Usage Examples

Calculate the electron difference 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 = HuckelCalculator()
molecule_configuration.setCalculator(calculator)

# Calculate and save the electron difference density
difference_density = ElectronDifferenceDensity(molecule_configuration)
nlsave('results.nc', difference_density)

nh3_density.py

Read in the electron difference density from a file and print

# import an electron difference density
difference_density = nlread('results.nc', ElectronDifferenceDensity)[0]
# Get the shape of the data.
shape = difference_density.shape()
# Find the volume elements.
dX, dY, dZ = difference_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 = [ difference_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]):
    print dx*i, n_x[i]
    sum+=n_x[i]*dx
print 'total density=', sum

nh3_density_integral.py

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

# import an electron difference density
difference_density = nlread('results.nc', ElectronDifferenceDensity)[0]

# Get the shape of the data.
ni, nj, nk = difference_density.shape()
# Find the volume elements.
dX, dY, dZ = difference_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 difference density.
grid_data = difference_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) = "
print m2

nh3_density_multipoles.py

For more examples on how to average and handle 3-d grids, see also the examples for ElectrostaticDifferencePotential.

Notes