Restrained ensemble molecular dynamics method implementation for small angle scattering data. This is a custom force plugin for OpenMM simulation package. Currently we provide only CUDA implementation of the plugin, so you’ll need CUDA-enebled device to be able to run the simulation. Besides current installation of OpenMM you’ll need NumPy and mpi4py.

Building The Plugin

This project uses CMake for its build system. To build it, follow these steps:

1. Create a directory in which to build the plugin.

2. Run the CMake GUI or ccmake, specifying your new directory as the build directory and the top level directory of this project as the source directory.

3. Press "Configure".

4. Set OPENMM_DIR to point to the directory where OpenMM is installed. This is needed to locate the OpenMM header files and libraries.

5. Set CMAKE_INSTALL_PREFIX to the directory where the plugin should be installed. Usually, this will be the same as OPENMM_DIR, so the plugin will be added to your OpenMM installation.

6. If you plan to build the CUDA platform, make sure that CUDA_TOOLKIT_ROOT_DIR is set correctly and that RE_BUILD_CUDA_LIB is on.

7. Press "Configure" again if necessary, then press "Generate".

8. Use the build system you selected to build and install the plugin. For example, if you selected Unix Makefiles, type make install.

Python API

OpenMM uses SWIG to generate its Python API. SWIG takes an "interface file", which is essentially a C++ header file with some extra annotations added, as its input. It then generates a Python extension module exposing the C++ API in Python.

When building OpenMM's Python API, the interface file is generated automatically from the C++ API. That guarantees the C++ and Python APIs are always synchronized with each other and avoids the potential bugs that would come from have duplicate definitions. It takes a lot of complex processing to do that, though, and for a single plugin it's far simpler to just write the interface file by hand. You will find it in the "python" directory.

To build and install the Python API, build the "PythonInstall" target, for example by typing make PythonInstall. (If you are installing into the system Python, you may need to use sudo.) This runs SWIG to generate the C++ and Python files for the extension module (REPluginWrapper.cpp and REplugin.py), then runs a setup.py script to build and install the module. Once you do that, you can use the plugin from your Python scripts:

from simtk.openmm import System
from REplugin import REForce
system = System()
force = REForce()
system.addForce(force)

You might also need to change the variables in python/setup.py to your local path to OpenMM and the plugin directory

Basic usage

Once you have a working installation of the plugin you can start running the simulations. An example of the python script that runs restrained ensemble simulation is provided with this tutorial (run_RE.py). The script is a modification of a standard script for running an OpenMM simulation generated by CHARMM-GUI and it utilize mpi4py for running multiple copies of the system. The number of MPI ranks that you use for running this script correspond to the number of systems used in the restrained ensemble simulation.

We will go through the basic steps you need to do for setting up a restrained ensemble simulation of a POPS bilayer system. First, you need to have a lipid bilayer system with an associated topology for running a simulation with OpenMM. The easiest way to get it is to use CHARMM-GUI input generator. It will provide you with all the structures and topologies you need to run a classic MD simulation. Once you’ve preequilibrated your lipid bilayer you can start running restrained ensemble simulation. We will not go into the details of how you can run classic MD simulation with OpenMM, instead we will focus on the details specific for the restrained ensemble simulations.

The tutorial folder contain all the files you need for running a restrained ensemble simulation of POPS bilayer: X-ray (.xff) and neutron (.nff) experimental data (POPS_ULV@25Cin*.xff/.nff), force field definition (toppar folder and toppar.str file), system topology (system_pops.psf), system coordinate file (system_pops.crd), input parameters file (production_param.inp), restart file with the coordinates of the preequilibrated POPS bilayer (restart.rst), the script to run the simulation (run_RE.py) and the modules required for the OpenMM simulation (omm_*.py).

The run_RE.py file contain lots of comments on a particular functions we use to set up the simulation. Besides the system topology and input parameters we use a couple of files that contain information on which atoms are subject of the restraining force or which atoms are water/ions/exchangeable hydrogens (e.g. in_*.dat). These files contain indices of the appropriate group of atoms. Here is one of the ways you can use to generate these files. We will use system_pops.crd file as a source of data, but we have to remove remove first 4 lines of this file so we have only the lines that start with an atom number. We can do that with the following shell command:

tail -n +5 system_pops.crd > system_pops_cut.crd

