X-Git-Url: http://mmka.chem.univ.gda.pl/gitweb/?a=blobdiff_plain;f=source%2Fcluster%2Fwham%2Fsrc-M-SAXS-homology%2Fenergy_p_new.F;h=c2d7f8556ef7caf5776bda05d7b006529b9b03d3;hb=5836ecdab5a8b95f079bbf6e07374dee3fce8a26;hp=f599f70476911880b993dd14f6d47603d1d3f888;hpb=a05e8996b6fb955dc21aef7a85db21ee73d0a9d3;p=unres.git

diff --git a/source/cluster/wham/src-M-SAXS-homology/energy_p_new.F b/source/cluster/wham/src-M-SAXS-homology/energy_p_new.F
index f599f70..c2d7f85 100644
--- a/source/cluster/wham/src-M-SAXS-homology/energy_p_new.F
+++ b/source/cluster/wham/src-M-SAXS-homology/energy_p_new.F
@@ -130,8 +130,10 @@ C
 
       if (wliptran.gt.0) then
         call Eliptransfer(eliptran)
+      else
+        eliptran=0.0d0
       endif
-
+#ifdef FOURBODY
 C 
 C 12/1/95 Multi-body terms
 C
@@ -153,6 +155,7 @@ c         write (iout,*) ecorr,ecorr5,ecorr6,eturn6
 c         write (iout,*) "Calling multibody_hbond"
          call multibody_hb(ecorr,ecorr5,ecorr6,n_corr,n_corr1)
       endif
+#endif
 c      write (iout,*) "NSAXS",nsaxs
       if (nsaxs.gt.0 .and. saxs_mode.eq.0) then
         call e_saxs(Esaxs_constr)
@@ -189,8 +192,10 @@ c      write(iout,*)'edfab is finished!', wdfa_beta,edfabet
      & +welec*fact(1)*ees
      & +fact(1)*wvdwpp*evdw1
      & +wang*ebe+wtor*fact(1)*etors+wscloc*escloc
-     & +wstrain*ehpb+wcorr*fact(3)*ecorr+wcorr5*fact(4)*ecorr5
-     & +wcorr6*fact(5)*ecorr6+wturn4*fact(3)*eello_turn4
+     & +wstrain*ehpb
+     & +wcorr*fact(3)*ecorr+wcorr5*fact(4)*ecorr5
+     & +wcorr6*fact(5)*ecorr6
+     & +wturn4*fact(3)*eello_turn4
      & +wturn3*fact(2)*eello_turn3+wturn6*fact(5)*eturn6
      & +wel_loc*fact(2)*eel_loc+edihcnstr+wtor_d*fact(2)*etors_d
      & +wbond*estr+wsccor*fact(1)*esccor+ethetacnstr
@@ -505,10 +510,17 @@ C     Bartek
 #ifdef SPLITELE
       write(iout,10) evdw,wsc,evdw2,wscp,ees,welec*fact(1),evdw1,wvdwpp,
      &  estr,wbond,ebe,wang,escloc,wscloc,etors,wtor*fact(1),
-     &  etors_d,wtor_d*fact(2),ehpb,wstrain,ecorr,wcorr*fact(3),
-     &  ecorr5,wcorr5*fact(4),ecorr6,wcorr6*fact(5),eel_loc,
+     &  etors_d,wtor_d*fact(2),ehpb,wstrain,
+#ifdef FOURBODY
+     &  ecorr,wcorr*fact(3),
+     &  ecorr5,wcorr5*fact(4),ecorr6,wcorr6*fact(5),
+#endif
+     &  eel_loc,
      &  wel_loc*fact(2),eello_turn3,wturn3*fact(2),
-     &  eello_turn4,wturn4*fact(3),eello_turn6,wturn6*fact(5),
+     &  eello_turn4,wturn4*fact(3),
+#ifdef FOURBODY
+     &  eello_turn6,wturn6*fact(5),
+#endif
      &  esccor,wsccor*fact(1),edihcnstr,
      &  ethetacnstr,ebr*nss,Uconst,wumb,eliptran,wliptran,Eafmforc,
      &  etube,wtube,esaxs,wsaxs,ehomology_constr,
@@ -527,13 +539,17 @@ C     Bartek
      & 'ETORSD=',1pE16.6,' WEIGHT=',1pE16.6,' (double torsional)'/
      & 'EHBP=  ',1pE16.6,' WEIGHT=',1pE16.6,
      & ' (SS bridges & dist. cnstr.)'/
