MDAnalysis contains the MDAnalysis.analysis.hole module that uses Oliver Smart's HOLE program. This notebook shows how to make use of the basic functionality in MDAnalysis.analysis.hole
. We will analyze the fluctuations in the pore radius profile of an ion channel (gramicidin A) from a simulation trajectory.
import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline
plt.style.use('ggplot')
import MDAnalysis as mda
import MDAnalysis.analysis.hole
mda.start_logging()
MDAnalysis : INFO MDAnalysis 0.17.1-dev STARTED logging to 'MDAnalysis.log'
I am testing HOLE on a simple "trajectory": I submitted gramicidin A (1grm, test file 1gmr_single.pdb
) to the ElNemo elastic network server. From an ENM it computed a number of modes and provided multi-frame PDB files showing these modes. I just picked the first mode (which is labelled mode 7), which shows a twist motion along the pore axis. Note that this "trajectory" is exaggerated and does not correspond to real protein motion.
The trajectory is included with the MDAnalysis test files.
from MDAnalysis.tests.datafiles import MULTIPDB_HOLE
u = mda.Universe(MULTIPDB_HOLE)
/Volumes/Data/oliver/Biop/Projects/Methods/MDAnalysis/mdanalysis/testsuite/MDAnalysisTests/__init__.py:109: UserWarning: This call to matplotlib.use() has no effect because the backend has already been chosen; matplotlib.use() must be called *before* pylab, matplotlib.pyplot, or matplotlib.backends is imported for the first time. The backend was *originally* set to 'module://ipykernel.pylab.backend_inline' by the following code: File "/Users/oliver/anaconda3/envs/mda_develop/lib/python3.6/runpy.py", line 193, in _run_module_as_main "__main__", mod_spec) File "/Users/oliver/anaconda3/envs/mda_develop/lib/python3.6/runpy.py", line 85, in _run_code exec(code, run_globals) File "/Users/oliver/anaconda3/envs/mda_develop/lib/python3.6/site-packages/ipykernel_launcher.py", line 16, in <module> app.launch_new_instance() File "/Users/oliver/anaconda3/envs/mda_develop/lib/python3.6/site-packages/traitlets/config/application.py", line 658, in launch_instance app.start() File "/Users/oliver/anaconda3/envs/mda_develop/lib/python3.6/site-packages/ipykernel/kernelapp.py", line 486, in start self.io_loop.start() File "/Users/oliver/anaconda3/envs/mda_develop/lib/python3.6/site-packages/tornado/ioloop.py", line 888, in start handler_func(fd_obj, events) File "/Users/oliver/anaconda3/envs/mda_develop/lib/python3.6/site-packages/tornado/stack_context.py", line 277, in null_wrapper return fn(*args, **kwargs) File "/Users/oliver/anaconda3/envs/mda_develop/lib/python3.6/site-packages/zmq/eventloop/zmqstream.py", line 450, in _handle_events self._handle_recv() File "/Users/oliver/anaconda3/envs/mda_develop/lib/python3.6/site-packages/zmq/eventloop/zmqstream.py", line 480, in _handle_recv self._run_callback(callback, msg) File "/Users/oliver/anaconda3/envs/mda_develop/lib/python3.6/site-packages/zmq/eventloop/zmqstream.py", line 432, in _run_callback callback(*args, **kwargs) File "/Users/oliver/anaconda3/envs/mda_develop/lib/python3.6/site-packages/tornado/stack_context.py", line 277, in null_wrapper return fn(*args, **kwargs) File "/Users/oliver/anaconda3/envs/mda_develop/lib/python3.6/site-packages/ipykernel/kernelbase.py", line 283, in dispatcher return self.dispatch_shell(stream, msg) File "/Users/oliver/anaconda3/envs/mda_develop/lib/python3.6/site-packages/ipykernel/kernelbase.py", line 233, in dispatch_shell handler(stream, idents, msg) File "/Users/oliver/anaconda3/envs/mda_develop/lib/python3.6/site-packages/ipykernel/kernelbase.py", line 399, in execute_request user_expressions, allow_stdin) File "/Users/oliver/anaconda3/envs/mda_develop/lib/python3.6/site-packages/ipykernel/ipkernel.py", line 208, in do_execute res = shell.run_cell(code, store_history=store_history, silent=silent) File "/Users/oliver/anaconda3/envs/mda_develop/lib/python3.6/site-packages/ipykernel/zmqshell.py", line 537, in run_cell return super(ZMQInteractiveShell, self).run_cell(*args, **kwargs) File "/Users/oliver/anaconda3/envs/mda_develop/lib/python3.6/site-packages/IPython/core/interactiveshell.py", line 2739, in run_cell self.events.trigger('post_run_cell') File "/Users/oliver/anaconda3/envs/mda_develop/lib/python3.6/site-packages/IPython/core/events.py", line 73, in trigger func(*args, **kwargs) File "/Users/oliver/anaconda3/envs/mda_develop/lib/python3.6/site-packages/ipykernel/pylab/backend_inline.py", line 160, in configure_once activate_matplotlib(backend) File "/Users/oliver/anaconda3/envs/mda_develop/lib/python3.6/site-packages/IPython/core/pylabtools.py", line 308, in activate_matplotlib matplotlib.pyplot.switch_backend(backend) File "/Users/oliver/anaconda3/envs/mda_develop/lib/python3.6/site-packages/matplotlib/pyplot.py", line 229, in switch_backend matplotlib.use(newbackend, warn=False, force=True) File "/Users/oliver/anaconda3/envs/mda_develop/lib/python3.6/site-packages/matplotlib/__init__.py", line 1305, in use reload(sys.modules['matplotlib.backends']) File "/Users/oliver/anaconda3/envs/mda_develop/lib/python3.6/importlib/__init__.py", line 166, in reload _bootstrap._exec(spec, module) File "/Users/oliver/anaconda3/envs/mda_develop/lib/python3.6/site-packages/matplotlib/backends/__init__.py", line 14, in <module> line for line in traceback.format_stack() matplotlib.use('agg')
You need to install HOLE
first, which is available from http://www.holeprogram.org/ under the Apache open source license.
Set up the HOLEtraj
class; provide a path to the hole
executable.
H = MDAnalysis.analysis.hole.HOLEtraj(u, executable="~/hole2/exe/hole")
Then run hole
on the individual frames (this takes a while...)
H.run()
MDAnalysis.analysis.hole: INFO HOLE analysis frame 0 (orderparameter 0) MDAnalysis.analysis.hole: INFO Setting up HOLE analysis for '/var/folders/ds/z7mbf7nx1gn31308g31t9tcc0000gp/T/tmpc54h20_6.pdb' MDAnalysis.analysis.hole: INFO Using radius file 'simple2.rad' MDAnalysis.analysis.hole: INFO HOLE will guess CPOINT MDAnalysis.analysis.hole: INFO HOLE will guess CVECT MDAnalysis.analysis.hole: INFO Starting HOLE on '/var/folders/ds/z7mbf7nx1gn31308g31t9tcc0000gp/T/tmpc54h20_6.pdb' (trajectory: None) MDAnalysis.analysis.hole: INFO HOLE finished: output file 'hole.out' MDAnalysis.analysis.hole: INFO Collecting HOLE profiles for run with id 1 MDAnalysis.analysis.hole: INFO Run 1: Reading 1 HOLE profiles from 'hole.out' MDAnalysis.analysis.hole: INFO Collected HOLE radius profiles for 1 frames MDAnalysis.analysis.hole: INFO HOLE analysis frame 1 (orderparameter 1) MDAnalysis.analysis.hole: INFO Setting up HOLE analysis for '/var/folders/ds/z7mbf7nx1gn31308g31t9tcc0000gp/T/tmpd2ksqofz.pdb' MDAnalysis.analysis.hole: INFO Using radius file 'simple2.rad' MDAnalysis.analysis.hole: INFO HOLE will guess CPOINT MDAnalysis.analysis.hole: INFO HOLE will guess CVECT MDAnalysis.analysis.hole: INFO Starting HOLE on '/var/folders/ds/z7mbf7nx1gn31308g31t9tcc0000gp/T/tmpd2ksqofz.pdb' (trajectory: None) MDAnalysis.analysis.hole: INFO HOLE finished: output file 'hole.out' MDAnalysis.analysis.hole: INFO Collecting HOLE profiles for run with id 1 MDAnalysis.analysis.hole: INFO Run 1: Reading 1 HOLE profiles from 'hole.out' MDAnalysis.analysis.hole: INFO Collected HOLE radius profiles for 1 frames MDAnalysis.analysis.hole: INFO HOLE analysis frame 2 (orderparameter 2) MDAnalysis.analysis.hole: INFO Setting up HOLE analysis for '/var/folders/ds/z7mbf7nx1gn31308g31t9tcc0000gp/T/tmpxkp9n4md.pdb' MDAnalysis.analysis.hole: INFO Using radius file 'simple2.rad' MDAnalysis.analysis.hole: INFO HOLE will guess CPOINT MDAnalysis.analysis.hole: INFO HOLE will guess CVECT MDAnalysis.analysis.hole: INFO Starting HOLE on '/var/folders/ds/z7mbf7nx1gn31308g31t9tcc0000gp/T/tmpxkp9n4md.pdb' (trajectory: None) MDAnalysis.analysis.hole: INFO HOLE finished: output file 'hole.out' MDAnalysis.analysis.hole: INFO Collecting HOLE profiles for run with id 1 MDAnalysis.analysis.hole: INFO Run 1: Reading 1 HOLE profiles from 'hole.out' MDAnalysis.analysis.hole: INFO Collected HOLE radius profiles for 1 frames MDAnalysis.analysis.hole: INFO HOLE analysis frame 3 (orderparameter 3) MDAnalysis.analysis.hole: INFO Setting up HOLE analysis for '/var/folders/ds/z7mbf7nx1gn31308g31t9tcc0000gp/T/tmp4rezxrfz.pdb' MDAnalysis.analysis.hole: INFO Using radius file 'simple2.rad' MDAnalysis.analysis.hole: INFO HOLE will guess CPOINT MDAnalysis.analysis.hole: INFO HOLE will guess CVECT MDAnalysis.analysis.hole: INFO Starting HOLE on '/var/folders/ds/z7mbf7nx1gn31308g31t9tcc0000gp/T/tmp4rezxrfz.pdb' (trajectory: None) MDAnalysis.analysis.hole: INFO HOLE finished: output file 'hole.out' MDAnalysis.analysis.hole: INFO Collecting HOLE profiles for run with id 1 MDAnalysis.analysis.hole: INFO Run 1: Reading 1 HOLE profiles from 'hole.out' MDAnalysis.analysis.hole: INFO Collected HOLE radius profiles for 1 frames MDAnalysis.analysis.hole: INFO HOLE analysis frame 4 (orderparameter 4) MDAnalysis.analysis.hole: INFO Setting up HOLE analysis for '/var/folders/ds/z7mbf7nx1gn31308g31t9tcc0000gp/T/tmpgy7fauh3.pdb' MDAnalysis.analysis.hole: INFO Using radius file 'simple2.rad' MDAnalysis.analysis.hole: INFO HOLE will guess CPOINT MDAnalysis.analysis.hole: INFO HOLE will guess CVECT MDAnalysis.analysis.hole: INFO Starting HOLE on '/var/folders/ds/z7mbf7nx1gn31308g31t9tcc0000gp/T/tmpgy7fauh3.pdb' (trajectory: None) MDAnalysis.analysis.hole: INFO HOLE finished: output file 'hole.out' MDAnalysis.analysis.hole: INFO Collecting HOLE profiles for run with id 1 MDAnalysis.analysis.hole: INFO Run 1: Reading 1 HOLE profiles from 'hole.out' MDAnalysis.analysis.hole: INFO Collected HOLE radius profiles for 1 frames MDAnalysis.analysis.hole: INFO HOLE analysis frame 5 (orderparameter 5) MDAnalysis.analysis.hole: INFO Setting up HOLE analysis for '/var/folders/ds/z7mbf7nx1gn31308g31t9tcc0000gp/T/tmpou9fm86u.pdb' MDAnalysis.analysis.hole: INFO Using radius file 'simple2.rad' MDAnalysis.analysis.hole: INFO HOLE will guess CPOINT MDAnalysis.analysis.hole: INFO HOLE will guess CVECT MDAnalysis.analysis.hole: INFO Starting HOLE on '/var/folders/ds/z7mbf7nx1gn31308g31t9tcc0000gp/T/tmpou9fm86u.pdb' (trajectory: None) MDAnalysis.analysis.hole: INFO HOLE finished: output file 'hole.out' MDAnalysis.analysis.hole: INFO Collecting HOLE profiles for run with id 1 MDAnalysis.analysis.hole: INFO Run 1: Reading 1 HOLE profiles from 'hole.out' MDAnalysis.analysis.hole: INFO Collected HOLE radius profiles for 1 frames MDAnalysis.analysis.hole: INFO HOLE analysis frame 6 (orderparameter 6) MDAnalysis.analysis.hole: INFO Setting up HOLE analysis for '/var/folders/ds/z7mbf7nx1gn31308g31t9tcc0000gp/T/tmph9gqfrqx.pdb' MDAnalysis.analysis.hole: INFO Using radius file 'simple2.rad' MDAnalysis.analysis.hole: INFO HOLE will guess CPOINT MDAnalysis.analysis.hole: INFO HOLE will guess CVECT MDAnalysis.analysis.hole: INFO Starting HOLE on '/var/folders/ds/z7mbf7nx1gn31308g31t9tcc0000gp/T/tmph9gqfrqx.pdb' (trajectory: None) MDAnalysis.analysis.hole: INFO HOLE finished: output file 'hole.out' MDAnalysis.analysis.hole: INFO Collecting HOLE profiles for run with id 1 MDAnalysis.analysis.hole: INFO Run 1: Reading 1 HOLE profiles from 'hole.out' MDAnalysis.analysis.hole: INFO Collected HOLE radius profiles for 1 frames MDAnalysis.analysis.hole: INFO HOLE analysis frame 7 (orderparameter 7) MDAnalysis.analysis.hole: INFO Setting up HOLE analysis for '/var/folders/ds/z7mbf7nx1gn31308g31t9tcc0000gp/T/tmpzx_kni6p.pdb' MDAnalysis.analysis.hole: INFO Using radius file 'simple2.rad' MDAnalysis.analysis.hole: INFO HOLE will guess CPOINT MDAnalysis.analysis.hole: INFO HOLE will guess CVECT MDAnalysis.analysis.hole: INFO Starting HOLE on '/var/folders/ds/z7mbf7nx1gn31308g31t9tcc0000gp/T/tmpzx_kni6p.pdb' (trajectory: None) MDAnalysis.analysis.hole: INFO HOLE finished: output file 'hole.out' MDAnalysis.analysis.hole: INFO Collecting HOLE profiles for run with id 1 MDAnalysis.analysis.hole: INFO Run 1: Reading 1 HOLE profiles from 'hole.out' MDAnalysis.analysis.hole: INFO Collected HOLE radius profiles for 1 frames MDAnalysis.analysis.hole: INFO HOLE analysis frame 8 (orderparameter 8) MDAnalysis.analysis.hole: INFO Setting up HOLE analysis for '/var/folders/ds/z7mbf7nx1gn31308g31t9tcc0000gp/T/tmpk5c9fu7b.pdb' MDAnalysis.analysis.hole: INFO Using radius file 'simple2.rad' MDAnalysis.analysis.hole: INFO HOLE will guess CPOINT MDAnalysis.analysis.hole: INFO HOLE will guess CVECT MDAnalysis.analysis.hole: INFO Starting HOLE on '/var/folders/ds/z7mbf7nx1gn31308g31t9tcc0000gp/T/tmpk5c9fu7b.pdb' (trajectory: None) MDAnalysis.analysis.hole: INFO HOLE finished: output file 'hole.out' MDAnalysis.analysis.hole: INFO Collecting HOLE profiles for run with id 1 MDAnalysis.analysis.hole: INFO Run 1: Reading 1 HOLE profiles from 'hole.out' MDAnalysis.analysis.hole: INFO Collected HOLE radius profiles for 1 frames MDAnalysis.analysis.hole: INFO HOLE analysis frame 9 (orderparameter 9) MDAnalysis.analysis.hole: INFO Setting up HOLE analysis for '/var/folders/ds/z7mbf7nx1gn31308g31t9tcc0000gp/T/tmp6u9jackb.pdb' MDAnalysis.analysis.hole: INFO Using radius file 'simple2.rad' MDAnalysis.analysis.hole: INFO HOLE will guess CPOINT MDAnalysis.analysis.hole: INFO HOLE will guess CVECT MDAnalysis.analysis.hole: INFO Starting HOLE on '/var/folders/ds/z7mbf7nx1gn31308g31t9tcc0000gp/T/tmp6u9jackb.pdb' (trajectory: None) MDAnalysis.analysis.hole: INFO HOLE finished: output file 'hole.out' MDAnalysis.analysis.hole: INFO Collecting HOLE profiles for run with id 1 MDAnalysis.analysis.hole: INFO Run 1: Reading 1 HOLE profiles from 'hole.out' MDAnalysis.analysis.hole: INFO Collected HOLE radius profiles for 1 frames MDAnalysis.analysis.hole: INFO HOLE analysis frame 10 (orderparameter 10) MDAnalysis.analysis.hole: INFO Setting up HOLE analysis for '/var/folders/ds/z7mbf7nx1gn31308g31t9tcc0000gp/T/tmpdwimmp5j.pdb' MDAnalysis.analysis.hole: INFO Using radius file 'simple2.rad' MDAnalysis.analysis.hole: INFO HOLE will guess CPOINT MDAnalysis.analysis.hole: INFO HOLE will guess CVECT MDAnalysis.analysis.hole: INFO Starting HOLE on '/var/folders/ds/z7mbf7nx1gn31308g31t9tcc0000gp/T/tmpdwimmp5j.pdb' (trajectory: None) MDAnalysis.analysis.hole: INFO HOLE finished: output file 'hole.out' MDAnalysis.analysis.hole: INFO Collecting HOLE profiles for run with id 1 MDAnalysis.analysis.hole: INFO Run 1: Reading 1 HOLE profiles from 'hole.out' MDAnalysis.analysis.hole: INFO Collected HOLE radius profiles for 1 frames
The trajectory contains 11 frames. We are typically interested in the HOLE profile, i.e., a plot of the pore radius as determined by HOLE against the reaction coordinate (often close to $z$ along the pore axis for simple channels).
All profiles are stored as
H.profiles
which is a dict
. The keys are the frame numbers (or a user supplied reaction coordinate for the conformation). The values are record arrays with columns frame
(the reaction coordinate for the conformation), rxncoord
(the HOLE pore reaction coordinate), and radius
(radius in Angstrom).
The H.plot()
method plots all the profiles together, the H.plot3D()
stacks them:
ax = H.plot()
ax.set_ylim(1, 4)
ax.set_xlim(-15, 25)
plt.tight_layout()
H.plot3D(rmax=4)
<matplotlib.axes._subplots.Axes3DSubplot at 0x11deb7668>
The output from HOLE is a pore profile $R(\zeta)$ (radius $R$ vs pore coordinate $\zeta$).
All pore profiles are stored in a dict
named H.profiles
that is indexed by the frame number. Each pore profile itself is a array containig rows with frame number (all identical), pore coordinate, and pore radius. For more details, see the docs on HOLE data structures.
H.profiles
To get at individual profiles in we can iterate over H.profiles
. However, as a convenience, just iterating over the HOLEtraj
instance itself will return the profiles in frame-sorted order (thanks to HOLEtraj.sorted_profiles_iter()).
For example, getting the minimum radius
for frame, profile in H:
print("frame {}: min(R) = {}".format(frame, profile.radius.min()))
frame 0: min(R) = 1.61259 frame 1: min(R) = 1.56638 frame 2: min(R) = 1.5328 frame 3: min(R) = 1.42548 frame 4: min(R) = 1.24337 frame 5: min(R) = 1.19749 frame 6: min(R) = 1.2947 frame 7: min(R) = 1.44292 frame 8: min(R) = 1.51135 frame 9: min(R) = 1.52775 frame 10: min(R) = 1.55973
(For a slightly more advanced analysis of minimum pore radius as function of a order parameter see the Advanced HOLE Hacking Notebook.)
TODO: Change these links to Binder notebooks.
with open("profiles.xvg", "w") as xvg:
for frame, profile in H:
for _, zeta, radius in profile:
xvg.write("{0} {1}\n".format(zeta, radius))
xvg.write("&\n")
The corresponding output file looks like this:
print("".join(open("profiles.xvg").readlines()[465:485]))
23.287 12.19109 23.387 12.30165 23.487 12.3973 23.587 12.49308 23.687 12.5897 23.787 12.68495 23.887 12.79422 23.987 13.12993 24.087 14.08415 24.187 15.57652 24.287 16.7467 24.387 18.19425 & -22.41283 20.68026 -22.31283 18.25184 -22.21283 16.22098 -22.11283 14.19682 -22.01283 12.62771 -21.91283 11.46761 -21.81283 10.71226
One can use the above code to write the profiles in any format; for instance, it would be easy to write each frame's profile to a separate file by moving the with
statement inside the outer for
loop.
At the moment (MDAnalysis 0.17.0) does not have a method to compute an average HOLE profile. However, in the following I demonstrate how this can be accomplished using the numkit.timeseries.regularized_function (a function found in GromacsWrapper) but we can also use it independently (I copied it from numkit/timeseries.py):
import numpy
def regularized_function(x, y, func, bins=100, range=None):
"""Compute *func()* over data aggregated in bins.
``(x,y) --> (x', func(Y'))`` with ``Y' = {y: y(x) where x in x' bin}``
First the data is collected in bins x' along x and then *func* is
applied to all data points Y' that have been collected in the bin.
.. function:: func(y) -> float
*func* takes exactly one argument, a numpy 1D array *y* (the
values in a single bin of the histogram), and reduces it to one
scalar float.
.. Note:: *x* and *y* must be 1D arrays.
:Arguments:
x
abscissa values (for binning)
y
ordinate values (func is applied)
func
a numpy ufunc that takes one argument, func(Y')
bins
number or array
range
limits (used with number of bins)
:Returns:
F,edges
function and edges (``midpoints = 0.5*(edges[:-1]+edges[1:])``)
(This function originated as
:func:`recsql.sqlfunctions.regularized_function`.)
"""
_x = numpy.asarray(x)
_y = numpy.asarray(y)
if len(_x.shape) != 1 or len(_y.shape) != 1:
raise TypeError("Can only deal with 1D arrays.")
# setup of bins (taken from numpy.histogram)
if (range is not None):
mn, mx = range
if (mn > mx):
raise AttributeError('max must be larger than min in range parameter.')
if not numpy.iterable(bins):
if range is None:
range = (_x.min(), _x.max())
mn, mx = [float(mi) for mi in range]
if mn == mx:
mn -= 0.5
mx += 0.5
bins = numpy.linspace(mn, mx, bins+1, endpoint=True)
else:
bins = numpy.asarray(bins)
if (numpy.diff(bins) < 0).any():
raise ValueError('bins must increase monotonically.')
sorting_index = numpy.argsort(_x)
sx = _x[sorting_index]
sy = _y[sorting_index]
# boundaries in SORTED data that demarcate bins; position in bin_index is the bin number
bin_index = numpy.r_[sx.searchsorted(bins[:-1], 'left'),
sx.searchsorted(bins[-1], 'right')]
# naive implementation: apply operator to each chunk = sy[start:stop] separately
#
# It's not clear to me how one could effectively block this procedure (cf
# block = 65536 in numpy.histogram) because there does not seem to be a
# general way to combine the chunks for different blocks, just think of
# func=median
F = numpy.zeros(len(bins)-1) # final function
F[:] = [func(sy[start:stop]) for start,stop in zip(bin_index[:-1],bin_index[1:])]
return F,bins
The strategy is to collect all data and the histogram the radii over the HOLE reactioncoordinate:
rxncoord = np.concatenate([profile.rxncoord for frame, profile in H.sorted_profiles_iter()])
radii = np.concatenate([profile.radius for frame, profile in H.sorted_profiles_iter()])
Now use regularized_function
.
Normally I would import it as
from numkit.timeseries import regularized_function
but for this notebook it is not necessary to install numkit (pip install numkit
...).
The function to reduce the data in each bin is either just the mean (np.mean
) or the standard deviation (np.std
). (For more ideas what else to use see the numkit notes on coarse graining timeseries.)
mean_r, q = regularized_function(rxncoord, radii, np.mean, bins=100)
std_r, q = regularized_function(rxncoord, radii, np.std, bins=100)
zeta = 0.5*(q[1:] + q[:-1])
Note that regularized_function()
returns the reduced data in each bin $i$ and the bin edges (the bin $i$ is defined as $\{q: q_i \leq q < q_{i+1}\}$). To plot over the bin midpoints, we compute them as $\zeta_i = \frac{1}{2}(q_i + q_{i+1})$.
ax = plt.subplot(111)
ax.fill_between(zeta, mean_r - std_r, mean_r + std_r, color="blue", alpha=0.3)
ax.plot(zeta, mean_r, color="blue", lw=2)
ax.set_xlabel(r"pore coordinate $\zeta$ ($\AA$)")
ax.set_ylabel(r"HOLE radius $R$ ($\AA$)");
This is a plot of the mean pore radius (solid line) with the bands showing the standard deviation of the radius over the trajectory.