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

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