This directory contains code for generating lattices using rational function
approximations for the fermion determinants.  You can compile either the
RHMC algorithm (Rational Hybrid Monte Carlo) or the "phi" algorithm", where
the difference is that the RHMC includes a Metropolis accept/reject step
at the end of each trajectory.

There are several possible targets for "make", among them are these

su3_rhmc
su3_rhmd
su3_leapfrog
su3_omelyan_rhmc
su3_rhmc_hisq
su3_rhmd_hisq
su3_spectrum_hisq

The others are for testing and development.

The RHMC codes (and the su3_leapfrog code) implement the Asqtad RHMC
algorithm (with the accept/reject decision).  The "leapfrog" and
"omelyan" codes call for the leapfrog or Omelyan integrator,
respectively.  The su3_rhmc code is generic and expects the
integration algorithm to be specified in the Makefile as the INT_ALG
parameter.  The rhmd codes do not reject any trajectories, but they
still report the change in the action at the end of the trajectory.)

The "hisq" codes implement the HISQ action.  Please note that this
code is still under development, so it is not recommended for
production.  The su3_spectrum_hisq code does RHMC but builds with
spectrum measurement procedures and takes extra input parameters to
specify the spectrum measurement.  The su3_spectrum_hisq code can also
do updating, and uses the RHMC algorithm.

The Omelyan integrator stepsize is normalized so that there is one
force computation in each step of size "epsilon".  This differs by a
factor of two from the normalization in Takaishi and deForcrand, but
makes algorithm comparison more sensible.  In particular, when
lambda=1, it reduces to the leapfrog algorithm.

The "rhmd" algorithms should work in single precision, but I recommend
double precision for the "rhmc" algorithms on larger lattices.

The code currently allows up to 10 rational function pseudofermion
fields.  For example, you may want to use two separate rational
function approximations for det(M_light)^(-1/2) and
det(M_strange)^(-1/4).  Alternatively (and probably better) you may
just want a single rational function approximation for
det(M_light)^(-1/2)*det(M_strange)^(-1/4)

The rational function parameters are constructed with the poly4
utility in the remez-milc subdirectory.  Instructions are given below.
That file is then specified in the load_rhmc_params line in the input
parameter file, and it becomes input to the code at run time.

To run the code, use the same input file as for our traditional
su3_hmc.  MAKE SURE that the masses, and numbers of flavors are
consistent with the choices you made in preparing the rational
function parameter file!!  NO CHECKING IS DONE FOR THIS.  The flavor
numbers are used only in the experimental gauge action
"symanzik_1loop_nf_action.h".

You need three rational function approximations for each pseudofermion
field.  For the fermion force in the updating, you need an
approximation for 1/(powers of determinants).  For example, for a
determinant for two light quarks, det(m)^(-2/4).  For the heatbath
updating of the phi field, you need (powers of determinants)^(1/2) For
the computation of the fermion action, you need (powers of
determinants)^(-1/2)

SAMPLE INPUT PARAMETER FILE

Please see in.sample.su3_leapfrog.1 for an annotated explanation of
the input parameters.

USING poly4 TO GENERATE THE RATIONAL FUNCTION FILE

To generate rational function approximations, use a modified version of
Mike Clark's code (requires the gnu multiple precision library)

At the moment, for each pseudofermion field, this code can handle a
product of up to four powers of determinants with different masses
(just in case you want to put in dynamical charm and isospin breaking
at the same time).  The code "poly4" takes a small input parameter
file that specifies the pseudofermion fields.  See examples in the
remez-milc directory.  The first line of the file specifies the number
of such fields.  The remaining lines contain the following parameters:

n1 m1  n2 m2   n3 m3  n4 m4  md_order action_order min_eig max_eig digits

where n1=number of flavors of first type quarks
      m1=mass of first type quark
      next 6 arguments are for the other three quark types
      md_order = order of approximation for molecular dynamics integration
      action_order = order of approximation for action evaluation
      min_eig = minimum eigenvalue of (D_slash)^2.  Zero doesn't work, but 1e-15 does.
      max_eig = maximum eigenvalue of (D_slash)^2.  Actually, the maximum you ever
	expect to encounter is OK, and for equilibrated lattices with u_0 tuned
        correctly, 90 seems to work.
      digits is precision used internally.  Use lots (50 or so).

an example may be useful.   For two degenerate quarks with mass 0.005:
  ./remez-milc/poly4   2 .005    0 .99    0 .99    0 .99   \
	8 10   1e-15  90   50
will generate the rational function approximation to
   (M^dagger M + 4*.005^2)^(-1/2) for the molecular dynamics
   (M^dagger M + 4*.005^2)^(+1/4) for heat bath fermion update.
   (M^dagger M + 4*.005^2)^(+1/4) for fermion action computation.

The output format is more or less what is needed for the "rationals.h"
include file

For a second example, to make a single rational function approximation including
two light quarks (mass .005) and one strange quark (mass .05), use
  ./remez-milc/poly4   2 .005    1 .05    0 .99    0 .99   \
	8 10   1e-15  90   50
this will generate the rational function approximation to
   (M^dagger M + 4*.005^2)^(-1/2) * (M^dagger M + 4*.05^2)^(-1/4)

(*) a special case is 4 flavors, where the expansion of
(M^dagger M + 4m^2) ^ (-nf/4) is trivial
Use A = {0.0,1.0}, B={0.0,4*m^2}

The output from poly4 is then an input file for the RHMC code.

-- C. DeTar 04/04/08.
