diff --git a/Grid/qcd/action/fermion/FermionOperator.h b/Grid/qcd/action/fermion/FermionOperator.h index 281c077e..55942dea 100644 --- a/Grid/qcd/action/fermion/FermionOperator.h +++ b/Grid/qcd/action/fermion/FermionOperator.h @@ -136,6 +136,7 @@ namespace Grid { ////////////////////////////////////////////////////////////////////// // Conserved currents, either contract at sink or insert sequentially. ////////////////////////////////////////////////////////////////////// + virtual void ContractConservedCurrent(PropagatorField &q_in_1, PropagatorField &q_in_2, PropagatorField &q_out, @@ -148,6 +149,12 @@ namespace Grid { unsigned int tmin, unsigned int tmax, ComplexField &lattice_cmplx)=0; + + // Only reimplemented in Wilson5D + // Default to just a zero correlation function + virtual void ContractJ5q(FermionField &q_in ,ComplexField &J5q) { J5q=zero; }; + virtual void ContractJ5q(PropagatorField &q_in,ComplexField &J5q) { J5q=zero; }; + /////////////////////////////////////////////// // Physical field import/export /////////////////////////////////////////////// diff --git a/Grid/qcd/action/fermion/WilsonFermion5D.cc b/Grid/qcd/action/fermion/WilsonFermion5D.cc index 47c9cbb5..fb3fd861 100644 --- a/Grid/qcd/action/fermion/WilsonFermion5D.cc +++ b/Grid/qcd/action/fermion/WilsonFermion5D.cc @@ -939,6 +939,75 @@ void WilsonFermion5D::MomentumSpacePropagatorHw(FermionField &out,const Fe merge(qSiteRev, qSiteVec); \ } +// psi = chiralProjectPlus(Result_s[Ls/2-1]); +// psi+= chiralProjectMinus(Result_s[Ls/2]); +// PJ5q+=localInnerProduct(psi,psi); + +template +Lattice spProj5p(const Lattice & in) +{ + GridBase *grid=in._grid; + Gamma G5(Gamma::Algebra::Gamma5); + Lattice ret(grid); + parallel_for(int ss=0;ssoSites();ss++){ + ret._odata[ss] = in._odata[ss] + G5*in._odata[ss]; + } + return ret; +} +template +Lattice spProj5m(const Lattice & in) +{ + Gamma G5(Gamma::Algebra::Gamma5); + GridBase *grid=in._grid; + Lattice ret(grid); + parallel_for(int ss=0;ssoSites();ss++){ + ret._odata[ss] = in._odata[ss] - G5*in._odata[ss]; + } + return ret; +} + +template +void WilsonFermion5D::ContractJ5q(FermionField &q_in,ComplexField &J5q) +{ + conformable(GaugeGrid(), J5q._grid); + conformable(q_in._grid, FermionGrid()); + + // 4d field + int Ls = this->Ls; + FermionField psi(GaugeGrid()); + FermionField p_plus (GaugeGrid()); + FermionField p_minus(GaugeGrid()); + FermionField p(GaugeGrid()); + + ExtractSlice(p_plus , q_in, Ls/2 , 0); + ExtractSlice(p_minus, q_in, Ls/2-1 , 0); + p_plus = spProj5p(p_plus ); + p_minus= spProj5m(p_minus); + p=p_plus+p_minus; + J5q = localInnerProduct(p,p); +} + +template +void WilsonFermion5D::ContractJ5q(PropagatorField &q_in,ComplexField &J5q) +{ + conformable(GaugeGrid(), J5q._grid); + conformable(q_in._grid, FermionGrid()); + + // 4d field + int Ls = this->Ls; + PropagatorField psi(GaugeGrid()); + PropagatorField p_plus (GaugeGrid()); + PropagatorField p_minus(GaugeGrid()); + PropagatorField p(GaugeGrid()); + + ExtractSlice(p_plus , q_in, Ls/2 , 0); + ExtractSlice(p_minus, q_in, Ls/2-1 , 0); + p_plus = spProj5p(p_plus ); + p_minus= spProj5m(p_minus); + p=p_plus+p_minus; + J5q = localInnerProduct(p,p); +} + template void WilsonFermion5D::ContractConservedCurrent(PropagatorField &q_in_1, PropagatorField &q_in_2, @@ -949,6 +1018,7 @@ void WilsonFermion5D::ContractConservedCurrent(PropagatorField &q_in_1, conformable(q_in_1._grid, FermionGrid()); conformable(q_in_1._grid, q_in_2._grid); conformable(_FourDimGrid, q_out._grid); + PropagatorField tmp1(FermionGrid()), tmp2(FermionGrid()); unsigned int LLs = q_in_1._grid->_rdimensions[0]; q_out = zero; @@ -995,7 +1065,6 @@ void WilsonFermion5D::ContractConservedCurrent(PropagatorField &q_in_1, } - template void WilsonFermion5D::SeqConservedCurrent(PropagatorField &q_in, PropagatorField &q_out, diff --git a/Grid/qcd/action/fermion/WilsonFermion5D.h b/Grid/qcd/action/fermion/WilsonFermion5D.h index d22b5d6f..4a31bb43 100644 --- a/Grid/qcd/action/fermion/WilsonFermion5D.h +++ b/Grid/qcd/action/fermion/WilsonFermion5D.h @@ -230,6 +230,10 @@ namespace QCD { unsigned int tmin, unsigned int tmax, ComplexField &lattice_cmplx); + + void ContractJ5q(PropagatorField &q_in,ComplexField &J5q); + void ContractJ5q(FermionField &q_in,ComplexField &J5q); + }; }} diff --git a/Grid/qcd/utils/GaugeFix.h b/Grid/qcd/utils/GaugeFix.h index c4ea31aa..2b0cd7f2 100644 --- a/Grid/qcd/utils/GaugeFix.h +++ b/Grid/qcd/utils/GaugeFix.h @@ -53,21 +53,32 @@ class FourierAcceleratedGaugeFixer : public Gimpl { } static void SteepestDescentGaugeFix(GaugeLorentz &Umu,Real & alpha,int maxiter,Real Omega_tol, Real Phi_tol,bool Fourier=false) { GridBase *grid = Umu._grid; + GaugeMat xform(grid); + SteepestDescentGaugeFix(Umu,xform,alpha,maxiter,Omega_tol,Phi_tol,Fourier); + } + static void SteepestDescentGaugeFix(GaugeLorentz &Umu,GaugeMat &xform,Real & alpha,int maxiter,Real Omega_tol, Real Phi_tol,bool Fourier=false) { + + GridBase *grid = Umu._grid; Real org_plaq =WilsonLoops::avgPlaquette(Umu); Real org_link_trace=WilsonLoops::linkTrace(Umu); Real old_trace = org_link_trace; Real trG; + + xform=1.0; std::vector U(Nd,grid); - GaugeMat dmuAmu(grid); + + GaugeMat dmuAmu(grid); for(int i=0;i(Umu,mu); + if ( Fourier==false ) { - trG = SteepestDescentStep(U,alpha,dmuAmu); + trG = SteepestDescentStep(U,xform,alpha,dmuAmu); } else { - trG = FourierAccelSteepestDescentStep(U,alpha,dmuAmu); + trG = FourierAccelSteepestDescentStep(U,xform,alpha,dmuAmu); } for(int mu=0;mu(Umu,U[mu],mu); // Monitor progress and convergence test @@ -84,7 +95,6 @@ class FourierAcceleratedGaugeFixer : public Gimpl { Real Phi = 1.0 - old_trace / link_trace ; Real Omega= 1.0 - trG; - std::cout << GridLogMessage << " Iteration "< &U,Real & alpha, GaugeMat & dmuAmu) { + static Real SteepestDescentStep(std::vector &U,GaugeMat &xform,Real & alpha, GaugeMat & dmuAmu) { GridBase *grid = U[0]._grid; std::vector A(Nd,grid); @@ -109,12 +119,13 @@ class FourierAcceleratedGaugeFixer : public Gimpl { Real vol = grid->gSites(); Real trG = TensorRemove(sum(trace(g))).real()/vol/Nc; + xform = g*xform ; SU::GaugeTransform(U,g); return trG; } - static Real FourierAccelSteepestDescentStep(std::vector &U,Real & alpha, GaugeMat & dmuAmu) { + static Real FourierAccelSteepestDescentStep(std::vector &U,GaugeMat &xform,Real & alpha, GaugeMat & dmuAmu) { GridBase *grid = U[0]._grid; @@ -153,13 +164,6 @@ class FourierAcceleratedGaugeFixer : public Gimpl { Complex psqMax(16.0); Fp = psqMax*one/psq; - /* - static int once; - if ( once == 0 ) { - std::cout << " Fp " << Fp <::GaugeTransform(U,g); return trG; diff --git a/Grid/qcd/utils/SUn.h b/Grid/qcd/utils/SUn.h index cdc6c961..5cb564bd 100644 --- a/Grid/qcd/utils/SUn.h +++ b/Grid/qcd/utils/SUn.h @@ -676,10 +676,18 @@ class SU { } } /* - add GaugeTrans -*/ - -template + * Fundamental rep gauge xform + */ + template + static void GaugeTransformFundamental( Fundamental &ferm, GaugeMat &g){ + GridBase *grid = ferm._grid; + conformable(grid,g._grid); + ferm = g*ferm; + } +/* + * Adjoint rep gauge xform + */ + template static void GaugeTransform( GaugeField &Umu, GaugeMat &g){ GridBase *grid = Umu._grid; conformable(grid,g._grid); diff --git a/tests/core/Test_fft_gfix.cc b/tests/core/Test_fft_gfix.cc index 916c4b0b..7fdf5f73 100644 --- a/tests/core/Test_fft_gfix.cc +++ b/tests/core/Test_fft_gfix.cc @@ -63,8 +63,11 @@ int main (int argc, char ** argv) LatticeGaugeField Umu(&GRID); LatticeGaugeField Urnd(&GRID); LatticeGaugeField Uorg(&GRID); + LatticeGaugeField Utmp(&GRID); LatticeColourMatrix g(&GRID); // Gauge xform + LatticeColourMatrix xform1(&GRID); // Gauge xform + LatticeColourMatrix xform2(&GRID); // Gauge xform SU3::ColdConfiguration(pRNG,Umu); // Unit gauge Uorg=Umu; @@ -78,7 +81,14 @@ int main (int argc, char ** argv) Real alpha=0.1; Umu = Urnd; - FourierAcceleratedGaugeFixer::SteepestDescentGaugeFix(Umu,alpha,10000,1.0e-12, 1.0e-12,false); + FourierAcceleratedGaugeFixer::SteepestDescentGaugeFix(Umu,xform1,alpha,10000,1.0e-12, 1.0e-12,false); + + // Check the gauge xform matrices + Utmp=Urnd; + SU::GaugeTransform(Utmp,xform1); + Utmp = Utmp - Umu; + std::cout << " Norm Difference of xformed gauge "<< norm2(Utmp) << std::endl; + plaq=WilsonLoops::avgPlaquette(Umu); std::cout << " Final plaquette "<::SteepestDescentGaugeFix(Umu,alpha,10000,1.0e-12, 1.0e-12,true); + FourierAcceleratedGaugeFixer::SteepestDescentGaugeFix(Umu,xform2,alpha,10000,1.0e-12, 1.0e-12,true); + + Utmp=Urnd; + SU::GaugeTransform(Utmp,xform2); + Utmp = Utmp - Umu; + std::cout << " Norm Difference of xformed gauge "<< norm2(Utmp) << std::endl; + plaq=WilsonLoops::avgPlaquette(Umu); std::cout << " Final plaquette "<