subroutine minimize(etot,x,iretcode,nfun)
      implicit real*8 (a-h,o-z)
      include 'DIMENSIONS'
      parameter (liv=60,lv=(77+maxvar*(maxvar+17)/2)) 
*********************************************************************
* OPTIMIZE sets up SUMSL or DFP and provides a simple interface for *
* the calling subprogram.                                           *     
* when d(i)=1.0, then v(35) is the length of the initial step,      *     
* calculated in the usual pythagorean way.                          *     
* absolute convergence occurs when the function is within v(31) of  *     
* zero. unless you know the minimum value in advance, abs convg     *     
* is probably not useful.                                           *     
* relative convergence is when the model predicts that the function *   
* will decrease by less than v(32)*abs(fun).                        *   
*********************************************************************
      include 'COMMON.IOUNITS'
      include 'COMMON.VAR'
      include 'COMMON.GEO'
      include 'COMMON.MINIM'
      include 'COMMON.CONTROL'
      common /srutu/ icall
      dimension iv(liv)                                               
      double precision minval,x(maxvar),d(maxvar),v(1:lv),xx(maxvar)
      double precision energia(0:n_ene)
      external func,gradient,fdum
      external func_restr,grad_restr
      logical not_done,change,reduce 
c      common /przechowalnia/ v

      icall = 1

      NOT_DONE=.TRUE.

c     DO WHILE (NOT_DONE)

      call deflt(2,iv,liv,lv,v)                                         
* 12 means fresh start, dont call deflt                                 
      iv(1)=12                                                          
* max num of fun calls                                                  
      if (maxfun.eq.0) maxfun=500
      iv(17)=maxfun
* max num of iterations                                                 
      if (maxmin.eq.0) maxmin=1000
      iv(18)=maxmin
* controls output                                                       
      iv(19)=2                                                          
* selects output unit                                                   
      iv(21)=0
      if (print_min_ini+print_min_stat+print_min_res.gt.0) iv(21)=iout
* 1 means to print out result                                           
      iv(22)=print_min_res
* 1 means to print out summary stats                                    
      iv(23)=print_min_stat
* 1 means to print initial x and d                                      
      iv(24)=print_min_ini
* min val for v(radfac) default is 0.1                                  
      v(24)=0.1D0                                                       
* max val for v(radfac) default is 4.0                                  
      v(25)=2.0D0                                                       
c     v(25)=4.0D0                                                       
* check false conv if (act fnctn decrease) .lt. v(26)*(exp decrease)    
* the sumsl default is 0.1                                              
      v(26)=0.1D0
* false conv if (act fnctn decrease) .lt. v(34)                         
* the sumsl default is 100*machep                                       
      v(34)=v(34)/100.0D0                                               
* absolute convergence                                                  
      if (tolf.eq.0.0D0) tolf=1.0D-4
      v(31)=tolf
* relative convergence                                                  
      if (rtolf.eq.0.0D0) rtolf=1.0D-4
      v(32)=rtolf
* controls initial step size                                            
       v(35)=1.0D-1                                                    
* large vals of d correspond to small components of step                
      do i=1,nphi
        d(i)=1.0D-1
      enddo
      do i=nphi+1,nvar
        d(i)=1.0D-1
      enddo
cd    print *,'Calling SUMSL'
c     call var_to_geom(nvar,x)
c     call chainbuild
c     call etotal(energia(0))
c     etot = energia(0)
c      icheckgrad=3
c      call exec_checkgrad
      IF (mask_r) THEN
       call x2xx(x,xx,nvar_restr)
       call sumsl(nvar_restr,d,xx,func_restr,grad_restr,
     &                    iv,liv,lv,v,idum,rdum,fdum)      
       call xx2x(x,xx)
      ELSE
       call sumsl(nvar,d,x,func,gradient,iv,liv,lv,v,idum,rdum,fdum)      
      ENDIF
c      icheckgrad=3
c      call exec_checkgrad
      etot=v(10)                                                      
      iretcode=iv(1)
cd    print *,'Exit SUMSL; return code:',iretcode,' energy:',etot
cd    write (iout,'(/a,i4/)') 'SUMSL return code:',iv(1)
c     call intout
c     change=reduce(x)
      call var_to_geom(nvar,x)
c     if (change) then
c       write (iout,'(a)') 'Reduction worked, minimizing again...'
c     else
c       not_done=.false.
c     endif
      call chainbuild
c     call etotal(energia(0))
c     etot=energia(0)
c     call enerprint(energia(0))
      nfun=iv(6)

