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=2a588bdab2630f8ae5e8c30dd3b51de03726a50d;hb=a30bd29e64da2aa47b84963fdd0bf4192ead2738;hp=65091c6df5c185a54ca8fff22958e69088399072;hpb=020e579626d686ec20ecd9f0cc4c8313f474e152;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 65091c6..2a588bd 100644
--- a/source/unres/src-HCD-5D/energy_p_new_barrier.F
+++ b/source/unres/src-HCD-5D/energy_p_new_barrier.F
@@ -1,5 +1,5 @@
       subroutine etotal(energia)
-      implicit real*8 (a-h,o-z)
+      implicit none
       include 'DIMENSIONS'
 #ifndef ISNAN
       external proc_proc
@@ -10,6 +10,8 @@ cMS$ATTRIBUTES C ::  proc_proc
 #ifdef MPI
       include "mpif.h"
       double precision weights_(n_ene)
+      double precision time00
+      integer ierror,ierr
 #endif
       include 'COMMON.SETUP'
       include 'COMMON.IOUNITS'
@@ -21,11 +23,19 @@ cMS$ATTRIBUTES C ::  proc_proc
       include 'COMMON.SBRIDGE'
       include 'COMMON.CHAIN'
       include 'COMMON.VAR'
-      include 'COMMON.MD'
+c      include 'COMMON.MD'
+      include 'COMMON.QRESTR'
       include 'COMMON.CONTROL'
       include 'COMMON.TIME1'
       include 'COMMON.SPLITELE'
       include 'COMMON.TORCNSTR'
+      include 'COMMON.SAXS'
+      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,
+     & eliptran,Eafmforce,Etube,
+     & esaxs_constr,ehomology_constr,edfator,edfanei,edfabet
+      integer n_corr,n_corr1
 #ifdef MPI      
 c      print*,"ETOTAL Processor",fg_rank," absolute rank",myrank,
 c     & " nfgtasks",nfgtasks
@@ -56,7 +66,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
@@ -88,7 +99,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)
@@ -314,6 +326,7 @@ C
       else
         esccor=0.0d0
       endif
+#ifdef FOURBODY
 C      print *,"PRZED MULIt"
 c      print *,"Processor",myrank," computed Usccorr"
 C 
@@ -342,6 +355,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
@@ -375,6 +389,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
@@ -457,7 +473,7 @@ c      print *," Processor",myrank," left SUM_ENERGY"
       end
 c-------------------------------------------------------------------------------
       subroutine sum_energy(energia,reduce)
-      implicit real*8 (a-h,o-z)
+      implicit none
       include 'DIMENSIONS'
 #ifndef ISNAN
       external proc_proc
@@ -467,6 +483,8 @@ cMS$ATTRIBUTES C ::  proc_proc
 #endif
 #ifdef MPI
       include "mpif.h"
+      integer ierr
+      double precision time00
 #endif
       include 'COMMON.SETUP'
       include 'COMMON.IOUNITS'
@@ -480,6 +498,13 @@ cMS$ATTRIBUTES C ::  proc_proc
       include 'COMMON.CONTROL'
       include 'COMMON.TIME1'
       logical reduce
+      integer i
+      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,
+     & eliptran,Eafmforce,Etube,
+     & esaxs_constr,ehomology_constr,edfator,edfanei,edfabet
+      double precision Uconst,etot
 #ifdef MPI
       if (nfgtasks.gt.1 .and. reduce) then
 #ifdef DEBUG
@@ -591,7 +616,7 @@ c detecting NaNQ
       end
 c-------------------------------------------------------------------------------
       subroutine sum_gradient
-      implicit real*8 (a-h,o-z)
+      implicit none
       include 'DIMENSIONS'
 #ifndef ISNAN
       external proc_proc
@@ -601,6 +626,8 @@ cMS$ATTRIBUTES C ::  proc_proc
 #endif
 #ifdef MPI
       include 'mpif.h'
+      integer ierror,ierr
+      double precision time00,time01
 #endif
       double precision gradbufc(3,-1:maxres),gradbufx(3,-1:maxres),
      & glocbuf(4*maxres),gradbufc_sum(3,-1:maxres)
@@ -617,7 +644,16 @@ cMS$ATTRIBUTES C ::  proc_proc
       include 'COMMON.TIME1'
       include 'COMMON.MAXGRAD'
       include 'COMMON.SCCOR'
