1) debug of Lorentzian like constrains in cluster 2) introduction of constrains on...

diff --git a/source/cluster/wham/src-M/COMMON.CONTROL b/source/cluster/wham/src-M/COMMON.CONTROL
index e65a5b2..f339ae9 100644 (file)
--- a/source/cluster/wham/src-M/COMMON.CONTROL
+++ b/source/cluster/wham/src-M/COMMON.CONTROL
@@ -3,9 +3,9 @@
      & constr_dist
       logical refstr,pdbref,punch_dist,print_dist,caonly,lside,
      & lprint_cart,lprint_int,from_cart,efree,from_bx,from_cx,
-     & with_dihed_constr
+     & with_dihed_constr,with_theta_constr
       common /cntrl/ betaT,iscode,indpdb,refstr,pdbref,outpdb,outmol2,
      & punch_dist,print_dist,caonly,lside,lprint_cart,lprint_int,
-     & from_cart,from_bx,from_cx, with_dihed_constr,
+     & from_cart,from_bx,from_cx, with_dihed_constr,with_theta_constr,
      & efree,iopt,nstart,nend,symetr,
      & constr_dist

diff --git a/source/cluster/wham/src-M/DIMENSIONS b/source/cluster/wham/src-M/DIMENSIONS
index 68ac673..29e8f24 100644 (file)
--- a/source/cluster/wham/src-M/DIMENSIONS
+++ b/source/cluster/wham/src-M/DIMENSIONS
@@ -60,7 +60,7 @@ C Max number of symetric chains
       parameter (maxperm=120)
 C Max. number of energy components
       integer max_ene
-      parameter (max_ene=21)
+      parameter (max_ene=24)
 C Max. number of temperatures
       integer maxt
       parameter (maxT=5)

diff --git a/source/cluster/wham/src-M/energy_p_new.F b/source/cluster/wham/src-M/energy_p_new.F
index 39d35dd..f78e2e9 100644 (file)
--- a/source/cluster/wham/src-M/energy_p_new.F
+++ b/source/cluster/wham/src-M/energy_p_new.F
@@ -67,7 +67,7 @@ cd    print *,'EHPB exitted succesfully.'
 C
 C Calculate the virtual-bond-angle energy.
 C
-      call ebend(ebe)
+      call ebend(ebe,ethetacnstr)
 cd    print *,'Bend energy finished.'
 C
 C Calculate the SC local energy.
@@ -111,7 +111,7 @@ c      write (iout,*) "ft(6)",fact(6)," evdw",evdw," evdw_t",evdw_t
      & +wcorr6*fact(5)*ecorr6+wturn4*fact(3)*eello_turn4
      & +wturn3*fact(2)*eello_turn3+wturn6*fact(5)*eturn6
      & +wel_loc*fact(2)*eel_loc+edihcnstr+wtor_d*fact(2)*etors_d
-     & +wbond*estr+wsccor*fact(1)*esccor
+     & +wbond*estr+wsccor*fact(1)*esccor+ethetacnstr
 #else
       etot=wsc*(evdw+fact(6)*evdw_t)+wscp*evdw2
      & +welec*fact(1)*(ees+evdw1)
@@ -120,7 +120,7 @@ c      write (iout,*) "ft(6)",fact(6)," evdw",evdw," evdw_t",evdw_t
      & +wcorr6*fact(5)*ecorr6+wturn4*fact(3)*eello_turn4
      & +wturn3*fact(2)*eello_turn3+wturn6*fact(5)*eturn6
      & +wel_loc*fact(2)*eel_loc+edihcnstr+wtor_d*fact(2)*etors_d
-     & +wbond*estr+wsccor*fact(1)*esccor
+     & +wbond*estr+wsccor*fact(1)*esccor+ethetacnstr
 #endif
       energia(0)=etot
       energia(1)=evdw
@@ -154,6 +154,7 @@ c      write (iout,*) "ft(6)",fact(6)," evdw",evdw," evdw_t",evdw_t
       energia(19)=esccor
       energia(20)=edihcnstr
       energia(21)=evdw_t
+      energia(24)=ethetacnstr
 c detecting NaNQ
 #ifdef ISNAN
 #ifdef AIX
@@ -270,6 +271,7 @@ C------------------------------------------------------------------------
       esccor=energia(19)
       edihcnstr=energia(20)
       estr=energia(18)
+      ethetacnstr=energia(24)
 #ifdef SPLITELE
       write (iout,10) evdw,wsc,evdw2,wscp,ees,welec*fact(1),evdw1,
      &  wvdwpp,
@@ -278,7 +280,7 @@ C------------------------------------------------------------------------
      &  ecorr,wcorr*fact(3),ecorr5,wcorr5*fact(4),ecorr6,wcorr6*fact(5),
      &  eel_loc,wel_loc*fact(2),eello_turn3,wturn3*fact(2),
      &  eello_turn4,wturn4*fact(3),eello_turn6,wturn6*fact(5),
-     &  esccor,wsccor*fact(1),edihcnstr,ebr*nss,etot
+     &  esccor,wsccor*fact(1),edihcnstr,ethetacnstr,ebr*nss,etot
    10 format (/'Virtual-chain energies:'//
      & 'EVDW=  ',1pE16.6,' WEIGHT=',1pD16.6,' (SC-SC)'/
      & 'EVDW2= ',1pE16.6,' WEIGHT=',1pD16.6,' (SC-p)'/
@@ -300,6 +302,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)'/ 
      & 'ETOT=  ',1pE16.6,' (total)')
 #else
