#ifndef __POISSON_H__
#define __POISSON_H__
#include "chromabase.h"
#include "tower_array.h"
namespace Chroma {
struct Poisson {
ComplexD sst;
ComplexD tst;
ComplexD ttsst;
ComplexD sttst;
ComplexD stsst;
ComplexD tttst;
Poisson()
{
sst = tst = ttsst = sttst = stsst = tttst = zero;
}
Poisson(const TowerArray<LatticeColorMatrix>& F,
const TowerArray<LatticeColorMatrix>& G,
const multi1d<LatticeColorMatrix>& P) {
computePoisson(F,G,P);
}
void computePoisson(const TowerArray<LatticeColorMatrix>& F,
const TowerArray<LatticeColorMatrix>& G,
const multi1d<LatticeColorMatrix>& P)
{
sst = tst = ttsst = sttst = stsst = tttst = zero;
for(int mu=0; mu < Nd; mu++) {
// {S,{S,T}} = tr(F1^2)
sst -= Real(2)*sum(trace( F[mu][0]*F[mu][0] ));
// {T,{S,T}} = - tr(F2 P)
tst += Real(2)*sum(trace( F[mu][1]*P[mu] ));
// {T,{T,{S,{S,T}}}} = 2( tr(F1 F3} + tr(F2^2) )
ttsst -= Real(4)*sum(trace( F[mu][0]*F[mu][2])); // F1 F3 term
ttsst -= Real(4)*sum(trace( F[mu][1]*F[mu][1])); // F2^2 term
// {{S,T},{T,{S,T}}} = -3 tr ( [F1,F2] P )
// - tr ( [F1,P] [F1,P] )
// + tr ( F1,F3 )
// -2 tr ( F2,F2 )
// Commutator 1 = [F1,F2] = F1 F2 - F2 F1
LatticeColorMatrix commut1 = F[mu][0]*F[mu][1];
commut1 -= F[mu][1]*F[mu][0];
// Commutator 2 = {F1, P] = F1 P - P F1
LatticeColorMatrix commut2 = F[mu][0]*P[mu];
commut2 -= P[mu]*F[mu][0];
sttst += Real(6)*sum(trace(commut1*P[mu])); // -3 tr ( [F1,F2] P )
sttst += Real(2)*sum(trace(commut2*commut2)); // - tr ( [F1,P] [F1,P] )
sttst -= Real(2)*sum(trace(F[mu][0]*F[mu][2])); // + tr ( F1,F3 )
sttst += Real(4)*sum(trace(F[mu][1]*F[mu][1])); // -2 tr ( F2,F2 )
// {{S,T}, {S,{S,T}}} = -2 tr (F1 G1)
stsst += Real(4)*sum(trace(F[mu][0]*G[mu][1]));
// {T,{T,{T,{S,T}}}} = -tr( F4 P )
tttst += Real(2)*sum(trace(F[mu][3]*P[mu]));
// {S,{S,{S,{S,T}}}} = 0 so don't compute
}
}
void print(void) {
QDPIO::cout << "PoissonBrackets: " << endl;
QDPIO::cout << "=================" << endl;
QDPIO::cout << " {S,{S,T}} = " << sst << endl;
QDPIO::cout << " {T,{S,T}} = " << tst << endl;
QDPIO::cout << " {T,{T,{S,{S,T}}}} = " << ttsst << endl;
QDPIO::cout << " {{S,T},{T,{S,T}}} = " << sttst << endl;
QDPIO::cout << " {{S,T},{S,{S,T}}} = " << stsst << endl;
QDPIO::cout << " {T,{T,{T,{S,T}}}} = " << tttst << endl;
QDPIO::cout << "=================" << endl;
}
};
}
#endif