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)
21 parameter (MaxHdim=200)
23 parameter (maxinde=200)
27 integer ierror,errcode,status(MPI_STATUS_SIZE)
29 include "COMMON.CONTROL"
30 include "COMMON.IOUNITS"
32 include "COMMON.ENERGIES"
33 include "COMMON.FFIELD"
34 include "COMMON.SBRIDGE"
36 include "COMMON.ENEPS"
37 integer MaxPoint,MaxPointProc
38 parameter (MaxPoint=MaxStr,
39 & MaxPointProc=MaxStr_Proc)
40 double precision finorm_max,potfac,entmin,entmax,expfac,vf
41 double precision entfac_min,entfac_min_t
42 parameter (finorm_max=1.0d0)
44 integer i,ii,j,jj,k,kk,l,m,ind,iter,t,tmax,ient,ientmax,iln
45 integer start,end,iharm,ib,iib,nbin1,nbin,nbin_rms,nbin_rgy,
46 & nbin_rmsrgy,liczba,iparm,nFi,indrgy,indrms
47 integer htot(0:MaxHdim),histent(0:2000)
48 double precision v(MaxPointProc,MaxR,MaxT_h,Max_Parm)
49 double precision energia(0:max_ene)
51 integer tmax_t,upindE_p
52 double precision fi_p(MaxR,MaxT_h,Max_Parm),
53 & fi_p_min(MaxR,MaxT_h,Max_Parm)
54 double precision sumW_p(0:Max_GridT,Max_Parm),
55 & sumE_p(0:Max_GridT,Max_Parm),sumEsq_p(0:Max_GridT,Max_Parm),
56 & sumQ_p(MaxQ1,0:Max_GridT,Max_Parm),
57 & sumQsq_p(MaxQ1,0:Max_GridT,Max_Parm),
58 & sumEQ_p(MaxQ1,0:Max_GridT,Max_Parm),
59 & sumEprim_p(MaxQ1,0:Max_GridT,Max_Parm),
60 & sumEbis_p(0:Max_GridT,Max_Parm)
61 double precision hfin_p(0:MaxHdim,maxT_h),
62 & hfin_ent_p(0:MaxHdim),histE_p(0:maxindE),sumH,
63 & hrmsrgy_p(0:MaxBinRgy,0:MaxBinRms,maxT_h)
64 double precision rgymin_t,rmsmin_t,rgymax_t,rmsmax_t
65 double precision potEmin_t_all(maxT_h,Max_Parm),entmin_p,entmax_p
66 integer histent_p(0:2000)
67 logical lprint /.true./
69 double precision rgymin,rmsmin,rgymax,rmsmax
70 double precision sumW(0:NGridT,Max_Parm),sumE(0:NGridT,Max_Parm),
71 & sumEsq(0:NGridT,Max_Parm),sumQ(MaxQ1,0:NGridT,Max_Parm),
72 & sumQsq(MaxQ1,0:NGridT,Max_Parm),sumEQ(MaxQ1,0:NGridT,Max_Parm),
73 & sumEprim(0:NGridT,Max_Parm),sumEbis(0:NGridT,Max_Parm),betaT,
75 double precision fi(MaxR,maxT_h,Max_Parm),
76 & fi_min(MaxR,maxT_h,Max_Parm),
77 & dd,dd1,dd2,hh,dmin,denom,finorm,avefi,pom,
78 & hfin(0:MaxHdim,maxT_h),histE(0:maxindE),
79 & hrmsrgy(0:MaxBinRgy,0:MaxBinRms,maxT_h),
80 & potEmin_all(maxT_h,Max_Parm),potEmin,potEmin_min,ent,
81 & hfin_ent(0:MaxHdim),vmax,aux
82 double precision fT(6),fTprim(6),fTbis(6),quot,quotl1,quotl,kfacl,
83 & eprim,ebis,temper,kfac/2.4d0/,T0/300.0d0/,
84 & eplus,eminus,logfac,tanhT,tt
85 double precision etot,evdw,evdw_t,evdw2,ees,evdw1,ebe,etors,
86 & escloc,ehpb,ecorr,ecorr5,ecorr6,eello_turn4,eello_turn3,
87 & eturn6,eel_loc,edihcnstr,etors_d,estr,evdw2_14,esccor,
88 & ehomology_constr,edfadis,edfator,edfanei,edfabet
91 integer ind_point(maxpoint),upindE,indE
99 write(licz2,'(bz,i2.2)') islice
101 write (iout,'(//80(1h-)/"Solving WHAM equations for slice",
102 & i2/80(1h-)//)') islice
103 write (iout,*) "delta",delta," nbin1",nbin1
104 write (iout,*) "MaxN",MaxN," MaxQ",MaxQ," MaHdim",MaxHdim
110 potEmin_all(j,i)=1.0d10
125 if (q(nQ+1,i).lt.rmsmin) rmsmin=q(nQ+1,i)
126 if (q(nQ+1,i).gt.rmsmax) rmsmax=q(nQ+1,i)
127 if (q(nQ+2,i).lt.rgymin) rgymin=q(nQ+2,i)
128 if (q(nQ+2,i).gt.rgymax) rgymax=q(nQ+2,i)
131 ind=(q(j,i)-dmin+1.0d-8)/delta
133 ind_point(i)=ind_point(i)+ind
135 ind_point(i)=ind_point(i)+nbin1**(j-1)*ind
137 if (ind_point(i).lt.0 .or. ind_point(i).gt.MaxHdim) then
138 write (iout,*) "Error - index exceeds range for point",i,
139 & " q=",q(j,i)," ind",ind_point(i)
141 write (iout,*) "Processor",me1
143 call MPI_Abort(MPI_COMM_WORLD, Ierror, Errcode )
148 if (ind_point(i).gt.tmax) tmax=ind_point(i)
149 htot(ind_point(i))=htot(ind_point(i))+1
151 write (iout,*) "i",i,"q",(q(j,i),j=1,nQ)," ind",ind_point(i),
152 & " htot",htot(ind_point(i))
159 write (iout,'(a)') "Numbers of counts in Q bins"
161 if (htot(t).gt.0) then
162 write (iout,'(i15,$)') t
165 jj = mod(liczba,nbin1)
167 write (iout,'(i5,$)') jj
169 write (iout,'(i8)') htot(t)
173 write (iout,'(a,i3)') "Number of data points for parameter set",
175 write (iout,'(i7,$)') ((snk(m,ib,iparm,islice),m=1,nr(ib,iparm)),
177 write (iout,'(i8)') stot(islice)
183 call MPI_AllReduce(tmax,tmax_t,1,MPI_INTEGER,MPI_MAX,
186 call MPI_AllReduce(rmsmin,rmsmin_t,1,MPI_DOUBLE_PRECISION,
187 & MPI_MIN,WHAM_COMM,IERROR)
188 call MPI_AllReduce(rmsmax,rmsmax_t,1,MPI_DOUBLE_PRECISION,
189 & MPI_MAX,WHAM_COMM,IERROR)
190 call MPI_AllReduce(rgymin,rgymin_t,1,MPI_DOUBLE_PRECISION,
191 & MPI_MIN,WHAM_COMM,IERROR)
192 call MPI_AllReduce(rgymax,rgymax_t,1,MPI_DOUBLE_PRECISION,
193 & MPI_MAX,WHAM_COMM,IERROR)
199 rmsmin=deltrms*dint(rmsmin/deltrms)
200 rmsmax=deltrms*dint(rmsmax/deltrms)
201 rgymin=deltrms*dint(rgymin/deltrgy)
202 rgymax=deltrms*dint(rgymax/deltrgy)
203 nbin_rms=(rmsmax-rmsmin)/deltrms
204 nbin_rgy=(rgymax-rgymin)/deltrgy
205 write (iout,*) "rmsmin",rmsmin," rmsmax",rmsmax," rgymin",rgymin,
206 & " rgymax",rgymax," nbin_rms",nbin_rms," nbin_rgy",nbin_rgy
213 write (iout,*) "nFi",nFi
214 ! Compute the Boltzmann factor corresponing to restrain potentials in different
221 c write (9,'(3i5,f10.5)') i,(iparm,potE(i,iparm),iparm=1,nParmSet)
224 write (iout,'(2i5,22f8.2)') i,iparm,
225 & (enetb(k,i,iparm),k=1,22)
227 call restore_parm(iparm)
229 write (iout,*) wsc,wscp,welec,wvdwpp,wang,wtor,wscloc,
230 & wcorr,wcorr5,wcorr6,wturn4,wturn3,wturn6,wel_loc,
231 & wtor_d,wsccor,wbond
234 if (rescale_mode.eq.1) then
235 quot=1.0d0/(beta_h(ib,iparm)*1.987D-3*T0)
242 fT(l)=kfacl/(kfacl-1.0d0+quotl)
245 tt = 1.0d0/(beta_h(ib,iparm)*1.987D-3)
246 ft(6)=(320.0d0+80.0d0*dtanh((tt-320.0d0)/80.0d0))/320.0d0
248 ft(6)=1.0d0/(beta_h(ib,iparm)*1.987D-3*T0)
252 else if (rescale_mode.eq.2) then
253 quot=1.0d0/(T0*beta_h(ib,iparm)*1.987D-3)
257 fT(l)=1.12692801104297249644d0/
258 & dlog(dexp(quotl)+dexp(-quotl))
261 tt = 1.0d0/(beta_h(ib,iparm)*1.987D-3)
262 ft(6)=(320.0d0+80.0d0*dtanh((tt-320.0d0)/80.0d0))/320.0d0
264 ft(6)=1.0d0/(beta_h(ib,iparm)*1.987D-3*T0)
268 c write (iout,*) 1.0d0/(beta_h(ib,iparm)*1.987D-3),ft
269 else if (rescale_mode.eq.0) then
274 write (iout,*) "Error in WHAM_CALC: wrong RESCALE_MODE",
279 evdw=enetb(1,i,iparm)
280 evdw_t=enetb(21,i,iparm)
282 evdw2_14=enetb(17,i,iparm)
283 evdw2=enetb(2,i,iparm)+evdw2_14
285 evdw2=enetb(2,i,iparm)
290 evdw1=enetb(16,i,iparm)
295 ecorr=enetb(4,i,iparm)
296 ecorr5=enetb(5,i,iparm)
297 ecorr6=enetb(6,i,iparm)
298 eel_loc=enetb(7,i,iparm)
299 eello_turn3=enetb(8,i,iparm)
300 eello_turn4=enetb(9,i,iparm)
301 eturn6=enetb(10,i,iparm)
302 ebe=enetb(11,i,iparm)
303 escloc=enetb(12,i,iparm)
304 etors=enetb(13,i,iparm)
305 etors_d=enetb(14,i,iparm)
306 ehpb=enetb(15,i,iparm)
307 estr=enetb(18,i,iparm)
308 esccor=enetb(19,i,iparm)
309 edihcnstr=enetb(20,i,iparm)
310 ehomology_constr=enetb(22,i,iparm)
311 edfadis=enetb(23,i,iparm)
312 edfator=enetb(24,i,iparm)
313 edfanei=enetb(25,i,iparm)
314 edfabet=enetb(26,i,iparm)
316 write (iout,'(3i5,6f5.2,16f12.3)') i,ib,iparm,(ft(l),l=1,6),
317 & evdw+evdw_t,evdw2,ees,evdw1,ecorr,eel_loc,estr,ebe,escloc,
318 & etors,etors_d,eello_turn3,eello_turn4,esccor,
319 & ehomology_constr,edfadis,edfator,edfanei,edfabet
323 etot=wsc*(evdw+ft(6)*evdw_t)+wscp*evdw2+ft(1)*welec*ees
325 & +wang*ebe+ft(1)*wtor*etors+wscloc*escloc
326 & +wstrain*ehpb+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5
327 & +ft(5)*wcorr6*ecorr6+ft(3)*wturn4*eello_turn4
328 & +ft(2)*wturn3*eello_turn3
329 & +ft(5)*wturn6*eturn6+ft(2)*wel_loc*eel_loc
330 & +edihcnstr+ft(2)*wtor_d*etors_d+ft(1)*wsccor*esccor
331 & +wbond*estr+wdfa_dist*edfadis
332 & +wdfa_tor*edfator+wdfa_nei*edfanei+wdfa_beta*edfabet
334 etot=wsc*(evdw+ft(6)*evdw_t)+wscp*evdw2
335 & +ft(1)*welec*(ees+evdw1)
336 & +wang*ebe+ft(1)*wtor*etors+wscloc*escloc
337 & +wstrain*ehpb+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5
338 & +ft(5)*wcorr6*ecorr6+ft(3)*wturn4*eello_turn4
339 & +ft(2)*wturn3*eello_turn3
340 & +ft(5)*wturn6*eturn6+ft(2)*wel_loc*eel_loc+edihcnstr
341 & +ft(2)*wtor_d*etors_d+ft(1)*wsccor*esccor
342 & +wbond*estr+wdfa_dist*edfadis
343 & +wdfa_tor*edfator+wdfa_nei*edfanei+wdfa_beta*edfabet
346 write (iout,*) i,iparm,1.0d0/(beta_h(ib,iparm)*1.987D-3),
350 if (iparm.eq.1 .and. ib.eq.1) then
351 write (iout,*)"Conformation",i
354 energia(k)=enetb(k,i,iparm)
356 call enerprint(energia(0),fT)
360 write (iout,*) "homol_nset",homol_nset,nR(ib,iparm)
362 if (homol_nset.gt.1) then
365 Econstr=waga_homology(kk)*ehomology_constr
367 & -beta_h(ib,iparm)*(etot+Econstr)
369 write (iout,'(4i5,4e15.5)') i,kk,ib,iparm,
370 & etot,Econstr,v(i,kk,ib,iparm)
376 etot=etot+ehomology_constr
381 Econstr=Econstr+Kh(j,kk,ib,iparm)
382 & *(dd-q0(j,kk,ib,iparm))**2
385 & -beta_h(ib,iparm)*(etot+Econstr)
387 write (iout,'(4i5,4e15.5)') i,kk,ib,iparm,
388 & etot,v(i,kk,ib,iparm)
396 ! Simple iteration to calculate free energies corresponding to all simulation
400 ! Compute new free-energy values corresponding to the righ-hand side of the
401 ! equation and their derivatives.
402 write (iout,*) "------------------------fi"
413 vf=v(t,l,k,i)+f(l,k,i)
414 if (vf.gt.vmax) vmax=vf
422 aux=f(l,k,i)+v(t,l,k,i)-vmax
424 & denom=denom+snk(l,k,i,islice)*dexp(aux)
428 entfac(t)=-dlog(denom)-vmax
429 if (entfac(t).lt.entfac_min) entfac_min=entfac(t)
431 write (iout,*) t,"vmax",vmax," denom",denom,"entfac",entfac(t)
435 c write (iout,*) "entfac_min before AllReduce",entfac_min
436 c call MPI_AllReduce(entfac_min,entfac_min_t,1,
437 c & MPI_DOUBLE_PRECISION,MPI_MAX,WHAM_COMM,IERROR)
438 c entfac_min=entfac_min_t
439 c write (iout,*) "entfac_min after AllReduce",entfac_min
443 c entfac(t)=entfac(t)-entfac_min
446 c do t=1,ntot(islice)
447 c entfac(t)=entfac(t)-entfac_min
452 do ii=1,nR(iib,iparm)
454 fi_p_min(ii,iib,iparm)=-1.0d10
456 aux=v(t,ii,iib,iparm)+entfac(t)
457 if (aux.gt.fi_p_min(ii,iib,iparm))
458 & fi_p_min(ii,iib,iparm)=aux
462 aux=v(t,ii,iib,iparm)+entfac(t)
463 if (aux.gt.fi_min(ii,iib,iparm))
464 & fi_min(ii,iib,iparm)=aux
472 write (iout,*) "fi_min before AllReduce"
475 write (iout,*) (i,j,k,fi_p_min(k,j,i),k=1,nR(j,i))
479 call MPI_AllReduce(fi_p_min,fi_min,MaxR*MaxT_h*nParmSet,
480 & MPI_DOUBLE_PRECISION,MPI_MAX,WHAM_COMM,IERROR)
482 write (iout,*) "fi_min after AllReduce"
485 write (iout,*) (i,j,k,fi_min(k,j,i),k=1,nR(j,i))
492 do ii=1,nR(iib,iparm)
494 fi_p(ii,iib,iparm)=0.0d0
496 fi_p(ii,iib,iparm)=fi_p(ii,iib,iparm)
497 & +dexp(v(t,ii,iib,iparm)+entfac(t)-fi_min(ii,iib,iparm))
499 write (iout,'(4i5,4e15.5)') t,ii,iib,iparm,
500 & v(t,ii,iib,iparm),entfac(t),fi_min(ii,iib,iparm),
505 fi(ii,iib,iparm)=0.0d0
507 fi(ii,iib,iparm)=fi(ii,iib,iparm)
508 & +dexp(v(t,ii,iib,iparm)+entfac(t)-fi_min(ii,iib,iparm))
517 write (iout,*) "fi before MPI_Reduce me",me,' master',master
520 write (iout,*) "iparm",iparm," ib",ib
521 write (iout,*) "beta=",beta_h(ib,iparm)
522 write (iout,'(8e15.5)') (fi_p(i,ib,iparm),i=1,nR(ib,iparm))
527 write (iout,*) "REDUCE size",maxR,MaxT_h,nParmSet,
528 & maxR*MaxT_h*nParmSet
529 write (iout,*) "MPI_COMM_WORLD",MPI_COMM_WORLD,
530 & " WHAM_COMM",WHAM_COMM
532 call MPI_Reduce(fi_p(1,1,1),fi(1,1,1),maxR*MaxT_h*nParmSet,
533 & MPI_DOUBLE_PRECISION,
534 & MPI_SUM,Master,WHAM_COMM,IERROR)
536 write (iout,*) "fi after MPI_Reduce nparmset",nparmset
538 write (iout,*) "iparm",iparm
540 write (iout,*) "beta=",beta_h(ib,iparm)
541 write (iout,'(8e15.5)') (fi(i,ib,iparm),i=1,nR(ib,iparm))
545 if (me1.eq.Master) then
551 fi(i,ib,iparm)=-dlog(fi(i,ib,iparm))-fi_min(i,ib,iparm)
552 avefi=avefi+fi(i,ib,iparm)
558 write (iout,*) "Parameter set",iparm
560 write (iout,*) "beta=",beta_h(ib,iparm)
562 fi(i,ib,iparm)=fi(i,ib,iparm)-avefi
564 write (iout,'(6f15.5)') (fi(i,ib,iparm),i=1,nR(ib,iparm))
565 write (iout,'(6f15.5)') (f(i,ib,iparm),i=1,nR(ib,iparm))
569 ! Compute the norm of free-energy increments.
574 finorm=finorm+dabs(fi(i,ib,iparm)-f(i,ib,iparm))
575 f(i,ib,iparm)=fi(i,ib,iparm)
580 write (iout,*) 'Iteration',iter,' finorm',finorm
584 call MPI_Bcast(f(1,1,1),MaxR*MaxT_h*nParmSet,
585 & MPI_DOUBLE_PRECISION,Master,
587 call MPI_Bcast(finorm,1,MPI_DOUBLE_PRECISION,Master,
590 ! Exit, if the increment norm is smaller than pre-assigned tolerance.
591 if (finorm.lt.fimin) then
592 write (iout,*) 'Iteration converged'
599 ! Now, put together the histograms from all simulations, in order to get the
600 ! unbiased total histogram.
602 C Determine the minimum free energies
608 c write (9,'(3i5,f10.5)') i,(iparm,potE(i,iparm),iparm=1,nParmSet)
611 write (iout,'(2i5,21f8.2)') i,iparm,
612 & (enetb(k,i,iparm),k=1,22)
614 call restore_parm(iparm)
616 write (iout,*) wsc,wscp,welec,wvdwpp,wang,wtor,wscloc,
617 & wcorr,wcorr5,wcorr6,wturn4,wturn3,wturn6,wel_loc,
618 & wtor_d,wsccor,wbond
621 if (rescale_mode.eq.1) then
622 quot=1.0d0/(beta_h(ib,iparm)*1.987D-3*T0)
629 fT(l)=kfacl/(kfacl-1.0d0+quotl)
632 tt = 1.0d0/(beta_h(ib,iparm)*1.987D-3)
633 ft(6)=(320.0d0+80.0d0*dtanh((tt-320.0d0)/80.0d0))/320.0d0
635 ft(6)=1.0d0/(beta_h(ib,iparm)*1.987D-3*T0)
639 else if (rescale_mode.eq.2) then
640 quot=1.0d0/(T0*beta_h(ib,iparm)*1.987D-3)
644 fT(l)=1.12692801104297249644d0/
645 & dlog(dexp(quotl)+dexp(-quotl))
648 tt = 1.0d0/(beta_h(ib,iparm)*1.987D-3)
649 ft(6)=(320.0d0+80.0d0*dtanh((tt-320.0d0)/80.0d0))/320.0d0
651 ft(6)=1.0d0/(beta_h(ib,iparm)*1.987D-3*T0)
655 c write (iout,*) 1.0d0/(beta_h(ib,iparm)*1.987D-3),ft
656 else if (rescale_mode.eq.0) then
661 write (iout,*) "Error in WHAM_CALC: wrong RESCALE_MODE",
666 evdw=enetb(1,i,iparm)
667 evdw_t=enetb(21,i,iparm)
669 evdw2_14=enetb(17,i,iparm)
670 evdw2=enetb(2,i,iparm)+evdw2_14
672 evdw2=enetb(2,i,iparm)
677 evdw1=enetb(16,i,iparm)
682 ecorr=enetb(4,i,iparm)
683 ecorr5=enetb(5,i,iparm)
684 ecorr6=enetb(6,i,iparm)
685 eel_loc=enetb(7,i,iparm)
686 eello_turn3=enetb(8,i,iparm)
687 eello_turn4=enetb(9,i,iparm)
688 eturn6=enetb(10,i,iparm)
689 ebe=enetb(11,i,iparm)
690 escloc=enetb(12,i,iparm)
691 etors=enetb(13,i,iparm)
692 etors_d=enetb(14,i,iparm)
693 ehpb=enetb(15,i,iparm)
694 estr=enetb(18,i,iparm)
695 esccor=enetb(19,i,iparm)
696 edihcnstr=enetb(20,i,iparm)
697 ehomology_constr=enetb(22,i,iparm)
699 & ehomology_constr=waga_homology(ihset)*ehomology_constr
700 edfadis=enetb(23,i,iparm)
701 edfator=enetb(24,i,iparm)
702 edfanei=enetb(25,i,iparm)
703 edfabet=enetb(26,i,iparm)
705 write (iout,'(3i5,6f5.2,17f12.3)') i,ib,iparm,(ft(l),l=1,6),
706 & evdw+evdw_t,evdw2,ees,evdw1,ecorr,eel_loc,estr,ebe,escloc,
707 & etors,etors_d,eello_turn3,eello_turn4,esccor,edihcnstr,
708 & ehomology_constr+wdfa_dist*edfadis+wdfa_tor*edfator+
709 & wdfa_nei*edfanei+wdfa_beta*edfabet
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+ehomology_constr+wdfa_dist*edfadis
722 & +wdfa_tor*edfator+wdfa_nei*edfanei+wdfa_beta*edfabet
724 etot=wsc*(evdw+ft(6)*evdw_t)+wscp*evdw2
725 & +ft(1)*welec*(ees+evdw1)
726 & +wang*ebe+ft(1)*wtor*etors+wscloc*escloc
727 & +wstrain*ehpb+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5
728 & +ft(5)*wcorr6*ecorr6+ft(3)*wturn4*eello_turn4
729 & +ft(2)*wturn3*eello_turn3
730 & +ft(5)*wturn6*eturn6+ft(2)*wel_loc*eel_loc+edihcnstr
731 & +ft(2)*wtor_d*etors_d+ft(1)*wsccor*esccor
732 & +wbond*estr+ehomology_constr+wdfa_dist*edfadis
733 & +wdfa_tor*edfator+wdfa_nei*edfanei+wdfa_beta*edfabet
736 write (iout,*) "i",i," ib",ib,
737 & " temp",1.0d0/(1.987d-3*beta_h(ib,iparm))," etot",etot,
738 & " entfac",entfac(i)," ecorr",etot-entfac(i)/beta_h(ib,iparm)
740 etot=etot-entfac(i)/beta_h(ib,iparm)
741 if(etot.lt.potEmin_all(ib,iparm)) potEmin_all(ib,iparm)=etot
743 write (iout,*) "efree",etot," potEmin",potEmin_all(ib,iparm)
749 write (iout,*) "The potEmin array before reduction"
751 write (iout,*) "Parameter set",i
753 write (iout,*) j,PotEmin_all(j,i)
756 write (iout,*) "potEmin_min",potEmin_min
759 C Determine the minimum energes for all parameter sets and temperatures
760 call MPI_AllReduce(potEmin_all(1,1),potEmin_t_all(1,1),
761 & maxT_h*nParmSet,MPI_DOUBLE_PRECISION,MPI_MIN,WHAM_COMM,IERROR)
764 potEmin_all(j,i)=potEmin_t_all(j,i)
768 potEmin_min=potEmin_all(1,1)
771 if (potEmin_all(j,i).lt.potEmin_min)
772 & potEmin_min=potEmin_all(j,i)
776 write (iout,*) "The potEmin array"
778 write (iout,*) "Parameter set",i
780 write (iout,*) j,1.0d0/(1.987d-3*beta_h(j,i)),
784 write (iout,*) "potEmin_min",potEmin_min
796 write (iout,*) "--------------hist"
800 sumW_p(i,iparm)=0.0d0
801 sumE_p(i,iparm)=0.0d0
802 sumEbis_p(i,iparm)=0.0d0
803 sumEsq_p(i,iparm)=0.0d0
805 sumQ_p(j,i,iparm)=0.0d0
806 sumQsq_p(j,i,iparm)=0.0d0
807 sumEQ_p(j,i,iparm)=0.0d0
817 sumEbis(i,iparm)=0.0d0
818 sumEsq(i,iparm)=0.0d0
820 sumQ(j,i,iparm)=0.0d0
821 sumQsq(j,i,iparm)=0.0d0
822 sumEQ(j,i,iparm)=0.0d0
828 c 8/26/05 entropy distribution
833 c ent=-dlog(entfac(t))
835 if (ent.lt.entmin_p) entmin_p=ent
836 if (ent.gt.entmax_p) entmax_p=ent
838 write (iout,*) "entmin",entmin_p," entmax",entmax_p
840 call MPI_Allreduce(entmin_p,entmin,1,MPI_DOUBLE_PRECISION,MPI_MIN,
842 call MPI_Allreduce(entmax_p,entmax,1,MPI_DOUBLE_PRECISION,MPI_MAX,
844 ientmax=entmax-entmin
845 if (ientmax.gt.2000) ientmax=2000
846 write (iout,*) "entmin",entmin," entmax",entmax," ientmax",ientmax
849 c ient=-dlog(entfac(t))-entmin
850 ient=entfac(t)-entmin
851 if (ient.le.2000) histent_p(ient)=histent_p(ient)+1
853 call MPI_Allreduce(histent_p(0),histent(0),ientmax+1,MPI_INTEGER,
854 & MPI_SUM,WHAM_COMM,IERROR)
855 if (me1.eq.Master) then
856 write (iout,*) "Entropy histogram"
858 write(iout,'(f15.4,i10)') entmin+i,histent(i)
866 if (ent.lt.entmin) entmin=ent
867 if (ent.gt.entmax) entmax=ent
869 ientmax=-dlog(entmax)-entmin
870 if (ientmax.gt.2000) ientmax=2000
872 ient=entfac(t)-entmin
873 if (ient.le.2000) histent(ient)=histent(ient)+1
875 write (iout,*) "Entropy histogram"
877 write(iout,'(2f15.4)') entmin+i,histent(i)
882 c write (iout,*) "me1",me1," scount",scount(me1)
908 hrmsrgy(j,i,ib)=0.0d0
910 hrmsrgy_p(j,i,ib)=0.0d0
922 hfin_ent_p(ind)=hfin_ent_p(ind)+dexp(entfac(t))
924 hfin_ent(ind)=hfin_ent(ind)+dexp(entfac(t))
926 call restore_parm(iparm)
927 evdw=enetb(21,t,iparm)
928 evdw_t=enetb(1,t,iparm)
930 evdw2_14=enetb(17,t,iparm)
931 evdw2=enetb(2,t,iparm)+evdw2_14
933 evdw2=enetb(2,t,iparm)
938 evdw1=enetb(16,t,iparm)
943 ecorr=enetb(4,t,iparm)
944 ecorr5=enetb(5,t,iparm)
945 ecorr6=enetb(6,t,iparm)
946 eel_loc=enetb(7,t,iparm)
947 eello_turn3=enetb(8,t,iparm)
948 eello_turn4=enetb(9,t,iparm)
949 eturn6=enetb(10,t,iparm)
950 ebe=enetb(11,t,iparm)
951 escloc=enetb(12,t,iparm)
952 etors=enetb(13,t,iparm)
953 etors_d=enetb(14,t,iparm)
954 ehpb=enetb(15,t,iparm)
955 estr=enetb(18,t,iparm)
956 esccor=enetb(19,t,iparm)
957 edihcnstr=enetb(20,t,iparm)
958 ehomology_constr=enetb(22,t,iparm)
960 & ehomology_constr=waga_homology(ihset)*ehomology_constr
961 edfadis=enetb(23,t,iparm)
962 edfator=enetb(24,t,iparm)
963 edfanei=enetb(25,t,iparm)
964 edfabet=enetb(26,t,iparm)
967 betaT=startGridT+k*delta_T
971 if (rescale_mode.eq.1) then
979 denom=kfacl-1.0d0+quotl
981 ftprim(l)=-l*ft(l)*quotl1/(T0*denom)
982 ftbis(l)=l*kfacl*quotl1*
983 & (2*l*quotl-(l-1)*denom)/(quot*t0*t0*denom**3)
986 ft(6)=(320.0d0+80.0d0*dtanh((betaT-320.0d0)/80.0d0))/
988 ftprim(6)=1.0d0/(320.0d0*dcosh((betaT-320.0d0)/80.0d0)**2)
989 ftbis(6)=-2.0d0*dtanh((betaT-320.0d0)/80.0d0)
990 & /(320.0d0*80.0d0*dcosh((betaT-320.0d0)/80.0d0)**3)
1000 else if (rescale_mode.eq.2) then
1008 logfac=1.0d0/dlog(eplus+eminus)
1009 tanhT=(eplus-eminus)/(eplus+eminus)
1010 fT(l)=1.12692801104297249644d0*logfac
1011 ftprim(l)=-l*quotl1*ft(l)*tanhT*logfac/T0
1012 ftbis(l)=(l-1)*ftprim(l)/(quot*T0)-
1013 & 2*l*quotl1/T0*logfac*
1014 & (2*l*quotl1*ft(l)/(T0*(eplus+eminus)**2)
1018 ft(6)=(320.0d0+80.0d0*dtanh((betaT-320.0d0)/80.0d0))/
1020 ftprim(6)=1.0d0/(320.0d0*dcosh((betaT-320.0d0)/80.0d0)**2)
1021 ftbis(6)=-2.0d0*dtanh((betaT-320.0d0)/80.0d0)
1022 & /(320.0d0*80.0d0*dcosh((betaT-320.0d0)/80.0d0)**3)
1023 #elif defined(FUNCT)
1032 else if (rescale_mode.eq.0) then
1038 write (iout,*) "Error in WHAM_CALC: wrong RESCALE_MODE",
1043 c write (iout,*) "ftprim",ftprim
1044 c write (iout,*) "ftbis",ftbis
1045 betaT=1.0d0/(1.987D-3*betaT)
1046 if (betaT.ge.beta_h(1,iparm)) then
1047 potEmin=potEmin_all(1,iparm)+
1048 & (potEmin_all(1,iparm)-potEmin_all(2,iparm))/
1049 & (1.0/beta_h(1,iparm)-1.0/beta_h(2,iparm))*
1050 & (1.0/betaT-1.0/beta_h(1,iparm))
1052 write(iout,*) "first",temper,potEmin
1054 else if (betaT.le.beta_h(nT_h(iparm),iparm)) then
1055 potEmin=potEmin_all(nT_h(iparm),iparm)+
1056 &(potEmin_all(nT_h(iparm),iparm)-potEmin_all(nT_h(iparm)-1,iparm))/
1057 &(1.0/beta_h(nT_h(iparm),iparm)-1.0/beta_h(nT_h(iparm)-1,iparm))*
1058 &(1.0/betaT-1.0/beta_h(nt_h(iparm),iparm))
1060 write (iout,*) "last",temper,potEmin
1063 do l=1,nT_h(iparm)-1
1064 if (betaT.le.beta_h(l,iparm) .and.
1065 & betaT.gt.beta_h(l+1,iparm)) then
1066 potEmin=potEmin_all(l,iparm)
1068 write (iout,*) "l",l,
1069 & betaT,1.0d0/(1.987D-3*beta_h(l,iparm)),
1070 & 1.0d0/(1.987D-3*beta_h(l+1,iparm)),temper,potEmin
1077 write (iout,*) "k",k," potEmin",potEmin
1080 etot=wsc*(evdw+ft(6)*evdw_t)+wscp*evdw2+ft(1)*welec*ees
1082 & +wang*ebe+ft(1)*wtor*etors+wscloc*escloc
1083 & +wstrain*ehpb+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5
1084 & +ft(5)*wcorr6*ecorr6+ft(3)*wturn4*eello_turn4
1085 & +ft(2)*wturn3*eello_turn3
1086 & +ft(5)*wturn6*eturn6+ft(2)*wel_loc*eel_loc
1087 & +edihcnstr+ft(2)*wtor_d*etors_d+ft(1)*wsccor*esccor
1088 & +wbond*estr+ehomology_constr
1089 eprim=ftprim(6)*evdw_t+ftprim(1)*welec*ees
1090 & +ftprim(1)*wtor*etors+
1091 & ftprim(3)*wcorr*ecorr+ftprim(4)*wcorr5*ecorr5+
1092 & ftprim(5)*wcorr6*ecorr6+ftprim(3)*wturn4*eello_turn4+
1093 & ftprim(2)*wturn3*eello_turn3+ftprim(5)*wturn6*eturn6+
1094 & ftprim(2)*wel_loc*eel_loc+ftprim(2)*wtor_d*etors_d+
1095 & ftprim(1)*wsccor*esccor
1096 ebis=ftbis(1)*welec*ees+ftbis(1)*wtor*etors+
1097 & ftbis(3)*wcorr*ecorr+ftbis(4)*wcorr5*ecorr5+
1098 & ftbis(5)*wcorr6*ecorr6+ftbis(3)*wturn4*eello_turn4+
1099 & ftbis(2)*wturn3*eello_turn3+ftbis(5)*wturn6*eturn6+
1100 & ftbis(2)*wel_loc*eel_loc+ftbis(2)*wtor_d*etors_d+
1101 & ftbis(1)*wsccor*esccor
1103 etot=wsc*(evdw+ft(6)*evdw_t)+wscp*evdw2
1104 & +ft(1)*welec*(ees+evdw1)
1105 & +wang*ebe+ft(1)*wtor*etors+wscloc*escloc
1106 & +wstrain*ehpb+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5
1107 & +ft(5)*wcorr6*ecorr6+ft(3)*wturn4*eello_turn4
1108 & +ft(2)*wturn3*eello_turn3
1109 & +ft(5)*wturn6*eturn6+ft(2)*wel_loc*eel_loc+edihcnstr
1110 & +ft(2)*wtor_d*etors_d+ft(1)*wsccor*esccor
1111 & +wbond*estr+ehomology_constr
1112 eprim=ftprim(6)*evdw_t+ftprim(1)*welec*(ees+evdw1)
1113 & +ftprim(1)*wtor*etors+
1114 & ftprim(3)*wcorr*ecorr+ftprim(4)*wcorr5*ecorr5+
1115 & ftprim(5)*wcorr6*ecorr6+ftprim(3)*wturn4*eello_turn4+
1116 & ftprim(2)*wturn3*eello_turn3+ftprim(5)*wturn6*eturn6+
1117 & ftprim(2)*wel_loc*eel_loc+ftprim(2)*wtor_d*etors_d+
1118 & ftprim(1)*wsccor*esccor
1119 ebis=ftbis(1)*welec*(ees+evdw1)+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 & ftprim(1)*wsccor*esccor
1126 weight=dexp(-betaT*(etot-potEmin)+entfac(t))
1128 write (iout,'(3i5,6f5.2,17f12.3)') i,ib,iparm,(ft(l),l=1,6),
1129 & evdw+evdw_t,evdw2,ees,evdw1,ecorr,eel_loc,estr,ebe,escloc,
1130 & etors,etors_d,eello_turn3,eello_turn4,esccor,edihcnstr,
1131 & ehomology_constr+wdfa_dist*edfadis+wdfa_tor*edfator+
1132 & wdfa_nei*edfanei+wdfa_beta*edfabet
1133 write (iout,*) "iparm",iparm," t",t," temper",temper,
1134 & " etot",etot," entfac",entfac(t),
1135 & " efree",etot-entfac(t)/betaT," potEmin",potEmin,
1136 & " boltz",-betaT*(etot-potEmin)+entfac(t),
1137 & " weight",weight," ebis",ebis
1139 etot=etot-temper*eprim
1141 sumW_p(k,iparm)=sumW_p(k,iparm)+weight
1142 sumE_p(k,iparm)=sumE_p(k,iparm)+etot*weight
1143 sumEbis_p(k,iparm)=sumEbis_p(k,iparm)+ebis*weight
1144 sumEsq_p(k,iparm)=sumEsq_p(k,iparm)+etot**2*weight
1146 sumQ_p(j,k,iparm)=sumQ_p(j,k,iparm)+q(j,t)*weight
1147 sumQsq_p(j,k,iparm)=sumQsq_p(j,k,iparm)+q(j,t)**2*weight
1148 sumEQ_p(j,k,iparm)=sumEQ_p(j,k,iparm)
1149 & +etot*q(j,t)*weight
1152 sumW(k,iparm)=sumW(k,iparm)+weight
1153 sumE(k,iparm)=sumE(k,iparm)+etot*weight
1154 sumEbis(k,iparm)=sumEbis(k,iparm)+ebis*weight
1155 sumEsq(k,iparm)=sumEsq(k,iparm)+etot**2*weight
1157 sumQ(j,k,iparm)=sumQ(j,k,iparm)+q(j,t)*weight
1158 sumQsq(j,k,iparm)=sumQsq(j,k,iparm)+q(j,t)**2*weight
1159 sumEQ(j,k,iparm)=sumEQ(j,k,iparm)
1160 & +etot*q(j,t)*weight
1164 indE = aint(potE(t,iparm)-aint(potEmin))
1165 if (indE.ge.0 .and. indE.le.maxinde) then
1166 if (indE.gt.upindE_p) upindE_p=indE
1167 histE_p(indE)=histE_p(indE)+dexp(-entfac(t))
1171 potEmin=potEmin_all(ib,iparm)
1172 expfac=dexp(-beta_h(ib,iparm)*(etot-potEmin)+entfac(t))
1173 hfin_p(ind,ib)=hfin_p(ind,ib)+
1174 & dexp(-beta_h(ib,iparm)*(etot-potEmin)+entfac(t))
1176 indrgy=dint((q(nQ+2,t)-rgymin)/deltrgy)
1177 indrms=dint((q(nQ+1,t)-rmsmin)/deltrms)
1178 hrmsrgy_p(indrgy,indrms,ib)=
1179 & hrmsrgy_p(indrgy,indrms,ib)+expfac
1184 potEmin=potEmin_all(ib,iparm)
1185 expfac=dexp(-beta_h(ib,iparm)*(etot-potEmin)+entfac(t))
1186 hfin(ind,ib)=hfin(ind,ib)+
1187 & dexp(-beta_h(ib,iparm)*(etot-potEmin)+entfac(t))
1189 indrgy=dint((q(nQ+2,t)-rgymin)/deltrgy)
1190 indrms=dint((q(nQ+1,t)-rmsmin)/deltrms)
1191 hrmsrgy(indrgy,indrms,ib)=
1192 & hrmsrgy(indrgy,indrms,ib)+expfac
1198 if (histout) call MPI_Reduce(hfin_p(0,ib),hfin(0,ib),nbin,
1199 & MPI_DOUBLE_PRECISION,MPI_SUM,Master,WHAM_COMM,IERROR)
1201 call MPI_Reduce(hrmsrgy_p(0,0,ib),hrmsrgy(0,0,ib),
1202 & (MaxBinRgy+1)*(nbin_rms+1),MPI_DOUBLE_PRECISION,MPI_SUM,Master,
1206 call MPI_Reduce(upindE_p,upindE,1,
1207 & MPI_INTEGER,MPI_MAX,Master,WHAM_COMM,IERROR)
1208 call MPI_Reduce(histE_p(0),histE(0),maxindE,
1209 & MPI_DOUBLE_PRECISION,MPI_SUM,Master,WHAM_COMM,IERROR)
1211 if (me1.eq.master) then
1215 write (iout,'(6x,$)')
1216 write (iout,'(f20.2,$)') (1.0d0/(1.987D-3*beta_h(ib,iparm)),
1220 write (iout,'(/a)') 'Final histograms'
1222 if (nslice.eq.1) then
1223 if (separate_parset) then
1224 write(licz3,"(bz,i3.3)") myparm
1225 histname=prefix(:ilen(prefix))//'_par'//licz3//'.hist'
1227 histname=prefix(:ilen(prefix))//'.hist'
1230 if (separate_parset) then
1231 write(licz3,"(bz,i3.3)") myparm
1232 histname=prefix(:ilen(prefix))//'_par'//licz3//
1233 & '_slice_'//licz2//'.hist'
1235 histname=prefix(:ilen(prefix))//'_slice_'//licz2//'.hist'
1238 #if defined(AIX) || defined(PGI)
1239 open (ihist,file=histname,position='append')
1241 open (ihist,file=histname,access='append')
1249 sumH=sumH+hfin(t,ib)
1251 if (sumH.gt.0.0d0) then
1253 jj = mod(liczba,nbin1)
1255 write (iout,'(f6.3,$)') dmin+(jj+0.5d0)*delta
1257 & write (ihist,'(f6.3,$)') dmin+(jj+0.5d0)*delta
1260 write (iout,'(e20.10,$)') hfin(t,ib)
1261 if (histfile) write (ihist,'(e20.10,$)') hfin(t,ib)
1263 write (iout,'(i5)') iparm
1264 if (histfile) write (ihist,'(i5)') iparm
1271 if (nslice.eq.1) then
1272 if (separate_parset) then
1273 write(licz3,"(bz,i3.3)") myparm
1274 histname=prefix(:ilen(prefix))//"_par"//licz3//'.ent'
1276 histname=prefix(:ilen(prefix))//'.ent'
1279 if (separate_parset) then
1280 write(licz3,"(bz,i3.3)") myparm
1281 histname=prefix(:ilen(prefix))//'par_'//licz3//
1282 & '_slice_'//licz2//'.ent'
1284 histname=prefix(:ilen(prefix))//'_slice_'//licz2//'.ent'
1287 #if defined(AIX) || defined(PGI)
1288 open (ihist,file=histname,position='append')
1290 open (ihist,file=histname,access='append')
1292 write (ihist,'(a)') "# Microcanonical entropy"
1294 write (ihist,'(f8.0,$)') dint(potEmin)+i
1295 if (histE(i).gt.0.0e0) then
1296 write (ihist,'(f15.5,$)') dlog(histE(i))
1298 write (ihist,'(f15.5,$)') 0.0d0
1304 write (iout,*) "Microcanonical entropy"
1306 write (iout,'(f8.0,$)') dint(potEmin)+i
1307 if (histE(i).gt.0.0e0) then
1308 write (iout,'(f15.5,$)') dlog(histE(i))
1310 write (iout,'(f15.5,$)') 0.0d0
1315 if (nslice.eq.1) then
1316 if (separate_parset) then
1317 write(licz3,"(bz,i3.3)") myparm
1318 histname=prefix(:ilen(prefix))//'_par'//licz3//'.rmsrgy'
1320 histname=prefix(:ilen(prefix))//'.rmsrgy'
1323 if (separate_parset) then
1324 write(licz3,"(bz,i3.3)") myparm
1325 histname=prefix(:ilen(prefix))//'_par'//licz3//
1326 & '_slice_'//licz2//'.rmsrgy'
1328 histname=prefix(:ilen(prefix))//'_slice_'//licz2//'.rmsrgy'
1331 #if defined(AIX) || defined(PGI)
1332 open (ihist,file=histname,position='append')
1334 open (ihist,file=histname,access='append')
1338 write(ihist,'(2f8.2,$)')
1339 & rgymin+deltrgy*j,rmsmin+deltrms*i
1341 if (hrmsrgy(j,i,ib).gt.0.0d0) then
1342 write(ihist,'(e14.5,$)')
1343 & -dlog(hrmsrgy(j,i,ib))/beta_h(ib,iparm)
1346 write(ihist,'(e14.5,$)') 1.0d6
1349 write (ihist,'(i2)') iparm
1357 call MPI_Reduce(hfin_ent_p(0),hfin_ent(0),nbin,
1358 & MPI_DOUBLE_PRECISION,MPI_SUM,Master,WHAM_COMM,IERROR)
1359 call MPI_Reduce(sumW_p(0,1),sumW(0,1),(nGridT+1)*nParmSet,
1360 & MPI_DOUBLE_PRECISION,MPI_SUM,Master,WHAM_COMM,IERROR)
1361 call MPI_Reduce(sumE_p(0,1),sumE(0,1),(nGridT+1)*nParmSet,
1362 & MPI_DOUBLE_PRECISION,MPI_SUM,Master,WHAM_COMM,IERROR)
1363 call MPI_Reduce(sumEbis_p(0,1),sumEbis(0,1),(nGridT+1)*nParmSet,
1364 & MPI_DOUBLE_PRECISION,MPI_SUM,Master,WHAM_COMM,IERROR)
1365 call MPI_Reduce(sumEsq_p(0,1),sumEsq(0,1),(nGridT+1)*nParmSet,
1366 & MPI_DOUBLE_PRECISION,MPI_SUM,Master,WHAM_COMM,IERROR)
1367 call MPI_Reduce(sumQ_p(1,0,1),sumQ(1,0,1),
1368 & MaxQ1*(nGridT+1)*nParmSet,MPI_DOUBLE_PRECISION,MPI_SUM,Master,
1370 call MPI_Reduce(sumQsq_p(1,0,1),sumQsq(1,0,1),
1371 & MaxQ1*(nGridT+1)*nParmSet,MPI_DOUBLE_PRECISION,MPI_SUM,Master,
1373 call MPI_Reduce(sumEQ_p(1,0,1),sumEQ(1,0,1),
1374 & MaxQ1*(nGridT+1)*nParmSet,MPI_DOUBLE_PRECISION,MPI_SUM,Master,
1376 if (me.eq.master) then
1378 write (iout,'(/a)') 'Thermal characteristics of folding'
1379 if (nslice.eq.1) then
1382 nazwa=prefix(:ilen(prefix))//"_slice_"//licz2
1385 if (nparmset.eq.1 .and. .not.separate_parset) then
1386 nazwa=nazwa(:iln)//".thermal"
1387 else if (nparmset.eq.1 .and. separate_parset) then
1388 write(licz3,"(bz,i3.3)") myparm
1389 nazwa=nazwa(:iln)//"_par_"//licz3//".thermal"
1392 if (nparmset.gt.1) then
1393 write(licz3,"(bz,i3.3)") iparm
1394 nazwa=nazwa(:iln)//"_par_"//licz3//".thermal"
1397 if (separate_parset) then
1398 write (iout,'(a,i3)') "Parameter set",myparm
1400 write (iout,'(a,i3)') "Parameter set",iparm
1404 betaT=1.0d0/(1.987D-3*(startGridT+i*delta_T))
1405 if (betaT.ge.beta_h(1,iparm)) then
1406 potEmin=potEmin_all(1,iparm)
1407 else if (betaT.lt.beta_h(nT_h(iparm),iparm)) then
1408 potEmin=potEmin_all(nT_h(iparm),iparm)
1410 do l=1,nT_h(iparm)-1
1411 if (betaT.le.beta_h(l,iparm) .and.
1412 & betaT.gt.beta_h(l+1,iparm)) then
1413 potEmin=potEmin_all(l,iparm)
1419 c write (iout,*) "i",i," potEmin",potEmin
1421 sumE(i,iparm)=sumE(i,iparm)/sumW(i,iparm)
1422 sumEbis(i,iparm)=(startGridT+i*delta_T)*sumEbis(i,iparm)/
1424 sumEsq(i,iparm)=(sumEsq(i,iparm)/sumW(i,iparm)
1425 & -sumE(i,iparm)**2)/(1.987D-3*(startGridT+i*delta_T)**2)
1427 sumQ(j,i,iparm)=sumQ(j,i,iparm)/sumW(i,iparm)
1428 sumQsq(j,i,iparm)=sumQsq(j,i,iparm)/sumW(i,iparm)
1429 & -sumQ(j,i,iparm)**2
1430 sumEQ(j,i,iparm)=sumEQ(j,i,iparm)/sumW(i,iparm)
1431 & -sumQ(j,i,iparm)*sumE(i,iparm)
1433 sumW(i,iparm)=-dlog(sumW(i,iparm))*(1.987D-3*
1434 & (startGridT+i*delta_T))+potEmin
1435 write (iout,'(f7.1,2f15.5,$)') startGridT+i*delta_T,
1436 & sumW(i,iparm),sumE(i,iparm)
1437 write (iout,'(f10.5,$)') (sumQ(j,i,iparm),j=1,nQ+2)
1438 write (iout,'(e15.5,$)') sumEsq(i,iparm)-sumEbis(i,iparm),
1439 & (sumQsq(j,i,iparm),j=1,nQ+2),(sumEQ(j,i,iparm),j=1,nQ+2)
1441 write (34,'(f7.1,2f15.5,$)') startGridT+i*delta_T,
1442 & sumW(i,iparm),sumE(i,iparm)
1443 write (34,'(f10.5,$)') (sumQ(j,i,iparm),j=1,nQ+2)
1444 write (34,'(e15.5,$)') sumEsq(i,iparm)-sumEbis(i,iparm),
1445 & (sumQsq(j,i,iparm),j=1,nQ+2),(sumEQ(j,i,iparm),j=1,nQ+2)
1452 if (hfin_ent(t).gt.0.0d0) then
1454 jj = mod(liczba,nbin1)
1455 write (iout,'(f6.3,e20.10," ent")') dmin+(jj+0.5d0)*delta,
1457 if (histfile) write (ihist,'(f6.3,e20.10," ent")')
1458 & dmin+(jj+0.5d0)*delta,
1462 if (histfile) close(ihist)
1466 ! Write data for zscore
1467 if (nslice.eq.1) then
1468 zscname=prefix(:ilen(prefix))//".zsc"
1470 zscname=prefix(:ilen(prefix))//"_slice_"//licz2//".zsc"
1472 #if defined(AIX) || defined(PGI)
1473 open (izsc,file=prefix(:ilen(prefix))//'.zsc',position='append')
1475 open (izsc,file=prefix(:ilen(prefix))//'.zsc',access='append')
1477 write (izsc,'("NQ=",i1," NPARM=",i1)') nQ,nParmSet
1479 write (izsc,'("NT=",i1)') nT_h(iparm)
1481 write (izsc,'("TEMP=",f6.1," NR=",i2," SNK=",$)')
1482 & 1.0d0/(beta_h(ib,iparm)*1.987D-3),nR(ib,iparm)
1483 jj = min0(nR(ib,iparm),7)
1484 write (izsc,'(i8,$)') (snk(i,ib,iparm,islice),i=1,jj)
1485 write (izsc,'(a1,$)') (" ",i=22+8*jj+1,79)
1486 write (izsc,'("&")')
1487 if (nR(ib,iparm).gt.7) then
1488 do ii=8,nR(ib,iparm),9
1489 jj = min0(nR(ib,iparm),ii+8)
1490 write (izsc,'(i8,$)') (snk(i,ib,iparm,islice),i=ii,jj)
1491 write (izsc,'(a1,$') (" ",i=(jj-ii+1)*8+1,79)
1492 write (izsc,'("&")')
1495 write (izsc,'("FI=",$)')
1496 jj=min0(nR(ib,iparm),7)
1497 write (izsc,'(f10.5,$)') (fi(i,ib,iparm),i=1,jj)
1498 write (izsc,'(a1,$)') (" ",i=3+10*jj+1,79)
1499 write (izsc,'("&")')
1500 if (nR(ib,iparm).gt.7) then
1501 do ii=8,nR(ib,iparm),9
1502 jj = min0(nR(ib,iparm),ii+8)
1503 write (izsc,'(f10.5,$)') (fi(i,ib,iparm),i=ii,jj)
1504 if (jj.eq.nR(ib,iparm)) then
1507 write (izsc,'(a1,$)') (" ",i=10*(jj-ii+1)+1,79)
1508 write (izsc,'(t80,"&")')
1513 write (izsc,'("KH=",$)')
1514 write (izsc,'(f7.2,$)') (Kh(j,i,ib,iparm),j=1,nQ)
1515 write (izsc,'(" Q0=",$)')
1516 write (izsc,'(f7.5,$)') (q0(j,i,ib,iparm),j=1,nQ)