X-Git-Url: http://mmka.chem.univ.gda.pl/gitweb/?a=blobdiff_plain;f=source%2Funres%2Fsrc-HCD-5D%2Fenergy_p_new_barrier.F;h=44023d05011ce004f21e1439b6f41d07cb0dd2e4;hb=946fab962642f90563ed751be8d845abc23f7457;hp=e9ec117a8e39dfc441d063bfaa32911a688a96fb;hpb=48ae9e01d2dd6571fa2cca6c704dc04f86e5fd7b;p=unres.git

diff --git a/source/unres/src-HCD-5D/energy_p_new_barrier.F b/source/unres/src-HCD-5D/energy_p_new_barrier.F
index e9ec117..44023d0 100644
--- a/source/unres/src-HCD-5D/energy_p_new_barrier.F
+++ b/source/unres/src-HCD-5D/energy_p_new_barrier.F
@@ -30,6 +30,7 @@ c      include 'COMMON.MD'
       include 'COMMON.SPLITELE'
       include 'COMMON.TORCNSTR'
       include 'COMMON.SAXS'
+      include 'COMMON.MD'
       double precision evdw,evdw1,evdw2,evdw2_14,ees,eel_loc,
      & eello_turn3,eello_turn4,edfadis,estr,ehpb,ebe,ethetacnstr,
      & escloc,etors,edihcnstr,etors_d,esccor,ecorr,ecorr5,ecorr6,eturn6,
@@ -66,7 +67,8 @@ C FG slaves as WEIGHTS array.
           weights_(17)=wbond
           weights_(18)=scal14
           weights_(21)=wsccor
-          weights_(22)=wtube
+          weights_(22)=wliptran
+          weights_(25)=wtube
           weights_(26)=wsaxs
           weights_(28)=wdfa_dist
           weights_(29)=wdfa_tor
@@ -98,7 +100,8 @@ C FG slaves receive the WEIGHTS array
           wbond=weights(17)
           scal14=weights(18)
           wsccor=weights(21)
-          wtube=weights(22)
+          wliptran=weights(22)
+          wtube=weights(25)
           wsaxs=weights(26)
           wdfa_dist=weights_(28)
           wdfa_tor=weights_(29)
@@ -109,12 +112,15 @@ C FG slaves receive the WEIGHTS array
         time_Bcastw=time_Bcastw+MPI_Wtime()-time00
 c        call chainbuild_cart
       endif
-#ifndef DFA
-      edfadis=0.0d0
-      edfator=0.0d0
-      edfanei=0.0d0
-      edfabet=0.0d0
-#endif
+      if (nfgtasks.gt.1) then
+        call MPI_Bcast(itime_mat,1,MPI_INT,king,FG_COMM,IERROR)
+      endif
+      if (mod(itime_mat,imatupdate).eq.0) then
+        call make_SCp_inter_list
+        call make_SCSC_inter_list
+        call make_pp_inter_list
+        call make_pp_vdw_inter_list
+      endif
 c      print *,'Processor',myrank,' calling etotal ipot=',ipot
 c      print *,'Processor',myrank,' nnt=',nnt,' nct=',nct
 #else
@@ -125,6 +131,13 @@ c      endif
 #ifdef TIMING
       time00=MPI_Wtime()
 #endif
+
+#ifndef DFA
+      edfadis=0.0d0
+      edfator=0.0d0
+      edfanei=0.0d0
+      edfabet=0.0d0
+#endif
 C 
 C Compute the side-chain and electrostatic interaction energy
 C
@@ -324,6 +337,7 @@ C
       else
         esccor=0.0d0
       endif
+#ifdef FOURBODY
 C      print *,"PRZED MULIt"
 c      print *,"Processor",myrank," computed Usccorr"
 C 
@@ -352,6 +366,7 @@ c         write (iout,*) "MULTIBODY_HB ecorr",ecorr,ecorr5,ecorr6,n_corr,
 c     &     n_corr1
 c         call flush(iout)
       endif
+#endif
 c      print *,"Processor",myrank," computed Ucorr"
 c      write (iout,*) "nsaxs",nsaxs," saxs_mode",saxs_mode
       if (nsaxs.gt.0 .and. saxs_mode.eq.0) then
@@ -385,6 +400,8 @@ C based on partition function
 C      print *,"przed lipidami"
       if (wliptran.gt.0) then
         call Eliptransfer(eliptran)
+      else
+        eliptran=0.0d0
       endif
 C      print *,"za lipidami"
       if (AFMlog.gt.0) then
@@ -1314,9 +1331,16 @@ C     Bartek
       write (iout,10) evdw,wsc,evdw2,wscp,ees,welec,evdw1,wvdwpp,
      &  estr,wbond,ebe,wang,
      &  escloc,wscloc,etors,wtor,etors_d,wtor_d,ehpb,wstrain,
+#ifdef FOURBODY
      &  ecorr,wcorr,
-     &  ecorr5,wcorr5,ecorr6,wcorr6,eel_loc,wel_loc,eello_turn3,wturn3,
-     &  eello_turn4,wturn4,eello_turn6,wturn6,esccor,wsccor,edihcnstr,
+     &  ecorr5,wcorr5,ecorr6,wcorr6,
+#endif
+     &  eel_loc,wel_loc,eello_turn3,wturn3,
+     &  eello_turn4,wturn4,
+#ifdef FOURBODY
+     &  eello_turn6,wturn6,
+#endif
+     &  esccor,wsccor,edihcnstr,
      &  ethetacnstr,ebr*nss,Uconst,wumb,eliptran,wliptran,Eafmforce,
      &  etube,wtube,esaxs,wsaxs,ehomology_constr,
      &  edfadis,wdfa_dist,edfator,wdfa_tor,edfanei,wdfa_nei,
@@ -1334,13 +1358,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)'/
@@ -1361,9 +1389,16 @@ C     Bartek
       write (iout,10) evdw,wsc,evdw2,wscp,ees,welec,
      &  estr,wbond,ebe,wang,
      &  escloc,wscloc,etors,wtor,etors_d,wtor_d,ehpb,wstrain,
+#ifdef FOURBODY
      &  ecorr,wcorr,
-     &  ecorr5,wcorr5,ecorr6,wcorr6,eel_loc,wel_loc,eello_turn3,wturn3,
-     &  eello_turn4,wturn4,eello_turn6,wturn6,esccor,wsccor,edihcnstr,
+     &  ecorr5,wcorr5,ecorr6,wcorr6,
+#endif
+     &  eel_loc,wel_loc,eello_turn3,wturn3,
+     &  eello_turn4,wturn4,
+#ifdef FOURBODY
+     &  eello_turn6,wturn6,
+#endif
+     &  esccor,wsccor,edihcnstr,
      &  ethetacnstr,ebr*nss,Uconst,wumb,eliptran,wliptran,Eafmforc,
      &  etube,wtube,esaxs,wsaxs,ehomology_constr,
      &  edfadis,wdfa_dist,edfator,wdfa_tor,edfanei,wdfa_nei,
@@ -1380,13 +1415,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)'/
@@ -1425,16 +1464,24 @@ C
       include 'COMMON.SBRIDGE'
       include 'COMMON.NAMES'
       include 'COMMON.IOUNITS'
+      include 'COMMON.SPLITELE'
+#ifdef FOURBODY
       include 'COMMON.CONTACTS'