-      include 'COMMON.MD'
+c      include 'COMMON.MD'
+      include 'COMMON.QRESTR'
+      integer i,j,k
+      double precision scalar
+      double precision gvdwc_norm,gvdwc_scp_norm,gelc_norm,gvdwpp_norm,
+     &gradb_norm,ghpbc_norm,gradcorr_norm,gel_loc_norm,gcorr3_turn_norm,
+     &gcorr4_turn_norm,gradcorr5_norm,gradcorr6_norm,
+     &gcorr6_turn_norm,gsccorrc_norm,gscloc_norm,gvdwx_norm,
+     &gradx_scp_norm,ghpbx_norm,gradxorr_norm,gsccorrx_norm,
+     &gsclocx_norm
 #ifdef TIMING
       time01=MPI_Wtime()
 #endif
@@ -1058,13 +1094,13 @@ c
       gradcorr5_max=0.0d0
       gradcorr6_max=0.0d0
       gcorr6_turn_max=0.0d0
-      gsccorc_max=0.0d0
+      gsccorrc_max=0.0d0
       gscloc_max=0.0d0
       gvdwx_max=0.0d0
       gradx_scp_max=0.0d0
       ghpbx_max=0.0d0
       gradxorr_max=0.0d0
-      gsccorx_max=0.0d0
+      gsccorrx_max=0.0d0
       gsclocx_max=0.0d0
       do i=1,nct
         gvdwc_norm=dsqrt(scalar(gvdwc(1,i),gvdwc(1,i)))
@@ -1096,13 +1132,13 @@ c
         if (gradcorr5_norm.gt.gradcorr5_max) 
      &    gradcorr5_max=gradcorr5_norm
         gradcorr6_norm=dsqrt(scalar(gradcorr6(1,i),gradcorr6(1,i)))
-        if (gradcorr6_norm.gt.gradcorr6_max) gcorr6_max=gradcorr6_norm
+        if (gradcorr6_norm.gt.gradcorr6_max)gradcorr6_max=gradcorr6_norm
         gcorr6_turn_norm=dsqrt(scalar(gcorr6_turn(1,i),
      &    gcorr6_turn(1,i)))
         if (gcorr6_turn_norm.gt.gcorr6_turn_max) 
      &    gcorr6_turn_max=gcorr6_turn_norm
-        gsccorr_norm=dsqrt(scalar(gsccorc(1,i),gsccorc(1,i)))
-        if (gsccorr_norm.gt.gsccorr_max) gsccorr_max=gsccorr_norm
+        gsccorrc_norm=dsqrt(scalar(gsccorc(1,i),gsccorc(1,i)))
+        if (gsccorrc_norm.gt.gsccorrc_max) gsccorrc_max=gsccorrc_norm
         gscloc_norm=dsqrt(scalar(gscloc(1,i),gscloc(1,i)))
         if (gscloc_norm.gt.gscloc_max) gscloc_max=gscloc_norm
         gvdwx_norm=dsqrt(scalar(gvdwx(1,i),gvdwx(1,i)))
@@ -1128,9 +1164,9 @@ c
         write (istat,'(1h#,21f10.2)') gvdwc_max,gvdwc_scp_max,
      &     gelc_max,gvdwpp_max,gradb_max,ghpbc_max,
      &     gradcorr_max,gel_loc_max,gcorr3_turn_max,gcorr4_turn_max,
-     &     gradcorr5_max,gradcorr6_max,gcorr6_turn_max,gsccorc_max,
+     &     gradcorr5_max,gradcorr6_max,gcorr6_turn_max,gsccorrc_max,
      &     gscloc_max,gvdwx_max,gradx_scp_max,ghpbx_max,gradxorr_max,
-     &     gsccorx_max,gsclocx_max
+     &     gsccorrx_max,gsclocx_max
         close(istat)
         if (gvdwc_max.gt.1.0d4) then
           write (iout,*) "gvdwc gvdwx gradb gradbx"
@@ -1157,12 +1193,18 @@ c
       end
 c-------------------------------------------------------------------------------
       subroutine rescale_weights(t_bath)
-      implicit real*8 (a-h,o-z)
+      implicit none
+#ifdef MPI
+      include 'mpif.h'
+      integer ierror
+#endif
       include 'DIMENSIONS'
       include 'COMMON.IOUNITS'
       include 'COMMON.FFIELD'
       include 'COMMON.SBRIDGE'
       include 'COMMON.CONTROL'
+      double precision t_bath
+      double precision facT,facT2,facT3,facT4,facT5
       double precision kfac /2.4d0/
       double precision x,x2,x3,x4,x5,licznik /1.12692801104297249644/
 c      facT=temp0/t_bath
@@ -1222,13 +1264,19 @@ c      write (iout,*) "t_bath",t_bath," temp0",temp0," wumb",wumb
       end
 C------------------------------------------------------------------------
       subroutine enerprint(energia)
-      implicit real*8 (a-h,o-z)
+      implicit none
       include 'DIMENSIONS'
       include 'COMMON.IOUNITS'
       include 'COMMON.FFIELD'
       include 'COMMON.SBRIDGE'
