LOOS  v2.3.2
cylindrical-density.py
1 #!/usr/bin/env python
2 """
3  cylindrical-density.py : compute the density of atoms ("targets") around
4  a centering selection, output a 2D histogram in
5  r,z
6  The expected use-case for this is looking at lipid packing around a membrane
7  protein, when for whatever region you want to average out the "polar"
8  degrees of freedom.
9 
10  For example, a command line could look like:
11 
12  cylindrical-density.py sim.psf sim.dcd 'segid == "PROT"' 'segid =~ "PE" && \!hydrogen' -25 25 50 10 30 20
13 
14  This would read the system info from sim.psf, and use the trajectory sim.dcd.
15  The system would be translated such that "PROT" is at the origin, and then
16  the density of all heavy atoms with segment names containing "PE" would be
17  computed. The z-range of the histogram would go from -25:25 with 50 bins,
18  and r-range would be 10:30 with 20 bins.
19 
20  I highly suggest including the "&& \!hydrogen" part of the target selection,
21  since that will make the program run significantly faster without
22  substantially changing the information content (the slash in front of the "!"
23  may or may not be necessary, depending on which shell you use).
24 
25  Alan Grossfield
26 """
27 """
28 
29  This file is part of LOOS.
30 
31  LOOS (Lightweight Object-Oriented Structure library)
32  Copyright (c) 2013 Tod Romo, Grossfield Lab
33  Department of Biochemistry and Biophysics
34  School of Medicine & Dentistry, University of Rochester
35 
36  This package (LOOS) is free software: you can redistribute it and/or modify
37  it under the terms of the GNU General Public License as published by
38  the Free Software Foundation under version 3 of the License.
39 
40  This package is distributed in the hope that it will be useful,
41  but WITHOUT ANY WARRANTY; without even the implied warranty of
42  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
43  GNU General Public License for more details.
44 
45  You should have received a copy of the GNU General Public License
46  along with this program. If not, see <http://www.gnu.org/licenses/>.
47 """
48 
49 
50 import sys
51 import loos
52 import loos.pyloos
53 import numpy
54 import math
55 
56 header = " ".join(sys.argv)
57 print "# ", header
58 
59 system_file = sys.argv[1]
60 traj_file = sys.argv[2]
61 centering_selection_string = sys.argv[3]
62 target_selection_string = sys.argv[4]
63 zmin = float(sys.argv[5])
64 zmax = float(sys.argv[6])
65 znum_bins = int(sys.argv[7])
66 rmin = float(sys.argv[8])
67 rmax = float(sys.argv[9])
68 rnum_bins = int(sys.argv[10])
69 
70 system = loos.createSystem(system_file)
71 traj = loos.pyloos.Trajectory(traj_file, system)
72 
73 centering = loos.selectAtoms(system, centering_selection_string)
74 target = loos.selectAtoms(system, target_selection_string)
75 
76 
77 zbin_width = (zmax - zmin) / znum_bins
78 
79 rbin_width = (rmax - rmin) / rnum_bins
80 rmin2 = rmin*rmin
81 rmax2 = rmax*rmax
82 
83 hist = numpy.zeros( [rnum_bins, znum_bins])
84 
85 for frame in traj:
86 
87  centroid = centering.centroid()
88  target.translate(-centroid)
89  target.reimageByAtom()
90 
91  for atom in target:
92  x = atom.coords().x()
93  y = atom.coords().y()
94  z = atom.coords().z()
95 
96  r2 = x*x + y*y
97 
98  #print r2, rmin2, rmax2, z, zmin, zmax
99  if (zmin < z < zmax) and (rmin2 < r2 < rmax2):
100  #print "got here"
101  r = math.sqrt(r2)
102 
103  rbin = int((r - rmin)/ rbin_width)
104  zbin = int((z - zmin)/ zbin_width)
105 
106  hist[rbin, zbin] += 1.0
107 
108 hist /= len(traj)
109 
110 for i in range(rnum_bins):
111  rinner = rmin + i*rbin_width
112  router = rinner + rbin_width
113  rval = rmin + (i+0.5)*rbin_width
114 
115  norm = math.pi * (router*router - rinner*rinner)
116  for j in range(znum_bins):
117  zval = zmin + (j+0.5)*zbin_width
118  print rval, zval, hist[i,j] / norm
119  print
Python-based wrapper for LOOS Trajectories This class turns a loos Trajectory into something more pyt...
Definition: trajectories.py:55
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