How to Prepare Your Data

Preparing your biomolecular system for Wepy simulation mostly follows the same steps you would use to run a normal molecular dynamics (MD) trajectory with OpenMM.

To prepare your data for Wepy you can follow these 3 steps:

  1. Use a system builder (e.g. CharmmGUI) to build your system.

  2. Perform a short MD simulation with OpenMM.

  3. Use the simulation output to build a restart file, as well as pickle files for the system and topology.

We describe a protocol using CharmmGUI below, but this can be easily adapted to accomodate other workflows as well.

Building your system with CharmmGUI

CharmmGUI is a powerful web-based interface for preparing input files for MD simulations. There are a number of resources on that page to assist you in choosing the appropriate options for your project. A minimal workflow is provided below:

  1. Visit CharmmGUI: Open your web browser and go to https://www.charmm-gui.org . You will need to create an account and/or log in.

  2. Select the Appropriate Module: Choose the module that best suits your needs (e.g., “Membrane Builder”, “Solution Builder”, etc.). For most MD simulations, the “Solution Builder” is commonly used.

  3. Upload Your Molecular Structure: Upload your molecular structure file (e.g., PDB file) to the CharmmGUI interface, or enter the PDB ID. If uploading, ensure your file is correctly formatted and contains all necessary components.

  4. Set Up Your Simulation Parameters: Follow the guided steps in CharmmGUI to set up the details of your simulation. This typically includes selecting the appropriate chains and entities, defining the type of solvent, adding ions, setting the box size, temperature and other relevant details.

  5. Generate Input Files: Once you have configured your system, CharmmGUI will generate the necessary input files for your MD simulation. For output type, be sure to select OpenMM. This will include all of the necessary topology files, coordinate files, and parameter files.

  6. Download the Input Files: Download the tgz-format archive file to your local machine and/or move them into your cluster or computing environment. Expand this using:

    tar -xzf <charmm-gui-file>.tgz.

    The files are now in the directory charmm-gui-XXXXXXXXXX/openmm where XXXXXXXXXX is your particular job number. You will use the files in this directory to perform the MD simulation in the next step.

MD Simulation

The output from CharmmGUI is not yet equilibrated. It should be energy minimized and then gradually heated to the production temperature. CharmmGUI has provided scripts for this purpose, which are contained in the README file.

A nice way to carry out this equilibration is to make this file executable (with chmod +x README) and then execute its contents:

./README

This can also be run in a SLURM script, if working on a compute cluster:

#!/bin/bash
#SBATCH --job-name=charmm-gui-equil
#SBATCH --gpus=1
#SBATCH --time=4:00:00
#SBATCH --mem=5G

# activate the conda environment with openmm
conda activate wepy-env

# load a CUDA module if necessary
module load CUDA/11.8.0

cd ~/charmm-gui-XXXXXXXXXX/openmm

./README >& t.log

Building System and Topology Pickle Files

Now that you have your equilibrated system, you can use the simulation output to build system and topology pickle files. To achieve this you will use the openmm_run.py script in the charmm-gui-XXXXXXXXXX/openmm folder and the corresponding README file generated by CHARMM-GUI.

To build the system and topology pickle files, follow these steps:

  1. Modify the openmmrun.py script: This script is generated by CHARMM-GUI. It contains the necessary code to generate the system and topology objects which we will save in pkl files, and then exit the script. An example of this modified file is given below, which is also included in examples/charmm-gui/openmm_make_pkls.py.

    """
Generated by CHARMM-GUI (http://www.charmm-gui.org)

openmm_run.py

This program is OpenMM running scripts written in python.

Correspondance: jul316@lehigh.edu or wonpil@lehigh.edu
Last update: June 18, 2021
"""

from __future__ import print_function
import argparse
import sys
import os
import pickle as pkl

from omm_readinputs import *
from omm_readparams import *
from omm_vfswitch import *
from omm_barostat import *
from omm_restraints import *
from omm_rewrap import *

from openmm.unit import *
from openmm import *
from openmm.app import *

