# Part 1: Spin-unpolarized Simulation

Now, let's put our hands on FHI-aims. This tutorial takes the H$$_2$$O molecule as an example for a structure optimization with "light" settings. The starting conformation assumes an unphysical bond angle of 90°. After structure optimization (here, carried out using the PBE exchange-correlation functional), a much more physical angle results.

Practical Recommendation

One calculation - one directory.

On the computer intended to run FHI-aims, you will need to manage files using a command line interface and in different folders ("directories") intended to organize your data.

We strongly recommend to create a new directory (with its own control.in and geometry.in files) for every new FHI-aims calculation.

Rerunning and/or continuing FHI-aims calculations in the same directory that was used before will overwrite output and input files. In any simulation (including non-trivial failed attempts), keeping the data can be essential to later understanding and/or reconstructing what happened in a particular calculation.

## Set up the geometry.in file

Relaxation/Optimization

The terms structure optimization and relaxation are used interchangeably. The same holds for structure and geometry. Strictly speaking, "structure optimization" means a directed search for a local minimum of the Born-Oppenheimer potential energy surface.

A good initial guess of the geometry is a first important step for any system one would like to simulate. For most cases, starting from a random guess will not result in finding the global minimum-energy structure (if that is the intention). Thus, looking for experimental data or existing numerical studies are always recommended.

However, the water molecule is simple enough to start from an educated guess. A possible FHI-aims geometry input file geometry.in for an initial guess for the water molecule (not the final geometry) might be:

atom    0.00000000    0.00000000    0.00000000    O
atom    0.70700000   -0.70700000    0.00000000    H
atom   -0.70700000   -0.70700000    0.00000000    H

In order to carry out this first tutorial, we recommend to use an ASCII text editor, such as, vim or emacs, enabling us to work as closely as possible to the input files (which are ASCII files). Create a new folder for this first exercise and generate the geometry.in file for the water molecule using the above example initial guess for the input geometry.

The syntax of the geometry.in file is pretty much self-explanatory. A more detailed definition of the format can be found here or in the manual.

## Set up the control.in file

Our the next step is to define two sets of parameters, i.e. runtime choices and numerical settings to be used. Together, these parameters control the various aspects of a simulation with FHI-aims. These parameters must be provided in a file called control.in. Our objective is to carry out a structure relaxation of the water molecule using the PBE generalized gradient density functional approximation.

### Runtime choices

In the first step, we define the run-time choices for the planned simulation. Please create a new file control.in, located in the same folder as the geometry.in file, and include the following information:

xc                 pbe
relativistic       atomic_zora scalar
relax_geometry     bfgs 5e-3
spin               none

The first line defines the the desired exchange-correlation (XC) functional; here, PBE. The second line specifies the level of approximation used to include relativistic effects; in practice, a scalar-relativistic approach (as used here) is always preferable over a non-relativistic calculation. Read more about these first two lines in the following information box.

Level of Theory: Choice of Density Functional and Choice of Scalar Relativity

In order to fully define the potential energy surface explored by FHI-aims, the choice of a density functional is needed. This choice is not always easy and does require some knowledge about the strengths and weaknesses of a given density functional for different problems.

A second important choice in any DFT calculation is the approximation used to represent the kinetic energy of the electrons. In fact, the Dirac equation, not the Schrödinger equation, is the equation that captures electronic properties correctly. The relativistic Dirac equation captures essentially all (electronic) effects relevant to chemistry and materials science. The non-relativistic Schrödinger form of the kinetic energy is acceptable only for the lightest chemical elements and becomes progressively wrong as elements grow heavier - even for valence electrons, not just for deep core states.

Yet, the numerical form of the Schrödinger equation is computationally much simpler and also less demanding than Dirac's equation. In practice, the solution adopted by practically every production code is to replace the Schrödinger kinetic energy with a "scalar relativistic" kinetic energy that captures the most important aspects of relativity, while retaining the same numerical shape as the Schrödinger Equation.

In our experience, there is never a reason not to use scalar relativity - by sticking to a single form of the kinetic energy operator in all calculations, consistent total energy differences, eigenvalues and other observables can be ensured.

Conversely, mixing total energies from scalar- and non-relativistic treatments in energy differences could be a catastrophic idea since the total energies of the same system can be very different numbers for different relativistic treatment.

The particular form of scalar relativity employed in FHI-aims can and must be invoked in control.in through the keywords

relativistic atomic_zora scalar