c     write (*,*) 'Processor',MyID,' leaves MINIMIZE.'

c     ENDDO ! NOT_DONE

      return  
      end  
#ifdef MPI
c----------------------------------------------------------------------------
      subroutine ergastulum
      implicit real*8 (a-h,o-z)
      include 'DIMENSIONS'
#ifdef MPI
      include "mpif.h"
#endif
      include 'COMMON.SETUP'
      include 'COMMON.DERIV'
      include 'COMMON.VAR'
      include 'COMMON.IOUNITS'
      include 'COMMON.FFIELD'
      include 'COMMON.INTERACT'
      include 'COMMON.MD'
      include 'COMMON.TIME1'
      double precision z(maxres6),d_a_tmp(maxres6)
      double precision edum(0:n_ene),time_order(0:10)
      double precision Gcopy(maxres2,maxres2)
      common /przechowalnia/ Gcopy
      integer icall /0/
C Workers wait for variables and NF, and NFL from the boss 
      iorder=0
      do while (iorder.ge.0)
c      write (*,*) 'Processor',fg_rank,' CG group',kolor,
c     & ' receives order from Master'
        time00=MPI_Wtime()
        call MPI_Bcast(iorder,1,MPI_INTEGER,king,FG_COMM,IERR)
        time_Bcast=time_Bcast+MPI_Wtime()-time00
        if (icall.gt.4 .and. iorder.ge.0) 
     &   time_order(iorder)=time_order(iorder)+MPI_Wtime()-time00
       icall=icall+1
c      write (*,*) 
c     & 'Processor',fg_rank,' completed receive MPI_BCAST order',iorder
        if (iorder.eq.0) then
          call zerograd
          call etotal(edum)
c          write (2,*) "After etotal"
c          write (2,*) "dimen",dimen," dimen3",dimen3
c          call flush(2)
        else if (iorder.eq.2) then
          call zerograd
          call etotal_short(edum)
c          write (2,*) "After etotal_short"
c          write (2,*) "dimen",dimen," dimen3",dimen3
c          call flush(2)
        else if (iorder.eq.3) then
          call zerograd
          call etotal_long(edum)
c          write (2,*) "After etotal_long"
c          write (2,*) "dimen",dimen," dimen3",dimen3
c          call flush(2)
        else if (iorder.eq.1) then
          call sum_gradient
c          write (2,*) "After sum_gradient"
c          write (2,*) "dimen",dimen," dimen3",dimen3
c          call flush(2)
        else if (iorder.eq.4) then
          call ginv_mult(z,d_a_tmp)
        else if (iorder.eq.5) then
c Setup MD things for a slave
          dimen=(nct-nnt+1)+nside
          dimen1=(nct-nnt)+(nct-nnt+1)
          dimen3=dimen*3
c          write (2,*) "dimen",dimen," dimen3",dimen3
c          call flush(2)
          call int_bounds(dimen,igmult_start,igmult_end)
          igmult_start=igmult_start-1
          call MPI_Allgather(3*igmult_start,1,MPI_INTEGER,
     &     ng_start(0),1,MPI_INTEGER,FG_COMM,IERROR)
           my_ng_count=igmult_end-igmult_start
          call MPI_Allgather(3*my_ng_count,1,MPI_INTEGER,ng_counts(0),1,
     &     MPI_INTEGER,FG_COMM,IERROR)
c          write (2,*) "ng_start",(ng_start(i),i=0,nfgtasks-1)
c          write (2,*) "ng_counts",(ng_counts(i),i=0,nfgtasks-1)
          myginv_ng_count=maxres2*my_ng_count
c          write (2,*) "igmult_start",igmult_start," igmult_end",
c     &     igmult_end," my_ng_count",my_ng_count
c          call flush(2)
          call MPI_Allgather(maxres2*igmult_start,1,MPI_INTEGER,
     &     nginv_start(0),1,MPI_INTEGER,FG_COMM,IERROR)
          call MPI_Allgather(myginv_ng_count,1,MPI_INTEGER,
     &     nginv_counts(0),1,MPI_INTEGER,FG_COMM,IERROR)
c          write (2,*) "nginv_start",(nginv_start(i),i=0,nfgtasks-1)
c          write (2,*) "nginv_counts",(nginv_counts(i),i=0,nfgtasks-1)
c          call flush(2)
c          call MPI_Barrier(FG_COMM,IERROR)
          time00=MPI_Wtime()
          call MPI_Scatterv(ginv(1,1),nginv_counts(0),
     &     nginv_start(0),MPI_DOUBLE_PRECISION,gcopy(1,1),
     &     myginv_ng_count,MPI_DOUBLE_PRECISION,king,FG_COMM,IERR)
