clean fix revert to previous vestion

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

diff --git a/source/wham/src-M/wham_calc1.F b/source/wham/src-M/wham_calc1.F
index a3d1658..dd76035 100644 (file)
--- a/source/wham/src-M/wham_calc1.F
+++ b/source/wham/src-M/wham_calc1.F
@@ -89,7 +89,7 @@ c      parameter (MaxHdim=200000)
       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,
-     &  eliptran
+     &  eliptran,etube
 
       integer ind_point(maxpoint),upindE,indE
       character*16 plik
@@ -239,7 +239,7 @@ 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,22)
+     &     (enetb(k,i,iparm),k=1,max_ene)
 #endif
           call restore_parm(iparm)
 #ifdef DEBUG
@@ -325,6 +325,8 @@ c              write (iout,*) 1.0d0/(beta_h(ib,iparm)*1.987D-3),ft
             esccor=enetb(19,i,iparm)
             edihcnstr=enetb(20,i,iparm)
             eliptran=enetb(22,i,iparm)
+            etube=enetb(25,i,iparm)
+
 
 #ifdef DEBUG
             write (iout,'(3i5,6f5.2,14f12.3)') i,ib,iparm,(ft(l),l=1,6),
@@ -343,7 +345,7 @@ c              write (iout,*) 1.0d0/(beta_h(ib,iparm)*1.987D-3),ft
      &      +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+wliptran*eliptran
+     &      +wbond*estr+wliptran*eliptran+wtube*Etube
              else
             etot=wsc*(evdw+ft(6)*evdw_t)+wscp*evdw2+ft(1)*welec*ees
      &      +wvdwpp*evdw1
@@ -353,7 +355,7 @@ c              write (iout,*) 1.0d0/(beta_h(ib,iparm)*1.987D-3),ft
      &      +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+wliptran*eliptran
+     &      +wbond*estr+wliptran*eliptran+wtube*Etube
              endif
 #else
       if (shield_mode.gt.0) then
@@ -365,7 +367,7 @@ c              write (iout,*) 1.0d0/(beta_h(ib,iparm)*1.987D-3),ft
      &      +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+wliptran*eliptran
+     &      +wbond*estr+wliptran*eliptran+wtube*Etube
             else
             etot=wsc*(evdw+ft(6)*evdw_t)+wscp*evdw2
      &      +ft(1)*welec*(ees+evdw1)
@@ -375,7 +377,7 @@ c              write (iout,*) 1.0d0/(beta_h(ib,iparm)*1.987D-3),ft
      &      +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+wliptran*eliptran
+     &      +wbond*estr+wliptran*eliptran+wtube*Etube
            endif
 
 #endif
@@ -693,6 +695,8 @@ c              write (iout,*) 1.0d0/(beta_h(ib,iparm)*1.987D-3),ft
           edihcnstr=enetb(20,i,iparm)
 C          edihcnstr=0.0d0
           eliptran=enetb(22,i,iparm)
+            etube=enetb(25,i,iparm)
+
 #ifdef SPLITELE
       if (shield_mode.gt.0) then
             etot=ft(1)*wsc*(evdw+ft(6)*evdw_t)+ft(1)*wscp*evdw2
@@ -704,7 +708,7 @@ C          edihcnstr=0.0d0
      &      +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+wliptran*eliptran
+     &      +wbond*estr+wliptran*eliptran+wtube*Etube
         else
             etot=wsc*(evdw+ft(6)*evdw_t)+wscp*evdw2+ft(1)*welec*ees
      &      +wvdwpp*evdw1
@@ -714,7 +718,7 @@ C          edihcnstr=0.0d0
      &      +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+wliptran*eliptran
+     &      +wbond*estr+wliptran*eliptran+wtube*Etube
       endif
 #else
       if (shield_mode.gt.0) then
@@ -726,7 +730,7 @@ C          edihcnstr=0.0d0
      &      +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+wliptran*eliptran
+     &      +wbond*estr+wliptran*eliptran+wtube*Etube
        else
             etot=wsc*(evdw+ft(6)*evdw_t)+wscp*evdw2
      &      +ft(1)*welec*(ees+evdw1)
@@ -736,7 +740,7 @@ C          edihcnstr=0.0d0
      &      +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+wliptran*eliptran
+     &      +wbond*estr+wliptran*eliptran+wtube*Etube
        endif
 
 #endif
@@ -847,23 +851,29 @@ c        ent=-dlog(entfac(t))
      &  WHAM_COMM,IERROR)
       call MPI_Allreduce(entmax_p,entmax,1,MPI_DOUBLE_PRECISION,MPI_MAX,
      &  WHAM_COMM,IERROR)
