inline_barspec_db_w.cc

Go to the documentation of this file.

/*! \file
 * \brief Inline construction of hadron spectrum
 *
 * Spectrum calculations
 */
 
 
#include "util/ferm/key_val_db.h"
#include "util/ferm/key_hadron_2pt_corr.h"
 
#include "inline_barspec_db_w.h"
#include "meas/inline/abs_inline_measurement_factory.h"
#include "meas/glue/mesplq.h"
#include "util/ft/sftmom.h"
#include "util/info/proginfo.h"
#include "io/param_io.h"
#include "io/qprop_io.h"
#include "meas/hadron/mesons2_w.h"
#include "meas/hadron/barhqlq_w.h"
#include "meas/hadron/curcor2_w.h"
#include "meas/inline/make_xml_file.h"
#include "meas/inline/io/named_objmap.h"
#include "meas/smear/no_quark_displacement.h"
 
#include <sstream>
 
 
namespace Chroma 
{ 
  namespace InlineBarSpecEnv 
  { 
    namespace
    {
      AbsInlineMeasurement* createMeasurement(XMLReader& xml_in, 
                                              const std::string& path) 
      {
        return new InlineMeas(Params(xml_in, path));
      }
 
      //! Local registration flag
      bool registered = false;
 
      const std::string name = "BARYON_DB_SPECTRUM";
    }
 
    //! Register all the factories
    bool registerAll() 
    {
      bool success = true; 
      if (! registered)
      {
        success &= TheInlineMeasurementFactory::Instance().registerObject(name, createMeasurement);
        registered = true;
      }
      return success;
    }
 
 
 
    void read(XMLReader& xml, const std::string& path, Params::SpinTerms_t& ter){
      XMLReader paramtop(xml, path);    
      read(paramtop, "spin"   ,ter.spin   );
      read(paramtop, "weight" ,ter.weight );
    }
 
    void write(XMLWriter& xml, const std::string& path, const Params::SpinTerms_t& ter){
      push(xml, path);    
      write(xml, "spin"   ,ter.spin   );
      write(xml, "weight" ,ter.weight );
      pop(xml);
    }
 
    void read(XMLReader& xml, const std::string& path, Params::SpinWF_t& spWF){
      XMLReader paramtop(xml, path);
      read(paramtop, "norm"     ,spWF.norm   );
      read(paramtop, "terms"    ,spWF.terms  );
    }
 
 
    void write(XMLWriter& xml,const std::string& path, const Params::SpinWF_t& spWF){
      push(xml, path);
      write(xml, "norm"     ,spWF.norm   );
      write(xml, "terms"    ,spWF.terms  );
      pop(xml);
    }
 
  
    void read(XMLReader& xml, const std::string& path, Params::Operators_t& op){
      XMLReader paramtop(xml, path);
      //QDPIO::cout<<"Reading Operators: "<<path<<std::endl ;
      //xml.print(std::cout);
      read(paramtop, "name"     ,op.name   );
      read(paramtop, "spinWF"   ,op.spinWF );
    }
 
    void write(XMLWriter& xml, const std::string& path, const Params::Operators_t& op){
      push(xml, path);
      write(xml, "name"     ,op.name   );
      write(xml, "spinWF"   ,op.spinWF );
      pop(xml);
    }
 
    void read(XMLReader& xml,const std::string& path, Params::State_t& s){
      XMLReader paramtop(xml, path);
      //QDPIO::cout<<"Reading Operators: "<<path<<std::endl ;
      //xml.print(std::cout);
      read(paramtop, "name"      ,s.name   );
      read(paramtop, "flavor"    ,s.flavor );
      read(paramtop, "db"        ,s.db     );
      read(paramtop, "spin"      ,s.spin   );
      read(paramtop, "Operators" ,s.ops    );
    }
 
    void write(XMLWriter& xml, const std::string& path, const Params::State_t& s){
      push(xml, path);
      write(xml, "name"      ,s.name   );
      write(xml, "flavor"    ,s.flavor );
      write(xml, "db"        ,s.db     );
      write(xml, "spin"      ,s.spin   );
      write(xml, "Operators" ,s.ops    );
      pop(xml);
    }
 
 
 
