Part 2: Relaxation of the initial structure
In this tutorial, we use the experimental structure of bulk Hematite Fe\(_2\)O\(_3\) as a starting point to determine the lowestenergy antiferromagnetic spin state, using the PBE density functional and light settings.
We also determine the DFTPBE optimized lattice positions for later use as a template to create Fe\(_2\)O\(_3\) cluster structures, which will also be handled (in this tutorial) using the same density functional and settings.
Finally, we use the DFTPBE prerelaxed bulk structure and approximate Hessian matrix as a starting point to calculate the optimized structure, as well as energy band structure and density of states, using a hybrid density functional, HSE06. The hybrid density functional would be expected to approximate certain aspects of the system better than the semilocal density functional PBE, particularly the energy band gap and spin moment of the Fe ions.
All calculations will here be based on "light" settings for the sake of computational feasibility within a tutorial. In a calculation for production purposes, we would recommend "intermediate" settings for a bulk oxide.
Here is an additional remark related to the choice of FHIaims settings specifically for dense solids and "tight" settings. In principle, it would seem reasonable to postconverge results for a bulk oxides with "tight" settings. This could be done; however, it turns out that the default "tight" basis set especially for oxygen is already rather large, as is the "tight" basis set for Fe. In a relatively sparse structure and for relatively weak interactions (e.g., a Fe ion complexated by a molecule), "tight" settings might still be a reasonable choice. For a dense oxide, however, the density of basis functions per volume becomes very high, i.e., the basis set is close to the (over)converged limit in this case. Such calculations will tend to become fairly expensive. In earlier simulations of bulk oxides, we have found that "intermediate" settings at least for oxygen are a perfectly acceptable choice, at significantly cheaper cost than a "tight" basis set.
Step 1: PreRelaxation with PBE, light settings
Within the primitive unit cell of hematite Fe\(_2\)O\(_3\), two inequivalent antiferromagnetic spin configurations can be envisioned. We next perform computational structure optimizations (atomic positions and lattice vectors) for input spin configurations pertaining to either one of them, aiming to find the lowerenergy spin configuration among them.
For either configuration, we repeat below only the relevant sections of the geometry.in
files pertaining to the spin state of the Fe atoms (see Part 1 for the full geometry.in
file derived from an experimental crystal structure):

upupdowndown (uudd) configuration:
... atom_frac 0.3549800000000001 0.3549800000000000 0.3549800000000000 Fe initial_moment 4. atom_frac 0.1450199999999998 0.1450199999999998 0.1450199999999998 Fe initial_moment 4. atom_frac 0.6450199999999998 0.6450200000000000 0.6450199999999996 Fe initial_moment 4. atom_frac 0.8549800000000003 0.8549799999999999 0.8549799999999997 Fe initial_moment 4. ...

downupupdown (duud) configuration:
... atom_frac 0.1450200000000000 0.1450200000000000 0.1450200000000000 Fe initial_moment 4. atom_frac 0.3549800000000000 0.3549800000000000 0.3549800000000000 Fe initial_moment 4. atom_frac 0.6450200000000000 0.6450200000000000 0.6450200000000000 Fe initial_moment 4. atom_frac 0.8549800000000000 0.8549800000000000 0.8549800000000000 Fe initial_moment 4. ...
Update:
Instead of the above input files (kept for consistency with the solutions currently included with the tutorial), it would be much better to initialize the Fe atoms with the physically expected spin and charge states of Fe\(^{3+}\) ions in their expected highspin state:
...
atom_frac 0.1450200000000000 0.1450200000000000 0.1450200000000000 Fe
initial_charge 3.
initial_moment 5.
atom_frac 0.3549800000000000 0.3549800000000000 0.3549800000000000 Fe
initial_charge 3.
initial_moment 5.
atom_frac 0.6450200000000000 0.6450200000000000 0.6450200000000000 Fe
initial_charge 3.
initial_moment 5.
atom_frac 0.8549800000000000 0.8549800000000000 0.8549800000000000 Fe
initial_charge 3.
initial_moment 5.
...
and, correspondingly, set initial charge 2.
keywords for each O atom in the structure. This initialization, which
corresponds to the physically expected state of the hematite structure, is much closer to the final selfconsistent
result and therefore leads to much faster convergence of the s.c.f. cycle.
Tasks:
 Prerelax the two configurations with light species_default settings. (Refer to Tutorial Basics of Running FHIaims regarding how to run the FHIaims code for given
geometry.in
andcontrol.in
files.)
You will need to prepare a control.in
file for this purpose. There are several options to do this  see Tutorial Basics of Running FHIaims for the meaning of some of the basic keywords.
Using a text editor, you could create a control.in file including the following keywords:
xc pbe
spin collinear
relativistic atomic_zora scalar
k_grid 10 10 10
relax_geometry trm 5e3
relax_unit_cell full
Additionally, append the "light" species defaults for O and for Fe as described in Tutorial Basics of Running FHIaims.
There are two important keywords here:
First, we do need to request specifically that the unit cell vectors be included as variable parameters in the totalenergy minimization with respect to the atomic positions, i.e.
relax_unit_cell full
Second, any periodic calculation will require a kspace grid to discretize reciprocal space and facilitate integrals in reciprocal space. The chosen discretization is
k_grid 10 10 10
meaning that the unit cell of reciprocal space will be discretized using a threedimensional grid of kvalues that is parallel to the three reciprocalspace lattice vectors, with ten grid spacings each along each reciprocal lattice vector.
How do we know that this discretization is sufficiently dense? In principle, we should test the kspace discretization for numerical convergence. As a starting point (used here), a useful rule of thumb is that the kspace discretization value n_i corresponding to a particular realspace unit vector a_i should fulfill an approximate criterion
All three realspace unit vectors of the primitive unit cell of Hematite Fe\(_2\)O\(_3\) have an approximate length of 5.42 Å, i.e., our chosen kspace grid fulfills this rule of thumb.
In a publishable simulation, we would likely want to use a kspace grid with better convergence for final results, or at least a kspace grid that is better tested.
We could, very similarly, also construct a control.in
file using GIMS. This choice has the additional advantage that the Brillouin zone of the material is visualized during the construction of the input files.
We could, finally, use the clims commandline tool climspreparerun relax
to obtain a template for the control.in file.

