From 04190ee7f30b181c851d40e01e1375ebee24cfc8 Mon Sep 17 00:00:00 2001 From: Vera Guelpers Date: Thu, 3 May 2018 12:31:36 +0100 Subject: [PATCH] 5D free propagator for DWF and boundary conditions for free propagators --- lib/qcd/action/fermion/DomainWallFermion.h | 55 ++++- lib/qcd/action/fermion/FermionOperator.h | 30 ++- .../fermion/OverlapWilsonCayleyTanhFermion.h | 4 +- lib/qcd/action/fermion/WilsonFermion.cc | 3 +- lib/qcd/action/fermion/WilsonFermion.h | 2 +- lib/qcd/action/fermion/WilsonFermion5D.cc | 227 +++++++++++++++++- lib/qcd/action/fermion/WilsonFermion5D.h | 5 +- tests/core/Test_fft.cc | 3 +- 8 files changed, 311 insertions(+), 18 deletions(-) diff --git a/lib/qcd/action/fermion/DomainWallFermion.h b/lib/qcd/action/fermion/DomainWallFermion.h index 72ce8f42..99f64865 100644 --- a/lib/qcd/action/fermion/DomainWallFermion.h +++ b/lib/qcd/action/fermion/DomainWallFermion.h @@ -8,6 +8,7 @@ Author: Peter Boyle Author: Peter Boyle +Author: Vera Guelpers 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 @@ -42,8 +43,58 @@ namespace Grid { INHERIT_IMPL_TYPES(Impl); public: - void MomentumSpacePropagator(FermionField &out,const FermionField &in,RealD _m) { - this->MomentumSpacePropagatorHt(out,in,_m); + void FreePropagator(const FermionField &in,FermionField &out,RealD mass, std::vector twist, bool fiveD) { + FermionField in_k(in._grid); + FermionField prop_k(in._grid); + + FFT theFFT((GridCartesian *) in._grid); + + //phase for boundary condition + ComplexField coor(in._grid); + ComplexField ph(in._grid); ph = zero; + FermionField in_buf(in._grid); in_buf = zero; + Complex ci(0.0,1.0); + assert(twist.size() == Nd);//check that twist 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]))); + } + in_buf = exp((RealD)(2.0*M_PI)*ci*ph*(-1.0))*in; + + if(fiveD){//FFT only on temporal and spatial dimensions + std::vector mask(Nd+1,1); mask[0] = 0; + theFFT.FFT_dim_mask(in_k,in_buf,mask,FFT::forward); + this->MomentumSpacePropagatorHt_5d(prop_k,in_k,mass,twist); + theFFT.FFT_dim_mask(out,prop_k,mask,FFT::backward); + } + else{ + theFFT.FFT_all_dim(in_k,in,FFT::forward); + this->MomentumSpacePropagatorHt(prop_k,in_k,mass,twist); + theFFT.FFT_all_dim(out,prop_k,FFT::backward); + } + + //phase for boundary condition + out = out * exp((RealD)(2.0*M_PI)*ci*ph); + }; + + virtual void FreePropagator(const FermionField &in,FermionField &out,RealD mass,std::vector twist) { + bool fiveD = true; //5d propagator by default + FreePropagator(in,out,mass,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); + }; + + virtual void FreePropagator(const FermionField &in,FermionField &out,RealD mass) { + bool fiveD = true; //5d propagator by default + std::vector twist(Nd,0.0); //default: periodic boundarys in all directions + FreePropagator(in,out,mass,twist,fiveD); }; virtual void Instantiatable(void) {}; diff --git a/lib/qcd/action/fermion/FermionOperator.h b/lib/qcd/action/fermion/FermionOperator.h index f81c9e2c..1ef99b85 100644 --- a/lib/qcd/action/fermion/FermionOperator.h +++ b/lib/qcd/action/fermion/FermionOperator.h @@ -9,6 +9,7 @@ Author: Peter Boyle Author: Peter Boyle Author: Peter Boyle +Author: Vera Guelpers 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 @@ -95,17 +96,38 @@ namespace Grid { virtual void Mdir (const FermionField &in, FermionField &out,int dir,int disp)=0; // case by case Wilson, Clover, Cayley, ContFrac, PartFrac - virtual void MomentumSpacePropagator(FermionField &out,const FermionField &in,RealD _m) { assert(0);}; + virtual void MomentumSpacePropagator(FermionField &out,const FermionField &in,RealD _m,std::vector twist) { assert(0);}; - virtual void FreePropagator(const FermionField &in,FermionField &out,RealD mass) { + virtual void FreePropagator(const FermionField &in,FermionField &out,RealD mass,std::vector twist) { FFT theFFT((GridCartesian *) in._grid); FermionField in_k(in._grid); FermionField prop_k(in._grid); - theFFT.FFT_all_dim(in_k,in,FFT::forward); - this->MomentumSpacePropagator(prop_k,in_k,mass); + //phase for boundary condition + ComplexField coor(in._grid); + ComplexField ph(in._grid); ph = zero; + FermionField in_buf(in._grid); in_buf = zero; + Complex ci(0.0,1.0); + assert(twist.size() == Nd);//check that twist is Nd + for(unsigned int nu = 0; nu < Nd; nu++) + { + LatticeCoordinate(coor, nu); + ph = ph + twist[nu]*coor*((1./(in._grid->_fdimensions[nu]))); + } + in_buf = exp((RealD)(2.0*M_PI)*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((RealD)(2.0*M_PI)*ci*ph); + + }; + virtual void FreePropagator(const FermionField &in,FermionField &out,RealD mass) { + std::vector twist(Nd,0.0); //default: periodic boundarys in all directions + FreePropagator(in,out,mass,twist); }; /////////////////////////////////////////////// diff --git a/lib/qcd/action/fermion/OverlapWilsonCayleyTanhFermion.h b/lib/qcd/action/fermion/OverlapWilsonCayleyTanhFermion.h index f516c5d0..fd7d74df 100644 --- a/lib/qcd/action/fermion/OverlapWilsonCayleyTanhFermion.h +++ b/lib/qcd/action/fermion/OverlapWilsonCayleyTanhFermion.h @@ -42,8 +42,8 @@ namespace Grid { INHERIT_IMPL_TYPES(Impl); public: - void MomentumSpacePropagator(FermionField &out,const FermionField &in,RealD _m) { - this->MomentumSpacePropagatorHw(out,in,_m); + void MomentumSpacePropagator(FermionField &out,const FermionField &in,RealD _m,std::vector twist) { + this->MomentumSpacePropagatorHw(out,in,_m,twist); }; // Constructors diff --git a/lib/qcd/action/fermion/WilsonFermion.cc b/lib/qcd/action/fermion/WilsonFermion.cc index 2d9cf22d..53e81c39 100644 --- a/lib/qcd/action/fermion/WilsonFermion.cc +++ b/lib/qcd/action/fermion/WilsonFermion.cc @@ -162,7 +162,7 @@ void WilsonFermion::MooeeInvDag(const FermionField &in, FermionField &out) MooeeInv(in,out); } template -void WilsonFermion::MomentumSpacePropagator(FermionField &out, const FermionField &in,RealD _m) +void WilsonFermion::MomentumSpacePropagator(FermionField &out, const FermionField &in,RealD _m,std::vector twist) { typedef typename FermionField::vector_type vector_type; typedef typename FermionField::scalar_type ScalComplex; @@ -195,6 +195,7 @@ void WilsonFermion::MomentumSpacePropagator(FermionField &out, const Fermi RealD TwoPiL = M_PI * 2.0/ latt_size[mu]; kmu = TwoPiL * kmu; + kmu = kmu + TwoPiL * one * twist[mu];//momentum for twisted boundary conditions wilson = wilson + 2.0*sin(kmu*0.5)*sin(kmu*0.5); // Wilson term diff --git a/lib/qcd/action/fermion/WilsonFermion.h b/lib/qcd/action/fermion/WilsonFermion.h index ea25ed7f..3985771c 100644 --- a/lib/qcd/action/fermion/WilsonFermion.h +++ b/lib/qcd/action/fermion/WilsonFermion.h @@ -96,7 +96,7 @@ class WilsonFermion : public WilsonKernels, public WilsonFermionStatic { virtual void MooeeInv(const FermionField &in, FermionField &out); virtual void MooeeInvDag(const FermionField &in, FermionField &out); - virtual void MomentumSpacePropagator(FermionField &out,const FermionField &in,RealD _mass) ; + virtual void MomentumSpacePropagator(FermionField &out,const FermionField &in,RealD _mass,std::vector twist) ; //////////////////////// // Derivative interface diff --git a/lib/qcd/action/fermion/WilsonFermion5D.cc b/lib/qcd/action/fermion/WilsonFermion5D.cc index 6f82aad2..6d4c6967 100644 --- a/lib/qcd/action/fermion/WilsonFermion5D.cc +++ b/lib/qcd/action/fermion/WilsonFermion5D.cc @@ -13,6 +13,7 @@ Author: Peter Boyle Author: paboyle Author: Guido Cossu Author: Andrew Lawson +Author: Vera Guelpers 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 @@ -563,7 +564,221 @@ void WilsonFermion5D::DW(const FermionField &in, FermionField &out,int dag } template -void WilsonFermion5D::MomentumSpacePropagatorHt(FermionField &out,const FermionField &in, RealD mass) +void WilsonFermion5D::MomentumSpacePropagatorHt_5d(FermionField &out,const FermionField &in, RealD mass,std::vector twist) +{ + // what type LatticeComplex + GridBase *_grid = _FourDimGrid; + GridBase *_5dgrid = _FiveDimGrid; + + conformable(_5dgrid,out._grid); + + FermionField PRsource(_5dgrid); + FermionField PLsource(_5dgrid); + FermionField buf1_4d(_grid); + FermionField buf2_4d(_grid); + FermionField GR(_5dgrid); + FermionField GL(_5dgrid); + FermionField bufL_4d(_grid); + FermionField bufR_4d(_grid); + + unsigned int Ls = in._grid->_rdimensions[0]; + + typedef typename FermionField::vector_type vector_type; + typedef typename FermionField::scalar_type ScalComplex; + typedef iSinglet Tcomplex; + typedef Lattice > LatComplex; + + Gamma::Algebra Gmu [] = { + Gamma::Algebra::GammaX, + Gamma::Algebra::GammaY, + Gamma::Algebra::GammaZ, + Gamma::Algebra::GammaT + }; + + Gamma g5(Gamma::Algebra::Gamma5); + + std::vector latt_size = _grid->_fdimensions; + + LatComplex sk(_grid); sk = zero; + LatComplex sk2(_grid); sk2= zero; + LatComplex W(_grid); W= zero; + LatComplex a(_grid); a= zero; + LatComplex one (_grid); one = ScalComplex(1.0,0.0); + LatComplex cosha(_grid); + LatComplex kmu(_grid); + LatComplex Wea(_grid); + LatComplex Wema(_grid); + LatComplex sinha(_grid); + LatComplex sinhaLs(_grid); + LatComplex coshaLs(_grid); + LatComplex A(_grid); + LatComplex F(_grid); + LatComplex App(_grid); + LatComplex Amm(_grid); + LatComplex Bpp(_grid); + LatComplex Bmm(_grid); + LatComplex ABpm(_grid); //Apm=Amp=Bpm=Bmp + LatComplex signW(_grid); + + ScalComplex ci(0.0,1.0); + + for(int mu=0;mu alpha + //////////////////////////////////////////// + cosha = (one + W*W + sk) / (abs(W)*2.0); + + // FIXME Need a Lattice acosh + for(int idx=0;idx<_grid->lSites();idx++){ + std::vector lcoor(Nd); + Tcomplex cc; + RealD sgn; + _grid->LocalIndexToLocalCoor(idx,lcoor); + peekLocalSite(cc,cosha,lcoor); + assert((double)real(cc)>=1.0); + assert(fabs((double)imag(cc))<=1.0e-15); + cc = ScalComplex(::acosh(real(cc)),0.0); + pokeLocalSite(cc,a,lcoor); + } + + Wea = ( exp( a) * abs(W) ); + Wema= ( exp(-a) * abs(W) ); + sinha = 0.5*(exp( a) - exp(-a)); + sinhaLs = 0.5*(exp( a*Ls) - exp(-a*Ls)); + coshaLs = 0.5*(exp( a*Ls) + exp(-a*Ls)); + + A = one / (abs(W) * sinha * 2.0) * one / (sinhaLs * 2.0); + F = exp( a*Ls) * (one - Wea + (Wema - one) * mass*mass); + F = F + exp(-a*Ls) * (Wema - one + (one - Wea) * mass*mass); + F = F - abs(W) * sinha * 4.0 * mass; + + Bpp = (A/F) * (exp(-a*Ls*2.0) - one) * (one - Wema) * (one - mass*mass * one); + Bmm = (A/F) * (one - exp(a*Ls*2.0)) * (one - Wea) * (one - mass*mass * one); + App = (A/F) * (exp(-a*Ls*2.0) - one) * exp(-a) * (exp(-a) - abs(W)) * (one - mass*mass * one); + Amm = (A/F) * (one - exp(a*Ls*2.0)) * exp(a) * (exp(a) - abs(W)) * (one - mass*mass * one); + ABpm = (A/F) * abs(W) * sinha * 2.0 * (one + mass * coshaLs * 2.0 + mass*mass * one); + + //P+ source, P- source + PRsource = (in + g5 * in) * 0.5; + PLsource = (in - g5 * in) * 0.5; + + //calculate GR, GL + for(unsigned int ss=1;ss<=Ls;ss++) + { + bufR_4d = zero; + bufL_4d = zero; + for(unsigned int tt=1;tt<=Ls;tt++) + { + //possible sign if W<0 + if((ss+tt)%2==1) signW = abs(W)/W; + else signW = one; + + unsigned int f = (ss > tt) ? ss-tt : tt-ss; //f = abs(ss-tt) + //GR + buf1_4d = zero; + ExtractSlice(buf1_4d, PRsource, (tt-1), 0); + //G(s,t) + bufR_4d = bufR_4d + A * exp(a*Ls) * exp(-a*f) * signW * buf1_4d + A * exp(-a*Ls) * exp(a*f) * signW * buf1_4d; + //A++*exp(a(s+t)) + bufR_4d = bufR_4d + App * exp(a*ss) * exp(a*tt) * signW * buf1_4d ; + //A+-*exp(a(s-t)) + bufR_4d = bufR_4d + ABpm * exp(a*ss) * exp(-a*tt) * signW * buf1_4d ; + //A-+*exp(a(-s+t)) + bufR_4d = bufR_4d + ABpm * exp(-a*ss) * exp(a*tt) * signW * buf1_4d ; + //A--*exp(a(-s-t)) + bufR_4d = bufR_4d + Amm * exp(-a*ss) * exp(-a*tt) * signW * buf1_4d ; + + //GL + buf2_4d = zero; + ExtractSlice(buf2_4d, PLsource, (tt-1), 0); + //G(s,t) + bufL_4d = bufL_4d + A * exp(a*Ls) * exp(-a*f) * signW * buf2_4d + A * exp(-a*Ls) * exp(a*f) * signW * buf2_4d; + //B++*exp(a(s+t)) + bufL_4d = bufL_4d + Bpp * exp(a*ss) * exp(a*tt) * signW * buf2_4d ; + //B+-*exp(a(s-t)) + bufL_4d = bufL_4d + ABpm * exp(a*ss) * exp(-a*tt) * signW * buf2_4d ; + //B-+*exp(a(-s+t)) + bufL_4d = bufL_4d + ABpm * exp(-a*ss) * exp(a*tt) * signW * buf2_4d ; + //B--*exp(a(-s-t)) + bufL_4d = bufL_4d + Bmm * exp(-a*ss) * exp(-a*tt) * signW * buf2_4d ; + } + InsertSlice(bufR_4d, GR, (ss-1), 0); + InsertSlice(bufL_4d, GL, (ss-1), 0); + } + +//calculate propagator + for(unsigned int ss=1;ss<=Ls;ss++) + { + bufR_4d = zero; + bufL_4d = zero; + + //(i*gamma_mu*sin(p_mu) - W)*(GL*P- source) + buf1_4d = zero; + ExtractSlice(buf1_4d, GL, (ss-1), 0); + buf2_4d = zero; + for(int mu=0;mu +void WilsonFermion5D::MomentumSpacePropagatorHt(FermionField &out,const FermionField &in, RealD mass,std::vector twist) { // what type LatticeComplex GridBase *_grid = _FourDimGrid; @@ -606,6 +821,7 @@ void WilsonFermion5D::MomentumSpacePropagatorHt(FermionField &out,const Fe RealD TwoPiL = M_PI * 2.0/ latt_size[mu]; kmu = TwoPiL * kmu; + kmu = kmu + TwoPiL * one * twist[mu];//momentum for twisted boundary conditions sk2 = sk2 + 2.0*sin(kmu*0.5)*sin(kmu*0.5); sk = sk + sin(kmu) *sin(kmu); @@ -619,7 +835,7 @@ void WilsonFermion5D::MomentumSpacePropagatorHt(FermionField &out,const Fe //////////////////////////////////////////// // Cosh alpha -> alpha //////////////////////////////////////////// - cosha = (one + W*W + sk) / (W*2.0); + cosha = (one + W*W + sk) / (abs(W)*2.0); // FIXME Need a Lattice acosh for(int idx=0;idx<_grid->lSites();idx++){ @@ -634,8 +850,8 @@ void WilsonFermion5D::MomentumSpacePropagatorHt(FermionField &out,const Fe pokeLocalSite(cc,a,lcoor); } - Wea = ( exp( a) * W ); - Wema= ( exp(-a) * W ); + Wea = ( exp( a) * abs(W) ); + Wema= ( exp(-a) * abs(W) ); num = num + ( one - Wema ) * mass * in; denom= ( Wea - one ) + mass*mass * (one - Wema); @@ -643,7 +859,7 @@ void WilsonFermion5D::MomentumSpacePropagatorHt(FermionField &out,const Fe } template -void WilsonFermion5D::MomentumSpacePropagatorHw(FermionField &out,const FermionField &in,RealD mass) +void WilsonFermion5D::MomentumSpacePropagatorHw(FermionField &out,const FermionField &in,RealD mass,std::vector twist) { Gamma::Algebra Gmu [] = { Gamma::Algebra::GammaX, @@ -683,6 +899,7 @@ void WilsonFermion5D::MomentumSpacePropagatorHw(FermionField &out,const Fe RealD TwoPiL = M_PI * 2.0/ latt_size[mu]; kmu = TwoPiL * kmu; + kmu = kmu + TwoPiL * one * twist[mu];//momentum for twisted boundary conditions sk2 = sk2 + 2.0*sin(kmu*0.5)*sin(kmu*0.5); sk = sk + sin(kmu)*sin(kmu); diff --git a/lib/qcd/action/fermion/WilsonFermion5D.h b/lib/qcd/action/fermion/WilsonFermion5D.h index 21da4c31..d22b5d6f 100644 --- a/lib/qcd/action/fermion/WilsonFermion5D.h +++ b/lib/qcd/action/fermion/WilsonFermion5D.h @@ -118,8 +118,9 @@ namespace QCD { virtual void DhopDerivEO(GaugeField &mat,const FermionField &U,const FermionField &V,int dag); virtual void DhopDerivOE(GaugeField &mat,const FermionField &U,const FermionField &V,int dag); - void MomentumSpacePropagatorHt(FermionField &out,const FermionField &in,RealD mass) ; - void MomentumSpacePropagatorHw(FermionField &out,const FermionField &in,RealD mass) ; + void MomentumSpacePropagatorHt_5d(FermionField &out,const FermionField &in,RealD mass,std::vector twist) ; + void MomentumSpacePropagatorHt(FermionField &out,const FermionField &in,RealD mass,std::vector twist) ; + void MomentumSpacePropagatorHw(FermionField &out,const FermionField &in,RealD mass,std::vector twist) ; // Implement hopping term non-hermitian hopping term; half cb or both // Implement s-diagonal DW diff --git a/tests/core/Test_fft.cc b/tests/core/Test_fft.cc index b2336cfa..b7f9f000 100644 --- a/tests/core/Test_fft.cc +++ b/tests/core/Test_fft.cc @@ -309,7 +309,8 @@ int main (int argc, char ** argv) // Momentum space prop std::cout << " Solving by FFT and Feynman rules" <