Rational Hybrid Monte Carlo

The rational hybrid Monte Carlo (RHMC) algorithm is a way of updating gauge fields when simulating dynamical staggered fermions. In particular, the system’s Hamiltonian is considered as a function of MC time \(\tau\). It then evolves according to Hamilton’s equations of motion, and integrating this differential equation along \(\tau\) is a strategy to suggest a new configuration. This is called molecular dynamics (MD), and its progress is measured in trajectories, i.e. in the number of integrations. After a trajectory, the configuration is accepted/rejected with a standard Metropolis step, which is an example of a hybrid Monte Carlo algorithm. When working with staggered fermions, one generally removes the four-fold degeneracy of fermion d.o.f. by taking the fourth root of the Dirac operator. In practice, this fourth root is implemented as a rational approximation. Using a rational approximation for the Hamiltonian of the MD is what is meant by RHMC. A thorough overview of HMC can be found in, e.g. Gattringer and Lang. A scheme for doing HMC with staggered fermions was first made explicit here.

The trajectory is pushed by a fictitious driving force. There come contributions from the gauge part of the action as well as the fermion part. Schematically, it comes out to something of the form \( F \sim -\phi^\dagger \left(D D^{\dagger}\right)^{-1} \left(\partial D D^{\dagger}\right) \left(D D^{\dagger}\right)^{-1}\phi, \) where \(\phi\) is a pseudo-fermion field, and \(\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, with possibilities for multiple RHS and multiple shifts to boost performance. The integrator 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

\( c_\text{plaq} = 5/4, \)

\( c_\text{rect} = -1/6. \)

The fermion force is described in more detail in the documentation here. There, we list the coefficients of the various link constructs needed for the force and smearing.

To use the RHMC class, the user will only have to call the constructor and two functions

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

ratApprox input.dat out.rational 

The input file input.dat should be structured as

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

\( f(x) = x^{1/4} (x + m_s^2 - m_l^2 )^{-1/4} \)

\( g(x) = x^{-1/4} (x + m_s^2 - m_l^2 )^{1/4} \)

\( h(x) = x^{-1/2} (x + m_s^2 + m_l^2 )^{-1/2} \)

while the strange are

\( f(x) = x^{3/8} \)

\( g(x) = x^{-3/8} \)

\( h(x) = x^{-3/4} \)

Update

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 force here.

Multiple pseudo-fermions

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.

Imaginary chemical potential

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.

In our code we implement the imaginary phase corresponding to the chemical potential chmp in staggeredPhases.h as:

img_chmp.cREAL = cos(chmp);
img_chmp.cIMAG = sin(chmp);