@@ -309,7 +312,7 @@ C------------------------------------------------------------------------
      &  ecorr6,wcorr6*fact(5),eel_loc,wel_loc*fact(2),
      &  eello_turn3,wturn3*fact(2),eello_turn4,wturn4*fact(3),
      &  eello_turn6,wturn6*fact(5),esccor*fact(1),wsccor,
-     &  edihcnstr,ebr*nss,etot
+     &  edihcnstr,ethetacnstr,ebr*nss,etot
    10 format (/'Virtual-chain energies:'//
      & 'EVDW=  ',1pE16.6,' WEIGHT=',1pD16.6,' (SC-SC)'/
      & 'EVDW2= ',1pE16.6,' WEIGHT=',1pD16.6,' (SC-p)'/
@@ -330,6 +333,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)'/ 
      & 'ETOT=  ',1pE16.6,' (total)')
 #endif
@@ -2909,6 +2913,7 @@ C
       include 'COMMON.DERIV'
       include 'COMMON.VAR'
       include 'COMMON.INTERACT'
+      include 'COMMON.CONTROL'
       dimension ggg(3)
       ehpb=0.0D0
 cd    print *,'edis: nhpb=',nhpb,' fbr=',fbr
@@ -2952,6 +2957,7 @@ C     &          *rlornmr1(dd,dhpb(i),dhpb1(i),forcon(i))
      &          *rlornmr1prim(dd,dhpb(i),dhpb1(i),forcon(i))/dd
 C          write (iout,'(a6,2i5,3f8.3)') "edisl",ii,jj,
 C     &    ehpb,fordepth(i),dd
+C             print *,"TUTU"
 C            write(iout,*) ehpb,"atu?"
 C            ehpb,"tu?"
 C            fac=fordepth(i)**4.0d0
@@ -3228,7 +3234,7 @@ c     &      AKSC(j,iti),abond0(j,iti),u(j),j=1,nbi)
       end
 #ifdef CRYST_THETA
 C--------------------------------------------------------------------------
-      subroutine ebend(etheta)
+      subroutine ebend(etheta,ethetacnstr)
 C
 C Evaluate the virtual-bond-angle energy given the virtual-bond dihedral
 C angles gamma and its derivatives in consecutive thetas and gammas.
@@ -3245,6 +3251,7 @@ C
       include 'COMMON.IOUNITS'
       include 'COMMON.NAMES'
       include 'COMMON.FFIELD'
+      include 'COMMON.TORCNSTR'
       common /calcthet/ term1,term2,termm,diffak,ratak,
      & ak,aktc,termpre,termexp,sigc,sig0i,time11,time12,sigcsq,
      & delthe0,sig0inv,sigtc,sigsqtc,delthec,it
@@ -3363,6 +3370,34 @@ c     &    rad2deg*phii,rad2deg*phii1,ethetai
 c 1215   continue
       enddo
 C Ufff.... We've done all this!!! 
+C now constrains
+      ethetacnstr=0.0d0
+C      print *,ithetaconstr_start,ithetaconstr_end,"TU"
+      do i=1,ntheta_constr
+        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
+C       if (energy_dec) then
+C        write (iout,'(a6,2i5,4f8.3,2e14.5)') "ethetc",
+C     &    i,itheta,rad2deg*thetiii,
+C     &    rad2deg*theta_constr0(i),  rad2deg*theta_drange(i),
+C     &    rad2deg*difi,0.25d0*for_thet_constr(i)*difi**4,
+C     &    gloc(itheta+nphi-2,icg)
+C        endif
+      enddo
       return
       end
 C---------------------------------------------------------------------------
@@ -3475,7 +3510,7 @@ C "Thank you" to MAPLE (probably spared one day of hand-differentiation).
       end
 #else
 C--------------------------------------------------------------------------
-      subroutine ebend(etheta)
+      subroutine ebend(etheta,ethetacnstr)
 C
 C Evaluate the virtual-bond-angle energy given the virtual-bond dihedral
 C angles gamma and its derivatives in consecutive thetas and gammas.
@@ -3495,6 +3530,7 @@ C
       include 'COMMON.NAMES'
       include 'COMMON.FFIELD'
       include 'COMMON.CONTROL'
+      include 'COMMON.TORCNSTR'
       double precision coskt(mmaxtheterm),sinkt(mmaxtheterm),
      & cosph1(maxsingle),sinph1(maxsingle),cosph2(maxsingle),
      & sinph2(maxsingle),cosph1ph2(maxdouble,maxdouble),
@@ -3675,6 +3711,34 @@ c        call flush(iout)
 c        gloc(nphi+i-2,icg)=wang*dethetai
         gloc(nphi+i-2,icg)=gloc(nphi+i-2,icg)+wang*dethetai
       enddo
+C now constrains
+      ethetacnstr=0.0d0
+C      print *,ithetaconstr_start,ithetaconstr_end,"TU"
+      do i=1,ntheta_constr
+        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
+C       if (energy_dec) then
+C        write (iout,'(a6,2i5,4f8.3,2e14.5)') "ethetc",
+C     &    i,itheta,rad2deg*thetiii,
+C     &    rad2deg*theta_constr0(i),  rad2deg*theta_drange(i),
+C     &    rad2deg*difi,0.25d0*for_thet_constr(i)*difi**4,
+C     &    gloc(itheta+nphi-2,icg)
+C        endif
+      enddo
       return
       end
 #endif

diff --git a/source/cluster/wham/src-M/gnmr1.f b/source/cluster/wham/src-M/gnmr1.f
index 905e746..2357e6d 100644 (file)
--- a/source/cluster/wham/src-M/gnmr1.f
+++ b/source/cluster/wham/src-M/gnmr1.f
@@ -41,3 +41,34 @@ c-------------------------------------------------------------------------------
       return
       end
 c---------------------------------------------------------------------------------
