# 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
```

*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
```

## 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
```

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

`[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
#
radial_base 36 5.0
radial_multiplier 1
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
#
radial_base 24 5.0
radial_multiplier 1
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
```

```
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.
#
trust_radius 0.2000000030
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:

- Run the pre-relaxation, which generates intermediate files
`geometry.in.next_step`

and`hessian.aims`

- After the pre-relaxation is finished, create a new directory for the followup relaxation
- Copy the
`geometry.in.next_step`

and`hessian.aims`

files from the pre-relaxation directory to the followup-relaxation directory. - Change the name of the copy of
`geometry.in.next_step`

to`geometry.in`

. - Create a new
`control.in`

file in the followup-relaxation directory - 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
```

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
```

```
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
```

*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.

## Additional Information

### Notes About Structure Relaxation in FHI-aims

The terms "local structure optimization" and "structure relaxation" are used interchangeably in the following. They mean the same thing.

- FHI-aims can calculate gradients of the total energy with respect to atomic positions ("forces") and lattice parameters ("stresses") for local-density, generalized-gradient and meta-generalized gradient density functionals, as well as hybrid density functionals.
- 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.