source/unres/src-HCD-5D/energy_p_new-sep_barrier.F

   1 C-----------------------------------------------------------------------
   2       double precision function sscalelip(r)
   3       implicit none
   4       double precision r,gamm
   5       include "COMMON.SPLITELE"
   6 C      if(r.lt.r_cut-rlamb) then
   7 C        sscale=1.0d0
   8 C      else if(r.le.r_cut.and.r.ge.r_cut-rlamb) then
   9 C        gamm=(r-(r_cut-rlamb))/rlamb
  10         sscalelip=1.0d0+r*r*(2*r-3.0d0)
  11 C      else
  12 C        sscale=0d0
  13 C      endif
  14       return
  15       end
  16 C-----------------------------------------------------------------------
  17       double precision function sscagradlip(r)
  18       implicit none
  19       double precision r,gamm
  20       include "COMMON.SPLITELE"
  21 C     if(r.lt.r_cut-rlamb) then
  22 C        sscagrad=0.0d0
  23 C      else if(r.le.r_cut.and.r.ge.r_cut-rlamb) then
  24 C        gamm=(r-(r_cut-rlamb))/rlamb
  25         sscagradlip=r*(6*r-6.0d0)
  26 C      else
  27 C        sscagrad=0.0d0
  28 C      endif
  29       return
  30       end
  31
  32 C-----------------------------------------------------------------------
  33       double precision function sscale(r,r_cut)
  34       implicit none
  35       double precision r,r_cut,gamm
  36       include "COMMON.SPLITELE"
  37       if(r.lt.r_cut-rlamb) then
  38         sscale=1.0d0
  39       else if(r.le.r_cut.and.r.ge.r_cut-rlamb) then
  40         gamm=(r-(r_cut-rlamb))/rlamb
  41         sscale=1.0d0+gamm*gamm*(2*gamm-3.0d0)
  42       else
  43         sscale=0d0
  44       endif
  45       return
  46       end
  47 C-----------------------------------------------------------------------
  48       double precision function sscagrad(r,r_cut)
  49       implicit none
  50       double precision r,r_cut,gamm
  51       include "COMMON.SPLITELE"
  52       if(r.lt.r_cut-rlamb) then
  53         sscagrad=0.0d0
  54       else if(r.le.r_cut.and.r.ge.r_cut-rlamb) then
  55         gamm=(r-(r_cut-rlamb))/rlamb
  56         sscagrad=gamm*(6*gamm-6.0d0)/rlamb
  57       else
  58         sscagrad=0.0d0
  59       endif
  60       return
  61       end
  62 C-----------------------------------------------------------------------
  63       subroutine elj_long(evdw)
  64 C
  65 C This subroutine calculates the interaction energy of nonbonded side chains
  66 C assuming the LJ potential of interaction.
  67 C
  68       implicit none
  69       include 'DIMENSIONS'
  70       include 'COMMON.GEO'
  71       include 'COMMON.VAR'
  72       include 'COMMON.LOCAL'
  73       include 'COMMON.CHAIN'
  74       include 'COMMON.DERIV'
  75       include 'COMMON.INTERACT'
  76       include 'COMMON.TORSION'
  77       include 'COMMON.SBRIDGE'
  78       include 'COMMON.NAMES'
  79       include 'COMMON.IOUNITS'
  80       include "COMMON.SPLITELE"
  81 c      include 'COMMON.CONTACTS'
  82       double precision gg(3)
  83       double precision evdw,evdwij
  84       integer i,j,k,itypi,itypj,itypi1,num_conti,iint,ikont
  85       double precision xi,yi,zi,xj,yj,zj,rij,eps0ij,fac,e1,e2,rrij,
  86      & sigij,r0ij,rcut,sss1,sssgrad1,sqrij
  87       double precision sscale,sscagrad
  88 c      write(iout,*)'Entering ELJ nnt=',nnt,' nct=',nct,' expon=',expon
  89       evdw=0.0D0
  90 c      do i=iatsc_s,iatsc_e
  91       do ikont=g_listscsc_start,g_listscsc_end
  92         i=newcontlisti(ikont)
  93         j=newcontlistj(ikont)
  94         itypi=iabs(itype(i))
  95         if (itypi.eq.ntyp1) cycle
  96         itypi1=iabs(itype(i+1))
  97         xi=c(1,nres+i)
  98         yi=c(2,nres+i)
  99         zi=c(3,nres+i)
 100 C
 101 C Calculate SC interaction energy.
 102 C
 103 c        do iint=1,nint_gr(i)
 104 cd        write (iout,*) 'i=',i,' iint=',iint,' istart=',istart(i,iint),
 105 cd   &                  'iend=',iend(i,iint)
 106 c          do j=istart(i,iint),iend(i,iint)
 107             itypj=iabs(itype(j))
 108             if (itypj.eq.ntyp1) cycle
 109             xj=c(1,nres+j)-xi
 110             yj=c(2,nres+j)-yi
 111             zj=c(3,nres+j)-zi
 112             rij=xj*xj+yj*yj+zj*zj
 113             sqrij=dsqrt(rrij)
 114             eps0ij=eps(itypi,itypj)
 115             sss1=sscale(sqrij,r_cut_int)
 116             if (sss1.eq.0.0d0) cycle
 117             sssgrad1=sscagrad(sqrij,r_cut_int)
 118             sssgrad=
 119      &        sscagrad(sqrij/sigma(itypi,itypj),r_cut_respa)
 120             sss=sscale(sqrij/sigma(itypi,itypj),r_cut_respa)
 121             if (sss.lt.1.0d0) then
 122               rrij=1.0D0/rij
 123               fac=rrij**expon2
 124               e1=fac*fac*aa
 125               e2=fac*bb
 126               evdwij=e1+e2
 127               evdw=evdw+(1.0d0-sss)*sss1*evdwij/sqrij/expon
 128 C
 129 C Calculate the components of the gradient in DC and X
 130 C
 131               fac=-rrij*(e1+evdwij)*(1.0d0-sss)*sss1
 132      &            +evdwij*(-sss1*sssgrad/sigma(itypi,itypj)
 133      &            +(1.0d0-sss)*sssgrad1)/sqrij
 134               gg(1)=xj*fac
 135               gg(2)=yj*fac
 136               gg(3)=zj*fac
 137               do k=1,3
 138                 gvdwx(k,i)=gvdwx(k,i)-gg(k)
 139                 gvdwx(k,j)=gvdwx(k,j)+gg(k)
 140                 gvdwc(k,i)=gvdwc(k,i)-gg(k)
 141                 gvdwc(k,j)=gvdwc(k,j)+gg(k)
 142               enddo
 143             endif
 144 c          enddo      ! j
 145 c        enddo        ! iint
 146       enddo          ! i
 147       do i=1,nct
 148         do j=1,3
 149           gvdwc(j,i)=expon*gvdwc(j,i)
 150           gvdwx(j,i)=expon*gvdwx(j,i)
 151         enddo
 152       enddo
 153 C******************************************************************************
 154 C
 155 C                              N O T E !!!
 156 C
 157 C To save time, the factor of EXPON has been extracted from ALL components
 158 C of GVDWC and GRADX. Remember to multiply them by this factor before further
 159 C use!
 160 C
 161 C******************************************************************************
 162       return
 163       end
 164 C-----------------------------------------------------------------------
 165       subroutine elj_short(evdw)
 166 C
 167 C This subroutine calculates the interaction energy of nonbonded side chains
 168 C assuming the LJ potential of interaction.
 169 C
 170       implicit none
 171       include 'DIMENSIONS'
 172       include 'COMMON.GEO'
 173       include 'COMMON.VAR'
 174       include 'COMMON.LOCAL'
 175       include 'COMMON.CHAIN'
 176       include 'COMMON.DERIV'
 177       include 'COMMON.INTERACT'
 178       include 'COMMON.TORSION'
 179       include 'COMMON.SBRIDGE'
 180       include 'COMMON.NAMES'
 181       include 'COMMON.IOUNITS'
 182       include "COMMON.SPLITELE"
 183 c      include 'COMMON.CONTACTS'
 184       double precision gg(3)
 185       double precision evdw,evdwij
 186       integer i,j,k,itypi,itypj,itypi1,num_conti,iint,ikont
 187       double precision xi,yi,zi,xj,yj,zj,rij,eps0ij,fac,e1,e2,rrij,
 188      & sigij,r0ij,rcut,sqrij,sss1,sssgrad1
 189       double precision sscale,sscagrad
 190 c      write(iout,*)'Entering ELJ nnt=',nnt,' nct=',nct,' expon=',expon
 191       evdw=0.0D0
 192 c      do i=iatsc_s,iatsc_e
 193       do ikont=g_listscsc_start,g_listscsc_end
 194         i=newcontlisti(ikont)
 195         j=newcontlistj(ikont)
 196         itypi=iabs(itype(i))
 197         if (itypi.eq.ntyp1) cycle
 198         itypi1=iabs(itype(i+1))
 199         xi=c(1,nres+i)
 200         yi=c(2,nres+i)
 201         zi=c(3,nres+i)
 202 C Change 12/1/95
 203         num_conti=0
 204 C
 205 C Calculate SC interaction energy.
 206 C
 207 c        do iint=1,nint_gr(i)
 208 cd        write (iout,*) 'i=',i,' iint=',iint,' istart=',istart(i,iint),
 209 cd   &                  'iend=',iend(i,iint)
 210 c          do j=istart(i,iint),iend(i,iint)
 211             itypj=iabs(itype(j))
 212             if (itypj.eq.ntyp1) cycle
 213             xj=c(1,nres+j)-xi
 214             yj=c(2,nres+j)-yi
 215             zj=c(3,nres+j)-zi
 216 C Change 12/1/95 to calculate four-body interactions
 217             rij=xj*xj+yj*yj+zj*zj
 218             sqrij=dsqrt(rij)
 219             sss=sscale(sqrij/sigma(itypi,itypj),r_cut_respa)
 220             if (sss.gt.0.0d0) then
 221               sssgrad=
 222      &          sscagrad(sqrij/sigma(itypi,itypj),r_cut_respa)
 223               rrij=1.0D0/rij
 224               eps0ij=eps(itypi,itypj)
 225               fac=rrij**expon2
 226               e1=fac*fac*aa
 227               e2=fac*bb
 228               evdwij=e1+e2
 229               evdw=evdw+sss*evdwij
 230 C
 231 C Calculate the components of the gradient in DC and X
 232 C
 233               fac=-rrij*(e1+evdwij)*sss+evdwij*sssgrad/sqrij/expon
 234               gg(1)=xj*fac
 235               gg(2)=yj*fac
 236               gg(3)=zj*fac
 237               do k=1,3
 238                 gvdwx(k,i)=gvdwx(k,i)-gg(k)
 239                 gvdwx(k,j)=gvdwx(k,j)+gg(k)
 240                 gvdwc(k,i)=gvdwc(k,i)-gg(k)
 241                 gvdwc(k,j)=gvdwc(k,j)+gg(k)
 242               enddo
 243             endif
 244 c          enddo      ! j
 245 c        enddo        ! iint
 246       enddo          ! i
 247       do i=1,nct
 248         do j=1,3
 249           gvdwc(j,i)=expon*gvdwc(j,i)
 250           gvdwx(j,i)=expon*gvdwx(j,i)
 251         enddo
 252       enddo
 253 C******************************************************************************
 254 C
 255 C                              N O T E !!!
 256 C
 257 C To save time, the factor of EXPON has been extracted from ALL components
 258 C of GVDWC and GRADX. Remember to multiply them by this factor before further
 259 C use!
 260 C
 261 C******************************************************************************
 262       return
 263       end
 264 C-----------------------------------------------------------------------------
 265       subroutine eljk_long(evdw)
 266 C
 267 C This subroutine calculates the interaction energy of nonbonded side chains
 268 C assuming the LJK potential of interaction.
 269 C
 270       implicit none
 271       include 'DIMENSIONS'
 272       include 'COMMON.GEO'
 273       include 'COMMON.VAR'
 274       include 'COMMON.LOCAL'
 275       include 'COMMON.CHAIN'
 276       include 'COMMON.DERIV'
 277       include 'COMMON.INTERACT'
 278       include 'COMMON.IOUNITS'
 279       include 'COMMON.NAMES'
 280       include "COMMON.SPLITELE"
 281       double precision gg(3)
 282       double precision evdw,evdwij
 283       integer i,j,k,itypi,itypj,itypi1,iint,ikont
 284       double precision xi,yi,zi,xj,yj,zj,rij,eps0ij,fac,e1,e2,rrij,
 285      & fac_augm,e_augm,r_inv_ij,r_shift_inv,sss1,sssgrad1
 286       logical scheck
 287       double precision sscale,sscagrad
 288 c     print *,'Entering ELJK nnt=',nnt,' nct=',nct,' expon=',expon
 289       evdw=0.0D0
 290 c      do i=iatsc_s,iatsc_e
 291       do ikont=g_listscsc_start,g_listscsc_end
 292         i=newcontlisti(ikont)
 293         j=newcontlistj(ikont)
 294         itypi=iabs(itype(i))
 295         if (itypi.eq.ntyp1) cycle
 296         itypi1=iabs(itype(i+1))
 297         xi=c(1,nres+i)
 298         yi=c(2,nres+i)
 299         zi=c(3,nres+i)
 300 C
 301 C Calculate SC interaction energy.
 302 C
 303 c        do iint=1,nint_gr(i)
 304 c          do j=istart(i,iint),iend(i,iint)
 305             itypj=iabs(itype(j))
 306             if (itypj.eq.ntyp1) cycle
 307             xj=c(1,nres+j)-xi
 308             yj=c(2,nres+j)-yi
 309             zj=c(3,nres+j)-zi
 310             rrij=1.0D0/(xj*xj+yj*yj+zj*zj)
 311             fac_augm=rrij**expon
 312             e_augm=augm(itypi,itypj)*fac_augm
 313             r_inv_ij=dsqrt(rrij)
 314             rij=1.0D0/r_inv_ij
 315             sss1=sscale(rij,r_cut_int)
 316             if (sss1.eq.0.0d0) cycle
 317             sssgrad1=sscagrad(rij,r_cut_int)
 318             sss=sscale(rij/sigma(itypi,itypj),r_cut_respa)
 319             if (sss.lt.1.0d0) then
 320               sssgrad=
 321      &          sscagrad(rij/sigma(itypi,itypj),r_cut_respa)
 322               r_shift_inv=1.0D0/(rij+r0(itypi,itypj)-sigma(itypi,itypj))
 323               fac=r_shift_inv**expon
 324               e1=fac*fac*aa
 325               e2=fac*bb
 326               evdwij=e_augm+e1+e2
 327 cd            sigm=dabs(aa(itypi,itypj)/bb(itypi,itypj))**(1.0D0/6.0D0)
 328 cd            epsi=bb(itypi,itypj)**2/aa(itypi,itypj)
 329 cd            write (iout,'(2(a3,i3,2x),8(1pd12.4)/2(3(1pd12.4),5x)/)')
 330 cd   &          restyp(itypi),i,restyp(itypj),j,aa(itypi,itypj),
 331 cd   &          bb(itypi,itypj),augm(itypi,itypj),epsi,sigm,
 332 cd   &          sigma(itypi,itypj),1.0D0/dsqrt(rrij),evdwij,
 333 cd   &          (c(k,i),k=1,3),(c(k,j),k=1,3)
 334               evdw=evdw+(1.0d0-sss)*sss1*evdwij
 335 C
 336 C Calculate the components of the gradient in DC and X
 337 C
 338               fac=-2.0D0*rrij*e_augm-r_inv_ij*r_shift_inv*(e1+e1+e2)
 339               fac=fac*(1.0d0-sss)*sss1
 340      &          +evdwij*(-sss1*sssgrad/sigma(itypi,itypj)
 341      &          +(1.0d0-sss)*sssgrad1)*r_inv_ij/expon
 342               gg(1)=xj*fac
 343               gg(2)=yj*fac
 344               gg(3)=zj*fac
 345               do k=1,3
 346                 gvdwx(k,i)=gvdwx(k,i)-gg(k)
 347                 gvdwx(k,j)=gvdwx(k,j)+gg(k)
 348                 gvdwc(k,i)=gvdwc(k,i)-gg(k)
 349                 gvdwc(k,j)=gvdwc(k,j)+gg(k)
 350               enddo
 351             endif
 352 c          enddo      ! j
 353 c       enddo        ! iint
 354       enddo          ! i
 355       do i=1,nct
 356         do j=1,3
 357           gvdwc(j,i)=expon*gvdwc(j,i)
 358           gvdwx(j,i)=expon*gvdwx(j,i)
 359         enddo
 360       enddo
 361       return
 362       end
 363 C-----------------------------------------------------------------------------
 364       subroutine eljk_short(evdw)
 365 C
 366 C This subroutine calculates the interaction energy of nonbonded side chains
 367 C assuming the LJK potential of interaction.
 368 C
 369       implicit real*8 (a-h,o-z)
 370       include 'DIMENSIONS'
 371       include 'COMMON.GEO'
 372       include 'COMMON.VAR'
 373       include 'COMMON.LOCAL'
 374       include 'COMMON.CHAIN'
 375       include 'COMMON.DERIV'
 376       include 'COMMON.INTERACT'
 377       include 'COMMON.IOUNITS'
 378       include 'COMMON.NAMES'
 379       include "COMMON.SPLITELE"
 380       double precision gg(3)
 381       double precision evdw,evdwij
 382       integer i,j,k,itypi,itypj,itypi1,iint,ikont
 383       double precision xi,yi,zi,xj,yj,zj,rij,eps0ij,fac,e1,e2,rrij,
 384      & fac_augm,e_augm,r_inv_ij,r_shift_inv,sss1,sssgrad1
 385       logical scheck
 386       double precision sscale,sscagrad
 387 c     print *,'Entering ELJK nnt=',nnt,' nct=',nct,' expon=',expon
 388       evdw=0.0D0
 389 c      do i=iatsc_s,iatsc_e
 390       do ikont=g_listscsc_start,g_listscsc_end
 391         i=newcontlisti(ikont)
 392         j=newcontlistj(ikont)
 393         itypi=iabs(itype(i))
 394         if (itypi.eq.ntyp1) cycle
 395         itypi1=iabs(itype(i+1))
 396         xi=c(1,nres+i)
 397         yi=c(2,nres+i)
 398         zi=c(3,nres+i)
 399 C
 400 C Calculate SC interaction energy.
 401 C
 402 c        do iint=1,nint_gr(i)
 403 c          do j=istart(i,iint),iend(i,iint)
 404             itypj=iabs(itype(j))
 405             if (itypj.eq.ntyp1) cycle
 406             xj=c(1,nres+j)-xi
 407             yj=c(2,nres+j)-yi
 408             zj=c(3,nres+j)-zi
 409             rrij=1.0D0/(xj*xj+yj*yj+zj*zj)
 410             fac_augm=rrij**expon
 411             e_augm=augm(itypi,itypj)*fac_augm
 412             r_inv_ij=dsqrt(rrij)
 413             rij=1.0D0/r_inv_ij
 414             sss=sscale(rij/sigma(itypi,itypj),r_cut_respa)
 415             if (sss.gt.0.0d0) then
 416               r_shift_inv=1.0D0/(rij+r0(itypi,itypj)-sigma(itypi,itypj))
 417               fac=r_shift_inv**expon
 418               e1=fac*fac*aa
 419               e2=fac*bb
 420               evdwij=e_augm+e1+e2
 421 cd            sigm=dabs(aa(itypi,itypj)/bb(itypi,itypj))**(1.0D0/6.0D0)
 422 cd            epsi=bb(itypi,itypj)**2/aa(itypi,itypj)
 423 cd            write (iout,'(2(a3,i3,2x),8(1pd12.4)/2(3(1pd12.4),5x)/)')
 424 cd   &          restyp(itypi),i,restyp(itypj),j,aa(itypi,itypj),
 425 cd   &          bb(itypi,itypj),augm(itypi,itypj),epsi,sigm,
 426 cd   &          sigma(itypi,itypj),1.0D0/dsqrt(rrij),evdwij,
 427 cd   &          (c(k,i),k=1,3),(c(k,j),k=1,3)
 428               evdw=evdw+sss*evdwij
 429 C
 430 C Calculate the components of the gradient in DC and X
 431 C
 432               fac=-2.0D0*rrij*e_augm-r_inv_ij*r_shift_inv*(e1+e1+e2)
 433      &            +evdwij*sssgrad/sigma(itypi,itypj)*r_inv_ij/expon
 434               fac=fac*sss
 435               gg(1)=xj*fac
 436               gg(2)=yj*fac
 437               gg(3)=zj*fac
 438               do k=1,3
 439                 gvdwx(k,i)=gvdwx(k,i)-gg(k)
 440                 gvdwx(k,j)=gvdwx(k,j)+gg(k)
 441                 gvdwc(k,i)=gvdwc(k,i)-gg(k)
 442                 gvdwc(k,j)=gvdwc(k,j)+gg(k)
 443               enddo
 444             endif
 445 c          enddo      ! j
 446 c        enddo        ! iint
 447       enddo          ! i
 448       do i=1,nct
 449         do j=1,3
 450           gvdwc(j,i)=expon*gvdwc(j,i)
 451           gvdwx(j,i)=expon*gvdwx(j,i)
 452         enddo
 453       enddo
 454       return
 455       end
 456 C-----------------------------------------------------------------------------
 457       subroutine ebp_long(evdw)
 458 C
 459 C This subroutine calculates the interaction energy of nonbonded side chains
 460 C assuming the Berne-Pechukas potential of interaction.
 461 C
 462       implicit none
 463       include 'DIMENSIONS'
 464       include 'COMMON.GEO'
 465       include 'COMMON.VAR'
 466       include 'COMMON.LOCAL'
 467       include 'COMMON.CHAIN'
 468       include 'COMMON.DERIV'
 469       include 'COMMON.NAMES'
 470       include 'COMMON.INTERACT'
 471       include 'COMMON.IOUNITS'
 472       include 'COMMON.CALC'
 473       include "COMMON.SPLITELE"
 474       integer icall
 475       common /srutu/ icall
 476       double precision evdw
 477       integer itypi,itypj,itypi1,iint,ind,ikont
 478       double precision eps0ij,epsi,sigm,fac,e1,e2,rrij,xi,yi,zi
 479       double precision sss1,sssgrad1
 480       double precision sscale,sscagrad
 481 c     double precision rrsave(maxdim)
 482       logical lprn
 483       evdw=0.0D0
 484 c     print *,'Entering EBP nnt=',nnt,' nct=',nct,' expon=',expon
 485       evdw=0.0D0
 486 c     if (icall.eq.0) then
 487 c       lprn=.true.
 488 c     else
 489         lprn=.false.
 490 c     endif
 491       ind=0
 492 c      do i=iatsc_s,iatsc_e
 493       do ikont=g_listscsc_start,g_listscsc_end
 494         i=newcontlisti(ikont)
 495         j=newcontlistj(ikont)
 496         itypi=iabs(itype(i))
 497         if (itypi.eq.ntyp1) cycle
 498         itypi1=iabs(itype(i+1))
 499         xi=c(1,nres+i)
 500         yi=c(2,nres+i)
 501         zi=c(3,nres+i)
 502         dxi=dc_norm(1,nres+i)
 503         dyi=dc_norm(2,nres+i)
 504         dzi=dc_norm(3,nres+i)
 505 c        dsci_inv=dsc_inv(itypi)
 506         dsci_inv=vbld_inv(i+nres)
 507 C
 508 C Calculate SC interaction energy.
 509 C
 510 c        do iint=1,nint_gr(i)
 511 c          do j=istart(i,iint),iend(i,iint)
 512             ind=ind+1
 513             itypj=iabs(itype(j))
 514             if (itypj.eq.ntyp1) cycle
 515 c            dscj_inv=dsc_inv(itypj)
 516             dscj_inv=vbld_inv(j+nres)
 517             chi1=chi(itypi,itypj)
 518             chi2=chi(itypj,itypi)
 519             chi12=chi1*chi2
 520             chip1=chip(itypi)
 521             chip2=chip(itypj)
 522             chip12=chip1*chip2
 523             alf1=alp(itypi)
 524             alf2=alp(itypj)
 525             alf12=0.5D0*(alf1+alf2)
 526             xj=c(1,nres+j)-xi
 527             yj=c(2,nres+j)-yi
 528             zj=c(3,nres+j)-zi
 529             dxj=dc_norm(1,nres+j)
 530             dyj=dc_norm(2,nres+j)
 531             dzj=dc_norm(3,nres+j)
 532             rrij=1.0D0/(xj*xj+yj*yj+zj*zj)
 533             rij=dsqrt(rrij)
 534             sss1=sscale(1.0d0/rij,r_cut_int)
 535             if (sss1.eq.0.0d0) cycle
 536             sss=sscale(1.0d0/(rij*sigmaii(itypi,itypj)),r_cut_respa)
 537             if (sss.lt.1.0d0) then
 538               sssgrad=
 539      &          sscagrad(1.0d0/(rij*sigmaii(itypi,itypj)),r_cut_respa)
 540               sssgrad1=sscagrad(1.0d0/rij,r_cut_int)
 541 C Calculate the angle-dependent terms of energy & contributions to derivatives.
 542               call sc_angular
 543 C Calculate whole angle-dependent part of epsilon and contributions
 544 C to its derivatives
 545               fac=(rrij*sigsq)**expon2
 546               e1=fac*fac*aa
 547               e2=fac*bb
 548               evdwij=eps1*eps2rt*eps3rt*(e1+e2)
 549               eps2der=evdwij*eps3rt
 550               eps3der=evdwij*eps2rt
 551               evdwij=evdwij*eps2rt*eps3rt
 552               evdw=evdw+evdwij*(1.0d0-sss)*sss1
 553               if (lprn) then
 554               sigm=dabs(aa/bb)**(1.0D0/6.0D0)
 555               epsi=bb**2/aa
 556 cd              write (iout,'(2(a3,i3,2x),15(0pf7.3))')
 557 cd     &          restyp(itypi),i,restyp(itypj),j,
 558 cd     &          epsi,sigm,chi1,chi2,chip1,chip2,
 559 cd     &          eps1,eps2rt**2,eps3rt**2,1.0D0/dsqrt(sigsq),
 560 cd     &          om1,om2,om12,1.0D0/dsqrt(rrij),
 561 cd     &          evdwij
 562               endif
 563 C Calculate gradient components.
 564               e1=e1*eps1*eps2rt**2*eps3rt**2
 565               fac=-expon*(e1+evdwij)
 566               sigder=fac/sigsq
 567               fac=(fac+evdwij*(sss1/(1.0d0-sss)*sssgrad/
 568      &            sigmaii(itypi,itypj)+(1.0d0-sss)/sss1*sssgrad1))*rij
 569 C Calculate radial part of the gradient
 570               gg(1)=xj*fac
 571               gg(2)=yj*fac
 572               gg(3)=zj*fac
 573 C Calculate the angular part of the gradient and sum add the contributions
 574 C to the appropriate components of the Cartesian gradient.
 575               call sc_grad_scale((1.0d0-sss)*sss1)
 576             endif
 577 c          enddo      ! j
 578 c        enddo        ! iint
 579       enddo          ! i
 580 c     stop
 581       return
 582       end
 583 C-----------------------------------------------------------------------------
 584       subroutine ebp_short(evdw)
 585 C
 586 C This subroutine calculates the interaction energy of nonbonded side chains
 587 C assuming the Berne-Pechukas potential of interaction.
 588 C
 589       implicit real*8 (a-h,o-z)
 590       include 'DIMENSIONS'
 591       include 'COMMON.GEO'
 592       include 'COMMON.VAR'
 593       include 'COMMON.LOCAL'
 594       include 'COMMON.CHAIN'
 595       include 'COMMON.DERIV'
 596       include 'COMMON.NAMES'
 597       include 'COMMON.INTERACT'
 598       include 'COMMON.IOUNITS'
 599       include 'COMMON.CALC'
 600       include "COMMON.SPLITELE"
 601       integer icall
 602       common /srutu/ icall
 603       double precision evdw
 604       integer itypi,itypj,itypi1,iint,ind,ikont
 605       double precision eps0ij,epsi,sigm,fac,e1,e2,rrij,xi,yi,zi
 606       double precision sscale,sscagrad
 607 c     double precision rrsave(maxdim)
 608       logical lprn
 609       evdw=0.0D0
 610 c     print *,'Entering EBP nnt=',nnt,' nct=',nct,' expon=',expon
 611       evdw=0.0D0
 612 c     if (icall.eq.0) then
 613 c       lprn=.true.
 614 c     else
 615         lprn=.false.
 616 c     endif
 617       ind=0
 618 c      do i=iatsc_s,iatsc_e
 619       do ikont=g_listscsc_start,g_listscsc_end
 620         i=newcontlisti(ikont)
 621         j=newcontlistj(ikont)
 622         itypi=iabs(itype(i))
 623         if (itypi.eq.ntyp1) cycle
 624         itypi1=iabs(itype(i+1))
 625         xi=c(1,nres+i)
 626         yi=c(2,nres+i)
 627         zi=c(3,nres+i)
 628         dxi=dc_norm(1,nres+i)
 629         dyi=dc_norm(2,nres+i)
 630         dzi=dc_norm(3,nres+i)
 631 c        dsci_inv=dsc_inv(itypi)
 632         dsci_inv=vbld_inv(i+nres)
 633 C
 634 C Calculate SC interaction energy.
 635 C
 636 c        do iint=1,nint_gr(i)
 637 c          do j=istart(i,iint),iend(i,iint)
 638             ind=ind+1
 639             itypj=iabs(itype(j))
 640             if (itypj.eq.ntyp1) cycle
 641 c            dscj_inv=dsc_inv(itypj)
 642             dscj_inv=vbld_inv(j+nres)
 643             chi1=chi(itypi,itypj)
 644             chi2=chi(itypj,itypi)
 645             chi12=chi1*chi2
 646             chip1=chip(itypi)
 647             chip2=chip(itypj)
 648             chip12=chip1*chip2
 649             alf1=alp(itypi)
 650             alf2=alp(itypj)
 651             alf12=0.5D0*(alf1+alf2)
 652             xj=c(1,nres+j)-xi
 653             yj=c(2,nres+j)-yi
 654             zj=c(3,nres+j)-zi
 655             dxj=dc_norm(1,nres+j)
 656             dyj=dc_norm(2,nres+j)
 657             dzj=dc_norm(3,nres+j)
 658             rrij=1.0D0/(xj*xj+yj*yj+zj*zj)
 659             rij=dsqrt(rrij)
 660             sss=sscale(1.0d0/(rij*sigmaii(itypi,itypj)),r_cut_respa)
 661
 662             if (sss.gt.0.0d0) then
 663
 664 C Calculate the angle-dependent terms of energy & contributions to derivatives.
 665               call sc_angular
 666 C Calculate whole angle-dependent part of epsilon and contributions
 667 C to its derivatives
 668               fac=(rrij*sigsq)**expon2
 669               e1=fac*fac*aa
 670               e2=fac*bb
 671               evdwij=eps1*eps2rt*eps3rt*(e1+e2)
 672               eps2der=evdwij*eps3rt
 673               eps3der=evdwij*eps2rt
 674               evdwij=evdwij*eps2rt*eps3rt
 675               evdw=evdw+evdwij*sss
 676               if (lprn) then
 677               sigm=dabs(aa/bb)**(1.0D0/6.0D0)
 678               epsi=bb**2/aa
 679 cd              write (iout,'(2(a3,i3,2x),15(0pf7.3))')
 680 cd     &          restyp(itypi),i,restyp(itypj),j,
 681 cd     &          epsi,sigm,chi1,chi2,chip1,chip2,
 682 cd     &          eps1,eps2rt**2,eps3rt**2,1.0D0/dsqrt(sigsq),
 683 cd     &          om1,om2,om12,1.0D0/dsqrt(rrij),
 684 cd     &          evdwij
 685               endif
 686 C Calculate gradient components.
 687               e1=e1*eps1*eps2rt**2*eps3rt**2
 688               fac=-expon*(e1+evdwij)
 689               sigder=fac/sigsq
 690               fac=(fac+evdwij*sssgrad/sss/sigmaii(itypi,itypj))*rrij
 691 C Calculate radial part of the gradient
 692               gg(1)=xj*fac
 693               gg(2)=yj*fac
 694               gg(3)=zj*fac
 695 C Calculate the angular part of the gradient and sum add the contributions
 696 C to the appropriate components of the Cartesian gradient.
 697               call sc_grad_scale(sss)
 698             endif
 699 c          enddo      ! j
 700 c        enddo        ! iint
 701       enddo          ! i
 702 c     stop
 703       return
 704       end
 705 C-----------------------------------------------------------------------------
 706       subroutine egb_long(evdw)
 707 C
 708 C This subroutine calculates the interaction energy of nonbonded side chains
 709 C assuming the Gay-Berne potential of interaction.
 710 C
 711       implicit none
 712       include 'DIMENSIONS'
 713       include 'COMMON.GEO'
 714       include 'COMMON.VAR'
 715       include 'COMMON.LOCAL'
 716       include 'COMMON.CHAIN'
 717       include 'COMMON.DERIV'
 718       include 'COMMON.NAMES'
 719       include 'COMMON.INTERACT'
 720       include 'COMMON.IOUNITS'
 721       include 'COMMON.CALC'
 722       include 'COMMON.CONTROL'
 723       include "COMMON.SPLITELE"
 724       logical lprn
 725       integer xshift,yshift,zshift
 726       double precision evdw
 727       integer itypi,itypj,itypi1,iint,ind,ikont
 728       double precision eps0ij,epsi,sigm,fac,e1,e2,rrij,xi,yi,zi
 729       double precision fracinbuf,sslipi,evdwij_przed_tri,sig0ij,
 730      & sslipj,ssgradlipj,ssgradlipi,dist_init,xj_safe,yj_safe,zj_safe,
 731      & xj_temp,yj_temp,zj_temp,dist_temp,sig,rij_shift,faclip
 732       double precision dist,sscale,sscagrad,sscagradlip,sscalelip
 733       double precision subchap,sss1,sssgrad1
 734       evdw=0.0D0
 735 ccccc      energy_dec=.false.
 736 c     print *,'Entering EGB nnt=',nnt,' nct=',nct,' expon=',expon
 737       evdw=0.0D0
 738       lprn=.false.
 739 c     if (icall.eq.0) lprn=.false.
 740       ind=0
 741 c      do i=iatsc_s,iatsc_e
 742       do ikont=g_listscsc_start,g_listscsc_end
 743         i=newcontlisti(ikont)
 744         j=newcontlistj(ikont)
 745         itypi=iabs(itype(i))
 746         if (itypi.eq.ntyp1) cycle
 747         itypi1=iabs(itype(i+1))
 748         xi=c(1,nres+i)
 749         yi=c(2,nres+i)
 750         zi=c(3,nres+i)
 751           xi=mod(xi,boxxsize)
 752           if (xi.lt.0) xi=xi+boxxsize
 753           yi=mod(yi,boxysize)
 754           if (yi.lt.0) yi=yi+boxysize
 755           zi=mod(zi,boxzsize)
 756           if (zi.lt.0) zi=zi+boxzsize
 757         dxi=dc_norm(1,nres+i)
 758         dyi=dc_norm(2,nres+i)
 759         dzi=dc_norm(3,nres+i)
 760 c        dsci_inv=dsc_inv(itypi)
 761         dsci_inv=vbld_inv(i+nres)
 762 c        write (iout,*) "i",i,dsc_inv(itypi),dsci_inv,1.0d0/vbld(i+nres)
 763 c        write (iout,*) "dcnori",dxi*dxi+dyi*dyi+dzi*dzi
 764 C
 765 C Calculate SC interaction energy.
 766 C
 767 c        do iint=1,nint_gr(i)
 768 c          do j=istart(i,iint),iend(i,iint)
 769             ind=ind+1
 770             itypj=iabs(itype(j))
 771             if (itypj.eq.ntyp1) cycle
 772 c            dscj_inv=dsc_inv(itypj)
 773             dscj_inv=vbld_inv(j+nres)
 774 c            write (iout,*) "j",j,dsc_inv(itypj),dscj_inv,
 775 c     &       1.0d0/vbld(j+nres)
 776 c            write (iout,*) "i",i," j", j," itype",itype(i),itype(j)
 777             sig0ij=sigma(itypi,itypj)
 778             chi1=chi(itypi,itypj)
 779             chi2=chi(itypj,itypi)
 780             chi12=chi1*chi2
 781             chip1=chip(itypi)
 782             chip2=chip(itypj)
 783             chip12=chip1*chip2
 784             alf1=alp(itypi)
 785             alf2=alp(itypj)
 786             alf12=0.5D0*(alf1+alf2)
 787             xj=c(1,nres+j)
 788             yj=c(2,nres+j)
 789             zj=c(3,nres+j)
 790             xj=mod(xj,boxxsize)
 791             if (xj.lt.0) xj=xj+boxxsize
 792             yj=mod(yj,boxysize)
 793             if (yj.lt.0) yj=yj+boxysize
 794             zj=mod(zj,boxzsize)
 795             if (zj.lt.0) zj=zj+boxzsize
 796             if ((zj.gt.bordlipbot).and.(zj.lt.bordliptop)) then
 797 C the energy transfer exist
 798               if (zj.lt.buflipbot) then
 799 C what fraction I am in
 800                 fracinbuf=1.0d0-((zi-bordlipbot)/lipbufthick)
 801 C lipbufthick is thickenes of lipid buffore
 802                 sslipj=sscalelip(fracinbuf)
 803                 ssgradlipj=-sscagradlip(fracinbuf)/lipbufthick
 804               else if (zi.gt.bufliptop) then
 805                 fracinbuf=1.0d0-((bordliptop-zi)/lipbufthick)
 806                 sslipj=sscalelip(fracinbuf)
 807                 ssgradlipj=sscagradlip(fracinbuf)/lipbufthick
 808              else
 809                sslipj=1.0d0
 810                ssgradlipj=0.0
 811              endif
 812            else
 813              sslipj=0.0d0
 814              ssgradlipj=0.0
 815            endif
 816            aa=aa_lip(itypi,itypj)*(sslipi+sslipj)/2.0d0
 817      &       +aa_aq(itypi,itypj)*(2.0d0-sslipi+sslipj)/2.0d0
 818            bb=bb_lip(itypi,itypj)*(sslipi+sslipj)/2.0d0
 819      &       +bb_aq(itypi,itypj)*(2.0d0-sslipi+sslipj)/2.0d0
 820            dist_init=(xj-xi)**2+(yj-yi)**2+(zj-zi)**2
 821            xj_safe=xj
 822            yj_safe=yj
 823            zj_safe=zj
 824            subchap=0
 825            do xshift=-1,1
 826              do yshift=-1,1
 827                do zshift=-1,1
 828                  xj=xj_safe+xshift*boxxsize
 829                  yj=yj_safe+yshift*boxysize
 830                  zj=zj_safe+zshift*boxzsize
 831                  dist_temp=(xj-xi)**2+(yj-yi)**2+(zj-zi)**2
 832                  if (dist_temp.lt.dist_init) then
 833                    dist_init=dist_temp
 834                    xj_temp=xj
 835                    yj_temp=yj
 836                    zj_temp=zj
 837                    subchap=1
 838                  endif
 839                enddo
 840              enddo
 841            enddo
 842            if (subchap.eq.1) then
 843              xj=xj_temp-xi
 844              yj=yj_temp-yi
 845              zj=zj_temp-zi
 846            else
 847              xj=xj_safe-xi
 848              yj=yj_safe-yi
 849              zj=zj_safe-zi
 850            endif
 851            dxj=dc_norm(1,nres+j)
 852            dyj=dc_norm(2,nres+j)
 853            dzj=dc_norm(3,nres+j)
 854            rrij=1.0D0/(xj*xj+yj*yj+zj*zj)
 855            rij=dsqrt(rrij)
 856            sss1=sscale(1.0d0/rij,r_cut_int)
 857            if (sss1.eq.0.0d0) cycle
 858            sss=sscale(1.0d0/(rij*sigmaii(itypi,itypj)),r_cut_respa)
 859            if (sss.lt.1.0d0) then
 860 C Calculate angle-dependent terms of energy and contributions to their
 861 C derivatives.
 862              sssgrad=
 863      &         sscagrad((1.0d0/rij)/sigmaii(itypi,itypj),r_cut_respa)
 864               sssgrad1=sscagrad(1.0d0/rij,r_cut_int)
 865               call sc_angular
 866               sigsq=1.0D0/sigsq
 867               sig=sig0ij*dsqrt(sigsq)
 868               rij_shift=1.0D0/rij-sig+sig0ij
 869 c for diagnostics; uncomment
 870 c              rij_shift=1.2*sig0ij
 871 C I hate to put IF's in the loops, but here don't have another choice!!!!
 872               if (rij_shift.le.0.0D0) then
 873                 evdw=1.0D20
 874 cd                write (iout,'(2(a3,i3,2x),17(0pf7.3))')
 875 cd     &          restyp(itypi),i,restyp(itypj),j,
 876 cd     &          rij_shift,1.0D0/rij,sig,sig0ij,sigsq,1-dsqrt(sigsq)
 877                 return
 878               endif
 879               sigder=-sig*sigsq
 880 c---------------------------------------------------------------
 881               rij_shift=1.0D0/rij_shift
 882               fac=rij_shift**expon
 883               e1=fac*fac*aa
 884               e2=fac*bb
 885               evdwij=eps1*eps2rt*eps3rt*(e1+e2)
 886               eps2der=evdwij*eps3rt
 887               eps3der=evdwij*eps2rt
 888 c              write (iout,*) "sigsq",sigsq," sig",sig," eps2rt",eps2rt,
 889 c     &        " eps3rt",eps3rt," eps1",eps1," e1",e1," e2",e2
 890               evdwij=evdwij*eps2rt*eps3rt
 891               evdw=evdw+evdwij*(1.0d0-sss)*sss1
 892               if (lprn) then
 893               sigm=dabs(aa/bb)**(1.0D0/6.0D0)
 894               epsi=bb**2/aa
 895               write (iout,'(2(a3,i3,2x),17(0pf7.3))')
 896      &          restyp(itypi),i,restyp(itypj),j,
 897      &          epsi,sigm,chi1,chi2,chip1,chip2,
 898      &          eps1,eps2rt**2,eps3rt**2,sig,sig0ij,
 899      &          om1,om2,om12,1.0D0/rij,1.0D0/rij_shift,
 900      &          evdwij
 901               endif
 902
 903               if (energy_dec) write (iout,'(a6,2i5,4f10.5)')
 904      &                        'evdw',i,j,rij,sss,sss1,evdwij
 905
 906 C Calculate gradient components.
 907               e1=e1*eps1*eps2rt**2*eps3rt**2
 908               fac=-expon*(e1+evdwij)*rij_shift
 909               sigder=fac*sigder
 910               fac=(fac+evdwij*(-sss1*sssgrad/(1.0d0-sss)
 911      &            /sigmaii(itypi,itypj)+(1.0d0-sss)*sssgrad1/sss1))*rij
 912 c              fac=0.0d0
 913 C Calculate the radial part of the gradient
 914               gg(1)=xj*fac
 915               gg(2)=yj*fac
 916               gg(3)=zj*fac
 917               gg_lipi(3)=ssgradlipi*evdwij
 918               gg_lipj(3)=ssgradlipj*evdwij
 919 C Calculate angular part of the gradient.
 920               call sc_grad_scale((1.0d0-sss)*sss1)
 921             endif
 922 c          enddo      ! j
 923 c        enddo        ! iint
 924       enddo          ! i
 925 c      write (iout,*) "Number of loop steps in EGB:",ind
 926 cccc      energy_dec=.false.
 927       return
 928       end
 929 C-----------------------------------------------------------------------------
 930       subroutine egb_short(evdw)
 931 C
 932 C This subroutine calculates the interaction energy of nonbonded side chains
 933 C assuming the Gay-Berne potential of interaction.
 934 C
 935       implicit none
 936       include 'DIMENSIONS'
 937       include 'COMMON.GEO'
 938       include 'COMMON.VAR'
 939       include 'COMMON.LOCAL'
 940       include 'COMMON.CHAIN'
 941       include 'COMMON.DERIV'
 942       include 'COMMON.NAMES'
 943       include 'COMMON.INTERACT'
 944       include 'COMMON.IOUNITS'
 945       include 'COMMON.CALC'
 946       include 'COMMON.CONTROL'
 947       include "COMMON.SPLITELE"
 948       logical lprn
 949       integer xshift,yshift,zshift
 950       double precision evdw
 951       integer itypi,itypj,itypi1,iint,ind,ikont
 952       double precision eps0ij,epsi,sigm,fac,e1,e2,rrij,xi,yi,zi
 953       double precision fracinbuf,sslipi,evdwij_przed_tri,sig0ij,
 954      & sslipj,ssgradlipj,ssgradlipi,dist_init,xj_safe,yj_safe,zj_safe,
 955      & xj_temp,yj_temp,zj_temp,dist_temp,sig,rij_shift,faclip
 956       double precision dist,sscale,sscagrad,sscagradlip,sscalelip
 957       double precision subchap
 958       evdw=0.0D0
 959 ccccc      energy_dec=.false.
 960 c     print *,'Entering EGB nnt=',nnt,' nct=',nct,' expon=',expon
 961       evdw=0.0D0
 962       lprn=.false.
 963 c     if (icall.eq.0) lprn=.false.
 964       ind=0
 965 c      do i=iatsc_s,iatsc_e
 966       do ikont=g_listscsc_start,g_listscsc_end
 967         i=newcontlisti(ikont)
 968         j=newcontlistj(ikont)
 969         itypi=iabs(itype(i))
 970         if (itypi.eq.ntyp1) cycle
 971         itypi1=iabs(itype(i+1))
 972         xi=c(1,nres+i)
 973         yi=c(2,nres+i)
 974         zi=c(3,nres+i)
 975           xi=mod(xi,boxxsize)
 976           if (xi.lt.0) xi=xi+boxxsize
 977           yi=mod(yi,boxysize)
 978           if (yi.lt.0) yi=yi+boxysize
 979           zi=mod(zi,boxzsize)
 980           if (zi.lt.0) zi=zi+boxzsize
 981         dxi=dc_norm(1,nres+i)
 982         dyi=dc_norm(2,nres+i)
 983         dzi=dc_norm(3,nres+i)
 984 c        dsci_inv=dsc_inv(itypi)
 985         dsci_inv=vbld_inv(i+nres)
 986 c        write (iout,*) "i",i,dsc_inv(itypi),dsci_inv,1.0d0/vbld(i+nres)
 987 c        write (iout,*) "dcnori",dxi*dxi+dyi*dyi+dzi*dzi
 988 C
 989 C Calculate SC interaction energy.
 990 C
 991 c        do iint=1,nint_gr(i)
 992 c          do j=istart(i,iint),iend(i,iint)
 993             ind=ind+1
 994             itypj=iabs(itype(j))
 995             if (itypj.eq.ntyp1) cycle
 996 c            dscj_inv=dsc_inv(itypj)
 997             dscj_inv=vbld_inv(j+nres)
 998 c            write (iout,*) "j",j,dsc_inv(itypj),dscj_inv,
 999 c     &       1.0d0/vbld(j+nres)