+#ifdef FOURBODY
      & 'ECORR4=',1pE16.6,' WEIGHT=',1pE16.6,' (multi-body)'/
      & 'ECORR5=',1pE16.6,' WEIGHT=',1pE16.6,' (multi-body)'/
      & 'ECORR6=',1pE16.6,' WEIGHT=',1pE16.6,' (multi-body)'/
+#endif
      & 'EELLO= ',1pE16.6,' WEIGHT=',1pE16.6,' (electrostatic-local)'/
      & 'ETURN3=',1pE16.6,' WEIGHT=',1pE16.6,' (turns, 3rd order)'/
      & 'ETURN4=',1pE16.6,' WEIGHT=',1pE16.6,' (turns, 4th order)'/
+#ifdef FOURBODY
      & 'ETURN6=',1pE16.6,' WEIGHT=',1pE16.6,' (turns, 6th order)'/
+#endif
      & 'ESCCOR=',1pE16.6,' WEIGHT=',1pE16.6,' (backbone-rotamer corr)'/
      & 'EDIHC= ',1pE16.6,' (virtual-bond dihedral angle restraints)'/
      & 'ETHETC=',1pE16.6,' (virtual-bond angle restraints)'/
@@ -553,10 +569,16 @@ C     Bartek
 #else
       write (iout,10) evdw,wsc,evdw2,wscp,ees,welec*fact(1),
      &  estr,wbond,ebe,wang,escloc,wscloc,etors,wtor*fact(1),
-     &  etors_d,wtor_d*fact(2),ehpb,wstrain,ecorr,wcorr*fact(3),
+     &  etors_d,wtor_d*fact(2),ehpb,
+#ifdef FOURBODY
+     &  wstrain,ecorr,wcorr*fact(3),
      &  ecorr5,wcorr5*fact(4),ecorr6,wcorr6*fact(5),
+#endif
      &  eel_loc,wel_loc*fact(2),eello_turn3,wturn3*fact(2),
-     &  eello_turn4,wturn4*fact(3),eello_turn6,wturn6*fact(5),
+     &  eello_turn4,wturn4*fact(3),
+#ifdef FOURBODY
+     &  eello_turn6,wturn6*fact(5),
+#endif
      &  esccor,wsccor*fact(1),edihcnstr,
      &  ethetacnstr,ebr*nss,Uconst,wumb,eliptran,wliptran,Eafmforc,
      &  etube,wtube,esaxs,wsaxs,ehomology_constr,
@@ -574,13 +596,17 @@ C     Bartek
      & 'ETORSD=',1pE16.6,' WEIGHT=',1pE16.6,' (double torsional)'/
      & 'EHBP=  ',1pE16.6,' WEIGHT=',1pE16.6,
      & ' (SS bridges & dist. restr.)'/
+#ifdef FOURBODY
      & 'ECORR4=',1pE16.6,' WEIGHT=',1pE16.6,' (multi-body)'/
      & 'ECORR5=',1pE16.6,' WEIGHT=',1pE16.6,' (multi-body)'/
      & 'ECORR6=',1pE16.6,' WEIGHT=',1pE16.6,' (multi-body)'/
+#endif
      & 'EELLO= ',1pE16.6,' WEIGHT=',1pE16.6,' (electrostatic-local)'/
      & 'ETURN3=',1pE16.6,' WEIGHT=',1pE16.6,' (turns, 3rd order)'/
      & 'ETURN4=',1pE16.6,' WEIGHT=',1pE16.6,' (turns, 4th order)'/
+#ifdef FOURBODY
      & 'ETURN6=',1pE16.6,' WEIGHT=',1pE16.6,' (turns, 6th order)'/
+#endif
      & 'ESCCOR=',1pE16.6,' WEIGHT=',1pE16.6,' (backbone-rotamer corr)'/
      & 'EDIHC= ',1pE16.6,' (virtual-bond dihedral angle restraints)'/
      & 'ETHETC=',1pE16.6,' (virtual-bond angle restraints)'/
@@ -619,7 +645,10 @@ C
       include 'COMMON.SBRIDGE'
       include 'COMMON.NAMES'
       include 'COMMON.IOUNITS'
+#ifdef FOURBODY
       include 'COMMON.CONTACTS'
+      include 'COMMON.CONTMAT'
+#endif
       dimension gg(3)
       integer icant
       external icant
@@ -658,6 +687,10 @@ cd   &                  'iend=',iend(i,iint)
 C Change 12/1/95 to calculate four-body interactions
             rij=xj*xj+yj*yj+zj*zj
             rrij=1.0D0/rij
