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