LOOS  v2.3.2
cylindrical-thickness.py
1 #!/usr/bin/env python
2 """
3  cylindrical-thickness.py : compute the thickness of the membrane around
4  a centering selection, output the average as
5  a function of lateral distance from origin
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-thickness.py sim.psf sim.dcd 'segid == "PROT"' 'segid =~ "PE" && name == "P"' 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 heights of phosphate with segment names containing "PE" would be
17  computed.
18 
19  NOTE: This code assumes the membrane is centered at z = 0. Atoms with
20  z>0 are assigned to the upper leaflet, z<0 the lower leaflet, and the
21  difference between the averages is computed as a function of r.
22 
23 
24  Alan Grossfield
25 """
26 """
27 
28  This file is part of LOOS.
29 
30  LOOS (Lightweight Object-Oriented Structure library)
31  Copyright (c) 2013 Tod Romo, Grossfield Lab
32  Department of Biochemistry and Biophysics
33  School of Medicine & Dentistry, University of Rochester
34 
35  This package (LOOS) is free software: you can redistribute it and/or modify
36  it under the terms of the GNU General Public License as published by
37  the Free Software Foundation under version 3 of the License.
38 
39  This package is distributed in the hope that it will be useful,
40  but WITHOUT ANY WARRANTY; without even the implied warranty of
41  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
42  GNU General Public License for more details.
43 
44  You should have received a copy of the GNU General Public License
45  along with this program. If not, see <http://www.gnu.org/licenses/>.
46 """
47 
48 
49 import sys
50 import loos
51 import loos.pyloos
52 import numpy
53 import math
54 
55 header = " ".join(sys.argv)
56 print "# ", header
57 
58 system_file = sys.argv[1]
59 traj_file = sys.argv[2]
60 centering_selection_string = sys.argv[3]
61 target_selection_string = sys.argv[4]
62 rmin = float(sys.argv[5])
63 rmax = float(sys.argv[6])
64 rnum_bins = int(sys.argv[7])
65 
66 system = loos.createSystem(system_file)
67 traj = loos.pyloos.Trajectory(traj_file, system)
68 
69 centering = loos.selectAtoms(system, centering_selection_string)
70 target = loos.selectAtoms(system, target_selection_string)
71 
72 
73 rbin_width = (rmax - rmin) / rnum_bins
74 rmin2 = rmin*rmin
75 rmax2 = rmax*rmax
76 
77 upper_sum = numpy.zeros(rnum_bins)
78 upper_count = numpy.zeros(rnum_bins)
79 lower_sum = numpy.zeros(rnum_bins)
80 lower_count = numpy.zeros(rnum_bins)
81 
82 for frame in traj:
83 
84  centroid = centering.centroid()
85  target.translate(-centroid)
86  target.reimageByAtom()
87 
88  for atom in target:
89  x = atom.coords().x()
90  y = atom.coords().y()
91  z = atom.coords().z()
92 
93  r2 = x*x + y*y
94 
95  if (rmin2 < r2 < rmax2):
96  r = math.sqrt(r2)
97 
98  rbin = int((r - rmin)/ rbin_width)
99  if z > 0:
100  upper_sum[rbin] += z
101  upper_count[rbin] += 1
102  else:
103  lower_sum[rbin] += z
104  lower_count[rbin] += 1
105 
106 
107 print "#r Thick UpperHeight LowerHeight"
108 for i in range(rnum_bins):
109  r = rmin + (i+0.5)*rbin_width
110  upper = lower = 0.0
111  if upper_count[rbin] > 0:
112  upper = upper_sum[i] / upper_count[i]
113  if lower_count[rbin] > 0:
114  lower = lower_sum[i] / lower_count[i]
115  print r, upper-lower, upper, lower
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