"""
Input/Output
"""
from __future__ import print_function
from dcma import tools
import os
import numpy as np
[docs]class Profiles(object):
"""
Diffusion and free energy profiles.
Diffusion is in 10^-5 cm^2/s, Free energy in kbT.
"""
def __init__(self, bin_edges, free_energy, diffusion):
"""
Args:
bin_edges (np.array): bin edges in Angstrom.
free_energy (np.array): Free energy profile in kbT.
diffusion (np.array): Diffusion profile in 10^-5 cm^2/s.
"""
self._diffusion = diffusion
self._free_energy = free_energy
self._bin_edges = bin_edges
self._bin_centers = tools.bin_edges_to_centers(bin_edges)
self._bin_widths = tools.bin_edges_to_widths(bin_edges)
self._n_bins = len(self.bin_centers)
assert len(self.bin_widths) == self.n_bins
assert len(diffusion) == self.n_bins
assert len(free_energy) == self.n_bins
@property
def diffusion(self):
""" Diffusion in 10^-5 cm^2/s. """
return self._diffusion
@property
def free_energy(self):
""" Free energy in kbT. """
return self._free_energy
@diffusion.setter
def diffusion(self, diffusion):
assert diffusion.shape == self._diffusion.shape
self._diffusion = diffusion
@free_energy.setter
def free_energy(self, free_energy):
assert free_energy.shape == self._free_energy.shape
self._free_energy = free_energy
@property
def bin_centers(self):
return self._bin_centers
@property
def bin_widths(self):
return self._bin_widths
@property
def bin_edges(self):
return self._bin_edges
@property
def n_bins(self):
return self._n_bins
def __len__(self):
return self._n_bins
def __eq__(self, other):
return (
(np.absolute(other.free_energy - self.free_energy) < 1e-10).all()
and (np.absolute(other.diffusion - self.diffusion) < 1e-10).all()
and (np.absolute(other.bin_centers - self.bin_centers) < 1e-10).all()
and (np.absolute(other.bin_widths - self.bin_widths) < 1e-10).all()
)
def __add__(self, other):
return Profiles(
bin_edges=np.average([self.bin_edges, other.bin_edges], axis=0),
free_energy=np.sum([self.free_energy, other.free_energy], axis=0),
diffusion=np.sum([self.diffusion, other.diffusion], axis=0),
)
def __radd__(self, other):
if other == 0:
return self
else:
return self + other
[docs] def find_membrane_bounds(self, threshold=0.05):
"""
Exclude water phase from the permeability calculation.
Args:
threshold:
start_index:
Returns: min_index, max_index
"""
d_gap = threshold*(self.diffusion.max() - self.diffusion.min())
fe_gap = threshold*(self.free_energy.max() - self.free_energy.min())
index = np.min(np.where(np.logical_or(
(np.absolute(self.diffusion - self.diffusion[0]) > d_gap),
(np.absolute(self.free_energy - self.free_energy[0]) > fe_gap)))
)
return index, self.n_bins - index
[docs] def concentration_in_water(self, num_permeants, area, water_width=None, threshold=0.05):
"""
Concentration of permeant in water in #molecules/A^3
Args:
num_permeants (int): Total number of permeants in the system
area (float): Cross section area of the membrane in Angstrom^2
water_width (int): Number of bins in the water phase on one side of the membrane
threshold (float): See calculate_permeability.
Returns:
float: The concentration of permeant in the water phase.
"""
integrand = np.exp(-(self.free_energy - self.reference_free_energy(water_width, threshold)))
integral = np.sum(self.bin_widths * integrand)
return num_permeants / area / integral
[docs] def concentration_in_membrane(self, num_permeants, area, skip=None, threshold=0.05):
"""
Concentration of permeant in the membrane in #molecules/A^3
Args:
num_permeants (int): Total number of permeants in the system
area (float): Cross section area of the membrane in Angstrom^2
skip (int): Number of bins in the water phase on one side of the membrane
threshold (float): See calculate_permeability.
Returns:
float: The average concentration of permeant in the membrane
"""
if skip is None:
imin, imax = self.find_membrane_bounds(threshold=0.05)
else:
imin, imax = skip, self.n_bins - skip
integrand = np.exp(-(self.free_energy - self.reference_free_energy(skip, threshold)))
integral_whole = np.sum(self.bin_widths * integrand)
integral_membrane = np.sum(self.bin_widths[imin:imax] * integrand[imin:imax])
length_membrane = np.sum(self.bin_widths[imin:imax])
volume_membrane = area * length_membrane
return num_permeants / volume_membrane * integral_membrane / integral_whole
[docs] def reference_free_energy(self, water_width=None, threshold=0.05):
"""
Reference free energy in the water phase.
Args:
water_width (int): Number of bins belonging to the water phase on one side of the membrane
threshold (float): see calculate_permeability
Returns:
The reference free energy.
"""
if water_width is None:
imin, imax = self.find_membrane_bounds(threshold)
else:
imin, imax = water_width, self.n_bins - water_width
water_fe = np.mean([self.free_energy[:imin]] + [self.free_energy[imax:]])
return water_fe
[docs] def calculate_permeability(self, skip=None, threshold=0.05):
"""
Permeability in cm/s
Args:
skip: Number of bins before membrane starts (overwrites threshold).
threshold: Relative deviation of free energy and diffusion profiles
that mark the end of the water phase and the start of
the membrane.
Returns:
A quadruple:
- Permeability in cm/s
- resistivity in s/cm,
- start index of membrane
- end index of membrane
"""
min_i, max_i = self.find_membrane_bounds(threshold)
reference_fe = self.free_energy[0]
if skip is not None:
min_i = skip
max_i = self.n_bins - min_i
integrand = (
np.exp(self.free_energy - reference_fe)
/ (self.diffusion * 1e-5)
* (1e-8 * self.bin_widths)
)
integral = np.sum(integrand[min_i:max_i])
return 1.0/integral, integrand, min_i, max_i
[docs] def save(self, filename):
p, integrand, min_i, max_i = self.calculate_permeability()
np.savetxt(
filename,
np.array([self.bin_centers,
self.free_energy,
self.diffusion,
self.bin_widths,
integrand]).transpose(),
header="Permeability/(cm/s): {} ".format(p)
+ os.linesep +
"( start_bin={} <= i < end_bin={} )".format(min_i, max_i)
+ os.linesep +
"bin_centers/(A) "
"free_energy/(kBT) "
"diffusion/(1e-5cm^2/s) "
"bin_widths/(A) "
"resistivity/(s/cm)",
fmt="%20.10f"
)
[docs] @staticmethod
def from_rate_matrix(rate):
raise NotImplementedError
[docs] def make_rate_matrix(self):
raise NotImplementedError
[docs] @staticmethod
def from_file(filename):
data = np.loadtxt(filename)
bin_centers = data[:,0]
free_energy = data[:,1]
diffusion = data[:,2]
bin_widths = data[:,3]
edges = bin_centers - bin_widths/2.0
bin_edges = np.append(edges, bin_centers[-1]+bin_widths[-1]/2.0)
return Profiles(bin_edges, free_energy, diffusion)
[docs]class SetOfProfiles(object):
"""
A set of profile objects to obtain statistics.
"""
def __init__(self, profiles_list, weights=None):
self.profiles = profiles_list
self.diffusions = np.array(
[p.diffusion for p in self.profiles]).transpose()
self.free_energies = np.array(
[p.free_energy for p in self.profiles]).transpose()
# unit converter to check consistency between bins and save bins
self.converter = tools.UnitConverter(
edges=[p.bin_edges for p in self.profiles],
weights=weights,
test_bin_consistency=True
)
self.weights = ([1.0 for p in self.profiles] if weights is None
else weights)
[docs] @staticmethod
def from_files(filenames, weights=None):
"""
Get profiles from files.
"""
profiles_list = []
for filename in filenames:
profiles_list += [Profiles.from_file(filename)]
return SetOfProfiles(profiles_list, weights)
def __len__(self):
return len(self.profiles)
def __getitem__(self, item):
return self.profiles[item]
[docs] def get_average(self):
"""
Returns:
An instance of Profiles containing the average
diffusion and free energy profiles.
"""
d = np.average(self.diffusions, weights=self.weights, axis=1)
e = np.average(self.free_energies, weights=self.weights, axis=1)
profiles = Profiles(self.converter.edges, e, d)
return profiles
[docs] def get_percentile(self, percent):
"""
Args:
percentile: A number between 0.0 and 100.0. Returns a set of profiles
so that the specified percentage of the data is below
the returned profile. E.g., for 95% confidence interval,
use percent=2.5 and percent=9.75.
Returns:
An instance of Profiles containing the quantiles.
"""
d = np.percentile(self.diffusions, percent, axis=1)
e = np.percentile(self.free_energies, percent, axis=1)
profiles = Profiles(self.converter.edges, e, d)
return profiles