#ifdef TIMING
          time_scatter_ginv=time_scatter_ginv+MPI_Wtime()-time00
#endif
          do i=1,dimen
            do j=1,2*my_ng_count
              ginv(j,i)=gcopy(i,j)
            enddo
          enddo
c          write (2,*) "dimen",dimen," dimen3",dimen3
c          write (2,*) "End MD setup"
c          call flush(2)
c           write (iout,*) "My chunk of ginv_block"
c           call MATOUT2(my_ng_count,dimen3,maxres2,maxers2,ginv_block)
        else if (iorder.eq.6) then
          call int_from_cart1(.false.)
        else if (iorder.eq.7) then
          call chainbuild_cart
        else if (iorder.eq.8) then
          call intcartderiv
        else if (iorder.eq.9) then
          call fricmat_mult(z,d_a_tmp)
        else if (iorder.eq.10) then
          call setup_fricmat
        endif
      enddo
      write (*,*) 'Processor',fg_rank,' CG group',kolor,
     &  ' absolute rank',myrank,' leves ERGASTULUM.'
      write(*,*)'Processor',fg_rank,' wait times for respective orders',
     & (' order[',i,']',time_order(i),i=0,10)
      return
      end
#endif
************************************************************************
      subroutine func(n,x,nf,f,uiparm,urparm,ufparm)  
      implicit real*8 (a-h,o-z)
      include 'DIMENSIONS'
      include 'COMMON.DERIV'
      include 'COMMON.IOUNITS'
      include 'COMMON.GEO'
      common /chuju/ jjj
      double precision energia(0:n_ene)
      integer jjj
      double precision ufparm
      external ufparm                                                   
      integer uiparm(1)                                                 
      real*8 urparm(1)                                                    
      dimension x(maxvar)
c     if (jjj.gt.0) then
c       write (iout,'(10f8.3)') (rad2deg*x(i),i=1,n)
c     endif
      nfl=nf
      icg=mod(nf,2)+1
cd      print *,'func',nf,nfl,icg
      call var_to_geom(n,x)
      call zerograd
      call chainbuild
cd    write (iout,*) 'ETOTAL called from FUNC'
      call etotal(energia(0))
      call sum_gradient
      f=energia(0)
c     if (jjj.gt.0) then
c       write (iout,'(10f8.3)') (rad2deg*x(i),i=1,n)
c       write (iout,*) 'f=',etot
c       jjj=0
c     endif
      return                                                            
      end                                                               
************************************************************************
      subroutine func_restr(n,x,nf,f,uiparm,urparm,ufparm)  
      implicit real*8 (a-h,o-z)
      include 'DIMENSIONS'
      include 'COMMON.DERIV'
      include 'COMMON.IOUNITS'
      include 'COMMON.GEO'
      common /chuju/ jjj
      double precision energia(0:n_ene)
      integer jjj
      double precision ufparm
      external ufparm                                                   
      integer uiparm(1)                                                 
      real*8 urparm(1)                                                    
      dimension x(maxvar)
c     if (jjj.gt.0) then
c       write (iout,'(10f8.3)') (rad2deg*x(i),i=1,n)
c     endif
      nfl=nf
      icg=mod(nf,2)+1
      call var_to_geom_restr(n,x)
      call zerograd
      call chainbuild
cd    write (iout,*) 'ETOTAL called from FUNC'
      call etotal(energia(0))
      call sum_gradient
      f=energia(0)
c     if (jjj.gt.0) then
c       write (iout,'(10f8.3)') (rad2deg*x(i),i=1,n)
c       write (iout,*) 'f=',etot
c       jjj=0
c     endif
      return                                                            
      end                                                               
