1 subroutine WHAM_CALC(islice,*)
2 ! Weighed Histogram Analysis Method (WHAM) code
3 ! Written by A. Liwo based on the work of Kumar et al.,
4 ! J.Comput.Chem., 13, 1011 (1992)
6 ! 2/1/05 Multiple temperatures allowed.
7 ! 2/2/05 Free energies calculated directly from data points
8 ! acc. to Eq. (21) of Kumar et al.; final histograms also
9 ! constructed based on this equation.
10 ! 2/12/05 Multiple parameter sets included
12 ! 2/2/05 Parallel version
15 include "DIMENSIONS.ZSCOPT"
16 include "DIMENSIONS.FREE"
18 parameter (NGridT=400)
19 integer MaxBinRms,MaxBinRgy
20 parameter (MaxBinRms=100,MaxBinRgy=100)
22 c parameter (MaxHdim=200)
24 parameter (maxinde=200)
28 integer ierror,errcode,status(MPI_STATUS_SIZE)
30 include "COMMON.CONTROL"
31 include "COMMON.IOUNITS"
33 include "COMMON.ENERGIES"
34 include "COMMON.HOMOLOGY"
35 include "COMMON.FFIELD"
36 include "COMMON.SBRIDGE"
38 include "COMMON.ENEPS"
39 include "COMMON.SHIELD"
40 integer MaxPoint,MaxPointProc
41 parameter (MaxPoint=MaxStr,
42 & MaxPointProc=MaxStr_Proc)
43 double precision finorm_max,potfac,entmin,entmax,expfac,vf
44 parameter (finorm_max=1.0d0)
46 integer i,ii,j,jj,k,kk,l,m,ind,iter,t,tmax,ient,ientmax,iln
47 integer start,end,iharm,ib,iib,nbin1,nbin,nbin_rms,nbin_rgy,
48 & nbin_rmsrgy,liczba,iparm,nFi,indrgy,indrms
49 integer htot(0:MaxHdim),histent(0:2000)
50 double precision v(MaxPointProc,MaxR,MaxT_h,Max_Parm)
51 double precision energia(0:max_ene)
53 integer tmax_t,upindE_p
54 double precision fi_p(MaxR,MaxT_h,Max_Parm),
55 & fimax_p(MaxR,MaxT_h,Max_Parm)
56 double precision sumW_p(0:nGridT,Max_Parm),
57 & sumE_p(0:nGridT,Max_Parm),sumEsq_p(0:nGridT,Max_Parm),
58 & sumQ_p(MaxQ1,0:nGridT,Max_Parm),
59 & sumQsq_p(MaxQ1,0:nGridT,Max_Parm),
60 & sumEQ_p(MaxQ1,0:nGridT,Max_Parm),
61 & sumEprim_p(MaxQ1,0:nGridT,Max_Parm),
62 & sumEbis_p(0:nGridT,Max_Parm)
63 double precision hfin_p(0:MaxHdim,maxT_h),
64 & hfin_ent_p(0:MaxHdim),histE_p(0:maxindE),sumH,
65 & hrmsrgy_p(0:MaxBinRgy,0:MaxBinRms,maxT_h)
66 double precision rgymin_t,rmsmin_t,rgymax_t,rmsmax_t
67 double precision potEmin_t,entmin_p,entmax_p
68 double precision ePMF,ePMF_q
69 double precision weimax_(0:ngridT)
70 integer histent_p(0:2000)
71 logical lprint /.true./
73 double precision delta_T /1.0d0/
74 double precision rgymin,rmsmin,rgymax,rmsmax
75 double precision sumW(0:NGridT,Max_Parm),sumE(0:NGridT,Max_Parm),
76 & sumEsq(0:NGridT,Max_Parm),sumQ(MaxQ1,0:NGridT,Max_Parm),
77 & sumQsq(MaxQ1,0:NGridT,Max_Parm),sumEQ(MaxQ1,0:NGridT,Max_Parm),
78 & sumEprim(0:NGridT,Max_Parm),sumEbis(0:NGridT,Max_Parm),betaT,
80 double precision fi(MaxR,maxT_h,Max_Parm),
81 & fimax(MaxR,maxT_h,Max_Parm),
82 & dd,dd1,dd2,hh,dmin,denom,finorm,avefi,pom,
83 & hfin(0:MaxHdim,maxT_h),histE(0:maxindE),
84 & hrmsrgy(0:MaxBinRgy,0:MaxBinRms,maxT_h),
86 & hfin_ent(0:MaxHdim),vmax,aux,weimax(0:nGridT,Max_Parm)
87 double precision fT(6),fTprim(6),fTbis(6),quot,quotl1,quotl,kfacl,
88 & eprim,ebis,temper,kfac/2.4d0/,T0/300.0d0/,startGridT/200.0d0/,
89 & eplus,eminus,logfac,tanhT,tt
90 double precision etot,evdw,evdw_t,evdw2,ees,evdw1,ebe,etors,
91 & escloc,ehpb,ecorr,ecorr5,ecorr6,eello_turn4,eello_turn3,
92 & eturn6,eel_loc,edihcnstr,etors_d,estr,evdw2_14,esccor,
94 & ehomology_constr,edfadis,edfator,edfanei,edfabet
95 integer ind_point(maxpoint),upindE,indE
104 write(licz2,'(bz,i2.2)') islice
106 write (iout,'(//80(1h-)/"Solving WHAM equations for slice",
107 & i2/80(1h-)//)') islice
108 write (iout,*) "delta",delta," nbin1",nbin1
109 write (iout,*) "MaxN",MaxN," MaxQ",MaxQ," MaHdim",MaxHdim
127 c write (iout,*) "i",i," potE",(potE(i,j),j=1,nParmset)
129 if (potE(i,j).le.potEmin) potEmin=potE(i,j)
131 if (q(nQ+1,i).lt.rmsmin) rmsmin=q(nQ+1,i)
132 if (q(nQ+1,i).gt.rmsmax) rmsmax=q(nQ+1,i)
133 if (q(nQ+2,i).lt.rgymin) rgymin=q(nQ+2,i)
134 if (q(nQ+2,i).gt.rgymax) rgymax=q(nQ+2,i)
137 ind=(q(j,i)-dmin+1.0d-8)/delta
139 ind_point(i)=ind_point(i)+ind
141 ind_point(i)=ind_point(i)+nbin1**(j-1)*ind
143 c write (iout,*) "i",i," j",j," q",q(j,i)," ind_point",
146 if (ind_point(i).lt.0 .or. ind_point(i).gt.MaxHdim) then
147 write (iout,*) "Error - index exceeds range for point",i,
148 & " q=",q(j,i)," ind",ind_point(i)
150 write (iout,*) "Processor",me1
152 call MPI_Abort(MPI_COMM_WORLD, Ierror, Errcode )
157 if (ind_point(i).gt.tmax) tmax=ind_point(i)
158 htot(ind_point(i))=htot(ind_point(i))+1
160 write (iout,*) "i",i,"q",(q(j,i),j=1,nQ)," ind",ind_point(i),
161 & " htot",htot(ind_point(i))
166 write (iout,*) "potEmin before reduce",potEmin
168 write (iout,'(a)') "Numbers of counts in Q bins"
170 if (htot(t).gt.0) then
171 write (iout,'(i15,$)') t
174 jj = mod(liczba,nbin1)
176 write (iout,'(i5,$)') jj
178 write (iout,'(i8)') htot(t)
182 write (iout,'(a,i3)') "Number of data points for parameter set",
184 write (iout,'(i7,$)') ((snk(m,ib,iparm,islice),m=1,nr(ib,iparm)),
186 write (iout,'(i8)') stot(islice)
192 call MPI_AllReduce(tmax,tmax_t,1,MPI_INTEGER,MPI_MAX,
195 call MPI_AllReduce(potEmin,potEmin_t,1,MPI_DOUBLE_PRECISION,
196 & MPI_MIN,WHAM_COMM,IERROR)
197 call MPI_AllReduce(rmsmin,rmsmin_t,1,MPI_DOUBLE_PRECISION,
198 & MPI_MIN,WHAM_COMM,IERROR)
199 call MPI_AllReduce(rmsmax,rmsmax_t,1,MPI_DOUBLE_PRECISION,
200 & MPI_MAX,WHAM_COMM,IERROR)
201 call MPI_AllReduce(rgymin,rgymin_t,1,MPI_DOUBLE_PRECISION,
202 & MPI_MIN,WHAM_COMM,IERROR)
203 call MPI_AllReduce(rgymax,rgymax_t,1,MPI_DOUBLE_PRECISION,
204 & MPI_MAX,WHAM_COMM,IERROR)
205 c potEmin=potEmin_t/2
211 write (iout,*) "potEmin",potEmin
213 rmsmin=deltrms*dint(rmsmin/deltrms)
214 rmsmax=deltrms*dint(rmsmax/deltrms)
215 rgymin=deltrms*dint(rgymin/deltrgy)
216 rgymax=deltrms*dint(rgymax/deltrgy)
217 nbin_rms=(rmsmax-rmsmin)/deltrms
218 nbin_rgy=(rgymax-rgymin)/deltrgy
219 write (iout,*) "rmsmin",rmsmin," rmsmax",rmsmax," rgymin",rgymin,
220 & " rgymax",rgymax," nbin_rms",nbin_rms," nbin_rgy",nbin_rgy
227 write (iout,*) "nFi",nFi
228 ! Compute the Boltzmann factor corresponing to restrain potentials in different
235 c write (9,'(3i5,f10.5)') i,(iparm,potE(i,iparm),iparm=1,nParmSet)
238 write (iout,'(2i5,21f8.2)') i,iparm,
239 & (enetb(k,i,iparm),k=1,22)
241 call restore_parm(iparm)
243 write (iout,*) wsc,wscp,welec,wvdwpp,wang,wtor,wscloc,
244 & wcorr,wcorr5,wcorr6,wturn4,wturn3,wturn6,wel_loc,
245 & wtor_d,wsccor,wbond
248 if (rescale_mode.eq.1) then
249 quot=1.0d0/(beta_h(ib,iparm)*1.987D-3*T0)
256 fT(l)=kfacl/(kfacl-1.0d0+quotl)
259 tt = 1.0d0/(beta_h(ib,iparm)*1.987D-3)
260 ft(6)=(320.0d0+80.0d0*dtanh((tt-320.0d0)/80.0d0))/320.0d0
262 ft(6)=1.0d0/(beta_h(ib,iparm)*1.987D-3*T0)
266 else if (rescale_mode.eq.2) then
267 quot=1.0d0/(T0*beta_h(ib,iparm)*1.987D-3)
271 fT(l)=1.12692801104297249644d0/
272 & dlog(dexp(quotl)+dexp(-quotl))
275 tt = 1.0d0/(beta_h(ib,iparm)*1.987D-3)
276 ft(6)=(320.0d0+80.0d0*dtanh((tt-320.0d0)/80.0d0))/320.0d0
278 ft(6)=1.0d0/(beta_h(ib,iparm)*1.987D-3*T0)
282 c write (iout,*) 1.0d0/(beta_h(ib,iparm)*1.987D-3),ft
283 else if (rescale_mode.eq.0) then
288 write (iout,*) "Error in WHAM_CALC: wrong RESCALE_MODE",
293 evdw=enetb(1,i,iparm)
294 evdw_t=enetb(21,i,iparm)
296 evdw2_14=enetb(17,i,iparm)
297 evdw2=enetb(2,i,iparm)+evdw2_14
299 evdw2=enetb(2,i,iparm)
304 evdw1=enetb(16,i,iparm)
309 ecorr=enetb(4,i,iparm)
310 ecorr5=enetb(5,i,iparm)
311 ecorr6=enetb(6,i,iparm)
312 eel_loc=enetb(7,i,iparm)
313 eello_turn3=enetb(8,i,iparm)
314 eello_turn4=enetb(9,i,iparm)
315 eturn6=enetb(10,i,iparm)
316 ebe=enetb(11,i,iparm)
317 escloc=enetb(12,i,iparm)
318 etors=enetb(13,i,iparm)
319 etors_d=enetb(14,i,iparm)
320 ehpb=enetb(15,i,iparm)
321 estr=enetb(18,i,iparm)
322 esccor=enetb(19,i,iparm)
323 edihcnstr=enetb(20,i,iparm)
324 eliptran=enetb(22,i,iparm)
325 esaxs=enetb(26,i,iparm)
326 ehomology_constr=enetb(27,i,iparm)
327 edfadis=enetb(28,i,iparm)
328 edfator=enetb(29,i,iparm)
329 edfanei=enetb(30,i,iparm)
330 edfabet=enetb(31,i,iparm)
333 write (iout,'(3i5,6f5.2,14f12.3)') i,ib,iparm,(ft(l),l=1,6),
334 & evdw+evdw_t,evdw2,ees,evdw1,ecorr,eel_loc,estr,ebe,escloc,
335 & etors,etors_d,eello_turn3,eello_turn4,esccor,esaxs,
336 & ehomology_constr,edfadis,edfator,edfanei,edfabet
340 if (shield_mode.gt.0) then
341 etot=ft(1)*wsc*(evdw+ft(6)*evdw_t)+ft(1)*wscp*evdw2
343 & +ft(1)*wvdwpp*evdw1
344 & +wang*ebe+ft(1)*wtor*etors+wscloc*escloc
345 ! & +wstrain*ehpb+nss*ebr+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5
346 & +wstrain*ehpb+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5
347 & +ft(5)*wcorr6*ecorr6+ft(3)*wturn4*eello_turn4
348 & +ft(2)*wturn3*eello_turn3
349 & +ft(5)*wturn6*eturn6+ft(2)*wel_loc*eel_loc
350 & +edihcnstr+ft(2)*wtor_d*etors_d+ft(1)*wsccor*esccor
351 & +wbond*estr+wliptran*eliptran+wsaxs*esaxs
353 & +wdfa_tor*edfator+wdfa_nei*edfanei+wdfa_beta*edfabet
355 etot=wsc*(evdw+ft(6)*evdw_t)+wscp*evdw2+ft(1)*welec*ees
357 & +wang*ebe+ft(1)*wtor*etors+wscloc*escloc
358 c & +wstrain*ehpb+nss*ebr+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5
359 & +wstrain*ehpb+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5
360 & +ft(5)*wcorr6*ecorr6+ft(3)*wturn4*eello_turn4
361 & +ft(2)*wturn3*eello_turn3
362 & +ft(5)*wturn6*eturn6+ft(2)*wel_loc*eel_loc
363 & +edihcnstr+ft(2)*wtor_d*etors_d+ft(1)*wsccor*esccor
364 & +wbond*estr+wliptran*eliptran+wsaxs*esaxs
366 & +wdfa_tor*edfator+wdfa_nei*edfanei+wdfa_beta*edfabet
369 if (shield_mode.gt.0) then
370 etot=ft(1)*wsc*(evdw+ft(6)*evdw_t)+ft(1)*wscp*evdw2
371 & +ft(1)*welec*(ees+evdw1)
372 & +wang*ebe+ft(1)*wtor*etors+wscloc*escloc
373 c & +wstrain*ehpb+nss*ebr+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5
374 & +wstrain*ehpb+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5
375 & +ft(5)*wcorr6*ecorr6+ft(3)*wturn4*eello_turn4
376 & +ft(2)*wturn3*eello_turn3
377 & +ft(5)*wturn6*eturn6+ft(2)*wel_loc*eel_loc+edihcnstr
378 & +ft(2)*wtor_d*etors_d+ft(1)*wsccor*esccor
379 & +wbond*estr+wliptran*eliptran+wsaxs*esaxs
381 & +wdfa_tor*edfator+wdfa_nei*edfanei+wdfa_beta*edfabet
383 etot=wsc*(evdw+ft(6)*evdw_t)+wscp*evdw2
384 & +ft(1)*welec*(ees+evdw1)
385 & +wang*ebe+ft(1)*wtor*etors+wscloc*escloc
386 c & +wstrain*ehpb+nss*ebr+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5
387 & +wstrain*ehpb+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5
388 & +ft(5)*wcorr6*ecorr6+ft(3)*wturn4*eello_turn4
389 & +ft(2)*wturn3*eello_turn3
390 & +ft(5)*wturn6*eturn6+ft(2)*wel_loc*eel_loc+edihcnstr
391 & +ft(2)*wtor_d*etors_d+ft(1)*wsccor*esccor
392 & +wbond*estr+wliptran*eliptran+wsaxs*esaxs
394 & +wdfa_tor*edfator+wdfa_nei*edfanei+wdfa_beta*edfabet
399 write (iout,*) i,iparm,1.0d0/(beta_h(ib,iparm)*1.987D-3),
403 if (iparm.eq.1 .and. ib.eq.1) then
404 write (iout,*)"Conformation",i
407 energia(k)=enetb(k,i,iparm)
409 call enerprint(energia(0),fT)
413 write (iout,*) "homol_nset",homol_nset,nR(ib,iparm)
415 if (homol_nset.gt.1) then
418 Econstr=waga_homology(kk)*ehomology_constr
420 & -beta_h(ib,iparm)*(etot+Econstr)
422 write (iout,'(4i5,4e15.5)') i,kk,ib,iparm,
423 & etot,Econstr,v(i,kk,ib,iparm)
429 etot=etot+ehomology_constr
435 Econstr=Econstr+Kh(j,kk,ib,iparm)
436 & *(dd-q0(j,kk,ib,iparm))**2
438 c Adaptive potential contribution
440 call PMF_energy(q(1,i),ib,kk,iparm,ePMF,ePMF_q)
444 & -beta_h(ib,iparm)*(etot-potEmin+Econstr)
446 write (iout,'(4i5,4e15.5)') i,kk,ib,iparm,
447 & etot,potEmin,etot-potEmin,v(i,kk,ib,iparm)
456 ! Simple iteration to calculate free energies corresponding to all simulation
460 ! Compute new free-energy values corresponding to the righ-hand side of the
461 ! equation and their derivatives.
462 write (iout,*) "------------------------fi"
472 vf=v(t,l,k,i)+f(l,k,i)
473 if (vf.gt.vmax) vmax=vf
481 aux=f(l,k,i)+v(t,l,k,i)-vmax
483 & denom=denom+snk(l,k,i,islice)*dexp(aux)
487 entfac(t)=-dlog(denom)-vmax
489 write (iout,*) t,"vmax",vmax," denom",denom,"entfac",entfac(t)
495 do ii=1,nR(iib,iparm)
497 fimax_p(ii,iib,iparm)=v(1,ii,iib,iparm)+entfac(1)
499 if(v(t,ii,iib,iparm)+entfac(t).gt.fimax_p(ii,iib,iparm))
500 & fimax_p(ii,iib,iparm)=v(t,ii,iib,iparm)+entfac(t)
503 fimax(ii,iib,iparm)=v(1,ii,iib,iparm)+entfac(1)
505 if(v(t,ii,iib,iparm)+entfac(t).gt.fimax(ii,iib,iparm))
506 & fimax(ii,iib,iparm)=v(t,ii,iib,iparm)+entfac(t)
513 call MPI_AllReduce(fimax_p(1,1,1),fimax(1,1,1),
514 & maxR*MaxT_h*nParmSet,MPI_DOUBLE_PRECISION,
515 & MPI_MAX,WHAM_COMM,IERROR)
519 do ii=1,nR(iib,iparm)
521 fi_p(ii,iib,iparm)=0.0d0
523 fi_p(ii,iib,iparm)=fi_p(ii,iib,iparm)
524 & +dexp(v(t,ii,iib,iparm)+entfac(t)-fimax(ii,iib,iparm))
526 write (iout,'(4i5,3e15.5)') t,ii,iib,iparm,
527 & v(t,ii,iib,iparm),entfac(t),fi_p(ii,iib,iparm)
531 fi(ii,iib,iparm)=0.0d0
533 fi(ii,iib,iparm)=fi(ii,iib,iparm)
534 & +dexp(v(t,ii,iib,iparm)+entfac(t)-fimax(ii,iib,iparm))
543 write (iout,*) "fi before MPI_Reduce me",me,' master',master
545 do ib=1,nT_h(nparmset)
546 write (iout,*) "iparm",iparm," ib",ib
547 write (iout,*) "beta=",beta_h(ib,iparm)
548 write (iout,'(8e15.5)') (fi_p(i,ib,iparm),i=1,nR(ib,iparm))
552 c write (iout,*) "REDUCE size",maxR,MaxT_h,nParmSet,
553 c & maxR*MaxT_h*nParmSet
554 c write (iout,*) "MPI_COMM_WORLD",MPI_COMM_WORLD,
555 c & " WHAM_COMM",WHAM_COMM
556 call MPI_Reduce(fi_p(1,1,1),fi(1,1,1),maxR*MaxT_h*nParmSet,
557 & MPI_DOUBLE_PRECISION,
558 & MPI_SUM,Master,WHAM_COMM,IERROR)
560 write (iout,*) "fi after MPI_Reduce nparmset",nparmset
562 write (iout,*) "iparm",iparm
564 write (iout,*) "beta=",beta_h(ib,iparm)
565 write (iout,'(8e15.5)') (fi(i,ib,iparm),i=1,nR(ib,iparm))
569 if (me1.eq.Master) then
575 fi(i,ib,iparm)=-dlog(fi(i,ib,iparm))-fimax(i,ib,iparm)
576 avefi=avefi+fi(i,ib,iparm)
582 write (iout,*) "Parameter set",iparm
584 write (iout,*) "beta=",beta_h(ib,iparm)
586 fi(i,ib,iparm)=fi(i,ib,iparm)-avefi
588 write (iout,'(8f10.5)') (fi(i,ib,iparm),i=1,nR(ib,iparm))
589 write (iout,'(8f10.5)') (f(i,ib,iparm),i=1,nR(ib,iparm))
593 ! Compute the norm of free-energy increments.
598 finorm=finorm+dabs(fi(i,ib,iparm)-f(i,ib,iparm))
599 f(i,ib,iparm)=fi(i,ib,iparm)
604 write (iout,*) 'Iteration',iter,' finorm',finorm
608 call MPI_Bcast(f(1,1,1),MaxR*MaxT_h*nParmSet,
609 & MPI_DOUBLE_PRECISION,Master,
611 call MPI_Bcast(finorm,1,MPI_DOUBLE_PRECISION,Master,
614 ! Exit, if the increment norm is smaller than pre-assigned tolerance.
615 if (finorm.lt.fimin) then
616 write (iout,*) 'Iteration converged'
623 ! Now, put together the histograms from all simulations, in order to get the
624 ! unbiased total histogram.
634 write (iout,*) "--------------hist"
638 sumW_p(i,iparm)=0.0d0
639 sumE_p(i,iparm)=0.0d0
640 sumEbis_p(i,iparm)=0.0d0
641 sumEsq_p(i,iparm)=0.0d0
643 sumQ_p(j,i,iparm)=0.0d0
644 sumQsq_p(j,i,iparm)=0.0d0
645 sumEQ_p(j,i,iparm)=0.0d0
655 sumEbis(i,iparm)=0.0d0
656 sumEsq(i,iparm)=0.0d0
658 sumQ(j,i,iparm)=0.0d0
659 sumQsq(j,i,iparm)=0.0d0
660 sumEQ(j,i,iparm)=0.0d0
666 c 8/26/05 entropy distribution
671 c ent=-dlog(entfac(t))
673 if (ent.lt.entmin_p) entmin_p=ent
674 if (ent.gt.entmax_p) entmax_p=ent
676 write (iout,*) "entmin",entmin_p," entmax",entmax_p
678 call MPI_Allreduce(entmin_p,entmin,1,MPI_DOUBLE_PRECISION,MPI_MIN,
680 call MPI_Allreduce(entmax_p,entmax,1,MPI_DOUBLE_PRECISION,MPI_MAX,
682 ientmax=entmax-entmin
683 if (ientmax.gt.2000) ientmax=2000
684 write (iout,*) "entmin",entmin," entmax",entmax," ientmax",ientmax
687 c ient=-dlog(entfac(t))-entmin
688 ient=entfac(t)-entmin
689 if (ient.le.2000) histent_p(ient)=histent_p(ient)+1
691 call MPI_Allreduce(histent_p(0),histent(0),ientmax+1,MPI_INTEGER,
692 & MPI_SUM,WHAM_COMM,IERROR)
693 if (me1.eq.Master) then
694 write (iout,*) "Entropy histogram"
696 write(iout,'(f15.4,i10)') entmin+i,histent(i)
704 if (ent.lt.entmin) entmin=ent
705 if (ent.gt.entmax) entmax=ent
707 ientmax=-dlog(entmax)-entmin
708 if (ientmax.gt.2000) ientmax=2000
710 ient=entfac(t)-entmin
711 if (ient.le.2000) histent(ient)=histent(ient)+1
713 write (iout,*) "Entropy histogram"
715 write(iout,'(2f15.4)') entmin+i,histent(i)
720 call restore_parm(iparm)
746 hrmsrgy(j,i,ib)=0.0d0
748 hrmsrgy_p(j,i,ib)=0.0d0
758 indE = aint(potE(t,iparm)-aint(potEmin))
759 if (indE.ge.0 .and. indE.le.maxinde) then
760 if (indE.gt.upindE_p) upindE_p=indE
761 histE_p(indE)=histE_p(indE)+dexp(-entfac(t))
765 expfac=dexp(-beta_h(ib,iparm)*(etot-potEmin)+entfac(t))
766 hfin_p(ind,ib)=hfin_p(ind,ib)+
767 & dexp(-beta_h(ib,iparm)*(etot-potEmin)+entfac(t))
769 indrgy=dint((q(nQ+2,t)-rgymin)/deltrgy)
770 indrms=dint((q(nQ+1,t)-rmsmin)/deltrms)
771 hrmsrgy_p(indrgy,indrms,ib)=
772 & hrmsrgy_p(indrgy,indrms,ib)+expfac
777 expfac=dexp(-beta_h(ib,iparm)*(etot-potEmin)+entfac(t))
778 hfin(ind,ib)=hfin(ind,ib)+
779 & dexp(-beta_h(ib,iparm)*(etot-potEmin)+entfac(t))
781 indrgy=dint((q(nQ+2,t)-rgymin)/deltrgy)
782 indrms=dint((q(nQ+1,t)-rmsmin)/deltrms)
783 hrmsrgy(indrgy,indrms,ib)=
784 & hrmsrgy(indrgy,indrms,ib)+expfac
790 c Thermo and ensemble averages
793 betaT=startGridT+k*delta_T
794 call temp_scalfac(betaT,ft,ftprim,ftbis,*10)
795 c write (iout,*) "ftprim",ftprim
796 c write (iout,*) "ftbis",ftbis
797 betaT=1.0d0/(1.987D-3*betaT)
798 c 7/10/18 AL Determine the max Botzmann weights for each temerature
799 call sum_ene(1,iparm,ft,etot)
800 weimax(k,iparm)=-betaT*(etot-potEmin)+entfac(1)
801 c write (iout,*) "k",k," t",1," weight",weimax(k,iparm)
807 call sum_ene(t,iparm,ft,etot)
808 weight=-betaT*(etot-potEmin)+entfac(t)
809 c write (iout,*) "k",k," t",t," weight",weight
810 if (weight.gt.weimax(k,iparm)) weimax(k,iparm)=weight
815 write (iout,*) "weimax before REDUCE"
816 write (iout,*) (weimax(k,iparm),k=0,ngridt)
819 weimax_(k)=weimax(k,iparm)
821 call MPI_Allreduce(weimax_(0),weimax(0,iparm),nGridT+1,
822 & MPI_DOUBLE_PRECISION,MPI_MAX,WHAM_COMM,IERROR)
824 write (iout,*) "weimax"
825 write (iout,*) (weimax(k,iparm),k=0,ngridt)
828 temper=startGridT+k*delta_T
829 betaT=1.0d0/(1.987D-3*temper)
830 call temp_scalfac(temper,ft,ftprim,ftbis,*10)
837 hfin_ent_p(ind)=hfin_ent_p(ind)+dexp(entfac(t))
839 hfin_ent(ind)=hfin_ent(ind)+dexp(entfac(t))
841 c write (iout,'(2i5,20f8.2)') t,t,(enetb(k,t,iparm),k=1,18)
842 c call restore_parm(iparm)
843 call sum_ene_deriv(t,iparm,ft,ftprim,ftbis,etot,eprim,ebis)
844 weight=dexp(-betaT*(etot-potEmin)+entfac(t)-weimax(k,iparm))
846 write (iout,*) "iparm",iparm," t",t," betaT",betaT,
847 & " etot",etot," entfac",entfac(t)," boltz",
848 & -betaT*(etot-potEmin)+entfac(t)," weimax",weimax(k,iparm),
849 & " weight",weight," ebis",ebis
851 etot=etot-temper*eprim
853 sumW_p(k,iparm)=sumW_p(k,iparm)+weight
854 sumE_p(k,iparm)=sumE_p(k,iparm)+etot*weight
855 sumEbis_p(k,iparm)=sumEbis_p(k,iparm)+ebis*weight
856 sumEsq_p(k,iparm)=sumEsq_p(k,iparm)+etot**2*weight
858 sumQ_p(j,k,iparm)=sumQ_p(j,k,iparm)+q(j,t)*weight
859 sumQsq_p(j,k,iparm)=sumQsq_p(j,k,iparm)+q(j,t)**2*weight
860 sumEQ_p(j,k,iparm)=sumEQ_p(j,k,iparm)
861 & +etot*q(j,t)*weight
864 sumW(k,iparm)=sumW(k,iparm)+weight
865 sumE(k,iparm)=sumE(k,iparm)+etot*weight
866 sumEbis(k,iparm)=sumEbis(k,iparm)+ebis*weight
867 sumEsq(k,iparm)=sumEsq(k,iparm)+etot**2*weight
869 sumQ(j,k,iparm)=sumQ(j,k,iparm)+q(j,t)*weight
870 sumQsq(j,k,iparm)=sumQsq(j,k,iparm)+q(j,t)**2*weight
871 sumEQ(j,k,iparm)=sumEQ(j,k,iparm)
872 & +etot*q(j,t)*weight
878 if (histout) call MPI_Reduce(hfin_p(0,ib),hfin(0,ib),nbin,
879 & MPI_DOUBLE_PRECISION,MPI_SUM,Master,WHAM_COMM,IERROR)
881 call MPI_Reduce(hrmsrgy_p(0,0,ib),hrmsrgy(0,0,ib),
882 & (MaxBinRgy+1)*(nbin_rms+1),MPI_DOUBLE_PRECISION,MPI_SUM,Master,
886 call MPI_Reduce(upindE_p,upindE,1,
887 & MPI_INTEGER,MPI_MAX,Master,WHAM_COMM,IERROR)
888 call MPI_Reduce(histE_p(0),histE(0),maxindE,
889 & MPI_DOUBLE_PRECISION,MPI_SUM,Master,WHAM_COMM,IERROR)
891 if (me1.eq.master) then
895 write (iout,'(6x,$)')
896 write (iout,'(f20.2,$)') (1.0d0/(1.987D-3*beta_h(ib,iparm)),
900 write (iout,'(/a)') 'Final histograms'
902 if (nslice.eq.1) then
903 if (separate_parset) then
904 write(licz3,"(bz,i3.3)") myparm
905 histname=prefix(:ilen(prefix))//'_par'//licz3//'.hist'
907 histname=prefix(:ilen(prefix))//'.hist'
910 if (separate_parset) then
911 write(licz3,"(bz,i3.3)") myparm
912 histname=prefix(:ilen(prefix))//'_par'//licz3//
913 & '_slice_'//licz2//'.hist'
915 histname=prefix(:ilen(prefix))//'_slice_'//licz2//'.hist'
918 #if defined(AIX) || defined(PGI)
919 open (ihist,file=histname,position='append')
921 open (ihist,file=histname,access='append')
931 if (sumH.gt.0.0d0) then
933 jj = mod(liczba,nbin1)
935 write (iout,'(f6.3,$)') dmin+(jj+0.5d0)*delta
937 & write (ihist,'(f6.3,$)') dmin+(jj+0.5d0)*delta
940 write (iout,'(e20.10,$)') hfin(t,ib)
941 if (histfile) write (ihist,'(e20.10,$)') hfin(t,ib)
943 write (iout,'(i5)') iparm
944 if (histfile) write (ihist,'(i5)') iparm
951 if (nslice.eq.1) then
952 if (separate_parset) then
953 write(licz3,"(bz,i3.3)") myparm
954 histname=prefix(:ilen(prefix))//"_par"//licz3//'.ent'
956 histname=prefix(:ilen(prefix))//'.ent'
959 if (separate_parset) then
960 write(licz3,"(bz,i3.3)") myparm
961 histname=prefix(:ilen(prefix))//'par_'//licz3//
962 & '_slice_'//licz2//'.ent'
964 histname=prefix(:ilen(prefix))//'_slice_'//licz2//'.ent'
967 #if defined(AIX) || defined(PGI)
968 open (ihist,file=histname,position='append')
970 open (ihist,file=histname,access='append')
972 write (ihist,'(a)') "# Microcanonical entropy"
974 write (ihist,'(f8.0,$)') dint(potEmin)+i
975 if (histE(i).gt.0.0e0) then
976 write (ihist,'(f15.5,$)') dlog(histE(i))
978 write (ihist,'(f15.5,$)') 0.0d0
984 write (iout,*) "Microcanonical entropy"
986 write (iout,'(f8.0,$)') dint(potEmin)+i
987 if (histE(i).gt.0.0e0) then
988 write (iout,'(f15.5,$)') dlog(histE(i))
990 write (iout,'(f15.5,$)') 0.0d0
995 if (nslice.eq.1) then
996 if (separate_parset) then
997 write(licz3,"(bz,i3.3)") myparm
998 histname=prefix(:ilen(prefix))//'_par'//licz3//'.rmsrgy'
1000 histname=prefix(:ilen(prefix))//'.rmsrgy'
1003 if (separate_parset) then
1004 write(licz3,"(bz,i3.3)") myparm
1005 histname=prefix(:ilen(prefix))//'_par'//licz3//
1006 & '_slice_'//licz2//'.rmsrgy'
1008 histname=prefix(:ilen(prefix))//'_slice_'//licz2//'.rmsrgy'
1011 #if defined(AIX) || defined(PGI)
1012 open (ihist,file=histname,position='append')
1014 open (ihist,file=histname,access='append')
1018 write(ihist,'(2f8.2,$)')
1019 & rgymin+deltrgy*j,rmsmin+deltrms*i
1021 if (hrmsrgy(j,i,ib).gt.0.0d0) then
1022 write(ihist,'(e14.5,$)')
1023 & -dlog(hrmsrgy(j,i,ib))/beta_h(ib,iparm)
1026 write(ihist,'(e14.5,$)') 1.0d6
1029 write (ihist,'(i2)') iparm
1037 call MPI_Reduce(hfin_ent_p(0),hfin_ent(0),nbin,
1038 & MPI_DOUBLE_PRECISION,MPI_SUM,Master,WHAM_COMM,IERROR)
1039 call MPI_Reduce(sumW_p(0,1),sumW(0,1),(nGridT+1)*nParmSet,
1040 & MPI_DOUBLE_PRECISION,MPI_SUM,Master,WHAM_COMM,IERROR)
1041 call MPI_Reduce(sumE_p(0,1),sumE(0,1),(nGridT+1)*nParmSet,
1042 & MPI_DOUBLE_PRECISION,MPI_SUM,Master,WHAM_COMM,IERROR)
1043 call MPI_Reduce(sumEbis_p(0,1),sumEbis(0,1),(nGridT+1)*nParmSet,
1044 & MPI_DOUBLE_PRECISION,MPI_SUM,Master,WHAM_COMM,IERROR)
1045 call MPI_Reduce(sumEsq_p(0,1),sumEsq(0,1),(nGridT+1)*nParmSet,
1046 & MPI_DOUBLE_PRECISION,MPI_SUM,Master,WHAM_COMM,IERROR)
1047 call MPI_Reduce(sumQ_p(1,0,1),sumQ(1,0,1),
1048 & MaxQ1*(nGridT+1)*nParmSet,MPI_DOUBLE_PRECISION,MPI_SUM,Master,
1050 call MPI_Reduce(sumQsq_p(1,0,1),sumQsq(1,0,1),
1051 & MaxQ1*(nGridT+1)*nParmSet,MPI_DOUBLE_PRECISION,MPI_SUM,Master,
1053 call MPI_Reduce(sumEQ_p(1,0,1),sumEQ(1,0,1),
1054 & MaxQ1*(nGridT+1)*nParmSet,MPI_DOUBLE_PRECISION,MPI_SUM,Master,
1056 if (me.eq.master) then
1058 write (iout,'(/a)') 'Thermal characteristics of folding'
1059 if (nslice.eq.1) then
1062 nazwa=prefix(:ilen(prefix))//"_slice_"//licz2
1065 if (nparmset.eq.1 .and. .not.separate_parset) then
1066 nazwa=nazwa(:iln)//".thermal"
1067 else if (nparmset.eq.1 .and. separate_parset) then
1068 write(licz3,"(bz,i3.3)") myparm
1069 nazwa=nazwa(:iln)//"_par_"//licz3//".thermal"
1072 if (nparmset.gt.1) then
1073 write(licz3,"(bz,i3.3)") iparm
1074 nazwa=nazwa(:iln)//"_par_"//licz3//".thermal"
1077 if (separate_parset) then
1078 write (iout,'(a,i3)') "Parameter set",myparm
1080 write (iout,'(a,i3)') "Parameter set",iparm
1083 sumE(i,iparm)=sumE(i,iparm)/sumW(i,iparm)
1084 sumEbis(i,iparm)=(startGridT+i*delta_T)*sumEbis(i,iparm)/
1086 sumEsq(i,iparm)=(sumEsq(i,iparm)/sumW(i,iparm)
1087 & -sumE(i,iparm)**2)/(1.987D-3*(startGridT+i*delta_T)**2)
1089 sumQ(j,i,iparm)=sumQ(j,i,iparm)/sumW(i,iparm)
1090 sumQsq(j,i,iparm)=sumQsq(j,i,iparm)/sumW(i,iparm)
1091 & -sumQ(j,i,iparm)**2
1092 sumEQ(j,i,iparm)=sumEQ(j,i,iparm)/sumW(i,iparm)
1093 & -sumQ(j,i,iparm)*sumE(i,iparm)
1095 sumW(i,iparm)=(-dlog(sumW(i,iparm))-weimax(i,iparm))*(1.987D-3*
1096 & (startGridT+i*delta_T))+potEmin
1097 write (iout,'(f7.1,2f15.5,$)') startGridT+i*delta_T,
1098 & sumW(i,iparm),sumE(i,iparm)
1099 write (iout,'(f10.5,$)') (sumQ(j,i,iparm),j=1,nQ+2)
1100 write (iout,'(e15.5,$)') sumEsq(i,iparm)-sumEbis(i,iparm),
1101 & (sumQsq(j,i,iparm),j=1,nQ+2),(sumEQ(j,i,iparm),j=1,nQ+2)
1103 write (34,'(f7.1,2f15.5,$)') startGridT+i*delta_T,
1104 & sumW(i,iparm),sumE(i,iparm)
1105 write (34,'(f10.5,$)') (sumQ(j,i,iparm),j=1,nQ+2)
1106 write (34,'(e15.5,$)') sumEsq(i,iparm)-sumEbis(i,iparm),
1107 & (sumQsq(j,i,iparm),j=1,nQ+2),(sumEQ(j,i,iparm),j=1,nQ+2)
1115 if (hfin_ent(t).gt.0.0d0) then
1117 jj = mod(liczba,nbin1)
1118 write (iout,'(f6.3,e20.10," ent")') dmin+(jj+0.5d0)*delta,
1120 if (histfile) write (ihist,'(f6.3,e20.10," ent")')
1121 & dmin+(jj+0.5d0)*delta,
1125 if (histfile) close(ihist)
1129 ! Write data for zscore
1130 if (nslice.eq.1) then
1131 zscname=prefix(:ilen(prefix))//".zsc"
1133 zscname=prefix(:ilen(prefix))//"_slice_"//licz2//".zsc"
1135 #if defined(AIX) || defined(PGI)
1136 open (izsc,file=prefix(:ilen(prefix))//'.zsc',position='append')
1138 open (izsc,file=prefix(:ilen(prefix))//'.zsc',access='append')
1140 write (izsc,'("NQ=",i1," NPARM=",i1)') nQ,nParmSet
1142 write (izsc,'("NT=",i1)') nT_h(iparm)
1144 write (izsc,'("TEMP=",f6.1," NR=",i2," SNK=",$)')
1145 & 1.0d0/(beta_h(ib,iparm)*1.987D-3),nR(ib,iparm)
1146 jj = min0(nR(ib,iparm),7)
1147 write (izsc,'(i8,$)') (snk(i,ib,iparm,islice),i=1,jj)
1148 write (izsc,'(a1,$)') (" ",i=22+8*jj+1,79)
1149 write (izsc,'("&")')
1150 if (nR(ib,iparm).gt.7) then
1151 do ii=8,nR(ib,iparm),9
1152 jj = min0(nR(ib,iparm),ii+8)
1153 write (izsc,'(i8,$)') (snk(i,ib,iparm,islice),i=ii,jj)
1154 write (izsc,'(a1,$') (" ",i=(jj-ii+1)*8+1,79)
1155 write (izsc,'("&")')
1158 write (izsc,'("FI=",$)')
1159 jj=min0(nR(ib,iparm),7)
1160 write (izsc,'(f10.5,$)') (fi(i,ib,iparm),i=1,jj)
1161 write (izsc,'(a1,$)') (" ",i=3+10*jj+1,79)
1162 write (izsc,'("&")')
1163 if (nR(ib,iparm).gt.7) then
1164 do ii=8,nR(ib,iparm),9
1165 jj = min0(nR(ib,iparm),ii+8)
1166 write (izsc,'(f10.5,$)') (fi(i,ib,iparm),i=ii,jj)
1167 if (jj.eq.nR(ib,iparm)) then
1170 write (izsc,'(a1,$)') (" ",i=10*(jj-ii+1)+1,79)
1171 write (izsc,'(t80,"&")')
1176 write (izsc,'("KH=",$)')
1177 write (izsc,'(f7.2,$)') (Kh(j,i,ib,iparm),j=1,nQ)
1178 write (izsc,'(" Q0=",$)')
1179 write (izsc,'(f7.5,$)') (q0(j,i,ib,iparm),j=1,nQ)
1194 c------------------------------------------------------------------------
1195 subroutine temp_scalfac(betaT,ft,ftprim,ftbis,*)
1197 include "DIMENSIONS"
1198 include "DIMENSIONS.FREE"
1199 include "COMMON.CONTROL"
1200 include "COMMON.FREE"
1201 include "COMMON.IOUNITS"
1202 double precision fT(6),fTprim(6),fTbis(6),quot,quotl1,quotl,kfacl,
1203 & kfac/2.4d0/,T0/300.0d0/,startGridT/200.0d0/,
1204 & logfac,tanhT,betaT,denom,eplus,eminus
1206 if (rescale_mode.eq.1) then
1214 denom=kfacl-1.0d0+quotl
1216 ftprim(l)=-l*ft(l)*quotl1/(T0*denom)
1217 ftbis(l)=l*kfacl*quotl1*
1218 & (2*l*quotl-(l-1)*denom)/(quot*t0*t0*denom**3)
1221 ft(6)=(320.0d0+80.0d0*dtanh((betaT-320.0d0)/80.0d0))/
1223 ftprim(6)=1.0d0/(320.0d0*dcosh((betaT-320.0d0)/80.0d0)**2)
1224 ftbis(6)=-2.0d0*dtanh((betaT-320.0d0)/80.0d0)
1225 & /(320.0d0*80.0d0*dcosh((betaT-320.0d0)/80.0d0)**3)
1226 #elif defined(FUNCT)
1235 else if (rescale_mode.eq.2) then
1243 logfac=1.0d0/dlog(eplus+eminus)
1244 tanhT=(eplus-eminus)/(eplus+eminus)
1245 fT(l)=1.12692801104297249644d0*logfac
1246 ftprim(l)=-l*quotl1*ft(l)*tanhT*logfac/T0
1247 ftbis(l)=(l-1)*ftprim(l)/(quot*T0)-
1248 & 2*l*quotl1/T0*logfac*
1249 & (2*l*quotl1*ft(l)/(T0*(eplus+eminus)**2)
1253 ft(6)=(320.0d0+80.0d0*dtanh((betaT-320.0d0)/80.0d0))/
1255 ftprim(6)=1.0d0/(320.0d0*dcosh((betaT-320.0d0)/80.0d0)**2)
1256 ftbis(6)=-2.0d0*dtanh((betaT-320.0d0)/80.0d0)
1257 & /(320.0d0*80.0d0*dcosh((betaT-320.0d0)/80.0d0)**3)
1258 #elif defined(FUNCT)
1267 else if (rescale_mode.eq.0) then
1273 write (iout,*) "Error in WHAM_CALC: wrong RESCALE_MODE",
1280 c--------------------------------------------------------------------
1281 subroutine sum_ene(t,iparm,ft,etot)
1283 include 'DIMENSIONS'
1284 include 'DIMENSIONS.ZSCOPT'
1285 include 'DIMENSIONS.FREE'
1286 include 'COMMON.CONTROL'
1287 include 'COMMON.FFIELD'
1288 include "COMMON.SBRIDGE"
1289 include "COMMON.ENERGIES"
1290 include "COMMON.IOUNITS"
1292 double precision fT(6)
1293 double precision etot,evdw,evdw_t,evdw2,ees,evdw1,ebe,etors,
1294 & escloc,ehpb,ecorr,ecorr5,ecorr6,eello_turn4,eello_turn3,
1295 & eturn6,eel_loc,edihcnstr,etors_d,estr,evdw2_14,esccor,
1297 evdw=enetb(21,t,iparm)
1298 evdw_t=enetb(1,t,iparm)
1300 evdw2_14=enetb(17,t,iparm)
1301 evdw2=enetb(2,t,iparm)+evdw2_14
1303 evdw2=enetb(2,t,iparm)
1307 ees=enetb(3,t,iparm)
1308 evdw1=enetb(16,t,iparm)
1310 ees=enetb(3,t,iparm)
1313 ecorr=enetb(4,t,iparm)
1314 ecorr5=enetb(5,t,iparm)
1315 ecorr6=enetb(6,t,iparm)
1316 eel_loc=enetb(7,t,iparm)
1317 eello_turn3=enetb(8,t,iparm)
1318 eello_turn4=enetb(9,t,iparm)
1319 eturn6=enetb(10,t,iparm)
1320 ebe=enetb(11,t,iparm)
1321 escloc=enetb(12,t,iparm)
1322 etors=enetb(13,t,iparm)
1323 etors_d=enetb(14,t,iparm)
1324 ehpb=enetb(15,t,iparm)
1325 estr=enetb(18,t,iparm)
1326 esccor=enetb(19,t,iparm)
1327 edihcnstr=enetb(20,t,iparm)
1328 eliptran=enetb(22,t,iparm)
1329 esaxs=enetb(26,t,iparm)
1331 if (shield_mode.gt.0) then
1332 etot=ft(1)*wsc*(evdw+ft(6)*evdw_t)+ft(1)*wscp*evdw2
1334 & +ft(1)*wvdwpp*evdw1
1335 & +wang*ebe+ft(1)*wtor*etors+wscloc*escloc
1336 c & +wstrain*ehpb+nss*ebr+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5
1337 & +wstrain*ehpb+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5
1338 & +ft(5)*wcorr6*ecorr6+ft(3)*wturn4*eello_turn4
1339 & +ft(2)*wturn3*eello_turn3
1340 & +ft(5)*wturn6*eturn6+ft(2)*wel_loc*eel_loc
1341 & +edihcnstr+ft(2)*wtor_d*etors_d+ft(1)*wsccor*esccor
1342 & +wbond*estr+wliptran*eliptran+wsaxs*esaxs
1344 etot=wsc*(evdw+ft(6)*evdw_t)+wscp*evdw2+ft(1)*welec*ees
1346 & +wang*ebe+ft(1)*wtor*etors+wscloc*escloc
1347 c & +wstrain*ehpb+nss*ebr+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5
1348 & +wstrain*ehpb+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5
1349 & +ft(5)*wcorr6*ecorr6+ft(3)*wturn4*eello_turn4
1350 & +ft(2)*wturn3*eello_turn3
1351 & +ft(5)*wturn6*eturn6+ft(2)*wel_loc*eel_loc
1352 & +edihcnstr+ft(2)*wtor_d*etors_d+ft(1)*wsccor*esccor
1353 & +wbond*estr+wliptran*eliptran+wsaxs*esaxs
1356 if (shield_mode.gt.0) then
1357 etot=ft(1)*wsc*(evdw+ft(6)*evdw_t)+ft(1)*wscp*evdw2
1358 & +ft(1)*welec*(ees+evdw1)
1359 & +wang*ebe+ft(1)*wtor*etors+wscloc*escloc
1360 c & +wstrain*ehpb+nss*ebr+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5
1361 & +wstrain*ehpb+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5
1362 & +ft(5)*wcorr6*ecorr6+ft(3)*wturn4*eello_turn4
1363 & +ft(2)*wturn3*eello_turn3
1364 & +ft(5)*wturn6*eturn6+ft(2)*wel_loc*eel_loc+edihcnstr
1365 & +ft(2)*wtor_d*etors_d+ft(1)*wsccor*esccor
1366 & +wbond*estr+wliptran*eliptran+wsaxs*esaxs
1368 etot=wsc*(evdw+ft(6)*evdw_t)+wscp*evdw2
1369 & +ft(1)*welec*(ees+evdw1)
1370 & +wang*ebe+ft(1)*wtor*etors+wscloc*escloc
1371 c & +wstrain*ehpb+nss*ebr+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5
1372 & +wstrain*ehpb+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5
1373 & +ft(5)*wcorr6*ecorr6+ft(3)*wturn4*eello_turn4
1374 & +ft(2)*wturn3*eello_turn3
1375 & +ft(5)*wturn6*eturn6+ft(2)*wel_loc*eel_loc+edihcnstr
1376 & +ft(2)*wtor_d*etors_d+ft(1)*wsccor*esccor
1377 & +wbond*estr+wliptran*eliptran+wsaxs*esaxs
1382 c--------------------------------------------------------------------
1383 subroutine sum_ene_deriv(t,iparm,ft,ftprim,ftbis,etot,eprim,ebis)
1385 include 'DIMENSIONS'
1386 include 'DIMENSIONS.ZSCOPT'
1387 include 'DIMENSIONS.FREE'
1388 include 'COMMON.CONTROL'
1389 include 'COMMON.FFIELD'
1390 include "COMMON.SBRIDGE"
1391 include 'COMMON.ENERGIES'
1392 include "COMMON.HOMOLOGY"
1393 include "COMMON.IOUNITS"
1395 double precision fT(6),fTprim(6),fTbis(6),
1397 double precision etot,evdw,evdw_t,evdw2,ees,evdw1,ebe,etors,
1398 & escloc,ehpb,ecorr,ecorr5,ecorr6,eello_turn4,eello_turn3,
1399 & eturn6,eel_loc,edihcnstr,etors_d,estr,evdw2_14,esccor,
1400 & eliptran,esaxs,ehomology_constr,edfadis,edfator,edfanei,edfabet
1401 evdw=enetb(21,t,iparm)
1402 evdw_t=enetb(1,t,iparm)
1404 evdw2_14=enetb(17,t,iparm)
1405 evdw2=enetb(2,t,iparm)+evdw2_14
1407 evdw2=enetb(2,t,iparm)
1411 ees=enetb(3,t,iparm)
1412 evdw1=enetb(16,t,iparm)
1414 ees=enetb(3,t,iparm)
1417 ecorr=enetb(4,t,iparm)
1418 ecorr5=enetb(5,t,iparm)
1419 ecorr6=enetb(6,t,iparm)
1420 eel_loc=enetb(7,t,iparm)
1421 eello_turn3=enetb(8,t,iparm)
1422 eello_turn4=enetb(9,t,iparm)
1423 eturn6=enetb(10,t,iparm)
1424 ebe=enetb(11,t,iparm)
1425 escloc=enetb(12,t,iparm)
1426 etors=enetb(13,t,iparm)
1427 etors_d=enetb(14,t,iparm)
1428 ehpb=enetb(15,t,iparm)
1429 estr=enetb(18,t,iparm)
1430 esccor=enetb(19,t,iparm)
1431 edihcnstr=enetb(20,t,iparm)
1432 eliptran=enetb(22,t,iparm)
1433 esaxs=enetb(26,t,iparm)
1434 ehomology_constr=enetb(27,t,iparm)
1435 edfadis=enetb(28,t,iparm)
1436 edfator=enetb(29,t,iparm)
1437 edfanei=enetb(30,t,iparm)
1438 edfabet=enetb(31,t,iparm)
1440 if (shield_mode.gt.0) then
1441 etot=ft(1)*wsc*(evdw+ft(6)*evdw_t)+ft(1)*wscp*evdw2
1443 & +ft(1)*wvdwpp*evdw1
1444 & +wang*ebe+ft(1)*wtor*etors+wscloc*escloc
1445 c & +wstrain*ehpb+nss*ebr+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5
1446 & +wstrain*ehpb+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5
1447 & +ft(5)*wcorr6*ecorr6+ft(3)*wturn4*eello_turn4
1448 & +ft(2)*wturn3*eello_turn3
1449 & +ft(5)*wturn6*eturn6+ft(2)*wel_loc*eel_loc
1450 & +edihcnstr+ft(2)*wtor_d*etors_d+ft(1)*wsccor*esccor
1451 & +wbond*estr+wliptran*eliptran+wsaxs*esaxs
1452 eprim=ftprim(1)*(ft(6)*evdw_t+evdw)
1453 C & +ftprim(6)*evdw_t
1454 & +ftprim(1)*wscp*evdw2
1455 & +ftprim(1)*welec*ees
1456 & +ftprim(1)*wvdwpp*evdw1
1457 & +ftprim(1)*wtor*etors+
1458 & ftprim(3)*wcorr*ecorr+ftprim(4)*wcorr5*ecorr5+
1459 & ftprim(5)*wcorr6*ecorr6+ftprim(3)*wturn4*eello_turn4+
1460 & ftprim(2)*wturn3*eello_turn3+ftprim(5)*wturn6*eturn6+
1461 & ftprim(2)*wel_loc*eel_loc+ftprim(2)*wtor_d*etors_d+
1462 & ftprim(1)*wsccor*esccor
1463 ebis=ftbis(1)*wsc*(evdw+ft(6)*evdw_t)
1464 & +ftbis(1)*wscp*evdw2+
1465 & ftbis(1)*welec*ees
1466 & +ftbis(1)*wvdwpp*evdw
1467 & +ftbis(1)*wtor*etors+
1468 & ftbis(3)*wcorr*ecorr+ftbis(4)*wcorr5*ecorr5+
1469 & ftbis(5)*wcorr6*ecorr6+ftbis(3)*wturn4*eello_turn4+
1470 & ftbis(2)*wturn3*eello_turn3+ftbis(5)*wturn6*eturn6+
1471 & ftbis(2)*wel_loc*eel_loc+ftbis(2)*wtor_d*etors_d+
1472 & ftbis(1)*wsccor*esccor
1474 etot=wsc*(evdw+ft(6)*evdw_t)+wscp*evdw2+ft(1)*welec*ees
1476 & +wang*ebe+ft(1)*wtor*etors+wscloc*escloc
1477 c & +wstrain*ehpb+nss*ebr+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5
1478 & +wstrain*ehpb+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5
1479 & +ft(5)*wcorr6*ecorr6+ft(3)*wturn4*eello_turn4
1480 & +ft(2)*wturn3*eello_turn3
1481 & +ft(5)*wturn6*eturn6+ft(2)*wel_loc*eel_loc
1482 & +edihcnstr+ft(2)*wtor_d*etors_d+ft(1)*wsccor*esccor
1483 & +wbond*estr+wliptran*eliptran+wsaxs*esaxs
1484 eprim=ftprim(6)*evdw_t+ftprim(1)*welec*ees
1485 & +ftprim(1)*wtor*etors+
1486 & ftprim(3)*wcorr*ecorr+ftprim(4)*wcorr5*ecorr5+
1487 & ftprim(5)*wcorr6*ecorr6+ftprim(3)*wturn4*eello_turn4+
1488 & ftprim(2)*wturn3*eello_turn3+ftprim(5)*wturn6*eturn6+
1489 & ftprim(2)*wel_loc*eel_loc+ftprim(2)*wtor_d*etors_d+
1490 & ftprim(1)*wsccor*esccor
1491 ebis=ftbis(1)*welec*ees+ftbis(1)*wtor*etors+
1492 & ftbis(3)*wcorr*ecorr+ftbis(4)*wcorr5*ecorr5+
1493 & ftbis(5)*wcorr6*ecorr6+ftbis(3)*wturn4*eello_turn4+
1494 & ftbis(2)*wturn3*eello_turn3+ftbis(5)*wturn6*eturn6+
1495 & ftbis(2)*wel_loc*eel_loc+ftbis(2)*wtor_d*etors_d+
1496 & ftbis(1)*wsccor*esccor
1499 if (shield_mode.gt.0) then
1500 etot=ft(1)*wsc*(evdw+ft(6)*evdw_t)+ft(1)*wscp*evdw2
1501 & +ft(1)*welec*(ees+evdw1)
1502 & +wang*ebe+ft(1)*wtor*etors+wscloc*escloc
1503 c & +wstrain*ehpb+nss*ebr+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5
1504 & +wstrain*ehpb+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5
1505 & +ft(5)*wcorr6*ecorr6+ft(3)*wturn4*eello_turn4
1506 & +ft(2)*wturn3*eello_turn3
1507 & +ft(5)*wturn6*eturn6+ft(2)*wel_loc*eel_loc+edihcnstr
1508 & +ft(2)*wtor_d*etors_d+ft(1)*wsccor*esccor
1509 & +wbond*estr+wliptran*eliptran+wsaxs*esaxs
1510 eprim=ftprim(1)*(evdw+ft(6)*evdw_t)
1511 & +ftprim(1)*welec*(ees+evdw1)
1512 & +ftprim(1)*wtor*etors+
1513 & ftprim(1)*wscp*evdw2+
1514 & ftprim(3)*wcorr*ecorr+ftprim(4)*wcorr5*ecorr5+
1515 & ftprim(5)*wcorr6*ecorr6+ftprim(3)*wturn4*eello_turn4+
1516 & ftprim(2)*wturn3*eello_turn3+ftprim(5)*wturn6*eturn6+
1517 & ftprim(2)*wel_loc*eel_loc+ftprim(2)*wtor_d*etors_d+
1518 & ftprim(1)*wsccor*esccor
1519 ebis= ftbis(1)*(evdw+ft(6)*evdw_t)
1520 & +ftbis(1)*wscp*evdw2
1521 & +ftbis(1)*welec*(ees+evdw1)+ftbis(1)*wtor*etors+
1522 & ftbis(3)*wcorr*ecorr+ftbis(4)*wcorr5*ecorr5+
1523 & ftbis(5)*wcorr6*ecorr6+ftbis(3)*wturn4*eello_turn4+
1524 & ftbis(2)*wturn3*eello_turn3+ftbis(5)*wturn6*eturn6+
1525 & ftbis(2)*wel_loc*eel_loc+ftbis(2)*wtor_d*etors_d+
1526 & ftprim(1)*wsccor*esccor
1528 etot=wsc*(evdw+ft(6)*evdw_t)+wscp*evdw2
1529 & +ft(1)*welec*(ees+evdw1)
1530 & +wang*ebe+ft(1)*wtor*etors+wscloc*escloc
1531 c & +wstrain*ehpb+nss*ebr+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5
1532 & +wstrain*ehpb+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5
1533 & +ft(5)*wcorr6*ecorr6+ft(3)*wturn4*eello_turn4
1534 & +ft(2)*wturn3*eello_turn3
1535 & +ft(5)*wturn6*eturn6+ft(2)*wel_loc*eel_loc+edihcnstr
1536 & +ft(2)*wtor_d*etors_d+ft(1)*wsccor*esccor
1537 & +wbond*estr+wliptran*eliptran+wsaxs*esaxs
1538 eprim=ftprim(6)*evdw_t+ftprim(1)*welec*(ees+evdw1)
1539 & +ftprim(1)*wtor*etors+
1540 & ftprim(3)*wcorr*ecorr+ftprim(4)*wcorr5*ecorr5+
1541 & ftprim(5)*wcorr6*ecorr6+ftprim(3)*wturn4*eello_turn4+
1542 & ftprim(2)*wturn3*eello_turn3+ftprim(5)*wturn6*eturn6+
1543 & ftprim(2)*wel_loc*eel_loc+ftprim(2)*wtor_d*etors_d+
1544 & ftprim(1)*wsccor*esccor
1545 ebis=ftbis(1)*welec*(ees+evdw1)+ftbis(1)*wtor*etors+
1546 & ftbis(3)*wcorr*ecorr+ftbis(4)*wcorr5*ecorr5+
1547 & ftbis(5)*wcorr6*ecorr6+ftbis(3)*wturn4*eello_turn4+
1548 & ftbis(2)*wturn3*eello_turn3+ftbis(5)*wturn6*eturn6+
1549 & ftbis(2)*wel_loc*eel_loc+ftbis(2)*wtor_d*etors_d+
1550 & ftprim(1)*wsccor*esccor