Note
Click here to download the full example code
Ideal component of the chemical potential#
import MDAnalysis as mda
import numpy as np
from scipy import constants as const
import maicos
def mu(rho, T, m):
"""Calculate the chemical potential.
The chemical potential is calculated from the
density: mu = k_B T log(rho. / m)
"""
# kT in KJ/mol
kT = T * const.Boltzmann * const.Avogadro / const.kilo
# De Broglie (converted to angstrom)
db = np.sqrt(
const.h ** 2 / (2 * np.pi * m * const.atomic_mass * const.Boltzmann * T)
) / const.angstrom
if np.all(rho > 0):
return kT * np.log(rho * db ** 3)
elif np.any(rho == 0):
return np.float64("-inf") * np.ones(rho.shape)
else:
return np.float64("nan") * np.ones(rho.shape)
def dmu(rho, drho, T):
"""Calculate the error of the chemical potential.
The error is calculated from the density using propagation of uncertainty.
"""
kT = T * const.Boltzmann * const.Avogadro / const.kilo
if np.all(rho > 0):
return kT * (drho / rho)
else:
return np.float64("nan") * np.ones(rho.shape)
water = mda.Universe('water.tpr', 'water.trr')
ana = maicos.DensityPlanar(water.atoms)
ana.run()
print('µ_id =', mu(ana.results.profile.mean(), 300, 18))
print('Δµ_id =', dmu(ana.results.profile.mean(),
ana.results.dprofile.mean(),
300))
Unwrapping in combination with atom grouping is superfluous. `unwrap` will be set to `False`.
µ_id = -12.05027764676635
Δµ_id = 0.041243491788978626