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.