debuging WHAM for constrains DEBUG 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 6e19387..f184eed 100644 (file)
--- a/source/wham/src-M/energy_p_new.F
+++ b/source/wham/src-M/energy_p_new.F
@@ -2947,10 +2947,13 @@ C
       include 'COMMON.DERIV'
       include 'COMMON.VAR'
       include 'COMMON.INTERACT'
+      include 'COMMON.CONTROL'
+      include 'COMMON.IOUNITS'
       dimension ggg(3)
       ehpb=0.0D0
 cd    print *,'edis: nhpb=',nhpb,' fbr=',fbr
 cd    print *,'link_start=',link_start,' link_end=',link_end
+      write(iout,*) link_end, "link_end"
       if (link_end.eq.0) return
       do i=link_start,link_end
 C If ihpb(i) and jhpb(i) > NRES, this is a SC-SC distance, otherwise a
@@ -2969,21 +2972,27 @@ C 24/11/03 AL: SS bridges handled separately because of introducing a specific
 C    distance and angle dependent SS bond potential.
 C        if (ii.gt.nres .and. iabs(itype(iii)).eq.1 .and. 
 C     & iabs(itype(jjj)).eq.1) then
-
+       write(iout,*) constr_dist,"const"
        if (.not.dyn_ss .and. i.le.nss) then
          if (ii.gt.nres .and. iabs(itype(iii)).eq.1 .and.
      & iabs(itype(jjj)).eq.1) then
           call ssbond_ene(iii,jjj,eij)
           ehpb=ehpb+2*eij
-           endif
+           endif !ii.gt.neres
         else if (ii.gt.nres .and. jj.gt.nres) then
 c Restraints from contact prediction
           dd=dist(ii,jj)
           if (constr_dist.eq.11) then
+C            ehpb=ehpb+fordepth(i)**4.0d0
+C     &          *rlornmr1(dd,dhpb(i),dhpb1(i),forcon(i))
             ehpb=ehpb+fordepth(i)**4.0d0
      &          *rlornmr1(dd,dhpb(i),dhpb1(i),forcon(i))
             fac=fordepth(i)**4.0d0
      &          *rlornmr1prim(dd,dhpb(i),dhpb1(i),forcon(i))/dd
+            write(iout,*) ehpb,"atu?"
+C            ehpb,"tu?"
+C            fac=fordepth(i)**4.0d0
+C     &          *rlornmr1prim(dd,dhpb(i),dhpb1(i),forcon(i))/dd
            else
           if (dhpb1(i).gt.0.0d0) then
             ehpb=ehpb+2*forcon(i)*gnmr1(dd,dhpb(i),dhpb1(i))
@@ -3002,8 +3011,8 @@ C
 C Evaluate gradient.
 C
             fac=waga*rdis/dd
-          endif
-          endif
+          endif !end dhpb1(i).gt.0
+          endif !end const_dist=11
           do j=1,3
             ggg(j)=fac*(c(j,jj)-c(j,ii))
           enddo
@@ -3015,7 +3024,26 @@ C
             ghpbc(k,jjj)=ghpbc(k,jjj)+ggg(k)
             ghpbc(k,iii)=ghpbc(k,iii)-ggg(k)
           enddo
-        else
+        else !ii.gt.nres
+          write(iout,*) "before"
+          dd=dist(ii,jj)
+          write(iout,*) "after",dd
+          if (constr_dist.eq.11) then
+            ehpb=ehpb+fordepth(i)**4.0d0
+     &          *rlornmr1(dd,dhpb(i),dhpb1(i),forcon(i))
+            fac=fordepth(i)**4.0d0
+     &          *rlornmr1prim(dd,dhpb(i),dhpb1(i),forcon(i))/dd
+C            ehpb=ehpb+fordepth(i)**4*rlornmr1(dd,dhpb(i),dhpb1(i))
+C            fac=fordepth(i)**4*rlornmr1prim(dd,dhpb(i),dhpb1(i))/dd
+C            print *,ehpb,"tu?"
+            write(iout,*) ehpb,"btu?",dd,dhpb(i),dhpb1(i),fordepth(i)
+           else   
+          if (dhpb1(i).gt.0.0d0) then
+            ehpb=ehpb+2*forcon(i)*gnmr1(dd,dhpb(i),dhpb1(i))
+            fac=forcon(i)*gnmr1prim(dd,dhpb(i),dhpb1(i))/dd
+c            write (iout,*) "alph nmr",
+c     &        dd,2*forcon(i)*gnmr1(dd,dhpb(i),dhpb1(i))
+          else
             rdis=dd-dhpb(i)
 C Get the force constant corresponding to this distance.
             waga=forcon(i)