+            sqrij=dsqrt(rij)
+            sss1=sscale(sqrij)
+            if (sss1.eq.0.0d0) cycle
+            sssgrad1=sscagrad(sqrij)
 c           write (iout,*)'i=',i,' j=',j,' itypi=',itypi,' itypj=',itypj
             eps0ij=eps(itypi,itypj)
             fac=rrij**expon2
@@ -677,15 +710,16 @@ cd   &        restyp(itypi),i,restyp(itypj),j,aa(itypi,itypj),
 cd   &        bb(itypi,itypj),1.0D0/dsqrt(rrij),evdwij,epsi,sigm,
 cd   &        (c(k,i),k=1,3),(c(k,j),k=1,3)
             if (bb.gt.0.0d0) then
-              evdw=evdw+evdwij
+              evdw=evdw+sss1*evdwij
             else
-              evdw_t=evdw_t+evdwij
+              evdw_t=evdw_t+sss1*evdwij
             endif
             if (calc_grad) then
 C 
 C Calculate the components of the gradient in DC and X
 C
-            fac=-rrij*(e1+evdwij)
+            fac=-rrij*(e1+evdwij)*sss1
+     &          +evdwij*sssgrad1/sqrij/expon
             gg(1)=xj*fac
             gg(2)=yj*fac
             gg(3)=zj*fac
@@ -699,6 +733,7 @@ C
               enddo
             enddo
             endif
+#ifdef FOURBODY
 C
 C 12/1/95, revised on 5/20/97
 C
@@ -755,10 +790,13 @@ cd              write (iout,'(2i3,3f10.5)')
 cd   &           i,j,(gacont(kk,num_conti,i),kk=1,3)
               endif
             endif
+#endif
           enddo      ! j
         enddo        ! iint
+#ifdef FOURBODY
 C Change 12/1/95
         num_cont(i)=num_conti
+#endif
       enddo          ! i
       if (calc_grad) then
       do i=1,nct
@@ -830,6 +868,9 @@ C
             e_augm=augm(itypi,itypj)*fac_augm
             r_inv_ij=dsqrt(rrij)
             rij=1.0D0/r_inv_ij 
+            sss1=sscale(rij)
+            if (sss1.eq.0.0d0) cycle
+            sssgrad1=sscagrad(rij)
             r_shift_inv=1.0D0/(rij+r0(itypi,itypj)-sigma(itypi,itypj))
             fac=r_shift_inv**expon
             e1=fac*fac*aa
@@ -847,15 +888,16 @@ cd   &        bb(itypi,itypj),augm(itypi,itypj),epsi,sigm,
 cd   &        sigma(itypi,itypj),1.0D0/dsqrt(rrij),evdwij,
 cd   &        (c(k,i),k=1,3),(c(k,j),k=1,3)
             if (bb.gt.0.0d0) then
-              evdw=evdw+evdwij
+              evdw=evdw+evdwij*sss1
             else 
-              evdw_t=evdw_t+evdwij
+              evdw_t=evdw_t+evdwij*sss1
             endif
             if (calc_grad) then
 C 
 C Calculate the components of the gradient in DC and X
 C
-            fac=-2.0D0*rrij*e_augm-r_inv_ij*r_shift_inv*(e1+e1+e2)
+           fac=(-2.0D0*rrij*e_augm-r_inv_ij*r_shift_inv*(e1+e1+e2))*sss1
+     &          +evdwij*sssgrad1*r_inv_ij/expon
             gg(1)=xj*fac
             gg(2)=yj*fac
             gg(3)=zj*fac
@@ -1230,8 +1272,8 @@ C finding the closest
 c            write (iout,*) i,j,xj,yj,zj
             rrij=1.0D0/(xj*xj+yj*yj+zj*zj)
             rij=dsqrt(rrij)
-            sss=sscale((1.0d0/rij)/sigma(itypi,itypj))
-            sssgrad=sscagrad((1.0d0/rij)/sigma(itypi,itypj))
+            sss=sscale(1.0d0/rij))
+            sssgrad=sscagrad(1.0d0/rij)
             if (sss.le.0.0) cycle
 C Calculate angle-dependent terms of energy and contributions to their
 C derivatives.
@@ -1387,6 +1429,9 @@ c           alf12=0.0D0
             dzj=dc_norm(3,nres+j)
             rrij=1.0D0/(xj*xj+yj*yj+zj*zj)
             rij=dsqrt(rrij)
+            sss=sscale(1.0d0/rij)
+            if (sss.eq.0.0d0) cycle
+            sssgrad=sscagrad(1.0d0/rij)
 C Calculate angle-dependent terms of energy and contributions to their
 C derivatives.
             call sc_angular
@@ -1411,9 +1456,9 @@ c---------------------------------------------------------------
             e_augm=augm(itypi,itypj)*fac_augm
             evdwij=evdwij*eps2rt*eps3rt
             if (bb.gt.0.0d0) then
