diff --git a/Hadrons/Modules/MScalarSUN/StochFreeField.hpp b/Hadrons/Modules/MScalarSUN/StochFreeField.hpp
index 1e4d2265..af55b115 100644
--- a/Hadrons/Modules/MScalarSUN/StochFreeField.hpp
+++ b/Hadrons/Modules/MScalarSUN/StochFreeField.hpp
@@ -44,7 +44,8 @@ class StochFreeFieldPar: Serializable
 public:
     GRID_SERIALIZABLE_CLASS_MEMBERS(StochFreeFieldPar,
                                     double, m2,
-                                    double, g);
+                                    double, g,
+                                    double, smearing);
 };
 
 template <typename SImpl>
@@ -111,6 +112,7 @@ void TStochFreeField<SImpl>::setup(void)
     if (!env().hasCreatedObject("_" + getName() + "_weight"))
     {
         envCacheLat(ComplexField, "_" + getName() + "_weight");
+        envTmpLat(ComplexField, "smear");
         create_weight = true;
     }
     envTmpLat(Field, "phift");
@@ -142,8 +144,11 @@ void TStochFreeField<SImpl>::execute(void)
     {
         LOG(Message) << "Caching momentum-space scalar action" << std::endl;
         
+        envGetTmp(ComplexField, smear);
+        SImpl::MomentaSquare(smear);
+        smear = exp(-par().smearing*smear);
         SImpl::MomentumSpacePropagator(w, sqrt(par().m2));
-        w *= par().g/N;
+        w *= par().g/N*smear;
         w  = sqrt(vol)*sqrt(w);
     }
     LOG(Message) << "Generating random momentum-space field" << std::endl;
diff --git a/Hadrons/Modules/MScalarSUN/TwoPointNPR.hpp b/Hadrons/Modules/MScalarSUN/TwoPointNPR.hpp
index 567e4662..83fa7f07 100644
--- a/Hadrons/Modules/MScalarSUN/TwoPointNPR.hpp
+++ b/Hadrons/Modules/MScalarSUN/TwoPointNPR.hpp
@@ -133,13 +133,14 @@ void TTwoPointNPR<SImpl>::setup(void)
 template <typename SImpl>
 void TTwoPointNPR<SImpl>::execute(void)
 {
-    const unsigned int             nd = env().getNd();
-    const unsigned int             nl = env().getDim(0);
+    const unsigned int             nd   = env().getNd();
+    const unsigned int             nl   = env().getDim(0);
+    const double                   invV = 1./env().getVolume();
     FFT                            fft(env().getGrid());
     std::vector<TwoPointNPRResult> result;
-    TwoPointNPRResult              twoPtp1, twoPtp2;
+    TwoPointNPRResult              twoPtp1, twoPtp2, twoPtDisc;
     auto                           &phi    = envGet(Field, par().field);
-    bool                           doTwoPt = true;
+    bool                           doAux = true;
 
     envGetTmp(ComplexField, ftBuf);
     envGetTmp(Field, ftMatBuf);
@@ -151,19 +152,23 @@ void TTwoPointNPR<SImpl>::execute(void)
         std::vector<int>  p1, p2, p;
         Site              phip1, phip2;
         TComplex          opp;
-        TwoPointNPRResult r;
+        TwoPointNPRResult r, rDisc;
 
         LOG(Message) << "FFT: operator '" << opName << "'" << std::endl;
         fft.FFT_all_dim(ftBuf, op, FFT::forward);
         LOG(Message) << "Generating vertex function" << std::endl;
         r.op = opName;
         r.data.resize(nl);
-        if (doTwoPt)
+        rDisc.op = opName + "_disc";
+        rDisc.data.resize(nl);
+        if (doAux)
         {
             twoPtp1.op = "phi_prop_p1";
             twoPtp1.data.resize(nl);
             twoPtp2.op = "phi_prop_p2";
             twoPtp2.data.resize(nl);
+            twoPtDisc.op = "phi_prop_disc";
+            twoPtDisc.data.resize(nl);
         }
         for (unsigned int n = 0; n < nl; ++n)
         {
@@ -184,20 +189,24 @@ void TTwoPointNPR<SImpl>::execute(void)
             peekSite(phip1, ftMatBuf, p1);
             peekSite(phip2, ftMatBuf, p2);
             peekSite(opp, ftBuf, p);
-            if (doTwoPt)
+            if (doAux)
             {
-                twoPtp1.data[n] = TensorRemove(trace(phip1*adj(phip1)));
-                twoPtp2.data[n] = TensorRemove(trace(phip2*adj(phip2)));
+                twoPtp1.data[n]   = invV*TensorRemove(trace(phip1*adj(phip1)));
+                twoPtp2.data[n]   = invV*TensorRemove(trace(phip2*adj(phip2)));
+                twoPtDisc.data[n] = invV*TensorRemove(trace(phip2*adj(phip1)));
             }
-            r.data[n]  = TensorRemove(trace(phip2*adj(phip1))*opp);
+            r.data[n]     = invV*TensorRemove(trace(phip2*adj(phip1))*opp);
+            rDisc.data[n] = invV*TensorRemove(trace(phip1*adj(phip1))*opp);
         }
-        if (doTwoPt)
+        if (doAux)
         {
             result.push_back(twoPtp1);
             result.push_back(twoPtp2);
+            result.push_back(twoPtDisc);
         }
         result.push_back(r);
-        doTwoPt = false;
+        result.push_back(rDisc);
+        doAux = false;
     }
     saveResult(par().output, "twoptnpr", result);
 }