molecular_simulations.simulate.omm_simulator module

OpenMM molecular dynamics simulation module.

This module provides classes for running molecular dynamics simulations using OpenMM with AMBER or CHARMM force fields. It supports both explicit and implicit solvent simulations, with built-in equilibration protocols and production MD.

Classes:

Simulator: Main class for explicit solvent OpenMM simulations. ImplicitSimulator: Simulator for implicit solvent (GB) simulations. CustomForcesSimulator: Simulator with user-defined custom forces. Minimizer: Simple energy minimization utility.

class molecular_simulations.simulate.omm_simulator.Simulator(path, top_name=None, coor_name=None, out_path=None, ff='amber', heat_steps=100000, equil_steps=1250000, prod_steps=250000000, prod_dt_in_ps=0.004, hydrogen_mass_repartitioning=1.5, n_equil_cycles=3, temperature=300.0, eq_reporter_frequency=1000, prod_reporter_frequency=10000, platform='CUDA', device_ids=[0], force_constant=10.0, params=None, membrane=False)[source]

Bases: object

Class for performing OpenMM simulations on AMBER FF inputs.

Inputs must conform to naming conventions found below in the init. Supports explicit solvent simulations with PME electrostatics and hydrogen mass repartitioning for longer timesteps.

Parameters:

  • path (Path | str) – Path to simulation inputs, same as output path.

  • top_name (str | None) – Optional topology file name. If not provided, assumes ‘system.prmtop’.

  • coor_name (str | None) – Optional coordinate file name. If not provided, assumes ‘system.inpcrd’.

  • out_path (Path | None) – Optional output path for simulation outputs. If not provided, uses the same path as inputs.

  • ff (str) – Force field to use, either ‘amber’ or ‘charmm’.

  • heat_steps (int) – Number of heating timesteps. Defaults to 100,000 (200 ps at 2 fs timestep).

  • equil_steps (int) – Number of equilibration timesteps. Defaults to 1,250,000 (2.5 ns at 2 fs timestep).

  • prod_steps (int) – Number of production timesteps. Defaults to 250,000,000 (1 µs at 4 fs timestep).

  • n_equil_cycles (int) – Number of unrestrained equilibration cycles after restraint relaxation. Defaults to 3.

  • temperature (float) – Simulation temperature in Kelvin. Defaults to 300.0.

  • eq_reporter_frequency (int) – Reporter frequency during equilibration in timesteps. Defaults to 1,000.

  • prod_reporter_frequency (int) – Reporter frequency during production in timesteps. Defaults to 10,000.

  • platform (str) – OpenMM platform to use. Options: ‘CUDA’, ‘CPU’, ‘OpenCL’. Defaults to ‘CUDA’.

  • device_ids (list[int]) – List of GPU device IDs to use. Defaults to [0].

  • force_constant (float) – Harmonic restraint force constant in kcal/mol*Ų. Defaults to 10.0.

  • params (str | list[str] | None) – Optional list of CHARMM parameter files for loading from psf/pdb file using CHARMM36m forcefield.

  • membrane (bool) – Whether this is a membrane system requiring anisotropic pressure coupling. Defaults to False.

path

Path object for simulation directory.

top_file

Path to topology file.

coor_file

Path to coordinate file.

temperature

Simulation temperature.

simulation

OpenMM Simulation object (available after production).

Example

>>> sim = Simulator(
...     path='./simulation',
...     equil_steps=500_000,
...     prod_steps=50_000_000,
...     device_ids=[0, 1],
... )
>>> sim.run()
__init__(path, top_name=None, coor_name=None, out_path=None, ff='amber', heat_steps=100000, equil_steps=1250000, prod_steps=250000000, prod_dt_in_ps=0.004, hydrogen_mass_repartitioning=1.5, n_equil_cycles=3, temperature=300.0, eq_reporter_frequency=1000, prod_reporter_frequency=10000, platform='CUDA', device_ids=[0], force_constant=10.0, params=None, membrane=False)[source]
setup_barostat(is_membrane_system)[source]