+      include 'COMMON.CONTMAT'
+#endif
       double precision gg(3)
       double precision evdw,evdwij
-      integer i,j,k,itypi,itypj,itypi1,num_conti,iint
+      integer i,j,k,itypi,itypj,itypi1,num_conti,iint,icont
       double precision xi,yi,zi,xj,yj,zj,rij,eps0ij,fac,e1,e2,rrij,
-     & sigij,r0ij,rcut
+     & sigij,r0ij,rcut,sqrij,sss1,sssgrad1
       double precision fcont,fprimcont
+      double precision sscale,sscagrad
 c      write(iout,*)'Entering ELJ nnt=',nnt,' nct=',nct,' expon=',expon
       evdw=0.0D0
-      do i=iatsc_s,iatsc_e
+c      do i=iatsc_s,iatsc_e
+      do icont=g_listscsc_start,g_listscsc_end
+        i=newcontlisti(icont)
+        j=newcontlistj(icont)
         itypi=iabs(itype(i))
         if (itypi.eq.ntyp1) cycle
         itypi1=iabs(itype(i+1))
@@ -1446,10 +1493,10 @@ C Change 12/1/95
 C
 C Calculate SC interaction energy.
 C
-        do iint=1,nint_gr(i)
+c        do iint=1,nint_gr(i)
 cd        write (iout,*) 'i=',i,' iint=',iint,' istart=',istart(i,iint),
 cd   &                  'iend=',iend(i,iint)
-          do j=istart(i,iint),iend(i,iint)
+c          do j=istart(i,iint),iend(i,iint)
             itypj=iabs(itype(j)) 
             if (itypj.eq.ntyp1) cycle
             xj=c(1,nres+j)-xi
@@ -1458,6 +1505,11 @@ 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,r_cut_int)
+            if (sss1.eq.0.0d0) cycle
+            sssgrad1=sscagrad(sqrij,r_cut_int)
+            
 c           write (iout,*)'i=',i,' j=',j,' itypi=',itypi,' itypj=',itypj
             eps0ij=eps(itypi,itypj)
             fac=rrij**expon2
@@ -1471,11 +1523,12 @@ cd          write (iout,'(2(a3,i3,2x),6(1pd12.4)/2(3(1pd12.4),5x)/)')
 cd   &        restyp(itypi),i,restyp(itypj),j,a(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)
-            evdw=evdw+evdwij
+            evdw=evdw+sss1*evdwij
 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
@@ -1491,6 +1544,7 @@ cgrad                gvdwc(l,k)=gvdwc(l,k)+gg(l)
 cgrad              enddo
 cgrad            enddo
 C
+#ifdef FOURBODY
 C 12/1/95, revised on 5/20/97
 C
 C Calculate the contact function. The ith column of the array JCONT will 
@@ -1546,10 +1600,13 @@ cd              write (iout,'(2i3,3f10.5)')
 cd   &           i,j,(gacont(kk,num_conti,i),kk=1,3)
               endif
             endif
-          enddo      ! j
-        enddo        ! iint
+#endif
+c          enddo      ! j
+c        enddo        ! iint
 C Change 12/1/95
+#ifdef FOURBODY
         num_cont(i)=num_conti
+#endif
       enddo          ! i
       do i=1,nct
         do j=1,3
@@ -1584,15 +1641,20 @@ C
       include 'COMMON.INTERACT'
       include 'COMMON.IOUNITS'
       include 'COMMON.NAMES'
+      include 'COMMON.SPLITELE'
       double precision gg(3)
       double precision evdw,evdwij
-      integer i,j,k,itypi,itypj,itypi1,iint
+      integer i,j,k,itypi,itypj,itypi1,iint,icont
       double precision xi,yi,zi,xj,yj,zj,rij,eps0ij,fac,e1,e2,rrij,
-     & fac_augm,e_augm,r_inv_ij,r_shift_inv
+     & fac_augm,e_augm,r_inv_ij,r_shift_inv,sss1,sssgrad1
       logical scheck
+      double precision sscale,sscagrad
 c     print *,'Entering ELJK nnt=',nnt,' nct=',nct,' expon=',expon
       evdw=0.0D0
-      do i=iatsc_s,iatsc_e
+c      do i=iatsc_s,iatsc_e
+      do icont=g_listscsc_start,g_listscsc_end
+        i=newcontlisti(icont)
+        j=newcontlistj(icont)
         itypi=iabs(itype(i))
         if (itypi.eq.ntyp1) cycle
         itypi1=iabs(itype(i+1))
@@ -1602,8 +1664,8 @@ c     print *,'Entering ELJK nnt=',nnt,' nct=',nct,' expon=',expon
 C
 C Calculate SC interaction energy.
 C
-        do iint=1,nint_gr(i)
-          do j=istart(i,iint),iend(i,iint)
+c        do iint=1,nint_gr(i)
+c          do j=istart(i,iint),iend(i,iint)
             itypj=iabs(itype(j))
             if (itypj.eq.ntyp1) cycle
             xj=c(1,nres+j)-xi
@@ -1614,6 +1676,9 @@ C
             e_augm=augm(itypi,itypj)*fac_augm
             r_inv_ij=dsqrt(rrij)
             rij=1.0D0/r_inv_ij 
+            sss1=sscale(rij,r_cut_int)
+            if (sss1.eq.0.0d0) cycle
+            sssgrad1=sscagrad(rij,r_cut_int)
             r_shift_inv=1.0D0/(rij+r0(itypi,itypj)-sigma(itypi,itypj))
             fac=r_shift_inv**expon
 C have you changed here?
@@ -1627,11 +1692,12 @@ cd   &        restyp(itypi),i,restyp(itypj),j,aa(itypi,itypj),
 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)
-            evdw=evdw+evdwij
+            evdw=evdw+evdwij*sss1
 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)
+     &          +evdwij*sssgrad1*r_inv_ij/expon
             gg(1)=xj*fac
             gg(2)=yj*fac
             gg(3)=zj*fac
@@ -1646,8 +1712,8 @@ cgrad              do l=1,3
 cgrad                gvdwc(l,k)=gvdwc(l,k)+gg(l)
 cgrad              enddo
 cgrad            enddo
-          enddo      ! j
-        enddo        ! iint
+c          enddo      ! j
+c        enddo        ! iint
       enddo          ! i
       do i=1,nct
         do j=1,3
@@ -1674,11 +1740,14 @@ C
       include 'COMMON.INTERACT'
       include 'COMMON.IOUNITS'
       include 'COMMON.CALC'
+      include 'COMMON.SPLITELE'
       integer icall
       common /srutu/ icall
       double precision evdw
-      integer itypi,itypj,itypi1,iint,ind
-      double precision eps0ij,epsi,sigm,fac,e1,e2,rrij,xi,yi,zi
+      integer itypi,itypj,itypi1,iint,ind,icont
+      double precision eps0ij,epsi,sigm,fac,e1,e2,rrij,xi,yi,zi,
+     & sss1,sssgrad1
+      double precision sscale,sscagrad
 c     double precision rrsave(maxdim)
       logical lprn
       evdw=0.0D0
@@ -1690,7 +1759,10 @@ c     else
         lprn=.false.
 c     endif
       ind=0
-      do i=iatsc_s,iatsc_e
+c      do i=iatsc_s,iatsc_e 
+      do icont=g_listscsc_start,g_listscsc_end
+        i=newcontlisti(icont)
+        j=newcontlistj(icont)
         itypi=iabs(itype(i))
         if (itypi.eq.ntyp1) cycle
         itypi1=iabs(itype(i+1))
