C Calculate the virtual-bond-angle energy.
C
if (wang.gt.0d0) then
+ if ((tor_mode.eq.0).or.(tor_mode.eq.2)) then
call ebend(ebe,ethetacnstr)
+ endif
+C ebend 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 ebend_kcc(ebe,ethetacnstr)
+ endif
else
ebe=0
ethetacnstr=0
C
cd print *,'nterm=',nterm
if (wtor.gt.0) then
+ if ((tor_mode.eq.0).or.(tor_mode.eq.2)) then
call etor(etors,edihcnstr)
+ endif
+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)
+ endif
else
etors=0
edihcnstr=0
C
C 6/23/01 Calculate double-torsional energy
C
- if (wtor_d.gt.0) then
+ if ((wtor_d.gt.0).and.((tor_mode.eq.0).or.(tor_mode.eq.2))) then
call etor_d(etors_d)
else
etors_d=0
return
end
#endif
+C----------------------------------------------------------------------------------
+C The rigorous attempt to derive energy function
+ subroutine etor_kcc(etors,edihcnstr)
+ 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),thybt2(maxtermkcc)
+C Set lprn=.true. for debugging
+ lprn=.false.
+c lprn=.true.
+ etors=0.0D0
+ 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.
+c & ((itype(i-2).eq.ntyp1).and.(itype(i-1).eq.ntyp1)) .or.
+c & ((itype(i-1).eq.ntyp1).and.(itype(i).eq.ntyp1))) cycle
+ if (itype(i-2).eq.ntyp1.or. itype(i-1).eq.ntyp1
+ & .or. itype(i).eq.ntyp1 .or. itype(i-3).eq.ntyp1) cycle
+ itori=itortyp_kcc(itype(i-2))
+ itori1=itortyp_kcc(itype(i-1))
+ phii=phi(i)
+ glocig=0.0D0
+ glocit1=0.0d0
+ glocit2=0.0d0
+ sumnonchebyshev=0.0d0
+ sumchebyshev=0.0d0
+C to avoid multiple devision by 2
+ theti22=0.5d0*theta(i)
+C theta 12 is the theta_1 /2
+C theta 22 is theta_2 /2
+ theti12=0.5d0*theta(i-1)
+C and appropriate sinus function
+ sinthet2=dsin(theta(i))
+ sinthet1=dsin(theta(i-1))
+ costhet1=dcos(theta(i-1))
+ costhet2=dcos(theta(i))
+C to speed up lets store its mutliplication
+ sint1t2=sinthet2*sinthet1
+C \sum_{i=1}^n (sin(theta_1) * sin(theta_2))^n * (c_n* cos(n*gamma)
+C +d_n*sin(n*gamma)) *
+C \sum_{i=1}^m (1+a_m*Tb_m(cos(theta_1 /2))+b_m*Tb_m(cos(theta_2 /2)))
+C we have two sum 1) Non-Chebyshev which is with n and gamma
+ do j=1,nterm_kcc(itori,itori1)
+
+ v1ij=v1_kcc(j,itori,itori1)
+ v2ij=v2_kcc(j,itori,itori1)
+C v1ij is c_n and d_n in euation above
+ cosphi=dcos(j*phii)
+ sinphi=dsin(j*phii)
+ sint1t2n=sint1t2**j
+ 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
+C now gradient over theta_1
+ glocit1=glocit1+actval/sinthet1*j*costhet1
+ glocit2=glocit2+actval/sinthet2*j*costhet2
+ enddo
+
+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)
+ enddo
+ sumth1thyb=tschebyshev
+ & (1,nterm_kcc_Tb(itori,itori1),thybt1(1),dcos(theta12))
+ gradthybt1=gradtschebyshev
+ & (0,nterm_kcc_Tb(itori,itori1)-1,thybt1(1),dcos(theta12))
+ & *(nterm_kcc_Tb(itori,itori1))*0.5*dsin(theta12)
+ sumth2thyb=tschebyshev
+ & (1,nterm_kcc_Tb(itori,itori1),thybt2(1),dcos(theta22))
+ gradthybt2=gradtschebyshev
+ & (0,nterm_kcc_Tb(itori,itori1)-1,thybt2(1),dcos(theta22))
+ & *(nterm_kcc_Tb(itori,itori1))*0.5*dsin(theta22)
+
+C now overal sumation
+ etors=etors+(1.0d0+sumth1thyb+sumth2thyb)*sumnonchebyshev
+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*
+ & (glocit1*(1.0d0+sumth1thyb+sumth2thyb)+
+ & sumnonchebyshev*gradthybt1)
+C now derivative over theta2
+ gloc(nphi+i-2,icg)=gloc(nphi+i-2,icg)+wang*
+ & (glocit2*(1.0d0+sumth1thyb+sumth2thyb)+
+ & sumnonchebyshev*gradthybt2)
+ enddo
+
+
+C gloc(i-3,icg)=gloc(i-3,icg)+wtor*gloci
+! 6/20/98 - dihedral angle constraints
+ if (tor_mode.ne.2) then
+ edihcnstr=0.0d0
+c do i=1,ndih_constr
+ do i=idihconstr_start,idihconstr_end
+ itori=idih_constr(i)
+ phii=phi(itori)
+ difi=pinorm(phii-phi0(i))
+ if (difi.gt.drange(i)) then
+ difi=difi-drange(i)
+ edihcnstr=edihcnstr+0.25d0*ftors(i)*difi**4
+ gloc(itori-3,icg)=gloc(itori-3,icg)+ftors(i)*difi**3
+ else if (difi.lt.-drange(i)) then
+ difi=difi+drange(i)
+ edihcnstr=edihcnstr+0.25d0*ftors(i)*difi**4
+ gloc(itori-3,icg)=gloc(itori-3,icg)+ftors(i)*difi**3
+ else
+ difi=0.0
+ endif
+ enddo
+ endif
+ return
+ end
+
+
+
c------------------------------------------------------------------------------
subroutine eback_sc_corr(esccor)
c 7/21/2007 Correlations between the backbone-local and side-chain-local
enddo
return
end
+C--------------------------------------------------------------------------
+ double precision function tschebyshev(m,n,x,y)
+ implicit none
+ include "DIMENSIONS"
+ integer i,m,n
+ double precision x(n),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)=y
+ do i=2,n
+ yy(i)=2*yy(1)*yy(i-1)-yy(i-2)
+ enddo
+ aux=0.0d0
+ do i=m,n
+ aux=aux+x(i)*yy(i)
+ enddo
+ tschebyshev=aux
+ return
+ end
+C--------------------------------------------------------------------------
+ double precision function gradtschebyshev(m,n,x,y)
+ implicit none
+ include "DIMENSIONS"
+ integer i,m,n
+ double precision x(n),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)
+ enddo
+ aux=0.0d0
+ do i=m,n
+ aux=aux+x(i)*yy(i)
+ enddo
+ gradtschebyshev=aux
+ return
+ end
projects / unres.git / blobdiff