"""
MembraneCurvature
=======================================
:Author: Estefania Barreto-Ojeda
:Year: 2021
:Copyright: GNU Public License v3
MembraneCurvature calculates the mean and Gaussian curvature of
surfaces derived from a selection of reference.
Mean curvature is calculated in units of Å :sup:`-1` and Gaussian curvature
in units of Å :sup:`-2`.
"""
import numpy as np
import warnings
from .binning_surface import get_z_surface
from .curvature import (
mean_curvature,
gaussian_curvature,
fourier_curvature,
)
from .fourier_surface import n_fourier_parameters
from .fourier_validators import validate_positive_bin_counts
import MDAnalysis
from MDAnalysis.analysis.base import AnalysisBase
import logging
MDAnalysis.start_logging()
logger = logging.getLogger('MDAnalysis.MDAKit.membrane_curvature')
[docs]
class MembraneCurvature(AnalysisBase):
r"""
MembraneCurvature is a tool to calculate membrane curvature.
Parameters
----------
universe : Universe or AtomGroup
An MDAnalysis Universe object.
select : str or iterable of str, optional.
The selection string of an atom selection to use as a
reference to derive a surface.
wrap : bool or None, optional
Apply coordinate wrapping with :meth:`~MDAnalysis.core.groups.AtomGroup.wrap`
Defaults to ``True`` for ``surface_method='binning'`` and must be omitted or
``False`` for ``surface_method='fourier'``.
n_x_bins : int, optional, default: '100'
Number of bins in grid in the x dimension.
n_y_bins : int, optional, default: '100'
Number of bins in grid in the y dimension.
x_range : tuple of (float, float), optional, default: (0, `universe.dimensions[0]`)
Range of coordinates (min, max) in the x dimension.
y_range : tuple of (float, float), optional, default: (0, `universe.dimensions[1]`)
Range of coordinates (min, max) in the y dimension.
surface_method : {'binning', 'fourier'}, optional
``fourier`` (default) fits a periodic Fourier sum to atom positions at
each frame and evaluates Monge-gauge curvature from analytic derivatives
on the same bin grid (bin centers); see :mod:`membrane_curvature.fourier_surface`.
``binning`` bins atoms and uses :func:`numpy.gradient` for derivatives.
fourier_m : int, optional
Maximum Fourier mode index in ``x`` when ``surface_method='fourier'``.
Default ``2``.
fourier_n : int, optional
Maximum Fourier mode index in ``y`` when ``surface_method='fourier'``.
Default ``2``.
fourier_rcond : float, optional
Singular-value cutoff for the Fourier fit via truncated SVD
with :func:`~membrane_curvature.fourier_surface._solve_design_least_squares_svd`
when ``surface_method='fourier'``. The cutoff is interpreted as a
relative threshold on singular values.
Attributes
----------
results.z_surface : ndarray
Surface derived from atom selection in every frame.
Array of shape (`n_frames`, `n_x_bins`, `n_y_bins`)
results.mean : ndarray
Mean curvature associated with the surface.
Array of shape (`n_frames`, `n_x_bins`, `n_y_bins`)
results.gaussian : ndarray
Gaussian curvature associated with the surface.
Array of shape (`n_frames`, `n_x_bins`, `n_y_bins`)
results.average_z_surface : ndarray
Average of the array elements in `z_surface`.
Each array has shape (`n_x_bins`, `n_y_bins`)
results.average_mean : ndarray
Average of the array elements in `mean_curvature`.
Each array has shape (`n_x_bins`, `n_y_bins`)
results.average_gaussian: ndarray
Average of the array elements in `gaussian_curvature`.
Each array has shape (`n_x_bins`, `n_y_bins`)
Raises
------
ValueError
If ``n_x_bins`` or ``n_y_bins`` is not a positive integer
(see :func:`~membrane_curvature.fourier_validators.validate_positive_bin_counts`),
if the selection is empty, if ``surface_method`` is not one of
``'binning'`` or ``'fourier'``, or, when ``surface_method='fourier'``,
if ``fourier_m`` / ``fourier_n`` are negative or the selection has fewer
atoms than Fourier parameters, or if ``wrap=True`` with
``surface_method='fourier'``.
See also
--------
:class:`~MDAnalysis.transformations.wrap.wrap`
Wrap/unwrap the atoms of a given AtomGroup in the unit cell.
Notes
-----
**Fourier surface method (default)**
``surface_method='fourier'`` uses ``fourier_m = fourier_n = 2`` as default.
Do not modify the default values unless you need shorter wavelengths.
Since the method performs periodic boundary conditions by itself, ``wrap`` defaults
to ``False`` and is not required.
**Binning mode**
The binning routine does not apply periodic wrapping itself; :class:`MembraneCurvature`
calls ``AtomGroup.wrap()`` when ``surface_method='binning'`` and ``wrap=True``.
When using binning, ``wrap`` defaults to ``True`` if not provided. Omit ``wrap`` or pass
``wrap=True`` for raw trajectories so atoms are packed into the unit cell before
binning. Run with ``wrap=False`` for preprocessed trajectories that have already applied
periodic boundary conditions. For membrane-protein systems without position restraints,
preprocessing should include rotational and translational fitting around the protein.
For more details on when to use ``wrap=True``, check the :ref:`usage` page.
For any method of choice, the derived surface and calculated curvatures are available
in the :attr:`results` attributes.
The attribute :attr:`~MembraneCurvature.results.average_z_surface` contains
the derived surface averaged over the `n_frames` of the trajectory.
The attributes :attr:`~MembraneCurvature.results.average_mean` and
:attr:`~MembraneCurvature.results.average_gaussian` contain the computed
values of mean and Gaussian curvature averaged over the `n_frames` of the
trajectory.
See also
--------
:class:`~MDAnalysis.transformations.wrap.wrap`
Wrap/unwrap the atoms of a given AtomGroup in the unit cell.
Notes
-----
**Fourier surface method (default)**
``surface_method='fourier'`` uses ``fourier_m = fourier_n = 2`` as default.
Do not modify the default values unless you need shorter wavelengths.
Since method performs periodic boundary contidions by itself, ``wrap`` is not required.
**Binning mode**
``surface_method='binning'`` does not apply periodic wrapping. Use ``wrap=True``
on raw trajectories so atoms fall inside the grid. Use ``wrap=False`` for
pre-processed trajectories where periodic boundary conditions have been applied.
In membrane-protein systems, pre-processing requires rotational and translational fits.
For more details on when to use ``wrap=True``, check the :ref:`usage` page.
For any method of choice, the derived surface and calculated curvatures are available
in the :attr:`results` attributes.
The attribute :attr:`~MembraneCurvature.results.average_z_surface` contains
the derived surface averaged over the `n_frames` of the trajectory.
The attributes :attr:`~MembraneCurvature.results.average_mean` and
:attr:`~MembraneCurvature.results.average_gaussian` contain the computed
values of mean and Gaussian curvature averaged over the `n_frames` of the
trajectory.
Example
-----------
You can pass a universe containing your selection of reference::
import MDAnalysis as mda
from membrane_curvature.base import MembraneCurvature
u = mda.Universe(coordinates, trajectory)
mc = MembraneCurvature(u).run()
surface = mc.results.average_z_surface
mean_curvature = mc.results.average_mean
gaussian_curvature = mc.results.average_gaussian
The respective 2D curvature plots can be obtained using the `matplotlib`
package for data visualization via :func:`~matplotlib.pyplot.contourf` or
:func:`~matplotlib.pyplot.imshow`.
For specific examples visit the :ref:`usage` page.
Check the :ref:`visualization` page for more examples to plot
MembraneCurvature results using :func:`~matplotlib.pyplot.contourf`
and :func:`~matplotlib.pyplot.imshow`.
"""
def __init__(
self,
universe,
select='all',
n_x_bins=100,
n_y_bins=100,
x_range=None,
y_range=None,
wrap=None,
surface_method='fourier',
fourier_m=2,
fourier_n=2,
fourier_rcond=None,
**kwargs,
):
super().__init__(universe.universe.trajectory, **kwargs)
self.ag = universe.select_atoms(select)
# Validate selection up front so an empty AtomGroup produces a clear
# message before any downstream check (bin counts, surface method,
# Fourier-parameter sizing) can mask it.
if len(self.ag) == 0:
raise ValueError('Invalid selection. Empty AtomGroup.')
validate_positive_bin_counts(n_x_bins, n_y_bins)
self.n_x_bins = n_x_bins
self.n_y_bins = n_y_bins
self.x_range = x_range if x_range else (0, universe.dimensions[0])
self.y_range = y_range if y_range else (0, universe.dimensions[1])
self.dx = (self.x_range[1] - self.x_range[0]) / n_x_bins
self.dy = (self.y_range[1] - self.y_range[0]) / n_y_bins
valid_methods = ('binning', 'fourier')
if surface_method not in valid_methods:
raise ValueError(f'surface_method must be one of {valid_methods}, got {surface_method!r}')
self.surface_method = surface_method
if self.surface_method == 'fourier':
if wrap is True:
raise ValueError("wrap=True is only valid when surface_method='binning'")
self.wrap = False
else:
self.wrap = True if wrap is None else wrap
self.fourier_m = int(fourier_m)
self.fourier_n = int(fourier_n)
self.fourier_rcond = fourier_rcond
if self.surface_method == 'fourier':
if self.fourier_m < 0 or self.fourier_n < 0:
raise ValueError('fourier_m and fourier_n must be non-negative')
n_param = n_fourier_parameters(self.fourier_m, self.fourier_n)
n_atom = len(self.ag)
if n_atom < n_param:
max_square_mode = max(int((np.sqrt(n_atom) - 1) // 2), 0)
# If too few atoms for the passed modes fourier_m / fourier_n, raise error.
# The message suggests the maximum fourier_m / fourier_n given the number
# of atoms in the selection n_atom = len(self.ag)
raise ValueError(
f"surface_method='fourier' needs at least {n_param} atoms in the selection "
f'(fourier_m={self.fourier_m}, fourier_n={self.fourier_n}), got {n_atom}. '
f'Suggested max modes for {n_atom} atoms is '
f'fourier_m = fourier_n = {max_square_mode}.'
)
# Only checks the first frame. NPT simulations not properly covered here.
# Warning message if range doesn't cover entire dimensions of simulation box
for dim_string, dim_range, num in [('x', self.x_range, 0), ('y', self.y_range, 1)]:
if dim_range[1] < universe.dimensions[num]:
msg = (
f'Grid range in {dim_string} does not cover entire '
'dimensions of simulation box.\n Minimum dimensions '
'must be equal to simulation box.'
)
warnings.warn(msg)
logger.warning(msg)
# Warn when binning without wrapping coordinates
if not self.wrap and self.surface_method == 'binning':
# Warning
msg = (
' `wrap == False` may result in inaccurate calculation '
"of membrane curvature with `surface_method='binning'`. "
'Surfaces will be derived from a reduced number of atoms. \n '
' Ignore this warning if your trajectory has '
' rotational/translational fit rotations! '
)
warnings.warn(msg)
logger.warning(msg)
def _prepare(self):
# Initialize empty np.array with results
self.results.z_surface = np.full((self.n_frames, self.n_x_bins, self.n_y_bins), np.nan)
self.results.mean = np.full((self.n_frames, self.n_x_bins, self.n_y_bins), np.nan)
self.results.gaussian = np.full((self.n_frames, self.n_x_bins, self.n_y_bins), np.nan)
def _single_frame(self):
if self.surface_method == 'binning':
# Check wrap in binning method only
if self.wrap:
# Apply wrapping coordinates
self.ag.wrap()
self.results.z_surface[self._frame_index] = get_z_surface(
self.ag.positions,
n_x_bins=self.n_x_bins,
n_y_bins=self.n_y_bins,
x_range=self.x_range,
y_range=self.y_range,
)
self.results.mean[self._frame_index] = mean_curvature(
self.results.z_surface[self._frame_index], self.dx, self.dy
)
self.results.gaussian[self._frame_index] = gaussian_curvature(
self.results.z_surface[self._frame_index], self.dx, self.dy
)
else:
z_f, mean_f, gauss_f = fourier_curvature(
self.ag.positions,
self.x_range,
self.y_range,
self.n_x_bins,
self.n_y_bins,
self.fourier_m,
self.fourier_n,
rcond=self.fourier_rcond,
)
self.results.z_surface[self._frame_index] = z_f
self.results.mean[self._frame_index] = mean_f
self.results.gaussian[self._frame_index] = gauss_f
def _conclude(self):
self.results.average_z_surface = np.nanmean(self.results.z_surface, axis=0)
self.results.average_mean = np.nanmean(self.results.mean, axis=0)
self.results.average_gaussian = np.nanmean(self.results.gaussian, axis=0)