cluster wham SACS cutoff

diff --git a/source/cluster/wham/src-M/COMMON.CONTROL b/source/cluster/wham/src-M/COMMON.CONTROL
index 1a6299f..9bd6fd4 100644 (file)
--- a/source/cluster/wham/src-M/COMMON.CONTROL
+++ b/source/cluster/wham/src-M/COMMON.CONTROL
@@ -6,7 +6,8 @@
       real*8 waga_homology
       real*8 waga_dist, waga_angle, waga_theta, waga_d, dist_cut,
      &  dist2_cut
-      real*8 Psaxs(maxsaxs),distsaxs(maxsaxs),CSAXS(3,maxsaxs),scal_rad
+      real*8 Psaxs(maxsaxs),distsaxs(maxsaxs),CSAXS(3,maxsaxs),
+     & scal_rad,wsaxs0,saxs_cutoff
       logical refstr,pdbref,punch_dist,print_dist,caonly,lside,
      & lprint_cart,lprint_int,from_cart,efree,from_bx,from_cx,
      & with_dihed_constr,with_theta_constr,out1file,
@@ -24,4 +25,6 @@
      & waga_dist,waga_angle,waga_theta,waga_d,dist_cut,dist2_cut,
      & iset,ihset,l_homo(max_template,maxdim),
      & print_homology_restraints,print_homology_models
-      common /saxsretr/ Psaxs,distsaxs,csaxs,scal_rad,nsaxs,saxs_mode
+      common /saxsretr/ Psaxs,distsaxs,csaxs,Wsaxs0,scal_rad,
+     & saxs_cutoff,nsaxs,saxs_mode
+

diff --git a/source/cluster/wham/src-M/energy_p_new.F b/source/cluster/wham/src-M/energy_p_new.F
index 55cc7a1..2807db0 100644 (file)
--- a/source/cluster/wham/src-M/energy_p_new.F
+++ b/source/cluster/wham/src-M/energy_p_new.F
@@ -9563,8 +9563,42 @@ C        sscagrad=0.0d0
 C      endif
       return
       end
-
+c----------------------------------------------------------------------------
+      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(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
+      endif
+      return
+      end
 C-----------------------------------------------------------------------
+      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(rr.le.r_cut.and.rr.ge.r_cut-rlamb) then
+        gamm=(rr-(r_cut-rlamb))/rlamb
+        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
+      return
+      end
 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
       subroutine Eliptransfer(eliptran)
       implicit real*8 (a-h,o-z)
@@ -9715,6 +9749,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
@@ -9738,8 +9773,10 @@ c  SAXS restraint penalty function
         enddo
       enddo
       do i=iatsc_s,iatsc_e
+       if (itype(i).eq.ntyp1) cycle
        do iint=1,nint_gr(i)
          do j=istart(i,iint),iend(i,iint)
+           if (itype(j).eq.ntyp1) cycle
 #ifdef ALLSAXS
            dijCACA=dist(i,j)
            dijCASC=dist(i,j+nres)
@@ -9808,21 +9845,42 @@ c SC SC
            dijCACA=dist(i,j)
            sigma2CACA=scal_rad**2*0.25d0/
      &        (restok(itype(j))**2+restok(itype(i))**2)
+
+           IF (saxs_cutoff.eq.0) THEN
            do k=1,nsaxs
              dk = distsaxs(k)
              expCACA = dexp(-0.5d0*sigma2CACA*(dijCACA-dk)**2)
              Pcalc(k) = Pcalc(k)+expCACA
+             CACAgrad = sigma2CACA*(dijCACA-dk)*expCACA
+             do l=1,3
+               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
+             enddo ! l
+           enddo ! k
+           ELSE
+           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.eq.0.0d0) cycle
+             ssgrad2 = sscalgrad2(dijCACA,rrr,dk,0.3d0)
+             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
            enddo ! k
+           ENDIF
 #endif
          enddo ! j
        enddo ! iint
@@ -9869,9 +9927,10 @@ c CA CA
       do k=1,nsaxs
         Cnorm = Cnorm + Pcalc(k)
       enddo
-      Esaxs_constr = dlog(Cnorm)
+      Esaxs_constr = dlog(Cnorm)-wsaxs0
       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
@@ -9886,7 +9945,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)

diff --git a/source/cluster/wham/src-M/readrtns.F b/source/cluster/wham/src-M/readrtns.F
index e5c31cf..60b1ec2 100644 (file)
--- a/source/cluster/wham/src-M/readrtns.F
+++ b/source/cluster/wham/src-M/readrtns.F
@@ -131,8 +131,9 @@ C long axis of side chain
       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
       if (min_var) iopt=1
       return
       end
@@ -997,6 +998,11 @@ c SAXS distance distribution
       do i=1,nsaxs
         write (iout,'(f8.2,e15.5)') distsaxs(i),Psaxs(i)
       enddo
+      Wsaxs0=0.0d0
+      do i=1,nsaxs
+        Wsaxs0=Wsaxs0-Psaxs(i)*dlog(Psaxs(i))
+      enddo
+      write (iout,*) "Wsaxs0",Wsaxs0
       else
 c SAXS "spheres".
       do i=1,nsaxs