diff --git a/Grid/qcd/action/momentum/DDHMCfilter.h b/Grid/qcd/action/momentum/DDHMCfilter.h index 27e0ac66..c96351e2 100644 --- a/Grid/qcd/action/momentum/DDHMCfilter.h +++ b/Grid/qcd/action/momentum/DDHMCfilter.h @@ -28,6 +28,7 @@ directory *************************************************************************************/ /* END LEGAL */ +NAMESPACE_BEGIN(Grid); //////////////////////////////////////////////////// // DDHMC filter with sub-block size B[mu] //////////////////////////////////////////////////// @@ -43,37 +44,60 @@ struct DDHMCFilter: public MomentumFilterBase typedef Lattice LatticeScalarType; Coordinate Block; + int Width; // Could also zero links in plane like Luscher advocates. - DDHMCFilter(const Coordinate &_Block): Block(_Block){} + DDHMCFilter(const Coordinate &_Block,int _Width): Block(_Block), Width(_Width) + { + assert( (Width==0) || (Width==1) || (Width==2) ); + } void applyFilter(MomentaField &P) const override { GridBase *grid = P.Grid(); - LatticeScalarType mask_mu(grid); - LatticeLorentzScalarType mask(grid); + + LatticeColourMatrix zz(grid); zz = Zero(); //////////////////////////////////////////////////// // Zero strictly links crossing between domains // Luscher also zeroes links in plane of domain boundaries // Keeping interior only. This prevents force from plaquettes // crossing domains and keeps whole MD trajectory local. - // Step further than I want to go. + // This is the case Width=1. + // Width = 2 should also decouple rectangle force. //////////////////////////////////////////////////// - ScalarType zz = scalar_type(0.0); - ScalarType one= scalar_type(1.0); - + Coordinate Global=grid->GlobalDimensions(); LatticeInteger coor(grid); for(int mu=0;mu(P, mu); + P_mu = where(mod(coor,Block[mu])==Integer(Block[mu]-1),zz,P_mu); + /* if(Width>=2) { + P_mu = where(mod(coor,Block[mu])==Integer(0),zz,P_mu); + P_mu = where(mod(coor,Block[mu])==Integer(Block[mu]-2),zz,P_mu); + }*/ + PokeIndex(P, P_mu, mu); - mask_mu = where(mod(coor,Block[mu])==Block[mu]-1,zz,one); - - PokeIndex(mask, mask_mu, mu); - + for(int nu=0;nu(P, nu); + if ( Width>=1 ){ + P_nu = where(mod(coor,Block[mu])==Integer(Block[mu]-1),zz,P_nu); + P_nu = where(mod(coor,Block[mu])==Integer(0) ,zz,P_nu); + } + /* if ( Width>=2 ){ + P_nu = where(mod(coor,Block[mu])==Integer(Block[mu]-2),zz,P_nu); + P_nu = where(mod(coor,Block[mu])==Integer(1) ,zz,P_nu); + }*/ + PokeIndex(P, P_nu, nu); + } + } + } } } }; +NAMESPACE_END(Grid);