That will write a new file (system_pops_cut.crd) without first four lines of system_pops.crd file. system_pops_cut.crd now contain all the atom indices as a first column, residue index, residue name as a third column, atom names as a fourth column and so on. First, we will generate a file with indices of the atoms that we want our restraining force to be able to act on. We want all the atoms in the system to be subject of the restraining force (if they are within the zcutoff distance from the bilayer center). So we just have to copy the first column of the system_pops_cut.crd to the in_atoms.dat file:

cat system_pops_cut.crd | awk '{print $1}' > in_atoms.dat

This command will feed the content of system_pops_cut.crd file to awk, which will print the first column of the input. Then we redirect the output to the in_atoms.dat file. Next, we want to get the file with indices of the atoms that we would use for the bilayer center of mass calculation. We would use all the POPS atoms for that:

cat system_pops_cut.crd | grep POPS | awk '{print $1}' > in_orig.dat

Again, we take the content of system_pops_cut.crd, keep only the lines that contain "POPS" in them (grep POPS) and then store the first column in the in_orig.dat file. Next, we generate the file with all the indices of the water:

cat system_pops_cut.crd | grep TIP3 | awk '{print $1}' > in_water.dat

And ions:

cat system_pops_cut.crd | grep POT | awk '{print $1}' >> in_ions.dat
cat system_pops_cut.crd | grep CLA | awk '{print $1}' >> in_ions.dat

POT and CLA are the atom names of potassium and chloride ions in .crd file. We also want all the atoms with the names HN1, HN2 and HN3 to be exchangeable hydrogens:

cat system_pops_cut.crd | grep HN1 | awk '{print $1}' >> in_exch_h.dat
cat system_pops_cut.crd | grep HN2 | awk '{print $1}' >> in_exch_h.dat
cat system_pops_cut.crd | grep HN3 | awk '{print $1}' >> in_exch_h.dat

We would need a file that we’ll use as a source of our atom names. We can use the system_pops_cut.crd file:

cat system_pops_cut.crd > for_names.dat

In .crd file we’ve had potassium and chloride atoms named as POT and CLA respectively, but when we determine scattering strength of an atom, we use either the first letter of the atom name (for non-ion or water atoms) or, for ions, the atom name should match one of the predefined ion names (which could be found in platforms/cuda/src/CudaFFMaps.cpp). For the potassium and chloride ions, their names should be ‘K+’ and ‘Cl-’ respectively. So, we should substitute POT with K+ and CLA with Cl- in for_names.dat file:

sed -i 's/'CLA'/'Cl-'/g' for_names.dat
sed -i 's/'POT'/'K+\ '/g' for_names.dat

Now, we have all the input files we need to run restrained ensemble simulation:

python -u run_RE.py -i production_param.inp -t toppar.str -p system_pops.psf -c system_pops.crd

You might also consider editing production_param.inp file to set classic MD-specific parameters. If we want to run multiple copies of the system, we should run the script as an MPI program. For example, if you want to run it on a cluster with SLURM scheduling system and save the output to the log.out file, you can do it with:

srun python -u run_RE.py -i production_param.inp -t toppar.str -p system_pops.psf -c system_pops.crd > log.out

srun can be substituted with mpirun or mpiexec commands if appropriate.

Once your simulation is finished, you can analyze the result. For example, you can have a look at the change in chi^2 and form factor. To do that, we first will filter the output files to get rid of the empty lines:

cat 0ff_xray.dat | grep -v '0 0 0 0' > 0ff_xray_tmp.dat
cat 0ff_xray_tmp.dat > 0ff_xray.dat
cat 0ff_neutron.dat | grep -v '0 0 0 0' > 0ff_neutron_tmp.dat
cat 0ff_neutron_tmp.dat > 0ff_neutron.dat
cat 0scale_xray.dat | sed -r '/^\s*$/d' > 0scale_xray_tmp.dat
cat 0scale_xray_tmp.dat > 0scale_xray.dat
cat 0scale_neutron.dat | sed -r '/^\s*$/d' > 0scale_neutron_tmp.dat
cat 0scale_neutron_tmp.dat > 0scale_neutron.dat
rm 0*_tmp.dat

Next, we can use get_chi.py script to calculate chi^2. It will write the chi^2 data and some other useful information:

python get_chi.py