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

   1 C-----------------------------------------------------------------------
   2       double precision function sscale(r)
   3       double precision r,gamm
   4       include "COMMON.SPLITELE"
   5       if(r.lt.r_cut-rlamb) then
   6         sscale=1.0d0
   7       else if(r.le.r_cut.and.r.ge.r_cut-rlamb) then
   8         gamm=(r-(r_cut-rlamb))/rlamb
   9         sscale=1.0d0+gamm*gamm*(2*gamm-3.0d0)
  10       else
  11         sscale=0d0
  12       endif
  13       return
  14       end
  15 C-----------------------------------------------------------------------
  16       subroutine elj_long(evdw)
  17 C
  18 C This subroutine calculates the interaction energy of nonbonded side chains
  19 C assuming the LJ potential of interaction.
  20 C
  21       implicit real*8 (a-h,o-z)
  22       include 'DIMENSIONS'
  23       parameter (accur=1.0d-10)
  24       include 'COMMON.GEO'
  25       include 'COMMON.VAR'
  26       include 'COMMON.LOCAL'
  27       include 'COMMON.CHAIN'
  28       include 'COMMON.DERIV'
  29       include 'COMMON.INTERACT'
  30       include 'COMMON.TORSION'
  31       include 'COMMON.SBRIDGE'
  32       include 'COMMON.NAMES'
  33       include 'COMMON.IOUNITS'
  34       include 'COMMON.CONTACTS'
  35       dimension gg(3)
  36 c      write(iout,*)'Entering ELJ nnt=',nnt,' nct=',nct,' expon=',expon
  37       evdw=0.0D0
  38       do i=iatsc_s,iatsc_e
  39         itypi=itype(i)
  40         itypi1=itype(i+1)
  41         xi=c(1,nres+i)
  42         yi=c(2,nres+i)
  43         zi=c(3,nres+i)
  44 C
  45 C Calculate SC interaction energy.
  46 C
  47         do iint=1,nint_gr(i)
  48 cd        write (iout,*) 'i=',i,' iint=',iint,' istart=',istart(i,iint),
  49 cd   &                  'iend=',iend(i,iint)
  50           do j=istart(i,iint),iend(i,iint)
  51             itypj=itype(j)
  52             xj=c(1,nres+j)-xi
  53             yj=c(2,nres+j)-yi
  54             zj=c(3,nres+j)-zi
  55             rij=xj*xj+yj*yj+zj*zj
  56             sss=sscale(dsqrt(rij)/sigma(itypi,itypj))
  57             if (sss.lt.1.0d0) then
  58               rrij=1.0D0/rij
  59               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 c            write(iout,*) "long",i,itypi,j,itypj," rij",1.0d0/rij,
 703 c     &          " 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 c              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 c            write(iout,*) "short",i,itypi,j,itypj," rij",1.0d0/rij,
 856 c     &          " 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 c              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 cd      do i=1,nres
