diff --git a/lib/qcd/action/Actions.h b/lib/qcd/action/Actions.h index e45f2e0c..ef2edc97 100644 --- a/lib/qcd/action/Actions.h +++ b/lib/qcd/action/Actions.h @@ -72,6 +72,10 @@ namespace QCD { typedef WilsonGaugeAction WilsonGaugeActionR; typedef WilsonGaugeAction WilsonGaugeActionF; typedef WilsonGaugeAction WilsonGaugeActionD; +typedef VariableWilsonGaugeAction VariableWilsonGaugeActionR; +typedef VariableWilsonGaugeAction VariableWilsonGaugeActionF; +typedef VariableWilsonGaugeAction VariableWilsonGaugeActionD; + typedef PlaqPlusRectangleAction PlaqPlusRectangleActionR; typedef PlaqPlusRectangleAction PlaqPlusRectangleActionF; typedef PlaqPlusRectangleAction PlaqPlusRectangleActionD; diff --git a/lib/qcd/action/gauge/WilsonGaugeAction.h b/lib/qcd/action/gauge/WilsonGaugeAction.h index aff67c67..686c5470 100644 --- a/lib/qcd/action/gauge/WilsonGaugeAction.h +++ b/lib/qcd/action/gauge/WilsonGaugeAction.h @@ -1,86 +1,260 @@ - /************************************************************************************* +/************************************************************************************* - Grid physics library, www.github.com/paboyle/Grid +Grid physics library, www.github.com/paboyle/Grid - Source file: ./lib/qcd/action/gauge/WilsonGaugeAction.h +Source file: ./lib/qcd/action/gauge/WilsonGaugeAction.h - Copyright (C) 2015 +Copyright (C) 2015 Author: Azusa Yamaguchi Author: Peter Boyle Author: neo Author: paboyle - This program is free software; you can redistribute it and/or modify - it under the terms of the GNU General Public License as published by - the Free Software Foundation; either version 2 of the License, or - (at your option) any later version. +This program is free software; you can redistribute it and/or modify +it under the terms of the GNU General Public License as published by +the Free Software Foundation; either version 2 of the License, or +(at your option) any later version. - This program is distributed in the hope that it will be useful, - but WITHOUT ANY WARRANTY; without even the implied warranty of - MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the - GNU General Public License for more details. +This program is distributed in the hope that it will be useful, +but WITHOUT ANY WARRANTY; without even the implied warranty of +MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +GNU General Public License for more details. - You should have received a copy of the GNU General Public License along - with this program; if not, write to the Free Software Foundation, Inc., - 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. +You should have received a copy of the GNU General Public License along +with this program; if not, write to the Free Software Foundation, Inc., +51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. - See the full license in the file "LICENSE" in the top level distribution directory - *************************************************************************************/ - /* END LEGAL */ +See the full license in the file "LICENSE" in the top level distribution +directory +*************************************************************************************/ +/* END LEGAL */ #ifndef QCD_WILSON_GAUGE_ACTION_H #define QCD_WILSON_GAUGE_ACTION_H -namespace Grid{ - namespace QCD{ - - //////////////////////////////////////////////////////////////////////// - // Wilson Gauge Action .. should I template the Nc etc.. - //////////////////////////////////////////////////////////////////////// - template - class WilsonGaugeAction : public Action { - public: +namespace Grid { +namespace QCD { - INHERIT_GIMPL_TYPES(Gimpl); +//////////////////////////////////////////////////////////////////////// +// Wilson Gauge Action .. should I template the Nc etc.. +//////////////////////////////////////////////////////////////////////// +template +class WilsonGaugeAction : public Action { + public: + INHERIT_GIMPL_TYPES(Gimpl); - // typedef LorentzScalar GaugeLinkField; + // typedef LorentzScalar GaugeLinkField; - private: - RealD beta; - public: - WilsonGaugeAction(RealD b):beta(b){}; + private: + RealD beta; + + public: + WilsonGaugeAction(RealD b) : beta(b){}; + + virtual void refresh(const GaugeField &U, + GridParallelRNG &pRNG){}; // noop as no pseudoferms + + virtual RealD S(const GaugeField &U) { + RealD plaq = WilsonLoops::avgPlaquette(U); + RealD vol = U._grid->gSites(); + RealD action = beta * (1.0 - plaq) * (Nd * (Nd - 1.0)) * vol * 0.5; + return action; + }; + + virtual void deriv(const GaugeField &U, GaugeField &dSdU) { + // not optimal implementation FIXME + // extend Ta to include Lorentz indexes + + // RealD factor = 0.5*beta/RealD(Nc); + RealD factor = 0.5 * beta / RealD(Nc); + + GaugeLinkField Umu(U._grid); + GaugeLinkField dSdU_mu(U._grid); + for (int mu = 0; mu < Nd; mu++) { + Umu = PeekIndex(U, mu); + + // Staple in direction mu + WilsonLoops::Staple(dSdU_mu, U, mu); + dSdU_mu = Ta(Umu * dSdU_mu) * factor; - virtual void refresh(const GaugeField &U, GridParallelRNG& pRNG) {}; // noop as no pseudoferms - - virtual RealD S(const GaugeField &U) { - RealD plaq = WilsonLoops::avgPlaquette(U); - RealD vol = U._grid->gSites(); - RealD action=beta*(1.0 -plaq)*(Nd*(Nd-1.0))*vol*0.5; - return action; - }; + PokeIndex(dSdU, dSdU_mu, mu); + } + }; +}; - virtual void deriv(const GaugeField &U,GaugeField & dSdU) { - //not optimal implementation FIXME - //extend Ta to include Lorentz indexes +template +class VariableWilsonGaugeAction : public Action { + public: + INHERIT_GIMPL_TYPES(Gimpl); - //RealD factor = 0.5*beta/RealD(Nc); - RealD factor = 0.5*beta/RealD(Nc); + private: + std::vector b_bulk;// bulk couplings + std::vector b_xdim;//extra dimension couplings + GridBase *grid; + LatticeComplex beta_xdim; + LatticeComplex beta_xdim_shifted; + LatticeComplex beta_bulk; - GaugeLinkField Umu(U._grid); - GaugeLinkField dSdU_mu(U._grid); - for (int mu=0; mu < Nd; mu++){ + public: + VariableWilsonGaugeAction(std::vector bulk, std::vector xdim, + GridBase *_grid) + : b_bulk(bulk), + b_xdim(xdim), + grid(_grid), + beta_xdim(grid), + beta_xdim_shifted(grid), + beta_bulk(grid) + { + //check that the grid is ok + //todo + int Ndim = Nd;//change later - Umu = PeekIndex(U,mu); + LatticeComplex temp(grid); - // Staple in direction mu - WilsonLoops::Staple(dSdU_mu,U,mu); - dSdU_mu = Ta(Umu*dSdU_mu)*factor; - PokeIndex(dSdU, dSdU_mu, mu); - } - }; - }; - + Lattice > coor(grid); + + LatticeCoordinate(coor, Ndim - 1); + + int Nex = grid->_fdimensions[Ndim - 1]; + assert(b_bulk.size() == Nex); + assert(b_xdim.size() == Nex); + + beta_xdim = zero; + for (int tau = 0; tau < Nex; tau++) { + temp = b_xdim[tau]; + beta_xdim = where(coor == tau, temp, beta_xdim); + } + + beta_xdim_shifted = Cshift(beta_xdim, Ndim - 1, -1); + + beta_bulk = zero; + for (int tau = 0; tau < Nex; tau++) { + temp = b_bulk[tau]; + beta_bulk = where(coor == tau, temp, beta_bulk); + } + }; + + virtual void refresh(const GaugeField &U, + GridParallelRNG &pRNG){}; // noop as no pseudoferms + + virtual RealD S(const GaugeField &Umu) { + int Ndim = Nd; // change later for generality + conformable(grid, Umu._grid); + + std::vector U(Ndim, grid); + + for (int mu = 0; mu < Ndim; mu++) { + U[mu] = PeekIndex(Umu, mu); + } + + LatticeComplex dirPlaq(grid); + LatticeComplex Plaq(grid); + + RealD OneOnNc = 1.0 / Real(Nc); + + ///////////// + // Lower dim plaquettes + ///////////// + Plaq = zero; + for (int mu = 1; mu < Ndim - 1; mu++) { + for (int nu = 0; nu < mu; nu++) { + WilsonLoops::traceDirPlaquette(dirPlaq, U, mu, nu); + Plaq = Plaq + (1.0 - dirPlaq * OneOnNc) * beta_bulk; + } + } + + ///////////// + // Extra dimension + ///////////// + { + int mu = Ndim - 1; + for (int nu = 0; nu < mu; nu++) { + WilsonLoops::traceDirPlaquette(dirPlaq, U, mu, nu); + Plaq = Plaq + (1.0 - dirPlaq * OneOnNc) * beta_xdim; + } + } + + TComplex Tp = sum(Plaq); + Complex p = TensorRemove(Tp); + RealD action = p.real(); + return action; + }; + + virtual void deriv(const GaugeField &U, GaugeField &dSdU) { + // not optimal implementation FIXME + // extend Ta to include Lorentz indexes + + // for the higher dimension plaquettes take the upper plaq of the + // 4d slice and multiply by beta[s] and the lower and multiply by beta[s-1] + // todo + // derivative of links mu = 0, ... Nd-1 inside plaq (mu,5) + // for these I need upper and lower staples separated + // each multiplied with their own beta + // derivative of links mu = 5 + // living on the same slice, share the same beta + + +conformable(grid,U._grid); +int Ndim = Nd; // change later +RealD factor = 0.5 / RealD(Nc); + +GaugeLinkField Umu(grid); +GaugeLinkField dSdU_mu(grid); +GaugeLinkField staple(grid); + +for (int mu = 0; mu < Ndim; mu++) { + Umu = PeekIndex(U, mu); + dSdU_mu = zero; + + for (int nu = 0; nu < Ndim; nu++) { + if (nu != mu) { + if ((mu < (Ndim - 1)) && (nu < (Ndim - 1))) { + // Spacelike case apply beta space + WilsonLoops::Staple(staple, U, mu, nu); + staple = staple * beta_bulk; + dSdU_mu += staple; + + } else if (mu == (Ndim - 1)) { + // nu space; mu time link + assert(nu < (Ndim - 1)); + assert(mu == (Ndim - 1)); + + // mu==tau dir link deriv, nu spatial + WilsonLoops::Staple(staple, U, mu, nu); + staple = staple * beta_xdim; + dSdU_mu += staple; + + } else { + assert(mu < (Ndim - 1)); + assert(nu == (Ndim - 1)); + + // nu time; mu space link + + // staple forwards in tau + WilsonLoops::StapleUpper(staple, U, mu, nu); + staple = staple * beta_xdim; + dSdU_mu += staple; + + // staple backwards in tau + WilsonLoops::StapleLower(staple, U, mu, nu); + staple = staple * beta_xdim_shifted; + dSdU_mu += staple; + } + } } + + dSdU_mu = Ta(Umu * dSdU_mu) * factor; + PokeIndex(dSdU, dSdU_mu, mu); + } + + + + }; +}; + + + +} } #endif diff --git a/lib/qcd/hmc/UsingHMC.md b/lib/qcd/hmc/UsingHMC.md new file mode 100644 index 00000000..da3c3c00 --- /dev/null +++ b/lib/qcd/hmc/UsingHMC.md @@ -0,0 +1,107 @@ +Using HMC in Grid version 0.5.1 + +These are the instructions to use the Generalised HMC on Grid version 0.5.1. +Disclaimer: GRID is still under active development so any information here can be changed in future releases. + + +Command line options +=================== +(relevant file GenericHMCrunner.h) +The initial configuration can be changed at the command line using +--StartType +valid choices, one among these +HotStart, ColdStart, TepidStart, CheckpointStart +default: HotStart + +example +./My_hmc_exec --StartType HotStart + +The CheckpointStart option uses the prefix for the configurations and rng seed files defined in your executable and the initial configuration is specified by +--StartTrajectory +default: 0 + +The number of trajectories for a specific run are specified at command line by +--Trajectories +default: 1 + +The number of thermalization steps (i.e. steps when the Metropolis acceptance check is turned off) is specified by +--Thermalizations +default: 10 + + +Any other parameter is defined in the source for the executable. + +HMC controls +=========== + +The lines + + std::vector SerSeed({1, 2, 3, 4, 5}); + std::vector ParSeed({6, 7, 8, 9, 10}); + +define the seeds for the serial and the parallel RNG. + +The line + + TheHMC.MDparameters.set(20, 1.0);// MDsteps, traj length + +declares the number of molecular dynamics steps and the total trajectory length. + + +Actions +====== + +Action names are defined in the file +lib/qcd/Actions.h + +Gauge actions list: + +WilsonGaugeActionR; +WilsonGaugeActionF; +WilsonGaugeActionD; +PlaqPlusRectangleActionR; +PlaqPlusRectangleActionF; +PlaqPlusRectangleActionD; +IwasakiGaugeActionR; +IwasakiGaugeActionF; +IwasakiGaugeActionD; +SymanzikGaugeActionR; +SymanzikGaugeActionF; +SymanzikGaugeActionD; + + +ConjugateWilsonGaugeActionR; +ConjugateWilsonGaugeActionF; +ConjugateWilsonGaugeActionD; +ConjugatePlaqPlusRectangleActionR; +ConjugatePlaqPlusRectangleActionF; +ConjugatePlaqPlusRectangleActionD; +ConjugateIwasakiGaugeActionR; +ConjugateIwasakiGaugeActionF; +ConjugateIwasakiGaugeActionD; +ConjugateSymanzikGaugeActionR; +ConjugateSymanzikGaugeActionF; +ConjugateSymanzikGaugeActionD; + + +ScalarActionR; +ScalarActionF; +ScalarActionD; + + +each of these action accept one single parameter at creation time (beta). +Example for creating a Symanzik action with beta=4.0 + + SymanzikGaugeActionR(4.0) + +The suffixes R,F,D in the action names refer to the Real +(the precision is defined at compile time by the --enable-precision flag in the configure), +Float and Double, that force the precision of the action to be 32, 64 bit respectively. + + + + + + + + diff --git a/lib/qcd/utils/WilsonLoops.h b/lib/qcd/utils/WilsonLoops.h index 10022f50..3051f830 100644 --- a/lib/qcd/utils/WilsonLoops.h +++ b/lib/qcd/utils/WilsonLoops.h @@ -86,8 +86,7 @@ public: // sum over all x,y,z,t and over all planes of plaquette ////////////////////////////////////////////////// static RealD sumPlaquette(const GaugeLorentz &Umu) { - std::vector U(4, Umu._grid); - + std::vector U(Nd, Umu._grid); for (int mu = 0; mu < Nd; mu++) { U[mu] = PeekIndex(Umu, mu); } @@ -95,11 +94,12 @@ public: LatticeComplex Plaq(Umu._grid); sitePlaquette(Plaq, U); - TComplex Tp = sum(Plaq); Complex p = TensorRemove(Tp); return p.real(); } + + ////////////////////////////////////////////////// // average over all x,y,z,t and over all planes of plaquette ////////////////////////////////////////////////// @@ -114,7 +114,7 @@ public: // average over traced single links ////////////////////////////////////////////////// static RealD linkTrace(const GaugeLorentz &Umu) { - std::vector U(4, Umu._grid); + std::vector U(Nd, Umu._grid); LatticeComplex Tr(Umu._grid); Tr = zero; @@ -139,7 +139,7 @@ public: GridBase *grid = Umu._grid; - std::vector U(4, grid); + std::vector U(Nd, grid); for (int d = 0; d < Nd; d++) { U[d] = PeekIndex(Umu, d); } @@ -233,7 +233,7 @@ public: if (nu != mu) { GridBase *grid = Umu._grid; - std::vector U(4, grid); + std::vector U(Nd, grid); for (int d = 0; d < Nd; d++) { U[d] = PeekIndex(Umu, d); } @@ -256,6 +256,37 @@ public: } } + ////////////////////////////////////////////////// + // the sum over all staples on each site in direction mu,nu, lower part + ////////////////////////////////////////////////// + static void StapleLower(GaugeMat &staple, const GaugeLorentz &Umu, int mu, + int nu) { + staple = zero; + + if (nu != mu) { + GridBase *grid = Umu._grid; + + std::vector U(Nd, grid); + for (int d = 0; d < Nd; d++) { + U[d] = PeekIndex(Umu, d); + } + + // mu + // ^ + // |__> nu + + // __ + // | + // |__ + // + // + staple += Gimpl::ShiftStaple( + Gimpl::CovShiftBackward(U[nu], nu, + Gimpl::CovShiftBackward(U[mu], mu, U[nu])), + mu); + } + } + ////////////////////////////////////////////////////// // Similar to above for rectangle is required ////////////////////////////////////////////////////// @@ -375,8 +406,8 @@ public: // |___ ___| // - // tmp= Staple2x1* Cshift(U[mu],mu,-2); - // Stap+= Cshift(tmp,mu,1) ; + // tmp= Staple2x1* Cshift(U[mu],mu,-2); + // Stap+= Cshift(tmp,mu,1) ; Stap += Cshift(Staple2x1, mu, 1) * Cshift(U[mu], mu, -1); ; diff --git a/tests/hmc/Make.inc b/tests/hmc/Make.inc index b37d797e..f4a3603f 100644 --- a/tests/hmc/Make.inc +++ b/tests/hmc/Make.inc @@ -1,5 +1,5 @@ -tests: Test_hmc_EODWFRatio_Binary Test_hmc_EODWFRatio Test_hmc_EODWFRatio_Gparity Test_hmc_EOWilsonFermionGauge Test_hmc_EOWilsonRatio Test_hmc_GparityIwasakiGauge Test_hmc_GparityWilsonGauge Test_hmc_IwasakiGauge Test_hmc_RectGauge Test_hmc_ScalarAction Test_hmc_WilsonAdjointFermionGauge Test_hmc_WilsonFermionGauge_Binary Test_hmc_WilsonFermionGauge Test_hmc_WilsonGauge Test_hmc_WilsonMixedRepresentationsFermionGauge Test_hmc_WilsonRatio Test_hmc_WilsonTwoIndexSymmetricFermionGauge Test_multishift_sqrt Test_remez Test_rhmc_EOWilson1p1 Test_rhmc_EOWilsonRatio Test_rhmc_Wilson1p1 Test_rhmc_WilsonRatio -EXTRA_PROGRAMS = Test_hmc_EODWFRatio_Binary Test_hmc_EODWFRatio Test_hmc_EODWFRatio_Gparity Test_hmc_EOWilsonFermionGauge Test_hmc_EOWilsonRatio Test_hmc_GparityIwasakiGauge Test_hmc_GparityWilsonGauge Test_hmc_IwasakiGauge Test_hmc_RectGauge Test_hmc_ScalarAction Test_hmc_WilsonAdjointFermionGauge Test_hmc_WilsonFermionGauge_Binary Test_hmc_WilsonFermionGauge Test_hmc_WilsonGauge Test_hmc_WilsonMixedRepresentationsFermionGauge Test_hmc_WilsonRatio Test_hmc_WilsonTwoIndexSymmetricFermionGauge Test_multishift_sqrt Test_remez Test_rhmc_EOWilson1p1 Test_rhmc_EOWilsonRatio Test_rhmc_Wilson1p1 Test_rhmc_WilsonRatio +tests: Test_hmc_EODWFRatio_Binary Test_hmc_EODWFRatio Test_hmc_EODWFRatio_Gparity Test_hmc_EOWilsonFermionGauge Test_hmc_EOWilsonRatio Test_hmc_GparityIwasakiGauge Test_hmc_GparityWilsonGauge Test_hmc_IwasakiGauge Test_hmc_RectGauge Test_hmc_ScalarAction Test_hmc_WilsonAdjointFermionGauge Test_hmc_WilsonFermionGauge_Binary Test_hmc_WilsonFermionGauge Test_hmc_WilsonGauge_Binary Test_hmc_WilsonGauge Test_hmc_WilsonMixedRepresentationsFermionGauge Test_hmc_WilsonRatio Test_hmc_WilsonTwoIndexSymmetricFermionGauge Test_multishift_sqrt Test_remez Test_rhmc_EOWilson1p1 Test_rhmc_EOWilsonRatio Test_rhmc_Wilson1p1 Test_rhmc_WilsonRatio +EXTRA_PROGRAMS = Test_hmc_EODWFRatio_Binary Test_hmc_EODWFRatio Test_hmc_EODWFRatio_Gparity Test_hmc_EOWilsonFermionGauge Test_hmc_EOWilsonRatio Test_hmc_GparityIwasakiGauge Test_hmc_GparityWilsonGauge Test_hmc_IwasakiGauge Test_hmc_RectGauge Test_hmc_ScalarAction Test_hmc_WilsonAdjointFermionGauge Test_hmc_WilsonFermionGauge_Binary Test_hmc_WilsonFermionGauge Test_hmc_WilsonGauge_Binary Test_hmc_WilsonGauge Test_hmc_WilsonMixedRepresentationsFermionGauge Test_hmc_WilsonRatio Test_hmc_WilsonTwoIndexSymmetricFermionGauge Test_multishift_sqrt Test_remez Test_rhmc_EOWilson1p1 Test_rhmc_EOWilsonRatio Test_rhmc_Wilson1p1 Test_rhmc_WilsonRatio Test_hmc_EODWFRatio_Binary_SOURCES=Test_hmc_EODWFRatio_Binary.cc Test_hmc_EODWFRatio_Binary_LDADD=-lGrid @@ -40,6 +40,9 @@ Test_hmc_WilsonFermionGauge_Binary_LDADD=-lGrid Test_hmc_WilsonFermionGauge_SOURCES=Test_hmc_WilsonFermionGauge.cc Test_hmc_WilsonFermionGauge_LDADD=-lGrid +Test_hmc_WilsonGauge_Binary_SOURCES=Test_hmc_WilsonGauge_Binary.cc +Test_hmc_WilsonGauge_Binary_LDADD=-lGrid + Test_hmc_WilsonGauge_SOURCES=Test_hmc_WilsonGauge.cc Test_hmc_WilsonGauge_LDADD=-lGrid diff --git a/tests/hmc/Test_hmc_EODWFRatio_Binary.cc b/tests/hmc/Test_hmc_EODWFRatio_Binary.cc index cd35968a..a7fb4db1 100644 --- a/tests/hmc/Test_hmc_EODWFRatio_Binary.cc +++ b/tests/hmc/Test_hmc_EODWFRatio_Binary.cc @@ -78,7 +78,8 @@ class HmcRunner : public BinaryHmcRunner { LatticeGaugeField U(UGrid); // Gauge action - IwasakiGaugeActionR Iaction(4.0); + double beta = 4.0; + IwasakiGaugeActionR Iaction(beta); Real mass = 0.04; Real pv = 1.0; @@ -87,7 +88,9 @@ class HmcRunner : public BinaryHmcRunner { FermionAction DenOp(U,*FGrid,*FrbGrid,*UGrid,*UrbGrid,mass,M5,scale); FermionAction NumOp(U,*FGrid,*FrbGrid,*UGrid,*UrbGrid,pv,M5,scale); - ConjugateGradient CG(1.0e-8,10000); + double StoppingCondition = 1.0e-8; + double MaxCGIterations = 10000; + ConjugateGradient CG(StoppingCondition,MaxCGIterations); TwoFlavourEvenOddRatioPseudoFermionAction Nf2(NumOp, DenOp,CG,CG); // Set smearing (true/false), default: false @@ -98,6 +101,9 @@ class HmcRunner : public BinaryHmcRunner { ActionLevel Level1(1); Level1.push_back(&Nf2); + // this level will integrate with a + // step that is 4 times finer + // than the previous level ActionLevel Level2(4); Level2.push_back(&Iaction); diff --git a/tests/hmc/Test_hmc_WilsonGauge_Binary.cc b/tests/hmc/Test_hmc_WilsonGauge_Binary.cc new file mode 100644 index 00000000..d6041679 --- /dev/null +++ b/tests/hmc/Test_hmc_WilsonGauge_Binary.cc @@ -0,0 +1,136 @@ +/************************************************************************************* + +Grid physics library, www.github.com/paboyle/Grid + +Source file: ./tests/Test_hmc_WilsonFermionGauge.cc + +Copyright (C) 2015 + +Author: Peter Boyle +Author: neo +Author: Guido Cossu + +This program is free software; you can redistribute it and/or modify +it under the terms of the GNU General Public License as published by +the Free Software Foundation; either version 2 of the License, or +(at your option) any later version. + +This program is distributed in the hope that it will be useful, +but WITHOUT ANY WARRANTY; without even the implied warranty of +MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +GNU General Public License for more details. + +You should have received a copy of the GNU General Public License along +with this program; if not, write to the Free Software Foundation, Inc., +51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. + +See the full license in the file "LICENSE" in the top level distribution +directory +*************************************************************************************/ +/* END LEGAL */ +#include + +using namespace std; +using namespace Grid; +using namespace Grid::QCD; + +namespace Grid { +namespace QCD { + +class HMCRunnerParameters : Serializable { + public: + GRID_SERIALIZABLE_CLASS_MEMBERS(HMCRunnerParameters, + double, beta, + double, mass, + int, MaxCGIterations, + double, StoppingCondition, + bool, smearedAction, + int, SaveInterval, + std::string, format, + std::string, conf_prefix, + std::string, rng_prefix, + double, rho, + int, SmearingLevels, + ); + + HMCRunnerParameters() {} +}; + +// Derive from the BinaryHmcRunner (templated for gauge fields) +class HmcRunner : public BinaryHmcRunner { + public: + void BuildTheAction(int argc, char **argv) + + { + typedef WilsonImplR ImplPolicy; + typedef WilsonFermionR FermionAction; + typedef typename FermionAction::FermionField FermionField; + + UGrid = SpaceTimeGrid::makeFourDimGrid( + GridDefaultLatt(), GridDefaultSimd(Nd, vComplex::Nsimd()), + GridDefaultMpi()); + UrbGrid = SpaceTimeGrid::makeFourDimRedBlackGrid(UGrid); + + FGrid = UGrid; + FrbGrid = UrbGrid; + + // temporarily need a gauge field + LatticeGaugeField U(UGrid); + + // Gauge action + int Ls = UGrid->_gdimensions[Nd - 1]; + std::vector betat(Ls,5); + std::vector betas(Ls); + betas={5,6,6,5}; + std:cout << "Betas:" << betas << std::endl; + VariableWilsonGaugeActionR Waction(betas, betat, UGrid); + //WilsonGaugeActionR Waction(5.6); + + // Collect actions + ActionLevel Level1(1); + Level1.push_back(&Waction); + TheAction.push_back(Level1); + + // Add observables + int SaveInterval = 1; + std::string format = std::string("IEEE64BIG"); + std::string conf_prefix = std::string("ckpoint_lat"); + std::string rng_prefix = std::string("ckpoint_rng"); + BinaryHmcCheckpointer Checkpoint( + conf_prefix, rng_prefix, SaveInterval, format); + // Can implement also a specific function in the hmcrunner + // AddCheckpoint (...) that takes the same parameters + a string/tag + // defining the type of the checkpointer + // with tags can be implemented by overloading and no ifs + // Then force all checkpoint to have few common functions + // return an object that is then passed to the Run function + + PlaquetteLogger PlaqLog( + std::string("Plaquette")); + ObservablesList.push_back(&PlaqLog); + ObservablesList.push_back(&Checkpoint); + + Run(argc, argv, Checkpoint); // no smearing + }; +}; +} +} + +int main(int argc, char **argv) { + Grid_init(&argc, &argv); + + int threads = GridThread::GetThreads(); + std::cout << GridLogMessage << "Grid is setup to use " << threads + << " threads" << std::endl; + + HmcRunner TheHMC; + + // Seeds for the random number generators + std::vector SerSeed({1, 2, 3, 4, 5}); + std::vector ParSeed({6, 7, 8, 9, 10}); + TheHMC.RNGSeeds(SerSeed, ParSeed); + + TheHMC.MDparameters.set(20, 1.0);// MDsteps, traj length + + TheHMC.BuildTheAction(argc, argv); +}