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),weimax_(0:ngridT)
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 weimax(0:nGridT,Max_Parm)
103 real(kind=8) :: fT(6),fTprim(6),fTbis(6),quot,quotl1,quotl,kfacl,&
104 eprim,ebis,temper,kfac=2.4d0,T0=300.0d0,startGridT=200.0d0,&
105 eplus,eminus,logfac,tanhT,tt
106 real(kind=8) :: etot,evdw,evdw_t,evdw2,ees,evdw1,ebe,etors,&
107 escloc,ehpb,ecorr,ecorr5,ecorr6,eello_turn4,eello_turn3,&
108 eello_turn6,eel_loc,edihcnstr,etors_d,estr,evdw2_14,esccor, &
109 ecationcation,ecation_prot, evdwpp,eespp ,evdwpsb,eelpsb, &
110 evdwsb, eelsb, estr_nucl,ebe_nucl,esbloc,etors_nucl,etors_d_nucl,&
111 ecorr_nucl, ecorr3_nucl,epeppho, escpho, epepbase,escbase,ecation_nucl,&
112 elipbond,elipang,eliplj,elipelec,ecat_prottran,ecation_protang,eliptran,&
116 integer :: ind_point(maxpoint),upindE,indE
117 character(len=16) :: plik
118 character(len=1) :: licz1
119 character(len=2) :: licz2
120 character(len=3) :: licz3
121 character(len=128) :: nazwa
126 write(licz2,'(bz,i2.2)') islice
128 write (iout,'(//80(1h-)/"Solving WHAM equations for slice",&
129 i2/80(1h-)//)') islice
130 write (iout,*) "delta",delta," nbin1",nbin1
131 write (iout,*) "MaxN",MaxN," MaxQ",MaxQ," MaHdim",MaxHdim
138 potEmin_all(j,i)=1.0d15
155 if (potE(i,j).le.potEmin) potEmin=potE(i,j)
157 if (q(nQ+1,i).lt.rmsmin) rmsmin=q(nQ+1,i)
158 if (q(nQ+1,i).gt.rmsmax) rmsmax=q(nQ+1,i)
159 if (q(nQ+2,i).lt.rgymin) rgymin=q(nQ+2,i)
160 if (q(nQ+2,i).gt.rgymax) rgymax=q(nQ+2,i)
163 ind=(q(j,i)-dmin+1.0d-8)/delta
165 ind_point(i)=ind_point(i)+ind
167 ind_point(i)=ind_point(i)+nbin1**(j-1)*ind
169 ! write (iout,*) "i",i," j",j," q",q(j,i)," ind_point",
172 if (ind_point(i).lt.0 .or. ind_point(i).gt.MaxHdim) then
173 write (iout,*) "Error - index exceeds range for point",i,&
174 " q=",q(j,i)," ind",ind_point(i)
176 write (iout,*) "Processor",me1
178 call MPI_Abort(MPI_COMM_WORLD, Ierror, Errcode )
183 if (ind_point(i).gt.tmax) tmax=ind_point(i)
184 htot(ind_point(i))=htot(ind_point(i))+1
186 write (iout,*) "i",i,"q",(q(j,i),j=1,nQ)," ind",ind_point(i),&
187 " htot",htot(ind_point(i))
194 write (iout,'(a)') "Numbers of counts in Q bins"
196 if (htot(t).gt.0) then
197 write (iout,'(i15,$)') t
200 jj = mod(liczbaW,nbin1)
201 liczbaW=liczbaW/nbin1
202 write (iout,'(i5,$)') jj
204 write (iout,'(i8)') htot(t)
208 write (iout,'(a,i3)') "Number of data points for parameter set",&
210 write (iout,'(i7,$)') ((snk(m,ib,iparm,islice),m=1,nr(ib,iparm)),&
212 write (iout,'(i8)') stot(islice)
218 call MPI_AllReduce(tmax,tmax_t,1,MPI_INTEGER,MPI_MAX,&
221 ! call MPI_AllReduce(potEmin,potEmin_t,1,MPI_DOUBLE_PRECISION,&
222 ! MPI_MIN,WHAM_COMM,IERROR) !????
223 call MPI_AllReduce(rmsmin,rmsmin_t,1,MPI_DOUBLE_PRECISION,&
224 MPI_MIN,WHAM_COMM,IERROR)
225 call MPI_AllReduce(rmsmax,rmsmax_t,1,MPI_DOUBLE_PRECISION,&
226 MPI_MAX,WHAM_COMM,IERROR)
227 call MPI_AllReduce(rgymin,rgymin_t,1,MPI_DOUBLE_PRECISION,&
228 MPI_MIN,WHAM_COMM,IERROR)
229 call MPI_AllReduce(rgymax,rgymax_t,1,MPI_DOUBLE_PRECISION,&
230 MPI_MAX,WHAM_COMM,IERROR)
231 ! potEmin=potEmin_t !/2 try now??
236 write (iout,*) "potEmin",potEmin
238 rmsmin=deltrms*dint(rmsmin/deltrms)
239 rmsmax=deltrms*dint(rmsmax/deltrms)
240 rgymin=deltrms*dint(rgymin/deltrgy)
241 rgymax=deltrms*dint(rgymax/deltrgy)
242 nbin_rms=(rmsmax-rmsmin)/deltrms
243 nbin_rgy=(rgymax-rgymin)/deltrgy
244 write (iout,*) "rmsmin",rmsmin," rmsmax",rmsmax," rgymin",rgymin,&
245 " rgymax",rgymax," nbin_rms",nbin_rms," nbin_rgy",nbin_rgy
252 write (iout,*) "nFi",nFi
253 ! Compute the Boltzmann factor corresponing to restrain potentials in different
260 ! write (9,'(3i5,f10.5)') i,(iparm,potE(i,iparm),iparm=1,nParmSet)
263 write (iout,'(2i5,21f8.2)') i,iparm,&
264 (enetb(k,i,iparm),k=1,21)
266 call restore_parm(iparm)
268 write (iout,*) wsc,wscp,welec,wvdwpp,wang,wtor,wscloc,&
269 wcorr,wcorr5,wcorr6,wturn4,wturn3,wturn6,wel_loc,&
270 wtor_d,wsccor,wbond,wcatcat
273 !el old rascale weights
275 ! if (rescale_modeW.eq.1) then
276 ! quot=1.0d0/(beta_h(ib,iparm)*1.987D-3*T0)
283 ! fT(l)=kfacl/(kfacl-1.0d0+quotl)
286 ! tt = 1.0d0/(beta_h(ib,iparm)*1.987D-3)
287 ! ft(6)=(320.0d0+80.0d0*dtanh((tt-320.0d0)/80.0d0))/320.0d0
288 !#elif defined(FUNCT)
289 ! ft(6)=1.0d0/(beta_h(ib,iparm)*1.987D-3*T0)
293 ! else if (rescale_modeW.eq.2) then
294 ! quot=1.0d0/(T0*beta_h(ib,iparm)*1.987D-3)
298 ! fT(l)=1.12692801104297249644d0/ &
299 ! dlog(dexp(quotl)+dexp(-quotl))
302 ! tt = 1.0d0/(beta_h(ib,iparm)*1.987D-3)
303 ! ft(6)=(320.0d0+80.0d0*dtanh((tt-320.0d0)/80.0d0))/320.0d0
304 !#elif defined(FUNCT)
305 ! ft(6)=1.0d0/(beta_h(ib,iparm)*1.987D-3*T0)
309 ! write (iout,*) 1.0d0/(beta_h(ib,iparm)*1.987D-3),ft
310 ! else if (rescale_modeW.eq.0) then
315 ! write (iout,*) "Error in WHAM_CALC: wrong RESCALE_MODE",&
320 ! el end old rescale weights
321 call rescale_weights(1.0d0/(beta_h(ib,iparm)*1.987D-3))
323 ! call etot(enetb(0,i,iparm))
324 evdw=enetb(1,i,iparm)
325 ! evdw_t=enetb(21,i,iparm)
326 evdw_t=enetb(20,i,iparm)
328 ! evdw2_14=enetb(17,i,iparm)
329 evdw2_14=enetb(18,i,iparm)
330 evdw2=enetb(2,i,iparm)+evdw2_14
332 evdw2=enetb(2,i,iparm)
337 evdw1=enetb(16,i,iparm)
342 ecorr=enetb(4,i,iparm)
343 ecorr5=enetb(5,i,iparm)
344 ecorr6=enetb(6,i,iparm)
345 eel_loc=enetb(7,i,iparm)
346 eello_turn3=enetb(8,i,iparm)
347 eello_turn4=enetb(9,i,iparm)
348 eello_turn6=enetb(10,i,iparm)
349 ebe=enetb(11,i,iparm)
350 escloc=enetb(12,i,iparm)
351 etors=enetb(13,i,iparm)
352 etors_d=enetb(14,i,iparm)
353 ehpb=enetb(15,i,iparm)
354 ! estr=enetb(18,i,iparm)
355 estr=enetb(17,i,iparm)
356 ! esccor=enetb(19,i,iparm)
357 esccor=enetb(21,i,iparm)
358 ! edihcnstr=enetb(20,i,iparm)
359 edihcnstr=enetb(19,i,iparm)
360 ecationcation=enetb(41,i,iparm)
361 ecation_prot=enetb(42,i,iparm)
362 evdwpp = enetb(26,i,iparm)
363 eespp = enetb(27,i,iparm)
364 evdwpsb = enetb(28,i,iparm)
365 eelpsb = enetb(29,i,iparm)
366 evdwsb = enetb(30,i,iparm)
367 eelsb = enetb(31,i,iparm)
368 estr_nucl = enetb(32,i,iparm)
369 ebe_nucl = enetb(33,i,iparm)
370 esbloc = enetb(34,i,iparm)
371 etors_nucl = enetb(35,i,iparm)
372 etors_d_nucl = enetb(36,i,iparm)
373 ecorr_nucl = enetb(37,i,iparm)
374 ecorr3_nucl = enetb(38,i,iparm)
375 epeppho= enetb(49,i,iparm)
376 escpho= enetb(48,i,iparm)
377 epepbase= enetb(47,i,iparm)
378 escbase= enetb(46,i,iparm)
379 ecation_nucl= enetb(50,i,iparm)
380 elipbond=enetb(52,i,iparm)
381 elipang=enetb(53,i,iparm)
382 eliplj=enetb(54,i,iparm)
383 elipelec=enetb(55,i,iparm)
384 ecat_prottran=enetb(56,i,iparm)
385 ecation_protang=enetb(57,i,iparm)
386 eliptran=enetb(22,i,iparm)
387 ehomology_constr=enetb(51,i,iparm)
389 write (iout,'(3i5,6f5.2,15f12.3)') i,ib,iparm,(ft(l),l=1,6),&
390 evdw+evdw_t,evdw2,ees,evdw1,ecorr,eel_loc,estr,ebe,escloc,&
391 etors,etors_d,eello_turn3,eello_turn4,esccor,ecationcation,&
396 ! etot=wsc*(evdw+ft(6)*evdw_t)+wscp*evdw2+ft(1)*welec*ees &
398 ! +wang*ebe+ft(1)*wtor*etors+wscloc*escloc &
399 ! +wstrain*ehpb+nss*ebr+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5 &
400 ! +ft(5)*wcorr6*ecorr6+ft(3)*wturn4*eello_turn4 &
401 ! +ft(2)*wturn3*eello_turn3 &
402 ! +ft(5)*wturn6*eello_turn6+ft(2)*wel_loc*eel_loc &
403 ! +edihcnstr+ft(2)*wtor_d*etors_d+ft(1)*wsccor*esccor &
406 ! etot=wsc*(evdw+ft(6)*evdw_t)+wscp*evdw2 &
407 ! +ft(1)*welec*(ees+evdw1) &
408 ! +wang*ebe+ft(1)*wtor*etors+wscloc*escloc &
409 ! +wstrain*ehpb+nss*ebr+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5 &
410 ! +ft(5)*wcorr6*ecorr6+ft(3)*wturn4*eello_turn4 &
411 ! +ft(2)*wturn3*eello_turn3 &
412 ! +ft(5)*wturn6*eello_turn6+ft(2)*wel_loc*eel_loc+edihcnstr &
413 ! +ft(2)*wtor_d*etors_d+ft(1)*wsccor*esccor &
418 etot=wsc*evdw+wscp*evdw2+welec*ees &
420 +wang*ebe+wtor*etors+wscloc*escloc &
421 +wstrain*ehpb+wcorr*ecorr+wcorr5*ecorr5 &
422 +wcorr6*ecorr6+wturn4*eello_turn4 &
423 +wturn3*eello_turn3 &
424 +wturn6*eello_turn6+wel_loc*eel_loc &
425 +edihcnstr+wtor_d*etors_d+wsccor*esccor &
426 +wbond*estr+wcatprot*ecation_prot+wcatcat*ecationcation&
427 +wbond_nucl*estr_nucl+wang_nucl*ebe_nucl&
428 +wvdwpp_nucl*evdwpp+welpp*eespp+wvdwpsb*evdwpsb+welpsb*eelpsb&
429 +wvdwsb*evdwsb+welsb*eelsb+wsbloc*esbloc+wtor_nucl*etors_nucl&
430 +wtor_d_nucl*etors_d_nucl+wcorr_nucl*ecorr_nucl+wcorr3_nucl*ecorr3_nucl&
432 +wpepbase*epepbase+wscpho*escpho+wpeppho*epeppho+wcatnucl*ecation_nucl&
433 +elipbond+elipang+eliplj+elipelec+wcat_tran*ecat_prottran+ecation_protang&
438 etot=wsc*evdw+wscp*evdw2 &
440 +wang*ebe+wtor*etors+wscloc*escloc &
441 +wstrain*ehpb+wcorr*ecorr+wcorr5*ecorr5 &
442 +wcorr6*ecorr6+wturn4*eello_turn4 &
443 +wturn3*eello_turn3 &
444 +wturn6*eello_turn6+wel_loc*eel_loc+edihcnstr &
445 +wtor_d*etors_d+wsccor*esccor &
446 +wbond*estr+wcatprot*ecation_prot+wcatcat*ecationcation&
447 +wbond_nucl*estr_nucl+wang_nucl*ebe_nucl&
448 +wvdwpp_nucl*evdwpp+welpp*eespp+wvdwpsb*evdwpsb+welpsb*eelpsb&
449 +wvdwsb*evdwsb+welsb*eelsb+wsbloc*esbloc+wtor_nucl*etors_nucl&
450 +wtor_d_nucl*etors_d_nucl+wcorr_nucl*ecorr_nucl+wcorr3_nucl*ecorr3_nucl&
452 +wpepbase*epepbase+wscpho*escpho+wpeppho*epeppho+wcatnucl*ecation_nucl&
453 +elipbond+elipang+eliplj+elipelec+wcat_tran*ecat_prottran+ecation_protang&
462 write (iout,*) i,iparm,1.0d0/(beta_h(ib,iparm)*1.987D-3),&
466 if (iparm.eq.1 .and. ib.eq.1) then
467 write (iout,*)"Conformation",i
470 energia(k)=enetb(k,i,iparm)
472 ! call enerprint(energia(0),fT)
473 call enerprint(energia(0))
476 write (iout,*) "homol_nset",homol_nset,nR(ib,iparm)
477 write (iout,*) "waga_homology", waga_homology
479 write (iout,*) "homol_nset",homol_nset,nR(ib,iparm)
481 if (homol_nset.gt.1) then
484 Econstr=waga_homology(kk)*ehomology_constr
486 -beta_h(ib,iparm)*(etot+Econstr)
488 write (iout,'(4i5,4e15.5)') i,kk,ib,iparm, &
489 etot,Econstr,v(i,kk,ib,iparm)
495 etot=etot+ehomology_constr
502 Econstr=Econstr+Kh(j,kk,ib,iparm) &
503 *(ddW-q0(j,kk,ib,iparm))**2
506 -beta_h(ib,iparm)*(etot-potEmin+Econstr)
508 write (iout,'(4i5,4e15.5)') i,kk,ib,iparm,&
509 etot,potEmin,etot-potEmin,v(i,kk,ib,iparm)
516 ! Simple iteration to calculate free energies corresponding to all simulation
520 ! Compute new free-energy values corresponding to the righ-hand side of the
521 ! equation and their derivatives.
522 write (iout,*) "------------------------fi"
532 vf=v(t,l,k,i)+f(l,k,i)
533 if (vf.gt.vmax) vmax=vf
541 aux=f(l,k,i)+v(t,l,k,i)-vmax
542 if (aux.gt.-200.0d0) &
543 denom=denom+snk(l,k,i,islice)*dexp(aux)
547 entfac(t)=-dlog(denom)-vmax
548 if (entfac(t).lt.entfac_min) entfac_min=entfac(t)
550 write (iout,*) t,"vmax",vmax," denom",denom,"entfac",entfac(t)
555 do ii=1,nR(iib,iparm)
557 fi_p_min(ii,iib,iparm)=-1.0d5
559 ! fi_p(ii,iib,iparm)=fi_p(ii,iib,iparm) &
560 ! +dexp(v(t,ii,iib,iparm)+entfac(t))
561 aux=v(t,ii,iib,iparm)+entfac(t)
562 if (aux.gt.fi_p_min(ii,iib,iparm)) fi_p_min(ii,iib,iparm)=aux
566 write (iout,'(4i5,4e15.5)') t,ii,iib,iparm,&
567 v(t,ii,iib,iparm),entfac(t),fi_p(ii,iib,iparm),fi_p_min(ii,iib,iparm)
572 ! fi(ii,iib,iparm)=0.0d0
574 ! fi(ii,iib,iparm)=fi(ii,iib,iparm) &
575 ! +dexp(v(t,ii,iib,iparm)+entfac(t))
576 aux=v(t,ii,iib,iparm)+entfac(t)
577 if (aux.gt.fi_min(ii,iib,iparm))
578 & fi_min(ii,iib,iparm)=aux
588 write (iout,*) "fi before MPI_Reduce me",me,' master',master
590 do ib=1,nT_h(nparmset)
591 write (iout,*) "iparm",iparm," ib",ib
592 write (iout,*) "beta=",beta_h(ib,iparm)
593 write (iout,'(8e15.5)') (fi_p(i,ib,iparm),i=1,nR(ib,iparm))
594 write (iout,'(8e15.5)') (fi_p_min(i,ib,iparm),i=1,nR(ib,iparm))
598 call MPI_AllReduce(fi_p_min,fi_min,MaxR*MaxT_h*nParmSet, &
599 MPI_DOUBLE_PRECISION,MPI_MAX,WHAM_COMM,IERROR)
602 write (iout,*) "fi_min after AllReduce"
605 write (iout,*) (i,j,k,fi_min(k,j,i),k=1,nR(j,i))
613 do ii=1,nR(iib,iparm)
615 fi_p(ii,iib,iparm)=0.0d0
617 fi_p(ii,iib,iparm)=fi_p(ii,iib,iparm) &
618 +dexp(v(t,ii,iib,iparm)+entfac(t)-fi_min(ii,iib,iparm))
620 write (iout,'(4i5,4e15.5)') t,ii,iib,iparm, &
621 v(t,ii,iib,iparm),entfac(t),fi_min(ii,iib,iparm), &
626 fi(ii,iib,iparm)=0.0d0
628 fi(ii,iib,iparm)=fi(ii,iib,iparm) &
629 +dexp(v(t,ii,iib,iparm)+entfac(t)-fi_min(ii,iib,iparm))
638 write (iout,*) "fi before MPI_Reduce me",me,' master',master
640 do ib=1,nT_h(nparmset)
641 write (iout,*) "iparm",iparm," ib",ib
642 write (iout,*) "beta=",beta_h(ib,iparm)
643 write (iout,'(8e15.5)') (fi_p(i,ib,iparm),i=1,nR(ib,iparm))
648 write (iout,*) "REDUCE size",maxR,MaxT_h,nParmSet, &
650 write (iout,*) "MPI_COMM_WORLD",MPI_COMM_WORLD, &
651 " WHAM_COMM",WHAM_COMM
654 write (iout,*) "REDUCE size",maxR,MaxT_h,nParmSet,&
656 write (iout,*) "MPI_COMM_WORLD",MPI_COMM_WORLD,&
657 " WHAM_COMM",WHAM_COMM
658 call MPI_Reduce(fi_p(1,1,1),fi(1,1,1),maxR*MaxT_h*nParmSet,&
659 MPI_DOUBLE_PRECISION,&
660 MPI_SUM,Master,WHAM_COMM,IERROR)
662 write (iout,*) "fi after MPI_Reduce nparmset",nparmset
664 write (iout,*) "iparm",iparm
666 write (iout,*) "beta=",beta_h(ib,iparm)
667 write (iout,'(8e15.5)') (fi(i,ib,iparm),i=1,nR(ib,iparm))
671 if (me1.eq.Master) then
677 fi(i,ib,iparm)=-dlog(fi(i,ib,iparm))-fi_min(i,ib,iparm)
678 avefi=avefi+fi(i,ib,iparm)
684 write (iout,*) "Parameter set",iparm
686 write (iout,*) "beta=",beta_h(ib,iparm)
688 fi(i,ib,iparm)=fi(i,ib,iparm)-avefi
690 write (iout,'(8f10.5)') (fi(i,ib,iparm),i=1,nR(ib,iparm))
691 write (iout,'(8f10.5)') (f(i,ib,iparm),i=1,nR(ib,iparm))
695 ! Compute the norm of free-energy increments.
700 finorm=finorm+dabs(fi(i,ib,iparm)-f(i,ib,iparm))
701 f(i,ib,iparm)=fi(i,ib,iparm)
706 write (iout,*) 'Iteration',iter,' finorm',finorm
710 call MPI_Bcast(f(1,1,1),MaxR*MaxT_h*nParmSet,&
711 MPI_DOUBLE_PRECISION,Master,&
713 call MPI_Bcast(finorm,1,MPI_DOUBLE_PRECISION,Master,&
716 ! Exit, if the increment norm is smaller than pre-assigned tolerance.
717 if (finorm.lt.fimin) then
718 write (iout,*) 'Iteration converged'
725 ! Now, put together the histograms from all simulations, in order to get the
726 ! unbiased total histogram.
727 !C Determine the minimum free energies
733 !c write (9,'(3i5,f10.5)') i,(iparm,potE(i,iparm),iparm=1,nParmSet)
736 write (iout,'(2i5,21f8.2)') i,iparm, &
737 (enetb(k,i,iparm),k=1,26)
739 call restore_parm(iparm)
741 write (iout,*) wsc,wscp,welec,wvdwpp,wang,wtor,wscloc, &
742 wcorr,wcorr5,wcorr6,wturn4,wturn3,wturn6,wel_loc, &
746 if (rescale_modeW.eq.1) then
747 quot=1.0d0/(beta_h(ib,iparm)*1.987D-3*T0)
754 fT(l)=kfacl/(kfacl-1.0d0+quotl)
757 tt = 1.0d0/(beta_h(ib,iparm)*1.987D-3)
758 ft(6)=(320.0d0+80.0d0*dtanh((tt-320.0d0)/80.0d0))/320.0d0
760 ft(6)=1.0d0/(beta_h(ib,iparm)*1.987D-3*T0)
764 else if (rescale_modeW.eq.2) then
765 quot=1.0d0/(T0*beta_h(ib,iparm)*1.987D-3)
769 fT(l)=1.12692801104297249644d0/ &
770 dlog(dexp(quotl)+dexp(-quotl))
773 tt = 1.0d0/(beta_h(ib,iparm)*1.987D-3)
774 ft(6)=(320.0d0+80.0d0*dtanh((tt-320.0d0)/80.0d0))/320.0d0
776 ft(6)=1.0d0/(beta_h(ib,iparm)*1.987D-3*T0)
780 !c write (iout,*) 1.0d0/(beta_h(ib,iparm)*1.987D-3),ft
781 else if (rescale_modeW.eq.0) then
786 write (iout,*) "Error in WHAM_CALC: wrong RESCALE_MODE", &
791 evdw=enetb(1,i,iparm)
792 ! evdw_t=enetb(21,i,iparm)
793 evdw_t=enetb(20,i,iparm)
795 ! evdw2_14=enetb(17,i,iparm)
796 evdw2_14=enetb(18,i,iparm)
797 evdw2=enetb(2,i,iparm)+evdw2_14
799 evdw2=enetb(2,i,iparm)
804 evdw1=enetb(16,i,iparm)
809 ecorr=enetb(4,i,iparm)
810 ecorr5=enetb(5,i,iparm)
811 ecorr6=enetb(6,i,iparm)
812 eel_loc=enetb(7,i,iparm)
813 eello_turn3=enetb(8,i,iparm)
814 eello_turn4=enetb(9,i,iparm)
815 eello_turn6=enetb(10,i,iparm)
816 ebe=enetb(11,i,iparm)
817 escloc=enetb(12,i,iparm)
818 etors=enetb(13,i,iparm)
819 etors_d=enetb(14,i,iparm)
820 ehpb=enetb(15,i,iparm)
821 ! estr=enetb(18,i,iparm)
822 estr=enetb(17,i,iparm)
823 ! esccor=enetb(19,i,iparm)
824 esccor=enetb(21,i,iparm)
825 ! edihcnstr=enetb(20,i,iparm)
826 edihcnstr=enetb(19,i,iparm)
827 ! ehomology_constr=enetb(22,i,iparm)
828 ! esaxs=enetb(26,i,iparm)
829 ecationcation=enetb(41,i,iparm)
830 ecation_prot=enetb(42,i,iparm)
831 evdwpp = enetb(26,i,iparm)
832 eespp = enetb(27,i,iparm)
833 evdwpsb = enetb(28,i,iparm)
834 eelpsb = enetb(29,i,iparm)
835 evdwsb = enetb(30,i,iparm)
836 eelsb = enetb(31,i,iparm)
837 estr_nucl = enetb(32,i,iparm)
838 ebe_nucl = enetb(33,i,iparm)
839 esbloc = enetb(34,i,iparm)
840 etors_nucl = enetb(35,i,iparm)
841 etors_d_nucl = enetb(36,i,iparm)
842 ecorr_nucl = enetb(37,i,iparm)
843 ecorr3_nucl = enetb(38,i,iparm)
844 epeppho= enetb(49,i,iparm)
845 escpho= enetb(48,i,iparm)
846 epepbase= enetb(47,i,iparm)
847 escbase= enetb(46,i,iparm)
848 ecation_nucl=enetb(50,i,iparm)
849 elipbond=enetb(52,i,iparm)
850 elipang=enetb(53,i,iparm)
851 eliplj=enetb(54,i,iparm)
852 elipelec=enetb(55,i,iparm)
853 ecat_prottran=enetb(56,i,iparm)
854 ecation_protang=enetb(57,i,iparm)
855 eliptran=enetb(22,i,iparm)
856 ehomology_constr=enetb(51,i,iparm)
858 etot=wsc*(evdw+ft(6)*evdw_t)+wscp*evdw2+ft(1)*welec*ees &
860 +wang*ebe+ft(1)*wtor*etors+wscloc*escloc &
861 +wstrain*ehpb+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5 &
862 +ft(5)*wcorr6*ecorr6+ft(3)*wturn4*eello_turn4 &
863 +ft(2)*wturn3*eello_turn3 &
864 +ft(5)*wturn6*eello_turn6+ft(2)*wel_loc*eel_loc &
865 +edihcnstr+ft(2)*wtor_d*etors_d+ft(1)*wsccor*esccor &
866 +wbond*estr+wcatprot*ecation_prot+wcatcat*ecationcation &
867 +wbond_nucl*estr_nucl+wang_nucl*ebe_nucl&
868 +wvdwpp_nucl*evdwpp+welpp*eespp+wvdwpsb*evdwpsb+welpsb*eelpsb&
869 +wvdwsb*evdwsb+welsb*eelsb+wsbloc*esbloc+wtor_nucl*etors_nucl&
870 *ft(1)+wtor_d_nucl*ft(2)*etors_d_nucl+wcorr_nucl*ecorr_nucl &
871 +wcorr3_nucl*ecorr3_nucl&
873 +ft(1)*wpepbase*epepbase+wscpho*escpho+wpeppho*epeppho+ wcatnucl*ecation_nucl&
874 +elipbond+elipang+eliplj+elipelec+wcat_tran*ecat_prottran+ecation_protang&
879 etot=wsc*(evdw+ft(6)*evdw_t)+wscp*evdw2 &
880 +ft(1)*welec*(ees+evdw1) &
881 +wang*ebe+ft(1)*wtor*etors+wscloc*escloc &
882 +wstrain*ehpb+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5 &
883 +ft(5)*wcorr6*ecorr6+ft(3)*wturn4*eello_turn4 &
884 +ft(2)*wturn3*eello_turn3 &
885 +ft(5)*wturn6*eello_turn6+ft(2)*wel_loc*eel_loc+edihcnstr &
886 +ft(2)*wtor_d*etors_d+ft(1)*wsccor*esccor &
887 +wbond*estr+wcatprot*ecation_prot+wcatcat*ecationcation&
888 +wbond_nucl*estr_nucl+wang_nucl*ebe_nucl&
889 +wvdwpp_nucl*evdwpp+welpp*eespp+wvdwpsb*evdwpsb+welpsb*eelpsb&
890 +wvdwsb*evdwsb+welsb*eelsb+wsbloc*esbloc+wtor_nucl*etors_nucl&
891 *ft(1)+wtor_d_nucl*ft(2)*etors_d_nucl+wcorr_nucl*ecorr_nucl &
892 +wcorr3_nucl*ecorr3_nucl&
894 +ft(1)*wpepbase*epepbase+wscpho*escpho+wpeppho*epeppho+ wcatnucl*ecation_nucl&
895 +elipbond+elipang+eliplj+elipelec+wcat_tran*ecat_prottran+ecation_protang&
901 write (iout,*) "WTF,",etot,potEmin_all(ib,iparm),entfac(i)/beta_h(ib,iparm)
902 etot=etot-entfac(i)/beta_h(ib,iparm)
903 if(etot.lt.potEmin_all(ib,iparm)) potEmin_all(ib,iparm)=etot
909 write (iout,*) "The potEmin array before reduction"
911 write (iout,*) "Parameter set",i
913 write (iout,*) j,PotEmin_all(j,i)
916 write (iout,*) "potEmin_min",potEmin_min
919 !C Determine the minimum energes for all parameter sets and temperatures
920 call MPI_AllReduce(potEmin_all(1,1),potEmin_t_all(1,1), &
921 maxT_h*nParmSet,MPI_DOUBLE_PRECISION,MPI_MIN,WHAM_COMM,IERROR)
924 potEmin_all(j,i)=potEmin_t_all(j,i)
928 potEmin_min=potEmin_all(1,1)
931 if (potEmin_all(j,i).lt.potEmin_min) &
932 potEmin_min=potEmin_all(j,i)
936 write (iout,*) "The potEmin array"
938 write (iout,*) "Parameter set",i
940 write (iout,*) j,1.0d0/(1.987d-3*beta_h(j,i)), &
944 write (iout,*) "potEmin_min",potEmin_min
956 write (iout,*) "--------------hist"
960 sumW_p(i,iparm)=0.0d0
961 sumE_p(i,iparm)=0.0d0
962 sumEbis_p(i,iparm)=0.0d0
963 sumEsq_p(i,iparm)=0.0d0
965 sumQ_p(j,i,iparm)=0.0d0
966 sumQsq_p(j,i,iparm)=0.0d0
967 sumEQ_p(j,i,iparm)=0.0d0
977 sumEbis(i,iparm)=0.0d0
978 sumEsq(i,iparm)=0.0d0
980 sumQ(j,i,iparm)=0.0d0
981 sumQsq(j,i,iparm)=0.0d0
982 sumEQ(j,i,iparm)=0.0d0
988 ! 8/26/05 entropy distribution
993 !! ent=-dlog(entfac(t))
995 ! if (ent.lt.entmin_p) entmin_p=ent
996 ! if (ent.gt.entmax_p) entmax_p=ent
998 ! write (iout,*) "entmin",entmin_p," entmax",entmax_p
999 !! write (iout,*) "entmin_p",entmin_p," entmax_p",entmax_p
1001 ! call MPI_Allreduce(entmin_p,entmin,1,MPI_DOUBLE_PRECISION,MPI_MIN,&
1003 ! call MPI_Allreduce(entmax_p,entmax,1,MPI_DOUBLE_PRECISION,MPI_MAX,&
1005 ! write (iout,*) "entmin",entmin," entmax",entmax
1006 ! write (iout,*) "entmin_p",entmin_p," entmax_p",entmax_p
1007 ! ientmax=entmax-entmin
1008 !iientmax=entmax-entmin !el
1009 !write (iout,*) "ientmax",ientmax,entmax,entmin
1010 !write (iout,*) "iientmax",iientmax
1011 ! if (ientmax.gt.2000) ientmax=2000
1012 ! write (iout,*) "entmin",entmin," entmax",entmax," ientmax",ientmax
1014 ! do t=1,scount(me1)
1015 !! ient=-dlog(entfac(t))-entmin
1016 ! ient=entfac(t)-entmin
1017 ! if (ient.le.2000) histent_p(ient)=histent_p(ient)+1
1019 ! call MPI_Allreduce(histent_p(0),histent(0),ientmax+1,MPI_INTEGER,&
1020 ! MPI_SUM,WHAM_COMM,IERROR)
1021 ! if (me1.eq.Master) then
1022 ! write (iout,*) "Entropy histogram"
1024 ! write(iout,'(f15.4,i10)') entmin+i,histent(i)
1030 ! do t=1,ntot(islice)
1032 ! if (ent.lt.entmin) entmin=ent
1033 ! if (ent.gt.entmax) entmax=ent
1035 ! ientmax=-dlog(entmax)-entmin
1036 ! if (ientmax.gt.2000) ientmax=2000
1037 ! do t=1,ntot(islice)
1038 ! ient=entfac(t)-entmin
1039 ! if (ient.le.2000) histent(ient)=histent(ient)+1
1041 ! write (iout,*) "Entropy histogram"
1043 ! write(iout,'(2f15.4)') entmin+i,histent(i)
1048 write (iout,*) "me1",me1," scount",scount(me1) !d
1074 hrmsrgy(j,i,ib)=0.0d0
1076 hrmsrgy_p(j,i,ib)=0.0d0
1088 hfin_ent_p(ind)=hfin_ent_p(ind)+dexp(entfac(t))
1090 hfin_ent(ind)=hfin_ent(ind)+dexp(entfac(t))
1092 ! write (iout,'(2i5,20f8.2)') "debug",t,t,(enetb(k,t,iparm),k=1,21)
1093 call restore_parm(iparm)
1094 ! evdw=enetb(21,t,iparm)
1095 evdw=enetb(20,t,iparm)
1096 evdw_t=enetb(1,t,iparm)
1098 ! evdw2_14=enetb(17,t,iparm)
1099 evdw2_14=enetb(18,t,iparm)
1100 evdw2=enetb(2,t,iparm)+evdw2_14
1102 evdw2=enetb(2,t,iparm)
1106 ees=enetb(3,t,iparm)
1107 evdw1=enetb(16,t,iparm)
1109 ees=enetb(3,t,iparm)
1112 ecorr=enetb(4,t,iparm)
1113 ecorr5=enetb(5,t,iparm)
1114 ecorr6=enetb(6,t,iparm)
1115 eel_loc=enetb(7,t,iparm)
1116 eello_turn3=enetb(8,t,iparm)
1117 eello_turn4=enetb(9,t,iparm)
1118 eello_turn6=enetb(10,t,iparm)
1119 ebe=enetb(11,t,iparm)
1120 escloc=enetb(12,t,iparm)
1121 etors=enetb(13,t,iparm)
1122 etors_d=enetb(14,t,iparm)
1123 ehpb=enetb(15,t,iparm)
1124 ! estr=enetb(18,t,iparm)
1125 estr=enetb(17,t,iparm)
1126 ! esccor=enetb(19,t,iparm)
1127 esccor=enetb(21,t,iparm)
1128 ! edihcnstr=enetb(20,t,iparm)
1129 edihcnstr=enetb(19,t,iparm)
1131 ecationcation=enetb(41,t,iparm)
1132 ecation_prot=enetb(42,t,iparm)
1133 evdwpp = enetb(26,t,iparm)
1134 eespp = enetb(27,t,iparm)
1135 evdwpsb = enetb(28,t,iparm)
1136 eelpsb = enetb(29,t,iparm)
1137 evdwsb = enetb(30,t,iparm)
1138 eelsb = enetb(31,t,iparm)
1139 estr_nucl = enetb(32,t,iparm)
1140 ebe_nucl = enetb(33,t,iparm)
1141 esbloc = enetb(34,t,iparm)
1142 etors_nucl = enetb(35,t,iparm)
1143 etors_d_nucl = enetb(36,t,iparm)
1144 ecorr_nucl = enetb(37,t,iparm)
1145 ecorr3_nucl = enetb(38,t,iparm)
1146 epeppho= enetb(49,t,iparm)
1147 escpho= enetb(48,t,iparm)
1148 epepbase= enetb(47,t,iparm)
1149 escbase= enetb(46,t,iparm)
1150 ecation_nucl=enetb(50,t,iparm)
1151 elipbond=enetb(52,t,iparm)
1152 elipang=enetb(53,t,iparm)
1153 eliplj=enetb(54,t,iparm)
1154 elipelec=enetb(55,t,iparm)
1155 ecat_prottran=enetb(56,t,iparm)
1156 ecation_protang=enetb(57,t,iparm)
1157 eliptran=enetb(22,t,iparm)
1158 ehomology_constr=enetb(51,t,iparm)
1161 betaT=startGridT+k*delta_T
1163 !write(iout,*)"kkkkkkkk",betaT,startGridT,k,delta_T
1165 !d ft=2*T0/(T0+betaT)
1166 if (rescale_modeW.eq.1) then
1174 denom=kfacl-1.0d0+quotl
1176 ftprim(l)=-l*ft(l)*quotl1/(T0*denom)
1177 ftbis(l)=l*kfacl*quotl1* &
1178 (2*l*quotl-(l-1)*denom)/(quot*t0*t0*denom**3)
1181 ft(6)=(320.0d0+80.0d0*dtanh((betaT-320.0d0)/80.0d0))/ &
1183 ftprim(6)=1.0d0/(320.0d0*dcosh((betaT-320.0d0)/80.0d0)**2)
1184 ftbis(6)=-2.0d0*dtanh((betaT-320.0d0)/80.0d0) &
1185 /(320.0d0*80.0d0*dcosh((betaT-320.0d0)/80.0d0)**3)
1186 #elif defined(FUNCT)
1195 else if (rescale_modeW.eq.2) then
1203 logfac=1.0d0/dlog(eplus+eminus)
1204 tanhT=(eplus-eminus)/(eplus+eminus)
1205 fT(l)=1.12692801104297249644d0*logfac
1206 ftprim(l)=-l*quotl1*ft(l)*tanhT*logfac/T0
1207 ftbis(l)=(l-1)*ftprim(l)/(quot*T0)- &
1208 2*l*quotl1/T0*logfac* &
1209 (2*l*quotl1*ft(l)/(T0*(eplus+eminus)**2) &
1213 ft(6)=(320.0d0+80.0d0*dtanh((betaT-320.0d0)/80.0d0))/ &
1215 ftprim(6)=1.0d0/(320.0d0*dcosh((betaT-320.0d0)/80.0d0)**2)
1216 ftbis(6)=-2.0d0*dtanh((betaT-320.0d0)/80.0d0) &
1217 /(320.0d0*80.0d0*dcosh((betaT-320.0d0)/80.0d0)**3)
1218 #elif defined(FUNCT)
1227 else if (rescale_modeW.eq.0) then
1233 write (iout,*) "Error in WHAM_CALC: wrong RESCALE_MODE",&
1238 ! write (iout,*) "ftprim",ftprim
1239 ! write (iout,*) "ftbis",ftbis
1240 betaT=1.0d0/(1.987D-3*betaT)
1242 etot=wsc*(evdw+ft(6)*evdw_t)+wscp*evdw2+ft(1)*welec*ees &
1244 +wang*ebe+ft(1)*wtor*etors+wscloc*escloc &
1245 +wstrain*ehpb+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5 &
1246 +ft(5)*wcorr6*ecorr6+ft(3)*wturn4*eello_turn4 &
1247 +ft(2)*wturn3*eello_turn3 &
1248 +ft(5)*wturn6*eello_turn6+ft(2)*wel_loc*eel_loc &
1249 +edihcnstr+ft(2)*wtor_d*etors_d+ft(1)*wsccor*esccor &
1250 +wbond*estr+wcatprot*ecation_prot+wcatcat*ecationcation &
1251 +wbond_nucl*estr_nucl+wang_nucl*ebe_nucl&
1252 +wvdwpp_nucl*evdwpp+welpp*eespp+wvdwpsb*evdwpsb+welpsb*eelpsb&
1253 +wvdwsb*evdwsb+welsb*eelsb+wsbloc*esbloc+wtor_nucl*etors_nucl&
1254 *ft(1)+wtor_d_nucl*ft(2)*etors_d_nucl+wcorr_nucl*ecorr_nucl &
1255 +wcorr3_nucl*ecorr3_nucl&
1257 +ft(1)*wpepbase*epepbase+wscpho*escpho+wpeppho*epeppho &
1258 + wcatnucl*ecation_nucl&
1259 +elipbond+elipang+eliplj+elipelec+wcat_tran*ecat_prottran+ecation_protang&
1264 eprim=ftprim(6)*evdw_t+ftprim(1)*welec*ees &
1265 +ftprim(1)*wtor*etors+ &
1266 ftprim(3)*wcorr*ecorr+ftprim(4)*wcorr5*ecorr5+ &
1267 ftprim(5)*wcorr6*ecorr6+ftprim(3)*wturn4*eello_turn4+ &
1268 ftprim(2)*wturn3*eello_turn3+ftprim(5)*wturn6*eello_turn6+ &
1269 ftprim(2)*wel_loc*eel_loc+ftprim(2)*wtor_d*etors_d+ &
1270 ftprim(1)*wsccor*esccor +ftprim(1)*wtor_nucl*etors_nucl&
1271 +wtor_d_nucl*ftprim(2)*etors_d_nucl+ftprim(1)*wpepbase*epepbase
1272 ebis=ftbis(1)*welec*ees+ftbis(1)*wtor*etors+ &
1273 ftbis(3)*wcorr*ecorr+ftbis(4)*wcorr5*ecorr5+ &
1274 ftbis(5)*wcorr6*ecorr6+ftbis(3)*wturn4*eello_turn4+ &
1275 ftbis(2)*wturn3*eello_turn3+ftbis(5)*wturn6*eello_turn6+ &
1276 ftbis(2)*wel_loc*eel_loc+ftbis(2)*wtor_d*etors_d+ &
1277 ftbis(1)*wsccor*esccor+ftbis(1)*wtor_nucl*etors_nucl&
1278 +wtor_d_nucl*ftbis(2)*etors_d_nucl+ftbis(1)*wpepbase*epepbase
1280 etot=wsc*(evdw+ft(6)*evdw_t)+wscp*evdw2 &
1281 +ft(1)*welec*(ees+evdw1) &
1282 +wang*ebe+ft(1)*wtor*etors+wscloc*escloc &
1283 +wstrain*ehpb+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5 &
1284 +ft(5)*wcorr6*ecorr6+ft(3)*wturn4*eello_turn4 &
1285 +ft(2)*wturn3*eello_turn3 &
1286 +ft(5)*wturn6*eello_turn6+ft(2)*wel_loc*eel_loc+edihcnstr &
1287 +ft(2)*wtor_d*etors_d+ft(1)*wsccor*esccor &
1288 +wbond*estr+wcatprot*ecation_prot+wcatcat*ecationcation&
1289 +wbond_nucl*estr_nucl+wang_nucl*ebe_nucl&
1290 +wvdwpp_nucl*evdwpp+welpp*eespp+wvdwpsb*evdwpsb+welpsb*eelpsb&
1291 +wvdwsb*evdwsb+welsb*eelsb+wsbloc*esbloc+wtor_nucl*etors_nucl&
1292 *ft(1)+wtor_d_nucl*ft(2)*etors_d_nucl+wcorr_nucl*ecorr_nucl &
1293 +wcorr3_nucl*ecorr3_nucl&
1295 +ft(1)*wpepbase*epepbase+wscpho*escpho+wpeppho*epeppho&
1296 + wcatnucl*ecation_nucl&
1297 +elipbond+elipang+eliplj+elipelec+wcat_tran*ecat_prottran+ecation_protang&
1302 eprim=ftprim(6)*evdw_t+ftprim(1)*welec*(ees+evdw1) &
1303 +ftprim(1)*wtor*etors+ &
1304 ftprim(3)*wcorr*ecorr+ftprim(4)*wcorr5*ecorr5+ &
1305 ftprim(5)*wcorr6*ecorr6+ftprim(3)*wturn4*eello_turn4+ &
1306 ftprim(2)*wturn3*eello_turn3+ftprim(5)*wturn6*eello_turn6+ &
1307 ftprim(2)*wel_loc*eel_loc+ftprim(2)*wtor_d*etors_d+ &
1308 ftprim(1)*wsccor*esccor+ftprim(1)*wtor_nucl*etors_nucl&
1309 +wtor_d_nucl*ftprim(2)*etors_d_nucl+ftprim(1)*wpepbase*epepbase
1310 ebis=ftbis(1)*welec*(ees+evdw1)+ftbis(1)*wtor*etors+ &
1311 ftbis(3)*wcorr*ecorr+ftbis(4)*wcorr5*ecorr5+ &
1312 ftbis(5)*wcorr6*ecorr6+ftbis(3)*wturn4*eello_turn4+ &
1313 ftbis(2)*wturn3*eello_turn3+ftbis(5)*wturn6*eello_turn6+ &
1314 ftbis(2)*wel_loc*eel_loc+ftbis(2)*wtor_d*etors_d+ &
1315 ftbis(1)*wsccor*esccor+ftbis(1)*wtor_nucl*etors_nucl&
1316 +wtor_d_nucl*ftbis(2)*etors_d_nucl*ftbis(1)*wpepbase*epepbase
1319 ! weight=dexp(-betaT*(etot-potEmin)+entfac(t))
1321 ! write (iout,*) "iparm",iparm," t",t," betaT",betaT,&
1322 ! " etot",etot," entfac",entfac(t),&
1323 ! " weight",weight," ebis",ebis
1325 ! etot=etot-temper*eprim
1327 ! sumW_p(k,iparm)=sumW_p(k,iparm)+weight
1328 ! sumE_p(k,iparm)=sumE_p(k,iparm)+etot*weight
1329 ! sumEbis_p(k,iparm)=sumEbis_p(k,iparm)+ebis*weight
1330 ! sumEsq_p(k,iparm)=sumEsq_p(k,iparm)+etot**2*weight
1332 ! sumQ_p(j,k,iparm)=sumQ_p(j,k,iparm)+q(j,t)*weight
1333 ! sumQsq_p(j,k,iparm)=sumQsq_p(j,k,iparm)+q(j,t)**2*weight
1334 ! sumEQ_p(j,k,iparm)=sumEQ_p(j,k,iparm) &
1335 ! +etot*q(j,t)*weight
1338 ! sumW(k,iparm)=sumW(k,iparm)+weight
1339 ! sumE(k,iparm)=sumE(k,iparm)+etot*weight
1340 ! sumEbis(k,iparm)=sumEbis(k,iparm)+ebis*weight
1341 ! sumEsq(k,iparm)=sumEsq(k,iparm)+etot**2*weight
1343 ! sumQ(j,k,iparm)=sumQ(j,k,iparm)+q(j,t)*weight
1344 ! sumQsq(j,k,iparm)=sumQsq(j,k,iparm)+q(j,t)**2*weight
1345 ! sumEQ(j,k,iparm)=sumEQ(j,k,iparm) &
1346 ! +etot*q(j,t)*weight
1350 indE = aint(potE(t,iparm)-aint(potEmin))
1351 if (indE.ge.0 .and. indE.le.maxinde) then
1352 if (indE.gt.upindE_p) upindE_p=indE
1353 histE_p(indE)=histE_p(indE)+dexp(-entfac(t))
1357 potEmin=potEmin_all(ib,iparm)
1358 expfac=dexp(-beta_h(ib,iparm)*(etot-potEmin)+entfac(t))
1359 hfin_p(ind,ib)=hfin_p(ind,ib)+ &
1360 dexp(-beta_h(ib,iparm)*(etot-potEmin)+entfac(t))
1362 indrgy=dint((q(nQ+2,t)-rgymin)/deltrgy)
1363 indrms=dint((q(nQ+1,t)-rmsmin)/deltrms)
1364 hrmsrgy_p(indrgy,indrms,ib)= &
1365 hrmsrgy_p(indrgy,indrms,ib)+expfac
1370 potEmin=potEmin_all(ib,iparm)
1371 expfac=dexp(-beta_h(ib,iparm)*(etot-potEmin)+entfac(t))
1372 hfin(ind,ib)=hfin(ind,ib)+ &
1373 dexp(-beta_h(ib,iparm)*(etot-potEmin)+entfac(t))
1375 indrgy=dint((q(nQ+2,t)-rgymin)/deltrgy)
1376 indrms=dint((q(nQ+1,t)-rmsmin)/deltrms)
1377 hrmsrgy(indrgy,indrms,ib)= &
1378 hrmsrgy(indrgy,indrms,ib)+expfac
1384 ! Thermo and ensemble averages
1387 betaT=startGridT+k*delta_T
1388 ! call temp_scalfac(betaT,ft,ftprim,ftbis,*10)
1389 if (rescale_mode.eq.1) then
1397 denom=kfacl-1.0d0+quotl
1399 ftprim(l)=-l*ft(l)*quotl1/(T0*denom)
1400 ftbis(l)=l*kfacl*quotl1* &
1401 (2*l*quotl-(l-1)*denom)/(quot*t0*t0*denom**3)
1404 ft(6)=(320.0d0+80.0d0*dtanh((betaT-320.0d0)/80.0d0))/ &
1406 ftprim(6)=1.0d0/(320.0d0*dcosh((betaT-320.0d0)/80.0d0)**2)
1407 ftbis(6)=-2.0d0*dtanh((betaT-320.0d0)/80.0d0) &
1408 /(320.0d0*80.0d0*dcosh((betaT-320.0d0)/80.0d0)**3)
1409 #elif defined(FUNCT)
1418 else if (rescale_mode.eq.2) then
1426 logfac=1.0d0/dlog(eplus+eminus)
1427 tanhT=(eplus-eminus)/(eplus+eminus)
1428 fT(l)=1.12692801104297249644d0*logfac
1429 ftprim(l)=-l*quotl1*ft(l)*tanhT*logfac/T0
1430 ftbis(l)=(l-1)*ftprim(l)/(quot*T0)- &
1431 2*l*quotl1/T0*logfac* &
1432 (2*l*quotl1*ft(l)/(T0*(eplus+eminus)**2) &
1436 ft(6)=(320.0d0+80.0d0*dtanh((betaT-320.0d0)/80.0d0))/ &
1438 ftprim(6)=1.0d0/(320.0d0*dcosh((betaT-320.0d0)/80.0d0)**2)
1439 ftbis(6)=-2.0d0*dtanh((betaT-320.0d0)/80.0d0) &
1440 /(320.0d0*80.0d0*dcosh((betaT-320.0d0)/80.0d0)**3)
1441 #elif defined(FUNCT)
1450 else if (rescale_mode.eq.0) then
1456 write (iout,*) "Error in WHAM_CALC: wrong RESCALE_MODE", &
1462 ! write (iout,*) "ftprim",ftprim
1463 ! write (iout,*) "ftbis",ftbis
1464 betaT=1.0d0/(1.987D-3*betaT)
1465 ! 7/10/18 AL Determine the max Botzmann weights for each temerature
1466 ! call sum_ene(1,iparm,ft,etot)
1467 evdw=enetb(1,1,iparm)
1468 evdw_t=enetb(20,1,iparm)
1470 evdw2_14=enetb(18,1,iparm)
1471 evdw2=enetb(2,1,iparm)+evdw2_14
1473 evdw2=enetb(2,1,iparm)
1477 ees=enetb(3,1,iparm)
1478 evdw1=enetb(16,1,iparm)
1480 ees=enetb(3,1,iparm)
1483 ecorr=enetb(4,1,iparm)
1484 ecorr5=enetb(5,1,iparm)
1485 ecorr6=enetb(6,1,iparm)
1486 eel_loc=enetb(7,1,iparm)
1487 eello_turn3=enetb(8,1,iparm)
1488 eello_turn4=enetb(9,1,iparm)
1489 eello_turn6=enetb(10,1,iparm)
1490 ebe=enetb(11,1,iparm)
1491 escloc=enetb(12,1,iparm)
1492 etors=enetb(13,1,iparm)
1493 etors_d=enetb(14,1,iparm)
1494 ehpb=enetb(15,1,iparm)
1495 estr=enetb(17,1,iparm)
1496 esccor=enetb(21,1,iparm)
1497 edihcnstr=enetb(19,1,iparm)
1498 ! eliptran=enetb(22,1,iparm)
1499 ! esaxs=enetb(26,1,iparm)
1500 ecationcation=enetb(41,1,iparm)
1501 ecation_prot=enetb(42,1,iparm)
1502 evdwpp = enetb(26,1,iparm)
1503 eespp = enetb(27,1,iparm)
1504 evdwpsb = enetb(28,1,iparm)
1505 eelpsb = enetb(29,1,iparm)
1506 evdwsb = enetb(30,1,iparm)
1507 eelsb = enetb(31,1,iparm)
1508 estr_nucl = enetb(32,1,iparm)
1509 ebe_nucl = enetb(33,1,iparm)
1510 esbloc = enetb(34,1,iparm)
1511 etors_nucl = enetb(35,1,iparm)
1512 etors_d_nucl = enetb(36,1,iparm)
1513 ecorr_nucl = enetb(37,1,iparm)
1514 ecorr3_nucl = enetb(38,1,iparm)
1515 epeppho= enetb(49,1,iparm)
1516 escpho= enetb(48,1,iparm)
1517 epepbase= enetb(47,1,iparm)
1518 escbase= enetb(46,1,iparm)
1519 ecation_nucl= enetb(50,1,iparm)
1520 ehomology_constr=enetb(51,1,iparm)
1523 if (shield_mode.gt.0) then
1524 etot=ft(1)*wsc*(evdw+ft(6)*evdw_t)+ft(1)*wscp*evdw2 &
1526 +ft(1)*wvdwpp*evdw1 &
1527 +wang*ebe+ft(1)*wtor*etors+wscloc*escloc &
1528 ! & +wstrain*ehpb+nss*ebr+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5
1529 +wstrain*ehpb+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5 &
1530 +ft(5)*wcorr6*ecorr6+ft(3)*wturn4*eello_turn4 &
1531 +ft(2)*wturn3*eello_turn3 &
1532 +ft(5)*wturn6*eello_turn6+ft(2)*wel_loc*eel_loc &
1533 +edihcnstr+ft(2)*wtor_d*etors_d+ft(1)*wsccor*esccor &
1534 +wbond*estr+wcatprot*ecation_prot+wcatcat*ecationcation &
1535 +wbond_nucl*estr_nucl+wang_nucl*ebe_nucl&
1536 +wvdwpp_nucl*evdwpp+welpp*eespp+wvdwpsb*evdwpsb+welpsb*eelpsb&
1537 +wvdwsb*evdwsb+welsb*eelsb+wsbloc*esbloc+wtor_nucl*etors_nucl&
1538 *ft(1)+wtor_d_nucl*ft(2)*etors_d_nucl+wcorr_nucl*ecorr_nucl &
1539 +wcorr3_nucl*ecorr3_nucl&
1541 +ft(1)*wpepbase*epepbase+wscpho*escpho+wpeppho*epeppho &
1542 + wcatnucl*ecation_nucl
1546 etot=wsc*(evdw+ft(6)*evdw_t)+wscp*evdw2+ft(1)*welec*ees &
1548 +wang*ebe+ft(1)*wtor*etors+wscloc*escloc &
1549 ! & +wstrain*ehpb+nss*ebr+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5
1550 +wstrain*ehpb+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5 &
1551 +ft(5)*wcorr6*ecorr6+ft(3)*wturn4*eello_turn4 &
1552 +ft(2)*wturn3*eello_turn3 &
1553 +ft(5)*wturn6*eello_turn6+ft(2)*wel_loc*eel_loc &
1554 +edihcnstr+ft(2)*wtor_d*etors_d+ft(1)*wsccor*esccor &
1555 +wbond*estr+wcatprot*ecation_prot+wcatcat*ecationcation &
1556 +wbond_nucl*estr_nucl+wang_nucl*ebe_nucl&
1557 +wvdwpp_nucl*evdwpp+welpp*eespp+wvdwpsb*evdwpsb+welpsb*eelpsb&
1558 +wvdwsb*evdwsb+welsb*eelsb+wsbloc*esbloc+wtor_nucl*etors_nucl&
1559 *ft(1)+wtor_d_nucl*ft(2)*etors_d_nucl+wcorr_nucl*ecorr_nucl &
1560 +wcorr3_nucl*ecorr3_nucl&
1562 +ft(1)*wpepbase*epepbase+wscpho*escpho+wpeppho*epeppho &
1563 + wcatnucl*ecation_nucl
1567 if (shield_mode.gt.0) then
1568 etot=ft(1)*wsc*(evdw+ft(6)*evdw_t)+ft(1)*wscp*evdw2 &
1569 +ft(1)*welec*(ees+evdw1) &
1570 +wang*ebe+ft(1)*wtor*etors+wscloc*escloc &
1571 ! & +wstrain*ehpb+nss*ebr+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5
1572 +wstrain*ehpb+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5 &
1573 +ft(5)*wcorr6*ecorr6+ft(3)*wturn4*eello_turn4 &
1574 +ft(2)*wturn3*eello_turn3 &
1575 +ft(5)*wturn6*eello_turn6+ft(2)*wel_loc*eel_loc+edihcnstr &
1576 +ft(2)*wtor_d*etors_d+ft(1)*wsccor*esccor &
1577 +wbond*estr+wcatprot*ecation_prot+wcatcat*ecationcation &
1578 +wbond_nucl*estr_nucl+wang_nucl*ebe_nucl&
1579 +wvdwpp_nucl*evdwpp+welpp*eespp+wvdwpsb*evdwpsb+welpsb*eelpsb&
1580 +wvdwsb*evdwsb+welsb*eelsb+wsbloc*esbloc+wtor_nucl*etors_nucl&
1581 *ft(1)+wtor_d_nucl*ft(2)*etors_d_nucl+wcorr_nucl*ecorr_nucl &
1582 +wcorr3_nucl*ecorr3_nucl&
1584 +ft(1)*wpepbase*epepbase+wscpho*escpho+wpeppho*epeppho &
1585 + wcatnucl*ecation_nucl
1588 etot=wsc*(evdw+ft(6)*evdw_t)+wscp*evdw2 &
1589 +ft(1)*welec*(ees+evdw1) &
1590 +wang*ebe+ft(1)*wtor*etors+wscloc*escloc &
1591 ! & +wstrain*ehpb+nss*ebr+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5
1592 +wstrain*ehpb+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5 &
1593 +ft(5)*wcorr6*ecorr6+ft(3)*wturn4*eello_turn4 &
1594 +ft(2)*wturn3*eello_turn3 &
1595 +ft(5)*wturn6*eello_turn6+ft(2)*wel_loc*eel_loc+edihcnstr &
1596 +ft(2)*wtor_d*etors_d+ft(1)*wsccor*esccor &
1597 +wbond*estr+wcatprot*ecation_prot+wcatcat*ecationcation &
1598 +wbond_nucl*estr_nucl+wang_nucl*ebe_nucl&
1599 +wvdwpp_nucl*evdwpp+welpp*eespp+wvdwpsb*evdwpsb+welpsb*eelpsb&
1600 +wvdwsb*evdwsb+welsb*eelsb+wsbloc*esbloc+wtor_nucl*etors_nucl&
1601 *ft(1)+wtor_d_nucl*ft(2)*etors_d_nucl+wcorr_nucl*ecorr_nucl &
1602 +wcorr3_nucl*ecorr3_nucl&
1604 +ft(1)*wpepbase*epepbase+wscpho*escpho+wpeppho*epeppho &
1605 + wcatnucl*ecation_nucl
1610 weimax(k,iparm)=-betaT*(etot-potEmin)+entfac(1)
1611 ! write (iout,*) "k",k," t",1," weight",weimax(k,iparm)
1617 ! call sum_ene(t,iparm,ft,etot)
1618 evdw=enetb(1,t,iparm)
1619 evdw_t=enetb(20,t,iparm)
1621 evdw2_14=enetb(18,t,iparm)
1622 evdw2=enetb(2,t,iparm)+evdw2_14
1624 evdw2=enetb(2,t,iparm)
1628 ees=enetb(3,t,iparm)
1629 evdw1=enetb(16,t,iparm)
1631 ees=enetb(3,t,iparm)
1634 ecorr=enetb(4,t,iparm)
1635 ecorr5=enetb(5,t,iparm)
1636 ecorr6=enetb(6,t,iparm)
1637 eel_loc=enetb(7,t,iparm)
1638 eello_turn3=enetb(8,t,iparm)
1639 eello_turn4=enetb(9,t,iparm)
1640 eello_turn6=enetb(10,t,iparm)
1641 ebe=enetb(11,t,iparm)
1642 escloc=enetb(12,t,iparm)
1643 etors=enetb(13,t,iparm)
1644 etors_d=enetb(14,t,iparm)
1645 ehpb=enetb(15,t,iparm)
1646 estr=enetb(17,t,iparm)
1647 esccor=enetb(21,t,iparm)
1648 edihcnstr=enetb(19,t,iparm)
1649 ! eliptran=enetb(22,t,iparm)
1650 ! esaxs=enetb(26,t,iparm)
1651 ecationcation=enetb(41,t,iparm)
1652 ecation_prot=enetb(42,t,iparm)
1653 evdwpp = enetb(26,t,iparm)
1654 eespp = enetb(27,t,iparm)
1655 evdwpsb = enetb(28,t,iparm)
1656 eelpsb = enetb(29,t,iparm)
1657 evdwsb = enetb(30,t,iparm)
1658 eelsb = enetb(31,t,iparm)
1659 estr_nucl = enetb(32,t,iparm)
1660 ebe_nucl = enetb(33,t,iparm)
1661 esbloc = enetb(34,t,iparm)
1662 etors_nucl = enetb(35,t,iparm)
1663 etors_d_nucl = enetb(36,t,iparm)
1664 ecorr_nucl = enetb(37,t,iparm)
1665 ecorr3_nucl = enetb(38,t,iparm)
1666 epeppho= enetb(49,t,iparm)
1667 escpho= enetb(48,t,iparm)
1668 epepbase= enetb(47,t,iparm)
1669 escbase= enetb(46,t,iparm)
1670 ecation_nucl=enetb(50,t,iparm)
1671 ehomology_constr=enetb(51,t,iparm)
1674 if (shield_mode.gt.0) then
1675 etot=ft(1)*wsc*(evdw+ft(6)*evdw_t)+ft(1)*wscp*evdw2 &
1677 +ft(1)*wvdwpp*evdw1 &
1678 +wang*ebe+ft(1)*wtor*etors+wscloc*escloc &
1679 ! & +wstrain*ehpb+nss*ebr+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5
1680 +wstrain*ehpb+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5 &
1681 +ft(5)*wcorr6*ecorr6+ft(3)*wturn4*eello_turn4 &
1682 +ft(2)*wturn3*eello_turn3 &
1683 +ft(5)*wturn6*eello_turn6+ft(2)*wel_loc*eel_loc &
1684 +edihcnstr+ft(2)*wtor_d*etors_d+ft(1)*wsccor*esccor &
1685 +wbond*estr+wcatprot*ecation_prot+wcatcat*ecationcation &
1686 +wbond_nucl*estr_nucl+wang_nucl*ebe_nucl&
1687 +wvdwpp_nucl*evdwpp+welpp*eespp+wvdwpsb*evdwpsb+welpsb*eelpsb&
1688 +wvdwsb*evdwsb+welsb*eelsb+wsbloc*esbloc+wtor_nucl*etors_nucl&
1689 *ft(1)+wtor_d_nucl*ft(2)*etors_d_nucl+wcorr_nucl*ecorr_nucl &
1690 +wcorr3_nucl*ecorr3_nucl&
1692 +ft(1)*wpepbase*epepbase+wscpho*escpho+wpeppho*epeppho &
1693 + wcatnucl*ecation_nucl
1696 etot=wsc*(evdw+ft(6)*evdw_t)+wscp*evdw2+ft(1)*welec*ees &
1698 +wang*ebe+ft(1)*wtor*etors+wscloc*escloc &
1699 ! & +wstrain*ehpb+nss*ebr+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5
1700 +wstrain*ehpb+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5 &
1701 +ft(5)*wcorr6*ecorr6+ft(3)*wturn4*eello_turn4 &
1702 +ft(2)*wturn3*eello_turn3 &
1703 +ft(5)*wturn6*eello_turn6+ft(2)*wel_loc*eel_loc &
1704 +edihcnstr+ft(2)*wtor_d*etors_d+ft(1)*wsccor*esccor &
1705 +wbond*estr+wcatprot*ecation_prot+wcatcat*ecationcation &
1706 +wbond_nucl*estr_nucl+wang_nucl*ebe_nucl&
1707 +wvdwpp_nucl*evdwpp+welpp*eespp+wvdwpsb*evdwpsb+welpsb*eelpsb&
1708 +wvdwsb*evdwsb+welsb*eelsb+wsbloc*esbloc+wtor_nucl*etors_nucl&
1709 *ft(1)+wtor_d_nucl*ft(2)*etors_d_nucl+wcorr_nucl*ecorr_nucl &
1710 +wcorr3_nucl*ecorr3_nucl&
1712 +ft(1)*wpepbase*epepbase+wscpho*escpho+wpeppho*epeppho &
1713 + wcatnucl*ecation_nucl
1716 if (shield_mode.gt.0) then
1717 etot=ft(1)*wsc*(evdw+ft(6)*evdw_t)+ft(1)*wscp*evdw2 &
1718 +ft(1)*welec*(ees+evdw1) &
1719 +wang*ebe+ft(1)*wtor*etors+wscloc*escloc &
1720 ! & +wstrain*ehpb+nss*ebr+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5
1721 +wstrain*ehpb+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5 &
1722 +ft(5)*wcorr6*ecorr6+ft(3)*wturn4*eello_turn4 &
1723 +ft(2)*wturn3*eello_turn3 &
1724 +ft(5)*wturn6*eello_turn6+ft(2)*wel_loc*eel_loc+edihcnstr &
1725 +ft(2)*wtor_d*etors_d+ft(1)*wsccor*esccor &
1726 +wbond*estr+wcatprot*ecation_prot+wcatcat*ecationcation &
1727 +wbond_nucl*estr_nucl+wang_nucl*ebe_nucl&
1728 +wvdwpp_nucl*evdwpp+welpp*eespp+wvdwpsb*evdwpsb+welpsb*eelpsb&
1729 +wvdwsb*evdwsb+welsb*eelsb+wsbloc*esbloc+wtor_nucl*etors_nucl&
1730 *ft(1)+wtor_d_nucl*ft(2)*etors_d_nucl+wcorr_nucl*ecorr_nucl &
1731 +wcorr3_nucl*ecorr3_nucl&
1733 +ft(1)*wpepbase*epepbase+wscpho*escpho+wpeppho*epeppho &
1734 + wcatnucl*ecation_nucl
1736 etot=wsc*(evdw+ft(6)*evdw_t)+wscp*evdw2 &
1737 +ft(1)*welec*(ees+evdw1) &
1738 +wang*ebe+ft(1)*wtor*etors+wscloc*escloc &
1739 ! & +wstrain*ehpb+nss*ebr+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5
1740 +wstrain*ehpb+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5 &
1741 +ft(5)*wcorr6*ecorr6+ft(3)*wturn4*eello_turn4 &
1742 +ft(2)*wturn3*eello_turn3 &
1743 +ft(5)*wturn6*eello_turn6+ft(2)*wel_loc*eel_loc+edihcnstr &
1744 +ft(2)*wtor_d*etors_d+ft(1)*wsccor*esccor &
1745 +wbond*estr+wcatprot*ecation_prot+wcatcat*ecationcation &
1746 +wbond_nucl*estr_nucl+wang_nucl*ebe_nucl&
1747 +wvdwpp_nucl*evdwpp+welpp*eespp+wvdwpsb*evdwpsb+welpsb*eelpsb&
1748 +wvdwsb*evdwsb+welsb*eelsb+wsbloc*esbloc+wtor_nucl*etors_nucl&
1749 *ft(1)+wtor_d_nucl*ft(2)*etors_d_nucl+wcorr_nucl*ecorr_nucl &
1750 +wcorr3_nucl*ecorr3_nucl&
1752 +ft(1)*wpepbase*epepbase+wscpho*escpho+wpeppho*epeppho &
1753 + wcatnucl*ecation_nucl
1757 weight=-betaT*(etot-potEmin)+entfac(t)
1758 ! write (iout,*) "k",k," t",t," weight",weight
1759 if (weight.gt.weimax(k,iparm)) weimax(k,iparm)=weight
1764 write (iout,*) "weimax before REDUCE"
1765 write (iout,*) (weimax(k,iparm),k=0,ngridt)
1768 weimax_(k)=weimax(k,iparm)
1770 call MPI_Allreduce(weimax_(0),weimax(0,iparm),nGridT+1, &
1771 MPI_DOUBLE_PRECISION,MPI_MAX,WHAM_COMM,IERROR)
1773 write (iout,*) "weimax"
1774 write (iout,*) (weimax(k,iparm),k=0,ngridt)
1777 weimax_(k)=weimax(k,iparm)
1779 call MPI_Allreduce(weimax_(0),weimax(0,iparm),nGridT+1, &
1780 MPI_DOUBLE_PRECISION,MPI_MAX,WHAM_COMM,IERROR)
1782 write (iout,*) "weimax"
1783 write (iout,*) (weimax(k,iparm),k=0,ngridt)
1786 temper=startGridT+k*delta_T
1787 betaT=1.0d0/(1.987D-3*temper)
1788 ! call temp_scalfac(temper,ft,ftprim,ftbis,*10)
1789 if (rescale_mode.eq.1) then
1797 denom=kfacl-1.0d0+quotl
1799 ftprim(l)=-l*ft(l)*quotl1/(T0*denom)
1800 ftbis(l)=l*kfacl*quotl1*&
1801 (2*l*quotl-(l-1)*denom)/(quot*t0*t0*denom**3)
1804 ft(6)=(320.0d0+80.0d0*dtanh((temper-320.0d0)/80.0d0))/ &
1806 ftprim(6)=1.0d0/(320.0d0*dcosh((temper-320.0d0)/80.0d0)**2)
1807 ftbis(6)=-2.0d0*dtanh((temper-320.0d0)/80.0d0) &
1808 /(320.0d0*80.0d0*dcosh((temper-320.0d0)/80.0d0)**3)
1809 #elif defined(FUNCT)
1818 else if (rescale_mode.eq.2) then
1826 logfac=1.0d0/dlog(eplus+eminus)
1827 tanhT=(eplus-eminus)/(eplus+eminus)
1828 fT(l)=1.12692801104297249644d0*logfac
1829 ftprim(l)=-l*quotl1*ft(l)*tanhT*logfac/T0
1830 ftbis(l)=(l-1)*ftprim(l)/(quot*T0)- &
1831 2*l*quotl1/T0*logfac* &
1832 (2*l*quotl1*ft(l)/(T0*(eplus+eminus)**2) &
1836 ft(6)=(320.0d0+80.0d0*dtanh((temper-320.0d0)/80.0d0))/ &
1838 ftprim(6)=1.0d0/(320.0d0*dcosh((temper-320.0d0)/80.0d0)**2)
1839 ftbis(6)=-2.0d0*dtanh((temper-320.0d0)/80.0d0) &
1840 /(320.0d0*80.0d0*dcosh((temper-320.0d0)/80.0d0)**3)
1841 #elif defined(FUNCT)
1850 else if (rescale_mode.eq.0) then
1856 write (iout,*) "Error in WHAM_CALC: wrong RESCALE_MODE", &
1870 hfin_ent_p(ind)=hfin_ent_p(ind)+dexp(entfac(t))
1872 hfin_ent(ind)=hfin_ent(ind)+dexp(entfac(t))
1874 ! write (iout,'(2i5,20f8.2)') t,t,(enetb(k,t,iparm),k=1,18)
1875 ! call restore_parm(iparm)
1876 ! call sum_ene_deriv(t,iparm,ft,ftprim,ftbis,etot,eprim,ebis)
1877 evdw=enetb(1,t,iparm)
1878 evdw_t=enetb(20,t,iparm)
1880 evdw2_14=enetb(18,t,iparm)
1881 evdw2=enetb(2,t,iparm)+evdw2_14
1883 evdw2=enetb(2,t,iparm)
1887 ees=enetb(3,t,iparm)
1888 evdw1=enetb(16,t,iparm)
1890 ees=enetb(3,t,iparm)
1893 ecorr=enetb(4,t,iparm)
1894 ecorr5=enetb(5,t,iparm)
1895 ecorr6=enetb(6,t,iparm)
1896 eel_loc=enetb(7,t,iparm)
1897 eello_turn3=enetb(8,t,iparm)
1898 eello_turn4=enetb(9,t,iparm)
1899 eello_turn6=enetb(10,t,iparm)
1900 ebe=enetb(11,t,iparm)
1901 escloc=enetb(12,t,iparm)
1902 etors=enetb(13,t,iparm)
1903 etors_d=enetb(14,t,iparm)
1904 ehpb=enetb(15,t,iparm)
1905 estr=enetb(17,t,iparm)
1906 esccor=enetb(21,t,iparm)
1907 edihcnstr=enetb(19,t,iparm)
1908 ! eliptran=enetb(22,t,iparm)
1909 ! esaxs=enetb(26,t,iparm)
1910 ! ehomology_constr=enetb(27,t,iparm)
1911 ! edfadis=enetb(28,t,iparm)
1912 ! edfator=enetb(29,t,iparm)
1913 ! edfanei=enetb(30,t,iparm)
1914 ! edfabet=enetb(31,t,iparm)
1915 ecationcation=enetb(41,i,iparm)
1916 ecation_prot=enetb(42,i,iparm)
1917 evdwpp = enetb(26,i,iparm)
1918 eespp = enetb(27,i,iparm)
1919 evdwpsb = enetb(28,i,iparm)
1920 eelpsb = enetb(29,i,iparm)
1921 evdwsb = enetb(30,i,iparm)
1922 eelsb = enetb(31,i,iparm)
1923 estr_nucl = enetb(32,i,iparm)
1924 ebe_nucl = enetb(33,i,iparm)
1925 esbloc = enetb(34,i,iparm)
1926 etors_nucl = enetb(35,i,iparm)
1927 etors_d_nucl = enetb(36,i,iparm)
1928 ecorr_nucl = enetb(37,i,iparm)
1929 ecorr3_nucl = enetb(38,i,iparm)
1930 epeppho= enetb(49,i,iparm)
1931 escpho= enetb(48,i,iparm)
1932 epepbase= enetb(47,i,iparm)
1933 escbase= enetb(46,i,iparm)
1934 ecation_nucl= enetb(50,i,iparm)
1935 ehomology_constr=enetb(51,i,iparm)
1938 if (shield_mode.gt.0) then
1939 etot=ft(1)*wsc*(evdw+ft(6)*evdw_t)+ft(1)*wscp*evdw2 &
1941 +ft(1)*wvdwpp*evdw1 &
1942 +wang*ebe+ft(1)*wtor*etors+wscloc*escloc &
1943 ! & +wstrain*ehpb+nss*ebr+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5
1944 +wstrain*ehpb+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5 &
1945 +ft(5)*wcorr6*ecorr6+ft(3)*wturn4*eello_turn4 &
1946 +ft(2)*wturn3*eello_turn3 &
1947 +ft(5)*wturn6*eello_turn6+ft(2)*wel_loc*eel_loc &
1948 +edihcnstr+ft(2)*wtor_d*etors_d+ft(1)*wsccor*esccor &
1949 +wbond*estr+wcatprot*ecation_prot+wcatcat*ecationcation &
1950 +wbond_nucl*estr_nucl+wang_nucl*ebe_nucl&
1951 +wvdwpp_nucl*evdwpp+welpp*eespp+wvdwpsb*evdwpsb+welpsb*eelpsb&
1952 +wvdwsb*evdwsb+welsb*eelsb+wsbloc*esbloc+wtor_nucl*etors_nucl&
1953 *ft(1)+wtor_d_nucl*ft(2)*etors_d_nucl+wcorr_nucl*ecorr_nucl &
1954 +wcorr3_nucl*ecorr3_nucl&
1956 +ft(1)*wpepbase*epepbase+wscpho*escpho+wpeppho*epeppho &
1957 + wcatnucl*ecation_nucl
1959 ! eprim=ftprim(1)*(ft(6)*evdw_t+evdw)
1960 !! & +ftprim(6)*evdw_t
1961 ! & +ftprim(1)*wscp*evdw2
1962 ! & +ftprim(1)*welec*ees
1963 ! & +ftprim(1)*wvdwpp*evdw1
1964 ! & +ftprim(1)*wtor*etors+
1965 ! & ftprim(3)*wcorr*ecorr+ftprim(4)*wcorr5*ecorr5+
1966 ! & ftprim(5)*wcorr6*ecorr6+ftprim(3)*wturn4*eello_turn4+
1967 ! & ftprim(2)*wturn3*eello_turn3+ftprim(5)*wturn6*eello_turn6+
1968 ! & ftprim(2)*wel_loc*eel_loc+ftprim(2)*wtor_d*etors_d+
1969 ! & ftprim(1)*wsccor*esccor
1971 ! ebis=ftbis(1)*wsc*(evdw+ft(6)*evdw_t)
1972 ! & +ftbis(1)*wscp*evdw2+
1973 ! & ftbis(1)*welec*ees
1974 ! & +ftbis(1)*wvdwpp*evdw
1975 ! & +ftbis(1)*wtor*etors+
1976 ! & ftbis(3)*wcorr*ecorr+ftbis(4)*wcorr5*ecorr5+
1977 ! & ftbis(5)*wcorr6*ecorr6+ftbis(3)*wturn4*eello_turn4+
1978 ! & ftbis(2)*wturn3*eello_turn3+ftbis(5)*wturn6*eello_turn6+
1979 ! & ftbis(2)*wel_loc*eel_loc+ftbis(2)*wtor_d*etors_d+
1980 ! & ftbis(1)*wsccor*esccor
1982 eprim=ftprim(6)*evdw_t+ftprim(1)*welec*(ees+evdw1) &
1983 +ftprim(1)*wtor*etors+ &
1984 ftprim(3)*wcorr*ecorr+ftprim(4)*wcorr5*ecorr5+ &
1985 ftprim(5)*wcorr6*ecorr6+ftprim(3)*wturn4*eello_turn4+ &
1986 ftprim(2)*wturn3*eello_turn3+ftprim(5)*wturn6*eello_turn6+ &
1987 ftprim(2)*wel_loc*eel_loc+ftprim(2)*wtor_d*etors_d+ &
1988 ftprim(1)*wsccor*esccor+ftprim(1)*wtor_nucl*etors_nucl&
1989 +wtor_d_nucl*ftprim(2)*etors_d_nucl+ftprim(1)*wpepbase*epepbase
1990 ebis=ftbis(1)*welec*(ees+evdw1)+ftbis(1)*wtor*etors+ &
1991 ftbis(3)*wcorr*ecorr+ftbis(4)*wcorr5*ecorr5+ &
1992 ftbis(5)*wcorr6*ecorr6+ftbis(3)*wturn4*eello_turn4+ &
1993 ftbis(2)*wturn3*eello_turn3+ftbis(5)*wturn6*eello_turn6+ &
1994 ftbis(2)*wel_loc*eel_loc+ftbis(2)*wtor_d*etors_d+ &
1995 ftbis(1)*wsccor*esccor+ftbis(1)*wtor_nucl*etors_nucl&
1996 +wtor_d_nucl*ftbis(2)*etors_d_nucl*ftbis(1)*wpepbase*epepbase
2001 etot=wsc*(evdw+ft(6)*evdw_t)+wscp*evdw2+ft(1)*welec*ees &
2003 +wang*ebe+ft(1)*wtor*etors+wscloc*escloc &
2004 !c & +wstrain*ehpb+nss*ebr+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5
2005 +wstrain*ehpb+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5 &
2006 +ft(5)*wcorr6*ecorr6+ft(3)*wturn4*eello_turn4 &
2007 +ft(2)*wturn3*eello_turn3 &
2008 +ft(5)*wturn6*eello_turn6+ft(2)*wel_loc*eel_loc &
2009 +edihcnstr+ft(2)*wtor_d*etors_d+ft(1)*wsccor*esccor &
2010 +wbond*estr+wcatprot*ecation_prot+wcatcat*ecationcation&
2011 +wbond_nucl*estr_nucl+wang_nucl*ebe_nucl&
2012 +wvdwpp_nucl*evdwpp+welpp*eespp+wvdwpsb*evdwpsb+welpsb*eelpsb&
2013 +wvdwsb*evdwsb+welsb*eelsb+wsbloc*esbloc+wtor_nucl*etors_nucl&
2014 *ft(1)+wtor_d_nucl*ft(2)*etors_d_nucl+wcorr_nucl*ecorr_nucl &
2015 +wcorr3_nucl*ecorr3_nucl&
2017 +ft(1)*wpepbase*epepbase+wscpho*escpho+wpeppho*epeppho&
2018 + wcatnucl*ecation_nucl
2021 eprim=ftprim(6)*evdw_t+ftprim(1)*welec*(ees+evdw1) &
2022 +ftprim(1)*wtor*etors+ &
2023 ftprim(3)*wcorr*ecorr+ftprim(4)*wcorr5*ecorr5+ &
2024 ftprim(5)*wcorr6*ecorr6+ftprim(3)*wturn4*eello_turn4+ &
2025 ftprim(2)*wturn3*eello_turn3+ftprim(5)*wturn6*eello_turn6+ &
2026 ftprim(2)*wel_loc*eel_loc+ftprim(2)*wtor_d*etors_d+ &
2027 ftprim(1)*wsccor*esccor+ftprim(1)*wtor_nucl*etors_nucl&
2028 +wtor_d_nucl*ftprim(2)*etors_d_nucl+ftprim(1)*wpepbase*epepbase
2029 ebis=ftbis(1)*welec*(ees+evdw1)+ftbis(1)*wtor*etors+ &
2030 ftbis(3)*wcorr*ecorr+ftbis(4)*wcorr5*ecorr5+ &
2031 ftbis(5)*wcorr6*ecorr6+ftbis(3)*wturn4*eello_turn4+ &
2032 ftbis(2)*wturn3*eello_turn3+ftbis(5)*wturn6*eello_turn6+ &
2033 ftbis(2)*wel_loc*eel_loc+ftbis(2)*wtor_d*etors_d+ &
2034 ftbis(1)*wsccor*esccor+ftbis(1)*wtor_nucl*etors_nucl&
2035 +wtor_d_nucl*ftbis(2)*etors_d_nucl*ftbis(1)*wpepbase*epepbase
2039 ! eprim=ftprim(6)*evdw_t+ftprim(1)*welec*ees
2040 ! & +ftprim(1)*wtor*etors+
2041 ! & ftprim(3)*wcorr*ecorr+ftprim(4)*wcorr5*ecorr5+
2042 ! & ftprim(5)*wcorr6*ecorr6+ftprim(3)*wturn4*eello_turn4+
2043 ! & ftprim(2)*wturn3*eello_turn3+ftprim(5)*wturn6*eello_turn6+
2044 ! & ftprim(2)*wel_loc*eel_loc+ftprim(2)*wtor_d*etors_d+
2045 ! & ftprim(1)*wsccor*esccor
2046 ! ebis=ftbis(1)*welec*ees+ftbis(1)*wtor*etors+
2047 ! & ftbis(3)*wcorr*ecorr+ftbis(4)*wcorr5*ecorr5+
2048 ! & ftbis(5)*wcorr6*ecorr6+ftbis(3)*wturn4*eello_turn4+
2049 ! & ftbis(2)*wturn3*eello_turn3+ftbis(5)*wturn6*eello_turn6+
2050 ! & ftbis(2)*wel_loc*eel_loc+ftbis(2)*wtor_d*etors_d+
2051 ! & ftbis(1)*wsccor*esccor
2054 if (shield_mode.gt.0) then
2055 etot=ft(1)*wsc*(evdw+ft(6)*evdw_t)+ft(1)*wscp*evdw2 &
2056 +ft(1)*welec*(ees+evdw1) &
2057 +wang*ebe+ft(1)*wtor*etors+wscloc*escloc &
2058 !c & +wstrain*ehpb+nss*ebr+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5
2059 +wstrain*ehpb+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5 &
2060 +ft(5)*wcorr6*ecorr6+ft(3)*wturn4*eello_turn4 &
2061 +ft(2)*wturn3*eello_turn3 &
2062 +ft(5)*wturn6*eello_turn6+ft(2)*wel_loc*eel_loc+edihcnstr &
2063 +ft(2)*wtor_d*etors_d+ft(1)*wsccor*esccor &
2064 +wbond*estr+wcatprot*ecation_prot+wcatcat*ecationcation&
2065 +wbond_nucl*estr_nucl+wang_nucl*ebe_nucl&
2066 +wvdwpp_nucl*evdwpp+welpp*eespp+wvdwpsb*evdwpsb+welpsb*eelpsb&
2067 +wvdwsb*evdwsb+welsb*eelsb+wsbloc*esbloc+wtor_nucl*etors_nucl&
2068 *ft(1)+wtor_d_nucl*ft(2)*etors_d_nucl+wcorr_nucl*ecorr_nucl &
2069 +wcorr3_nucl*ecorr3_nucl&
2071 +ft(1)*wpepbase*epepbase+wscpho*escpho+wpeppho*epeppho&
2072 + wcatnucl*ecation_nucl
2075 eprim=ftprim(6)*evdw_t+ftprim(1)*welec*(ees+evdw1) &
2076 +ftprim(1)*wtor*etors+ &
2077 ftprim(3)*wcorr*ecorr+ftprim(4)*wcorr5*ecorr5+ &
2078 ftprim(5)*wcorr6*ecorr6+ftprim(3)*wturn4*eello_turn4+ &
2079 ftprim(2)*wturn3*eello_turn3+ftprim(5)*wturn6*eello_turn6+ &
2080 ftprim(2)*wel_loc*eel_loc+ftprim(2)*wtor_d*etors_d+ &
2081 ftprim(1)*wsccor*esccor+ftprim(1)*wtor_nucl*etors_nucl&
2082 +wtor_d_nucl*ftprim(2)*etors_d_nucl+ftprim(1)*wpepbase*epepbase
2083 ebis=ftbis(1)*welec*(ees+evdw1)+ftbis(1)*wtor*etors+ &
2084 ftbis(3)*wcorr*ecorr+ftbis(4)*wcorr5*ecorr5+ &
2085 ftbis(5)*wcorr6*ecorr6+ftbis(3)*wturn4*eello_turn4+ &
2086 ftbis(2)*wturn3*eello_turn3+ftbis(5)*wturn6*eello_turn6+ &
2087 ftbis(2)*wel_loc*eel_loc+ftbis(2)*wtor_d*etors_d+ &
2088 ftbis(1)*wsccor*esccor+ftbis(1)*wtor_nucl*etors_nucl&
2089 +wtor_d_nucl*ftbis(2)*etors_d_nucl*ftbis(1)*wpepbase*epepbase
2094 ! eprim=ftprim(1)*(evdw+ft(6)*evdw_t)
2095 ! & +ftprim(1)*welec*(ees+evdw1)
2096 ! & +ftprim(1)*wtor*etors+
2097 ! & ftprim(1)*wscp*evdw2+
2098 ! & ftprim(3)*wcorr*ecorr+ftprim(4)*wcorr5*ecorr5+
2099 ! & ftprim(5)*wcorr6*ecorr6+ftprim(3)*wturn4*eello_turn4+
2100 ! & ftprim(2)*wturn3*eello_turn3+ftprim(5)*wturn6*eello_turn6+
2101 ! & ftprim(2)*wel_loc*eel_loc+ftprim(2)*wtor_d*etors_d+
2102 ! & ftprim(1)*wsccor*esccor
2103 ! ebis= ftbis(1)*(evdw+ft(6)*evdw_t)
2104 ! & +ftbis(1)*wscp*evdw2
2105 ! & +ftbis(1)*welec*(ees+evdw1)+ftbis(1)*wtor*etors+
2106 ! & ftbis(3)*wcorr*ecorr+ftbis(4)*wcorr5*ecorr5+
2107 ! & ftbis(5)*wcorr6*ecorr6+ftbis(3)*wturn4*eello_turn4+
2108 ! & ftbis(2)*wturn3*eello_turn3+ftbis(5)*wturn6*eello_turn6+
2109 ! & ftbis(2)*wel_loc*eel_loc+ftbis(2)*wtor_d*etors_d+
2110 ! & ftprim(1)*wsccor*esccor
2112 etot=wsc*(evdw+ft(6)*evdw_t)+wscp*evdw2 &
2113 +ft(1)*welec*(ees+evdw1) &
2114 +wang*ebe+ft(1)*wtor*etors+wscloc*escloc &
2115 !c & +wstrain*ehpb+nss*ebr+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5
2116 +wstrain*ehpb+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5 &
2117 +ft(5)*wcorr6*ecorr6+ft(3)*wturn4*eello_turn4 &
2118 +ft(2)*wturn3*eello_turn3 &
2119 +ft(5)*wturn6*eello_turn6+ft(2)*wel_loc*eel_loc+edihcnstr &
2120 +ft(2)*wtor_d*etors_d+ft(1)*wsccor*esccor &
2121 +wbond*estr+wcatprot*ecation_prot+wcatcat*ecationcation&
2122 +wbond_nucl*estr_nucl+wang_nucl*ebe_nucl&
2123 +wvdwpp_nucl*evdwpp+welpp*eespp+wvdwpsb*evdwpsb+welpsb*eelpsb&
2124 +wvdwsb*evdwsb+welsb*eelsb+wsbloc*esbloc+wtor_nucl*etors_nucl&
2125 *ft(1)+wtor_d_nucl*ft(2)*etors_d_nucl+wcorr_nucl*ecorr_nucl &
2126 +wcorr3_nucl*ecorr3_nucl&
2128 +ft(1)*wpepbase*epepbase+wscpho*escpho+wpeppho*epeppho&
2129 + wcatnucl*ecation_nucl
2132 eprim=ftprim(6)*evdw_t+ftprim(1)*welec*(ees+evdw1) &
2133 +ftprim(1)*wtor*etors+ &
2134 ftprim(3)*wcorr*ecorr+ftprim(4)*wcorr5*ecorr5+ &
2135 ftprim(5)*wcorr6*ecorr6+ftprim(3)*wturn4*eello_turn4+ &
2136 ftprim(2)*wturn3*eello_turn3+ftprim(5)*wturn6*eello_turn6+ &
2137 ftprim(2)*wel_loc*eel_loc+ftprim(2)*wtor_d*etors_d+ &
2138 ftprim(1)*wsccor*esccor+ftprim(1)*wtor_nucl*etors_nucl&
2139 +wtor_d_nucl*ftprim(2)*etors_d_nucl+ftprim(1)*wpepbase*epepbase
2140 ebis=ftbis(1)*welec*(ees+evdw1)+ftbis(1)*wtor*etors+ &
2141 ftbis(3)*wcorr*ecorr+ftbis(4)*wcorr5*ecorr5+ &
2142 ftbis(5)*wcorr6*ecorr6+ftbis(3)*wturn4*eello_turn4+ &
2143 ftbis(2)*wturn3*eello_turn3+ftbis(5)*wturn6*eello_turn6+ &
2144 ftbis(2)*wel_loc*eel_loc+ftbis(2)*wtor_d*etors_d+ &
2145 ftbis(1)*wsccor*esccor+ftbis(1)*wtor_nucl*etors_nucl&
2146 +wtor_d_nucl*ftbis(2)*etors_d_nucl*ftbis(1)*wpepbase*epepbase
2152 ! eprim=ftprim(6)*evdw_t+ftprim(1)*welec*(ees+evdw1)
2153 ! & +ftprim(1)*wtor*etors+
2154 ! & ftprim(3)*wcorr*ecorr+ftprim(4)*wcorr5*ecorr5+
2155 ! & ftprim(5)*wcorr6*ecorr6+ftprim(3)*wturn4*eello_turn4+
2156 ! & ftprim(2)*wturn3*eello_turn3+ftprim(5)*wturn6*eello_turn6+
2157 ! & ftprim(2)*wel_loc*eel_loc+ftprim(2)*wtor_d*etors_d+
2158 ! & ftprim(1)*wsccor*esccor
2159 ! ebis=ftbis(1)*welec*(ees+evdw1)+ftbis(1)*wtor*etors+
2160 ! & ftbis(3)*wcorr*ecorr+ftbis(4)*wcorr5*ecorr5+
2161 ! & ftbis(5)*wcorr6*ecorr6+ftbis(3)*wturn4*eello_turn4+
2162 ! & ftbis(2)*wturn3*eello_turn3+ftbis(5)*wturn6*eello_turn6+
2163 ! & ftbis(2)*wel_loc*eel_loc+ftbis(2)*wtor_d*etors_d+
2164 ! & ftprim(1)*wsccor*esccor
2168 weight=dexp(-betaT*(etot-potEmin)+entfac(t)-weimax(k,iparm))
2170 write (iout,*) "iparm",iparm," t",t," betaT",betaT, &
2171 " etot",etot," entfac",entfac(t)," boltz", &
2172 -betaT*(etot-potEmin)+entfac(t)," weimax",weimax(k,iparm), &
2173 " weight",weight," ebis",ebis
2175 etot=etot-temper*eprim
2177 sumW_p(k,iparm)=sumW_p(k,iparm)+weight
2178 sumE_p(k,iparm)=sumE_p(k,iparm)+etot*weight
2179 sumEbis_p(k,iparm)=sumEbis_p(k,iparm)+ebis*weight
2180 sumEsq_p(k,iparm)=sumEsq_p(k,iparm)+etot**2*weight
2182 sumQ_p(j,k,iparm)=sumQ_p(j,k,iparm)+q(j,t)*weight
2183 sumQsq_p(j,k,iparm)=sumQsq_p(j,k,iparm)+q(j,t)**2*weight
2184 sumEQ_p(j,k,iparm)=sumEQ_p(j,k,iparm) &
2188 sumW(k,iparm)=sumW(k,iparm)+weight
2189 sumE(k,iparm)=sumE(k,iparm)+etot*weight
2190 sumEbis(k,iparm)=sumEbis(k,iparm)+ebis*weight
2191 sumEsq(k,iparm)=sumEsq(k,iparm)+etot**2*weight
2193 sumQ(j,k,iparm)=sumQ(j,k,iparm)+q(j,t)*weight
2194 sumQsq(j,k,iparm)=sumQsq(j,k,iparm)+q(j,t)**2*weight
2195 sumEQ(j,k,iparm)=sumEQ(j,k,iparm) &
2203 if (histout) call MPI_Reduce(hfin_p(0,ib),hfin(0,ib),nbin,&
2204 MPI_DOUBLE_PRECISION,MPI_SUM,Master,WHAM_COMM,IERROR)
2206 call MPI_Reduce(hrmsrgy_p(0,0,ib),hrmsrgy(0,0,ib),&
2207 (MaxBinRgy+1)*(nbin_rms+1),MPI_DOUBLE_PRECISION,MPI_SUM,Master,&
2211 call MPI_Reduce(upindE_p,upindE,1,&
2212 MPI_INTEGER,MPI_MAX,Master,WHAM_COMM,IERROR)
2213 call MPI_Reduce(histE_p(0),histE(0),maxindE,&
2214 MPI_DOUBLE_PRECISION,MPI_SUM,Master,WHAM_COMM,IERROR)
2216 if (me1.eq.master) then
2220 write (iout,'(6x,$)')
2221 write (iout,'(f20.2,$)') (1.0d0/(1.987D-3*beta_h(ib,iparm)),&
2225 write (iout,'(/a)') 'Final histograms'
2227 if (nslice.eq.1) then
2228 if (separate_parset) then
2229 write(licz3,"(bz,i3.3)") myparm
2230 histname=prefix(:ilen(prefix))//'_par'//licz3//'.hist'
2232 histname=prefix(:ilen(prefix))//'.hist'
2235 if (separate_parset) then
2236 write(licz3,"(bz,i3.3)") myparm
2237 histname=prefix(:ilen(prefix))//'_par'//licz3// &
2238 '_slice_'//licz2//'.hist'
2240 histname=prefix(:ilen(prefix))//'_slice_'//licz2//'.hist'
2243 #if defined(AIX) || defined(PGI)
2244 open (ihist,file=histname,position='append')
2246 open (ihist,file=histname,access='append')
2254 sumH=sumH+hfin(t,ib)
2256 if (sumH.gt.0.0d0) then
2258 jj = mod(liczbaW,nbin1)
2259 liczbaW=liczbaW/nbin1
2260 write (iout,'(f6.3,$)') dmin+(jj+0.5d0)*delta
2262 write (ihist,'(f6.3,$)') dmin+(jj+0.5d0)*delta
2265 write (iout,'(e20.10,$)') hfin(t,ib)
2266 if (histfile) write (ihist,'(e20.10,$)') hfin(t,ib)
2268 write (iout,'(i5)') iparm
2269 if (histfile) write (ihist,'(i5)') iparm
2276 if (nslice.eq.1) then
2277 if (separate_parset) then
2278 write(licz3,"(bz,i3.3)") myparm
2279 histname=prefix(:ilen(prefix))//"_par"//licz3//'.ent'
2281 histname=prefix(:ilen(prefix))//'.ent'
2284 if (separate_parset) then
2285 write(licz3,"(bz,i3.3)") myparm
2286 histname=prefix(:ilen(prefix))//'par_'//licz3// &
2287 '_slice_'//licz2//'.ent'
2289 histname=prefix(:ilen(prefix))//'_slice_'//licz2//'.ent'
2292 #if defined(AIX) || defined(PGI)
2293 open (ihist,file=histname,position='append')
2295 open (ihist,file=histname,access='append')
2297 write (ihist,'(a)') "# Microcanonical entropy"
2299 write (ihist,'(f8.0,$)') dint(potEmin)+i
2300 if (histE(i).gt.0.0e0) then
2301 write (ihist,'(f15.5,$)') dlog(histE(i))
2303 write (ihist,'(f15.5,$)') 0.0d0
2309 write (iout,*) "Microcanonical entropy"
2311 write (iout,'(f8.0,$)') dint(potEmin)+i
2312 if (histE(i).gt.0.0e0) then
2313 write (iout,'(f15.5,$)') dlog(histE(i))
2315 write (iout,'(f15.5,$)') 0.0d0
2320 if (nslice.eq.1) then
2321 if (separate_parset) then
2322 write(licz3,"(bz,i3.3)") myparm
2323 histname=prefix(:ilen(prefix))//'_par'//licz3//'.rmsrgy'
2325 histname=prefix(:ilen(prefix))//'.rmsrgy'
2328 if (separate_parset) then
2329 write(licz3,"(bz,i3.3)") myparm
2330 histname=prefix(:ilen(prefix))//'_par'//licz3// &
2331 '_slice_'//licz2//'.rmsrgy'
2333 histname=prefix(:ilen(prefix))//'_slice_'//licz2//'.rmsrgy'
2336 #if defined(AIX) || defined(PGI)
2337 open (ihist,file=histname,position='append')
2339 open (ihist,file=histname,access='append')
2343 write(ihist,'(2f8.2,$)') &
2344 rgymin+deltrgy*j,rmsmin+deltrms*i
2346 if (hrmsrgy(j,i,ib).gt.0.0d0) then
2347 write(ihist,'(e14.5,$)') &
2348 -dlog(hrmsrgy(j,i,ib))/beta_h(ib,iparm) &
2351 write(ihist,'(e14.5,$)') 1.0d6
2354 write (ihist,'(i2)') iparm
2362 call MPI_Reduce(hfin_ent_p(0),hfin_ent(0),nbin,&
2363 MPI_DOUBLE_PRECISION,MPI_SUM,Master,WHAM_COMM,IERROR)
2364 call MPI_Reduce(sumW_p(0,1),sumW(0,1),(nGridT+1)*nParmSet,&
2365 MPI_DOUBLE_PRECISION,MPI_SUM,Master,WHAM_COMM,IERROR)
2366 call MPI_Reduce(sumE_p(0,1),sumE(0,1),(nGridT+1)*nParmSet,&
2367 MPI_DOUBLE_PRECISION,MPI_SUM,Master,WHAM_COMM,IERROR)
2368 call MPI_Reduce(sumEbis_p(0,1),sumEbis(0,1),(nGridT+1)*nParmSet,&
2369 MPI_DOUBLE_PRECISION,MPI_SUM,Master,WHAM_COMM,IERROR)
2370 call MPI_Reduce(sumEsq_p(0,1),sumEsq(0,1),(nGridT+1)*nParmSet,&
2371 MPI_DOUBLE_PRECISION,MPI_SUM,Master,WHAM_COMM,IERROR)
2372 call MPI_Reduce(sumQ_p(1,0,1),sumQ(1,0,1),&
2373 MaxQ1*(nGridT+1)*nParmSet,MPI_DOUBLE_PRECISION,MPI_SUM,Master,&
2375 call MPI_Reduce(sumQsq_p(1,0,1),sumQsq(1,0,1),&
2376 MaxQ1*(nGridT+1)*nParmSet,MPI_DOUBLE_PRECISION,MPI_SUM,Master,&
2378 call MPI_Reduce(sumEQ_p(1,0,1),sumEQ(1,0,1),&
2379 MaxQ1*(nGridT+1)*nParmSet,MPI_DOUBLE_PRECISION,MPI_SUM,Master,&
2381 if (me.eq.master) then
2383 write (iout,'(/a)') 'Thermal characteristics of folding'
2384 if (nslice.eq.1) then
2387 nazwa=prefix(:ilen(prefix))//"_slice_"//licz2
2390 if (nparmset.eq.1 .and. .not.separate_parset) then
2391 nazwa=nazwa(:iln)//".thermal"
2392 else if (nparmset.eq.1 .and. separate_parset) then
2393 write(licz3,"(bz,i3.3)") myparm
2394 nazwa=nazwa(:iln)//"_par_"//licz3//".thermal"
2397 if (nparmset.gt.1) then
2398 write(licz3,"(bz,i3.3)") iparm
2399 nazwa=nazwa(:iln)//"_par_"//licz3//".thermal"
2402 if (separate_parset) then
2403 write (iout,'(a,i3)') "Parameter set",myparm
2405 write (iout,'(a,i3)') "Parameter set",iparm
2408 betaT=1.0d0/(1.987D-3*(startGridT+i*delta_T))
2409 if (betaT.ge.beta_h(1,iparm)) then
2410 potEmin=potEmin_all(1,iparm)
2411 else if (betaT.lt.beta_h(nT_h(iparm),iparm)) then
2412 potEmin=potEmin_all(nT_h(iparm),iparm)
2414 do l=1,nT_h(iparm)-1
2415 if (betaT.le.beta_h(l,iparm) .and. &
2416 betaT.gt.beta_h(l+1,iparm)) then
2417 potEmin=potEmin_all(l,iparm)
2423 sumE(i,iparm)=sumE(i,iparm)/sumW(i,iparm)
2424 sumEbis(i,iparm)=(startGridT+i*delta_T)*sumEbis(i,iparm)/ &
2426 sumEsq(i,iparm)=(sumEsq(i,iparm)/sumW(i,iparm) &
2427 -sumE(i,iparm)**2)/(1.987D-3*(startGridT+i*delta_T)**2)
2429 sumQ(j,i,iparm)=sumQ(j,i,iparm)/sumW(i,iparm)
2430 sumQsq(j,i,iparm)=sumQsq(j,i,iparm)/sumW(i,iparm) &
2432 sumEQ(j,i,iparm)=sumEQ(j,i,iparm)/sumW(i,iparm) &
2433 -sumQ(j,i,iparm)*sumE(i,iparm)
2435 sumW(i,iparm)=-dlog(sumW(i,iparm))*(1.987D-3* &
2436 (startGridT+i*delta_T))+potEmin
2437 write (iout,'(f7.1,2f15.5,$)') startGridT+i*delta_T,&
2438 sumW(i,iparm),sumE(i,iparm)
2439 write (iout,'(f10.5,$)') (sumQ(j,i,iparm),j=1,nQ+2)
2440 write (iout,'(e15.5,$)') sumEsq(i,iparm)-sumEbis(i,iparm),&
2441 (sumQsq(j,i,iparm),j=1,nQ+2),(sumEQ(j,i,iparm),j=1,nQ+2)
2443 write (34,'(f7.1,2f15.5,$)') startGridT+i*delta_T,&
2444 sumW(i,iparm),sumE(i,iparm)
2445 write (34,'(f10.5,$)') (sumQ(j,i,iparm),j=1,nQ+2)
2446 write (34,'(e15.5,$)') sumEsq(i,iparm)-sumEbis(i,iparm),&
2447 (sumQsq(j,i,iparm),j=1,nQ+2),(sumEQ(j,i,iparm),j=1,nQ+2)
2455 if (hfin_ent(t).gt.0.0d0) then
2457 jj = mod(liczbaW,nbin1)
2458 write (iout,'(f6.3,e20.10," ent")') dmin+(jj+0.5d0)*delta,&
2460 if (histfile) write (ihist,'(f6.3,e20.10," ent")') &
2461 dmin+(jj+0.5d0)*delta,&
2465 if (histfile) close(ihist)
2469 ! Write data for zscore
2470 if (nslice.eq.1) then
2471 zscname=prefix(:ilen(prefix))//".zsc"
2473 zscname=prefix(:ilen(prefix))//"_slice_"//licz2//".zsc"
2475 #if defined(AIX) || defined(PGI)
2476 open (izsc,file=prefix(:ilen(prefix))//'.zsc',position='append')
2478 open (izsc,file=prefix(:ilen(prefix))//'.zsc',access='append')
2480 write (izsc,'("NQ=",i1," NPARM=",i1)') nQ,nParmSet
2482 write (izsc,'("NT=",i1)') nT_h(iparm)
2484 write (izsc,'("TEMP=",f6.1," NR=",i2," SNK=",$)') &
2485 1.0d0/(beta_h(ib,iparm)*1.987D-3),nR(ib,iparm)
2486 jj = min0(nR(ib,iparm),7)
2487 write (izsc,'(i8,$)') (snk(i,ib,iparm,islice),i=1,jj)
2488 write (izsc,'(a1,$)') (" ",i=22+8*jj+1,79)
2489 write (izsc,'("&")')
2490 if (nR(ib,iparm).gt.7) then
2491 do ii=8,nR(ib,iparm),9
2492 jj = min0(nR(ib,iparm),ii+8)
2493 write (izsc,'(i8,$)') (snk(i,ib,iparm,islice),i=ii,jj)
2494 write (izsc,'(a1,$') (" ",i=(jj-ii+1)*8+1,79)
2495 write (izsc,'("&")')
2498 write (izsc,'("FI=",$)')
2499 jj=min0(nR(ib,iparm),7)
2500 write (izsc,'(f10.5,$)') (fi(i,ib,iparm),i=1,jj)
2501 write (izsc,'(a1,$)') (" ",i=3+10*jj+1,79)
2502 write (izsc,'("&")')
2503 if (nR(ib,iparm).gt.7) then
2504 do ii=8,nR(ib,iparm),9
2505 jj = min0(nR(ib,iparm),ii+8)
2506 write (izsc,'(f10.5,$)') (fi(i,ib,iparm),i=ii,jj)
2507 if (jj.eq.nR(ib,iparm)) then
2510 write (izsc,'(a1,$)') (" ",i=10*(jj-ii+1)+1,79)
2511 write (izsc,'(t80,"&")')
2516 write (izsc,'("KH=",$)')
2517 write (izsc,'(f7.2,$)') (Kh(j,i,ib,iparm),j=1,nQ)
2518 write (izsc,'(" Q0=",$)')
2519 write (izsc,'(f7.5,$)') (q0(j,i,ib,iparm),j=1,nQ)
2530 end subroutine WHAMCALC
2531 !-----------------------------------------------------------------------------
2532 end module wham_calc