mirror of
https://github.com/paboyle/Grid.git
synced 2024-11-10 07:55:35 +00:00
Rearranging files in hmc
This commit is contained in:
neo
parent
68fe0769a1
commit
fa42b652e5
@ -1,4 +1,4 @@
|
||||
|
||||
HFILES=./Cshift.h ./simd/Grid_avx.h ./simd/Grid_vector_types.h ./simd/Grid_sse4.h ./simd/Grid_avx512.h ./simd/Grid_empty.h ./simd/Grid_vector_unops.h ./simd/Grid_neon.h ./simd/Grid_qpx.h ./Tensors.h ./Algorithms.h ./communicator/Communicator_base.h ./lattice/Lattice_rng.h ./lattice/Lattice_reduction.h ./lattice/Lattice_transfer.h ./lattice/Lattice_unary.h ./lattice/Lattice_peekpoke.h ./lattice/Lattice_coordinate.h ./lattice/Lattice_comparison.h ./lattice/Lattice_overload.h ./lattice/Lattice_reality.h ./lattice/Lattice_local.h ./lattice/Lattice_conformable.h ./lattice/Lattice_where.h ./lattice/Lattice_comparison_utils.h ./lattice/Lattice_arith.h ./lattice/Lattice_base.h ./lattice/Lattice_ET.h ./lattice/Lattice_transpose.h ./lattice/Lattice_trace.h ./Stencil.h ./tensors/Tensor_arith_sub.h ./tensors/Tensor_exp.h ./tensors/Tensor_arith_mul.h ./tensors/Tensor_class.h ./tensors/Tensor_logical.h ./tensors/Tensor_transpose.h ./tensors/Tensor_arith_mac.h ./tensors/Tensor_arith_scalar.h ./tensors/Tensor_reality.h ./tensors/Tensor_trace.h ./tensors/Tensor_index.h ./tensors/Tensor_arith_add.h ./tensors/Tensor_outer.h ./tensors/Tensor_inner.h ./tensors/Tensor_traits.h ./tensors/Tensor_Ta.h ./tensors/Tensor_unary.h ./tensors/Tensor_determinant.h ./tensors/Tensor_arith.h ./tensors/Tensor_extract_merge.h ./Communicator.h ./Cartesian.h ./parallelIO/NerscIO.h ./qcd/QCD.h ./qcd/hmc/integrators/Integrator_base.h ./qcd/hmc/HMC.h ./qcd/utils/SpaceTimeGrid.h ./qcd/utils/SUn.h ./qcd/utils/LinalgUtils.h ./qcd/utils/CovariantCshift.h ./qcd/utils/WilsonLoops.h ./qcd/action/ActionBase.h ./qcd/action/gauge/WilsonGaugeAction.h ./qcd/action/Actions.h ./qcd/action/fermion/CayleyFermion5D.h ./qcd/action/fermion/ScaledShamirFermion.h ./qcd/action/fermion/MobiusFermion.h ./qcd/action/fermion/OverlapWilsonContfracTanhFermion.h ./qcd/action/fermion/PartialFractionFermion5D.h ./qcd/action/fermion/ShamirZolotarevFermion.h ./qcd/action/fermion/FermionOperator.h ./qcd/action/fermion/WilsonFermion5D.h ./qcd/action/fermion/WilsonCompressor.h ./qcd/action/fermion/OverlapWilsonPartialFractionZolotarevFermion.h ./qcd/action/fermion/WilsonKernels.h ./qcd/action/fermion/DomainWallFermion.h ./qcd/action/fermion/OverlapWilsonPartialFractionTanhFermion.h ./qcd/action/fermion/OverlapWilsonContfracZolotarevFermion.h ./qcd/action/fermion/MobiusZolotarevFermion.h ./qcd/action/fermion/g5HermitianLinop.h ./qcd/action/fermion/OverlapWilsonCayleyTanhFermion.h ./qcd/action/fermion/WilsonFermion.h ./qcd/action/fermion/ContinuedFractionFermion5D.h ./qcd/action/fermion/OverlapWilsonCayleyZolotarevFermion.h ./qcd/spin/TwoSpinor.h ./qcd/spin/Dirac.h ./cshift/Cshift_common.h ./cshift/Cshift_none.h ./cshift/Cshift_mpi.h ./Simd.h ./GridConfig.h ./cartesian/Cartesian_base.h ./cartesian/Cartesian_red_black.h ./cartesian/Cartesian_full.h ./AlignedAllocator.h ./Lattice.h ./Old/Tensor_poke.h ./Old/Tensor_peek.h ./Threads.h ./Grid.h ./algorithms/Preconditioner.h ./algorithms/iterative/ConjugateResidual.h ./algorithms/iterative/PrecGeneralisedConjugateResidual.h ./algorithms/iterative/ConjugateGradientMultiShift.h ./algorithms/iterative/SchurRedBlack.h ./algorithms/iterative/NormalEquations.h ./algorithms/iterative/ConjugateGradient.h ./algorithms/iterative/AdefGeneric.h ./algorithms/approx/Chebyshev.h ./algorithms/approx/Zolotarev.h ./algorithms/approx/MultiShiftFunction.h ./algorithms/approx/bigfloat.h ./algorithms/approx/bigfloat_double.h ./algorithms/approx/Remez.h ./algorithms/LinearOperator.h ./algorithms/SparseMatrix.h ./algorithms/CoarsenedMatrix.h ./stencil/Lebesgue.h
|
||||
HFILES=./Cshift.h ./simd/Grid_avx.h ./simd/Grid_vector_types.h ./simd/Grid_sse4.h ./simd/Grid_avx512.h ./simd/Grid_empty.h ./simd/Grid_vector_unops.h ./simd/Grid_neon.h ./simd/Grid_qpx.h ./Tensors.h ./Algorithms.h ./communicator/Communicator_base.h ./lattice/Lattice_rng.h ./lattice/Lattice_reduction.h ./lattice/Lattice_transfer.h ./lattice/Lattice_unary.h ./lattice/Lattice_peekpoke.h ./lattice/Lattice_coordinate.h ./lattice/Lattice_comparison.h ./lattice/Lattice_overload.h ./lattice/Lattice_reality.h ./lattice/Lattice_local.h ./lattice/Lattice_conformable.h ./lattice/Lattice_where.h ./lattice/Lattice_comparison_utils.h ./lattice/Lattice_arith.h ./lattice/Lattice_base.h ./lattice/Lattice_ET.h ./lattice/Lattice_transpose.h ./lattice/Lattice_trace.h ./Stencil.h ./tensors/Tensor_arith_sub.h ./tensors/Tensor_exp.h ./tensors/Tensor_arith_mul.h ./tensors/Tensor_class.h ./tensors/Tensor_logical.h ./tensors/Tensor_transpose.h ./tensors/Tensor_arith_mac.h ./tensors/Tensor_arith_scalar.h ./tensors/Tensor_reality.h ./tensors/Tensor_trace.h ./tensors/Tensor_index.h ./tensors/Tensor_arith_add.h ./tensors/Tensor_outer.h ./tensors/Tensor_inner.h ./tensors/Tensor_traits.h ./tensors/Tensor_Ta.h ./tensors/Tensor_unary.h ./tensors/Tensor_determinant.h ./tensors/Tensor_arith.h ./tensors/Tensor_extract_merge.h ./Communicator.h ./Cartesian.h ./parallelIO/NerscIO.h ./qcd/QCD.h ./qcd/hmc/integrators/Integrator.h ./qcd/hmc/integrators/Integrator_algorithm.h ./qcd/hmc/HMC.h ./qcd/utils/SpaceTimeGrid.h ./qcd/utils/SUn.h ./qcd/utils/LinalgUtils.h ./qcd/utils/CovariantCshift.h ./qcd/utils/WilsonLoops.h ./qcd/action/ActionBase.h ./qcd/action/gauge/WilsonGaugeAction.h ./qcd/action/Actions.h ./qcd/action/fermion/CayleyFermion5D.h ./qcd/action/fermion/ScaledShamirFermion.h ./qcd/action/fermion/MobiusFermion.h ./qcd/action/fermion/OverlapWilsonContfracTanhFermion.h ./qcd/action/fermion/PartialFractionFermion5D.h ./qcd/action/fermion/ShamirZolotarevFermion.h ./qcd/action/fermion/FermionOperator.h ./qcd/action/fermion/WilsonFermion5D.h ./qcd/action/fermion/WilsonCompressor.h ./qcd/action/fermion/OverlapWilsonPartialFractionZolotarevFermion.h ./qcd/action/fermion/WilsonKernels.h ./qcd/action/fermion/DomainWallFermion.h ./qcd/action/fermion/OverlapWilsonPartialFractionTanhFermion.h ./qcd/action/fermion/OverlapWilsonContfracZolotarevFermion.h ./qcd/action/fermion/MobiusZolotarevFermion.h ./qcd/action/fermion/g5HermitianLinop.h ./qcd/action/fermion/OverlapWilsonCayleyTanhFermion.h ./qcd/action/fermion/WilsonFermion.h ./qcd/action/fermion/ContinuedFractionFermion5D.h ./qcd/action/fermion/OverlapWilsonCayleyZolotarevFermion.h ./qcd/spin/TwoSpinor.h ./qcd/spin/Dirac.h ./cshift/Cshift_common.h ./cshift/Cshift_none.h ./cshift/Cshift_mpi.h ./Simd.h ./GridConfig.h ./cartesian/Cartesian_base.h ./cartesian/Cartesian_red_black.h ./cartesian/Cartesian_full.h ./AlignedAllocator.h ./Lattice.h ./Old/Tensor_poke.h ./Old/Tensor_peek.h ./Threads.h ./Grid.h ./algorithms/Preconditioner.h ./algorithms/iterative/ConjugateResidual.h ./algorithms/iterative/PrecGeneralisedConjugateResidual.h ./algorithms/iterative/ConjugateGradientMultiShift.h ./algorithms/iterative/SchurRedBlack.h ./algorithms/iterative/NormalEquations.h ./algorithms/iterative/ConjugateGradient.h ./algorithms/iterative/AdefGeneric.h ./algorithms/approx/Chebyshev.h ./algorithms/approx/Zolotarev.h ./algorithms/approx/MultiShiftFunction.h ./algorithms/approx/bigfloat.h ./algorithms/approx/bigfloat_double.h ./algorithms/approx/Remez.h ./algorithms/LinearOperator.h ./algorithms/SparseMatrix.h ./algorithms/CoarsenedMatrix.h ./stencil/Lebesgue.h
|
||||
|
||||
CCFILES=./qcd/hmc/integrators/Integrator_base.cc ./qcd/hmc/HMC.cc ./qcd/utils/SpaceTimeGrid.cc ./qcd/action/fermion/WilsonKernels.cc ./qcd/action/fermion/PartialFractionFermion5D.cc ./qcd/action/fermion/CayleyFermion5D.cc ./qcd/action/fermion/WilsonKernelsHand.cc ./qcd/action/fermion/WilsonFermion.cc ./qcd/action/fermion/ContinuedFractionFermion5D.cc ./qcd/action/fermion/WilsonFermion5D.cc ./qcd/spin/Dirac.cc ./GridInit.cc ./algorithms/approx/MultiShiftFunction.cc ./algorithms/approx/Remez.cc ./algorithms/approx/Zolotarev.cc ./stencil/Lebesgue.cc ./stencil/Stencil_common.cc
|
||||
CCFILES=./qcd/hmc/integrators/Integrator.cc ./qcd/hmc/HMC.cc ./qcd/utils/SpaceTimeGrid.cc ./qcd/action/fermion/WilsonKernels.cc ./qcd/action/fermion/PartialFractionFermion5D.cc ./qcd/action/fermion/CayleyFermion5D.cc ./qcd/action/fermion/WilsonKernelsHand.cc ./qcd/action/fermion/WilsonFermion.cc ./qcd/action/fermion/ContinuedFractionFermion5D.cc ./qcd/action/fermion/WilsonFermion5D.cc ./qcd/spin/Dirac.cc ./GridInit.cc ./algorithms/approx/MultiShiftFunction.cc ./algorithms/approx/Remez.cc ./algorithms/approx/Zolotarev.cc ./stencil/Lebesgue.cc ./stencil/Stencil_common.cc
|
||||
|
||||
@ -413,7 +413,8 @@ namespace QCD {
|
||||
#include <qcd/utils/WilsonLoops.h>
|
||||
#include <qcd/utils/SUn.h>
|
||||
#include <qcd/action/Actions.h>
|
||||
#include <qcd/hmc/integrators/Integrator_base.h>
|
||||
#include <qcd/hmc/integrators/Integrator.h>
|
||||
#include <qcd/hmc/integrators/Integrator_algorithm.h>
|
||||
#include <qcd/hmc/HMC.h>
|
||||
|
||||
|
||||
|
||||
@ -25,11 +25,11 @@ namespace Grid{
|
||||
HMCparameters();
|
||||
};
|
||||
|
||||
template <class IntegType>
|
||||
template <class Algorithm>
|
||||
class HybridMonteCarlo{
|
||||
const HMCparameters Params;
|
||||
GridSerialRNG sRNG;
|
||||
Integrator<IntegType>& MD;
|
||||
Integrator<Algorithm>& MD;
|
||||
|
||||
bool metropolis_test(const RealD DeltaH){
|
||||
RealD rn_test;
|
||||
@ -51,23 +51,21 @@ namespace Grid{
|
||||
|
||||
RealD evolve_step(LatticeLorentzColourMatrix& U){
|
||||
|
||||
MD.init(U); // set U and initialize P and phi's
|
||||
RealD H0 = MD.S(U); // current state
|
||||
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";
|
||||
|
||||
MD.integrate(U,0);
|
||||
MD.integrate(U);
|
||||
|
||||
RealD H1 = MD.S(U); // updated state
|
||||
RealD H1 = MD.S(U); // updated state action
|
||||
std::cout<<"Total H after = "<< H1 << "\n";
|
||||
|
||||
return (H1-H0);
|
||||
}
|
||||
|
||||
|
||||
|
||||
public:
|
||||
HybridMonteCarlo(HMCparameters Pms,
|
||||
Integrator<IntegType>& MolDyn):
|
||||
Integrator<Algorithm>& MolDyn):
|
||||
Params(Pms),MD(MolDyn){
|
||||
//FIXME
|
||||
|
||||
@ -87,7 +85,6 @@ namespace Grid{
|
||||
|
||||
DeltaH = evolve_step(Uin);
|
||||
std::cout<< " dH = "<< DeltaH << "\n";
|
||||
|
||||
}
|
||||
|
||||
// Actual updates (evolve a copy Ucopy then copy back eventually)
|
||||
|
||||
148
lib/qcd/hmc/integrators/Integrator.h
Normal file
148
lib/qcd/hmc/integrators/Integrator.h
Normal file
@ -0,0 +1,148 @@
|
||||
//--------------------------------------------------------------------
|
||||
/*! @file Integrator.h
|
||||
* @brief Declaration of classes for the Molecular Dynamics integrator
|
||||
*
|
||||
* @author Guido Cossu
|
||||
*/
|
||||
//--------------------------------------------------------------------
|
||||
|
||||
#ifndef INTEGRATOR_INCLUDED
|
||||
#define INTEGRATOR_INCLUDED
|
||||
|
||||
class Observer;
|
||||
|
||||
#include <memory>
|
||||
|
||||
namespace Grid{
|
||||
namespace QCD{
|
||||
|
||||
typedef Action<LatticeLorentzColourMatrix>* ActPtr; // now force the same colours as the rest of the code
|
||||
typedef std::vector<ActPtr> ActionLevel;
|
||||
typedef std::vector<ActionLevel> ActionSet;
|
||||
typedef std::vector<Observer*> ObserverList;
|
||||
|
||||
struct IntegratorParameters{
|
||||
int Nexp;
|
||||
int MDsteps; // number of outer steps
|
||||
RealD trajL; // trajectory length
|
||||
RealD stepsize;
|
||||
|
||||
IntegratorParameters(int Nexp_,
|
||||
int MDsteps_,
|
||||
RealD trajL_):
|
||||
Nexp(Nexp_),MDsteps(MDsteps_),trajL(trajL_),stepsize(trajL/MDsteps){};
|
||||
};
|
||||
|
||||
|
||||
namespace MDutils{
|
||||
void generate_momenta(LatticeLorentzColourMatrix&,GridParallelRNG&);
|
||||
void generate_momenta_su3(LatticeLorentzColourMatrix&,GridParallelRNG&);
|
||||
}
|
||||
|
||||
/*! @brief Class for Molecular Dynamics management */
|
||||
template< class IntegratorAlgorithm >
|
||||
class Integrator{
|
||||
private:
|
||||
IntegratorParameters Params;
|
||||
const ActionSet as;
|
||||
const std::vector<int> Nrel; //relative step size per level
|
||||
std::unique_ptr<LatticeLorentzColourMatrix> P;
|
||||
GridParallelRNG pRNG;
|
||||
//ObserverList observers; // not yet
|
||||
|
||||
IntegratorAlgorithm TheIntegrator;
|
||||
|
||||
void register_observers();
|
||||
void notify_observers();
|
||||
|
||||
void update_P(LatticeLorentzColourMatrix&U, int level,double ep){
|
||||
for(int a=0; a<as[level].size(); ++a){
|
||||
LatticeLorentzColourMatrix force(U._grid);
|
||||
as[level].at(a)->deriv(U,force);
|
||||
*P -= force*ep;
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
void update_U(LatticeLorentzColourMatrix&U, double ep){
|
||||
//rewrite exponential to deal automatically with the lorentz index?
|
||||
LatticeColourMatrix Umu(U._grid);
|
||||
LatticeColourMatrix Pmu(U._grid);
|
||||
for (int mu = 0; mu < Nd; mu++){
|
||||
Umu=peekLorentz(U, mu);
|
||||
Pmu=peekLorentz(*P, mu);
|
||||
Umu = expMat(Pmu, ep, Params.Nexp)*Umu;
|
||||
pokeLorentz(U, Umu, mu);
|
||||
}
|
||||
|
||||
}
|
||||
|
||||
|
||||
|
||||
friend void IntegratorAlgorithm::step (LatticeLorentzColourMatrix& U,
|
||||
int level, std::vector<int>& clock,
|
||||
Integrator<IntegratorAlgorithm>* Integ);
|
||||
public:
|
||||
Integrator(GridBase* grid, IntegratorParameters Par,
|
||||
ActionSet& Aset, std::vector<int> Nrel_):
|
||||
Params(Par),as(Aset),Nrel(Nrel_),P(new LatticeLorentzColourMatrix(grid)),pRNG(grid){
|
||||
assert(as.size() == Nrel.size());
|
||||
pRNG.SeedRandomDevice();
|
||||
};
|
||||
|
||||
~Integrator(){}
|
||||
|
||||
|
||||
//Initialization of momenta and actions
|
||||
void init(LatticeLorentzColourMatrix& U){
|
||||
std::cout<< "Integrator init\n";
|
||||
|
||||
MDutils::generate_momenta(*P,pRNG);
|
||||
for(int level=0; level< as.size(); ++level){
|
||||
for(int actionID=0; actionID<as.at(level).size(); ++actionID){
|
||||
as[level].at(actionID)->init(U, pRNG);
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
// Calculate action
|
||||
RealD S(LatticeLorentzColourMatrix& U){
|
||||
LatticeComplex Hloc(U._grid);
|
||||
Hloc = zero;
|
||||
// Momenta
|
||||
for (int mu=0; mu <Nd; mu++){
|
||||
LatticeColourMatrix Pmu = peekLorentz(*P, mu);
|
||||
Hloc -= trace(Pmu*Pmu);
|
||||
}
|
||||
Complex Hsum = sum(Hloc);
|
||||
|
||||
RealD H = Hsum.real();
|
||||
|
||||
std::cout << "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);
|
||||
|
||||
return H;
|
||||
}
|
||||
|
||||
void integrate(LatticeLorentzColourMatrix& U){
|
||||
std::vector<int> clock;
|
||||
clock.resize(as.size(),0);
|
||||
for(int step=0; step< Params.MDsteps; ++step) // MD step
|
||||
TheIntegrator.step(U,0,clock, (this));
|
||||
}
|
||||
};
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
}
|
||||
}
|
||||
#endif//INTEGRATOR_INCLUDED
|
||||
152
lib/qcd/hmc/integrators/Integrator_algorithm.h
Normal file
152
lib/qcd/hmc/integrators/Integrator_algorithm.h
Normal file
@ -0,0 +1,152 @@
|
||||
//--------------------------------------------------------------------
|
||||
/*! @file Integrator_algorithm.h
|
||||
* @brief Declaration of classes for the Molecular Dynamics algorithms
|
||||
*
|
||||
* @author Guido Cossu
|
||||
*/
|
||||
//--------------------------------------------------------------------
|
||||
|
||||
#ifndef INTEGRATOR_ALG_INCLUDED
|
||||
#define INTEGRATOR_ALG_INCLUDED
|
||||
|
||||
|
||||
namespace Grid{
|
||||
namespace QCD{
|
||||
|
||||
|
||||
class MinimumNorm2{
|
||||
const double lambda = 0.1931833275037836;
|
||||
public:
|
||||
void step (LatticeLorentzColourMatrix& U,
|
||||
int level, std::vector<int>& clock,
|
||||
Integrator<MinimumNorm2>* Integ){
|
||||
// level : current level
|
||||
// fl : final level
|
||||
// eps : current step size
|
||||
|
||||
int fl = Integ->as.size() -1;
|
||||
double eps = Integ->Params.stepsize;
|
||||
|
||||
for(int l=0; l<=level; ++l) eps/= 2.0*Integ->Nrel[l];
|
||||
|
||||
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
|
||||
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;
|
||||
}
|
||||
|
||||
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;
|
||||
}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;
|
||||
|
||||
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;
|
||||
}else{ // recursive function call
|
||||
step(U,level+1,clock, Integ);
|
||||
}
|
||||
|
||||
|
||||
if(clock[level] == fin){ // final 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;
|
||||
}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;
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
|
||||
}
|
||||
|
||||
};
|
||||
|
||||
class LeapFrog{
|
||||
public:
|
||||
void step (LatticeLorentzColourMatrix& U,
|
||||
int level, std::vector<int>& clock,
|
||||
Integrator<LeapFrog>* Integ){
|
||||
// level : current level
|
||||
// fl : final level
|
||||
// eps : current step size
|
||||
|
||||
int fl = Integ->as.size() -1;
|
||||
double eps = Integ->Params.stepsize;
|
||||
|
||||
// Get current level step size
|
||||
for(int l=0; l<=level; ++l) eps/= Integ->Nrel[l];
|
||||
|
||||
int fin = 1;
|
||||
for(int l=0; l<=level; ++l) fin*= Integ->Nrel[l];
|
||||
fin = 2*Integ->Params.MDsteps*fin - 1;
|
||||
|
||||
for(int e=0; e<Integ->Nrel[level]; ++e){
|
||||
|
||||
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;
|
||||
}
|
||||
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;
|
||||
}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;
|
||||
}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;
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
|
||||
|
||||
}
|
||||
};
|
||||
|
||||
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
|
||||
#endif//INTEGRATOR_INCLUDED
|
||||
@ -1,280 +0,0 @@
|
||||
//--------------------------------------------------------------------
|
||||
/*! @file Integrator_base.h
|
||||
* @brief Declaration of classes for the abstract Molecular Dynamics integrator
|
||||
*
|
||||
* @author Guido Cossu
|
||||
*/
|
||||
//--------------------------------------------------------------------
|
||||
|
||||
#ifndef INTEGRATOR_INCLUDED
|
||||
#define INTEGRATOR_INCLUDED
|
||||
|
||||
#include <memory>
|
||||
|
||||
class Observer;
|
||||
|
||||
|
||||
/*! @brief Abstract base class for Molecular Dynamics management */
|
||||
|
||||
namespace Grid{
|
||||
namespace QCD{
|
||||
|
||||
typedef Action<LatticeLorentzColourMatrix>* ActPtr; // now force the same size as the rest of the code
|
||||
typedef std::vector<ActPtr> ActionLevel;
|
||||
typedef std::vector<ActionLevel> ActionSet;
|
||||
typedef std::vector<Observer*> ObserverList;
|
||||
|
||||
struct IntegratorParameters{
|
||||
int Nexp;
|
||||
int MDsteps; // number of outer steps
|
||||
RealD trajL; // trajectory length
|
||||
RealD stepsize;
|
||||
|
||||
IntegratorParameters(int Nexp_,
|
||||
int MDsteps_,
|
||||
RealD trajL_):
|
||||
Nexp(Nexp_),MDsteps(MDsteps_),trajL(trajL_),stepsize(trajL/MDsteps){};
|
||||
};
|
||||
|
||||
|
||||
namespace MDutils{
|
||||
void generate_momenta(LatticeLorentzColourMatrix&,GridParallelRNG&);
|
||||
void generate_momenta_su3(LatticeLorentzColourMatrix&,GridParallelRNG&);
|
||||
}
|
||||
|
||||
|
||||
template< class IntegratorPolicy >
|
||||
class Integrator{
|
||||
private:
|
||||
IntegratorParameters Params;
|
||||
const ActionSet as;
|
||||
const std::vector<int> Nrel; //relative step size per level
|
||||
std::unique_ptr<LatticeLorentzColourMatrix> P;
|
||||
GridParallelRNG pRNG;
|
||||
//ObserverList observers; // not yet
|
||||
|
||||
IntegratorPolicy TheIntegrator;
|
||||
|
||||
void register_observers();
|
||||
void notify_observers();
|
||||
|
||||
void update_P(LatticeLorentzColourMatrix&U, int level,double ep){
|
||||
for(int a=0; a<as[level].size(); ++a){
|
||||
LatticeLorentzColourMatrix force(U._grid);
|
||||
as[level].at(a)->deriv(U,force);
|
||||
*P -= force*ep;
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
void update_U(LatticeLorentzColourMatrix&U, double ep){
|
||||
//rewrite exponential to deal automatically with the lorentz index?
|
||||
LatticeColourMatrix Umu(U._grid);
|
||||
LatticeColourMatrix Pmu(U._grid);
|
||||
for (int mu = 0; mu < Nd; mu++){
|
||||
Umu=peekLorentz(U, mu);
|
||||
Pmu=peekLorentz(*P, mu);
|
||||
Umu = expMat(Pmu, ep, Params.Nexp)*Umu;
|
||||
pokeLorentz(U, Umu, mu);
|
||||
}
|
||||
|
||||
}
|
||||
|
||||
|
||||
|
||||
friend void IntegratorPolicy::step (LatticeLorentzColourMatrix& U,
|
||||
int level, std::vector<int>& clock,
|
||||
Integrator<IntegratorPolicy>* Integ);
|
||||
public:
|
||||
Integrator(GridBase* grid, IntegratorParameters Par,
|
||||
ActionSet& Aset, std::vector<int> Nrel_):
|
||||
Params(Par),as(Aset),Nrel(Nrel_),P(new LatticeLorentzColourMatrix(grid)),pRNG(grid){
|
||||
assert(as.size() == Nrel.size());
|
||||
pRNG.SeedRandomDevice();
|
||||
};
|
||||
|
||||
~Integrator(){}
|
||||
|
||||
|
||||
//Initialization of momenta and actions
|
||||
void init(LatticeLorentzColourMatrix& U){
|
||||
std::cout<< "Integrator init\n";
|
||||
|
||||
MDutils::generate_momenta(*P,pRNG);
|
||||
for(int level=0; level< as.size(); ++level){
|
||||
for(int actionID=0; actionID<as.at(level).size(); ++actionID){
|
||||
as[level].at(actionID)->init(U, pRNG);
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
// Calculate action
|
||||
RealD S(LatticeLorentzColourMatrix& U){
|
||||
LatticeComplex Hloc(U._grid);
|
||||
Hloc = zero;
|
||||
// Momenta
|
||||
for (int mu=0; mu <Nd; mu++){
|
||||
LatticeColourMatrix Pmu = peekLorentz(*P, mu);
|
||||
Hloc -= trace(Pmu*Pmu);
|
||||
}
|
||||
Complex Hsum = sum(Hloc);
|
||||
|
||||
RealD H = Hsum.real();
|
||||
|
||||
std::cout << "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);
|
||||
|
||||
return H;
|
||||
}
|
||||
|
||||
|
||||
|
||||
void integrate(LatticeLorentzColourMatrix& U, int level){
|
||||
std::vector<int> clock;
|
||||
clock.resize(as.size(),0);
|
||||
for(int step=0; step< Params.MDsteps; ++step) // MD step
|
||||
TheIntegrator.step(U,0,clock, (this));
|
||||
}
|
||||
};
|
||||
|
||||
|
||||
class MinimumNorm2{
|
||||
const double lambda = 0.1931833275037836;
|
||||
public:
|
||||
void step (LatticeLorentzColourMatrix& U,
|
||||
int level, std::vector<int>& clock,
|
||||
Integrator<MinimumNorm2>* Integ){
|
||||
// level : current level
|
||||
// fl : final level
|
||||
// eps : current step size
|
||||
|
||||
int fl = Integ->as.size() -1;
|
||||
double eps = Integ->Params.stepsize;
|
||||
|
||||
for(int l=0; l<=level; ++l) eps/= 2.0*Integ->Nrel[l];
|
||||
|
||||
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
|
||||
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;
|
||||
}
|
||||
|
||||
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;
|
||||
}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;
|
||||
|
||||
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;
|
||||
}else{ // recursive function call
|
||||
step(U,level+1,clock, Integ);
|
||||
}
|
||||
|
||||
|
||||
if(clock[level] == fin){ // final 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;
|
||||
}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;
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
|
||||
}
|
||||
|
||||
};
|
||||
|
||||
class LeapFrog{
|
||||
public:
|
||||
void step (LatticeLorentzColourMatrix& U,
|
||||
int level, std::vector<int>& clock,
|
||||
Integrator<LeapFrog>* Integ){
|
||||
// fl : final level
|
||||
// eps : current step size
|
||||
|
||||
int fl = Integ->as.size() -1;
|
||||
double eps = Integ->Params.stepsize;
|
||||
|
||||
// Get current level step size
|
||||
for(int l=0; l<=level; ++l) eps/= Integ->Nrel[l];
|
||||
|
||||
int fin = 1;
|
||||
for(int l=0; l<=level; ++l) fin*= Integ->Nrel[l];
|
||||
fin = 2*Integ->Params.MDsteps*fin - 1;
|
||||
|
||||
for(int e=0; e<Integ->Nrel[level]; ++e){
|
||||
|
||||
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;
|
||||
}
|
||||
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;
|
||||
}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;
|
||||
}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;
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
|
||||
|
||||
}
|
||||
};
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
}
|
||||
}
|
||||
#endif//INTEGRATOR_INCLUDED
|
||||
@ -40,7 +40,7 @@ int main (int argc, char ** argv)
|
||||
|
||||
// Create integrator
|
||||
typedef MinimumNorm2 IntegratorAlgorithm;// change here to change the algorithm
|
||||
IntegratorParameters MDpar(12,30,1.0);
|
||||
IntegratorParameters MDpar(12,5,1.0);
|
||||
std::vector<int> rel ={1};
|
||||
Integrator<IntegratorAlgorithm> MDynamics(&Fine,MDpar, FullSet,rel);
|
||||
|
||||
|
||||
Loading…
Reference in New Issue
Block a user