From: Cezary Czaplewski <czarek@chem.univ.gda.pl>
Date: Wed, 3 Jan 2018 14:25:01 +0000 (+0100)
Subject: SAXS cutoff correct gradiend and cutoff scaling
X-Git-Url: http://mmka.chem.univ.gda.pl/gitweb/?p=unres.git;a=commitdiff_plain;h=024c09f85d96670bfddf91ec2305f8a5e2532adc

SAXS cutoff correct gradiend and cutoff scaling
---

diff --git a/source/unres/src_MD-M/COMMON.CONTROL b/source/unres/src_MD-M/COMMON.CONTROL
index 4de9697..6405131 100644
--- a/source/unres/src_MD-M/COMMON.CONTROL
+++ b/source/unres/src_MD-M/COMMON.CONTROL
@@ -3,7 +3,7 @@
      & constr_homology,homol_nset,nsaxs,saxs_mode
       real*8 waga_homology
       real*8 waga_dist, waga_angle, waga_theta, waga_d, dist_cut,
-     &  dist2_cut, scal_rad
+     &  dist2_cut, scal_rad, saxs_cutoff
       real*8 Psaxs(maxsaxs),distsaxs(maxsaxs),CSAXS(3,maxsaxs)
       logical minim,refstr,pdbref,outpdb,outmol2,overlapsc,energy_dec,
      &                 sideadd,lsecondary,read_cart,unres_pdb,
@@ -19,6 +19,7 @@
      & constr_homology,homol_nset,read2sigma,start_from_model
       common /homol/ waga_homology(maxprocs/20),
      & waga_dist, waga_angle, waga_theta, waga_d, dist_cut,dist2_cut
-      common /saxsretr/ Psaxs,distsaxs,csaxs,scal_rad,nsaxs,saxs_mode
+      common /saxsretr/ Psaxs,distsaxs,csaxs,scal_rad,saxs_cutoff,
+     & nsaxs,saxs_mode
 C... minim = .true. means DO minimization.
 C... energy_dec = .true. means print energy decomposition matrix
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 1517e8f..f95ca0c 100644
--- a/source/unres/src_MD-M/energy_p_new_barrier.F
+++ b/source/unres/src_MD-M/energy_p_new_barrier.F
@@ -11075,8 +11075,11 @@ C-----------------------------------------------------------------------
         sscalgrad2=0.0d0
       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
+        if (r.ge.r0) then
+          sscalgrad2=gamm*(6*gamm-6.0d0)/rlamb
+        else
+          sscalgrad2=-gamm*(6*gamm-6.0d0)/rlamb
+        endif
       else
         sscalgrad2=0.0d0
       endif
@@ -11213,27 +11216,27 @@ 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)
+           rrr = saxs_cutoff*2.0d0/dsqrt(sigma2CACA)
            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
+             if (energy_dec) write(iout,'(a4,3i5,5f10.4)') 
+     &          'saxs',i,j,k,dijCACA,rrr,dk,sss2,ssgrad2
              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
-               PgradC(k,l,i) = PgradC(k,l,i)-aux
-               PgradC(k,l,j) = PgradC(k,l,j)+aux
+               PgradC(k,l,i) = PgradC(k,l,i)+aux
+               PgradC(k,l,j) = PgradC(k,l,j)-aux
              enddo ! l
              endif ! sss
            enddo ! k
diff --git a/source/unres/src_MD-M/readrtns_CSA.F b/source/unres/src_MD-M/readrtns_CSA.F
index 5b71db2..70bb647 100644
--- a/source/unres/src_MD-M/readrtns_CSA.F
+++ b/source/unres/src_MD-M/readrtns_CSA.F
@@ -102,8 +102,9 @@ C Set up the time limit (caution! The time must be input in minutes!)
       call readi(controlcard,'NSAXS',nsaxs,0)
       call readi(controlcard,'SAXS_MODE',saxs_mode,0)
       call reada(controlcard,'SCAL_RAD',scal_rad,1.0d0)
+      call reada(controlcard,'SAXS_CUTOFF',saxs_cutoff,1.0d0)
       write (iout,*) "Number of SAXS restraints",NSAXS," SAXS_MODE",
-     &   SAXS_MODE," SCAL_RAD",scal_rad
+     &   SAXS_MODE," SCAL_RAD",scal_rad,"SAXS_CUTOFF",saxs_cutoff
       call readi(controlcard,'CONSTR_HOMOL',constr_homology,0)
       call readi(controlcard,'SYM',symetr,1)
       call reada(controlcard,'TIMLIM',timlim,960.0D0) ! default 16 hours