Configure the barostat based on system type.

Chooses the correct barostat for the system based on whether or not there is a membrane present. Membrane systems use MonteCarloMembraneBarostat with anisotropic pressure coupling.

Parameters:

is_membrane_system (bool) – True if this is a membrane-containing system.

Return type:

None

load_system()[source]

Load the molecular system based on force field type.

Dispatches to the appropriate file loader based on the specified force field (AMBER or CHARMM).

Return type:

System

Returns:

OpenMM System object configured with the appropriate force field.

Raises:

AttributeError – If an invalid force field type is specified.

load_amber_files()[source]

Build an OpenMM system from AMBER prmtop/inpcrd files.

Uses PME for electrostatics with a 1 nm non-bonded cutoff and 1.5 amu hydrogen mass repartitioning for longer timesteps.

Return type:

System

Returns:

OpenMM System configured for AMBER force field.

load_charmm_files()[source]

Build an OpenMM system from CHARMM psf/pdb files.

Uses PME for electrostatics with a 1.2 nm non-bonded cutoff.

Return type:

System

Returns:

OpenMM System configured for CHARMM force field.

setup_sim(system, dt)[source]

Build OpenMM Simulation and Integrator objects.

Creates a LangevinMiddleIntegrator with the specified timestep and builds the Simulation object with the configured platform.

Parameters:

  • system (System) – OpenMM System object to simulate.

  • dt (float) – Integration timestep in picoseconds.

Return type:

tuple[Simulation, LangevinMiddleIntegrator]

Returns:

Tuple containing (Simulation, LangevinMiddleIntegrator) objects.

run()[source]

Execute the full simulation workflow.

Determines whether to restart based on existing equilibration outputs. Checks production log to resume from checkpoints if available. Runs equilibration if needed, then production MD.

Return type:

None

equilibrate()[source]

Run the equilibration protocol.

Performs energy minimization, slow heating, and gradual restraint relaxation followed by unrestrained NPT equilibration.

Return type:

Simulation

Returns:

Equilibrated OpenMM Simulation object.

production(chkpt, restart=False)[source]

Run production molecular dynamics.

Loads a new system with barostat, loads checkpoint, attaches reporters, and runs production MD for the specified number of steps.

Parameters:

  • chkpt (Path | str) – Path to checkpoint file (equilibration or previous production).

  • restart (bool) – If True, append to existing log/DCD files instead of overwriting. Defaults to False.

Return type:

None

load_checkpoint(simulation, checkpoint)[source]

Load a checkpoint into the simulation.

Parameters:

  • simulation (Simulation) – OpenMM Simulation object to load checkpoint into.

  • checkpoint (Path | str) – Path to OpenMM checkpoint file.

Return type:

Simulation

Returns:

Simulation with restored positions, velocities, and step count.

attach_reporters(simulation, dcd_file, log_file, rst_file, restart=False)[source]

Attach trajectory, logging, and checkpoint reporters.

Parameters:

  • simulation (Simulation) – OpenMM Simulation object.

  • dcd_file (Path | str) – Path for DCD trajectory output.

  • log_file (Path | str | TextIOWrapper) – Path for state data log output.

  • rst_file (Path | str) – Path for checkpoint file output.

  • restart (bool) – If True, append to existing DCD file. Defaults to False.

Return type:

Simulation

Returns:

Simulation with reporters attached.

get_restraint_indices(addtl_selection='')[source]

Get atom indices for harmonic restraints.

Uses MDAnalysis to select protein/nucleic acid backbone atoms. Additional selections can be included for ligand restraints.

Parameters:

addtl_selection (str) – Additional MDAnalysis selection string to include in restraints. Defaults to empty string.

Return type:

list[int]

Returns:

List of atom indices to restrain.

check_num_steps_left()[source]

Check production log to determine remaining simulation steps.

Reads the log file to find the last completed step, then calculates the remaining steps needed. Also handles duplicate frames that may occur when restarting from checkpoints (since checkpoint frequency may not align with reporter frequency).

