# 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.](https://doi.org/10.1021/OM950966X)

## 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](../etot/tutorial_etot.md))

```python
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`

```python
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](#initialize-the-mulliken-calculator) from a ground state calculation (`TotalEnergy`).
2. Perform a [Mulliken population calculation](#mulliken-population-calculation) (non-self-consistently).
3. [Analyze the population data](#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`.

```python
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

```python
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

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