parser = argparse.ArgumentParser()
parser.add_argument('--platform', nargs=1, help='OpenMM platform (default: CUDA or OpenCL)')
parser.add_argument('-i', dest='inpfile', help='Input parameter file', required=True)
parser.add_argument('-p', dest='topfile', help='Input topology file', required=True)
parser.add_argument('-c', dest='crdfile', help='Input coordinate file', required=True)
parser.add_argument('-t', dest='toppar', help='Input CHARMM-GUI toppar stream file (optional)')
parser.add_argument('-b', dest='sysinfo', help='Input CHARMM-GUI sysinfo stream file (optional)')
parser.add_argument('-ff', dest='fftype', help='Input force field type (default: CHARMM)', default='CHARMM')
parser.add_argument('-icrst', metavar='RSTFILE', dest='icrst', help='Input CHARMM RST file (optional)')
parser.add_argument('-irst', metavar='RSTFILE', dest='irst', help='Input restart file (optional)')
parser.add_argument('-ichk', metavar='CHKFILE', dest='ichk', help='Input checkpoint file (optional)')
parser.add_argument('-opdb', metavar='PDBFILE', dest='opdb', help='Output PDB file (optional)')
parser.add_argument('-orst', metavar='RSTFILE', dest='orst', help='Output restart file (optional)')
parser.add_argument('-ochk', metavar='CHKFILE', dest='ochk', help='Output checkpoint file (optional)')
parser.add_argument('-odcd', metavar='DCDFILE', dest='odcd', help='Output trajectory file (optional)')
parser.add_argument('-rewrap', dest='rewrap', help='Re-wrap the coordinates in a molecular basis (optional)', action='store_true', default=False)
args = parser.parse_args()

# Load parameters
print("Loading parameters")
inputs = read_inputs(args.inpfile)

top = read_top(args.topfile, args.fftype.upper())
crd = read_crd(args.crdfile, args.fftype.upper())
if args.fftype.upper() == 'CHARMM':
params = read_params(args.toppar)
top = read_box(top, args.sysinfo) if args.sysinfo else gen_box(top, crd)

# Build system
nboptions = dict(nonbondedMethod=inputs.coulomb,
               nonbondedCutoff=inputs.r_off*nanometers,
               constraints=inputs.cons,
               ewaldErrorTolerance=inputs.ewald_Tol)
if inputs.vdw == 'Switch': nboptions['switchDistance'] = inputs.r_on*nanometers
if inputs.vdw == 'LJPME':  nboptions['nonbondedMethod'] = LJPME
if inputs.implicitSolvent:
nboptions['implicitSolvent'] = inputs.implicitSolvent
nboptions['implicitSolventSaltConc'] = inputs.implicit_salt*(moles/liter)
nboptions['temperature'] = inputs.temp*kelvin
nboptions['soluteDielectric'] = inputs.solut_diele
nboptions['solventDielectric'] = inputs.solve_diele
nboptions['gbsaModel'] = inputs.gbsamodel

if   args.fftype.upper() == 'CHARMM': system = top.createSystem(params, **nboptions)
elif args.fftype.upper() == 'AMBER':  system = top.createSystem(**nboptions)

if inputs.vdw == 'Force-switch': system = vfswitch(system, top, inputs)
if inputs.lj_lrc == 'yes':
for force in system.getForces():
      if isinstance(force, NonbondedForce): force.setUseDispersionCorrection(True)
         if isinstance(force, CustomNonbondedForce) and force.getNumTabulatedFunctions() != 1:
            force.setUseLongRangeCorrection(True)
if inputs.e14scale != 1.0:
for force in system.getForces():
      if isinstance(force, NonbondedForce): nonbonded = force; break
for i in range(nonbonded.getNumExceptions()):
         atom1, atom2, chg, sig, eps = nonbonded.getExceptionParameters(i)
         nonbonded.setExceptionParameters(i, atom1, atom2, chg*inputs.e14scale, sig, eps)

if inputs.pcouple == 'yes':      system = barostat(system, inputs)
if inputs.rest == 'yes':         system = restraints(system, crd, inputs)

# system and topology are made!
# write the pkl files
pkl.dump(system,open('system.pkl','wb'))
pkl.dump(top.topology,open('topology.pkl','wb'))

    NOTE again that this is just a truncated version of the openmm_run.py script with two pkl.dump commands at the end. It will create system and topology objects that correspond to the options specified in the CharmmGUI configuration files.

  2. Call the new script as done in the README: The README file, also generated by CHARMM-GUI, provides the call signature for running this script with the “production” parameters. Check your README script for the appropriate call signature. A sample command is as follows:

    python -u openmm_make_pkls.py -i step5_production.inp -t toppar.str -p step3_input.psf -c step3_input.crd

    This script will generate two pickle files: system.pkl and topology.pkl. These files will be used in the next step to run the simulation with the generated system and topology. Now you are all set to start setting up the Wepy simulation.

​** For rebinding simulations, you will need to generate system and topology files for the bonded state as well. **