In [1]:
from openforcefield.topology import Molecule, Topology
from openforcefield.typing.engines.smirnoff.forcefield import ForceField
from openforcefield.utils import get_data_file_path
from simtk import openmm, unit
import numpy as np
Define utility function we'll use to get energy of an OpenMM system
In [2]:
def get_energy(system, positions):
"""
Return the potential energy.
Parameters
----------
system : simtk.openmm.System
The system to check
positions : simtk.unit.Quantity of dimension (natoms,3) with units of length
The positions to use
Returns
---------
energy
"""
integrator = openmm.VerletIntegrator(1.0 * unit.femtoseconds)
context = openmm.Context(system, integrator)
context.setPositions(positions)
state = context.getState(getEnergy=True)
energy = state.getPotentialEnergy().in_units_of(unit.kilocalories_per_mole)
return energy
In this example, we load a single ethanol molecule with geometry information, parameterize it using the smirnoff99Frosst forcefield, and evaluate its energy. We then modify the parameter that is applied to the C-O-H angle and re-evaluate the energy.
Load a molecule
In [3]:
molecule = Molecule.from_file(get_data_file_path('molecules/ethanol.sdf'))
Get positions for use below
In [4]:
positions = molecule.conformers[0]
Load the smirnoff99Frosst forcefield file
In [5]:
ff = ForceField('openff-1.0.0.offxml')
Generate an openforcefield Topology containing only this molecule
In [6]:
topology = molecule.to_topology()
Parameterize the molecule, creating an OpenMM system
In [7]:
orig_system = ff.create_openmm_system(topology)
Calculate energy before parameter modification
In [8]:
orig_energy = get_energy(orig_system, positions)
print(f"Original energy: {orig_energy}")
Get parameters for the C-O-H angle
In [9]:
smirks = '[*:1]-[#8:2]-[*:3]' # SMIRKS for the parameter to retrieve
parameter = ff.get_parameter_handler('Angles').parameters[smirks]
Modify the parameters
In [10]:
parameter.k *= 0.9
parameter.angle *= 1.1
Evaluate energy after parameter modification
In [11]:
new_system = ff.create_openmm_system(topology)
new_energy = get_energy(new_system, positions)
Print out energy
In [12]:
print(f"Original energy: {orig_energy}. New energy: {new_energy}")
The SMIRNOFF spec aims to specify all aspects of a system that contribute to the energy within the ForceField object. This includes the treatment of long-range electrostatics and van der Waals interactions. This may be different for some users, as other packages set these parameters at runtime, for example in an AMBER mdin file or GROMACS MDP file.
This example evaluates the energy of a periodic box of solvent molecules using the "standard" smirnoff99Frosst settings (PME for electrostatics, 9 Angstrom cutoff for vdW interactions). It then changes the ForceField's vdW treatment method to "PME" and re-evaluates the energy.
In [13]:
# Create a new ForceField containing the smirnoff99Frosst parameter set:
forcefield = ForceField('openff-1.0.0.offxml')
# Inspect the long-range van der Waals method:
vdw_handler = forcefield.get_parameter_handler('vdW')
print(f"The vdW method is currently set to: {vdw_handler.method}")
Select a solvent box to parameterize.
In [14]:
from simtk.openmm import app
# A 239-molecule mixture of cyclohexane and ethanol
pdbfile = app.PDBFile(get_data_file_path('systems/packmol_boxes/cyclohexane_ethanol_0.4_0.6.pdb'))
# A 340-molecule mixture of propane, methane, and butanol.
#pdbfile = app.PDBFile(get_data_file_path('systems/packmol_boxes/propane_methane_butanol_0.2_0.3_0.5.pdb'))
# One cyclohexane in a box of roughly 1400 waters
#pdbfile = app.PDBFile(get_data_file_path('systems/packmol_boxes/cyclohexane_water.pdb'))
# One ethanol in a box of roughly 1300 waters
#pdbfile = app.PDBFile(get_data_file_path('systems/packmol_boxes/ethanol_water.pdb'))
Provide a "complete" (including bond orders, charges, and stereochemistry) representation of each molecule that might be in the PDB. This is necessary because a PDB representation of a molecule does not contain sufficient information for parameterization.
In [15]:
molecules = [Molecule.from_smiles('C'), # methane
Molecule.from_smiles('CCC'),# propane
Molecule.from_smiles('CCCCO'), # butanol
Molecule.from_smiles('O'), # water
Molecule.from_smiles('CCO'), #ethanol
Molecule.from_smiles('C1CCCCC1'), #cyclohexane
]
Create an Open Force Field Topology object by matching the Open Force Field molecules defined above to those in the PDB
In [16]:
top = Topology.from_openmm(pdbfile.topology, unique_molecules=molecules)
orig_system = forcefield.create_openmm_system(top)
orig_energy = get_energy(orig_system, pdbfile.getPositions())
print(f"Original energy: {orig_energy}")
Change the long-range van der Waals method to be PME
In [17]:
vdw_handler.method = 'PME'
print(f"The vdW method is currently set to: {vdw_handler.method}")
The Open Force Field toolkit applies vdW parameters using a SMIRKS-based typing scheme. Inspect the first few vdW parameters. These can be changed programmatically, as shown in example 1.
In [18]:
for vdw_param in forcefield.get_parameter_handler('vdW').parameters[0:3]:
print(vdw_param)
Now recompute the energy of the system using PME for long-range vdW interactions
In [19]:
new_system = forcefield.create_openmm_system(top)
new_energy = get_energy(new_system, pdbfile.getPositions())
print(f"Original energy (with LJ cutoff): {orig_energy}")
print(f"New energy (using LJ PME): {new_energy}")