Example

We start by providing an example. Here is a python function script:

import openmm as mm
import openmm.app as app
import openmm.unit as unit

parameters = {
    "low_dG_and_low_barrier": (-7812.5, 7265.625, -2214.84375, 220.1171875, 4000.0, -4800.0, 1890.0, -245.),
    "low_dG_and_mid_barrier": (-15625, 14531.25, -4429.6875, 440.234375, 8000.0, -9600.0, 3780.0, -490.),
    "mid_dG_and_mid_barrier": (-5859.375, 5449.21875, -1661.1328125, 165.08789062, 8000.0, -9600.0, 3780.0, -490.),
    "low_dG_and_high_barrier": (-23437.5, 21796.875, -6644.53125, 660.3515625, 12000.0, -14400.0, 5670.0, -735.),
    "mid_dG_and_high_barrier": (-13671.875, 12714.84375, -3875.9765625, 385.20507812, 12000.0, -14400.0, 5670.0, -735.),
    "high_dG_and_high_barrier": (-3906.25, 3632.8125, -1107.421875, 110.05859375, 12000.0, -14400.0, 5670.0, -735.)
}

def poly(prms, xx):
    return (prms[0] * xx ** 3 + prms[1] * xx ** 2 + prms[2] * xx + prms[3]) * np.heaviside(0.35 - xx, 1) + (prms[4] * xx ** 3 + prms[5] * xx ** 2 + prms[6] * xx + prms[7]) * np.heaviside(xx - 0.35, 1)

def generate_efunc(prms):
    return f"4.184*(step(0.35-r)*({prms[0]}*r^3+{prms[1]}*r^2+{prms[2]}*r+{prms[3]})+step(r-0.35)*({prms[4]}*r^3+{prms[5]}*r^2+{prms[6]}*r+{prms[7]}))"

def sample(path):
    total_steps = 50000 # in ps

    syst = mm.System()
    for npair in range(100):
        syst.addParticle(12.0 * unit.amu)
        syst.addParticle(12.0 * unit.amu)
    force = mm.CustomBondForce(generate_efunc(parameters[path]))
    syst.addForce(force)
    for npair in range(100):
        force.addBond(2*npair+0, 2*npair+1, [])

    integrator = mm.LangevinIntegrator(300 * unit.kelvin, 1.0 / unit.picosecond, 0.002 * unit.picoseconds)
    topology = app.Topology()
    chain = topology.addChain()
    for npair in range(100):
        residue = topology.addResidue("MOL", chain)
        atom1 = topology.addAtom("C1", app.Element.getBySymbol('C'), residue)
        atom2 = topology.addAtom("C2", app.Element.getBySymbol('C'), residue)
    topology.addBond(atom1, atom2)

    simulation = app.Simulation(topology, syst, integrator)
    pos = []
    for npair in range(100):
        pos.append([0.0, 0.0, 0.0])
        pos.append([0.35 * np.random.random()*0.1, 0.0, 0.0])
    simulation.context.setPositions(np.array(pos))
    simulation.minimizeEnergy()
    simulation.context.setVelocitiesToTemperature(300 * unit.kelvin)
    simulation.reporters.append(app.StateDataReporter(f'{path}/output.txt', 100, step=True, potentialEnergy=True, temperature=True, speed=True, remainingTime=True, totalSteps=total_steps*500))
    simulation.reporters.append(app.DCDReporter(f'{path}/output.dcd', 100))
    for nstep in trange(total_steps):
        simulation.step(100)

Now we explain this function line by line, the first line:

syst = mm.System()

Creates an instance of the OpenMM System class, which represents the entire system to be simulated. Then the lines:

for npair in range(100):
    syst.addParticle(12.0 * unit.amu)
    syst.addParticle(12.0 * unit.amu)

Adds 200 particles (100 pairs) to the system, each with a mass of 12 atomic mass units (amu). Then the lines:

force = mm.CustomBondForce(generate_efunc(parameters[path]))
syst.addForce(force)
for npair in range(100):
    force.addBond(2*npair+0, 2*npair+1, [])

Creates a custom bond force using mm.CustomBondForce. Adds the force to the system and Defines bonds between pairs of particles. Then the line:

