source/unres/src_MD/src/old_F/energy_p_new-sep.F.org

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