parameter (MaxPoint=MaxStr,
& MaxPointProc=MaxStr_Proc)
double precision finorm_max,potfac,entmin,entmax,expfac,vf
+ double precision entfac_min,entfac_min_t
parameter (finorm_max=1.0d0)
integer islice
integer i,ii,j,jj,k,kk,l,m,ind,iter,t,tmax,ient,ientmax,iln
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:Max_GridT,Max_Parm),
& sumE_p(0:Max_GridT,Max_Parm),sumEsq_p(0:Max_GridT,Max_Parm),
& sumQ_p(MaxQ1,0:Max_GridT,Max_Parm),
& 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),
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",
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
+ & +wstrain*ehpb+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
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
+ & +wstrain*ehpb+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
! 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
+c#ifdef MPI
+c write (iout,*) "entfac_min before AllReduce",entfac_min
+c call MPI_AllReduce(entfac_min,entfac_min_t,1,
+c & MPI_DOUBLE_PRECISION,MPI_MAX,WHAM_COMM,IERROR)
+c entfac_min=entfac_min_t
+c write (iout,*) "entfac_min after AllReduce",entfac_min
+c#endif
+c#ifdef MPI
+c do t=1,scount(me)
+c entfac(t)=entfac(t)-entfac_min
+c enddo
+c#else
+c do t=1,ntot(islice)
+c entfac(t)=entfac(t)-entfac_min
+c enddo
+c#endif
+ 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
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))
+ write (iout,'(6f15.5)') (fi(i,ib,iparm),i=1,nR(ib,iparm))
+ write (iout,'(6f15.5)') (f(i,ib,iparm),i=1,nR(ib,iparm))
enddo
enddo
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
+ & +wstrain*ehpb+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
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
+ & +wstrain*ehpb+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
do i=1,nParmSet
write (iout,*) "Parameter set",i
do j=1,nT_h(i)
- write (iout,*) j,PotEmin_all(j,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
-#undef DEBUG
#ifdef MPI
do t=0,tmax
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)
-c write(iout,*) "first",temper,potEmin
- else if (betaT.lt.beta_h(nT_h(iparm),iparm)) then
- potEmin=potEmin_all(nT_h(iparm),iparm)
-c write (iout,*) "last",temper,potEmin
+ 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)
-c write (iout,*) "l",l,
-c & betaT,1.0d0/(1.987D-3*beta_h(l,iparm)),
-c & 1.0d0/(1.987D-3*beta_h(l+1,iparm)),temper,potEmin
+#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
-c write (iout,*) ib," PotEmin",potEmin
+#ifdef DEBUG
+ write (iout,*) "k",k," potEmin",potEmin
+#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
+ & +wstrain*ehpb+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
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
+ & +wstrain*ehpb+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
endif
enddo
endif
+
+c write (iout,*) "i",i," potEmin",potEmin
+
sumE(i,iparm)=sumE(i,iparm)/sumW(i,iparm)
sumEbis(i,iparm)=(startGridT+i*delta_T)*sumEbis(i,iparm)/
& sumW(i,iparm)
projects / unres.git / blobdiff