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=200000)
23 parameter (MaxHdim=200)
25 parameter (maxinde=200)
29 integer ierror,errcode,status(MPI_STATUS_SIZE)
31 include "COMMON.CONTROL"
32 include "COMMON.IOUNITS"
34 include "COMMON.ENERGIES"
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 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 integer histent_p(0:2000)
68 logical lprint /.true./
70 double precision delta_T /1.0d0/
71 double precision rgymin,rmsmin,rgymax,rmsmax
72 double precision sumW(0:NGridT,Max_Parm),sumE(0:NGridT,Max_Parm),
73 & sumEsq(0:NGridT,Max_Parm),sumQ(MaxQ1,0:NGridT,Max_Parm),
74 & sumQsq(MaxQ1,0:NGridT,Max_Parm),sumEQ(MaxQ1,0:NGridT,Max_Parm),
75 & sumEprim(0:NGridT,Max_Parm),sumEbis(0:NGridT,Max_Parm),betaT,
77 double precision fi(MaxR,maxT_h,Max_Parm),
78 & dd,dd1,dd2,hh,dmin,denom,finorm,avefi,pom,
79 & hfin(0:MaxHdim,maxT_h),histE(0:maxindE),
80 & hrmsrgy(0:MaxBinRgy,0:MaxBinRms,maxT_h),
82 & hfin_ent(0:MaxHdim),vmax,aux
83 double precision fT(6),fTprim(6),fTbis(6),quot,quotl1,quotl,kfacl,
84 & eprim,ebis,temper,kfac/2.4d0/,T0/300.0d0/,startGridT/200.0d0/,
85 & eplus,eminus,logfac,tanhT,tt
86 double precision etot,evdw,evdw_t,evdw2,ees,evdw1,ebe,etors,
87 & escloc,ehpb,ecorr,ecorr5,ecorr6,eello_turn4,eello_turn3,
88 & eturn6,eel_loc,edihcnstr,etors_d,estr,evdw2_14,esccor,
91 integer ind_point(maxpoint),upindE,indE
100 write(licz2,'(bz,i2.2)') islice
102 write (iout,'(//80(1h-)/"Solving WHAM equations for slice",
103 & i2/80(1h-)//)') islice
104 write (iout,*) "delta",delta," nbin1",nbin1
105 write (iout,*) "MaxN",MaxN," MaxQ",MaxQ," MaHdim",MaxHdim
124 if (potE(i,j).le.potEmin) potEmin=potE(i,j)
126 if (q(nQ+1,i).lt.rmsmin) rmsmin=q(nQ+1,i)
127 if (q(nQ+1,i).gt.rmsmax) rmsmax=q(nQ+1,i)
128 if (q(nQ+2,i).lt.rgymin) rgymin=q(nQ+2,i)
129 if (q(nQ+2,i).gt.rgymax) rgymax=q(nQ+2,i)
132 ind=(q(j,i)-dmin+1.0d-8)/delta
134 ind_point(i)=ind_point(i)+ind
136 ind_point(i)=ind_point(i)+nbin1**(j-1)*ind
138 c write (iout,*) "i",i," j",j," q",q(j,i)," ind_point",
141 if (ind_point(i).lt.0 .or. ind_point(i).gt.MaxHdim) then
142 write (iout,*) "Error - index exceeds range for point",i,
143 & " q=",q(j,i)," ind",ind_point(i)
145 write (iout,*) "Processor",me1
147 call MPI_Abort(MPI_COMM_WORLD, Ierror, Errcode )
152 if (ind_point(i).gt.tmax) tmax=ind_point(i)
153 htot(ind_point(i))=htot(ind_point(i))+1
155 write (iout,*) "i",i,"q",(q(j,i),j=1,nQ)," ind",ind_point(i),
156 & " htot",htot(ind_point(i))
163 write (iout,'(a)') "Numbers of counts in Q bins"
165 if (htot(t).gt.0) then
166 write (iout,'(i15,$)') t
169 jj = mod(liczba,nbin1)
171 write (iout,'(i5,$)') jj
173 write (iout,'(i8)') htot(t)
177 write (iout,'(a,i3)') "Number of data points for parameter set",
179 write (iout,'(i7,$)') ((snk(m,ib,iparm,islice),m=1,nr(ib,iparm)),
181 write (iout,'(i8)') stot(islice)
187 call MPI_AllReduce(tmax,tmax_t,1,MPI_INTEGER,MPI_MAX,
190 call MPI_AllReduce(potEmin,potEmin_t,1,MPI_DOUBLE_PRECISION,
191 & MPI_MIN,WHAM_COMM,IERROR)
192 call MPI_AllReduce(rmsmin,rmsmin_t,1,MPI_DOUBLE_PRECISION,
193 & MPI_MIN,WHAM_COMM,IERROR)
194 call MPI_AllReduce(rmsmax,rmsmax_t,1,MPI_DOUBLE_PRECISION,
195 & MPI_MAX,WHAM_COMM,IERROR)
196 call MPI_AllReduce(rgymin,rgymin_t,1,MPI_DOUBLE_PRECISION,
197 & MPI_MIN,WHAM_COMM,IERROR)
198 call MPI_AllReduce(rgymax,rgymax_t,1,MPI_DOUBLE_PRECISION,
199 & MPI_MAX,WHAM_COMM,IERROR)
205 write (iout,*) "potEmin",potEmin
207 rmsmin=deltrms*dint(rmsmin/deltrms)
208 rmsmax=deltrms*dint(rmsmax/deltrms)
209 rgymin=deltrms*dint(rgymin/deltrgy)
210 rgymax=deltrms*dint(rgymax/deltrgy)
211 nbin_rms=(rmsmax-rmsmin)/deltrms
212 nbin_rgy=(rgymax-rgymin)/deltrgy
213 write (iout,*) "rmsmin",rmsmin," rmsmax",rmsmax," rgymin",rgymin,
214 & " rgymax",rgymax," nbin_rms",nbin_rms," nbin_rgy",nbin_rgy
221 write (iout,*) "nFi",nFi
222 ! Compute the Boltzmann factor corresponing to restrain potentials in different
229 c write (9,'(3i5,f10.5)') i,(iparm,potE(i,iparm),iparm=1,nParmSet)
232 write (iout,'(2i5,21f8.2)') i,iparm,
233 & (enetb(k,i,iparm),k=1,22)
235 call restore_parm(iparm)
237 write (iout,*) wsc,wscp,welec,wvdwpp,wang,wtor,wscloc,
238 & wcorr,wcorr5,wcorr6,wturn4,wturn3,wturn6,wel_loc,
239 & wtor_d,wsccor,wbond
242 if (rescale_mode.eq.1) then
243 quot=1.0d0/(beta_h(ib,iparm)*1.987D-3*T0)
250 fT(l)=kfacl/(kfacl-1.0d0+quotl)
253 tt = 1.0d0/(beta_h(ib,iparm)*1.987D-3)
254 ft(6)=(320.0d0+80.0d0*dtanh((tt-320.0d0)/80.0d0))/320.0d0
256 ft(6)=1.0d0/(beta_h(ib,iparm)*1.987D-3*T0)
260 else if (rescale_mode.eq.2) then
261 quot=1.0d0/(T0*beta_h(ib,iparm)*1.987D-3)
265 fT(l)=1.12692801104297249644d0/
266 & dlog(dexp(quotl)+dexp(-quotl))
269 tt = 1.0d0/(beta_h(ib,iparm)*1.987D-3)
270 ft(6)=(320.0d0+80.0d0*dtanh((tt-320.0d0)/80.0d0))/320.0d0
272 ft(6)=1.0d0/(beta_h(ib,iparm)*1.987D-3*T0)
276 c write (iout,*) 1.0d0/(beta_h(ib,iparm)*1.987D-3),ft
277 else if (rescale_mode.eq.0) then
282 write (iout,*) "Error in WHAM_CALC: wrong RESCALE_MODE",
287 evdw=enetb(1,i,iparm)
288 evdw_t=enetb(21,i,iparm)
290 evdw2_14=enetb(17,i,iparm)
291 evdw2=enetb(2,i,iparm)+evdw2_14
293 evdw2=enetb(2,i,iparm)
298 evdw1=enetb(16,i,iparm)
303 ecorr=enetb(4,i,iparm)
304 ecorr5=enetb(5,i,iparm)
305 ecorr6=enetb(6,i,iparm)
306 eel_loc=enetb(7,i,iparm)
307 eello_turn3=enetb(8,i,iparm)
308 eello_turn4=enetb(9,i,iparm)
309 eturn6=enetb(10,i,iparm)
310 ebe=enetb(11,i,iparm)
311 escloc=enetb(12,i,iparm)
312 etors=enetb(13,i,iparm)
313 etors_d=enetb(14,i,iparm)
314 ehpb=enetb(15,i,iparm)
315 estr=enetb(18,i,iparm)
316 esccor=enetb(19,i,iparm)
317 edihcnstr=enetb(20,i,iparm)
318 eliptran=enetb(22,i,iparm)
321 write (iout,'(3i5,6f5.2,14f12.3)') i,ib,iparm,(ft(l),l=1,6),
322 & evdw+evdw_t,evdw2,ees,evdw1,ecorr,eel_loc,estr,ebe,escloc,
323 & etors,etors_d,eello_turn3,eello_turn4,esccor
327 if (shield_mode.gt.0) then
328 etot=ft(1)*wsc*(evdw+ft(6)*evdw_t)+ft(1)*wscp*evdw2
330 & +ft(1)*wvdwpp*evdw1
331 & +wang*ebe+ft(1)*wtor*etors+wscloc*escloc
332 & +wstrain*ehpb+nss*ebr+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5
333 & +ft(5)*wcorr6*ecorr6+ft(3)*wturn4*eello_turn4
334 & +ft(2)*wturn3*eello_turn3
335 & +ft(5)*wturn6*eturn6+ft(2)*wel_loc*eel_loc
336 & +edihcnstr+ft(2)*wtor_d*etors_d+ft(1)*wsccor*esccor
337 & +wbond*estr+wliptran*eliptran
339 etot=wsc*(evdw+ft(6)*evdw_t)+wscp*evdw2+ft(1)*welec*ees
341 & +wang*ebe+ft(1)*wtor*etors+wscloc*escloc
342 & +wstrain*ehpb+nss*ebr+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5
343 & +ft(5)*wcorr6*ecorr6+ft(3)*wturn4*eello_turn4
344 & +ft(2)*wturn3*eello_turn3
345 & +ft(5)*wturn6*eturn6+ft(2)*wel_loc*eel_loc
346 & +edihcnstr+ft(2)*wtor_d*etors_d+ft(1)*wsccor*esccor
347 & +wbond*estr+wliptran*eliptran
350 if (shield_mode.gt.0) then
351 etot=ft(1)*wsc*(evdw+ft(6)*evdw_t)+ft(1)*wscp*evdw2
352 & +ft(1)*welec*(ees+evdw1)
353 & +wang*ebe+ft(1)*wtor*etors+wscloc*escloc
354 & +wstrain*ehpb+nss*ebr+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5
355 & +ft(5)*wcorr6*ecorr6+ft(3)*wturn4*eello_turn4
356 & +ft(2)*wturn3*eello_turn3
357 & +ft(5)*wturn6*eturn6+ft(2)*wel_loc*eel_loc+edihcnstr
358 & +ft(2)*wtor_d*etors_d+ft(1)*wsccor*esccor
359 & +wbond*estr+wliptran*eliptran
361 etot=wsc*(evdw+ft(6)*evdw_t)+wscp*evdw2
362 & +ft(1)*welec*(ees+evdw1)
363 & +wang*ebe+ft(1)*wtor*etors+wscloc*escloc
364 & +wstrain*ehpb+nss*ebr+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5
365 & +ft(5)*wcorr6*ecorr6+ft(3)*wturn4*eello_turn4
366 & +ft(2)*wturn3*eello_turn3
367 & +ft(5)*wturn6*eturn6+ft(2)*wel_loc*eel_loc+edihcnstr
368 & +ft(2)*wtor_d*etors_d+ft(1)*wsccor*esccor
369 & +wbond*estr+wliptran*eliptran
374 write (iout,*) i,iparm,1.0d0/(beta_h(ib,iparm)*1.987D-3),
378 if (iparm.eq.1 .and. ib.eq.1) then
379 write (iout,*)"Conformation",i
382 energia(k)=enetb(k,i,iparm)
384 call enerprint(energia(0),fT)
391 Econstr=Econstr+Kh(j,kk,ib,iparm)
392 & *(dd-q0(j,kk,ib,iparm))**2
395 & -beta_h(ib,iparm)*(etot-potEmin+Econstr)
397 write (iout,'(4i5,4e15.5)') i,kk,ib,iparm,
398 & etot,potEmin,etot-potEmin,v(i,kk,ib,iparm)
404 ! Simple iteration to calculate free energies corresponding to all simulation
408 ! Compute new free-energy values corresponding to the righ-hand side of the
409 ! equation and their derivatives.
410 write (iout,*) "------------------------fi"
420 vf=v(t,l,k,i)+f(l,k,i)
421 if (vf.gt.vmax) vmax=vf
429 aux=f(l,k,i)+v(t,l,k,i)-vmax
431 & denom=denom+snk(l,k,i,islice)*dexp(aux)
435 entfac(t)=-dlog(denom)-vmax
437 write (iout,*) t,"vmax",vmax," denom",denom,"entfac",entfac(t)
442 do ii=1,nR(iib,iparm)
444 fi_p(ii,iib,iparm)=0.0d0
446 fi_p(ii,iib,iparm)=fi_p(ii,iib,iparm)
447 & +dexp(v(t,ii,iib,iparm)+entfac(t))
449 write (iout,'(4i5,3e15.5)') t,ii,iib,iparm,
450 & v(t,ii,iib,iparm),entfac(t),fi_p(ii,iib,iparm)
454 fi(ii,iib,iparm)=0.0d0
456 fi(ii,iib,iparm)=fi(ii,iib,iparm)
457 & +dexp(v(t,ii,iib,iparm)+entfac(t))
466 write (iout,*) "fi before MPI_Reduce me",me,' master',master
468 do ib=1,nT_h(nparmset)
469 write (iout,*) "iparm",iparm," ib",ib
470 write (iout,*) "beta=",beta_h(ib,iparm)
471 write (iout,'(8e15.5)') (fi_p(i,ib,iparm),i=1,nR(ib,iparm))
475 write (iout,*) "REDUCE size",maxR,MaxT_h,nParmSet,
476 & maxR*MaxT_h*nParmSet
477 write (iout,*) "MPI_COMM_WORLD",MPI_COMM_WORLD,
478 & " WHAM_COMM",WHAM_COMM
479 call MPI_Reduce(fi_p(1,1,1),fi(1,1,1),maxR*MaxT_h*nParmSet,
480 & MPI_DOUBLE_PRECISION,
481 & MPI_SUM,Master,WHAM_COMM,IERROR)
483 write (iout,*) "fi after MPI_Reduce nparmset",nparmset
485 write (iout,*) "iparm",iparm
487 write (iout,*) "beta=",beta_h(ib,iparm)
488 write (iout,'(8e15.5)') (fi(i,ib,iparm),i=1,nR(ib,iparm))
492 if (me1.eq.Master) then
498 fi(i,ib,iparm)=-dlog(fi(i,ib,iparm))
499 avefi=avefi+fi(i,ib,iparm)
505 write (iout,*) "Parameter set",iparm
507 write (iout,*) "beta=",beta_h(ib,iparm)
509 fi(i,ib,iparm)=fi(i,ib,iparm)-avefi
511 write (iout,'(8f10.5)') (fi(i,ib,iparm),i=1,nR(ib,iparm))
512 write (iout,'(8f10.5)') (f(i,ib,iparm),i=1,nR(ib,iparm))
516 ! Compute the norm of free-energy increments.
521 finorm=finorm+dabs(fi(i,ib,iparm)-f(i,ib,iparm))
522 f(i,ib,iparm)=fi(i,ib,iparm)
527 write (iout,*) 'Iteration',iter,' finorm',finorm
531 call MPI_Bcast(f(1,1,1),MaxR*MaxT_h*nParmSet,
532 & MPI_DOUBLE_PRECISION,Master,
534 call MPI_Bcast(finorm,1,MPI_DOUBLE_PRECISION,Master,
537 ! Exit, if the increment norm is smaller than pre-assigned tolerance.
538 if (finorm.lt.fimin) then
539 write (iout,*) 'Iteration converged'
546 ! Now, put together the histograms from all simulations, in order to get the
547 ! unbiased total histogram.
557 write (iout,*) "--------------hist"
561 sumW_p(i,iparm)=0.0d0
562 sumE_p(i,iparm)=0.0d0
563 sumEbis_p(i,iparm)=0.0d0
564 sumEsq_p(i,iparm)=0.0d0
566 sumQ_p(j,i,iparm)=0.0d0
567 sumQsq_p(j,i,iparm)=0.0d0
568 sumEQ_p(j,i,iparm)=0.0d0
578 sumEbis(i,iparm)=0.0d0
579 sumEsq(i,iparm)=0.0d0
581 sumQ(j,i,iparm)=0.0d0
582 sumQsq(j,i,iparm)=0.0d0
583 sumEQ(j,i,iparm)=0.0d0
589 c 8/26/05 entropy distribution
594 c ent=-dlog(entfac(t))
596 if (ent.lt.entmin_p) entmin_p=ent
597 if (ent.gt.entmax_p) entmax_p=ent
599 write (iout,*) "entmin",entmin_p," entmax",entmax_p
601 call MPI_Allreduce(entmin_p,entmin,1,MPI_DOUBLE_PRECISION,MPI_MIN,
603 call MPI_Allreduce(entmax_p,entmax,1,MPI_DOUBLE_PRECISION,MPI_MAX,
605 ientmax=entmax-entmin
606 if (ientmax.gt.2000) ientmax=2000
607 write (iout,*) "entmin",entmin," entmax",entmax," ientmax",ientmax
610 c ient=-dlog(entfac(t))-entmin
611 ient=entfac(t)-entmin
612 if (ient.le.2000) histent_p(ient)=histent_p(ient)+1
614 call MPI_Allreduce(histent_p(0),histent(0),ientmax+1,MPI_INTEGER,
615 & MPI_SUM,WHAM_COMM,IERROR)
616 if (me1.eq.Master) then
617 write (iout,*) "Entropy histogram"
619 write(iout,'(f15.4,i10)') entmin+i,histent(i)
627 if (ent.lt.entmin) entmin=ent
628 if (ent.gt.entmax) entmax=ent
630 ientmax=-dlog(entmax)-entmin
631 if (ientmax.gt.2000) ientmax=2000
633 ient=entfac(t)-entmin
634 if (ient.le.2000) histent(ient)=histent(ient)+1
636 write (iout,*) "Entropy histogram"
638 write(iout,'(2f15.4)') entmin+i,histent(i)
643 c write (iout,*) "me1",me1," scount",scount(me1)
669 hrmsrgy(j,i,ib)=0.0d0
671 hrmsrgy_p(j,i,ib)=0.0d0
683 hfin_ent_p(ind)=hfin_ent_p(ind)+dexp(entfac(t))
685 hfin_ent(ind)=hfin_ent(ind)+dexp(entfac(t))
687 c write (iout,'(2i5,20f8.2)') t,t,(enetb(k,t,iparm),k=1,18)
688 call restore_parm(iparm)
689 evdw=enetb(21,t,iparm)
690 evdw_t=enetb(1,t,iparm)
692 evdw2_14=enetb(17,t,iparm)
693 evdw2=enetb(2,t,iparm)+evdw2_14
695 evdw2=enetb(2,t,iparm)
700 evdw1=enetb(16,t,iparm)
705 ecorr=enetb(4,t,iparm)
706 ecorr5=enetb(5,t,iparm)
707 ecorr6=enetb(6,t,iparm)
708 eel_loc=enetb(7,t,iparm)
709 eello_turn3=enetb(8,t,iparm)
710 eello_turn4=enetb(9,t,iparm)
711 eturn6=enetb(10,t,iparm)
712 ebe=enetb(11,t,iparm)
713 escloc=enetb(12,t,iparm)
714 etors=enetb(13,t,iparm)
715 etors_d=enetb(14,t,iparm)
716 ehpb=enetb(15,t,iparm)
717 estr=enetb(18,t,iparm)
718 esccor=enetb(19,t,iparm)
719 edihcnstr=enetb(20,t,iparm)
722 betaT=startGridT+k*delta_T
726 if (rescale_mode.eq.1) then
734 denom=kfacl-1.0d0+quotl
736 ftprim(l)=-l*ft(l)*quotl1/(T0*denom)
737 ftbis(l)=l*kfacl*quotl1*
738 & (2*l*quotl-(l-1)*denom)/(quot*t0*t0*denom**3)
741 ft(6)=(320.0d0+80.0d0*dtanh((betaT-320.0d0)/80.0d0))/
743 ftprim(6)=1.0d0/(320.0d0*dcosh((betaT-320.0d0)/80.0d0)**2)
744 ftbis(6)=-2.0d0*dtanh((betaT-320.0d0)/80.0d0)
745 & /(320.0d0*80.0d0*dcosh((betaT-320.0d0)/80.0d0)**3)
755 else if (rescale_mode.eq.2) then
763 logfac=1.0d0/dlog(eplus+eminus)
764 tanhT=(eplus-eminus)/(eplus+eminus)
765 fT(l)=1.12692801104297249644d0*logfac
766 ftprim(l)=-l*quotl1*ft(l)*tanhT*logfac/T0
767 ftbis(l)=(l-1)*ftprim(l)/(quot*T0)-
768 & 2*l*quotl1/T0*logfac*
769 & (2*l*quotl1*ft(l)/(T0*(eplus+eminus)**2)
773 ft(6)=(320.0d0+80.0d0*dtanh((betaT-320.0d0)/80.0d0))/
775 ftprim(6)=1.0d0/(320.0d0*dcosh((betaT-320.0d0)/80.0d0)**2)
776 ftbis(6)=-2.0d0*dtanh((betaT-320.0d0)/80.0d0)
777 & /(320.0d0*80.0d0*dcosh((betaT-320.0d0)/80.0d0)**3)
787 else if (rescale_mode.eq.0) then
793 write (iout,*) "Error in WHAM_CALC: wrong RESCALE_MODE",
798 c write (iout,*) "ftprim",ftprim
799 c write (iout,*) "ftbis",ftbis
800 betaT=1.0d0/(1.987D-3*betaT)
802 if (shield_mode.gt.0) then
803 etot=ft(1)*wsc*(evdw+ft(6)*evdw_t)+ft(1)*wscp*evdw2
805 & +ft(1)*wvdwpp*evdw1
806 & +wang*ebe+ft(1)*wtor*etors+wscloc*escloc
807 & +wstrain*ehpb+nss*ebr+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5
808 & +ft(5)*wcorr6*ecorr6+ft(3)*wturn4*eello_turn4
809 & +ft(2)*wturn3*eello_turn3
810 & +ft(5)*wturn6*eturn6+ft(2)*wel_loc*eel_loc
811 & +edihcnstr+ft(2)*wtor_d*etors_d+ft(1)*wsccor*esccor
812 & +wbond*estr+wliptran*eliptran
813 eprim=ftprim(1)*(ft(6)*evdw_t+evdw)
814 C & +ftprim(6)*evdw_t
815 & +ftprim(1)*wscp*evdw2
816 & +ftprim(1)*welec*ees
817 & +ftprim(1)*wvdwpp*evdw1
818 & +ftprim(1)*wtor*etors+
819 & ftprim(3)*wcorr*ecorr+ftprim(4)*wcorr5*ecorr5+
820 & ftprim(5)*wcorr6*ecorr6+ftprim(3)*wturn4*eello_turn4+
821 & ftprim(2)*wturn3*eello_turn3+ftprim(5)*wturn6*eturn6+
822 & ftprim(2)*wel_loc*eel_loc+ftprim(2)*wtor_d*etors_d+
823 & ftprim(1)*wsccor*esccor
824 ebis=ftbis(1)*wsc*(evdw+ft(6)*evdw_t)
825 & +ftbis(1)*wscp*evdw2+
827 & +ftbis(1)*wvdwpp*evdw
828 & +ftbis(1)*wtor*etors+
829 & ftbis(3)*wcorr*ecorr+ftbis(4)*wcorr5*ecorr5+
830 & ftbis(5)*wcorr6*ecorr6+ftbis(3)*wturn4*eello_turn4+
831 & ftbis(2)*wturn3*eello_turn3+ftbis(5)*wturn6*eturn6+
832 & ftbis(2)*wel_loc*eel_loc+ftbis(2)*wtor_d*etors_d+
833 & ftbis(1)*wsccor*esccor
835 etot=wsc*(evdw+ft(6)*evdw_t)+wscp*evdw2+ft(1)*welec*ees
837 & +wang*ebe+ft(1)*wtor*etors+wscloc*escloc
838 & +wstrain*ehpb+nss*ebr+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5
839 & +ft(5)*wcorr6*ecorr6+ft(3)*wturn4*eello_turn4
840 & +ft(2)*wturn3*eello_turn3
841 & +ft(5)*wturn6*eturn6+ft(2)*wel_loc*eel_loc
842 & +edihcnstr+ft(2)*wtor_d*etors_d+ft(1)*wsccor*esccor
843 & +wbond*estr+wliptran*eliptran
844 eprim=ftprim(6)*evdw_t+ftprim(1)*welec*ees
845 & +ftprim(1)*wtor*etors+
846 & ftprim(3)*wcorr*ecorr+ftprim(4)*wcorr5*ecorr5+
847 & ftprim(5)*wcorr6*ecorr6+ftprim(3)*wturn4*eello_turn4+
848 & ftprim(2)*wturn3*eello_turn3+ftprim(5)*wturn6*eturn6+
849 & ftprim(2)*wel_loc*eel_loc+ftprim(2)*wtor_d*etors_d+
850 & ftprim(1)*wsccor*esccor
851 ebis=ftbis(1)*welec*ees+ftbis(1)*wtor*etors+
852 & ftbis(3)*wcorr*ecorr+ftbis(4)*wcorr5*ecorr5+
853 & ftbis(5)*wcorr6*ecorr6+ftbis(3)*wturn4*eello_turn4+
854 & ftbis(2)*wturn3*eello_turn3+ftbis(5)*wturn6*eturn6+
855 & ftbis(2)*wel_loc*eel_loc+ftbis(2)*wtor_d*etors_d+
856 & ftbis(1)*wsccor*esccor
859 if (shield_mode.gt.0) then
860 etot=ft(1)*wsc*(evdw+ft(6)*evdw_t)+ft(1)*wscp*evdw2
861 & +ft(1)*welec*(ees+evdw1)
862 & +wang*ebe+ft(1)*wtor*etors+wscloc*escloc
863 & +wstrain*ehpb+nss*ebr+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5
864 & +ft(5)*wcorr6*ecorr6+ft(3)*wturn4*eello_turn4
865 & +ft(2)*wturn3*eello_turn3
866 & +ft(5)*wturn6*eturn6+ft(2)*wel_loc*eel_loc+edihcnstr
867 & +ft(2)*wtor_d*etors_d+ft(1)*wsccor*esccor
868 & +wbond*estr+wliptran*eliptran
869 eprim=ftprim(1)*(evdw+ft(6)*evdw_t)
870 & +ftprim(1)*welec*(ees+evdw1)
871 & +ftprim(1)*wtor*etors+
872 & ftprim(1)*wscp*evdw2+
873 & ftprim(3)*wcorr*ecorr+ftprim(4)*wcorr5*ecorr5+
874 & ftprim(5)*wcorr6*ecorr6+ftprim(3)*wturn4*eello_turn4+
875 & ftprim(2)*wturn3*eello_turn3+ftprim(5)*wturn6*eturn6+
876 & ftprim(2)*wel_loc*eel_loc+ftprim(2)*wtor_d*etors_d+
877 & ftprim(1)*wsccor*esccor
878 ebis= ftbis(1)*(evdw+ft(6)*evdw_t)
879 & +ftbis(1)*wscp*evdw2
880 & +ftbis(1)*welec*(ees+evdw1)+ftbis(1)*wtor*etors+
881 & ftbis(3)*wcorr*ecorr+ftbis(4)*wcorr5*ecorr5+
882 & ftbis(5)*wcorr6*ecorr6+ftbis(3)*wturn4*eello_turn4+
883 & ftbis(2)*wturn3*eello_turn3+ftbis(5)*wturn6*eturn6+
884 & ftbis(2)*wel_loc*eel_loc+ftbis(2)*wtor_d*etors_d+
885 & ftprim(1)*wsccor*esccor
887 etot=wsc*(evdw+ft(6)*evdw_t)+wscp*evdw2
888 & +ft(1)*welec*(ees+evdw1)
889 & +wang*ebe+ft(1)*wtor*etors+wscloc*escloc
890 & +wstrain*ehpb+nss*ebr+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5
891 & +ft(5)*wcorr6*ecorr6+ft(3)*wturn4*eello_turn4
892 & +ft(2)*wturn3*eello_turn3
893 & +ft(5)*wturn6*eturn6+ft(2)*wel_loc*eel_loc+edihcnstr
894 & +ft(2)*wtor_d*etors_d+ft(1)*wsccor*esccor
895 & +wbond*estr+wliptran*eliptran
896 eprim=ftprim(6)*evdw_t+ftprim(1)*welec*(ees+evdw1)
897 & +ftprim(1)*wtor*etors+
898 & ftprim(3)*wcorr*ecorr+ftprim(4)*wcorr5*ecorr5+
899 & ftprim(5)*wcorr6*ecorr6+ftprim(3)*wturn4*eello_turn4+
900 & ftprim(2)*wturn3*eello_turn3+ftprim(5)*wturn6*eturn6+
901 & ftprim(2)*wel_loc*eel_loc+ftprim(2)*wtor_d*etors_d+
902 & ftprim(1)*wsccor*esccor
903 ebis=ftbis(1)*welec*(ees+evdw1)+ftbis(1)*wtor*etors+
904 & ftbis(3)*wcorr*ecorr+ftbis(4)*wcorr5*ecorr5+
905 & ftbis(5)*wcorr6*ecorr6+ftbis(3)*wturn4*eello_turn4+
906 & ftbis(2)*wturn3*eello_turn3+ftbis(5)*wturn6*eturn6+
907 & ftbis(2)*wel_loc*eel_loc+ftbis(2)*wtor_d*etors_d+
908 & ftprim(1)*wsccor*esccor
913 weight=dexp(-betaT*(etot-potEmin)+entfac(t))
915 write (iout,*) "iparm",iparm," t",t," betaT",betaT,
916 & " etot",etot," entfac",entfac(t),
917 & " weight",weight," ebis",ebis
919 etot=etot-temper*eprim
921 sumW_p(k,iparm)=sumW_p(k,iparm)+weight
922 sumE_p(k,iparm)=sumE_p(k,iparm)+etot*weight
923 sumEbis_p(k,iparm)=sumEbis_p(k,iparm)+ebis*weight
924 sumEsq_p(k,iparm)=sumEsq_p(k,iparm)+etot**2*weight
926 sumQ_p(j,k,iparm)=sumQ_p(j,k,iparm)+q(j,t)*weight
927 sumQsq_p(j,k,iparm)=sumQsq_p(j,k,iparm)+q(j,t)**2*weight
928 sumEQ_p(j,k,iparm)=sumEQ_p(j,k,iparm)
929 & +etot*q(j,t)*weight
932 sumW(k,iparm)=sumW(k,iparm)+weight
933 sumE(k,iparm)=sumE(k,iparm)+etot*weight
934 sumEbis(k,iparm)=sumEbis(k,iparm)+ebis*weight
935 sumEsq(k,iparm)=sumEsq(k,iparm)+etot**2*weight
937 sumQ(j,k,iparm)=sumQ(j,k,iparm)+q(j,t)*weight
938 sumQsq(j,k,iparm)=sumQsq(j,k,iparm)+q(j,t)**2*weight
939 sumEQ(j,k,iparm)=sumEQ(j,k,iparm)
940 & +etot*q(j,t)*weight
944 indE = aint(potE(t,iparm)-aint(potEmin))
945 if (indE.ge.0 .and. indE.le.maxinde) then
946 if (indE.gt.upindE_p) upindE_p=indE
947 histE_p(indE)=histE_p(indE)+dexp(-entfac(t))
951 expfac=dexp(-beta_h(ib,iparm)*(etot-potEmin)+entfac(t))
952 hfin_p(ind,ib)=hfin_p(ind,ib)+
953 & dexp(-beta_h(ib,iparm)*(etot-potEmin)+entfac(t))
955 indrgy=dint((q(nQ+2,t)-rgymin)/deltrgy)
956 indrms=dint((q(nQ+1,t)-rmsmin)/deltrms)
957 hrmsrgy_p(indrgy,indrms,ib)=
958 & hrmsrgy_p(indrgy,indrms,ib)+expfac
963 expfac=dexp(-beta_h(ib,iparm)*(etot-potEmin)+entfac(t))
964 hfin(ind,ib)=hfin(ind,ib)+
965 & dexp(-beta_h(ib,iparm)*(etot-potEmin)+entfac(t))
967 indrgy=dint((q(nQ+2,t)-rgymin)/deltrgy)
968 indrms=dint((q(nQ+1,t)-rmsmin)/deltrms)
969 hrmsrgy(indrgy,indrms,ib)=
970 & hrmsrgy(indrgy,indrms,ib)+expfac
976 if (histout) call MPI_Reduce(hfin_p(0,ib),hfin(0,ib),nbin,
977 & MPI_DOUBLE_PRECISION,MPI_SUM,Master,WHAM_COMM,IERROR)
979 call MPI_Reduce(hrmsrgy_p(0,0,ib),hrmsrgy(0,0,ib),
980 & (MaxBinRgy+1)*(nbin_rms+1),MPI_DOUBLE_PRECISION,MPI_SUM,Master,
984 call MPI_Reduce(upindE_p,upindE,1,
985 & MPI_INTEGER,MPI_MAX,Master,WHAM_COMM,IERROR)
986 call MPI_Reduce(histE_p(0),histE(0),maxindE,
987 & MPI_DOUBLE_PRECISION,MPI_SUM,Master,WHAM_COMM,IERROR)
989 if (me1.eq.master) then
993 write (iout,'(6x,$)')
994 write (iout,'(f20.2,$)') (1.0d0/(1.987D-3*beta_h(ib,iparm)),
998 write (iout,'(/a)') 'Final histograms'
1000 if (nslice.eq.1) then
1001 if (separate_parset) then
1002 write(licz3,"(bz,i3.3)") myparm
1003 histname=prefix(:ilen(prefix))//'_par'//licz3//'.hist'
1005 histname=prefix(:ilen(prefix))//'.hist'
1008 if (separate_parset) then
1009 write(licz3,"(bz,i3.3)") myparm
1010 histname=prefix(:ilen(prefix))//'_par'//licz3//
1011 & '_slice_'//licz2//'.hist'
1013 histname=prefix(:ilen(prefix))//'_slice_'//licz2//'.hist'
1016 #if defined(AIX) || defined(PGI)
1017 open (ihist,file=histname,position='append')
1019 open (ihist,file=histname,access='append')
1027 sumH=sumH+hfin(t,ib)
1029 if (sumH.gt.0.0d0) then
1031 jj = mod(liczba,nbin1)
1033 write (iout,'(f6.3,$)') dmin+(jj+0.5d0)*delta
1035 & write (ihist,'(f6.3,$)') dmin+(jj+0.5d0)*delta
1038 write (iout,'(e20.10,$)') hfin(t,ib)
1039 if (histfile) write (ihist,'(e20.10,$)') hfin(t,ib)
1041 write (iout,'(i5)') iparm
1042 if (histfile) write (ihist,'(i5)') iparm
1049 if (nslice.eq.1) then
1050 if (separate_parset) then
1051 write(licz3,"(bz,i3.3)") myparm
1052 histname=prefix(:ilen(prefix))//"_par"//licz3//'.ent'
1054 histname=prefix(:ilen(prefix))//'.ent'
1057 if (separate_parset) then
1058 write(licz3,"(bz,i3.3)") myparm
1059 histname=prefix(:ilen(prefix))//'par_'//licz3//
1060 & '_slice_'//licz2//'.ent'
1062 histname=prefix(:ilen(prefix))//'_slice_'//licz2//'.ent'
1065 #if defined(AIX) || defined(PGI)
1066 open (ihist,file=histname,position='append')
1068 open (ihist,file=histname,access='append')
1070 write (ihist,'(a)') "# Microcanonical entropy"
1072 write (ihist,'(f8.0,$)') dint(potEmin)+i
1073 if (histE(i).gt.0.0e0) then
1074 write (ihist,'(f15.5,$)') dlog(histE(i))
1076 write (ihist,'(f15.5,$)') 0.0d0
1082 write (iout,*) "Microcanonical entropy"
1084 write (iout,'(f8.0,$)') dint(potEmin)+i
1085 if (histE(i).gt.0.0e0) then
1086 write (iout,'(f15.5,$)') dlog(histE(i))
1088 write (iout,'(f15.5,$)') 0.0d0
1093 if (nslice.eq.1) then
1094 if (separate_parset) then
1095 write(licz3,"(bz,i3.3)") myparm
1096 histname=prefix(:ilen(prefix))//'_par'//licz3//'.rmsrgy'
1098 histname=prefix(:ilen(prefix))//'.rmsrgy'
1101 if (separate_parset) then
1102 write(licz3,"(bz,i3.3)") myparm
1103 histname=prefix(:ilen(prefix))//'_par'//licz3//
1104 & '_slice_'//licz2//'.rmsrgy'
1106 histname=prefix(:ilen(prefix))//'_slice_'//licz2//'.rmsrgy'
1109 #if defined(AIX) || defined(PGI)
1110 open (ihist,file=histname,position='append')
1112 open (ihist,file=histname,access='append')
1116 write(ihist,'(2f8.2,$)')
1117 & rgymin+deltrgy*j,rmsmin+deltrms*i
1119 if (hrmsrgy(j,i,ib).gt.0.0d0) then
1120 write(ihist,'(e14.5,$)')
1121 & -dlog(hrmsrgy(j,i,ib))/beta_h(ib,iparm)
1124 write(ihist,'(e14.5,$)') 1.0d6
1127 write (ihist,'(i2)') iparm
1135 call MPI_Reduce(hfin_ent_p(0),hfin_ent(0),nbin,
1136 & MPI_DOUBLE_PRECISION,MPI_SUM,Master,WHAM_COMM,IERROR)
1137 call MPI_Reduce(sumW_p(0,1),sumW(0,1),(nGridT+1)*nParmSet,
1138 & MPI_DOUBLE_PRECISION,MPI_SUM,Master,WHAM_COMM,IERROR)
1139 call MPI_Reduce(sumE_p(0,1),sumE(0,1),(nGridT+1)*nParmSet,
1140 & MPI_DOUBLE_PRECISION,MPI_SUM,Master,WHAM_COMM,IERROR)
1141 call MPI_Reduce(sumEbis_p(0,1),sumEbis(0,1),(nGridT+1)*nParmSet,
1142 & MPI_DOUBLE_PRECISION,MPI_SUM,Master,WHAM_COMM,IERROR)
1143 call MPI_Reduce(sumEsq_p(0,1),sumEsq(0,1),(nGridT+1)*nParmSet,
1144 & MPI_DOUBLE_PRECISION,MPI_SUM,Master,WHAM_COMM,IERROR)
1145 call MPI_Reduce(sumQ_p(1,0,1),sumQ(1,0,1),
1146 & MaxQ1*(nGridT+1)*nParmSet,MPI_DOUBLE_PRECISION,MPI_SUM,Master,
1148 call MPI_Reduce(sumQsq_p(1,0,1),sumQsq(1,0,1),
1149 & MaxQ1*(nGridT+1)*nParmSet,MPI_DOUBLE_PRECISION,MPI_SUM,Master,
1151 call MPI_Reduce(sumEQ_p(1,0,1),sumEQ(1,0,1),
1152 & MaxQ1*(nGridT+1)*nParmSet,MPI_DOUBLE_PRECISION,MPI_SUM,Master,
1154 if (me.eq.master) then
1156 write (iout,'(/a)') 'Thermal characteristics of folding'
1157 if (nslice.eq.1) then
1160 nazwa=prefix(:ilen(prefix))//"_slice_"//licz2
1163 if (nparmset.eq.1 .and. .not.separate_parset) then
1164 nazwa=nazwa(:iln)//".thermal"
1165 else if (nparmset.eq.1 .and. separate_parset) then
1166 write(licz3,"(bz,i3.3)") myparm
1167 nazwa=nazwa(:iln)//"_par_"//licz3//".thermal"
1170 if (nparmset.gt.1) then
1171 write(licz3,"(bz,i3.3)") iparm
1172 nazwa=nazwa(:iln)//"_par_"//licz3//".thermal"
1175 if (separate_parset) then
1176 write (iout,'(a,i3)') "Parameter set",myparm
1178 write (iout,'(a,i3)') "Parameter set",iparm
1181 sumE(i,iparm)=sumE(i,iparm)/sumW(i,iparm)
1182 sumEbis(i,iparm)=(startGridT+i*delta_T)*sumEbis(i,iparm)/
1184 sumEsq(i,iparm)=(sumEsq(i,iparm)/sumW(i,iparm)
1185 & -sumE(i,iparm)**2)/(1.987D-3*(startGridT+i*delta_T)**2)
1187 sumQ(j,i,iparm)=sumQ(j,i,iparm)/sumW(i,iparm)
1188 sumQsq(j,i,iparm)=sumQsq(j,i,iparm)/sumW(i,iparm)
1189 & -sumQ(j,i,iparm)**2
1190 sumEQ(j,i,iparm)=sumEQ(j,i,iparm)/sumW(i,iparm)
1191 & -sumQ(j,i,iparm)*sumE(i,iparm)
1193 sumW(i,iparm)=-dlog(sumW(i,iparm))*(1.987D-3*
1194 & (startGridT+i*delta_T))+potEmin
1195 write (iout,'(f7.1,2f15.5,$)') startGridT+i*delta_T,
1196 & sumW(i,iparm),sumE(i,iparm)
1197 write (iout,'(f10.5,$)') (sumQ(j,i,iparm),j=1,nQ+2)
1198 write (iout,'(e15.5,$)') sumEsq(i,iparm)-sumEbis(i,iparm),
1199 & (sumQsq(j,i,iparm),j=1,nQ+2),(sumEQ(j,i,iparm),j=1,nQ+2)
1201 write (34,'(f7.1,2f15.5,$)') startGridT+i*delta_T,
1202 & sumW(i,iparm),sumE(i,iparm)
1203 write (34,'(f10.5,$)') (sumQ(j,i,iparm),j=1,nQ+2)
1204 write (34,'(e15.5,$)') sumEsq(i,iparm)-sumEbis(i,iparm),
1205 & (sumQsq(j,i,iparm),j=1,nQ+2),(sumEQ(j,i,iparm),j=1,nQ+2)
1213 if (hfin_ent(t).gt.0.0d0) then
1215 jj = mod(liczba,nbin1)
1216 write (iout,'(f6.3,e20.10," ent")') dmin+(jj+0.5d0)*delta,
1218 if (histfile) write (ihist,'(f6.3,e20.10," ent")')
1219 & dmin+(jj+0.5d0)*delta,
1223 if (histfile) close(ihist)
1227 ! Write data for zscore
1228 if (nslice.eq.1) then
1229 zscname=prefix(:ilen(prefix))//".zsc"
1231 zscname=prefix(:ilen(prefix))//"_slice_"//licz2//".zsc"
1233 #if defined(AIX) || defined(PGI)
1234 open (izsc,file=prefix(:ilen(prefix))//'.zsc',position='append')
1236 open (izsc,file=prefix(:ilen(prefix))//'.zsc',access='append')
1238 write (izsc,'("NQ=",i1," NPARM=",i1)') nQ,nParmSet
1240 write (izsc,'("NT=",i1)') nT_h(iparm)
1242 write (izsc,'("TEMP=",f6.1," NR=",i2," SNK=",$)')
1243 & 1.0d0/(beta_h(ib,iparm)*1.987D-3),nR(ib,iparm)
1244 jj = min0(nR(ib,iparm),7)
1245 write (izsc,'(i8,$)') (snk(i,ib,iparm,islice),i=1,jj)
1246 write (izsc,'(a1,$)') (" ",i=22+8*jj+1,79)
1247 write (izsc,'("&")')
1248 if (nR(ib,iparm).gt.7) then
1249 do ii=8,nR(ib,iparm),9
1250 jj = min0(nR(ib,iparm),ii+8)
1251 write (izsc,'(i8,$)') (snk(i,ib,iparm,islice),i=ii,jj)
1252 write (izsc,'(a1,$') (" ",i=(jj-ii+1)*8+1,79)
1253 write (izsc,'("&")')
1256 write (izsc,'("FI=",$)')
1257 jj=min0(nR(ib,iparm),7)
1258 write (izsc,'(f10.5,$)') (fi(i,ib,iparm),i=1,jj)
1259 write (izsc,'(a1,$)') (" ",i=3+10*jj+1,79)
1260 write (izsc,'("&")')
1261 if (nR(ib,iparm).gt.7) then
1262 do ii=8,nR(ib,iparm),9
1263 jj = min0(nR(ib,iparm),ii+8)
1264 write (izsc,'(f10.5,$)') (fi(i,ib,iparm),i=ii,jj)
1265 if (jj.eq.nR(ib,iparm)) then
1268 write (izsc,'(a1,$)') (" ",i=10*(jj-ii+1)+1,79)
1269 write (izsc,'(t80,"&")')
1274 write (izsc,'("KH=",$)')
1275 write (izsc,'(f7.2,$)') (Kh(j,i,ib,iparm),j=1,nQ)
1276 write (izsc,'(" Q0=",$)')
1277 write (izsc,'(f7.5,$)') (q0(j,i,ib,iparm),j=1,nQ)