omm_remake_forces

remaking OpenMM forces with Custom*Forces

This repo demonstrates how to replace the 'native' forces in OpenMM - HarmonicBondForce, HarmonicAngleForce, PeriodicTorsionForce, and NonbondedForce - with Custom*Force versions. To run one of these, simply run python remake_bonds_force.py (or any of the other force types). This generates System objects with or without a Custom*Force and prints the energies to demonstrate they are equivalent. For example:

in:

python remake_bonds_force.py; python remake_angles_force.py; python remake_torsions_force.py; python remake_nonbonded_force.py

out: (see 🚨 below for why remaking the NonbondedForces compares a different total energy to the other systems)

Total force using OpenMM-derived forces:
	 -107896.91075643156 kJ/mol
Total force using custom forces:
	 -107896.91075643156 kJ/mol
Total force using OpenMM-derived forces:
	 -107896.91075643156 kJ/mol
Total force using custom forces:
	 -107896.91075643156 kJ/mol
Total force using OpenMM-derived forces:
	 -107896.91075643156 kJ/mol
Total force using custom forces:
	 -107896.91075643156 kJ/mol
Total force using OpenMM-derived forces:
	 -107055.8852928183 kJ/mol
Total force using custom forces:
	 -107055.88443283795 kJ/mol

The input file is a PDB of protein in solvent+ions from the share installed with OpenMM (typically run using simulatePdb.py).

Some of the code is hidden away in boilerplate.py. This is just for readability. The code is just to instantiate the System and Simulation objects and so it's relatively uninteresting:

def makeSystem(pdb):
    """
    Generate a `System` object
    Note:
        - I use CutoffPeriodic here - this uses Reaction Field for
          long range forces, not PME. Reason for this is that it's easier
          to write a CustomNonbondedForce for reaction field.
    """
    forcefield = ForceField('amber14-all.xml', 'amber14/tip3pfb.xml')
    system = forcefield.createSystem(
        pdb.topology,
        nonbondedMethod=CutoffPeriodic,
        nonbondedCutoff=0.9*nanometer
    )
    return system

def makeSimulation(system, pdb):
    """
    Generate an OpenMM `Simulation` object
    Note:
       - I use the `Reference` implementation here for simplicity.
         OpenCL or CUDA are preferred, particularly when replacing
         the nonbonded forces - CPU impl. is slow in this case.
    """
    integrator = openmm.LangevinIntegrator(
        298*kelvin,
        1/picosecond,
        1*femtosecond
    )
    platform = openmm.Platform.getPlatformByName('Reference')
    simulation = app.Simulation(pdb.topology, system, integrator, platform)
    simulation.context.setPositions(pdb.positions)
    return simulation

Notes:

• 🚨 remake_nonbonded_force.py currently turns OFF the long-range dispersion correction. For the reason why, see: openmm/openmm#3281 (comment) .
• All System objects are generted with app.CutoffPeriodic, which applies reaction field, rather than PME, for long-range electrostatics. For the reason why, see this comment: openmm/openmm#3281 (comment)
• Use of CustomNonbondedForce is slow on the CPU platform. For the reason why, see this comment: openmm/openmm#3162 (comment)