source/unres/src_MD/src/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         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,evdw_p,evdw_m)
 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 ccccc      energy_dec=.false.
 590 c     print *,'Entering EGB nnt=',nnt,' nct=',nct,' expon=',expon
 591       evdw=0.0D0
 592       evdw_p=0.0D0
 593       evdw_m=0.0D0
 594       lprn=.false.
 595 c     if (icall.eq.0) lprn=.false.
 596       ind=0
 597       do i=iatsc_s,iatsc_e
 598         itypi=itype(i)
 599         itypi1=itype(i+1)
 600         xi=c(1,nres+i)
 601         yi=c(2,nres+i)
 602         zi=c(3,nres+i)
 603         dxi=dc_norm(1,nres+i)
 604         dyi=dc_norm(2,nres+i)
 605         dzi=dc_norm(3,nres+i)
 606 c        dsci_inv=dsc_inv(itypi)
 607         dsci_inv=vbld_inv(i+nres)
 608 c        write (iout,*) "i",i,dsc_inv(itypi),dsci_inv,1.0d0/vbld(i+nres)
 609 c        write (iout,*) "dcnori",dxi*dxi+dyi*dyi+dzi*dzi
 610 C
 611 C Calculate SC interaction energy.
 612 C
 613         do iint=1,nint_gr(i)
 614           do j=istart(i,iint),iend(i,iint)
 615             ind=ind+1
 616             itypj=itype(j)
 617 c            dscj_inv=dsc_inv(itypj)
 618             dscj_inv=vbld_inv(j+nres)
 619 c            write (iout,*) "j",j,dsc_inv(itypj),dscj_inv,
 620 c     &       1.0d0/vbld(j+nres)
 621 c            write (iout,*) "i",i," j", j," itype",itype(i),itype(j)
 622             sig0ij=sigma(itypi,itypj)
 623             chi1=chi(itypi,itypj)
 624             chi2=chi(itypj,itypi)
 625             chi12=chi1*chi2
 626             chip1=chip(itypi)
 627             chip2=chip(itypj)
 628             chip12=chip1*chip2
 629             alf1=alp(itypi)
 630             alf2=alp(itypj)
 631             alf12=0.5D0*(alf1+alf2)
 632             xj=c(1,nres+j)-xi
 633             yj=c(2,nres+j)-yi
 634             zj=c(3,nres+j)-zi
 635             dxj=dc_norm(1,nres+j)
 636             dyj=dc_norm(2,nres+j)
 637             dzj=dc_norm(3,nres+j)
 638             rrij=1.0D0/(xj*xj+yj*yj+zj*zj)
 639             rij=dsqrt(rrij)
 640             sss=sscale(1.0d0/(rij*sigmaii(itypi,itypj)))
 641
 642             if (sss.lt.1.0d0) then
 643
 644 C Calculate angle-dependent terms of energy and contributions to their
 645 C derivatives.
 646               call sc_angular
 647               sigsq=1.0D0/sigsq
 648               sig=sig0ij*dsqrt(sigsq)
 649               rij_shift=1.0D0/rij-sig+sig0ij
 650 c for diagnostics; uncomment
 651 c              rij_shift=1.2*sig0ij
 652 C I hate to put IF's in the loops, but here don't have another choice!!!!
 653               if (rij_shift.le.0.0D0) then
 654                 evdw=1.0D20
 655 cd                write (iout,'(2(a3,i3,2x),17(0pf7.3))')
 656 cd     &          restyp(itypi),i,restyp(itypj),j,
 657 cd     &          rij_shift,1.0D0/rij,sig,sig0ij,sigsq,1-dsqrt(sigsq)
 658                 return
 659               endif
 660               sigder=-sig*sigsq
 661 c---------------------------------------------------------------
 662               rij_shift=1.0D0/rij_shift
 663               fac=rij_shift**expon
 664               e1=fac*fac*aa(itypi,itypj)
 665               e2=fac*bb(itypi,itypj)
 666               evdwij=eps1*eps2rt*eps3rt*(e1+e2)
 667               eps2der=evdwij*eps3rt
 668               eps3der=evdwij*eps2rt
 669 c              write (iout,*) "sigsq",sigsq," sig",sig," eps2rt",eps2rt,
 670 c     &        " eps3rt",eps3rt," eps1",eps1," e1",e1," e2",e2
 671               evdwij=evdwij*eps2rt*eps3rt
 672 #ifdef TSCSC
 673               if (bb(itypi,itypj).gt.0) then
 674                  evdw_p=evdw_p+evdwij*(1.0d0-sss)
 675               else
 676                  evdw_m=evdw_m+evdwij*(1.0d0-sss)
 677               endif
 678 #else
 679               evdw=evdw+evdwij*(1.0d0-sss)
 680 #endif
 681               if (lprn) then
 682               sigm=dabs(aa(itypi,itypj)/bb(itypi,itypj))**(1.0D0/6.0D0)
 683               epsi=bb(itypi,itypj)**2/aa(itypi,itypj)
 684               write (iout,'(2(a3,i3,2x),17(0pf7.3))')
 685      &          restyp(itypi),i,restyp(itypj),j,
 686      &          epsi,sigm,chi1,chi2,chip1,chip2,
 687      &          eps1,eps2rt**2,eps3rt**2,sig,sig0ij,
 688      &          om1,om2,om12,1.0D0/rij,1.0D0/rij_shift,
 689      &          evdwij
 690               endif
 691
 692               if (energy_dec) write (iout,'(a6,2i5,0pf7.3)')
 693      &                        'evdw',i,j,evdwij
 694
 695 C Calculate gradient components.
 696               e1=e1*eps1*eps2rt**2*eps3rt**2
 697               fac=-expon*(e1+evdwij)*rij_shift
 698               sigder=fac*sigder
 699               fac=rij*fac
 700 c              fac=0.0d0
 701 C Calculate the radial part of the gradient
 702               gg(1)=xj*fac
 703               gg(2)=yj*fac
 704               gg(3)=zj*fac
 705 C Calculate angular part of the gradient.
 706 #ifdef TSCSC
 707               if (bb(itypi,itypj).gt.0) then
 708                call sc_grad_scale_T(1.0d0-sss)
 709               else
 710                call sc_grad_scale(1.0d0-sss)
 711               endif
 712 #else
 713               call sc_grad_scale(1.0d0-sss)
 714 #endif
 715             endif
 716           enddo      ! j
 717         enddo        ! iint
 718       enddo          ! i
 719 c      write (iout,*) "Number of loop steps in EGB:",ind
 720 cccc      energy_dec=.false.
 721       return
 722       end
 723 C-----------------------------------------------------------------------------
 724       subroutine egb_short(evdw,evdw_p,evdw_m)
 725 C
 726 C This subroutine calculates the interaction energy of nonbonded side chains
 727 C assuming the Gay-Berne potential of interaction.
 728 C
 729       implicit real*8 (a-h,o-z)
 730       include 'DIMENSIONS'
 731       include 'COMMON.GEO'
 732       include 'COMMON.VAR'
 733       include 'COMMON.LOCAL'
 734       include 'COMMON.CHAIN'
 735       include 'COMMON.DERIV'
 736       include 'COMMON.NAMES'
 737       include 'COMMON.INTERACT'
 738       include 'COMMON.IOUNITS'
 739       include 'COMMON.CALC'
 740       include 'COMMON.CONTROL'
 741       logical lprn
 742       evdw=0.0D0
 743       evdw_p=0.0D0
 744       evdw_m=0.0D0
 745 ccccc      energy_dec=.false.
 746 c     print *,'Entering EGB nnt=',nnt,' nct=',nct,' expon=',expon
 747       evdw=0.0D0
 748       lprn=.false.
 749 c     if (icall.eq.0) lprn=.false.
 750       ind=0
 751       do i=iatsc_s,iatsc_e
 752         itypi=itype(i)
 753         itypi1=itype(i+1)
 754         xi=c(1,nres+i)
 755         yi=c(2,nres+i)
 756         zi=c(3,nres+i)
 757         dxi=dc_norm(1,nres+i)
 758         dyi=dc_norm(2,nres+i)
 759         dzi=dc_norm(3,nres+i)
 760 c        dsci_inv=dsc_inv(itypi)
 761         dsci_inv=vbld_inv(i+nres)
 762 c        write (iout,*) "i",i,dsc_inv(itypi),dsci_inv,1.0d0/vbld(i+nres)
 763 c        write (iout,*) "dcnori",dxi*dxi+dyi*dyi+dzi*dzi
 764 C
 765 C Calculate SC interaction energy.
 766 C
 767         do iint=1,nint_gr(i)
 768           do j=istart(i,iint),iend(i,iint)
 769             ind=ind+1
 770             itypj=itype(j)
 771 c            dscj_inv=dsc_inv(itypj)
 772             dscj_inv=vbld_inv(j+nres)
 773 c            write (iout,*) "j",j,dsc_inv(itypj),dscj_inv,
 774 c     &       1.0d0/vbld(j+nres)
 775 c            write (iout,*) "i",i," j", j," itype",itype(i),itype(j)
 776             sig0ij=sigma(itypi,itypj)
 777             chi1=chi(itypi,itypj)
 778             chi2=chi(itypj,itypi)
 779             chi12=chi1*chi2
 780             chip1=chip(itypi)
 781             chip2=chip(itypj)
 782             chip12=chip1*chip2
 783             alf1=alp(itypi)
 784             alf2=alp(itypj)
 785             alf12=0.5D0*(alf1+alf2)
 786             xj=c(1,nres+j)-xi
 787             yj=c(2,nres+j)-yi
 788             zj=c(3,nres+j)-zi
 789             dxj=dc_norm(1,nres+j)
 790             dyj=dc_norm(2,nres+j)
 791             dzj=dc_norm(3,nres+j)
 792             rrij=1.0D0/(xj*xj+yj*yj+zj*zj)
 793             rij=dsqrt(rrij)
 794             sss=sscale(1.0d0/(rij*sigmaii(itypi,itypj)))
 795
 796             if (sss.gt.0.0d0) then
 797
 798 C Calculate angle-dependent terms of energy and contributions to their
 799 C derivatives.
 800               call sc_angular
 801               sigsq=1.0D0/sigsq
 802               sig=sig0ij*dsqrt(sigsq)
 803               rij_shift=1.0D0/rij-sig+sig0ij
 804 c for diagnostics; uncomment
 805 c              rij_shift=1.2*sig0ij
 806 C I hate to put IF's in the loops, but here don't have another choice!!!!
 807               if (rij_shift.le.0.0D0) then
 808                 evdw=1.0D20
 809 cd                write (iout,'(2(a3,i3,2x),17(0pf7.3))')
 810 cd     &          restyp(itypi),i,restyp(itypj),j,
 811 cd     &          rij_shift,1.0D0/rij,sig,sig0ij,sigsq,1-dsqrt(sigsq)
 812                 return
 813               endif
 814               sigder=-sig*sigsq
 815 c---------------------------------------------------------------
 816               rij_shift=1.0D0/rij_shift
 817               fac=rij_shift**expon
 818               e1=fac*fac*aa(itypi,itypj)
 819               e2=fac*bb(itypi,itypj)
 820               evdwij=eps1*eps2rt*eps3rt*(e1+e2)
 821               eps2der=evdwij*eps3rt
 822               eps3der=evdwij*eps2rt
 823 c              write (iout,*) "sigsq",sigsq," sig",sig," eps2rt",eps2rt,
 824 c     &        " eps3rt",eps3rt," eps1",eps1," e1",e1," e2",e2
 825               evdwij=evdwij*eps2rt*eps3rt
 826 #ifdef TSCSC
 827               if (bb(itypi,itypj).gt.0) then
 828                  evdw_p=evdw_p+evdwij*sss
 829               else
 830                  evdw_m=evdw_m+evdwij*sss
 831               endif
 832 #else
 833               evdw=evdw+evdwij*sss
 834 #endif
 835               if (lprn) then
 836               sigm=dabs(aa(itypi,itypj)/bb(itypi,itypj))**(1.0D0/6.0D0)
 837               epsi=bb(itypi,itypj)**2/aa(itypi,itypj)
 838               write (iout,'(2(a3,i3,2x),17(0pf7.3))')
 839      &          restyp(itypi),i,restyp(itypj),j,
 840      &          epsi,sigm,chi1,chi2,chip1,chip2,
 841      &          eps1,eps2rt**2,eps3rt**2,sig,sig0ij,
 842      &          om1,om2,om12,1.0D0/rij,1.0D0/rij_shift,
 843      &          evdwij
 844               endif
 845
 846               if (energy_dec) write (iout,'(a6,2i5,0pf7.3)')
 847      &                        'evdw',i,j,evdwij
 848
 849 C Calculate gradient components.
 850               e1=e1*eps1*eps2rt**2*eps3rt**2
 851               fac=-expon*(e1+evdwij)*rij_shift
 852               sigder=fac*sigder
 853               fac=rij*fac
 854 c              fac=0.0d0
 855 C Calculate the radial part of the gradient
 856               gg(1)=xj*fac
 857               gg(2)=yj*fac
 858               gg(3)=zj*fac
 859 C Calculate angular part of the gradient.
 860 #ifdef TSCSC
 861               if (bb(itypi,itypj).gt.0) then
 862                call sc_grad_scale_T(sss)
 863               else
 864                call sc_grad_scale(sss)
 865               endif
 866 #else
 867               call sc_grad_scale(sss)
 868 #endif
 869             endif
 870           enddo      ! j
 871         enddo        ! iint
 872       enddo          ! i
 873 c      write (iout,*) "Number of loop steps in EGB:",ind
 874 cccc      energy_dec=.false.
 875       return
 876       end
 877 C-----------------------------------------------------------------------------
 878       subroutine egbv_long(evdw)
 879 C
 880 C This subroutine calculates the interaction energy of nonbonded side chains
 881 C assuming the Gay-Berne-Vorobjev potential of interaction.
 882 C
 883       implicit real*8 (a-h,o-z)
 884       include 'DIMENSIONS'
 885       include 'COMMON.GEO'
 886       include 'COMMON.VAR'
 887       include 'COMMON.LOCAL'
 888       include 'COMMON.CHAIN'
 889       include 'COMMON.DERIV'
 890       include 'COMMON.NAMES'
 891       include 'COMMON.INTERACT'
 892       include 'COMMON.IOUNITS'
 893       include 'COMMON.CALC'
 894       common /srutu/ icall
 895       logical lprn
 896       evdw=0.0D0
 897 c     print *,'Entering EGB nnt=',nnt,' nct=',nct,' expon=',expon
 898       evdw=0.0D0
 899       lprn=.false.
 900 c     if (icall.eq.0) lprn=.true.
 901       ind=0
 902       do i=iatsc_s,iatsc_e
 903         itypi=itype(i)
 904         itypi1=itype(i+1)
 905         xi=c(1,nres+i)
 906         yi=c(2,nres+i)
 907         zi=c(3,nres+i)
 908         dxi=dc_norm(1,nres+i)
 909         dyi=dc_norm(2,nres+i)
 910         dzi=dc_norm(3,nres+i)
 911 c        dsci_inv=dsc_inv(itypi)
 912         dsci_inv=vbld_inv(i+nres)
 913 C
 914 C Calculate SC interaction energy.
 915 C
 916         do iint=1,nint_gr(i)
 917           do j=istart(i,iint),iend(i,iint)
 918             ind=ind+1
 919             itypj=itype(j)
 920 c            dscj_inv=dsc_inv(itypj)
 921             dscj_inv=vbld_inv(j+nres)
 922             sig0ij=sigma(itypi,itypj)
 923             r0ij=r0(itypi,itypj)
 924             chi1=chi(itypi,itypj)
 925             chi2=chi(itypj,itypi)
 926             chi12=chi1*chi2
 927             chip1=chip(itypi)
 928             chip2=chip(itypj)
 929             chip12=chip1*chip2
 930             alf1=alp(itypi)
 931             alf2=alp(itypj)
 932             alf12=0.5D0*(alf1+alf2)
 933             xj=c(1,nres+j)-xi
 934             yj=c(2,nres+j)-yi
 935             zj=c(3,nres+j)-zi
 936             dxj=dc_norm(1,nres+j)
 937             dyj=dc_norm(2,nres+j)
 938             dzj=dc_norm(3,nres+j)
 939             rrij=1.0D0/(xj*xj+yj*yj+zj*zj)
 940             rij=dsqrt(rrij)
 941
 942             sss=sscale(1.0d0/(rij*sigmaii(itypi,itypj)))
 943
 944             if (sss.lt.1.0d0) then
 945
 946 C Calculate angle-dependent terms of energy and contributions to their
 947 C derivatives.
 948               call sc_angular
 949               sigsq=1.0D0/sigsq
 950               sig=sig0ij*dsqrt(sigsq)
 951               rij_shift=1.0D0/rij-sig+r0ij
 952 C I hate to put IF's in the loops, but here don't have another choice!!!!
 953               if (rij_shift.le.0.0D0) then
 954                 evdw=1.0D20
 955                 return
 956               endif
 957               sigder=-sig*sigsq
 958 c---------------------------------------------------------------
 959               rij_shift=1.0D0/rij_shift
 960               fac=rij_shift**expon
 961               e1=fac*fac*aa(itypi,itypj)
 962               e2=fac*bb(itypi,itypj)
 963               evdwij=eps1*eps2rt*eps3rt*(e1+e2)
 964               eps2der=evdwij*eps3rt
 965               eps3der=evdwij*eps2rt
 966               fac_augm=rrij**expon
 967               e_augm=augm(itypi,itypj)*fac_augm
 968               evdwij=evdwij*eps2rt*eps3rt
 969               evdw=evdw+(evdwij+e_augm)*(1.0d0-sss)
 970               if (lprn) then
 971               sigm=dabs(aa(itypi,itypj)/bb(itypi,itypj))**(1.0D0/6.0D0)
 972               epsi=bb(itypi,itypj)**2/aa(itypi,itypj)
 973               write (iout,'(2(a3,i3,2x),17(0pf7.3))')
 974      &          restyp(itypi),i,restyp(itypj),j,
 975      &          epsi,sigm,sig,(augm(itypi,itypj)/epsi)**(1.0D0/12.0D0),
 976      &          chi1,chi2,chip1,chip2,
 977      &          eps1,eps2rt**2,eps3rt**2,
 978      &          om1,om2,om12,1.0D0/rij,1.0D0/rij_shift,
 979      &          evdwij+e_augm
 980               endif
 981 C Calculate gradient components.
 982               e1=e1*eps1*eps2rt**2*eps3rt**2
 983               fac=-expon*(e1+evdwij)*rij_shift
 984               sigder=fac*sigder
 985               fac=rij*fac-2*expon*rrij*e_augm
 986 C Calculate the radial part of the gradient
 987               gg(1)=xj*fac
 988               gg(2)=yj*fac
 989               gg(3)=zj*fac
 990 C Calculate angular part of the gradient.
 991               call sc_grad_scale(1.0d0-sss)
 992             endif
 993           enddo      ! j
 994         enddo        ! iint
 995       enddo          ! i
 996       end
 997 C-----------------------------------------------------------------------------
 998       subroutine egbv_short(evdw)
 999 C
