PartialElectronDensity¶

class PartialElectronDensity(configuration, kpoints=None, band_indices=None, energy_range=None, spectrum_method=None, density_mesh_cutoff=None)¶

Class for calculating the partial electron density.

Parameters:

configuration (BulkConfiguration | MoleculeConfiguration) – The configuration with an attached calculator for which to calculate the partial electron density.
kpoints – The k-points for which to perform the calculation. This can either be given as a k-point grid or a list of fractional k-points, e.g., [[0.0, 0.0, 0.0], [0.0, 0.0, 0.1], ...].
Default: The Monkhorst-Pack grid used for the self-consistent calculation.
band_indices (list of int) – The band indices of the Bloch states to include in the partial density.
Default: All occupied bands
energy_range (list of two PhysicalQuantity of type energy.) – Specifies the band energy range, relative to the Fermi level, of the bands to include in the partial density. Needs to be an interval or None.
Default: All energies below the Fermi level.
spectrum_method (GaussianBroadening) – The method to use for calculating the partial electron density.
Default: GaussianBroadening(0.1*eV)
density_mesh_cutoff (PhysicalQuantity of type energy | GridSampling | OptimizedFFTGridSampling) – The mesh cutoff to be used to determine the grid sampling. The mesh cutoff must be a positive energy or a GridSampling object.
Default: Specific for each calculator.

absolute()¶

Returns:: A new grid containing the absolute values (or modulus) of the current field.
Return type:: GridValues

axisProjection(projection_type='sum', axis='c', spin=None, projection_point=None, coordinate_type=<class 'NL.ComputerScienceUtilities.NLFlag._NLFlag.Fractional'>)¶

Get the values projected on one of the grid axes.

Parameters:

projection_type (str) –
The type of projection to perform. Should be either
- ’sum’ for the sum over the plane spanned by the two other axes.
- ’average’ or ‘avg’ for the average value over the plane spanned by the two other axes.
- ’line’ for the value along a line parallel to the axis and through a point specified by the projection_point parameter.
Default: ‘sum’
axis (str) – The axis to project the data onto. Should be either ‘a’, ‘b’ or ‘c’.
Default: ‘c’
spin (Spin.Sum | Spin.Z | Spin.X | Spin.Y | Spin.Up | Spin.Down | Spin.RealUpDown | Spin.ImagUpDown) – Which spin component to project on.
Default: Spin.All
projection_point (sequence, PhysicalQuantity) – Axis coordinates of the point through which to take a line if projection_type is ‘projection_point’. Must be given as a sequence of three coordinates [a, b, c]. It the numbers have units of length, they are first divided by the length of the respective primitive vectors [A, B, C], and then interpreted as fractional coordinates. Unitless coordinates are immidiately interpreted as fractional.
coordinate_type (Fractional | Cartesian) – Flag to toggle if the returned axis values should be given in units of Angstrom (NLFlag.Cartesian) or in units of the norm of the axis primitive vector (NLFlag.Fractional).
Default: Fractional

Returns:

A 2-tuple of 1D numpy.arrays containing the axis values and the projected data. For Cartesian coordinate type the grid offset is added to the axis values.

Return type:

tuple.

bandIndices()¶

Returns:: The band indicies to include in the partial electron density.
Return type:: list of int.

derivatives(x, y, z, spin=None)¶

Calculate the derivative in the point (x, y, z).

Parameters:

x (PhysicalQuantity with type length) – The Cartesian x coordinate.
y (PhysicalQuantity with type length) – The Cartesian y coordinate.
z (PhysicalQuantity with type length) – The Cartesian z coordinate.
spin (Spin.All | Spin.Sum | Spin.Up | Spin.Down | Spin.X | Spin.Y | Spin.Z) – The spin component to project on.
Default: Spin.All

Returns:

The gradient at the specified point for the given spin. For Spin.All, a tuple with (Spin.Sum, Spin.X, Spin.Y, Spin.Z) components is returned.

Return type:

PhysicalQuantity of type length^-4

downsample(downsampling_a=None, downsampling_b=None, downsampling_c=None)¶

Generate a new GridValues object where the grid is downsampled. Along periodic directions an FFT downsampling is performed. Along non-periodic directions antialiasing and downsampling is performed.

Parameters:

downsampling_a (int) – The new number of grid points along the A direction.
Default: No downsampling.
downsampling_b (int) – The new number of grid points along the B direction.
Default: No downsampling.
downsampling_c (int) – The new number of grid points along the C direction.
Default: No downsampling.

energyRange()¶

Returns:: The band energy range, relative to the Fermi level, of the bands to include in the partial density.
Return type:: list of two PhysicalQuantity of type energy.

evaluate(x, y, z, spin=None)¶

