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 double precision energia(0:max_ene),epskl,eneps_num
404 double precision f,fplus,fminus,xi,delta /5.0d-5/,
405 & g(max_paropt),gnum(max_paropt),
406 & rdum,zscder,pom,sumlikder1(maxvar,maxT,maxprot)
408 double precision x(n)
409 integer uiparm,maxfun,i,ibatch,iprot,j,nf,k,kk,l,ii,iti,itj,ib,
411 double precision nuave_orig(maxdimnat,maxT,maxnatlike,maxprot)
413 double precision enetb_orig1(maxstr_proc,max_ene,maxprot)
417 double precision urparm,aux,fac,eoldsc
420 & rder_all(max_paropt,maxdimnat,maxT,maxnatlike,maxprot),
421 & zvtot_orig(maxT,maxprot),sumlik_orig(maxT,maxprot)
423 integer ind_shield /25/
425 write (iout,*) "Entered CHECKGRAD"
428 call targetfun(n,x,nf,f,uiparm,fdum)
429 call targetgrad(n,x,nf,g,uiparm,rdum,fdum)
430 write (iout,*) "Initial function & gradient evaluation completed."
439 enetb_orig1(jj,j,iprot)=enetb(jj,j,iprot)
441 c write (iout,*) i,enetb_orig1(jj,1,iprot)
446 write (iout,*) "Derivatives in epsilons"
449 call restore_molinfo(iprot)
454 c write (iout,*) "iprot",iprot," i",i
455 call etotal(energia(0))
460 eps(k,l)=eps(k,l)+delta
462 c write (iout,*) "***",k,l,epskl
465 call restore_ccoords(iprot,i)
466 call etotal(energia(0))
467 eneps_num=(energia(1)-eoldsc)
469 c write (iout,*) energia(1),enetb_orig1(jj,1,iprot)
472 c write (iout,*) "---",k,l,epskl
473 write(iout,'(2a4,3e20.10)') restyp(k),restyp(l),
474 & eneps(kk,i,iprot),eneps_num,
475 & (eneps_num-eneps(kk,i,iprot))*100
487 write (iout,*) "calling func1 0"
491 do inat=1,natlike(iprot)
492 do ib=1,nbeta(inat+2,iprot)
493 do i=1,natdim(inat,iprot)
494 nuave_orig(i,ib,inat,iprot)=nuave(i,ib,inat,iprot)
500 do ib=1,nbeta(1,iprot)
501 sumlik_orig(ib,iprot)=sumlik(ib,iprot)
505 do ib=1,nbeta(2,iprot)
506 zvtot_orig(ib,iprot)=zvtot(ib,iprot)
509 write(iout,*) "exit func1 0"
511 do ib=1,nbeta(2,iprot)
513 zvdertot(i,ib,iprot)=0.0d0
520 if (imask(k).gt.0 .and. k.ne.ind_shield) ii=ii+1
523 do ib=1,nbeta(1,iprot)
525 sumlikder1(k,ib,iprot)=sumlikder(k,ib,iprot)
526 c write (iout,*) "ib",ib," k",k," sumlikder",
527 c & sumlikder1(k,ib,iprot)
532 c Maximum likelihood terms
533 do ib=1,nbeta(1,iprot)
538 sumlikder1(kk,ib,iprot)=
539 & sumlikeps(icant(iti,iti),ib,iprot)
541 sumlikder1(kk,ib,iprot)=
542 & sumlikder1(kk,ib,iprot)+
543 & sumlikeps(icant(iti,itj),ib,iprot)
550 sumlikder1(kk,ib,iprot)=
551 & sumlikeps(icant(iti,itj),ib,iprot)
552 c write (iout,*) "kk",kk," iprot",iprot," ib",ib,
553 c & " iti",iti," itj",itj,
554 c & " ind",icant(iti,itj)," sumlikeps",
555 c & sumlikeps(icant(iti,itj),ib,iprot)
559 sumlikder1(kk,ib,iprot)=sumliktor(i,ib,iprot)
563 sumlikder1(kk,ib,iprot)=sumlikang(i,ib,iprot)
566 c Heat capacity terms
567 do ib=1,nbeta(2,iprot)
572 zvdertot(kk,ib,iprot)=
573 & zvepstot(icant(iti,iti),ib,iprot)
575 zvdertot(kk,ib,iprot)=
576 & zvdertot(kk,ib,iprot)+
577 & zvepstot(icant(iti,j),ib,iprot)
584 zvdertot(kk,ib,iprot)=
585 & zvepstot(icant(iti,itj),ib,iprot)
589 zvdertot(kk,ib,iprot)=zvtortot(i,ib,iprot)
593 zvdertot(kk,ib,iprot)=zvangtot(i,ib,iprot)
597 c Molecular properties
599 do inat=1,natlike(iprot)
600 do ib=1,nbeta(inat+2,iprot)
601 call propderiv(ib,inat,iprot)
602 do k=1,natdim(inat,iprot)
604 rder_all(kk,k,ib,inat,iprot)=rder(kk,k)
610 rder_all(kk,k,ib,inat,iprot)=
611 & reps(icant(iti,iti),k)
613 rder_all(kk,k,ib,inat,iprot)=
614 & rder_all(kk,k,ib,inat,iprot)+
615 & reps(icant(iti,j),k)
622 rder_all(kk,k,ib,inat,iprot)=
623 & reps(icant(iti,itj),k)
628 rder_all(kk,k,ib,inat,iprot)=rtor(i,k)
632 rder_all(kk,k,ib,inat,iprot)=rang(i,k)
638 c Calculate numerical derivatives and compare with analytical derivatives
642 c write (iout,*) (x(k),k=1,kk)
644 write (iout,*) "Variable",i
647 write (iout,*) "Derivatives of maximum likelihood"
648 do ib=1,nbeta(1,iprot)
649 c write (iout,*) sumlik_orig(ib,iprot),sumlik(ib,iprot)
650 pom=(sumlik(ib,iprot)-sumlik_orig(ib,iprot))/(xi*delta)
651 write (iout,'(2i5,3e14.5)') iprot,ib,
652 & sumlikder1(i,ib,iprot),pom,(pom-sumlikder1(i,ib,iprot))
653 & /(dabs(sumlikder1(i,ib,iprot))+1.0d-10)
654 c write (iout,*) sumlik(ib,iprot),
655 c & sumlik_orig(ib,iprot)
660 write (iout,*) "Derivatives of heat capacity"
661 do ib=1,nbeta(2,iprot)
662 pom=(zvtot(ib,iprot)-zvtot_orig(ib,iprot))/(xi*delta)
663 write (iout,'(2i5,3e14.5)') iprot,ib,
664 & zvdertot(i,ib,iprot),pom,(pom-zvdertot(i,ib,iprot))
665 & /(dabs(zvdertot(i,ib,iprot))+1.0d-10)
666 c write (iout,*) zvtot(ib,ibatch,iprot),
667 c & zvtot_orig(ib,ibatch,iprot)
670 c Numerical derivatives of natlike properties
671 write (iout,*) "Derivatives of molecular properties"
673 do inat=1,natlike(iprot)
674 do ib=1,nbeta(inat+2,iprot)
675 do k=1,natdim(inat,iprot)
676 rder_num=(nuave(k,ib,inat,iprot)
677 & -nuave_orig(k,ib,inat,iprot))/(xi*delta)
678 write (iout,'(4i3,2e14.5,f5.1)') iprot,inat,ib,k,
679 & rder_all(kk,k,ib,inat,iprot),
680 & rder_num,(rder_num-rder_all(kk,k,ib,inat,iprot))/
681 & (dabs(rder_all(kk,k,ib,inat,iprot))+1.0d-10)*100
692 c write (iout,*) "i",i," xplus ",x(i),xi
693 call targetfun(n,x,nf,fplus,uiparm,fdum)
695 c write (iout,*) "i",i," xminus",x(i),xi
696 call targetfun(n,x,nf,fminus,uiparm,fdum)
698 c write (iout,*) "x",i,x(i)," fplus",fplus," fminus",fminus,
699 c & "obfplus",obfplus," obfminus",obfminus
700 gnum(i)=0.5d0*(fplus-fminus)/(xi*delta)
702 write (iout,*) "Comparison of analytical and numerical gradient"
704 write (iout,'(i5,2e15.7,f10.2)') i,g(i),gnum(i),
705 & 1.0d2*(g(i)-gnum(i))/(dabs(g(i))+1.0d-10)
710 c-----------------------------------------------------------------------
711 subroutine prepare_sc
714 include "DIMENSIONS.ZSCOPT"
715 include "COMMON.NAMES"
716 include "COMMON.WEIGHTS"
717 include "COMMON.INTERACT"
718 include "COMMON.ENERGIES"
719 include "COMMON.FFIELD"
720 include "COMMON.TORSION"
721 include "COMMON.IOUNITS"
723 integer i,j,ii,nvarr,iti,itj
725 double precision epsij,rrij,sigeps,sigt1sq,sigt2sq,sigii1,sigii2,
726 & ratsig1,ratsig2,rsum_max,dwa16,r_augm
730 if (lprint) write(iout,'(/3a,2i3)') 'Potential is ',
731 & potname(ipot),', exponents are ',expon,2*expon
732 C-----------------------------------------------------------------------
733 C Calculate the "working" parameters of SC interactions.
734 dwa16=2.0d0**(1.0d0/6.0d0)
737 sigma(i,j)=dsqrt(sigma0(i)**2+sigma0(j)**2)
738 sigma(j,i)=sigma(i,j)
739 rs0(i,j)=dwa16*sigma(i,j)
743 if (lprint) write (iout,'(/a/10x,7a/72(1h-))')
744 & 'Working parameters of the SC interactions:',
745 & ' a ',' b ',' augm ',' sigma ',' r0 ',
750 if (ipot.eq.1 .or. ipot.eq.3 .or. ipot.eq.4) then
759 sigeps=dsign(1.0D0,epsij)
761 aa(i,j)=epsij*rrij*rrij
762 bb(i,j)=-sigeps*epsij*rrij
770 ratsig1=sigt2sq/sigt1sq
771 ratsig2=1.0D0/ratsig1
772 chi(i,j)=(sigii1-1.0D0)/(sigii1+ratsig1)
773 if (j.gt.i) chi(j,i)=(sigii2-1.0D0)/(sigii2+ratsig2)
774 rsum_max=dsqrt(sigii1*sigt1sq+sigii2*sigt2sq)
778 sigmaii(i,j)=rsum_max
779 sigmaii(j,i)=rsum_max
780 if ((ipot.eq.2 .or. ipot.eq.5) .and. r0(i,j).gt.rsum_max)
782 r_augm=sigma(i,j)*(rrij-sigma(i,j))/rrij
783 augm(i,j)=epsij*r_augm**(2*expon)
784 c augm(i,j)=0.5D0**(2*expon)*aa(i,j)
791 write (iout,'(2(a3,2x),3(1pe10.3),5(0pf8.3))')
792 & restyp(i),restyp(j),aa(i,j),bb(i,j),augm(i,j),
793 & sigma(i,j),r0(i,j),chi(i,j),chi(j,i)