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           xi=mod(xi,boxxsize)
 631           if (xi.lt.0) xi=xi+boxxsize
 632           yi=mod(yi,boxysize)
 633           if (yi.lt.0) yi=yi+boxysize
 634           zi=mod(zi,boxzsize)
 635           if (zi.lt.0) zi=zi+boxzsize
 636         dxi=dc_norm(1,nres+i)
 637         dyi=dc_norm(2,nres+i)
 638         dzi=dc_norm(3,nres+i)
 639 c        dsci_inv=dsc_inv(itypi)
 640         dsci_inv=vbld_inv(i+nres)
 641 c        write (iout,*) "i",i,dsc_inv(itypi),dsci_inv,1.0d0/vbld(i+nres)
 642 c        write (iout,*) "dcnori",dxi*dxi+dyi*dyi+dzi*dzi
 643 C
 644 C Calculate SC interaction energy.
 645 C
 646         do iint=1,nint_gr(i)
 647           do j=istart(i,iint),iend(i,iint)
 648             ind=ind+1
 649             itypj=itype(j)
 650             if (itypj.eq.ntyp1) cycle
 651 c            dscj_inv=dsc_inv(itypj)
 652             dscj_inv=vbld_inv(j+nres)
 653 c            write (iout,*) "j",j,dsc_inv(itypj),dscj_inv,
 654 c     &       1.0d0/vbld(j+nres)
 655 c            write (iout,*) "i",i," j", j," itype",itype(i),itype(j)
 656             sig0ij=sigma(itypi,itypj)
 657             chi1=chi(itypi,itypj)
 658             chi2=chi(itypj,itypi)
 659             chi12=chi1*chi2
 660             chip1=chip(itypi)
 661             chip2=chip(itypj)
 662             chip12=chip1*chip2
 663             alf1=alp(itypi)
 664             alf2=alp(itypj)
 665             alf12=0.5D0*(alf1+alf2)
 666             xj=c(1,nres+j)
 667             yj=c(2,nres+j)
 668             zj=c(3,nres+j)
 669           xj=mod(xj,boxxsize)
 670           if (xj.lt.0) xj=xj+boxxsize
 671           yj=mod(yj,boxysize)
 672           if (yj.lt.0) yj=yj+boxysize
 673           zj=mod(zj,boxzsize)
 674           if (zj.lt.0) zj=zj+boxzsize
 675       dist_init=(xj-xi)**2+(yj-yi)**2+(zj-zi)**2
 676       xj_safe=xj
 677       yj_safe=yj
 678       zj_safe=zj
 679       subchap=0
 680       do xshift=-1,1
 681       do yshift=-1,1
 682       do zshift=-1,1
 683           xj=xj_safe+xshift*boxxsize
 684           yj=yj_safe+yshift*boxysize
 685           zj=zj_safe+zshift*boxzsize
 686           dist_temp=(xj-xi)**2+(yj-yi)**2+(zj-zi)**2
 687           if(dist_temp.lt.dist_init) then
 688             dist_init=dist_temp
 689             xj_temp=xj
 690             yj_temp=yj
 691             zj_temp=zj
 692             subchap=1
 693           endif
 694        enddo
 695        enddo
 696        enddo
 697        if (subchap.eq.1) then
 698           xj=xj_temp-xi
 699           yj=yj_temp-yi
 700           zj=zj_temp-zi
 701        else
 702           xj=xj_safe-xi
 703           yj=yj_safe-yi
 704           zj=zj_safe-zi
 705        endif
 706             dxj=dc_norm(1,nres+j)
 707             dyj=dc_norm(2,nres+j)
 708             dzj=dc_norm(3,nres+j)
 709             rrij=1.0D0/(xj*xj+yj*yj+zj*zj)
 710             rij=dsqrt(rrij)
 711             sss=sscale(1.0d0/(rij*sigmaii(itypi,itypj)))
 712             sssgrad=sscagrad((1.0d0/rij)/sigmaii(itypi,itypj))
 713             if (sss.lt.1.0d0) then
 714
 715 C Calculate angle-dependent terms of energy and contributions to their
 716 C derivatives.
 717               call sc_angular
 718               sigsq=1.0D0/sigsq
 719               sig=sig0ij*dsqrt(sigsq)
 720               rij_shift=1.0D0/rij-sig+sig0ij
 721 c for diagnostics; uncomment
 722 c              rij_shift=1.2*sig0ij
 723 C I hate to put IF's in the loops, but here don't have another choice!!!!
 724               if (rij_shift.le.0.0D0) then
 725                 evdw=1.0D20
 726 cd                write (iout,'(2(a3,i3,2x),17(0pf7.3))')
 727 cd     &          restyp(itypi),i,restyp(itypj),j,
 728 cd     &          rij_shift,1.0D0/rij,sig,sig0ij,sigsq,1-dsqrt(sigsq)
 729                 return
 730               endif
 731               sigder=-sig*sigsq
 732 c---------------------------------------------------------------
 733               rij_shift=1.0D0/rij_shift
 734               fac=rij_shift**expon
 735               e1=fac*fac*aa(itypi,itypj)
 736               e2=fac*bb(itypi,itypj)
 737               evdwij=eps1*eps2rt*eps3rt*(e1+e2)
 738               eps2der=evdwij*eps3rt
 739               eps3der=evdwij*eps2rt
 740 c              write (iout,*) "sigsq",sigsq," sig",sig," eps2rt",eps2rt,
 741 c     &        " eps3rt",eps3rt," eps1",eps1," e1",e1," e2",e2
 742               evdwij=evdwij*eps2rt*eps3rt
 743               evdw=evdw+evdwij*(1.0d0-sss)
 744               if (lprn) then
 745               sigm=dabs(aa(itypi,itypj)/bb(itypi,itypj))**(1.0D0/6.0D0)
 746               epsi=bb(itypi,itypj)**2/aa(itypi,itypj)
 747               write (iout,'(2(a3,i3,2x),17(0pf7.3))')
 748      &          restyp(itypi),i,restyp(itypj),j,
 749      &          epsi,sigm,chi1,chi2,chip1,chip2,
 750      &          eps1,eps2rt**2,eps3rt**2,sig,sig0ij,
 751      &          om1,om2,om12,1.0D0/rij,1.0D0/rij_shift,
 752      &          evdwij
 753               endif
 754
 755               if (energy_dec) write (iout,'(a6,2i5,0pf7.3)')
 756      &                        'evdw',i,j,evdwij
 757
 758 C Calculate gradient components.
 759               e1=e1*eps1*eps2rt**2*eps3rt**2
 760               fac=-expon*(e1+evdwij)*rij_shift
 761               sigder=fac*sigder
 762               fac=rij*fac
 763             fac=fac+evdwij/(1.0-sss)*(-sssgrad)/sigmaii(itypi,itypj)*rij
 764 c              fac=0.0d0
 765 C Calculate the radial part of the gradient
 766               gg(1)=xj*fac
 767               gg(2)=yj*fac
 768               gg(3)=zj*fac
 769 C Calculate angular part of the gradient.
 770               call sc_grad_scale(1.0d0-sss)
 771             endif
 772           enddo      ! j
 773         enddo        ! iint
 774       enddo          ! i
 775 c      write (iout,*) "Number of loop steps in EGB:",ind
 776 cccc      energy_dec=.false.
 777       return
 778       end
 779 C-----------------------------------------------------------------------------
 780       subroutine egb_short(evdw)
 781 C
 782 C This subroutine calculates the interaction energy of nonbonded side chains
 783 C assuming the Gay-Berne potential of interaction.
 784 C
 785       implicit real*8 (a-h,o-z)
 786       include 'DIMENSIONS'
 787       include 'COMMON.GEO'
 788       include 'COMMON.VAR'
 789       include 'COMMON.LOCAL'
 790       include 'COMMON.CHAIN'
 791       include 'COMMON.DERIV'
 792       include 'COMMON.NAMES'
 793       include 'COMMON.INTERACT'
 794       include 'COMMON.IOUNITS'
 795       include 'COMMON.CALC'
 796       include 'COMMON.CONTROL'
 797       logical lprn
 798       evdw=0.0D0
 799 ccccc      energy_dec=.false.
 800 c     print *,'Entering EGB nnt=',nnt,' nct=',nct,' expon=',expon
 801       evdw=0.0D0
 802       lprn=.false.
 803 c     if (icall.eq.0) lprn=.false.
 804       ind=0
 805       do i=iatsc_s,iatsc_e
 806         itypi=itype(i)
 807         if (itypi.eq.ntyp1) cycle
 808         itypi1=itype(i+1)
 809         xi=c(1,nres+i)
 810         yi=c(2,nres+i)
 811         zi=c(3,nres+i)
 812           xi=mod(xi,boxxsize)
 813           if (xi.lt.0) xi=xi+boxxsize
 814           yi=mod(yi,boxysize)
 815           if (yi.lt.0) yi=yi+boxysize
 816           zi=mod(zi,boxzsize)
 817           if (zi.lt.0) zi=zi+boxzsize
 818         dxi=dc_norm(1,nres+i)
 819         dyi=dc_norm(2,nres+i)
 820         dzi=dc_norm(3,nres+i)
 821 c        dsci_inv=dsc_inv(itypi)
 822         dsci_inv=vbld_inv(i+nres)
 823 c        write (iout,*) "i",i,dsc_inv(itypi),dsci_inv,1.0d0/vbld(i+nres)
 824 c        write (iout,*) "dcnori",dxi*dxi+dyi*dyi+dzi*dzi
 825 C
 826 C Calculate SC interaction energy.
 827 C
 828         do iint=1,nint_gr(i)
 829           do j=istart(i,iint),iend(i,iint)
 830             ind=ind+1
 831             itypj=itype(j)
 832             if (itypj.eq.ntyp1) cycle
 833 c            dscj_inv=dsc_inv(itypj)
 834             dscj_inv=vbld_inv(j+nres)
 835 c            write (iout,*) "j",j,dsc_inv(itypj),dscj_inv,
 836 c     &       1.0d0/vbld(j+nres)
 837 c            write (iout,*) "i",i," j", j," itype",itype(i),itype(j)
 838             sig0ij=sigma(itypi,itypj)
 839             chi1=chi(itypi,itypj)
 840             chi2=chi(itypj,itypi)
 841             chi12=chi1*chi2
 842             chip1=chip(itypi)
 843             chip2=chip(itypj)
 844             chip12=chip1*chip2
 845             alf1=alp(itypi)
 846             alf2=alp(itypj)
 847             alf12=0.5D0*(alf1+alf2)
 848             xj=c(1,nres+j)
 849             yj=c(2,nres+j)
 850             zj=c(3,nres+j)
 851           xj=mod(xj,boxxsize)
 852           if (xj.lt.0) xj=xj+boxxsize
 853           yj=mod(yj,boxysize)
 854           if (yj.lt.0) yj=yj+boxysize
 855           zj=mod(zj,boxzsize)
 856           if (zj.lt.0) zj=zj+boxzsize
 857       dist_init=(xj-xi)**2+(yj-yi)**2+(zj-zi)**2
 858       xj_safe=xj
 859       yj_safe=yj
 860       zj_safe=zj
 861       subchap=0
 862       do xshift=-1,1
 863       do yshift=-1,1
 864       do zshift=-1,1
 865           xj=xj_safe+xshift*boxxsize
 866           yj=yj_safe+yshift*boxysize
 867           zj=zj_safe+zshift*boxzsize
 868           dist_temp=(xj-xi)**2+(yj-yi)**2+(zj-zi)**2
 869           if(dist_temp.lt.dist_init) then
 870             dist_init=dist_temp
 871             xj_temp=xj
 872             yj_temp=yj
 873             zj_temp=zj
 874             subchap=1
 875           endif
 876        enddo
 877        enddo
 878        enddo
 879        if (subchap.eq.1) then
 880           xj=xj_temp-xi
 881           yj=yj_temp-yi
 882           zj=zj_temp-zi
 883        else
 884           xj=xj_safe-xi
 885           yj=yj_safe-yi
 886           zj=zj_safe-zi
 887        endif
 888             dxj=dc_norm(1,nres+j)
 889             dyj=dc_norm(2,nres+j)
 890             dzj=dc_norm(3,nres+j)
 891             rrij=1.0D0/(xj*xj+yj*yj+zj*zj)
 892             rij=dsqrt(rrij)
 893             sss=sscale(1.0d0/(rij*sigmaii(itypi,itypj)))
 894             sssgrad=sscagrad((1.0d0/rij)/sigmaii(itypi,itypj))
 895             if (sss.gt.0.0d0) then
 896
 897 C Calculate angle-dependent terms of energy and contributions to their
 898 C derivatives.
 899               call sc_angular
 900               sigsq=1.0D0/sigsq
 901               sig=sig0ij*dsqrt(sigsq)
 902               rij_shift=1.0D0/rij-sig+sig0ij
 903 c for diagnostics; uncomment
 904 c              rij_shift=1.2*sig0ij
 905 C I hate to put IF's in the loops, but here don't have another choice!!!!
 906               if (rij_shift.le.0.0D0) then
 907                 evdw=1.0D20
 908 cd                write (iout,'(2(a3,i3,2x),17(0pf7.3))')
 909 cd     &          restyp(itypi),i,restyp(itypj),j,
 910 cd     &          rij_shift,1.0D0/rij,sig,sig0ij,sigsq,1-dsqrt(sigsq)
 911                 return
 912               endif
 913               sigder=-sig*sigsq
 914 c---------------------------------------------------------------
 915               rij_shift=1.0D0/rij_shift
 916               fac=rij_shift**expon
 917               e1=fac*fac*aa(itypi,itypj)
 918               e2=fac*bb(itypi,itypj)
 919               evdwij=eps1*eps2rt*eps3rt*(e1+e2)
 920               eps2der=evdwij*eps3rt
 921               eps3der=evdwij*eps2rt
 922 c              write (iout,*) "sigsq",sigsq," sig",sig," eps2rt",eps2rt,
 923 c     &        " eps3rt",eps3rt," eps1",eps1," e1",e1," e2",e2
 924               evdwij=evdwij*eps2rt*eps3rt
 925               evdw=evdw+evdwij*sss
 926               if (lprn) then
 927               sigm=dabs(aa(itypi,itypj)/bb(itypi,itypj))**(1.0D0/6.0D0)
 928               epsi=bb(itypi,itypj)**2/aa(itypi,itypj)
 929               write (iout,'(2(a3,i3,2x),17(0pf7.3))')
 930      &          restyp(itypi),i,restyp(itypj),j,
 931      &          epsi,sigm,chi1,chi2,chip1,chip2,
 932      &          eps1,eps2rt**2,eps3rt**2,sig,sig0ij,
 933      &          om1,om2,om12,1.0D0/rij,1.0D0/rij_shift,
 934      &          evdwij
 935               endif
 936
 937               if (energy_dec) write (iout,'(a6,2i5,0pf7.3)')
 938      &                        'evdw',i,j,evdwij
 939
 940 C Calculate gradient components.
 941               e1=e1*eps1*eps2rt**2*eps3rt**2
 942               fac=-expon*(e1+evdwij)*rij_shift
 943               sigder=fac*sigder
 944               fac=rij*fac
 945             fac=fac+evdwij/sss*sssgrad/sigmaii(itypi,itypj)*rij
 946 c              fac=0.0d0
 947 C Calculate the radial part of the gradient
 948               gg(1)=xj*fac
 949               gg(2)=yj*fac
 950               gg(3)=zj*fac
 951 C Calculate angular part of the gradient.
 952               call sc_grad_scale(sss)
 953             endif
 954           enddo      ! j
 955         enddo        ! iint
 956       enddo          ! i
 957 c      write (iout,*) "Number of loop steps in EGB:",ind
 958 cccc      energy_dec=.false.
 959       return
 960       end
 961 C-----------------------------------------------------------------------------
 962       subroutine egbv_long(evdw)
 963 C
 964 C This subroutine calculates the interaction energy of nonbonded side chains
 965 C assuming the Gay-Berne-Vorobjev potential of interaction.
 966 C
 967       implicit real*8 (a-h,o-z)
 968       include 'DIMENSIONS'
 969       include 'COMMON.GEO'
 970       include 'COMMON.VAR'
 971       include 'COMMON.LOCAL'
 972       include 'COMMON.CHAIN'
 973       include 'COMMON.DERIV'
 974       include 'COMMON.NAMES'
 975       include 'COMMON.INTERACT'
 976       include 'COMMON.IOUNITS'
 977       include 'COMMON.CALC'
 978       common /srutu/ icall
 979       logical lprn
 980       evdw=0.0D0
 981 c     print *,'Entering EGB nnt=',nnt,' nct=',nct,' expon=',expon
 982       evdw=0.0D0
 983       lprn=.false.
 984 c     if (icall.eq.0) lprn=.true.
 985       ind=0
 986       do i=iatsc_s,iatsc_e
 987         itypi=itype(i)
 988         if (itypi.eq.ntyp1) cycle
 989         itypi1=itype(i+1)
 990         xi=c(1,nres+i)
 991         yi=c(2,nres+i)
 992         zi=c(3,nres+i)
 993         dxi=dc_norm(1,nres+i)
 994         dyi=dc_norm(2,nres+i)
 995         dzi=dc_norm(3,nres+i)
 996 c        dsci_inv=dsc_inv(itypi)
 997         dsci_inv=vbld_inv(i+nres)
 998 C
 999 C Calculate SC interaction energy.
