portelli/Grid
1
0
Fork 0
mirror of https://github.com/paboyle/Grid.git synced 2025-06-13 20:57:06 +01:00
Code Releases Activity

Small change in the HMC interface.

Browse Source 
Example of multiple levels in the WilsonFermion hmc test.

Merge remote-tracking branch 'upstream/master'

Conflicts:
	lib/qcd/hmc/HMC.h
	lib/qcd/hmc/integrators/Integrator.h
	lib/qcd/hmc/integrators/Integrator_algorithm.h
	tests/Test_simd.cc
This commit is contained in:
neo
2015-07-30 17:16:04 +09:00
parent 479912a5ed 0b603225d1
commit bcdc67b152
100 changed files with 16360 additions and 1284 deletions
Download Patch File Download Diff File Expand all files Collapse all files

24
lib/qcd/action/ActionBase.h
View File

@ -7,14 +7,30 @@ template<class GaugeField>
class Action {
public:
virtual void init(const GaugeField &U, GridParallelRNG& pRNG) = 0;
virtual RealD S(const GaugeField &U) = 0; // evaluate the action
virtual void deriv(const GaugeField &U,GaugeField & dSdU ) = 0; // evaluate the action derivative
//virtual void refresh(const GaugeField & ) {} ;
virtual void init (const GaugeField &U, GridParallelRNG& pRNG) = 0; //
virtual RealD S (const GaugeField &U) = 0; // evaluate the action
virtual void deriv(const GaugeField &U,GaugeField & dSdU ) = 0; // evaluate the action derivative
virtual void refresh(const GaugeField & ) {}; // Default to no-op for actions with no internal fields
// Boundary conditions?
// Heatbath?
virtual ~Action() {};
};
// Could derive PseudoFermion action with a PF field, FermionField, and a Grid; implement refresh
template<class GaugeField, class FermionField>
class PseudoFermionAction : public Action<GaugeField> {
public:
FermionField Phi;
GridParallelRNG &pRNG;
GridBase &Grid;
PseudoFermionAction(GridBase &_Grid,GridParallelRNG &_pRNG) : Grid(_Grid), Phi(&_Grid), pRNG(_pRNG) {
};
virtual void refresh(const GaugeField &gauge) {
gaussian(Phi,pRNG);
};
};
}}
#endif

6
lib/qcd/action/Actions.h
View File

@ -79,4 +79,10 @@
///////////////////////////////////////////////////////////////////////////////
#include <qcd/action/fermion/g5HermitianLinop.h>
////////////////////////////////////////
// Pseudo fermion combinations
////////////////////////////////////////
#include <qcd/action/pseudofermion/TwoFlavour.h>
#endif

14
lib/qcd/action/fermion/ContinuedFractionFermion5D.cc
View File

