calculate_potential_profile

mdhelper.analysis.profile.calculate_potential_profile(bins: ndarray[float], charge_density: ndarray[float], L: float, dielectric: float = 1, *, sigma_q: float = None, dV: float = None, threshold: float = 1e-05, V0: float = 0, method: str = 'integral', pbc: bool = False, reduced: bool = False) → ndarray[float][source]

Calculates the potential profile \(\Psi(z)\) using the charge density profile by numerically solving Poisson’s equation for electrostatics.

Poisson’s equation is given by

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

where \(\varepsilon_0\) is the vacuum permittivity, \(\varepsilon_\mathrm{r}\) is the relative permittivity, \(\rho_q\) is the charge density, and \(\Psi\) is the potential.

The boundary conditions (BCs) are

\[\begin{split}\left.\frac{\partial\Psi}{\partial z}\right\rvert_{z=0}& =-\frac{\sigma_q}{\varepsilon_0\varepsilon_\mathrm{r}}\\ \left.\Psi\right\rvert_{z=0}&=\Psi_0\end{split}\]

The first BC is used to ensure that the electric field in the bulk of the solution is zero, while the second BC is used to set the potential on the left electrode.

Poisson’s equation can be evaluated by using the trapezoidal rule to numerically integrate the charge density profile twice:

  1. Integrate the charge density profile.

  2. Apply the first BC by adding \(\sigma_q\) to all points.

  3. Integrate the profile from Step 2.

  4. Divide by \(-\varepsilon_0\varepsilon_\mathrm{r}\).

  5. Apply the second BC by adding \(\Psi_0\) to all points.

This method is fast but requires many histogram bins to accurately resolve the potential profile.

Alternatively, Poisson’s equation can be written as a system of linear equations

\[A\mathbf{\Psi}=\mathbf{b}\]

using second-order finite differences. \(A\) is a tridiagonal matrix and \(b\) is a vector containing the charge density profile, with boundary conditions applied in the first and last rows.

The inner equations are given by

\[\Psi_i^{''}=\frac{\Psi_{i-1}-2\Psi_i+\Psi_{i+1}}{h^2} =-\frac{1}{\varepsilon_0\varepsilon_\mathrm{r}}\rho_{q,\,i}\]

where \(i\) is the bin index and \(h\) is the bin width.

In the case of periodic boundary conditions, the first and last equations are given by

\[\begin{split}\Psi_0^{''}&=\frac{\Psi_{N-1}-2\Psi_0+\Psi_1}{h^2} =-\frac{1}{\varepsilon_0\varepsilon_\mathrm{r}}\rho_{q,\,0}\\ \Psi_{N-1}^{''}&=\frac{\Psi_{N-2}-2\Psi_{N-1}+\Psi_0}{h^2} =-\frac{1}{\varepsilon_0\varepsilon_\mathrm{r}}\rho_{q,\,N-1}\end{split}\]

When the system has slab geometry, the boundary conditions are implemented via

\[\begin{split}\Psi_0&=\Psi_0\\ \Psi_0^\prime&=\frac{-3\Psi_0+4\Psi_1-\Psi_2}{2h} =-\frac{\sigma_q}{\varepsilon_0\varepsilon_\mathrm{r}}\end{split}\]

This method is slower but can be more accurate even with fewer histogram bins for bulk systems with periodic boundary conditions.

Parameters:
binsarray-like

Histogram bin centers corresponding to the charge density profile in charge_density.

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

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

charge_densityarray-like

Array containing the charge density profile.

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

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

Lfloat

System size in the dimension that bins and charge_density were calculated in.

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

dielectricfloat, default: 1

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

sigma_qfloat, 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 if method="integral".

Note

This value should be negative if the potential difference is positive and vice versa.

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

dVfloat, 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, 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 axis has periodic boundary conditions. Only used when method="matrix".

reducedbool, keyword-only, default: False

Specifies whether the data is in reduced units.

Returns:
potentialnumpy.ndarray

Potential profile \(\Psi(z)\).

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

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