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