-      include 'COMMON.MD'
+      include 'COMMON.QRESTR'
       double precision energia(0:n_ene)
+      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,
+     & eello_turn6,
+     & eliptran,Eafmforce,Etube,
+     & esaxs,ehomology_constr,edfator,edfanei,edfabet,etot
       etot=energia(0)
       evdw=energia(1)
       evdw2=energia(2)
@@ -1272,10 +1320,17 @@ 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,
-     &  ethetacnstr,ebr*nss,Uconst,wumb,eliptran,wliptran,Eafmforc,
+     &  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,
      &  edfabet,wdfa_beta,
@@ -1292,13 +1347,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)'/
@@ -1319,9 +1378,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,
@@ -1338,13 +1404,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)'/
@@ -1369,7 +1439,8 @@ C
 C This subroutine calculates the interaction energy of nonbonded side chains
 C assuming the LJ potential of interaction.
 C
-      implicit real*8 (a-h,o-z)
+      implicit none
+      double precision accur
       include 'DIMENSIONS'
       parameter (accur=1.0d-10)
       include 'COMMON.GEO'
@@ -1382,8 +1453,18 @@ C
       include 'COMMON.SBRIDGE'
       include 'COMMON.NAMES'
       include 'COMMON.IOUNITS'
+      include 'COMMON.SPLITELE'
+#ifdef FOURBODY
       include 'COMMON.CONTACTS'
-      dimension gg(3)
+      include 'COMMON.CONTMAT'
+#endif
+      double precision gg(3)
+      double precision evdw,evdwij
+      integer i,j,k,itypi,itypj,itypi1,num_conti,iint
+      double precision xi,yi,zi,xj,yj,zj,rij,eps0ij,fac,e1,e2,rrij,
+     & 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
@@ -1410,6 +1491,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
@@ -1423,11 +1509,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
@@ -1443,6 +1530,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 
@@ -1498,10 +1586,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
 C Change 12/1/95
+#ifdef FOURBODY
         num_cont(i)=num_conti
+#endif
       enddo          ! i
       do i=1,nct
         do j=1,3
@@ -1526,7 +1617,7 @@ C
 C This subroutine calculates the interaction energy of nonbonded side chains
 C assuming the LJK potential of interaction.
 C
-      implicit real*8 (a-h,o-z)
+      implicit none
       include 'DIMENSIONS'
       include 'COMMON.GEO'
       include 'COMMON.VAR'
@@ -1536,8 +1627,14 @@ C
       include 'COMMON.INTERACT'
       include 'COMMON.IOUNITS'
       include 'COMMON.NAMES'
-      dimension gg(3)
+      include 'COMMON.SPLITELE'
+      double precision gg(3)
+      double precision evdw,evdwij
+      integer i,j,k,itypi,itypj,itypi1,iint
+      double precision xi,yi,zi,xj,yj,zj,rij,eps0ij,fac,e1,e2,rrij,
+     & 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
@@ -1562,6 +1659,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?
@@ -1575,11 +1675,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
@@ -1611,7 +1712,7 @@ C
 C This subroutine calculates the interaction energy of nonbonded side chains
 C assuming the Berne-Pechukas potential of interaction.
 C
-      implicit real*8 (a-h,o-z)
+      implicit none
       include 'DIMENSIONS'
       include 'COMMON.GEO'
       include 'COMMON.VAR'
@@ -1622,7 +1723,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,
+     & sss1,sssgrad1
+      double precision sscale,sscagrad
 c     double precision rrsave(maxdim)
       logical lprn
       evdw=0.0D0
@@ -1688,6 +1796,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
@@ -1700,7 +1811,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
@@ -1716,6 +1827,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
@@ -1735,7 +1847,7 @@ C
 C This subroutine calculates the interaction energy of nonbonded side chains
 C assuming the Gay-Berne potential of interaction.
 C
-      implicit real*8 (a-h,o-z)
+      implicit none
       include 'DIMENSIONS'
       include 'COMMON.GEO'
       include 'COMMON.VAR'
@@ -1750,8 +1862,14 @@ C
       include 'COMMON.SPLITELE'
       include 'COMMON.SBRIDGE'
       logical lprn
-      integer xshift,yshift,zshift
-
+      integer xshift,yshift,zshift,subchap
+      double precision evdw
+      integer itypi,itypj,itypi1,iint,ind
+      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,
+     & xj_temp,yj_temp,zj_temp,dist_temp,sig,rij_shift,faclip
+      double precision dist,sscale,sscagrad,sscagradlip,sscalelip
       evdw=0.0D0
 ccccc      energy_dec=.false.
 C      print *,'Entering EGB nnt=',nnt,' nct=',nct,' expon=',expon
