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

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