Part 5: Converging clusters

We next attempt a full geometry optimization of the initial cluster structure generated in Part 4, aiming to find a local minimum of the potential energy surface for a given density functional (here, PBE).

We refer again to our discussion of initial structure, charge and spin state in Tutorial Basics of Running FHI-aims. A given structure optimization might eventually converge to a particular geometry and electronic and/or spin state. However, there is no guarantee whatsoever that such an outcome will be the only possible state (it won't) or even a state of the system that is in any way physically attainable. A more refined set of hypotheses for possible cluster structures and spin states would gave to be considered in order to explore whether a given result is physically meaningful.

Naive Electronic Structure Calculation

We have an initial structure, charge and spin state from Part 4 of this tutorial. While this structure has not been deeply vetted for physical consistency, we might as well go ahead and try to optimize its structure with a given density functional. The control.in file for this non-periodic structure optimization using "light" settings looks surprisingly simple, easily constructed using GIMS or following the same steps as described in Tutorial Basics of Running FHI-aims.

What could possibly go wrong?

If this were a simple, non-spinpolarized structure, close to the chemical equilibrium of a given material, not much would normally go wrong. FHI-aims performs the task proposed here routinely and with competitive efficiency for molecules, semiconductors, organic-inorganic hybrid materials of significant complexity etc., up to structure sizes of a few thousand atoms.

However, we are here dealing with an initial structure that is not chemically reasonable (it has near-surface atoms in awkward positions, having lost their bonding partners). Additionally, we have a potentially complex spin state to sort out.

Just trying to relax this structure without thinking much is not a practical way forward. An example can be found in solution/relaxation-naive-try.

The s.c.f. convergence of this structure is painfully slow with default settings (which are, themselves, fairly conservative). We might be able to find a way to make it work and eventually arrive at a chemically more reasonable, local minimum-energy geometry.

However, at this stage, it is preferable to do some quick, chemically motivated assessment of the problem. Looking at the input structure, it turns out that the surface has some oxygen atoms that have lost their bonding partners. They look, loosely speaking, "singly bonded". In any chemically reasonable scenario, the highly reactive oxygen atom would not be in such a place. In a solution, those atoms would at least be saturated by H atoms or by a similar bonding partner. However, on their own, these O atoms are very high in energy and chemically unhappy.

It is a fundamental fact of electronic structure theory that a structure that is chemically unreasonable often does not lead to efficient or sensible convergence of the s.c.f. cycle.

Rather than battling an electronic structure code, however well-intentioned, to find a way to do non-linear optimization for a badly converging problem, the appropriate path forward is usually to pursue a better problem definition (in scientific terms, a better hypothesis for the structure under consideration).

It follows that we might benefit from a more "chemically motivated" strategy to guide our cluster towards a chemically more sensible configuration. In the following, we pursue two such approaches.

1. Hydrogen saturation

Keep the cluster size - however, pursue a two-step strategy (here, manual and somewhat crude) to guide the system towards a better structure:

First, saturate all singly bonded O atoms with H to place them in a chemically more acceptable, bonded state. Relax this structure.
Second, after the initial structure optimization, remove the H atoms. This re-creates some singly bonded O atoms. However, the remaining atoms are now in more acceptable positions. Take the remaining singly bonded O atoms, move them individually (by editing the structure in GIMS) to a more plausible position close by, and post-relax the resulting structure.

This process is documented in solutions/relaxation-H-saturated

The structure optimization with H-saturated O atoms is a very long trajectory, but it finds a local minimum-energy structure without significant problems, as documented in solutions/relaxation-H-saturated/step-1-H-saturated.

The cluster looks quite different after this first step, indicating that this cluster was too small and does not remain stable in a hematite-like structure.

The process can be visualized in GIMS.

However, this is quite a lengthy structure optimization and so some more targeted analysis may also be in order. As an alternative to GIMS, the script call

[FHI-aims-directory]/utilities/create_geometry_zip.pl aims.out

will create an archive geometries.zip that contains all intermediate geometries created during the structure optimization, in order. This file geometries.zip can be opened with the Jmol viewing program. Specifically, and for the example of the command line of the authors' Macbook:

/Users/blum/software/jmol-14.19.1/jmol.sh geometries.zip

This viewer enables some further flexibility in visual analysis of the structure and relaxation in question.

Next, in GIMS, remove the Hydrogen atoms and move the remaining, single-bonded oxygen atoms to plausible positions. For our example, we edited atom 8 and atom 68::

# atom 8
atom 2.56675 4.53067 -3.70961 O
->
atom 0.40518 2.24984 -5.51587 O

and

# atom 68
atom 0.76001 -4.07773 -5.79183 O
->
atom 0.54959 -4.07795 -3.58703 O

Begin a new calculation using this adjusted, H-free geometry (cf. solutions/relaxation-H-saturated/step-2-H-removed).

We finally end up with a stoichometric, 80-atom Fe\(_2\)O\(_3\) cluster in a local minimum-energy configuration.

2. Removing single-bonded atoms

We removed the single-bonded oxygen atoms from the initial structure, but we kept the stoichiometry of the cluster right (i.e., we also removed additionally O and Fe atoms to achieve a multiple of the Fe2O3 stoichiometry). See solutions/no-single-bonds for the result.

This is no longer an 80-atom cluster but it does yield one more (out of many possible) hypothetical Fe\(_2\)O\(_3\) nanostructures as model systems to simulate catalytic processes.

3. Other options

In a solvated environment and using physical principles, we could have created many other structure variants of Fe\(_2\)O\(_3\) clusters. One obvious pathway would be to assess the surface energies of different surface orientations and terminations as a function of variables such as the experimental pH value and use a Wulff-type minimum surface free energy construction to obtain plausible cluster shapes.

In all likelihood, such geometries would include H- or OH-terminations of some kind. Also, an implicit or explicit solvation environment might help.

We have not pursued these routes here since they would lead well beyond this already lengthy tutorial. Nevertheless, in a more realistic simulations, the steps above might be one way to formulate a more sensible hypothesis than we have done.