@@ -1705,8 +1777,8 @@ c        dsci_inv=dsc_inv(itypi)
 C
 C Calculate SC interaction energy.
 C
-        do iint=1,nint_gr(i)
-          do j=istart(i,iint),iend(i,iint)
+c        do iint=1,nint_gr(i)
+c          do j=istart(i,iint),iend(i,iint)
             ind=ind+1
             itypj=iabs(itype(j))
             if (itypj.eq.ntyp1) cycle
@@ -1744,6 +1816,9 @@ cd          else
 cd            rrij=rrsave(ind)
 cd          endif
             rij=dsqrt(rrij)
+            sss1=sscale(1.0d0/rij,r_cut_int)
+            if (sss1.eq.0.0d0) cycle
+            sssgrad1=sscagrad(1.0d0/rij,r_cut_int)
 C Calculate the angle-dependent terms of energy & contributions to derivatives.
             call sc_angular
 C Calculate whole angle-dependent part of epsilon and contributions
@@ -1756,7 +1831,7 @@ C have you changed here?
             eps2der=evdwij*eps3rt
             eps3der=evdwij*eps2rt
             evdwij=evdwij*eps2rt*eps3rt
-            evdw=evdw+evdwij
+            evdw=evdw+sss1*evdwij
             if (lprn) then
             sigm=dabs(aa/bb)**(1.0D0/6.0D0)
             epsi=bb**2/aa
@@ -1772,6 +1847,7 @@ C Calculate gradient components.
             fac=-expon*(e1+evdwij)
             sigder=fac/sigsq
             fac=rrij*fac
+     &          +evdwij*sssgrad1/sss1*rij
 C Calculate radial part of the gradient
             gg(1)=xj*fac
             gg(2)=yj*fac
@@ -1779,8 +1855,8 @@ C Calculate radial part of the gradient
 C Calculate the angular part of the gradient and sum add the contributions
 C to the appropriate components of the Cartesian gradient.
             call sc_grad
-          enddo      ! j
-        enddo        ! iint
+!          enddo      ! j
+!        enddo        ! iint
       enddo          ! i
 c     stop
       return
@@ -1808,7 +1884,7 @@ C
       logical lprn
       integer xshift,yshift,zshift,subchap
       double precision evdw
-      integer itypi,itypj,itypi1,iint,ind
+      integer itypi,itypj,itypi1,iint,ind,icont
       double precision eps0ij,epsi,sigm,fac,e1,e2,rrij,xi,yi,zi
       double precision fracinbuf,sslipi,evdwij_przed_tri,sig0ij,
      & sslipj,ssgradlipj,ssgradlipi,dist_init,xj_safe,yj_safe,zj_safe,
@@ -1826,7 +1902,10 @@ C we have the original box)
 C      do xshift=-1,1
 C      do yshift=-1,1
 C      do zshift=-1,1
-      do i=iatsc_s,iatsc_e
+c      do i=iatsc_s,iatsc_e
+      do icont=g_listscsc_start,g_listscsc_end
+        i=newcontlisti(icont)
+        j=newcontlistj(icont)
         itypi=iabs(itype(i))
         if (itypi.eq.ntyp1) cycle
         itypi1=iabs(itype(i+1))
@@ -1907,8 +1986,8 @@ c        write (iout,*) "dcnori",dxi*dxi+dyi*dyi+dzi*dzi
 C
 C Calculate SC interaction energy.
 C
-        do iint=1,nint_gr(i)
-          do j=istart(i,iint),iend(i,iint)
+c        do iint=1,nint_gr(i)
+c          do j=istart(i,iint),iend(i,iint)
             IF (dyn_ss_mask(i).and.dyn_ss_mask(j)) THEN
 
 c              write(iout,*) "PRZED ZWYKLE", evdwij
@@ -2076,12 +2155,11 @@ c            write (iout,*) "j",j," dc_norm",
 c     &       dc_norm(1,nres+j),dc_norm(2,nres+j),dc_norm(3,nres+j)
             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,r_cut_int)
 c            write (iout,'(a7,4f8.3)') 
 c    &      "ssscale",sss,((1.0d0/rij)/sigma(itypi,itypj)),r_cut,rlamb
-            if (sss.gt.0.0d0) then
+            if (sss.eq.0.0d0) cycle
+            sssgrad=sscagrad(1.0d0/rij,r_cut_int)
 C Calculate angle-dependent terms of energy and contributions to their
 C derivatives.
             call sc_angular
@@ -2128,8 +2206,8 @@ c     &        " eps3rt",eps3rt," eps1",eps1," e1",e1," e2",e2
      &        evdwij
             endif
 
-            if (energy_dec) write (iout,'(a6,2i5,0pf7.3)') 
-     &                        'evdw',i,j,evdwij
+            if (energy_dec) write (iout,'(a,2i5,3f10.5)') 
+     &                    'r sss evdw',i,j,rij,sss,evdwij
 
 C Calculate gradient components.
             e1=e1*eps1*eps2rt**2*eps3rt**2
@@ -2138,13 +2216,13 @@ C Calculate gradient components.
             fac=rij*fac
 c            print '(2i4,6f8.4)',i,j,sss,sssgrad*
 c     &      evdwij,fac,sigma(itypi,itypj),expon
-            fac=fac+evdwij/sss*sssgrad/sigma(itypi,itypj)*rij
+            fac=fac+evdwij*sssgrad/sss*rij
 c            fac=0.0d0
 C Calculate the radial part of the gradient
             gg_lipi(3)=eps1*(eps2rt*eps2rt)
-     &*(eps3rt*eps3rt)*sss/2.0d0*(faclip*faclip*
-     & (aa_lip(itypi,itypj)-aa_aq(itypi,itypj))
-     &+faclip*(bb_lip(itypi,itypj)-bb_aq(itypi,itypj)))
+     &       *(eps3rt*eps3rt)*sss/2.0d0*(faclip*faclip*
+     &        (aa_lip(itypi,itypj)-aa_aq(itypi,itypj))
+     &       +faclip*(bb_lip(itypi,itypj)-bb_aq(itypi,itypj)))
             gg_lipj(3)=ssgradlipj*gg_lipi(3)
             gg_lipi(3)=gg_lipi(3)*ssgradlipi
 C            gg_lipi(3)=0.0d0
@@ -2153,11 +2231,11 @@ C            gg_lipj(3)=0.0d0
             gg(2)=yj*fac
             gg(3)=zj*fac
 C Calculate angular part of the gradient.
+c            call sc_grad_scale(sss)
             call sc_grad
-            endif
             ENDIF    ! dyn_ss            
-          enddo      ! j
-        enddo        ! iint
+c          enddo      ! j
+c        enddo        ! iint
       enddo          ! i
 C      enddo          ! zshift
 C      enddo          ! yshift
@@ -2183,17 +2261,18 @@ C
       include 'COMMON.INTERACT'
       include 'COMMON.IOUNITS'
       include 'COMMON.CALC'
+      include 'COMMON.SPLITELE'
       integer xshift,yshift,zshift,subchap
       integer icall
       common /srutu/ icall
       logical lprn
       double precision evdw
