DensityProfile¶

Bases: DynamicAnalysisBase

Serial and parallel implementations to calculate the number and charge density profiles \(\rho_i(z)\) and \(\rho_q(z)\) of a system along the specified axes.

The microscopic number density profile of species \(i\) in a constant-volume system is calculated by binning particle positions along an axis \(z\) using

\[\rho_i(z)=\frac{V}{N_\mathrm{bin}}\left\langle \sum_\alpha\delta(z-z_\alpha)\right\rangle\]

where \(V\) is the system volume and \(N_\mathrm{bins}\) is the number of bins. The angular brackets denote an ensemble average.

If the species carry charges, the charge density profile can be obtained using

\[\rho_q(z)=\sum_i z_ie\rho_i(z)\]

where \(z_i\) is the charge number of species \(i\) and \(e\) is the elementary charge.

With the charge density profile, the potential profile can be computed by numerically solving Poisson’s equation for electrostatics:

\[\varepsilon_0\varepsilon_\mathrm{r}\nabla^2\Psi(z)=-\rho_q(z)\]

Parameters:

groupsMDAnalysis.AtomGroup or array-like

Groups of atoms for which density profiles are calculated.

groupingsstr or array-like, default: "atoms"

Determines whether the centers of mass are used in lieu of individual atom positions. If groupings is a str, the same value is used for all groups.

Note

If the desired grouping is not "atoms",

the trajectory file should have segments (or chains) containing residues (or molecules) and residues containing atoms, and
residues and segments should be locally unwrapped at the simulation box edges, if not already, using MDAnalysis.transformations.wrap.unwrap, MDAnalysis.core.groups.AtomGroup.unwrap(), or MDAnalysis.lib.mdamath.make_whole().

Valid values:

"atoms": Atom positions (generally or for coarse-grained simulations).
"residues": Residues’ centers of mass (for atomistic simulations).
"segments": Segments’ centers of mass (for atomistic polymer simulations).

axesint, str, or array-like, default: "xyz"

Axes along which to compute the density profiles.

Examples:

2 for the \(z\)-direction.
"xy" for the \(x\)- and \(y\)-directions.
(0, 1) for the \(x\)- and \(y\)-directions.

n_binsint or array-like

Number of bins for each axis. If an int is provided, the same value is used for all axes. If parallel=True, the number of bins in all axes must be the same.

chargesarray-like, keyword-only, optional

Charge numbers \(z_i\) for the specified groupings in the \(N_\mathrm{g}\) groups. If not provided, it will be retrieved from the main MDAnalysis.core.universe.Universe object if available.

Note

Depending on the grouping for a specific group, all atoms, residues, or segments should have the same charge since the charge density profile for the group would not make sense otherwise. If this condition does not hold, change how the particles are grouped in grouping such that all entities share the same charge.

Shape: \((N_\mathrm{g})\).

Reference unit: \(\mathrm{e}\).

dimensionsarray-like, keyword-only, optional

System dimensions. If the MDAnalysis.core.universe.Universe object that the groups in groups belong to does not contain dimensionality information, provide it here. Affected by scales.

Shape: \((3,)\).

Reference unit: \(\mathrm{Å}\).

dtfloat, openmm.unit.Quantity, or pint.Quantity, keyword-only, optional

Time between frames \(\Delta t\). While this is normally determined from the trajectory, the trajectory may not have the correct information if the data is in reduced units. For example, if your reduced timestep is \(0.01\) and you output trajectory data every \(10,000\) timesteps, then \(\Delta t = 100\).

Reference unit: \(\mathrm{ps}\).

scalesarray-like, keyword-only, optional

Scaling factors for each system dimension. If an int is provided, the same value is used for all axes.

Shape: \((3,)\).

averagebool, keyword-only, default: True

Determines whether the density profiles are averaged over the specified frames.

recenterint, list, MDAnalysis.AtomGroup, or tuple, keyword-only, optional

Constrains the center of mass of an atom group by adjusting the particle coordinates every analysis frame. Either specify an MDAnalysis.core.groups.AtomGroup, its index within groups, or a tuple containing the aforementioned information and the fixed center of mass coordinates, in that order. To avoid recentering in a specific dimension, set the coordinate to numpy.nan. If the center of mass is not specified, the center of the simulation box is used.

