C Calculate the virtual-bond torsional energy.
C
cd print *,'nterm=',nterm
+C print *,"tor",tor_mode
if (wtor.gt.0) then
if ((tor_mode.eq.0).or.(tor_mode.eq.2)) then
call etor(etors,edihcnstr)
C etor kcc is Kubo cumulant clustered rigorous attemp to derive the
C energy function
if ((tor_mode.eq.1).or.(tor_mode.eq.2)) then
- call etor(etors,edihcnstr)
+ call etor_kcc(etors,edihcnstr)
endif
else
etors=0
C Set lprn=.true. for debugging
lprn=.false.
c lprn=.true.
+C print *,"wchodze kcc"
+ if (tor_mode.ne.2) then
etors=0.0D0
+ endif
do i=iphi_start,iphi_end
C ANY TWO ARE DUMMY ATOMS in row CYCLE
c if (((itype(i-3).eq.ntyp1).and.(itype(i-2).eq.ntyp1)).or.
cosphi=dcos(j*phii)
sinphi=dsin(j*phii)
sint1t2n=sint1t2**j
- sumnonchebyshev=sumnonchebyshev+
+ sumnonchebyshev=
& sint1t2n*(v1ij*cosphi+v2ij*sinphi)
actval=sint1t2n*(v1ij*cosphi+v2ij*sinphi)
C etors=etors+sint1t2n*(v1ij*cosphi+v2ij*sinphi)
C if (energy_dec) etors_ii=etors_ii+
C & v1ij*cosphi+v2ij*sinphi
C glocig is the gradient local i site in gamma
- glocig=glocig+j*(v2ij*cosphi-v1ij*sinphi)*sint1t2n
+ glocig=j*(v2ij*cosphi-v1ij*sinphi)*sint1t2n
C now gradient over theta_1
- glocit1=glocit1+actval/sinthet1*j*costhet1
- glocit2=glocit2+actval/sinthet2*j*costhet2
- enddo
+ glocit1=actval/sinthet1*j*costhet1
+ glocit2=actval/sinthet2*j*costhet2
C now the Czebyshev polinominal sum
- do j=1,nterm_kcc_Tb(itori,itori1)
- thybt1(j)=v1_chyb(j,itori,itori1)
- thybt2(j)=v2_chyb(j,itori,itori1)
+ do k=1,nterm_kcc_Tb(itori,itori1)
+ thybt1(k)=v1_chyb(k,j,itori,itori1)
+ thybt2(k)=v2_chyb(k,j,itori,itori1)
+C thybt1(k)=0.0
+C thybt2(k)=0.0
enddo
sumth1thyb=tschebyshev
- & (1,nterm_kcc_Tb(itori,itori1),thybt1(1),dcos(theta12))
+ & (1,nterm_kcc_Tb(itori,itori1),thybt1(1),dcos(theti12)**2)
gradthybt1=gradtschebyshev
- & (0,nterm_kcc_Tb(itori,itori1)-1,thybt1(1),dcos(theta12))
- & *(nterm_kcc_Tb(itori,itori1))*0.5*dsin(theta12)
+ & (0,nterm_kcc_Tb(itori,itori1)-1,thybt1(1),
+ & dcos(theti12)**2)
+ & *dcos(theti12)*(-dsin(theti12))
sumth2thyb=tschebyshev
- & (1,nterm_kcc_Tb(itori,itori1),thybt2(1),dcos(theta22))
+ & (1,nterm_kcc_Tb(itori,itori1),thybt2(1),dcos(theti22)**2)
gradthybt2=gradtschebyshev
- & (0,nterm_kcc_Tb(itori,itori1)-1,thybt2(1),dcos(theta22))
- & *(nterm_kcc_Tb(itori,itori1))*0.5*dsin(theta22)
+ & (0,nterm_kcc_Tb(itori,itori1)-1,thybt2(1),
+ & dcos(theti22)**2)
+ & *dcos(theti22)*(-dsin(theti22))
+C print *, sumth1thyb, gradthybt1, sumth2thyb, gradthybt2,
+C & gradtschebyshev
+C & (0,nterm_kcc_Tb(itori,itori1)-1,thybt2(1),
+C & dcos(theti22)**2),
+C & dsin(theti22)
C now overal sumation
etors=etors+(1.0d0+sumth1thyb+sumth2thyb)*sumnonchebyshev
+C print *,"sumnon", sumnonchebyshev,sumth1thyb+sumth2thyb
C derivative over gamma
gloc(i-3,icg)=gloc(i-3,icg)+wtor*glocig
& *(1.0d0+sumth1thyb+sumth2thyb)
C derivative over theta1
- gloc(nphi+i-3,icg)=gloc(nphi+i-3,icg)+wang*
+ gloc(nphi+i-3,icg)=gloc(nphi+i-3,icg)+wtor*
& (glocit1*(1.0d0+sumth1thyb+sumth2thyb)+
& sumnonchebyshev*gradthybt1)
C now derivative over theta2
- gloc(nphi+i-2,icg)=gloc(nphi+i-2,icg)+wang*
+ gloc(nphi+i-2,icg)=gloc(nphi+i-2,icg)+wtor*
& (glocit2*(1.0d0+sumth1thyb+sumth2thyb)+
& sumnonchebyshev*gradthybt2)
enddo
-
+ enddo
C gloc(i-3,icg)=gloc(i-3,icg)+wtor*gloci
! 6/20/98 - dihedral angle constraints
return
end
+C The rigorous attempt to derive energy function
+ subroutine ebend_kcc(etheta,ethetacnstr)
-
+ implicit real*8 (a-h,o-z)
+ include 'DIMENSIONS'
+ include 'COMMON.VAR'
+ include 'COMMON.GEO'
+ include 'COMMON.LOCAL'
+ include 'COMMON.TORSION'
+ include 'COMMON.INTERACT'
+ include 'COMMON.DERIV'
+ include 'COMMON.CHAIN'
+ include 'COMMON.NAMES'
+ include 'COMMON.IOUNITS'
+ include 'COMMON.FFIELD'
+ include 'COMMON.TORCNSTR'
+ include 'COMMON.CONTROL'
+ logical lprn
+ double precision thybt1(maxtermkcc)
+C Set lprn=.true. for debugging
+ lprn=.false.
+c lprn=.true.
+C print *,"wchodze kcc"
+ if (tormode.ne.2) etheta=0.0D0
+ do i=ithet_start,ithet_end
+c print *,i,itype(i-1),itype(i),itype(i-2)
+ if ((itype(i-1).eq.ntyp1).or.itype(i-2).eq.ntyp1
+ & .or.itype(i).eq.ntyp1) cycle
+ iti=itortyp_kcc(itype(i-1))
+ sinthet=dsin(theta(i)/2.0d0)
+ costhet=dcos(theta(i)/2.0d0)
+ do j=1,nbend_kcc_Tb(iti)
+ thybt1(j)=v1bend_chyb(j,iti)
+ enddo
+ sumth1thyb=tschebyshev
+ & (1,nbend_kcc_Tb(iti),thybt1(1),costhet)
+ ihelp=nbend_kcc_Tb(iti)-1
+ gradthybt1=gradtschebyshev
+ & (0,ihelp,thybt1(1),costhet)
+ etheta=etheta+sumth1thyb
+C print *,sumth1thyb,gradthybt1,sinthet*(-0.5d0)
+ gloc(nphi+i-2,icg)=gloc(nphi+i-2,icg)+wang*
+ & gradthybt1*sinthet*(-0.5d0)
+ enddo
+ if (tormode.ne.2) then
+ ethetacnstr=0.0d0
+C print *,ithetaconstr_start,ithetaconstr_end,"TU"
+ do i=ithetaconstr_start,ithetaconstr_end
+ itheta=itheta_constr(i)
+ thetiii=theta(itheta)
+ difi=pinorm(thetiii-theta_constr0(i))
+ if (difi.gt.theta_drange(i)) then
+ difi=difi-theta_drange(i)
+ ethetacnstr=ethetacnstr+0.25d0*for_thet_constr(i)*difi**4
+ gloc(itheta+nphi-2,icg)=gloc(itheta+nphi-2,icg)
+ & +for_thet_constr(i)*difi**3
+ else if (difi.lt.-drange(i)) then
+ difi=difi+drange(i)
+ ethetacnstr=ethetacnstr+0.25d0*for_thet_constr(i)*difi**4
+ gloc(itheta+nphi-2,icg)=gloc(itheta+nphi-2,icg)
+ & +for_thet_constr(i)*difi**3
+ else
+ difi=0.0
+ endif
+ if (energy_dec) then
+ write (iout,'(a6,2i5,4f8.3,2e14.5)') "ethetc",
+ & i,itheta,rad2deg*thetiii,
+ & rad2deg*theta_constr0(i), rad2deg*theta_drange(i),
+ & rad2deg*difi,0.25d0*for_thet_constr(i)*difi**4,
+ & gloc(itheta+nphi-2,icg)
+ endif
+ enddo
+ endif
+ return
+ end
c------------------------------------------------------------------------------
subroutine eback_sc_corr(esccor)
c 7/21/2007 Correlations between the backbone-local and side-chain-local
implicit none
include "DIMENSIONS"
integer i,m,n
- double precision x(n),y,yy(0:maxvar),aux
+ double precision x(n+1),y,yy(0:maxvar),aux
c Tschebyshev polynomial. Note that the first term is omitted
c m=0: the constant term is included
c m=1: the constant term is not included
yy(0)=1.0d0
yy(1)=2.0d0*y
do i=2,n
- yy(i)=2*yy(1)*yy(i-1)-yy(i-2)
+ yy(i)=2*y*yy(i-1)-yy(i-2)
enddo
aux=0.0d0
do i=m,n
- aux=aux+x(i)*yy(i)
+ aux=aux+x(i+1)*yy(i)*(i+1)
+C print *, x(i+1),yy(i),i
enddo
gradtschebyshev=aux
return
projects / unres.git / blobdiff