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