+c---------------------------------------------------------------------------------
+      double precision function rlornmr1(y,ymin,ymax,sigma)
+      implicit none
+      double precision y,ymin,ymax,sigma
+      double precision wykl /4.0d0/
+      if (y.lt.ymin) then
+        rlornmr1=(ymin-y)**wykl/((ymin-y)**wykl+sigma**wykl)
+      else if (y.gt.ymax) then
+        rlornmr1=(y-ymax)**wykl/((y-ymax)**wykl+sigma**wykl)
+      else
+        rlornmr1=0.0d0
+      endif
+      return
+      end
+c------------------------------------------------------------------------------
+      double precision function rlornmr1prim(y,ymin,ymax,sigma)
+      implicit none
+      double precision y,ymin,ymax,sigma
+      double precision wykl /4.0d0/
+      if (y.lt.ymin) then
+        rlornmr1prim=-(ymin-y)**(wykl-1)*sigma**wykl*wykl/
+     &   ((ymin-y)**wykl+sigma**wykl)**2
+      else if (y.gt.ymax) then
+        rlornmr1prim=(y-ymax)**(wykl-1)*sigma**wykl*wykl/
+     & ((y-ymax)**wykl+sigma**wykl)**2
+      else
+        rlornmr1prim=0.0d0
+      endif
+      return
+      end
+

diff --git a/source/cluster/wham/src-M/include_unres/COMMON.TORCNSTR b/source/cluster/wham/src-M/include_unres/COMMON.TORCNSTR
index 5c26487..845c44d 100644 (file)
--- a/source/cluster/wham/src-M/include_unres/COMMON.TORCNSTR
+++ b/source/cluster/wham/src-M/include_unres/COMMON.TORCNSTR
@@ -1,6 +1,14 @@
-      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,ithetaconstr_start,
+     & ithetaconstr_end
       double precision phi0(maxdih_constr),drange(maxdih_constr),
-     &   ftors(maxdih_constr)
-      common /torcnstr/ phi0,drange,ftors,ndih_constr,idih_constr,
-     &                  ndih_nconstr,idih_nconstr
+     & ftors(maxdih_constr),theta_constr0(maxdih_constr),
+     & theta_drange(maxdih_constr),for_thet_constr(maxdih_constr)
+      common /torcnstr/ phi0,drange,ftors,theta_constr0,theta_drange,
+     & for_thet_constr,
+     &  ndih_constr,idih_constr,
+     &  ndih_nconstr,idih_nconstr,idihconstr_start,idihconstr_end,
+     & ntheta_constr,itheta_constr,ithetaconstr_start,
+     & ithetaconstr_end

diff --git a/source/cluster/wham/src-M/initialize_p.F b/source/cluster/wham/src-M/initialize_p.F
index 4ebced6..f385f1a 100644 (file)
--- a/source/cluster/wham/src-M/initialize_p.F
+++ b/source/cluster/wham/src-M/initialize_p.F
@@ -256,14 +256,16 @@ c-------------------------------------------------------------------------
      &   "EVDW SC-SC","EVDW2 SC-p","EES p-p","ECORR4 ","ECORR5 ",
      &   "ECORR6 ","EELLO ","ETURN3 ","ETURN4 ","ETURN6 ",
      &   "EBE bend","ESC SCloc","ETORS ","ETORSD ","EHPB","EVDWPP",
-     &   "ESTR","ESCCOR","EVDW2_14","","EVDW_T"/
+     &   "EVDW2_14","ESTR","ESCCOR","EDIHC","EVDW_T","ELIPTRAN",
+     &   "EAFM","ETHETC"/
       data wname /
      &   "WSC","WSCP","WELEC","WCORR","WCORR5","WCORR6","WEL_LOC",
      &   "WTURN3","WTURN4","WTURN6","WANG","WSCLOC","WTOR","WTORD",
-     &   "WHPB","WVDWPP","WBOND","WSCCOR","WSCP14","","WSC"/
-      data nprint_ene /21/
+     &   "WHPB","WVDWPP","WSCP14","WBOND","WSCCOR","WDIHC","WSC",
+     &   "WLIPTRAN","WAFM","WTHETC"/
+      data nprint_ene /22/
       data print_order /1,2,3,18,11,12,13,14,4,5,6,7,8,9,10,19,
-     &  16,15,17,20,21/
+     &  16,15,17,20,21,24,22,23/
       end 
 c---------------------------------------------------------------------------
       subroutine init_int_table

diff --git a/source/cluster/wham/src-M/readrtns.F b/source/cluster/wham/src-M/readrtns.F
index 8624355..aa85ca1 100644 (file)
--- a/source/cluster/wham/src-M/readrtns.F
+++ b/source/cluster/wham/src-M/readrtns.F
@@ -39,6 +39,8 @@ C
       call readi(controlcard,'CONSTR_DIST',constr_dist,0)
       write (iout,*) "with_dihed_constr ",with_dihed_constr,
      & " CONSTR_DIST",constr_dist
+      with_theta_constr = index(controlcard,"WITH_THETA_CONSTR").gt.0
+      write (iout,*) "with_theta_constr ",with_theta_constr
       call flush(iout)
       min_var=(index(controlcard,'MINVAR').gt.0)
       plot_tree=(index(controlcard,'PLOT_TREE').gt.0)
@@ -306,6 +308,44 @@ C ftors is the force constant for torsional quartic constrains
         enddo
       endif ! endif ndif_constr.gt.0
       endif ! with_dihed_constr
