Molecular dynamics for full QCD simulations with an improved action

TL;DR

This paper develops a molecular dynamics approach for full lattice QCD simulations using clover quarks, incorporating an improved action and advanced preconditioning to enhance computational efficiency.

Contribution

It introduces a new derivation of the equation of motion and extends even-odd preconditioning techniques for more efficient full QCD simulations.

Findings

Derived the equation of motion for molecular dynamics in lattice QCD with clover quarks.

Enhanced even-odd preconditioning methods to reduce computational effort.

Facilitated more efficient full lattice QCD simulations.

Abstract

I derive the equation of motion in molecular dynamics for doing full lattice QCD simulations with clover quarks. The even-odd preconditioning technique, expected to significantly reduce the computational effort, is further developed for the simulations.

Equations

$$S_G=-\frac{\beta }{N_c}\sum _{p}Retr(U_p)$$

$$=-\frac{\beta }{N_c}\sum _{x,\mu >\nu }Retr[U_{\mu }(x)U_{\nu }(x+\mu )U_{\mu }(x+\nu {)}^{†}{U}_{\nu }(x{)}^{†}]$$

$$S_F=\sum _{x,y}\bar{\psi }(x)M_{xy}\psi (y)=\sum _{x,y}\bar{\psi }(x)(A_{x,y}-\kappa B_{x,y})\psi (y),$$

$$A_{x,y}={\delta }_{x,y}[1-\frac{\kappa C}{2}\sum _{\mu ,\nu }{\sigma }_{\mu \nu }{\mathcal{F}}_{\mu \nu }(x)],$$

$$B_{x,y}=\sum _{\mu =1}^4(1-{\gamma }_{\mu })U_{\mu }(x){\delta }_{x,y-\mu }+(1+{\gamma }_{\mu })U_{\mu }^{†}(x-\mu ){\delta }_{x,y+\mu }.$$

$$Q_{\mu \nu }(x)=\frac{1}{4}[U_{\mu }(x)U_{\nu }(x+\mu )U_{\mu }(x+\nu {)}^{†}{U}_{\nu }(x{)}^{†}$$

$$+U_{\nu }(x)U_{\mu }(x-\mu +\nu {)}^{†}{U}_{\nu }(x-\mu {)}^{†}{U}_{\mu }(x-\mu )$$

$$+U_{\mu }(x-\mu {)}^{†}{U}_{\nu }(x-\mu -\nu {)}^{†}{U}_{\mu }(x-\mu -\nu )U_{\nu }(x-\nu )$$

$$+U_{\nu }(x-\nu {)}^{†}{U}_{\mu }(x-\nu )U_{\nu }(x+\mu -\nu )U_{\mu }(x{)}^{†}]$$

$$\displaystyle M=\left(\begin{array}[]{cc}A_{ee}&-\kappa B_{eo}\\ -\kappa B_{oe}&A_{oo}\end{array}\right),$$

$$\det (M)=\det (A_{oo})\det (A_{ee}-{\kappa }^2{B}_{eo}{A}_{oo}^{-1}{B}_{oe})=\det A_{oo}\det M_{ee}$$

$$M_{ee}=A_{ee}-{\kappa }^2{B}_{eo}{A}_{oo}^{-1}{B}_{oe},$$

$$\det (M^{†}M)=\det (A_{oo}^{†}{A}_{oo})\det (M_{ee}^{†}{M}_{ee})=\int d{\eta }_o^{†}d{\eta }_o\int d{\varphi }_e^{†}d{\varphi }_{e}exp(-S_{pf})$$

$$=\int d{\eta }_o^{†}d{\eta }_{o}exp[-{\eta }_o^{†}(A_{oo}^{†}{A}_{oo}{)}^{-1}{\eta }_o]\int d{\varphi }_e^{†}d{\varphi }_{e}exp[-{\varphi }_e^{†}(M_{ee}^{†}{M}_{ee}{)}^{-1}{\varphi }_e],$$

$$S_{pf}={\eta }_o^{†}(A_{oo}^{†}{A}_{oo}{)}^{-1}{\eta }_o+{\varphi }_e^{†}(M_{ee}^{†}{M}_{ee}{)}^{-1}{\varphi }_e,$$

