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.FFIELD"
35 include "COMMON.SBRIDGE"
37 include "COMMON.ENEPS"
38 include "COMMON.SHIELD"
39 integer MaxPoint,MaxPointProc
40 parameter (MaxPoint=MaxStr,
41 & MaxPointProc=MaxStr_Proc)
42 double precision finorm_max,potfac,entmin,entmax,expfac,vf
43 parameter (finorm_max=1.0d0)
45 integer i,ii,j,jj,k,kk,l,m,ind,iter,t,tmax,ient,ientmax,iln
46 integer start,end,iharm,ib,iib,nbin1,nbin,nbin_rms,nbin_rgy,
47 & nbin_rmsrgy,liczba,iparm,nFi,indrgy,indrms
48 integer htot(0:MaxHdim),histent(0:2000)
49 double precision v(MaxPointProc,MaxR,MaxT_h,Max_Parm)
50 double precision energia(0:max_ene)
52 integer tmax_t,upindE_p
53 double precision fi_p(MaxR,MaxT_h,Max_Parm),
54 & fimax_p(MaxR,MaxT_h,Max_Parm)
55 double precision sumW_p(0:nGridT,Max_Parm),
56 & sumE_p(0:nGridT,Max_Parm),sumEsq_p(0:nGridT,Max_Parm),
57 & sumQ_p(MaxQ1,0:nGridT,Max_Parm),
58 & sumQsq_p(MaxQ1,0:nGridT,Max_Parm),
59 & sumEQ_p(MaxQ1,0:nGridT,Max_Parm),
60 & sumEprim_p(MaxQ1,0:nGridT,Max_Parm),
61 & sumEbis_p(0:nGridT,Max_Parm)
62 double precision hfin_p(0:MaxHdim,maxT_h),
63 & hfin_ent_p(0:MaxHdim),histE_p(0:maxindE),sumH,
64 & hrmsrgy_p(0:MaxBinRgy,0:MaxBinRms,maxT_h)
65 double precision rgymin_t,rmsmin_t,rgymax_t,rmsmax_t
66 double precision potEmin_t,entmin_p,entmax_p
67 double precision ePMF,ePMF_q
68 double precision weimax_(0:ngridT)
69 integer histent_p(0:2000)
70 logical lprint /.true./
72 double precision delta_T /1.0d0/
73 double precision rgymin,rmsmin,rgymax,rmsmax
74 double precision sumW(0:NGridT,Max_Parm),sumE(0:NGridT,Max_Parm),
75 & sumEsq(0:NGridT,Max_Parm),sumQ(MaxQ1,0:NGridT,Max_Parm),
76 & sumQsq(MaxQ1,0:NGridT,Max_Parm),sumEQ(MaxQ1,0:NGridT,Max_Parm),
77 & sumEprim(0:NGridT,Max_Parm),sumEbis(0:NGridT,Max_Parm),betaT,
79 double precision fi(MaxR,maxT_h,Max_Parm),
80 & fimax(MaxR,maxT_h,Max_Parm),
81 & dd,dd1,dd2,hh,dmin,denom,finorm,avefi,pom,
82 & hfin(0:MaxHdim,maxT_h),histE(0:maxindE),
83 & hrmsrgy(0:MaxBinRgy,0:MaxBinRms,maxT_h),
85 & hfin_ent(0:MaxHdim),vmax,aux,weimax(0:nGridT,Max_Parm)
86 double precision fT(6),fTprim(6),fTbis(6),quot,quotl1,quotl,kfacl,
87 & eprim,ebis,temper,kfac/2.4d0/,T0/300.0d0/,startGridT/200.0d0/,
88 & eplus,eminus,logfac,tanhT,tt
89 double precision etot,evdw,evdw_t,evdw2,ees,evdw1,ebe,etors,
90 & escloc,ehpb,ecorr,ecorr5,ecorr6,eello_turn4,eello_turn3,
91 & eturn6,eel_loc,edihcnstr,etors_d,estr,evdw2_14,esccor,
94 integer ind_point(maxpoint),upindE,indE
103 write(licz2,'(bz,i2.2)') islice
105 write (iout,'(//80(1h-)/"Solving WHAM equations for slice",
106 & i2/80(1h-)//)') islice
107 write (iout,*) "delta",delta," nbin1",nbin1
108 write (iout,*) "MaxN",MaxN," MaxQ",MaxQ," MaHdim",MaxHdim
126 c write (iout,*) "i",i," potE",(potE(i,j),j=1,nParmset)
128 if (potE(i,j).le.potEmin) potEmin=potE(i,j)
130 if (q(nQ+1,i).lt.rmsmin) rmsmin=q(nQ+1,i)
131 if (q(nQ+1,i).gt.rmsmax) rmsmax=q(nQ+1,i)
132 if (q(nQ+2,i).lt.rgymin) rgymin=q(nQ+2,i)
133 if (q(nQ+2,i).gt.rgymax) rgymax=q(nQ+2,i)
136 ind=(q(j,i)-dmin+1.0d-8)/delta
138 ind_point(i)=ind_point(i)+ind
140 ind_point(i)=ind_point(i)+nbin1**(j-1)*ind
142 c write (iout,*) "i",i," j",j," q",q(j,i)," ind_point",
145 if (ind_point(i).lt.0 .or. ind_point(i).gt.MaxHdim) then
146 write (iout,*) "Error - index exceeds range for point",i,
147 & " q=",q(j,i)," ind",ind_point(i)
149 write (iout,*) "Processor",me1
151 call MPI_Abort(MPI_COMM_WORLD, Ierror, Errcode )
156 if (ind_point(i).gt.tmax) tmax=ind_point(i)
157 htot(ind_point(i))=htot(ind_point(i))+1
159 write (iout,*) "i",i,"q",(q(j,i),j=1,nQ)," ind",ind_point(i),
160 & " htot",htot(ind_point(i))
165 write (iout,*) "potEmin before reduce",potEmin
167 write (iout,'(a)') "Numbers of counts in Q bins"
169 if (htot(t).gt.0) then
170 write (iout,'(i15,$)') t
173 jj = mod(liczba,nbin1)
175 write (iout,'(i5,$)') jj
177 write (iout,'(i8)') htot(t)
181 write (iout,'(a,i3)') "Number of data points for parameter set",
183 write (iout,'(i7,$)') ((snk(m,ib,iparm,islice),m=1,nr(ib,iparm)),
185 write (iout,'(i8)') stot(islice)
191 call MPI_AllReduce(tmax,tmax_t,1,MPI_INTEGER,MPI_MAX,
194 call MPI_AllReduce(potEmin,potEmin_t,1,MPI_DOUBLE_PRECISION,
195 & MPI_MIN,WHAM_COMM,IERROR)
196 call MPI_AllReduce(rmsmin,rmsmin_t,1,MPI_DOUBLE_PRECISION,
197 & MPI_MIN,WHAM_COMM,IERROR)
198 call MPI_AllReduce(rmsmax,rmsmax_t,1,MPI_DOUBLE_PRECISION,
199 & MPI_MAX,WHAM_COMM,IERROR)
200 call MPI_AllReduce(rgymin,rgymin_t,1,MPI_DOUBLE_PRECISION,
201 & MPI_MIN,WHAM_COMM,IERROR)
202 call MPI_AllReduce(rgymax,rgymax_t,1,MPI_DOUBLE_PRECISION,
203 & MPI_MAX,WHAM_COMM,IERROR)
204 c potEmin=potEmin_t/2
210 write (iout,*) "potEmin",potEmin
212 rmsmin=deltrms*dint(rmsmin/deltrms)
213 rmsmax=deltrms*dint(rmsmax/deltrms)
214 rgymin=deltrms*dint(rgymin/deltrgy)
215 rgymax=deltrms*dint(rgymax/deltrgy)
216 nbin_rms=(rmsmax-rmsmin)/deltrms
217 nbin_rgy=(rgymax-rgymin)/deltrgy
218 write (iout,*) "rmsmin",rmsmin," rmsmax",rmsmax," rgymin",rgymin,
219 & " rgymax",rgymax," nbin_rms",nbin_rms," nbin_rgy",nbin_rgy
226 write (iout,*) "nFi",nFi
227 ! Compute the Boltzmann factor corresponing to restrain potentials in different
234 c write (9,'(3i5,f10.5)') i,(iparm,potE(i,iparm),iparm=1,nParmSet)
237 write (iout,'(2i5,21f8.2)') i,iparm,
238 & (enetb(k,i,iparm),k=1,22)
240 call restore_parm(iparm)
242 write (iout,*) wsc,wscp,welec,wvdwpp,wang,wtor,wscloc,
243 & wcorr,wcorr5,wcorr6,wturn4,wturn3,wturn6,wel_loc,
244 & wtor_d,wsccor,wbond
247 if (rescale_mode.eq.1) then
248 quot=1.0d0/(beta_h(ib,iparm)*1.987D-3*T0)
255 fT(l)=kfacl/(kfacl-1.0d0+quotl)
258 tt = 1.0d0/(beta_h(ib,iparm)*1.987D-3)
259 ft(6)=(320.0d0+80.0d0*dtanh((tt-320.0d0)/80.0d0))/320.0d0
261 ft(6)=1.0d0/(beta_h(ib,iparm)*1.987D-3*T0)
265 else if (rescale_mode.eq.2) then
266 quot=1.0d0/(T0*beta_h(ib,iparm)*1.987D-3)
270 fT(l)=1.12692801104297249644d0/
271 & dlog(dexp(quotl)+dexp(-quotl))
274 tt = 1.0d0/(beta_h(ib,iparm)*1.987D-3)
275 ft(6)=(320.0d0+80.0d0*dtanh((tt-320.0d0)/80.0d0))/320.0d0
277 ft(6)=1.0d0/(beta_h(ib,iparm)*1.987D-3*T0)
281 c write (iout,*) 1.0d0/(beta_h(ib,iparm)*1.987D-3),ft
282 else if (rescale_mode.eq.0) then
287 write (iout,*) "Error in WHAM_CALC: wrong RESCALE_MODE",
292 evdw=enetb(1,i,iparm)
293 evdw_t=enetb(21,i,iparm)
295 evdw2_14=enetb(17,i,iparm)
296 evdw2=enetb(2,i,iparm)+evdw2_14
298 evdw2=enetb(2,i,iparm)
303 evdw1=enetb(16,i,iparm)
308 ecorr=enetb(4,i,iparm)
309 ecorr5=enetb(5,i,iparm)
310 ecorr6=enetb(6,i,iparm)
311 eel_loc=enetb(7,i,iparm)
312 eello_turn3=enetb(8,i,iparm)
313 eello_turn4=enetb(9,i,iparm)
314 eturn6=enetb(10,i,iparm)
315 ebe=enetb(11,i,iparm)
316 escloc=enetb(12,i,iparm)
317 etors=enetb(13,i,iparm)
318 etors_d=enetb(14,i,iparm)
319 ehpb=enetb(15,i,iparm)
320 estr=enetb(18,i,iparm)
321 esccor=enetb(19,i,iparm)
322 edihcnstr=enetb(20,i,iparm)
323 eliptran=enetb(22,i,iparm)
324 esaxs=enetb(26,i,iparm)
327 write (iout,'(3i5,6f5.2,14f12.3)') i,ib,iparm,(ft(l),l=1,6),
328 & evdw+evdw_t,evdw2,ees,evdw1,ecorr,eel_loc,estr,ebe,escloc,
329 & etors,etors_d,eello_turn3,eello_turn4,esccor,esaxs
333 if (shield_mode.gt.0) then
334 etot=ft(1)*wsc*(evdw+ft(6)*evdw_t)+ft(1)*wscp*evdw2
336 & +ft(1)*wvdwpp*evdw1
337 & +wang*ebe+ft(1)*wtor*etors+wscloc*escloc
338 ! & +wstrain*ehpb+nss*ebr+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5
339 & +wstrain*ehpb+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5
340 & +ft(5)*wcorr6*ecorr6+ft(3)*wturn4*eello_turn4
341 & +ft(2)*wturn3*eello_turn3
342 & +ft(5)*wturn6*eturn6+ft(2)*wel_loc*eel_loc
343 & +edihcnstr+ft(2)*wtor_d*etors_d+ft(1)*wsccor*esccor
344 & +wbond*estr+wliptran*eliptran+wsaxs*esaxs
346 etot=wsc*(evdw+ft(6)*evdw_t)+wscp*evdw2+ft(1)*welec*ees
348 & +wang*ebe+ft(1)*wtor*etors+wscloc*escloc
349 c & +wstrain*ehpb+nss*ebr+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5
350 & +wstrain*ehpb+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5
351 & +ft(5)*wcorr6*ecorr6+ft(3)*wturn4*eello_turn4
352 & +ft(2)*wturn3*eello_turn3
353 & +ft(5)*wturn6*eturn6+ft(2)*wel_loc*eel_loc
354 & +edihcnstr+ft(2)*wtor_d*etors_d+ft(1)*wsccor*esccor
355 & +wbond*estr+wliptran*eliptran+wsaxs*esaxs
358 if (shield_mode.gt.0) then
359 etot=ft(1)*wsc*(evdw+ft(6)*evdw_t)+ft(1)*wscp*evdw2
360 & +ft(1)*welec*(ees+evdw1)
361 & +wang*ebe+ft(1)*wtor*etors+wscloc*escloc
362 c & +wstrain*ehpb+nss*ebr+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5
363 & +wstrain*ehpb+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5
364 & +ft(5)*wcorr6*ecorr6+ft(3)*wturn4*eello_turn4
365 & +ft(2)*wturn3*eello_turn3
366 & +ft(5)*wturn6*eturn6+ft(2)*wel_loc*eel_loc+edihcnstr
367 & +ft(2)*wtor_d*etors_d+ft(1)*wsccor*esccor
368 & +wbond*estr+wliptran*eliptran+wsaxs*esaxs
370 etot=wsc*(evdw+ft(6)*evdw_t)+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
384 write (iout,*) i,iparm,1.0d0/(beta_h(ib,iparm)*1.987D-3),
388 if (iparm.eq.1 .and. ib.eq.1) then
389 write (iout,*)"Conformation",i
392 energia(k)=enetb(k,i,iparm)
394 call enerprint(energia(0),fT)
401 Econstr=Econstr+Kh(j,kk,ib,iparm)
402 & *(dd-q0(j,kk,ib,iparm))**2
404 c Adaptive potential contribution
406 call PMF_energy(q(1,i),ib,kk,iparm,ePMF,ePMF_q)
410 & -beta_h(ib,iparm)*(etot-potEmin+Econstr)
412 write (iout,'(4i5,4e15.5)') i,kk,ib,iparm,
413 & etot,potEmin,etot-potEmin,v(i,kk,ib,iparm)
419 ! Simple iteration to calculate free energies corresponding to all simulation
423 ! Compute new free-energy values corresponding to the righ-hand side of the
424 ! equation and their derivatives.
425 write (iout,*) "------------------------fi"
435 vf=v(t,l,k,i)+f(l,k,i)
436 if (vf.gt.vmax) vmax=vf
444 aux=f(l,k,i)+v(t,l,k,i)-vmax
446 & denom=denom+snk(l,k,i,islice)*dexp(aux)
450 entfac(t)=-dlog(denom)-vmax
452 write (iout,*) t,"vmax",vmax," denom",denom,"entfac",entfac(t)
458 do ii=1,nR(iib,iparm)
460 fimax_p(ii,iib,iparm)=v(1,ii,iib,iparm)+entfac(1)
462 if(v(t,ii,iib,iparm)+entfac(t).gt.fimax_p(ii,iib,iparm))
463 & fimax_p(ii,iib,iparm)=v(t,ii,iib,iparm)+entfac(t)
466 fimax(ii,iib,iparm)=v(1,ii,iib,iparm)+entfac(1)
468 if(v(t,ii,iib,iparm)+entfac(t).gt.fimax(ii,iib,iparm))
469 & fimax(ii,iib,iparm)=v(t,ii,iib,iparm)+entfac(t)
476 call MPI_AllReduce(fimax_p(1,1,1),fimax(1,1,1),
477 & maxR*MaxT_h*nParmSet,MPI_DOUBLE_PRECISION,
478 & MPI_MAX,WHAM_COMM,IERROR)
482 do ii=1,nR(iib,iparm)
484 fi_p(ii,iib,iparm)=0.0d0
486 fi_p(ii,iib,iparm)=fi_p(ii,iib,iparm)
487 & +dexp(v(t,ii,iib,iparm)+entfac(t)-fimax(ii,iib,iparm))
489 write (iout,'(4i5,3e15.5)') t,ii,iib,iparm,
490 & v(t,ii,iib,iparm),entfac(t),fi_p(ii,iib,iparm)
494 fi(ii,iib,iparm)=0.0d0
496 fi(ii,iib,iparm)=fi(ii,iib,iparm)
497 & +dexp(v(t,ii,iib,iparm)+entfac(t)-fimax(ii,iib,iparm))
506 write (iout,*) "fi before MPI_Reduce me",me,' master',master
508 do ib=1,nT_h(nparmset)
509 write (iout,*) "iparm",iparm," ib",ib
510 write (iout,*) "beta=",beta_h(ib,iparm)
511 write (iout,'(8e15.5)') (fi_p(i,ib,iparm),i=1,nR(ib,iparm))
515 c write (iout,*) "REDUCE size",maxR,MaxT_h,nParmSet,
516 c & maxR*MaxT_h*nParmSet
517 c write (iout,*) "MPI_COMM_WORLD",MPI_COMM_WORLD,
518 c & " WHAM_COMM",WHAM_COMM
519 call MPI_Reduce(fi_p(1,1,1),fi(1,1,1),maxR*MaxT_h*nParmSet,
520 & MPI_DOUBLE_PRECISION,
521 & MPI_SUM,Master,WHAM_COMM,IERROR)
523 write (iout,*) "fi after MPI_Reduce nparmset",nparmset
525 write (iout,*) "iparm",iparm
527 write (iout,*) "beta=",beta_h(ib,iparm)
528 write (iout,'(8e15.5)') (fi(i,ib,iparm),i=1,nR(ib,iparm))
532 if (me1.eq.Master) then
538 fi(i,ib,iparm)=-dlog(fi(i,ib,iparm))-fimax(i,ib,iparm)
539 avefi=avefi+fi(i,ib,iparm)
545 write (iout,*) "Parameter set",iparm
547 write (iout,*) "beta=",beta_h(ib,iparm)
549 fi(i,ib,iparm)=fi(i,ib,iparm)-avefi
551 write (iout,'(8f10.5)') (fi(i,ib,iparm),i=1,nR(ib,iparm))
552 write (iout,'(8f10.5)') (f(i,ib,iparm),i=1,nR(ib,iparm))
556 ! Compute the norm of free-energy increments.
561 finorm=finorm+dabs(fi(i,ib,iparm)-f(i,ib,iparm))
562 f(i,ib,iparm)=fi(i,ib,iparm)
567 write (iout,*) 'Iteration',iter,' finorm',finorm
571 call MPI_Bcast(f(1,1,1),MaxR*MaxT_h*nParmSet,
572 & MPI_DOUBLE_PRECISION,Master,
574 call MPI_Bcast(finorm,1,MPI_DOUBLE_PRECISION,Master,
577 ! Exit, if the increment norm is smaller than pre-assigned tolerance.
578 if (finorm.lt.fimin) then
579 write (iout,*) 'Iteration converged'
586 ! Now, put together the histograms from all simulations, in order to get the
587 ! unbiased total histogram.
597 write (iout,*) "--------------hist"
601 sumW_p(i,iparm)=0.0d0
602 sumE_p(i,iparm)=0.0d0
603 sumEbis_p(i,iparm)=0.0d0
604 sumEsq_p(i,iparm)=0.0d0
606 sumQ_p(j,i,iparm)=0.0d0
607 sumQsq_p(j,i,iparm)=0.0d0
608 sumEQ_p(j,i,iparm)=0.0d0
618 sumEbis(i,iparm)=0.0d0
619 sumEsq(i,iparm)=0.0d0
621 sumQ(j,i,iparm)=0.0d0
622 sumQsq(j,i,iparm)=0.0d0
623 sumEQ(j,i,iparm)=0.0d0
629 c 8/26/05 entropy distribution
634 c ent=-dlog(entfac(t))
636 if (ent.lt.entmin_p) entmin_p=ent
637 if (ent.gt.entmax_p) entmax_p=ent
639 write (iout,*) "entmin",entmin_p," entmax",entmax_p
641 call MPI_Allreduce(entmin_p,entmin,1,MPI_DOUBLE_PRECISION,MPI_MIN,
643 call MPI_Allreduce(entmax_p,entmax,1,MPI_DOUBLE_PRECISION,MPI_MAX,
645 ientmax=entmax-entmin
646 if (ientmax.gt.2000) ientmax=2000
647 write (iout,*) "entmin",entmin," entmax",entmax," ientmax",ientmax
650 c ient=-dlog(entfac(t))-entmin
651 ient=entfac(t)-entmin
652 if (ient.le.2000) histent_p(ient)=histent_p(ient)+1
654 call MPI_Allreduce(histent_p(0),histent(0),ientmax+1,MPI_INTEGER,
655 & MPI_SUM,WHAM_COMM,IERROR)
656 if (me1.eq.Master) then
657 write (iout,*) "Entropy histogram"
659 write(iout,'(f15.4,i10)') entmin+i,histent(i)
667 if (ent.lt.entmin) entmin=ent
668 if (ent.gt.entmax) entmax=ent
670 ientmax=-dlog(entmax)-entmin
671 if (ientmax.gt.2000) ientmax=2000
673 ient=entfac(t)-entmin
674 if (ient.le.2000) histent(ient)=histent(ient)+1
676 write (iout,*) "Entropy histogram"
678 write(iout,'(2f15.4)') entmin+i,histent(i)
683 call restore_parm(iparm)
709 hrmsrgy(j,i,ib)=0.0d0
711 hrmsrgy_p(j,i,ib)=0.0d0
721 indE = aint(potE(t,iparm)-aint(potEmin))
722 if (indE.ge.0 .and. indE.le.maxinde) then
723 if (indE.gt.upindE_p) upindE_p=indE
724 histE_p(indE)=histE_p(indE)+dexp(-entfac(t))
728 expfac=dexp(-beta_h(ib,iparm)*(etot-potEmin)+entfac(t))
729 hfin_p(ind,ib)=hfin_p(ind,ib)+
730 & dexp(-beta_h(ib,iparm)*(etot-potEmin)+entfac(t))
732 indrgy=dint((q(nQ+2,t)-rgymin)/deltrgy)
733 indrms=dint((q(nQ+1,t)-rmsmin)/deltrms)
734 hrmsrgy_p(indrgy,indrms,ib)=
735 & hrmsrgy_p(indrgy,indrms,ib)+expfac
740 expfac=dexp(-beta_h(ib,iparm)*(etot-potEmin)+entfac(t))
741 hfin(ind,ib)=hfin(ind,ib)+
742 & dexp(-beta_h(ib,iparm)*(etot-potEmin)+entfac(t))
744 indrgy=dint((q(nQ+2,t)-rgymin)/deltrgy)
745 indrms=dint((q(nQ+1,t)-rmsmin)/deltrms)
746 hrmsrgy(indrgy,indrms,ib)=
747 & hrmsrgy(indrgy,indrms,ib)+expfac
753 c Thermo and ensemble averages
756 betaT=startGridT+k*delta_T
757 call temp_scalfac(betaT,ft,ftprim,ftbis,*10)
758 c write (iout,*) "ftprim",ftprim
759 c write (iout,*) "ftbis",ftbis
760 betaT=1.0d0/(1.987D-3*betaT)
761 c 7/10/18 AL Determine the max Botzmann weights for each temerature
762 call sum_ene(1,iparm,ft,etot)
763 weimax(k,iparm)=-betaT*(etot-potEmin)+entfac(1)
764 c write (iout,*) "k",k," t",1," weight",weimax(k,iparm)
770 call sum_ene(t,iparm,ft,etot)
771 weight=-betaT*(etot-potEmin)+entfac(t)
772 c write (iout,*) "k",k," t",t," weight",weight
773 if (weight.gt.weimax(k,iparm)) weimax(k,iparm)=weight
778 write (iout,*) "weimax before REDUCE"
779 write (iout,*) (weimax(k,iparm),k=0,ngridt)
782 weimax_(k)=weimax(k,iparm)
784 call MPI_Allreduce(weimax_(0),weimax(0,iparm),nGridT+1,
785 & MPI_DOUBLE_PRECISION,MPI_MAX,WHAM_COMM,IERROR)
787 write (iout,*) "weimax"
788 write (iout,*) (weimax(k,iparm),k=0,ngridt)
791 temper=startGridT+k*delta_T
792 betaT=1.0d0/(1.987D-3*temper)
793 call temp_scalfac(temper,ft,ftprim,ftbis,*10)
800 hfin_ent_p(ind)=hfin_ent_p(ind)+dexp(entfac(t))
802 hfin_ent(ind)=hfin_ent(ind)+dexp(entfac(t))
804 c write (iout,'(2i5,20f8.2)') t,t,(enetb(k,t,iparm),k=1,18)
805 c call restore_parm(iparm)
806 call sum_ene_deriv(t,iparm,ft,ftprim,ftbis,etot,eprim,ebis)
807 weight=dexp(-betaT*(etot-potEmin)+entfac(t)-weimax(k,iparm))
809 write (iout,*) "iparm",iparm," t",t," betaT",betaT,
810 & " etot",etot," entfac",entfac(t)," boltz",
811 & -betaT*(etot-potEmin)+entfac(t)," weimax",weimax(k,iparm),
812 & " weight",weight," ebis",ebis
814 etot=etot-temper*eprim
816 sumW_p(k,iparm)=sumW_p(k,iparm)+weight
817 sumE_p(k,iparm)=sumE_p(k,iparm)+etot*weight
818 sumEbis_p(k,iparm)=sumEbis_p(k,iparm)+ebis*weight
819 sumEsq_p(k,iparm)=sumEsq_p(k,iparm)+etot**2*weight
821 sumQ_p(j,k,iparm)=sumQ_p(j,k,iparm)+q(j,t)*weight
822 sumQsq_p(j,k,iparm)=sumQsq_p(j,k,iparm)+q(j,t)**2*weight
823 sumEQ_p(j,k,iparm)=sumEQ_p(j,k,iparm)
824 & +etot*q(j,t)*weight
827 sumW(k,iparm)=sumW(k,iparm)+weight
828 sumE(k,iparm)=sumE(k,iparm)+etot*weight
829 sumEbis(k,iparm)=sumEbis(k,iparm)+ebis*weight
830 sumEsq(k,iparm)=sumEsq(k,iparm)+etot**2*weight
832 sumQ(j,k,iparm)=sumQ(j,k,iparm)+q(j,t)*weight
833 sumQsq(j,k,iparm)=sumQsq(j,k,iparm)+q(j,t)**2*weight
834 sumEQ(j,k,iparm)=sumEQ(j,k,iparm)
835 & +etot*q(j,t)*weight
841 if (histout) call MPI_Reduce(hfin_p(0,ib),hfin(0,ib),nbin,
842 & MPI_DOUBLE_PRECISION,MPI_SUM,Master,WHAM_COMM,IERROR)
844 call MPI_Reduce(hrmsrgy_p(0,0,ib),hrmsrgy(0,0,ib),
845 & (MaxBinRgy+1)*(nbin_rms+1),MPI_DOUBLE_PRECISION,MPI_SUM,Master,
849 call MPI_Reduce(upindE_p,upindE,1,
850 & MPI_INTEGER,MPI_MAX,Master,WHAM_COMM,IERROR)
851 call MPI_Reduce(histE_p(0),histE(0),maxindE,
852 & MPI_DOUBLE_PRECISION,MPI_SUM,Master,WHAM_COMM,IERROR)
854 if (me1.eq.master) then
858 write (iout,'(6x,$)')
859 write (iout,'(f20.2,$)') (1.0d0/(1.987D-3*beta_h(ib,iparm)),
863 write (iout,'(/a)') 'Final histograms'
865 if (nslice.eq.1) then
866 if (separate_parset) then
867 write(licz3,"(bz,i3.3)") myparm
868 histname=prefix(:ilen(prefix))//'_par'//licz3//'.hist'
870 histname=prefix(:ilen(prefix))//'.hist'
873 if (separate_parset) then
874 write(licz3,"(bz,i3.3)") myparm
875 histname=prefix(:ilen(prefix))//'_par'//licz3//
876 & '_slice_'//licz2//'.hist'
878 histname=prefix(:ilen(prefix))//'_slice_'//licz2//'.hist'
881 #if defined(AIX) || defined(PGI)
882 open (ihist,file=histname,position='append')
884 open (ihist,file=histname,access='append')
894 if (sumH.gt.0.0d0) then
896 jj = mod(liczba,nbin1)
898 write (iout,'(f6.3,$)') dmin+(jj+0.5d0)*delta
900 & write (ihist,'(f6.3,$)') dmin+(jj+0.5d0)*delta
903 write (iout,'(e20.10,$)') hfin(t,ib)
904 if (histfile) write (ihist,'(e20.10,$)') hfin(t,ib)
906 write (iout,'(i5)') iparm
907 if (histfile) write (ihist,'(i5)') iparm
914 if (nslice.eq.1) then
915 if (separate_parset) then
916 write(licz3,"(bz,i3.3)") myparm
917 histname=prefix(:ilen(prefix))//"_par"//licz3//'.ent'
919 histname=prefix(:ilen(prefix))//'.ent'
922 if (separate_parset) then
923 write(licz3,"(bz,i3.3)") myparm
924 histname=prefix(:ilen(prefix))//'par_'//licz3//
925 & '_slice_'//licz2//'.ent'
927 histname=prefix(:ilen(prefix))//'_slice_'//licz2//'.ent'
930 #if defined(AIX) || defined(PGI)
931 open (ihist,file=histname,position='append')
933 open (ihist,file=histname,access='append')
935 write (ihist,'(a)') "# Microcanonical entropy"
937 write (ihist,'(f8.0,$)') dint(potEmin)+i
938 if (histE(i).gt.0.0e0) then
939 write (ihist,'(f15.5,$)') dlog(histE(i))
941 write (ihist,'(f15.5,$)') 0.0d0
947 write (iout,*) "Microcanonical entropy"
949 write (iout,'(f8.0,$)') dint(potEmin)+i
950 if (histE(i).gt.0.0e0) then
951 write (iout,'(f15.5,$)') dlog(histE(i))
953 write (iout,'(f15.5,$)') 0.0d0
958 if (nslice.eq.1) then
959 if (separate_parset) then
960 write(licz3,"(bz,i3.3)") myparm
961 histname=prefix(:ilen(prefix))//'_par'//licz3//'.rmsrgy'
963 histname=prefix(:ilen(prefix))//'.rmsrgy'
966 if (separate_parset) then
967 write(licz3,"(bz,i3.3)") myparm
968 histname=prefix(:ilen(prefix))//'_par'//licz3//
969 & '_slice_'//licz2//'.rmsrgy'
971 histname=prefix(:ilen(prefix))//'_slice_'//licz2//'.rmsrgy'
974 #if defined(AIX) || defined(PGI)
975 open (ihist,file=histname,position='append')
977 open (ihist,file=histname,access='append')
981 write(ihist,'(2f8.2,$)')
982 & rgymin+deltrgy*j,rmsmin+deltrms*i
984 if (hrmsrgy(j,i,ib).gt.0.0d0) then
985 write(ihist,'(e14.5,$)')
986 & -dlog(hrmsrgy(j,i,ib))/beta_h(ib,iparm)
989 write(ihist,'(e14.5,$)') 1.0d6
992 write (ihist,'(i2)') iparm
1000 call MPI_Reduce(hfin_ent_p(0),hfin_ent(0),nbin,
1001 & MPI_DOUBLE_PRECISION,MPI_SUM,Master,WHAM_COMM,IERROR)
1002 call MPI_Reduce(sumW_p(0,1),sumW(0,1),(nGridT+1)*nParmSet,
1003 & MPI_DOUBLE_PRECISION,MPI_SUM,Master,WHAM_COMM,IERROR)
1004 call MPI_Reduce(sumE_p(0,1),sumE(0,1),(nGridT+1)*nParmSet,
1005 & MPI_DOUBLE_PRECISION,MPI_SUM,Master,WHAM_COMM,IERROR)
1006 call MPI_Reduce(sumEbis_p(0,1),sumEbis(0,1),(nGridT+1)*nParmSet,
1007 & MPI_DOUBLE_PRECISION,MPI_SUM,Master,WHAM_COMM,IERROR)
1008 call MPI_Reduce(sumEsq_p(0,1),sumEsq(0,1),(nGridT+1)*nParmSet,
1009 & MPI_DOUBLE_PRECISION,MPI_SUM,Master,WHAM_COMM,IERROR)
1010 call MPI_Reduce(sumQ_p(1,0,1),sumQ(1,0,1),
1011 & MaxQ1*(nGridT+1)*nParmSet,MPI_DOUBLE_PRECISION,MPI_SUM,Master,
1013 call MPI_Reduce(sumQsq_p(1,0,1),sumQsq(1,0,1),
1014 & MaxQ1*(nGridT+1)*nParmSet,MPI_DOUBLE_PRECISION,MPI_SUM,Master,
1016 call MPI_Reduce(sumEQ_p(1,0,1),sumEQ(1,0,1),
1017 & MaxQ1*(nGridT+1)*nParmSet,MPI_DOUBLE_PRECISION,MPI_SUM,Master,
1019 if (me.eq.master) then
1021 write (iout,'(/a)') 'Thermal characteristics of folding'
1022 if (nslice.eq.1) then
1025 nazwa=prefix(:ilen(prefix))//"_slice_"//licz2
1028 if (nparmset.eq.1 .and. .not.separate_parset) then
1029 nazwa=nazwa(:iln)//".thermal"
1030 else if (nparmset.eq.1 .and. separate_parset) then
1031 write(licz3,"(bz,i3.3)") myparm
1032 nazwa=nazwa(:iln)//"_par_"//licz3//".thermal"
1035 if (nparmset.gt.1) then
1036 write(licz3,"(bz,i3.3)") iparm
1037 nazwa=nazwa(:iln)//"_par_"//licz3//".thermal"
1040 if (separate_parset) then
1041 write (iout,'(a,i3)') "Parameter set",myparm
1043 write (iout,'(a,i3)') "Parameter set",iparm
1046 sumE(i,iparm)=sumE(i,iparm)/sumW(i,iparm)
1047 sumEbis(i,iparm)=(startGridT+i*delta_T)*sumEbis(i,iparm)/
1049 sumEsq(i,iparm)=(sumEsq(i,iparm)/sumW(i,iparm)
1050 & -sumE(i,iparm)**2)/(1.987D-3*(startGridT+i*delta_T)**2)
1052 sumQ(j,i,iparm)=sumQ(j,i,iparm)/sumW(i,iparm)
1053 sumQsq(j,i,iparm)=sumQsq(j,i,iparm)/sumW(i,iparm)
1054 & -sumQ(j,i,iparm)**2
1055 sumEQ(j,i,iparm)=sumEQ(j,i,iparm)/sumW(i,iparm)
1056 & -sumQ(j,i,iparm)*sumE(i,iparm)
1058 sumW(i,iparm)=(-dlog(sumW(i,iparm))-weimax(i,iparm))*(1.987D-3*
1059 & (startGridT+i*delta_T))+potEmin
1060 write (iout,'(f7.1,2f15.5,$)') startGridT+i*delta_T,
1061 & sumW(i,iparm),sumE(i,iparm)
1062 write (iout,'(f10.5,$)') (sumQ(j,i,iparm),j=1,nQ+2)
1063 write (iout,'(e15.5,$)') sumEsq(i,iparm)-sumEbis(i,iparm),
1064 & (sumQsq(j,i,iparm),j=1,nQ+2),(sumEQ(j,i,iparm),j=1,nQ+2)
1066 write (34,'(f7.1,2f15.5,$)') startGridT+i*delta_T,
1067 & sumW(i,iparm),sumE(i,iparm)
1068 write (34,'(f10.5,$)') (sumQ(j,i,iparm),j=1,nQ+2)
1069 write (34,'(e15.5,$)') sumEsq(i,iparm)-sumEbis(i,iparm),
1070 & (sumQsq(j,i,iparm),j=1,nQ+2),(sumEQ(j,i,iparm),j=1,nQ+2)
1078 if (hfin_ent(t).gt.0.0d0) then
1080 jj = mod(liczba,nbin1)
1081 write (iout,'(f6.3,e20.10," ent")') dmin+(jj+0.5d0)*delta,
1083 if (histfile) write (ihist,'(f6.3,e20.10," ent")')
1084 & dmin+(jj+0.5d0)*delta,
1088 if (histfile) close(ihist)
1092 ! Write data for zscore
1093 if (nslice.eq.1) then
1094 zscname=prefix(:ilen(prefix))//".zsc"
1096 zscname=prefix(:ilen(prefix))//"_slice_"//licz2//".zsc"
1098 #if defined(AIX) || defined(PGI)
1099 open (izsc,file=prefix(:ilen(prefix))//'.zsc',position='append')
1101 open (izsc,file=prefix(:ilen(prefix))//'.zsc',access='append')
1103 write (izsc,'("NQ=",i1," NPARM=",i1)') nQ,nParmSet
1105 write (izsc,'("NT=",i1)') nT_h(iparm)
1107 write (izsc,'("TEMP=",f6.1," NR=",i2," SNK=",$)')
1108 & 1.0d0/(beta_h(ib,iparm)*1.987D-3),nR(ib,iparm)
1109 jj = min0(nR(ib,iparm),7)
1110 write (izsc,'(i8,$)') (snk(i,ib,iparm,islice),i=1,jj)
1111 write (izsc,'(a1,$)') (" ",i=22+8*jj+1,79)
1112 write (izsc,'("&")')
1113 if (nR(ib,iparm).gt.7) then
1114 do ii=8,nR(ib,iparm),9
1115 jj = min0(nR(ib,iparm),ii+8)
1116 write (izsc,'(i8,$)') (snk(i,ib,iparm,islice),i=ii,jj)
1117 write (izsc,'(a1,$') (" ",i=(jj-ii+1)*8+1,79)
1118 write (izsc,'("&")')
1121 write (izsc,'("FI=",$)')
1122 jj=min0(nR(ib,iparm),7)
1123 write (izsc,'(f10.5,$)') (fi(i,ib,iparm),i=1,jj)
1124 write (izsc,'(a1,$)') (" ",i=3+10*jj+1,79)
1125 write (izsc,'("&")')
1126 if (nR(ib,iparm).gt.7) then
1127 do ii=8,nR(ib,iparm),9
1128 jj = min0(nR(ib,iparm),ii+8)
1129 write (izsc,'(f10.5,$)') (fi(i,ib,iparm),i=ii,jj)
1130 if (jj.eq.nR(ib,iparm)) then
1133 write (izsc,'(a1,$)') (" ",i=10*(jj-ii+1)+1,79)
1134 write (izsc,'(t80,"&")')
1139 write (izsc,'("KH=",$)')
1140 write (izsc,'(f7.2,$)') (Kh(j,i,ib,iparm),j=1,nQ)
1141 write (izsc,'(" Q0=",$)')
1142 write (izsc,'(f7.5,$)') (q0(j,i,ib,iparm),j=1,nQ)
1157 c------------------------------------------------------------------------
1158 subroutine temp_scalfac(betaT,ft,ftprim,ftbis,*)
1160 include "DIMENSIONS"
1161 include "DIMENSIONS.FREE"
1162 include "COMMON.CONTROL"
1163 include "COMMON.FREE"
1164 include "COMMON.IOUNITS"
1165 double precision fT(6),fTprim(6),fTbis(6),quot,quotl1,quotl,kfacl,
1166 & kfac/2.4d0/,T0/300.0d0/,startGridT/200.0d0/,
1167 & logfac,tanhT,betaT,denom,eplus,eminus
1169 if (rescale_mode.eq.1) then
1177 denom=kfacl-1.0d0+quotl
1179 ftprim(l)=-l*ft(l)*quotl1/(T0*denom)
1180 ftbis(l)=l*kfacl*quotl1*
1181 & (2*l*quotl-(l-1)*denom)/(quot*t0*t0*denom**3)
1184 ft(6)=(320.0d0+80.0d0*dtanh((betaT-320.0d0)/80.0d0))/
1186 ftprim(6)=1.0d0/(320.0d0*dcosh((betaT-320.0d0)/80.0d0)**2)
1187 ftbis(6)=-2.0d0*dtanh((betaT-320.0d0)/80.0d0)
1188 & /(320.0d0*80.0d0*dcosh((betaT-320.0d0)/80.0d0)**3)
1189 #elif defined(FUNCT)
1198 else if (rescale_mode.eq.2) then
1206 logfac=1.0d0/dlog(eplus+eminus)
1207 tanhT=(eplus-eminus)/(eplus+eminus)
1208 fT(l)=1.12692801104297249644d0*logfac
1209 ftprim(l)=-l*quotl1*ft(l)*tanhT*logfac/T0
1210 ftbis(l)=(l-1)*ftprim(l)/(quot*T0)-
1211 & 2*l*quotl1/T0*logfac*
1212 & (2*l*quotl1*ft(l)/(T0*(eplus+eminus)**2)
1216 ft(6)=(320.0d0+80.0d0*dtanh((betaT-320.0d0)/80.0d0))/
1218 ftprim(6)=1.0d0/(320.0d0*dcosh((betaT-320.0d0)/80.0d0)**2)
1219 ftbis(6)=-2.0d0*dtanh((betaT-320.0d0)/80.0d0)
1220 & /(320.0d0*80.0d0*dcosh((betaT-320.0d0)/80.0d0)**3)
1221 #elif defined(FUNCT)
1230 else if (rescale_mode.eq.0) then
1236 write (iout,*) "Error in WHAM_CALC: wrong RESCALE_MODE",
1243 c--------------------------------------------------------------------
1244 subroutine sum_ene(t,iparm,ft,etot)
1246 include 'DIMENSIONS'
1247 include 'DIMENSIONS.ZSCOPT'
1248 include 'DIMENSIONS.FREE'
1249 include 'COMMON.CONTROL'
1250 include 'COMMON.FFIELD'
1251 include "COMMON.SBRIDGE"
1252 include "COMMON.ENERGIES"
1253 include "COMMON.IOUNITS"
1255 double precision fT(6)
1256 double precision etot,evdw,evdw_t,evdw2,ees,evdw1,ebe,etors,
1257 & escloc,ehpb,ecorr,ecorr5,ecorr6,eello_turn4,eello_turn3,
1258 & eturn6,eel_loc,edihcnstr,etors_d,estr,evdw2_14,esccor,
1260 evdw=enetb(21,t,iparm)
1261 evdw_t=enetb(1,t,iparm)
1263 evdw2_14=enetb(17,t,iparm)
1264 evdw2=enetb(2,t,iparm)+evdw2_14
1266 evdw2=enetb(2,t,iparm)
1270 ees=enetb(3,t,iparm)
1271 evdw1=enetb(16,t,iparm)
1273 ees=enetb(3,t,iparm)
1276 ecorr=enetb(4,t,iparm)
1277 ecorr5=enetb(5,t,iparm)
1278 ecorr6=enetb(6,t,iparm)
1279 eel_loc=enetb(7,t,iparm)
1280 eello_turn3=enetb(8,t,iparm)
1281 eello_turn4=enetb(9,t,iparm)
1282 eturn6=enetb(10,t,iparm)
1283 ebe=enetb(11,t,iparm)
1284 escloc=enetb(12,t,iparm)
1285 etors=enetb(13,t,iparm)
1286 etors_d=enetb(14,t,iparm)
1287 ehpb=enetb(15,t,iparm)
1288 estr=enetb(18,t,iparm)
1289 esccor=enetb(19,t,iparm)
1290 edihcnstr=enetb(20,t,iparm)
1291 eliptran=enetb(22,t,iparm)
1292 esaxs=enetb(26,t,iparm)
1294 if (shield_mode.gt.0) then
1295 etot=ft(1)*wsc*(evdw+ft(6)*evdw_t)+ft(1)*wscp*evdw2
1297 & +ft(1)*wvdwpp*evdw1
1298 & +wang*ebe+ft(1)*wtor*etors+wscloc*escloc
1299 c & +wstrain*ehpb+nss*ebr+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5
1300 & +wstrain*ehpb+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5
1301 & +ft(5)*wcorr6*ecorr6+ft(3)*wturn4*eello_turn4
1302 & +ft(2)*wturn3*eello_turn3
1303 & +ft(5)*wturn6*eturn6+ft(2)*wel_loc*eel_loc
1304 & +edihcnstr+ft(2)*wtor_d*etors_d+ft(1)*wsccor*esccor
1305 & +wbond*estr+wliptran*eliptran+wsaxs*esaxs
1307 etot=wsc*(evdw+ft(6)*evdw_t)+wscp*evdw2+ft(1)*welec*ees
1309 & +wang*ebe+ft(1)*wtor*etors+wscloc*escloc
1310 c & +wstrain*ehpb+nss*ebr+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5
1311 & +wstrain*ehpb+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5
1312 & +ft(5)*wcorr6*ecorr6+ft(3)*wturn4*eello_turn4
1313 & +ft(2)*wturn3*eello_turn3
1314 & +ft(5)*wturn6*eturn6+ft(2)*wel_loc*eel_loc
1315 & +edihcnstr+ft(2)*wtor_d*etors_d+ft(1)*wsccor*esccor
1316 & +wbond*estr+wliptran*eliptran+wsaxs*esaxs
1319 if (shield_mode.gt.0) then
1320 etot=ft(1)*wsc*(evdw+ft(6)*evdw_t)+ft(1)*wscp*evdw2
1321 & +ft(1)*welec*(ees+evdw1)
1322 & +wang*ebe+ft(1)*wtor*etors+wscloc*escloc
1323 c & +wstrain*ehpb+nss*ebr+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5
1324 & +wstrain*ehpb+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5
1325 & +ft(5)*wcorr6*ecorr6+ft(3)*wturn4*eello_turn4
1326 & +ft(2)*wturn3*eello_turn3
1327 & +ft(5)*wturn6*eturn6+ft(2)*wel_loc*eel_loc+edihcnstr
1328 & +ft(2)*wtor_d*etors_d+ft(1)*wsccor*esccor
1329 & +wbond*estr+wliptran*eliptran+wsaxs*esaxs
1331 etot=wsc*(evdw+ft(6)*evdw_t)+wscp*evdw2
1332 & +ft(1)*welec*(ees+evdw1)
1333 & +wang*ebe+ft(1)*wtor*etors+wscloc*escloc
1334 c & +wstrain*ehpb+nss*ebr+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5
1335 & +wstrain*ehpb+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5
1336 & +ft(5)*wcorr6*ecorr6+ft(3)*wturn4*eello_turn4
1337 & +ft(2)*wturn3*eello_turn3
1338 & +ft(5)*wturn6*eturn6+ft(2)*wel_loc*eel_loc+edihcnstr
1339 & +ft(2)*wtor_d*etors_d+ft(1)*wsccor*esccor
1340 & +wbond*estr+wliptran*eliptran+wsaxs*esaxs
1345 c--------------------------------------------------------------------
1346 subroutine sum_ene_deriv(t,iparm,ft,ftprim,ftbis,etot,eprim,ebis)
1348 include 'DIMENSIONS'
1349 include 'DIMENSIONS.ZSCOPT'
1350 include 'DIMENSIONS.FREE'
1351 include 'COMMON.CONTROL'
1352 include 'COMMON.FFIELD'
1353 include "COMMON.SBRIDGE"
1354 include 'COMMON.ENERGIES'
1355 include "COMMON.IOUNITS"
1357 double precision fT(6),fTprim(6),fTbis(6),
1359 double precision etot,evdw,evdw_t,evdw2,ees,evdw1,ebe,etors,
1360 & escloc,ehpb,ecorr,ecorr5,ecorr6,eello_turn4,eello_turn3,
1361 & eturn6,eel_loc,edihcnstr,etors_d,estr,evdw2_14,esccor,
1363 evdw=enetb(21,t,iparm)
1364 evdw_t=enetb(1,t,iparm)
1366 evdw2_14=enetb(17,t,iparm)
1367 evdw2=enetb(2,t,iparm)+evdw2_14
1369 evdw2=enetb(2,t,iparm)
1373 ees=enetb(3,t,iparm)
1374 evdw1=enetb(16,t,iparm)
1376 ees=enetb(3,t,iparm)
1379 ecorr=enetb(4,t,iparm)
1380 ecorr5=enetb(5,t,iparm)
1381 ecorr6=enetb(6,t,iparm)
1382 eel_loc=enetb(7,t,iparm)
1383 eello_turn3=enetb(8,t,iparm)
1384 eello_turn4=enetb(9,t,iparm)
1385 eturn6=enetb(10,t,iparm)
1386 ebe=enetb(11,t,iparm)
1387 escloc=enetb(12,t,iparm)
1388 etors=enetb(13,t,iparm)
1389 etors_d=enetb(14,t,iparm)
1390 ehpb=enetb(15,t,iparm)
1391 estr=enetb(18,t,iparm)
1392 esccor=enetb(19,t,iparm)
1393 edihcnstr=enetb(20,t,iparm)
1394 eliptran=enetb(22,t,iparm)
1395 esaxs=enetb(26,t,iparm)
1397 if (shield_mode.gt.0) then
1398 etot=ft(1)*wsc*(evdw+ft(6)*evdw_t)+ft(1)*wscp*evdw2
1400 & +ft(1)*wvdwpp*evdw1
1401 & +wang*ebe+ft(1)*wtor*etors+wscloc*escloc
1402 c & +wstrain*ehpb+nss*ebr+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5
1403 & +wstrain*ehpb+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5
1404 & +ft(5)*wcorr6*ecorr6+ft(3)*wturn4*eello_turn4
1405 & +ft(2)*wturn3*eello_turn3
1406 & +ft(5)*wturn6*eturn6+ft(2)*wel_loc*eel_loc
1407 & +edihcnstr+ft(2)*wtor_d*etors_d+ft(1)*wsccor*esccor
1408 & +wbond*estr+wliptran*eliptran+wsaxs*esaxs
1409 eprim=ftprim(1)*(ft(6)*evdw_t+evdw)
1410 C & +ftprim(6)*evdw_t
1411 & +ftprim(1)*wscp*evdw2
1412 & +ftprim(1)*welec*ees
1413 & +ftprim(1)*wvdwpp*evdw1
1414 & +ftprim(1)*wtor*etors+
1415 & ftprim(3)*wcorr*ecorr+ftprim(4)*wcorr5*ecorr5+
1416 & ftprim(5)*wcorr6*ecorr6+ftprim(3)*wturn4*eello_turn4+
1417 & ftprim(2)*wturn3*eello_turn3+ftprim(5)*wturn6*eturn6+
1418 & ftprim(2)*wel_loc*eel_loc+ftprim(2)*wtor_d*etors_d+
1419 & ftprim(1)*wsccor*esccor
1420 ebis=ftbis(1)*wsc*(evdw+ft(6)*evdw_t)
1421 & +ftbis(1)*wscp*evdw2+
1422 & ftbis(1)*welec*ees
1423 & +ftbis(1)*wvdwpp*evdw
1424 & +ftbis(1)*wtor*etors+
1425 & ftbis(3)*wcorr*ecorr+ftbis(4)*wcorr5*ecorr5+
1426 & ftbis(5)*wcorr6*ecorr6+ftbis(3)*wturn4*eello_turn4+
1427 & ftbis(2)*wturn3*eello_turn3+ftbis(5)*wturn6*eturn6+
1428 & ftbis(2)*wel_loc*eel_loc+ftbis(2)*wtor_d*etors_d+
1429 & ftbis(1)*wsccor*esccor
1431 etot=wsc*(evdw+ft(6)*evdw_t)+wscp*evdw2+ft(1)*welec*ees
1433 & +wang*ebe+ft(1)*wtor*etors+wscloc*escloc
1434 c & +wstrain*ehpb+nss*ebr+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5
1435 & +wstrain*ehpb+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5
1436 & +ft(5)*wcorr6*ecorr6+ft(3)*wturn4*eello_turn4
1437 & +ft(2)*wturn3*eello_turn3
1438 & +ft(5)*wturn6*eturn6+ft(2)*wel_loc*eel_loc
1439 & +edihcnstr+ft(2)*wtor_d*etors_d+ft(1)*wsccor*esccor
1440 & +wbond*estr+wliptran*eliptran+wsaxs*esaxs
1441 eprim=ftprim(6)*evdw_t+ftprim(1)*welec*ees
1442 & +ftprim(1)*wtor*etors+
1443 & ftprim(3)*wcorr*ecorr+ftprim(4)*wcorr5*ecorr5+
1444 & ftprim(5)*wcorr6*ecorr6+ftprim(3)*wturn4*eello_turn4+
1445 & ftprim(2)*wturn3*eello_turn3+ftprim(5)*wturn6*eturn6+
1446 & ftprim(2)*wel_loc*eel_loc+ftprim(2)*wtor_d*etors_d+
1447 & ftprim(1)*wsccor*esccor
1448 ebis=ftbis(1)*welec*ees+ftbis(1)*wtor*etors+
1449 & ftbis(3)*wcorr*ecorr+ftbis(4)*wcorr5*ecorr5+
1450 & ftbis(5)*wcorr6*ecorr6+ftbis(3)*wturn4*eello_turn4+
1451 & ftbis(2)*wturn3*eello_turn3+ftbis(5)*wturn6*eturn6+
1452 & ftbis(2)*wel_loc*eel_loc+ftbis(2)*wtor_d*etors_d+
1453 & ftbis(1)*wsccor*esccor
1456 if (shield_mode.gt.0) then
1457 etot=ft(1)*wsc*(evdw+ft(6)*evdw_t)+ft(1)*wscp*evdw2
1458 & +ft(1)*welec*(ees+evdw1)
1459 & +wang*ebe+ft(1)*wtor*etors+wscloc*escloc
1460 c & +wstrain*ehpb+nss*ebr+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5
1461 & +wstrain*ehpb+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5
1462 & +ft(5)*wcorr6*ecorr6+ft(3)*wturn4*eello_turn4
1463 & +ft(2)*wturn3*eello_turn3
1464 & +ft(5)*wturn6*eturn6+ft(2)*wel_loc*eel_loc+edihcnstr
1465 & +ft(2)*wtor_d*etors_d+ft(1)*wsccor*esccor
1466 & +wbond*estr+wliptran*eliptran+wsaxs*esaxs
1467 eprim=ftprim(1)*(evdw+ft(6)*evdw_t)
1468 & +ftprim(1)*welec*(ees+evdw1)
1469 & +ftprim(1)*wtor*etors+
1470 & ftprim(1)*wscp*evdw2+
1471 & ftprim(3)*wcorr*ecorr+ftprim(4)*wcorr5*ecorr5+
1472 & ftprim(5)*wcorr6*ecorr6+ftprim(3)*wturn4*eello_turn4+
1473 & ftprim(2)*wturn3*eello_turn3+ftprim(5)*wturn6*eturn6+
1474 & ftprim(2)*wel_loc*eel_loc+ftprim(2)*wtor_d*etors_d+
1475 & ftprim(1)*wsccor*esccor
1476 ebis= ftbis(1)*(evdw+ft(6)*evdw_t)
1477 & +ftbis(1)*wscp*evdw2
1478 & +ftbis(1)*welec*(ees+evdw1)+ftbis(1)*wtor*etors+
1479 & ftbis(3)*wcorr*ecorr+ftbis(4)*wcorr5*ecorr5+
1480 & ftbis(5)*wcorr6*ecorr6+ftbis(3)*wturn4*eello_turn4+
1481 & ftbis(2)*wturn3*eello_turn3+ftbis(5)*wturn6*eturn6+
1482 & ftbis(2)*wel_loc*eel_loc+ftbis(2)*wtor_d*etors_d+
1483 & ftprim(1)*wsccor*esccor
1485 etot=wsc*(evdw+ft(6)*evdw_t)+wscp*evdw2
1486 & +ft(1)*welec*(ees+evdw1)
1487 & +wang*ebe+ft(1)*wtor*etors+wscloc*escloc
1488 c & +wstrain*ehpb+nss*ebr+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5
1489 & +wstrain*ehpb+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5
1490 & +ft(5)*wcorr6*ecorr6+ft(3)*wturn4*eello_turn4
1491 & +ft(2)*wturn3*eello_turn3
1492 & +ft(5)*wturn6*eturn6+ft(2)*wel_loc*eel_loc+edihcnstr
1493 & +ft(2)*wtor_d*etors_d+ft(1)*wsccor*esccor
1494 & +wbond*estr+wliptran*eliptran+wsaxs*esaxs
1495 eprim=ftprim(6)*evdw_t+ftprim(1)*welec*(ees+evdw1)
1496 & +ftprim(1)*wtor*etors+
1497 & ftprim(3)*wcorr*ecorr+ftprim(4)*wcorr5*ecorr5+
1498 & ftprim(5)*wcorr6*ecorr6+ftprim(3)*wturn4*eello_turn4+
1499 & ftprim(2)*wturn3*eello_turn3+ftprim(5)*wturn6*eturn6+
1500 & ftprim(2)*wel_loc*eel_loc+ftprim(2)*wtor_d*etors_d+
1501 & ftprim(1)*wsccor*esccor
1502 ebis=ftbis(1)*welec*(ees+evdw1)+ftbis(1)*wtor*etors+
1503 & ftbis(3)*wcorr*ecorr+ftbis(4)*wcorr5*ecorr5+
1504 & ftbis(5)*wcorr6*ecorr6+ftbis(3)*wturn4*eello_turn4+
1505 & ftbis(2)*wturn3*eello_turn3+ftbis(5)*wturn6*eturn6+
1506 & ftbis(2)*wel_loc*eel_loc+ftbis(2)*wtor_d*etors_d+
1507 & ftprim(1)*wsccor*esccor