Get least RMSD#

Getting the least RMSD of a molecular system from a reference structure after optimal superimposition.

The function molsysmt.structure.get_least_rmsd() calculates the root-mean-square deviation (RMSD) between structures after finding the optimal translation and rotation that minimizes coordinate differences (using the quaternion-based algorithm by Coutsias et al.).

Added in version 1.0.0.

API documentation

Follow this link for a detailed description of the input arguments, raised errors, and returned objects of this function: molsysmt.structure.get_least_rmsd().

Basic usage#

Let’s show how to calculate the least RMSD along a pentalanine simulation trajectory with respect to the initial configuration:

import molsysmt as msm
import numpy as np
import matplotlib.pyplot as plt

molsys = msm.convert(msm.systems['pentalanine']['traj_pentalanine.h5'], to_form='molsysmt.MolSys')

Note

Demo Systems Catalog

This tutorial uses demonstration datasets provided by MolSysMT. To explore the full catalog of bundled systems, forms, and file paths, visit the Demo Systems guide.

We compute both the unaligned standard RMSD and the optimal least RMSD (LRMSD) for the backbone atoms:

rmsd = msm.structure.get_rmsd(molsys, selection='backbone', reference_structure_index=0)
lrmsd = msm.structure.get_least_rmsd(molsys, selection='backbone', reference_structure_index=0)
time = msm.get(molsys, element='system', time=True)
time_ns = msm.pyunitwizard.convert(time, to_unit='ns')

We plot the unaligned RMSD versus the optimal least RMSD across the entire simulation:

plt.figure(figsize=(6, 4))
plt.plot(time_ns, rmsd, label='RMSD (unaligned)', color='lightblue')
plt.plot(time_ns, lrmsd, label='Least RMSD (optimal)', color='blue')
plt.xlabel('Time [ns]')
plt.ylabel('RMSD [nm]')
plt.title('Backbone RMSD vs. Least RMSD')
plt.legend()
plt.show()

../../../../_images/50eab50931ced5babe067371a70bac13a3aee0d503cc93300cdc19d82ff53f7b.png

Comparing against an external reference structure#

We can compute the least RMSD with respect to an external reference molecular system using reference_molecular_system:

native_ref = msm.extract(molsys, structure_indices=0)
lrmsd_ref = msm.structure.get_least_rmsd(molsys, selection='backbone', reference_molecular_system=native_ref)
print('Mean LRMSD to native reference:', np.mean(lrmsd_ref))

Mean LRMSD to native reference: 0.2224203571949468 nanometer