Evaluate in the point (x, y, z).

Parameters:

x (PhysicalQuantity with type length) – The Cartesian x coordinate.
y (PhysicalQuantity with type length) – The Cartesian y coordinate.
z (PhysicalQuantity with type length) – The Cartesian z coordinate.
spin (Spin.All | Spin.Sum | Spin.Up | Spin.Down | Spin.X | Spin.Y | Spin.Z) – The spin component to project on.
Default: Spin.All

Returns:

The value at the specified point for the given spin. For Spin.All, a tuple with (Spin.Sum, Spin.X, Spin.Y, Spin.Z) components is returned.

Return type:

PhysicalQuantity of type length^-3

gridCoordinate(i, j, k)¶

Return the coordinate for a given grid index.

Parameters:

i (int) – The grid index in the A direction.
j (int) – The grid index in the B direction.
k (int) – The grid index in the C direction.

Returns:

The Cartesian coordinate of the given grid index.

Return type:

PhysicalQuantity of type length.

kpoints()¶

Returns:: The kpoints given as input.
Return type:: MonkhorstPackGrid | RegularKpointGrid | numpy.ndarray

metatext()¶

Returns:: The metatext of the object or None if no metatext is present.
Return type:: str | None

nlprint(stream=<_io.TextIOWrapper name='<stdout>' mode='w' encoding='utf-8'>)¶

Print a string containing an ASCII table useful for plotting the partial electron density.

Parameters:: stream (file-like) – The stream to write to. This should be an object that supports strings being written to it using a write method.
Default: sys.stdout

primitiveVectors()¶

Returns:: The primitive vectors of the grid.
Return type:: PhysicalQuantity of type length.

scale(scale)¶

Scale the field with a float.

Parameters:: scale (float) – The parameter to scale with.

setMetatext(metatext)¶

Set a given metatext string on the object.

Parameters:: metatext (str | None) – The metatext string that should be set. A value of “None” can be given to remove the current metatext.

shape()¶

Returns:: The number of grid points in each direction.
Return type:: tuple of three int.

spinProjection(spin=None)¶

Construct a new GridValues object with the values of this object projected on a given spin component.

Parameters:: spin (Spin.All | Spin.Sum | Spin.X | Spin.Y | Spin.Z) – The spin component to project on.
Default: Spin.All
Returns:: A new GridValues object for the specified spin.
Return type:: GridValues

toArray()¶

Returns:: The values of the grid as a numpy array slicing off any units.
Return type:: numpy.array

uniqueString()¶: Return a unique string representing the state of the object.

unit()¶

Returns:: The unit of the data in the grid.
Return type:: A physical unit.

unitCell()¶

Returns:: The unit cell of the grid.
Return type:: PhysicalQuantity of type length.

volumeElement()¶

Returns:: The volume element of the grid represented by three vectors.
Return type:: PhysicalQuantity of type length.

Usage Examples¶

Analyze the partial electron density of a silver (111) surface. This structure is characterized by a surface state around the \(\Gamma\)-point. In the example, the kpoint sampling is specified with a RegularKpointGrid in a kpoint region around the \(\Gamma\)-point.

# -------------------------------------------------------------
# Bulk Configuration
# -------------------------------------------------------------

# Set up lattice
vector_a = [2.88903, 0.0, 0.0]*Angstrom
vector_b = [-1.44451, 2.50197, 0.0]*Angstrom
vector_c = [0.0, 0.0, 40.0]*Angstrom
lattice = UnitCell(vector_a, vector_b, vector_c)

# Define elements
elements = [Silver, Silver, Silver, Silver, Silver, Silver, Silver, Silver,
            Silver, Silver, Silver, Silver]

# Define coordinates
fractional_coordinates = [[ 0.749999669084,  0.5           ,  0.175654000711],
                          [ 0.749999669084,  0.5           ,  0.352570000323],
                          [ 0.749999669084,  0.5           ,  0.529485999935],
                          [ 0.749999669084,  0.5           ,  0.706401999547],
                          [ 0.416666424719,  0.833333341397,  0.234626000582],
                          [ 0.416666424719,  0.833333341397,  0.411542000194],
                          [ 0.416666424719,  0.833333341397,  0.588457999806],
                          [ 0.416666424719,  0.833333341397,  0.765373999418],
                          [ 0.083334237113,  0.166666658603,  0.293598000453],
                          [ 0.083334237113,  0.166666658603,  0.470514000065],
                          [ 0.083334237113,  0.166666658603,  0.647429999677],
                          [ 0.083334237113,  0.166666658603,  0.824345999289]]

# Set up configuration
bulk_configuration = BulkConfiguration(
    bravais_lattice=lattice,
    elements=elements,
    fractional_coordinates=fractional_coordinates
    )