integrator = mm.LangevinIntegrator(300 * unit.kelvin, 1.0 / unit.picosecond, 0.002 * unit.picoseconds)

Creates a Langevin integrator, which integrates the equations of motion using Langevin dynamics: Temperature: 300 Kelvin, Friction coefficient: 1.0 / picoseconds, Time step: 0.002 picoseconds. Then the lines:

topology = app.Topology()
chain = topology.addChain()
for npair in range(100):
    residue = topology.addResidue("MOL", chain)
    atom1 = topology.addAtom("C1", app.Element.getBySymbol('C'), residue)
    atom2 = topology.addAtom("C2", app.Element.getBySymbol('C'), residue)
    topology.addBond(atom1, atom2)

Initializes the topology of the system, which represents the connectivity between atoms. Adds a single chain to the topology. For each particle pair, adds a residue, and two carbon atoms (C1 and C2), and defines a bond between them. Then the line:

simulation = app.Simulation(topology, syst, integrator)

Creates a Simulation object that ties together the system, topology, and integrator. Then the lines:

pos = []
for npair in range(100):
    pos.append([0.0, 0.0, 0.0])
    pos.append([0.35 * np.random.random()*0.1, 0.0, 0.0])
simulation.context.setPositions(np.array(pos))

Defines initial positions for the particles and sets these positions in the simulation context. Then the lines:

simulation.minimizeEnergy()
simulation.context.setVelocitiesToTemperature(300 * unit.kelvin)

Minimizes the energy of the system to remove any initial bad contacts, and Sets the initial velocities according to the specified temperature (300 Kelvin). Then the lines:

simulation.reporters.append(app.StateDataReporter(f'{path}/output.txt', 100, step=True, potentialEnergy=True, temperature=True, speed=True, remainingTime=True, totalSteps=total_steps*500))
simulation.reporters.append(app.DCDReporter(f'{path}/output.dcd', 100))

Adds reporters to record the state of the simulation. StateDataReporter: Writes various data (step, potential energy, temperature, speed, etc.) to a text file every 100 steps. DCDReporter: Writes the trajectory to a DCD file every 100 steps. Then the lines:

for nstep in trange(total_steps):
    simulation.step(100)

Runs the simulation for the specified total number of steps (total_steps), performing 100 steps in each iteration of the loop. Let’s discuss the form of the energy function, it can be written as:

\[\begin{equation} E(r) = \begin{cases} 4.184 \left( a_1 r^3 + a_2 r^2 + a_3 r + a_4 \right) & \text{for } r \leq 0.35 \\ 4.184 \left( b_1 r^3 + b_2 r^2 + b_3 r + b_4 \right) & \text{for } r > 0.35 \end{cases} \end{equation}\]

Different parameters define different values of a and b. (1kcal = 4.184kj)

Units

OpenMM uses the following units:

\[\begin{equation} \begin{array}{|l|l|} \hline \text { Quantity } & \text { Units } \\ \hline \text { distance } & \mathrm{nm} \\ \hline \text { time } & \mathrm{ps} \\ \hline \text { mass } & \text { atomic mass units } \\ \hline \text { charge } & \text { proton charge } \\ \hline \text { temperature } & \text { Kelvin } \\ \hline \text { angle } & \text { radians } \\ \hline \text { energy } & \mathrm{kJ} / \mathrm{mol} \\ \hline \end{array} \end{equation}\]

An atomic mass unit is equal to 1/12 the mass of a single atom of carbon-12.

Integrators

The langevin integrator is:

\[\begin{equation} m_i \frac{d \mathbf{v}\_i}{d t}=\mathbf{f}\_i-\gamma m_i \mathbf{v}\_i+\mathbf{R}\_i \end{equation}\]

In OpenMM the langevin integrator can be defined in the class LangevinIntegrator, where three parameters can be specified in the initialization (temperature, frictionCoeff and stepsize) like:

LangevinIntegrator(300 * unit.kelvin, 1.0 / unit.picosecond, 0.002 * unit.picoseconds)

OpenMM python package

Example

Units

Integrators

Enjoy Reading This Article?