The particular expressions that define this level of scalar relativity can be found in Eqs. (55) and (56) of Computer Physics Communications 180, 2175-2196 (2009) http://dx.doi.org/10.1016/j.cpc.2009.06.022.

More importantly, this level of scalar relativity has proven to be extremely robust and accurate for valence electrons related properties over time, consistently matching benchmark results of other high-precision codes in the community: see. e.g., http://dx.doi.org/10.1126/science.aad3000 for bulk cohesive properties and https://doi.org/10.1103/PhysRevMaterials.1.033803 for a broad range of energy band structures.

In short, we recommend to simply always use relativistic atomic_zora scalar as outlined above.

We will use a variant of the bfgs algorithm (BFGS: Broyden–Fletcher–Goldfarb–Shanno; more specifically, a trust radius based version) for geometry optimization, as described below. The keyword trm (trust radius method) can be used synonymously with bfgs. We will consider a local minimum of the total energy to be reached when the magnitude of all residual forces on the atomic nuclei is smaller than a predefined threshold, i.e., $$5\cdot10^{-3}$$ eV/Å. For a smaller threshold, the simulation can take rather long. On the other hand, a larger threshold may be insufficient to find the actual (local) minimum structure. Please find more details about the default relaxation algorithm in FHI-aims at the end of this first part.

Finding a Local Minimum-Energy Structure of a Molecule, Solid, or Nanostructure

Finding atomic coordinates and (optionally) lattice parameters that correspond to a local minimum of the potential-energy surface at a given level of theory is usually the first step of an electronic structure based simulation. This step ensures that the input structure used to derive observables later is mathematically well defined and that the structure does not contain any accidentally misconfigured chemical configurations or bonds that are far from thermodynamic equilibrium conditions.

It is important to remember that, for a given set of atoms, there is usually more than one such local minimum. Finding the global minimum and/or a physically meaningful structure cannot usually be accomplished by just entering a random initial structure and hoping for the best.

Setting up an initial structure verifying that it is physically plausible is a critical task that can be the key step in a scientific project. Often, this step amounts to formulating the overarching scientific question correctly in the first place.

Likewise, visualizing an input structure before using it for expensive follow-up simulations is a critical task that should never be skipped. Check for physically implausible bond lengths, bond angles, misplaced atoms, etc.

There are many good visualization applications. In this tutorial, we use GIMS. However, Jmol, Vesta, XCrysDen and others are all excellent and frequently used tools for visualization tasks related to FHI-aims.

### Species Defaults

Species Defaults: Choice of the numerical settings, including the basis set

FHI-aims uses numerically tabulated, atom-centered orbital (NAO) basis sets to represent the quantum mechanical orbitals and wave functions. These basis functions are specific to each chemical element ("species") used in the simulation.

FHI-aims also employs a variety of other numerical parameters that are specific to each atom and/or chemical element: Angular-momentum expansion order of the electrostatic (Hartree) potential, atom-centered real-space integration grids, etc.

Predefined files that contain these numerical definitions for each chemical element are located in species_defaults directory of the FHI-aims code distribution. At least four different levels of precision are provided for each chemical element: light, intermediate, tight, and really_tight. A more accurate choice of defaults means a higher computational cost of the simulation.

In order to obtain the complete control.in file, the "light" species defaults for the H and O atoms will still need to be appended to control.in.

At the command line, the appropriate command to accomplish this task is

cat [FHI-aims-directory]/species_defaults/defaults_2020/light/01_H_default >> control.in

and
cat [FHI-aims-directory]/species_defaults/defaults_2020/light/08_O_default >> control.in

where [FHI-aims-directory] should be replaced with the proper path to your FHI-aims distribution on the computer to be used.

A note on GIMS

The initial conformation geometry.in of the molecule and the initial control.in file can be built using GIMS and the "Simple Calculation" workflow app.

Using GIMS for this task does have the advantage that the "species defaults" for H and O will be added to control.in automatically. When using a text editor and the command line, the species defaults need to be appended to control.in manually using an appropriate command (see above).

A valid control.in for the water molecule

Your final control.in file should look as follows:

xc                                 pbe
relativistic                       atomic_zora  scalar
relax_geometry                     bfgs         5e-3

