first introduction of valence constrains - not working yet

diff --git a/source/unres/src_MD-M/COMMON.TORCNSTR b/source/unres/src_MD-M/COMMON.TORCNSTR
index ec35b4f..4d77fc6 100644 (file)
--- a/source/unres/src_MD-M/COMMON.TORCNSTR
+++ b/source/unres/src_MD-M/COMMON.TORCNSTR
@@ -1,7 +1,10 @@
-      integer ndih_constr,idih_constr(maxdih_constr)
+      integer ndih_constr,idih_constr(maxdih_constr),ntheta_constr,
+     & itheta_constr(maxdih_constr)
       integer ndih_nconstr,idih_nconstr(maxdih_constr)
-      integer idihconstr_start,idihconstr_end
+      integer idihconstr_start,idihconstr_end,ithetaconstr_start,
+     & ithetaconstr_end
       double precision phi0(maxdih_constr),drange(maxdih_constr),
-     & ftors(maxdih_constr)
+     & ftors(maxdih_constr),theta_constr0(maxdih_constr),
+     & theta_drange(maxdih_constr),for_thet_constr(maxdih_constr)
       common /torcnstr/ phi0,drange,ftors,ndih_constr,idih_constr,
      &  ndih_nconstr,idih_nconstr,idihconstr_start,idihconstr_end

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 766d2d5..ec0aef9 100644 (file)
--- a/source/unres/src_MD-M/energy_p_new_barrier.F
+++ b/source/unres/src_MD-M/energy_p_new_barrier.F
@@ -190,7 +190,7 @@ C
 C Calculate the virtual-bond-angle energy.
 C
       if (wang.gt.0d0) then
-        call ebend(ebe)
+        call ebend(ebe,ethetacnstr)
       else
         ebe=0
       endif
@@ -301,6 +301,7 @@ C
       energia(17)=estr
       energia(20)=Uconst+Uconst_back
       energia(21)=esccor
+C      energia(22)=
 c    Here are the energies showed per procesor if the are more processors 
 c    per molecule then we sum it up in sum_energy subroutine 
 c      print *," Processor",myrank," calls SUM_ENERGY"
@@ -977,7 +978,8 @@ C------------------------------------------------------------------------
      &  ecorr,wcorr,
      &  ecorr5,wcorr5,ecorr6,wcorr6,eel_loc,wel_loc,eello_turn3,wturn3,
      &  eello_turn4,wturn4,eello_turn6,wturn6,esccor,wsccor,
-     &  edihcnstr,ebr*nss,
+     &  edihcnstr,
+     &  ethetacnstr,ebr*nss,
      &  Uconst,etot
    10 format (/'Virtual-chain energies:'//
      & 'EVDW=  ',1pE16.6,' WEIGHT=',1pD16.6,' (SC-SC)'/
@@ -1000,6 +1002,7 @@ C------------------------------------------------------------------------
      & 'ETURN6=',1pE16.6,' WEIGHT=',1pD16.6,' (turns, 6th order)'/
      & 'ESCCOR=',1pE16.6,' WEIGHT=',1pD16.6,' (backbone-rotamer corr)'/
      & 'EDIHC= ',1pE16.6,' (dihedral angle constraints)'/
+     & 'ETHETC= ',1pE16.6,' (valence angle constraints)'/
      & 'ESS=   ',1pE16.6,' (disulfide-bridge intrinsic energy)'/
      & 'UCONST= ',1pE16.6,' (Constraint energy)'/ 
      & 'ETOT=  ',1pE16.6,' (total)')
@@ -1010,6 +1013,7 @@ C------------------------------------------------------------------------
      &  ecorr,wcorr,
      &  ecorr5,wcorr5,ecorr6,wcorr6,eel_loc,wel_loc,eello_turn3,wturn3,
      &  eello_turn4,wturn4,eello_turn6,wturn6,esccor,wsccro,edihcnstr,