-      integer itypi,itypj,itypi1,iint,ind
+      integer itypi,itypj,itypi1,iint,ind,icont
       double precision eps0ij,epsi,sigm,fac,e1,e2,rrij,r0ij,
      & xi,yi,zi,fac_augm,e_augm
       double precision fracinbuf,sslipi,evdwij_przed_tri,sig0ij,
      & sslipj,ssgradlipj,ssgradlipi,dist_init,xj_safe,yj_safe,zj_safe,
-     & xj_temp,yj_temp,zj_temp,dist_temp,sig,rij_shift,faclip
+     & xj_temp,yj_temp,zj_temp,dist_temp,sig,rij_shift,faclip,sssgrad1
       double precision dist,sscale,sscagrad,sscagradlip,sscalelip
       evdw=0.0D0
 c     print *,'Entering EGB nnt=',nnt,' nct=',nct,' expon=',expon
@@ -2201,7 +2280,10 @@ c     print *,'Entering EGB nnt=',nnt,' nct=',nct,' expon=',expon
       lprn=.false.
 c     if (icall.eq.0) lprn=.true.
       ind=0
-      do i=iatsc_s,iatsc_e
+c      do i=iatsc_s,iatsc_e
+      do icont=g_listscsc_start,g_listscsc_end
+        i=newcontlisti(icont)
+        j=newcontlistj(icont)
         itypi=iabs(itype(i))
         if (itypi.eq.ntyp1) cycle
         itypi1=iabs(itype(i+1))
@@ -2251,8 +2333,8 @@ c        dsci_inv=dsc_inv(itypi)
 C
 C Calculate SC interaction energy.
 C
-        do iint=1,nint_gr(i)
-          do j=istart(i,iint),iend(i,iint)
+c        do iint=1,nint_gr(i)
+c          do j=istart(i,iint),iend(i,iint)
             ind=ind+1
             itypj=iabs(itype(j))
             if (itypj.eq.ntyp1) cycle
@@ -2353,6 +2435,9 @@ C      write(iout,*) "tu,", i,j,aa,bb,aa_lip(itypi,itypj),sslipi,sslipj
             dzj=dc_norm(3,nres+j)
             rrij=1.0D0/(xj*xj+yj*yj+zj*zj)
             rij=dsqrt(rrij)
+            sss=sscale(1.0d0/rij,r_cut_int)
+            if (sss.eq.0.0d0) cycle
+            sssgrad=sscagrad(1.0d0/rij,r_cut_int)
 C Calculate angle-dependent terms of energy and contributions to their
 C derivatives.
             call sc_angular
@@ -2393,15 +2478,16 @@ C Calculate gradient components.
             fac=-expon*(e1+evdwij)*rij_shift
             sigder=fac*sigder
             fac=rij*fac-2*expon*rrij*e_augm
-            fac=fac+evdwij/sss*sssgrad/sigma(itypi,itypj)*rij
+            fac=fac+(evdwij+e_augm)*sssgrad/sss*rij
 C Calculate the radial part of the gradient
             gg(1)=xj*fac
             gg(2)=yj*fac
             gg(3)=zj*fac
 C Calculate angular part of the gradient.
+c            call sc_grad_scale(sss)
             call sc_grad
-          enddo      ! j
-        enddo        ! iint
+c          enddo      ! j
+c        enddo        ! iint
       enddo          ! i
       end
 C-----------------------------------------------------------------------------
@@ -2548,11 +2634,14 @@ C
       include 'COMMON.SBRIDGE'
       include 'COMMON.NAMES'
       include 'COMMON.IOUNITS'
-      include 'COMMON.CONTACTS'
+c      include 'COMMON.CONTACTS'
       dimension gg(3)
 cd    print *,'Entering Esoft_sphere nnt=',nnt,' nct=',nct
       evdw=0.0D0
-      do i=iatsc_s,iatsc_e
+c      do i=iatsc_s,iatsc_e
+      do icont=g_listscsc_start,g_listscsc_end
+        i=newcontlisti(icont)
+        j=newcontlistj(icont)
         itypi=iabs(itype(i))
         if (itypi.eq.ntyp1) cycle
         itypi1=iabs(itype(i+1))
@@ -2562,10 +2651,10 @@ cd    print *,'Entering Esoft_sphere nnt=',nnt,' nct=',nct
 C
 C Calculate SC interaction energy.
 C
-        do iint=1,nint_gr(i)
+c        do iint=1,nint_gr(i)
 cd        write (iout,*) 'i=',i,' iint=',iint,' istart=',istart(i,iint),
 cd   &                  'iend=',iend(i,iint)
-          do j=istart(i,iint),iend(i,iint)
+c          do j=istart(i,iint),iend(i,iint)
             itypj=iabs(itype(j))
             if (itypj.eq.ntyp1) cycle
             xj=c(1,nres+j)-xi
@@ -2601,8 +2690,8 @@ cgrad              do l=1,3
 cgrad                gvdwc(l,k)=gvdwc(l,k)+gg(l)
 cgrad              enddo
 cgrad            enddo
-          enddo ! j
-        enddo ! iint
+c          enddo ! j
+c        enddo ! iint
       enddo ! i
       return
       end
@@ -2622,7 +2711,7 @@ C
       include 'COMMON.CHAIN'
       include 'COMMON.DERIV'
       include 'COMMON.INTERACT'
-      include 'COMMON.CONTACTS'
+c      include 'COMMON.CONTACTS'
       include 'COMMON.TORSION'
       include 'COMMON.VECTORS'
       include 'COMMON.FFIELD'
@@ -2703,8 +2792,8 @@ c        write (iout,*) 'i',i,' ielstart',ielstart(i),' ielend',ielend(i)
           zj=zj_safe-zmedi
        endif
           rij=xj*xj+yj*yj+zj*zj
-            sss=sscale(sqrt(rij))
-            sssgrad=sscagrad(sqrt(rij))
+            sss=sscale(sqrt(rij),r_cut_int)
+            sssgrad=sscagrad(sqrt(rij),r_cut_int)
           if (rij.lt.r0ijsq) then
             evdw1ij=0.25d0*(rij-r0ijsq)**2
             fac=rij-r0ijsq
@@ -2946,7 +3035,7 @@ C--------------------------------------------------------------------------
       include 'COMMON.CHAIN'
       include 'COMMON.DERIV'
       include 'COMMON.INTERACT'
-      include 'COMMON.CONTACTS'
+      include 'COMMON.CORRMAT'
       include 'COMMON.TORSION'
       include 'COMMON.VECTORS'
       include 'COMMON.FFIELD'
@@ -3197,6 +3286,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))
@@ -3205,6 +3295,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
@@ -3249,6 +3340,7 @@ c          mu(k,i-2)=Ub2(k,i-2)
 cd        write (iout,*) 'mu1',mu1(:,i-2)
 cd        write (iout,*) 'mu2',mu2(:,i-2)
 cd        write (iout,*) 'mu',i-2,mu(:,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-1),Ugder(1,1,i-2),CUgder(1,1,i-2))
@@ -3267,7 +3359,9 @@ C Vectors and matrices dependent on a single virtual-bond dihedral.
         call matmat2(EUg(1,1,i-2),DD(1,1,i-1),EUgD(1,1,i-2))
         call matmat2(EUgder(1,1,i-2),DD(1,1,i-1),EUgDder(1,1,i-2))
         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)
@@ -3284,6 +3378,7 @@ c      do i=max0(ivec_start,2),ivec_end
         call matmat2(auxmat(1,1),EUg(1,1,i),Ug2DtEUgder(1,1,1,i))
       enddo
       endif