Shape: \((3,)\) for the fixed center of mass.

reducedbool, keyword-only, default: False

Specifies whether the data is in reduced units. Affects results.number_densities, results.charge_densities, etc.

parallelbool, keyword-only, default: False

Determines whether the analysis is performed in parallel.

verbosebool, keyword-only, default: True

Determines whether detailed progress is shown.

**kwargs

Additional keyword arguments to pass to MDAnalysis.analysis.base.AnalysisBase.

Attributes:

universeMDAnalysis.Universe

MDAnalysis.core.universe.Universe object containing all information describing the simulation system.

results.unitsdict

Reference units for the results. For example, to get the reference units for results.bins, call results.units["results.bins"].

results.timesnumpy.ndarray

Times at which the density profiles are calculated.

Shape: \((N_\mathrm{frames},)\).

Reference unit: \(\mathrm{ps}\).

results.binslist

Bin centers corresponding to the density profiles in each dimension.

Shape: \((N_\mathrm{axes},)\) list of \((N_\mathrm{bins},)\) arrays.

Reference unit: \(\mathrm{Å}\).

results.number_densitieslist

Number density profiles.

Shape: \((N_\mathrm{axes},)\) list of \((N_\mathrm{bins},)\) arrays.

Reference unit: \(\mathrm{Å}^{-3}\).

results.charge_densitieslist

Charge density profiles, if charge information is available.

Shape: \((N_\mathrm{axes},)\) list of \((N_\mathrm{bins},)\) arrays.

Reference unit: \(\mathrm{e/Å}^{-3}\).

results.potentialsdict

Potential profiles, if charge information is available, with the key being the axis index (e.g., 1 for the \(z\)- direction if axes="yz"). Only available after running calculate_potential_profile().

Shape: \((N_\mathrm{bins},)\) for the potential profiles.

Reference unit: \(\mathrm{V}\).

Methods

`calculate_potential_profile`	Calculates the average potential profile in the given dimension using the charge density profile by numerically solving Poisson's equation for electrostatics.
`run`	Performs the calculation.
`save`	Saves results to a binary or archive file in NumPy format.

Calculates the average potential profile in the given dimension using the charge density profile by numerically solving Poisson’s equation for electrostatics.

Parameters:

dielectricfloat

Relative permittivity or dielectric constant \(\varepsilon_\mathrm{r}\).

axisint or str

Axis along which to compute the potential profiles.

Examples:

2 for the \(z\)-direction.
"x" for the \(x\)-direction.

sigma_qfloat, openmm.unit.Quantity, or pint.Quantity, keyword-only, optional

Total surface charge density \(\sigma_q\). Used to ensure that the electric field in the bulk of the solution is zero. If not provided, it is determined using dV and the charge density profile, or the average value in the center of the integrated charge density profile.

Reference unit: \(\mathrm{e/Å^2}\).

dVfloat, openmm.unit.Quantity, or pint.Quantity, keyword-only, optional

Potential difference \(\Delta\Psi\) across the system dimension specified in axis. Has no effect if sigma_q is provided since this value is used solely to calculate sigma_q.

Reference unit: \(\mathrm{V}\).

thresholdfloat, keyword-only, default: 1e-5

Threshold for determining the plateau region of the first integral of the charge density profile to calculate sigma_q. Has no effect if sigma_q is provided, or if sigma_q can be calculated using dV and dielectric.

V0float, openmm.unit.Quantity, or pint.Quantity, keyword-only, default: 0

Potential \(\Psi_0\) at the left boundary.

Reference unit: \(\mathrm{V}\).

methodstr, keyword-only, default: "integral"

Method to use to calculate the potential profile.

Valid values: "integral", "matrix".

pbcbool, keyword-only, default: False

Specifies whether the system has periodic boundary conditions. Only used when method="matrix".

run(start: int = None, stop: int = None, step: int = None, frames: slice | ndarray[int] = None, verbose: bool = None, **kwargs) → SerialAnalysisBase | ParallelAnalysisBase¶

Performs the calculation.