source/unres/src_MD-M-SAXS.safe/energy_p_new-sep_barrier.F

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