@@ -2014,12 +2132,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
@@ -2066,8 +2183,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
@@ -2076,13 +2193,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
@@ -2091,8 +2208,8 @@ 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
@@ -2110,7 +2227,7 @@ C
 C This subroutine calculates the interaction energy of nonbonded side chains
 C assuming the Gay-Berne-Vorobjev potential of interaction.
 C
-      implicit real*8 (a-h,o-z)
+      implicit none
       include 'DIMENSIONS'
       include 'COMMON.GEO'
       include 'COMMON.VAR'
@@ -2121,9 +2238,19 @@ C
       include 'COMMON.INTERACT'
       include 'COMMON.IOUNITS'
       include 'COMMON.CALC'
-      integer xshift,yshift,zshift
+      include 'COMMON.SPLITELE'
+      integer xshift,yshift,zshift,subchap
+      integer icall
       common /srutu/ icall
       logical lprn
+      double precision evdw
+      integer itypi,itypj,itypi1,iint,ind
+      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,sssgrad1
+      double precision dist,sscale,sscagrad,sscagradlip,sscalelip
       evdw=0.0D0
 c     print *,'Entering EGB nnt=',nnt,' nct=',nct,' expon=',expon
       evdw=0.0D0
@@ -2282,6 +2409,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
@@ -2322,12 +2452,13 @@ 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
@@ -2477,7 +2608,7 @@ 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
@@ -2551,7 +2682,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'
@@ -2632,8 +2763,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
@@ -2858,90 +2989,6 @@ c     &   " ivec_count",(ivec_count(i),i=0,nfgtasks1-1)
 #endif
       return
       end
-C-----------------------------------------------------------------------------
-      subroutine check_vecgrad
-      implicit real*8 (a-h,o-z)
-      include 'DIMENSIONS'
-      include 'COMMON.IOUNITS'
-      include 'COMMON.GEO'
-      include 'COMMON.VAR'
-      include 'COMMON.LOCAL'
-      include 'COMMON.CHAIN'
-      include 'COMMON.VECTORS'
-      dimension uygradt(3,3,2,maxres),uzgradt(3,3,2,maxres)
-      dimension uyt(3,maxres),uzt(3,maxres)
-      dimension uygradn(3,3,2),uzgradn(3,3,2),erij(3)
-      double precision delta /1.0d-7/
-      call vec_and_deriv
-cd      do i=1,nres
-crc          write(iout,'(2i5,2(3f10.5,5x))') i,1,dc_norm(:,i)
-crc          write(iout,'(2i5,2(3f10.5,5x))') i,2,uy(:,i)
-crc          write(iout,'(2i5,2(3f10.5,5x)/)')i,3,uz(:,i)
-cd          write(iout,'(2i5,2(3f10.5,5x))') i,1,
-cd     &     (dc_norm(if90,i),if90=1,3)
-cd          write(iout,'(2i5,2(3f10.5,5x))') i,2,(uy(if90,i),if90=1,3)
-cd          write(iout,'(2i5,2(3f10.5,5x)/)')i,3,(uz(if90,i),if90=1,3)
-cd          write(iout,'(a)')
-cd      enddo
-      do i=1,nres
-        do j=1,2
-          do k=1,3
-            do l=1,3
-              uygradt(l,k,j,i)=uygrad(l,k,j,i)
-              uzgradt(l,k,j,i)=uzgrad(l,k,j,i)
-            enddo
-          enddo
-        enddo
-      enddo
-      call vec_and_deriv
-      do i=1,nres
-        do j=1,3
-          uyt(j,i)=uy(j,i)
-          uzt(j,i)=uz(j,i)
-        enddo
-      enddo
-      do i=1,nres
-cd        write (iout,*) 'i=',i
-        do k=1,3
-          erij(k)=dc_norm(k,i)
-        enddo
-        do j=1,3
-          do k=1,3
-            dc_norm(k,i)=erij(k)
-          enddo
-          dc_norm(j,i)=dc_norm(j,i)+delta
-c          fac=dsqrt(scalar(dc_norm(1,i),dc_norm(1,i)))
-c          do k=1,3
-c            dc_norm(k,i)=dc_norm(k,i)/fac
-c          enddo
-c          write (iout,*) (dc_norm(k,i),k=1,3)
-c          write (iout,*) (erij(k),k=1,3)
-          call vec_and_deriv
-          do k=1,3
-            uygradn(k,j,1)=(uy(k,i)-uyt(k,i))/delta
-            uygradn(k,j,2)=(uy(k,i-1)-uyt(k,i-1))/delta
-            uzgradn(k,j,1)=(uz(k,i)-uzt(k,i))/delta
-            uzgradn(k,j,2)=(uz(k,i-1)-uzt(k,i-1))/delta
-          enddo 
-c          write (iout,'(i5,3f8.5,3x,3f8.5,5x,3f8.5,3x,3f8.5)') 
-c     &      j,(uzgradt(k,j,1,i),k=1,3),(uzgradn(k,j,1),k=1,3),
-c     &      (uzgradt(k,j,2,i-1),k=1,3),(uzgradn(k,j,2),k=1,3)
-        enddo
-        do k=1,3
-          dc_norm(k,i)=erij(k)
-        enddo
-cd        do k=1,3
-cd          write (iout,'(i5,3f8.5,3x,3f8.5,5x,3f8.5,3x,3f8.5)') 
-cd     &      k,(uygradt(k,l,1,i),l=1,3),(uygradn(k,l,1),l=1,3),
-cd     &      (uygradt(k,l,2,i-1),l=1,3),(uygradn(k,l,2),l=1,3)
-cd          write (iout,'(i5,3f8.5,3x,3f8.5,5x,3f8.5,3x,3f8.5)') 
-cd     &      k,(uzgradt(k,l,1,i),l=1,3),(uzgradn(k,l,1),l=1,3),
-cd     &      (uzgradt(k,l,2,i-1),l=1,3),(uzgradn(k,l,2),l=1,3)
-cd          write (iout,'(a)')
-cd        enddo
-      enddo
-      return
-      end
 C--------------------------------------------------------------------------
       subroutine set_matrices
       implicit real*8 (a-h,o-z)
