5/29/2012 by Adam

[unres.git] / source / wham / src / wham_calc1.F

diff --git a/source/wham/src/wham_calc1.F b/source/wham/src/wham_calc1.F
index 5783bd4..a6044cd 100644 (file)
--- a/source/wham/src/wham_calc1.F
+++ b/source/wham/src/wham_calc1.F
@@ -38,6 +38,7 @@ c      parameter (MaxHdim=200000)
       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
@@ -48,7 +49,8 @@ c      parameter (MaxHdim=200000)
       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),
@@ -71,6 +73,7 @@ c      parameter (MaxHdim=200000)
      & 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),
@@ -368,6 +371,7 @@ c              write (iout,*) 1.0d0/(beta_h(ib,iparm)*1.987D-3),ft
 ! 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
@@ -393,10 +397,67 @@ c              write (iout,*) 1.0d0/(beta_h(ib,iparm)*1.987D-3),ft
             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)
@@ -404,17 +465,18 @@ c              write (iout,*) 1.0d0/(beta_h(ib,iparm)*1.987D-3),ft
               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
@@ -432,10 +494,12 @@ c              write (iout,*) 1.0d0/(beta_h(ib,iparm)*1.987D-3),ft
           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)
@@ -455,7 +519,7 @@ c              write (iout,*) 1.0d0/(beta_h(ib,iparm)*1.987D-3),ft
         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
@@ -674,7 +738,6 @@ C Determine the minimum energes for all parameter sets and temperatures
       enddo
       write (iout,*) "potEmin_min",potEmin_min
 #endif
-#undef DEBUG
 
 #ifdef MPI
       do t=0,tmax