+      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 
+C        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
+C        endif
+        do i=1,ntheta_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
+C        do i=1,ntheta_constr
+C          ii = itheta_constr(i)
+C          thetabound(1,ii) = phi0(i)-drange(i)
+C          thetabound(2,ii) = phi0(i)+drange(i)
+C        enddo
+      endif ! ntheta_constr.gt.0
+      endif! with_theta_constr
+
       nnt=1
       nct=nres
       print *,'NNT=',NNT,' NCT=',NCT

diff --git a/source/unres/src_MD-M/COMMON.CONTROL b/source/unres/src_MD-M/COMMON.CONTROL
index 40346c7..4d05560 100644 (file)
--- a/source/unres/src_MD-M/COMMON.CONTROL
+++ b/source/unres/src_MD-M/COMMON.CONTROL
@@ -3,11 +3,12 @@
       logical minim,refstr,pdbref,outpdb,outmol2,overlapsc,energy_dec,
      &                 sideadd,lsecondary,read_cart,unres_pdb,
      &                 vdisulf,searchsc,lmuca,dccart,extconf,out1file,
-     &                 gnorm_check,gradout,split_ene
+     &                 gnorm_check,gradout,split_ene,with_theta_constr
       common /cntrl/ modecalc,iscode,indpdb,indback,indphi,iranconf,
      & icheckgrad,minim,i2ndstr,refstr,pdbref,outpdb,outmol2,iprint,
      & overlapsc,energy_dec,sideadd,lsecondary,read_cart,unres_pdb
      & ,vdisulf,searchsc,lmuca,dccart,mucadyn,extconf,out1file,
-     & constr_dist,gnorm_check,gradout,split_ene,symetr
+     & constr_dist,gnorm_check,gradout,split_ene,with_theta_constr,
+     & symetr
 C... minim = .true. means DO minimization.
 C... energy_dec = .true. means print energy decomposition matrix

diff --git a/source/unres/src_MD-M/COMMON.TORCNSTR b/source/unres/src_MD-M/COMMON.TORCNSTR
index 4d77fc6..845c44d 100644 (file)
--- a/source/unres/src_MD-M/COMMON.TORCNSTR
+++ b/source/unres/src_MD-M/COMMON.TORCNSTR
@@ -6,5 +6,9 @@
       double precision phi0(maxdih_constr),drange(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
+      common /torcnstr/ phi0,drange,ftors,theta_constr0,theta_drange,
+     & for_thet_constr,
+     &  ndih_constr,idih_constr,
+     &  ndih_nconstr,idih_nconstr,idihconstr_start,idihconstr_end,
+     & ntheta_constr,itheta_constr,ithetaconstr_start,
+     & ithetaconstr_end

diff --git a/source/unres/src_MD-M/DIMENSIONS b/source/unres/src_MD-M/DIMENSIONS
index c8a95df..bda06b4 100644 (file)
--- a/source/unres/src_MD-M/DIMENSIONS
+++ b/source/unres/src_MD-M/DIMENSIONS
@@ -95,7 +95,7 @@ C Max. number of conformations in the pool
       parameter (max_pool=10)
 C Number of energy components
       integer n_ene,n_ene2
-      parameter (n_ene=21,n_ene2=2*n_ene)
+      parameter (n_ene=25,n_ene2=2*n_ene)
 C Number of threads in deformation
       integer max_thread,max_thread2
       parameter (max_thread=4,max_thread2=2*max_thread)     

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 ec0aef9..9414f1c 100644 (file)
--- a/source/unres/src_MD-M/energy_p_new_barrier.F
+++ b/source/unres/src_MD-M/energy_p_new_barrier.F
@@ -193,6 +193,7 @@ C
         call ebend(ebe,ethetacnstr)
       else
         ebe=0
+        ethetacnstr=0
       endif
 c      print *,"Processor",myrank," computed UB"
 C
@@ -301,7 +302,9 @@ C
       energia(17)=estr
       energia(20)=Uconst+Uconst_back
       energia(21)=esccor
-C      energia(22)=
+C      energia(22)=eliptrans (the energy for lipid transfere implemented in lipid branch)
+C      energia(23)= ... (energy for AFM, steered molecular dynamics)
+      energia(24)=ethetacnstr
 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"
@@ -393,20 +396,22 @@ cMS$ATTRIBUTES C ::  proc_proc
       estr=energia(17)
       Uconst=energia(20)
       esccor=energia(21)
+      ethetacnstr=energia(24)
+
 #ifdef SPLITELE
       etot=wsc*evdw+wscp*evdw2+welec*ees+wvdwpp*evdw1
      & +wang*ebe+wtor*etors+wscloc*escloc
      & +wstrain*ehpb+wcorr*ecorr+wcorr5*ecorr5
      & +wcorr6*ecorr6+wturn4*eello_turn4+wturn3*eello_turn3
      & +wturn6*eturn6+wel_loc*eel_loc+edihcnstr+wtor_d*etors_d
-     & +wbond*estr+Uconst+wsccor*esccor
+     & +wbond*estr+Uconst+wsccor*esccor+ethetacnstr
 #else
       etot=wsc*evdw+wscp*evdw2+welec*(ees+evdw1)
      & +wang*ebe+wtor*etors+wscloc*escloc
      & +wstrain*ehpb+wcorr*ecorr+wcorr5*ecorr5
      & +wcorr6*ecorr6+wturn4*eello_turn4+wturn3*eello_turn3
      & +wturn6*eturn6+wel_loc*eel_loc+edihcnstr+wtor_d*etors_d
-     & +wbond*estr+Uconst+wsccor*esccor
+     & +wbond*estr+Uconst+wsccor*esccor+ethetacnstr
 #endif
       energia(0)=etot
 c detecting NaNQ
@@ -971,6 +976,7 @@ C------------------------------------------------------------------------
       estr=energia(17)
       Uconst=energia(20)
       esccor=energia(21)
+      ethetacnstr=energia(24)
 #ifdef SPLITELE
       write (iout,10) evdw,wsc,evdw2,wscp,ees,welec,evdw1,wvdwpp,
      &  estr,wbond,ebe,wang,
@@ -4405,7 +4411,7 @@ c
       end 
 #ifdef CRYST_THETA
 C--------------------------------------------------------------------------
-      subroutine ebend(etheta)
+      subroutine ebend(etheta,ethetacnstr)
 C
 C Evaluate the virtual-bond-angle energy given the virtual-bond dihedral
 C angles gamma and its derivatives in consecutive thetas and gammas.
@@ -4422,6 +4428,7 @@ C
       include 'COMMON.NAMES'
       include 'COMMON.FFIELD'
       include 'COMMON.CONTROL'
+      include 'COMMON.TORCNSTR'
       common /calcthet/ term1,term2,termm,diffak,ratak,
      & ak,aktc,termpre,termexp,sigc,sig0i,time11,time12,sigcsq,
      & delthe0,sig0inv,sigtc,sigsqtc,delthec,it
@@ -4525,18 +4532,26 @@ C Derivatives of the "mean" values in gamma1 and gamma2.
      &        E_theta,E_tc)
         endif
         etheta=etheta+ethetai
