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