    //! Reader for parameters
    void read(XMLReader& xml, const std::string& path, Params::Param_t& param)
    {
      XMLReader paramtop(xml, path);
 
      int version;
      read(paramtop, "version", version);
 
      switch (version) 
      {
      case 1:
        break;
 
      default:
        QDPIO::cerr << "Input parameter version " << version << " unsupported." << std::endl;
        QDP_abort(1);
      }
 
      //read(paramtop, "MesonP", param.MesonP);
      //read(paramtop, "CurrentP", param.CurrentP);
      //read(paramtop, "BaryonP", param.BaryonP);
 
      read(paramtop, "time_rev", param.time_rev);
      read(paramtop, "mom2_max", param.mom2_max);
      read(paramtop, "avg_equiv_mom", param.avg_equiv_mom);
      read(paramtop, "ensemble", param.ensemble);
      read(paramtop, "States", param.states);
    }
 
 
    //! Writer for parameters
    void write(XMLWriter& xml, const std::string& path, const Params::Param_t& param)
    {
      push(xml, path);
 
      int version = 1;
      write(xml, "version", version);
 
      //write(xml, "MesonP", param.MesonP);
      //write(xml, "CurrentP", param.CurrentP);
      //write(xml, "BaryonP", param.BaryonP);
 
      write(xml, "time_rev", param.time_rev);
 
      write(xml, "mom2_max", param.mom2_max);
      write(xml, "avg_equiv_mom", param.avg_equiv_mom);
      write(xml, "ensemble", param.ensemble);
      write(xml, "States", param.states);
 
      pop(xml);
    }
 
  
 
    //! Propagator input
    void read(XMLReader& xml, const std::string& path, Params::NamedObject_t::Props_t& input)
    {
      XMLReader inputtop(xml, path);
 
      read(inputtop, "up_id", input.up_id);
      if(inputtop.count("down_id")!=0)
        read(inputtop, "down_id", input.down_id);
      else
        input.down_id = "NULL";
      if(inputtop.count("strange_id")!=0)
        read(inputtop, "strange_id", input.strange_id);
      else
        input.strange_id = "NULL";
      if(inputtop.count("charm_id")!=0)
        read(inputtop, "charm_id", input.charm_id);
      else
        input.charm_id = "NULL";
    }
 
    //! Propagator output
    void write(XMLWriter& xml, const std::string& path, const Params::NamedObject_t::Props_t& input)
    {
      push(xml, path);
 
      write(xml, "up_id", input.up_id);
      write(xml, "down_id", input.down_id);
      write(xml, "strange_id", input.strange_id);
      write(xml, "charm_id", input.charm_id);
 
      pop(xml);
    }
 
 
    //! Propagator input
    void read(XMLReader& xml, const std::string& path, Params::NamedObject_t& input)
    {
      XMLReader inputtop(xml, path);
 
      read(inputtop, "gauge_id" , input.gauge_id);
      read(inputtop, "props"    , input.props);
    }
 
    //! Propagator output
    void write(XMLWriter& xml, const std::string& path, const Params::NamedObject_t& input)
    {
      push(xml, path);
      write(xml, "gauge_id" , input.gauge_id);
      write(xml, "props"    , input.props);
      pop(xml);
    }
 
 
    // Param stuff
    Params::Params()
    { 
      frequency = 0; 
    }
 