-              evdw=evdw+evdwij+e_augm
+              evdw=evdw+(evdwij+e_augm)*sss
             else
-              evdw_t=evdw_t+evdwij+e_augm
+              evdw_t=evdw_t+(evdwij+e_augm)*sss
             endif
             ij=icant(itypi,itypj)
             aux=eps1*eps2rt**2*eps3rt**2
@@ -1439,6 +1484,7 @@ C Calculate gradient components.
             fac=-expon*(e1+evdwij)*rij_shift
             sigder=fac*sigder
             fac=rij*fac-2*expon*rrij*e_augm
+            fac=fac+(evdwij+e_augm)*sssgrad/sss*rij
 C Calculate the radial part of the gradient
             gg(1)=xj*fac
             gg(2)=yj*fac
@@ -1717,6 +1763,7 @@ C--------------------------------------------------------------------------
       include 'COMMON.TORSION'
       include 'COMMON.VECTORS'
       include 'COMMON.FFIELD'
+      include 'COMMON.CORRMAT'
       double precision auxvec(2),auxmat(2,2)
 C
 C Compute the virtual-bond-torsional-angle dependent quantities needed
@@ -1943,6 +1990,7 @@ c     &    EE(1,2,iti),EE(2,2,i)
 c          write(iout,*) "Macierz EUG",
 c     &    eug(1,1,i-2),eug(1,2,i-2),eug(2,1,i-2),
 c     &    eug(2,2,i-2)
+#ifdef FOURBODY
           if (wcorr4.gt.0.0d0 .or. wcorr5.gt.0.0d0 .or. wcorr6.gt.0.0d0) 
      &    then
           call matmat2(CC(1,1,i-2),Ug(1,1,i-2),CUg(1,1,i-2))
@@ -1951,6 +1999,7 @@ c     &    eug(2,2,i-2)
           call matvec2(Ctilde(1,1,i-1),obrot(1,i-2),Ctobr(1,i-2))
           call matvec2(Dtilde(1,1,i-2),obrot2(1,i-2),Dtobr2(1,i-2))
           endif
+#endif
         else
           do k=1,2
             Ub2(k,i-2)=0.0d0
@@ -1992,6 +2041,7 @@ c        if (i.gt. iatel_s+1 .and. i.lt.iatel_e+4) then
 #endif
 cd        write (iout,*) 'mu1',mu1(:,i-2)
 cd        write (iout,*) 'mu2',mu2(:,i-2)
+#ifdef FOURBODY
         if (wcorr4.gt.0.0d0 .or. wcorr5.gt.0.0d0 .or.wcorr6.gt.0.0d0)
      &  then  
         if (calc_grad) then
@@ -2014,7 +2064,9 @@ C Vectors and matrices dependent on a single virtual-bond dihedral.
         call matmat2(EUgder(1,1,i-2),DD(1,1,i-1),EUgDder(1,1,i-2))
         endif
         endif
+#endif
       enddo
+#ifdef FOURBODY
 C Matrices dependent on two consecutive virtual-bond dihedrals.
 C The order of matrices is from left to right.
       if (wcorr4.gt.0.0d0 .or. wcorr5.gt.0.0d0 .or.wcorr6.gt.0.0d0)
@@ -2032,6 +2084,7 @@ C The order of matrices is from left to right.
         call transpose2(DtUg2der(1,1,i-1),auxmat(1,1))
         call matmat2(auxmat(1,1),EUg(1,1,i),Ug2DtEUgder(1,1,1,i))
         endif
+#endif
       enddo
       endif
       return
@@ -2058,7 +2111,11 @@ C
       include 'COMMON.CHAIN'
       include 'COMMON.DERIV'
       include 'COMMON.INTERACT'
+#ifdef FOURBODY
       include 'COMMON.CONTACTS'
+      include 'COMMON.CONTMAP'
+#endif
+      include 'COMMON.CORRMAT'
       include 'COMMON.TORSION'
       include 'COMMON.VECTORS'
       include 'COMMON.FFIELD'
@@ -2131,9 +2188,11 @@ cd      enddo
       eello_turn3=0.0d0
       eello_turn4=0.0d0
       ind=0
+#ifdef FOURBODY
       do i=1,nres
         num_cont_hb(i)=0
       enddo
+#endif
 cd      print '(a)','Enter EELEC'
 c      write (iout,*) 'iatel_s=',iatel_s,' iatel_e=',iatel_e
 c      call flush(iout)
