working wham and cluster for NARES_UNRES

diff --git a/source/cluster/io_clust.F90 b/source/cluster/io_clust.F90
index 197a1a7..a75c27e 100644 (file)
--- a/source/cluster/io_clust.F90
+++ b/source/cluster/io_clust.F90
@@ -1448,8 +1448,28 @@
       call reada(weightcard,'WCATPROT',wcatprot,0.0d0)
       call reada(weightcard,'TEMP0',temp0,300.0d0)   !!! el
       if (index(weightcard,'SOFT').gt.0) ipot=6
+      call reada(weightcard,'WVDWPP_NUCL',wvdwpp_nucl,0.0D0)
+      call reada(weightcard,'WELPP',welpp,0.0d0)
+      call reada(weightcard,'WVDWPSB',wvdwpsb,0.0d0)
+      call reada(weightcard,'WELPSB',welpsb,0.0D0)
+      call reada(weightcard,'WVDWSB',wvdwsb,0.0d0)
+      call reada(weightcard,'WELSB',welsb,0.0D0)
+      call reada(weightcard,'WBOND_NUCL',wbond_nucl,0.0D0)
+      call reada(weightcard,'WANG_NUCL',wang_nucl,0.0D0)
+      call reada(weightcard,'WSBLOC',wsbloc,0.0D0)
+      call reada(weightcard,'WTOR_NUCL',wtor_nucl,0.0D0)
+!      print *,"KUR..",wtor_nucl
+      call reada(weightcard,'WTORD_NUCL',wtor_d_nucl,0.0D0)
+      call reada(weightcard,'WCORR_NUCL',wcorr_nucl,0.0D0)
+      call reada(weightcard,'WCORR3_NUC',wcorr3_nucl,0.0D0)
 ! 12/1/95 Added weight for the multi-body term WCORR
       call reada(weightcard,'WCORRH',wcorr,1.0D0)
+      call reada(weightcard,'WSCBASE',wscbase,0.0D0)
+      call reada(weightcard,'WPEPBASE',wpepbase,1.0d0)
+      call reada(weightcard,'WSCPHO',wscpho,0.0d0)
+      call reada(weightcard,'WPEPPHO',wpeppho,0.0d0)
+
+
       if (wcorr4.gt.0.0d0) wcorr=wcorr4
       weights(1)=wsc
       weights(2)=wscp
@@ -1471,6 +1491,26 @@
       weights(18)=scal14
 !el      weights(19)=wsccor !!!!!!!!!!!!!!!!
       weights(21)=wsccor
+          weights(26)=wvdwpp_nucl
+          weights(27)=welpp
+          weights(28)=wvdwpsb
+          weights(29)=welpsb
+          weights(30)=wvdwsb
+          weights(31)=welsb
+          weights(32)=wbond_nucl
+          weights(33)=wang_nucl
+          weights(34)=wsbloc
+          weights(35)=wtor_nucl
+          weights(36)=wtor_d_nucl
+          weights(37)=wcorr_nucl
+          weights(38)=wcorr3_nucl
+          weights(41)=wcatcat
+          weights(42)=wcatprot
+          weights(46)=wscbase
+          weights(47)=wscpho
+          weights(48)=wpeppho
+
+
       write (iout,10) wsc,wscp,welec,wvdwpp,wbond,wang,wscloc,wtor,&
         wtor_d,wstrain,wel_loc,wcorr,wcorr5,wcorr6,wturn3,&
         wturn4,wturn6,wsccor
@@ -1729,6 +1769,8 @@ write(iout,*)"po setup var"
       open (isidep1,file=sidepname,status="old")
       call getenv('SCCORPAR',sccorname)
       open (isccor,file=sccorname,status="old")
+      call getenv('BONDPAR_NUCL',bondname_nucl)
+      open (ibond_nucl,file=bondname_nucl,status='old')
       call getenv('THETPAR_NUCL',thetname_nucl)
       open (ithep_nucl,file=thetname_nucl,status='old')
       call getenv('ROTPAR_NUCL',rotname_nucl)
@@ -1763,6 +1805,9 @@ write(iout,*)"po setup var"
 !
       call getenv('SCPPAR',scpname)
       open (iscpp,file=scpname,status='old')
+      call getenv('SCPPAR_NUCL',scpname_nucl)
+      open (iscpp_nucl,file=scpname_nucl,status='old')
+
 #endif
       return
       end subroutine openunits

