LOOS  v2.3.2
Voronoi.py
1 #!/usr/bin/env python
2 # Do voronoi analysis on slices of a membrane
3 # 2D only
4 # This code is a wrapper around the scipy Voronoi class, which is itself
5 # a wrapper around QHull (http://www.qhull.org/)
6 #
7 # Alan Grossfield, University of Rochester Medical Center, Dept
8 # of Biochemistry and Biophysics
9 # Copyright 2014
10 #
11 
12 import sys
13 import loos
14 try:
15  import numpy
16 except ImportError as e:
17  print "Error importing numpy: ({0}): {1}".format(e.errno, e.strerror)
18 try:
19  from scipy.spatial import Voronoi
20 except ImportError as e:
21  print "Error importing Voronoi from scipy: ({0}): {1}".format(e.errno, e.strerror)
22 
24  def __init__(self, min_z, max_z):
25  self.min_z = float(min_z)
26  self.max_z = float(max_z)
27 
28  def __call__(self, atomicgroup):
29  ag = loos.AtomicGroup()
30  ag.periodicBox(atomicgroup.periodicBox())
31  for a in atomicgroup:
32  if self.min_z < a.coords().z() < self.max_z:
33  ag.append(a)
34  return ag
35 
36 
37 class VoronoiError(Exception):
38  """Base class for Voronoi package exceptions"""
39 
40  def __init__(self, msg):
41  self.msg = msg
42 
43  def __str__(self):
44  return(repr(self.msg))
45 
46 
48  """
49  Wrap the scipy Voronoi class, which in turn is a wrapper around QHull.
50  So, it's a wrapper wrapper, and if this program is called from a script,
51  then we'd have a wrapper wrapper wrapper. Now I just need m4 and Tod will
52  approve.
53 
54  atomicgroup is a LOOS AtomicGroup
55  pad is a float specifying how far out we will generate padding atoms
56  (fake "image" atoms used to emulate periodicity). 15 ang is a good
57  default if you're using all atoms or all heavy atoms, but you may
58  need to go farther if you're using a sparser selection (e.g. just
59  lipid phosphates)
60  """
61 
62  def __init__(self, atomicGroup, pad=15.0):
63  self.atoms = atomicGroup
64  self.edges = []
65  self.pad = pad
66  self.padding_atoms = []
67  self.voronoi = None
68  self.regions = []
69  self.superRegions = []
70  self.atoms_to_regions = {}
71 
72  def isPeriodic(self):
73  return self.atoms.isPeriodic()
74 
75  def num_atoms(self):
76  return len(self.atoms)
77 
78  def update_atoms(self, atomicGroup):
79  self.atoms = atomicGroup
80 
82  """
83  Build the list of atoms in the periodic images surrounding the central
84  image.
85 
86  This is necessary because QHull doesn't know anything about periodicity,
87  so we need to do fake it; otherwise, you'd have bizarre or infinite
88  areas for atoms near the edge of the box
89  """
90  if not self.isPeriodic():
91  raise VoronoiError("Periodic boundaries are required")
92 
93  box = self.atoms.periodicBox()
94  half_box = 0.5 * box
95  xbox = loos.GCoord(box.x(), 0.0, 0.0)
96  ybox = loos.GCoord(0.0, box.y(), 0.0)
97 
98  for a in self.atoms:
99  c = a.coords()
100  # generate the square neighbors
101  if (c.x() < -half_box.x() + self.pad):
102  self.padding_atoms.append(c+xbox)
103  if (c.x() > half_box.x() - self.pad):
104  self.padding_atoms.append(c-xbox)
105  if (c.y() < -half_box.y() + self.pad):
106  self.padding_atoms.append(c+ybox)
107  if (c.y() > half_box.y() - self.pad):
108  self.padding_atoms.append(c-ybox)
109 
110  # generate the diagonal neighbors
111  if ( (c.x() < -half_box.x() + self.pad) and
112  (c.y() < -half_box.y() + self.pad) ):
113  self.padding_atoms.append(c+xbox+ybox)
114 
115  if ( (c.x() > half_box.x() - self.pad) and
116  (c.y() > half_box.y() - self.pad) ):
117  self.padding_atoms.append(c-xbox-ybox)
118 
119  if ( (c.x() < -half_box.x() + self.pad) and
120  (c.y() > half_box.y() - self.pad) ):
121  self.padding_atoms.append(c+xbox-ybox)
122 
123  if ( (c.x() > half_box.x() - self.pad) and
124  (c.y() < -half_box.y() + self.pad) ):
125  self.padding_atoms.append(c-xbox+ybox)
126 
127  return len(self.padding_atoms)
128 
129  def num_padding_atoms(self):
130  return len(self.padding_atoms)
131 
132  def get_region_from_atomid(self, atomid):
133  return self.atoms_to_regions[atomid]
134 
135  def generate_voronoi(self):
136  # make the 2D list of points
137  points = []
138  for a in self.atoms:
139  points.append([a.coords().x(), a.coords().y()])
140 
141  numpad = self.generate_padding_atoms()
142  for p in self.padding_atoms:
143  points.append([p.x(), p.y()])
144 
145  points = numpy.array(points)
146  self.voronoi = Voronoi(points)
147 
148  # read in all of the ridges
149  for r in self.voronoi.ridge_vertices:
150  self.edges.append(Edge(r[0], r[1]))
151 
152  # now build the regions
153  # the scipy voronoi object has an array point_region, which
154  # maps the index of the input point to its associated Voronoi region
155  v = []
156  for vert in self.voronoi.vertices:
157  c = loos.GCoord(vert[0], vert[1], 0.0)
158  v.append(c)
159  for i in range(self.num_atoms()):
160  index = self.voronoi.point_region[i]
161  if index == -1:
162  raise ValueError, "point %d (atomId = %d) from voronoi decomposition isn't associated with a voronoi region; you may need to increase the padding value" % i, self.atoms[i].id()
163  r = self.voronoi.regions[index]
164  self.regions.append(Region(v, r, self.atoms[i]))
165  self.atoms_to_regions[self.atoms[i].id()] = self.regions[i]
166 
167 
168 
169 class Edge:
170  def __init__(self, index1, index2):
171  self.ind1 = index1
172  self.ind2 = index2
173 
174  def __eq__(self, other):
175  if ( (self.ind1 == other.ind1) and (self.ind2 == other.ind2) or
176  (self.ind2 == other.ind1) and (self.ind1 == other.ind2) ):
177  return True
178  else:
179  return False
180 
181 class Region:
182  def __init__(self, vert_array, indices, atom):
183  if len(indices) == 0:
184  raise ValueError, "can't have 0-length list of indices"
185  self.vertices = vert_array
186  self.indices = indices
187  self.atom = atom
188  self.edges = []
189 
190  # build list of indices
191  for i in range(len(indices)-1):
192  self.edges.append(Edge(indices[i], indices[i+1]))
193  self.edges.append(Edge(indices[-1], indices[0]))
194 
195  def num_indices(self):
196  return len(self.indices)
197 
198  def num_edges(self):
199  return len(self.edges)
200 
201  def area(self):
202  area = 0.0
203  for i in range(self.num_indices()):
204  if ( (i+1) == self.num_indices()):
205  j = 0
206  else:
207  j = i+1
208 
209  p1 = self.vertices[self.indices[i]]
210  p2 = self.vertices[self.indices[j]]
211  area += p1.x()*p2.y() - p2.x()*p1.y();
212 
213 
214  area = 0.5*abs(area)
215  return(area)
216 
217  def is_neighbor(self, other):
218  """
219  Two regions are neighbors if they have an edge in common
220  """
221  for e1 in self.edges:
222  for e2 in other.edges:
223  if (e1 == e2):
224  return(True)
225  return(False)
226 
227  def print_indices(self):
228  for i in range(self.num_indices()):
229  print self.indices[i], " ", self.vertices[self.indices[i]]
230 
231  def atomId(self):
232  return self.atom.id()
233 
234  def coords(self):
235  return self.atom.coords()
236 
237  def __str__(self):
238  string = ""
239  for i in range(self.num_indices()):
240  g = self.vertices[self.indices[i]]
241  string += "%f\t%f\n" % (g.x(), g.y())
242  g = self.vertices[self.indices[0]]
243  string += "%f\t%f\n" % (g.x(), g.y())
244  return string
245 
246 
248 
249  def __init__(self, regions=None):
250  if regions:
251  self.regions = regions
252  else:
253  self.regions = []
254 
255  def add_region(self, region):
256  self.regions.append(region)
257 
258  def buildFromAtoms(self, atomicGroup, voronoi):
259  for atom in atomicGroup:
260  self.add_region(voronoi.get_region_from_atomid(atom.id()))
261 
262  def area(self):
263  a = 0.0
264  for r in self.regions:
265  a += r.area()
266  return(a)
267 
268  def is_neighborRegion(self, other):
269  for r in self.regions:
270  if r.is_neighbor(other):
271  return True
272  return False
273 
274  def is_neighborSuperRegion(self, other):
275  for r in self.regions:
276  for r2 in other.regions:
277  if r.is_neighbor(r2):
278  return True
279  return False
280 
281  def print_indices(self):
282  for r in self.regions:
283  r.print_indices()
284  print
285 
286 if __name__ == '__main__':
287 
288  import sys
289 
290  structure = loos.createSystem("trj_1.pdb")
291  #structure = loos.createSystem("example.pdb")
292  #structure = loos.createSystem("b2ar.pdb")
293 
294 
295  #box = loos.GCoord(55., 77, 100)
296  #box = loos.GCoord(55., 77, 100)
297  phos = loos.selectAtoms(structure, 'name == "P"')
298  upper = loos.AtomicGroup()
299  for p in phos:
300  if p.coords().z() > 0:
301  upper.append(p)
302  upper.periodicBox(structure.periodicBox())
303  print upper.isPeriodic()
304  """
305  slice = loos.AtomicGroup()
306  for a in structure:
307  if a.coords().z() > 20 and a.coords().z() < 21:
308  slice.append(a)
309  slice.periodicBox(box)
310  """
311 
312 
313  v = VoronoiWrapper(upper)
314  #v = VoronoiWrapper(slice)
315  #print v.isPeriodic()
316  #print v.num_atoms()
317  v.generate_voronoi()
318  print v.num_atoms(), v.num_padding_atoms()
319  for i in range(v.num_atoms()):
320  print i, v.atoms[i].coords()
321  print v.regions[i].area()
322 
323  for i in range(v.num_padding_atoms()):
324  print i, v.padding_atoms[i]
325 
326  s = SuperRegion(v.regions[:1])
327  s2 = SuperRegion(v.regions[:2])
328  print s.area(), v.regions[0].area()
329  print s2.area(), v.regions[0].area()+v.regions[1].area()
330  s.add_region(v.regions[1])
331  print s.area()
332  s.print_indices()
333 
334  s3 = SuperRegion()
335  s3.buildFromAtoms(upper, v)
336  print s3.area()
337  total = 0.0
338  for r in v.regions:
339  total += r.area()
340  print r.coords()
341  print r
342  print total
343 
344 
def add_region(self, region)
Definition: Voronoi.py:255
def is_neighbor(self, other)
Definition: Voronoi.py:217
def num_indices(self)
Definition: Voronoi.py:195
def isPeriodic(self)
Definition: Voronoi.py:72
AtomicGroup selectAtoms(const AtomicGroup &source, const std::string selection)
Applies a string-based selection to an atomic group...
Definition: utils.cpp:195
def generate_padding_atoms(self)
Definition: Voronoi.py:81
Class for handling groups of Atoms (pAtoms, actually)
Definition: AtomicGroup.hpp:87
def num_atoms(self)
Definition: Voronoi.py:75
AtomicGroup createSystem(const std::string &filename, const std::string &filetype)
Definition: sfactories.cpp:119