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

   1 C-----------------------------------------------------------------------
   2       double precision function sscale(r)
   3       double precision r,gamm
   4       include "COMMON.SPLITELE"
   5       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         num_conti=0
2132 c        write (iout,*) 'i',i,' ielstart',ielstart(i),' ielend',ielend(i)
2133         do j=ielstart(i),ielend(i)
2134           ind=ind+1
2135           iteli=itel(i)
2136           itelj=itel(j)
2137           if (j.eq.i+2 .and. itelj.eq.2) iteli=2
2138           aaa=app(iteli,itelj)
2139           bbb=bpp(iteli,itelj)
2140           dxj=dc(1,j)
2141           dyj=dc(2,j)
2142           dzj=dc(3,j)
2143           dx_normj=dc_norm(1,j)
2144           dy_normj=dc_norm(2,j)
2145           dz_normj=dc_norm(3,j)
2146           xj=c(1,j)+0.5D0*dxj-xmedi
2147           yj=c(2,j)+0.5D0*dyj-ymedi
2148           zj=c(3,j)+0.5D0*dzj-zmedi
2149           rij=xj*xj+yj*yj+zj*zj
2150           rrmij=1.0D0/rij
2151           rij=dsqrt(rij)
2152           sss=sscale(rij/rpp(iteli,itelj))
2153           if (sss.lt.1.0d0) then
2154             rmij=1.0D0/rij
2155             r3ij=rrmij*rmij
2156             r6ij=r3ij*r3ij
2157             ev1=aaa*r6ij*r6ij
2158 c 4/26/02 - AL scaling down 1,4 repulsive VDW interactions
2159             if (j.eq.i+2) ev1=scal_el*ev1
2160             ev2=bbb*r6ij
2161             evdwij=ev1+ev2
2162             if (energy_dec) then
2163               write (iout,'(a6,2i5,0pf7.3)') 'evdw1',i,j,evdwij
2164               write (iout,'(a6,2i5,0pf7.3)') 'ees',i,j,eesij
2165             endif
2166             evdw1=evdw1+evdwij*(1.0d0-sss)
2167 C
2168 C Calculate contributions to the Cartesian gradient.
2169 C
2170             facvdw=-6*rrmij*(ev1+evdwij)*(1.0d0-sss)
2171             ggg(1)=facvdw*xj
2172             ggg(2)=facvdw*yj
2173             ggg(3)=facvdw*zj
2174
2175             do k=1,3
2176               ghalf=0.5D0*ggg(k)
2177               gvdwpp(k,i)=gvdwpp(k,i)+ghalf
2178               gvdwpp(k,j)=gvdwpp(k,j)+ghalf
2179             enddo
2180 *
2181 * Loop over residues i+1 thru j-1.
2182 *
2183             do k=i+1,j-1
2184               do l=1,3
2185                 gvdwpp(l,k)=gvdwpp(l,k)+ggg(l)
2186               enddo
2187             enddo
2188           endif
2189         enddo ! j
2190       enddo   ! i
2191       return
2192       end
2193 C-----------------------------------------------------------------------
2194       subroutine evdwpp_short(evdw1)
2195 C
2196 C Compute Evdwpp
2197 C
2198       implicit real*8 (a-h,o-z)
2199       include 'DIMENSIONS'
2200       include 'COMMON.CONTROL'
2201       include 'COMMON.IOUNITS'
2202       include 'COMMON.GEO'
2203       include 'COMMON.VAR'
2204       include 'COMMON.LOCAL'
2205       include 'COMMON.CHAIN'
2206       include 'COMMON.DERIV'
2207       include 'COMMON.INTERACT'
2208       include 'COMMON.CONTACTS'
2209       include 'COMMON.TORSION'
2210       include 'COMMON.VECTORS'
2211       include 'COMMON.FFIELD'
2212       dimension ggg(3)
2213 c 4/26/02 - AL scaling factor for 1,4 repulsive VDW interactions
2214 #ifdef MOMENT
2215       double precision scal_el /1.0d0/
2216 #else
2217       double precision scal_el /0.5d0/
2218 #endif
2219       evdw1=0.0D0
2220       do i=iatel_s,iatel_e
2221         dxi=dc(1,i)
2222         dyi=dc(2,i)
2223         dzi=dc(3,i)
2224         dx_normi=dc_norm(1,i)
2225         dy_normi=dc_norm(2,i)
2226         dz_normi=dc_norm(3,i)
2227         xmedi=c(1,i)+0.5d0*dxi
2228         ymedi=c(2,i)+0.5d0*dyi
2229         zmedi=c(3,i)+0.5d0*dzi
2230         num_conti=0
2231 c        write (iout,*) 'i',i,' ielstart',ielstart(i),' ielend',ielend(i)
2232         do j=ielstart(i),ielend(i)
2233           ind=ind+1
2234           iteli=itel(i)
2235           itelj=itel(j)
2236           if (j.eq.i+2 .and. itelj.eq.2) iteli=2
2237           aaa=app(iteli,itelj)
2238           bbb=bpp(iteli,itelj)
2239           dxj=dc(1,j)
2240           dyj=dc(2,j)
2241           dzj=dc(3,j)
2242           dx_normj=dc_norm(1,j)
2243           dy_normj=dc_norm(2,j)
2244           dz_normj=dc_norm(3,j)
2245           xj=c(1,j)+0.5D0*dxj-xmedi
2246           yj=c(2,j)+0.5D0*dyj-ymedi
2247           zj=c(3,j)+0.5D0*dzj-zmedi
2248           rij=xj*xj+yj*yj+zj*zj
2249           rrmij=1.0D0/rij
2250           rij=dsqrt(rij)
2251           sss=sscale(rij/rpp(iteli,itelj))
2252           if (sss.gt.0.0d0) then
2253             rmij=1.0D0/rij
2254             r3ij=rrmij*rmij
2255             r6ij=r3ij*r3ij
2256             ev1=aaa*r6ij*r6ij
2257 c 4/26/02 - AL scaling down 1,4 repulsive VDW interactions
2258             if (j.eq.i+2) ev1=scal_el*ev1
2259             ev2=bbb*r6ij
2260             evdwij=ev1+ev2
2261             if (energy_dec) then
2262               write (iout,'(a6,2i5,0pf7.3)') 'evdw1',i,j,evdwij
2263               write (iout,'(a6,2i5,0pf7.3)') 'ees',i,j,eesij
2264             endif
2265             evdw1=evdw1+evdwij*sss
2266 C
2267 C Calculate contributions to the Cartesian gradient.
2268 C
2269             facvdw=-6*rrmij*(ev1+evdwij)*sss
2270             ggg(1)=facvdw*xj
2271             ggg(2)=facvdw*yj
2272             ggg(3)=facvdw*zj
2273
2274             do k=1,3
2275               ghalf=0.5D0*ggg(k)
2276               gvdwpp(k,i)=gvdwpp(k,i)+ghalf
2277               gvdwpp(k,j)=gvdwpp(k,j)+ghalf
2278             enddo
2279 *
2280 * Loop over residues i+1 thru j-1.
2281 *
2282             do k=i+1,j-1
2283               do l=1,3
2284                 gvdwpp(l,k)=gvdwpp(l,k)+ggg(l)
2285               enddo
2286             enddo
2287           endif
2288         enddo ! j
2289       enddo   ! i
2290       return
2291       end
2292 C-----------------------------------------------------------------------------
2293       subroutine escp_long(evdw2,evdw2_14)
2294 C
2295 C This subroutine calculates the excluded-volume interaction energy between
2296 C peptide-group centers and side chains and its gradient in virtual-bond and
2297 C side-chain vectors.
2298 C
2299       implicit real*8 (a-h,o-z)
2300       include 'DIMENSIONS'
2301       include 'COMMON.GEO'
2302       include 'COMMON.VAR'
2303       include 'COMMON.LOCAL'
2304       include 'COMMON.CHAIN'
2305       include 'COMMON.DERIV'
2306       include 'COMMON.INTERACT'
2307       include 'COMMON.FFIELD'
2308       include 'COMMON.IOUNITS'
2309       include 'COMMON.CONTROL'
2310       dimension ggg(3)
2311       evdw2=0.0D0
2312       evdw2_14=0.0d0
2313 cd    print '(a)','Enter ESCP'
2314 cd    write (iout,*) 'iatscp_s=',iatscp_s,' iatscp_e=',iatscp_e
2315       do i=iatscp_s,iatscp_e
2316         iteli=itel(i)
2317         xi=0.5D0*(c(1,i)+c(1,i+1))
2318         yi=0.5D0*(c(2,i)+c(2,i+1))
2319         zi=0.5D0*(c(3,i)+c(3,i+1))
2320
2321         do iint=1,nscp_gr(i)
2322
2323         do j=iscpstart(i,iint),iscpend(i,iint)
2324           itypj=itype(j)
2325 C Uncomment following three lines for SC-p interactions
2326 c         xj=c(1,nres+j)-xi
2327 c         yj=c(2,nres+j)-yi
2328 c         zj=c(3,nres+j)-zi
2329 C Uncomment following three lines for Ca-p interactions
2330           xj=c(1,j)-xi
2331           yj=c(2,j)-yi
2332           zj=c(3,j)-zi
2333           rrij=1.0D0/(xj*xj+yj*yj+zj*zj)
2334
2335           sss=sscale(1.0d0/(dsqrt(rrij)*rscp(itypj,iteli)))
2336
2337           if (sss.lt.1.0d0) then
2338
2339             fac=rrij**expon2
2340             e1=fac*fac*aad(itypj,iteli)
2341             e2=fac*bad(itypj,iteli)
2342             if (iabs(j-i) .le. 2) then
2343               e1=scal14*e1
2344               e2=scal14*e2
2345               evdw2_14=evdw2_14+(e1+e2)*(1.0d0-sss)
2346             endif
2347             evdwij=e1+e2
2348             evdw2=evdw2+evdwij*(1.0d0-sss)
2349             if (energy_dec) write (iout,'(a6,2i5,0pf7.3)')
2350      &          'evdw2',i,j,evdwij
2351 C
2352 C Calculate contributions to the gradient in the virtual-bond and SC vectors.
2353 C
2354             fac=-(evdwij+e1)*rrij*(1.0d0-sss)
2355             ggg(1)=xj*fac
2356             ggg(2)=yj*fac
2357             ggg(3)=zj*fac
2358             if (j.lt.i) then
2359 cd            write (iout,*) 'j<i'
2360 C Uncomment following three lines for SC-p interactions
2361 c             do k=1,3
2362 c             gradx_scp(k,j)=gradx_scp(k,j)+ggg(k)
2363 c             enddo
2364             else
2365 cd            write (iout,*) 'j>i'
2366               do k=1,3
2367                 ggg(k)=-ggg(k)
2368 C Uncomment following line for SC-p interactions
2369 c               gradx_scp(k,j)=gradx_scp(k,j)-ggg(k)
2370               enddo
2371             endif
2372             do k=1,3
2373               gvdwc_scp(k,i)=gvdwc_scp(k,i)-0.5D0*ggg(k)
2374             enddo
2375             kstart=min0(i+1,j)
2376             kend=max0(i-1,j-1)
2377 cd        write (iout,*) 'i=',i,' j=',j,' kstart=',kstart,' kend=',kend
2378 cd        write (iout,*) ggg(1),ggg(2),ggg(3)
2379             do k=kstart,kend
2380               do l=1,3
2381                 gvdwc_scp(l,k)=gvdwc_scp(l,k)-ggg(l)
2382               enddo
2383             enddo
2384
2385           endif
2386
2387         enddo
2388
2389         enddo ! iint
2390       enddo ! i
2391       do i=1,nct
2392         do j=1,3
2393           gvdwc_scp(j,i)=expon*gvdwc_scp(j,i)
2394           gradx_scp(j,i)=expon*gradx_scp(j,i)
2395         enddo
2396       enddo
2397 C******************************************************************************
2398 C
2399 C                              N O T E !!!
2400 C
2401 C To save time the factor EXPON has been extracted from ALL components
2402 C of GVDWC and GRADX. Remember to multiply them by this factor before further
2403 C use!
2404 C
2405 C******************************************************************************
2406       return
2407       end
2408 C-----------------------------------------------------------------------------
2409       subroutine escp_short(evdw2,evdw2_14)
2410 C
2411 C This subroutine calculates the excluded-volume interaction energy between
2412 C peptide-group centers and side chains and its gradient in virtual-bond and
2413 C side-chain vectors.
2414 C
2415       implicit real*8 (a-h,o-z)
2416       include 'DIMENSIONS'
2417       include 'COMMON.GEO'
2418       include 'COMMON.VAR'
2419       include 'COMMON.LOCAL'
2420       include 'COMMON.CHAIN'
2421       include 'COMMON.DERIV'
2422       include 'COMMON.INTERACT'
2423       include 'COMMON.FFIELD'
2424       include 'COMMON.IOUNITS'
2425       include 'COMMON.CONTROL'
2426       dimension ggg(3)
2427       evdw2=0.0D0
2428       evdw2_14=0.0d0
2429 cd    print '(a)','Enter ESCP'
2430 cd    write (iout,*) 'iatscp_s=',iatscp_s,' iatscp_e=',iatscp_e
2431       do i=iatscp_s,iatscp_e
2432         iteli=itel(i)
2433         xi=0.5D0*(c(1,i)+c(1,i+1))
2434         yi=0.5D0*(c(2,i)+c(2,i+1))
2435         zi=0.5D0*(c(3,i)+c(3,i+1))
2436
2437         do iint=1,nscp_gr(i)
2438
2439         do j=iscpstart(i,iint),iscpend(i,iint)
2440           itypj=itype(j)
2441 C Uncomment following three lines for SC-p interactions
2442 c         xj=c(1,nres+j)-xi
2443 c         yj=c(2,nres+j)-yi
2444 c         zj=c(3,nres+j)-zi
2445 C Uncomment following three lines for Ca-p interactions
2446           xj=c(1,j)-xi
2447           yj=c(2,j)-yi
2448           zj=c(3,j)-zi
2449           rrij=1.0D0/(xj*xj+yj*yj+zj*zj)
2450
2451           sss=sscale(1.0d0/(dsqrt(rrij)*rscp(itypj,iteli)))
2452
2453           if (sss.gt.0.0d0) then
2454
2455             fac=rrij**expon2
2456             e1=fac*fac*aad(itypj,iteli)
2457             e2=fac*bad(itypj,iteli)
2458             if (iabs(j-i) .le. 2) then
2459               e1=scal14*e1
2460               e2=scal14*e2
2461               evdw2_14=evdw2_14+(e1+e2)*sss
2462             endif
2463             evdwij=e1+e2
2464             evdw2=evdw2+evdwij*sss
2465             if (energy_dec) write (iout,'(a6,2i5,0pf7.3)')
2466      &          'evdw2',i,j,evdwij
2467 C
2468 C Calculate contributions to the gradient in the virtual-bond and SC vectors.
2469 C
2470             fac=-(evdwij+e1)*rrij*sss
2471             ggg(1)=xj*fac
2472             ggg(2)=yj*fac
2473             ggg(3)=zj*fac
2474             if (j.lt.i) then
2475 cd            write (iout,*) 'j<i'
2476 C Uncomment following three lines for SC-p interactions
2477 c             do k=1,3
2478 c             gradx_scp(k,j)=gradx_scp(k,j)+ggg(k)
2479 c             enddo
2480             else
2481 cd            write (iout,*) 'j>i'
2482               do k=1,3
2483                 ggg(k)=-ggg(k)
2484 C Uncomment following line for SC-p interactions
2485 c               gradx_scp(k,j)=gradx_scp(k,j)-ggg(k)
2486               enddo
2487             endif
2488             do k=1,3
2489               gvdwc_scp(k,i)=gvdwc_scp(k,i)-0.5D0*ggg(k)
2490             enddo
2491             kstart=min0(i+1,j)
2492             kend=max0(i-1,j-1)
2493 cd        write (iout,*) 'i=',i,' j=',j,' kstart=',kstart,' kend=',kend
2494 cd        write (iout,*) ggg(1),ggg(2),ggg(3)
2495             do k=kstart,kend
2496               do l=1,3
2497                 gvdwc_scp(l,k)=gvdwc_scp(l,k)-ggg(l)
2498               enddo
2499             enddo
2500
2501           endif
2502
2503         enddo
2504
2505         enddo ! iint
2506       enddo ! i
2507       do i=1,nct
2508         do j=1,3
2509           gvdwc_scp(j,i)=expon*gvdwc_scp(j,i)
2510           gradx_scp(j,i)=expon*gradx_scp(j,i)
2511         enddo
2512       enddo
2513 C******************************************************************************
2514 C
2515 C                              N O T E !!!
2516 C
2517 C To save time the factor EXPON has been extracted from ALL components
2518 C of GVDWC and GRADX. Remember to multiply them by this factor before further
2519 C use!
2520 C
2521 C******************************************************************************
2522       return
2523       end