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