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 lowest-energy antiferromagnetic spin state, using the PBE density functional and light settings.

We also determine the DFT-PBE 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 DFT-PBE pre-relaxed 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 FHI-aims settings specifically for dense solids and "tight" settings. In principle, it would seem reasonable to post-converge 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: Pre-Relaxation 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 lower-energy 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):

up-up-down-down (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.
...

down-up-up-down (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 high-spin 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 self-consistent result and therefore leads to much faster convergence of the s.c.f. cycle.

Tasks:

Pre-relax the two configurations with light species_default settings. (Refer to Tutorial Basics of Running FHI-aims regarding how to run the FHI-aims code for given geometry.in and control.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 FHI-aims 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 5e-3
relax_unit_cell                    full

Additionally, append the "light" species defaults for O and for Fe as described in Tutorial Basics of Running FHI-aims.

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 total-energy minimization with respect to the atomic positions, i.e.

relax_unit_cell                    full

Second, any periodic calculation will require a k-space 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 three-dimensional grid of k-values that is parallel to the three reciprocal-space 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 k-space discretization for numerical convergence. As a starting point (used here), a useful rule of thumb is that the k-space discretization value n_i corresponding to a particular real-space unit vector a_i should fulfill an approximate criterion

\[ n_i * a_i > 50 Å \]

All three real-space unit vectors of the primitive unit cell of Hematite Fe\(_2\)O\(_3\) have an approximate length of 5.42 Å, i.e., our chosen k-space grid fulfills this rule of thumb.

In a publishable simulation, we would likely want to use a k-space grid with better convergence for final results, or at least a k-space 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 command-line tool clims-prepare-run --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 post-relaxing 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, self-consistent spin state of the system, in an atom-resolved 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

\[ \int_{uc} d^3r |n^\uparrow(r) - n^\downarrow(r)| \]

This quantity is printed in every s.c.f. iteration of a spin-polarized FHI-aims calculation. To obtain the final value from the converged FHI-aims standard output file at the command line, type

grep rho_up-rho_dn aims.out | tail -1

This command will reveal the overall integrated polarization of the system.

A qualitative atom-resolved 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 atomic-like 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" DFT-PBE pre-relaxation, 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 higher-level calculation.

Step 2: Final-Relaxation 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 high-spin 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 DFT-PBE pre-relaxation.

FHI-aims actually does include a fairly capable, linear-scaling hybrid density functional theory implementation for non-periodic and periodic systems. This implementation has been used to compute spin-orbit 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 DFT-HSE06 post-relaxation. 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 the geometry.in file in the new folder and initialize the atom with the correct initial_moment (they are not in the geometry.in.next_step file as there is no exact atomic spin decomposition possible with the converged density). Also, copy over the earlier hessian.aims file.
Prepare the input file control.in for the final relaxation with HSE06. The control.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) clims-prepare-run --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 k-points 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!