\section{Introduction}

In two recent papers\cite{nvu1,nvu2} (henceforth: Papers I and II) molecular dynamics at constant potential energy was introduced, tested, and
compared to well-known molecular dynamics algorithms. This new molecular dynamics is entitled \textit{NVU} dynamics in analogy to standard energy-conserving Newtonian \textit{NVE} dynamics.
The conclusion was that \textit{NVU} dynamics is a fully valid molecular dynamics, which for sufficiently large systems can be used
interchangeably with \textit{NVE} dynamics for calculating most quantities of interest. \textit{NVU} dynamics
is not faster than standard \textit{NVE} or \textit{NVT} dynamics,
but introduces a new way of thinking about molecular dynamics.
Molecular dynamics at constant potential energy was previously considered by Cotterill and co-workers\cite{cotterill1,cotterill4,cotterill3,cotterill2}, by
Scala \textit{et al.}\cite{scala}, and most recently by
Stratt and co-workers\cite{stratt1,stratt2,stratt3,stratt4}, who actually allowed also lower potential energy values. Our motivation for studying \textit{NVU} dynamics derived from recent work on
strongly correlating liquids and their isomorphs\cite{paper1,paper2,paper3,paper4,paper5,moleculesisomorphs,thermoscl,PRX} (see the Introduction of Paper I).

\textit{NVU} dynamics is defined by geodesic motion on the constant-potential-energy hypersurface $\Omega$ defined by

\begin{equation}\label{omega_def}
  \Omega = \{ \textbf{R} \in R^{3N} \,\ | \,\ U(\textbf{R})=U_0 \}.
\end{equation}
Here $\textbf{R} \equiv \{\textbf{r}^{(1)},...,\textbf{r}^{(N)}\}$ in which $\textbf{r}^{(k)}$ is the position vector of the $k$'te particle (we follow here the notation
of the Appendix of Paper II), and $U$ is the potential-energy function of
an $N$-particle classical system. A geodesic on $\Omega$ is a curve that satisfies the condition of stationary length for fixed endpoints $\textbf{R}_{A}$ and $\textbf{R}_{B}$, i.e.,

\begin{equation}\label{startcond}
  \delta \int_{\textbf{R}_{A}}^{\textbf{R}_{B}} dl\,\ \Bigg|_{\Omega} = 0,
\end{equation}
where $dl$ is the line element of the metric. The shortest path between any two points is a geodesic. On a sphere geodesics are great circles, the ''straightest lines'' of
the surface. Traversing a geodesic at constant velocity thus corresponds to a generalization of Newton's first law to a curved space (the
surface itself).

In Paper I the \textit{NVU} algorithm was developed via a discretization of Eq. (\ref{startcond}), subsequently carrying out the variation. This technique, which is known as variational
integration\cite{var1,var2,var3,var4}, resulted in a ''basic'' \textit{NVU} algorithm that is similar to the well-known Verlet
algorithm ${\bf R}_{i+1} \,=\, 2{\bf R}_i - {\bf R}_{i-1} + (\Delta t)^{2}{\bf F}_i/m$ for
Newtonian (\textit{NVE}) dynamics ($m$ is the particle mass which is assumed identical in this section, and $\textbf{F}_{i} \equiv -\nabla_{\textbf{R}_{i}} U$ is the $3N$-dimensional
force vector); the index $i$ refers to step $i$ of the
integration sequence. In the Verlet algorithm $\Delta t$ is a fixed time step length. In comparison, the basic \textit{NVU} algorithm is given by (Paper I)

\begin{equation}\label{numgeo}
  {\bf R}_{i+1} \,=\, 2{\bf R}_i - {\bf R}_{i-1} + \frac{-2{\bf F}_{i}\cdot({\bf R}_{i} - {\bf R}_{i-1})}{{\bf F}_{i}^{2}}\,{\bf F}_i.
\end{equation}
If the number of particles $N$ increases, the relative variation of the term $-2{\bf F}_{i}\cdot({\bf R}_{i} - {\bf R}_{i-1})/{\bf F}_{i}^{2}$
decreases, and this is why equivalence with Newtonian \textit{NVE} dynamics is established in the thermodynamic limit. This equivalence should be understood in the sense that the relative deviations
between, for instance, \textit{NVE} and \textit{NVU} auto-correlation functions go to zero as $N \to \infty$.

Paper I additionally developed a ''stabilized'' version
of the basic \textit{NVU} algorithm to prevent accumulation of numerical errors. This version of the algorithm is given by (defining the position changes
$\boldsymbol{\Delta}_{i+1/2} \equiv \textbf{R}_{i+1} - \textbf{R}_{i}$)

\begin{align}
  \boldsymbol{\Delta}_{i+1/2} & = l_{0} \frac{\textbf{A}_{i-1/2}}{||\textbf{A}_{i-1/2}||},\label{stabil1} \\
  {\bf R}_{i+1} & = \textbf{R}_{i} + \boldsymbol{\Delta}_{i+1/2},\label{stabil2}
\end{align}
where $l_0$ is the step length and

\begin{equation}
  \textbf{A}_{i-1/2}= \frac{\boldsymbol{\Delta}_{i-1/2}+(-2{\bf F}_i\cdot\boldsymbol{\Delta}_{i-1/2}+U_{i-1}-U_0)}{{\bf F}_i^2}\,{\bf F}_i.\label{stabil3}
\end{equation}
All simulations in Papers I and II were performed with the stabilized algorithm. The basic algorithm was used, however, for
theoretical considerations.

In this article we extend the stabilized \textit{NVU} algorithm to deal
with simulations of molecular systems. Molecular systems are simulated by introducing rigid and/or flexible
bonds between the atoms in the modelling. Flexible bonds introduce merely an additional contribution to $U$, for instance, harmonic spring potentials. The basic \textit{NVU} algorithm
conserves the total potential energy and can readily simulate flexible bonds. The focus in this paper is thus on implementing rigid bonds
in the framework of \textit{NVU} dynamics.