    Params::Params(XMLReader& xml_in, const std::string& path) 
    {
      try 
      {
        XMLReader paramtop(xml_in, path);
 
        if (paramtop.count("Frequency") == 1)
          read(paramtop, "Frequency", frequency);
        else
          frequency = 1;
 
        // Parameters for source construction
        read(paramtop, "Param", param);
 
        // Read in the output propagator/source configuration info
        read(paramtop, "NamedObject", named_obj);
 
        // Possible alternate XML file pattern
        if (paramtop.count("xml_file") != 0) 
        {
          read(paramtop, "xml_file", xml_file);
        }
      }
      catch(const std::string& e) 
      {
        QDPIO::cerr << __func__ << ": Caught Exception reading XML: " << e << std::endl;
        QDP_abort(1);
      }
    }
 
 
    void
    Params::writeXML(XMLWriter& xml_out, const std::string& path) 
    {
      push(xml_out, path);
    
      write(xml_out, "Param", param);
      write(xml_out, "NamedObject", named_obj);
      write(xml_out, "xml_file", xml_file);
 
      pop(xml_out);
    }
 
 
    void epsilon_contract(LatticeComplex& res,
                          const multi2d<LatticeComplex>& l,
                          const multi2d<LatticeComplex>& m,
                          const multi2d<LatticeComplex>& r){
      //LatticeComplex res = 0.0;
    
      //generated by : make_epsilon_contr.pl | sort                
 
      res  = l(0,0)*m(1,1)*r(2,2);
      res -= l(0,0)*m(1,2)*r(2,1);
      res -= l(0,0)*m(2,1)*r(1,2);
      res += l(0,0)*m(2,2)*r(1,1);
      res -= l(0,1)*m(1,0)*r(2,2);
      res += l(0,1)*m(1,2)*r(2,0);
      res += l(0,1)*m(2,0)*r(1,2);
      res -= l(0,1)*m(2,2)*r(1,0);
      res += l(0,2)*m(1,0)*r(2,1);
      res -= l(0,2)*m(1,1)*r(2,0);
      res -= l(0,2)*m(2,0)*r(1,1);
      res += l(0,2)*m(2,1)*r(1,0);
      res -= l(1,0)*m(0,1)*r(2,2);
      res += l(1,0)*m(0,2)*r(2,1);
      res += l(1,0)*m(2,1)*r(0,2);
      res -= l(1,0)*m(2,2)*r(0,1);
      res += l(1,1)*m(0,0)*r(2,2);
      res -= l(1,1)*m(0,2)*r(2,0);
      res -= l(1,1)*m(2,0)*r(0,2);
      res += l(1,1)*m(2,2)*r(0,0);
      res -= l(1,2)*m(0,0)*r(2,1);
      res += l(1,2)*m(0,1)*r(2,0);
      res += l(1,2)*m(2,0)*r(0,1);
      res -= l(1,2)*m(2,1)*r(0,0);
      res += l(2,0)*m(0,1)*r(1,2);
      res -= l(2,0)*m(0,2)*r(1,1);
      res -= l(2,0)*m(1,1)*r(0,2);
      res += l(2,0)*m(1,2)*r(0,1);
      res -= l(2,1)*m(0,0)*r(1,2);
      res += l(2,1)*m(0,2)*r(1,0);
      res += l(2,1)*m(1,0)*r(0,2);
      res -= l(2,1)*m(1,2)*r(0,0);
      res += l(2,2)*m(0,0)*r(1,1);
      res -= l(2,2)*m(0,1)*r(1,0);
      res -= l(2,2)*m(1,0)*r(0,1);
      res += l(2,2)*m(1,1)*r(0,0);
 
      //return res ;
    }
 
    // Anonymous namespace
    namespace 
    {
      //! Useful structure holding sink props
      class SinkPropContainer_t{
      public:
        ForwardProp_t prop_header;
        std::string quark_propagator_id;
        Real Mass;
 
        bool exists ;
 
        multi1d<int> bc; 
    
        std::string source_type;
        std::string source_disp_type;
        std::string sink_type;
        std::string sink_disp_type;
 
