From: Cezary Czaplewski <czarek@chem.univ.gda.pl>
Date: Tue, 2 Jan 2018 11:28:18 +0000 (+0100)
Subject: SAXS cutoff
X-Git-Url: http://mmka.chem.univ.gda.pl/gitweb/?p=unres.git;a=commitdiff_plain;h=59f6a39e24764975c960c93afbc2b1cb930c5a35

SAXS cutoff
---

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 7a158ed..1be4376 100644
--- a/source/unres/src_MD-M/energy_p_new_barrier.F
+++ b/source/unres/src_MD-M/energy_p_new_barrier.F
@@ -11050,6 +11050,32 @@ 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
+        sscale2=1.0d0
+      else if(r.le.r_cut.and.r.ge.r_cut-rlamb) then
+        gamm=(r-(r_cut-rlamb))/rlamb
+        sscale2=1.0d0+gamm*gamm*(2*gamm-3.0d0)
+      else
+        sscale2=0d0
+      endif
+      return
+      end
+C-----------------------------------------------------------------------
+      double precision function sscalgrad2(r,r_cut,rlamb)
+      double precision r,gamm
+      if(r.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
+        sscalgrad2=gamm*(6*gamm-6.0d0)/rlamb
+      else
+        sscalgrad2=0.0d0
+      endif
+      return
+      end
+c----------------------------------------------------------------------------
       subroutine e_saxs(Esaxs_constr)
       implicit none
       include 'DIMENSIONS'
@@ -11084,6 +11110,7 @@ c
      & sigma2CACA,sigma2CASC,sigma2SCCA,sigma2SCSC,expCACA,expCASC,
      & expSCCA,expSCSC,CASCgrad,SCCAgrad,SCSCgrad,aux,auxC,auxC1,
      & auxX,auxX1,CACAgrad,Cnorm
+      double precision sss2,ssgrad2,rrr,sscalgrad2,sscale2
       double precision dist
       external dist
 c  SAXS restraint penalty function
@@ -11179,14 +11206,19 @@ c SC SC
            dijCACA=dist(i,j)
            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)
-             expCACA = dexp(-0.5d0*sigma2CACA*(dijCACA-dk)**2)
+             expCACA = dexp(-0.5d0*sigma2CACA*(dijCACA-dk)**2)*sss2
              Pcalc(k) = Pcalc(k)+expCACA
 #ifdef DEBUG
              write(iout,*) "i j k Pcalc",i,j,Pcalc(k)
 #endif
-             CACAgrad = sigma2CACA*(dijCACA-dk)*expCACA
+             CACAgrad = sigma2CACA*(dijCACA-dk)*expCACA+
+     &             ssgrad2*expCACA/sss2
              do l=1,3
 c CA CA 
                aux = CACAgrad*(C(l,j)-C(l,i))/dijCACA
@@ -11194,6 +11226,7 @@ c CA CA
                PgradC(k,l,j) = PgradC(k,l,j)+aux
              enddo ! l
            enddo ! k
+           endif ! sss
 #endif
          enddo ! j
        enddo ! iint