+        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
+        if (i.lt.nres) gloc(i-2,icg)=gloc(i-2,icg)+wang*E_tc*dthetg2
+        gloc(nphi+i-2,icg)=wang*(E_theta+E_tc*dthett)
+      enddo
+      ethetacnstr=0.0d0
+C      print *,ithetaconstr_start,ithetaconstr_end,"TU"
       do i=ithetaconstr_start,ithetaconstr_end
         itheta=itheta_constr(i)
-        thetiii=theta(itori)
+        thetiii=theta(itheta)
         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
+          ethetacnstr=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
+          ethetacnstr=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
@@ -4550,12 +4565,7 @@ C Derivatives of the "mean" values in gamma1 and gamma2.
      &    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
-        if (i.lt.nres) gloc(i-2,icg)=gloc(i-2,icg)+wang*E_tc*dthetg2
-        gloc(nphi+i-2,icg)=wang*(E_theta+E_tc*dthett)
-      enddo
+
 C Ufff.... We've done all this!!! 
       return
       end
@@ -4669,7 +4679,7 @@ C "Thank you" to MAPLE (probably spared one day of hand-differentiation).
       end
 #else
 C--------------------------------------------------------------------------
-      subroutine ebend(etheta)
+      subroutine ebend(etheta,ethetacnstr)
 C
 C Evaluate the virtual-bond-angle energy given the virtual-bond dihedral
 C angles gamma and its derivatives in consecutive thetas and gammas.
@@ -4688,6 +4698,7 @@ C
       include 'COMMON.NAMES'
       include 'COMMON.FFIELD'
       include 'COMMON.CONTROL'
+      include 'COMMON.TORCNSTR'
       double precision coskt(mmaxtheterm),sinkt(mmaxtheterm),
      & cosph1(maxsingle),sinph1(maxsingle),cosph2(maxsingle),
      & sinph2(maxsingle),cosph1ph2(maxdouble,maxdouble),
@@ -4862,19 +4873,33 @@ C        print *,"tu",cosph1(k),sinph1(k),cosph2(k),sinph2(k)
           enddo
         enddo
 10      continue
-C now we have the theta_constrains
+c        lprn1=.true.
+C        print *,ethetai
+        if (lprn1) 
+     &    write (iout,'(i2,3f8.1,9h ethetai ,f10.5)') 
+     &   i,theta(i)*rad2deg,phii*rad2deg,
+     &   phii1*rad2deg,ethetai
+c        lprn1=.false.
+        etheta=etheta+ethetai
+        if (i.gt.3) gloc(i-3,icg)=gloc(i-3,icg)+wang*dephii
+        if (i.lt.nres) gloc(i-2,icg)=gloc(i-2,icg)+wang*dephii1
+        gloc(nphi+i-2,icg)=gloc(nphi+i-2,icg)+wang*dethetai
+      enddo
+C now constrains
+      ethetacnstr=0.0d0
+C      print *,ithetaconstr_start,ithetaconstr_end,"TU"
       do i=ithetaconstr_start,ithetaconstr_end
         itheta=itheta_constr(i)
-        thetiii=theta(itori)
+        thetiii=theta(itheta)
         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
+          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)
-          ethetcnstr=ethetcnstr+0.25d0*for_thet_constr(i)*difi**4
+          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
@@ -4889,18 +4914,6 @@ C now we have the theta_constrains
         endif
       enddo
 
-c        lprn1=.true.
-C        print *,ethetai
-        if (lprn1) 
-     &    write (iout,'(i2,3f8.1,9h ethetai ,f10.5)') 
-     &   i,theta(i)*rad2deg,phii*rad2deg,
-     &   phii1*rad2deg,ethetai
-c        lprn1=.false.
-        etheta=etheta+ethetai
-        if (i.gt.3) gloc(i-3,icg)=gloc(i-3,icg)+wang*dephii
-        if (i.lt.nres) gloc(i-2,icg)=gloc(i-2,icg)+wang*dephii1
-        gloc(nphi+i-2,icg)=gloc(nphi+i-2,icg)+wang*dethetai
-      enddo
       return
       end
 #endif

