# Total energy tutorial

## Requirements

• 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

```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 and specify calculation parameters like resolution, k-sampling, tolerance, etc.

2. Create a `TotalEnergy` calculator. In the process, all unspecified parameters will be initialized.

3. Solve the Kohn-Sham equation self-consistently calling `solve`.

4. 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

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

```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")
```

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

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

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

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

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

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