small changes in cation cation interactions and wham correction
[unres4.git] / source / wham / enecalc.F90
index b8fd48f..f1546ad 100644 (file)
@@ -2,6 +2,7 @@
 !-----------------------------------------------------------------------------
       use io_units
       use wham_data
+      use control_data, only: tor_mode
 !
       use geometry_data, only:nres,boxxsize,boxysize,boxzsize
       use energy_data
@@ -555,7 +556,7 @@ write(iout,*)"enecalc_ i ntot",i,ntot
           (vbld(j).lt.2.0d0 .or. vbld(j).gt.5.0d0)) then
           if (iprint.gt.0) &
           write (iout,*) "Bad CA-CA bond length",j," ",vbld(j),&
-            " for conformation",ii
+            " for conformation",ii,mnum
           if (iprint.gt.1) then
             write (iout,*) "The Cartesian geometry is:"
             write (iout,'(8f10.5)') ((c(l,k),l=1,3),k=1,nres)
@@ -675,6 +676,8 @@ write(iout,*)"enecalc_ i ntot",i,ntot
       ww_all(16,iparm)=wvdwpp
       ww_all(17,iparm)=wbond
       ww_all(19,iparm)=wsccor
+      ww_all(42,iparm)=wcatprot
+      ww_all(41,iparm)=wcatcat
 ! Store bond parameters
       vbldp0_all(iparm)=vbldp0
       akp_all(iparm)=akp
@@ -709,6 +712,7 @@ write(iout,*)"enecalc_ i ntot",i,ntot
         sigc0_all(i,iparm)=sigc0(i)
       enddo
 #else
+      IF (tor_mode.eq.0) THEN
       nthetyp_all(iparm)=nthetyp
       ntheterm_all(iparm)=ntheterm
       ntheterm2_all(iparm)=ntheterm2
@@ -760,6 +764,9 @@ write(iout,*)"enecalc_ i ntot",i,ntot
         enddo
       enddo
       enddo
+      ELSE
+      write(iout,*) "Need storing parameters"
+      ENDIF
 #endif
 #ifdef CRYST_SC
 ! Store the sidechain rotamer parameters
@@ -787,6 +794,7 @@ write(iout,*)"enecalc_ i ntot",i,ntot
         enddo
       enddo
 #endif
+      IF (tor_mode.eq.0) THEN
 ! Store the torsional parameters
       do iblock=1,2
       do i=-ntortyp+1,ntortyp-1
@@ -829,7 +837,7 @@ write(iout,*)"enecalc_ i ntot",i,ntot
       enddo
       enddo
 ! Store parameters of the cumulants
-      do i=-nloctyp,nloctyp
+      do i=1,nloctyp
         do j=1,2
           b1_all(j,i,iparm)=b1(j,i)
           b1tilde_all(j,i,iparm)=b1tilde(j,i)
@@ -845,6 +853,9 @@ write(iout,*)"enecalc_ i ntot",i,ntot
           enddo
         enddo
       enddo
+      ELSE
+       write(iout,*) "NEED storing parameters"
+      ENDIF
 ! Store the parameters of electrostatic interactions
       do i=1,2
         do j=1,2
@@ -947,6 +958,8 @@ write(iout,*)"end of store_parm"
       wvdwpp=ww_all(16,iparm)
       wbond=ww_all(17,iparm)
       wsccor=ww_all(19,iparm)
+      wcatprot=ww_all(42,iparm)
+      wcatcat=ww_all(41,iparm)
 ! Restore bond parameters
       vbldp0=vbldp0_all(iparm)
       akp=akp_all(iparm)
@@ -981,6 +994,7 @@ write(iout,*)"end of store_parm"
         sigc0(i)=sigc0_all(i,iparm)
       enddo
 #else
+      IF (tor_mode.eq.0) THEN
       nthetyp=nthetyp_all(iparm)
       ntheterm=ntheterm_all(iparm)
       ntheterm2=ntheterm2_all(iparm)
@@ -1032,6 +1046,9 @@ write(iout,*)"end of store_parm"
         enddo
       enddo
       enddo
+      ELSE
+      write (iout,*) "Need storing parameters"
+      ENDIF
 #endif
 ! Restore the sidechain rotamer parameters
 #ifdef CRYST_SC
@@ -1058,6 +1075,7 @@ write(iout,*)"end of store_parm"
         enddo
       enddo
 #endif
+      IF (tor_mode.eq.0) THEN
 ! Restore the torsional parameters
       do iblock=1,2
       do i=-ntortyp+1,ntortyp-1
@@ -1099,8 +1117,9 @@ write(iout,*)"end of store_parm"
         enddo
       enddo
       enddo
+     
 ! Restore parameters of the cumulants
-      do i=-nloctyp,nloctyp
+      do i=1,nloctyp
         do j=1,2
           b1(j,i)=b1_all(j,i,iparm)
           b1tilde(j,i)=b1tilde_all(j,i,iparm)
@@ -1116,6 +1135,9 @@ write(iout,*)"end of store_parm"
           enddo
         enddo
       enddo
+      ELSE
+      write(iout,*) "need storing parameters"
+      ENDIF
 ! Restore the parameters of electrostatic interactions
       do i=1,2
         do j=1,2
@@ -1207,7 +1229,10 @@ write(iout,*)"end of store_parm"
       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,tt, &
-            ecation_prot, ecationcation
+            ecation_prot, ecationcation,evdwpp,eespp,evdwpsb,eelpsb,evdwsb, &
+            eelsb,estr_nucl,ebe_nucl,esbloc,etors_nucl,etors_d_nucl,&
+            ecorr_nucl,ecorr3_nucl,escbase, epepbase,escpho, epeppho
+
       integer :: i,ii,ik,iproc,iscor,j,k,l,ib,iparm,iprot,nlist
       real(kind=8) :: qfree,sumprob,eini,efree,rmsdev
       character(len=80) :: bxname
@@ -1358,6 +1383,27 @@ write(iout,*)"end of store_parm"
             edihcnstr=enetb(19,i,iparm)
             ecationcation=enetb(42,i,iparm)
             ecation_prot=enetb(41,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)
+!      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 &
@@ -1390,7 +1436,14 @@ write(iout,*)"end of store_parm"
             +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) &
@@ -1400,7 +1453,14 @@ write(iout,*)"end of store_parm"
             +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 MPI
@@ -1418,13 +1478,22 @@ write(iout,*)"end of store_parm"
 #endif
           enddo   ! i
 #ifdef MPI
+          scount_(:)=0
           do i=1,scount(me1)
             Fdimless_(i)=Fdimless(i)
           enddo
+          write(iout,*) "before gather Fdimless"
+          write(iout,*) scount(me),scount_(0),idispl(0)
+          write (iout,*) "added update of scount_"
+         call MPI_AllGather( scount(me), 1, MPI_INTEGER, scount_(0), 1, &
+        MPI_INTEGER, WHAM_COMM, IERROR)
+
+
           call MPI_Gatherv(Fdimless_(1),scount(me),&
            MPI_REAL,Fdimless(1),&
            scount_(0),idispl(0),MPI_REAL,Master,&
            WHAM_COMM, IERROR)
+          write(iout,*) "after gather Fdimless"
 #ifdef DEBUG
           call MPI_Gatherv(potE(1,iparm),scount_(me),&
            MPI_DOUBLE_PRECISION,potE(1,iparm),&