diff --git a/source/unres/src_MD-M/initialize_p.F b/source/unres/src_MD-M/initialize_p.F
index 6c5f4eb..bdc3269 100644 (file)
--- a/source/unres/src_MD-M/initialize_p.F
+++ b/source/unres/src_MD-M/initialize_p.F
@@ -670,7 +670,7 @@ c      nlen=nres-nnt+1
      & ' ivec_start',ivec_start,' ivec_end',ivec_end,
      & ' iset_start',iset_start,' iset_end',iset_end,
      & ' idihconstr_start',idihconstr_start,' idihconstr_end',
-     &   idihconstr_end
+     &   idihconstr_end,
      & ' ithetaconstr_start',ithetaconstr_start,' ithetaconstr_end',
      &   ithetaconstr_end
 

diff --git a/source/unres/src_MD-M/readrtns_CSA.F b/source/unres/src_MD-M/readrtns_CSA.F
index 286cc57..e10b3ba 100644 (file)
--- a/source/unres/src_MD-M/readrtns_CSA.F
+++ b/source/unres/src_MD-M/readrtns_CSA.F
@@ -856,11 +856,12 @@ C     &   write (iout,*) 'FTORS',ftors
         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 this mean that there is no energy penalty for any angle occuring this can be applied 
+C for generate random conformation but is not implemented in this way
+C      do i=1,nres
+C        thetabound(1,i)=0
+C        thetabound(2,i)=pi
+C      enddo
 C begin reading theta constrains this is quartic constrains allowing to 
 C have smooth second derivative 
       if (with_theta_constr) then
@@ -887,17 +888,17 @@ C E=k*x**4
      &    for_thet_constr(i)
          enddo
         endif
-        do i=1,nthet_constr
+        do i=1,ntheta_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
+C        do i=1,ntheta_constr
+C          ii = itheta_constr(i)
+C          thetabound(1,ii) = phi0(i)-drange(i)
+C          thetabound(2,ii) = phi0(i)+drange(i)
+C        enddo
       endif ! ntheta_constr.gt.0
       endif! with_theta_constr
 C

diff --git a/source/wham/src-M/COMMON.CONTROL b/source/wham/src-M/COMMON.CONTROL
index 6c87389..122ff8d 100644 (file)
--- a/source/wham/src-M/COMMON.CONTROL
+++ b/source/wham/src-M/COMMON.CONTROL
@@ -2,9 +2,10 @@
      & ensembles,constr_dist,symetr
       logical refstr,pdbref,punch_dist,print_rms,caonly,verbose,
      & merge_helices,bxfile,cxfile,histfile,entfile,zscfile,
-     & rmsrgymap,with_dihed_constr,check_conf,histout
+     & rmsrgymap,with_dihed_constr,check_conf,histout,with_theta_constr
       common /cntrl/ iscode,indpdb,refstr,pdbref,outpdb,outmol2,
      & punch_dist,print_rms,caonly,verbose,icomparfunc,pdbint,
      & merge_helices,bxfile,cxfile,histfile,entfile,zscfile,rmsrgymap,
      & ensembles,with_dihed_constr,constr_dist,check_conf,histout,
+     & with_theta_constr,
      &symetr

diff --git a/source/wham/src-M/DIMENSIONS.ZSCOPT b/source/wham/src-M/DIMENSIONS.ZSCOPT
index cfc74a4..dbeffcd 100644 (file)
--- a/source/wham/src-M/DIMENSIONS.ZSCOPT
+++ b/source/wham/src-M/DIMENSIONS.ZSCOPT
@@ -3,7 +3,7 @@
 c Maximum number of structures in the database, energy components, proteins,
 c and structural classes
 c#ifdef JUBL
-      parameter (maxstr=200000,max_ene=21,maxprot=7,maxclass=5000)
+      parameter (maxstr=200000,max_ene=24,maxprot=7,maxclass=5000)
       parameter (maxclass1=10)
 c Maximum number of structures to be dealt with by one processor
       parameter (maxstr_proc=10000)

diff --git a/source/wham/src-M/energy_p_new.F b/source/wham/src-M/energy_p_new.F
index 33e670a..44e5cba 100644 (file)
--- a/source/wham/src-M/energy_p_new.F
+++ b/source/wham/src-M/energy_p_new.F
@@ -67,7 +67,7 @@ cd    print *,'EHPB exitted succesfully.'
 C
 C Calculate the virtual-bond-angle energy.
 C
-      call ebend(ebe)
+      call ebend(ebe,ethetacnstr)
 cd    print *,'Bend energy finished.'
 C
 C Calculate the SC local energy.
@@ -111,7 +111,7 @@ c      write (iout,*) "ft(6)",fact(6)," evdw",evdw," evdw_t",evdw_t
      & +wcorr6*fact(5)*ecorr6+wturn4*fact(3)*eello_turn4
      & +wturn3*fact(2)*eello_turn3+wturn6*fact(5)*eturn6
      & +wel_loc*fact(2)*eel_loc+edihcnstr+wtor_d*fact(2)*etors_d
-     & +wbond*estr+wsccor*fact(1)*esccor
+     & +wbond*estr+wsccor*fact(1)*esccor+ethetacnstr
 #else
       etot=wsc*(evdw+fact(6)*evdw_t)+wscp*evdw2
      & +welec*fact(1)*(ees+evdw1)
