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