c-------------------------------------------------------
      subroutine x2xx(x,xx,n)
      implicit real*8 (a-h,o-z)
      include 'DIMENSIONS'
      include 'COMMON.VAR'
      include 'COMMON.CHAIN'
      include 'COMMON.INTERACT'
      double precision xx(maxvar),x(maxvar)

      do i=1,nvar
        varall(i)=x(i)
      enddo

      ig=0                                                                      
      igall=0                                                                   
      do i=4,nres                                                               
        igall=igall+1                                                           
        if (mask_phi(i).eq.1) then                                              
          ig=ig+1                                                               
          xx(ig)=x(igall)                       
        endif                                                                   
      enddo                                                                     
                                                                                
      do i=3,nres                                                               
        igall=igall+1                                                           
        if (mask_theta(i).eq.1) then                                            
          ig=ig+1                                                               
          xx(ig)=x(igall)
        endif                                                                   
      enddo                                          

      do ij=1,2                                                                 
      do i=2,nres-1                                                             
        if (itype(i).ne.10) then                                                
          igall=igall+1                                                         
          if (mask_side(i).eq.1) then                                           
            ig=ig+1                                                             
            xx(ig)=x(igall)
          endif                                                                 
        endif                                                                   
      enddo                                                                     
      enddo                              
 
      n=ig

      return
      end

      subroutine xx2x(x,xx)
      implicit real*8 (a-h,o-z)
      include 'DIMENSIONS'
      include 'COMMON.VAR'
      include 'COMMON.CHAIN'
      include 'COMMON.INTERACT'
      double precision xx(maxvar),x(maxvar)

      do i=1,nvar
        x(i)=varall(i)
      enddo

      ig=0                                                                      
      igall=0                                                                   
      do i=4,nres                                                               
        igall=igall+1                                                           
        if (mask_phi(i).eq.1) then                                              
          ig=ig+1                                                               
          x(igall)=xx(ig)
        endif                                                                   
      enddo                                                                     
                                                                                
      do i=3,nres                                                               
        igall=igall+1                                                           
        if (mask_theta(i).eq.1) then                                            
          ig=ig+1                                                               
          x(igall)=xx(ig)
        endif                                                                   
      enddo                                          

      do ij=1,2                                                                 
      do i=2,nres-1                                                             
        if (itype(i).ne.10) then                                                
          igall=igall+1                                                         
          if (mask_side(i).eq.1) then                                           
            ig=ig+1                                                             
            x(igall)=xx(ig)
          endif                                                                 
        endif                                                                   
      enddo                                                             
      enddo                              

      return
      end
  
c---------------------------------------------------------- 
      subroutine minim_dc(etot,iretcode,nfun)
      implicit real*8 (a-h,o-z)
      include 'DIMENSIONS'
      parameter (liv=60,lv=(77+maxvar*(maxvar+17)/2)) 
#ifdef MPI
      include 'mpif.h'
#endif
      include 'COMMON.SETUP'
      include 'COMMON.IOUNITS'
      include 'COMMON.VAR'
      include 'COMMON.GEO'
      include 'COMMON.MINIM'
      include 'COMMON.CHAIN'
      dimension iv(liv)                                               
      double precision minval,x(maxvar),d(maxvar),v(1:lv),xx(maxvar)
c      common /przechowalnia/ v

      double precision energia(0:n_ene)
      external func_dc,grad_dc,fdum
      logical not_done,change,reduce 
      double precision g(maxvar),f1

      call deflt(2,iv,liv,lv,v)                                         
* 12 means fresh start, dont call deflt                                 
      iv(1)=12                                                          
* max num of fun calls                                                  
      if (maxfun.eq.0) maxfun=500
      iv(17)=maxfun
* max num of iterations                                                 
      if (maxmin.eq.0) maxmin=1000
      iv(18)=maxmin
* controls output                                                       
      iv(19)=2                                                          
* selects output unit                                                   
      iv(21)=0
      if (print_min_ini+print_min_stat+print_min_res.gt.0) iv(21)=iout 
* 1 means to print out result                                           
      iv(22)=print_min_res
* 1 means to print out summary stats                                    
      iv(23)=print_min_stat
* 1 means to print initial x and d                                      
      iv(24)=print_min_ini
* min val for v(radfac) default is 0.1                                  
      v(24)=0.1D0                                                       
* max val for v(radfac) default is 4.0                                  
      v(25)=2.0D0                                                       
c     v(25)=4.0D0                                                       
* check false conv if (act fnctn decrease) .lt. v(26)*(exp decrease)    
* the sumsl default is 0.1                                              
      v(26)=0.1D0
* false conv if (act fnctn decrease) .lt. v(34)                         
* the sumsl default is 100*machep                                       
      v(34)=v(34)/100.0D0                                               
* absolute convergence                                                  
      if (tolf.eq.0.0D0) tolf=1.0D-4
      v(31)=tolf
* relative convergence                                                  
      if (rtolf.eq.0.0D0) rtolf=1.0D-4
      v(32)=rtolf
* controls initial step size                                            
       v(35)=1.0D-1                                                    
* large vals of d correspond to small components of step                
      do i=1,6*nres
        d(i)=1.0D-1
      enddo

      k=0
      do i=1,nres-1
        do j=1,3
          k=k+1
          x(k)=dc(j,i)
        enddo
      enddo
      do i=2,nres-1
        if (ialph(i,1).gt.0) then
        do j=1,3
          k=k+1
          x(k)=dc(j,i+nres)
        enddo
        endif
      enddo
      print *,"check_ecart before sumsl"
