--- /dev/null
+C-----------------------------------------------------------------------
+ double precision function sscale(r)
+ double precision r,gamm
+ include "COMMON.SPLITELE"
+ if(r.lt.r_cut-rlamb) then
+ sscale=1.0d0
+ else if(r.le.r_cut.and.r.ge.r_cut-rlamb) then
+ gamm=(r-(r_cut-rlamb))/rlamb
+ sscale=1.0d0+gamm*gamm*(2*gamm-3.0d0)
+ else
+ sscale=0d0
+ endif
+ return
+ end
+C-----------------------------------------------------------------------
+ subroutine elj_long(evdw)
+C
+C This subroutine calculates the interaction energy of nonbonded side chains
+C assuming the LJ potential of interaction.
+C
+ implicit real*8 (a-h,o-z)
+ include 'DIMENSIONS'
+ parameter (accur=1.0d-10)
+ include 'COMMON.GEO'
+ include 'COMMON.VAR'
+ include 'COMMON.LOCAL'
+ include 'COMMON.CHAIN'
+ include 'COMMON.DERIV'
+ include 'COMMON.INTERACT'
+ include 'COMMON.TORSION'
+ include 'COMMON.SBRIDGE'
+ include 'COMMON.NAMES'
+ include 'COMMON.IOUNITS'
+ include 'COMMON.CONTACTS'
+ dimension gg(3)
+c write(iout,*)'Entering ELJ nnt=',nnt,' nct=',nct,' expon=',expon
+ evdw=0.0D0
+ do i=iatsc_s,iatsc_e
+ itypi=itype(i)
+ if (itypi.eq.ntyp1) cycle
+ itypi1=itype(i+1)
+ xi=c(1,nres+i)
+ yi=c(2,nres+i)
+ zi=c(3,nres+i)
+C
+C Calculate SC interaction energy.
+C
+ do iint=1,nint_gr(i)
+cd write (iout,*) 'i=',i,' iint=',iint,' istart=',istart(i,iint),
+cd & 'iend=',iend(i,iint)
+ do j=istart(i,iint),iend(i,iint)
+ itypj=itype(j)
+ if (itypj.eq.ntyp1) cycle
+ xj=c(1,nres+j)-xi
+ yj=c(2,nres+j)-yi
+ zj=c(3,nres+j)-zi
+ rij=xj*xj+yj*yj+zj*zj
+ sss=sscale(dsqrt(rij)/sigma(itypi,itypj))
+ if (sss.lt.1.0d0) then
+ rrij=1.0D0/rij
+ eps0ij=eps(itypi,itypj)
+ fac=rrij**expon2
+ e1=fac*fac*aa(itypi,itypj)
+ e2=fac*bb(itypi,itypj)
+ evdwij=e1+e2
+ evdw=evdw+(1.0d0-sss)*evdwij
+C
+C Calculate the components of the gradient in DC and X
+C
+ fac=-rrij*(e1+evdwij)*(1.0d0-sss)
+ gg(1)=xj*fac
+ gg(2)=yj*fac
+ gg(3)=zj*fac
+ do k=1,3
+ gvdwx(k,i)=gvdwx(k,i)-gg(k)
+ gvdwx(k,j)=gvdwx(k,j)+gg(k)
+ gvdwc(k,i)=gvdwc(k,i)-gg(k)
+ gvdwc(k,j)=gvdwc(k,j)+gg(k)
+ enddo
+ endif
+ enddo ! j
+ enddo ! iint
+ enddo ! i
+ do i=1,nct
+ do j=1,3
+ gvdwc(j,i)=expon*gvdwc(j,i)
+ gvdwx(j,i)=expon*gvdwx(j,i)
+ enddo
+ enddo
+C******************************************************************************
+C
+C N O T E !!!
+C
+C To save time, the factor of EXPON has been extracted from ALL components
+C of GVDWC and GRADX. Remember to multiply them by this factor before further
+C use!
+C
+C******************************************************************************
+ return
+ end
+C-----------------------------------------------------------------------
+ subroutine elj_short(evdw)
+C
+C This subroutine calculates the interaction energy of nonbonded side chains
+C assuming the LJ potential of interaction.
+C
+ implicit real*8 (a-h,o-z)
+ include 'DIMENSIONS'
+ parameter (accur=1.0d-10)
+ include 'COMMON.GEO'
+ include 'COMMON.VAR'
+ include 'COMMON.LOCAL'
+ include 'COMMON.CHAIN'
+ include 'COMMON.DERIV'
+ include 'COMMON.INTERACT'
+ include 'COMMON.TORSION'
+ include 'COMMON.SBRIDGE'
+ include 'COMMON.NAMES'
+ include 'COMMON.IOUNITS'
+ include 'COMMON.CONTACTS'
+ dimension gg(3)
+c write(iout,*)'Entering ELJ nnt=',nnt,' nct=',nct,' expon=',expon
+ evdw=0.0D0
+ do i=iatsc_s,iatsc_e
+ itypi=itype(i)
+ if (itypi.eq.ntyp1) cycle
+ itypi1=itype(i+1)
+ xi=c(1,nres+i)
+ yi=c(2,nres+i)
+ zi=c(3,nres+i)
+C Change 12/1/95
+ num_conti=0
+C
+C Calculate SC interaction energy.
+C
+ do iint=1,nint_gr(i)
+cd write (iout,*) 'i=',i,' iint=',iint,' istart=',istart(i,iint),
+cd & 'iend=',iend(i,iint)
+ do j=istart(i,iint),iend(i,iint)
+ itypj=itype(j)
+ if (itypj.eq.ntyp1) cycle
+ xj=c(1,nres+j)-xi
+ yj=c(2,nres+j)-yi
+ zj=c(3,nres+j)-zi
+C Change 12/1/95 to calculate four-body interactions
+ rij=xj*xj+yj*yj+zj*zj
+ sss=sscale(dsqrt(rij)/sigma(itypi,itypj))
+ if (sss.gt.0.0d0) then
+ rrij=1.0D0/rij
+ eps0ij=eps(itypi,itypj)
+ fac=rrij**expon2
+ e1=fac*fac*aa(itypi,itypj)
+ e2=fac*bb(itypi,itypj)
+ evdwij=e1+e2
+ evdw=evdw+sss*evdwij
+C
+C Calculate the components of the gradient in DC and X
+C
+ fac=-rrij*(e1+evdwij)*sss
+ gg(1)=xj*fac
+ gg(2)=yj*fac
+ gg(3)=zj*fac
+ do k=1,3
+ gvdwx(k,i)=gvdwx(k,i)-gg(k)
+ gvdwx(k,j)=gvdwx(k,j)+gg(k)
+ gvdwc(k,i)=gvdwc(k,i)-gg(k)
+ gvdwc(k,j)=gvdwc(k,j)+gg(k)
+ enddo
+ endif
+ enddo ! j
+ enddo ! iint
+ enddo ! i
+ do i=1,nct
+ do j=1,3
+ gvdwc(j,i)=expon*gvdwc(j,i)
+ gvdwx(j,i)=expon*gvdwx(j,i)
+ enddo
+ enddo
+C******************************************************************************
+C
+C N O T E !!!
+C
+C To save time, the factor of EXPON has been extracted from ALL components
+C of GVDWC and GRADX. Remember to multiply them by this factor before further
+C use!
+C
+C******************************************************************************
+ return
+ end
+C-----------------------------------------------------------------------------
+ subroutine eljk_long(evdw)
+C
+C This subroutine calculates the interaction energy of nonbonded side chains
+C assuming the LJK potential of interaction.
+C
+ implicit real*8 (a-h,o-z)
+ include 'DIMENSIONS'
+ include 'COMMON.GEO'
+ include 'COMMON.VAR'
+ include 'COMMON.LOCAL'
+ include 'COMMON.CHAIN'
+ include 'COMMON.DERIV'
+ include 'COMMON.INTERACT'
+ include 'COMMON.IOUNITS'
+ include 'COMMON.NAMES'
+ dimension gg(3)
+ logical scheck
+c print *,'Entering ELJK nnt=',nnt,' nct=',nct,' expon=',expon
+ evdw=0.0D0
+ do i=iatsc_s,iatsc_e
+ itypi=itype(i)
+ if (itypi.eq.ntyp1) cycle
+ itypi1=itype(i+1)
+ xi=c(1,nres+i)
+ yi=c(2,nres+i)
+ zi=c(3,nres+i)
+C
+C Calculate SC interaction energy.
+C
+ do iint=1,nint_gr(i)
+ do j=istart(i,iint),iend(i,iint)
+ itypj=itype(j)
+ if (itypj.eq.ntyp1) cycle
+ xj=c(1,nres+j)-xi
+ yj=c(2,nres+j)-yi
+ zj=c(3,nres+j)-zi
+ rrij=1.0D0/(xj*xj+yj*yj+zj*zj)
+ fac_augm=rrij**expon
+ e_augm=augm(itypi,itypj)*fac_augm
+ r_inv_ij=dsqrt(rrij)
+ rij=1.0D0/r_inv_ij
+ sss=sscale(rij/sigma(itypi,itypj))
+ if (sss.lt.1.0d0) then
+ r_shift_inv=1.0D0/(rij+r0(itypi,itypj)-sigma(itypi,itypj))
+ fac=r_shift_inv**expon
+ e1=fac*fac*aa(itypi,itypj)
+ e2=fac*bb(itypi,itypj)
+ evdwij=e_augm+e1+e2
+cd sigm=dabs(aa(itypi,itypj)/bb(itypi,itypj))**(1.0D0/6.0D0)
+cd epsi=bb(itypi,itypj)**2/aa(itypi,itypj)
+cd write (iout,'(2(a3,i3,2x),8(1pd12.4)/2(3(1pd12.4),5x)/)')
+cd & restyp(itypi),i,restyp(itypj),j,aa(itypi,itypj),
+cd & bb(itypi,itypj),augm(itypi,itypj),epsi,sigm,
+cd & sigma(itypi,itypj),1.0D0/dsqrt(rrij),evdwij,
+cd & (c(k,i),k=1,3),(c(k,j),k=1,3)
+ evdw=evdw+(1.0d0-sss)*evdwij
+C
+C Calculate the components of the gradient in DC and X
+C
+ fac=-2.0D0*rrij*e_augm-r_inv_ij*r_shift_inv*(e1+e1+e2)
+ fac=fac*(1.0d0-sss)
+ gg(1)=xj*fac
+ gg(2)=yj*fac
+ gg(3)=zj*fac
+ do k=1,3
+ gvdwx(k,i)=gvdwx(k,i)-gg(k)
+ gvdwx(k,j)=gvdwx(k,j)+gg(k)
+ gvdwc(k,i)=gvdwc(k,i)-gg(k)
+ gvdwc(k,j)=gvdwc(k,j)+gg(k)
+ enddo
+ endif
+ enddo ! j
+ enddo ! iint
+ enddo ! i
+ do i=1,nct
+ do j=1,3
+ gvdwc(j,i)=expon*gvdwc(j,i)
+ gvdwx(j,i)=expon*gvdwx(j,i)
+ enddo
+ enddo
+ return
+ end
+C-----------------------------------------------------------------------------
+ subroutine eljk_short(evdw)
+C
+C This subroutine calculates the interaction energy of nonbonded side chains
+C assuming the LJK potential of interaction.
+C
+ implicit real*8 (a-h,o-z)
+ include 'DIMENSIONS'
+ include 'COMMON.GEO'
+ include 'COMMON.VAR'
+ include 'COMMON.LOCAL'
+ include 'COMMON.CHAIN'
+ include 'COMMON.DERIV'
+ include 'COMMON.INTERACT'
+ include 'COMMON.IOUNITS'
+ include 'COMMON.NAMES'
+ dimension gg(3)
+ logical scheck
+c print *,'Entering ELJK nnt=',nnt,' nct=',nct,' expon=',expon
+ evdw=0.0D0
+ do i=iatsc_s,iatsc_e
+ itypi=itype(i)
+ if (itypi.eq.ntyp1) cycle
+ itypi1=itype(i+1)
+ xi=c(1,nres+i)
+ yi=c(2,nres+i)
+ zi=c(3,nres+i)
+C
+C Calculate SC interaction energy.
+C
+ do iint=1,nint_gr(i)
+ do j=istart(i,iint),iend(i,iint)
+ itypj=itype(j)
+ if (itypj.eq.ntyp1) cycle
+ xj=c(1,nres+j)-xi
+ yj=c(2,nres+j)-yi
+ zj=c(3,nres+j)-zi
+ rrij=1.0D0/(xj*xj+yj*yj+zj*zj)
+ fac_augm=rrij**expon
+ e_augm=augm(itypi,itypj)*fac_augm
+ r_inv_ij=dsqrt(rrij)
+ rij=1.0D0/r_inv_ij
+ sss=sscale(rij/sigma(itypi,itypj))
+ if (sss.gt.0.0d0) then
+ r_shift_inv=1.0D0/(rij+r0(itypi,itypj)-sigma(itypi,itypj))
+ fac=r_shift_inv**expon
+ e1=fac*fac*aa(itypi,itypj)
+ e2=fac*bb(itypi,itypj)
+ evdwij=e_augm+e1+e2
+cd sigm=dabs(aa(itypi,itypj)/bb(itypi,itypj))**(1.0D0/6.0D0)
+cd epsi=bb(itypi,itypj)**2/aa(itypi,itypj)
+cd write (iout,'(2(a3,i3,2x),8(1pd12.4)/2(3(1pd12.4),5x)/)')
+cd & restyp(itypi),i,restyp(itypj),j,aa(itypi,itypj),
+cd & bb(itypi,itypj),augm(itypi,itypj),epsi,sigm,
+cd & sigma(itypi,itypj),1.0D0/dsqrt(rrij),evdwij,
+cd & (c(k,i),k=1,3),(c(k,j),k=1,3)
+ evdw=evdw+sss*evdwij
+C
+C Calculate the components of the gradient in DC and X
+C
+ fac=-2.0D0*rrij*e_augm-r_inv_ij*r_shift_inv*(e1+e1+e2)
+ fac=fac*sss
+ gg(1)=xj*fac
+ gg(2)=yj*fac
+ gg(3)=zj*fac
+ do k=1,3
+ gvdwx(k,i)=gvdwx(k,i)-gg(k)
+ gvdwx(k,j)=gvdwx(k,j)+gg(k)
+ gvdwc(k,i)=gvdwc(k,i)-gg(k)
+ gvdwc(k,j)=gvdwc(k,j)+gg(k)
+ enddo
+ endif
+ enddo ! j
+ enddo ! iint
+ enddo ! i
+ do i=1,nct
+ do j=1,3
+ gvdwc(j,i)=expon*gvdwc(j,i)
+ gvdwx(j,i)=expon*gvdwx(j,i)
+ enddo
+ enddo
+ return
+ end
+C-----------------------------------------------------------------------------
+ subroutine ebp_long(evdw)
+C
+C This subroutine calculates the interaction energy of nonbonded side chains
+C assuming the Berne-Pechukas potential of interaction.
+C
+ implicit real*8 (a-h,o-z)
+ include 'DIMENSIONS'
+ include 'COMMON.GEO'
+ include 'COMMON.VAR'
+ include 'COMMON.LOCAL'
+ include 'COMMON.CHAIN'
+ include 'COMMON.DERIV'
+ include 'COMMON.NAMES'
+ include 'COMMON.INTERACT'
+ include 'COMMON.IOUNITS'
+ include 'COMMON.CALC'
+ common /srutu/ icall
+c double precision rrsave(maxdim)
+ logical lprn
+ evdw=0.0D0
+c print *,'Entering EBP nnt=',nnt,' nct=',nct,' expon=',expon
+ evdw=0.0D0
+c if (icall.eq.0) then
+c lprn=.true.
+c else
+ lprn=.false.
+c endif
+ ind=0
+ do i=iatsc_s,iatsc_e
+ itypi=itype(i)
+ if (itypi.eq.ntyp1) cycle
+ itypi1=itype(i+1)
+ xi=c(1,nres+i)
+ yi=c(2,nres+i)
+ zi=c(3,nres+i)
+ dxi=dc_norm(1,nres+i)
+ dyi=dc_norm(2,nres+i)
+ dzi=dc_norm(3,nres+i)
+c dsci_inv=dsc_inv(itypi)
+ dsci_inv=vbld_inv(i+nres)
+C
+C Calculate SC interaction energy.
+C
+ do iint=1,nint_gr(i)
+ do j=istart(i,iint),iend(i,iint)
+ ind=ind+1
+ itypj=itype(j)
+ if (itypj.eq.ntyp1) cycle
+c dscj_inv=dsc_inv(itypj)
+ dscj_inv=vbld_inv(j+nres)
+ chi1=chi(itypi,itypj)
+ chi2=chi(itypj,itypi)
+ chi12=chi1*chi2
+ chip1=chip(itypi)
+ chip2=chip(itypj)
+ chip12=chip1*chip2
+ alf1=alp(itypi)
+ alf2=alp(itypj)
+ alf12=0.5D0*(alf1+alf2)
+ xj=c(1,nres+j)-xi
+ yj=c(2,nres+j)-yi
+ zj=c(3,nres+j)-zi
+ dxj=dc_norm(1,nres+j)
+ dyj=dc_norm(2,nres+j)
+ dzj=dc_norm(3,nres+j)
+ rrij=1.0D0/(xj*xj+yj*yj+zj*zj)
+ rij=dsqrt(rrij)
+ sss=sscale(1.0d0/(rij*sigmaii(itypi,itypj)))
+
+ if (sss.lt.1.0d0) then
+
+C Calculate the angle-dependent terms of energy & contributions to derivatives.
+ call sc_angular
+C Calculate whole angle-dependent part of epsilon and contributions
+C to its derivatives
+ fac=(rrij*sigsq)**expon2
+ e1=fac*fac*aa(itypi,itypj)
+ e2=fac*bb(itypi,itypj)
+ evdwij=eps1*eps2rt*eps3rt*(e1+e2)
+ eps2der=evdwij*eps3rt
+ eps3der=evdwij*eps2rt
+ evdwij=evdwij*eps2rt*eps3rt
+ evdw=evdw+evdwij*(1.0d0-sss)
+ if (lprn) then
+ sigm=dabs(aa(itypi,itypj)/bb(itypi,itypj))**(1.0D0/6.0D0)
+ epsi=bb(itypi,itypj)**2/aa(itypi,itypj)
+cd write (iout,'(2(a3,i3,2x),15(0pf7.3))')
+cd & restyp(itypi),i,restyp(itypj),j,
+cd & epsi,sigm,chi1,chi2,chip1,chip2,
+cd & eps1,eps2rt**2,eps3rt**2,1.0D0/dsqrt(sigsq),
+cd & om1,om2,om12,1.0D0/dsqrt(rrij),
+cd & evdwij
+ endif
+C Calculate gradient components.
+ e1=e1*eps1*eps2rt**2*eps3rt**2
+ fac=-expon*(e1+evdwij)
+ sigder=fac/sigsq
+ fac=rrij*fac
+C Calculate radial part of the gradient
+ gg(1)=xj*fac
+ gg(2)=yj*fac
+ gg(3)=zj*fac
+C Calculate the angular part of the gradient and sum add the contributions
+C to the appropriate components of the Cartesian gradient.
+ call sc_grad_scale(1.0d0-sss)
+ endif
+ enddo ! j
+ enddo ! iint
+ enddo ! i
+c stop
+ return
+ end
+C-----------------------------------------------------------------------------
+ subroutine ebp_short(evdw)
+C
+C This subroutine calculates the interaction energy of nonbonded side chains
+C assuming the Berne-Pechukas potential of interaction.
+C
+ implicit real*8 (a-h,o-z)
+ include 'DIMENSIONS'
+ include 'COMMON.GEO'
+ include 'COMMON.VAR'
+ include 'COMMON.LOCAL'
+ include 'COMMON.CHAIN'
+ include 'COMMON.DERIV'
+ include 'COMMON.NAMES'
+ include 'COMMON.INTERACT'
+ include 'COMMON.IOUNITS'
+ include 'COMMON.CALC'
+ common /srutu/ icall
+c double precision rrsave(maxdim)
+ logical lprn
+ evdw=0.0D0
+c print *,'Entering EBP nnt=',nnt,' nct=',nct,' expon=',expon
+ evdw=0.0D0
+c if (icall.eq.0) then
+c lprn=.true.
+c else
+ lprn=.false.
+c endif
+ ind=0
+ do i=iatsc_s,iatsc_e
+ itypi=itype(i)
+ if (itypi.eq.ntyp1) cycle
+ itypi1=itype(i+1)
+ xi=c(1,nres+i)
+ yi=c(2,nres+i)
+ zi=c(3,nres+i)
+ dxi=dc_norm(1,nres+i)
+ dyi=dc_norm(2,nres+i)
+ dzi=dc_norm(3,nres+i)
+c dsci_inv=dsc_inv(itypi)
+ dsci_inv=vbld_inv(i+nres)
+C
+C Calculate SC interaction energy.
+C
+ do iint=1,nint_gr(i)
+ do j=istart(i,iint),iend(i,iint)
+ ind=ind+1
+ itypj=itype(j)
+ if (itypj.eq.ntyp1) cycle
+c dscj_inv=dsc_inv(itypj)
+ dscj_inv=vbld_inv(j+nres)
+ chi1=chi(itypi,itypj)
+ chi2=chi(itypj,itypi)
+ chi12=chi1*chi2
+ chip1=chip(itypi)
+ chip2=chip(itypj)
+ chip12=chip1*chip2
+ alf1=alp(itypi)
+ alf2=alp(itypj)
+ alf12=0.5D0*(alf1+alf2)
+ xj=c(1,nres+j)-xi
+ yj=c(2,nres+j)-yi
+ zj=c(3,nres+j)-zi
+ dxj=dc_norm(1,nres+j)
+ dyj=dc_norm(2,nres+j)
+ dzj=dc_norm(3,nres+j)
+ rrij=1.0D0/(xj*xj+yj*yj+zj*zj)
+ rij=dsqrt(rrij)
+ sss=sscale(1.0d0/(rij*sigmaii(itypi,itypj)))
+
+ if (sss.gt.0.0d0) then
+
+C Calculate the angle-dependent terms of energy & contributions to derivatives.
+ call sc_angular
+C Calculate whole angle-dependent part of epsilon and contributions
+C to its derivatives
+ fac=(rrij*sigsq)**expon2
+ e1=fac*fac*aa(itypi,itypj)
+ e2=fac*bb(itypi,itypj)
+ evdwij=eps1*eps2rt*eps3rt*(e1+e2)
+ eps2der=evdwij*eps3rt
+ eps3der=evdwij*eps2rt
+ evdwij=evdwij*eps2rt*eps3rt
+ evdw=evdw+evdwij*sss
+ if (lprn) then
+ sigm=dabs(aa(itypi,itypj)/bb(itypi,itypj))**(1.0D0/6.0D0)
+ epsi=bb(itypi,itypj)**2/aa(itypi,itypj)
+cd write (iout,'(2(a3,i3,2x),15(0pf7.3))')
+cd & restyp(itypi),i,restyp(itypj),j,
+cd & epsi,sigm,chi1,chi2,chip1,chip2,
+cd & eps1,eps2rt**2,eps3rt**2,1.0D0/dsqrt(sigsq),
+cd & om1,om2,om12,1.0D0/dsqrt(rrij),
+cd & evdwij
+ endif
+C Calculate gradient components.
+ e1=e1*eps1*eps2rt**2*eps3rt**2
+ fac=-expon*(e1+evdwij)
+ sigder=fac/sigsq
+ fac=rrij*fac
+C Calculate radial part of the gradient
+ gg(1)=xj*fac
+ gg(2)=yj*fac
+ gg(3)=zj*fac
+C Calculate the angular part of the gradient and sum add the contributions
+C to the appropriate components of the Cartesian gradient.
+ call sc_grad_scale(sss)
+ endif
+ enddo ! j
+ enddo ! iint
+ enddo ! i
+c stop
+ return
+ end
+C-----------------------------------------------------------------------------
+ subroutine egb_long(evdw)
+C
+C This subroutine calculates the interaction energy of nonbonded side chains
+C assuming the Gay-Berne potential of interaction.
+C
+ implicit real*8 (a-h,o-z)
+ include 'DIMENSIONS'
+ include 'COMMON.GEO'
+ include 'COMMON.VAR'
+ include 'COMMON.LOCAL'
+ include 'COMMON.CHAIN'
+ include 'COMMON.DERIV'
+ include 'COMMON.NAMES'
+ include 'COMMON.INTERACT'
+ include 'COMMON.IOUNITS'
+ include 'COMMON.CALC'
+ include 'COMMON.CONTROL'
+ logical lprn
+ evdw=0.0D0
+ccccc energy_dec=.false.
+c print *,'Entering EGB nnt=',nnt,' nct=',nct,' expon=',expon
+ evdw=0.0D0
+ lprn=.false.
+c if (icall.eq.0) lprn=.false.
+ ind=0
+ do i=iatsc_s,iatsc_e
+ itypi=itype(i)
+ if (itypi.eq.ntyp1) cycle
+ itypi1=itype(i+1)
+ xi=c(1,nres+i)
+ yi=c(2,nres+i)
+ zi=c(3,nres+i)
+ dxi=dc_norm(1,nres+i)
+ dyi=dc_norm(2,nres+i)
+ dzi=dc_norm(3,nres+i)
+c dsci_inv=dsc_inv(itypi)
+ dsci_inv=vbld_inv(i+nres)
+c write (iout,*) "i",i,dsc_inv(itypi),dsci_inv,1.0d0/vbld(i+nres)
+c write (iout,*) "dcnori",dxi*dxi+dyi*dyi+dzi*dzi
+C
+C Calculate SC interaction energy.
+C
+ do iint=1,nint_gr(i)
+ do j=istart(i,iint),iend(i,iint)
+ ind=ind+1
+ itypj=itype(j)
+ if (itypj.eq.ntyp1) cycle
+c dscj_inv=dsc_inv(itypj)
+ dscj_inv=vbld_inv(j+nres)
+c write (iout,*) "j",j,dsc_inv(itypj),dscj_inv,
+c & 1.0d0/vbld(j+nres)
+c write (iout,*) "i",i," j", j," itype",itype(i),itype(j)
+ sig0ij=sigma(itypi,itypj)
+ chi1=chi(itypi,itypj)
+ chi2=chi(itypj,itypi)
+ chi12=chi1*chi2
+ chip1=chip(itypi)
+ chip2=chip(itypj)
+ chip12=chip1*chip2
+ alf1=alp(itypi)
+ alf2=alp(itypj)
+ alf12=0.5D0*(alf1+alf2)
+ xj=c(1,nres+j)-xi
+ yj=c(2,nres+j)-yi
+ zj=c(3,nres+j)-zi
+ dxj=dc_norm(1,nres+j)
+ dyj=dc_norm(2,nres+j)
+ dzj=dc_norm(3,nres+j)
+ rrij=1.0D0/(xj*xj+yj*yj+zj*zj)
+ rij=dsqrt(rrij)
+ sss=sscale(1.0d0/(rij*sigmaii(itypi,itypj)))
+
+ if (sss.lt.1.0d0) then
+
+C Calculate angle-dependent terms of energy and contributions to their
+C derivatives.
+ call sc_angular
+ sigsq=1.0D0/sigsq
+ sig=sig0ij*dsqrt(sigsq)
+ rij_shift=1.0D0/rij-sig+sig0ij
+c for diagnostics; uncomment
+c rij_shift=1.2*sig0ij
+C I hate to put IF's in the loops, but here don't have another choice!!!!
+ if (rij_shift.le.0.0D0) then
+ evdw=1.0D20
+cd write (iout,'(2(a3,i3,2x),17(0pf7.3))')
+cd & restyp(itypi),i,restyp(itypj),j,
+cd & rij_shift,1.0D0/rij,sig,sig0ij,sigsq,1-dsqrt(sigsq)
+ return
+ endif
+ sigder=-sig*sigsq
+c---------------------------------------------------------------
+ rij_shift=1.0D0/rij_shift
+ fac=rij_shift**expon
+ e1=fac*fac*aa(itypi,itypj)
+ e2=fac*bb(itypi,itypj)
+ evdwij=eps1*eps2rt*eps3rt*(e1+e2)
+ eps2der=evdwij*eps3rt
+ eps3der=evdwij*eps2rt
+c write (iout,*) "sigsq",sigsq," sig",sig," eps2rt",eps2rt,
+c & " eps3rt",eps3rt," eps1",eps1," e1",e1," e2",e2
+ evdwij=evdwij*eps2rt*eps3rt
+ evdw=evdw+evdwij*(1.0d0-sss)
+ if (lprn) then
+ sigm=dabs(aa(itypi,itypj)/bb(itypi,itypj))**(1.0D0/6.0D0)
+ epsi=bb(itypi,itypj)**2/aa(itypi,itypj)
+ write (iout,'(2(a3,i3,2x),17(0pf7.3))')
+ & restyp(itypi),i,restyp(itypj),j,
+ & epsi,sigm,chi1,chi2,chip1,chip2,
+ & eps1,eps2rt**2,eps3rt**2,sig,sig0ij,
+ & om1,om2,om12,1.0D0/rij,1.0D0/rij_shift,
+ & evdwij
+ endif
+
+ if (energy_dec) write (iout,'(a6,2i5,0pf7.3)')
+ & 'evdw',i,j,evdwij
+
+C Calculate gradient components.
+ e1=e1*eps1*eps2rt**2*eps3rt**2
+ fac=-expon*(e1+evdwij)*rij_shift
+ sigder=fac*sigder
+ fac=rij*fac
+c fac=0.0d0
+C Calculate the radial part of the gradient
+ gg(1)=xj*fac
+ gg(2)=yj*fac
+ gg(3)=zj*fac
+C Calculate angular part of the gradient.
+ call sc_grad_scale(1.0d0-sss)
+ endif
+ enddo ! j
+ enddo ! iint
+ enddo ! i
+c write (iout,*) "Number of loop steps in EGB:",ind
+cccc energy_dec=.false.
+ return
+ end
+C-----------------------------------------------------------------------------
+ subroutine egb_short(evdw)
+C
+C This subroutine calculates the interaction energy of nonbonded side chains
+C assuming the Gay-Berne potential of interaction.
+C
+ implicit real*8 (a-h,o-z)
+ include 'DIMENSIONS'
+ include 'COMMON.GEO'
+ include 'COMMON.VAR'
+ include 'COMMON.LOCAL'
+ include 'COMMON.CHAIN'
+ include 'COMMON.DERIV'
+ include 'COMMON.NAMES'
+ include 'COMMON.INTERACT'
+ include 'COMMON.IOUNITS'
+ include 'COMMON.CALC'
+ include 'COMMON.CONTROL'
+ logical lprn
+ evdw=0.0D0
+ccccc energy_dec=.false.
+c print *,'Entering EGB nnt=',nnt,' nct=',nct,' expon=',expon
+ evdw=0.0D0
+ lprn=.false.
+c if (icall.eq.0) lprn=.false.
+ ind=0
+ do i=iatsc_s,iatsc_e
+ itypi=itype(i)
+ if (itypi.eq.ntyp1) cycle
+ itypi1=itype(i+1)
+ xi=c(1,nres+i)
+ yi=c(2,nres+i)
+ zi=c(3,nres+i)
+ dxi=dc_norm(1,nres+i)
+ dyi=dc_norm(2,nres+i)
+ dzi=dc_norm(3,nres+i)
+c dsci_inv=dsc_inv(itypi)
+ dsci_inv=vbld_inv(i+nres)
+c write (iout,*) "i",i,dsc_inv(itypi),dsci_inv,1.0d0/vbld(i+nres)
+c write (iout,*) "dcnori",dxi*dxi+dyi*dyi+dzi*dzi
+C
+C Calculate SC interaction energy.
+C
+ do iint=1,nint_gr(i)
+ do j=istart(i,iint),iend(i,iint)
+ ind=ind+1
+ itypj=itype(j)
+ if (itypj.eq.ntyp1) cycle
+c dscj_inv=dsc_inv(itypj)
+ dscj_inv=vbld_inv(j+nres)
+c write (iout,*) "j",j,dsc_inv(itypj),dscj_inv,
+c & 1.0d0/vbld(j+nres)
+c write (iout,*) "i",i," j", j," itype",itype(i),itype(j)
+ sig0ij=sigma(itypi,itypj)
+ chi1=chi(itypi,itypj)
+ chi2=chi(itypj,itypi)
+ chi12=chi1*chi2
+ chip1=chip(itypi)
+ chip2=chip(itypj)
+ chip12=chip1*chip2
+ alf1=alp(itypi)
+ alf2=alp(itypj)
+ alf12=0.5D0*(alf1+alf2)
+ xj=c(1,nres+j)-xi
+ yj=c(2,nres+j)-yi
+ zj=c(3,nres+j)-zi
+ dxj=dc_norm(1,nres+j)
+ dyj=dc_norm(2,nres+j)
+ dzj=dc_norm(3,nres+j)
+ rrij=1.0D0/(xj*xj+yj*yj+zj*zj)
+ rij=dsqrt(rrij)
+ sss=sscale(1.0d0/(rij*sigmaii(itypi,itypj)))
+
+ if (sss.gt.0.0d0) then
+
+C Calculate angle-dependent terms of energy and contributions to their
+C derivatives.
+ call sc_angular
+ sigsq=1.0D0/sigsq
+ sig=sig0ij*dsqrt(sigsq)
+ rij_shift=1.0D0/rij-sig+sig0ij
+c for diagnostics; uncomment
+c rij_shift=1.2*sig0ij
+C I hate to put IF's in the loops, but here don't have another choice!!!!
+ if (rij_shift.le.0.0D0) then
+ evdw=1.0D20
+cd write (iout,'(2(a3,i3,2x),17(0pf7.3))')
+cd & restyp(itypi),i,restyp(itypj),j,
+cd & rij_shift,1.0D0/rij,sig,sig0ij,sigsq,1-dsqrt(sigsq)
+ return
+ endif
+ sigder=-sig*sigsq
+c---------------------------------------------------------------
+ rij_shift=1.0D0/rij_shift
+ fac=rij_shift**expon
+ e1=fac*fac*aa(itypi,itypj)
+ e2=fac*bb(itypi,itypj)
+ evdwij=eps1*eps2rt*eps3rt*(e1+e2)
+ eps2der=evdwij*eps3rt
+ eps3der=evdwij*eps2rt
+c write (iout,*) "sigsq",sigsq," sig",sig," eps2rt",eps2rt,
+c & " eps3rt",eps3rt," eps1",eps1," e1",e1," e2",e2
+ evdwij=evdwij*eps2rt*eps3rt
+ evdw=evdw+evdwij*sss
+ if (lprn) then
+ sigm=dabs(aa(itypi,itypj)/bb(itypi,itypj))**(1.0D0/6.0D0)
+ epsi=bb(itypi,itypj)**2/aa(itypi,itypj)
+ write (iout,'(2(a3,i3,2x),17(0pf7.3))')
+ & restyp(itypi),i,restyp(itypj),j,
+ & epsi,sigm,chi1,chi2,chip1,chip2,
+ & eps1,eps2rt**2,eps3rt**2,sig,sig0ij,
+ & om1,om2,om12,1.0D0/rij,1.0D0/rij_shift,
+ & evdwij
+ endif
+
+ if (energy_dec) write (iout,'(a6,2i5,0pf7.3)')
+ & 'evdw',i,j,evdwij
+
+C Calculate gradient components.
+ e1=e1*eps1*eps2rt**2*eps3rt**2
+ fac=-expon*(e1+evdwij)*rij_shift
+ sigder=fac*sigder
+ fac=rij*fac
+c fac=0.0d0
+C Calculate the radial part of the gradient
+ gg(1)=xj*fac
+ gg(2)=yj*fac
+ gg(3)=zj*fac
+C Calculate angular part of the gradient.
+ call sc_grad_scale(sss)
+ endif
+ enddo ! j
+ enddo ! iint
+ enddo ! i
+c write (iout,*) "Number of loop steps in EGB:",ind
+cccc energy_dec=.false.
+ return
+ end
+C-----------------------------------------------------------------------------
+ subroutine egbv_long(evdw)
+C
+C This subroutine calculates the interaction energy of nonbonded side chains
+C assuming the Gay-Berne-Vorobjev potential of interaction.
+C
+ implicit real*8 (a-h,o-z)
+ include 'DIMENSIONS'
+ include 'COMMON.GEO'
+ include 'COMMON.VAR'
+ include 'COMMON.LOCAL'
+ include 'COMMON.CHAIN'
+ include 'COMMON.DERIV'
+ include 'COMMON.NAMES'
+ include 'COMMON.INTERACT'
+ include 'COMMON.IOUNITS'
+ include 'COMMON.CALC'
+ common /srutu/ icall
+ logical lprn
+ evdw=0.0D0
+c print *,'Entering EGB nnt=',nnt,' nct=',nct,' expon=',expon
+ evdw=0.0D0
+ lprn=.false.
+c if (icall.eq.0) lprn=.true.
+ ind=0
+ do i=iatsc_s,iatsc_e
+ itypi=itype(i)
+ if (itypi.eq.ntyp1) cycle
+ itypi1=itype(i+1)
+ xi=c(1,nres+i)
+ yi=c(2,nres+i)
+ zi=c(3,nres+i)
+ dxi=dc_norm(1,nres+i)
+ dyi=dc_norm(2,nres+i)
+ dzi=dc_norm(3,nres+i)
+c dsci_inv=dsc_inv(itypi)
+ dsci_inv=vbld_inv(i+nres)
+C
+C Calculate SC interaction energy.
+C
+ do iint=1,nint_gr(i)
+ do j=istart(i,iint),iend(i,iint)
+ ind=ind+1
+ itypj=itype(j)
+ if (itypj.eq.ntyp1) cycle
+c dscj_inv=dsc_inv(itypj)
+ dscj_inv=vbld_inv(j+nres)
+ sig0ij=sigma(itypi,itypj)
+ r0ij=r0(itypi,itypj)
+ chi1=chi(itypi,itypj)
+ chi2=chi(itypj,itypi)
+ chi12=chi1*chi2
+ chip1=chip(itypi)
+ chip2=chip(itypj)
+ chip12=chip1*chip2
+ alf1=alp(itypi)
+ alf2=alp(itypj)
+ alf12=0.5D0*(alf1+alf2)
+ xj=c(1,nres+j)-xi
+ yj=c(2,nres+j)-yi
+ zj=c(3,nres+j)-zi
+ dxj=dc_norm(1,nres+j)
+ dyj=dc_norm(2,nres+j)
+ dzj=dc_norm(3,nres+j)
+ rrij=1.0D0/(xj*xj+yj*yj+zj*zj)
+ rij=dsqrt(rrij)
+
+ sss=sscale(1.0d0/(rij*sigmaii(itypi,itypj)))
+
+ if (sss.lt.1.0d0) then
+
+C Calculate angle-dependent terms of energy and contributions to their
+C derivatives.
+ call sc_angular
+ sigsq=1.0D0/sigsq
+ sig=sig0ij*dsqrt(sigsq)
+ rij_shift=1.0D0/rij-sig+r0ij
+C I hate to put IF's in the loops, but here don't have another choice!!!!
+ if (rij_shift.le.0.0D0) then
+ evdw=1.0D20
+ return
+ endif
+ sigder=-sig*sigsq
+c---------------------------------------------------------------
+ rij_shift=1.0D0/rij_shift
+ fac=rij_shift**expon
+ e1=fac*fac*aa(itypi,itypj)
+ e2=fac*bb(itypi,itypj)
+ evdwij=eps1*eps2rt*eps3rt*(e1+e2)
+ eps2der=evdwij*eps3rt
+ eps3der=evdwij*eps2rt
+ fac_augm=rrij**expon
+ e_augm=augm(itypi,itypj)*fac_augm
+ evdwij=evdwij*eps2rt*eps3rt
+ evdw=evdw+(evdwij+e_augm)*(1.0d0-sss)
+ if (lprn) then
+ sigm=dabs(aa(itypi,itypj)/bb(itypi,itypj))**(1.0D0/6.0D0)
+ epsi=bb(itypi,itypj)**2/aa(itypi,itypj)
+ write (iout,'(2(a3,i3,2x),17(0pf7.3))')
+ & restyp(itypi),i,restyp(itypj),j,
+ & epsi,sigm,sig,(augm(itypi,itypj)/epsi)**(1.0D0/12.0D0),
+ & chi1,chi2,chip1,chip2,
+ & eps1,eps2rt**2,eps3rt**2,
+ & om1,om2,om12,1.0D0/rij,1.0D0/rij_shift,
+ & evdwij+e_augm
+ endif
+C Calculate gradient components.
+ e1=e1*eps1*eps2rt**2*eps3rt**2
+ fac=-expon*(e1+evdwij)*rij_shift
+ sigder=fac*sigder
+ fac=rij*fac-2*expon*rrij*e_augm
+C Calculate the radial part of the gradient
+ gg(1)=xj*fac
+ gg(2)=yj*fac
+ gg(3)=zj*fac
+C Calculate angular part of the gradient.
+ call sc_grad_scale(1.0d0-sss)
+ endif
+ enddo ! j
+ enddo ! iint
+ enddo ! i
+ end
+C-----------------------------------------------------------------------------
+ subroutine egbv_short(evdw)
+C
+C This subroutine calculates the interaction energy of nonbonded side chains
+C assuming the Gay-Berne-Vorobjev potential of interaction.
+C
+ implicit real*8 (a-h,o-z)
+ include 'DIMENSIONS'
+ include 'COMMON.GEO'
+ include 'COMMON.VAR'
+ include 'COMMON.LOCAL'
+ include 'COMMON.CHAIN'
+ include 'COMMON.DERIV'
+ include 'COMMON.NAMES'
+ include 'COMMON.INTERACT'
+ include 'COMMON.IOUNITS'
+ include 'COMMON.CALC'
+ common /srutu/ icall
+ logical lprn
+ evdw=0.0D0
+c print *,'Entering EGB nnt=',nnt,' nct=',nct,' expon=',expon
+ evdw=0.0D0
+ lprn=.false.
+c if (icall.eq.0) lprn=.true.
+ ind=0
+ do i=iatsc_s,iatsc_e
+ itypi=itype(i)
+ if (itypi.eq.ntyp1) cycle
+ itypi1=itype(i+1)
+ xi=c(1,nres+i)
+ yi=c(2,nres+i)
+ zi=c(3,nres+i)
+ dxi=dc_norm(1,nres+i)
+ dyi=dc_norm(2,nres+i)
+ dzi=dc_norm(3,nres+i)
+c dsci_inv=dsc_inv(itypi)
+ dsci_inv=vbld_inv(i+nres)
+C
+C Calculate SC interaction energy.
+C
+ do iint=1,nint_gr(i)
+ do j=istart(i,iint),iend(i,iint)
+ ind=ind+1
+ itypj=itype(j)
+ if (itypj.eq.ntyp1) cycle
+c dscj_inv=dsc_inv(itypj)
+ dscj_inv=vbld_inv(j+nres)
+ sig0ij=sigma(itypi,itypj)
+ r0ij=r0(itypi,itypj)
+ chi1=chi(itypi,itypj)
+ chi2=chi(itypj,itypi)
+ chi12=chi1*chi2
+ chip1=chip(itypi)
+ chip2=chip(itypj)
+ chip12=chip1*chip2
+ alf1=alp(itypi)
+ alf2=alp(itypj)
+ alf12=0.5D0*(alf1+alf2)
+ xj=c(1,nres+j)-xi
+ yj=c(2,nres+j)-yi
+ zj=c(3,nres+j)-zi
+ dxj=dc_norm(1,nres+j)
+ dyj=dc_norm(2,nres+j)
+ dzj=dc_norm(3,nres+j)
+ rrij=1.0D0/(xj*xj+yj*yj+zj*zj)
+ rij=dsqrt(rrij)
+
+ sss=sscale(1.0d0/(rij*sigmaii(itypi,itypj)))
+
+ if (sss.gt.0.0d0) then
+
+C Calculate angle-dependent terms of energy and contributions to their
+C derivatives.
+ call sc_angular
+ sigsq=1.0D0/sigsq
+ sig=sig0ij*dsqrt(sigsq)
+ rij_shift=1.0D0/rij-sig+r0ij
+C I hate to put IF's in the loops, but here don't have another choice!!!!
+ if (rij_shift.le.0.0D0) then
+ evdw=1.0D20
+ return
+ endif
+ sigder=-sig*sigsq
+c---------------------------------------------------------------
+ rij_shift=1.0D0/rij_shift
+ fac=rij_shift**expon
+ e1=fac*fac*aa(itypi,itypj)
+ e2=fac*bb(itypi,itypj)
+ evdwij=eps1*eps2rt*eps3rt*(e1+e2)
+ eps2der=evdwij*eps3rt
+ eps3der=evdwij*eps2rt
+ fac_augm=rrij**expon
+ e_augm=augm(itypi,itypj)*fac_augm
+ evdwij=evdwij*eps2rt*eps3rt
+ evdw=evdw+(evdwij+e_augm)*sss
+ if (lprn) then
+ sigm=dabs(aa(itypi,itypj)/bb(itypi,itypj))**(1.0D0/6.0D0)
+ epsi=bb(itypi,itypj)**2/aa(itypi,itypj)
+ write (iout,'(2(a3,i3,2x),17(0pf7.3))')
+ & restyp(itypi),i,restyp(itypj),j,
+ & epsi,sigm,sig,(augm(itypi,itypj)/epsi)**(1.0D0/12.0D0),
+ & chi1,chi2,chip1,chip2,
+ & eps1,eps2rt**2,eps3rt**2,
+ & om1,om2,om12,1.0D0/rij,1.0D0/rij_shift,
+ & evdwij+e_augm
+ endif
+C Calculate gradient components.
+ e1=e1*eps1*eps2rt**2*eps3rt**2
+ fac=-expon*(e1+evdwij)*rij_shift
+ sigder=fac*sigder
+ fac=rij*fac-2*expon*rrij*e_augm
+C Calculate the radial part of the gradient
+ gg(1)=xj*fac
+ gg(2)=yj*fac
+ gg(3)=zj*fac
+C Calculate angular part of the gradient.
+ call sc_grad_scale(sss)
+ endif
+ enddo ! j
+ enddo ! iint
+ enddo ! i
+ end
+C----------------------------------------------------------------------------
+ subroutine sc_grad_scale(scalfac)
+ implicit real*8 (a-h,o-z)
+ include 'DIMENSIONS'
+ include 'COMMON.CHAIN'
+ include 'COMMON.DERIV'
+ include 'COMMON.CALC'
+ include 'COMMON.IOUNITS'
+ double precision dcosom1(3),dcosom2(3)
+ double precision scalfac
+ eom1=eps2der*eps2rt_om1-2.0D0*alf1*eps3der+sigder*sigsq_om1
+ eom2=eps2der*eps2rt_om2+2.0D0*alf2*eps3der+sigder*sigsq_om2
+ eom12=evdwij*eps1_om12+eps2der*eps2rt_om12
+ & -2.0D0*alf12*eps3der+sigder*sigsq_om12
+c diagnostics only
+c eom1=0.0d0
+c eom2=0.0d0
+c eom12=evdwij*eps1_om12
+c end diagnostics
+c write (iout,*) "eps2der",eps2der," eps3der",eps3der,
+c & " sigder",sigder
+c write (iout,*) "eps1_om12",eps1_om12," eps2rt_om12",eps2rt_om12
+c write (iout,*) "eom1",eom1," eom2",eom2," eom12",eom12
+ do k=1,3
+ dcosom1(k)=rij*(dc_norm(k,nres+i)-om1*erij(k))
+ dcosom2(k)=rij*(dc_norm(k,nres+j)-om2*erij(k))
+ enddo
+ do k=1,3
+ gg(k)=(gg(k)+eom1*dcosom1(k)+eom2*dcosom2(k))*scalfac
+ enddo
+c write (iout,*) "gg",(gg(k),k=1,3)
+ do k=1,3
+ gvdwx(k,i)=gvdwx(k,i)-gg(k)
+ & +(eom12*(dc_norm(k,nres+j)-om12*dc_norm(k,nres+i))
+ & +eom1*(erij(k)-om1*dc_norm(k,nres+i)))*dsci_inv*scalfac
+ gvdwx(k,j)=gvdwx(k,j)+gg(k)
+ & +(eom12*(dc_norm(k,nres+i)-om12*dc_norm(k,nres+j))
+ & +eom2*(erij(k)-om2*dc_norm(k,nres+j)))*dscj_inv*scalfac
+c write (iout,*)(eom12*(dc_norm(k,nres+j)-om12*dc_norm(k,nres+i))
+c & +eom1*(erij(k)-om1*dc_norm(k,nres+i)))*dsci_inv
+c write (iout,*)(eom12*(dc_norm(k,nres+i)-om12*dc_norm(k,nres+j))
+c & +eom2*(erij(k)-om2*dc_norm(k,nres+j)))*dscj_inv
+ enddo
+C
+C Calculate the components of the gradient in DC and X
+C
+ do l=1,3
+ gvdwc(l,i)=gvdwc(l,i)-gg(l)
+ gvdwc(l,j)=gvdwc(l,j)+gg(l)
+ enddo
+ return
+ end
+C--------------------------------------------------------------------------
+ subroutine eelec_scale(ees,evdw1,eel_loc,eello_turn3,eello_turn4)
+C
+C This subroutine calculates the average interaction energy and its gradient
+C in the virtual-bond vectors between non-adjacent peptide groups, based on
+C the potential described in Liwo et al., Protein Sci., 1993, 2, 1715.
+C The potential depends both on the distance of peptide-group centers and on
+C the orientation of the CA-CA virtual bonds.
+C
+ implicit real*8 (a-h,o-z)
+#ifdef MPI
+ include 'mpif.h'
+#endif
+ include 'DIMENSIONS'
+ include 'COMMON.CONTROL'
+ include 'COMMON.SETUP'
+ include 'COMMON.IOUNITS'
+ include 'COMMON.GEO'
+ include 'COMMON.VAR'
+ include 'COMMON.LOCAL'
+ include 'COMMON.CHAIN'
+ include 'COMMON.DERIV'
+ include 'COMMON.INTERACT'
+ include 'COMMON.CONTACTS'
+ include 'COMMON.TORSION'
+ include 'COMMON.VECTORS'
+ include 'COMMON.FFIELD'
+ include 'COMMON.TIME1'
+ dimension ggg(3),gggp(3),gggm(3),erij(3),dcosb(3),dcosg(3),
+ & erder(3,3),uryg(3,3),urzg(3,3),vryg(3,3),vrzg(3,3)
+ double precision acipa(2,2),agg(3,4),aggi(3,4),aggi1(3,4),
+ & aggj(3,4),aggj1(3,4),a_temp(2,2),muij(4)
+ common /locel/ a_temp,agg,aggi,aggi1,aggj,aggj1,a22,a23,a32,a33,
+ & dxi,dyi,dzi,dx_normi,dy_normi,dz_normi,xmedi,ymedi,zmedi,
+ & num_conti,j1,j2
+c 4/26/02 - AL scaling factor for 1,4 repulsive VDW interactions
+#ifdef MOMENT
+ double precision scal_el /1.0d0/
+#else
+ double precision scal_el /0.5d0/
+#endif
+C 12/13/98
+C 13-go grudnia roku pamietnego...
+ double precision unmat(3,3) /1.0d0,0.0d0,0.0d0,
+ & 0.0d0,1.0d0,0.0d0,
+ & 0.0d0,0.0d0,1.0d0/
+cd write(iout,*) 'In EELEC'
+cd do i=1,nloctyp
+cd write(iout,*) 'Type',i
+cd write(iout,*) 'B1',B1(:,i)
+cd write(iout,*) 'B2',B2(:,i)
+cd write(iout,*) 'CC',CC(:,:,i)
+cd write(iout,*) 'DD',DD(:,:,i)
+cd write(iout,*) 'EE',EE(:,:,i)
+cd enddo
+cd call check_vecgrad
+cd stop
+ if (icheckgrad.eq.1) then
+ do i=1,nres-1
+ fac=1.0d0/dsqrt(scalar(dc(1,i),dc(1,i)))
+ do k=1,3
+ dc_norm(k,i)=dc(k,i)*fac
+ enddo
+c write (iout,*) 'i',i,' fac',fac
+ enddo
+ endif
+ if (wel_loc.gt.0.0d0 .or. wcorr4.gt.0.0d0 .or. wcorr5.gt.0.0d0
+ & .or. wcorr6.gt.0.0d0 .or. wturn3.gt.0.0d0 .or.
+ & wturn4.gt.0.0d0 .or. wturn6.gt.0.0d0) then
+c call vec_and_deriv
+#ifdef TIMING
+ time01=MPI_Wtime()
+#endif
+ call set_matrices
+#ifdef TIMING
+ time_mat=time_mat+MPI_Wtime()-time01
+#endif
+ endif
+cd do i=1,nres-1
+cd write (iout,*) 'i=',i
+cd do k=1,3
+cd write (iout,'(i5,2f10.5)') k,uy(k,i),uz(k,i)
+cd enddo
+cd do k=1,3
+cd write (iout,'(f10.5,2x,3f10.5,2x,3f10.5)')
+cd & uz(k,i),(uzgrad(k,l,1,i),l=1,3),(uzgrad(k,l,2,i),l=1,3)
+cd enddo
+cd enddo
+ t_eelecij=0.0d0
+ ees=0.0D0
+ evdw1=0.0D0
+ eel_loc=0.0d0
+ eello_turn3=0.0d0
+ eello_turn4=0.0d0
+ ind=0
+ do i=1,nres
+ num_cont_hb(i)=0
+ enddo
+cd print '(a)','Enter EELEC'
+cd write (iout,*) 'iatel_s=',iatel_s,' iatel_e=',iatel_e
+ do i=1,nres
+ gel_loc_loc(i)=0.0d0
+ gcorr_loc(i)=0.0d0
+ enddo
+c
+c
+c 9/27/08 AL Split the interaction loop to ensure load balancing of turn terms
+C
+C Loop over i,i+2 and i,i+3 pairs of the peptide groups
+C
+ do i=iturn3_start,iturn3_end
+ if (itype(i).eq.ntyp1.or. itype(i+1).eq.ntyp1
+ & .or. itype(i+2).eq.ntyp1 .or. itype(i+3).eq.ntyp1) cycle
+ dxi=dc(1,i)
+ dyi=dc(2,i)
+ dzi=dc(3,i)
+ dx_normi=dc_norm(1,i)
+ dy_normi=dc_norm(2,i)
+ dz_normi=dc_norm(3,i)
+ xmedi=c(1,i)+0.5d0*dxi
+ ymedi=c(2,i)+0.5d0*dyi
+ zmedi=c(3,i)+0.5d0*dzi
+ num_conti=0
+ call eelecij_scale(i,i+2,ees,evdw1,eel_loc)
+ if (wturn3.gt.0.0d0) call eturn3(i,eello_turn3)
+ num_cont_hb(i)=num_conti
+ enddo
+ do i=iturn4_start,iturn4_end
+ if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1
+ & .or. itype(i+3).eq.ntyp1
+ & .or. itype(i+4).eq.ntyp1) cycle
+ dxi=dc(1,i)
+ dyi=dc(2,i)
+ dzi=dc(3,i)
+ dx_normi=dc_norm(1,i)
+ dy_normi=dc_norm(2,i)
+ dz_normi=dc_norm(3,i)
+ xmedi=c(1,i)+0.5d0*dxi
+ ymedi=c(2,i)+0.5d0*dyi
+ zmedi=c(3,i)+0.5d0*dzi
+ num_conti=num_cont_hb(i)
+ call eelecij_scale(i,i+3,ees,evdw1,eel_loc)
+ if (wturn4.gt.0.0d0 .and. itype(i+2).ne.ntyp1)
+ & call eturn4(i,eello_turn4)
+ num_cont_hb(i)=num_conti
+ enddo ! i
+c
+c Loop over all pairs of interacting peptide groups except i,i+2 and i,i+3
+c
+ do i=iatel_s,iatel_e
+ if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1) cycle
+ dxi=dc(1,i)
+ dyi=dc(2,i)
+ dzi=dc(3,i)
+ dx_normi=dc_norm(1,i)
+ dy_normi=dc_norm(2,i)
+ dz_normi=dc_norm(3,i)
+ xmedi=c(1,i)+0.5d0*dxi
+ ymedi=c(2,i)+0.5d0*dyi
+ zmedi=c(3,i)+0.5d0*dzi
+c write (iout,*) 'i',i,' ielstart',ielstart(i),' ielend',ielend(i)
+ num_conti=num_cont_hb(i)
+ do j=ielstart(i),ielend(i)
+ if (itype(j).eq.ntyp1 .or. itype(j+1).eq.ntyp1) cycle
+ call eelecij_scale(i,j,ees,evdw1,eel_loc)
+ enddo ! j
+ num_cont_hb(i)=num_conti
+ enddo ! i
+c write (iout,*) "Number of loop steps in EELEC:",ind
+cd do i=1,nres
+cd write (iout,'(i3,3f10.5,5x,3f10.5)')
+cd & i,(gel_loc(k,i),k=1,3),gel_loc_loc(i)
+cd enddo
+c 12/7/99 Adam eello_turn3 will be considered as a separate energy term
+ccc eel_loc=eel_loc+eello_turn3
+cd print *,"Processor",fg_rank," t_eelecij",t_eelecij
+ return
+ end
+C-------------------------------------------------------------------------------
+ subroutine eelecij_scale(i,j,ees,evdw1,eel_loc)
+ implicit real*8 (a-h,o-z)
+ include 'DIMENSIONS'
+#ifdef MPI
+ include "mpif.h"
+#endif
+ include 'COMMON.CONTROL'
+ include 'COMMON.IOUNITS'
+ include 'COMMON.GEO'
+ include 'COMMON.VAR'
+ include 'COMMON.LOCAL'
+ include 'COMMON.CHAIN'
+ include 'COMMON.DERIV'
+ include 'COMMON.INTERACT'
+ include 'COMMON.CONTACTS'
+ include 'COMMON.TORSION'
+ include 'COMMON.VECTORS'
+ include 'COMMON.FFIELD'
+ include 'COMMON.TIME1'
+ dimension ggg(3),gggp(3),gggm(3),erij(3),dcosb(3),dcosg(3),
+ & erder(3,3),uryg(3,3),urzg(3,3),vryg(3,3),vrzg(3,3)
+ double precision acipa(2,2),agg(3,4),aggi(3,4),aggi1(3,4),
+ & aggj(3,4),aggj1(3,4),a_temp(2,2),muij(4)
+ common /locel/ a_temp,agg,aggi,aggi1,aggj,aggj1,a22,a23,a32,a33,
+ & dxi,dyi,dzi,dx_normi,dy_normi,dz_normi,xmedi,ymedi,zmedi,
+ & num_conti,j1,j2
+c 4/26/02 - AL scaling factor for 1,4 repulsive VDW interactions
+#ifdef MOMENT
+ double precision scal_el /1.0d0/
+#else
+ double precision scal_el /0.5d0/
+#endif
+C 12/13/98
+C 13-go grudnia roku pamietnego...
+ double precision unmat(3,3) /1.0d0,0.0d0,0.0d0,
+ & 0.0d0,1.0d0,0.0d0,
+ & 0.0d0,0.0d0,1.0d0/
+c time00=MPI_Wtime()
+cd write (iout,*) "eelecij",i,j
+ ind=ind+1
+ iteli=itel(i)
+ itelj=itel(j)
+ if (j.eq.i+2 .and. itelj.eq.2) iteli=2
+ aaa=app(iteli,itelj)
+ bbb=bpp(iteli,itelj)
+ ael6i=ael6(iteli,itelj)
+ ael3i=ael3(iteli,itelj)
+ dxj=dc(1,j)
+ dyj=dc(2,j)
+ dzj=dc(3,j)
+ dx_normj=dc_norm(1,j)
+ dy_normj=dc_norm(2,j)
+ dz_normj=dc_norm(3,j)
+ xj=c(1,j)+0.5D0*dxj-xmedi
+ yj=c(2,j)+0.5D0*dyj-ymedi
+ zj=c(3,j)+0.5D0*dzj-zmedi
+ rij=xj*xj+yj*yj+zj*zj
+ rrmij=1.0D0/rij
+ rij=dsqrt(rij)
+ rmij=1.0D0/rij
+c For extracting the short-range part of Evdwpp
+ sss=sscale(rij/rpp(iteli,itelj))
+
+ r3ij=rrmij*rmij
+ r6ij=r3ij*r3ij
+ cosa=dx_normi*dx_normj+dy_normi*dy_normj+dz_normi*dz_normj
+ cosb=(xj*dx_normi+yj*dy_normi+zj*dz_normi)*rmij
+ cosg=(xj*dx_normj+yj*dy_normj+zj*dz_normj)*rmij
+ fac=cosa-3.0D0*cosb*cosg
+ ev1=aaa*r6ij*r6ij
+c 4/26/02 - AL scaling down 1,4 repulsive VDW interactions
+ if (j.eq.i+2) ev1=scal_el*ev1
+ ev2=bbb*r6ij
+ fac3=ael6i*r6ij
+ fac4=ael3i*r3ij
+ evdwij=ev1+ev2
+ el1=fac3*(4.0D0+fac*fac-3.0D0*(cosb*cosb+cosg*cosg))
+ el2=fac4*fac
+ eesij=el1+el2
+C 12/26/95 - for the evaluation of multi-body H-bonding interactions
+ ees0ij=4.0D0+fac*fac-3.0D0*(cosb*cosb+cosg*cosg)
+ ees=ees+eesij
+ evdw1=evdw1+evdwij*(1.0d0-sss)
+cd write(iout,'(2(2i3,2x),7(1pd12.4)/2(3(1pd12.4),5x)/)')
+cd & iteli,i,itelj,j,aaa,bbb,ael6i,ael3i,
+cd & 1.0D0/dsqrt(rrmij),evdwij,eesij,
+cd & xmedi,ymedi,zmedi,xj,yj,zj
+
+ if (energy_dec) then
+ write (iout,'(a6,2i5,0pf7.3,f7.3)') 'evdw1',i,j,evdwij,sss
+ write (iout,'(a6,2i5,0pf7.3)') 'ees',i,j,eesij
+ endif
+
+C
+C Calculate contributions to the Cartesian gradient.
+C
+#ifdef SPLITELE
+ facvdw=-6*rrmij*(ev1+evdwij)*(1.0d0-sss)
+ facel=-3*rrmij*(el1+eesij)
+ fac1=fac
+ erij(1)=xj*rmij
+ erij(2)=yj*rmij
+ erij(3)=zj*rmij
+*
+* Radial derivatives. First process both termini of the fragment (i,j)
+*
+ ggg(1)=facel*xj
+ ggg(2)=facel*yj
+ ggg(3)=facel*zj
+c do k=1,3
+c ghalf=0.5D0*ggg(k)
+c gelc(k,i)=gelc(k,i)+ghalf
+c gelc(k,j)=gelc(k,j)+ghalf
+c enddo
+c 9/28/08 AL Gradient compotents will be summed only at the end
+ do k=1,3
+ gelc_long(k,j)=gelc_long(k,j)+ggg(k)
+ gelc_long(k,i)=gelc_long(k,i)-ggg(k)
+ enddo
+*
+* Loop over residues i+1 thru j-1.
+*
+cgrad do k=i+1,j-1
+cgrad do l=1,3
+cgrad gelc(l,k)=gelc(l,k)+ggg(l)
+cgrad enddo
+cgrad enddo
+ ggg(1)=facvdw*xj
+ ggg(2)=facvdw*yj
+ ggg(3)=facvdw*zj
+c do k=1,3
+c ghalf=0.5D0*ggg(k)
+c gvdwpp(k,i)=gvdwpp(k,i)+ghalf
+c gvdwpp(k,j)=gvdwpp(k,j)+ghalf
+c enddo
+c 9/28/08 AL Gradient compotents will be summed only at the end
+ do k=1,3
+ gvdwpp(k,j)=gvdwpp(k,j)+ggg(k)
+ gvdwpp(k,i)=gvdwpp(k,i)-ggg(k)
+ enddo
+*
+* Loop over residues i+1 thru j-1.
+*
+cgrad do k=i+1,j-1
+cgrad do l=1,3
+cgrad gvdwpp(l,k)=gvdwpp(l,k)+ggg(l)
+cgrad enddo
+cgrad enddo
+#else
+ facvdw=ev1+evdwij*(1.0d0-sss)
+ facel=el1+eesij
+ fac1=fac
+ fac=-3*rrmij*(facvdw+facvdw+facel)
+ erij(1)=xj*rmij
+ erij(2)=yj*rmij
+ erij(3)=zj*rmij
+*
+* Radial derivatives. First process both termini of the fragment (i,j)
+*
+ ggg(1)=fac*xj
+ ggg(2)=fac*yj
+ ggg(3)=fac*zj
+c do k=1,3
+c ghalf=0.5D0*ggg(k)
+c gelc(k,i)=gelc(k,i)+ghalf
+c gelc(k,j)=gelc(k,j)+ghalf
+c enddo
+c 9/28/08 AL Gradient compotents will be summed only at the end
+ do k=1,3
+ gelc_long(k,j)=gelc(k,j)+ggg(k)
+ gelc_long(k,i)=gelc(k,i)-ggg(k)
+ enddo
+*
+* Loop over residues i+1 thru j-1.
+*
+cgrad do k=i+1,j-1
+cgrad do l=1,3
+cgrad gelc(l,k)=gelc(l,k)+ggg(l)
+cgrad enddo
+cgrad enddo
+c 9/28/08 AL Gradient compotents will be summed only at the end
+ ggg(1)=facvdw*xj
+ ggg(2)=facvdw*yj
+ ggg(3)=facvdw*zj
+ do k=1,3
+ gvdwpp(k,j)=gvdwpp(k,j)+ggg(k)
+ gvdwpp(k,i)=gvdwpp(k,i)-ggg(k)
+ enddo
+#endif
+*
+* Angular part
+*
+ ecosa=2.0D0*fac3*fac1+fac4
+ fac4=-3.0D0*fac4
+ fac3=-6.0D0*fac3
+ ecosb=(fac3*(fac1*cosg+cosb)+cosg*fac4)
+ ecosg=(fac3*(fac1*cosb+cosg)+cosb*fac4)
+ do k=1,3
+ dcosb(k)=rmij*(dc_norm(k,i)-erij(k)*cosb)
+ dcosg(k)=rmij*(dc_norm(k,j)-erij(k)*cosg)
+ enddo
+cd print '(2i3,2(3(1pd14.5),3x))',i,j,(dcosb(k),k=1,3),
+cd & (dcosg(k),k=1,3)
+ do k=1,3
+ ggg(k)=ecosb*dcosb(k)+ecosg*dcosg(k)
+ enddo
+c do k=1,3
+c ghalf=0.5D0*ggg(k)
+c gelc(k,i)=gelc(k,i)+ghalf
+c & +(ecosa*(dc_norm(k,j)-cosa*dc_norm(k,i))
+c & + ecosb*(erij(k)-cosb*dc_norm(k,i)))*vbld_inv(i+1)
+c gelc(k,j)=gelc(k,j)+ghalf
+c & +(ecosa*(dc_norm(k,i)-cosa*dc_norm(k,j))
+c & + ecosg*(erij(k)-cosg*dc_norm(k,j)))*vbld_inv(j+1)
+c enddo
+cgrad do k=i+1,j-1
+cgrad do l=1,3
+cgrad gelc(l,k)=gelc(l,k)+ggg(l)
+cgrad enddo
+cgrad enddo
+ do k=1,3
+ gelc(k,i)=gelc(k,i)
+ & +(ecosa*(dc_norm(k,j)-cosa*dc_norm(k,i))
+ & + ecosb*(erij(k)-cosb*dc_norm(k,i)))*vbld_inv(i+1)
+ gelc(k,j)=gelc(k,j)
+ & +(ecosa*(dc_norm(k,i)-cosa*dc_norm(k,j))
+ & + ecosg*(erij(k)-cosg*dc_norm(k,j)))*vbld_inv(j+1)
+ gelc_long(k,j)=gelc_long(k,j)+ggg(k)
+ gelc_long(k,i)=gelc_long(k,i)-ggg(k)
+ enddo
+ IF (wel_loc.gt.0.0d0 .or. wcorr4.gt.0.0d0 .or. wcorr5.gt.0.0d0
+ & .or. wcorr6.gt.0.0d0 .or. wturn3.gt.0.0d0
+ & .or. wturn4.gt.0.0d0 .or. wturn6.gt.0.0d0) THEN
+C
+C 9/25/99 Mixed third-order local-electrostatic terms. The local-interaction
+C energy of a peptide unit is assumed in the form of a second-order
+C Fourier series in the angles lambda1 and lambda2 (see Nishikawa et al.
+C Macromolecules, 1974, 7, 797-806 for definition). This correlation terms
+C are computed for EVERY pair of non-contiguous peptide groups.
+C
+ if (j.lt.nres-1) then
+ j1=j+1
+ j2=j-1
+ else
+ j1=j-1
+ j2=j-2
+ endif
+ kkk=0
+ do k=1,2
+ do l=1,2
+ kkk=kkk+1
+ muij(kkk)=mu(k,i)*mu(l,j)
+ enddo
+ enddo
+cd write (iout,*) 'EELEC: i',i,' j',j
+cd write (iout,*) 'j',j,' j1',j1,' j2',j2
+cd write(iout,*) 'muij',muij
+ ury=scalar(uy(1,i),erij)
+ urz=scalar(uz(1,i),erij)
+ vry=scalar(uy(1,j),erij)
+ vrz=scalar(uz(1,j),erij)
+ a22=scalar(uy(1,i),uy(1,j))-3*ury*vry
+ a23=scalar(uy(1,i),uz(1,j))-3*ury*vrz
+ a32=scalar(uz(1,i),uy(1,j))-3*urz*vry
+ a33=scalar(uz(1,i),uz(1,j))-3*urz*vrz
+ fac=dsqrt(-ael6i)*r3ij
+ a22=a22*fac
+ a23=a23*fac
+ a32=a32*fac
+ a33=a33*fac
+cd write (iout,'(4i5,4f10.5)')
+cd & i,itortyp(itype(i)),j,itortyp(itype(j)),a22,a23,a32,a33
+cd write (iout,'(6f10.5)') (muij(k),k=1,4),fac,eel_loc_ij
+cd write (iout,'(2(3f10.5,5x)/2(3f10.5,5x))') uy(:,i),uz(:,i),
+cd & uy(:,j),uz(:,j)
+cd write (iout,'(4f10.5)')
+cd & scalar(uy(1,i),uy(1,j)),scalar(uy(1,i),uz(1,j)),
+cd & scalar(uz(1,i),uy(1,j)),scalar(uz(1,i),uz(1,j))
+cd write (iout,'(4f10.5)') ury,urz,vry,vrz
+cd write (iout,'(9f10.5/)')
+cd & fac22,a22,fac23,a23,fac32,a32,fac33,a33,eel_loc_ij
+C Derivatives of the elements of A in virtual-bond vectors
+ call unormderiv(erij(1),unmat(1,1),rmij,erder(1,1))
+ do k=1,3
+ uryg(k,1)=scalar(erder(1,k),uy(1,i))
+ uryg(k,2)=scalar(uygrad(1,k,1,i),erij(1))
+ uryg(k,3)=scalar(uygrad(1,k,2,i),erij(1))
+ urzg(k,1)=scalar(erder(1,k),uz(1,i))
+ urzg(k,2)=scalar(uzgrad(1,k,1,i),erij(1))
+ urzg(k,3)=scalar(uzgrad(1,k,2,i),erij(1))
+ vryg(k,1)=scalar(erder(1,k),uy(1,j))
+ vryg(k,2)=scalar(uygrad(1,k,1,j),erij(1))
+ vryg(k,3)=scalar(uygrad(1,k,2,j),erij(1))
+ vrzg(k,1)=scalar(erder(1,k),uz(1,j))
+ vrzg(k,2)=scalar(uzgrad(1,k,1,j),erij(1))
+ vrzg(k,3)=scalar(uzgrad(1,k,2,j),erij(1))
+ enddo
+C Compute radial contributions to the gradient
+ facr=-3.0d0*rrmij
+ a22der=a22*facr
+ a23der=a23*facr
+ a32der=a32*facr
+ a33der=a33*facr
+ agg(1,1)=a22der*xj
+ agg(2,1)=a22der*yj
+ agg(3,1)=a22der*zj
+ agg(1,2)=a23der*xj
+ agg(2,2)=a23der*yj
+ agg(3,2)=a23der*zj
+ agg(1,3)=a32der*xj
+ agg(2,3)=a32der*yj
+ agg(3,3)=a32der*zj
+ agg(1,4)=a33der*xj
+ agg(2,4)=a33der*yj
+ agg(3,4)=a33der*zj
+C Add the contributions coming from er
+ fac3=-3.0d0*fac
+ do k=1,3
+ agg(k,1)=agg(k,1)+fac3*(uryg(k,1)*vry+vryg(k,1)*ury)
+ agg(k,2)=agg(k,2)+fac3*(uryg(k,1)*vrz+vrzg(k,1)*ury)
+ agg(k,3)=agg(k,3)+fac3*(urzg(k,1)*vry+vryg(k,1)*urz)
+ agg(k,4)=agg(k,4)+fac3*(urzg(k,1)*vrz+vrzg(k,1)*urz)
+ enddo
+ do k=1,3
+C Derivatives in DC(i)
+cgrad ghalf1=0.5d0*agg(k,1)
+cgrad ghalf2=0.5d0*agg(k,2)
+cgrad ghalf3=0.5d0*agg(k,3)
+cgrad ghalf4=0.5d0*agg(k,4)
+ aggi(k,1)=fac*(scalar(uygrad(1,k,1,i),uy(1,j))
+ & -3.0d0*uryg(k,2)*vry)!+ghalf1
+ aggi(k,2)=fac*(scalar(uygrad(1,k,1,i),uz(1,j))
+ & -3.0d0*uryg(k,2)*vrz)!+ghalf2
+ aggi(k,3)=fac*(scalar(uzgrad(1,k,1,i),uy(1,j))
+ & -3.0d0*urzg(k,2)*vry)!+ghalf3
+ aggi(k,4)=fac*(scalar(uzgrad(1,k,1,i),uz(1,j))
+ & -3.0d0*urzg(k,2)*vrz)!+ghalf4
+C Derivatives in DC(i+1)
+ aggi1(k,1)=fac*(scalar(uygrad(1,k,2,i),uy(1,j))
+ & -3.0d0*uryg(k,3)*vry)!+agg(k,1)
+ aggi1(k,2)=fac*(scalar(uygrad(1,k,2,i),uz(1,j))
+ & -3.0d0*uryg(k,3)*vrz)!+agg(k,2)
+ aggi1(k,3)=fac*(scalar(uzgrad(1,k,2,i),uy(1,j))
+ & -3.0d0*urzg(k,3)*vry)!+agg(k,3)
+ aggi1(k,4)=fac*(scalar(uzgrad(1,k,2,i),uz(1,j))
+ & -3.0d0*urzg(k,3)*vrz)!+agg(k,4)
+C Derivatives in DC(j)
+ aggj(k,1)=fac*(scalar(uygrad(1,k,1,j),uy(1,i))
+ & -3.0d0*vryg(k,2)*ury)!+ghalf1
+ aggj(k,2)=fac*(scalar(uzgrad(1,k,1,j),uy(1,i))
+ & -3.0d0*vrzg(k,2)*ury)!+ghalf2
+ aggj(k,3)=fac*(scalar(uygrad(1,k,1,j),uz(1,i))
+ & -3.0d0*vryg(k,2)*urz)!+ghalf3
+ aggj(k,4)=fac*(scalar(uzgrad(1,k,1,j),uz(1,i))
+ & -3.0d0*vrzg(k,2)*urz)!+ghalf4
+C Derivatives in DC(j+1) or DC(nres-1)
+ aggj1(k,1)=fac*(scalar(uygrad(1,k,2,j),uy(1,i))
+ & -3.0d0*vryg(k,3)*ury)
+ aggj1(k,2)=fac*(scalar(uzgrad(1,k,2,j),uy(1,i))
+ & -3.0d0*vrzg(k,3)*ury)
+ aggj1(k,3)=fac*(scalar(uygrad(1,k,2,j),uz(1,i))
+ & -3.0d0*vryg(k,3)*urz)
+ aggj1(k,4)=fac*(scalar(uzgrad(1,k,2,j),uz(1,i))
+ & -3.0d0*vrzg(k,3)*urz)
+cgrad if (j.eq.nres-1 .and. i.lt.j-2) then
+cgrad do l=1,4
+cgrad aggj1(k,l)=aggj1(k,l)+agg(k,l)
+cgrad enddo
+cgrad endif
+ enddo
+ acipa(1,1)=a22
+ acipa(1,2)=a23
+ acipa(2,1)=a32
+ acipa(2,2)=a33
+ a22=-a22
+ a23=-a23
+ do l=1,2
+ do k=1,3
+ agg(k,l)=-agg(k,l)
+ aggi(k,l)=-aggi(k,l)
+ aggi1(k,l)=-aggi1(k,l)
+ aggj(k,l)=-aggj(k,l)
+ aggj1(k,l)=-aggj1(k,l)
+ enddo
+ enddo
+ if (j.lt.nres-1) then
+ a22=-a22
+ a32=-a32
+ do l=1,3,2
+ do k=1,3
+ agg(k,l)=-agg(k,l)
+ aggi(k,l)=-aggi(k,l)
+ aggi1(k,l)=-aggi1(k,l)
+ aggj(k,l)=-aggj(k,l)
+ aggj1(k,l)=-aggj1(k,l)
+ enddo
+ enddo
+ else
+ a22=-a22
+ a23=-a23
+ a32=-a32
+ a33=-a33
+ do l=1,4
+ do k=1,3
+ agg(k,l)=-agg(k,l)
+ aggi(k,l)=-aggi(k,l)
+ aggi1(k,l)=-aggi1(k,l)
+ aggj(k,l)=-aggj(k,l)
+ aggj1(k,l)=-aggj1(k,l)
+ enddo
+ enddo
+ endif
+ ENDIF ! WCORR
+ IF (wel_loc.gt.0.0d0) THEN
+C Contribution to the local-electrostatic energy coming from the i-j pair
+ eel_loc_ij=a22*muij(1)+a23*muij(2)+a32*muij(3)
+ & +a33*muij(4)
+cd write (iout,*) 'i',i,' j',j,' eel_loc_ij',eel_loc_ij
+
+ if (energy_dec) write (iout,'(a6,2i5,0pf7.3)')
+ & 'eelloc',i,j,eel_loc_ij
+
+ eel_loc=eel_loc+eel_loc_ij
+C Partial derivatives in virtual-bond dihedral angles gamma
+ if (i.gt.1)
+ & gel_loc_loc(i-1)=gel_loc_loc(i-1)+
+ & a22*muder(1,i)*mu(1,j)+a23*muder(1,i)*mu(2,j)
+ & +a32*muder(2,i)*mu(1,j)+a33*muder(2,i)*mu(2,j)
+ gel_loc_loc(j-1)=gel_loc_loc(j-1)+
+ & a22*mu(1,i)*muder(1,j)+a23*mu(1,i)*muder(2,j)
+ & +a32*mu(2,i)*muder(1,j)+a33*mu(2,i)*muder(2,j)
+C Derivatives of eello in DC(i+1) thru DC(j-1) or DC(nres-2)
+ do l=1,3
+ ggg(l)=agg(l,1)*muij(1)+
+ & agg(l,2)*muij(2)+agg(l,3)*muij(3)+agg(l,4)*muij(4)
+ gel_loc_long(l,j)=gel_loc_long(l,j)+ggg(l)
+ gel_loc_long(l,i)=gel_loc_long(l,i)-ggg(l)
+cgrad ghalf=0.5d0*ggg(l)
+cgrad gel_loc(l,i)=gel_loc(l,i)+ghalf
+cgrad gel_loc(l,j)=gel_loc(l,j)+ghalf
+ enddo
+cgrad do k=i+1,j2
+cgrad do l=1,3
+cgrad gel_loc(l,k)=gel_loc(l,k)+ggg(l)
+cgrad enddo
+cgrad enddo
+C Remaining derivatives of eello
+ do l=1,3
+ gel_loc(l,i)=gel_loc(l,i)+aggi(l,1)*muij(1)+
+ & aggi(l,2)*muij(2)+aggi(l,3)*muij(3)+aggi(l,4)*muij(4)
+ gel_loc(l,i+1)=gel_loc(l,i+1)+aggi1(l,1)*muij(1)+
+ & aggi1(l,2)*muij(2)+aggi1(l,3)*muij(3)+aggi1(l,4)*muij(4)
+ gel_loc(l,j)=gel_loc(l,j)+aggj(l,1)*muij(1)+
+ & aggj(l,2)*muij(2)+aggj(l,3)*muij(3)+aggj(l,4)*muij(4)
+ gel_loc(l,j1)=gel_loc(l,j1)+aggj1(l,1)*muij(1)+
+ & aggj1(l,2)*muij(2)+aggj1(l,3)*muij(3)+aggj1(l,4)*muij(4)
+ enddo
+ ENDIF
+C Change 12/26/95 to calculate four-body contributions to H-bonding energy
+c if (j.gt.i+1 .and. num_conti.le.maxconts) then
+ if (wcorr+wcorr4+wcorr5+wcorr6.gt.0.0d0
+ & .and. num_conti.le.maxconts) then
+c write (iout,*) i,j," entered corr"
+C
+C Calculate the contact function. The ith column of the array JCONT will
+C contain the numbers of atoms that make contacts with the atom I (of numbers
+C greater than I). The arrays FACONT and GACONT will contain the values of
+C the contact function and its derivative.
+c r0ij=1.02D0*rpp(iteli,itelj)
+c r0ij=1.11D0*rpp(iteli,itelj)
+ r0ij=2.20D0*rpp(iteli,itelj)
+c r0ij=1.55D0*rpp(iteli,itelj)
+ call gcont(rij,r0ij,1.0D0,0.2d0*r0ij,fcont,fprimcont)
+ if (fcont.gt.0.0D0) then
+ num_conti=num_conti+1
+ if (num_conti.gt.maxconts) then
+ write (iout,*) 'WARNING - max. # of contacts exceeded;',
+ & ' will skip next contacts for this conf.'
+ else
+ jcont_hb(num_conti,i)=j
+cd write (iout,*) "i",i," j",j," num_conti",num_conti,
+cd & " jcont_hb",jcont_hb(num_conti,i)
+ IF (wcorr4.gt.0.0d0 .or. wcorr5.gt.0.0d0 .or.
+ & wcorr6.gt.0.0d0 .or. wturn6.gt.0.0d0) THEN
+C 9/30/99 (AL) - store components necessary to evaluate higher-order loc-el
+C terms.
+ d_cont(num_conti,i)=rij
+cd write (2,'(3e15.5)') rij,r0ij+0.2d0*r0ij,rij
+C --- Electrostatic-interaction matrix ---
+ a_chuj(1,1,num_conti,i)=a22
+ a_chuj(1,2,num_conti,i)=a23
+ a_chuj(2,1,num_conti,i)=a32
+ a_chuj(2,2,num_conti,i)=a33
+C --- Gradient of rij
+ do kkk=1,3
+ grij_hb_cont(kkk,num_conti,i)=erij(kkk)
+ enddo
+ kkll=0
+ do k=1,2
+ do l=1,2
+ kkll=kkll+1
+ do m=1,3
+ a_chuj_der(k,l,m,1,num_conti,i)=agg(m,kkll)
+ a_chuj_der(k,l,m,2,num_conti,i)=aggi(m,kkll)
+ a_chuj_der(k,l,m,3,num_conti,i)=aggi1(m,kkll)
+ a_chuj_der(k,l,m,4,num_conti,i)=aggj(m,kkll)
+ a_chuj_der(k,l,m,5,num_conti,i)=aggj1(m,kkll)
+ enddo
+ enddo
+ enddo
+ ENDIF
+ IF (wcorr4.eq.0.0d0 .and. wcorr.gt.0.0d0) THEN
+C Calculate contact energies
+ cosa4=4.0D0*cosa
+ wij=cosa-3.0D0*cosb*cosg
+ cosbg1=cosb+cosg
+ cosbg2=cosb-cosg
+c fac3=dsqrt(-ael6i)/r0ij**3
+ fac3=dsqrt(-ael6i)*r3ij
+c ees0pij=dsqrt(4.0D0+cosa4+wij*wij-3.0D0*cosbg1*cosbg1)
+ ees0tmp=4.0D0+cosa4+wij*wij-3.0D0*cosbg1*cosbg1
+ if (ees0tmp.gt.0) then
+ ees0pij=dsqrt(ees0tmp)
+ else
+ ees0pij=0
+ endif
+c ees0mij=dsqrt(4.0D0-cosa4+wij*wij-3.0D0*cosbg2*cosbg2)
+ ees0tmp=4.0D0-cosa4+wij*wij-3.0D0*cosbg2*cosbg2
+ if (ees0tmp.gt.0) then
+ ees0mij=dsqrt(ees0tmp)
+ else
+ ees0mij=0
+ endif
+c ees0mij=0.0D0
+ ees0p(num_conti,i)=0.5D0*fac3*(ees0pij+ees0mij)
+ ees0m(num_conti,i)=0.5D0*fac3*(ees0pij-ees0mij)
+C Diagnostics. Comment out or remove after debugging!
+c ees0p(num_conti,i)=0.5D0*fac3*ees0pij
+c ees0m(num_conti,i)=0.5D0*fac3*ees0mij
+c ees0m(num_conti,i)=0.0D0
+C End diagnostics.
+c write (iout,*) 'i=',i,' j=',j,' rij=',rij,' r0ij=',r0ij,
+c & ' ees0ij=',ees0p(num_conti,i),ees0m(num_conti,i),' fcont=',fcont
+C Angular derivatives of the contact function
+ ees0pij1=fac3/ees0pij
+ ees0mij1=fac3/ees0mij
+ fac3p=-3.0D0*fac3*rrmij
+ ees0pijp=0.5D0*fac3p*(ees0pij+ees0mij)
+ ees0mijp=0.5D0*fac3p*(ees0pij-ees0mij)
+c ees0mij1=0.0D0
+ ecosa1= ees0pij1*( 1.0D0+0.5D0*wij)
+ ecosb1=-1.5D0*ees0pij1*(wij*cosg+cosbg1)
+ ecosg1=-1.5D0*ees0pij1*(wij*cosb+cosbg1)
+ ecosa2= ees0mij1*(-1.0D0+0.5D0*wij)
+ ecosb2=-1.5D0*ees0mij1*(wij*cosg+cosbg2)
+ ecosg2=-1.5D0*ees0mij1*(wij*cosb-cosbg2)
+ ecosap=ecosa1+ecosa2
+ ecosbp=ecosb1+ecosb2
+ ecosgp=ecosg1+ecosg2
+ ecosam=ecosa1-ecosa2
+ ecosbm=ecosb1-ecosb2
+ ecosgm=ecosg1-ecosg2
+C Diagnostics
+c ecosap=ecosa1
+c ecosbp=ecosb1
+c ecosgp=ecosg1
+c ecosam=0.0D0
+c ecosbm=0.0D0
+c ecosgm=0.0D0
+C End diagnostics
+ facont_hb(num_conti,i)=fcont
+ fprimcont=fprimcont/rij
+cd facont_hb(num_conti,i)=1.0D0
+C Following line is for diagnostics.
+cd fprimcont=0.0D0
+ do k=1,3
+ dcosb(k)=rmij*(dc_norm(k,i)-erij(k)*cosb)
+ dcosg(k)=rmij*(dc_norm(k,j)-erij(k)*cosg)
+ enddo
+ do k=1,3
+ gggp(k)=ecosbp*dcosb(k)+ecosgp*dcosg(k)
+ gggm(k)=ecosbm*dcosb(k)+ecosgm*dcosg(k)
+ enddo
+ gggp(1)=gggp(1)+ees0pijp*xj
+ gggp(2)=gggp(2)+ees0pijp*yj
+ gggp(3)=gggp(3)+ees0pijp*zj
+ gggm(1)=gggm(1)+ees0mijp*xj
+ gggm(2)=gggm(2)+ees0mijp*yj
+ gggm(3)=gggm(3)+ees0mijp*zj
+C Derivatives due to the contact function
+ gacont_hbr(1,num_conti,i)=fprimcont*xj
+ gacont_hbr(2,num_conti,i)=fprimcont*yj
+ gacont_hbr(3,num_conti,i)=fprimcont*zj
+ do k=1,3
+c
+c 10/24/08 cgrad and ! comments indicate the parts of the code removed
+c following the change of gradient-summation algorithm.
+c
+cgrad ghalfp=0.5D0*gggp(k)
+cgrad ghalfm=0.5D0*gggm(k)
+ gacontp_hb1(k,num_conti,i)=!ghalfp
+ & +(ecosap*(dc_norm(k,j)-cosa*dc_norm(k,i))
+ & + ecosbp*(erij(k)-cosb*dc_norm(k,i)))*vbld_inv(i+1)
+ gacontp_hb2(k,num_conti,i)=!ghalfp
+ & +(ecosap*(dc_norm(k,i)-cosa*dc_norm(k,j))
+ & + ecosgp*(erij(k)-cosg*dc_norm(k,j)))*vbld_inv(j+1)
+ gacontp_hb3(k,num_conti,i)=gggp(k)
+ gacontm_hb1(k,num_conti,i)=!ghalfm
+ & +(ecosam*(dc_norm(k,j)-cosa*dc_norm(k,i))
+ & + ecosbm*(erij(k)-cosb*dc_norm(k,i)))*vbld_inv(i+1)
+ gacontm_hb2(k,num_conti,i)=!ghalfm
+ & +(ecosam*(dc_norm(k,i)-cosa*dc_norm(k,j))
+ & + ecosgm*(erij(k)-cosg*dc_norm(k,j)))*vbld_inv(j+1)
+ gacontm_hb3(k,num_conti,i)=gggm(k)
+ enddo
+ ENDIF ! wcorr
+ endif ! num_conti.le.maxconts
+ endif ! fcont.gt.0
+ endif ! j.gt.i+1
+ if (wturn3.gt.0.0d0 .or. wturn4.gt.0.0d0) then
+ do k=1,4
+ do l=1,3
+ ghalf=0.5d0*agg(l,k)
+ aggi(l,k)=aggi(l,k)+ghalf
+ aggi1(l,k)=aggi1(l,k)+agg(l,k)
+ aggj(l,k)=aggj(l,k)+ghalf
+ enddo
+ enddo
+ if (j.eq.nres-1 .and. i.lt.j-2) then
+ do k=1,4
+ do l=1,3
+ aggj1(l,k)=aggj1(l,k)+agg(l,k)
+ enddo
+ enddo
+ endif
+ endif
+c t_eelecij=t_eelecij+MPI_Wtime()-time00
+ return
+ end
+C-----------------------------------------------------------------------
+ subroutine evdwpp_short(evdw1)
+C
+C Compute Evdwpp
+C
+ implicit real*8 (a-h,o-z)
+ include 'DIMENSIONS'
+ include 'COMMON.CONTROL'
+ include 'COMMON.IOUNITS'
+ include 'COMMON.GEO'
+ include 'COMMON.VAR'
+ include 'COMMON.LOCAL'
+ include 'COMMON.CHAIN'
+ include 'COMMON.DERIV'
+ include 'COMMON.INTERACT'
+ include 'COMMON.CONTACTS'
+ include 'COMMON.TORSION'
+ include 'COMMON.VECTORS'
+ include 'COMMON.FFIELD'
+ dimension ggg(3)
+c 4/26/02 - AL scaling factor for 1,4 repulsive VDW interactions
+#ifdef MOMENT
+ double precision scal_el /1.0d0/
+#else
+ double precision scal_el /0.5d0/
+#endif
+ evdw1=0.0D0
+c write (iout,*) "iatel_s_vdw",iatel_s_vdw,
+c & " iatel_e_vdw",iatel_e_vdw
+ call flush(iout)
+ do i=iatel_s_vdw,iatel_e_vdw
+ if (itype(i).eq.ntyp1.or. itype(i+1).eq.ntyp1) cycle
+ dxi=dc(1,i)
+ dyi=dc(2,i)
+ dzi=dc(3,i)
+ dx_normi=dc_norm(1,i)
+ dy_normi=dc_norm(2,i)
+ dz_normi=dc_norm(3,i)
+ xmedi=c(1,i)+0.5d0*dxi
+ ymedi=c(2,i)+0.5d0*dyi
+ zmedi=c(3,i)+0.5d0*dzi
+ num_conti=0
+c write (iout,*) 'i',i,' ielstart',ielstart_vdw(i),
+c & ' ielend',ielend_vdw(i)
+ call flush(iout)
+ do j=ielstart_vdw(i),ielend_vdw(i)
+ if (itype(j).eq.ntyp1 .or. itype(j+1).eq.ntyp1) cycle
+ ind=ind+1
+ iteli=itel(i)
+ itelj=itel(j)
+ if (j.eq.i+2 .and. itelj.eq.2) iteli=2
+ aaa=app(iteli,itelj)
+ bbb=bpp(iteli,itelj)
+ dxj=dc(1,j)
+ dyj=dc(2,j)
+ dzj=dc(3,j)
+ dx_normj=dc_norm(1,j)
+ dy_normj=dc_norm(2,j)
+ dz_normj=dc_norm(3,j)
+ xj=c(1,j)+0.5D0*dxj-xmedi
+ yj=c(2,j)+0.5D0*dyj-ymedi
+ zj=c(3,j)+0.5D0*dzj-zmedi
+ rij=xj*xj+yj*yj+zj*zj
+ rrmij=1.0D0/rij
+ rij=dsqrt(rij)
+ sss=sscale(rij/rpp(iteli,itelj))
+ if (sss.gt.0.0d0) then
+ rmij=1.0D0/rij
+ r3ij=rrmij*rmij
+ r6ij=r3ij*r3ij
+ ev1=aaa*r6ij*r6ij
+c 4/26/02 - AL scaling down 1,4 repulsive VDW interactions
+ if (j.eq.i+2) ev1=scal_el*ev1
+ ev2=bbb*r6ij
+ evdwij=ev1+ev2
+ if (energy_dec) then
+ write (iout,'(a6,2i5,0pf7.3,f7.3)') 'evdw1',i,j,evdwij,sss
+ endif
+ evdw1=evdw1+evdwij*sss
+C
+C Calculate contributions to the Cartesian gradient.
+C
+ facvdw=-6*rrmij*(ev1+evdwij)*sss
+ ggg(1)=facvdw*xj
+ ggg(2)=facvdw*yj
+ ggg(3)=facvdw*zj
+ do k=1,3
+ gvdwpp(k,j)=gvdwpp(k,j)+ggg(k)
+ gvdwpp(k,i)=gvdwpp(k,i)-ggg(k)
+ enddo
+ endif
+ enddo ! j
+ enddo ! i
+ return
+ end
+C-----------------------------------------------------------------------------
+ subroutine escp_long(evdw2,evdw2_14)
+C
+C This subroutine calculates the excluded-volume interaction energy between
+C peptide-group centers and side chains and its gradient in virtual-bond and
+C side-chain vectors.
+C
+ implicit real*8 (a-h,o-z)
+ include 'DIMENSIONS'
+ include 'COMMON.GEO'
+ include 'COMMON.VAR'
+ include 'COMMON.LOCAL'
+ include 'COMMON.CHAIN'
+ include 'COMMON.DERIV'
+ include 'COMMON.INTERACT'
+ include 'COMMON.FFIELD'
+ include 'COMMON.IOUNITS'
+ include 'COMMON.CONTROL'
+ dimension ggg(3)
+ evdw2=0.0D0
+ evdw2_14=0.0d0
+cd print '(a)','Enter ESCP'
+cd write (iout,*) 'iatscp_s=',iatscp_s,' iatscp_e=',iatscp_e
+ do i=iatscp_s,iatscp_e
+ if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1) cycle
+ iteli=itel(i)
+ xi=0.5D0*(c(1,i)+c(1,i+1))
+ yi=0.5D0*(c(2,i)+c(2,i+1))
+ zi=0.5D0*(c(3,i)+c(3,i+1))
+
+ do iint=1,nscp_gr(i)
+
+ do j=iscpstart(i,iint),iscpend(i,iint)
+ itypj=itype(j)
+ if (itypj.eq.ntyp1) cycle
+C Uncomment following three lines for SC-p interactions
+c xj=c(1,nres+j)-xi
+c yj=c(2,nres+j)-yi
+c zj=c(3,nres+j)-zi
+C Uncomment following three lines for Ca-p interactions
+ xj=c(1,j)-xi
+ yj=c(2,j)-yi
+ zj=c(3,j)-zi
+ rrij=1.0D0/(xj*xj+yj*yj+zj*zj)
+
+ sss=sscale(1.0d0/(dsqrt(rrij)*rscp(itypj,iteli)))
+
+ if (sss.lt.1.0d0) then
+
+ fac=rrij**expon2
+ e1=fac*fac*aad(itypj,iteli)
+ e2=fac*bad(itypj,iteli)
+ if (iabs(j-i) .le. 2) then
+ e1=scal14*e1
+ e2=scal14*e2
+ evdw2_14=evdw2_14+(e1+e2)*(1.0d0-sss)
+ endif
+ evdwij=e1+e2
+ evdw2=evdw2+evdwij*(1.0d0-sss)
+ if (energy_dec) write (iout,'(a6,2i5,2(0pf7.3))')
+ & 'evdw2',i,j,sss,evdwij
+C
+C Calculate contributions to the gradient in the virtual-bond and SC vectors.
+C
+ fac=-(evdwij+e1)*rrij*(1.0d0-sss)
+ ggg(1)=xj*fac
+ ggg(2)=yj*fac
+ ggg(3)=zj*fac
+C Uncomment following three lines for SC-p interactions
+c do k=1,3
+c gradx_scp(k,j)=gradx_scp(k,j)+ggg(k)
+c enddo
+C Uncomment following line for SC-p interactions
+c gradx_scp(k,j)=gradx_scp(k,j)+ggg(k)
+ do k=1,3
+ gvdwc_scpp(k,i)=gvdwc_scpp(k,i)-ggg(k)
+ gvdwc_scp(k,j)=gvdwc_scp(k,j)+ggg(k)
+ enddo
+ endif
+ enddo
+
+ enddo ! iint
+ enddo ! i
+ do i=1,nct
+ do j=1,3
+ gvdwc_scp(j,i)=expon*gvdwc_scp(j,i)
+ gvdwc_scpp(j,i)=expon*gvdwc_scpp(j,i)
+ gradx_scp(j,i)=expon*gradx_scp(j,i)
+ enddo
+ enddo
+C******************************************************************************
+C
+C N O T E !!!
+C
+C To save time the factor EXPON has been extracted from ALL components
+C of GVDWC and GRADX. Remember to multiply them by this factor before further
+C use!
+C
+C******************************************************************************
+ return
+ end
+C-----------------------------------------------------------------------------
+ subroutine escp_short(evdw2,evdw2_14)
+C
+C This subroutine calculates the excluded-volume interaction energy between
+C peptide-group centers and side chains and its gradient in virtual-bond and
+C side-chain vectors.
+C
+ implicit real*8 (a-h,o-z)
+ include 'DIMENSIONS'
+ include 'COMMON.GEO'
+ include 'COMMON.VAR'
+ include 'COMMON.LOCAL'
+ include 'COMMON.CHAIN'
+ include 'COMMON.DERIV'
+ include 'COMMON.INTERACT'
+ include 'COMMON.FFIELD'
+ include 'COMMON.IOUNITS'
+ include 'COMMON.CONTROL'
+ dimension ggg(3)
+ evdw2=0.0D0
+ evdw2_14=0.0d0
+cd print '(a)','Enter ESCP'
+cd write (iout,*) 'iatscp_s=',iatscp_s,' iatscp_e=',iatscp_e
+ do i=iatscp_s,iatscp_e
+ if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1) cycle
+ iteli=itel(i)
+ xi=0.5D0*(c(1,i)+c(1,i+1))
+ yi=0.5D0*(c(2,i)+c(2,i+1))
+ zi=0.5D0*(c(3,i)+c(3,i+1))
+
+ do iint=1,nscp_gr(i)
+
+ do j=iscpstart(i,iint),iscpend(i,iint)
+ itypj=itype(j)
+ if (itypj.eq.ntyp1) cycle
+C Uncomment following three lines for SC-p interactions
+c xj=c(1,nres+j)-xi
+c yj=c(2,nres+j)-yi
+c zj=c(3,nres+j)-zi
+C Uncomment following three lines for Ca-p interactions
+ xj=c(1,j)-xi
+ yj=c(2,j)-yi
+ zj=c(3,j)-zi
+ rrij=1.0D0/(xj*xj+yj*yj+zj*zj)
+
+ sss=sscale(1.0d0/(dsqrt(rrij)*rscp(itypj,iteli)))
+
+ if (sss.gt.0.0d0) then
+
+ fac=rrij**expon2
+ e1=fac*fac*aad(itypj,iteli)
+ e2=fac*bad(itypj,iteli)
+ if (iabs(j-i) .le. 2) then
+ e1=scal14*e1
+ e2=scal14*e2
+ evdw2_14=evdw2_14+(e1+e2)*sss
+ endif
+ evdwij=e1+e2
+ evdw2=evdw2+evdwij*sss
+ if (energy_dec) write (iout,'(a6,2i5,2(0pf7.3))')
+ & 'evdw2',i,j,sss,evdwij
+C
+C Calculate contributions to the gradient in the virtual-bond and SC vectors.
+C
+ fac=-(evdwij+e1)*rrij*sss
+ ggg(1)=xj*fac
+ ggg(2)=yj*fac
+ ggg(3)=zj*fac
+C Uncomment following three lines for SC-p interactions
+c do k=1,3
+c gradx_scp(k,j)=gradx_scp(k,j)+ggg(k)
+c enddo
+C Uncomment following line for SC-p interactions
+c gradx_scp(k,j)=gradx_scp(k,j)+ggg(k)
+ do k=1,3
+ gvdwc_scpp(k,i)=gvdwc_scpp(k,i)-ggg(k)
+ gvdwc_scp(k,j)=gvdwc_scp(k,j)+ggg(k)
+ enddo
+ endif
+ enddo
+
+ enddo ! iint
+ enddo ! i
+ do i=1,nct
+ do j=1,3
+ gvdwc_scp(j,i)=expon*gvdwc_scp(j,i)
+ gvdwc_scpp(j,i)=expon*gvdwc_scpp(j,i)
+ gradx_scp(j,i)=expon*gradx_scp(j,i)
+ enddo
+ enddo
+C******************************************************************************
+C
+C N O T E !!!
+C
+C To save time the factor EXPON has been extracted from ALL components
+C of GVDWC and GRADX. Remember to multiply them by this factor before further
+C use!
+C
+C******************************************************************************
+ return
+ end
projects / unres.git / blobdiff