LOOS  v2.3.2
area_per_molecule.py
1 #!/usr/bin/env python
2 """
3 Compute distribution of areas/molecule for a z-slice
4 
5 Alan Grossfield, University of Rochester Medical Center
6 Copyright 2014
7 """
8 
9 import sys
10 import loos
11 import numpy
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 [selection-string3 ...]
22 
23 system: system file (e.g. psf) -- MUST CONTAIN CONNECTIVITY INFORMATION
24 trajectory: trajectory file (must have periodic boundary information)
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 
41  This program uses scipy Voronoi implementation, which is itself a wrapper
42  around Qhull. Qhull doesn't know anything about periodic boundary
43  conditions as far as I can tell, so we fake it by generating extra atoms so
44  provide an environment around the real ones; otherwise, you get insane
45  areas for the atoms at the "edge" of the box. You need to pick this value
46  such that there are ~2-3 layers of atoms around the real ones; in my
47  experience, a value of 15 ang works well if you're selecting all atoms or
48  all heavy atoms. If your selection is more sparse (e.g. just phosphate
49  atoms), you will need to make it significantly larger, perhaps 30 ang).
50  The best test is to run with different pad values on a short trajectory and
51  find the value where the answer doesn't change. Since a larger pad value
52  increases the number of atoms used in the Voronoi decomposition, it does
53  get more expensive, but it's never too bad, and there's always a faster way
54  to get a wrong answer.
55 
56 
57 
58 Example selection choices:
59  To compute the distribution of areas/lipid molecule, where there are
60  2 kinds of lipids, one with segids L1, L2, ... L120, and the other with
61  segids L121, L121, ... you could do 3 selections like
62 
63  'segid =~ "^L\d+" && \!hydrogen' 'segid -> "^L(\d+)" <= 120' 'segid -> "^L(\d+)" > 120'
64 
65 Note that for a bilayer you have to think carefully about the z-range you use.
66 If you use a z-range of [0, 20], you'll get the upper leaflet accurately, but
67 you'll also pick up some junk at very lower area due to the other leaflet
68 "leaking" across
69 
70  """
71  sys.exit(0)
72  elif len(sys.argv) < 11 or sys.argv[1] == "-h" or sys.argv[1] == "--h":
73  print sys.argv[0], " system trajectory skip zmin zmax padding selection-string1 selection-string2 [selection-string3 ...]"
74  sys.exit(0)
75 
76 
77  system_filename = sys.argv[1]
78  traj_filename = sys.argv[2]
79  skip = int(sys.argv[3])
80  zmin = float(sys.argv[4])
81  zmax = float(sys.argv[5])
82  padding = float(sys.argv[6])
83  min_area = float(sys.argv[7])
84  max_area = float(sys.argv[8])
85  num_bins = int(sys.argv[9])
86  selection_strings = sys.argv[10:] # first selection gives the all atoms to use
87  # in area calculations, all others tell you
88  # how to group the areas
89  bin_width = (max_area - min_area) / num_bins
90 
91  histograms = numpy.zeros([len(selection_strings[1:]), num_bins], numpy.float)
92 
93  print "# ", " ".join(sys.argv)
94 
95  system = loos.createSystem(system_filename)
96  traj = loos.createTrajectory(traj_filename, system)
97 
98  traj.readFrame(skip)
99  traj.updateGroupCoords(system)
100 
101  slicer = ZSliceSelector(zmin, zmax)
102 
103  selections = []
104  selections.append(loos.selectAtoms(system, selection_strings[0]))
105  for s in selection_strings[1:]:
106  selections.append(loos.selectAtoms(selections[0], s))
107  for i in range(1,len(selections)):
108  selections[i] = selections[i].splitByMolecule()
109 
110 
111  frame = skip
112  string = ""
113  for i in range(len(selections)):
114  string += "\tArea" + str(i)
115 
116  print "# Frame", string
117  while (traj.readFrame()):
118  traj.updateGroupCoords(system)
119  system.reimageByAtom()
120 
121  slice_atoms = slicer(selections[0])
122 
123  # run voronoi
124  v = VoronoiWrapper(slice_atoms, padding)
125  v.generate_voronoi()
126 
127  # generate the areas for the selections
128  for i in range(len(selections[1:])):
129  s = selections[i+1]
130  for j in range(len(s)):
131  sr = SuperRegion()
132  gr = slicer(s[j])
133  if (len(gr) == 0): # skip if there are no atoms selected from mol
134  continue
135  sr.buildFromAtoms(gr, v)
136  a = sr.area()
137  index = int((a-min_area) / bin_width)
138  try:
139  histograms[i][index] += 1
140  except IndexError:
141  print "# ", i, j, a, index
142 
143  # normalize the histograms
144  for i in range(len(histograms)):
145  histograms[i] /= numpy.add.reduce(histograms[i])
146 
147 
148  # final output
149  string = ""
150  for i in range(len(selections[1:])):
151  string += "\tSel" + str(i)
152  print "# Area", string
153  for i in range(num_bins):
154  a = min_area + (i+0.5)*bin_width
155  string = map(str, histograms[:,i])
156  print a, "\t".join(string)
157 
158 
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