[Chimera-users] HBond on every frames of a large trajectory
pett at cgl.ucsf.edu
Fri Apr 22 17:29:14 PDT 2016
On Apr 22, 2016, at 9:30 AM, Jérémie KNOOPS  wrote:
> What is the easiest and fastest way to get the number of H bonds on every frame of a large trajectory (>10 000 frames) and print it to a file?
> The hbonds command does not seem to work with the perframe command and displaying the whole trajectory to run a per-frame script is quite slow.
To do this efficiently, we need to use a Python script. In the script the important things to avoid are: 1) a redraw for every trajectory frame, and 2) creating and destroying graphics objects representing the H-bonds every frame. This means we are going to avoid using the handy-but-non-efficient runCommand() Python call, which typically causes a redraw, and we are going to call the functions underlying the FindHBonds tool, in order to find the h-bonded pairs without creating pseudobonds and so on. We will also step through the coordinate sets (frames) in Python.
Here's a script:
# the script assumes you've already opened the trajectory in Chimera
from chimera import openModels, Molecule
# we assume the trajectory is the only open model
traj = openModels.list(modelTypes=[Molecule])
# load all the frames
# open the output file (change as desired)
f = open(os.path.expanduser("~/hbonds.out"), "w")
# get the functions and parameters we need for finding H-bonds
from FindHBond import findHBonds, recDistSlop, recAngleSlop
# loop through the coordinate sets (frames)
for i, cs in enumerate(traj.coordSets):
traj.activeCoordSet = cs
# find and print the number of H-bonds, and cache which atoms
# are donors/acceptors, so they're not computed each frame
print>>f, "Frame", i+1, len(findHBonds([traj], distSlop=recDistSlop,
# close/flush the output file
Put the above in a file whose name ends in ".py" and open it in Chimera (after opening your trajectory) and it should write the number of hydrogens bonds for each frame to a file named "hbonds.out" in your home directory. I didn't test it but I'm fairly certain the logic is correct. I'm also pretty sure the syntax is correct, thought it's more possible I missed a parenthesis or quote somewhere. If you know Python it should be easy to fix any such errors . Otherwise, let me know the error and I'll provide a fix.
UCSF Computer Graphics Lab
-------------- next part --------------
An HTML attachment was scrubbed...
More information about the Chimera-users