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.XBOUND"
84 include "COMMON.WEIGHTDER"
85 include "COMMON.OPTIM"
86 include "COMMON.CLASSES"
87 include "COMMON.COMPAR"
88 include "COMMON.MINPAR"
89 include "COMMON.IOUNITS"
90 include "COMMON.VMCPAR"
91 include "COMMON.TIME1"
92 include "COMMON.ENERGIES"
93 include "COMMON.PROTNAME"
94 include "COMMON.THERMAL"
100 double precision x(max_paropt),d(max_paropt),v(1:lv),g(max_paropt)
101 real p(max_paropt+1,max_paropt),y(max_paropt+1),pb(max_paropt),
102 & pom(max_paropt),yb,y1,y2,pomi,ftol,temptr
105 external targetfun,targetgrad,fdum
106 integer uiparm,i,ii,j,k,nf,nfun,ngrad,iretcode,ianneal,maxanneal
107 double precision urparm,penalty
108 double precision f,f0,obf,obg(max_paropt),tcpu
109 double precision Evarmax(maxprot),aux
110 integer iscan,iprot,ibatch,niter
114 write (iout,*) "Minimizing with ",MINIMIZER
115 IF (MINIMIZER.EQ."AMEBSA") THEN
131 pom(i-1)=p(i,i-1)-0.2
133 if (pom(i-1).lt.x_low(i-1)) pom(i-1)=x_low(i-1)
135 pom(i-1)=p(i,i-1)+0.2
136 if (pom(i-1).gt.x_up(i-1)) pom(i-1)=x_up(i-1)
147 write (iout,*) "p and y"
149 write (iout,'(20e15.5)') (p(i,j),j=1,n),y(i)
152 do ianneal=1,maxanneal
155 write (iout,*) "ianneal",ianneal," temptr",temptr
156 call amebsa(p,y,max_paropt+1,max_paropt,n,pb,yb,ftol,funk,
158 write (iout,*) "Optimized P (PB)"
160 write (iout,'(i5,f10.5)') i,pb(i)
162 write (iout,*) "Optimal function value",yb
163 write(iout,*) "AMEBSA return code:",iretcode
168 call targetfun(n,x,nf,f,uiparm,fdum)
170 write(iout,*) "Final function value:",f,f
171 write(iout,*) "Final value of the object function:",obf
172 write(iout,*) "Number of target function evaluations:",nfun
178 open (88,file=restartname,status="unknown")
183 ELSE IF (MINIMIZER.EQ."SUMSL") THEN
188 call minimize_init(n,iv,v,d)
190 write (iout,*) "Initial checkgrad:"
192 write (iout,*) "Second checkgrad"
195 write (iout,*) "Calling SUMSL"
197 call targetfun(n,x,nf,f,uiparm,fdum)
198 write (iout,*) "Initial function value",f
201 do while ( iscan.lt.maxscan .and. (f0-f)/(f0+1.0d0).gt.0.2d0 )
203 call scan(n,x,f,.false.)
209 if (maxscan.gt.0) then
211 do while ( iscan.lt.maxscan .and. (f0-f)/(f0+1.0d0).gt.0.2d0 )
214 call sumsl(n,d,x,targetfun,targetgrad,iv,liv,lv,v,uiparm,
216 write (iout,*) "Left SUMSL"
221 call targetfun(n,x,nf,f,uiparm,fdum)
223 call scan(n,x,f,.false.)
225 if (iv(1).eq.11 .and. cutoffviol) then
228 &"Revising the list of conformations because of cut-off violation."
230 write (iout,'(a,f7.1)')
231 & protname(iprot)(:ilen(protname(iprot))),enecut(iprot)
232 do j=1,natlike(iprot)+2
233 write (iout,'(5(f7.1,f9.1,2x))')
234 & (1.0d0/(Rgas*betaT(k,ibatch,iprot)),
235 & Evar(j,k,iprot),k=1,nbeta(j,iprot))
242 call make_list(.true.,.false.,n,x)
245 else if (iv(1).eq.11 .and. WhatsUp.eq.1) then
246 write (iout,*) "Time up, terminating."
254 write (iout,*) "MAXMIN",maxmin
256 & call sumsl(n,d,x,targetfun,targetgrad,iv,liv,lv,v,uiparm,
258 write (iout,*) "Left SUMSL"
264 write (iout,*) "iretcode",iv(1)," cutoffviol",cutoffviol
265 if (iv(1).eq.11 .and. cutoffviol) then
267 &"Revising the list of conformations because of cut-off violation."
269 write (iout,'(a,f7.1)')
270 & protname(iprot)(:ilen(protname(iprot))),enecut(iprot)
271 do j=1,natlike(iprot)+2
272 write (iout,'(5(f7.1,f9.1,2x))')
273 & (1.0d0/(Rgas*betaT(k,ibatch,iprot)),
274 & Evar(j,k,iprot),k=1,nbeta(j,iprot))
279 Evarmax(iprot)=Evar(1,1,iprot)*betaT(1,1,iprot)
280 do ibatch=2,natlike(iprot)+2
281 do k=1,nbeta(ibatch,iprot)
282 aux = Evar(k,ibatch,iprot)*betaT(k,ibatch,iprot)
283 if (aux.gt.Evarmax(iprot)) Evarmax(iprot)=aux
288 aux=dmin1(0.5d0*enecut_min(iprot)+1.25d0*Evarmax(iprot)
289 & ,enecut_max(iprot))
290 enecut(iprot)=dmax1(aux,enecut_min(iprot))
292 write (iout,*) "Revised energy cutoffs"
294 write (iout,*) "Protein",iprot,":",enecut(iprot)
299 c Send the order to the second master to revise the list of conformations.
304 call make_list(.true.,.false.,n,x)
306 call minimize_init(n,iv,v,d)
309 write (iout,*) "Diminishing IV(17) to",iv(17)
310 write (iout,*) "Diminishing IV(18) to",iv(18)
312 else if (iv(1).eq.11 .and. WhatsUp.eq.1) then
313 write (iout,*) "Time up, terminating."
318 call targetfun(n,x,nf,f,uiparm,fdum)
320 write(iout,*) "Final function value:",f,v(10)
321 write(iout,*) "Final value of the object function:",obf
322 write(iout,*) "Number of target function evaluations:",nfun
323 write(iout,*) "Number of target gradient evaluations:",ngrad
324 write(iout,*) "SUMSL return code:",iretcode
325 write(*,*) "Final function value:",f,v(10)
326 write(*,*) "Final value of the object function:",obf
327 write(*,*) "Number of target function evaluations:",nfun
328 write(*,*) "Number of target gradient evaluations:",ngrad
329 write(*,*) "SUMSL return code:",iretcode
335 open (88,file=restartname,status="unknown")
341 write (iout,*) "Final checkgrad:"
345 write (iout,*) "Unknown minimizer."
349 c-----------------------------------------------------------------------------
350 real function funk(p)
353 include "DIMENSIONS.ZSCOPT"
355 double precision x(max_paropt)
359 double precision ufparm
360 include "COMMON.IOUNITS"
366 write (iout,*) "funk n",n," p",(p(i),i=1,n)
370 c write (iout,*) "funk n",n," x",(x(i),i=1,n)
371 call targetfun(n,x,nf,f,uiparm,fdum)
375 c-----------------------------------------------------------------------------
376 c logical function stopx(idum)
381 c-----------------------------------------------------------------------------
382 subroutine checkgrad(n,x)
385 include "DIMENSIONS.ZSCOPT"
390 integer status(MPI_STATUS_SIZE)
392 include "COMMON.IOUNITS"
393 include "COMMON.VMCPAR"
394 include "COMMON.OPTIM"
395 include "COMMON.TIME1"
396 include "COMMON.ENERGIES"
397 include "COMMON.WEIGHTDER"
398 include "COMMON.WEIGHTS"
399 include "COMMON.CLASSES"
400 include "COMMON.INTERACT"
401 include "COMMON.NAMES"
403 include "COMMON.COMPAR"
404 double precision energia(0:max_ene),epskl,eneps_num
405 double precision f,fplus,fminus,xi,delta /5.0d-7/,
406 & g(max_paropt),gnum(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 urparm,aux,fac,eoldsc
421 & rder_all(max_paropt,maxdimnat,maxT,maxnatlike,maxprot),
422 & zvtot_orig(maxT,maxprot),sumlik_orig(maxT,maxprot)
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) 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)
562 c Heat capacity terms
563 do ib=1,nbeta(2,iprot)
568 zvdertot(kk,ib,iprot)=
569 & zvepstot(icant(iti,iti),ib,iprot)
571 zvdertot(kk,ib,iprot)=
572 & zvdertot(kk,ib,iprot)+
573 & zvepstot(icant(iti,j),ib,iprot)
580 zvdertot(kk,ib,iprot)=
581 & zvepstot(icant(iti,itj),ib,iprot)
585 zvdertot(kk,ib,iprot)=zvtortot(i,ib,iprot)
589 c Molecular properties
591 do inat=1,natlike(iprot)
592 do ib=1,nbeta(inat+2,iprot)
593 call propderiv(ib,inat,iprot)
594 do k=1,natdim(inat,iprot)
596 rder_all(kk,k,ib,inat,iprot)=rder(kk,k)
602 rder_all(kk,k,ib,inat,iprot)=
603 & reps(icant(iti,iti),k)
605 rder_all(kk,k,ib,inat,iprot)=
606 & rder_all(kk,k,ib,inat,iprot)+
607 & reps(icant(iti,j),k)
614 rder_all(kk,k,ib,inat,iprot)=
615 & reps(icant(iti,itj),k)
620 rder_all(kk,k,ib,inat,iprot)=rtor(i,k)
626 c Calculate numerical derivatives and compare with analytical derivatives
630 c write (iout,*) (x(k),k=1,kk)
632 write (iout,*) "Variable",i
635 write (iout,*) "Derivatives of maximum likelihood"
636 do ib=1,nbeta(1,iprot)
637 c write (iout,*) sumlik_orig(ib,iprot),sumlik(ib,iprot)
638 pom=(sumlik(ib,iprot)-sumlik_orig(ib,iprot))/delta
639 write (iout,'(2i5,3e14.5)') iprot,ib,
640 & sumlikder1(i,ib,iprot),pom,(pom-sumlikder1(i,ib,iprot))
641 & /(dabs(sumlikder1(i,ib,iprot))+1.0d-10)
642 c write (iout,*) sumlik(ib,iprot),
643 c & sumlik_orig(ib,iprot)
648 write (iout,*) "Derivatives of heat capacity"
649 do ib=1,nbeta(2,iprot)
650 pom=(zvtot(ib,iprot)-zvtot_orig(ib,iprot))/delta
651 write (iout,'(2i5,3e14.5)') iprot,ib,
652 & zvdertot(i,ib,iprot),pom,(pom-zvdertot(i,ib,iprot))
653 & /(dabs(zvdertot(i,ib,iprot))+1.0d-10)
654 c write (iout,*) zvtot(ib,ibatch,iprot),
655 c & zvtot_orig(ib,ibatch,iprot)
658 c Numerical derivatives of natlike properties
659 write (iout,*) "Derivatives of molecular properties"
661 do inat=1,natlike(iprot)
662 do ib=1,nbeta(inat+2,iprot)
663 do k=1,natdim(inat,iprot)
664 rder_num=(nuave(k,ib,inat,iprot)
665 & -nuave_orig(k,ib,inat,iprot))/delta
666 write (iout,'(4i3,2e14.5,f5.1)') iprot,inat,ib,k,
667 & rder_all(kk,k,ib,inat,iprot),
668 & rder_num,(rder_num-rder_all(kk,k,ib,inat,iprot))/
669 & (dabs(rder_all(kk,k,ib,inat,iprot))+1.0d-10)*100
680 c write (iout,*) "i",i," xplus ",x(i),xi
681 call targetfun(n,x,nf,fplus,uiparm,fdum)
683 c write (iout,*) "i",i," xminus",x(i),xi
684 call targetfun(n,x,nf,fminus,uiparm,fdum)
686 c write (iout,*) "x",i,x(i)," fplus",fplus," fminus",fminus,
687 c & "obfplus",obfplus," obfminus",obfminus
688 gnum(i)=0.5d0*(fplus-fminus)/delta
690 write (iout,*) "Comparison of analytical and numerical gradient"
692 write (iout,'(i5,2e15.7,f10.2)') i,g(i),gnum(i),
693 & 1.0d2*(g(i)-gnum(i))/(dabs(g(i))+1.0d-10)
698 c-----------------------------------------------------------------------
699 subroutine prepare_sc
702 include "DIMENSIONS.ZSCOPT"
703 include "COMMON.NAMES"
704 include "COMMON.WEIGHTS"
705 include "COMMON.INTERACT"
706 include "COMMON.ENERGIES"
707 include "COMMON.FFIELD"
708 include "COMMON.TORSION"
709 include "COMMON.IOUNITS"
711 integer i,j,ii,nvarr,iti,itj
713 double precision epsij,rrij,sigeps,sigt1sq,sigt2sq,sigii1,sigii2,
714 & ratsig1,ratsig2,rsum_max,dwa16,r_augm
718 if (lprint) write(iout,'(/3a,2i3)') 'Potential is ',
719 & potname(ipot),', exponents are ',expon,2*expon
720 C-----------------------------------------------------------------------
721 C Calculate the "working" parameters of SC interactions.
722 dwa16=2.0d0**(1.0d0/6.0d0)
725 sigma(i,j)=dsqrt(sigma0(i)**2+sigma0(j)**2)
726 sigma(j,i)=sigma(i,j)
727 rs0(i,j)=dwa16*sigma(i,j)
731 if (lprint) write (iout,'(/a/10x,7a/72(1h-))')
732 & 'Working parameters of the SC interactions:',
733 & ' a ',' b ',' augm ',' sigma ',' r0 ',
738 if (ipot.eq.1 .or. ipot.eq.3 .or. ipot.eq.4) then
747 sigeps=dsign(1.0D0,epsij)
749 aa(i,j)=epsij*rrij*rrij
750 bb(i,j)=-sigeps*epsij*rrij
758 ratsig1=sigt2sq/sigt1sq
759 ratsig2=1.0D0/ratsig1
760 chi(i,j)=(sigii1-1.0D0)/(sigii1+ratsig1)
761 if (j.gt.i) chi(j,i)=(sigii2-1.0D0)/(sigii2+ratsig2)
762 rsum_max=dsqrt(sigii1*sigt1sq+sigii2*sigt2sq)
766 sigmaii(i,j)=rsum_max
767 sigmaii(j,i)=rsum_max
768 if ((ipot.eq.2 .or. ipot.eq.5) .and. r0(i,j).gt.rsum_max)
770 r_augm=sigma(i,j)*(rrij-sigma(i,j))/rrij
771 augm(i,j)=epsij*r_augm**(2*expon)
772 c augm(i,j)=0.5D0**(2*expon)*aa(i,j)
779 write (iout,'(2(a3,2x),3(1pe10.3),5(0pf8.3))')
780 & restyp(i),restyp(j),aa(i,j),bb(i,j),augm(i,j),
781 & sigma(i,j),r0(i,j),chi(i,j),chi(j,i)