SAXS cutoff correction

[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 1be4376..1517e8f 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
@@ -11275,7 +11285,8 @@ c CA CA
       enddo
       Esaxs_constr = dlog(Cnorm)
       do k=1,nsaxs
-        Esaxs_constr = Esaxs_constr - Psaxs(k)*dlog(Pcalc(k)) 
+        if (Pcalc(k).gt.0.0d0) 
+     &  Esaxs_constr = Esaxs_constr - Psaxs(k)*dlog(Pcalc(k)) 
 #ifdef DEBUG
         write (iout,*) "k",k," Esaxs_constr",Esaxs_constr
 #endif
@@ -11290,7 +11301,8 @@ c CA CA
           auxX=0.0d0
           auxX1=0.d0 
           do k=1,nsaxs
-            auxC  = auxC +Psaxs(k)*PgradC(k,l,i)/Pcalc(k)
+            if (Pcalc(k).gt.0) 
+     &      auxC  = auxC +Psaxs(k)*PgradC(k,l,i)/Pcalc(k)
             auxC1 = auxC1+PgradC(k,l,i)
 #ifdef ALLSAXS
             auxX  = auxX +Psaxs(k)*PgradX(k,l,i)/Pcalc(k)