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