Analyze the results using GIMS. Compare the final total energies of both configurations. Using "light" settings, at least, the 'duud' configuration championed in Ref. https://doi.org/10.3390/min4010089 is lower in energy by a significant amount (0.85 eV per unit cell). It would be interesting to confirm this result by postrelaxing the respective structures with "intermediate" or better settings and a better density functional (e.g., HSE06).

One thing that we are missing is a deeper analysis of the final, selfconsistent spin state of the system, in an atomresolved manner. This state does not necessarily have to be the same as the initially configured spin state.
A look into the aims.out
file would show that the overall spin moment of the system is exactly zero, not unexpected for an antiferromagnetic system.
A better measure for the overall spin state is the integral over the modulus of the difference of the spin densities over the entire unit cell, that is
This quantity is printed in every s.c.f. iteration of a spinpolarized FHIaims calculation. To obtain the final value from the converged FHIaims standard output file at the command line, type
grep rho_uprho_dn aims.out  tail 1
This command will reveal the overall integrated polarization of the system.
A qualitative atomresolved value for the spin configuration of the Fe ions could also be obtained using a Mulliken analysis, i.e., a projection of atomic states and densities onto individual atomiclike components, using basis functions as a projector.
A Mulliken analysis, or better, analysis in terms of a projected density of states (which we did not include here), could be activated using the keyword
output mulliken
in control.in
.
The analysis of spin state and of the projected density of states is provided for a different, likely better density functional (HSE06) in Part 3.
The corresponding output is there obtained as a byproduct of a calculation of the projected density of states. For a "light" DFTPBE prerelaxation, a Mulliken analysis would still be useful as a consistency check and diagnostic output, but it would have to be verified by a more precise higherlevel calculation.
Step 2: FinalRelaxation with HSE06, light settings
Semilocal density functionals (such as PBE) can have severe problems in transition metal oxides, particularly in cases in which the individual transition metal ions are not formally in a highspin state. We skirted this problem with our choice of Fe\(_2\)O\(_3\), which places the Fe(III) ion in its highest possible spin state in the ground state.
Nevertheless, for more precise calculations (including band structures and densities of states), a hybrid density functional would be the better choice.
We here use the HSE06 density functional to follow up on our DFTPBE prerelaxation.
FHIaims actually does include a fairly capable, linearscaling hybrid density functional theory implementation for nonperiodic and periodic systems. This implementation has been used to compute spinorbit coupled band structures (intermediate settings) for periodic structure sizes above 1,500 atoms including heavy elements, without making any significant approximations to the HSE06 functional or other settings along the way.
Nevertheless, hybrid DFT is still significantly more expensive than semilocal DFT, both in CPU time and in memory consumption. This increase in resources is particularly significant for dense solids, in which a high density of basis functions per volume exists and overlaps. For Fe\(_2\)O\(_3\), this creates somewhat of a computational challenge.
For the purposes of this tutorial, we therefore use "light" settings for the DFTHSE06 postrelaxation. The output file available in the solutions
folder shows that this calculation (full relaxation including unit cell) completes in a somewhat reasonable wall clock time (just over half an hour) on a reasonably large computer (here, 216 current Intel CPU cores).
Usually, we would want to perform this simulation with "intermediate" settings.
During this tutorial, you may want to try the calculation or you may simply wish to look at the existing output file. The steps to perform the simulation are:
 Create a new folder for the HSE06 relaxation. Copy the
geometry.in.next_step
to thegeometry.in
file in the new folder and initialize the atom with the correctinitial_moment
(they are not in thegeometry.in.next_step
file as there is no exact atomic spin decomposition possible with the converged density). Also, copy over the earlierhessian.aims
file.  Prepare the input file
control.in
for the final relaxation with HSE06. Thecontrol.in
should contain the following, important lines to change the density functional to HSE06 with a screening parameter of 0.11 (Bohr radii)\(^{1}\) and its (default) exchange mixing parameter of 0.25:xc hse06 0.11 hse_unit bohr hybrid_xc_coeff 0.25
Alternatively, the modification can be accomplished using GIMS.
Finally, you could use (again) climspreparerun relax hse06
to create the template for the control.in
file.

Please also (manually) add the following line to the
control.in
file:elsi_restart read_and_write 100
This keyword will write the converged density matrices for all kpoints to disk (unfortunately, a large number of such files for a dense solid).
The parameter 100
implies that the density matrices will only be written after every 100 s.c.f. cycles or for the fully s.c.f.converged density matrix. We can later reuse this converged density matrix to restart the calculation for the converged atomic positions and save some time for postprocessing tasks, i.e., to generate output such as energy band structures, full or partial densities of states, etc.
All solutions to this tutorial can be again found in the solution
folder. Have fun!