$$D_i=A_{ii}-\sum _{k=1}^{i-1}{L}_{ik}{D}_k{L}_{ik}^{\ast },$$

$$L_{ij}{D}_j=A_{ij}-\sum _{k=1}^{j-1}{L}_{ik}{D}_k{L}_{jk}^{\ast },(j=1,...,i-1).$$

$$Y_i=b_i-\sum _{k=1}^{i-1}{L}_{ik}{Y}_k,(i=1,...,n),$$

$$X_i=Y_i/D_i-\sum _{k=i+1}^n{L}_{ki}^{\ast }{X}_k,(i=n,...,1),$$

$$\frac{d{\mathcal{A}}_{\mu }(x)}{d\tau }=P_{{\mathcal{A}}_{\mu }(x)},\frac{d{P}_{{\mathcal{A}}_{\mu }(x)}}{d\tau }=\frac{\partial H}{\partial {\mathcal{A}}_{\mu }(x)}=-\frac{\partial S_G}{\partial {\mathcal{A}}_{\mu }(x)}-\frac{\partial S_{pf}}{\partial {\mathcal{A}}_{\mu }(x)}.$$

$$-\frac{\partial S_G}{\partial {\mathcal{A}}_{\mu }(x)}=\frac{\beta }{2N_c}[\frac{\partial U_{\mu }(x)}{\partial {\mathcal{A}}_{\mu }(x)}STAPL{E}_{\mu }(x)+STAPL{E}_{\mu }(x{)}^{†}\frac{\partial U_{\mu }(x{)}^{†}}{\partial {\mathcal{A}}_{\mu }(x)}]$$

$$=\frac{i\beta }{2N_c}[U_{\mu }(x)STAPL{E}_{\mu }(x)-h.c.],$$

$$-\frac{\partial S_{pf}}{\partial {\mathcal{A}}_{\mu }(x)}=X_o^{A†}\frac{\partial (A_{oo}^{†}{A}_{oo})}{\partial {\mathcal{A}}_{\mu }(x)}{X}_o^A+X_e^{†}\frac{\partial (M_{ee}^{†}{M}_{ee})}{\partial {\mathcal{A}}_{\mu }(x)}{X}_e,$$

$$X_o^A=(A_{oo}^{†}{A}_{oo}{)}^{-1}{\eta }_o,X_e=(M_{ee}^{†}{M}_{ee}{)}^{-1}{\varphi }_e.$$

$$Y_o^A=A_{oo}{X}_o^A,Y_e=M_{ee}{X}_e,$$

$$-\frac{\partial S_{pf}}{\partial {\mathcal{A}}_{\mu }(x)}=X_o^{A†}\frac{\partial A_{oo}^{†}}{\partial {\mathcal{A}}_{\mu }(x)}{Y}_o^A+Y_o^{A†}\frac{\partial A_{oo}}{\partial {\mathcal{A}}_{\mu }(x)}{X}_o^A$$

$$+X_e^{†}\frac{\partial M_{ee}^{†}}{\partial {\mathcal{A}}_{\mu }(x)}{Y}_e+Y_e^{†}\frac{\partial M_{ee}}{\partial {\mathcal{A}}_{\mu }(x)}{X}_e.$$

$$\frac{\partial M_{ee}}{\partial {\mathcal{A}}_{\mu }(x)}=\frac{\partial A_{ee}}{\partial {\mathcal{A}}_{\mu }(x)}-{\kappa }^2\frac{\partial B_{eo}}{\partial {\mathcal{A}}_{\mu }(x)}{A}_{oo}^{-1}{B}_{oe}$$

$$+{\kappa }^2{B}_{eo}{A}_{oo}^{-1}\frac{\partial A_{oo}}{\partial {\mathcal{A}}_{\mu }(x)}{A}_{oo}^{-1}{B}_{oe}-{\kappa }^2{B}_{oe}{A}_{oo}^{-1}\frac{\partial B_{oe}}{\partial {\mathcal{A}}_{\mu }(x)},$$

$$\frac{\partial M_{ee}^{†}}{\partial {\mathcal{A}}_{\mu }(x)}=\frac{\partial A_{ee}}{\partial {\mathcal{A}}_{\mu }(x)}-{\kappa }^2\frac{\partial B_{eo}^{†}}{\partial {\mathcal{A}}_{\mu }(x)}{A}_{oo}^{-1}{B}_{oe}^{†}$$

