Source code of RHMC (from SIMULATeQCD)
/*
* main_rhmc.cpp
*
* This is the main to use RHMC to generate Nf=2+1 HISQ configurations. By default it will also measure
* the chiral condensate.
*
*/
#include "../simulateqcd.h"
#include "../modules/rhmc/rhmc.h"
#include "../modules/observables/polyakovLoop.h"
#include "../modules/dslash/condensate.h"
int main(int argc, char *argv[]) {
stdLogger.setVerbosity(INFO);
CommunicationBase commBase(&argc, &argv);
RhmcParameters param;
param.readfile(commBase, "../parameter/tests/rhmcTest.param", argc, argv);
if (not param.confnumber.isSet()){
param.confnumber.set(0);
}
const int HaloDepth = 2;
const size_t Nmeas = 10; // Number of RHS in Condensate measurements
std::string rand_file;
std::string gauge_file;
RationalCoeff rat;
rat.readfile(commBase, param.rat_file());
rat.check_rat(param);
commBase.init(param.nodeDim(), param.gpuTopo());
StopWatch<true> timer, totaltime;
typedef float floatT; // Define the precision here
rootLogger.info("STARTING RHMC Update:");
if (sizeof(floatT)==4) {
rootLogger.info("update done in single precision");
} else if(sizeof(floatT)==8) {
rootLogger.info("update done in double precision");
} else {
rootLogger.info("update done in unknown precision");
}
initIndexer(4,param, commBase);
Gaugefield<floatT, true, HaloDepth> gauge(commBase);
grnd_state<true> d_rand;
if (param.load_rand()) {
grnd_state<false> h_rand;
rand_file = param.rand_file() + std::to_string(param.confnumber());
rootLogger.info("With random numbers from file: " , rand_file);
h_rand.read_from_file(rand_file, commBase);
d_rand=h_rand;
} else {
rootLogger.info("With new random numbers generated from seed: " , param.seed());
initialize_rng(param.seed(), d_rand);
}
if (param.load_conf() == 0) {
rootLogger.info("Starting from unit configuration");
gauge.one();
} else if(param.load_conf() == 1) {
rootLogger.info("Starting from random configuration");
gauge.random(d_rand.state);
} else if(param.load_conf() == 2) {
gauge_file = param.gauge_file() + std::to_string(param.confnumber());
rootLogger.info("Starting from configuration: " , gauge_file);
gauge.readconf_nersc(gauge_file);
}
gauge.updateAll();
gauge.su3latunitarize();
if(param.always_acc()) {
rootLogger.warn("Skipping Metropolis step!");
}
rhmc<floatT, true, HaloDepth> HMC(param, rat, gauge, d_rand.state);
rootLogger.info("constructed the HMC");
HMC.init_ratapprox();
int acc = 0;
floatT acceptance = 0.0;
PolyakovLoop<floatT, true, HaloDepth, R18> ploop(gauge);
GaugeAction<floatT, true, HaloDepth, R18> gaugeaction(gauge);
for (int i = 1; i <= param.no_updates(); ++i) {
timer.start();
acc += HMC.update(!param.always_acc());
acceptance = floatT(acc)/floatT(i);
rootLogger.info("current acceptance = ", acceptance);
rootLogger.info("MEASUREMENT: ", param.confnumber()+i);
rootLogger.info("Polyakov Loop = ", ploop.getPolyakovLoop());
rootLogger.info("Plaquette = ", gaugeaction.plaquette());
rootLogger.info("Rectangle = ", gaugeaction.rectangle());
SimpleArray<double,Nmeas> chi_l = measure_condensate<floatT, true, 2, 4, Nmeas>(commBase, param, param.m_ud(), gauge, d_rand);
for (int j = 0; j < Nmeas; ++j) {
rootLogger.info("CHI_UD = ", chi_l[j]);
}
SimpleArray<double,Nmeas> chi_s = measure_condensate<floatT, true, 2, 4, Nmeas>(commBase, param, param.m_s(), gauge, d_rand);
for (int j = 0; j < Nmeas; ++j) {
rootLogger.info("CHI_S = ", chi_s[j]);
}
timer.stop();
totaltime += timer;
rootLogger.info("Time (TTRAJ) for trajectory without IO: ", timer, " | avg traj. time : " , totaltime/i);
timer.reset();
if (i % param.write_every()==0) {
gauge_file = param.gauge_file() + std::to_string(param.confnumber()+i);
rand_file = param.rand_file() + std::to_string(param.confnumber()+i);
gauge.writeconf_nersc(gauge_file);
grnd_state<false> h_rand;
h_rand=d_rand;
h_rand.write_to_file(rand_file, commBase);
}
}
rootLogger.info("Run has ended! acceptance = " , acceptance);
return 0;
}
[The rational hybrid Monte Carlo (rhmc)](https://doi.org/10.1016/S0920-5632(99)85217-7)
algorithm is a way of updating gauge fields when simulating dynamical fermions.
At the end of the rhmc process, usually called a trajectory, there is a
typical Metropolis accept/reject step. During the trajectory, a proposal gauge
field is generated by integrating some ficitious, Hamiltonian equations of motion
in Monte Carlo time.
The trajectory is pushed by a fictious driving force. There comes contributions from
the gauge part of the action as well as the fermion part. Schematically, it comes
out to something of the form
where \(\phi\) is a pseudo-fermion field, and where the \(\partial\) indicates a
Lie group derivative. Hence we need to spend some effort finding
the vector \(X=\left(D D^{\dagger}\right)^{-1}\phi\). Since \(D\) depends on the gauge
field, we have to solve for \(X\) at each step in the rhmc trajectory, which
makes it an important bottleneck in the rhmc algorithm. Hence it is important
that the inversion carried out in that equation is very fast.
The inverter is a [conjugate gradient](../05_modules/inverter.md), with possibilities
for multiple RHS and multiple shifts to boost performance.
The [integrator](../05_modules/integrator.md) uses a leapfrog by default, but it
can use an Omelyan on the largest scale.
We use the HISQ/tree action, which is a tree-level improved
Lüscher-Weisz action in the gauge sector. The relative
weights of the plaquette and rectangle terms are
To use the rhmc class, the user will only have to call the constructor and two functions
```C++
rhmc(RhmcParameters rhmc_param, Gaugefield<floatT,onDevice,All,HaloDepth> &gaugeField, uint4* rand_state)
void init_ratapprox(RationalCoeff rat)
int update(bool metro=true, bool reverse=false)
```
The constructor has to be called with the usual template arguments and passed
an instance of RhmcParameters, the gauge field, and an uint4 array with
the RNG state. The function init_ratapprox will set the coefficients for
all the rational approximations and has to be called before updating.
The function update will generate one molecular dynamics (MD) trajectory.
If no arguments are passed to update() it will also perform a Metropolis
step after the trajectory. The Metropolis step can
also be omitted by passing the argument false to update. This is handy in
the beginning of thermalization. The second argument is false by default;
passing true to update will make the updater run the trajectory forward
and backwards for testing if the integration is reversible.
## Generating rational approximation input files
A tool to generate rational approximation files written by K. Clark
can be found in SIMULATeQCD/src/tools/rational_approx. One can find a document
explainRatApprox.pdf by Q. Yuan explaining the idea behind the rational
approximation for the fermion determinant, as well as some of the following
notation. The makefile makeRatApprox will compile the executable ratApprox,
which will generate you a rational approximation file for use with the
rhmc of SIMULATeQCD. You can call it with
```shell
ratApprox input.dat out.rational
```
The input file input.dat should be structured as
```C
npff // Number of pseudo-fermion flavors
y1
y2
mprec // Pre-conditioner mass (reduces the condition number in CG)
mq
order1
order2
lambda_low
lambda_high
precision
```
One block will generate three rational approximations according to
\(f(x) = x^{y_1/8} (x+ m_\text{prec}^2 -m_q^2 )^{y_2/8}\)
\(g(x) = x^{-y_1/8} (x+ m_\text{prec}^2 -m_q^2 )^{-y_2/8}\)
\(h(x) = x^{-y_1/4} (x+ m_\text{prec}^2 -m_q^2 )^{-y_2/4}\)
with \(m^2 = m_\text{prec}^2 - m_q^2\).
An example input file for \(N_f=2+1\) with standard Hasenbusch preconditioning
for the light flavors is given in exampleInput.dat.
This input file has a light block and a strange block, which together
generate six rational approximations. The light approximations are
The update routine saves a copy of the gauge field, applies the smearing to
the gauge field, builds the pseudo-fermion fields, generates the conjugate
momenta and calculates the Hamiltonian.
Then it starts the MD evolution by calling integrate() from the integrator
class (the integrator object is instantiated by the rhmc constructor). After
the MD trajectory the new Hamiltonian is calculated and - depending on the
arguments - the Metropolis step is done.
You can learn more about the smearing [here](../05_modules/gaugeSmearing.md).
When you want to use multiple pseudo-fermion fields, set no_pf in the rhmc
input file to the respective number. Be aware that this changes the way you
have to construct your ratapprox: In the remez input.dat, if you want to
generate \(N_f\) flavors using \(N_{pf}\) pseudo-fermion fields, you have to use \(N_f/N_{pf}\)
as an input, (which is then used \(N_{pf}\) times). Note that \(N_f/N_{pf}\) must be < 4.
no_pf is 1 per default.
The rhmc has the option to generate HISQ lattices with \(\mu_B=i\mu_I\). This can be accomplished
by setting the rhmc parameter mu0 to your desired value. (The default value is 0.)
This can be accomplished straightforwardly in lattice calculations by multiplying time-facing
staggered phases by an appropriate function of \(\mu_I\); see for instance
[this work](https://doi.org/10.1016/0370-2693(83)91290-X).
In our code we implement the imaginary phase corresponding to the chemical potential
chmp in staggeredPhases.h as:
```C++
img_chmp.cREAL = cos(chmp);
img_chmp.cIMAG = sin(chmp);
`