        SinkPropContainer_t(){
          exists=false ;
        }
        //! Read a sink prop
        void readSinkProp(const std::string& id)
          {
            if(id=="NULL"){
              return ;
            }
            try
            {
              // Try a cast to see if it succeeds
              const LatticePropagator& foo = 
                TheNamedObjMap::Instance().getData<LatticePropagator>(id);
            
              // Snarf the data into a copy
              quark_propagator_id = id;
            
              exists = true ;
              // Snarf the prop info. This is will throw if the prop_id is not there
              XMLReader prop_file_xml, prop_record_xml;
              TheNamedObjMap::Instance().get(id).getFileXML(prop_file_xml);
              TheNamedObjMap::Instance().get(id).getRecordXML(prop_record_xml);
            
              // Try to invert this record XML into a ChromaProp struct
              // Also pull out the id of this source
              {
                std::string xpath;
                read(prop_record_xml, "/SinkSmear", prop_header);
              
                read(prop_record_xml, "/SinkSmear/PropSource/Source/SourceType", source_type);
                xpath = "/SinkSmear/PropSource/Source/Displacement/DisplacementType";
                if (prop_record_xml.count(xpath) != 0)
                  read(prop_record_xml, xpath, source_disp_type);
                else
                  source_disp_type = NoQuarkDisplacementEnv::getName();
              
                read(prop_record_xml, "/SinkSmear/PropSink/Sink/SinkType", sink_type);
                xpath = "/SinkSmear/PropSink/Sink/Displacement/DisplacementType";
                if (prop_record_xml.count(xpath) != 0)
                  read(prop_record_xml, xpath, sink_disp_type);
                else
                  sink_disp_type = NoQuarkDisplacementEnv::getName();
              }
            }
            catch( std::bad_cast ) 
            {
              QDPIO::cerr << name << ": caught dynamic cast error" 
                          << std::endl;
              QDP_abort(1);
            }
            catch (const std::string& e) 
            {
              QDPIO::cerr << name << ": error message: " << e 
                          << std::endl;
              QDP_abort(1);
            }
        
        
            // Derived from input prop
            // Hunt around to find the mass
            // NOTE: this may be problematic in the future if actions are used with no
            // clear def. of a Mass
            QDPIO::cout << "Try action and mass" << std::endl;
            Mass = getMass(prop_header.prop_header.fermact);
        
            // Only baryons care about boundaries
            // Try to find them. If not present, assume dirichlet.
            // This turns off any attempt to time reverse which is the
            // only thing that the BC are affecting.
            bc.resize(Nd);
            bc = 0;
        
            try
            {
              bc = getFermActBoundary(prop_header.prop_header.fermact);
            }
            catch (const std::string& e) 
            {
              QDPIO::cerr << name 
                          << ": caught exception. No BC found in these headers. Will assume dirichlet: " << e 
                          << std::endl;
            }
        
            QDPIO::cout << "FermAct = " << prop_header.prop_header.fermact.id << std::endl;
            QDPIO::cout << "Mass = " << Mass << std::endl;
          }
      
      };
 
 
 
 
 
      void sanity_check_props(const SinkPropContainer_t& p1, 
                              const SinkPropContainer_t& p2){
      
        int j_decay = p1.prop_header.source_header.j_decay;
        if (p2.prop_header.source_header.j_decay != 
            p1.prop_header.source_header.j_decay){
          QDPIO::cerr << "Error!! j_decay must be the same for all propagators " << std::endl;
          QDP_abort(1);
        }
        if (p2.prop_header.source_header.t_source != 
            p1.prop_header.source_header.t_source)
        {
          QDPIO::cerr << "Error!! t_source must be the same for all propagators " << std::endl;
          QDP_abort(1);
        }
        if (p1.source_type != p2.source_type)
        {
          QDPIO::cerr << "Error!! source_type must be the same in a pair " << std::endl;
          QDP_abort(1);
        }
        if (p1.sink_type != p2.sink_type)
        {
          QDPIO::cerr << "Error!! source_type must be the same in a pair " << std::endl;
          QDP_abort(1);
        }
      
        if (p2.bc[j_decay] != p1.bc[j_decay])
        {
          QDPIO::cerr << "Error!! bc must be the same for all propagators " << std::endl;
          QDP_abort(1);
        }
      }
 
 
 
 
      //! Useful structure holding sink props
      struct AllSinkProps_t
      {
        int j_decay;
        int t0 ; 
        multi1d<int> t_srce ;
        int bc_spec ;
        std::map<std::string,SinkPropContainer_t>  prop;
        std::map<std::string,BarSpec::RPropagator> rprop;
 
