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=c91990790254af7ebbd469338db9a82d77a64bd5;hb=88af2feafd40050318384e0a14aac77a4f9fffa4;hp=8fdffe3b05a3697357cc22f8b838c1d4037a4127;hpb=aedd2f72960d9cec69fcc7f0b6fd6555c7e22312;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 8fdffe3..c919907 100644 --- a/source/unres/src_MD-M/energy_p_new_barrier.F +++ b/source/unres/src_MD-M/energy_p_new_barrier.F @@ -224,6 +224,7 @@ C 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) @@ -231,7 +232,7 @@ cd print *,'nterm=',nterm 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 @@ -7468,7 +7469,10 @@ C The rigorous attempt to derive energy function 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. @@ -7508,50 +7512,59 @@ C v1ij is c_n and d_n in euation above 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 @@ -7578,8 +7591,81 @@ c do i=1,ndih_constr 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 @@ -11345,18 +11431,19 @@ C-------------------------------------------------------------------------- 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