\documentclass{article}
\usepackage{fullpage}
\parindent=0pt
\parskip=5pt
\usepackage[stable]{footmisc}
\usepackage{amsmath}
\usepackage{graphicx}
\usepackage{algorithmic}

\newcommand{\bd}{\begin{equation}}
\newcommand{\ed}{\end{equation}}
\newcommand{\bea}{\begin{eqnarray}}
\newcommand{\eea}{\end{eqnarray}}

% The mchem package at CTAN looks pretty useful if one
% were going to write a lot of chemistry...
\newcommand{\AHN}{\mathrm{AH}_n}
\newcommand{\AHthree}{\mathrm{AH}_3}
\newcommand{\AHtwo}{\mathrm{AH}_2}
\newcommand{\AHone}{\mathrm{AH}_1}
\newcommand{\HtO}{\mathrm{H}_2\mathrm{O}}

\newcommand{\sa}{\sigma^\alpha}
\newcommand{\tsa}{{\tilde\sigma}^\alpha}
\newcommand{\sbeta}{\sigma^\beta}
\newcommand{\tsbeta}{{\tilde \sigma}^\beta}
\newcommand{\la}{\lambda_\alpha}
\newcommand{\lb}{\lambda_\beta}
\newcommand{\disa}{\partial_i \sa}
\newcommand{\djsa}{\partial_j \sa}
\newcommand{\djsb}{\partial_j \sbeta}
\newcommand{\minvij}{{R}^{ij}}
\newcommand{\Zab}{Z^{\alpha\beta}}
\newcommand{\tZab}{{\tilde Z}^{\alpha\beta}}
\newcommand{\tZ}{\tilde Z}
\newcommand{\Zs}{Z_s}