        //! Read all sinks
        AllSinkProps_t(const Params::NamedObject_t::Props_t& p){
 
          QDPIO::cout<<"Attempt to parse forward propagator= "<<p.up_id<<std::endl;
          prop["up"].readSinkProp(p.up_id);
          rprop[p.up_id].ConvertProp(TheNamedObjMap::Instance().getData<LatticePropagator>(p.up_id));
          QDPIO::cout<<"up quark  propagator successfully parsed" << std::endl;
          j_decay = prop["up"].prop_header.source_header.j_decay;
          t0      = prop["up"].prop_header.source_header.t_source;
          t_srce  = prop["up"].prop_header.source_header.getTSrce();
          bc_spec = prop["up"].bc[j_decay];
 
          QDPIO::cout<<"Attempt to parse forward propagator= " <<p.down_id<<std::endl;
          //Always need a down quark 
          prop["down"].readSinkProp(p.down_id);
          QDPIO::cout << "down quark propagator successfully parsed" << std::endl;
          if(rprop.find(p.down_id) == rprop.end()){
            QDPIO::cout<<__func__<<": Need to convert prop id: "
                       <<p.down_id<<std::endl;
            rprop[p.down_id].ConvertProp(TheNamedObjMap::Instance().getData<LatticePropagator>(p.down_id));
          }
          QDPIO::cout<<"Attempt to parse forward propagator= "<<p.strange_id<<std::endl;
          prop["strange"].readSinkProp(p.strange_id);
          if(p.strange_id != "NULL"){
            QDPIO::cout <<"strange quark propagator successfully parsed" << std::endl;
            if(rprop.find(p.strange_id) == rprop.end()){
              QDPIO::cout<<__func__<<": Need to convert prop id: "
                         <<p.strange_id<<std::endl;
              rprop[p.strange_id].ConvertProp(TheNamedObjMap::Instance().getData<LatticePropagator>(p.strange_id));
            }
          }
 
          QDPIO::cout<<"Attempt to parse forward propagator= "<<p.charm_id<<std::endl;
          prop["charm"].readSinkProp(p.charm_id);
          if(p.charm_id != "NULL"){
            QDPIO::cout << "charm quark propagator successfully parsed" << std::endl;
            if(rprop.find(p.charm_id) == rprop.end()){
              QDPIO::cout<<__func__<<": Need to convert prop id: "
                         <<p.charm_id<<std::endl;
              rprop[p.charm_id].ConvertProp(TheNamedObjMap::Instance().getData<LatticePropagator>(p.charm_id));
            }
          }
        
         
          sanity_check_props(prop["up"], prop["down"]) ;
          if(prop["strange"].exists)
            sanity_check_props(prop["up"], prop["strange"]) ;
          if(prop["charm"].exists)
            sanity_check_props(prop["up"], prop["charm"]) ;
 
        }
 
 
        std::string sink(const std::string& flavor){
          return prop[flavor].sink_type ;
        }
 
        std::string source(const std::string& flavor){
          return prop[flavor].source_type ;
        }
 
        const BarSpec::RPropagator& prop_ref(const std::string& flavor){
 
          return  rprop[prop[flavor].quark_propagator_id] ;
      
        }
 
 
      };
 
 
    
 
      void xml_print_prop_info(XMLWriter& xml_out, 
                               const SinkPropContainer_t& p,
                               const Set& s){
      
        // Sanity check - write out the norm2 of the forward prop 
        // Use this for any possible verification
        if(p.exists){
          push(xml_out, p.quark_propagator_id);
          {
            const LatticePropagator& prop = 
              TheNamedObjMap::Instance().getData<LatticePropagator>(p.quark_propagator_id);
          
            QDPIO::cout << "propagator_id = " << p.quark_propagator_id << std::endl;
            QDPIO::cout << "Source_type   = " << p.source_type << std::endl;
            QDPIO::cout << "Sink_type     = " << p.sink_type << std::endl;
          
            write(xml_out, "correlator",sumMulti(localNorm2(prop),s));
            write(xml_out, "source_type", p.source_type);
            write(xml_out, "source_disp_type", p.source_disp_type);
            write(xml_out, "sink_type", p.sink_type);
            write(xml_out, "sink_disp_type", p.sink_disp_type);
            write(xml_out, "Mass", p.Mass);
 
          }
          pop(xml_out);
        }
      }
 
 
 
    } // namespace anonymous
 
