double precision energia(0:max_ene)
#ifdef MPI
integer tmax_t,upindE_p
- double precision fi_p(MaxR,MaxT_h,Max_Parm)
+ double precision fi_p(MaxR,MaxT_h,Max_Parm),
+ & fi_p_min(MaxR,MaxT_h,Max_Parm)
double precision sumW_p(0:nGridT,Max_Parm),
& sumE_p(0:nGridT,Max_Parm),sumEsq_p(0:nGridT,Max_Parm),
& sumQ_p(MaxQ1,0:nGridT,Max_Parm),
& hfin_ent_p(0:MaxHdim),histE_p(0:maxindE),sumH,
& hrmsrgy_p(0:MaxBinRgy,0:MaxBinRms,maxT_h)
double precision rgymin_t,rmsmin_t,rgymax_t,rmsmax_t
- double precision potEmin_t,entmin_p,entmax_p
+ double precision potEmin_t,entmin_p,entmax_p,
+ & potEmin_t_all(maxT_h,Max_Parm)
integer histent_p(0:2000)
logical lprint /.true./
#endif
& sumEprim(0:NGridT,Max_Parm),sumEbis(0:NGridT,Max_Parm),betaT,
& weight,econstr
double precision fi(MaxR,maxT_h,Max_Parm),
+ & fi_min(MaxR,maxT_h,Max_Parm),
& dd,dd1,dd2,hh,dmin,denom,finorm,avefi,pom,
& hfin(0:MaxHdim,maxT_h),histE(0:maxindE),
& hrmsrgy(0:MaxBinRgy,0:MaxBinRms,maxT_h),
- & potEmin,ent,
- & hfin_ent(0:MaxHdim),vmax,aux
+ & potEmin,ent,potEmin_all(maxT_h,Max_Parm),potEmin_min,
+ & hfin_ent(0:MaxHdim),vmax,aux,entfac_min
double precision fT(6),fTprim(6),fTbis(6),quot,quotl1,quotl,kfacl,
& eprim,ebis,temper,kfac/2.4d0/,T0/300.0d0/,startGridT/200.0d0/,
& eplus,eminus,logfac,tanhT,tt
double precision etot,evdw,evdw_t,evdw2,ees,evdw1,ebe,etors,
& escloc,ehpb,ecorr,ecorr5,ecorr6,eello_turn4,eello_turn3,
- & eturn6,eel_loc,edihcnstr,etors_d,estr,evdw2_14,esccor
+ & eturn6,eel_loc,edihcnstr,etors_d,estr,evdw2_14,esccor,
+ & ehomology_constr,esaxs
+
integer ind_point(maxpoint),upindE,indE
character*16 plik
character*128 nazwa
integer ilen
external ilen
-
write(licz2,'(bz,i2.2)') islice
nbin1 = 1.0d0/delta
write (iout,'(//80(1h-)/"Solving WHAM equations for slice",
dmin=0.0d0
tmax=0
potEmin=1.0d10
+ do i=1,nParmset
+ do j=1,nT_h(i)
+ potEmin_all(j,i)=1.0d10
+ enddo
+ enddo
rgymin=1.0d10
rmsmin=1.0d10
rgymax=0.0d0
do t=0,MaxN
htot(t)=0
enddo
+C#define DEBUG
#ifdef MPI
do i=1,scount(me1)
#else
do i=1,ntot(islice)
#endif
- do j=1,nParmSet
- if (potE(i,j).le.potEmin) potEmin=potE(i,j)
- enddo
+C do j=1,nParmSet
+C if (potE(i,j).le.potEmin) potEmin=potE(i,j)
+C enddo
if (q(nQ+1,i).lt.rmsmin) rmsmin=q(nQ+1,i)
if (q(nQ+1,i).gt.rmsmax) rmsmax=q(nQ+1,i)
if (q(nQ+2,i).lt.rgymin) rgymin=q(nQ+2,i)
call MPI_AllReduce(tmax,tmax_t,1,MPI_INTEGER,MPI_MAX,
& WHAM_COMM,IERROR)
tmax=tmax_t
- call MPI_AllReduce(potEmin,potEmin_t,1,MPI_DOUBLE_PRECISION,
- & MPI_MIN,WHAM_COMM,IERROR)
+C call MPI_AllReduce(potEmin,potEmin_t,1,MPI_DOUBLE_PRECISION,
+C & MPI_MIN,WHAM_COMM,IERROR)
call MPI_AllReduce(rmsmin,rmsmin_t,1,MPI_DOUBLE_PRECISION,
& MPI_MIN,WHAM_COMM,IERROR)
call MPI_AllReduce(rmsmax,rmsmax_t,1,MPI_DOUBLE_PRECISION,
& MPI_MIN,WHAM_COMM,IERROR)
call MPI_AllReduce(rgymax,rgymax_t,1,MPI_DOUBLE_PRECISION,
& MPI_MAX,WHAM_COMM,IERROR)
- potEmin=potEmin_t/2
+C potEmin=potEmin_t/2
+
rgymin=rgymin_t
rgymax=rgymax_t
rmsmin=rmsmin_t
do iparm=1,nParmSet
#ifdef DEBUG
write (iout,'(2i5,21f8.2)') i,iparm,
- & (enetb(k,i,iparm),k=1,21)
+ & (enetb(k,i,iparm),k=1,26)
#endif
call restore_parm(iparm)
#ifdef DEBUG
write (iout,*) wsc,wscp,welec,wvdwpp,wang,wtor,wscloc,
& wcorr,wcorr5,wcorr6,wturn4,wturn3,wturn6,wel_loc,
- & wtor_d,wsccor,wbond
+ & wtor_d,wsccor,wbond,wsaxs
#endif
do ib=1,nT_h(iparm)
if (rescale_mode.eq.1) then
estr=enetb(18,i,iparm)
esccor=enetb(19,i,iparm)
edihcnstr=enetb(20,i,iparm)
+ ehomology_constr=enetb(22,i,iparm)
+ esaxs=enetb(26,i,iparm)
#ifdef DEBUG
write (iout,'(3i5,6f5.2,14f12.3)') i,ib,iparm,(ft(l),l=1,6),
& evdw+evdw_t,evdw2,ees,evdw1,ecorr,eel_loc,estr,ebe,escloc,
- & etors,etors_d,eello_turn3,eello_turn4,esccor
+ & etors,etors_d,eello_turn3,eello_turn4,esccor,esaxs
#endif
#ifdef SPLITELE
& +ft(2)*wturn3*eello_turn3
& +ft(5)*wturn6*eturn6+ft(2)*wel_loc*eel_loc
& +edihcnstr+ft(2)*wtor_d*etors_d+ft(1)*wsccor*esccor
- & +wbond*estr
+ & +wbond*estr+wsaxs*esaxs
#else
etot=wsc*(evdw+ft(6)*evdw_t)+wscp*evdw2
& +ft(1)*welec*(ees+evdw1)
& +ft(2)*wturn3*eello_turn3
& +ft(5)*wturn6*eturn6+ft(2)*wel_loc*eel_loc+edihcnstr
& +ft(2)*wtor_d*etors_d+ft(1)*wsccor*esccor
- & +wbond*estr
+ & +wbond*estr+wsaxs*esaxs
#endif
#ifdef DEBUG
write (iout,*) i,iparm,1.0d0/(beta_h(ib,iparm)*1.987D-3),
call enerprint(energia(0),fT)
endif
#endif
+#ifdef DEBUG
+ write (iout,*) "homol_nset",homol_nset,nR(ib,iparm)
+#endif
+ if (homol_nset.gt.1) then
+
+ do kk=1,nR(ib,iparm)
+ Econstr=waga_homology(kk)*ehomology_constr
+ v(i,kk,ib,iparm)=
+ & -beta_h(ib,iparm)*(etot+Econstr)
+#ifdef DEBUG
+ write (iout,'(4i5,4e15.5)') i,kk,ib,iparm,
+ & etot,Econstr,v(i,kk,ib,iparm)
+#endif
+ enddo ! kk
+
+ else
+
+ etot=etot+ehomology_constr
do kk=1,nR(ib,iparm)
Econstr=0.0d0
do j=1,nQ
& *(dd-q0(j,kk,ib,iparm))**2
enddo
v(i,kk,ib,iparm)=
- & -beta_h(ib,iparm)*(etot-potEmin+Econstr)
+ & -beta_h(ib,iparm)*(etot+Econstr)
#ifdef DEBUG
write (iout,'(4i5,4e15.5)') i,kk,ib,iparm,
& etot,potEmin,etot-potEmin,v(i,kk,ib,iparm)
#endif
enddo ! kk
+
+ endif
enddo ! ib
enddo ! iparm
enddo ! i
! Compute new free-energy values corresponding to the righ-hand side of the
! equation and their derivatives.
write (iout,*) "------------------------fi"
+ entfac_min=1.0d10
#ifdef MPI
do t=1,scount(me1)
#else
enddo
enddo
entfac(t)=-dlog(denom)-vmax
+ if (entfac(t).lt.entfac_min) entfac_min=entfac(t)
#ifdef DEBUG
write (iout,*) t,"vmax",vmax," denom",denom,"entfac",entfac(t)
#endif
enddo
+
+ do iparm=1,nParmSet
+ do iib=1,nT_h(iparm)
+ do ii=1,nR(iib,iparm)
+#ifdef MPI
+ fi_p_min(ii,iib,iparm)=-1.0d10
+ do t=1,scount(me)
+ aux=v(t,ii,iib,iparm)+entfac(t)
+ if (aux.gt.fi_p_min(ii,iib,iparm))
+ & fi_p_min(ii,iib,iparm)=aux
+ enddo
+#else
+ do t=1,ntot(islice)
+ aux=v(t,ii,iib,iparm)+entfac(t)
+ if (aux.gt.fi_min(ii,iib,iparm))
+ & fi_min(ii,iib,iparm)=aux
+ enddo
+#endif
+ enddo ! ii
+ enddo ! iib
+ enddo ! iparm
+#ifdef MPI
+#ifdef DEBUG
+ write (iout,*) "fi_min before AllReduce"
+ do i=1,nParmSet
+ do j=1,nT_h(i)
+ write (iout,*) (i,j,k,fi_p_min(k,j,i),k=1,nR(j,i))
+ enddo
+ enddo
+#endif
+ call MPI_AllReduce(fi_p_min,fi_min,MaxR*MaxT_h*nParmSet,
+ & MPI_DOUBLE_PRECISION,MPI_MAX,WHAM_COMM,IERROR)
+#ifdef DEBUG
+ write (iout,*) "fi_min after AllReduce"
+ do i=1,nParmSet
+ do j=1,nT_h(i)
+ write (iout,*) (i,j,k,fi_min(k,j,i),k=1,nR(j,i))
+ enddo
+ enddo
+#endif
+#endif
do iparm=1,nParmSet
do iib=1,nT_h(iparm)
do ii=1,nR(iib,iparm)
fi_p(ii,iib,iparm)=0.0d0
do t=1,scount(me)
fi_p(ii,iib,iparm)=fi_p(ii,iib,iparm)
- & +dexp(v(t,ii,iib,iparm)+entfac(t))
+ & +dexp(v(t,ii,iib,iparm)+entfac(t)-fi_min(ii,iib,iparm))
#ifdef DEBUG
- write (iout,'(4i5,3e15.5)') t,ii,iib,iparm,
- & v(t,ii,iib,iparm),entfac(t),fi_p(ii,iib,iparm)
+ write (iout,'(4i5,4e15.5)') t,ii,iib,iparm,
+ & v(t,ii,iib,iparm),entfac(t),fi_min(ii,iib,iparm),
+ & fi_p(ii,iib,iparm)
#endif
enddo
#else
fi(ii,iib,iparm)=0.0d0
do t=1,ntot(islice)
fi(ii,iib,iparm)=fi(ii,iib,iparm)
- & +dexp(v(t,ii,iib,iparm)+entfac(t))
+ & +dexp(v(t,ii,iib,iparm)+entfac(t)-fi_min(ii,iib,iparm))
enddo
#endif
enddo ! ii
enddo
enddo
#endif
+#ifdef DEBUG
write (iout,*) "REDUCE size",maxR,MaxT_h,nParmSet,
& maxR*MaxT_h*nParmSet
write (iout,*) "MPI_COMM_WORLD",MPI_COMM_WORLD,
& " WHAM_COMM",WHAM_COMM
+#endif
call MPI_Reduce(fi_p(1,1,1),fi(1,1,1),maxR*MaxT_h*nParmSet,
& MPI_DOUBLE_PRECISION,
& MPI_SUM,Master,WHAM_COMM,IERROR)
do iparm=1,nParmSet
do ib=1,nT_h(iparm)
do i=1,nR(ib,iparm)
- fi(i,ib,iparm)=-dlog(fi(i,ib,iparm))
+ fi(i,ib,iparm)=-dlog(fi(i,ib,iparm))-fi_min(i,ib,iparm)
avefi=avefi+fi(i,ib,iparm)
enddo
enddo
20 continue
! Now, put together the histograms from all simulations, in order to get the
! unbiased total histogram.
+
+C Determine the minimum free energies
+#ifdef MPI
+ do i=1,scount(me1)
+#else
+ do i=1,ntot(islice)
+#endif
+c write (9,'(3i5,f10.5)') i,(iparm,potE(i,iparm),iparm=1,nParmSet)
+ do iparm=1,nParmSet
+#ifdef DEBUG
+ write (iout,'(2i5,21f8.2)') i,iparm,
+ & (enetb(k,i,iparm),k=1,26)
+#endif
+ call restore_parm(iparm)
+#ifdef DEBUG
+ write (iout,*) wsc,wscp,welec,wvdwpp,wang,wtor,wscloc,
+ & wcorr,wcorr5,wcorr6,wturn4,wturn3,wturn6,wel_loc,
+ & wtor_d,wsccor,wbond,wsaxs
+#endif
+ do ib=1,nT_h(iparm)
+ if (rescale_mode.eq.1) then
+ quot=1.0d0/(beta_h(ib,iparm)*1.987D-3*T0)
+ quotl=1.0d0
+ kfacl=1.0d0
+ do l=1,5
+ quotl1=quotl
+ quotl=quotl*quot
+ kfacl=kfacl*kfac
+ fT(l)=kfacl/(kfacl-1.0d0+quotl)
+ enddo
+#if defined(FUNCTH)
+ tt = 1.0d0/(beta_h(ib,iparm)*1.987D-3)
+ ft(6)=(320.0d0+80.0d0*dtanh((tt-320.0d0)/80.0d0))/320.0d0
+#elif defined(FUNCT)
+ ft(6)=1.0d0/(beta_h(ib,iparm)*1.987D-3*T0)
+#else
+ ft(6)=1.0d0
+#endif
+ else if (rescale_mode.eq.2) then
+ quot=1.0d0/(T0*beta_h(ib,iparm)*1.987D-3)
+ quotl=1.0d0
+ do l=1,5
+ quotl=quotl*quot
+ fT(l)=1.12692801104297249644d0/
+ & dlog(dexp(quotl)+dexp(-quotl))
+ enddo
+#if defined(FUNCTH)
+ tt = 1.0d0/(beta_h(ib,iparm)*1.987D-3)
+ ft(6)=(320.0d0+80.0d0*dtanh((tt-320.0d0)/80.0d0))/320.0d0
+#elif defined(FUNCT)
+ ft(6)=1.0d0/(beta_h(ib,iparm)*1.987D-3*T0)
+#else
+ ft(6)=1.0d0
+#endif
+c write (iout,*) 1.0d0/(beta_h(ib,iparm)*1.987D-3),ft
+ else if (rescale_mode.eq.0) then
+ do l=1,6
+ fT(l)=1.0d0
+ enddo
+ else
+ write (iout,*) "Error in WHAM_CALC: wrong RESCALE_MODE",
+ & rescale_mode
+ call flush(iout)
+ return1
+ endif
+ evdw=enetb(21,i,iparm)
+ evdw_t=enetb(1,i,iparm)
+#ifdef SCP14
+ evdw2_14=enetb(17,i,iparm)
+ evdw2=enetb(2,i,iparm)+evdw2_14
+#else
+ evdw2=enetb(2,i,iparm)
+ evdw2_14=0.0d0
+#endif
+#ifdef SPLITELE
+ ees=enetb(3,i,iparm)
+ evdw1=enetb(16,i,iparm)
+#else
+ ees=enetb(3,i,iparm)
+ evdw1=0.0d0
+#endif
+ ecorr=enetb(4,i,iparm)
+ ecorr5=enetb(5,i,iparm)
+ ecorr6=enetb(6,i,iparm)
+ eel_loc=enetb(7,i,iparm)
+ eello_turn3=enetb(8,i,iparm)
+ eello_turn4=enetb(9,i,iparm)
+ eturn6=enetb(10,i,iparm)
+ ebe=enetb(11,i,iparm)
+ escloc=enetb(12,i,iparm)
+ etors=enetb(13,i,iparm)
+ etors_d=enetb(14,i,iparm)
+ ehpb=enetb(15,i,iparm)
+ estr=enetb(18,i,iparm)
+ esccor=enetb(19,i,iparm)
+ edihcnstr=enetb(20,i,iparm)
+ ehomology_constr=enetb(22,i,iparm)
+ esaxs=enetb(26,i,iparm)
+#ifdef SPLITELE
+ etot=wsc*(evdw+ft(6)*evdw_t)+wscp*evdw2+ft(1)*welec*ees
+ & +wvdwpp*evdw1
+ & +wang*ebe+ft(1)*wtor*etors+wscloc*escloc
+ & +wstrain*ehpb+nss*ebr+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5
+ & +ft(5)*wcorr6*ecorr6+ft(3)*wturn4*eello_turn4
+ & +ft(2)*wturn3*eello_turn3
+ & +ft(5)*wturn6*eturn6+ft(2)*wel_loc*eel_loc
+ & +edihcnstr+ft(2)*wtor_d*etors_d+ft(1)*wsccor*esccor
+ & +wbond*estr+ehomology_constr+wsaxs*esaxs
+#else
+ etot=wsc*(evdw+ft(6)*evdw_t)+wscp*evdw2
+ & +ft(1)*welec*(ees+evdw1)
+ & +wang*ebe+ft(1)*wtor*etors+wscloc*escloc
+ & +wstrain*ehpb+nss*ebr+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5
+ & +ft(5)*wcorr6*ecorr6+ft(3)*wturn4*eello_turn4
+ & +ft(2)*wturn3*eello_turn3
+ & +ft(5)*wturn6*eturn6+ft(2)*wel_loc*eel_loc+edihcnstr
+ & +ft(2)*wtor_d*etors_d+ft(1)*wsccor*esccor
+ & +wbond*estr+ehomology_constr+wsaxs*esaxs
+
+#endif
+ etot=etot-entfac(i)/beta_h(ib,iparm)
+ if(etot.lt.potEmin_all(ib,iparm)) potEmin_all(ib,iparm)=etot
+
+ enddo ! ib
+ enddo ! iparm
+ enddo ! i
+#ifdef DEBUG
+ write (iout,*) "The potEmin array before reduction"
+ do i=1,nParmSet
+ write (iout,*) "Parameter set",i
+ do j=1,nT_h(i)
+ write (iout,*) j,PotEmin_all(j,i)
+ enddo
+ enddo
+ write (iout,*) "potEmin_min",potEmin_min
+#endif
+#ifdef MPI
+C Determine the minimum energes for all parameter sets and temperatures
+ call MPI_AllReduce(potEmin_all(1,1),potEmin_t_all(1,1),
+ & maxT_h*nParmSet,MPI_DOUBLE_PRECISION,MPI_MIN,WHAM_COMM,IERROR)
+ do i=1,nParmSet
+ do j=1,nT_h(i)
+ potEmin_all(j,i)=potEmin_t_all(j,i)
+ enddo
+ enddo
+#endif
+ potEmin_min=potEmin_all(1,1)
+ do i=1,nParmSet
+ do j=1,nT_h(i)
+ if (potEmin_all(j,i).lt.potEmin_min)
+ & potEmin_min=potEmin_all(j,i)
+ enddo
+ enddo
+#ifdef DEBUG
+ write (iout,*) "The potEmin array"
+ do i=1,nParmSet
+ write (iout,*) "Parameter set",i
+ do j=1,nT_h(i)
+ write (iout,*) j,1.0d0/(1.987d-3*beta_h(j,i)),
+ & PotEmin_all(j,i)
+ enddo
+ enddo
+ write (iout,*) "potEmin_min",potEmin_min
+#endif
+
+
+! Now, put together the histograms from all simulations, in order to get the
+! unbiased total histogram.
#ifdef MPI
do t=0,tmax
hfin_ent_p(t)=0.0d0
estr=enetb(18,t,iparm)
esccor=enetb(19,t,iparm)
edihcnstr=enetb(20,t,iparm)
- edihcnstr=0.0d0
+ ehomology_constr=enetb(22,t,iparm)
+ esaxs=enetb(26,t,iparm)
+ if (homol_nset.gt.1)
+ & ehomology_constr=waga_homology(ihset)*ehomology_constr
do k=0,nGridT
betaT=startGridT+k*delta_T
temper=betaT
c write (iout,*) "ftprim",ftprim
c write (iout,*) "ftbis",ftbis
betaT=1.0d0/(1.987D-3*betaT)
+ if (betaT.ge.beta_h(1,iparm)) then
+ potEmin=potEmin_all(1,iparm)+
+ & (potEmin_all(1,iparm)-potEmin_all(2,iparm))/
+ & (1.0/beta_h(1,iparm)-1.0/beta_h(2,iparm))*
+ & (1.0/betaT-1.0/beta_h(1,iparm))
+#ifdef DEBUG
+ write(iout,*) "first",temper,potEmin
+#endif
+ else if (betaT.le.beta_h(nT_h(iparm),iparm)) then
+ potEmin=potEmin_all(nT_h(iparm),iparm)+
+ &(potEmin_all(nT_h(iparm),iparm)-potEmin_all(nT_h(iparm)-1,iparm))/
+ &(1.0/beta_h(nT_h(iparm),iparm)-1.0/beta_h(nT_h(iparm)-1,iparm))*
+ &(1.0/betaT-1.0/beta_h(nt_h(iparm),iparm))
+#ifdef DEBUG
+ write (iout,*) "last",temper,potEmin
+#endif
+ else
+ do l=1,nT_h(iparm)-1
+ if (betaT.le.beta_h(l,iparm) .and.
+ & betaT.gt.beta_h(l+1,iparm)) then
+ potEmin=potEmin_all(l,iparm)
+#ifdef DEBUG
+ write (iout,*) "l",l,
+ & betaT,1.0d0/(1.987D-3*beta_h(l,iparm)),
+ & 1.0d0/(1.987D-3*beta_h(l+1,iparm)),temper,potEmin
+#endif
+ exit
+ endif
+ enddo
+ endif
#ifdef SPLITELE
etot=wsc*(evdw+ft(6)*evdw_t)+wscp*evdw2+ft(1)*welec*ees
& +wvdwpp*evdw1
& +ft(2)*wturn3*eello_turn3
& +ft(5)*wturn6*eturn6+ft(2)*wel_loc*eel_loc
& +edihcnstr+ft(2)*wtor_d*etors_d+ft(1)*wsccor*esccor
- & +wbond*estr
+ & +wbond*estr+ehomology_constr+wsaxs*esaxs
eprim=ftprim(6)*evdw_t+ftprim(1)*welec*ees
& +ftprim(1)*wtor*etors+
& ftprim(3)*wcorr*ecorr+ftprim(4)*wcorr5*ecorr5+
& +ft(2)*wturn3*eello_turn3
& +ft(5)*wturn6*eturn6+ft(2)*wel_loc*eel_loc+edihcnstr
& +ft(2)*wtor_d*etors_d+ft(1)*wsccor*esccor
- & +wbond*estr
+ & +wbond*estr+ehomology_constr+wsaxs*esaxs
eprim=ftprim(6)*evdw_t+ftprim(1)*welec*(ees+evdw1)
& +ftprim(1)*wtor*etors+
& ftprim(3)*wcorr*ecorr+ftprim(4)*wcorr5*ecorr5+
& ftbis(5)*wcorr6*ecorr6+ftbis(3)*wturn4*eello_turn4+
& ftbis(2)*wturn3*eello_turn3+ftbis(5)*wturn6*eturn6+
& ftbis(2)*wel_loc*eel_loc+ftbis(2)*wtor_d*etors_d+
- & ftprim(1)*wsccor*esccor
+ & ftbis(1)*wsccor*esccor
+
+ endif
+
#endif
weight=dexp(-betaT*(etot-potEmin)+entfac(t))
#ifdef DEBUG
endif
#ifdef MPI
do ib=1,nT_h(iparm)
+ potEmin=potEmin_all(ib,iparm)
expfac=dexp(-beta_h(ib,iparm)*(etot-potEmin)+entfac(t))
hfin_p(ind,ib)=hfin_p(ind,ib)+
& dexp(-beta_h(ib,iparm)*(etot-potEmin)+entfac(t))
enddo
#else
do ib=1,nT_h(iparm)
+ potEmin=potEmin_all(ib,iparm)
expfac=dexp(-beta_h(ib,iparm)*(etot-potEmin)+entfac(t))
hfin(ind,ib)=hfin(ind,ib)+
& dexp(-beta_h(ib,iparm)*(etot-potEmin)+entfac(t))
write (iout,'(a,i3)') "Parameter set",iparm
endif
do i=0,NGridT
+ betaT=1.0d0/(1.987D-3*(startGridT+i*delta_T))
+ if (betaT.ge.beta_h(1,iparm)) then
+ potEmin=potEmin_all(1,iparm)
+ else if (betaT.lt.beta_h(nT_h(iparm),iparm)) then
+ potEmin=potEmin_all(nT_h(iparm),iparm)
+ else
+ do l=1,nT_h(iparm)-1
+ if (betaT.le.beta_h(l,iparm) .and.
+ & betaT.gt.beta_h(l+1,iparm)) then
+ potEmin=potEmin_all(l,iparm)
+ exit
+ endif
+ enddo
+ endif
+
sumE(i,iparm)=sumE(i,iparm)/sumW(i,iparm)
sumEbis(i,iparm)=(startGridT+i*delta_T)*sumEbis(i,iparm)/
& sumW(i,iparm)
write (iout,*)
write (34,'(f7.1,2f15.5,$)') startGridT+i*delta_T,
& sumW(i,iparm),sumE(i,iparm)
- write (34,'(f10.5,$)') (sumQ(j,i,iparm),j=1,nQ+2)
- write (34,'(e15.5,$)') sumEsq(i,iparm)-sumEbis(i,iparm),
- & (sumQsq(j,i,iparm),j=1,nQ+2),(sumEQ(j,i,iparm),j=1,nQ+2)
+ do j=1,nQ+2
+ write (34,'(f10.5,$)') sumQ(j,i,iparm)
+ enddo
+ write (34,'(e15.5,$)') sumEsq(i,iparm)-sumEbis(i,iparm)
+ do j=1,nQ+2
+ write (34,'(e15.5,$)') sumQsq(j,i,iparm)
+ enddo
+ do j=1,nQ+2
+ write (34,'(e15.5,$)') sumEQ(j,i,iparm)
+ enddo
write (34,*)
call flush(34)
enddo