LOOS  v2.3.2
area_profile.py
1 #!/usr/bin/env python
2 """
3 Compute the voronoi cross-sectional area for something (e.g. a protein) through the membrane.
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  if len(sys.argv) > 1 and sys.argv[1] == "--fullhelp":
17  print """
18 Usage:
19 area_profile.py system trajectory skip stride zmin zmax num_slices padding all-selection-string target-selection-string
20 
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
24  "stride"
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
28  (see below)
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
31 
32 Notes
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
38  is centered at x=y=0.
39 
40 Padding:
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.
53 
54 
55 
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)
62 
63 
64  """
65  sys.exit(0)
66  elif len(sys.argv) != 11 or sys.argv[1] == "-h" or sys.argv[1] == "--h":
67  print "Usage: area_profile.py system trajectory skip stride zmin zmax num_slices padding all-selection-string target-selection-string"
68  sys.exit(0)
69 
70  system_filename = sys.argv[1]
71  traj_filename = sys.argv[2]
72  skip = int(sys.argv[3])
73  stride = int(sys.argv[4])
74  zmin = float(sys.argv[5])
75  zmax = float(sys.argv[6])
76  num_slices = int(sys.argv[7])
77  padding = float(sys.argv[8])
78 
79  all_selection = sys.argv[9]
80  target_selection = sys.argv[10]
81 
82  print "# ", " ".join(sys.argv)
83 
84  system = loos.createSystem(system_filename)
85  traj = loos.createTrajectory(traj_filename, system)
86  pytraj = loos.PyTraj(traj, system, skip=skip, stride=stride)
87 
88  #traj.readFrame(skip)
89  #traj.updateGroupCoords(system)
90 
91  all_atoms = loos.selectAtoms(system, all_selection)
92  target_atoms = loos.selectAtoms(all_atoms, target_selection)
93 
94  slicers = []
95  slice_width = (zmax - zmin) / num_slices
96  for i in range(num_slices):
97  low = zmin + i*slice_width
98  high = zmin + (i+1)*slice_width
99  slicers.append(ZSliceSelector(low, high))
100 
101  areas = numpy.zeros([num_slices])
102  areas2 = numpy.zeros([num_slices])
103 
104  for snap in pytraj:
105  system.reimageByAtom()
106 
107  for i in range(num_slices):
108  all_slice_atoms = slicers[i](all_atoms)
109  target_slice_atoms = slicers[i](target_atoms)
110 
111  # run voronoi
112  v = VoronoiWrapper(all_slice_atoms, padding)
113  v.generate_voronoi()
114 
115  sr = SuperRegion()
116  sr.buildFromAtoms(target_slice_atoms, v)
117  a = sr.area()
118  areas[i] += a
119  areas2[i] += a*a
120 
121  # normalize
122  areas /= len(pytraj)
123  areas2 /= len(pytraj)
124  areas2 = numpy.sqrt(areas2 - areas*areas)
125  print "# ZSlice\tArea\tDev"
126  for i in range(num_slices):
127  z = zmin + (i+0.5)*slice_width
128  print z, areas[i], areas2[i]
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