C----------------------------------------------------------------------- double precision function sscalelip(r) implicit none double precision r,gamm include "COMMON.SPLITELE" C if(r.lt.r_cut-rlamb) then C sscale=1.0d0 C else if(r.le.r_cut.and.r.ge.r_cut-rlamb) then C gamm=(r-(r_cut-rlamb))/rlamb sscalelip=1.0d0+r*r*(2*r-3.0d0) C else C sscale=0d0 C endif return end C----------------------------------------------------------------------- double precision function sscagradlip(r) implicit none double precision r,gamm include "COMMON.SPLITELE" C if(r.lt.r_cut-rlamb) then C sscagrad=0.0d0 C else if(r.le.r_cut.and.r.ge.r_cut-rlamb) then C gamm=(r-(r_cut-rlamb))/rlamb sscagradlip=r*(6*r-6.0d0) C else C sscagrad=0.0d0 C endif return end C----------------------------------------------------------------------- double precision function sscale(r,r_cut) implicit none double precision r,r_cut,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----------------------------------------------------------------------- double precision function sscagrad(r,r_cut) implicit none double precision r,r_cut,gamm include "COMMON.SPLITELE" if(r.lt.r_cut-rlamb) then sscagrad=0.0d0 else if(r.le.r_cut.and.r.ge.r_cut-rlamb) then gamm=(r-(r_cut-rlamb))/rlamb sscagrad=gamm*(6*gamm-6.0d0)/rlamb else sscagrad=0.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 none include 'DIMENSIONS' 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.SPLITELE" c include 'COMMON.CONTACTS' double precision gg(3) double precision evdw,evdwij integer i,j,k,itypi,itypj,itypi1,num_conti,iint double precision xi,yi,zi,xj,yj,zj,rij,eps0ij,fac,e1,e2,rrij, & sigij,r0ij,rcut,sss1,sssgrad1,sqrij double precision sscale,sscagrad c write(iout,*)'Entering ELJ nnt=',nnt,' nct=',nct,' expon=',expon evdw=0.0D0 do i=iatsc_s,iatsc_e itypi=iabs(itype(i)) if (itypi.eq.ntyp1) cycle itypi1=iabs(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=iabs(itype(j)) if (itypj.eq.ntyp1) cycle xj=c(1,nres+j)-xi yj=c(2,nres+j)-yi zj=c(3,nres+j)-zi rij=xj*xj+yj*yj+zj*zj sqrij=dsqrt(rrij) eps0ij=eps(itypi,itypj) sss1=sscale(sqrij,r_cut_int) if (sss1.eq.0.0d0) cycle sssgrad1=sscagrad(sqrij,r_cut_int) sssgrad= & sscagrad(sqrij/sigma(itypi,itypj),r_cut_respa) sss=sscale(sqrij/sigma(itypi,itypj),r_cut_respa) if (sss.lt.1.0d0) then rrij=1.0D0/rij fac=rrij**expon2 e1=fac*fac*aa e2=fac*bb evdwij=e1+e2 evdw=evdw+(1.0d0-sss)*sss1*evdwij/sqrij/expon C C Calculate the components of the gradient in DC and X C fac=-rrij*(e1+evdwij)*(1.0d0-sss)*sss1 & +evdwij*(-sss1*sssgrad/sigma(itypi,itypj) & +(1.0d0-sss)*sssgrad1)/sqrij gg(1)=xj*fac gg(2)=yj*fac gg(3)=zj*fac do k=1,3 gvdwx(k,i)=gvdwx(k,i)-gg(k) gvdwx(k,j)=gvdwx(k,j)+gg(k) gvdwc(k,i)=gvdwc(k,i)-gg(k) gvdwc(k,j)=gvdwc(k,j)+gg(k) enddo endif enddo ! j enddo ! iint enddo ! i do i=1,nct do j=1,3 gvdwc(j,i)=expon*gvdwc(j,i) gvdwx(j,i)=expon*gvdwx(j,i) enddo enddo C****************************************************************************** C C N O T E !!! C C To save time, the factor of EXPON has been extracted from ALL components C of GVDWC and GRADX. Remember to multiply them by this factor before further C use! C C****************************************************************************** return end C----------------------------------------------------------------------- subroutine elj_short(evdw) C C This subroutine calculates the interaction energy of nonbonded side chains C assuming the LJ potential of interaction. C implicit none include 'DIMENSIONS' 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.SPLITELE" c include 'COMMON.CONTACTS' double precision gg(3) double precision evdw,evdwij integer i,j,k,itypi,itypj,itypi1,num_conti,iint double precision xi,yi,zi,xj,yj,zj,rij,eps0ij,fac,e1,e2,rrij, & sigij,r0ij,rcut,sqrij,sss1,sssgrad1 double precision sscale,sscagrad c write(iout,*)'Entering ELJ nnt=',nnt,' nct=',nct,' expon=',expon evdw=0.0D0 do i=iatsc_s,iatsc_e itypi=iabs(itype(i)) if (itypi.eq.ntyp1) cycle itypi1=iabs(itype(i+1)) xi=c(1,nres+i) yi=c(2,nres+i) zi=c(3,nres+i) C Change 12/1/95 num_conti=0 C C Calculate SC interaction energy. C do iint=1,nint_gr(i) cd write (iout,*) 'i=',i,' iint=',iint,' istart=',istart(i,iint), cd & 'iend=',iend(i,iint) do j=istart(i,iint),iend(i,iint) itypj=iabs(itype(j)) if (itypj.eq.ntyp1) cycle xj=c(1,nres+j)-xi yj=c(2,nres+j)-yi zj=c(3,nres+j)-zi C Change 12/1/95 to calculate four-body interactions rij=xj*xj+yj*yj+zj*zj sqrij=dsqrt(rij) sss=sscale(sqrij/sigma(itypi,itypj),r_cut_respa) if (sss.gt.0.0d0) then sssgrad= & sscagrad(sqrij/sigma(itypi,itypj),r_cut_respa) rrij=1.0D0/rij eps0ij=eps(itypi,itypj) fac=rrij**expon2 e1=fac*fac*aa e2=fac*bb 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+evdwij*sssgrad/sqrij/expon gg(1)=xj*fac gg(2)=yj*fac gg(3)=zj*fac do k=1,3 gvdwx(k,i)=gvdwx(k,i)-gg(k) gvdwx(k,j)=gvdwx(k,j)+gg(k) gvdwc(k,i)=gvdwc(k,i)-gg(k) gvdwc(k,j)=gvdwc(k,j)+gg(k) enddo endif enddo ! j enddo ! iint enddo ! i do i=1,nct do j=1,3 gvdwc(j,i)=expon*gvdwc(j,i) gvdwx(j,i)=expon*gvdwx(j,i) enddo enddo C****************************************************************************** C C N O T E !!! C C To save time, the factor of EXPON has been extracted from ALL components C of GVDWC and GRADX. Remember to multiply them by this factor before further C use! C C****************************************************************************** return end C----------------------------------------------------------------------------- subroutine eljk_long(evdw) C C This subroutine calculates the interaction energy of nonbonded side chains C assuming the LJK potential of interaction. C implicit none 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' include "COMMON.SPLITELE" double precision gg(3) double precision evdw,evdwij integer i,j,k,itypi,itypj,itypi1,iint double precision xi,yi,zi,xj,yj,zj,rij,eps0ij,fac,e1,e2,rrij, & fac_augm,e_augm,r_inv_ij,r_shift_inv,sss1,sssgrad1 logical scheck double precision sscale,sscagrad c print *,'Entering ELJK nnt=',nnt,' nct=',nct,' expon=',expon evdw=0.0D0 do i=iatsc_s,iatsc_e itypi=iabs(itype(i)) if (itypi.eq.ntyp1) cycle itypi1=iabs(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=iabs(itype(j)) if (itypj.eq.ntyp1) cycle xj=c(1,nres+j)-xi yj=c(2,nres+j)-yi zj=c(3,nres+j)-zi rrij=1.0D0/(xj*xj+yj*yj+zj*zj) fac_augm=rrij**expon e_augm=augm(itypi,itypj)*fac_augm r_inv_ij=dsqrt(rrij) rij=1.0D0/r_inv_ij sss1=sscale(rij,r_cut_int) if (sss1.eq.0.0d0) cycle sssgrad1=sscagrad(rij,r_cut_int) sss=sscale(rij/sigma(itypi,itypj),r_cut_respa) if (sss.lt.1.0d0) then sssgrad= & sscagrad(rij/sigma(itypi,itypj),r_cut_respa) r_shift_inv=1.0D0/(rij+r0(itypi,itypj)-sigma(itypi,itypj)) fac=r_shift_inv**expon e1=fac*fac*aa e2=fac*bb evdwij=e_augm+e1+e2 cd sigm=dabs(aa(itypi,itypj)/bb(itypi,itypj))**(1.0D0/6.0D0) cd epsi=bb(itypi,itypj)**2/aa(itypi,itypj) cd write (iout,'(2(a3,i3,2x),8(1pd12.4)/2(3(1pd12.4),5x)/)') cd & restyp(itypi),i,restyp(itypj),j,aa(itypi,itypj), cd & bb(itypi,itypj),augm(itypi,itypj),epsi,sigm, cd & sigma(itypi,itypj),1.0D0/dsqrt(rrij),evdwij, cd & (c(k,i),k=1,3),(c(k,j),k=1,3) evdw=evdw+(1.0d0-sss)*sss1*evdwij C C Calculate the components of the gradient in DC and X C fac=-2.0D0*rrij*e_augm-r_inv_ij*r_shift_inv*(e1+e1+e2) fac=fac*(1.0d0-sss)*sss1 & +evdwij*(-sss1*sssgrad/sigma(itypi,itypj) & +(1.0d0-sss)*sssgrad1)*r_inv_ij/expon gg(1)=xj*fac gg(2)=yj*fac gg(3)=zj*fac do k=1,3 gvdwx(k,i)=gvdwx(k,i)-gg(k) gvdwx(k,j)=gvdwx(k,j)+gg(k) gvdwc(k,i)=gvdwc(k,i)-gg(k) gvdwc(k,j)=gvdwc(k,j)+gg(k) enddo endif enddo ! j enddo ! iint enddo ! i do i=1,nct do j=1,3 gvdwc(j,i)=expon*gvdwc(j,i) gvdwx(j,i)=expon*gvdwx(j,i) enddo enddo return end C----------------------------------------------------------------------------- subroutine eljk_short(evdw) C C This subroutine calculates the interaction energy of nonbonded side chains C assuming the LJK potential of interaction. C implicit real*8 (a-h,o-z) include 'DIMENSIONS' include 'COMMON.GEO' include 'COMMON.VAR' include 'COMMON.LOCAL' include 'COMMON.CHAIN' include 'COMMON.DERIV' include 'COMMON.INTERACT' include 'COMMON.IOUNITS' include 'COMMON.NAMES' include "COMMON.SPLITELE" double precision gg(3) double precision evdw,evdwij integer i,j,k,itypi,itypj,itypi1,iint double precision xi,yi,zi,xj,yj,zj,rij,eps0ij,fac,e1,e2,rrij, & fac_augm,e_augm,r_inv_ij,r_shift_inv,sss1,sssgrad1 logical scheck double precision sscale,sscagrad c print *,'Entering ELJK nnt=',nnt,' nct=',nct,' expon=',expon evdw=0.0D0 do i=iatsc_s,iatsc_e itypi=iabs(itype(i)) if (itypi.eq.ntyp1) cycle itypi1=iabs(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=iabs(itype(j)) if (itypj.eq.ntyp1) cycle xj=c(1,nres+j)-xi yj=c(2,nres+j)-yi zj=c(3,nres+j)-zi rrij=1.0D0/(xj*xj+yj*yj+zj*zj) fac_augm=rrij**expon e_augm=augm(itypi,itypj)*fac_augm r_inv_ij=dsqrt(rrij) rij=1.0D0/r_inv_ij sss=sscale(rij/sigma(itypi,itypj),r_cut_respa) 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 e2=fac*bb evdwij=e_augm+e1+e2 cd sigm=dabs(aa(itypi,itypj)/bb(itypi,itypj))**(1.0D0/6.0D0) cd epsi=bb(itypi,itypj)**2/aa(itypi,itypj) cd write (iout,'(2(a3,i3,2x),8(1pd12.4)/2(3(1pd12.4),5x)/)') cd & restyp(itypi),i,restyp(itypj),j,aa(itypi,itypj), cd & bb(itypi,itypj),augm(itypi,itypj),epsi,sigm, cd & sigma(itypi,itypj),1.0D0/dsqrt(rrij),evdwij, cd & (c(k,i),k=1,3),(c(k,j),k=1,3) evdw=evdw+sss*evdwij C C Calculate the components of the gradient in DC and X C fac=-2.0D0*rrij*e_augm-r_inv_ij*r_shift_inv*(e1+e1+e2) & +evdwij*sssgrad/sigma(itypi,itypj)*r_inv_ij/expon fac=fac*sss gg(1)=xj*fac gg(2)=yj*fac gg(3)=zj*fac do k=1,3 gvdwx(k,i)=gvdwx(k,i)-gg(k) gvdwx(k,j)=gvdwx(k,j)+gg(k) gvdwc(k,i)=gvdwc(k,i)-gg(k) gvdwc(k,j)=gvdwc(k,j)+gg(k) enddo endif enddo ! j enddo ! iint enddo ! i do i=1,nct do j=1,3 gvdwc(j,i)=expon*gvdwc(j,i) gvdwx(j,i)=expon*gvdwx(j,i) enddo enddo return end C----------------------------------------------------------------------------- subroutine ebp_long(evdw) C C This subroutine calculates the interaction energy of nonbonded side chains C assuming the Berne-Pechukas potential of interaction. C implicit none 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.SPLITELE" integer icall common /srutu/ icall double precision evdw integer itypi,itypj,itypi1,iint,ind double precision eps0ij,epsi,sigm,fac,e1,e2,rrij,xi,yi,zi double precision sss1,sssgrad1 double precision sscale,sscagrad 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=iabs(itype(i)) if (itypi.eq.ntyp1) cycle itypi1=iabs(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=iabs(itype(j)) if (itypj.eq.ntyp1) cycle c dscj_inv=dsc_inv(itypj) dscj_inv=vbld_inv(j+nres) chi1=chi(itypi,itypj) chi2=chi(itypj,itypi) chi12=chi1*chi2 chip1=chip(itypi) chip2=chip(itypj) chip12=chip1*chip2 alf1=alp(itypi) alf2=alp(itypj) alf12=0.5D0*(alf1+alf2) xj=c(1,nres+j)-xi yj=c(2,nres+j)-yi zj=c(3,nres+j)-zi dxj=dc_norm(1,nres+j) dyj=dc_norm(2,nres+j) dzj=dc_norm(3,nres+j) rrij=1.0D0/(xj*xj+yj*yj+zj*zj) rij=dsqrt(rrij) sss1=sscale(1.0d0/rij,r_cut_int) if (sss1.eq.0.0d0) cycle sss=sscale(1.0d0/(rij*sigmaii(itypi,itypj)),r_cut_respa) if (sss.lt.1.0d0) then sssgrad= & sscagrad(1.0d0/(rij*sigmaii(itypi,itypj)),r_cut_respa) sssgrad1=sscagrad(1.0d0/rij,r_cut_int) 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 e2=fac*bb evdwij=eps1*eps2rt*eps3rt*(e1+e2) eps2der=evdwij*eps3rt eps3der=evdwij*eps2rt evdwij=evdwij*eps2rt*eps3rt evdw=evdw+evdwij*(1.0d0-sss)*sss1 if (lprn) then sigm=dabs(aa/bb)**(1.0D0/6.0D0) epsi=bb**2/aa 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=(fac+evdwij*(sss1/(1.0d0-sss)*sssgrad/ & sigmaii(itypi,itypj)+(1.0d0-sss)/sss1*sssgrad1))*rij 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)*sss1) 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' include "COMMON.SPLITELE" integer icall common /srutu/ icall double precision evdw integer itypi,itypj,itypi1,iint,ind double precision eps0ij,epsi,sigm,fac,e1,e2,rrij,xi,yi,zi double precision sscale,sscagrad 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=iabs(itype(i)) if (itypi.eq.ntyp1) cycle itypi1=iabs(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=iabs(itype(j)) if (itypj.eq.ntyp1) cycle c dscj_inv=dsc_inv(itypj) dscj_inv=vbld_inv(j+nres) chi1=chi(itypi,itypj) chi2=chi(itypj,itypi) chi12=chi1*chi2 chip1=chip(itypi) chip2=chip(itypj) chip12=chip1*chip2 alf1=alp(itypi) alf2=alp(itypj) alf12=0.5D0*(alf1+alf2) xj=c(1,nres+j)-xi yj=c(2,nres+j)-yi zj=c(3,nres+j)-zi dxj=dc_norm(1,nres+j) dyj=dc_norm(2,nres+j) dzj=dc_norm(3,nres+j) rrij=1.0D0/(xj*xj+yj*yj+zj*zj) rij=dsqrt(rrij) sss=sscale(1.0d0/(rij*sigmaii(itypi,itypj)),r_cut_respa) 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 e2=fac*bb 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/bb)**(1.0D0/6.0D0) epsi=bb**2/aa 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=(fac+evdwij*sssgrad/sss/sigmaii(itypi,itypj))*rrij 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 none 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' include "COMMON.SPLITELE" logical lprn integer xshift,yshift,zshift double precision evdw integer itypi,itypj,itypi1,iint,ind double precision eps0ij,epsi,sigm,fac,e1,e2,rrij,xi,yi,zi double precision fracinbuf,sslipi,evdwij_przed_tri,sig0ij, & sslipj,ssgradlipj,ssgradlipi,dist_init,xj_safe,yj_safe,zj_safe, & xj_temp,yj_temp,zj_temp,dist_temp,sig,rij_shift,faclip double precision dist,sscale,sscagrad,sscagradlip,sscalelip double precision subchap,sss1,sssgrad1 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=iabs(itype(i)) if (itypi.eq.ntyp1) cycle itypi1=iabs(itype(i+1)) xi=c(1,nres+i) yi=c(2,nres+i) zi=c(3,nres+i) xi=mod(xi,boxxsize) if (xi.lt.0) xi=xi+boxxsize yi=mod(yi,boxysize) if (yi.lt.0) yi=yi+boxysize zi=mod(zi,boxzsize) if (zi.lt.0) zi=zi+boxzsize 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=iabs(itype(j)) if (itypj.eq.ntyp1) cycle c dscj_inv=dsc_inv(itypj) dscj_inv=vbld_inv(j+nres) c write (iout,*) "j",j,dsc_inv(itypj),dscj_inv, c & 1.0d0/vbld(j+nres) c write (iout,*) "i",i," j", j," itype",itype(i),itype(j) sig0ij=sigma(itypi,itypj) chi1=chi(itypi,itypj) chi2=chi(itypj,itypi) chi12=chi1*chi2 chip1=chip(itypi) chip2=chip(itypj) chip12=chip1*chip2 alf1=alp(itypi) alf2=alp(itypj) alf12=0.5D0*(alf1+alf2) xj=c(1,nres+j) yj=c(2,nres+j) zj=c(3,nres+j) xj=mod(xj,boxxsize) if (xj.lt.0) xj=xj+boxxsize yj=mod(yj,boxysize) if (yj.lt.0) yj=yj+boxysize zj=mod(zj,boxzsize) if (zj.lt.0) zj=zj+boxzsize if ((zj.gt.bordlipbot).and.(zj.lt.bordliptop)) then C the energy transfer exist if (zj.lt.buflipbot) then C what fraction I am in fracinbuf=1.0d0-((zi-bordlipbot)/lipbufthick) C lipbufthick is thickenes of lipid buffore sslipj=sscalelip(fracinbuf) ssgradlipj=-sscagradlip(fracinbuf)/lipbufthick else if (zi.gt.bufliptop) then fracinbuf=1.0d0-((bordliptop-zi)/lipbufthick) sslipj=sscalelip(fracinbuf) ssgradlipj=sscagradlip(fracinbuf)/lipbufthick else sslipj=1.0d0 ssgradlipj=0.0 endif else sslipj=0.0d0 ssgradlipj=0.0 endif aa=aa_lip(itypi,itypj)*(sslipi+sslipj)/2.0d0 & +aa_aq(itypi,itypj)*(2.0d0-sslipi+sslipj)/2.0d0 bb=bb_lip(itypi,itypj)*(sslipi+sslipj)/2.0d0 & +bb_aq(itypi,itypj)*(2.0d0-sslipi+sslipj)/2.0d0 dist_init=(xj-xi)**2+(yj-yi)**2+(zj-zi)**2 xj_safe=xj yj_safe=yj zj_safe=zj subchap=0 do xshift=-1,1 do yshift=-1,1 do zshift=-1,1 xj=xj_safe+xshift*boxxsize yj=yj_safe+yshift*boxysize zj=zj_safe+zshift*boxzsize dist_temp=(xj-xi)**2+(yj-yi)**2+(zj-zi)**2 if (dist_temp.lt.dist_init) then dist_init=dist_temp xj_temp=xj yj_temp=yj zj_temp=zj subchap=1 endif enddo enddo enddo if (subchap.eq.1) then xj=xj_temp-xi yj=yj_temp-yi zj=zj_temp-zi else xj=xj_safe-xi yj=yj_safe-yi zj=zj_safe-zi endif 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) sss1=sscale(1.0d0/rij,r_cut_int) if (sss1.eq.0.0d0) cycle sss=sscale(1.0d0/(rij*sigmaii(itypi,itypj)),r_cut_respa) if (sss.lt.1.0d0) then C Calculate angle-dependent terms of energy and contributions to their C derivatives. sssgrad= & sscagrad((1.0d0/rij)/sigmaii(itypi,itypj),r_cut_respa) sssgrad1=sscagrad(1.0d0/rij,r_cut_int) 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 e2=fac*bb 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)*sss1 if (lprn) then sigm=dabs(aa/bb)**(1.0D0/6.0D0) epsi=bb**2/aa 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,4f10.5)') & 'evdw',i,j,rij,sss,sss1,evdwij C Calculate gradient components. e1=e1*eps1*eps2rt**2*eps3rt**2 fac=-expon*(e1+evdwij)*rij_shift sigder=fac*sigder fac=(fac+evdwij*(-sss1*sssgrad/(1.0d0-sss) & /sigmaii(itypi,itypj)+(1.0d0-sss)*sssgrad1/sss1))*rij c fac=0.0d0 C Calculate the radial part of the gradient gg(1)=xj*fac gg(2)=yj*fac gg(3)=zj*fac gg_lipi(3)=ssgradlipi*evdwij gg_lipj(3)=ssgradlipj*evdwij C Calculate angular part of the gradient. call sc_grad_scale((1.0d0-sss)*sss1) endif enddo ! j enddo ! iint enddo ! i c write (iout,*) "Number of loop steps in EGB:",ind cccc energy_dec=.false. return end C----------------------------------------------------------------------------- subroutine egb_short(evdw) C C This subroutine calculates the interaction energy of nonbonded side chains C assuming the Gay-Berne potential of interaction. C implicit none 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' include "COMMON.SPLITELE" logical lprn integer xshift,yshift,zshift double precision evdw integer itypi,itypj,itypi1,iint,ind double precision eps0ij,epsi,sigm,fac,e1,e2,rrij,xi,yi,zi double precision fracinbuf,sslipi,evdwij_przed_tri,sig0ij, & sslipj,ssgradlipj,ssgradlipi,dist_init,xj_safe,yj_safe,zj_safe, & xj_temp,yj_temp,zj_temp,dist_temp,sig,rij_shift,faclip double precision dist,sscale,sscagrad,sscagradlip,sscalelip double precision subchap 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=iabs(itype(i)) if (itypi.eq.ntyp1) cycle itypi1=iabs(itype(i+1)) xi=c(1,nres+i) yi=c(2,nres+i) zi=c(3,nres+i) xi=mod(xi,boxxsize) if (xi.lt.0) xi=xi+boxxsize yi=mod(yi,boxysize) if (yi.lt.0) yi=yi+boxysize zi=mod(zi,boxzsize) if (zi.lt.0) zi=zi+boxzsize 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=iabs(itype(j)) if (itypj.eq.ntyp1) cycle c dscj_inv=dsc_inv(itypj) dscj_inv=vbld_inv(j+nres) c write (iout,*) "j",j,dsc_inv(itypj),dscj_inv, c & 1.0d0/vbld(j+nres) c write (iout,*) "i",i," j", j," itype",itype(i),itype(j) sig0ij=sigma(itypi,itypj) chi1=chi(itypi,itypj) chi2=chi(itypj,itypi) chi12=chi1*chi2 chip1=chip(itypi) chip2=chip(itypj) chip12=chip1*chip2 alf1=alp(itypi) alf2=alp(itypj) alf12=0.5D0*(alf1+alf2) xj=c(1,nres+j) yj=c(2,nres+j) zj=c(3,nres+j) xj=mod(xj,boxxsize) if (xj.lt.0) xj=xj+boxxsize yj=mod(yj,boxysize) if (yj.lt.0) yj=yj+boxysize zj=mod(zj,boxzsize) if (zj.lt.0) zj=zj+boxzsize if ((zj.gt.bordlipbot).and.(zj.lt.bordliptop)) then C the energy transfer exist if (zj.lt.buflipbot) then C what fraction I am in fracinbuf=1.0d0-((zi-bordlipbot)/lipbufthick) C lipbufthick is thickenes of lipid buffore sslipj=sscalelip(fracinbuf) ssgradlipj=-sscagradlip(fracinbuf)/lipbufthick elseif (zi.gt.bufliptop) then fracinbuf=1.0d0-((bordliptop-zi)/lipbufthick) sslipj=sscalelip(fracinbuf) ssgradlipj=sscagradlip(fracinbuf)/lipbufthick else sslipj=1.0d0 ssgradlipj=0.0 endif else sslipj=0.0d0 ssgradlipj=0.0 endif aa=aa_lip(itypi,itypj)*(sslipi+sslipj)/2.0d0 & +aa_aq(itypi,itypj)*(2.0d0-sslipi+sslipj)/2.0d0 bb=bb_lip(itypi,itypj)*(sslipi+sslipj)/2.0d0 & +bb_aq(itypi,itypj)*(2.0d0-sslipi+sslipj)/2.0d0 dist_init=(xj-xi)**2+(yj-yi)**2+(zj-zi)**2 xj_safe=xj yj_safe=yj zj_safe=zj subchap=0 do xshift=-1,1 do yshift=-1,1 do zshift=-1,1 xj=xj_safe+xshift*boxxsize yj=yj_safe+yshift*boxysize zj=zj_safe+zshift*boxzsize dist_temp=(xj-xi)**2+(yj-yi)**2+(zj-zi)**2 if(dist_temp.lt.dist_init) then dist_init=dist_temp xj_temp=xj yj_temp=yj zj_temp=zj subchap=1 endif enddo enddo enddo if (subchap.eq.1) then xj=xj_temp-xi yj=yj_temp-yi zj=zj_temp-zi else xj=xj_safe-xi yj=yj_safe-yi zj=zj_safe-zi endif 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)),r_cut_respa) sssgrad=sscagrad((1.0d0/rij)/sigmaii(itypi,itypj),r_cut_respa) 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 e2=fac*bb evdwij=eps1*eps2rt*eps3rt*(e1+e2) eps2der=evdwij*eps3rt eps3der=evdwij*eps2rt c write (iout,*) "sigsq",sigsq," sig",sig," eps2rt",eps2rt, c & " eps3rt",eps3rt," eps1",eps1," e1",e1," e2",e2 evdwij=evdwij*eps2rt*eps3rt evdw=evdw+evdwij*sss if (lprn) then sigm=dabs(aa/bb)**(1.0D0/6.0D0) epsi=bb**2/aa 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=(fac+evdwij*sssgrad/sss/sigmaii(itypi,itypj))*rij c fac=0.0d0 C Calculate the radial part of the gradient gg(1)=xj*fac gg(2)=yj*fac gg(3)=zj*fac gg_lipi(3)=ssgradlipi*evdwij gg_lipj(3)=ssgradlipj*evdwij C Calculate angular part of the gradient. call sc_grad_scale(sss) endif enddo ! j enddo ! iint enddo ! i c write (iout,*) "Number of loop steps in EGB:",ind cccc energy_dec=.false. return end C----------------------------------------------------------------------------- subroutine egbv_long(evdw) C C This subroutine calculates the interaction energy of nonbonded side chains C assuming the Gay-Berne-Vorobjev potential of interaction. C implicit real*8 (a-h,o-z) include 'DIMENSIONS' include 'COMMON.GEO' include 'COMMON.VAR' include 'COMMON.LOCAL' include 'COMMON.CHAIN' include 'COMMON.DERIV' include 'COMMON.NAMES' include 'COMMON.INTERACT' include 'COMMON.IOUNITS' include 'COMMON.CALC' include "COMMON.SPLITELE" integer icall common /srutu/ icall logical lprn integer itypi,itypj,itypi1,iint,ind double precision eps0ij,epsi,sigm,fac,e1,e2,rrij,r0ij, & xi,yi,zi,fac_augm,e_augm double precision fracinbuf,sslipi,evdwij_przed_tri,sig0ij, & sslipj,ssgradlipj,ssgradlipi,dist_init,xj_safe,yj_safe,zj_safe, & xj_temp,yj_temp,zj_temp,dist_temp,sig,rij_shift,faclip double precision dist,sscale,sscagrad,sscagradlip,sscalelip double precision sss1,sssgrad1 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=iabs(itype(i)) if (itypi.eq.ntyp1) cycle itypi1=iabs(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=iabs(itype(j)) if (itypj.eq.ntyp1) cycle c dscj_inv=dsc_inv(itypj) dscj_inv=vbld_inv(j+nres) sig0ij=sigma(itypi,itypj) r0ij=r0(itypi,itypj) chi1=chi(itypi,itypj) chi2=chi(itypj,itypi) chi12=chi1*chi2 chip1=chip(itypi) chip2=chip(itypj) chip12=chip1*chip2 alf1=alp(itypi) alf2=alp(itypj) alf12=0.5D0*(alf1+alf2) xj=c(1,nres+j)-xi yj=c(2,nres+j)-yi zj=c(3,nres+j)-zi dxj=dc_norm(1,nres+j) dyj=dc_norm(2,nres+j) dzj=dc_norm(3,nres+j) rrij=1.0D0/(xj*xj+yj*yj+zj*zj) rij=dsqrt(rrij) sss1=sscale(1.0d0/rij,r_cut_int) if (sss1.eq.0.0d0) cycle if (sss.lt.1.0d0) then sssgrad= & sscagrad((1.0d0/rij)/sigmaii(itypi,itypj),r_cut_respa) sssgrad1=sscagrad(1.0d0/rij,r_cut_int) 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 e2=fac*bb 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)*sss1*(1.0d0-sss) if (lprn) then sigm=dabs(aa/bb)**(1.0D0/6.0D0) epsi=bb**2/aa 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 fac=fac+(evdwij+e_augm)* & (-sss1*sssgrad/(1.0d0-sss)/sigmaii(itypi,itypj) & +(1.0d0-sss)*sssgrad1/sss1)*rij 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)*sss1) 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 none 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.SPLITELE" integer icall common /srutu/ icall logical lprn integer itypi,itypj,itypi1,iint,ind double precision eps0ij,epsi,sigm,fac,e1,e2,rrij,r0ij, & xi,yi,zi,fac_augm,e_augm double precision evdw double precision fracinbuf,sslipi,evdwij_przed_tri,sig0ij, & sslipj,ssgradlipj,ssgradlipi,dist_init,xj_safe,yj_safe,zj_safe, & xj_temp,yj_temp,zj_temp,dist_temp,sig,rij_shift,faclip double precision dist,sscale,sscagrad,sscagradlip,sscalelip 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=iabs(itype(i)) if (itypi.eq.ntyp1) cycle itypi1=iabs(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=iabs(itype(j)) if (itypj.eq.ntyp1) cycle c dscj_inv=dsc_inv(itypj) dscj_inv=vbld_inv(j+nres) sig0ij=sigma(itypi,itypj) r0ij=r0(itypi,itypj) chi1=chi(itypi,itypj) chi2=chi(itypj,itypi) chi12=chi1*chi2 chip1=chip(itypi) chip2=chip(itypj) chip12=chip1*chip2 alf1=alp(itypi) alf2=alp(itypj) alf12=0.5D0*(alf1+alf2) xj=c(1,nres+j)-xi yj=c(2,nres+j)-yi zj=c(3,nres+j)-zi dxj=dc_norm(1,nres+j) dyj=dc_norm(2,nres+j) dzj=dc_norm(3,nres+j) rrij=1.0D0/(xj*xj+yj*yj+zj*zj) rij=dsqrt(rrij) sss=sscale(1.0d0/(rij*sigmaii(itypi,itypj)),r_cut_respa) 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 e2=fac*bb 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/bb)**(1.0D0/6.0D0) epsi=bb**2/aa 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+ & (evdwij+e_augm)*sssgrad/sigmaii(itypi,itypj)/sss*rij 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' include "COMMON.SPLITELE" 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)+gg_lipi(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)+gg_lipj(k) & +(eom12*(dc_norm(k,nres+i)-om12*dc_norm(k,nres+j)) & +eom2*(erij(k)-om2*dc_norm(k,nres+j)))*dscj_inv*scalfac c write (iout,*)(eom12*(dc_norm(k,nres+j)-om12*dc_norm(k,nres+i)) c & +eom1*(erij(k)-om1*dc_norm(k,nres+i)))*dsci_inv c write (iout,*)(eom12*(dc_norm(k,nres+i)-om12*dc_norm(k,nres+j)) c & +eom2*(erij(k)-om2*dc_norm(k,nres+j)))*dscj_inv enddo C C Calculate the components of the gradient in DC and X C do l=1,3 gvdwc(l,i)=gvdwc(l,i)-gg(l)+gg_lipi(l) gvdwc(l,j)=gvdwc(l,j)+gg(l)+gg_lipj(l) enddo return end C-------------------------------------------------------------------------- subroutine eelec_scale(ees,evdw1,eel_loc,eello_turn3,eello_turn4) C C This subroutine calculates the average interaction energy and its gradient C in the virtual-bond vectors between non-adjacent peptide groups, based on C the potential described in Liwo et al., Protein Sci., 1993, 2, 1715. C The potential depends both on the distance of peptide-group centers and on C the orientation of the CA-CA virtual bonds. C implicit real*8 (a-h,o-z) #ifdef MPI include 'mpif.h' #endif include 'DIMENSIONS' include 'COMMON.CONTROL' include 'COMMON.SETUP' include 'COMMON.IOUNITS' include 'COMMON.GEO' include 'COMMON.VAR' include 'COMMON.LOCAL' include 'COMMON.CHAIN' include 'COMMON.DERIV' include 'COMMON.INTERACT' #ifdef FOURBODY include 'COMMON.CONTACTS' include 'COMMON.CONTMAT' #endif include 'COMMON.CORRMAT' include 'COMMON.TORSION' include 'COMMON.VECTORS' include 'COMMON.FFIELD' include 'COMMON.TIME1' include 'COMMON.SHIELD' include "COMMON.SPLITELE" dimension ggg(3),gggp(3),gggm(3),erij(3),dcosb(3),dcosg(3), & erder(3,3),uryg(3,3),urzg(3,3),vryg(3,3),vrzg(3,3) double precision acipa(2,2),agg(3,4),aggi(3,4),aggi1(3,4), & aggj(3,4),aggj1(3,4),a_temp(2,2),muij(4) common /locel/ a_temp,agg,aggi,aggi1,aggj,aggj1,a22,a23,a32,a33, & dxi,dyi,dzi,dx_normi,dy_normi,dz_normi,xmedi,ymedi,zmedi, & num_conti,j1,j2 c 4/26/02 - AL scaling factor for 1,4 repulsive VDW interactions #ifdef MOMENT double precision scal_el /1.0d0/ #else double precision scal_el /0.5d0/ #endif C 12/13/98 C 13-go grudnia roku pamietnego... double precision unmat(3,3) /1.0d0,0.0d0,0.0d0, & 0.0d0,1.0d0,0.0d0, & 0.0d0,0.0d0,1.0d0/ cd write(iout,*) 'In EELEC' cd do i=1,nloctyp cd write(iout,*) 'Type',i cd write(iout,*) 'B1',B1(:,i) cd write(iout,*) 'B2',B2(:,i) cd write(iout,*) 'CC',CC(:,:,i) cd write(iout,*) 'DD',DD(:,:,i) cd write(iout,*) 'EE',EE(:,:,i) cd enddo cd call check_vecgrad cd stop C print *,"WCHODZE3" if (icheckgrad.eq.1) then do i=1,nres-1 fac=1.0d0/dsqrt(scalar(dc(1,i),dc(1,i))) do k=1,3 dc_norm(k,i)=dc(k,i)*fac enddo c write (iout,*) 'i',i,' fac',fac enddo endif if (wel_loc.gt.0.0d0 .or. wcorr4.gt.0.0d0 .or. wcorr5.gt.0.0d0 & .or. wcorr6.gt.0.0d0 .or. wturn3.gt.0.0d0 .or. & wturn4.gt.0.0d0 .or. wturn6.gt.0.0d0) then c call vec_and_deriv #ifdef TIMING time01=MPI_Wtime() #endif call set_matrices #ifdef TIMING time_mat=time_mat+MPI_Wtime()-time01 #endif endif cd do i=1,nres-1 cd write (iout,*) 'i=',i cd do k=1,3 cd write (iout,'(i5,2f10.5)') k,uy(k,i),uz(k,i) cd enddo cd do k=1,3 cd write (iout,'(f10.5,2x,3f10.5,2x,3f10.5)') cd & uz(k,i),(uzgrad(k,l,1,i),l=1,3),(uzgrad(k,l,2,i),l=1,3) cd enddo cd enddo t_eelecij=0.0d0 ees=0.0D0 evdw1=0.0D0 eel_loc=0.0d0 eello_turn3=0.0d0 eello_turn4=0.0d0 ind=0 #ifdef FOURBODY do i=1,nres num_cont_hb(i)=0 enddo #endif cd print '(a)','Enter EELEC' cd write (iout,*) 'iatel_s=',iatel_s,' iatel_e=',iatel_e do i=1,nres gel_loc_loc(i)=0.0d0 gcorr_loc(i)=0.0d0 enddo c c c 9/27/08 AL Split the interaction loop to ensure load balancing of turn terms C C Loop over i,i+2 and i,i+3 pairs of the peptide groups C do i=iturn3_start,iturn3_end if (itype(i).eq.ntyp1.or. itype(i+1).eq.ntyp1 & .or. itype(i+2).eq.ntyp1 .or. itype(i+3).eq.ntyp1 C & .or. itype(i-1).eq.ntyp1 C & .or. itype(i+4).eq.ntyp1 & ) cycle dxi=dc(1,i) dyi=dc(2,i) dzi=dc(3,i) dx_normi=dc_norm(1,i) dy_normi=dc_norm(2,i) dz_normi=dc_norm(3,i) xmedi=c(1,i)+0.5d0*dxi ymedi=c(2,i)+0.5d0*dyi zmedi=c(3,i)+0.5d0*dzi xmedi=mod(xmedi,boxxsize) if (xmedi.lt.0) xmedi=xmedi+boxxsize ymedi=mod(ymedi,boxysize) if (ymedi.lt.0) ymedi=ymedi+boxysize zmedi=mod(zmedi,boxzsize) if (zmedi.lt.0) zmedi=zmedi+boxzsize num_conti=0 call eelecij_scale(i,i+2,ees,evdw1,eel_loc) if (wturn3.gt.0.0d0) call eturn3(i,eello_turn3) #ifdef FOURBODY num_cont_hb(i)=num_conti #endif enddo do i=iturn4_start,iturn4_end if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1 & .or. itype(i+3).eq.ntyp1 & .or. itype(i+4).eq.ntyp1 C & .or. itype(i+5).eq.ntyp1 C & .or. itype(i-1).eq.ntyp1 & ) cycle dxi=dc(1,i) dyi=dc(2,i) dzi=dc(3,i) dx_normi=dc_norm(1,i) dy_normi=dc_norm(2,i) dz_normi=dc_norm(3,i) xmedi=c(1,i)+0.5d0*dxi ymedi=c(2,i)+0.5d0*dyi zmedi=c(3,i)+0.5d0*dzi xmedi=mod(xmedi,boxxsize) if (xmedi.lt.0) xmedi=xmedi+boxxsize ymedi=mod(ymedi,boxysize) if (ymedi.lt.0) ymedi=ymedi+boxysize zmedi=mod(zmedi,boxzsize) if (zmedi.lt.0) zmedi=zmedi+boxzsize #ifdef FOURBODY num_conti=num_cont_hb(i) #endif call eelecij_scale(i,i+3,ees,evdw1,eel_loc) if (wturn4.gt.0.0d0 .and. itype(i+2).ne.ntyp1) & call eturn4(i,eello_turn4) #ifdef FOURBODY num_cont_hb(i)=num_conti #endif enddo ! i c c Loop over all pairs of interacting peptide groups except i,i+2 and i,i+3 c do i=iatel_s,iatel_e if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1 C & .or. itype(i+2).eq.ntyp1 C & .or. itype(i-1).eq.ntyp1 &) cycle dxi=dc(1,i) dyi=dc(2,i) dzi=dc(3,i) dx_normi=dc_norm(1,i) dy_normi=dc_norm(2,i) dz_normi=dc_norm(3,i) xmedi=c(1,i)+0.5d0*dxi ymedi=c(2,i)+0.5d0*dyi zmedi=c(3,i)+0.5d0*dzi xmedi=mod(xmedi,boxxsize) if (xmedi.lt.0) xmedi=xmedi+boxxsize ymedi=mod(ymedi,boxysize) if (ymedi.lt.0) ymedi=ymedi+boxysize zmedi=mod(zmedi,boxzsize) if (zmedi.lt.0) zmedi=zmedi+boxzsize c write (iout,*) 'i',i,' ielstart',ielstart(i),' ielend',ielend #ifdef FOURBODY num_conti=num_cont_hb(i) #endif do j=ielstart(i),ielend(i) if (itype(j).eq.ntyp1 .or. itype(j+1).eq.ntyp1 C & .or.itype(j+2).eq.ntyp1 C & .or.itype(j-1).eq.ntyp1 &) cycle call eelecij_scale(i,j,ees,evdw1,eel_loc) enddo ! j #ifdef FOURBODY num_cont_hb(i)=num_conti #endif enddo ! i c write (iout,*) "Number of loop steps in EELEC:",ind cd do i=1,nres cd write (iout,'(i3,3f10.5,5x,3f10.5)') cd & i,(gel_loc(k,i),k=1,3),gel_loc_loc(i) cd enddo c 12/7/99 Adam eello_turn3 will be considered as a separate energy term ccc eel_loc=eel_loc+eello_turn3 cd print *,"Processor",fg_rank," t_eelecij",t_eelecij return end C------------------------------------------------------------------------------- subroutine eelecij_scale(i,j,ees,evdw1,eel_loc) implicit none include 'DIMENSIONS' #ifdef MPI include "mpif.h" #endif include 'COMMON.CONTROL' include 'COMMON.IOUNITS' include 'COMMON.GEO' include 'COMMON.VAR' include 'COMMON.LOCAL' include 'COMMON.CHAIN' include 'COMMON.DERIV' include 'COMMON.INTERACT' #ifdef FOURBODY include 'COMMON.CONTACTS' include 'COMMON.CONTMAT' #endif include 'COMMON.CORRMAT' include 'COMMON.TORSION' include 'COMMON.VECTORS' include 'COMMON.FFIELD' include 'COMMON.TIME1' include 'COMMON.SHIELD' include "COMMON.SPLITELE" integer xshift,yshift,zshift double precision 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),gmuij1(4),gmuji1(4), & gmuij2(4),gmuji2(4) integer j1,j2,num_conti 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 integer k,i,j,iteli,itelj,kkk,l,kkll,m,isubchap,ind,itypi,itypj integer ilist,iresshield double precision rlocshield double precision ael6i,rrmij,rmij,r0ij,fcont,fprimcont,ees0tmp double precision ees,evdw1,eel_loc,aaa,bbb,ael3i double precision dxj,dyj,dzj,dx_normj,dy_normj,dz_normj,xj,yj,zj, & rij,r3ij,r6ij,cosa,cosb,cosg,fac,ev1,ev2,fac3,fac4, & evdwij,el1,el2,eesij,ees0ij,facvdw,facel,fac1,ecosa, & ecosb,ecosg,ury,urz,vry,vrz,facr,a22der,a23der,a32der, & a33der,eel_loc_ij,cosa4,wij,cosbg1,cosbg2,ees0pij, & ees0pij1,ees0mij,ees0mij1,fac3p,ees0mijp,ees0pijp, & ecosa1,ecosb1,ecosg1,ecosa2,ecosb2,ecosg2,ecosap,ecosbp, & ecosgp,ecosam,ecosbm,ecosgm,ghalf,geel_loc_ij,geel_loc_ji, & dxi,dyi,dzi,a22,a23,a32,a33 double precision dist_init,xmedi,ymedi,zmedi,xj_safe,yj_safe, & zj_safe,xj_temp,yj_temp,zj_temp,dist_temp,dx_normi,dy_normi, & dz_normi,aux double precision sss1,sssgrad1 double precision sscale,sscagrad double precision scalar c 4/26/02 - AL scaling factor for 1,4 repulsive VDW interactions #ifdef MOMENT double precision scal_el /1.0d0/ #else double precision scal_el /0.5d0/ #endif C 12/13/98 C 13-go grudnia roku pamietnego... double precision unmat(3,3) /1.0d0,0.0d0,0.0d0, & 0.0d0,1.0d0,0.0d0, & 0.0d0,0.0d0,1.0d0/ c time00=MPI_Wtime() cd write (iout,*) "eelecij",i,j C print *,"WCHODZE2" ind=ind+1 iteli=itel(i) itelj=itel(j) if (j.eq.i+2 .and. itelj.eq.2) iteli=2 aaa=app(iteli,itelj) bbb=bpp(iteli,itelj) ael6i=ael6(iteli,itelj) ael3i=ael3(iteli,itelj) dxj=dc(1,j) dyj=dc(2,j) dzj=dc(3,j) dx_normj=dc_norm(1,j) dy_normj=dc_norm(2,j) dz_normj=dc_norm(3,j) xj=c(1,j)+0.5D0*dxj yj=c(2,j)+0.5D0*dyj zj=c(3,j)+0.5D0*dzj xj=mod(xj,boxxsize) if (xj.lt.0) xj=xj+boxxsize yj=mod(yj,boxysize) if (yj.lt.0) yj=yj+boxysize zj=mod(zj,boxzsize) if (zj.lt.0) zj=zj+boxzsize dist_init=(xj-xmedi)**2+(yj-ymedi)**2+(zj-zmedi)**2 xj_safe=xj yj_safe=yj zj_safe=zj isubchap=0 do xshift=-1,1 do yshift=-1,1 do zshift=-1,1 xj=xj_safe+xshift*boxxsize yj=yj_safe+yshift*boxysize zj=zj_safe+zshift*boxzsize dist_temp=(xj-xmedi)**2+(yj-ymedi)**2+(zj-zmedi)**2 if(dist_temp.lt.dist_init) then dist_init=dist_temp xj_temp=xj yj_temp=yj zj_temp=zj isubchap=1 endif enddo enddo enddo if (isubchap.eq.1) then xj=xj_temp-xmedi yj=yj_temp-ymedi zj=zj_temp-zmedi else xj=xj_safe-xmedi yj=yj_safe-ymedi zj=zj_safe-zmedi endif 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 sss1=sscale(rij,r_cut_int) if (sss1.eq.0.0d0) return sss=sscale(rij/rpp(iteli,itelj),r_cut_respa) sssgrad=sscagrad(rij/rpp(iteli,itelj),r_cut_respa) sssgrad1=sscagrad(rij,r_cut_int) 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 if (shield_mode.eq.0) then fac_shield(i)=1.0 fac_shield(j)=1.0 endif el1=fac3*(4.0D0+fac*fac-3.0D0*(cosb*cosb+cosg*cosg)) el2=fac4*fac el1=el1*fac_shield(i)**2*fac_shield(j)**2 el2=el2*fac_shield(i)**2*fac_shield(j)**2 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*sss1 evdw1=evdw1+evdwij*(1.0d0-sss)*sss1 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,2f7.3)') & 'evdw1',i,j,evdwij,sss,sss1 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)*sss1 c & *((sslipi+sslipj)/2.0d0*lipscale**2+1.0d0) facel=-3*rrmij*(el1+eesij)*sss1 c & *((sslipi+sslipj)/2.0d0*lipscale**2+1.0d0) 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) * aux=facel+sssgrad1*(1.0d0-sss)*eesij*rmij c & *((sslipi+sslipj)/2.0d0*lipscale**2+1.0d0) ggg(1)=aux*xj ggg(2)=aux*yj ggg(3)=aux*zj c ggg(1)=facel*xj c ggg(2)=facel*yj c ggg(3)=facel*zj if ((fac_shield(i).gt.0).and.(fac_shield(j).gt.0).and. & (shield_mode.gt.0)) then C print *,i,j do ilist=1,ishield_list(i) iresshield=shield_list(ilist,i) do k=1,3 rlocshield=grad_shield_side(k,ilist,i)*eesij*sss1 & /fac_shield(i)*2.0*sss1 gshieldx(k,iresshield)=gshieldx(k,iresshield)+ & rlocshield & +grad_shield_loc(k,ilist,i)*eesij*sss1/fac_shield(i)*2.0 & *sss1 gshieldc(k,iresshield-1)=gshieldc(k,iresshield-1)+rlocshield C gshieldc_loc(k,iresshield)=gshieldc_loc(k,iresshield) C & +grad_shield_loc(k,ilist,i)*eesij/fac_shield(i) C if (iresshield.gt.i) then C do ishi=i+1,iresshield-1 C gshieldc(k,ishi)=gshieldc(k,ishi)+rlocshield C & +grad_shield_loc(k,ilist,i)*eesij/fac_shield(i) C C enddo C else C do ishi=iresshield,i C gshieldc(k,ishi)=gshieldc(k,ishi)-rlocshield C & -grad_shield_loc(k,ilist,i)*eesij/fac_shield(i) C C enddo C endif enddo enddo do ilist=1,ishield_list(j) iresshield=shield_list(ilist,j) do k=1,3 rlocshield=grad_shield_side(k,ilist,j)*eesij/fac_shield(j) & *2.0*sss1 gshieldx(k,iresshield)=gshieldx(k,iresshield)+ & rlocshield & +grad_shield_loc(k,ilist,j)*eesij/fac_shield(j)*2.0*sss1 gshieldc(k,iresshield-1)=gshieldc(k,iresshield-1)+rlocshield enddo enddo do k=1,3 gshieldc(k,i)=gshieldc(k,i)+ & grad_shield(k,i)*eesij/fac_shield(i)*2.0*sss1 gshieldc(k,j)=gshieldc(k,j)+ & grad_shield(k,j)*eesij/fac_shield(j)*2.0*sss1 gshieldc(k,i-1)=gshieldc(k,i-1)+ & grad_shield(k,i)*eesij/fac_shield(i)*2.0*sss1 gshieldc(k,j-1)=gshieldc(k,j-1)+ & grad_shield(k,j)*eesij/fac_shield(j)*2.0*sss1 enddo endif c do k=1,3 c ghalf=0.5D0*ggg(k) c gelc(k,i)=gelc(k,i)+ghalf c gelc(k,j)=gelc(k,j)+ghalf c enddo c 9/28/08 AL Gradient compotents will be summed only at the end do k=1,3 gelc_long(k,j)=gelc_long(k,j)+ggg(k) gelc_long(k,i)=gelc_long(k,i)-ggg(k) enddo c gelc_long(3,i)=gelc_long(3,i)+ c ssgradlipi*eesij/2.0d0*lipscale**2*sss1 * * Loop over residues i+1 thru j-1. * cgrad do k=i+1,j-1 cgrad do l=1,3 cgrad gelc(l,k)=gelc(l,k)+ggg(l) cgrad enddo cgrad enddo facvdw=facvdw+ & (-sss1*sssgrad/rpp(iteli,itelj)+(1.0d0-sss)*sssgrad1)*rmij*evdwij c & *((sslipi+sslipj)/2.0d0*lipscale**2+1.0d0) ggg(1)=facvdw*xj ggg(2)=facvdw*yj ggg(3)=facvdw*zj c do k=1,3 c ghalf=0.5D0*ggg(k) c gvdwpp(k,i)=gvdwpp(k,i)+ghalf c gvdwpp(k,j)=gvdwpp(k,j)+ghalf c enddo c 9/28/08 AL Gradient compotents will be summed only at the end do k=1,3 gvdwpp(k,j)=gvdwpp(k,j)+ggg(k) gvdwpp(k,i)=gvdwpp(k,i)-ggg(k) enddo * * Loop over residues i+1 thru j-1. * cgrad do k=i+1,j-1 cgrad do l=1,3 cgrad gvdwpp(l,k)=gvdwpp(l,k)+ggg(l) cgrad enddo cgrad enddo #else facvdw=-6*rrmij*(ev1+evdwij)*(1.0d0-sss)*sss1 c & *((sslipi+sslipj)/2.0d0*lipscale**2+1.0d0) facel=-3*rrmij*(el1+eesij)*sss1 c & *((sslipi+sslipj)/2.0d0*lipscale**2+1.0d0) c facvdw=ev1+evdwij*(1.0d0-sss)*sss1 c 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) * aux=fac+(sssgrad1*(1.0d0-sss)-sssgrad*sss1/rpp(iteli,itelj)) & *eesij*rmij c & *((sslipi+sslipj)/2.0d0*lipscale**2+1.0d0) ggg(1)=aux*xj ggg(2)=aux*yj ggg(3)=axu*zj c ggg(1)=fac*xj c ggg(2)=fac*yj c ggg(3)=fac*zj c do k=1,3 c ghalf=0.5D0*ggg(k) c gelc(k,i)=gelc(k,i)+ghalf c gelc(k,j)=gelc(k,j)+ghalf c enddo c 9/28/08 AL Gradient compotents will be summed only at the end do k=1,3 gelc_long(k,j)=gelc(k,j)+ggg(k) gelc_long(k,i)=gelc(k,i)-ggg(k) enddo * * Loop over residues i+1 thru j-1. * cgrad do k=i+1,j-1 cgrad do l=1,3 cgrad gelc(l,k)=gelc(l,k)+ggg(l) cgrad enddo cgrad enddo c 9/28/08 AL Gradient compotents will be summed only at the end C ggg(1)=facvdw*xj C ggg(2)=facvdw*yj C ggg(3)=facvdw*zj facvdw=facvdw & (-sssgrad*sss1/rpp(iteli,itelj)+sssgrad1*(1.0d0-sss))*rmij*evdwij ggg(1)=facvdw*xj ggg(2)=facvdw*yj ggg(3)=facvdw*zj do k=1,3 gvdwpp(k,j)=gvdwpp(k,j)+ggg(k) gvdwpp(k,i)=gvdwpp(k,i)-ggg(k) enddo #endif * * Angular part * ecosa=2.0D0*fac3*fac1+fac4 fac4=-3.0D0*fac4 fac3=-6.0D0*fac3 ecosb=(fac3*(fac1*cosg+cosb)+cosg*fac4) ecosg=(fac3*(fac1*cosb+cosg)+cosb*fac4) do k=1,3 dcosb(k)=rmij*(dc_norm(k,i)-erij(k)*cosb) dcosg(k)=rmij*(dc_norm(k,j)-erij(k)*cosg) enddo cd print '(2i3,2(3(1pd14.5),3x))',i,j,(dcosb(k),k=1,3), cd & (dcosg(k),k=1,3) do k=1,3 ggg(k)=(ecosb*dcosb(k)+ecosg*dcosg(k))*sss1 & *fac_shield(i)**2*fac_shield(j)**2 c & *((sslipi+sslipj)/2.0d0*lipscale**2+1.0d0) enddo c do k=1,3 c ghalf=0.5D0*ggg(k) c gelc(k,i)=gelc(k,i)+ghalf c & +(ecosa*(dc_norm(k,j)-cosa*dc_norm(k,i)) c & + ecosb*(erij(k)-cosb*dc_norm(k,i)))*vbld_inv(i+1) c gelc(k,j)=gelc(k,j)+ghalf c & +(ecosa*(dc_norm(k,i)-cosa*dc_norm(k,j)) c & + ecosg*(erij(k)-cosg*dc_norm(k,j)))*vbld_inv(j+1) c enddo cgrad do k=i+1,j-1 cgrad do l=1,3 cgrad gelc(l,k)=gelc(l,k)+ggg(l) cgrad enddo cgrad enddo do k=1,3 gelc(k,i)=gelc(k,i) & +((ecosa*(dc_norm(k,j)-cosa*dc_norm(k,i)) & +ecosb*(erij(k)-cosb*dc_norm(k,i)))*vbld_inv(i+1))*sss1 & *fac_shield(i)**2*fac_shield(j)**2 c & *((sslipi+sslipj)/2.0d0*lipscale**2+1.0d0) gelc(k,j)=gelc(k,j) & +((ecosa*(dc_norm(k,i)-cosa*dc_norm(k,j)) & +ecosg*(erij(k)-cosg*dc_norm(k,j)))*vbld_inv(j+1))*sss1 & *fac_shield(i)**2*fac_shield(j)**2 c & *((sslipi+sslipj)/2.0d0*lipscale**2+1.0d0) gelc_long(k,j)=gelc_long(k,j)+ggg(k) gelc_long(k,i)=gelc_long(k,i)-ggg(k) enddo IF (wel_loc.gt.0.0d0 .or. wcorr4.gt.0.0d0 .or. wcorr5.gt.0.0d0 & .or. wcorr6.gt.0.0d0 .or. wturn3.gt.0.0d0 & .or. wturn4.gt.0.0d0 .or. wturn6.gt.0.0d0) THEN C C 9/25/99 Mixed third-order local-electrostatic terms. The local-interaction C energy of a peptide unit is assumed in the form of a second-order C Fourier series in the angles lambda1 and lambda2 (see Nishikawa et al. C Macromolecules, 1974, 7, 797-806 for definition). This correlation terms C are computed for EVERY pair of non-contiguous peptide groups. C if (j.lt.nres-1) then j1=j+1 j2=j-1 else j1=j-1 j2=j-2 endif kkk=0 do k=1,2 do l=1,2 kkk=kkk+1 muij(kkk)=mu(k,i)*mu(l,j) #ifdef NEWCORR gmuij1(kkk)=gtb1(k,i+1)*mu(l,j) c write(iout,*) 'k=',k,i,gtb1(k,i+1),gtb1(k,i+1)*mu(l,j) gmuij2(kkk)=gUb2(k,i)*mu(l,j) gmuji1(kkk)=mu(k,i)*gtb1(l,j+1) c write(iout,*) 'l=',l,j,gtb1(l,j+1),gtb1(l,j+1)*mu(k,i) gmuji2(kkk)=mu(k,i)*gUb2(l,j) #endif enddo enddo cd write (iout,*) 'EELEC: i',i,' j',j cd write (iout,*) 'j',j,' j1',j1,' j2',j2 cd write(iout,*) 'muij',muij ury=scalar(uy(1,i),erij) urz=scalar(uz(1,i),erij) vry=scalar(uy(1,j),erij) vrz=scalar(uz(1,j),erij) a22=scalar(uy(1,i),uy(1,j))-3*ury*vry a23=scalar(uy(1,i),uz(1,j))-3*ury*vrz a32=scalar(uz(1,i),uy(1,j))-3*urz*vry a33=scalar(uz(1,i),uz(1,j))-3*urz*vrz fac=dsqrt(-ael6i)*r3ij a22=a22*fac a23=a23*fac a32=a32*fac a33=a33*fac cd write (iout,'(4i5,4f10.5)') cd & i,itortyp(itype(i)),j,itortyp(itype(j)),a22,a23,a32,a33 cd write (iout,'(6f10.5)') (muij(k),k=1,4),fac,eel_loc_ij cd write (iout,'(2(3f10.5,5x)/2(3f10.5,5x))') uy(:,i),uz(:,i), cd & uy(:,j),uz(:,j) cd write (iout,'(4f10.5)') cd & scalar(uy(1,i),uy(1,j)),scalar(uy(1,i),uz(1,j)), cd & scalar(uz(1,i),uy(1,j)),scalar(uz(1,i),uz(1,j)) cd write (iout,'(4f10.5)') ury,urz,vry,vrz cd write (iout,'(9f10.5/)') cd & fac22,a22,fac23,a23,fac32,a32,fac33,a33,eel_loc_ij C Derivatives of the elements of A in virtual-bond vectors call unormderiv(erij(1),unmat(1,1),rmij,erder(1,1)) do k=1,3 uryg(k,1)=scalar(erder(1,k),uy(1,i)) uryg(k,2)=scalar(uygrad(1,k,1,i),erij(1)) uryg(k,3)=scalar(uygrad(1,k,2,i),erij(1)) urzg(k,1)=scalar(erder(1,k),uz(1,i)) urzg(k,2)=scalar(uzgrad(1,k,1,i),erij(1)) urzg(k,3)=scalar(uzgrad(1,k,2,i),erij(1)) vryg(k,1)=scalar(erder(1,k),uy(1,j)) vryg(k,2)=scalar(uygrad(1,k,1,j),erij(1)) vryg(k,3)=scalar(uygrad(1,k,2,j),erij(1)) vrzg(k,1)=scalar(erder(1,k),uz(1,j)) vrzg(k,2)=scalar(uzgrad(1,k,1,j),erij(1)) vrzg(k,3)=scalar(uzgrad(1,k,2,j),erij(1)) enddo C Compute radial contributions to the gradient facr=-3.0d0*rrmij a22der=a22*facr a23der=a23*facr a32der=a32*facr a33der=a33*facr agg(1,1)=a22der*xj agg(2,1)=a22der*yj agg(3,1)=a22der*zj agg(1,2)=a23der*xj agg(2,2)=a23der*yj agg(3,2)=a23der*zj agg(1,3)=a32der*xj agg(2,3)=a32der*yj agg(3,3)=a32der*zj agg(1,4)=a33der*xj agg(2,4)=a33der*yj agg(3,4)=a33der*zj C Add the contributions coming from er fac3=-3.0d0*fac do k=1,3 agg(k,1)=agg(k,1)+fac3*(uryg(k,1)*vry+vryg(k,1)*ury) agg(k,2)=agg(k,2)+fac3*(uryg(k,1)*vrz+vrzg(k,1)*ury) agg(k,3)=agg(k,3)+fac3*(urzg(k,1)*vry+vryg(k,1)*urz) agg(k,4)=agg(k,4)+fac3*(urzg(k,1)*vrz+vrzg(k,1)*urz) enddo do k=1,3 C Derivatives in DC(i) cgrad ghalf1=0.5d0*agg(k,1) cgrad ghalf2=0.5d0*agg(k,2) cgrad ghalf3=0.5d0*agg(k,3) cgrad ghalf4=0.5d0*agg(k,4) aggi(k,1)=fac*(scalar(uygrad(1,k,1,i),uy(1,j)) & -3.0d0*uryg(k,2)*vry)!+ghalf1 aggi(k,2)=fac*(scalar(uygrad(1,k,1,i),uz(1,j)) & -3.0d0*uryg(k,2)*vrz)!+ghalf2 aggi(k,3)=fac*(scalar(uzgrad(1,k,1,i),uy(1,j)) & -3.0d0*urzg(k,2)*vry)!+ghalf3 aggi(k,4)=fac*(scalar(uzgrad(1,k,1,i),uz(1,j)) & -3.0d0*urzg(k,2)*vrz)!+ghalf4 C Derivatives in DC(i+1) aggi1(k,1)=fac*(scalar(uygrad(1,k,2,i),uy(1,j)) & -3.0d0*uryg(k,3)*vry)!+agg(k,1) aggi1(k,2)=fac*(scalar(uygrad(1,k,2,i),uz(1,j)) & -3.0d0*uryg(k,3)*vrz)!+agg(k,2) aggi1(k,3)=fac*(scalar(uzgrad(1,k,2,i),uy(1,j)) & -3.0d0*urzg(k,3)*vry)!+agg(k,3) aggi1(k,4)=fac*(scalar(uzgrad(1,k,2,i),uz(1,j)) & -3.0d0*urzg(k,3)*vrz)!+agg(k,4) C Derivatives in DC(j) aggj(k,1)=fac*(scalar(uygrad(1,k,1,j),uy(1,i)) & -3.0d0*vryg(k,2)*ury)!+ghalf1 aggj(k,2)=fac*(scalar(uzgrad(1,k,1,j),uy(1,i)) & -3.0d0*vrzg(k,2)*ury)!+ghalf2 aggj(k,3)=fac*(scalar(uygrad(1,k,1,j),uz(1,i)) & -3.0d0*vryg(k,2)*urz)!+ghalf3 aggj(k,4)=fac*(scalar(uzgrad(1,k,1,j),uz(1,i)) & -3.0d0*vrzg(k,2)*urz)!+ghalf4 C Derivatives in DC(j+1) or DC(nres-1) aggj1(k,1)=fac*(scalar(uygrad(1,k,2,j),uy(1,i)) & -3.0d0*vryg(k,3)*ury) aggj1(k,2)=fac*(scalar(uzgrad(1,k,2,j),uy(1,i)) & -3.0d0*vrzg(k,3)*ury) aggj1(k,3)=fac*(scalar(uygrad(1,k,2,j),uz(1,i)) & -3.0d0*vryg(k,3)*urz) aggj1(k,4)=fac*(scalar(uzgrad(1,k,2,j),uz(1,i)) & -3.0d0*vrzg(k,3)*urz) cgrad if (j.eq.nres-1 .and. i.lt.j-2) then cgrad do l=1,4 cgrad aggj1(k,l)=aggj1(k,l)+agg(k,l) cgrad enddo cgrad endif enddo acipa(1,1)=a22 acipa(1,2)=a23 acipa(2,1)=a32 acipa(2,2)=a33 a22=-a22 a23=-a23 do l=1,2 do k=1,3 agg(k,l)=-agg(k,l) aggi(k,l)=-aggi(k,l) aggi1(k,l)=-aggi1(k,l) aggj(k,l)=-aggj(k,l) aggj1(k,l)=-aggj1(k,l) enddo enddo if (j.lt.nres-1) then a22=-a22 a32=-a32 do l=1,3,2 do k=1,3 agg(k,l)=-agg(k,l) aggi(k,l)=-aggi(k,l) aggi1(k,l)=-aggi1(k,l) aggj(k,l)=-aggj(k,l) aggj1(k,l)=-aggj1(k,l) enddo enddo else a22=-a22 a23=-a23 a32=-a32 a33=-a33 do l=1,4 do k=1,3 agg(k,l)=-agg(k,l) aggi(k,l)=-aggi(k,l) aggi1(k,l)=-aggi1(k,l) aggj(k,l)=-aggj(k,l) aggj1(k,l)=-aggj1(k,l) enddo enddo endif ENDIF ! WCORR IF (wel_loc.gt.0.0d0) THEN C Contribution to the local-electrostatic energy coming from the i-j pair eel_loc_ij=a22*muij(1)+a23*muij(2)+a32*muij(3) & +a33*muij(4) cd write (iout,*) 'i',i,' j',j,' eel_loc_ij',eel_loc_ij if (energy_dec) write (iout,'(a6,2i5,0pf7.3)') & 'eelloc',i,j,eel_loc_ij if (shield_mode.eq.0) then fac_shield(i)=1.0 fac_shield(j)=1.0 C else C fac_shield(i)=0.4 C fac_shield(j)=0.6 endif eel_loc_ij=eel_loc_ij & *fac_shield(i)*fac_shield(j)*sss1 eel_loc=eel_loc+eel_loc_ij if ((fac_shield(i).gt.0).and.(fac_shield(j).gt.0).and. & (shield_mode.gt.0)) then C print *,i,j do ilist=1,ishield_list(i) iresshield=shield_list(ilist,i) do k=1,3 rlocshield=grad_shield_side(k,ilist,i)*eel_loc_ij & /fac_shield(i) C & *2.0 gshieldx_ll(k,iresshield)=gshieldx_ll(k,iresshield)+ & rlocshield & +grad_shield_loc(k,ilist,i)*eel_loc_ij/fac_shield(i) gshieldc_ll(k,iresshield-1)=gshieldc_ll(k,iresshield-1) & +rlocshield enddo enddo do ilist=1,ishield_list(j) iresshield=shield_list(ilist,j) do k=1,3 rlocshield=grad_shield_side(k,ilist,j)*eel_loc_ij & /fac_shield(j) C & *2.0 gshieldx_ll(k,iresshield)=gshieldx_ll(k,iresshield)+ & rlocshield & +grad_shield_loc(k,ilist,j)*eel_loc_ij/fac_shield(j) gshieldc_ll(k,iresshield-1)=gshieldc_ll(k,iresshield-1) & +rlocshield enddo enddo do k=1,3 gshieldc_ll(k,i)=gshieldc_ll(k,i)+ & grad_shield(k,i)*eel_loc_ij/fac_shield(i) gshieldc_ll(k,j)=gshieldc_ll(k,j)+ & grad_shield(k,j)*eel_loc_ij/fac_shield(j) gshieldc_ll(k,i-1)=gshieldc_ll(k,i-1)+ & grad_shield(k,i)*eel_loc_ij/fac_shield(i) gshieldc_ll(k,j-1)=gshieldc_ll(k,j-1)+ & grad_shield(k,j)*eel_loc_ij/fac_shield(j) enddo endif #ifdef NEWCORR geel_loc_ij=(a22*gmuij1(1) & +a23*gmuij1(2) & +a32*gmuij1(3) & +a33*gmuij1(4)) & *fac_shield(i)*fac_shield(j)*sss1 c write(iout,*) "derivative over thatai" c write(iout,*) a22*gmuij1(1), a23*gmuij1(2) ,a32*gmuij1(3), c & a33*gmuij1(4) gloc(nphi+i,icg)=gloc(nphi+i,icg)+ & geel_loc_ij*wel_loc c write(iout,*) "derivative over thatai-1" c write(iout,*) a22*gmuij2(1), a23*gmuij2(2) ,a32*gmuij2(3), c & a33*gmuij2(4) geel_loc_ij= & a22*gmuij2(1) & +a23*gmuij2(2) & +a32*gmuij2(3) & +a33*gmuij2(4) gloc(nphi+i-1,icg)=gloc(nphi+i-1,icg)+ & geel_loc_ij*wel_loc & *fac_shield(i)*fac_shield(j)*sss1 c Derivative over j residue geel_loc_ji=a22*gmuji1(1) & +a23*gmuji1(2) & +a32*gmuji1(3) & +a33*gmuji1(4) c write(iout,*) "derivative over thataj" c write(iout,*) a22*gmuji1(1), a23*gmuji1(2) ,a32*gmuji1(3), c & a33*gmuji1(4) gloc(nphi+j,icg)=gloc(nphi+j,icg)+ & geel_loc_ji*wel_loc & *fac_shield(i)*fac_shield(j)*sss1 geel_loc_ji= & +a22*gmuji2(1) & +a23*gmuji2(2) & +a32*gmuji2(3) & +a33*gmuji2(4) c write(iout,*) "derivative over thataj-1" c write(iout,*) a22*gmuji2(1), a23*gmuji2(2) ,a32*gmuji2(3), c & a33*gmuji2(4) gloc(nphi+j-1,icg)=gloc(nphi+j-1,icg)+ & geel_loc_ji*wel_loc & *fac_shield(i)*fac_shield(j)*sss1 #endif cC Paral 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)) & *fac_shield(i)*fac_shield(j)*sss1 c & *fac_shield(i)*fac_shield(j) c & *((sslipi+sslipj)/2.0d0*lipscale+1.0d0) 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)) & *fac_shield(i)*fac_shield(j)*sss1 c & *fac_shield(i)*fac_shield(j) c & *((sslipi+sslipj)/2.0d0*lipscale+1.0d0) C Derivatives of eello in DC(i+1) thru DC(j-1) or DC(nres-2) aux=eel_loc_ij/sss1*sssgrad1*rmij ggg(1)=aux*xj ggg(2)=aux*yj ggg(3)=aux*zj do l=1,3 ggg(l)=ggg(l)+(agg(l,1)*muij(1)+ & agg(l,2)*muij(2)+agg(l,3)*muij(3)+agg(l,4)*muij(4)) & *fac_shield(i)*fac_shield(j)*sss1 c & *fac_shield(i)*fac_shield(j) c & *((sslipi+sslipj)/2.0d0*lipscale+1.0d0) gel_loc_long(l,j)=gel_loc_long(l,j)+ggg(l) gel_loc_long(l,i)=gel_loc_long(l,i)-ggg(l) cgrad ghalf=0.5d0*ggg(l) cgrad gel_loc(l,i)=gel_loc(l,i)+ghalf cgrad gel_loc(l,j)=gel_loc(l,j)+ghalf enddo cgrad do k=i+1,j2 cgrad do l=1,3 cgrad gel_loc(l,k)=gel_loc(l,k)+ggg(l) cgrad enddo cgrad enddo C Remaining derivatives of eello c gel_loc_long(3,j)=gel_loc_long(3,j)+ & c ssgradlipj*eel_loc_ij/2.0d0*lipscale/ & c ((sslipi+sslipj)/2.0d0*lipscale+1.0d0)*sss_ele_cut c c gel_loc_long(3,i)=gel_loc_long(3,i)+ & c ssgradlipi*eel_loc_ij/2.0d0*lipscale/ & c ((sslipi+sslipj)/2.0d0*lipscale+1.0d0)*sss_ele_cut 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)) & *fac_shield(i)*fac_shield(j)*sss1 c & *((sslipi+sslipj)/2.0d0*lipscale+1.0d0) 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)) & *fac_shield(i)*fac_shield(j)*sss1 c & *((sslipi+sslipj)/2.0d0*lipscale+1.0d0) 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)) & *fac_shield(i)*fac_shield(j)*sss1 c & *((sslipi+sslipj)/2.0d0*lipscale+1.0d0) 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)) & *fac_shield(i)*fac_shield(j)*sss1 c & *((sslipi+sslipj)/2.0d0*lipscale+1.0d0) enddo ENDIF #ifdef FOURBODY C Change 12/26/95 to calculate four-body contributions to H-bonding energy c if (j.gt.i+1 .and. num_conti.le.maxconts) then if (wcorr+wcorr4+wcorr5+wcorr6.gt.0.0d0 & .and. num_conti.le.maxconts) then c write (iout,*) i,j," entered corr" C C Calculate the contact function. The ith column of the array JCONT will C contain the numbers of atoms that make contacts with the atom I (of numbers C greater than I). The arrays FACONT and GACONT will contain the values of C the contact function and its derivative. c r0ij=1.02D0*rpp(iteli,itelj) c r0ij=1.11D0*rpp(iteli,itelj) r0ij=2.20D0*rpp(iteli,itelj) c r0ij=1.55D0*rpp(iteli,itelj) call gcont(rij,r0ij,1.0D0,0.2d0*r0ij,fcont,fprimcont) if (fcont.gt.0.0D0) then num_conti=num_conti+1 if (num_conti.gt.maxconts) then write (iout,*) 'WARNING - max. # of contacts exceeded;', & ' will skip next contacts for this conf.' else jcont_hb(num_conti,i)=j cd write (iout,*) "i",i," j",j," num_conti",num_conti, cd " jcont_hb",jcont_hb(num_conti,i) IF (wcorr4.gt.0.0d0 .or. wcorr5.gt.0.0d0 .or. & wcorr6.gt.0.0d0 .or. wturn6.gt.0.0d0) THEN C 9/30/99 (AL) - store components necessary to evaluate higher-order loc-el C terms. d_cont(num_conti,i)=rij cd write (2,'(3e15.5)') rij,r0ij+0.2d0*r0ij,rij C --- Electrostatic-interaction matrix --- a_chuj(1,1,num_conti,i)=a22 a_chuj(1,2,num_conti,i)=a23 a_chuj(2,1,num_conti,i)=a32 a_chuj(2,2,num_conti,i)=a33 C --- Gradient of rij do kkk=1,3 grij_hb_cont(kkk,num_conti,i)=erij(kkk) enddo kkll=0 do k=1,2 do l=1,2 kkll=kkll+1 do m=1,3 a_chuj_der(k,l,m,1,num_conti,i)=agg(m,kkll) a_chuj_der(k,l,m,2,num_conti,i)=aggi(m,kkll) a_chuj_der(k,l,m,3,num_conti,i)=aggi1(m,kkll) a_chuj_der(k,l,m,4,num_conti,i)=aggj(m,kkll) a_chuj_der(k,l,m,5,num_conti,i)=aggj1(m,kkll) enddo enddo enddo ENDIF IF (wcorr4.eq.0.0d0 .and. wcorr.gt.0.0d0) THEN C Calculate contact energies cosa4=4.0D0*cosa wij=cosa-3.0D0*cosb*cosg cosbg1=cosb+cosg cosbg2=cosb-cosg c fac3=dsqrt(-ael6i)/r0ij**3 fac3=dsqrt(-ael6i)*r3ij c ees0pij=dsqrt(4.0D0+cosa4+wij*wij-3.0D0*cosbg1*cosbg1) ees0tmp=4.0D0+cosa4+wij*wij-3.0D0*cosbg1*cosbg1 if (ees0tmp.gt.0) then ees0pij=dsqrt(ees0tmp) else ees0pij=0 endif c ees0mij=dsqrt(4.0D0-cosa4+wij*wij-3.0D0*cosbg2*cosbg2) ees0tmp=4.0D0-cosa4+wij*wij-3.0D0*cosbg2*cosbg2 if (ees0tmp.gt.0) then ees0mij=dsqrt(ees0tmp) else ees0mij=0 endif c ees0mij=0.0D0 if (shield_mode.eq.0) then fac_shield(i)=1.0d0 fac_shield(j)=1.0d0 else ees0plist(num_conti,i)=j C fac_shield(i)=0.4d0 C fac_shield(j)=0.6d0 endif ees0p(num_conti,i)=0.5D0*fac3*(ees0pij+ees0mij) & *fac_shield(i)*fac_shield(j)*sss1 ees0m(num_conti,i)=0.5D0*fac3*(ees0pij-ees0mij) & *fac_shield(i)*fac_shield(j)*sss1 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 & +ees0p(num_conti,i)/sss1*rmij*xj*sssgrad1 gggp(2)=gggp(2)+ees0pijp*yj & +ees0p(num_conti,i)/sss1*rmij*yj*sssgrad1 gggp(3)=gggp(3)+ees0pijp*zj & +ees0p(num_conti,i)/sss1*rmij*zj*sssgrad1 gggm(1)=gggm(1)+ees0mijp*xj & +ees0m(num_conti,i)/sss1*rmij*xj*sssgrad1 gggm(2)=gggm(2)+ees0mijp*yj & +ees0m(num_conti,i)/sss1*rmij*yj*sssgrad1 gggm(3)=gggm(3)+ees0mijp*zj & +ees0m(num_conti,i)/sss1*rmij*zj*sssgrad1 C Derivatives due to the contact function gacont_hbr(1,num_conti,i)=fprimcont*xj gacont_hbr(2,num_conti,i)=fprimcont*yj gacont_hbr(3,num_conti,i)=fprimcont*zj do k=1,3 c c 10/24/08 cgrad and ! comments indicate the parts of the code removed c following the change of gradient-summation algorithm. c cgrad ghalfp=0.5D0*gggp(k) cgrad ghalfm=0.5D0*gggm(k) gacontp_hb1(k,num_conti,i)=!ghalfp & +(ecosap*(dc_norm(k,j)-cosa*dc_norm(k,i)) & + ecosbp*(erij(k)-cosb*dc_norm(k,i)))*vbld_inv(i+1) & *sss1*fac_shield(i)*fac_shield(j) 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) & *sss1*fac_shield(i)*fac_shield(j) gacontp_hb3(k,num_conti,i)=gggp(k) & *sss1*fac_shield(i)*fac_shield(j) 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) & *sss1*fac_shield(i)*fac_shield(j) 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) & *sss1*fac_shield(i)*fac_shield(j) gacontm_hb3(k,num_conti,i)=gggm(k) & *sss1*fac_shield(i)*fac_shield(j) enddo ENDIF ! wcorr endif ! num_conti.le.maxconts endif ! fcont.gt.0 endif ! j.gt.i+1 #endif if (wturn3.gt.0.0d0 .or. wturn4.gt.0.0d0) then do k=1,4 do l=1,3 ghalf=0.5d0*agg(l,k) aggi(l,k)=aggi(l,k)+ghalf aggi1(l,k)=aggi1(l,k)+agg(l,k) aggj(l,k)=aggj(l,k)+ghalf enddo enddo if (j.eq.nres-1 .and. i.lt.j-2) then do k=1,4 do l=1,3 aggj1(l,k)=aggj1(l,k)+agg(l,k) enddo enddo endif endif c t_eelecij=t_eelecij+MPI_Wtime()-time00 return end C----------------------------------------------------------------------- subroutine evdwpp_short(evdw1) C C Compute Evdwpp C implicit none 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' c include 'COMMON.CONTACTS' include 'COMMON.TORSION' include 'COMMON.VECTORS' include 'COMMON.FFIELD' include "COMMON.SPLITELE" double precision ggg(3) integer xshift,yshift,zshift 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 write (iout,*) "evdwpp_short" integer i,j,k,iteli,itelj,num_conti,ind,isubchap double precision dxi,dyi,dzi,dxj,dyj,dzj,aaa,bbb double precision xj,yj,zj,rij,rrmij,r3ij,r6ij,evdw1, & dx_normi,dy_normi,dz_normi,xmedi,ymedi,zmedi, & dx_normj,dy_normj,dz_normj,rmij,ev1,ev2,evdwij,facvdw double precision xj_safe,yj_safe,zj_safe,xj_temp,yj_temp,zj_temp, & dist_temp, dist_init,sss_grad double precision sscale,sscagrad evdw1=0.0D0 C print *,"WCHODZE" c write (iout,*) "iatel_s_vdw",iatel_s_vdw, c & " iatel_e_vdw",iatel_e_vdw c call flush(iout) do i=iatel_s_vdw,iatel_e_vdw if (itype(i).eq.ntyp1.or. itype(i+1).eq.ntyp1) cycle dxi=dc(1,i) dyi=dc(2,i) dzi=dc(3,i) dx_normi=dc_norm(1,i) dy_normi=dc_norm(2,i) dz_normi=dc_norm(3,i) xmedi=c(1,i)+0.5d0*dxi ymedi=c(2,i)+0.5d0*dyi zmedi=c(3,i)+0.5d0*dzi xmedi=mod(xmedi,boxxsize) if (xmedi.lt.0.0d0) xmedi=xmedi+boxxsize ymedi=mod(ymedi,boxysize) if (ymedi.lt.0.0d0) ymedi=ymedi+boxysize zmedi=mod(zmedi,boxzsize) if (zmedi.lt.0.0d0) zmedi=zmedi+boxzsize num_conti=0 c write (iout,*) 'i',i,' ielstart',ielstart_vdw(i), c & ' ielend',ielend_vdw(i) c call flush(iout) do j=ielstart_vdw(i),ielend_vdw(i) if (itype(j).eq.ntyp1 .or. itype(j+1).eq.ntyp1) cycle ind=ind+1 iteli=itel(i) itelj=itel(j) if (j.eq.i+2 .and. itelj.eq.2) iteli=2 aaa=app(iteli,itelj) bbb=bpp(iteli,itelj) dxj=dc(1,j) dyj=dc(2,j) dzj=dc(3,j) dx_normj=dc_norm(1,j) dy_normj=dc_norm(2,j) dz_normj=dc_norm(3,j) xj=c(1,j)+0.5D0*dxj yj=c(2,j)+0.5D0*dyj zj=c(3,j)+0.5D0*dzj xj=mod(xj,boxxsize) if (xj.lt.0) xj=xj+boxxsize yj=mod(yj,boxysize) if (yj.lt.0) yj=yj+boxysize zj=mod(zj,boxzsize) if (zj.lt.0) zj=zj+boxzsize dist_init=(xj-xmedi)**2+(yj-ymedi)**2+(zj-zmedi)**2 xj_safe=xj yj_safe=yj zj_safe=zj isubchap=0 do xshift=-1,1 do yshift=-1,1 do zshift=-1,1 xj=xj_safe+xshift*boxxsize yj=yj_safe+yshift*boxysize zj=zj_safe+zshift*boxzsize dist_temp=(xj-xmedi)**2+(yj-ymedi)**2+(zj-zmedi)**2 if(dist_temp.lt.dist_init) then dist_init=dist_temp xj_temp=xj yj_temp=yj zj_temp=zj isubchap=1 endif enddo enddo enddo if (isubchap.eq.1) then xj=xj_temp-xmedi yj=yj_temp-ymedi zj=zj_temp-zmedi else xj=xj_safe-xmedi yj=yj_safe-ymedi zj=zj_safe-zmedi endif rij=xj*xj+yj*yj+zj*zj rrmij=1.0D0/rij rij=dsqrt(rij) c sss=sscale(rij/rpp(iteli,itelj)) c sssgrad=sscagrad(rij/rpp(iteli,itelj)) sss=sscale(rij/rpp(iteli,itelj),r_cut_respa) sssgrad=sscagrad(rij/rpp(iteli,itelj),r_cut_respa) if (sss.gt.0.0d0) then rmij=1.0D0/rij r3ij=rrmij*rmij r6ij=r3ij*r3ij ev1=aaa*r6ij*r6ij c 4/26/02 - AL scaling down 1,4 repulsive VDW interactions if (j.eq.i+2) ev1=scal_el*ev1 ev2=bbb*r6ij evdwij=ev1+ev2 if (energy_dec) then write (iout,'(a6,2i5,0pf7.3,f7.3)') 'evdw1',i,j,evdwij,sss endif evdw1=evdw1+evdwij*sss if (energy_dec) write (iout,'(a10,2i5,0pf7.3)') & 'evdw1_sum',i,j,evdw1 C C Calculate contributions to the Cartesian gradient. C facvdw=-6*rrmij*(ev1+evdwij)*sss ggg(1)=facvdw*xj+sssgrad*rmij*evdwij*xj/rpp(iteli,itelj) ggg(2)=facvdw*yj+sssgrad*rmij*evdwij*yj/rpp(iteli,itelj) ggg(3)=facvdw*zj+sssgrad*rmij*evdwij*zj/rpp(iteli,itelj) C ggg(1)=facvdw*xj C ggg(2)=facvdw*yj C ggg(3)=facvdw*zj do k=1,3 gvdwpp(k,j)=gvdwpp(k,j)+ggg(k) gvdwpp(k,i)=gvdwpp(k,i)-ggg(k) enddo endif enddo ! j enddo ! i return end C----------------------------------------------------------------------------- subroutine escp_long(evdw2,evdw2_14) C C This subroutine calculates the excluded-volume interaction energy between C peptide-group centers and side chains and its gradient in virtual-bond and C side-chain vectors. C implicit real*8 (a-h,o-z) include 'DIMENSIONS' include 'COMMON.GEO' include 'COMMON.VAR' include 'COMMON.LOCAL' include 'COMMON.CHAIN' include 'COMMON.DERIV' include 'COMMON.INTERACT' include 'COMMON.FFIELD' include 'COMMON.IOUNITS' include 'COMMON.CONTROL' include "COMMON.SPLITELE" logical lprint_short common /shortcheck/ lprint_short double precision ggg(3) integer xshift,yshift,zshift integer i,iint,j,k,iteli,itypj,subchap double precision xi,yi,zi,xj,yj,zj,rrij,sss1,sssgrad1, & fac,e1,e2,rij double precision evdw2,evdw2_14,evdwij double precision xj_safe,yj_safe,zj_safe,xj_temp,yj_temp,zj_temp, & dist_temp, dist_init double precision sscale,sscagrad if (energy_dec) write (iout,*) "escp_long:",r_cut,rlamb evdw2=0.0D0 evdw2_14=0.0d0 CD print '(a)','Enter ESCP KURWA' cd write (iout,*) 'iatscp_s=',iatscp_s,' iatscp_e=',iatscp_e c if (lprint_short) c & write (iout,*) 'ESCP_LONG iatscp_s=',iatscp_s, c & ' iatscp_e=',iatscp_e do i=iatscp_s,iatscp_e if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1) cycle iteli=itel(i) xi=0.5D0*(c(1,i)+c(1,i+1)) yi=0.5D0*(c(2,i)+c(2,i+1)) zi=0.5D0*(c(3,i)+c(3,i+1)) xi=mod(xi,boxxsize) if (xi.lt.0) xi=xi+boxxsize yi=mod(yi,boxysize) if (yi.lt.0) yi=yi+boxysize zi=mod(zi,boxzsize) if (zi.lt.0) zi=zi+boxzsize do iint=1,nscp_gr(i) do j=iscpstart(i,iint),iscpend(i,iint) itypj=iabs(itype(j)) if (itypj.eq.ntyp1) cycle C Uncomment following three lines for SC-p interactions c xj=c(1,nres+j)-xi c yj=c(2,nres+j)-yi c zj=c(3,nres+j)-zi C Uncomment following three lines for Ca-p interactions xj=c(1,j) yj=c(2,j) zj=c(3,j) c corrected by AL xj=mod(xj,boxxsize) if (xj.lt.0) xj=xj+boxxsize yj=mod(yj,boxysize) if (yj.lt.0) yj=yj+boxysize zj=mod(zj,boxzsize) if (zj.lt.0) zj=zj+boxzsize c end correction dist_init=(xj-xi)**2+(yj-yi)**2+(zj-zi)**2 xj_safe=xj yj_safe=yj zj_safe=zj subchap=0 do xshift=-1,1 do yshift=-1,1 do zshift=-1,1 xj=xj_safe+xshift*boxxsize yj=yj_safe+yshift*boxysize zj=zj_safe+zshift*boxzsize dist_temp=(xj-xi)**2+(yj-yi)**2+(zj-zi)**2 if(dist_temp.lt.dist_init) then dist_init=dist_temp xj_temp=xj yj_temp=yj zj_temp=zj subchap=1 endif enddo enddo enddo if (subchap.eq.1) then xj=xj_temp-xi yj=yj_temp-yi zj=zj_temp-zi else xj=xj_safe-xi yj=yj_safe-yi zj=zj_safe-zi endif rrij=1.0D0/(xj*xj+yj*yj+zj*zj) sss1=sscale(1.0d0/(dsqrt(rrij)),r_cut_int) if (sss1.eq.0) cycle sss=sscale(1.0d0/(dsqrt(rrij)*rscp(itypj,iteli)),r_cut_respa) sssgrad= & sscagrad(1.0d0/(dsqrt(rrij)*rscp(itypj,iteli)),r_cut_respa) sssgrad1=sscagrad(1.0d0/dsqrt(rrij),r_cut_int) if (energy_dec) write (iout,*) "rrij",1.0d0/dsqrt(rrij), & " rscp",rscp(itypj,iteli)," subchap",subchap," sss",sss 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)*sss1 endif evdwij=e1+e2 evdw2=evdw2+evdwij*(1.0d0-sss)*sss1 if (energy_dec) write (iout,'(a6,2i5,2(0pf7.3))') & 'evdw2',i,j,sss,evdwij C C Calculate contributions to the gradient in the virtual-bond and SC vectors. C fac=-(evdwij+e1)*rrij*(1.0d0-sss)*sss1 fac=fac+evdwij*dsqrt(rrij)*(-sssgrad/rscp(itypj,iteli) & +sssgrad1)/expon ggg(1)=xj*fac ggg(2)=yj*fac ggg(3)=zj*fac C Uncomment following three lines for SC-p interactions c do k=1,3 c gradx_scp(k,j)=gradx_scp(k,j)+ggg(k) c enddo C Uncomment following line for SC-p interactions c gradx_scp(k,j)=gradx_scp(k,j)+ggg(k) do k=1,3 gvdwc_scpp(k,i)=gvdwc_scpp(k,i)-ggg(k) gvdwc_scp(k,j)=gvdwc_scp(k,j)+ggg(k) enddo endif enddo enddo ! iint enddo ! i do i=1,nct do j=1,3 gvdwc_scp(j,i)=expon*gvdwc_scp(j,i) gvdwc_scpp(j,i)=expon*gvdwc_scpp(j,i) gradx_scp(j,i)=expon*gradx_scp(j,i) enddo enddo C****************************************************************************** C C N O T E !!! C C To save time the factor EXPON has been extracted from ALL components C of GVDWC and GRADX. Remember to multiply them by this factor before further C use! C C****************************************************************************** return end C----------------------------------------------------------------------------- subroutine escp_short(evdw2,evdw2_14) C C This subroutine calculates the excluded-volume interaction energy between C peptide-group centers and side chains and its gradient in virtual-bond and C side-chain vectors. C implicit none 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' include "COMMON.SPLITELE" integer xshift,yshift,zshift logical lprint_short common /shortcheck/ lprint_short integer i,iint,j,k,iteli,itypj,subchap double precision xi,yi,zi,xj,yj,zj,rrij,sss1,sssgrad1, & fac,e1,e2,rij double precision evdw2,evdw2_14,evdwij double precision xj_safe,yj_safe,zj_safe,xj_temp,yj_temp,zj_temp, & dist_temp, dist_init double precision ggg(3) double precision sscale,sscagrad evdw2=0.0D0 evdw2_14=0.0d0 cd print '(a)','Enter ESCP' c if (lprint_short) c & write (iout,*) 'ESCP_SHORT iatscp_s=',iatscp_s, c & ' iatscp_e=',iatscp_e if (energy_dec) write (iout,*) "escp_short:",r_cut_int,rlamb do i=iatscp_s,iatscp_e if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1) cycle iteli=itel(i) xi=0.5D0*(c(1,i)+c(1,i+1)) yi=0.5D0*(c(2,i)+c(2,i+1)) zi=0.5D0*(c(3,i)+c(3,i+1)) xi=mod(xi,boxxsize) if (xi.lt.0) xi=xi+boxxsize yi=mod(yi,boxysize) if (yi.lt.0) yi=yi+boxysize zi=mod(zi,boxzsize) if (zi.lt.0) zi=zi+boxzsize c if (lprint_short) c & write (iout,*) "i",i," itype",itype(i),itype(i+1), c & " nscp_gr",nscp_gr(i) do iint=1,nscp_gr(i) do j=iscpstart(i,iint),iscpend(i,iint) itypj=iabs(itype(j)) c if (lprint_short) c & write (iout,*) "j",j," itypj",itypj if (itypj.eq.ntyp1) cycle C Uncomment following three lines for SC-p interactions c xj=c(1,nres+j)-xi c yj=c(2,nres+j)-yi c zj=c(3,nres+j)-zi C Uncomment following three lines for Ca-p interactions xj=c(1,j) yj=c(2,j) zj=c(3,j) c corrected by AL xj=mod(xj,boxxsize) if (xj.lt.0) xj=xj+boxxsize yj=mod(yj,boxysize) if (yj.lt.0) yj=yj+boxysize zj=mod(zj,boxzsize) if (zj.lt.0) zj=zj+boxzsize c end correction dist_init=(xj-xi)**2+(yj-yi)**2+(zj-zi)**2 c if (lprint_short) then c write (iout,*) i,j,xi,yi,zi,xj,yj,zj c write (iout,*) "dist_init",dsqrt(dist_init) c endif xj_safe=xj yj_safe=yj zj_safe=zj subchap=0 do xshift=-1,1 do yshift=-1,1 do zshift=-1,1 xj=xj_safe+xshift*boxxsize yj=yj_safe+yshift*boxysize zj=zj_safe+zshift*boxzsize dist_temp=(xj-xi)**2+(yj-yi)**2+(zj-zi)**2 if(dist_temp.lt.dist_init) then dist_init=dist_temp xj_temp=xj yj_temp=yj zj_temp=zj subchap=1 endif enddo enddo enddo c if (lprint_short) write (iout,*) "dist_temp",dsqrt(dist_temp) if (subchap.eq.1) then xj=xj_temp-xi yj=yj_temp-yi zj=zj_temp-zi else xj=xj_safe-xi yj=yj_safe-yi zj=zj_safe-zi endif rrij=1.0D0/(xj*xj+yj*yj+zj*zj) c sss=sscale(1.0d0/(dsqrt(rrij)*rscp(itypj,iteli))) c sssgrad=sscagrad(1.0d0/(dsqrt(rrij)*rscp(itypj,iteli))) sss=sscale(1.0d0/(dsqrt(rrij)*rscp(itypj,iteli)),r_cut_respa) sssgrad=sscagrad(1.0d0/(dsqrt(rrij)*rscp(itypj,iteli)), & r_cut_respa) if (energy_dec) write (iout,*) "rrij",1.0d0/dsqrt(rrij), & " rscp",rscp(itypj,iteli)," subchap",subchap," sss",sss c if (lprint_short) write (iout,*) "rij",1.0/dsqrt(rrij), c & " subchap",subchap," sss",sss if (sss.gt.0.0d0) then fac=rrij**expon2 e1=fac*fac*aad(itypj,iteli) e2=fac*bad(itypj,iteli) if (iabs(j-i) .le. 2) then e1=scal14*e1 e2=scal14*e2 evdw2_14=evdw2_14+(e1+e2)*sss endif evdwij=e1+e2 evdw2=evdw2+evdwij*sss if (energy_dec) write (iout,'(a6,2i5,2(0pf7.3))') & 'evdw2',i,j,sss,evdwij C C Calculate contributions to the gradient in the virtual-bond and SC vectors. C fac=-(evdwij+e1)*rrij*sss fac=fac+evdwij*sssgrad*dsqrt(rrij)/rscp(itypj,iteli)/expon ggg(1)=xj*fac ggg(2)=yj*fac ggg(3)=zj*fac C Uncomment following three lines for SC-p interactions c do k=1,3 c gradx_scp(k,j)=gradx_scp(k,j)+ggg(k) c enddo C Uncomment following line for SC-p interactions c gradx_scp(k,j)=gradx_scp(k,j)+ggg(k) do k=1,3 gvdwc_scpp(k,i)=gvdwc_scpp(k,i)-ggg(k) gvdwc_scp(k,j)=gvdwc_scp(k,j)+ggg(k) enddo endif enddo enddo ! iint enddo ! i do i=1,nct do j=1,3 gvdwc_scp(j,i)=expon*gvdwc_scp(j,i) gvdwc_scpp(j,i)=expon*gvdwc_scpp(j,i) gradx_scp(j,i)=expon*gradx_scp(j,i) enddo enddo C****************************************************************************** C C N O T E !!! C C To save time the factor EXPON has been extracted from ALL components C of GVDWC and GRADX. Remember to multiply them by this factor before further C use! C C****************************************************************************** return end