################################################################################
#
#  FHI-aims code project
#  VB, Fritz-Haber Institut, 2009
#
#  Suggested "light" defaults for O atom (to be pasted into control.in file)
#  Be sure to double-check any results obtained with these settings for post-processing,
#  e.g., with the "tight" defaults and larger basis sets.
#
################################################################################
species        O
#     global species definitions
nucleus             8
mass                15.9994
#
l_hartree           4
#
cut_pot             3.5  1.5  1.0
basis_dep_cutoff    1e-4
#
angular_grids specified
division   0.2659   50
division   0.4451  110
division   0.6052  194
division   0.7543  302
#      division   0.8014  434
#      division   0.8507  590
#      division   0.8762  770
#      division   0.9023  974
#      division   1.2339 1202
#      outer_grid 974
outer_grid 302
################################################################################
#
#  Definition of "minimal" basis
#
################################################################################
#     valence basis states
valence      2  s   2.
valence      2  p   4.
#     ion occupancy
ion_occ      2  s   1.
ion_occ      2  p   3.
################################################################################
#
#  Suggested additional basis functions. For production calculations,
#  uncomment them one after another (the most important basis functions are
#  listed first).
#
#  Constructed for dimers: 1.0 A, 1.208 A, 1.5 A, 2.0 A, 3.0 A
#
################################################################################
#  "First tier" - improvements: -699.05 meV to -159.38 meV
hydro 2 p 1.8
hydro 3 d 7.6
hydro 3 s 6.4
#  "Second tier" - improvements: -49.91 meV to -5.39 meV
#     hydro 4 f 11.6
#     hydro 3 p 6.2
#     hydro 3 d 5.6
#     hydro 5 g 17.6
#     hydro 1 s 0.75
#  "Third tier" - improvements: -2.83 meV to -0.50 meV
#     ionic 2 p auto
#     hydro 4 f 10.8
#     hydro 4 d 4.7
#     hydro 2 s 6.8
#  "Fourth tier" - improvements: -0.40 meV to -0.12 meV
#     hydro 3 p 5
#     hydro 3 s 3.3
#     hydro 5 g 15.6
#     hydro 4 f 17.6
#     hydro 4 d 14
# Further basis functions - -0.08 meV and below
#     hydro 3 s 2.1
#     hydro 4 d 11.6
#     hydro 3 p 16
#     hydro 2 s 17.2
################################################################################
#
#  FHI-aims code project
#  VB, Fritz-Haber Institut, 2009
#
#  Suggested "light" defaults for H atom (to be pasted into control.in file)
#  Be sure to double-check any results obtained with these settings for post-processing,
#  e.g., with the "tight" defaults and larger basis sets.
#
################################################################################
species        H
#     global species definitions
nucleus             1
mass                1.00794
#
l_hartree           4
#
cut_pot             3.5  1.5  1.0
basis_dep_cutoff    1e-4
#
angular_grids       specified
division   0.2421   50
division   0.3822  110
division   0.4799  194
division   0.5341  302
#      division   0.5626  434
#      division   0.5922  590
#      division   0.6542  770
#      division   0.6868 1202
#      outer_grid  770
outer_grid  302
################################################################################
#
#  Definition of "minimal" basis
#
################################################################################
#     valence basis states
valence      1  s   1.
#     ion occupancy
ion_occ      1  s   0.5
################################################################################
#
#  Suggested additional basis functions. For production calculations,
#  uncomment them one after another (the most important basis functions are
#  listed first).
#
#  Basis constructed for dimers: 0.5 A, 0.7 A, 1.0 A, 1.5 A, 2.5 A
#
################################################################################
#  "First tier" - improvements: -1014.90 meV to -62.69 meV
hydro 2 s 2.1
hydro 2 p 3.5
#  "Second tier" - improvements: -12.89 meV to -1.83 meV
#     hydro 1 s 0.85
#     hydro 2 p 3.7
#     hydro 2 s 1.2
#     hydro 3 d 7
#  "Third tier" - improvements: -0.25 meV to -0.12 meV
#     hydro 4 f 11.2
#     hydro 3 p 4.8
#     hydro 4 d 9
#     hydro 3 s 3.2


## Running the calculation

A parallel simulation with N processes can be launched by typing

mpirun -n N aims.x > aims.out 2>&1

or
mpirun -n N aims.x | tee aims.out


Here we write the output of the simulation to a file called aims.out file but you can choose any other name that you prefer.

The binary name aims.x should be replaced with whatever is the name of the FHI-aims binary file compiled by you (including the corresponding path, i.e., the location of the directory in which that file is located). For example, on this writer's laptop, the current binary name is /Users/blum/codes/fhi-aims/bin/aims.210716_1.scalapack.mpi.x.

