diff --git a/extras/Hadrons/Modules/MScalar/ScalarVP.cc b/extras/Hadrons/Modules/MScalar/ScalarVP.cc
index 7745c40e..880b3ffd 100644
--- a/extras/Hadrons/Modules/MScalar/ScalarVP.cc
+++ b/extras/Hadrons/Modules/MScalar/ScalarVP.cc
@@ -205,15 +205,15 @@ void TScalarVP::execute(void)
     }
 
     // CONTRACTIONS
-    vpTensor_.clear();
+    std::vector<std::vector<ScalarField> > vpTensor;
     for (unsigned int mu = 0; mu < env().getNd(); ++mu)
     {
-        std::vector<ScalarField *> vpTensor_mu;
+        std::vector<ScalarField> vpTensor_mu;
         for (unsigned int nu = 0; nu < env().getNd(); ++nu)
         {
-            vpTensor_mu.push_back(env().createLattice<ScalarField>(vpTensorName_[mu][nu]));
+            vpTensor_mu.push_back(*env().createLattice<ScalarField>(vpTensorName_[mu][nu]));
         }
-        vpTensor_.push_back(vpTensor_mu);
+        vpTensor.push_back(vpTensor_mu);
     }
     ScalarField prop1(env().getGrid()), prop2(env().getGrid());
     EmField     &A = *env().getObject<EmField>(par().emField);
@@ -221,7 +221,6 @@ void TScalarVP::execute(void)
     TComplex    Anu0;
     std::vector<int> coor0 = {0, 0, 0, 0};
 
-    // Position-space implementation
     prop1 = *GFSrc_ + q*propQ + q*q*propSun + q*q*propTad;
     fft.FFT_all_dim(prop1, prop1, FFT::backward);
     for (unsigned int nu = 0; nu < env().getNd(); ++nu)
@@ -231,57 +230,24 @@ void TScalarVP::execute(void)
                 + q*q*(*(muPropSun_[nu]) + *(muPropTad_[nu]));
         fft.FFT_all_dim(prop2, prop2, FFT::backward);
 
+        std::vector<ScalarField> pi_nu;
         for (unsigned int mu = 0; mu < env().getNd(); ++mu)
         {
             LOG(Message) << "Computing Pi[" << mu << "][" << nu << "]..."
                          << std::endl;
             Amu = peekLorentz(A, mu);
-            ScalarField &pi_mu_nu = *(vpTensor_[mu][nu]);
-            pi_mu_nu = adj(prop2)
-                       * (1.0 + ci*q*Amu - 0.5*q*q*Amu*Amu)
-                       * Cshift(prop1, mu, 1)
-                       * (1.0 + ci*q*Anu0 - 0.5*q*q*Anu0*Anu0);
-            pi_mu_nu -= Cshift(adj(prop2), mu, 1)
-                        * (1.0 - ci*q*Amu - 0.5*q*q*Amu*Amu)
-                        * prop1
-                        * (1.0 + ci*q*Anu0 - 0.5*q*q*Anu0*Anu0);
-            pi_mu_nu = 2.0*real(pi_mu_nu);
+            vpTensor[mu][nu] = adj(prop2)
+                               * (1.0 + ci*q*Amu - 0.5*q*q*Amu*Amu)
+                               * Cshift(prop1, mu, 1)
+                               * (1.0 + ci*q*Anu0 - 0.5*q*q*Anu0*Anu0);
+            vpTensor[mu][nu] -= Cshift(adj(prop2), mu, 1)
+                                * (1.0 - ci*q*Amu - 0.5*q*q*Amu*Amu)
+                                * prop1
+                                * (1.0 + ci*q*Anu0 - 0.5*q*q*Anu0*Anu0);
+            vpTensor[mu][nu] = 2.0*real(vpTensor[mu][nu]);
         }
     }
 
-    // // Momentum-space implementation
-    // ScalarField propbuf1(env().getGrid()), propbuf2(env().getGrid());
-    // prop1 = *GFSrc_ + q*propQ + q*q*propSun + q*q*propTad;
-    // for (unsigned int nu = 0; nu < env().getNd(); ++nu)
-    // {
-    //     peekSite(Anu0, peekLorentz(A, nu), coor0);
-    //     prop2 = adj(*phase_[nu])*(*GFSrc_) + q*(*(muPropQ_[nu]))
-    //             + q*q*(*(muPropSun_[nu]) + *(muPropTad_[nu]));
-
-    //     for (unsigned int mu = 0; mu < env().getNd(); ++mu)
-    //     {
-    //         LOG(Message) << "Computing Pi[" << mu << "][" << nu << "]..."
-    //                      << std::endl;
-    //         Amu = peekLorentz(A, mu);
-    //         ScalarField &pi_mu_nu = *(vpTensor_[mu][nu]);
-    //         propbuf1 = (*phase_[mu])*prop1;
-    //         fft.FFT_all_dim(propbuf1, propbuf1, FFT::backward);
-    //         fft.FFT_all_dim(propbuf2, prop2, FFT::backward);
-    //         pi_mu_nu = adj(propbuf2)
-    //                    * (1.0 + ci*q*Amu - 0.5*q*q*Amu*Amu)
-    //                    * propbuf1
-    //                    * (1.0 + ci*q*Anu0 - 0.5*q*q*Anu0*Anu0);
-    //         propbuf2 = (*phase_[mu])*prop2;
-    //         fft.FFT_all_dim(propbuf1, prop1, FFT::backward);
-    //         fft.FFT_all_dim(propbuf2, propbuf2, FFT::backward);
-    //         pi_mu_nu -= adj(propbuf2)
-    //                     * (1.0 - ci*q*Amu - 0.5*q*q*Amu*Amu)
-    //                     * propbuf1
-    //                     * (1.0 + ci*q*Anu0 - 0.5*q*q*Anu0*Anu0);
-    //         pi_mu_nu = 2.0*real(pi_mu_nu);
-    //     }
-    // }
-
     // OUTPUT IF NECESSARY
     if (!par().output.empty())
     {
@@ -302,7 +268,7 @@ void TScalarVP::execute(void)
         {
             for (unsigned int nu = 0; nu < env().getNd(); ++nu)
             {
-                sliceSum(*(vpTensor_[mu][nu]), vecBuf, Tp);
+                sliceSum(vpTensor[mu][nu], vecBuf, Tp);
                 result.resize(vecBuf.size());
                 for (unsigned int t = 0; t < vecBuf.size(); ++t)
                 {
diff --git a/extras/Hadrons/Modules/MScalar/ScalarVP.hpp b/extras/Hadrons/Modules/MScalar/ScalarVP.hpp
index fbe73d85..9d884575 100644
--- a/extras/Hadrons/Modules/MScalar/ScalarVP.hpp
+++ b/extras/Hadrons/Modules/MScalar/ScalarVP.hpp
@@ -57,7 +57,6 @@ private:
     ScalarField                                 *freeMomProp_, *GFSrc_,
                                                 *prop0_;
     std::vector<ScalarField *>                  phase_;
-    std::vector<std::vector<ScalarField *> >    vpTensor_;
     EmField                                     *A;
 };