1
0
mirror of https://github.com/paboyle/Grid.git synced 2024-11-09 23:45:36 +00:00

Added two spinor functionality required to support the Wilson hopping term.

This commit is contained in:
Peter Boyle 2015-04-25 12:54:06 +01:00
parent 8b4073d84c
commit c5fa18eb20
5 changed files with 1227 additions and 59 deletions

22
TODO
View File

@ -13,7 +13,26 @@
* Consider switch std::vector to boost arrays or something lighter weight
boost::multi_array<type, 3> A()... to replace multi1d, multi2d etc..
* How to define simple matrix operations, such as flavour matrices?
* TensorRemove is a hack, come up with a long term rationalised approach to Complex vs. Scalar<Scalar<Scalar<Complex > > >
QDP forces use of "toDouble" to get back to non tensor scalar. This role is presently taken TensorRemove, but I
want to introduce a syntax that does not require this.
* norm2l is a hack. figure out syntax error and make this norm2. c.f. tests/Grid_gamma.cc
* std::vector replacement;
Had to change "reserve" back to "resize" on std::vector in Lattice class.
This forces the constructor call on EVERY element of the array with negative
performance effects on temporaries.
The reversion was required because copy constructur has to work.
CONCLUSION: I must implement a similar to vector without construction/fill on
resize. Find out if valarray or alternative works differently prior to
doing this since there may still be something I can use..
* Flavour matrices?
* Make the Tensor types and Complex etc... play more nicely.
@ -36,7 +55,6 @@
* Conformable test in Cshift routines.
* QDP++ regression suite and comparative benchmark
AUDITS:

View File

