# 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.
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 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. 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 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. 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 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. 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. The Cartesian coordinate of the given grid index. PhysicalQuantity of type length.
metatext()
Returns: The metatext of the object or None if no metatext is present. 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. 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. tuple of three int.
spin()
Returns: The spin the electron density is calculated for, always Spin.All. 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 A new GridValues object for the specified spin. GridValues
toArray()
Returns: The values of the grid as a numpy array slicing off any units. numpy.array
unit()
Returns: The unit of the data in the grid. A physical unit.
unitCell()
Returns: The unit cell of the grid. PhysicalQuantity of type length.
volumeElement()
Returns: The volume element of the grid represented by three vectors. 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

# 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

# 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)

nh3_density_multipoles.py