@@ -120,7 +120,7 @@ c      write (iout,*) "ft(6)",fact(6)," evdw",evdw," evdw_t",evdw_t
      & +wcorr6*fact(5)*ecorr6+wturn4*fact(3)*eello_turn4
      & +wturn3*fact(2)*eello_turn3+wturn6*fact(5)*eturn6
      & +wel_loc*fact(2)*eel_loc+edihcnstr+wtor_d*fact(2)*etors_d
-     & +wbond*estr+wsccor*fact(1)*esccor
+     & +wbond*estr+wsccor*fact(1)*esccor+ethetacnstr
 #endif
       energia(0)=etot
       energia(1)=evdw
@@ -154,6 +154,7 @@ c      write (iout,*) "ft(6)",fact(6)," evdw",evdw," evdw_t",evdw_t
       energia(19)=esccor
       energia(20)=edihcnstr
       energia(21)=evdw_t
+      energia(24)=ethetacnstr
 c detecting NaNQ
 #ifdef ISNAN
 #ifdef AIX
@@ -270,6 +271,7 @@ C------------------------------------------------------------------------
       esccor=energia(19)
       edihcnstr=energia(20)
       estr=energia(18)
+      ethetacnstr=energia(24)
 #ifdef SPLITELE
       write (iout,10) evdw,wsc,evdw2,wscp,ees,welec*fact(1),evdw1,
      &  wvdwpp,
@@ -278,7 +280,7 @@ C------------------------------------------------------------------------
      &  ecorr,wcorr*fact(3),ecorr5,wcorr5*fact(4),ecorr6,wcorr6*fact(5),
      &  eel_loc,wel_loc*fact(2),eello_turn3,wturn3*fact(2),
      &  eello_turn4,wturn4*fact(3),eello_turn6,wturn6*fact(5),
-     &  esccor,wsccor*fact(1),edihcnstr,ebr*nss,etot
+     &  esccor,wsccor*fact(1),edihcnstr,ethetacnstr,ebr*nss,etot
    10 format (/'Virtual-chain energies:'//
      & 'EVDW=  ',1pE16.6,' WEIGHT=',1pD16.6,' (SC-SC)'/
      & 'EVDW2= ',1pE16.6,' WEIGHT=',1pD16.6,' (SC-p)'/
@@ -300,6 +302,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)'/ 
      & 'ETOT=  ',1pE16.6,' (total)')
 #else
@@ -309,7 +312,7 @@ C------------------------------------------------------------------------
      &  ecorr6,wcorr6*fact(5),eel_loc,wel_loc*fact(2),
      &  eello_turn3,wturn3*fact(2),eello_turn4,wturn4*fact(3),
      &  eello_turn6,wturn6*fact(5),esccor*fact(1),wsccor,
-     &  edihcnstr,ebr*nss,etot
+     &  edihcnstr,ethetacnstr,ebr*nss,etot
    10 format (/'Virtual-chain energies:'//
      & 'EVDW=  ',1pE16.6,' WEIGHT=',1pD16.6,' (SC-SC)'/
      & 'EVDW2= ',1pE16.6,' WEIGHT=',1pD16.6,' (SC-p)'/
@@ -330,6 +333,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)'/ 
      & 'ETOT=  ',1pE16.6,' (total)')
 #endif
@@ -3265,7 +3269,7 @@ c     &      AKSC(j,iti),abond0(j,iti),u(j),j=1,nbi)
       end
 #ifdef CRYST_THETA
 C--------------------------------------------------------------------------
-      subroutine ebend(etheta)
+      subroutine ebend(etheta,ethetacnstr)
 C
 C Evaluate the virtual-bond-angle energy given the virtual-bond dihedral
 C angles gamma and its derivatives in consecutive thetas and gammas.
@@ -3282,6 +3286,7 @@ C
       include 'COMMON.IOUNITS'
       include 'COMMON.NAMES'
       include 'COMMON.FFIELD'
+      include 'COMMON.TORCNSTR'
       common /calcthet/ term1,term2,termm,diffak,ratak,
      & ak,aktc,termpre,termexp,sigc,sig0i,time11,time12,sigcsq,
      & delthe0,sig0inv,sigtc,sigsqtc,delthec,it
@@ -3400,6 +3405,33 @@ c     &    rad2deg*phii,rad2deg*phii1,ethetai
         gloc(nphi+i-2,icg)=wang*(E_theta+E_tc*dthett)
 c 1215   continue
       enddo
+      ethetacnstr=0.0d0
+C      print *,ithetaconstr_start,ithetaconstr_end,"TU"
+      do i=1,ntheta_constr
+        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
+C       if (energy_dec) then
+C        write (iout,'(a6,2i5,4f8.3,2e14.5)') "ethetc",
+C     &    i,itheta,rad2deg*thetiii,
+C     &    rad2deg*theta_constr0(i),  rad2deg*theta_drange(i),
+C     &    rad2deg*difi,0.25d0*for_thet_constr(i)*difi**4,
+C     &    gloc(itheta+nphi-2,icg)
+        endif
+      enddo
 C Ufff.... We've done all this!!! 
       return
       end
@@ -3513,7 +3545,7 @@ C "Thank you" to MAPLE (probably spared one day of hand-differentiation).
       end
 #else
 C--------------------------------------------------------------------------
-      subroutine ebend(etheta)
+      subroutine ebend(etheta,ethetacnstr)
 C
 C Evaluate the virtual-bond-angle energy given the virtual-bond dihedral
 C angles gamma and its derivatives in consecutive thetas and gammas.
@@ -3533,6 +3565,7 @@ C
       include 'COMMON.NAMES'
       include 'COMMON.FFIELD'
       include 'COMMON.CONTROL'
