PLUMED-TUTORIALS by plumed-tutorials

Input files

All files and data required to perform this tutorial are on the GitHub repository. They can either be accessed by cloning the repository to your machine;

git clone https://github.com/blake-armstrong/binding-tutorial-PLUMED

or by following the download instructions in the README.md page on the repository to download the basic ZIP file.

Once downloaded, the input files required to run a simulation with OpenMM and PLUMED are provided in the simulation_files directory. Included in the scripts directory are bash scripts that will automate setting up the multiple-walker simulation files and will analyse the output for you. These scripts are not strictly necessary to perform this tutorial, as the actions they perform can be carried out manually by the user. All the data has already been produced and analysed, and is presented in the ‘data’ directory. We will discuss how that data was produced and analysed so you can do it all yourself.

OpennMM & PLUMED Installation

The first thing we need to do is get a working installation of OpenMM and openmmplumed. Comprehensive installation instructions for OpenMM. We recommend using conda, as it is the most straight forward way to install the software without needing to worry about dependencies. See instructions for installing conda. After testing the installation with:

python -m openmm.testInstallation

All platforms available to you will be shown. In this tutorial we use the CPU platform, and the input files provided set up a simulation to use the CPU platform, which every user should have available by default. The instructions to install the OpenMM PLUMED plugin are also very simple when using conda. This will also install the command line plumed tool which will be used for analysis later on.

Running a simulation

Running the bash script setup_walkers.sh

bash scripts/setup_walkers.sh

bash scripts/setup_walkers.sh NUMWALKERS

will create a set of directory trees within the data directory labelled with cylinder radii ranging from 0.0 nm to 0.5 nm. NUMWALKERS is an integer for the number of walkers desired with the default being 4 (NB this can be changed but must be greater than 1 for multiple walker metadynamics). These directories are where the multiple walker simulations will be run for systems with various radii for the flat-bottom cylindrical restraining potential. This setup script has been included so that the user can play around with different numbers of walkers in an easy way, if they so desire. By default, it will create directories and setups for four walkers. The details on how to run a single walker are discussed further later.

If we navigate to the directory for a radius of 0.0 nm,

cd data/radius_0.0

here we will find a run.sh script here that will handle submitting the run.py job for each walker. The run.py file (which is located in the simulation_files directory) will allow us to run our simulation using OpenMM through its Python API. The plumed.inp file which contains the PLUMED specifics. There is an ‘analysis.sh’ script here as well, which will handle the analysis of our simulation output, which is discussed here.

There is extensive documentation for using the OpenMM Python API which we direct you to if you cannot understand some component of the run.py file while we briefly discuss it here.

run.py

import openmm as mm
import openmm.app as app
import openmm.unit as unit
from openmmplumed import PlumedForce

These lines give us access to the OpenMM code, and the OpenMM PLUMED plugin that can parse a PLUMED input file, interpret it, and provide the instructions to OpenMM.

coordinates = "coord.pdb"           # Name of our coordinates file.    
forcefield = "argon.xml"            # Name of our force field file.
plumed_file = "plumed.inp"          # Name of our PLUMED input file.
temperature = 300 * unit.kelvin     # Simulation temperature.
timestep = 1e-3 * unit.picosecond   # Simulation timestep.
trel = 1 / unit.picosecond          # Thermostat relaxation constant for the Langevin thermostat.
thermo_file = "output.out"          # Name of the file to write our periodic simulation output to.
nthermo = 10000                     # How often to periodicly write simulation output (in steps).
traj_file = "trajectory.dcd"        # Name of the file to periodicly write our system coordinates to.
ntraj = 10000                       # How often to periodicly write our system coordinates (in steps).
nsteps = 2000000                    # Total number of steps to simulate for (this corresponds to 2 ns).

Simulation parameter definitions that will be used further in the file. When providing a unit with a dimension, we make use of the OpenMM unit object. By default, each walker will run for 2 ns which should take 5 - 10 minutes to complete. This should be enough to produce a PMF that is sufficiently converged to analyse for our purposes due to the simplicity of the model system.

pdb = app.PDBFile(coordinates)
forcefield = app.ForceField(forcefield)
platform = mm.Platform.getPlatformByName("CPU")
properties = {}

system = forcefield.createSystem(
    pdb.topology,
    nonbondedMethod=app.CutoffNonPeriodic,      # Run a non-periodic simulation.
    nonbondedCutoff=0.9 * unit.nanometer,
    useDispersionCorrection=False,
    constraints=None,
    rigidWater=False,
)

The details of this section are better explained in the OpenMM docs.

for n, force in enumerate(system.getForces()):
    if force.getName() == "CMMotionRemover":
        system.removeForce(n)
        break

This section is needed to prevent the thermal energy from being removed from the system as we have only one particle. It is not important to understand, as it would not generally be necessary for a real simulation.

with open(plumed_file, "r") as f:
    script = f.read()
system.addForce(PlumedForce(script))

This is where a PLUMED input script is read and our forces will be created and added to the system.

integrator = mm.LangevinMiddleIntegrator(temperature, trel, timestep)

simulation = app.Simulation(pdb.topology, system, integrator, platform, properties)

simulation.context.setPositions(pdb.positions)
simulation.context.setVelocitiesToTemperature(temperature, random.randrange(99999))

simulation.reporters.append(
    app.StateDataReporter(
        thermo_file,
        nthermo,
        separator="\t",
        step=False,
        time=True,
        potentialEnergy=True,
        kineticEnergy=False,
        totalEnergy=False,
        temperature=True,
        volume=True,
        density=True,
        progress=False,
        remainingTime=False,
        speed=True,
        elapsedTime=False,
    )
)

simulation.reporters.append(app.DCDReporter(traj_file, ntraj, enforcePeriodicBox=False))

simulation.reporters.append(
    app.StateDataReporter(
        stdout,
        int(nsteps / 50),
        separator="\t",
        step=False,
        time=False,
        potentialEnergy=True,
        kineticEnergy=False,
        totalEnergy=False,
        temperature=True,
        volume=True,
        density=False,
        progress=True,
        remainingTime=True,
        speed=True,
        elapsedTime=True,
        totalSteps=nsteps,
    )
)

simulation.step(nsteps)

The details of this section are better explained in the OpenMM docs.