The mpirun command facilitates the parallel execution of the code and can have different names that depend on the particular computer system and/or queueing system used. The actual name and usage are usually documented by the computing center in question.

## Checking the results

After the FHI-aims calculation is complete (see above for the command to execute FHI-aims), the fully relaxed structure ("light") is written to geometry.in.next_step and looks like this:

#
# This is the geometry file that corresponds to the current relaxation step.
# If you do not want this file to be written, set the "write_restart_geometry" flag to .false.
#  aims_uuid : 9496E7B2-D8E2-4C5F-AA87-D4B66DF85A6B
#
atom       0.00000000     -0.07320102      0.00000000 O
atom       0.76731112     -0.67039949      0.00000000 H
atom      -0.76731112     -0.67039949      0.00000000 H
#
# What follows is the current estimated Hessian matrix constructed by the BFGS algorithm.
# This is NOT the true Hessian matrix of the system.
# If you do not want this information here, switch it off using the "hessian_to_restart_geometry" keyword.
#
hessian_file


A file hessian.aims is also written. As outlined below, the geometry.in.next_step and hessian.aims files could be used to continue the relaxation with "intermediate" or "tight" species defaults and obtain fully converged results.

It is furthermore worthwhile to take a look at the actual FHI-aims output stream (written into a file aims.out in our example).

This stream begins with some important information regarding how the calculation was executed, then includes complete copies of control.in and geometry.in, details all self-consistent field iterations, as well as relaxation steps. The aims.out ends with the line:

Have a nice day.


which is used as FHI-aims' indicator that the calculation converged as expected. (Also, this line was among the first, if not the first line written in FHI-aims in 2004 and the sentiment is absolutely genuine.)

A bit before the end of the standard output aims.out, one can identify the final total energy delivered by the s.c.f. cycle, as well as the final geometry reached by the geometry optimization.

This information, and more, can also be obtained by opening the "Output Analyzer" App of GIMS and importing the respective aims.out file there. One advantage of GIMS is that it does provide a graphical analysis of the calculation, particularly of the self-consistent field cycle and of successive geometry steps until a converged geometry is found.

## Restarting a Structure Relaxation

The files geometry.in.next_step and hessian.aims are the two key pieces needed to restart a structure relaxation from the results of an earlier relaxation. Such a restart can be desirable for several reasons, e.g, in order to:

• Continue an unfinished relaxation that has stopped (for example, due to a queue wall time limit or because the earlier FHI-aims relaxation run encountered a problem and stopped with an error message).
• Use the data created by a pre-relaxation run using computationally cheaper "light" settings and follow up with a post-relaxation using more accurate "intermediate" or "tight" settings.
• Use the data generated by pre-relaxing with a cheaper density functional (e.g., PBE) as a starting point to initialize a post-relaxation with a more expensive density functional (e.g., HSE06). Read more details here.

The method to restart a structure relaxation from an earlier one is simple:

1. Run the pre-relaxation, which generates intermediate files geometry.in.next_step and hessian.aims
2. After the pre-relaxation is finished, create a new directory for the followup relaxation
3. Copy the geometry.in.next_step and hessian.aims files from the pre-relaxation directory to the followup-relaxation directory.
4. Change the name of the copy of geometry.in.next_step to geometry.in.
5. Create a new control.in file in the followup-relaxation directory
6. Rerun FHI-aims in the follow-up directory.

We will demonstrate the restart function for a followup relaxation using the hybrid-density functional HSE06 with intermediate species defaults.

Create a new folder and copy geometry.in.next_step and hessian.aims from the above PBE-based relaxation into it. Rename geometry.in.next_step to geometry.in. The top of your control.in file for the followup relaxation should look like this:

xc                                 hse06 0.11
hse_unit                           bohr
relativistic                       atomic_zora  scalar
relax_geometry                     bfgs         5e-3
spin                               none

The first two lines indicate that the HSE06 functional will be used for the calculation.

Please attach the intermediate species defaults for O and H to the end of this file. Again, this can be done with the following commands:

cat [FHI-aims-directory]/species_defaults/defaults_2020/intermediate/01_H_default >> control.in

and
cat [FHI-aims-directory]/species_defaults/defaults_2020/intermediate/08_O_default >> control.in


Now, you are ready to start your followup relaxation (again, use the mpirun ... given above).

