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