SAXS cutoff correction

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 1be4376..ff2d579 100644 (file)
--- a/source/unres/src_MD-M/energy_p_new_barrier.F
+++ b/source/unres/src_MD-M/energy_p_new_barrier.F
@@ -11050,12 +11050,16 @@ C      print *,'AFM',Eafmforce,totTafm*velAFMconst,dist
       end
 
 c----------------------------------------------------------------------------
-      double precision function sscale2(r,r_cut,rlamb)
-      double precision r,gamm
-      if(r.lt.r_cut-rlamb) then
+      double precision function sscale2(r,r_cut,r0,rlamb)
+      implicit none
+      double precision r,gamm,r_cut,r0,rlamb,rr
+      rr = dabs(r-r0)
+c      write (2,*) "r",r," r_cut",r_cut," r0",r0," rlamb",rlamb
+c      write (2,*) "rr",rr
+      if(rr.lt.r_cut-rlamb) then
         sscale2=1.0d0
-      else if(r.le.r_cut.and.r.ge.r_cut-rlamb) then
-        gamm=(r-(r_cut-rlamb))/rlamb
+      else if(rr.le.r_cut.and.rr.ge.r_cut-rlamb) then
+        gamm=(rr-(r_cut-rlamb))/rlamb
         sscale2=1.0d0+gamm*gamm*(2*gamm-3.0d0)
       else
         sscale2=0d0
@@ -11063,13 +11067,16 @@ c----------------------------------------------------------------------------
       return
       end
 C-----------------------------------------------------------------------
-      double precision function sscalgrad2(r,r_cut,rlamb)
-      double precision r,gamm
-      if(r.lt.r_cut-rlamb) then
+      double precision function sscalgrad2(r,r_cut,r0,rlamb)
+      implicit none
+      double precision r,gamm,r_cut,r0,rlamb,rr
+      rr = dabs(r-r0)
+      if(rr.lt.r_cut-rlamb) then
         sscalgrad2=0.0d0
-      else if(r.le.r_cut.and.r.ge.r_cut-rlamb) then
-        gamm=(r-(r_cut-rlamb))/rlamb
+      else if(rr.le.r_cut.and.rr.ge.r_cut-rlamb) then
+        gamm=(rr-(r_cut-rlamb))/rlamb
         sscalgrad2=gamm*(6*gamm-6.0d0)/rlamb
+        if (r.lt.r0) sscalgrad2=-sscalgrad2
       else
         sscalgrad2=0.0d0
       endif
@@ -11207,11 +11214,14 @@ c SC SC
            sigma2CACA=scal_rad**2*0.25d0/
      &        (restok(itype(j))**2+restok(itype(i))**2)
            rrr = 2.0d0/dsqrt(sigma2CACA)
-           sss2 = sscale2(dijCACA,rrr,0.3d0)
-           if (sss2.ne.0.0d0) then
-           ssgrad2 = sscalgrad2(dijCACA,rrr,0.3d0)
            do k=1,nsaxs
              dk = distsaxs(k)
+c             write (2,*) "ijk",i,j,k
+             sss2 = sscale2(dijCACA,rrr,dk,0.3d0)
+             if (sss2.ne.0.0d0) then
+             ssgrad2 = sscalgrad2(dijCACA,rrr,dk,0.3d0)
+             if (energy_dec) write(iout,'(a4,3i5,4f10.4)') 
+     &          'saxs',i,j,k,dijCACA,rrr,dk,sss2
              expCACA = dexp(-0.5d0*sigma2CACA*(dijCACA-dk)**2)*sss2
              Pcalc(k) = Pcalc(k)+expCACA
 #ifdef DEBUG
@@ -11225,8 +11235,8 @@ c CA CA
                PgradC(k,l,i) = PgradC(k,l,i)-aux
                PgradC(k,l,j) = PgradC(k,l,j)+aux
              enddo ! l
+             endif ! sss
            enddo ! k
-           endif ! sss
 #endif
          enddo ! j
        enddo ! iint