source/unres/src_MD-M-newcorr/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.ntyp1) 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.ntyp1) 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.ntyp1) 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.ntyp1) 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.ntyp1) 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.ntyp1) 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.ntyp1) 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.ntyp1) cycle
 307             xj=c(1,nres+j)-xi
 308             yj=c(2,nres+j)-yi
 309             zj=c(3,nres+j)-zi
 310             rrij=1.0D0/(xj*xj+yj*yj+zj*zj)
 311             fac_augm=rrij**expon
 312             e_augm=augm(itypi,itypj)*fac_augm
 313             r_inv_ij=dsqrt(rrij)
 314             rij=1.0D0/r_inv_ij
 315             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.ntyp1) 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.ntyp1) 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.ntyp1) 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.ntyp1) 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.ntyp1) 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.ntyp1) 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.ntyp1) 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.ntyp1) 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.ntyp1) 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.ntyp1) 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.ntyp1) 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.ntyp1) 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.ntyp1.or. itype(i+1).eq.ntyp1
1266      &  .or. itype(i+2).eq.ntyp1 .or. itype(i+3).eq.ntyp1) 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.ntyp1 .or. itype(i+1).eq.ntyp1
1283      &    .or. itype(i+3).eq.ntyp1
1284      &    .or. itype(i+4).eq.ntyp1) 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.ntyp1)
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.ntyp1 .or. itype(i+1).eq.ntyp1) 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.ntyp1 .or. itype(j+1).eq.ntyp1) 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),gmuij1(4),gmuji1(4),
1356      &    gmuij2(4),gmuji2(4)
1357       common /locel/ a_temp,agg,aggi,aggi1,aggj,aggj1,a22,a23,a32,a33,
1358      &    dxi,dyi,dzi,dx_normi,dy_normi,dz_normi,xmedi,ymedi,zmedi,
1359      &    num_conti,j1,j2
1360 c 4/26/02 - AL scaling factor for 1,4 repulsive VDW interactions
1361 #ifdef MOMENT
1362       double precision scal_el /1.0d0/
1363 #else
1364       double precision scal_el /0.5d0/
1365 #endif
1366 C 12/13/98
1367 C 13-go grudnia roku pamietnego...
1368       double precision unmat(3,3) /1.0d0,0.0d0,0.0d0,
1369      &                   0.0d0,1.0d0,0.0d0,
1370      &                   0.0d0,0.0d0,1.0d0/
1371 c          time00=MPI_Wtime()
1372 cd      write (iout,*) "eelecij",i,j
1373           ind=ind+1
1374           iteli=itel(i)
1375           itelj=itel(j)
1376           if (j.eq.i+2 .and. itelj.eq.2) iteli=2
1377           aaa=app(iteli,itelj)
1378           bbb=bpp(iteli,itelj)
1379           ael6i=ael6(iteli,itelj)
1380           ael3i=ael3(iteli,itelj)
1381           dxj=dc(1,j)
1382           dyj=dc(2,j)
1383           dzj=dc(3,j)
1384           dx_normj=dc_norm(1,j)
1385           dy_normj=dc_norm(2,j)
1386           dz_normj=dc_norm(3,j)
1387           xj=c(1,j)+0.5D0*dxj-xmedi
1388           yj=c(2,j)+0.5D0*dyj-ymedi
1389           zj=c(3,j)+0.5D0*dzj-zmedi
1390           rij=xj*xj+yj*yj+zj*zj
1391           rrmij=1.0D0/rij
1392           rij=dsqrt(rij)
1393           rmij=1.0D0/rij
1394 c For extracting the short-range part of Evdwpp
1395           sss=sscale(rij/rpp(iteli,itelj))
1396
1397           r3ij=rrmij*rmij
1398           r6ij=r3ij*r3ij
1399           cosa=dx_normi*dx_normj+dy_normi*dy_normj+dz_normi*dz_normj
1400           cosb=(xj*dx_normi+yj*dy_normi+zj*dz_normi)*rmij
1401           cosg=(xj*dx_normj+yj*dy_normj+zj*dz_normj)*rmij
1402           fac=cosa-3.0D0*cosb*cosg
1403           ev1=aaa*r6ij*r6ij
1404 c 4/26/02 - AL scaling down 1,4 repulsive VDW interactions
1405           if (j.eq.i+2) ev1=scal_el*ev1
1406           ev2=bbb*r6ij
1407           fac3=ael6i*r6ij
1408           fac4=ael3i*r3ij
1409           evdwij=ev1+ev2
1410           el1=fac3*(4.0D0+fac*fac-3.0D0*(cosb*cosb+cosg*cosg))
1411           el2=fac4*fac
1412           eesij=el1+el2
1413 C 12/26/95 - for the evaluation of multi-body H-bonding interactions
1414           ees0ij=4.0D0+fac*fac-3.0D0*(cosb*cosb+cosg*cosg)
1415           ees=ees+eesij
1416           evdw1=evdw1+evdwij*(1.0d0-sss)
1417 cd          write(iout,'(2(2i3,2x),7(1pd12.4)/2(3(1pd12.4),5x)/)')
1418 cd     &      iteli,i,itelj,j,aaa,bbb,ael6i,ael3i,
1419 cd     &      1.0D0/dsqrt(rrmij),evdwij,eesij,
1420 cd     &      xmedi,ymedi,zmedi,xj,yj,zj
1421
1422           if (energy_dec) then
1423               write (iout,'(a6,2i5,0pf7.3,f7.3)') 'evdw1',i,j,evdwij,sss
1424               write (iout,'(a6,2i5,0pf7.3)') 'ees',i,j,eesij
1425           endif
1426
1427 C
1428 C Calculate contributions to the Cartesian gradient.
1429 C
1430 #ifdef SPLITELE
1431           facvdw=-6*rrmij*(ev1+evdwij)*(1.0d0-sss)
1432           facel=-3*rrmij*(el1+eesij)
1433           fac1=fac
1434           erij(1)=xj*rmij
1435           erij(2)=yj*rmij
1436           erij(3)=zj*rmij
1437 *
1438 * Radial derivatives. First process both termini of the fragment (i,j)
1439 *
1440           ggg(1)=facel*xj
1441           ggg(2)=facel*yj
1442           ggg(3)=facel*zj
1443 c          do k=1,3
1444 c            ghalf=0.5D0*ggg(k)
1445 c            gelc(k,i)=gelc(k,i)+ghalf
1446 c            gelc(k,j)=gelc(k,j)+ghalf
1447 c          enddo
1448 c 9/28/08 AL Gradient compotents will be summed only at the end
1449           do k=1,3
1450             gelc_long(k,j)=gelc_long(k,j)+ggg(k)
1451             gelc_long(k,i)=gelc_long(k,i)-ggg(k)
1452           enddo
1453 *
1454 * Loop over residues i+1 thru j-1.
1455 *
1456 cgrad          do k=i+1,j-1
1457 cgrad            do l=1,3
1458 cgrad              gelc(l,k)=gelc(l,k)+ggg(l)
1459 cgrad            enddo
1460 cgrad          enddo
1461           ggg(1)=facvdw*xj
1462           ggg(2)=facvdw*yj
1463           ggg(3)=facvdw*zj
1464 c          do k=1,3
1465 c            ghalf=0.5D0*ggg(k)
1466 c            gvdwpp(k,i)=gvdwpp(k,i)+ghalf
1467 c            gvdwpp(k,j)=gvdwpp(k,j)+ghalf
1468 c          enddo
1469 c 9/28/08 AL Gradient compotents will be summed only at the end
1470           do k=1,3
1471             gvdwpp(k,j)=gvdwpp(k,j)+ggg(k)
1472             gvdwpp(k,i)=gvdwpp(k,i)-ggg(k)
1473           enddo
1474 *
1475 * Loop over residues i+1 thru j-1.
1476 *
1477 cgrad          do k=i+1,j-1
1478 cgrad            do l=1,3
1479 cgrad              gvdwpp(l,k)=gvdwpp(l,k)+ggg(l)
1480 cgrad            enddo
1481 cgrad          enddo
1482 #else
1483           facvdw=ev1+evdwij*(1.0d0-sss)
1484           facel=el1+eesij
1485           fac1=fac
1486           fac=-3*rrmij*(facvdw+facvdw+facel)
1487           erij(1)=xj*rmij
1488           erij(2)=yj*rmij
1489           erij(3)=zj*rmij
1490 *
1491 * Radial derivatives. First process both termini of the fragment (i,j)
1492 *
1493           ggg(1)=fac*xj
1494           ggg(2)=fac*yj
1495           ggg(3)=fac*zj
1496 c          do k=1,3
1497 c            ghalf=0.5D0*ggg(k)
1498 c            gelc(k,i)=gelc(k,i)+ghalf
1499 c            gelc(k,j)=gelc(k,j)+ghalf
1500 c          enddo
1501 c 9/28/08 AL Gradient compotents will be summed only at the end
1502           do k=1,3
1503             gelc_long(k,j)=gelc(k,j)+ggg(k)
1504             gelc_long(k,i)=gelc(k,i)-ggg(k)
1505           enddo
1506 *
1507 * Loop over residues i+1 thru j-1.
1508 *
1509 cgrad          do k=i+1,j-1
1510 cgrad            do l=1,3
1511 cgrad              gelc(l,k)=gelc(l,k)+ggg(l)
1512 cgrad            enddo
1513 cgrad          enddo
1514 c 9/28/08 AL Gradient compotents will be summed only at the end
1515           ggg(1)=facvdw*xj
1516           ggg(2)=facvdw*yj
1517           ggg(3)=facvdw*zj
1518           do k=1,3
1519             gvdwpp(k,j)=gvdwpp(k,j)+ggg(k)
1520             gvdwpp(k,i)=gvdwpp(k,i)-ggg(k)
1521           enddo
1522 #endif
1523 *
1524 * Angular part
1525 *
1526           ecosa=2.0D0*fac3*fac1+fac4
1527           fac4=-3.0D0*fac4
1528           fac3=-6.0D0*fac3
1529           ecosb=(fac3*(fac1*cosg+cosb)+cosg*fac4)
1530           ecosg=(fac3*(fac1*cosb+cosg)+cosb*fac4)
1531           do k=1,3
1532             dcosb(k)=rmij*(dc_norm(k,i)-erij(k)*cosb)
1533             dcosg(k)=rmij*(dc_norm(k,j)-erij(k)*cosg)
1534           enddo
1535 cd        print '(2i3,2(3(1pd14.5),3x))',i,j,(dcosb(k),k=1,3),
1536 cd   &          (dcosg(k),k=1,3)
1537           do k=1,3
1538             ggg(k)=ecosb*dcosb(k)+ecosg*dcosg(k)
1539           enddo
1540 c          do k=1,3
1541 c            ghalf=0.5D0*ggg(k)
1542 c            gelc(k,i)=gelc(k,i)+ghalf
1543 c     &               +(ecosa*(dc_norm(k,j)-cosa*dc_norm(k,i))
1544 c     &               + ecosb*(erij(k)-cosb*dc_norm(k,i)))*vbld_inv(i+1)
1545 c            gelc(k,j)=gelc(k,j)+ghalf
1546 c     &               +(ecosa*(dc_norm(k,i)-cosa*dc_norm(k,j))
1547 c     &               + ecosg*(erij(k)-cosg*dc_norm(k,j)))*vbld_inv(j+1)
1548 c          enddo
1549 cgrad          do k=i+1,j-1
1550 cgrad            do l=1,3
1551 cgrad              gelc(l,k)=gelc(l,k)+ggg(l)
1552 cgrad            enddo
1553 cgrad          enddo
1554           do k=1,3
1555             gelc(k,i)=gelc(k,i)
1556      &               +(ecosa*(dc_norm(k,j)-cosa*dc_norm(k,i))
1557      &               + ecosb*(erij(k)-cosb*dc_norm(k,i)))*vbld_inv(i+1)
1558             gelc(k,j)=gelc(k,j)
1559      &               +(ecosa*(dc_norm(k,i)-cosa*dc_norm(k,j))
1560      &               + ecosg*(erij(k)-cosg*dc_norm(k,j)))*vbld_inv(j+1)
1561             gelc_long(k,j)=gelc_long(k,j)+ggg(k)
1562             gelc_long(k,i)=gelc_long(k,i)-ggg(k)
1563           enddo
1564           IF (wel_loc.gt.0.0d0 .or. wcorr4.gt.0.0d0 .or. wcorr5.gt.0.0d0
1565      &        .or. wcorr6.gt.0.0d0 .or. wturn3.gt.0.0d0
1566      &        .or. wturn4.gt.0.0d0 .or. wturn6.gt.0.0d0) THEN
1567 C
1568 C 9/25/99 Mixed third-order local-electrostatic terms. The local-interaction
1569 C   energy of a peptide unit is assumed in the form of a second-order
1570 C   Fourier series in the angles lambda1 and lambda2 (see Nishikawa et al.
1571 C   Macromolecules, 1974, 7, 797-806 for definition). This correlation terms
1572 C   are computed for EVERY pair of non-contiguous peptide groups.
1573 C
1574           if (j.lt.nres-1) then
1575             j1=j+1
1576             j2=j-1
1577           else
1578             j1=j-1
1579             j2=j-2
1580           endif
1581           kkk=0
1582           do k=1,2
1583             do l=1,2
1584               kkk=kkk+1
1585               muij(kkk)=mu(k,i)*mu(l,j)
1586 #ifdef NEWCORR
1587              gmuij1(kkk)=gtb1(k,i+1)*mu(l,j)
1588 c             write(iout,*) 'kkk=', gtb1(k,i)*mu(l,j),gtb1(k,i),k,i
1589              gmuij2(kkk)=gUb2(k,i)*mu(l,j)
1590              gmuji1(kkk)=mu(k,i)*gtb1(l,j+1)
1591 c             write(iout,*) 'kkk=', gtb1(k,i)*mu(l,j),gtb1(l,j),l,j
1592              gmuji2(kkk)=mu(k,i)*gUb2(l,j)
1593 #endif
1594             enddo
1595           enddo
1596 cd         write (iout,*) 'EELEC: i',i,' j',j
1597 cd          write (iout,*) 'j',j,' j1',j1,' j2',j2
1598 cd          write(iout,*) 'muij',muij
1599           ury=scalar(uy(1,i),erij)
1600           urz=scalar(uz(1,i),erij)
1601           vry=scalar(uy(1,j),erij)
1602           vrz=scalar(uz(1,j),erij)
1603           a22=scalar(uy(1,i),uy(1,j))-3*ury*vry
1604           a23=scalar(uy(1,i),uz(1,j))-3*ury*vrz
1605           a32=scalar(uz(1,i),uy(1,j))-3*urz*vry
1606           a33=scalar(uz(1,i),uz(1,j))-3*urz*vrz
1607           fac=dsqrt(-ael6i)*r3ij
1608           a22=a22*fac
1609           a23=a23*fac
1610           a32=a32*fac
1611           a33=a33*fac
1612 cd          write (iout,'(4i5,4f10.5)')
1613 cd     &     i,itortyp(itype(i)),j,itortyp(itype(j)),a22,a23,a32,a33
1614 cd          write (iout,'(6f10.5)') (muij(k),k=1,4),fac,eel_loc_ij
1615 cd          write (iout,'(2(3f10.5,5x)/2(3f10.5,5x))') uy(:,i),uz(:,i),
1616 cd     &      uy(:,j),uz(:,j)
1617 cd          write (iout,'(4f10.5)')
1618 cd     &      scalar(uy(1,i),uy(1,j)),scalar(uy(1,i),uz(1,j)),
1619 cd     &      scalar(uz(1,i),uy(1,j)),scalar(uz(1,i),uz(1,j))
1620 cd          write (iout,'(4f10.5)') ury,urz,vry,vrz
1621 cd           write (iout,'(9f10.5/)')
1622 cd     &      fac22,a22,fac23,a23,fac32,a32,fac33,a33,eel_loc_ij
1623 C Derivatives of the elements of A in virtual-bond vectors
1624           call unormderiv(erij(1),unmat(1,1),rmij,erder(1,1))
1625           do k=1,3
1626             uryg(k,1)=scalar(erder(1,k),uy(1,i))
1627             uryg(k,2)=scalar(uygrad(1,k,1,i),erij(1))
1628             uryg(k,3)=scalar(uygrad(1,k,2,i),erij(1))
1629             urzg(k,1)=scalar(erder(1,k),uz(1,i))
1630             urzg(k,2)=scalar(uzgrad(1,k,1,i),erij(1))
1631             urzg(k,3)=scalar(uzgrad(1,k,2,i),erij(1))
1632             vryg(k,1)=scalar(erder(1,k),uy(1,j))
1633             vryg(k,2)=scalar(uygrad(1,k,1,j),erij(1))
1634             vryg(k,3)=scalar(uygrad(1,k,2,j),erij(1))
1635             vrzg(k,1)=scalar(erder(1,k),uz(1,j))
1636             vrzg(k,2)=scalar(uzgrad(1,k,1,j),erij(1))
1637             vrzg(k,3)=scalar(uzgrad(1,k,2,j),erij(1))
1638           enddo
1639 C Compute radial contributions to the gradient
1640           facr=-3.0d0*rrmij
1641           a22der=a22*facr
1642           a23der=a23*facr
1643           a32der=a32*facr
1644           a33der=a33*facr
1645           agg(1,1)=a22der*xj
1646           agg(2,1)=a22der*yj
1647           agg(3,1)=a22der*zj
1648           agg(1,2)=a23der*xj
1649           agg(2,2)=a23der*yj
1650           agg(3,2)=a23der*zj
1651           agg(1,3)=a32der*xj
1652           agg(2,3)=a32der*yj
1653           agg(3,3)=a32der*zj
1654           agg(1,4)=a33der*xj
1655           agg(2,4)=a33der*yj
1656           agg(3,4)=a33der*zj
1657 C Add the contributions coming from er
1658           fac3=-3.0d0*fac
1659           do k=1,3
1660             agg(k,1)=agg(k,1)+fac3*(uryg(k,1)*vry+vryg(k,1)*ury)
1661             agg(k,2)=agg(k,2)+fac3*(uryg(k,1)*vrz+vrzg(k,1)*ury)
1662             agg(k,3)=agg(k,3)+fac3*(urzg(k,1)*vry+vryg(k,1)*urz)
1663             agg(k,4)=agg(k,4)+fac3*(urzg(k,1)*vrz+vrzg(k,1)*urz)
1664           enddo
1665           do k=1,3
1666 C Derivatives in DC(i)
1667 cgrad            ghalf1=0.5d0*agg(k,1)
1668 cgrad            ghalf2=0.5d0*agg(k,2)
1669 cgrad            ghalf3=0.5d0*agg(k,3)
1670 cgrad            ghalf4=0.5d0*agg(k,4)
1671             aggi(k,1)=fac*(scalar(uygrad(1,k,1,i),uy(1,j))
1672      &      -3.0d0*uryg(k,2)*vry)!+ghalf1
1673             aggi(k,2)=fac*(scalar(uygrad(1,k,1,i),uz(1,j))
1674      &      -3.0d0*uryg(k,2)*vrz)!+ghalf2
1675             aggi(k,3)=fac*(scalar(uzgrad(1,k,1,i),uy(1,j))
1676      &      -3.0d0*urzg(k,2)*vry)!+ghalf3
1677             aggi(k,4)=fac*(scalar(uzgrad(1,k,1,i),uz(1,j))
1678      &      -3.0d0*urzg(k,2)*vrz)!+ghalf4
1679 C Derivatives in DC(i+1)
1680             aggi1(k,1)=fac*(scalar(uygrad(1,k,2,i),uy(1,j))
1681      &      -3.0d0*uryg(k,3)*vry)!+agg(k,1)
1682             aggi1(k,2)=fac*(scalar(uygrad(1,k,2,i),uz(1,j))
1683      &      -3.0d0*uryg(k,3)*vrz)!+agg(k,2)
1684             aggi1(k,3)=fac*(scalar(uzgrad(1,k,2,i),uy(1,j))
1685      &      -3.0d0*urzg(k,3)*vry)!+agg(k,3)
1686             aggi1(k,4)=fac*(scalar(uzgrad(1,k,2,i),uz(1,j))
1687      &      -3.0d0*urzg(k,3)*vrz)!+agg(k,4)
1688 C Derivatives in DC(j)
1689             aggj(k,1)=fac*(scalar(uygrad(1,k,1,j),uy(1,i))
1690      &      -3.0d0*vryg(k,2)*ury)!+ghalf1
1691             aggj(k,2)=fac*(scalar(uzgrad(1,k,1,j),uy(1,i))
1692      &      -3.0d0*vrzg(k,2)*ury)!+ghalf2
1693             aggj(k,3)=fac*(scalar(uygrad(1,k,1,j),uz(1,i))
1694      &      -3.0d0*vryg(k,2)*urz)!+ghalf3
1695             aggj(k,4)=fac*(scalar(uzgrad(1,k,1,j),uz(1,i))
1696      &      -3.0d0*vrzg(k,2)*urz)!+ghalf4
1697 C Derivatives in DC(j+1) or DC(nres-1)
1698             aggj1(k,1)=fac*(scalar(uygrad(1,k,2,j),uy(1,i))
1699      &      -3.0d0*vryg(k,3)*ury)
1700             aggj1(k,2)=fac*(scalar(uzgrad(1,k,2,j),uy(1,i))
1701      &      -3.0d0*vrzg(k,3)*ury)
1702             aggj1(k,3)=fac*(scalar(uygrad(1,k,2,j),uz(1,i))
1703      &      -3.0d0*vryg(k,3)*urz)
1704             aggj1(k,4)=fac*(scalar(uzgrad(1,k,2,j),uz(1,i))
1705      &      -3.0d0*vrzg(k,3)*urz)
1706 cgrad            if (j.eq.nres-1 .and. i.lt.j-2) then
1707 cgrad              do l=1,4
1708 cgrad                aggj1(k,l)=aggj1(k,l)+agg(k,l)
1709 cgrad              enddo
1710 cgrad            endif
1711           enddo
1712           acipa(1,1)=a22
1713           acipa(1,2)=a23
1714           acipa(2,1)=a32
1715           acipa(2,2)=a33
1716           a22=-a22
1717           a23=-a23
1718           do l=1,2
1719             do k=1,3
1720               agg(k,l)=-agg(k,l)
1721               aggi(k,l)=-aggi(k,l)
1722               aggi1(k,l)=-aggi1(k,l)
1723               aggj(k,l)=-aggj(k,l)
1724               aggj1(k,l)=-aggj1(k,l)
1725             enddo
1726           enddo
1727           if (j.lt.nres-1) then
1728             a22=-a22
1729             a32=-a32
1730             do l=1,3,2
1731               do k=1,3
1732                 agg(k,l)=-agg(k,l)
1733                 aggi(k,l)=-aggi(k,l)
1734                 aggi1(k,l)=-aggi1(k,l)
1735                 aggj(k,l)=-aggj(k,l)
1736                 aggj1(k,l)=-aggj1(k,l)
1737               enddo
1738             enddo
1739           else
1740             a22=-a22
1741             a23=-a23
1742             a32=-a32
1743             a33=-a33
1744             do l=1,4
1745               do k=1,3
1746                 agg(k,l)=-agg(k,l)
1747                 aggi(k,l)=-aggi(k,l)
1748                 aggi1(k,l)=-aggi1(k,l)
1749                 aggj(k,l)=-aggj(k,l)
1750                 aggj1(k,l)=-aggj1(k,l)
1751               enddo
1752             enddo
1753           endif
1754           ENDIF ! WCORR
1755           IF (wel_loc.gt.0.0d0) THEN
1756 C Contribution to the local-electrostatic energy coming from the i-j pair
1757           eel_loc_ij=a22*muij(1)+a23*muij(2)+a32*muij(3)
1758      &     +a33*muij(4)
1759 cd          write (iout,*) 'i',i,' j',j,' eel_loc_ij',eel_loc_ij
1760 C Calculate patrial derivative for theta angle
1761 #ifdef NEWCORR
1762          geel_loc_ij=a22*gmuij1(1)
1763      &     +a23*gmuij1(2)
1764      &     +a32*gmuij1(3)
1765      &     +a33*gmuij1(4)
1766 c         write(iout,*) "derivative over thatai"
1767 c         write(iout,*) a22*gmuij1(1), a23*gmuij1(2) ,a32*gmuij1(3),
1768 c     &   a33*gmuij1(4)
1769          gloc(nphi+i,icg)=gloc(nphi+i,icg)+
1770      &      geel_loc_ij*wel_loc
1771 c         write(iout,*) "derivative over thatai-1"
1772 c         write(iout,*) a22*gmuij2(1), a23*gmuij2(2) ,a32*gmuij2(3),
1773 c     &   a33*gmuij2(4)
1774          geel_loc_ij=a22*gmuij2(1)+a23*gmuij2(2)+a32*gmuij2(3)
1775      &     +a33*gmuij2(4)
1776          gloc(nphi+i-1,icg)=gloc(nphi+i-1,icg)+
1777      &      geel_loc_ij*wel_loc
1778          geel_loc_ji=a22*gmuji1(1)+a23*gmuji1(2)+a32*gmuji1(3)
1779      &     +a33*gmuji1(4)
1780 c         write(iout,*) "derivative over thataj"
1781 c         write(iout,*) a22*gmuji1(1), a23*gmuji1(2) ,a32*gmuji1(3),
1782 c     &   a33*gmuji1(4)
1783
1784          gloc(nphi+j,icg)=gloc(nphi+j,icg)+
1785      &      geel_loc_ji*wel_loc
1786          geel_loc_ji=a22*gmuji2(1)+a23*gmuji2(2)+a32*gmuji2(3)
1787      &     +a33*gmuji2(4)
1788 c         write(iout,*) "derivative over thataj-1"
1789 c         write(iout,*) a22*gmuji2(1), a23*gmuji2(2) ,a32*gmuji2(3),
1790 c     &   a33*gmuji2(4)
1791          gloc(nphi+j-1,icg)=gloc(nphi+j-1,icg)+
1792      &      geel_loc_ji*wel_loc
1793 #endif
1794           if (energy_dec) write (iout,'(a6,2i5,0pf7.3)')
1795      &            'eelloc',i,j,eel_loc_ij
1796
1797           eel_loc=eel_loc+eel_loc_ij
1798 C Partial derivatives in virtual-bond dihedral angles gamma
1799           if (i.gt.1)
1800      &    gel_loc_loc(i-1)=gel_loc_loc(i-1)+
1801      &            a22*muder(1,i)*mu(1,j)+a23*muder(1,i)*mu(2,j)
1802      &           +a32*muder(2,i)*mu(1,j)+a33*muder(2,i)*mu(2,j)
1803           gel_loc_loc(j-1)=gel_loc_loc(j-1)+
1804      &            a22*mu(1,i)*muder(1,j)+a23*mu(1,i)*muder(2,j)
1805      &           +a32*mu(2,i)*muder(1,j)+a33*mu(2,i)*muder(2,j)
1806 C Derivatives of eello in DC(i+1) thru DC(j-1) or DC(nres-2)
1807           do l=1,3
1808             ggg(l)=agg(l,1)*muij(1)+
1809      &          agg(l,2)*muij(2)+agg(l,3)*muij(3)+agg(l,4)*muij(4)
1810             gel_loc_long(l,j)=gel_loc_long(l,j)+ggg(l)
1811             gel_loc_long(l,i)=gel_loc_long(l,i)-ggg(l)
1812 cgrad            ghalf=0.5d0*ggg(l)
1813 cgrad            gel_loc(l,i)=gel_loc(l,i)+ghalf
1814 cgrad            gel_loc(l,j)=gel_loc(l,j)+ghalf
1815           enddo
1816 cgrad          do k=i+1,j2
1817 cgrad            do l=1,3
1818 cgrad              gel_loc(l,k)=gel_loc(l,k)+ggg(l)
1819 cgrad            enddo
1820 cgrad          enddo
1821 C Remaining derivatives of eello
1822           do l=1,3
1823             gel_loc(l,i)=gel_loc(l,i)+aggi(l,1)*muij(1)+
1824      &          aggi(l,2)*muij(2)+aggi(l,3)*muij(3)+aggi(l,4)*muij(4)
1825             gel_loc(l,i+1)=gel_loc(l,i+1)+aggi1(l,1)*muij(1)+
1826      &          aggi1(l,2)*muij(2)+aggi1(l,3)*muij(3)+aggi1(l,4)*muij(4)
1827             gel_loc(l,j)=gel_loc(l,j)+aggj(l,1)*muij(1)+
1828      &          aggj(l,2)*muij(2)+aggj(l,3)*muij(3)+aggj(l,4)*muij(4)
1829             gel_loc(l,j1)=gel_loc(l,j1)+aggj1(l,1)*muij(1)+
1830      &          aggj1(l,2)*muij(2)+aggj1(l,3)*muij(3)+aggj1(l,4)*muij(4)
1831           enddo
1832           ENDIF
1833 C Change 12/26/95 to calculate four-body contributions to H-bonding energy
1834 c          if (j.gt.i+1 .and. num_conti.le.maxconts) then
1835           if (wcorr+wcorr4+wcorr5+wcorr6.gt.0.0d0
1836      &       .and. num_conti.le.maxconts) then
1837 c            write (iout,*) i,j," entered corr"
1838 C
1839 C Calculate the contact function. The ith column of the array JCONT will
1840 C contain the numbers of atoms that make contacts with the atom I (of numbers
1841 C greater than I). The arrays FACONT and GACONT will contain the values of
1842 C the contact function and its derivative.
1843 c           r0ij=1.02D0*rpp(iteli,itelj)
1844 c           r0ij=1.11D0*rpp(iteli,itelj)
1845             r0ij=2.20D0*rpp(iteli,itelj)
1846 c           r0ij=1.55D0*rpp(iteli,itelj)
1847             call gcont(rij,r0ij,1.0D0,0.2d0*r0ij,fcont,fprimcont)
1848             if (fcont.gt.0.0D0) then
1849               num_conti=num_conti+1
1850               if (num_conti.gt.maxconts) then
1851                 write (iout,*) 'WARNING - max. # of contacts exceeded;',
1852      &                         ' will skip next contacts for this conf.'
1853               else
1854                 jcont_hb(num_conti,i)=j
1855 cd                write (iout,*) "i",i," j",j," num_conti",num_conti,
1856 cd     &           " jcont_hb",jcont_hb(num_conti,i)
1857                 IF (wcorr4.gt.0.0d0 .or. wcorr5.gt.0.0d0 .or.
1858      &          wcorr6.gt.0.0d0 .or. wturn6.gt.0.0d0) THEN
1859 C 9/30/99 (AL) - store components necessary to evaluate higher-order loc-el
1860 C  terms.
1861                 d_cont(num_conti,i)=rij
1862 cd                write (2,'(3e15.5)') rij,r0ij+0.2d0*r0ij,rij
1863 C     --- Electrostatic-interaction matrix ---
1864                 a_chuj(1,1,num_conti,i)=a22
1865                 a_chuj(1,2,num_conti,i)=a23
1866                 a_chuj(2,1,num_conti,i)=a32
1867                 a_chuj(2,2,num_conti,i)=a33
1868 C     --- Gradient of rij
1869                 do kkk=1,3
1870                   grij_hb_cont(kkk,num_conti,i)=erij(kkk)
1871                 enddo
1872                 kkll=0
1873                 do k=1,2
1874                   do l=1,2
1875                     kkll=kkll+1
1876                     do m=1,3
1877                       a_chuj_der(k,l,m,1,num_conti,i)=agg(m,kkll)
1878                       a_chuj_der(k,l,m,2,num_conti,i)=aggi(m,kkll)
1879                       a_chuj_der(k,l,m,3,num_conti,i)=aggi1(m,kkll)
1880                       a_chuj_der(k,l,m,4,num_conti,i)=aggj(m,kkll)
1881                       a_chuj_der(k,l,m,5,num_conti,i)=aggj1(m,kkll)
1882                     enddo
1883                   enddo
1884                 enddo
1885                 ENDIF
1886                 IF (wcorr4.eq.0.0d0 .and. wcorr.gt.0.0d0) THEN
1887 C Calculate contact energies
1888                 cosa4=4.0D0*cosa
1889                 wij=cosa-3.0D0*cosb*cosg
1890                 cosbg1=cosb+cosg
1891                 cosbg2=cosb-cosg
1892 c               fac3=dsqrt(-ael6i)/r0ij**3
1893                 fac3=dsqrt(-ael6i)*r3ij
1894 c                 ees0pij=dsqrt(4.0D0+cosa4+wij*wij-3.0D0*cosbg1*cosbg1)
1895                 ees0tmp=4.0D0+cosa4+wij*wij-3.0D0*cosbg1*cosbg1
1896                 if (ees0tmp.gt.0) then
1897                   ees0pij=dsqrt(ees0tmp)
1898                 else
1899                   ees0pij=0
1900                 endif
1901 c                ees0mij=dsqrt(4.0D0-cosa4+wij*wij-3.0D0*cosbg2*cosbg2)
1902                 ees0tmp=4.0D0-cosa4+wij*wij-3.0D0*cosbg2*cosbg2
1903                 if (ees0tmp.gt.0) then
1904                   ees0mij=dsqrt(ees0tmp)
1905                 else
1906                   ees0mij=0
1907                 endif
1908 c               ees0mij=0.0D0
1909                 ees0p(num_conti,i)=0.5D0*fac3*(ees0pij+ees0mij)
1910                 ees0m(num_conti,i)=0.5D0*fac3*(ees0pij-ees0mij)
1911 C Diagnostics. Comment out or remove after debugging!
1912 c               ees0p(num_conti,i)=0.5D0*fac3*ees0pij
1913 c               ees0m(num_conti,i)=0.5D0*fac3*ees0mij
1914 c               ees0m(num_conti,i)=0.0D0
1915 C End diagnostics.
1916 c               write (iout,*) 'i=',i,' j=',j,' rij=',rij,' r0ij=',r0ij,
1917 c    & ' ees0ij=',ees0p(num_conti,i),ees0m(num_conti,i),' fcont=',fcont
1918 C Angular derivatives of the contact function
1919                 ees0pij1=fac3/ees0pij
1920                 ees0mij1=fac3/ees0mij
1921                 fac3p=-3.0D0*fac3*rrmij
1922                 ees0pijp=0.5D0*fac3p*(ees0pij+ees0mij)
1923                 ees0mijp=0.5D0*fac3p*(ees0pij-ees0mij)
1924 c               ees0mij1=0.0D0
1925                 ecosa1=       ees0pij1*( 1.0D0+0.5D0*wij)
1926                 ecosb1=-1.5D0*ees0pij1*(wij*cosg+cosbg1)
1927                 ecosg1=-1.5D0*ees0pij1*(wij*cosb+cosbg1)
1928                 ecosa2=       ees0mij1*(-1.0D0+0.5D0*wij)
1929                 ecosb2=-1.5D0*ees0mij1*(wij*cosg+cosbg2)
1930                 ecosg2=-1.5D0*ees0mij1*(wij*cosb-cosbg2)
1931                 ecosap=ecosa1+ecosa2
1932                 ecosbp=ecosb1+ecosb2
1933                 ecosgp=ecosg1+ecosg2
1934                 ecosam=ecosa1-ecosa2
1935                 ecosbm=ecosb1-ecosb2
1936                 ecosgm=ecosg1-ecosg2
1937 C Diagnostics
1938 c               ecosap=ecosa1
1939 c               ecosbp=ecosb1
1940 c               ecosgp=ecosg1
1941 c               ecosam=0.0D0
1942 c               ecosbm=0.0D0
1943 c               ecosgm=0.0D0
1944 C End diagnostics
1945                 facont_hb(num_conti,i)=fcont
1946                 fprimcont=fprimcont/rij
1947 cd              facont_hb(num_conti,i)=1.0D0
1948 C Following line is for diagnostics.
1949 cd              fprimcont=0.0D0
1950                 do k=1,3
1951                   dcosb(k)=rmij*(dc_norm(k,i)-erij(k)*cosb)
1952                   dcosg(k)=rmij*(dc_norm(k,j)-erij(k)*cosg)
1953                 enddo
1954                 do k=1,3
1955                   gggp(k)=ecosbp*dcosb(k)+ecosgp*dcosg(k)
1956                   gggm(k)=ecosbm*dcosb(k)+ecosgm*dcosg(k)
1957                 enddo
1958                 gggp(1)=gggp(1)+ees0pijp*xj
1959                 gggp(2)=gggp(2)+ees0pijp*yj
1960                 gggp(3)=gggp(3)+ees0pijp*zj
1961                 gggm(1)=gggm(1)+ees0mijp*xj
1962                 gggm(2)=gggm(2)+ees0mijp*yj
1963                 gggm(3)=gggm(3)+ees0mijp*zj
1964 C Derivatives due to the contact function
1965                 gacont_hbr(1,num_conti,i)=fprimcont*xj
1966                 gacont_hbr(2,num_conti,i)=fprimcont*yj
1967                 gacont_hbr(3,num_conti,i)=fprimcont*zj
1968                 do k=1,3
1969 c
1970 c 10/24/08 cgrad and ! comments indicate the parts of the code removed
1971 c          following the change of gradient-summation algorithm.
1972 c
1973 cgrad                  ghalfp=0.5D0*gggp(k)
1974 cgrad                  ghalfm=0.5D0*gggm(k)
1975                   gacontp_hb1(k,num_conti,i)=!ghalfp
1976      &              +(ecosap*(dc_norm(k,j)-cosa*dc_norm(k,i))
1977      &              + ecosbp*(erij(k)-cosb*dc_norm(k,i)))*vbld_inv(i+1)
1978                   gacontp_hb2(k,num_conti,i)=!ghalfp
1979      &              +(ecosap*(dc_norm(k,i)-cosa*dc_norm(k,j))
1980      &              + ecosgp*(erij(k)-cosg*dc_norm(k,j)))*vbld_inv(j+1)
1981                   gacontp_hb3(k,num_conti,i)=gggp(k)
1982                   gacontm_hb1(k,num_conti,i)=!ghalfm
1983      &              +(ecosam*(dc_norm(k,j)-cosa*dc_norm(k,i))
1984      &              + ecosbm*(erij(k)-cosb*dc_norm(k,i)))*vbld_inv(i+1)
1985                   gacontm_hb2(k,num_conti,i)=!ghalfm
1986      &              +(ecosam*(dc_norm(k,i)-cosa*dc_norm(k,j))
1987      &              + ecosgm*(erij(k)-cosg*dc_norm(k,j)))*vbld_inv(j+1)
1988                   gacontm_hb3(k,num_conti,i)=gggm(k)
1989                 enddo
1990               ENDIF ! wcorr
1991               endif  ! num_conti.le.maxconts
1992             endif  ! fcont.gt.0
1993           endif    ! j.gt.i+1
1994           if (wturn3.gt.0.0d0 .or. wturn4.gt.0.0d0) then
1995             do k=1,4
1996               do l=1,3
1997                 ghalf=0.5d0*agg(l,k)
1998                 aggi(l,k)=aggi(l,k)+ghalf
1999                 aggi1(l,k)=aggi1(l,k)+agg(l,k)
2000                 aggj(l,k)=aggj(l,k)+ghalf
2001               enddo
2002             enddo
2003             if (j.eq.nres-1 .and. i.lt.j-2) then
2004               do k=1,4
2005                 do l=1,3
2006                   aggj1(l,k)=aggj1(l,k)+agg(l,k)
2007                 enddo
2008               enddo
2009             endif
2010           endif
2011 c          t_eelecij=t_eelecij+MPI_Wtime()-time00
2012       return
2013       end
2014 C-----------------------------------------------------------------------
2015       subroutine evdwpp_short(evdw1)
2016 C
2017 C Compute Evdwpp
2018 C
2019       implicit real*8 (a-h,o-z)
2020       include 'DIMENSIONS'
2021       include 'COMMON.CONTROL'
2022       include 'COMMON.IOUNITS'
2023       include 'COMMON.GEO'
2024       include 'COMMON.VAR'
2025       include 'COMMON.LOCAL'
2026       include 'COMMON.CHAIN'
2027       include 'COMMON.DERIV'
2028       include 'COMMON.INTERACT'
2029       include 'COMMON.CONTACTS'
2030       include 'COMMON.TORSION'
2031       include 'COMMON.VECTORS'
2032       include 'COMMON.FFIELD'
2033       dimension ggg(3)
2034 c 4/26/02 - AL scaling factor for 1,4 repulsive VDW interactions
2035 #ifdef MOMENT
2036       double precision scal_el /1.0d0/
2037 #else
2038       double precision scal_el /0.5d0/
2039 #endif
2040       evdw1=0.0D0
2041 c      write (iout,*) "iatel_s_vdw",iatel_s_vdw,
2042 c     & " iatel_e_vdw",iatel_e_vdw
2043       call flush(iout)
2044       do i=iatel_s_vdw,iatel_e_vdw
2045         if (itype(i).eq.ntyp1.or. itype(i+1).eq.ntyp1) cycle
2046         dxi=dc(1,i)
2047         dyi=dc(2,i)
2048         dzi=dc(3,i)
2049         dx_normi=dc_norm(1,i)
2050         dy_normi=dc_norm(2,i)
2051         dz_normi=dc_norm(3,i)
2052         xmedi=c(1,i)+0.5d0*dxi
2053         ymedi=c(2,i)+0.5d0*dyi
2054         zmedi=c(3,i)+0.5d0*dzi
2055         num_conti=0
2056 c        write (iout,*) 'i',i,' ielstart',ielstart_vdw(i),
2057 c     &   ' ielend',ielend_vdw(i)
2058         call flush(iout)
2059         do j=ielstart_vdw(i),ielend_vdw(i)
2060           if (itype(j).eq.ntyp1 .or. itype(j+1).eq.ntyp1) cycle
2061           ind=ind+1
2062           iteli=itel(i)
2063           itelj=itel(j)
2064           if (j.eq.i+2 .and. itelj.eq.2) iteli=2
2065           aaa=app(iteli,itelj)
2066           bbb=bpp(iteli,itelj)
2067           dxj=dc(1,j)
2068           dyj=dc(2,j)
2069           dzj=dc(3,j)
2070           dx_normj=dc_norm(1,j)
2071           dy_normj=dc_norm(2,j)
2072           dz_normj=dc_norm(3,j)
2073           xj=c(1,j)+0.5D0*dxj-xmedi
2074           yj=c(2,j)+0.5D0*dyj-ymedi
2075           zj=c(3,j)+0.5D0*dzj-zmedi
2076           rij=xj*xj+yj*yj+zj*zj
2077           rrmij=1.0D0/rij
2078           rij=dsqrt(rij)
2079           sss=sscale(rij/rpp(iteli,itelj))
2080           if (sss.gt.0.0d0) then
2081             rmij=1.0D0/rij
2082             r3ij=rrmij*rmij
2083             r6ij=r3ij*r3ij
2084             ev1=aaa*r6ij*r6ij
2085 c 4/26/02 - AL scaling down 1,4 repulsive VDW interactions
2086             if (j.eq.i+2) ev1=scal_el*ev1
2087             ev2=bbb*r6ij
2088             evdwij=ev1+ev2
2089             if (energy_dec) then
2090               write (iout,'(a6,2i5,0pf7.3,f7.3)') 'evdw1',i,j,evdwij,sss
2091             endif
2092             evdw1=evdw1+evdwij*sss
2093 C
2094 C Calculate contributions to the Cartesian gradient.
2095 C
2096             facvdw=-6*rrmij*(ev1+evdwij)*sss
2097             ggg(1)=facvdw*xj
2098             ggg(2)=facvdw*yj
2099             ggg(3)=facvdw*zj
2100             do k=1,3
2101               gvdwpp(k,j)=gvdwpp(k,j)+ggg(k)
2102               gvdwpp(k,i)=gvdwpp(k,i)-ggg(k)
2103             enddo
2104           endif
2105         enddo ! j
2106       enddo   ! i
2107       return
2108       end
2109 C-----------------------------------------------------------------------------
2110       subroutine escp_long(evdw2,evdw2_14)
2111 C
2112 C This subroutine calculates the excluded-volume interaction energy between
2113 C peptide-group centers and side chains and its gradient in virtual-bond and
2114 C side-chain vectors.
2115 C
2116       implicit real*8 (a-h,o-z)
2117       include 'DIMENSIONS'
2118       include 'COMMON.GEO'
2119       include 'COMMON.VAR'
2120       include 'COMMON.LOCAL'
2121       include 'COMMON.CHAIN'
2122       include 'COMMON.DERIV'
2123       include 'COMMON.INTERACT'
2124       include 'COMMON.FFIELD'
2125       include 'COMMON.IOUNITS'
2126       include 'COMMON.CONTROL'
2127       dimension ggg(3)
2128       evdw2=0.0D0
2129       evdw2_14=0.0d0
2130 cd    print '(a)','Enter ESCP'
2131 cd    write (iout,*) 'iatscp_s=',iatscp_s,' iatscp_e=',iatscp_e
2132       do i=iatscp_s,iatscp_e
2133         if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1) cycle
2134         iteli=itel(i)
2135         xi=0.5D0*(c(1,i)+c(1,i+1))
2136         yi=0.5D0*(c(2,i)+c(2,i+1))
2137         zi=0.5D0*(c(3,i)+c(3,i+1))
2138
2139         do iint=1,nscp_gr(i)
2140
2141         do j=iscpstart(i,iint),iscpend(i,iint)
2142           itypj=itype(j)
2143           if (itypj.eq.ntyp1) cycle
2144 C Uncomment following three lines for SC-p interactions
2145 c         xj=c(1,nres+j)-xi
2146 c         yj=c(2,nres+j)-yi
2147 c         zj=c(3,nres+j)-zi
2148 C Uncomment following three lines for Ca-p interactions
2149           xj=c(1,j)-xi
2150           yj=c(2,j)-yi
2151           zj=c(3,j)-zi
2152           rrij=1.0D0/(xj*xj+yj*yj+zj*zj)
2153
2154           sss=sscale(1.0d0/(dsqrt(rrij)*rscp(itypj,iteli)))
2155
2156           if (sss.lt.1.0d0) then
2157
2158             fac=rrij**expon2
2159             e1=fac*fac*aad(itypj,iteli)
2160             e2=fac*bad(itypj,iteli)
2161             if (iabs(j-i) .le. 2) then
2162               e1=scal14*e1
2163               e2=scal14*e2
2164               evdw2_14=evdw2_14+(e1+e2)*(1.0d0-sss)
2165             endif
2166             evdwij=e1+e2
2167             evdw2=evdw2+evdwij*(1.0d0-sss)
2168             if (energy_dec) write (iout,'(a6,2i5,2(0pf7.3))')
2169      &          'evdw2',i,j,sss,evdwij
2170 C
2171 C Calculate contributions to the gradient in the virtual-bond and SC vectors.
2172 C
2173             fac=-(evdwij+e1)*rrij*(1.0d0-sss)
2174             ggg(1)=xj*fac
2175             ggg(2)=yj*fac
2176             ggg(3)=zj*fac
2177 C Uncomment following three lines for SC-p interactions
2178 c           do k=1,3
2179 c             gradx_scp(k,j)=gradx_scp(k,j)+ggg(k)
2180 c           enddo
2181 C Uncomment following line for SC-p interactions
2182 c             gradx_scp(k,j)=gradx_scp(k,j)+ggg(k)
2183             do k=1,3
2184               gvdwc_scpp(k,i)=gvdwc_scpp(k,i)-ggg(k)
2185               gvdwc_scp(k,j)=gvdwc_scp(k,j)+ggg(k)
2186             enddo
2187           endif
2188         enddo
2189
2190         enddo ! iint
2191       enddo ! i
2192       do i=1,nct
2193         do j=1,3
2194           gvdwc_scp(j,i)=expon*gvdwc_scp(j,i)
2195           gvdwc_scpp(j,i)=expon*gvdwc_scpp(j,i)
2196           gradx_scp(j,i)=expon*gradx_scp(j,i)
2197         enddo
2198       enddo
2199 C******************************************************************************
2200 C
2201 C                              N O T E !!!
2202 C
2203 C To save time the factor EXPON has been extracted from ALL components
2204 C of GVDWC and GRADX. Remember to multiply them by this factor before further
2205 C use!
2206 C
2207 C******************************************************************************
2208       return
2209       end
2210 C-----------------------------------------------------------------------------
2211       subroutine escp_short(evdw2,evdw2_14)
2212 C
2213 C This subroutine calculates the excluded-volume interaction energy between
2214 C peptide-group centers and side chains and its gradient in virtual-bond and
2215 C side-chain vectors.
2216 C
2217       implicit real*8 (a-h,o-z)
2218       include 'DIMENSIONS'
2219       include 'COMMON.GEO'
2220       include 'COMMON.VAR'
2221       include 'COMMON.LOCAL'
2222       include 'COMMON.CHAIN'
2223       include 'COMMON.DERIV'
2224       include 'COMMON.INTERACT'
2225       include 'COMMON.FFIELD'
2226       include 'COMMON.IOUNITS'
2227       include 'COMMON.CONTROL'
2228       dimension ggg(3)
2229       evdw2=0.0D0
2230       evdw2_14=0.0d0
2231 cd    print '(a)','Enter ESCP'
2232 cd    write (iout,*) 'iatscp_s=',iatscp_s,' iatscp_e=',iatscp_e
2233       do i=iatscp_s,iatscp_e
2234         if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1) cycle
2235         iteli=itel(i)
2236         xi=0.5D0*(c(1,i)+c(1,i+1))
2237         yi=0.5D0*(c(2,i)+c(2,i+1))
2238         zi=0.5D0*(c(3,i)+c(3,i+1))
2239
2240         do iint=1,nscp_gr(i)
2241
2242         do j=iscpstart(i,iint),iscpend(i,iint)
2243           itypj=itype(j)
2244           if (itypj.eq.ntyp1) cycle
2245 C Uncomment following three lines for SC-p interactions
2246 c         xj=c(1,nres+j)-xi
2247 c         yj=c(2,nres+j)-yi
2248 c         zj=c(3,nres+j)-zi
2249 C Uncomment following three lines for Ca-p interactions
2250           xj=c(1,j)-xi
2251           yj=c(2,j)-yi
2252           zj=c(3,j)-zi
2253           rrij=1.0D0/(xj*xj+yj*yj+zj*zj)
2254
2255           sss=sscale(1.0d0/(dsqrt(rrij)*rscp(itypj,iteli)))
2256
2257           if (sss.gt.0.0d0) then
2258
2259             fac=rrij**expon2
2260             e1=fac*fac*aad(itypj,iteli)
2261             e2=fac*bad(itypj,iteli)
2262             if (iabs(j-i) .le. 2) then
2263               e1=scal14*e1
2264               e2=scal14*e2
2265               evdw2_14=evdw2_14+(e1+e2)*sss
2266             endif
2267             evdwij=e1+e2
2268             evdw2=evdw2+evdwij*sss
2269             if (energy_dec) write (iout,'(a6,2i5,2(0pf7.3))')
2270      &          'evdw2',i,j,sss,evdwij
2271 C
2272 C Calculate contributions to the gradient in the virtual-bond and SC vectors.
2273 C
2274             fac=-(evdwij+e1)*rrij*sss
2275             ggg(1)=xj*fac
2276             ggg(2)=yj*fac
2277             ggg(3)=zj*fac
2278 C Uncomment following three lines for SC-p interactions
2279 c           do k=1,3
2280 c             gradx_scp(k,j)=gradx_scp(k,j)+ggg(k)
2281 c           enddo
2282 C Uncomment following line for SC-p interactions
2283 c             gradx_scp(k,j)=gradx_scp(k,j)+ggg(k)
2284             do k=1,3
2285               gvdwc_scpp(k,i)=gvdwc_scpp(k,i)-ggg(k)
2286               gvdwc_scp(k,j)=gvdwc_scp(k,j)+ggg(k)
2287             enddo
2288           endif
2289         enddo
2290
2291         enddo ! iint
2292       enddo ! i
2293       do i=1,nct
2294         do j=1,3
2295           gvdwc_scp(j,i)=expon*gvdwc_scp(j,i)
2296           gvdwc_scpp(j,i)=expon*gvdwc_scpp(j,i)
2297           gradx_scp(j,i)=expon*gradx_scp(j,i)
2298         enddo
2299       enddo
2300 C******************************************************************************
2301 C
2302 C                              N O T E !!!
2303 C
2304 C To save time the factor EXPON has been extracted from ALL components
2305 C of GVDWC and GRADX. Remember to multiply them by this factor before further
2306 C use!
2307 C
2308 C******************************************************************************
2309       return
2310       end