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