wsaxs0 and map

[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 ff2d579..526b4ab 100644 (file)
--- 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,30 +11216,44 @@ 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)
+
+           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.ne.0.0d0) then
+             if (sss2.eq.0.0d0) cycle
              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
+           ENDIF
 #endif
          enddo ! j
        enddo ! iint
@@ -11283,9 +11300,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
@@ -11300,7 +11318,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)