energy_p_new-sep.F

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