Full text

††thanks: Present Address. E-Mail: stslxq@zsu.edu.cn

Published in Computer Physics Communications 94 (1996) 119-127.

1 Introduction

The effects of dynamical quarks are important in QCD at finite temperature as well as in some phenomenological aspects at zero temperature. Unfortunately, the inclusion of dynamical quarks is the most demanding task in computer simulations of lattice QCD. The hybrid molecular dynamics or Hybrid Monte Carlo (HMC) methods have been developed into very efficient algorithms (maybe the most popular) for dynamical quarks. In these algorithms, the equation of motion is the essential ingredient. One has to derive the relevant equations before writing the programs for molecular dynamical simulations. For lattice QCD with staggered or Wilson fermions, these equations have already been available in the literature [1, 2, 3].

Lattice QCD has discretization errors due to the lattice spacing $a$. At intermediate bare coupling, corresponding to relatively large $a$, these systematic errors might sometimes be very severe for Wilson fermions due to the chiral symmetry breaking term. The current computers do not allow the calculations done for very small $a$, because to reduce $a$ implies to use a much larger lattice. Another way out is to use the improved fermionic actions. Recently, it has been shown that the use of the clover action [4, 5] can significantly reduce these finite cut-off errors. However, the calculations of the clover action are much more complicated than the standard Wilson action. To my knowledge, there has not been a simulation of lattice QCD with the clover action in the literature.

The purpose of this paper is to derive the equation of motion for full QCD simulations with the clover action. Because the fermionic matrix has to be inverted in each step of the molecular dynamics step, it is also challenging to devise efficient algorithms for preconditioning [2, 6] the fermionic matrix so that the inverse is easier to compute. For this reason, I also extend the even-odd preconditioning technique, previously used for quark propagator measurements, to the case of dynamical clover fermions.

2 Preconditioning

2.1 The action

The action of the theory is $S=S_{G}+S_{F}$, where

[TABLE]

[TABLE]

is the gauge action. The clover action for the quarks [4, 5] is

[TABLE]

[TABLE]

[TABLE]

$A$ is local and hermitian, and $B$ connects only the nearest neighbor sites. The field strength tensor on the lattice is defined by ${\mathcal{F}_{\mu\nu}}(x)=[Q_{\mu\nu}(x)-Q_{\mu\nu}^{\dagger}(x)]/2i$, where $Q_{\mu\nu}$ is the averaged sum of four plaquettes on the $\mu\nu$ plane with the lattice site $x$ as one corner. Each plaquette is the product of four link variables in the counterclockwise sense and begins with the link variable directed away from the site $x$ and ends with the link variable directed towards site $x$, i.e.,

[TABLE]

[TABLE]

[TABLE]

[TABLE]

as shown in Fig. 1. This operator is so chosen as for the maximal symmetry on the lattice. Most symbols in above equations are conventional, while the coefficient $C$ in (2) depends on the choice of improvement strategy: $C=1$ for tree level improvement, and $C=[Re~{}tr(U_{p})/N_{c}]^{-3/4}$ for the tadpole improvement [7].

2.2 Even-odd splitting

The lattice sites can be organized in an even-odd checkerboard and the even sites are numbered before the odd sites such that the fermionic matrix can be written as [8, 9]

[TABLE]

where $e$ or $o$ denotes even or odd site on the lattice.

Using such an arrangement, we obtain

[TABLE]

where

[TABLE]

which couples only to even sites of the lattice. Now $\det(M)$ on the whole lattice has been factorized as a product of the determinant of the local matrix $A$ on the odd lattice and that of $M_{ee}$ on the even lattice.

To calculate the fermionic determinant, it is useful to introduce the pseudo-scalar variables $\eta_{o}$ and $\phi_{e}$, so that

[TABLE]

[TABLE]

where $S_{pf}$ is the pseudo-fermionic action

[TABLE]

describing two flavor quarks with the same bare mass. Notice that $\eta_{o}=A^{\dagger}_{oo}\xi_{o}$ and $\phi_{e}=M^{\dagger}_{ee}\theta_{e}$ have no direct coupling, where $\xi_{o}$ and $\theta_{e}$ are Gaussian noises injected at beginning of each molecular dynamics trajectory and held fixed during each trajectory.

