source/unres/src_MD-M-homology/energy_p_new-sep_barrier.F

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