1 subroutine minimize_init(n,iv,v,d)
4 include "DIMENSIONS.ZSCOPT"
6 parameter (liv=60,lv=(77+max_paropt*(max_paropt+17)/2))
8 include "COMMON.MINPAR"
9 include "COMMON.IOUNITS"
10 include "COMMON.TIME1"
12 double precision d(max_paropt),v(1:lv)
16 call deflt(2,iv,liv,lv,v)
17 * 12 means fresh start, dont call deflt
19 * max num of fun calls
21 * max num of iterations
27 * 1 means to print out result
29 * 1 means to print out summary stats
31 * 1 means to print initial x and d
33 * min val for v(radfac) default is 0.1
35 * max val for v(radfac) default is 4.0
38 * check false conv if (act fnctn decrease) .lt. v(26)*(exp decrease)
39 * the sumsl default is 0.1
41 * false conv if (act fnctn decrease) .lt. v(34)
42 * the sumsl default is 100*machep
44 * absolute convergence
45 if (tolf.eq.0.0D0) tolf=1.0D-6
47 * relative convergence
48 if (rtolf.eq.0.0D0) rtolf=1.0D-6
50 * controls initial step size
52 * large vals of d correspond to small components of step
58 c------------------------------------------------------------------------------
59 subroutine minimize(n,x,f)
62 include "DIMENSIONS.ZSCOPT"
64 parameter (liv=60,lv=(77+max_paropt*(max_paropt+17)/2))
65 *********************************************************************
66 * OPTIMIZE sets up SUMSL or DFP and provides a simple interface for *
67 * the calling subprogram. *
68 * when d(i)=1.0, then v(35) is the length of the initial step, *
69 * calculated in the usual pythagorean way. *
70 * absolute convergence occurs when the function is within v(31) of *
71 * zero. unless you know the minimum value in advance, abs convg *
72 * is probably not useful. *
73 * relative convergence is when the model predicts that the function *
74 * will decrease by less than v(32)*abs(fun). *
75 *********************************************************************
80 integer status(MPI_STATUS_SIZE)
82 include "COMMON.WEIGHTS"
83 include "COMMON.WEIGHTDER"
84 include "COMMON.OPTIM"
85 include "COMMON.CLASSES"
86 include "COMMON.COMPAR"
87 include "COMMON.MINPAR"
88 include "COMMON.IOUNITS"
89 include "COMMON.VMCPAR"
90 include "COMMON.TIME1"
91 include "COMMON.ENERGIES"
92 include "COMMON.PROTNAME"
93 include "COMMON.THERMAL"
99 double precision x(max_paropt),d(max_paropt),v(1:lv),g(max_paropt)
100 real p(max_paropt+1,max_paropt),y(max_paropt+1),pb(max_paropt),
101 & pom(max_paropt),yb,y1,y2,pomi,ftol,temptr
104 external targetfun,targetgrad,fdum
105 integer uiparm,i,ii,j,k,nf,nfun,ngrad,iretcode,ianneal,maxanneal
106 double precision urparm,penalty
107 double precision f,f0,obf,obg(max_paropt),tcpu
108 double precision Evarmax(maxprot),aux
109 integer iscan,iprot,ibatch,niter
113 write (iout,*) "Minimizing with ",MINIMIZER
114 IF (MINIMIZER.EQ."AMEBSA") THEN
130 pom(i-1)=p(i,i-1)-0.2
132 if (pom(i-1).lt.x_low(i-1)) pom(i-1)=x_low(i-1)
134 pom(i-1)=p(i,i-1)+0.2
135 if (pom(i-1).gt.x_up(i-1)) pom(i-1)=x_up(i-1)
146 write (iout,*) "p and y"
148 write (iout,'(20e15.5)') (p(i,j),j=1,n),y(i)
151 do ianneal=1,maxanneal
154 write (iout,*) "ianneal",ianneal," temptr",temptr
155 call amebsa(p,y,max_paropt+1,max_paropt,n,pb,yb,ftol,funk,
157 write (iout,*) "Optimized P (PB)"
159 write (iout,'(i5,f10.5)') i,pb(i)
161 write (iout,*) "Optimal function value",yb
162 write(iout,*) "AMEBSA return code:",iretcode
167 call targetfun(n,x,nf,f,uiparm,fdum)
169 write(iout,*) "Final function value:",f,f
170 write(iout,*) "Final value of the object function:",obf
171 write(iout,*) "Number of target function evaluations:",nfun
177 open (88,file=restartname,status="unknown")
182 ELSE IF (MINIMIZER.EQ."SUMSL") THEN
187 call minimize_init(n,iv,v,d)
189 write (iout,*) "Initial checkgrad:"
191 c write (iout,*) "Second checkgrad"
192 c call checkgrad(n,x)
194 write (iout,*) "Calling SUMSL"
196 call targetfun(n,x,nf,f,uiparm,fdum)
197 write (iout,*) "Initial function value",f
200 do while ( iscan.lt.maxscan .and. (f0-f)/(f0+1.0d0).gt.0.2d0 )
202 call scan(n,x,f,.false.)
208 if (maxscan.gt.0) then
210 do while ( iscan.lt.maxscan .and. (f0-f)/(f0+1.0d0).gt.0.2d0 )
213 call sumsl(n,d,x,targetfun,targetgrad,iv,liv,lv,v,uiparm,
215 write (iout,*) "Left SUMSL"
220 call targetfun(n,x,nf,f,uiparm,fdum)
222 call scan(n,x,f,.false.)
224 if (iv(1).eq.11 .and. cutoffviol) then
227 &"Revising the list of conformations because of cut-off violation."
229 write (iout,'(a,f7.1)')
230 & protname(iprot)(:ilen(protname(iprot))),enecut(iprot)
231 do j=1,natlike(iprot)+2
232 write (iout,'(5(f7.1,f9.1,2x))')
233 & (1.0d0/(Rgas*betaT(k,ibatch,iprot)),
234 & Evar(j,k,iprot),k=1,nbeta(j,iprot))
241 call make_list(.true.,.false.,n,x)
244 else if (iv(1).eq.11 .and. WhatsUp.eq.1) then
245 write (iout,*) "Time up, terminating."
253 write (iout,*) "MAXMIN",maxmin
255 & call sumsl(n,d,x,targetfun,targetgrad,iv,liv,lv,v,uiparm,
257 write (iout,*) "Left SUMSL"
263 write (iout,*) "iretcode",iv(1)," cutoffviol",cutoffviol
264 if (iv(1).eq.11 .and. cutoffviol) then
266 &"Revising the list of conformations because of cut-off violation."
268 write (iout,'(a,f7.1)')
269 & protname(iprot)(:ilen(protname(iprot))),enecut(iprot)
270 do j=1,natlike(iprot)+2
271 write (iout,'(5(f7.1,f9.1,2x))')
272 & (1.0d0/(Rgas*betaT(k,ibatch,iprot)),
273 & Evar(j,k,iprot),k=1,nbeta(j,iprot))
278 Evarmax(iprot)=Evar(1,1,iprot)*betaT(1,1,iprot)
279 do ibatch=2,natlike(iprot)+2
280 do k=1,nbeta(ibatch,iprot)
281 aux = Evar(k,ibatch,iprot)*betaT(k,ibatch,iprot)
282 if (aux.gt.Evarmax(iprot)) Evarmax(iprot)=aux
287 aux=dmin1(0.5d0*enecut_min(iprot)+1.25d0*Evarmax(iprot)
288 & ,enecut_max(iprot))
289 enecut(iprot)=dmax1(aux,enecut_min(iprot))
291 write (iout,*) "Revised energy cutoffs"
293 write (iout,*) "Protein",iprot,":",enecut(iprot)
298 c Send the order to the second master to revise the list of conformations.
303 call make_list(.true.,.false.,n,x)
305 call minimize_init(n,iv,v,d)
308 write (iout,*) "Diminishing IV(17) to",iv(17)
309 write (iout,*) "Diminishing IV(18) to",iv(18)
311 else if (iv(1).eq.11 .and. WhatsUp.eq.1) then
312 write (iout,*) "Time up, terminating."
317 call targetfun(n,x,nf,f,uiparm,fdum)
319 write(iout,*) "Final function value:",f,v(10)
320 write(iout,*) "Final value of the object function:",obf
321 write(iout,*) "Number of target function evaluations:",nfun
322 write(iout,*) "Number of target gradient evaluations:",ngrad
323 write(iout,*) "SUMSL return code:",iretcode
324 write(*,*) "Final function value:",f,v(10)
325 write(*,*) "Final value of the object function:",obf
326 write(*,*) "Number of target function evaluations:",nfun
327 write(*,*) "Number of target gradient evaluations:",ngrad
328 write(*,*) "SUMSL return code:",iretcode
334 open (88,file=restartname,status="unknown")
340 write (iout,*) "Final checkgrad:"
344 write (iout,*) "Unknown minimizer."
348 c-----------------------------------------------------------------------------
349 real function funk(p)
352 include "DIMENSIONS.ZSCOPT"
354 double precision x(max_paropt)
358 double precision ufparm
359 include "COMMON.IOUNITS"
365 write (iout,*) "funk n",n," p",(p(i),i=1,n)
369 c write (iout,*) "funk n",n," x",(x(i),i=1,n)
370 call targetfun(n,x,nf,f,uiparm,fdum)
374 c-----------------------------------------------------------------------------
375 c logical function stopx(idum)
380 c-----------------------------------------------------------------------------
381 subroutine checkgrad(n,x)
384 include "DIMENSIONS.ZSCOPT"
389 integer status(MPI_STATUS_SIZE)
391 include "COMMON.IOUNITS"
392 include "COMMON.VMCPAR"
393 include "COMMON.OPTIM"
394 include "COMMON.TIME1"
395 include "COMMON.ENERGIES"
396 include "COMMON.WEIGHTDER"
397 include "COMMON.WEIGHTS"
398 include "COMMON.CLASSES"
399 include "COMMON.INTERACT"
400 include "COMMON.NAMES"
402 include "COMMON.COMPAR"
403 include "COMMON.TORSION"
404 double precision energia(0:max_ene),epskl,eneps_num
405 double precision f,fplus,fminus,fpot,fpotplus,xi,delta /5.0d-5/,
406 & g(max_paropt),gnum(max_paropt),gfpot(max_paropt),gg(max_paropt),
407 & rdum,zscder,pom,sumlikder1(maxvar,maxT,maxprot)
409 double precision x(n)
410 integer uiparm,maxfun,i,ibatch,iprot,j,nf,k,kk,l,ii,iti,itj,ib,
412 double precision nuave_orig(maxdimnat,maxT,maxnatlike,maxprot)
414 double precision enetb_orig1(maxstr_proc,max_ene,maxprot)
418 double precision rdif
419 double precision urparm,aux,fac,eoldsc
422 & rder_all(max_paropt,maxdimnat,maxT,maxnatlike,maxprot),
423 & zvtot_orig(maxT,maxprot),sumlik_orig(maxT,maxprot)
425 integer ind_shield /25/
427 write (iout,*) "Entered CHECKGRAD"
430 call targetfun(n,x,nf,f,uiparm,fdum)
431 call targetgrad(n,x,nf,g,uiparm,rdum,fdum)
432 write (iout,*) "Initial function & gradient evaluation completed."
441 enetb_orig1(jj,j,iprot)=enetb(jj,j,iprot)
443 c write (iout,*) i,enetb_orig1(jj,1,iprot)
448 write (iout,*) "Derivatives in epsilons"
451 call restore_molinfo(iprot)
456 c write (iout,*) "iprot",iprot," i",i
457 call etotal(energia(0))
462 eps(k,l)=eps(k,l)+delta
464 c write (iout,*) "***",k,l,epskl
467 call restore_ccoords(iprot,i)
468 call etotal(energia(0))
469 eneps_num=(energia(1)-eoldsc)
471 c write (iout,*) energia(1),enetb_orig1(jj,1,iprot)
474 c write (iout,*) "---",k,l,epskl
475 write(iout,'(2a4,3e20.10)') restyp(k),restyp(l),
476 & eneps(kk,i,iprot),eneps_num,
477 & (eneps_num-eneps(kk,i,iprot))*100
489 write (iout,*) "calling func1 0"
493 c AL 2/10/2018 Add PMF gradient
494 if (mod_fourier(nloctyp).gt.0) then
495 fpot = rdif(n,x,gfpot,2)
499 do inat=1,natlike(iprot)
500 do ib=1,nbeta(inat+2,iprot)
501 do i=1,natdim(inat,iprot)
502 nuave_orig(i,ib,inat,iprot)=nuave(i,ib,inat,iprot)
508 do ib=1,nbeta(1,iprot)
509 sumlik_orig(ib,iprot)=sumlik(ib,iprot)
513 do ib=1,nbeta(2,iprot)
514 zvtot_orig(ib,iprot)=zvtot(ib,iprot)
517 write(iout,*) "exit func1 0"
519 do ib=1,nbeta(2,iprot)
521 zvdertot(i,ib,iprot)=0.0d0
528 if (imask(k).gt.0 .and. k.ne.ind_shield) ii=ii+1
531 do ib=1,nbeta(1,iprot)
533 sumlikder1(k,ib,iprot)=sumlikder(k,ib,iprot)
534 c write (iout,*) "ib",ib," k",k," sumlikder",
535 c & sumlikder1(k,ib,iprot)
540 c AL 2/10/2018 Add PMF gradient
541 if (mod_fourier(nloctyp).gt.0 .and. torsion_pmf) then
542 do i=1,nloctyp*(2*nloctyp-1)
544 sumlikder1(ii,ib,iprot)=0.0d0
549 c Maximum likelihood terms
550 do ib=1,nbeta(1,iprot)
555 sumlikder1(kk,ib,iprot)=
556 & sumlikeps(icant(iti,iti),ib,iprot)
558 sumlikder1(kk,ib,iprot)=
559 & sumlikder1(kk,ib,iprot)+
560 & sumlikeps(icant(iti,itj),ib,iprot)
567 sumlikder1(kk,ib,iprot)=
568 & sumlikeps(icant(iti,itj),ib,iprot)
569 c write (iout,*) "kk",kk," iprot",iprot," ib",ib,
570 c & " iti",iti," itj",itj,
571 c & " ind",icant(iti,itj)," sumlikeps",
572 c & sumlikeps(icant(iti,itj),ib,iprot)
576 sumlikder1(kk,ib,iprot)=sumliktor(i,ib,iprot)
580 sumlikder1(kk,ib,iprot)=sumlikang(i,ib,iprot)
583 c Heat capacity terms
584 do ib=1,nbeta(2,iprot)
589 zvdertot(kk,ib,iprot)=
590 & zvepstot(icant(iti,iti),ib,iprot)
592 zvdertot(kk,ib,iprot)=
593 & zvdertot(kk,ib,iprot)+
594 & zvepstot(icant(iti,j),ib,iprot)
601 zvdertot(kk,ib,iprot)=
602 & zvepstot(icant(iti,itj),ib,iprot)
606 zvdertot(kk,ib,iprot)=zvtortot(i,ib,iprot)
610 zvdertot(kk,ib,iprot)=zvangtot(i,ib,iprot)
614 c Molecular properties
616 do inat=1,natlike(iprot)
617 do ib=1,nbeta(inat+2,iprot)
618 call propderiv(ib,inat,iprot)
619 do k=1,natdim(inat,iprot)
621 rder_all(kk,k,ib,inat,iprot)=rder(kk,k)
627 rder_all(kk,k,ib,inat,iprot)=
628 & reps(icant(iti,iti),k)
630 rder_all(kk,k,ib,inat,iprot)=
631 & rder_all(kk,k,ib,inat,iprot)+
632 & reps(icant(iti,j),k)
639 rder_all(kk,k,ib,inat,iprot)=
640 & reps(icant(iti,itj),k)
645 rder_all(kk,k,ib,inat,iprot)=rtor(i,k)
649 rder_all(kk,k,ib,inat,iprot)=rang(i,k)
655 c Calculate numerical derivatives and compare with analytical derivatives
659 c write (iout,*) (x(k),k=1,kk)
661 c write (iout,*) "after func1"
664 c AL 2/10/2018 Add PMF gradient
665 if (mod_fourier(nloctyp).gt.0) then
666 fpotplus = rdif(n,x,gfpot,1)
667 write (iout,*) "after rdif"
671 write (iout,*) "Variable",i
674 write (iout,*) "Derivatives of maximum likelihood"
675 do ib=1,nbeta(1,iprot)
676 c write (iout,*) sumlik_orig(ib,iprot),sumlik(ib,iprot)
677 pom=(sumlik(ib,iprot)-sumlik_orig(ib,iprot))/(xi*delta)
678 write (iout,'(2i5,3e14.5)') iprot,ib,
679 & sumlikder1(i,ib,iprot),pom,(pom-sumlikder1(i,ib,iprot))
680 & /(dabs(sumlikder1(i,ib,iprot))+1.0d-10)
681 c write (iout,*) sumlik(ib,iprot),
682 c & sumlik_orig(ib,iprot)
687 write (iout,*) "Derivatives of heat capacity"
688 do ib=1,nbeta(2,iprot)
689 pom=(zvtot(ib,iprot)-zvtot_orig(ib,iprot))/(xi*delta)
690 write (iout,'(2i5,3e14.5)') iprot,ib,
691 & zvdertot(i,ib,iprot),pom,(pom-zvdertot(i,ib,iprot))
692 & /(dabs(zvdertot(i,ib,iprot))+1.0d-10)
693 c write (iout,*) zvtot(ib,ibatch,iprot),
694 c & zvtot_orig(ib,ibatch,iprot)
697 c Numerical derivatives of natlike properties
698 write (iout,*) "Derivatives of molecular properties"
700 do inat=1,natlike(iprot)
701 do ib=1,nbeta(inat+2,iprot)
702 do k=1,natdim(inat,iprot)
703 rder_num=(nuave(k,ib,inat,iprot)
704 & -nuave_orig(k,ib,inat,iprot))/(xi*delta)
705 write (iout,'(4i3,2e14.5,f5.1)') iprot,inat,ib,k,
706 & rder_all(kk,k,ib,inat,iprot),
707 & rder_num,(rder_num-rder_all(kk,k,ib,inat,iprot))/
708 & (dabs(rder_all(kk,k,ib,inat,iprot))+1.0d-10)*100
714 c AL 2/10/2018 Add PMF gradient
715 if (mod_fourier(nloctyp).gt.0) then
716 write(iout,*)"Derivative of fpot",(-fpotplus+fpot)/(xi*delta),
717 & gfpot(i),(gfpot(i)-(-fpotplus+fpot)/(xi*delta))/gfpot(i)
726 c write (iout,*) "i",i," xplus ",x(i),xi
727 call targetfun(n,x,nf,fplus,uiparm,fdum)
729 c write (iout,*) "i",i," xminus",x(i),xi
730 call targetfun(n,x,nf,fminus,uiparm,fdum)
732 c write (iout,*) "x",i,x(i)," fplus",fplus," fminus",fminus,
733 c & "obfplus",obfplus," obfminus",obfminus
734 gnum(i)=0.5d0*(fplus-fminus)/(xi*delta)
736 write (iout,*) "Comparison of analytical and numerical gradient"
738 write (iout,'(i5,2e15.7,f10.2)') i,g(i),gnum(i),
739 & 1.0d2*(g(i)-gnum(i))/(dabs(g(i))+1.0d-10)
744 c-----------------------------------------------------------------------
745 subroutine prepare_sc
748 include "DIMENSIONS.ZSCOPT"
749 include "COMMON.NAMES"
750 include "COMMON.WEIGHTS"
751 include "COMMON.INTERACT"
752 include "COMMON.ENERGIES"
753 include "COMMON.FFIELD"
754 include "COMMON.TORSION"
755 include "COMMON.IOUNITS"
757 integer i,j,ii,nvarr,iti,itj
759 double precision epsij,rrij,sigeps,sigt1sq,sigt2sq,sigii1,sigii2,
760 & ratsig1,ratsig2,rsum_max,dwa16,r_augm
764 if (lprint) write(iout,'(/3a,2i3)') 'Potential is ',
765 & potname(ipot),', exponents are ',expon,2*expon
766 C-----------------------------------------------------------------------
767 C Calculate the "working" parameters of SC interactions.
768 dwa16=2.0d0**(1.0d0/6.0d0)
771 sigma(i,j)=dsqrt(sigma0(i)**2+sigma0(j)**2)
772 sigma(j,i)=sigma(i,j)
773 rs0(i,j)=dwa16*sigma(i,j)
777 if (lprint) write (iout,'(/a/10x,7a/72(1h-))')
778 & 'Working parameters of the SC interactions:',
779 & ' a ',' b ',' augm ',' sigma ',' r0 ',
784 if (ipot.eq.1 .or. ipot.eq.3 .or. ipot.eq.4) then
793 sigeps=dsign(1.0D0,epsij)
795 aa(i,j)=epsij*rrij*rrij
796 bb(i,j)=-sigeps*epsij*rrij
804 ratsig1=sigt2sq/sigt1sq
805 ratsig2=1.0D0/ratsig1
806 chi(i,j)=(sigii1-1.0D0)/(sigii1+ratsig1)
807 if (j.gt.i) chi(j,i)=(sigii2-1.0D0)/(sigii2+ratsig2)
808 rsum_max=dsqrt(sigii1*sigt1sq+sigii2*sigt2sq)
812 sigmaii(i,j)=rsum_max
813 sigmaii(j,i)=rsum_max
814 if ((ipot.eq.2 .or. ipot.eq.5) .and. r0(i,j).gt.rsum_max)
816 r_augm=sigma(i,j)*(rrij-sigma(i,j))/rrij
817 augm(i,j)=epsij*r_augm**(2*expon)
818 c augm(i,j)=0.5D0**(2*expon)*aa(i,j)
825 write (iout,'(2(a3,2x),3(1pe10.3),5(0pf8.3))')
826 & restyp(i),restyp(j),aa(i,j),bb(i,j),augm(i,j),
827 & sigma(i,j),r0(i,j),chi(i,j),chi(j,i)