\newcommand{\qi}{q^i}
\newcommand{\tqi}{{\tilde q}^i}
\newcommand{\qj}{q^j}
\newcommand{\rprimea}{r'^{\alpha}}
\newcommand{\trprimea}{{\tilde r}'^{\alpha}}
\newcommand{\qprimei}{q'^i}
\newcommand{\tqprimei}{{\tilde q}'^i}

\newcommand{\Cai}{S^{\alpha}_i}
\newcommand{\Caj}{S^{\alpha}_j}
\newcommand{\Cbj}{S^{\beta}_j}
\newcommand{\da}{d_{\alpha}}
\newcommand{\half}{\frac{1}{2}}
\newcommand{\halftsq}{\frac{t^2}{2}}
\newcommand{\halftsqcsq}{\frac{t^2c^2}{2}}
\newcommand{\ra}{r^\alpha}
\newcommand{\rb}{r^\beta}
\newcommand{\Fi}{F_i}
\newcommand{\Fj}{F_j}
\newcommand{\Fia}{F_i^\alpha}
\newcommand{\Mab}{M^{\alpha\beta}}
\newcommand{\nosum}{\qquad\mbox{(no sum over $\alpha$)}}
\newcommand{\nosumab}{\qquad\mbox{(no sum over $\alpha$ or $\beta$)}}

\newcommand{\rma}{1/m_A}
\newcommand{\rmh}{1/m_H}
\newcommand{\rmo}{1/m_O}

\newcommand{\trma}{1/{\tilde m}_A}
\newcommand{\trmh}{1/{\tilde m}_H}
\newcommand{\trmo}{1/{\tilde m}_O}

\newcommand{\tFi}{{\tilde F}_i}
\newcommand{\tra}{{\tilde r}^\alpha}
\newcommand{\tda}{{\tilde d}^\alpha}
\newcommand{\trb}{{\tilde r}^\beta}
\newcommand{\tla}{{\tilde \lambda}_{\alpha}}
\newcommand{\tlb}{{\tilde \lambda}_{\beta}}
\newcommand{\tMab}{{\tilde M}^{\alpha\beta}}

\newcommand{\Rs}{R_s}
\newcommand{\Fs}{F_s}
\newcommand{\ls}{\lambda_s}
\newcommand{\Bs}{B_s}
\newcommand{\Ms}{M_s}
\newcommand{\ms}{m_s}

\newcommand{\Length}{\mbox{Length}}
\newcommand{\Force}{\mbox{Force}}
\newcommand{\Mass}{\mbox{Mass}}
\newcommand{\Time}{\mbox{Time}}

\newcommand{\bLs}{b}
\newcommand{\pshift}{\mbox{integrator\_momentum\_shift}}
\newcommand{\qshift}{\mbox{integrator\_position\_shift}}

\title{Notes on Mshake}
\author{John Salmon}
%\date{Feb 6, 2007}
\begin{document}
\maketitle

\section{Introduction}
The constraints that we typically solve for in MD, i.e., the $\HtO$ and $\AHN$ constraints,
are special cases of so-called ``holonomic'' constraints.  General formulas
for the equations of motion in the presence of holonomic constraints are well
known.
For the purposes of saving code-space in the GC, we'd like to abstract the
constraint-solving code as much as possible.  Therefore, it's helpful to
have the constraint algorithm written in a general form.  The most general
possible form is interesting theoretically, but is not particularly useful
for computation or programming.  Therefore this note works through the details
of the position-constraint algorithm for a class of constraints that includes
the $\AHN$ and $\HtO$ constraints as special cases, but that are less general
than fully holonomic constraints.  Furthermore, the GCs must work in
fixed-point and ``prefer'' to work with numbers whose magnitudes are
in the range $[-1, 1)$, so this note also works through the necessary
scalings to guarantee this.

\subsection{Notation}
The notational convention will be that greek indices ($\alpha, \beta, \gamma$, etc)
identify constraints, while roman indices ($i, j, k$, etc) represent coordinates.
Repeated indices follow the ``Einstein summation convention'' and imply
that the index should be summed.  I've tried as much as possible to
keep ``covariant'' and ``contravariant'' indices up and down respectively.
This is a ``no-op'' for Cartesian coordinate systems, but is a useful
convention to adopt in curvilinear coordinates as it tends to remind one
where the ``metric tensor'' needs to go.

\section{Holonomic constraints (the general case)}

``Holonomic constraints'' are written in the form:
\bd
\sigma^\alpha(q) = 0
\ed
where $\alpha$ is an index that identifies the constraint and q are the
coordinates in the state space.  The functions $\sigma^\alpha$ are arbitrary.
In fact, the constraints are specified only up to a constant scale factor.
The dynamics of two systems whose constraint functions
differ only by a constant non-zero factor are exactly the same.  

The standard approach is to augment the lagrangian, $L$, with a set of
lagrange multipliers, $\la$ i.e.,
\bd
L'(q, \dot q) = L + \lambda_\alpha \sigma^\alpha(q)
\ed

Discretizing in time, let $q^i(t)$ be the result of propagating the
unconstrained equations of motion forward in time by t.  Then the
discretized constrained equations result in:

\bd
q'^i(t) = q^i(t) + \half t^2 \minvij \la \djsa|_{q(t=0)} \label{eq:qprime}
\ed

Where $\minvij$ is the ``inverse mass matrix'' which is found
in the momentum part of the Hamiltonian, i.e.,
\bd
H = \half \minvij p_i p_j + V(q)
\ed

These equations behave exactly as if a ``constraint force'',
\bd
\Fi = \la \disa|_{q(t=0)}
\ed
were applied to each atom.  The total constraint force, $\Fi$ is a linear
combination of terms associated with each constraint $\disa$.  The magnitudes
of these constraint forces are adjusted by the lagrange multipliers, but
the directions are fixed by the functional form of the constraint and the
coordinates at t=0.

For a typical MD system, the inverse mass matrix is diagonal with
elements given by the inverse of the individual atomic masses, i.e.,
\bd
\minvij = diag(1/m_i)
\ed

The constraints are satisfied by solving for $\la$ such that
\bd
\sa(q'(t)) = 0 \label{eq:constraints}
\ed

Remarkably, if the unconstrained propagator for $q^i(t)$ is
symplectic, then the constrained propagator defined by
Eqs~\ref{eq:qprime} and~\ref{eq:constraints} is also symplectic.  Among
other things, it preserves ``phase space volume'', and is expected to have very
good long-term stability properties.  Symplecticity is fragile.  In
particular, it is crucial that the ``direction'' of each of the constraint
forces be defined by the partial derivatives of the constraint functions
with respect to the coordinates at t=0, and it is also crucial that the
constraint be solved as exactly as possible.  Allowing the state to drift
even slightly from the constraint surface affects symplecticity, and
hence the long-term stability of the trajectory -- not to mention the fact
that it violates the constraint!

Eq~\ref{eq:constraints} is a set of coupled 
nonlinear equations.  One way to solve
them is by the Newton-Raphson method, i.e., to start with a guess, $\la$ and
to refine it iteratively according to:
\bd
\la' = \la - \left(\frac{\partial \sbeta(q'(t))}{\partial \la}\right)^{-1} 
   \sbeta(q(t)) \label{eq:nr}
\ed

It is straightforward to apply the chain rule to simplify the partial
derivatives in Eq~\ref{eq:nr}, resulting in:

\bd
\la' = \la - (Z^{-1})_{\alpha\beta}\sbeta \label{eq:laprime}
\ed
where

\bd
\Zab = \half t^2 \disa(q'(t)) \minvij \djsb(q(0))
\ed

It is worth noting that $\Zab$ combines derivatives of the constraint
functions at the ``new new'' coordinates $q'(t)$ and at the ``old coordinates''
$q(0)$.  It is important to compute $\Zab$ in exactly this way.  Otherwise,
one will lose either symplecticity or the convergence properties of the
Newton-Raphson iteration, or both.

\section{Quadratic Holonomic Constraints}

In practice, constraints are not completely arbitrary functions, $\sa(q)$.
The constraints that we're interested in for MD are generally designed to
fix one or more bond-lengths to specific values.  We can write such
constraints as:

\bd
\ra = \Cai \qi \label{eq:radef}
\ed

\bd
\sa(q) = \half \left( (\ra)^2 - (\da)^2 \right) \label{eq:saq}
\ed

I.e., we define ``separation vectors'', $\ra$ as linear combinations
of the coordinates $\qi$, and we define the constraints to enforce the
condition that the magnitudes of the separation vectors are fixed at
predetermined values, $\da$.

The ``separation matrices'', $\Cai$ are very simple matrices
consisting of +1, -1 and 0.  For a constraint that fixes the bond
length between two atoms, $i_1$ and $i_2$, $\Cai$ is +1 when i is
equal to $i_1$ and -1 when i is equal to $i_2$, i.e.,
\bd
\ra = q_{i_1} - q_{i_2}
\ed
\bd
\Cai = \delta_{ii_1} - \delta_{ii_2}
\ed

For constraints defined this way, the general equations above simplify
considerably:

\bd
\disa = \Cai \ra \nosum
\ed
\bd
\Fi = \la \Cai \ra \label{eq:Fi}
\ed

In other words, the $\alpha$ constraint introduces a force directed
along $\ra$, the separation vector between the atoms connected by
$\alpha$ at t=0.  The force is felt with equal and opposite magnitude
by the two atoms connected by the constraint, and its magnitude is
given by the lagrange multiplier, $\la$, times the separation vector.

Solving the constraint equations by Newton's method 
requires us to evaluate the constraint
functions for trial values of $\la$.  I.e., we will need to
compute $\sa(q'(t))$, which obviously requires us to
compute:
\bea
\rprimea(t) &=& c \Cai \qprimei(t) \\
&=& c \Cai ( \qi(t) + \halftsq \minvij \Fj ) \\
&=& c \Cai ( \qi(t) + \halftsq \minvij \Cbj \lb \rb(0) ) \\
&=& \ra(t) + \halftsq \Cai \minvij \Cbj \lb \rb(0) 
\eea

Since Newton's method is iterative, it is useful to pre-compute
as much as possible outside the loop.  Define:
\bd
\Mab = \half \Cai \minvij \Cbj
\ed

Then we have:
\bd
\rprimea(t) = \ra(t) + t^2 \Mab \lb \rb(0) \label{eq:rprimea}
\ed

The $Z$ matrix used in the Newton-Raphson iteration, Eq~\ref{eq:laprime}, takes the form:
\bd
\Zab = t^2 {r'^{\alpha}}(t) \Mab \rb(0) \nosumab \label{eq:Zab}
\ed

A physical analogy is apt.  Just as we think of a bond as ``balls on a
spring'', we can think of a constrained bond as ``balls on a stick''.
If one pushes and pulls on the ends of the stick, the separation
between the ends never changes.  The stick provides a compressive or
tensile force that exactly balances the applied forces along the
direction of the stick.  Solving the constraint equations amounts
to solving for the compressive or tensile force supplied by the stick.
Solving for the force in a single stick is
easy - one just projects out the net applied forces along the
direction of the stick.  On the other hand, solving for the internal
forces in a network of sticks (think of a ball-and-stick molecular
model) requires solving a set of equations that result in no net
acceleration between the ends of any of the sticks.  These coupled
equations are exactly those solved by ``mshake'', and the results,
i.e., the lagrange multipliers times the separation vectors, $\la \ra$ are
exactly the compressive and tensile forces in the physical sticks.

The M matrix (and hence the Z matrix and its inverse) is blocked.
Non-zero matrix elements only connect constraints that share atoms.
In other words, the constraints form a graph, and the blocks
correspond to the maximally connected subgraphs.
Because the matrix is blocked, one can solve the equations
one block, i.e., one connected subgraph, at a time.
For the cases that arise in MD, the connected subgraphs
correspond to $\AHN$ and $\HtO$ constraints.  I.e., the N individual
constraints that form each $\AHN$ constraint are solved together for N
lagrange coefficients by iteratively constructing an NxN Z matrix,
inverting it, and applying the Newton-Raphson iteration.  Similarly,
the three lagrange multipliers corresponding to each $\HtO$ are found by
iteratively constructing a 3x3 Z matrix, inverting it and applying the
Newton-Raphson iteration.

\section{Specific cases}
It is helpful to see the equations ``unrolled'' for the specific
cases of interest to Anton, i.e., $\AHN$ and $\HtO$ constraints.

In contrast to Desmond, we will incorporate the assumptions that
all hydrogen atoms have the same mass and that all AH or OH bonds in
a constraint group have the same length.  The more general case
obviously follows from the equations above, but since we are pressed for
space on the GCs, it is worth exploiting obvious symmetries to save
storage where possible.

\subsection{$\AHthree$}
The connected subgraph consists of four atoms, with three bond-length constraints.  So $\Cai$ is
a 3x4 matrix:
\bd
\Cai = \left(
\begin{array}{cccc}
1 &-1 & 0 & 0 \\
1 & 0 &-1 & 0 \\
1 & 0 & 0 &-1 \\
\end{array}
\right)
\ed
The columns of $\Cai$ correspond to the A, $H_1$, $H_2$ and $H_3$ atoms,
respectively.  The rows of $\Cai$ correspond to the $A-H_1$, $A-H_2$ and
$A-H_3$ constraints, respectively.

$\Mab$ is a 3x3 matrix:
\bd
\Mab = \half\left(
\begin{array}{cccc}
\rma +\rmh & \rma & \rma \\
\rma & \rma + \rmh & \rma \\
\rma & \rma & \rma+\rmh \\
\end{array}
\right) \label{eq:Mabahn}
\ed

It's worth noting that the mass of the
A atom is typically 12 (Carbon), 14 (Nitrogen) or 16 (Oxygen) times
larger than the mass of a Hydrogen.  Therefore, the M matrix (and
hence the Z matrix) is diagonally dominant because the large $\rmh$ terms
appear only on the diagonal.  We can expect all our favorite formulas 
and methods from
Linear Algebra to be well-behaved, and we can expect that
an initial guess that ignores $\rma$ completely will provide us with
a good starting point for our Newton-Raphson iteration.

\subsection{$\HtO$}
The water constraint consists 
of three bond-length constraints connecting three atoms.
Therefore, $\Cai$ is a 3x3 matrix:
\bd
\Cai = \left(
\begin{array}{ccc}
1 &-1 & 0 \\ % H1-O
0 & 1 &-1 \\ % O-H2
-1& 0 &1 \\ %  H2-H1
\end{array}
\right)
\ed

The columns of $\Cai$ correspond to the $H_1$, $O$, and $H_2$ atoms,
respectively.  The rows of $\Cai$ correspond to the $H_1-O$, $O-H_2$ and
$H_2-H_1$ constraints, respectively.

$\Mab$ is a 3x3 matrix:
\bd
\Mab = \half \left(
\begin{array}{ccc}
\rmo +\rmh & -\rmo        & -\rmh \\
-\rmo      & \rmo + \rmh  & -\rmh \\
-\rmh & -\rmh & 2/m_H \\
\end{array}
\right)
\label{eq:Mabh2o}
\ed

The Z matrix is again diagonally dominant, but the argument now relies
on the separation vectors not being colinear.  This is obviously true
for actual water molecules (none of the $H_1-O$, $O-H_2$ nor $H_2-H_1$ bonds are
even close to colinear), so all our favorite formulas and methods from
linear algebra are again safe to use.

\subsection{$\AHtwo$ and $\AHone$}
The $\Cai$ and $\Mab$ matrices for $\AHtwo$ and $\AHone$ constraints are
exactly the same as for the $\AHthree$ case, with one or two of the right-most
columns and bottom-most rows removed.  However, on Anton, in the
interest of code-space, we will implement the Newton-Raphson iteration
only once for 3x3 Z matrices.  For $\AHtwo$ and $\AHone$, we'll fill out the 3x3
matrix and the ``right-hand side'' with zeros and ones so that the ``extra'' 
rows and columns don't change the results for the real lagrange multipliers.

\section{Scaling for Anton}
As usual, some thought is required to prescale values so that Anton's
fixed-point arithmetic neither overflows nor loses precision unnecessarily.

\subsection{Units}

We'll refer to the quantities in the equations above as having
``textbook units'', i.e., they're dimensional quantities with units of
Length, inverse Mass, Force, etc.  For implementation in Anton, we
need to find a rescaling of the equations so that the quantities
manipulated in GC registers are dimensionless
and of bounded magnitude in the range [-1, 1). \footnote{Boundedness
is not strictly necessary.  Anton can, if necessary, represent a wider range
of values by keeping an explicit exponent in a separate register.
We prefer to avoid this because it is consumes both registers
and instructions.
}

First, let's review the ``textbook units'' in the equations above.
The notation ``$:=$'' should be read as ``has dimensions of''.  Then:

\bea
\qi &:=& \Length \\
\ra &:=& \Length \\
\sa &:=& \Length^2 \\
\da &:=& \Length \\
\la &:=& \Force / \Length \\
\Mab &:=& 1 /{\Mass } \\
\Zab &:=& \Length^3 /{\Force } := \Time^2 \Length^2/\Mass
\eea

\subsection{Scale factors}
When working with scaled quantities the notational convention
is that a subscript 's' indicates a dimensional scale factor,
a superscript tilde indicates a normalized, dimensionless
quantity suitable for manipulation directly by the GCs and an
unadorned symbol represents a dimensional ``textbook'' quantity
like force, energy, distance, inverse mass matrix, etc.  I.e.,
% \bd
% X = \tilde X X_s \label{eq:tildescale}
% \ed

Rewriting ``textbook'' equations in terms of normalized quantities
is straightforward.  One simply replaces $X$ with $\tilde X X_s$ 
and rearranges terms to isolate normalized quantities.  Then one
specifies the scale factors to make the resulting computation most
convenient.  Multiplicative factors that evaluate to 1 are best.
Powers of 2 are almost as good.

\subsubsection{Normalizing the constraints:  $r_s$, $d_s$, and $\sigma_s$}
Let's start with Eq~\ref{eq:saq}, rewritten as:
\bd
\tsa = \half \left( (\tra)^2 (r_s/\sigma_s)^2 - (\tda)^2 (d_s^2/\sigma_s^2) \right)
\ed

The above equation suggests that we define:
\bd
r_s = \sigma_s = d_s
\ed
in which case:
\bd
\tsa = \half \left( (\tra)^2  - (\tda)^2 \right) \label{eq:tsa}
\ed

The bond calculation in the GCs makes frequent use of the ``bond scale'',
$\Bs$ defined in terms of the global 
position scale, $\Rs$ and the ``bondLshift'', $\bLs$.
\bd
\Bs = \Rs 2^{-\bLs}
\ed
The bondLshift, $\bLs$, is the smallest integer that guarantees
that the separation between all bonded pairs of atoms is less
than the bond scale.  The same considerations that make the
bond scale useful in bond calculations make it useful for constraints
as well.  In particular, if we set $r_s = \Bs$, then the $\tra$ and
$\tsa$ will be representable with the most possible significant
bits.  Therefore, we set:
\bd
r_s = \sigma_s = d_s = \Bs = \Rs 2^{-\bLs}\label{eq:rs}
\ed

Note that Eqs~\ref{eq:rs} and~\ref{eq:radef} imply:
\bd 
\tra = 2^\bLs \Cai \tqi \label{eq:tra}
\ed
I.e., to
compute the normalized constraint vectors, $\tra$, we subtract atomic
coordinates (the $\pm 1$ terms in $\Cai$) and shift the result left
by $\bLs$, just as we do when calculating bond terms.

\subsubsection{Normalizing the lagrange multipliers: $\ls$}
We can follow the same approach for Eq~\ref{eq:Fi}, giving:
\bd
\tFi = \Cai \tla \tra ( \ls \Bs / \Fs )
\ed
which simplifies to:
\bd
\tFi = \Cai \tla \tra \label{eq:tFi}
\ed
if we define:
\bd
\ls = \Fs/\Bs \label{eq:ls}
\ed

\subsubsection{Normalizing the $Z$ matrix:  $Z_s$}
Equations~\ref{eq:laprime} and~\ref{eq:ls} give:
\bd
\tla' = \tla - (\tZ^{-1})_{\alpha\beta}\tsbeta (\sigma_s / Z_s \ls)
      = \tla - (\tZ^{-1})_{\alpha\beta}\tsbeta (\Bs^3 / Z_s F_s)
\ed
simplifying to:
\bd
\tla' 
      = \tla - (\tZ^{-1})_{\alpha\beta}\tsbeta \label{eq:tlaprime}
\ed
if we define:
\bd
Z_s = \Bs^3/F_s
\ed

\subsubsection{Normalizing the $M$ matrix: $M_s$}
Equation~\ref{eq:Zab} may be rewritten as:
\bd
\tZab = {{\tilde r}'^{\alpha}}(t) \tMab \trb(0) (t^2 M_s \Bs^2 / Z_s)\nosumab
\ed
which simplifies to:
\bd
\tZab = {{\tilde r}'^{\alpha}}(t) \tMab \trb(0) \nosumab \label{eq:tZab}
\ed
if we define:
\bd
M_s = Z_s / t^2 \Bs^2 = B_s / F_s t^2 \label{eq:Ms}
\ed

The force, bond, and mass scales and the timestep, t, are not independent.
In order to replace multiplicative factors with shifts
in the integrator, the clink infrastructure guarantees
that:
\begin{align}
t \Fs / P_s &= 2^{\pshift} \\
t P_s / R_s m_s &= 2^{\qshift}
\end{align}
where $P_s$ is the momentum scale and $m_s$ is the mass scale.  It follows that:
\begin{align}
M_s &= 2^{-k} / m_s \\
k &= \pshift + \qshift + \bLs \\
\end{align}
A typical system has $\qshift=-9$,
${\pshift=-3}$, and $\bLs=8$, in which case, $k=-4$.

For the $\HtO$ constraint, Eq~\ref{eq:Mabh2o} gives:
\bd
\tMab = \Mab / M_s = 2^{k-1}
\left(
\begin{array}{ccc}
\trmo +\trmh & -\trmo        & -\trmh \\
-\trmo      & \trmo + \trmh  & -\trmh \\
-\trmh & -\trmh & 2/{\tilde m}_H \\
\end{array}
\right)
\label{eq:tMabh2o}
\ed

For the $\AHN$ constraint, Eq~\ref{eq:Mabahn} gives:
\bd
\tMab = 2^{k-1}\left(
\begin{array}{cccc}
\trma +\trmh & \trma & \trma \\
\trma & \trma + \trmh & \trma \\
\trma & \trma & \trma+\trmh \\
\end{array}
\right) \label{eq:tMabahn}
\ed

In both cases, the $\Mab$ matrix is a trivial combination of shifted
normalized reciprocal masses already known to the integrator.

Finally, equation~\ref{eq:rprimea} may now be rewritten using
Eqs~\ref{eq:ls} and~\ref{eq:Ms}:
\bd
\trprimea(t) = \tra(t) + \tMab \tlb \trb(0) \label{eq:trprimea}
\ed

\section{Putting it all together}

Constraints depend on two sets of coordinates: $\tqi(0)$ and
$\tqi(t)$.  $\tqi(t)$ are obtained by using the unconstrained
integrator, but are {\em only} used as input to the constraint solver.
They are never used as a bona fide coordinates of the system.  Similarly,
any intermediate values used in the computation of $\tqi(t)$, e.g.,
$\tilde p_i(t/2)$ should not to be used as an intermediate
state of the system and should be discarded.

The output of the constraint
calculation is a set of normalized constraint forces, $\tFi$ to be added to
the total force on each atom.  With this force incorporated, the
position and momentum updates can be carried out, 
resulting in a new state, $(\tqprimei(t), {\tilde p}'_i(t))$ 
that satisfies the constraint.  Global quantities (e.g.,
kinetic energy) should be computed from the new state resulting
from the integration of the constrained forces.  Intermediate
states obtained from unconstrained forces should be discarded.

\begin{algorithmic}
\REQUIRE Inputs: $\tqi(0), \tqi(t), \tMab, (\tda)^2, \Cai$
\STATE Compute $\tra(0)$ // Eq~\ref{eq:tra}
\STATE Compute $\tra(t)$ // Eq~\ref{eq:tra}
\STATE Set $\tla = 0$
\STATE Set $\trprimea(t) = \tra(t)$ // follows from $\tla=0$
\STATE Compute $\tsa(\trprimea(t))$ // Eq~\ref{eq:tsa}
\WHILE{ $\tsa(\trprimea(t))$ excedes threshold }
\STATE Compute $\tZab$        // Eq~\ref{eq:tZab}
\STATE Compute $\tla'$        // Eq~\ref{eq:tlaprime}, Cramer's rule
\STATE Compute $\trprimea(t)$ // Eq~\ref{eq:trprimea}
\STATE Compute $\tsa(\trprimea(t))$ // Eq~\ref{eq:tsa}
\ENDWHILE
\STATE Compute $\tFi$         // Eq~\ref{eq:tFi}
\end{algorithmic}

\section{Virials and NPT ensembles}
More work is needed here, but my intuition is that if we compute the
``atomic'' virial with the constraint force included in the total
force, the result will be the same as if we had computed the
molecular virial.  I.e., the constraint force exactly balances the
external forces to make the net intramolecular acceleration zero.  This
is exactly (?) the same as the correction to the atomic virial based
on intra-molecular forces.

\end{document}
