3 Compute the voronoi cross-sectional area for something (e.g. a protein) through the membrane.
5 Alan Grossfield, University of Rochester Medical Center
__name__ == '__main__'
len(sys.argv) > 1 and
sys.argv == "--fullhelp"
19 area_profile.py system trajectory skip stride zmin zmax num_slices padding all-selection-string target-selection-string
21 system: system file (e.g. pdb, psf, parmtop, gro)
22 trajectory: trajectory file (periodic boundary information required)
23 skip, stride: how to move through the trajectory, skipping the first "skip" values and stepping by
25 zmin, zmax: range of z-values to consider for area calculation
26 num_slices: number of slices to break the vertical span into
27 padding: floating point number specifying how many extra layers of atoms are generated
29 all-selection-string: the set of atoms used to compute the voronoi decomposition
30 target-selection-string: set of atoms whose areas we report
33 1) all selections are forced to be subsets of the initial selection. This is
34 necessary for the mapping of areas to work correctly.
35 2) this program assumes that the system has already been centered such that
36 the z location of the membrane isn't drifting (z-slices are absolute, not
37 relative to the membrane center) and such that the periodic box
41 This program uses scipy Voronoi implementation, which is itself a wrapper around
42 Qhull. Qhull doesn't know anything about periodic boundary conditions as far as I
43 can tell, so we fake it by generating extra atoms so provide an environment around
44 the real ones; otherwise, you get insane areas for the atoms at the "edge" of the
45 box. You need to pick this value such that there are ~2-3 layers of atoms around
46 the real ones; in my experience, a value of 15 ang works well if you're selecting
47 all atoms or all heavy atoms. If your selection is more sparse (e.g. just phosphate
48 atoms), you will need to make it significantly larger, perhaps 30 ang). The best
49 test is to run with different pad values on a short trajectory and find the value
50 where the answer doesn't change. Since a larger pad value increases the number
51 of atoms used in the Voronoi decomposition, it does get more expensive, but it's
52 never too bad, and there's always a faster way to get a wrong answer.
56 Example selection choice:
57 '\!hydrogen' 'segid == "RHOD"'
58 use all heavy atoms for the voronoi decomposition, then pick out the rhodopsin
59 molecule to get its area
60 (as a rule, you should get virtually identical answers with all atoms and all
61 heavy atoms, but the latter will be dramatically faster)
len(sys.argv) != 11 or
sys.argv == "-h" or
sys.argv == "--h"
67 print "Usage: area_profile.py system trajectory skip stride zmin zmax num_slices padding all-selection-string target-selection-string"
system_filename = sys.argv
traj_filename = sys.argv
skip = int(sys.argv)
stride = int(sys.argv)
zmin = float(sys.argv)
zmax = float(sys.argv)
num_slices = int(sys.argv)
padding = float(sys.argv)
all_selection = sys.argv
target_selection = sys.argv
82 print "# "
, " "
pytraj = loos.PyTraj(traj, system, skip=skip, stride=stride)
slice_width = (zmax - zmin) / num_slices
low = zmin + i*slice_width
high = zmin + (i+1)*slice_width
areas = numpy.zeros([num_slices])
areas2 = numpy.zeros([num_slices])
all_slice_atoms = slicers[i](all_atoms)
target_slice_atoms = slicers[i](target_atoms)
v = VoronoiWrapper(all_slice_atoms, padding)
areas2 /= len(pytraj)
areas2 = numpy.sqrt(areas2 - areas*areas)
125 print "# ZSlice\tArea\tDev"
z = zmin + (i+0.5)*slice_width
z, areas[i], areas2[i]
pTraj createTrajectory(const std::string &filename, const AtomicGroup &g)
Factory function for reading in a trajectory file.
AtomicGroup selectAtoms(const AtomicGroup &source, const std::string selection)
Applies a string-based selection to an atomic group...
AtomicGroup createSystem(const std::string &filename, const std::string &filetype)