diff --git a/source/wham/enecalc.F90 b/source/wham/enecalc.F90
index 1568a9a..9adc19a 100644 (file)
--- a/source/wham/enecalc.F90
+++ b/source/wham/enecalc.F90
@@ -1396,6 +1396,14 @@ write(iout,*)"end of store_parm"
             etors_d_nucl =      enetb(36,i,iparm)
             ecorr_nucl =      enetb(37,i,iparm)
             ecorr3_nucl =      enetb(38,i,iparm)
+            epeppho=   enetb(49,i,iparm)
+            escpho=    enetb(48,i,iparm)
+            epepbase=  enetb(47,i,iparm)
+            escbase=   enetb(46,i,iparm)
+!      wscbase=ww(46)
+!      wpepbase=ww(47)
+!      wscpho=ww(48)
+!      wpeppho=ww(49)
 
 !#ifdef SPLITELE
 !            etot=wsc*(evdw+ft(6)*evdw_t)+wscp*evdw2+ft(1)*welec*ees &

diff --git a/source/wham/wham_calc.F90 b/source/wham/wham_calc.F90
index d1ed2f7..e2bb757 100644 (file)
--- a/source/wham/wham_calc.F90
+++ b/source/wham/wham_calc.F90
@@ -103,7 +103,10 @@
       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, &
-        ecationcation,ecation_prot
+        ecationcation,ecation_prot, evdwpp,eespp ,evdwpsb,eelpsb, &
+        evdwsb, eelsb, estr_nucl,ebe_nucl,esbloc,etors_nucl,etors_d_nucl,&
+        ecorr_nucl, ecorr3_nucl,epeppho, escpho, epepbase,escbase
+
 
       integer :: ind_point(maxpoint),upindE,indE
       character(len=16) :: plik
@@ -345,6 +348,24 @@
             edihcnstr=enetb(19,i,iparm)
           ecationcation=enetb(41,i,iparm)
           ecation_prot=enetb(42,i,iparm)
+            evdwpp =      enetb(26,i,iparm)
+            eespp =      enetb(27,i,iparm)
+            evdwpsb =      enetb(28,i,iparm)
+            eelpsb =      enetb(29,i,iparm)
+            evdwsb =      enetb(30,i,iparm)
+            eelsb =      enetb(31,i,iparm)
+            estr_nucl =      enetb(32,i,iparm)
+            ebe_nucl =      enetb(33,i,iparm)
+            esbloc =      enetb(34,i,iparm)
+            etors_nucl =      enetb(35,i,iparm)
+            etors_d_nucl =      enetb(36,i,iparm)
+            ecorr_nucl =      enetb(37,i,iparm)
+            ecorr3_nucl =      enetb(38,i,iparm)
+            epeppho=   enetb(49,i,iparm)
+            escpho=    enetb(48,i,iparm)
+            epepbase=  enetb(47,i,iparm)
+            escbase=   enetb(46,i,iparm)
+
 #ifdef DEBUG
             write (iout,'(3i5,6f5.2,15f12.3)') i,ib,iparm,(ft(l),l=1,6),&
              evdw+evdw_t,evdw2,ees,evdw1,ecorr,eel_loc,estr,ebe,escloc,&
@@ -382,7 +403,15 @@
             +wturn3*eello_turn3 &
             +wturn6*eello_turn6+wel_loc*eel_loc &
             +edihcnstr+wtor_d*etors_d+wsccor*esccor &
-            +wbond*estr+wcatprot*ecation_prot+wcatcat*ecationcation
+            +wbond*estr+wcatprot*ecation_prot+wcatcat*ecationcation&
+            +wbond_nucl*estr_nucl+wang_nucl*ebe_nucl&
+            +wvdwpp_nucl*evdwpp+welpp*eespp+wvdwpsb*evdwpsb+welpsb*eelpsb&
+            +wvdwsb*evdwsb+welsb*eelsb+wsbloc*esbloc+wtor_nucl*etors_nucl&
+            +wtor_d_nucl*etors_d_nucl+wcorr_nucl*ecorr_nucl+wcorr3_nucl*ecorr3_nucl&
+            +wscbase*escbase&
+            +wpepbase*epepbase+wscpho*escpho+wpeppho*epeppho
+
+
 #else
             etot=wsc*evdw+wscp*evdw2 &
             +welec*(ees+evdw1) &
@@ -392,7 +421,15 @@
             +wturn3*eello_turn3 &
             +wturn6*eello_turn6+wel_loc*eel_loc+edihcnstr &
             +wtor_d*etors_d+wsccor*esccor &