1000 c            write (iout,*) "i",i," j", j," itype",itype(i),itype(j)
1001             sig0ij=sigma(itypi,itypj)
1002             chi1=chi(itypi,itypj)
1003             chi2=chi(itypj,itypi)
1004             chi12=chi1*chi2
1005             chip1=chip(itypi)
1006             chip2=chip(itypj)
1007             chip12=chip1*chip2
1008             alf1=alp(itypi)
1009             alf2=alp(itypj)
1010             alf12=0.5D0*(alf1+alf2)
1011             xj=c(1,nres+j)
1012             yj=c(2,nres+j)
1013             zj=c(3,nres+j)
1014             xj=mod(xj,boxxsize)
1015             if (xj.lt.0) xj=xj+boxxsize
1016             yj=mod(yj,boxysize)
1017             if (yj.lt.0) yj=yj+boxysize
1018             zj=mod(zj,boxzsize)
1019             if (zj.lt.0) zj=zj+boxzsize
1020             if ((zj.gt.bordlipbot).and.(zj.lt.bordliptop)) then
1021 C the energy transfer exist
1022               if (zj.lt.buflipbot) then
1023 C what fraction I am in
1024                 fracinbuf=1.0d0-((zi-bordlipbot)/lipbufthick)
1025 C lipbufthick is thickenes of lipid buffore
1026                 sslipj=sscalelip(fracinbuf)
1027                 ssgradlipj=-sscagradlip(fracinbuf)/lipbufthick
1028               elseif (zi.gt.bufliptop) then
1029                 fracinbuf=1.0d0-((bordliptop-zi)/lipbufthick)
1030                 sslipj=sscalelip(fracinbuf)
1031                 ssgradlipj=sscagradlip(fracinbuf)/lipbufthick
1032               else
1033                 sslipj=1.0d0
1034                 ssgradlipj=0.0
1035               endif
1036             else
1037               sslipj=0.0d0
1038               ssgradlipj=0.0
1039             endif
1040             aa=aa_lip(itypi,itypj)*(sslipi+sslipj)/2.0d0
1041      &        +aa_aq(itypi,itypj)*(2.0d0-sslipi+sslipj)/2.0d0
1042             bb=bb_lip(itypi,itypj)*(sslipi+sslipj)/2.0d0
1043      &        +bb_aq(itypi,itypj)*(2.0d0-sslipi+sslipj)/2.0d0
1044             dist_init=(xj-xi)**2+(yj-yi)**2+(zj-zi)**2
1045             xj_safe=xj
1046             yj_safe=yj
1047             zj_safe=zj
1048             subchap=0
1049             do xshift=-1,1
1050               do yshift=-1,1
1051                 do zshift=-1,1
1052                   xj=xj_safe+xshift*boxxsize
1053                   yj=yj_safe+yshift*boxysize
1054                   zj=zj_safe+zshift*boxzsize
1055                   dist_temp=(xj-xi)**2+(yj-yi)**2+(zj-zi)**2
1056                   if(dist_temp.lt.dist_init) then
1057                     dist_init=dist_temp
1058                     xj_temp=xj
1059                     yj_temp=yj
1060                     zj_temp=zj
1061                     subchap=1
1062                   endif
1063                 enddo
1064               enddo
1065             enddo
1066             if (subchap.eq.1) then
1067               xj=xj_temp-xi
1068               yj=yj_temp-yi
1069               zj=zj_temp-zi
1070             else
1071               xj=xj_safe-xi
1072               yj=yj_safe-yi
1073               zj=zj_safe-zi
1074             endif
1075             dxj=dc_norm(1,nres+j)
1076             dyj=dc_norm(2,nres+j)
1077             dzj=dc_norm(3,nres+j)
1078             rrij=1.0D0/(xj*xj+yj*yj+zj*zj)
1079             rij=dsqrt(rrij)
1080             sss=sscale(1.0d0/(rij*sigmaii(itypi,itypj)),r_cut_respa)
1081           sssgrad=sscagrad((1.0d0/rij)/sigmaii(itypi,itypj),r_cut_respa)
1082             if (sss.gt.0.0d0) then
1083
1084 C Calculate angle-dependent terms of energy and contributions to their
1085 C derivatives.
1086               call sc_angular
1087               sigsq=1.0D0/sigsq
1088               sig=sig0ij*dsqrt(sigsq)
1089               rij_shift=1.0D0/rij-sig+sig0ij
1090 c for diagnostics; uncomment
1091 c              rij_shift=1.2*sig0ij
1092 C I hate to put IF's in the loops, but here don't have another choice!!!!
1093               if (rij_shift.le.0.0D0) then
1094                 evdw=1.0D20
1095 cd                write (iout,'(2(a3,i3,2x),17(0pf7.3))')
1096 cd     &          restyp(itypi),i,restyp(itypj),j,
1097 cd     &          rij_shift,1.0D0/rij,sig,sig0ij,sigsq,1-dsqrt(sigsq)
1098                 return
1099               endif
1100               sigder=-sig*sigsq
1101 c---------------------------------------------------------------
1102               rij_shift=1.0D0/rij_shift
1103               fac=rij_shift**expon
1104               e1=fac*fac*aa
1105               e2=fac*bb
1106               evdwij=eps1*eps2rt*eps3rt*(e1+e2)
1107               eps2der=evdwij*eps3rt
1108               eps3der=evdwij*eps2rt
1109 c              write (iout,*) "sigsq",sigsq," sig",sig," eps2rt",eps2rt,
1110 c     &        " eps3rt",eps3rt," eps1",eps1," e1",e1," e2",e2
1111               evdwij=evdwij*eps2rt*eps3rt
1112               evdw=evdw+evdwij*sss
1113               if (lprn) then
1114               sigm=dabs(aa/bb)**(1.0D0/6.0D0)
1115               epsi=bb**2/aa
1116               write (iout,'(2(a3,i3,2x),17(0pf7.3))')
1117      &          restyp(itypi),i,restyp(itypj),j,
1118      &          epsi,sigm,chi1,chi2,chip1,chip2,
1119      &          eps1,eps2rt**2,eps3rt**2,sig,sig0ij,
1120      &          om1,om2,om12,1.0D0/rij,1.0D0/rij_shift,
1121      &          evdwij
1122               endif
1123
1124               if (energy_dec) write (iout,'(a6,2i5,0pf7.3)')
1125      &                        'evdw',i,j,evdwij
1126
1127 C Calculate gradient components.
1128               e1=e1*eps1*eps2rt**2*eps3rt**2
1129               fac=-expon*(e1+evdwij)*rij_shift
1130               sigder=fac*sigder
1131               fac=(fac+evdwij*sssgrad/sss/sigmaii(itypi,itypj))*rij
1132 c              fac=0.0d0
1133 C Calculate the radial part of the gradient
1134               gg(1)=xj*fac
1135               gg(2)=yj*fac
1136               gg(3)=zj*fac
1137               gg_lipi(3)=ssgradlipi*evdwij
1138               gg_lipj(3)=ssgradlipj*evdwij
1139 C Calculate angular part of the gradient.
1140               call sc_grad_scale(sss)
1141             endif
1142 c          enddo      ! j
1143 c        enddo        ! iint
1144       enddo          ! i
1145 c      write (iout,*) "Number of loop steps in EGB:",ind
1146 cccc      energy_dec=.false.
1147       return
1148       end
1149 C-----------------------------------------------------------------------------
1150       subroutine egbv_long(evdw)
1151 C
1152 C This subroutine calculates the interaction energy of nonbonded side chains
1153 C assuming the Gay-Berne-Vorobjev potential of interaction.
1154 C
1155       implicit real*8 (a-h,o-z)
1156       include 'DIMENSIONS'
1157       include 'COMMON.GEO'
1158       include 'COMMON.VAR'
1159       include 'COMMON.LOCAL'
1160       include 'COMMON.CHAIN'
1161       include 'COMMON.DERIV'
1162       include 'COMMON.NAMES'
1163       include 'COMMON.INTERACT'
1164       include 'COMMON.IOUNITS'
1165       include 'COMMON.CALC'
1166       include "COMMON.SPLITELE"
1167       integer icall
1168       common /srutu/ icall
1169       logical lprn
1170       integer itypi,itypj,itypi1,iint,ind,ikont
1171       double precision eps0ij,epsi,sigm,fac,e1,e2,rrij,r0ij,
1172      & xi,yi,zi,fac_augm,e_augm
1173       double precision fracinbuf,sslipi,evdwij_przed_tri,sig0ij,
1174      & sslipj,ssgradlipj,ssgradlipi,dist_init,xj_safe,yj_safe,zj_safe,
1175      & xj_temp,yj_temp,zj_temp,dist_temp,sig,rij_shift,faclip
1176       double precision dist,sscale,sscagrad,sscagradlip,sscalelip
1177       double precision sss1,sssgrad1
1178       evdw=0.0D0
1179 c     print *,'Entering EGB nnt=',nnt,' nct=',nct,' expon=',expon
1180       evdw=0.0D0
1181       lprn=.false.
1182 c     if (icall.eq.0) lprn=.true.
1183       ind=0
1184 c      do i=iatsc_s,iatsc_e
1185       do ikont=g_listscsc_start,g_listscsc_end
1186         i=newcontlisti(ikont)
1187         j=newcontlistj(ikont)
1188         itypi=iabs(itype(i))
1189         if (itypi.eq.ntyp1) cycle
1190         itypi1=iabs(itype(i+1))
1191         xi=c(1,nres+i)
1192         yi=c(2,nres+i)
1193         zi=c(3,nres+i)
1194         dxi=dc_norm(1,nres+i)
1195         dyi=dc_norm(2,nres+i)
1196         dzi=dc_norm(3,nres+i)
1197 c        dsci_inv=dsc_inv(itypi)
1198         dsci_inv=vbld_inv(i+nres)
1199 C
1200 C Calculate SC interaction energy.
1201 C
1202 c        do iint=1,nint_gr(i)
1203 c          do j=istart(i,iint),iend(i,iint)
1204             ind=ind+1
1205             itypj=iabs(itype(j))
1206             if (itypj.eq.ntyp1) cycle
1207 c            dscj_inv=dsc_inv(itypj)
1208             dscj_inv=vbld_inv(j+nres)
1209             sig0ij=sigma(itypi,itypj)
1210             r0ij=r0(itypi,itypj)
1211             chi1=chi(itypi,itypj)
1212             chi2=chi(itypj,itypi)
1213             chi12=chi1*chi2
1214             chip1=chip(itypi)
1215             chip2=chip(itypj)
1216             chip12=chip1*chip2
1217             alf1=alp(itypi)
1218             alf2=alp(itypj)
1219             alf12=0.5D0*(alf1+alf2)
1220             xj=c(1,nres+j)-xi
1221             yj=c(2,nres+j)-yi
1222             zj=c(3,nres+j)-zi
1223             dxj=dc_norm(1,nres+j)
1224             dyj=dc_norm(2,nres+j)
1225             dzj=dc_norm(3,nres+j)
1226             rrij=1.0D0/(xj*xj+yj*yj+zj*zj)
1227             rij=dsqrt(rrij)
1228
1229             sss1=sscale(1.0d0/rij,r_cut_int)
1230             if (sss1.eq.0.0d0) cycle
1231
1232             if (sss.lt.1.0d0) then
1233               sssgrad=
1234      &         sscagrad((1.0d0/rij)/sigmaii(itypi,itypj),r_cut_respa)
1235               sssgrad1=sscagrad(1.0d0/rij,r_cut_int)
1236
1237 C Calculate angle-dependent terms of energy and contributions to their
1238 C derivatives.
1239               call sc_angular
1240               sigsq=1.0D0/sigsq
1241               sig=sig0ij*dsqrt(sigsq)
1242               rij_shift=1.0D0/rij-sig+r0ij
1243 C I hate to put IF's in the loops, but here don't have another choice!!!!
1244               if (rij_shift.le.0.0D0) then
1245                 evdw=1.0D20
1246                 return
1247               endif
1248               sigder=-sig*sigsq
1249 c---------------------------------------------------------------
1250               rij_shift=1.0D0/rij_shift
1251               fac=rij_shift**expon
1252               e1=fac*fac*aa
1253               e2=fac*bb
1254               evdwij=eps1*eps2rt*eps3rt*(e1+e2)
1255               eps2der=evdwij*eps3rt
1256               eps3der=evdwij*eps2rt
1257               fac_augm=rrij**expon
1258               e_augm=augm(itypi,itypj)*fac_augm
1259               evdwij=evdwij*eps2rt*eps3rt
1260               evdw=evdw+(evdwij+e_augm)*sss1*(1.0d0-sss)
1261               if (lprn) then
1262               sigm=dabs(aa/bb)**(1.0D0/6.0D0)
1263               epsi=bb**2/aa
1264               write (iout,'(2(a3,i3,2x),17(0pf7.3))')
1265      &          restyp(itypi),i,restyp(itypj),j,
1266      &          epsi,sigm,sig,(augm(itypi,itypj)/epsi)**(1.0D0/12.0D0),
1267      &          chi1,chi2,chip1,chip2,
1268      &          eps1,eps2rt**2,eps3rt**2,
1269      &          om1,om2,om12,1.0D0/rij,1.0D0/rij_shift,
1270      &          evdwij+e_augm
1271               endif
1272 C Calculate gradient components.
1273               e1=e1*eps1*eps2rt**2*eps3rt**2
1274               fac=-expon*(e1+evdwij)*rij_shift
1275               sigder=fac*sigder
1276               fac=rij*fac-2*expon*rrij*e_augm
1277               fac=fac+(evdwij+e_augm)*
1278      &           (-sss1*sssgrad/(1.0d0-sss)/sigmaii(itypi,itypj)
1279      &            +(1.0d0-sss)*sssgrad1/sss1)*rij
1280 C Calculate the radial part of the gradient
1281               gg(1)=xj*fac
1282               gg(2)=yj*fac
1283               gg(3)=zj*fac
1284 C Calculate angular part of the gradient.
1285               call sc_grad_scale((1.0d0-sss)*sss1)
1286             endif
1287 c          enddo      ! j
1288 c        enddo        ! iint
1289       enddo          ! i
1290       end
1291 C-----------------------------------------------------------------------------
1292       subroutine egbv_short(evdw)
1293 C
1294 C This subroutine calculates the interaction energy of nonbonded side chains
1295 C assuming the Gay-Berne-Vorobjev potential of interaction.
1296 C
1297       implicit none
1298       include 'DIMENSIONS'
1299       include 'COMMON.GEO'
1300       include 'COMMON.VAR'
1301       include 'COMMON.LOCAL'
1302       include 'COMMON.CHAIN'
1303       include 'COMMON.DERIV'
1304       include 'COMMON.NAMES'
1305       include 'COMMON.INTERACT'
1306       include 'COMMON.IOUNITS'
1307       include 'COMMON.CALC'
1308       include "COMMON.SPLITELE"
1309       integer icall
1310       common /srutu/ icall
1311       logical lprn
1312       integer itypi,itypj,itypi1,iint,ind,ikont
1313       double precision eps0ij,epsi,sigm,fac,e1,e2,rrij,r0ij,
1314      & xi,yi,zi,fac_augm,e_augm
1315       double precision evdw
1316       double precision fracinbuf,sslipi,evdwij_przed_tri,sig0ij,
1317      & sslipj,ssgradlipj,ssgradlipi,dist_init,xj_safe,yj_safe,zj_safe,
1318      & xj_temp,yj_temp,zj_temp,dist_temp,sig,rij_shift,faclip
1319       double precision dist,sscale,sscagrad,sscagradlip,sscalelip
1320       evdw=0.0D0
1321 c     print *,'Entering EGB nnt=',nnt,' nct=',nct,' expon=',expon
1322       evdw=0.0D0
1323       lprn=.false.
1324 c     if (icall.eq.0) lprn=.true.
1325       ind=0
1326 c      do i=iatsc_s,iatsc_e
1327       do ikont=g_listscsc_start,g_listscsc_end
1328         i=newcontlisti(ikont)
1329         j=newcontlistj(ikont)
1330         itypi=iabs(itype(i))
1331         if (itypi.eq.ntyp1) cycle
1332         itypi1=iabs(itype(i+1))
1333         xi=c(1,nres+i)
1334         yi=c(2,nres+i)
1335         zi=c(3,nres+i)
1336         dxi=dc_norm(1,nres+i)
1337         dyi=dc_norm(2,nres+i)
1338         dzi=dc_norm(3,nres+i)
1339 c        dsci_inv=dsc_inv(itypi)
1340         dsci_inv=vbld_inv(i+nres)
1341 C
1342 C Calculate SC interaction energy.
1343 C
1344 c        do iint=1,nint_gr(i)
1345 c          do j=istart(i,iint),iend(i,iint)
1346             ind=ind+1
1347             itypj=iabs(itype(j))
1348             if (itypj.eq.ntyp1) cycle
1349 c            dscj_inv=dsc_inv(itypj)
1350             dscj_inv=vbld_inv(j+nres)
1351             sig0ij=sigma(itypi,itypj)
1352             r0ij=r0(itypi,itypj)
1353             chi1=chi(itypi,itypj)
1354             chi2=chi(itypj,itypi)
1355             chi12=chi1*chi2
1356             chip1=chip(itypi)
1357             chip2=chip(itypj)
1358             chip12=chip1*chip2
1359             alf1=alp(itypi)
1360             alf2=alp(itypj)
1361             alf12=0.5D0*(alf1+alf2)
1362             xj=c(1,nres+j)-xi
1363             yj=c(2,nres+j)-yi
1364             zj=c(3,nres+j)-zi
1365             dxj=dc_norm(1,nres+j)
1366             dyj=dc_norm(2,nres+j)
1367             dzj=dc_norm(3,nres+j)
1368             rrij=1.0D0/(xj*xj+yj*yj+zj*zj)
1369             rij=dsqrt(rrij)
1370
1371             sss=sscale(1.0d0/(rij*sigmaii(itypi,itypj)),r_cut_respa)
1372
1373             if (sss.gt.0.0d0) then
1374
1375 C Calculate angle-dependent terms of energy and contributions to their
1376 C derivatives.
1377               call sc_angular
1378               sigsq=1.0D0/sigsq
1379               sig=sig0ij*dsqrt(sigsq)
1380               rij_shift=1.0D0/rij-sig+r0ij
1381 C I hate to put IF's in the loops, but here don't have another choice!!!!
1382               if (rij_shift.le.0.0D0) then
1383                 evdw=1.0D20
1384                 return
1385               endif
1386               sigder=-sig*sigsq
1387 c---------------------------------------------------------------
1388               rij_shift=1.0D0/rij_shift
1389               fac=rij_shift**expon
1390               e1=fac*fac*aa
1391               e2=fac*bb
1392               evdwij=eps1*eps2rt*eps3rt*(e1+e2)
1393               eps2der=evdwij*eps3rt
1394               eps3der=evdwij*eps2rt
1395               fac_augm=rrij**expon
1396               e_augm=augm(itypi,itypj)*fac_augm
1397               evdwij=evdwij*eps2rt*eps3rt
1398               evdw=evdw+(evdwij+e_augm)*sss
1399               if (lprn) then
1400               sigm=dabs(aa/bb)**(1.0D0/6.0D0)
1401               epsi=bb**2/aa
1402               write (iout,'(2(a3,i3,2x),17(0pf7.3))')
1403      &          restyp(itypi),i,restyp(itypj),j,
1404      &          epsi,sigm,sig,(augm(itypi,itypj)/epsi)**(1.0D0/12.0D0),
1405      &          chi1,chi2,chip1,chip2,
1406      &          eps1,eps2rt**2,eps3rt**2,
1407      &          om1,om2,om12,1.0D0/rij,1.0D0/rij_shift,
1408      &          evdwij+e_augm
1409               endif
1410 C Calculate gradient components.
1411               e1=e1*eps1*eps2rt**2*eps3rt**2
1412               fac=-expon*(e1+evdwij)*rij_shift
1413               sigder=fac*sigder
1414               fac=rij*fac-2*expon*rrij*e_augm+
1415      &          (evdwij+e_augm)*sssgrad/sigmaii(itypi,itypj)/sss*rij
1416 C Calculate the radial part of the gradient
1417               gg(1)=xj*fac
1418               gg(2)=yj*fac
1419               gg(3)=zj*fac
1420 C Calculate angular part of the gradient.
1421               call sc_grad_scale(sss)
1422             endif
1423 c          enddo      ! j
1424 c        enddo        ! iint
1425       enddo          ! i
1426       end
1427 C----------------------------------------------------------------------------
1428       subroutine sc_grad_scale(scalfac)
1429       implicit real*8 (a-h,o-z)
1430       include 'DIMENSIONS'
1431       include 'COMMON.CHAIN'
1432       include 'COMMON.DERIV'
1433       include 'COMMON.CALC'
1434       include 'COMMON.IOUNITS'
1435       include "COMMON.SPLITELE"
1436       double precision dcosom1(3),dcosom2(3)
1437       double precision scalfac
1438       eom1=eps2der*eps2rt_om1-2.0D0*alf1*eps3der+sigder*sigsq_om1
1439       eom2=eps2der*eps2rt_om2+2.0D0*alf2*eps3der+sigder*sigsq_om2
1440       eom12=evdwij*eps1_om12+eps2der*eps2rt_om12
1441      &     -2.0D0*alf12*eps3der+sigder*sigsq_om12
1442 c diagnostics only
1443 c      eom1=0.0d0
1444 c      eom2=0.0d0
1445 c      eom12=evdwij*eps1_om12
1446 c end diagnostics
1447 c      write (iout,*) "eps2der",eps2der," eps3der",eps3der,
1448 c     &  " sigder",sigder
1449 c      write (iout,*) "eps1_om12",eps1_om12," eps2rt_om12",eps2rt_om12
1450 c      write (iout,*) "eom1",eom1," eom2",eom2," eom12",eom12
1451       do k=1,3
1452         dcosom1(k)=rij*(dc_norm(k,nres+i)-om1*erij(k))
1453         dcosom2(k)=rij*(dc_norm(k,nres+j)-om2*erij(k))
1454       enddo
1455       do k=1,3
1456         gg(k)=(gg(k)+eom1*dcosom1(k)+eom2*dcosom2(k))*scalfac
1457       enddo
1458 c      write (iout,*) "gg",(gg(k),k=1,3)
1459       do k=1,3
1460         gvdwx(k,i)=gvdwx(k,i)-gg(k)+gg_lipi(k)
1461      &            +(eom12*(dc_norm(k,nres+j)-om12*dc_norm(k,nres+i))
1462      &          +eom1*(erij(k)-om1*dc_norm(k,nres+i)))*dsci_inv*scalfac
1463         gvdwx(k,j)=gvdwx(k,j)+gg(k)+gg_lipj(k)
1464      &            +(eom12*(dc_norm(k,nres+i)-om12*dc_norm(k,nres+j))
1465      &          +eom2*(erij(k)-om2*dc_norm(k,nres+j)))*dscj_inv*scalfac
1466 c        write (iout,*)(eom12*(dc_norm(k,nres+j)-om12*dc_norm(k,nres+i))
1467 c     &            +eom1*(erij(k)-om1*dc_norm(k,nres+i)))*dsci_inv
1468 c        write (iout,*)(eom12*(dc_norm(k,nres+i)-om12*dc_norm(k,nres+j))
1469 c     &            +eom2*(erij(k)-om2*dc_norm(k,nres+j)))*dscj_inv
1470       enddo
1471 C
1472 C Calculate the components of the gradient in DC and X
1473 C
1474       do l=1,3
1475         gvdwc(l,i)=gvdwc(l,i)-gg(l)+gg_lipi(l)
1476         gvdwc(l,j)=gvdwc(l,j)+gg(l)+gg_lipj(l)
1477       enddo
1478       return
1479       end
1480 C--------------------------------------------------------------------------
1481       subroutine eelec_scale(ees,evdw1,eel_loc,eello_turn3,eello_turn4)
1482 C
1483 C This subroutine calculates the average interaction energy and its gradient
1484 C in the virtual-bond vectors between non-adjacent peptide groups, based on
1485 C the potential described in Liwo et al., Protein Sci., 1993, 2, 1715.
1486 C The potential depends both on the distance of peptide-group centers and on
1487 C the orientation of the CA-CA virtual bonds.
1488 C
1489       implicit real*8 (a-h,o-z)
1490 #ifdef MPI
1491       include 'mpif.h'
1492 #endif
1493       include 'DIMENSIONS'
1494       include 'COMMON.CONTROL'
1495       include 'COMMON.SETUP'
1496       include 'COMMON.IOUNITS'
1497       include 'COMMON.GEO'
1498       include 'COMMON.VAR'
1499       include 'COMMON.LOCAL'
1500       include 'COMMON.CHAIN'
1501       include 'COMMON.DERIV'
1502       include 'COMMON.INTERACT'
1503 #ifdef FOURBODY
1504       include 'COMMON.CONTACTS'
1505       include 'COMMON.CONTMAT'
1506 #endif
1507       include 'COMMON.CORRMAT'
1508       include 'COMMON.TORSION'
1509       include 'COMMON.VECTORS'
1510       include 'COMMON.FFIELD'
1511       include 'COMMON.TIME1'
1512       include 'COMMON.SHIELD'
1513       include "COMMON.SPLITELE"
1514       integer ikont
1515       dimension ggg(3),gggp(3),gggm(3),erij(3),dcosb(3),dcosg(3),
1516      &          erder(3,3),uryg(3,3),urzg(3,3),vryg(3,3),vrzg(3,3)
1517       double precision acipa(2,2),agg(3,4),aggi(3,4),aggi1(3,4),
1518      &    aggj(3,4),aggj1(3,4),a_temp(2,2),muij(4)
1519       common /locel/ a_temp,agg,aggi,aggi1,aggj,aggj1,a22,a23,a32,a33,
1520      &    dxi,dyi,dzi,dx_normi,dy_normi,dz_normi,xmedi,ymedi,zmedi,
1521      &    num_conti,j1,j2
1522 c 4/26/02 - AL scaling factor for 1,4 repulsive VDW interactions
1523 #ifdef MOMENT
1524       double precision scal_el /1.0d0/
1525 #else
1526       double precision scal_el /0.5d0/
1527 #endif
1528 C 12/13/98
1529 C 13-go grudnia roku pamietnego...
1530       double precision unmat(3,3) /1.0d0,0.0d0,0.0d0,
1531      &                   0.0d0,1.0d0,0.0d0,
1532      &                   0.0d0,0.0d0,1.0d0/
1533 cd      write(iout,*) 'In EELEC'
1534 cd      do i=1,nloctyp
1535 cd        write(iout,*) 'Type',i
1536 cd        write(iout,*) 'B1',B1(:,i)
1537 cd        write(iout,*) 'B2',B2(:,i)
1538 cd        write(iout,*) 'CC',CC(:,:,i)
1539 cd        write(iout,*) 'DD',DD(:,:,i)
1540 cd        write(iout,*) 'EE',EE(:,:,i)
1541 cd      enddo
1542 cd      call check_vecgrad
1543 cd      stop
1544 C      print *,"WCHODZE3"
1545       if (icheckgrad.eq.1) then
1546         do i=1,nres-1
1547           fac=1.0d0/dsqrt(scalar(dc(1,i),dc(1,i)))
1548           do k=1,3
1549             dc_norm(k,i)=dc(k,i)*fac
1550           enddo
1551 c          write (iout,*) 'i',i,' fac',fac
1552         enddo
1553       endif
1554       if (wel_loc.gt.0.0d0 .or. wcorr4.gt.0.0d0 .or. wcorr5.gt.0.0d0
1555      &    .or. wcorr6.gt.0.0d0 .or. wturn3.gt.0.0d0 .or.
1556      &    wturn4.gt.0.0d0 .or. wturn6.gt.0.0d0) then
1557 c        call vec_and_deriv
1558 #ifdef TIMING
1559         time01=MPI_Wtime()
1560 #endif
1561         call set_matrices
1562 #ifdef TIMING
1563         time_mat=time_mat+MPI_Wtime()-time01
1564 #endif
1565       endif
1566 cd      do i=1,nres-1
1567 cd        write (iout,*) 'i=',i
1568 cd        do k=1,3
1569 cd        write (iout,'(i5,2f10.5)') k,uy(k,i),uz(k,i)
1570 cd        enddo
1571 cd        do k=1,3
1572 cd          write (iout,'(f10.5,2x,3f10.5,2x,3f10.5)')
1573 cd     &     uz(k,i),(uzgrad(k,l,1,i),l=1,3),(uzgrad(k,l,2,i),l=1,3)
1574 cd        enddo
1575 cd      enddo
1576       t_eelecij=0.0d0
1577       ees=0.0D0
1578       evdw1=0.0D0
1579       eel_loc=0.0d0
1580       eello_turn3=0.0d0
1581       eello_turn4=0.0d0
1582       ind=0
1583 #ifdef FOURBODY
1584       do i=1,nres
1585         num_cont_hb(i)=0
1586       enddo
1587 #endif
1588 cd      print '(a)','Enter EELEC'
1589 cd      write (iout,*) 'iatel_s=',iatel_s,' iatel_e=',iatel_e
1590       do i=1,nres
1591         gel_loc_loc(i)=0.0d0
1592         gcorr_loc(i)=0.0d0
1593       enddo
1594 c
1595 c
1596 c 9/27/08 AL Split the interaction loop to ensure load balancing of turn terms
1597 C
1598 C Loop over i,i+2 and i,i+3 pairs of the peptide groups
1599 C
1600       do i=iturn3_start,iturn3_end
1601         if (itype(i).eq.ntyp1.or. itype(i+1).eq.ntyp1
1602      &  .or. itype(i+2).eq.ntyp1 .or. itype(i+3).eq.ntyp1
1603 C     &  .or. itype(i-1).eq.ntyp1
1604 C     &  .or. itype(i+4).eq.ntyp1
1605      &   ) cycle
1606         dxi=dc(1,i)
1607         dyi=dc(2,i)
1608         dzi=dc(3,i)
1609         dx_normi=dc_norm(1,i)
1610         dy_normi=dc_norm(2,i)
1611         dz_normi=dc_norm(3,i)
1612         xmedi=c(1,i)+0.5d0*dxi
1613         ymedi=c(2,i)+0.5d0*dyi
1614         zmedi=c(3,i)+0.5d0*dzi
1615           xmedi=mod(xmedi,boxxsize)
1616           if (xmedi.lt.0) xmedi=xmedi+boxxsize
1617           ymedi=mod(ymedi,boxysize)
1618           if (ymedi.lt.0) ymedi=ymedi+boxysize
1619           zmedi=mod(zmedi,boxzsize)
1620           if (zmedi.lt.0) zmedi=zmedi+boxzsize
1621         num_conti=0
1622         call eelecij_scale(i,i+2,ees,evdw1,eel_loc)
1623         if (wturn3.gt.0.0d0) call eturn3(i,eello_turn3)
1624 #ifdef FOURBODY
1625         num_cont_hb(i)=num_conti
1626 #endif
1627       enddo
1628       do i=iturn4_start,iturn4_end
1629         if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1
1630      &    .or. itype(i+3).eq.ntyp1
1631      &    .or. itype(i+4).eq.ntyp1
1632 C     &    .or. itype(i+5).eq.ntyp1
1633 C     &    .or. itype(i-1).eq.ntyp1
1634      &    ) cycle
1635         dxi=dc(1,i)
1636         dyi=dc(2,i)
1637         dzi=dc(3,i)
1638         dx_normi=dc_norm(1,i)
1639         dy_normi=dc_norm(2,i)
1640         dz_normi=dc_norm(3,i)
1641         xmedi=c(1,i)+0.5d0*dxi
1642         ymedi=c(2,i)+0.5d0*dyi
1643         zmedi=c(3,i)+0.5d0*dzi
1644           xmedi=mod(xmedi,boxxsize)
1645           if (xmedi.lt.0) xmedi=xmedi+boxxsize
1646           ymedi=mod(ymedi,boxysize)
1647           if (ymedi.lt.0) ymedi=ymedi+boxysize
1648           zmedi=mod(zmedi,boxzsize)
1649           if (zmedi.lt.0) zmedi=zmedi+boxzsize
1650 #ifdef FOURBODY
1651         num_conti=num_cont_hb(i)
1652 #endif
1653         call eelecij_scale(i,i+3,ees,evdw1,eel_loc)
1654         if (wturn4.gt.0.0d0 .and. itype(i+2).ne.ntyp1)
1655      &    call eturn4(i,eello_turn4)
1656 #ifdef FOURBODY
1657         num_cont_hb(i)=num_conti
1658 #endif
1659       enddo   ! i
1660 c
1661 c Loop over all pairs of interacting peptide groups except i,i+2 and i,i+3
1662 c
1663 c      do i=iatel_s,iatel_e
1664       do ikont=g_listpp_start,g_listpp_end
1665         i=newcontlistppi(ikont)
1666         j=newcontlistppj(ikont)
1667         if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1
1668 C     &  .or. itype(i+2).eq.ntyp1
1669 C     &  .or. itype(i-1).eq.ntyp1
1670      &) cycle
1671         dxi=dc(1,i)
1672         dyi=dc(2,i)
1673         dzi=dc(3,i)
1674         dx_normi=dc_norm(1,i)
1675         dy_normi=dc_norm(2,i)
1676         dz_normi=dc_norm(3,i)
1677         xmedi=c(1,i)+0.5d0*dxi
1678         ymedi=c(2,i)+0.5d0*dyi
1679         zmedi=c(3,i)+0.5d0*dzi
1680           xmedi=mod(xmedi,boxxsize)
1681           if (xmedi.lt.0) xmedi=xmedi+boxxsize
1682           ymedi=mod(ymedi,boxysize)
1683           if (ymedi.lt.0) ymedi=ymedi+boxysize
1684           zmedi=mod(zmedi,boxzsize)
1685           if (zmedi.lt.0) zmedi=zmedi+boxzsize
1686 c        write (iout,*) 'i',i,' ielstart',ielstart(i),' ielend',ielend
1687 #ifdef FOURBODY
1688         num_conti=num_cont_hb(i)
1689 #endif
1690 c        do j=ielstart(i),ielend(i)
1691           if (itype(j).eq.ntyp1 .or. itype(j+1).eq.ntyp1
1692 C     & .or.itype(j+2).eq.ntyp1
1693 C     & .or.itype(j-1).eq.ntyp1
1694      &) cycle
1695           call eelecij_scale(i,j,ees,evdw1,eel_loc)
1696 c        enddo ! j
1697 #ifdef FOURBODY
1698         num_cont_hb(i)=num_conti
1699 #endif
1700       enddo   ! i
1701 c      write (iout,*) "Number of loop steps in EELEC:",ind
1702 cd      do i=1,nres
1703 cd        write (iout,'(i3,3f10.5,5x,3f10.5)')
1704 cd     &     i,(gel_loc(k,i),k=1,3),gel_loc_loc(i)
1705 cd      enddo
1706 c 12/7/99 Adam eello_turn3 will be considered as a separate energy term
1707 ccc      eel_loc=eel_loc+eello_turn3
1708 cd      print *,"Processor",fg_rank," t_eelecij",t_eelecij
1709       return
1710       end
1711 C-------------------------------------------------------------------------------
1712       subroutine eelecij_scale(i,j,ees,evdw1,eel_loc)
1713       implicit none
1714       include 'DIMENSIONS'
1715 #ifdef MPI
1716       include "mpif.h"
1717 #endif
1718       include 'COMMON.CONTROL'
1719       include 'COMMON.IOUNITS'
1720       include 'COMMON.GEO'
1721       include 'COMMON.VAR'
1722       include 'COMMON.LOCAL'
1723       include 'COMMON.CHAIN'
1724       include 'COMMON.DERIV'
1725       include 'COMMON.INTERACT'
1726 #ifdef FOURBODY
1727       include 'COMMON.CONTACTS'
1728       include 'COMMON.CONTMAT'
1729 #endif
1730       include 'COMMON.CORRMAT'
1731       include 'COMMON.TORSION'
1732       include 'COMMON.VECTORS'
1733       include 'COMMON.FFIELD'
1734       include 'COMMON.TIME1'
1735       include 'COMMON.SHIELD'
1736       include "COMMON.SPLITELE"
1737       integer xshift,yshift,zshift
1738       double precision ggg(3),gggp(3),gggm(3),erij(3),dcosb(3),dcosg(3),
1739      &          erder(3,3),uryg(3,3),urzg(3,3),vryg(3,3),vrzg(3,3)
1740       double precision acipa(2,2),agg(3,4),aggi(3,4),aggi1(3,4),
1741      &    aggj(3,4),aggj1(3,4),a_temp(2,2),muij(4),gmuij1(4),gmuji1(4),
1742      &    gmuij2(4),gmuji2(4)
1743       integer j1,j2,num_conti
1744       common /locel/ a_temp,agg,aggi,aggi1,aggj,aggj1,a22,a23,a32,a33,
1745      &    dxi,dyi,dzi,dx_normi,dy_normi,dz_normi,xmedi,ymedi,zmedi,
1746      &    num_conti,j1,j2
1747       integer k,i,j,iteli,itelj,kkk,l,kkll,m,isubchap,ind,itypi,itypj
1748       integer ilist,iresshield
1749       double precision rlocshield
1750       double precision ael6i,rrmij,rmij,r0ij,fcont,fprimcont,ees0tmp
1751       double precision ees,evdw1,eel_loc,aaa,bbb,ael3i
1752       double precision dxj,dyj,dzj,dx_normj,dy_normj,dz_normj,xj,yj,zj,
1753      &  rij,r3ij,r6ij,cosa,cosb,cosg,fac,ev1,ev2,fac3,fac4,
1754      &  evdwij,el1,el2,eesij,ees0ij,facvdw,facel,fac1,ecosa,
1755      &  ecosb,ecosg,ury,urz,vry,vrz,facr,a22der,a23der,a32der,
1756      &  a33der,eel_loc_ij,cosa4,wij,cosbg1,cosbg2,ees0pij,
1757      &  ees0pij1,ees0mij,ees0mij1,fac3p,ees0mijp,ees0pijp,
1758      &  ecosa1,ecosb1,ecosg1,ecosa2,ecosb2,ecosg2,ecosap,ecosbp,
1759      &  ecosgp,ecosam,ecosbm,ecosgm,ghalf,geel_loc_ij,geel_loc_ji,
1760      &  dxi,dyi,dzi,a22,a23,a32,a33
1761       double precision dist_init,xmedi,ymedi,zmedi,xj_safe,yj_safe,
1762      &  zj_safe,xj_temp,yj_temp,zj_temp,dist_temp,dx_normi,dy_normi,
1763      &  dz_normi,aux
1764       double precision sss1,sssgrad1
1765       double precision sscale,sscagrad
1766       double precision scalar
1767
1768 c 4/26/02 - AL scaling factor for 1,4 repulsive VDW interactions
1769 #ifdef MOMENT
1770       double precision scal_el /1.0d0/
1771 #else
1772       double precision scal_el /0.5d0/
1773 #endif
1774 C 12/13/98
1775 C 13-go grudnia roku pamietnego...
1776       double precision unmat(3,3) /1.0d0,0.0d0,0.0d0,
1777      &                   0.0d0,1.0d0,0.0d0,
1778      &                   0.0d0,0.0d0,1.0d0/
1779 c          time00=MPI_Wtime()
1780 cd      write (iout,*) "eelecij",i,j
1781 C      print *,"WCHODZE2"
1782       ind=ind+1
1783       iteli=itel(i)
1784       itelj=itel(j)
1785       if (j.eq.i+2 .and. itelj.eq.2) iteli=2
1786       aaa=app(iteli,itelj)
1787       bbb=bpp(iteli,itelj)
1788       ael6i=ael6(iteli,itelj)
1789       ael3i=ael3(iteli,itelj)
1790       dxj=dc(1,j)
1791       dyj=dc(2,j)
1792       dzj=dc(3,j)
1793       dx_normj=dc_norm(1,j)
1794       dy_normj=dc_norm(2,j)
1795       dz_normj=dc_norm(3,j)
1796       xj=c(1,j)+0.5D0*dxj
1797       yj=c(2,j)+0.5D0*dyj
1798       zj=c(3,j)+0.5D0*dzj
1799       xj=mod(xj,boxxsize)
1800       if (xj.lt.0) xj=xj+boxxsize
1801       yj=mod(yj,boxysize)
1802       if (yj.lt.0) yj=yj+boxysize
1803       zj=mod(zj,boxzsize)
1804       if (zj.lt.0) zj=zj+boxzsize
1805       dist_init=(xj-xmedi)**2+(yj-ymedi)**2+(zj-zmedi)**2
1806       xj_safe=xj
1807       yj_safe=yj
1808       zj_safe=zj
1809       isubchap=0
1810       do xshift=-1,1
1811       do yshift=-1,1
1812       do zshift=-1,1
1813           xj=xj_safe+xshift*boxxsize
1814           yj=yj_safe+yshift*boxysize
1815           zj=zj_safe+zshift*boxzsize
1816           dist_temp=(xj-xmedi)**2+(yj-ymedi)**2+(zj-zmedi)**2
1817           if(dist_temp.lt.dist_init) then
1818             dist_init=dist_temp
1819             xj_temp=xj
1820             yj_temp=yj
1821             zj_temp=zj
1822             isubchap=1
1823           endif
1824       enddo
1825       enddo
1826       enddo
1827       if (isubchap.eq.1) then
1828          xj=xj_temp-xmedi
1829          yj=yj_temp-ymedi
1830          zj=zj_temp-zmedi
1831       else
1832          xj=xj_safe-xmedi
1833          yj=yj_safe-ymedi
1834          zj=zj_safe-zmedi
1835       endif
1836
1837       rij=xj*xj+yj*yj+zj*zj
1838       rrmij=1.0D0/rij
1839       rij=dsqrt(rij)
1840       rmij=1.0D0/rij
1841 c For extracting the short-range part of Evdwpp
1842       sss1=sscale(rij,r_cut_int)
1843       if (sss1.eq.0.0d0) return
1844       sss=sscale(rij/rpp(iteli,itelj),r_cut_respa)
1845       sssgrad=sscagrad(rij/rpp(iteli,itelj),r_cut_respa)
1846       sssgrad1=sscagrad(rij,r_cut_int)
1847       r3ij=rrmij*rmij
1848       r6ij=r3ij*r3ij
1849       cosa=dx_normi*dx_normj+dy_normi*dy_normj+dz_normi*dz_normj
1850       cosb=(xj*dx_normi+yj*dy_normi+zj*dz_normi)*rmij
1851       cosg=(xj*dx_normj+yj*dy_normj+zj*dz_normj)*rmij
1852       fac=cosa-3.0D0*cosb*cosg
1853       ev1=aaa*r6ij*r6ij
1854 c 4/26/02 - AL scaling down 1,4 repulsive VDW interactions
1855       if (j.eq.i+2) ev1=scal_el*ev1
1856       ev2=bbb*r6ij
1857       fac3=ael6i*r6ij
1858       fac4=ael3i*r3ij
1859       evdwij=ev1+ev2
1860       if (shield_mode.eq.0) then
1861         fac_shield(i)=1.0
1862         fac_shield(j)=1.0
1863       endif
1864       el1=fac3*(4.0D0+fac*fac-3.0D0*(cosb*cosb+cosg*cosg))
1865       el2=fac4*fac
1866       el1=el1*fac_shield(i)**2*fac_shield(j)**2
1867       el2=el2*fac_shield(i)**2*fac_shield(j)**2
1868       eesij=el1+el2
1869 C 12/26/95 - for the evaluation of multi-body H-bonding interactions
1870       ees0ij=4.0D0+fac*fac-3.0D0*(cosb*cosb+cosg*cosg)
1871       ees=ees+eesij*sss1
1872       evdw1=evdw1+evdwij*(1.0d0-sss)*sss1
1873 cd          write(iout,'(2(2i3,2x),7(1pd12.4)/2(3(1pd12.4),5x)/)')
1874 cd     &      iteli,i,itelj,j,aaa,bbb,ael6i,ael3i,
1875 cd     &      1.0D0/dsqrt(rrmij),evdwij,eesij,
1876 cd     &      xmedi,ymedi,zmedi,xj,yj,zj
1877
1878       if (energy_dec) then
1879           write (iout,'(a6,2i5,0pf7.3,2f7.3)')
1880      &              'evdw1',i,j,evdwij,sss,sss1
1881           write (iout,'(a6,2i5,0pf7.3)') 'ees',i,j,eesij
1882       endif
1883
1884 C
1885 C Calculate contributions to the Cartesian gradient.
1886 C
1887 #ifdef SPLITELE
1888       facvdw=-6*rrmij*(ev1+evdwij)*(1.0d0-sss)*sss1
1889 c     &  *((sslipi+sslipj)/2.0d0*lipscale**2+1.0d0)
1890       facel=-3*rrmij*(el1+eesij)*sss1
1891 c     &  *((sslipi+sslipj)/2.0d0*lipscale**2+1.0d0)
1892       fac1=fac
1893       erij(1)=xj*rmij
1894       erij(2)=yj*rmij
1895       erij(3)=zj*rmij
1896 *
1897 * Radial derivatives. First process both termini of the fragment (i,j)
1898 *
1899       aux=facel+sssgrad1*(1.0d0-sss)*eesij*rmij
1900 c     & *((sslipi+sslipj)/2.0d0*lipscale**2+1.0d0)
1901       ggg(1)=aux*xj
1902       ggg(2)=aux*yj
1903       ggg(3)=aux*zj
1904 c      ggg(1)=facel*xj
1905 c      ggg(2)=facel*yj
1906 c      ggg(3)=facel*zj
1907       if ((fac_shield(i).gt.0).and.(fac_shield(j).gt.0).and.
1908      &  (shield_mode.gt.0)) then
1909 C          print *,i,j
1910         do ilist=1,ishield_list(i)
1911           iresshield=shield_list(ilist,i)
1912           do k=1,3
1913            rlocshield=grad_shield_side(k,ilist,i)*eesij*sss1
1914      &      /fac_shield(i)*2.0*sss1
1915            gshieldx(k,iresshield)=gshieldx(k,iresshield)+
1916      &              rlocshield
1917      &     +grad_shield_loc(k,ilist,i)*eesij*sss1/fac_shield(i)*2.0
1918      &      *sss1
1919             gshieldc(k,iresshield-1)=gshieldc(k,iresshield-1)+rlocshield
1920 C           gshieldc_loc(k,iresshield)=gshieldc_loc(k,iresshield)
1921 C     & +grad_shield_loc(k,ilist,i)*eesij/fac_shield(i)
1922 C             if (iresshield.gt.i) then
1923 C               do ishi=i+1,iresshield-1
1924 C                gshieldc(k,ishi)=gshieldc(k,ishi)+rlocshield
1925 C     & +grad_shield_loc(k,ilist,i)*eesij/fac_shield(i)
1926 C
1927 C              enddo
1928 C             else
1929 C               do ishi=iresshield,i
1930 C                gshieldc(k,ishi)=gshieldc(k,ishi)-rlocshield
1931 C     & -grad_shield_loc(k,ilist,i)*eesij/fac_shield(i)
1932 C
1933 C               enddo
1934 C              endif
1935           enddo
1936         enddo
1937         do ilist=1,ishield_list(j)
1938           iresshield=shield_list(ilist,j)
1939           do k=1,3
1940            rlocshield=grad_shield_side(k,ilist,j)*eesij/fac_shield(j)
1941      &     *2.0*sss1
1942            gshieldx(k,iresshield)=gshieldx(k,iresshield)+
1943      &        rlocshield
1944      &        +grad_shield_loc(k,ilist,j)*eesij/fac_shield(j)*2.0*sss1
1945            gshieldc(k,iresshield-1)=gshieldc(k,iresshield-1)+rlocshield
1946           enddo
1947         enddo
1948
1949         do k=1,3
1950            gshieldc(k,i)=gshieldc(k,i)+
1951      &             grad_shield(k,i)*eesij/fac_shield(i)*2.0*sss1
1952            gshieldc(k,j)=gshieldc(k,j)+
1953      &             grad_shield(k,j)*eesij/fac_shield(j)*2.0*sss1
1954            gshieldc(k,i-1)=gshieldc(k,i-1)+
1955      &             grad_shield(k,i)*eesij/fac_shield(i)*2.0*sss1
1956            gshieldc(k,j-1)=gshieldc(k,j-1)+
1957      &             grad_shield(k,j)*eesij/fac_shield(j)*2.0*sss1
1958
1959         enddo
1960       endif
1961
1962 c          do k=1,3
1963 c            ghalf=0.5D0*ggg(k)
1964 c            gelc(k,i)=gelc(k,i)+ghalf
1965 c            gelc(k,j)=gelc(k,j)+ghalf
1966 c          enddo
1967 c 9/28/08 AL Gradient compotents will be summed only at the end
1968       do k=1,3
1969         gelc_long(k,j)=gelc_long(k,j)+ggg(k)
1970         gelc_long(k,i)=gelc_long(k,i)-ggg(k)
1971       enddo
1972 c      gelc_long(3,i)=gelc_long(3,i)+
1973 c        ssgradlipi*eesij/2.0d0*lipscale**2*sss1
1974 *
1975 * Loop over residues i+1 thru j-1.
1976 *
1977 cgrad          do k=i+1,j-1
1978 cgrad            do l=1,3
1979 cgrad              gelc(l,k)=gelc(l,k)+ggg(l)
1980 cgrad            enddo
1981 cgrad          enddo
1982       facvdw=facvdw+
1983      & (-sss1*sssgrad/rpp(iteli,itelj)+(1.0d0-sss)*sssgrad1)*rmij*evdwij
1984 c     &   *((sslipi+sslipj)/2.0d0*lipscale**2+1.0d0)
1985       ggg(1)=facvdw*xj
1986       ggg(2)=facvdw*yj
1987       ggg(3)=facvdw*zj
1988 c          do k=1,3
1989 c            ghalf=0.5D0*ggg(k)
1990 c            gvdwpp(k,i)=gvdwpp(k,i)+ghalf
1991 c            gvdwpp(k,j)=gvdwpp(k,j)+ghalf
1992 c          enddo
1993 c 9/28/08 AL Gradient compotents will be summed only at the end
1994       do k=1,3
1995         gvdwpp(k,j)=gvdwpp(k,j)+ggg(k)
1996         gvdwpp(k,i)=gvdwpp(k,i)-ggg(k)
1997       enddo
1998 *
1999 * Loop over residues i+1 thru j-1.
2000 *
2001 cgrad          do k=i+1,j-1
2002 cgrad            do l=1,3
2003 cgrad              gvdwpp(l,k)=gvdwpp(l,k)+ggg(l)
2004 cgrad            enddo
2005 cgrad          enddo
2006 #else
2007       facvdw=-6*rrmij*(ev1+evdwij)*(1.0d0-sss)*sss1
2008 c     &  *((sslipi+sslipj)/2.0d0*lipscale**2+1.0d0)
2009       facel=-3*rrmij*(el1+eesij)*sss1
2010 c     &  *((sslipi+sslipj)/2.0d0*lipscale**2+1.0d0)
2011
2012 c      facvdw=ev1+evdwij*(1.0d0-sss)*sss1
2013 c      facel=el1+eesij
2014       fac1=fac
2015       fac=-3*rrmij*(facvdw+facvdw+facel)
2016       erij(1)=xj*rmij
2017       erij(2)=yj*rmij
2018       erij(3)=zj*rmij
2019 *
2020 * Radial derivatives. First process both termini of the fragment (i,j)
2021 *
2022       aux=fac+(sssgrad1*(1.0d0-sss)-sssgrad*sss1/rpp(iteli,itelj))
2023      &  *eesij*rmij
2024 c     & *((sslipi+sslipj)/2.0d0*lipscale**2+1.0d0)
2025       ggg(1)=aux*xj
2026       ggg(2)=aux*yj
2027       ggg(3)=axu*zj
2028 c      ggg(1)=fac*xj
2029 c      ggg(2)=fac*yj
2030 c      ggg(3)=fac*zj
2031 c          do k=1,3
2032 c            ghalf=0.5D0*ggg(k)
2033 c            gelc(k,i)=gelc(k,i)+ghalf
2034 c            gelc(k,j)=gelc(k,j)+ghalf
2035 c          enddo
2036 c 9/28/08 AL Gradient compotents will be summed only at the end
2037       do k=1,3
2038         gelc_long(k,j)=gelc(k,j)+ggg(k)
2039         gelc_long(k,i)=gelc(k,i)-ggg(k)
2040       enddo
2041 *
2042 * Loop over residues i+1 thru j-1.
2043 *
2044 cgrad          do k=i+1,j-1
2045 cgrad            do l=1,3
2046 cgrad              gelc(l,k)=gelc(l,k)+ggg(l)
2047 cgrad            enddo
2048 cgrad          enddo
2049 c 9/28/08 AL Gradient compotents will be summed only at the end
2050 C          ggg(1)=facvdw*xj
2051 C          ggg(2)=facvdw*yj
2052 C          ggg(3)=facvdw*zj
2053       facvdw=facvdw
2054      & (-sssgrad*sss1/rpp(iteli,itelj)+sssgrad1*(1.0d0-sss))*rmij*evdwij
2055       ggg(1)=facvdw*xj
2056       ggg(2)=facvdw*yj
2057       ggg(3)=facvdw*zj
2058       do k=1,3
2059         gvdwpp(k,j)=gvdwpp(k,j)+ggg(k)
2060         gvdwpp(k,i)=gvdwpp(k,i)-ggg(k)
2061       enddo
2062 #endif
2063 *
2064 * Angular part
2065 *
2066       ecosa=2.0D0*fac3*fac1+fac4
2067       fac4=-3.0D0*fac4
2068       fac3=-6.0D0*fac3
2069       ecosb=(fac3*(fac1*cosg+cosb)+cosg*fac4)
2070       ecosg=(fac3*(fac1*cosb+cosg)+cosb*fac4)
2071       do k=1,3
2072         dcosb(k)=rmij*(dc_norm(k,i)-erij(k)*cosb)
2073         dcosg(k)=rmij*(dc_norm(k,j)-erij(k)*cosg)
2074       enddo
2075 cd        print '(2i3,2(3(1pd14.5),3x))',i,j,(dcosb(k),k=1,3),
2076 cd   &          (dcosg(k),k=1,3)
2077       do k=1,3
2078         ggg(k)=(ecosb*dcosb(k)+ecosg*dcosg(k))*sss1
2079      &  *fac_shield(i)**2*fac_shield(j)**2
2080 c     &  *((sslipi+sslipj)/2.0d0*lipscale**2+1.0d0)
2081
2082       enddo
2083 c          do k=1,3
2084 c            ghalf=0.5D0*ggg(k)
2085 c            gelc(k,i)=gelc(k,i)+ghalf
2086 c     &               +(ecosa*(dc_norm(k,j)-cosa*dc_norm(k,i))
2087 c     &               + ecosb*(erij(k)-cosb*dc_norm(k,i)))*vbld_inv(i+1)
2088 c            gelc(k,j)=gelc(k,j)+ghalf
2089 c     &               +(ecosa*(dc_norm(k,i)-cosa*dc_norm(k,j))
2090 c     &               + ecosg*(erij(k)-cosg*dc_norm(k,j)))*vbld_inv(j+1)
2091 c          enddo
2092 cgrad          do k=i+1,j-1
2093 cgrad            do l=1,3
2094 cgrad              gelc(l,k)=gelc(l,k)+ggg(l)
2095 cgrad            enddo
2096 cgrad          enddo
2097       do k=1,3
2098         gelc(k,i)=gelc(k,i)
2099      &           +((ecosa*(dc_norm(k,j)-cosa*dc_norm(k,i))
2100      &           +ecosb*(erij(k)-cosb*dc_norm(k,i)))*vbld_inv(i+1))*sss1
2101      &       *fac_shield(i)**2*fac_shield(j)**2
2102 c     &       *((sslipi+sslipj)/2.0d0*lipscale**2+1.0d0)
2103
2104         gelc(k,j)=gelc(k,j)
2105      &           +((ecosa*(dc_norm(k,i)-cosa*dc_norm(k,j))
2106      &           +ecosg*(erij(k)-cosg*dc_norm(k,j)))*vbld_inv(j+1))*sss1
2107      &       *fac_shield(i)**2*fac_shield(j)**2
2108 c     &       *((sslipi+sslipj)/2.0d0*lipscale**2+1.0d0)
2109         gelc_long(k,j)=gelc_long(k,j)+ggg(k)
2110         gelc_long(k,i)=gelc_long(k,i)-ggg(k)
2111       enddo
2112       IF (wel_loc.gt.0.0d0 .or. wcorr4.gt.0.0d0 .or. wcorr5.gt.0.0d0
2113      &    .or. wcorr6.gt.0.0d0 .or. wturn3.gt.0.0d0
2114      &    .or. wturn4.gt.0.0d0 .or. wturn6.gt.0.0d0) THEN
2115 C
2116 C 9/25/99 Mixed third-order local-electrostatic terms. The local-interaction
2117 C   energy of a peptide unit is assumed in the form of a second-order
2118 C   Fourier series in the angles lambda1 and lambda2 (see Nishikawa et al.
2119 C   Macromolecules, 1974, 7, 797-806 for definition). This correlation terms
2120 C   are computed for EVERY pair of non-contiguous peptide groups.
2121 C
2122         if (j.lt.nres-1) then
2123           j1=j+1
2124           j2=j-1
2125         else
2126           j1=j-1
2127           j2=j-2
2128         endif
2129         kkk=0
2130         do k=1,2
2131           do l=1,2
2132             kkk=kkk+1
2133             muij(kkk)=mu(k,i)*mu(l,j)
2134 #ifdef NEWCORR
2135             gmuij1(kkk)=gtb1(k,i+1)*mu(l,j)
2136 c             write(iout,*) 'k=',k,i,gtb1(k,i+1),gtb1(k,i+1)*mu(l,j)
2137             gmuij2(kkk)=gUb2(k,i)*mu(l,j)
2138             gmuji1(kkk)=mu(k,i)*gtb1(l,j+1)
2139 c             write(iout,*) 'l=',l,j,gtb1(l,j+1),gtb1(l,j+1)*mu(k,i)
2140             gmuji2(kkk)=mu(k,i)*gUb2(l,j)
2141 #endif
2142           enddo
2143         enddo
2144 cd         write (iout,*) 'EELEC: i',i,' j',j
2145 cd          write (iout,*) 'j',j,' j1',j1,' j2',j2
2146 cd          write(iout,*) 'muij',muij
2147         ury=scalar(uy(1,i),erij)
2148         urz=scalar(uz(1,i),erij)
2149         vry=scalar(uy(1,j),erij)
2150         vrz=scalar(uz(1,j),erij)
2151         a22=scalar(uy(1,i),uy(1,j))-3*ury*vry
2152         a23=scalar(uy(1,i),uz(1,j))-3*ury*vrz
2153         a32=scalar(uz(1,i),uy(1,j))-3*urz*vry
2154         a33=scalar(uz(1,i),uz(1,j))-3*urz*vrz
2155         fac=dsqrt(-ael6i)*r3ij
2156         a22=a22*fac
2157         a23=a23*fac
2158         a32=a32*fac
2159         a33=a33*fac
2160 cd          write (iout,'(4i5,4f10.5)')
2161 cd     &     i,itortyp(itype(i)),j,itortyp(itype(j)),a22,a23,a32,a33
2162 cd          write (iout,'(6f10.5)') (muij(k),k=1,4),fac,eel_loc_ij
2163 cd          write (iout,'(2(3f10.5,5x)/2(3f10.5,5x))') uy(:,i),uz(:,i),
2164 cd     &      uy(:,j),uz(:,j)
2165 cd          write (iout,'(4f10.5)')
2166 cd     &      scalar(uy(1,i),uy(1,j)),scalar(uy(1,i),uz(1,j)),
2167 cd     &      scalar(uz(1,i),uy(1,j)),scalar(uz(1,i),uz(1,j))
2168 cd          write (iout,'(4f10.5)') ury,urz,vry,vrz
2169 cd           write (iout,'(9f10.5/)')
2170 cd     &      fac22,a22,fac23,a23,fac32,a32,fac33,a33,eel_loc_ij
2171 C Derivatives of the elements of A in virtual-bond vectors
2172         call unormderiv(erij(1),unmat(1,1),rmij,erder(1,1))
2173         do k=1,3
2174           uryg(k,1)=scalar(erder(1,k),uy(1,i))
2175           uryg(k,2)=scalar(uygrad(1,k,1,i),erij(1))
2176           uryg(k,3)=scalar(uygrad(1,k,2,i),erij(1))
2177           urzg(k,1)=scalar(erder(1,k),uz(1,i))
2178           urzg(k,2)=scalar(uzgrad(1,k,1,i),erij(1))
2179           urzg(k,3)=scalar(uzgrad(1,k,2,i),erij(1))
2180           vryg(k,1)=scalar(erder(1,k),uy(1,j))
2181           vryg(k,2)=scalar(uygrad(1,k,1,j),erij(1))
2182           vryg(k,3)=scalar(uygrad(1,k,2,j),erij(1))
2183           vrzg(k,1)=scalar(erder(1,k),uz(1,j))
2184           vrzg(k,2)=scalar(uzgrad(1,k,1,j),erij(1))
2185           vrzg(k,3)=scalar(uzgrad(1,k,2,j),erij(1))
2186         enddo
2187 C Compute radial contributions to the gradient
2188         facr=-3.0d0*rrmij
2189         a22der=a22*facr
2190         a23der=a23*facr
2191         a32der=a32*facr
2192         a33der=a33*facr
2193         agg(1,1)=a22der*xj
2194         agg(2,1)=a22der*yj
2195         agg(3,1)=a22der*zj
2196         agg(1,2)=a23der*xj
2197         agg(2,2)=a23der*yj
2198         agg(3,2)=a23der*zj
2199         agg(1,3)=a32der*xj
2200         agg(2,3)=a32der*yj
2201         agg(3,3)=a32der*zj
2202         agg(1,4)=a33der*xj
2203         agg(2,4)=a33der*yj
2204         agg(3,4)=a33der*zj
2205 C Add the contributions coming from er
2206         fac3=-3.0d0*fac
2207         do k=1,3
2208           agg(k,1)=agg(k,1)+fac3*(uryg(k,1)*vry+vryg(k,1)*ury)
2209           agg(k,2)=agg(k,2)+fac3*(uryg(k,1)*vrz+vrzg(k,1)*ury)
2210           agg(k,3)=agg(k,3)+fac3*(urzg(k,1)*vry+vryg(k,1)*urz)
2211           agg(k,4)=agg(k,4)+fac3*(urzg(k,1)*vrz+vrzg(k,1)*urz)
2212         enddo
2213         do k=1,3
2214 C Derivatives in DC(i)
2215 cgrad            ghalf1=0.5d0*agg(k,1)
2216 cgrad            ghalf2=0.5d0*agg(k,2)
2217 cgrad            ghalf3=0.5d0*agg(k,3)
2218 cgrad            ghalf4=0.5d0*agg(k,4)
2219           aggi(k,1)=fac*(scalar(uygrad(1,k,1,i),uy(1,j))
2220      &    -3.0d0*uryg(k,2)*vry)!+ghalf1
2221           aggi(k,2)=fac*(scalar(uygrad(1,k,1,i),uz(1,j))
2222      &    -3.0d0*uryg(k,2)*vrz)!+ghalf2
2223           aggi(k,3)=fac*(scalar(uzgrad(1,k,1,i),uy(1,j))
2224      &    -3.0d0*urzg(k,2)*vry)!+ghalf3
2225           aggi(k,4)=fac*(scalar(uzgrad(1,k,1,i),uz(1,j))
2226      &    -3.0d0*urzg(k,2)*vrz)!+ghalf4
2227 C Derivatives in DC(i+1)
2228           aggi1(k,1)=fac*(scalar(uygrad(1,k,2,i),uy(1,j))
2229      &    -3.0d0*uryg(k,3)*vry)!+agg(k,1)
2230           aggi1(k,2)=fac*(scalar(uygrad(1,k,2,i),uz(1,j))
2231      &    -3.0d0*uryg(k,3)*vrz)!+agg(k,2)
2232           aggi1(k,3)=fac*(scalar(uzgrad(1,k,2,i),uy(1,j))
2233      &    -3.0d0*urzg(k,3)*vry)!+agg(k,3)
2234           aggi1(k,4)=fac*(scalar(uzgrad(1,k,2,i),uz(1,j))
2235      &    -3.0d0*urzg(k,3)*vrz)!+agg(k,4)
2236 C Derivatives in DC(j)
2237           aggj(k,1)=fac*(scalar(uygrad(1,k,1,j),uy(1,i))
2238      &    -3.0d0*vryg(k,2)*ury)!+ghalf1
2239           aggj(k,2)=fac*(scalar(uzgrad(1,k,1,j),uy(1,i))
2240      &    -3.0d0*vrzg(k,2)*ury)!+ghalf2
2241           aggj(k,3)=fac*(scalar(uygrad(1,k,1,j),uz(1,i))
2242      &    -3.0d0*vryg(k,2)*urz)!+ghalf3
2243           aggj(k,4)=fac*(scalar(uzgrad(1,k,1,j),uz(1,i))
2244      &    -3.0d0*vrzg(k,2)*urz)!+ghalf4
2245 C Derivatives in DC(j+1) or DC(nres-1)
2246           aggj1(k,1)=fac*(scalar(uygrad(1,k,2,j),uy(1,i))
2247      &    -3.0d0*vryg(k,3)*ury)
2248           aggj1(k,2)=fac*(scalar(uzgrad(1,k,2,j),uy(1,i))
2249      &    -3.0d0*vrzg(k,3)*ury)
2250           aggj1(k,3)=fac*(scalar(uygrad(1,k,2,j),uz(1,i))
2251      &    -3.0d0*vryg(k,3)*urz)
2252           aggj1(k,4)=fac*(scalar(uzgrad(1,k,2,j),uz(1,i))
2253      &    -3.0d0*vrzg(k,3)*urz)
2254 cgrad            if (j.eq.nres-1 .and. i.lt.j-2) then
2255 cgrad              do l=1,4
2256 cgrad                aggj1(k,l)=aggj1(k,l)+agg(k,l)
2257 cgrad              enddo
2258 cgrad            endif
2259         enddo
2260         acipa(1,1)=a22
2261         acipa(1,2)=a23
2262         acipa(2,1)=a32
2263         acipa(2,2)=a33
2264         a22=-a22
2265         a23=-a23
2266         do l=1,2
2267           do k=1,3
2268             agg(k,l)=-agg(k,l)
2269             aggi(k,l)=-aggi(k,l)
2270             aggi1(k,l)=-aggi1(k,l)
2271             aggj(k,l)=-aggj(k,l)
2272             aggj1(k,l)=-aggj1(k,l)
2273           enddo
2274         enddo
2275         if (j.lt.nres-1) then
2276           a22=-a22
2277           a32=-a32
2278           do l=1,3,2
2279             do k=1,3
2280               agg(k,l)=-agg(k,l)
2281               aggi(k,l)=-aggi(k,l)
2282               aggi1(k,l)=-aggi1(k,l)
2283               aggj(k,l)=-aggj(k,l)
2284               aggj1(k,l)=-aggj1(k,l)
2285             enddo
2286           enddo
2287         else
2288           a22=-a22
2289           a23=-a23
2290           a32=-a32
2291           a33=-a33
2292           do l=1,4
2293             do k=1,3
2294               agg(k,l)=-agg(k,l)
2295               aggi(k,l)=-aggi(k,l)
2296               aggi1(k,l)=-aggi1(k,l)
2297               aggj(k,l)=-aggj(k,l)
2298               aggj1(k,l)=-aggj1(k,l)
2299             enddo
2300           enddo
2301         endif
2302       ENDIF ! WCORR
2303       IF (wel_loc.gt.0.0d0) THEN
2304 C Contribution to the local-electrostatic energy coming from the i-j pair
2305         eel_loc_ij=a22*muij(1)+a23*muij(2)+a32*muij(3)
2306      &   +a33*muij(4)
2307 cd        write (iout,*) 'i',i,' j',j,' eel_loc_ij',eel_loc_ij
2308
2309         if (energy_dec) write (iout,'(a6,2i5,0pf7.3)')
2310      &          'eelloc',i,j,eel_loc_ij
2311
2312
2313         if (shield_mode.eq.0) then
2314           fac_shield(i)=1.0
2315           fac_shield(j)=1.0
2316 C        else
2317 C         fac_shield(i)=0.4
2318 C         fac_shield(j)=0.6
2319         endif
2320         eel_loc_ij=eel_loc_ij
2321      &  *fac_shield(i)*fac_shield(j)*sss1
2322         eel_loc=eel_loc+eel_loc_ij
2323
2324         if ((fac_shield(i).gt.0).and.(fac_shield(j).gt.0).and.
2325      &  (shield_mode.gt.0)) then
2326 C          print *,i,j
2327
2328           do ilist=1,ishield_list(i)
2329            iresshield=shield_list(ilist,i)
2330            do k=1,3
2331            rlocshield=grad_shield_side(k,ilist,i)*eel_loc_ij
2332      &                                          /fac_shield(i)
2333 C     &      *2.0
2334            gshieldx_ll(k,iresshield)=gshieldx_ll(k,iresshield)+
2335      &              rlocshield
2336      &      +grad_shield_loc(k,ilist,i)*eel_loc_ij/fac_shield(i)
2337             gshieldc_ll(k,iresshield-1)=gshieldc_ll(k,iresshield-1)
2338      &      +rlocshield
2339            enddo
2340           enddo
2341           do ilist=1,ishield_list(j)
2342            iresshield=shield_list(ilist,j)
2343            do k=1,3
2344            rlocshield=grad_shield_side(k,ilist,j)*eel_loc_ij
2345      &                                       /fac_shield(j)
2346 C     &     *2.0
2347            gshieldx_ll(k,iresshield)=gshieldx_ll(k,iresshield)+
2348      &              rlocshield
2349      &     +grad_shield_loc(k,ilist,j)*eel_loc_ij/fac_shield(j)
2350            gshieldc_ll(k,iresshield-1)=gshieldc_ll(k,iresshield-1)
2351      &             +rlocshield
2352
2353            enddo
2354           enddo
2355
2356           do k=1,3
2357             gshieldc_ll(k,i)=gshieldc_ll(k,i)+
2358      &              grad_shield(k,i)*eel_loc_ij/fac_shield(i)
2359             gshieldc_ll(k,j)=gshieldc_ll(k,j)+
2360      &              grad_shield(k,j)*eel_loc_ij/fac_shield(j)
2361             gshieldc_ll(k,i-1)=gshieldc_ll(k,i-1)+
2362      &              grad_shield(k,i)*eel_loc_ij/fac_shield(i)
2363             gshieldc_ll(k,j-1)=gshieldc_ll(k,j-1)+
2364      &              grad_shield(k,j)*eel_loc_ij/fac_shield(j)
2365           enddo
2366         endif
2367
2368 #ifdef NEWCORR
2369         geel_loc_ij=(a22*gmuij1(1)
2370      &     +a23*gmuij1(2)
2371      &     +a32*gmuij1(3)
2372      &     +a33*gmuij1(4))
2373      &    *fac_shield(i)*fac_shield(j)*sss1
2374 c         write(iout,*) "derivative over thatai"
2375 c         write(iout,*) a22*gmuij1(1), a23*gmuij1(2) ,a32*gmuij1(3),
2376 c     &   a33*gmuij1(4)
2377         gloc(nphi+i,icg)=gloc(nphi+i,icg)+
2378      &      geel_loc_ij*wel_loc
2379 c         write(iout,*) "derivative over thatai-1"
2380 c         write(iout,*) a22*gmuij2(1), a23*gmuij2(2) ,a32*gmuij2(3),
2381 c     &   a33*gmuij2(4)
2382         geel_loc_ij=
2383      &     a22*gmuij2(1)
2384      &     +a23*gmuij2(2)
2385      &     +a32*gmuij2(3)
2386      &     +a33*gmuij2(4)
2387         gloc(nphi+i-1,icg)=gloc(nphi+i-1,icg)+
2388      &      geel_loc_ij*wel_loc
2389      &    *fac_shield(i)*fac_shield(j)*sss1
2390
2391 c  Derivative over j residue
2392         geel_loc_ji=a22*gmuji1(1)
2393      &     +a23*gmuji1(2)
2394      &     +a32*gmuji1(3)
2395      &     +a33*gmuji1(4)
2396 c         write(iout,*) "derivative over thataj"
2397 c         write(iout,*) a22*gmuji1(1), a23*gmuji1(2) ,a32*gmuji1(3),
2398 c     &   a33*gmuji1(4)
2399
2400         gloc(nphi+j,icg)=gloc(nphi+j,icg)+
2401      &      geel_loc_ji*wel_loc
2402      &    *fac_shield(i)*fac_shield(j)*sss1
2403
2404         geel_loc_ji=
2405      &     +a22*gmuji2(1)
2406      &     +a23*gmuji2(2)
2407      &     +a32*gmuji2(3)
2408      &     +a33*gmuji2(4)
2409 c         write(iout,*) "derivative over thataj-1"
2410 c         write(iout,*) a22*gmuji2(1), a23*gmuji2(2) ,a32*gmuji2(3),
2411 c     &   a33*gmuji2(4)
2412         gloc(nphi+j-1,icg)=gloc(nphi+j-1,icg)+
2413      &    geel_loc_ji*wel_loc
2414      &    *fac_shield(i)*fac_shield(j)*sss1
2415 #endif
2416 cC Paral derivatives in virtual-bond dihedral angles gamma
2417         if (i.gt.1)
2418      &    gel_loc_loc(i-1)=gel_loc_loc(i-1)+
2419      &          (a22*muder(1,i)*mu(1,j)+a23*muder(1,i)*mu(2,j)
2420      &         +a32*muder(2,i)*mu(1,j)+a33*muder(2,i)*mu(2,j))
2421      &         *fac_shield(i)*fac_shield(j)*sss1
2422 c     &         *fac_shield(i)*fac_shield(j)
2423 c     &         *((sslipi+sslipj)/2.0d0*lipscale+1.0d0)
2424
2425
2426         gel_loc_loc(j-1)=gel_loc_loc(j-1)+
2427      &          (a22*mu(1,i)*muder(1,j)+a23*mu(1,i)*muder(2,j)
2428      &         +a32*mu(2,i)*muder(1,j)+a33*mu(2,i)*muder(2,j))
2429      &         *fac_shield(i)*fac_shield(j)*sss1
2430 c     &         *fac_shield(i)*fac_shield(j)
2431 c     &         *((sslipi+sslipj)/2.0d0*lipscale+1.0d0)
2432
2433 C Derivatives of eello in DC(i+1) thru DC(j-1) or DC(nres-2)
2434           aux=eel_loc_ij/sss1*sssgrad1*rmij
2435           ggg(1)=aux*xj
2436           ggg(2)=aux*yj
2437           ggg(3)=aux*zj
2438         do l=1,3
2439           ggg(l)=ggg(l)+(agg(l,1)*muij(1)+
2440      &       agg(l,2)*muij(2)+agg(l,3)*muij(3)+agg(l,4)*muij(4))
2441      &       *fac_shield(i)*fac_shield(j)*sss1
2442 c     &         *fac_shield(i)*fac_shield(j)
2443 c     &         *((sslipi+sslipj)/2.0d0*lipscale+1.0d0)
2444
2445           gel_loc_long(l,j)=gel_loc_long(l,j)+ggg(l)
2446           gel_loc_long(l,i)=gel_loc_long(l,i)-ggg(l)
2447 cgrad          ghalf=0.5d0*ggg(l)
2448 cgrad          gel_loc(l,i)=gel_loc(l,i)+ghalf
2449 cgrad          gel_loc(l,j)=gel_loc(l,j)+ghalf
2450         enddo
2451 cgrad          do k=i+1,j2
2452 cgrad            do l=1,3
2453 cgrad              gel_loc(l,k)=gel_loc(l,k)+ggg(l)
2454 cgrad            enddo
2455 cgrad          enddo
2456 C Remaining derivatives of eello
2457 c          gel_loc_long(3,j)=gel_loc_long(3,j)+ &
2458 c          ssgradlipj*eel_loc_ij/2.0d0*lipscale/  &
2459 c          ((sslipi+sslipj)/2.0d0*lipscale+1.0d0)*sss_ele_cut
2460 c
2461 c          gel_loc_long(3,i)=gel_loc_long(3,i)+ &
2462 c          ssgradlipi*eel_loc_ij/2.0d0*lipscale/  &
2463 c          ((sslipi+sslipj)/2.0d0*lipscale+1.0d0)*sss_ele_cut
2464
2465         do l=1,3
2466           gel_loc(l,i)=gel_loc(l,i)+(aggi(l,1)*muij(1)+
2467      &        aggi(l,2)*muij(2)+aggi(l,3)*muij(3)+aggi(l,4)*muij(4))
2468      &       *fac_shield(i)*fac_shield(j)*sss1
2469 c     &       *((sslipi+sslipj)/2.0d0*lipscale+1.0d0)
2470
2471           gel_loc(l,i+1)=gel_loc(l,i+1)+(aggi1(l,1)*muij(1)+
2472      &       aggi1(l,2)*muij(2)+aggi1(l,3)*muij(3)+aggi1(l,4)*muij(4))
2473      &      *fac_shield(i)*fac_shield(j)*sss1
2474 c     &       *((sslipi+sslipj)/2.0d0*lipscale+1.0d0)
2475
2476           gel_loc(l,j)=gel_loc(l,j)+(aggj(l,1)*muij(1)+
2477      &        aggj(l,2)*muij(2)+aggj(l,3)*muij(3)+aggj(l,4)*muij(4))
2478      &       *fac_shield(i)*fac_shield(j)*sss1
2479 c     &       *((sslipi+sslipj)/2.0d0*lipscale+1.0d0)
2480
2481           gel_loc(l,j1)=gel_loc(l,j1)+(aggj1(l,1)*muij(1)+
2482      &       aggj1(l,2)*muij(2)+aggj1(l,3)*muij(3)+aggj1(l,4)*muij(4))
2483      &      *fac_shield(i)*fac_shield(j)*sss1
2484 c     &       *((sslipi+sslipj)/2.0d0*lipscale+1.0d0)
2485
2486         enddo
2487       ENDIF
2488 #ifdef FOURBODY
2489 C Change 12/26/95 to calculate four-body contributions to H-bonding energy
2490 c          if (j.gt.i+1 .and. num_conti.le.maxconts) then
2491       if (wcorr+wcorr4+wcorr5+wcorr6.gt.0.0d0
2492      &   .and. num_conti.le.maxconts) then
2493 c            write (iout,*) i,j," entered corr"
2494 C
2495 C Calculate the contact function. The ith column of the array JCONT will
2496 C contain the numbers of atoms that make contacts with the atom I (of numbers
2497 C greater than I). The arrays FACONT and GACONT will contain the values of
2498 C the contact function and its derivative.
2499 c       r0ij=1.02D0*rpp(iteli,itelj)
2500 c       r0ij=1.11D0*rpp(iteli,itelj)
2501         r0ij=2.20D0*rpp(iteli,itelj)
2502 c       r0ij=1.55D0*rpp(iteli,itelj)
2503         call gcont(rij,r0ij,1.0D0,0.2d0*r0ij,fcont,fprimcont)
2504         if (fcont.gt.0.0D0) then
2505           num_conti=num_conti+1
2506           if (num_conti.gt.maxconts) then
2507             write (iout,*) 'WARNING - max. # of contacts exceeded;',
2508      &                     ' will skip next contacts for this conf.'
2509           else
2510             jcont_hb(num_conti,i)=j
2511 cd         write (iout,*) "i",i," j",j," num_conti",num_conti,
2512 cd          " jcont_hb",jcont_hb(num_conti,i)
2513             IF (wcorr4.gt.0.0d0 .or. wcorr5.gt.0.0d0 .or.
2514      &      wcorr6.gt.0.0d0 .or. wturn6.gt.0.0d0) THEN
2515 C 9/30/99 (AL) - store components necessary to evaluate higher-order loc-el
2516 C  terms.
2517               d_cont(num_conti,i)=rij
2518 cd                write (2,'(3e15.5)') rij,r0ij+0.2d0*r0ij,rij
2519 C     --- Electrostatic-interaction matrix ---
2520               a_chuj(1,1,num_conti,i)=a22
2521               a_chuj(1,2,num_conti,i)=a23
2522               a_chuj(2,1,num_conti,i)=a32
2523               a_chuj(2,2,num_conti,i)=a33
2524 C     --- Gradient of rij
2525               do kkk=1,3
2526                 grij_hb_cont(kkk,num_conti,i)=erij(kkk)
2527               enddo
2528               kkll=0
2529               do k=1,2
2530                 do l=1,2
2531                   kkll=kkll+1
2532                   do m=1,3
2533                     a_chuj_der(k,l,m,1,num_conti,i)=agg(m,kkll)
2534                     a_chuj_der(k,l,m,2,num_conti,i)=aggi(m,kkll)
2535                     a_chuj_der(k,l,m,3,num_conti,i)=aggi1(m,kkll)
2536                     a_chuj_der(k,l,m,4,num_conti,i)=aggj(m,kkll)
2537                     a_chuj_der(k,l,m,5,num_conti,i)=aggj1(m,kkll)
2538                   enddo
2539                 enddo
2540               enddo
2541             ENDIF
2542             IF (wcorr4.eq.0.0d0 .and. wcorr.gt.0.0d0) THEN
2543 C Calculate contact energies
2544               cosa4=4.0D0*cosa
2545               wij=cosa-3.0D0*cosb*cosg
2546               cosbg1=cosb+cosg
2547               cosbg2=cosb-cosg
2548 c             fac3=dsqrt(-ael6i)/r0ij**3
2549               fac3=dsqrt(-ael6i)*r3ij
2550 c               ees0pij=dsqrt(4.0D0+cosa4+wij*wij-3.0D0*cosbg1*cosbg1)
2551               ees0tmp=4.0D0+cosa4+wij*wij-3.0D0*cosbg1*cosbg1
2552               if (ees0tmp.gt.0) then
2553                 ees0pij=dsqrt(ees0tmp)
2554               else
2555                 ees0pij=0
2556               endif
2557 c              ees0mij=dsqrt(4.0D0-cosa4+wij*wij-3.0D0*cosbg2*cosbg2)
2558               ees0tmp=4.0D0-cosa4+wij*wij-3.0D0*cosbg2*cosbg2
2559               if (ees0tmp.gt.0) then
2560                 ees0mij=dsqrt(ees0tmp)
2561               else
2562                 ees0mij=0
2563               endif
2564 c             ees0mij=0.0D0
2565               if (shield_mode.eq.0) then
2566                 fac_shield(i)=1.0d0
2567                 fac_shield(j)=1.0d0
2568               else
2569                 ees0plist(num_conti,i)=j
2570 C                fac_shield(i)=0.4d0
2571 C                fac_shield(j)=0.6d0
2572               endif
2573               ees0p(num_conti,i)=0.5D0*fac3*(ees0pij+ees0mij)
2574      &        *fac_shield(i)*fac_shield(j)*sss1
2575               ees0m(num_conti,i)=0.5D0*fac3*(ees0pij-ees0mij)
2576      &        *fac_shield(i)*fac_shield(j)*sss1
2577
2578 C Diagnostics. Comment out or remove after debugging!
2579 c               ees0p(num_conti,i)=0.5D0*fac3*ees0pij
2580 c               ees0m(num_conti,i)=0.5D0*fac3*ees0mij
2581 c               ees0m(num_conti,i)=0.0D0
2582 C End diagnostics.
2583 c               write (iout,*) 'i=',i,' j=',j,' rij=',rij,' r0ij=',r0ij,
2584 c    & ' ees0ij=',ees0p(num_conti,i),ees0m(num_conti,i),' fcont=',fcont
2585 C Angular derivatives of the contact function
2586               ees0pij1=fac3/ees0pij
2587               ees0mij1=fac3/ees0mij
2588               fac3p=-3.0D0*fac3*rrmij
2589               ees0pijp=0.5D0*fac3p*(ees0pij+ees0mij)
2590               ees0mijp=0.5D0*fac3p*(ees0pij-ees0mij)
2591 c             ees0mij1=0.0D0
2592               ecosa1=       ees0pij1*( 1.0D0+0.5D0*wij)
2593               ecosb1=-1.5D0*ees0pij1*(wij*cosg+cosbg1)
2594               ecosg1=-1.5D0*ees0pij1*(wij*cosb+cosbg1)
2595               ecosa2=       ees0mij1*(-1.0D0+0.5D0*wij)
2596               ecosb2=-1.5D0*ees0mij1*(wij*cosg+cosbg2)
2597               ecosg2=-1.5D0*ees0mij1*(wij*cosb-cosbg2)
2598               ecosap=ecosa1+ecosa2
2599               ecosbp=ecosb1+ecosb2
2600               ecosgp=ecosg1+ecosg2
2601               ecosam=ecosa1-ecosa2
2602               ecosbm=ecosb1-ecosb2
2603               ecosgm=ecosg1-ecosg2
2604 C Diagnostics
2605 c               ecosap=ecosa1
2606 c               ecosbp=ecosb1
2607 c               ecosgp=ecosg1
2608 c               ecosam=0.0D0
2609 c               ecosbm=0.0D0
2610 c               ecosgm=0.0D0
2611 C End diagnostics
2612               facont_hb(num_conti,i)=fcont
2613               fprimcont=fprimcont/rij
2614 cd              facont_hb(num_conti,i)=1.0D0
2615 C Following line is for diagnostics.
2616 cd              fprimcont=0.0D0
2617               do k=1,3
2618                 dcosb(k)=rmij*(dc_norm(k,i)-erij(k)*cosb)
2619                 dcosg(k)=rmij*(dc_norm(k,j)-erij(k)*cosg)
2620               enddo
2621               do k=1,3
2622                 gggp(k)=ecosbp*dcosb(k)+ecosgp*dcosg(k)
2623                 gggm(k)=ecosbm*dcosb(k)+ecosgm*dcosg(k)
2624               enddo
2625               gggp(1)=gggp(1)+ees0pijp*xj
2626      &          +ees0p(num_conti,i)/sss1*rmij*xj*sssgrad1
2627               gggp(2)=gggp(2)+ees0pijp*yj
2628      &          +ees0p(num_conti,i)/sss1*rmij*yj*sssgrad1
2629               gggp(3)=gggp(3)+ees0pijp*zj
2630      &          +ees0p(num_conti,i)/sss1*rmij*zj*sssgrad1
2631               gggm(1)=gggm(1)+ees0mijp*xj
2632      &          +ees0m(num_conti,i)/sss1*rmij*xj*sssgrad1
2633               gggm(2)=gggm(2)+ees0mijp*yj
2634      &          +ees0m(num_conti,i)/sss1*rmij*yj*sssgrad1
2635               gggm(3)=gggm(3)+ees0mijp*zj
2636      &          +ees0m(num_conti,i)/sss1*rmij*zj*sssgrad1
2637 C Derivatives due to the contact function
2638               gacont_hbr(1,num_conti,i)=fprimcont*xj
2639               gacont_hbr(2,num_conti,i)=fprimcont*yj
2640               gacont_hbr(3,num_conti,i)=fprimcont*zj
2641               do k=1,3
2642 c
2643 c 10/24/08 cgrad and ! comments indicate the parts of the code removed
2644 c          following the change of gradient-summation algorithm.
2645 c
2646 cgrad                  ghalfp=0.5D0*gggp(k)
2647 cgrad                  ghalfm=0.5D0*gggm(k)
2648                  gacontp_hb1(k,num_conti,i)=!ghalfp
2649      &             +(ecosap*(dc_norm(k,j)-cosa*dc_norm(k,i))
2650      &             + ecosbp*(erij(k)-cosb*dc_norm(k,i)))*vbld_inv(i+1)
2651      &               *sss1*fac_shield(i)*fac_shield(j)
2652
2653                  gacontp_hb2(k,num_conti,i)=!ghalfp
2654      &             +(ecosap*(dc_norm(k,i)-cosa*dc_norm(k,j))
2655      &             + ecosgp*(erij(k)-cosg*dc_norm(k,j)))*vbld_inv(j+1)
2656      &               *sss1*fac_shield(i)*fac_shield(j)
2657
2658                  gacontp_hb3(k,num_conti,i)=gggp(k)
2659      &               *sss1*fac_shield(i)*fac_shield(j)
2660
2661                  gacontm_hb1(k,num_conti,i)=!ghalfm
2662      &             +(ecosam*(dc_norm(k,j)-cosa*dc_norm(k,i))
2663      &             + ecosbm*(erij(k)-cosb*dc_norm(k,i)))*vbld_inv(i+1)
2664      &              *sss1*fac_shield(i)*fac_shield(j)
2665
2666                  gacontm_hb2(k,num_conti,i)=!ghalfm
2667      &             +(ecosam*(dc_norm(k,i)-cosa*dc_norm(k,j))
2668      &             + ecosgm*(erij(k)-cosg*dc_norm(k,j)))*vbld_inv(j+1)
2669      &              *sss1*fac_shield(i)*fac_shield(j)
2670
2671                  gacontm_hb3(k,num_conti,i)=gggm(k)
2672      &              *sss1*fac_shield(i)*fac_shield(j)
2673
2674               enddo
2675             ENDIF ! wcorr
2676           endif  ! num_conti.le.maxconts
2677         endif  ! fcont.gt.0
2678       endif    ! j.gt.i+1
2679 #endif
2680       if (wturn3.gt.0.0d0 .or. wturn4.gt.0.0d0) then
2681         do k=1,4
2682           do l=1,3
2683             ghalf=0.5d0*agg(l,k)
2684             aggi(l,k)=aggi(l,k)+ghalf
2685             aggi1(l,k)=aggi1(l,k)+agg(l,k)
2686             aggj(l,k)=aggj(l,k)+ghalf
2687           enddo
2688         enddo
2689         if (j.eq.nres-1 .and. i.lt.j-2) then
2690           do k=1,4
2691             do l=1,3
2692               aggj1(l,k)=aggj1(l,k)+agg(l,k)
2693             enddo
2694           enddo
2695         endif
2696       endif
2697 c          t_eelecij=t_eelecij+MPI_Wtime()-time00
2698       return
2699       end
2700 C-----------------------------------------------------------------------
2701       subroutine evdwpp_short(evdw1)
2702 C
2703 C Compute Evdwpp
2704 C
2705       implicit none
2706       include 'DIMENSIONS'
2707       include 'COMMON.CONTROL'
2708       include 'COMMON.IOUNITS'
2709       include 'COMMON.GEO'
2710       include 'COMMON.VAR'
2711       include 'COMMON.LOCAL'
2712       include 'COMMON.CHAIN'
2713       include 'COMMON.DERIV'
2714       include 'COMMON.INTERACT'
2715 c      include 'COMMON.CONTACTS'
2716       include 'COMMON.TORSION'
2717       include 'COMMON.VECTORS'
2718       include 'COMMON.FFIELD'
2719       include "COMMON.SPLITELE"
2720       double precision ggg(3)
2721       integer xshift,yshift,zshift
2722 c 4/26/02 - AL scaling factor for 1,4 repulsive VDW interactions
2723 #ifdef MOMENT
2724       double precision scal_el /1.0d0/
2725 #else
2726       double precision scal_el /0.5d0/
2727 #endif
2728 c      write (iout,*) "evdwpp_short"
2729       integer i,j,k,iteli,itelj,num_conti,ind,isubchap
2730       double precision dxi,dyi,dzi,dxj,dyj,dzj,aaa,bbb
2731       double precision xj,yj,zj,rij,rrmij,r3ij,r6ij,evdw1,
2732      & dx_normi,dy_normi,dz_normi,xmedi,ymedi,zmedi,
2733      & dx_normj,dy_normj,dz_normj,rmij,ev1,ev2,evdwij,facvdw
2734       double precision xj_safe,yj_safe,zj_safe,xj_temp,yj_temp,zj_temp,
2735      & dist_temp, dist_init,sss_grad
2736       double precision sscale,sscagrad
2737       integer ikont
2738       evdw1=0.0D0
2739 C      print *,"WCHODZE"
2740 c      write (iout,*) "iatel_s_vdw",iatel_s_vdw,
2741 c     & " iatel_e_vdw",iatel_e_vdw
2742 c      call flush(iout)
2743 c      do i=iatel_s_vdw,iatel_e_vdw
2744       do ikont=g_listpp_vdw_start,g_listpp_vdw_end
2745         i=newcontlistpp_vdwi(ikont)
2746         j=newcontlistpp_vdwj(ikont)
2747         if (itype(i).eq.ntyp1.or. itype(i+1).eq.ntyp1) cycle
2748         dxi=dc(1,i)
2749         dyi=dc(2,i)
2750         dzi=dc(3,i)
2751         dx_normi=dc_norm(1,i)
2752         dy_normi=dc_norm(2,i)
2753         dz_normi=dc_norm(3,i)
2754         xmedi=c(1,i)+0.5d0*dxi
2755         ymedi=c(2,i)+0.5d0*dyi
2756         zmedi=c(3,i)+0.5d0*dzi
2757         xmedi=mod(xmedi,boxxsize)
2758         if (xmedi.lt.0.0d0) xmedi=xmedi+boxxsize
2759         ymedi=mod(ymedi,boxysize)
2760         if (ymedi.lt.0.0d0) ymedi=ymedi+boxysize
2761         zmedi=mod(zmedi,boxzsize)
2762         if (zmedi.lt.0.0d0) zmedi=zmedi+boxzsize
2763         num_conti=0
2764 c        write (iout,*) 'i',i,' ielstart',ielstart_vdw(i),
2765 c     &   ' ielend',ielend_vdw(i)
2766 c        call flush(iout)
2767 c        do j=ielstart_vdw(i),ielend_vdw(i)
2768           if (itype(j).eq.ntyp1 .or. itype(j+1).eq.ntyp1) cycle
2769           ind=ind+1
2770           iteli=itel(i)
2771           itelj=itel(j)
2772           if (j.eq.i+2 .and. itelj.eq.2) iteli=2
2773           aaa=app(iteli,itelj)
2774           bbb=bpp(iteli,itelj)
2775           dxj=dc(1,j)
2776           dyj=dc(2,j)
2777           dzj=dc(3,j)
2778           dx_normj=dc_norm(1,j)
2779           dy_normj=dc_norm(2,j)
2780           dz_normj=dc_norm(3,j)
2781           xj=c(1,j)+0.5D0*dxj
2782           yj=c(2,j)+0.5D0*dyj
2783           zj=c(3,j)+0.5D0*dzj
2784           xj=mod(xj,boxxsize)
2785           if (xj.lt.0) xj=xj+boxxsize
2786           yj=mod(yj,boxysize)
2787           if (yj.lt.0) yj=yj+boxysize
2788           zj=mod(zj,boxzsize)
2789           if (zj.lt.0) zj=zj+boxzsize
2790           dist_init=(xj-xmedi)**2+(yj-ymedi)**2+(zj-zmedi)**2
2791           xj_safe=xj
2792           yj_safe=yj
2793           zj_safe=zj
2794           isubchap=0
2795           do xshift=-1,1
2796           do yshift=-1,1
2797           do zshift=-1,1
2798               xj=xj_safe+xshift*boxxsize
2799               yj=yj_safe+yshift*boxysize
2800               zj=zj_safe+zshift*boxzsize
2801               dist_temp=(xj-xmedi)**2+(yj-ymedi)**2+(zj-zmedi)**2
2802               if(dist_temp.lt.dist_init) then
2803                 dist_init=dist_temp
2804                 xj_temp=xj
2805                 yj_temp=yj
2806                 zj_temp=zj
2807                 isubchap=1
2808               endif
2809            enddo
2810            enddo
2811            enddo
2812            if (isubchap.eq.1) then
2813               xj=xj_temp-xmedi
2814               yj=yj_temp-ymedi
2815               zj=zj_temp-zmedi
2816            else
2817               xj=xj_safe-xmedi
2818               yj=yj_safe-ymedi
2819               zj=zj_safe-zmedi
2820            endif
2821           rij=xj*xj+yj*yj+zj*zj
2822           rrmij=1.0D0/rij
2823           rij=dsqrt(rij)
2824 c            sss=sscale(rij/rpp(iteli,itelj))
2825 c            sssgrad=sscagrad(rij/rpp(iteli,itelj))
2826           sss=sscale(rij/rpp(iteli,itelj),r_cut_respa)
2827           sssgrad=sscagrad(rij/rpp(iteli,itelj),r_cut_respa)
2828           if (sss.gt.0.0d0) then
2829             rmij=1.0D0/rij
2830             r3ij=rrmij*rmij
2831             r6ij=r3ij*r3ij
2832             ev1=aaa*r6ij*r6ij
2833 c 4/26/02 - AL scaling down 1,4 repulsive VDW interactions
2834             if (j.eq.i+2) ev1=scal_el*ev1
2835             ev2=bbb*r6ij
2836             evdwij=ev1+ev2
2837             if (energy_dec) then
2838               write (iout,'(a6,2i5,0pf7.3,f7.3)') 'evdw1',i,j,evdwij,sss
2839             endif
2840             evdw1=evdw1+evdwij*sss
2841             if (energy_dec) write (iout,'(a10,2i5,0pf7.3)')
2842      &        'evdw1_sum',i,j,evdw1
2843 C
2844 C Calculate contributions to the Cartesian gradient.
2845 C
2846             facvdw=-6*rrmij*(ev1+evdwij)*sss
2847             ggg(1)=facvdw*xj+sssgrad*rmij*evdwij*xj/rpp(iteli,itelj)
2848             ggg(2)=facvdw*yj+sssgrad*rmij*evdwij*yj/rpp(iteli,itelj)
2849             ggg(3)=facvdw*zj+sssgrad*rmij*evdwij*zj/rpp(iteli,itelj)
2850 C            ggg(1)=facvdw*xj
2851 C            ggg(2)=facvdw*yj
2852 C            ggg(3)=facvdw*zj
2853             do k=1,3
2854               gvdwpp(k,j)=gvdwpp(k,j)+ggg(k)
2855               gvdwpp(k,i)=gvdwpp(k,i)-ggg(k)
2856             enddo
2857           endif
2858 c        enddo ! j
2859       enddo   ! i
2860       return
2861       end
2862 C-----------------------------------------------------------------------------
2863       subroutine escp_long(evdw2,evdw2_14)
2864 C
2865 C This subroutine calculates the excluded-volume interaction energy between
2866 C peptide-group centers and side chains and its gradient in virtual-bond and
2867 C side-chain vectors.
2868 C
2869       implicit real*8 (a-h,o-z)
2870       include 'DIMENSIONS'
2871       include 'COMMON.GEO'
2872       include 'COMMON.VAR'
2873       include 'COMMON.LOCAL'
2874       include 'COMMON.CHAIN'
2875       include 'COMMON.DERIV'
2876       include 'COMMON.INTERACT'
2877       include 'COMMON.FFIELD'
2878       include 'COMMON.IOUNITS'
2879       include 'COMMON.CONTROL'
2880       include "COMMON.SPLITELE"
2881       logical lprint_short
2882       common /shortcheck/ lprint_short
2883       double precision ggg(3)
2884       integer xshift,yshift,zshift
2885       integer i,iint,j,k,iteli,itypj,subchap
2886       double precision xi,yi,zi,xj,yj,zj,rrij,sss1,sssgrad1,
2887      & fac,e1,e2,rij
2888       double precision evdw2,evdw2_14,evdwij
2889       double precision xj_safe,yj_safe,zj_safe,xj_temp,yj_temp,zj_temp,
2890      & dist_temp, dist_init
2891       double precision sscale,sscagrad
2892       integer ikont
2893       if (energy_dec) write (iout,*) "escp_long:",r_cut,rlamb
2894       evdw2=0.0D0
2895       evdw2_14=0.0d0
2896 CD        print '(a)','Enter ESCP KURWA'
2897 cd    write (iout,*) 'iatscp_s=',iatscp_s,' iatscp_e=',iatscp_e
2898 c      if (lprint_short)
2899 c     &  write (iout,*) 'ESCP_LONG iatscp_s=',iatscp_s,
2900 c     & ' iatscp_e=',iatscp_e
2901 c      do i=iatscp_s,iatscp_e
2902       do ikont=g_listscp_start,g_listscp_end
2903         i=newcontlistscpi(ikont)
2904         j=newcontlistscpj(ikont)
2905         if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1) cycle
2906         iteli=itel(i)
2907         xi=0.5D0*(c(1,i)+c(1,i+1))
2908         yi=0.5D0*(c(2,i)+c(2,i+1))
2909         zi=0.5D0*(c(3,i)+c(3,i+1))
2910         xi=mod(xi,boxxsize)
2911         if (xi.lt.0) xi=xi+boxxsize
2912         yi=mod(yi,boxysize)
2913         if (yi.lt.0) yi=yi+boxysize
2914         zi=mod(zi,boxzsize)
2915         if (zi.lt.0) zi=zi+boxzsize
2916
2917 c        do iint=1,nscp_gr(i)
2918
2919 c        do j=iscpstart(i,iint),iscpend(i,iint)
2920           itypj=iabs(itype(j))
2921           if (itypj.eq.ntyp1) cycle
2922 C Uncomment following three lines for SC-p interactions
2923 c         xj=c(1,nres+j)-xi
2924 c         yj=c(2,nres+j)-yi
2925 c         zj=c(3,nres+j)-zi
2926 C Uncomment following three lines for Ca-p interactions
2927           xj=c(1,j)
2928           yj=c(2,j)
2929           zj=c(3,j)
2930 c corrected by AL
2931           xj=mod(xj,boxxsize)
2932           if (xj.lt.0) xj=xj+boxxsize
2933           yj=mod(yj,boxysize)
2934           if (yj.lt.0) yj=yj+boxysize
2935           zj=mod(zj,boxzsize)
2936           if (zj.lt.0) zj=zj+boxzsize
2937 c end correction
2938           dist_init=(xj-xi)**2+(yj-yi)**2+(zj-zi)**2
2939           xj_safe=xj
2940           yj_safe=yj
2941           zj_safe=zj
2942           subchap=0
2943           do xshift=-1,1
2944           do yshift=-1,1
2945           do zshift=-1,1
2946           xj=xj_safe+xshift*boxxsize
2947           yj=yj_safe+yshift*boxysize
2948           zj=zj_safe+zshift*boxzsize
2949           dist_temp=(xj-xi)**2+(yj-yi)**2+(zj-zi)**2
2950           if(dist_temp.lt.dist_init) then
2951             dist_init=dist_temp
2952             xj_temp=xj
2953             yj_temp=yj
2954             zj_temp=zj
2955             subchap=1
2956           endif
2957           enddo
2958           enddo
2959           enddo
2960           if (subchap.eq.1) then
2961             xj=xj_temp-xi
2962             yj=yj_temp-yi
2963             zj=zj_temp-zi
2964           else
2965             xj=xj_safe-xi
2966             yj=yj_safe-yi
2967             zj=zj_safe-zi
2968           endif
2969
2970           rrij=1.0D0/(xj*xj+yj*yj+zj*zj)
2971
2972           sss1=sscale(1.0d0/(dsqrt(rrij)),r_cut_int)
2973           if (sss1.eq.0) cycle
2974           sss=sscale(1.0d0/(dsqrt(rrij)*rscp(itypj,iteli)),r_cut_respa)
2975           sssgrad=
2976      &      sscagrad(1.0d0/(dsqrt(rrij)*rscp(itypj,iteli)),r_cut_respa)
2977           sssgrad1=sscagrad(1.0d0/dsqrt(rrij),r_cut_int)
2978           if (energy_dec) write (iout,*) "rrij",1.0d0/dsqrt(rrij),
2979      &     " rscp",rscp(itypj,iteli)," subchap",subchap," sss",sss
2980           if (sss.lt.1.0d0) then
2981             fac=rrij**expon2
2982             e1=fac*fac*aad(itypj,iteli)
2983             e2=fac*bad(itypj,iteli)
2984             if (iabs(j-i) .le. 2) then
2985               e1=scal14*e1
2986               e2=scal14*e2
2987               evdw2_14=evdw2_14+(e1+e2)*(1.0d0-sss)*sss1
2988             endif
2989             evdwij=e1+e2
2990             evdw2=evdw2+evdwij*(1.0d0-sss)*sss1
2991             if (energy_dec) write (iout,'(a6,2i5,2(0pf7.3))')
2992      &          'evdw2',i,j,sss,evdwij
2993 C
2994 C Calculate contributions to the gradient in the virtual-bond and SC vectors.
2995 C
2996
2997             fac=-(evdwij+e1)*rrij*(1.0d0-sss)*sss1
2998             fac=fac+evdwij*dsqrt(rrij)*(-sssgrad/rscp(itypj,iteli)
2999      &        +sssgrad1)/expon
3000             ggg(1)=xj*fac
3001             ggg(2)=yj*fac
3002             ggg(3)=zj*fac
3003 C Uncomment following three lines for SC-p interactions
3004 c           do k=1,3
3005 c             gradx_scp(k,j)=gradx_scp(k,j)+ggg(k)
3006 c           enddo
3007 C Uncomment following line for SC-p interactions
3008 c             gradx_scp(k,j)=gradx_scp(k,j)+ggg(k)
3009             do k=1,3
3010               gvdwc_scpp(k,i)=gvdwc_scpp(k,i)-ggg(k)
3011               gvdwc_scp(k,j)=gvdwc_scp(k,j)+ggg(k)
3012             enddo
3013           endif
3014 c        enddo
3015
3016 c        enddo ! iint
3017       enddo ! i
3018       do i=1,nct
3019         do j=1,3
3020           gvdwc_scp(j,i)=expon*gvdwc_scp(j,i)
3021           gvdwc_scpp(j,i)=expon*gvdwc_scpp(j,i)
3022           gradx_scp(j,i)=expon*gradx_scp(j,i)
3023         enddo
3024       enddo
3025 C******************************************************************************
3026 C
3027 C                              N O T E !!!
3028 C
3029 C To save time the factor EXPON has been extracted from ALL components
3030 C of GVDWC and GRADX. Remember to multiply them by this factor before further
3031 C use!
3032 C
3033 C******************************************************************************
3034       return
3035       end
3036 C-----------------------------------------------------------------------------
3037       subroutine escp_short(evdw2,evdw2_14)
3038 C
3039 C This subroutine calculates the excluded-volume interaction energy between
3040 C peptide-group centers and side chains and its gradient in virtual-bond and
3041 C side-chain vectors.
3042 C
3043       implicit none
3044       include 'DIMENSIONS'
3045       include 'COMMON.GEO'
3046       include 'COMMON.VAR'
3047       include 'COMMON.LOCAL'
3048       include 'COMMON.CHAIN'
3049       include 'COMMON.DERIV'
3050       include 'COMMON.INTERACT'
3051       include 'COMMON.FFIELD'
3052       include 'COMMON.IOUNITS'
3053       include 'COMMON.CONTROL'
3054       include "COMMON.SPLITELE"
3055       integer xshift,yshift,zshift
3056       logical lprint_short
3057       common /shortcheck/ lprint_short
3058       integer i,iint,j,k,iteli,itypj,subchap
3059       double precision xi,yi,zi,xj,yj,zj,rrij,sss1,sssgrad1,
3060      & fac,e1,e2,rij
3061       double precision evdw2,evdw2_14,evdwij
3062       double precision xj_safe,yj_safe,zj_safe,xj_temp,yj_temp,zj_temp,
3063      & dist_temp, dist_init
3064       double precision ggg(3)
3065       double precision sscale,sscagrad
3066       evdw2=0.0D0
3067       evdw2_14=0.0d0
3068 cd    print '(a)','Enter ESCP'
3069 c      if (lprint_short)
3070 c     &  write (iout,*) 'ESCP_SHORT iatscp_s=',iatscp_s,
3071 c     & ' iatscp_e=',iatscp_e
3072       if (energy_dec) write (iout,*) "escp_short:",r_cut_int,rlamb
3073       do i=iatscp_s,iatscp_e
3074         if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1) cycle
3075         iteli=itel(i)
3076         xi=0.5D0*(c(1,i)+c(1,i+1))
3077         yi=0.5D0*(c(2,i)+c(2,i+1))
3078         zi=0.5D0*(c(3,i)+c(3,i+1))
3079         xi=mod(xi,boxxsize)
3080         if (xi.lt.0) xi=xi+boxxsize
3081         yi=mod(yi,boxysize)
3082         if (yi.lt.0) yi=yi+boxysize
3083         zi=mod(zi,boxzsize)
3084         if (zi.lt.0) zi=zi+boxzsize
3085
3086 c        if (lprint_short)
3087 c     &    write (iout,*) "i",i," itype",itype(i),itype(i+1),
3088 c     &     " nscp_gr",nscp_gr(i)
3089         do iint=1,nscp_gr(i)
3090
3091         do j=iscpstart(i,iint),iscpend(i,iint)
3092           itypj=iabs(itype(j))
3093 c        if (lprint_short)
3094 c     &    write (iout,*) "j",j," itypj",itypj
3095           if (itypj.eq.ntyp1) cycle
3096 C Uncomment following three lines for SC-p interactions
3097 c         xj=c(1,nres+j)-xi
3098 c         yj=c(2,nres+j)-yi
3099 c         zj=c(3,nres+j)-zi
3100 C Uncomment following three lines for Ca-p interactions
3101           xj=c(1,j)
3102           yj=c(2,j)
3103           zj=c(3,j)
3104 c corrected by AL
3105           xj=mod(xj,boxxsize)
3106           if (xj.lt.0) xj=xj+boxxsize
3107           yj=mod(yj,boxysize)
3108           if (yj.lt.0) yj=yj+boxysize
3109           zj=mod(zj,boxzsize)
3110           if (zj.lt.0) zj=zj+boxzsize
3111 c end correction
3112           dist_init=(xj-xi)**2+(yj-yi)**2+(zj-zi)**2
3113 c          if (lprint_short) then
3114 c            write (iout,*) i,j,xi,yi,zi,xj,yj,zj
3115 c            write (iout,*) "dist_init",dsqrt(dist_init)
3116 c          endif
3117           xj_safe=xj
3118           yj_safe=yj
3119           zj_safe=zj
3120           subchap=0
3121           do xshift=-1,1
3122           do yshift=-1,1
3123           do zshift=-1,1
3124           xj=xj_safe+xshift*boxxsize
3125           yj=yj_safe+yshift*boxysize
3126           zj=zj_safe+zshift*boxzsize
3127           dist_temp=(xj-xi)**2+(yj-yi)**2+(zj-zi)**2
3128           if(dist_temp.lt.dist_init) then
3129             dist_init=dist_temp
3130             xj_temp=xj
3131             yj_temp=yj
3132             zj_temp=zj
3133             subchap=1
3134           endif
3135           enddo
3136           enddo
3137           enddo
3138 c          if (lprint_short) write (iout,*) "dist_temp",dsqrt(dist_temp)
3139           if (subchap.eq.1) then
3140              xj=xj_temp-xi
3141              yj=yj_temp-yi
3142              zj=zj_temp-zi
3143           else
3144              xj=xj_safe-xi
3145              yj=yj_safe-yi
3146              zj=zj_safe-zi
3147           endif
3148           rrij=1.0D0/(xj*xj+yj*yj+zj*zj)
3149 c          sss=sscale(1.0d0/(dsqrt(rrij)*rscp(itypj,iteli)))
3150 c          sssgrad=sscagrad(1.0d0/(dsqrt(rrij)*rscp(itypj,iteli)))
3151           sss=sscale(1.0d0/(dsqrt(rrij)*rscp(itypj,iteli)),r_cut_respa)
3152           sssgrad=sscagrad(1.0d0/(dsqrt(rrij)*rscp(itypj,iteli)),
3153      &        r_cut_respa)
3154           if (energy_dec) write (iout,*) "rrij",1.0d0/dsqrt(rrij),
3155      &     " rscp",rscp(itypj,iteli)," subchap",subchap," sss",sss
3156 c          if (lprint_short) write (iout,*) "rij",1.0/dsqrt(rrij),
3157 c     &     " subchap",subchap," sss",sss
3158           if (sss.gt.0.0d0) then
3159
3160             fac=rrij**expon2
3161             e1=fac*fac*aad(itypj,iteli)
3162             e2=fac*bad(itypj,iteli)
3163             if (iabs(j-i) .le. 2) then
3164               e1=scal14*e1
3165               e2=scal14*e2
3166               evdw2_14=evdw2_14+(e1+e2)*sss
3167             endif
3168             evdwij=e1+e2
3169             evdw2=evdw2+evdwij*sss
3170             if (energy_dec) write (iout,'(a6,2i5,2(0pf7.3))')
3171      &          'evdw2',i,j,sss,evdwij
3172 C
3173 C Calculate contributions to the gradient in the virtual-bond and SC vectors.
3174 C
3175             fac=-(evdwij+e1)*rrij*sss
3176             fac=fac+evdwij*sssgrad*dsqrt(rrij)/rscp(itypj,iteli)/expon
3177             ggg(1)=xj*fac
3178             ggg(2)=yj*fac
3179             ggg(3)=zj*fac
3180 C Uncomment following three lines for SC-p interactions
3181 c           do k=1,3
3182 c             gradx_scp(k,j)=gradx_scp(k,j)+ggg(k)
3183 c           enddo
3184 C Uncomment following line for SC-p interactions
3185 c             gradx_scp(k,j)=gradx_scp(k,j)+ggg(k)
3186             do k=1,3
3187               gvdwc_scpp(k,i)=gvdwc_scpp(k,i)-ggg(k)
3188               gvdwc_scp(k,j)=gvdwc_scp(k,j)+ggg(k)
3189             enddo
3190           endif
3191         enddo
3192
3193         enddo ! iint
3194       enddo ! i
3195       do i=1,nct
3196         do j=1,3
3197           gvdwc_scp(j,i)=expon*gvdwc_scp(j,i)
3198           gvdwc_scpp(j,i)=expon*gvdwc_scpp(j,i)
3199           gradx_scp(j,i)=expon*gradx_scp(j,i)
3200         enddo
3201       enddo
3202 C******************************************************************************
3203 C
3204 C                              N O T E !!!
3205 C
3206 C To save time the factor EXPON has been extracted from ALL components
3207 C of GVDWC and GRADX. Remember to multiply them by this factor before further
3208 C use!
3209 C
3210 C******************************************************************************
3211       return
3212       end