@ -10,12 +10,12 @@ namespace Grid {
void ContinuedFractionFermion5D::SetCoefficientsZolotarev(RealD zolo_hi,Approx::zolotarev_data *zdata)
{
// How to check Ls matches??
// std::cout << Ls << " Ls"<<std::endl;
// std::cout << zdata->n << " - n"<<std::endl;
// std::cout << zdata->da << " -da "<<std::endl;
// std::cout << zdata->db << " -db"<<std::endl;
// std::cout << zdata->dn << " -dn"<<std::endl;
// std::cout << zdata->dd << " -dd"<<std::endl;
// std::cout<<GridLogMessage << Ls << " Ls"<<std::endl;
// std::cout<<GridLogMessage << zdata->n << " - n"<<std::endl;
// std::cout<<GridLogMessage << zdata->da << " -da "<<std::endl;
// std::cout<<GridLogMessage << zdata->db << " -db"<<std::endl;
// std::cout<<GridLogMessage << zdata->dn << " -dn"<<std::endl;
// std::cout<<GridLogMessage << zdata->dd << " -dd"<<std::endl;
assert(zdata->db==Ls);// Beta has Ls coeffs
@ -55,7 +55,7 @@ namespace Grid {
See[s] = Aee[s] - 1.0/See[s-1];
}
for(int s=0;s<Ls;s++){
std::cout <<"s = "<<s<<" Beta "<<Beta[s]<<" Aee "<<Aee[s] <<" See "<<See[s] <<std::endl;
std::cout<<GridLogMessage <<"s = "<<s<<" Beta "<<Beta[s]<<" Aee "<<Aee[s] <<" See "<<See[s] <<std::endl;
}
}

2
lib/qcd/action/fermion/DomainWallFermion.h
View File

@ -32,7 +32,7 @@ namespace Grid {
Approx::zolotarev_data *zdata = Approx::higham(eps,this->Ls);// eps is ignored for higham
assert(zdata->n==this->Ls);
std::cout << "DomainWallFermion with Ls="<<Ls<<std::endl;
std::cout<<GridLogMessage << "DomainWallFermion with Ls="<<Ls<<std::endl;
// Call base setter
this->CayleyFermion5D::SetCoefficientsTanh(zdata,1.0,0.0);

23
lib/qcd/action/fermion/FermionOperator.h
View File

@ -39,10 +39,27 @@ namespace Grid {
virtual void Dhop (const FermionField &in, FermionField &out,int dag)=0;
virtual void DhopOE(const FermionField &in, FermionField &out,int dag)=0;
virtual void DhopEO(const FermionField &in, FermionField &out,int dag)=0;
virtual void DhopDir(const FermionField &in, FermionField &out,int dir,int disp)=0; // implemented by WilsonFermion and WilsonFermion5D
virtual void Mdiag(const FermionField &in, FermionField &out) { Mooee(in,out);}; // Same as Mooee applied to both CB's
virtual void Mdir (const FermionField &in, FermionField &out,int dir,int disp)=0; // case by case Wilson, Clover, Cayley, ContFrac, PartFrac
virtual void DhopDir(const FermionField &in, FermionField &out,int dir,int disp)=0; // implemented by WilsonFermion and WilsonFermion5D
// force terms; five routines; default to Dhop on diagonal
virtual void MDeriv (LatticeGaugeField &mat,const FermionField &U,const FermionField &V,int dag){DhopDeriv(mat,U,V,dag);};
virtual void MoeDeriv(LatticeGaugeField &mat,const FermionField &U,const FermionField &V,int dag){DhopDerivOE(mat,U,V,dag);};
virtual void MeoDeriv(LatticeGaugeField &mat,const FermionField &U,const FermionField &V,int dag){DhopDerivEO(mat,U,V,dag);};
virtual void MooDeriv(LatticeGaugeField &mat,const FermionField &U,const FermionField &V,int dag){mat=zero;};
virtual void MeeDeriv(LatticeGaugeField &mat,const FermionField &U,const FermionField &V,int dag){mat=zero;};
virtual void DhopDeriv (LatticeGaugeField &mat,const FermionField &U,const FermionField &V,int dag)=0;
virtual void DhopDerivEO(LatticeGaugeField &mat,const FermionField &U,const FermionField &V,int dag)=0;
virtual void DhopDerivOE(LatticeGaugeField &mat,const FermionField &U,const FermionField &V,int dag)=0;
virtual void Mdiag (const FermionField &in, FermionField &out) { Mooee(in,out);}; // Same as Mooee applied to both CB's
virtual void Mdir (const FermionField &in, FermionField &out,int dir,int disp)=0; // case by case Wilson, Clover, Cayley, ContFrac, PartFrac
///////////////////////////////////////////////
// Updates gauge field during HMC
///////////////////////////////////////////////
virtual void ImportGauge(const GaugeField & _U);
};

2
lib/qcd/action/fermion/MobiusFermion.h
View File

@ -30,7 +30,7 @@ namespace Grid {
{
RealD eps = 1.0;
std::cout << "MobiusFermion (b="<<b<<",c="<<c<<") with Ls= "<<Ls<<" Tanh approx"<<std::endl;
std::cout<<GridLogMessage << "MobiusFermion (b="<<b<<",c="<<c<<") with Ls= "<<Ls<<" Tanh approx"<<std::endl;
Approx::zolotarev_data *zdata = Approx::higham(eps,this->Ls);// eps is ignored for higham
assert(zdata->n==this->Ls);

2
lib/qcd/action/fermion/MobiusZolotarevFermion.h
View File

@ -34,7 +34,7 @@ namespace Grid {
Approx::zolotarev_data *zdata = Approx::zolotarev(eps,this->Ls,0);
assert(zdata->n==this->Ls);
std::cout << "MobiusZolotarevFermion (b="<<b<<",c="<<c<<") with Ls= "<<Ls<<" Zolotarev range ["<<lo<<","<<hi<<"]"<<std::endl;
std::cout<<GridLogMessage << "MobiusZolotarevFermion (b="<<b<<",c="<<c<<") with Ls= "<<Ls<<" Zolotarev range ["<<lo<<","<<hi<<"]"<<std::endl;
// Call base setter
this->CayleyFermion5D::SetCoefficientsZolotarev(hi,zdata,b,c);

12
lib/qcd/action/fermion/PartialFractionFermion5D.cc
View File

@ -260,12 +260,12 @@ namespace Grid {
void PartialFractionFermion5D::SetCoefficientsZolotarev(RealD zolo_hi,Approx::zolotarev_data *zdata){
// check on degree matching
// std::cout << Ls << " Ls"<<std::endl;
// std::cout << zdata->n << " - n"<<std::endl;
// std::cout << zdata->da << " -da "<<std::endl;
// std::cout << zdata->db << " -db"<<std::endl;
// std::cout << zdata->dn << " -dn"<<std::endl;
// std::cout << zdata->dd << " -dd"<<std::endl;
// std::cout<<GridLogMessage << Ls << " Ls"<<std::endl;
// std::cout<<GridLogMessage << zdata->n << " - n"<<std::endl;
// std::cout<<GridLogMessage << zdata->da << " -da "<<std::endl;
// std::cout<<GridLogMessage << zdata->db << " -db"<<std::endl;
// std::cout<<GridLogMessage << zdata->dn << " -dn"<<std::endl;
// std::cout<<GridLogMessage << zdata->dd << " -dd"<<std::endl;
assert(Ls == (2*zdata->da -1) );
// Part frac

92
lib/qcd/action/fermion/WilsonFermion.cc
View File

@ -24,6 +24,10 @@ WilsonFermion::WilsonFermion(LatticeGaugeField &_Umu,
{
// Allocate the required comms buffer
comm_buf.resize(Stencil._unified_buffer_size); // this is always big enough to contain EO
ImportGauge(_Umu);
}
void WilsonFermion::ImportGauge(const LatticeGaugeField &_Umu)
{
DoubleStore(Umu,_Umu);
pickCheckerboard(Even,UmuEven,Umu);
pickCheckerboard(Odd ,UmuOdd,Umu);
@ -98,7 +102,9 @@ void WilsonFermion::Mdir (const LatticeFermion &in, LatticeFermion &out,int dir,
DhopDir(in,out,dir,disp);
}
void WilsonFermion::DhopDir(const LatticeFermion &in, LatticeFermion &out,int dir,int disp){
WilsonCompressor compressor(DaggerNo);
Stencil.HaloExchange<vSpinColourVector,vHalfSpinColourVector,WilsonCompressor>(in,comm_buf,compressor);
assert( (disp==1)||(disp==-1) );
@ -109,9 +115,22 @@ void WilsonFermion::DhopDir(const LatticeFermion &in, LatticeFermion &out,int di
PARALLEL_FOR_LOOP
for(int sss=0;sss<in._grid->oSites();sss++){
DiracOptDhopDir(Stencil,Umu,comm_buf,sss,sss,in,out,dirdisp);
DiracOptDhopDir(Stencil,Umu,comm_buf,sss,sss,in,out,dirdisp,dirdisp);
}
};
void WilsonFermion::DhopDirDisp(const LatticeFermion &in, LatticeFermion &out,int dirdisp,int gamma,int dag)
{
WilsonCompressor compressor(dag);
Stencil.HaloExchange<vSpinColourVector,vHalfSpinColourVector,WilsonCompressor>(in,comm_buf,compressor);
PARALLEL_FOR_LOOP
for(int sss=0;sss<in._grid->oSites();sss++){
DiracOptDhopDir(Stencil,Umu,comm_buf,sss,sss,in,out,dirdisp,gamma);
}
};
void WilsonFermion::DhopInternal(CartesianStencil & st,LatticeDoubledGaugeField & U,
@ -177,6 +196,77 @@ void WilsonFermion::Dhop(const LatticeFermion &in, LatticeFermion &out,int dag)
DhopInternal(Stencil,Umu,in,out,dag);
}
void WilsonFermion::DerivInternal(CartesianStencil & st,LatticeDoubledGaugeField & U,
LatticeGaugeField &mat,const LatticeFermion &A,const LatticeFermion &B,int dag)
{
assert((dag==DaggerNo) ||(dag==DaggerYes));
WilsonCompressor compressor(dag);
LatticeColourMatrix tmp(B._grid);
LatticeFermion Btilde(B._grid);
st.HaloExchange<vSpinColourVector,vHalfSpinColourVector,WilsonCompressor>(B,comm_buf,compressor);
for(int mu=0;mu<Nd;mu++){
////////////////////////////////////////////////////////////////////////
// Flip gamma (1+g)<->(1-g) if dag
////////////////////////////////////////////////////////////////////////
int gamma = mu;
if ( dag ) gamma+= Nd;
////////////////////////
// Call the single hop
////////////////////////
PARALLEL_FOR_LOOP
for(int sss=0;sss<B._grid->oSites();sss++){
DiracOptDhopDir(st,U,comm_buf,sss,sss,B,Btilde,mu,gamma);
}
//////////////////////////////////////////////////
// spin trace outer product
//////////////////////////////////////////////////
tmp = TraceIndex<SpinIndex>(outerProduct(Btilde,A));
PokeIndex<LorentzIndex>(mat,tmp,mu);
}
}
void WilsonFermion::DhopDeriv(LatticeGaugeField &mat,const LatticeFermion &U,const LatticeFermion &V,int dag)
{
conformable(U._grid,_grid);
conformable(U._grid,V._grid);
conformable(U._grid,mat._grid);
mat.checkerboard = U.checkerboard;
DerivInternal(Stencil,Umu,mat,U,V,dag);
}
void WilsonFermion::DhopDerivOE(LatticeGaugeField &mat,const LatticeFermion &U,const LatticeFermion &V,int dag)
{
conformable(U._grid,_cbgrid);
conformable(U._grid,V._grid);
conformable(U._grid,mat._grid);
assert(V.checkerboard==Even);
assert(U.checkerboard==Odd);
mat.checkerboard = Odd;
DerivInternal(StencilEven,UmuOdd,mat,U,V,dag);
}
void WilsonFermion::DhopDerivEO(LatticeGaugeField &mat,const LatticeFermion &U,const LatticeFermion &V,int dag)
{
conformable(U._grid,_cbgrid);
conformable(U._grid,V._grid);
conformable(U._grid,mat._grid);
assert(V.checkerboard==Odd);
assert(U.checkerboard==Even);
mat.checkerboard = Even;
DerivInternal(StencilOdd,UmuEven,mat,U,V,dag);
}
}}

89
lib/qcd/action/fermion/WilsonFermion.h
View File

@ -24,11 +24,95 @@ namespace Grid {
// half checkerboard operaions
void Meooe (const LatticeFermion &in, LatticeFermion &out);
void MeooeDag (const LatticeFermion &in, LatticeFermion &out);
virtual void Mooee (const LatticeFermion &in, LatticeFermion &out); // remain virtual so we
virtual void MooeeDag (const LatticeFermion &in, LatticeFermion &out); // can derive Clover
virtual void MooeeInv (const LatticeFermion &in, LatticeFermion &out); // from Wilson base
virtual void MooeeInvDag (const LatticeFermion &in, LatticeFermion &out);
////////////////////////
//
// Force term: d/dtau S = 0
//
// It is simplest to consider the two flavour force term
//
// S[U,phi] = phidag (MdagM)^-1 phi
//
// But simplify even this to
//
// S[U,phi] = phidag MdagM phi
//
// (other options exist depending on nature of action fragment.)
//
// Require momentum be traceless anti-hermitian to move within group manifold [ P = i P^a T^a ]
//
// Define the HMC hamiltonian
//
// H = 1/2 Tr P^2 + S(U,phi)
//
// .
// U = P U (lorentz & color indices multiplied)
//
// Hence
//
// .c c c c
// U = U P = - U P (c == dagger)
//
// So, taking some liberty with implicit indices
// . . .c c
// dH/dt = 0 = Tr P P +Tr[ U dS/dU + U dS/dU ]
//
// . c c
// = Tr P P + i Tr[ P U dS/dU - U P dS/dU ]
//
// . c c
// = Tr P (P + i ( U dS/dU - P dS/dU U ]
//
// . c c
// => P = -i [ U dS/dU - dS/dU U ] generates HMC EoM
//
// Simple case work this out using S = phi^dag MdagM phi for wilson:
// c c
// dSdt = dU_xdt dSdUx + dUxdt dSdUx
//
// = Tr i P U_x [ (\phi^\dag)_x (1+g) (M \phi)_x+\mu +(\phi^\dag M^\dag)_x (1-g) \phi_{x+\mu} ]
// c
// - i U_x P [ (\phi^\dag)_x+mu (1-g) (M \phi)_x +(\phi^\dag M^\dag)_(x+\mu) (1+g) \phi_{x} ]
//
// = i [(\phi^\dag)_x ]_j P_jk [U_x(1+g) (M \phi)_x+\mu]_k (1)
// + i [(\phi^\dagM^\dag)_x]_j P_jk [U_x(1-g) (\phi)_x+\mu]_k (2)
// - i [(\phi^\dag)_x+mu (1-g) U^dag_x]_j P_jk [(M \phi)_xk (3)
// - i [(\phi^\dagM^\dag)_x+mu (1+g) U^dag_x]_j P_jk [ \phi]_xk (4)
//
// Observe that (1)* = (4)
// (2)* = (3)
//
// Write as .
// P_{kj} = - i ( [U_x(1+g) (M \phi)_x+\mu] (x) [(\phi^\dag)_x] + [U_x(1-g) (\phi)_x+\mu] (x) [(\phi^\dagM^\dag)_x] - h.c )
//
// where (x) denotes outer product in colour and spins are traced.
//
// Need only evaluate (1) and (2) [Chroma] or (2) and (4) [IroIro] and take the
// traceless anti hermitian part (of term in brackets w/o the "i")
//
// Generalisation to S=phi^dag (MdagM)^{-1} phi is simple:
//
// For more complicated DWF etc... apply product rule in differentiation
//
////////////////////////
void DhopDeriv (LatticeGaugeField &mat,const LatticeFermion &U,const LatticeFermion &V,int dag);
void DhopDerivEO(LatticeGaugeField &mat,const LatticeFermion &U,const LatticeFermion &V,int dag);
void DhopDerivOE(LatticeGaugeField &mat,const LatticeFermion &U,const LatticeFermion &V,int dag);
// Extra support internal
void DerivInternal(CartesianStencil & st,
LatticeDoubledGaugeField & U,
LatticeGaugeField &mat,
const LatticeFermion &A,
const LatticeFermion &B,
int dag);
// non-hermitian hopping term; half cb or both
void Dhop (const LatticeFermion &in, LatticeFermion &out,int dag);
void DhopOE(const LatticeFermion &in, LatticeFermion &out,int dag);
@ -37,6 +121,7 @@ namespace Grid {
// Multigrid assistance
void Mdir (const LatticeFermion &in, LatticeFermion &out,int dir,int disp);
void DhopDir(const LatticeFermion &in, LatticeFermion &out,int dir,int disp);
void DhopDirDisp(const LatticeFermion &in, LatticeFermion &out,int dirdisp,int gamma,int dag);
///////////////////////////////////////////////////////////////
// Extra methods added by derived
@ -51,6 +136,7 @@ namespace Grid {
WilsonFermion(LatticeGaugeField &_Umu,GridCartesian &Fgrid,GridRedBlackCartesian &Hgrid,RealD _mass);
// DoubleStore
virtual void ImportGauge(const LatticeGaugeField &_Umu);
void DoubleStore(LatticeDoubledGaugeField &Uds,const LatticeGaugeField &Umu);
///////////////////////////////////////////////////////////////
@ -59,7 +145,8 @@ namespace Grid {
static int HandOptDslash; // these are a temporary hack
static int MortonOrder;
protected:
// protected:
public:
RealD mass;

103
lib/qcd/action/fermion/WilsonFermion5D.cc
View File

@ -65,7 +65,10 @@ namespace QCD {
// Allocate the required comms buffer
comm_buf.resize(Stencil._unified_buffer_size); // this is always big enough to contain EO
ImportGauge(_Umu);
}
void WilsonFermion5D::ImportGauge(const LatticeGaugeField &_Umu)
{
DoubleStore(Umu,_Umu);
pickCheckerboard(Even,UmuEven,Umu);
pickCheckerboard(Odd ,UmuOdd,Umu);
@ -100,19 +103,111 @@ void WilsonFermion5D::DhopDir(const LatticeFermion &in, LatticeFermion &out,int
assert(dirdisp<=7);
assert(dirdisp>=0);
//PARALLEL_FOR_LOOP
PARALLEL_FOR_LOOP
for(int ss=0;ss<Umu._grid->oSites();ss++){
for(int s=0;s<Ls;s++){
int sU=ss;
int sF = s+Ls*sU;
DiracOptDhopDir(Stencil,Umu,comm_buf,sF,sU,in,out,dirdisp);
DiracOptDhopDir(Stencil,Umu,comm_buf,sF,sU,in,out,dirdisp,dirdisp);
}
}
};
void WilsonFermion5D::DerivInternal(CartesianStencil & st,
LatticeDoubledGaugeField & U,
LatticeGaugeField &mat,
const LatticeFermion &A,
const LatticeFermion &B,
int dag)
{
assert((dag==DaggerNo) ||(dag==DaggerYes));
WilsonCompressor compressor(dag);
LatticeColourMatrix tmp(B._grid);
LatticeFermion Btilde(B._grid);
st.HaloExchange<vSpinColourVector,vHalfSpinColourVector,WilsonCompressor>(B,comm_buf,compressor);
for(int mu=0;mu<Nd;mu++){
////////////////////////////////////////////////////////////////////////
// Flip gamma if dag
////////////////////////////////////////////////////////////////////////
int gamma = mu;
if ( dag ) gamma+= Nd;
////////////////////////
// Call the single hop
////////////////////////
PARALLEL_FOR_LOOP
for(int sss=0;sss<B._grid->oSites();sss++){
for(int s=0;s<Ls;s++){
int sU=sss;
int sF = s+Ls*sU;
DiracOptDhopDir(st,U,comm_buf,sF,sU,B,Btilde,mu,gamma);
}
}
////////////////////////////
// spin trace outer product
////////////////////////////
// FIXME : need to sum over fifth direction.
tmp = TraceIndex<SpinIndex>(outerProduct(Btilde,A)); // ordering here
PokeIndex<LorentzIndex>(mat,tmp,mu);
}
}
void WilsonFermion5D::DhopDeriv( LatticeGaugeField &mat,
const LatticeFermion &A,
const LatticeFermion &B,
int dag)
{
conformable(A._grid,FermionGrid());
conformable(A._grid,B._grid);
conformable(GaugeGrid(),mat._grid);
mat.checkerboard = A.checkerboard;
DerivInternal(Stencil,Umu,mat,A,B,dag);
}
void WilsonFermion5D::DhopDerivEO(LatticeGaugeField &mat,
const LatticeFermion &A,
const LatticeFermion &B,
int dag)
{
conformable(A._grid,FermionRedBlackGrid());
conformable(GaugeRedBlackGrid(),mat._grid);
conformable(A._grid,B._grid);
assert(B.checkerboard==Odd);
assert(A.checkerboard==Even);
mat.checkerboard = Even;
DerivInternal(StencilOdd,UmuEven,mat,A,B,dag);
}
void WilsonFermion5D::DhopDerivOE(LatticeGaugeField &mat,
const LatticeFermion &A,
const LatticeFermion &B,
int dag)
{
conformable(A._grid,FermionRedBlackGrid());
conformable(GaugeRedBlackGrid(),mat._grid);
conformable(A._grid,B._grid);
assert(B.checkerboard==Even);
assert(A.checkerboard==Odd);
mat.checkerboard = Odd;
DerivInternal(StencilEven,UmuOdd,mat,A,B,dag);
}
void WilsonFermion5D::DhopInternal(CartesianStencil & st, LebesgueOrder &lo,
LatticeDoubledGaugeField & U,
const LatticeFermion &in, LatticeFermion &out,int dag)
const LatticeFermion &in, LatticeFermion &out,int dag)
{
// assert((dag==DaggerNo) ||(dag==DaggerYes));

19
lib/qcd/action/fermion/WilsonFermion5D.h
View File

@ -44,12 +44,18 @@ namespace Grid {
// half checkerboard operations; leave unimplemented as abstract for now
virtual void Meooe (const LatticeFermion &in, LatticeFermion &out){assert(0);};
virtual void MeooeDag (const LatticeFermion &in, LatticeFermion &out){assert(0);};
virtual void Mooee (const LatticeFermion &in, LatticeFermion &out){assert(0);};
virtual void MooeeDag (const LatticeFermion &in, LatticeFermion &out){assert(0);};
virtual void MooeeInv (const LatticeFermion &in, LatticeFermion &out){assert(0);};
virtual void MeooeDag (const LatticeFermion &in, LatticeFermion &out){assert(0);};
virtual void MooeeDag (const LatticeFermion &in, LatticeFermion &out){assert(0);};
virtual void MooeeInvDag (const LatticeFermion &in, LatticeFermion &out){assert(0);};
// These can be overridden by fancy 5d chiral actions
virtual void DhopDeriv (LatticeGaugeField &mat,const LatticeFermion &U,const LatticeFermion &V,int dag);
virtual void DhopDerivEO(LatticeGaugeField &mat,const LatticeFermion &U,const LatticeFermion &V,int dag);
virtual void DhopDerivOE(LatticeGaugeField &mat,const LatticeFermion &U,const LatticeFermion &V,int dag);
// Implement hopping term non-hermitian hopping term; half cb or both
// Implement s-diagonal DW
void DW (const LatticeFermion &in, LatticeFermion &out,int dag);
@ -64,6 +70,14 @@ namespace Grid {
///////////////////////////////////////////////////////////////
// New methods added
///////////////////////////////////////////////////////////////
void DerivInternal(CartesianStencil & st,
LatticeDoubledGaugeField & U,
LatticeGaugeField &mat,
const LatticeFermion &A,
const LatticeFermion &B,
int dag);
void DhopInternal(CartesianStencil & st,
LebesgueOrder &lo,
LatticeDoubledGaugeField &U,
@ -80,6 +94,7 @@ namespace Grid {
double _M5);
// DoubleStore
virtual void ImportGauge(const LatticeGaugeField &_Umu);
void DoubleStore(LatticeDoubledGaugeField &Uds,const LatticeGaugeField &Umu);
///////////////////////////////////////////////////////////////

44
lib/qcd/action/fermion/WilsonKernels.cc
View File

@ -294,8 +294,8 @@ void DiracOptDhopSiteDag(CartesianStencil &st,LatticeDoubledGaugeField &U,
}
void DiracOptDhopDir(CartesianStencil &st,LatticeDoubledGaugeField &U,
std::vector<vHalfSpinColourVector,alignedAllocator<vHalfSpinColourVector> > &buf,
int sF,int sU,const LatticeFermion &in, LatticeFermion &out,int dirdisp)
std::vector<vHalfSpinColourVector,alignedAllocator<vHalfSpinColourVector> > &buf,
int sF,int sU,const LatticeFermion &in, LatticeFermion &out,int dir,int gamma)
{
vHalfSpinColourVector tmp;
vHalfSpinColourVector chi;
@ -304,13 +304,13 @@ void DiracOptDhopDir(CartesianStencil &st,LatticeDoubledGaugeField &U,
int offset,local,perm, ptype;
int ss=sF;
offset = st._offsets [dirdisp][ss];
local = st._is_local[dirdisp][ss];
perm = st._permute[dirdisp][ss];
ptype = st._permute_type[dirdisp];
offset = st._offsets [dir][ss];
local = st._is_local[dir][ss];
perm = st._permute[dir][ss];
ptype = st._permute_type[dir];
// Xp
if(dirdisp==Xp){
if(gamma==Xp){
if ( local && perm ) {
spProjXp(tmp,in._odata[offset]);
permute(chi,tmp,ptype);
@ -319,12 +319,12 @@ void DiracOptDhopDir(CartesianStencil &st,LatticeDoubledGaugeField &U,
} else {
chi=buf[offset];
}
mult(&Uchi(),&U._odata[sU](Xp),&chi());
mult(&Uchi(),&U._odata[sU](dir),&chi());
spReconXp(result,Uchi);
}
// Yp
if ( dirdisp==Yp ){
if ( gamma==Yp ){
if ( local && perm ) {
spProjYp(tmp,in._odata[offset]);
permute(chi,tmp,ptype);
@ -333,12 +333,12 @@ void DiracOptDhopDir(CartesianStencil &st,LatticeDoubledGaugeField &U,
} else {
chi=buf[offset];
}
mult(&Uchi(),&U._odata[sU](Yp),&chi());
mult(&Uchi(),&U._odata[sU](dir),&chi());
spReconYp(result,Uchi);
}
// Zp
if ( dirdisp ==Zp ){
if ( gamma ==Zp ){
if ( local && perm ) {
spProjZp(tmp,in._odata[offset]);
permute(chi,tmp,ptype);
@ -347,12 +347,12 @@ void DiracOptDhopDir(CartesianStencil &st,LatticeDoubledGaugeField &U,
} else {
chi=buf[offset];
}
mult(&Uchi(),&U._odata[sU](Zp),&chi());
mult(&Uchi(),&U._odata[sU](dir),&chi());
spReconZp(result,Uchi);
}
// Tp
if ( dirdisp ==Tp ){
if ( gamma ==Tp ){
if ( local && perm ) {
spProjTp(tmp,in._odata[offset]);
permute(chi,tmp,ptype);
@ -361,12 +361,12 @@ void DiracOptDhopDir(CartesianStencil &st,LatticeDoubledGaugeField &U,
} else {
chi=buf[offset];
}
mult(&Uchi(),&U._odata[sU](Tp),&chi());
mult(&Uchi(),&U._odata[sU](dir),&chi());
spReconTp(result,Uchi);
}
// Xm
if ( dirdisp==Xm ){
if ( gamma==Xm ){
if ( local && perm ) {
spProjXm(tmp,in._odata[offset]);
permute(chi,tmp,ptype);
@ -375,12 +375,12 @@ void DiracOptDhopDir(CartesianStencil &st,LatticeDoubledGaugeField &U,
} else {
chi=buf[offset];
}
mult(&Uchi(),&U._odata[sU](Xm),&chi());
mult(&Uchi(),&U._odata[sU](dir),&chi());
spReconXm(result,Uchi);
}
// Ym
if ( dirdisp == Ym ){
if ( gamma == Ym ){
if ( local && perm ) {
spProjYm(tmp,in._odata[offset]);
permute(chi,tmp,ptype);
@ -389,12 +389,12 @@ void DiracOptDhopDir(CartesianStencil &st,LatticeDoubledGaugeField &U,
} else {
chi=buf[offset];
}
mult(&Uchi(),&U._odata[sU](Ym),&chi());
mult(&Uchi(),&U._odata[sU](dir),&chi());
spReconYm(result,Uchi);
}
// Zm
if ( dirdisp == Zm ){
if ( gamma == Zm ){
if ( local && perm ) {
spProjZm(tmp,in._odata[offset]);
permute(chi,tmp,ptype);
@ -403,12 +403,12 @@ void DiracOptDhopDir(CartesianStencil &st,LatticeDoubledGaugeField &U,
} else {
chi=buf[offset];
}
mult(&Uchi(),&U._odata[sU](Zm),&chi());
mult(&Uchi(),&U._odata[sU](dir),&chi());
spReconZm(result,Uchi);
}
// Tm
if ( dirdisp==Tm ) {
if ( gamma==Tm ) {
if ( local && perm ) {
spProjTm(tmp,in._odata[offset]);
permute(chi,tmp,ptype);
@ -417,7 +417,7 @@ void DiracOptDhopDir(CartesianStencil &st,LatticeDoubledGaugeField &U,
} else {
chi=buf[offset];
}
mult(&Uchi(),&U._odata[sU](Tm),&chi());
mult(&Uchi(),&U._odata[sU](dir),&chi());
spReconTm(result,Uchi);
}

2
lib/qcd/action/fermion/WilsonKernels.h
View File

@ -22,7 +22,7 @@ namespace Grid {
int sF,int sU,const LatticeFermion &in, LatticeFermion &out);
void DiracOptDhopDir(CartesianStencil &st,LatticeDoubledGaugeField &U,
std::vector<vHalfSpinColourVector,alignedAllocator<vHalfSpinColourVector> > &buf,
int sF,int sU,const LatticeFermion &in, LatticeFermion &out,int dirdisp);
int sF,int sU,const LatticeFermion &in, LatticeFermion &out,int dirdisp,int gamma);
// };

20
lib/qcd/action/fermion/WilsonKernelsHand.cc
View File

@ -360,11 +360,11 @@ void DiracOptHandDhopSite(CartesianStencil &st,LatticeDoubledGaugeField &U,
MULT_2SPIN(Xp);
}
XP_RECON;
// std::cout << "XP_RECON"<<std::endl;
// std::cout << result_00 <<" "<<result_01 <<" "<<result_02 <<std::endl;
// std::cout << result_10 <<" "<<result_11 <<" "<<result_12 <<std::endl;
// std::cout << result_20 <<" "<<result_21 <<" "<<result_22 <<std::endl;
// std::cout << result_30 <<" "<<result_31 <<" "<<result_32 <<std::endl;
// std::cout<<GridLogMessage << "XP_RECON"<<std::endl;
// std::cout<<GridLogMessage << result_00 <<" "<<result_01 <<" "<<result_02 <<std::endl;
// std::cout<<GridLogMessage << result_10 <<" "<<result_11 <<" "<<result_12 <<std::endl;
// std::cout<<GridLogMessage << result_20 <<" "<<result_21 <<" "<<result_22 <<std::endl;
// std::cout<<GridLogMessage << result_30 <<" "<<result_31 <<" "<<result_32 <<std::endl;
// Yp
offset = st._offsets [Yp][ss];
@ -446,11 +446,11 @@ void DiracOptHandDhopSite(CartesianStencil &st,LatticeDoubledGaugeField &U,
MULT_2SPIN(Xm);
}
XM_RECON_ACCUM;
// std::cout << "XM_RECON_ACCUM"<<std::endl;
// std::cout << result_00 <<" "<<result_01 <<" "<<result_02 <<std::endl;
// std::cout << result_10 <<" "<<result_11 <<" "<<result_12 <<std::endl;
// std::cout << result_20 <<" "<<result_21 <<" "<<result_22 <<std::endl;
// std::cout << result_30 <<" "<<result_31 <<" "<<result_32 <<std::endl;
// std::cout<<GridLogMessage << "XM_RECON_ACCUM"<<std::endl;
// std::cout<<GridLogMessage << result_00 <<" "<<result_01 <<" "<<result_02 <<std::endl;
// std::cout<<GridLogMessage << result_10 <<" "<<result_11 <<" "<<result_12 <<std::endl;
// std::cout<<GridLogMessage << result_20 <<" "<<result_21 <<" "<<result_22 <<std::endl;
// std::cout<<GridLogMessage << result_30 <<" "<<result_31 <<" "<<result_32 <<std::endl;
// Ym

8
lib/qcd/action/gauge/WilsonGaugeAction.h
View File

@ -18,9 +18,11 @@ namespace Grid{
virtual RealD S(const GaugeField &U) {
RealD plaq = WilsonLoops<MatrixField,GaugeField>::avgPlaquette(U);
std::cout << "Plaq : "<<plaq << "\n";
double vol = U._grid->gSites();
return beta*(1.0 -plaq)*(Nd*(Nd-1.0))*vol*0.5;
std::cout<<GridLogMessage << "Plaq : "<<plaq << "\n";
RealD vol = U._grid->gSites();
RealD action=beta*(1.0 -plaq)*(Nd*(Nd-1.0))*vol*0.5;
std::cout << GridLogMessage << "WilsonGauge action "<<action<<std::endl;
return action;
};
virtual void deriv(const GaugeField &U,GaugeField & dSdU) {
//not optimal implementation FIXME

209
lib/qcd/action/pseudofermion/TwoFlavour.h Normal file
View File

@ -0,0 +1,209 @@
#ifndef QCD_PSEUDOFERMION_TWO_FLAVOUR_H
#define QCD_PSEUDOFERMION_TWO_FLAVOUR_H
namespace Grid{
namespace QCD{
// Placeholder comments:
///////////////////////////////////////
// Two flavour ratio
///////////////////////////////////////
// S = phi^dag V (Mdag M)^-1 V^dag phi
// dS/du = phi^dag dV (Mdag M)^-1 V^dag phi
// - phi^dag V (Mdag M)^-1 [ Mdag dM + dMdag M ] (Mdag M)^-1 V^dag phi
// + phi^dag V (Mdag M)^-1 dV^dag phi
///////////////////////////////////////
// One flavour rational
///////////////////////////////////////
// S_f = chi^dag * N(M^dag*M)/D(M^dag*M) * chi
//
// Here, M is some operator
// N and D makeup the rat. poly
//
// Need
// dS_f/dU = chi^dag P/Q d[N/D] P/Q chi
//
// Here N/D \sim R_{-1/2} ~ (M^dagM)^{-1/2}
//
// N/D is expressed as partial fraction expansion:
//
// a0 + \sum_k ak/(M^dagM + bk)
//
// d[N/D] is then
//
// \sum_k -ak [M^dagM+bk]^{-1} [ dM^dag M + M^dag dM ] [M^dag M + bk]^{-1}
//
// Need
//
// Mf Phi_k = [MdagM+bk]^{-1} Phi
// Mf Phi = \sum_k ak [MdagM+bk]^{-1} Phi
//
// With these building blocks
//
// dS/dU = \sum_k -ak Mf Phi_k^dag [ dM^dag M + M^dag dM ] Mf Phi_k
// S = innerprodReal(Phi,Mf Phi);
///////////////////////////////////////
// One flavour rational ratio
///////////////////////////////////////
// S_f = chi^dag* P(V^dag*V)/Q(V^dag*V)* N(M^dag*M)/D(M^dag*M)* P(V^dag*V)/Q(V^dag*V)* chi
//
// Here, M is some 5D operator and V is the Pauli-Villars field
// N and D makeup the rat. poly of the M term and P and & makeup the rat.poly of the denom term
//
// Need
// dS_f/dU = chi^dag d[P/Q] N/D P/Q chi
// + chi^dag P/Q d[N/D] P/Q chi
// + chi^dag P/Q N/D d[P/Q] chi
//
// Here P/Q \sim R_{1/4} ~ (V^dagV)^{1/4}
// Here N/D \sim R_{-1/2} ~ (M^dagM)^{-1/2}
//
// P/Q is expressed as partial fraction expansion:
//
// a0 + \sum_k ak/(V^dagV + bk)
//
// d[P/Q] is then
//
// \sum_k -ak [V^dagV+bk]^{-1} [ dV^dag V + V^dag dV ] [V^dag V + bk]^{-1}
//
// and similar for N/D.
//
// Need
// MpvPhi_k = [Vdag V + bk]^{-1} chi
//
// MpvPhi = {a0 + \sum_k ak [Vdag V + bk]^{-1} }chi
//
// MfMpvPhi_k = [MdagM+bk]^{-1} MpvPhi
//
// MfMpvPhi = {a0 + \sum_k ak [Mdag M + bk]^{-1} } MpvPhi
//
// MpvMfMpvPhi_k = [Vdag V + bk]^{-1} MfMpvchi
//
// With these building blocks
//
// dS/dU =
// \sum_k -ak MpvPhi_k^dag [ dV^dag V + V^dag dV ] MpvMfMpvPhi_k <- deriv on P left
// + \sum_k -ak MpvMfMpvPhi_k^\dag [ dV^dag V + V^dag dV ] MpvPhi_k
// + \sum_k -ak MfMpvPhi_k^dag [ dM^dag M + M^dag dM ] MfMpvPhi_k
////////////////////////////////////////////////////////////////////////
// Two flavour pseudofermion action for any dop
////////////////////////////////////////////////////////////////////////
template<class GaugeField,class MatrixField,class FermionField>
class TwoFlavourPseudoFermionAction : public Action<GaugeField> {
private:
FermionOperator<FermionField,GaugeField> & FermOp;// the basic operator
OperatorFunction<FermionField> &DerivativeSolver;
OperatorFunction<FermionField> &ActionSolver;
GridBase &Grid;
FermionField Phi; // the pseudo fermion field for this trajectory
public:
/////////////////////////////////////////////////
// Pass in required objects.
/////////////////////////////////////////////////
TwoFlavourPseudoFermionAction(FermionOperator<FermionField,GaugeField> &Op,
OperatorFunction<FermionField> & DS,
OperatorFunction<FermionField> & AS,
GridBase &_Grid
) : FermOp(Op), DerivativeSolver(DS), ActionSolver(AS), Phi(&_Grid), Grid(_Grid) {
};
//////////////////////////////////////////////////////////////////////////////////////
// Push the gauge field in to the dops. Assume any BC's and smearing already applied
//////////////////////////////////////////////////////////////////////////////////////
virtual void init(const GaugeField &U, GridParallelRNG& pRNG) {
// P(phi) = e^{- phi^dag (MdagM)^-1 phi}
// Phi = Mdag eta
// P(eta) = e^{- eta^dag eta}
//
// e^{x^2/2 sig^2} => sig^2 = 0.5.
//
// So eta should be of width sig = 1/sqrt(2).
// and must multiply by 0.707....
//
// Chroma has this scale factor: two_flavor_monomial_w.h
// IroIro: does not use this scale. It is absorbed by a change of vars
// in the Phi integral, and thus is only an irrelevant prefactor for the partition function.
//
RealD scale = std::sqrt(0.5);
FermionField eta(&Grid);
gaussian(pRNG,eta);
FermOp.Mdag(eta,Phi);
Phi=Phi*scale;
};
//////////////////////////////////////////////////////
// S = phi^dag (Mdag M)^-1 phi
//////////////////////////////////////////////////////
virtual RealD S(const GaugeField &U) {
FermOp.ImportGauge(U);
FermionField X(&Grid);
FermionField Y(&Grid);
MdagMLinearOperator<FermionOperator<FermionField,GaugeField> ,FermionField> MdagMOp(FermOp);
X=zero;
ActionSolver(MdagMOp,Phi,X);
MdagMOp.Op(X,Y);
RealD action = norm2(Y);
std::cout << GridLogMessage << "Pseudofermion action "<<action<<std::endl;
return action;
};
//////////////////////////////////////////////////////
// dS/du = - phi^dag (Mdag M)^-1 [ Mdag dM + dMdag M ] (Mdag M)^-1 phi
// = - phi^dag M^-1 dM (MdagM)^-1 phi - phi^dag (MdagM)^-1 dMdag dM (Mdag)^-1 phi
//
// = - Ydag dM X - Xdag dMdag Y
//
//////////////////////////////////////////////////////
virtual void deriv(const GaugeField &U,GaugeField & dSdU) {
FermOp.ImportGauge(U);
FermionField X(&Grid);
FermionField Y(&Grid);
GaugeField tmp(&Grid);
MdagMLinearOperator<FermionOperator<FermionField,GaugeField> ,FermionField> MdagMOp(FermOp);
X=zero;
DerivativeSolver(MdagMOp,Phi,X);
MdagMOp.Op(X,Y);
// Our conventions really make this UdSdU; We do not differentiate wrt Udag here.
// So must take dSdU - adj(dSdU) and left multiply by mom to get dS/dt.
FermOp.MDeriv(tmp , Y, X,DaggerNo ); dSdU=tmp;
FermOp.MDeriv(tmp , X, Y,DaggerYes); dSdU=dSdU+tmp;
dSdU = Ta(dSdU);
};
};
}
}
#endif

7
lib/qcd/hmc/HMC.cc
View File

@ -7,8 +7,8 @@ namespace Grid{
// FIXME fill this constructor now just default values
////////////////////////////// Default values
Nsweeps = 100;
TotalSweeps = 20;
Nsweeps = 200;
TotalSweeps = 220;
ThermalizationSteps = 20;
StartingConfig = 0;
SaveInterval = 1;
@ -17,8 +17,5 @@ namespace Grid{
}
}
}

58
lib/qcd/hmc/HMC.h
View File

@ -3,7 +3,7 @@
* @brief Classes for Hybrid Monte Carlo update
*
* @author Guido Cossu
* Time-stamp: <2015-07-07 14:58:13 neo>
* Time-stamp: <2015-07-30 16:58:26 neo>
*/
//--------------------------------------------------------------------
#ifndef HMC_INCLUDED
@ -28,75 +28,89 @@ namespace Grid{
template <class Algorithm>
class HybridMonteCarlo{
const HMCparameters Params;
GridSerialRNG sRNG;
GridSerialRNG sRNG; // Fixme: need a RNG management strategy.
Integrator<Algorithm>& MD;
/////////////////////////////////////////////////////////
// Metropolis step
/////////////////////////////////////////////////////////
bool metropolis_test(const RealD DeltaH){
RealD rn_test;
RealD prob = std::exp(-DeltaH);
random(sRNG,rn_test);
std::cout<< "--------------------------------------------\n";
std::cout<< "dH = "<<DeltaH << " Random = "<< rn_test
<< "\nAcc. Probability = " << ((prob<1.0)? prob: 1.0)<< " ";
std::cout<<GridLogMessage<< "--------------------------------------------\n";
std::cout<<GridLogMessage<< "dH = "<<DeltaH << " Random = "<< rn_test <<"\n";
std::cout<<GridLogMessage<< "Acc. Probability = " << ((prob<1.0)? prob: 1.0)<< " ";
if((prob >1.0) || (rn_test <= prob)){ // accepted
std::cout <<"-- ACCEPTED\n";
std::cout<<GridLogMessage <<"-- ACCEPTED\n";
return true;
} else { // rejected
std::cout <<"-- REJECTED\n";
std::cout<<GridLogMessage <<"-- REJECTED\n";
return false;
}
}
/////////////////////////////////////////////////////////
// Evolution
/////////////////////////////////////////////////////////
RealD evolve_step(LatticeGaugeField& U){
MD.init(U); // set U and initialize P and phi's
RealD H0 = MD.S(U); // initial state action
std::cout<<"Total H before = "<< H0 << "\n";
std::cout<<GridLogMessage<<"Total H before = "<< H0 << "\n";
MD.integrate(U);
RealD H1 = MD.S(U); // updated state action
std::cout<<"Total H after = "<< H1 << "\n";
std::cout<<GridLogMessage<<"Total H after = "<< H1 << "\n";
return (H1-H0);
}
public:
HybridMonteCarlo(HMCparameters Pms,
Integrator<Algorithm>& MolDyn):
Params(Pms),MD(MolDyn){
//FIXME
// initialize RNGs also with seed
/////////////////////////////////////////
// Constructor
/////////////////////////////////////////
HybridMonteCarlo(HMCparameters Pms, Integrator<Algorithm>& MolDyn): Params(Pms),MD(MolDyn) {
//FIXME... initialize RNGs also with seed ; RNG management strategy
sRNG.SeedRandomDevice();
}
~HybridMonteCarlo(){};
void evolve(LatticeGaugeField& Uin){
Real DeltaH;
// Thermalizations
for(int iter=1; iter <= Params.ThermalizationSteps; ++iter){
std::cout << "-- # Thermalization step = "<< iter << "\n";
std::cout<<GridLogMessage << "-- # Thermalization step = "<< iter << "\n";
DeltaH = evolve_step(Uin);
std::cout<< " dH = "<< DeltaH << "\n";
std::cout<<GridLogMessage<< "dH = "<< DeltaH << "\n";
}
// Actual updates (evolve a copy Ucopy then copy back eventually)
LatticeGaugeField Ucopy(Uin._grid);
for(int iter=Params.StartingConfig;
iter < Params.Nsweeps+Params.StartingConfig; ++iter){
std::cout << "-- # Sweep = "<< iter << "\n";
std::cout<<GridLogMessage << "-- # Sweep = "<< iter << "\n";
Ucopy = Uin;
DeltaH = evolve_step(Ucopy);
if(metropolis_test(DeltaH)) Uin = Ucopy;
// here save config and RNG seed
}
}
};

13
lib/qcd/hmc/integrators/Integrator.h
View File

@ -3,7 +3,7 @@
* @brief Classes for the Molecular Dynamics integrator
*
* @author Guido Cossu
* Time-stamp: <2015-07-07 14:58:40 neo>
* Time-stamp: <2015-07-30 16:21:29 neo>
*/
//--------------------------------------------------------------------
@ -71,7 +71,6 @@ namespace Grid{
}
}
void update_U(LatticeGaugeField&U, double ep){
//rewrite exponential to deal automatically with the lorentz index?
LatticeColourMatrix Umu(U._grid);
@ -85,7 +84,6 @@ namespace Grid{
}
friend void IntegratorAlgorithm::step (LatticeGaugeField& U,
int level, std::vector<int>& clock,
@ -99,11 +97,9 @@ namespace Grid{
~Integrator(){}
//Initialization of momenta and actions
void init(LatticeGaugeField& U){
std::cout<< "Integrator init\n";
std::cout<<GridLogMessage<< "Integrator init\n";
MDutils::generate_momenta(*P,pRNG);
for(int level=0; level< as.size(); ++level){
for(int actionID=0; actionID<as[level].actions.size(); ++actionID){
@ -112,7 +108,6 @@ namespace Grid{
}
}
// Calculate action
RealD S(LatticeGaugeField& U){
LatticeComplex Hloc(U._grid);
@ -126,12 +121,14 @@ namespace Grid{
RealD H = Hsum.real();
std::cout << "H_p = "<< H << "\n";
std::cout<<GridLogMessage << "Momentum action H_p = "<< H << "\n";
// Actions
for(int level=0; level<as.size(); ++level)
for(int actionID=0; actionID<as[level].actions.size(); ++actionID)
H += as[level].actions.at(actionID)->S(U);
std::cout<<GridLogMessage << "Total action H = "<< H << "\n";
return H;
}

56
lib/qcd/hmc/integrators/Integrator_algorithm.h
View File

@ -33,35 +33,32 @@ namespace Grid{
for(int l=1; l<=level; ++l) fin*= 2.0*Integ->as[l].multiplier;
fin = 3*Integ->Params.MDsteps*fin -1;
for(int e=0; e<Integ->as[level].multiplier; ++e){
if(clock[level] == 0){ // initial half step
Integ->update_P(U,level,lambda*eps);
++clock[level];
for(int l=0; l<level;++l) std::cout<<" ";
std::cout<<"P "<< clock[level] <<std::endl;
for(int l=0; l<level;++l) std::cout<<GridLogMessage<<" ";
std::cout<<GridLogMessage<<"P "<< clock[level] <<std::endl;
}
if(level == fl){ // lowest level
Integ->update_U(U,0.5*eps);
for(int l=0; l<level;++l) std::cout<<" ";
std::cout<<"U "<< (clock[level]+1) <<std::endl;
for(int l=0; l<level;++l) std::cout<<GridLogMessage<<" ";
std::cout<<GridLogMessage<<"U "<< (clock[level]+1) <<std::endl;
}else{ // recursive function call
step(U,level+1,clock, Integ);
}
Integ->update_P(U,level,(1.0-2.0*lambda)*eps);
++clock[level];
for(int l=0; l<level;++l) std::cout<<" ";
std::cout<<"P "<< (clock[level]) <<std::endl;
for(int l=0; l<level;++l) std::cout<<GridLogMessage<<" ";
std::cout<<GridLogMessage<<"P "<< (clock[level]) <<std::endl;
if(level == fl){ // lowest level
Integ->update_U(U,0.5*eps);
for(int l=0; l<level;++l) std::cout<<" ";
std::cout<<"U "<< (clock[level]+1) <<std::endl;
for(int l=0; l<level;++l) std::cout<<GridLogMessage<<" ";
std::cout<<GridLogMessage<<"U "<< (clock[level]+1) <<std::endl;
}else{ // recursive function call
step(U,level+1,clock, Integ);
}
@ -71,19 +68,17 @@ namespace Grid{
Integ->update_P(U,level,lambda*eps);
++clock[level];
for(int l=0; l<level;++l) std::cout<<" ";
std::cout<<"P "<< clock[level] <<std::endl;
for(int l=0; l<level;++l) std::cout<<GridLogMessage<<" ";
std::cout<<GridLogMessage<<"P "<< clock[level] <<std::endl;
}else{ // bulk step
Integ->update_P(U,level,lambda*2.0*eps);
clock[level]+=2;
for(int l=0; l<level;++l) std::cout<<" ";
std::cout<<"P "<< clock[level] <<std::endl;
for(int l=0; l<level;++l) std::cout<<GridLogMessage<<" ";
std::cout<<GridLogMessage<<"P "<< clock[level] <<std::endl;
}
}
}
};
@ -93,6 +88,7 @@ namespace Grid{
void step (LatticeLorentzColourMatrix& U,
int level, std::vector<int>& clock,
Integrator<LeapFrog>* Integ){
// level : current level
// fl : final level
// eps : current step size
@ -112,34 +108,32 @@ namespace Grid{
if(clock[level] == 0){ // initial half step
Integ->update_P(U, level,eps/2.0);
++clock[level];
for(int l=0; l<level;++l) std::cout<<" ";
std::cout<<"P "<< 0.5*clock[level] <<std::endl;
for(int l=0; l<level;++l) std::cout<<GridLogMessage<<" ";
std::cout<<GridLogMessage<<"P "<< 0.5*clock[level] <<std::endl;
}
if(level == fl){ // lowest level
Integ->update_U(U, eps);
for(int l=0; l<level;++l) std::cout<<" ";
std::cout<<"U "<< 0.5*(clock[level]+1) <<std::endl;
for(int l=0; l<level;++l) std::cout<<GridLogMessage<<" ";
std::cout<<GridLogMessage<<"U "<< 0.5*(clock[level]+1) <<std::endl;
}else{ // recursive function call
step(U, level+1,clock, Integ);
}
if(clock[level] == fin){ // final half step
Integ->update_P(U, level,eps/2.0);
++clock[level];
for(int l=0; l<level;++l) std::cout<<" ";
std::cout<<"P "<< 0.5*clock[level] <<std::endl;
for(int l=0; l<level;++l) std::cout<<GridLogMessage<<" ";
std::cout<<GridLogMessage<<"P "<< 0.5*clock[level] <<std::endl;
}else{ // bulk step
Integ->update_P(U, level,eps);
clock[level]+=2;
for(int l=0; l<level;++l) std::cout<<" ";
std::cout<<"P "<< 0.5*clock[level] <<std::endl;
for(int l=0; l<level;++l) std::cout<<GridLogMessage<<" ";
std::cout<<GridLogMessage<<"P "<< 0.5*clock[level] <<std::endl;
}
}
}
};

30
lib/qcd/hmc/integrators/Integrator_base.h
View File

@ -65,6 +65,18 @@ namespace Grid{
for(int a=0; a<as[level].size(); ++a){
LatticeLorentzColourMatrix force(U._grid);
as[level].at(a)->deriv(U,force);
Complex dSdt=0.0;
for(int mu=0;mu<Nd;mu++){
LatticeColourMatrix forcemu(U._grid);
LatticeColourMatrix mommu(U._grid);
forcemu=PeekIndex<LorentzIndex>(force,mu);
mommu=PeekIndex<LorentzIndex>(*P,mu);
dSdt += sum(trace(forcemu*(*P)));
}
std::cout << GridLogMessage << " action "<<level<<","<<a<<" dSdt "<< dSdt << " dt "<<ep <<std::endl;
*P -= force*ep;
}
}
@ -101,7 +113,7 @@ namespace Grid{
//Initialization of momenta and actions
void init(LatticeLorentzColourMatrix& U,
GridParallelRNG& pRNG){
std::cout<< "Integrator init\n";
std::cout<<GridLogMessage<< "Integrator init\n";
if (!P)
P = new LatticeLorentzColourMatrix(U._grid);
MDutils::generate_momenta(*P,pRNG);
@ -172,13 +184,13 @@ namespace Grid{
if(clock[level] == 0){ // initial half step
Integ->update_P(U, level,eps/2);
++clock[level];
for(int l=0; l<level;++l) std::cout<<" ";
std::cout<<"P "<< 0.5*clock[level] <<std::endl;
for(int l=0; l<level;++l) std::cout<<GridLogMessage<<" ";
std::cout<<GridLogMessage<<"P "<< 0.5*clock[level] <<std::endl;
}
if(level == fl){ // lowest level
Integ->update_U(U, eps);
for(int l=0; l<level;++l) std::cout<<" ";
std::cout<<"U "<< 0.5*(clock[level]+1) <<std::endl;
for(int l=0; l<level;++l) std::cout<<GridLogMessage<<" ";
std::cout<<GridLogMessage<<"U "<< 0.5*(clock[level]+1) <<std::endl;
}else{ // recursive function call
step(U, level+1,clock, Integ);
}
@ -186,14 +198,14 @@ namespace Grid{
Integ->update_P(U, level,eps/2);
++clock[level];
for(int l=0; l<level;++l) std::cout<<" ";
std::cout<<"P "<< 0.5*clock[level] <<std::endl;
for(int l=0; l<level;++l) std::cout<<GridLogMessage<<" ";
std::cout<<GridLogMessage<<"P "<< 0.5*clock[level] <<std::endl;
}else{ // bulk step
Integ->update_P(U, level,eps);
clock[level]+=2;
for(int l=0; l<level;++l) std::cout<<" ";
std::cout<<"P "<< 0.5*clock[level] <<std::endl;
for(int l=0; l<level;++l) std::cout<<GridLogMessage<<" ";
std::cout<<GridLogMessage<<"P "<< 0.5*clock[level] <<std::endl;
}
}

30
lib/qcd/utils/SUn.h
View File

@ -372,7 +372,7 @@ Note that in step D setting B ~ X - A and using B in place of A in step E will g
LatticeReal d(grid); d=zero;
LatticeReal alpha(grid);
// std::cout<<"xi "<<xi <<std::endl;
// std::cout<<GridLogMessage<<"xi "<<xi <<std::endl;
alpha = toReal(2.0*xi);
do {
@ -468,11 +468,11 @@ Note that in step D setting B ~ X - A and using B in place of A in step E will g
LatticeMatrix Vcheck(grid);
Vcheck = zero;
Vcheck = where(Accepted,V*adj(V) - 1.0,Vcheck);
// std::cout << "SU3 check " <<norm2(Vcheck)<<std::endl;
// std::cout<<GridLogMessage << "SU3 check " <<norm2(Vcheck)<<std::endl;
assert(norm2(Vcheck)<1.0e-4);
// Verify the link stays in SU(3)
// std::cout <<"Checking the modified link"<<std::endl;
// std::cout<<GridLogMessage <<"Checking the modified link"<<std::endl;
Vcheck = link*adj(link) - 1.0;
assert(norm2(Vcheck)<1.0e-4);
/////////////////////////////////
@ -483,42 +483,42 @@ Note that in step D setting B ~ X - A and using B in place of A in step E will g
for(int gen=0;gen<generators();gen++){
Matrix ta;
generator(gen,ta);
std::cout<< "Nc = "<<ncolour<<" t_"<<gen<<std::endl;
std::cout<<ta<<std::endl;
std::cout<<GridLogMessage<< "Nc = "<<ncolour<<" t_"<<gen<<std::endl;
std::cout<<GridLogMessage<<ta<<std::endl;
}
}
static void testGenerators(void){
Matrix ta;
Matrix tb;
std::cout<<"Checking trace ta tb is 0.5 delta_ab"<<std::endl;
std::cout<<GridLogMessage<<"Checking trace ta tb is 0.5 delta_ab"<<std::endl;
for(int a=0;a<generators();a++){
for(int b=0;b<generators();b++){
generator(a,ta);
generator(b,tb);
Complex tr =TensorRemove(trace(ta*tb));
std::cout<<tr<<" ";
std::cout<<GridLogMessage<<tr<<" ";
if(a==b) assert(abs(tr-Complex(0.5))<1.0e-6);
if(a!=b) assert(abs(tr)<1.0e-6);
}
std::cout<<std::endl;
std::cout<<GridLogMessage<<std::endl;
}
std::cout<<"Checking hermitian"<<std::endl;
std::cout<<GridLogMessage<<"Checking hermitian"<<std::endl;
for(int a=0;a<generators();a++){
generator(a,ta);
std::cout<<a<<" ";
std::cout<<GridLogMessage<<a<<" ";
assert(norm2(ta-adj(ta))<1.0e-6);
}
std::cout<<std::endl;
std::cout<<GridLogMessage<<std::endl;
std::cout<<"Checking traceless"<<std::endl;
std::cout<<GridLogMessage<<"Checking traceless"<<std::endl;
for(int a=0;a<generators();a++){
generator(a,ta);
Complex tr =TensorRemove(trace(ta));
std::cout<<a<<" ";
std::cout<<GridLogMessage<<a<<" ";
assert(abs(tr)<1.0e-6);
}
std::cout<<std::endl;
std::cout<<GridLogMessage<<std::endl;
}
// reunitarise??
@ -554,9 +554,7 @@ Note that in step D setting B ~ X - A and using B in place of A in step E will g
for(int a=0;a<generators();a++){
gaussian(pRNG,ca);
generator(a,ta);
la=toComplex(ca)*ci*ta;
out += la;
}