1000 C This subroutine calculates the interaction energy of nonbonded side chains
1001 C assuming the Gay-Berne-Vorobjev potential of interaction.
1002 C
1003       implicit real*8 (a-h,o-z)
1004       include 'DIMENSIONS'
1005       include 'COMMON.GEO'
1006       include 'COMMON.VAR'
1007       include 'COMMON.LOCAL'
1008       include 'COMMON.CHAIN'
1009       include 'COMMON.DERIV'
1010       include 'COMMON.NAMES'
1011       include 'COMMON.INTERACT'
1012       include 'COMMON.IOUNITS'
1013       include 'COMMON.CALC'
1014       common /srutu/ icall
1015       logical lprn
1016       evdw=0.0D0
1017 c     print *,'Entering EGB nnt=',nnt,' nct=',nct,' expon=',expon
1018       evdw=0.0D0
1019       lprn=.false.
1020 c     if (icall.eq.0) lprn=.true.
1021       ind=0
1022       do i=iatsc_s,iatsc_e
1023         itypi=itype(i)
1024         itypi1=itype(i+1)
1025         xi=c(1,nres+i)
1026         yi=c(2,nres+i)
1027         zi=c(3,nres+i)
1028         dxi=dc_norm(1,nres+i)
1029         dyi=dc_norm(2,nres+i)
1030         dzi=dc_norm(3,nres+i)
1031 c        dsci_inv=dsc_inv(itypi)
1032         dsci_inv=vbld_inv(i+nres)
1033 C
1034 C Calculate SC interaction energy.
1035 C
1036         do iint=1,nint_gr(i)
1037           do j=istart(i,iint),iend(i,iint)
1038             ind=ind+1
1039             itypj=itype(j)
1040 c            dscj_inv=dsc_inv(itypj)
1041             dscj_inv=vbld_inv(j+nres)
1042             sig0ij=sigma(itypi,itypj)
1043             r0ij=r0(itypi,itypj)
1044             chi1=chi(itypi,itypj)
1045             chi2=chi(itypj,itypi)
1046             chi12=chi1*chi2
1047             chip1=chip(itypi)
1048             chip2=chip(itypj)
1049             chip12=chip1*chip2
1050             alf1=alp(itypi)
1051             alf2=alp(itypj)
1052             alf12=0.5D0*(alf1+alf2)
1053             xj=c(1,nres+j)-xi
1054             yj=c(2,nres+j)-yi
1055             zj=c(3,nres+j)-zi
1056             dxj=dc_norm(1,nres+j)
1057             dyj=dc_norm(2,nres+j)
1058             dzj=dc_norm(3,nres+j)
1059             rrij=1.0D0/(xj*xj+yj*yj+zj*zj)
1060             rij=dsqrt(rrij)
1061
1062             sss=sscale(1.0d0/(rij*sigmaii(itypi,itypj)))
1063
1064             if (sss.gt.0.0d0) then
1065
1066 C Calculate angle-dependent terms of energy and contributions to their
1067 C derivatives.
1068               call sc_angular
1069               sigsq=1.0D0/sigsq
1070               sig=sig0ij*dsqrt(sigsq)
1071               rij_shift=1.0D0/rij-sig+r0ij
1072 C I hate to put IF's in the loops, but here don't have another choice!!!!
1073               if (rij_shift.le.0.0D0) then
1074                 evdw=1.0D20
1075                 return
1076               endif
1077               sigder=-sig*sigsq
1078 c---------------------------------------------------------------
1079               rij_shift=1.0D0/rij_shift
1080               fac=rij_shift**expon
1081               e1=fac*fac*aa(itypi,itypj)
1082               e2=fac*bb(itypi,itypj)
1083               evdwij=eps1*eps2rt*eps3rt*(e1+e2)
1084               eps2der=evdwij*eps3rt
1085               eps3der=evdwij*eps2rt
1086               fac_augm=rrij**expon
1087               e_augm=augm(itypi,itypj)*fac_augm
1088               evdwij=evdwij*eps2rt*eps3rt
1089               evdw=evdw+(evdwij+e_augm)*sss
1090               if (lprn) then
1091               sigm=dabs(aa(itypi,itypj)/bb(itypi,itypj))**(1.0D0/6.0D0)
1092               epsi=bb(itypi,itypj)**2/aa(itypi,itypj)
1093               write (iout,'(2(a3,i3,2x),17(0pf7.3))')
1094      &          restyp(itypi),i,restyp(itypj),j,
1095      &          epsi,sigm,sig,(augm(itypi,itypj)/epsi)**(1.0D0/12.0D0),
1096      &          chi1,chi2,chip1,chip2,
1097      &          eps1,eps2rt**2,eps3rt**2,
1098      &          om1,om2,om12,1.0D0/rij,1.0D0/rij_shift,
1099      &          evdwij+e_augm
1100               endif
1101 C Calculate gradient components.
1102               e1=e1*eps1*eps2rt**2*eps3rt**2
1103               fac=-expon*(e1+evdwij)*rij_shift
1104               sigder=fac*sigder
1105               fac=rij*fac-2*expon*rrij*e_augm
1106 C Calculate the radial part of the gradient
1107               gg(1)=xj*fac
1108               gg(2)=yj*fac
1109               gg(3)=zj*fac
1110 C Calculate angular part of the gradient.
1111               call sc_grad_scale(sss)
1112             endif
1113           enddo      ! j
1114         enddo        ! iint
1115       enddo          ! i
1116       end
1117 C----------------------------------------------------------------------------
1118       subroutine sc_grad_scale(scalfac)
1119       implicit real*8 (a-h,o-z)
1120       include 'DIMENSIONS'
1121       include 'COMMON.CHAIN'
1122       include 'COMMON.DERIV'
1123       include 'COMMON.CALC'
1124       include 'COMMON.IOUNITS'
1125       double precision dcosom1(3),dcosom2(3)
1126       double precision scalfac
1127       eom1=eps2der*eps2rt_om1-2.0D0*alf1*eps3der+sigder*sigsq_om1
1128       eom2=eps2der*eps2rt_om2+2.0D0*alf2*eps3der+sigder*sigsq_om2
1129       eom12=evdwij*eps1_om12+eps2der*eps2rt_om12
1130      &     -2.0D0*alf12*eps3der+sigder*sigsq_om12
1131 c diagnostics only
1132 c      eom1=0.0d0
1133 c      eom2=0.0d0
1134 c      eom12=evdwij*eps1_om12
1135 c end diagnostics
1136 c      write (iout,*) "eps2der",eps2der," eps3der",eps3der,
1137 c     &  " sigder",sigder
1138 c      write (iout,*) "eps1_om12",eps1_om12," eps2rt_om12",eps2rt_om12
1139 c      write (iout,*) "eom1",eom1," eom2",eom2," eom12",eom12
1140       do k=1,3
1141         dcosom1(k)=rij*(dc_norm(k,nres+i)-om1*erij(k))
1142         dcosom2(k)=rij*(dc_norm(k,nres+j)-om2*erij(k))
1143       enddo
1144       do k=1,3
1145         gg(k)=(gg(k)+eom1*dcosom1(k)+eom2*dcosom2(k))*scalfac
1146       enddo
1147 c      write (iout,*) "gg",(gg(k),k=1,3)
1148       do k=1,3
1149         gvdwx(k,i)=gvdwx(k,i)-gg(k)
1150      &            +(eom12*(dc_norm(k,nres+j)-om12*dc_norm(k,nres+i))
1151      &          +eom1*(erij(k)-om1*dc_norm(k,nres+i)))*dsci_inv*scalfac
1152         gvdwx(k,j)=gvdwx(k,j)+gg(k)
1153      &            +(eom12*(dc_norm(k,nres+i)-om12*dc_norm(k,nres+j))
1154      &          +eom2*(erij(k)-om2*dc_norm(k,nres+j)))*dscj_inv*scalfac
1155 c        write (iout,*)(eom12*(dc_norm(k,nres+j)-om12*dc_norm(k,nres+i))
1156 c     &            +eom1*(erij(k)-om1*dc_norm(k,nres+i)))*dsci_inv
1157 c        write (iout,*)(eom12*(dc_norm(k,nres+i)-om12*dc_norm(k,nres+j))
1158 c     &            +eom2*(erij(k)-om2*dc_norm(k,nres+j)))*dscj_inv
1159       enddo
1160 C
1161 C Calculate the components of the gradient in DC and X
1162 C
1163       do l=1,3
1164         gvdwc(l,i)=gvdwc(l,i)-gg(l)
1165         gvdwc(l,j)=gvdwc(l,j)+gg(l)
1166       enddo
1167       return
1168       end
1169 C----------------------------------------------------------------------------
1170       subroutine sc_grad_scale_T(scalfac)
1171       implicit real*8 (a-h,o-z)
1172       include 'DIMENSIONS'
1173       include 'COMMON.CHAIN'
1174       include 'COMMON.DERIV'
1175       include 'COMMON.CALC'
1176       include 'COMMON.IOUNITS'
1177       double precision dcosom1(3),dcosom2(3)
1178       double precision scalfac
1179       eom1=eps2der*eps2rt_om1-2.0D0*alf1*eps3der+sigder*sigsq_om1
1180       eom2=eps2der*eps2rt_om2+2.0D0*alf2*eps3der+sigder*sigsq_om2
1181       eom12=evdwij*eps1_om12+eps2der*eps2rt_om12
1182      &     -2.0D0*alf12*eps3der+sigder*sigsq_om12
1183 c diagnostics only
1184 c      eom1=0.0d0
1185 c      eom2=0.0d0
1186 c      eom12=evdwij*eps1_om12
1187 c end diagnostics
1188 c      write (iout,*) "eps2der",eps2der," eps3der",eps3der,
1189 c     &  " sigder",sigder
1190 c      write (iout,*) "eps1_om12",eps1_om12," eps2rt_om12",eps2rt_om12
1191 c      write (iout,*) "eom1",eom1," eom2",eom2," eom12",eom12
1192       do k=1,3
1193         dcosom1(k)=rij*(dc_norm(k,nres+i)-om1*erij(k))
1194         dcosom2(k)=rij*(dc_norm(k,nres+j)-om2*erij(k))
1195       enddo
1196       do k=1,3
1197         gg(k)=(gg(k)+eom1*dcosom1(k)+eom2*dcosom2(k))*scalfac
1198       enddo
1199 c      write (iout,*) "gg",(gg(k),k=1,3)
1200       do k=1,3
1201         gvdwxT(k,i)=gvdwxT(k,i)-gg(k)
1202      &            +(eom12*(dc_norm(k,nres+j)-om12*dc_norm(k,nres+i))
1203      &          +eom1*(erij(k)-om1*dc_norm(k,nres+i)))*dsci_inv*scalfac
1204         gvdwxT(k,j)=gvdwxT(k,j)+gg(k)
1205      &            +(eom12*(dc_norm(k,nres+i)-om12*dc_norm(k,nres+j))
1206      &          +eom2*(erij(k)-om2*dc_norm(k,nres+j)))*dscj_inv*scalfac
1207 c        write (iout,*)(eom12*(dc_norm(k,nres+j)-om12*dc_norm(k,nres+i))
1208 c     &            +eom1*(erij(k)-om1*dc_norm(k,nres+i)))*dsci_inv
1209 c        write (iout,*)(eom12*(dc_norm(k,nres+i)-om12*dc_norm(k,nres+j))
1210 c     &            +eom2*(erij(k)-om2*dc_norm(k,nres+j)))*dscj_inv
1211       enddo
1212 C
1213 C Calculate the components of the gradient in DC and X
1214 C
1215       do l=1,3
1216         gvdwcT(l,i)=gvdwcT(l,i)-gg(l)
1217         gvdwcT(l,j)=gvdwcT(l,j)+gg(l)
1218       enddo
1219       return
1220       end
1221
1222 C--------------------------------------------------------------------------
1223       subroutine eelec_scale(ees,evdw1,eel_loc,eello_turn3,eello_turn4)
1224 C
1225 C This subroutine calculates the average interaction energy and its gradient
1226 C in the virtual-bond vectors between non-adjacent peptide groups, based on
1227 C the potential described in Liwo et al., Protein Sci., 1993, 2, 1715.
1228 C The potential depends both on the distance of peptide-group centers and on
1229 C the orientation of the CA-CA virtual bonds.
1230 C
1231       implicit real*8 (a-h,o-z)
1232 #ifdef MPI
1233       include 'mpif.h'
1234 #endif
1235       include 'DIMENSIONS'
1236       include 'COMMON.CONTROL'
1237       include 'COMMON.SETUP'
1238       include 'COMMON.IOUNITS'
1239       include 'COMMON.GEO'
1240       include 'COMMON.VAR'
1241       include 'COMMON.LOCAL'
1242       include 'COMMON.CHAIN'
1243       include 'COMMON.DERIV'
1244       include 'COMMON.INTERACT'
1245       include 'COMMON.CONTACTS'
1246       include 'COMMON.TORSION'
1247       include 'COMMON.VECTORS'
1248       include 'COMMON.FFIELD'
1249       include 'COMMON.TIME1'
1250       dimension ggg(3),gggp(3),gggm(3),erij(3),dcosb(3),dcosg(3),
1251      &          erder(3,3),uryg(3,3),urzg(3,3),vryg(3,3),vrzg(3,3)
1252       double precision acipa(2,2),agg(3,4),aggi(3,4),aggi1(3,4),
1253      &    aggj(3,4),aggj1(3,4),a_temp(2,2),muij(4)
1254       common /locel/ a_temp,agg,aggi,aggi1,aggj,aggj1,a22,a23,a32,a33,
1255      &    dxi,dyi,dzi,dx_normi,dy_normi,dz_normi,xmedi,ymedi,zmedi,
1256      &    num_conti,j1,j2
1257 c 4/26/02 - AL scaling factor for 1,4 repulsive VDW interactions
1258 #ifdef MOMENT
1259       double precision scal_el /1.0d0/
1260 #else
1261       double precision scal_el /0.5d0/
1262 #endif
1263 C 12/13/98
1264 C 13-go grudnia roku pamietnego...
1265       double precision unmat(3,3) /1.0d0,0.0d0,0.0d0,
1266      &                   0.0d0,1.0d0,0.0d0,
1267      &                   0.0d0,0.0d0,1.0d0/
1268 cd      write(iout,*) 'In EELEC'
1269 cd      do i=1,nloctyp
1270 cd        write(iout,*) 'Type',i
1271 cd        write(iout,*) 'B1',B1(:,i)
1272 cd        write(iout,*) 'B2',B2(:,i)
1273 cd        write(iout,*) 'CC',CC(:,:,i)
1274 cd        write(iout,*) 'DD',DD(:,:,i)
1275 cd        write(iout,*) 'EE',EE(:,:,i)
1276 cd      enddo
1277 cd      call check_vecgrad
1278 cd      stop
1279       if (icheckgrad.eq.1) then
1280         do i=1,nres-1
1281           fac=1.0d0/dsqrt(scalar(dc(1,i),dc(1,i)))
1282           do k=1,3
1283             dc_norm(k,i)=dc(k,i)*fac
1284           enddo
1285 c          write (iout,*) 'i',i,' fac',fac
1286         enddo
1287       endif
1288       if (wel_loc.gt.0.0d0 .or. wcorr4.gt.0.0d0 .or. wcorr5.gt.0.0d0
1289      &    .or. wcorr6.gt.0.0d0 .or. wturn3.gt.0.0d0 .or.
1290      &    wturn4.gt.0.0d0 .or. wturn6.gt.0.0d0) then
1291 c        call vec_and_deriv
1292 #ifdef TIMING
1293         time01=MPI_Wtime()
1294 #endif
1295         call set_matrices
1296 #ifdef TIMING
1297         time_mat=time_mat+MPI_Wtime()-time01
1298 #endif
1299       endif
1300 cd      do i=1,nres-1
1301 cd        write (iout,*) 'i=',i
1302 cd        do k=1,3
1303 cd        write (iout,'(i5,2f10.5)') k,uy(k,i),uz(k,i)
1304 cd        enddo
1305 cd        do k=1,3
1306 cd          write (iout,'(f10.5,2x,3f10.5,2x,3f10.5)')
1307 cd     &     uz(k,i),(uzgrad(k,l,1,i),l=1,3),(uzgrad(k,l,2,i),l=1,3)
1308 cd        enddo
1309 cd      enddo
1310       t_eelecij=0.0d0
1311       ees=0.0D0
1312       evdw1=0.0D0
1313       eel_loc=0.0d0
1314       eello_turn3=0.0d0
1315       eello_turn4=0.0d0
1316       ind=0
1317       do i=1,nres
1318         num_cont_hb(i)=0
1319       enddo
1320 cd      print '(a)','Enter EELEC'
1321 cd      write (iout,*) 'iatel_s=',iatel_s,' iatel_e=',iatel_e
1322       do i=1,nres
1323         gel_loc_loc(i)=0.0d0
1324         gcorr_loc(i)=0.0d0
1325       enddo
1326 c
1327 c
1328 c 9/27/08 AL Split the interaction loop to ensure load balancing of turn terms
1329 C
1330 C Loop over i,i+2 and i,i+3 pairs of the peptide groups
1331 C
1332       do i=iturn3_start,iturn3_end
1333         dxi=dc(1,i)
1334         dyi=dc(2,i)
1335         dzi=dc(3,i)
1336         dx_normi=dc_norm(1,i)
1337         dy_normi=dc_norm(2,i)
1338         dz_normi=dc_norm(3,i)
1339         xmedi=c(1,i)+0.5d0*dxi
1340         ymedi=c(2,i)+0.5d0*dyi
1341         zmedi=c(3,i)+0.5d0*dzi
1342         num_conti=0
1343         call eelecij_scale(i,i+2,ees,evdw1,eel_loc)
1344         if (wturn3.gt.0.0d0) call eturn3(i,eello_turn3)
1345         num_cont_hb(i)=num_conti
1346       enddo
1347       do i=iturn4_start,iturn4_end
1348         dxi=dc(1,i)
1349         dyi=dc(2,i)
1350         dzi=dc(3,i)
1351         dx_normi=dc_norm(1,i)
1352         dy_normi=dc_norm(2,i)
1353         dz_normi=dc_norm(3,i)
1354         xmedi=c(1,i)+0.5d0*dxi
1355         ymedi=c(2,i)+0.5d0*dyi
1356         zmedi=c(3,i)+0.5d0*dzi
1357         num_conti=num_cont_hb(i)
1358         call eelecij_scale(i,i+3,ees,evdw1,eel_loc)
1359         if (wturn4.gt.0.0d0) call eturn4(i,eello_turn4)
1360         num_cont_hb(i)=num_conti
1361       enddo   ! i
1362 c
1363 c Loop over all pairs of interacting peptide groups except i,i+2 and i,i+3
1364 c
1365       do i=iatel_s,iatel_e
1366         dxi=dc(1,i)
1367         dyi=dc(2,i)
1368         dzi=dc(3,i)
1369         dx_normi=dc_norm(1,i)
1370         dy_normi=dc_norm(2,i)
1371         dz_normi=dc_norm(3,i)
1372         xmedi=c(1,i)+0.5d0*dxi
1373         ymedi=c(2,i)+0.5d0*dyi
1374         zmedi=c(3,i)+0.5d0*dzi
1375 c        write (iout,*) 'i',i,' ielstart',ielstart(i),' ielend',ielend(i)
1376         num_conti=num_cont_hb(i)
1377         do j=ielstart(i),ielend(i)
1378           call eelecij_scale(i,j,ees,evdw1,eel_loc)
1379         enddo ! j
1380         num_cont_hb(i)=num_conti
1381       enddo   ! i
1382 c      write (iout,*) "Number of loop steps in EELEC:",ind
1383 cd      do i=1,nres
1384 cd        write (iout,'(i3,3f10.5,5x,3f10.5)')
1385 cd     &     i,(gel_loc(k,i),k=1,3),gel_loc_loc(i)
1386 cd      enddo
1387 c 12/7/99 Adam eello_turn3 will be considered as a separate energy term
1388 ccc      eel_loc=eel_loc+eello_turn3
1389 cd      print *,"Processor",fg_rank," t_eelecij",t_eelecij
1390       return
1391       end
1392 C-------------------------------------------------------------------------------
1393       subroutine eelecij_scale(i,j,ees,evdw1,eel_loc)
1394       implicit real*8 (a-h,o-z)
1395       include 'DIMENSIONS'
1396 #ifdef MPI
1397       include "mpif.h"
1398 #endif
1399       include 'COMMON.CONTROL'
1400       include 'COMMON.IOUNITS'
1401       include 'COMMON.GEO'
1402       include 'COMMON.VAR'
1403       include 'COMMON.LOCAL'
1404       include 'COMMON.CHAIN'
1405       include 'COMMON.DERIV'
1406       include 'COMMON.INTERACT'
1407       include 'COMMON.CONTACTS'
1408       include 'COMMON.TORSION'
1409       include 'COMMON.VECTORS'
1410       include 'COMMON.FFIELD'
1411       include 'COMMON.TIME1'
1412       dimension ggg(3),gggp(3),gggm(3),erij(3),dcosb(3),dcosg(3),
1413      &          erder(3,3),uryg(3,3),urzg(3,3),vryg(3,3),vrzg(3,3)
1414       double precision acipa(2,2),agg(3,4),aggi(3,4),aggi1(3,4),
1415      &    aggj(3,4),aggj1(3,4),a_temp(2,2),muij(4)
1416       common /locel/ a_temp,agg,aggi,aggi1,aggj,aggj1,a22,a23,a32,a33,
1417      &    dxi,dyi,dzi,dx_normi,dy_normi,dz_normi,xmedi,ymedi,zmedi,
1418      &    num_conti,j1,j2
1419 c 4/26/02 - AL scaling factor for 1,4 repulsive VDW interactions
1420 #ifdef MOMENT
1421       double precision scal_el /1.0d0/
1422 #else
1423       double precision scal_el /0.5d0/
1424 #endif
1425 C 12/13/98
1426 C 13-go grudnia roku pamietnego...
1427       double precision unmat(3,3) /1.0d0,0.0d0,0.0d0,
1428      &                   0.0d0,1.0d0,0.0d0,
1429      &                   0.0d0,0.0d0,1.0d0/
1430 c          time00=MPI_Wtime()
1431 cd      write (iout,*) "eelecij",i,j
1432           ind=ind+1
1433           iteli=itel(i)
1434           itelj=itel(j)
1435           if (j.eq.i+2 .and. itelj.eq.2) iteli=2
1436           aaa=app(iteli,itelj)
1437           bbb=bpp(iteli,itelj)
1438           ael6i=ael6(iteli,itelj)
1439           ael3i=ael3(iteli,itelj)
1440           dxj=dc(1,j)
1441           dyj=dc(2,j)
1442           dzj=dc(3,j)
1443           dx_normj=dc_norm(1,j)
1444           dy_normj=dc_norm(2,j)
1445           dz_normj=dc_norm(3,j)
1446           xj=c(1,j)+0.5D0*dxj-xmedi
1447           yj=c(2,j)+0.5D0*dyj-ymedi
1448           zj=c(3,j)+0.5D0*dzj-zmedi
1449           rij=xj*xj+yj*yj+zj*zj
1450           rrmij=1.0D0/rij
1451           rij=dsqrt(rij)
1452           rmij=1.0D0/rij
1453 c For extracting the short-range part of Evdwpp
1454           sss=sscale(rij/rpp(iteli,itelj))
1455
1456           r3ij=rrmij*rmij
1457           r6ij=r3ij*r3ij
1458           cosa=dx_normi*dx_normj+dy_normi*dy_normj+dz_normi*dz_normj
1459           cosb=(xj*dx_normi+yj*dy_normi+zj*dz_normi)*rmij
1460           cosg=(xj*dx_normj+yj*dy_normj+zj*dz_normj)*rmij
1461           fac=cosa-3.0D0*cosb*cosg
1462           ev1=aaa*r6ij*r6ij
1463 c 4/26/02 - AL scaling down 1,4 repulsive VDW interactions
1464           if (j.eq.i+2) ev1=scal_el*ev1
1465           ev2=bbb*r6ij
1466           fac3=ael6i*r6ij
1467           fac4=ael3i*r3ij
1468           evdwij=ev1+ev2
1469           el1=fac3*(4.0D0+fac*fac-3.0D0*(cosb*cosb+cosg*cosg))
1470           el2=fac4*fac
1471           eesij=el1+el2
1472 C 12/26/95 - for the evaluation of multi-body H-bonding interactions
1473           ees0ij=4.0D0+fac*fac-3.0D0*(cosb*cosb+cosg*cosg)
1474           ees=ees+eesij
1475           evdw1=evdw1+evdwij*(1.0d0-sss)
1476 cd          write(iout,'(2(2i3,2x),7(1pd12.4)/2(3(1pd12.4),5x)/)')
1477 cd     &      iteli,i,itelj,j,aaa,bbb,ael6i,ael3i,
1478 cd     &      1.0D0/dsqrt(rrmij),evdwij,eesij,
1479 cd     &      xmedi,ymedi,zmedi,xj,yj,zj
1480
1481           if (energy_dec) then
1482               write (iout,'(a6,2i5,0pf7.3,f7.3)') 'evdw1',i,j,evdwij,sss
1483               write (iout,'(a6,2i5,0pf7.3)') 'ees',i,j,eesij
1484           endif
1485
1486 C
1487 C Calculate contributions to the Cartesian gradient.
1488 C
1489 #ifdef SPLITELE
1490           facvdw=-6*rrmij*(ev1+evdwij)*(1.0d0-sss)
1491           facel=-3*rrmij*(el1+eesij)
1492           fac1=fac
1493           erij(1)=xj*rmij
1494           erij(2)=yj*rmij
1495           erij(3)=zj*rmij
1496 *
1497 * Radial derivatives. First process both termini of the fragment (i,j)
1498 *
1499           ggg(1)=facel*xj
1500           ggg(2)=facel*yj
1501           ggg(3)=facel*zj
1502 c          do k=1,3
1503 c            ghalf=0.5D0*ggg(k)
1504 c            gelc(k,i)=gelc(k,i)+ghalf
1505 c            gelc(k,j)=gelc(k,j)+ghalf
1506 c          enddo
1507 c 9/28/08 AL Gradient compotents will be summed only at the end
1508           do k=1,3
1509             gelc_long(k,j)=gelc_long(k,j)+ggg(k)
1510             gelc_long(k,i)=gelc_long(k,i)-ggg(k)
1511           enddo
1512 *
1513 * Loop over residues i+1 thru j-1.
1514 *
1515 cgrad          do k=i+1,j-1
1516 cgrad            do l=1,3
1517 cgrad              gelc(l,k)=gelc(l,k)+ggg(l)
1518 cgrad            enddo
1519 cgrad          enddo
1520           ggg(1)=facvdw*xj
1521           ggg(2)=facvdw*yj
1522           ggg(3)=facvdw*zj
1523 c          do k=1,3
1524 c            ghalf=0.5D0*ggg(k)
1525 c            gvdwpp(k,i)=gvdwpp(k,i)+ghalf
1526 c            gvdwpp(k,j)=gvdwpp(k,j)+ghalf
1527 c          enddo
1528 c 9/28/08 AL Gradient compotents will be summed only at the end
1529           do k=1,3
1530             gvdwpp(k,j)=gvdwpp(k,j)+ggg(k)
1531             gvdwpp(k,i)=gvdwpp(k,i)-ggg(k)
1532           enddo
1533 *
1534 * Loop over residues i+1 thru j-1.
1535 *
1536 cgrad          do k=i+1,j-1
1537 cgrad            do l=1,3
1538 cgrad              gvdwpp(l,k)=gvdwpp(l,k)+ggg(l)
1539 cgrad            enddo
1540 cgrad          enddo
1541 #else
1542           facvdw=ev1+evdwij*(1.0d0-sss)
1543           facel=el1+eesij
1544           fac1=fac
1545           fac=-3*rrmij*(facvdw+facvdw+facel)
1546           erij(1)=xj*rmij
1547           erij(2)=yj*rmij
1548           erij(3)=zj*rmij
1549 *
1550 * Radial derivatives. First process both termini of the fragment (i,j)
1551 *
1552           ggg(1)=fac*xj
1553           ggg(2)=fac*yj
1554           ggg(3)=fac*zj
1555 c          do k=1,3
1556 c            ghalf=0.5D0*ggg(k)
1557 c            gelc(k,i)=gelc(k,i)+ghalf
1558 c            gelc(k,j)=gelc(k,j)+ghalf
1559 c          enddo
1560 c 9/28/08 AL Gradient compotents will be summed only at the end
1561           do k=1,3
1562             gelc_long(k,j)=gelc(k,j)+ggg(k)
1563             gelc_long(k,i)=gelc(k,i)-ggg(k)
1564           enddo
1565 *
1566 * Loop over residues i+1 thru j-1.
1567 *
1568 cgrad          do k=i+1,j-1
1569 cgrad            do l=1,3
1570 cgrad              gelc(l,k)=gelc(l,k)+ggg(l)
1571 cgrad            enddo
1572 cgrad          enddo
1573 c 9/28/08 AL Gradient compotents will be summed only at the end
1574           ggg(1)=facvdw*xj
1575           ggg(2)=facvdw*yj
1576           ggg(3)=facvdw*zj
1577           do k=1,3
1578             gvdwpp(k,j)=gvdwpp(k,j)+ggg(k)
1579             gvdwpp(k,i)=gvdwpp(k,i)-ggg(k)
1580           enddo
1581 #endif
1582 *
1583 * Angular part
1584 *
1585           ecosa=2.0D0*fac3*fac1+fac4
1586           fac4=-3.0D0*fac4
1587           fac3=-6.0D0*fac3
1588           ecosb=(fac3*(fac1*cosg+cosb)+cosg*fac4)
1589           ecosg=(fac3*(fac1*cosb+cosg)+cosb*fac4)
1590           do k=1,3
1591             dcosb(k)=rmij*(dc_norm(k,i)-erij(k)*cosb)
1592             dcosg(k)=rmij*(dc_norm(k,j)-erij(k)*cosg)
1593           enddo
1594 cd        print '(2i3,2(3(1pd14.5),3x))',i,j,(dcosb(k),k=1,3),
1595 cd   &          (dcosg(k),k=1,3)
1596           do k=1,3
1597             ggg(k)=ecosb*dcosb(k)+ecosg*dcosg(k)
1598           enddo
1599 c          do k=1,3
1600 c            ghalf=0.5D0*ggg(k)
1601 c            gelc(k,i)=gelc(k,i)+ghalf
1602 c     &               +(ecosa*(dc_norm(k,j)-cosa*dc_norm(k,i))
1603 c     &               + ecosb*(erij(k)-cosb*dc_norm(k,i)))*vbld_inv(i+1)
1604 c            gelc(k,j)=gelc(k,j)+ghalf
1605 c     &               +(ecosa*(dc_norm(k,i)-cosa*dc_norm(k,j))
1606 c     &               + ecosg*(erij(k)-cosg*dc_norm(k,j)))*vbld_inv(j+1)
1607 c          enddo
1608 cgrad          do k=i+1,j-1
1609 cgrad            do l=1,3
1610 cgrad              gelc(l,k)=gelc(l,k)+ggg(l)
1611 cgrad            enddo
1612 cgrad          enddo
1613           do k=1,3
1614             gelc(k,i)=gelc(k,i)
1615      &               +(ecosa*(dc_norm(k,j)-cosa*dc_norm(k,i))
1616      &               + ecosb*(erij(k)-cosb*dc_norm(k,i)))*vbld_inv(i+1)
1617             gelc(k,j)=gelc(k,j)
1618      &               +(ecosa*(dc_norm(k,i)-cosa*dc_norm(k,j))
1619      &               + ecosg*(erij(k)-cosg*dc_norm(k,j)))*vbld_inv(j+1)
1620             gelc_long(k,j)=gelc_long(k,j)+ggg(k)
1621             gelc_long(k,i)=gelc_long(k,i)-ggg(k)
1622           enddo
1623           IF (wel_loc.gt.0.0d0 .or. wcorr4.gt.0.0d0 .or. wcorr5.gt.0.0d0
1624      &        .or. wcorr6.gt.0.0d0 .or. wturn3.gt.0.0d0
1625      &        .or. wturn4.gt.0.0d0 .or. wturn6.gt.0.0d0) THEN
1626 C
1627 C 9/25/99 Mixed third-order local-electrostatic terms. The local-interaction
1628 C   energy of a peptide unit is assumed in the form of a second-order
1629 C   Fourier series in the angles lambda1 and lambda2 (see Nishikawa et al.
1630 C   Macromolecules, 1974, 7, 797-806 for definition). This correlation terms
1631 C   are computed for EVERY pair of non-contiguous peptide groups.
1632 C
1633           if (j.lt.nres-1) then
1634             j1=j+1
1635             j2=j-1
1636           else
1637             j1=j-1
1638             j2=j-2
1639           endif
1640           kkk=0
1641           do k=1,2
1642             do l=1,2
1643               kkk=kkk+1
1644               muij(kkk)=mu(k,i)*mu(l,j)
1645             enddo
1646           enddo
1647 cd         write (iout,*) 'EELEC: i',i,' j',j
1648 cd          write (iout,*) 'j',j,' j1',j1,' j2',j2
1649 cd          write(iout,*) 'muij',muij
1650           ury=scalar(uy(1,i),erij)
1651           urz=scalar(uz(1,i),erij)
1652           vry=scalar(uy(1,j),erij)
1653           vrz=scalar(uz(1,j),erij)
1654           a22=scalar(uy(1,i),uy(1,j))-3*ury*vry
1655           a23=scalar(uy(1,i),uz(1,j))-3*ury*vrz
1656           a32=scalar(uz(1,i),uy(1,j))-3*urz*vry
1657           a33=scalar(uz(1,i),uz(1,j))-3*urz*vrz
1658           fac=dsqrt(-ael6i)*r3ij
1659           a22=a22*fac
1660           a23=a23*fac
1661           a32=a32*fac
1662           a33=a33*fac
1663 cd          write (iout,'(4i5,4f10.5)')
1664 cd     &     i,itortyp(itype(i)),j,itortyp(itype(j)),a22,a23,a32,a33
1665 cd          write (iout,'(6f10.5)') (muij(k),k=1,4),fac,eel_loc_ij
1666 cd          write (iout,'(2(3f10.5,5x)/2(3f10.5,5x))') uy(:,i),uz(:,i),
1667 cd     &      uy(:,j),uz(:,j)
1668 cd          write (iout,'(4f10.5)')
1669 cd     &      scalar(uy(1,i),uy(1,j)),scalar(uy(1,i),uz(1,j)),
1670 cd     &      scalar(uz(1,i),uy(1,j)),scalar(uz(1,i),uz(1,j))
1671 cd          write (iout,'(4f10.5)') ury,urz,vry,vrz
1672 cd           write (iout,'(9f10.5/)')
1673 cd     &      fac22,a22,fac23,a23,fac32,a32,fac33,a33,eel_loc_ij
1674 C Derivatives of the elements of A in virtual-bond vectors
1675           call unormderiv(erij(1),unmat(1,1),rmij,erder(1,1))
1676           do k=1,3
1677             uryg(k,1)=scalar(erder(1,k),uy(1,i))
1678             uryg(k,2)=scalar(uygrad(1,k,1,i),erij(1))
1679             uryg(k,3)=scalar(uygrad(1,k,2,i),erij(1))
1680             urzg(k,1)=scalar(erder(1,k),uz(1,i))
1681             urzg(k,2)=scalar(uzgrad(1,k,1,i),erij(1))
1682             urzg(k,3)=scalar(uzgrad(1,k,2,i),erij(1))
1683             vryg(k,1)=scalar(erder(1,k),uy(1,j))
1684             vryg(k,2)=scalar(uygrad(1,k,1,j),erij(1))
1685             vryg(k,3)=scalar(uygrad(1,k,2,j),erij(1))
1686             vrzg(k,1)=scalar(erder(1,k),uz(1,j))
1687             vrzg(k,2)=scalar(uzgrad(1,k,1,j),erij(1))
1688             vrzg(k,3)=scalar(uzgrad(1,k,2,j),erij(1))
1689           enddo
1690 C Compute radial contributions to the gradient
1691           facr=-3.0d0*rrmij
1692           a22der=a22*facr
1693           a23der=a23*facr
1694           a32der=a32*facr
1695           a33der=a33*facr
1696           agg(1,1)=a22der*xj
1697           agg(2,1)=a22der*yj
1698           agg(3,1)=a22der*zj
1699           agg(1,2)=a23der*xj
1700           agg(2,2)=a23der*yj
1701           agg(3,2)=a23der*zj
1702           agg(1,3)=a32der*xj
1703           agg(2,3)=a32der*yj
1704           agg(3,3)=a32der*zj
1705           agg(1,4)=a33der*xj
1706           agg(2,4)=a33der*yj
1707           agg(3,4)=a33der*zj
1708 C Add the contributions coming from er
1709           fac3=-3.0d0*fac
1710           do k=1,3
1711             agg(k,1)=agg(k,1)+fac3*(uryg(k,1)*vry+vryg(k,1)*ury)
1712             agg(k,2)=agg(k,2)+fac3*(uryg(k,1)*vrz+vrzg(k,1)*ury)
1713             agg(k,3)=agg(k,3)+fac3*(urzg(k,1)*vry+vryg(k,1)*urz)
1714             agg(k,4)=agg(k,4)+fac3*(urzg(k,1)*vrz+vrzg(k,1)*urz)
1715           enddo
1716           do k=1,3
1717 C Derivatives in DC(i)
1718 cgrad            ghalf1=0.5d0*agg(k,1)
1719 cgrad            ghalf2=0.5d0*agg(k,2)
1720 cgrad            ghalf3=0.5d0*agg(k,3)
1721 cgrad            ghalf4=0.5d0*agg(k,4)
1722             aggi(k,1)=fac*(scalar(uygrad(1,k,1,i),uy(1,j))
1723      &      -3.0d0*uryg(k,2)*vry)!+ghalf1
1724             aggi(k,2)=fac*(scalar(uygrad(1,k,1,i),uz(1,j))
1725      &      -3.0d0*uryg(k,2)*vrz)!+ghalf2
1726             aggi(k,3)=fac*(scalar(uzgrad(1,k,1,i),uy(1,j))
1727      &      -3.0d0*urzg(k,2)*vry)!+ghalf3
1728             aggi(k,4)=fac*(scalar(uzgrad(1,k,1,i),uz(1,j))
1729      &      -3.0d0*urzg(k,2)*vrz)!+ghalf4
1730 C Derivatives in DC(i+1)
1731             aggi1(k,1)=fac*(scalar(uygrad(1,k,2,i),uy(1,j))
1732      &      -3.0d0*uryg(k,3)*vry)!+agg(k,1)
1733             aggi1(k,2)=fac*(scalar(uygrad(1,k,2,i),uz(1,j))
1734      &      -3.0d0*uryg(k,3)*vrz)!+agg(k,2)
1735             aggi1(k,3)=fac*(scalar(uzgrad(1,k,2,i),uy(1,j))
1736      &      -3.0d0*urzg(k,3)*vry)!+agg(k,3)
1737             aggi1(k,4)=fac*(scalar(uzgrad(1,k,2,i),uz(1,j))
1738      &      -3.0d0*urzg(k,3)*vrz)!+agg(k,4)
1739 C Derivatives in DC(j)
1740             aggj(k,1)=fac*(scalar(uygrad(1,k,1,j),uy(1,i))
1741      &      -3.0d0*vryg(k,2)*ury)!+ghalf1
1742             aggj(k,2)=fac*(scalar(uzgrad(1,k,1,j),uy(1,i))
1743      &      -3.0d0*vrzg(k,2)*ury)!+ghalf2
1744             aggj(k,3)=fac*(scalar(uygrad(1,k,1,j),uz(1,i))
1745      &      -3.0d0*vryg(k,2)*urz)!+ghalf3
1746             aggj(k,4)=fac*(scalar(uzgrad(1,k,1,j),uz(1,i))
1747      &      -3.0d0*vrzg(k,2)*urz)!+ghalf4
1748 C Derivatives in DC(j+1) or DC(nres-1)
1749             aggj1(k,1)=fac*(scalar(uygrad(1,k,2,j),uy(1,i))
1750      &      -3.0d0*vryg(k,3)*ury)
1751             aggj1(k,2)=fac*(scalar(uzgrad(1,k,2,j),uy(1,i))
1752      &      -3.0d0*vrzg(k,3)*ury)
1753             aggj1(k,3)=fac*(scalar(uygrad(1,k,2,j),uz(1,i))
1754      &      -3.0d0*vryg(k,3)*urz)
1755             aggj1(k,4)=fac*(scalar(uzgrad(1,k,2,j),uz(1,i))
1756      &      -3.0d0*vrzg(k,3)*urz)
1757 cgrad            if (j.eq.nres-1 .and. i.lt.j-2) then
1758 cgrad              do l=1,4
1759 cgrad                aggj1(k,l)=aggj1(k,l)+agg(k,l)
1760 cgrad              enddo
1761 cgrad            endif
1762           enddo
1763           acipa(1,1)=a22
1764           acipa(1,2)=a23
1765           acipa(2,1)=a32
1766           acipa(2,2)=a33
1767           a22=-a22
1768           a23=-a23
1769           do l=1,2
1770             do k=1,3
1771               agg(k,l)=-agg(k,l)
1772               aggi(k,l)=-aggi(k,l)
1773               aggi1(k,l)=-aggi1(k,l)
1774               aggj(k,l)=-aggj(k,l)
1775               aggj1(k,l)=-aggj1(k,l)
1776             enddo
1777           enddo
1778           if (j.lt.nres-1) then
1779             a22=-a22
1780             a32=-a32
1781             do l=1,3,2
1782               do k=1,3
1783                 agg(k,l)=-agg(k,l)
1784                 aggi(k,l)=-aggi(k,l)
1785                 aggi1(k,l)=-aggi1(k,l)
1786                 aggj(k,l)=-aggj(k,l)
1787                 aggj1(k,l)=-aggj1(k,l)
1788               enddo
1789             enddo
1790           else
1791             a22=-a22
1792             a23=-a23
1793             a32=-a32
1794             a33=-a33
1795             do l=1,4
1796               do k=1,3
1797                 agg(k,l)=-agg(k,l)
1798                 aggi(k,l)=-aggi(k,l)
1799                 aggi1(k,l)=-aggi1(k,l)
1800                 aggj(k,l)=-aggj(k,l)
1801                 aggj1(k,l)=-aggj1(k,l)
1802               enddo
1803             enddo
1804           endif
1805           ENDIF ! WCORR
1806           IF (wel_loc.gt.0.0d0) THEN
1807 C Contribution to the local-electrostatic energy coming from the i-j pair
1808           eel_loc_ij=a22*muij(1)+a23*muij(2)+a32*muij(3)
1809      &     +a33*muij(4)
1810 cd          write (iout,*) 'i',i,' j',j,' eel_loc_ij',eel_loc_ij
1811
1812           if (energy_dec) write (iout,'(a6,2i5,0pf7.3)')
1813      &            'eelloc',i,j,eel_loc_ij
1814
1815           eel_loc=eel_loc+eel_loc_ij
1816 C Partial derivatives in virtual-bond dihedral angles gamma
1817           if (i.gt.1)
1818      &    gel_loc_loc(i-1)=gel_loc_loc(i-1)+
1819      &            a22*muder(1,i)*mu(1,j)+a23*muder(1,i)*mu(2,j)
1820      &           +a32*muder(2,i)*mu(1,j)+a33*muder(2,i)*mu(2,j)
1821           gel_loc_loc(j-1)=gel_loc_loc(j-1)+
1822      &            a22*mu(1,i)*muder(1,j)+a23*mu(1,i)*muder(2,j)
1823      &           +a32*mu(2,i)*muder(1,j)+a33*mu(2,i)*muder(2,j)
1824 C Derivatives of eello in DC(i+1) thru DC(j-1) or DC(nres-2)
1825           do l=1,3
1826             ggg(l)=agg(l,1)*muij(1)+
1827      &          agg(l,2)*muij(2)+agg(l,3)*muij(3)+agg(l,4)*muij(4)
1828             gel_loc_long(l,j)=gel_loc_long(l,j)+ggg(l)
1829             gel_loc_long(l,i)=gel_loc_long(l,i)-ggg(l)
1830 cgrad            ghalf=0.5d0*ggg(l)
1831 cgrad            gel_loc(l,i)=gel_loc(l,i)+ghalf
1832 cgrad            gel_loc(l,j)=gel_loc(l,j)+ghalf
1833           enddo
1834 cgrad          do k=i+1,j2
1835 cgrad            do l=1,3
1836 cgrad              gel_loc(l,k)=gel_loc(l,k)+ggg(l)
1837 cgrad            enddo
1838 cgrad          enddo
1839 C Remaining derivatives of eello
1840           do l=1,3
1841             gel_loc(l,i)=gel_loc(l,i)+aggi(l,1)*muij(1)+
1842      &          aggi(l,2)*muij(2)+aggi(l,3)*muij(3)+aggi(l,4)*muij(4)
1843             gel_loc(l,i+1)=gel_loc(l,i+1)+aggi1(l,1)*muij(1)+
1844      &          aggi1(l,2)*muij(2)+aggi1(l,3)*muij(3)+aggi1(l,4)*muij(4)
1845             gel_loc(l,j)=gel_loc(l,j)+aggj(l,1)*muij(1)+
1846      &          aggj(l,2)*muij(2)+aggj(l,3)*muij(3)+aggj(l,4)*muij(4)
1847             gel_loc(l,j1)=gel_loc(l,j1)+aggj1(l,1)*muij(1)+
1848      &          aggj1(l,2)*muij(2)+aggj1(l,3)*muij(3)+aggj1(l,4)*muij(4)
1849           enddo
1850           ENDIF
1851 C Change 12/26/95 to calculate four-body contributions to H-bonding energy
1852 c          if (j.gt.i+1 .and. num_conti.le.maxconts) then
1853           if (wcorr+wcorr4+wcorr5+wcorr6.gt.0.0d0
1854      &       .and. num_conti.le.maxconts) then
1855 c            write (iout,*) i,j," entered corr"
1856 C
1857 C Calculate the contact function. The ith column of the array JCONT will
1858 C contain the numbers of atoms that make contacts with the atom I (of numbers
1859 C greater than I). The arrays FACONT and GACONT will contain the values of
1860 C the contact function and its derivative.
1861 c           r0ij=1.02D0*rpp(iteli,itelj)
1862 c           r0ij=1.11D0*rpp(iteli,itelj)
1863             r0ij=2.20D0*rpp(iteli,itelj)
1864 c           r0ij=1.55D0*rpp(iteli,itelj)
1865             call gcont(rij,r0ij,1.0D0,0.2d0*r0ij,fcont,fprimcont)
1866             if (fcont.gt.0.0D0) then
1867               num_conti=num_conti+1
1868               if (num_conti.gt.maxconts) then
1869                 write (iout,*) 'WARNING - max. # of contacts exceeded;',
1870      &                         ' will skip next contacts for this conf.'
1871               else
1872                 jcont_hb(num_conti,i)=j
1873 cd                write (iout,*) "i",i," j",j," num_conti",num_conti,
1874 cd     &           " jcont_hb",jcont_hb(num_conti,i)
1875                 IF (wcorr4.gt.0.0d0 .or. wcorr5.gt.0.0d0 .or.
1876      &          wcorr6.gt.0.0d0 .or. wturn6.gt.0.0d0) THEN
1877 C 9/30/99 (AL) - store components necessary to evaluate higher-order loc-el
1878 C  terms.
1879                 d_cont(num_conti,i)=rij
1880 cd                write (2,'(3e15.5)') rij,r0ij+0.2d0*r0ij,rij
1881 C     --- Electrostatic-interaction matrix ---
1882                 a_chuj(1,1,num_conti,i)=a22
1883                 a_chuj(1,2,num_conti,i)=a23
1884                 a_chuj(2,1,num_conti,i)=a32
1885                 a_chuj(2,2,num_conti,i)=a33
1886 C     --- Gradient of rij
1887                 do kkk=1,3
1888                   grij_hb_cont(kkk,num_conti,i)=erij(kkk)
1889                 enddo
1890                 kkll=0
1891                 do k=1,2
1892                   do l=1,2
1893                     kkll=kkll+1
1894                     do m=1,3
1895                       a_chuj_der(k,l,m,1,num_conti,i)=agg(m,kkll)
1896                       a_chuj_der(k,l,m,2,num_conti,i)=aggi(m,kkll)
1897                       a_chuj_der(k,l,m,3,num_conti,i)=aggi1(m,kkll)
1898                       a_chuj_der(k,l,m,4,num_conti,i)=aggj(m,kkll)
1899                       a_chuj_der(k,l,m,5,num_conti,i)=aggj1(m,kkll)
1900                     enddo
1901                   enddo
1902                 enddo
1903                 ENDIF
1904                 IF (wcorr4.eq.0.0d0 .and. wcorr.gt.0.0d0) THEN
1905 C Calculate contact energies
1906                 cosa4=4.0D0*cosa
1907                 wij=cosa-3.0D0*cosb*cosg
1908                 cosbg1=cosb+cosg
1909                 cosbg2=cosb-cosg
1910 c               fac3=dsqrt(-ael6i)/r0ij**3
1911                 fac3=dsqrt(-ael6i)*r3ij
1912 c                 ees0pij=dsqrt(4.0D0+cosa4+wij*wij-3.0D0*cosbg1*cosbg1)
1913                 ees0tmp=4.0D0+cosa4+wij*wij-3.0D0*cosbg1*cosbg1
1914                 if (ees0tmp.gt.0) then
1915                   ees0pij=dsqrt(ees0tmp)
1916                 else
1917                   ees0pij=0
1918                 endif
1919 c                ees0mij=dsqrt(4.0D0-cosa4+wij*wij-3.0D0*cosbg2*cosbg2)
1920                 ees0tmp=4.0D0-cosa4+wij*wij-3.0D0*cosbg2*cosbg2
1921                 if (ees0tmp.gt.0) then
1922                   ees0mij=dsqrt(ees0tmp)
1923                 else
1924                   ees0mij=0
1925                 endif
1926 c               ees0mij=0.0D0
1927                 ees0p(num_conti,i)=0.5D0*fac3*(ees0pij+ees0mij)
1928                 ees0m(num_conti,i)=0.5D0*fac3*(ees0pij-ees0mij)
1929 C Diagnostics. Comment out or remove after debugging!
1930 c               ees0p(num_conti,i)=0.5D0*fac3*ees0pij
1931 c               ees0m(num_conti,i)=0.5D0*fac3*ees0mij
1932 c               ees0m(num_conti,i)=0.0D0
1933 C End diagnostics.
1934 c               write (iout,*) 'i=',i,' j=',j,' rij=',rij,' r0ij=',r0ij,
1935 c    & ' ees0ij=',ees0p(num_conti,i),ees0m(num_conti,i),' fcont=',fcont
1936 C Angular derivatives of the contact function
1937                 ees0pij1=fac3/ees0pij
1938                 ees0mij1=fac3/ees0mij
1939                 fac3p=-3.0D0*fac3*rrmij
1940                 ees0pijp=0.5D0*fac3p*(ees0pij+ees0mij)
1941                 ees0mijp=0.5D0*fac3p*(ees0pij-ees0mij)
1942 c               ees0mij1=0.0D0
1943                 ecosa1=       ees0pij1*( 1.0D0+0.5D0*wij)
1944                 ecosb1=-1.5D0*ees0pij1*(wij*cosg+cosbg1)
1945                 ecosg1=-1.5D0*ees0pij1*(wij*cosb+cosbg1)
1946                 ecosa2=       ees0mij1*(-1.0D0+0.5D0*wij)
1947                 ecosb2=-1.5D0*ees0mij1*(wij*cosg+cosbg2)
1948                 ecosg2=-1.5D0*ees0mij1*(wij*cosb-cosbg2)
1949                 ecosap=ecosa1+ecosa2
1950                 ecosbp=ecosb1+ecosb2
1951                 ecosgp=ecosg1+ecosg2
1952                 ecosam=ecosa1-ecosa2
1953                 ecosbm=ecosb1-ecosb2
1954                 ecosgm=ecosg1-ecosg2
1955 C Diagnostics
1956 c               ecosap=ecosa1
1957 c               ecosbp=ecosb1
1958 c               ecosgp=ecosg1
1959 c               ecosam=0.0D0
1960 c               ecosbm=0.0D0
1961 c               ecosgm=0.0D0
1962 C End diagnostics
1963                 facont_hb(num_conti,i)=fcont
1964                 fprimcont=fprimcont/rij
1965 cd              facont_hb(num_conti,i)=1.0D0
1966 C Following line is for diagnostics.
1967 cd              fprimcont=0.0D0
1968                 do k=1,3
1969                   dcosb(k)=rmij*(dc_norm(k,i)-erij(k)*cosb)
1970                   dcosg(k)=rmij*(dc_norm(k,j)-erij(k)*cosg)
1971                 enddo
1972                 do k=1,3
1973                   gggp(k)=ecosbp*dcosb(k)+ecosgp*dcosg(k)
1974                   gggm(k)=ecosbm*dcosb(k)+ecosgm*dcosg(k)
1975                 enddo
1976                 gggp(1)=gggp(1)+ees0pijp*xj
1977                 gggp(2)=gggp(2)+ees0pijp*yj
1978                 gggp(3)=gggp(3)+ees0pijp*zj
1979                 gggm(1)=gggm(1)+ees0mijp*xj
1980                 gggm(2)=gggm(2)+ees0mijp*yj
1981                 gggm(3)=gggm(3)+ees0mijp*zj
1982 C Derivatives due to the contact function
1983                 gacont_hbr(1,num_conti,i)=fprimcont*xj
1984                 gacont_hbr(2,num_conti,i)=fprimcont*yj
1985                 gacont_hbr(3,num_conti,i)=fprimcont*zj
1986                 do k=1,3
1987 c
1988 c 10/24/08 cgrad and ! comments indicate the parts of the code removed
1989 c          following the change of gradient-summation algorithm.
1990 c
1991 cgrad                  ghalfp=0.5D0*gggp(k)
1992 cgrad                  ghalfm=0.5D0*gggm(k)
1993                   gacontp_hb1(k,num_conti,i)=!ghalfp
1994      &              +(ecosap*(dc_norm(k,j)-cosa*dc_norm(k,i))
1995      &              + ecosbp*(erij(k)-cosb*dc_norm(k,i)))*vbld_inv(i+1)
1996                   gacontp_hb2(k,num_conti,i)=!ghalfp
1997      &              +(ecosap*(dc_norm(k,i)-cosa*dc_norm(k,j))
1998      &              + ecosgp*(erij(k)-cosg*dc_norm(k,j)))*vbld_inv(j+1)
1999                   gacontp_hb3(k,num_conti,i)=gggp(k)
2000                   gacontm_hb1(k,num_conti,i)=!ghalfm
2001      &              +(ecosam*(dc_norm(k,j)-cosa*dc_norm(k,i))
2002      &              + ecosbm*(erij(k)-cosb*dc_norm(k,i)))*vbld_inv(i+1)
2003                   gacontm_hb2(k,num_conti,i)=!ghalfm
2004      &              +(ecosam*(dc_norm(k,i)-cosa*dc_norm(k,j))
2005      &              + ecosgm*(erij(k)-cosg*dc_norm(k,j)))*vbld_inv(j+1)
2006                   gacontm_hb3(k,num_conti,i)=gggm(k)
2007                 enddo
2008               ENDIF ! wcorr
2009               endif  ! num_conti.le.maxconts
2010             endif  ! fcont.gt.0
2011           endif    ! j.gt.i+1
2012           if (wturn3.gt.0.0d0 .or. wturn4.gt.0.0d0) then
2013             do k=1,4
2014               do l=1,3
2015                 ghalf=0.5d0*agg(l,k)
2016                 aggi(l,k)=aggi(l,k)+ghalf
2017                 aggi1(l,k)=aggi1(l,k)+agg(l,k)
2018                 aggj(l,k)=aggj(l,k)+ghalf
2019               enddo
2020             enddo
2021             if (j.eq.nres-1 .and. i.lt.j-2) then
2022               do k=1,4
2023                 do l=1,3
2024                   aggj1(l,k)=aggj1(l,k)+agg(l,k)
2025                 enddo
2026               enddo
2027             endif
2028           endif
2029 c          t_eelecij=t_eelecij+MPI_Wtime()-time00
2030       return
2031       end
2032 C-----------------------------------------------------------------------
2033       subroutine evdwpp_short(evdw1)
2034 C
2035 C Compute Evdwpp
2036 C
2037       implicit real*8 (a-h,o-z)
2038       include 'DIMENSIONS'
2039       include 'COMMON.CONTROL'
2040       include 'COMMON.IOUNITS'
2041       include 'COMMON.GEO'
2042       include 'COMMON.VAR'
2043       include 'COMMON.LOCAL'
2044       include 'COMMON.CHAIN'
2045       include 'COMMON.DERIV'
2046       include 'COMMON.INTERACT'
2047       include 'COMMON.CONTACTS'
2048       include 'COMMON.TORSION'
2049       include 'COMMON.VECTORS'
2050       include 'COMMON.FFIELD'
2051       dimension ggg(3)
2052 c 4/26/02 - AL scaling factor for 1,4 repulsive VDW interactions
2053 #ifdef MOMENT
2054       double precision scal_el /1.0d0/
2055 #else
2056       double precision scal_el /0.5d0/
2057 #endif
2058       evdw1=0.0D0
2059 c      write (iout,*) "iatel_s_vdw",iatel_s_vdw,
2060 c     & " iatel_e_vdw",iatel_e_vdw
2061       call flush(iout)
2062       do i=iatel_s_vdw,iatel_e_vdw
2063         dxi=dc(1,i)
2064         dyi=dc(2,i)
2065         dzi=dc(3,i)
2066         dx_normi=dc_norm(1,i)
2067         dy_normi=dc_norm(2,i)
2068         dz_normi=dc_norm(3,i)
2069         xmedi=c(1,i)+0.5d0*dxi
2070         ymedi=c(2,i)+0.5d0*dyi
2071         zmedi=c(3,i)+0.5d0*dzi
2072         num_conti=0
2073 c        write (iout,*) 'i',i,' ielstart',ielstart_vdw(i),
2074 c     &   ' ielend',ielend_vdw(i)
2075         call flush(iout)
2076         do j=ielstart_vdw(i),ielend_vdw(i)
2077           ind=ind+1
2078           iteli=itel(i)
2079           itelj=itel(j)
2080           if (j.eq.i+2 .and. itelj.eq.2) iteli=2
2081           aaa=app(iteli,itelj)
2082           bbb=bpp(iteli,itelj)
2083           dxj=dc(1,j)
2084           dyj=dc(2,j)
2085           dzj=dc(3,j)
2086           dx_normj=dc_norm(1,j)
2087           dy_normj=dc_norm(2,j)
2088           dz_normj=dc_norm(3,j)
2089           xj=c(1,j)+0.5D0*dxj-xmedi
2090           yj=c(2,j)+0.5D0*dyj-ymedi
2091           zj=c(3,j)+0.5D0*dzj-zmedi
2092           rij=xj*xj+yj*yj+zj*zj
2093           rrmij=1.0D0/rij
2094           rij=dsqrt(rij)
2095           sss=sscale(rij/rpp(iteli,itelj))
2096           if (sss.gt.0.0d0) then
2097             rmij=1.0D0/rij
2098             r3ij=rrmij*rmij
2099             r6ij=r3ij*r3ij
2100             ev1=aaa*r6ij*r6ij
2101 c 4/26/02 - AL scaling down 1,4 repulsive VDW interactions
2102             if (j.eq.i+2) ev1=scal_el*ev1
2103             ev2=bbb*r6ij
2104             evdwij=ev1+ev2
2105             if (energy_dec) then
2106               write (iout,'(a6,2i5,0pf7.3,f7.3)') 'evdw1',i,j,evdwij,sss
2107             endif
2108             evdw1=evdw1+evdwij*sss
2109 C
2110 C Calculate contributions to the Cartesian gradient.
2111 C
2112             facvdw=-6*rrmij*(ev1+evdwij)*sss
2113             ggg(1)=facvdw*xj
2114             ggg(2)=facvdw*yj
2115             ggg(3)=facvdw*zj
2116             do k=1,3
2117               gvdwpp(k,j)=gvdwpp(k,j)+ggg(k)
2118               gvdwpp(k,i)=gvdwpp(k,i)-ggg(k)
2119             enddo
2120           endif
2121         enddo ! j
2122       enddo   ! i
2123       return
2124       end
2125 C-----------------------------------------------------------------------------
2126       subroutine escp_long(evdw2,evdw2_14)
2127 C
2128 C This subroutine calculates the excluded-volume interaction energy between
2129 C peptide-group centers and side chains and its gradient in virtual-bond and
2130 C side-chain vectors.
2131 C
2132       implicit real*8 (a-h,o-z)
2133       include 'DIMENSIONS'
2134       include 'COMMON.GEO'
2135       include 'COMMON.VAR'
2136       include 'COMMON.LOCAL'
2137       include 'COMMON.CHAIN'
2138       include 'COMMON.DERIV'
2139       include 'COMMON.INTERACT'
2140       include 'COMMON.FFIELD'
2141       include 'COMMON.IOUNITS'
2142       include 'COMMON.CONTROL'
2143       dimension ggg(3)
2144       evdw2=0.0D0
2145       evdw2_14=0.0d0
2146 cd    print '(a)','Enter ESCP'
2147 cd    write (iout,*) 'iatscp_s=',iatscp_s,' iatscp_e=',iatscp_e
2148       do i=iatscp_s,iatscp_e
2149         iteli=itel(i)
2150         xi=0.5D0*(c(1,i)+c(1,i+1))
2151         yi=0.5D0*(c(2,i)+c(2,i+1))
2152         zi=0.5D0*(c(3,i)+c(3,i+1))
2153
2154         do iint=1,nscp_gr(i)
2155
2156         do j=iscpstart(i,iint),iscpend(i,iint)
2157           itypj=itype(j)
2158 C Uncomment following three lines for SC-p interactions
2159 c         xj=c(1,nres+j)-xi
2160 c         yj=c(2,nres+j)-yi
2161 c         zj=c(3,nres+j)-zi
2162 C Uncomment following three lines for Ca-p interactions
2163           xj=c(1,j)-xi
2164           yj=c(2,j)-yi
2165           zj=c(3,j)-zi
2166           rrij=1.0D0/(xj*xj+yj*yj+zj*zj)
2167
2168           sss=sscale(1.0d0/(dsqrt(rrij)*rscp(itypj,iteli)))
2169
2170           if (sss.lt.1.0d0) then
2171
2172             fac=rrij**expon2
2173             e1=fac*fac*aad(itypj,iteli)
2174             e2=fac*bad(itypj,iteli)
2175             if (iabs(j-i) .le. 2) then
2176               e1=scal14*e1
2177               e2=scal14*e2
2178               evdw2_14=evdw2_14+(e1+e2)*(1.0d0-sss)
2179             endif
2180             evdwij=e1+e2
2181             evdw2=evdw2+evdwij*(1.0d0-sss)
2182             if (energy_dec) write (iout,'(a6,2i5,0pf7.3)')
2183      &          'evdw2',i,j,evdwij
2184 C
2185 C Calculate contributions to the gradient in the virtual-bond and SC vectors.
2186 C
2187             fac=-(evdwij+e1)*rrij*(1.0d0-sss)
2188             ggg(1)=xj*fac
2189             ggg(2)=yj*fac
2190             ggg(3)=zj*fac
2191 C Uncomment following three lines for SC-p interactions
2192 c           do k=1,3
2193 c             gradx_scp(k,j)=gradx_scp(k,j)+ggg(k)
2194 c           enddo
2195 C Uncomment following line for SC-p interactions
2196 c             gradx_scp(k,j)=gradx_scp(k,j)+ggg(k)
2197             do k=1,3
2198               gvdwc_scpp(k,i)=gvdwc_scpp(k,i)-ggg(k)
2199               gvdwc_scp(k,j)=gvdwc_scp(k,j)+ggg(k)
2200             enddo
2201           endif
2202         enddo
2203
2204         enddo ! iint
2205       enddo ! i
2206       do i=1,nct
2207         do j=1,3
2208           gvdwc_scp(j,i)=expon*gvdwc_scp(j,i)
2209           gvdwc_scpp(j,i)=expon*gvdwc_scpp(j,i)
2210           gradx_scp(j,i)=expon*gradx_scp(j,i)
2211         enddo
2212       enddo
2213 C******************************************************************************
2214 C
2215 C                              N O T E !!!
2216 C
2217 C To save time the factor EXPON has been extracted from ALL components
2218 C of GVDWC and GRADX. Remember to multiply them by this factor before further
2219 C use!
2220 C
2221 C******************************************************************************
2222       return
2223       end
2224 C-----------------------------------------------------------------------------
2225       subroutine escp_short(evdw2,evdw2_14)
2226 C
2227 C This subroutine calculates the excluded-volume interaction energy between
2228 C peptide-group centers and side chains and its gradient in virtual-bond and
2229 C side-chain vectors.
2230 C
2231       implicit real*8 (a-h,o-z)
2232       include 'DIMENSIONS'
2233       include 'COMMON.GEO'
2234       include 'COMMON.VAR'
2235       include 'COMMON.LOCAL'
2236       include 'COMMON.CHAIN'
2237       include 'COMMON.DERIV'
2238       include 'COMMON.INTERACT'
2239       include 'COMMON.FFIELD'
2240       include 'COMMON.IOUNITS'
2241       include 'COMMON.CONTROL'
2242       dimension ggg(3)
2243       evdw2=0.0D0
2244       evdw2_14=0.0d0
2245 cd    print '(a)','Enter ESCP'
2246 cd    write (iout,*) 'iatscp_s=',iatscp_s,' iatscp_e=',iatscp_e
2247       do i=iatscp_s,iatscp_e
2248         iteli=itel(i)
2249         xi=0.5D0*(c(1,i)+c(1,i+1))
2250         yi=0.5D0*(c(2,i)+c(2,i+1))
2251         zi=0.5D0*(c(3,i)+c(3,i+1))
2252
2253         do iint=1,nscp_gr(i)
2254
2255         do j=iscpstart(i,iint),iscpend(i,iint)
2256           itypj=itype(j)
2257 C Uncomment following three lines for SC-p interactions
2258 c         xj=c(1,nres+j)-xi
2259 c         yj=c(2,nres+j)-yi
2260 c         zj=c(3,nres+j)-zi
2261 C Uncomment following three lines for Ca-p interactions
2262           xj=c(1,j)-xi
2263           yj=c(2,j)-yi
2264           zj=c(3,j)-zi
2265           rrij=1.0D0/(xj*xj+yj*yj+zj*zj)
2266
2267           sss=sscale(1.0d0/(dsqrt(rrij)*rscp(itypj,iteli)))
2268
2269           if (sss.gt.0.0d0) then
2270
2271             fac=rrij**expon2
2272             e1=fac*fac*aad(itypj,iteli)
2273             e2=fac*bad(itypj,iteli)
2274             if (iabs(j-i) .le. 2) then
2275               e1=scal14*e1
2276               e2=scal14*e2
2277               evdw2_14=evdw2_14+(e1+e2)*sss
2278             endif
2279             evdwij=e1+e2
2280             evdw2=evdw2+evdwij*sss
2281             if (energy_dec) write (iout,'(a6,2i5,0pf7.3)')
2282      &          'evdw2',i,j,evdwij
2283 C
2284 C Calculate contributions to the gradient in the virtual-bond and SC vectors.
2285 C
2286             fac=-(evdwij+e1)*rrij*sss
2287             ggg(1)=xj*fac
2288             ggg(2)=yj*fac
2289             ggg(3)=zj*fac
2290 C Uncomment following three lines for SC-p interactions
2291 c           do k=1,3
2292 c             gradx_scp(k,j)=gradx_scp(k,j)+ggg(k)
2293 c           enddo
2294 C Uncomment following line for SC-p interactions
2295 c             gradx_scp(k,j)=gradx_scp(k,j)+ggg(k)
2296             do k=1,3
2297               gvdwc_scpp(k,i)=gvdwc_scpp(k,i)-ggg(k)
2298               gvdwc_scp(k,j)=gvdwc_scp(k,j)+ggg(k)
2299             enddo
2300           endif
2301         enddo
2302
2303         enddo ! iint
2304       enddo ! i
2305       do i=1,nct
2306         do j=1,3
2307           gvdwc_scp(j,i)=expon*gvdwc_scp(j,i)
2308           gvdwc_scpp(j,i)=expon*gvdwc_scpp(j,i)
2309           gradx_scp(j,i)=expon*gradx_scp(j,i)
2310         enddo
2311       enddo
2312 C******************************************************************************
2313 C
2314 C                              N O T E !!!
2315 C
2316 C To save time the factor EXPON has been extracted from ALL components
2317 C of GVDWC and GRADX. Remember to multiply them by this factor before further
2318 C use!
2319 C
2320 C******************************************************************************
2321       return
2322       end