3 Compute areas for different sets of atoms within a particular slice along the membrane
6 Alan Grossfield, University of Rochester Medical Center
__name__ == '__main__'
len(sys.argv) > 1 and
sys.argv == "--fullhelp"
21 program system trajectory skip zmin zmax padding selection-string1 [selection-string2 ...]
23 system: system file (e.g. pdb, psf, parmtop, gro)
24 trajectory: trajectory file (Periodic boundary information required)
25 zmin, zmax: floating point numbers used to select a particular slice of the system
26 padding: floating point number specifying how many extra layers of atoms are generated
28 selection-string1: the set of atoms used to compute the voronoi decomposition
29 selection-string2, etc: sets of atoms for which areas are reported
32 1) all selections are forced to be subsets of the initial selection. This is
33 necessary for the mapping of areas to work correctly.
34 2) this program assumes that the system has already been centered such that
35 the z location of the membrane isn't drifting (z-slices are absolute, not
36 relative to the membrane center) and such that the periodic box
40 This program uses scipy Voronoi implementation, which is itself a wrapper around
41 Qhull. Qhull doesn't know anything about periodic boundary conditions as far as I
42 can tell, so we fake it by generating extra atoms so provide an environment around
43 the real ones; otherwise, you get insane areas for the atoms at the "edge" of the
44 box. You need to pick this value such that there are ~2-3 layers of atoms around
45 the real ones; in my experience, a value of 15 ang works well if you're selecting
46 all atoms or all heavy atoms. If your selection is more sparse (e.g. just phosphate
47 atoms), you will need to make it significantly larger, perhaps 30 ang). The best
48 test is to run with different pad values on a short trajectory and find the value
49 where the answer doesn't change. Since a larger pad value increases the number
50 of atoms used in the Voronoi decomposition, it does get more expensive, but it's
51 never too bad, and there's always a faster way to get a wrong answer.
55 Example selection choices:
56 To map the areas for different membrane constituents, one could do something like
57 '\!hydrogen' 'resname == "PCGL"' 'resname == "PALM"' 'resname == "OLEO"' 'resname == "TIP3"'
58 The first selection is all heavy atoms, which we then break into contributions from
59 the headgroup, palmitate chain, oleate chain, and water (this is the old charmm27
60 format). Note that skipping hydrogens is a good idea -- it won't change the answer
61 noticeably, but it greatly reduces the number of atoms, which will make things
64 To track "lipid areas" in a bilayer with POPE and POPG lipids, one could do
65 'name == "P"' 'resname == "POPE"' 'resname == "POPG"'
len(sys.argv) < 8 or
sys.argv == "-h" or
sys.argv == "--h"
sys.argv, " system trajectory skip zmin zmax padding selection-string1 [selection-string2 ...]"
system_filename = sys.argv
traj_filename = sys.argv
skip = int(sys.argv)
zmin = float(sys.argv)
zmax = float(sys.argv)
padding = float(sys.argv)
selection_strings = sys.argv[7:]
86 print "# "
, " "
slicer = ZSliceSelector(zmin, zmax)
string += "\tArea"
106 print "# Frame"
slice_atoms = slicer(selections)
v = VoronoiWrapper(slice_atoms, padding)
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)