update HCD-5D
[unres.git] / source / unres / src-HCD-5D / energy_p_new_barrier.F.safe
index ae8e449..2a588bd 100644 (file)
@@ -66,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
@@ -98,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)
@@ -387,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
@@ -1449,6 +1453,7 @@ C
       include 'COMMON.SBRIDGE'
       include 'COMMON.NAMES'
       include 'COMMON.IOUNITS'
+      include 'COMMON.SPLITELE'
 #ifdef FOURBODY
       include 'COMMON.CONTACTS'
       include 'COMMON.CONTMAT'
@@ -1457,8 +1462,9 @@ C
       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
+     & 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
@@ -1485,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
@@ -1498,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
       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
       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
@@ -1645,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?
@@ -1658,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
       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
+      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
@@ -1775,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
@@ -1787,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
@@ -1803,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
@@ -2108,11 +2133,10 @@ 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,r_cut_int)
-            sssgrad=sscagrad(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
@@ -2159,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
@@ -2169,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*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
@@ -2184,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
@@ -2214,6 +2238,7 @@ C
       include 'COMMON.INTERACT'
       include 'COMMON.IOUNITS'
       include 'COMMON.CALC'
+      include 'COMMON.SPLITELE'
       integer xshift,yshift,zshift,subchap
       integer icall
       common /srutu/ icall
@@ -2224,7 +2249,7 @@ C
      & 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
@@ -2384,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
@@ -2424,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
@@ -3909,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"
@@ -3933,14 +3962,33 @@ C-------------------------------------------------------------------------------
       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/
@@ -4044,8 +4092,9 @@ C        yj=yj-ymedi
 C        zj=zj-zmedi
           rij=xj*xj+yj*yj+zj*zj
 
-            sss=sscale(sqrt(rij),r_cut_int)
-            sssgrad=sscagrad(sqrt(rij),r_cut_int)
+          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)
@@ -4080,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)/)')
@@ -4089,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
@@ -4110,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     
@@ -4146,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)
@@ -4172,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
@@ -4209,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
@@ -4238,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
@@ -4297,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)
@@ -4317,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)
@@ -4557,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
@@ -4615,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) 
@@ -4631,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)
@@ -4644,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)
@@ -4656,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
 
@@ -4672,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)
@@ -4698,19 +4747,19 @@ 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
@@ -4800,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
@@ -4851,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
@@ -4870,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!
@@ -5672,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'
@@ -5685,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'
@@ -5831,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.