@@ -2959,7 +3006,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'
@@ -2975,18 +3022,26 @@ c      write(iout,*) "itype2loc",itype2loc
 #else
       do i=3,nres+1
 #endif
-        if (i.gt. nnt+2 .and. i.lt.nct+2) then
+        ii=ireschain(i-2)
+c        write (iout,*) "i",i,i-2," ii",ii
+        if (ii.eq.0) cycle
+        innt=chain_border(1,ii)
+        inct=chain_border(2,ii)
+c        write (iout,*) "i",i,i-2," ii",ii," innt",innt," inct",inct
+c        if (i.gt. nnt+2 .and. i.lt.nct+2) then 
+        if (i.gt. innt+2 .and. i.lt.inct+2) then 
           iti = itype2loc(itype(i-2))
         else
           iti=nloctyp
         endif
 c        if (i.gt. iatel_s+1 .and. i.lt.iatel_e+4) then
-        if (i.gt. nnt+1 .and. i.lt.nct+1) then
+        if (i.gt. innt+1 .and. i.lt.inct+1) then 
           iti1 = itype2loc(itype(i-1))
         else
           iti1=nloctyp
         endif
-c        write(iout,*),i
+c        write(iout,*),"i",i,i-2," iti",itype(i-2),iti,
+c     &  " iti1",itype(i-1),iti1
 #ifdef NEWCORR
         cost1=dcos(theta(i-1))
         sint1=dsin(theta(i-1))
@@ -3052,7 +3107,8 @@ c        b2tilde(2,i-2)=-b2(2,i-2)
         write (iout,*) 'theta=', theta(i-1)
 #endif
 #else
-        if (i.gt. nnt+2 .and. i.lt.nct+2) then
+        if (i.gt. innt+2 .and. i.lt.inct+2) then 
+c        if (i.gt. nnt+2 .and. i.lt.nct+2) then
           iti = itype2loc(itype(i-2))
         else
           iti=nloctyp
@@ -3097,12 +3153,14 @@ c
         write(iout,*)  'b2=',(b2(k,i-2),k=1,2)
 #endif
       enddo
+      mu=0.0d0
 #ifdef PARMAT
       do i=ivec_start+2,ivec_end+2
 #else
       do i=3,nres+1
 #endif
-        if (i .lt. nres+1) then
+c        if (itype(i-1).eq.ntyp1 .and. itype(i).eq.ntyp1) cycle
+        if (i .lt. nres+1 .and. itype(i-1).lt.ntyp1) then
           sin1=dsin(phi(i))
           cos1=dcos(phi(i))
           sintab(i-2)=sin1
@@ -3139,7 +3197,7 @@ c
           Ug2(2,1,i-2)=0.0d0
           Ug2(2,2,i-2)=0.0d0
         endif
-        if (i .gt. 3 .and. i .lt. nres+1) then
+        if (i .gt. 3) then
           obrot_der(1,i-2)=-sin1
           obrot_der(2,i-2)= cos1
           Ugder(1,1,i-2)= sin1
@@ -3169,7 +3227,8 @@ c
           Ug2der(2,2,i-2)=0.0d0
         endif
 c        if (i.gt. iatel_s+2 .and. i.lt.iatel_e+5) then
-        if (i.gt. nnt+2 .and. i.lt.nct+2) then
+c        if (i.gt. nnt+2 .and. i.lt.nct+2) then
+        if (i.gt.nnt+2 .and.i.lt.nct+2) then
           iti = itype2loc(itype(i-2))
         else
           iti=nloctyp
