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