calculate_potential_profile

mdcraft.analysis.profile.calculate_potential_profile(bins: ndarray[float], charge_density_profile: ndarray[float], dielectric: float, *, L: float = None, sigma_q: float | ndarray[float] = None, dV: float | ndarray[float] = None, threshold: float = 1e-05, V0: float | ndarray[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 \(\rho_q(z)\) 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\) and \(\varepsilon_\mathrm{r}\) are the vacuum and relative permittivities, respectively.

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 subtracting \(\sigma_q\) from 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 \(\mathbf{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 \(z\) for the charge density profile in charge_density_profile.

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

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

charge_density_profilearray-like

Charge density profile \(\rho_q(z)\).

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

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

dielectricfloat

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

Lfloat, keyword-only, optional

System size in the dimension that bins and charge_density_profile were calculated in. If not specified, it is determined using the first and last values in bins, assuming that the bin centers are uniformly spaced.

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

sigma_qfloat or array-like, keyword-only, optional

Surface charge density \(\sigma_q\). Used as a boundary condition. If not provided, it is determined using dV and charge_density_profile. If dV is also not provided, the average value in the center of the integrated charge density profile is used if method="integral".

Note

\(\sigma_q\) and \(\Delta\Psi\) should have the same sign.

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

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

dVfloat or array-like, keyword-only, optional

Potential difference \(\Delta\Psi\) across the system dimension specified in axis. Only used if sigma_q was not provided.

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

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

V0float or array-like, keyword-only, default: 0

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

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

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

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

Threshold for determining the plateau regions in the centers of the integrated charge density profiles to be used as estimates of sigma_q. Only used if sigma_q was not provided and cannot be calculated using dV and charge_density_profile.

methodstr, keyword-only, default: "integral"

Method used to calculate the potential profiles.

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: Same as charge_density_profile.

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