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