diff --git a/Grid/qcd/action/fermion/DomainWallFermion.h b/Grid/qcd/action/fermion/DomainWallFermion.h index 56179f26..51c2f0ba 100644 --- a/Grid/qcd/action/fermion/DomainWallFermion.h +++ b/Grid/qcd/action/fermion/DomainWallFermion.h @@ -43,7 +43,7 @@ namespace Grid { INHERIT_IMPL_TYPES(Impl); public: - void FreePropagator(const FermionField &in,FermionField &out,RealD mass, std::vector twist, bool fiveD) { + void FreePropagator(const FermionField &in,FermionField &out,RealD mass,std::vector boundary, std::vector twist, bool fiveD) { FermionField in_k(in._grid); FermionField prop_k(in._grid); @@ -55,15 +55,19 @@ namespace Grid { FermionField in_buf(in._grid); in_buf = zero; Complex ci(0.0,1.0); assert(twist.size() == Nd);//check that twist is Nd + assert(boundary.size() == Nd);//check that boundary conditions is Nd int shift = 0; if(fiveD) shift = 1; for(unsigned int nu = 0; nu < Nd; nu++) { // Shift coordinate lattice index by 1 to account for 5th dimension. LatticeCoordinate(coor, nu + shift); - ph = ph + twist[nu]*coor*((1./(in._grid->_fdimensions[nu+shift]))); + double boundary_phase = ::acos(real(boundary[nu])); + ph = ph + boundary_phase*coor*((1./(in._grid->_fdimensions[nu+shift]))); + //momenta for propagator shifted by twist+boundary + twist[nu] = twist[nu] + boundary_phase/((2.0*M_PI)); } - in_buf = exp((Real)(2.0*M_PI)*ci*ph*(-1.0))*in; + in_buf = exp(ci*ph*(-1.0))*in; if(fiveD){//FFT only on temporal and spatial dimensions std::vector mask(Nd+1,1); mask[0] = 0; @@ -76,25 +80,28 @@ namespace Grid { this->MomentumSpacePropagatorHt(prop_k,in_k,mass,twist); theFFT.FFT_all_dim(out,prop_k,FFT::backward); } - //phase for boundary condition - out = out * exp((Real)(2.0*M_PI)*ci*ph); + out = out * exp(ci*ph); }; - virtual void FreePropagator(const FermionField &in,FermionField &out,RealD mass,std::vector twist) { + virtual void FreePropagator(const FermionField &in,FermionField &out,RealD mass,std::vector boundary,std::vector twist) { bool fiveD = true; //5d propagator by default - FreePropagator(in,out,mass,twist,fiveD); + FreePropagator(in,out,mass,boundary,twist,fiveD); }; virtual void FreePropagator(const FermionField &in,FermionField &out,RealD mass, bool fiveD) { std::vector twist(Nd,0.0); //default: periodic boundarys in all directions - FreePropagator(in,out,mass,twist,fiveD); + std::vector boundary; + for(int i=0;i twist(Nd,0.0); //default: periodic boundarys in all directions - FreePropagator(in,out,mass,twist,fiveD); + std::vector twist(Nd,0.0); //default: twist angle 0 + std::vector boundary; + for(int i=0;i twist) { assert(0);}; - virtual void FreePropagator(const FermionField &in,FermionField &out,RealD mass,std::vector twist) { + virtual void FreePropagator(const FermionField &in,FermionField &out,RealD mass,std::vector boundary,std::vector twist) { FFT theFFT((GridCartesian *) in._grid); FermionField in_k(in._grid); @@ -108,24 +108,31 @@ namespace Grid { FermionField in_buf(in._grid); in_buf = zero; Complex ci(0.0,1.0); assert(twist.size() == Nd);//check that twist is Nd + assert(boundary.size() == Nd);//check that boundary conditions is Nd for(unsigned int nu = 0; nu < Nd; nu++) { LatticeCoordinate(coor, nu); - ph = ph + twist[nu]*coor*((1./(in._grid->_fdimensions[nu]))); + double boundary_phase = ::acos(real(boundary[nu])); + ph = ph + boundary_phase*coor*((1./(in._grid->_fdimensions[nu]))); + //momenta for propagator shifted by twist+boundary + twist[nu] = twist[nu] + boundary_phase/((2.0*M_PI)); } - in_buf = exp((Real)(2.0*M_PI)*ci*ph*(-1.0))*in; + in_buf = exp(ci*ph*(-1.0))*in; theFFT.FFT_all_dim(in_k,in_buf,FFT::forward); this->MomentumSpacePropagator(prop_k,in_k,mass,twist); theFFT.FFT_all_dim(out,prop_k,FFT::backward); //phase for boundary condition - out = out * exp((Real)(2.0*M_PI)*ci*ph); + out = out * exp(ci*ph); }; + virtual void FreePropagator(const FermionField &in,FermionField &out,RealD mass) { + std::vector boundary; + for(int i=0;i twist(Nd,0.0); //default: periodic boundarys in all directions - FreePropagator(in,out,mass,twist); + FreePropagator(in,out,mass,boundary,twist); }; /////////////////////////////////////////////// diff --git a/Hadrons/Modules/MFermion/EMLepton.hpp b/Hadrons/Modules/MFermion/EMLepton.hpp index a3f5fa64..a67f7495 100644 --- a/Hadrons/Modules/MFermion/EMLepton.hpp +++ b/Hadrons/Modules/MFermion/EMLepton.hpp @@ -72,6 +72,7 @@ public: std::string, action, std::string, emField, double, mass, + std::string , boundary, std::string, twist, unsigned int, deltat); }; @@ -163,8 +164,7 @@ void TEMLepton::execute(void) envGetTmp(FermionField, sol); envGetTmp(FermionField, tmp); LOG(Message) << "Calculating a lepton Propagator with sequential Aslash insertion with lepton mass " - << mass << " and twist (" - << par().twist << ") using the action '" << par().action + << mass << " using the action '" << par().action << "' for fixed source-sink separation of " << par().deltat << std::endl; envGetTmp(Lattice>, tlat); @@ -176,6 +176,11 @@ void TEMLepton::execute(void) { HADRONS_ERROR(Size, "number of twist angles does not match number of dimensions"); } + std::vector boundary = strToVec(par().boundary); + if(boundary.size() != Nd) + { + HADRONS_ERROR(Size, "number of boundary conditions does not match number of dimensions"); + } auto &stoch_photon = envGet(EmField, par().emField); unsigned int nt = env().getDim(Tp); @@ -212,7 +217,7 @@ void TEMLepton::execute(void) mat.ImportPhysicalFermionSource(tmp, source); } sol = zero; - mat.FreePropagator(source,sol,mass,twist); + mat.FreePropagator(source,sol,mass,boundary,twist); if (Ls_ == 1) { FermToProp(freetmp, sol, s, 0); @@ -229,11 +234,7 @@ void TEMLepton::execute(void) //shift free propagator to different source positions proptmp = Cshift(freetmp,Tp, -tl); - //take anti-periodic boundary conditions into account, if used - if(twist[Tp]==0.5) - { - proptmp = where( tlat < tl, (-1.0)*proptmp, proptmp); - } + proptmp = where( tlat < tl, boundary[Tp]*proptmp, proptmp); // free propagator for fixed source-sink separation lep = where(tlat == (tl-par().deltat+nt)%nt, proptmp, lep); @@ -264,7 +265,7 @@ void TEMLepton::execute(void) mat.ImportPhysicalFermionSource(tmp, source); } sol = zero; - mat.FreePropagator(source,sol,mass,twist); + mat.FreePropagator(source,sol,mass,boundary,twist); if (Ls_ == 1) { FermToProp(proptmp, sol, s, 0); diff --git a/Hadrons/Modules/MFermion/FreeProp.hpp b/Hadrons/Modules/MFermion/FreeProp.hpp index 4485b6f9..baccba6a 100644 --- a/Hadrons/Modules/MFermion/FreeProp.hpp +++ b/Hadrons/Modules/MFermion/FreeProp.hpp @@ -49,6 +49,7 @@ public: std::string, source, std::string, action, double, mass, + std::string , boundary, std::string, twist); }; @@ -172,7 +173,12 @@ void TFreeProp::execute(void) { HADRONS_ERROR(Size, "number of twist angles does not match number of dimensions"); } - mat.FreePropagator(source,sol,mass,twist); + std::vector boundary = strToVec(par().boundary); + if(boundary.size() != Nd) + { + HADRONS_ERROR(Size, "number of boundary conditions does not match number of dimensions"); + } + mat.FreePropagator(source,sol,mass,boundary,twist); FermToProp(prop, sol, s, c); // create 4D propagators from 5D one if necessary if (Ls_ > 1)