Return type:

None

static add_backbone_posres(system, positions, atoms, indices, restraint_force=10.0)[source]

Add harmonic position restraints to selected atoms.

Parameters:

  • system (System) – OpenMM System object.

  • positions (Any) – Position array for all atoms (with units).

  • atoms (Any) – List of OpenMM Atom objects from topology.

  • indices (list[int]) – List of atom indices to restrain.

  • restraint_force (float) – Force constant in kcal/mol*Ų. Defaults to 10.0.

Return type:

System

Returns:

Copy of System with harmonic restraints added.

class molecular_simulations.simulate.omm_simulator.ImplicitSimulator(path, top_name=None, coor_name=None, out_path=None, ff='amber', equil_steps=1250000, prod_steps=250000000, n_equil_cycles=3, temperature=300.0, eq_reporter_frequency=1000, prod_reporter_frequency=10000, platform='CUDA', device_ids=[0], force_constant=10.0, implicit_solvent=openmm.app.GBn2, solute_dielectric=1.0, solvent_dielectric=78.5, **kwargs)[source]

Bases: Simulator

Simulator for implicit solvent (Generalized Born) simulations.

Inherits from Simulator and overloads methods for implicit solvent compatibility. Uses GBn2 model by default with ionic screening.

Parameters:

  • path (Path | str) – Path to simulation inputs, same as output path.

  • top_name (str | None) – Optional topology file name. If not provided, assumes ‘system.prmtop’.

  • coor_name (str | None) – Optional coordinate file name. If not provided, assumes ‘system.inpcrd’.

  • out_path (Path | None) – Optional output path for simulation outputs.

  • ff (str) – Force field to use, either ‘amber’ or ‘charmm’.

  • equil_steps (int) – Number of equilibration timesteps. Defaults to 1,250,000.

  • prod_steps (int) – Number of production timesteps. Defaults to 250,000,000.

  • n_equil_cycles (int) – Number of unrestrained equilibration cycles.

  • temperature (float) – Simulation temperature in Kelvin. Defaults to 300.0.

  • eq_reporter_frequency (int) – Reporter frequency during equilibration.

  • prod_reporter_frequency (int) – Reporter frequency during production.

  • platform (str) – OpenMM platform to use. Defaults to ‘CUDA’.

  • device_ids (list[int]) – List of GPU device IDs to use. Defaults to [0].

  • force_constant (float) – Harmonic restraint force constant in kcal/mol*Ų.

  • implicit_solvent (Any) – GB model to use. Defaults to GBn2.

  • solute_dielectric (float) – Solute dielectric constant. Defaults to 1.0.

  • solvent_dielectric (float) – Solvent dielectric constant. Defaults to 78.5.

Example

>>> sim = ImplicitSimulator(
...     path='./simulation', implicit_solvent=GBn2, prod_steps=100_000_000
... )
>>> sim.run()
__init__(path, top_name=None, coor_name=None, out_path=None, ff='amber', equil_steps=1250000, prod_steps=250000000, n_equil_cycles=3, temperature=300.0, eq_reporter_frequency=1000, prod_reporter_frequency=10000, platform='CUDA', device_ids=[0], force_constant=10.0, implicit_solvent=openmm.app.GBn2, solute_dielectric=1.0, solvent_dielectric=78.5, **kwargs)[source]
load_amber_files()[source]

Build an OpenMM system with implicit solvent.

Uses the specified GB model with ionic screening for 150mM salt.

Return type:

System

Returns:

OpenMM System configured for implicit solvent simulation.

equilibrate()[source]

Run reduced equilibration protocol for implicit solvent.

Due to faster convergence with implicit solvent, uses a simplified equilibration protocol compared to explicit solvent.

Return type:

Simulation

Returns:

Equilibrated OpenMM Simulation object.

production(chkpt, restart=False)[source]

Run production MD for implicit solvent (no barostat).

Parameters:

  • chkpt (Path | str) – Path to checkpoint file.

  • restart (bool) – If True, append to existing files. Defaults to False.

