first introduction of valence constrains - not working yet

[unres.git] / source / unres / src_MD-M / readrtns_CSA.F

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