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

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

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 using ASE.

  2. Create a relax calculator.

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

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. Next, proceed as follows.

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.

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

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.