1000 C
1001         do iint=1,nint_gr(i)
1002           do j=istart(i,iint),iend(i,iint)
1003             ind=ind+1
1004             itypj=itype(j)
1005             if (itypj.eq.ntyp1) cycle
1006 c            dscj_inv=dsc_inv(itypj)
1007             dscj_inv=vbld_inv(j+nres)
1008             sig0ij=sigma(itypi,itypj)
1009             r0ij=r0(itypi,itypj)
1010             chi1=chi(itypi,itypj)
1011             chi2=chi(itypj,itypi)
1012             chi12=chi1*chi2
1013             chip1=chip(itypi)
1014             chip2=chip(itypj)
1015             chip12=chip1*chip2
1016             alf1=alp(itypi)
1017             alf2=alp(itypj)
1018             alf12=0.5D0*(alf1+alf2)
1019             xj=c(1,nres+j)-xi
1020             yj=c(2,nres+j)-yi
1021             zj=c(3,nres+j)-zi
1022             dxj=dc_norm(1,nres+j)
1023             dyj=dc_norm(2,nres+j)
1024             dzj=dc_norm(3,nres+j)
1025             rrij=1.0D0/(xj*xj+yj*yj+zj*zj)
1026             rij=dsqrt(rrij)
1027
1028             sss=sscale(1.0d0/(rij*sigmaii(itypi,itypj)))
1029
1030             if (sss.lt.1.0d0) then
1031
1032 C Calculate angle-dependent terms of energy and contributions to their
1033 C derivatives.
1034               call sc_angular
1035               sigsq=1.0D0/sigsq
1036               sig=sig0ij*dsqrt(sigsq)
1037               rij_shift=1.0D0/rij-sig+r0ij
1038 C I hate to put IF's in the loops, but here don't have another choice!!!!
1039               if (rij_shift.le.0.0D0) then
1040                 evdw=1.0D20
1041                 return
1042               endif
1043               sigder=-sig*sigsq
1044 c---------------------------------------------------------------
1045               rij_shift=1.0D0/rij_shift
1046               fac=rij_shift**expon
1047               e1=fac*fac*aa(itypi,itypj)
1048               e2=fac*bb(itypi,itypj)
1049               evdwij=eps1*eps2rt*eps3rt*(e1+e2)
1050               eps2der=evdwij*eps3rt
1051               eps3der=evdwij*eps2rt
1052               fac_augm=rrij**expon
1053               e_augm=augm(itypi,itypj)*fac_augm
1054               evdwij=evdwij*eps2rt*eps3rt
1055               evdw=evdw+(evdwij+e_augm)*(1.0d0-sss)
1056               if (lprn) then
1057               sigm=dabs(aa(itypi,itypj)/bb(itypi,itypj))**(1.0D0/6.0D0)
1058               epsi=bb(itypi,itypj)**2/aa(itypi,itypj)
1059               write (iout,'(2(a3,i3,2x),17(0pf7.3))')
1060      &          restyp(itypi),i,restyp(itypj),j,
1061      &          epsi,sigm,sig,(augm(itypi,itypj)/epsi)**(1.0D0/12.0D0),
1062      &          chi1,chi2,chip1,chip2,
1063      &          eps1,eps2rt**2,eps3rt**2,
1064      &          om1,om2,om12,1.0D0/rij,1.0D0/rij_shift,
1065      &          evdwij+e_augm
1066               endif
1067 C Calculate gradient components.
1068               e1=e1*eps1*eps2rt**2*eps3rt**2
1069               fac=-expon*(e1+evdwij)*rij_shift
1070               sigder=fac*sigder
1071               fac=rij*fac-2*expon*rrij*e_augm
1072 C Calculate the radial part of the gradient
1073               gg(1)=xj*fac
1074               gg(2)=yj*fac
1075               gg(3)=zj*fac
1076 C Calculate angular part of the gradient.
1077               call sc_grad_scale(1.0d0-sss)
1078             endif
1079           enddo      ! j
1080         enddo        ! iint
1081       enddo          ! i
1082       end
1083 C-----------------------------------------------------------------------------
1084       subroutine egbv_short(evdw)
1085 C
1086 C This subroutine calculates the interaction energy of nonbonded side chains
1087 C assuming the Gay-Berne-Vorobjev potential of interaction.
1088 C
1089       implicit real*8 (a-h,o-z)
1090       include 'DIMENSIONS'
1091       include 'COMMON.GEO'
1092       include 'COMMON.VAR'
1093       include 'COMMON.LOCAL'
1094       include 'COMMON.CHAIN'
1095       include 'COMMON.DERIV'
1096       include 'COMMON.NAMES'
1097       include 'COMMON.INTERACT'
1098       include 'COMMON.IOUNITS'
1099       include 'COMMON.CALC'
1100       common /srutu/ icall
1101       logical lprn
1102       evdw=0.0D0
1103 c     print *,'Entering EGB nnt=',nnt,' nct=',nct,' expon=',expon
1104       evdw=0.0D0
1105       lprn=.false.
1106 c     if (icall.eq.0) lprn=.true.
1107       ind=0
1108       do i=iatsc_s,iatsc_e
1109         itypi=itype(i)
1110         if (itypi.eq.ntyp1) cycle
1111         itypi1=itype(i+1)
1112         xi=c(1,nres+i)
1113         yi=c(2,nres+i)
1114         zi=c(3,nres+i)
1115         dxi=dc_norm(1,nres+i)
1116         dyi=dc_norm(2,nres+i)
1117         dzi=dc_norm(3,nres+i)
1118 c        dsci_inv=dsc_inv(itypi)
1119         dsci_inv=vbld_inv(i+nres)
1120 C
1121 C Calculate SC interaction energy.
1122 C
1123         do iint=1,nint_gr(i)
1124           do j=istart(i,iint),iend(i,iint)
1125             ind=ind+1
1126             itypj=itype(j)
1127             if (itypj.eq.ntyp1) cycle
1128 c            dscj_inv=dsc_inv(itypj)
1129             dscj_inv=vbld_inv(j+nres)
1130             sig0ij=sigma(itypi,itypj)
1131             r0ij=r0(itypi,itypj)
1132             chi1=chi(itypi,itypj)
1133             chi2=chi(itypj,itypi)
1134             chi12=chi1*chi2
1135             chip1=chip(itypi)
1136             chip2=chip(itypj)
1137             chip12=chip1*chip2
1138             alf1=alp(itypi)
1139             alf2=alp(itypj)
1140             alf12=0.5D0*(alf1+alf2)
1141             xj=c(1,nres+j)-xi
1142             yj=c(2,nres+j)-yi
1143             zj=c(3,nres+j)-zi
1144             dxj=dc_norm(1,nres+j)
1145             dyj=dc_norm(2,nres+j)
1146             dzj=dc_norm(3,nres+j)
1147             rrij=1.0D0/(xj*xj+yj*yj+zj*zj)
1148             rij=dsqrt(rrij)
1149
1150             sss=sscale(1.0d0/(rij*sigmaii(itypi,itypj)))
1151
1152             if (sss.gt.0.0d0) then
1153
1154 C Calculate angle-dependent terms of energy and contributions to their
1155 C derivatives.
1156               call sc_angular
1157               sigsq=1.0D0/sigsq
1158               sig=sig0ij*dsqrt(sigsq)
1159               rij_shift=1.0D0/rij-sig+r0ij
1160 C I hate to put IF's in the loops, but here don't have another choice!!!!
1161               if (rij_shift.le.0.0D0) then
1162                 evdw=1.0D20
1163                 return
1164               endif
1165               sigder=-sig*sigsq
1166 c---------------------------------------------------------------
1167               rij_shift=1.0D0/rij_shift
1168               fac=rij_shift**expon
1169               e1=fac*fac*aa(itypi,itypj)
1170               e2=fac*bb(itypi,itypj)
1171               evdwij=eps1*eps2rt*eps3rt*(e1+e2)
1172               eps2der=evdwij*eps3rt
1173               eps3der=evdwij*eps2rt
1174               fac_augm=rrij**expon
1175               e_augm=augm(itypi,itypj)*fac_augm
1176               evdwij=evdwij*eps2rt*eps3rt
1177               evdw=evdw+(evdwij+e_augm)*sss
1178               if (lprn) then
1179               sigm=dabs(aa(itypi,itypj)/bb(itypi,itypj))**(1.0D0/6.0D0)
1180               epsi=bb(itypi,itypj)**2/aa(itypi,itypj)
1181               write (iout,'(2(a3,i3,2x),17(0pf7.3))')
1182      &          restyp(itypi),i,restyp(itypj),j,
1183      &          epsi,sigm,sig,(augm(itypi,itypj)/epsi)**(1.0D0/12.0D0),
1184      &          chi1,chi2,chip1,chip2,
1185      &          eps1,eps2rt**2,eps3rt**2,
1186      &          om1,om2,om12,1.0D0/rij,1.0D0/rij_shift,
1187      &          evdwij+e_augm
1188               endif
1189 C Calculate gradient components.
1190               e1=e1*eps1*eps2rt**2*eps3rt**2
1191               fac=-expon*(e1+evdwij)*rij_shift
1192               sigder=fac*sigder
1193               fac=rij*fac-2*expon*rrij*e_augm
1194 C Calculate the radial part of the gradient
1195               gg(1)=xj*fac
1196               gg(2)=yj*fac
1197               gg(3)=zj*fac
1198 C Calculate angular part of the gradient.
1199               call sc_grad_scale(sss)
1200             endif
1201           enddo      ! j
1202         enddo        ! iint
1203       enddo          ! i
1204       end
1205 C----------------------------------------------------------------------------
1206       subroutine sc_grad_scale(scalfac)
1207       implicit real*8 (a-h,o-z)
1208       include 'DIMENSIONS'
1209       include 'COMMON.CHAIN'
1210       include 'COMMON.DERIV'
1211       include 'COMMON.CALC'
1212       include 'COMMON.IOUNITS'
1213       double precision dcosom1(3),dcosom2(3)
1214       double precision scalfac
1215       eom1=eps2der*eps2rt_om1-2.0D0*alf1*eps3der+sigder*sigsq_om1
1216       eom2=eps2der*eps2rt_om2+2.0D0*alf2*eps3der+sigder*sigsq_om2
1217       eom12=evdwij*eps1_om12+eps2der*eps2rt_om12
1218      &     -2.0D0*alf12*eps3der+sigder*sigsq_om12
1219 c diagnostics only
1220 c      eom1=0.0d0
1221 c      eom2=0.0d0
1222 c      eom12=evdwij*eps1_om12
1223 c end diagnostics
1224 c      write (iout,*) "eps2der",eps2der," eps3der",eps3der,
1225 c     &  " sigder",sigder
1226 c      write (iout,*) "eps1_om12",eps1_om12," eps2rt_om12",eps2rt_om12
1227 c      write (iout,*) "eom1",eom1," eom2",eom2," eom12",eom12
1228       do k=1,3
1229         dcosom1(k)=rij*(dc_norm(k,nres+i)-om1*erij(k))
1230         dcosom2(k)=rij*(dc_norm(k,nres+j)-om2*erij(k))
1231       enddo
1232       do k=1,3
1233         gg(k)=(gg(k)+eom1*dcosom1(k)+eom2*dcosom2(k))*scalfac
1234       enddo
1235 c      write (iout,*) "gg",(gg(k),k=1,3)
1236       do k=1,3
1237         gvdwx(k,i)=gvdwx(k,i)-gg(k)
1238      &            +(eom12*(dc_norm(k,nres+j)-om12*dc_norm(k,nres+i))
1239      &          +eom1*(erij(k)-om1*dc_norm(k,nres+i)))*dsci_inv*scalfac
1240         gvdwx(k,j)=gvdwx(k,j)+gg(k)
1241      &            +(eom12*(dc_norm(k,nres+i)-om12*dc_norm(k,nres+j))
1242      &          +eom2*(erij(k)-om2*dc_norm(k,nres+j)))*dscj_inv*scalfac
1243 c        write (iout,*)(eom12*(dc_norm(k,nres+j)-om12*dc_norm(k,nres+i))
1244 c     &            +eom1*(erij(k)-om1*dc_norm(k,nres+i)))*dsci_inv
1245 c        write (iout,*)(eom12*(dc_norm(k,nres+i)-om12*dc_norm(k,nres+j))
1246 c     &            +eom2*(erij(k)-om2*dc_norm(k,nres+j)))*dscj_inv
1247       enddo
1248 C
1249 C Calculate the components of the gradient in DC and X
1250 C
1251       do l=1,3
1252         gvdwc(l,i)=gvdwc(l,i)-gg(l)
1253         gvdwc(l,j)=gvdwc(l,j)+gg(l)
1254       enddo
1255       return
1256       end
1257 C--------------------------------------------------------------------------
1258       subroutine eelec_scale(ees,evdw1,eel_loc,eello_turn3,eello_turn4)
1259 C
1260 C This subroutine calculates the average interaction energy and its gradient
1261 C in the virtual-bond vectors between non-adjacent peptide groups, based on
1262 C the potential described in Liwo et al., Protein Sci., 1993, 2, 1715.
1263 C The potential depends both on the distance of peptide-group centers and on
1264 C the orientation of the CA-CA virtual bonds.
1265 C
1266       implicit real*8 (a-h,o-z)
1267 #ifdef MPI
1268       include 'mpif.h'
1269 #endif
1270       include 'DIMENSIONS'
1271       include 'COMMON.CONTROL'
1272       include 'COMMON.SETUP'
1273       include 'COMMON.IOUNITS'
1274       include 'COMMON.GEO'
1275       include 'COMMON.VAR'
1276       include 'COMMON.LOCAL'
1277       include 'COMMON.CHAIN'
1278       include 'COMMON.DERIV'
1279       include 'COMMON.INTERACT'
1280       include 'COMMON.CONTACTS'
1281       include 'COMMON.TORSION'
1282       include 'COMMON.VECTORS'
1283       include 'COMMON.FFIELD'
1284       include 'COMMON.TIME1'
1285       dimension ggg(3),gggp(3),gggm(3),erij(3),dcosb(3),dcosg(3),
1286      &          erder(3,3),uryg(3,3),urzg(3,3),vryg(3,3),vrzg(3,3)
1287       double precision acipa(2,2),agg(3,4),aggi(3,4),aggi1(3,4),
1288      &    aggj(3,4),aggj1(3,4),a_temp(2,2),muij(4)
1289       common /locel/ a_temp,agg,aggi,aggi1,aggj,aggj1,a22,a23,a32,a33,
1290      &    dxi,dyi,dzi,dx_normi,dy_normi,dz_normi,xmedi,ymedi,zmedi,
1291      &    num_conti,j1,j2
1292 c 4/26/02 - AL scaling factor for 1,4 repulsive VDW interactions
1293 #ifdef MOMENT
1294       double precision scal_el /1.0d0/
1295 #else
1296       double precision scal_el /0.5d0/
1297 #endif
1298 C 12/13/98
1299 C 13-go grudnia roku pamietnego...
1300       double precision unmat(3,3) /1.0d0,0.0d0,0.0d0,
1301      &                   0.0d0,1.0d0,0.0d0,
1302      &                   0.0d0,0.0d0,1.0d0/
1303 cd      write(iout,*) 'In EELEC'
1304 cd      do i=1,nloctyp
1305 cd        write(iout,*) 'Type',i
1306 cd        write(iout,*) 'B1',B1(:,i)
1307 cd        write(iout,*) 'B2',B2(:,i)
1308 cd        write(iout,*) 'CC',CC(:,:,i)
1309 cd        write(iout,*) 'DD',DD(:,:,i)
1310 cd        write(iout,*) 'EE',EE(:,:,i)
1311 cd      enddo
1312 cd      call check_vecgrad
1313 cd      stop
1314 C      print *,"WCHODZE3"
1315       if (icheckgrad.eq.1) then
1316         do i=1,nres-1
1317           fac=1.0d0/dsqrt(scalar(dc(1,i),dc(1,i)))
1318           do k=1,3
1319             dc_norm(k,i)=dc(k,i)*fac
1320           enddo
1321 c          write (iout,*) 'i',i,' fac',fac
1322         enddo
1323       endif
1324       if (wel_loc.gt.0.0d0 .or. wcorr4.gt.0.0d0 .or. wcorr5.gt.0.0d0
1325      &    .or. wcorr6.gt.0.0d0 .or. wturn3.gt.0.0d0 .or.
1326      &    wturn4.gt.0.0d0 .or. wturn6.gt.0.0d0) then
1327 c        call vec_and_deriv
1328 #ifdef TIMING
1329         time01=MPI_Wtime()
1330 #endif
1331         call set_matrices
1332 #ifdef TIMING
1333         time_mat=time_mat+MPI_Wtime()-time01
1334 #endif
1335       endif
1336 cd      do i=1,nres-1
1337 cd        write (iout,*) 'i=',i
1338 cd        do k=1,3
1339 cd        write (iout,'(i5,2f10.5)') k,uy(k,i),uz(k,i)
1340 cd        enddo
1341 cd        do k=1,3
1342 cd          write (iout,'(f10.5,2x,3f10.5,2x,3f10.5)')
1343 cd     &     uz(k,i),(uzgrad(k,l,1,i),l=1,3),(uzgrad(k,l,2,i),l=1,3)
1344 cd        enddo
1345 cd      enddo
1346       t_eelecij=0.0d0
1347       ees=0.0D0
1348       evdw1=0.0D0
1349       eel_loc=0.0d0
1350       eello_turn3=0.0d0
1351       eello_turn4=0.0d0
1352       ind=0
1353       do i=1,nres
1354         num_cont_hb(i)=0
1355       enddo
1356 cd      print '(a)','Enter EELEC'
1357 cd      write (iout,*) 'iatel_s=',iatel_s,' iatel_e=',iatel_e
1358       do i=1,nres
1359         gel_loc_loc(i)=0.0d0
1360         gcorr_loc(i)=0.0d0
1361       enddo
1362 c
1363 c
1364 c 9/27/08 AL Split the interaction loop to ensure load balancing of turn terms
1365 C
1366 C Loop over i,i+2 and i,i+3 pairs of the peptide groups
1367 C
1368       do i=iturn3_start,iturn3_end
1369         if (itype(i).eq.ntyp1.or. itype(i+1).eq.ntyp1
1370      &  .or. itype(i+2).eq.ntyp1 .or. itype(i+3).eq.ntyp1
1371      &  .or. itype(i-1).eq.ntyp1
1372      &  .or. itype(i+4).eq.ntyp1
1373      &   ) cycle
1374         dxi=dc(1,i)
1375         dyi=dc(2,i)
1376         dzi=dc(3,i)
1377         dx_normi=dc_norm(1,i)
1378         dy_normi=dc_norm(2,i)
1379         dz_normi=dc_norm(3,i)
1380         xmedi=c(1,i)+0.5d0*dxi
1381         ymedi=c(2,i)+0.5d0*dyi
1382         zmedi=c(3,i)+0.5d0*dzi
1383           xmedi=mod(xmedi,boxxsize)
1384           if (xmedi.lt.0) xmedi=xmedi+boxxsize
1385           ymedi=mod(ymedi,boxysize)
1386           if (ymedi.lt.0) ymedi=ymedi+boxysize
1387           zmedi=mod(zmedi,boxzsize)
1388           if (zmedi.lt.0) zmedi=zmedi+boxzsize
1389         num_conti=0
1390         call eelecij_scale(i,i+2,ees,evdw1,eel_loc)
1391         if (wturn3.gt.0.0d0) call eturn3(i,eello_turn3)
1392         num_cont_hb(i)=num_conti
1393       enddo
1394       do i=iturn4_start,iturn4_end
1395         if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1
1396      &    .or. itype(i+3).eq.ntyp1
1397      &    .or. itype(i+4).eq.ntyp1
1398      &    .or. itype(i+5).eq.ntyp1
1399      &    .or. itype(i-1).eq.ntyp1
1400      &    ) cycle
1401         dxi=dc(1,i)
1402         dyi=dc(2,i)
1403         dzi=dc(3,i)
1404         dx_normi=dc_norm(1,i)
1405         dy_normi=dc_norm(2,i)
1406         dz_normi=dc_norm(3,i)
1407         xmedi=c(1,i)+0.5d0*dxi
1408         ymedi=c(2,i)+0.5d0*dyi
1409         zmedi=c(3,i)+0.5d0*dzi
1410           xmedi=mod(xmedi,boxxsize)
1411           if (xmedi.lt.0) xmedi=xmedi+boxxsize
1412           ymedi=mod(ymedi,boxysize)
1413           if (ymedi.lt.0) ymedi=ymedi+boxysize
1414           zmedi=mod(zmedi,boxzsize)
1415           if (zmedi.lt.0) zmedi=zmedi+boxzsize
1416         num_conti=num_cont_hb(i)
1417         call eelecij_scale(i,i+3,ees,evdw1,eel_loc)
1418         if (wturn4.gt.0.0d0 .and. itype(i+2).ne.ntyp1)
1419      &    call eturn4(i,eello_turn4)
1420         num_cont_hb(i)=num_conti
1421       enddo   ! i
1422 c
1423 c Loop over all pairs of interacting peptide groups except i,i+2 and i,i+3
1424 c
1425       do i=iatel_s,iatel_e
1426         if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1
1427      &  .or. itype(i+2).eq.ntyp1
1428      &  .or. itype(i-1).eq.ntyp1
1429      &) cycle
1430         dxi=dc(1,i)
1431         dyi=dc(2,i)
1432         dzi=dc(3,i)
1433         dx_normi=dc_norm(1,i)
1434         dy_normi=dc_norm(2,i)
1435         dz_normi=dc_norm(3,i)
1436         xmedi=c(1,i)+0.5d0*dxi
1437         ymedi=c(2,i)+0.5d0*dyi
1438         zmedi=c(3,i)+0.5d0*dzi
1439           xmedi=mod(xmedi,boxxsize)
1440           if (xmedi.lt.0) xmedi=xmedi+boxxsize
1441           ymedi=mod(ymedi,boxysize)
1442           if (ymedi.lt.0) ymedi=ymedi+boxysize
1443           zmedi=mod(zmedi,boxzsize)
1444           if (zmedi.lt.0) zmedi=zmedi+boxzsize
1445 c        write (iout,*) 'i',i,' ielstart',ielstart(i),' ielend',ielend(i)
1446         num_conti=num_cont_hb(i)
1447         do j=ielstart(i),ielend(i)
1448           if (itype(j).eq.ntyp1 .or. itype(j+1).eq.ntyp1
1449      & .or.itype(j+2).eq.ntyp1
1450      & .or.itype(j-1).eq.ntyp1
1451      &) cycle
1452           call eelecij_scale(i,j,ees,evdw1,eel_loc)
1453         enddo ! j
1454         num_cont_hb(i)=num_conti
1455       enddo   ! i
1456 c      write (iout,*) "Number of loop steps in EELEC:",ind
1457 cd      do i=1,nres
1458 cd        write (iout,'(i3,3f10.5,5x,3f10.5)')
1459 cd     &     i,(gel_loc(k,i),k=1,3),gel_loc_loc(i)
1460 cd      enddo
1461 c 12/7/99 Adam eello_turn3 will be considered as a separate energy term
1462 ccc      eel_loc=eel_loc+eello_turn3
1463 cd      print *,"Processor",fg_rank," t_eelecij",t_eelecij
1464       return
1465       end
1466 C-------------------------------------------------------------------------------
1467       subroutine eelecij_scale(i,j,ees,evdw1,eel_loc)
1468       implicit real*8 (a-h,o-z)
1469       include 'DIMENSIONS'
1470 #ifdef MPI
1471       include "mpif.h"
1472 #endif
1473       include 'COMMON.CONTROL'
1474       include 'COMMON.IOUNITS'
1475       include 'COMMON.GEO'
1476       include 'COMMON.VAR'
1477       include 'COMMON.LOCAL'
1478       include 'COMMON.CHAIN'
1479       include 'COMMON.DERIV'
1480       include 'COMMON.INTERACT'
1481       include 'COMMON.CONTACTS'
1482       include 'COMMON.TORSION'
1483       include 'COMMON.VECTORS'
1484       include 'COMMON.FFIELD'
1485       include 'COMMON.TIME1'
1486       dimension ggg(3),gggp(3),gggm(3),erij(3),dcosb(3),dcosg(3),
1487      &          erder(3,3),uryg(3,3),urzg(3,3),vryg(3,3),vrzg(3,3)
1488       double precision acipa(2,2),agg(3,4),aggi(3,4),aggi1(3,4),
1489      &    aggj(3,4),aggj1(3,4),a_temp(2,2),muij(4)
1490       common /locel/ a_temp,agg,aggi,aggi1,aggj,aggj1,a22,a23,a32,a33,
1491      &    dxi,dyi,dzi,dx_normi,dy_normi,dz_normi,xmedi,ymedi,zmedi,
1492      &    num_conti,j1,j2
1493 c 4/26/02 - AL scaling factor for 1,4 repulsive VDW interactions
1494 #ifdef MOMENT
1495       double precision scal_el /1.0d0/
1496 #else
1497       double precision scal_el /0.5d0/
1498 #endif
1499 C 12/13/98
1500 C 13-go grudnia roku pamietnego...
1501       double precision unmat(3,3) /1.0d0,0.0d0,0.0d0,
1502      &                   0.0d0,1.0d0,0.0d0,
1503      &                   0.0d0,0.0d0,1.0d0/
1504 c          time00=MPI_Wtime()
1505 cd      write (iout,*) "eelecij",i,j
1506 C      print *,"WCHODZE2"
1507           ind=ind+1
1508           iteli=itel(i)
1509           itelj=itel(j)
1510           if (j.eq.i+2 .and. itelj.eq.2) iteli=2
1511           aaa=app(iteli,itelj)
1512           bbb=bpp(iteli,itelj)
1513           ael6i=ael6(iteli,itelj)
1514           ael3i=ael3(iteli,itelj)
1515           dxj=dc(1,j)
1516           dyj=dc(2,j)
1517           dzj=dc(3,j)
1518           dx_normj=dc_norm(1,j)
1519           dy_normj=dc_norm(2,j)
1520           dz_normj=dc_norm(3,j)
1521           xj=c(1,j)+0.5D0*dxj
1522           yj=c(2,j)+0.5D0*dyj
1523           zj=c(3,j)+0.5D0*dzj
1524          xj=mod(xj,boxxsize)
1525           if (xj.lt.0) xj=xj+boxxsize
1526           yj=mod(yj,boxysize)
1527           if (yj.lt.0) yj=yj+boxysize
1528           zj=mod(zj,boxzsize)
1529           if (zj.lt.0) zj=zj+boxzsize
1530       dist_init=(xj-xmedi)**2+(yj-ymedi)**2+(zj-zmedi)**2
1531       xj_safe=xj
1532       yj_safe=yj
1533       zj_safe=zj
1534       isubchap=0
1535       do xshift=-1,1
1536       do yshift=-1,1
1537       do zshift=-1,1
1538           xj=xj_safe+xshift*boxxsize
1539           yj=yj_safe+yshift*boxysize
1540           zj=zj_safe+zshift*boxzsize
1541           dist_temp=(xj-xmedi)**2+(yj-ymedi)**2+(zj-zmedi)**2
1542           if(dist_temp.lt.dist_init) then
1543             dist_init=dist_temp
1544             xj_temp=xj
1545             yj_temp=yj
1546             zj_temp=zj
1547             isubchap=1
1548           endif
1549        enddo
1550        enddo
1551        enddo
1552        if (isubchap.eq.1) then
1553           xj=xj_temp-xmedi
1554           yj=yj_temp-ymedi
1555           zj=zj_temp-zmedi
1556        else
1557           xj=xj_safe-xmedi
1558           yj=yj_safe-ymedi
1559           zj=zj_safe-zmedi
1560        endif
1561
1562           rij=xj*xj+yj*yj+zj*zj
1563           rrmij=1.0D0/rij
1564           rij=dsqrt(rij)
1565           rmij=1.0D0/rij
1566 c For extracting the short-range part of Evdwpp
1567           sss=sscale(rij/rpp(iteli,itelj))
1568           sssgrad=sscagrad(rij/rpp(iteli,itelj))
1569           r3ij=rrmij*rmij
1570           r6ij=r3ij*r3ij
1571           cosa=dx_normi*dx_normj+dy_normi*dy_normj+dz_normi*dz_normj
1572           cosb=(xj*dx_normi+yj*dy_normi+zj*dz_normi)*rmij
1573           cosg=(xj*dx_normj+yj*dy_normj+zj*dz_normj)*rmij
1574           fac=cosa-3.0D0*cosb*cosg
1575           ev1=aaa*r6ij*r6ij
1576 c 4/26/02 - AL scaling down 1,4 repulsive VDW interactions
1577           if (j.eq.i+2) ev1=scal_el*ev1
1578           ev2=bbb*r6ij
1579           fac3=ael6i*r6ij
1580           fac4=ael3i*r3ij
1581           evdwij=ev1+ev2
1582           el1=fac3*(4.0D0+fac*fac-3.0D0*(cosb*cosb+cosg*cosg))
1583           el2=fac4*fac
1584           eesij=el1+el2
1585 C 12/26/95 - for the evaluation of multi-body H-bonding interactions
1586           ees0ij=4.0D0+fac*fac-3.0D0*(cosb*cosb+cosg*cosg)
1587           ees=ees+eesij
1588           evdw1=evdw1+evdwij*(1.0d0-sss)
1589 cd          write(iout,'(2(2i3,2x),7(1pd12.4)/2(3(1pd12.4),5x)/)')
1590 cd     &      iteli,i,itelj,j,aaa,bbb,ael6i,ael3i,
1591 cd     &      1.0D0/dsqrt(rrmij),evdwij,eesij,
1592 cd     &      xmedi,ymedi,zmedi,xj,yj,zj
1593
1594           if (energy_dec) then
1595               write (iout,'(a6,2i5,0pf7.3,f7.3)') 'evdw1',i,j,evdwij,sss
1596               write (iout,'(a6,2i5,0pf7.3)') 'ees',i,j,eesij
1597           endif
1598
1599 C
1600 C Calculate contributions to the Cartesian gradient.
1601 C
1602 #ifdef SPLITELE
1603           facvdw=-6*rrmij*(ev1+evdwij)*(1.0d0-sss)
1604           facel=-3*rrmij*(el1+eesij)
1605           fac1=fac
1606           erij(1)=xj*rmij
1607           erij(2)=yj*rmij
1608           erij(3)=zj*rmij
1609 *
1610 * Radial derivatives. First process both termini of the fragment (i,j)
1611 *
1612           ggg(1)=facel*xj
1613           ggg(2)=facel*yj
1614           ggg(3)=facel*zj
1615 c          do k=1,3
1616 c            ghalf=0.5D0*ggg(k)
1617 c            gelc(k,i)=gelc(k,i)+ghalf
1618 c            gelc(k,j)=gelc(k,j)+ghalf
1619 c          enddo
1620 c 9/28/08 AL Gradient compotents will be summed only at the end
1621           do k=1,3
1622             gelc_long(k,j)=gelc_long(k,j)+ggg(k)
1623             gelc_long(k,i)=gelc_long(k,i)-ggg(k)
1624           enddo
1625 *
1626 * Loop over residues i+1 thru j-1.
1627 *
1628 cgrad          do k=i+1,j-1
1629 cgrad            do l=1,3
1630 cgrad              gelc(l,k)=gelc(l,k)+ggg(l)
1631 cgrad            enddo
1632 cgrad          enddo
1633           ggg(1)=facvdw*xj-sssgrad*rmij*evdwij*xj/rpp(iteli,itelj)
1634           ggg(2)=facvdw*yj-sssgrad*rmij*evdwij*yj/rpp(iteli,itelj)
1635           ggg(3)=facvdw*zj-sssgrad*rmij*evdwij*zj/rpp(iteli,itelj)
1636 c          do k=1,3
1637 c            ghalf=0.5D0*ggg(k)
1638 c            gvdwpp(k,i)=gvdwpp(k,i)+ghalf
1639 c            gvdwpp(k,j)=gvdwpp(k,j)+ghalf
1640 c          enddo
1641 c 9/28/08 AL Gradient compotents will be summed only at the end
1642           do k=1,3
1643             gvdwpp(k,j)=gvdwpp(k,j)+ggg(k)
1644             gvdwpp(k,i)=gvdwpp(k,i)-ggg(k)
1645           enddo
1646 *
1647 * Loop over residues i+1 thru j-1.
1648 *
1649 cgrad          do k=i+1,j-1
1650 cgrad            do l=1,3
1651 cgrad              gvdwpp(l,k)=gvdwpp(l,k)+ggg(l)
1652 cgrad            enddo
1653 cgrad          enddo
1654 #else
1655           facvdw=ev1+evdwij*(1.0d0-sss)
1656           facel=el1+eesij
1657           fac1=fac
1658           fac=-3*rrmij*(facvdw+facvdw+facel)
1659           erij(1)=xj*rmij
1660           erij(2)=yj*rmij
1661           erij(3)=zj*rmij
1662 *
1663 * Radial derivatives. First process both termini of the fragment (i,j)
1664 *
1665           ggg(1)=fac*xj
1666           ggg(2)=fac*yj
1667           ggg(3)=fac*zj
1668 c          do k=1,3
1669 c            ghalf=0.5D0*ggg(k)
1670 c            gelc(k,i)=gelc(k,i)+ghalf
1671 c            gelc(k,j)=gelc(k,j)+ghalf
1672 c          enddo
1673 c 9/28/08 AL Gradient compotents will be summed only at the end
1674           do k=1,3
1675             gelc_long(k,j)=gelc(k,j)+ggg(k)
1676             gelc_long(k,i)=gelc(k,i)-ggg(k)
1677           enddo
1678 *
1679 * Loop over residues i+1 thru j-1.
1680 *
1681 cgrad          do k=i+1,j-1
1682 cgrad            do l=1,3
1683 cgrad              gelc(l,k)=gelc(l,k)+ggg(l)
1684 cgrad            enddo
1685 cgrad          enddo
1686 c 9/28/08 AL Gradient compotents will be summed only at the end
1687 C          ggg(1)=facvdw*xj
1688 C          ggg(2)=facvdw*yj
1689 C          ggg(3)=facvdw*zj
1690           ggg(1)=facvdw*xj-sssgrad*rmij*evdwij*xj/rpp(iteli,itelj)
1691           ggg(2)=facvdw*yj-sssgrad*rmij*evdwij*yj/rpp(iteli,itelj)
1692           ggg(3)=facvdw*zj-sssgrad*rmij*evdwij*zj/rpp(iteli,itelj)
1693           do k=1,3
1694             gvdwpp(k,j)=gvdwpp(k,j)+ggg(k)
1695             gvdwpp(k,i)=gvdwpp(k,i)-ggg(k)
1696           enddo
1697 #endif
1698 *
1699 * Angular part
1700 *
1701           ecosa=2.0D0*fac3*fac1+fac4
1702           fac4=-3.0D0*fac4
1703           fac3=-6.0D0*fac3
1704           ecosb=(fac3*(fac1*cosg+cosb)+cosg*fac4)
1705           ecosg=(fac3*(fac1*cosb+cosg)+cosb*fac4)
1706           do k=1,3
1707             dcosb(k)=rmij*(dc_norm(k,i)-erij(k)*cosb)
1708             dcosg(k)=rmij*(dc_norm(k,j)-erij(k)*cosg)
1709           enddo
1710 cd        print '(2i3,2(3(1pd14.5),3x))',i,j,(dcosb(k),k=1,3),
1711 cd   &          (dcosg(k),k=1,3)
1712           do k=1,3
1713             ggg(k)=ecosb*dcosb(k)+ecosg*dcosg(k)
1714           enddo
1715 c          do k=1,3
1716 c            ghalf=0.5D0*ggg(k)
1717 c            gelc(k,i)=gelc(k,i)+ghalf
1718 c     &               +(ecosa*(dc_norm(k,j)-cosa*dc_norm(k,i))
1719 c     &               + ecosb*(erij(k)-cosb*dc_norm(k,i)))*vbld_inv(i+1)
1720 c            gelc(k,j)=gelc(k,j)+ghalf
1721 c     &               +(ecosa*(dc_norm(k,i)-cosa*dc_norm(k,j))
1722 c     &               + ecosg*(erij(k)-cosg*dc_norm(k,j)))*vbld_inv(j+1)
1723 c          enddo
1724 cgrad          do k=i+1,j-1
1725 cgrad            do l=1,3
1726 cgrad              gelc(l,k)=gelc(l,k)+ggg(l)
1727 cgrad            enddo
1728 cgrad          enddo
1729           do k=1,3
1730             gelc(k,i)=gelc(k,i)
1731      &               +(ecosa*(dc_norm(k,j)-cosa*dc_norm(k,i))
1732      &               + ecosb*(erij(k)-cosb*dc_norm(k,i)))*vbld_inv(i+1)
1733             gelc(k,j)=gelc(k,j)
1734      &               +(ecosa*(dc_norm(k,i)-cosa*dc_norm(k,j))
1735      &               + ecosg*(erij(k)-cosg*dc_norm(k,j)))*vbld_inv(j+1)
1736             gelc_long(k,j)=gelc_long(k,j)+ggg(k)
1737             gelc_long(k,i)=gelc_long(k,i)-ggg(k)
1738           enddo
1739           IF (wel_loc.gt.0.0d0 .or. wcorr4.gt.0.0d0 .or. wcorr5.gt.0.0d0
1740      &        .or. wcorr6.gt.0.0d0 .or. wturn3.gt.0.0d0
1741      &        .or. wturn4.gt.0.0d0 .or. wturn6.gt.0.0d0) THEN
1742 C
1743 C 9/25/99 Mixed third-order local-electrostatic terms. The local-interaction
1744 C   energy of a peptide unit is assumed in the form of a second-order
1745 C   Fourier series in the angles lambda1 and lambda2 (see Nishikawa et al.
1746 C   Macromolecules, 1974, 7, 797-806 for definition). This correlation terms
1747 C   are computed for EVERY pair of non-contiguous peptide groups.
1748 C
1749           if (j.lt.nres-1) then
1750             j1=j+1
1751             j2=j-1
1752           else
1753             j1=j-1
1754             j2=j-2
1755           endif
1756           kkk=0
1757           do k=1,2
1758             do l=1,2
1759               kkk=kkk+1
1760               muij(kkk)=mu(k,i)*mu(l,j)
1761             enddo
1762           enddo
1763 cd         write (iout,*) 'EELEC: i',i,' j',j
1764 cd          write (iout,*) 'j',j,' j1',j1,' j2',j2
1765 cd          write(iout,*) 'muij',muij
1766           ury=scalar(uy(1,i),erij)
1767           urz=scalar(uz(1,i),erij)
1768           vry=scalar(uy(1,j),erij)
1769           vrz=scalar(uz(1,j),erij)
1770           a22=scalar(uy(1,i),uy(1,j))-3*ury*vry
1771           a23=scalar(uy(1,i),uz(1,j))-3*ury*vrz
1772           a32=scalar(uz(1,i),uy(1,j))-3*urz*vry
1773           a33=scalar(uz(1,i),uz(1,j))-3*urz*vrz
1774           fac=dsqrt(-ael6i)*r3ij
1775           a22=a22*fac
1776           a23=a23*fac
1777           a32=a32*fac
1778           a33=a33*fac
1779 cd          write (iout,'(4i5,4f10.5)')
1780 cd     &     i,itortyp(itype(i)),j,itortyp(itype(j)),a22,a23,a32,a33
1781 cd          write (iout,'(6f10.5)') (muij(k),k=1,4),fac,eel_loc_ij
1782 cd          write (iout,'(2(3f10.5,5x)/2(3f10.5,5x))') uy(:,i),uz(:,i),
1783 cd     &      uy(:,j),uz(:,j)
1784 cd          write (iout,'(4f10.5)')
1785 cd     &      scalar(uy(1,i),uy(1,j)),scalar(uy(1,i),uz(1,j)),
1786 cd     &      scalar(uz(1,i),uy(1,j)),scalar(uz(1,i),uz(1,j))
1787 cd          write (iout,'(4f10.5)') ury,urz,vry,vrz
1788 cd           write (iout,'(9f10.5/)')
1789 cd     &      fac22,a22,fac23,a23,fac32,a32,fac33,a33,eel_loc_ij
1790 C Derivatives of the elements of A in virtual-bond vectors
1791           call unormderiv(erij(1),unmat(1,1),rmij,erder(1,1))
1792           do k=1,3
1793             uryg(k,1)=scalar(erder(1,k),uy(1,i))
1794             uryg(k,2)=scalar(uygrad(1,k,1,i),erij(1))
1795             uryg(k,3)=scalar(uygrad(1,k,2,i),erij(1))
1796             urzg(k,1)=scalar(erder(1,k),uz(1,i))
1797             urzg(k,2)=scalar(uzgrad(1,k,1,i),erij(1))
1798             urzg(k,3)=scalar(uzgrad(1,k,2,i),erij(1))
1799             vryg(k,1)=scalar(erder(1,k),uy(1,j))
1800             vryg(k,2)=scalar(uygrad(1,k,1,j),erij(1))
1801             vryg(k,3)=scalar(uygrad(1,k,2,j),erij(1))
1802             vrzg(k,1)=scalar(erder(1,k),uz(1,j))
1803             vrzg(k,2)=scalar(uzgrad(1,k,1,j),erij(1))
1804             vrzg(k,3)=scalar(uzgrad(1,k,2,j),erij(1))
1805           enddo
1806 C Compute radial contributions to the gradient
1807           facr=-3.0d0*rrmij
1808           a22der=a22*facr
1809           a23der=a23*facr
1810           a32der=a32*facr
1811           a33der=a33*facr
1812           agg(1,1)=a22der*xj
1813           agg(2,1)=a22der*yj
1814           agg(3,1)=a22der*zj
1815           agg(1,2)=a23der*xj
1816           agg(2,2)=a23der*yj
1817           agg(3,2)=a23der*zj
1818           agg(1,3)=a32der*xj
1819           agg(2,3)=a32der*yj
1820           agg(3,3)=a32der*zj
1821           agg(1,4)=a33der*xj
1822           agg(2,4)=a33der*yj
1823           agg(3,4)=a33der*zj
1824 C Add the contributions coming from er
1825           fac3=-3.0d0*fac
1826           do k=1,3
1827             agg(k,1)=agg(k,1)+fac3*(uryg(k,1)*vry+vryg(k,1)*ury)
1828             agg(k,2)=agg(k,2)+fac3*(uryg(k,1)*vrz+vrzg(k,1)*ury)
1829             agg(k,3)=agg(k,3)+fac3*(urzg(k,1)*vry+vryg(k,1)*urz)
1830             agg(k,4)=agg(k,4)+fac3*(urzg(k,1)*vrz+vrzg(k,1)*urz)
1831           enddo
1832           do k=1,3
1833 C Derivatives in DC(i)
1834 cgrad            ghalf1=0.5d0*agg(k,1)
1835 cgrad            ghalf2=0.5d0*agg(k,2)
1836 cgrad            ghalf3=0.5d0*agg(k,3)
1837 cgrad            ghalf4=0.5d0*agg(k,4)
1838             aggi(k,1)=fac*(scalar(uygrad(1,k,1,i),uy(1,j))
1839      &      -3.0d0*uryg(k,2)*vry)!+ghalf1
1840             aggi(k,2)=fac*(scalar(uygrad(1,k,1,i),uz(1,j))
1841      &      -3.0d0*uryg(k,2)*vrz)!+ghalf2
1842             aggi(k,3)=fac*(scalar(uzgrad(1,k,1,i),uy(1,j))
1843      &      -3.0d0*urzg(k,2)*vry)!+ghalf3
1844             aggi(k,4)=fac*(scalar(uzgrad(1,k,1,i),uz(1,j))
1845      &      -3.0d0*urzg(k,2)*vrz)!+ghalf4
1846 C Derivatives in DC(i+1)
1847             aggi1(k,1)=fac*(scalar(uygrad(1,k,2,i),uy(1,j))
1848      &      -3.0d0*uryg(k,3)*vry)!+agg(k,1)
1849             aggi1(k,2)=fac*(scalar(uygrad(1,k,2,i),uz(1,j))
1850      &      -3.0d0*uryg(k,3)*vrz)!+agg(k,2)
1851             aggi1(k,3)=fac*(scalar(uzgrad(1,k,2,i),uy(1,j))
1852      &      -3.0d0*urzg(k,3)*vry)!+agg(k,3)
1853             aggi1(k,4)=fac*(scalar(uzgrad(1,k,2,i),uz(1,j))
1854      &      -3.0d0*urzg(k,3)*vrz)!+agg(k,4)
1855 C Derivatives in DC(j)
1856             aggj(k,1)=fac*(scalar(uygrad(1,k,1,j),uy(1,i))
1857      &      -3.0d0*vryg(k,2)*ury)!+ghalf1
1858             aggj(k,2)=fac*(scalar(uzgrad(1,k,1,j),uy(1,i))
1859      &      -3.0d0*vrzg(k,2)*ury)!+ghalf2
1860             aggj(k,3)=fac*(scalar(uygrad(1,k,1,j),uz(1,i))
1861      &      -3.0d0*vryg(k,2)*urz)!+ghalf3
1862             aggj(k,4)=fac*(scalar(uzgrad(1,k,1,j),uz(1,i))
1863      &      -3.0d0*vrzg(k,2)*urz)!+ghalf4
1864 C Derivatives in DC(j+1) or DC(nres-1)
1865             aggj1(k,1)=fac*(scalar(uygrad(1,k,2,j),uy(1,i))
1866      &      -3.0d0*vryg(k,3)*ury)
1867             aggj1(k,2)=fac*(scalar(uzgrad(1,k,2,j),uy(1,i))
1868      &      -3.0d0*vrzg(k,3)*ury)
1869             aggj1(k,3)=fac*(scalar(uygrad(1,k,2,j),uz(1,i))
1870      &      -3.0d0*vryg(k,3)*urz)
1871             aggj1(k,4)=fac*(scalar(uzgrad(1,k,2,j),uz(1,i))
1872      &      -3.0d0*vrzg(k,3)*urz)
1873 cgrad            if (j.eq.nres-1 .and. i.lt.j-2) then
1874 cgrad              do l=1,4
1875 cgrad                aggj1(k,l)=aggj1(k,l)+agg(k,l)
1876 cgrad              enddo
1877 cgrad            endif
1878           enddo
1879           acipa(1,1)=a22
1880           acipa(1,2)=a23
1881           acipa(2,1)=a32
1882           acipa(2,2)=a33
1883           a22=-a22
1884           a23=-a23
1885           do l=1,2
1886             do k=1,3
1887               agg(k,l)=-agg(k,l)
1888               aggi(k,l)=-aggi(k,l)
1889               aggi1(k,l)=-aggi1(k,l)
1890               aggj(k,l)=-aggj(k,l)
1891               aggj1(k,l)=-aggj1(k,l)
1892             enddo
1893           enddo
1894           if (j.lt.nres-1) then
1895             a22=-a22
1896             a32=-a32
1897             do l=1,3,2
1898               do k=1,3
1899                 agg(k,l)=-agg(k,l)
1900                 aggi(k,l)=-aggi(k,l)
1901                 aggi1(k,l)=-aggi1(k,l)
1902                 aggj(k,l)=-aggj(k,l)
1903                 aggj1(k,l)=-aggj1(k,l)
1904               enddo
1905             enddo
1906           else
1907             a22=-a22
1908             a23=-a23
1909             a32=-a32
1910             a33=-a33
1911             do l=1,4
1912               do k=1,3
1913                 agg(k,l)=-agg(k,l)
1914                 aggi(k,l)=-aggi(k,l)
1915                 aggi1(k,l)=-aggi1(k,l)
1916                 aggj(k,l)=-aggj(k,l)
1917                 aggj1(k,l)=-aggj1(k,l)
1918               enddo
1919             enddo
1920           endif
1921           ENDIF ! WCORR
1922           IF (wel_loc.gt.0.0d0) THEN
1923 C Contribution to the local-electrostatic energy coming from the i-j pair
1924           eel_loc_ij=a22*muij(1)+a23*muij(2)+a32*muij(3)
1925      &     +a33*muij(4)
1926 cd          write (iout,*) 'i',i,' j',j,' eel_loc_ij',eel_loc_ij
1927
1928           if (energy_dec) write (iout,'(a6,2i5,0pf7.3)')
1929      &            'eelloc',i,j,eel_loc_ij
1930
1931           eel_loc=eel_loc+eel_loc_ij
1932 C Partial derivatives in virtual-bond dihedral angles gamma
1933           if (i.gt.1)
1934      &    gel_loc_loc(i-1)=gel_loc_loc(i-1)+
1935      &            a22*muder(1,i)*mu(1,j)+a23*muder(1,i)*mu(2,j)
1936      &           +a32*muder(2,i)*mu(1,j)+a33*muder(2,i)*mu(2,j)
1937           gel_loc_loc(j-1)=gel_loc_loc(j-1)+
1938      &            a22*mu(1,i)*muder(1,j)+a23*mu(1,i)*muder(2,j)
1939      &           +a32*mu(2,i)*muder(1,j)+a33*mu(2,i)*muder(2,j)
1940 C Derivatives of eello in DC(i+1) thru DC(j-1) or DC(nres-2)
1941           do l=1,3
1942             ggg(l)=agg(l,1)*muij(1)+
1943      &          agg(l,2)*muij(2)+agg(l,3)*muij(3)+agg(l,4)*muij(4)
1944             gel_loc_long(l,j)=gel_loc_long(l,j)+ggg(l)
1945             gel_loc_long(l,i)=gel_loc_long(l,i)-ggg(l)
1946 cgrad            ghalf=0.5d0*ggg(l)
1947 cgrad            gel_loc(l,i)=gel_loc(l,i)+ghalf
1948 cgrad            gel_loc(l,j)=gel_loc(l,j)+ghalf
1949           enddo
1950 cgrad          do k=i+1,j2
1951 cgrad            do l=1,3
1952 cgrad              gel_loc(l,k)=gel_loc(l,k)+ggg(l)
1953 cgrad            enddo
1954 cgrad          enddo
1955 C Remaining derivatives of eello
1956           do l=1,3
1957             gel_loc(l,i)=gel_loc(l,i)+aggi(l,1)*muij(1)+
1958      &          aggi(l,2)*muij(2)+aggi(l,3)*muij(3)+aggi(l,4)*muij(4)
1959             gel_loc(l,i+1)=gel_loc(l,i+1)+aggi1(l,1)*muij(1)+
1960      &          aggi1(l,2)*muij(2)+aggi1(l,3)*muij(3)+aggi1(l,4)*muij(4)
1961             gel_loc(l,j)=gel_loc(l,j)+aggj(l,1)*muij(1)+
1962      &          aggj(l,2)*muij(2)+aggj(l,3)*muij(3)+aggj(l,4)*muij(4)
1963             gel_loc(l,j1)=gel_loc(l,j1)+aggj1(l,1)*muij(1)+
1964      &          aggj1(l,2)*muij(2)+aggj1(l,3)*muij(3)+aggj1(l,4)*muij(4)
1965           enddo
1966           ENDIF
1967 C Change 12/26/95 to calculate four-body contributions to H-bonding energy
1968 c          if (j.gt.i+1 .and. num_conti.le.maxconts) then
1969           if (wcorr+wcorr4+wcorr5+wcorr6.gt.0.0d0
1970      &       .and. num_conti.le.maxconts) then
1971 c            write (iout,*) i,j," entered corr"
1972 C
1973 C Calculate the contact function. The ith column of the array JCONT will
1974 C contain the numbers of atoms that make contacts with the atom I (of numbers
1975 C greater than I). The arrays FACONT and GACONT will contain the values of
1976 C the contact function and its derivative.
1977 c           r0ij=1.02D0*rpp(iteli,itelj)
1978 c           r0ij=1.11D0*rpp(iteli,itelj)
1979             r0ij=2.20D0*rpp(iteli,itelj)
1980 c           r0ij=1.55D0*rpp(iteli,itelj)
1981             call gcont(rij,r0ij,1.0D0,0.2d0*r0ij,fcont,fprimcont)
1982             if (fcont.gt.0.0D0) then
1983               num_conti=num_conti+1
1984               if (num_conti.gt.maxconts) then
1985                 write (iout,*) 'WARNING - max. # of contacts exceeded;',
1986      &                         ' will skip next contacts for this conf.'
1987               else
1988                 jcont_hb(num_conti,i)=j
1989 cd                write (iout,*) "i",i," j",j," num_conti",num_conti,
1990 cd     &           " jcont_hb",jcont_hb(num_conti,i)
1991                 IF (wcorr4.gt.0.0d0 .or. wcorr5.gt.0.0d0 .or.
1992      &          wcorr6.gt.0.0d0 .or. wturn6.gt.0.0d0) THEN
1993 C 9/30/99 (AL) - store components necessary to evaluate higher-order loc-el
1994 C  terms.
1995                 d_cont(num_conti,i)=rij
1996 cd                write (2,'(3e15.5)') rij,r0ij+0.2d0*r0ij,rij
1997 C     --- Electrostatic-interaction matrix ---
1998                 a_chuj(1,1,num_conti,i)=a22
1999                 a_chuj(1,2,num_conti,i)=a23
2000                 a_chuj(2,1,num_conti,i)=a32
2001                 a_chuj(2,2,num_conti,i)=a33
2002 C     --- Gradient of rij
2003                 do kkk=1,3
2004                   grij_hb_cont(kkk,num_conti,i)=erij(kkk)
2005                 enddo
2006                 kkll=0
2007                 do k=1,2
2008                   do l=1,2
2009                     kkll=kkll+1
2010                     do m=1,3
2011                       a_chuj_der(k,l,m,1,num_conti,i)=agg(m,kkll)
2012                       a_chuj_der(k,l,m,2,num_conti,i)=aggi(m,kkll)
2013                       a_chuj_der(k,l,m,3,num_conti,i)=aggi1(m,kkll)
2014                       a_chuj_der(k,l,m,4,num_conti,i)=aggj(m,kkll)
2015                       a_chuj_der(k,l,m,5,num_conti,i)=aggj1(m,kkll)
2016                     enddo
2017                   enddo
2018                 enddo
2019                 ENDIF
2020                 IF (wcorr4.eq.0.0d0 .and. wcorr.gt.0.0d0) THEN
2021 C Calculate contact energies
2022                 cosa4=4.0D0*cosa
2023                 wij=cosa-3.0D0*cosb*cosg
2024                 cosbg1=cosb+cosg
2025                 cosbg2=cosb-cosg
2026 c               fac3=dsqrt(-ael6i)/r0ij**3
2027                 fac3=dsqrt(-ael6i)*r3ij
2028 c                 ees0pij=dsqrt(4.0D0+cosa4+wij*wij-3.0D0*cosbg1*cosbg1)
2029                 ees0tmp=4.0D0+cosa4+wij*wij-3.0D0*cosbg1*cosbg1
2030                 if (ees0tmp.gt.0) then
2031                   ees0pij=dsqrt(ees0tmp)
2032                 else
2033                   ees0pij=0
2034                 endif
2035 c                ees0mij=dsqrt(4.0D0-cosa4+wij*wij-3.0D0*cosbg2*cosbg2)
2036                 ees0tmp=4.0D0-cosa4+wij*wij-3.0D0*cosbg2*cosbg2
2037                 if (ees0tmp.gt.0) then
2038                   ees0mij=dsqrt(ees0tmp)
2039                 else
2040                   ees0mij=0
2041                 endif
2042 c               ees0mij=0.0D0
2043                 ees0p(num_conti,i)=0.5D0*fac3*(ees0pij+ees0mij)
2044                 ees0m(num_conti,i)=0.5D0*fac3*(ees0pij-ees0mij)
2045 C Diagnostics. Comment out or remove after debugging!
2046 c               ees0p(num_conti,i)=0.5D0*fac3*ees0pij
2047 c               ees0m(num_conti,i)=0.5D0*fac3*ees0mij
2048 c               ees0m(num_conti,i)=0.0D0
2049 C End diagnostics.
2050 c               write (iout,*) 'i=',i,' j=',j,' rij=',rij,' r0ij=',r0ij,
2051 c    & ' ees0ij=',ees0p(num_conti,i),ees0m(num_conti,i),' fcont=',fcont
2052 C Angular derivatives of the contact function
2053                 ees0pij1=fac3/ees0pij
2054                 ees0mij1=fac3/ees0mij
2055                 fac3p=-3.0D0*fac3*rrmij
2056                 ees0pijp=0.5D0*fac3p*(ees0pij+ees0mij)
2057                 ees0mijp=0.5D0*fac3p*(ees0pij-ees0mij)
2058 c               ees0mij1=0.0D0
2059                 ecosa1=       ees0pij1*( 1.0D0+0.5D0*wij)
2060                 ecosb1=-1.5D0*ees0pij1*(wij*cosg+cosbg1)
2061                 ecosg1=-1.5D0*ees0pij1*(wij*cosb+cosbg1)
2062                 ecosa2=       ees0mij1*(-1.0D0+0.5D0*wij)
2063                 ecosb2=-1.5D0*ees0mij1*(wij*cosg+cosbg2)
2064                 ecosg2=-1.5D0*ees0mij1*(wij*cosb-cosbg2)
2065                 ecosap=ecosa1+ecosa2
2066                 ecosbp=ecosb1+ecosb2
2067                 ecosgp=ecosg1+ecosg2
2068                 ecosam=ecosa1-ecosa2
2069                 ecosbm=ecosb1-ecosb2
2070                 ecosgm=ecosg1-ecosg2
2071 C Diagnostics
2072 c               ecosap=ecosa1
2073 c               ecosbp=ecosb1
2074 c               ecosgp=ecosg1
2075 c               ecosam=0.0D0
2076 c               ecosbm=0.0D0
2077 c               ecosgm=0.0D0
2078 C End diagnostics
2079                 facont_hb(num_conti,i)=fcont
2080                 fprimcont=fprimcont/rij
2081 cd              facont_hb(num_conti,i)=1.0D0
2082 C Following line is for diagnostics.
2083 cd              fprimcont=0.0D0
2084                 do k=1,3
2085                   dcosb(k)=rmij*(dc_norm(k,i)-erij(k)*cosb)
2086                   dcosg(k)=rmij*(dc_norm(k,j)-erij(k)*cosg)
2087                 enddo
2088                 do k=1,3
2089                   gggp(k)=ecosbp*dcosb(k)+ecosgp*dcosg(k)
2090                   gggm(k)=ecosbm*dcosb(k)+ecosgm*dcosg(k)
2091                 enddo
2092                 gggp(1)=gggp(1)+ees0pijp*xj
2093                 gggp(2)=gggp(2)+ees0pijp*yj
2094                 gggp(3)=gggp(3)+ees0pijp*zj
2095                 gggm(1)=gggm(1)+ees0mijp*xj
2096                 gggm(2)=gggm(2)+ees0mijp*yj
2097                 gggm(3)=gggm(3)+ees0mijp*zj
2098 C Derivatives due to the contact function
2099                 gacont_hbr(1,num_conti,i)=fprimcont*xj
2100                 gacont_hbr(2,num_conti,i)=fprimcont*yj
2101                 gacont_hbr(3,num_conti,i)=fprimcont*zj
2102                 do k=1,3
2103 c
2104 c 10/24/08 cgrad and ! comments indicate the parts of the code removed
2105 c          following the change of gradient-summation algorithm.
2106 c
2107 cgrad                  ghalfp=0.5D0*gggp(k)
2108 cgrad                  ghalfm=0.5D0*gggm(k)
2109                   gacontp_hb1(k,num_conti,i)=!ghalfp
2110      &              +(ecosap*(dc_norm(k,j)-cosa*dc_norm(k,i))
2111      &              + ecosbp*(erij(k)-cosb*dc_norm(k,i)))*vbld_inv(i+1)
2112                   gacontp_hb2(k,num_conti,i)=!ghalfp
2113      &              +(ecosap*(dc_norm(k,i)-cosa*dc_norm(k,j))
2114      &              + ecosgp*(erij(k)-cosg*dc_norm(k,j)))*vbld_inv(j+1)
2115                   gacontp_hb3(k,num_conti,i)=gggp(k)
2116                   gacontm_hb1(k,num_conti,i)=!ghalfm
2117      &              +(ecosam*(dc_norm(k,j)-cosa*dc_norm(k,i))
2118      &              + ecosbm*(erij(k)-cosb*dc_norm(k,i)))*vbld_inv(i+1)
2119                   gacontm_hb2(k,num_conti,i)=!ghalfm
2120      &              +(ecosam*(dc_norm(k,i)-cosa*dc_norm(k,j))
2121      &              + ecosgm*(erij(k)-cosg*dc_norm(k,j)))*vbld_inv(j+1)
2122                   gacontm_hb3(k,num_conti,i)=gggm(k)
2123                 enddo
2124               ENDIF ! wcorr
2125               endif  ! num_conti.le.maxconts
2126             endif  ! fcont.gt.0
2127           endif    ! j.gt.i+1
2128           if (wturn3.gt.0.0d0 .or. wturn4.gt.0.0d0) then
2129             do k=1,4
2130               do l=1,3
2131                 ghalf=0.5d0*agg(l,k)
2132                 aggi(l,k)=aggi(l,k)+ghalf
2133                 aggi1(l,k)=aggi1(l,k)+agg(l,k)
2134                 aggj(l,k)=aggj(l,k)+ghalf
2135               enddo
2136             enddo
2137             if (j.eq.nres-1 .and. i.lt.j-2) then
2138               do k=1,4
2139                 do l=1,3
2140                   aggj1(l,k)=aggj1(l,k)+agg(l,k)
2141                 enddo
2142               enddo
2143             endif
2144           endif
2145 c          t_eelecij=t_eelecij+MPI_Wtime()-time00
2146       return
2147       end
2148 C-----------------------------------------------------------------------
2149       subroutine evdwpp_short(evdw1)
2150 C
2151 C Compute Evdwpp
2152 C
2153       implicit real*8 (a-h,o-z)
2154       include 'DIMENSIONS'
2155       include 'COMMON.CONTROL'
2156       include 'COMMON.IOUNITS'
2157       include 'COMMON.GEO'
2158       include 'COMMON.VAR'
2159       include 'COMMON.LOCAL'
2160       include 'COMMON.CHAIN'
2161       include 'COMMON.DERIV'
2162       include 'COMMON.INTERACT'
2163       include 'COMMON.CONTACTS'
2164       include 'COMMON.TORSION'
2165       include 'COMMON.VECTORS'
2166       include 'COMMON.FFIELD'
2167       dimension ggg(3)
2168 c 4/26/02 - AL scaling factor for 1,4 repulsive VDW interactions
2169 #ifdef MOMENT
2170       double precision scal_el /1.0d0/
2171 #else
2172       double precision scal_el /0.5d0/
2173 #endif
2174       evdw1=0.0D0
2175 C      print *,"WCHODZE"
2176 c      write (iout,*) "iatel_s_vdw",iatel_s_vdw,
2177 c     & " iatel_e_vdw",iatel_e_vdw
2178       call flush(iout)
2179       do i=iatel_s_vdw,iatel_e_vdw
2180         if (itype(i).eq.ntyp1.or. itype(i+1).eq.ntyp1) cycle
2181         dxi=dc(1,i)
2182         dyi=dc(2,i)
2183         dzi=dc(3,i)
2184         dx_normi=dc_norm(1,i)
2185         dy_normi=dc_norm(2,i)
2186         dz_normi=dc_norm(3,i)
2187         xmedi=c(1,i)+0.5d0*dxi
2188         ymedi=c(2,i)+0.5d0*dyi
2189         zmedi=c(3,i)+0.5d0*dzi
2190           xmedi=mod(xmedi,boxxsize)
2191           if (xmedi.lt.0) xmedi=xmedi+boxxsize
2192           ymedi=mod(ymedi,boxysize)
2193           if (ymedi.lt.0) ymedi=ymedi+boxysize
2194           zmedi=mod(zmedi,boxzsize)
2195           if (zmedi.lt.0) zmedi=zmedi+boxzsize
2196         num_conti=0
2197 c        write (iout,*) 'i',i,' ielstart',ielstart_vdw(i),
2198 c     &   ' ielend',ielend_vdw(i)
2199         call flush(iout)
2200         do j=ielstart_vdw(i),ielend_vdw(i)
2201           if (itype(j).eq.ntyp1 .or. itype(j+1).eq.ntyp1) cycle
2202           ind=ind+1
2203           iteli=itel(i)
2204           itelj=itel(j)
2205           if (j.eq.i+2 .and. itelj.eq.2) iteli=2
2206           aaa=app(iteli,itelj)
2207           bbb=bpp(iteli,itelj)
2208           dxj=dc(1,j)
2209           dyj=dc(2,j)
2210           dzj=dc(3,j)
2211           dx_normj=dc_norm(1,j)
2212           dy_normj=dc_norm(2,j)
2213           dz_normj=dc_norm(3,j)
2214           xj=c(1,j)+0.5D0*dxj
2215           yj=c(2,j)+0.5D0*dyj
2216           zj=c(3,j)+0.5D0*dzj
2217           xj=mod(xj,boxxsize)
2218           if (xj.lt.0) xj=xj+boxxsize
2219           yj=mod(yj,boxysize)
2220           if (yj.lt.0) yj=yj+boxysize
2221           zj=mod(zj,boxzsize)
2222           if (zj.lt.0) zj=zj+boxzsize
2223       dist_init=(xj-xmedi)**2+(yj-ymedi)**2+(zj-zmedi)**2
2224       xj_safe=xj
2225       yj_safe=yj
2226       zj_safe=zj
2227       isubchap=0
2228       do xshift=-1,1
2229       do yshift=-1,1
2230       do zshift=-1,1
2231           xj=xj_safe+xshift*boxxsize
2232           yj=yj_safe+yshift*boxysize
2233           zj=zj_safe+zshift*boxzsize
2234           dist_temp=(xj-xmedi)**2+(yj-ymedi)**2+(zj-zmedi)**2
2235           if(dist_temp.lt.dist_init) then
2236             dist_init=dist_temp
2237             xj_temp=xj
2238             yj_temp=yj
2239             zj_temp=zj
2240             isubchap=1
2241           endif
2242        enddo
2243        enddo
2244        enddo
2245        if (isubchap.eq.1) then
2246           xj=xj_temp-xmedi
2247           yj=yj_temp-ymedi
2248           zj=zj_temp-zmedi
2249        else
2250           xj=xj_safe-xmedi
2251           yj=yj_safe-ymedi
2252           zj=zj_safe-zmedi
2253        endif
2254           rij=xj*xj+yj*yj+zj*zj
2255           rrmij=1.0D0/rij
2256           rij=dsqrt(rij)
2257           sss=sscale(rij/rpp(iteli,itelj))
2258             sssgrad=sscagrad(rij/rpp(iteli,itelj))
2259           if (sss.gt.0.0d0) then
2260             rmij=1.0D0/rij
2261             r3ij=rrmij*rmij
2262             r6ij=r3ij*r3ij
2263             ev1=aaa*r6ij*r6ij
2264 c 4/26/02 - AL scaling down 1,4 repulsive VDW interactions
2265             if (j.eq.i+2) ev1=scal_el*ev1
2266             ev2=bbb*r6ij
2267             evdwij=ev1+ev2
2268             if (energy_dec) then
2269               write (iout,'(a6,2i5,0pf7.3,f7.3)') 'evdw1',i,j,evdwij,sss
2270             endif
2271             evdw1=evdw1+evdwij*sss
2272 C
2273 C Calculate contributions to the Cartesian gradient.
2274 C
2275             facvdw=-6*rrmij*(ev1+evdwij)*sss
2276           ggg(1)=facvdw*xj+sssgrad*rmij*evdwij*xj/rpp(iteli,itelj)
2277           ggg(2)=facvdw*yj+sssgrad*rmij*evdwij*yj/rpp(iteli,itelj)
2278           ggg(3)=facvdw*zj+sssgrad*rmij*evdwij*zj/rpp(iteli,itelj)
2279 C            ggg(1)=facvdw*xj
2280 C            ggg(2)=facvdw*yj
2281 C            ggg(3)=facvdw*zj
2282             do k=1,3
2283               gvdwpp(k,j)=gvdwpp(k,j)+ggg(k)
2284               gvdwpp(k,i)=gvdwpp(k,i)-ggg(k)
2285             enddo
2286           endif
2287         enddo ! j
2288       enddo   ! i
2289       return
2290       end
2291 C-----------------------------------------------------------------------------
2292       subroutine escp_long(evdw2,evdw2_14)
2293 C
2294 C This subroutine calculates the excluded-volume interaction energy between
2295 C peptide-group centers and side chains and its gradient in virtual-bond and
2296 C side-chain vectors.
2297 C
2298       implicit real*8 (a-h,o-z)
2299       include 'DIMENSIONS'
2300       include 'COMMON.GEO'
2301       include 'COMMON.VAR'
2302       include 'COMMON.LOCAL'
2303       include 'COMMON.CHAIN'
2304       include 'COMMON.DERIV'
2305       include 'COMMON.INTERACT'
2306       include 'COMMON.FFIELD'
2307       include 'COMMON.IOUNITS'
2308       include 'COMMON.CONTROL'
2309       dimension ggg(3)
2310       evdw2=0.0D0
2311       evdw2_14=0.0d0
2312 CD        print '(a)','Enter ESCP KURWA'
2313 cd    write (iout,*) 'iatscp_s=',iatscp_s,' iatscp_e=',iatscp_e
2314       do i=iatscp_s,iatscp_e
2315         if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1) cycle
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          xi=mod(xi,boxxsize)
2321           if (xi.lt.0) xi=xi+boxxsize
2322           yi=mod(yi,boxysize)
2323           if (yi.lt.0) yi=yi+boxysize
2324           zi=mod(zi,boxzsize)
2325           if (zi.lt.0) zi=zi+boxzsize
2326         do iint=1,nscp_gr(i)
2327
2328         do j=iscpstart(i,iint),iscpend(i,iint)
2329           itypj=itype(j)
2330           if (itypj.eq.ntyp1) cycle
2331 C Uncomment following three lines for SC-p interactions
2332 c         xj=c(1,nres+j)-xi
2333 c         yj=c(2,nres+j)-yi
2334 c         zj=c(3,nres+j)-zi
2335 C Uncomment following three lines for Ca-p interactions
2336           xj=c(1,j)
2337           yj=c(2,j)
2338           zj=c(3,j)
2339       dist_init=(xj-xi)**2+(yj-yi)**2+(zj-zi)**2
2340       xj_safe=xj
2341       yj_safe=yj
2342       zj_safe=zj
2343       subchap=0
2344       do xshift=-1,1
2345       do yshift=-1,1
2346       do zshift=-1,1
2347           xj=xj_safe+xshift*boxxsize
2348           yj=yj_safe+yshift*boxysize
2349           zj=zj_safe+zshift*boxzsize
2350           dist_temp=(xj-xi)**2+(yj-yi)**2+(zj-zi)**2
2351           if(dist_temp.lt.dist_init) then
2352             dist_init=dist_temp
2353             xj_temp=xj
2354             yj_temp=yj
2355             zj_temp=zj
2356             subchap=1
2357           endif
2358        enddo
2359        enddo
2360        enddo
2361        if (subchap.eq.1) then
2362           xj=xj_temp-xi
2363           yj=yj_temp-yi
2364           zj=zj_temp-zi
2365        else
2366           xj=xj_safe-xi
2367           yj=yj_safe-yi
2368           zj=zj_safe-zi
2369        endif
2370
2371           rrij=1.0D0/(xj*xj+yj*yj+zj*zj)
2372
2373           sss=sscale(1.0d0/(dsqrt(rrij)*rscp(itypj,iteli)))
2374           sssgrad=sscagrad(1.0d0/(dsqrt(rrij)*rscp(itypj,iteli)))
2375           if (sss.lt.1.0d0) then
2376             fac=rrij**expon2
2377             e1=fac*fac*aad(itypj,iteli)
2378             e2=fac*bad(itypj,iteli)
2379             if (iabs(j-i) .le. 2) then
2380               e1=scal14*e1
2381               e2=scal14*e2
2382               evdw2_14=evdw2_14+(e1+e2)*(1.0d0-sss)
2383             endif
2384             evdwij=e1+e2
2385             evdw2=evdw2+evdwij*(1.0d0-sss)
2386             if (energy_dec) write (iout,'(a6,2i5,2(0pf7.3))')
2387      &          'evdw2',i,j,sss,evdwij
2388 C
2389 C Calculate contributions to the gradient in the virtual-bond and SC vectors.
2390 C
2391
2392             fac=-(evdwij+e1)*rrij*(1.0d0-sss)
2393             fac=fac-(evdwij)*sssgrad*dsqrt(rrij)/rscp(itypj,iteli)
2394             ggg(1)=xj*fac
2395             ggg(2)=yj*fac
2396             ggg(3)=zj*fac
2397 C Uncomment following three lines for SC-p interactions
2398 c           do k=1,3
2399 c             gradx_scp(k,j)=gradx_scp(k,j)+ggg(k)
2400 c           enddo
2401 C Uncomment following line for SC-p interactions
2402 c             gradx_scp(k,j)=gradx_scp(k,j)+ggg(k)
2403             do k=1,3
2404               gvdwc_scpp(k,i)=gvdwc_scpp(k,i)-ggg(k)
2405               gvdwc_scp(k,j)=gvdwc_scp(k,j)+ggg(k)
2406             enddo
2407           endif
2408         enddo
2409
2410         enddo ! iint
2411       enddo ! i
2412       do i=1,nct
2413         do j=1,3
2414           gvdwc_scp(j,i)=expon*gvdwc_scp(j,i)
2415           gvdwc_scpp(j,i)=expon*gvdwc_scpp(j,i)
2416           gradx_scp(j,i)=expon*gradx_scp(j,i)
2417         enddo
2418       enddo
2419 C******************************************************************************
2420 C
2421 C                              N O T E !!!
2422 C
2423 C To save time the factor EXPON has been extracted from ALL components
2424 C of GVDWC and GRADX. Remember to multiply them by this factor before further
2425 C use!
2426 C
2427 C******************************************************************************
2428       return
2429       end
2430 C-----------------------------------------------------------------------------
2431       subroutine escp_short(evdw2,evdw2_14)
2432 C
2433 C This subroutine calculates the excluded-volume interaction energy between
2434 C peptide-group centers and side chains and its gradient in virtual-bond and
2435 C side-chain vectors.
2436 C
2437       implicit real*8 (a-h,o-z)
2438       include 'DIMENSIONS'
2439       include 'COMMON.GEO'
2440       include 'COMMON.VAR'
2441       include 'COMMON.LOCAL'
2442       include 'COMMON.CHAIN'
2443       include 'COMMON.DERIV'
2444       include 'COMMON.INTERACT'
2445       include 'COMMON.FFIELD'
2446       include 'COMMON.IOUNITS'
2447       include 'COMMON.CONTROL'
2448       dimension ggg(3)
2449       evdw2=0.0D0
2450       evdw2_14=0.0d0
2451 cd    print '(a)','Enter ESCP'
2452 cd    write (iout,*) 'iatscp_s=',iatscp_s,' iatscp_e=',iatscp_e
2453       do i=iatscp_s,iatscp_e
2454         if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1) cycle
2455         iteli=itel(i)
2456         xi=0.5D0*(c(1,i)+c(1,i+1))
2457         yi=0.5D0*(c(2,i)+c(2,i+1))
2458         zi=0.5D0*(c(3,i)+c(3,i+1))
2459          xi=mod(xi,boxxsize)
2460           if (xi.lt.0) xi=xi+boxxsize
2461           yi=mod(yi,boxysize)
2462           if (yi.lt.0) yi=yi+boxysize
2463           zi=mod(zi,boxzsize)
2464           if (zi.lt.0) zi=zi+boxzsize
2465
2466         do iint=1,nscp_gr(i)
2467
2468         do j=iscpstart(i,iint),iscpend(i,iint)
2469           itypj=itype(j)
2470           if (itypj.eq.ntyp1) cycle
2471 C Uncomment following three lines for SC-p interactions
2472 c         xj=c(1,nres+j)-xi
2473 c         yj=c(2,nres+j)-yi
2474 c         zj=c(3,nres+j)-zi
2475 C Uncomment following three lines for Ca-p interactions
2476           xj=c(1,j)
2477           yj=c(2,j)
2478           zj=c(3,j)
2479       dist_init=(xj-xi)**2+(yj-yi)**2+(zj-zi)**2
2480       xj_safe=xj
2481       yj_safe=yj
2482       zj_safe=zj
2483       subchap=0
2484       do xshift=-1,1
2485       do yshift=-1,1
2486       do zshift=-1,1
2487           xj=xj_safe+xshift*boxxsize
2488           yj=yj_safe+yshift*boxysize
2489           zj=zj_safe+zshift*boxzsize
2490           dist_temp=(xj-xi)**2+(yj-yi)**2+(zj-zi)**2
2491           if(dist_temp.lt.dist_init) then
2492             dist_init=dist_temp
2493             xj_temp=xj
2494             yj_temp=yj
2495             zj_temp=zj
2496             subchap=1
2497           endif
2498        enddo
2499        enddo
2500        enddo
2501        if (subchap.eq.1) then
2502           xj=xj_temp-xi
2503           yj=yj_temp-yi
2504           zj=zj_temp-zi
2505        else
2506           xj=xj_safe-xi
2507           yj=yj_safe-yi
2508           zj=zj_safe-zi
2509        endif
2510           rrij=1.0D0/(xj*xj+yj*yj+zj*zj)
2511           sss=sscale(1.0d0/(dsqrt(rrij)*rscp(itypj,iteli)))
2512           sssgrad=sscagrad(1.0d0/(dsqrt(rrij)*rscp(itypj,iteli)))
2513           if (sss.gt.0.0d0) then
2514
2515             fac=rrij**expon2
2516             e1=fac*fac*aad(itypj,iteli)
2517             e2=fac*bad(itypj,iteli)
2518             if (iabs(j-i) .le. 2) then
2519               e1=scal14*e1
2520               e2=scal14*e2
2521               evdw2_14=evdw2_14+(e1+e2)*sss
2522             endif
2523             evdwij=e1+e2
2524             evdw2=evdw2+evdwij*sss
2525             if (energy_dec) write (iout,'(a6,2i5,2(0pf7.3))')
2526      &          'evdw2',i,j,sss,evdwij
2527 C
2528 C Calculate contributions to the gradient in the virtual-bond and SC vectors.
2529 C
2530             fac=-(evdwij+e1)*rrij*sss
2531             fac=fac+(evdwij)*sssgrad*dsqrt(rrij)/rscp(itypj,iteli)
2532             ggg(1)=xj*fac
2533             ggg(2)=yj*fac
2534             ggg(3)=zj*fac
2535 C Uncomment following three lines for SC-p interactions
2536 c           do k=1,3
2537 c             gradx_scp(k,j)=gradx_scp(k,j)+ggg(k)
2538 c           enddo
2539 C Uncomment following line for SC-p interactions
2540 c             gradx_scp(k,j)=gradx_scp(k,j)+ggg(k)
2541             do k=1,3
2542               gvdwc_scpp(k,i)=gvdwc_scpp(k,i)-ggg(k)
2543               gvdwc_scp(k,j)=gvdwc_scp(k,j)+ggg(k)
2544             enddo
2545           endif
2546         enddo
2547
2548         enddo ! iint
2549       enddo ! i
2550       do i=1,nct
2551         do j=1,3
2552           gvdwc_scp(j,i)=expon*gvdwc_scp(j,i)
2553           gvdwc_scpp(j,i)=expon*gvdwc_scpp(j,i)
2554           gradx_scp(j,i)=expon*gradx_scp(j,i)
2555         enddo
2556       enddo
2557 C******************************************************************************
2558 C
2559 C                              N O T E !!!
2560 C
2561 C To save time the factor EXPON has been extracted from ALL components
2562 C of GVDWC and GRADX. Remember to multiply them by this factor before further
2563 C use!
2564 C
2565 C******************************************************************************
2566       return
2567       end