# Structure relaxation tutorial

## Requirements

### Software components

- ASE
- nanotools
- RESCU+

### Pseudopotentials

I will need the following pseudopotentials.

- Si_AtomicData.mat

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

### References

- [ASE website](https://wiki.fysik.dtu.dk/ase/index.html)
- [tutorial_etot.md](../etot/tutorial_etot.md)

## Briefing

In this tutorial, I explain how to find the equilibrium positions of the atoms of a silicon unit cell using the `Relax` class.

## Code

```python
from ase.build import bulk
from nanotools import System, TotalEnergy, Relax

a = 5.43 # lattice constant in ang
atoms = bulk("Si", "diamond", a=a, cubic=True)
atoms.rattle(stdev=0.05, seed=1)
atoms.positions[0,:] = 0.

sys = System.from_ase_atoms(atoms)
sys.cell.set_resolution(0.2)
sys.kpoint.set_grid([5,5,5])
ecalc = TotalEnergy(sys)
ecalc.solver.mix.alpha = 0.5
ecalc.solver.set_mpi_command("mpiexec --bind-to core -n 16")
rlx = Relax.from_totalenergy(ecalc)
rlx.solve()
```

## Explanations

Here is a high level view of the calculation work flow:

1. [Build the system](#build-the-system) using ASE.
2. Create a [relax calculator](#create-a-relax-calculator).
3. [Relax the structure](#relax-the-structure) using RESCU+ as a force calculator and ASE as an optimization engine.

### Build the system

Let's first load some Python modules and class definitions, and then build a silicon primitive cell with lattice parameter of 5.43 Ang.

```python
from ase.build import bulk
from nanotools import System, TotalEnergy, Relax

a = 5.43 # lattice constant in ang
atoms = bulk("Si", "diamond", a=a, cubic=True)
atoms.rattle(stdev=0.05, seed=1)
atoms.positions[0,:] = 0.
```

I used the `bulk` function imported from the module `ase.build`.
This function is quite handy to build various basic crystals.
The `Atoms` class has a method `rattle` which moves the atoms in random directions.

### Create a relax calculator

The following step is creating and linking a `Relax` calculator.
If you haven't done so already, you should go through the total energy tutorial [tutorial_etot.md](../etot/tutorial_etot.md).
Next, proceed as follows.

```python
sys = System.from_ase_atoms(atoms)
sys.cell.set_resolution(0.2)
sys.kpoint.set_grid([5,5,5])
ecalc = TotalEnergy(sys)
ecalc.solver.mix.alpha = 0.5
ecalc.solver.set_mpi_command("mpiexec --bind-to core -n 16")
rlx = Relax.from_totalenergy(ecalc)
```

I first create a reference `TotalEnergy` calculator.
I then pass it to `Relax.from_totalenergy` to create the `Relax` calculator.

### Relax the structure

Finally, I call `solve` which will internally use RESCU for the forces and stresses, and ASE's relaxation engine.

```python
rlx.solve()
```

During the relaxation, the results are saved to various `nano_rlx` files (e.g. the energy and max force in `nano_rlx_log.out`).
The last results are saved in `nano_rlx_out.json` and `nano_rlx_out.h5`.
One can then look at the equilibrium positions as follows

```python
from nanotools import TotalEnergy
import numpy as np
ecalc = TotalEnergy.read("nano_rlx_out.json")
fxyz = ecalc.system.atoms.fractional_positions
print(np.round(fxyz-fxyz[0,:],4))
```

and recognized the expected atomic positions.