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