+     &  ethetacnstr,
      &  ebr*nss,Uconst,etot
    10 format (/'Virtual-chain energies:'//
      & 'EVDW=  ',1pE16.6,' WEIGHT=',1pD16.6,' (SC-SC)'/
@@ -1031,6 +1035,7 @@ C------------------------------------------------------------------------
      & 'ETURN6=',1pE16.6,' WEIGHT=',1pD16.6,' (turns, 6th order)'/
      & 'ESCCOR=',1pE16.6,' WEIGHT=',1pD16.6,' (backbone-rotamer corr)'/
      & 'EDIHC= ',1pE16.6,' (dihedral angle constraints)'/
+     & 'ETHETC= ',1pE16.6,' (valence angle constraints)'/
      & 'ESS=   ',1pE16.6,' (disulfide-bridge intrinsic energy)'/
      & 'UCONST=',1pE16.6,' (Constraint energy)'/ 
      & 'ETOT=  ',1pE16.6,' (total)')
@@ -4520,6 +4525,31 @@ C Derivatives of the "mean" values in gamma1 and gamma2.
      &        E_theta,E_tc)
         endif
         etheta=etheta+ethetai
+      do i=ithetaconstr_start,ithetaconstr_end
+        itheta=itheta_constr(i)
+        thetiii=theta(itori)
+        difi=pinorm(thetiii-theta_constr0(i))
+        if (difi.gt.theta_drange(i)) then
+          difi=difi-theta_drange(i)
+          ethetcnstr=ethetcnstr+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)
+          ethetcnstr=ethetcnstr+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
         if (energy_dec) write (iout,'(a6,i5,0pf7.3)')
      &      'ebend',i,ethetai
         if (i.gt.3) gloc(i-3,icg)=gloc(i-3,icg)+wang*E_tc*dthetg1
@@ -4832,6 +4862,33 @@ C        print *,"tu",cosph1(k),sinph1(k),cosph2(k),sinph2(k)
           enddo
         enddo
 10      continue
+C now we have the theta_constrains
+      do i=ithetaconstr_start,ithetaconstr_end
+        itheta=itheta_constr(i)
+        thetiii=theta(itori)
+        difi=pinorm(thetiii-theta_constr0(i))
+        if (difi.gt.theta_drange(i)) then
+          difi=difi-theta_drange(i)
+          ethetcnstr=ethetcnstr+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)
+          ethetcnstr=ethetcnstr+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
+
 c        lprn1=.true.
 C        print *,ethetai
         if (lprn1) 

diff --git a/source/unres/src_MD-M/initialize_p.F b/source/unres/src_MD-M/initialize_p.F
index f520972..6c5f4eb 100644 (file)
--- a/source/unres/src_MD-M/initialize_p.F
+++ b/source/unres/src_MD-M/initialize_p.F
@@ -629,6 +629,13 @@ c     &  " ivec_start",ivec_start," ivec_end",ivec_end
       else
         call int_bounds(ndih_constr,idihconstr_start,idihconstr_end)
       endif
+      if (ntheta_constr.eq.0) then
+        idihconstr_start=1
+        idihconstr_end=0
+      else
+        call int_bounds
+     &  (ntheta_constr,ithetaconstr_start,ithetaconstr_end)
+      endif
 c      nsumgrad=(nres-nnt)*(nres-nnt+1)/2
 c      nlen=nres-nnt+1
       nsumgrad=(nres-nnt)*(nres-nnt+1)/2
@@ -664,6 +671,9 @@ c      nlen=nres-nnt+1
      & ' iset_start',iset_start,' iset_end',iset_end,
      & ' idihconstr_start',idihconstr_start,' idihconstr_end',
      &   idihconstr_end
+     & ' ithetaconstr_start',ithetaconstr_start,' ithetaconstr_end',
+     &   ithetaconstr_end
+
        write (*,*) 'Processor:',fg_rank,myrank,' igrad_start',
      &   igrad_start,' igrad_end',igrad_end,' ngrad_start',ngrad_start,
      &   ' ngrad_end',ngrad_end
@@ -1128,6 +1138,8 @@ c        write (iout,*) "MPI_ROTAT2",MPI_ROTAT2
       iphi1_end=nres
       idihconstr_start=1
       idihconstr_end=ndih_constr
+      ithetaconstr_start=1
+      ithetaconstr_end=ntheta_constr
       iphid_start=iphi_start
       iphid_end=iphi_end-1
       itau_start=4