+#endif
 #if defined(MPI) && defined(PARMAT)
 #ifdef DEBUG
 c      if (fg_rank.eq.0) then
@@ -3352,6 +3447,7 @@ c     &   " ivec_count",(ivec_count(i),i=0,nfgtasks-1)
         call MPI_Allgatherv(sintab2(ivec_start),ivec_count(fg_rank1),
      &   MPI_DOUBLE_PRECISION,sintab2(1),ivec_count(0),ivec_displ(0),
      &   MPI_DOUBLE_PRECISION,FG_COMM1,IERR)
+#ifdef FOURBODY
         if (wcorr4.gt.0.0d0 .or. wcorr5.gt.0.0d0 .or. wcorr6.gt.0.0d0)
      &  then
         call MPI_Allgatherv(Ctobr(1,ivec_start),ivec_count(fg_rank1),
@@ -3427,6 +3523,7 @@ c     &   " ivec_count",(ivec_count(i),i=0,nfgtasks-1)
      &   MPI_MAT2,Ug2DtEUgder(1,1,1,1),ivec_count(0),ivec_displ(0),
      &   MPI_MAT2,FG_COMM1,IERR)
         endif
+#endif
 #else
 c Passes matrix info through the ring
       isend=fg_rank1
@@ -3471,6 +3568,7 @@ c        call flush(iout)
      &   iprev,6600+irecv,FG_COMM,status,IERR)
 c        write (iout,*) "Gather PRECOMP12"
 c        call flush(iout)
+#ifdef FOURBODY
         if (wcorr4.gt.0.0d0 .or. wcorr5.gt.0.0d0 .or. wcorr6.gt.0.0d0) 
      &  then
         call MPI_SENDRECV(ug2db1t(1,ivec_displ(isend)+1),1,
@@ -3490,6 +3588,7 @@ c        call flush(iout)
      &   Ug2DtEUgder(1,1,1,ivec_displ(irecv)+1),1,
      &   MPI_PRECOMP23(lenrecv),
      &   iprev,9900+irecv,FG_COMM,status,IERR)
+#endif
 c        write (iout,*) "Gather PRECOMP23"
 c        call flush(iout)
         endif
@@ -3565,7 +3664,11 @@ C
       include 'COMMON.CHAIN'
       include 'COMMON.DERIV'
       include 'COMMON.INTERACT'
+#ifdef FOURBODY
       include 'COMMON.CONTACTS'
+      include 'COMMON.CONTMAT'
+#endif
+      include 'COMMON.CORRMAT'
       include 'COMMON.TORSION'
       include 'COMMON.VECTORS'
       include 'COMMON.FFIELD'
@@ -3638,9 +3741,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'
 cd      write (iout,*) 'iatel_s=',iatel_s,' iatel_e=',iatel_e
       do i=1,nres
@@ -3690,7 +3795,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
@@ -3746,12 +3853,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
@@ -3761,7 +3872,10 @@ c
 c Loop over all pairs of interacting peptide groups except i,i+2 and i,i+3
 c
 CTU KURWA
-      do i=iatel_s,iatel_e
+c      do i=iatel_s,iatel_e
+      do icont=g_listpp_start,g_listpp_end
+        i=newcontlistppi(icont)
+        j=newcontlistppj(icont)
 C        do i=75,75
 c        if (i.le.1) cycle
         if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1
@@ -3818,9 +3932,11 @@ 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=ielstart(i),ielend(i)
 C          do j=16,17
 C          write (iout,*) i,j
 C         if (j.le.1) cycle
@@ -3833,8 +3949,10 @@ c     & .or.itype(j+2).eq.ntyp1
 c     & .or.itype(j-1).eq.ntyp1
      &) cycle
           call eelecij(i,j,ees,evdw1,eel_loc)
-        enddo ! j
+c        enddo ! j
+#ifdef FOURBODY
         num_cont_hb(i)=num_conti
+#endif
       enddo   ! i
 C     enddo   ! zshift
 C      enddo   ! yshift
@@ -3852,7 +3970,7 @@ cd      print *,"Processor",fg_rank," t_eelecij",t_eelecij
       end
 C-------------------------------------------------------------------------------
       subroutine eelecij(i,j,ees,evdw1,eel_loc)
-      implicit real*8 (a-h,o-z)
+      implicit none
       include 'DIMENSIONS'
 #ifdef MPI
       include "mpif.h"
@@ -3865,21 +3983,44 @@ C-------------------------------------------------------------------------------
       include 'COMMON.CHAIN'
       include 'COMMON.DERIV'
       include 'COMMON.INTERACT'
+#ifdef FOURBODY
       include 'COMMON.CONTACTS'
+      include 'COMMON.CONTMAT'
+#endif
+      include 'COMMON.CORRMAT'
       include 'COMMON.TORSION'
       include 'COMMON.VECTORS'
       include 'COMMON.FFIELD'
       include 'COMMON.TIME1'
       include 'COMMON.SPLITELE'
       include 'COMMON.SHIELD'
-      dimension ggg(3),gggp(3),gggm(3),erij(3),dcosb(3),dcosg(3),
+      double precision ggg(3),gggp(3),gggm(3),erij(3),dcosb(3),dcosg(3),
      &          erder(3,3),uryg(3,3),urzg(3,3),vryg(3,3),vrzg(3,3)
       double precision acipa(2,2),agg(3,4),aggi(3,4),aggi1(3,4),
      &    aggj(3,4),aggj1(3,4),a_temp(2,2),muij(4),gmuij1(4),gmuji1(4),
      &    gmuij2(4),gmuji2(4)
+      double precision dxi,dyi,dzi
+      double precision dx_normi,dy_normi,dz_normi,aux
+      integer j1,j2,lll,num_conti
       common /locel/ a_temp,agg,aggi,aggi1,aggj,aggj1,a22,a23,a32,a33,
      &    dxi,dyi,dzi,dx_normi,dy_normi,dz_normi,xmedi,ymedi,zmedi,
      &    num_conti,j1,j2
+      integer k,i,j,iteli,itelj,kkk,l,kkll,m,isubchap,ilist,iresshield
+      double precision ael6i,rrmij,rmij,r0ij,fcont,fprimcont,ees0tmp
+      double precision ees,evdw1,eel_loc,aaa,bbb,ael3i
+      double precision dxj,dyj,dzj,dx_normj,dy_normj,dz_normj,xj,yj,zj,
+     &  rij,r3ij,r6ij,cosa,cosb,cosg,fac,ev1,ev2,fac3,fac4,
+     &  evdwij,el1,el2,eesij,ees0ij,facvdw,facel,fac1,ecosa,
+     &  ecosb,ecosg,ury,urz,vry,vrz,facr,a22der,a23der,a32der,
+     &  a33der,eel_loc_ij,cosa4,wij,cosbg1,cosbg2,ees0pij,
+     &  ees0pij1,ees0mij,ees0mij1,fac3p,ees0mijp,ees0pijp,
+     &  ecosa1,ecosb1,ecosg1,ecosa2,ecosb2,ecosg2,ecosap,ecosbp,
+     &  ecosgp,ecosam,ecosbm,ecosgm,ghalf,rlocshield
+      double precision a22,a23,a32,a33,geel_loc_ij,geel_loc_ji
+      double precision dist_init,xj_safe,yj_safe,zj_safe,
+     &  xj_temp,yj_temp,zj_temp,dist_temp,xmedi,ymedi,zmedi
+      double precision sscale,sscagrad,scalar
+
 c 4/26/02 - AL scaling factor for 1,4 repulsive VDW interactions
 #ifdef MOMENT
       double precision scal_el /1.0d0/
