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), &
75 fi_p_min(MaxR,MaxT_h,nParmSet) !(MaxR,MaxT_h,Max_Parm)
76 real(kind=8),dimension(0:nGridT,nParmSet) :: sumW_p,sumE_p,&
77 sumEbis_p,sumEsq_p !(0:nGridT,Max_Parm)
78 real(kind=8),dimension(MaxQ1,0:nGridT,nParmSet) :: sumQ_p,&
79 sumQsq_p,sumEQ_p,sumEprim_p !(MaxQ1,0:nGridT,Max_Parm)
80 real(kind=8) :: hfin_p(0:MaxHdim,maxT_h),&
81 hfin_ent_p(0:MaxHdim),histE_p(0:maxindE),sumH,&
82 hrmsrgy_p(0:MaxBinRgy,0:MaxBinRms,maxT_h)
83 real(kind=8) :: rgymin_t,rmsmin_t,rgymax_t,rmsmax_t
84 real(kind=8) :: potEmin_t,potEmin_t_all(maxT_h,Max_Parm)!,entmin_p,entmax_p
85 ! integer :: histent_p(0:2000)
86 logical :: lprint=.true.
88 real(kind=8) :: delta_T=1.0d0,iientmax
89 real(kind=8) :: rgymin,rmsmin,rgymax,rmsmax
90 real(kind=8),dimension(0:nGridT,nParmSet) :: sumW,sumE,&
91 sumEsq,sumEprim,sumEbis !(0:NGridT,Max_Parm)
92 real(kind=8),dimension(MaxQ1,0:nGridT,nParmSet) :: sumQ,&
93 sumQsq,sumEQ !(MaxQ1,0:NGridT,Max_Parm)
94 real(kind=8) :: betaT,weight,econstr
95 real(kind=8) :: fi(MaxR,MaxT_h,nParmSet),& !(MaxR,maxT_h,Max_Parm)
96 ddW,dd1,dd2,hh,dmin,denom,finorm,avefi,pom,&
97 hfin(0:MaxHdim,maxT_h),histE(0:maxindE),&
98 hrmsrgy(0:MaxBinRgy,0:MaxBinRms,maxT_h),&
100 hfin_ent(0:MaxHdim),vmax,aux,fi_min(MaxR,maxT_h,nParmSet), &
101 potEmin_all(maxT_h,Max_Parm),potEmin_min,entfac_min
102 real(kind=8) :: fT(6),fTprim(6),fTbis(6),quot,quotl1,quotl,kfacl,&
103 eprim,ebis,temper,kfac=2.4d0,T0=300.0d0,startGridT=200.0d0,&
104 eplus,eminus,logfac,tanhT,tt
105 real(kind=8) :: etot,evdw,evdw_t,evdw2,ees,evdw1,ebe,etors,&
106 escloc,ehpb,ecorr,ecorr5,ecorr6,eello_turn4,eello_turn3,&
107 eello_turn6,eel_loc,edihcnstr,etors_d,estr,evdw2_14,esccor, &
108 ecationcation,ecation_prot, evdwpp,eespp ,evdwpsb,eelpsb, &
109 evdwsb, eelsb, estr_nucl,ebe_nucl,esbloc,etors_nucl,etors_d_nucl,&
110 ecorr_nucl, ecorr3_nucl,epeppho, escpho, epepbase,escbase,ecation_nucl,&
111 elipbond,elipang,eliplj,elipelec,ecat_prottran,ecation_protang,eliptran
115 integer :: ind_point(maxpoint),upindE,indE
116 character(len=16) :: plik
117 character(len=1) :: licz1
118 character(len=2) :: licz2
119 character(len=3) :: licz3
120 character(len=128) :: nazwa
125 write(licz2,'(bz,i2.2)') islice
127 write (iout,'(//80(1h-)/"Solving WHAM equations for slice",&
128 i2/80(1h-)//)') islice
129 write (iout,*) "delta",delta," nbin1",nbin1
130 write (iout,*) "MaxN",MaxN," MaxQ",MaxQ," MaHdim",MaxHdim
137 potEmin_all(j,i)=1.0d15
154 if (potE(i,j).le.potEmin) potEmin=potE(i,j)
156 if (q(nQ+1,i).lt.rmsmin) rmsmin=q(nQ+1,i)
157 if (q(nQ+1,i).gt.rmsmax) rmsmax=q(nQ+1,i)
158 if (q(nQ+2,i).lt.rgymin) rgymin=q(nQ+2,i)
159 if (q(nQ+2,i).gt.rgymax) rgymax=q(nQ+2,i)
162 ind=(q(j,i)-dmin+1.0d-8)/delta
164 ind_point(i)=ind_point(i)+ind
166 ind_point(i)=ind_point(i)+nbin1**(j-1)*ind
168 ! write (iout,*) "i",i," j",j," q",q(j,i)," ind_point",
171 if (ind_point(i).lt.0 .or. ind_point(i).gt.MaxHdim) then
172 write (iout,*) "Error - index exceeds range for point",i,&
173 " q=",q(j,i)," ind",ind_point(i)
175 write (iout,*) "Processor",me1
177 call MPI_Abort(MPI_COMM_WORLD, Ierror, Errcode )
182 if (ind_point(i).gt.tmax) tmax=ind_point(i)
183 htot(ind_point(i))=htot(ind_point(i))+1
185 write (iout,*) "i",i,"q",(q(j,i),j=1,nQ)," ind",ind_point(i),&
186 " htot",htot(ind_point(i))
193 write (iout,'(a)') "Numbers of counts in Q bins"
195 if (htot(t).gt.0) then
196 write (iout,'(i15,$)') t
199 jj = mod(liczbaW,nbin1)
200 liczbaW=liczbaW/nbin1
201 write (iout,'(i5,$)') jj
203 write (iout,'(i8)') htot(t)
207 write (iout,'(a,i3)') "Number of data points for parameter set",&
209 write (iout,'(i7,$)') ((snk(m,ib,iparm,islice),m=1,nr(ib,iparm)),&
211 write (iout,'(i8)') stot(islice)
217 call MPI_AllReduce(tmax,tmax_t,1,MPI_INTEGER,MPI_MAX,&
220 ! call MPI_AllReduce(potEmin,potEmin_t,1,MPI_DOUBLE_PRECISION,&
221 ! MPI_MIN,WHAM_COMM,IERROR) !????
222 call MPI_AllReduce(rmsmin,rmsmin_t,1,MPI_DOUBLE_PRECISION,&
223 MPI_MIN,WHAM_COMM,IERROR)
224 call MPI_AllReduce(rmsmax,rmsmax_t,1,MPI_DOUBLE_PRECISION,&
225 MPI_MAX,WHAM_COMM,IERROR)
226 call MPI_AllReduce(rgymin,rgymin_t,1,MPI_DOUBLE_PRECISION,&
227 MPI_MIN,WHAM_COMM,IERROR)
228 call MPI_AllReduce(rgymax,rgymax_t,1,MPI_DOUBLE_PRECISION,&
229 MPI_MAX,WHAM_COMM,IERROR)
230 ! potEmin=potEmin_t !/2 try now??
235 write (iout,*) "potEmin",potEmin
237 rmsmin=deltrms*dint(rmsmin/deltrms)
238 rmsmax=deltrms*dint(rmsmax/deltrms)
239 rgymin=deltrms*dint(rgymin/deltrgy)
240 rgymax=deltrms*dint(rgymax/deltrgy)
241 nbin_rms=(rmsmax-rmsmin)/deltrms
242 nbin_rgy=(rgymax-rgymin)/deltrgy
243 write (iout,*) "rmsmin",rmsmin," rmsmax",rmsmax," rgymin",rgymin,&
244 " rgymax",rgymax," nbin_rms",nbin_rms," nbin_rgy",nbin_rgy
251 write (iout,*) "nFi",nFi
252 ! Compute the Boltzmann factor corresponing to restrain potentials in different
259 ! write (9,'(3i5,f10.5)') i,(iparm,potE(i,iparm),iparm=1,nParmSet)
262 write (iout,'(2i5,21f8.2)') i,iparm,&
263 (enetb(k,i,iparm),k=1,21)
265 call restore_parm(iparm)
267 write (iout,*) wsc,wscp,welec,wvdwpp,wang,wtor,wscloc,&
268 wcorr,wcorr5,wcorr6,wturn4,wturn3,wturn6,wel_loc,&
269 wtor_d,wsccor,wbond,wcatcat
272 !el old rascale weights
274 ! if (rescale_modeW.eq.1) then
275 ! quot=1.0d0/(beta_h(ib,iparm)*1.987D-3*T0)
282 ! fT(l)=kfacl/(kfacl-1.0d0+quotl)
285 ! tt = 1.0d0/(beta_h(ib,iparm)*1.987D-3)
286 ! ft(6)=(320.0d0+80.0d0*dtanh((tt-320.0d0)/80.0d0))/320.0d0
287 !#elif defined(FUNCT)
288 ! ft(6)=1.0d0/(beta_h(ib,iparm)*1.987D-3*T0)
292 ! else if (rescale_modeW.eq.2) then
293 ! quot=1.0d0/(T0*beta_h(ib,iparm)*1.987D-3)
297 ! fT(l)=1.12692801104297249644d0/ &
298 ! dlog(dexp(quotl)+dexp(-quotl))
301 ! tt = 1.0d0/(beta_h(ib,iparm)*1.987D-3)
302 ! ft(6)=(320.0d0+80.0d0*dtanh((tt-320.0d0)/80.0d0))/320.0d0
303 !#elif defined(FUNCT)
304 ! ft(6)=1.0d0/(beta_h(ib,iparm)*1.987D-3*T0)
308 ! write (iout,*) 1.0d0/(beta_h(ib,iparm)*1.987D-3),ft
309 ! else if (rescale_modeW.eq.0) then
314 ! write (iout,*) "Error in WHAM_CALC: wrong RESCALE_MODE",&
319 ! el end old rescale weights
320 call rescale_weights(1.0d0/(beta_h(ib,iparm)*1.987D-3))
322 ! call etot(enetb(0,i,iparm))
323 evdw=enetb(1,i,iparm)
324 ! evdw_t=enetb(21,i,iparm)
325 evdw_t=enetb(20,i,iparm)
327 ! evdw2_14=enetb(17,i,iparm)
328 evdw2_14=enetb(18,i,iparm)
329 evdw2=enetb(2,i,iparm)+evdw2_14
331 evdw2=enetb(2,i,iparm)
336 evdw1=enetb(16,i,iparm)
341 ecorr=enetb(4,i,iparm)
342 ecorr5=enetb(5,i,iparm)
343 ecorr6=enetb(6,i,iparm)
344 eel_loc=enetb(7,i,iparm)
345 eello_turn3=enetb(8,i,iparm)
346 eello_turn4=enetb(9,i,iparm)
347 eello_turn6=enetb(10,i,iparm)
348 ebe=enetb(11,i,iparm)
349 escloc=enetb(12,i,iparm)
350 etors=enetb(13,i,iparm)
351 etors_d=enetb(14,i,iparm)
352 ehpb=enetb(15,i,iparm)
353 ! estr=enetb(18,i,iparm)
354 estr=enetb(17,i,iparm)
355 ! esccor=enetb(19,i,iparm)
356 esccor=enetb(21,i,iparm)
357 ! edihcnstr=enetb(20,i,iparm)
358 edihcnstr=enetb(19,i,iparm)
359 ecationcation=enetb(41,i,iparm)
360 ecation_prot=enetb(42,i,iparm)
361 evdwpp = enetb(26,i,iparm)
362 eespp = enetb(27,i,iparm)
363 evdwpsb = enetb(28,i,iparm)
364 eelpsb = enetb(29,i,iparm)
365 evdwsb = enetb(30,i,iparm)
366 eelsb = enetb(31,i,iparm)
367 estr_nucl = enetb(32,i,iparm)
368 ebe_nucl = enetb(33,i,iparm)
369 esbloc = enetb(34,i,iparm)
370 etors_nucl = enetb(35,i,iparm)
371 etors_d_nucl = enetb(36,i,iparm)
372 ecorr_nucl = enetb(37,i,iparm)
373 ecorr3_nucl = enetb(38,i,iparm)
374 epeppho= enetb(49,i,iparm)
375 escpho= enetb(48,i,iparm)
376 epepbase= enetb(47,i,iparm)
377 escbase= enetb(46,i,iparm)
378 ecation_nucl= enetb(50,i,iparm)
379 elipbond=enetb(52,i,iparm)
380 elipang=enetb(53,i,iparm)
381 eliplj=enetb(54,i,iparm)
382 elipelec=enetb(55,i,iparm)
383 ecat_prottran=enetb(56,i,iparm)
384 ecation_protang=enetb(57,i,iparm)
385 eliptran=enetb(22,i,iparm)
387 write (iout,'(3i5,6f5.2,15f12.3)') i,ib,iparm,(ft(l),l=1,6),&
388 evdw+evdw_t,evdw2,ees,evdw1,ecorr,eel_loc,estr,ebe,escloc,&
389 etors,etors_d,eello_turn3,eello_turn4,esccor,ecationcation
393 ! etot=wsc*(evdw+ft(6)*evdw_t)+wscp*evdw2+ft(1)*welec*ees &
395 ! +wang*ebe+ft(1)*wtor*etors+wscloc*escloc &
396 ! +wstrain*ehpb+nss*ebr+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5 &
397 ! +ft(5)*wcorr6*ecorr6+ft(3)*wturn4*eello_turn4 &
398 ! +ft(2)*wturn3*eello_turn3 &
399 ! +ft(5)*wturn6*eturn6+ft(2)*wel_loc*eel_loc &
400 ! +edihcnstr+ft(2)*wtor_d*etors_d+ft(1)*wsccor*esccor &
403 ! etot=wsc*(evdw+ft(6)*evdw_t)+wscp*evdw2 &
404 ! +ft(1)*welec*(ees+evdw1) &
405 ! +wang*ebe+ft(1)*wtor*etors+wscloc*escloc &
406 ! +wstrain*ehpb+nss*ebr+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5 &
407 ! +ft(5)*wcorr6*ecorr6+ft(3)*wturn4*eello_turn4 &
408 ! +ft(2)*wturn3*eello_turn3 &
409 ! +ft(5)*wturn6*eturn6+ft(2)*wel_loc*eel_loc+edihcnstr &
410 ! +ft(2)*wtor_d*etors_d+ft(1)*wsccor*esccor &
415 etot=wsc*evdw+wscp*evdw2+welec*ees &
417 +wang*ebe+wtor*etors+wscloc*escloc &
418 +wstrain*ehpb+wcorr*ecorr+wcorr5*ecorr5 &
419 +wcorr6*ecorr6+wturn4*eello_turn4 &
420 +wturn3*eello_turn3 &
421 +wturn6*eello_turn6+wel_loc*eel_loc &
422 +edihcnstr+wtor_d*etors_d+wsccor*esccor &
423 +wbond*estr+wcatprot*ecation_prot+wcatcat*ecationcation&
424 +wbond_nucl*estr_nucl+wang_nucl*ebe_nucl&
425 +wvdwpp_nucl*evdwpp+welpp*eespp+wvdwpsb*evdwpsb+welpsb*eelpsb&
426 +wvdwsb*evdwsb+welsb*eelsb+wsbloc*esbloc+wtor_nucl*etors_nucl&
427 +wtor_d_nucl*etors_d_nucl+wcorr_nucl*ecorr_nucl+wcorr3_nucl*ecorr3_nucl&
429 +wpepbase*epepbase+wscpho*escpho+wpeppho*epeppho+wcatnucl*ecation_nucl&
430 +elipbond+elipang+eliplj+elipelec+wcat_tran*ecat_prottran+ecation_protang&
435 etot=wsc*evdw+wscp*evdw2 &
437 +wang*ebe+wtor*etors+wscloc*escloc &
438 +wstrain*ehpb+wcorr*ecorr+wcorr5*ecorr5 &
439 +wcorr6*ecorr6+wturn4*eello_turn4 &
440 +wturn3*eello_turn3 &
441 +wturn6*eello_turn6+wel_loc*eel_loc+edihcnstr &
442 +wtor_d*etors_d+wsccor*esccor &
443 +wbond*estr+wcatprot*ecation_prot+wcatcat*ecationcation&
444 +wbond_nucl*estr_nucl+wang_nucl*ebe_nucl&
445 +wvdwpp_nucl*evdwpp+welpp*eespp+wvdwpsb*evdwpsb+welpsb*eelpsb&
446 +wvdwsb*evdwsb+welsb*eelsb+wsbloc*esbloc+wtor_nucl*etors_nucl&
447 +wtor_d_nucl*etors_d_nucl+wcorr_nucl*ecorr_nucl+wcorr3_nucl*ecorr3_nucl&
449 +wpepbase*epepbase+wscpho*escpho+wpeppho*epeppho+wcatnucl*ecation_nucl&
450 +elipbond+elipang+eliplj+elipelec+wcat_tran*ecat_prottran+ecation_protang&
459 write (iout,*) i,iparm,1.0d0/(beta_h(ib,iparm)*1.987D-3),&
463 if (iparm.eq.1 .and. ib.eq.1) then
464 write (iout,*)"Conformation",i
467 energia(k)=enetb(k,i,iparm)
469 ! call enerprint(energia(0),fT)
470 call enerprint(energia(0))
477 Econstr=Econstr+Kh(j,kk,ib,iparm) &
478 *(ddW-q0(j,kk,ib,iparm))**2
481 -beta_h(ib,iparm)*(etot-potEmin+Econstr)
483 write (iout,'(4i5,4e15.5)') i,kk,ib,iparm,&
484 etot,potEmin,etot-potEmin,v(i,kk,ib,iparm)
490 ! Simple iteration to calculate free energies corresponding to all simulation
494 ! Compute new free-energy values corresponding to the righ-hand side of the
495 ! equation and their derivatives.
496 write (iout,*) "------------------------fi"
506 vf=v(t,l,k,i)+f(l,k,i)
507 if (vf.gt.vmax) vmax=vf
515 aux=f(l,k,i)+v(t,l,k,i)-vmax
516 if (aux.gt.-200.0d0) &
517 denom=denom+snk(l,k,i,islice)*dexp(aux)
521 entfac(t)=-dlog(denom)-vmax
522 if (entfac(t).lt.entfac_min) entfac_min=entfac(t)
524 write (iout,*) t,"vmax",vmax," denom",denom,"entfac",entfac(t)
529 do ii=1,nR(iib,iparm)
531 fi_p_min(ii,iib,iparm)=-1.0d5
533 ! fi_p(ii,iib,iparm)=fi_p(ii,iib,iparm) &
534 ! +dexp(v(t,ii,iib,iparm)+entfac(t))
535 aux=v(t,ii,iib,iparm)+entfac(t)
536 if (aux.gt.fi_p_min(ii,iib,iparm)) fi_p_min(ii,iib,iparm)=aux
540 write (iout,'(4i5,4e15.5)') t,ii,iib,iparm,&
541 v(t,ii,iib,iparm),entfac(t),fi_p(ii,iib,iparm),fi_p_min(ii,iib,iparm)
546 ! fi(ii,iib,iparm)=0.0d0
548 ! fi(ii,iib,iparm)=fi(ii,iib,iparm) &
549 ! +dexp(v(t,ii,iib,iparm)+entfac(t))
550 aux=v(t,ii,iib,iparm)+entfac(t)
551 if (aux.gt.fi_min(ii,iib,iparm))
552 & fi_min(ii,iib,iparm)=aux
562 write (iout,*) "fi before MPI_Reduce me",me,' master',master
564 do ib=1,nT_h(nparmset)
565 write (iout,*) "iparm",iparm," ib",ib
566 write (iout,*) "beta=",beta_h(ib,iparm)
567 write (iout,'(8e15.5)') (fi_p(i,ib,iparm),i=1,nR(ib,iparm))
568 write (iout,'(8e15.5)') (fi_p_min(i,ib,iparm),i=1,nR(ib,iparm))
572 call MPI_AllReduce(fi_p_min,fi_min,MaxR*MaxT_h*nParmSet, &
573 MPI_DOUBLE_PRECISION,MPI_MAX,WHAM_COMM,IERROR)
576 write (iout,*) "fi_min after AllReduce"
579 write (iout,*) (i,j,k,fi_min(k,j,i),k=1,nR(j,i))
587 do ii=1,nR(iib,iparm)
589 fi_p(ii,iib,iparm)=0.0d0
591 fi_p(ii,iib,iparm)=fi_p(ii,iib,iparm) &
592 +dexp(v(t,ii,iib,iparm)+entfac(t)-fi_min(ii,iib,iparm))
594 write (iout,'(4i5,4e15.5)') t,ii,iib,iparm, &
595 v(t,ii,iib,iparm),entfac(t),fi_min(ii,iib,iparm), &
600 fi(ii,iib,iparm)=0.0d0
602 fi(ii,iib,iparm)=fi(ii,iib,iparm) &
603 +dexp(v(t,ii,iib,iparm)+entfac(t)-fi_min(ii,iib,iparm))
612 write (iout,*) "fi before MPI_Reduce me",me,' master',master
614 do ib=1,nT_h(nparmset)
615 write (iout,*) "iparm",iparm," ib",ib
616 write (iout,*) "beta=",beta_h(ib,iparm)
617 write (iout,'(8e15.5)') (fi_p(i,ib,iparm),i=1,nR(ib,iparm))
622 write (iout,*) "REDUCE size",maxR,MaxT_h,nParmSet, &
624 write (iout,*) "MPI_COMM_WORLD",MPI_COMM_WORLD, &
625 " WHAM_COMM",WHAM_COMM
628 write (iout,*) "REDUCE size",maxR,MaxT_h,nParmSet,&
630 write (iout,*) "MPI_COMM_WORLD",MPI_COMM_WORLD,&
631 " WHAM_COMM",WHAM_COMM
632 call MPI_Reduce(fi_p(1,1,1),fi(1,1,1),maxR*MaxT_h*nParmSet,&
633 MPI_DOUBLE_PRECISION,&
634 MPI_SUM,Master,WHAM_COMM,IERROR)
636 write (iout,*) "fi after MPI_Reduce nparmset",nparmset
638 write (iout,*) "iparm",iparm
640 write (iout,*) "beta=",beta_h(ib,iparm)
641 write (iout,'(8e15.5)') (fi(i,ib,iparm),i=1,nR(ib,iparm))
645 if (me1.eq.Master) then
651 fi(i,ib,iparm)=-dlog(fi(i,ib,iparm))-fi_min(i,ib,iparm)
652 avefi=avefi+fi(i,ib,iparm)
658 write (iout,*) "Parameter set",iparm
660 write (iout,*) "beta=",beta_h(ib,iparm)
662 fi(i,ib,iparm)=fi(i,ib,iparm)-avefi
664 write (iout,'(8f10.5)') (fi(i,ib,iparm),i=1,nR(ib,iparm))
665 write (iout,'(8f10.5)') (f(i,ib,iparm),i=1,nR(ib,iparm))
669 ! Compute the norm of free-energy increments.
674 finorm=finorm+dabs(fi(i,ib,iparm)-f(i,ib,iparm))
675 f(i,ib,iparm)=fi(i,ib,iparm)
680 write (iout,*) 'Iteration',iter,' finorm',finorm
684 call MPI_Bcast(f(1,1,1),MaxR*MaxT_h*nParmSet,&
685 MPI_DOUBLE_PRECISION,Master,&
687 call MPI_Bcast(finorm,1,MPI_DOUBLE_PRECISION,Master,&
690 ! Exit, if the increment norm is smaller than pre-assigned tolerance.
691 if (finorm.lt.fimin) then
692 write (iout,*) 'Iteration converged'
699 ! Now, put together the histograms from all simulations, in order to get the
700 ! unbiased total histogram.
701 !C Determine the minimum free energies
707 !c write (9,'(3i5,f10.5)') i,(iparm,potE(i,iparm),iparm=1,nParmSet)
710 write (iout,'(2i5,21f8.2)') i,iparm, &
711 (enetb(k,i,iparm),k=1,26)
713 call restore_parm(iparm)
715 write (iout,*) wsc,wscp,welec,wvdwpp,wang,wtor,wscloc, &
716 wcorr,wcorr5,wcorr6,wturn4,wturn3,wturn6,wel_loc, &
720 if (rescale_modeW.eq.1) then
721 quot=1.0d0/(beta_h(ib,iparm)*1.987D-3*T0)
728 fT(l)=kfacl/(kfacl-1.0d0+quotl)
731 tt = 1.0d0/(beta_h(ib,iparm)*1.987D-3)
732 ft(6)=(320.0d0+80.0d0*dtanh((tt-320.0d0)/80.0d0))/320.0d0
734 ft(6)=1.0d0/(beta_h(ib,iparm)*1.987D-3*T0)
738 else if (rescale_modeW.eq.2) then
739 quot=1.0d0/(T0*beta_h(ib,iparm)*1.987D-3)
743 fT(l)=1.12692801104297249644d0/ &
744 dlog(dexp(quotl)+dexp(-quotl))
747 tt = 1.0d0/(beta_h(ib,iparm)*1.987D-3)
748 ft(6)=(320.0d0+80.0d0*dtanh((tt-320.0d0)/80.0d0))/320.0d0
750 ft(6)=1.0d0/(beta_h(ib,iparm)*1.987D-3*T0)
754 !c write (iout,*) 1.0d0/(beta_h(ib,iparm)*1.987D-3),ft
755 else if (rescale_modeW.eq.0) then
760 write (iout,*) "Error in WHAM_CALC: wrong RESCALE_MODE", &
765 evdw=enetb(1,i,iparm)
766 ! evdw_t=enetb(21,i,iparm)
767 evdw_t=enetb(20,i,iparm)
769 ! evdw2_14=enetb(17,i,iparm)
770 evdw2_14=enetb(18,i,iparm)
771 evdw2=enetb(2,i,iparm)+evdw2_14
773 evdw2=enetb(2,i,iparm)
778 evdw1=enetb(16,i,iparm)
783 ecorr=enetb(4,i,iparm)
784 ecorr5=enetb(5,i,iparm)
785 ecorr6=enetb(6,i,iparm)
786 eel_loc=enetb(7,i,iparm)
787 eello_turn3=enetb(8,i,iparm)
788 eello_turn4=enetb(9,i,iparm)
789 eello_turn6=enetb(10,i,iparm)
790 ebe=enetb(11,i,iparm)
791 escloc=enetb(12,i,iparm)
792 etors=enetb(13,i,iparm)
793 etors_d=enetb(14,i,iparm)
794 ehpb=enetb(15,i,iparm)
795 ! estr=enetb(18,i,iparm)
796 estr=enetb(17,i,iparm)
797 ! esccor=enetb(19,i,iparm)
798 esccor=enetb(21,i,iparm)
799 ! edihcnstr=enetb(20,i,iparm)
800 edihcnstr=enetb(19,i,iparm)
801 ! ehomology_constr=enetb(22,i,iparm)
802 ! esaxs=enetb(26,i,iparm)
803 ecationcation=enetb(41,i,iparm)
804 ecation_prot=enetb(42,i,iparm)
805 evdwpp = enetb(26,i,iparm)
806 eespp = enetb(27,i,iparm)
807 evdwpsb = enetb(28,i,iparm)
808 eelpsb = enetb(29,i,iparm)
809 evdwsb = enetb(30,i,iparm)
810 eelsb = enetb(31,i,iparm)
811 estr_nucl = enetb(32,i,iparm)
812 ebe_nucl = enetb(33,i,iparm)
813 esbloc = enetb(34,i,iparm)
814 etors_nucl = enetb(35,i,iparm)
815 etors_d_nucl = enetb(36,i,iparm)
816 ecorr_nucl = enetb(37,i,iparm)
817 ecorr3_nucl = enetb(38,i,iparm)
818 epeppho= enetb(49,i,iparm)
819 escpho= enetb(48,i,iparm)
820 epepbase= enetb(47,i,iparm)
821 escbase= enetb(46,i,iparm)
822 ecation_nucl=enetb(50,i,iparm)
823 elipbond=enetb(52,i,iparm)
824 elipang=enetb(53,i,iparm)
825 eliplj=enetb(54,i,iparm)
826 elipelec=enetb(55,i,iparm)
827 ecat_prottran=enetb(56,i,iparm)
828 ecation_protang=enetb(57,i,iparm)
829 eliptran=enetb(22,i,iparm)
831 etot=wsc*(evdw+ft(6)*evdw_t)+wscp*evdw2+ft(1)*welec*ees &
833 +wang*ebe+ft(1)*wtor*etors+wscloc*escloc &
834 +wstrain*ehpb+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5 &
835 +ft(5)*wcorr6*ecorr6+ft(3)*wturn4*eello_turn4 &
836 +ft(2)*wturn3*eello_turn3 &
837 +ft(5)*wturn6*eello_turn6+ft(2)*wel_loc*eel_loc &
838 +edihcnstr+ft(2)*wtor_d*etors_d+ft(1)*wsccor*esccor &
839 +wbond*estr+wcatprot*ecation_prot+wcatcat*ecationcation &
840 +wbond_nucl*estr_nucl+wang_nucl*ebe_nucl&
841 +wvdwpp_nucl*evdwpp+welpp*eespp+wvdwpsb*evdwpsb+welpsb*eelpsb&
842 +wvdwsb*evdwsb+welsb*eelsb+wsbloc*esbloc+wtor_nucl*etors_nucl&
843 *ft(1)+wtor_d_nucl*ft(2)*etors_d_nucl+wcorr_nucl*ecorr_nucl &
844 +wcorr3_nucl*ecorr3_nucl&
846 +ft(1)*wpepbase*epepbase+wscpho*escpho+wpeppho*epeppho+ wcatnucl*ecation_nucl&
847 +elipbond+elipang+eliplj+elipelec+wcat_tran*ecat_prottran+ecation_protang&
852 etot=wsc*(evdw+ft(6)*evdw_t)+wscp*evdw2 &
853 +ft(1)*welec*(ees+evdw1) &
854 +wang*ebe+ft(1)*wtor*etors+wscloc*escloc &
855 +wstrain*ehpb+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5 &
856 +ft(5)*wcorr6*ecorr6+ft(3)*wturn4*eello_turn4 &
857 +ft(2)*wturn3*eello_turn3 &
858 +ft(5)*wturn6*eello_turn6+ft(2)*wel_loc*eel_loc+edihcnstr &
859 +ft(2)*wtor_d*etors_d+ft(1)*wsccor*esccor &
860 +wbond*estr+wcatprot*ecation_prot+wcatcat*ecationcation&
861 +wbond_nucl*estr_nucl+wang_nucl*ebe_nucl&
862 +wvdwpp_nucl*evdwpp+welpp*eespp+wvdwpsb*evdwpsb+welpsb*eelpsb&
863 +wvdwsb*evdwsb+welsb*eelsb+wsbloc*esbloc+wtor_nucl*etors_nucl&
864 *ft(1)+wtor_d_nucl*ft(2)*etors_d_nucl+wcorr_nucl*ecorr_nucl &
865 +wcorr3_nucl*ecorr3_nucl&
867 +ft(1)*wpepbase*epepbase+wscpho*escpho+wpeppho*epeppho+ wcatnucl*ecation_nucl&
868 +elipbond+elipang+eliplj+elipelec+wcat_tran*ecat_prottran+ecation_protang&
874 write (iout,*) "WTF,",etot,potEmin_all(ib,iparm),entfac(i)/beta_h(ib,iparm)
875 etot=etot-entfac(i)/beta_h(ib,iparm)
876 if(etot.lt.potEmin_all(ib,iparm)) potEmin_all(ib,iparm)=etot
882 write (iout,*) "The potEmin array before reduction"
884 write (iout,*) "Parameter set",i
886 write (iout,*) j,PotEmin_all(j,i)
889 write (iout,*) "potEmin_min",potEmin_min
892 !C Determine the minimum energes for all parameter sets and temperatures
893 call MPI_AllReduce(potEmin_all(1,1),potEmin_t_all(1,1), &
894 maxT_h*nParmSet,MPI_DOUBLE_PRECISION,MPI_MIN,WHAM_COMM,IERROR)
897 potEmin_all(j,i)=potEmin_t_all(j,i)
901 potEmin_min=potEmin_all(1,1)
904 if (potEmin_all(j,i).lt.potEmin_min) &
905 potEmin_min=potEmin_all(j,i)
909 write (iout,*) "The potEmin array"
911 write (iout,*) "Parameter set",i
913 write (iout,*) j,1.0d0/(1.987d-3*beta_h(j,i)), &
917 write (iout,*) "potEmin_min",potEmin_min
929 write (iout,*) "--------------hist"
933 sumW_p(i,iparm)=0.0d0
934 sumE_p(i,iparm)=0.0d0
935 sumEbis_p(i,iparm)=0.0d0
936 sumEsq_p(i,iparm)=0.0d0
938 sumQ_p(j,i,iparm)=0.0d0
939 sumQsq_p(j,i,iparm)=0.0d0
940 sumEQ_p(j,i,iparm)=0.0d0
950 sumEbis(i,iparm)=0.0d0
951 sumEsq(i,iparm)=0.0d0
953 sumQ(j,i,iparm)=0.0d0
954 sumQsq(j,i,iparm)=0.0d0
955 sumEQ(j,i,iparm)=0.0d0
961 ! 8/26/05 entropy distribution
966 !! ent=-dlog(entfac(t))
968 ! if (ent.lt.entmin_p) entmin_p=ent
969 ! if (ent.gt.entmax_p) entmax_p=ent
971 ! write (iout,*) "entmin",entmin_p," entmax",entmax_p
972 !! write (iout,*) "entmin_p",entmin_p," entmax_p",entmax_p
974 ! call MPI_Allreduce(entmin_p,entmin,1,MPI_DOUBLE_PRECISION,MPI_MIN,&
976 ! call MPI_Allreduce(entmax_p,entmax,1,MPI_DOUBLE_PRECISION,MPI_MAX,&
978 ! write (iout,*) "entmin",entmin," entmax",entmax
979 ! write (iout,*) "entmin_p",entmin_p," entmax_p",entmax_p
980 ! ientmax=entmax-entmin
981 !iientmax=entmax-entmin !el
982 !write (iout,*) "ientmax",ientmax,entmax,entmin
983 !write (iout,*) "iientmax",iientmax
984 ! if (ientmax.gt.2000) ientmax=2000
985 ! write (iout,*) "entmin",entmin," entmax",entmax," ientmax",ientmax
988 !! ient=-dlog(entfac(t))-entmin
989 ! ient=entfac(t)-entmin
990 ! if (ient.le.2000) histent_p(ient)=histent_p(ient)+1
992 ! call MPI_Allreduce(histent_p(0),histent(0),ientmax+1,MPI_INTEGER,&
993 ! MPI_SUM,WHAM_COMM,IERROR)
994 ! if (me1.eq.Master) then
995 ! write (iout,*) "Entropy histogram"
997 ! write(iout,'(f15.4,i10)') entmin+i,histent(i)
1003 ! do t=1,ntot(islice)
1005 ! if (ent.lt.entmin) entmin=ent
1006 ! if (ent.gt.entmax) entmax=ent
1008 ! ientmax=-dlog(entmax)-entmin
1009 ! if (ientmax.gt.2000) ientmax=2000
1010 ! do t=1,ntot(islice)
1011 ! ient=entfac(t)-entmin
1012 ! if (ient.le.2000) histent(ient)=histent(ient)+1
1014 ! write (iout,*) "Entropy histogram"
1016 ! write(iout,'(2f15.4)') entmin+i,histent(i)
1021 write (iout,*) "me1",me1," scount",scount(me1) !d
1047 hrmsrgy(j,i,ib)=0.0d0
1049 hrmsrgy_p(j,i,ib)=0.0d0
1061 hfin_ent_p(ind)=hfin_ent_p(ind)+dexp(entfac(t))
1063 hfin_ent(ind)=hfin_ent(ind)+dexp(entfac(t))
1065 ! write (iout,'(2i5,20f8.2)') "debug",t,t,(enetb(k,t,iparm),k=1,21)
1066 call restore_parm(iparm)
1067 ! evdw=enetb(21,t,iparm)
1068 evdw=enetb(20,t,iparm)
1069 evdw_t=enetb(1,t,iparm)
1071 ! evdw2_14=enetb(17,t,iparm)
1072 evdw2_14=enetb(18,t,iparm)
1073 evdw2=enetb(2,t,iparm)+evdw2_14
1075 evdw2=enetb(2,t,iparm)
1079 ees=enetb(3,t,iparm)
1080 evdw1=enetb(16,t,iparm)
1082 ees=enetb(3,t,iparm)
1085 ecorr=enetb(4,t,iparm)
1086 ecorr5=enetb(5,t,iparm)
1087 ecorr6=enetb(6,t,iparm)
1088 eel_loc=enetb(7,t,iparm)
1089 eello_turn3=enetb(8,t,iparm)
1090 eello_turn4=enetb(9,t,iparm)
1091 eello_turn6=enetb(10,t,iparm)
1092 ebe=enetb(11,t,iparm)
1093 escloc=enetb(12,t,iparm)
1094 etors=enetb(13,t,iparm)
1095 etors_d=enetb(14,t,iparm)
1096 ehpb=enetb(15,t,iparm)
1097 ! estr=enetb(18,t,iparm)
1098 estr=enetb(17,t,iparm)
1099 ! esccor=enetb(19,t,iparm)
1100 esccor=enetb(21,t,iparm)
1101 ! edihcnstr=enetb(20,t,iparm)
1102 edihcnstr=enetb(19,t,iparm)
1104 ecationcation=enetb(41,t,iparm)
1105 ecation_prot=enetb(42,t,iparm)
1106 evdwpp = enetb(26,t,iparm)
1107 eespp = enetb(27,t,iparm)
1108 evdwpsb = enetb(28,t,iparm)
1109 eelpsb = enetb(29,t,iparm)
1110 evdwsb = enetb(30,t,iparm)
1111 eelsb = enetb(31,t,iparm)
1112 estr_nucl = enetb(32,t,iparm)
1113 ebe_nucl = enetb(33,t,iparm)
1114 esbloc = enetb(34,t,iparm)
1115 etors_nucl = enetb(35,t,iparm)
1116 etors_d_nucl = enetb(36,t,iparm)
1117 ecorr_nucl = enetb(37,t,iparm)
1118 ecorr3_nucl = enetb(38,t,iparm)
1119 epeppho= enetb(49,t,iparm)
1120 escpho= enetb(48,t,iparm)
1121 epepbase= enetb(47,t,iparm)
1122 escbase= enetb(46,t,iparm)
1123 ecation_nucl=enetb(50,t,iparm)
1124 elipbond=enetb(52,t,iparm)
1125 elipang=enetb(53,t,iparm)
1126 eliplj=enetb(54,t,iparm)
1127 elipelec=enetb(55,t,iparm)
1128 ecat_prottran=enetb(56,t,iparm)
1129 ecation_protang=enetb(57,t,iparm)
1130 eliptran=enetb(22,t,iparm)
1132 betaT=startGridT+k*delta_T
1134 !write(iout,*)"kkkkkkkk",betaT,startGridT,k,delta_T
1136 !d ft=2*T0/(T0+betaT)
1137 if (rescale_modeW.eq.1) then
1145 denom=kfacl-1.0d0+quotl
1147 ftprim(l)=-l*ft(l)*quotl1/(T0*denom)
1148 ftbis(l)=l*kfacl*quotl1* &
1149 (2*l*quotl-(l-1)*denom)/(quot*t0*t0*denom**3)
1152 ft(6)=(320.0d0+80.0d0*dtanh((betaT-320.0d0)/80.0d0))/ &
1154 ftprim(6)=1.0d0/(320.0d0*dcosh((betaT-320.0d0)/80.0d0)**2)
1155 ftbis(6)=-2.0d0*dtanh((betaT-320.0d0)/80.0d0) &
1156 /(320.0d0*80.0d0*dcosh((betaT-320.0d0)/80.0d0)**3)
1157 #elif defined(FUNCT)
1166 else if (rescale_modeW.eq.2) then
1174 logfac=1.0d0/dlog(eplus+eminus)
1175 tanhT=(eplus-eminus)/(eplus+eminus)
1176 fT(l)=1.12692801104297249644d0*logfac
1177 ftprim(l)=-l*quotl1*ft(l)*tanhT*logfac/T0
1178 ftbis(l)=(l-1)*ftprim(l)/(quot*T0)- &
1179 2*l*quotl1/T0*logfac* &
1180 (2*l*quotl1*ft(l)/(T0*(eplus+eminus)**2) &
1184 ft(6)=(320.0d0+80.0d0*dtanh((betaT-320.0d0)/80.0d0))/ &
1186 ftprim(6)=1.0d0/(320.0d0*dcosh((betaT-320.0d0)/80.0d0)**2)
1187 ftbis(6)=-2.0d0*dtanh((betaT-320.0d0)/80.0d0) &
1188 /(320.0d0*80.0d0*dcosh((betaT-320.0d0)/80.0d0)**3)
1189 #elif defined(FUNCT)
1198 else if (rescale_modeW.eq.0) then
1204 write (iout,*) "Error in WHAM_CALC: wrong RESCALE_MODE",&
1209 ! write (iout,*) "ftprim",ftprim
1210 ! write (iout,*) "ftbis",ftbis
1211 betaT=1.0d0/(1.987D-3*betaT)
1213 etot=wsc*(evdw+ft(6)*evdw_t)+wscp*evdw2+ft(1)*welec*ees &
1215 +wang*ebe+ft(1)*wtor*etors+wscloc*escloc &
1216 +wstrain*ehpb+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5 &
1217 +ft(5)*wcorr6*ecorr6+ft(3)*wturn4*eello_turn4 &
1218 +ft(2)*wturn3*eello_turn3 &
1219 +ft(5)*wturn6*eello_turn6+ft(2)*wel_loc*eel_loc &
1220 +edihcnstr+ft(2)*wtor_d*etors_d+ft(1)*wsccor*esccor &
1221 +wbond*estr+wcatprot*ecation_prot+wcatcat*ecationcation &
1222 +wbond_nucl*estr_nucl+wang_nucl*ebe_nucl&
1223 +wvdwpp_nucl*evdwpp+welpp*eespp+wvdwpsb*evdwpsb+welpsb*eelpsb&
1224 +wvdwsb*evdwsb+welsb*eelsb+wsbloc*esbloc+wtor_nucl*etors_nucl&
1225 *ft(1)+wtor_d_nucl*ft(2)*etors_d_nucl+wcorr_nucl*ecorr_nucl &
1226 +wcorr3_nucl*ecorr3_nucl&
1228 +ft(1)*wpepbase*epepbase+wscpho*escpho+wpeppho*epeppho &
1229 + wcatnucl*ecation_nucl&
1230 +elipbond+elipang+eliplj+elipelec+wcat_tran*ecat_prottran+ecation_protang&
1235 eprim=ftprim(6)*evdw_t+ftprim(1)*welec*ees &
1236 +ftprim(1)*wtor*etors+ &
1237 ftprim(3)*wcorr*ecorr+ftprim(4)*wcorr5*ecorr5+ &
1238 ftprim(5)*wcorr6*ecorr6+ftprim(3)*wturn4*eello_turn4+ &
1239 ftprim(2)*wturn3*eello_turn3+ftprim(5)*wturn6*eello_turn6+ &
1240 ftprim(2)*wel_loc*eel_loc+ftprim(2)*wtor_d*etors_d+ &
1241 ftprim(1)*wsccor*esccor +ftprim(1)*wtor_nucl*etors_nucl&
1242 +wtor_d_nucl*ftprim(2)*etors_d_nucl+ftprim(1)*wpepbase*epepbase
1243 ebis=ftbis(1)*welec*ees+ftbis(1)*wtor*etors+ &
1244 ftbis(3)*wcorr*ecorr+ftbis(4)*wcorr5*ecorr5+ &
1245 ftbis(5)*wcorr6*ecorr6+ftbis(3)*wturn4*eello_turn4+ &
1246 ftbis(2)*wturn3*eello_turn3+ftbis(5)*wturn6*eello_turn6+ &
1247 ftbis(2)*wel_loc*eel_loc+ftbis(2)*wtor_d*etors_d+ &
1248 ftbis(1)*wsccor*esccor+ftbis(1)*wtor_nucl*etors_nucl&
1249 +wtor_d_nucl*ftbis(2)*etors_d_nucl+ftbis(1)*wpepbase*epepbase
1251 etot=wsc*(evdw+ft(6)*evdw_t)+wscp*evdw2 &
1252 +ft(1)*welec*(ees+evdw1) &
1253 +wang*ebe+ft(1)*wtor*etors+wscloc*escloc &
1254 +wstrain*ehpb+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5 &
1255 +ft(5)*wcorr6*ecorr6+ft(3)*wturn4*eello_turn4 &
1256 +ft(2)*wturn3*eello_turn3 &
1257 +ft(5)*wturn6*eello_turn6+ft(2)*wel_loc*eel_loc+edihcnstr &
1258 +ft(2)*wtor_d*etors_d+ft(1)*wsccor*esccor &
1259 +wbond*estr+wcatprot*ecation_prot+wcatcat*ecationcation&
1260 +wbond_nucl*estr_nucl+wang_nucl*ebe_nucl&
1261 +wvdwpp_nucl*evdwpp+welpp*eespp+wvdwpsb*evdwpsb+welpsb*eelpsb&
1262 +wvdwsb*evdwsb+welsb*eelsb+wsbloc*esbloc+wtor_nucl*etors_nucl&
1263 *ft(1)+wtor_d_nucl*ft(2)*etors_d_nucl+wcorr_nucl*ecorr_nucl &
1264 +wcorr3_nucl*ecorr3_nucl&
1266 +ft(1)*wpepbase*epepbase+wscpho*escpho+wpeppho*epeppho&
1267 + wcatnucl*ecation_nucl&
1268 +elipbond+elipang+eliplj+elipelec+wcat_tran*ecat_prottran+ecation_protang&
1273 eprim=ftprim(6)*evdw_t+ftprim(1)*welec*(ees+evdw1) &
1274 +ftprim(1)*wtor*etors+ &
1275 ftprim(3)*wcorr*ecorr+ftprim(4)*wcorr5*ecorr5+ &
1276 ftprim(5)*wcorr6*ecorr6+ftprim(3)*wturn4*eello_turn4+ &
1277 ftprim(2)*wturn3*eello_turn3+ftprim(5)*wturn6*eello_turn6+ &
1278 ftprim(2)*wel_loc*eel_loc+ftprim(2)*wtor_d*etors_d+ &
1279 ftprim(1)*wsccor*esccor+ftprim(1)*wtor_nucl*etors_nucl&
1280 +wtor_d_nucl*ftprim(2)*etors_d_nucl+ftprim(1)*wpepbase*epepbase
1281 ebis=ftbis(1)*welec*(ees+evdw1)+ftbis(1)*wtor*etors+ &
1282 ftbis(3)*wcorr*ecorr+ftbis(4)*wcorr5*ecorr5+ &
1283 ftbis(5)*wcorr6*ecorr6+ftbis(3)*wturn4*eello_turn4+ &
1284 ftbis(2)*wturn3*eello_turn3+ftbis(5)*wturn6*eello_turn6+ &
1285 ftbis(2)*wel_loc*eel_loc+ftbis(2)*wtor_d*etors_d+ &
1286 ftbis(1)*wsccor*esccor+ftbis(1)*wtor_nucl*etors_nucl&
1287 +wtor_d_nucl*ftbis(2)*etors_d_nucl*ftbis(1)*wpepbase*epepbase
1290 weight=dexp(-betaT*(etot-potEmin)+entfac(t))
1292 write (iout,*) "iparm",iparm," t",t," betaT",betaT,&
1293 " etot",etot," entfac",entfac(t),&
1294 " weight",weight," ebis",ebis
1296 etot=etot-temper*eprim
1298 sumW_p(k,iparm)=sumW_p(k,iparm)+weight
1299 sumE_p(k,iparm)=sumE_p(k,iparm)+etot*weight
1300 sumEbis_p(k,iparm)=sumEbis_p(k,iparm)+ebis*weight
1301 sumEsq_p(k,iparm)=sumEsq_p(k,iparm)+etot**2*weight
1303 sumQ_p(j,k,iparm)=sumQ_p(j,k,iparm)+q(j,t)*weight
1304 sumQsq_p(j,k,iparm)=sumQsq_p(j,k,iparm)+q(j,t)**2*weight
1305 sumEQ_p(j,k,iparm)=sumEQ_p(j,k,iparm) &
1309 sumW(k,iparm)=sumW(k,iparm)+weight
1310 sumE(k,iparm)=sumE(k,iparm)+etot*weight
1311 sumEbis(k,iparm)=sumEbis(k,iparm)+ebis*weight
1312 sumEsq(k,iparm)=sumEsq(k,iparm)+etot**2*weight
1314 sumQ(j,k,iparm)=sumQ(j,k,iparm)+q(j,t)*weight
1315 sumQsq(j,k,iparm)=sumQsq(j,k,iparm)+q(j,t)**2*weight
1316 sumEQ(j,k,iparm)=sumEQ(j,k,iparm) &
1321 indE = aint(potE(t,iparm)-aint(potEmin))
1322 if (indE.ge.0 .and. indE.le.maxinde) then
1323 if (indE.gt.upindE_p) upindE_p=indE
1324 histE_p(indE)=histE_p(indE)+dexp(-entfac(t))
1328 potEmin=potEmin_all(ib,iparm)
1329 expfac=dexp(-beta_h(ib,iparm)*(etot-potEmin)+entfac(t))
1330 hfin_p(ind,ib)=hfin_p(ind,ib)+ &
1331 dexp(-beta_h(ib,iparm)*(etot-potEmin)+entfac(t))
1333 indrgy=dint((q(nQ+2,t)-rgymin)/deltrgy)
1334 indrms=dint((q(nQ+1,t)-rmsmin)/deltrms)
1335 hrmsrgy_p(indrgy,indrms,ib)= &
1336 hrmsrgy_p(indrgy,indrms,ib)+expfac
1341 potEmin=potEmin_all(ib,iparm)
1342 expfac=dexp(-beta_h(ib,iparm)*(etot-potEmin)+entfac(t))
1343 hfin(ind,ib)=hfin(ind,ib)+ &
1344 dexp(-beta_h(ib,iparm)*(etot-potEmin)+entfac(t))
1346 indrgy=dint((q(nQ+2,t)-rgymin)/deltrgy)
1347 indrms=dint((q(nQ+1,t)-rmsmin)/deltrms)
1348 hrmsrgy(indrgy,indrms,ib)= &
1349 hrmsrgy(indrgy,indrms,ib)+expfac
1355 if (histout) call MPI_Reduce(hfin_p(0,ib),hfin(0,ib),nbin,&
1356 MPI_DOUBLE_PRECISION,MPI_SUM,Master,WHAM_COMM,IERROR)
1358 call MPI_Reduce(hrmsrgy_p(0,0,ib),hrmsrgy(0,0,ib),&
1359 (MaxBinRgy+1)*(nbin_rms+1),MPI_DOUBLE_PRECISION,MPI_SUM,Master,&
1363 call MPI_Reduce(upindE_p,upindE,1,&
1364 MPI_INTEGER,MPI_MAX,Master,WHAM_COMM,IERROR)
1365 call MPI_Reduce(histE_p(0),histE(0),maxindE,&
1366 MPI_DOUBLE_PRECISION,MPI_SUM,Master,WHAM_COMM,IERROR)
1368 if (me1.eq.master) then
1372 write (iout,'(6x,$)')
1373 write (iout,'(f20.2,$)') (1.0d0/(1.987D-3*beta_h(ib,iparm)),&
1377 write (iout,'(/a)') 'Final histograms'
1379 if (nslice.eq.1) then
1380 if (separate_parset) then
1381 write(licz3,"(bz,i3.3)") myparm
1382 histname=prefix(:ilen(prefix))//'_par'//licz3//'.hist'
1384 histname=prefix(:ilen(prefix))//'.hist'
1387 if (separate_parset) then
1388 write(licz3,"(bz,i3.3)") myparm
1389 histname=prefix(:ilen(prefix))//'_par'//licz3// &
1390 '_slice_'//licz2//'.hist'
1392 histname=prefix(:ilen(prefix))//'_slice_'//licz2//'.hist'
1395 #if defined(AIX) || defined(PGI)
1396 open (ihist,file=histname,position='append')
1398 open (ihist,file=histname,access='append')
1406 sumH=sumH+hfin(t,ib)
1408 if (sumH.gt.0.0d0) then
1410 jj = mod(liczbaW,nbin1)
1411 liczbaW=liczbaW/nbin1
1412 write (iout,'(f6.3,$)') dmin+(jj+0.5d0)*delta
1414 write (ihist,'(f6.3,$)') dmin+(jj+0.5d0)*delta
1417 write (iout,'(e20.10,$)') hfin(t,ib)
1418 if (histfile) write (ihist,'(e20.10,$)') hfin(t,ib)
1420 write (iout,'(i5)') iparm
1421 if (histfile) write (ihist,'(i5)') iparm
1428 if (nslice.eq.1) then
1429 if (separate_parset) then
1430 write(licz3,"(bz,i3.3)") myparm
1431 histname=prefix(:ilen(prefix))//"_par"//licz3//'.ent'
1433 histname=prefix(:ilen(prefix))//'.ent'
1436 if (separate_parset) then
1437 write(licz3,"(bz,i3.3)") myparm
1438 histname=prefix(:ilen(prefix))//'par_'//licz3// &
1439 '_slice_'//licz2//'.ent'
1441 histname=prefix(:ilen(prefix))//'_slice_'//licz2//'.ent'
1444 #if defined(AIX) || defined(PGI)
1445 open (ihist,file=histname,position='append')
1447 open (ihist,file=histname,access='append')
1449 write (ihist,'(a)') "# Microcanonical entropy"
1451 write (ihist,'(f8.0,$)') dint(potEmin)+i
1452 if (histE(i).gt.0.0e0) then
1453 write (ihist,'(f15.5,$)') dlog(histE(i))
1455 write (ihist,'(f15.5,$)') 0.0d0
1461 write (iout,*) "Microcanonical entropy"
1463 write (iout,'(f8.0,$)') dint(potEmin)+i
1464 if (histE(i).gt.0.0e0) then
1465 write (iout,'(f15.5,$)') dlog(histE(i))
1467 write (iout,'(f15.5,$)') 0.0d0
1472 if (nslice.eq.1) then
1473 if (separate_parset) then
1474 write(licz3,"(bz,i3.3)") myparm
1475 histname=prefix(:ilen(prefix))//'_par'//licz3//'.rmsrgy'
1477 histname=prefix(:ilen(prefix))//'.rmsrgy'
1480 if (separate_parset) then
1481 write(licz3,"(bz,i3.3)") myparm
1482 histname=prefix(:ilen(prefix))//'_par'//licz3// &
1483 '_slice_'//licz2//'.rmsrgy'
1485 histname=prefix(:ilen(prefix))//'_slice_'//licz2//'.rmsrgy'
1488 #if defined(AIX) || defined(PGI)
1489 open (ihist,file=histname,position='append')
1491 open (ihist,file=histname,access='append')
1495 write(ihist,'(2f8.2,$)') &
1496 rgymin+deltrgy*j,rmsmin+deltrms*i
1498 if (hrmsrgy(j,i,ib).gt.0.0d0) then
1499 write(ihist,'(e14.5,$)') &
1500 -dlog(hrmsrgy(j,i,ib))/beta_h(ib,iparm) &
1503 write(ihist,'(e14.5,$)') 1.0d6
1506 write (ihist,'(i2)') iparm
1514 call MPI_Reduce(hfin_ent_p(0),hfin_ent(0),nbin,&
1515 MPI_DOUBLE_PRECISION,MPI_SUM,Master,WHAM_COMM,IERROR)
1516 call MPI_Reduce(sumW_p(0,1),sumW(0,1),(nGridT+1)*nParmSet,&
1517 MPI_DOUBLE_PRECISION,MPI_SUM,Master,WHAM_COMM,IERROR)
1518 call MPI_Reduce(sumE_p(0,1),sumE(0,1),(nGridT+1)*nParmSet,&
1519 MPI_DOUBLE_PRECISION,MPI_SUM,Master,WHAM_COMM,IERROR)
1520 call MPI_Reduce(sumEbis_p(0,1),sumEbis(0,1),(nGridT+1)*nParmSet,&
1521 MPI_DOUBLE_PRECISION,MPI_SUM,Master,WHAM_COMM,IERROR)
1522 call MPI_Reduce(sumEsq_p(0,1),sumEsq(0,1),(nGridT+1)*nParmSet,&
1523 MPI_DOUBLE_PRECISION,MPI_SUM,Master,WHAM_COMM,IERROR)
1524 call MPI_Reduce(sumQ_p(1,0,1),sumQ(1,0,1),&
1525 MaxQ1*(nGridT+1)*nParmSet,MPI_DOUBLE_PRECISION,MPI_SUM,Master,&
1527 call MPI_Reduce(sumQsq_p(1,0,1),sumQsq(1,0,1),&
1528 MaxQ1*(nGridT+1)*nParmSet,MPI_DOUBLE_PRECISION,MPI_SUM,Master,&
1530 call MPI_Reduce(sumEQ_p(1,0,1),sumEQ(1,0,1),&
1531 MaxQ1*(nGridT+1)*nParmSet,MPI_DOUBLE_PRECISION,MPI_SUM,Master,&
1533 if (me.eq.master) then
1535 write (iout,'(/a)') 'Thermal characteristics of folding'
1536 if (nslice.eq.1) then
1539 nazwa=prefix(:ilen(prefix))//"_slice_"//licz2
1542 if (nparmset.eq.1 .and. .not.separate_parset) then
1543 nazwa=nazwa(:iln)//".thermal"
1544 else if (nparmset.eq.1 .and. separate_parset) then
1545 write(licz3,"(bz,i3.3)") myparm
1546 nazwa=nazwa(:iln)//"_par_"//licz3//".thermal"
1549 if (nparmset.gt.1) then
1550 write(licz3,"(bz,i3.3)") iparm
1551 nazwa=nazwa(:iln)//"_par_"//licz3//".thermal"
1554 if (separate_parset) then
1555 write (iout,'(a,i3)') "Parameter set",myparm
1557 write (iout,'(a,i3)') "Parameter set",iparm
1560 betaT=1.0d0/(1.987D-3*(startGridT+i*delta_T))
1561 if (betaT.ge.beta_h(1,iparm)) then
1562 potEmin=potEmin_all(1,iparm)
1563 else if (betaT.lt.beta_h(nT_h(iparm),iparm)) then
1564 potEmin=potEmin_all(nT_h(iparm),iparm)
1566 do l=1,nT_h(iparm)-1
1567 if (betaT.le.beta_h(l,iparm) .and. &
1568 betaT.gt.beta_h(l+1,iparm)) then
1569 potEmin=potEmin_all(l,iparm)
1575 sumE(i,iparm)=sumE(i,iparm)/sumW(i,iparm)
1576 sumEbis(i,iparm)=(startGridT+i*delta_T)*sumEbis(i,iparm)/ &
1578 sumEsq(i,iparm)=(sumEsq(i,iparm)/sumW(i,iparm) &
1579 -sumE(i,iparm)**2)/(1.987D-3*(startGridT+i*delta_T)**2)
1581 sumQ(j,i,iparm)=sumQ(j,i,iparm)/sumW(i,iparm)
1582 sumQsq(j,i,iparm)=sumQsq(j,i,iparm)/sumW(i,iparm) &
1584 sumEQ(j,i,iparm)=sumEQ(j,i,iparm)/sumW(i,iparm) &
1585 -sumQ(j,i,iparm)*sumE(i,iparm)
1587 sumW(i,iparm)=-dlog(sumW(i,iparm))*(1.987D-3* &
1588 (startGridT+i*delta_T))+potEmin
1589 write (iout,'(f7.1,2f15.5,$)') startGridT+i*delta_T,&
1590 sumW(i,iparm),sumE(i,iparm)
1591 write (iout,'(f10.5,$)') (sumQ(j,i,iparm),j=1,nQ+2)
1592 write (iout,'(e15.5,$)') sumEsq(i,iparm)-sumEbis(i,iparm),&
1593 (sumQsq(j,i,iparm),j=1,nQ+2),(sumEQ(j,i,iparm),j=1,nQ+2)
1595 write (34,'(f7.1,2f15.5,$)') startGridT+i*delta_T,&
1596 sumW(i,iparm),sumE(i,iparm)
1597 write (34,'(f10.5,$)') (sumQ(j,i,iparm),j=1,nQ+2)
1598 write (34,'(e15.5,$)') sumEsq(i,iparm)-sumEbis(i,iparm),&
1599 (sumQsq(j,i,iparm),j=1,nQ+2),(sumEQ(j,i,iparm),j=1,nQ+2)
1607 if (hfin_ent(t).gt.0.0d0) then
1609 jj = mod(liczbaW,nbin1)
1610 write (iout,'(f6.3,e20.10," ent")') dmin+(jj+0.5d0)*delta,&
1612 if (histfile) write (ihist,'(f6.3,e20.10," ent")') &
1613 dmin+(jj+0.5d0)*delta,&
1617 if (histfile) close(ihist)
1621 ! Write data for zscore
1622 if (nslice.eq.1) then
1623 zscname=prefix(:ilen(prefix))//".zsc"
1625 zscname=prefix(:ilen(prefix))//"_slice_"//licz2//".zsc"
1627 #if defined(AIX) || defined(PGI)
1628 open (izsc,file=prefix(:ilen(prefix))//'.zsc',position='append')
1630 open (izsc,file=prefix(:ilen(prefix))//'.zsc',access='append')
1632 write (izsc,'("NQ=",i1," NPARM=",i1)') nQ,nParmSet
1634 write (izsc,'("NT=",i1)') nT_h(iparm)
1636 write (izsc,'("TEMP=",f6.1," NR=",i2," SNK=",$)') &
1637 1.0d0/(beta_h(ib,iparm)*1.987D-3),nR(ib,iparm)
1638 jj = min0(nR(ib,iparm),7)
1639 write (izsc,'(i8,$)') (snk(i,ib,iparm,islice),i=1,jj)
1640 write (izsc,'(a1,$)') (" ",i=22+8*jj+1,79)
1641 write (izsc,'("&")')
1642 if (nR(ib,iparm).gt.7) then
1643 do ii=8,nR(ib,iparm),9
1644 jj = min0(nR(ib,iparm),ii+8)
1645 write (izsc,'(i8,$)') (snk(i,ib,iparm,islice),i=ii,jj)
1646 write (izsc,'(a1,$') (" ",i=(jj-ii+1)*8+1,79)
1647 write (izsc,'("&")')
1650 write (izsc,'("FI=",$)')
1651 jj=min0(nR(ib,iparm),7)
1652 write (izsc,'(f10.5,$)') (fi(i,ib,iparm),i=1,jj)
1653 write (izsc,'(a1,$)') (" ",i=3+10*jj+1,79)
1654 write (izsc,'("&")')
1655 if (nR(ib,iparm).gt.7) then
1656 do ii=8,nR(ib,iparm),9
1657 jj = min0(nR(ib,iparm),ii+8)
1658 write (izsc,'(f10.5,$)') (fi(i,ib,iparm),i=ii,jj)
1659 if (jj.eq.nR(ib,iparm)) then
1662 write (izsc,'(a1,$)') (" ",i=10*(jj-ii+1)+1,79)
1663 write (izsc,'(t80,"&")')
1668 write (izsc,'("KH=",$)')
1669 write (izsc,'(f7.2,$)') (Kh(j,i,ib,iparm),j=1,nQ)
1670 write (izsc,'(" Q0=",$)')
1671 write (izsc,'(f7.5,$)') (q0(j,i,ib,iparm),j=1,nQ)
1682 end subroutine WHAMCALC
1683 !-----------------------------------------------------------------------------
1684 end module wham_calc