-      ientmax=entmax-entmin 
-      if (ientmax.gt.2000) ientmax=2000
+C      ientmax=entmax-entmin 
+C      if (ientmax.gt.2000) ientmax=2000
+      if ((-dlog(entmax)-entmin).lt.2000.0d0) then
+      ientmax=-dlog(entmax)-entmin
+      else
+       ientmax=2000
+      endif
       write (iout,*) "entmin",entmin," entmax",entmax," ientmax",ientmax
       call flush(iout)
       do t=1,scount(me1)
 c        ient=-dlog(entfac(t))-entmin
         ient=entfac(t)-entmin
-        if (ient.le.2000) histent_p(ient)=histent_p(ient)+1
+      write (iout,*) "ient",ient,entfac(t),entmin
+C        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
+C      call MPI_Allreduce(histent_p(0),histent(0),ientmax+1,MPI_INTEGER,
+C     &  MPI_SUM,WHAM_COMM,IERROR)
+C      if (me1.eq.Master) then
+C        write (iout,*) "Entropy histogram"
+C        do i=0,ientmax
+C          write(iout,'(f15.4,i10)') entmin+i,histent(i)
+C        enddo
+C      endif
 #else
       entmin=1.0d10
       entmax=-1.0d10
@@ -872,16 +882,19 @@ c        ient=-dlog(entfac(t))-entmin
         if (ent.lt.entmin) entmin=ent
         if (ent.gt.entmax) entmax=ent
       enddo
+      if ((-dlog(entmax)-entmin).lt.2000.0d0) then
       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
+      else
+       ientmax=2000
+      endif
+C      do t=1,ntot(islice)
+C        ient=entfac(t)-entmin
+C        if (ient.le.2000) histent(ient)=histent(ient)+1
+C      enddo
+C      write (iout,*) "Entropy histogram"
+C      do i=0,ientmax
+C        write(iout,'(2f15.4)') entmin+i,histent(i)
+C      enddo
 #endif
       
 #ifdef MPI
@@ -963,7 +976,8 @@ c          write (iout,'(2i5,20f8.2)') t,t,(enetb(k,t,iparm),k=1,18)
           esccor=enetb(19,t,iparm)
           edihcnstr=enetb(20,t,iparm)
 C          edihcnstr=0.0d0
-          eliptran=enetb(22,i,iparm)
+          eliptran=enetb(22,t,iparm)
+            etube=enetb(25,t,iparm)
 
           do k=0,nGridT
             betaT=startGridT+k*delta_T
@@ -1086,7 +1100,7 @@ c            write (iout,*) "ftbis",ftbis
      &      +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+wliptran*eliptran
+     &      +wbond*estr+wliptran*eliptran+wtube*Etube
             eprim=ftprim(1)*(ft(6)*evdw_t+evdw)
 C     &            +ftprim(6)*evdw_t
      &            +ftprim(1)*wscp*evdw2
@@ -1117,7 +1131,7 @@ C     &            +ftprim(6)*evdw_t
      &      +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+wliptran*eliptran
+     &      +wbond*estr+wliptran*eliptran+wtube*Etube
             eprim=ftprim(6)*evdw_t+ftprim(1)*welec*ees
      &            +ftprim(1)*wtor*etors+
      &            ftprim(3)*wcorr*ecorr+ftprim(4)*wcorr5*ecorr5+
@@ -1142,7 +1156,7 @@ C     &            +ftprim(6)*evdw_t
      &      +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+wliptran*eliptran
+     &      +wbond*estr+wliptran*eliptran+wtube*Etube
             eprim=ftprim(1)*(evdw+ft(6)*evdw_t)
      &           +ftprim(1)*welec*(ees+evdw1)
      &           +ftprim(1)*wtor*etors+
@@ -1169,7 +1183,7 @@ C     &            +ftprim(6)*evdw_t
      &      +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+wliptran*eliptran
+     &      +wbond*estr+wliptran*eliptran+wtube*Etube
             eprim=ftprim(6)*evdw_t+ftprim(1)*welec*(ees+evdw1)
      &           +ftprim(1)*wtor*etors+
      &            ftprim(3)*wcorr*ecorr+ftprim(4)*wcorr5*ecorr5+
@@ -1187,6 +1201,7 @@ C     &            +ftprim(6)*evdw_t
        endif
 
 #endif
+
             weight=dexp(-betaT*(etot-potEmin)+entfac(t))
 #ifdef DEBUG
             write (iout,*) "iparm",iparm," t",t," betaT",betaT,