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

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