debuging WHAM for constrains DEBUG ON

diff --git a/source/wham/src-M/DIMENSIONS b/source/wham/src-M/DIMENSIONS
index 154f3df..e9dac01 100644 (file)
--- a/source/wham/src-M/DIMENSIONS
+++ b/source/wham/src-M/DIMENSIONS
@@ -13,8 +13,8 @@ C Max. number of coarse-grain processors
 c      parameter (max_cg_procs=maxprocs)
 C Max. number of AA residues
       integer maxres
-c      parameter (maxres=250)
-      parameter (maxres=100)
+      parameter (maxres=130)
+C      parameter (maxres=100)
 C Appr. max. number of interaction sites
       integer maxres2
       parameter (maxres2=2*maxres)

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)

diff --git a/source/wham/src-M/gnmr1.f b/source/wham/src-M/gnmr1.f
index 905e746..8bfc43a 100644 (file)
--- a/source/wham/src-M/gnmr1.f
+++ b/source/wham/src-M/gnmr1.f
@@ -41,3 +41,33 @@ c-------------------------------------------------------------------------------
       return
       end
 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/wham/src-M/include_unres/COMMON.SBRIDGE b/source/wham/src-M/include_unres/COMMON.SBRIDGE
index bea57e4..028f9ae 100644 (file)
--- a/source/wham/src-M/include_unres/COMMON.SBRIDGE
+++ b/source/wham/src-M/include_unres/COMMON.SBRIDGE
@@ -2,9 +2,10 @@
       integer ns,nss,nfree,iss
       common /sbridge/ ss_depth,ebr,d0cm,akcm,akth,akct,v1ss,v2ss,v3ss,
      & ns,nss,nfree,iss(maxss)
-      double precision dhpb,dhpb1,forcon
+      double precision dhpb,dhpb1,forcon,fordepth
       integer ihpb,jhpb,nhpb,idssb,jdssb,ibecarb
       common /links/ dhpb(maxdim),dhpb1(maxdim),forcon(maxdim),
+     & fordepth(maxdim),
      & ihpb(maxdim),jhpb(maxdim),nhpb
       double precision weidis
       common /restraints/ weidis

diff --git a/source/wham/src-M/read_dist_constr.F b/source/wham/src-M/read_dist_constr.F
index 2e9ba3b..a97f797 100644 (file)
--- a/source/wham/src-M/read_dist_constr.F
+++ b/source/wham/src-M/read_dist_constr.F
@@ -13,9 +13,12 @@
       character*500 controlcard
       logical lprn /.true./
       write (iout,*) "Calling read_dist_constr"
-      write (iout,*) "nres",nres," nstart_sup",nstart_sup," nsup",nsup
+C      write (iout,*) "nres",nres," nstart_sup",nstart_sup," nsup",nsup
+C      call flush(iout)
+      write(iout,*) "TU sie wywalam?"
+      call card_concat(controlcard,.false.)
+      write (iout,*) controlcard
       call flush(iout)
-      call card_concat(controlcard)
       call readi(controlcard,"NFRAG",nfrag_,0)
       call readi(controlcard,"NPAIR",npair_,0)
       call readi(controlcard,"NDIST",ndist_,0)
@@ -100,16 +103,20 @@ c        write (iout,*) i,ifrag_(1,i),ifrag_(2,i),wfrag_(i)
         read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),dhpb(i),dhpb1(i),
      &     ibecarb(i),forcon(nhpb+1),fordepth(nhpb+1)
         fordepth(nhpb+1)=fordepth(nhpb+1)/forcon(nhpb+1)
+C        write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",
+C     &     nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
         else
         read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),forcon(nhpb+1)
+        endif
         if (forcon(nhpb+1).gt.0.0d0) then
           nhpb=nhpb+1
           dhpb(nhpb)=dist(ihpb(nhpb),jhpb(nhpb))
           write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",
      &     nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
         endif
-      endif
+C      endif
       enddo
+      call hpb_partition
       call flush(iout)
       return
       end

diff --git a/source/wham/src-M/wham_multparm.F b/source/wham/src-M/wham_multparm.F
index e80a39d..de73250 100644 (file)
--- a/source/wham/src-M/wham_multparm.F
+++ b/source/wham/src-M/wham_multparm.F
@@ -58,6 +58,7 @@ c NaNQ initialization
       call flush(iout)
       call molread(*10)
       call flush(iout)
+      if (constr_dist.gt.0) call read_dist_constr
 #ifdef MPI 
       write (iout,*) "Calling proc_groups"
       call proc_groups
@@ -88,8 +89,8 @@ c NaNQ initialization
         call read_ref_structure(*10)
         call proc_cont
         call fragment_list
-        if (constr_dist.gt.0) call read_dist_constr
       endif
+C      if (constr_dist.gt.0) call read_dist_constr
       write (iout,*) "Begin read_database"
       call flush(iout)
       call read_database(*10)