@@ -2185,7 +2244,9 @@ c        end if
         num_conti=0
         call eelecij(i,i+2,ees,evdw1,eel_loc)
         if (wturn3.gt.0.0d0) call eturn3(i,eello_turn3)
+#ifdef FOURBODY
         num_cont_hb(i)=num_conti
+#endif
       enddo
       do i=iturn4_start,iturn4_end
         if (i.lt.1) cycle
@@ -2241,12 +2302,16 @@ c        endif
           zmedi=mod(zmedi,boxzsize)
           if (zmedi.lt.0) zmedi=zmedi+boxzsize
 
+#ifdef FOURBODY
         num_conti=num_cont_hb(i)
+#endif
 c        write(iout,*) "JESTEM W PETLI"
         call eelecij(i,i+3,ees,evdw1,eel_loc)
         if (wturn4.gt.0.0d0 .and. itype(i+2).ne.ntyp1) 
      &   call eturn4(i,eello_turn4)
+#ifdef FOURBODY
         num_cont_hb(i)=num_conti
+#endif
       enddo   ! i
 C Loop over all neighbouring boxes
 C      do xshift=-1,1
@@ -2313,7 +2378,9 @@ c        go to 166
 c        endif
 
 c        write (iout,*) 'i',i,' ielstart',ielstart(i),' ielend',ielend(i)
+#ifdef FOURBODY
         num_conti=num_cont_hb(i)
+#endif
 C I TU KURWA
         do j=ielstart(i),ielend(i)
 C          do j=16,17
@@ -2329,7 +2396,9 @@ c     & .or.itype(j-1).eq.ntyp1
      &) cycle
           call eelecij(i,j,ees,evdw1,eel_loc)
         enddo ! j
+#ifdef FOURBODY
         num_cont_hb(i)=num_conti
+#endif
       enddo   ! i
 C     enddo   ! zshift
 C      enddo   ! yshift
@@ -2360,7 +2429,11 @@ C-------------------------------------------------------------------------------
       include 'COMMON.CHAIN'
       include 'COMMON.DERIV'
       include 'COMMON.INTERACT'
+#ifdef FOURBODY
       include 'COMMON.CONTACTS'
+      include 'COMMON.CONTMAP'
+#endif
+      include 'COMMON.CORRMAT'
       include 'COMMON.TORSION'
       include 'COMMON.VECTORS'
       include 'COMMON.FFIELD'
@@ -2478,8 +2551,9 @@ C        yj=yj-ymedi
 C        zj=zj-zmedi
           rij=xj*xj+yj*yj+zj*zj
 
-            sss=sscale(sqrt(rij))
-            sssgrad=sscagrad(sqrt(rij))
+          sss=sscale(sqrt(rij))
+          if (sss.eq.0.0d0) return
+          sssgrad=sscagrad(sqrt(rij))
 c            write (iout,*) "ij",i,j," rij",sqrt(rij)," r_cut",r_cut,
 c     &       " rlamb",rlamb," sss",sss
 c            if (sss.gt.0.0d0) then  
@@ -2647,9 +2721,10 @@ cgrad              gelc(l,k)=gelc(l,k)+ggg(l)
 cgrad            enddo
 cgrad          enddo
           if (sss.gt.0.0) then
-          ggg(1)=facvdw*xj+sssgrad*rmij*evdwij*xj
-          ggg(2)=facvdw*yj+sssgrad*rmij*evdwij*yj
-          ggg(3)=facvdw*zj+sssgrad*rmij*evdwij*zj
+          facvdw=facvdw+sssgrad*rmij*evdwij
+          ggg(1)=facvdw*xj
+          ggg(2)=facvdw*yj
+          ggg(3)=facvdw*zj
           else
           ggg(1)=0.0
           ggg(2)=0.0
@@ -2676,10 +2751,11 @@ cgrad          enddo
           endif ! calc_grad
 #else
 C MARYSIA
-          facvdw=(ev1+evdwij)*sss
+          facvdw=(ev1+evdwij)
           facel=(el1+eesij)
           fac1=fac
-          fac=-3*rrmij*(facvdw+facvdw+facel)
+          fac=-3*rrmij*(facvdw+facvdw+facel)*sss
+     &       +(evdwij+eesij)*sssgrad*rrmij
           erij(1)=xj*rmij
           erij(2)=yj*rmij
           erij(3)=zj*rmij
@@ -3002,7 +3078,7 @@ c           if (eel_loc_ij.ne.0)
 c     &      write (iout,'(a4,2i4,8f9.5)')'chuj',
 c     &     i,j,a22,muij(1),a23,muij(2),a32,muij(3),a33,muij(4)
 