-            +wbond*estr+wcatprot*ecation_prot+wcatcat*ecationcation
+            +wbond*estr+wcatprot*ecation_prot+wcatcat*ecationcation&
+            +wbond_nucl*estr_nucl+wang_nucl*ebe_nucl&
+            +wvdwpp_nucl*evdwpp+welpp*eespp+wvdwpsb*evdwpsb+welpsb*eelpsb&
+            +wvdwsb*evdwsb+welsb*eelsb+wsbloc*esbloc+wtor_nucl*etors_nucl&
+            +wtor_d_nucl*etors_d_nucl+wcorr_nucl*ecorr_nucl+wcorr3_nucl*ecorr3_nucl&
+            +wscbase*escbase&
+            +wpepbase*epepbase+wscpho*escpho+wpeppho*epeppho
+
+
 #endif
 
 #ifdef DEBUG
@@ -757,6 +794,24 @@
           edihcnstr=0.0d0
           ecationcation=enetb(41,t,iparm)
           ecation_prot=enetb(42,t,iparm)
+            evdwpp =      enetb(26,t,iparm)
+            eespp =      enetb(27,t,iparm)
+            evdwpsb =      enetb(28,t,iparm)
+            eelpsb =      enetb(29,t,iparm)
+            evdwsb =      enetb(30,t,iparm)
+            eelsb =      enetb(31,t,iparm)
+            estr_nucl =      enetb(32,t,iparm)
+            ebe_nucl =      enetb(33,t,iparm)
+            esbloc =      enetb(34,t,iparm)
+            etors_nucl =      enetb(35,t,iparm)
+            etors_d_nucl =      enetb(36,t,iparm)
+            ecorr_nucl =      enetb(37,t,iparm)
+            ecorr3_nucl =      enetb(38,t,iparm)
+            epeppho=   enetb(49,t,iparm)
+            escpho=    enetb(48,t,iparm)
+            epepbase=  enetb(47,t,iparm)
+            escbase=   enetb(46,t,iparm)
+
 
           do k=0,nGridT
             betaT=startGridT+k*delta_T
@@ -848,20 +903,30 @@
             +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+wcatprot*ecation_prot+wcatcat*ecationcation
+            +wbond*estr+wcatprot*ecation_prot+wcatcat*ecationcation &
+            +wbond_nucl*estr_nucl+wang_nucl*ebe_nucl&
+            +wvdwpp_nucl*evdwpp+welpp*eespp+wvdwpsb*evdwpsb+welpsb*eelpsb&
+            +wvdwsb*evdwsb+welsb*eelsb+wsbloc*esbloc+wtor_nucl*etors_nucl&
+            *ft(1)+wtor_d_nucl*ft(2)*etors_d_nucl+wcorr_nucl*ecorr_nucl &
+            +wcorr3_nucl*ecorr3_nucl&
+            +wscbase*escbase&
+            +wpepbase*epepbase+wscpho*escpho+wpeppho*epeppho
+
             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
+                 ftprim(1)*wsccor*esccor +ftprim(1)*wtor_nucl*etors_nucl&
+                 +wtor_d_nucl*ftprim(2)*etors_d_nucl
             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
+                 ftbis(1)*wsccor*esccor+ftbis(1)*wtor_nucl*etors_nucl&
+                 +wtor_d_nucl*ftbis(2)*etors_d_nucl
 #else
             etot=wsc*(evdw+ft(6)*evdw_t)+wscp*evdw2 &
             +ft(1)*welec*(ees+evdw1) &
@@ -871,20 +936,31 @@
             +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+wcatprot*ecation_prot+wcatcat*ecationcation
+            +wbond*estr+wcatprot*ecation_prot+wcatcat*ecationcation&
+            +wbond_nucl*estr_nucl+wang_nucl*ebe_nucl&
+            +wvdwpp_nucl*evdwpp+welpp*eespp+wvdwpsb*evdwpsb+welpsb*eelpsb&
+            +wvdwsb*evdwsb+welsb*eelsb+wsbloc*esbloc+wtor_nucl*etors_nucl&
+            *ft(1)+wtor_d_nucl*ft(2)*etors_d_nucl+wcorr_nucl*ecorr_nucl &
+            +wcorr3_nucl*ecorr3_nucl&
+            +wscbase*escbase&
+            +wpepbase*epepbase+wscpho*escpho+wpeppho*epeppho
+
             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
+                 ftprim(1)*wsccor*esccor+ftprim(1)*wtor_nucl*etors_nucl&
+                 +wtor_d_nucl*ftprim(2)*etors_d_nucl
             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
+                 ftprim(1)*wsccor*esccor+ftbis(1)*wtor_nucl*etors_nucl&
+                 +wtor_d_nucl*ftbis(2)*etors_d_nucl
+
 #endif
             weight=dexp(-betaT*(etot-potEmin)+entfac(t))
 #ifdef DEBUG