Grid physics library, www.github.com/paboyle/Grid 

Source file: lib/qcd/spin/gamma-gen/gamma-gen.nb

Copyright (C) 2015-2018

Author: Antonin Portelli <antonin.portelli@me.com>

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
the Free Software Foundation; either version 2 of the License, or
(at your option) any later version.

This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
GNU General Public License for more details.

You should have received a copy of the GNU General Public License along
with this program; if not, write to the Free Software Foundation, Inc.,
51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.

See the full license in the file “LICENSE” in the top level distribution directory

ClearAll["Global`*"]
Unprotect["Dot"];Unprotect["ConjugateTranspose"];
SetDirectory[NotebookDirectory[]]

Clifford algebra generation

com[a_,b_]:=a.b-b.a
acom[a_,b_]:=a.b+b.a

Definition of the matrix representation of the algebra:

mat[id]=IdentityMatrix[4];
mat[gx]=(({{0, 0, 0, I}, {0, 0, I, 0}, {0, -I, 0, 0}, {-I, 0, 0, 0}}));mat[gy]=(({{0, 0, 0, -1}, {0, 0, 1, 0}, {0, 1, 0, 0}, {-1, 0, 0, 0}}));mat[gz]=(({{0, 0, I, 0}, {0, 0, 0, -I}, {-I, 0, 0, 0}, {0, I, 0, 0}}));mat[gt]=(({{0, 0, 1, 0}, {0, 0, 0, 1}, {1, 0, 0, 0}, {0, 1, 0, 0}}));
mat[g5]=mat[gx].mat[gy].mat[gz].mat[gt];
mat[gx5]=mat[gx].mat[g5];
mat[gy5]=mat[gy].mat[g5];
mat[gz5]=mat[gz].mat[g5];
mat[gt5]=mat[gt].mat[g5];
mat[σxy]=com[mat[gx],mat[gy]]/2;
mat[σxz]=com[mat[gx],mat[gz]]/2;
mat[σxt]=com[mat[gx],mat[gt]]/2;
mat[σyz]=com[mat[gy],mat[gz]]/2;
mat[σyt]=com[mat[gy],mat[gt]]/2;
mat[σzt]=com[mat[gz],mat[gt]]/2;
palg={id,gx,gy,gz,gt,g5,gx5,gy5,gz5,gt5,σxy,σxz,σxt,σyz,σyt,σzt};
malg=Table[-i,{i,palg}];
alg=Union[palg,malg];
Do[mat[-i]=-mat[i],{i,palg}]