@@ -3198,6 +3257,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))
@@ -3206,6 +3266,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
@@ -3250,6 +3311,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))
@@ -3268,7 +3330,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)
@@ -3285,6 +3349,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
@@ -3353,6 +3418,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),
@@ -3428,6 +3494,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
@@ -3472,6 +3539,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,
@@ -3491,6 +3559,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
@@ -3543,7 +3612,7 @@ cd        enddo
 cd      enddo
       return
       end
-C--------------------------------------------------------------------------
+C-----------------------------------------------------------------------------
       subroutine eelec(ees,evdw1,eel_loc,eello_turn3,eello_turn4)
 C
 C This subroutine calculates the average interaction energy and its gradient
@@ -3566,7 +3635,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'
@@ -3639,9 +3712,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
@@ -3691,7 +3766,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
@@ -3747,12 +3824,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
@@ -3819,7 +3900,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
@@ -3835,7 +3918,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
@@ -3853,7 +3938,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"
@@ -3866,21 +3951,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/
@@ -3984,8 +4092,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)
@@ -4020,7 +4129,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)/)')
@@ -4029,11 +4138,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
@@ -4050,9 +4158,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     
@@ -4086,10 +4195,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)
@@ -4112,13 +4221,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
@@ -4149,15 +4258,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
@@ -4178,10 +4282,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
@@ -4237,7 +4342,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)
@@ -4257,11 +4362,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)
@@ -4497,7 +4602,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
@@ -4555,7 +4660,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) 
@@ -4571,7 +4676,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)
@@ -4584,7 +4689,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)
@@ -4596,7 +4701,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
 
@@ -4612,17 +4717,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)
@@ -4638,24 +4747,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"
@@ -4739,9 +4849,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
@@ -4790,11 +4900,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
@@ -4809,28 +4925,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!
@@ -4846,6 +4962,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
@@ -4878,7 +4995,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'
@@ -5061,7 +5178,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'
@@ -5610,7 +5727,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'
@@ -5623,7 +5740,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
+      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'
@@ -5632,7 +5756,7 @@ 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
+      if (energy_dec) write (iout,*) "escp:",r_cut_int,rlamb
       do i=iatscp_s,iatscp_e
         if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1) cycle
         iteli=itel(i)
@@ -5754,11 +5878,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)
@@ -5769,8 +5893,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.
@@ -6185,6 +6310,12 @@ c
       estr=0.0d0
       estr1=0.0d0
       do i=ibondp_start,ibondp_end
+c  3/4/2020 Adam: removed dummy bond graient if Calpha and SC coords are
+c      used
+#ifdef FIVEDIAG
+        if (itype(i-1).eq.ntyp1 .or. itype(i).eq.ntyp1) cycle
+        diff = vbld(i)-vbldp0
+#else
         if (itype(i-1).eq.ntyp1 .and. itype(i).eq.ntyp1) cycle
 c          estr1=estr1+gnmr1(vbld(i),-1.0d0,distchainmax)
 c          do j=1,3
@@ -6195,15 +6326,16 @@ c          if (energy_dec) write(iout,*)
 c     &       "estr1",i,gnmr1(vbld(i),-1.0d0,distchainmax)
 c        else
 C       Checking if it involves dummy (NH3+ or COO-) group
-         if (itype(i-1).eq.ntyp1 .or. itype(i).eq.ntyp1) then
+        if (itype(i-1).eq.ntyp1 .or. itype(i).eq.ntyp1) then
 C YES   vbldpDUM is the equlibrium length of spring for Dummy atom
-        diff = vbld(i)-vbldpDUM
-        if (energy_dec) write(iout,*) "dum_bond",i,diff 
-         else
-C NO    vbldp0 is the equlibrium lenght of spring for peptide group
-        diff = vbld(i)-vbldp0
-         endif 
-        if (energy_dec)    write (iout,'(a7,i5,4f7.3)') 
+          diff = vbld(i)-vbldpDUM
+          if (energy_dec) write(iout,*) "dum_bond",i,diff 
+        else
+C NO    vbldp0 is the equlibrium length of spring for peptide group
+          diff = vbld(i)-vbldp0
+        endif 
+#endif
+        if (energy_dec) write (iout,'(a7,i5,4f7.3)') 
      &     "estr bb",i,vbld(i),vbldp0,diff,AKP*diff*diff
         estr=estr+diff*diff
         do j=1,3
@@ -6713,8 +6845,6 @@ C        print *,ethetai
      &   phii1*rad2deg,ethetai
 c        lprn1=.false.
         etheta=etheta+ethetai
