Source code for MDAnalysis.analysis.dielectric

# -*- Mode: python; tab-width: 4; indent-tabs-mode:nil; coding: utf-8 -*-
# vim: tabstop=4 expandtab shiftwidth=4 softtabstop=4
#
# MDAnalysis --- https://www.mdanalysis.org
# Copyright (c) 2006-2017 The MDAnalysis Development Team and contributors
# (see the file AUTHORS for the full list of names)
#
# Released under the GNU Public Licence, v2 or any higher version
#
# Please cite your use of MDAnalysis in published work:
#
# R. J. Gowers, M. Linke, J. Barnoud, T. J. E. Reddy, M. N. Melo, S. L. Seyler,
# D. L. Dotson, J. Domanski, S. Buchoux, I. M. Kenney, and O. Beckstein.
# MDAnalysis: A Python package for the rapid analysis of molecular dynamics
# simulations. In S. Benthall and S. Rostrup editors, Proceedings of the 15th
# Python in Science Conference, pages 102-109, Austin, TX, 2016. SciPy.
# doi: 10.25080/majora-629e541a-00e
#
# N. Michaud-Agrawal, E. J. Denning, T. B. Woolf, and O. Beckstein.
# MDAnalysis: A Toolkit for the Analysis of Molecular Dynamics Simulations.
# J. Comput. Chem. 32 (2011), 2319--2327, doi:10.1002/jcc.21787
#

r"""
Dielectric --- :mod:`MDAnalysis.analysis.dielectric`
===========================================================

:Authors: Mattia Felice Palermo, Philip Loche
:Year: 2022
:Copyright: GNU Public License v3
"""

import numpy as np

from MDAnalysis.units import constants, convert
from MDAnalysis.analysis.base import AnalysisBase
from MDAnalysis.due import due, Doi
from MDAnalysis.exceptions import NoDataError

due.cite(Doi("10.1080/00268978300102721"),
         description="Dielectric analysis",
         path="MDAnalysis.analysis.dielectric",
         cite_module=True)
del Doi


[docs]class DielectricConstant(AnalysisBase): r""" Computes the average dipole moment .. math:: \boldsymbol M = \sum_i q_i \boldsymbol r_i where :math:`q_i` is the charge and :math:`\boldsymbol r_i` the position of atom :math:`i` in the given :class:`MDAnalysis.core.groups.AtomGroup`. Also, the static dielectric constant .. math:: \varepsilon = 1 + \frac{\langle M^2 \rangle - \langle M \rangle^2} {3 \varepsilon_ 0 V k_B T} is calculated for a system in tin foil boundary conditions, which is the usual case if electrostatics are handled with a Ewald summation technique. See [Neumann1983]_ for details on the derivation. Parameters ---------- atomgroup : MDAnalysis.core.groups.AtomGroup Atomgroup on which the analysis is executed temperature : float Temperature (Kelvin) at which the system has been simulated make_whole : bool Make molecules whole; If the input already contains whole molecules this can be disabled to gain speedup verbose : bool Show detailed progress of the calculation Attributes ---------- results.M : numpy.ndarray Directional dependant dipole moment :math:`\langle \boldsymbol M \rangle` in :math:`eÅ`. results.M2 : numpy.ndarray Directional dependant squared dipole moment :math:`\langle \boldsymbol M^2 \rangle` in :math:`(eÅ)^2` results.fluct : float Directional dependant dipole moment fluctuation :math:`\langle \boldsymbol M^2 \rangle - \langle \boldsymbol M \rangle^2` in :math:`(eÅ)^2` results.eps : numpy.ndarray Directional dependant static dielectric constant results.eps_mean : float Static dielectric constant Example ------- Create a :class:`DielectricConstant` instance by supplying an :class:`~MDAnalysis.core.groups.AtomGroup`, then use the :meth:`run` method:: import MDAnalysis as mda from MDAnalysis.analysis.dielectric import DielectricConstant from MDAnalysisTests.datafiles import PSF_TRICLINIC, DCD_TRICLINIC # Load a pure water system universe = mda.Universe(PSF_TRICLINIC, DCD_TRICLINIC) diel = DielectricConstant(universe.atoms) diel.run() print(diel.results) The static dielectric constant of the provided atomgroup is saved within the :class:`~MDAnalysis.analysis.base.Results` attribute. .. versionadded:: 2.1.0 """ def __init__(self, atomgroup, temperature=300, make_whole=True, **kwargs): super(DielectricConstant, self).__init__(atomgroup.universe.trajectory, **kwargs) self.atomgroup = atomgroup self.temperature = temperature self.make_whole = make_whole def _prepare(self): if not hasattr(self.atomgroup, "charges"): raise NoDataError("No charges defined given atomgroup.") if not np.allclose(self.atomgroup.total_charge(compound='fragments'), 0.0, atol=1E-5): raise NotImplementedError("Analysis for non-neutral systems or" " systems with free charges are not" " available.") self.volume = 0 self.results.M = np.zeros(3) self.results.M2 = np.zeros(3) self.results.fluct = np.zeros(3) self.results.eps = np.zeros(3) self.results.eps_mean = 0 def _single_frame(self): if self.make_whole: self.atomgroup.unwrap() self.volume += self.atomgroup.universe.trajectory.ts.volume M = np.dot(self.atomgroup.charges, self.atomgroup.positions) self.results.M += M self.results.M2 += M * M def _conclude(self): self.results.M /= self.n_frames self.results.M2 /= self.n_frames self.volume /= self.n_frames self.results.fluct = self.results.M2 - self.results.M * self.results.M self.results.eps = self.results.fluct / ( convert(constants["Boltzmann_constant"], "kJ/mol", "eV") * self.temperature * self.volume * constants["electric_constant"]) self.results.eps_mean = self.results.eps.mean() self.results.eps += 1 self.results.eps_mean += 1