HBond on every frames of a large trajectory
Hi, 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. Jérémie Knoops PhD Student Laboratory for Chemistry of Novel Materials<http://morris.umh.ac.be/default.aspx> University of Mons<http://portail.umons.ac.be/en2/> 20, Place du Parc B-7000 Mons (Belgium) Tel : +32(0)65 37 38 65<tel:%2B32%280%2965%2037%2038%2065>
Hi Jérémie, I suspect the problems may be due to the large data size rather than using the hbonds command in per-frame scripting. It is common to combine hbonds with per-frame scripting of trajectories... it’s even in a tutorial: <http://www.rbvi.ucsf.edu/chimera/docs/UsersGuide/tutorials/ensembles2.html#part1> However, if you had a specific example of what you mean by not working, and the problem is reproducible (say on a smaller dataset that you could easily share with us), you could use menu: Help… Report a Bug. I’ll have to leave it to others with no-GUI scripting experience to say whether there is a way to get these results without displaying the trajectory frames. Sorry about the difficulties, Elaine ----- Elaine C. Meng, Ph.D. UCSF Computer Graphics Lab (Chimera team) and Babbitt Lab Department of Pharmaceutical Chemistry University of California, San Francisco On Apr 22, 2016, at 9:30 AM, Jérémie KNOOPS [531802] <Jeremie.KNOOPS@umons.ac.be> wrote:
Hi,
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.
Jérémie Knoops PhD Student
On Apr 22, 2016, at 9:30 AM, Jérémie KNOOPS [531802] wrote:
Hi,
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.
Hi Jérémie, 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])[0] # load all the frames traj.loadAllFrames() # open the output file (change as desired) import os.path 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, angleSlop=recAngleSlop, cacheDA=True)) # close/flush the output file f.close() 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. --Eric Eric Pettersen UCSF Computer Graphics Lab
Hi Eric, Thank you very much, I didn't expect an (almost) ready to use script. I only change the loop, the one with enumerate didn't work. But now it works perfectly. This is the script I use : # 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])[0] # load all the frames traj.loadAllFrames() # open the output file (change as desired) import os.path f = open(os.path.expanduser('~/Documents/hbonds.txt'), "w") # get the functions and parameters we need for finding H-bonds from FindHBond import findHBonds, recDistSlop, recAngleSlop # loop through the coordinate sets (frames) c = traj.coordSets for frame in range(len(c)): traj.activeCoordSet = c[frame+1] # 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", frame+1, len(findHBonds([traj], distSlop=recDistSlop, angleSlop=recAngleSlop, cacheDA=True)) # close/flush the output file f.close() I have an additional question : what's the way to restrict the residues to search for hbond donors & acceptors with python commands? Jérémie ________________________________ De : Eric Pettersen [pett@cgl.ucsf.edu] Envoyé : samedi 23 avril 2016 2:29 À : Jérémie KNOOPS [531802] Cc : chimera-users@cgl.ucsf.edu Objet : Re: [Chimera-users] HBond on every frames of a large trajectory On Apr 22, 2016, at 9:30 AM, Jérémie KNOOPS [531802] wrote: Hi, 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. Hi Jérémie, 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])[0] # load all the frames traj.loadAllFrames() # open the output file (change as desired) import os.path 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, angleSlop=recAngleSlop, cacheDA=True)) # close/flush the output file f.close() 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. --Eric Eric Pettersen UCSF Computer Graphics Lab
Great! Yeah, I forgot that traj.coordSets is a dictionary (frame# -> coordinate set) rather than a list. Good work! I guess I should point out that your version of the script will only work if the frame numbers are continuous and start from 1, which should be the case for almost all trajectories. —Eric
On Apr 26, 2016, at 11:30 AM, Jérémie KNOOPS [531802] <Jeremie.KNOOPS@umons.ac.be> wrote:
Hi Eric,
Thank you very much, I didn't expect an (almost) ready to use script. I only change the loop, the one with enumerate didn't work. But now it works perfectly. This is the script I use :
# 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])[0]
# load all the frames traj.loadAllFrames()
# open the output file (change as desired) import os.path f = open(os.path.expanduser('~/Documents/hbonds.txt'), "w")
# get the functions and parameters we need for finding H-bonds from FindHBond import findHBonds, recDistSlop, recAngleSlop
# loop through the coordinate sets (frames) c = traj.coordSets for frame in range(len(c)): traj.activeCoordSet = c[frame+1] # 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", frame+1, len(findHBonds([traj], distSlop=recDistSlop, angleSlop=recAngleSlop, cacheDA=True))
# close/flush the output file f.close()
I have an additional question : what's the way to restrict the residues to search for hbond donors & acceptors with python commands?
Jérémie De : Eric Pettersen [pett@cgl.ucsf.edu <mailto:pett@cgl.ucsf.edu>] Envoyé : samedi 23 avril 2016 2:29 À : Jérémie KNOOPS [531802] Cc : chimera-users@cgl.ucsf.edu <mailto:chimera-users@cgl.ucsf.edu> Objet : Re: [Chimera-users] HBond on every frames of a large trajectory
On Apr 22, 2016, at 9:30 AM, Jérémie KNOOPS [531802] wrote:
Hi,
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.
Hi Jérémie, 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])[0]
# load all the frames traj.loadAllFrames()
# open the output file (change as desired) import os.path 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, angleSlop=recAngleSlop, cacheDA=True))
# close/flush the output file f.close()
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.
--Eric
Eric Pettersen UCSF Computer Graphics Lab
On Apr 26, 2016, at 11:30 AM, Jérémie KNOOPS [531802] <Jeremie.KNOOPS@umons.ac.be> wrote:
I have an additional question : what's the way to restrict the residues to search for hbond donors & acceptors with python commands?
I overlooked this question on my first read. You do it by specifying the donors/acceptors keywords in the findHBonds() call. The simplest way to do that is to select the residues you want to restrict the search to, by whatever means you find convenient, and then in the script: from chimera.selection import currentAtoms donors = acceptors = currentAtoms() …and the end of the FindHBonds call in the script would be modified to look like: …, cacheDA=True, donors=donors, acceptors=acceptors)) —Eric Eric Pettersen UCSF Computer Graphics Lab
participants (3)
-
Elaine Meng -
Eric Pettersen -
Jérémie KNOOPS [531802]