ElectronDensity¶

class ElectronDensity(configuration, density_mesh_cutoff=None)¶

A class for calculating the electron density for a configuration.

Parameters:

configuration (MoleculeConfiguration | BulkConfiguration | DeviceConfiguration | SurfaceConfiguration) – The configuration for which the density should be calculated.
density_mesh_cutoff (PhysicalQuantity of type energy | GridSampling | OptimizedFFTGridSampling) – The mesh cutoff to be used to determine the density 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.

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.

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.

metatext()¶

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

nlprint(stream=None)¶

Print a string containing an ASCII table useful for plotting the AnalysisSpin object.

Parameters:: stream (python stream) – The stream the table should be written to.
Default: NLPrintLogger()

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.

spin()¶

Returns:: The spin the electron density is calculated for, always Spin.All.
Return type:: Spin.All

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¶

Calculate the electron density for a NH₃ molecule and save it to a file:

# Set up configuration
molecule_configuration = MoleculeConfiguration(
    elements=[Nitrogen, Hydrogen, Hydrogen, Hydrogen],
    cartesian_coordinates=[[ 0.     ,  0.      ,  0.124001],
                           [ 0.     ,  0.941173, -0.289336],
                           [ 0.81508, -0.470587, -0.289336],
                           [-0.81508, -0.470587, -0.289336]]*Angstrom
    )

# Define the calculator
calculator = LCAOCalculator()
molecule_configuration.setCalculator(calculator)

# Calculate and save the electron density together with the configration.
electron_density = ElectronDensity(molecule_configuration)
nlsave('results.hdf5', electron_density)
nlsave('results.hdf5', molecule_configuration)

nh3_density.py

Read in the electron density from a file and calculate the integrated density, which should equal the number of valence electrons in the calculation: \(\int n (x, y, z) dy dz dx = N_{elec}\). In this case this equals \(N_{elec}=3\cdot n_H + n_N=3 + 5=8\) since there are 5 (1) valence electrons included in the default PseudoDojo pseudo-potential for Nitrogen (Hydrogen).

# Import an electron density
density = nlread('results.hdf5', ElectronDensity)[0]

# Get the shape of the data.
shape = density.shape()

# Find the volume elements.
dX, dY, dZ = density.volumeElement()

# Calculate the unit area in the y-z plane.
dAYZ = numpy.cross(dY,dZ)[0]

# calculate density along x integrated over y,z
n_x = [ density[i,:,:].sum() * dAYZ for i in range(shape[0]) ]

# Integrate the density along x to get the number of electrons.
dx = dX.norm()
number_of_electrons = 0

for i in range(shape[0]):
    number_of_electrons += n_x[i]*dx

print('Total number of electrons = {:.5f}'.format(number_of_electrons))

 # Alternative to calculating the integral
dX, dY, dZ = density.volumeElement()

# Volume of a single grid point
dV = numpy.dot(numpy.cross(dX, dY), dZ)

density_unit = density.unit()
density_array = density.toArray() * density_unit
integral = density_array.sum() * dV

print('Alternative method:')
print('Total number of electrons= {:.5f}'.format(integral))

nh3_density_integral.py

Read in the electron density from a file and calculate the first moment of both the negatively charged electron density and the positively charged ions. The sum of the two contributions give the dipole moment for NH₃:

# Import an electron density and molecule configuration
density = nlread('results.hdf5', ElectronDensity)[0]
molecule_configuration = nlread('results.hdf5', MoleculeConfiguration)[0]

# Volume of a single grid point
dX, dY, dZ = density.volumeElement()
dV = numpy.dot(numpy.cross(dX, dY), dZ)

# Setup array for the first moment of the electron density.
electron_moment = numpy.zeros(3) * Ang * elementary_charge

# Loop over unit cell vectors:
for axis_index, axis in enumerate(['a', 'b', 'c']):
    (axis_values, density_values) = density.axisProjection(
        projection_type='sum',axis=axis, coordinate_type=Cartesian)
    # The minus sign accounts for the negative charge of electrons.
    moment = numpy.sum(axis_values * density_values) * dV * (-1*elementary_charge)
    electron_moment[axis_index] += moment

# Calculate the dipole moment for the positive nuclei
elements = molecule_configuration.elements()
cartesian_coordinates = molecule_configuration.cartesianCoordinates()

# Get the valence charge for each element from the basis set.
basis_set = molecule_configuration.calculator().basisSet()
ionic_charges = {b.element() : b.numberOfValenceElectrons() * elementary_charge for b in basis_set}

# Calculate the first moment of the positive ions:
ion_moment = numpy.zeros(3) * Ang * elementary_charge
for element, coordinate in zip(elements, cartesian_coordinates):
    ion_moment += ionic_charges[element] * coordinate

# Calculate the total dipole moment taking both negative electrons and positive ions into account.
total_dipole_moment = ion_moment + electron_moment

# Get the norm of the dipole moment in units of Debye (1 D = 0.2081943 * e*Ang)
dipole_moment = numpy.linalg.norm(total_dipole_moment.inUnitsOf(Angstrom*elementary_charge)) / 0.2081943

nlprint('Dipole moment = {:.3f} D'.format(dipole_moment))

For more examples on working with 3D grids, see HartreePotential.

Notes¶

This class inherits from the GridValues class.
Returns the electron density \(n ({\bf r})\), as defined in The Hartree potential and the electrostatic potential.
When PAW potentials are used, \(n(\mathbf{r})\) is the smooth pseudo valence density.