source/unres/src_MD-M/energy_p_new-sep_barrier.F

   1 C-----------------------------------------------------------------------
   2       double precision function sscale(r)
   3       double precision r,gamm
   4       include "COMMON.SPLITELE"
   5       if(r.lt.r_cut-rlamb) then
   6         sscale=1.0d0
   7       else if(r.le.r_cut.and.r.ge.r_cut-rlamb) then
   8         gamm=(r-(r_cut-rlamb))/rlamb
   9         sscale=1.0d0+gamm*gamm*(2*gamm-3.0d0)
  10       else
  11         sscale=0d0
  12       endif
  13       return
  14       end
  15 C-----------------------------------------------------------------------
  16       subroutine elj_long(evdw)
  17 C
  18 C This subroutine calculates the interaction energy of nonbonded side chains
  19 C assuming the LJ potential of interaction.
  20 C
  21       implicit real*8 (a-h,o-z)
  22       include 'DIMENSIONS'
  23       parameter (accur=1.0d-10)
  24       include 'COMMON.GEO'
  25       include 'COMMON.VAR'
  26       include 'COMMON.LOCAL'
  27       include 'COMMON.CHAIN'
  28       include 'COMMON.DERIV'
  29       include 'COMMON.INTERACT'
  30       include 'COMMON.TORSION'
  31       include 'COMMON.SBRIDGE'
  32       include 'COMMON.NAMES'
  33       include 'COMMON.IOUNITS'
  34       include 'COMMON.CONTACTS'
  35       dimension gg(3)
  36 c      write(iout,*)'Entering ELJ nnt=',nnt,' nct=',nct,' expon=',expon
  37       evdw=0.0D0
  38       do i=iatsc_s,iatsc_e
  39         itypi=itype(i)
  40         if (itypi.eq.21) cycle
  41         itypi1=itype(i+1)
  42         xi=c(1,nres+i)
  43         yi=c(2,nres+i)
  44         zi=c(3,nres+i)
  45 C
  46 C Calculate SC interaction energy.
  47 C
  48         do iint=1,nint_gr(i)
  49 cd        write (iout,*) 'i=',i,' iint=',iint,' istart=',istart(i,iint),
  50 cd   &                  'iend=',iend(i,iint)
  51           do j=istart(i,iint),iend(i,iint)
  52             itypj=itype(j)
  53             if (itypj.eq.21) cycle
  54             xj=c(1,nres+j)-xi
  55             yj=c(2,nres+j)-yi
  56             zj=c(3,nres+j)-zi
  57             rij=xj*xj+yj*yj+zj*zj
  58             sss=sscale(dsqrt(rij)/sigma(itypi,itypj))
  59             if (sss.lt.1.0d0) then
  60               rrij=1.0D0/rij
  61               eps0ij=eps(itypi,itypj)
  62               fac=rrij**expon2
  63               e1=fac*fac*aa(itypi,itypj)
  64               e2=fac*bb(itypi,itypj)
  65               evdwij=e1+e2
  66               evdw=evdw+(1.0d0-sss)*evdwij
  67 C
  68 C Calculate the components of the gradient in DC and X
  69 C
  70               fac=-rrij*(e1+evdwij)*(1.0d0-sss)
  71               gg(1)=xj*fac
  72               gg(2)=yj*fac
  73               gg(3)=zj*fac
  74               do k=1,3
  75                 gvdwx(k,i)=gvdwx(k,i)-gg(k)
  76                 gvdwx(k,j)=gvdwx(k,j)+gg(k)
  77                 gvdwc(k,i)=gvdwc(k,i)-gg(k)
  78                 gvdwc(k,j)=gvdwc(k,j)+gg(k)
  79               enddo
  80             endif
  81           enddo      ! j
  82         enddo        ! iint
  83       enddo          ! i
  84       do i=1,nct
  85         do j=1,3
  86           gvdwc(j,i)=expon*gvdwc(j,i)
  87           gvdwx(j,i)=expon*gvdwx(j,i)
  88         enddo
  89       enddo
  90 C******************************************************************************
  91 C
  92 C                              N O T E !!!
  93 C
  94 C To save time, the factor of EXPON has been extracted from ALL components
  95 C of GVDWC and GRADX. Remember to multiply them by this factor before further
  96 C use!
  97 C
  98 C******************************************************************************
  99       return
 100       end
 101 C-----------------------------------------------------------------------
 102       subroutine elj_short(evdw)
 103 C
 104 C This subroutine calculates the interaction energy of nonbonded side chains
 105 C assuming the LJ potential of interaction.
 106 C
 107       implicit real*8 (a-h,o-z)
 108       include 'DIMENSIONS'
 109       parameter (accur=1.0d-10)
 110       include 'COMMON.GEO'
 111       include 'COMMON.VAR'
 112       include 'COMMON.LOCAL'
 113       include 'COMMON.CHAIN'
 114       include 'COMMON.DERIV'
 115       include 'COMMON.INTERACT'
 116       include 'COMMON.TORSION'
 117       include 'COMMON.SBRIDGE'
 118       include 'COMMON.NAMES'
 119       include 'COMMON.IOUNITS'
 120       include 'COMMON.CONTACTS'
 121       dimension gg(3)
 122 c      write(iout,*)'Entering ELJ nnt=',nnt,' nct=',nct,' expon=',expon
 123       evdw=0.0D0
 124       do i=iatsc_s,iatsc_e
 125         itypi=itype(i)
 126         if (itypi.eq.21) cycle
 127         itypi1=itype(i+1)
 128         xi=c(1,nres+i)
 129         yi=c(2,nres+i)
 130         zi=c(3,nres+i)
 131 C Change 12/1/95
 132         num_conti=0
 133 C
 134 C Calculate SC interaction energy.
 135 C
 136         do iint=1,nint_gr(i)
 137 cd        write (iout,*) 'i=',i,' iint=',iint,' istart=',istart(i,iint),
 138 cd   &                  'iend=',iend(i,iint)
 139           do j=istart(i,iint),iend(i,iint)
 140             itypj=itype(j)
 141             if (itypj.eq.21) cycle
 142             xj=c(1,nres+j)-xi
 143             yj=c(2,nres+j)-yi
 144             zj=c(3,nres+j)-zi
 145 C Change 12/1/95 to calculate four-body interactions
 146             rij=xj*xj+yj*yj+zj*zj
 147             sss=sscale(dsqrt(rij)/sigma(itypi,itypj))
 148             if (sss.gt.0.0d0) then
 149               rrij=1.0D0/rij
 150               eps0ij=eps(itypi,itypj)
 151               fac=rrij**expon2
 152               e1=fac*fac*aa(itypi,itypj)
 153               e2=fac*bb(itypi,itypj)
 154               evdwij=e1+e2
 155               evdw=evdw+sss*evdwij
 156 C
 157 C Calculate the components of the gradient in DC and X
 158 C
 159               fac=-rrij*(e1+evdwij)*sss
 160               gg(1)=xj*fac
 161               gg(2)=yj*fac
 162               gg(3)=zj*fac
 163               do k=1,3
 164                 gvdwx(k,i)=gvdwx(k,i)-gg(k)
 165                 gvdwx(k,j)=gvdwx(k,j)+gg(k)
 166                 gvdwc(k,i)=gvdwc(k,i)-gg(k)
 167                 gvdwc(k,j)=gvdwc(k,j)+gg(k)
 168               enddo
 169             endif
 170           enddo      ! j
 171         enddo        ! iint
 172       enddo          ! i
 173       do i=1,nct
 174         do j=1,3
 175           gvdwc(j,i)=expon*gvdwc(j,i)
 176           gvdwx(j,i)=expon*gvdwx(j,i)
 177         enddo
 178       enddo
 179 C******************************************************************************
 180 C
 181 C                              N O T E !!!
 182 C
 183 C To save time, the factor of EXPON has been extracted from ALL components
 184 C of GVDWC and GRADX. Remember to multiply them by this factor before further
 185 C use!
 186 C
 187 C******************************************************************************
 188       return
 189       end
 190 C-----------------------------------------------------------------------------
 191       subroutine eljk_long(evdw)
 192 C
 193 C This subroutine calculates the interaction energy of nonbonded side chains
 194 C assuming the LJK potential of interaction.
 195 C
 196       implicit real*8 (a-h,o-z)
 197       include 'DIMENSIONS'
 198       include 'COMMON.GEO'
 199       include 'COMMON.VAR'
 200       include 'COMMON.LOCAL'
 201       include 'COMMON.CHAIN'
 202       include 'COMMON.DERIV'
 203       include 'COMMON.INTERACT'
 204       include 'COMMON.IOUNITS'
 205       include 'COMMON.NAMES'
 206       dimension gg(3)
 207       logical scheck
 208 c     print *,'Entering ELJK nnt=',nnt,' nct=',nct,' expon=',expon
 209       evdw=0.0D0
 210       do i=iatsc_s,iatsc_e
 211         itypi=itype(i)
 212         if (itypi.eq.21) cycle
 213         itypi1=itype(i+1)
 214         xi=c(1,nres+i)
 215         yi=c(2,nres+i)
 216         zi=c(3,nres+i)
 217 C
 218 C Calculate SC interaction energy.
 219 C
 220         do iint=1,nint_gr(i)
 221           do j=istart(i,iint),iend(i,iint)
 222             itypj=itype(j)
 223             if (itypj.eq.21) cycle
 224             xj=c(1,nres+j)-xi
 225             yj=c(2,nres+j)-yi
 226             zj=c(3,nres+j)-zi
 227             rrij=1.0D0/(xj*xj+yj*yj+zj*zj)
 228             fac_augm=rrij**expon
 229             e_augm=augm(itypi,itypj)*fac_augm
 230             r_inv_ij=dsqrt(rrij)
 231             rij=1.0D0/r_inv_ij
 232             sss=sscale(rij/sigma(itypi,itypj))
 233             if (sss.lt.1.0d0) then
 234               r_shift_inv=1.0D0/(rij+r0(itypi,itypj)-sigma(itypi,itypj))
 235               fac=r_shift_inv**expon
 236               e1=fac*fac*aa(itypi,itypj)
 237               e2=fac*bb(itypi,itypj)
 238               evdwij=e_augm+e1+e2
 239 cd            sigm=dabs(aa(itypi,itypj)/bb(itypi,itypj))**(1.0D0/6.0D0)
 240 cd            epsi=bb(itypi,itypj)**2/aa(itypi,itypj)
 241 cd            write (iout,'(2(a3,i3,2x),8(1pd12.4)/2(3(1pd12.4),5x)/)')
 242 cd   &          restyp(itypi),i,restyp(itypj),j,aa(itypi,itypj),
 243 cd   &          bb(itypi,itypj),augm(itypi,itypj),epsi,sigm,
 244 cd   &          sigma(itypi,itypj),1.0D0/dsqrt(rrij),evdwij,
 245 cd   &          (c(k,i),k=1,3),(c(k,j),k=1,3)
 246               evdw=evdw+(1.0d0-sss)*evdwij
 247 C
 248 C Calculate the components of the gradient in DC and X
 249 C
 250               fac=-2.0D0*rrij*e_augm-r_inv_ij*r_shift_inv*(e1+e1+e2)
 251               fac=fac*(1.0d0-sss)
 252               gg(1)=xj*fac
 253               gg(2)=yj*fac
 254               gg(3)=zj*fac
 255               do k=1,3
 256                 gvdwx(k,i)=gvdwx(k,i)-gg(k)
 257                 gvdwx(k,j)=gvdwx(k,j)+gg(k)
 258                 gvdwc(k,i)=gvdwc(k,i)-gg(k)
 259                 gvdwc(k,j)=gvdwc(k,j)+gg(k)
 260               enddo
 261             endif
 262           enddo      ! j
 263         enddo        ! iint
 264       enddo          ! i
 265       do i=1,nct
 266         do j=1,3
 267           gvdwc(j,i)=expon*gvdwc(j,i)
 268           gvdwx(j,i)=expon*gvdwx(j,i)
 269         enddo
 270       enddo
 271       return
 272       end
 273 C-----------------------------------------------------------------------------
 274       subroutine eljk_short(evdw)
 275 C
 276 C This subroutine calculates the interaction energy of nonbonded side chains
 277 C assuming the LJK potential of interaction.
 278 C
 279       implicit real*8 (a-h,o-z)
 280       include 'DIMENSIONS'
 281       include 'COMMON.GEO'
 282       include 'COMMON.VAR'
 283       include 'COMMON.LOCAL'
 284       include 'COMMON.CHAIN'
 285       include 'COMMON.DERIV'
 286       include 'COMMON.INTERACT'
 287       include 'COMMON.IOUNITS'
 288       include 'COMMON.NAMES'
 289       dimension gg(3)
 290       logical scheck
 291 c     print *,'Entering ELJK nnt=',nnt,' nct=',nct,' expon=',expon
 292       evdw=0.0D0
 293       do i=iatsc_s,iatsc_e
 294         itypi=itype(i)
 295         if (itypi.eq.21) cycle
 296         itypi1=itype(i+1)
 297         xi=c(1,nres+i)
 298         yi=c(2,nres+i)
 299         zi=c(3,nres+i)
 300 C
 301 C Calculate SC interaction energy.
 302 C
 303         do iint=1,nint_gr(i)
 304           do j=istart(i,iint),iend(i,iint)
 305             itypj=itype(j)
 306             if (itypj.eq.21) cycle
 307             xj=c(1,nres+j)-xi
 308             yj=c(2,nres+j)-yi
 309             zj=c(3,nres+j)-zi
 310             rrij=1.0D0/(xj*xj+yj*yj+zj*zj)
 311             fac_augm=rrij**expon
 312             e_augm=augm(itypi,itypj)*fac_augm
 313             r_inv_ij=dsqrt(rrij)
 314             rij=1.0D0/r_inv_ij
 315             sss=sscale(rij/sigma(itypi,itypj))
 316             if (sss.gt.0.0d0) then
 317               r_shift_inv=1.0D0/(rij+r0(itypi,itypj)-sigma(itypi,itypj))
 318               fac=r_shift_inv**expon
 319               e1=fac*fac*aa(itypi,itypj)
 320               e2=fac*bb(itypi,itypj)
 321               evdwij=e_augm+e1+e2
 322 cd            sigm=dabs(aa(itypi,itypj)/bb(itypi,itypj))**(1.0D0/6.0D0)
 323 cd            epsi=bb(itypi,itypj)**2/aa(itypi,itypj)
 324 cd            write (iout,'(2(a3,i3,2x),8(1pd12.4)/2(3(1pd12.4),5x)/)')
 325 cd   &          restyp(itypi),i,restyp(itypj),j,aa(itypi,itypj),
 326 cd   &          bb(itypi,itypj),augm(itypi,itypj),epsi,sigm,
 327 cd   &          sigma(itypi,itypj),1.0D0/dsqrt(rrij),evdwij,
 328 cd   &          (c(k,i),k=1,3),(c(k,j),k=1,3)
 329               evdw=evdw+sss*evdwij
 330 C
 331 C Calculate the components of the gradient in DC and X
 332 C
 333               fac=-2.0D0*rrij*e_augm-r_inv_ij*r_shift_inv*(e1+e1+e2)
 334               fac=fac*sss
 335               gg(1)=xj*fac
 336               gg(2)=yj*fac
 337               gg(3)=zj*fac
 338               do k=1,3
 339                 gvdwx(k,i)=gvdwx(k,i)-gg(k)
 340                 gvdwx(k,j)=gvdwx(k,j)+gg(k)
 341                 gvdwc(k,i)=gvdwc(k,i)-gg(k)
 342                 gvdwc(k,j)=gvdwc(k,j)+gg(k)
 343               enddo
 344             endif
 345           enddo      ! j
 346         enddo        ! iint
 347       enddo          ! i
 348       do i=1,nct
 349         do j=1,3
 350           gvdwc(j,i)=expon*gvdwc(j,i)
 351           gvdwx(j,i)=expon*gvdwx(j,i)
 352         enddo
 353       enddo
 354       return
 355       end
 356 C-----------------------------------------------------------------------------
 357       subroutine ebp_long(evdw)
 358 C
 359 C This subroutine calculates the interaction energy of nonbonded side chains
 360 C assuming the Berne-Pechukas potential of interaction.
 361 C
 362       implicit real*8 (a-h,o-z)
 363       include 'DIMENSIONS'
 364       include 'COMMON.GEO'
 365       include 'COMMON.VAR'
 366       include 'COMMON.LOCAL'
 367       include 'COMMON.CHAIN'
 368       include 'COMMON.DERIV'
 369       include 'COMMON.NAMES'
 370       include 'COMMON.INTERACT'
 371       include 'COMMON.IOUNITS'
 372       include 'COMMON.CALC'
 373       common /srutu/ icall
 374 c     double precision rrsave(maxdim)
 375       logical lprn
 376       evdw=0.0D0
 377 c     print *,'Entering EBP nnt=',nnt,' nct=',nct,' expon=',expon
 378       evdw=0.0D0
 379 c     if (icall.eq.0) then
 380 c       lprn=.true.
 381 c     else
 382         lprn=.false.
 383 c     endif
 384       ind=0
 385       do i=iatsc_s,iatsc_e
 386         itypi=itype(i)
 387         if (itypi.eq.21) cycle
 388         itypi1=itype(i+1)
 389         xi=c(1,nres+i)
 390         yi=c(2,nres+i)
 391         zi=c(3,nres+i)
 392         dxi=dc_norm(1,nres+i)
 393         dyi=dc_norm(2,nres+i)
 394         dzi=dc_norm(3,nres+i)
 395 c        dsci_inv=dsc_inv(itypi)
 396         dsci_inv=vbld_inv(i+nres)
 397 C
 398 C Calculate SC interaction energy.
 399 C
 400         do iint=1,nint_gr(i)
 401           do j=istart(i,iint),iend(i,iint)
 402             ind=ind+1
 403             itypj=itype(j)
 404             if (itypj.eq.21) cycle
 405 c            dscj_inv=dsc_inv(itypj)
 406             dscj_inv=vbld_inv(j+nres)
 407             chi1=chi(itypi,itypj)
 408             chi2=chi(itypj,itypi)
 409             chi12=chi1*chi2
 410             chip1=chip(itypi)
 411             chip2=chip(itypj)
 412             chip12=chip1*chip2
 413             alf1=alp(itypi)
 414             alf2=alp(itypj)
 415             alf12=0.5D0*(alf1+alf2)
 416             xj=c(1,nres+j)-xi
 417             yj=c(2,nres+j)-yi
 418             zj=c(3,nres+j)-zi
 419             dxj=dc_norm(1,nres+j)
 420             dyj=dc_norm(2,nres+j)
 421             dzj=dc_norm(3,nres+j)
 422             rrij=1.0D0/(xj*xj+yj*yj+zj*zj)
 423             rij=dsqrt(rrij)
 424             sss=sscale(1.0d0/(rij*sigmaii(itypi,itypj)))
 425
 426             if (sss.lt.1.0d0) then
 427
 428 C Calculate the angle-dependent terms of energy & contributions to derivatives.
 429               call sc_angular
 430 C Calculate whole angle-dependent part of epsilon and contributions
 431 C to its derivatives
 432               fac=(rrij*sigsq)**expon2
 433               e1=fac*fac*aa(itypi,itypj)
 434               e2=fac*bb(itypi,itypj)
 435               evdwij=eps1*eps2rt*eps3rt*(e1+e2)
 436               eps2der=evdwij*eps3rt
 437               eps3der=evdwij*eps2rt
 438               evdwij=evdwij*eps2rt*eps3rt
 439               evdw=evdw+evdwij*(1.0d0-sss)
 440               if (lprn) then
 441               sigm=dabs(aa(itypi,itypj)/bb(itypi,itypj))**(1.0D0/6.0D0)
 442               epsi=bb(itypi,itypj)**2/aa(itypi,itypj)
 443 cd              write (iout,'(2(a3,i3,2x),15(0pf7.3))')
 444 cd     &          restyp(itypi),i,restyp(itypj),j,
 445 cd     &          epsi,sigm,chi1,chi2,chip1,chip2,
 446 cd     &          eps1,eps2rt**2,eps3rt**2,1.0D0/dsqrt(sigsq),
 447 cd     &          om1,om2,om12,1.0D0/dsqrt(rrij),
 448 cd     &          evdwij
 449               endif
 450 C Calculate gradient components.
 451               e1=e1*eps1*eps2rt**2*eps3rt**2
 452               fac=-expon*(e1+evdwij)
 453               sigder=fac/sigsq
 454               fac=rrij*fac
 455 C Calculate radial part of the gradient
 456               gg(1)=xj*fac
 457               gg(2)=yj*fac
 458               gg(3)=zj*fac
 459 C Calculate the angular part of the gradient and sum add the contributions
 460 C to the appropriate components of the Cartesian gradient.
 461               call sc_grad_scale(1.0d0-sss)
 462             endif
 463           enddo      ! j
 464         enddo        ! iint
 465       enddo          ! i
 466 c     stop
 467       return
 468       end
 469 C-----------------------------------------------------------------------------
 470       subroutine ebp_short(evdw)
 471 C
 472 C This subroutine calculates the interaction energy of nonbonded side chains
 473 C assuming the Berne-Pechukas potential of interaction.
 474 C
 475       implicit real*8 (a-h,o-z)
 476       include 'DIMENSIONS'
 477       include 'COMMON.GEO'
 478       include 'COMMON.VAR'
 479       include 'COMMON.LOCAL'
 480       include 'COMMON.CHAIN'
 481       include 'COMMON.DERIV'
 482       include 'COMMON.NAMES'
 483       include 'COMMON.INTERACT'
 484       include 'COMMON.IOUNITS'
 485       include 'COMMON.CALC'
 486       common /srutu/ icall
 487 c     double precision rrsave(maxdim)
 488       logical lprn
 489       evdw=0.0D0
 490 c     print *,'Entering EBP nnt=',nnt,' nct=',nct,' expon=',expon
 491       evdw=0.0D0
 492 c     if (icall.eq.0) then
 493 c       lprn=.true.
 494 c     else
 495         lprn=.false.
 496 c     endif
 497       ind=0
 498       do i=iatsc_s,iatsc_e
 499         itypi=itype(i)
 500         if (itypi.eq.21) cycle
 501         itypi1=itype(i+1)
 502         xi=c(1,nres+i)
 503         yi=c(2,nres+i)
 504         zi=c(3,nres+i)
 505         dxi=dc_norm(1,nres+i)
 506         dyi=dc_norm(2,nres+i)
 507         dzi=dc_norm(3,nres+i)
 508 c        dsci_inv=dsc_inv(itypi)
 509         dsci_inv=vbld_inv(i+nres)
 510 C
 511 C Calculate SC interaction energy.
 512 C
 513         do iint=1,nint_gr(i)
 514           do j=istart(i,iint),iend(i,iint)
 515             ind=ind+1
 516             itypj=itype(j)
 517             if (itypj.eq.21) cycle
 518 c            dscj_inv=dsc_inv(itypj)
 519             dscj_inv=vbld_inv(j+nres)
 520             chi1=chi(itypi,itypj)
 521             chi2=chi(itypj,itypi)
 522             chi12=chi1*chi2
 523             chip1=chip(itypi)
 524             chip2=chip(itypj)
 525             chip12=chip1*chip2
 526             alf1=alp(itypi)
 527             alf2=alp(itypj)
 528             alf12=0.5D0*(alf1+alf2)
 529             xj=c(1,nres+j)-xi
 530             yj=c(2,nres+j)-yi
 531             zj=c(3,nres+j)-zi
 532             dxj=dc_norm(1,nres+j)
 533             dyj=dc_norm(2,nres+j)
 534             dzj=dc_norm(3,nres+j)
 535             rrij=1.0D0/(xj*xj+yj*yj+zj*zj)
 536             rij=dsqrt(rrij)
 537             sss=sscale(1.0d0/(rij*sigmaii(itypi,itypj)))
 538
 539             if (sss.gt.0.0d0) then
 540
 541 C Calculate the angle-dependent terms of energy & contributions to derivatives.
 542               call sc_angular
 543 C Calculate whole angle-dependent part of epsilon and contributions
 544 C to its derivatives
 545               fac=(rrij*sigsq)**expon2
 546               e1=fac*fac*aa(itypi,itypj)
 547               e2=fac*bb(itypi,itypj)
 548               evdwij=eps1*eps2rt*eps3rt*(e1+e2)
 549               eps2der=evdwij*eps3rt
 550               eps3der=evdwij*eps2rt
 551               evdwij=evdwij*eps2rt*eps3rt
 552               evdw=evdw+evdwij*sss
 553               if (lprn) then
 554               sigm=dabs(aa(itypi,itypj)/bb(itypi,itypj))**(1.0D0/6.0D0)
 555               epsi=bb(itypi,itypj)**2/aa(itypi,itypj)
 556 cd              write (iout,'(2(a3,i3,2x),15(0pf7.3))')
 557 cd     &          restyp(itypi),i,restyp(itypj),j,
 558 cd     &          epsi,sigm,chi1,chi2,chip1,chip2,
 559 cd     &          eps1,eps2rt**2,eps3rt**2,1.0D0/dsqrt(sigsq),
 560 cd     &          om1,om2,om12,1.0D0/dsqrt(rrij),
 561 cd     &          evdwij
 562               endif
 563 C Calculate gradient components.
 564               e1=e1*eps1*eps2rt**2*eps3rt**2
 565               fac=-expon*(e1+evdwij)
 566               sigder=fac/sigsq
 567               fac=rrij*fac
 568 C Calculate radial part of the gradient
 569               gg(1)=xj*fac
 570               gg(2)=yj*fac
 571               gg(3)=zj*fac
 572 C Calculate the angular part of the gradient and sum add the contributions
 573 C to the appropriate components of the Cartesian gradient.
 574               call sc_grad_scale(sss)
 575             endif
 576           enddo      ! j
 577         enddo        ! iint
 578       enddo          ! i
 579 c     stop
 580       return
 581       end
 582 C-----------------------------------------------------------------------------
 583       subroutine egb_long(evdw)
 584 C
 585 C This subroutine calculates the interaction energy of nonbonded side chains
 586 C assuming the Gay-Berne potential of interaction.
 587 C
 588       implicit real*8 (a-h,o-z)
 589       include 'DIMENSIONS'
 590       include 'COMMON.GEO'
 591       include 'COMMON.VAR'
 592       include 'COMMON.LOCAL'
 593       include 'COMMON.CHAIN'
 594       include 'COMMON.DERIV'
 595       include 'COMMON.NAMES'
 596       include 'COMMON.INTERACT'
 597       include 'COMMON.IOUNITS'
 598       include 'COMMON.CALC'
 599       include 'COMMON.CONTROL'
 600       logical lprn
 601       evdw=0.0D0
 602 ccccc      energy_dec=.false.
 603 c     print *,'Entering EGB nnt=',nnt,' nct=',nct,' expon=',expon
 604       evdw=0.0D0
 605       lprn=.false.
 606 c     if (icall.eq.0) lprn=.false.
 607       ind=0
 608       do i=iatsc_s,iatsc_e
 609         itypi=itype(i)
 610         if (itypi.eq.21) cycle
 611         itypi1=itype(i+1)
 612         xi=c(1,nres+i)
 613         yi=c(2,nres+i)
 614         zi=c(3,nres+i)
 615         dxi=dc_norm(1,nres+i)
 616         dyi=dc_norm(2,nres+i)
 617         dzi=dc_norm(3,nres+i)
 618 c        dsci_inv=dsc_inv(itypi)
 619         dsci_inv=vbld_inv(i+nres)
 620 c        write (iout,*) "i",i,dsc_inv(itypi),dsci_inv,1.0d0/vbld(i+nres)
 621 c        write (iout,*) "dcnori",dxi*dxi+dyi*dyi+dzi*dzi
 622 C
 623 C Calculate SC interaction energy.
 624 C
 625         do iint=1,nint_gr(i)
 626           do j=istart(i,iint),iend(i,iint)
 627             ind=ind+1
 628             itypj=itype(j)
 629             if (itypj.eq.21) cycle
 630 c            dscj_inv=dsc_inv(itypj)
 631             dscj_inv=vbld_inv(j+nres)
 632 c            write (iout,*) "j",j,dsc_inv(itypj),dscj_inv,
 633 c     &       1.0d0/vbld(j+nres)
 634 c            write (iout,*) "i",i," j", j," itype",itype(i),itype(j)
 635             sig0ij=sigma(itypi,itypj)
 636             chi1=chi(itypi,itypj)
 637             chi2=chi(itypj,itypi)
 638             chi12=chi1*chi2
 639             chip1=chip(itypi)
 640             chip2=chip(itypj)
 641             chip12=chip1*chip2
 642             alf1=alp(itypi)
 643             alf2=alp(itypj)
 644             alf12=0.5D0*(alf1+alf2)
 645             xj=c(1,nres+j)-xi
 646             yj=c(2,nres+j)-yi
 647             zj=c(3,nres+j)-zi
 648             dxj=dc_norm(1,nres+j)
 649             dyj=dc_norm(2,nres+j)
 650             dzj=dc_norm(3,nres+j)
 651             rrij=1.0D0/(xj*xj+yj*yj+zj*zj)
 652             rij=dsqrt(rrij)
 653             sss=sscale(1.0d0/(rij*sigmaii(itypi,itypj)))
 654
 655             if (sss.lt.1.0d0) then
 656
 657 C Calculate angle-dependent terms of energy and contributions to their
 658 C derivatives.
 659               call sc_angular
 660               sigsq=1.0D0/sigsq
 661               sig=sig0ij*dsqrt(sigsq)
 662               rij_shift=1.0D0/rij-sig+sig0ij
 663 c for diagnostics; uncomment
 664 c              rij_shift=1.2*sig0ij
 665 C I hate to put IF's in the loops, but here don't have another choice!!!!
 666               if (rij_shift.le.0.0D0) then
 667                 evdw=1.0D20
 668 cd                write (iout,'(2(a3,i3,2x),17(0pf7.3))')
 669 cd     &          restyp(itypi),i,restyp(itypj),j,
 670 cd     &          rij_shift,1.0D0/rij,sig,sig0ij,sigsq,1-dsqrt(sigsq)
 671                 return
 672               endif
 673               sigder=-sig*sigsq
 674 c---------------------------------------------------------------
 675               rij_shift=1.0D0/rij_shift
 676               fac=rij_shift**expon
 677               e1=fac*fac*aa(itypi,itypj)
 678               e2=fac*bb(itypi,itypj)
 679               evdwij=eps1*eps2rt*eps3rt*(e1+e2)
 680               eps2der=evdwij*eps3rt
 681               eps3der=evdwij*eps2rt
 682 c              write (iout,*) "sigsq",sigsq," sig",sig," eps2rt",eps2rt,
 683 c     &        " eps3rt",eps3rt," eps1",eps1," e1",e1," e2",e2
 684               evdwij=evdwij*eps2rt*eps3rt
 685               evdw=evdw+evdwij*(1.0d0-sss)
 686               if (lprn) then
 687               sigm=dabs(aa(itypi,itypj)/bb(itypi,itypj))**(1.0D0/6.0D0)
 688               epsi=bb(itypi,itypj)**2/aa(itypi,itypj)
 689               write (iout,'(2(a3,i3,2x),17(0pf7.3))')
 690      &          restyp(itypi),i,restyp(itypj),j,
 691      &          epsi,sigm,chi1,chi2,chip1,chip2,
 692      &          eps1,eps2rt**2,eps3rt**2,sig,sig0ij,
 693      &          om1,om2,om12,1.0D0/rij,1.0D0/rij_shift,
 694      &          evdwij
 695               endif
 696
 697               if (energy_dec) write (iout,'(a6,2i5,0pf7.3)')
 698      &                        'evdw',i,j,evdwij
 699
 700 C Calculate gradient components.
 701               e1=e1*eps1*eps2rt**2*eps3rt**2
 702               fac=-expon*(e1+evdwij)*rij_shift
 703               sigder=fac*sigder
 704               fac=rij*fac
 705 c              fac=0.0d0
 706 C Calculate the radial part of the gradient
 707               gg(1)=xj*fac
 708               gg(2)=yj*fac
 709               gg(3)=zj*fac
 710 C Calculate angular part of the gradient.
 711               call sc_grad_scale(1.0d0-sss)
 712             endif
 713           enddo      ! j
 714         enddo        ! iint
 715       enddo          ! i
 716 c      write (iout,*) "Number of loop steps in EGB:",ind
 717 cccc      energy_dec=.false.
 718       return
 719       end
 720 C-----------------------------------------------------------------------------
 721       subroutine egb_short(evdw)
 722 C
 723 C This subroutine calculates the interaction energy of nonbonded side chains
 724 C assuming the Gay-Berne potential of interaction.
 725 C
 726       implicit real*8 (a-h,o-z)
 727       include 'DIMENSIONS'
 728       include 'COMMON.GEO'
 729       include 'COMMON.VAR'
 730       include 'COMMON.LOCAL'
 731       include 'COMMON.CHAIN'
 732       include 'COMMON.DERIV'
 733       include 'COMMON.NAMES'
 734       include 'COMMON.INTERACT'
 735       include 'COMMON.IOUNITS'
 736       include 'COMMON.CALC'
 737       include 'COMMON.CONTROL'
 738       logical lprn
 739       evdw=0.0D0
 740 ccccc      energy_dec=.false.
 741 c     print *,'Entering EGB nnt=',nnt,' nct=',nct,' expon=',expon
 742       evdw=0.0D0
 743       lprn=.false.
 744 c     if (icall.eq.0) lprn=.false.
 745       ind=0
 746       do i=iatsc_s,iatsc_e
 747         itypi=itype(i)
 748         if (itypi.eq.21) cycle
 749         itypi1=itype(i+1)
 750         xi=c(1,nres+i)
 751         yi=c(2,nres+i)
 752         zi=c(3,nres+i)
 753         dxi=dc_norm(1,nres+i)
 754         dyi=dc_norm(2,nres+i)
 755         dzi=dc_norm(3,nres+i)
 756 c        dsci_inv=dsc_inv(itypi)
 757         dsci_inv=vbld_inv(i+nres)
 758 c        write (iout,*) "i",i,dsc_inv(itypi),dsci_inv,1.0d0/vbld(i+nres)
 759 c        write (iout,*) "dcnori",dxi*dxi+dyi*dyi+dzi*dzi
 760 C
 761 C Calculate SC interaction energy.
 762 C
 763         do iint=1,nint_gr(i)
 764           do j=istart(i,iint),iend(i,iint)
 765             ind=ind+1
 766             itypj=itype(j)
 767             if (itypj.eq.21) cycle
 768 c            dscj_inv=dsc_inv(itypj)
 769             dscj_inv=vbld_inv(j+nres)
 770 c            write (iout,*) "j",j,dsc_inv(itypj),dscj_inv,
 771 c     &       1.0d0/vbld(j+nres)
 772 c            write (iout,*) "i",i," j", j," itype",itype(i),itype(j)
 773             sig0ij=sigma(itypi,itypj)
 774             chi1=chi(itypi,itypj)
 775             chi2=chi(itypj,itypi)
 776             chi12=chi1*chi2
 777             chip1=chip(itypi)
 778             chip2=chip(itypj)
 779             chip12=chip1*chip2
 780             alf1=alp(itypi)
 781             alf2=alp(itypj)
 782             alf12=0.5D0*(alf1+alf2)
 783             xj=c(1,nres+j)-xi
 784             yj=c(2,nres+j)-yi
 785             zj=c(3,nres+j)-zi
 786             dxj=dc_norm(1,nres+j)
 787             dyj=dc_norm(2,nres+j)
 788             dzj=dc_norm(3,nres+j)
 789             rrij=1.0D0/(xj*xj+yj*yj+zj*zj)
 790             rij=dsqrt(rrij)
 791             sss=sscale(1.0d0/(rij*sigmaii(itypi,itypj)))
 792
 793             if (sss.gt.0.0d0) then
 794
 795 C Calculate angle-dependent terms of energy and contributions to their
 796 C derivatives.
 797               call sc_angular
 798               sigsq=1.0D0/sigsq
 799               sig=sig0ij*dsqrt(sigsq)
 800               rij_shift=1.0D0/rij-sig+sig0ij
 801 c for diagnostics; uncomment
 802 c              rij_shift=1.2*sig0ij
 803 C I hate to put IF's in the loops, but here don't have another choice!!!!
 804               if (rij_shift.le.0.0D0) then
 805                 evdw=1.0D20
 806 cd                write (iout,'(2(a3,i3,2x),17(0pf7.3))')
 807 cd     &          restyp(itypi),i,restyp(itypj),j,
 808 cd     &          rij_shift,1.0D0/rij,sig,sig0ij,sigsq,1-dsqrt(sigsq)
 809                 return
 810               endif
 811               sigder=-sig*sigsq
 812 c---------------------------------------------------------------
 813               rij_shift=1.0D0/rij_shift
 814               fac=rij_shift**expon
 815               e1=fac*fac*aa(itypi,itypj)
 816               e2=fac*bb(itypi,itypj)
 817               evdwij=eps1*eps2rt*eps3rt*(e1+e2)
 818               eps2der=evdwij*eps3rt
 819               eps3der=evdwij*eps2rt
 820 c              write (iout,*) "sigsq",sigsq," sig",sig," eps2rt",eps2rt,
 821 c     &        " eps3rt",eps3rt," eps1",eps1," e1",e1," e2",e2
 822               evdwij=evdwij*eps2rt*eps3rt
 823               evdw=evdw+evdwij*sss
 824               if (lprn) then
 825               sigm=dabs(aa(itypi,itypj)/bb(itypi,itypj))**(1.0D0/6.0D0)
 826               epsi=bb(itypi,itypj)**2/aa(itypi,itypj)
 827               write (iout,'(2(a3,i3,2x),17(0pf7.3))')
 828      &          restyp(itypi),i,restyp(itypj),j,
 829      &          epsi,sigm,chi1,chi2,chip1,chip2,
 830      &          eps1,eps2rt**2,eps3rt**2,sig,sig0ij,
 831      &          om1,om2,om12,1.0D0/rij,1.0D0/rij_shift,
 832      &          evdwij
 833               endif
 834
 835               if (energy_dec) write (iout,'(a6,2i5,0pf7.3)')
 836      &                        'evdw',i,j,evdwij
 837
 838 C Calculate gradient components.
 839               e1=e1*eps1*eps2rt**2*eps3rt**2
 840               fac=-expon*(e1+evdwij)*rij_shift
 841               sigder=fac*sigder
 842               fac=rij*fac
 843 c              fac=0.0d0
 844 C Calculate the radial part of the gradient
 845               gg(1)=xj*fac
 846               gg(2)=yj*fac
 847               gg(3)=zj*fac
 848 C Calculate angular part of the gradient.
 849               call sc_grad_scale(sss)
 850             endif
 851           enddo      ! j
 852         enddo        ! iint
 853       enddo          ! i
 854 c      write (iout,*) "Number of loop steps in EGB:",ind
 855 cccc      energy_dec=.false.
 856       return
 857       end
 858 C-----------------------------------------------------------------------------
 859       subroutine egbv_long(evdw)
 860 C
 861 C This subroutine calculates the interaction energy of nonbonded side chains
 862 C assuming the Gay-Berne-Vorobjev potential of interaction.
 863 C
 864       implicit real*8 (a-h,o-z)
 865       include 'DIMENSIONS'
 866       include 'COMMON.GEO'
 867       include 'COMMON.VAR'
 868       include 'COMMON.LOCAL'
 869       include 'COMMON.CHAIN'
 870       include 'COMMON.DERIV'
 871       include 'COMMON.NAMES'
 872       include 'COMMON.INTERACT'
 873       include 'COMMON.IOUNITS'
 874       include 'COMMON.CALC'
 875       common /srutu/ icall
 876       logical lprn
 877       evdw=0.0D0
 878 c     print *,'Entering EGB nnt=',nnt,' nct=',nct,' expon=',expon
 879       evdw=0.0D0
 880       lprn=.false.
 881 c     if (icall.eq.0) lprn=.true.
 882       ind=0
 883       do i=iatsc_s,iatsc_e
 884         itypi=itype(i)
 885         if (itypi.eq.21) cycle
 886         itypi1=itype(i+1)
 887         xi=c(1,nres+i)
 888         yi=c(2,nres+i)
 889         zi=c(3,nres+i)
 890         dxi=dc_norm(1,nres+i)
 891         dyi=dc_norm(2,nres+i)
 892         dzi=dc_norm(3,nres+i)
 893 c        dsci_inv=dsc_inv(itypi)
 894         dsci_inv=vbld_inv(i+nres)
 895 C
 896 C Calculate SC interaction energy.
 897 C
 898         do iint=1,nint_gr(i)
 899           do j=istart(i,iint),iend(i,iint)
 900             ind=ind+1
 901             itypj=itype(j)
 902             if (itypj.eq.21) cycle
 903 c            dscj_inv=dsc_inv(itypj)
 904             dscj_inv=vbld_inv(j+nres)
 905             sig0ij=sigma(itypi,itypj)
 906             r0ij=r0(itypi,itypj)
 907             chi1=chi(itypi,itypj)
 908             chi2=chi(itypj,itypi)
 909             chi12=chi1*chi2
 910             chip1=chip(itypi)
 911             chip2=chip(itypj)
 912             chip12=chip1*chip2
 913             alf1=alp(itypi)
 914             alf2=alp(itypj)
 915             alf12=0.5D0*(alf1+alf2)
 916             xj=c(1,nres+j)-xi
 917             yj=c(2,nres+j)-yi
 918             zj=c(3,nres+j)-zi
 919             dxj=dc_norm(1,nres+j)
 920             dyj=dc_norm(2,nres+j)
 921             dzj=dc_norm(3,nres+j)
 922             rrij=1.0D0/(xj*xj+yj*yj+zj*zj)
 923             rij=dsqrt(rrij)
 924
 925             sss=sscale(1.0d0/(rij*sigmaii(itypi,itypj)))
 926
 927             if (sss.lt.1.0d0) then
 928
 929 C Calculate angle-dependent terms of energy and contributions to their
 930 C derivatives.
 931               call sc_angular
 932               sigsq=1.0D0/sigsq
 933               sig=sig0ij*dsqrt(sigsq)
 934               rij_shift=1.0D0/rij-sig+r0ij
 935 C I hate to put IF's in the loops, but here don't have another choice!!!!
 936               if (rij_shift.le.0.0D0) then
 937                 evdw=1.0D20
 938                 return
 939               endif
 940               sigder=-sig*sigsq
 941 c---------------------------------------------------------------
 942               rij_shift=1.0D0/rij_shift
 943               fac=rij_shift**expon
 944               e1=fac*fac*aa(itypi,itypj)
 945               e2=fac*bb(itypi,itypj)
 946               evdwij=eps1*eps2rt*eps3rt*(e1+e2)
 947               eps2der=evdwij*eps3rt
 948               eps3der=evdwij*eps2rt
 949               fac_augm=rrij**expon
 950               e_augm=augm(itypi,itypj)*fac_augm
 951               evdwij=evdwij*eps2rt*eps3rt
 952               evdw=evdw+(evdwij+e_augm)*(1.0d0-sss)
 953               if (lprn) then
 954               sigm=dabs(aa(itypi,itypj)/bb(itypi,itypj))**(1.0D0/6.0D0)
 955               epsi=bb(itypi,itypj)**2/aa(itypi,itypj)
 956               write (iout,'(2(a3,i3,2x),17(0pf7.3))')
 957      &          restyp(itypi),i,restyp(itypj),j,
 958      &          epsi,sigm,sig,(augm(itypi,itypj)/epsi)**(1.0D0/12.0D0),
 959      &          chi1,chi2,chip1,chip2,
 960      &          eps1,eps2rt**2,eps3rt**2,
 961      &          om1,om2,om12,1.0D0/rij,1.0D0/rij_shift,
 962      &          evdwij+e_augm
 963               endif
 964 C Calculate gradient components.
 965               e1=e1*eps1*eps2rt**2*eps3rt**2
 966               fac=-expon*(e1+evdwij)*rij_shift
 967               sigder=fac*sigder
 968               fac=rij*fac-2*expon*rrij*e_augm
 969 C Calculate the radial part of the gradient
 970               gg(1)=xj*fac
 971               gg(2)=yj*fac
 972               gg(3)=zj*fac
 973 C Calculate angular part of the gradient.
 974               call sc_grad_scale(1.0d0-sss)
 975             endif
 976           enddo      ! j
 977         enddo        ! iint
 978       enddo          ! i
 979       end
 980 C-----------------------------------------------------------------------------
 981       subroutine egbv_short(evdw)
 982 C
 983 C This subroutine calculates the interaction energy of nonbonded side chains
 984 C assuming the Gay-Berne-Vorobjev potential of interaction.
 985 C
 986       implicit real*8 (a-h,o-z)
 987       include 'DIMENSIONS'
 988       include 'COMMON.GEO'
 989       include 'COMMON.VAR'
 990       include 'COMMON.LOCAL'
 991       include 'COMMON.CHAIN'
 992       include 'COMMON.DERIV'
 993       include 'COMMON.NAMES'
 994       include 'COMMON.INTERACT'
 995       include 'COMMON.IOUNITS'
 996       include 'COMMON.CALC'
 997       common /srutu/ icall
 998       logical lprn
 999       evdw=0.0D0
