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       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       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       include 'COMMON.CONTACTS'
1372       include 'COMMON.TORSION'
1373       include 'COMMON.VECTORS'
1374       include 'COMMON.FFIELD'
1375       include 'COMMON.TIME1'
1376       include 'COMMON.SHIELD'
1377       dimension ggg(3),gggp(3),gggm(3),erij(3),dcosb(3),dcosg(3),
1378      &          erder(3,3),uryg(3,3),urzg(3,3),vryg(3,3),vrzg(3,3)
1379       double precision acipa(2,2),agg(3,4),aggi(3,4),aggi1(3,4),
1380      &    aggj(3,4),aggj1(3,4),a_temp(2,2),muij(4)
1381       common /locel/ a_temp,agg,aggi,aggi1,aggj,aggj1,a22,a23,a32,a33,
1382      &    dxi,dyi,dzi,dx_normi,dy_normi,dz_normi,xmedi,ymedi,zmedi,
1383      &    num_conti,j1,j2
1384 c 4/26/02 - AL scaling factor for 1,4 repulsive VDW interactions
1385 #ifdef MOMENT
1386       double precision scal_el /1.0d0/
1387 #else
1388       double precision scal_el /0.5d0/
1389 #endif
1390 C 12/13/98
1391 C 13-go grudnia roku pamietnego...
1392       double precision unmat(3,3) /1.0d0,0.0d0,0.0d0,
1393      &                   0.0d0,1.0d0,0.0d0,
1394      &                   0.0d0,0.0d0,1.0d0/
1395 cd      write(iout,*) 'In EELEC'
1396 cd      do i=1,nloctyp
1397 cd        write(iout,*) 'Type',i
1398 cd        write(iout,*) 'B1',B1(:,i)
1399 cd        write(iout,*) 'B2',B2(:,i)
1400 cd        write(iout,*) 'CC',CC(:,:,i)
1401 cd        write(iout,*) 'DD',DD(:,:,i)
1402 cd        write(iout,*) 'EE',EE(:,:,i)
1403 cd      enddo
1404 cd      call check_vecgrad
1405 cd      stop
1406 C      print *,"WCHODZE3"
1407       if (icheckgrad.eq.1) then
1408         do i=1,nres-1
1409           fac=1.0d0/dsqrt(scalar(dc(1,i),dc(1,i)))
1410           do k=1,3
1411             dc_norm(k,i)=dc(k,i)*fac
1412           enddo
1413 c          write (iout,*) 'i',i,' fac',fac
1414         enddo
1415       endif
1416       if (wel_loc.gt.0.0d0 .or. wcorr4.gt.0.0d0 .or. wcorr5.gt.0.0d0
1417      &    .or. wcorr6.gt.0.0d0 .or. wturn3.gt.0.0d0 .or.
1418      &    wturn4.gt.0.0d0 .or. wturn6.gt.0.0d0) then
1419 c        call vec_and_deriv
1420 #ifdef TIMING
1421         time01=MPI_Wtime()
1422 #endif
1423         call set_matrices
1424 #ifdef TIMING
1425         time_mat=time_mat+MPI_Wtime()-time01
1426 #endif
1427       endif
1428 cd      do i=1,nres-1
1429 cd        write (iout,*) 'i=',i
1430 cd        do k=1,3
1431 cd        write (iout,'(i5,2f10.5)') k,uy(k,i),uz(k,i)
1432 cd        enddo
1433 cd        do k=1,3
1434 cd          write (iout,'(f10.5,2x,3f10.5,2x,3f10.5)')
1435 cd     &     uz(k,i),(uzgrad(k,l,1,i),l=1,3),(uzgrad(k,l,2,i),l=1,3)
1436 cd        enddo
1437 cd      enddo
1438       t_eelecij=0.0d0
1439       ees=0.0D0
1440       evdw1=0.0D0
1441       eel_loc=0.0d0
1442       eello_turn3=0.0d0
1443       eello_turn4=0.0d0
1444       ind=0
1445       do i=1,nres
1446         num_cont_hb(i)=0
1447       enddo
1448 cd      print '(a)','Enter EELEC'
1449 cd      write (iout,*) 'iatel_s=',iatel_s,' iatel_e=',iatel_e
1450       do i=1,nres
1451         gel_loc_loc(i)=0.0d0
1452         gcorr_loc(i)=0.0d0
1453       enddo
1454 c
1455 c
1456 c 9/27/08 AL Split the interaction loop to ensure load balancing of turn terms
1457 C
1458 C Loop over i,i+2 and i,i+3 pairs of the peptide groups
1459 C
1460       do i=iturn3_start,iturn3_end
1461         if (itype(i).eq.ntyp1.or. itype(i+1).eq.ntyp1
1462      &  .or. itype(i+2).eq.ntyp1 .or. itype(i+3).eq.ntyp1
1463 C     &  .or. itype(i-1).eq.ntyp1
1464 C     &  .or. itype(i+4).eq.ntyp1
1465      &   ) cycle
1466         dxi=dc(1,i)
1467         dyi=dc(2,i)
1468         dzi=dc(3,i)
1469         dx_normi=dc_norm(1,i)
1470         dy_normi=dc_norm(2,i)
1471         dz_normi=dc_norm(3,i)
1472         xmedi=c(1,i)+0.5d0*dxi
1473         ymedi=c(2,i)+0.5d0*dyi
1474         zmedi=c(3,i)+0.5d0*dzi
1475           xmedi=mod(xmedi,boxxsize)
1476           if (xmedi.lt.0) xmedi=xmedi+boxxsize
1477           ymedi=mod(ymedi,boxysize)
1478           if (ymedi.lt.0) ymedi=ymedi+boxysize
1479           zmedi=mod(zmedi,boxzsize)
1480           if (zmedi.lt.0) zmedi=zmedi+boxzsize
1481         num_conti=0
1482         call eelecij_scale(i,i+2,ees,evdw1,eel_loc)
1483         if (wturn3.gt.0.0d0) call eturn3(i,eello_turn3)
1484         num_cont_hb(i)=num_conti
1485       enddo
1486       do i=iturn4_start,iturn4_end
1487         if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1
1488      &    .or. itype(i+3).eq.ntyp1
1489      &    .or. itype(i+4).eq.ntyp1
1490 C     &    .or. itype(i+5).eq.ntyp1
1491 C     &    .or. itype(i-1).eq.ntyp1
1492      &    ) cycle
1493         dxi=dc(1,i)
1494         dyi=dc(2,i)
1495         dzi=dc(3,i)
1496         dx_normi=dc_norm(1,i)
1497         dy_normi=dc_norm(2,i)
1498         dz_normi=dc_norm(3,i)
1499         xmedi=c(1,i)+0.5d0*dxi
1500         ymedi=c(2,i)+0.5d0*dyi
1501         zmedi=c(3,i)+0.5d0*dzi
1502           xmedi=mod(xmedi,boxxsize)
1503           if (xmedi.lt.0) xmedi=xmedi+boxxsize
1504           ymedi=mod(ymedi,boxysize)
1505           if (ymedi.lt.0) ymedi=ymedi+boxysize
1506           zmedi=mod(zmedi,boxzsize)
1507           if (zmedi.lt.0) zmedi=zmedi+boxzsize
1508         num_conti=num_cont_hb(i)
1509         call eelecij_scale(i,i+3,ees,evdw1,eel_loc)
1510         if (wturn4.gt.0.0d0 .and. itype(i+2).ne.ntyp1)
1511      &    call eturn4(i,eello_turn4)
1512         num_cont_hb(i)=num_conti
1513       enddo   ! i
1514 c
1515 c Loop over all pairs of interacting peptide groups except i,i+2 and i,i+3
1516 c
1517       do i=iatel_s,iatel_e
1518         if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1
1519 C     &  .or. itype(i+2).eq.ntyp1
1520 C     &  .or. itype(i-1).eq.ntyp1
1521      &) cycle
1522         dxi=dc(1,i)
1523         dyi=dc(2,i)
1524         dzi=dc(3,i)
1525         dx_normi=dc_norm(1,i)
1526         dy_normi=dc_norm(2,i)
1527         dz_normi=dc_norm(3,i)
1528         xmedi=c(1,i)+0.5d0*dxi
1529         ymedi=c(2,i)+0.5d0*dyi
1530         zmedi=c(3,i)+0.5d0*dzi
1531           xmedi=mod(xmedi,boxxsize)
1532           if (xmedi.lt.0) xmedi=xmedi+boxxsize
1533           ymedi=mod(ymedi,boxysize)
1534           if (ymedi.lt.0) ymedi=ymedi+boxysize
1535           zmedi=mod(zmedi,boxzsize)
1536           if (zmedi.lt.0) zmedi=zmedi+boxzsize
1537 c        write (iout,*) 'i',i,' ielstart',ielstart(i),' ielend',ielend(i)
1538         num_conti=num_cont_hb(i)
1539         do j=ielstart(i),ielend(i)
1540           if (itype(j).eq.ntyp1 .or. itype(j+1).eq.ntyp1
1541 C     & .or.itype(j+2).eq.ntyp1
1542 C     & .or.itype(j-1).eq.ntyp1
1543      &) cycle
1544           call eelecij_scale(i,j,ees,evdw1,eel_loc)
1545         enddo ! j
1546         num_cont_hb(i)=num_conti
1547       enddo   ! i
1548 c      write (iout,*) "Number of loop steps in EELEC:",ind
1549 cd      do i=1,nres
1550 cd        write (iout,'(i3,3f10.5,5x,3f10.5)')
1551 cd     &     i,(gel_loc(k,i),k=1,3),gel_loc_loc(i)
1552 cd      enddo
1553 c 12/7/99 Adam eello_turn3 will be considered as a separate energy term
1554 ccc      eel_loc=eel_loc+eello_turn3
1555 cd      print *,"Processor",fg_rank," t_eelecij",t_eelecij
1556       return
1557       end
1558 C-------------------------------------------------------------------------------
1559       subroutine eelecij_scale(i,j,ees,evdw1,eel_loc)
1560       implicit real*8 (a-h,o-z)
1561       include 'DIMENSIONS'
1562 #ifdef MPI
1563       include "mpif.h"
1564 #endif
1565       include 'COMMON.CONTROL'
1566       include 'COMMON.IOUNITS'
1567       include 'COMMON.GEO'
1568       include 'COMMON.VAR'
1569       include 'COMMON.LOCAL'
1570       include 'COMMON.CHAIN'
1571       include 'COMMON.DERIV'
1572       include 'COMMON.INTERACT'
1573       include 'COMMON.CONTACTS'
1574       include 'COMMON.TORSION'
1575       include 'COMMON.VECTORS'
1576       include 'COMMON.FFIELD'
1577       include 'COMMON.TIME1'
1578       include 'COMMON.SHIELD'
1579       integer xshift,yshift,zshift
1580       dimension ggg(3),gggp(3),gggm(3),erij(3),dcosb(3),dcosg(3),
1581      &          erder(3,3),uryg(3,3),urzg(3,3),vryg(3,3),vrzg(3,3)
1582       double precision acipa(2,2),agg(3,4),aggi(3,4),aggi1(3,4),
1583      &    aggj(3,4),aggj1(3,4),a_temp(2,2),muij(4),gmuij1(4),gmuji1(4),
1584      &    gmuij2(4),gmuji2(4)
1585       common /locel/ a_temp,agg,aggi,aggi1,aggj,aggj1,a22,a23,a32,a33,
1586      &    dxi,dyi,dzi,dx_normi,dy_normi,dz_normi,xmedi,ymedi,zmedi,
1587      &    num_conti,j1,j2
1588 c 4/26/02 - AL scaling factor for 1,4 repulsive VDW interactions
1589 #ifdef MOMENT
1590       double precision scal_el /1.0d0/
1591 #else
1592       double precision scal_el /0.5d0/
1593 #endif
1594 C 12/13/98
1595 C 13-go grudnia roku pamietnego...
1596       double precision unmat(3,3) /1.0d0,0.0d0,0.0d0,
1597      &                   0.0d0,1.0d0,0.0d0,
1598      &                   0.0d0,0.0d0,1.0d0/
1599 c          time00=MPI_Wtime()
1600 cd      write (iout,*) "eelecij",i,j
1601 C      print *,"WCHODZE2"
1602           ind=ind+1
1603           iteli=itel(i)
1604           itelj=itel(j)
1605           if (j.eq.i+2 .and. itelj.eq.2) iteli=2
1606           aaa=app(iteli,itelj)
1607           bbb=bpp(iteli,itelj)
1608           ael6i=ael6(iteli,itelj)
1609           ael3i=ael3(iteli,itelj)
1610           dxj=dc(1,j)
1611           dyj=dc(2,j)
1612           dzj=dc(3,j)
1613           dx_normj=dc_norm(1,j)
1614           dy_normj=dc_norm(2,j)
1615           dz_normj=dc_norm(3,j)
1616           xj=c(1,j)+0.5D0*dxj
1617           yj=c(2,j)+0.5D0*dyj
1618           zj=c(3,j)+0.5D0*dzj
1619          xj=mod(xj,boxxsize)
1620           if (xj.lt.0) xj=xj+boxxsize
1621           yj=mod(yj,boxysize)
1622           if (yj.lt.0) yj=yj+boxysize
1623           zj=mod(zj,boxzsize)
1624           if (zj.lt.0) zj=zj+boxzsize
1625       dist_init=(xj-xmedi)**2+(yj-ymedi)**2+(zj-zmedi)**2
1626       xj_safe=xj
1627       yj_safe=yj
1628       zj_safe=zj
1629       isubchap=0
1630       do xshift=-1,1
1631       do yshift=-1,1
1632       do zshift=-1,1
1633           xj=xj_safe+xshift*boxxsize
1634           yj=yj_safe+yshift*boxysize
1635           zj=zj_safe+zshift*boxzsize
1636           dist_temp=(xj-xmedi)**2+(yj-ymedi)**2+(zj-zmedi)**2
1637           if(dist_temp.lt.dist_init) then
1638             dist_init=dist_temp
1639             xj_temp=xj
1640             yj_temp=yj
1641             zj_temp=zj
1642             isubchap=1
1643           endif
1644        enddo
1645        enddo
1646        enddo
1647        if (isubchap.eq.1) then
1648           xj=xj_temp-xmedi
1649           yj=yj_temp-ymedi
1650           zj=zj_temp-zmedi
1651        else
1652           xj=xj_safe-xmedi
1653           yj=yj_safe-ymedi
1654           zj=zj_safe-zmedi
1655        endif
1656
1657           rij=xj*xj+yj*yj+zj*zj
1658           rrmij=1.0D0/rij
1659           rij=dsqrt(rij)
1660           rmij=1.0D0/rij
1661 c For extracting the short-range part of Evdwpp
1662           sss=sscale(rij/rpp(iteli,itelj))
1663           sssgrad=sscagrad(rij/rpp(iteli,itelj))
1664           r3ij=rrmij*rmij
1665           r6ij=r3ij*r3ij
1666           cosa=dx_normi*dx_normj+dy_normi*dy_normj+dz_normi*dz_normj
1667           cosb=(xj*dx_normi+yj*dy_normi+zj*dz_normi)*rmij
1668           cosg=(xj*dx_normj+yj*dy_normj+zj*dz_normj)*rmij
1669           fac=cosa-3.0D0*cosb*cosg
1670           ev1=aaa*r6ij*r6ij
1671 c 4/26/02 - AL scaling down 1,4 repulsive VDW interactions
1672           if (j.eq.i+2) ev1=scal_el*ev1
1673           ev2=bbb*r6ij
1674           fac3=ael6i*r6ij
1675           fac4=ael3i*r3ij
1676           evdwij=ev1+ev2
1677           if (shield_mode.eq.0) then
1678           fac_shield(i)=1.0
1679           fac_shield(j)=1.0
1680           endif
1681           el1=fac3*(4.0D0+fac*fac-3.0D0*(cosb*cosb+cosg*cosg))
1682           el2=fac4*fac
1683           el1=el1*fac_shield(i)**2*fac_shield(j)**2
1684           el2=el2*fac_shield(i)**2*fac_shield(j)**2
1685           eesij=el1+el2
1686 C 12/26/95 - for the evaluation of multi-body H-bonding interactions
1687           ees0ij=4.0D0+fac*fac-3.0D0*(cosb*cosb+cosg*cosg)
1688           ees=ees+eesij
1689           evdw1=evdw1+evdwij*(1.0d0-sss)
1690 cd          write(iout,'(2(2i3,2x),7(1pd12.4)/2(3(1pd12.4),5x)/)')
1691 cd     &      iteli,i,itelj,j,aaa,bbb,ael6i,ael3i,
1692 cd     &      1.0D0/dsqrt(rrmij),evdwij,eesij,
1693 cd     &      xmedi,ymedi,zmedi,xj,yj,zj
1694
1695           if (energy_dec) then
1696               write (iout,'(a6,2i5,0pf7.3,f7.3)') 'evdw1',i,j,evdwij,sss
1697               write (iout,'(a6,2i5,0pf7.3)') 'ees',i,j,eesij
1698           endif
1699
1700 C
1701 C Calculate contributions to the Cartesian gradient.
1702 C
1703 #ifdef SPLITELE
1704           facvdw=-6*rrmij*(ev1+evdwij)*(1.0d0-sss)
1705           facel=-3*rrmij*(el1+eesij)
1706           fac1=fac
1707           erij(1)=xj*rmij
1708           erij(2)=yj*rmij
1709           erij(3)=zj*rmij
1710 *
1711 * Radial derivatives. First process both termini of the fragment (i,j)
1712 *
1713           ggg(1)=facel*xj
1714           ggg(2)=facel*yj
1715           ggg(3)=facel*zj
1716           if ((fac_shield(i).gt.0).and.(fac_shield(j).gt.0).and.
1717      &  (shield_mode.gt.0)) then
1718 C          print *,i,j
1719           do ilist=1,ishield_list(i)
1720            iresshield=shield_list(ilist,i)
1721            do k=1,3
1722            rlocshield=grad_shield_side(k,ilist,i)*eesij/fac_shield(i)
1723      &      *2.0
1724            gshieldx(k,iresshield)=gshieldx(k,iresshield)+
1725      &              rlocshield
1726      & +grad_shield_loc(k,ilist,i)*eesij/fac_shield(i)*2.0
1727             gshieldc(k,iresshield-1)=gshieldc(k,iresshield-1)+rlocshield
1728 C           gshieldc_loc(k,iresshield)=gshieldc_loc(k,iresshield)
1729 C     & +grad_shield_loc(k,ilist,i)*eesij/fac_shield(i)
1730 C             if (iresshield.gt.i) then
1731 C               do ishi=i+1,iresshield-1
1732 C                gshieldc(k,ishi)=gshieldc(k,ishi)+rlocshield
1733 C     & +grad_shield_loc(k,ilist,i)*eesij/fac_shield(i)
1734 C
1735 C              enddo
1736 C             else
1737 C               do ishi=iresshield,i
1738 C                gshieldc(k,ishi)=gshieldc(k,ishi)-rlocshield
1739 C     & -grad_shield_loc(k,ilist,i)*eesij/fac_shield(i)
1740 C
1741 C               enddo
1742 C              endif
1743            enddo
1744           enddo
1745           do ilist=1,ishield_list(j)
1746            iresshield=shield_list(ilist,j)
1747            do k=1,3
1748            rlocshield=grad_shield_side(k,ilist,j)*eesij/fac_shield(j)
1749      &     *2.0
1750            gshieldx(k,iresshield)=gshieldx(k,iresshield)+
1751      &              rlocshield
1752      & +grad_shield_loc(k,ilist,j)*eesij/fac_shield(j)*2.0
1753            gshieldc(k,iresshield-1)=gshieldc(k,iresshield-1)+rlocshield
1754            enddo
1755           enddo
1756
1757           do k=1,3
1758             gshieldc(k,i)=gshieldc(k,i)+
1759      &              grad_shield(k,i)*eesij/fac_shield(i)*2.0
1760             gshieldc(k,j)=gshieldc(k,j)+
1761      &              grad_shield(k,j)*eesij/fac_shield(j)*2.0
1762             gshieldc(k,i-1)=gshieldc(k,i-1)+
1763      &              grad_shield(k,i)*eesij/fac_shield(i)*2.0
1764             gshieldc(k,j-1)=gshieldc(k,j-1)+
1765      &              grad_shield(k,j)*eesij/fac_shield(j)*2.0
1766
1767            enddo
1768            endif
1769
1770 c          do k=1,3
1771 c            ghalf=0.5D0*ggg(k)
1772 c            gelc(k,i)=gelc(k,i)+ghalf
1773 c            gelc(k,j)=gelc(k,j)+ghalf
1774 c          enddo
1775 c 9/28/08 AL Gradient compotents will be summed only at the end
1776           do k=1,3
1777             gelc_long(k,j)=gelc_long(k,j)+ggg(k)
1778             gelc_long(k,i)=gelc_long(k,i)-ggg(k)
1779           enddo
1780 *
1781 * Loop over residues i+1 thru j-1.
1782 *
1783 cgrad          do k=i+1,j-1
1784 cgrad            do l=1,3
1785 cgrad              gelc(l,k)=gelc(l,k)+ggg(l)
1786 cgrad            enddo
1787 cgrad          enddo
1788           ggg(1)=facvdw*xj-sssgrad*rmij*evdwij*xj/rpp(iteli,itelj)
1789           ggg(2)=facvdw*yj-sssgrad*rmij*evdwij*yj/rpp(iteli,itelj)
1790           ggg(3)=facvdw*zj-sssgrad*rmij*evdwij*zj/rpp(iteli,itelj)
1791 c          do k=1,3
1792 c            ghalf=0.5D0*ggg(k)
1793 c            gvdwpp(k,i)=gvdwpp(k,i)+ghalf
1794 c            gvdwpp(k,j)=gvdwpp(k,j)+ghalf
1795 c          enddo
1796 c 9/28/08 AL Gradient compotents will be summed only at the end
1797           do k=1,3
1798             gvdwpp(k,j)=gvdwpp(k,j)+ggg(k)
1799             gvdwpp(k,i)=gvdwpp(k,i)-ggg(k)
1800           enddo
1801 *
1802 * Loop over residues i+1 thru j-1.
1803 *
1804 cgrad          do k=i+1,j-1
1805 cgrad            do l=1,3
1806 cgrad              gvdwpp(l,k)=gvdwpp(l,k)+ggg(l)
1807 cgrad            enddo
1808 cgrad          enddo
1809 #else
1810           facvdw=ev1+evdwij*(1.0d0-sss)
1811           facel=el1+eesij
1812           fac1=fac
1813           fac=-3*rrmij*(facvdw+facvdw+facel)
1814           erij(1)=xj*rmij
1815           erij(2)=yj*rmij
1816           erij(3)=zj*rmij
1817 *
1818 * Radial derivatives. First process both termini of the fragment (i,j)
1819 *
1820           ggg(1)=fac*xj
1821           ggg(2)=fac*yj
1822           ggg(3)=fac*zj
1823 c          do k=1,3
1824 c            ghalf=0.5D0*ggg(k)
1825 c            gelc(k,i)=gelc(k,i)+ghalf
1826 c            gelc(k,j)=gelc(k,j)+ghalf
1827 c          enddo
1828 c 9/28/08 AL Gradient compotents will be summed only at the end
1829           do k=1,3
1830             gelc_long(k,j)=gelc(k,j)+ggg(k)
1831             gelc_long(k,i)=gelc(k,i)-ggg(k)
1832           enddo
1833 *
1834 * Loop over residues i+1 thru j-1.
1835 *
1836 cgrad          do k=i+1,j-1
1837 cgrad            do l=1,3
1838 cgrad              gelc(l,k)=gelc(l,k)+ggg(l)
1839 cgrad            enddo
1840 cgrad          enddo
1841 c 9/28/08 AL Gradient compotents will be summed only at the end
1842 C          ggg(1)=facvdw*xj
1843 C          ggg(2)=facvdw*yj
1844 C          ggg(3)=facvdw*zj
1845           ggg(1)=facvdw*xj-sssgrad*rmij*evdwij*xj/rpp(iteli,itelj)
1846           ggg(2)=facvdw*yj-sssgrad*rmij*evdwij*yj/rpp(iteli,itelj)
1847           ggg(3)=facvdw*zj-sssgrad*rmij*evdwij*zj/rpp(iteli,itelj)
1848           do k=1,3
1849             gvdwpp(k,j)=gvdwpp(k,j)+ggg(k)
1850             gvdwpp(k,i)=gvdwpp(k,i)-ggg(k)
1851           enddo
1852 #endif
1853 *
1854 * Angular part
1855 *
1856           ecosa=2.0D0*fac3*fac1+fac4
1857           fac4=-3.0D0*fac4
1858           fac3=-6.0D0*fac3
1859           ecosb=(fac3*(fac1*cosg+cosb)+cosg*fac4)
1860           ecosg=(fac3*(fac1*cosb+cosg)+cosb*fac4)
1861           do k=1,3
1862             dcosb(k)=rmij*(dc_norm(k,i)-erij(k)*cosb)
1863             dcosg(k)=rmij*(dc_norm(k,j)-erij(k)*cosg)
1864           enddo
1865 cd        print '(2i3,2(3(1pd14.5),3x))',i,j,(dcosb(k),k=1,3),
1866 cd   &          (dcosg(k),k=1,3)
1867           do k=1,3
1868             ggg(k)=(ecosb*dcosb(k)+ecosg*dcosg(k))*
1869      &      fac_shield(i)**2*fac_shield(j)**2
1870
1871           enddo
1872 c          do k=1,3
1873 c            ghalf=0.5D0*ggg(k)
1874 c            gelc(k,i)=gelc(k,i)+ghalf
1875 c     &               +(ecosa*(dc_norm(k,j)-cosa*dc_norm(k,i))
1876 c     &               + ecosb*(erij(k)-cosb*dc_norm(k,i)))*vbld_inv(i+1)
1877 c            gelc(k,j)=gelc(k,j)+ghalf
1878 c     &               +(ecosa*(dc_norm(k,i)-cosa*dc_norm(k,j))
1879 c     &               + ecosg*(erij(k)-cosg*dc_norm(k,j)))*vbld_inv(j+1)
1880 c          enddo
1881 cgrad          do k=i+1,j-1
1882 cgrad            do l=1,3
1883 cgrad              gelc(l,k)=gelc(l,k)+ggg(l)
1884 cgrad            enddo
1885 cgrad          enddo
1886           do k=1,3
1887             gelc(k,i)=gelc(k,i)
1888      &               +((ecosa*(dc_norm(k,j)-cosa*dc_norm(k,i))
1889      &               + ecosb*(erij(k)-cosb*dc_norm(k,i)))*vbld_inv(i+1))
1890      &           *fac_shield(i)**2*fac_shield(j)**2
1891
1892             gelc(k,j)=gelc(k,j)
1893      &               +((ecosa*(dc_norm(k,i)-cosa*dc_norm(k,j))
1894      &               + ecosg*(erij(k)-cosg*dc_norm(k,j)))*vbld_inv(j+1))
1895      &           *fac_shield(i)**2*fac_shield(j)**2
1896             gelc_long(k,j)=gelc_long(k,j)+ggg(k)
1897             gelc_long(k,i)=gelc_long(k,i)-ggg(k)
1898           enddo
1899           IF (wel_loc.gt.0.0d0 .or. wcorr4.gt.0.0d0 .or. wcorr5.gt.0.0d0
1900      &        .or. wcorr6.gt.0.0d0 .or. wturn3.gt.0.0d0
1901      &        .or. wturn4.gt.0.0d0 .or. wturn6.gt.0.0d0) THEN
1902 C
1903 C 9/25/99 Mixed third-order local-electrostatic terms. The local-interaction
1904 C   energy of a peptide unit is assumed in the form of a second-order
1905 C   Fourier series in the angles lambda1 and lambda2 (see Nishikawa et al.
1906 C   Macromolecules, 1974, 7, 797-806 for definition). This correlation terms
1907 C   are computed for EVERY pair of non-contiguous peptide groups.
1908 C
1909           if (j.lt.nres-1) then
1910             j1=j+1
1911             j2=j-1
1912           else
1913             j1=j-1
1914             j2=j-2
1915           endif
1916           kkk=0
1917           do k=1,2
1918             do l=1,2
1919               kkk=kkk+1
1920               muij(kkk)=mu(k,i)*mu(l,j)
1921 #ifdef NEWCORR
1922              gmuij1(kkk)=gtb1(k,i+1)*mu(l,j)
1923 c             write(iout,*) 'k=',k,i,gtb1(k,i+1),gtb1(k,i+1)*mu(l,j)
1924              gmuij2(kkk)=gUb2(k,i)*mu(l,j)
1925              gmuji1(kkk)=mu(k,i)*gtb1(l,j+1)
1926 c             write(iout,*) 'l=',l,j,gtb1(l,j+1),gtb1(l,j+1)*mu(k,i)
1927              gmuji2(kkk)=mu(k,i)*gUb2(l,j)
1928 #endif
1929             enddo
1930           enddo
1931 cd         write (iout,*) 'EELEC: i',i,' j',j
1932 cd          write (iout,*) 'j',j,' j1',j1,' j2',j2
1933 cd          write(iout,*) 'muij',muij
1934           ury=scalar(uy(1,i),erij)
1935           urz=scalar(uz(1,i),erij)
1936           vry=scalar(uy(1,j),erij)
1937           vrz=scalar(uz(1,j),erij)
1938           a22=scalar(uy(1,i),uy(1,j))-3*ury*vry
1939           a23=scalar(uy(1,i),uz(1,j))-3*ury*vrz
1940           a32=scalar(uz(1,i),uy(1,j))-3*urz*vry
1941           a33=scalar(uz(1,i),uz(1,j))-3*urz*vrz
1942           fac=dsqrt(-ael6i)*r3ij
1943           a22=a22*fac
1944           a23=a23*fac
1945           a32=a32*fac
1946           a33=a33*fac
1947 cd          write (iout,'(4i5,4f10.5)')
1948 cd     &     i,itortyp(itype(i)),j,itortyp(itype(j)),a22,a23,a32,a33
1949 cd          write (iout,'(6f10.5)') (muij(k),k=1,4),fac,eel_loc_ij
1950 cd          write (iout,'(2(3f10.5,5x)/2(3f10.5,5x))') uy(:,i),uz(:,i),
1951 cd     &      uy(:,j),uz(:,j)
1952 cd          write (iout,'(4f10.5)')
1953 cd     &      scalar(uy(1,i),uy(1,j)),scalar(uy(1,i),uz(1,j)),
1954 cd     &      scalar(uz(1,i),uy(1,j)),scalar(uz(1,i),uz(1,j))
1955 cd          write (iout,'(4f10.5)') ury,urz,vry,vrz
1956 cd           write (iout,'(9f10.5/)')
1957 cd     &      fac22,a22,fac23,a23,fac32,a32,fac33,a33,eel_loc_ij
1958 C Derivatives of the elements of A in virtual-bond vectors
1959           call unormderiv(erij(1),unmat(1,1),rmij,erder(1,1))
1960           do k=1,3
1961             uryg(k,1)=scalar(erder(1,k),uy(1,i))
1962             uryg(k,2)=scalar(uygrad(1,k,1,i),erij(1))
1963             uryg(k,3)=scalar(uygrad(1,k,2,i),erij(1))
1964             urzg(k,1)=scalar(erder(1,k),uz(1,i))
1965             urzg(k,2)=scalar(uzgrad(1,k,1,i),erij(1))
1966             urzg(k,3)=scalar(uzgrad(1,k,2,i),erij(1))
1967             vryg(k,1)=scalar(erder(1,k),uy(1,j))
1968             vryg(k,2)=scalar(uygrad(1,k,1,j),erij(1))
1969             vryg(k,3)=scalar(uygrad(1,k,2,j),erij(1))
1970             vrzg(k,1)=scalar(erder(1,k),uz(1,j))
1971             vrzg(k,2)=scalar(uzgrad(1,k,1,j),erij(1))
1972             vrzg(k,3)=scalar(uzgrad(1,k,2,j),erij(1))
1973           enddo
1974 C Compute radial contributions to the gradient
1975           facr=-3.0d0*rrmij
1976           a22der=a22*facr
1977           a23der=a23*facr
1978           a32der=a32*facr
1979           a33der=a33*facr
1980           agg(1,1)=a22der*xj
1981           agg(2,1)=a22der*yj
1982           agg(3,1)=a22der*zj
1983           agg(1,2)=a23der*xj
1984           agg(2,2)=a23der*yj
1985           agg(3,2)=a23der*zj
1986           agg(1,3)=a32der*xj
1987           agg(2,3)=a32der*yj
1988           agg(3,3)=a32der*zj
1989           agg(1,4)=a33der*xj
1990           agg(2,4)=a33der*yj
1991           agg(3,4)=a33der*zj
1992 C Add the contributions coming from er
1993           fac3=-3.0d0*fac
1994           do k=1,3
1995             agg(k,1)=agg(k,1)+fac3*(uryg(k,1)*vry+vryg(k,1)*ury)
1996             agg(k,2)=agg(k,2)+fac3*(uryg(k,1)*vrz+vrzg(k,1)*ury)
1997             agg(k,3)=agg(k,3)+fac3*(urzg(k,1)*vry+vryg(k,1)*urz)
1998             agg(k,4)=agg(k,4)+fac3*(urzg(k,1)*vrz+vrzg(k,1)*urz)
1999           enddo
2000           do k=1,3
2001 C Derivatives in DC(i)
2002 cgrad            ghalf1=0.5d0*agg(k,1)
2003 cgrad            ghalf2=0.5d0*agg(k,2)
2004 cgrad            ghalf3=0.5d0*agg(k,3)
2005 cgrad            ghalf4=0.5d0*agg(k,4)
2006             aggi(k,1)=fac*(scalar(uygrad(1,k,1,i),uy(1,j))
2007      &      -3.0d0*uryg(k,2)*vry)!+ghalf1
2008             aggi(k,2)=fac*(scalar(uygrad(1,k,1,i),uz(1,j))
2009      &      -3.0d0*uryg(k,2)*vrz)!+ghalf2
2010             aggi(k,3)=fac*(scalar(uzgrad(1,k,1,i),uy(1,j))
2011      &      -3.0d0*urzg(k,2)*vry)!+ghalf3
2012             aggi(k,4)=fac*(scalar(uzgrad(1,k,1,i),uz(1,j))
2013      &      -3.0d0*urzg(k,2)*vrz)!+ghalf4
2014 C Derivatives in DC(i+1)
2015             aggi1(k,1)=fac*(scalar(uygrad(1,k,2,i),uy(1,j))
2016      &      -3.0d0*uryg(k,3)*vry)!+agg(k,1)
2017             aggi1(k,2)=fac*(scalar(uygrad(1,k,2,i),uz(1,j))
2018      &      -3.0d0*uryg(k,3)*vrz)!+agg(k,2)
2019             aggi1(k,3)=fac*(scalar(uzgrad(1,k,2,i),uy(1,j))
2020      &      -3.0d0*urzg(k,3)*vry)!+agg(k,3)
2021             aggi1(k,4)=fac*(scalar(uzgrad(1,k,2,i),uz(1,j))
2022      &      -3.0d0*urzg(k,3)*vrz)!+agg(k,4)
2023 C Derivatives in DC(j)
2024             aggj(k,1)=fac*(scalar(uygrad(1,k,1,j),uy(1,i))
2025      &      -3.0d0*vryg(k,2)*ury)!+ghalf1
2026             aggj(k,2)=fac*(scalar(uzgrad(1,k,1,j),uy(1,i))
2027      &      -3.0d0*vrzg(k,2)*ury)!+ghalf2
2028             aggj(k,3)=fac*(scalar(uygrad(1,k,1,j),uz(1,i))
2029      &      -3.0d0*vryg(k,2)*urz)!+ghalf3
2030             aggj(k,4)=fac*(scalar(uzgrad(1,k,1,j),uz(1,i))
2031      &      -3.0d0*vrzg(k,2)*urz)!+ghalf4
2032 C Derivatives in DC(j+1) or DC(nres-1)
2033             aggj1(k,1)=fac*(scalar(uygrad(1,k,2,j),uy(1,i))
2034      &      -3.0d0*vryg(k,3)*ury)
2035             aggj1(k,2)=fac*(scalar(uzgrad(1,k,2,j),uy(1,i))
2036      &      -3.0d0*vrzg(k,3)*ury)
2037             aggj1(k,3)=fac*(scalar(uygrad(1,k,2,j),uz(1,i))
2038      &      -3.0d0*vryg(k,3)*urz)
2039             aggj1(k,4)=fac*(scalar(uzgrad(1,k,2,j),uz(1,i))
2040      &      -3.0d0*vrzg(k,3)*urz)
2041 cgrad            if (j.eq.nres-1 .and. i.lt.j-2) then
2042 cgrad              do l=1,4
2043 cgrad                aggj1(k,l)=aggj1(k,l)+agg(k,l)
2044 cgrad              enddo
2045 cgrad            endif
2046           enddo
2047           acipa(1,1)=a22
2048           acipa(1,2)=a23
2049           acipa(2,1)=a32
2050           acipa(2,2)=a33
2051           a22=-a22
2052           a23=-a23
2053           do l=1,2
2054             do k=1,3
2055               agg(k,l)=-agg(k,l)
2056               aggi(k,l)=-aggi(k,l)
2057               aggi1(k,l)=-aggi1(k,l)
2058               aggj(k,l)=-aggj(k,l)
2059               aggj1(k,l)=-aggj1(k,l)
2060             enddo
2061           enddo
2062           if (j.lt.nres-1) then
2063             a22=-a22
2064             a32=-a32
2065             do l=1,3,2
2066               do k=1,3
2067                 agg(k,l)=-agg(k,l)
2068                 aggi(k,l)=-aggi(k,l)
2069                 aggi1(k,l)=-aggi1(k,l)
2070                 aggj(k,l)=-aggj(k,l)
2071                 aggj1(k,l)=-aggj1(k,l)
2072               enddo
2073             enddo
2074           else
2075             a22=-a22
2076             a23=-a23
2077             a32=-a32
2078             a33=-a33
2079             do l=1,4
2080               do k=1,3
2081                 agg(k,l)=-agg(k,l)
2082                 aggi(k,l)=-aggi(k,l)
2083                 aggi1(k,l)=-aggi1(k,l)
2084                 aggj(k,l)=-aggj(k,l)
2085                 aggj1(k,l)=-aggj1(k,l)
2086               enddo
2087             enddo
2088           endif
2089           ENDIF ! WCORR
2090           IF (wel_loc.gt.0.0d0) THEN
2091 C Contribution to the local-electrostatic energy coming from the i-j pair
2092           eel_loc_ij=a22*muij(1)+a23*muij(2)+a32*muij(3)
2093      &     +a33*muij(4)
2094 cd          write (iout,*) 'i',i,' j',j,' eel_loc_ij',eel_loc_ij
2095
2096           if (energy_dec) write (iout,'(a6,2i5,0pf7.3)')
2097      &            'eelloc',i,j,eel_loc_ij
2098
2099
2100           if (shield_mode.eq.0) then
2101            fac_shield(i)=1.0
2102            fac_shield(j)=1.0
2103 C          else
2104 C           fac_shield(i)=0.4
2105 C           fac_shield(j)=0.6
2106           endif
2107           eel_loc_ij=eel_loc_ij
2108      &    *fac_shield(i)*fac_shield(j)
2109           eel_loc=eel_loc+eel_loc_ij
2110
2111           if ((fac_shield(i).gt.0).and.(fac_shield(j).gt.0).and.
2112      &  (shield_mode.gt.0)) then
2113 C          print *,i,j
2114
2115           do ilist=1,ishield_list(i)
2116            iresshield=shield_list(ilist,i)
2117            do k=1,3
2118            rlocshield=grad_shield_side(k,ilist,i)*eel_loc_ij
2119      &                                          /fac_shield(i)
2120 C     &      *2.0
2121            gshieldx_ll(k,iresshield)=gshieldx_ll(k,iresshield)+
2122      &              rlocshield
2123      & +grad_shield_loc(k,ilist,i)*eel_loc_ij/fac_shield(i)
2124             gshieldc_ll(k,iresshield-1)=gshieldc_ll(k,iresshield-1)
2125      &      +rlocshield
2126            enddo
2127           enddo
2128           do ilist=1,ishield_list(j)
2129            iresshield=shield_list(ilist,j)
2130            do k=1,3
2131            rlocshield=grad_shield_side(k,ilist,j)*eel_loc_ij
2132      &                                       /fac_shield(j)
2133 C     &     *2.0
2134            gshieldx_ll(k,iresshield)=gshieldx_ll(k,iresshield)+
2135      &              rlocshield
2136      & +grad_shield_loc(k,ilist,j)*eel_loc_ij/fac_shield(j)
2137            gshieldc_ll(k,iresshield-1)=gshieldc_ll(k,iresshield-1)
2138      &             +rlocshield
2139
2140            enddo
2141           enddo
2142
2143           do k=1,3
2144             gshieldc_ll(k,i)=gshieldc_ll(k,i)+
2145      &              grad_shield(k,i)*eel_loc_ij/fac_shield(i)
2146             gshieldc_ll(k,j)=gshieldc_ll(k,j)+
2147      &              grad_shield(k,j)*eel_loc_ij/fac_shield(j)
2148             gshieldc_ll(k,i-1)=gshieldc_ll(k,i-1)+
2149      &              grad_shield(k,i)*eel_loc_ij/fac_shield(i)
2150             gshieldc_ll(k,j-1)=gshieldc_ll(k,j-1)+
2151      &              grad_shield(k,j)*eel_loc_ij/fac_shield(j)
2152            enddo
2153            endif
2154
2155 #ifdef NEWCORR
2156          geel_loc_ij=(a22*gmuij1(1)
2157      &     +a23*gmuij1(2)
2158      &     +a32*gmuij1(3)
2159      &     +a33*gmuij1(4))
2160      &    *fac_shield(i)*fac_shield(j)
2161 c         write(iout,*) "derivative over thatai"
2162 c         write(iout,*) a22*gmuij1(1), a23*gmuij1(2) ,a32*gmuij1(3),
2163 c     &   a33*gmuij1(4)
2164          gloc(nphi+i,icg)=gloc(nphi+i,icg)+
2165      &      geel_loc_ij*wel_loc
2166 c         write(iout,*) "derivative over thatai-1"
2167 c         write(iout,*) a22*gmuij2(1), a23*gmuij2(2) ,a32*gmuij2(3),
2168 c     &   a33*gmuij2(4)
2169          geel_loc_ij=
2170      &     a22*gmuij2(1)
2171      &     +a23*gmuij2(2)
2172      &     +a32*gmuij2(3)
2173      &     +a33*gmuij2(4)
2174          gloc(nphi+i-1,icg)=gloc(nphi+i-1,icg)+
2175      &      geel_loc_ij*wel_loc
2176      &    *fac_shield(i)*fac_shield(j)
2177
2178 c  Derivative over j residue
2179          geel_loc_ji=a22*gmuji1(1)
2180      &     +a23*gmuji1(2)
2181      &     +a32*gmuji1(3)
2182      &     +a33*gmuji1(4)
2183 c         write(iout,*) "derivative over thataj"
2184 c         write(iout,*) a22*gmuji1(1), a23*gmuji1(2) ,a32*gmuji1(3),
2185 c     &   a33*gmuji1(4)
2186
2187         gloc(nphi+j,icg)=gloc(nphi+j,icg)+
2188      &      geel_loc_ji*wel_loc
2189      &    *fac_shield(i)*fac_shield(j)
2190
2191          geel_loc_ji=
2192      &     +a22*gmuji2(1)
2193      &     +a23*gmuji2(2)
2194      &     +a32*gmuji2(3)
2195      &     +a33*gmuji2(4)
2196 c         write(iout,*) "derivative over thataj-1"
2197 c         write(iout,*) a22*gmuji2(1), a23*gmuji2(2) ,a32*gmuji2(3),
2198 c     &   a33*gmuji2(4)
2199          gloc(nphi+j-1,icg)=gloc(nphi+j-1,icg)+
2200      &      geel_loc_ji*wel_loc
2201      &    *fac_shield(i)*fac_shield(j)
2202 #endif
2203 cC Partial derivatives in virtual-bond dihedral angles gamma
2204           if (i.gt.1)
2205      &    gel_loc_loc(i-1)=gel_loc_loc(i-1)+
2206      &            (a22*muder(1,i)*mu(1,j)+a23*muder(1,i)*mu(2,j)
2207      &           +a32*muder(2,i)*mu(1,j)+a33*muder(2,i)*mu(2,j))
2208      &    *fac_shield(i)*fac_shield(j)
2209
2210           gel_loc_loc(j-1)=gel_loc_loc(j-1)+
2211      &            (a22*mu(1,i)*muder(1,j)+a23*mu(1,i)*muder(2,j)
2212      &           +a32*mu(2,i)*muder(1,j)+a33*mu(2,i)*muder(2,j))
2213      &    *fac_shield(i)*fac_shield(j)
2214
2215 C Derivatives of eello in DC(i+1) thru DC(j-1) or DC(nres-2)
2216           do l=1,3
2217             ggg(l)=(agg(l,1)*muij(1)+
2218      &          agg(l,2)*muij(2)+agg(l,3)*muij(3)+agg(l,4)*muij(4))
2219      &    *fac_shield(i)*fac_shield(j)
2220
2221             gel_loc_long(l,j)=gel_loc_long(l,j)+ggg(l)
2222             gel_loc_long(l,i)=gel_loc_long(l,i)-ggg(l)
2223 cgrad            ghalf=0.5d0*ggg(l)
2224 cgrad            gel_loc(l,i)=gel_loc(l,i)+ghalf
2225 cgrad            gel_loc(l,j)=gel_loc(l,j)+ghalf
2226           enddo
2227 cgrad          do k=i+1,j2
2228 cgrad            do l=1,3
2229 cgrad              gel_loc(l,k)=gel_loc(l,k)+ggg(l)
2230 cgrad            enddo
2231 cgrad          enddo
2232 C Remaining derivatives of eello
2233           do l=1,3
2234             gel_loc(l,i)=gel_loc(l,i)+(aggi(l,1)*muij(1)+
2235      &          aggi(l,2)*muij(2)+aggi(l,3)*muij(3)+aggi(l,4)*muij(4))
2236      &    *fac_shield(i)*fac_shield(j)
2237
2238             gel_loc(l,i+1)=gel_loc(l,i+1)+(aggi1(l,1)*muij(1)+
2239      &         aggi1(l,2)*muij(2)+aggi1(l,3)*muij(3)+aggi1(l,4)*muij(4))
2240      &    *fac_shield(i)*fac_shield(j)
2241
2242             gel_loc(l,j)=gel_loc(l,j)+(aggj(l,1)*muij(1)+
2243      &          aggj(l,2)*muij(2)+aggj(l,3)*muij(3)+aggj(l,4)*muij(4))
2244      &    *fac_shield(i)*fac_shield(j)
2245
2246             gel_loc(l,j1)=gel_loc(l,j1)+(aggj1(l,1)*muij(1)+
2247      &         aggj1(l,2)*muij(2)+aggj1(l,3)*muij(3)+aggj1(l,4)*muij(4))
2248      &    *fac_shield(i)*fac_shield(j)
2249
2250           enddo
2251           ENDIF
2252 C Change 12/26/95 to calculate four-body contributions to H-bonding energy
2253 c          if (j.gt.i+1 .and. num_conti.le.maxconts) then
2254           if (wcorr+wcorr4+wcorr5+wcorr6.gt.0.0d0
2255      &       .and. num_conti.le.maxconts) then
2256 c            write (iout,*) i,j," entered corr"
2257 C
2258 C Calculate the contact function. The ith column of the array JCONT will
2259 C contain the numbers of atoms that make contacts with the atom I (of numbers
2260 C greater than I). The arrays FACONT and GACONT will contain the values of
2261 C the contact function and its derivative.
2262 c           r0ij=1.02D0*rpp(iteli,itelj)
2263 c           r0ij=1.11D0*rpp(iteli,itelj)
2264             r0ij=2.20D0*rpp(iteli,itelj)
2265 c           r0ij=1.55D0*rpp(iteli,itelj)
2266             call gcont(rij,r0ij,1.0D0,0.2d0*r0ij,fcont,fprimcont)
2267             if (fcont.gt.0.0D0) then
2268               num_conti=num_conti+1
2269               if (num_conti.gt.maxconts) then
2270                 write (iout,*) 'WARNING - max. # of contacts exceeded;',
2271      &                         ' will skip next contacts for this conf.'
2272               else
2273                 jcont_hb(num_conti,i)=j
2274 cd                write (iout,*) "i",i," j",j," num_conti",num_conti,
2275 cd     &           " jcont_hb",jcont_hb(num_conti,i)
2276                 IF (wcorr4.gt.0.0d0 .or. wcorr5.gt.0.0d0 .or.
2277      &          wcorr6.gt.0.0d0 .or. wturn6.gt.0.0d0) THEN
2278 C 9/30/99 (AL) - store components necessary to evaluate higher-order loc-el
2279 C  terms.
2280                 d_cont(num_conti,i)=rij
2281 cd                write (2,'(3e15.5)') rij,r0ij+0.2d0*r0ij,rij
2282 C     --- Electrostatic-interaction matrix ---
2283                 a_chuj(1,1,num_conti,i)=a22
2284                 a_chuj(1,2,num_conti,i)=a23
2285                 a_chuj(2,1,num_conti,i)=a32
2286                 a_chuj(2,2,num_conti,i)=a33
2287 C     --- Gradient of rij
2288                 do kkk=1,3
2289                   grij_hb_cont(kkk,num_conti,i)=erij(kkk)
2290                 enddo
2291                 kkll=0
2292                 do k=1,2
2293                   do l=1,2
2294                     kkll=kkll+1
2295                     do m=1,3
2296                       a_chuj_der(k,l,m,1,num_conti,i)=agg(m,kkll)
2297                       a_chuj_der(k,l,m,2,num_conti,i)=aggi(m,kkll)
2298                       a_chuj_der(k,l,m,3,num_conti,i)=aggi1(m,kkll)
2299                       a_chuj_der(k,l,m,4,num_conti,i)=aggj(m,kkll)
2300                       a_chuj_der(k,l,m,5,num_conti,i)=aggj1(m,kkll)
2301                     enddo
2302                   enddo
2303                 enddo
2304                 ENDIF
2305                 IF (wcorr4.eq.0.0d0 .and. wcorr.gt.0.0d0) THEN
2306 C Calculate contact energies
2307                 cosa4=4.0D0*cosa
2308                 wij=cosa-3.0D0*cosb*cosg
2309                 cosbg1=cosb+cosg
2310                 cosbg2=cosb-cosg
2311 c               fac3=dsqrt(-ael6i)/r0ij**3
2312                 fac3=dsqrt(-ael6i)*r3ij
2313 c                 ees0pij=dsqrt(4.0D0+cosa4+wij*wij-3.0D0*cosbg1*cosbg1)
2314                 ees0tmp=4.0D0+cosa4+wij*wij-3.0D0*cosbg1*cosbg1
2315                 if (ees0tmp.gt.0) then
2316                   ees0pij=dsqrt(ees0tmp)
2317                 else
2318                   ees0pij=0
2319                 endif
2320 c                ees0mij=dsqrt(4.0D0-cosa4+wij*wij-3.0D0*cosbg2*cosbg2)
2321                 ees0tmp=4.0D0-cosa4+wij*wij-3.0D0*cosbg2*cosbg2
2322                 if (ees0tmp.gt.0) then
2323                   ees0mij=dsqrt(ees0tmp)
2324                 else
2325                   ees0mij=0
2326                 endif
2327 c               ees0mij=0.0D0
2328                 if (shield_mode.eq.0) then
2329                 fac_shield(i)=1.0d0
2330                 fac_shield(j)=1.0d0
2331                 else
2332                 ees0plist(num_conti,i)=j
2333 C                fac_shield(i)=0.4d0
2334 C                fac_shield(j)=0.6d0
2335                 endif
2336                 ees0p(num_conti,i)=0.5D0*fac3*(ees0pij+ees0mij)
2337      &          *fac_shield(i)*fac_shield(j)
2338                 ees0m(num_conti,i)=0.5D0*fac3*(ees0pij-ees0mij)
2339      &          *fac_shield(i)*fac_shield(j)
2340
2341 C Diagnostics. Comment out or remove after debugging!
2342 c               ees0p(num_conti,i)=0.5D0*fac3*ees0pij
2343 c               ees0m(num_conti,i)=0.5D0*fac3*ees0mij
2344 c               ees0m(num_conti,i)=0.0D0
2345 C End diagnostics.
2346 c               write (iout,*) 'i=',i,' j=',j,' rij=',rij,' r0ij=',r0ij,
2347 c    & ' ees0ij=',ees0p(num_conti,i),ees0m(num_conti,i),' fcont=',fcont
2348 C Angular derivatives of the contact function
2349                 ees0pij1=fac3/ees0pij
2350                 ees0mij1=fac3/ees0mij
2351                 fac3p=-3.0D0*fac3*rrmij
2352                 ees0pijp=0.5D0*fac3p*(ees0pij+ees0mij)
2353                 ees0mijp=0.5D0*fac3p*(ees0pij-ees0mij)
2354 c               ees0mij1=0.0D0
2355                 ecosa1=       ees0pij1*( 1.0D0+0.5D0*wij)
2356                 ecosb1=-1.5D0*ees0pij1*(wij*cosg+cosbg1)
2357                 ecosg1=-1.5D0*ees0pij1*(wij*cosb+cosbg1)
2358                 ecosa2=       ees0mij1*(-1.0D0+0.5D0*wij)
2359                 ecosb2=-1.5D0*ees0mij1*(wij*cosg+cosbg2)
2360                 ecosg2=-1.5D0*ees0mij1*(wij*cosb-cosbg2)
2361                 ecosap=ecosa1+ecosa2
2362                 ecosbp=ecosb1+ecosb2
2363                 ecosgp=ecosg1+ecosg2
2364                 ecosam=ecosa1-ecosa2
2365                 ecosbm=ecosb1-ecosb2
2366                 ecosgm=ecosg1-ecosg2
2367 C Diagnostics
2368 c               ecosap=ecosa1
2369 c               ecosbp=ecosb1
2370 c               ecosgp=ecosg1
2371 c               ecosam=0.0D0
2372 c               ecosbm=0.0D0
2373 c               ecosgm=0.0D0
2374 C End diagnostics
2375                 facont_hb(num_conti,i)=fcont
2376                 fprimcont=fprimcont/rij
2377 cd              facont_hb(num_conti,i)=1.0D0
2378 C Following line is for diagnostics.
2379 cd              fprimcont=0.0D0
2380                 do k=1,3
2381                   dcosb(k)=rmij*(dc_norm(k,i)-erij(k)*cosb)
2382                   dcosg(k)=rmij*(dc_norm(k,j)-erij(k)*cosg)
2383                 enddo
2384                 do k=1,3
2385                   gggp(k)=ecosbp*dcosb(k)+ecosgp*dcosg(k)
2386                   gggm(k)=ecosbm*dcosb(k)+ecosgm*dcosg(k)
2387                 enddo
2388                 gggp(1)=gggp(1)+ees0pijp*xj
2389                 gggp(2)=gggp(2)+ees0pijp*yj
2390                 gggp(3)=gggp(3)+ees0pijp*zj
2391                 gggm(1)=gggm(1)+ees0mijp*xj
2392                 gggm(2)=gggm(2)+ees0mijp*yj
2393                 gggm(3)=gggm(3)+ees0mijp*zj
2394 C Derivatives due to the contact function
2395                 gacont_hbr(1,num_conti,i)=fprimcont*xj
2396                 gacont_hbr(2,num_conti,i)=fprimcont*yj
2397                 gacont_hbr(3,num_conti,i)=fprimcont*zj
2398                 do k=1,3
2399 c
2400 c 10/24/08 cgrad and ! comments indicate the parts of the code removed
2401 c          following the change of gradient-summation algorithm.
2402 c
2403 cgrad                  ghalfp=0.5D0*gggp(k)
2404 cgrad                  ghalfm=0.5D0*gggm(k)
2405                   gacontp_hb1(k,num_conti,i)=!ghalfp
2406      &              +(ecosap*(dc_norm(k,j)-cosa*dc_norm(k,i))
2407      &              + ecosbp*(erij(k)-cosb*dc_norm(k,i)))*vbld_inv(i+1)
2408      &          *fac_shield(i)*fac_shield(j)
2409
2410                   gacontp_hb2(k,num_conti,i)=!ghalfp
2411      &              +(ecosap*(dc_norm(k,i)-cosa*dc_norm(k,j))
2412      &              + ecosgp*(erij(k)-cosg*dc_norm(k,j)))*vbld_inv(j+1)
2413      &          *fac_shield(i)*fac_shield(j)
2414
2415                   gacontp_hb3(k,num_conti,i)=gggp(k)
2416      &          *fac_shield(i)*fac_shield(j)
2417
2418                   gacontm_hb1(k,num_conti,i)=!ghalfm
2419      &              +(ecosam*(dc_norm(k,j)-cosa*dc_norm(k,i))
2420      &              + ecosbm*(erij(k)-cosb*dc_norm(k,i)))*vbld_inv(i+1)
2421      &          *fac_shield(i)*fac_shield(j)
2422
2423                   gacontm_hb2(k,num_conti,i)=!ghalfm
2424      &              +(ecosam*(dc_norm(k,i)-cosa*dc_norm(k,j))
2425      &              + ecosgm*(erij(k)-cosg*dc_norm(k,j)))*vbld_inv(j+1)
2426      &          *fac_shield(i)*fac_shield(j)
2427
2428                   gacontm_hb3(k,num_conti,i)=gggm(k)
2429      &          *fac_shield(i)*fac_shield(j)
2430
2431                 enddo
2432               ENDIF ! wcorr
2433               endif  ! num_conti.le.maxconts
2434             endif  ! fcont.gt.0
2435           endif    ! j.gt.i+1
2436           if (wturn3.gt.0.0d0 .or. wturn4.gt.0.0d0) then
2437             do k=1,4
2438               do l=1,3
2439                 ghalf=0.5d0*agg(l,k)
2440                 aggi(l,k)=aggi(l,k)+ghalf
2441                 aggi1(l,k)=aggi1(l,k)+agg(l,k)
2442                 aggj(l,k)=aggj(l,k)+ghalf
2443               enddo
2444             enddo
2445             if (j.eq.nres-1 .and. i.lt.j-2) then
2446               do k=1,4
2447                 do l=1,3
2448                   aggj1(l,k)=aggj1(l,k)+agg(l,k)
2449                 enddo
2450               enddo
2451             endif
2452           endif
2453 c          t_eelecij=t_eelecij+MPI_Wtime()-time00
2454       return
2455       end
2456 C-----------------------------------------------------------------------
2457       subroutine evdwpp_short(evdw1)
2458 C
2459 C Compute Evdwpp
2460 C
2461       implicit real*8 (a-h,o-z)
2462       include 'DIMENSIONS'
2463       include 'COMMON.CONTROL'
2464       include 'COMMON.IOUNITS'
2465       include 'COMMON.GEO'
2466       include 'COMMON.VAR'
2467       include 'COMMON.LOCAL'
2468       include 'COMMON.CHAIN'
2469       include 'COMMON.DERIV'
2470       include 'COMMON.INTERACT'
2471       include 'COMMON.CONTACTS'
2472       include 'COMMON.TORSION'
2473       include 'COMMON.VECTORS'
2474       include 'COMMON.FFIELD'
2475       dimension ggg(3)
2476       integer xshift,yshift,zshift
2477 c 4/26/02 - AL scaling factor for 1,4 repulsive VDW interactions
2478 #ifdef MOMENT
2479       double precision scal_el /1.0d0/
2480 #else
2481       double precision scal_el /0.5d0/
2482 #endif
2483 c      write (iout,*) "evdwpp_short"
2484       evdw1=0.0D0
2485 C      print *,"WCHODZE"
2486 c      write (iout,*) "iatel_s_vdw",iatel_s_vdw,
2487 c     & " iatel_e_vdw",iatel_e_vdw
2488 c      call flush(iout)
2489       do i=iatel_s_vdw,iatel_e_vdw
2490         if (itype(i).eq.ntyp1.or. itype(i+1).eq.ntyp1) cycle
2491         dxi=dc(1,i)
2492         dyi=dc(2,i)
2493         dzi=dc(3,i)
2494         dx_normi=dc_norm(1,i)
2495         dy_normi=dc_norm(2,i)
2496         dz_normi=dc_norm(3,i)
2497         xmedi=c(1,i)+0.5d0*dxi
2498         ymedi=c(2,i)+0.5d0*dyi
2499         zmedi=c(3,i)+0.5d0*dzi
2500           xmedi=mod(xmedi,boxxsize)
2501           if (xmedi.lt.0.0d0) xmedi=xmedi+boxxsize
2502           ymedi=mod(ymedi,boxysize)
2503           if (ymedi.lt.0.0d0) ymedi=ymedi+boxysize
2504           zmedi=mod(zmedi,boxzsize)
2505           if (zmedi.lt.0.0d0) zmedi=zmedi+boxzsize
2506         num_conti=0
2507 c        write (iout,*) 'i',i,' ielstart',ielstart_vdw(i),
2508 c     &   ' ielend',ielend_vdw(i)
2509 c        call flush(iout)
2510         do j=ielstart_vdw(i),ielend_vdw(i)
2511           if (itype(j).eq.ntyp1 .or. itype(j+1).eq.ntyp1) cycle
2512           ind=ind+1
2513           iteli=itel(i)
2514           itelj=itel(j)
2515           if (j.eq.i+2 .and. itelj.eq.2) iteli=2
2516           aaa=app(iteli,itelj)
2517           bbb=bpp(iteli,itelj)
2518           dxj=dc(1,j)
2519           dyj=dc(2,j)
2520           dzj=dc(3,j)
2521           dx_normj=dc_norm(1,j)
2522           dy_normj=dc_norm(2,j)
2523           dz_normj=dc_norm(3,j)
2524           xj=c(1,j)+0.5D0*dxj
2525           yj=c(2,j)+0.5D0*dyj
2526           zj=c(3,j)+0.5D0*dzj
2527           xj=mod(xj,boxxsize)
2528           if (xj.lt.0) xj=xj+boxxsize
2529           yj=mod(yj,boxysize)
2530           if (yj.lt.0) yj=yj+boxysize
2531           zj=mod(zj,boxzsize)
2532           if (zj.lt.0) zj=zj+boxzsize
2533       dist_init=(xj-xmedi)**2+(yj-ymedi)**2+(zj-zmedi)**2
2534       xj_safe=xj
2535       yj_safe=yj
2536       zj_safe=zj
2537       isubchap=0
2538       do xshift=-1,1
2539       do yshift=-1,1
2540       do zshift=-1,1
2541           xj=xj_safe+xshift*boxxsize
2542           yj=yj_safe+yshift*boxysize
2543           zj=zj_safe+zshift*boxzsize
2544           dist_temp=(xj-xmedi)**2+(yj-ymedi)**2+(zj-zmedi)**2
2545           if(dist_temp.lt.dist_init) then
2546             dist_init=dist_temp
2547             xj_temp=xj
2548             yj_temp=yj
2549             zj_temp=zj
2550             isubchap=1
2551           endif
2552        enddo
2553        enddo
2554        enddo
2555        if (isubchap.eq.1) then
2556           xj=xj_temp-xmedi
2557           yj=yj_temp-ymedi
2558           zj=zj_temp-zmedi
2559        else
2560           xj=xj_safe-xmedi
2561           yj=yj_safe-ymedi
2562           zj=zj_safe-zmedi
2563        endif
2564           rij=xj*xj+yj*yj+zj*zj
2565           rrmij=1.0D0/rij
2566           rij=dsqrt(rij)
2567 c            sss=sscale(rij/rpp(iteli,itelj))
2568 c            sssgrad=sscagrad(rij/rpp(iteli,itelj))
2569           sss=sscale(rij)
2570           sssgrad=sscagrad(rij)
2571           if (sss.gt.0.0d0) then
2572             rmij=1.0D0/rij
2573             r3ij=rrmij*rmij
2574             r6ij=r3ij*r3ij
2575             ev1=aaa*r6ij*r6ij
2576 c 4/26/02 - AL scaling down 1,4 repulsive VDW interactions
2577             if (j.eq.i+2) ev1=scal_el*ev1
2578             ev2=bbb*r6ij
2579             evdwij=ev1+ev2
2580             if (energy_dec) then
2581               write (iout,'(a6,2i5,0pf7.3,f7.3)') 'evdw1',i,j,evdwij,sss
2582             endif
2583             evdw1=evdw1+evdwij*sss
2584             if (energy_dec) write (iout,'(a10,2i5,0pf7.3)')
2585      &        'evdw1_sum',i,j,evdw1
2586 C
2587 C Calculate contributions to the Cartesian gradient.
2588 C
2589             facvdw=-6*rrmij*(ev1+evdwij)*sss
2590           ggg(1)=facvdw*xj+sssgrad*rmij*evdwij*xj/rpp(iteli,itelj)
2591           ggg(2)=facvdw*yj+sssgrad*rmij*evdwij*yj/rpp(iteli,itelj)
2592           ggg(3)=facvdw*zj+sssgrad*rmij*evdwij*zj/rpp(iteli,itelj)
2593 C            ggg(1)=facvdw*xj
2594 C            ggg(2)=facvdw*yj
2595 C            ggg(3)=facvdw*zj
2596             do k=1,3
2597               gvdwpp(k,j)=gvdwpp(k,j)+ggg(k)
2598               gvdwpp(k,i)=gvdwpp(k,i)-ggg(k)
2599             enddo
2600           endif
2601         enddo ! j
2602       enddo   ! i
2603       return
2604       end
2605 C-----------------------------------------------------------------------------
2606       subroutine escp_long(evdw2,evdw2_14)
2607 C
2608 C This subroutine calculates the excluded-volume interaction energy between
2609 C peptide-group centers and side chains and its gradient in virtual-bond and
2610 C side-chain vectors.
2611 C
2612       implicit real*8 (a-h,o-z)
2613       include 'DIMENSIONS'
2614       include 'COMMON.GEO'
2615       include 'COMMON.VAR'
2616       include 'COMMON.LOCAL'
2617       include 'COMMON.CHAIN'
2618       include 'COMMON.DERIV'
2619       include 'COMMON.INTERACT'
2620       include 'COMMON.FFIELD'
2621       include 'COMMON.IOUNITS'
2622       include 'COMMON.CONTROL'
2623       logical lprint_short
2624       common /shortcheck/ lprint_short
2625       dimension ggg(3)
2626       integer xshift,yshift,zshift
2627       if (energy_dec) write (iout,*) "escp_long:",r_cut,rlamb
2628       evdw2=0.0D0
2629       evdw2_14=0.0d0
2630 CD        print '(a)','Enter ESCP KURWA'
2631 cd    write (iout,*) 'iatscp_s=',iatscp_s,' iatscp_e=',iatscp_e
2632 c      if (lprint_short)
2633 c     &  write (iout,*) 'ESCP_LONG iatscp_s=',iatscp_s,
2634 c     & ' iatscp_e=',iatscp_e
2635       do i=iatscp_s,iatscp_e
2636         if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1) cycle
2637         iteli=itel(i)
2638         xi=0.5D0*(c(1,i)+c(1,i+1))
2639         yi=0.5D0*(c(2,i)+c(2,i+1))
2640         zi=0.5D0*(c(3,i)+c(3,i+1))
2641          xi=mod(xi,boxxsize)
2642           if (xi.lt.0) xi=xi+boxxsize
2643           yi=mod(yi,boxysize)
2644           if (yi.lt.0) yi=yi+boxysize
2645           zi=mod(zi,boxzsize)
2646           if (zi.lt.0) zi=zi+boxzsize
2647         do iint=1,nscp_gr(i)
2648
2649         do j=iscpstart(i,iint),iscpend(i,iint)
2650           itypj=iabs(itype(j))
2651           if (itypj.eq.ntyp1) cycle
2652 C Uncomment following three lines for SC-p interactions
2653 c         xj=c(1,nres+j)-xi
2654 c         yj=c(2,nres+j)-yi
2655 c         zj=c(3,nres+j)-zi
2656 C Uncomment following three lines for Ca-p interactions
2657           xj=c(1,j)
2658           yj=c(2,j)
2659           zj=c(3,j)
2660 c corrected by AL
2661           xj=mod(xj,boxxsize)
2662           if (xj.lt.0) xj=xj+boxxsize
2663           yj=mod(yj,boxysize)
2664           if (yj.lt.0) yj=yj+boxysize
2665           zj=mod(zj,boxzsize)
2666           if (zj.lt.0) zj=zj+boxzsize
2667 c end correction
2668       dist_init=(xj-xi)**2+(yj-yi)**2+(zj-zi)**2
2669       xj_safe=xj
2670       yj_safe=yj
2671       zj_safe=zj
2672       subchap=0
2673       do xshift=-1,1
2674       do yshift=-1,1
2675       do zshift=-1,1
2676           xj=xj_safe+xshift*boxxsize
2677           yj=yj_safe+yshift*boxysize
2678           zj=zj_safe+zshift*boxzsize
2679           dist_temp=(xj-xi)**2+(yj-yi)**2+(zj-zi)**2
2680           if(dist_temp.lt.dist_init) then
2681             dist_init=dist_temp
2682             xj_temp=xj
2683             yj_temp=yj
2684             zj_temp=zj
2685             subchap=1
2686           endif
2687        enddo
2688        enddo
2689        enddo
2690        if (subchap.eq.1) then
2691           xj=xj_temp-xi
2692           yj=yj_temp-yi
2693           zj=zj_temp-zi
2694        else
2695           xj=xj_safe-xi
2696           yj=yj_safe-yi
2697           zj=zj_safe-zi
2698        endif
2699
2700           rrij=1.0D0/(xj*xj+yj*yj+zj*zj)
2701
2702           sss=sscale(1.0d0/(dsqrt(rrij)*rscp(itypj,iteli)))
2703           sssgrad=sscagrad(1.0d0/(dsqrt(rrij)*rscp(itypj,iteli)))
2704           if (energy_dec) write (iout,*) "rrij",1.0d0/dsqrt(rrij),
2705      &     " rscp",rscp(itypj,iteli)," subchap",subchap," sss",sss
2706           if (sss.lt.1.0d0) then
2707             fac=rrij**expon2
2708             e1=fac*fac*aad(itypj,iteli)
2709             e2=fac*bad(itypj,iteli)
2710             if (iabs(j-i) .le. 2) then
2711               e1=scal14*e1
2712               e2=scal14*e2
2713               evdw2_14=evdw2_14+(e1+e2)*(1.0d0-sss)
2714             endif
2715             evdwij=e1+e2
2716             evdw2=evdw2+evdwij*(1.0d0-sss)
2717             if (energy_dec) write (iout,'(a6,2i5,2(0pf7.3))')
2718      &          'evdw2',i,j,sss,evdwij
2719 C
2720 C Calculate contributions to the gradient in the virtual-bond and SC vectors.
2721 C
2722
2723             fac=-(evdwij+e1)*rrij*(1.0d0-sss)
2724             fac=fac-(evdwij)*sssgrad*dsqrt(rrij)/rscp(itypj,iteli)
2725             ggg(1)=xj*fac
2726             ggg(2)=yj*fac
2727             ggg(3)=zj*fac
2728 C Uncomment following three lines for SC-p interactions
2729 c           do k=1,3
2730 c             gradx_scp(k,j)=gradx_scp(k,j)+ggg(k)
2731 c           enddo
2732 C Uncomment following line for SC-p interactions
2733 c             gradx_scp(k,j)=gradx_scp(k,j)+ggg(k)
2734             do k=1,3
2735               gvdwc_scpp(k,i)=gvdwc_scpp(k,i)-ggg(k)
2736               gvdwc_scp(k,j)=gvdwc_scp(k,j)+ggg(k)
2737             enddo
2738           endif
2739         enddo
2740
2741         enddo ! iint
2742       enddo ! i
2743       do i=1,nct
2744         do j=1,3
2745           gvdwc_scp(j,i)=expon*gvdwc_scp(j,i)
2746           gvdwc_scpp(j,i)=expon*gvdwc_scpp(j,i)
2747           gradx_scp(j,i)=expon*gradx_scp(j,i)
2748         enddo
2749       enddo
2750 C******************************************************************************
2751 C
2752 C                              N O T E !!!
2753 C
2754 C To save time the factor EXPON has been extracted from ALL components
2755 C of GVDWC and GRADX. Remember to multiply them by this factor before further
2756 C use!
2757 C
2758 C******************************************************************************
2759       return
2760       end
2761 C-----------------------------------------------------------------------------
2762       subroutine escp_short(evdw2,evdw2_14)
2763 C
2764 C This subroutine calculates the excluded-volume interaction energy between
2765 C peptide-group centers and side chains and its gradient in virtual-bond and
2766 C side-chain vectors.
2767 C
2768       implicit real*8 (a-h,o-z)
2769       include 'DIMENSIONS'
2770       include 'COMMON.GEO'
2771       include 'COMMON.VAR'
2772       include 'COMMON.LOCAL'
2773       include 'COMMON.CHAIN'
2774       include 'COMMON.DERIV'
2775       include 'COMMON.INTERACT'
2776       include 'COMMON.FFIELD'
2777       include 'COMMON.IOUNITS'
2778       include 'COMMON.CONTROL'
2779       integer xshift,yshift,zshift
2780       logical lprint_short
2781       common /shortcheck/ lprint_short
2782       dimension ggg(3)
2783       evdw2=0.0D0
2784       evdw2_14=0.0d0
2785 cd    print '(a)','Enter ESCP'
2786 c      if (lprint_short)
2787 c     &  write (iout,*) 'ESCP_SHORT iatscp_s=',iatscp_s,
2788 c     & ' iatscp_e=',iatscp_e
2789       if (energy_dec) write (iout,*) "escp_short:",r_cut,rlamb
2790       do i=iatscp_s,iatscp_e
2791         if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1) cycle
2792         iteli=itel(i)
2793         xi=0.5D0*(c(1,i)+c(1,i+1))
2794         yi=0.5D0*(c(2,i)+c(2,i+1))
2795         zi=0.5D0*(c(3,i)+c(3,i+1))
2796          xi=mod(xi,boxxsize)
2797           if (xi.lt.0) xi=xi+boxxsize
2798           yi=mod(yi,boxysize)
2799           if (yi.lt.0) yi=yi+boxysize
2800           zi=mod(zi,boxzsize)
2801           if (zi.lt.0) zi=zi+boxzsize
2802
2803 c        if (lprint_short)
2804 c     &    write (iout,*) "i",i," itype",itype(i),itype(i+1),
2805 c     &     " nscp_gr",nscp_gr(i)
2806         do iint=1,nscp_gr(i)
2807
2808         do j=iscpstart(i,iint),iscpend(i,iint)
2809           itypj=iabs(itype(j))
2810 c        if (lprint_short)
2811 c     &    write (iout,*) "j",j," itypj",itypj
2812           if (itypj.eq.ntyp1) cycle
2813 C Uncomment following three lines for SC-p interactions
2814 c         xj=c(1,nres+j)-xi
2815 c         yj=c(2,nres+j)-yi
2816 c         zj=c(3,nres+j)-zi
2817 C Uncomment following three lines for Ca-p interactions
2818           xj=c(1,j)
2819           yj=c(2,j)
2820           zj=c(3,j)
2821 c corrected by AL
2822           xj=mod(xj,boxxsize)
2823           if (xj.lt.0) xj=xj+boxxsize
2824           yj=mod(yj,boxysize)
2825           if (yj.lt.0) yj=yj+boxysize
2826           zj=mod(zj,boxzsize)
2827           if (zj.lt.0) zj=zj+boxzsize
2828 c end correction
2829       dist_init=(xj-xi)**2+(yj-yi)**2+(zj-zi)**2
2830 c          if (lprint_short) then
2831 c            write (iout,*) i,j,xi,yi,zi,xj,yj,zj
2832 c            write (iout,*) "dist_init",dsqrt(dist_init)
2833 c          endif
2834       xj_safe=xj
2835       yj_safe=yj
2836       zj_safe=zj
2837       subchap=0
2838       do xshift=-1,1
2839       do yshift=-1,1
2840       do zshift=-1,1
2841           xj=xj_safe+xshift*boxxsize
2842           yj=yj_safe+yshift*boxysize
2843           zj=zj_safe+zshift*boxzsize
2844           dist_temp=(xj-xi)**2+(yj-yi)**2+(zj-zi)**2
2845           if(dist_temp.lt.dist_init) then
2846             dist_init=dist_temp
2847             xj_temp=xj
2848             yj_temp=yj
2849             zj_temp=zj
2850             subchap=1
2851           endif
2852        enddo
2853        enddo
2854        enddo
2855 c          if (lprint_short) write (iout,*) "dist_temp",dsqrt(dist_temp)
2856        if (subchap.eq.1) then
2857           xj=xj_temp-xi
2858           yj=yj_temp-yi
2859           zj=zj_temp-zi
2860        else
2861           xj=xj_safe-xi
2862           yj=yj_safe-yi
2863           zj=zj_safe-zi
2864        endif
2865           rrij=1.0D0/(xj*xj+yj*yj+zj*zj)
2866 c          sss=sscale(1.0d0/(dsqrt(rrij)*rscp(itypj,iteli)))
2867 c          sssgrad=sscagrad(1.0d0/(dsqrt(rrij)*rscp(itypj,iteli)))
2868           sss=sscale(1.0d0/(dsqrt(rrij)))
2869           sssgrad=sscagrad(1.0d0/(dsqrt(rrij)))
2870           if (energy_dec) write (iout,*) "rrij",1.0d0/dsqrt(rrij),
2871      &     " rscp",rscp(itypj,iteli)," subchap",subchap," sss",sss
2872 c          if (lprint_short) write (iout,*) "rij",1.0/dsqrt(rrij),
2873 c     &     " subchap",subchap," sss",sss
2874           if (sss.gt.0.0d0) then
2875
2876             fac=rrij**expon2
2877             e1=fac*fac*aad(itypj,iteli)
2878             e2=fac*bad(itypj,iteli)
2879             if (iabs(j-i) .le. 2) then
2880               e1=scal14*e1
2881               e2=scal14*e2
2882               evdw2_14=evdw2_14+(e1+e2)*sss
2883             endif
2884             evdwij=e1+e2
2885             evdw2=evdw2+evdwij*sss
2886             if (energy_dec) write (iout,'(a6,2i5,2(0pf7.3))')
2887      &          'evdw2',i,j,sss,evdwij
2888 C
2889 C Calculate contributions to the gradient in the virtual-bond and SC vectors.
2890 C
2891             fac=-(evdwij+e1)*rrij*sss
2892             fac=fac+(evdwij)*sssgrad*dsqrt(rrij)/rscp(itypj,iteli)
2893             ggg(1)=xj*fac
2894             ggg(2)=yj*fac
2895             ggg(3)=zj*fac
2896 C Uncomment following three lines for SC-p interactions
2897 c           do k=1,3
2898 c             gradx_scp(k,j)=gradx_scp(k,j)+ggg(k)
2899 c           enddo
2900 C Uncomment following line for SC-p interactions
2901 c             gradx_scp(k,j)=gradx_scp(k,j)+ggg(k)
2902             do k=1,3
2903               gvdwc_scpp(k,i)=gvdwc_scpp(k,i)-ggg(k)
2904               gvdwc_scp(k,j)=gvdwc_scp(k,j)+ggg(k)
2905             enddo
2906           endif
2907         enddo
2908
2909         enddo ! iint
2910       enddo ! i
2911       do i=1,nct
2912         do j=1,3
2913           gvdwc_scp(j,i)=expon*gvdwc_scp(j,i)
2914           gvdwc_scpp(j,i)=expon*gvdwc_scpp(j,i)
2915           gradx_scp(j,i)=expon*gradx_scp(j,i)
2916         enddo
2917       enddo
2918 C******************************************************************************
2919 C
2920 C                              N O T E !!!
2921 C
2922 C To save time the factor EXPON has been extracted from ALL components
2923 C of GVDWC and GRADX. Remember to multiply them by this factor before further
2924 C use!
2925 C
2926 C******************************************************************************
2927       return
2928       end