Mulliken population tutorial

Requirements

Software components

  • nanotools

  • RESCU+

Pseudopotentials

I will need the following pseudopotentials.

  • Al_AtomicData.mat

  • P_AtomicData.mat

Let’s copy them from the pseudo database to the current working directory and export RESCUPLUS_PSEUDO=$PWD.

References

  • F. Matthias Bickelhaupt, Nicolaas J. R. van Eikema Hommes, Célia Fonseca Guerra, and Baerends, E. J. (1996). The Carbon−Lithium Electron Pair Bond in (CH3Li)n (n = 1, 2, 4). Organometallics, 15(13), 2923–2931.

Briefing

In this tutorial, I show how to calculate the Mulliken populations of AlP.

The Mulliken population matrix \(\mathbf {P_{\mu\nu}}\) is defined as

\[ \mathbf {P_{\mu\nu}} = \mathbf {D_{\mu\nu}} \mathbf {S_{\mu\nu}} \]

where the density matrix \(\mathbf {D_{\mu\nu}}\) is defined as

\[ \mathbf {D_{\mu\nu}} = \mathbf{2}\sum_{i} \mathbf {C_{\mu i}} \mathbf {C_{\nu i}^*} \]

and \(\mathbf {C_{\mu i}}\) is the matrix of atomic orbital coefficients. The above equation is generalized to account for k-sampling and spin in our implementation. The Mulliken orbital populations are defined as

\[ p_{\mu} = P_{\mu\mu} + \sum\limits_{\nu \neq \mu}\frac{P_{\mu\mu}}{(P_{\mu\mu} + P_{\mu\nu})}\mathbf {P_{\mu\nu}} \]

and the Mulliken atomic populations \(p_A\) are the sum over all orbitals \(\mu\) belonging to atom A.

Setup

Create the Python script etot.py with the following content (see the tutorial on TotalEnergy)

from nanotools import Atoms, Cell, System, TotalEnergy
import numpy as np
a = 5.4510 # lattice constant (ang)
avec = a * np.array([[0.0,0.5,0.5],[0.5,0.0,0.5],[0.5,0.5,0.0]])
cell = Cell(avec=avec, resolution=0.14)
fxyz = [[0.00,0.00,0.00],[0.25,0.25,0.25]]
atoms = Atoms(fractional_positions=fxyz, formula="AlP")
sys = System(cell=cell, atoms=atoms)
sys.kpoint.set_grid([4,4,4])
calc = TotalEnergy(sys)
calc.solve()

and mul.py

from nanotools import Mulliken
import numpy as np
calc = Mulliken.from_totalenergy("nano_scf_out.json")
calc.solve()

Explanations

Here is a high level view of the calculation workflow:

  1. Initialize the Mulliken calculator from a ground state calculation (TotalEnergy).

  2. Perform a Mulliken population calculation (non-self-consistently).

  3. Analyze the population data.

Initialize the Mulliken calculator

Upon completing the total energy calculation, the results are saved in nano_scf_out.json. I will initialize a Mulliken calculator using the class method from_totalenergy.

calc = Mulliken.from_totalenergy("nano_scf_out.json")
calc.solve()

The method from_totalenergy will copy the system information and initialize the Mulliken calculator. I could increase the k-sampling resetting the kpoint grid with set_grid but I’ll maintain it here.

Mulliken population calculation

RESCU+’s solvers are invoked calling the solve method

calc.solve()

The method writes all parameters to a JSON file, then calls the relevant (Fortran) program, then loads the data back into the calculator. The output of rescuplus goes to nano_mul_out.h5 and nano_mul_out.json. The population matrix can grow quite large, and hence it is stored in nano_mul_out.h5 under /mulliken/total. In the collinear spin framework, the population matrices are stored in /mulliken/spin-up, /mulliken/spin-down and in the non-collinear spin framework under /mulliken/total, /mulliken/spin-x, /mulliken/spin-y, /mulliken/spin-z.

Analyze the population data

The Mulliken population data is stored in the calc.mulliken after completing the calculation. For convenience, the Mulliken atomic and orbital populations are calculated and returned in calc.mulliken.atom_populations and calc.mulliken.orbital_populations respectively. In the present example, I can print the Mulliken population of each atoms as follows

print(calc.mulliken.atom_populations)

and get \([[1.9413278], [6.0586722]]\), which means P is grabbing an electron from Al according to this tool, more or less. If spin is included, there is one column per spin channel, corresponding with the HDF5 paths above.

To further analyze the Mulliken population data, use the following attributes of calc.mulliken:

  • orbital_atom: atomic index;

  • orbital_index: orbital index;

  • orbital_l: orbital angular momentum;

  • orbital_m: z-axis angular momentum;

  • orbital_species: species index.