c      icheckgrad=2
c      call exec_checkgrad

      call sumsl(k,d,x,func_dc,grad_dc,iv,liv,lv,v,idum,rdum,fdum)      

      k=0
      do i=1,nres-1
        do j=1,3
          k=k+1
          dc(j,i)=x(k)
        enddo
      enddo
      do i=2,nres-1
        if (ialph(i,1).gt.0) then
        do j=1,3
          k=k+1
          dc(j,i+nres)=x(k)
        enddo
        endif
      enddo
      call chainbuild_cart

cd      call zerograd
cd      nf=0
cd      call func_dc(k,x,nf,f,idum,rdum,fdum)
cd      call grad_dc(k,x,nf,g,idum,rdum,fdum)
cd
cd      do i=1,k
cd       x(i)=x(i)+1.0D-5
cd       call func_dc(k,x,nf,f1,idum,rdum,fdum)
cd       x(i)=x(i)-1.0D-5
cd       print '(i5,2f15.5)',i,g(i),(f1-f)/1.0D-5
cd      enddo

      etot=v(10)                                                      
      iretcode=iv(1)
      nfun=iv(6)
      print *,"check_ecart"
c      icheckgrad=2
c      call exec_checkgrad
      return  
      end  

      subroutine func_dc(n,x,nf,f,uiparm,urparm,ufparm)  
      implicit real*8 (a-h,o-z)
      include 'DIMENSIONS'
#ifdef MPI
      include 'mpif.h'
#endif
      include 'COMMON.SETUP'
      include 'COMMON.DERIV'
      include 'COMMON.IOUNITS'
      include 'COMMON.GEO'
      include 'COMMON.CHAIN'
      include 'COMMON.VAR'
      double precision energia(0:n_ene)
      double precision ufparm
      external ufparm                                                   
      integer uiparm(1)                                                 
      real*8 urparm(1)                                                    
      dimension x(maxvar)
      nfl=nf
cbad      icg=mod(nf,2)+1
      icg=1

      k=0
      do i=1,nres-1
        do j=1,3
          k=k+1
          dc(j,i)=x(k)
        enddo
      enddo
      do i=2,nres-1
        if (ialph(i,1).gt.0) then
        do j=1,3
          k=k+1
          dc(j,i+nres)=x(k)
        enddo
        endif
      enddo
      call chainbuild_cart

      call zerograd
      call etotal(energia(0))
      f=energia(0)

cd      print *,'func_dc ',nf,nfl,f

      return                                                            
      end                                                               

      subroutine grad_dc(n,x,nf,g,uiparm,urparm,ufparm)
      implicit real*8 (a-h,o-z)
      include 'DIMENSIONS'
#ifdef MPI
      include 'mpif.h'
#endif
      include 'COMMON.SETUP'
      include 'COMMON.CHAIN'
      include 'COMMON.DERIV'
      include 'COMMON.VAR'
      include 'COMMON.INTERACT'
      include 'COMMON.FFIELD'
      include 'COMMON.MD'
      include 'COMMON.IOUNITS'
      external ufparm
      integer uiparm(1),k
      double precision urparm(1)
      dimension x(maxvar),g(maxvar)
c
c
c
cbad      icg=mod(nf,2)+1
      icg=1
cd      print *,'grad_dc ',nf,nfl,nf-nfl+1,icg
      if (nf-nfl+1) 20,30,40
   20 call func_dc(n,x,nf,f,uiparm,urparm,ufparm)
cd      print *,20
      if (nf.eq.0) return
      goto 40
   30 continue
cd      print *,30
      k=0
      do i=1,nres-1
        do j=1,3
          k=k+1
          dc(j,i)=x(k)
        enddo
      enddo
      do i=2,nres-1
        if (ialph(i,1).gt.0) then
        do j=1,3
          k=k+1
          dc(j,i+nres)=x(k)
        enddo
        endif
      enddo
      call chainbuild_cart

C
C Evaluate the derivatives of virtual bond lengths and SC vectors in variables.
C
   40 call cartgrad
cd      print *,40
      k=0
      do i=1,nres-1
        do j=1,3
          k=k+1
          g(k)=gcart(j,i)
        enddo
      enddo
      do i=2,nres-1
        if (ialph(i,1).gt.0) then
        do j=1,3
          k=k+1
          g(k)=gxcart(j,i)
        enddo
        endif
      enddo       

      return
      end