@@ -3983,8 +4124,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(dsqrt(rij),r_cut_int)
+          if (sss.eq.0.0d0) return
+          sssgrad=sscagrad(dsqrt(rij),r_cut_int)
 c            if (sss.gt.0.0d0) then  
           rrmij=1.0D0/rij
           rij=dsqrt(rij)
@@ -4019,7 +4161,7 @@ C          fac_shield(j)=0.6
           fac_shield(i)=1.0
           fac_shield(j)=1.0
           eesij=(el1+el2)
-          ees=ees+eesij
+          ees=ees+eesij*sss
           endif
           evdw1=evdw1+evdwij*sss
 cd          write(iout,'(2(2i3,2x),7(1pd12.4)/2(3(1pd12.4),5x)/)')
@@ -4028,11 +4170,10 @@ cd     &      1.0D0/dsqrt(rrmij),evdwij,eesij,
 cd     &      xmedi,ymedi,zmedi,xj,yj,zj
 
           if (energy_dec) then 
-              write (iout,'(a6,2i5,0pf7.3,2i5,3e11.3)') 
-     &'evdw1',i,j,evdwij
-     &,iteli,itelj,aaa,evdw1,sss
-              write (iout,'(a6,2i5,0pf7.3,2f8.3)') 'ees',i,j,eesij,
-     &fac_shield(i),fac_shield(j)
+            write (iout,'(a6,2i5,0pf7.3,2i5,e11.3,3f10.5)') 
+     &        'evdw1',i,j,evdwij,iteli,itelj,aaa,evdw1,sss,rij
+            write (iout,'(a6,2i5,0pf7.3,2f8.3)') 'ees',i,j,eesij,
+     &        fac_shield(i),fac_shield(j)
           endif
 
 C
@@ -4049,9 +4190,10 @@ C
 *
 * Radial derivatives. First process both termini of the fragment (i,j)
 *
-          ggg(1)=facel*xj
-          ggg(2)=facel*yj
-          ggg(3)=facel*zj
+          aux=facel*sss+rmij*sssgrad*eesij
+          ggg(1)=aux*xj
+          ggg(2)=aux*yj
+          ggg(3)=aux*zj
           if ((fac_shield(i).gt.0).and.(fac_shield(j).gt.0).and.
      &  (shield_mode.gt.0)) then
 C          print *,i,j     
@@ -4085,10 +4227,10 @@ C              endif
            iresshield=shield_list(ilist,j)
            do k=1,3
            rlocshield=grad_shield_side(k,ilist,j)*eesij/fac_shield(j)
-     &     *2.0
+     &     *2.0*sss
            gshieldx(k,iresshield)=gshieldx(k,iresshield)+
      &              rlocshield
-     & +grad_shield_loc(k,ilist,j)*eesij/fac_shield(j)*2.0
+     & +grad_shield_loc(k,ilist,j)*eesij/fac_shield(j)*2.0*sss
            gshieldc(k,iresshield-1)=gshieldc(k,iresshield-1)+rlocshield
 
 C     & +grad_shield_loc(k,ilist,j)*eesij/fac_shield(j)
@@ -4111,13 +4253,13 @@ C              endif
 
           do k=1,3
             gshieldc(k,i)=gshieldc(k,i)+
-     &              grad_shield(k,i)*eesij/fac_shield(i)*2.0
+     &              grad_shield(k,i)*eesij/fac_shield(i)*2.0*sss
             gshieldc(k,j)=gshieldc(k,j)+
-     &              grad_shield(k,j)*eesij/fac_shield(j)*2.0
+     &              grad_shield(k,j)*eesij/fac_shield(j)*2.0*sss
             gshieldc(k,i-1)=gshieldc(k,i-1)+
-     &              grad_shield(k,i)*eesij/fac_shield(i)*2.0
+     &              grad_shield(k,i)*eesij/fac_shield(i)*2.0*sss
             gshieldc(k,j-1)=gshieldc(k,j-1)+
-     &              grad_shield(k,j)*eesij/fac_shield(j)*2.0
+     &              grad_shield(k,j)*eesij/fac_shield(j)*2.0*sss
 
            enddo
            endif
@@ -4148,15 +4290,10 @@ cgrad            do l=1,3
 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
-          else
-          ggg(1)=0.0
-          ggg(2)=0.0
-          ggg(3)=0.0
-          endif
+          facvdw=facvdw+sssgrad*rmij*evdwij
+          ggg(1)=facvdw*xj
+          ggg(2)=facvdw*yj
+          ggg(3)=facvdw*zj
 c          do k=1,3
 c            ghalf=0.5D0*ggg(k)
 c            gvdwpp(k,i)=gvdwpp(k,i)+ghalf
@@ -4177,10 +4314,11 @@ cgrad            enddo
 cgrad          enddo
 #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
@@ -4236,7 +4374,7 @@ cd        print '(2i3,2(3(1pd14.5),3x))',i,j,(dcosb(k),k=1,3),
 cd   &          (dcosg(k),k=1,3)
           do k=1,3
             ggg(k)=(ecosb*dcosb(k)+ecosg*dcosg(k))*
-     &      fac_shield(i)**2*fac_shield(j)**2
+     &      fac_shield(i)**2*fac_shield(j)**2*sss
           enddo
 c          do k=1,3
 c            ghalf=0.5D0*ggg(k)
@@ -4256,11 +4394,11 @@ C                     print *,"before22", gelc_long(1,i), gelc_long(1,j)
           do k=1,3
             gelc(k,i)=gelc(k,i)
      &           +((ecosa*(dc_norm(k,j)-cosa*dc_norm(k,i))
-     &           + ecosb*(erij(k)-cosb*dc_norm(k,i)))*vbld_inv(i+1))
+     &           + ecosb*(erij(k)-cosb*dc_norm(k,i)))*vbld_inv(i+1))*sss
      &           *fac_shield(i)**2*fac_shield(j)**2   
             gelc(k,j)=gelc(k,j)
      &           +((ecosa*(dc_norm(k,i)-cosa*dc_norm(k,j))
-     &           + ecosg*(erij(k)-cosg*dc_norm(k,j)))*vbld_inv(j+1))
+     &           + ecosg*(erij(k)-cosg*dc_norm(k,j)))*vbld_inv(j+1))*sss
      &           *fac_shield(i)**2*fac_shield(j)**2
             gelc_long(k,j)=gelc_long(k,j)+ggg(k)
             gelc_long(k,i)=gelc_long(k,i)-ggg(k)
@@ -4496,7 +4634,7 @@ C           fac_shield(i)=0.4
 C           fac_shield(j)=0.6
           endif
           eel_loc_ij=eel_loc_ij
-     &    *fac_shield(i)*fac_shield(j)
+     &    *fac_shield(i)*fac_shield(j)*sss
 c          if (energy_dec) write (iout,'(a6,2i5,0pf7.3)')
 c     &            'eelloc',i,j,eel_loc_ij
 C Now derivative over eel_loc
