LOOS  v2.3.2
run_areas.py
1 #!/usr/bin/env python
2 """
3 Compute areas for different sets of atoms within a particular slice along the membrane
4 normal.
5 
6 Alan Grossfield, University of Rochester Medical Center
7 Copyright 2014
8 """
9 
10 import sys
11 import loos
12 from Voronoi import *
13 
14 if __name__ == '__main__':
15 
16 
17  if len(sys.argv) > 1 and sys.argv[1] == "--fullhelp":
18  print """
19 
20 Usage:
21 program system trajectory skip zmin zmax padding selection-string1 [selection-string2 ...]
22 
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
27  (see below)
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
30 
31 Notes
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
37  is centered at x=y=0.
38 
39 Padding:
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.
52 
53 
54 
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
62  run much faster.
63 
64  To track "lipid areas" in a bilayer with POPE and POPG lipids, one could do
65  'name == "P"' 'resname == "POPE"' 'resname == "POPG"'
66 
67 
68 
69  """
70  sys.exit(0)
71  elif len(sys.argv) < 8 or sys.argv[1] == "-h" or sys.argv[1] == "--h":
72  print sys.argv[0], " system trajectory skip zmin zmax padding selection-string1 [selection-string2 ...]"
73  sys.exit(0)
74 
75 
76  system_filename = sys.argv[1]
77  traj_filename = sys.argv[2]
78  skip = int(sys.argv[3])
79  zmin = float(sys.argv[4])
80  zmax = float(sys.argv[5])
81  padding = float(sys.argv[6])
82  selection_strings = sys.argv[7:] # first selection gives the all atoms to use
83  # in area calculations, all others tell you
84  # how to group the areas
85 
86  print "# ", " ".join(sys.argv)
87 
88  system = loos.createSystem(system_filename)
89  traj = loos.createTrajectory(traj_filename, system)
90 
91  traj.readFrame(skip)
92  traj.updateGroupCoords(system)
93 
94  slicer = ZSliceSelector(zmin, zmax)
95 
96  selections = []
97  selections.append(loos.selectAtoms(system, selection_strings[0]))
98  for s in selection_strings[1:]:
99  selections.append(loos.selectAtoms(selections[0], s))
100 
101  frame = skip
102  string = ""
103  for i in range(len(selections)):
104  string += "\tArea" + str(i)
105 
106  print "# Frame", string
107  while (traj.readFrame()):
108  traj.updateGroupCoords(system)
109  system.reimageByAtom()
110 
111  slice_atoms = slicer(selections[0])
112 
113  # run voronoi
114  v = VoronoiWrapper(slice_atoms, padding)
115  v.generate_voronoi()
116 
117  areas = []
118  # generate the areas for the selections
119  for s in selections: # really should be selections[1:]
120  sr = SuperRegion()
121  sr.buildFromAtoms(slicer(s), v)
122  areas.append(sr.area())
123  print frame, "\t".join(map(str,areas))
124  frame += 1
pTraj createTrajectory(const std::string &filename, const AtomicGroup &g)
Factory function for reading in a trajectory file.
Definition: sfactories.cpp:184
AtomicGroup selectAtoms(const AtomicGroup &source, const std::string selection)
Applies a string-based selection to an atomic group...
Definition: utils.cpp:195
AtomicGroup createSystem(const std::string &filename, const std::string &filetype)
Definition: sfactories.cpp:119