CustomTorsionForce: How to obtain the dihedral angle?

[Justin Chan]

I am using the CustomTorsionForce with the potential: "0.5*k*(1-cos(theta-theta0))" where theta is the current dihedral angle and theta0 is the equilibrium dihedral angle.

I would like to compute the force of this potential with respect to theta therefore I need the current theta. Is there a way to obtain the current theta given the current state of the system? Or do I have to compute it from the coordinates?

Thanks!

[Peter Eastman]

MDTraj has a routine to calculate dihedrals from atomic coordinates: http://mdtraj.org/development/api/gener ... _dihedrals

[Gideon Simpson]

I had a similar question to this regarding the various functions that are a part of MDTraj, including dihedral angles. While it's clear to me how to calculate these quantities in a post-processing step, by dumping to disk and then loading in as an mdtraj Trajectory object, it's not at all clear how to do things online, such as the original question about computing the dihedral angle for a force computation. Is there a way to put a wrapper around the OpenMM data structures (i.e., Simulation), to be able to call the MDTraj functions on the current state of the system?

[Peter Eastman]

The following code might have some errors in it, but it should give you the idea.

First, you'll need a MDTraj Topology object:

Code: Select all

topology = mdtraj.from_openmm(simulation.topology)

Now query the current positions from OpenMM and create a Trajectory from them.

Code: Select all

state = simulation.context.getState(getPositions=True)
positions = state.getPositions().value_in_unit(nanometers)
traj = mdtraj.Trajectory([positions], topology)

Now you can use any of the MDTraj analysis functions.

[Gideon Simpson]

Yes, this works:

Code: Select all

md_topology = mdtraj.Topology.from_openmm(simulation.topology)
positions = state.getPositions().value_in_unit(unit.meter)
traj = mdtraj.Trajectory([positions], md_topology)
print('center of mass', mdtraj.compute_center_of_mass(traj)[0])

A few follow up questions:
1. Is there a way to integrate quantities computed from the mdtraj in a custom integrator? That is to say, perhaps I'd have some additional variable I had hadded, either global or per dof, and it was based on some quantity that mdtraj had a function for. Is it possible to set up (within python) the custom integrator to update this value at each step?

2. Is it possible to set up a reporter which would call one of the mdtraj analysis functions?

3. Am I adding a lot of computational overhead with mdtraj by including these computations? Right now, I am just playing with small clusters of atoms, so maybe there's a penalty I'm not yet seeing.

[Peter Eastman]

Is there a way to integrate quantities computed from the mdtraj in a custom integrator?

You could do something like this:

Code: Select all

for i in range(steps):
  # ... compute the quantity with MDTraj
  integrator.setGlobalVariable(...)
  integrator.step(1)

but that could easily end up being really slow, especially if you're using a GPU.

Is it possible to set up a reporter which would call one of the mdtraj analysis functions?

Sure. Just follow the example at http://docs.openmm.org/latest/userguide ... other-data.

Am I adding a lot of computational overhead with mdtraj by including these computations?

That depends how expensive the computation is, and how often you do it. Also, if you're running on a GPU then any interruption to bring data back from the GPU will have some cost.