diff --git a/extras/Hadrons/Modules/MScalar/ChargedProp.cc b/extras/Hadrons/Modules/MScalar/ChargedProp.cc
index 40d4504c..b76ea8d2 100644
--- a/extras/Hadrons/Modules/MScalar/ChargedProp.cc
+++ b/extras/Hadrons/Modules/MScalar/ChargedProp.cc
@@ -23,7 +23,8 @@ std::vector<std::string> TChargedProp::getInput(void)
 
 std::vector<std::string> TChargedProp::getOutput(void)
 {
-    std::vector<std::string> out = {getName()};
+    std::vector<std::string> out = {getName(), getName()+"_0", getName()+"_D1",
+                                    getName()+"_D1D1", getName()+"_D2"};
     
     return out;
 }
@@ -38,6 +39,10 @@ void TChargedProp::setup(void)
         phaseName_.push_back("_shiftphase_" + std::to_string(mu));
     }
     GFSrcName_ = "_" + getName() + "_DinvSrc";
+    prop0Name_ = getName() + "_0";
+    propD1Name_ = getName() + "_D1";
+    propD1D1Name_ = getName() + "_D1D1";
+    propD2Name_ = getName() + "_D2";
     if (!env().hasRegisteredObject(freeMomPropName_))
     {
         env().registerLattice<ScalarField>(freeMomPropName_);
@@ -53,7 +58,14 @@ void TChargedProp::setup(void)
     {
         env().registerLattice<ScalarField>(GFSrcName_);
     }
+    if (!env().hasRegisteredObject(prop0Name_))
+    {
+        env().registerLattice<ScalarField>(prop0Name_);
+    }
     env().registerLattice<ScalarField>(getName());
+    env().registerLattice<ScalarField>(propD1Name_);
+    env().registerLattice<ScalarField>(propD1D1Name_);
+    env().registerLattice<ScalarField>(propD2Name_);
 }
 
 // execution ///////////////////////////////////////////////////////////////////
@@ -64,7 +76,7 @@ void TChargedProp::execute(void)
     Complex     ci(0.0,1.0);
     FFT         fft(env().getGrid());
     
