ElectronDensity

class ElectronDensity(configuration)

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

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

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 | unicode | 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 | 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.
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
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 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
electron_density = ElectronDensity(molecule_configuration)
nlsave('results.nc', electron_density)

nh3_density.py

Read in the electron density from a file and print \(\int n (x, y, z) dy dz dx = 0\):

# 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.linalg.norm( numpy.cross(dY,dZ) )
# calculate density along x integrated over y,z
n_x = [ density[i,:,:].sum() * dAYZ for i in range(shape[0]) ]

#print out the result
dx = dX.norm()

sum = 0 * Units.Bohr**-2

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

print('Total electron density=', sum)

nh3_density_integral.py

Read in the electron density from a file and calculate multipoles of the density:

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

# Get the shape of the data.
ni, nj, nk = density.shape()
# Find the volume elements.
dX, dY, dZ = density.volumeElement()
length_unit = dX.unit()
# Calculate the volume of the volume element.
dV = numpy.dot(dX, numpy.cross(dY,dZ)) * length_unit**3
# Get a reference for the data of the density.
grid_data = density[:,:,:]*Units.e

# Calculate the absoute density sum.
density_abs_sum = abs(grid_data).sum()

# Calculate center of mass
Mi = sum([ i*abs(grid_data[i,:,:]).sum() for i in range(ni)])/density_abs_sum
Mj = sum([ j*abs(grid_data[:,j,:]).sum() for j in range(nj)])/density_abs_sum
Mk = sum([ k*abs(grid_data[:,:,k]).sum() for k in range(nk)])/density_abs_sum
center_of_mass = Mi*dX + Mj*dY+ Mk*dZ

# Calculate monopole
density_sum = grid_data.sum()
m0 = density_sum*dV

#calculate dipole
dens_i = sum([ (i-Mi)*grid_data[i,:,:].sum() for i in range(ni)])
dens_j = sum([ (j-Mj)*grid_data[:,j,:].sum() for j in range(nj)])
dens_k = sum([ (k-Mk)*grid_data[:,:,k].sum() for k in range(nk)])
m1 = (dX*dens_i + dY*dens_j + dZ*dens_k)*dV

# Calculate quadropoles (3 x_i *x_j-r*r delta_ij)*dens
dens_ii = sum([ (i-Mi)**2 * grid_data[i,:,:].sum() for i in range(ni)])
dens_jj = sum([ (j-Mj)**2 * grid_data[:,j,:].sum() for j in range(nj)])
dens_kk = sum([ (k-Mk)**2 * grid_data[:,:,k].sum() for k in range(nk)])
dens_ij = sum([ (i-Mi)*(j-Mj)*grid_data[i,j,:].sum()
	          for i in range(ni) for j in range(nj)])
dens_ik = sum([ (i-Mi)*(k-Mk)*grid_data[i,:,k].sum()
	          for i in range(ni) for k in range(nk)])
dens_jk = sum([ (j-Mj)*(k-Mk)*grid_data[:,j,k].sum()
	          for j in range(nj) for k in range(nk)])

m2 = (dX.outer(dX)*dens_ii + dY.outer(dY)*dens_jj + dZ.outer(dZ)*dens_kk\
   + (dX.outer(dY)+dY.outer(dX))*dens_ij + (dX.outer(dZ)+dZ.outer(dX))*dens_ik\
   + (dY.outer(dZ)+dZ.outer(dY))*dens_jk) * dV

r2 = sum(m2)
m2 = 3*m2 - numpy.identity(3)*r2

print("center of mass (bohr)  = ", center_of_mass)

print("monopole   (e)         = ", m0)
print("dipole     (e*bohr)    = ", m1)
print("quadropole (e*bohr**2) = ", m2)

nh3_density_multipoles.py

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

Notes