Section \ref{rigidalgo} considers \textit{NVU} dynamics with rigid bonds. Introducing rigid bonds in the simulations leads to Lagrangian multipliers
in addition to those introduced in order to keep the potential energy constant (Paper I). Section \ref{rigidalgo} is fairly technical and easiest to read after reading Paper I.
Section \ref{rigidmodel} gives simulation and model details. Section \ref{testing} tests the rigid-bond \textit{NVU} algorithm, and Sec. \ref{sampling}
investigates the \textit{NVU} sampling properties by comparing the \textit{NVU}
results to Nosé-Hoover \textit{NVT} results\cite{nose,hoover} on three different systems:
 the asymmetric dumbbell model\cite{moleculeshidden}, Lewis-Wahnstr{\"o}m OTP\cite{otp2}, and rigid SPC/E water\cite{spce}.
Nosé-Hoover \textit{NVT} dynamics is known to give results equivalent to \textit{NVE} dynamics in the thermodynamic limit\cite{evans1985}, and we refer to these dynamics interchangeably
in the forthcoming sections. Finally, Sec. \ref{sum} concludes.

\section{Rigid-bond $NVU$ algorithm}\label{rigidalgo}

The rigid bonds\cite{shake,goldstein} introduce constraints among the particle coordinates of the system. Each constraint $\alpha = 1,...,G$ is of the form

\begin{equation}
  \sigma_{\alpha}(\textbf{R}) \equiv (\textbf{r}^{(k_{\alpha})}-\textbf{r}^{(l_{\alpha})})^{2} \equiv (\textbf{r}^{\alpha})^{2} = C_{\alpha}^{2};
\end{equation}
it expresses that the distance between particles $k_{\alpha}$ and $l_{\alpha}$ is a constant, $C_{\alpha}$.
In Papers I and II the integral of Eq. (\ref{startcond}) was merely restricted to the constant-potential-energy
hypersurface $\Omega$. Each rigid bond constraint introduces a function $\sigma_{\alpha}$ to be kept constant, and thus
the integral of Eq. (\ref{startcond}) is now further restricted to the sub-manifold $\omega$ of $\Omega$ where the bond constraints are satisfied,

\begin{equation}
  \omega = \{ \textbf{R} \in \Omega \,\ | \,\ \sigma_{\alpha}(\textbf{R}) = C_{\alpha}\,\ , \,\ \alpha = 1,...,G \}.
\end{equation}
If the bond constraints are independent, as assumed throughout the paper, $\omega$ is a (3N - G - 1)-dimensional compact Riemannian manifold.
The variational principle defining \textit{NVU} dynamics with rigid bonds is given by

\begin{equation}\label{condrigid}
  \delta \int_{\textbf{R}_{A}}^{\textbf{R}_{B}} dl = 0 \,\ \Bigg|_{\omega}.
\end{equation}
Most of Papers I and II dealt with the case of identical particle masses, but we wish here to develop a completely general molecular \textit{NVU} algorithm.
The line element $dl$ is defined by

\begin{equation}\label{metric}
  dl^{2} \equiv \sum_{k} \tilde m _k \big(d{\bf r}^{(k)}\big)^{2},
\end{equation}
where $\tilde{m}_{k}$ = $m_{k}/\langle m \rangle$ is the ''reduced'' mass of particle $k$.
Equation (\ref{metric}) is not the standard Euclidean line element, but a mass-weighted line element that goes back to Hertz\cite{hertz1,hertz2}.
We shall refer to this metric as the ''Hertzian'' metric. The point of this particular
metric is that it ensures equivalence between \textit{NVU} and \textit{NVE} dynamics for systems of atoms and molecules of varying mass.
In appendix \ref{hertz} we derive the variable-mass atomic \textit{NVU} algorithm applying the
Hertzian metric (correcting also a typo of the Appendix of Paper II).

Applying the variational integration technique to Eq. (\ref{condrigid}) gives

\begin{equation}\label{nvu_cond}
\delta \left( \sum_i \sqrt{\sum_k \tilde{m}_{k} \left( \textbf{r}_{i}^{(k)} - \textbf{r}_{i-1}^{(k)} \right )^2 } -\sum_i \lambda_i U( {\bf R}_i)
  + \sum_{i,\alpha}\Lambda_{\alpha i}\sigma_{\alpha}(\textbf{R}_{i})\right) \,=\,0\,.
\end{equation}
In Eq. (\ref{nvu_cond}) the path is divided into a number of discrete points and one Lagrangian multiplier $\Lambda_{\alpha i}$ is introduced for each constraint $\alpha$
at every point $i$.
Following standard notation for constraint molecular dynamics\cite{goldstein,shake}, the Lagrangian multipliers of the bond constraints are chosen with a positive sign.
As in Papers I and II we now make the Ansatz of constant step length $l_{0}$, i.e.,

\begin{equation}\label{ansatz}
  \sum_k \tilde{m}_{k} \left( \textbf{r}_{i}^{(k)} - \textbf{r}_{i-1}^{(k)} \right )^2 \equiv l_{0}^{2}.
\end{equation}
Carrying out the variation of Eq. (\ref{nvu_cond}) using Eq. (\ref{ansatz}) leads to (compare the derivation in Paper I)

\begin{equation}\label{some}
  \textbf{r}_{i+1}^{(k)} = 2\textbf{r}_{i}^{(k)} - \textbf{r}_{i-1}^{(k)} + \frac{l_{0}}{\tilde{m}_{k}}\lambda_{i}\textbf{f}_{i}^{(k)}
  + \frac{l_{0}}{\tilde{m}_{k}}\nabla_{\textbf{r}_{i}^{(k)}}\sum_{\alpha}\Lambda_{\alpha i}\sigma_{\alpha},
\end{equation}
where $\textbf{f}_{i}^{(k)} = -\nabla_{\textbf{r}_{i}^{(k)}}U$ is the force on particle $k$ at step $i$.
This equation constitutes the \textit{NVU} algorithm with rigid bonds. It has a close
resemblance to the Lagrangian equations of motion with holonomic constraints\cite{goldstein}, i.e., rigid-bond \textit{NVE} dynamics\cite{shake}. Equation (\ref{some}) contains $G+1$
Lagrangian multipliers for each integration step, which must be determined to complete the algorithm.

\subsection{Determining the \textit{NVU} Lagrangian multipliers}

This section shows how to calculate the Lagrangian multipliers. Since the algorithm is to be implemented on a computer (with finite-precision), we shall
proceed directly to a ''stabilized'' algorithm conserving for indefinitely long runs potential energy, bond lengths, and step length (in $3N$-dimensions).
The resulting algorithm reduces to the stabilized atomic \textit{NVU} algorithm of Eqs. (\ref{stabil1})-(\ref{stabil3}) in the case of no bonds constraints.

Some notation used in the following derivation is now introduced (the nomenclature of text is summarized in Table \ref{notationtab}).

\begin{table}[H]
  \centering
  \begin{tabular}
    {|| p{2.0cm} || p{14.0cm} ||}
    \hline
    \hline
    Symbol & Definition \\
    \hline
    \hline
    &  \\
    \hline
    $\sigma_{\alpha}(\textbf{R})$ & The $\alpha$'te bond constraint between particles $k_{\alpha}$ and $l_{\alpha}$ with $\alpha = 1, ..., G$. ($\sigma_{\alpha} = (\textbf{r}^{\alpha})^{2} = C_{\alpha}^{2}$). \\
    $\tilde{m}_{k}$ & The mass of particle $k$ divided by the average mass of the system. ($\tilde{m}_{k} = m_{k}/\langle m \rangle$). \\
    \hline
    & 3-dimensional vectors \\
    \hline
    $\textbf{r}^{(k)}_{i}$ & Position of particle $k$ at step $i$. \\
    $\boldsymbol{\delta}_{i+1/2}^{(k)}$ & Displacement of the position of particle $k$ between step $i$ and $i+1$. ($\boldsymbol{\delta}_{i+1/2}^{(k)} = \textbf{r}_{i+1}^{(k)}
    - \textbf{r}_{i}^{(k)}$).\\
    $\textbf{f}^{(k)}_{i}$ & Force on particle $k$ at step $i$. ($\textbf{f}^{(k)}_{i} = - \nabla_{\textbf{r}^{(k)}_{i}}U$). \\
    $\textbf{g}^{(k)}_{i}$ & Constraint force on particle $k$ at step $i$. ($\textbf{g}_{i}^{(k)} = \nabla_{\textbf{r}_{i}^{(k)}}\sum_{\alpha}\Lambda_{\alpha i}\sigma_{\alpha}$). \\
    $\textbf{r}_{i}^{\alpha}$ & Displacement of the positions of particles $k_{\alpha}$ and $l_{\alpha}$ at step $i$.
    ($\textbf{r}_{i}^{\alpha} = \textbf{r}_{i}^{(k_{\alpha})} - \textbf{r}_{i}^{(l_{\alpha})}$). \\
    $\boldsymbol{\delta}_{i-1/2}^{\alpha}$ & Displacement of the velocities of particles $k_{\alpha}$ and $l_{\alpha}$ at step $i-1/2$.
    ($\boldsymbol{\delta}_{i-1/2}^{\alpha} = \boldsymbol{\delta}^{(k_{\alpha})}_{i-1/2} - \boldsymbol{\delta}^{(l_{\alpha})}_{i-1/2}$).  \\
    $\textbf{s}_{i}^{\alpha}$ & Sum of displacements of positions and velocities of particles $k_{\alpha}$ and $l_{\alpha}$ at, respectively, step $i$ and $i-1/2$.
    $(\textbf{s}_{i}^{\alpha} = \textbf{r}_{i}^{\alpha} +  \boldsymbol{\delta}_{i-1/2}^{\alpha}$).  \\
    $\tilde{\textbf{f}}^{\alpha}_{i}$ & Displacement of the forces on particles $k_{\alpha}$ and $l_{\alpha}$ at step $i$ divided by their reduced particle mass.
    ($\tilde{\textbf{f}}_{i}^{\alpha} = \textbf{f}_{i}^{(k_{\alpha})}/\tilde{m}_{k_{\alpha}} - {\textbf{f}}_{i}^{(l_{\alpha})}/\tilde{m}_{l_{\alpha}}$). \\
    $\tilde{\textbf{g}}^{\alpha}_{i}$ & Displacement of the constraint forces on particles $k_{\alpha}$ and $l_{\alpha}$ at step $i$ divided by their reduced particle mass.
    ($\tilde{\boldsymbol{g}}^{\alpha}_{i} = \textbf{g}_{i}^{(k_{\alpha})}/\tilde{m}_{k_{\alpha}} - {\textbf{g}}_{i}^{(l_{\alpha})}/\tilde{m}_{l_{\alpha}}$). \\
    \hline
    & $3N$-dimensional vectors \\
    \hline
    $\textbf{R}_{i}$ & Position of all particles at step $i$. ($\textbf{R}_{i} = \{\textbf{r}_{i}^{(1)},...,\textbf{r}_{i}^{(N)}\}$). \\
    $\boldsymbol{\Delta}_{i+1/2}$  & Displacement  of the positions between step $i$ and $i+1$. ($\boldsymbol{\Delta}_{i+1/2} = \textbf{R}_{i+1} - \textbf{R}_{i}$). \\
    $\textbf{F}_{i}$ & Force on all particles at step $i$. ($\textbf{F}_{i} = - \nabla_{\textbf{R}_{i}}U$). \\
    $\tilde{\textbf{F}}_{i}$ & Force on all particles at step $i$ divided by the reduced particle mass. ($\tilde{\textbf{F}}_{i}
    = \{\textbf{f}_{i}^{(1)}/\tilde{m}_{1},...,\textbf{f}_{i}^{(N)}/\tilde{m}_{N}\}$).  \\
    $\tilde{\textbf{G}}_{i}$ & Constraint force on all particles at step $i$ divided by the reduced particle mass.
    ($\tilde{\textbf{G}}_{i} = \{\textbf{g}_{i}^{(1)}/\tilde{m}_{1},...,\textbf{g}_{i}^{(N)}/\tilde{m}_{N}\}$). \\
    \hline
    \hline
  \end{tabular}
  \caption{Definitions and nomenclature of the text.}
  \label{notationtab}
\end{table}
Defining $\boldsymbol{\delta}_{i+1/2}^{(k)} \equiv \textbf{r}_{i+1}^{(k)} - \textbf{r}_{i}^{(k)}$  and
$\textbf{g}_{i}^{(k)} \equiv \nabla_{\textbf{r}_{i}^{(k)}}\sum_{\alpha}\Lambda_{\alpha i}\sigma_{\alpha}$ the
''Leap-frog''\cite{tildesley} version of the rigid-bond \textit{NVU} algorithm Eq. (\ref{some}) reads

\begin{align}
  \boldsymbol{\delta}_{i+1/2}^{(k)} & = \boldsymbol{\delta}_{i-1/2}^{(k)}  + \frac{l_{0}}{\tilde{m}_{k}}\lambda_{i}\textbf{f}_{i}^{(k)} +
  \frac{l_{0}}{\tilde{m}_{k}}\textbf{g}_{i}^{(k)},\label{LF1} \\
  \textbf{r}_{i+1}^{(k)} & = \textbf{r}_{i}^{(k)} + \boldsymbol{\delta}_{i+1/2}^{(k)}. \label{LF2}
\end{align}
In analogy to rigid-bond \textit{NVE} dynamics we call $\textbf{g}^{(k)}_{i}$ the ''constraint force'' on particle $k$ at step $i$. Introducing the
notation $\tilde{\textbf{F}}_{i} \equiv \{\textbf{f}_{i}^{(1)}/\tilde{m}_{1}, ..., \textbf{f}_{i}^{(N)}/\tilde{m}_{N}\}$ and
$\tilde{\textbf{G}}_{i} \equiv \{\textbf{g}^{(1)}_{i}/\tilde{m}_{1}, ..., \textbf{g}^{(N)}_{i}/\tilde{m}_{N}\}$,
the \textit{NVU} algorithm in the full $3N$-dimensional coordinate space reads

\begin{align}
  \boldsymbol{\Delta}_{i+1/2} & =  \boldsymbol{\Delta}_{i-1/2} + l_{0}\lambda_{i}\tilde{\textbf{F}}_{i} + l_{0}\tilde{\textbf{G}}_{i} \label{NOT1}, \\
  \textbf{R}_{i+1} & = \textbf{R}_{i} + \boldsymbol{\Delta}_{i+1/2},\label{NOT2}.
\end{align}
The Lagrangian multipliers are calculated by combining a result derived in Paper I with the method applied in the \textit{SHAKE} algorithm\cite{shake} for
rigid bonds in \textit{NVE} dynamics\cite{shake,toxconstraintnve,toxconstraintnph}. The \textit{SHAKE} algorithm
calculates the Lagrangian multipliers from the equations $(\textbf{r}_{i+1}^{\alpha})^{2} = C^{2}_{\alpha}$. In doing so, the target value of the constraints $C_{\alpha}$
appears explicitly in the algorithm, making the bond lengths insensitive to numerical error.
The expression for $\textbf{r}_{i+1}^{\alpha}$ is supplied by
the integration algorithm containing herein the Lagrangian multipliers. In our case, this gives $G$ equations with $G+1$ unknowns. The missing equation is supplied
by an expression derived in Paper I, namely that $U_{i+1}$ = $U_{i-1} - \textbf{F}_{i}\cdot(\textbf{R}_{i+1} - \textbf{R}_{i-1})$ to third order in the step length. In the discrete sequence of
points $U_{i+1}$ is set equal to $U_{0}$ (the constant defining $\Omega$), making the constraint of constant potential energy also insensitive to numerical error.
We thus have the following $G+1$ equations for calculating the Lagrangian multipliers

\begin{align}
  U_{i-1} - \textbf{F}_{i}\cdot(\textbf{R}_{i+1} - \textbf{R}_{i-1}) - U_{0}  & = 0 \label{trick}, \\
  (\textbf{r}_{i+1}^{\alpha})^{2} - C^{2}_{\alpha} & = 0,\,\ (\alpha = 1,...,G).\label{trick1}
\end{align}
By Eqs. (\ref{NOT1}) and (\ref{NOT2}); $\textbf{R}_{i+1} - \textbf{R}_{i-1}$ = $\boldsymbol{\Delta}_{i+1/2} + \boldsymbol{\Delta}_{i-1/2}$
= $2\boldsymbol{\Delta}_{i-1/2} + l_{0}\lambda_{i}\tilde{\textbf{F}}_{i} + l_{0}\tilde{\textbf{G}}_{i}$.
Defining $\boldsymbol{\delta}^{\alpha}_{i-1/2} \equiv \boldsymbol{\delta}^{(k_{\alpha})}_{i-1/2} - \boldsymbol{\delta}^{(l_{\alpha})}_{i-1/2}$,
$\tilde{\textbf{f}}_{i}^{\alpha} \equiv \textbf{f}_{i}^{(k_{\alpha})}/\tilde{m}_{k_{\alpha}} - {\textbf{f}}_{i}^{(l_{\alpha})}/\tilde{m}_{l_{\alpha}}$, and
$\tilde{\boldsymbol{g}}^{\alpha}_{i} \equiv \textbf{g}_{i}^{(k_{\alpha})}/\tilde{m}_{k_{\alpha}} - {\textbf{g}}_{i}^{(l_{\alpha})}/\tilde{m}_{l_{\alpha}}$,
since by Eqs. (\ref{LF1}) and (\ref{LF2}); $\textbf{r}^{\alpha}_{i+1}$
= $\textbf{r}^{(k_{\alpha})}_{i+1} - \textbf{r}^{(l_{\alpha})}_{i+1}$
= $\textbf{r}_{i}^{(k_{\alpha})} - \textbf{r}_{i}^{(l_{\alpha})} + \boldsymbol{\delta}_{i+1/2}^{(k_{\alpha})} - \boldsymbol{\delta}_{i+1/2}^{(l_{\alpha})}$
= $\textbf{r}^{\alpha}_{i} + \boldsymbol{\delta}^{\alpha}_{i-1/2} + l_{0}\lambda_{i}\tilde{\textbf{f}}_{i}^{\alpha}
  + l_{0}\tilde{\textbf{g}}_{i}^{\alpha}$, it follows that

\begin{align}
  U_{i-1} - \textbf{F}_{i}\cdot\Big[2\boldsymbol{\Delta}_{i-1/2} + l_{0}\lambda_{i}\tilde{\textbf{F}}_{i} + l_{0}\tilde{\textbf{G}}_{i} \Big] - U_{0} & = 0,\label{taylor1} \\
  \Big[ \textbf{r}^{\alpha}_{i} + \boldsymbol{\delta}^{\alpha}_{i-1/2} + l_{0}\lambda_{i}\tilde{\textbf{f}}_{i}^{\alpha}
  + l_{0}\tilde{\textbf{g}}_{i}^{\alpha} \Big]^{2} - C^{2}_{\alpha} & = 0,\,\ (\alpha = 1,...,G).\label{taylor2}
\end{align}
The above coupled quadratic equations for the Lagrangian multipliers are now solved following the produce of the \textit{MILC-SHAKE} algorithm\cite{milcshake}, which
starts by neglecting the second order terms in the Lagrangian multipliers and solving the resulting linear equations. Afterwards, the second order
terms are taken into account in an iterative manner - the details of which are described below.

For each integration step $i$, the linearized equations are given by

\begin{equation}\label{linear}
  \textbf{A}_{i}\boldsymbol{\lambda}_{i} = \textbf{b}_{i},
\end{equation}
where $\textbf{A}_{i}$ is a $(G+1) \times (G+1)$ matrix, $\boldsymbol{\lambda}_{i} \equiv \{\lambda_{i}, \Lambda_{1 i}, ..., \Lambda_{G i} \}$, and
$\textbf{b}_{i}$ a $G+1$ column vector. We start by calculating explicitly the first few elements of the matrix $\textbf{A}_{i}$.
$A_{11}$ consists merely of the factor infront of $\lambda_{i}$ in Eq. (\ref{taylor1}), i.e., $A_{11} = -l_{0}\tilde{\textbf{F}}_{i}\cdot\textbf{F}_{i}$.
The second element $A_{12}$ appears after expansion of the dot product $\textbf{F}_{i} \cdot \tilde{\textbf{G}}_{i}$.
Noting that $\nabla_{\textbf{r}_{i}^{(k_{\alpha})}}\sigma_{\alpha} = 2\textbf{r}_{i}^{\alpha}$, we have $\textbf{F}_{i} \cdot \tilde{\textbf{G}}_{i}$ =
$\textbf{f}_{i}^{(1)}\cdot \textbf{g}_{i}^{(1)}/\tilde{m}_{1} + ... + \textbf{f}_{i}^{(N)}\cdot \textbf{g}_{i}^{(N)}/\tilde{m}_{N}$
= 2$\Lambda_{1i}(\tilde{\textbf{f}}^{1}_{i}\cdot\textbf{r}_{i}^{1}) + ... + 2\Lambda_{Gi}(\tilde{\textbf{f}}^{G}_{i}\cdot\textbf{r}_{i}^{G})$. The last equation follows as the Lagrangian multipliers
appear in pairs, differing only by the sign from $\nabla_{\textbf{r}_{i}^{(k_{\alpha})}}\sigma_{\alpha}$ and the term $\textbf{f}_{i}^{(k_{\alpha})}/\tilde{m}_{k_{\alpha}}$. We thus
find $A_{12}$ = $-2l_{0}\tilde{\textbf{f}}^{1}_{i}\cdot\textbf{r}_{i}^{1}$, $A_{13} = -2l_{0}\tilde{\textbf{f}}^{2}_{i}\cdot\textbf{r}_{i}^{2}$, etc.
In the second row of $\textbf{A}_{i}$, the short-hand notation $\textbf{s}^{\alpha}_{i} \equiv \textbf{r}_{i}^{\alpha} + \boldsymbol{\delta}^{\alpha}_{i-1/2}$ is introduced, making
 $A_{21}$ = $2l_{0}(\textbf{s}^{1}_{i} \cdot \tilde{\textbf{f}}^{1}_{i})$, i.e., the factor infront of $\lambda_{i}$ after squaring of the parentheses. The next element $A_{22}$
appears after expanding $\textbf{s}_{i}^{1} \cdot \tilde{\textbf{g}}^{1}_{i}$
= $\boldsymbol{s}^{1}_{i}\cdot\sum_{\beta}\Lambda_{\beta i}(\frac{1}{\tilde{m}_{k_{1}}}\nabla_{\textbf{r}_{i}^{(k_{1})}}\sigma_{\beta} - \frac{1}{\tilde{m}_{l_{1}}}\nabla_{\textbf{r}_{i}^{(l_{1})}}\sigma_{\beta})$.
In this sum, we identify the factor in front of $\Lambda_{1i}$, giving $A_{22}$ = $2l_{0}\boldsymbol{s}^{1}_{i}\cdot(\frac{1}{\tilde{m}_{k_{1}}}\nabla_{\textbf{r}_{i}^{(k_{1})}}\sigma_{1}
- \frac{1}{\tilde{m}_{l_{1}}}\nabla_{\textbf{r}_{i}^{(l_{1})}}\sigma_{1})$, and similarly for the remaining elements of the second row.

Altogether, the elements of $\textbf{A}_{i}$ are thus given by
\begin{align}
  \textbf{A}_{i} & = 2l_{0}
  \begin{pmatrix}
    -\tilde{\textbf{F}}_{i}\cdot\textbf{F}_{i} / 2 & -\tilde{\textbf{f}}_{i}^{1}\cdot \textbf{r}^{1}_{i} & \cdots &  -\tilde{\textbf{f}}_{i}^{G}\cdot \textbf{r}^{G}_{i} \\
    \boldsymbol{s}^{1}_{i} \cdot \tilde{\textbf{f}}_{i}^{1} & \boldsymbol{s}^{1}_{i}\cdot(\frac{1}{\tilde{m}_{k_{1}}}\nabla_{\textbf{r}_{i}^{(k_{1})}}\sigma_{1}
    - \frac{1}{\tilde{m}_{l_{1}}}\nabla_{\textbf{r}_{i}^{(l_{1})}}\sigma_{1}) & \cdots &
    \boldsymbol{s}^{1}_{i}\cdot(\frac{1}{\tilde{m}_{k_{1}}}\nabla_{\textbf{r}_{i}^{(k_{1})}}\sigma_{G}
    - \frac{1}{\tilde{m}_{l_{1}}}\nabla_{\textbf{r}_{i}^{(l_{1})}}\sigma_{G})  \\
    \vdots  & \vdots  & \vdots & \vdots \\
    \boldsymbol{s}^{G}_{i} \cdot \tilde{\textbf{f}}_{i}^{G} & \boldsymbol{s}^{G}_{i}\cdot(\frac{1}{\tilde{m}_{k_{G}}}\nabla_{\textbf{r}_{i}^{(k_{G})}}\sigma_{1}
    - \frac{1}{\tilde{m}_{l_{G}}}\nabla_{\textbf{r}_{i}^{(l_{G})}}\sigma_{1}) & \cdots &
    \boldsymbol{s}^{G}_{i}\cdot(\frac{1}{\tilde{m}_{k_{G}}}\nabla_{\textbf{r}_{i}^{(k_{G})}}\sigma_{G}
    - \frac{1}{\tilde{m}_{l_{G}}}\nabla_{\textbf{r}_{i}^{(l_{G})}}\sigma_{G})
  \end{pmatrix}.
\end{align}
The column vector $\textbf{b}_{i}$ consists of all  zeroth-order terms in Eqs. (\ref{taylor1}) and (\ref{taylor2})

\begin{align}
  \textbf{b}_{i} & =
  \begin{pmatrix}
    U_{0} - U_{i-1} + 2\textbf{F}_{i}\cdot\boldsymbol{\Delta}_{i-1/2} \\
    C^{2}_{1} - (\textbf{s}^{1}_{i})^{2} \\
    \vdots \\
    C^{2}_{G} - (\textbf{s}^{G}_{i})^{2}
  \end{pmatrix}.
\end{align}
Turning now to the iteration procedure, the second-order terms in the Lagrangian multipliers (Eq. (\ref{taylor2})) are taken into account by iterating
the right-hand side of Eq. (\ref{linear}) via the scheme ($\alpha = 1, ..., G$)

\begin{equation}\label{ite}
  b_{\alpha}^{j+1} = b_{\alpha}^{j} + \Big[C_{\alpha}^{2} - \big((\textbf{r}_{i+1}^{\alpha})^{2} \big)^{j}\Big].
\end{equation}
The superscript $j$ refers here to iteration $j$, and $\big((\textbf{r}_{i+1}^{\alpha})^{2}\big)^{j}$ are the positions associated with
iteration $j$. The element $b_{0}$ is not updated as it derives from the
constraint of constant potential energy. For each iteration $j$ the term $C_{\alpha}^{2} - \big((\textbf{r}_{i+1}^{\alpha})^{2}\big)^{j}$ is expected to become smaller as the bonds
are satisfied better and better, and convergence was achieved within a few iterations\cite{milcshake}.

For each integration step $i$, the algorithm for determining the \textit{NVU} Lagrangian multipliers thus proceeds as follows

\begin{enumerate}
\item The Lagrangian multipliers of iteration $j$, $(\boldsymbol{\lambda}_{i})^{j}$, are calculated from Eq. (\ref{linear}).
\item $\big((\textbf{r}_{i+1}^{\alpha})^{2} \big)^{j}$ is calculated via Eqs. (\ref{LF1}) and (\ref{LF2}) using $(\boldsymbol{\lambda}_{i})^{j}$.
\item $\textbf{b}_{i}$ is updated via Eq. (\ref{ite}) from $\big((\textbf{r}_{i+1}^{\alpha})^{2} \big)^{j}$.
\item The above steps are repeated (starting iteration $j+1$) until convergence is established (we used a preset number of iterations, typically 3-5).
\end{enumerate}
How is constant step length $l_{0}$ ensured numerically? Generalizing the approach of Paper I we introduce a normalizing factor such that

\begin{align}
  \boldsymbol{\delta}_{i+1/2}^{(k)} & = l_{0}\frac{\boldsymbol{\chi}^{(k)}_{i-1/2}}{\sqrt{\sum_{k}\tilde{m}_{k}(\boldsymbol{\chi}_{i-1/2}^{(k)})^{2}}}, \\
  \textbf{r}_{i+1}^{(k)} & = \textbf{r}_{i}^{(k)} + \boldsymbol{\delta}_{i+1/2}^{(k)},
\end{align}
where

\begin{equation}
  \boldsymbol{\chi}_{i-1/2}^{(k)} \equiv \boldsymbol{\delta}_{i-1/2}^{(k)}  + \frac{l_{0}}{\tilde{m}_{k}}\lambda_{i}\textbf{f}_{i}^{(k)}
  +  \frac{l_{0}}{\tilde{m}_{k}}\textbf{g}_{i}^{(k)}.
\end{equation}
The normalizing factor is close to unity\cite{nvu1} and ensures trivially
$\sum_{k}\tilde{m}_{k}(\boldsymbol{\delta}_{i+1/2}^{(k)})^{2}$ = $l^{2}_{0}$, i.e., that the step length is conserved. The algorithm is now absolutely stable, conserving
potential energy, bond lengths, and step length for indefinitely long runs. The stability of the \textit{NVU} algorithm is tested numerically in Sec. \ref{testing}.

\subsection{Alternative determination of the \textit{NVU} Lagrangian multipliers}\label{alternative}

The previous section followed the traditional way of calculating the Lagrangian multipliers.
The \textit{NVU} Lagrangian multipliers may also be calculated by Taylor expanding the constraints $\sigma_{\alpha}$ in analogy to the method sketched above for the
potential energy. In this way, the constraints of constant potential energy and constant bond lengths are treated on equal footing.
The set of equations to be solved is the following (recall that
$\textbf{R}_{i+1} - \textbf{R}_{i-1}$ = $2\boldsymbol{\Delta}_{i-1/2} + l_{0}\lambda_{i}\tilde{\textbf{F}}_{i} + l_{0}\tilde{\textbf{G}}_{i}$)

\begin{align}
  U_{i-1} - \textbf{F}_{i} \cdot (\textbf{R}_{i+1} - \textbf{R}_{i-1}) - U_{0}  & = 0, \label{linear1} \\
  \sigma_{\alpha(i-1)} + \nabla_{\textbf{R}_{i}}\sigma_{\alpha i} \cdot (\textbf{R}_{i+1} - \textbf{R}_{i-1}) - C^{2}_{\alpha}  & = 0,\,\ (\alpha = 1,...,G). \label{linear2}
\end{align}
The determination of the Lagrangian multipliers is linear and thus no iterations are needed. The bond
constraints $\sigma_{\alpha}$ are obeyed to the same order $O(l_{0}^{3})$ as the constraint of constant potential energy. The sampling properties of this novel, alternative
determination method is tested briefly in Sec. \ref{sampling}. It appears to be a promising new way of determining the Lagrangian multipliers in connection with rigid bonds,
which might also be useful for standard bond-constraint \textit{NVE} or \textit{NVT} simulations.

\section{Simulation details and model systems}\label{rigidmodel}

We investigated three systems: The asymmetric dumbbell model, the Lewis-Wahnstr{\"o}m OTP model, and rigid SPC/E water. For all simulated pair potentials the shifted-force
truncation scheme was applied at a cut-off radius $r_{c}$. If the pair potential is $v(r)$ and the pair force is $f(r)=-v'(r)$, the shifted force is given by\cite{tildesley,FCS1}

\begin{equation}
  f_{\rm SF}(r)\,=\,
  \begin{cases}  f(r)-f(r_c) & \text{if}\,\, r<r_c\,, \\
    0 &\text{if}\,\, r>r_c\,.
  \end{cases}
\end{equation}
This corresponds to using the following pair potential below $r_c$: $v_{\rm SF}(r) = v(r) - v'(r_c) ( r-r_c) - v(r_c)$. All simulations were
performed with the \textit{NVT} and \textit{NVU} algorithms. Recall that \textit{NVE} and \textit{NVT} dynamics give equivalent results\cite{evans1985}; for this reason no
simulations are presented for \textit{NVE} dynamics.
The RUMD code\cite{rumd} was used for molecular dynamics simulations (an optimized open-source GPU code). The \textit{NVT} ensemble is
generated via the Nosé-Hoover algorithm\cite{nose,hoover,nvttoxvaerd}, and the bonds held fixed using the time-reversible constraint
algorithm of Refs. \onlinecite{toxconstraintnve,toxconstraintnph}. The \textit{NVU} algorithm is described in Sec. \ref{rigidalgo}. The starting files for \textit{NVU} dynamics
were taken from an equilibrated \textit{NVT} simulation. The positions and velocities of the \textit{NVT} configuration do not correspond perfectly
to motion on $\omega$, since the potential energy and step length are not those of $U_{0}$ and $l_{0}$, respectively. As all three constraints are to be satisfied
simultaneously, this results in numerical problems when starting the simulation from the particular \textit{NVT} configuration.
A more gentle procedure is thus applied, where the atomic \textit{NVU} algorithm is used for a couple of integration steps to ensure the values of $U_{0}$ and $l_{0}$. Afterwards, the
rigid-bond \textit{NVU} algorithm is used.

\subsection{\textit{NVU} iteration procedure}

The quadratic equations (Eq. (\ref{ite})) were iterated with a fixed number of iterations (between 3 and 5).
The linear systems were solved utilizing CUSP\cite{cusp}, a library for solving systems of linear equations on the \textit{GPU}. More specifically,
the stabilized biconjugate gradient algorithm with a Jacobi preconditioner\cite{recipes} was used with the initial value $\boldsymbol{\lambda}_{i} = 0$.
The relative tolerance $\tau$ for the asymmetric dumbbell and Lewis-Wahnstr{\"o}m OTP models was chosen as $\tau = 10^{-7}$ and for rigid SPC/E water as $\tau = 3 \cdot 10^{-7}$.
A larger tolerance was chosen for rigid SPC/E water due to convergence issues in connection with shifted-force Coulomb interactions (see below).

The maximum number of allowed iterations was 50. A restart scheme was applied when the solver did not converge within the chosen tolerance. In this case the solver (and quadratic iteration) was
restarted from the partially estimated ''solution'' adding $2 \cdot 10^{-7}$ to the tolerance. It should be noted that the stabilized biconjugate
gradient algorithm may get ''trapped'', resulting in a break-down of the CUSP linear solver. If this happens, it is
detected by our program, and the solver and quadratic iteration are restarted, with a
smaller number (10) of maximum allowed iterations for the solver.

\subsection{The asymmetric dumbbell}

The asymmetric dumbbell model\cite{moleculeshidden} consists of a large ($A$) and a small ($B$) Lennard-Jones (LJ) particle, rigidly bonded  with bond distance of $r_{AB} = 0.29/0.4963$ (here and
henceforth units are given in LJ units referring to the $A$ particle such that $\sigma_{AA}$ = 1, $\epsilon_{AA}$ = 1, and $m_{A}$ = 1). The asymmetric dumbbell
model has $\sigma_{BB}=0.3910/0.4963$, $\epsilon_{BB}=0.66944/5.726$, and $m_{B}=15.035/77.106$. The $AB$ interaction between different
molecules is determined by the Lorentz-Berthelot mixing rule\cite{tildesley}.
$n$ = 500 molecules (here and henceforth $n$ denotes the number of molecules and $N$ the number of atoms) were used in the simulations
with a pair-potential cut-off of $r_{c} = 2.5$. The step length $l_{0}$ was fixed in the range 0.125-0.138 depending on the state point.

Simulations were also performed where the rigid bonds were replaced with stiff harmonic springs. The spring constant was $k = 3000$, while all other model parameters remained unchanged.

\subsection{Lewis-Wahnstr{\"o}m OTP}

The Lewis-Wahnstr{\"o}m OTP model\cite{otp2} consists of three identical LJ particles rigidly bonded in an isosceles triangle with sides of $r_{AA} = 1$ and
top angle of $75^\circ$. All parameters (including the masses) are unity for the OTP model. $n$ = 320 molecules were simulated and a pair-potential cut-off of $r_{c} = 2.5$ was used.
The step length was 0.100.

\subsection{SPC/E Water}

The SPC/E water model\cite{spce} is an isosceles triangle with sides $r_{OH} = 1/3.166$ and top angle $109.47^{\circ}$. The $OO$ intermolecular interactions are given by the
LJ pair potential ($\epsilon_{OO}=1$, $\sigma_{OO}=1$, and $m_{O}=15.9994/1.00794$). The three particles are
charged with $q_{O}=-22.0$ and $q_{H}=|q_{O}|/2$. $n$ = 2000 molecules were simulated and a pair-potential cut-off
of $r_{c} = 6.28$ for both LJ and Coulomb interactions was applied\cite{shiftedforceCoulomb,shiftedforceHansen}. The step length was fixed in the range 0.06-0.07
depending on the state point. For this system the numerical
stability is surprisingly sensitive to the cut-off used in the Coulomb interactions, but a larger shifted-force cut-off improves this behavior\cite{shiftedforceHansen}.

\section{Testing the stability of the rigid-bond $NVU$ algorithm}\label{testing}

This section tests the conservation properties of the rigid-bond \textit{NVU} algorithm. Table \ref{ustepconservation} shows the potential energy, the deviation of bond lengths, and step length
as functions of integration step number for Lewis-Wahnstr{\"o}m OTP at $\rho$ = 0.329 and T = 0.700. It is clear that these quantities are conserved by the algorithm
and that no drift occurs. The step length is conserved to the highest accuracy since it is not prone to numerical error in determining the Lagrangian multipliers.

\begin{table}[H]
  \centering
  \begin{tabular}
    {|| p{4.0cm} || p{4.0cm} | p{4.0cm} | p{4.0cm} || }
    \hline
    \hline
    Integration steps & $U$ / $N$ & $(1/G\sum_{\alpha}(r^{\alpha}$ $-$ $C_{\alpha})^{2})^{1/2}$ & $\sum_{k}\tilde{m}_{k}(\boldsymbol{\delta}_{i+1/2}^{(k)})^{2}$ \\
    \hline
    \hline
    $10^{1}$  & -4.42550 & 2.81207 $\cdot 10^{-7}$ & 0.0999999 \\
    $10^{2}$  & -4.42552 & 3.03535 $\cdot 10^{-7}$ & 0.1000000 \\
    $10^{3}$  & -4.42552 & 2.81128 $\cdot 10^{-7}$ & 0.1000000 \\
    $10^{4}$  & -4.42552 & 2.95078 $\cdot 10^{-7}$ & 0.1000000 \\
    $10^{5}$  & -4.42550 & 3.08793 $\cdot 10^{-7}$ & 0.1000000 \\
    $10^{6}$  & -4.42551 & 2.90477 $\cdot 10^{-7}$ & 0.1000000 \\
    \hline
    \hline
  \end{tabular}
  \caption{Potential energy, deviation of bond lengths and step length as functions of integration step number in the \textit{NVU} algorithm for Lewis-Wahnstr{\"o}m OTP ($\rho$ = 0.329, T = 0.700).
    Single-precision floating-point arithmetic was used for the simulations.}
  \label{ustepconservation}
\end{table}
Figure \ref{lambdadist} shows the distribution of the term $l_{0}\lambda_{i}\langle m \rangle$ in Eq. (\ref{some}) (recall $\tilde{m}_{k} = m_{k}/\langle m \rangle$).
In \textit{NVU} dynamics there is, as such, no notation of time;
a geodesic on the manifold can be traversed with any velocity. Comparing the \textit{NVU} algorithm of Eq. (\ref{some})
to the rigid-bond Verlet algorithm\cite{shake} ${\bf r}^{(k)}_{i+1} \,=\, 2{\bf r}^{(k)}_i - {\bf r}^{(k)}_{i-1} + ((\Delta t)^{2}/m_{k})[\textbf{f}^{(k)}_i + \textbf{g}^{(k)}_i]$,
we can define the term $l_{0}\lambda_{i}\langle m \rangle$ as
a varying ''time step'' length of the \textit{NVU} algorithm (see also Paper II), i.e.,

\begin{equation}\label{time}
  (\Delta t_{i,NVU})^{2} \equiv l_{0}\lambda_{i}\langle m \rangle.
\end{equation}
The integration steps of the \textit{NVU} algorithm are thus henceforth referred to as ''time steps''. The average of
Eq. (\ref{time}) is used in Sec. \ref{sampling} when comparing to \textit{NVT} dynamics. As was the case for the atomic \textit{NVU}
algorithm (Paper I), $l_{0}\lambda_{i}\langle m \rangle$ is Gaussian distributed for large systems and its relative variation decreases as
the number of particles increases. It thus becomes a better and better approximation to treat this term as constant, implying equivalent
sampling properties  of \textit{NVU} and \textit{NVE} dynamics also when rigid bonds are included in the simulations.\newline \newline

\begin{figure}[H]
  \centering
  \includegraphics[width=70mm]{lambda_dist_otp}
  \caption{The probability density of the ''time step'' length $(\Delta t_{i,NVU})^{2} \equiv l_{0}\lambda_{i}\langle m \rangle$  of the rigid-bond \textit{NVU} algorithm for
    Lewis-Wahnstr{\"o}m OTP at $\rho$ = 0.329 and T = 0.700. $n = 320$ molecules were used in the simulations.}
  \label{lambdadist}
\end{figure}

\section{Sampling properties of the rigid-bond $NVU$ algorithm}\label{sampling}

The \textit{NVU} algorithm is now compared to \textit{NVT} dynamics for the three different models.
First, we consider the asymmetric dumbbell model\cite{moleculeshidden}, both rigid and flexible. Afterwards, the {Lewis-Wahnstr{\"o}m OTP model\cite{otp2}, and
finally rigid SPC/E water\cite{spce}.

\subsection{The asymmetric dumbbell model}

In Figs. \ref{dumbrigid}(a) and (b) are shown, respectively, the molecular center-of-mass (CM) radial distribution functions and the CM incoherent intermediate scattering functions for the rigid
asymmetric dumbbell model\cite{moleculeshidden} for different temperatures at $\rho = 0.932$. The black circles and curves
give \textit{NVT} simulation results while the red crosses give the \textit{NVU} simulation results. The two radial distribution functions in Fig. \ref{dumbrigid}(a)
agree very well, and this is also the case for the dynamics in Fig. \ref{dumbrigid}(b). \newline \newline

\begin{figure}[H]
  \begin{minipage}[t]{0.45\linewidth}
    \centering
    \includegraphics[width=70mm]{db_rigid_rdf}
  \end{minipage}
  \hspace{0.5cm}
  \begin{minipage}[t]{0.45\linewidth}
    \includegraphics[width=70mm]{db_rigid_fs}
    \centering
  \end{minipage}
  \caption{Comparison of structure and dynamics in \textit{NVU} and \textit{NVT} simulations of the rigid asymmetric dumbbell model.
    The black circles and curves give \textit{NVT}, the red crosses \textit{NVU} simulation results. (a) The molecular CM radial
    distribution functions at $\rho$ = 0.932 and T = 0.500. (b) The molecular CM incoherent intermediate scattering functions
    at $\rho$ = 0.932 and T = 0.500, 0.600, 0.700, 0.800, 0.900. }
  \label{dumbrigid}
\end{figure}
For reference, we also simulated (Fig. \ref{dumbflex}) the corresponding quantities for the flexible-bond asymmetric dumbbell model at the state points of Fig. \ref{dumbrigid}. Again, there
is a very good agreement between \textit{NVU} and \textit{NVT} dynamics. \newline \newline

\begin{figure}[H]
  \begin{minipage}[t]{0.45\linewidth}
    \centering
    \includegraphics[width=70mm]{db_flex_rdf}
  \end{minipage}
  \hspace{0.5cm}
  \begin{minipage}[t]{0.45\linewidth}
    \includegraphics[width=70mm]{db_flex_fs}
    \centering
  \end{minipage}
  \caption{Comparison of structure and dynamics in \textit{NVU} and \textit{NVT} simulations of the flexible-bond asymmetric dumbbell model.
    The black circles and curves give \textit{NVT}, the red crosses \textit{NVU} simulation results. The same state points as in Fig. \ref{dumbrigid} were simulated.
    (a) The molecular CM radial distribution functions at $\rho$ = 0.932 and T = 0.500. (b) The molecular CM incoherent intermediate scattering functions
    at $\rho$ = 0.932 and T = 0.500, 0.600, 0.700, 0.800, 0.900. }
  \label{dumbflex}
\end{figure}

\subsection{Lewis-Wahnstr{\"o}m OTP}

We show in Figs. \ref{otp}(a) and (b), respectively, the molecular CM radial distribution functions and CM incoherent intermediate scattering functions for
the Lewis-Wahnstr{\"o}m OTP model\cite{otp2}. The same symbols and meanings as in the preceding section are used. Again, the \textit{NVU} and \textit{NVT} simulations agree very well for both
structure and dynamics. \newline \newline

\begin{figure}[H]
  \begin{minipage}[t]{0.45\linewidth}
    \centering
    \includegraphics[width=70mm]{otp_rdf_cm}
  \end{minipage}
  \hspace{0.5cm}
  \begin{minipage}[t]{0.45\linewidth}
    \includegraphics[width=70mm]{otp_Fs_cm}
    \centering
  \end{minipage}
  \caption{Comparison of center-of-mass structure and dynamics in \textit{NVU} and \textit{NVT} simulations of the Lewis-Wahnstr{\"o}m OTP model.
    The black circles and curves give \textit{NVT}, the red crosses \textit{NVU} simulation results. (a) The molecular CM radial
    distribution functions at $\rho$ = 0.329 and T = 0.700. (b) The molecular CM incoherent intermediate scattering functions
    at $\rho$ = 0.329 and T = 0.700, 0.800, 0.900, 1.000. }
  \label{otp}
\end{figure}
For comparison, we also show in Fig. \ref{otp1} the corresponding particle quantities for the OTP model. \newline \newline

\begin{figure}[H]
  \begin{minipage}[t]{0.45\linewidth}
    \centering
    \includegraphics[width=70mm]{otp_rdf_particle}
  \end{minipage}
  \hspace{0.5cm}
  \begin{minipage}[t]{0.45\linewidth}
    \includegraphics[width=70mm]{otp_Fs_particle}
    \centering
  \end{minipage}
  \caption{Comparison of particle structure and dynamics in \textit{NVU} and \textit{NVT} simulations for the Lewis-Wahnstr{\"o}m OTP model.
    The black circles and curves give \textit{NVT}, the red crosses \textit{NVU} simulation results. (a) The particle radial
    distribution function at $\rho$ = 0.329 and T = 0.700. (b) The particle incoherent intermediate scattering functions
    at $\rho$ = 0.329 and T = 0.700, 0.800, 0.900, 1.000. }
  \label{otp1}
\end{figure}

\subsection{SPC/E Water}

Finally, we consider in Fig. \ref{spce} the same quantities as above for the rigid SPC/E water model\cite{spce}. Again,
full equivalence between \textit{NVU} and \textit{NVT} dynamics is found. \newline \newline

\begin{figure}[H]
  \begin{minipage}[t]{0.45\linewidth}
    \centering
    \includegraphics[width=70mm]{spce_rdf_cm}
  \end{minipage}
  \hspace{0.5cm}
  \begin{minipage}[t]{0.45\linewidth}
    \includegraphics[width=70mm]{spce_Fs_cm}
    \centering
  \end{minipage}
  \caption{Comparison of structure and dynamics in \textit{NVU} and \textit{NVT} simulations of rigid SPC/E water.
    The black circles and curves give \textit{NVT}, the red crosses \textit{NVU} simulation results. (a) The molecular CM radial
    distribution functions at $\rho$ = 1.000 and T = 3.800. (b) The molecular CM incoherent intermediate scattering functions
    at $\rho$ = 1.000 and T = 3.800, 4.200, 5.000. }
  \label{spce}
\end{figure}
The linear algorithm for determining the Lagrangian multipliers presented in Sec. \ref{alternative} (Eqs. (\ref{linear1}) and (\ref{linear2})) is
tested in Fig. \ref{spcelinear} by probing the molecular CM radial distribution functions. \textit{NVU} and \textit{NVT} dynamics
also here give identical results. \newline \newline

\begin{figure}[H]
  \centering
  \includegraphics[width=70mm]{spce_rdf_cm_linear}
  \caption{Comparison of structure in \textit{NVU} and \textit{NVT} simulations of rigid SPC/E water at $\rho$ = 1.000 and T = 3.800 applying the
    linear method to determine the Lagrangian multipliers (Eqs. (\ref{linear1}) and (\ref{linear2})). The bond lengths are here
conserved to order $10^{-6}$ in the standard deviation of the bonds (using single-precision).}
  \label{spcelinear}
\end{figure}
We conclude from the presented results that for sufficiently large molecular systems with
flexible and/or rigid bonds, \textit{NVU} dynamics is equivalent to Nosé-Hoover \textit{NVT} dynamics (and, by implication, Newtonian \textit{NVE} dynamics).

\section{Summary}\label{sum}

\textit{NVU} dynamics is molecular dynamics at constant potential energy realized by tracing out a geodesic on the constant-potential-energy
hypersurface $\Omega$ (Eq. (\ref{omega_def})). In Papers I and II\cite{nvu1,nvu2}, a ''basic'' and a ''stabilized'' atomic \textit{NVU} algorithm for simulating geodesics on $\Omega$ was developed. The
basic \textit{NVU} algorithm has excellent stability; it is time-reversible and symplectic, and the stabilized algorithm was developed only to prevent accumulation of
numerical error. It was found that atomic \textit{NVU} dynamics becomes equivalent to atomic \textit{NVE} dynamics in the thermodynamic limit.

In this paper the  stabilized \textit{NVU} algorithm has been extended to simulate molecules at constant potential energy. Molecules are simulated by introducing rigid and/or flexible bonds in the models.
The atomic \textit{NVU} algorithm keeps the potential energy constant and can thus right away simulate flexible bonds.
The focus here was on incorporating rigid bonds in the framework of \textit{NVU} dynamics, which leads to the introduction of additional Lagrangian multipliers beyond
those of the constraint of constant potential energy. This is completely analogous to the approach for simulating rigid bonds in
standard Newtonian \textit{NVE} dynamics\cite{shake,toxconstraintnve,toxconstraintnph}. In the \textit{NVU} algorithm, a set of coupled quadratic equations
was constructed for calculating the Lagrangian multipliers and solved in an iterative manner as a linear system, a procedure developed for rigid-bond \textit{NVE} dynamics
in the \textit{MILC-SHAKE} algorithm\cite{milcshake}. In addition, a set of linear equations was presented for calculating the Lagrangian multipliers, and appears to be a
promising new way of simulating rigid bonds.

The rigid-bond \textit{NVU} algorithm reduces to the atomic
\textit{NVU} algorithm when there are no rigid bonds. The algorithm was tested on three different model systems: the asymmetric dumbbell model, Lewis-Wahnstr{\"o}m OTP, and rigid SPC/E water.
The probed quantities in the simulation gave
identical results to those of Nos$\acute{e}$-Hoover \textit{NVT} dynamics.
We conclude that also for molecular systems do \textit{NVU} dynamics become equivalent to \textit{NVE} dynamics in the thermodynamic limit (since \textit{NVE} and \textit{NVT}
dynamics are known to give equivalent results\cite{evans1985}).

\acknowledgments

The centre for viscous liquid dynamics ``Glass and Time'' is sponsored by the Danish National Research Foundation (DNRF). The authors are grateful to Ole J. Heilmann for pointing
out the alternative method for determining the Lagrangian multipliers (Sec. \ref{alternative}).


\section{Introduction}
    Fusion devices contain fast particle populations due to external plasma
    heating and (eventually) fusion-borne $\alpha$-particles.
    Fast particle populations can interact with global electromagnetic waves,
    leading to the growth of MHD-like and kinetic instabilities -- e.g., Toroidicity
    Induced Eigenmodes (TAE) \cite{Cheng85, Cheng86}, Reversed Shear
    Alfv{\'e}n Eigenmodes (RSAE) \cite{Berk01}
    or Beta Induced Alfv{\'e}n Eigenmodes (BAE) \cite{Heidbrink93, Turnbull93}.\\\\
    In this work, drift-kinetic fast particle simulations performed
    with the \textsc{Hagis} code (\cite{Pinches98,sip_phd}, shortly introduced in \sref{hagis}) are
    carried out to obtain a deeper understanding of
    the processes in phase space that occur due to fast particle interaction
    with MHD modes. Of interest is the dynamics of wave-fast particle interaction in
    scenarios with two modes of different frequencies: what are the mode coupling mechanisms, and
    what is the dependence on radial mode distance?
    \Sref{theory} gives an overview of the theoretical basis for
    these coupling mechanisms.
    \Sref{hagis} explains in brief, which ASDEX Upgrade reference scenario is selected and
    how the MHD equilibrium is processed to give the basis for the \textsc{Hagis} simulations.
    In \sref{numstudy} numerical studies are described, investigating the stochastic threshold and the influence of $q$
    on the mode drive. It is found that the presence of a second mode can enhance growth rates as
    well as mode amplitudes significantly. As this effect is based on particles that have resonance regions in phase
    space with both modes, it is referred to as ``double-resonance''. The double-resonant effect was expected
    to decrease with the radial mode distance, but it turned out, that the picture is more complicated.

\section{Theoretical Picture of Double-Resonance}\label{theory}
    Theory (e.g.\ Ref.\ \cite{Berk92}) predicts that conversion of free energy
    to wave energy is enhanced in a multiple-mode scenario, i.e., the interaction of multiple modes produces
    energy conversion rates higher than that which would be achieved with each mode acting independently.
    This can be partially explained by the principle of gradient (of the radial particle distribution) driven
    mode growth -- according to\footnote{Throughout this work, $s$ refers to the radial coordinate as
    the square root of the normalized poloidal flux:
    $s=(\psi_{\mathrm{pol}}/\psi_{\mathrm{pol, edge}})^{1/2}$.}
    $\gamma \propto \nabla f(s)$ \cite{Fu89} -- which can be extended to
    multiple modes \cite{Berk90-I, Berk90-II, Berk90-III,Berk92}.
    This picture of \emph{\textbf{gradient driven double-resonance}} is based on the
    precondition that modes share resonances in the same phase space area.
    Through the resulting redistribution by each mode, a steeper gradient
    is produced at the other mode's position, enhancing its drive. The
    overlapping of modes leads then to a much larger conversion of free energy to wave
    energy.\\
    However, this mechanism can only work if there is also \emph{spatial} mode overlap in
    the radial direction. In Ref.\ \cite{mwb_phd} simulations were carried
    out, finding a double-resonant effect also without this precondition.
    Furthermore, a superimposed oscillation on the modes' amplitudes was observed,
    clearly indicating mode-mode interaction. The modes without radial overlap
    are then coupled radially through the particles' trajectories: A population
    of particles that shares resonances in phase space with both modes and passes
    both modes' location at once, can transfer energy from one mode to the other \cite{mwb_phd}.
    Since the particle orbits are characterized by a certain width, it is not necessary that
    both modes have a radial overlap.
    In the following, this mechanism is called \emph{\textbf{inter-mode energy transfer}}:
    By damping one mode, particles gain energy $E$ and also toroidal momentum $P_{\zeta}$
    due to \cite{sip_phd}
    \begin{equation}\label{EPzeta}
        (E - \frac{\omega}{n}P_{\zeta}) = \mathrm{const.}
    \end{equation}
    As Alfv{\'e}nic mode frequencies are very low (compared to the particles' cyclotron frequency),
    momentum transfer dominates over  energy transfer.
    Since $P_{\zeta} \propto -\Psi$, particles gaining energy from the wave (i.e.\ damping it),
    are redistributed radially inwards, whereas particles losing energy to the
    wave (i.e.\ driving it), are redistributed radially outwards. The latter is the
    dominant process and is caused by the negative
    slope of the particles' radial distribution function.
    When passing through the second mode, the particles lose energy and toroidal momentum by
    driving the mode.
    As there is no radial net drift in this mechanism, it can continue as long as the dominant
    mode is strong, making this the dominant process over other possible
    combinations of mode-mode energy exchange. Particles that gain energy from
    both modes or lose energy to both modes soon leave the resonant phase space
    area.
    It is the exchange of energy between the modes that leads to the observed
    oscillation of their amplitudes: The mode receives energy in the rhythm of the
    particles' bounce frequency $\omega_{\mathrm{b}}$,
    which equals the beat frequency of both modes $\Delta \omega$ as shown in the following:
    Trapped particles with a bounce frequency $\omega_{\mathrm{b}}$ and toroidal
    precession frequency $\omega_{\mathrm{tp}}$
    interact with MHD modes of a certain frequency $\omega$,
    if the \emph{resonance condition} \cite{porcelli94}
    \begin{equation}\label{rescond}
        \omega - n\omega_{\mathrm{tp}}-p\omega_{\mathrm{b}} \approx 0
    \end{equation}
    is fulfilled, where $n$ is the toroidal mode
    harmonic and $p$ the particles' bounce harmonic.
    For double-resonance, this resonance condition has to be fulfilled for both modes '1' and '2'
    simultaneously, leading to
    \begin{equation}\label{doublerescond}
      \omega_{\mathrm{tp}}(n_{1}-n_{2})+\omega_{\mathrm{b}}(p_{1} - p_{2})
      = \omega_{\mathrm{1}}-\omega_{\mathrm{2}} \equiv \Delta \omega
    \end{equation}
    For the simplest case, the toroidal mode numbers $n$ are considered equal,
    and $p_{1} = 1$, $p_{2} = 0$.
    This leads to $\omega_{\mathrm{b}} = \Delta \omega$.
    Indeed, the lowest bounce harmonics $p=0, \pm 1$ are the most relevant ones;
    however, the issue of different toroidal mode numbers $n$ will be discussed later.

\section{Using the Hagis Code with an ASDEX Upgrade Plasma Equilibrium}\label{hagis}
    The numerical investigations presented in this work are performed with
    the \textsc{Hagis} Code \cite{Pinches98, sip_phd},
    a nonlinear, drift-kinetic, perturbative Particle-in-Cell code (current release 12.05).
    \textsc{Hagis} models the interaction between a distribution of energetic particles and a set
    of Alfv{\'e}n Eigenmodes. It calculates the linear growth rates as well as the nonlinear behavior
    of the mode amplitudes and the fast ion distribution function that are determined by
    kinetic wave-particle nonlinearities. It is fully
    updated to work with MHD equilibria given by the recent \textsc{Helena} \cite{Huysmans91} version.\\
    \begin{minipage}{1\textwidth}
      \begin{minipage}{0.45\textwidth}
        \begin{figure}[H]
           \hspace{1cm}\includegraphics[width=0.7\textwidth]{./figure1.eps}
           \caption{\itshape Hagis equilibrium with vacuum region (red),
             representative banana orbit (blue) and FILD position (black square).}
           \label{hagis_equilibrium}
        \end{figure}
        \vfill
      \end{minipage}
      \hfill
      \begin{minipage}{0.45\textwidth}
        \begin{figure}[H]
          \includegraphics[width=0.9\textwidth, height=4.7cm]{./figure2.eps}
          \caption{\itshape $q$ profile\ according to experimental measurements in
            AUG discharge \#23824 at different time
          points: $t=1.16$\ s (black solid line)  and  $t=1.51$\ s (blue dashed line):
          the $q$ profile s differ in the shape (the earlier one is inverted, the later one monotonic), but
          also in the absolute values (the inverted $q$ profile\ has higher absolute $q$ values).}
          \label{q23824}
        \end{figure}
      \end{minipage}
    \end{minipage}
    The plasma equilibrium (\fref{hagis_equilibrium}) for \textsc{Hagis} is based
    on the \textsc{Cliste} code \cite{Carthy99, Carthy12},
    then transformed via \textsc{Helena} to straight field line coordinates and
    to Boozer coordinates by \textsc{Hagis}.\\
    The data for the MHD equilibria of all simulations presented in the following originate from the
    ICRH minority heated ASDEX Upgrade (AUG) discharge \#23824, at time  $t=1.16$\ s
    or $t=1.51$\ s.
    At the earlier time point, the $q$ profile\ is slightly inverted (\fref{q23824}, black solid line),
    whereas at the later time point, it is monotonic, with lower absolute values (\fref{q23824},
    blue dashed line).\\
    The plasma equilibrium, and in particular the $q$ profile\ is determined by Alfv{\'e}n spectroscopy of
    the RSAEs:
    magnetic pick-up coil data and soft X-ray emission measurements are used to determine the
    $q$ profile's minimum value and location.\\
    This particular discharge was chosen due to the availability of especially detailed
    experimental data concerning fast particle-mode interaction. However, a detailed comparison between
    numerical and experimental results goes beyond the scope of this work, but is the
    object of current and future
    investigation. In this work, a numerical study is presented, that is based on an AUG reference case:
    the $q$ profile\ as well as the mode frequencies, harmonics and widths are chosen according to
    experimental observation. However, significant abstractions have also been made:
    the mode structures are chosen analytically, as well
    as the particle distribution functions. Furthermore, mode saturation is based completely on
    the depletion of the particle distribution function, no collisions, nor turbulence are taken
    into account, and there is no source term for energetic particles.

\section{Numerical Study on Double Mode Resonance}\label{numstudy}
  In this section, a numerical study is presented that examines
  under which circumstances mode-particle
  interaction can take place. A special focus is set on the interaction
  with multiple modes of different frequencies.

  \subsection{Simulation Conditions}
    As the question of understanding double mode resonance is very fundamental,
    the simulations were performed under quite simple, but still realistic physical conditions:
    the volume averaged fast particle beta is chosen to be $\beta_{\mathrm{f}}$\ $= 1\%$.
    To avoid different mode drive at different radial mode positions
    only due to a steeper gradient in the distribution function, a radial particle distribution $f(\psi)$
    with constant gradient is chosen.
    As energy distribution function $f(E)$ a slowing down function \cite{Gaffey76} is used:
    \begin{equation}\label{slowdown}
      f(E) = \frac{1}{E^{3/2}+E^{3/2}_{\mathrm{c}}}\mathrm{erfc}\left(\frac{E-E_0}{\Delta E}\right)
    \end{equation}
    ($\Delta E=0.1499$\ MeV, $E_{\mathrm{c}}=0.1934$\ MeV, $E_0=1.0$\ MeV).
    The particles are distributed isotropically in pitch angle (as
    e.g., fusion-borne $\alpha$ particles would be).
    As MHD perturbations, analytic, gauss shaped functions are used (see \fref{pert_exam}),
    without background damping. The mode frequencies are chosen to
    match experimental data: at both time points within the considered discharge \#23824,
    a high frequency TAE with $120$\ kHz is found, as well as
    a lower frequency mode at $55$\ kHz. Ref.\ \cite{Lauber09} describes
    their mode as a RSAE at $t= 1.16$\ s
    and as BAE at $t= 1.51$\ s. The widths of the modes are based on the MHD eigenfunctions, calculated
    numerically with the linear eigenvalue solver \textsc{Ligka} \cite{Lauber07}.
     \begin{figure}[H]
        \centering
        \includegraphics[width=0.4\textwidth, height=3.5cm]{./figure3.eps}
        \caption{\itshape Analytical, gauss-shaped perturbation as used in the simulations.}
        \label{pert_exam}
     \end{figure}
    Convergence tests indicate that simulations with 120000 markers are
    sufficient, leading to computational costs of a few hundred CPUh for each simulation.

  \subsection{Simulation Results}
    Multi-mode simulations are carried out with different radial distances $\Delta s$ between
    the Alfv{\'e}nic modes.
    It is helpful to look at the different stages within the simulation individually:
    although all simulations are performed including all nonlinearities, the nonlinear
    effects are small at the beginning of each simulation, leading to an amplitude
    growing according to $\propto \mathrm{exp}(\gamma t)$ with $\gamma$ the so-called
    linear growth rate.
    \paragraph{The linear regime}
    The modes' amplitudes in this linear regime, as depicted in
    \fref{run0161+0162+0163+0165+0166+0175_linamp}
    confirm already some of the results presented in \cite{mwb_phd}:
    at least one mode grows faster in the double mode scenario compared
    to the single mode simulation (the TAEs in \fref{run0161+0162+0163+0165+0166+0175_linamp}
    dark colors).
    In most cases, both modes have larger growth rates in the double mode scenario than
    in the single mode simulation. Furthermore, it is usually the radially outer mode that
    is driven most strongly due to the double mode resonance.
    This is because \textit{gradient driven double-resonance} is most effective in
    driving the outer mode: through resonant particle redistribution, one mode produces steeper
    gradients at the other mode's position. Through the inner mode, particles are
    redistributed towards the outer position and supply it with more (resonant) particles,
    that can transfer their energy to the wave. At the inner mode's position, in contrast,
    the particle population decreases, and with it, the energy transfer to the wave.\\
    However, in simulations with low growth rates
    (see linear phase of \fref{run0186+0187+0188_amplitude}),
    the stronger -- and that is mostly the outer mode -- was enhanced
    in mode drive less effectively than the subdominant one, or even weakened
    compared to the single mode case. In these cases, the \textit{inter-mode energy exchange}
    is the prevailing mechanism, conducting energy from the stronger
    to the weaker mode.\\
    A superimposed oscillation is visible, strongest if the
    modes' growth rates are relatively
    different from each other. However, if they differ too much,
    the oscillation becomes jagged, due to
    the large impact of the dominant mode on the subdominant one.
    The oscillation frequency matches the beat frequency of both
    modes $\Delta \omega = \omega(\mathrm{TAE}) - \omega(\mathrm{RSAE}) = 65~$kHz,
    a hint for working inter-mode energy transfer.\\
    As stated in Ref.\ \cite{mwb_phd}, both mode numbers $n$ have to be equal to lead
    to the oscillatory behavior, i.e.\ to double mode resonance, as
    will be discussed later.\\\\
    In the course of each simulation, nonlinear effects become more and more important,
    leading to a mode amplitude evolution different from $\propto \mathrm{exp}(\gamma t)$.
    The beginning of the nonlinear regime is observed, when the mode amplitude changes
    from evolving linearly (in a semilogarithmic plot) to saturating at a certain level.\\
    \begin{figure}[H]
      \centering
      \includegraphics[width=0.6\textwidth]{./figure4.eps}
      \caption{\itshape Mode amplitude evolution in linear growth phase for different
        multi-mode (solid lines) scenarios (for the MHD equilibrium at $t=1.16$\ s),
        differing in radial mode distances $\Delta s$. For comparison, the single mode
        simulations are depicted as dashed lines. In all cases, the TAE
        (120\ kHz at $s=0.4$, dark colors)
        grows stronger in the double mode scenario, the RSAE
        (55\ kHz at various radial positions, reddish colors) in contrast,
        benefits from the double mode case better, if $\Delta s$ is small enough.}
      \label{run0161+0162+0163+0165+0166+0175_linamp}
    \end{figure}

    \paragraph{The nonlinear regime}
    The effect of the double-resonance in the nonlinear regime is clearly visible
    in \fref{run0163+0166+0167_amplitude} ($\Delta s =0.2$):
    \begin{figure}[H]
      \centering
      \includegraphics[width=0.6\textwidth, height=5.5cm]{./figure5.eps}
      \caption{\itshape Mode amplitudes evolution in double mode (solid lines)
        versus single mode (dashed) simulation,
        with a RSAE (pink) at $s=0.2$ and a TAE (blue) at $s=0.4$ for the MHD
        equilibrium at $t=1.16$\ s. In the double mode case, one
        can clearly see a superimposed oscillation and higher saturation levels.}
      \label{run0163+0166+0167_amplitude}
    \end{figure}
    in the double mode scenario, the modes not only grow faster (in the case of the TAE),
    but also saturate at a higher level compared to the single mode case (by about a factor of five to ten).
    This is due to the orbit \textit{stochasticization} at
    $\delta B/B \approx 2 \cdot 10^{-3}$
    already discovered in single mode simulations: The stochastic
    threshold \cite{Berk92} is only reached in the double mode scenario, not in the single mode case.
    To visualize the effect of stochasticization, an representative particle orbit and the particle's energy
    evolution \textit{before} and \textit{after} the stochasticization sets in is shown
    in \fref{run0175_orbit}.
    \begin{figure}[H]
      \centering
      \subfigure[\itshape during resonant phase]{\includegraphics[width=0.42\textwidth, height=3.0cm]{./figure6a.eps}}
      \hspace{0.8cm}\subfigure[\itshape during stochastic phase]{\includegraphics[width=0.42\textwidth, height=3.0cm]{./figure6b.eps}}
      \caption{\itshape Representative particle orbit (left, blue line) and energy evolution (right)
        before (a) and after (b) the stochastic threshold is
        reached: During the early saturation phase, the orbit is banana-shaped and shifted outwards, with
        the particle losing its energy. In the stochastic regime, the orbit becomes distorted
        and the particle's energy loss occurs in bursts. The black line in the left plots indicates radial
        envelope of the MHD mode.}
      \label{run0175_orbit}
    \end{figure}
    The picture is not as simple in the second example considered
    ($\Delta s=0.1$, \fref{run0162+0165+0167_ampl}):
    at the beginning of the nonlinear regime,
    double-resonance is effective -- mode amplitudes are higher by a factor of five compared
    to the single mode simulation. At a later time point in the simulation
    however, the situation changes
    completely: the TAE amplitude vanishes and the RSAE saturates at a lower level compared to
    the amplitude it reached in the single mode scenario during stochasticization. This scenario
    was found frequently for close radial mode distances, and shows that a linearly dominant
    mode can be stabilized nonlinearly in a multi-mode scenario.\\
    \begin{figure}[H]
      \centering
      \includegraphics[width=0.6\textwidth, height=5.5cm]{./figure7.eps}
      \caption{\itshape Mode amplitudes evolution in double mode (straight lines)
        versus single mode (dashed) simulation,
        with a RSAE (red) at $s=0.3$ and a TAE (green) at $s=0.4$ for
        the MHD equilibrium at $t=1.16$\ s. Here, the TAE is destroyed
        in the double mode case.}
      \label{run0162+0165+0167_ampl}
    \end{figure}
    The reason for this behavior can be found in the radial particle redistribution.
    To find out the redistribution caused solely by the RSAE, the single mode
    simulation of that mode has to be considered: the RSAE leads to a radial particle redistribution
    as depicted in \fref{run0165_disb}a). This redistribution leads
    to a flattened gradient (\fref{run0165_disb}b) at $s=0.4$,
    the position of the TAE in the
    the double mode case of \fref{run0162+0165+0167_ampl}. The reason for the
    destruction of the TAE in this simulation is therefore its radial location too close to the RSAE.
    However, at all positions $s > 0.45$, the gradient becomes steeper in the
    course of the single mode simulation. Indeed, simulating the RSAE with a TAE in this
    radial range leads to strong TAE drive.
    \begin{figure}[H]
      \centering
      \includegraphics[width=0.4\textwidth]{./figure8a.eps}\\
      \includegraphics[width=0.4\textwidth]{./figure8b.eps}
      \caption{\itshape Temporal evolution (early: solid, late: dashed) of the distribution
        function $f(s)$ in the single mode simulation
       ($t=1.16$\ s equilibrium) with a RSAE at $s=0.3$ (upper).
       Lower: the gradient $\mathrm{d}f(s)/\mathrm{d}s$ at $s=0.4$
       shows a transition from negative values (straight) to positive values (dashed),
       leading to a very low amplitude of the
       TAE mode at that position in the double mode simulation.}
      \label{run0165_disb}
    \end{figure}
    A scan over the radial mode distance reveals an effective double-resonance as shown in
    \fref{deltas-scan}: depicted are the ratios of the linear growth rates (a) and the amplitudes
    (b) in the double mode case vs.\ the single
    mode case over the radial mode distance $\Delta s$.
    The amplitude level was compared after $\approx 300$ TAE periods ($=~2.5$\ ms) of
    simulated time. This time is sufficient for the single amplitudes to saturate,
    but still significantly below energy slowing down time. Therefore, the fact that
    no source term is implemented in the code does not disturb the physical picture.
    One can see that the growth rates of
    both the TAE and the RSAE are enhanced in all double mode cases compared to the single mode ones.
    However, the growth rate of the outer TAE is enhanced most strongly and independently
    of the radial
    mode distance -- i.e.\ gradient driven double-resonance works even if there is no radial mode overlap.
    In contrast, the enhancement of the inner and weaker RSAE decreases with the radial
    mode distance for small $\Delta s$. Then it increases again for $\Delta s > 0.15$.
    These larger mode distances match the double-resonant
    particle orbits and therefore enable inter-mode energy exchange, driving the weak mode. For
    higher $\Delta s$, the larger, i.e.\ higher energetic orbits fit the mode distance
    and lead to even more energy exchange.
    Furthermore, with larger radial mode distances, the modes are able to tap energy from a wider
    gradient region. The amplitude ratios, however, are even lowered in the double-resonant case
    compared to the respective single mode levels, if the radial mode distance is small. This happens
    due to the mutual gradient depletion at the
    other mode's radial position. If the modes feature a larger radial distance ($\Delta s > 0.15$), the
    double mode scenario amplitudes are much higher compared to the single mode amplitudes, both for the
    TAE and the RSAE. Both modes benefit from each other -- most for a radial distance of about
    $\Delta s \approx 0.25$.
    It is important to note that the distance $\Delta s=0.25$ giving maximum amplitude ratios depends
    strongly on the absolute mode positions with respect to the radial distribution function,
    and especially on the amplitude regime (stochastic or non-stochastic) of each mode.
    The same applies for the value $\Delta s=0.15$ at which the transition towards
    double-resonant amplitude enhancement takes place.\\
    To summarize, one can say, that growth rates are
    generally enhanced by the presence of another mode, whereas for the amplitudes, this is only
    the case, if the modes feature a sufficient radial distance. For small distances, modes at radial
    positions, where the initial distribution function is already relatively flat, double-resonance
    leads to strong mode stabilization. If the amplitudes are enhanced, their amplification level is,
    however, mainly determined by whether the mode reaches the stochastic regime.
    The probability of reaching the stochastic threshold rises strongly if a second mode is present.
    \begin{figure}[H]
      \centering
      \subfigure[\itshape Growth rates enhancement]{\includegraphics[width=0.37\textwidth,height=4cm]{./figure9a.eps}}
      \hspace{1.0cm}\subfigure[\itshape Amplitudes enhancement (at $\approx 300$ TAE periods)]{\includegraphics[width=0.37\textwidth,height=4cm]{./figure9b.eps}}
      \caption{\itshape Mode evolution scanned over the radial mode distance $\Delta s$ in the
        inverted $q$ profile\ case. Depicted is the $\Delta s$ dependence of the ratios
        of growth rates (a) and amplitudes (b) in double mode
        simulations over those from single mode simulations. Red: RSAE, blue: TAE.
        One can clearly see that larger radial mode distances lead to higher amplitudes,
        whereas amplitudes are even lower than in the single mode case for $\Delta s \le 0.15$.
        The linear growth rates, however, are higher than in the single mode simulation throughout the
        $\Delta s$ range. The RSAE growth rate experiences a small drop at $\Delta s \approx 0.15$.}
      \label{deltas-scan}
    \end{figure}
    However, if the growth rates are relatively low for some reason
    (e.g.\ small mode width, small fast particle beta), the particle redistribution
    is not strong enough to lead to a dominant gradient driven double-resonance.
    In these cases, the \textit{inter-mode energy transfer} mechanism can prevail
    (even in a later phase). The dominant mode is weakened through inter-mode energy transfer
    to the subdominant mode,
    as depicted in \fref{run0186+0187+0188_amplitude}. The process saturates, as
    the modes' amplitudes converge towards comparable levels.
    In this simulation, there is no strong fast particle redistribution.
    \begin{figure}[H]
      \centering
      \includegraphics[width=0.6\textwidth,height=5.5cm]{./figure10.eps}
      \caption{\itshape Mode amplitudes over time in double mode (solid lines) versus single
        mode (dashed) simulation with low growth rates (for the MHD equilibrium
        at $t=1.51$\ s,  $\beta_{\mathrm{f}}$\ $= 2.5\%$ with
        BAE (pink) at $s=0.3$ and TAE (blue) at $s=0.4$). Here, the
        redistribution is too weak to lead to significant
        gradient driven double-resonance. Inter-mode energy
        transfer (from TAE to BAE) dominates.}
      \label{run0186+0187+0188_amplitude}
    \end{figure}
    In the monotonic $q$ profile\ equilibrium (shown in \fref{q23824}),
    the Alfv{\'e}nic modes as simulated within this study do not grow
    up to the stochastic regime. In order to also be able to investigate
    the stochastic phase in the monotonic profile case, one can shift the modes further outward radially,
    where the mode drive is higher. Such simulations revealed interesting possibilities for interplay
    between the two double-resonance mechanisms: e.g., it is possible that the mode in the
    flattened-gradient region is not weakened, but rather the \emph{other} mode.
    In this scenario the mode at the position
    of the flattened gradient is driven solely through energy transfer
    from the other mode, as there is no gradient-caused mode drive remaining.
    Eventually, this leads to a depletion of the other (stronger) mode's amplitude.

    \subsubsection{The role of trapped particles}
    Another question regards the role of trapped particles in the double-resonance mechanisms:
    Especially for the inter-mode energy exchange, the particles' orbit
    widths are crucial. Therefore, the
    question arises if the major part of the double-resonance is carried
    by the trapped particle population or by the passing particles.\\
    In \textsc{Hagis} it is possible not to load any trapped particles at all. In
    the simulations presented in the following,
    the code was modified to even neglect particles that become trapped
    in the course of the simulation.\\
    It was found that for the linear phase, the absence of trapped particles
    leads to significantly lower growth rates. This effect is slightly stronger
    for inner mode positions, although there would be more trapped particles further
    outside. This is due to the fact that the orbit width broadens with radial position,
    and therefore also the loss region. But especially the resonance with a broader TAE
    is effective only with broader orbit particles that are more easily lost.\\
    However, except for the effect of generally lower amplitudes in the scenario without trapped particles,
    the double-resonance also works with passing particles alone:
    the effect of double-resonance
    remains, which is reflected by the ratio of growth rates and amplitude levels in double mode simulations
    over those from single mode simulations. This means, passing as well as banana orbits
    can be adequate orbits for the double-resonance mechanism of inter-mode energy transfer:
    As for passing particles, the orbit does not close entirely on the low field side,
    the double-resonance mechanism must work with the higher poloidal harmonics $m$ of one mode,
    located at the high field side. For the TAE, due to strong ballooning effects, there are no
    significant maxima at the high field side. Core-localized modes, such as
    BAE and RSAE do not couple to the next harmonic $m+1$ and are cylindrical modes that
    feature higher harmonics at the high field side. There, double-resonant passing particles
    can transit the mode. Compared to the banana width $w$ of a trapped particle, the
    drift displacement $\Delta r$ of a passing particle's orbit (at the same energy) is
    significantly smaller. But, as the resonance condition is different for passing particles,
    higher energetic particles are responsible for double-resonance. These feature a larger
    drift displacement, that fits again the simulated radial mode distances $\Delta s$.\\
    However, for the gradient driven (double as well as single) resonance, the trapped particles'
    resonances are very important. It was not possible to quantify this further.
    Without trapped particles, the saturation levels are much lower and do not reach the stochastic regime.
    This is due to the importance of the trapped particles' (single) resonances, leading to gradient driven
    mode-particle interaction. The importance of the trapped particles is consistent
    with the fact, that the TAE is an even TAE, which is known to interact very strongly with trapped particles.

    \subsubsection{The importance of equal toroidal mode numbers}
    Especially for the inter-mode energy transfer, there has to be a population of particles that
    is resonant with both modes at once. For two modes with different toroidal mode numbers $n$,
    the overlapping volume of resonances in phase space becomes smaller, and it becomes
    less probable to find particles that fulfill the resonance condition eq.\ (\ref{doublerescond}) for
    two modes at once, and it is expected that less double mode resonance effects occur,
    if the modes differ in their toroidal mode number $n$.
    \Fref{run0162+0165+0167+0227+0228_amplitude} shows the resulting amplitude
    evolution in a double mode simulation
    with different $n$ for the RSAE ($n=4$) and TAE ($n=5$).
    In the linear phase, one can see clearly that the amplitudes of the single TAE
    simulations are almost exactly equal, independent on $n$. However, in the double mode scenario,
    the RSAE is enhanced (compared to the single mode amplitude), and the TAE
    even more strongly -- but only
    if both modes have equal $n$. If this is not the case, only the TAE is enhanced, whereas the
    RSAE actually decreases. Furthermore, the superimposed oscillation
    vanishes to a slightly jagged curve. This shows that gradient driven
    double-resonance still works with different toroidal $n$ -- the outer mode amplitude is still
    enhanced -- but the inter-mode energy transfer mechanism breaks down, as the resonance
    condition cannot be easily fulfilled for both modes at once for particles with trajectories
    through both modes: the inner RSAE amplitude is not enhanced, but weakened.
    When reaching saturation, the effect of inter-mode energy transfer is overlapped by
    the dominance of the gradient driven double-resonance, and therefore, the saturation levels do
    not depend strongly on equal $n$. Simulations with different $n$ show almost identical nonlinear behavior,
    although saturation is reached later. The stochastic threshold does not seem to be influenced by
    the differing $n$, nor the effect that one mode can be damped by the
    redistribution of the other one. However, the differences again depend on the exact
    scenario.
    \begin{figure}[H]
      \centering
      \includegraphics[width=0.6\textwidth,height=5.5cm]{./figure11.eps}
      \caption{\itshape Mode amplitudes evolution in double mode simulations with equal
        toroidal mode numbers $n = 4$ (pink solid line: RSAE, blue solid line: TAE)
        and with different $n$: RSAE (red solid, $n=4$) and TAE (green solid $n=5$).
        For comparison, the single mode simulations are given (dashed). The RSAE is
        located at $s=0.3$, the TAE at $s=0.4$.}
          \label{run0162+0165+0167+0227+0228_amplitude}
    \end{figure}

\section{Conclusions}
    The interaction of fast particles with Alfv{\'e}n Eigenmodes of different frequencies was
    studied numerically with the \textsc{Hagis} Code
    in a simple, but physically realistic picture for two different MHD equilibria occurring
    during the ASDEX Upgrade discharge \#23824.
    Double-resonant mode drive was compared to single mode scenarios, verifying previous
    findings of double-resonance mechanisms:
    gradient driven double-resonance \cite{Berk92} and inter-mode energy transfer \cite{mwb_phd}.
    The latter was observed to prevail only in low amplitude cases,
    enhancing the weaker mode at the expense of the dominant one.
    A superimposed beat frequency oscillation on both modes' amplitudes was found,
    as well as higher linear growth rates -- at least for one mode -- compared to the single mode
    reference cases.
    The growth rate enhancement was observed to have only a very weak dependence
    on the radial mode distance.
    Concerning the amplitudes of the Alfv{\'e}n Eigenmodes,
    double-resonance can enhance modes into the stochastic regime and therefore lead to much higher
    saturation levels compared to the single mode scenarios, even if there
    is no radial mode overlap.
    However, close radial mode distances were found to destroy
    linearly dominant modes when the nonlinear regime sets in.
    As a consequence, linearly weaker modes may become nonlinearly dominant.
    These results reveal a complex nonlinear evolution of multi-mode scenarios, rather
    than merely a simple continuation of the linear multi-mode behavior.\\\\
    The authors wish to thank D.R.\ Hatch for proof reading.
\bibliographystyle{iopart-num}