projects / unres.git / blobdiff
? search:
summary | shortlog | log | commit | commitdiff | tree
raw | inline | side by side
first introduction of valence constrains - not working yet
[unres.git] / source / unres / src_MD-M / energy_p_new_barrier.F
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 745f7d5..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)')
@@ -4090,6 +4095,7 @@ C
       do i=1,3
        ggg(i)=0.0d0
       enddo
+C      write (iout,*) ,"link_end",link_end,constr_dist
 cd      write(iout,*)'edis: nhpb=',nhpb,' fbr=',fbr
 cd      write(iout,*)'link_start=',link_start,' link_end=',link_end
       if (link_end.eq.0) return
@@ -4129,6 +4135,8 @@ c Restraints from contact prediction
      &          *rlornmr1(dd,dhpb(i),dhpb1(i),forcon(i))
             fac=fordepth(i)**4.0d0
      &          *rlornmr1prim(dd,dhpb(i),dhpb1(i),forcon(i))/dd
+          if (energy_dec) write (iout,'(a6,2i5,3f8.3)') "edisl",ii,jj,
+     &    ehpb,fordepth(i),dd
            else
           if (dhpb1(i).gt.0.0d0) then
             ehpb=ehpb+2*forcon(i)*gnmr1(dd,dhpb(i),dhpb1(i))
@@ -4169,6 +4177,8 @@ C target distance.
      &           *rlornmr1(dd,dhpb(i),dhpb1(i),forcon(i))
             fac=fordepth(i)**4.0d0
      &           *rlornmr1prim(dd,dhpb(i),dhpb1(i),forcon(i))/dd
+          if (energy_dec) write (iout,'(a6,2i5,3f8.3)') "edisl",ii,jj,
+     &    ehpb,fordepth(i),dd
            else   
           if (dhpb1(i).gt.0.0d0) then
             ehpb=ehpb+2*forcon(i)*gnmr1(dd,dhpb(i),dhpb1(i))
@@ -4515,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
@@ -4827,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) 
@@ -5658,12 +5720,12 @@ c       write (iout,*) 'i=',i,' gloc=',gloc(i-3,icg)
         difi=phii-phi0(i)
         if (difi.gt.drange(i)) then
           difi=difi-drange(i)
-          edihcnstr=edihcnstr+0.25d0*ftors*difi**4
-          gloc(itori-3,icg)=gloc(itori-3,icg)+ftors*difi**3
+          edihcnstr=edihcnstr+0.25d0*ftors(i)*difi**4
+          gloc(itori-3,icg)=gloc(itori-3,icg)+ftors(i)*difi**3
         else if (difi.lt.-drange(i)) then
           difi=difi+drange(i)
-          edihcnstr=edihcnstr+0.25d0*ftors*difi**4
-          gloc(itori-3,icg)=gloc(itori-3,icg)+ftors*difi**3
+          edihcnstr=edihcnstr+0.25d0*ftors(i)**difi**4
+          gloc(itori-3,icg)=gloc(itori-3,icg)+ftors(i)*difi**3
         endif
 !        write (iout,'(2i5,2f8.3,2e14.5)') i,itori,rad2deg*phii,
 !     &    rad2deg*difi,0.25d0*ftors*difi**4,gloc(itori-3,icg)
@@ -5762,18 +5824,21 @@ c      do i=1,ndih_constr
         difi=pinorm(phii-phi0(i))
         if (difi.gt.drange(i)) then
           difi=difi-drange(i)
-          edihcnstr=edihcnstr+0.25d0*ftors*difi**4
-          gloc(itori-3,icg)=gloc(itori-3,icg)+ftors*difi**3
+          edihcnstr=edihcnstr+0.25d0*ftors(i)*difi**4
+          gloc(itori-3,icg)=gloc(itori-3,icg)+ftors(i)*difi**3
         else if (difi.lt.-drange(i)) then
           difi=difi+drange(i)
-          edihcnstr=edihcnstr+0.25d0*ftors*difi**4
-          gloc(itori-3,icg)=gloc(itori-3,icg)+ftors*difi**3
+          edihcnstr=edihcnstr+0.25d0*ftors(i)*difi**4
+          gloc(itori-3,icg)=gloc(itori-3,icg)+ftors(i)*difi**3
         else
           difi=0.0
         endif
-cd        write (iout,'(2i5,4f8.3,2e14.5)') i,itori,rad2deg*phii,
-cd     &    rad2deg*phi0(i),  rad2deg*drange(i),
-cd     &    rad2deg*difi,0.25d0*ftors*difi**4,gloc(itori-3,icg)
+       if (energy_dec) then
+        write (iout,'(a6,2i5,4f8.3,2e14.5)') "edihc",
+     &    i,itori,rad2deg*phii,
+     &    rad2deg*phi0(i),  rad2deg*drange(i),
+     &    rad2deg*difi,0.25d0*ftors(i)*difi**4,gloc(itori-3,icg)
+        endif
       enddo
 cd       write (iout,*) 'edihcnstr',edihcnstr
       return
UNRESPACK 3.x