"""Virial isotherm model."""
import numpy
from scipy import optimize
from pygaps import logger
from pygaps.graphing.calc_graphs import virial_plot
from pygaps.modelling.base_model import IsothermBaseModel
from pygaps.utilities.exceptions import CalculationError
[docs]class Virial(IsothermBaseModel):
r"""
A virial isotherm model with 3 factors.
.. math::
p(n) = n \exp{(-\ln{K_H} + An + Bn^2 + Cn^3)}
Notes
-----
A virial isotherm model attempts to fit the measured data to a factorized
exponent relationship between loading and pressure.
.. math::
p = n \exp{(K_1n^0 + K_2n^1 + K_3n^2 + K_4n^3 + ... + K_i n^{i-1})}
It has been applied with success to describe the behaviour of standard as
well as supercritical isotherms. The factors are usually empirical,
although some relationship with physical can be determined:
the first constant is related to the Henry constant at zero loading, while
the second constant is a measure of the interaction strength with the surface.
.. math::
K_1 = -\ln{K_{H,0}}
In practice, besides the first constant, only 2-3 factors are used.
"""
# Model parameters
name = 'Virial'
formula = r"p(n) = n \exp{(-\ln{K_H} + An + Bn^2 + Cn^3)}"
calculates = 'pressure'
param_names = ("K", "A", "B", "C")
param_default_bounds = (
(0, numpy.inf),
(-numpy.inf, numpy.inf),
(-numpy.inf, numpy.inf),
(-numpy.inf, numpy.inf),
)
[docs] def loading(self, pressure):
"""
Calculate loading at specified pressure.
Careful!
For the Virial model, the loading has to
be computed numerically.
Parameters
----------
pressure : float
The pressure at which to calculate the loading.
Returns
-------
float
Loading at specified pressure.
"""
def fun(x):
return (self.pressure(x) - pressure)**2
opt_res = optimize.minimize(fun, pressure, method='Nelder-Mead')
if not opt_res.success:
raise CalculationError(f"Root finding failed. Error: \n\t{opt_res.message}")
return opt_res.x
[docs] def pressure(self, loading):
"""
Calculate pressure at specified loading.
The Virial model calculates the pressure directly.
Parameters
----------
loading : float
The loading at which to calculate the pressure.
Returns
-------
float
Pressure at specified loading.
"""
return loading * numpy.exp(
-numpy.log(self.params['K']) + self.params['A'] * loading +
self.params['B'] * loading**2 + self.params['C'] * loading**3
)
[docs] def spreading_pressure(self, pressure):
r"""
Calculate spreading pressure at specified gas pressure.
Function that calculates spreading pressure by solving the
following integral at each point i.
.. math::
\pi = \int_{0}^{p_i} \frac{n_i(p_i)}{p_i} dp_i
The integral for the Virial model cannot be solved analytically
and must be calculated numerically.
Parameters
----------
pressure : float
The pressure at which to calculate the spreading pressure.
Returns
-------
float
Spreading pressure at specified pressure.
"""
raise NotImplementedError
[docs] def initial_guess(self, pressure, loading):
"""
Return initial guess for fitting.
Parameters
----------
pressure : ndarray
Pressure data.
loading : ndarray
Loading data.
Returns
-------
dict
Dictionary of initial guesses for the parameters.
"""
saturation_loading, langmuir_k = super().initial_guess(pressure, loading)
guess = {"K": saturation_loading * langmuir_k, "A": 0, "B": 0, "C": 0}
guess = self.initial_guess_bounds(guess)
return guess
[docs] def fit(self, pressure, loading, param_guess, optimization_params=None, verbose=False):
"""
Fit model to data using nonlinear optimization with least squares loss function.
Resulting parameters are assigned to self.
Parameters
----------
pressure : ndarray
The pressures of each point.
loading : ndarray
The loading for each point.
optimization_params : dict
Custom parameters to pass to SciPy.optimize.least_squares.
verbose : bool, optional
Prints out extra information about steps taken.
"""
if verbose:
logger.info(f"Attempting to model using {self.name}")
# parameter names (cannot rely on order in Dict)
param_names = [param for param in self.params]
guess = numpy.array([param_guess[param] for param in param_names])
bounds = [[self.param_bounds[param][0] for param in param_names],
[self.param_bounds[param][1] for param in param_names]]
# remove invalid values in function
zero_values = ~numpy.logical_and(pressure > 0, loading > 0)
if any(zero_values):
logger.warning('Removed points which are equal to 0.')
pressure = pressure[~zero_values]
loading = loading[~zero_values]
# define fitting function as polynomial transformed input
ln_p_over_n = numpy.log(numpy.divide(pressure, loading))
# add point
add_point = False
added_point = False
if optimization_params:
add_point = optimization_params.pop('add_point', None)
fractional_loading = loading / max(loading)
if len(fractional_loading[fractional_loading < 0.5]) < 3:
if not add_point:
raise CalculationError(
"""
The isotherm recorded has very few points below 0.5
fractional loading. If a virial model fit is attempted
the resulting polynomial will likely be unstable in the
low loading region.
You can pass ``add_point=True`` in ``optimization_params``
to attempt to add a point in the low pressure region or
record better isotherms.
"""
)
added_point = True
ln_p_over_n = numpy.hstack([ln_p_over_n[0], ln_p_over_n])
loading = numpy.hstack([1e-1, loading])
def fit_func(x, L, ln_p_over_n):
for i, _ in enumerate(param_names):
self.params[param_names[i]] = x[i]
return self.params['C'] * L**3 + self.params['B'] * L**2 \
+ self.params['A'] * L - numpy.log(self.params['K']) - ln_p_over_n
kwargs = dict(
fun=fit_func, # fitting function
x0=guess, # initial guess
bounds=bounds, # bounds of the parameters
args=(loading, ln_p_over_n), # extra arguments to the fit function
# loss='huber', # use a loss function against outliers
# f_scale=0.1, # scale of outliers
)
if optimization_params:
kwargs.update(optimization_params)
opt_res = self.fit_leastsq(kwargs)
# assign params
for index, _ in enumerate(param_names):
self.params[param_names[index]] = opt_res.x[index]
self.rmse = numpy.sqrt(numpy.sum((opt_res.fun)**2) / len(loading))
if verbose:
logger.info(f"Model {self.name} success, RMSE is {self.rmse:.4g}")
n_load = numpy.linspace(1e-2, numpy.amax(loading), 100)
virial_plot(
loading, ln_p_over_n, n_load,
numpy.log(numpy.divide(self.pressure(n_load), n_load)), added_point
)