    namespace BarSpec{
    
      std::vector<int> permutation(int k, const std::vector<int>& s){
        std::vector<int> p = s ;
        int ss ;
        int jj ;
        for(int j(2); j<=s.size(); j++){
          jj=k%j ;
          ss=p[jj]; p[jj]=p[j-1]; p[j-1]=ss;
          k=k/j ;
        }
        return p ;
      }
    
      int permutation_sign(int k, int n){
        int jj ;
        int sign_(1) ;
        for(int j(2); j<=n; j++){
          jj=k%j ;
          if(jj!=(j-1))
            sign_ *= (-1) ;
          k=k/j ;
        }
        return sign_ ;
      }
 
      void contract(LatticeComplex& latC, 
                    const RPropagator& q1,
                    const RPropagator& q2,
                    const RPropagator& q3,
                    const SpinWF_t& snk,
                    const SpinWF_t& src){
        LatticeComplex cc ;
        latC = 0.0;
        for(int s(0);s<src.terms.size();s++)
          for(int ss(0);ss<snk.terms.size();ss++){
            epsilon_contract(cc,
                             q1.p(snk.terms[ss].spin[0],src.terms[s].spin[0]),
                             q2.p(snk.terms[ss].spin[1],src.terms[s].spin[1]),
                             q3.p(snk.terms[ss].spin[2],src.terms[s].spin[2]) );
            latC += (snk.terms[ss].weight*src.terms[s].weight)*cc ;
          }// loop over source sink wavefunction components
        latC *= (snk.norm*src.norm) ;
      }
    }//barspec name space
 
 
    // Function call
    void 
    InlineMeas::operator()(unsigned long update_no,
                           XMLWriter& xml_out) 
    {
      // If xml file not empty, then use alternate
      if (params.xml_file != "")
      {
        std::string xml_file = makeXMLFileName(params.xml_file, update_no);
 
        push(xml_out, "barspec");
        write(xml_out, "update_no", update_no);
        write(xml_out, "xml_file", xml_file);
        pop(xml_out);
 
        XMLFileWriter xml(xml_file);
        func(update_no, xml);
      }
      else
      {
        func(update_no, xml_out);
      }
    }
 
 
    // Real work done here
    void 
    InlineMeas::func(unsigned long update_no,
                     XMLWriter& xml_out) 
    {
      START_CODE();
 
      StopWatch snoop;
      snoop.reset();
      snoop.start();
 
      // Test and grab a reference to the gauge field
      XMLBufferWriter gauge_xml;
      try
      {
        TheNamedObjMap::Instance().getData< multi1d<LatticeColorMatrix> >(params.named_obj.gauge_id);
        TheNamedObjMap::Instance().get(params.named_obj.gauge_id).getRecordXML(gauge_xml);
      }
      catch( std::bad_cast ) 
      {
        QDPIO::cerr << name << ": caught dynamic cast error" 
                    << std::endl;
        QDP_abort(1);
      }
      catch (const std::string& e) 
      {
        QDPIO::cerr << name << ": std::map call failed: " << e 
                    << std::endl;
        QDP_abort(1);
      }
      const multi1d<LatticeColorMatrix>& u = 
        TheNamedObjMap::Instance().getData< multi1d<LatticeColorMatrix> >(params.named_obj.gauge_id);
 
      push(xml_out, "barspec");
      write(xml_out, "update_no", update_no);
 
      QDPIO::cout << " BARSPEC: Spectroscopy for Wilson-like fermions" << std::endl;
      QDPIO::cout << std::endl << "     Gauge group: SU(" << Nc << ")" << std::endl;
      QDPIO::cout << "     volume: " << Layout::lattSize()[0];
      for (int i=1; i<Nd; ++i) {
        QDPIO::cout << " x " << Layout::lattSize()[i];
      }
      QDPIO::cout << std::endl;
 
      proginfo(xml_out);    // Print out basic program info
 
      // Write out the input
      params.writeXML(xml_out, "Input");
 
      // Write out the config info
      write(xml_out, "Config_info", gauge_xml);
 
      push(xml_out, "Output_version");
      write(xml_out, "out_version", 1);
      pop(xml_out);
 
 
      // First calculate some gauge invariant observables just for info.
      MesPlq(xml_out, "Observables", u);
 
      //open the database
      //
      BinaryStoreDB< SerialDBKey<KeyHadron2PtCorr_t>, SerialDBData<multi1d<ComplexD> > > qdp_db;
 
 
      // Now loop over the various fermion pairs
      for(int p=0; p < params.named_obj.props.size(); ++p)
      {
        const Params::NamedObject_t::Props_t named_obj = params.named_obj.props[p];
 
        push(xml_out, "elem");
 
        AllSinkProps_t all_sinks(params.named_obj.props[p]);
 
        // Derived from input prop
        const multi1d<int>& t_srce = all_sinks.t_srce ;
        const int& j_decay         = all_sinks.j_decay ;
        const int& t0              = all_sinks.t0 ;
 
        const int& bc_spec         = all_sinks.bc_spec ;
     
 
        // Initialize the slow Fourier transform phases
        SftMom phases(params.param.mom2_max, t_srce, params.param.avg_equiv_mom,
                      j_decay);
 
        // Keep a copy of the phases with NO momenta
        SftMom phases_nomom(0, true, j_decay);
 
        write(xml_out, "t0", t0);
 
        // Save prop input
        push(xml_out, "Forward_prop_headers");
        write(xml_out, "up_prop", all_sinks.prop["up"].prop_header);
        write(xml_out, "down_prop", all_sinks.prop["down"].prop_header);
        if(all_sinks.prop["strange"].exists)
          write(xml_out, "strange_prop", all_sinks.prop["strange"].prop_header);
        if(all_sinks.prop["charm"].exists)
          write(xml_out, "charm_prop", all_sinks.prop["charm"].prop_header);
 
        pop(xml_out);
 
        push(xml_out, "Forward_Propagator_Properties"                     );
        std::map<std::string,SinkPropContainer_t>::iterator it;
        for(it=all_sinks.prop.begin();it!=all_sinks.prop.end();it++)
          xml_print_prop_info(xml_out, it->second      , phases.getSet()  );
      
        pop(xml_out);
 
        int Nt = Layout::lattSize()[j_decay];
 
        LatticeComplex latC ;
        StopWatch tictoc;
        tictoc.reset();
        tictoc.start();
        for(int s(0);s<params.param.states.size();s++){
          QDPIO::cout<<"Doing state "<<params.param.states[s].name<<std::endl;
          QDPIO::cout<<"   Flavor structure: " ;
          for(int k(0);k<params.param.states[s].flavor.size();k++)
            QDPIO::cout<<"  "<<params.param.states[s].flavor[k] ;
          QDPIO::cout<<std::endl;
          multi1d<std::string> prop_id(params.param.states[s].flavor.size()) ;
          for(int k(0);k<params.param.states[s].flavor.size();k++){
            prop_id[k] = params.param.states[s].flavor[k] ;
          }
          // Open the file, and write the meta-data and the binary for this state
          if (! qdp_db.fileExists(params.param.states[s].db)){
            XMLBufferWriter file_xml;
          
            push(file_xml, "DBMetaData");
            write(file_xml, "id", std::string("hadron2Pt"));
            write(file_xml, "lattSize", QDP::Layout::lattSize());
            write(file_xml, "decay_dir", j_decay);
            write(file_xml, "State", params.param.states[s].name);
            proginfo(file_xml);    // Print out basic program info
            write(file_xml, "Params", params.param);
            write(file_xml, "Config_info", gauge_xml);
            pop(file_xml);
          
            std::string file_str(file_xml.str());
            qdp_db.setMaxUserInfoLen(file_str.size());    
            qdp_db.open(params.param.states[s].db, O_RDWR | O_CREAT, 0664);
            qdp_db.insertUserdata(file_str);
          }
          else
          {
            qdp_db.open(params.param.states[s].db, O_RDWR, 0664);
          }
 
 
          // References for use later 
          const BarSpec::RPropagator& q1 = all_sinks.prop_ref(prop_id[0]) ;
          const BarSpec::RPropagator& q2 = all_sinks.prop_ref(prop_id[1]) ;
          const BarSpec::RPropagator& q3 = all_sinks.prop_ref(prop_id[2]) ;
 
          KeyHadron2PtCorr_t key  ;
 
          std::ostringstream os ;
          os<<all_sinks.prop[prop_id[0]].Mass<<" ";
          os<<all_sinks.prop[prop_id[1]].Mass<<" ";
          os<<all_sinks.prop[prop_id[2]].Mass ;
          key.mass        = os.str();
          key.ensemble    = params.param.ensemble;
          key.num_vecs = Nc;  /*!< Number of color vectors */
 
          key.src_smear    = all_sinks.source(prop_id[0])+" ";
          key.src_smear   += all_sinks.source(prop_id[1])+" ";
          key.src_smear   += all_sinks.source(prop_id[2]) ;         
          key.src_spin  =   params.param.states[s].spin ;
          key.src_lorentz.resize(0);
          //key.src_lorentz =  key.src_spin  ;
 
          key.snk_smear    = all_sinks.sink(prop_id[0])+" ";
          key.snk_smear   += all_sinks.sink(prop_id[1])+" ";
          key.snk_smear   += all_sinks.sink(prop_id[2]) ;
          key.snk_spin  =  params.param.states[s].spin ;
          key.snk_lorentz.resize(0);
          //key.snk_lorentz =  key.snk_spin  ;
 
          //loop over momenta goes here
          for(int oi(0);oi<params.param.states[s].ops.size();oi++){ //sink
            BarSpec::SpinWF_t snk(params.param.states[s].ops[oi].spinWF);
            snk.permutations(prop_id);
            for(int oj(0);oj<params.param.states[s].ops.size();oj++){//source
 
              BarSpec::SpinWF_t src(params.param.states[s].ops[oj].spinWF);
 
              key.src_name    = params.param.states[s].ops[oj].name;
              key.snk_name    = params.param.states[s].ops[oi].name;
            
              BarSpec::contract(latC,q1,q2,q3,snk,src);
 
              multi2d<ComplexD> hsum;
              hsum = phases.sft(latC) ;
            
              for(int mom(0);mom<phases.numMom();mom++){
                key.mom = phases.numToMom(mom);    /*<! Momentum  */
                SerialDBKey<KeyHadron2PtCorr_t> K;
                K.key() = key ;
                SerialDBData<multi1d<ComplexD> >  V ;
                V.data().resize(Nt);
                for(int t(0);t<Nt;t++){
                  int t_eff = (t - t0 + Nt) % Nt;
                  if ( bc_spec < 0 && (t_eff+t0) >= Nt)
                    V.data()[t_eff] = -hsum[mom][t];
                  else
                    V.data()[t_eff] =  hsum[mom][t];
                }//loop over time
                qdp_db.insert(K,V);
              }// loop over momenta
 
            }// loop  over source ops
          }// loop over sink ops
          qdp_db.close() ;
        }// loop over states
        tictoc.stop();
        QDPIO::cout << name << ": contraction time = "
                    << tictoc.getTimeInSeconds() 
                    << " secs" << std::endl;
 
        pop(xml_out); //elem
      }// loop over propagator groups
      pop(xml_out);  // barspec
 
      snoop.stop();
      QDPIO::cout << name << ": total time = "
                  << snoop.getTimeInSeconds() 
                  << " secs" << std::endl;
 
      QDPIO::cout << name << ": ran successfully" << std::endl;
 
      END_CODE();
    } 
 
  }
 
}