@@ -4554,7 +4692,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) 
@@ -4570,7 +4708,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)
@@ -4583,7 +4721,7 @@ c     &   a33*gmuji1(4)
 
         gloc(nphi+j,icg)=gloc(nphi+j,icg)+
      &      geel_loc_ji*wel_loc
-     &    *fac_shield(i)*fac_shield(j)
+     &    *fac_shield(i)*fac_shield(j)*sss
 
          geel_loc_ji=
      &     +a22*gmuji2(1)
@@ -4595,7 +4733,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
 
@@ -4611,17 +4749,21 @@ C Partial derivatives in virtual-bond dihedral angles gamma
      &    gel_loc_loc(i-1)=gel_loc_loc(i-1)+ 
      &            (a22*muder(1,i)*mu(1,j)+a23*muder(1,i)*mu(2,j)
      &           +a32*muder(2,i)*mu(1,j)+a33*muder(2,i)*mu(2,j))
-     &    *fac_shield(i)*fac_shield(j)
+     &    *fac_shield(i)*fac_shield(j)*sss
 
           gel_loc_loc(j-1)=gel_loc_loc(j-1)+ 
      &           (a22*mu(1,i)*muder(1,j)+a23*mu(1,i)*muder(2,j)
      &           +a32*mu(2,i)*muder(1,j)+a33*mu(2,i)*muder(2,j))
-     &    *fac_shield(i)*fac_shield(j)
+     &    *fac_shield(i)*fac_shield(j)*sss
 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)
@@ -4637,24 +4779,25 @@ C Remaining derivatives of eello
           do l=1,3
             gel_loc(l,i)=gel_loc(l,i)+(aggi(l,1)*muij(1)+
      &        aggi(l,2)*muij(2)+aggi(l,3)*muij(3)+aggi(l,4)*muij(4))
-     &    *fac_shield(i)*fac_shield(j)
+     &    *fac_shield(i)*fac_shield(j)*sss
 
             gel_loc(l,i+1)=gel_loc(l,i+1)+(aggi1(l,1)*muij(1)+
      &     aggi1(l,2)*muij(2)+aggi1(l,3)*muij(3)+aggi1(l,4)*muij(4))
-     &    *fac_shield(i)*fac_shield(j)
+     &    *fac_shield(i)*fac_shield(j)*sss
 
             gel_loc(l,j)=gel_loc(l,j)+(aggj(l,1)*muij(1)+
      &       aggj(l,2)*muij(2)+aggj(l,3)*muij(3)+aggj(l,4)*muij(4))
-     &    *fac_shield(i)*fac_shield(j)
+     &    *fac_shield(i)*fac_shield(j)*sss
 
             gel_loc(l,j1)=gel_loc(l,j1)+(aggj1(l,1)*muij(1)+
      &     aggj1(l,2)*muij(2)+aggj1(l,3)*muij(3)+aggj1(l,4)*muij(4))
-     &    *fac_shield(i)*fac_shield(j)
+     &    *fac_shield(i)*fac_shield(j)*sss
 
           enddo
           ENDIF
 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"
@@ -4738,9 +4881,9 @@ C                fac_shield(i)=0.4d0
 C                fac_shield(j)=0.6d0
                 endif
                 ees0p(num_conti,i)=0.5D0*fac3*(ees0pij+ees0mij)
-     &          *fac_shield(i)*fac_shield(j) 
+     &          *fac_shield(i)*fac_shield(j)*sss
                 ees0m(num_conti,i)=0.5D0*fac3*(ees0pij-ees0mij)
-     &          *fac_shield(i)*fac_shield(j)
+     &          *fac_shield(i)*fac_shield(j)*sss
 C Diagnostics. Comment out or remove after debugging!
 c               ees0p(num_conti,i)=0.5D0*fac3*ees0pij
 c               ees0m(num_conti,i)=0.5D0*fac3*ees0mij
@@ -4789,11 +4932,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
@@ -4808,28 +4957,28 @@ 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)
+     &          *fac_shield(i)*fac_shield(j)*sss
 
                 enddo
 C Diagnostics. Comment out or remove after debugging!
@@ -4845,6 +4994,7 @@ cdiag           enddo
               endif  ! num_conti.le.maxconts
             endif  ! fcont.gt.0
           endif    ! j.gt.i+1
+#endif
           if (wturn3.gt.0.0d0 .or. wturn4.gt.0.0d0) then
             do k=1,4
               do l=1,3
@@ -4877,7 +5027,7 @@ C Third- and fourth-order contributions from turns
       include 'COMMON.CHAIN'
       include 'COMMON.DERIV'
       include 'COMMON.INTERACT'
-      include 'COMMON.CONTACTS'
+      include 'COMMON.CORRMAT'
       include 'COMMON.TORSION'
       include 'COMMON.VECTORS'
       include 'COMMON.FFIELD'
@@ -5060,7 +5210,7 @@ C Third- and fourth-order contributions from turns
       include 'COMMON.CHAIN'
       include 'COMMON.DERIV'
       include 'COMMON.INTERACT'
-      include 'COMMON.CONTACTS'
+      include 'COMMON.CORRMAT'
       include 'COMMON.TORSION'
       include 'COMMON.VECTORS'
       include 'COMMON.FFIELD'
@@ -5428,7 +5578,10 @@ cd    write (iout,*) 'iatscp_s=',iatscp_s,' iatscp_e=',iatscp_e
 C      do xshift=-1,1
 C      do yshift=-1,1
 C      do zshift=-1,1
-      do i=iatscp_s,iatscp_e
+c      do i=iatscp_s,iatscp_e
+      do icont=g_listscp_start,g_listscp_end
+        i=newcontlistscpi(icont)
+        j=newcontlistscpj(icont)
         if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1) cycle
         iteli=itel(i)
         xi=0.5D0*(c(1,i)+c(1,i+1))
@@ -5468,9 +5621,9 @@ c        endif
 C          xi=xi+xshift*boxxsize
 C          yi=yi+yshift*boxysize
 C          zi=zi+zshift*boxzsize
-        do iint=1,nscp_gr(i)
+c        do iint=1,nscp_gr(i)
 
-        do j=iscpstart(i,iint),iscpend(i,iint)
+c        do j=iscpstart(i,iint),iscpend(i,iint)
           if (itype(j).eq.ntyp1) cycle
           itypj=iabs(itype(j))
 C Uncomment following three lines for SC-p interactions
@@ -5593,9 +5746,9 @@ cgrad          enddo
             gvdwc_scpp(k,i)=gvdwc_scpp(k,i)-ggg(k)
             gvdwc_scp(k,j)=gvdwc_scp(k,j)+ggg(k)
           enddo
-        enddo
+c        enddo
 
-        enddo ! iint
+c        enddo ! iint
       enddo ! i
 C      enddo !zshift
 C      enddo !yshift
@@ -5609,7 +5762,7 @@ C This subroutine calculates the excluded-volume interaction energy between
 C peptide-group centers and side chains and its gradient in virtual-bond and
 C side-chain vectors.
 C
-      implicit real*8 (a-h,o-z)
+      implicit none
       include 'DIMENSIONS'
       include 'COMMON.GEO'
       include 'COMMON.VAR'
@@ -5622,7 +5775,14 @@ C
       include 'COMMON.CONTROL'
       include 'COMMON.SPLITELE'
       integer xshift,yshift,zshift