After the calculation has finished, take a look into the output file. First, FHI-aims found the restart information (after parsing the geometry.in file):

Found keyword 'hessian_file', indicating that this is the continuation of a structure optimization.
Verifying the existence of the file 'hessian.aims' on all processors

Next, look at the total calculation time at the end of the output file. You may note that the relaxation took longer. This is due to a more expensive hybrid XC functional and the use of more precise species defaults (intermediate), which are also computationally more demanding.

To understand the difference between restarting from a pre-relaxed calculation and starting a relaxation from scratch, we will now repeat the same relaxation without restart files.

For this purpose, create a new folder and copy geometry.in and control.in from the first HSE06 relaxation into the new folder. Open geometry.in and remove the last two lines. The file should look like this:

#
# This is the geometry file that corresponds to the current relaxation step.
# If you do not want this file to be written, set the "write_restart_geometry" flag to .false.
#  aims_uuid : 9496E7B2-D8E2-4C5F-AA87-D4B66DF85A6B
#
atom       0.00000000     -0.07320102      0.00000000 O
atom       0.76731112     -0.67039949      0.00000000 H
atom      -0.76731112     -0.67039949      0.00000000 H
#
# What follows is the current estimated Hessian matrix constructed by the BFGS algorithm.
# This is NOT the true Hessian matrix of the system.
# If you do not want this information here, switch it off using the "hessian_to_restart_geometry" keyword.
#


Now, start this calculation.

We can now compare the two relaxations from their output files.

• HSE06, intermediate, with restart files:
Computational steps:
| Number of self-consistency cycles          :           34
| Number of SCF (re)initializations          :            3
| Number of relaxation steps                 :            2

• HSE06, intermediate, without restart files:
Computational steps:
| Number of self-consistency cycles          :           49
| Number of SCF (re)initializations          :            5
| Number of relaxation steps                 :            4


One can immediately see the effect of using restart files. The relaxation using the restart files from a pre-relaxation needs fewer relaxation steps to converge.

We always recommend to begin a relaxation run with "light" settings and, in some cases, even with a cheaper density functional (assuming a reliable cheaper density functional can be found). The resulting pre-optimized geometry (geometry.in.next_step) and Hessian matrix (hessian.aims) can be MUCH better than their corresponding initial guesses.

In contrast, beginning a relaxation from a bad geometry guess while using highly accurate numerical settings throughout the entire structure optimization can be very expensive and entirely unnecessary (most of the intermediate structures found during the relaxation process will likely never be used again). Thus, always:

• Pre-relax with appropriate, relatively cheap settings (usually, "light")
• Post-relax using tighter settings in a second run, using the information generated during the pre-relaxation as a starting point.

We are done with the first part of the tutorial! Below, we provide some additional information, which can greatly help with later workflows/research with FHI-aims.

• A variant of the Broyden–Fletcher–Goldfarb–Shanno (bfgs) algorithm, enhanced by a trust-region method (trm), is employed in FHI-aims to find a local minimum of the potential energy surface. Such a minimum is considered to be found if the moduli of all forces should be smaller than a small, configurable threshold value, typically smaller than $$5\cdot10^{-3}$$ eV/Å. Within FHI-aims, the terms bfgs and trm are used to denote the same algorithm.
• It is highly worthwhile to understand at least the basics of the optimization algorithms mentioned above - follow the links. Importantly, this shows that the bfgs algorithm relies on a good approximation of the Hessian matrix, that is, the second derivatives of the total energy with respect to all combinations of atomic and/or lattice coordinates for a given structure. A good initial guess for the Hessian matrix at the outset of a structure optimization can be immensely helpful to speed up the calculation. We use a variant of the initial guess suggested by Lindh and coworkers but if this does not work, a simple diagonal initial guess for the Hessian can also be used.
• Even more importantly, the bfgs algorithm updates and improves its own guess of the Hessian matrix as the structure optimization progresses. Structure optimization is a step-wise process. At each optimization step, FHI-aims writes the current atomic/lattice geometry into a file geometry.in.next_step. It also writes a file hessian.aims, which contains the current guess of the Hessian matrix determined by the bfgs algorithm.
• To start a new simulation with the relaxed structure, geometry.in.next_step should be used as the geometry.in file of the new simulation, and hessian.aims should be present in the new simulation folder. For example, one can copy old-simulation/geometry.in.next_step into new-simulation/geometry.in, prepare appropriate control.in, and then start the new simulation.