Merged latest changes from develop, in preparation for release.

Grid/algorithms/iterative/ImplicitlyRestartedLanczos.h

@@ -64,7 +64,7 @@ void basisRotate(std::vector<Field> &basis,Eigen::MatrixXd& Qt,int j0, int j1, i
 
   thread_region
   {
-    std::vector < vobj > B(Nm); // Thread private
+    std::vector < vobj , commAllocator<vobj> > B(Nm); // Thread private
     thread_for_in_region(ss, grid->oSites(),{
       for(int j=j0; j<j1; ++j) B[j]=0.;
       

Grid/qcd/smearing/StoutSmearing.h

@@ -40,7 +40,7 @@ template <class Gimpl>
 class Smear_Stout : public Smear<Gimpl> {
  private:
   int OrthogDim = -1;
-  const std::vector<double> SmearRho{};
+  const std::vector<double> SmearRho;
   // Smear<Gimpl>* ownership semantics:
   //    Smear<Gimpl>* passed in to constructor are owned by caller, so we don't delete them here
   //    Smear<Gimpl>* created within constructor need to be deleted as part of the destructor
@@ -61,7 +61,6 @@ public:
 
   /*! Stout smearing with base explicitly specified */
   Smear_Stout(Smear<Gimpl>* base) : SmearBase{base} {
-    std::cout << GridLogDebug << "Stout smearing constructor : Smear_Stout(Smear<Gimpl>* base)" << std::endl
     assert(Nc == 3 && "Stout smearing currently implemented only for Nc==3");
   }
 
@@ -73,32 +72,25 @@ public:
     }
 
   /*! Default constructor. rho is constant in all directions, optionally except for orthogonal dimension */
-  Smear_Stout(double rho, int orthogdim = -1)
+  Smear_Stout(double rho = 1.0, int orthogdim = -1)
-  //: OwnedBase{(orthogdim<0 || orthogdim>=Nd) ? new Smear_APE<Gimpl>(rho) : new Smear_APE<Gimpl>(rho3D(rho,orthogdim))},
   : OrthogDim{orthogdim}, SmearRho{ rho3D(rho,orthogdim) }, OwnedBase{ new Smear_APE<Gimpl>(SmearRho) }, SmearBase{OwnedBase.get()} {
-    std::cout << GridLogDebug << "Stout smearing constructor : Smear_StoutSmear_Stout(double " << rho << ", int " << OrthogDim << " )\nrho3d=" << SmearRho << std::endl;
     assert(Nc == 3 && "Stout smearing currently implemented only for Nc==3");
   }
 
-  /*
-  Smear_Stout(double rho = 1.0) : SmearBase(new Smear_APE<Gimpl>(rho)) {
-    assert(Nc == 3 && "Stout smearing currently implemented only for Nc==3");
-  }*/
-
   ~Smear_Stout() {}  // delete SmearBase...
 
   void smear(GaugeField& u_smr, const GaugeField& U) const {
     GaugeField C(U.Grid());
     GaugeLinkField tmp(U.Grid()), iq_mu(U.Grid()), Umu(U.Grid());
 
-    std::cout << GridLogDebug << "Stout smearing started" << std::endl;
+    std::cout << GridLogDebug << "Stout smearing started\n";
 
     // Smear the configurations
     SmearBase->smear(C, U);
 
     for (int mu = 0; mu < Nd; mu++) {
       if( mu == OrthogDim )
-        tmp = 1.0;
+        tmp = 1.0;  // Don't smear in the orthogonal direction
       else {
         tmp = peekLorentz(C, mu);
         Umu = peekLorentz(U, mu);
@@ -109,7 +101,7 @@ public:
       }
       pokeLorentz(u_smr, tmp * Umu, mu);  // u_smr = exp(iQ_mu)*U_mu
     }
-    std::cout << GridLogDebug << "Stout smearing completed" << std::endl;
+    std::cout << GridLogDebug << "Stout smearing completed\n";
   };
 
   void derivative(GaugeField& SigmaTerm, const GaugeField& iLambda,
@@ -169,7 +161,7 @@ public:
     c0max = 2.0 * pow(tmp, 1.5);
 
     theta = acos(c0 / c0max) *
-            one_over_three;  // divide by three here, now leave as it is
+      one_over_three;  // divide by three here, now leave as it is
     u = sqrt(tmp) * cos(theta);
     w = sqrt(c1) * sin(theta);
   }
@@ -194,7 +186,7 @@ public:
     e2iu = cos(2.0 * u) + timesI(sin(2.0 * u));
 
     h0 = e2iu * (u2 - w2) +
-         emiu * ((8.0 * u2 * cosw) + (2.0 * u * (3.0 * u2 + w2) * ixi0));
+      emiu * ((8.0 * u2 * cosw) + (2.0 * u * (3.0 * u2 + w2) * ixi0));
     h1 = e2iu * (2.0 * u) - emiu * ((2.0 * u * cosw) - (3.0 * u2 - w2) * ixi0);
     h2 = e2iu - emiu * (cosw + (3.0 * u) * ixi0);
 

Hadrons/Modules/MDistil/BContraction.cc

@@ -1,36 +0,0 @@
-/*************************************************************************************
- 
- Grid physics library, www.github.com/paboyle/Grid
- 
- Source file: Hadrons/Modules/MDistil/BContraction.cc
- 
- Copyright (C) 2019
- 
- Author: Felix Erben <ferben@ed.ac.uk>
- Author: Michael Marshall <Michael.Marshall@ed.ac.uk>
- 
- This program is free software; you can redistribute it and/or modify
- it under the terms of the GNU General Public License as published by
- the Free Software Foundation; either version 2 of the License, or
- (at your option) any later version.
- 
- This program is distributed in the hope that it will be useful,
- but WITHOUT ANY WARRANTY; without even the implied warranty of
- MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
- GNU General Public License for more details.
- 
- You should have received a copy of the GNU General Public License along
- with this program; if not, write to the Free Software Foundation, Inc.,
- 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
- 
- See the full license in the file "LICENSE" in the top level distribution directory
- *************************************************************************************/
-/*  END LEGAL */
-
-#include <Hadrons/Modules/MDistil/BContraction.hpp>
-
-using namespace Grid;
-using namespace Hadrons;
-using namespace MDistil;
-
-template class Grid::Hadrons::MDistil::TBContraction<FIMPL>;

Hadrons/Modules/MDistil/BContraction.hpp

@@ -1,289 +0,0 @@
-/*************************************************************************************
- 
- Grid physics library, www.github.com/paboyle/Grid
- 
- Source file: Hadrons/Modules/MDistil/BContraction.hpp
- 
- Copyright (C) 2019
- 
- Author: Felix Erben <ferben@ed.ac.uk>
- Author: Michael Marshall <Michael.Marshall@ed.ac.uk>
- 
- This program is free software; you can redistribute it and/or modify
- it under the terms of the GNU General Public License as published by
- the Free Software Foundation; either version 2 of the License, or
- (at your option) any later version.
- 
- This program is distributed in the hope that it will be useful,
- but WITHOUT ANY WARRANTY; without even the implied warranty of
- MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
- GNU General Public License for more details.
- 
- You should have received a copy of the GNU General Public License along
- with this program; if not, write to the Free Software Foundation, Inc.,
- 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
- 
- See the full license in the file "LICENSE" in the top level distribution directory
- *************************************************************************************/
-/*  END LEGAL */
-
-#ifndef Hadrons_MDistil_BContraction_hpp_
-#define Hadrons_MDistil_BContraction_hpp_
-
-#include <Hadrons/Global.hpp>
-#include <Hadrons/Module.hpp>
-#include <Hadrons/ModuleFactory.hpp>
-#include <Hadrons/Solver.hpp>
-#include <Hadrons/EigenPack.hpp>
-#include <Hadrons/A2AVectors.hpp>
-#include <Hadrons/DilutedNoise.hpp>
-
-// These are members of Distillation
-#include <Hadrons/Distil.hpp>
-BEGIN_HADRONS_NAMESPACE
-
-/******************************************************************************
- *                         BContraction                                 *
- ******************************************************************************/
-BEGIN_MODULE_NAMESPACE(MDistil)
-
-    // general baryon tensor set based on Eigen tensors and Grid-allocated memory
-    // Dimensions:
-    //   0 - ext - external field (momentum, EM field, ...)
-    //   1 - str - dirac structure
-    //   2 - t   - timeslice
-    //   3 - s   - free spin index
-    //   4 - i   - left  distillation mode index
-    //   5 - j   - middle  distillation mode index
-    //   6 - k   - left  distillation mode index
-    // template <typename T>
-    // using BaryonTensorSet = Eigen::TensorMap<Eigen::Tensor<T, 7, Eigen::RowMajor>>;
-
-
-class BContractionPar: Serializable
-{
-public:
-    GRID_SERIALIZABLE_CLASS_MEMBERS(BContractionPar,
-		                    std::string, one,
-		                    std::string, two,
-		                    std::string, three,
-		                    std::string, output,
-				    int, parity,
-                                    std::vector<std::string>, mom);
-};
-
-template <typename FImpl>
-class TBContraction: public Module<BContractionPar>
-{
-public:
-  using W = WilsonImplR; // Debug so I can see type info for default FImpl
-    FERM_TYPE_ALIASES(FImpl,);
-public:
-    // constructor
-    TBContraction(const std::string name);
-    // destructor
-    virtual ~TBContraction(void) {};
-    // dependency relation
-    virtual std::vector<std::string> getInput(void);
-    virtual std::vector<std::string> getOutput(void);
-    // setup
-    virtual void setup(void);
-    // execution
-    virtual void execute(void);
-private:
-    bool                               hasPhase_{false};
-    std::string                        momphName_;
-    std::vector<Gamma::Algebra>        gamma12_;
-    std::vector<Gamma::Algebra>        gamma23_;
-    std::vector<std::vector<Real>>     mom_;
-protected:
-    GridCartesian * grid4d;
-    GridCartesian * grid3d;
-};
-
-/*class BFieldIO: Serializable{
-public:
-  using BaryonTensorSet = Eigen::Tensor<Complex, 7>;
-  GRID_SERIALIZABLE_CLASS_MEMBERS(BFieldIO,
-                                  BaryonTensorSet, BField
-		                  );
-};*/
-
-MODULE_REGISTER_TMP(BContraction, TBContraction<FIMPL>, MDistil);
-
-/******************************************************************************
- *                 TBContraction implementation                             *
- ******************************************************************************/
-// constructor /////////////////////////////////////////////////////////////////
-template <typename FImpl>
-TBContraction<FImpl>::TBContraction(const std::string name)
-: Module<BContractionPar>(name)
-{}
-
-// dependencies/products ///////////////////////////////////////////////////////
-template <typename FImpl>
-std::vector<std::string> TBContraction<FImpl>::getInput(void)
-{
-    std::vector<std::string> in = {par().one, par().two, par().three};
-    
-    return in;
-}
-
-template <typename FImpl>
-std::vector<std::string> TBContraction<FImpl>::getOutput(void)
-{
-    std::vector<std::string> out = {};
-    
-    return out;
-}
-
-// setup ///////////////////////////////////////////////////////////////////////
-template <typename FImpl>
-void TBContraction<FImpl>::setup(void)
-{
-
-}
-
-// execution ///////////////////////////////////////////////////////////////////
-template <typename FImpl>
-void TBContraction<FImpl>::execute(void)
-{
-    auto &one   = envGet(std::vector<FermionField>, par().one);
-    auto &two   = envGet(std::vector<FermionField>, par().two);
-    auto &three = envGet(std::vector<FermionField>, par().three);
-
-    int N_1     = one.size();
-    int N_2     = two.size();
-    int N_3     = three.size();
-
-    int parity = par().parity;
-    const std::string &output{par().output};
-
-    LOG(Message) << "Computing distillation baryon fields" << std::endl;
-    LOG(Message) << "One: '" << par().one << "' Two: '" << par().two  << "' Three: '" << par().three << "'" << std::endl;
-    LOG(Message) << "Momenta:" << std::endl;
-    for (auto &p: mom_)
-    {
-        LOG(Message) << "  " << p << std::endl;
-    }
-
-    grid4d = env().getGrid();
-    grid3d = MakeLowerDimGrid(grid4d);    
-    int Nmom=1;
-    int Nt=64;
-   // std::vector<Complex> BField(Nmom*Nt*N_1*N_2*N_3);    
-    int Bindex;
-    int Nc=3; //Num colours
-
-    FermionField ftmp1(grid3d);
-    FermionField ftmp2(grid3d);
-    FermionField ftmp3(grid3d);
-    LatticeView<typename FImpl::SiteSpinor> tmp1{ ftmp1 };
-    LatticeView<typename FImpl::SiteSpinor> tmp2{ ftmp2 };
-    LatticeView<typename FImpl::SiteSpinor> tmp3{ ftmp3 };
-    //std::complex<double> * tmp33 = reinterpret_cast<std::complex<double> *>(&(tmp3[0]()(0)(0)));
-
-#ifdef THIS_IS_NAUGHTY_TODO_FIXME
-  // The reinterpret_cast gets rid of SIMD attributes
-  // Plus other badness - e.g. we perhaps shouldn't explicitly say SpinColourVector
-  // ... but rather use some correct FIMPL types
-  // REVIEW WITH PETER
-#endif
-  
-    SpinColourVector * tmp11 = reinterpret_cast<SpinColourVector *>(&(tmp1[0]()(0)(0)));
-    SpinColourVector * tmp22 = reinterpret_cast<SpinColourVector *>(&(tmp2[0]()(0)(0)));
-    SpinColourVector * tmp33 = reinterpret_cast<SpinColourVector *>(&(tmp3[0]()(0)(0)));
-
-    SpinVector tmp11s; 
-    SpinVector tmp22s; 
-    SpinVector tmp33s; 
-    SpinVector tmp333;
-    SpinMatrix diquark;
-    SpinMatrix g_diquark;
-    SpinVector tmp222;
-    SpinVector tmp111;
-
-
-    assert(parity == 1 || parity == -1);
-
-    std::vector<std::vector<int>> epsilon = {{0,1,2},{1,2,0},{2,0,1},{0,2,1},{2,1,0},{1,0,2}};
-    std::vector<int> epsilon_sgn = {1,1,1,-1,-1,-1};
-
-    Gamma g4(Gamma::Algebra::GammaT);
-
-    gamma12_ = {
-       Gamma::Algebra::Identity, // I
-       Gamma::Algebra::Gamma5, // gamma_5
-       Gamma::Algebra::Identity, // I
-    };
-    gamma23_ = { // C = i gamma_2 gamma_4
-       Gamma::Algebra::SigmaXZ, // C gamma_5 = -i gamma_1 gamma_3
-       Gamma::Algebra::SigmaYT, // C = i gamma_2 gamma_4
-       Gamma::Algebra::GammaYGamma5, // i gamma_4 C gamma_5 = i gamma_2 gamma_5
-    };
-  std::vector<Complex> factor23{{0.,-1.},{0.,1.},{0.,1.}}; 
-  using BaryonTensorSet = Eigen::Tensor<Complex, 7>;
-  int Ngamma=3;
-  BaryonTensorSet BField3(Nmom,Ngamma,Nt,4,N_1,N_2,N_3); 
-
-  Eigen::Tensor<Complex, 3> corr(Nmom,4,Nt); 
-
-   
-    Complex diquark2;
-    for (int i1=0 ; i1 < N_1 ; i1++){
-      for (int i2=0 ; i2 < N_2 ; i2++){
-        for (int i3=0 ; i3 < N_3 ; i3++){
-          for (int imom=0 ; imom < Nmom ; imom++){
-            for (int t=0 ; t < Nt ; t++){
-              Bindex = i1 + N_1*(i2 + N_2*(i3 + N_3*(imom+Nmom*t)));
-	      ExtractSliceLocal(ftmp1,one[i1],0,t,3);
-	      ExtractSliceLocal(ftmp2,two[i2],0,t,3);
-	      ExtractSliceLocal(ftmp3,three[i3],0,t,3);
-              accelerator_for(sU, grid3d->oSites(), grid3d->Nsimd(),
-              {
-                for (int ie=0 ; ie < 6 ; ie++){
-		  // Why does peekColour not work????
-                  for (int is=0 ; is < 4 ; is++){
-	            tmp11s()(is)() = tmp11[sU]()(is)(epsilon[ie][0]);
-	            tmp22s()(is)() = tmp22[sU]()(is)(epsilon[ie][1]);
-	            tmp33s()(is)() = tmp33[sU]()(is)(epsilon[ie][2]);
-		  }
-                  for (int ig=0 ; ig < Ngamma ; ig++){
-		    tmp333 = Gamma(gamma23_[ig])*tmp33s;
-		    tmp111 = Gamma(gamma12_[ig])*tmp11s;
-		    tmp222 = g4*tmp111;
-		    tmp111 = 0.5*(double)parity*(tmp111 + tmp222); // P_\pm * ...
-                    diquark2 = factor23[0]*innerProduct(tmp22s,tmp333);
-                    for (int is=0 ; is < 4 ; is++){
-                      BField3(imom,ig,t,is,i1,i2,i3)+=static_cast<Real>(epsilon_sgn[ie])*tmp111()(is)()*diquark2;
-                    }
-  		  }
-		}
-              } );
-            }
-	  }
-	}
-      }
-    }
-    for (int is=0 ; is < 4 ; is++){
-      for (int t=0 ; t < Nt ; t++){
-        std::cout << "BaryonField(is=" << is << ",t=" << t << ") = " << BField3(0,0,t,is,0,0,0) << std::endl;
-      }
-    }
-
-    BFieldIO BField_save;
-#ifdef FELIX_ISSUE
-    BField_save.BField = BField3;
-#endif
-  std::string filename ="./" + output + ".h5"; 
-  std::cout << "Writing to file " << filename << std::endl;
-  Hdf5Writer writer(filename);
-  write(writer,"BaryonField",BField_save.BField);
-
-}
-
-END_MODULE_NAMESPACE
-
-END_HADRONS_NAMESPACE
-
-#endif // Hadrons_MDistil_BContraction_hpp_

Hadrons/Modules/MDistil/Baryon2pt.cc

@@ -1,36 +0,0 @@
-/*************************************************************************************
- 
- Grid physics library, www.github.com/paboyle/Grid
- 
- Source file: Hadrons/Modules/MDistil/Baryon2pt.cc
- 
- Copyright (C) 2019
- 
- Author: Felix Erben <ferben@ed.ac.uk>
- Author: Michael Marshall <Michael.Marshall@ed.ac.uk>
- 
- This program is free software; you can redistribute it and/or modify
- it under the terms of the GNU General Public License as published by
- the Free Software Foundation; either version 2 of the License, or
- (at your option) any later version.
- 
- This program is distributed in the hope that it will be useful,
- but WITHOUT ANY WARRANTY; without even the implied warranty of
- MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
- GNU General Public License for more details.
- 
- You should have received a copy of the GNU General Public License along
- with this program; if not, write to the Free Software Foundation, Inc.,
- 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
- 
- See the full license in the file "LICENSE" in the top level distribution directory
- *************************************************************************************/
-/*  END LEGAL */
-
-#include <Hadrons/Modules/MDistil/Baryon2pt.hpp>
-
-using namespace Grid;
-using namespace Hadrons;
-using namespace MDistil;
-
-template class Grid::Hadrons::MDistil::TBaryon2pt<FIMPL>;

Hadrons/Modules/MDistil/Baryon2pt.hpp

@@ -1,221 +0,0 @@
-/*************************************************************************************
- 
- Grid physics library, www.github.com/paboyle/Grid
- 
- Source file: Hadrons/Modules/MDistil/Baryon2pt.hpp
- 
- Copyright (C) 2019
- 
- Author: Felix Erben <ferben@ed.ac.uk>
- Author: Michael Marshall <Michael.Marshall@ed.ac.uk>
- 
- This program is free software; you can redistribute it and/or modify
- it under the terms of the GNU General Public License as published by
- the Free Software Foundation; either version 2 of the License, or
- (at your option) any later version.
- 
- This program is distributed in the hope that it will be useful,
- but WITHOUT ANY WARRANTY; without even the implied warranty of
- MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
- GNU General Public License for more details.
- 
- You should have received a copy of the GNU General Public License along
- with this program; if not, write to the Free Software Foundation, Inc.,
- 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
- 
- See the full license in the file "LICENSE" in the top level distribution directory
- *************************************************************************************/
-/*  END LEGAL */
-
-#ifndef Hadrons_MDistil_Baryon2pt_hpp_
-#define Hadrons_MDistil_Baryon2pt_hpp_
-
-#include <Hadrons/Global.hpp>
-#include <Hadrons/Module.hpp>
-#include <Hadrons/ModuleFactory.hpp>
-#include <Hadrons/Solver.hpp>
-#include <Hadrons/EigenPack.hpp>
-#include <Hadrons/A2AVectors.hpp>
-#include <Hadrons/DilutedNoise.hpp>
-
-// These are members of Distillation
-#include <Hadrons/Distil.hpp>
-BEGIN_HADRONS_NAMESPACE
-
-/******************************************************************************
- *                         Baryon2pt                                 *
- ******************************************************************************/
-BEGIN_MODULE_NAMESPACE(MDistil)
-
-class Baryon2ptPar: Serializable
-{
-public:
-    GRID_SERIALIZABLE_CLASS_MEMBERS(Baryon2ptPar,
-		                    std::string, inputL,
-		                    std::string, inputR,
-		                    std::string, quarksL,
-		                    std::string, quarksR,
-		                    std::string, output
-                                    );
-};
-
-template <typename FImpl>
-class TBaryon2pt: public Module<Baryon2ptPar>
-{
-public:
-    // constructor
-    TBaryon2pt(const std::string name);
-    // destructor
-    virtual ~TBaryon2pt(void) {};
-    // dependency relation
-    virtual std::vector<std::string> getInput(void);
-    virtual std::vector<std::string> getOutput(void);
-    // setup
-    virtual void setup(void);
-    // execution
-    virtual void execute(void);
-};
-
-class C2IO: Serializable{
-public:
-using C2Set = Eigen::Tensor<Complex, 2>;
-  GRID_SERIALIZABLE_CLASS_MEMBERS(C2IO,
-                                  C2Set, C2
-		                  );
-};
-
-MODULE_REGISTER_TMP(Baryon2pt, TBaryon2pt<FIMPL>, MDistil);
-
-/******************************************************************************
- *                 TBaryon2pt implementation                             *
- ******************************************************************************/
-// constructor /////////////////////////////////////////////////////////////////
-template <typename FImpl>
-TBaryon2pt<FImpl>::TBaryon2pt(const std::string name)
-: Module<Baryon2ptPar>(name)
-{}
-
-// dependencies/products ///////////////////////////////////////////////////////
-template <typename FImpl>
-std::vector<std::string> TBaryon2pt<FImpl>::getInput(void)
-{
-    std::vector<std::string> in;
-    
-    return in;
-}
-
-template <typename FImpl>
-std::vector<std::string> TBaryon2pt<FImpl>::getOutput(void)
-{
-    std::vector<std::string> out = {getName()};
-    
-    return out;
-}
-
-// setup ///////////////////////////////////////////////////////////////////////
-template <typename FImpl>
-void TBaryon2pt<FImpl>::setup(void)
-{
-    
-}
-
-// execution ///////////////////////////////////////////////////////////////////
-template <typename FImpl>
-void TBaryon2pt<FImpl>::execute(void)
-{
-
-    const std::string &inputL{par().inputL};
-    const std::string &inputR{par().inputR};
-    const std::string &quarksL{par().quarksL};
-    const std::string &quarksR{par().quarksR};
-    const std::string &output{par().output};
-
-    int Nmom=1;
-    int Nt=32;
-    int Nc=3; //Num colours
-    std::vector<std::vector<int>> epsilon = {{0,1,2},{1,2,0},{2,0,1},{0,2,1},{2,1,0},{1,0,2}};
-    std::vector<int> epsilon_sgn = {1,1,1,-1,-1,-1};
-    int Ngamma=3;
-
-    int N_1=12;
-    int N_2=12;
-    int N_3=12;
-
-  //  using BaryonTensorSet = Eigen::Tensor<Complex, 7>;
-
-    BFieldIO BFieldL;
-    BFieldL.BField.resize(Nmom,Nt,N_1,N_2,N_3,4);
-
-    std::string filenameL ="./" + inputL + ".h5"; 
-    std::cout << "Reading from file " << filenameL << std::endl;
-    Hdf5Reader readerL(filenameL);
-    read(readerL,"BaryonField",BFieldL.BField);
-
-    BFieldIO BFieldR;
-    BFieldR.BField.resize(Nmom,Nt,N_1,N_2,N_3,4);
-
-    std::string filenameR ="./" + inputR + ".h5"; 
-    std::cout << "Reading from file " << filenameR << std::endl;
-    Hdf5Reader readerR(filenameR);
-    read(readerR,"BaryonField",BFieldR.BField);
-
-    Eigen::Tensor<Complex, 2> corr(Nmom,Nt); 
-
-    int Npairs = 0;
-    char left[] = "uud";
-    char right[] = "uud";
-    std::vector<int> pairs(6);
-    for (int ie=0, i=0 ; ie < 6 ; ie++){
-      if (left[0] == right[epsilon[ie][0]] && left[1] == right[epsilon[ie][1]] && left[2] == right[epsilon[ie][2]]){
-        pairs[i] = ie;
-        i++;
-        Npairs++;
-      }
-    }
-    pairs.resize(Npairs);
-    std::cout << Npairs << " pairs: " << pairs << std::endl;
-    for (int imom=0 ; imom < Nmom ; imom++){
-      for (int t=0 ; t < Nt ; t++){
-        corr(imom,t) = 0.;
-      }
-    }
-
-    int tsrc=0;
-
-    for (int ipair=0 ; ipair < Npairs ; ipair++){
-      Eigen::array<Eigen::IndexPair<int>, 3> product_dims = { Eigen::IndexPair<int>(0,epsilon[pairs[ipair]][0]),Eigen::IndexPair<int>(1,epsilon[pairs[ipair]][1]) ,Eigen::IndexPair<int>(2,epsilon[pairs[ipair]][2])  };
-      for (int imom=0 ; imom < Nmom ; imom++){
-        std::cout << imom << std::endl;
-        Eigen::Tensor<Complex,5> B5L = BFieldL.BField.chip(imom,0);
-        Eigen::Tensor<Complex,5> B5R = BFieldR.BField.chip(imom,0);
-        for (int t=0 ; t < Nt ; t++){
-          Eigen::Tensor<Complex,4> B4L = B5L.chip(t,0);
-          Eigen::Tensor<Complex,4> B4R = B5R.chip(tsrc,0);
-          for (int is=0 ; is < 4 ; is++){
-            Eigen::Tensor<Complex,3> B3L = B4L.chip(is,3);
-            Eigen::Tensor<Complex,3> B3R = B4R.chip(is,3);
-            Eigen::Tensor<Complex,0> C2 = B3L.contract(B3R,product_dims);
-            corr(imom,t) += static_cast<Real>(epsilon_sgn[pairs[ipair]])*C2(0);
-          }
-        }
-      }
-    }
-    for (int t=0 ; t < Nt ; t++){
-      std::cout << "C2(t=" << t << ") = " << corr(0,t) << std::endl;
-    }    
-
-  /*  C2IO C2_save;
-    C2_save.C2 = corr;
-
-    std::string filename ="./" + output + ".h5"; 
-    std::cout << "Writing to file " << filename << std::endl;
-    Hdf5Writer writer(filename);
-    write(writer,"C2",C2_save.C2);
-  */
-}
-
-END_MODULE_NAMESPACE
-
-END_HADRONS_NAMESPACE
-
-#endif // Hadrons_MDistil_Baryon2pt_hpp_