[Chimera-users] HBond on every frames of a large trajectory

[Eric Pettersen]

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 at 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 at cgl.ucsf.edu <mailto:pett at cgl.ucsf.edu>]
> Envoyé : samedi 23 avril 2016 2:29
> À : Jérémie KNOOPS [531802]
> Cc : chimera-users at cgl.ucsf.edu <mailto:chimera-users at 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

-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://plato.cgl.ucsf.edu/pipermail/chimera-users/attachments/20160426/29905214/attachment.html>