+++ /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)
- 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)
- 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
- 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)
- enddo
- do k=i,j-1
- do l=1,3
- gvdwc(l,k)=gvdwc(l,k)+gg(l)
- enddo
- 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)
- 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)
- 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.gt.0.0d0) then
- rrij=1.0D0/rij
- 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)
- enddo
- do k=i,j-1
- do l=1,3
- gvdwc(l,k)=gvdwc(l,k)+gg(l)
- enddo
- 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)
- 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)
- 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+evdwij*(1.0d0-sss)
-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)
- enddo
- do k=i,j-1
- do l=1,3
- gvdwc(l,k)=gvdwc(l,k)+gg(l)
- enddo
- 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)
- 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)
- 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+evdwij*sss
-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)
- enddo
- do k=i,j-1
- do l=1,3
- gvdwc(l,k)=gvdwc(l,k)+gg(l)
- enddo
- 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)
- 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)
-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)
-C For diagnostics only!!!
-c chi1=0.0D0
-c chi2=0.0D0
-c chi12=0.0D0
-c chip1=0.0D0
-c chip2=0.0D0
-c chip12=0.0D0
-c alf1=0.0D0
-c alf2=0.0D0
-c alf12=0.0D0
- 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)
-cd if (icall.eq.0) then
-cd rrsave(ind)=rrij
-cd else
-cd rrij=rrsave(ind)
-cd endif
- 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)
- 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)
-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)
-C For diagnostics only!!!
-c chi1=0.0D0
-c chi2=0.0D0
-c chi12=0.0D0
-c chip1=0.0D0
-c chip2=0.0D0
-c chip12=0.0D0
-c alf1=0.0D0
-c alf2=0.0D0
-c alf12=0.0D0
- 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)
-cd if (icall.eq.0) then
-cd rrsave(ind)=rrij
-cd else
-cd rrij=rrsave(ind)
-cd endif
- 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)
- 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)
-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)
-C For diagnostics only!!!
-c chi1=0.0D0
-c chi2=0.0D0
-c chi12=0.0D0
-c chip1=0.0D0
-c chip2=0.0D0
-c chip12=0.0D0
-c alf1=0.0D0
-c alf2=0.0D0
-c alf12=0.0D0
- 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)
-c write (iout,*) "dcnorj",dxi*dxi+dyi*dyi+dzi*dzi
-c write (iout,*) "j",j," dc_norm",
-c & dc_norm(1,nres+j),dc_norm(2,nres+j),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)))
- write(iout,*) "long",i,itypi,j,itypj," rij",1.0d0/rij,
- & " sigmaii",sigmaii(itypi,itypj)," sss",sss
-
- 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)
- write (iout,*) "evdwij",evdwij," evdw",evdw
- 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
-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)
- 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)
-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)
-C For diagnostics only!!!
-c chi1=0.0D0
-c chi2=0.0D0
-c chi12=0.0D0
-c chip1=0.0D0
-c chip2=0.0D0
-c chip12=0.0D0
-c alf1=0.0D0
-c alf2=0.0D0
-c alf12=0.0D0
- 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)
-c write (iout,*) "dcnorj",dxi*dxi+dyi*dyi+dzi*dzi
-c write (iout,*) "j",j," dc_norm",
-c & dc_norm(1,nres+j),dc_norm(2,nres+j),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)))
- write(iout,*) "short",i,itypi,j,itypj," rij",1.0d0/rij,
- & " sigmaii",sigmaii(itypi,itypj)," sss",sss
- 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
- write (iout,*) "evdwij",evdwij," evdw",evdw
- 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
-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)
- 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)
-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)
-C For diagnostics only!!!
-c chi1=0.0D0
-c chi2=0.0D0
-c chi12=0.0D0
-c chip1=0.0D0
-c chip2=0.0D0
-c chip12=0.0D0
-c alf1=0.0D0
-c alf2=0.0D0
-c alf12=0.0D0
- 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)
- 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)
-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)
-C For diagnostics only!!!
-c chi1=0.0D0
-c chi2=0.0D0
-c chi12=0.0D0
-c chip1=0.0D0
-c chip2=0.0D0
-c chip12=0.0D0
-c alf1=0.0D0
-c alf2=0.0D0
-c alf12=0.0D0
- 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 k=i,j-1
- do l=1,3
- gvdwc(l,k)=gvdwc(l,k)+gg(l)
- enddo
- 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)
- 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),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,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
- call set_matrices
- 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
- num_conti_hb=0
- 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
- do i=iatel_s,iatel_e
- 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(i),' ielend',ielend(i)
- do j=ielstart(i),ielend(i)
- 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)
-C Diagnostics only!!!
-c aaa=0.0D0
-c bbb=0.0D0
-c ael6i=0.0D0
-c ael3i=0.0D0
-C End diagnostics
- 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))
-c
- 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)') 'evdw1',i,j,evdwij
- 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
- do k=1,3
- ghalf=0.5D0*ggg(k)
- gelc(k,i)=gelc(k,i)+ghalf
- gelc(k,j)=gelc(k,j)+ghalf
- enddo
-*
-* Loop over residues i+1 thru j-1.
-*
- do k=i+1,j-1
- do l=1,3
- gelc(l,k)=gelc(l,k)+ggg(l)
- enddo
- enddo
- ggg(1)=facvdw*xj
- ggg(2)=facvdw*yj
- ggg(3)=facvdw*zj
- do k=1,3
- ghalf=0.5D0*ggg(k)
- gvdwpp(k,i)=gvdwpp(k,i)+ghalf
- gvdwpp(k,j)=gvdwpp(k,j)+ghalf
- enddo
-*
-* Loop over residues i+1 thru j-1.
-*
- do k=i+1,j-1
- do l=1,3
- gvdwpp(l,k)=gvdwpp(l,k)+ggg(l)
- enddo
- 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
- do k=1,3
- ghalf=0.5D0*ggg(k)
- gelc(k,i)=gelc(k,i)+ghalf
- gelc(k,j)=gelc(k,j)+ghalf
- enddo
-*
-* Loop over residues i+1 thru j-1.
-*
- do k=i+1,j-1
- do l=1,3
- gelc(l,k)=gelc(l,k)+ggg(l)
- enddo
- 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
- do k=1,3
- ghalf=0.5D0*ggg(k)
- gelc(k,i)=gelc(k,i)+ghalf
- & +(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)+ghalf
- & +(ecosa*(dc_norm(k,i)-cosa*dc_norm(k,j))
- & + ecosg*(erij(k)-cosg*dc_norm(k,j)))*vbld_inv(j+1)
- enddo
- do k=i+1,j-1
- do l=1,3
- gelc(l,k)=gelc(l,k)+ggg(l)
- enddo
- 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
-C For diagnostics only
-cd a22=1.0d0
-cd a23=1.0d0
-cd a32=1.0d0
-cd a33=1.0d0
- fac=dsqrt(-ael6i)*r3ij
-cd write (2,*) 'fac=',fac
-C For diagnostics only
-cd fac=1.0d0
- 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))
-cd do k=1,3
-cd do l=1,3
-cd erder(k,l)=0.0d0
-cd enddo
-cd enddo
- 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
-cd do k=1,3
-cd do l=1,3
-cd uryg(k,l)=0.0d0
-cd urzg(k,l)=0.0d0
-cd vryg(k,l)=0.0d0
-cd vrzg(k,l)=0.0d0
-cd enddo
-cd enddo
-C Compute radial contributions to the gradient
- facr=-3.0d0*rrmij
- a22der=a22*facr
- a23der=a23*facr
- a32der=a32*facr
- a33der=a33*facr
-cd a22der=0.0d0
-cd a23der=0.0d0
-cd a32der=0.0d0
-cd a33der=0.0d0
- 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)
- ghalf1=0.5d0*agg(k,1)
- ghalf2=0.5d0*agg(k,2)
- ghalf3=0.5d0*agg(k,3)
- 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)
-cd aggi(k,1)=ghalf1
-cd aggi(k,2)=ghalf2
-cd aggi(k,3)=ghalf3
-cd aggi(k,4)=ghalf4
-C Derivatives in DC(i+1)
-cd aggi1(k,1)=agg(k,1)
-cd aggi1(k,2)=agg(k,2)
-cd aggi1(k,3)=agg(k,3)
-cd aggi1(k,4)=agg(k,4)
-C Derivatives in DC(j)
-cd aggj(k,1)=ghalf1
-cd aggj(k,2)=ghalf2
-cd aggj(k,3)=ghalf3
-cd aggj(k,4)=ghalf4
-C Derivatives in DC(j+1)
-cd aggj1(k,1)=0.0d0
-cd aggj1(k,2)=0.0d0
-cd aggj1(k,3)=0.0d0
-cd aggj1(k,4)=0.0d0
- if (j.eq.nres-1 .and. i.lt.j-2) then
- do l=1,4
- aggj1(k,l)=aggj1(k,l)+agg(k,l)
-cd aggj1(k,l)=agg(k,l)
- enddo
- endif
- enddo
-c goto 11111
-C Check the loc-el terms by numerical integration
- 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
-11111 continue
- 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)
-cd call checkint3(i,j,mu1,mu2,a22,a23,a32,a33,acipa,eel_loc_ij)
-cd write(iout,*) 'agg ',agg
-cd write(iout,*) 'aggi ',aggi
-cd write(iout,*) 'aggi1',aggi1
-cd write(iout,*) 'aggj ',aggj
-cd write(iout,*) 'aggj1',aggj1
-
-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)
- enddo
- do k=i+2,j2
- do l=1,3
- gel_loc(l,k)=gel_loc(l,k)+ggg(l)
- enddo
- 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
- if (wturn3.gt.0.0d0 .or. wturn4.gt.0.0d0) then
-C Contributions from turns
- a_temp(1,1)=a22
- a_temp(1,2)=a23
- a_temp(2,1)=a32
- a_temp(2,2)=a33
- call eturn34(i,j,eello_turn3,eello_turn4)
- endif
-C Change 12/26/95 to calculate four-body contributions to H-bonding energy
- if (j.gt.i+1 .and. num_conti.le.maxconts) then
-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
- 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
-c if (i.eq.1) then
-c a_chuj(1,1,num_conti,i)=-0.61d0
-c a_chuj(1,2,num_conti,i)= 0.4d0
-c a_chuj(2,1,num_conti,i)= 0.65d0
-c a_chuj(2,2,num_conti,i)= 0.50d0
-c else if (i.eq.2) then
-c a_chuj(1,1,num_conti,i)= 0.0d0
-c a_chuj(1,2,num_conti,i)= 0.0d0
-c a_chuj(2,1,num_conti,i)= 0.0d0
-c a_chuj(2,2,num_conti,i)= 0.0d0
-c endif
-C --- and its gradients
-cd write (iout,*) 'i',i,' j',j
-cd do kkk=1,3
-cd write (iout,*) 'iii 1 kkk',kkk
-cd write (iout,*) agg(kkk,:)
-cd enddo
-cd do kkk=1,3
-cd write (iout,*) 'iii 2 kkk',kkk
-cd write (iout,*) aggi(kkk,:)
-cd enddo
-cd do kkk=1,3
-cd write (iout,*) 'iii 3 kkk',kkk
-cd write (iout,*) aggi1(kkk,:)
-cd enddo
-cd do kkk=1,3
-cd write (iout,*) 'iii 4 kkk',kkk
-cd write (iout,*) aggj(kkk,:)
-cd enddo
-cd do kkk=1,3
-cd write (iout,*) 'iii 5 kkk',kkk
-cd write (iout,*) aggj1(kkk,:)
-cd 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)
-c do mm=1,5
-c a_chuj_der(k,l,m,mm,num_conti,i)=0.0d0
-c enddo
- 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
- ghalfp=0.5D0*gggp(k)
- 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
-C Diagnostics. Comment out or remove after debugging!
-cdiag do k=1,3
-cdiag gacontp_hb1(k,num_conti,i)=0.0D0
-cdiag gacontp_hb2(k,num_conti,i)=0.0D0
-cdiag gacontp_hb3(k,num_conti,i)=0.0D0
-cdiag gacontm_hb1(k,num_conti,i)=0.0D0
-cdiag gacontm_hb2(k,num_conti,i)=0.0D0
-cdiag gacontm_hb3(k,num_conti,i)=0.0D0
-cdiag enddo
- ENDIF ! wcorr
- endif ! num_conti.le.maxconts
- endif ! fcont.gt.0
- endif ! j.gt.i+1
- enddo ! j
- num_cont_hb(i)=num_conti
- enddo ! i
-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
- return
- end
-C-----------------------------------------------------------------------
- subroutine evdwpp_long(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
- do i=iatel_s,iatel_e
- 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(i),' ielend',ielend(i)
- do j=ielstart(i),ielend(i)
- 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.lt.1.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)') 'evdw1',i,j,evdwij
- write (iout,'(a6,2i5,0pf7.3)') 'ees',i,j,eesij
- endif
- evdw1=evdw1+evdwij*(1.0d0-sss)
-C
-C Calculate contributions to the Cartesian gradient.
-C
- facvdw=-6*rrmij*(ev1+evdwij)*(1.0d0-sss)
- ggg(1)=facvdw*xj
- ggg(2)=facvdw*yj
- ggg(3)=facvdw*zj
-
- do k=1,3
- ghalf=0.5D0*ggg(k)
- gvdwpp(k,i)=gvdwpp(k,i)+ghalf
- gvdwpp(k,j)=gvdwpp(k,j)+ghalf
- enddo
-*
-* Loop over residues i+1 thru j-1.
-*
- do k=i+1,j-1
- do l=1,3
- gvdwpp(l,k)=gvdwpp(l,k)+ggg(l)
- enddo
- enddo
- endif
- enddo ! j
- enddo ! i
- 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
- do i=iatel_s,iatel_e
- 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(i),' ielend',ielend(i)
- do j=ielstart(i),ielend(i)
- 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)') 'evdw1',i,j,evdwij
- write (iout,'(a6,2i5,0pf7.3)') 'ees',i,j,eesij
- 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
- ghalf=0.5D0*ggg(k)
- gvdwpp(k,i)=gvdwpp(k,i)+ghalf
- gvdwpp(k,j)=gvdwpp(k,j)+ghalf
- enddo
-*
-* Loop over residues i+1 thru j-1.
-*
- do k=i+1,j-1
- do l=1,3
- gvdwpp(l,k)=gvdwpp(l,k)+ggg(l)
- enddo
- 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
- 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)
-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,0pf7.3)')
- & 'evdw2',i,j,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
- if (j.lt.i) then
-cd write (iout,*) 'j<i'
-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
- else
-cd write (iout,*) 'j>i'
- do k=1,3
- ggg(k)=-ggg(k)
-C Uncomment following line for SC-p interactions
-c gradx_scp(k,j)=gradx_scp(k,j)-ggg(k)
- enddo
- endif
- do k=1,3
- gvdwc_scp(k,i)=gvdwc_scp(k,i)-0.5D0*ggg(k)
- enddo
- kstart=min0(i+1,j)
- kend=max0(i-1,j-1)
-cd write (iout,*) 'i=',i,' j=',j,' kstart=',kstart,' kend=',kend
-cd write (iout,*) ggg(1),ggg(2),ggg(3)
- do k=kstart,kend
- do l=1,3
- gvdwc_scp(l,k)=gvdwc_scp(l,k)-ggg(l)
- enddo
- enddo
-
- endif
-
- enddo
-
- enddo ! iint
- enddo ! i
- do i=1,nct
- do j=1,3
- gvdwc_scp(j,i)=expon*gvdwc_scp(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
- 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)
-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,0pf7.3)')
- & 'evdw2',i,j,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
- if (j.lt.i) then
-cd write (iout,*) 'j<i'
-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
- else
-cd write (iout,*) 'j>i'
- do k=1,3
- ggg(k)=-ggg(k)
-C Uncomment following line for SC-p interactions
-c gradx_scp(k,j)=gradx_scp(k,j)-ggg(k)
- enddo
- endif
- do k=1,3
- gvdwc_scp(k,i)=gvdwc_scp(k,i)-0.5D0*ggg(k)
- enddo
- kstart=min0(i+1,j)
- kend=max0(i-1,j-1)
-cd write (iout,*) 'i=',i,' j=',j,' kstart=',kstart,' kend=',kend
-cd write (iout,*) ggg(1),ggg(2),ggg(3)
- do k=kstart,kend
- do l=1,3
- gvdwc_scp(l,k)=gvdwc_scp(l,k)-ggg(l)
- enddo
- enddo
-
- endif
-
- enddo
-
- enddo ! iint
- enddo ! i
- do i=1,nct
- do j=1,3
- gvdwc_scp(j,i)=expon*gvdwc_scp(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