mirror of
https://github.com/paboyle/Grid.git
synced 2024-11-10 07:55:35 +00:00
Two flavour HMC for Wilson/Wilson is conserving energy.
Still to check plaq and <e(-dH)>, but nevertheless this is progress
This commit is contained in:
Peter Boyle
parent
4cc2ef84d3
commit
4fe110bd07
@ -1,4 +1,4 @@
|
||||
|
||||
HFILES=./algorithms/approx/bigfloat.h ./algorithms/approx/bigfloat_double.h ./algorithms/approx/Chebyshev.h ./algorithms/approx/MultiShiftFunction.h ./algorithms/approx/Remez.h ./algorithms/approx/Zolotarev.h ./algorithms/CoarsenedMatrix.h ./algorithms/iterative/AdefGeneric.h ./algorithms/iterative/BfmHDCG.h ./algorithms/iterative/ConjugateGradient.h ./algorithms/iterative/ConjugateGradientMultiShift.h ./algorithms/iterative/ConjugateResidual.h ./algorithms/iterative/GeneralisedMinimumResidual.h ./algorithms/iterative/gmres.h ./algorithms/iterative/NormalEquations.h ./algorithms/iterative/PrecConjugateResidual.h ./algorithms/iterative/PrecGeneralisedConjugateResidual.h ./algorithms/iterative/SchurRedBlack.h ./algorithms/LinearOperator.h ./algorithms/Preconditioner.h ./algorithms/SparseMatrix.h ./Algorithms.h ./AlignedAllocator.h ./cartesian/Cartesian_base.h ./cartesian/Cartesian_full.h ./cartesian/Cartesian_red_black.h ./Cartesian.h ./communicator/Communicator_base.h ./Communicator.h ./Config.h ./cshift/Cshift_common.h ./cshift/Cshift_mpi.h ./cshift/Cshift_none.h ./Cshift.h ./Grid.h ./Init.h ./lattice/Lattice_arith.h ./lattice/Lattice_base.h ./lattice/Lattice_comparison.h ./lattice/Lattice_comparison_utils.h ./lattice/Lattice_conformable.h ./lattice/Lattice_coordinate.h ./lattice/Lattice_ET.h ./lattice/Lattice_local.h ./lattice/Lattice_overload.h ./lattice/Lattice_peekpoke.h ./lattice/Lattice_reality.h ./lattice/Lattice_reduction.h ./lattice/Lattice_rng.h ./lattice/Lattice_trace.h ./lattice/Lattice_transfer.h ./lattice/Lattice_transpose.h ./lattice/Lattice_unary.h ./lattice/Lattice_where.h ./Lattice.h ./Log.h ./MacroMagic.h ./Old/Tensor_peek.h ./Old/Tensor_poke.h ./parallelIO/NerscIO.h ./qcd/action/ActionBase.h ./qcd/action/Actions.h ./qcd/action/DiffAction.h ./qcd/action/fermion/CayleyFermion5D.h ./qcd/action/fermion/ContinuedFractionFermion5D.h ./qcd/action/fermion/DomainWallFermion.h ./qcd/action/fermion/FermionOperator.h ./qcd/action/fermion/g5HermitianLinop.h ./qcd/action/fermion/MobiusFermion.h ./qcd/action/fermion/MobiusZolotarevFermion.h ./qcd/action/fermion/OverlapWilsonCayleyTanhFermion.h ./qcd/action/fermion/OverlapWilsonCayleyZolotarevFermion.h ./qcd/action/fermion/OverlapWilsonContfracTanhFermion.h ./qcd/action/fermion/OverlapWilsonContfracZolotarevFermion.h ./qcd/action/fermion/OverlapWilsonPartialFractionTanhFermion.h ./qcd/action/fermion/OverlapWilsonPartialFractionZolotarevFermion.h ./qcd/action/fermion/PartialFractionFermion5D.h ./qcd/action/fermion/ScaledShamirFermion.h ./qcd/action/fermion/ShamirZolotarevFermion.h ./qcd/action/fermion/WilsonCompressor.h ./qcd/action/fermion/WilsonFermion.h ./qcd/action/fermion/WilsonFermion5D.h ./qcd/action/fermion/WilsonKernels.h ./qcd/action/gauge/WilsonGaugeAction.h ./qcd/action/pseudofermion/TwoFlavour.h ./qcd/hmc/HMC.h ./qcd/hmc/integrators/Integrator.h ./qcd/hmc/integrators/Integrator_algorithm.h ./qcd/hmc/integrators/Integrator_base.h ./qcd/QCD.h ./qcd/spin/Dirac.h ./qcd/spin/TwoSpinor.h ./qcd/utils/CovariantCshift.h ./qcd/utils/LinalgUtils.h ./qcd/utils/SpaceTimeGrid.h ./qcd/utils/SUn.h ./qcd/utils/WilsonLoops.h ./simd/Grid_avx.h ./simd/Grid_avx512.h ./simd/Grid_empty.h ./simd/Grid_neon.h ./simd/Grid_qpx.h ./simd/Grid_sse4.h ./simd/Grid_vector_types.h ./simd/Grid_vector_unops.h ./Simd.h ./stencil/Lebesgue.h ./Stencil.h ./tensors/Tensor_arith.h ./tensors/Tensor_arith_add.h ./tensors/Tensor_arith_mac.h ./tensors/Tensor_arith_mul.h ./tensors/Tensor_arith_scalar.h ./tensors/Tensor_arith_sub.h ./tensors/Tensor_class.h ./tensors/Tensor_determinant.h ./tensors/Tensor_exp.h ./tensors/Tensor_extract_merge.h ./tensors/Tensor_index.h ./tensors/Tensor_inner.h ./tensors/Tensor_logical.h ./tensors/Tensor_outer.h ./tensors/Tensor_reality.h ./tensors/Tensor_Ta.h ./tensors/Tensor_trace.h ./tensors/Tensor_traits.h ./tensors/Tensor_transpose.h ./tensors/Tensor_unary.h ./Tensors.h ./Threads.h ./Timer.h
|
||||
HFILES=./algorithms/approx/bigfloat.h ./algorithms/approx/bigfloat_double.h ./algorithms/approx/Chebyshev.h ./algorithms/approx/MultiShiftFunction.h ./algorithms/approx/Remez.h ./algorithms/approx/Zolotarev.h ./algorithms/CoarsenedMatrix.h ./algorithms/iterative/AdefGeneric.h ./algorithms/iterative/BfmHDCG.h ./algorithms/iterative/ConjugateGradient.h ./algorithms/iterative/ConjugateGradientMultiShift.h ./algorithms/iterative/ConjugateResidual.h ./algorithms/iterative/GeneralisedMinimumResidual.h ./algorithms/iterative/gmres.h ./algorithms/iterative/NormalEquations.h ./algorithms/iterative/PrecConjugateResidual.h ./algorithms/iterative/PrecGeneralisedConjugateResidual.h ./algorithms/iterative/SchurRedBlack.h ./algorithms/LinearOperator.h ./algorithms/Preconditioner.h ./algorithms/SparseMatrix.h ./Algorithms.h ./AlignedAllocator.h ./cartesian/Cartesian_base.h ./cartesian/Cartesian_full.h ./cartesian/Cartesian_red_black.h ./Cartesian.h ./communicator/Communicator_base.h ./Communicator.h ./Config.h ./cshift/Cshift_common.h ./cshift/Cshift_mpi.h ./cshift/Cshift_none.h ./Cshift.h ./Grid.h ./Init.h ./lattice/Lattice_arith.h ./lattice/Lattice_base.h ./lattice/Lattice_comparison.h ./lattice/Lattice_comparison_utils.h ./lattice/Lattice_conformable.h ./lattice/Lattice_coordinate.h ./lattice/Lattice_ET.h ./lattice/Lattice_local.h ./lattice/Lattice_overload.h ./lattice/Lattice_peekpoke.h ./lattice/Lattice_reality.h ./lattice/Lattice_reduction.h ./lattice/Lattice_rng.h ./lattice/Lattice_trace.h ./lattice/Lattice_transfer.h ./lattice/Lattice_transpose.h ./lattice/Lattice_unary.h ./lattice/Lattice_where.h ./Lattice.h ./Log.h ./MacroMagic.h ./Old/Tensor_peek.h ./Old/Tensor_poke.h ./parallelIO/NerscIO.h ./parameterIO/XMLReader.h ./qcd/action/ActionBase.h ./qcd/action/Actions.h ./qcd/action/DiffAction.h ./qcd/action/fermion/CayleyFermion5D.h ./qcd/action/fermion/ContinuedFractionFermion5D.h ./qcd/action/fermion/DomainWallFermion.h ./qcd/action/fermion/FermionOperator.h ./qcd/action/fermion/g5HermitianLinop.h ./qcd/action/fermion/MobiusFermion.h ./qcd/action/fermion/MobiusZolotarevFermion.h ./qcd/action/fermion/OverlapWilsonCayleyTanhFermion.h ./qcd/action/fermion/OverlapWilsonCayleyZolotarevFermion.h ./qcd/action/fermion/OverlapWilsonContfracTanhFermion.h ./qcd/action/fermion/OverlapWilsonContfracZolotarevFermion.h ./qcd/action/fermion/OverlapWilsonPartialFractionTanhFermion.h ./qcd/action/fermion/OverlapWilsonPartialFractionZolotarevFermion.h ./qcd/action/fermion/PartialFractionFermion5D.h ./qcd/action/fermion/ScaledShamirFermion.h ./qcd/action/fermion/ShamirZolotarevFermion.h ./qcd/action/fermion/WilsonCompressor.h ./qcd/action/fermion/WilsonFermion.h ./qcd/action/fermion/WilsonFermion5D.h ./qcd/action/fermion/WilsonKernels.h ./qcd/action/gauge/WilsonGaugeAction.h ./qcd/action/pseudofermion/TwoFlavour.h ./qcd/hmc/HMC.h ./qcd/hmc/integrators/Integrator.h ./qcd/hmc/integrators/Integrator_algorithm.h ./qcd/hmc/integrators/Integrator_base.h ./qcd/QCD.h ./qcd/spin/Dirac.h ./qcd/spin/TwoSpinor.h ./qcd/utils/CovariantCshift.h ./qcd/utils/LinalgUtils.h ./qcd/utils/SpaceTimeGrid.h ./qcd/utils/SUn.h ./qcd/utils/WilsonLoops.h ./simd/Grid_avx.h ./simd/Grid_avx512.h ./simd/Grid_empty.h ./simd/Grid_neon.h ./simd/Grid_qpx.h ./simd/Grid_sse4.h ./simd/Grid_vector_types.h ./simd/Grid_vector_unops.h ./Simd.h ./stencil/Lebesgue.h ./Stencil.h ./tensors/Tensor_arith.h ./tensors/Tensor_arith_add.h ./tensors/Tensor_arith_mac.h ./tensors/Tensor_arith_mul.h ./tensors/Tensor_arith_scalar.h ./tensors/Tensor_arith_sub.h ./tensors/Tensor_class.h ./tensors/Tensor_determinant.h ./tensors/Tensor_exp.h ./tensors/Tensor_extract_merge.h ./tensors/Tensor_index.h ./tensors/Tensor_inner.h ./tensors/Tensor_logical.h ./tensors/Tensor_outer.h ./tensors/Tensor_reality.h ./tensors/Tensor_Ta.h ./tensors/Tensor_trace.h ./tensors/Tensor_traits.h ./tensors/Tensor_transpose.h ./tensors/Tensor_unary.h ./Tensors.h ./Threads.h ./Timer.h
|
||||
|
||||
CCFILES=./algorithms/approx/MultiShiftFunction.cc ./algorithms/approx/Remez.cc ./algorithms/approx/Zolotarev.cc ./Init.cc ./Log.cc ./qcd/action/fermion/CayleyFermion5D.cc ./qcd/action/fermion/ContinuedFractionFermion5D.cc ./qcd/action/fermion/PartialFractionFermion5D.cc ./qcd/action/fermion/WilsonFermion.cc ./qcd/action/fermion/WilsonFermion5D.cc ./qcd/action/fermion/WilsonKernels.cc ./qcd/action/fermion/WilsonKernelsHand.cc ./qcd/hmc/HMC.cc ./qcd/hmc/integrators/Integrator.cc ./qcd/spin/Dirac.cc ./qcd/utils/SpaceTimeGrid.cc ./stencil/Lebesgue.cc ./stencil/Stencil_common.cc
|
||||
|
||||
@ -13,9 +13,7 @@ namespace Grid {
|
||||
public:
|
||||
RealD Tolerance;
|
||||
Integer MaxIterations;
|
||||
int verbose;
|
||||
ConjugateGradient(RealD tol,Integer maxit) : Tolerance(tol), MaxIterations(maxit) {
|
||||
verbose=0;
|
||||
};
|
||||
|
||||
|
||||
@ -42,14 +40,12 @@ public:
|
||||
cp =a;
|
||||
ssq=norm2(src);
|
||||
|
||||
if ( verbose ) {
|
||||
std::cout<<GridLogMessage <<std::setprecision(4)<< "ConjugateGradient: guess "<<guess<<std::endl;
|
||||
std::cout<<GridLogMessage <<std::setprecision(4)<< "ConjugateGradient: src "<<ssq <<std::endl;
|
||||
std::cout<<GridLogMessage <<std::setprecision(4)<< "ConjugateGradient: mp "<<d <<std::endl;
|
||||
std::cout<<GridLogMessage <<std::setprecision(4)<< "ConjugateGradient: mmp "<<b <<std::endl;
|
||||
std::cout<<GridLogMessage <<std::setprecision(4)<< "ConjugateGradient: cp,r "<<cp <<std::endl;
|
||||
std::cout<<GridLogMessage <<std::setprecision(4)<< "ConjugateGradient: p "<<a <<std::endl;
|
||||
}
|
||||
std::cout<<GridLogIterative <<std::setprecision(4)<< "ConjugateGradient: guess "<<guess<<std::endl;
|
||||
std::cout<<GridLogIterative <<std::setprecision(4)<< "ConjugateGradient: src "<<ssq <<std::endl;
|
||||
std::cout<<GridLogIterative <<std::setprecision(4)<< "ConjugateGradient: mp "<<d <<std::endl;
|
||||
std::cout<<GridLogIterative <<std::setprecision(4)<< "ConjugateGradient: mmp "<<b <<std::endl;
|
||||
std::cout<<GridLogIterative <<std::setprecision(4)<< "ConjugateGradient: cp,r "<<cp <<std::endl;
|
||||
std::cout<<GridLogIterative <<std::setprecision(4)<< "ConjugateGradient: p "<<a <<std::endl;
|
||||
|
||||
RealD rsq = Tolerance* Tolerance*ssq;
|
||||
|
||||
@ -58,7 +54,7 @@ public:
|
||||
return;
|
||||
}
|
||||
|
||||
if(verbose) std::cout<<GridLogMessage << std::setprecision(4)<< "ConjugateGradient: k=0 residual "<<cp<<" rsq"<<rsq<<std::endl;
|
||||
std::cout<<GridLogIterative << std::setprecision(4)<< "ConjugateGradient: k=0 residual "<<cp<<" rsq"<<rsq<<std::endl;
|
||||
|
||||
int k;
|
||||
for (k=1;k<=MaxIterations;k++){
|
||||
@ -80,7 +76,7 @@ public:
|
||||
psi= a*p+psi;
|
||||
p = p*b+r;
|
||||
|
||||
if (verbose) std::cout<<GridLogMessage<<"ConjugateGradient: Iteration " <<k<<" residual "<<cp<< " target"<< rsq<<std::endl;
|
||||
std::cout<<GridLogIterative<<"ConjugateGradient: Iteration " <<k<<" residual "<<cp<< " target"<< rsq<<std::endl;
|
||||
|
||||
// Stopping condition
|
||||
if ( cp <= rsq ) {
|
||||
|
||||
@ -79,4 +79,10 @@
|
||||
///////////////////////////////////////////////////////////////////////////////
|
||||
#include <qcd/action/fermion/g5HermitianLinop.h>
|
||||
|
||||
////////////////////////////////////////
|
||||
// Pseudo fermion combinations
|
||||
////////////////////////////////////////
|
||||
#include <qcd/action/pseudofermion/TwoFlavour.h>
|
||||
|
||||
|
||||
#endif
|
||||
|
||||
@ -56,6 +56,10 @@ namespace Grid {
|
||||
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);
|
||||
|
||||
};
|
||||
|
||||
|
||||
@ -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);
|
||||
|
||||
@ -136,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);
|
||||
|
||||
///////////////////////////////////////////////////////////////
|
||||
|
||||
@ -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);
|
||||
|
||||
@ -94,6 +94,7 @@ namespace Grid {
|
||||
double _M5);
|
||||
|
||||
// DoubleStore
|
||||
virtual void ImportGauge(const LatticeGaugeField &_Umu);
|
||||
void DoubleStore(LatticeDoubledGaugeField &Uds,const LatticeGaugeField &Umu);
|
||||
|
||||
///////////////////////////////////////////////////////////////
|
||||
|
||||
@ -19,8 +19,10 @@ namespace Grid{
|
||||
virtual RealD S(const GaugeField &U) {
|
||||
RealD plaq = WilsonLoops<MatrixField,GaugeField>::avgPlaquette(U);
|
||||
std::cout<<GridLogMessage << "Plaq : "<<plaq << "\n";
|
||||
double vol = U._grid->gSites();
|
||||
return beta*(1.0 -plaq)*(Nd*(Nd-1.0))*vol*0.5;
|
||||
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) {
|
||||
|
||||
|
||||
@ -95,18 +95,18 @@ namespace Grid{
|
||||
////////////////////////////////////////////////////////////////////////
|
||||
// Two flavour pseudofermion action for any dop
|
||||
////////////////////////////////////////////////////////////////////////
|
||||
template<class GaugeField,class MatrixField,class FermionField,class FermionOperator>
|
||||
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;
|
||||
GridBase &Grid;
|
||||
|
||||
FermionField Phi; // the pseudo fermion field for this trajectory
|
||||
|
||||
@ -114,10 +114,11 @@ namespace Grid{
|
||||
/////////////////////////////////////////////////
|
||||
// Pass in required objects.
|
||||
/////////////////////////////////////////////////
|
||||
TwoFlavourPseudoFermionAction(FermionOperator &Op,
|
||||
OperatorFunction<FermionField> & DS,
|
||||
OperatorFunction<FermionField> & AS
|
||||
) : FermOp(Op), DerivativeSolver(DS), ActionSolver(AS) {
|
||||
TwoFlavourPseudoFermionAction(FermionOperator<FermionField,GaugeField> &Op,
|
||||
OperatorFunction<FermionField> & DS,
|
||||
OperatorFunction<FermionField> & AS,
|
||||
GridBase &_Grid
|
||||
) : FermOp(Op), DerivativeSolver(DS), ActionSolver(AS), Phi(&_Grid), Grid(_Grid) {
|
||||
};
|
||||
|
||||
//////////////////////////////////////////////////////////////////////////////////////
|
||||
@ -125,9 +126,28 @@ namespace Grid{
|
||||
//////////////////////////////////////////////////////////////////////////////////////
|
||||
virtual void init(const GaugeField &U, GridParallelRNG& pRNG) {
|
||||
|
||||
// width? Must check
|
||||
gaussian(Phi,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;
|
||||
|
||||
};
|
||||
|
||||
//////////////////////////////////////////////////////
|
||||
@ -135,38 +155,49 @@ namespace Grid{
|
||||
//////////////////////////////////////////////////////
|
||||
virtual RealD S(const GaugeField &U) {
|
||||
|
||||
FermionField X(Grid);
|
||||
FermionField Y(Grid);
|
||||
|
||||
MdagMLinearOperator<FermionOperator<FermionField,GaugeField>,FermionField> MdagMOp(FermOp);
|
||||
FermOp.ImportGauge(U);
|
||||
|
||||
ActionSolver(MdagMop,Phi,X);
|
||||
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) {
|
||||
|
||||
FermionField X(Grid);
|
||||
FermionField Y(Grid);
|
||||
GaugeField tmp(Grid);
|
||||
FermOp.ImportGauge(U);
|
||||
|
||||
MdagMLinearOperator<FermionOperator<FermionField,GaugeField>,FermionField> MdagMOp(FermOp);
|
||||
FermionField X(&Grid);
|
||||
FermionField Y(&Grid);
|
||||
GaugeField tmp(&Grid);
|
||||
|
||||
DerivativeSolver(MdagMop,Phi,X);
|
||||
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=-UdSdU-tmp;
|
||||
FermOp.MDeriv(tmp , X, Y,DaggerYes); dSdU=dSdU+tmp;
|
||||
|
||||
dSdU = Ta(dSdU);
|
||||
|
||||
};
|
||||
|
||||
|
||||
@ -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;
|
||||
|
||||
@ -73,6 +73,7 @@ namespace Grid{
|
||||
|
||||
RealD H1 = MD.S(U); // updated state action
|
||||
std::cout<<GridLogMessage<<"Total H after = "<< H1 << "\n";
|
||||
std::cout<<GridLogMessage<<"DeltaH is "<< H1-H0 << "\n";
|
||||
|
||||
return (H1-H0);
|
||||
}
|
||||
|
||||
@ -63,7 +63,6 @@ namespace Grid{
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
void update_U(LatticeLorentzColourMatrix&U, double ep){
|
||||
//rewrite exponential to deal automatically with the lorentz index?
|
||||
LatticeColourMatrix Umu(U._grid);
|
||||
@ -77,7 +76,6 @@ namespace Grid{
|
||||
|
||||
}
|
||||
|
||||
|
||||
|
||||
friend void IntegratorAlgorithm::step (LatticeLorentzColourMatrix& U,
|
||||
int level, std::vector<int>& clock,
|
||||
@ -92,9 +90,9 @@ namespace Grid{
|
||||
|
||||
~Integrator(){}
|
||||
|
||||
|
||||
//Initialization of momenta and actions
|
||||
void init(LatticeLorentzColourMatrix& U){
|
||||
|
||||
std::cout<<GridLogMessage<< "Integrator init\n";
|
||||
|
||||
MDutils::generate_momenta(*P,pRNG);
|
||||
@ -105,7 +103,6 @@ namespace Grid{
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
// Calculate action
|
||||
RealD S(LatticeLorentzColourMatrix& U){
|
||||
LatticeComplex Hloc(U._grid);
|
||||
@ -119,12 +116,14 @@ namespace Grid{
|
||||
|
||||
RealD H = Hsum.real();
|
||||
|
||||
std::cout<<GridLogMessage << "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.at(level).size(); ++actionID)
|
||||
H += as[level].at(actionID)->S(U);
|
||||
|
||||
std::cout<<GridLogMessage << "Total action H = "<< H << "\n";
|
||||
|
||||
return H;
|
||||
}
|
||||
|
||||
@ -32,8 +32,7 @@ namespace Grid{
|
||||
int fin = Integ->Nrel[0];
|
||||
for(int l=1; l<=level; ++l) fin*= 2.0*Integ->Nrel[l];
|
||||
fin = 3*Integ->Params.MDsteps*fin -1;
|
||||
|
||||
|
||||
|
||||
for(int e=0; e<Integ->Nrel[level]; ++e){
|
||||
|
||||
if(clock[level] == 0){ // initial half step
|
||||
@ -45,7 +44,6 @@ namespace Grid{
|
||||
|
||||
if(level == fl){ // lowest level
|
||||
Integ->update_U(U,0.5*eps);
|
||||
|
||||
for(int l=0; l<level;++l) std::cout<<GridLogMessage<<" ";
|
||||
std::cout<<GridLogMessage<<"U "<< (clock[level]+1) <<std::endl;
|
||||
}else{ // recursive function call
|
||||
@ -81,9 +79,7 @@ namespace Grid{
|
||||
std::cout<<GridLogMessage<<"P "<< clock[level] <<std::endl;
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
|
||||
|
||||
}
|
||||
|
||||
};
|
||||
@ -93,6 +89,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
|
||||
@ -115,6 +112,7 @@ namespace Grid{
|
||||
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<<GridLogMessage<<" ";
|
||||
@ -122,24 +120,21 @@ namespace Grid{
|
||||
}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<<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<<GridLogMessage<<" ";
|
||||
std::cout<<GridLogMessage<<"P "<< 0.5*clock[level] <<std::endl;
|
||||
}
|
||||
|
||||
}
|
||||
|
||||
|
||||
|
||||
|
||||
}
|
||||
};
|
||||
|
||||
|
||||
@ -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;
|
||||
}
|
||||
}
|
||||
|
||||
@ -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;
|
||||
}
|
||||
|
||||
|
||||
@ -7,11 +7,12 @@ namespace Grid {
|
||||
///////////////////////////////////////////////
|
||||
// Ta function for scalar, vector, matrix
|
||||
///////////////////////////////////////////////
|
||||
/*
|
||||
inline ComplexF Ta( const ComplexF &arg){ return arg;}
|
||||
inline ComplexD Ta( const ComplexD &arg){ return arg;}
|
||||
inline RealF Ta( const RealF &arg){ return arg;}
|
||||
inline RealD Ta( const RealD &arg){ return arg;}
|
||||
|
||||
*/
|
||||
|
||||
template<class vtype> inline iScalar<vtype> Ta(const iScalar<vtype>&r)
|
||||
{
|
||||
@ -29,10 +30,11 @@ namespace Grid {
|
||||
}
|
||||
template<class vtype,int N> inline iMatrix<vtype,N> Ta(const iMatrix<vtype,N> &arg)
|
||||
{
|
||||
iMatrix<vtype,N> ret(arg);
|
||||
double factor = (1/(double)N);
|
||||
ret = (ret - adj(arg))*0.5;
|
||||
ret -= trace(ret)*factor;
|
||||
iMatrix<vtype,N> ret;
|
||||
|
||||
double factor = (1.0/(double)N);
|
||||
ret= (arg - adj(arg))*0.5;
|
||||
ret=ret - (trace(ret)*factor);
|
||||
return ret;
|
||||
}
|
||||
|
||||
|
||||
@ -1,5 +1,5 @@
|
||||
|
||||
bin_PROGRAMS = Test_GaugeAction Test_cayley_cg Test_cayley_coarsen_support Test_cayley_even_odd Test_cayley_ldop_cr Test_cf_coarsen_support Test_cf_cr_unprec Test_contfrac_cg Test_contfrac_even_odd Test_cshift Test_cshift_red_black Test_dwf_cg_prec Test_dwf_cg_schur Test_dwf_cg_unprec Test_dwf_cr_unprec Test_dwf_even_odd Test_dwf_fpgcr Test_dwf_hdcr Test_gamma Test_hmc_WilsonGauge Test_lie_generators Test_main Test_multishift_sqrt Test_nersc_io Test_quenched_update Test_remez Test_rng Test_rng_fixed Test_simd Test_stencil Test_wilson_cg_prec Test_wilson_cg_schur Test_wilson_cg_unprec Test_wilson_cr_unprec Test_wilson_even_odd Test_wilson_force Test_wilson_force_phiMdagMphi Test_wilson_force_phiMphi
|
||||
bin_PROGRAMS = Test_GaugeAction Test_cayley_cg Test_cayley_coarsen_support Test_cayley_even_odd Test_cayley_ldop_cr Test_cf_coarsen_support Test_cf_cr_unprec Test_contfrac_cg Test_contfrac_even_odd Test_cshift Test_cshift_red_black Test_dwf_cg_prec Test_dwf_cg_schur Test_dwf_cg_unprec Test_dwf_cr_unprec Test_dwf_even_odd Test_dwf_fpgcr Test_dwf_hdcr Test_gamma Test_hmc_WilsonFermionGauge Test_hmc_WilsonGauge Test_lie_generators Test_main Test_multishift_sqrt Test_nersc_io Test_quenched_update Test_remez Test_rng Test_rng_fixed Test_simd Test_stencil Test_wilson_cg_prec Test_wilson_cg_schur Test_wilson_cg_unprec Test_wilson_cr_unprec Test_wilson_even_odd Test_wilson_force Test_wilson_force_phiMdagMphi Test_wilson_force_phiMphi
|
||||
|
||||
|
||||
Test_GaugeAction_SOURCES=Test_GaugeAction.cc
|
||||
@ -78,6 +78,10 @@ Test_gamma_SOURCES=Test_gamma.cc
|
||||
Test_gamma_LDADD=-lGrid
|
||||
|
||||
|
||||
Test_hmc_WilsonFermionGauge_SOURCES=Test_hmc_WilsonFermionGauge.cc
|
||||
Test_hmc_WilsonFermionGauge_LDADD=-lGrid
|
||||
|
||||
|
||||
Test_hmc_WilsonGauge_SOURCES=Test_hmc_WilsonGauge.cc
|
||||
Test_hmc_WilsonGauge_LDADD=-lGrid
|
||||
|
||||
|
||||
55
tests/Test_hmc_WilsonFermionGauge.cc
Normal file
55
tests/Test_hmc_WilsonFermionGauge.cc
Normal file
@ -0,0 +1,55 @@
|
||||
#include "Grid.h"
|
||||
|
||||
|
||||
using namespace std;
|
||||
using namespace Grid;
|
||||
using namespace Grid::QCD;
|
||||
|
||||
int main (int argc, char ** argv)
|
||||
{
|
||||
Grid_init(&argc,&argv);
|
||||
|
||||
std::vector<int> latt_size = GridDefaultLatt();
|
||||
std::vector<int> simd_layout = GridDefaultSimd(4,vComplex::Nsimd());
|
||||
std::vector<int> mpi_layout = GridDefaultMpi();
|
||||
|
||||
GridCartesian Fine(latt_size,simd_layout,mpi_layout);
|
||||
GridRedBlackCartesian RBFine(latt_size,simd_layout,mpi_layout);
|
||||
GridParallelRNG pRNG(&Fine);
|
||||
pRNG.SeedRandomDevice();
|
||||
LatticeLorentzColourMatrix U(&Fine);
|
||||
|
||||
SU3::HotConfiguration(pRNG, U);
|
||||
|
||||
// simplify template declaration? Strip the lorentz from the second template
|
||||
WilsonGaugeAction<LatticeLorentzColourMatrix, LatticeColourMatrix> Waction(5.6);
|
||||
|
||||
Real mass=0.01;
|
||||
WilsonFermion FermOp(U,Fine,RBFine,mass);
|
||||
|
||||
ConjugateGradient<LatticeFermion> CG(1.0e-8,10000);
|
||||
|
||||
TwoFlavourPseudoFermionAction<LatticeLorentzColourMatrix, LatticeColourMatrix,LatticeFermion>
|
||||
Pseudofermion(FermOp,CG,CG,Fine);
|
||||
|
||||
|
||||
//Collect actions
|
||||
ActionLevel Level1;
|
||||
Level1.push_back(&Waction);
|
||||
Level1.push_back(&Pseudofermion);
|
||||
ActionSet FullSet;
|
||||
FullSet.push_back(Level1);
|
||||
|
||||
// Create integrator
|
||||
typedef LeapFrog IntegratorAlgorithm;// change here to change the algorithm
|
||||
IntegratorParameters MDpar(12,40,1.0);
|
||||
std::vector<int> rel ={1};
|
||||
Integrator<IntegratorAlgorithm> MDynamics(&Fine,MDpar, FullSet,rel);
|
||||
|
||||
// Create HMC
|
||||
HMCparameters HMCpar;
|
||||
HybridMonteCarlo<IntegratorAlgorithm> HMC(HMCpar, MDynamics);
|
||||
|
||||
HMC.evolve(U);
|
||||
|
||||
}
|
||||
@ -48,17 +48,19 @@ int main (int argc, char ** argv)
|
||||
LatticeGaugeField tmp(&Grid);
|
||||
|
||||
Dw.MDeriv(tmp , Mphi, phi,DaggerNo ); UdSdU=tmp;
|
||||
Dw.MDeriv(tmp , phi, Mphi,DaggerYes ); UdSdU=UdSdU+tmp;
|
||||
|
||||
|
||||
Dw.MDeriv(tmp , phi, Mphi,DaggerYes ); UdSdU=(UdSdU+tmp);
|
||||
|
||||
LatticeFermion Ftmp (&Grid);
|
||||
|
||||
////////////////////////////////////
|
||||
// Modify the gauge field a little
|
||||
////////////////////////////////////
|
||||
RealD dt = 1.0e-6;
|
||||
|
||||
RealD dt = 0.0001;
|
||||
RealD Hmom = 0.0;
|
||||
RealD Hmomprime = 0.0;
|
||||
RealD Hmompp = 0.0;
|
||||
LatticeColourMatrix mommu(&Grid);
|
||||
LatticeColourMatrix forcemu(&Grid);
|
||||
LatticeGaugeField mom(&Grid);
|
||||
LatticeGaugeField Uprime(&Grid);
|
||||
|
||||
@ -66,13 +68,26 @@ int main (int argc, char ** argv)
|
||||
|
||||
SU3::GaussianLieAlgebraMatrix(pRNG, mommu); // Traceless antihermitian momentum; gaussian in lie alg
|
||||
|
||||
Hmom -= real(sum(trace(mommu*mommu)));
|
||||
|
||||
PokeIndex<LorentzIndex>(mom,mommu,mu);
|
||||
|
||||
// fourth order exponential approx
|
||||
parallel_for(auto i=mom.begin();i<mom.end();i++){
|
||||
Uprime[i](mu) =U[i](mu)+ mom[i](mu)*U[i](mu)*dt;
|
||||
Uprime[i](mu) =
|
||||
U[i](mu)
|
||||
+ mom[i](mu)*U[i](mu)*dt
|
||||
+ mom[i](mu) *mom[i](mu) *U[i](mu)*(dt*dt/2.0)
|
||||
+ mom[i](mu) *mom[i](mu) *mom[i](mu) *U[i](mu)*(dt*dt*dt/6.0)
|
||||
+ mom[i](mu) *mom[i](mu) *mom[i](mu) *mom[i](mu) *U[i](mu)*(dt*dt*dt*dt/24.0)
|
||||
+ mom[i](mu) *mom[i](mu) *mom[i](mu) *mom[i](mu) *mom[i](mu) *U[i](mu)*(dt*dt*dt*dt*dt/120.0)
|
||||
+ mom[i](mu) *mom[i](mu) *mom[i](mu) *mom[i](mu) *mom[i](mu) *mom[i](mu) *U[i](mu)*(dt*dt*dt*dt*dt*dt/720.0)
|
||||
;
|
||||
}
|
||||
|
||||
}
|
||||
|
||||
std::cout << GridLogMessage <<"Initial mom hamiltonian is "<< Hmom <<std::endl;
|
||||
Dw.DoubleStore(Dw.Umu,Uprime);
|
||||
Dw.M (phi,MphiPrime);
|
||||
|
||||
@ -81,21 +96,75 @@ int main (int argc, char ** argv)
|
||||
//////////////////////////////////////////////
|
||||
// Use derivative to estimate dS
|
||||
//////////////////////////////////////////////
|
||||
LatticeComplex dS(&Grid); dS = zero;
|
||||
|
||||
parallel_for(auto i=mom.begin();i<mom.end();i++){
|
||||
for(int mu=0;mu<Nd;mu++){
|
||||
// dS[i]() = dS[i]()+trace(mom[i](mu) * UdSdU[i](mu) - mom[i](mu)* adj( UdSdU[i](mu)) )*dt;
|
||||
dS[i]() = dS[i]()+trace(mom[i](mu) * (UdSdU[i](mu)))*dt;
|
||||
dS[i]() = dS[i]()-trace(mom[i](mu) * adj(UdSdU[i](mu)))*dt;
|
||||
}
|
||||
|
||||
for(int mu=0;mu<Nd;mu++){
|
||||
std::cout << "" <<std::endl;
|
||||
mommu = PeekIndex<LorentzIndex>(mom,mu);
|
||||
std::cout << GridLogMessage<< " Mommu " << norm2(mommu)<<std::endl;
|
||||
mommu = mommu+adj(mommu);
|
||||
std::cout << GridLogMessage<< " Mommu + Mommudag " << norm2(mommu)<<std::endl;
|
||||
mommu = PeekIndex<LorentzIndex>(UdSdU,mu);
|
||||
std::cout << GridLogMessage<< " dsdumu " << norm2(mommu)<<std::endl;
|
||||
mommu = mommu+adj(mommu);
|
||||
std::cout << GridLogMessage<< " dsdumu + dag " << norm2(mommu)<<std::endl;
|
||||
}
|
||||
|
||||
LatticeComplex dS(&Grid); dS = zero;
|
||||
LatticeComplex dSmom(&Grid); dSmom = zero;
|
||||
LatticeComplex dSmom2(&Grid); dSmom2 = zero;
|
||||
for(int mu=0;mu<Nd;mu++){
|
||||
mommu = PeekIndex<LorentzIndex>(UdSdU,mu);
|
||||
mommu=Ta(mommu)*2.0;
|
||||
PokeIndex<LorentzIndex>(UdSdU,mommu,mu);
|
||||
}
|
||||
|
||||
for(int mu=0;mu<Nd;mu++){
|
||||
mommu = PeekIndex<LorentzIndex>(mom,mu);
|
||||
std::cout << GridLogMessage<< " Mommu " << norm2(mommu)<<std::endl;
|
||||
mommu = mommu+adj(mommu);
|
||||
std::cout << GridLogMessage<< " Mommu + Mommudag " << norm2(mommu)<<std::endl;
|
||||
mommu = PeekIndex<LorentzIndex>(UdSdU,mu);
|
||||
std::cout << GridLogMessage<< " dsdumu " << norm2(mommu)<<std::endl;
|
||||
mommu = mommu+adj(mommu);
|
||||
std::cout << GridLogMessage<< " dsdumu + dag " << norm2(mommu)<<std::endl;
|
||||
}
|
||||
|
||||
for(int mu=0;mu<Nd;mu++){
|
||||
forcemu = PeekIndex<LorentzIndex>(UdSdU,mu);
|
||||
mommu = PeekIndex<LorentzIndex>(mom,mu);
|
||||
|
||||
// Update PF action density
|
||||
dS = dS+trace(mommu*forcemu)*dt;
|
||||
|
||||
dSmom = dSmom - trace(mommu*forcemu) * dt;
|
||||
dSmom2 = dSmom2 - trace(forcemu*forcemu) *(0.25* dt*dt);
|
||||
|
||||
// Update mom action density
|
||||
mommu = mommu + forcemu*(dt*0.5);
|
||||
|
||||
Hmomprime -= real(sum(trace(mommu*mommu)));
|
||||
|
||||
}
|
||||
|
||||
Complex dSpred = sum(dS);
|
||||
Complex dSm = sum(dSmom);
|
||||
Complex dSm2 = sum(dSmom2);
|
||||
|
||||
|
||||
std::cout << GridLogMessage <<"Initial mom hamiltonian is "<< Hmom <<std::endl;
|
||||
std::cout << GridLogMessage <<"Final mom hamiltonian is "<< Hmomprime <<std::endl;
|
||||
std::cout << GridLogMessage <<"Delta mom hamiltonian is "<< Hmomprime-Hmom <<std::endl;
|
||||
|
||||
std::cout << GridLogMessage << " S "<<S<<std::endl;
|
||||
std::cout << GridLogMessage << " Sprime "<<Sprime<<std::endl;
|
||||
std::cout << GridLogMessage << "dS "<<Sprime-S<<std::endl;
|
||||
std::cout << GridLogMessage << "predict dS "<< dSpred <<std::endl;
|
||||
std::cout << GridLogMessage <<"dSm "<< dSm<<std::endl;
|
||||
std::cout << GridLogMessage <<"dSm2"<< dSm2<<std::endl;
|
||||
|
||||
std::cout << GridLogMessage << "Total dS "<< Hmomprime - Hmom + Sprime - S <<std::endl;
|
||||
|
||||
|
||||
std::cout<< GridLogMessage << "Done" <<std::endl;
|
||||
Grid_finalize();
|
||||
|
||||
Loading…
Reference in New Issue
Block a user