+      include 'COMMON.TORCNSTR'
       double precision coskt(mmaxtheterm),sinkt(mmaxtheterm),
      & cosph1(maxsingle),sinph1(maxsingle),cosph2(maxsingle),
      & sinph2(maxsingle),cosph1ph2(maxdouble,maxdouble),
@@ -3712,6 +3745,34 @@ c        call flush(iout)
 c        gloc(nphi+i-2,icg)=wang*dethetai
         gloc(nphi+i-2,icg)=gloc(nphi+i-2,icg)+wang*dethetai
       enddo
+C now constrains
+      ethetacnstr=0.0d0
+C      print *,ithetaconstr_start,ithetaconstr_end,"TU"
+      do i=1,ntheta_constr
+        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
+C       if (energy_dec) then
+C        write (iout,'(a6,2i5,4f8.3,2e14.5)') "ethetc",
+C     &    i,itheta,rad2deg*thetiii,
+C     &    rad2deg*theta_constr0(i),  rad2deg*theta_drange(i),
+C     &    rad2deg*difi,0.25d0*for_thet_constr(i)*difi**4,
+C     &    gloc(itheta+nphi-2,icg)
+C        endif
+      enddo
       return
       end
 #endif

diff --git a/source/wham/src-M/include_unres/COMMON.TORCNSTR b/source/wham/src-M/include_unres/COMMON.TORCNSTR
index bef2264..845c44d 100644 (file)
--- a/source/wham/src-M/include_unres/COMMON.TORCNSTR
+++ b/source/wham/src-M/include_unres/COMMON.TORCNSTR
@@ -1,6 +1,14 @@
-      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,ithetaconstr_start,
+     & ithetaconstr_end
       double precision phi0(maxdih_constr),drange(maxdih_constr),
-     & ftors(maxdih_constr)
-      common /torcnstr/ phi0,drange,ftors,ndih_constr,idih_constr,
-     &                  ndih_nconstr,idih_nconstr
+     & ftors(maxdih_constr),theta_constr0(maxdih_constr),
+     & theta_drange(maxdih_constr),for_thet_constr(maxdih_constr)
+      common /torcnstr/ phi0,drange,ftors,theta_constr0,theta_drange,
+     & for_thet_constr,
+     &  ndih_constr,idih_constr,
+     &  ndih_nconstr,idih_nconstr,idihconstr_start,idihconstr_end,
+     & ntheta_constr,itheta_constr,ithetaconstr_start,
+     & ithetaconstr_end

diff --git a/source/wham/src-M/initialize_p.F b/source/wham/src-M/initialize_p.F
index c890e09..aea09dd 100644 (file)
--- a/source/wham/src-M/initialize_p.F
+++ b/source/wham/src-M/initialize_p.F
@@ -262,17 +262,19 @@ c-------------------------------------------------------------------------
      &   "EVDW SC-SC","EVDW2 SC-p","EES p-p","ECORR4 ","ECORR5 ",
      &   "ECORR6 ","EELLO ","ETURN3 ","ETURN4 ","ETURN6 ",
      &   "EBE bend","ESC SCloc","ETORS ","ETORSD ","EHPB","EVDWPP",
-     &   "EVDW2_14","ESTR","ESCCOR","EDIHC","EVDW_T"/
+     &   "EVDW2_14","ESTR","ESCCOR","EDIHC","EVDW_T","ELIPTRAN",
+     &   "EAFM","ETHETC"/
       data wname /
      &   "WSC","WSCP","WELEC","WCORR","WCORR5","WCORR6","WEL_LOC",
      &   "WTURN3","WTURN4","WTURN6","WANG","WSCLOC","WTOR","WTORD",
-     &   "WHPB","WVDWPP","WSCP14","WBOND","WSCCOR","WDIHC","WSC"/
+     &   "WHPB","WVDWPP","WSCP14","WBOND","WSCCOR","WDIHC","WSC",
+     &   "WLIPTRAN","WAFM","WTHETC"/
       data ww0 /1.0d0,1.0d0,1.0d0,1.0d0,1.0d0,1.0d0,1.0d0,1.0d0,1.0d0,
      &    1.0d0,1.0d0,1.0d0,1.0d0,1.0d0,1.0d0,1.0d0,0.4d0,1.0d0,1.0d0,
-     &    0.0d0,0.0/
-      data nprint_ene /21/
+     &    0.0d0,0.0,0.0d0,0.0d0,0.0d0/
+      data nprint_ene /22/
       data print_order /1,2,3,18,11,12,13,14,4,5,6,7,8,9,10,19,
-     &  16,15,17,20,21/
+     &  16,15,17,20,21,24,22,23/
       end 
 c---------------------------------------------------------------------------
       subroutine init_int_table

diff --git a/source/wham/src-M/molread_zs.F b/source/wham/src-M/molread_zs.F
index 061ae05..2b9fe8c 100644 (file)
--- a/source/wham/src-M/molread_zs.F
+++ b/source/wham/src-M/molread_zs.F
@@ -94,7 +94,43 @@ C        write (iout,*) 'FTORS',ftors
       endif
 
       endif
-
+      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 
+C        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
+C        endif
+        do i=1,ntheta_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
+C        do i=1,ntheta_constr
+C          ii = itheta_constr(i)
+C          thetabound(1,ii) = phi0(i)-drange(i)
+C          thetabound(2,ii) = phi0(i)+drange(i)
+C        enddo
+      endif ! ntheta_constr.gt.0
+      endif! with_theta_constr
       nnt=1
       nct=nres
       if (itype(1).eq.ntyp1) nnt=2