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) implicit none 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----------------------------------------------------------------------- double precision function sscagrad(r) implicit none double precision r,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 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' c 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=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 sss=sscale(dsqrt(rij)/sigma(itypi,itypj)) if (sss.lt.1.0d0) then rrij=1.0D0/rij eps0ij=eps(itypi,itypj) fac=rrij**expon2 e1=fac*fac*aa e2=fac*bb evdwij=e1+e2 evdw=evdw+(1.0d0-sss)*evdwij C C Calculate the components of the gradient in DC and X C fac=-rrij*(e1+evdwij)*(1.0d0-sss) gg(1)=xj*fac gg(2)=yj*fac gg(3)=zj*fac do k=1,3 gvdwx(k,i)=gvdwx(k,i)-gg(k) gvdwx(k,j)=gvdwx(k,j)+gg(k) gvdwc(k,i)=gvdwc(k,i)-gg(k) gvdwc(k,j)=gvdwc(k,j)+gg(k) enddo endif enddo ! j enddo ! iint enddo ! i do i=1,nct do j=1,3 gvdwc(j,i)=expon*gvdwc(j,i) gvdwx(j,i)=expon*gvdwx(j,i) enddo enddo C****************************************************************************** C C N O T E !!! C C To save time, the factor of EXPON has been extracted from ALL components C of GVDWC and GRADX. Remember to multiply them by this factor before further C use! C C****************************************************************************** return end C----------------------------------------------------------------------- subroutine elj_short(evdw) C C This subroutine calculates the interaction energy of nonbonded side chains C assuming the LJ potential of interaction. C implicit real*8 (a-h,o-z) include 'DIMENSIONS' parameter (accur=1.0d-10) include 'COMMON.GEO' include 'COMMON.VAR' include 'COMMON.LOCAL' include 'COMMON.CHAIN' include 'COMMON.DERIV' include 'COMMON.INTERACT' include 'COMMON.TORSION' include 'COMMON.SBRIDGE' include 'COMMON.NAMES' include 'COMMON.IOUNITS' c 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=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 sss=sscale(dsqrt(rij)/sigma(itypi,itypj)) if (sss.gt.0.0d0) then rrij=1.0D0/rij eps0ij=eps(itypi,itypj) fac=rrij**expon2 e1=fac*fac*aa 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 gg(1)=xj*fac gg(2)=yj*fac gg(3)=zj*fac do k=1,3 gvdwx(k,i)=gvdwx(k,i)-gg(k) gvdwx(k,j)=gvdwx(k,j)+gg(k) gvdwc(k,i)=gvdwc(k,i)-gg(k) gvdwc(k,j)=gvdwc(k,j)+gg(k) enddo endif enddo ! j enddo ! iint enddo ! i do i=1,nct do j=1,3 gvdwc(j,i)=expon*gvdwc(j,i) gvdwx(j,i)=expon*gvdwx(j,i) enddo enddo C****************************************************************************** C C N O T E !!! C C To save time, the factor of EXPON has been extracted from ALL components C of GVDWC and GRADX. Remember to multiply them by this factor before further C use! C C****************************************************************************** return end C----------------------------------------------------------------------------- subroutine eljk_long(evdw) C C This subroutine calculates the interaction energy of nonbonded side chains C assuming the LJK potential of interaction. C implicit real*8 (a-h,o-z) include 'DIMENSIONS' include 'COMMON.GEO' include 'COMMON.VAR' include 'COMMON.LOCAL' include 'COMMON.CHAIN' include 'COMMON.DERIV' include 'COMMON.INTERACT' include 'COMMON.IOUNITS' include 'COMMON.NAMES' dimension gg(3) logical scheck c print *,'Entering ELJK nnt=',nnt,' nct=',nct,' expon=',expon evdw=0.0D0 do i=iatsc_s,iatsc_e itypi=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)) 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 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)*evdwij C C Calculate the components of the gradient in DC and X C fac=-2.0D0*rrij*e_augm-r_inv_ij*r_shift_inv*(e1+e1+e2) fac=fac*(1.0d0-sss) gg(1)=xj*fac gg(2)=yj*fac gg(3)=zj*fac do k=1,3 gvdwx(k,i)=gvdwx(k,i)-gg(k) gvdwx(k,j)=gvdwx(k,j)+gg(k) gvdwc(k,i)=gvdwc(k,i)-gg(k) gvdwc(k,j)=gvdwc(k,j)+gg(k) enddo endif enddo ! j enddo ! iint enddo ! i do i=1,nct do j=1,3 gvdwc(j,i)=expon*gvdwc(j,i) gvdwx(j,i)=expon*gvdwx(j,i) enddo enddo return end C----------------------------------------------------------------------------- subroutine eljk_short(evdw) C C This subroutine calculates the interaction energy of nonbonded side chains C assuming the LJK potential of interaction. C implicit real*8 (a-h,o-z) include 'DIMENSIONS' include 'COMMON.GEO' include 'COMMON.VAR' include 'COMMON.LOCAL' include 'COMMON.CHAIN' include 'COMMON.DERIV' include 'COMMON.INTERACT' include 'COMMON.IOUNITS' include 'COMMON.NAMES' dimension gg(3) logical scheck c print *,'Entering ELJK nnt=',nnt,' nct=',nct,' expon=',expon evdw=0.0D0 do i=iatsc_s,iatsc_e itypi=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)) 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) fac=fac*sss gg(1)=xj*fac gg(2)=yj*fac gg(3)=zj*fac do k=1,3 gvdwx(k,i)=gvdwx(k,i)-gg(k) gvdwx(k,j)=gvdwx(k,j)+gg(k) gvdwc(k,i)=gvdwc(k,i)-gg(k) gvdwc(k,j)=gvdwc(k,j)+gg(k) enddo endif enddo ! j enddo ! iint enddo ! i do i=1,nct do j=1,3 gvdwc(j,i)=expon*gvdwc(j,i) gvdwx(j,i)=expon*gvdwx(j,i) enddo enddo return end C----------------------------------------------------------------------------- subroutine ebp_long(evdw) C C This subroutine calculates the interaction energy of nonbonded side chains C assuming the Berne-Pechukas potential of interaction. C implicit real*8 (a-h,o-z) include 'DIMENSIONS' include 'COMMON.GEO' include 'COMMON.VAR' include 'COMMON.LOCAL' include 'COMMON.CHAIN' include 'COMMON.DERIV' include 'COMMON.NAMES' include 'COMMON.INTERACT' include 'COMMON.IOUNITS' include 'COMMON.CALC' common /srutu/ icall c double precision rrsave(maxdim) logical lprn evdw=0.0D0 c print *,'Entering EBP nnt=',nnt,' nct=',nct,' expon=',expon evdw=0.0D0 c if (icall.eq.0) then c lprn=.true. c else lprn=.false. c endif ind=0 do i=iatsc_s,iatsc_e itypi=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))) 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 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) 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=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=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))) 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=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 integer xshift,yshift,zshift 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- & ((positi-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-positi)/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))) sssgrad=sscagrad((1.0d0/rij)/sigmaii(itypi,itypj)) if (sss.lt.1.0d0) then C Calculate angle-dependent terms of energy and contributions to their C derivatives. call sc_angular sigsq=1.0D0/sigsq sig=sig0ij*dsqrt(sigsq) rij_shift=1.0D0/rij-sig+sig0ij c for diagnostics; uncomment c rij_shift=1.2*sig0ij C I hate to put IF's in the loops, but here don't have another choice!!!! if (rij_shift.le.0.0D0) then evdw=1.0D20 cd write (iout,'(2(a3,i3,2x),17(0pf7.3))') cd & restyp(itypi),i,restyp(itypj),j, cd & rij_shift,1.0D0/rij,sig,sig0ij,sigsq,1-dsqrt(sigsq) return endif sigder=-sig*sigsq c--------------------------------------------------------------- rij_shift=1.0D0/rij_shift fac=rij_shift**expon e1=fac*fac*aa 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) 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=rij*fac fac=fac+evdwij/(1.0-sss)*(-sssgrad)/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(1.0d0-sss) endif enddo ! j enddo ! iint enddo ! i c write (iout,*) "Number of loop steps in EGB:",ind cccc energy_dec=.false. return end C----------------------------------------------------------------------------- subroutine egb_short(evdw) C C This subroutine calculates the interaction energy of nonbonded side chains C assuming the Gay-Berne potential of interaction. C implicit real*8 (a-h,o-z) include 'DIMENSIONS' include 'COMMON.GEO' include 'COMMON.VAR' include 'COMMON.LOCAL' include 'COMMON.CHAIN' include 'COMMON.DERIV' include 'COMMON.NAMES' include 'COMMON.INTERACT' include 'COMMON.IOUNITS' include 'COMMON.CALC' include 'COMMON.CONTROL' logical lprn integer xshift,yshift,zshift 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- & ((positi-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-positi)/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))) sssgrad=sscagrad((1.0d0/rij)/sigmaii(itypi,itypj)) if (sss.gt.0.0d0) then C Calculate angle-dependent terms of energy and contributions to their C derivatives. call sc_angular sigsq=1.0D0/sigsq sig=sig0ij*dsqrt(sigsq) rij_shift=1.0D0/rij-sig+sig0ij c for diagnostics; uncomment c rij_shift=1.2*sig0ij C I hate to put IF's in the loops, but here don't have another choice!!!! if (rij_shift.le.0.0D0) then evdw=1.0D20 cd write (iout,'(2(a3,i3,2x),17(0pf7.3))') cd & restyp(itypi),i,restyp(itypj),j, cd & rij_shift,1.0D0/rij,sig,sig0ij,sigsq,1-dsqrt(sigsq) return endif sigder=-sig*sigsq c--------------------------------------------------------------- rij_shift=1.0D0/rij_shift fac=rij_shift**expon e1=fac*fac*aa 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=rij*fac fac=fac+evdwij/sss*sssgrad/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' 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=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))) 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 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)*(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 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=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))) 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 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)+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' 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 real*8 (a-h,o-z) include 'DIMENSIONS' #ifdef MPI include "mpif.h" #endif include 'COMMON.CONTROL' include 'COMMON.IOUNITS' include 'COMMON.GEO' include 'COMMON.VAR' include 'COMMON.LOCAL' include 'COMMON.CHAIN' include 'COMMON.DERIV' include 'COMMON.INTERACT' #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' integer xshift,yshift,zshift 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),gmuij1(4),gmuji1(4), & gmuij2(4),gmuji2(4) common /locel/ a_temp,agg,aggi,aggi1,aggj,aggj1,a22,a23,a32,a33, & dxi,dyi,dzi,dx_normi,dy_normi,dz_normi,xmedi,ymedi,zmedi, & num_conti,j1,j2 c 4/26/02 - AL scaling factor for 1,4 repulsive VDW interactions #ifdef MOMENT double precision scal_el /1.0d0/ #else double precision scal_el /0.5d0/ #endif C 12/13/98 C 13-go grudnia roku pamietnego... double precision unmat(3,3) /1.0d0,0.0d0,0.0d0, & 0.0d0,1.0d0,0.0d0, & 0.0d0,0.0d0,1.0d0/ c time00=MPI_Wtime() cd write (iout,*) "eelecij",i,j 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 sss=sscale(rij/rpp(iteli,itelj)) sssgrad=sscagrad(rij/rpp(iteli,itelj)) r3ij=rrmij*rmij r6ij=r3ij*r3ij cosa=dx_normi*dx_normj+dy_normi*dy_normj+dz_normi*dz_normj cosb=(xj*dx_normi+yj*dy_normi+zj*dz_normi)*rmij cosg=(xj*dx_normj+yj*dy_normj+zj*dz_normj)*rmij fac=cosa-3.0D0*cosb*cosg ev1=aaa*r6ij*r6ij c 4/26/02 - AL scaling down 1,4 repulsive VDW interactions if (j.eq.i+2) ev1=scal_el*ev1 ev2=bbb*r6ij fac3=ael6i*r6ij fac4=ael3i*r3ij evdwij=ev1+ev2 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 evdw1=evdw1+evdwij*(1.0d0-sss) cd write(iout,'(2(2i3,2x),7(1pd12.4)/2(3(1pd12.4),5x)/)') cd & iteli,i,itelj,j,aaa,bbb,ael6i,ael3i, cd & 1.0D0/dsqrt(rrmij),evdwij,eesij, cd & xmedi,ymedi,zmedi,xj,yj,zj if (energy_dec) then write (iout,'(a6,2i5,0pf7.3,f7.3)') 'evdw1',i,j,evdwij,sss write (iout,'(a6,2i5,0pf7.3)') 'ees',i,j,eesij endif C C Calculate contributions to the Cartesian gradient. C #ifdef SPLITELE facvdw=-6*rrmij*(ev1+evdwij)*(1.0d0-sss) facel=-3*rrmij*(el1+eesij) fac1=fac erij(1)=xj*rmij erij(2)=yj*rmij erij(3)=zj*rmij * * Radial derivatives. First process both termini of the fragment (i,j) * ggg(1)=facel*xj ggg(2)=facel*yj ggg(3)=facel*zj 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/fac_shield(i) & *2.0 gshieldx(k,iresshield)=gshieldx(k,iresshield)+ & rlocshield & +grad_shield_loc(k,ilist,i)*eesij/fac_shield(i)*2.0 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 gshieldx(k,iresshield)=gshieldx(k,iresshield)+ & rlocshield & +grad_shield_loc(k,ilist,j)*eesij/fac_shield(j)*2.0 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 gshieldc(k,j)=gshieldc(k,j)+ & grad_shield(k,j)*eesij/fac_shield(j)*2.0 gshieldc(k,i-1)=gshieldc(k,i-1)+ & grad_shield(k,i)*eesij/fac_shield(i)*2.0 gshieldc(k,j-1)=gshieldc(k,j-1)+ & grad_shield(k,j)*eesij/fac_shield(j)*2.0 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 * * Loop over residues i+1 thru j-1. * cgrad do k=i+1,j-1 cgrad do l=1,3 cgrad gelc(l,k)=gelc(l,k)+ggg(l) cgrad enddo cgrad enddo ggg(1)=facvdw*xj-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 do k=1,3 c ghalf=0.5D0*ggg(k) c gvdwpp(k,i)=gvdwpp(k,i)+ghalf c gvdwpp(k,j)=gvdwpp(k,j)+ghalf c enddo c 9/28/08 AL Gradient compotents will be summed only at the end do k=1,3 gvdwpp(k,j)=gvdwpp(k,j)+ggg(k) gvdwpp(k,i)=gvdwpp(k,i)-ggg(k) enddo * * Loop over residues i+1 thru j-1. * cgrad do k=i+1,j-1 cgrad do l=1,3 cgrad gvdwpp(l,k)=gvdwpp(l,k)+ggg(l) cgrad enddo cgrad enddo #else facvdw=ev1+evdwij*(1.0d0-sss) facel=el1+eesij fac1=fac fac=-3*rrmij*(facvdw+facvdw+facel) erij(1)=xj*rmij erij(2)=yj*rmij erij(3)=zj*rmij * * Radial derivatives. First process both termini of the fragment (i,j) * ggg(1)=fac*xj ggg(2)=fac*yj ggg(3)=fac*zj c do k=1,3 c ghalf=0.5D0*ggg(k) c gelc(k,i)=gelc(k,i)+ghalf c gelc(k,j)=gelc(k,j)+ghalf c enddo c 9/28/08 AL Gradient compotents will be summed only at the end do k=1,3 gelc_long(k,j)=gelc(k,j)+ggg(k) gelc_long(k,i)=gelc(k,i)-ggg(k) enddo * * Loop over residues i+1 thru j-1. * cgrad do k=i+1,j-1 cgrad do l=1,3 cgrad gelc(l,k)=gelc(l,k)+ggg(l) cgrad enddo cgrad enddo c 9/28/08 AL Gradient compotents will be summed only at the end C ggg(1)=facvdw*xj C ggg(2)=facvdw*yj C ggg(3)=facvdw*zj 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) 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))* & fac_shield(i)**2*fac_shield(j)**2 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)) & *fac_shield(i)**2*fac_shield(j)**2 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)) & *fac_shield(i)**2*fac_shield(j)**2 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) 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) 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) 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) 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) #endif cC 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)) & *fac_shield(i)*fac_shield(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)) & *fac_shield(i)*fac_shield(j) C Derivatives of eello in DC(i+1) thru DC(j-1) or DC(nres-2) do l=1,3 ggg(l)=(agg(l,1)*muij(1)+ & agg(l,2)*muij(2)+agg(l,3)*muij(3)+agg(l,4)*muij(4)) & *fac_shield(i)*fac_shield(j) gel_loc_long(l,j)=gel_loc_long(l,j)+ggg(l) gel_loc_long(l,i)=gel_loc_long(l,i)-ggg(l) cgrad ghalf=0.5d0*ggg(l) cgrad gel_loc(l,i)=gel_loc(l,i)+ghalf cgrad gel_loc(l,j)=gel_loc(l,j)+ghalf enddo cgrad do k=i+1,j2 cgrad do l=1,3 cgrad gel_loc(l,k)=gel_loc(l,k)+ggg(l) cgrad enddo cgrad enddo C Remaining derivatives of eello do l=1,3 gel_loc(l,i)=gel_loc(l,i)+(aggi(l,1)*muij(1)+ & aggi(l,2)*muij(2)+aggi(l,3)*muij(3)+aggi(l,4)*muij(4)) & *fac_shield(i)*fac_shield(j) 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) 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) 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) 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) ees0m(num_conti,i)=0.5D0*fac3*(ees0pij-ees0mij) & *fac_shield(i)*fac_shield(j) C Diagnostics. Comment out or remove after debugging! c ees0p(num_conti,i)=0.5D0*fac3*ees0pij c ees0m(num_conti,i)=0.5D0*fac3*ees0mij c ees0m(num_conti,i)=0.0D0 C End diagnostics. c write (iout,*) 'i=',i,' j=',j,' rij=',rij,' r0ij=',r0ij, c & ' ees0ij=',ees0p(num_conti,i),ees0m(num_conti,i),' fcont=',fcont C Angular derivatives of the contact function ees0pij1=fac3/ees0pij ees0mij1=fac3/ees0mij fac3p=-3.0D0*fac3*rrmij ees0pijp=0.5D0*fac3p*(ees0pij+ees0mij) ees0mijp=0.5D0*fac3p*(ees0pij-ees0mij) c ees0mij1=0.0D0 ecosa1= ees0pij1*( 1.0D0+0.5D0*wij) ecosb1=-1.5D0*ees0pij1*(wij*cosg+cosbg1) ecosg1=-1.5D0*ees0pij1*(wij*cosb+cosbg1) ecosa2= ees0mij1*(-1.0D0+0.5D0*wij) ecosb2=-1.5D0*ees0mij1*(wij*cosg+cosbg2) ecosg2=-1.5D0*ees0mij1*(wij*cosb-cosbg2) ecosap=ecosa1+ecosa2 ecosbp=ecosb1+ecosb2 ecosgp=ecosg1+ecosg2 ecosam=ecosa1-ecosa2 ecosbm=ecosb1-ecosb2 ecosgm=ecosg1-ecosg2 C Diagnostics c ecosap=ecosa1 c ecosbp=ecosb1 c ecosgp=ecosg1 c ecosam=0.0D0 c ecosbm=0.0D0 c ecosgm=0.0D0 C End diagnostics facont_hb(num_conti,i)=fcont fprimcont=fprimcont/rij cd facont_hb(num_conti,i)=1.0D0 C Following line is for diagnostics. cd fprimcont=0.0D0 do k=1,3 dcosb(k)=rmij*(dc_norm(k,i)-erij(k)*cosb) dcosg(k)=rmij*(dc_norm(k,j)-erij(k)*cosg) enddo do k=1,3 gggp(k)=ecosbp*dcosb(k)+ecosgp*dcosg(k) gggm(k)=ecosbm*dcosb(k)+ecosgm*dcosg(k) enddo gggp(1)=gggp(1)+ees0pijp*xj gggp(2)=gggp(2)+ees0pijp*yj gggp(3)=gggp(3)+ees0pijp*zj gggm(1)=gggm(1)+ees0mijp*xj gggm(2)=gggm(2)+ees0mijp*yj gggm(3)=gggm(3)+ees0mijp*zj C Derivatives due to the contact function gacont_hbr(1,num_conti,i)=fprimcont*xj gacont_hbr(2,num_conti,i)=fprimcont*yj gacont_hbr(3,num_conti,i)=fprimcont*zj do k=1,3 c c 10/24/08 cgrad and ! comments indicate the parts of the code removed c following the change of gradient-summation algorithm. c cgrad ghalfp=0.5D0*gggp(k) cgrad ghalfm=0.5D0*gggm(k) gacontp_hb1(k,num_conti,i)=!ghalfp & +(ecosap*(dc_norm(k,j)-cosa*dc_norm(k,i)) & + ecosbp*(erij(k)-cosb*dc_norm(k,i)))*vbld_inv(i+1) & *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) & *fac_shield(i)*fac_shield(j) gacontp_hb3(k,num_conti,i)=gggp(k) & *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) & *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) & *fac_shield(i)*fac_shield(j) gacontm_hb3(k,num_conti,i)=gggm(k) & *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 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' c include 'COMMON.CONTACTS' include 'COMMON.TORSION' include 'COMMON.VECTORS' include 'COMMON.FFIELD' dimension 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" 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) sssgrad=sscagrad(rij) 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' logical lprint_short common /shortcheck/ lprint_short dimension ggg(3) integer xshift,yshift,zshift 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) sss=sscale(1.0d0/(dsqrt(rrij)*rscp(itypj,iteli))) sssgrad=sscagrad(1.0d0/(dsqrt(rrij)*rscp(itypj,iteli))) 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) endif evdwij=e1+e2 evdw2=evdw2+evdwij*(1.0d0-sss) if (energy_dec) write (iout,'(a6,2i5,2(0pf7.3))') & 'evdw2',i,j,sss,evdwij C C Calculate contributions to the gradient in the virtual-bond and SC vectors. C fac=-(evdwij+e1)*rrij*(1.0d0-sss) fac=fac-(evdwij)*sssgrad*dsqrt(rrij)/rscp(itypj,iteli) ggg(1)=xj*fac ggg(2)=yj*fac ggg(3)=zj*fac C Uncomment following three lines for SC-p interactions c do k=1,3 c gradx_scp(k,j)=gradx_scp(k,j)+ggg(k) c enddo C Uncomment following line for SC-p interactions c gradx_scp(k,j)=gradx_scp(k,j)+ggg(k) do k=1,3 gvdwc_scpp(k,i)=gvdwc_scpp(k,i)-ggg(k) gvdwc_scp(k,j)=gvdwc_scp(k,j)+ggg(k) enddo endif enddo enddo ! iint enddo ! i do i=1,nct do j=1,3 gvdwc_scp(j,i)=expon*gvdwc_scp(j,i) gvdwc_scpp(j,i)=expon*gvdwc_scpp(j,i) gradx_scp(j,i)=expon*gradx_scp(j,i) enddo enddo C****************************************************************************** C C N O T E !!! C C To save time the factor EXPON has been extracted from ALL components C of GVDWC and GRADX. Remember to multiply them by this factor before further C use! C C****************************************************************************** return end C----------------------------------------------------------------------------- subroutine escp_short(evdw2,evdw2_14) C C This subroutine calculates the excluded-volume interaction energy between C peptide-group centers and side chains and its gradient in virtual-bond and C side-chain vectors. C implicit real*8 (a-h,o-z) include 'DIMENSIONS' include 'COMMON.GEO' include 'COMMON.VAR' include 'COMMON.LOCAL' include 'COMMON.CHAIN' include 'COMMON.DERIV' include 'COMMON.INTERACT' include 'COMMON.FFIELD' include 'COMMON.IOUNITS' include 'COMMON.CONTROL' integer xshift,yshift,zshift logical lprint_short common /shortcheck/ lprint_short dimension ggg(3) 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,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))) sssgrad=sscagrad(1.0d0/(dsqrt(rrij))) 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) 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