projects / unres.git / blobdiff
? search:
summary | shortlog | log | commit | commitdiff | tree
raw | inline | side by side
dfa wham and cluster corrections
[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 986f6a3..41f42bc 100644 (file)
--- a/source/wham/src/wham_calc1.F
+++ b/source/wham/src/wham_calc1.F
@@ -356,7 +356,10 @@ c              write (iout,*) 1.0d0/(beta_h(ib,iparm)*1.987D-3),ft
             call enerprint(energia(0),fT)
             endif
 #endif
-            if (homol_nset.le.1) then
+#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
@@ -364,7 +367,7 @@ c              write (iout,*) 1.0d0/(beta_h(ib,iparm)*1.987D-3),ft
      &          -beta_h(ib,iparm)*(etot+Econstr)
 #ifdef DEBUG
               write (iout,'(4i5,4e15.5)') i,kk,ib,iparm,
-     &         etot,v(i,kk,ib,iparm)
+     &         etot,Econstr,v(i,kk,ib,iparm)
 #endif
             enddo ! kk
 
@@ -558,8 +561,8 @@ c#endif
             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
 
@@ -692,6 +695,8 @@ c              write (iout,*) 1.0d0/(beta_h(ib,iparm)*1.987D-3),ft
             esccor=enetb(19,i,iparm)
             edihcnstr=enetb(20,i,iparm)
             ehomology_constr=enetb(22,i,iparm)
+            if (homol_nset.gt.1) 
+     &       ehomology_constr=waga_homology(ihset)*ehomology_constr
             edfadis=enetb(23,i,iparm)
             edfator=enetb(24,i,iparm)
             edfanei=enetb(25,i,iparm)
@@ -731,8 +736,8 @@ c              write (iout,*) 1.0d0/(beta_h(ib,iparm)*1.987D-3),ft
             write (iout,*) "i",i," ib",ib,
      &      " temp",1.0d0/(1.987d-3*beta_h(ib,iparm))," etot",etot,
      &      " entfac",entfac(i)," ecorr",etot-entfac(i)/beta_h(ib,iparm)
-            etot=etot-entfac(i)/beta_h(ib,iparm)
 #endif
+            etot=etot-entfac(i)/beta_h(ib,iparm)
             if(etot.lt.potEmin_all(ib,iparm)) potEmin_all(ib,iparm)=etot
 #ifdef DEBUG
             write (iout,*) "efree",etot," potEmin",potEmin_all(ib,iparm)
@@ -772,7 +777,8 @@ C Determine the minimum energes for all parameter sets and temperatures
       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
@@ -950,10 +956,13 @@ c      write (iout,*) "me1",me1," scount",scount(me1)
           esccor=enetb(19,t,iparm)
           edihcnstr=enetb(20,t,iparm)
           ehomology_constr=enetb(22,t,iparm)
+          if (homol_nset.gt.1)
+     &       ehomology_constr=waga_homology(ihset)*ehomology_constr
           edfadis=enetb(23,t,iparm)
           edfator=enetb(24,t,iparm)
           edfanei=enetb(25,t,iparm)
           edfabet=enetb(26,t,iparm)
+c          do k=0,nGridT
           do k=0,nGridT
             betaT=startGridT+k*delta_T
             temper=betaT
@@ -1035,24 +1044,38 @@ 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)
-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
@@ -1063,6 +1086,8 @@ c            write (iout,*) ib," PotEmin",potEmin
      &      +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
+     &      +wdfa_dist*edfadis
+     &      +wdfa_tor*edfator+wdfa_nei*edfanei+wdfa_beta*edfabet
             eprim=ftprim(6)*evdw_t+ftprim(1)*welec*ees
      &            +ftprim(1)*wtor*etors+
      &            ftprim(3)*wcorr*ecorr+ftprim(4)*wcorr5*ecorr5+
@@ -1086,6 +1111,8 @@ c            write (iout,*) ib," PotEmin",potEmin
      &      +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
+     &      +wdfa_dist*edfadis
+     &      +wdfa_tor*edfator+wdfa_nei*edfanei+wdfa_beta*edfabet
             eprim=ftprim(6)*evdw_t+ftprim(1)*welec*(ees+evdw1)
      &           +ftprim(1)*wtor*etors+
      &            ftprim(3)*wcorr*ecorr+ftprim(4)*wcorr5*ecorr5+
@@ -1376,6 +1403,7 @@ c            write (iout,*) ib," PotEmin",potEmin
       else
         write (iout,'(a,i3)') "Parameter set",iparm
       endif
+c      do i=0,NGridT
       do i=0,NGridT
         betaT=1.0d0/(1.987D-3*(startGridT+i*delta_T))
         if (betaT.ge.beta_h(1,iparm)) then
@@ -1391,6 +1419,9 @@ c            write (iout,*) ib," PotEmin",potEmin
             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)
UNRESPACK 3.x