diff --git a/source/unres/src_MD-M/readrtns_CSA.F b/source/unres/src_MD-M/readrtns_CSA.F
index 1cda5ea..286cc57 100644 (file)
--- a/source/unres/src_MD-M/readrtns_CSA.F
+++ b/source/unres/src_MD-M/readrtns_CSA.F
@@ -96,6 +96,10 @@ c      print *,"Processor",me," fg_rank",fg_rank," out1file",out1file
 C Set up the time limit (caution! The time must be input in minutes!)
       read_cart=index(controlcard,'READ_CART').gt.0
       call readi(controlcard,'CONSTR_DIST',constr_dist,0)
+C this variable with_theta_constr is the variable which allow to read and execute the
+C constrains on theta angles WITH_THETA_CONSTR is the keyword
+      with_theta_constr = index(controlcard,"WITH_THETA_CONSTR").gt.0
+      write (iout,*) "with_theta_constr ",with_theta_constr
       call readi(controlcard,'SYM',symetr,1)
       call reada(controlcard,'TIMLIM',timlim,960.0D0) ! default 16 hours
       unres_pdb = index(controlcard,'UNRES_PDB') .gt. 0
@@ -851,6 +855,54 @@ C     &   write (iout,*) 'FTORS',ftors
           phibound(2,ii) = phi0(i)+drange(i)
         enddo 
       endif
+C first setting the theta boundaries to 0 to pi
+C this mean that there is no energy penalty for any angle occuring
+      do i=1,nres
+        thetabound(1,i)=0
+        thetabound(2,i)=pi
+      enddo
+C begin reading theta constrains this is quartic constrains allowing to 
+C have smooth second derivative 
+      if (with_theta_constr) then
+C with_theta_constr is keyword allowing for occurance of theta constrains
+      read (inp,*) ntheta_constr
+C ntheta_constr is the number of theta constrains
+      if (ntheta_constr.gt.0) then
+C        read (inp,*) ftors
+        read (inp,*) (itheta_constr(i),theta_constr0(i),
+     &  theta_drange(i),for_thet_constr(i),
+     &  i=1,ntheta_constr)
+C the above code reads from 1 to ntheta_constr 
+C itheta_constr(i) residue i for which is theta_constr
+C theta_constr0 the global minimum value
+C theta_drange is range for which there is no energy penalty
+C for_thet_constr is the force constant for quartic energy penalty
+C E=k*x**4 
+        if(me.eq.king.or..not.out1file)then
+         write (iout,*)
+     &   'There are',ntheta_constr,' constraints on phi angles.'
+         do i=1,ntheta_constr
+          write (iout,'(i5,3f8.3)') itheta_constr(i),theta_constr0(i),
+     &    theta_drange(i),
+     &    for_thet_constr(i)
+         enddo
+        endif
+        do i=1,nthet_constr
+          theta_constr0(i)=deg2rad*theta_constr0(i)
+          theta_drange(i)=deg2rad*theta_drange(i)
+        enddo
+C        if(me.eq.king.or..not.out1file)
+C     &   write (iout,*) 'FTORS',ftors
+        do i=1,ntheta_constr
+          ii = itheta_constr(i)
+          thetabound(1,ii) = phi0(i)-drange(i)
+          thetabound(2,ii) = phi0(i)+drange(i)
+        enddo
+      endif ! ntheta_constr.gt.0
+      endif! with_theta_constr
+C
+C      with_dihed_constr = index(controlcard,"WITH_DIHED_CONSTR").gt.0
+C      write (iout,*) "with_dihed_constr ",with_dihed_constr
       nnt=1
 #ifdef MPI
       if (me.eq.king) then

diff --git a/source/wham/src-M/readrtns.F b/source/wham/src-M/readrtns.F
index 33924b5..5dd092e 100644 (file)
--- a/source/wham/src-M/readrtns.F
+++ b/source/wham/src-M/readrtns.F
@@ -91,6 +91,8 @@
       zscfile=index(controlcard,"ZSCFILE").gt.0
       with_dihed_constr = index(controlcard,"WITH_DIHED_CONSTR").gt.0
       write (iout,*) "with_dihed_constr ",with_dihed_constr
+      with_theta_constr = index(controlcard,"WITH_THETA_CONSTR").gt.0
+      write (iout,*) "with_theta_constr ",with_theta_constr
       call readi(controlcard,'CONSTR_DIST',constr_dist,0)
       dyn_ss=index(controlcard,"DYN_SS").gt.0
       return