This article explains how to run hybrid simulations (PFP/MM) by using OpenMM with PFP. In these simulations, a specific region (such as a ligand or a particular substrate) is calculated using a universal machine learning potential (PFP), while the surrounding environment is calculated using a molecular mechanics (MM) force field.
We will use OpenMM to execute PFP/MM. For simulations using OpenMM, please first refer to the following article:
Additionally, AmberTools is useful for preparing classical force fields. Please see the article below for installing AmberTools on Matlantis:
Step 0. Required Input Files
To perform a simulation using PFP from OpenMM, the following input files are required. Please prepare them with appropriate modeling before starting the simulation.
-
topology.pdb: The full-system topology + coordinate file loaded by OpenMM. -
system.xml: The system file defining OpenMM's MM force field. An Amber format leap.parm7 or Gromacs format files can be used as an alternative. -
topology.xyz: The coordinate file for constructing Atoms in ASE. The atom ordering must match 1:1 perfectly withtopology.pdbandsystem.xml.
Step 1. Environment Setup and Verification
Step 1.1 Importing Required Libraries
Import OpenMM, ASE (Atomistic Simulation Environment), and the PFP API client required for the simulation.
import os
import sys
import numpy as np
# OpenMM-related imports
import openmm as mm
from openmm import app, unit
from openmm.app import PDBFile, Modeller, DCDReporter, StateDataReporter
from openmmml.mlpotential import MLPotential
# ASE and PFP-related imports
from ase.io import read
from pfp_api_client.pfp.estimator import Estimator, EstimatorCalcMode, EstimatorMethodType
from pfp_api_client.pfp.calculators.ase_calculator import ASECalculator
Step 1.2 Verifying Available Computing Platforms
Check the computing backends (Platforms) recognized by OpenMM. You will specify one of these verified backends to be used in the simulation later.
# Get and display available platform names
platform_names = [mm.Platform.getPlatform(i).getName() for i in range(mm.Platform.getNumPlatforms())]
print("Available Platform:")
for name in platform_names:
print(f" - {name}")
Step 2. System Construction and Configuration
Step 2.1 Initializing the PFP Calculator
Create a Calculator by specifying the PFP model version, method type, and calculation mode. Here, we use the R2SCAN_PLUS_D3 mode of version v8.0.0.
# PFP Estimator configuration
estimator = Estimator(
model_version="v8.0.0",
method_type=EstimatorMethodType.PFVM_D3_PFVM,
calc_mode=EstimatorCalcMode.R2SCAN_PLUS_D3
)
# Creating the ASE Calculator
pfp_calculator = ASECalculator(estimator)
Step 2.2 Setting Up the Topology and MM Force Field
Load the structure from the PDB file. Also, construct the system for MM force field calculations (mm_system) from OpenMM's XML-formatted topology and parameter file.
# 1. Load coordinates and topology
pdb_path = 'topology.pdb'
pdb = PDBFile(pdb_path)
modeller = Modeller(pdb.topology, pdb.positions)
# 2. Load the MM system (XML)
xml_path = 'system.xml'
with open(xml_path, 'r') as f:
xml_data = f.read()
mm_system = mm.XmlSerializer.deserialize(xml_data)
Those who are familiar with MM force field simulations might be more accustomed to Amber or GROMACS format topology and parameter files. In OpenMM, besides the OpenMM-specified XML format, you can also load MM force field topology and parameter information from formats like Amber or GROMACS to construct the mm_system.
For example, in the case of the Amber format, you can construct the mm_system as follows:
# Load Amber input files
prmtop = app.AmberPrmtopFile('leap.parm7') # Amber format parameter file
inpcrd = app.AmberInpcrdFile('leap.rst7') # Amber format coordinate/velocity file
# Prepare mm_system
mm_system = prmtop.createSystem(nonbondedMethod=app.PME, #app.CutoffNonPeriodic,
nonbondedCutoff = 1.0 * unit.nanometers,
constraints = app.HBonds,
rigidWater = True,
ewaldErrorTolerance=0.0005)
Step 2.3 Specifying the Region to Apply the Machine Learning Potential (PFP)
Determine which atoms among all atoms will be calculated with PFP. Here, we extract the indices of atoms whose residue name is "LIG" (generally the ligand).
Note: The PDB file loaded for the MM force field (
topology.pdb) and the XYZ file loaded for ASE (topology.xyz) must have the exact same atom ordering and match 1:1. If they do not match, PFP will not be applied to the correct region.
# 1. Load the initial structure (XYZ file) and set up the calculator
atoms = read('topology.xyz')
atoms.calc = pfp_calculator
potential = MLPotential('ase')
# 2. Extract atoms with residue name 'LIG' as the ML region
ml_atoms = [
atom.index for atom in modeller.topology.atoms()
if atom.residue.name == "LIG"
]
print(f"# of atoms in ML region: {len(ml_atoms)}")
print(f"Atoms in ML region: {ml_atoms}")
Step 2.4 Creating the Coupled PFP/MM System
Generate a mixed system (Mixed System) that incorporates the MLIP region (PFP) into the MM system using MLPotential.createMixedSystem.
# Constructing the coupled PFP/MM system
pfpmm_system = potential.createMixedSystem(
pdb.topology,
mm_system,
ml_atoms, # Atom indices of the ML region
calculator=pfp_calculator, # PFP Calculator
info={'charge': 0}, # Always set to 0 as PFP assumes charge neutrality
interpolate=False, # Set to False if you do not use a potential that mixes MM force field and PFP at a certain ratio
)
Step 2.5 Verifying the Force Field and ML Region
To confirm whether the MM force field and the machine learning potential are correctly configured, output the potential energy functions (Forces) set in the system and the details of the atoms specified in the ML region.
print("--- List of Applied Forces ---")
for force in pfpmm_system.getForces():
print(force)
force.setForceGroup(0)
print("\n--- Atom Details in the ML Region ---")
for iatom in ml_atoms:
atom = list(pdb.topology.atoms())[iatom]
print(f"Index: {iatom:4d} | Name: {atom.name:4s} | Element: {atom.element.symbol}")
Step 3. Defining Simulation Conditions
Step 3.1 Defining Output and Simulation Parameters
Create the output directory and define the conditions for the MD simulation (number of steps, saving intervals) as variables
# Create output directory
output_dir = 'output'
os.makedirs(output_dir, exist_ok=True)
# Simulation parameters
temperature = 300.0 * unit.kelvin # Temperature (Kelvin)
pressure = 1.0 * unit.bar # Pressure (bar)
dt = 1.0 * unit.femtoseconds # Timestep (fs)
fric_coeff = 1.0 / unit.picoseconds # Friction coefficient (ps^{-1})
md_steps = 200_000 # Total number of steps
dcd_period = 100 # Output interval for the trajectory file (DCD)
log_period = 100 # Output interval for the log file
Step 3.2 Initializing the Simulation Environment
Set up the platform to execute the calculation (CPU), the integrator for time evolution (LangevinMiddleIntegrator), and the simulation object.
Make sure to specify one of the available platforms verified in Step 1.2. If nothing is specified (leaving it to auto-selection), you can set platform=None or omit the platform argument when constructing the Simulation object.
# 1. Specify the platform
platform = mm.Platform.getPlatformByName('CPU')
# 2. Configure the integrator
integrator = mm.LangevinMiddleIntegrator(
temperature,
fric_coeff,
dt
)
integrator.setConstraintTolerance(0.00001)
# 3. Add a Monte Carlo Barostat
barostat = mm.MonteCarloBarostat(pressure, temperature)
pfpmm_system.addForce(barostat)
# 4. Construct the Simulation object
simulation = app.Simulation(modeller.topology, pfpmm_system, integrator, platform)
# 5. Set initial coordinates (Note: convert from ASE's Å units to OpenMM's nm units)
simulation.context.setPositions(atoms.get_positions() / 10.0)
Step 4. Energy Minimization
Before starting the simulation, perform energy minimization to remove unphysical atomic clashes or strains in the system. Here, we define a custom reporter to record and output the progress of the optimization.
from openmm import MinimizationReporter
class MinLogger(MinimizationReporter):
"""Custom reporter to monitor and record the progress of energy minimization"""
def __init__(self, log_path=None):
super().__init__()
self.history = []
self.f = open(log_path, 'w') if log_path else None
header = "iter\tsystem_E(kJ/mol)\trestraint_E\trestraint_k\tmax_constraint_err\t|g|max"
print(header)
if self.f:
self.f.write(header + "\n")
def report(self, iteration, x, grad, args):
gmax = float(np.max(np.abs(grad)))
line = (f"{iteration}\t{args['system energy']:.4f}\t"
f"{args['restraint energy']:.4f}\t"
f"{args['restraint strength']:.4f}\t"
f"{args['max constraint error']:.6f}\t"
f"|g|max={gmax:.4f}")
print(line, flush=True)
if self.f:
self.f.write(line + "\n")
self.f.flush()
self.history.append((iteration, args['system energy'], gmax))
return False # Returning True stops minimization process
# Initialize the reporter and execute energy minimization
minimize_log_path = os.path.join(output_dir, 'minimize.log')
reporter = MinLogger(minimize_log_path)
# Execute energy minimization
print("Start Energy Minimization")
simulation.minimizeEnergy(
tolerance=10 * unit.kilojoule_per_mole / unit.nanometer,
maxIterations=1000,
reporter=reporter,
)
reporter.f.close()
print("Complete Energy Minimization")
Step 5. Production MD Simulation
Step 5.1 Setting Up the Ensemble and Assigning Initial Velocities
Set the initial velocities for the system according to the specified temperature.
# Generate initial velocities based on the specified temperature
simulation.context.setVelocitiesToTemperature(temperature)
Step 5.2 Configuring Reporters
Configure the settings to write the trajectory file (DCD file) and thermodynamic quantities (energy, temperature, pressure, simulation speed, etc.) at each step into a log file and standard output (stdout).
# 1. Configuration to save the trajectory file (DCD)
dcd_path = os.path.join(output_dir, 'mdtraj.dcd')
simulation.reporters.append(app.DCDReporter(dcd_path, dcd_period))
# 2. Configuration for log file output
log_path = os.path.join(output_dir, 'mdtraj.log')
simulation.reporters.append(app.StateDataReporter(
log_path, log_period, step=True, potentialEnergy=True,
kineticEnergy=True, totalEnergy=True, temperature=True,
volume=True, progress=True, remainingTime=True,
speed=True, totalSteps=md_steps, separator='\t'
))
# 3. Configuration for standard output display
simulation.reporters.append(app.StateDataReporter(
sys.stdout, 100, step=True, potentialEnergy=True,
kineticEnergy=True, totalEnergy=True, temperature=True,
volume=True, progress=True, remainingTime=True,
speed=True, totalSteps=md_steps, separator='\t'
))
Step 5.3 Running the Simulation
Once everything is ready, execute the molecular dynamics simulation for the specified number of steps. Finally, save the simulation state (State) for restarting purposes.
print("Simulation Run")
# Run the simulation for the specified number of steps
simulation.step(md_steps)
# Save the final state
state_path = os.path.join(output_dir, 'mdtraj.xml')
simulation.saveState(state_path)
print("Simulation Done")
Limitations
Limitation 1: Handling of Covalent Bonds and LinkAtoms Across Regional Boundaries
If you want to specify a part of a protein's residues or a substrate bound to the main chain as the PFP region, you will cross "covalent bonds" at the boundary of the regions. If you try to run the current PFP/MM system by simply severing the bonds, you will face technical limitations.
Specifically, because the system is constructed by simply cutting the covalent bonds at the boundary between the PFP region and the MM region, dangling bonds (unbonded electrons) will occur in the region where PFP performs calculations. This makes the electronic state and local atomic environment recognized by PFP unnatural.
To address these issues, the LinkAtom method has been proposed as a way to cap severed covalent bonds with virtual atoms. However, the OpenMM code explained in this article (PFP/MM using MLPotential.createMixedSystem) does not support link atoms and cannot handle the "severing of covalent bonds across regions" and the "redistribution of forces associated with link atoms."
Therefore, when performing PFP/MM, it is strongly recommended to divide the PFP and MM regions using completely independent molecules that are not connected by covalent bonds.
Limitation 2: Mechanical Embedding
In the current basic PFP/MM implementation using OpenMM, "Mechanical Embedding" is adopted as the coupling scheme for interactions between PFP and MM regions.
In mechanical embedding, the PFP region is evaluated in a virtually isolated state (in vacuo), and the hybrid potential is constructed by adding up the interactions occurring with atoms in the MM region using only non-bonded interactions at the MM force field level (electrostatic Coulomb interactions and Lennard-Jones potential).
The limitation of this approach is that it cannot directly account for the electronic polarization effects exerted on the PFP region by the surrounding MM environment (such as proteins or polar solvents). Therefore, in systems where chemical reactions or charge transfer reactions involve severe charge redistribution, or where the environmental polarity or hydrogen-bonding network plays a decisive role in stabilizing the transition state, there is a risk of significant divergence from experimental values or all-atom QM calculations.
As a countermeasure to this limitation, a method is used where not only the reacting solute (substrate) but also several molecules of the nearby surrounding solvent are explicitly included in the PFP region for calculation. Since PFP calculation is highly scalable, it is easible to set a relatively large PFP region.