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)
326 write (iout,'(3i5,6f5.2,14f12.3)') i,ib,iparm,(ft(l),l=1,6),
327 & evdw+evdw_t,evdw2,ees,evdw1,ecorr,eel_loc,estr,ebe,escloc,
328 & etors,etors_d,eello_turn3,eello_turn4,esccor
332 if (shield_mode.gt.0) then
333 etot=ft(1)*wsc*(evdw+ft(6)*evdw_t)+ft(1)*wscp*evdw2
335 & +ft(1)*wvdwpp*evdw1
336 & +wang*ebe+ft(1)*wtor*etors+wscloc*escloc
337 ! & +wstrain*ehpb+nss*ebr+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5
338 & +wstrain*ehpb+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5
339 & +ft(5)*wcorr6*ecorr6+ft(3)*wturn4*eello_turn4
340 & +ft(2)*wturn3*eello_turn3
341 & +ft(5)*wturn6*eturn6+ft(2)*wel_loc*eel_loc
342 & +edihcnstr+ft(2)*wtor_d*etors_d+ft(1)*wsccor*esccor
343 & +wbond*estr+wliptran*eliptran
345 etot=wsc*(evdw+ft(6)*evdw_t)+wscp*evdw2+ft(1)*welec*ees
347 & +wang*ebe+ft(1)*wtor*etors+wscloc*escloc
348 c & +wstrain*ehpb+nss*ebr+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5
349 & +wstrain*ehpb+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5
350 & +ft(5)*wcorr6*ecorr6+ft(3)*wturn4*eello_turn4
351 & +ft(2)*wturn3*eello_turn3
352 & +ft(5)*wturn6*eturn6+ft(2)*wel_loc*eel_loc
353 & +edihcnstr+ft(2)*wtor_d*etors_d+ft(1)*wsccor*esccor
354 & +wbond*estr+wliptran*eliptran
357 if (shield_mode.gt.0) then
358 etot=ft(1)*wsc*(evdw+ft(6)*evdw_t)+ft(1)*wscp*evdw2
359 & +ft(1)*welec*(ees+evdw1)
360 & +wang*ebe+ft(1)*wtor*etors+wscloc*escloc
361 c & +wstrain*ehpb+nss*ebr+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5
362 & +wstrain*ehpb+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5
363 & +ft(5)*wcorr6*ecorr6+ft(3)*wturn4*eello_turn4
364 & +ft(2)*wturn3*eello_turn3
365 & +ft(5)*wturn6*eturn6+ft(2)*wel_loc*eel_loc+edihcnstr
366 & +ft(2)*wtor_d*etors_d+ft(1)*wsccor*esccor
367 & +wbond*estr+wliptran*eliptran
369 etot=wsc*(evdw+ft(6)*evdw_t)+wscp*evdw2
370 & +ft(1)*welec*(ees+evdw1)
371 & +wang*ebe+ft(1)*wtor*etors+wscloc*escloc
372 c & +wstrain*ehpb+nss*ebr+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5
373 & +wstrain*ehpb+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5
374 & +ft(5)*wcorr6*ecorr6+ft(3)*wturn4*eello_turn4
375 & +ft(2)*wturn3*eello_turn3
376 & +ft(5)*wturn6*eturn6+ft(2)*wel_loc*eel_loc+edihcnstr
377 & +ft(2)*wtor_d*etors_d+ft(1)*wsccor*esccor
378 & +wbond*estr+wliptran*eliptran
383 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)
402 Econstr=Econstr+Kh(j,kk,ib,iparm)
403 & *(dd-q0(j,kk,ib,iparm))**2
405 c Adaptive potential contribution
407 call PMF_energy(q(1,i),ib,kk,iparm,ePMF,ePMF_q)
411 & -beta_h(ib,iparm)*(etot-potEmin+Econstr)
413 write (iout,'(4i5,4e15.5)') i,kk,ib,iparm,
414 & etot,potEmin,etot-potEmin,v(i,kk,ib,iparm)
420 ! Simple iteration to calculate free energies corresponding to all simulation
424 ! Compute new free-energy values corresponding to the righ-hand side of the
425 ! equation and their derivatives.
426 write (iout,*) "------------------------fi"
436 vf=v(t,l,k,i)+f(l,k,i)
437 if (vf.gt.vmax) vmax=vf
445 aux=f(l,k,i)+v(t,l,k,i)-vmax
447 & denom=denom+snk(l,k,i,islice)*dexp(aux)
451 entfac(t)=-dlog(denom)-vmax
453 write (iout,*) t,"vmax",vmax," denom",denom,"entfac",entfac(t)
459 do ii=1,nR(iib,iparm)
461 fimax_p(ii,iib,iparm)=v(1,ii,iib,iparm)+entfac(1)
463 if(v(t,ii,iib,iparm)+entfac(t).gt.fimax_p(ii,iib,iparm))
464 & fimax_p(ii,iib,iparm)=v(t,ii,iib,iparm)+entfac(t)
467 fimax(ii,iib,iparm)=v(1,ii,iib,iparm)+entfac(1)
469 if(v(t,ii,iib,iparm)+entfac(t).gt.fimax(ii,iib,iparm))
470 & fimax(ii,iib,iparm)=v(t,ii,iib,iparm)+entfac(t)
477 call MPI_AllReduce(fimax_p(1,1,1),fimax(1,1,1),
478 & maxR*MaxT_h*nParmSet,MPI_DOUBLE_PRECISION,
479 & MPI_MAX,WHAM_COMM,IERROR)
483 do ii=1,nR(iib,iparm)
485 fi_p(ii,iib,iparm)=0.0d0
487 fi_p(ii,iib,iparm)=fi_p(ii,iib,iparm)
488 & +dexp(v(t,ii,iib,iparm)+entfac(t)-fimax(ii,iib,iparm))
490 write (iout,'(4i5,3e15.5)') t,ii,iib,iparm,
491 & v(t,ii,iib,iparm),entfac(t),fi_p(ii,iib,iparm)
495 fi(ii,iib,iparm)=0.0d0
497 fi(ii,iib,iparm)=fi(ii,iib,iparm)
498 & +dexp(v(t,ii,iib,iparm)+entfac(t)-fimax(ii,iib,iparm))
507 write (iout,*) "fi before MPI_Reduce me",me,' master',master
509 do ib=1,nT_h(nparmset)
510 write (iout,*) "iparm",iparm," ib",ib
511 write (iout,*) "beta=",beta_h(ib,iparm)
512 write (iout,'(8e15.5)') (fi_p(i,ib,iparm),i=1,nR(ib,iparm))
516 c write (iout,*) "REDUCE size",maxR,MaxT_h,nParmSet,
517 c & maxR*MaxT_h*nParmSet
518 c write (iout,*) "MPI_COMM_WORLD",MPI_COMM_WORLD,
519 c & " WHAM_COMM",WHAM_COMM
520 call MPI_Reduce(fi_p(1,1,1),fi(1,1,1),maxR*MaxT_h*nParmSet,
521 & MPI_DOUBLE_PRECISION,
522 & MPI_SUM,Master,WHAM_COMM,IERROR)
524 write (iout,*) "fi after MPI_Reduce nparmset",nparmset
526 write (iout,*) "iparm",iparm
528 write (iout,*) "beta=",beta_h(ib,iparm)
529 write (iout,'(8e15.5)') (fi(i,ib,iparm),i=1,nR(ib,iparm))
533 if (me1.eq.Master) then
539 fi(i,ib,iparm)=-dlog(fi(i,ib,iparm))-fimax(i,ib,iparm)
540 avefi=avefi+fi(i,ib,iparm)
546 write (iout,*) "Parameter set",iparm
548 write (iout,*) "beta=",beta_h(ib,iparm)
550 fi(i,ib,iparm)=fi(i,ib,iparm)-avefi
552 write (iout,'(8f10.5)') (fi(i,ib,iparm),i=1,nR(ib,iparm))
553 write (iout,'(8f10.5)') (f(i,ib,iparm),i=1,nR(ib,iparm))
557 ! Compute the norm of free-energy increments.
562 finorm=finorm+dabs(fi(i,ib,iparm)-f(i,ib,iparm))
563 f(i,ib,iparm)=fi(i,ib,iparm)
568 write (iout,*) 'Iteration',iter,' finorm',finorm
572 call MPI_Bcast(f(1,1,1),MaxR*MaxT_h*nParmSet,
573 & MPI_DOUBLE_PRECISION,Master,
575 call MPI_Bcast(finorm,1,MPI_DOUBLE_PRECISION,Master,
578 ! Exit, if the increment norm is smaller than pre-assigned tolerance.
579 if (finorm.lt.fimin) then
580 write (iout,*) 'Iteration converged'
587 ! Now, put together the histograms from all simulations, in order to get the
588 ! unbiased total histogram.
598 write (iout,*) "--------------hist"
602 sumW_p(i,iparm)=0.0d0
603 sumE_p(i,iparm)=0.0d0
604 sumEbis_p(i,iparm)=0.0d0
605 sumEsq_p(i,iparm)=0.0d0
607 sumQ_p(j,i,iparm)=0.0d0
608 sumQsq_p(j,i,iparm)=0.0d0
609 sumEQ_p(j,i,iparm)=0.0d0
619 sumEbis(i,iparm)=0.0d0
620 sumEsq(i,iparm)=0.0d0
622 sumQ(j,i,iparm)=0.0d0
623 sumQsq(j,i,iparm)=0.0d0
624 sumEQ(j,i,iparm)=0.0d0
630 c 8/26/05 entropy distribution
635 c ent=-dlog(entfac(t))
637 if (ent.lt.entmin_p) entmin_p=ent
638 if (ent.gt.entmax_p) entmax_p=ent
640 write (iout,*) "entmin",entmin_p," entmax",entmax_p
642 call MPI_Allreduce(entmin_p,entmin,1,MPI_DOUBLE_PRECISION,MPI_MIN,
644 call MPI_Allreduce(entmax_p,entmax,1,MPI_DOUBLE_PRECISION,MPI_MAX,
646 ientmax=entmax-entmin
647 if (ientmax.gt.2000) ientmax=2000
648 write (iout,*) "entmin",entmin," entmax",entmax," ientmax",ientmax
651 c ient=-dlog(entfac(t))-entmin
652 ient=entfac(t)-entmin
653 if (ient.le.2000) histent_p(ient)=histent_p(ient)+1
655 call MPI_Allreduce(histent_p(0),histent(0),ientmax+1,MPI_INTEGER,
656 & MPI_SUM,WHAM_COMM,IERROR)
657 if (me1.eq.Master) then
658 write (iout,*) "Entropy histogram"
660 write(iout,'(f15.4,i10)') entmin+i,histent(i)
668 if (ent.lt.entmin) entmin=ent
669 if (ent.gt.entmax) entmax=ent
671 ientmax=-dlog(entmax)-entmin
672 if (ientmax.gt.2000) ientmax=2000
674 ient=entfac(t)-entmin
675 if (ient.le.2000) histent(ient)=histent(ient)+1
677 write (iout,*) "Entropy histogram"
679 write(iout,'(2f15.4)') entmin+i,histent(i)
684 call restore_parm(iparm)
710 hrmsrgy(j,i,ib)=0.0d0
712 hrmsrgy_p(j,i,ib)=0.0d0
722 indE = aint(potE(t,iparm)-aint(potEmin))
723 if (indE.ge.0 .and. indE.le.maxinde) then
724 if (indE.gt.upindE_p) upindE_p=indE
725 histE_p(indE)=histE_p(indE)+dexp(-entfac(t))
729 expfac=dexp(-beta_h(ib,iparm)*(etot-potEmin)+entfac(t))
730 hfin_p(ind,ib)=hfin_p(ind,ib)+
731 & dexp(-beta_h(ib,iparm)*(etot-potEmin)+entfac(t))
733 indrgy=dint((q(nQ+2,t)-rgymin)/deltrgy)
734 indrms=dint((q(nQ+1,t)-rmsmin)/deltrms)
735 hrmsrgy_p(indrgy,indrms,ib)=
736 & hrmsrgy_p(indrgy,indrms,ib)+expfac
741 expfac=dexp(-beta_h(ib,iparm)*(etot-potEmin)+entfac(t))
742 hfin(ind,ib)=hfin(ind,ib)+
743 & dexp(-beta_h(ib,iparm)*(etot-potEmin)+entfac(t))
745 indrgy=dint((q(nQ+2,t)-rgymin)/deltrgy)
746 indrms=dint((q(nQ+1,t)-rmsmin)/deltrms)
747 hrmsrgy(indrgy,indrms,ib)=
748 & hrmsrgy(indrgy,indrms,ib)+expfac
754 c Thermo and ensemble averages
757 betaT=startGridT+k*delta_T
758 call temp_scalfac(betaT,ft,ftprim,ftbis,*10)
759 c write (iout,*) "ftprim",ftprim
760 c write (iout,*) "ftbis",ftbis
761 betaT=1.0d0/(1.987D-3*betaT)
762 c 7/10/18 AL Determine the max Botzmann weights for each temerature
763 call sum_ene(1,iparm,ft,etot)
764 weimax(k,iparm)=-betaT*(etot-potEmin)+entfac(1)
765 c write (iout,*) "k",k," t",1," weight",weimax(k,iparm)
771 call sum_ene(t,iparm,ft,etot)
772 weight=-betaT*(etot-potEmin)+entfac(t)
773 c write (iout,*) "k",k," t",t," weight",weight
774 if (weight.gt.weimax(k,iparm)) weimax(k,iparm)=weight
779 write (iout,*) "weimax before REDUCE"
780 write (iout,*) (weimax(k,iparm),k=0,ngridt)
783 weimax_(k)=weimax(k,iparm)
785 call MPI_Allreduce(weimax_(0),weimax(0,iparm),nGridT+1,
786 & MPI_DOUBLE_PRECISION,MPI_MAX,WHAM_COMM,IERROR)
788 write (iout,*) "weimax"
789 write (iout,*) (weimax(k,iparm),k=0,ngridt)
792 temper=startGridT+k*delta_T
793 betaT=1.0d0/(1.987D-3*temper)
794 call temp_scalfac(temper,ft,ftprim,ftbis,*10)
801 hfin_ent_p(ind)=hfin_ent_p(ind)+dexp(entfac(t))
803 hfin_ent(ind)=hfin_ent(ind)+dexp(entfac(t))
805 c write (iout,'(2i5,20f8.2)') t,t,(enetb(k,t,iparm),k=1,18)
806 c call restore_parm(iparm)
807 call sum_ene_deriv(t,iparm,ft,ftprim,ftbis,etot,eprim,ebis)
808 weight=dexp(-betaT*(etot-potEmin)+entfac(t)-weimax(k,iparm))
810 write (iout,*) "iparm",iparm," t",t," betaT",betaT,
811 & " etot",etot," entfac",entfac(t)," boltz",
812 & -betaT*(etot-potEmin)+entfac(t)," weimax",weimax(k,iparm),
813 & " weight",weight," ebis",ebis
815 etot=etot-temper*eprim
817 sumW_p(k,iparm)=sumW_p(k,iparm)+weight
818 sumE_p(k,iparm)=sumE_p(k,iparm)+etot*weight
819 sumEbis_p(k,iparm)=sumEbis_p(k,iparm)+ebis*weight
820 sumEsq_p(k,iparm)=sumEsq_p(k,iparm)+etot**2*weight
822 sumQ_p(j,k,iparm)=sumQ_p(j,k,iparm)+q(j,t)*weight
823 sumQsq_p(j,k,iparm)=sumQsq_p(j,k,iparm)+q(j,t)**2*weight
824 sumEQ_p(j,k,iparm)=sumEQ_p(j,k,iparm)
825 & +etot*q(j,t)*weight
828 sumW(k,iparm)=sumW(k,iparm)+weight
829 sumE(k,iparm)=sumE(k,iparm)+etot*weight
830 sumEbis(k,iparm)=sumEbis(k,iparm)+ebis*weight
831 sumEsq(k,iparm)=sumEsq(k,iparm)+etot**2*weight
833 sumQ(j,k,iparm)=sumQ(j,k,iparm)+q(j,t)*weight
834 sumQsq(j,k,iparm)=sumQsq(j,k,iparm)+q(j,t)**2*weight
835 sumEQ(j,k,iparm)=sumEQ(j,k,iparm)
836 & +etot*q(j,t)*weight
842 if (histout) call MPI_Reduce(hfin_p(0,ib),hfin(0,ib),nbin,
843 & MPI_DOUBLE_PRECISION,MPI_SUM,Master,WHAM_COMM,IERROR)
845 call MPI_Reduce(hrmsrgy_p(0,0,ib),hrmsrgy(0,0,ib),
846 & (MaxBinRgy+1)*(nbin_rms+1),MPI_DOUBLE_PRECISION,MPI_SUM,Master,
850 call MPI_Reduce(upindE_p,upindE,1,
851 & MPI_INTEGER,MPI_MAX,Master,WHAM_COMM,IERROR)
852 call MPI_Reduce(histE_p(0),histE(0),maxindE,
853 & MPI_DOUBLE_PRECISION,MPI_SUM,Master,WHAM_COMM,IERROR)
855 if (me1.eq.master) then
859 write (iout,'(6x,$)')
860 write (iout,'(f20.2,$)') (1.0d0/(1.987D-3*beta_h(ib,iparm)),
864 write (iout,'(/a)') 'Final histograms'
866 if (nslice.eq.1) then
867 if (separate_parset) then
868 write(licz3,"(bz,i3.3)") myparm
869 histname=prefix(:ilen(prefix))//'_par'//licz3//'.hist'
871 histname=prefix(:ilen(prefix))//'.hist'
874 if (separate_parset) then
875 write(licz3,"(bz,i3.3)") myparm
876 histname=prefix(:ilen(prefix))//'_par'//licz3//
877 & '_slice_'//licz2//'.hist'
879 histname=prefix(:ilen(prefix))//'_slice_'//licz2//'.hist'
882 #if defined(AIX) || defined(PGI)
883 open (ihist,file=histname,position='append')
885 open (ihist,file=histname,access='append')
895 if (sumH.gt.0.0d0) then
897 jj = mod(liczba,nbin1)
899 write (iout,'(f6.3,$)') dmin+(jj+0.5d0)*delta
901 & write (ihist,'(f6.3,$)') dmin+(jj+0.5d0)*delta
904 write (iout,'(e20.10,$)') hfin(t,ib)
905 if (histfile) write (ihist,'(e20.10,$)') hfin(t,ib)
907 write (iout,'(i5)') iparm
908 if (histfile) write (ihist,'(i5)') iparm
915 if (nslice.eq.1) then
916 if (separate_parset) then
917 write(licz3,"(bz,i3.3)") myparm
918 histname=prefix(:ilen(prefix))//"_par"//licz3//'.ent'
920 histname=prefix(:ilen(prefix))//'.ent'
923 if (separate_parset) then
924 write(licz3,"(bz,i3.3)") myparm
925 histname=prefix(:ilen(prefix))//'par_'//licz3//
926 & '_slice_'//licz2//'.ent'
928 histname=prefix(:ilen(prefix))//'_slice_'//licz2//'.ent'
931 #if defined(AIX) || defined(PGI)
932 open (ihist,file=histname,position='append')
934 open (ihist,file=histname,access='append')
936 write (ihist,'(a)') "# Microcanonical entropy"
938 write (ihist,'(f8.0,$)') dint(potEmin)+i
939 if (histE(i).gt.0.0e0) then
940 write (ihist,'(f15.5,$)') dlog(histE(i))
942 write (ihist,'(f15.5,$)') 0.0d0
948 write (iout,*) "Microcanonical entropy"
950 write (iout,'(f8.0,$)') dint(potEmin)+i
951 if (histE(i).gt.0.0e0) then
952 write (iout,'(f15.5,$)') dlog(histE(i))
954 write (iout,'(f15.5,$)') 0.0d0
959 if (nslice.eq.1) then
960 if (separate_parset) then
961 write(licz3,"(bz,i3.3)") myparm
962 histname=prefix(:ilen(prefix))//'_par'//licz3//'.rmsrgy'
964 histname=prefix(:ilen(prefix))//'.rmsrgy'
967 if (separate_parset) then
968 write(licz3,"(bz,i3.3)") myparm
969 histname=prefix(:ilen(prefix))//'_par'//licz3//
970 & '_slice_'//licz2//'.rmsrgy'
972 histname=prefix(:ilen(prefix))//'_slice_'//licz2//'.rmsrgy'
975 #if defined(AIX) || defined(PGI)
976 open (ihist,file=histname,position='append')
978 open (ihist,file=histname,access='append')
982 write(ihist,'(2f8.2,$)')
983 & rgymin+deltrgy*j,rmsmin+deltrms*i
985 if (hrmsrgy(j,i,ib).gt.0.0d0) then
986 write(ihist,'(e14.5,$)')
987 & -dlog(hrmsrgy(j,i,ib))/beta_h(ib,iparm)
990 write(ihist,'(e14.5,$)') 1.0d6
993 write (ihist,'(i2)') iparm
1001 call MPI_Reduce(hfin_ent_p(0),hfin_ent(0),nbin,
1002 & MPI_DOUBLE_PRECISION,MPI_SUM,Master,WHAM_COMM,IERROR)
1003 call MPI_Reduce(sumW_p(0,1),sumW(0,1),(nGridT+1)*nParmSet,
1004 & MPI_DOUBLE_PRECISION,MPI_SUM,Master,WHAM_COMM,IERROR)
1005 call MPI_Reduce(sumE_p(0,1),sumE(0,1),(nGridT+1)*nParmSet,
1006 & MPI_DOUBLE_PRECISION,MPI_SUM,Master,WHAM_COMM,IERROR)
1007 call MPI_Reduce(sumEbis_p(0,1),sumEbis(0,1),(nGridT+1)*nParmSet,
1008 & MPI_DOUBLE_PRECISION,MPI_SUM,Master,WHAM_COMM,IERROR)
1009 call MPI_Reduce(sumEsq_p(0,1),sumEsq(0,1),(nGridT+1)*nParmSet,
1010 & MPI_DOUBLE_PRECISION,MPI_SUM,Master,WHAM_COMM,IERROR)
1011 call MPI_Reduce(sumQ_p(1,0,1),sumQ(1,0,1),
1012 & MaxQ1*(nGridT+1)*nParmSet,MPI_DOUBLE_PRECISION,MPI_SUM,Master,
1014 call MPI_Reduce(sumQsq_p(1,0,1),sumQsq(1,0,1),
1015 & MaxQ1*(nGridT+1)*nParmSet,MPI_DOUBLE_PRECISION,MPI_SUM,Master,
1017 call MPI_Reduce(sumEQ_p(1,0,1),sumEQ(1,0,1),
1018 & MaxQ1*(nGridT+1)*nParmSet,MPI_DOUBLE_PRECISION,MPI_SUM,Master,
1020 if (me.eq.master) then
1022 write (iout,'(/a)') 'Thermal characteristics of folding'
1023 if (nslice.eq.1) then
1026 nazwa=prefix(:ilen(prefix))//"_slice_"//licz2
1029 if (nparmset.eq.1 .and. .not.separate_parset) then
1030 nazwa=nazwa(:iln)//".thermal"
1031 else if (nparmset.eq.1 .and. separate_parset) then
1032 write(licz3,"(bz,i3.3)") myparm
1033 nazwa=nazwa(:iln)//"_par_"//licz3//".thermal"
1036 if (nparmset.gt.1) then
1037 write(licz3,"(bz,i3.3)") iparm
1038 nazwa=nazwa(:iln)//"_par_"//licz3//".thermal"
1041 if (separate_parset) then
1042 write (iout,'(a,i3)') "Parameter set",myparm
1044 write (iout,'(a,i3)') "Parameter set",iparm
1047 sumE(i,iparm)=sumE(i,iparm)/sumW(i,iparm)
1048 sumEbis(i,iparm)=(startGridT+i*delta_T)*sumEbis(i,iparm)/
1050 sumEsq(i,iparm)=(sumEsq(i,iparm)/sumW(i,iparm)
1051 & -sumE(i,iparm)**2)/(1.987D-3*(startGridT+i*delta_T)**2)
1053 sumQ(j,i,iparm)=sumQ(j,i,iparm)/sumW(i,iparm)
1054 sumQsq(j,i,iparm)=sumQsq(j,i,iparm)/sumW(i,iparm)
1055 & -sumQ(j,i,iparm)**2
1056 sumEQ(j,i,iparm)=sumEQ(j,i,iparm)/sumW(i,iparm)
1057 & -sumQ(j,i,iparm)*sumE(i,iparm)
1059 sumW(i,iparm)=(-dlog(sumW(i,iparm))-weimax(i,iparm))*(1.987D-3*
1060 & (startGridT+i*delta_T))+potEmin
1061 write (iout,'(f7.1,2f15.5,$)') startGridT+i*delta_T,
1062 & sumW(i,iparm),sumE(i,iparm)
1063 write (iout,'(f10.5,$)') (sumQ(j,i,iparm),j=1,nQ+2)
1064 write (iout,'(e15.5,$)') sumEsq(i,iparm)-sumEbis(i,iparm),
1065 & (sumQsq(j,i,iparm),j=1,nQ+2),(sumEQ(j,i,iparm),j=1,nQ+2)
1067 write (34,'(f7.1,2f15.5,$)') startGridT+i*delta_T,
1068 & sumW(i,iparm),sumE(i,iparm)
1069 write (34,'(f10.5,$)') (sumQ(j,i,iparm),j=1,nQ+2)
1070 write (34,'(e15.5,$)') sumEsq(i,iparm)-sumEbis(i,iparm),
1071 & (sumQsq(j,i,iparm),j=1,nQ+2),(sumEQ(j,i,iparm),j=1,nQ+2)
1079 if (hfin_ent(t).gt.0.0d0) then
1081 jj = mod(liczba,nbin1)
1082 write (iout,'(f6.3,e20.10," ent")') dmin+(jj+0.5d0)*delta,
1084 if (histfile) write (ihist,'(f6.3,e20.10," ent")')
1085 & dmin+(jj+0.5d0)*delta,
1089 if (histfile) close(ihist)
1093 ! Write data for zscore
1094 if (nslice.eq.1) then
1095 zscname=prefix(:ilen(prefix))//".zsc"
1097 zscname=prefix(:ilen(prefix))//"_slice_"//licz2//".zsc"
1099 #if defined(AIX) || defined(PGI)
1100 open (izsc,file=prefix(:ilen(prefix))//'.zsc',position='append')
1102 open (izsc,file=prefix(:ilen(prefix))//'.zsc',access='append')
1104 write (izsc,'("NQ=",i1," NPARM=",i1)') nQ,nParmSet
1106 write (izsc,'("NT=",i1)') nT_h(iparm)
1108 write (izsc,'("TEMP=",f6.1," NR=",i2," SNK=",$)')
1109 & 1.0d0/(beta_h(ib,iparm)*1.987D-3),nR(ib,iparm)
1110 jj = min0(nR(ib,iparm),7)
1111 write (izsc,'(i8,$)') (snk(i,ib,iparm,islice),i=1,jj)
1112 write (izsc,'(a1,$)') (" ",i=22+8*jj+1,79)
1113 write (izsc,'("&")')
1114 if (nR(ib,iparm).gt.7) then
1115 do ii=8,nR(ib,iparm),9
1116 jj = min0(nR(ib,iparm),ii+8)
1117 write (izsc,'(i8,$)') (snk(i,ib,iparm,islice),i=ii,jj)
1118 write (izsc,'(a1,$') (" ",i=(jj-ii+1)*8+1,79)
1119 write (izsc,'("&")')
1122 write (izsc,'("FI=",$)')
1123 jj=min0(nR(ib,iparm),7)
1124 write (izsc,'(f10.5,$)') (fi(i,ib,iparm),i=1,jj)
1125 write (izsc,'(a1,$)') (" ",i=3+10*jj+1,79)
1126 write (izsc,'("&")')
1127 if (nR(ib,iparm).gt.7) then
1128 do ii=8,nR(ib,iparm),9
1129 jj = min0(nR(ib,iparm),ii+8)
1130 write (izsc,'(f10.5,$)') (fi(i,ib,iparm),i=ii,jj)
1131 if (jj.eq.nR(ib,iparm)) then
1134 write (izsc,'(a1,$)') (" ",i=10*(jj-ii+1)+1,79)
1135 write (izsc,'(t80,"&")')
1140 write (izsc,'("KH=",$)')
1141 write (izsc,'(f7.2,$)') (Kh(j,i,ib,iparm),j=1,nQ)
1142 write (izsc,'(" Q0=",$)')
1143 write (izsc,'(f7.5,$)') (q0(j,i,ib,iparm),j=1,nQ)
1158 c------------------------------------------------------------------------
1159 subroutine temp_scalfac(betaT,ft,ftprim,ftbis,*)
1161 include "DIMENSIONS"
1162 include "DIMENSIONS.FREE"
1163 include "COMMON.CONTROL"
1164 include "COMMON.FREE"
1165 include "COMMON.IOUNITS"
1166 double precision fT(6),fTprim(6),fTbis(6),quot,quotl1,quotl,kfacl,
1167 & kfac/2.4d0/,T0/300.0d0/,startGridT/200.0d0/,
1168 & logfac,tanhT,betaT,denom,eplus,eminus
1170 if (rescale_mode.eq.1) then
1178 denom=kfacl-1.0d0+quotl
1180 ftprim(l)=-l*ft(l)*quotl1/(T0*denom)
1181 ftbis(l)=l*kfacl*quotl1*
1182 & (2*l*quotl-(l-1)*denom)/(quot*t0*t0*denom**3)
1185 ft(6)=(320.0d0+80.0d0*dtanh((betaT-320.0d0)/80.0d0))/
1187 ftprim(6)=1.0d0/(320.0d0*dcosh((betaT-320.0d0)/80.0d0)**2)
1188 ftbis(6)=-2.0d0*dtanh((betaT-320.0d0)/80.0d0)
1189 & /(320.0d0*80.0d0*dcosh((betaT-320.0d0)/80.0d0)**3)
1190 #elif defined(FUNCT)
1199 else if (rescale_mode.eq.2) then
1207 logfac=1.0d0/dlog(eplus+eminus)
1208 tanhT=(eplus-eminus)/(eplus+eminus)
1209 fT(l)=1.12692801104297249644d0*logfac
1210 ftprim(l)=-l*quotl1*ft(l)*tanhT*logfac/T0
1211 ftbis(l)=(l-1)*ftprim(l)/(quot*T0)-
1212 & 2*l*quotl1/T0*logfac*
1213 & (2*l*quotl1*ft(l)/(T0*(eplus+eminus)**2)
1217 ft(6)=(320.0d0+80.0d0*dtanh((betaT-320.0d0)/80.0d0))/
1219 ftprim(6)=1.0d0/(320.0d0*dcosh((betaT-320.0d0)/80.0d0)**2)
1220 ftbis(6)=-2.0d0*dtanh((betaT-320.0d0)/80.0d0)
1221 & /(320.0d0*80.0d0*dcosh((betaT-320.0d0)/80.0d0)**3)
1222 #elif defined(FUNCT)
1231 else if (rescale_mode.eq.0) then
1237 write (iout,*) "Error in WHAM_CALC: wrong RESCALE_MODE",
1244 c--------------------------------------------------------------------
1245 subroutine sum_ene(t,iparm,ft,etot)
1247 include 'DIMENSIONS'
1248 include 'DIMENSIONS.ZSCOPT'
1249 include 'DIMENSIONS.FREE'
1250 include 'COMMON.CONTROL'
1251 include 'COMMON.FFIELD'
1252 include "COMMON.SBRIDGE"
1253 include "COMMON.ENERGIES"
1254 include "COMMON.IOUNITS"
1256 double precision fT(6)
1257 double precision etot,evdw,evdw_t,evdw2,ees,evdw1,ebe,etors,
1258 & escloc,ehpb,ecorr,ecorr5,ecorr6,eello_turn4,eello_turn3,
1259 & eturn6,eel_loc,edihcnstr,etors_d,estr,evdw2_14,esccor,
1261 evdw=enetb(21,t,iparm)
1262 evdw_t=enetb(1,t,iparm)
1264 evdw2_14=enetb(17,t,iparm)
1265 evdw2=enetb(2,t,iparm)+evdw2_14
1267 evdw2=enetb(2,t,iparm)
1271 ees=enetb(3,t,iparm)
1272 evdw1=enetb(16,t,iparm)
1274 ees=enetb(3,t,iparm)
1277 ecorr=enetb(4,t,iparm)
1278 ecorr5=enetb(5,t,iparm)
1279 ecorr6=enetb(6,t,iparm)
1280 eel_loc=enetb(7,t,iparm)
1281 eello_turn3=enetb(8,t,iparm)
1282 eello_turn4=enetb(9,t,iparm)
1283 eturn6=enetb(10,t,iparm)
1284 ebe=enetb(11,t,iparm)
1285 escloc=enetb(12,t,iparm)
1286 etors=enetb(13,t,iparm)
1287 etors_d=enetb(14,t,iparm)
1288 ehpb=enetb(15,t,iparm)
1289 estr=enetb(18,t,iparm)
1290 esccor=enetb(19,t,iparm)
1291 edihcnstr=enetb(20,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
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
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
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
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)
1396 if (shield_mode.gt.0) then
1397 etot=ft(1)*wsc*(evdw+ft(6)*evdw_t)+ft(1)*wscp*evdw2
1399 & +ft(1)*wvdwpp*evdw1
1400 & +wang*ebe+ft(1)*wtor*etors+wscloc*escloc
1401 c & +wstrain*ehpb+nss*ebr+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5
1402 & +wstrain*ehpb+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5
1403 & +ft(5)*wcorr6*ecorr6+ft(3)*wturn4*eello_turn4
1404 & +ft(2)*wturn3*eello_turn3
1405 & +ft(5)*wturn6*eturn6+ft(2)*wel_loc*eel_loc
1406 & +edihcnstr+ft(2)*wtor_d*etors_d+ft(1)*wsccor*esccor
1407 & +wbond*estr+wliptran*eliptran
1408 eprim=ftprim(1)*(ft(6)*evdw_t+evdw)
1409 C & +ftprim(6)*evdw_t
1410 & +ftprim(1)*wscp*evdw2
1411 & +ftprim(1)*welec*ees
1412 & +ftprim(1)*wvdwpp*evdw1
1413 & +ftprim(1)*wtor*etors+
1414 & ftprim(3)*wcorr*ecorr+ftprim(4)*wcorr5*ecorr5+
1415 & ftprim(5)*wcorr6*ecorr6+ftprim(3)*wturn4*eello_turn4+
1416 & ftprim(2)*wturn3*eello_turn3+ftprim(5)*wturn6*eturn6+
1417 & ftprim(2)*wel_loc*eel_loc+ftprim(2)*wtor_d*etors_d+
1418 & ftprim(1)*wsccor*esccor
1419 ebis=ftbis(1)*wsc*(evdw+ft(6)*evdw_t)
1420 & +ftbis(1)*wscp*evdw2+
1421 & ftbis(1)*welec*ees
1422 & +ftbis(1)*wvdwpp*evdw
1423 & +ftbis(1)*wtor*etors+
1424 & ftbis(3)*wcorr*ecorr+ftbis(4)*wcorr5*ecorr5+
1425 & ftbis(5)*wcorr6*ecorr6+ftbis(3)*wturn4*eello_turn4+
1426 & ftbis(2)*wturn3*eello_turn3+ftbis(5)*wturn6*eturn6+
1427 & ftbis(2)*wel_loc*eel_loc+ftbis(2)*wtor_d*etors_d+
1428 & ftbis(1)*wsccor*esccor
1430 etot=wsc*(evdw+ft(6)*evdw_t)+wscp*evdw2+ft(1)*welec*ees
1432 & +wang*ebe+ft(1)*wtor*etors+wscloc*escloc
1433 c & +wstrain*ehpb+nss*ebr+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5
1434 & +wstrain*ehpb+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5
1435 & +ft(5)*wcorr6*ecorr6+ft(3)*wturn4*eello_turn4
1436 & +ft(2)*wturn3*eello_turn3
1437 & +ft(5)*wturn6*eturn6+ft(2)*wel_loc*eel_loc
1438 & +edihcnstr+ft(2)*wtor_d*etors_d+ft(1)*wsccor*esccor
1439 & +wbond*estr+wliptran*eliptran
1440 eprim=ftprim(6)*evdw_t+ftprim(1)*welec*ees
1441 & +ftprim(1)*wtor*etors+
1442 & ftprim(3)*wcorr*ecorr+ftprim(4)*wcorr5*ecorr5+
1443 & ftprim(5)*wcorr6*ecorr6+ftprim(3)*wturn4*eello_turn4+
1444 & ftprim(2)*wturn3*eello_turn3+ftprim(5)*wturn6*eturn6+
1445 & ftprim(2)*wel_loc*eel_loc+ftprim(2)*wtor_d*etors_d+
1446 & ftprim(1)*wsccor*esccor
1447 ebis=ftbis(1)*welec*ees+ftbis(1)*wtor*etors+
1448 & ftbis(3)*wcorr*ecorr+ftbis(4)*wcorr5*ecorr5+
1449 & ftbis(5)*wcorr6*ecorr6+ftbis(3)*wturn4*eello_turn4+
1450 & ftbis(2)*wturn3*eello_turn3+ftbis(5)*wturn6*eturn6+
1451 & ftbis(2)*wel_loc*eel_loc+ftbis(2)*wtor_d*etors_d+
1452 & ftbis(1)*wsccor*esccor
1455 if (shield_mode.gt.0) then
1456 etot=ft(1)*wsc*(evdw+ft(6)*evdw_t)+ft(1)*wscp*evdw2
1457 & +ft(1)*welec*(ees+evdw1)
1458 & +wang*ebe+ft(1)*wtor*etors+wscloc*escloc
1459 c & +wstrain*ehpb+nss*ebr+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5
1460 & +wstrain*ehpb+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5
1461 & +ft(5)*wcorr6*ecorr6+ft(3)*wturn4*eello_turn4
1462 & +ft(2)*wturn3*eello_turn3
1463 & +ft(5)*wturn6*eturn6+ft(2)*wel_loc*eel_loc+edihcnstr
1464 & +ft(2)*wtor_d*etors_d+ft(1)*wsccor*esccor
1465 & +wbond*estr+wliptran*eliptran
1466 eprim=ftprim(1)*(evdw+ft(6)*evdw_t)
1467 & +ftprim(1)*welec*(ees+evdw1)
1468 & +ftprim(1)*wtor*etors+
1469 & ftprim(1)*wscp*evdw2+
1470 & ftprim(3)*wcorr*ecorr+ftprim(4)*wcorr5*ecorr5+
1471 & ftprim(5)*wcorr6*ecorr6+ftprim(3)*wturn4*eello_turn4+
1472 & ftprim(2)*wturn3*eello_turn3+ftprim(5)*wturn6*eturn6+
1473 & ftprim(2)*wel_loc*eel_loc+ftprim(2)*wtor_d*etors_d+
1474 & ftprim(1)*wsccor*esccor
1475 ebis= ftbis(1)*(evdw+ft(6)*evdw_t)
1476 & +ftbis(1)*wscp*evdw2
1477 & +ftbis(1)*welec*(ees+evdw1)+ftbis(1)*wtor*etors+
1478 & ftbis(3)*wcorr*ecorr+ftbis(4)*wcorr5*ecorr5+
1479 & ftbis(5)*wcorr6*ecorr6+ftbis(3)*wturn4*eello_turn4+
1480 & ftbis(2)*wturn3*eello_turn3+ftbis(5)*wturn6*eturn6+
1481 & ftbis(2)*wel_loc*eel_loc+ftbis(2)*wtor_d*etors_d+
1482 & ftprim(1)*wsccor*esccor
1484 etot=wsc*(evdw+ft(6)*evdw_t)+wscp*evdw2
1485 & +ft(1)*welec*(ees+evdw1)
1486 & +wang*ebe+ft(1)*wtor*etors+wscloc*escloc
1487 c & +wstrain*ehpb+nss*ebr+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5
1488 & +wstrain*ehpb+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5
1489 & +ft(5)*wcorr6*ecorr6+ft(3)*wturn4*eello_turn4
1490 & +ft(2)*wturn3*eello_turn3
1491 & +ft(5)*wturn6*eturn6+ft(2)*wel_loc*eel_loc+edihcnstr
1492 & +ft(2)*wtor_d*etors_d+ft(1)*wsccor*esccor
1493 & +wbond*estr+wliptran*eliptran
1494 eprim=ftprim(6)*evdw_t+ftprim(1)*welec*(ees+evdw1)
1495 & +ftprim(1)*wtor*etors+
1496 & ftprim(3)*wcorr*ecorr+ftprim(4)*wcorr5*ecorr5+
1497 & ftprim(5)*wcorr6*ecorr6+ftprim(3)*wturn4*eello_turn4+
1498 & ftprim(2)*wturn3*eello_turn3+ftprim(5)*wturn6*eturn6+
1499 & ftprim(2)*wel_loc*eel_loc+ftprim(2)*wtor_d*etors_d+
1500 & ftprim(1)*wsccor*esccor
1501 ebis=ftbis(1)*welec*(ees+evdw1)+ftbis(1)*wtor*etors+
1502 & ftbis(3)*wcorr*ecorr+ftbis(4)*wcorr5*ecorr5+
1503 & ftbis(5)*wcorr6*ecorr6+ftbis(3)*wturn4*eello_turn4+
1504 & ftbis(2)*wturn3*eello_turn3+ftbis(5)*wturn6*eturn6+
1505 & ftbis(2)*wel_loc*eel_loc+ftbis(2)*wtor_d*etors_d+
1506 & ftprim(1)*wsccor*esccor