-        if (energy_dec) write (iout,'(a6,i5,0pf7.3)')
-     &      'ebend',i,ethetai
         if (i.gt.3) gloc(i-3,icg)=gloc(i-3,icg)+wang*dephii
         if (i.lt.nres) gloc(i-2,icg)=gloc(i-2,icg)+wang*dephii1
         gloc(nphi+i-2,icg)=gloc(nphi+i-2,icg)+wang*dethetai
@@ -7157,9 +7287,8 @@ c     &   sumene4,
 c     &   dscp1,dscp2,sumene
 c        sumene = enesc(x,xx,yy,zz,cost2tab(i+1),sint2tab(i+1))
         escloc = escloc + sumene
-        if (energy_dec) write (iout,'(a6,i5,0pf7.3)')
-     &     'escloc',i,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
@@ -7656,7 +7785,6 @@ C 6/23/01 Compute double torsional energy
       include 'COMMON.IOUNITS'
       include 'COMMON.FFIELD'
       include 'COMMON.TORCNSTR'
-      include 'COMMON.CONTROL'
       logical lprn
 C Set lprn=.true. for debugging
       lprn=.false.
@@ -7673,7 +7801,6 @@ C     &      ((itype(i).eq.ntyp1).and.(itype(i+1).eq.ntyp1))) cycle
      &  (itype(i-1).eq.ntyp1).or.(itype(i).eq.ntyp1).or.
      &  (itype(i+1).eq.ntyp1)) cycle
 C In current verion the ALL DUMMY ATOM POTENTIALS ARE OFF
-        etors_d_ii=0.0D0
         itori=itortyp(itype(i-2))
         itori1=itortyp(itype(i-1))
         itori2=itortyp(itype(i))
@@ -7708,8 +7835,6 @@ C          v1cij=v1c(1,j,itori,itori1,itori2,iblock,ntblock)
           sinphi2=dsin(j*phii1)
           etors_d=etors_d+v1cij*cosphi1+v1sij*sinphi1+
      &     v2cij*cosphi2+v2sij*sinphi2
-          if (energy_dec) etors_d_ii=etors_d_ii+
-     &     v1cij*cosphi1+v1sij*sinphi1+v2cij*cosphi2+v2sij*sinphi2
           gloci1=gloci1+j*(v1sij*cosphi1-v1cij*sinphi1)
           gloci2=gloci2+j*(v2sij*cosphi2-v2cij*sinphi2)
         enddo
@@ -7725,17 +7850,12 @@ C          v1cij=v1c(1,j,itori,itori1,itori2,iblock,ntblock)
             sinphi1m2=dsin(l*phii-(k-l)*phii1)
             etors_d=etors_d+v1cdij*cosphi1p2+v2cdij*cosphi1m2+
      &        v1sdij*sinphi1p2+v2sdij*sinphi1m2
-            if (energy_dec) etors_d_ii=etors_d_ii+
-     &        v1cdij*cosphi1p2+v2cdij*cosphi1m2+
-     &        v1sdij*sinphi1p2+v2sdij*sinphi1m2
             gloci1=gloci1+l*(v1sdij*cosphi1p2+v2sdij*cosphi1m2
      &        -v1cdij*sinphi1p2-v2cdij*sinphi1m2)
             gloci2=gloci2+(k-l)*(v1sdij*cosphi1p2-v2sdij*cosphi1m2
      &        -v1cdij*sinphi1p2+v2cdij*sinphi1m2) 
           enddo
         enddo
-          if (energy_dec) write (iout,'(a6,i5,0pf7.3)')
-     &         'etor_d',i,etors_d_ii
         gloc(i-3,icg)=gloc(i-3,icg)+wtor_d*gloci1
         gloc(i-2,icg)=gloc(i-2,icg)+wtor_d*gloci2
       enddo
@@ -7936,10 +8056,11 @@ c      do i=1,ndih_constr
 c----------------------------------------------------------------------------
 c MODELLER restraint function
       subroutine e_modeller(ehomology_constr)
-      implicit real*8 (a-h,o-z)
+      implicit none
       include 'DIMENSIONS'
 
-      integer nnn, i, j, k, ki, irec, l
+      double precision ehomology_constr
+      integer nnn,i,ii,j,k,ijk,jik,ki,kk,nexl,irec,l
       integer katy, odleglosci, test7
       real*8 odleg, odleg2, odleg3, kat, kat2, kat3, gdih(max_template)
       real*8 Eval,Erot
@@ -7955,8 +8076,11 @@ c
       double precision, dimension (max_template) ::  
      &           gtheta,dscdiff,uscdiffk,guscdiff2,guscdiff3,
      &           theta_diff
+      double precision sum_godl,sgodl,grad_odl3,ggodl,sum_gdih,
+     & sum_guscdiff,sum_sgdih,sgdih,grad_dih3,usc_diff_i,dxx,dyy,dzz,
+     & betai,sum_sgodl,dij
+      double precision dist,pinorm
 c
