X-Git-Url: http://mmka.chem.univ.gda.pl/gitweb/?a=blobdiff_plain;f=source%2Funres%2Fsrc_MD-M%2Fenergy_p_new_barrier.F;h=8fdffe3b05a3697357cc22f8b838c1d4037a4127;hb=aedd2f72960d9cec69fcc7f0b6fd6555c7e22312;hp=a889c75dfd9a6bdde55732701d63794cfe729a6c;hpb=9111b59335b350bc16fcac69159ec3d3cbb62ba4;p=unres.git diff --git a/source/unres/src_MD-M/energy_p_new_barrier.F b/source/unres/src_MD-M/energy_p_new_barrier.F index a889c75..8fdffe3 100644 --- a/source/unres/src_MD-M/energy_p_new_barrier.F +++ b/source/unres/src_MD-M/energy_p_new_barrier.F @@ -201,7 +201,14 @@ C 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 @@ -218,7 +225,14 @@ C Calculate the virtual-bond torsional energy. 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 @@ -227,7 +241,7 @@ c print *,"Processor",myrank," computed Utor" 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 @@ -7432,6 +7446,140 @@ C v1cij=v1c(1,j,itori,itori1,itori2,iblock,ntblock) 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 @@ -11171,4 +11319,46 @@ C write(2,*) "TOTAL VOLUME",i,VolumeTotal,fac_shield(i) 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