7
$\begingroup$

I'm trying to use openpathsampling to simulate the binding of a ligand to my enzyme. I'm using Gromacs for my normal simulation on the HPC system provided by my university. I created an initial trajectory using some pull code for Gromacs. I tried to load this trajectory (in trr format) as described here:

from openpathsampling.engines.openmm.tools import trajectory_from_mdtraj

init_traj = md.load("md_pull_conf2_1st.trr", top="npt_conf2.gro")

trr = md.formats.TRRTrajectoryFile("md_pull_conf2_1st.trr")
vel = trr._read(n_frames=len(init_traj), atom_indices=None, get_velocities=True)[5]

OPS_init_traj = trajectory_from_mdtraj(init_traj, velocities=vel)

However, then I get the following error for the last step:

Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
  File "/[...]/openpathsampling/engines/openmm/tools.py", line 237, in trajectory_from_mdtraj
    
NameError: name 'MDTrajTopology' is not defined

From what I understand, in the tools.py script MDTrajTopology is only imported if OpenMM in installed. Since I can't install OpenMM on the HPC system, I tried to add the following import statement outside of if HAS_OPENMM:

from openpathsampling.engines import Topology, MDTrajTopology

If I run the import code like this, I get the following error:

Traceback (most recent call last):
  File "<stdin>", line 2, in <module>
  File "/[...]openpathsampling/engines/openmm/tools.py", line 244, in trajectory_from_mdtraj
    unit.nanometers)
  File "/[...]/simtk/unit/quantity.py", line 170, in __init__
    if value == None:
ValueError: The truth value of an array with more than one element is ambiguous. Use a.any() or a.all()

I also tried to load trajectory_from_mdtraj from openpathsampling.engines.openmm instead of openpathsampling.engines.openmm.tools. Then I get the following error:

Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
  File "/[...]/openpathsampling/engines/openmm/__init__.py", line 2, in missing_openmm
    raise RuntimeError("Install OpenMM to use this feature")
RuntimeError: Install OpenMM to use this feature

Is OpenMM required just to load the trajectory? Or is there a workaround which is not dependent on OpenMM?

$\endgroup$
0

2 Answers 2

6
$\begingroup$

Snapshot types in OpenPathSampling depend on the type of engine you're using. To some extent, snapshot types can be interchanged (if you had OpenMM available, I think the approach you're trying would work), but the OPS Gromacs engine has some different internals.

The Gromacs engine in OpenPathSampling mainly uses externally-defined trajectories (that is, trajectories that exist in external files, not in the internal format used by OPS). In this way, it is very different from the OpenMM engine, which stores everything in the internal format. Because the Gromacs engine only saves a reference to a file, if you need to recalculate something from an existing TRR file, it needs to be at the same relative path to the original when you use it later (example below).

Here's the overall process I would recommend to load a TRR file for use with Gromacs (code follows):

  1. Create the engine that you'll use for the dynamics.
  2. Use the engine's read_frame_from_file method to load the snapshots and wrap them in an OPS Trajectory object.

(Side note: In the following, I'll use list comprehensions (things that look like[f(x) for x in a_list]). If you (OP, or possibly a future reader) aren't familiar with list comprehensions, I recommend Trey Hunner's description.)

import openpathsampling as paths

# MDTraj only used in this script to get the number of frames
# (but OPS uses it internally, so you need it on your system)
import mdtraj as md

# change these filenames and the engine definition to meet your needs
trr_file = "example.trr"
pdb_file = "example.pdb"
engine = paths.engines.gromacs.Engine(
    gro=pdb_file,  # using pdb gives better MDTraj integration than gro
    mdp="md.mdp",
    top="topol.top",
    options={
        'snapshot_timestep': 0.02,
        'n_frames_max': 10000,
        'grompp_args': "-maxwarn 1",  # may be specific to my example.pdb
    }
).named("engine")
n_frames = len(md.load(trr_file, top=pdb_file))

# external_traj uses externally-stored snapshots
external_traj = paths.Trajectory([engine.read_frame_from_file(trr_file, num)
                                  for num in range(n_frames)])

From here, you can even save that external_traj and reuse it:

st = paths.Storage("setup.nc", mode='w')
st.save(external_traj)
st.close()

About keeping the original trajectory with the file: Now lets say I run the OPS command line from a subdirectory. From your shell (if you've installed the OpenPathSampling CLI) what happens if you run:

mkdir subdir && cd subdir
openpathsampling md ../setup.nc -e engine --init-frame -1 --nsteps 10 -o output.nc

(where --init-frame -1 will use the file's last frame as the initial frame.) Well, right now, you'll get a mess of an error that includes OSError: The file 'example.trr' doesn't exist. This is because OPS expects a file called example.trr in your current working directory. Fix it with:

cp ../example.trr ./
openpathsampling md ../setup.nc -e engine --init-frame -1 --nsteps 10 -o output.nc

And you'll see it runs Gromacs! Point of this digression: If you need the data from the TRR file, you'll need the TRR file. You'll also need to ensure it has the same path relative to your working directory. Be very careful not put a different file there with the same name, or you could get nonsense results!

(This isn't a concern for other Gromacs files, such as the options for mdp, gro, and top. Those files are actually stored with the OPS engine, and will be recreated if the directory is missing them.)

$\endgroup$
5
$\begingroup$

Note: not surprisingly, this answer doesn't resolve the main issue. See dwhswenson's answer for a full fix.

trajectory_from_mdtraj((mdtrajectory, simple_topology=False, velocities=None) is the full signature for the method you are calling. Within the method, there is this small segment, which controls whether MDTrajTopology is used:

if simple_topology:
    topology = Topology(*mdtrajectory.xyz[0].shape)
else:
    topology = MDTrajTopology(mdtrajectory.topology)

So specifying simple_topology=True should let you use OPS's native topology checker. I haven't tested this to ensure your example will work, but from looking at the code it seems like this would enough to skip around the OpenMM dependency.

$\endgroup$
5
  • $\begingroup$ Wow! Do you use OPS? I thought this would be answered by dwhswensen! $\endgroup$ Commented Oct 29, 2021 at 18:15
  • $\begingroup$ @NikeDattani Haven't ever used it, but the error message pointed to a specific spot in the code and the function was small/simple enough where I could at least pick this out. Not sure if it will fix everything for the OP, but it should at least fix the initial error message. $\endgroup$
    – Tyberius
    Commented Oct 29, 2021 at 18:18
  • 1
    $\begingroup$ This won't quite do it -- OPS treats snapshots differently for different MD engines. OpenMM requires explicit units, and this function (within engines.openmm) uses OpenMM's unit system, which won't be available on the OP's HPC machine. I've got an answer on the way. I was going to recommend a recent enhancement, but in testing out this specific use case I found a bug in that -- so I guess I'll only recommend the less-elegant-but-actually-working solution! $\endgroup$
    – dwhswenson
    Commented Oct 29, 2021 at 18:49
  • 1
    $\begingroup$ ... and just in case anyone ever stumbles across my previous comment and wonders whether the bug was fixed, here's the PR for it: github.com/openpathsampling/openpathsampling/pull/1086 $\endgroup$
    – dwhswenson
    Commented Oct 29, 2021 at 21:51
  • $\begingroup$ @dwhswenson at least the quick answer got this into the Hot Network Questions list, which is bringing in much more traffic than the average question here gets! $\endgroup$ Commented Oct 30, 2021 at 0:14

You must log in to answer this question.

Not the answer you're looking for? Browse other questions tagged .