2 !-----------------------------------------------------------------------------
12 !-----------------------------------------------------------------------------
15 !-----------------------------------------------------------------------------
17 !-----------------------------------------------------------------------------
19 subroutine WHAMCALC(islice,*)
20 ! Weighed Histogram Analysis Method (WHAM) code
21 ! Written by A. Liwo based on the work of Kumar et al.,
22 ! J.Comput.Chem., 13, 1011 (1992)
24 ! 2/1/05 Multiple temperatures allowed.
25 ! 2/2/05 Free energies calculated directly from data points
26 ! acc. to Eq. (21) of Kumar et al.; final histograms also
27 ! constructed based on this equation.
28 ! 2/12/05 Multiple parameter sets included
30 ! 2/2/05 Parallel version
34 use io_base, only:ilen
39 ! include "DIMENSIONS"
40 ! include "DIMENSIONS.ZSCOPT"
41 ! include "DIMENSIONS.FREE"
42 integer,parameter :: NGridT=400
43 integer,parameter :: MaxBinRms=100,MaxBinRgy=100
44 integer,parameter :: MaxHdim=200
45 ! parameter (MaxHdim=200000)
46 integer,parameter :: maxinde=200
48 integer :: ierror,errcode,status(MPI_STATUS_SIZE)
50 ! include "COMMON.CONTROL"
51 ! include "COMMON.IOUNITS"
52 ! include "COMMON.FREE"
53 ! include "COMMON.ENERGIES"
54 ! include "COMMON.FFIELD"
55 ! include "COMMON.SBRIDGE"
56 ! include "COMMON.PROT"
57 ! include "COMMON.ENEPS"
58 integer,parameter :: MaxPoint=MaxStr,&
59 MaxPointProc=MaxStr_Proc
60 real(kind=8),parameter :: finorm_max=1.0d0
61 real(kind=8) :: potfac,expfac,vf
62 ! real(kind=8) :: potfac,entmin,entmax,expfac,vf
64 integer :: i,ii,j,jj,k,kk,l,m,ind,iter,t,tmax,ient,ientmax,iln
65 integer :: start,end,iharm,ib,iib,nbin1,nbin,nbin_rms,nbin_rgy,&
66 nbin_rmsrgy,liczbaW,iparm,nFi,indrgy,indrms
67 ! 4/17/17 AKS & AL: histent is obsolete
68 integer :: htot(0:MaxHdim)!,histent(0:2000)
69 real(kind=8) :: v(MaxPointProc,MaxR,MaxT_h,nParmSet) !(MaxPointProc,MaxR,MaxT_h,Max_Parm)
70 real(kind=8) :: energia(0:n_ene)
71 !el real(kind=8) :: energia(0:max_ene)
73 integer :: tmax_t,upindE_p
74 real(kind=8) :: fi_p(MaxR,MaxT_h,nParmSet) !(MaxR,MaxT_h,Max_Parm)
75 real(kind=8),dimension(0:nGridT,nParmSet) :: sumW_p,sumE_p,&
76 sumEbis_p,sumEsq_p !(0:nGridT,Max_Parm)
77 real(kind=8),dimension(MaxQ1,0:nGridT,nParmSet) :: sumQ_p,&
78 sumQsq_p,sumEQ_p,sumEprim_p !(MaxQ1,0:nGridT,Max_Parm)
79 real(kind=8) :: hfin_p(0:MaxHdim,maxT_h),&
80 hfin_ent_p(0:MaxHdim),histE_p(0:maxindE),sumH,&
81 hrmsrgy_p(0:MaxBinRgy,0:MaxBinRms,maxT_h)
82 real(kind=8) :: rgymin_t,rmsmin_t,rgymax_t,rmsmax_t
83 real(kind=8) :: potEmin_t!,entmin_p,entmax_p
84 ! integer :: histent_p(0:2000)
85 logical :: lprint=.true.
87 real(kind=8) :: delta_T=1.0d0,iientmax
88 real(kind=8) :: rgymin,rmsmin,rgymax,rmsmax
89 real(kind=8),dimension(0:nGridT,nParmSet) :: sumW,sumE,&
90 sumEsq,sumEprim,sumEbis !(0:NGridT,Max_Parm)
91 real(kind=8),dimension(MaxQ1,0:nGridT,nParmSet) :: sumQ,&
92 sumQsq,sumEQ !(MaxQ1,0:NGridT,Max_Parm)
93 real(kind=8) :: betaT,weight,econstr
94 real(kind=8) :: fi(MaxR,MaxT_h,nParmSet),& !(MaxR,maxT_h,Max_Parm)
95 ddW,dd1,dd2,hh,dmin,denom,finorm,avefi,pom,&
96 hfin(0:MaxHdim,maxT_h),histE(0:maxindE),&
97 hrmsrgy(0:MaxBinRgy,0:MaxBinRms,maxT_h),&
99 hfin_ent(0:MaxHdim),vmax,aux
100 real(kind=8) :: fT(6),fTprim(6),fTbis(6),quot,quotl1,quotl,kfacl,&
101 eprim,ebis,temper,kfac=2.4d0,T0=300.0d0,startGridT=200.0d0,&
102 eplus,eminus,logfac,tanhT,tt
103 real(kind=8) :: etot,evdw,evdw_t,evdw2,ees,evdw1,ebe,etors,&
104 escloc,ehpb,ecorr,ecorr5,ecorr6,eello_turn4,eello_turn3,&
105 eello_turn6,eel_loc,edihcnstr,etors_d,estr,evdw2_14,esccor, &
106 ecationcation,ecation_prot
108 integer :: ind_point(maxpoint),upindE,indE
109 character(len=16) :: plik
110 character(len=1) :: licz1
111 character(len=2) :: licz2
112 character(len=3) :: licz3
113 character(len=128) :: nazwa
118 write(licz2,'(bz,i2.2)') islice
120 write (iout,'(//80(1h-)/"Solving WHAM equations for slice",&
121 i2/80(1h-)//)') islice
122 write (iout,*) "delta",delta," nbin1",nbin1
123 write (iout,*) "MaxN",MaxN," MaxQ",MaxQ," MaHdim",MaxHdim
141 if (potE(i,j).le.potEmin) potEmin=potE(i,j)
143 if (q(nQ+1,i).lt.rmsmin) rmsmin=q(nQ+1,i)
144 if (q(nQ+1,i).gt.rmsmax) rmsmax=q(nQ+1,i)
145 if (q(nQ+2,i).lt.rgymin) rgymin=q(nQ+2,i)
146 if (q(nQ+2,i).gt.rgymax) rgymax=q(nQ+2,i)
149 ind=(q(j,i)-dmin+1.0d-8)/delta
151 ind_point(i)=ind_point(i)+ind
153 ind_point(i)=ind_point(i)+nbin1**(j-1)*ind
155 ! write (iout,*) "i",i," j",j," q",q(j,i)," ind_point",
158 if (ind_point(i).lt.0 .or. ind_point(i).gt.MaxHdim) then
159 write (iout,*) "Error - index exceeds range for point",i,&
160 " q=",q(j,i)," ind",ind_point(i)
162 write (iout,*) "Processor",me1
164 call MPI_Abort(MPI_COMM_WORLD, Ierror, Errcode )
169 if (ind_point(i).gt.tmax) tmax=ind_point(i)
170 htot(ind_point(i))=htot(ind_point(i))+1
172 write (iout,*) "i",i,"q",(q(j,i),j=1,nQ)," ind",ind_point(i),&
173 " htot",htot(ind_point(i))
180 write (iout,'(a)') "Numbers of counts in Q bins"
182 if (htot(t).gt.0) then
183 write (iout,'(i15,$)') t
186 jj = mod(liczbaW,nbin1)
187 liczbaW=liczbaW/nbin1
188 write (iout,'(i5,$)') jj
190 write (iout,'(i8)') htot(t)
194 write (iout,'(a,i3)') "Number of data points for parameter set",&
196 write (iout,'(i7,$)') ((snk(m,ib,iparm,islice),m=1,nr(ib,iparm)),&
198 write (iout,'(i8)') stot(islice)
204 call MPI_AllReduce(tmax,tmax_t,1,MPI_INTEGER,MPI_MAX,&
207 call MPI_AllReduce(potEmin,potEmin_t,1,MPI_DOUBLE_PRECISION,&
208 MPI_MIN,WHAM_COMM,IERROR)
209 call MPI_AllReduce(rmsmin,rmsmin_t,1,MPI_DOUBLE_PRECISION,&
210 MPI_MIN,WHAM_COMM,IERROR)
211 call MPI_AllReduce(rmsmax,rmsmax_t,1,MPI_DOUBLE_PRECISION,&
212 MPI_MAX,WHAM_COMM,IERROR)
213 call MPI_AllReduce(rgymin,rgymin_t,1,MPI_DOUBLE_PRECISION,&
214 MPI_MIN,WHAM_COMM,IERROR)
215 call MPI_AllReduce(rgymax,rgymax_t,1,MPI_DOUBLE_PRECISION,&
216 MPI_MAX,WHAM_COMM,IERROR)
222 write (iout,*) "potEmin",potEmin
224 rmsmin=deltrms*dint(rmsmin/deltrms)
225 rmsmax=deltrms*dint(rmsmax/deltrms)
226 rgymin=deltrms*dint(rgymin/deltrgy)
227 rgymax=deltrms*dint(rgymax/deltrgy)
228 nbin_rms=(rmsmax-rmsmin)/deltrms
229 nbin_rgy=(rgymax-rgymin)/deltrgy
230 write (iout,*) "rmsmin",rmsmin," rmsmax",rmsmax," rgymin",rgymin,&
231 " rgymax",rgymax," nbin_rms",nbin_rms," nbin_rgy",nbin_rgy
238 write (iout,*) "nFi",nFi
239 ! Compute the Boltzmann factor corresponing to restrain potentials in different
246 ! write (9,'(3i5,f10.5)') i,(iparm,potE(i,iparm),iparm=1,nParmSet)
249 write (iout,'(2i5,21f8.2)') i,iparm,&
250 (enetb(k,i,iparm),k=1,21)
252 call restore_parm(iparm)
254 write (iout,*) wsc,wscp,welec,wvdwpp,wang,wtor,wscloc,&
255 wcorr,wcorr5,wcorr6,wturn4,wturn3,wturn6,wel_loc,&
259 !el old rascale weights
261 ! if (rescale_modeW.eq.1) then
262 ! quot=1.0d0/(beta_h(ib,iparm)*1.987D-3*T0)
269 ! fT(l)=kfacl/(kfacl-1.0d0+quotl)
272 ! tt = 1.0d0/(beta_h(ib,iparm)*1.987D-3)
273 ! ft(6)=(320.0d0+80.0d0*dtanh((tt-320.0d0)/80.0d0))/320.0d0
274 !#elif defined(FUNCT)
275 ! ft(6)=1.0d0/(beta_h(ib,iparm)*1.987D-3*T0)
279 ! else if (rescale_modeW.eq.2) then
280 ! quot=1.0d0/(T0*beta_h(ib,iparm)*1.987D-3)
284 ! fT(l)=1.12692801104297249644d0/ &
285 ! dlog(dexp(quotl)+dexp(-quotl))
288 ! tt = 1.0d0/(beta_h(ib,iparm)*1.987D-3)
289 ! ft(6)=(320.0d0+80.0d0*dtanh((tt-320.0d0)/80.0d0))/320.0d0
290 !#elif defined(FUNCT)
291 ! ft(6)=1.0d0/(beta_h(ib,iparm)*1.987D-3*T0)
295 ! write (iout,*) 1.0d0/(beta_h(ib,iparm)*1.987D-3),ft
296 ! else if (rescale_modeW.eq.0) then
301 ! write (iout,*) "Error in WHAM_CALC: wrong RESCALE_MODE",&
306 ! el end old rescale weights
307 call rescale_weights(1.0d0/(beta_h(ib,iparm)*1.987D-3))
309 ! call etot(enetb(0,i,iparm))
310 evdw=enetb(1,i,iparm)
311 ! evdw_t=enetb(21,i,iparm)
312 evdw_t=enetb(20,i,iparm)
314 ! evdw2_14=enetb(17,i,iparm)
315 evdw2_14=enetb(18,i,iparm)
316 evdw2=enetb(2,i,iparm)+evdw2_14
318 evdw2=enetb(2,i,iparm)
323 evdw1=enetb(16,i,iparm)
328 ecorr=enetb(4,i,iparm)
329 ecorr5=enetb(5,i,iparm)
330 ecorr6=enetb(6,i,iparm)
331 eel_loc=enetb(7,i,iparm)
332 eello_turn3=enetb(8,i,iparm)
333 eello_turn4=enetb(9,i,iparm)
334 eello_turn6=enetb(10,i,iparm)
335 ebe=enetb(11,i,iparm)
336 escloc=enetb(12,i,iparm)
337 etors=enetb(13,i,iparm)
338 etors_d=enetb(14,i,iparm)
339 ehpb=enetb(15,i,iparm)
340 ! estr=enetb(18,i,iparm)
341 estr=enetb(17,i,iparm)
342 ! esccor=enetb(19,i,iparm)
343 esccor=enetb(21,i,iparm)
344 ! edihcnstr=enetb(20,i,iparm)
345 edihcnstr=enetb(19,i,iparm)
346 ecationcation=enetb(42,i,iparm)
347 ecation_prot=enetb(41,i,iparm)
349 write (iout,'(3i5,6f5.2,14f12.3)') i,ib,iparm,(ft(l),l=1,6),&
350 evdw+evdw_t,evdw2,ees,evdw1,ecorr,eel_loc,estr,ebe,escloc,&
351 etors,etors_d,eello_turn3,eello_turn4,esccor
355 ! etot=wsc*(evdw+ft(6)*evdw_t)+wscp*evdw2+ft(1)*welec*ees &
357 ! +wang*ebe+ft(1)*wtor*etors+wscloc*escloc &
358 ! +wstrain*ehpb+nss*ebr+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5 &
359 ! +ft(5)*wcorr6*ecorr6+ft(3)*wturn4*eello_turn4 &
360 ! +ft(2)*wturn3*eello_turn3 &
361 ! +ft(5)*wturn6*eturn6+ft(2)*wel_loc*eel_loc &
362 ! +edihcnstr+ft(2)*wtor_d*etors_d+ft(1)*wsccor*esccor &
365 ! etot=wsc*(evdw+ft(6)*evdw_t)+wscp*evdw2 &
366 ! +ft(1)*welec*(ees+evdw1) &
367 ! +wang*ebe+ft(1)*wtor*etors+wscloc*escloc &
368 ! +wstrain*ehpb+nss*ebr+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5 &
369 ! +ft(5)*wcorr6*ecorr6+ft(3)*wturn4*eello_turn4 &
370 ! +ft(2)*wturn3*eello_turn3 &
371 ! +ft(5)*wturn6*eturn6+ft(2)*wel_loc*eel_loc+edihcnstr &
372 ! +ft(2)*wtor_d*etors_d+ft(1)*wsccor*esccor &
377 etot=wsc*evdw+wscp*evdw2+welec*ees &
379 +wang*ebe+wtor*etors+wscloc*escloc &
380 +wstrain*ehpb+nss*ebr+wcorr*ecorr+wcorr5*ecorr5 &
381 +wcorr6*ecorr6+wturn4*eello_turn4 &
382 +wturn3*eello_turn3 &
383 +wturn6*eello_turn6+wel_loc*eel_loc &
384 +edihcnstr+wtor_d*etors_d+wsccor*esccor &
385 +wbond*estr+wcatprot*ecation_prot+wcatcat*ecationcation
387 etot=wsc*evdw+wscp*evdw2 &
389 +wang*ebe+wtor*etors+wscloc*escloc &
390 +wstrain*ehpb+nss*ebr+wcorr*ecorr+wcorr5*ecorr5 &
391 +wcorr6*ecorr6+wturn4*eello_turn4 &
392 +wturn3*eello_turn3 &
393 +wturn6*eello_turn6+wel_loc*eel_loc+edihcnstr &
394 +wtor_d*etors_d+wsccor*esccor &
395 +wbond*estr+wcatprot*ecation_prot+wcatcat*ecationcation
399 write (iout,*) i,iparm,1.0d0/(beta_h(ib,iparm)*1.987D-3),&
403 if (iparm.eq.1 .and. ib.eq.1) then
404 write (iout,*)"Conformation",i
407 energia(k)=enetb(k,i,iparm)
409 ! call enerprint(energia(0),fT)
410 call enerprint(energia(0))
417 Econstr=Econstr+Kh(j,kk,ib,iparm) &
418 *(ddW-q0(j,kk,ib,iparm))**2
421 -beta_h(ib,iparm)*(etot-potEmin+Econstr)
423 write (iout,'(4i5,4e15.5)') i,kk,ib,iparm,&
424 etot,potEmin,etot-potEmin,v(i,kk,ib,iparm)
430 ! Simple iteration to calculate free energies corresponding to all simulation
434 ! Compute new free-energy values corresponding to the righ-hand side of the
435 ! equation and their derivatives.
436 write (iout,*) "------------------------fi"
446 vf=v(t,l,k,i)+f(l,k,i)
447 if (vf.gt.vmax) vmax=vf
455 aux=f(l,k,i)+v(t,l,k,i)-vmax
456 if (aux.gt.-200.0d0) &
457 denom=denom+snk(l,k,i,islice)*dexp(aux)
461 entfac(t)=-dlog(denom)-vmax
463 write (iout,*) t,"vmax",vmax," denom",denom,"entfac",entfac(t)
468 do ii=1,nR(iib,iparm)
470 fi_p(ii,iib,iparm)=0.0d0
472 fi_p(ii,iib,iparm)=fi_p(ii,iib,iparm) &
473 +dexp(v(t,ii,iib,iparm)+entfac(t))
475 write (iout,'(4i5,3e15.5)') t,ii,iib,iparm,&
476 v(t,ii,iib,iparm),entfac(t),fi_p(ii,iib,iparm)
480 fi(ii,iib,iparm)=0.0d0
482 fi(ii,iib,iparm)=fi(ii,iib,iparm) &
483 +dexp(v(t,ii,iib,iparm)+entfac(t))
492 write (iout,*) "fi before MPI_Reduce me",me,' master',master
494 do ib=1,nT_h(nparmset)
495 write (iout,*) "iparm",iparm," ib",ib
496 write (iout,*) "beta=",beta_h(ib,iparm)
497 write (iout,'(8e15.5)') (fi_p(i,ib,iparm),i=1,nR(ib,iparm))
501 write (iout,*) "REDUCE size",maxR,MaxT_h,nParmSet,&
503 write (iout,*) "MPI_COMM_WORLD",MPI_COMM_WORLD,&
504 " WHAM_COMM",WHAM_COMM
505 call MPI_Reduce(fi_p(1,1,1),fi(1,1,1),maxR*MaxT_h*nParmSet,&
506 MPI_DOUBLE_PRECISION,&
507 MPI_SUM,Master,WHAM_COMM,IERROR)
509 write (iout,*) "fi after MPI_Reduce nparmset",nparmset
511 write (iout,*) "iparm",iparm
513 write (iout,*) "beta=",beta_h(ib,iparm)
514 write (iout,'(8e15.5)') (fi(i,ib,iparm),i=1,nR(ib,iparm))
518 if (me1.eq.Master) then
524 fi(i,ib,iparm)=-dlog(fi(i,ib,iparm))
525 avefi=avefi+fi(i,ib,iparm)
531 write (iout,*) "Parameter set",iparm
533 write (iout,*) "beta=",beta_h(ib,iparm)
535 fi(i,ib,iparm)=fi(i,ib,iparm)-avefi
537 write (iout,'(8f10.5)') (fi(i,ib,iparm),i=1,nR(ib,iparm))
538 write (iout,'(8f10.5)') (f(i,ib,iparm),i=1,nR(ib,iparm))
542 ! Compute the norm of free-energy increments.
547 finorm=finorm+dabs(fi(i,ib,iparm)-f(i,ib,iparm))
548 f(i,ib,iparm)=fi(i,ib,iparm)
553 write (iout,*) 'Iteration',iter,' finorm',finorm
557 call MPI_Bcast(f(1,1,1),MaxR*MaxT_h*nParmSet,&
558 MPI_DOUBLE_PRECISION,Master,&
560 call MPI_Bcast(finorm,1,MPI_DOUBLE_PRECISION,Master,&
563 ! Exit, if the increment norm is smaller than pre-assigned tolerance.
564 if (finorm.lt.fimin) then
565 write (iout,*) 'Iteration converged'
572 ! Now, put together the histograms from all simulations, in order to get the
573 ! unbiased total histogram.
583 write (iout,*) "--------------hist"
587 sumW_p(i,iparm)=0.0d0
588 sumE_p(i,iparm)=0.0d0
589 sumEbis_p(i,iparm)=0.0d0
590 sumEsq_p(i,iparm)=0.0d0
592 sumQ_p(j,i,iparm)=0.0d0
593 sumQsq_p(j,i,iparm)=0.0d0
594 sumEQ_p(j,i,iparm)=0.0d0
604 sumEbis(i,iparm)=0.0d0
605 sumEsq(i,iparm)=0.0d0
607 sumQ(j,i,iparm)=0.0d0
608 sumQsq(j,i,iparm)=0.0d0
609 sumEQ(j,i,iparm)=0.0d0
615 ! 8/26/05 entropy distribution
620 !! ent=-dlog(entfac(t))
622 ! if (ent.lt.entmin_p) entmin_p=ent
623 ! if (ent.gt.entmax_p) entmax_p=ent
625 ! write (iout,*) "entmin",entmin_p," entmax",entmax_p
626 !! write (iout,*) "entmin_p",entmin_p," entmax_p",entmax_p
628 ! call MPI_Allreduce(entmin_p,entmin,1,MPI_DOUBLE_PRECISION,MPI_MIN,&
630 ! call MPI_Allreduce(entmax_p,entmax,1,MPI_DOUBLE_PRECISION,MPI_MAX,&
632 ! write (iout,*) "entmin",entmin," entmax",entmax
633 ! write (iout,*) "entmin_p",entmin_p," entmax_p",entmax_p
634 ! ientmax=entmax-entmin
635 !iientmax=entmax-entmin !el
636 !write (iout,*) "ientmax",ientmax,entmax,entmin
637 !write (iout,*) "iientmax",iientmax
638 ! if (ientmax.gt.2000) ientmax=2000
639 ! write (iout,*) "entmin",entmin," entmax",entmax," ientmax",ientmax
642 !! ient=-dlog(entfac(t))-entmin
643 ! ient=entfac(t)-entmin
644 ! if (ient.le.2000) histent_p(ient)=histent_p(ient)+1
646 ! call MPI_Allreduce(histent_p(0),histent(0),ientmax+1,MPI_INTEGER,&
647 ! MPI_SUM,WHAM_COMM,IERROR)
648 ! if (me1.eq.Master) then
649 ! write (iout,*) "Entropy histogram"
651 ! write(iout,'(f15.4,i10)') entmin+i,histent(i)
657 ! do t=1,ntot(islice)
659 ! if (ent.lt.entmin) entmin=ent
660 ! if (ent.gt.entmax) entmax=ent
662 ! ientmax=-dlog(entmax)-entmin
663 ! if (ientmax.gt.2000) ientmax=2000
664 ! do t=1,ntot(islice)
665 ! ient=entfac(t)-entmin
666 ! if (ient.le.2000) histent(ient)=histent(ient)+1
668 ! write (iout,*) "Entropy histogram"
670 ! write(iout,'(2f15.4)') entmin+i,histent(i)
675 write (iout,*) "me1",me1," scount",scount(me1) !d
701 hrmsrgy(j,i,ib)=0.0d0
703 hrmsrgy_p(j,i,ib)=0.0d0
715 hfin_ent_p(ind)=hfin_ent_p(ind)+dexp(entfac(t))
717 hfin_ent(ind)=hfin_ent(ind)+dexp(entfac(t))
719 ! write (iout,'(2i5,20f8.2)') "debug",t,t,(enetb(k,t,iparm),k=1,21)
720 call restore_parm(iparm)
721 ! evdw=enetb(21,t,iparm)
722 evdw=enetb(20,t,iparm)
723 evdw_t=enetb(1,t,iparm)
725 ! evdw2_14=enetb(17,t,iparm)
726 evdw2_14=enetb(18,t,iparm)
727 evdw2=enetb(2,t,iparm)+evdw2_14
729 evdw2=enetb(2,t,iparm)
734 evdw1=enetb(16,t,iparm)
739 ecorr=enetb(4,t,iparm)
740 ecorr5=enetb(5,t,iparm)
741 ecorr6=enetb(6,t,iparm)
742 eel_loc=enetb(7,t,iparm)
743 eello_turn3=enetb(8,t,iparm)
744 eello_turn4=enetb(9,t,iparm)
745 eello_turn6=enetb(10,t,iparm)
746 ebe=enetb(11,t,iparm)
747 escloc=enetb(12,t,iparm)
748 etors=enetb(13,t,iparm)
749 etors_d=enetb(14,t,iparm)
750 ehpb=enetb(15,t,iparm)
751 ! estr=enetb(18,t,iparm)
752 estr=enetb(17,t,iparm)
753 ! esccor=enetb(19,t,iparm)
754 esccor=enetb(21,t,iparm)
755 ! edihcnstr=enetb(20,t,iparm)
756 edihcnstr=enetb(19,t,iparm)
758 ecationcation=enetb(42,t,iparm)
759 ecation_prot=enetb(41,t,iparm)
762 betaT=startGridT+k*delta_T
764 !write(iout,*)"kkkkkkkk",betaT,startGridT,k,delta_T
766 !d ft=2*T0/(T0+betaT)
767 if (rescale_modeW.eq.1) then
775 denom=kfacl-1.0d0+quotl
777 ftprim(l)=-l*ft(l)*quotl1/(T0*denom)
778 ftbis(l)=l*kfacl*quotl1* &
779 (2*l*quotl-(l-1)*denom)/(quot*t0*t0*denom**3)
782 ft(6)=(320.0d0+80.0d0*dtanh((betaT-320.0d0)/80.0d0))/ &
784 ftprim(6)=1.0d0/(320.0d0*dcosh((betaT-320.0d0)/80.0d0)**2)
785 ftbis(6)=-2.0d0*dtanh((betaT-320.0d0)/80.0d0) &
786 /(320.0d0*80.0d0*dcosh((betaT-320.0d0)/80.0d0)**3)
796 else if (rescale_modeW.eq.2) then
804 logfac=1.0d0/dlog(eplus+eminus)
805 tanhT=(eplus-eminus)/(eplus+eminus)
806 fT(l)=1.12692801104297249644d0*logfac
807 ftprim(l)=-l*quotl1*ft(l)*tanhT*logfac/T0
808 ftbis(l)=(l-1)*ftprim(l)/(quot*T0)- &
809 2*l*quotl1/T0*logfac* &
810 (2*l*quotl1*ft(l)/(T0*(eplus+eminus)**2) &
814 ft(6)=(320.0d0+80.0d0*dtanh((betaT-320.0d0)/80.0d0))/ &
816 ftprim(6)=1.0d0/(320.0d0*dcosh((betaT-320.0d0)/80.0d0)**2)
817 ftbis(6)=-2.0d0*dtanh((betaT-320.0d0)/80.0d0) &
818 /(320.0d0*80.0d0*dcosh((betaT-320.0d0)/80.0d0)**3)
828 else if (rescale_modeW.eq.0) then
834 write (iout,*) "Error in WHAM_CALC: wrong RESCALE_MODE",&
839 ! write (iout,*) "ftprim",ftprim
840 ! write (iout,*) "ftbis",ftbis
841 betaT=1.0d0/(1.987D-3*betaT)
843 etot=wsc*(evdw+ft(6)*evdw_t)+wscp*evdw2+ft(1)*welec*ees &
845 +wang*ebe+ft(1)*wtor*etors+wscloc*escloc &
846 +wstrain*ehpb+nss*ebr+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5 &
847 +ft(5)*wcorr6*ecorr6+ft(3)*wturn4*eello_turn4 &
848 +ft(2)*wturn3*eello_turn3 &
849 +ft(5)*wturn6*eello_turn6+ft(2)*wel_loc*eel_loc &
850 +edihcnstr+ft(2)*wtor_d*etors_d+ft(1)*wsccor*esccor &
851 +wbond*estr+wcatprot*ecation_prot+wcatcat*ecationcation
852 eprim=ftprim(6)*evdw_t+ftprim(1)*welec*ees &
853 +ftprim(1)*wtor*etors+ &
854 ftprim(3)*wcorr*ecorr+ftprim(4)*wcorr5*ecorr5+ &
855 ftprim(5)*wcorr6*ecorr6+ftprim(3)*wturn4*eello_turn4+ &
856 ftprim(2)*wturn3*eello_turn3+ftprim(5)*wturn6*eello_turn6+ &
857 ftprim(2)*wel_loc*eel_loc+ftprim(2)*wtor_d*etors_d+ &
858 ftprim(1)*wsccor*esccor
859 ebis=ftbis(1)*welec*ees+ftbis(1)*wtor*etors+ &
860 ftbis(3)*wcorr*ecorr+ftbis(4)*wcorr5*ecorr5+ &
861 ftbis(5)*wcorr6*ecorr6+ftbis(3)*wturn4*eello_turn4+ &
862 ftbis(2)*wturn3*eello_turn3+ftbis(5)*wturn6*eello_turn6+ &
863 ftbis(2)*wel_loc*eel_loc+ftbis(2)*wtor_d*etors_d+ &
864 ftbis(1)*wsccor*esccor
866 etot=wsc*(evdw+ft(6)*evdw_t)+wscp*evdw2 &
867 +ft(1)*welec*(ees+evdw1) &
868 +wang*ebe+ft(1)*wtor*etors+wscloc*escloc &
869 +wstrain*ehpb+nss*ebr+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5 &
870 +ft(5)*wcorr6*ecorr6+ft(3)*wturn4*eello_turn4 &
871 +ft(2)*wturn3*eello_turn3 &
872 +ft(5)*wturn6*eello_turn6+ft(2)*wel_loc*eel_loc+edihcnstr &
873 +ft(2)*wtor_d*etors_d+ft(1)*wsccor*esccor &
874 +wbond*estr+wcatprot*ecation_prot+wcatcat*ecationcation
875 eprim=ftprim(6)*evdw_t+ftprim(1)*welec*(ees+evdw1) &
876 +ftprim(1)*wtor*etors+ &
877 ftprim(3)*wcorr*ecorr+ftprim(4)*wcorr5*ecorr5+ &
878 ftprim(5)*wcorr6*ecorr6+ftprim(3)*wturn4*eello_turn4+ &
879 ftprim(2)*wturn3*eello_turn3+ftprim(5)*wturn6*eello_turn6+ &
880 ftprim(2)*wel_loc*eel_loc+ftprim(2)*wtor_d*etors_d+ &
881 ftprim(1)*wsccor*esccor
882 ebis=ftbis(1)*welec*(ees+evdw1)+ftbis(1)*wtor*etors+ &
883 ftbis(3)*wcorr*ecorr+ftbis(4)*wcorr5*ecorr5+ &
884 ftbis(5)*wcorr6*ecorr6+ftbis(3)*wturn4*eello_turn4+ &
885 ftbis(2)*wturn3*eello_turn3+ftbis(5)*wturn6*eello_turn6+ &
886 ftbis(2)*wel_loc*eel_loc+ftbis(2)*wtor_d*etors_d+ &
887 ftprim(1)*wsccor*esccor
889 weight=dexp(-betaT*(etot-potEmin)+entfac(t))
891 write (iout,*) "iparm",iparm," t",t," betaT",betaT,&
892 " etot",etot," entfac",entfac(t),&
893 " weight",weight," ebis",ebis
895 etot=etot-temper*eprim
897 sumW_p(k,iparm)=sumW_p(k,iparm)+weight
898 sumE_p(k,iparm)=sumE_p(k,iparm)+etot*weight
899 sumEbis_p(k,iparm)=sumEbis_p(k,iparm)+ebis*weight
900 sumEsq_p(k,iparm)=sumEsq_p(k,iparm)+etot**2*weight
902 sumQ_p(j,k,iparm)=sumQ_p(j,k,iparm)+q(j,t)*weight
903 sumQsq_p(j,k,iparm)=sumQsq_p(j,k,iparm)+q(j,t)**2*weight
904 sumEQ_p(j,k,iparm)=sumEQ_p(j,k,iparm) &
908 sumW(k,iparm)=sumW(k,iparm)+weight
909 sumE(k,iparm)=sumE(k,iparm)+etot*weight
910 sumEbis(k,iparm)=sumEbis(k,iparm)+ebis*weight
911 sumEsq(k,iparm)=sumEsq(k,iparm)+etot**2*weight
913 sumQ(j,k,iparm)=sumQ(j,k,iparm)+q(j,t)*weight
914 sumQsq(j,k,iparm)=sumQsq(j,k,iparm)+q(j,t)**2*weight
915 sumEQ(j,k,iparm)=sumEQ(j,k,iparm) &
920 indE = aint(potE(t,iparm)-aint(potEmin))
921 if (indE.ge.0 .and. indE.le.maxinde) then
922 if (indE.gt.upindE_p) upindE_p=indE
923 histE_p(indE)=histE_p(indE)+dexp(-entfac(t))
927 expfac=dexp(-beta_h(ib,iparm)*(etot-potEmin)+entfac(t))
928 hfin_p(ind,ib)=hfin_p(ind,ib)+ &
929 dexp(-beta_h(ib,iparm)*(etot-potEmin)+entfac(t))
931 indrgy=dint((q(nQ+2,t)-rgymin)/deltrgy)
932 indrms=dint((q(nQ+1,t)-rmsmin)/deltrms)
933 hrmsrgy_p(indrgy,indrms,ib)= &
934 hrmsrgy_p(indrgy,indrms,ib)+expfac
939 expfac=dexp(-beta_h(ib,iparm)*(etot-potEmin)+entfac(t))
940 hfin(ind,ib)=hfin(ind,ib)+ &
941 dexp(-beta_h(ib,iparm)*(etot-potEmin)+entfac(t))
943 indrgy=dint((q(nQ+2,t)-rgymin)/deltrgy)
944 indrms=dint((q(nQ+1,t)-rmsmin)/deltrms)
945 hrmsrgy(indrgy,indrms,ib)= &
946 hrmsrgy(indrgy,indrms,ib)+expfac
952 if (histout) call MPI_Reduce(hfin_p(0,ib),hfin(0,ib),nbin,&
953 MPI_DOUBLE_PRECISION,MPI_SUM,Master,WHAM_COMM,IERROR)
955 call MPI_Reduce(hrmsrgy_p(0,0,ib),hrmsrgy(0,0,ib),&
956 (MaxBinRgy+1)*(nbin_rms+1),MPI_DOUBLE_PRECISION,MPI_SUM,Master,&
960 call MPI_Reduce(upindE_p,upindE,1,&
961 MPI_INTEGER,MPI_MAX,Master,WHAM_COMM,IERROR)
962 call MPI_Reduce(histE_p(0),histE(0),maxindE,&
963 MPI_DOUBLE_PRECISION,MPI_SUM,Master,WHAM_COMM,IERROR)
965 if (me1.eq.master) then
969 write (iout,'(6x,$)')
970 write (iout,'(f20.2,$)') (1.0d0/(1.987D-3*beta_h(ib,iparm)),&
974 write (iout,'(/a)') 'Final histograms'
976 if (nslice.eq.1) then
977 if (separate_parset) then
978 write(licz3,"(bz,i3.3)") myparm
979 histname=prefix(:ilen(prefix))//'_par'//licz3//'.hist'
981 histname=prefix(:ilen(prefix))//'.hist'
984 if (separate_parset) then
985 write(licz3,"(bz,i3.3)") myparm
986 histname=prefix(:ilen(prefix))//'_par'//licz3// &
987 '_slice_'//licz2//'.hist'
989 histname=prefix(:ilen(prefix))//'_slice_'//licz2//'.hist'
992 #if defined(AIX) || defined(PGI)
993 open (ihist,file=histname,position='append')
995 open (ihist,file=histname,access='append')
1003 sumH=sumH+hfin(t,ib)
1005 if (sumH.gt.0.0d0) then
1007 jj = mod(liczbaW,nbin1)
1008 liczbaW=liczbaW/nbin1
1009 write (iout,'(f6.3,$)') dmin+(jj+0.5d0)*delta
1011 write (ihist,'(f6.3,$)') dmin+(jj+0.5d0)*delta
1014 write (iout,'(e20.10,$)') hfin(t,ib)
1015 if (histfile) write (ihist,'(e20.10,$)') hfin(t,ib)
1017 write (iout,'(i5)') iparm
1018 if (histfile) write (ihist,'(i5)') iparm
1025 if (nslice.eq.1) then
1026 if (separate_parset) then
1027 write(licz3,"(bz,i3.3)") myparm
1028 histname=prefix(:ilen(prefix))//"_par"//licz3//'.ent'
1030 histname=prefix(:ilen(prefix))//'.ent'
1033 if (separate_parset) then
1034 write(licz3,"(bz,i3.3)") myparm
1035 histname=prefix(:ilen(prefix))//'par_'//licz3// &
1036 '_slice_'//licz2//'.ent'
1038 histname=prefix(:ilen(prefix))//'_slice_'//licz2//'.ent'
1041 #if defined(AIX) || defined(PGI)
1042 open (ihist,file=histname,position='append')
1044 open (ihist,file=histname,access='append')
1046 write (ihist,'(a)') "# Microcanonical entropy"
1048 write (ihist,'(f8.0,$)') dint(potEmin)+i
1049 if (histE(i).gt.0.0e0) then
1050 write (ihist,'(f15.5,$)') dlog(histE(i))
1052 write (ihist,'(f15.5,$)') 0.0d0
1058 write (iout,*) "Microcanonical entropy"
1060 write (iout,'(f8.0,$)') dint(potEmin)+i
1061 if (histE(i).gt.0.0e0) then
1062 write (iout,'(f15.5,$)') dlog(histE(i))
1064 write (iout,'(f15.5,$)') 0.0d0
1069 if (nslice.eq.1) then
1070 if (separate_parset) then
1071 write(licz3,"(bz,i3.3)") myparm
1072 histname=prefix(:ilen(prefix))//'_par'//licz3//'.rmsrgy'
1074 histname=prefix(:ilen(prefix))//'.rmsrgy'
1077 if (separate_parset) then
1078 write(licz3,"(bz,i3.3)") myparm
1079 histname=prefix(:ilen(prefix))//'_par'//licz3// &
1080 '_slice_'//licz2//'.rmsrgy'
1082 histname=prefix(:ilen(prefix))//'_slice_'//licz2//'.rmsrgy'
1085 #if defined(AIX) || defined(PGI)
1086 open (ihist,file=histname,position='append')
1088 open (ihist,file=histname,access='append')
1092 write(ihist,'(2f8.2,$)') &
1093 rgymin+deltrgy*j,rmsmin+deltrms*i
1095 if (hrmsrgy(j,i,ib).gt.0.0d0) then
1096 write(ihist,'(e14.5,$)') &
1097 -dlog(hrmsrgy(j,i,ib))/beta_h(ib,iparm) &
1100 write(ihist,'(e14.5,$)') 1.0d6
1103 write (ihist,'(i2)') iparm
1111 call MPI_Reduce(hfin_ent_p(0),hfin_ent(0),nbin,&
1112 MPI_DOUBLE_PRECISION,MPI_SUM,Master,WHAM_COMM,IERROR)
1113 call MPI_Reduce(sumW_p(0,1),sumW(0,1),(nGridT+1)*nParmSet,&
1114 MPI_DOUBLE_PRECISION,MPI_SUM,Master,WHAM_COMM,IERROR)
1115 call MPI_Reduce(sumE_p(0,1),sumE(0,1),(nGridT+1)*nParmSet,&
1116 MPI_DOUBLE_PRECISION,MPI_SUM,Master,WHAM_COMM,IERROR)
1117 call MPI_Reduce(sumEbis_p(0,1),sumEbis(0,1),(nGridT+1)*nParmSet,&
1118 MPI_DOUBLE_PRECISION,MPI_SUM,Master,WHAM_COMM,IERROR)
1119 call MPI_Reduce(sumEsq_p(0,1),sumEsq(0,1),(nGridT+1)*nParmSet,&
1120 MPI_DOUBLE_PRECISION,MPI_SUM,Master,WHAM_COMM,IERROR)
1121 call MPI_Reduce(sumQ_p(1,0,1),sumQ(1,0,1),&
1122 MaxQ1*(nGridT+1)*nParmSet,MPI_DOUBLE_PRECISION,MPI_SUM,Master,&
1124 call MPI_Reduce(sumQsq_p(1,0,1),sumQsq(1,0,1),&
1125 MaxQ1*(nGridT+1)*nParmSet,MPI_DOUBLE_PRECISION,MPI_SUM,Master,&
1127 call MPI_Reduce(sumEQ_p(1,0,1),sumEQ(1,0,1),&
1128 MaxQ1*(nGridT+1)*nParmSet,MPI_DOUBLE_PRECISION,MPI_SUM,Master,&
1130 if (me.eq.master) then
1132 write (iout,'(/a)') 'Thermal characteristics of folding'
1133 if (nslice.eq.1) then
1136 nazwa=prefix(:ilen(prefix))//"_slice_"//licz2
1139 if (nparmset.eq.1 .and. .not.separate_parset) then
1140 nazwa=nazwa(:iln)//".thermal"
1141 else if (nparmset.eq.1 .and. separate_parset) then
1142 write(licz3,"(bz,i3.3)") myparm
1143 nazwa=nazwa(:iln)//"_par_"//licz3//".thermal"
1146 if (nparmset.gt.1) then
1147 write(licz3,"(bz,i3.3)") iparm
1148 nazwa=nazwa(:iln)//"_par_"//licz3//".thermal"
1151 if (separate_parset) then
1152 write (iout,'(a,i3)') "Parameter set",myparm
1154 write (iout,'(a,i3)') "Parameter set",iparm
1157 sumE(i,iparm)=sumE(i,iparm)/sumW(i,iparm)
1158 sumEbis(i,iparm)=(startGridT+i*delta_T)*sumEbis(i,iparm)/ &
1160 sumEsq(i,iparm)=(sumEsq(i,iparm)/sumW(i,iparm) &
1161 -sumE(i,iparm)**2)/(1.987D-3*(startGridT+i*delta_T)**2)
1163 sumQ(j,i,iparm)=sumQ(j,i,iparm)/sumW(i,iparm)
1164 sumQsq(j,i,iparm)=sumQsq(j,i,iparm)/sumW(i,iparm) &
1166 sumEQ(j,i,iparm)=sumEQ(j,i,iparm)/sumW(i,iparm) &
1167 -sumQ(j,i,iparm)*sumE(i,iparm)
1169 sumW(i,iparm)=-dlog(sumW(i,iparm))*(1.987D-3* &
1170 (startGridT+i*delta_T))+potEmin
1171 write (iout,'(f7.1,2f15.5,$)') startGridT+i*delta_T,&
1172 sumW(i,iparm),sumE(i,iparm)
1173 write (iout,'(f10.5,$)') (sumQ(j,i,iparm),j=1,nQ+2)
1174 write (iout,'(e15.5,$)') sumEsq(i,iparm)-sumEbis(i,iparm),&
1175 (sumQsq(j,i,iparm),j=1,nQ+2),(sumEQ(j,i,iparm),j=1,nQ+2)
1177 write (34,'(f7.1,2f15.5,$)') startGridT+i*delta_T,&
1178 sumW(i,iparm),sumE(i,iparm)
1179 write (34,'(f10.5,$)') (sumQ(j,i,iparm),j=1,nQ+2)
1180 write (34,'(e15.5,$)') sumEsq(i,iparm)-sumEbis(i,iparm),&
1181 (sumQsq(j,i,iparm),j=1,nQ+2),(sumEQ(j,i,iparm),j=1,nQ+2)
1189 if (hfin_ent(t).gt.0.0d0) then
1191 jj = mod(liczbaW,nbin1)
1192 write (iout,'(f6.3,e20.10," ent")') dmin+(jj+0.5d0)*delta,&
1194 if (histfile) write (ihist,'(f6.3,e20.10," ent")') &
1195 dmin+(jj+0.5d0)*delta,&
1199 if (histfile) close(ihist)
1203 ! Write data for zscore
1204 if (nslice.eq.1) then
1205 zscname=prefix(:ilen(prefix))//".zsc"
1207 zscname=prefix(:ilen(prefix))//"_slice_"//licz2//".zsc"
1209 #if defined(AIX) || defined(PGI)
1210 open (izsc,file=prefix(:ilen(prefix))//'.zsc',position='append')
1212 open (izsc,file=prefix(:ilen(prefix))//'.zsc',access='append')
1214 write (izsc,'("NQ=",i1," NPARM=",i1)') nQ,nParmSet
1216 write (izsc,'("NT=",i1)') nT_h(iparm)
1218 write (izsc,'("TEMP=",f6.1," NR=",i2," SNK=",$)') &
1219 1.0d0/(beta_h(ib,iparm)*1.987D-3),nR(ib,iparm)
1220 jj = min0(nR(ib,iparm),7)
1221 write (izsc,'(i8,$)') (snk(i,ib,iparm,islice),i=1,jj)
1222 write (izsc,'(a1,$)') (" ",i=22+8*jj+1,79)
1223 write (izsc,'("&")')
1224 if (nR(ib,iparm).gt.7) then
1225 do ii=8,nR(ib,iparm),9
1226 jj = min0(nR(ib,iparm),ii+8)
1227 write (izsc,'(i8,$)') (snk(i,ib,iparm,islice),i=ii,jj)
1228 write (izsc,'(a1,$') (" ",i=(jj-ii+1)*8+1,79)
1229 write (izsc,'("&")')
1232 write (izsc,'("FI=",$)')
1233 jj=min0(nR(ib,iparm),7)
1234 write (izsc,'(f10.5,$)') (fi(i,ib,iparm),i=1,jj)
1235 write (izsc,'(a1,$)') (" ",i=3+10*jj+1,79)
1236 write (izsc,'("&")')
1237 if (nR(ib,iparm).gt.7) then
1238 do ii=8,nR(ib,iparm),9
1239 jj = min0(nR(ib,iparm),ii+8)
1240 write (izsc,'(f10.5,$)') (fi(i,ib,iparm),i=ii,jj)
1241 if (jj.eq.nR(ib,iparm)) then
1244 write (izsc,'(a1,$)') (" ",i=10*(jj-ii+1)+1,79)
1245 write (izsc,'(t80,"&")')
1250 write (izsc,'("KH=",$)')
1251 write (izsc,'(f7.2,$)') (Kh(j,i,ib,iparm),j=1,nQ)
1252 write (izsc,'(" Q0=",$)')
1253 write (izsc,'(f7.5,$)') (q0(j,i,ib,iparm),j=1,nQ)
1264 end subroutine WHAMCALC
1265 !-----------------------------------------------------------------------------
1266 end module wham_calc