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