1000 c     print *,'Entering EGB nnt=',nnt,' nct=',nct,' expon=',expon
1001       evdw=0.0D0
1002       lprn=.false.
1003 c     if (icall.eq.0) lprn=.true.
1004       ind=0
1005       do i=iatsc_s,iatsc_e
1006         itypi=itype(i)
1007         if (itypi.eq.21) cycle
1008         itypi1=itype(i+1)
1009         xi=c(1,nres+i)
1010         yi=c(2,nres+i)
1011         zi=c(3,nres+i)
1012         dxi=dc_norm(1,nres+i)
1013         dyi=dc_norm(2,nres+i)
1014         dzi=dc_norm(3,nres+i)
1015 c        dsci_inv=dsc_inv(itypi)
1016         dsci_inv=vbld_inv(i+nres)
1017 C
1018 C Calculate SC interaction energy.
1019 C
1020         do iint=1,nint_gr(i)
1021           do j=istart(i,iint),iend(i,iint)
1022             ind=ind+1
1023             itypj=itype(j)
1024             if (itypj.eq.21) cycle
1025 c            dscj_inv=dsc_inv(itypj)
1026             dscj_inv=vbld_inv(j+nres)
1027             sig0ij=sigma(itypi,itypj)
1028             r0ij=r0(itypi,itypj)
1029             chi1=chi(itypi,itypj)
1030             chi2=chi(itypj,itypi)
1031             chi12=chi1*chi2
1032             chip1=chip(itypi)
1033             chip2=chip(itypj)
1034             chip12=chip1*chip2
1035             alf1=alp(itypi)
1036             alf2=alp(itypj)
1037             alf12=0.5D0*(alf1+alf2)
1038             xj=c(1,nres+j)-xi
1039             yj=c(2,nres+j)-yi
1040             zj=c(3,nres+j)-zi
1041             dxj=dc_norm(1,nres+j)
1042             dyj=dc_norm(2,nres+j)
1043             dzj=dc_norm(3,nres+j)
1044             rrij=1.0D0/(xj*xj+yj*yj+zj*zj)
1045             rij=dsqrt(rrij)
1046
1047             sss=sscale(1.0d0/(rij*sigmaii(itypi,itypj)))
1048
1049             if (sss.gt.0.0d0) then
1050
1051 C Calculate angle-dependent terms of energy and contributions to their
1052 C derivatives.
1053               call sc_angular
1054               sigsq=1.0D0/sigsq
1055               sig=sig0ij*dsqrt(sigsq)
1056               rij_shift=1.0D0/rij-sig+r0ij
1057 C I hate to put IF's in the loops, but here don't have another choice!!!!
1058               if (rij_shift.le.0.0D0) then
1059                 evdw=1.0D20
1060                 return
1061               endif
1062               sigder=-sig*sigsq
1063 c---------------------------------------------------------------
1064               rij_shift=1.0D0/rij_shift
1065               fac=rij_shift**expon
1066               e1=fac*fac*aa(itypi,itypj)
1067               e2=fac*bb(itypi,itypj)
1068               evdwij=eps1*eps2rt*eps3rt*(e1+e2)
1069               eps2der=evdwij*eps3rt
1070               eps3der=evdwij*eps2rt
1071               fac_augm=rrij**expon
1072               e_augm=augm(itypi,itypj)*fac_augm
1073               evdwij=evdwij*eps2rt*eps3rt
1074               evdw=evdw+(evdwij+e_augm)*sss
1075               if (lprn) then
1076               sigm=dabs(aa(itypi,itypj)/bb(itypi,itypj))**(1.0D0/6.0D0)
1077               epsi=bb(itypi,itypj)**2/aa(itypi,itypj)
1078               write (iout,'(2(a3,i3,2x),17(0pf7.3))')
1079      &          restyp(itypi),i,restyp(itypj),j,
1080      &          epsi,sigm,sig,(augm(itypi,itypj)/epsi)**(1.0D0/12.0D0),
1081      &          chi1,chi2,chip1,chip2,
1082      &          eps1,eps2rt**2,eps3rt**2,
1083      &          om1,om2,om12,1.0D0/rij,1.0D0/rij_shift,
1084      &          evdwij+e_augm
1085               endif
1086 C Calculate gradient components.
1087               e1=e1*eps1*eps2rt**2*eps3rt**2
1088               fac=-expon*(e1+evdwij)*rij_shift
1089               sigder=fac*sigder
1090               fac=rij*fac-2*expon*rrij*e_augm
1091 C Calculate the radial part of the gradient
1092               gg(1)=xj*fac
1093               gg(2)=yj*fac
1094               gg(3)=zj*fac
1095 C Calculate angular part of the gradient.
1096               call sc_grad_scale(sss)
1097             endif
1098           enddo      ! j
1099         enddo        ! iint
1100       enddo          ! i
1101       end
1102 C----------------------------------------------------------------------------
1103       subroutine sc_grad_scale(scalfac)
1104       implicit real*8 (a-h,o-z)
1105       include 'DIMENSIONS'
1106       include 'COMMON.CHAIN'
1107       include 'COMMON.DERIV'
1108       include 'COMMON.CALC'
1109       include 'COMMON.IOUNITS'
1110       double precision dcosom1(3),dcosom2(3)
1111       double precision scalfac
1112       eom1=eps2der*eps2rt_om1-2.0D0*alf1*eps3der+sigder*sigsq_om1
1113       eom2=eps2der*eps2rt_om2+2.0D0*alf2*eps3der+sigder*sigsq_om2
1114       eom12=evdwij*eps1_om12+eps2der*eps2rt_om12
1115      &     -2.0D0*alf12*eps3der+sigder*sigsq_om12
1116 c diagnostics only
1117 c      eom1=0.0d0
1118 c      eom2=0.0d0
1119 c      eom12=evdwij*eps1_om12
1120 c end diagnostics
1121 c      write (iout,*) "eps2der",eps2der," eps3der",eps3der,
1122 c     &  " sigder",sigder
1123 c      write (iout,*) "eps1_om12",eps1_om12," eps2rt_om12",eps2rt_om12
1124 c      write (iout,*) "eom1",eom1," eom2",eom2," eom12",eom12
1125       do k=1,3
1126         dcosom1(k)=rij*(dc_norm(k,nres+i)-om1*erij(k))
1127         dcosom2(k)=rij*(dc_norm(k,nres+j)-om2*erij(k))
1128       enddo
1129       do k=1,3
1130         gg(k)=(gg(k)+eom1*dcosom1(k)+eom2*dcosom2(k))*scalfac
1131       enddo
1132 c      write (iout,*) "gg",(gg(k),k=1,3)
1133       do k=1,3
1134         gvdwx(k,i)=gvdwx(k,i)-gg(k)
1135      &            +(eom12*(dc_norm(k,nres+j)-om12*dc_norm(k,nres+i))
1136      &          +eom1*(erij(k)-om1*dc_norm(k,nres+i)))*dsci_inv*scalfac
1137         gvdwx(k,j)=gvdwx(k,j)+gg(k)
1138      &            +(eom12*(dc_norm(k,nres+i)-om12*dc_norm(k,nres+j))
1139      &          +eom2*(erij(k)-om2*dc_norm(k,nres+j)))*dscj_inv*scalfac
1140 c        write (iout,*)(eom12*(dc_norm(k,nres+j)-om12*dc_norm(k,nres+i))
1141 c     &            +eom1*(erij(k)-om1*dc_norm(k,nres+i)))*dsci_inv
1142 c        write (iout,*)(eom12*(dc_norm(k,nres+i)-om12*dc_norm(k,nres+j))
1143 c     &            +eom2*(erij(k)-om2*dc_norm(k,nres+j)))*dscj_inv
1144       enddo
1145 C
1146 C Calculate the components of the gradient in DC and X
1147 C
1148       do l=1,3
1149         gvdwc(l,i)=gvdwc(l,i)-gg(l)
1150         gvdwc(l,j)=gvdwc(l,j)+gg(l)
1151       enddo
1152       return
1153       end
1154 C--------------------------------------------------------------------------
1155       subroutine eelec_scale(ees,evdw1,eel_loc,eello_turn3,eello_turn4)
1156 C
1157 C This subroutine calculates the average interaction energy and its gradient
1158 C in the virtual-bond vectors between non-adjacent peptide groups, based on
1159 C the potential described in Liwo et al., Protein Sci., 1993, 2, 1715.
1160 C The potential depends both on the distance of peptide-group centers and on
1161 C the orientation of the CA-CA virtual bonds.
1162 C
1163       implicit real*8 (a-h,o-z)
1164 #ifdef MPI
1165       include 'mpif.h'
1166 #endif
1167       include 'DIMENSIONS'
1168       include 'COMMON.CONTROL'
1169       include 'COMMON.SETUP'
1170       include 'COMMON.IOUNITS'
1171       include 'COMMON.GEO'
1172       include 'COMMON.VAR'
1173       include 'COMMON.LOCAL'
1174       include 'COMMON.CHAIN'
1175       include 'COMMON.DERIV'
1176       include 'COMMON.INTERACT'
1177       include 'COMMON.CONTACTS'
1178       include 'COMMON.TORSION'
1179       include 'COMMON.VECTORS'
1180       include 'COMMON.FFIELD'
1181       include 'COMMON.TIME1'
1182       dimension ggg(3),gggp(3),gggm(3),erij(3),dcosb(3),dcosg(3),
1183      &          erder(3,3),uryg(3,3),urzg(3,3),vryg(3,3),vrzg(3,3)
1184       double precision acipa(2,2),agg(3,4),aggi(3,4),aggi1(3,4),
1185      &    aggj(3,4),aggj1(3,4),a_temp(2,2),muij(4)
1186       common /locel/ a_temp,agg,aggi,aggi1,aggj,aggj1,a22,a23,a32,a33,
1187      &    dxi,dyi,dzi,dx_normi,dy_normi,dz_normi,xmedi,ymedi,zmedi,
1188      &    num_conti,j1,j2
1189 c 4/26/02 - AL scaling factor for 1,4 repulsive VDW interactions
1190 #ifdef MOMENT
1191       double precision scal_el /1.0d0/
1192 #else
1193       double precision scal_el /0.5d0/
1194 #endif
1195 C 12/13/98
1196 C 13-go grudnia roku pamietnego...
1197       double precision unmat(3,3) /1.0d0,0.0d0,0.0d0,
1198      &                   0.0d0,1.0d0,0.0d0,
1199      &                   0.0d0,0.0d0,1.0d0/
1200 cd      write(iout,*) 'In EELEC'
1201 cd      do i=1,nloctyp
1202 cd        write(iout,*) 'Type',i
1203 cd        write(iout,*) 'B1',B1(:,i)
1204 cd        write(iout,*) 'B2',B2(:,i)
1205 cd        write(iout,*) 'CC',CC(:,:,i)
1206 cd        write(iout,*) 'DD',DD(:,:,i)
1207 cd        write(iout,*) 'EE',EE(:,:,i)
1208 cd      enddo
1209 cd      call check_vecgrad
1210 cd      stop
1211       if (icheckgrad.eq.1) then
1212         do i=1,nres-1
1213           fac=1.0d0/dsqrt(scalar(dc(1,i),dc(1,i)))
1214           do k=1,3
1215             dc_norm(k,i)=dc(k,i)*fac
1216           enddo
1217 c          write (iout,*) 'i',i,' fac',fac
1218         enddo
1219       endif
1220       if (wel_loc.gt.0.0d0 .or. wcorr4.gt.0.0d0 .or. wcorr5.gt.0.0d0
1221      &    .or. wcorr6.gt.0.0d0 .or. wturn3.gt.0.0d0 .or.
1222      &    wturn4.gt.0.0d0 .or. wturn6.gt.0.0d0) then
1223 c        call vec_and_deriv
1224 #ifdef TIMING
1225         time01=MPI_Wtime()
1226 #endif
1227         call set_matrices
1228 #ifdef TIMING
1229         time_mat=time_mat+MPI_Wtime()-time01
1230 #endif
1231       endif
1232 cd      do i=1,nres-1
1233 cd        write (iout,*) 'i=',i
1234 cd        do k=1,3
1235 cd        write (iout,'(i5,2f10.5)') k,uy(k,i),uz(k,i)
1236 cd        enddo
1237 cd        do k=1,3
1238 cd          write (iout,'(f10.5,2x,3f10.5,2x,3f10.5)')
1239 cd     &     uz(k,i),(uzgrad(k,l,1,i),l=1,3),(uzgrad(k,l,2,i),l=1,3)
1240 cd        enddo
1241 cd      enddo
1242       t_eelecij=0.0d0
1243       ees=0.0D0
1244       evdw1=0.0D0
1245       eel_loc=0.0d0
1246       eello_turn3=0.0d0
1247       eello_turn4=0.0d0
1248       ind=0
1249       do i=1,nres
1250         num_cont_hb(i)=0
1251       enddo
1252 cd      print '(a)','Enter EELEC'
1253 cd      write (iout,*) 'iatel_s=',iatel_s,' iatel_e=',iatel_e
1254       do i=1,nres
1255         gel_loc_loc(i)=0.0d0
1256         gcorr_loc(i)=0.0d0
1257       enddo
1258 c
1259 c
1260 c 9/27/08 AL Split the interaction loop to ensure load balancing of turn terms
1261 C
1262 C Loop over i,i+2 and i,i+3 pairs of the peptide groups
1263 C
1264       do i=iturn3_start,iturn3_end
1265         if (itype(i).eq.21 .or. itype(i+1).eq.21
1266      &  .or. itype(i+2).eq.21 .or. itype(i+3).eq.21) cycle
1267         dxi=dc(1,i)
1268         dyi=dc(2,i)
1269         dzi=dc(3,i)
1270         dx_normi=dc_norm(1,i)
1271         dy_normi=dc_norm(2,i)
1272         dz_normi=dc_norm(3,i)
1273         xmedi=c(1,i)+0.5d0*dxi
1274         ymedi=c(2,i)+0.5d0*dyi
1275         zmedi=c(3,i)+0.5d0*dzi
1276         num_conti=0
1277         call eelecij_scale(i,i+2,ees,evdw1,eel_loc)
1278         if (wturn3.gt.0.0d0) call eturn3(i,eello_turn3)
1279         num_cont_hb(i)=num_conti
1280       enddo
1281       do i=iturn4_start,iturn4_end
1282         if (itype(i).eq.21 .or. itype(i+1).eq.21
1283      &    .or. itype(i+3).eq.21
1284      &    .or. itype(i+4).eq.21) cycle
1285         dxi=dc(1,i)
1286         dyi=dc(2,i)
1287         dzi=dc(3,i)
1288         dx_normi=dc_norm(1,i)
1289         dy_normi=dc_norm(2,i)
1290         dz_normi=dc_norm(3,i)
1291         xmedi=c(1,i)+0.5d0*dxi
1292         ymedi=c(2,i)+0.5d0*dyi
1293         zmedi=c(3,i)+0.5d0*dzi
1294         num_conti=num_cont_hb(i)
1295         call eelecij_scale(i,i+3,ees,evdw1,eel_loc)
1296         if (wturn4.gt.0.0d0 .and. itype(i+2).ne.21)
1297      &    call eturn4(i,eello_turn4)
1298         num_cont_hb(i)=num_conti
1299       enddo   ! i
1300 c
1301 c Loop over all pairs of interacting peptide groups except i,i+2 and i,i+3
1302 c
1303       do i=iatel_s,iatel_e
1304         if (itype(i).eq.21 .or. itype(i+1).eq.21) cycle
1305         dxi=dc(1,i)
1306         dyi=dc(2,i)
1307         dzi=dc(3,i)
1308         dx_normi=dc_norm(1,i)
1309         dy_normi=dc_norm(2,i)
1310         dz_normi=dc_norm(3,i)
1311         xmedi=c(1,i)+0.5d0*dxi
1312         ymedi=c(2,i)+0.5d0*dyi
1313         zmedi=c(3,i)+0.5d0*dzi
1314 c        write (iout,*) 'i',i,' ielstart',ielstart(i),' ielend',ielend(i)
1315         num_conti=num_cont_hb(i)
1316         do j=ielstart(i),ielend(i)
1317           if (itype(j).eq.21 .or. itype(j+1).eq.21) cycle
1318           call eelecij_scale(i,j,ees,evdw1,eel_loc)
1319         enddo ! j
1320         num_cont_hb(i)=num_conti
1321       enddo   ! i
1322 c      write (iout,*) "Number of loop steps in EELEC:",ind
1323 cd      do i=1,nres
1324 cd        write (iout,'(i3,3f10.5,5x,3f10.5)')
1325 cd     &     i,(gel_loc(k,i),k=1,3),gel_loc_loc(i)
1326 cd      enddo
1327 c 12/7/99 Adam eello_turn3 will be considered as a separate energy term
1328 ccc      eel_loc=eel_loc+eello_turn3
1329 cd      print *,"Processor",fg_rank," t_eelecij",t_eelecij
1330       return
1331       end
1332 C-------------------------------------------------------------------------------
1333       subroutine eelecij_scale(i,j,ees,evdw1,eel_loc)
1334       implicit real*8 (a-h,o-z)
1335       include 'DIMENSIONS'
1336 #ifdef MPI
1337       include "mpif.h"
1338 #endif
1339       include 'COMMON.CONTROL'
1340       include 'COMMON.IOUNITS'
1341       include 'COMMON.GEO'
1342       include 'COMMON.VAR'
1343       include 'COMMON.LOCAL'
1344       include 'COMMON.CHAIN'
1345       include 'COMMON.DERIV'
1346       include 'COMMON.INTERACT'
1347       include 'COMMON.CONTACTS'
1348       include 'COMMON.TORSION'
1349       include 'COMMON.VECTORS'
1350       include 'COMMON.FFIELD'
1351       include 'COMMON.TIME1'
1352       dimension ggg(3),gggp(3),gggm(3),erij(3),dcosb(3),dcosg(3),
1353      &          erder(3,3),uryg(3,3),urzg(3,3),vryg(3,3),vrzg(3,3)
1354       double precision acipa(2,2),agg(3,4),aggi(3,4),aggi1(3,4),
1355      &    aggj(3,4),aggj1(3,4),a_temp(2,2),muij(4)
1356       common /locel/ a_temp,agg,aggi,aggi1,aggj,aggj1,a22,a23,a32,a33,
1357      &    dxi,dyi,dzi,dx_normi,dy_normi,dz_normi,xmedi,ymedi,zmedi,
1358      &    num_conti,j1,j2
1359 c 4/26/02 - AL scaling factor for 1,4 repulsive VDW interactions
1360 #ifdef MOMENT
1361       double precision scal_el /1.0d0/
1362 #else
1363       double precision scal_el /0.5d0/
1364 #endif
1365 C 12/13/98
1366 C 13-go grudnia roku pamietnego...
1367       double precision unmat(3,3) /1.0d0,0.0d0,0.0d0,
1368      &                   0.0d0,1.0d0,0.0d0,
1369      &                   0.0d0,0.0d0,1.0d0/
1370 c          time00=MPI_Wtime()
1371 cd      write (iout,*) "eelecij",i,j
1372           ind=ind+1
1373           iteli=itel(i)
1374           itelj=itel(j)
1375           if (j.eq.i+2 .and. itelj.eq.2) iteli=2
1376           aaa=app(iteli,itelj)
1377           bbb=bpp(iteli,itelj)
1378           ael6i=ael6(iteli,itelj)
1379           ael3i=ael3(iteli,itelj)
1380           dxj=dc(1,j)
1381           dyj=dc(2,j)
1382           dzj=dc(3,j)
1383           dx_normj=dc_norm(1,j)
1384           dy_normj=dc_norm(2,j)
1385           dz_normj=dc_norm(3,j)
1386           xj=c(1,j)+0.5D0*dxj-xmedi
1387           yj=c(2,j)+0.5D0*dyj-ymedi
1388           zj=c(3,j)+0.5D0*dzj-zmedi
1389           rij=xj*xj+yj*yj+zj*zj
1390           rrmij=1.0D0/rij
1391           rij=dsqrt(rij)
1392           rmij=1.0D0/rij
1393 c For extracting the short-range part of Evdwpp
1394           sss=sscale(rij/rpp(iteli,itelj))
1395
1396           r3ij=rrmij*rmij
1397           r6ij=r3ij*r3ij
1398           cosa=dx_normi*dx_normj+dy_normi*dy_normj+dz_normi*dz_normj
1399           cosb=(xj*dx_normi+yj*dy_normi+zj*dz_normi)*rmij
1400           cosg=(xj*dx_normj+yj*dy_normj+zj*dz_normj)*rmij
1401           fac=cosa-3.0D0*cosb*cosg
1402           ev1=aaa*r6ij*r6ij
1403 c 4/26/02 - AL scaling down 1,4 repulsive VDW interactions
1404           if (j.eq.i+2) ev1=scal_el*ev1
1405           ev2=bbb*r6ij
1406           fac3=ael6i*r6ij
1407           fac4=ael3i*r3ij
1408           evdwij=ev1+ev2
1409           el1=fac3*(4.0D0+fac*fac-3.0D0*(cosb*cosb+cosg*cosg))
1410           el2=fac4*fac
1411           eesij=el1+el2
1412 C 12/26/95 - for the evaluation of multi-body H-bonding interactions
1413           ees0ij=4.0D0+fac*fac-3.0D0*(cosb*cosb+cosg*cosg)
1414           ees=ees+eesij
1415           evdw1=evdw1+evdwij*(1.0d0-sss)
1416 cd          write(iout,'(2(2i3,2x),7(1pd12.4)/2(3(1pd12.4),5x)/)')
1417 cd     &      iteli,i,itelj,j,aaa,bbb,ael6i,ael3i,
1418 cd     &      1.0D0/dsqrt(rrmij),evdwij,eesij,
1419 cd     &      xmedi,ymedi,zmedi,xj,yj,zj
1420
1421           if (energy_dec) then
1422               write (iout,'(a6,2i5,0pf7.3,f7.3)') 'evdw1',i,j,evdwij,sss
1423               write (iout,'(a6,2i5,0pf7.3)') 'ees',i,j,eesij
1424           endif
1425
1426 C
1427 C Calculate contributions to the Cartesian gradient.
1428 C
1429 #ifdef SPLITELE
1430           facvdw=-6*rrmij*(ev1+evdwij)*(1.0d0-sss)
1431           facel=-3*rrmij*(el1+eesij)
1432           fac1=fac
1433           erij(1)=xj*rmij
1434           erij(2)=yj*rmij
1435           erij(3)=zj*rmij
1436 *
1437 * Radial derivatives. First process both termini of the fragment (i,j)
1438 *
1439           ggg(1)=facel*xj
1440           ggg(2)=facel*yj
1441           ggg(3)=facel*zj
1442 c          do k=1,3
1443 c            ghalf=0.5D0*ggg(k)
1444 c            gelc(k,i)=gelc(k,i)+ghalf
1445 c            gelc(k,j)=gelc(k,j)+ghalf
1446 c          enddo
1447 c 9/28/08 AL Gradient compotents will be summed only at the end
1448           do k=1,3
1449             gelc_long(k,j)=gelc_long(k,j)+ggg(k)
1450             gelc_long(k,i)=gelc_long(k,i)-ggg(k)
1451           enddo
1452 *
1453 * Loop over residues i+1 thru j-1.
1454 *
1455 cgrad          do k=i+1,j-1
1456 cgrad            do l=1,3
1457 cgrad              gelc(l,k)=gelc(l,k)+ggg(l)
1458 cgrad            enddo
1459 cgrad          enddo
1460           ggg(1)=facvdw*xj
1461           ggg(2)=facvdw*yj
1462           ggg(3)=facvdw*zj
1463 c          do k=1,3
1464 c            ghalf=0.5D0*ggg(k)
1465 c            gvdwpp(k,i)=gvdwpp(k,i)+ghalf
1466 c            gvdwpp(k,j)=gvdwpp(k,j)+ghalf
1467 c          enddo
1468 c 9/28/08 AL Gradient compotents will be summed only at the end
1469           do k=1,3
1470             gvdwpp(k,j)=gvdwpp(k,j)+ggg(k)
1471             gvdwpp(k,i)=gvdwpp(k,i)-ggg(k)
1472           enddo
1473 *
1474 * Loop over residues i+1 thru j-1.
1475 *
1476 cgrad          do k=i+1,j-1
1477 cgrad            do l=1,3
1478 cgrad              gvdwpp(l,k)=gvdwpp(l,k)+ggg(l)
1479 cgrad            enddo
1480 cgrad          enddo
1481 #else
1482           facvdw=ev1+evdwij*(1.0d0-sss)
1483           facel=el1+eesij
1484           fac1=fac
1485           fac=-3*rrmij*(facvdw+facvdw+facel)
1486           erij(1)=xj*rmij
1487           erij(2)=yj*rmij
1488           erij(3)=zj*rmij
1489 *
1490 * Radial derivatives. First process both termini of the fragment (i,j)
1491 *
1492           ggg(1)=fac*xj
1493           ggg(2)=fac*yj
1494           ggg(3)=fac*zj
1495 c          do k=1,3
1496 c            ghalf=0.5D0*ggg(k)
1497 c            gelc(k,i)=gelc(k,i)+ghalf
1498 c            gelc(k,j)=gelc(k,j)+ghalf
1499 c          enddo
1500 c 9/28/08 AL Gradient compotents will be summed only at the end
1501           do k=1,3
1502             gelc_long(k,j)=gelc(k,j)+ggg(k)
1503             gelc_long(k,i)=gelc(k,i)-ggg(k)
1504           enddo
1505 *
1506 * Loop over residues i+1 thru j-1.
1507 *
1508 cgrad          do k=i+1,j-1
1509 cgrad            do l=1,3
1510 cgrad              gelc(l,k)=gelc(l,k)+ggg(l)
1511 cgrad            enddo
1512 cgrad          enddo
1513 c 9/28/08 AL Gradient compotents will be summed only at the end
1514           ggg(1)=facvdw*xj
1515           ggg(2)=facvdw*yj
1516           ggg(3)=facvdw*zj
1517           do k=1,3
1518             gvdwpp(k,j)=gvdwpp(k,j)+ggg(k)
1519             gvdwpp(k,i)=gvdwpp(k,i)-ggg(k)
1520           enddo
1521 #endif
1522 *
1523 * Angular part
1524 *
1525           ecosa=2.0D0*fac3*fac1+fac4
1526           fac4=-3.0D0*fac4
1527           fac3=-6.0D0*fac3
1528           ecosb=(fac3*(fac1*cosg+cosb)+cosg*fac4)
1529           ecosg=(fac3*(fac1*cosb+cosg)+cosb*fac4)
1530           do k=1,3
1531             dcosb(k)=rmij*(dc_norm(k,i)-erij(k)*cosb)
1532             dcosg(k)=rmij*(dc_norm(k,j)-erij(k)*cosg)
1533           enddo
1534 cd        print '(2i3,2(3(1pd14.5),3x))',i,j,(dcosb(k),k=1,3),
1535 cd   &          (dcosg(k),k=1,3)
1536           do k=1,3
1537             ggg(k)=ecosb*dcosb(k)+ecosg*dcosg(k)
1538           enddo
1539 c          do k=1,3
1540 c            ghalf=0.5D0*ggg(k)
1541 c            gelc(k,i)=gelc(k,i)+ghalf
1542 c     &               +(ecosa*(dc_norm(k,j)-cosa*dc_norm(k,i))
1543 c     &               + ecosb*(erij(k)-cosb*dc_norm(k,i)))*vbld_inv(i+1)
1544 c            gelc(k,j)=gelc(k,j)+ghalf
1545 c     &               +(ecosa*(dc_norm(k,i)-cosa*dc_norm(k,j))
1546 c     &               + ecosg*(erij(k)-cosg*dc_norm(k,j)))*vbld_inv(j+1)
1547 c          enddo
1548 cgrad          do k=i+1,j-1
1549 cgrad            do l=1,3
1550 cgrad              gelc(l,k)=gelc(l,k)+ggg(l)
1551 cgrad            enddo
1552 cgrad          enddo
1553           do k=1,3
1554             gelc(k,i)=gelc(k,i)
1555      &               +(ecosa*(dc_norm(k,j)-cosa*dc_norm(k,i))
1556      &               + ecosb*(erij(k)-cosb*dc_norm(k,i)))*vbld_inv(i+1)
1557             gelc(k,j)=gelc(k,j)
1558      &               +(ecosa*(dc_norm(k,i)-cosa*dc_norm(k,j))
1559      &               + ecosg*(erij(k)-cosg*dc_norm(k,j)))*vbld_inv(j+1)
1560             gelc_long(k,j)=gelc_long(k,j)+ggg(k)
1561             gelc_long(k,i)=gelc_long(k,i)-ggg(k)
1562           enddo
1563           IF (wel_loc.gt.0.0d0 .or. wcorr4.gt.0.0d0 .or. wcorr5.gt.0.0d0
1564      &        .or. wcorr6.gt.0.0d0 .or. wturn3.gt.0.0d0
1565      &        .or. wturn4.gt.0.0d0 .or. wturn6.gt.0.0d0) THEN
1566 C
1567 C 9/25/99 Mixed third-order local-electrostatic terms. The local-interaction
1568 C   energy of a peptide unit is assumed in the form of a second-order
1569 C   Fourier series in the angles lambda1 and lambda2 (see Nishikawa et al.
1570 C   Macromolecules, 1974, 7, 797-806 for definition). This correlation terms
1571 C   are computed for EVERY pair of non-contiguous peptide groups.
1572 C
1573           if (j.lt.nres-1) then
1574             j1=j+1
1575             j2=j-1
1576           else
1577             j1=j-1
1578             j2=j-2
1579           endif
1580           kkk=0
1581           do k=1,2
1582             do l=1,2
1583               kkk=kkk+1
1584               muij(kkk)=mu(k,i)*mu(l,j)
1585             enddo
1586           enddo
1587 cd         write (iout,*) 'EELEC: i',i,' j',j
1588 cd          write (iout,*) 'j',j,' j1',j1,' j2',j2
1589 cd          write(iout,*) 'muij',muij
1590           ury=scalar(uy(1,i),erij)
1591           urz=scalar(uz(1,i),erij)
1592           vry=scalar(uy(1,j),erij)
1593           vrz=scalar(uz(1,j),erij)
1594           a22=scalar(uy(1,i),uy(1,j))-3*ury*vry
1595           a23=scalar(uy(1,i),uz(1,j))-3*ury*vrz
1596           a32=scalar(uz(1,i),uy(1,j))-3*urz*vry
1597           a33=scalar(uz(1,i),uz(1,j))-3*urz*vrz
1598           fac=dsqrt(-ael6i)*r3ij
1599           a22=a22*fac
1600           a23=a23*fac
1601           a32=a32*fac
1602           a33=a33*fac
1603 cd          write (iout,'(4i5,4f10.5)')
1604 cd     &     i,itortyp(itype(i)),j,itortyp(itype(j)),a22,a23,a32,a33
1605 cd          write (iout,'(6f10.5)') (muij(k),k=1,4),fac,eel_loc_ij
1606 cd          write (iout,'(2(3f10.5,5x)/2(3f10.5,5x))') uy(:,i),uz(:,i),
1607 cd     &      uy(:,j),uz(:,j)
1608 cd          write (iout,'(4f10.5)')
1609 cd     &      scalar(uy(1,i),uy(1,j)),scalar(uy(1,i),uz(1,j)),
1610 cd     &      scalar(uz(1,i),uy(1,j)),scalar(uz(1,i),uz(1,j))
1611 cd          write (iout,'(4f10.5)') ury,urz,vry,vrz
1612 cd           write (iout,'(9f10.5/)')
1613 cd     &      fac22,a22,fac23,a23,fac32,a32,fac33,a33,eel_loc_ij
1614 C Derivatives of the elements of A in virtual-bond vectors
1615           call unormderiv(erij(1),unmat(1,1),rmij,erder(1,1))
1616           do k=1,3
1617             uryg(k,1)=scalar(erder(1,k),uy(1,i))
1618             uryg(k,2)=scalar(uygrad(1,k,1,i),erij(1))
1619             uryg(k,3)=scalar(uygrad(1,k,2,i),erij(1))
1620             urzg(k,1)=scalar(erder(1,k),uz(1,i))
1621             urzg(k,2)=scalar(uzgrad(1,k,1,i),erij(1))
1622             urzg(k,3)=scalar(uzgrad(1,k,2,i),erij(1))
1623             vryg(k,1)=scalar(erder(1,k),uy(1,j))
1624             vryg(k,2)=scalar(uygrad(1,k,1,j),erij(1))
1625             vryg(k,3)=scalar(uygrad(1,k,2,j),erij(1))
1626             vrzg(k,1)=scalar(erder(1,k),uz(1,j))
1627             vrzg(k,2)=scalar(uzgrad(1,k,1,j),erij(1))
1628             vrzg(k,3)=scalar(uzgrad(1,k,2,j),erij(1))
1629           enddo
1630 C Compute radial contributions to the gradient
1631           facr=-3.0d0*rrmij
1632           a22der=a22*facr
1633           a23der=a23*facr
1634           a32der=a32*facr
1635           a33der=a33*facr
1636           agg(1,1)=a22der*xj
1637           agg(2,1)=a22der*yj
1638           agg(3,1)=a22der*zj
1639           agg(1,2)=a23der*xj
1640           agg(2,2)=a23der*yj
1641           agg(3,2)=a23der*zj
1642           agg(1,3)=a32der*xj
1643           agg(2,3)=a32der*yj
1644           agg(3,3)=a32der*zj
1645           agg(1,4)=a33der*xj
1646           agg(2,4)=a33der*yj
1647           agg(3,4)=a33der*zj
1648 C Add the contributions coming from er
1649           fac3=-3.0d0*fac
1650           do k=1,3
1651             agg(k,1)=agg(k,1)+fac3*(uryg(k,1)*vry+vryg(k,1)*ury)
1652             agg(k,2)=agg(k,2)+fac3*(uryg(k,1)*vrz+vrzg(k,1)*ury)
1653             agg(k,3)=agg(k,3)+fac3*(urzg(k,1)*vry+vryg(k,1)*urz)
1654             agg(k,4)=agg(k,4)+fac3*(urzg(k,1)*vrz+vrzg(k,1)*urz)
1655           enddo
1656           do k=1,3
1657 C Derivatives in DC(i)
1658 cgrad            ghalf1=0.5d0*agg(k,1)
1659 cgrad            ghalf2=0.5d0*agg(k,2)
1660 cgrad            ghalf3=0.5d0*agg(k,3)
1661 cgrad            ghalf4=0.5d0*agg(k,4)
1662             aggi(k,1)=fac*(scalar(uygrad(1,k,1,i),uy(1,j))
1663      &      -3.0d0*uryg(k,2)*vry)!+ghalf1
1664             aggi(k,2)=fac*(scalar(uygrad(1,k,1,i),uz(1,j))
1665      &      -3.0d0*uryg(k,2)*vrz)!+ghalf2
1666             aggi(k,3)=fac*(scalar(uzgrad(1,k,1,i),uy(1,j))
1667      &      -3.0d0*urzg(k,2)*vry)!+ghalf3
1668             aggi(k,4)=fac*(scalar(uzgrad(1,k,1,i),uz(1,j))
1669      &      -3.0d0*urzg(k,2)*vrz)!+ghalf4
1670 C Derivatives in DC(i+1)
1671             aggi1(k,1)=fac*(scalar(uygrad(1,k,2,i),uy(1,j))
1672      &      -3.0d0*uryg(k,3)*vry)!+agg(k,1)
1673             aggi1(k,2)=fac*(scalar(uygrad(1,k,2,i),uz(1,j))
1674      &      -3.0d0*uryg(k,3)*vrz)!+agg(k,2)
1675             aggi1(k,3)=fac*(scalar(uzgrad(1,k,2,i),uy(1,j))
1676      &      -3.0d0*urzg(k,3)*vry)!+agg(k,3)
1677             aggi1(k,4)=fac*(scalar(uzgrad(1,k,2,i),uz(1,j))
1678      &      -3.0d0*urzg(k,3)*vrz)!+agg(k,4)
1679 C Derivatives in DC(j)
1680             aggj(k,1)=fac*(scalar(uygrad(1,k,1,j),uy(1,i))
1681      &      -3.0d0*vryg(k,2)*ury)!+ghalf1
1682             aggj(k,2)=fac*(scalar(uzgrad(1,k,1,j),uy(1,i))
1683      &      -3.0d0*vrzg(k,2)*ury)!+ghalf2
1684             aggj(k,3)=fac*(scalar(uygrad(1,k,1,j),uz(1,i))
1685      &      -3.0d0*vryg(k,2)*urz)!+ghalf3
1686             aggj(k,4)=fac*(scalar(uzgrad(1,k,1,j),uz(1,i))
1687      &      -3.0d0*vrzg(k,2)*urz)!+ghalf4
1688 C Derivatives in DC(j+1) or DC(nres-1)
1689             aggj1(k,1)=fac*(scalar(uygrad(1,k,2,j),uy(1,i))
1690      &      -3.0d0*vryg(k,3)*ury)
1691             aggj1(k,2)=fac*(scalar(uzgrad(1,k,2,j),uy(1,i))
1692      &      -3.0d0*vrzg(k,3)*ury)
1693             aggj1(k,3)=fac*(scalar(uygrad(1,k,2,j),uz(1,i))
1694      &      -3.0d0*vryg(k,3)*urz)
1695             aggj1(k,4)=fac*(scalar(uzgrad(1,k,2,j),uz(1,i))
1696      &      -3.0d0*vrzg(k,3)*urz)
1697 cgrad            if (j.eq.nres-1 .and. i.lt.j-2) then
1698 cgrad              do l=1,4
1699 cgrad                aggj1(k,l)=aggj1(k,l)+agg(k,l)
1700 cgrad              enddo
1701 cgrad            endif
1702           enddo
1703           acipa(1,1)=a22
1704           acipa(1,2)=a23
1705           acipa(2,1)=a32
1706           acipa(2,2)=a33
1707           a22=-a22
1708           a23=-a23
1709           do l=1,2
1710             do k=1,3
1711               agg(k,l)=-agg(k,l)
1712               aggi(k,l)=-aggi(k,l)
1713               aggi1(k,l)=-aggi1(k,l)
1714               aggj(k,l)=-aggj(k,l)
1715               aggj1(k,l)=-aggj1(k,l)
1716             enddo
1717           enddo
1718           if (j.lt.nres-1) then
1719             a22=-a22
1720             a32=-a32
1721             do l=1,3,2
1722               do k=1,3
1723                 agg(k,l)=-agg(k,l)
1724                 aggi(k,l)=-aggi(k,l)
1725                 aggi1(k,l)=-aggi1(k,l)
1726                 aggj(k,l)=-aggj(k,l)
1727                 aggj1(k,l)=-aggj1(k,l)
1728               enddo
1729             enddo
1730           else
1731             a22=-a22
1732             a23=-a23
1733             a32=-a32
1734             a33=-a33
1735             do l=1,4
1736               do k=1,3
1737                 agg(k,l)=-agg(k,l)
1738                 aggi(k,l)=-aggi(k,l)
1739                 aggi1(k,l)=-aggi1(k,l)
1740                 aggj(k,l)=-aggj(k,l)
1741                 aggj1(k,l)=-aggj1(k,l)
1742               enddo
1743             enddo
1744           endif
1745           ENDIF ! WCORR
1746           IF (wel_loc.gt.0.0d0) THEN
1747 C Contribution to the local-electrostatic energy coming from the i-j pair
1748           eel_loc_ij=a22*muij(1)+a23*muij(2)+a32*muij(3)
1749      &     +a33*muij(4)
1750 cd          write (iout,*) 'i',i,' j',j,' eel_loc_ij',eel_loc_ij
1751
1752           if (energy_dec) write (iout,'(a6,2i5,0pf7.3)')
1753      &            'eelloc',i,j,eel_loc_ij
1754
1755           eel_loc=eel_loc+eel_loc_ij
1756 C Partial derivatives in virtual-bond dihedral angles gamma
1757           if (i.gt.1)
1758      &    gel_loc_loc(i-1)=gel_loc_loc(i-1)+
1759      &            a22*muder(1,i)*mu(1,j)+a23*muder(1,i)*mu(2,j)
1760      &           +a32*muder(2,i)*mu(1,j)+a33*muder(2,i)*mu(2,j)
1761           gel_loc_loc(j-1)=gel_loc_loc(j-1)+
1762      &            a22*mu(1,i)*muder(1,j)+a23*mu(1,i)*muder(2,j)
1763      &           +a32*mu(2,i)*muder(1,j)+a33*mu(2,i)*muder(2,j)
1764 C Derivatives of eello in DC(i+1) thru DC(j-1) or DC(nres-2)
1765           do l=1,3
1766             ggg(l)=agg(l,1)*muij(1)+
1767      &          agg(l,2)*muij(2)+agg(l,3)*muij(3)+agg(l,4)*muij(4)
1768             gel_loc_long(l,j)=gel_loc_long(l,j)+ggg(l)
1769             gel_loc_long(l,i)=gel_loc_long(l,i)-ggg(l)
1770 cgrad            ghalf=0.5d0*ggg(l)
1771 cgrad            gel_loc(l,i)=gel_loc(l,i)+ghalf
1772 cgrad            gel_loc(l,j)=gel_loc(l,j)+ghalf
1773           enddo
1774 cgrad          do k=i+1,j2
1775 cgrad            do l=1,3
1776 cgrad              gel_loc(l,k)=gel_loc(l,k)+ggg(l)
1777 cgrad            enddo
1778 cgrad          enddo
1779 C Remaining derivatives of eello
1780           do l=1,3
1781             gel_loc(l,i)=gel_loc(l,i)+aggi(l,1)*muij(1)+
1782      &          aggi(l,2)*muij(2)+aggi(l,3)*muij(3)+aggi(l,4)*muij(4)
1783             gel_loc(l,i+1)=gel_loc(l,i+1)+aggi1(l,1)*muij(1)+
1784      &          aggi1(l,2)*muij(2)+aggi1(l,3)*muij(3)+aggi1(l,4)*muij(4)
1785             gel_loc(l,j)=gel_loc(l,j)+aggj(l,1)*muij(1)+
1786      &          aggj(l,2)*muij(2)+aggj(l,3)*muij(3)+aggj(l,4)*muij(4)
1787             gel_loc(l,j1)=gel_loc(l,j1)+aggj1(l,1)*muij(1)+
1788      &          aggj1(l,2)*muij(2)+aggj1(l,3)*muij(3)+aggj1(l,4)*muij(4)
1789           enddo
1790           ENDIF
1791 C Change 12/26/95 to calculate four-body contributions to H-bonding energy
1792 c          if (j.gt.i+1 .and. num_conti.le.maxconts) then
1793           if (wcorr+wcorr4+wcorr5+wcorr6.gt.0.0d0
1794      &       .and. num_conti.le.maxconts) then
1795 c            write (iout,*) i,j," entered corr"
1796 C
1797 C Calculate the contact function. The ith column of the array JCONT will
1798 C contain the numbers of atoms that make contacts with the atom I (of numbers
1799 C greater than I). The arrays FACONT and GACONT will contain the values of
1800 C the contact function and its derivative.
1801 c           r0ij=1.02D0*rpp(iteli,itelj)
1802 c           r0ij=1.11D0*rpp(iteli,itelj)
1803             r0ij=2.20D0*rpp(iteli,itelj)
1804 c           r0ij=1.55D0*rpp(iteli,itelj)
1805             call gcont(rij,r0ij,1.0D0,0.2d0*r0ij,fcont,fprimcont)
1806             if (fcont.gt.0.0D0) then
1807               num_conti=num_conti+1
1808               if (num_conti.gt.maxconts) then
1809                 write (iout,*) 'WARNING - max. # of contacts exceeded;',
1810      &                         ' will skip next contacts for this conf.'
1811               else
1812                 jcont_hb(num_conti,i)=j
1813 cd                write (iout,*) "i",i," j",j," num_conti",num_conti,
1814 cd     &           " jcont_hb",jcont_hb(num_conti,i)
1815                 IF (wcorr4.gt.0.0d0 .or. wcorr5.gt.0.0d0 .or.
1816      &          wcorr6.gt.0.0d0 .or. wturn6.gt.0.0d0) THEN
1817 C 9/30/99 (AL) - store components necessary to evaluate higher-order loc-el
1818 C  terms.
1819                 d_cont(num_conti,i)=rij
1820 cd                write (2,'(3e15.5)') rij,r0ij+0.2d0*r0ij,rij
1821 C     --- Electrostatic-interaction matrix ---
1822                 a_chuj(1,1,num_conti,i)=a22
1823                 a_chuj(1,2,num_conti,i)=a23
1824                 a_chuj(2,1,num_conti,i)=a32
1825                 a_chuj(2,2,num_conti,i)=a33
1826 C     --- Gradient of rij
1827                 do kkk=1,3
1828                   grij_hb_cont(kkk,num_conti,i)=erij(kkk)
1829                 enddo
1830                 kkll=0
1831                 do k=1,2
1832                   do l=1,2
1833                     kkll=kkll+1
1834                     do m=1,3
1835                       a_chuj_der(k,l,m,1,num_conti,i)=agg(m,kkll)
1836                       a_chuj_der(k,l,m,2,num_conti,i)=aggi(m,kkll)
1837                       a_chuj_der(k,l,m,3,num_conti,i)=aggi1(m,kkll)
1838                       a_chuj_der(k,l,m,4,num_conti,i)=aggj(m,kkll)
1839                       a_chuj_der(k,l,m,5,num_conti,i)=aggj1(m,kkll)
1840                     enddo
1841                   enddo
1842                 enddo
1843                 ENDIF
1844                 IF (wcorr4.eq.0.0d0 .and. wcorr.gt.0.0d0) THEN
1845 C Calculate contact energies
1846                 cosa4=4.0D0*cosa
1847                 wij=cosa-3.0D0*cosb*cosg
1848                 cosbg1=cosb+cosg
1849                 cosbg2=cosb-cosg
1850 c               fac3=dsqrt(-ael6i)/r0ij**3
1851                 fac3=dsqrt(-ael6i)*r3ij
1852 c                 ees0pij=dsqrt(4.0D0+cosa4+wij*wij-3.0D0*cosbg1*cosbg1)
1853                 ees0tmp=4.0D0+cosa4+wij*wij-3.0D0*cosbg1*cosbg1
1854                 if (ees0tmp.gt.0) then
1855                   ees0pij=dsqrt(ees0tmp)
1856                 else
1857                   ees0pij=0
1858                 endif
1859 c                ees0mij=dsqrt(4.0D0-cosa4+wij*wij-3.0D0*cosbg2*cosbg2)
1860                 ees0tmp=4.0D0-cosa4+wij*wij-3.0D0*cosbg2*cosbg2
1861                 if (ees0tmp.gt.0) then
1862                   ees0mij=dsqrt(ees0tmp)
1863                 else
1864                   ees0mij=0
1865                 endif
1866 c               ees0mij=0.0D0
1867                 ees0p(num_conti,i)=0.5D0*fac3*(ees0pij+ees0mij)
1868                 ees0m(num_conti,i)=0.5D0*fac3*(ees0pij-ees0mij)
1869 C Diagnostics. Comment out or remove after debugging!
1870 c               ees0p(num_conti,i)=0.5D0*fac3*ees0pij
1871 c               ees0m(num_conti,i)=0.5D0*fac3*ees0mij
1872 c               ees0m(num_conti,i)=0.0D0
1873 C End diagnostics.
1874 c               write (iout,*) 'i=',i,' j=',j,' rij=',rij,' r0ij=',r0ij,
1875 c    & ' ees0ij=',ees0p(num_conti,i),ees0m(num_conti,i),' fcont=',fcont
1876 C Angular derivatives of the contact function
1877                 ees0pij1=fac3/ees0pij
1878                 ees0mij1=fac3/ees0mij
1879                 fac3p=-3.0D0*fac3*rrmij
1880                 ees0pijp=0.5D0*fac3p*(ees0pij+ees0mij)
1881                 ees0mijp=0.5D0*fac3p*(ees0pij-ees0mij)
1882 c               ees0mij1=0.0D0
1883                 ecosa1=       ees0pij1*( 1.0D0+0.5D0*wij)
1884                 ecosb1=-1.5D0*ees0pij1*(wij*cosg+cosbg1)
1885                 ecosg1=-1.5D0*ees0pij1*(wij*cosb+cosbg1)
1886                 ecosa2=       ees0mij1*(-1.0D0+0.5D0*wij)
1887                 ecosb2=-1.5D0*ees0mij1*(wij*cosg+cosbg2)
1888                 ecosg2=-1.5D0*ees0mij1*(wij*cosb-cosbg2)
1889                 ecosap=ecosa1+ecosa2
1890                 ecosbp=ecosb1+ecosb2
1891                 ecosgp=ecosg1+ecosg2
1892                 ecosam=ecosa1-ecosa2
1893                 ecosbm=ecosb1-ecosb2
1894                 ecosgm=ecosg1-ecosg2
1895 C Diagnostics
1896 c               ecosap=ecosa1
1897 c               ecosbp=ecosb1
1898 c               ecosgp=ecosg1
1899 c               ecosam=0.0D0
1900 c               ecosbm=0.0D0
1901 c               ecosgm=0.0D0
1902 C End diagnostics
1903                 facont_hb(num_conti,i)=fcont
1904                 fprimcont=fprimcont/rij
1905 cd              facont_hb(num_conti,i)=1.0D0
1906 C Following line is for diagnostics.
1907 cd              fprimcont=0.0D0
1908                 do k=1,3
1909                   dcosb(k)=rmij*(dc_norm(k,i)-erij(k)*cosb)
1910                   dcosg(k)=rmij*(dc_norm(k,j)-erij(k)*cosg)
1911                 enddo
1912                 do k=1,3
1913                   gggp(k)=ecosbp*dcosb(k)+ecosgp*dcosg(k)
1914                   gggm(k)=ecosbm*dcosb(k)+ecosgm*dcosg(k)
1915                 enddo
1916                 gggp(1)=gggp(1)+ees0pijp*xj
1917                 gggp(2)=gggp(2)+ees0pijp*yj
1918                 gggp(3)=gggp(3)+ees0pijp*zj
1919                 gggm(1)=gggm(1)+ees0mijp*xj
1920                 gggm(2)=gggm(2)+ees0mijp*yj
1921                 gggm(3)=gggm(3)+ees0mijp*zj
1922 C Derivatives due to the contact function
1923                 gacont_hbr(1,num_conti,i)=fprimcont*xj
1924                 gacont_hbr(2,num_conti,i)=fprimcont*yj
1925                 gacont_hbr(3,num_conti,i)=fprimcont*zj
1926                 do k=1,3
1927 c
1928 c 10/24/08 cgrad and ! comments indicate the parts of the code removed
1929 c          following the change of gradient-summation algorithm.
1930 c
1931 cgrad                  ghalfp=0.5D0*gggp(k)
1932 cgrad                  ghalfm=0.5D0*gggm(k)
1933                   gacontp_hb1(k,num_conti,i)=!ghalfp
1934      &              +(ecosap*(dc_norm(k,j)-cosa*dc_norm(k,i))
1935      &              + ecosbp*(erij(k)-cosb*dc_norm(k,i)))*vbld_inv(i+1)
1936                   gacontp_hb2(k,num_conti,i)=!ghalfp
1937      &              +(ecosap*(dc_norm(k,i)-cosa*dc_norm(k,j))
1938      &              + ecosgp*(erij(k)-cosg*dc_norm(k,j)))*vbld_inv(j+1)
1939                   gacontp_hb3(k,num_conti,i)=gggp(k)
1940                   gacontm_hb1(k,num_conti,i)=!ghalfm
1941      &              +(ecosam*(dc_norm(k,j)-cosa*dc_norm(k,i))
1942      &              + ecosbm*(erij(k)-cosb*dc_norm(k,i)))*vbld_inv(i+1)
1943                   gacontm_hb2(k,num_conti,i)=!ghalfm
1944      &              +(ecosam*(dc_norm(k,i)-cosa*dc_norm(k,j))
1945      &              + ecosgm*(erij(k)-cosg*dc_norm(k,j)))*vbld_inv(j+1)
1946                   gacontm_hb3(k,num_conti,i)=gggm(k)
1947                 enddo
1948               ENDIF ! wcorr
1949               endif  ! num_conti.le.maxconts
1950             endif  ! fcont.gt.0
1951           endif    ! j.gt.i+1
1952           if (wturn3.gt.0.0d0 .or. wturn4.gt.0.0d0) then
1953             do k=1,4
1954               do l=1,3
1955                 ghalf=0.5d0*agg(l,k)
1956                 aggi(l,k)=aggi(l,k)+ghalf
1957                 aggi1(l,k)=aggi1(l,k)+agg(l,k)
1958                 aggj(l,k)=aggj(l,k)+ghalf
1959               enddo
1960             enddo
1961             if (j.eq.nres-1 .and. i.lt.j-2) then
1962               do k=1,4
1963                 do l=1,3
1964                   aggj1(l,k)=aggj1(l,k)+agg(l,k)
1965                 enddo
1966               enddo
1967             endif
1968           endif
1969 c          t_eelecij=t_eelecij+MPI_Wtime()-time00
1970       return
1971       end
1972 C-----------------------------------------------------------------------
1973       subroutine evdwpp_short(evdw1)
1974 C
1975 C Compute Evdwpp
1976 C
1977       implicit real*8 (a-h,o-z)
1978       include 'DIMENSIONS'
1979       include 'COMMON.CONTROL'
1980       include 'COMMON.IOUNITS'
1981       include 'COMMON.GEO'
1982       include 'COMMON.VAR'
1983       include 'COMMON.LOCAL'
1984       include 'COMMON.CHAIN'
1985       include 'COMMON.DERIV'
1986       include 'COMMON.INTERACT'
1987       include 'COMMON.CONTACTS'
1988       include 'COMMON.TORSION'
1989       include 'COMMON.VECTORS'
1990       include 'COMMON.FFIELD'
1991       dimension ggg(3)
1992 c 4/26/02 - AL scaling factor for 1,4 repulsive VDW interactions
1993 #ifdef MOMENT
1994       double precision scal_el /1.0d0/
1995 #else
1996       double precision scal_el /0.5d0/
1997 #endif
1998       evdw1=0.0D0
1999 c      write (iout,*) "iatel_s_vdw",iatel_s_vdw,
2000 c     & " iatel_e_vdw",iatel_e_vdw
2001       call flush(iout)
2002       do i=iatel_s_vdw,iatel_e_vdw
2003         if (itype(i).eq.21 .or. itype(i+1).eq.21) cycle
2004         dxi=dc(1,i)
2005         dyi=dc(2,i)
2006         dzi=dc(3,i)
2007         dx_normi=dc_norm(1,i)
2008         dy_normi=dc_norm(2,i)
2009         dz_normi=dc_norm(3,i)
2010         xmedi=c(1,i)+0.5d0*dxi
2011         ymedi=c(2,i)+0.5d0*dyi
2012         zmedi=c(3,i)+0.5d0*dzi
2013         num_conti=0
2014 c        write (iout,*) 'i',i,' ielstart',ielstart_vdw(i),
2015 c     &   ' ielend',ielend_vdw(i)
2016         call flush(iout)
2017         do j=ielstart_vdw(i),ielend_vdw(i)
2018           if (itype(j).eq.21 .or. itype(j+1).eq.21) cycle
2019           ind=ind+1
2020           iteli=itel(i)
2021           itelj=itel(j)
2022           if (j.eq.i+2 .and. itelj.eq.2) iteli=2
2023           aaa=app(iteli,itelj)
2024           bbb=bpp(iteli,itelj)
2025           dxj=dc(1,j)
2026           dyj=dc(2,j)
2027           dzj=dc(3,j)
2028           dx_normj=dc_norm(1,j)
2029           dy_normj=dc_norm(2,j)
2030           dz_normj=dc_norm(3,j)
2031           xj=c(1,j)+0.5D0*dxj-xmedi
2032           yj=c(2,j)+0.5D0*dyj-ymedi
2033           zj=c(3,j)+0.5D0*dzj-zmedi
2034           rij=xj*xj+yj*yj+zj*zj
2035           rrmij=1.0D0/rij
2036           rij=dsqrt(rij)
2037           sss=sscale(rij/rpp(iteli,itelj))
2038           if (sss.gt.0.0d0) then
2039             rmij=1.0D0/rij
2040             r3ij=rrmij*rmij
2041             r6ij=r3ij*r3ij
2042             ev1=aaa*r6ij*r6ij
2043 c 4/26/02 - AL scaling down 1,4 repulsive VDW interactions
2044             if (j.eq.i+2) ev1=scal_el*ev1
2045             ev2=bbb*r6ij
2046             evdwij=ev1+ev2
2047             if (energy_dec) then
2048               write (iout,'(a6,2i5,0pf7.3,f7.3)') 'evdw1',i,j,evdwij,sss
2049             endif
2050             evdw1=evdw1+evdwij*sss
2051 C
2052 C Calculate contributions to the Cartesian gradient.
2053 C
2054             facvdw=-6*rrmij*(ev1+evdwij)*sss
2055             ggg(1)=facvdw*xj
2056             ggg(2)=facvdw*yj
2057             ggg(3)=facvdw*zj
2058             do k=1,3
2059               gvdwpp(k,j)=gvdwpp(k,j)+ggg(k)
2060               gvdwpp(k,i)=gvdwpp(k,i)-ggg(k)
2061             enddo
2062           endif
2063         enddo ! j
2064       enddo   ! i
2065       return
2066       end
2067 C-----------------------------------------------------------------------------
2068       subroutine escp_long(evdw2,evdw2_14)
2069 C
2070 C This subroutine calculates the excluded-volume interaction energy between
2071 C peptide-group centers and side chains and its gradient in virtual-bond and
2072 C side-chain vectors.
2073 C
2074       implicit real*8 (a-h,o-z)
2075       include 'DIMENSIONS'
2076       include 'COMMON.GEO'
2077       include 'COMMON.VAR'
2078       include 'COMMON.LOCAL'
2079       include 'COMMON.CHAIN'
2080       include 'COMMON.DERIV'
2081       include 'COMMON.INTERACT'
2082       include 'COMMON.FFIELD'
2083       include 'COMMON.IOUNITS'
2084       include 'COMMON.CONTROL'
2085       dimension ggg(3)
2086       evdw2=0.0D0
2087       evdw2_14=0.0d0
2088 cd    print '(a)','Enter ESCP'
2089 cd    write (iout,*) 'iatscp_s=',iatscp_s,' iatscp_e=',iatscp_e
2090       do i=iatscp_s,iatscp_e
2091         if (itype(i).eq.21 .or. itype(i+1).eq.21) cycle
2092         iteli=itel(i)
2093         xi=0.5D0*(c(1,i)+c(1,i+1))
2094         yi=0.5D0*(c(2,i)+c(2,i+1))
2095         zi=0.5D0*(c(3,i)+c(3,i+1))
2096
2097         do iint=1,nscp_gr(i)
2098
2099         do j=iscpstart(i,iint),iscpend(i,iint)
2100           itypj=itype(j)
2101           if (itypj.eq.21) cycle
2102 C Uncomment following three lines for SC-p interactions
2103 c         xj=c(1,nres+j)-xi
2104 c         yj=c(2,nres+j)-yi
2105 c         zj=c(3,nres+j)-zi
2106 C Uncomment following three lines for Ca-p interactions
2107           xj=c(1,j)-xi
2108           yj=c(2,j)-yi
2109           zj=c(3,j)-zi
2110           rrij=1.0D0/(xj*xj+yj*yj+zj*zj)
2111
2112           sss=sscale(1.0d0/(dsqrt(rrij)*rscp(itypj,iteli)))
2113
2114           if (sss.lt.1.0d0) then
2115
2116             fac=rrij**expon2
2117             e1=fac*fac*aad(itypj,iteli)
2118             e2=fac*bad(itypj,iteli)
2119             if (iabs(j-i) .le. 2) then
2120               e1=scal14*e1
2121               e2=scal14*e2
2122               evdw2_14=evdw2_14+(e1+e2)*(1.0d0-sss)
2123             endif
2124             evdwij=e1+e2
2125             evdw2=evdw2+evdwij*(1.0d0-sss)
2126             if (energy_dec) write (iout,'(a6,2i5,2(0pf7.3))')
2127      &          'evdw2',i,j,sss,evdwij
2128 C
2129 C Calculate contributions to the gradient in the virtual-bond and SC vectors.
2130 C
2131             fac=-(evdwij+e1)*rrij*(1.0d0-sss)
2132             ggg(1)=xj*fac
2133             ggg(2)=yj*fac
2134             ggg(3)=zj*fac
2135 C Uncomment following three lines for SC-p interactions
2136 c           do k=1,3
2137 c             gradx_scp(k,j)=gradx_scp(k,j)+ggg(k)
2138 c           enddo
2139 C Uncomment following line for SC-p interactions
2140 c             gradx_scp(k,j)=gradx_scp(k,j)+ggg(k)
2141             do k=1,3
2142               gvdwc_scpp(k,i)=gvdwc_scpp(k,i)-ggg(k)
2143               gvdwc_scp(k,j)=gvdwc_scp(k,j)+ggg(k)
2144             enddo
2145           endif
2146         enddo
2147
2148         enddo ! iint
2149       enddo ! i
2150       do i=1,nct
2151         do j=1,3
2152           gvdwc_scp(j,i)=expon*gvdwc_scp(j,i)
2153           gvdwc_scpp(j,i)=expon*gvdwc_scpp(j,i)
2154           gradx_scp(j,i)=expon*gradx_scp(j,i)
2155         enddo
2156       enddo
2157 C******************************************************************************
2158 C
2159 C                              N O T E !!!
2160 C
2161 C To save time the factor EXPON has been extracted from ALL components
2162 C of GVDWC and GRADX. Remember to multiply them by this factor before further
2163 C use!
2164 C
2165 C******************************************************************************
2166       return
2167       end
2168 C-----------------------------------------------------------------------------
2169       subroutine escp_short(evdw2,evdw2_14)
2170 C
2171 C This subroutine calculates the excluded-volume interaction energy between
2172 C peptide-group centers and side chains and its gradient in virtual-bond and
2173 C side-chain vectors.
2174 C
2175       implicit real*8 (a-h,o-z)
2176       include 'DIMENSIONS'
2177       include 'COMMON.GEO'
2178       include 'COMMON.VAR'
2179       include 'COMMON.LOCAL'
2180       include 'COMMON.CHAIN'
2181       include 'COMMON.DERIV'
2182       include 'COMMON.INTERACT'
2183       include 'COMMON.FFIELD'
2184       include 'COMMON.IOUNITS'
2185       include 'COMMON.CONTROL'
2186       dimension ggg(3)
2187       evdw2=0.0D0
2188       evdw2_14=0.0d0
2189 cd    print '(a)','Enter ESCP'
2190 cd    write (iout,*) 'iatscp_s=',iatscp_s,' iatscp_e=',iatscp_e
2191       do i=iatscp_s,iatscp_e
2192         if (itype(i).eq.21 .or. itype(i+1).eq.21) cycle
2193         iteli=itel(i)
2194         xi=0.5D0*(c(1,i)+c(1,i+1))
2195         yi=0.5D0*(c(2,i)+c(2,i+1))
2196         zi=0.5D0*(c(3,i)+c(3,i+1))
2197
2198         do iint=1,nscp_gr(i)
2199
2200         do j=iscpstart(i,iint),iscpend(i,iint)
2201           itypj=itype(j)
2202           if (itypj.eq.21) cycle
2203 C Uncomment following three lines for SC-p interactions
2204 c         xj=c(1,nres+j)-xi
2205 c         yj=c(2,nres+j)-yi
2206 c         zj=c(3,nres+j)-zi
2207 C Uncomment following three lines for Ca-p interactions
2208           xj=c(1,j)-xi
2209           yj=c(2,j)-yi
2210           zj=c(3,j)-zi
2211           rrij=1.0D0/(xj*xj+yj*yj+zj*zj)
2212
2213           sss=sscale(1.0d0/(dsqrt(rrij)*rscp(itypj,iteli)))
2214
2215           if (sss.gt.0.0d0) then
2216
2217             fac=rrij**expon2
2218             e1=fac*fac*aad(itypj,iteli)
2219             e2=fac*bad(itypj,iteli)
2220             if (iabs(j-i) .le. 2) then
2221               e1=scal14*e1
2222               e2=scal14*e2
2223               evdw2_14=evdw2_14+(e1+e2)*sss
2224             endif
2225             evdwij=e1+e2
2226             evdw2=evdw2+evdwij*sss
2227             if (energy_dec) write (iout,'(a6,2i5,2(0pf7.3))')
2228      &          'evdw2',i,j,sss,evdwij
2229 C
2230 C Calculate contributions to the gradient in the virtual-bond and SC vectors.
2231 C
2232             fac=-(evdwij+e1)*rrij*sss
2233             ggg(1)=xj*fac
2234             ggg(2)=yj*fac
2235             ggg(3)=zj*fac
2236 C Uncomment following three lines for SC-p interactions
2237 c           do k=1,3
2238 c             gradx_scp(k,j)=gradx_scp(k,j)+ggg(k)
2239 c           enddo
2240 C Uncomment following line for SC-p interactions
2241 c             gradx_scp(k,j)=gradx_scp(k,j)+ggg(k)
2242             do k=1,3
2243               gvdwc_scpp(k,i)=gvdwc_scpp(k,i)-ggg(k)
2244               gvdwc_scp(k,j)=gvdwc_scp(k,j)+ggg(k)
2245             enddo
2246           endif
2247         enddo
2248
2249         enddo ! iint
2250       enddo ! i
2251       do i=1,nct
2252         do j=1,3
2253           gvdwc_scp(j,i)=expon*gvdwc_scp(j,i)
2254           gvdwc_scpp(j,i)=expon*gvdwc_scpp(j,i)
2255           gradx_scp(j,i)=expon*gradx_scp(j,i)
2256         enddo
2257       enddo
2258 C******************************************************************************
2259 C
2260 C                              N O T E !!!
2261 C
2262 C To save time the factor EXPON has been extracted from ALL components
2263 C of GVDWC and GRADX. Remember to multiply them by this factor before further
2264 C use!
2265 C
2266 C******************************************************************************
2267       return
2268       end