--- /dev/null
+ subroutine minim_jlee
+c controls minimization and sorting routines
+ implicit real*8 (a-h,o-z)
+ include 'DIMENSIONS'
+ parameter (liv=60,lv=(77+maxvar*(maxvar+17)/2))
+ include 'COMMON.VAR'
+ include 'COMMON.IOUNITS'
+ include 'COMMON.MINIM'
+ include 'COMMON.CONTROL'
+ include 'mpif.h'
+ external func,gradient,fdum
+ real ran1,ran2,ran3
+ include 'COMMON.SETUP'
+ include 'COMMON.GEO'
+ include 'COMMON.FFIELD'
+ include 'COMMON.SBRIDGE'
+ include 'COMMON.DISTFIT'
+ include 'COMMON.CHAIN'
+ dimension muster(mpi_status_size)
+ dimension var(maxvar),erg(mxch*(mxch+1)/2+1)
+ dimension var2(maxvar)
+ integer iffr(maxres),ihpbt(maxdim),jhpbt(maxdim)
+ double precision d(maxvar),v(1:lv+1),garbage(maxvar)
+ double precision energia(0:n_ene),time0s,time1s
+ dimension indx(9),info(12)
+ dimension iv(liv)
+ dimension idum(1),rdum(1)
+ dimension icont(2,maxcont)
+ logical check_var,fail
+ integer iloop(2)
+ common /przechowalnia/ v
+ data rad /1.745329252d-2/
+c receive # of start
+! print *,'Processor',me,' calling MINIM_JLEE maxfun',maxfun,
+! & ' maxmin',maxmin,' tolf',tolf,' rtolf',rtolf
+ nhpb0=nhpb
+ 10 continue
+ time0s=MPI_WTIME()
+c print *, 'MINIM_JLEE: ',me,' is waiting'
+ call mpi_recv(info,12,mpi_integer,king,idint,CG_COMM,
+ * muster,ierr)
+ time1s=MPI_WTIME()
+ write (iout,'(a12,f10.4,a4)')'Waiting for ',time1s-time0s,' sec'
+ call flush(iout)
+ n=info(1)
+c print *, 'MINIM_JLEE: ',me,' received: ',n
+
+crc if (ierr.ne.0) go to 100
+c if # = 0, return
+ if (n.eq.0) then
+ write (iout,*) 'Finishing minim_jlee - signal',n,' from master'
+ call flush(iout)
+ return
+ endif
+
+ nfun=0
+ IF (n.lt.0) THEN
+ call mpi_recv(var,nvar,mpi_double_precision,
+ * king,idreal,CG_COMM,muster,ierr)
+ call mpi_recv(iffr,nres,mpi_integer,
+ * king,idint,CG_COMM,muster,ierr)
+ call mpi_recv(var2,nvar,mpi_double_precision,
+ * king,idreal,CG_COMM,muster,ierr)
+ ELSE
+c receive initial values of variables
+ call mpi_recv(var,nvar,mpi_double_precision,
+ * king,idreal,CG_COMM,muster,ierr)
+crc if (ierr.ne.0) go to 100
+ ENDIF
+
+ if(vdisulf.and.info(2).ne.-1) then
+ if(info(4).ne.0)then
+ call mpi_recv(ihpbt,info(4),mpi_integer,
+ * king,idint,CG_COMM,muster,ierr)
+ call mpi_recv(jhpbt,info(4),mpi_integer,
+ * king,idint,CG_COMM,muster,ierr)
+ endif
+ endif
+
+ IF (n.lt.0) THEN
+ n=-n
+ nhpb=nhpb0
+ link_start=1
+ link_end=nhpb
+ call init_int_table
+ call contact_cp(var,var2,iffr,nfun,n)
+ ENDIF
+
+ if(vdisulf.and.info(2).ne.-1) then
+ nss=0
+ if(info(4).ne.0)then
+cd write(iout,*) 'SS=',info(4),'N=',info(1),'IT=',info(2)
+ call var_to_geom(nvar,var)
+ call chainbuild
+ do i=1,info(4)
+ if (dist(ihpbt(i),jhpbt(i)).lt.7.0) then
+ nss=nss+1
+ ihpb(nss)=ihpbt(i)
+ jhpb(nss)=jhpbt(i)
+cd write(iout,*) 'SS mv=',info(3),
+cd & ihpb(nss)-nres,jhpb(nss)-nres,
+cd & dist(ihpb(nss),jhpb(nss))
+ dhpb(nss)=dbr
+ forcon(nss)=fbr
+ else
+cd write(iout,*) 'rm SS mv=',info(3),
+cd & ihpbt(i)-nres,jhpbt(i)-nres,dist(ihpbt(i),jhpbt(i))
+ endif
+ enddo
+ endif
+ nhpb=nss
+ link_start=1
+ link_end=nhpb
+ call init_int_table
+ endif
+
+ if (info(3).eq.14) then
+ write(iout,*) 'calling local_move',info(7),info(8)
+ call local_move_init(.false.)
+ call var_to_geom(nvar,var)
+ call local_move(info(7),info(8),20d0,50d0)
+ call geom_to_var(nvar,var)
+ endif
+
+
+ if (info(3).eq.16) then
+ write(iout,*) 'calling beta_slide',info(7),info(8),
+ & info(10), info(11), info(12)
+ call var_to_geom(nvar,var)
+ call beta_slide(info(7),info(8),info(10),info(11),info(12)
+ & ,nfun,n)
+ call geom_to_var(nvar,var)
+ endif
+
+
+ if (info(3).eq.17) then
+ write(iout,*) 'calling beta_zip',info(7),info(8)
+ call var_to_geom(nvar,var)
+ call beta_zip(info(7),info(8),nfun,n)
+ call geom_to_var(nvar,var)
+ endif
+
+
+crc overlap test
+
+ if (overlapsc) then
+
+ call var_to_geom(nvar,var)
+ call chainbuild
+ call etotal(energia(0))
+ nfun=nfun+1
+ if (energia(1).eq.1.0d20) then
+ info(3)=-info(3)
+ write (iout,'(a,1pe14.5)')'#OVERLAP evdw=1d20',energia(1)
+ call overlap_sc(fail)
+ if(.not.fail) then
+ call geom_to_var(nvar,var)
+ call etotal(energia(0))
+ nfun=nfun+1
+ write (iout,'(a,1pe14.5)')'#OVERLAP evdw after',energia(1)
+ else
+ v(10)=1.0d20
+ iv(1)=-1
+ goto 201
+ endif
+ endif
+ endif
+
+ if (searchsc) then
+ call var_to_geom(nvar,var)
+ call sc_move(2,nres-1,1,10d0,nft_sc,etot)
+ call geom_to_var(nvar,var)
+cd write(iout,*) 'sc_move',nft_sc,etot
+ endif
+
+ if (check_var(var,info)) then
+ v(10)=1.0d21
+ iv(1)=6
+ goto 201
+ endif
+
+
+crc
+
+! write (iout,*) 'MINIM_JLEE: Processor',me,' nvar',nvar
+! write (iout,'(8f10.4)') (var(i),i=1,nvar)
+! write (*,*) 'MINIM_JLEE: Processor',me,' received nvar',nvar
+! write (*,'(8f10.4)') (var(i),i=1,nvar)
+
+ do i=1,nvar
+ garbage(i)=var(i)
+ enddo
+
+ 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
+cd iv(21)=iout
+ iv(21)=0
+* 1 means to print out result
+ iv(22)=0
+cd iv(22)=1
+* 1 means to print out summary stats
+ iv(23)=0
+* 1 means to print initial x and d
+ iv(24)=0
+
+c if(me.eq.3.and.n.eq.255) then
+c print *,' CHUJ: stoi'
+c iv(21)=6
+c iv(22)=1
+c iv(23)=1
+c iv(24)=1
+c endif
+
+* 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
+c minimize energy
+! write (iout,*) 'Processor',me,' nvar',nvar
+! write (iout,*) 'Variables BEFORE minimization:'
+! write (iout,'(8f10.4)') (rad2deg*var(i),i=1,nvar)
+
+c print *, 'MINIM_JLEE: ',me,' before SUMSL '
+
+ call func(nvar,var,nf,eee,idum,rdum,fdum)
+ nfun=nfun+1
+ if(eee.ge.1.0d20) then
+c print *,'MINIM_JLEE: ',me,' CHUJ NASTAPIL'
+c print *,' energy before SUMSL =',eee
+c print *,' aborting local minimization'
+ iv(1)=-1
+ v(10)=eee
+ go to 201
+ endif
+
+ct time0s=MPI_WTIME()
+ call sumsl(nvar,d,var,func,gradient,iv,liv,lv,v,idum,rdum,fdum)
+ct write(iout,*) 'sumsl time=',MPI_WTIME()-time0s,iv(7),v(10)
+c print *, 'MINIM_JLEE: ',me,' after SUMSL '
+
+c find which conformation was returned from sumsl
+ nfun=nfun+iv(7)
+! print *,'Processor',me,' iv(17)',iv(17),' iv(18)',iv(18),' nf',nf,
+! & ' retcode',iv(1),' energy',v(10),' tolf',v(31),' rtolf',v(32)
+c if (iv(1).ne.4 .or. nf.le.1) then
+c write (*,*) 'Processor',me,' something bad in SUMSL',iv(1),nf
+c write (*,*) 'Initial Variables'
+c write (*,'(8f10.4)') (rad2deg*garbage(i),i=1,nvar)
+c write (*,*) 'Variables'
+c write (*,'(8f10.4)') (rad2deg*var(i),i=1,nvar)
+c write (*,*) 'Vector d'
+c write (*,'(8f10.4)') (d(i),i=1,nvar)
+c write (iout,*) 'Processor',me,' something bad in SUMSL',
+c & iv(1),nf
+c write (iout,*) 'Initial Variables'
+c write (iout,'(8f10.4)') (rad2deg*garbage(i),i=1,nvar)
+c write (iout,*) 'Variables'
+c write (iout,'(8f10.4)') (rad2deg*var(i),i=1,nvar)
+c write (iout,*) 'Vector d'
+c write (iout,'(8f10.4)') (d(i),i=1,nvar)
+c endif
+c if (nf.lt.iv(6)-1) then
+c recalculate intra- and interchain energies
+c call func(nvar,var,nf,v(10),iv,v,fdum)
+c else if (nf.eq.iv(6)-1) then
+c regenerate conformation
+c call var_to_geom(nvar,var)
+c call chainbuild
+c endif
+c change origin and axes to standard ECEPP format
+c call var_to_geom(nvar,var)
+! write (iout,*) 'MINIM_JLEE after minim: Processor',me,' nvar',nvar
+! write (iout,'(8f10.4)') (var(i),i=1,nvar)
+! write (iout,*) 'Energy:',v(10)
+c send back output
+c print *, 'MINIM_JLEE: ',me,' minimized: ',n
+ 201 continue
+ indx(1)=n
+c return code: 6-gradient 9-number of ftn evaluation, etc
+ indx(2)=iv(1)
+c total # of ftn evaluations (for iwf=0, it includes all minimizations).
+ indx(3)=nfun
+ indx(4)=info(2)
+ indx(5)=info(3)
+ indx(6)=nss
+ indx(7)=info(5)
+ indx(8)=info(6)
+ indx(9)=info(9)
+ call mpi_send(indx,9,mpi_integer,king,idint,CG_COMM,
+ * ierr)
+c send back energies
+c al & cc
+c calculate contact order
+#ifdef CO_BIAS
+ call contact(.false.,ncont,icont,co)
+ erg(1)=v(10)-1.0d2*co
+#else
+ erg(1)=v(10)
+#endif
+ j=1
+ call mpi_send(erg,j,mpi_double_precision,king,idreal,
+ * CG_COMM,ierr)
+#ifdef CO_BIAS
+ call mpi_send(co,j,mpi_double_precision,king,idreal,
+ * CG_COMM,ierr)
+#endif
+c send back values of variables
+ call mpi_send(var,nvar,mpi_double_precision,
+ * king,idreal,CG_COMM,ierr)
+! print * , 'MINIM_JLEE: Processor',me,' send erg and var '
+
+ if(vdisulf.and.info(2).ne.-1.and.nss.ne.0) then
+cd call intout
+cd call chainbuild
+cd call etotal(energia(0))
+cd etot=energia(0)
+cd call enerprint(energia(0))
+ call mpi_send(ihpb,nss,mpi_integer,
+ * king,idint,CG_COMM,ierr)
+ call mpi_send(jhpb,nss,mpi_integer,
+ * king,idint,CG_COMM,ierr)
+ endif
+
+ go to 10
+ 100 print *, ' error in receiving message from emperor', me
+ call mpi_abort(mpi_comm_world,ierror,ierrcode)
+ return
+ 200 print *, ' error in sending message to emperor'
+ call mpi_abort(mpi_comm_world,ierror,ierrcode)
+ return
+ 300 print *, ' error in communicating with emperor'
+ call mpi_abort(mpi_comm_world,ierror,ierrcode)
+ return
+ 956 format (' initial energy could not be calculated',41x)
+ 957 format (80x)
+ 965 format (' convergence code ',i2,' # of function calls ',
+ * i4,' # of gradient calls ',i4,10x)
+ 975 format (' energy ',1p,e12.4,' scaled gradient ',e11.3,32x)
+ end
+
+ logical function check_var(var,info)
+ implicit real*8 (a-h,o-z)
+ include 'DIMENSIONS'
+ include 'COMMON.VAR'
+ include 'COMMON.IOUNITS'
+ include 'COMMON.GEO'
+ include 'COMMON.SETUP'
+ dimension var(maxvar)
+ dimension info(3)
+C AL -------
+ check_var=.false.
+ do i=nphi+ntheta+1,nphi+ntheta+nside
+! Check the side chain "valence" angles alpha
+ if (var(i).lt.1.0d-7) then
+ write (iout,*) 'CHUJ NASTAPIL ABSOLUTNY!!!!!!!!!!!!'
+ write (iout,*) 'Processor',me,'received bad variables!!!!'
+ write (iout,*) 'Variables'
+ write (iout,'(8f10.4)') (rad2deg*var(j),j=1,nvar)
+ write (iout,*) 'Continuing calculations at this point',
+ & ' could destroy the results obtained so far... ABORTING!!!!!!'
+ write (iout,'(a19,i5,f10.4,a4,2i4,a3,i3)')
+ & 'valence angle alpha',i-nphi-ntheta,var(i),
+ & 'n it',info(1),info(2),'mv ',info(3)
+ write (*,*) 'CHUJ NASTAPIL ABSOLUTNY!!!!!!!!!!!!'
+ write (*,*) 'Processor',me,'received bad variables!!!!'
+ write (*,*) 'Variables'
+ write (*,'(8f10.4)') (rad2deg*var(j),j=1,nvar)
+ write (*,*) 'Continuing calculations at this point',
+ & ' could destroy the results obtained so far... ABORTING!!!!!!'
+ write (*,'(a19,i5,f10.4,a4,2i4,a3,i3)')
+ & 'valence angle alpha',i-nphi-ntheta,var(i),
+ & 'n it',info(1),info(2),'mv ',info(3)
+ check_var=.true.
+ return
+ endif
+ enddo
+! Check the backbone "valence" angles theta
+ do i=nphi+1,nphi+ntheta
+ if (var(i).lt.1.0d-7) then
+ write (iout,*) 'CHUJ NASTAPIL ABSOLUTNY!!!!!!!!!!!!'
+ write (iout,*) 'Processor',me,'received bad variables!!!!'
+ write (iout,*) 'Variables'
+ write (iout,'(8f10.4)') (rad2deg*var(j),j=1,nvar)
+ write (iout,*) 'Continuing calculations at this point',
+ & ' could destroy the results obtained so far... ABORTING!!!!!!'
+ write (iout,'(a19,i5,f10.4,a4,2i4,a3,i3)')
+ & 'valence angle theta',i-nphi,var(i),
+ & 'n it',info(1),info(2),'mv ',info(3)
+ write (*,*) 'CHUJ NASTAPIL ABSOLUTNY!!!!!!!!!!!!'
+ write (*,*) 'Processor',me,'received bad variables!!!!'
+ write (*,*) 'Variables'
+ write (*,'(8f10.4)') (rad2deg*var(j),j=1,nvar)
+ write (*,*) 'Continuing calculations at this point',
+ & ' could destroy the results obtained so far... ABORTING!!!!!!'
+ write (*,'(a19,i5,f10.4,a4,2i4,a3,i3)')
+ & 'valence angle theta',i-nphi,var(i),
+ & 'n it',info(1),info(2),'mv ',info(3)
+ check_var=.true.
+ return
+ endif
+ enddo
+ return
+ end
projects / unres.git / blobdiff