Rational Hybrid Monte Carlo

The rational hybrid Monte Carlo (rhmc) 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

\( 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 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, 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. \)

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 smearing 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);