implementation of czybyszeb grad ok

diff --git a/source/unres/src_MD-M/COMMON.CONTROL b/source/unres/src_MD-M/COMMON.CONTROL
index a6cb608..d03828e 100644 (file)
--- 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 (file)
--- 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 (file)
--- 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 (file)
--- 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 (file)
--- 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 (file)
--- 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