# -------------------------------------------------------------
# Calculator
# -------------------------------------------------------------
#----------------------------------------
# Basis Set
#----------------------------------------
basis_set = [
    LDABasis.Silver_SingleZetaPolarized,
    ]

#----------------------------------------
# Exchange-Correlation
#----------------------------------------
exchange_correlation = LDA.PZ

k_point_sampling = MonkhorstPackGrid(
    na=11,
    nb=11,
    )
numerical_accuracy_parameters = NumericalAccuracyParameters(
    density_mesh_cutoff=45.0*Hartree,
    k_point_sampling=k_point_sampling,
    )

calculator = LCAOCalculator(
    basis_set=basis_set,
    exchange_correlation=exchange_correlation,
    numerical_accuracy_parameters=numerical_accuracy_parameters,
    )

bulk_configuration.setCalculator(calculator)
nlprint(bulk_configuration)
bulk_configuration.update()

# -------------------------------------------------------------
# Partial Electron Density
# -------------------------------------------------------------
kpoints = RegularKpointGrid(
    ka_range=[-0.05, 0.05],
    kb_range=[-0.05, 0.05],
    kc_range=[0.0, 0.0],
    na=5,
    nb=5,
    nc=1,
    )

partial_electron_density = PartialElectronDensity(
    configuration=bulk_configuration,
    kpoints=kpoints,
    energy_range=[-0.5, 0.5]*eV,
    band_indices=All,
    spectrum_method=GaussianBroadening(0.1*eV),
    )

partial_electron_density.py

Notes¶

The PartialElectronDensity can be calculated by specifying the

energy_range. Only electronic states with energies in the specified energy range \([E_{min}; E_{max}]\) will contribute to the partial electron density. The Fermi occupation factor used to construct the full electron density is not taken into account for the partial electron density. If \(E_{max}>0\) unoccupied states will therefore contribute to the partial density.

band_indices. Only states from the specified band indices will contribute to the partial electron density. By default, all bands are taken into account and the selection is determined by the energy range.

kpoints. The partial electron density is calculated from the eigenstates at the specified k-points.

spectrum_method. The method used to integrate the density. Currently only GaussianBroadening is supported. Each electronic state is assigned a weight which is computed by integrating a Gaussian centered at the state’s energy over the specified energy range.

The partial electron density is defined as:

\[n(\mathbf{r}) = \sum_{i}\sum_\mathbf{k} w_{\mathbf{k}} |\psi_{i,\mathbf{k}}(\mathbf{r})|^2 w(\epsilon_{i,\mathbf{k}}),\]

where the sum over band indices \(i\) includes only the specified bands (by default all bands). The weights \(w_{\mathbf{k}}\) and \(w(\epsilon_{i,\mathbf{k}})\) are the k-point weights and the Gaussian weights, respectively. The energy dependent Gaussian weight is computed as the integral of a normalized Gaussian centered at \(\epsilon_{i,\mathbf{k}}\) with standard deviation \(\sigma\) over the specified energy range, i.e.

\[\begin{split}\begin{align*} w(\epsilon_{i,\mathbf{k}}) &= \int_{E_{min}}^{E_{max}} \frac{1}{\sqrt{2 \pi \sigma^2}} \exp\left(\frac{\left(E - \epsilon_{i,\mathbf{k}}\right)^2}{2 \sigma^2}\right)\,dE \\ &= {\frac {1}{2}}\left[\mathrm{erf}\left(\frac{E_{max}-\epsilon_{i,\mathbf{k}}}{\sqrt{2} \sigma}\right)-\mathrm{erf}\left(\frac{E_{min}-\epsilon_{i,\mathbf{k}}}{\sqrt{2} \sigma}\right)\right], \end{align*}\end{split}\]

with \(\operatorname {erf}\) being the error function.

Note that in combination with the projector-augmented wave (PAW) method, \(n(\mathbf{r})\) is the smooth pseudo valence density, which does not necessarily integrate to the correct number of electrons. However, outside of the augmentation spheres, i.e. in the region where the bonds are located, \(n(\mathbf{r})\) does correspond to the all electron valence density.

Tersoff-Hamann approximation of simulated STM¶

The PartialElectronDensity can be used to simulate scanning tunneling microscopy (STM) images within the Tersoff-Hamann approximation [1] assuming that the tip wave function has s-orbital character. When simulating an STM image, the energy_range can be interpreted as the bias difference applied between the tip and the surface.

Notice that the Tersoff-Hamann approximation is the simplest approach to simulated STM images and it is not always applicable. Also bare in mind that far away from a surface the partial density will be zero for an LCAO calculation since the basis orbitals have a finite range.