Simplify[Conjugate[{w0,w1,w2,w3}].((a0 mat[gx]+a1 mat[gy]+a2 mat[gz]+a3 mat[gt])/.{a0->(1/(2I))(b0-Conjugate[b0]),a1->(1/2)(b0+Conjugate[b0]),a2->(1/(2I))(b1-Conjugate[b1]),a3->(1/2)(b1+Conjugate[b1])}).{v0,v1,v2,v3}//MatrixForm]

Simplify[((a0 mat[gx]+a1 mat[gy]+a2 mat[gz]+a3 mat[gt])/.{a0->(1/(2I))(b0-Conjugate[b0]),a1->(1/2)(b0+Conjugate[b0]),a2->(1/(2I))(b1-Conjugate[b1]),a3->(1/2)(b1+Conjugate[b1])})//MatrixForm]

Generation of the abstract algebra:

Do[Dot[i1,i2]=Select[alg,mat[i1].mat[i2]==mat[#]&][[1]],{i1,palg},{i2,palg}]
Dot[-i1_,i2_]:=-i1.i2;Dot[i1_,-i2_]:=-i1.i2;Dot[-i1_,-i2_]:=i1.i2;
Do[ConjugateTranspose[i]=Select[alg,ConjugateTranspose[mat[i]]==mat[#]&][[1]],{i,palg}];
ConjugateTranspose[-i_]:=-ConjugateTranspose[i]

Check that we can reconstruct the Euclidean metric:

Table[acom[i1,i2]/2,{i1,{gx,gy,gz,gt}},{i2,{gx,gy,gz,gt}}]//MatrixForm

Full multiplication table:

Table[i1.i2,{i1,palg},{i2,palg}]//MatrixForm

Header file Gamma.h generation

File skeleton:

makeHeaderCode[enumCode_,funcCode_]:=Module[{out=""},
out = out<>
"#ifndef GRID_QCD_GAMMA_H
#define GRID_QCD_GAMMA_H

// code generated by the Mathematica notebook gamma-gen/gamma-gen.nb

#include <array>

namespace Grid {
namespace QCD {

class Gamma {
  public:
";
out = out<>enumCode;
out=out<>"    static const std::array<const char *, nGamma>                name;
    static const std::array<std::array<Algebra, nGamma>, nGamma> mul;
    static const std::array<Algebra, nGamma>                     adj;
    static const std::array<const Gamma, 4>                      gmu;
    static const std::array<const Gamma, 16>                     gall;
    Algebra                                                      g;
  public:
    Gamma(Algebra initg): g(initg) {}  
};

";
out =out<>funcCode;
out=out<>"
}}

#endif // GRID_QCD_GAMMA_H
";
Return[out]]

Gamma enum generation:

enum[id]="Algebra::Identity";
enum[gx]="Algebra::GammaX";
enum[gy]="Algebra::GammaY";
enum[gz]="Algebra::GammaZ";
enum[gt]="Algebra::GammaT";
enum[g5]="Algebra::Gamma5";
enum[gx5]="Algebra::GammaXGamma5";
enum[gy5]="Algebra::GammaYGamma5";
enum[gz5]="Algebra::GammaZGamma5";
enum[gt5]="Algebra::GammaTGamma5";
enum[σxy]="Algebra::SigmaXY";
enum[σxz]="Algebra::SigmaXZ";
enum[σxt]="Algebra::SigmaXT";
enum[σyz]="Algebra::SigmaYZ";
enum[σyt]="Algebra::SigmaYT";
enum[σzt]="Algebra::SigmaZT";
enum[-id]="Algebra::MinusIdentity";
enum[-gx]="Algebra::MinusGammaX";
enum[-gy]="Algebra::MinusGammaY";
enum[-gz]="Algebra::MinusGammaZ";
enum[-gt]="Algebra::MinusGammaT";
enum[-g5]="Algebra::MinusGamma5";
enum[-gx5]="Algebra::MinusGammaXGamma5";
enum[-gy5]="Algebra::MinusGammaYGamma5";
enum[-gz5]="Algebra::MinusGammaZGamma5";
enum[-gt5]="Algebra::MinusGammaTGamma5";
enum[-σxy]="Algebra::MinusSigmaXY";
enum[-σxz]="Algebra::MinusSigmaXZ";
enum[-σxt]="Algebra::MinusSigmaXT";
enum[-σyz]="Algebra::MinusSigmaYZ";
enum[-σyt]="Algebra::MinusSigmaYT";
enum[-σzt]="Algebra::MinusSigmaZT";

makeEnumCode[indent_]:=Module[{out},
out=indent<>"GRID_SERIALIZABLE_ENUM(Algebra, undef,";
ind=0;
Do[
out=out<>"\n"<>indent<> "                       "<>StringPadRight[StringTrim[enum[i],"Algebra::"],17]<>", "<>ToString[ind];
out = out<>If[i==Last[alg],");\n"<>indent<>"static constexpr unsigned int nGamma = "<>ToString[ind+1]<>";\n",",",","];
ind=ind+1,
{i,alg}];
Return[out]]

Multiplication functions generation:

rhsv={rhs[0],rhs[1],rhs[2],rhs[3]};
vecCoef[v_[i_]]:=ToString[v]<>"("<>ToString[i]<>")"
lmatCoef[v_[i_]]:=ToString[v]<>"("<>ToString[i]<>", i)"
rmatCoef[v_[i_]]:=ToString[v]<>"(i, "<>ToString[i]<>")"
expr[ir_,rhs[i_],strFn_]:=strFn[ret[ir]]<>" = "<> strFn[rhs[i]]<>";"
expr[ir_,Times[-1,rhs[i_]],strFn_]:=strFn[ret[ir]]<>" = -"<> strFn[rhs[i]]<>";"
expr[ir_,Times[Complex [0,1],rhs[i_]],strFn_]:=strFn[ret[ir]]<>" = timesI("<> strFn[rhs[i]]<>");"
expr[ir_,Times[Complex [0,-1],rhs[i_]],strFn_]:=strFn[ret[ir]]<>" = timesMinusI("<> strFn[rhs[i]]<>");"
vecFunc[g_]:=Module[{out="",ir=0},
out="template<class vtype>
inline void mult"<>
StringTrim[enum[g],"Algebra::"]<>
"(iVector<vtype, Ns> &ret, const iVector<vtype, Ns> &rhs)\n{\n";
Do[out=out<>"  "<>expr[ir,c,vecCoef]<>"\n";ir=ir+1,{c,mat[g].rhsv}];
out =out<>"};\n";
Return[out]];
lmatFunc[g_]:=Module[{out="",ir=0},
out="template<class vtype>
inline void lmult"<>
StringTrim[enum[g],"Algebra::"]<>
"(iMatrix<vtype, Ns> &ret, const iMatrix<vtype, Ns> &rhs)\n{\n"<>
"  for(unsigned int i = 0; i < Ns; ++i)\n  {\n";
Do[out=out<>"    "<>expr[ir,c,lmatCoef]<>"\n";ir=ir+1,{c,mat[g].rhsv}];
out =out<>"  }\n};\n";
Return[out]]
rmatFunc[g_]:=Module[{out="",ir=0},
out="template<class vtype>
inline void rmult"<>
StringTrim[enum[g],"Algebra::"]<>
"(iMatrix<vtype, Ns> &ret, const iMatrix<vtype, Ns> &rhs)\n{\n"<>
"  for(unsigned int i = 0; i < Ns; ++i)\n  {\n";
Do[out=out<>"    "<>expr[ir,c,rmatCoef]<>"\n";ir=ir+1,{c,rhsv.mat[g]}];
out =out<>"  }\n};\n";
Return[out]]
allFunc[g_]:=Module[{out=""},
out=out<>vecFunc[g]<>"\n";
out=out<>lmatFunc[g]<>"\n";
out=out<>rmatFunc[g]<>"\n";
Return[out]]
vecFinal[]:=Module[{out=""},
out="template<class vtype> 
inline auto operator*(const Gamma &G, const iVector<vtype, Ns> &arg)
->typename std::enable_if<matchGridTensorIndex<iVector<vtype, Ns>, SpinorIndex>::value, iVector<vtype, Ns>>::type
{
  iVector<vtype, Ns> ret;

  switch (G.g) 
  {
";
Do[out =out<>"  case Gamma::"<>enum[g]<>":\n";
out=out<>"    mult"<>StringTrim[enum[g],"Algebra::"]<>"(ret, arg); break;\n",
{g,alg}];
out=out<>"  default: assert(0);
  }
  
  return ret;
}\n";
Return[out]]
lmatFinal[]:=Module[{out=""},
out="template<class vtype> 
inline auto operator*(const Gamma &G, const iMatrix<vtype, Ns> &arg)
->typename std::enable_if<matchGridTensorIndex<iMatrix<vtype, Ns>, SpinorIndex>::value, iMatrix<vtype, Ns>>::type
{
  iMatrix<vtype, Ns> ret;

  switch (G.g) 
  {
";
Do[out =out<>"  case Gamma::"<>enum[g]<>":\n";
out=out<>"    lmult"<>StringTrim[enum[g],"Algebra::"]<>"(ret, arg); break;\n",
{g,alg}];
out=out<>"  default: assert(0);
  }
  
  return ret;
}\n";
Return[out]]
makeFuncCode[]:=Module[{out=""},Do[out=out<>allFunc[g],{g,alg}];Return[out]]
rmatFinal[]:=Module[{out=""},
out="template<class vtype> 
inline auto operator*(const iMatrix<vtype, Ns> &arg, const Gamma &G)
->typename std::enable_if<matchGridTensorIndex<iMatrix<vtype, Ns>, SpinorIndex>::value, iMatrix<vtype, Ns>>::type
{
  iMatrix<vtype, Ns> ret;

  switch (G.g) 
  {
";
Do[out =out<>"  case Gamma::"<>enum[g]<>":\n";
out=out<>"    rmult"<>StringTrim[enum[g],"Algebra::"]<>"(ret, arg); break;\n",
{g,alg}];
out=out<>"  default: assert(0);
  }
  
  return ret;
}\n";
Return[out]]
makeFuncCode[]:=Module[{out=""},
Do[out=out<>allFunc[g],{g,alg}];
out=out<>vecFinal[]<>"\n"<>lmatFinal[]<>"\n"<>rmatFinal[];
Return[out]]

Header file generation:

file=OpenWrite["Gamma.h"];WriteString[file,makeHeaderCode[makeEnumCode["    "],makeFuncCode[]]];Close[file];

Gamma.cc source file generation

makeSourceCode[]:=Module[{out=""},
out = out<>
"// code generated by the Mathematica notebook gamma-gen/gamma-gen.nb

#include <Grid/Grid.h>

namespace Grid {
namespace QCD {";
out = out <>"

const std::array<const Gamma, 4> Gamma::gmu = {{
  Gamma(Gamma::Algebra::GammaX),
  Gamma(Gamma::Algebra::GammaY),
  Gamma(Gamma::Algebra::GammaZ),
  Gamma(Gamma::Algebra::GammaT)}};

const std::array<const Gamma, 16> Gamma::gall = {{
  Gamma(Gamma::Algebra::Identity),
  Gamma(Gamma::Algebra::Gamma5),
  Gamma(Gamma::Algebra::GammaX),
  Gamma(Gamma::Algebra::GammaY),
  Gamma(Gamma::Algebra::GammaZ),
  Gamma(Gamma::Algebra::GammaT),
  Gamma(Gamma::Algebra::GammaXGamma5),
  Gamma(Gamma::Algebra::GammaYGamma5),
  Gamma(Gamma::Algebra::GammaZGamma5),
  Gamma(Gamma::Algebra::GammaTGamma5),
  Gamma(Gamma::Algebra::SigmaXT),      
  Gamma(Gamma::Algebra::SigmaXY),      
  Gamma(Gamma::Algebra::SigmaXZ),      
  Gamma(Gamma::Algebra::SigmaYT),
  Gamma(Gamma::Algebra::SigmaYZ),
  Gamma(Gamma::Algebra::SigmaZT)}};

const std::array<const char *, Gamma::nGamma> Gamma::name = {{
";
Do[
out = out<>"  \""<>
StringPadRight[StringTrim[StringReplace[enum[g],"Minus"->"-"],"Algebra::"],13]<>
"\""<>
If[g==Last[alg],"}};\n\n",",\n",",\n"],
{g,alg}];
out = out<>"const std::array<Gamma::Algebra, Gamma::nGamma> Gamma::adj = {{\n";
Do[
out=out<>"  Gamma::"<>enum[g\[ConjugateTranspose]]<>If[g==Last[alg],"}};\n\n",",\n",",\n"],
{g,alg}];
out = out <>"const std::array<std::array<Gamma::Algebra, Gamma::nGamma>, Gamma::nGamma> Gamma::mul = {{\n";
Do[
Do[
out=out<>If[g2==First[alg],"  {{","   ","   "]
<>"Gamma::"<>enum[g1.g2]<>If[g2==Last[alg],"}}","",""]
<>If[g2==Last[alg],If[g1==Last[alg],"","\n  ,\n","\n  ,\n"],",\n",",\n"],
{g2,alg}],
{g1,alg}];
out=out<>"\n}};\n\n";
out =out<>"}}\n";
Return[out]]

file=OpenWrite["Gamma.cc"];WriteString[file,makeSourceCode[]];Close[file];