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-7/,
406 & g(max_paropt),gnum(max_paropt),gfpot(max_paropt),gg(max_paropt),
407 & gforce(max_paropt),rdum,zscder,pom,fchi,fchiplus,
408 & sumlikder1(maxvar,maxT,maxprot)
409 double precision chisquare_force
411 double precision x(n)
412 integer uiparm,maxfun,i,ibatch,iprot,j,nf,k,kk,l,ii,iti,itj,ib,
414 double precision nuave_orig(maxdimnat,maxT,maxnatlike,maxprot)
416 double precision enetb_orig1(maxstr_proc,max_ene,maxprot)
420 double precision rdif
421 double precision urparm,aux,fac,eoldsc
424 & rder_all(max_paropt,maxdimnat,maxT,maxnatlike,maxprot),
425 & zvtot_orig(maxT,maxprot),sumlik_orig(maxT,maxprot)
427 integer ind_shield /25/
429 c write (iout,*) "Entered CHECKGRAD"
432 call targetfun(n,x,nf,f,uiparm,fdum)
433 c write (iout,*) "f",f
434 call targetgrad(n,x,nf,g,uiparm,rdum,fdum)
435 c write (iout,*) "g",g(:n)
436 write (iout,*) "Initial function & gradient evaluation completed."
445 enetb_orig1(jj,j,iprot)=enetb(jj,j,iprot)
447 c write (iout,*) i,enetb_orig1(jj,1,iprot)
452 write (iout,*) "Derivatives in epsilons"
455 call restore_molinfo(iprot)
460 c write (iout,*) "iprot",iprot," i",i
461 call etotal(energia(0))
466 eps(k,l)=eps(k,l)+delta
468 c write (iout,*) "***",k,l,epskl
471 call restore_ccoords(iprot,i)
472 call etotal(energia(0))
473 eneps_num=(energia(1)-eoldsc)
475 c write (iout,*) energia(1),enetb_orig1(jj,1,iprot)
478 c write (iout,*) "---",k,l,epskl
479 write(iout,'(2a4,3e20.10)') restyp(k),restyp(l),
480 & eneps(kk,i,iprot),eneps_num,
481 & (eneps_num-eneps(kk,i,iprot))*100
493 write (iout,*) "calling func1 0"
497 c AL 2/10/2018 Add PMF gradient
498 if (mod_fourier(nloctyp).gt.0) then
499 fpot = rdif(n,x,gfpot,2)
502 c write (iout,*) "-----------------------------------"
503 c write (iout,*) "Computing force chisquare"
504 fchi = chisquare_force(n,x,gforce,2)
505 c write (iout,*) "After first call to chisquare_force"
506 c write (iout,*) "fchi",fchi
507 c write (iout,*) "gforce",gforce(:n)
508 c write (iout,*) "-----------------------------------"
510 do inat=1,natlike(iprot)
511 do ib=1,nbeta(inat+2,iprot)
512 do i=1,natdim(inat,iprot)
513 nuave_orig(i,ib,inat,iprot)=nuave(i,ib,inat,iprot)
519 do ib=1,nbeta(1,iprot)
520 sumlik_orig(ib,iprot)=sumlik(ib,iprot)
524 do ib=1,nbeta(2,iprot)
525 zvtot_orig(ib,iprot)=zvtot(ib,iprot)
528 c write(iout,*) "exit func1 0"
530 do ib=1,nbeta(2,iprot)
532 zvdertot(i,ib,iprot)=0.0d0
539 if (imask(k).gt.0 .and. k.ne.ind_shield) ii=ii+1
542 do ib=1,nbeta(1,iprot)
544 sumlikder1(k,ib,iprot)=sumlikder(k,ib,iprot)
545 c write (iout,*) "ib",ib," k",k," sumlikder",
546 c & sumlikder1(k,ib,iprot)
551 c AL 2/10/2018 Add PMF gradient
552 if (mod_fourier(nloctyp).gt.0 .and. torsion_pmf) then
553 do i=1,nloctyp*(2*nloctyp-1)
555 sumlikder1(ii,ib,iprot)=0.0d0
561 c Maximum likelihood terms
562 do ib=1,nbeta(1,iprot)
567 sumlikder1(kk,ib,iprot)=
568 & sumlikeps(icant(iti,iti),ib,iprot)
570 sumlikder1(kk,ib,iprot)=
571 & sumlikder1(kk,ib,iprot)+
572 & sumlikeps(icant(iti,itj),ib,iprot)
579 sumlikder1(kk,ib,iprot)=
580 & sumlikeps(icant(iti,itj),ib,iprot)
581 c write (iout,*) "kk",kk," iprot",iprot," ib",ib,
582 c & " iti",iti," itj",itj,
583 c & " ind",icant(iti,itj)," sumlikeps",
584 c & sumlikeps(icant(iti,itj),ib,iprot)
588 sumlikder1(kk,ib,iprot)=sumliktor(i,ib,iprot)
592 sumlikder1(kk,ib,iprot)=sumlikang(i,ib,iprot)
595 c Heat capacity terms
596 do ib=1,nbeta(2,iprot)
601 zvdertot(kk,ib,iprot)=
602 & zvepstot(icant(iti,iti),ib,iprot)
604 zvdertot(kk,ib,iprot)=
605 & zvdertot(kk,ib,iprot)+
606 & zvepstot(icant(iti,j),ib,iprot)
613 zvdertot(kk,ib,iprot)=
614 & zvepstot(icant(iti,itj),ib,iprot)
618 zvdertot(kk,ib,iprot)=zvtortot(i,ib,iprot)
622 zvdertot(kk,ib,iprot)=zvangtot(i,ib,iprot)
626 c Molecular properties
628 do inat=1,natlike(iprot)
629 do ib=1,nbeta(inat+2,iprot)
630 call propderiv(ib,inat,iprot)
631 do k=1,natdim(inat,iprot)
633 rder_all(kk,k,ib,inat,iprot)=rder(kk,k)
639 rder_all(kk,k,ib,inat,iprot)=
640 & reps(icant(iti,iti),k)
642 rder_all(kk,k,ib,inat,iprot)=
643 & rder_all(kk,k,ib,inat,iprot)+
644 & reps(icant(iti,j),k)
651 rder_all(kk,k,ib,inat,iprot)=
652 & reps(icant(iti,itj),k)
657 rder_all(kk,k,ib,inat,iprot)=rtor(i,k)
661 rder_all(kk,k,ib,inat,iprot)=rang(i,k)
667 c Calculate numerical derivatives and compare with analytical derivatives
669 c write (iout,*) "gforce1:",gforce(:n)
672 c write (iout,*) (x(k),k=1,kk)
675 c write (iout,*) "after func1"
678 c AL 2/10/2018 Add PMF gradient
679 if (mod_fourier(nloctyp).gt.0) then
680 fpotplus = rdif(n,x,gfpot,1)
681 c write (iout,*) "after rdif"
685 fchiplus = chisquare_force(n,x,gg,1)
686 c write (iout,*) "gforce2:",gforce(:n)
687 c write (iout,*) "Variable",i," fchiplus",fchiplus
690 write (iout,*) "Derivatives of maximum likelihood"
691 do ib=1,nbeta(1,iprot)
692 c write (iout,*) sumlik_orig(ib,iprot),sumlik(ib,iprot)
693 pom=(sumlik(ib,iprot)-sumlik_orig(ib,iprot))/(xi*delta)
694 write (iout,'(2i5,3e14.5)') iprot,ib,
695 & sumlikder1(i,ib,iprot),pom,(pom-sumlikder1(i,ib,iprot))
696 & /(dabs(sumlikder1(i,ib,iprot))+1.0d-10)
697 c write (iout,*) sumlik(ib,iprot),
698 c & sumlik_orig(ib,iprot)
703 write (iout,*) "Derivatives of heat capacity"
704 do ib=1,nbeta(2,iprot)
705 pom=(zvtot(ib,iprot)-zvtot_orig(ib,iprot))/(xi*delta)
706 write (iout,'(2i5,3e14.5)') iprot,ib,
707 & zvdertot(i,ib,iprot),pom,(pom-zvdertot(i,ib,iprot))
708 & /(dabs(zvdertot(i,ib,iprot))+1.0d-10)
709 c write (iout,*) zvtot(ib,ibatch,iprot),
710 c & zvtot_orig(ib,ibatch,iprot)
713 c Numerical derivatives of natlike properties
714 write (iout,*) "Derivatives of molecular properties"
716 do inat=1,natlike(iprot)
717 do ib=1,nbeta(inat+2,iprot)
718 do k=1,natdim(inat,iprot)
719 rder_num=(nuave(k,ib,inat,iprot)
720 & -nuave_orig(k,ib,inat,iprot))/(xi*delta)
721 write (iout,'(4i3,2e14.5,f5.1)') iprot,inat,ib,k,
722 & rder_all(kk,k,ib,inat,iprot),
723 & rder_num,(rder_num-rder_all(kk,k,ib,inat,iprot))/
724 & (dabs(rder_all(kk,k,ib,inat,iprot))+1.0d-10)*100
730 c AL 2/10/2018 Add PMF gradient
731 if (mod_fourier(nloctyp).gt.0) then
732 write(iout,*)"Derivative of fpot",(-fpotplus+fpot)/(xi*delta),
733 & gfpot(i),(gfpot(i)-(-fpotplus+fpot)/(xi*delta))/gfpot(i)
736 c write (iout,*) "fchiplus",fchiplus," fchi",fchi," xi",xi,
738 write (iout,*)"Derivatives of fchi",(-fchiplus+fchi)/(xi*delta),
739 & gforce(i),(gforce(i)-(-fchiplus+fchi)/(xi*delta))/gforce(i)
746 c write (iout,*) "i",i," xplus ",x(i),xi
747 call targetfun(n,x,nf,fplus,uiparm,fdum)
749 c write (iout,*) "i",i," xminus",x(i),xi
750 call targetfun(n,x,nf,fminus,uiparm,fdum)
752 write (iout,*) "x",i,x(i)," fplus",fplus," fminus",fminus
753 gnum(i)=0.5d0*(fplus-fminus)/(xi*delta)
755 write (iout,*) "Comparison of analytical and numerical gradient"
757 write (iout,'(i5,2e15.7,f10.2)') i,g(i),gnum(i),
758 & 1.0d2*(g(i)-gnum(i))/(dabs(g(i))+1.0d-10)
763 c-----------------------------------------------------------------------
764 subroutine prepare_sc
767 include "DIMENSIONS.ZSCOPT"
768 include "COMMON.NAMES"
769 include "COMMON.WEIGHTS"
770 include "COMMON.INTERACT"
771 include "COMMON.ENERGIES"
772 include "COMMON.FFIELD"
773 include "COMMON.TORSION"
774 include "COMMON.IOUNITS"
776 integer i,j,ii,nvarr,iti,itj
778 double precision epsij,rrij,sigeps,sigt1sq,sigt2sq,sigii1,sigii2,
779 & ratsig1,ratsig2,rsum_max,dwa16,r_augm
783 if (lprint) write(iout,'(/3a,2i3)') 'Potential is ',
784 & potname(ipot),', exponents are ',expon,2*expon
785 C-----------------------------------------------------------------------
786 C Calculate the "working" parameters of SC interactions.
787 dwa16=2.0d0**(1.0d0/6.0d0)
790 sigma(i,j)=dsqrt(sigma0(i)**2+sigma0(j)**2)
791 sigma(j,i)=sigma(i,j)
792 rs0(i,j)=dwa16*sigma(i,j)
796 if (lprint) write (iout,'(/a/10x,7a/72(1h-))')
797 & 'Working parameters of the SC interactions:',
798 & ' a ',' b ',' augm ',' sigma ',' r0 ',
803 if (ipot.eq.1 .or. ipot.eq.3 .or. ipot.eq.4) then
812 sigeps=dsign(1.0D0,epsij)
814 aa(i,j)=epsij*rrij*rrij
815 bb(i,j)=-sigeps*epsij*rrij
823 ratsig1=sigt2sq/sigt1sq
824 ratsig2=1.0D0/ratsig1
825 chi(i,j)=(sigii1-1.0D0)/(sigii1+ratsig1)
826 if (j.gt.i) chi(j,i)=(sigii2-1.0D0)/(sigii2+ratsig2)
827 rsum_max=dsqrt(sigii1*sigt1sq+sigii2*sigt2sq)
831 sigmaii(i,j)=rsum_max
832 sigmaii(j,i)=rsum_max
833 if ((ipot.eq.2 .or. ipot.eq.5) .and. r0(i,j).gt.rsum_max)
835 r_augm=sigma(i,j)*(rrij-sigma(i,j))/rrij
836 augm(i,j)=epsij*r_augm**(2*expon)
837 c augm(i,j)=0.5D0**(2*expon)*aa(i,j)
844 write (iout,'(2(a3,2x),3(1pe10.3),5(0pf8.3))')
845 & restyp(i),restyp(j),aa(i,j),bb(i,j),augm(i,j),
846 & sigma(i,j),r0(i,j),chi(i,j),chi(j,i)