Return type:

None

class molecular_simulations.simulate.omm_simulator.CustomForcesSimulator(path, custom_force_objects, equil_steps=1250000, prod_steps=250000000, n_equil_cycles=3, reporter_frequency=1000, platform='CUDA', device_ids=[0], equilibration_force_constant=10.0)[source]

Bases: Simulator

Simulator with user-defined custom forces.

Inherits from Simulator and provides injection of custom OpenMM forces for enhanced sampling or specialized simulations.

Parameters:

  • path (Path | str) – Path to simulation inputs.

  • custom_force_objects (list) – List of OpenMM Force objects to add to system.

  • equil_steps (int) – Number of equilibration timesteps. Defaults to 1,250,000.

  • prod_steps (int) – Number of production timesteps. Defaults to 250,000,000.

  • n_equil_cycles (int) – Number of unrestrained equilibration cycles.

  • reporter_frequency (int) – Reporter frequency in timesteps. Defaults to 1,000.

  • platform (str) – OpenMM platform to use. Defaults to ‘CUDA’.

  • device_ids (list[int]) – List of GPU device IDs. Defaults to [0].

  • equilibration_force_constant (float) – Restraint force constant in kcal/mol*Ų.

Example

>>> from openmm import CustomBondForce
>>> force = CustomBondForce('0.5*k*(r-r0)^2')
>>> force.addGlobalParameter('k', 1000)
>>> force.addGlobalParameter('r0', 0.3)
>>> sim = CustomForcesSimulator('./sim', [force])
>>> sim.run()
__init__(path, custom_force_objects, equil_steps=1250000, prod_steps=250000000, n_equil_cycles=3, reporter_frequency=1000, platform='CUDA', device_ids=[0], equilibration_force_constant=10.0)[source]
load_amber_files()[source]

Load OpenMM system and add custom forces.

Return type:

System

Returns:

OpenMM System with custom forces added.

add_forces(system)[source]

Add all custom forces to the system.

Parameters:

system (System) – OpenMM System object.

Return type:

System

Returns:

System with all custom forces added.

class molecular_simulations.simulate.omm_simulator.Minimizer(topology, coordinates, out='min.pdb', platform='OpenCL', device_ids=[0])[source]

Bases: object

Simple energy minimization utility for molecular structures.

Supports AMBER, GROMACS, and PDB input formats. Performs energy minimization and writes out the minimized structure as a PDB file.

Parameters:

  • topology (Path | str) – Path to topology file (prmtop, top, or pdb).

  • coordinates (Path | str) – Path to coordinate file (inpcrd, gro, or same pdb).

  • out (Path | str) – Output filename for minimized structure. Defaults to ‘min.pdb’.

  • platform (str) – OpenMM platform to use. Defaults to ‘OpenCL’.

  • device_ids (list[int] | None) – List of GPU device IDs, or None for CPU. Defaults to [0].

Example

>>> minimizer = Minimizer(
...     topology='system.prmtop',
...     coordinates='system.inpcrd',
...     out='minimized.pdb',
... )
>>> minimizer.minimize()
__init__(topology, coordinates, out='min.pdb', platform='OpenCL', device_ids=[0])[source]
minimize()[source]

Perform energy minimization and save structure.

Loads the system, runs OpenMM energy minimization, and writes the final coordinates to a PDB file.

Return type:

None

load_files()[source]

Load system based on topology file extension.

Return type:

System

Returns:

OpenMM System object.

Raises:

FileNotFoundError – If no valid input files are found.

load_amber()[source]

Load AMBER input files into OpenMM System.

Return type:

System

Returns:

OpenMM System object from AMBER files.

load_gromacs()[source]

Load GROMACS input files into OpenMM System.

Note

This method is untested and may require adjustments.

Return type:

System

Returns:

OpenMM System object from GROMACS files.

load_pdb()[source]

Load PDB file into OpenMM System.

Uses AMBER14 force field with automatic topology generation.

Return type:

System

Returns:

OpenMM System object from PDB file.