@ -6,7 +6,7 @@ namespace Grid {
// innerProduct Vector x Vector -> Scalar
// innerProduct Matrix x Matrix -> Scalar
///////////////////////////////////////////////////////////////////////////////////////
template<class sobj> inline RealD norm2l(sobj &arg){
template<class sobj> inline RealD norm2l(const sobj &arg){
typedef typename sobj::scalar_type scalar;
decltype(innerProduct(arg,arg)) nrm;
nrm = innerProduct(arg,arg);

View File

@ -7,6 +7,7 @@ namespace QCD {
static const int Nc=3;
static const int Ns=4;
static const int Nd=4;
static const int Nhs=2; // half spinor
static const int CbRed =0;
static const int CbBlack=1;
@ -38,16 +39,21 @@ namespace QCD {
template<typename vtype> using iColourVector = iScalar<iScalar<iVector<vtype, Nc> > >;
template<typename vtype> using iSpinColourVector = iScalar<iVector<iVector<vtype, Nc>, Ns> >;
typedef iSpinMatrix<Complex > SpinMatrix;
typedef iColourMatrix<Complex > ColourMatrix;
typedef iSpinColourMatrix<Complex > SpinColourMatrix;
typedef iLorentzColourMatrix<Complex > LorentzColourMatrix;
template<typename vtype> using iHalfSpinVector = iScalar<iVector<iScalar<vtype>, Nhs> >;
template<typename vtype> using iHalfSpinColourVector = iScalar<iVector<iVector<vtype, Nc>, Nhs> >;
typedef iSpinMatrix<Complex > SpinMatrix;
typedef iColourMatrix<Complex > ColourMatrix;
typedef iSpinColourMatrix<Complex > SpinColourMatrix;
typedef iLorentzColourMatrix<Complex > LorentzColourMatrix;
typedef iLorentzColourMatrix<ComplexF > LorentzColourMatrixF;
typedef iLorentzColourMatrix<ComplexD > LorentzColourMatrixD;
typedef iSpinVector<Complex > SpinVector;
typedef iColourVector<Complex > ColourVector;
typedef iSpinColourVector<Complex > SpinColourVector;
typedef iSpinVector<Complex > SpinVector;
typedef iColourVector<Complex > ColourVector;
typedef iSpinColourVector<Complex > SpinColourVector;
typedef iHalfSpinVector<Complex > HalfSpinVector;
typedef iHalfSpinColourVector<Complex > HalfSpinColourVector;
typedef iSpinMatrix<vComplex > vSpinMatrix;
@ -55,36 +61,46 @@ namespace QCD {
typedef iSpinColourMatrix<vComplex > vSpinColourMatrix;
typedef iLorentzColourMatrix<vComplex > vLorentzColourMatrix;
typedef iSpinVector<vComplex > vSpinVector;
typedef iColourVector<vComplex > vColourVector;
typedef iSpinColourVector<vComplex > vSpinColourVector;
typedef iSpinVector<vComplex > vSpinVector;
typedef iColourVector<vComplex > vColourVector;
typedef iSpinColourVector<vComplex > vSpinColourVector;
typedef iHalfSpinVector<vComplex > vHalfSpinVector;
typedef iHalfSpinColourVector<vComplex > vHalfSpinColourVector;
typedef iSinglet<Complex > TComplex; // This is painful. Tensor singlet complex type.
typedef iSinglet<vComplex > vTComplex; // what if we don't know the tensor structure
typedef iSinglet<Real > TReal; // Shouldn't need these; can I make it work without?
typedef iSinglet<vReal > vTReal;
typedef iSinglet<vInteger > vTInteger;
typedef iSinglet<Integer > TInteger;
typedef iSinglet<Complex > TComplex; // FIXME This is painful. Tensor singlet complex type.
typedef iSinglet<vComplex > vTComplex; // what if we don't know the tensor structure
typedef iSinglet<Real > TReal; // Shouldn't need these; can I make it work without?
typedef iSinglet<vReal > vTReal;
typedef iSinglet<vInteger > vTInteger;
typedef iSinglet<Integer > TInteger;
typedef Lattice<vTReal> LatticeReal;
typedef Lattice<vTComplex> LatticeComplex;
typedef Lattice<vTInteger> LatticeInteger; // Predicates for "where"
typedef Lattice<vTReal> LatticeReal;
typedef Lattice<vTComplex> LatticeComplex;
typedef Lattice<vTInteger> LatticeInteger; // Predicates for "where"
typedef Lattice<vColourMatrix> LatticeColourMatrix;
typedef Lattice<vSpinMatrix> LatticeSpinMatrix;
typedef Lattice<vSpinColourMatrix> LatticeSpinColourMatrix;
typedef Lattice<vSpinColourVector> LatticeSpinColourVector;
typedef Lattice<vSpinVector> LatticeSpinVector;
typedef Lattice<vColourVector> LatticeColourVector;
typedef Lattice<vSpinVector> LatticeSpinVector;
typedef Lattice<vColourVector> LatticeColourVector;
typedef Lattice<vSpinColourVector> LatticeSpinColourVector;
typedef Lattice<vHalfSpinVector> LatticeHalfSpinVector;
typedef Lattice<vHalfSpinColourVector> LatticeHalfSpinColourVector;
///////////////////////////////////////////
// Physical names for things
///////////////////////////////////////////
typedef Lattice<vSpinColourVector> LatticeFermion;
typedef Lattice<vSpinColourMatrix> LatticePropagator;
typedef Lattice<vLorentzColourMatrix> LatticeGaugeField;
typedef Lattice<vHalfSpinColourVector> LatticeHalfFermion;
typedef Lattice<vSpinColourVector> LatticeFermion;
typedef Lattice<vSpinColourMatrix> LatticePropagator;
typedef Lattice<vLorentzColourMatrix> LatticeGaugeField;
// Uhgg... typing this hurt ;)
// (my keyboard got burning hot when I typed this, must be the anti-Fermion)
typedef Lattice<vColourVector> LatticeStaggeredFermion;
typedef Lattice<vColourMatrix> LatticeStaggeredPropagator;
//////////////////////////////////////////////////////////////////////////////
// Peek and Poke named after physics attributes
@ -145,5 +161,7 @@ namespace QCD {
} // Grid
#include <qcd/Grid_qcd_dirac.h>
#include <qcd/Grid_qcd_2spinor.h>
//#include <qcd/Grid_qcd_pauli.h>
#endif

1082
lib/qcd/Grid_qcd_2spinor.h Normal file

File diff suppressed because it is too large Load Diff

View File

@ -29,42 +29,16 @@ int main (int argc, char ** argv)
SpinMatrix rr=zero;
SpinMatrix result;
SpinVector lv=zero;
SpinVector rv=zero;
SpinVector lv; random(sRNG,lv);
SpinVector rv; random(sRNG,rv);
for(int a=0;a<Ns;a++){
ident()(a,a) = 1.0;
}
Gamma::GammaMatrix g [] = {
Gamma::Identity,
Gamma::GammaX,
Gamma::GammaY,
Gamma::GammaZ,
Gamma::GammaT,
Gamma::Gamma5,
Gamma::MinusIdentity,
Gamma::MinusGammaX,
Gamma::MinusGammaY,
Gamma::MinusGammaZ,
Gamma::MinusGammaT,
Gamma::MinusGamma5
};
const char *list[] = {
"Identity ",
"GammaX ",
"GammaY ",
"GammaZ ",
"GammaT ",
"Gamma5 ",
"-Identity",
"-GammaX ",
"-GammaY ",
"-GammaZ ",
"-GammaT ",
"-Gamma5 ",
" "
};
const Gamma::GammaMatrix *g = Gamma::GammaMatrices;
const char **list = Gamma::GammaMatrixNames;
result =ll*Gamma(g[0])*rr;
result =ll*Gamma(g[0]);
rv = Gamma(g[0])*lv;
@ -124,5 +98,81 @@ int main (int argc, char ** argv)
std::cout << list[mu]<<" " << mag<<std::endl;
}
// Testing spins and reconstructs
SpinVector recon; random(sRNG,rv);
SpinVector full;
HalfSpinVector hsp,hsm;
// Xp
double mag;
spProjXp(hsm,rv);
spReconXp(recon,hsm);
full = rv + Gamma(Gamma::GammaX) *rv;
mag = TensorRemove(norm2l(full-recon));
std::cout << "Xp "<< mag<<std::endl;
// Xm
spProjXm(hsm,rv);
spReconXm(recon,hsm);
full = rv - Gamma(Gamma::GammaX) *rv;
mag = TensorRemove(norm2l(full-recon));
std::cout << "Xm "<< mag<<std::endl;
// Yp
spProjYp(hsm,rv);
spReconYp(recon,hsm);
full = rv + Gamma(Gamma::GammaY) *rv;
mag = TensorRemove(norm2l(full-recon));
std::cout << "Yp "<< mag<<std::endl;
// Ym
spProjYm(hsm,rv);
spReconYm(recon,hsm);
full = rv - Gamma(Gamma::GammaY) *rv;
mag = TensorRemove(norm2l(full-recon));
std::cout << "Ym "<< mag<<std::endl;
// Zp
spProjZp(hsm,rv);
spReconZp(recon,hsm);
full = rv + Gamma(Gamma::GammaZ) *rv;
mag = TensorRemove(norm2l(full-recon));
std::cout << "Zp "<< mag<<std::endl;
// Zm
spProjZm(hsm,rv);
spReconZm(recon,hsm);
full = rv - Gamma(Gamma::GammaZ) *rv;
mag = TensorRemove(norm2l(full-recon));
std::cout << "Zm "<< mag<<std::endl;
// Tp
spProjTp(hsm,rv);
spReconTp(recon,hsm);
full = rv + Gamma(Gamma::GammaT) *rv;
mag = TensorRemove(norm2l(full-recon));
std::cout << "Tp "<< mag<<std::endl;
// Tm
spProjTm(hsm,rv);
spReconTm(recon,hsm);
full = rv - Gamma(Gamma::GammaT) *rv;
mag = TensorRemove(norm2l(full-recon));
std::cout << "Tm "<< mag<<std::endl;
// 5p
spProj5p(hsm,rv);
spRecon5p(recon,hsm);
full = rv + Gamma(Gamma::Gamma5) *rv;
mag = TensorRemove(norm2l(full-recon));
std::cout << "5p "<< mag<<std::endl;
// 5m
spProj5m(hsm,rv);
spRecon5m(recon,hsm);
full = rv - Gamma(Gamma::Gamma5) *rv;
mag = TensorRemove(norm2l(full-recon));
std::cout << "5m "<< mag<<std::endl;
Grid_finalize();
}