-          eel_loc=eel_loc+eel_loc_ij
+          eel_loc=eel_loc+eel_loc_ij*sss
 C Now derivative over eel_loc
           if (calc_grad) then
           if ((fac_shield(i).gt.0).and.(fac_shield(j).gt.0).and.
@@ -3059,7 +3135,7 @@ C Calculate patrial derivative for theta angle
      &     +a23*gmuij1(2)
      &     +a32*gmuij1(3)
      &     +a33*gmuij1(4))
-     &    *fac_shield(i)*fac_shield(j)
+     &    *fac_shield(i)*fac_shield(j)*sss
 c         write(iout,*) "derivative over thatai"
 c         write(iout,*) a22*gmuij1(1), a23*gmuij1(2) ,a32*gmuij1(3),
 c     &   a33*gmuij1(4) 
@@ -3075,7 +3151,7 @@ c     &   a33*gmuij2(4)
      &     +a33*gmuij2(4)
          gloc(nphi+i-1,icg)=gloc(nphi+i-1,icg)+
      &      geel_loc_ij*wel_loc
-     &    *fac_shield(i)*fac_shield(j)
+     &    *fac_shield(i)*fac_shield(j)*sss
 
 c  Derivative over j residue
          geel_loc_ji=a22*gmuji1(1)
@@ -3100,7 +3176,7 @@ c         write(iout,*) a22*gmuji2(1), a23*gmuji2(2) ,a32*gmuji2(3),
 c     &   a33*gmuji2(4)
          gloc(nphi+j-1,icg)=gloc(nphi+j-1,icg)+
      &      geel_loc_ji*wel_loc
-     &    *fac_shield(i)*fac_shield(j)
+     &    *fac_shield(i)*fac_shield(j)*sss
 #endif
 cd          write (iout,*) 'i',i,' j',j,' eel_loc_ij',eel_loc_ij
 
@@ -3116,10 +3192,14 @@ C Partial derivatives in virtual-bond dihedral angles gamma
      &           +a32*mu(2,i)*muder(1,j)+a33*mu(2,i)*muder(2,j))
      &    *fac_shield(i)*fac_shield(j)
 C Derivatives of eello in DC(i+1) thru DC(j-1) or DC(nres-2)
+          aux=eel_loc_ij/sss*sssgrad*rmij
+          ggg(1)=aux*xj
+          ggg(2)=aux*yj
+          ggg(3)=aux*zj
           do l=1,3