In the remaining text, I will use the above even-odd preconditioning to discuss the molecular dynamics.

2.3 Fermionic inversion

For quark propagator measurements and also in each molecular dynamics step, one has to calculate $M_{ee}^{-1}$ or $(M^{\dagger}_{ee}M_{ee})^{-1}$, which can be implemented using the standard techniques like minimum residue, conjugate gradient or stabilized biconjugate gradient algorithms. The advantage of the even-odd splitting, as can also be seen later, is that such inversion is implemented only on the even lattice. Furthermore, due to the factor $\kappa^{2}$ in (8), $M_{ee}$ is better conditioned than $M$.

For the inversion $A_{oo}^{-1}$ on each odd site, because it is completely local, we can use the $LDL^{\dagger}$ decomposition [9, 10] to solve it. Since it is a hermitian matrix, there exists a diagonal matrix $D$ and lower-triangular matrix $L$ such that $A=LDL^{\dagger}$. Denoting $i$ and $j$ as the color-spin indexes of $A$, then

[TABLE]

[TABLE]

We can also compute the solution of $AX=b$ by $y=L^{-1}b$ and $X=(L^{\dagger})^{-1}D^{-1}Y$, i.e.,

[TABLE]

[TABLE]

with $n=12$, which is the number of colors times the number of spins. The calculation is quite easy because there are $n^{2}$ multiplications and only $n$ divisions.

3 Molecular dynamics

3.1 Equation of motion

To develop the equation of motion for the gauge field ${\mathcal{A}_{\mu}}(x)$, one has to introduce a Hamiltonian $H=\sum_{x,\mu}P^{2}_{{\mathcal{A}_{\mu}}(x)}/2+S_{g}+S_{pf}$, with $P_{{\mathcal{A}_{\mu}}(x)}$ the canonical conjugate momentum defined by $P_{{\mathcal{A}_{\mu}}(x)}=\partial H/\partial(d{\mathcal{A}_{\mu}}(x)/d\tau)$, and $\tau$ being the fictious molecular dynamics time.

The gauge configurations are generated by solving the Hamiltonian equation of motion:

[TABLE]

For the gauge action, it is quite easy to show that

[TABLE]

[TABLE]

where $STAPLE_{\mu}(x)$ is the sum over six staples surrounding the link $U_{\mu}(x)$. For the fermionic part,

[TABLE]

where

[TABLE]

If we define two more variables

[TABLE]

(15) can be simply rewritten as

[TABLE]

[TABLE]

A straightforward computation leads to

[TABLE]

[TABLE]

[TABLE]

[TABLE]

By defining the following variables on the odd sites

[TABLE]

we have

[TABLE]

[TABLE]

[TABLE]

[TABLE]

We can further demonstrate that the last two terms of these two equations in (21) are summarized as

[TABLE]

[TABLE]

where

[TABLE]

the same form as the fermionic force in the Wilson fermion case. As can also be seen later, the last two terms in (18) have exactly the same form for even and odd sites. Therefore, the introduction of the variables (16), (17) and (20) has a great advantage.

A remark has to be made: to keep the conjugate momentum traceless, the right hand side (r.h.s.) of (13) should finally be subtracted by a term being the trace of the r.h.s. divided by $N_{c}$.

3.2 Practical implementation

The simulations should be carried out in the first two steps as follows.

1. Generating the full configurations. In numerical integration of the equation of motion, one has to Taylor expand $U(\tau+d\tau)=exp[id\tau P(\tau)]U(\tau)=U(\tau)+id\tau P(\tau)U(\tau)+...$, $P(\tau+d\tau)=P(\tau)+d\tau{dP\over d\tau}(\tau)+...$, with finite order truncation. The leapfrog scheme can reduce the truncation errors to $N_{md}O(d\tau^{3})=O(d\tau^{2})$ at $N_{md}$ molecular dynamical steps. These errors can be canceled by a Metropolis test at the end of the trajectory. Of course, this $d\tau$ has to be fine tuned to maintain high acceptance rate and small auto-correlation time.