1333 cd        write (iout,'(i3,3f10.5,5x,3f10.5)')
1334 cd     &     i,(gel_loc(k,i),k=1,3),gel_loc_loc(i)
1335 cd      enddo
1336 c
1337 c 9/27/08 AL Split the interaction loop to ensure load balancing of turn terms
1338 C
1339 C Loop over i,i+2 and i,i+3 pairs of the peptide groups
1340 C
1341       do i=iturn3_start,iturn3_end
1342         dxi=dc(1,i)
1343         dyi=dc(2,i)
1344         dzi=dc(3,i)
1345         dx_normi=dc_norm(1,i)
1346         dy_normi=dc_norm(2,i)
1347         dz_normi=dc_norm(3,i)
1348         xmedi=c(1,i)+0.5d0*dxi
1349         ymedi=c(2,i)+0.5d0*dyi
1350         zmedi=c(3,i)+0.5d0*dzi
1351         num_conti=0
1352         call eelecij_scale(i,i+2,ees,evdw1,eel_loc)
1353         if (wturn3.gt.0.0d0) call eturn3(i,eello_turn3)
1354         num_cont_hb(i)=num_conti
1355       enddo
1356       do i=iturn4_start,iturn4_end
1357         dxi=dc(1,i)
1358         dyi=dc(2,i)
1359         dzi=dc(3,i)
1360         dx_normi=dc_norm(1,i)
1361         dy_normi=dc_norm(2,i)
1362         dz_normi=dc_norm(3,i)
1363         xmedi=c(1,i)+0.5d0*dxi
1364         ymedi=c(2,i)+0.5d0*dyi
1365         zmedi=c(3,i)+0.5d0*dzi
1366         num_conti=0
1367         call eelecij_scale(i,i+3,ees,evdw1,eel_loc)
1368         if (wturn4.gt.0.0d0) call eturn4(i,eello_turn4)
1369         num_cont_hb(i)=num_cont_hb(i)+num_conti
1370       enddo   ! i
1371 c
1372 c Loop over all pairs of interacting peptide groups except i,i+2 and i,i+3
1373 c
1374       do i=iatel_s,iatel_e
1375         dxi=dc(1,i)
1376         dyi=dc(2,i)
1377         dzi=dc(3,i)
1378         dx_normi=dc_norm(1,i)
1379         dy_normi=dc_norm(2,i)
1380         dz_normi=dc_norm(3,i)
1381         xmedi=c(1,i)+0.5d0*dxi
1382         ymedi=c(2,i)+0.5d0*dyi
1383         zmedi=c(3,i)+0.5d0*dzi
1384         num_conti=0
1385 c        write (iout,*) 'i',i,' ielstart',ielstart(i),' ielend',ielend(i)
1386         do j=ielstart(i),ielend(i)
1387           call eelecij_scale(i,j,ees,evdw1,eel_loc)
1388         enddo ! j
1389         num_cont_hb(i)=num_cont_hb(i)+num_conti
1390       enddo   ! i
1391       return
1392       end
1393 C-------------------------------------------------------------------------------
1394       subroutine eelecij_scale(i,j,ees,evdw1,eel_loc)
1395 C
1396 C This subroutine calculates the average interaction energy and its gradient
1397 C in the virtual-bond vectors between non-adjacent peptide groups, based on
1398 C the potential described in Liwo et al., Protein Sci., 1993, 2, 1715.
1399 C The potential depends both on the distance of peptide-group centers and on
1400 C the orientation of the CA-CA virtual bonds.
1401 C
1402       implicit real*8 (a-h,o-z)
1403       include 'DIMENSIONS'
1404       include 'COMMON.CONTROL'
1405       include 'COMMON.IOUNITS'
1406       include 'COMMON.GEO'
1407       include 'COMMON.VAR'
1408       include 'COMMON.LOCAL'
1409       include 'COMMON.CHAIN'
1410       include 'COMMON.DERIV'
1411       include 'COMMON.INTERACT'
1412       include 'COMMON.CONTACTS'
1413       include 'COMMON.TORSION'
1414       include 'COMMON.VECTORS'
1415       include 'COMMON.FFIELD'
1416       dimension ggg(3),gggp(3),gggm(3),erij(3),dcosb(3),dcosg(3),
1417      &          erder(3,3),uryg(3,3),urzg(3,3),vryg(3,3),vrzg(3,3)
1418       double precision acipa(2,2),agg(3,4),aggi(3,4),aggi1(3,4),
1419      &    aggj(3,4),aggj1(3,4),a_temp(2,2),muij(4)
1420       common /locel/ a_temp,agg,aggi,aggi1,aggj,aggj1,a22,a23,a32,a33,
1421      &    dxi,dyi,dzi,dx_normi,dy_normi,dz_normi,xmedi,ymedi,zmedi,
1422      &    num_conti,j1,j2
1423 c 4/26/02 - AL scaling factor for 1,4 repulsive VDW interactions
1424 #ifdef MOMENT
1425       double precision scal_el /1.0d0/
1426 #else
1427       double precision scal_el /0.5d0/
1428 #endif
1429 C 12/13/98
1430 C 13-go grudnia roku pamietnego...
1431       double precision unmat(3,3) /1.0d0,0.0d0,0.0d0,
1432      &                   0.0d0,1.0d0,0.0d0,
1433      &                   0.0d0,0.0d0,1.0d0/
1434           ind=ind+1
1435           iteli=itel(i)
1436           itelj=itel(j)
1437           if (j.eq.i+2 .and. itelj.eq.2) iteli=2
1438           aaa=app(iteli,itelj)
1439           bbb=bpp(iteli,itelj)
1440           ael6i=ael6(iteli,itelj)
1441           ael3i=ael3(iteli,itelj)
1442 C Diagnostics only!!!
1443 c         aaa=0.0D0
1444 c         bbb=0.0D0
1445 c         ael6i=0.0D0
1446 c         ael3i=0.0D0
1447 C End diagnostics
1448           dxj=dc(1,j)
1449           dyj=dc(2,j)
1450           dzj=dc(3,j)
1451           dx_normj=dc_norm(1,j)
1452           dy_normj=dc_norm(2,j)
1453           dz_normj=dc_norm(3,j)
1454           xj=c(1,j)+0.5D0*dxj-xmedi
1455           yj=c(2,j)+0.5D0*dyj-ymedi
1456           zj=c(3,j)+0.5D0*dzj-zmedi
1457           rij=xj*xj+yj*yj+zj*zj
1458           rrmij=1.0D0/rij
1459           rij=dsqrt(rij)
1460           rmij=1.0D0/rij
1461 c For extracting the short-range part of Evdwpp
1462           sss=sscale(rij/rpp(iteli,itelj))
1463 c
1464           r3ij=rrmij*rmij
1465           r6ij=r3ij*r3ij
1466           cosa=dx_normi*dx_normj+dy_normi*dy_normj+dz_normi*dz_normj
1467           cosb=(xj*dx_normi+yj*dy_normi+zj*dz_normi)*rmij
1468           cosg=(xj*dx_normj+yj*dy_normj+zj*dz_normj)*rmij
1469           fac=cosa-3.0D0*cosb*cosg
1470           ev1=aaa*r6ij*r6ij
1471 c 4/26/02 - AL scaling down 1,4 repulsive VDW interactions
1472           if (j.eq.i+2) ev1=scal_el*ev1
1473           ev2=bbb*r6ij
1474           fac3=ael6i*r6ij
1475           fac4=ael3i*r3ij
1476           evdwij=ev1+ev2
1477           el1=fac3*(4.0D0+fac*fac-3.0D0*(cosb*cosb+cosg*cosg))
1478           el2=fac4*fac
1479           eesij=el1+el2
1480 C 12/26/95 - for the evaluation of multi-body H-bonding interactions
1481           ees0ij=4.0D0+fac*fac-3.0D0*(cosb*cosb+cosg*cosg)
1482           ees=ees+eesij
1483           evdw1=evdw1+evdwij*(1.0d0-sss)
1484 cd          write(iout,'(2(2i3,2x),7(1pd12.4)/2(3(1pd12.4),5x)/)')
1485 cd     &      iteli,i,itelj,j,aaa,bbb,ael6i,ael3i,
1486 cd     &      1.0D0/dsqrt(rrmij),evdwij,eesij,
1487 cd     &      xmedi,ymedi,zmedi,xj,yj,zj
1488
1489           if (energy_dec) then
1490               write (iout,'(a6,2i5,0pf7.3)') 'evdw1',i,j,evdwij
1491               write (iout,'(a6,2i5,0pf7.3)') 'ees',i,j,eesij
1492           endif
1493
1494 C
1495 C Calculate contributions to the Cartesian gradient.
1496 C
1497 #ifdef SPLITELE
1498           facvdw=-6*rrmij*(ev1+evdwij)*(1.0d0-sss)
1499           facel=-3*rrmij*(el1+eesij)
1500           fac1=fac
1501           erij(1)=xj*rmij
1502           erij(2)=yj*rmij
1503           erij(3)=zj*rmij
1504 *
1505 * Radial derivatives. First process both termini of the fragment (i,j)
1506 *
1507           ggg(1)=facel*xj
1508           ggg(2)=facel*yj
1509           ggg(3)=facel*zj
1510           do k=1,3
1511             ghalf=0.5D0*ggg(k)
1512             gelc(k,i)=gelc(k,i)+ghalf
1513             gelc(k,j)=gelc(k,j)+ghalf
1514           enddo
1515 *
1516 * Loop over residues i+1 thru j-1.
1517 *
1518           do k=i+1,j-1
1519             do l=1,3
1520               gelc(l,k)=gelc(l,k)+ggg(l)
1521             enddo
1522           enddo
1523           ggg(1)=facvdw*xj
1524           ggg(2)=facvdw*yj
1525           ggg(3)=facvdw*zj
1526           do k=1,3
1527             ghalf=0.5D0*ggg(k)
1528             gvdwpp(k,i)=gvdwpp(k,i)+ghalf
1529             gvdwpp(k,j)=gvdwpp(k,j)+ghalf
1530           enddo
1531 *
1532 * Loop over residues i+1 thru j-1.
1533 *
1534           do k=i+1,j-1
1535             do l=1,3
1536               gvdwpp(l,k)=gvdwpp(l,k)+ggg(l)
1537             enddo
1538           enddo
1539 #else
1540           facvdw=(ev1+evdwij)*(1.0d0-sss)
1541           facel=el1+eesij
1542           fac1=fac
1543           fac=-3*rrmij*(facvdw+facvdw+facel)
1544           erij(1)=xj*rmij
1545           erij(2)=yj*rmij
1546           erij(3)=zj*rmij
1547 *
1548 * Radial derivatives. First process both termini of the fragment (i,j)
1549 *
1550           ggg(1)=fac*xj
1551           ggg(2)=fac*yj
1552           ggg(3)=fac*zj
1553           do k=1,3
1554             ghalf=0.5D0*ggg(k)
1555             gelc(k,i)=gelc(k,i)+ghalf
1556             gelc(k,j)=gelc(k,j)+ghalf
1557           enddo
1558 *
1559 * Loop over residues i+1 thru j-1.
1560 *
1561           do k=i+1,j-1
1562             do l=1,3
1563               gelc(l,k)=gelc(l,k)+ggg(l)
1564             enddo
1565           enddo
1566 #endif
1567 *
1568 * Angular part
1569 *
1570           ecosa=2.0D0*fac3*fac1+fac4
1571           fac4=-3.0D0*fac4
1572           fac3=-6.0D0*fac3
1573           ecosb=(fac3*(fac1*cosg+cosb)+cosg*fac4)
1574           ecosg=(fac3*(fac1*cosb+cosg)+cosb*fac4)
1575           do k=1,3
1576             dcosb(k)=rmij*(dc_norm(k,i)-erij(k)*cosb)
1577             dcosg(k)=rmij*(dc_norm(k,j)-erij(k)*cosg)
1578           enddo
1579 cd        print '(2i3,2(3(1pd14.5),3x))',i,j,(dcosb(k),k=1,3),
1580 cd   &          (dcosg(k),k=1,3)
1581           do k=1,3
1582             ggg(k)=ecosb*dcosb(k)+ecosg*dcosg(k)
1583           enddo
1584           do k=1,3
1585             ghalf=0.5D0*ggg(k)
1586             gelc(k,i)=gelc(k,i)+ghalf
1587      &               +(ecosa*(dc_norm(k,j)-cosa*dc_norm(k,i))
1588      &               + ecosb*(erij(k)-cosb*dc_norm(k,i)))*vbld_inv(i+1)
1589             gelc(k,j)=gelc(k,j)+ghalf
1590      &               +(ecosa*(dc_norm(k,i)-cosa*dc_norm(k,j))
1591      &               + ecosg*(erij(k)-cosg*dc_norm(k,j)))*vbld_inv(j+1)
1592           enddo
1593           do k=i+1,j-1
1594             do l=1,3
1595               gelc(l,k)=gelc(l,k)+ggg(l)
1596             enddo
1597           enddo
1598
1599           IF (wel_loc.gt.0.0d0 .or. wcorr4.gt.0.0d0 .or. wcorr5.gt.0.0d0
1600      &        .or. wcorr6.gt.0.0d0 .or. wturn3.gt.0.0d0
1601      &        .or. wturn4.gt.0.0d0 .or. wturn6.gt.0.0d0) THEN
1602 C
1603 C 9/25/99 Mixed third-order local-electrostatic terms. The local-interaction
1604 C   energy of a peptide unit is assumed in the form of a second-order
1605 C   Fourier series in the angles lambda1 and lambda2 (see Nishikawa et al.
1606 C   Macromolecules, 1974, 7, 797-806 for definition). This correlation terms
1607 C   are computed for EVERY pair of non-contiguous peptide groups.
1608 C
1609           if (j.lt.nres-1) then
1610             j1=j+1
1611             j2=j-1
1612           else
1613             j1=j-1
1614             j2=j-2
1615           endif
1616           kkk=0
1617           do k=1,2
1618             do l=1,2
1619               kkk=kkk+1
1620               muij(kkk)=mu(k,i)*mu(l,j)
1621             enddo
1622           enddo
1623 cd         write (iout,*) 'EELEC: i',i,' j',j
1624 cd          write (iout,*) 'j',j,' j1',j1,' j2',j2
1625 cd          write(iout,*) 'muij',muij
1626           ury=scalar(uy(1,i),erij)
1627           urz=scalar(uz(1,i),erij)
1628           vry=scalar(uy(1,j),erij)
1629           vrz=scalar(uz(1,j),erij)
1630           a22=scalar(uy(1,i),uy(1,j))-3*ury*vry
1631           a23=scalar(uy(1,i),uz(1,j))-3*ury*vrz
1632           a32=scalar(uz(1,i),uy(1,j))-3*urz*vry
1633           a33=scalar(uz(1,i),uz(1,j))-3*urz*vrz
1634 C For diagnostics only
1635 cd          a22=1.0d0
1636 cd          a23=1.0d0
1637 cd          a32=1.0d0
1638 cd          a33=1.0d0
1639           fac=dsqrt(-ael6i)*r3ij
1640 cd          write (2,*) 'fac=',fac
1641 C For diagnostics only
1642 cd          fac=1.0d0
1643           a22=a22*fac
1644           a23=a23*fac
1645           a32=a32*fac
1646           a33=a33*fac
1647 cd          write (iout,'(4i5,4f10.5)')
1648 cd     &     i,itortyp(itype(i)),j,itortyp(itype(j)),a22,a23,a32,a33
1649 cd          write (iout,'(6f10.5)') (muij(k),k=1,4),fac,eel_loc_ij
1650 cd          write (iout,'(2(3f10.5,5x)/2(3f10.5,5x))') uy(:,i),uz(:,i),
1651 cd     &      uy(:,j),uz(:,j)
1652 cd          write (iout,'(4f10.5)')
1653 cd     &      scalar(uy(1,i),uy(1,j)),scalar(uy(1,i),uz(1,j)),
1654 cd     &      scalar(uz(1,i),uy(1,j)),scalar(uz(1,i),uz(1,j))
1655 cd          write (iout,'(4f10.5)') ury,urz,vry,vrz
1656 cd           write (iout,'(9f10.5/)')
1657 cd     &      fac22,a22,fac23,a23,fac32,a32,fac33,a33,eel_loc_ij
1658 C Derivatives of the elements of A in virtual-bond vectors
1659           call unormderiv(erij(1),unmat(1,1),rmij,erder(1,1))
1660 cd          do k=1,3
1661 cd            do l=1,3
1662 cd              erder(k,l)=0.0d0
1663 cd            enddo
1664 cd          enddo
1665           do k=1,3
1666             uryg(k,1)=scalar(erder(1,k),uy(1,i))
1667             uryg(k,2)=scalar(uygrad(1,k,1,i),erij(1))
1668             uryg(k,3)=scalar(uygrad(1,k,2,i),erij(1))
1669             urzg(k,1)=scalar(erder(1,k),uz(1,i))
1670             urzg(k,2)=scalar(uzgrad(1,k,1,i),erij(1))
1671             urzg(k,3)=scalar(uzgrad(1,k,2,i),erij(1))
1672             vryg(k,1)=scalar(erder(1,k),uy(1,j))
1673             vryg(k,2)=scalar(uygrad(1,k,1,j),erij(1))
1674             vryg(k,3)=scalar(uygrad(1,k,2,j),erij(1))
1675             vrzg(k,1)=scalar(erder(1,k),uz(1,j))
1676             vrzg(k,2)=scalar(uzgrad(1,k,1,j),erij(1))
1677             vrzg(k,3)=scalar(uzgrad(1,k,2,j),erij(1))
1678           enddo
1679 cd          do k=1,3
1680 cd            do l=1,3
1681 cd              uryg(k,l)=0.0d0
1682 cd              urzg(k,l)=0.0d0
1683 cd              vryg(k,l)=0.0d0
1684 cd              vrzg(k,l)=0.0d0
1685 cd            enddo
1686 cd          enddo
1687 C Compute radial contributions to the gradient
1688           facr=-3.0d0*rrmij
1689           a22der=a22*facr
1690           a23der=a23*facr
1691           a32der=a32*facr
1692           a33der=a33*facr
1693 cd          a22der=0.0d0
1694 cd          a23der=0.0d0
1695 cd          a32der=0.0d0
1696 cd          a33der=0.0d0
1697           agg(1,1)=a22der*xj
1698           agg(2,1)=a22der*yj
1699           agg(3,1)=a22der*zj
1700           agg(1,2)=a23der*xj
1701           agg(2,2)=a23der*yj
1702           agg(3,2)=a23der*zj
1703           agg(1,3)=a32der*xj
1704           agg(2,3)=a32der*yj
1705           agg(3,3)=a32der*zj
1706           agg(1,4)=a33der*xj
1707           agg(2,4)=a33der*yj
1708           agg(3,4)=a33der*zj
1709 C Add the contributions coming from er
1710           fac3=-3.0d0*fac
1711           do k=1,3
1712             agg(k,1)=agg(k,1)+fac3*(uryg(k,1)*vry+vryg(k,1)*ury)
1713             agg(k,2)=agg(k,2)+fac3*(uryg(k,1)*vrz+vrzg(k,1)*ury)
1714             agg(k,3)=agg(k,3)+fac3*(urzg(k,1)*vry+vryg(k,1)*urz)
1715             agg(k,4)=agg(k,4)+fac3*(urzg(k,1)*vrz+vrzg(k,1)*urz)
1716           enddo
1717           do k=1,3
1718 C Derivatives in DC(i)
1719             ghalf1=0.5d0*agg(k,1)
1720             ghalf2=0.5d0*agg(k,2)
1721             ghalf3=0.5d0*agg(k,3)
1722             ghalf4=0.5d0*agg(k,4)
1723             aggi(k,1)=fac*(scalar(uygrad(1,k,1,i),uy(1,j))
1724      &      -3.0d0*uryg(k,2)*vry)+ghalf1
1725             aggi(k,2)=fac*(scalar(uygrad(1,k,1,i),uz(1,j))
1726      &      -3.0d0*uryg(k,2)*vrz)+ghalf2
1727             aggi(k,3)=fac*(scalar(uzgrad(1,k,1,i),uy(1,j))
1728      &      -3.0d0*urzg(k,2)*vry)+ghalf3
1729             aggi(k,4)=fac*(scalar(uzgrad(1,k,1,i),uz(1,j))
1730      &      -3.0d0*urzg(k,2)*vrz)+ghalf4
1731 C Derivatives in DC(i+1)
1732             aggi1(k,1)=fac*(scalar(uygrad(1,k,2,i),uy(1,j))
1733      &      -3.0d0*uryg(k,3)*vry)+agg(k,1)
1734             aggi1(k,2)=fac*(scalar(uygrad(1,k,2,i),uz(1,j))
1735      &      -3.0d0*uryg(k,3)*vrz)+agg(k,2)
1736             aggi1(k,3)=fac*(scalar(uzgrad(1,k,2,i),uy(1,j))
1737      &      -3.0d0*urzg(k,3)*vry)+agg(k,3)
1738             aggi1(k,4)=fac*(scalar(uzgrad(1,k,2,i),uz(1,j))
1739      &      -3.0d0*urzg(k,3)*vrz)+agg(k,4)
1740 C Derivatives in DC(j)
1741             aggj(k,1)=fac*(scalar(uygrad(1,k,1,j),uy(1,i))
1742      &      -3.0d0*vryg(k,2)*ury)+ghalf1
1743             aggj(k,2)=fac*(scalar(uzgrad(1,k,1,j),uy(1,i))
1744      &      -3.0d0*vrzg(k,2)*ury)+ghalf2
1745             aggj(k,3)=fac*(scalar(uygrad(1,k,1,j),uz(1,i))
1746      &      -3.0d0*vryg(k,2)*urz)+ghalf3
1747             aggj(k,4)=fac*(scalar(uzgrad(1,k,1,j),uz(1,i))
1748      &      -3.0d0*vrzg(k,2)*urz)+ghalf4
1749 C Derivatives in DC(j+1) or DC(nres-1)
1750             aggj1(k,1)=fac*(scalar(uygrad(1,k,2,j),uy(1,i))
1751      &      -3.0d0*vryg(k,3)*ury)
1752             aggj1(k,2)=fac*(scalar(uzgrad(1,k,2,j),uy(1,i))
1753      &      -3.0d0*vrzg(k,3)*ury)
1754             aggj1(k,3)=fac*(scalar(uygrad(1,k,2,j),uz(1,i))
1755      &      -3.0d0*vryg(k,3)*urz)
1756             aggj1(k,4)=fac*(scalar(uzgrad(1,k,2,j),uz(1,i))
1757      &      -3.0d0*vrzg(k,3)*urz)
1758 cd            aggi(k,1)=ghalf1
1759 cd            aggi(k,2)=ghalf2
1760 cd            aggi(k,3)=ghalf3
1761 cd            aggi(k,4)=ghalf4
1762 C Derivatives in DC(i+1)
1763 cd            aggi1(k,1)=agg(k,1)
1764 cd            aggi1(k,2)=agg(k,2)
1765 cd            aggi1(k,3)=agg(k,3)
1766 cd            aggi1(k,4)=agg(k,4)
1767 C Derivatives in DC(j)
1768 cd            aggj(k,1)=ghalf1
1769 cd            aggj(k,2)=ghalf2
1770 cd            aggj(k,3)=ghalf3
1771 cd            aggj(k,4)=ghalf4
1772 C Derivatives in DC(j+1)
1773 cd            aggj1(k,1)=0.0d0
1774 cd            aggj1(k,2)=0.0d0
1775 cd            aggj1(k,3)=0.0d0
1776 cd            aggj1(k,4)=0.0d0
1777             if (j.eq.nres-1 .and. i.lt.j-2) then
1778               do l=1,4
1779                 aggj1(k,l)=aggj1(k,l)+agg(k,l)
1780 cd                aggj1(k,l)=agg(k,l)
1781               enddo
1782             endif
1783           enddo
1784 c          goto 11111
1785 C Check the loc-el terms by numerical integration
1786           acipa(1,1)=a22
1787           acipa(1,2)=a23
1788           acipa(2,1)=a32
1789           acipa(2,2)=a33
1790           a22=-a22
1791           a23=-a23
1792           do l=1,2
1793             do k=1,3
1794               agg(k,l)=-agg(k,l)
1795               aggi(k,l)=-aggi(k,l)
1796               aggi1(k,l)=-aggi1(k,l)
1797               aggj(k,l)=-aggj(k,l)
1798               aggj1(k,l)=-aggj1(k,l)
1799             enddo
1800           enddo
1801           if (j.lt.nres-1) then
1802             a22=-a22
1803             a32=-a32
1804             do l=1,3,2
1805               do k=1,3
1806                 agg(k,l)=-agg(k,l)
1807                 aggi(k,l)=-aggi(k,l)
1808                 aggi1(k,l)=-aggi1(k,l)
1809                 aggj(k,l)=-aggj(k,l)
1810                 aggj1(k,l)=-aggj1(k,l)
1811               enddo
1812             enddo
1813           else
1814             a22=-a22
1815             a23=-a23
1816             a32=-a32
1817             a33=-a33
1818             do l=1,4
1819               do k=1,3
1820                 agg(k,l)=-agg(k,l)
1821                 aggi(k,l)=-aggi(k,l)
1822                 aggi1(k,l)=-aggi1(k,l)
1823                 aggj(k,l)=-aggj(k,l)
1824                 aggj1(k,l)=-aggj1(k,l)
1825               enddo
1826             enddo
1827           endif
1828           ENDIF ! WCORR
1829 11111     continue
1830           IF (wel_loc.gt.0.0d0) THEN
1831 C Contribution to the local-electrostatic energy coming from the i-j pair
1832           eel_loc_ij=a22*muij(1)+a23*muij(2)+a32*muij(3)
1833      &     +a33*muij(4)
1834 cd          write (iout,*) 'i',i,' j',j,' eel_loc_ij',eel_loc_ij
1835
1836           if (energy_dec) write (iout,'(a6,2i5,0pf7.3)')
1837      &            'eelloc',i,j,eel_loc_ij
1838
1839           eel_loc=eel_loc+eel_loc_ij
1840 C Partial derivatives in virtual-bond dihedral angles gamma
1841           if (i.gt.1)
1842      &    gel_loc_loc(i-1)=gel_loc_loc(i-1)+
1843      &            a22*muder(1,i)*mu(1,j)+a23*muder(1,i)*mu(2,j)
1844      &           +a32*muder(2,i)*mu(1,j)+a33*muder(2,i)*mu(2,j)
1845           gel_loc_loc(j-1)=gel_loc_loc(j-1)+
1846      &            a22*mu(1,i)*muder(1,j)+a23*mu(1,i)*muder(2,j)
1847      &           +a32*mu(2,i)*muder(1,j)+a33*mu(2,i)*muder(2,j)
1848 cd          call checkint3(i,j,mu1,mu2,a22,a23,a32,a33,acipa,eel_loc_ij)
1849 cd          write(iout,*) 'agg  ',agg
1850 cd          write(iout,*) 'aggi ',aggi
1851 cd          write(iout,*) 'aggi1',aggi1
1852 cd          write(iout,*) 'aggj ',aggj
1853 cd          write(iout,*) 'aggj1',aggj1
1854
1855 C Derivatives of eello in DC(i+1) thru DC(j-1) or DC(nres-2)
1856           do l=1,3
1857             ggg(l)=agg(l,1)*muij(1)+
1858      &          agg(l,2)*muij(2)+agg(l,3)*muij(3)+agg(l,4)*muij(4)
1859           enddo
1860           do k=i+2,j2
1861             do l=1,3
1862               gel_loc(l,k)=gel_loc(l,k)+ggg(l)
1863             enddo
1864           enddo
1865 C Remaining derivatives of eello
1866           do l=1,3
1867             gel_loc(l,i)=gel_loc(l,i)+aggi(l,1)*muij(1)+
1868      &          aggi(l,2)*muij(2)+aggi(l,3)*muij(3)+aggi(l,4)*muij(4)
1869             gel_loc(l,i+1)=gel_loc(l,i+1)+aggi1(l,1)*muij(1)+
1870      &          aggi1(l,2)*muij(2)+aggi1(l,3)*muij(3)+aggi1(l,4)*muij(4)
1871             gel_loc(l,j)=gel_loc(l,j)+aggj(l,1)*muij(1)+
1872      &          aggj(l,2)*muij(2)+aggj(l,3)*muij(3)+aggj(l,4)*muij(4)
1873             gel_loc(l,j1)=gel_loc(l,j1)+aggj1(l,1)*muij(1)+
1874      &          aggj1(l,2)*muij(2)+aggj1(l,3)*muij(3)+aggj1(l,4)*muij(4)
1875           enddo
1876           ENDIF
1877 C Change 12/26/95 to calculate four-body contributions to H-bonding energy
1878           if (j.gt.i+1 .and. num_conti.le.maxconts) then
1879 C
1880 C Calculate the contact function. The ith column of the array JCONT will
1881 C contain the numbers of atoms that make contacts with the atom I (of numbers
1882 C greater than I). The arrays FACONT and GACONT will contain the values of
1883 C the contact function and its derivative.
1884 c           r0ij=1.02D0*rpp(iteli,itelj)
1885 c           r0ij=1.11D0*rpp(iteli,itelj)
1886             r0ij=2.20D0*rpp(iteli,itelj)
1887 c           r0ij=1.55D0*rpp(iteli,itelj)
1888             call gcont(rij,r0ij,1.0D0,0.2d0*r0ij,fcont,fprimcont)
1889             if (fcont.gt.0.0D0) then
1890               num_conti=num_conti+1
1891               if (num_conti.gt.maxconts) then
1892                 write (iout,*) 'WARNING - max. # of contacts exceeded;',
1893      &                         ' will skip next contacts for this conf.'
1894               else
1895                 jcont_hb(num_conti,i)=j
1896                 IF (wcorr4.gt.0.0d0 .or. wcorr5.gt.0.0d0 .or.
1897      &          wcorr6.gt.0.0d0 .or. wturn6.gt.0.0d0) THEN
1898 C 9/30/99 (AL) - store components necessary to evaluate higher-order loc-el
1899 C  terms.
1900                 d_cont(num_conti,i)=rij
1901 cd                write (2,'(3e15.5)') rij,r0ij+0.2d0*r0ij,rij
1902 C     --- Electrostatic-interaction matrix ---
1903                 a_chuj(1,1,num_conti,i)=a22
1904                 a_chuj(1,2,num_conti,i)=a23
1905                 a_chuj(2,1,num_conti,i)=a32
1906                 a_chuj(2,2,num_conti,i)=a33
1907 C     --- Gradient of rij
1908                 do kkk=1,3
1909                   grij_hb_cont(kkk,num_conti,i)=erij(kkk)
1910                 enddo
1911 c             if (i.eq.1) then
1912 c                a_chuj(1,1,num_conti,i)=-0.61d0
1913 c                a_chuj(1,2,num_conti,i)= 0.4d0
1914 c                a_chuj(2,1,num_conti,i)= 0.65d0
1915 c                a_chuj(2,2,num_conti,i)= 0.50d0
1916 c             else if (i.eq.2) then
1917 c                a_chuj(1,1,num_conti,i)= 0.0d0
1918 c                a_chuj(1,2,num_conti,i)= 0.0d0
1919 c                a_chuj(2,1,num_conti,i)= 0.0d0
1920 c                a_chuj(2,2,num_conti,i)= 0.0d0
1921 c             endif
1922 C     --- and its gradients
1923 cd                write (iout,*) 'i',i,' j',j
1924 cd                do kkk=1,3
1925 cd                write (iout,*) 'iii 1 kkk',kkk
1926 cd                write (iout,*) agg(kkk,:)
1927 cd                enddo
1928 cd                do kkk=1,3
1929 cd                write (iout,*) 'iii 2 kkk',kkk
1930 cd                write (iout,*) aggi(kkk,:)
1931 cd                enddo
1932 cd                do kkk=1,3
1933 cd                write (iout,*) 'iii 3 kkk',kkk
1934 cd                write (iout,*) aggi1(kkk,:)
1935 cd                enddo
1936 cd                do kkk=1,3
1937 cd                write (iout,*) 'iii 4 kkk',kkk
1938 cd                write (iout,*) aggj(kkk,:)
1939 cd                enddo
1940 cd                do kkk=1,3
1941 cd                write (iout,*) 'iii 5 kkk',kkk
1942 cd                write (iout,*) aggj1(kkk,:)
1943 cd                enddo
1944                 kkll=0
1945                 do k=1,2
1946                   do l=1,2
1947                     kkll=kkll+1
1948                     do m=1,3
1949                       a_chuj_der(k,l,m,1,num_conti,i)=agg(m,kkll)
1950                       a_chuj_der(k,l,m,2,num_conti,i)=aggi(m,kkll)
1951                       a_chuj_der(k,l,m,3,num_conti,i)=aggi1(m,kkll)
1952                       a_chuj_der(k,l,m,4,num_conti,i)=aggj(m,kkll)
1953                       a_chuj_der(k,l,m,5,num_conti,i)=aggj1(m,kkll)
1954 c                      do mm=1,5
1955 c                      a_chuj_der(k,l,m,mm,num_conti,i)=0.0d0
1956 c                      enddo
1957                     enddo
1958                   enddo
1959                 enddo
1960                 ENDIF
1961                 IF (wcorr4.eq.0.0d0 .and. wcorr.gt.0.0d0) THEN
1962 C Calculate contact energies
1963                 cosa4=4.0D0*cosa
1964                 wij=cosa-3.0D0*cosb*cosg
1965                 cosbg1=cosb+cosg
1966                 cosbg2=cosb-cosg
1967 c               fac3=dsqrt(-ael6i)/r0ij**3
1968                 fac3=dsqrt(-ael6i)*r3ij
1969 c                 ees0pij=dsqrt(4.0D0+cosa4+wij*wij-3.0D0*cosbg1*cosbg1)
1970                 ees0tmp=4.0D0+cosa4+wij*wij-3.0D0*cosbg1*cosbg1
1971                 if (ees0tmp.gt.0) then
1972                   ees0pij=dsqrt(ees0tmp)
1973                 else
1974                   ees0pij=0
1975                 endif
1976 c                ees0mij=dsqrt(4.0D0-cosa4+wij*wij-3.0D0*cosbg2*cosbg2)
1977                 ees0tmp=4.0D0-cosa4+wij*wij-3.0D0*cosbg2*cosbg2
1978                 if (ees0tmp.gt.0) then
1979                   ees0mij=dsqrt(ees0tmp)
1980                 else
1981                   ees0mij=0
1982                 endif
1983 c               ees0mij=0.0D0
1984                 ees0p(num_conti,i)=0.5D0*fac3*(ees0pij+ees0mij)
1985                 ees0m(num_conti,i)=0.5D0*fac3*(ees0pij-ees0mij)
1986 C Diagnostics. Comment out or remove after debugging!
1987 c               ees0p(num_conti,i)=0.5D0*fac3*ees0pij
1988 c               ees0m(num_conti,i)=0.5D0*fac3*ees0mij
1989 c               ees0m(num_conti,i)=0.0D0
1990 C End diagnostics.
1991 c               write (iout,*) 'i=',i,' j=',j,' rij=',rij,' r0ij=',r0ij,
1992 c    & ' ees0ij=',ees0p(num_conti,i),ees0m(num_conti,i),' fcont=',fcont
1993 C Angular derivatives of the contact function
1994                 ees0pij1=fac3/ees0pij
1995                 ees0mij1=fac3/ees0mij
1996                 fac3p=-3.0D0*fac3*rrmij
1997                 ees0pijp=0.5D0*fac3p*(ees0pij+ees0mij)
1998                 ees0mijp=0.5D0*fac3p*(ees0pij-ees0mij)
1999 c               ees0mij1=0.0D0
2000                 ecosa1=       ees0pij1*( 1.0D0+0.5D0*wij)
2001                 ecosb1=-1.5D0*ees0pij1*(wij*cosg+cosbg1)
2002                 ecosg1=-1.5D0*ees0pij1*(wij*cosb+cosbg1)
2003                 ecosa2=       ees0mij1*(-1.0D0+0.5D0*wij)
2004                 ecosb2=-1.5D0*ees0mij1*(wij*cosg+cosbg2)
2005                 ecosg2=-1.5D0*ees0mij1*(wij*cosb-cosbg2)
2006                 ecosap=ecosa1+ecosa2
2007                 ecosbp=ecosb1+ecosb2
2008                 ecosgp=ecosg1+ecosg2
2009                 ecosam=ecosa1-ecosa2
2010                 ecosbm=ecosb1-ecosb2
2011                 ecosgm=ecosg1-ecosg2
2012 C Diagnostics
2013 c               ecosap=ecosa1
2014 c               ecosbp=ecosb1
2015 c               ecosgp=ecosg1
2016 c               ecosam=0.0D0
2017 c               ecosbm=0.0D0
2018 c               ecosgm=0.0D0
2019 C End diagnostics
2020                 facont_hb(num_conti,i)=fcont
2021                 fprimcont=fprimcont/rij
2022 cd              facont_hb(num_conti,i)=1.0D0
2023 C Following line is for diagnostics.
2024 cd              fprimcont=0.0D0
2025                 do k=1,3
2026                   dcosb(k)=rmij*(dc_norm(k,i)-erij(k)*cosb)
2027                   dcosg(k)=rmij*(dc_norm(k,j)-erij(k)*cosg)
2028                 enddo
2029                 do k=1,3
2030                   gggp(k)=ecosbp*dcosb(k)+ecosgp*dcosg(k)
2031                   gggm(k)=ecosbm*dcosb(k)+ecosgm*dcosg(k)
2032                 enddo
2033                 gggp(1)=gggp(1)+ees0pijp*xj
2034                 gggp(2)=gggp(2)+ees0pijp*yj
2035                 gggp(3)=gggp(3)+ees0pijp*zj
2036                 gggm(1)=gggm(1)+ees0mijp*xj
2037                 gggm(2)=gggm(2)+ees0mijp*yj
2038                 gggm(3)=gggm(3)+ees0mijp*zj
2039 C Derivatives due to the contact function
2040                 gacont_hbr(1,num_conti,i)=fprimcont*xj
2041                 gacont_hbr(2,num_conti,i)=fprimcont*yj
2042                 gacont_hbr(3,num_conti,i)=fprimcont*zj
2043                 do k=1,3
2044                   ghalfp=0.5D0*gggp(k)
2045                   ghalfm=0.5D0*gggm(k)
2046                   gacontp_hb1(k,num_conti,i)=ghalfp
2047      &              +(ecosap*(dc_norm(k,j)-cosa*dc_norm(k,i))
2048      &              + ecosbp*(erij(k)-cosb*dc_norm(k,i)))*vbld_inv(i+1)
2049                   gacontp_hb2(k,num_conti,i)=ghalfp
2050      &              +(ecosap*(dc_norm(k,i)-cosa*dc_norm(k,j))
2051      &              + ecosgp*(erij(k)-cosg*dc_norm(k,j)))*vbld_inv(j+1)
2052                   gacontp_hb3(k,num_conti,i)=gggp(k)
2053                   gacontm_hb1(k,num_conti,i)=ghalfm
2054      &              +(ecosam*(dc_norm(k,j)-cosa*dc_norm(k,i))
2055      &              + ecosbm*(erij(k)-cosb*dc_norm(k,i)))*vbld_inv(i+1)
2056                   gacontm_hb2(k,num_conti,i)=ghalfm
2057      &              +(ecosam*(dc_norm(k,i)-cosa*dc_norm(k,j))
2058      &              + ecosgm*(erij(k)-cosg*dc_norm(k,j)))*vbld_inv(j+1)
2059                   gacontm_hb3(k,num_conti,i)=gggm(k)
2060                 enddo
2061 C Diagnostics. Comment out or remove after debugging!
2062 cdiag           do k=1,3
2063 cdiag             gacontp_hb1(k,num_conti,i)=0.0D0
2064 cdiag             gacontp_hb2(k,num_conti,i)=0.0D0
2065 cdiag             gacontp_hb3(k,num_conti,i)=0.0D0
2066 cdiag             gacontm_hb1(k,num_conti,i)=0.0D0
2067 cdiag             gacontm_hb2(k,num_conti,i)=0.0D0
2068 cdiag             gacontm_hb3(k,num_conti,i)=0.0D0
2069 cdiag           enddo
2070               ENDIF ! wcorr
2071               endif  ! num_conti.le.maxconts
2072             endif  ! fcont.gt.0
2073           endif    ! j.gt.i+1
2074       return
2075       end
2076 C-----------------------------------------------------------------------
2077       subroutine evdwpp_long(evdw1)
2078 C
2079 C Compute Evdwpp
2080 C
2081       implicit real*8 (a-h,o-z)
2082       include 'DIMENSIONS'
2083       include 'COMMON.CONTROL'
2084       include 'COMMON.IOUNITS'
2085       include 'COMMON.GEO'
2086       include 'COMMON.VAR'
2087       include 'COMMON.LOCAL'
2088       include 'COMMON.CHAIN'
2089       include 'COMMON.DERIV'
2090       include 'COMMON.INTERACT'
2091       include 'COMMON.CONTACTS'
2092       include 'COMMON.TORSION'
2093       include 'COMMON.VECTORS'
2094       include 'COMMON.FFIELD'
2095       dimension ggg(3)
2096 c 4/26/02 - AL scaling factor for 1,4 repulsive VDW interactions
2097 #ifdef MOMENT
2098       double precision scal_el /1.0d0/
2099 #else
2100       double precision scal_el /0.5d0/
2101 #endif
2102       evdw1=0.0D0
2103       do i=iatel_s,iatel_e
2104         dxi=dc(1,i)
2105         dyi=dc(2,i)
2106         dzi=dc(3,i)
2107         dx_normi=dc_norm(1,i)
2108         dy_normi=dc_norm(2,i)
2109         dz_normi=dc_norm(3,i)
2110         xmedi=c(1,i)+0.5d0*dxi
2111         ymedi=c(2,i)+0.5d0*dyi
2112         zmedi=c(3,i)+0.5d0*dzi
2113         num_conti=0
2114 c        write (iout,*) 'i',i,' ielstart',ielstart(i),' ielend',ielend(i)
2115         do j=ielstart(i),ielend(i)
2116           ind=ind+1
2117           iteli=itel(i)
2118           itelj=itel(j)
2119           if (j.eq.i+2 .and. itelj.eq.2) iteli=2
2120           aaa=app(iteli,itelj)
2121           bbb=bpp(iteli,itelj)
2122           dxj=dc(1,j)
2123           dyj=dc(2,j)
2124           dzj=dc(3,j)
2125           dx_normj=dc_norm(1,j)
2126           dy_normj=dc_norm(2,j)
2127           dz_normj=dc_norm(3,j)
2128           xj=c(1,j)+0.5D0*dxj-xmedi
2129           yj=c(2,j)+0.5D0*dyj-ymedi
2130           zj=c(3,j)+0.5D0*dzj-zmedi
2131           rij=xj*xj+yj*yj+zj*zj
2132           rrmij=1.0D0/rij
2133           rij=dsqrt(rij)
2134           sss=sscale(rij/rpp(iteli,itelj))
2135           if (sss.lt.1.0d0) then
2136             rmij=1.0D0/rij
2137             r3ij=rrmij*rmij
2138             r6ij=r3ij*r3ij
2139             ev1=aaa*r6ij*r6ij
2140 c 4/26/02 - AL scaling down 1,4 repulsive VDW interactions
2141             if (j.eq.i+2) ev1=scal_el*ev1
2142             ev2=bbb*r6ij
2143             evdwij=ev1+ev2
2144             if (energy_dec) then
2145               write (iout,'(a6,2i5,0pf7.3)') 'evdw1',i,j,evdwij
2146               write (iout,'(a6,2i5,0pf7.3)') 'ees',i,j,eesij
2147             endif
2148             evdw1=evdw1+evdwij*(1.0d0-sss)
2149 C
2150 C Calculate contributions to the Cartesian gradient.
2151 C
2152             facvdw=-6*rrmij*(ev1+evdwij)*(1.0d0-sss)
2153             ggg(1)=facvdw*xj
2154             ggg(2)=facvdw*yj
2155             ggg(3)=facvdw*zj
2156
2157             do k=1,3
2158               ghalf=0.5D0*ggg(k)
2159               gvdwpp(k,i)=gvdwpp(k,i)+ghalf
2160               gvdwpp(k,j)=gvdwpp(k,j)+ghalf
2161             enddo
2162 *
2163 * Loop over residues i+1 thru j-1.
2164 *
2165             do k=i+1,j-1
2166               do l=1,3
2167                 gvdwpp(l,k)=gvdwpp(l,k)+ggg(l)
2168               enddo
2169             enddo
2170           endif
2171         enddo ! j
2172       enddo   ! i
2173       return
2174       end
2175 C-----------------------------------------------------------------------
2176       subroutine evdwpp_short(evdw1)
2177 C
2178 C Compute Evdwpp
2179 C
2180       implicit real*8 (a-h,o-z)
2181       include 'DIMENSIONS'
2182       include 'COMMON.CONTROL'
2183       include 'COMMON.IOUNITS'
2184       include 'COMMON.GEO'
2185       include 'COMMON.VAR'
2186       include 'COMMON.LOCAL'
2187       include 'COMMON.CHAIN'
2188       include 'COMMON.DERIV'
2189       include 'COMMON.INTERACT'
2190       include 'COMMON.CONTACTS'
2191       include 'COMMON.TORSION'
2192       include 'COMMON.VECTORS'
2193       include 'COMMON.FFIELD'
2194       dimension ggg(3)
2195 c 4/26/02 - AL scaling factor for 1,4 repulsive VDW interactions
2196 #ifdef MOMENT
2197       double precision scal_el /1.0d0/
2198 #else
2199       double precision scal_el /0.5d0/
2200 #endif
2201       evdw1=0.0D0
2202       do i=iatel_s,iatel_e
2203         dxi=dc(1,i)
2204         dyi=dc(2,i)
2205         dzi=dc(3,i)
2206         dx_normi=dc_norm(1,i)
2207         dy_normi=dc_norm(2,i)
2208         dz_normi=dc_norm(3,i)
2209         xmedi=c(1,i)+0.5d0*dxi
2210         ymedi=c(2,i)+0.5d0*dyi
2211         zmedi=c(3,i)+0.5d0*dzi
2212         num_conti=0
2213 c        write (iout,*) 'i',i,' ielstart',ielstart(i),' ielend',ielend(i)
2214         do j=ielstart(i),ielend(i)
2215           ind=ind+1
2216           iteli=itel(i)
2217           itelj=itel(j)
2218           if (j.eq.i+2 .and. itelj.eq.2) iteli=2
2219           aaa=app(iteli,itelj)
2220           bbb=bpp(iteli,itelj)
2221           dxj=dc(1,j)
2222           dyj=dc(2,j)
2223           dzj=dc(3,j)
2224           dx_normj=dc_norm(1,j)
2225           dy_normj=dc_norm(2,j)
2226           dz_normj=dc_norm(3,j)
2227           xj=c(1,j)+0.5D0*dxj-xmedi
2228           yj=c(2,j)+0.5D0*dyj-ymedi
2229           zj=c(3,j)+0.5D0*dzj-zmedi
2230           rij=xj*xj+yj*yj+zj*zj
2231           rrmij=1.0D0/rij
2232           rij=dsqrt(rij)
2233           sss=sscale(rij/rpp(iteli,itelj))
2234           if (sss.gt.0.0d0) then
2235             rmij=1.0D0/rij
2236             r3ij=rrmij*rmij
2237             r6ij=r3ij*r3ij
2238             ev1=aaa*r6ij*r6ij
2239 c 4/26/02 - AL scaling down 1,4 repulsive VDW interactions
2240             if (j.eq.i+2) ev1=scal_el*ev1
2241             ev2=bbb*r6ij
2242             evdwij=ev1+ev2
2243             if (energy_dec) then
2244               write (iout,'(a6,2i5,0pf7.3)') 'evdw1',i,j,evdwij
2245               write (iout,'(a6,2i5,0pf7.3)') 'ees',i,j,eesij
2246             endif
2247             evdw1=evdw1+evdwij*sss
2248 C
2249 C Calculate contributions to the Cartesian gradient.
2250 C
2251             facvdw=-6*rrmij*(ev1+evdwij)*sss
2252             ggg(1)=facvdw*xj
2253             ggg(2)=facvdw*yj
2254             ggg(3)=facvdw*zj
2255
2256             do k=1,3
2257               ghalf=0.5D0*ggg(k)
2258               gvdwpp(k,i)=gvdwpp(k,i)+ghalf
2259               gvdwpp(k,j)=gvdwpp(k,j)+ghalf
2260             enddo
2261 *
2262 * Loop over residues i+1 thru j-1.
2263 *
2264             do k=i+1,j-1
2265               do l=1,3
2266                 gvdwpp(l,k)=gvdwpp(l,k)+ggg(l)
2267               enddo
2268             enddo
2269           endif
2270         enddo ! j
2271       enddo   ! i
2272       return
2273       end
2274 C-----------------------------------------------------------------------------
2275       subroutine escp_long(evdw2,evdw2_14)
2276 C
2277 C This subroutine calculates the excluded-volume interaction energy between
2278 C peptide-group centers and side chains and its gradient in virtual-bond and
2279 C side-chain vectors.
2280 C
2281       implicit real*8 (a-h,o-z)
2282       include 'DIMENSIONS'
2283       include 'COMMON.GEO'
2284       include 'COMMON.VAR'
2285       include 'COMMON.LOCAL'
2286       include 'COMMON.CHAIN'
2287       include 'COMMON.DERIV'
2288       include 'COMMON.INTERACT'
2289       include 'COMMON.FFIELD'
2290       include 'COMMON.IOUNITS'
2291       include 'COMMON.CONTROL'
2292       dimension ggg(3)
2293       evdw2=0.0D0
2294       evdw2_14=0.0d0
2295 cd    print '(a)','Enter ESCP'
2296 cd    write (iout,*) 'iatscp_s=',iatscp_s,' iatscp_e=',iatscp_e
2297       do i=iatscp_s,iatscp_e
2298         iteli=itel(i)
2299         xi=0.5D0*(c(1,i)+c(1,i+1))
2300         yi=0.5D0*(c(2,i)+c(2,i+1))
2301         zi=0.5D0*(c(3,i)+c(3,i+1))
2302
2303         do iint=1,nscp_gr(i)
2304
2305         do j=iscpstart(i,iint),iscpend(i,iint)
2306           itypj=itype(j)
2307 C Uncomment following three lines for SC-p interactions
2308 c         xj=c(1,nres+j)-xi
2309 c         yj=c(2,nres+j)-yi
2310 c         zj=c(3,nres+j)-zi
2311 C Uncomment following three lines for Ca-p interactions
2312           xj=c(1,j)-xi
2313           yj=c(2,j)-yi
2314           zj=c(3,j)-zi
2315           rrij=1.0D0/(xj*xj+yj*yj+zj*zj)
2316
2317           sss=sscale(1.0d0/(dsqrt(rrij)*rscp(itypj,iteli)))
2318
2319           if (sss.lt.1.0d0) then
2320
2321             fac=rrij**expon2
2322             e1=fac*fac*aad(itypj,iteli)
2323             e2=fac*bad(itypj,iteli)
2324             if (iabs(j-i) .le. 2) then
2325               e1=scal14*e1
2326               e2=scal14*e2
2327               evdw2_14=evdw2_14+(e1+e2)*(1.0d0-sss)
2328             endif
2329             evdwij=e1+e2
2330             evdw2=evdw2+evdwij*(1.0d0-sss)
2331             if (energy_dec) write (iout,'(a6,2i5,0pf7.3)')
2332      &          'evdw2',i,j,evdwij
2333 C
2334 C Calculate contributions to the gradient in the virtual-bond and SC vectors.
2335 C
2336             fac=-(evdwij+e1)*rrij*(1.0d0-sss)
2337             ggg(1)=xj*fac
2338             ggg(2)=yj*fac
2339             ggg(3)=zj*fac
2340             if (j.lt.i) then
2341 cd            write (iout,*) 'j<i'
2342 C Uncomment following three lines for SC-p interactions
2343 c             do k=1,3
2344 c             gradx_scp(k,j)=gradx_scp(k,j)+ggg(k)
2345 c             enddo
2346             else
2347 cd            write (iout,*) 'j>i'
2348               do k=1,3
2349                 ggg(k)=-ggg(k)
2350 C Uncomment following line for SC-p interactions
2351 c               gradx_scp(k,j)=gradx_scp(k,j)-ggg(k)
2352               enddo
2353             endif
2354             do k=1,3
2355               gvdwc_scp(k,i)=gvdwc_scp(k,i)-0.5D0*ggg(k)
2356             enddo
2357             kstart=min0(i+1,j)
2358             kend=max0(i-1,j-1)
2359 cd        write (iout,*) 'i=',i,' j=',j,' kstart=',kstart,' kend=',kend
2360 cd        write (iout,*) ggg(1),ggg(2),ggg(3)
2361             do k=kstart,kend
2362               do l=1,3
2363                 gvdwc_scp(l,k)=gvdwc_scp(l,k)-ggg(l)
2364               enddo
2365             enddo
2366
2367           endif
2368
2369         enddo
2370
2371         enddo ! iint
2372       enddo ! i
2373       do i=1,nct
2374         do j=1,3
2375           gvdwc_scp(j,i)=expon*gvdwc_scp(j,i)
2376           gradx_scp(j,i)=expon*gradx_scp(j,i)
2377         enddo
2378       enddo
2379 C******************************************************************************
2380 C
2381 C                              N O T E !!!
2382 C
2383 C To save time the factor EXPON has been extracted from ALL components
2384 C of GVDWC and GRADX. Remember to multiply them by this factor before further
2385 C use!
2386 C
2387 C******************************************************************************
2388       return
2389       end
2390 C-----------------------------------------------------------------------------
2391       subroutine escp_short(evdw2,evdw2_14)
2392 C
2393 C This subroutine calculates the excluded-volume interaction energy between
2394 C peptide-group centers and side chains and its gradient in virtual-bond and
2395 C side-chain vectors.
2396 C
2397       implicit real*8 (a-h,o-z)
2398       include 'DIMENSIONS'
2399       include 'COMMON.GEO'
2400       include 'COMMON.VAR'
2401       include 'COMMON.LOCAL'
2402       include 'COMMON.CHAIN'
2403       include 'COMMON.DERIV'
2404       include 'COMMON.INTERACT'
2405       include 'COMMON.FFIELD'
2406       include 'COMMON.IOUNITS'
2407       include 'COMMON.CONTROL'
2408       dimension ggg(3)
2409       evdw2=0.0D0
2410       evdw2_14=0.0d0
2411 cd    print '(a)','Enter ESCP'
2412 cd    write (iout,*) 'iatscp_s=',iatscp_s,' iatscp_e=',iatscp_e
2413       do i=iatscp_s,iatscp_e
2414         iteli=itel(i)
2415         xi=0.5D0*(c(1,i)+c(1,i+1))
2416         yi=0.5D0*(c(2,i)+c(2,i+1))
2417         zi=0.5D0*(c(3,i)+c(3,i+1))
2418
2419         do iint=1,nscp_gr(i)
2420
2421         do j=iscpstart(i,iint),iscpend(i,iint)
2422           itypj=itype(j)
2423 C Uncomment following three lines for SC-p interactions
2424 c         xj=c(1,nres+j)-xi
2425 c         yj=c(2,nres+j)-yi
2426 c         zj=c(3,nres+j)-zi
2427 C Uncomment following three lines for Ca-p interactions
2428           xj=c(1,j)-xi
2429           yj=c(2,j)-yi
2430           zj=c(3,j)-zi
2431           rrij=1.0D0/(xj*xj+yj*yj+zj*zj)
2432
2433           sss=sscale(1.0d0/(dsqrt(rrij)*rscp(itypj,iteli)))
2434
2435           if (sss.gt.0.0d0) then
2436
2437             fac=rrij**expon2
2438             e1=fac*fac*aad(itypj,iteli)
2439             e2=fac*bad(itypj,iteli)
2440             if (iabs(j-i) .le. 2) then
2441               e1=scal14*e1
2442               e2=scal14*e2
2443               evdw2_14=evdw2_14+(e1+e2)*sss
2444             endif
2445             evdwij=e1+e2
2446             evdw2=evdw2+evdwij*sss
2447             if (energy_dec) write (iout,'(a6,2i5,0pf7.3)')
2448      &          'evdw2',i,j,evdwij
2449 C
2450 C Calculate contributions to the gradient in the virtual-bond and SC vectors.
2451 C
2452             fac=-(evdwij+e1)*rrij*sss
2453             ggg(1)=xj*fac
2454             ggg(2)=yj*fac
2455             ggg(3)=zj*fac
2456             if (j.lt.i) then
2457 cd            write (iout,*) 'j<i'
2458 C Uncomment following three lines for SC-p interactions
2459 c             do k=1,3
2460 c             gradx_scp(k,j)=gradx_scp(k,j)+ggg(k)
2461 c             enddo
2462             else
2463 cd            write (iout,*) 'j>i'
2464               do k=1,3
2465                 ggg(k)=-ggg(k)
2466 C Uncomment following line for SC-p interactions
2467 c               gradx_scp(k,j)=gradx_scp(k,j)-ggg(k)
2468               enddo
2469             endif
2470             do k=1,3
2471               gvdwc_scp(k,i)=gvdwc_scp(k,i)-0.5D0*ggg(k)
2472             enddo
2473             kstart=min0(i+1,j)
2474             kend=max0(i-1,j-1)
2475 cd        write (iout,*) 'i=',i,' j=',j,' kstart=',kstart,' kend=',kend
2476 cd        write (iout,*) ggg(1),ggg(2),ggg(3)
2477             do k=kstart,kend
2478               do l=1,3
2479                 gvdwc_scp(l,k)=gvdwc_scp(l,k)-ggg(l)
2480               enddo
2481             enddo
2482
2483           endif
2484
2485         enddo
2486
2487         enddo ! iint
2488       enddo ! i
2489       do i=1,nct
2490         do j=1,3
2491           gvdwc_scp(j,i)=expon*gvdwc_scp(j,i)
2492           gradx_scp(j,i)=expon*gradx_scp(j,i)
2493         enddo
2494       enddo
2495 C******************************************************************************
2496 C
2497 C                              N O T E !!!
2498 C
2499 C To save time the factor EXPON has been extracted from ALL components
2500 C of GVDWC and GRADX. Remember to multiply them by this factor before further
2501 C use!
2502 C
2503 C******************************************************************************
2504       return
2505       end