-            ggg(l)=(agg(l,1)*muij(1)+
+            ggg(l)=ggg(l)+(agg(l,1)*muij(1)+
      &          agg(l,2)*muij(2)+agg(l,3)*muij(3)+agg(l,4)*muij(4))
-     &    *fac_shield(i)*fac_shield(j)
+     &    *fac_shield(i)*fac_shield(j)*sss
             gel_loc_long(l,j)=gel_loc_long(l,j)+ggg(l)
             gel_loc_long(l,i)=gel_loc_long(l,i)-ggg(l)
 cgrad            ghalf=0.5d0*ggg(l)
@@ -3156,6 +3236,7 @@ C Remaining derivatives of eello
 
 C Change 12/26/95 to calculate four-body contributions to H-bonding energy
 c          if (j.gt.i+1 .and. num_conti.le.maxconts) then
+#ifdef FOURBODY
           if (wcorr+wcorr4+wcorr5+wcorr6.gt.0.0d0
      &       .and. num_conti.le.maxconts) then
 c            write (iout,*) i,j," entered corr"
@@ -3295,11 +3376,17 @@ cd              fprimcont=0.0D0
                   gggm(k)=ecosbm*dcosb(k)+ecosgm*dcosg(k)
                 enddo
                 gggp(1)=gggp(1)+ees0pijp*xj
+     &          +ees0p(num_conti,i)/sss*rmij*xj*sssgrad
                 gggp(2)=gggp(2)+ees0pijp*yj
+     &          +ees0p(num_conti,i)/sss*rmij*yj*sssgrad
                 gggp(3)=gggp(3)+ees0pijp*zj
+     &          +ees0p(num_conti,i)/sss*rmij*zj*sssgrad
                 gggm(1)=gggm(1)+ees0mijp*xj
+     &          +ees0m(num_conti,i)/sss*rmij*xj*sssgrad
                 gggm(2)=gggm(2)+ees0mijp*yj
+     &          +ees0m(num_conti,i)/sss*rmij*yj*sssgrad
                 gggm(3)=gggm(3)+ees0mijp*zj
+     &          +ees0m(num_conti,i)/sss*rmij*zj*sssgrad
 C Derivatives due to the contact function
                 gacont_hbr(1,num_conti,i)=fprimcont*xj
                 gacont_hbr(2,num_conti,i)=fprimcont*yj
@@ -3314,29 +3401,29 @@ cgrad                  ghalfm=0.5D0*gggm(k)
                   gacontp_hb1(k,num_conti,i)=!ghalfp
      &              +(ecosap*(dc_norm(k,j)-cosa*dc_norm(k,i))
      &              + ecosbp*(erij(k)-cosb*dc_norm(k,i)))*vbld_inv(i+1)
-     &          *fac_shield(i)*fac_shield(j)
+     &          *fac_shield(i)*fac_shield(j)*sss
 
                   gacontp_hb2(k,num_conti,i)=!ghalfp
      &              +(ecosap*(dc_norm(k,i)-cosa*dc_norm(k,j))
      &              + ecosgp*(erij(k)-cosg*dc_norm(k,j)))*vbld_inv(j+1)
-     &          *fac_shield(i)*fac_shield(j)
+     &          *fac_shield(i)*fac_shield(j)*sss
 
                   gacontp_hb3(k,num_conti,i)=gggp(k)
-     &          *fac_shield(i)*fac_shield(j)
+     &          *fac_shield(i)*fac_shield(j)*sss
 
                   gacontm_hb1(k,num_conti,i)=!ghalfm
      &              +(ecosam*(dc_norm(k,j)-cosa*dc_norm(k,i))
      &              + ecosbm*(erij(k)-cosb*dc_norm(k,i)))*vbld_inv(i+1)
-     &          *fac_shield(i)*fac_shield(j)
+     &          *fac_shield(i)*fac_shield(j)*sss
 
                   gacontm_hb2(k,num_conti,i)=!ghalfm
      &              +(ecosam*(dc_norm(k,i)-cosa*dc_norm(k,j))
      &              + ecosgm*(erij(k)-cosg*dc_norm(k,j)))*vbld_inv(j+1)
-     &          *fac_shield(i)*fac_shield(j)
+     &          *fac_shield(i)*fac_shield(j)*sss
 
                   gacontm_hb3(k,num_conti,i)=gggm(k)
      &          *fac_shield(i)*fac_shield(j)
-
+*sss
                 enddo
 C Diagnostics. Comment out or remove after debugging!
 cdiag           do k=1,3
@@ -3354,6 +3441,7 @@ cdiag           enddo
               endif  ! num_conti.le.maxconts
             endif  ! fcont.gt.0
           endif    ! j.gt.i+1
+#endif
           if (calc_grad) then
           if (wturn3.gt.0.0d0 .or. wturn4.gt.0.0d0) then
             do k=1,4
@@ -6270,6 +6358,7 @@ c        gsccor_loc(i-3)=gloci
       enddo
       return
       end
+#ifdef FOURBODY
 c------------------------------------------------------------------------------
       subroutine multibody(ecorr)
 C This subroutine calculates multi-body contributions to energy following
@@ -6282,6 +6371,8 @@ C contribution equal to sqrt(eps(i,j)*eps(i+1,j+1)) is added.
       include 'COMMON.DERIV'
       include 'COMMON.INTERACT'
       include 'COMMON.CONTACTS'
+      include 'COMMON.CONTMAT'
+      include 'COMMON.CORRMAT'
       double precision gx(3),gx1(3)
       logical lprn
 
@@ -6336,6 +6427,8 @@ c------------------------------------------------------------------------------
       include 'COMMON.DERIV'
       include 'COMMON.INTERACT'
       include 'COMMON.CONTACTS'
+      include 'COMMON.CONTMAT'
+      include 'COMMON.CORRMAT'
       double precision gx(3),gx1(3)
       logical lprn
       lprn=.false.
@@ -6377,6 +6470,8 @@ C This subroutine calculates multi-body contributions to hydrogen-bonding
       include 'COMMON.DERIV'
       include 'COMMON.INTERACT'
       include 'COMMON.CONTACTS'
+      include 'COMMON.CONTMAT'
+      include 'COMMON.CORRMAT'
       double precision gx(3),gx1(3)
       logical lprn,ldone
 
@@ -6449,6 +6544,8 @@ C This subroutine calculates multi-body contributions to hydrogen-bonding
       include 'COMMON.LOCAL'
       include 'COMMON.INTERACT'
       include 'COMMON.CONTACTS'
+      include 'COMMON.CONTMAT'
+      include 'COMMON.CORRMAT'
       include 'COMMON.CHAIN'
       include 'COMMON.CONTROL'
       include 'COMMON.SHIELD'
@@ -6605,6 +6702,8 @@ c------------------------------------------------------------------------------
       include 'COMMON.DERIV'
       include 'COMMON.INTERACT'
       include 'COMMON.CONTACTS'
+      include 'COMMON.CONTMAT'
+      include 'COMMON.CORRMAT'
       include 'COMMON.SHIELD'
       include 'COMMON.CONTROL'
       double precision gx(3),gx1(3)
@@ -6780,6 +6879,8 @@ C---------------------------------------------------------------------------
       include 'COMMON.DERIV'
       include 'COMMON.INTERACT'
       include 'COMMON.CONTACTS'
+      include 'COMMON.CONTMAT'
+      include 'COMMON.CORRMAT'
       include 'COMMON.TORSION'
       include 'COMMON.VAR'
       include 'COMMON.GEO'
@@ -6845,6 +6946,8 @@ C
       include 'COMMON.DERIV'
       include 'COMMON.INTERACT'
       include 'COMMON.CONTACTS'
+      include 'COMMON.CONTMAT'
+      include 'COMMON.CORRMAT'
       include 'COMMON.TORSION'
       include 'COMMON.VAR'
       include 'COMMON.GEO'
@@ -7223,6 +7326,8 @@ C---------------------------------------------------------------------------
       include 'COMMON.DERIV'
       include 'COMMON.INTERACT'
       include 'COMMON.CONTACTS'
+      include 'COMMON.CONTMAT'
+      include 'COMMON.CORRMAT'
       include 'COMMON.TORSION'
       include 'COMMON.VAR'
       include 'COMMON.GEO'
@@ -7337,6 +7442,8 @@ C---------------------------------------------------------------------------
       include 'COMMON.DERIV'
       include 'COMMON.INTERACT'
       include 'COMMON.CONTACTS'
+      include 'COMMON.CONTMAT'
+      include 'COMMON.CORRMAT'
       include 'COMMON.TORSION'
       include 'COMMON.VAR'
       include 'COMMON.GEO'
@@ -7753,6 +7860,8 @@ c--------------------------------------------------------------------------
       include 'COMMON.DERIV'
       include 'COMMON.INTERACT'
       include 'COMMON.CONTACTS'
+      include 'COMMON.CONTMAT'
+      include 'COMMON.CORRMAT'
       include 'COMMON.TORSION'
       include 'COMMON.VAR'
       include 'COMMON.GEO'
@@ -7895,6 +8004,8 @@ c--------------------------------------------------------------------------
       include 'COMMON.DERIV'
       include 'COMMON.INTERACT'
       include 'COMMON.CONTACTS'
+      include 'COMMON.CONTMAT'
+      include 'COMMON.CORRMAT'
       include 'COMMON.TORSION'
       include 'COMMON.VAR'
       include 'COMMON.GEO'
@@ -8001,6 +8112,8 @@ c----------------------------------------------------------------------------
       include 'COMMON.DERIV'
       include 'COMMON.INTERACT'
       include 'COMMON.CONTACTS'
+      include 'COMMON.CONTMAT'
+      include 'COMMON.CORRMAT'
       include 'COMMON.TORSION'
       include 'COMMON.VAR'
       include 'COMMON.GEO'
@@ -8188,6 +8301,8 @@ c----------------------------------------------------------------------------
       include 'COMMON.DERIV'
       include 'COMMON.INTERACT'
       include 'COMMON.CONTACTS'
+      include 'COMMON.CONTMAT'
+      include 'COMMON.CORRMAT'
       include 'COMMON.TORSION'
       include 'COMMON.VAR'
       include 'COMMON.GEO'
@@ -8305,6 +8420,8 @@ c----------------------------------------------------------------------------
       include 'COMMON.DERIV'
       include 'COMMON.INTERACT'
       include 'COMMON.CONTACTS'
+      include 'COMMON.CONTMAT'
+      include 'COMMON.CORRMAT'
       include 'COMMON.TORSION'
       include 'COMMON.VAR'
       include 'COMMON.GEO'
@@ -8551,6 +8668,8 @@ c----------------------------------------------------------------------------
       include 'COMMON.DERIV'
       include 'COMMON.INTERACT'
       include 'COMMON.CONTACTS'
+      include 'COMMON.CONTMAT'
+      include 'COMMON.CORRMAT'
       include 'COMMON.TORSION'
       include 'COMMON.VAR'
       include 'COMMON.GEO'
@@ -8871,7 +8990,7 @@ cd      write (2,*) 'ekont',ekont
 cd      write (2,*) 'eel_turn6',ekont*eel_turn6
       return
       end
-
+#endif
 crc-------------------------------------------------
 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
       subroutine Eliptransfer(eliptran)