diff --git a/lib/Make.inc b/lib/Make.inc index 92699dbe..409bf754 100644 --- a/lib/Make.inc +++ b/lib/Make.inc @@ -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 diff --git a/lib/algorithms/iterative/ConjugateGradient.h b/lib/algorithms/iterative/ConjugateGradient.h index 70b89f90..b919dd91 100644 --- a/lib/algorithms/iterative/ConjugateGradient.h +++ b/lib/algorithms/iterative/ConjugateGradient.h @@ -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< +//////////////////////////////////////// +// Pseudo fermion combinations +//////////////////////////////////////// +#include + + #endif diff --git a/lib/qcd/action/fermion/FermionOperator.h b/lib/qcd/action/fermion/FermionOperator.h index a8f10380..b6ba2d08 100644 --- a/lib/qcd/action/fermion/FermionOperator.h +++ b/lib/qcd/action/fermion/FermionOperator.h @@ -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); }; diff --git a/lib/qcd/action/fermion/WilsonFermion.cc b/lib/qcd/action/fermion/WilsonFermion.cc index afd7edd7..2c2132d1 100644 --- a/lib/qcd/action/fermion/WilsonFermion.cc +++ b/lib/qcd/action/fermion/WilsonFermion.cc @@ -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); diff --git a/lib/qcd/action/fermion/WilsonFermion.h b/lib/qcd/action/fermion/WilsonFermion.h index 5422ca52..79d3e626 100644 --- a/lib/qcd/action/fermion/WilsonFermion.h +++ b/lib/qcd/action/fermion/WilsonFermion.h @@ -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); /////////////////////////////////////////////////////////////// diff --git a/lib/qcd/action/fermion/WilsonFermion5D.cc b/lib/qcd/action/fermion/WilsonFermion5D.cc index 0387a26e..c3329503 100644 --- a/lib/qcd/action/fermion/WilsonFermion5D.cc +++ b/lib/qcd/action/fermion/WilsonFermion5D.cc @@ -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); diff --git a/lib/qcd/action/fermion/WilsonFermion5D.h b/lib/qcd/action/fermion/WilsonFermion5D.h index 7ff5d1c3..d5cff3e6 100644 --- a/lib/qcd/action/fermion/WilsonFermion5D.h +++ b/lib/qcd/action/fermion/WilsonFermion5D.h @@ -94,6 +94,7 @@ namespace Grid { double _M5); // DoubleStore + virtual void ImportGauge(const LatticeGaugeField &_Umu); void DoubleStore(LatticeDoubledGaugeField &Uds,const LatticeGaugeField &Umu); /////////////////////////////////////////////////////////////// diff --git a/lib/qcd/action/gauge/WilsonGaugeAction.h b/lib/qcd/action/gauge/WilsonGaugeAction.h index fc31a5c3..f5f46070 100644 --- a/lib/qcd/action/gauge/WilsonGaugeAction.h +++ b/lib/qcd/action/gauge/WilsonGaugeAction.h @@ -19,8 +19,10 @@ namespace Grid{ virtual RealD S(const GaugeField &U) { RealD plaq = WilsonLoops::avgPlaquette(U); std::cout<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 "< + template class TwoFlavourPseudoFermionAction : public Action { private: - + FermionOperator & FermOp;// the basic operator OperatorFunction &DerivativeSolver; OperatorFunction &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 & DS, - OperatorFunction & AS - ) : FermOp(Op), DerivativeSolver(DS), ActionSolver(AS) { + TwoFlavourPseudoFermionAction(FermionOperator &Op, + OperatorFunction & DS, + OperatorFunction & 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,FermionField> MdagMOp(FermOp); + FermOp.ImportGauge(U); - ActionSolver(MdagMop,Phi,X); + FermionField X(&Grid); + FermionField Y(&Grid); + + MdagMLinearOperator ,FermionField> MdagMOp(FermOp); + X=zero; + ActionSolver(MdagMOp,Phi,X); MdagMOp.Op(X,Y); RealD action = norm2(Y); - + std::cout << GridLogMessage << "Pseudofermion action "<,FermionField> MdagMOp(FermOp); + FermionField X(&Grid); + FermionField Y(&Grid); + GaugeField tmp(&Grid); - DerivativeSolver(MdagMop,Phi,X); + MdagMLinearOperator ,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); }; diff --git a/lib/qcd/hmc/HMC.cc b/lib/qcd/hmc/HMC.cc index a6130224..f7a5d0b8 100644 --- a/lib/qcd/hmc/HMC.cc +++ b/lib/qcd/hmc/HMC.cc @@ -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; diff --git a/lib/qcd/hmc/HMC.h b/lib/qcd/hmc/HMC.h index 10947719..bc4f47ff 100644 --- a/lib/qcd/hmc/HMC.h +++ b/lib/qcd/hmc/HMC.h @@ -73,6 +73,7 @@ namespace Grid{ RealD H1 = MD.S(U); // updated state action std::cout<& clock, @@ -92,9 +90,9 @@ namespace Grid{ ~Integrator(){} - //Initialization of momenta and actions void init(LatticeLorentzColourMatrix& U){ + std::cout<S(U); + + std::cout<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; eNrel[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& clock, Integrator* Integ){ + // level : current level // fl : final level // eps : current step size @@ -115,6 +112,7 @@ namespace Grid{ for(int l=0; lupdate_U(U, eps); for(int l=0; lupdate_P(U, level,eps/2.0); - ++clock[level]; for(int l=0; lupdate_P(U, level,eps); - clock[level]+=2; for(int l=0; lderiv(U,force); + + Complex dSdt=0.0; + for(int mu=0;mu(force,mu); + mommu=PeekIndex(*P,mu); + + dSdt += sum(trace(forcemu*(*P))); + + } + std::cout << GridLogMessage << " action "< inline iScalar Ta(const iScalar&r) { @@ -29,10 +30,11 @@ namespace Grid { } template inline iMatrix Ta(const iMatrix &arg) { - iMatrix ret(arg); - double factor = (1/(double)N); - ret = (ret - adj(arg))*0.5; - ret -= trace(ret)*factor; + iMatrix ret; + + double factor = (1.0/(double)N); + ret= (arg - adj(arg))*0.5; + ret=ret - (trace(ret)*factor); return ret; } diff --git a/tests/Make.inc b/tests/Make.inc index 8b313ccc..38714c2f 100644 --- a/tests/Make.inc +++ b/tests/Make.inc @@ -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 diff --git a/tests/Test_hmc_WilsonFermionGauge.cc b/tests/Test_hmc_WilsonFermionGauge.cc new file mode 100644 index 00000000..3a10b191 --- /dev/null +++ b/tests/Test_hmc_WilsonFermionGauge.cc @@ -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 latt_size = GridDefaultLatt(); + std::vector simd_layout = GridDefaultSimd(4,vComplex::Nsimd()); + std::vector 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 Waction(5.6); + + Real mass=0.01; + WilsonFermion FermOp(U,Fine,RBFine,mass); + + ConjugateGradient CG(1.0e-8,10000); + + TwoFlavourPseudoFermionAction + 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 rel ={1}; + Integrator MDynamics(&Fine,MDpar, FullSet,rel); + + // Create HMC + HMCparameters HMCpar; + HybridMonteCarlo HMC(HMCpar, MDynamics); + + HMC.evolve(U); + +} diff --git a/tests/Test_wilson_force.cc b/tests/Test_wilson_force.cc index 6a8d7611..b08ccd45 100644 --- a/tests/Test_wilson_force.cc +++ b/tests/Test_wilson_force.cc @@ -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(mom,mommu,mu); + + // fourth order exponential approx parallel_for(auto i=mom.begin();i(mom,mu); + std::cout << GridLogMessage<< " Mommu " << norm2(mommu)<(UdSdU,mu); + std::cout << GridLogMessage<< " dsdumu " << norm2(mommu)<(UdSdU,mu); + mommu=Ta(mommu)*2.0; + PokeIndex(UdSdU,mommu,mu); + } + + for(int mu=0;mu(mom,mu); + std::cout << GridLogMessage<< " Mommu " << norm2(mommu)<(UdSdU,mu); + std::cout << GridLogMessage<< " dsdumu " << norm2(mommu)<(UdSdU,mu); + mommu = PeekIndex(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 <