-      dimension ggg(3)
+      double precision ggg(3)
+      integer i,iint,j,k,iteli,itypj,subchap,icont
+      double precision xi,yi,zi,xj,yj,zj,rrij,sss1,sssgrad1,
+     & fac,e1,e2,rij
+      double precision evdw2,evdw2_14,evdwij
+      double precision xj_safe,yj_safe,zj_safe,xj_temp,yj_temp,zj_temp,
+     & dist_temp, dist_init
+      double precision sscale,sscagrad
       evdw2=0.0D0
       evdw2_14=0.0d0
 c        print *,boxxsize,boxysize,boxzsize,'wymiary pudla'
@@ -5631,8 +5791,11 @@ cd    write (iout,*) 'iatscp_s=',iatscp_s,' iatscp_e=',iatscp_e
 C      do xshift=-1,1
 C      do yshift=-1,1
 C      do zshift=-1,1
-      if (energy_dec) write (iout,*) "escp:",r_cut,rlamb
-      do i=iatscp_s,iatscp_e
+      if (energy_dec) write (iout,*) "escp:",r_cut_int,rlamb
+c      do i=iatscp_s,iatscp_e
+      do icont=g_listscp_start,g_listscp_end
+        i=newcontlistscpi(icont)
+        j=newcontlistscpj(icont)
         if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1) cycle
         iteli=itel(i)
         xi=0.5D0*(c(1,i)+c(1,i+1))
@@ -5675,9 +5838,9 @@ c        if ((zi.gt.((zshift+0.5d0)*boxzsize)).or.
 c     &       (zi.lt.((zshift-0.5d0)*boxzsize))) then
 c        go to 136
 c        endif
-        do iint=1,nscp_gr(i)
+c        do iint=1,nscp_gr(i)
 
-        do j=iscpstart(i,iint),iscpend(i,iint)
+c        do j=iscpstart(i,iint),iscpend(i,iint)
           itypj=iabs(itype(j))
           if (itypj.eq.ntyp1) cycle
 C Uncomment following three lines for SC-p interactions
@@ -5753,11 +5916,11 @@ CHERE IS THE CALCULATION WHICH MIRROR IMAGE IS THE CLOSEST ONE
 c          print *,xj,yj,zj,'polozenie j'
           rrij=1.0D0/(xj*xj+yj*yj+zj*zj)
 c          print *,rrij
-          sss=sscale(1.0d0/(dsqrt(rrij)))
+          sss=sscale(1.0d0/(dsqrt(rrij)),r_cut_int)
 c          print *,r_cut,1.0d0/dsqrt(rrij),sss,'tu patrz'
 c          if (sss.eq.0) print *,'czasem jest OK'
           if (sss.le.0.0d0) cycle
-          sssgrad=sscagrad(1.0d0/(dsqrt(rrij)))
+          sssgrad=sscagrad(1.0d0/(dsqrt(rrij)),r_cut_int)
           fac=rrij**expon2
           e1=fac*fac*aad(itypj,iteli)
           e2=fac*bad(itypj,iteli)
@@ -5768,8 +5931,9 @@ c          if (sss.eq.0) print *,'czasem jest OK'
           endif
           evdwij=e1+e2
           evdw2=evdw2+evdwij*sss
-          if (energy_dec) write (iout,'(a6,2i5,0pf7.3,2i3,3e11.3)')
-     &        'evdw2',i,j,evdwij,iteli,itypj,fac,aad(itypj,iteli),
+          if (energy_dec) write (iout,'(a6,2i5,3f7.3,2i3,3e11.3)')
+     &        'evdw2',i,j,1.0d0/dsqrt(rrij),sss,
+     &       evdwij,iteli,itypj,fac,aad(itypj,iteli),
      &       bad(itypj,iteli)
 C
 C Calculate contributions to the gradient in the virtual-bond and SC vectors.
@@ -5811,9 +5975,9 @@ cgrad          enddo
             gvdwc_scp(k,j)=gvdwc_scp(k,j)+ggg(k)
           enddo
 c        endif !endif for sscale cutoff
-        enddo ! j
+c        enddo ! j
 
-        enddo ! iint
+c        enddo ! iint
       enddo ! i
 c      enddo !zshift
 c      enddo !yshift
@@ -7161,7 +7325,8 @@ c     &   sumene4,
 c     &   dscp1,dscp2,sumene
 c        sumene = enesc(x,xx,yy,zz,cost2tab(i+1),sint2tab(i+1))
         escloc = escloc + sumene
-c        write (2,*) "i",i," escloc",sumene,escloc,it,itype(i)
+        if (energy_dec) write (2,*) "i",i," itype",itype(i)," it",it,
+     &   " escloc",sumene,escloc,it,itype(i)
 c     & ,zz,xx,yy
 c#define DEBUG
 #ifdef DEBUG
@@ -8726,6 +8891,7 @@ c      write (iout,*) "EBACK_SC_COR",i,v1ij*cosphi+v2ij*sinphi,intertyp
 
       return
       end
+#ifdef FOURBODY
 c----------------------------------------------------------------------------
       subroutine multibody(ecorr)
 C This subroutine calculates multi-body contributions to energy following
@@ -8738,6 +8904,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
 
@@ -8792,6 +8960,8 @@ c------------------------------------------------------------------------------
       include 'COMMON.DERIV'
       include 'COMMON.INTERACT'
       include 'COMMON.CONTACTS'
+      include 'COMMON.CONTMAT'
+      include 'COMMON.CORRMAT'
       include 'COMMON.SHIELD'
       double precision gx(3),gx1(3)
       logical lprn
@@ -8846,6 +9016,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'
       include 'COMMON.CONTROL'
       include 'COMMON.LOCAL'
       double precision gx(3),gx1(3),time00
@@ -9139,6 +9311,8 @@ c------------------------------------------------------------------------------
       parameter (max_cont=maxconts)
       parameter (max_dim=26)
       include "COMMON.CONTACTS"
+      include 'COMMON.CONTMAT'
+      include 'COMMON.CORRMAT'
       double precision zapas(max_dim,maxconts,max_fg_procs),
      &  zapas_recv(max_dim,maxconts,max_fg_procs)
       common /przechowalnia/ zapas
@@ -9210,6 +9384,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'
@@ -9580,6 +9756,8 @@ c------------------------------------------------------------------------------
       parameter (max_cont=maxconts)
       parameter (max_dim=70)
       include "COMMON.CONTACTS"
+      include 'COMMON.CONTMAT'
+      include 'COMMON.CORRMAT'
       double precision zapas(max_dim,maxconts,max_fg_procs),
      &  zapas_recv(max_dim,maxconts,max_fg_procs)
       common /przechowalnia/ zapas
@@ -9633,6 +9811,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)
@@ -9808,6 +9988,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'
@@ -9873,6 +10055,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'
@@ -10259,6 +10443,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'
@@ -10380,6 +10566,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'
@@ -10784,6 +10972,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'
@@ -10924,6 +11114,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'
@@ -11028,6 +11220,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'
@@ -11213,6 +11407,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'
@@ -11328,6 +11524,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'
@@ -11572,6 +11770,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'
@@ -11890,8 +12090,8 @@ cd      write (2,*) 'ekont',ekont
 cd      write (2,*) 'eel_turn6',ekont*eel_turn6
       return
       end
-
 C-----------------------------------------------------------------------------
+#endif
       double precision function scalar(u,v)
 !DIR$ INLINEALWAYS scalar
 #ifndef OSF