-
       include 'COMMON.SBRIDGE'
       include 'COMMON.CHAIN'
       include 'COMMON.GEO'
@@ -7965,8 +8089,10 @@ c
       include 'COMMON.INTERACT'
       include 'COMMON.VAR'
       include 'COMMON.IOUNITS'
-      include 'COMMON.MD'
+c      include 'COMMON.MD'
       include 'COMMON.CONTROL'
+      include 'COMMON.HOMOLOGY'
+      include 'COMMON.QRESTR'
 c
 c     From subroutine Econstr_back
 c
@@ -8682,12 +8808,12 @@ c      write (iout,*) "EBACK_SC_COR",itau_start,itau_end
       esccor=0.0D0
       do i=itau_start,itau_end
         if ((itype(i-2).eq.ntyp1).or.(itype(i-1).eq.ntyp1)) cycle
+        esccor_ii=0.0D0
         isccori=isccortyp(itype(i-2))
         isccori1=isccortyp(itype(i-1))
 c      write (iout,*) "EBACK_SC_COR",i,nterm_sccor(isccori,isccori1)
         phii=phi(i)
         do intertyp=1,3 !intertyp
-         esccor_ii=0.0D0
 cc Added 09 May 2012 (Adasko)
 cc  Intertyp means interaction type of backbone mainchain correlation: 
 c   1 = SC...Ca...Ca...Ca
@@ -8711,12 +8837,9 @@ c   3 = SC...Ca...Ca...SCi
           v2ij=v2sccor(j,intertyp,isccori,isccori1)
           cosphi=dcos(j*tauangle(intertyp,i))
           sinphi=dsin(j*tauangle(intertyp,i))
-          if (energy_dec) esccor_ii=esccor_ii+v1ij*cosphi+v2ij*sinphi
           esccor=esccor+v1ij*cosphi+v2ij*sinphi
           gloci=gloci+j*(v2ij*cosphi-v1ij*sinphi)
         enddo
-         if (energy_dec) write (iout,'(a6,i5,i2,0pf7.3)')
-     &         'esccor',i,intertyp,esccor_ii
 c      write (iout,*) "EBACK_SC_COR",i,v1ij*cosphi+v2ij*sinphi,intertyp
         gloc_sc(intertyp,i-3,icg)=gloc_sc(intertyp,i-3,icg)+wsccor*gloci
         if (lprn)
@@ -8730,6 +8853,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
@@ -8742,6 +8866,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
 
@@ -8796,6 +8922,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
@@ -8850,6 +8978,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
@@ -9143,6 +9273,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
@@ -9214,6 +9346,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'
@@ -9584,6 +9718,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
@@ -9637,6 +9773,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)
@@ -9812,6 +9950,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'
@@ -9877,6 +10017,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'
@@ -10263,6 +10405,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'
@@ -10384,6 +10528,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'
@@ -10788,6 +10934,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'
@@ -10928,6 +11076,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'
@@ -11032,6 +11182,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'
@@ -11217,6 +11369,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'
@@ -11332,6 +11486,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'
@@ -11576,6 +11732,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'
@@ -11894,8 +12052,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
@@ -12969,8 +13127,18 @@ c----------------------------------------------------------------------------
       include 'COMMON.INTERACT'
       include 'COMMON.VAR'
       include 'COMMON.IOUNITS'
-      include 'COMMON.MD'
+c      include 'COMMON.MD'
+#ifdef LANG0
+#ifdef FIVEDIAG
+      include 'COMMON.LANGEVIN.lang0.5diag'
+#else
+      include 'COMMON.LANGEVIN.lang0'
+#endif
+#else
+      include 'COMMON.LANGEVIN'
+#endif
       include 'COMMON.CONTROL'
+      include 'COMMON.SAXS'
       include 'COMMON.NAMES'
       include 'COMMON.TIME1'
       include 'COMMON.FFIELD'
@@ -13279,8 +13447,18 @@ c----------------------------------------------------------------------------
       include 'COMMON.INTERACT'
       include 'COMMON.VAR'
       include 'COMMON.IOUNITS'
-      include 'COMMON.MD'
+c      include 'COMMON.MD'
+#ifdef LANG0
+#ifdef FIVEDIAG
+      include 'COMMON.LANGEVIN.lang0.5diag'
+#else
+      include 'COMMON.LANGEVIN.lang0'
+#endif
+#else
+      include 'COMMON.LANGEVIN'
+#endif
       include 'COMMON.CONTROL'
+      include 'COMMON.SAXS'
       include 'COMMON.NAMES'
       include 'COMMON.TIME1'
       include 'COMMON.FFIELD'