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!
CustomTorsionForce: How to obtain the dihedral angle?
- Justin Chan
- Posts: 3
- Joined: Wed Jan 27, 2016 12:19 am
- Peter Eastman
- Posts: 2564
- Joined: Thu Aug 09, 2007 1:25 pm
Re: CustomTorsionForce: How to obtain the dihedral angle?
MDTraj has a routine to calculate dihedrals from atomic coordinates: http://mdtraj.org/development/api/gener ... _dihedrals
- Gideon Simpson
- Posts: 28
- Joined: Tue Sep 26, 2017 11:15 pm
Re: CustomTorsionForce: How to obtain the dihedral angle?
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
- Posts: 2564
- Joined: Thu Aug 09, 2007 1:25 pm
Re: CustomTorsionForce: How to obtain the dihedral angle?
The following code might have some errors in it, but it should give you the idea.
First, you'll need a MDTraj Topology object:
Now query the current positions from OpenMM and create a Trajectory from them.
Now you can use any of the MDTraj analysis functions.
First, you'll need a MDTraj Topology object:
Code: Select all
topology = mdtraj.from_openmm(simulation.topology)
Code: Select all
state = simulation.context.getState(getPositions=True)
positions = state.getPositions().value_in_unit(nanometers)
traj = mdtraj.Trajectory([positions], topology)
- Gideon Simpson
- Posts: 28
- Joined: Tue Sep 26, 2017 11:15 pm
Re: CustomTorsionForce: How to obtain the dihedral angle?
Yes, this works:
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.
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])
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
- Posts: 2564
- Joined: Thu Aug 09, 2007 1:25 pm
Re: CustomTorsionForce: How to obtain the dihedral angle?
You could do something like this:Is there a way to integrate quantities computed from the mdtraj in a custom integrator?
Code: Select all
for i in range(steps):
# ... compute the quantity with MDTraj
integrator.setGlobalVariable(...)
integrator.step(1)
Sure. Just follow the example at http://docs.openmm.org/latest/userguide ... other-data.Is it possible to set up a reporter which would call one of the mdtraj analysis functions?
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.Am I adding a lot of computational overhead with mdtraj by including these computations?