+++ /dev/null
- module wham_calc
-!-----------------------------------------------------------------------------
- use io_units
- use wham_data
-!
- use ene_calc
-#ifdef MPI
- use MPI_data
-! include "COMMON.MPI"
-#endif
- implicit none
-!-----------------------------------------------------------------------------
-!
-!
-!-----------------------------------------------------------------------------
- contains
-!-----------------------------------------------------------------------------
-
- subroutine WHAMCALC(islice,*)
-! Weighed Histogram Analysis Method (WHAM) code
-! Written by A. Liwo based on the work of Kumar et al.,
-! J.Comput.Chem., 13, 1011 (1992)
-!
-! 2/1/05 Multiple temperatures allowed.
-! 2/2/05 Free energies calculated directly from data points
-! acc. to Eq. (21) of Kumar et al.; final histograms also
-! constructed based on this equation.
-! 2/12/05 Multiple parameter sets included
-!
-! 2/2/05 Parallel version
-! use wham_data
-! use io_units
- use names
- use io_base, only:ilen
- use energy_data
-#ifdef MPI
- include "mpif.h"
-#endif
-! include "DIMENSIONS"
-! include "DIMENSIONS.ZSCOPT"
-! include "DIMENSIONS.FREE"
- integer,parameter :: NGridT=400
- integer,parameter :: MaxBinRms=100,MaxBinRgy=100
- integer,parameter :: MaxHdim=200
-! parameter (MaxHdim=200000)
- integer,parameter :: maxinde=200
-#ifdef MPI
- integer :: ierror,errcode,status(MPI_STATUS_SIZE)
-#endif
-! include "COMMON.CONTROL"
-! include "COMMON.IOUNITS"
-! include "COMMON.FREE"
-! include "COMMON.ENERGIES"
-! include "COMMON.FFIELD"
-! include "COMMON.SBRIDGE"
-! include "COMMON.PROT"
-! include "COMMON.ENEPS"
- integer,parameter :: MaxPoint=MaxStr,&
- MaxPointProc=MaxStr_Proc
- real(kind=8),parameter :: finorm_max=1.0d0
- real(kind=8) :: potfac,entmin,entmax,expfac,vf
- integer :: islice
- integer :: i,ii,j,jj,k,kk,l,m,ind,iter,t,tmax,ient,ientmax,iln
- integer :: start,end,iharm,ib,iib,nbin1,nbin,nbin_rms,nbin_rgy,&
- nbin_rmsrgy,liczbaW,iparm,nFi,indrgy,indrms
- integer :: htot(0:MaxHdim),histent(0:2000)
- real(kind=8) :: v(MaxPointProc,MaxR,MaxT_h,nParmSet) !(MaxPointProc,MaxR,MaxT_h,Max_Parm)
- real(kind=8) :: energia(0:n_ene)
-!el real(kind=8) :: energia(0:max_ene)
-#ifdef MPI
- integer :: tmax_t,upindE_p
- real(kind=8) :: fi_p(MaxR,MaxT_h,nParmSet) !(MaxR,MaxT_h,Max_Parm)
- real(kind=8),dimension(0:nGridT,nParmSet) :: sumW_p,sumE_p,&
- sumEbis_p,sumEsq_p !(0:nGridT,Max_Parm)
- real(kind=8),dimension(MaxQ1,0:nGridT,nParmSet) :: sumQ_p,&
- sumQsq_p,sumEQ_p,sumEprim_p !(MaxQ1,0:nGridT,Max_Parm)
- real(kind=8) :: hfin_p(0:MaxHdim,maxT_h),&
- hfin_ent_p(0:MaxHdim),histE_p(0:maxindE),sumH,&
- hrmsrgy_p(0:MaxBinRgy,0:MaxBinRms,maxT_h)
- real(kind=8) :: rgymin_t,rmsmin_t,rgymax_t,rmsmax_t
- real(kind=8) :: potEmin_t,entmin_p,entmax_p
- integer :: histent_p(0:2000)
- logical :: lprint=.true.
-#endif
- real(kind=8) :: delta_T=1.0d0,iientmax
- real(kind=8) :: rgymin,rmsmin,rgymax,rmsmax
- real(kind=8),dimension(0:nGridT,nParmSet) :: sumW,sumE,&
- sumEsq,sumEprim,sumEbis !(0:NGridT,Max_Parm)
- real(kind=8),dimension(MaxQ1,0:nGridT,nParmSet) :: sumQ,&
- sumQsq,sumEQ !(MaxQ1,0:NGridT,Max_Parm)
- real(kind=8) :: betaT,weight,econstr
- real(kind=8) :: fi(MaxR,MaxT_h,nParmSet),& !(MaxR,maxT_h,Max_Parm)
- ddW,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
- real(kind=8) :: 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
- real(kind=8) :: etot,evdw,evdw_t,evdw2,ees,evdw1,ebe,etors,&
- escloc,ehpb,ecorr,ecorr5,ecorr6,eello_turn4,eello_turn3,&
- eello_turn6,eel_loc,edihcnstr,etors_d,estr,evdw2_14,esccor
-
- integer :: ind_point(maxpoint),upindE,indE
- character(len=16) :: plik
- character(len=1) :: licz1
- character(len=2) :: licz2
- character(len=3) :: licz3
- character(len=128) :: nazwa
-! integer ilen
-! external ilen
-!el ientmax=0
-!el ent=0.0d0
- write(licz2,'(bz,i2.2)') islice
- nbin1 = 1.0d0/delta
- write (iout,'(//80(1h-)/"Solving WHAM equations for slice",&
- i2/80(1h-)//)') islice
- write (iout,*) "delta",delta," nbin1",nbin1
- write (iout,*) "MaxN",MaxN," MaxQ",MaxQ," MaHdim",MaxHdim
- call flush(iout)
- dmin=0.0d0
- tmax=0
- potEmin=1.0d10
- rgymin=1.0d10
- rmsmin=1.0d10
- rgymax=0.0d0
- rmsmax=0.0d0
- do t=0,MaxN
- htot(t)=0
- enddo
-#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
- 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)
- if (q(nQ+2,i).gt.rgymax) rgymax=q(nQ+2,i)
- ind_point(i)=0
- do j=nQ,1,-1
- ind=(q(j,i)-dmin+1.0d-8)/delta
- if (j.eq.1) then
- ind_point(i)=ind_point(i)+ind
- else
- ind_point(i)=ind_point(i)+nbin1**(j-1)*ind
- endif
-! write (iout,*) "i",i," j",j," q",q(j,i)," ind_point",
-! & ind_point(i)
- call flush(iout)
- if (ind_point(i).lt.0 .or. ind_point(i).gt.MaxHdim) then
- write (iout,*) "Error - index exceeds range for point",i,&
- " q=",q(j,i)," ind",ind_point(i)
-#ifdef MPI
- write (iout,*) "Processor",me1
- call flush(iout)
- call MPI_Abort(MPI_COMM_WORLD, Ierror, Errcode )
-#endif
- stop
- endif
- enddo ! j
- if (ind_point(i).gt.tmax) tmax=ind_point(i)
- htot(ind_point(i))=htot(ind_point(i))+1
-#ifdef DEBUG
- write (iout,*) "i",i,"q",(q(j,i),j=1,nQ)," ind",ind_point(i),&
- " htot",htot(ind_point(i))
- call flush(iout)
-#endif
- enddo ! i
- call flush(iout)
-
- nbin=nbin1**nQ-1
- write (iout,'(a)') "Numbers of counts in Q bins"
- do t=0,tmax
- if (htot(t).gt.0) then
- write (iout,'(i15,$)') t
- liczbaW=t
- do j=1,nQ
- jj = mod(liczbaW,nbin1)
- liczbaW=liczbaW/nbin1
- write (iout,'(i5,$)') jj
- enddo
- write (iout,'(i8)') htot(t)
- endif
- enddo
- do iparm=1,nParmSet
- write (iout,'(a,i3)') "Number of data points for parameter set",&
- iparm
- write (iout,'(i7,$)') ((snk(m,ib,iparm,islice),m=1,nr(ib,iparm)),&
- ib=1,nT_h(iparm))
- write (iout,'(i8)') stot(islice)
- write (iout,'(a)')
- enddo
- call flush(iout)
-
-#ifdef MPI
- 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)
- 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_MAX,WHAM_COMM,IERROR)
- call MPI_AllReduce(rgymin,rgymin_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
- rgymin=rgymin_t
- rgymax=rgymax_t
- rmsmin=rmsmin_t
- rmsmax=rmsmax_t
- write (iout,*) "potEmin",potEmin
-#endif
- rmsmin=deltrms*dint(rmsmin/deltrms)
- rmsmax=deltrms*dint(rmsmax/deltrms)
- rgymin=deltrms*dint(rgymin/deltrgy)
- rgymax=deltrms*dint(rgymax/deltrgy)
- nbin_rms=(rmsmax-rmsmin)/deltrms
- nbin_rgy=(rgymax-rgymin)/deltrgy
- write (iout,*) "rmsmin",rmsmin," rmsmax",rmsmax," rgymin",rgymin,&
- " rgymax",rgymax," nbin_rms",nbin_rms," nbin_rgy",nbin_rgy
- nFi=0
- do i=1,nParmSet
- do j=1,nT_h(i)
- nFi=nFi+nR(j,i)
- enddo
- enddo
- write (iout,*) "nFi",nFi
-! Compute the Boltzmann factor corresponing to restrain potentials in different
-! simulations.
-#ifdef MPI
- do i=1,scount(me1)
-#else
- do i=1,ntot(islice)
-#endif
-! 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,21)
-!#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
-!#endif
- do ib=1,nT_h(iparm)
-!el old rascale weights
-!
-! if (rescale_modeW.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_modeW.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
-! write (iout,*) 1.0d0/(beta_h(ib,iparm)*1.987D-3),ft
-! else if (rescale_modeW.eq.0) then
-! do l=1,6
-! fT(l)=1.0d0
-! enddo
-! else
-! write (iout,*) "Error in WHAM_CALC: wrong RESCALE_MODE",&
-! rescale_modeW
-! call flush(iout)
-! return 1
-! endif
-! el end old rescale weights
- call rescale_weights(1.0d0/(beta_h(ib,iparm)*1.987D-3))
-
-! call etot(enetb(0,i,iparm))
- evdw=enetb(1,i,iparm)
-! evdw_t=enetb(21,i,iparm)
- evdw_t=enetb(20,i,iparm)
-#ifdef SCP14
-! evdw2_14=enetb(17,i,iparm)
- evdw2_14=enetb(18,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)
- eello_turn6=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)
- estr=enetb(17,i,iparm)
-! esccor=enetb(19,i,iparm)
- esccor=enetb(21,i,iparm)
-! edihcnstr=enetb(20,i,iparm)
- edihcnstr=enetb(19,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
-#endif
-
-!#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
-!#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
-!#endif
-
-#ifdef SPLITELE
- etot=wsc*evdw+wscp*evdw2+welec*ees &
- +wvdwpp*evdw1 &
- +wang*ebe+wtor*etors+wscloc*escloc &
- +wstrain*ehpb+nss*ebr+wcorr*ecorr+wcorr5*ecorr5 &
- +wcorr6*ecorr6+wturn4*eello_turn4 &
- +wturn3*eello_turn3 &
- +wturn6*eello_turn6+wel_loc*eel_loc &
- +edihcnstr+wtor_d*etors_d+wsccor*esccor &
- +wbond*estr
-#else
- etot=wsc*evdw+wscp*evdw2 &
- +welec*(ees+evdw1) &
- +wang*ebe+wtor*etors+wscloc*escloc &
- +wstrain*ehpb+nss*ebr+wcorr*ecorr+wcorr5*ecorr5 &
- +wcorr6*ecorr6+wturn4*eello_turn4 &
- +wturn3*eello_turn3 &
- +wturn6*eello_turn6+wel_loc*eel_loc+edihcnstr &
- +wtor_d*etors_d+wsccor*esccor &
- +wbond*estr
-#endif
-
-#ifdef DEBUG
- write (iout,*) i,iparm,1.0d0/(beta_h(ib,iparm)*1.987D-3),&
- etot,potEmin
-#endif
-#ifdef DEBUG
- if (iparm.eq.1 .and. ib.eq.1) then
- write (iout,*)"Conformation",i
- energia(0)=etot
- do k=1,max_ene
- energia(k)=enetb(k,i,iparm)
- enddo
-! call enerprint(energia(0),fT)
- call enerprint(energia(0))
- endif
-#endif
- do kk=1,nR(ib,iparm)
- Econstr=0.0d0
- do j=1,nQ
- ddW = q(j,i)
- Econstr=Econstr+Kh(j,kk,ib,iparm) &
- *(ddW-q0(j,kk,ib,iparm))**2
- enddo
- v(i,kk,ib,iparm)= &
- -beta_h(ib,iparm)*(etot-potEmin+Econstr)
-#ifdef DEBUG
- write (iout,'(4i5,4e15.5)') i,kk,ib,iparm,&
- etot,potEmin,etot-potEmin,v(i,kk,ib,iparm)
-#endif
- enddo ! kk
- enddo ! ib
- enddo ! iparm
- enddo ! i
-! Simple iteration to calculate free energies corresponding to all simulation
-! runs.
- do iter=1,maxit
-
-! Compute new free-energy values corresponding to the righ-hand side of the
-! equation and their derivatives.
- write (iout,*) "------------------------fi"
-#ifdef MPI
- do t=1,scount(me1)
-#else
- do t=1,ntot(islice)
-#endif
- vmax=-1.0d+20
- do i=1,nParmSet
- do k=1,nT_h(i)
- do l=1,nR(k,i)
- vf=v(t,l,k,i)+f(l,k,i)
- if (vf.gt.vmax) vmax=vf
- enddo
- enddo
- enddo
- denom=0.0d0
- do i=1,nParmSet
- do k=1,nT_h(i)
- do l=1,nR(k,i)
- aux=f(l,k,i)+v(t,l,k,i)-vmax
- if (aux.gt.-200.0d0) &
- denom=denom+snk(l,k,i,islice)*dexp(aux)
- enddo
- enddo
- enddo
- entfac(t)=-dlog(denom)-vmax
-#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(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))
-#ifdef DEBUG
- write (iout,'(4i5,3e15.5)') t,ii,iib,iparm,&
- v(t,ii,iib,iparm),entfac(t),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))
- enddo
-#endif
- enddo ! ii
- enddo ! iib
- enddo ! iparm
-
-#ifdef MPI
-#ifdef DEBUG
- write (iout,*) "fi before MPI_Reduce me",me,' master',master
- do iparm=1,nParmSet
- do ib=1,nT_h(nparmset)
- write (iout,*) "iparm",iparm," ib",ib
- write (iout,*) "beta=",beta_h(ib,iparm)
- write (iout,'(8e15.5)') (fi_p(i,ib,iparm),i=1,nR(ib,iparm))
- enddo
- enddo
-#endif
- write (iout,*) "REDUCE size",maxR,MaxT_h,nParmSet,&
- maxR*MaxT_h*nParmSet
- write (iout,*) "MPI_COMM_WORLD",MPI_COMM_WORLD,&
- " WHAM_COMM",WHAM_COMM
- 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)
-#ifdef DEBUG
- write (iout,*) "fi after MPI_Reduce nparmset",nparmset
- do iparm=1,nParmSet
- write (iout,*) "iparm",iparm
- do ib=1,nT_h(iparm)
- write (iout,*) "beta=",beta_h(ib,iparm)
- write (iout,'(8e15.5)') (fi(i,ib,iparm),i=1,nR(ib,iparm))
- enddo
- enddo
-#endif
- if (me1.eq.Master) then
-#endif
- avefi=0.0d0
- 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))
- avefi=avefi+fi(i,ib,iparm)
- enddo
- enddo
- enddo
- avefi=avefi/nFi
- do iparm=1,nParmSet
- write (iout,*) "Parameter set",iparm
- do ib =1,nT_h(iparm)
- write (iout,*) "beta=",beta_h(ib,iparm)
- do i=1,nR(ib,iparm)
- fi(i,ib,iparm)=fi(i,ib,iparm)-avefi
- enddo
- write (iout,'(8f10.5)') (fi(i,ib,iparm),i=1,nR(ib,iparm))
- write (iout,'(8f10.5)') (f(i,ib,iparm),i=1,nR(ib,iparm))
- enddo
- enddo
-
-! Compute the norm of free-energy increments.
- finorm=0.0d0
- do iparm=1,nParmSet
- do ib=1,nT_h(iparm)
- do i=1,nR(ib,iparm)
- finorm=finorm+dabs(fi(i,ib,iparm)-f(i,ib,iparm))
- f(i,ib,iparm)=fi(i,ib,iparm)
- enddo
- enddo
- enddo
-
- write (iout,*) 'Iteration',iter,' finorm',finorm
-
-#ifdef MPI
- endif
- call MPI_Bcast(f(1,1,1),MaxR*MaxT_h*nParmSet,&
- MPI_DOUBLE_PRECISION,Master,&
- WHAM_COMM,IERROR)
- call MPI_Bcast(finorm,1,MPI_DOUBLE_PRECISION,Master,&
- WHAM_COMM,IERROR)
-#endif
-! Exit, if the increment norm is smaller than pre-assigned tolerance.
- if (finorm.lt.fimin) then
- write (iout,*) 'Iteration converged'
- goto 20
- endif
-
- enddo ! iter
-
- 20 continue
-! 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
- enddo
-#else
- do t=0,tmax
- hfin_ent(t)=0.0d0
- enddo
-#endif
- write (iout,*) "--------------hist"
-#ifdef MPI
- do iparm=1,nParmSet
- do i=0,nGridT
- sumW_p(i,iparm)=0.0d0
- sumE_p(i,iparm)=0.0d0
- sumEbis_p(i,iparm)=0.0d0
- sumEsq_p(i,iparm)=0.0d0
- do j=1,nQ+2
- sumQ_p(j,i,iparm)=0.0d0
- sumQsq_p(j,i,iparm)=0.0d0
- sumEQ_p(j,i,iparm)=0.0d0
- enddo
- enddo
- enddo
- upindE_p=0
-#else
- do iparm=1,nParmSet
- do i=0,nGridT
- sumW(i,iparm)=0.0d0
- sumE(i,iparm)=0.0d0
- sumEbis(i,iparm)=0.0d0
- sumEsq(i,iparm)=0.0d0
- do j=1,nQ+2
- sumQ(j,i,iparm)=0.0d0
- sumQsq(j,i,iparm)=0.0d0
- sumEQ(j,i,iparm)=0.0d0
- enddo
- enddo
- enddo
- upindE=0
-#endif
-! 8/26/05 entropy distribution
-#ifdef MPI
- entmin_p=1.0d10
- entmax_p=-1.0d10
- do t=1,scount(me1)
-! ent=-dlog(entfac(t))
- ent=entfac(t)
- if (ent.lt.entmin_p) entmin_p=ent
- if (ent.gt.entmax_p) entmax_p=ent
- enddo
- write (iout,*) "entmin",entmin_p," entmax",entmax_p
-! write (iout,*) "entmin_p",entmin_p," entmax_p",entmax_p
- call flush(iout)
- call MPI_Allreduce(entmin_p,entmin,1,MPI_DOUBLE_PRECISION,MPI_MIN,&
- WHAM_COMM,IERROR)
- call MPI_Allreduce(entmax_p,entmax,1,MPI_DOUBLE_PRECISION,MPI_MAX,&
- WHAM_COMM,IERROR)
- write (iout,*) "entmin",entmin_p," entmax",entmax_p
-! write (iout,*) "entmin_p",entmin_p," entmax_p",entmax_p
- ientmax=entmax-entmin
-!iientmax=entmax-entmin !el
-!write (iout,*) "ientmax",ientmax,entmax,entmin
-!write (iout,*) "iientmax",iientmax
- if (ientmax.gt.2000) ientmax=2000
- write (iout,*) "entmin",entmin," entmax",entmax," ientmax",ientmax
- call flush(iout)
- do t=1,scount(me1)
-! ient=-dlog(entfac(t))-entmin
- ient=entfac(t)-entmin
- if (ient.le.2000) histent_p(ient)=histent_p(ient)+1
- enddo
- call MPI_Allreduce(histent_p(0),histent(0),ientmax+1,MPI_INTEGER,&
- MPI_SUM,WHAM_COMM,IERROR)
- if (me1.eq.Master) then
- write (iout,*) "Entropy histogram"
- do i=0,ientmax
- write(iout,'(f15.4,i10)') entmin+i,histent(i)
- enddo
- endif
-#else
- entmin=1.0d10
- entmax=-1.0d10
- do t=1,ntot(islice)
- ent=entfac(t)
- if (ent.lt.entmin) entmin=ent
- if (ent.gt.entmax) entmax=ent
- enddo
- ientmax=-dlog(entmax)-entmin
- if (ientmax.gt.2000) ientmax=2000
- do t=1,ntot(islice)
- ient=entfac(t)-entmin
- if (ient.le.2000) histent(ient)=histent(ient)+1
- enddo
- write (iout,*) "Entropy histogram"
- do i=0,ientmax
- write(iout,'(2f15.4)') entmin+i,histent(i)
- enddo
-#endif
-
-#ifdef MPI
- write (iout,*) "me1",me1," scount",scount(me1) !d
-
- do iparm=1,nParmSet
-
-#ifdef MPI
- do ib=1,nT_h(iparm)
- do t=0,tmax
- hfin_p(t,ib)=0.0d0
- enddo
- enddo
- do i=1,maxindE
- histE_p(i)=0.0d0
- enddo
-#else
- do ib=1,nT_h(iparm)
- do t=0,tmax
- hfin(t,ib)=0.0d0
- enddo
- enddo
- do i=1,maxindE
- histE(i)=0.0d0
- enddo
-#endif
- do ib=1,nT_h(iparm)
- do i=0,MaxBinRms
- do j=0,MaxBinRgy
- hrmsrgy(j,i,ib)=0.0d0
-#ifdef MPI
- hrmsrgy_p(j,i,ib)=0.0d0
-#endif
- enddo
- enddo
- enddo
-
- do t=1,scount(me1)
-#else
- do t=1,ntot(islice)
-#endif
- ind = ind_point(t)
-#ifdef MPI
- hfin_ent_p(ind)=hfin_ent_p(ind)+dexp(entfac(t))
-#else
- hfin_ent(ind)=hfin_ent(ind)+dexp(entfac(t))
-#endif
-! write (iout,'(2i5,20f8.2)') "debug",t,t,(enetb(k,t,iparm),k=1,21)
- call restore_parm(iparm)
-! evdw=enetb(21,t,iparm)
- evdw=enetb(20,t,iparm)
- evdw_t=enetb(1,t,iparm)
-#ifdef SCP14
-! evdw2_14=enetb(17,t,iparm)
- evdw2_14=enetb(18,t,iparm)
- evdw2=enetb(2,t,iparm)+evdw2_14
-#else
- evdw2=enetb(2,t,iparm)
- evdw2_14=0.0d0
-#endif
-#ifdef SPLITELE
- ees=enetb(3,t,iparm)
- evdw1=enetb(16,t,iparm)
-#else
- ees=enetb(3,t,iparm)
- evdw1=0.0d0
-#endif
- ecorr=enetb(4,t,iparm)
- ecorr5=enetb(5,t,iparm)
- ecorr6=enetb(6,t,iparm)
- eel_loc=enetb(7,t,iparm)
- eello_turn3=enetb(8,t,iparm)
- eello_turn4=enetb(9,t,iparm)
- eello_turn6=enetb(10,t,iparm)
- ebe=enetb(11,t,iparm)
- escloc=enetb(12,t,iparm)
- etors=enetb(13,t,iparm)
- etors_d=enetb(14,t,iparm)
- ehpb=enetb(15,t,iparm)
-! estr=enetb(18,t,iparm)
- estr=enetb(17,t,iparm)
-! esccor=enetb(19,t,iparm)
- esccor=enetb(21,t,iparm)
-! edihcnstr=enetb(20,t,iparm)
- edihcnstr=enetb(19,t,iparm)
- edihcnstr=0.0d0
- do k=0,nGridT
- betaT=startGridT+k*delta_T
- temper=betaT
-!write(iout,*)"kkkkkkkk",betaT,startGridT,k,delta_T
-!d fT=T0/betaT
-!d ft=2*T0/(T0+betaT)
- if (rescale_modeW.eq.1) then
- quot=betaT/T0
- quotl=1.0d0
- kfacl=1.0d0
- do l=1,5
- quotl1=quotl
- quotl=quotl*quot
- kfacl=kfacl*kfac
- denom=kfacl-1.0d0+quotl
- fT(l)=kfacl/denom
- ftprim(l)=-l*ft(l)*quotl1/(T0*denom)
- ftbis(l)=l*kfacl*quotl1* &
- (2*l*quotl-(l-1)*denom)/(quot*t0*t0*denom**3)
- enddo
-#if defined(FUNCTH)
- ft(6)=(320.0d0+80.0d0*dtanh((betaT-320.0d0)/80.0d0))/ &
- 320.0d0
- ftprim(6)=1.0d0/(320.0d0*dcosh((betaT-320.0d0)/80.0d0)**2)
- ftbis(6)=-2.0d0*dtanh((betaT-320.0d0)/80.0d0) &
- /(320.0d0*80.0d0*dcosh((betaT-320.0d0)/80.0d0)**3)
-#elif defined(FUNCT)
- fT(6)=betaT/T0
- ftprim(6)=1.0d0/T0
- ftbis(6)=0.0d0
-#else
- fT(6)=1.0d0
- ftprim(6)=0.0d0
- ftbis(6)=0.0d0
-#endif
- else if (rescale_modeW.eq.2) then
- quot=betaT/T0
- quotl=1.0d0
- do l=1,5
- quotl1=quotl
- quotl=quotl*quot
- eplus=dexp(quotl)
- eminus=dexp(-quotl)
- logfac=1.0d0/dlog(eplus+eminus)
- tanhT=(eplus-eminus)/(eplus+eminus)
- fT(l)=1.12692801104297249644d0*logfac
- ftprim(l)=-l*quotl1*ft(l)*tanhT*logfac/T0
- ftbis(l)=(l-1)*ftprim(l)/(quot*T0)- &
- 2*l*quotl1/T0*logfac* &
- (2*l*quotl1*ft(l)/(T0*(eplus+eminus)**2) &
- +ftprim(l)*tanhT)
- enddo
-#if defined(FUNCTH)
- ft(6)=(320.0d0+80.0d0*dtanh((betaT-320.0d0)/80.0d0))/ &
- 320.0d0
- ftprim(6)=1.0d0/(320.0d0*dcosh((betaT-320.0d0)/80.0d0)**2)
- ftbis(6)=-2.0d0*dtanh((betaT-320.0d0)/80.0d0) &
- /(320.0d0*80.0d0*dcosh((betaT-320.0d0)/80.0d0)**3)
-#elif defined(FUNCT)
- fT(6)=betaT/T0
- ftprim(6)=1.0d0/T0
- ftbis(6)=0.0d0
-#else
- fT(6)=1.0d0
- ftprim(6)=0.0d0
- ftbis(6)=0.0d0
-#endif
- else if (rescale_modeW.eq.0) then
- do l=1,5
- fT(l)=1.0d0
- ftprim(l)=0.0d0
- enddo
- else
- write (iout,*) "Error in WHAM_CALC: wrong RESCALE_MODE",&
- rescale_modeW
- call flush(iout)
- return 1
- endif
-! write (iout,*) "ftprim",ftprim
-! write (iout,*) "ftbis",ftbis
- betaT=1.0d0/(1.987D-3*betaT)
-#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*eello_turn6+ft(2)*wel_loc*eel_loc &
- +edihcnstr+ft(2)*wtor_d*etors_d+ft(1)*wsccor*esccor &
- +wbond*estr
- eprim=ftprim(6)*evdw_t+ftprim(1)*welec*ees &
- +ftprim(1)*wtor*etors+ &
- ftprim(3)*wcorr*ecorr+ftprim(4)*wcorr5*ecorr5+ &
- ftprim(5)*wcorr6*ecorr6+ftprim(3)*wturn4*eello_turn4+ &
- ftprim(2)*wturn3*eello_turn3+ftprim(5)*wturn6*eello_turn6+ &
- ftprim(2)*wel_loc*eel_loc+ftprim(2)*wtor_d*etors_d+ &
- ftprim(1)*wsccor*esccor
- ebis=ftbis(1)*welec*ees+ftbis(1)*wtor*etors+ &
- ftbis(3)*wcorr*ecorr+ftbis(4)*wcorr5*ecorr5+ &
- ftbis(5)*wcorr6*ecorr6+ftbis(3)*wturn4*eello_turn4+ &
- ftbis(2)*wturn3*eello_turn3+ftbis(5)*wturn6*eello_turn6+ &
- ftbis(2)*wel_loc*eel_loc+ftbis(2)*wtor_d*etors_d+ &
- ftbis(1)*wsccor*esccor
-#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*eello_turn6+ft(2)*wel_loc*eel_loc+edihcnstr &
- +ft(2)*wtor_d*etors_d+ft(1)*wsccor*esccor &
- +wbond*estr
- eprim=ftprim(6)*evdw_t+ftprim(1)*welec*(ees+evdw1) &
- +ftprim(1)*wtor*etors+ &
- ftprim(3)*wcorr*ecorr+ftprim(4)*wcorr5*ecorr5+ &
- ftprim(5)*wcorr6*ecorr6+ftprim(3)*wturn4*eello_turn4+ &
- ftprim(2)*wturn3*eello_turn3+ftprim(5)*wturn6*eello_turn6+ &
- ftprim(2)*wel_loc*eel_loc+ftprim(2)*wtor_d*etors_d+ &
- ftprim(1)*wsccor*esccor
- ebis=ftbis(1)*welec*(ees+evdw1)+ftbis(1)*wtor*etors+ &
- ftbis(3)*wcorr*ecorr+ftbis(4)*wcorr5*ecorr5+ &
- ftbis(5)*wcorr6*ecorr6+ftbis(3)*wturn4*eello_turn4+ &
- ftbis(2)*wturn3*eello_turn3+ftbis(5)*wturn6*eello_turn6+ &
- ftbis(2)*wel_loc*eel_loc+ftbis(2)*wtor_d*etors_d+ &
- ftprim(1)*wsccor*esccor
-#endif
- weight=dexp(-betaT*(etot-potEmin)+entfac(t))
-#ifdef DEBUG
- write (iout,*) "iparm",iparm," t",t," betaT",betaT,&
- " etot",etot," entfac",entfac(t),&
- " weight",weight," ebis",ebis
-#endif
- etot=etot-temper*eprim
-#ifdef MPI
- sumW_p(k,iparm)=sumW_p(k,iparm)+weight
- sumE_p(k,iparm)=sumE_p(k,iparm)+etot*weight
- sumEbis_p(k,iparm)=sumEbis_p(k,iparm)+ebis*weight
- sumEsq_p(k,iparm)=sumEsq_p(k,iparm)+etot**2*weight
- do j=1,nQ+2
- sumQ_p(j,k,iparm)=sumQ_p(j,k,iparm)+q(j,t)*weight
- sumQsq_p(j,k,iparm)=sumQsq_p(j,k,iparm)+q(j,t)**2*weight
- sumEQ_p(j,k,iparm)=sumEQ_p(j,k,iparm) &
- +etot*q(j,t)*weight
- enddo
-#else
- sumW(k,iparm)=sumW(k,iparm)+weight
- sumE(k,iparm)=sumE(k,iparm)+etot*weight
- sumEbis(k,iparm)=sumEbis(k,iparm)+ebis*weight
- sumEsq(k,iparm)=sumEsq(k,iparm)+etot**2*weight
- do j=1,nQ+2
- sumQ(j,k,iparm)=sumQ(j,k,iparm)+q(j,t)*weight
- sumQsq(j,k,iparm)=sumQsq(j,k,iparm)+q(j,t)**2*weight
- sumEQ(j,k,iparm)=sumEQ(j,k,iparm) &
- +etot*q(j,t)*weight
- enddo
-#endif
- enddo
- indE = aint(potE(t,iparm)-aint(potEmin))
- if (indE.ge.0 .and. indE.le.maxinde) then
- if (indE.gt.upindE_p) upindE_p=indE
- histE_p(indE)=histE_p(indE)+dexp(-entfac(t))
- endif
-#ifdef MPI
- do ib=1,nT_h(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))
- if (rmsrgymap) then
- indrgy=dint((q(nQ+2,t)-rgymin)/deltrgy)
- indrms=dint((q(nQ+1,t)-rmsmin)/deltrms)
- hrmsrgy_p(indrgy,indrms,ib)= &
- hrmsrgy_p(indrgy,indrms,ib)+expfac
- endif
- enddo
-#else
- do ib=1,nT_h(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))
- if (rmsrgymap) then
- indrgy=dint((q(nQ+2,t)-rgymin)/deltrgy)
- indrms=dint((q(nQ+1,t)-rmsmin)/deltrms)
- hrmsrgy(indrgy,indrms,ib)= &
- hrmsrgy(indrgy,indrms,ib)+expfac
- endif
- enddo
-#endif
- enddo ! t
- do ib=1,nT_h(iparm)
- if (histout) call MPI_Reduce(hfin_p(0,ib),hfin(0,ib),nbin,&
- MPI_DOUBLE_PRECISION,MPI_SUM,Master,WHAM_COMM,IERROR)
- if (rmsrgymap) then
- call MPI_Reduce(hrmsrgy_p(0,0,ib),hrmsrgy(0,0,ib),&
- (MaxBinRgy+1)*(nbin_rms+1),MPI_DOUBLE_PRECISION,MPI_SUM,Master,&
- WHAM_COMM,IERROR)
- endif
- enddo
- call MPI_Reduce(upindE_p,upindE,1,&
- MPI_INTEGER,MPI_MAX,Master,WHAM_COMM,IERROR)
- call MPI_Reduce(histE_p(0),histE(0),maxindE,&
- MPI_DOUBLE_PRECISION,MPI_SUM,Master,WHAM_COMM,IERROR)
-
- if (me1.eq.master) then
-
- if (histout) then
-
- write (iout,'(6x,$)')
- write (iout,'(f20.2,$)') (1.0d0/(1.987D-3*beta_h(ib,iparm)),&
- ib=1,nT_h(iparm))
- write (iout,*)
-
- write (iout,'(/a)') 'Final histograms'
- if (histfile) then
- if (nslice.eq.1) then
- if (separate_parset) then
- write(licz3,"(bz,i3.3)") myparm
- histname=prefix(:ilen(prefix))//'_par'//licz3//'.hist'
- else
- histname=prefix(:ilen(prefix))//'.hist'
- endif
- else
- if (separate_parset) then
- write(licz3,"(bz,i3.3)") myparm
- histname=prefix(:ilen(prefix))//'_par'//licz3// &
- '_slice_'//licz2//'.hist'
- else
- histname=prefix(:ilen(prefix))//'_slice_'//licz2//'.hist'
- endif
- endif
-#if defined(AIX) || defined(PGI)
- open (ihist,file=histname,position='append')
-#else
- open (ihist,file=histname,access='append')
-#endif
- endif
-
- do t=0,tmax
- liczbaW=t
- sumH=0.0d0
- do ib=1,nT_h(iparm)
- sumH=sumH+hfin(t,ib)
- enddo
- if (sumH.gt.0.0d0) then
- do j=1,nQ
- jj = mod(liczbaW,nbin1)
- liczbaW=liczbaW/nbin1
- write (iout,'(f6.3,$)') dmin+(jj+0.5d0)*delta
- if (histfile) &
- write (ihist,'(f6.3,$)') dmin+(jj+0.5d0)*delta
- enddo
- do ib=1,nT_h(iparm)
- write (iout,'(e20.10,$)') hfin(t,ib)
- if (histfile) write (ihist,'(e20.10,$)') hfin(t,ib)
- enddo
- write (iout,'(i5)') iparm
- if (histfile) write (ihist,'(i5)') iparm
- endif
- enddo
-
- endif
-
- if (entfile) then
- if (nslice.eq.1) then
- if (separate_parset) then
- write(licz3,"(bz,i3.3)") myparm
- histname=prefix(:ilen(prefix))//"_par"//licz3//'.ent'
- else
- histname=prefix(:ilen(prefix))//'.ent'
- endif
- else
- if (separate_parset) then
- write(licz3,"(bz,i3.3)") myparm
- histname=prefix(:ilen(prefix))//'par_'//licz3// &
- '_slice_'//licz2//'.ent'
- else
- histname=prefix(:ilen(prefix))//'_slice_'//licz2//'.ent'
- endif
- endif
-#if defined(AIX) || defined(PGI)
- open (ihist,file=histname,position='append')
-#else
- open (ihist,file=histname,access='append')
-#endif
- write (ihist,'(a)') "# Microcanonical entropy"
- do i=0,upindE
- write (ihist,'(f8.0,$)') dint(potEmin)+i
- if (histE(i).gt.0.0e0) then
- write (ihist,'(f15.5,$)') dlog(histE(i))
- else
- write (ihist,'(f15.5,$)') 0.0d0
- endif
- enddo
- write (ihist,*)
- close(ihist)
- endif
- write (iout,*) "Microcanonical entropy"
- do i=0,upindE
- write (iout,'(f8.0,$)') dint(potEmin)+i
- if (histE(i).gt.0.0e0) then
- write (iout,'(f15.5,$)') dlog(histE(i))
- else
- write (iout,'(f15.5,$)') 0.0d0
- endif
- write (iout,*)
- enddo
- if (rmsrgymap) then
- if (nslice.eq.1) then
- if (separate_parset) then
- write(licz3,"(bz,i3.3)") myparm
- histname=prefix(:ilen(prefix))//'_par'//licz3//'.rmsrgy'
- else
- histname=prefix(:ilen(prefix))//'.rmsrgy'
- endif
- else
- if (separate_parset) then
- write(licz3,"(bz,i3.3)") myparm
- histname=prefix(:ilen(prefix))//'_par'//licz3// &
- '_slice_'//licz2//'.rmsrgy'
- else
- histname=prefix(:ilen(prefix))//'_slice_'//licz2//'.rmsrgy'
- endif
- endif
-#if defined(AIX) || defined(PGI)
- open (ihist,file=histname,position='append')
-#else
- open (ihist,file=histname,access='append')
-#endif
- do i=0,nbin_rms
- do j=0,nbin_rgy
- write(ihist,'(2f8.2,$)') &
- rgymin+deltrgy*j,rmsmin+deltrms*i
- do ib=1,nT_h(iparm)
- if (hrmsrgy(j,i,ib).gt.0.0d0) then
- write(ihist,'(e14.5,$)') &
- -dlog(hrmsrgy(j,i,ib))/beta_h(ib,iparm) &
- +potEmin
- else
- write(ihist,'(e14.5,$)') 1.0d6
- endif
- enddo
- write (ihist,'(i2)') iparm
- enddo
- enddo
- close(ihist)
- endif
- endif
- enddo ! iparm
-#ifdef MPI
- call MPI_Reduce(hfin_ent_p(0),hfin_ent(0),nbin,&
- MPI_DOUBLE_PRECISION,MPI_SUM,Master,WHAM_COMM,IERROR)
- call MPI_Reduce(sumW_p(0,1),sumW(0,1),(nGridT+1)*nParmSet,&
- MPI_DOUBLE_PRECISION,MPI_SUM,Master,WHAM_COMM,IERROR)
- call MPI_Reduce(sumE_p(0,1),sumE(0,1),(nGridT+1)*nParmSet,&
- MPI_DOUBLE_PRECISION,MPI_SUM,Master,WHAM_COMM,IERROR)
- call MPI_Reduce(sumEbis_p(0,1),sumEbis(0,1),(nGridT+1)*nParmSet,&
- MPI_DOUBLE_PRECISION,MPI_SUM,Master,WHAM_COMM,IERROR)
- call MPI_Reduce(sumEsq_p(0,1),sumEsq(0,1),(nGridT+1)*nParmSet,&
- MPI_DOUBLE_PRECISION,MPI_SUM,Master,WHAM_COMM,IERROR)
- call MPI_Reduce(sumQ_p(1,0,1),sumQ(1,0,1),&
- MaxQ1*(nGridT+1)*nParmSet,MPI_DOUBLE_PRECISION,MPI_SUM,Master,&
- WHAM_COMM,IERROR)
- call MPI_Reduce(sumQsq_p(1,0,1),sumQsq(1,0,1),&
- MaxQ1*(nGridT+1)*nParmSet,MPI_DOUBLE_PRECISION,MPI_SUM,Master,&
- WHAM_COMM,IERROR)
- call MPI_Reduce(sumEQ_p(1,0,1),sumEQ(1,0,1),&
- MaxQ1*(nGridT+1)*nParmSet,MPI_DOUBLE_PRECISION,MPI_SUM,Master,&
- WHAM_COMM,IERROR)
- if (me.eq.master) then
-#endif
- write (iout,'(/a)') 'Thermal characteristics of folding'
- if (nslice.eq.1) then
- nazwa=prefix
- else
- nazwa=prefix(:ilen(prefix))//"_slice_"//licz2
- endif
- iln=ilen(nazwa)
- if (nparmset.eq.1 .and. .not.separate_parset) then
- nazwa=nazwa(:iln)//".thermal"
- else if (nparmset.eq.1 .and. separate_parset) then
- write(licz3,"(bz,i3.3)") myparm
- nazwa=nazwa(:iln)//"_par_"//licz3//".thermal"
- endif
- do iparm=1,nParmSet
- if (nparmset.gt.1) then
- write(licz3,"(bz,i3.3)") iparm
- nazwa=nazwa(:iln)//"_par_"//licz3//".thermal"
- endif
- open(34,file=nazwa)
- if (separate_parset) then
- write (iout,'(a,i3)') "Parameter set",myparm
- else
- write (iout,'(a,i3)') "Parameter set",iparm
- endif
- do i=0,NGridT
- sumE(i,iparm)=sumE(i,iparm)/sumW(i,iparm)
- sumEbis(i,iparm)=(startGridT+i*delta_T)*sumEbis(i,iparm)/ &
- sumW(i,iparm)
- sumEsq(i,iparm)=(sumEsq(i,iparm)/sumW(i,iparm) &
- -sumE(i,iparm)**2)/(1.987D-3*(startGridT+i*delta_T)**2)
- do j=1,nQ+2
- sumQ(j,i,iparm)=sumQ(j,i,iparm)/sumW(i,iparm)
- sumQsq(j,i,iparm)=sumQsq(j,i,iparm)/sumW(i,iparm) &
- -sumQ(j,i,iparm)**2
- sumEQ(j,i,iparm)=sumEQ(j,i,iparm)/sumW(i,iparm) &
- -sumQ(j,i,iparm)*sumE(i,iparm)
- enddo
- sumW(i,iparm)=-dlog(sumW(i,iparm))*(1.987D-3* &
- (startGridT+i*delta_T))+potEmin
- write (iout,'(f7.1,2f15.5,$)') startGridT+i*delta_T,&
- sumW(i,iparm),sumE(i,iparm)
- write (iout,'(f10.5,$)') (sumQ(j,i,iparm),j=1,nQ+2)
- write (iout,'(e15.5,$)') sumEsq(i,iparm)-sumEbis(i,iparm),&
- (sumQsq(j,i,iparm),j=1,nQ+2),(sumEQ(j,i,iparm),j=1,nQ+2)
- 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)
- write (34,*)
- call flush(34)
- enddo
- close(34)
- enddo
- if (histout) then
- do t=0,tmax
- if (hfin_ent(t).gt.0.0d0) then
- liczbaW=t
- jj = mod(liczbaW,nbin1)
- write (iout,'(f6.3,e20.10," ent")') dmin+(jj+0.5d0)*delta,&
- hfin_ent(t)
- if (histfile) write (ihist,'(f6.3,e20.10," ent")') &
- dmin+(jj+0.5d0)*delta,&
- hfin_ent(t)
- endif
- enddo
- if (histfile) close(ihist)
- endif
-
-#ifdef ZSCORE
-! Write data for zscore
- if (nslice.eq.1) then
- zscname=prefix(:ilen(prefix))//".zsc"
- else
- zscname=prefix(:ilen(prefix))//"_slice_"//licz2//".zsc"
- endif
-#if defined(AIX) || defined(PGI)
- open (izsc,file=prefix(:ilen(prefix))//'.zsc',position='append')
-#else
- open (izsc,file=prefix(:ilen(prefix))//'.zsc',access='append')
-#endif
- write (izsc,'("NQ=",i1," NPARM=",i1)') nQ,nParmSet
- do iparm=1,nParmSet
- write (izsc,'("NT=",i1)') nT_h(iparm)
- do ib=1,nT_h(iparm)
- write (izsc,'("TEMP=",f6.1," NR=",i2," SNK=",$)') &
- 1.0d0/(beta_h(ib,iparm)*1.987D-3),nR(ib,iparm)
- jj = min0(nR(ib,iparm),7)
- write (izsc,'(i8,$)') (snk(i,ib,iparm,islice),i=1,jj)
- write (izsc,'(a1,$)') (" ",i=22+8*jj+1,79)
- write (izsc,'("&")')
- if (nR(ib,iparm).gt.7) then
- do ii=8,nR(ib,iparm),9
- jj = min0(nR(ib,iparm),ii+8)
- write (izsc,'(i8,$)') (snk(i,ib,iparm,islice),i=ii,jj)
- write (izsc,'(a1,$') (" ",i=(jj-ii+1)*8+1,79)
- write (izsc,'("&")')
- enddo
- endif
- write (izsc,'("FI=",$)')
- jj=min0(nR(ib,iparm),7)
- write (izsc,'(f10.5,$)') (fi(i,ib,iparm),i=1,jj)
- write (izsc,'(a1,$)') (" ",i=3+10*jj+1,79)
- write (izsc,'("&")')
- if (nR(ib,iparm).gt.7) then
- do ii=8,nR(ib,iparm),9
- jj = min0(nR(ib,iparm),ii+8)
- write (izsc,'(f10.5,$)') (fi(i,ib,iparm),i=ii,jj)
- if (jj.eq.nR(ib,iparm)) then
- write (izsc,*)
- else
- write (izsc,'(a1,$)') (" ",i=10*(jj-ii+1)+1,79)
- write (izsc,'(t80,"&")')
- endif
- enddo
- endif
- do i=1,nR(ib,iparm)
- write (izsc,'("KH=",$)')
- write (izsc,'(f7.2,$)') (Kh(j,i,ib,iparm),j=1,nQ)
- write (izsc,'(" Q0=",$)')
- write (izsc,'(f7.5,$)') (q0(j,i,ib,iparm),j=1,nQ)
- write (izsc,*)
- enddo
- enddo
- enddo
- close(izsc)
-#endif
-#ifdef MPI
- endif
-#endif
- return
- end subroutine WHAMCALC
-!-----------------------------------------------------------------------------
- end module wham_calc
-
projects / unres4.git / blobdiff