Get forces#

Extracting atomic force vectors and scalar magnitudes from molecular mechanics force fields.

The function molsysmt.molecular_mechanics.get_forces() evaluates the negative gradient of the potential energy with respect to atomic coordinates (\(-\nabla U\)).

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.molecular_mechanics.get_forces().

Note

Molecular Mechanics Backend Engine

Molecular mechanics and force evaluations in MolSysMT are executed using backend simulation engines (currently OpenMM by default). MolSysMT intentionally unifies the Pythonic interface and data structures while delegating energy evaluations, parameterizations, and gradient calculations to high-performance, battle-tested molecular simulation engines.

Basic usage#

Let’s evaluate the residual atomic force vectors on a relaxed Trp-Cage miniprotein:

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

molsys = msm.convert(msm.systems['Trp-Cage']['1l2y.h5msm'], structure_indices=0)
molsys_min = msm.molecular_mechanics.potential_energy_minimization(molsys)
context = msm.convert(molsys_min, to_form='openmm.Context')

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 extract the 3D atomic force vectors acting on the system:

forces = msm.molecular_mechanics.get_forces(context)
print('Forces array shape (n_atoms, 3):', forces.shape)
print('First 5 atomic force vectors (kJ/(mol*nm)):\n', forces[:5])

Forces array shape (n_atoms, 3): (304, 3)
First 5 atomic force vectors (kJ/(mol*nm)):
 [[-203.52891418538445 -59.48664780872133 -41.75353447795948] [5.266354881566144 19.56653812680109 12.627775171915474] [-12.145142046770047 -11.768486409033017 5.129915645151726] [-1.461332942163267 -4.805862072016566 9.03546569932422] [0.06741908959327958 12.002086706279599 22.31171933490519]] kilojoule / mole / nanometer

Force magnitudes and residue profiles#

Setting magnitude=True calculates the scalar Euclidean norm of the forces (\(\|\mathbf{F}_i\|\)) for each atom:

force_mags = msm.molecular_mechanics.get_forces(context, magnitude=True)
mag_vals = puw.get_value(force_mags, to_unit='kJ/(mol*nm)')

fig, ax = plt.subplots(figsize=(8, 3.5))
ax.plot(range(len(mag_vals)), mag_vals, color='firebrick', marker='o', markersize=3, linewidth=1.2)
ax.set_title('Atomic Force Magnitudes on Minimized Trp-Cage')
ax.set_xlabel('Atom Index')
ax.set_ylabel('Force Magnitude (kJ/(mol$\cdot$nm))')
ax.grid(True, linestyle='--', alpha=0.4)
plt.tight_layout()
plt.show()

../../../../_images/4df72da060cb4dae52efbdcb26f6e99df0377015ca91130d873d8d0f9390f3c0.png

Visualizing atomic force magnitudes in 3D#

We can color the 3D molecular structure by atomic force magnitudes to inspect the spatial distribution of mechanical forces:

# view = msm.view(molsys_min)
# view.whole.set_color_by_values(puw.get_value(force_mags), element='atom', palette='viridis')
# view