2. About the clover coefficient. For the tree level improvement, $C=1$. For the tadpole improvement scheme, the value of $C$ depends dynamically on the configurations. One has to determine $C$ self-consistently from the simulation. For example, one may first have an initial guess for it, then generate a gauge configuration. From the plaquette, we can get a new $C$ value. After some iteration, $C$ might converge to some stable value for some given $\beta$ and $\kappa$. (This could be done very quickly since the plaquette can be accurately measured with a small number of configurations and on small lattices, provided the system is far away from a phase transition).

3. Measuring the physical quantities. To obtain the improved hadronic matrix elements, rotation of quark fields [5, 8] is necessary.

4 More details about the fermionic force

I have described that how the introduction of the variables (16), (17) and (20) leads to the equation of motion similar to that of the Wilson case. What is different is the terms with matrix $A$, which makes the equation of motion much more complicate. Therefore, it deserves further discussions. Note the term in the pseudo-fermion action $Y_{o}^{A\dagger}A_{oo}X_{o}^{A}$ (also the second term) is placed only on odd sites of the lattice, then for $x$ being odd sites, there are terms only on $x$, $x+\mu+\nu$ and $x+\mu-\nu$ relevant for $Y_{o}^{A\dagger}(\partial A_{oo}/\partial{\mathcal{A}_{\mu}}(x))X_{o}^{A}$ as shown in Fig. 2, i.e.,

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

Here $\sum_{\nu}^{\prime}$ means the sum over $\nu\neq\mu$. For $x$ being even sites, there are terms only on $x+\mu$, $x+\nu$ and $x-\nu$ as shown in Fig. 3, i.e.,

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

Therefore for odd-site $x$, the first two terms in (18) read

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

Similarly, for even $x$, the first two terms in (18) are

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

Equations (26) and (27) can be generalized to terms with $y$ being odd or even:

[TABLE]

[TABLE]

[TABLE]

These relations are also quite useful when deriving the first two terms in (21).

Summarizing (22) and (28), (18) becomes

[TABLE]

[TABLE]

[TABLE]

Note that the difference in the form of the molecular dynamics equation on even and odd sites is only in the first term.

5 Discussions

In this paper, I have derived (14) and (29), relevant for equation of motion (13) in molecular dynamics (or HMC) simulations with clover fermions. I have also extended the even-odd precondition technique, previously introduced for the quark propagator measurements, to the case of simulations with dynamical clover fermions. With the preconditioning technique, the number of iterations required would be reduced by a factor of $3$ according to experience, and the most expensive part of the fermionic inversion is performed only on half lattice size. Therefore, it is expected that the preconditioned equation of motion would lead to considerable improvement over the unpreconditioned one. This scheme is vectorizable and has been parallelized. Currently, the simulations of QCD at finite temperature using the clover action are being carrying out on the Quadrics-APE100 parallel computers. The above work might lay a foundation of further computer simulations using dynamical clover fermions.

Of course, to obtain physical results from the numerical simulations, there are still of lot of work to do. For instance, because the clover constant $C$ depends dynamically on the gauge configuration, then there would be delicate interplay between $C$ and $\beta$, $\kappa$, in particular when the system is at criticality. This is a new subject beyond the scope of the paper.

It has been mentioned in [8] that even for the quenched clover propagator calculations, each minimum residue iteration took $35\%$ longer than for the Wilson action. One has to choose a more efficient algorithm for fermionic inversion because a good algorithm is critical for full simulations. It is known that the stabilized biconjugate gradient is more efficient than the minimum residue or conjugate gradient, reducing the CPU time by at least a factor of 1.5.

Even if with the clover action there is an improvement of the finite cut-off error, the simulations with this action require larger statistical samples, more arithmetic operations and much memory. Concerning the feasibility of a full QCD clover simulation on supercomputers, it is not easy to quantify, because it is machine and code dependent.

I hope to discuss these problems and report the physical results in the near future.

I am grateful to R. Horsley for valuable discussions on TAO (the programming language of Quadrics-APE100), P. Lepage on his improved perturbation theory, H. Shanahan for reference [9], and D. Richards and some members of UKQCD collaboration for useful conversations about $A=LDL^{\dagger}$ decomposition at Lattice 95. This work is sponsored by DESY.