projects / unres.git / blobdiff
? search:
summary | shortlog | log | commit | commitdiff | tree
raw | inline | side by side
1) debug of Lorentzian like constrains in cluster 2) introduction of constrains on...
[unres.git] / source / wham / src-M / energy_p_new.F
diff --git a/source/wham/src-M/energy_p_new.F b/source/wham/src-M/energy_p_new.F
index 8cf29b2..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
@@ -4478,15 +4539,16 @@ 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)
+C        write (iout,'(a6,2i5,2f8.3,2e14.5)') "edih",
+C     &    i,itori,rad2deg*phii,
+C     &    rad2deg*difi,0.25d0*ftors(i)*difi**4,gloc(itori-3,icg)
       enddo
 !      write (iout,*) 'edihcnstr',edihcnstr
       return
@@ -4573,21 +4635,24 @@ c       write (iout,*) 'i=',i,' gloc=',gloc(i-3,icg)
         edihi=0.0d0
         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
-          edihi=0.25d0*ftors*difi**4
+          edihcnstr=edihcnstr+0.25d0*ftors(i)*difi**4
+          gloc(itori-3,icg)=gloc(itori-3,icg)+ftors(i)*difi**3
+          edihi=0.25d0*ftors(i)*difi**4
         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
-          edihi=0.25d0*ftors*difi**4
+          edihcnstr=edihcnstr+0.25d0*ftors(i)*difi**4
+          gloc(itori-3,icg)=gloc(itori-3,icg)+ftors(i)*difi**3
+          edihi=0.25d0*ftors(i)*difi**4
         else
           difi=0.0d0
         endif
+        write (iout,'(a6,2i5,2f8.3,2e14.5)') "edih",
+     &    i,itori,rad2deg*phii,
+     &    rad2deg*difi,0.25d0*ftors(i)*difi**4
 c        write (iout,'(2i5,4f10.5,e15.5)') i,itori,phii,phi0(i),difi,
 c     &    drange(i),edihi
 !        write (iout,'(2i5,2f8.3,2e14.5)') i,itori,rad2deg*phii,
-!     &    rad2deg*difi,0.25d0*ftors*difi**4,gloc(itori-3,icg)
+!     &    rad2deg*difi,0.25d0*ftors(i)*difi**4,gloc(itori-3,icg)
       enddo
 !      write (iout,*) 'edihcnstr',edihcnstr
       return
UNRESPACK 3.x