# Total energy tutorial

## Requirements

### Software components

- nanotools
- RESCU+

### Pseudopotentials

I will need the following pseudopotentials.

- Ga_AtomicData.mat
- As_AtomicData.mat

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

## Briefing

In this tutorial, I show how to calculate the total energy of the GaAs crystal.

## Code

```python
from nanotools import Atoms, Cell, System, TotalEnergy
a = 2.818
cell = Cell(avec=[[0.,a,a],[a,0.,a],[a,a,0.]], resolution=0.12)
fxyz = [[0.00,0.00,0.00],[0.25,0.25,0.25]]
atoms = Atoms(fractional_positions=fxyz, formula="GaAs")
sys = System(cell=cell, atoms=atoms)
calc = TotalEnergy(sys)
calc.solve()
print(calc.energy.etot)
```

## Explanations

Here is a high level view of the calculation workflow:

1. [Define a system](#system-definition) and specify calculation parameters like resolution, k-sampling, tolerance, etc.
2. [Create a `TotalEnergy` calculator](#total-energy-calculator). In the process, all unspecified parameters will be initialized.
3. [Solve the Kohn-Sham equation](#solve-the-kohn-sham-equation) self-consistently calling `solve`.
4. [Analyze the data](#analyze-the-data) and look at the results.

### System definition

I first need to specify the system configuration and parameters, which is done creating a `System` object.
All I need is specify what the lattice is, which means

- what is the shape of the cell;
- what are the species involved;
- what are the atoms' positions.

The cell can be specified with a `Cell` object, and the atomic species and coordinates with an `Atoms` object as follows

```python
cell = Cell(avec=[[0.,a,a],[a,0.,a],[a,a,0.]], resolution=0.12)
fxyz = [[0.00,0.00,0.00],[0.25,0.25,0.25]]
atoms = Atoms(fractional_positions=fxyz, formula="GaAs")
```

The `avec` attribute holds the lattice vectors as three row-vectors.
Here we set the `resolution` to 0.12 Ang, the default length unit.
The fractional coordinates are passed via `fractional_positions` for the gallium and arsenic atoms respectively, as specified by `formula`.

### Total energy calculator

Next, we create a `System` object passing `cell` and `atoms` to the class constructor; and we initialize a `TotalEnergy` calculator with `sys`.
We could then adjust additional parameters, for instance

```python
sys = System(cell=cell, atoms=atoms)
calc = TotalEnergy(sys)
calc.system.kpoint.set_grid([5,5,5])
calc.solver.set_mpi_command("mpiexec -n 16")
```

#### Note about syntax

If one prefers a more traditional, _flat text_ input style, the above can be written as follows

```python
from nanotools import TotalEnergy
a = 2.818 # lattice constant (ang)
inp = {"system": {"atoms": {}, "cell": {}}}
inp["system"]["cell"]["avec"] = [[0.,a,a],[a,0.,a],[a,a,0.]]
inp["system"]["cell"]["resolution"] = 0.12
inp["system"]["atoms"]["formula"] = "GaAs"
inp["system"]["atoms"]["fractional_positions"] = [[0.00,0.00,0.00],[0.25,0.25,0.25]]
calc = TotalEnergy(**inp)
```

The `**` operator unpacks the `dict` `inp` into keyword arguments.
We advise against it at the beginning since one looses the contextual documentation of a Python shell.

### Solve the Kohn-Sham equation

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_scf_out.h5` and `nano_scf_out.json`, the former stores large variables like matrices (e.g. the Hamiltonian) and discretized fields (e.g. density, wavefunctions) while the later stores the rest.
`nano_scf_out.json` reproduces all parameters of the input file, but calculated data like the total energy is updated.
As suggested by the extensions, the format are HDF5 and JSON respectively.
You may thus parse and read the data using your favorite HDF5 and JSON toolboxes.
If for some reason you would like to move the output to `gaas_scf_out` (with both extensions), you may call

```python
calc.solve(input="gaas_scf_in", output="gaas_scf_out")
```

Important data is stored in

- `calc.energy.eigenvalues`, the Kohn-Sham eigenvalues.
- `calc.energy.etot`, the total energy.
- `calc.energy.forces`, the atomic forces.
- `calc.energy.stress`, the unit cell stress tensor.

### Analyze the data

As mentioned above, `calc` will load the results after `rescuplus` is finished.
The total energy is printed typing

```python
print(calc.energy.etot)
-2724.9067779871266 electron_volt
```

Note that if I look at `nano_scf_out.json`, the total energy is in atomic units, but after constructing the object the total energy is in eV.
This is because `TotalEnergy` will convert the units after reading `nano_scf_out.json`.
You can always change it back to atomic units as follows

```python
calc.energy.set_units("atomic")
print(calc.energy.etot)
-100.13847719270561 hartree
```