-    // cache free scalar propagator
+    // cache momentum-space free scalar propagator
     if (!env().hasCreatedObject(freeMomPropName_))
     {
         LOG(Message) << "Caching momentum space free scalar propagator"
@@ -88,6 +100,17 @@ void TChargedProp::execute(void)
     {
         GFSrc_ = env().getObject<ScalarField>(GFSrcName_);
     }
+    // cache free scalar propagator
+    if (!env().hasCreatedObject(prop0Name_))
+    {
+        prop0_ = env().createLattice<ScalarField>(prop0Name_);
+        *prop0_ = *GFSrc_;
+        fft.FFT_all_dim(*prop0_, *prop0_, FFT::backward);
+    }
+    else
+    {
+        prop0_ = env().getObject<ScalarField>(prop0Name_);
+    }
     // cache phases
     if (!env().hasCreatedObject(phaseName_[0]))
     {
@@ -117,30 +140,33 @@ void TChargedProp::execute(void)
                  << ", charge= " << par().charge << ")..." << std::endl;
     
     ScalarField &prop   = *env().createLattice<ScalarField>(getName());
+    ScalarField &propD1   = *env().createLattice<ScalarField>(propD1Name_);
+    ScalarField &propD1D1   = *env().createLattice<ScalarField>(propD1D1Name_);
+    ScalarField &propD2   = *env().createLattice<ScalarField>(propD2Name_);
     ScalarField buf(env().getGrid());
     ScalarField &GFSrc = *GFSrc_, &G = *freeMomProp_;
     double      q = par().charge;
     
-    // G*F*Src
-    prop = GFSrc;
-
-    // - q*G*momD1*G*F*Src (momD1 = F*D1*Finv)
+    // -G*momD1*G*F*Src (momD1 = F*D1*Finv)
     buf = GFSrc;
     momD1(buf, fft);
-    buf = G*buf;
-    prop = prop - q*buf;
+    buf = -G*buf;
+    fft.FFT_all_dim(propD1, buf, FFT::backward);
 
-    // + q^2*G*momD1*G*momD1*G*F*Src (here buf = G*momD1*G*F*Src)
+    // G*momD1*G*momD1*G*F*Src (here buf = G*momD1*G*F*Src)
+    buf = -buf;
     momD1(buf, fft);
-    prop = prop + q*q*G*buf;
+    propD1D1 = G*buf;
+    fft.FFT_all_dim(propD1D1, propD1D1, FFT::backward);
 
-    // - q^2*G*momD2*G*F*Src (momD2 = F*D2*Finv)
+    // -G*momD2*G*F*Src (momD2 = F*D2*Finv)
     buf = GFSrc;
     momD2(buf, fft);
-    prop = prop - q*q*G*buf;
+    buf = -G*buf;
+    fft.FFT_all_dim(propD2, buf, FFT::backward);
 
-    // final FT
-    fft.FFT_all_dim(prop, prop, FFT::backward);
+    // full charged scalar propagator
+    prop = (*prop0_) + q*propD1 + q*q*propD1D1 + q*q*propD2;
     
     // OUTPUT IF NECESSARY
     if (!par().output.empty())
@@ -155,14 +181,48 @@ void TChargedProp::execute(void)
         std::vector<TComplex> vecBuf;
         std::vector<Complex>  result;
         
+        write(writer, "charge", q);
+
+        // Write full propagator
         sliceSum(prop, vecBuf, Tp);
         result.resize(vecBuf.size());
         for (unsigned int t = 0; t < vecBuf.size(); ++t)
         {
             result[t] = TensorRemove(vecBuf[t]);
         }
-        write(writer, "charge", q);
         write(writer, "prop", result);
+
+        // Write free propagator
+        sliceSum(*prop0_, vecBuf, Tp);
+        for (unsigned int t = 0; t < vecBuf.size(); ++t)
+        {
+            result[t] = TensorRemove(vecBuf[t]);
+        }
+        write(writer, "prop_0", result);
+
+        // Write propagator D1 term
+        sliceSum(propD1, vecBuf, Tp);
+        for (unsigned int t = 0; t < vecBuf.size(); ++t)
+        {
+            result[t] = TensorRemove(vecBuf[t]);
+        }
+        write(writer, "prop_D1", result);
+
+        // Write propagator D1D1 term
+        sliceSum(propD1D1, vecBuf, Tp);
+        for (unsigned int t = 0; t < vecBuf.size(); ++t)
+        {
+            result[t] = TensorRemove(vecBuf[t]);
+        }
+        write(writer, "prop_D1D1", result);
+
+        // Write propagator D2 term
+        sliceSum(propD2, vecBuf, Tp);
+        for (unsigned int t = 0; t < vecBuf.size(); ++t)
+        {
+            result[t] = TensorRemove(vecBuf[t]);
+        }
+        write(writer, "prop_D2", result);
     }
 }
 
diff --git a/extras/Hadrons/Modules/MScalar/ChargedProp.hpp b/extras/Hadrons/Modules/MScalar/ChargedProp.hpp
index 8bb5faa0..6a6c6c39 100644
--- a/extras/Hadrons/Modules/MScalar/ChargedProp.hpp
+++ b/extras/Hadrons/Modules/MScalar/ChargedProp.hpp
@@ -45,9 +45,10 @@ private:
     void momD1(ScalarField &s, FFT &fft);
     void momD2(ScalarField &s, FFT &fft);
 private:
-    std::string                freeMomPropName_, GFSrcName_;
+    std::string                freeMomPropName_, GFSrcName_, prop0Name_,
+                               propD1Name_, propD1D1Name_, propD2Name_;
     std::vector<std::string>   phaseName_;
-    ScalarField                *freeMomProp_, *GFSrc_;
+    ScalarField                *freeMomProp_, *GFSrc_, *prop0_;
     std::vector<ScalarField *> phase_;
     EmField                    *A;
 };