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
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 : C138C83D-E5EA-45F0-A93E-9EA40FD4A8F9
#
atom 0.00000000 -0.07320103 -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_stepandhessian.aims - After the pre-relaxation is finished, create a new directory for the followup relaxation
- Copy the
geometry.in.next_stepandhessian.aimsfiles from the pre-relaxation directory to the followup-relaxation directory. - Change the name of the copy of
geometry.in.next_steptogeometry.in. - Create a new
control.infile 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
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 : C138C83D-E5EA-45F0-A93E-9EA40FD4A8F9
#
atom 0.00000000 -0.07320103 -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 : 30 | Number of SCF (re)initializations : 3 | Number of relaxation steps : 2 - HSE06, intermediate, without restart files:
Computational steps: | Number of self-consistency cycles : 46 | 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.
Solutions
You find all the solution to all the above exercises by clicking on the button below.
Additional Information
Checking the input files for trivial errors
This is about the dry_run keyword.
What good is a "dry run" (i.e., a test of the input files before running the actual calculations)?
Upon starting the FHI-aims calculation, the code will first check the input
keywords in both control.in and geometry.in for consistency and (hopefully)
provide an error message in case of any typographic errors or known inconsistent
keyword combinations that would prevent a correct run of FHI-aims.
We could just run FHI-aims at the command line (see above) and see what happens.
In many real world scenarios, however, FHI-aims will be run through a queueing system, i.e., everything will be set up, a job will be queued, the user goes away and waits for the supercomputer to execute the code. In this case, stopping with a trivial error can cost serious amounts of real-world time since the error would need to be corrected and the task re-queued.
FHI-aims does allow one to check the consistency of input files at the command line without executing the full run if the keyword
dry_run
is included in the control.in file. In that case, the code will just go through all motions
to set up the calculation but not actually execute the time consuming parts. If you're so inclined,
modify your control.in file to include the dry_run keyword and see what happens if you execute
FHI-aims at the command line. This should tell you if the input files would be accepted for a
later, real-world run, helping you to identify incorrect keywords or typos ahead of time.
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 termsbfgsandtrmare 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
bfgsalgorithm 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
bfgsalgorithm 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 filegeometry.in.next_step. It also writes a filehessian.aims, which contains the current guess of the Hessian matrix determined by thebfgsalgorithm. - To start a new simulation with the relaxed structure,
geometry.in.next_stepshould be used as thegeometry.infile of the new simulation, andhessian.aimsshould be present in the new simulation folder. For example, one can copyold-simulation/geometry.in.next_stepintonew-simulation/geometry.in, prepare appropriatecontrol.in, and then start the new simulation.