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"
17 integer MaxBinRms,MaxBinRgy
18 parameter (MaxBinRms=100,MaxBinRgy=100)
20 c parameter (MaxHdim=200000)
22 parameter (MaxHdim=200)
24 parameter (maxinde=200)
26 parameter (MaxHdim=100)
28 parameter (maxinde=100)
29 >>>>>>> e183793... Added src_MD-M-newcorr (Adasko's source) and src-NEWSC of WHAM (with Momo's SCSC potentials)
33 integer ierror,errcode,status(MPI_STATUS_SIZE)
35 include "COMMON.CONTROL"
36 include "COMMON.IOUNITS"
38 include "COMMON.ENERGIES"
39 include "COMMON.FFIELD"
40 include "COMMON.SBRIDGE"
42 include "COMMON.ENEPS"
43 integer MaxPoint,MaxPointProc
44 parameter (MaxPoint=MaxStr,
45 & MaxPointProc=MaxStr_Proc)
46 double precision finorm_max,potfac,entmin,entmax,expfac,vf
47 double precision entfac_min,entfac_min_t
48 parameter (finorm_max=1.0d0)
50 integer i,ii,j,jj,k,kk,l,m,ind,iter,t,tmax,ient,ientmax,iln
51 integer start,end,iharm,ib,iib,nbin1,nbin,nbin_rms,nbin_rgy,
52 & nbin_rmsrgy,liczba,iparm,nFi,indrgy,indrms
53 integer htot(0:MaxHdim),histent(0:2000)
54 double precision v(MaxPointProc,MaxR,MaxT_h,Max_Parm)
55 double precision energia(0:max_ene)
57 integer tmax_t,upindE_p
58 double precision fi_p(MaxR,MaxT_h,Max_Parm),
59 & fi_p_min(MaxR,MaxT_h,Max_Parm)
60 double precision sumW_p(0:Max_GridT,Max_Parm),
61 & sumE_p(0:Max_GridT,Max_Parm),sumEsq_p(0:Max_GridT,Max_Parm),
62 & sumQ_p(MaxQ1,0:Max_GridT,Max_Parm),
63 & sumQsq_p(MaxQ1,0:Max_GridT,Max_Parm),
64 & sumEQ_p(MaxQ1,0:Max_GridT,Max_Parm),
65 & sumEprim_p(MaxQ1,0:Max_GridT,Max_Parm),
66 & sumEbis_p(0:Max_GridT,Max_Parm)
67 double precision hfin_p(0:MaxHdim,maxT_h),
68 & hfin_ent_p(0:MaxHdim),histE_p(0:maxindE),sumH,
69 & hrmsrgy_p(0:MaxBinRgy,0:MaxBinRms,maxT_h)
70 double precision rgymin_t,rmsmin_t,rgymax_t,rmsmax_t
71 double precision potEmin_t_all(maxT_h,Max_Parm),entmin_p,entmax_p
72 integer histent_p(0:2000)
73 logical lprint /.true./
75 double precision rgymin,rmsmin,rgymax,rmsmax
76 double precision sumW(0:NGridT,Max_Parm),sumE(0:NGridT,Max_Parm),
77 & sumEsq(0:NGridT,Max_Parm),sumQ(MaxQ1,0:NGridT,Max_Parm),
78 & sumQsq(MaxQ1,0:NGridT,Max_Parm),sumEQ(MaxQ1,0:NGridT,Max_Parm),
79 & sumEprim(0:NGridT,Max_Parm),sumEbis(0:NGridT,Max_Parm),betaT,
81 double precision fi(MaxR,maxT_h,Max_Parm),
82 & fi_min(MaxR,maxT_h,Max_Parm),
83 & dd,dd1,dd2,hh,dmin,denom,finorm,avefi,pom,
84 & hfin(0:MaxHdim,maxT_h),histE(0:maxindE),
85 & hrmsrgy(0:MaxBinRgy,0:MaxBinRms,maxT_h),
86 & potEmin_all(maxT_h,Max_Parm),potEmin,potEmin_min,ent,
87 & hfin_ent(0:MaxHdim),vmax,aux
88 double precision fT(6),fTprim(6),fTbis(6),quot,quotl1,quotl,kfacl,
89 & eprim,ebis,temper,kfac/2.4d0/,T0/300.0d0/,
90 & eplus,eminus,logfac,tanhT,tt
91 double precision etot,evdw,evdw_t,evdw2,ees,evdw1,ebe,etors,
92 & escloc,ehpb,ecorr,ecorr5,ecorr6,eello_turn4,eello_turn3,
93 & eturn6,eel_loc,edihcnstr,etors_d,estr,evdw2_14,esccor
95 integer ind_point(maxpoint),upindE,indE
106 write (iout,*) "Enter WHAM_calc"
108 >>>>>>> e183793... Added src_MD-M-newcorr (Adasko's source) and src-NEWSC of WHAM (with Momo's SCSC potentials)
109 write(licz2,'(bz,i2.2)') islice
111 write (iout,'(//80(1h-)/"Solving WHAM equations for slice",
112 & i2/80(1h-)//)') islice
113 write (iout,*) "delta",delta," nbin1",nbin1
114 write (iout,*) "MaxN",MaxN," MaxQ",MaxQ," MaHdim",MaxHdim
120 potEmin_all(j,i)=1.0d10
135 if (q(nQ+1,i).lt.rmsmin) rmsmin=q(nQ+1,i)
136 if (q(nQ+1,i).gt.rmsmax) rmsmax=q(nQ+1,i)
137 if (q(nQ+2,i).lt.rgymin) rgymin=q(nQ+2,i)
138 if (q(nQ+2,i).gt.rgymax) rgymax=q(nQ+2,i)
141 ind=(q(j,i)-dmin+1.0d-8)/delta
143 ind_point(i)=ind_point(i)+ind
145 ind_point(i)=ind_point(i)+nbin1**(j-1)*ind
147 if (ind_point(i).lt.0 .or. ind_point(i).gt.MaxHdim) then
148 write (iout,*) "Error - index exceeds range for point",i,
149 & " q=",q(j,i)," ind",ind_point(i)
151 write (iout,*) "Processor",me1
153 call MPI_Abort(MPI_COMM_WORLD, Ierror, Errcode )
158 if (ind_point(i).gt.tmax) tmax=ind_point(i)
159 htot(ind_point(i))=htot(ind_point(i))+1
161 write (iout,*) "i",i,"q",(q(j,i),j=1,nQ)," ind",ind_point(i),
162 & " htot",htot(ind_point(i))
169 write (iout,'(a)') "Numbers of counts in Q bins"
171 if (htot(t).gt.0) then
172 write (iout,'(i15,$)') t
175 jj = mod(liczba,nbin1)
177 write (iout,'(i5,$)') jj
179 write (iout,'(i8)') htot(t)
183 write (iout,'(a,i3)') "Number of data points for parameter set",
185 write (iout,'(i7,$)') ((snk(m,ib,iparm,islice),m=1,nr(ib,iparm)),
187 write (iout,'(i8)') stot(islice)
193 call MPI_AllReduce(tmax,tmax_t,1,MPI_INTEGER,MPI_MAX,
196 call MPI_AllReduce(rmsmin,rmsmin_t,1,MPI_DOUBLE_PRECISION,
197 & MPI_MIN,WHAM_COMM,IERROR)
198 call MPI_AllReduce(rmsmax,rmsmax_t,1,MPI_DOUBLE_PRECISION,
199 & MPI_MAX,WHAM_COMM,IERROR)
200 call MPI_AllReduce(rgymin,rgymin_t,1,MPI_DOUBLE_PRECISION,
201 & MPI_MIN,WHAM_COMM,IERROR)
202 call MPI_AllReduce(rgymax,rgymax_t,1,MPI_DOUBLE_PRECISION,
203 & MPI_MAX,WHAM_COMM,IERROR)
209 rmsmin=deltrms*dint(rmsmin/deltrms)
210 rmsmax=deltrms*dint(rmsmax/deltrms)
211 rgymin=deltrms*dint(rgymin/deltrgy)
212 rgymax=deltrms*dint(rgymax/deltrgy)
213 nbin_rms=(rmsmax-rmsmin)/deltrms
214 nbin_rgy=(rgymax-rgymin)/deltrgy
215 write (iout,*) "rmsmin",rmsmin," rmsmax",rmsmax," rgymin",rgymin,
216 & " rgymax",rgymax," nbin_rms",nbin_rms," nbin_rgy",nbin_rgy
223 write (iout,*) "nFi",nFi
224 ! Compute the Boltzmann factor corresponing to restrain potentials in different
231 c write (9,'(3i5,f10.5)') i,(iparm,potE(i,iparm),iparm=1,nParmSet)
234 write (iout,'(2i5,21f8.2)') i,iparm,
235 & (enetb(k,i,iparm),k=1,21)
237 call restore_parm(iparm)
239 write (iout,*) wsc,wscp,welec,wvdwpp,wang,wtor,wscloc,
240 & wcorr,wcorr5,wcorr6,wturn4,wturn3,wturn6,wel_loc,
241 & wtor_d,wsccor,wbond
244 if (rescale_mode.eq.1) then
245 quot=1.0d0/(beta_h(ib,iparm)*1.987D-3*T0)
252 fT(l)=kfacl/(kfacl-1.0d0+quotl)
255 tt = 1.0d0/(beta_h(ib,iparm)*1.987D-3)
256 ft(6)=(320.0d0+80.0d0*dtanh((tt-320.0d0)/80.0d0))/320.0d0
258 ft(6)=1.0d0/(beta_h(ib,iparm)*1.987D-3*T0)
262 else if (rescale_mode.eq.2) then
263 quot=1.0d0/(T0*beta_h(ib,iparm)*1.987D-3)
267 fT(l)=1.12692801104297249644d0/
268 & dlog(dexp(quotl)+dexp(-quotl))
271 tt = 1.0d0/(beta_h(ib,iparm)*1.987D-3)
272 ft(6)=(320.0d0+80.0d0*dtanh((tt-320.0d0)/80.0d0))/320.0d0
274 ft(6)=1.0d0/(beta_h(ib,iparm)*1.987D-3*T0)
278 c write (iout,*) 1.0d0/(beta_h(ib,iparm)*1.987D-3),ft
279 else if (rescale_mode.eq.0) then
284 write (iout,*) "Error in WHAM_CALC: wrong RESCALE_MODE",
289 evdw=enetb(1,i,iparm)
290 evdw_t=enetb(21,i,iparm)
292 evdw2_14=enetb(17,i,iparm)
293 evdw2=enetb(2,i,iparm)+evdw2_14
295 evdw2=enetb(2,i,iparm)
300 evdw1=enetb(16,i,iparm)
305 ecorr=enetb(4,i,iparm)
306 ecorr5=enetb(5,i,iparm)
307 ecorr6=enetb(6,i,iparm)
308 eel_loc=enetb(7,i,iparm)
309 eello_turn3=enetb(8,i,iparm)
310 eello_turn4=enetb(9,i,iparm)
311 eturn6=enetb(10,i,iparm)
312 ebe=enetb(11,i,iparm)
313 escloc=enetb(12,i,iparm)
314 etors=enetb(13,i,iparm)
315 etors_d=enetb(14,i,iparm)
316 ehpb=enetb(15,i,iparm)
317 estr=enetb(18,i,iparm)
318 esccor=enetb(19,i,iparm)
319 edihcnstr=enetb(20,i,iparm)
321 write (iout,'(3i5,6f5.2,14f12.3)') i,ib,iparm,(ft(l),l=1,6),
322 & evdw+evdw_t,evdw2,ees,evdw1,ecorr,eel_loc,estr,ebe,escloc,
323 & etors,etors_d,eello_turn3,eello_turn4,esccor
327 etot=wsc*(evdw+ft(6)*evdw_t)+wscp*evdw2+ft(1)*welec*ees
329 & +wang*ebe+ft(1)*wtor*etors+wscloc*escloc
330 & +wstrain*ehpb+nss*ebr+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5
331 & +ft(5)*wcorr6*ecorr6+ft(3)*wturn4*eello_turn4
332 & +ft(2)*wturn3*eello_turn3
333 & +ft(5)*wturn6*eturn6+ft(2)*wel_loc*eel_loc
334 & +edihcnstr+ft(2)*wtor_d*etors_d+ft(1)*wsccor*esccor
337 etot=wsc*(evdw+ft(6)*evdw_t)+wscp*evdw2
338 & +ft(1)*welec*(ees+evdw1)
339 & +wang*ebe+ft(1)*wtor*etors+wscloc*escloc
340 & +wstrain*ehpb+nss*ebr+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5
341 & +ft(5)*wcorr6*ecorr6+ft(3)*wturn4*eello_turn4
342 & +ft(2)*wturn3*eello_turn3
343 & +ft(5)*wturn6*eturn6+ft(2)*wel_loc*eel_loc+edihcnstr
344 & +ft(2)*wtor_d*etors_d+ft(1)*wsccor*esccor
348 write (iout,*) i,iparm,1.0d0/(beta_h(ib,iparm)*1.987D-3),
352 if (iparm.eq.1 .and. ib.eq.1) then
353 write (iout,*)"Conformation",i
356 energia(k)=enetb(k,i,iparm)
358 call enerprint(energia(0),fT)
365 Econstr=Econstr+Kh(j,kk,ib,iparm)
366 & *(dd-q0(j,kk,ib,iparm))**2
369 & -beta_h(ib,iparm)*(etot+Econstr)
371 write (iout,'(4i5,4e15.5)') i,kk,ib,iparm,
372 & etot,v(i,kk,ib,iparm)
378 ! Simple iteration to calculate free energies corresponding to all simulation
382 ! Compute new free-energy values corresponding to the righ-hand side of the
383 ! equation and their derivatives.
384 write (iout,*) "------------------------fi"
395 vf=v(t,l,k,i)+f(l,k,i)
396 if (vf.gt.vmax) vmax=vf
404 aux=f(l,k,i)+v(t,l,k,i)-vmax
406 & denom=denom+snk(l,k,i,islice)*dexp(aux)
410 entfac(t)=-dlog(denom)-vmax
411 if (entfac(t).lt.entfac_min) entfac_min=entfac(t)
413 write (iout,*) t,"vmax",vmax," denom",denom,"entfac",entfac(t)
417 c write (iout,*) "entfac_min before AllReduce",entfac_min
418 c call MPI_AllReduce(entfac_min,entfac_min_t,1,
419 c & MPI_DOUBLE_PRECISION,MPI_MAX,WHAM_COMM,IERROR)
420 c entfac_min=entfac_min_t
421 c write (iout,*) "entfac_min after AllReduce",entfac_min
425 c entfac(t)=entfac(t)-entfac_min
428 c do t=1,ntot(islice)
429 c entfac(t)=entfac(t)-entfac_min
434 do ii=1,nR(iib,iparm)
436 fi_p_min(ii,iib,iparm)=-1.0d10
438 aux=v(t,ii,iib,iparm)+entfac(t)
439 if (aux.gt.fi_p_min(ii,iib,iparm))
440 & fi_p_min(ii,iib,iparm)=aux
444 aux=v(t,ii,iib,iparm)+entfac(t)
445 if (aux.gt.fi_min(ii,iib,iparm))
446 & fi_min(ii,iib,iparm)=aux
454 write (iout,*) "fi_min before AllReduce"
457 write (iout,*) (i,j,k,fi_p_min(k,j,i),k=1,nR(j,i))
461 call MPI_AllReduce(fi_p_min,fi_min,MaxR*MaxT_h*nParmSet,
462 & MPI_DOUBLE_PRECISION,MPI_MAX,WHAM_COMM,IERROR)
464 write (iout,*) "fi_min after AllReduce"
467 write (iout,*) (i,j,k,fi_min(k,j,i),k=1,nR(j,i))
474 do ii=1,nR(iib,iparm)
476 fi_p(ii,iib,iparm)=0.0d0
478 fi_p(ii,iib,iparm)=fi_p(ii,iib,iparm)
479 & +dexp(v(t,ii,iib,iparm)+entfac(t)-fi_min(ii,iib,iparm))
481 write (iout,'(4i5,4e15.5)') t,ii,iib,iparm,
482 & v(t,ii,iib,iparm),entfac(t),fi_min(ii,iib,iparm),
487 fi(ii,iib,iparm)=0.0d0
489 fi(ii,iib,iparm)=fi(ii,iib,iparm)
490 & +dexp(v(t,ii,iib,iparm)+entfac(t)-fi_min(ii,iib,iparm))
499 write (iout,*) "fi before MPI_Reduce me",me,' master',master
502 write (iout,*) "iparm",iparm," ib",ib
503 write (iout,*) "beta=",beta_h(ib,iparm)
504 write (iout,'(8e15.5)') (fi_p(i,ib,iparm),i=1,nR(ib,iparm))
509 write (iout,*) "REDUCE size",maxR,MaxT_h,nParmSet,
510 & maxR*MaxT_h*nParmSet
511 write (iout,*) "MPI_COMM_WORLD",MPI_COMM_WORLD,
512 & " WHAM_COMM",WHAM_COMM
514 call MPI_Reduce(fi_p(1,1,1),fi(1,1,1),maxR*MaxT_h*nParmSet,
515 & MPI_DOUBLE_PRECISION,
516 & MPI_SUM,Master,WHAM_COMM,IERROR)
518 write (iout,*) "fi after MPI_Reduce nparmset",nparmset
520 write (iout,*) "iparm",iparm
522 write (iout,*) "beta=",beta_h(ib,iparm)
523 write (iout,'(8e15.5)') (fi(i,ib,iparm),i=1,nR(ib,iparm))
527 if (me1.eq.Master) then
533 fi(i,ib,iparm)=-dlog(fi(i,ib,iparm))-fi_min(i,ib,iparm)
534 avefi=avefi+fi(i,ib,iparm)
540 write (iout,*) "Parameter set",iparm
542 write (iout,*) "beta=",beta_h(ib,iparm)
544 fi(i,ib,iparm)=fi(i,ib,iparm)-avefi
546 write (iout,'(8f10.5)') (fi(i,ib,iparm),i=1,nR(ib,iparm))
547 write (iout,'(8f10.5)') (f(i,ib,iparm),i=1,nR(ib,iparm))
551 ! Compute the norm of free-energy increments.
556 finorm=finorm+dabs(fi(i,ib,iparm)-f(i,ib,iparm))
557 f(i,ib,iparm)=fi(i,ib,iparm)
562 write (iout,*) 'Iteration',iter,' finorm',finorm
566 call MPI_Bcast(f(1,1,1),MaxR*MaxT_h*nParmSet,
567 & MPI_DOUBLE_PRECISION,Master,
569 call MPI_Bcast(finorm,1,MPI_DOUBLE_PRECISION,Master,
572 ! Exit, if the increment norm is smaller than pre-assigned tolerance.
573 if (finorm.lt.fimin) then
574 write (iout,*) 'Iteration converged'
581 ! Now, put together the histograms from all simulations, in order to get the
582 ! unbiased total histogram.
584 C Determine the minimum free energies
590 c write (9,'(3i5,f10.5)') i,(iparm,potE(i,iparm),iparm=1,nParmSet)
593 write (iout,'(2i5,21f8.2)') i,iparm,
594 & (enetb(k,i,iparm),k=1,21)
596 call restore_parm(iparm)
598 write (iout,*) wsc,wscp,welec,wvdwpp,wang,wtor,wscloc,
599 & wcorr,wcorr5,wcorr6,wturn4,wturn3,wturn6,wel_loc,
600 & wtor_d,wsccor,wbond
603 if (rescale_mode.eq.1) then
604 quot=1.0d0/(beta_h(ib,iparm)*1.987D-3*T0)
611 fT(l)=kfacl/(kfacl-1.0d0+quotl)
614 tt = 1.0d0/(beta_h(ib,iparm)*1.987D-3)
615 ft(6)=(320.0d0+80.0d0*dtanh((tt-320.0d0)/80.0d0))/320.0d0
617 ft(6)=1.0d0/(beta_h(ib,iparm)*1.987D-3*T0)
621 else if (rescale_mode.eq.2) then
622 quot=1.0d0/(T0*beta_h(ib,iparm)*1.987D-3)
626 fT(l)=1.12692801104297249644d0/
627 & dlog(dexp(quotl)+dexp(-quotl))
630 tt = 1.0d0/(beta_h(ib,iparm)*1.987D-3)
631 ft(6)=(320.0d0+80.0d0*dtanh((tt-320.0d0)/80.0d0))/320.0d0
633 ft(6)=1.0d0/(beta_h(ib,iparm)*1.987D-3*T0)
637 c write (iout,*) 1.0d0/(beta_h(ib,iparm)*1.987D-3),ft
638 else if (rescale_mode.eq.0) then
643 write (iout,*) "Error in WHAM_CALC: wrong RESCALE_MODE",
648 evdw=enetb(1,i,iparm)
649 evdw_t=enetb(21,i,iparm)
651 evdw2_14=enetb(17,i,iparm)
652 evdw2=enetb(2,i,iparm)+evdw2_14
654 evdw2=enetb(2,i,iparm)
659 evdw1=enetb(16,i,iparm)
664 ecorr=enetb(4,i,iparm)
665 ecorr5=enetb(5,i,iparm)
666 ecorr6=enetb(6,i,iparm)
667 eel_loc=enetb(7,i,iparm)
668 eello_turn3=enetb(8,i,iparm)
669 eello_turn4=enetb(9,i,iparm)
670 eturn6=enetb(10,i,iparm)
671 ebe=enetb(11,i,iparm)
672 escloc=enetb(12,i,iparm)
673 etors=enetb(13,i,iparm)
674 etors_d=enetb(14,i,iparm)
675 ehpb=enetb(15,i,iparm)
676 estr=enetb(18,i,iparm)
677 esccor=enetb(19,i,iparm)
678 edihcnstr=enetb(20,i,iparm)
680 write (iout,'(3i5,6f5.2,14f12.3)') i,ib,iparm,(ft(l),l=1,6),
681 & evdw+evdw_t,evdw2,ees,evdw1,ecorr,eel_loc,estr,ebe,escloc,
682 & etors,etors_d,eello_turn3,eello_turn4,esccor,edihcnstr
686 etot=wsc*(evdw+ft(6)*evdw_t)+wscp*evdw2+ft(1)*welec*ees
688 & +wang*ebe+ft(1)*wtor*etors+wscloc*escloc
689 & +wstrain*ehpb+nss*ebr+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5
690 & +ft(5)*wcorr6*ecorr6+ft(3)*wturn4*eello_turn4
691 & +ft(2)*wturn3*eello_turn3
692 & +ft(5)*wturn6*eturn6+ft(2)*wel_loc*eel_loc
693 & +edihcnstr+ft(2)*wtor_d*etors_d+ft(1)*wsccor*esccor
696 etot=wsc*(evdw+ft(6)*evdw_t)+wscp*evdw2
697 & +ft(1)*welec*(ees+evdw1)
698 & +wang*ebe+ft(1)*wtor*etors+wscloc*escloc
699 & +wstrain*ehpb+nss*ebr+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5
700 & +ft(5)*wcorr6*ecorr6+ft(3)*wturn4*eello_turn4
701 & +ft(2)*wturn3*eello_turn3
702 & +ft(5)*wturn6*eturn6+ft(2)*wel_loc*eel_loc+edihcnstr
703 & +ft(2)*wtor_d*etors_d+ft(1)*wsccor*esccor
706 c write (iout,*) "i",i," ib",ib,
707 c & " temp",1.0d0/(1.987d-3*beta_h(ib,iparm))," etot",etot,
708 c & " entfac",entfac(i)
709 etot=etot-entfac(i)/beta_h(ib,iparm)
710 if(etot.lt.potEmin_all(ib,iparm)) potEmin_all(ib,iparm)=etot
711 c write (iout,*) "efree",etot," potEmin",potEmin_all(ib,iparm)
716 write (iout,*) "The potEmin array before reduction"
718 write (iout,*) "Parameter set",i
720 write (iout,*) j,PotEmin_all(j,i)
723 write (iout,*) "potEmin_min",potEmin_min
726 C Determine the minimum energes for all parameter sets and temperatures
727 call MPI_AllReduce(potEmin_all(1,1),potEmin_t_all(1,1),
728 & maxT_h*nParmSet,MPI_DOUBLE_PRECISION,MPI_MIN,WHAM_COMM,IERROR)
731 potEmin_all(j,i)=potEmin_t_all(j,i)
735 potEmin_min=potEmin_all(1,1)
738 if (potEmin_all(j,i).lt.potEmin_min)
739 & potEmin_min=potEmin_all(j,i)
743 write (iout,*) "The potEmin array"
745 write (iout,*) "Parameter set",i
747 write (iout,*) j,PotEmin_all(j,i)
750 write (iout,*) "potEmin_min",potEmin_min
762 write (iout,*) "--------------hist"
766 sumW_p(i,iparm)=0.0d0
767 sumE_p(i,iparm)=0.0d0
768 sumEbis_p(i,iparm)=0.0d0
769 sumEsq_p(i,iparm)=0.0d0
771 sumQ_p(j,i,iparm)=0.0d0
772 sumQsq_p(j,i,iparm)=0.0d0
773 sumEQ_p(j,i,iparm)=0.0d0
783 sumEbis(i,iparm)=0.0d0
784 sumEsq(i,iparm)=0.0d0
786 sumQ(j,i,iparm)=0.0d0
787 sumQsq(j,i,iparm)=0.0d0
788 sumEQ(j,i,iparm)=0.0d0
794 c 8/26/05 entropy distribution
799 c ent=-dlog(entfac(t))
801 if (ent.lt.entmin_p) entmin_p=ent
802 if (ent.gt.entmax_p) entmax_p=ent
804 write (iout,*) "entmin",entmin_p," entmax",entmax_p
806 call MPI_Allreduce(entmin_p,entmin,1,MPI_DOUBLE_PRECISION,MPI_MIN,
808 call MPI_Allreduce(entmax_p,entmax,1,MPI_DOUBLE_PRECISION,MPI_MAX,
810 ientmax=entmax-entmin
811 if (ientmax.gt.2000) ientmax=2000
812 write (iout,*) "entmin",entmin," entmax",entmax," ientmax",ientmax
815 c ient=-dlog(entfac(t))-entmin
816 ient=entfac(t)-entmin
817 if (ient.le.2000) histent_p(ient)=histent_p(ient)+1
819 call MPI_Allreduce(histent_p(0),histent(0),ientmax+1,MPI_INTEGER,
820 & MPI_SUM,WHAM_COMM,IERROR)
821 if (me1.eq.Master) then
822 write (iout,*) "Entropy histogram"
824 write(iout,'(f15.4,i10)') entmin+i,histent(i)
832 if (ent.lt.entmin) entmin=ent
833 if (ent.gt.entmax) entmax=ent
835 ientmax=-dlog(entmax)-entmin
836 if (ientmax.gt.2000) ientmax=2000
838 ient=entfac(t)-entmin
839 if (ient.le.2000) histent(ient)=histent(ient)+1
841 write (iout,*) "Entropy histogram"
843 write(iout,'(2f15.4)') entmin+i,histent(i)
848 c write (iout,*) "me1",me1," scount",scount(me1)
874 hrmsrgy(j,i,ib)=0.0d0
876 hrmsrgy_p(j,i,ib)=0.0d0
888 hfin_ent_p(ind)=hfin_ent_p(ind)+dexp(entfac(t))
890 hfin_ent(ind)=hfin_ent(ind)+dexp(entfac(t))
892 call restore_parm(iparm)
893 evdw=enetb(21,t,iparm)
894 evdw_t=enetb(1,t,iparm)
896 evdw2_14=enetb(17,t,iparm)
897 evdw2=enetb(2,t,iparm)+evdw2_14
899 evdw2=enetb(2,t,iparm)
904 evdw1=enetb(16,t,iparm)
909 ecorr=enetb(4,t,iparm)
910 ecorr5=enetb(5,t,iparm)
911 ecorr6=enetb(6,t,iparm)
912 eel_loc=enetb(7,t,iparm)
913 eello_turn3=enetb(8,t,iparm)
914 eello_turn4=enetb(9,t,iparm)
915 eturn6=enetb(10,t,iparm)
916 ebe=enetb(11,t,iparm)
917 escloc=enetb(12,t,iparm)
918 etors=enetb(13,t,iparm)
919 etors_d=enetb(14,t,iparm)
920 ehpb=enetb(15,t,iparm)
921 estr=enetb(18,t,iparm)
922 esccor=enetb(19,t,iparm)
923 edihcnstr=enetb(20,t,iparm)
925 betaT=startGridT+k*delta_T
929 if (rescale_mode.eq.1) then
937 denom=kfacl-1.0d0+quotl
939 ftprim(l)=-l*ft(l)*quotl1/(T0*denom)
940 ftbis(l)=l*kfacl*quotl1*
941 & (2*l*quotl-(l-1)*denom)/(quot*t0*t0*denom**3)
944 ft(6)=(320.0d0+80.0d0*dtanh((betaT-320.0d0)/80.0d0))/
946 ftprim(6)=1.0d0/(320.0d0*dcosh((betaT-320.0d0)/80.0d0)**2)
947 ftbis(6)=-2.0d0*dtanh((betaT-320.0d0)/80.0d0)
948 & /(320.0d0*80.0d0*dcosh((betaT-320.0d0)/80.0d0)**3)
958 else if (rescale_mode.eq.2) then
966 logfac=1.0d0/dlog(eplus+eminus)
967 tanhT=(eplus-eminus)/(eplus+eminus)
968 fT(l)=1.12692801104297249644d0*logfac
969 ftprim(l)=-l*quotl1*ft(l)*tanhT*logfac/T0
970 ftbis(l)=(l-1)*ftprim(l)/(quot*T0)-
971 & 2*l*quotl1/T0*logfac*
972 & (2*l*quotl1*ft(l)/(T0*(eplus+eminus)**2)
976 ft(6)=(320.0d0+80.0d0*dtanh((betaT-320.0d0)/80.0d0))/
978 ftprim(6)=1.0d0/(320.0d0*dcosh((betaT-320.0d0)/80.0d0)**2)
979 ftbis(6)=-2.0d0*dtanh((betaT-320.0d0)/80.0d0)
980 & /(320.0d0*80.0d0*dcosh((betaT-320.0d0)/80.0d0)**3)
990 else if (rescale_mode.eq.0) then
996 write (iout,*) "Error in WHAM_CALC: wrong RESCALE_MODE",
1001 c write (iout,*) "ftprim",ftprim
1002 c write (iout,*) "ftbis",ftbis
1003 betaT=1.0d0/(1.987D-3*betaT)
1004 if (betaT.ge.beta_h(1,iparm)) then
1005 potEmin=potEmin_all(1,iparm)
1006 c write(iout,*) "first",temper,potEmin
1007 else if (betaT.lt.beta_h(nT_h(iparm),iparm)) then
1008 potEmin=potEmin_all(nT_h(iparm),iparm)
1009 c write (iout,*) "last",temper,potEmin
1011 do l=1,nT_h(iparm)-1
1012 if (betaT.le.beta_h(l,iparm) .and.
1013 & betaT.gt.beta_h(l+1,iparm)) then
1014 potEmin=potEmin_all(l,iparm)
1015 c write (iout,*) "l",l,
1016 c & betaT,1.0d0/(1.987D-3*beta_h(l,iparm)),
1017 c & 1.0d0/(1.987D-3*beta_h(l+1,iparm)),temper,potEmin
1022 c write (iout,*) ib," PotEmin",potEmin
1024 etot=wsc*(evdw+ft(6)*evdw_t)+wscp*evdw2+ft(1)*welec*ees
1026 & +wang*ebe+ft(1)*wtor*etors+wscloc*escloc
1027 & +wstrain*ehpb+nss*ebr+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5
1028 & +ft(5)*wcorr6*ecorr6+ft(3)*wturn4*eello_turn4
1029 & +ft(2)*wturn3*eello_turn3
1030 & +ft(5)*wturn6*eturn6+ft(2)*wel_loc*eel_loc
1031 & +edihcnstr+ft(2)*wtor_d*etors_d+ft(1)*wsccor*esccor
1033 eprim=ftprim(6)*evdw_t+ftprim(1)*welec*ees
1034 & +ftprim(1)*wtor*etors+
1035 & ftprim(3)*wcorr*ecorr+ftprim(4)*wcorr5*ecorr5+
1036 & ftprim(5)*wcorr6*ecorr6+ftprim(3)*wturn4*eello_turn4+
1037 & ftprim(2)*wturn3*eello_turn3+ftprim(5)*wturn6*eturn6+
1038 & ftprim(2)*wel_loc*eel_loc+ftprim(2)*wtor_d*etors_d+
1039 & ftprim(1)*wsccor*esccor
1040 ebis=ftbis(1)*welec*ees+ftbis(1)*wtor*etors+
1041 & ftbis(3)*wcorr*ecorr+ftbis(4)*wcorr5*ecorr5+
1042 & ftbis(5)*wcorr6*ecorr6+ftbis(3)*wturn4*eello_turn4+
1043 & ftbis(2)*wturn3*eello_turn3+ftbis(5)*wturn6*eturn6+
1044 & ftbis(2)*wel_loc*eel_loc+ftbis(2)*wtor_d*etors_d+
1045 & ftbis(1)*wsccor*esccor
1047 etot=wsc*(evdw+ft(6)*evdw_t)+wscp*evdw2
1048 & +ft(1)*welec*(ees+evdw1)
1049 & +wang*ebe+ft(1)*wtor*etors+wscloc*escloc
1050 & +wstrain*ehpb+nss*ebr+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5
1051 & +ft(5)*wcorr6*ecorr6+ft(3)*wturn4*eello_turn4
1052 & +ft(2)*wturn3*eello_turn3
1053 & +ft(5)*wturn6*eturn6+ft(2)*wel_loc*eel_loc+edihcnstr
1054 & +ft(2)*wtor_d*etors_d+ft(1)*wsccor*esccor
1056 eprim=ftprim(6)*evdw_t+ftprim(1)*welec*(ees+evdw1)
1057 & +ftprim(1)*wtor*etors+
1058 & ftprim(3)*wcorr*ecorr+ftprim(4)*wcorr5*ecorr5+
1059 & ftprim(5)*wcorr6*ecorr6+ftprim(3)*wturn4*eello_turn4+
1060 & ftprim(2)*wturn3*eello_turn3+ftprim(5)*wturn6*eturn6+
1061 & ftprim(2)*wel_loc*eel_loc+ftprim(2)*wtor_d*etors_d+
1062 & ftprim(1)*wsccor*esccor
1063 ebis=ftbis(1)*welec*(ees+evdw1)+ftbis(1)*wtor*etors+
1064 & ftbis(3)*wcorr*ecorr+ftbis(4)*wcorr5*ecorr5+
1065 & ftbis(5)*wcorr6*ecorr6+ftbis(3)*wturn4*eello_turn4+
1066 & ftbis(2)*wturn3*eello_turn3+ftbis(5)*wturn6*eturn6+
1067 & ftbis(2)*wel_loc*eel_loc+ftbis(2)*wtor_d*etors_d+
1068 & ftprim(1)*wsccor*esccor
1070 weight=dexp(-betaT*(etot-potEmin)+entfac(t))
1072 write (iout,*) "iparm",iparm," t",t," temper",temper,
1073 & " etot",etot," entfac",entfac(t),
1074 & " efree",etot-entfac(t)/betaT," potEmin",potEmin,
1075 & " boltz",-betaT*(etot-potEmin)+entfac(t),
1076 & " weight",weight," ebis",ebis
1078 etot=etot-temper*eprim
1080 sumW_p(k,iparm)=sumW_p(k,iparm)+weight
1081 sumE_p(k,iparm)=sumE_p(k,iparm)+etot*weight
1082 sumEbis_p(k,iparm)=sumEbis_p(k,iparm)+ebis*weight
1083 sumEsq_p(k,iparm)=sumEsq_p(k,iparm)+etot**2*weight
1085 sumQ_p(j,k,iparm)=sumQ_p(j,k,iparm)+q(j,t)*weight
1086 sumQsq_p(j,k,iparm)=sumQsq_p(j,k,iparm)+q(j,t)**2*weight
1087 sumEQ_p(j,k,iparm)=sumEQ_p(j,k,iparm)
1088 & +etot*q(j,t)*weight
1091 sumW(k,iparm)=sumW(k,iparm)+weight
1092 sumE(k,iparm)=sumE(k,iparm)+etot*weight
1093 sumEbis(k,iparm)=sumEbis(k,iparm)+ebis*weight
1094 sumEsq(k,iparm)=sumEsq(k,iparm)+etot**2*weight
1096 sumQ(j,k,iparm)=sumQ(j,k,iparm)+q(j,t)*weight
1097 sumQsq(j,k,iparm)=sumQsq(j,k,iparm)+q(j,t)**2*weight
1098 sumEQ(j,k,iparm)=sumEQ(j,k,iparm)
1099 & +etot*q(j,t)*weight
1103 indE = aint(potE(t,iparm)-aint(potEmin))
1104 if (indE.ge.0 .and. indE.le.maxinde) then
1105 if (indE.gt.upindE_p) upindE_p=indE
1106 histE_p(indE)=histE_p(indE)+dexp(-entfac(t))
1110 potEmin=potEmin_all(ib,iparm)
1111 expfac=dexp(-beta_h(ib,iparm)*(etot-potEmin)+entfac(t))
1112 hfin_p(ind,ib)=hfin_p(ind,ib)+
1113 & dexp(-beta_h(ib,iparm)*(etot-potEmin)+entfac(t))
1115 indrgy=dint((q(nQ+2,t)-rgymin)/deltrgy)
1116 indrms=dint((q(nQ+1,t)-rmsmin)/deltrms)
1117 hrmsrgy_p(indrgy,indrms,ib)=
1118 & hrmsrgy_p(indrgy,indrms,ib)+expfac
1123 potEmin=potEmin_all(ib,iparm)
1124 expfac=dexp(-beta_h(ib,iparm)*(etot-potEmin)+entfac(t))
1125 hfin(ind,ib)=hfin(ind,ib)+
1126 & dexp(-beta_h(ib,iparm)*(etot-potEmin)+entfac(t))
1128 indrgy=dint((q(nQ+2,t)-rgymin)/deltrgy)
1129 indrms=dint((q(nQ+1,t)-rmsmin)/deltrms)
1130 hrmsrgy(indrgy,indrms,ib)=
1131 & hrmsrgy(indrgy,indrms,ib)+expfac
1137 if (histout) call MPI_Reduce(hfin_p(0,ib),hfin(0,ib),nbin,
1138 & MPI_DOUBLE_PRECISION,MPI_SUM,Master,WHAM_COMM,IERROR)
1140 call MPI_Reduce(hrmsrgy_p(0,0,ib),hrmsrgy(0,0,ib),
1141 & (MaxBinRgy+1)*(nbin_rms+1),MPI_DOUBLE_PRECISION,MPI_SUM,Master,
1145 call MPI_Reduce(upindE_p,upindE,1,
1146 & MPI_INTEGER,MPI_MAX,Master,WHAM_COMM,IERROR)
1147 call MPI_Reduce(histE_p(0),histE(0),maxindE,
1148 & MPI_DOUBLE_PRECISION,MPI_SUM,Master,WHAM_COMM,IERROR)
1150 if (me1.eq.master) then
1154 write (iout,'(6x,$)')
1155 write (iout,'(f20.2,$)') (1.0d0/(1.987D-3*beta_h(ib,iparm)),
1159 write (iout,'(/a)') 'Final histograms'
1161 if (nslice.eq.1) then
1162 if (separate_parset) then
1163 write(licz3,"(bz,i3.3)") myparm
1164 histname=prefix(:ilen(prefix))//'_par'//licz3//'.hist'
1166 histname=prefix(:ilen(prefix))//'.hist'
1169 if (separate_parset) then
1170 write(licz3,"(bz,i3.3)") myparm
1171 histname=prefix(:ilen(prefix))//'_par'//licz3//
1172 & '_slice_'//licz2//'.hist'
1174 histname=prefix(:ilen(prefix))//'_slice_'//licz2//'.hist'
1177 #if defined(AIX) || defined(PGI)
1178 open (ihist,file=histname,position='append')
1180 open (ihist,file=histname,access='append')
1188 sumH=sumH+hfin(t,ib)
1190 if (sumH.gt.0.0d0) then
1192 jj = mod(liczba,nbin1)
1194 write (iout,'(f6.3,$)') dmin+(jj+0.5d0)*delta
1196 & write (ihist,'(f6.3,$)') dmin+(jj+0.5d0)*delta
1199 write (iout,'(e20.10,$)') hfin(t,ib)
1200 if (histfile) write (ihist,'(e20.10,$)') hfin(t,ib)
1202 write (iout,'(i5)') iparm
1203 if (histfile) write (ihist,'(i5)') iparm
1210 if (nslice.eq.1) then
1211 if (separate_parset) then
1212 write(licz3,"(bz,i3.3)") myparm
1213 histname=prefix(:ilen(prefix))//"_par"//licz3//'.ent'
1215 histname=prefix(:ilen(prefix))//'.ent'
1218 if (separate_parset) then
1219 write(licz3,"(bz,i3.3)") myparm
1220 histname=prefix(:ilen(prefix))//'par_'//licz3//
1221 & '_slice_'//licz2//'.ent'
1223 histname=prefix(:ilen(prefix))//'_slice_'//licz2//'.ent'
1226 #if defined(AIX) || defined(PGI)
1227 open (ihist,file=histname,position='append')
1229 open (ihist,file=histname,access='append')
1231 write (ihist,'(a)') "# Microcanonical entropy"
1233 write (ihist,'(f8.0,$)') dint(potEmin)+i
1234 if (histE(i).gt.0.0e0) then
1235 write (ihist,'(f15.5,$)') dlog(histE(i))
1237 write (ihist,'(f15.5,$)') 0.0d0
1243 write (iout,*) "Microcanonical entropy"
1245 write (iout,'(f8.0,$)') dint(potEmin)+i
1246 if (histE(i).gt.0.0e0) then
1247 write (iout,'(f15.5,$)') dlog(histE(i))
1249 write (iout,'(f15.5,$)') 0.0d0
1254 if (nslice.eq.1) then
1255 if (separate_parset) then
1256 write(licz3,"(bz,i3.3)") myparm
1257 histname=prefix(:ilen(prefix))//'_par'//licz3//'.rmsrgy'
1259 histname=prefix(:ilen(prefix))//'.rmsrgy'
1262 if (separate_parset) then
1263 write(licz3,"(bz,i3.3)") myparm
1264 histname=prefix(:ilen(prefix))//'_par'//licz3//
1265 & '_slice_'//licz2//'.rmsrgy'
1267 histname=prefix(:ilen(prefix))//'_slice_'//licz2//'.rmsrgy'
1270 #if defined(AIX) || defined(PGI)
1271 open (ihist,file=histname,position='append')
1273 open (ihist,file=histname,access='append')
1277 write(ihist,'(2f8.2,$)')
1278 & rgymin+deltrgy*j,rmsmin+deltrms*i
1280 if (hrmsrgy(j,i,ib).gt.0.0d0) then
1281 write(ihist,'(e14.5,$)')
1282 & -dlog(hrmsrgy(j,i,ib))/beta_h(ib,iparm)
1285 write(ihist,'(e14.5,$)') 1.0d6
1288 write (ihist,'(i2)') iparm
1296 call MPI_Reduce(hfin_ent_p(0),hfin_ent(0),nbin,
1297 & MPI_DOUBLE_PRECISION,MPI_SUM,Master,WHAM_COMM,IERROR)
1298 call MPI_Reduce(sumW_p(0,1),sumW(0,1),(nGridT+1)*nParmSet,
1299 & MPI_DOUBLE_PRECISION,MPI_SUM,Master,WHAM_COMM,IERROR)
1300 call MPI_Reduce(sumE_p(0,1),sumE(0,1),(nGridT+1)*nParmSet,
1301 & MPI_DOUBLE_PRECISION,MPI_SUM,Master,WHAM_COMM,IERROR)
1302 call MPI_Reduce(sumEbis_p(0,1),sumEbis(0,1),(nGridT+1)*nParmSet,
1303 & MPI_DOUBLE_PRECISION,MPI_SUM,Master,WHAM_COMM,IERROR)
1304 call MPI_Reduce(sumEsq_p(0,1),sumEsq(0,1),(nGridT+1)*nParmSet,
1305 & MPI_DOUBLE_PRECISION,MPI_SUM,Master,WHAM_COMM,IERROR)
1306 call MPI_Reduce(sumQ_p(1,0,1),sumQ(1,0,1),
1307 & MaxQ1*(nGridT+1)*nParmSet,MPI_DOUBLE_PRECISION,MPI_SUM,Master,
1309 call MPI_Reduce(sumQsq_p(1,0,1),sumQsq(1,0,1),
1310 & MaxQ1*(nGridT+1)*nParmSet,MPI_DOUBLE_PRECISION,MPI_SUM,Master,
1312 call MPI_Reduce(sumEQ_p(1,0,1),sumEQ(1,0,1),
1313 & MaxQ1*(nGridT+1)*nParmSet,MPI_DOUBLE_PRECISION,MPI_SUM,Master,
1315 if (me.eq.master) then
1317 write (iout,'(/a)') 'Thermal characteristics of folding'
1318 if (nslice.eq.1) then
1321 nazwa=prefix(:ilen(prefix))//"_slice_"//licz2
1324 if (nparmset.eq.1 .and. .not.separate_parset) then
1325 nazwa=nazwa(:iln)//".thermal"
1326 else if (nparmset.eq.1 .and. separate_parset) then
1327 write(licz3,"(bz,i3.3)") myparm
1328 nazwa=nazwa(:iln)//"_par_"//licz3//".thermal"
1331 if (nparmset.gt.1) then
1332 write(licz3,"(bz,i3.3)") iparm
1333 nazwa=nazwa(:iln)//"_par_"//licz3//".thermal"
1336 if (separate_parset) then
1337 write (iout,'(a,i3)') "Parameter set",myparm
1339 write (iout,'(a,i3)') "Parameter set",iparm
1342 betaT=1.0d0/(1.987D-3*(startGridT+i*delta_T))
1343 if (betaT.ge.beta_h(1,iparm)) then
1344 potEmin=potEmin_all(1,iparm)
1345 else if (betaT.lt.beta_h(nT_h(iparm),iparm)) then
1346 potEmin=potEmin_all(nT_h(iparm),iparm)
1348 do l=1,nT_h(iparm)-1
1349 if (betaT.le.beta_h(l,iparm) .and.
1350 & betaT.gt.beta_h(l+1,iparm)) then
1351 potEmin=potEmin_all(l,iparm)
1356 sumE(i,iparm)=sumE(i,iparm)/sumW(i,iparm)
1357 sumEbis(i,iparm)=(startGridT+i*delta_T)*sumEbis(i,iparm)/
1359 sumEsq(i,iparm)=(sumEsq(i,iparm)/sumW(i,iparm)
1360 & -sumE(i,iparm)**2)/(1.987D-3*(startGridT+i*delta_T)**2)
1362 sumQ(j,i,iparm)=sumQ(j,i,iparm)/sumW(i,iparm)
1363 sumQsq(j,i,iparm)=sumQsq(j,i,iparm)/sumW(i,iparm)
1364 & -sumQ(j,i,iparm)**2
1365 sumEQ(j,i,iparm)=sumEQ(j,i,iparm)/sumW(i,iparm)
1366 & -sumQ(j,i,iparm)*sumE(i,iparm)
1368 sumW(i,iparm)=-dlog(sumW(i,iparm))*(1.987D-3*
1369 & (startGridT+i*delta_T))+potEmin
1370 write (iout,'(f7.1,2f15.5,$)') startGridT+i*delta_T,
1371 & sumW(i,iparm),sumE(i,iparm)
1372 write (iout,'(f10.5,$)') (sumQ(j,i,iparm),j=1,nQ+2)
1373 write (iout,'(e15.5,$)') sumEsq(i,iparm)-sumEbis(i,iparm),
1374 & (sumQsq(j,i,iparm),j=1,nQ+2),(sumEQ(j,i,iparm),j=1,nQ+2)
1376 write (34,'(f7.1,2f15.5,$)') startGridT+i*delta_T,
1377 & sumW(i,iparm),sumE(i,iparm)
1378 write (34,'(f10.5,$)') (sumQ(j,i,iparm),j=1,nQ+2)
1379 write (34,'(e15.5,$)') sumEsq(i,iparm)-sumEbis(i,iparm),
1380 & (sumQsq(j,i,iparm),j=1,nQ+2),(sumEQ(j,i,iparm),j=1,nQ+2)
1387 if (hfin_ent(t).gt.0.0d0) then
1389 jj = mod(liczba,nbin1)
1390 write (iout,'(f6.3,e20.10," ent")') dmin+(jj+0.5d0)*delta,
1392 if (histfile) write (ihist,'(f6.3,e20.10," ent")')
1393 & dmin+(jj+0.5d0)*delta,
1397 if (histfile) close(ihist)
1401 ! Write data for zscore
1402 if (nslice.eq.1) then
1403 zscname=prefix(:ilen(prefix))//".zsc"
1405 zscname=prefix(:ilen(prefix))//"_slice_"//licz2//".zsc"
1407 #if defined(AIX) || defined(PGI)
1408 open (izsc,file=prefix(:ilen(prefix))//'.zsc',position='append')
1410 open (izsc,file=prefix(:ilen(prefix))//'.zsc',access='append')
1412 write (izsc,'("NQ=",i1," NPARM=",i1)') nQ,nParmSet
1414 write (izsc,'("NT=",i1)') nT_h(iparm)
1416 write (izsc,'("TEMP=",f6.1," NR=",i2," SNK=",$)')
1417 & 1.0d0/(beta_h(ib,iparm)*1.987D-3),nR(ib,iparm)
1418 jj = min0(nR(ib,iparm),7)
1419 write (izsc,'(i8,$)') (snk(i,ib,iparm,islice),i=1,jj)
1420 write (izsc,'(a1,$)') (" ",i=22+8*jj+1,79)
1421 write (izsc,'("&")')
1422 if (nR(ib,iparm).gt.7) then
1423 do ii=8,nR(ib,iparm),9
1424 jj = min0(nR(ib,iparm),ii+8)
1425 write (izsc,'(i8,$)') (snk(i,ib,iparm,islice),i=ii,jj)
1426 write (izsc,'(a1,$') (" ",i=(jj-ii+1)*8+1,79)
1427 write (izsc,'("&")')
1430 write (izsc,'("FI=",$)')
1431 jj=min0(nR(ib,iparm),7)
1432 write (izsc,'(f10.5,$)') (fi(i,ib,iparm),i=1,jj)
1433 write (izsc,'(a1,$)') (" ",i=3+10*jj+1,79)
1434 write (izsc,'("&")')
1435 if (nR(ib,iparm).gt.7) then
1436 do ii=8,nR(ib,iparm),9
1437 jj = min0(nR(ib,iparm),ii+8)
1438 write (izsc,'(f10.5,$)') (fi(i,ib,iparm),i=ii,jj)
1439 if (jj.eq.nR(ib,iparm)) then
1442 write (izsc,'(a1,$)') (" ",i=10*(jj-ii+1)+1,79)
1443 write (izsc,'(t80,"&")')
1448 write (izsc,'("KH=",$)')
1449 write (izsc,'(f7.2,$)') (Kh(j,i,ib,iparm),j=1,nQ)
1450 write (izsc,'(" Q0=",$)')
1451 write (izsc,'(f7.5,$)') (q0(j,i,ib,iparm),j=1,nQ)
1466 >>>>>>> e183793... Added src_MD-M-newcorr (Adasko's source) and src-NEWSC of WHAM (with Momo's SCSC potentials)