X-Git-Url: http://mmka.chem.univ.gda.pl/gitweb/?a=blobdiff_plain;f=source%2Funres%2Fsrc_MD-M-newcorr%2Fenergy_p_new-sep.F;fp=source%2Funres%2Fsrc_MD-M-newcorr%2Fenergy_p_new-sep.F;h=0b8f27b07646e94854c1853316a020b691decf2c;hb=d101c97dea752458d76055fdbae49c26fff03c1f;hp=0000000000000000000000000000000000000000;hpb=325eda160c9ad2982501e091ca40606a29043712;p=unres.git diff --git a/source/unres/src_MD-M-newcorr/energy_p_new-sep.F b/source/unres/src_MD-M-newcorr/energy_p_new-sep.F new file mode 100644 index 0000000..0b8f27b --- /dev/null +++ b/source/unres/src_MD-M-newcorr/energy_p_new-sep.F @@ -0,0 +1,2505 @@ +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))) +c write(iout,*) "long",i,itypi,j,itypj," rij",1.0d0/rij, +c & " 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) +c 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))) +c write(iout,*) "short",i,itypi,j,itypj," rij",1.0d0/rij, +c & " 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 +c 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 +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 +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 + 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 + 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+3,ees,evdw1,eel_loc) + if (wturn4.gt.0.0d0) call eturn4(i,eello_turn4) + num_cont_hb(i)=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 + 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) + call eelecij_scale(i,j,ees,evdw1,eel_loc) + enddo ! j + num_cont_hb(i)=num_cont_hb(i)+num_conti + enddo ! i + return + end +C------------------------------------------------------------------------------- + subroutine eelecij_scale(i,j,ees,evdw1,eel_loc) +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,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/ + 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 +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 + 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,*) 'ji' + 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,*) 'ji' + 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