From 88af2feafd40050318384e0a14aac77a4f9fffa4 Mon Sep 17 00:00:00 2001 From: Adam Sieradzan Date: Thu, 26 Nov 2015 11:43:25 +0100 Subject: [PATCH] implementation of czybyszeb grad ok --- source/unres/src_MD-M/COMMON.CONTROL | 2 +- source/unres/src_MD-M/COMMON.TORSION | 11 ++- source/unres/src_MD-M/DIMENSIONS | 2 +- source/unres/src_MD-M/eigen.f | 2 +- source/unres/src_MD-M/energy_p_new_barrier.F | 131 +++++++++++++++++++++----- source/unres/src_MD-M/parmread.F | 30 ++++-- 6 files changed, 140 insertions(+), 38 deletions(-) diff --git a/source/unres/src_MD-M/COMMON.CONTROL b/source/unres/src_MD-M/COMMON.CONTROL index a6cb608..d03828e 100644 --- a/source/unres/src_MD-M/COMMON.CONTROL +++ b/source/unres/src_MD-M/COMMON.CONTROL @@ -9,7 +9,7 @@ & icheckgrad,minim,i2ndstr,refstr,pdbref,outpdb,outmol2,iprint, & overlapsc,energy_dec,sideadd,lsecondary,read_cart,unres_pdb & ,vdisulf,searchsc,lmuca,dccart,mucadyn,extconf,out1file, - & selfguide,AFMlog,shield_mode, + & selfguide,AFMlog,shield_mode,tor_mode, & constr_dist,gnorm_check,gradout,split_ene,with_theta_constr, & symetr C... minim = .true. means DO minimization. diff --git a/source/unres/src_MD-M/COMMON.TORSION b/source/unres/src_MD-M/COMMON.TORSION index d9a2fea..56c86ac 100644 --- a/source/unres/src_MD-M/COMMON.TORSION +++ b/source/unres/src_MD-M/COMMON.TORSION @@ -1,7 +1,8 @@ C Torsional constants of the rotation about virtual-bond dihedral angles - double precision v1,v2,vlor1,vlor2,vlor3,v0 + double precision v1,v2,vlor1,vlor2,vlor3,v0,v1_kcc,v2_kcc, + & v1_chyb,v2_chyb,v1bend_chyb integer itortyp,ntortyp,nterm,nlor,nterm_old,nterm_kcc_Tb, - & nterm_kcc,itortyp_kcc + & nterm_kcc,itortyp_kcc,nbend_kcc_Tb common/torsion/v0(-maxtor:maxtor,-maxtor:maxtor,2), & v1(maxterm,-maxtor:maxtor,-maxtor:maxtor,2), & v2(maxterm,-maxtor:maxtor,-maxtor:maxtor,2), @@ -9,14 +10,16 @@ C Torsional constants of the rotation about virtual-bond dihedral angles & vlor2(maxlor,maxtor,maxtor),vlor3(maxlor,maxtor,maxtor), & v1_kcc(maxtermkcc,-maxtor:maxtor,-maxtor:maxtor), & v2_kcc(maxtermkcc,-maxtor:maxtor,-maxtor:maxtor), - & v1_chyb(maxtermkcc,-maxtor:maxtor,-maxtor:maxtor), - & v2_chyb(maxtermkcc,-maxtor:maxtor,-maxtor:maxtor), + & v1_chyb(maxtermkcc,maxtermkcc,-maxtor:maxtor,-maxtor:maxtor), + & v2_chyb(maxtermkcc,maxtermkcc,-maxtor:maxtor,-maxtor:maxtor), + & v1bend_chyb(maxtermkcc,-maxtor:maxtor), & itortyp(-ntyp1:ntyp1),ntortyp, & itortyp_kcc(-ntyp1:ntyp1), & nterm(-maxtor:maxtor,-maxtor:maxtor,2), & nlor(-maxtor:maxtor,-maxtor:maxtor,2), & nterm_kcc_Tb(-maxtor:maxtor,-maxtor:maxtor), & nterm_kcc(-maxtor:maxtor,-maxtor:maxtor), + & nbend_kcc_Tb(-maxtor:maxtor), & nterm_old C 6/23/01 - constants for double torsionals double precision v1c,v1s,v2c,v2s diff --git a/source/unres/src_MD-M/DIMENSIONS b/source/unres/src_MD-M/DIMENSIONS index 02deda8..4c9a85c 100644 --- a/source/unres/src_MD-M/DIMENSIONS +++ b/source/unres/src_MD-M/DIMENSIONS @@ -50,7 +50,7 @@ C Max. number of types of dihedral angles & multiplicity of torsional barriers C and the number of terms in double torsionals integer maxtor,maxterm,maxlor,maxtermd_1,maxtermd_2,maxtermkcc parameter (maxtor=4,maxterm=10,maxlor=3,maxtermd_1=8,maxtermd_2=8) - parameter (maxtermkcc=20) + parameter (maxtermkcc=40) C Max. number of residue types and parameters in expressions for C virtual-bond angle bending potentials integer maxthetyp,maxthetyp1,maxtheterm,maxtheterm2,maxtheterm3, diff --git a/source/unres/src_MD-M/eigen.f b/source/unres/src_MD-M/eigen.f index e4088ee..287792f 100644 --- a/source/unres/src_MD-M/eigen.f +++ b/source/unres/src_MD-M/eigen.f @@ -127,7 +127,7 @@ C INTEGER GROUP,I,IERR,ITS,J,JJ,M,N,NM,P,Q,R,S,SUBMAT,TAG INTEGER IND(M) C - DOUBLE PRECISION D(N),E(N),E2(N),W(M),Z(NM,M) + DOUBLE PRECISION D(N),E(3*N),E2(N),W(M),Z(NM,M) DOUBLE PRECISION RV1(N),RV2(N),RV3(N),RV4(N),RV6(N) DOUBLE PRECISION ANORM,EPS2,EPS3,EPS4,NORM,ORDER,RHO,U,UK,V DOUBLE PRECISION X0,X1,XU 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 diff --git a/source/unres/src_MD-M/parmread.F b/source/unres/src_MD-M/parmread.F index e156c97..438f544 100644 --- a/source/unres/src_MD-M/parmread.F +++ b/source/unres/src_MD-M/parmread.F @@ -784,22 +784,34 @@ C read Czybyshev torsional parameters do i=1,nkcctyp do j=1,nkcctyp C first we read the cos and sin gamma parameters - read (itorkcc,*,end=121,err=121) nterm_kcc(j,i) + read (itorkcc,*,end=121,err=121) + & nterm_kcc(j,i),nterm_kcc_Tb(j,i) +C read (itorkcc,*,end=121,err=121) nterm_kcc_Tb(j,i) do k=1,nterm_kcc(j,i) - read (itorkcc,*,end=121,err=121) - & v1_kcc(k,j,i),v2_kcc(k,j,i) + do l=1,nterm_kcc_Tb(j,i) + read (itorkcc,*,end=121,err=121) + & v1_chyb(l,k,j,i) + enddo + do l=1,nterm_kcc_Tb(j,i) + read (itorkcc,*,end=121,err=121) + & v2_chyb(l,k,j,i) enddo -C now Czybyshev parameters - read (itorkcc,*,end=121,err=121) nterm_kcc_Tb(j,i) - do k=1,nterm_kcc_Tb(j,i) read (itorkcc,*,end=121,err=121) - & v1_chyb(k,j,i),v2_chyb(k,j,i) + & v1_kcc(k,j,i) + read (itorkcc,*,end=121,err=121) + & v2_kcc(k,j,i) enddo enddo enddo C here will be the apropriate recalibrating for D-aminoacid - - +C read (ithetkcc,*,end=121,err=121) nkcctyp + do i=1,nkcctyp + read (ithetkcc,*,end=121,err=121) nbend_kcc_Tb(i) + do j=1,nbend_kcc_Tb(i) + read (ithetkcc,*,end=121,err=121) + & v1bend_chyb(j,i) + enddo + enddo C Read of Side-chain backbone correlation parameters C Modified 11 May 2012 by Adasko CCC -- 1.7.9.5