RMSD (root mean square deviation) and RMSF (root mean square fluctuation) are common measures of biomolecules’ spatial variations in a molecular dynamics (MD) simulation. RMSD describes the molecule’s overall discrepancy with respect to a reference conformation. RMSF, on the other hand, is a per atom quantity describing the atom’s variation over the whole trajectory. Their definitions are

\[\begin{align} \text{RMSD}(t) \equiv& \sqrt{\frac{1}{N_a}\sum_i \|\mathbf r_i(t) - \mathbf r_i'\|^2} \\ \text{RMSF}_i\equiv&\sqrt{\left<\|\mathbf r_i(t) - \left< \mathbf r_i\right> \|^2 \right>} \end{align}\]

where \(N_a\) is the number of atoms of interest, \(\mathbf r_i\) is the location of the \(i\)‘th atom, \(\mathbf r_i'\) is its reference location, \(tr\) is trace operation, \(\otimes\) is the outer product, and \(\left<\cdot\right>\) denotes temporal average. In the case of discrete time dynamics, e.g., MD simulation,

\[\left<f(t)\right>\equiv \frac{1}{N_t}\sum_k f(t_k)\]

where \(N_t\) is the number of time points.

Often times one is interested in the conformational changes after aligning a group of atoms. For example, the side chain’s RMSD and RMSF with the protein backbone aligned.

Such alignment is naturally formulated as a minimization problem

\[\min f(R, \mathbf t) \equiv \sum_{i=1}^{N_a} w_i\|(R\mathbf r_i + \mathbf t) - \mathbf r_i'\|^2\]

where \(w_i\) is a weighting factor (for example, mass or charge), \(R\) and \(\mathbf t\) are the rotation matrix and translation vector to be determined. The solution to the unweighted version is commonly known as Kabsch algorithm.

To get \(R\) and \(\mathbf t\), we solve the two stationary conditions

\[\begin{align} \frac{\partial{f}}{\partial\mathbf t} =& 0 \\ \frac{\partial{f}}{\partial R} =& 0 \end{align}\]

The following vector equality will be handy

\[\nabla_A tr AB = B^T\]

The first condition gives rise to

\[\mathbf t = \frac{\sum_i w_i \mathbf r_i'}{\sum_i w_i} - R\frac{\sum_i w_i \mathbf r_i}{\sum_i w_i}\equiv \overline{\mathbf r}' - R\overline{\mathbf r}\]

Plugging into \(f(R, \mathbf t)\), we get

\[f(R) = \sum_i w_i \|R(\mathbf r_i - \overline{\mathbf r}) - (\mathbf r_i' -\overline{\mathbf r}')\|^2\equiv \sum_i w_i \|R\mathbf p_i - \mathbf q_i \|^2\]

Here we define new sets of vectors \(\mathbf p_i\) and \(\mathbf q_i\) to simplify the notation. The two groups of atoms have their geometric centers (defined with respect to the weights \(w_i\)) removed.

Then \(f(R)\) can be written as

\[f(R) = tr WP^TP + tr WQ^TQ - 2tr WQ^TRP\]

where \(P\) and \(Q\) are both \(3\times N_a\) matrices whose columns are the atom positions. Here \(R^TR=1\) is used.

Since only the last term depends on \(R\), the second stationary condition won’t give us \(R\) directly. But we can see that the original minimization problem is equivalent to the maximization problem of

\[\max tr R P W Q^T\equiv \max tr R S^T\]

It can be further simplified by applying singular decomposition on the covariance matrix

\[S \equiv QWP^T = U\Sigma V^T\]

We have

\[\max tr \left(U^T R V\right) \Sigma\]

It’s easy to see that the maximum value is attained when

\[R = UV^T\]

Note that this \(R\) matrix is not guaranteed to be a rotation. It could be an improper rotation. Depending on the application, one may or may not want to allow reflection in the solution. I will leave it as an exercise to you to restrict the answer to proper rotations.