1 subroutine WHAM_CALC(islice,*)
2 ! Weighed Histogram Analysis Method (WHAM) code
3 ! Written by A. Liwo based on the work of Kumar et al.,
4 ! J.Comput.Chem., 13, 1011 (1992)
6 ! 2/1/05 Multiple temperatures allowed.
7 ! 2/2/05 Free energies calculated directly from data points
8 ! acc. to Eq. (21) of Kumar et al.; final histograms also
9 ! constructed based on this equation.
10 ! 2/12/05 Multiple parameter sets included
12 ! 2/2/05 Parallel version
15 include "DIMENSIONS.ZSCOPT"
16 include "DIMENSIONS.FREE"
18 parameter (NGridT=400)
19 integer MaxBinRms,MaxBinRgy
20 parameter (MaxBinRms=100,MaxBinRgy=100)
22 c parameter (MaxHdim=200000)
23 parameter (MaxHdim=200)
25 parameter (maxinde=200)
29 integer ierror,errcode,status(MPI_STATUS_SIZE)
31 include "COMMON.CONTROL"
32 include "COMMON.IOUNITS"
34 include "COMMON.ENERGIES"
35 include "COMMON.FFIELD"
36 include "COMMON.SBRIDGE"
38 include "COMMON.ENEPS"
39 include "COMMON.SHIELD"
40 integer MaxPoint,MaxPointProc
41 parameter (MaxPoint=MaxStr,
42 & MaxPointProc=MaxStr_Proc)
43 double precision finorm_max,potfac,entmin,entmax,expfac,vf
44 parameter (finorm_max=1.0d0)
46 integer i,ii,j,jj,k,kk,l,m,ind,iter,t,tmax,ient,ientmax,iln
47 integer start,end,iharm,ib,iib,nbin1,nbin,nbin_rms,nbin_rgy,
48 & nbin_rmsrgy,liczba,iparm,nFi,indrgy,indrms
49 integer htot(0:MaxHdim),histent(0:2000)
50 double precision v(MaxPointProc,MaxR,MaxT_h,Max_Parm)
51 double precision energia(0:max_ene)
53 integer tmax_t,upindE_p
54 double precision fi_p(MaxR,MaxT_h,Max_Parm),
55 & fi_p_min(MaxR,MaxT_h,Max_Parm)
56 double precision sumW_p(0:nGridT,Max_Parm),
57 & sumE_p(0:nGridT,Max_Parm),sumEsq_p(0:nGridT,Max_Parm),
58 & sumQ_p(MaxQ1,0:nGridT,Max_Parm),
59 & sumQsq_p(MaxQ1,0:nGridT,Max_Parm),
60 & sumEQ_p(MaxQ1,0:nGridT,Max_Parm),
61 & sumEprim_p(MaxQ1,0:nGridT,Max_Parm),
62 & sumEbis_p(0:nGridT,Max_Parm)
63 double precision hfin_p(0:MaxHdim,maxT_h),
64 & hfin_ent_p(0:MaxHdim),histE_p(0:maxindE),sumH,
65 & hrmsrgy_p(0:MaxBinRgy,0:MaxBinRms,maxT_h)
66 double precision rgymin_t,rmsmin_t,rgymax_t,rmsmax_t
67 double precision potEmin_t,entmin_p,entmax_p,
68 & potEmin_t_all(maxT_h,Max_Parm)
69 integer histent_p(0:2000)
70 logical lprint /.true./
72 double precision delta_T /1.0d0/
73 double precision rgymin,rmsmin,rgymax,rmsmax
74 double precision sumW(0:NGridT,Max_Parm),sumE(0:NGridT,Max_Parm),
75 & sumEsq(0:NGridT,Max_Parm),sumQ(MaxQ1,0:NGridT,Max_Parm),
76 & sumQsq(MaxQ1,0:NGridT,Max_Parm),sumEQ(MaxQ1,0:NGridT,Max_Parm),
77 & sumEprim(0:NGridT,Max_Parm),sumEbis(0:NGridT,Max_Parm),betaT,
79 double precision fi(MaxR,maxT_h,Max_Parm),
80 & fi_min(MaxR,maxT_h,Max_Parm),
81 & dd,dd1,dd2,hh,dmin,denom,finorm,avefi,pom,
82 & hfin(0:MaxHdim,maxT_h),histE(0:maxindE),
83 & hrmsrgy(0:MaxBinRgy,0:MaxBinRms,maxT_h),
84 & potEmin,ent,potEmin_all(maxT_h,Max_Parm),potEmin_min,
85 & hfin_ent(0:MaxHdim),vmax,aux,entfac_min
86 double precision fT(6),fTprim(6),fTbis(6),quot,quotl1,quotl,kfacl,
87 & eprim,ebis,temper,kfac/2.4d0/,T0/300.0d0/,startGridT/200.0d0/,
88 & eplus,eminus,logfac,tanhT,tt
89 double precision etot,evdw,evdw_t,evdw2,ees,evdw1,ebe,etors,
90 & escloc,ehpb,ecorr,ecorr5,ecorr6,eello_turn4,eello_turn3,
91 & eturn6,eel_loc,edihcnstr,etors_d,estr,evdw2_14,esccor,
94 integer ind_point(maxpoint),upindE,indE
103 write(licz2,'(bz,i2.2)') islice
105 write (iout,'(//80(1h-)/"Solving WHAM equations for slice",
106 & i2/80(1h-)//)') islice
107 write (iout,*) "delta",delta," nbin1",nbin1
108 write (iout,*) "MaxN",MaxN," MaxQ",MaxQ," MaHdim",MaxHdim
115 potEmin_all(j,i)=1.0d10
132 C if (potE(i,j).le.potEmin) potEmin=potE(i,j)
134 if (q(nQ+1,i).lt.rmsmin) rmsmin=q(nQ+1,i)
135 if (q(nQ+1,i).gt.rmsmax) rmsmax=q(nQ+1,i)
136 if (q(nQ+2,i).lt.rgymin) rgymin=q(nQ+2,i)
137 if (q(nQ+2,i).gt.rgymax) rgymax=q(nQ+2,i)
140 ind=(q(j,i)-dmin+1.0d-8)/delta
142 ind_point(i)=ind_point(i)+ind
144 ind_point(i)=ind_point(i)+nbin1**(j-1)*ind
146 c write (iout,*) "i",i," j",j," q",q(j,i)," ind_point",
149 if (ind_point(i).lt.0 .or. ind_point(i).gt.MaxHdim) then
150 write (iout,*) "Error - index exceeds range for point",i,
151 & " q=",q(j,i)," ind",ind_point(i)
153 write (iout,*) "Processor",me1
155 call MPI_Abort(MPI_COMM_WORLD, Ierror, Errcode )
160 if (ind_point(i).gt.tmax) tmax=ind_point(i)
161 htot(ind_point(i))=htot(ind_point(i))+1
163 write (iout,*) "i",i,"q",(q(j,i),j=1,nQ)," ind",ind_point(i),
164 & " htot",htot(ind_point(i))
171 write (iout,'(a)') "Numbers of counts in Q bins"
173 if (htot(t).gt.0) then
174 write (iout,'(i15,$)') t
177 jj = mod(liczba,nbin1)
179 write (iout,'(i5,$)') jj
181 write (iout,'(i8)') htot(t)
185 write (iout,'(a,i3)') "Number of data points for parameter set",
187 write (iout,'(i7,$)') ((snk(m,ib,iparm,islice),m=1,nr(ib,iparm)),
189 write (iout,'(i8)') stot(islice)
195 call MPI_AllReduce(tmax,tmax_t,1,MPI_INTEGER,MPI_MAX,
198 C call MPI_AllReduce(potEmin,potEmin_t,1,MPI_DOUBLE_PRECISION,
199 C & MPI_MIN,WHAM_COMM,IERROR)
200 call MPI_AllReduce(rmsmin,rmsmin_t,1,MPI_DOUBLE_PRECISION,
201 & MPI_MIN,WHAM_COMM,IERROR)
202 call MPI_AllReduce(rmsmax,rmsmax_t,1,MPI_DOUBLE_PRECISION,
203 & MPI_MAX,WHAM_COMM,IERROR)
204 call MPI_AllReduce(rgymin,rgymin_t,1,MPI_DOUBLE_PRECISION,
205 & MPI_MIN,WHAM_COMM,IERROR)
206 call MPI_AllReduce(rgymax,rgymax_t,1,MPI_DOUBLE_PRECISION,
207 & MPI_MAX,WHAM_COMM,IERROR)
208 C potEmin=potEmin_t/2
214 write (iout,*) "potEmin",potEmin
216 rmsmin=deltrms*dint(rmsmin/deltrms)
217 rmsmax=deltrms*dint(rmsmax/deltrms)
218 rgymin=deltrms*dint(rgymin/deltrgy)
219 rgymax=deltrms*dint(rgymax/deltrgy)
220 nbin_rms=(rmsmax-rmsmin)/deltrms
221 nbin_rgy=(rgymax-rgymin)/deltrgy
222 write (iout,*) "rmsmin",rmsmin," rmsmax",rmsmax," rgymin",rgymin,
223 & " rgymax",rgymax," nbin_rms",nbin_rms," nbin_rgy",nbin_rgy
230 write (iout,*) "nFi",nFi
231 ! Compute the Boltzmann factor corresponing to restrain potentials in different
238 c write (9,'(3i5,f10.5)') i,(iparm,potE(i,iparm),iparm=1,nParmSet)
241 write (iout,'(2i5,21f8.2)') i,iparm,
242 & (enetb(k,i,iparm),k=1,max_ene)
244 call restore_parm(iparm)
246 write (iout,*) wsc,wscp,welec,wvdwpp,wang,wtor,wscloc,
247 & wcorr,wcorr5,wcorr6,wturn4,wturn3,wturn6,wel_loc,
248 & wtor_d,wsccor,wbond
251 if (rescale_mode.eq.1) then
252 quot=1.0d0/(beta_h(ib,iparm)*1.987D-3*T0)
259 fT(l)=kfacl/(kfacl-1.0d0+quotl)
262 tt = 1.0d0/(beta_h(ib,iparm)*1.987D-3)
263 ft(6)=(320.0d0+80.0d0*dtanh((tt-320.0d0)/80.0d0))/320.0d0
265 ft(6)=1.0d0/(beta_h(ib,iparm)*1.987D-3*T0)
269 else if (rescale_mode.eq.2) then
270 quot=1.0d0/(T0*beta_h(ib,iparm)*1.987D-3)
274 fT(l)=1.12692801104297249644d0/
275 & dlog(dexp(quotl)+dexp(-quotl))
278 tt = 1.0d0/(beta_h(ib,iparm)*1.987D-3)
279 ft(6)=(320.0d0+80.0d0*dtanh((tt-320.0d0)/80.0d0))/320.0d0
281 ft(6)=1.0d0/(beta_h(ib,iparm)*1.987D-3*T0)
285 c write (iout,*) 1.0d0/(beta_h(ib,iparm)*1.987D-3),ft
286 else if (rescale_mode.eq.0) then
291 write (iout,*) "Error in WHAM_CALC: wrong RESCALE_MODE",
296 evdw=enetb(1,i,iparm)
297 evdw_t=enetb(21,i,iparm)
299 evdw2_14=enetb(17,i,iparm)
300 evdw2=enetb(2,i,iparm)+evdw2_14
302 evdw2=enetb(2,i,iparm)
307 evdw1=enetb(16,i,iparm)
312 ecorr=enetb(4,i,iparm)
313 ecorr5=enetb(5,i,iparm)
314 ecorr6=enetb(6,i,iparm)
315 eel_loc=enetb(7,i,iparm)
316 eello_turn3=enetb(8,i,iparm)
317 eello_turn4=enetb(9,i,iparm)
318 eturn6=enetb(10,i,iparm)
319 ebe=enetb(11,i,iparm)
320 escloc=enetb(12,i,iparm)
321 etors=enetb(13,i,iparm)
322 etors_d=enetb(14,i,iparm)
323 ehpb=enetb(15,i,iparm)
324 estr=enetb(18,i,iparm)
325 esccor=enetb(19,i,iparm)
326 edihcnstr=enetb(20,i,iparm)
327 eliptran=enetb(22,i,iparm)
328 etube=enetb(25,i,iparm)
332 write (iout,'(3i5,6f5.2,14f12.3)') i,ib,iparm,(ft(l),l=1,6),
333 & evdw+evdw_t,evdw2,ees,evdw1,ecorr,eel_loc,estr,ebe,escloc,
334 & etors,etors_d,eello_turn3,eello_turn4,esccor
338 if (shield_mode.gt.0) then
339 etot=ft(1)*wsc*(evdw+ft(6)*evdw_t)+ft(1)*wscp*evdw2
341 & +ft(1)*wvdwpp*evdw1
342 & +wang*ebe+ft(1)*wtor*etors+wscloc*escloc
343 & +wstrain*ehpb+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5
344 & +ft(5)*wcorr6*ecorr6+ft(3)*wturn4*eello_turn4
345 & +ft(2)*wturn3*eello_turn3
346 & +ft(5)*wturn6*eturn6+ft(2)*wel_loc*eel_loc
347 & +edihcnstr+ft(2)*wtor_d*etors_d+ft(1)*wsccor*esccor
348 & +wbond*estr+wliptran*eliptran+wtube*Etube
350 etot=wsc*(evdw+ft(6)*evdw_t)+wscp*evdw2+ft(1)*welec*ees
352 & +wang*ebe+ft(1)*wtor*etors+wscloc*escloc
353 & +wstrain*ehpb+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5
354 & +ft(5)*wcorr6*ecorr6+ft(3)*wturn4*eello_turn4
355 & +ft(2)*wturn3*eello_turn3
356 & +ft(5)*wturn6*eturn6+ft(2)*wel_loc*eel_loc
357 & +edihcnstr+ft(2)*wtor_d*etors_d+ft(1)*wsccor*esccor
358 & +wbond*estr+wliptran*eliptran+wtube*Etube
361 if (shield_mode.gt.0) then
362 etot=ft(1)*wsc*(evdw+ft(6)*evdw_t)+ft(1)*wscp*evdw2
363 & +ft(1)*welec*(ees+evdw1)
364 & +wang*ebe+ft(1)*wtor*etors+wscloc*escloc
365 & +wstrain*ehpb+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5
366 & +ft(5)*wcorr6*ecorr6+ft(3)*wturn4*eello_turn4
367 & +ft(2)*wturn3*eello_turn3
368 & +ft(5)*wturn6*eturn6+ft(2)*wel_loc*eel_loc+edihcnstr
369 & +ft(2)*wtor_d*etors_d+ft(1)*wsccor*esccor
370 & +wbond*estr+wliptran*eliptran+wtube*Etube
372 etot=wsc*(evdw+ft(6)*evdw_t)+wscp*evdw2
373 & +ft(1)*welec*(ees+evdw1)
374 & +wang*ebe+ft(1)*wtor*etors+wscloc*escloc
375 & +wstrain*ehpb+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5
376 & +ft(5)*wcorr6*ecorr6+ft(3)*wturn4*eello_turn4
377 & +ft(2)*wturn3*eello_turn3
378 & +ft(5)*wturn6*eturn6+ft(2)*wel_loc*eel_loc+edihcnstr
379 & +ft(2)*wtor_d*etors_d+ft(1)*wsccor*esccor
380 & +wbond*estr+wliptran*eliptran+wtube*Etube
385 write (iout,*) i,iparm,1.0d0/(beta_h(ib,iparm)*1.987D-3),
389 if (iparm.eq.1 .and. ib.eq.1) then
390 write (iout,*)"Conformation",i
393 energia(k)=enetb(k,i,iparm)
395 call enerprint(energia(0),fT)
402 Econstr=Econstr+Kh(j,kk,ib,iparm)
403 & *(dd-q0(j,kk,ib,iparm))**2
406 & -beta_h(ib,iparm)*(etot+Econstr)
408 write (iout,'(4i5,4e15.5)') i,kk,ib,iparm,
409 & etot,potEmin,etot-potEmin,v(i,kk,ib,iparm)
415 ! Simple iteration to calculate free energies corresponding to all simulation
419 ! Compute new free-energy values corresponding to the righ-hand side of the
420 ! equation and their derivatives.
421 write (iout,*) "------------------------fi"
433 vf=v(t,l,k,i)+f(l,k,i)
434 if (vf.gt.vmax) vmax=vf
442 aux=f(l,k,i)+v(t,l,k,i)-vmax
444 & denom=denom+snk(l,k,i,islice)*dexp(aux)
448 entfac(t)=-dlog(denom)-vmax
449 if (entfac(t).lt.entfac_min) entfac_min=entfac(t)
451 write (iout,*) t,"vmax",vmax," denom",denom,"entfac",entfac(t)
457 do ii=1,nR(iib,iparm)
459 fi_p_min(ii,iib,iparm)=-1.0d10
461 aux=v(t,ii,iib,iparm)+entfac(t)
462 if (aux.gt.fi_p_min(ii,iib,iparm))
463 & fi_p_min(ii,iib,iparm)=aux
467 aux=v(t,ii,iib,iparm)+entfac(t)
468 if (aux.gt.fi_min(ii,iib,iparm))
469 & fi_min(ii,iib,iparm)=aux
477 write (iout,*) "fi_min before AllReduce"
480 write (iout,*) (i,j,k,fi_p_min(k,j,i),k=1,nR(j,i))
484 call MPI_AllReduce(fi_p_min,fi_min,MaxR*MaxT_h*nParmSet,
485 & MPI_DOUBLE_PRECISION,MPI_MAX,WHAM_COMM,IERROR)
487 write (iout,*) "fi_min after AllReduce"
490 write (iout,*) (i,j,k,fi_min(k,j,i),k=1,nR(j,i))
497 do ii=1,nR(iib,iparm)
499 fi_p(ii,iib,iparm)=0.0d0
501 fi_p(ii,iib,iparm)=fi_p(ii,iib,iparm)
502 & +dexp(v(t,ii,iib,iparm)+entfac(t)-fi_min(ii,iib,iparm))
504 write (iout,'(4i5,4e15.5)') t,ii,iib,iparm,
505 & v(t,ii,iib,iparm),entfac(t),fi_min(ii,iib,iparm),
510 fi(ii,iib,iparm)=0.0d0
512 fi(ii,iib,iparm)=fi(ii,iib,iparm)
513 & +dexp(v(t,ii,iib,iparm)+entfac(t)-fi_min(ii,iib,iparm))
521 write (iout,*) "fi before MPI_Reduce me",me,' master',master
523 do ib=1,nT_h(nparmset)
524 write (iout,*) "iparm",iparm," ib",ib
525 write (iout,*) "beta=",beta_h(ib,iparm)
526 write (iout,'(8e15.5)') (fi_p(i,ib,iparm),i=1,nR(ib,iparm))
530 write (iout,*) "REDUCE size",maxR,MaxT_h,nParmSet,
531 & maxR*MaxT_h*nParmSet
532 write (iout,*) "MPI_COMM_WORLD",MPI_COMM_WORLD,
533 & " WHAM_COMM",WHAM_COMM
534 call MPI_Reduce(fi_p(1,1,1),fi(1,1,1),maxR*MaxT_h*nParmSet,
535 & MPI_DOUBLE_PRECISION,
536 & MPI_SUM,Master,WHAM_COMM,IERROR)
538 write (iout,*) "fi after MPI_Reduce nparmset",nparmset
540 write (iout,*) "iparm",iparm
542 write (iout,*) "beta=",beta_h(ib,iparm)
543 write (iout,'(8e15.5)') (fi(i,ib,iparm),i=1,nR(ib,iparm))
547 if (me1.eq.Master) then
553 fi(i,ib,iparm)=-dlog(fi(i,ib,iparm))-fi_min(i,ib,iparm)
554 avefi=avefi+fi(i,ib,iparm)
560 write (iout,*) "Parameter set",iparm
562 write (iout,*) "beta=",beta_h(ib,iparm)
564 fi(i,ib,iparm)=fi(i,ib,iparm)-avefi
566 write (iout,'(8f10.5)') (fi(i,ib,iparm),i=1,nR(ib,iparm))
567 write (iout,'(8f10.5)') (f(i,ib,iparm),i=1,nR(ib,iparm))
571 ! Compute the norm of free-energy increments.
576 finorm=finorm+dabs(fi(i,ib,iparm)-f(i,ib,iparm))
577 f(i,ib,iparm)=fi(i,ib,iparm)
582 write (iout,*) 'Iteration',iter,' finorm',finorm
586 call MPI_Bcast(f(1,1,1),MaxR*MaxT_h*nParmSet,
587 & MPI_DOUBLE_PRECISION,Master,
589 call MPI_Bcast(finorm,1,MPI_DOUBLE_PRECISION,Master,
592 ! Exit, if the increment norm is smaller than pre-assigned tolerance.
593 if (finorm.lt.fimin) then
594 write (iout,*) 'Iteration converged'
607 c write (9,'(3i5,f10.5)') i,(iparm,potE(i,iparm),iparm=1,nParmSet)
610 write (iout,'(2i5,21f8.2)') i,iparm,
611 & (enetb(k,i,iparm),k=1,21)
613 call restore_parm(iparm)
615 write (iout,*) wsc,wscp,welec,wvdwpp,wang,wtor,wscloc,
616 & wcorr,wcorr5,wcorr6,wturn4,wturn3,wturn6,wel_loc,
617 & wtor_d,wsccor,wbond
620 if (rescale_mode.eq.1) then
621 quot=1.0d0/(beta_h(ib,iparm)*1.987D-3*T0)
628 fT(l)=kfacl/(kfacl-1.0d0+quotl)
631 tt = 1.0d0/(beta_h(ib,iparm)*1.987D-3)
632 ft(6)=(320.0d0+80.0d0*dtanh((tt-320.0d0)/80.0d0))/320.0d0
634 ft(6)=1.0d0/(beta_h(ib,iparm)*1.987D-3*T0)
638 else if (rescale_mode.eq.2) then
639 quot=1.0d0/(T0*beta_h(ib,iparm)*1.987D-3)
643 fT(l)=1.12692801104297249644d0/
644 & dlog(dexp(quotl)+dexp(-quotl))
647 tt = 1.0d0/(beta_h(ib,iparm)*1.987D-3)
648 ft(6)=(320.0d0+80.0d0*dtanh((tt-320.0d0)/80.0d0))/320.0d0
650 ft(6)=1.0d0/(beta_h(ib,iparm)*1.987D-3*T0)
654 c write (iout,*) 1.0d0/(beta_h(ib,iparm)*1.987D-3),ft
655 else if (rescale_mode.eq.0) then
660 write (iout,*) "Error in WHAM_CALC: wrong RESCALE_MODE",
665 evdw=enetb(21,i,iparm)
666 evdw_t=enetb(1,i,iparm)
668 evdw2_14=enetb(17,i,iparm)
669 evdw2=enetb(2,i,iparm)+evdw2_14
671 evdw2=enetb(2,i,iparm)
676 evdw1=enetb(16,i,iparm)
681 ecorr=enetb(4,i,iparm)
682 ecorr5=enetb(5,i,iparm)
683 ecorr6=enetb(6,i,iparm)
684 eel_loc=enetb(7,i,iparm)
685 eello_turn3=enetb(8,i,iparm)
686 eello_turn4=enetb(9,i,iparm)
687 eturn6=enetb(10,i,iparm)
688 ebe=enetb(11,i,iparm)
689 escloc=enetb(12,i,iparm)
690 etors=enetb(13,i,iparm)
691 etors_d=enetb(14,i,iparm)
692 ehpb=enetb(15,i,iparm)
693 estr=enetb(18,i,iparm)
694 esccor=enetb(19,i,iparm)
695 edihcnstr=enetb(20,i,iparm)
697 eliptran=enetb(22,i,iparm)
698 etube=enetb(25,i,iparm)
701 if (shield_mode.gt.0) then
702 etot=ft(1)*wsc*(evdw+ft(6)*evdw_t)+ft(1)*wscp*evdw2
704 & +ft(1)*wvdwpp*evdw1
705 & +wang*ebe+ft(1)*wtor*etors+wscloc*escloc
706 & +wstrain*ehpb+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5
707 & +ft(5)*wcorr6*ecorr6+ft(3)*wturn4*eello_turn4
708 & +ft(2)*wturn3*eello_turn3
709 & +ft(5)*wturn6*eturn6+ft(2)*wel_loc*eel_loc
710 & +edihcnstr+ft(2)*wtor_d*etors_d+ft(1)*wsccor*esccor
711 & +wbond*estr+wliptran*eliptran+wtube*Etube
713 etot=wsc*(evdw+ft(6)*evdw_t)+wscp*evdw2+ft(1)*welec*ees
715 & +wang*ebe+ft(1)*wtor*etors+wscloc*escloc
716 & +wstrain*ehpb+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5
717 & +ft(5)*wcorr6*ecorr6+ft(3)*wturn4*eello_turn4
718 & +ft(2)*wturn3*eello_turn3
719 & +ft(5)*wturn6*eturn6+ft(2)*wel_loc*eel_loc
720 & +edihcnstr+ft(2)*wtor_d*etors_d+ft(1)*wsccor*esccor
721 & +wbond*estr+wliptran*eliptran+wtube*Etube
724 if (shield_mode.gt.0) then
725 etot=ft(1)*wsc*(evdw+ft(6)*evdw_t)+ft(1)*wscp*evdw2
726 & +ft(1)*welec*(ees+evdw1)
727 & +wang*ebe+ft(1)*wtor*etors+wscloc*escloc
728 & +wstrain*ehpb+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5
729 & +ft(5)*wcorr6*ecorr6+ft(3)*wturn4*eello_turn4
730 & +ft(2)*wturn3*eello_turn3
731 & +ft(5)*wturn6*eturn6+ft(2)*wel_loc*eel_loc+edihcnstr
732 & +ft(2)*wtor_d*etors_d+ft(1)*wsccor*esccor
733 & +wbond*estr+wliptran*eliptran+wtube*Etube
735 etot=wsc*(evdw+ft(6)*evdw_t)+wscp*evdw2
736 & +ft(1)*welec*(ees+evdw1)
737 & +wang*ebe+ft(1)*wtor*etors+wscloc*escloc
738 & +wstrain*ehpb+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5
739 & +ft(5)*wcorr6*ecorr6+ft(3)*wturn4*eello_turn4
740 & +ft(2)*wturn3*eello_turn3
741 & +ft(5)*wturn6*eturn6+ft(2)*wel_loc*eel_loc+edihcnstr
742 & +ft(2)*wtor_d*etors_d+ft(1)*wsccor*esccor
743 & +wbond*estr+wliptran*eliptran+wtube*Etube
747 etot=etot-entfac(i)/beta_h(ib,iparm)
748 if(etot.lt.potEmin_all(ib,iparm)) potEmin_all(ib,iparm)=etot
756 write (iout,*) "The potEmin array before reduction"
758 write (iout,*) "Parameter set",i
760 write (iout,*) j,PotEmin_all(j,i)
763 write (iout,*) "potEmin_min",potEmin_min
766 C Determine the minimum energes for all parameter sets and temperatures
767 call MPI_AllReduce(potEmin_all(1,1),potEmin_t_all(1,1),
768 & maxT_h*nParmSet,MPI_DOUBLE_PRECISION,MPI_MIN,WHAM_COMM,IERROR)
771 potEmin_all(j,i)=potEmin_t_all(j,i)
775 potEmin_min=potEmin_all(1,1)
778 if (potEmin_all(j,i).lt.potEmin_min)
779 & potEmin_min=potEmin_all(j,i)
783 write (iout,*) "The potEmin array"
785 write (iout,*) "Parameter set",i
787 write (iout,*) j,1.0d0/(1.987d-3*beta_h(j,i)),
791 write (iout,*) "potEmin_min",potEmin_min
795 ! Now, put together the histograms from all simulations, in order to get the
796 ! unbiased total histogram.
806 write (iout,*) "--------------hist"
810 sumW_p(i,iparm)=0.0d0
811 sumE_p(i,iparm)=0.0d0
812 sumEbis_p(i,iparm)=0.0d0
813 sumEsq_p(i,iparm)=0.0d0
815 sumQ_p(j,i,iparm)=0.0d0
816 sumQsq_p(j,i,iparm)=0.0d0
817 sumEQ_p(j,i,iparm)=0.0d0
827 sumEbis(i,iparm)=0.0d0
828 sumEsq(i,iparm)=0.0d0
830 sumQ(j,i,iparm)=0.0d0
831 sumQsq(j,i,iparm)=0.0d0
832 sumEQ(j,i,iparm)=0.0d0
838 c 8/26/05 entropy distribution
843 c ent=-dlog(entfac(t))
845 if (ent.lt.entmin_p) entmin_p=ent
846 if (ent.gt.entmax_p) entmax_p=ent
848 write (iout,*) "entmin",entmin_p," entmax",entmax_p
850 call MPI_Allreduce(entmin_p,entmin,1,MPI_DOUBLE_PRECISION,MPI_MIN,
852 call MPI_Allreduce(entmax_p,entmax,1,MPI_DOUBLE_PRECISION,MPI_MAX,
854 C ientmax=entmax-entmin
855 C if (ientmax.gt.2000) ientmax=2000
856 if ((-dlog(entmax)-entmin).lt.2000.0d0) then
857 ientmax=-dlog(entmax)-entmin
861 write (iout,*) "entmin",entmin," entmax",entmax," ientmax",ientmax
864 c ient=-dlog(entfac(t))-entmin
865 ient=entfac(t)-entmin
866 write (iout,*) "ient",ient,entfac(t),entmin
867 C if (ient.le.2000) histent_p(ient)=histent_p(ient)+1
869 C call MPI_Allreduce(histent_p(0),histent(0),ientmax+1,MPI_INTEGER,
870 C & MPI_SUM,WHAM_COMM,IERROR)
871 C if (me1.eq.Master) then
872 C write (iout,*) "Entropy histogram"
874 C write(iout,'(f15.4,i10)') entmin+i,histent(i)
882 if (ent.lt.entmin) entmin=ent
883 if (ent.gt.entmax) entmax=ent
885 if ((-dlog(entmax)-entmin).lt.2000.0d0) then
886 ientmax=-dlog(entmax)-entmin
890 C do t=1,ntot(islice)
891 C ient=entfac(t)-entmin
892 C if (ient.le.2000) histent(ient)=histent(ient)+1
894 C write (iout,*) "Entropy histogram"
896 C write(iout,'(2f15.4)') entmin+i,histent(i)
901 c write (iout,*) "me1",me1," scount",scount(me1)
927 hrmsrgy(j,i,ib)=0.0d0
929 hrmsrgy_p(j,i,ib)=0.0d0
941 hfin_ent_p(ind)=hfin_ent_p(ind)+dexp(entfac(t))
943 hfin_ent(ind)=hfin_ent(ind)+dexp(entfac(t))
945 c write (iout,'(2i5,20f8.2)') t,t,(enetb(k,t,iparm),k=1,18)
946 call restore_parm(iparm)
947 evdw=enetb(21,t,iparm)
948 evdw_t=enetb(1,t,iparm)
950 evdw2_14=enetb(17,t,iparm)
951 evdw2=enetb(2,t,iparm)+evdw2_14
953 evdw2=enetb(2,t,iparm)
958 evdw1=enetb(16,t,iparm)
963 ecorr=enetb(4,t,iparm)
964 ecorr5=enetb(5,t,iparm)
965 ecorr6=enetb(6,t,iparm)
966 eel_loc=enetb(7,t,iparm)
967 eello_turn3=enetb(8,t,iparm)
968 eello_turn4=enetb(9,t,iparm)
969 eturn6=enetb(10,t,iparm)
970 ebe=enetb(11,t,iparm)
971 escloc=enetb(12,t,iparm)
972 etors=enetb(13,t,iparm)
973 etors_d=enetb(14,t,iparm)
974 ehpb=enetb(15,t,iparm)
975 estr=enetb(18,t,iparm)
976 esccor=enetb(19,t,iparm)
977 edihcnstr=enetb(20,t,iparm)
979 eliptran=enetb(22,t,iparm)
980 etube=enetb(25,t,iparm)
983 betaT=startGridT+k*delta_T
987 if (rescale_mode.eq.1) then
995 denom=kfacl-1.0d0+quotl
997 ftprim(l)=-l*ft(l)*quotl1/(T0*denom)
998 ftbis(l)=l*kfacl*quotl1*
999 & (2*l*quotl-(l-1)*denom)/(quot*t0*t0*denom**3)
1002 ft(6)=(320.0d0+80.0d0*dtanh((betaT-320.0d0)/80.0d0))/
1004 ftprim(6)=1.0d0/(320.0d0*dcosh((betaT-320.0d0)/80.0d0)**2)
1005 ftbis(6)=-2.0d0*dtanh((betaT-320.0d0)/80.0d0)
1006 & /(320.0d0*80.0d0*dcosh((betaT-320.0d0)/80.0d0)**3)
1007 #elif defined(FUNCT)
1016 else if (rescale_mode.eq.2) then
1024 logfac=1.0d0/dlog(eplus+eminus)
1025 tanhT=(eplus-eminus)/(eplus+eminus)
1026 fT(l)=1.12692801104297249644d0*logfac
1027 ftprim(l)=-l*quotl1*ft(l)*tanhT*logfac/T0
1028 ftbis(l)=(l-1)*ftprim(l)/(quot*T0)-
1029 & 2*l*quotl1/T0*logfac*
1030 & (2*l*quotl1*ft(l)/(T0*(eplus+eminus)**2)
1034 ft(6)=(320.0d0+80.0d0*dtanh((betaT-320.0d0)/80.0d0))/
1036 ftprim(6)=1.0d0/(320.0d0*dcosh((betaT-320.0d0)/80.0d0)**2)
1037 ftbis(6)=-2.0d0*dtanh((betaT-320.0d0)/80.0d0)
1038 & /(320.0d0*80.0d0*dcosh((betaT-320.0d0)/80.0d0)**3)
1039 #elif defined(FUNCT)
1048 else if (rescale_mode.eq.0) then
1054 write (iout,*) "Error in WHAM_CALC: wrong RESCALE_MODE",
1059 c write (iout,*) "ftprim",ftprim
1060 c write (iout,*) "ftbis",ftbis
1061 betaT=1.0d0/(1.987D-3*betaT)
1062 if (betaT.ge.beta_h(1,iparm)) then
1063 potEmin=potEmin_all(1,iparm)+
1064 & (potEmin_all(1,iparm)-potEmin_all(2,iparm))/
1065 & (1.0/beta_h(1,iparm)-1.0/beta_h(2,iparm))*
1066 & (1.0/betaT-1.0/beta_h(1,iparm))
1068 write(iout,*) "first",temper,potEmin
1070 else if (betaT.le.beta_h(nT_h(iparm),iparm)) then
1071 potEmin=potEmin_all(nT_h(iparm),iparm)+
1072 &(potEmin_all(nT_h(iparm),iparm)-potEmin_all(nT_h(iparm)-1,iparm))/
1073 &(1.0/beta_h(nT_h(iparm),iparm)-1.0/beta_h(nT_h(iparm)-1,iparm))*
1074 &(1.0/betaT-1.0/beta_h(nt_h(iparm),iparm))
1076 write (iout,*) "last",temper,potEmin
1079 do l=1,nT_h(iparm)-1
1080 if (betaT.le.beta_h(l,iparm) .and.
1081 & betaT.gt.beta_h(l+1,iparm)) then
1082 potEmin=potEmin_all(l,iparm)
1084 write (iout,*) "l",l,
1085 & betaT,1.0d0/(1.987D-3*beta_h(l,iparm)),
1086 & 1.0d0/(1.987D-3*beta_h(l+1,iparm)),temper,potEmin
1093 if (shield_mode.gt.0) then
1094 etot=ft(1)*wsc*(evdw+ft(6)*evdw_t)+ft(1)*wscp*evdw2
1096 & +ft(1)*wvdwpp*evdw1
1097 & +wang*ebe+ft(1)*wtor*etors+wscloc*escloc
1098 & +wstrain*ehpb+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5
1099 & +ft(5)*wcorr6*ecorr6+ft(3)*wturn4*eello_turn4
1100 & +ft(2)*wturn3*eello_turn3
1101 & +ft(5)*wturn6*eturn6+ft(2)*wel_loc*eel_loc
1102 & +edihcnstr+ft(2)*wtor_d*etors_d+ft(1)*wsccor*esccor
1103 & +wbond*estr+wliptran*eliptran+wtube*Etube
1104 eprim=ftprim(1)*(ft(6)*evdw_t+evdw)
1105 C & +ftprim(6)*evdw_t
1106 & +ftprim(1)*wscp*evdw2
1107 & +ftprim(1)*welec*ees
1108 & +ftprim(1)*wvdwpp*evdw1
1109 & +ftprim(1)*wtor*etors+
1110 & ftprim(3)*wcorr*ecorr+ftprim(4)*wcorr5*ecorr5+
1111 & ftprim(5)*wcorr6*ecorr6+ftprim(3)*wturn4*eello_turn4+
1112 & ftprim(2)*wturn3*eello_turn3+ftprim(5)*wturn6*eturn6+
1113 & ftprim(2)*wel_loc*eel_loc+ftprim(2)*wtor_d*etors_d+
1114 & ftprim(1)*wsccor*esccor
1115 ebis=ftbis(1)*wsc*(evdw+ft(6)*evdw_t)
1116 & +ftbis(1)*wscp*evdw2+
1117 & ftbis(1)*welec*ees
1118 & +ftbis(1)*wvdwpp*evdw
1119 & +ftbis(1)*wtor*etors+
1120 & ftbis(3)*wcorr*ecorr+ftbis(4)*wcorr5*ecorr5+
1121 & ftbis(5)*wcorr6*ecorr6+ftbis(3)*wturn4*eello_turn4+
1122 & ftbis(2)*wturn3*eello_turn3+ftbis(5)*wturn6*eturn6+
1123 & ftbis(2)*wel_loc*eel_loc+ftbis(2)*wtor_d*etors_d+
1124 & ftbis(1)*wsccor*esccor
1126 etot=wsc*(evdw+ft(6)*evdw_t)+wscp*evdw2+ft(1)*welec*ees
1128 & +wang*ebe+ft(1)*wtor*etors+wscloc*escloc
1129 & +wstrain*ehpb+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5
1130 & +ft(5)*wcorr6*ecorr6+ft(3)*wturn4*eello_turn4
1131 & +ft(2)*wturn3*eello_turn3
1132 & +ft(5)*wturn6*eturn6+ft(2)*wel_loc*eel_loc
1133 & +edihcnstr+ft(2)*wtor_d*etors_d+ft(1)*wsccor*esccor
1134 & +wbond*estr+wliptran*eliptran+wtube*Etube
1135 eprim=ftprim(6)*evdw_t+ftprim(1)*welec*ees
1136 & +ftprim(1)*wtor*etors+
1137 & ftprim(3)*wcorr*ecorr+ftprim(4)*wcorr5*ecorr5+
1138 & ftprim(5)*wcorr6*ecorr6+ftprim(3)*wturn4*eello_turn4+
1139 & ftprim(2)*wturn3*eello_turn3+ftprim(5)*wturn6*eturn6+
1140 & ftprim(2)*wel_loc*eel_loc+ftprim(2)*wtor_d*etors_d+
1141 & ftprim(1)*wsccor*esccor
1142 ebis=ftbis(1)*welec*ees+ftbis(1)*wtor*etors+
1143 & ftbis(3)*wcorr*ecorr+ftbis(4)*wcorr5*ecorr5+
1144 & ftbis(5)*wcorr6*ecorr6+ftbis(3)*wturn4*eello_turn4+
1145 & ftbis(2)*wturn3*eello_turn3+ftbis(5)*wturn6*eturn6+
1146 & ftbis(2)*wel_loc*eel_loc+ftbis(2)*wtor_d*etors_d+
1147 & ftbis(1)*wsccor*esccor
1150 if (shield_mode.gt.0) then
1151 etot=ft(1)*wsc*(evdw+ft(6)*evdw_t)+ft(1)*wscp*evdw2
1152 & +ft(1)*welec*(ees+evdw1)
1153 & +wang*ebe+ft(1)*wtor*etors+wscloc*escloc
1154 & +wstrain*ehpb+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5
1155 & +ft(5)*wcorr6*ecorr6+ft(3)*wturn4*eello_turn4
1156 & +ft(2)*wturn3*eello_turn3
1157 & +ft(5)*wturn6*eturn6+ft(2)*wel_loc*eel_loc+edihcnstr
1158 & +ft(2)*wtor_d*etors_d+ft(1)*wsccor*esccor
1159 & +wbond*estr+wliptran*eliptran+wtube*Etube
1160 eprim=ftprim(1)*(evdw+ft(6)*evdw_t)
1161 & +ftprim(1)*welec*(ees+evdw1)
1162 & +ftprim(1)*wtor*etors+
1163 & ftprim(1)*wscp*evdw2+
1164 & ftprim(3)*wcorr*ecorr+ftprim(4)*wcorr5*ecorr5+
1165 & ftprim(5)*wcorr6*ecorr6+ftprim(3)*wturn4*eello_turn4+
1166 & ftprim(2)*wturn3*eello_turn3+ftprim(5)*wturn6*eturn6+
1167 & ftprim(2)*wel_loc*eel_loc+ftprim(2)*wtor_d*etors_d+
1168 & ftprim(1)*wsccor*esccor
1169 ebis= ftbis(1)*(evdw+ft(6)*evdw_t)
1170 & +ftbis(1)*wscp*evdw2
1171 & +ftbis(1)*welec*(ees+evdw1)+ftbis(1)*wtor*etors+
1172 & ftbis(3)*wcorr*ecorr+ftbis(4)*wcorr5*ecorr5+
1173 & ftbis(5)*wcorr6*ecorr6+ftbis(3)*wturn4*eello_turn4+
1174 & ftbis(2)*wturn3*eello_turn3+ftbis(5)*wturn6*eturn6+
1175 & ftbis(2)*wel_loc*eel_loc+ftbis(2)*wtor_d*etors_d+
1176 & ftbis(1)*wsccor*esccor
1178 etot=wsc*(evdw+ft(6)*evdw_t)+wscp*evdw2
1179 & +ft(1)*welec*(ees+evdw1)
1180 & +wang*ebe+ft(1)*wtor*etors+wscloc*escloc
1181 & +wstrain*ehpb+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5
1182 & +ft(5)*wcorr6*ecorr6+ft(3)*wturn4*eello_turn4
1183 & +ft(2)*wturn3*eello_turn3
1184 & +ft(5)*wturn6*eturn6+ft(2)*wel_loc*eel_loc+edihcnstr
1185 & +ft(2)*wtor_d*etors_d+ft(1)*wsccor*esccor
1186 & +wbond*estr+wliptran*eliptran+wtube*Etube
1187 eprim=ftprim(6)*evdw_t+ftprim(1)*welec*(ees+evdw1)
1188 & +ftprim(1)*wtor*etors+
1189 & ftprim(3)*wcorr*ecorr+ftprim(4)*wcorr5*ecorr5+
1190 & ftprim(5)*wcorr6*ecorr6+ftprim(3)*wturn4*eello_turn4+
1191 & ftprim(2)*wturn3*eello_turn3+ftprim(5)*wturn6*eturn6+
1192 & ftprim(2)*wel_loc*eel_loc+ftprim(2)*wtor_d*etors_d+
1193 & ftprim(1)*wsccor*esccor
1194 ebis=ftbis(1)*welec*(ees+evdw1)+ftbis(1)*wtor*etors+
1195 & ftbis(3)*wcorr*ecorr+ftbis(4)*wcorr5*ecorr5+
1196 & ftbis(5)*wcorr6*ecorr6+ftbis(3)*wturn4*eello_turn4+
1197 & ftbis(2)*wturn3*eello_turn3+ftbis(5)*wturn6*eturn6+
1198 & ftbis(2)*wel_loc*eel_loc+ftbis(2)*wtor_d*etors_d+
1199 & ftbis(1)*wsccor*esccor
1205 weight=dexp(-betaT*(etot-potEmin)+entfac(t))
1207 write (iout,*) "iparm",iparm," t",t," betaT",betaT,
1208 & " etot",etot," entfac",entfac(t),
1209 & " weight",weight," ebis",ebis
1211 etot=etot-temper*eprim
1213 sumW_p(k,iparm)=sumW_p(k,iparm)+weight
1214 sumE_p(k,iparm)=sumE_p(k,iparm)+etot*weight
1215 sumEbis_p(k,iparm)=sumEbis_p(k,iparm)+ebis*weight
1216 sumEsq_p(k,iparm)=sumEsq_p(k,iparm)+etot**2*weight
1218 sumQ_p(j,k,iparm)=sumQ_p(j,k,iparm)+q(j,t)*weight
1219 sumQsq_p(j,k,iparm)=sumQsq_p(j,k,iparm)+q(j,t)**2*weight
1220 sumEQ_p(j,k,iparm)=sumEQ_p(j,k,iparm)
1221 & +etot*q(j,t)*weight
1224 sumW(k,iparm)=sumW(k,iparm)+weight
1225 sumE(k,iparm)=sumE(k,iparm)+etot*weight
1226 sumEbis(k,iparm)=sumEbis(k,iparm)+ebis*weight
1227 sumEsq(k,iparm)=sumEsq(k,iparm)+etot**2*weight
1229 sumQ(j,k,iparm)=sumQ(j,k,iparm)+q(j,t)*weight
1230 sumQsq(j,k,iparm)=sumQsq(j,k,iparm)+q(j,t)**2*weight
1231 sumEQ(j,k,iparm)=sumEQ(j,k,iparm)
1232 & +etot*q(j,t)*weight
1236 indE = aint(potE(t,iparm)-aint(potEmin))
1237 if (indE.ge.0 .and. indE.le.maxinde) then
1238 if (indE.gt.upindE_p) upindE_p=indE
1239 histE_p(indE)=histE_p(indE)+dexp(-entfac(t))
1243 potEmin=potEmin_all(ib,iparm)
1244 expfac=dexp(-beta_h(ib,iparm)*(etot-potEmin)+entfac(t))
1245 hfin_p(ind,ib)=hfin_p(ind,ib)+
1246 & dexp(-beta_h(ib,iparm)*(etot-potEmin)+entfac(t))
1248 indrgy=dint((q(nQ+2,t)-rgymin)/deltrgy)
1249 indrms=dint((q(nQ+1,t)-rmsmin)/deltrms)
1250 hrmsrgy_p(indrgy,indrms,ib)=
1251 & hrmsrgy_p(indrgy,indrms,ib)+expfac
1256 potEmin=potEmin_all(ib,iparm)
1257 expfac=dexp(-beta_h(ib,iparm)*(etot-potEmin)+entfac(t))
1258 hfin(ind,ib)=hfin(ind,ib)+
1259 & dexp(-beta_h(ib,iparm)*(etot-potEmin)+entfac(t))
1261 indrgy=dint((q(nQ+2,t)-rgymin)/deltrgy)
1262 indrms=dint((q(nQ+1,t)-rmsmin)/deltrms)
1263 hrmsrgy(indrgy,indrms,ib)=
1264 & hrmsrgy(indrgy,indrms,ib)+expfac
1270 if (histout) call MPI_Reduce(hfin_p(0,ib),hfin(0,ib),nbin,
1271 & MPI_DOUBLE_PRECISION,MPI_SUM,Master,WHAM_COMM,IERROR)
1273 call MPI_Reduce(hrmsrgy_p(0,0,ib),hrmsrgy(0,0,ib),
1274 & (MaxBinRgy+1)*(nbin_rms+1),MPI_DOUBLE_PRECISION,MPI_SUM,Master,
1278 call MPI_Reduce(upindE_p,upindE,1,
1279 & MPI_INTEGER,MPI_MAX,Master,WHAM_COMM,IERROR)
1280 call MPI_Reduce(histE_p(0),histE(0),maxindE,
1281 & MPI_DOUBLE_PRECISION,MPI_SUM,Master,WHAM_COMM,IERROR)
1283 if (me1.eq.master) then
1287 write (iout,'(6x,$)')
1288 write (iout,'(f20.2,$)') (1.0d0/(1.987D-3*beta_h(ib,iparm)),
1292 write (iout,'(/a)') 'Final histograms'
1294 if (nslice.eq.1) then
1295 if (separate_parset) then
1296 write(licz3,"(bz,i3.3)") myparm
1297 histname=prefix(:ilen(prefix))//'_par'//licz3//'.hist'
1299 histname=prefix(:ilen(prefix))//'.hist'
1302 if (separate_parset) then
1303 write(licz3,"(bz,i3.3)") myparm
1304 histname=prefix(:ilen(prefix))//'_par'//licz3//
1305 & '_slice_'//licz2//'.hist'
1307 histname=prefix(:ilen(prefix))//'_slice_'//licz2//'.hist'
1310 #if defined(AIX) || defined(PGI)
1311 open (ihist,file=histname,position='append')
1313 open (ihist,file=histname,access='append')
1321 sumH=sumH+hfin(t,ib)
1323 if (sumH.gt.0.0d0) then
1325 jj = mod(liczba,nbin1)
1327 write (iout,'(f6.3,$)') dmin+(jj+0.5d0)*delta
1329 & write (ihist,'(f6.3,$)') dmin+(jj+0.5d0)*delta
1332 write (iout,'(e20.10,$)') hfin(t,ib)
1333 if (histfile) write (ihist,'(e20.10,$)') hfin(t,ib)
1335 write (iout,'(i5)') iparm
1336 if (histfile) write (ihist,'(i5)') iparm
1343 if (nslice.eq.1) then
1344 if (separate_parset) then
1345 write(licz3,"(bz,i3.3)") myparm
1346 histname=prefix(:ilen(prefix))//"_par"//licz3//'.ent'
1348 histname=prefix(:ilen(prefix))//'.ent'
1351 if (separate_parset) then
1352 write(licz3,"(bz,i3.3)") myparm
1353 histname=prefix(:ilen(prefix))//'par_'//licz3//
1354 & '_slice_'//licz2//'.ent'
1356 histname=prefix(:ilen(prefix))//'_slice_'//licz2//'.ent'
1359 #if defined(AIX) || defined(PGI)
1360 open (ihist,file=histname,position='append')
1362 open (ihist,file=histname,access='append')
1364 write (ihist,'(a)') "# Microcanonical entropy"
1366 write (ihist,'(f8.0,$)') dint(potEmin)+i
1367 if (histE(i).gt.0.0e0) then
1368 write (ihist,'(f15.5,$)') dlog(histE(i))
1370 write (ihist,'(f15.5,$)') 0.0d0
1376 write (iout,*) "Microcanonical entropy"
1378 write (iout,'(f8.0,$)') dint(potEmin)+i
1379 if (histE(i).gt.0.0e0) then
1380 write (iout,'(f15.5,$)') dlog(histE(i))
1382 write (iout,'(f15.5,$)') 0.0d0
1387 if (nslice.eq.1) then
1388 if (separate_parset) then
1389 write(licz3,"(bz,i3.3)") myparm
1390 histname=prefix(:ilen(prefix))//'_par'//licz3//'.rmsrgy'
1392 histname=prefix(:ilen(prefix))//'.rmsrgy'
1395 if (separate_parset) then
1396 write(licz3,"(bz,i3.3)") myparm
1397 histname=prefix(:ilen(prefix))//'_par'//licz3//
1398 & '_slice_'//licz2//'.rmsrgy'
1400 histname=prefix(:ilen(prefix))//'_slice_'//licz2//'.rmsrgy'
1403 #if defined(AIX) || defined(PGI)
1404 open (ihist,file=histname,position='append')
1406 open (ihist,file=histname,access='append')
1410 write(ihist,'(2f8.2,$)')
1411 & rgymin+deltrgy*j,rmsmin+deltrms*i
1413 if (hrmsrgy(j,i,ib).gt.0.0d0) then
1414 write(ihist,'(e14.5,$)')
1415 & -dlog(hrmsrgy(j,i,ib))/beta_h(ib,iparm)
1418 write(ihist,'(e14.5,$)') 1.0d6
1421 write (ihist,'(i2)') iparm
1429 call MPI_Reduce(hfin_ent_p(0),hfin_ent(0),nbin,
1430 & MPI_DOUBLE_PRECISION,MPI_SUM,Master,WHAM_COMM,IERROR)
1431 call MPI_Reduce(sumW_p(0,1),sumW(0,1),(nGridT+1)*nParmSet,
1432 & MPI_DOUBLE_PRECISION,MPI_SUM,Master,WHAM_COMM,IERROR)
1433 call MPI_Reduce(sumE_p(0,1),sumE(0,1),(nGridT+1)*nParmSet,
1434 & MPI_DOUBLE_PRECISION,MPI_SUM,Master,WHAM_COMM,IERROR)
1435 call MPI_Reduce(sumEbis_p(0,1),sumEbis(0,1),(nGridT+1)*nParmSet,
1436 & MPI_DOUBLE_PRECISION,MPI_SUM,Master,WHAM_COMM,IERROR)
1437 call MPI_Reduce(sumEsq_p(0,1),sumEsq(0,1),(nGridT+1)*nParmSet,
1438 & MPI_DOUBLE_PRECISION,MPI_SUM,Master,WHAM_COMM,IERROR)
1439 call MPI_Reduce(sumQ_p(1,0,1),sumQ(1,0,1),
1440 & MaxQ1*(nGridT+1)*nParmSet,MPI_DOUBLE_PRECISION,MPI_SUM,Master,
1442 call MPI_Reduce(sumQsq_p(1,0,1),sumQsq(1,0,1),
1443 & MaxQ1*(nGridT+1)*nParmSet,MPI_DOUBLE_PRECISION,MPI_SUM,Master,
1445 call MPI_Reduce(sumEQ_p(1,0,1),sumEQ(1,0,1),
1446 & MaxQ1*(nGridT+1)*nParmSet,MPI_DOUBLE_PRECISION,MPI_SUM,Master,
1448 if (me.eq.master) then
1450 write (iout,'(/a)') 'Thermal characteristics of folding'
1451 if (nslice.eq.1) then
1454 nazwa=prefix(:ilen(prefix))//"_slice_"//licz2
1457 if (nparmset.eq.1 .and. .not.separate_parset) then
1458 nazwa=nazwa(:iln)//".thermal"
1459 else if (nparmset.eq.1 .and. separate_parset) then
1460 write(licz3,"(bz,i3.3)") myparm
1461 nazwa=nazwa(:iln)//"_par_"//licz3//".thermal"
1464 if (nparmset.gt.1) then
1465 write(licz3,"(bz,i3.3)") iparm
1466 nazwa=nazwa(:iln)//"_par_"//licz3//".thermal"
1469 if (separate_parset) then
1470 write (iout,'(a,i3)') "Parameter set",myparm
1472 write (iout,'(a,i3)') "Parameter set",iparm
1475 betaT=1.0d0/(1.987D-3*(startGridT+i*delta_T))
1476 if (betaT.ge.beta_h(1,iparm)) then
1477 potEmin=potEmin_all(1,iparm)
1478 else if (betaT.lt.beta_h(nT_h(iparm),iparm)) then
1479 potEmin=potEmin_all(nT_h(iparm),iparm)
1481 do l=1,nT_h(iparm)-1
1482 if (betaT.le.beta_h(l,iparm) .and.
1483 & betaT.gt.beta_h(l+1,iparm)) then
1484 potEmin=potEmin_all(l,iparm)
1490 sumE(i,iparm)=sumE(i,iparm)/sumW(i,iparm)
1491 sumEbis(i,iparm)=(startGridT+i*delta_T)*sumEbis(i,iparm)/
1493 sumEsq(i,iparm)=(sumEsq(i,iparm)/sumW(i,iparm)
1494 & -sumE(i,iparm)**2)/(1.987D-3*(startGridT+i*delta_T)**2)
1496 sumQ(j,i,iparm)=sumQ(j,i,iparm)/sumW(i,iparm)
1497 sumQsq(j,i,iparm)=sumQsq(j,i,iparm)/sumW(i,iparm)
1498 & -sumQ(j,i,iparm)**2
1499 sumEQ(j,i,iparm)=sumEQ(j,i,iparm)/sumW(i,iparm)
1500 & -sumQ(j,i,iparm)*sumE(i,iparm)
1502 sumW(i,iparm)=-dlog(sumW(i,iparm))*(1.987D-3*
1503 & (startGridT+i*delta_T))+potEmin
1504 write (iout,'(f7.1,2f15.5,$)') startGridT+i*delta_T,
1505 & sumW(i,iparm),sumE(i,iparm)
1506 write (iout,'(f10.5,$)') (sumQ(j,i,iparm),j=1,nQ+2)
1507 write (iout,'(e15.5,$)') sumEsq(i,iparm)-sumEbis(i,iparm),
1508 & (sumQsq(j,i,iparm),j=1,nQ+2),(sumEQ(j,i,iparm),j=1,nQ+2)
1510 write (34,'(f7.1,2f15.5,$)') startGridT+i*delta_T,
1511 & sumW(i,iparm),sumE(i,iparm)
1512 write (34,'(f10.5,$)') (sumQ(j,i,iparm),j=1,nQ+2)
1513 write (34,'(e15.5,$)') sumEsq(i,iparm)-sumEbis(i,iparm),
1514 & (sumQsq(j,i,iparm),j=1,nQ+2),(sumEQ(j,i,iparm),j=1,nQ+2)
1522 if (hfin_ent(t).gt.0.0d0) then
1524 jj = mod(liczba,nbin1)
1525 write (iout,'(f6.3,e20.10," ent")') dmin+(jj+0.5d0)*delta,
1527 if (histfile) write (ihist,'(f6.3,e20.10," ent")')
1528 & dmin+(jj+0.5d0)*delta,
1532 if (histfile) close(ihist)
1536 ! Write data for zscore
1537 if (nslice.eq.1) then
1538 zscname=prefix(:ilen(prefix))//".zsc"
1540 zscname=prefix(:ilen(prefix))//"_slice_"//licz2//".zsc"
1542 #if defined(AIX) || defined(PGI)
1543 open (izsc,file=prefix(:ilen(prefix))//'.zsc',position='append')
1545 open (izsc,file=prefix(:ilen(prefix))//'.zsc',access='append')
1547 write (izsc,'("NQ=",i1," NPARM=",i1)') nQ,nParmSet
1549 write (izsc,'("NT=",i1)') nT_h(iparm)
1551 write (izsc,'("TEMP=",f6.1," NR=",i2," SNK=",$)')
1552 & 1.0d0/(beta_h(ib,iparm)*1.987D-3),nR(ib,iparm)
1553 jj = min0(nR(ib,iparm),7)
1554 write (izsc,'(i8,$)') (snk(i,ib,iparm,islice),i=1,jj)
1555 write (izsc,'(a1,$)') (" ",i=22+8*jj+1,79)
1556 write (izsc,'("&")')
1557 if (nR(ib,iparm).gt.7) then
1558 do ii=8,nR(ib,iparm),9
1559 jj = min0(nR(ib,iparm),ii+8)
1560 write (izsc,'(i8,$)') (snk(i,ib,iparm,islice),i=ii,jj)
1561 write (izsc,'(a1,$') (" ",i=(jj-ii+1)*8+1,79)
1562 write (izsc,'("&")')
1565 write (izsc,'("FI=",$)')
1566 jj=min0(nR(ib,iparm),7)
1567 write (izsc,'(f10.5,$)') (fi(i,ib,iparm),i=1,jj)
1568 write (izsc,'(a1,$)') (" ",i=3+10*jj+1,79)
1569 write (izsc,'("&")')
1570 if (nR(ib,iparm).gt.7) then
1571 do ii=8,nR(ib,iparm),9
1572 jj = min0(nR(ib,iparm),ii+8)
1573 write (izsc,'(f10.5,$)') (fi(i,ib,iparm),i=ii,jj)
1574 if (jj.eq.nR(ib,iparm)) then
1577 write (izsc,'(a1,$)') (" ",i=10*(jj-ii+1)+1,79)
1578 write (izsc,'(t80,"&")')
1583 write (izsc,'("KH=",$)')
1584 write (izsc,'(f7.2,$)') (Kh(j,i,ib,iparm),j=1,nQ)
1585 write (izsc,'(" Q0=",$)')
1586 write (izsc,'(f7.5,$)') (q0(j,i,ib,iparm),j=1,nQ)