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

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

Note about syntax

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