source/unres/src-HCD-5D/energy_p_new-sep_barrier.F

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