PartialElectronDensity

class PartialElectronDensity(configuration, kpoints=None, band_indices=None, energy_range=None, spectrum_method=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 calculate the eigenvalues. This can either be a MonkhorstPackGrid, RegularKpointGrid, or a list of 3-dimensional fractional k-points e.g. [[0.0, 0.0, 0.0], [0.0, 0.0, 0.1], ...].
    Default: MonkhorstPackGrid(nx,ny,nz), where nx, ny, nz is the sampling 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)
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 Cartesian 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.

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

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 | unicode | 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 | unicode | 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
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

  1. 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.
  2. 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.
  3. kpoints. The partial electron density is calculated from the eigenstates at the specified k-points.
  4. 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.

Tersoff-Hamann approximation of simulated STM

The PartialElectronDensity can be used to simulate scanning tunneling microscopy (STM) images within the Tersoff-Hamann approximation [TH85] 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.

[TH85]J. Tersoff and D. R. Hamann. Theory of the scanning tunneling microscopy. Phys. Rev. B, 31:805, 1985. doi:10.1103/PhysRevB.31.805.