3.52 Output options

Tag: verbatim_writeout (geometry.in)
        Usage: verbatim_writeout flag
Purpose: Enables or suppresses the writing of geometry.in to the FHI-aims standard output stream exactly as it is read the first time.
flag is a logical variable (.true. or .false.). Default: .true. .

By default, geometry.in is now written (copied) verbatim into the FHI-aims standard output as it is parsed for the first time, allowing to reproduce exactly any FHI-aims calculation simply by copy-pasting that part to a new geometry.in file.

If verbatim_writeout is set to false anywhere in geometry.in, no writing will occur from that point on forward.

The exact same option (same keyword / syntax) can also be used in control.in, producing the same effect there.

Note that the keyword has the same name in geometry.in and control.in, and is therefore only documented as a clickable link for control.in. Apologies for this omission.

When the output cube functionality is requested in order to print a three-dimensional array on an even-spaced grid, FHI-aims can use internal defaults to determine the dimensions of this cube file. Each cube file is based on an even-spaced grid with $n$ grid points, where $n=n_1\times n_2\times n_3$ and $n_i$ is the number of grid segments used along each of the grid directions $i$=1, 2, 3. The output cube edge subtag can be used to specify the grid spacing. However, if that subtag is not set, a default value of 0.1 Å is used for discretization of the three cube edges. For large structures, this could lead to very large cube file sizes. As a safeguard, FHI-aims stops when the number of points in a cube, $n$, exceeds the number set by cube_default_size_safeguard. Its default value, 5$\cdot$10⁷, corresponds to a cube file size of about 800 MB (i.e., 16 bytes of information are stored for each point of a cube grid). If the edges of any cube files are set explicitly using the output cube edge subtag, the criterion set by cube_default_size_safeguard is not checked and the code does noty stop.

Only useful in conjunction with the keywords output dos, output dos_tetrahedron or output postscf_eigenvalues, and only for periodic systems.

In a periodic calculation, one usually specifies the basic k_grid used to obtain the self-consistent electron density, total energy etc. Such $k$-space grids are usually fairly sparse, and if a density of states (DOS) is calculated directly from the eigenvalues stored at these $k$-points only, the DOS will either look choppy, or (after significant broadening), smooth, broad, and blurred.

A simple remedy is to use the original k_grid while approaching self-consistency as usual, but then compute the DOS using an auxiliary $k$-grid that is made denser by factors $n_1$, $n_2$, $n_3$, respectively. For example, the settings
k_grid 10 10 10
output dos […]
dos_kgrid_factors 8 8 8
mean that the s.c.f.-cycle itself is run with a 10$\times$10$\times$10 $k$-point grid, but subsequently, a density of states is computed with an 80$\times$80$\times$80 $k$-point grid.

Available options for verbosity are:

• •

  none : No output from ELSI. This is the default if the overall output level of FHI-aims is MD_light.

• •

  light : Enables output from the ELSI interface, but no output from the solvers.

• •

  detail : Enables output from the ELSI interface as well as the solvers. When using libOMM or PEXSI, additional output will be written to an separate log file.

• •

  debug : Enables the same output as does the detail option, with additional memory usage information. Creates large output files and thus should not be chosen in production runs.

• •

  json : Enables the output of the runtime parameters used in ELSI in a separate JSON file, powered by the FortJSON library. May be used on top of the above options.

Available options for matrix are:

• •

  hamiltonian : Outputs the SCF converged Hamiltonian matrix.

• •

  overlap : Outputs the overlap matrix.

• •

  density_matrix : Outputs the SCF converged density matrix.

• •

  eigenvectors : Outputs the SCF converged KS eigenvectors.

Note that this keyword outputs matrices in the $k$-space instead of the real space. Therefore, the number of matrix files would be equal to the number of $k$-points, multiplied by the number of spin channels. The output files are in the ELSI CSC binary format (see the documentation of ELSI). The Python scripts in the “utilities/elsi_matrix” directory may be used to post-process an ELSI matrix file, e.g., converting it to a human-readable format.

Available options for format are:

• •

  plain : Outputs the matrices in human-readable format.

• •

  HDF5 : Outputs the matrices using HDF5 format.

Note that this keyword outputs matrices in the real-space instead of the $k$-space. Therefore, the number of Hamiltonian matrix files is number of spin channels. Python scripts in the “utilities/ will soon be added to parse the output files. If these scripts are desirable please contact the FHI-aims developers.

This keyword will only work for periodic systems and it does not cover meta-GGA or hybrid density functionals.

This option requires that a reference $z$ coordinate for the electrostatic potential evaluation deep in the vacuum be provided by hand, through the keyword set_vacuum_level. The surface must be parallel to the $xy$ plane.

The output for the “upper’ and “lower” surface of the slab (larger and smaller $z$ value, respectively), will be printed at the end of the first SCF step and at the end of the last SCF step where the chemical potential is included too. For ease of evaluation, the self-consistent work function will be repeated again at the end of the output file in the summary block. In the case of a semiconducting slab (band gap > 0.1 eV), the output also includes the work function referenced to the VBM and CBM.

This behavior is shown in Fig 3.16, the work function is defined as the difference between the Fermi surface and vacuum level which are different on the saturated and unsaturated surfaces. Accordingly, the work function exhibits the same asymmetry. In practice, this behavior is correctly reproduced only if the use_dipole_correction is additionally specified.

If the include_spin_orbit is requested, the work function output, both neglecting and including the post-SCF SOC corrections, will be printed directly after the output of the SOC calculation.

If the work function output is requested, keyword compensate_multipole_errors is now automatically switched on by default. This will change total energies slightly compared to the uncompensating case, and – we believe – even for the better. It will certainly lead to a better description of the long-range Hartree potential.

However, it must be possible to find a vacuum plane $z$, where the surface dipole is compensated, that is further than 6 Å away from the nearest atom. Otherwise, the calculation will stop and alert the user.

Specifically, the reference Hartree potential component for the work function evaluation is only the long-range (reciprocal-space) Hartree potential term of the Ewald sum, not the full electrostatic potential. Thus, the vacuum level must be specified in a region where all real-space components of the electrostatic potential have safely died away to zero. One may achieve this by increasing the vacuum layer thickness, which can be done at very small overhead cost in FHI-aims.

The list of additional output types is given as a separate subsubsection below.

The maximization procedure follows JCP 135, 134107 (2011). Currently only the cartesian position of the centers is outputted in xyz format in geometry_boys.xyz. The transformation matrix is also calculated, but only applied to the orbitals when using the keyword apply_boys.

By default, cube files are written once, after the SCF cycle has converged. With this keyword,all cube files specified in control.in will be written each $n^{th}$ iteration, where n is an integer greater than zero. This keyword should be helpful to analyse what is going on during subsequent SCF cycles. However, the output of cubes is usually quite slow, so choosing this option will slow down the calculation a lot.

This keyword does not support spin-orbit coupling, as spin-orbit coupling is applied after the SCF cycle has converged.

This keyword is only applicable when the output cube keyword(s) are being used; please see the manual entry for output cube for more information.

In some, atoms in FHI-aims can unexpectedly “switch” unit cells during relaxation. This has no effect on the calculation, but the output geometry coordinates (written to the standard output stream) will look strangely detached when visualized. By default, FHI-aims maps its coordinates back to the first unit cell anyway, but this behavior can be forcibly switched off if so requested (makes nicer movies).

If increased to full, the Kohn-Sham eigenvalues of every s.c.f. iteration are written to the standard output file. For single-point calculations, this may be quite desirable, but leads to unmanageable file sizes for long relaxation or molecular dymnamics runs.
Another option, useful for long molecular dynamics (MD) runs, is MD_light. It writes standard output only in the initialization part and at the end of each MD step, while a minimal output is written for the single scf cycle.

If set to .false., FHI-aims will check during the output of cube files whether a file with the selected file name already exists. If such a file is found, the file name will be changed by adding a number (1,2,3…) to the end of the file name. This is very useful when relying on default names or when plotting cubes during the SCF cycle.

If set to .true. FHI-aims will not check during output whether a file with the selected name already exists. If it does exist, it will be simply overwritten!

This keyword is only applicable when the output cube keyword(s) are being used; please see the manual entry for output cube for more information.

By default, control.in is now written (copied) verbatim into the FHI-aims standard output as it is parsed for the first time, allowing to reproduce exactly any FHI-aims calculation simply by copy-pasting that part to a new control.in file.

If verbatim_writeout is set to false anywhere in control.in, no writing will occur from that point on forward.

The exact same option (same keyword / syntax) can also be used in geometry.in, producing the same effect there.

TREXIO [253] is a powerful file format for storing wave functions in quantum chemistry. Its goal is to allow the exchange of wave functions between different codes. For example, a Hartree-Fock calculation can be performed in FHI-aims using a NAO basis set. Through TREXIO, this wave function is then read in by other compatible programs so that post-Hartree-Fock methods can be performed using a NAO basis set usually not available to such programs. The information exported by this module does not only allow to do post-Hartree-Fock methods, but also real space methods like Diffusion Monte Carlo (DMC), among many other options.

Whenever a TREXIO export is performed, the "nucleus" and "electron" groups are written regardless of further keywords. For the definition of these groups, please refer to the TREXIO documentation and GitHub repository.

This sets the information in the "basis" and "ao" groups of the TREXIO file. The basis functions are exported as the splines used within FHI-aims. A TREXIO file generated with this argument can e.g. be used for further analysis of the basis functions or for real-space methods like diffusion Monte Carlo.

This sets the information in the "mo" group of the TREXIO file. This includes the coefficient matrix to convert between the AO and MO bases, occupation numbers, eigenvalues and spin information in the UHF case. Depending on the kind parameter, the overlap matrix, core Hamiltonian and electron-repulsion integrals are exported for the desired bases. The generated TREXIO file can be used for post-Hartree-Fock methods, such as variants of configuration interaction (CI), coupled cluster (CC) or other methods such as Full-Configuration-Interaction Quantum Monte Carlo (FCIQMC). Having the integrals in both bases is necessary for orbital rotations in some programs.

Please, look at Chapter 5 for a comprehensive description on how to perform transport simulations across nanoscale objects.

This option is based on a Mulliken Analysis and shares its syntax with output dos and output species_proj_dos. See also section 4.4 for more details.

There are two types of output files for each atom:

• •

  atom_proj_dos_number_raw.dat, where number denotes the atom number in the order of geometry.in. This file contains the total and angular-momentum resolved DOS components as a function of eigenvalue energy (first column) as used internally in FHI-aims. The energy zero is then given by the vacuum level (non-periodic systems) or by the $𝐆$=0 component of the long-range Hartree potential (periodic systems).

• •

  atom_projected_dos_number.dat, which gives the same information, except that the energy zero is shifted to the Fermi energy (metallic systems) or valence band maximum (insulators), respectively.

Note that projected densities of states such as given here must be based on some kind of projection orbitals, the choice of which is somewhat arbitrary by necessity. This is thus a tool for qualitative analyses.

In FHI-aims, we project directly on the atom-centered angular-momentum components as defined by the overlapping basis set. This definition becomes the more arbitrary the larger the basis set, just like a mulliken analysis. The closer the full basis comes to completeness, the more ambiguous will a mulliken-like analysis become, since it may not be a priori clear which electrons should be counted towards the basis functions of one atom vs. those of another atom. Thus, do not expect a pDOS to simply converge as the basis set size is increased; use it as a qualitative indicator of trends, but nothing more.

This keyword supports spin-orbit coupling. When spin-orbit coupling is enabled, the file(s) containing the spin-orbit-coupled DOS will have the default filename(s) and the file(s) containing the scalar-relativistic (i.e. no SOC) DOS will have an additional suffix “.no_soc”. Note that if you requested spin collinear in the control.in file, there will be only one file per atom (and equivalently for all _raw.dat files) containing the projected DOS of all spin-coupled states:

• •

  atom_proj_dos_number.dat

The reason is that the separate spin channels of scalar relativity no longer exist after the SOC operator is applied – the states now form a single set. However, there are two DOS files for the output without spin-orbit coupling; one for each spin channel:

• •

  atom_proj_dos_spin_upnumber.dat.no_soc

• •

  atom_proj_dos_spin_dnnumber.dat.no_soc.

As of version 240605 this keyword can now be used in combination with
dos_kgrid_factors. The n_broadening parameter was also introduced in this version.

When this keyword is combined with atom_proj_dos_tetrahedron,
species_proj_dos_tetrahedron, or mulliken, and dos_kgrid_factors is not used, then the n_broadening parameter will be ignored. This is because a Mulliken decomposition over all states is performed for each of these other keywords, which is then re-used for the projected DOS with atom_proj_dos.

When this keyword is combined with species_proj_dos, and different values of Estart, broadening, n_broadening, or Eend are used for each, the code will take the most extreme limit from the values of Estart-broadening*n_broadening and Eend + broadening*n_broadening for the atom and species projected DOS, to ensure that a single Mulliken decomposition is done that utilizes states to encompass both projections, such that two separate decompositions do not need to be carried out.

This is the keyword that should be used to obtain atom-resolved densities of states with high integration resolution.

This option is based on a Mulliken Analysis and shares its syntax with output dos_tetrahedron and output species_proj_dos_tetrahedron. See also section 4.4 for more details.

There are two types of output files for each atom:

• •

  atom_proj_dos_number_tetrahedron_raw.dat, where number denotes the atom number in the order of geometry.in. This file contains the total and angular-momentum resolved DOS components as a function of eigenvalue energy (first column) as used internally in FHI-aims. The energy zero is then given by the vacuum level (non-periodic systems) or by the $𝐆$=0 component of the long-range Hartree potential (periodic systems).

• •

  atom_projected_dos_number_tetrahedron.dat, which gives the same information, except that the energy zero is shifted to the Fermi energy (metallic systems) or valence band maximum (insulators), respectively.

Note that projected densities of states such as given here must be based on some kind of projection orbitals, the choice of which is somewhat arbitrary by necessity. This is thus a tool for qualitative analyses.

In FHI-aims, we project directly on the atom-centered angular-momentum components as defined by the overlapping basis set. This definition becomes the more arbitrary the larger the basis set, just like a mulliken analysis. The closer the full basis comes to completeness, the more ambiguous will a mulliken-like analysis become, since it may not be a priori clear which electrons should be counted towards the basis functions of one atom vs. those of another atom. Thus, do not expect a pDOS to simply converge as the basis set size is increased; use it as a qualitative indicator of trends, but nothing more.

This keyword supports spin-orbit coupling. When spin-orbit coupling is enabled, the file(s) containing the spin-orbit-coupled DOS will have the default filename(s) and the file(s) containing the scalar-relativistic (i.e. no SOC) DOS will have an additional suffix “.no_soc”. Note that if you requested spin collinear in the control.in file, there will be only one file per atom (and equivalently for all tetrahedron_raw.dat files) containing the projected DOS of all spin-coupled states:

• •

  atom_proj_dos_tetrahedron_number.dat

The reason is that the separate spin channels of scalar relativity no longer exist after the SOC operator is applied – the states now form a single set. However, there are two DOS files for the output without spin-orbit coupling; one for each spin channel:

• •

  atom_proj_dos_tetrahedron_spin_upnumber.dat.no_soc

• •

  atom_proj_dos_tetrahedron_spin_dnnumber.dat.no_soc.

Several bands can be plotted; FHI-aims outputs one file per specified output band line.

The band structure output files are named bandXYYY.out, where the letters X and YYY are replaced with numbers in the actual output file names:

• •

  The letter X encode the spin channel. In a non-spinpolarized or in a spin-orbit coupled calculation, X will always be 1. In a spin-polarized calculation, X=1 indicates the first spin channel, X=2 indicates the second spin channel.

• •

  The letters YYY are the consecutive numbers of the bands (starting from 001) requested in the control.in file, that is, the bands given in the order of output lines in control.in.

The files bandXYYY.out have the format
ipoint k1 k2 k3 occ1 E1 occ2 E2 ... occN EN,
i.e. they specify not only the bands for each $k$-point but also the occupation number for this particular point.

A safe starting value for n_points when performing semi-local calculations is 21. We have found that this generally samples the fine features of the k-path reasonably well, even for small unit cells with correspondingly large Brillouin zones. For hybrid-functional calculations, due to the computational expense one should consider using a smaller value. We also note that n_points includes the end-points, i.e. n_points=21 will give 20 intervals for a given branch. For comparison with results calculated by other DFT codes, it’s recommended to use values of form 1, 6, 11, 16, 21, … to ensure that the reciprocal coordinates are nice, simple fractions.

Note that a fully occupied band has the occupation number 2.0 in a non-spinpolarized calculation. In a spin-polarized or spin-orbit-coupled calculation, the maximum occupation number is 1.0.

Since this format contains all the important information, but is not particularly useful for actually plotting the band structure, we provide a small script which is described in section 4.4.

Note: the last two input options are technically not needed by the FHI-aims main program, but they are seriously helpful when turning this data into a plot and are used by the band plotting script provided along with this distribution, see Section 4.4 for details.

For periodic calculations, the eigenvectors, overlap matrices, and hamiltonian matrices at each k-point requested by output band can be written out using the output eigenvectors, output overlap_matrix, and output hamiltonian_matrix keywords, respectively.

This keyword supports spin-orbit coupling. When spin-orbit coupling is enabled, the file(s) containing the spin-orbit-coupled band structures will have the default filename(s) and the file(s) containing the scalar-relativistic (i.e. no SOC) band structures will have an additional suffix “.no_soc”. Note that if you requested spin collinear in the control.in file in addition, there will be only one file per band segment containing all spin-coupled states (output with spin-orbit coupling):

• •

  band1YYY.dat

The reason is that the separate spin channels of scalar relativity no longer exist after the SOC operator is applied – the states now form a single set. However, there are two band segment files for the output without spin-orbit coupling; one for each spin channel:

• •

  band1YYY.dat.no_soc

• •

  band2YYY.dat.no_soc

This keyword allows to get the band structure along a certain reciprocal-space direction, but only at those $k$-points that are already used during the s.c.f. calculation. If there are no appropriate $k$-points, no band structure is printed. If there are appropriate $k$-points, they are printed in the same format as the normal band structure from outputband, although some additional editing may be required to get a clean plot.

The purpose of this keyword is to allow to extract a band structure even in cases when the normal outputband functionality is experimental or, for some reason, not available.

This keyword allows to calculate the mulliken charge analysis at all K points along the band K path. The file name is named as bandmlk1001.out, bandmlk1002.out, …, etc. In the output file, the mulliken charge data is written first by K point, then by eigenstates. In each eigenstate, the mulliken charge on all atoms is written line by line. In each line, the following information is written: eigenstate ID, eigenvalue, occupation number, atom ID, spin ID, total mulliken charge, mulliken charge for l = 0, 1, 2,.. etc. For each state at a given point, the total mulliken charge for all atoms should sum up to 1 or 2, for non-SOC and SOC, respectively.

There is a keyword called band_mulliken_orbit_num specifying how many orbitals (states) to be written out. If the number following the keyword band_mulliken_orbit_num is I, then the orbitals in the range of HOMO - I + 1 and LUMO + I would be written out in the output file. The default value of band_mulliken_orbit_num is 50 for non-SOC and 100 for SOC.

To plot the band structures with Mulliken decomposition, two python script in the utilities directory in FHI-aims distribution can be used, i.e., band_mlk.py and band_mlk_soc.py for non-SOC and SOC, respectively. Running of these scripts is instructed in the first lines of these python files.

This output option allows to visualize the basis functions used, as well as some of the other defining pieces of the basis set. Note that the output is written for each grid point of the dense 1-dimensional logarithmic grid, with the radius given in bohr.

Specifically, this option produces the following types of files:

• •

  A$i$_$j$_$nl$_base.dat : Atomic (minimal basis) radial function $u(r)$ after the basis-confining potential was applied, for species number $i$, atomic-like (minimal) radial function number $j$, radial and angular quantum numbers $nl$.

• •

  A$i$_$j$_$nl$_base_kin.dat : Corresponding kinetic energy expression $[ϵ-v(r)]\cdot u(r)$

• •

  C$i$_$j$_$nl$_base.dat : confined free-atom like radial function $u(r)$ number $j$ for species number $i$, radial and angular quantum numbers $nl$.

• •

  C$i$_$j$_$nl$_base_pot.dat : Corresponding basis-defining potential including confining potential

• •

  El_base_$n$_$l$.dat : Radial function $u(r)$ after the basis-confining potential was applied for the species named El, radial and angular quantum numbers $nl$ (same as A$i$_$j$_$nl$_base.dat).

• •

  El_base_pot.dat : Basis-defining potential for atomic (minimal) radial functions of the species named El, after addition of the confining potential (as defined by cut_pot).

• •

  El_base_pot.dat Free-atom density of species named El (same as El_base_pot.dat).

• •

  El_free_$n$_$l$.dat: Radial function $u(r)$ before the basis-confining potential was applied for the species named El, radial and angular quantum numbers $nl$

• •

  El_free_pot.dat Spherical self-consistent free-atom potential of the species named El (implicitly confined by the cut_free_atom potential, but this artificial part is here not included)

• •

  El_free_rho.dat Free-atom density of the species named El.

• •

  H$i$_$j$_$nl$_base.dat : hydro radial function $u(r)$ number $j$ for species number $i$, radial and angular quantum numbers $nl$.

• •

  H$i$_$j$_$nl$_base_kin.dat : Corresponding kinetic energy expression $[ϵ-v(r)]\cdot u(r)$

• •

  I$i$_$j$_$nl$_base.dat : ionic radial function $u(r)$ number $j$ for species number $i$, radial and angular quantum numbers $nl$.

• •

  I$i$_$j$_$nl$_base_pot.dat : Corresponding basis-defining potential including confining potential

• •

  ty_$i$_$j$_$n$_$l$.dat : After on-site orthonormalization, radial function $u(r)$ of type ty (atomic, confined, hydro, ionic, …), for species number $i$, radial function number $j$, radial and angular quantum numbers $n,l$.

• •

  kin_ty_$i$_$j$_$n$_$l$.dat : Corresponding kinetic energy expression after on-site orthonormalization.

• •

  S$i$_$j$_$nl$_base.dat : Slater-type orbital radial function $u(r)$ for species number $i$, radial function number $j$, radial and angular quantum numbers $nl$.

• •

  S$i$_$j$_$nl$_base_kin.dat : Corresponding kinetic energy expression $[ϵ-v(r)]\cdot u(r)$

This keyword outputs information about the batch distribution used by FHI-aims to evaluate real-space quantities (real-space Hamiltonian, charge density update, etc.) Every MPI task creates a batch_statistics_task_###.dat file, where ### is the MPI task’s rank, and statistics for each batch are output sequentially to file. Statistics output for each batch include:

• •

  Number of points in the batch

• •

  Minimum and maximum number of radial basis functions evaluated on batch

• •

  Maximum number of basis functions evaluated on batch

• •

  Minmum and maximum number of atoms whose basis functions are evaluated on the batch

• •

  Minimum and maximum values for the integration weights for points in batch

Note: As of this writing, the output files will be rewritten with every SCF restart, including SCF reinitialization, geometry relaxation steps, and MD steps.

This keyword can be specified in combination with a SCF-calculation and will calculate and/or write to file quantities CC4S requires after the SCF has been completed.

If the xc-functional used is Hartree-Fock the output files generated by CC-aims can be used by CC4S to perform a canonical Coupled cluster calculation. In principal, one can also specify different xc-functionals and the quantities for a non-canonical Coupled cluster calculation will be calculated correctly. However, at this point in time, the non-canonical Coupled cluster algorithms are not yet implemented.

Note, that the CC-aims interface relies on ScaLAPACK, so that KS_method parallel must be used.

A number of additional keywords have been added to FHI-aims to give the user more control over the calculation of the CC4S-quantities. These can be found in Section CC4S.

The “cube” file format originates from the Gaussian code, but publically available descriptions exist, for example here:
http://paulbourke.net/dataformats/cube/
It is accepted by multiple visualization tools.

Visualization tools which may be used to plot cube file and are available for all major operating systems are Avogadro (https://avogadro.cc), jmol (http://www.jmol.org), and VMD (http://www.ks.uiuc.edu/Research/vmd/). For periodic systems, excellent success has been reported by multiple users using Vesta (https://jp-minerals.org/vesta/en/). Please see the documentation of those programs for more information on plotting the resulting cube files.

By default, the cube files will be output once, after the SCF cycle has converged, and FHI-aims will avoid overwriting pre-existing cube files it finds by appending numbers to the end of new file names. To output the cube files at regular intervals during the SCF cycle, use the output_cube_nth_iteration keyword, and to overwrite pre-existing cube files with new files, use the overwrite_existing_cube_files keyword.

This keyword supports spin-orbit coupling, but only when the type is either eigenstate_density or eigenstate. The large-scale use_local_index and load_balancing keywords are only supported when the type is either eigenstate_density or eigenstate.

Before we move on to the supported options for type, as well as other keywords related to cube output, there is an important note about specifying the dimensions of the cube that all users should know.

If no cube grid spacings are specified using the output cube edge tag, FHI-aims will use its own internal default for this grid spacing. In some cases, such as separated molecules or surfaces with large amounts of vacuum, the cube dimensions that FHI-aims would silently default to may not be ideal and may result in excessively large files.

In order to prevent uncontrolled damage (such as, filling up a file system or quota to beyond any reasonable limits) the code will therefore stop if a default number of cube grid points would be written to a single file in excess of limiting value defined by keyword cube_default_size_safeguard. This limit is configurable (see the description of that keyword).

If any output cube edge tag is specified in control.in, the limit given by cube_default_size_safeguard does not apply.

Also, many viewers do not implement non-rectangular cube edges, which would result from non-rectangular unit cells by default.

In short: Users are always strongly encouraged to specify cube geometries directly.

A list of all keywords related to cube plotting is given below. After this list of keywords, we have provided an example set of lines for plotting the total density of a molecule as well as the densities for individual eigenstates. This example should be adaptable for other types of cube files. Units for densities and eigenstates are Å^-3 and Å^-3/2, respectively. The unit for the long range and hartree potential is Hartree [Ha]. However, for electronic density response the unit is VÅ^-2.

• •

  output cube type This is the only mandatory line, specifying which type of cube file should be produced. FHI-aims presently allows the following options for type:

  1. 1.

     delta_density : Writes the difference between the initial (superposition of free atoms) and the final self-consistent density to a file.

  2. 2.

     eigenstate_density $n$. Writes the electron density of the $n-th$ eigenstate to a file. One may put string $homo$ or $lumo$ here instead of an int. In this case, the eigenstate index and kpt will be calculated based on the SCF occupation result. The eigenstate density is obtained as the square of the wavefunction. In periodic calculations and spin-orbit coupled calculations, this includes the contribution from the imaginary part of the wave function, and thus the output of this type is more physically relevant than the eigenstate type, which outputs only the real part of the wave function. By default, the first spin channel and the first k-point is printed out in non-spin-orbit coupled calculations, see also cube spinstate. In the spin-orbit coupled case, electron density contributed from both spin channels and the first k-point is printed out.

  3. 3.

     eigenstate $n$ : Writes the real part of the wave function of the $n$-th eigenstate to a file. $n$ must be an integer number. For non-periodic non-spin-orbit coupled calculations, the wave function has no imaginary part, so this type is sufficient. However, periodic and/or spin-orbit-coupled calculations produce wave functions with both a real part and an imaginary part. The corresponding imaginary part can be printed by requesting the cube type eigenstate_imag (see below) for the same state, in addition to the eigenstate type. Alternatively, one may eigenstate_density type for periodic calculations (and non-periodic calculations with spin-orbit coupling) instead. By default, the first spin channel and the first k-point is printed out, see also cube spinstate and cube kpoint.

  4. 4.

     eigenstate_imag $n$ : If applicable (for example, in periodic or spin-orbit coupled calculations), writes the imaginary part of the of the $n$-th eigenstate to a file. $n$ must be an integer number. Specifically for a spin-orbit coupled orbital, it is important to remember that such an orbital is a vector of two complex-valued functions:

     $${\psi }_{n,\mathrm{k}}(𝐫)=\left(\begin{array}{c}{\psi }_{n,𝐤}^{(1)}(𝐫)\\ {\psi }_{n,𝐤}^{(2)}(𝐫)\end{array}\right)$$

     (3.278)

     The two functions ${\psi }_{n,𝐤}^{(1)}(𝐫)$ and ${\psi }_{n,𝐤}^{(2)}(𝐫)$ are complex-valued, i.e., they each have a real part and an imaginary part. Thus, a total of four scalar functions $\mathrm{Re}({\psi }_{n,𝐤}^{(1)})$, $\mathrm{Im}({\psi }_{n,𝐤}^{(1)})$, $\mathrm{Re}({\psi }_{n,𝐤}^{(2)})$, $\mathrm{Im}({\psi }_{n,𝐤}^{(2)})$ need to be plotted to get the full orbital ($n$,$\mathrm{k}$). The two vector components (1) and (2) are sometimes loosely called “spin channels” although they are spin channels only in the non-relativistic limit; in the actual spin-orbit coupled case, the expectation value of the Pauli matrices would have to be calculated to get the spin direction. In any case, to get a full spin-orbit coupled orbital (say, $n$=1568 and the default $k$-point, usually – for unshifted $k$-point grids – $\mathrm{\Gamma }$) as cube file output, the syntax to use in control.in is this:

     Copied!

     
           output cube eigenstate 1568
           cube spinstate 1
     
           output cube eigenstate_imag 1568
           cube spinstate 1
     
           output cube eigenstate 1568
           cube spinstate 2
     
           output cube eigenstate_imag 1568
           cube spinstate 2
     
         
     

  5. 5.

     spin_density : The spin density $n^{\uparrow }(𝐫)-n^{\downarrow }(𝐫)$ is written to a file named spin_density.cube. Only available for spin collinear.

  6. 6.

     basis_function $n$ : Writes the real part of the $n$-th basis function to a file. $n$ must be an integer number and no larger than the number of basis functions.

  7. 7.

     stm : Must be followed by a real number $V$. Calculates 3D tunneling current map (more precisely, the tunneling current is proportional to the printed values) which can be used to plot STM images for a given voltage $V$ (in Volts) in the frame of the Tersoff-Hamann model. This is done by summing up eigenstate densities for all eigenstates between the Fermi level and $V$ (in eV), and the result is multiplied by $V$. In addition, a file cube_xxx_stm_z_map.cube will be printed. It contains values of the z-coordinate at the vertices of the cube, and can be used along with the tunneling current map to color the constant current isosurfaces according to their extent in the z-direction (to mimic the constant current image contrast in STM imaging). The output of stm-cubes works only for periodic systems.

  8. 8.

     ri_density : Resolution-of-the-Identity (RI) coefficients are read from file "ri_restart_coeffs.out" and used to evaluate the corresponding electronic density on the cube grid using the RI basis functions. Can be used in conjunction with keywords ri_restart and ri_full_output to visualize the RI approximation of the density. See Lewis et. al, JCTC, 17, 11 (2021), DOI: 10.1021/acs.jctc.1c00576 for more details.

  9. 9.

     total_density : The full electron density is recomputed and written to the a. In case of a periodic calculation, electron density from all unit cells that overlap with the cube output region will be printed.

  10. 10.

      total_density_integrable : The full electron density is recomputed and written to the a. In case of a periodic calculation, electron density from all unit cells that overlap with the cube output region will be printed. The output value is the integration in the voxel instead of the value at the center of the voxel. The total density here is more suitable for functions like Bader analysis.

  11. 11.

      long_range_potential : Prints the long range electrostatic potential of the Ewald summation. This result is useful in regions where no electron density is found and is much faster than the output of the full potential.

  12. 12.

      hartree_potential: The whole (i.e, short-range and long-range) electrostatic potential is recomputed on a cube grid and written out. Be careful with keyword potential. Please report errors.

  13. 13.

      xc_potential: The xc potential is recomputed on a cube grid and written out. PBE only, spin unpolarized only.

  14. 14.

      potential: Legacy. The whole (i.e, short-range and long-range) electrostatic potential is recomputed on a cube grid and written out. The keyword is considered broken/experimental. Use hartree_potential.

  15. 15.

      delta_v: Output $\delta v=v-v^{\mathrm{free}}$. This is especially tested for MPB solvation effects and is also an experimental feature especially for vacuum calculations.

  16. 16.

      ion_dens: Ionic charge density $n_{\mathrm{ion}}^{\mathrm{MPB}}$ as obtained from an MPB-DFT calculation. Still experimental.

  17. 17.

      dielec_func: Dielectric function $\epsilon [n_{\mathrm{el}}]$ as used in the MPB-DFT calculation.

  18. 18.

      elf: Electron localization function. Different options are available for spin-polarized systems, see keyword cube elf_type. Currently under testing, please report any errors. The implementation is not yet compatible with spin-orbit calculations.

  19. 19.

      first_order_density $n$: Density response with respect to an applied homogeneous electric field is printed in a cube file ($\frac{\partial \rho }{\partial {ϵ}_n}$). So that, n represents the direction of the field (1,2 or 3). a. For nonperiodic systems, DFPT polarizability should be specified in control.in before the cube file commands. b. For periodic system, DFPT dielectric should be specified in control.in before the cube file commands. Note that as of version stamp 241126 this keyword is only possible with the old DFPT interface, by setting dfpt_centralised .false.

      An example outputting density response with respect to the three Cartesian direction in control.in is this:

      Copied!

      
      output cube first_order_density 1
          cube origin       7.3408      7.6288     52.3975
          cube edge 100      0.1468      0.0000      0.0000
          cube edge 100      0.0000      0.1526      0.0000
          cube edge 300      0.0000      0.0000      0.1480
      output cube first_order_density 2
          cube origin       7.3408      7.6288     52.3975
          cube edge 100      0.1468      0.0000      0.0000
          cube edge 100      0.0000      0.1526      0.0000
          cube edge 300      0.0000      0.0000      0.1480
      output cube first_order_density 3
          cube origin       7.3408      7.6288     52.3975
          cube edge 100      0.1468      0.0000      0.0000
          cube edge 100      0.0000      0.1526      0.0000
          cube edge 300      0.0000      0.0000      0.1480
      
      
          
      

      Hence, three cube file will be outputted for each direction.

• •

  cube spinstate spin
  This keyword allows the user to choose whether to print spin-channel 1 or 2. The default value is 1. An example for how to use the cube spinstate keyword to print eigenstates in both spin channels is given below. This keyword is only useful for spin-polarized calculations (keywordspin collinear) and for spin-orbit coupled calculations. Otherwise, scalar-relativistic spin-non-polarized calculations (keywordspin none) have degenerate spin channels by definition. See the description of eigenstate_imag for an example of how to write all parts of a spin-orbit coupled orbital to a cube file.

• •

  cube kpoint kpoint
  This keyword allows to choose the k-point to be printed. $kpoint$ is an integer number following the same ordering as the k-points within the SCF-cycle. It is presently not possible to output cube files at a k-point not included in the scf. Keyword output k_point_list may be used to print out the entire list of k-points used in the s.c.f. cycle. This will help identify which k-point is printed in a cube file. Default: 1

• •

  cube state spin k-point
  Deprecated keyword to choose spin and k-point. The cube spinstate keyword should be used instead as given in the example below. If nothing else is specified, this keyword defaults to cube state $1$ $1$. Note that for cluster calculations, $k-point$ must always be 1.

• •

  cube filename name_of_the_file
  Allows to customize the name of the cubefile. If this line is not given, FHI-AIMS will default to a file name which contains the number of the cube requested, its type, and, if applicable, the corresponding spin and k-point of the data.

• •

  cube format format
  Apart from the default cube format, FHI-AIMS also supports output in the formats of the gOpenMol and XCrysden software. This is requested by the line cube format $format$. The options for $format$ are cube, gOpenMol, and xsf (the XCrysden format). Default: cube

• •

  cube divisor number
  This is a technical settings which governs the paralellization of the cube output. The whole cube is divided into smaller, so-called minicubes, which are then treated one after the other. The value governs the number of points in each directed to be used for each minicube, i.e., its size. A larger setting therefore results in less minicubes (and it thus potentially faster), but also a larger demand for memory. Unless there are problems with memory, this setting usually does not need to be touched, with the exception of output potential, where it should be set to its maximal value (45), independent of the number of processors used. Default: 10

• •

  cube spinmask $i$ $j$
  For total_density cube files, the line cube spinmask $i$ $j$ with integer $i$ and $j$ allows to manipulate the spin channels independently according to the following formula: $i\cdot n^{\uparrow }(𝐫)-j\cdot n^{\downarrow }(𝐫)$. If $i$ = 1 and $j$ = 1 (which is the default), the total density will be computed as normal. This keyword replaces the earlier keyword cube spin .

• •

  cube elf_type $i$
  Specifies the type of electron localization function (ELF) to be calculated. The default $i=0$ (and the only available option for spin-unpolarized calculations) is the Savin et al. formula [269]. For spin-polarized systems, $i=1$ and $i=2$ correspond to the original formulation by Becke and Edgecombe [25] for spin channels 1 and 2, respectively. If $i$ is not 0, 1, or 2, and the system is spin-polarized, the Kohout-Savin variant of ELF [173] will be calculated.

• •

  cube origin $x$ $y$ $z$
  Single line which specifies the origin , i.e., the center of the region to be plotted. Values are given in ÅIf omitted, the same origin as for the previous cube file is used. If no origin has been given yet, it defaults to the geometric center of the molecule in the cluster case or (0,0,0) for periodic calculations. For slab type calculations (when using dipole_correction .true.), the origin is set to the center of the slab.

• •

  cube edge_density $den$
  Speciefies the edge_density Å for the cube file, the number of edges in each direction is then calculated by $\frac{latticelength}{edge\mathrm{\_}density}$.

• •

  cube edge $n$ $dx$ $dy$ $dz$
  Specifies the edges of the volumetric data to be plotted. Separate lines have to be given for each of the three edges of the cube. In each line, $n$ indicates the number of steps a particular edge (voxel), and $dx$, $dy$, $dz$ indicate the length of each individual step [i.e., the full cube edge length is ($n\cdot dx$, $n\cdot dy$, $n\cdot dz$). If omitted, the same edges as for the previous cube file are used. If no edges have been specified yet, FHI-aims defaults (in the cluster case) to orthogonal grids of 0.1 Å length, which span the whole molecule plus 14 Bohr beyond the outermost nucleii. For periodic calculations, the default edges are the same as the lattice vectors, again with 0.1 Å step length.

  Note: In some cases, such as separated molecules or surfaces with large amounts of vacuum, the defaults might be far from ideal and result in excessively large files. In order to prevent uncontrolled damage (such as, filling up a file system or quota to beyond any reasonable limits) the code will therefore stop if a default number of cube grid points would be written to a single file in excess of limiting value defined by keyword cube_default_size_safeguard. This limit is configurable (see the description of that keyword).

  Users are always strongly encouraged to specify cube geometries themselves.

NOTE: A word of warning here: FHI-aims outputs the density and the wave function in units of 1/ų and 1/Å^3/2, respectively, by default. In contrast, the length units of the cube voxels are in atomic units. The choice of units must be accounted for when doing any kind of postprocessing.

The keyword cube_content_unit allows one to switch the output units of the printed array.
cube_content_unit legacy is the default and gives the output in SI units, while
cube_content_unit bohr switches to output of densities and orbitals in atomic units. (Note that other codes’ definition of cube file usually is in atomic units, which is why we do not provide a switch to Å).

As promised, an example set of lines for the control.in file showing how to plot the total density and densities for individual eigenstates of a hypothetical spin-polarized non-periodic system. This example will not exactly correspond to your system and should not be blindly copy-pasted. Nevertheless, it should give a good idea of how the keywords presented previously work together.

    output cube total_density
    cube origin 1.59 9.85 12.80
    cube edge 101 0.15 0.0 0.0
    cube edge 101 0.0 0.15 0.0
    cube edge 101 0.0 0.0 0.15
    output cube eigenstate_density 151
    cube spinstate 1
    output cube eigenstate_density 151
    cube spinstate 2
    output cube eigenstate_density 152
    cube spinstate 1
    output cube eigenstate_density 152
    cube spinstate 2
    output cube eigenstate_density 153
    cube spinstate 1
    output cube eigenstate_density 153
    cube spinstate 2
    output cube eigenstate_density 154
    cube spinstate 1
    output cube eigenstate_density 154
    cube spinstate 2

The first line requests cube output for the total density. See the details above on how to get individual versions of the spin density. Then, the center of the cube is specified at (1.59, 9.5, 12.80) Å. Each cube direction has 101 points that are spaced apart by 0.15 Å, giving a total cube edge length of 15 Å in each direction. Finally, output is requested for the densities correspond to the Kohn-Sham wave functions associated with eigenstates number 151-154. For each eigenstate, the additional cube spinstate i line requests first the spin-up channel ($i$=1), then spin-down ($i$=2).

NOTE that for periodic systems, the cube spinstate lines may have to be followed by cube kpoint lines to specify the k-point at which an eigenstate is printed. By default, only $k$-point number 1 is printed. For unshifted $k$ grids, that is the $\mathrm{\Gamma }$ point.

Thus ends our brief treatise on cube plotting. We return you to your regularly scheduled programming.

Note that this density output is not given on a uniform grid, but simply on the overlapping atom-centered grid used for all internal operations of FHI-aims. For a density on a uniform grid, see the output cube subkeyword.

output density additionally writes the difference between the current density and the superposition-of-free-atom reference density to a file diff-density.dat

This option serves as interface to the DGrid program (written by Miroslav Kohout), which is available free of charge at https://www2.cpfs.mpg.de/~kohout/dgrid.html. DGrid allows to employ various electronic structure and chemical bonding analysis alogorithms in real space, e.g., the QTAIM (Quantum Theory of Atoms in Molecules of R.F.W. Bader) or the ELI-D/ELF method. After installation of DGrid (Version $\ge 5.0$) the command dgrid dgrid_aims.dat converts the interface wave function file dgrid_aims.dat into a new file dgrid_aims.fhi with native DGrid file format. This file provides the basis to all capabilities of the DGrid program described in detail in its manual at https://www2.cpfs.mpg.de/~kohout/Documents/dgrid-html/dgrid.html.

This is generally useful, but the dipole moment is particularly needed to compute molecular oscillator strengths for individual vibrational frequencies.

Note that, for charged systems, the electrical dipole is defined with respect to the (possibly arbitrary) origin of the file geometry.in, rather than an origin within the system itself. Thus, charged systems will yield different dipole moments for different choices of origin, but the important dipole differences needed to compute, e.g., oscillator strengths via finite diffences remain well-defined.

This keyword shares its syntax with output atom_proj_dos and output species_proj_dos. See also section 4.4 for more details.

Two output files emerge from this option:

• •

  KS_DOS_total_raw.dat contains the total DOS components as a function of the eigenvalue energy (first column) as used internally in FHI-aims. The energy zero is then given by the vacuum level (non-periodic systems) or by the $𝐆$=0 component of the long-range Hartree potential (periodic systems).

• •

  KS_DOS_total.dat contains the same information, except that the energy zero is shifted to the Fermi energy (metallic systems) or valence band maximum (insulators), respectively.

When followed by the following keyword

dos_kgrid_factors $n_1$ $n_2$ $n_3$

where $n_1$, $n_2$ and $n_3$ are integers, the dimensions of the k-point grid along the first, second and third lattice vectors are multiplied by $n_1$, $n_2$ and $n_3$, respectively. New eigenvalues are re-calculated non-selfconsistently on the new denser k-point grid. The new eigenvalues are then used to plot an improved (so-called perturbative) density of states.

If no dos_kgrid_factors are specified, the original k-point grid is used.

The unit of output for the density of states is ${(eV\cdot V_{\text{unit cell}})}^{-1}$, i.e., number of states per energy unit (eV) and unit cell volume.

This keyword supports spin-orbit coupling. When spin-orbit coupling is enabled, the file(s) containing the spin-orbit-coupled DOS will have the default filename(s) and the file(s) containing the scalar-relativistic (i.e. no SOC) band structure will have an additional suffix “.no_soc”. Note that if you requested spin collinear in the control.in in addition, the file with spin-orbit coupling included (suffix “.dat”) has only one column containing the dos of all spin-coupled states. The reason is that the separate spin channels of scalar relativity no longer exist after the SOC operator is applied – the states now form a single set. The output of the total DOS without spin-orbit coupling (suffix “.dat.no_soc”) has two columns - one for each spin channel.

This keyword shares its syntax with output atom_proj_dos_tetrahedron and output species_proj_dos_tetrahedron. See also section 4.4 for more details.

Two output files emerge from this option:

• •

  KS_DOS_total_raw_tetrahedron.dat contains the total DOS components as a function of the eigenvalue energy (first column) as used internally in FHI-aims. The energy zero is then given by the vacuum level (non-periodic systems) or by the $𝐆$=0 component of the long-range Hartree potential (periodic systems).

• •

  KS_DOS_total_tetrahedron.dat contains the same information, except that the energy zero is shifted to the Fermi energy (metallic systems) or valence band maximum (insulators), respectively.

The keyword dos_kgrid_factors is now supported by the tetrahedron method.

Two output files emerge from this option:

• •

  KS_DOS_total_raw.dat contains the total DOS components as a function of the eigenvalue energy (first column) as used internally in FHI-aims. The energy zero is then given by the vacuum level (non-periodic systems) or by the $𝐆$=0 component of the long-range Hartree potential (periodic systems).

• •

  KS_DOS_total.dat contains the same information, except that the energy zero is shifted to the Fermi energy (metallic systems) or valence band maximum (insulators), respectively.

The unit of output for the density of states is ${(eV\cdot V_{\text{unit cell}})}^{-1}$, i.e., number of states per energy unit (eV) and unit cell volume.

This keyword supports spin-orbit coupling. When spin-orbit coupling is enabled, the file(s) containing the spin-orbit-coupled DOS will have the default filename(s) and the file(s) containing the scalar-relativistic (i.e. no SOC) band structure will have an additional suffix “.no_soc”. Note that if you requested spin collinear in the control.in file in addition, the file with spin-orbit coupling included (suffix “.dat”) has only one column containing the dos of all spin-coupled states. The reason is that the separate spin channels of scalar relativity no longer exist after the SOC operator is applied – the states now form a single set. The output of the total DOS without spin-orbit coupling (suffix “.dat.no_soc”) has two columns - one for each spin channel.

If this keyword is set, the following information for each eigenstate is written to file: index, Kohn-Sham state, spin, k-point, k-weight, electronic occupation, eigenvalue (in eV). This data is printed one row per eigenstate. Crucially, data is only printed for eigenstates in the k-point-symmetrised set, if applicable (for instance when the system is periodic and k-point symmetry reduction is being used).

For non-periodic systems, this option causes the Kohn-Sham eigenvectors $c_{il}$ (for basis functions $i$, eigenstates $l$) to be written out for each state, whenever the Kohn-Sham eigenvalues are written. However, at present the non-periodic version is mostly of ’debug’ character. In particular, KS_method parallel is not supported under all circumstances. Please test carefully.

For periodic systems, eigenvector output must be requested together with the output band functionality described above. If output eigenvectors is set in addition, the Kohn-Sham eigenvectors $c_{il}(𝒌)$ will be written into separate files for each spin channel, only for each $k$-point for which band output was requested.

For each state $l$, the real and imaginary part of $c_{il}(𝒌)$ will then be written out for each basis function $i$. Output file names are assigned automatically based on the number of the output band, and to the specific k-point in that band, to which they pertain, e.g.:
KS_eigenvectors.band_$\ast$.kpt_$\ast$.out
In addition to the eigenvectors themselves, these files also contain as header information:

• •

  the relative coordinates of the $k$-point in question (in units of the reciprocal lattice vectors of the structure in question)
  

• •

  Information on the identity (atom number and angular momentum) for each basis function
  

• •

  The Kohn-Sham eigenvalue and occupation number of each state.
  

The listed eigenvectors pertain to the superimposed, Bloch-like basis functions (with phase factors!) ${\chi }_{i,k}(𝒓)$ as defined through Eq. (22) of the FHI-aims Computer Physics Communications description, Ref. [36].

Note that the output aitranss keyword provides one further option to print eigenvectors and other information for any non-periodic system, serial or parallel (ScaLAPACK, ELPA, ELSI).

This keyword will output scalar-relativistic eigenvectors. For spin-orbit-coupled eigenvectors, please see the soc_eigenvectors keyword.

When this flag is enabled, the ELectronic Structure Infrastructure (ELSI) will output information every time it is used to solve the Kohn-Sham eigenvalue problem, including timings for the ELSI invocation, a list of values for relevant variables, and versioning information.

This information is stored in a JSON format and is written to the file elsi_log.json using the FortJSON library, which is distributed with FHI-aims and ELSI and is built automatically. Only task 0 will output this file.

To output information from an FHI-aims calculation in a JSON format, please use the output json_log keyword.

For periodic systems, this option allows to write out the k-dependent (Bloch) hamiltonian matrices that correspond to the set of k-points requested with the output band keyword. Both spin channels are written into the same file.

Defining “atoms-in-molecules” is a classic, intuition based problem; one would like to associate individual (partial) charges or moments with individual atoms in a bonded structure. This process is by necessity not uniquely definable (molecules are not atoms, and the are no rigorously defined boundaries between atoms). Nonetheless, much chemical intuition is based on such a concept.

Hirshfeld’s [144] “atoms-in-molecules” partitioning relies on the same idea as the partitioning of our charge density for the electrostatic potential (Eq. 3.30), using the free-atom electron density $n_{\text{at}}^{\text{free}}(r)$ as the weight function $g_{\text{at}}(𝐫)$ in Eq. (3.26). Since (for a given functional), we know the spherical $n_{\text{at}}^{\text{free}}(r)$ exactly, this analysis remains well-defined even as external circumstances such as the basis set change. However, the resulting values are still qualitative in the sense that Hirshfeld’s [144] “atoms-in-molecules” partitioning is just one among many other prescriptions that have been suggested in the literature. Any ghost atoms in the system are skipped for this analysis. This could lead to “missing” charge if the ghost atoms carried any nonnegligible charge, but that should happen only in abnormal systems.

Note that the output hirshfeld keyword itself only writes a Hirshfeld analysis for the final geometry of an FHI-aims run, not, for instance, for intermediate molecular dynamics steps. This ensures that the Hirshfeld analysis does not accidentally clutter an output file with large amounts of data if the output_level MD_light output level is set. Output for every geometry can be accomplished with the output hirshfeld_always keyword below.

This keyword ensures that a Hirshfeld analysis is written at every geometry step of a FHI-aims run, not just after the final step. If output hirshfeld-I is requested together with output hirshfeld_always, results from both the normal and the iterative Hirshfeld analysis are written at every step.

Similar functionality to the normal output hirshfeld Hirshfeld analysis – see the comments for that keyword – except that in the “iterative Hirshfeld” [47] analysis the partition weights are changed. Here, the partitioning densities are not those of neutral atoms but rather those of ions with the same formal charge as the formal charge determined by the “iterative Hirshfeld” analysis.

This keyword is implemented but has not seen much production use. It is therefore not guaranteed that it will always work or that the results will always make sense or even be in line with the original “iterative Hirshfeld” publication.[47] All may be well, but if you do use the functionality, please check very carefully that the results appear to be correct.

When this flag is enabled, FHI-aims will output information in a JSON format, which may be easily parsed by your favorite post-processing language (or Python).

This feature was added in 2018, and much of the functionality in FHI-aims outside of the main SCF cycle will not write out any information. A partial list of quantities which will be written includes:

• •

  Initial/final geometries

• •

  Versioning information

• •

  Runtime settings: number of basis functions, number of k-points, etc.

• •

  SCF convergence evolution

• •

  Mulliken decompositions (when calculated)

• •

  Domain decomposition (a.k.a. batch partitioning) statistics

• •

  Total energies

• •

  HOMO/LUMO levels (on the SCF k-grid)

• •

  Fundamental gap (on the SCF k-grid)

• •

  Timings

This information is written to the file aims.json using the FortJSON library, which is distributed with FHI-aims and built automatically. Only task 0 will output this file. Unlike the ELSI JSON log (which is written out using the output elsi_log keyword), this JSON log will persist through SCF reinitialization, geometry relaxation steps, and MD steps; it will only be overwritten when a new FHI-aims calculation is performed.

To output information from ELSI directly in a JSON format, please use the output elsi_log keyword.

Eigenvalues will only be written for the first number $k$-points by default. For dense $𝐤$-grids, the sheer number of $k$-points simply gets too large to allow for a full output.

The $k$-point coordinates are written in units of the reciprocal lattice vectors of the structure. This option is also the default when using output_level full and KS_method serial.

This subkeyword provides the file-based interface to the LibRPA package (https://srlive1201.github.io/LibRPA) for efficient Green-function-based calculations. Note that it is available only for periodic calculations. For usage examples of this keyword, please see Sec. 3.34.

Setting this keyword is equivalent to setting elsi_output_matrix hamiltonian and elsi_output_matrix overlap. Note that this keyword outputs matrices in the $k$-space, not in real space. Therefore, the number of Hamiltonian matrix files is equal to the number of $k$-points, multiplied by the number of spin channels. The overlap matrix, however, is only written out for the first spin channel to save disk space. Additionally, the Hamiltonian is written out in every SCF iteration, while the overlap is only written out in the first SCF iteration. The output files are in the ELSI CSC binary format (see the documentation of ELSI). The Python scripts in the “utilities/elsi_matrix” directory may be used to post-process an ELSI matrix file, e.g., converting it to a human-readable format.

Depending on option “types”, the upper part of up to three different matrices is written into separate files. If types is “n” (without quotes), no matrices are written. If types is set to “h”, then the Hamiltonian is written. The choice “o” causes the overlap matrix to be output. Finally, “s” refers to the system matrix from which the eigenvalues are calculated. The last three options can be combined. For example, “ho” means Hamiltonian and overlap.

The format of the files can be controlled with the optional parameter format. Possible values are “asc” for ASCII output (default) and “bin” for binary output.

Note that the output is written as binary files in the ELSI CSC format (see the documentation of ELSI). The Python script in the “utilities/elsi_matrix” directory may be used to convert an ELSI matrix file to a human-readable format.

Defining “atoms-in-molecules” is a classic, intuition based problem; one would like to associate individual (partial) charges or moments with individual atoms in a bonded structure. This process is by necessity not uniquely definable (molecules are not atoms, and the are no rigorously defined boundaries between atoms). Nonetheless, much chemical intuition is based on such a concept.

A classic “atoms-in-molecules” concept is the Mulliken analysis [225], which defines electronic occupations of atoms by projected occupations into the localized basis functions associated with them (see the standard literature for exact definitions and use).

In short, when so requested, FHI-aims will provide a decomposition of the electronic density per atom, angular momentum channel, and possibly spin channel, allowing to deduce approximate partial charges. The summarized Mulliken analysis is written into the standard output stream, while a separate file Mulliken.out contains detailed state-by-state projected electron occupations.

Note that a Mulliken analysis is somewhat ill-defined because of basis function overlap; thus electrons can be counted to one atom or another at will. For small basis sets, a Mulliken analysis may still yield qualitatively reasonable numbers, but it becomes increasingly ill-defined as the atom-centered basis sets approach basis set completeness.

This keyword supports spin-orbit coupling. When spin-orbit coupling is enabled, the file(s) containing the spin-orbit-coupled Mulliken analysis will have the default filename(s) and the file(s) containing the scalar-relativistic (i.e. no SOC) Mulliken analysis will have an additional suffix ".no_soc".

As above, except here the supplementary Mulliken.out is not generated, which can result in substantial savings on filespace.

This can be useful if the FHI-aims basis functions are needed for a further, separate purpose (Quantum Monte Carlo etc) but please note that the integration accuracy for the Coulomb singularity near the nuclei must be higher than in our standard calculations (where the singularity is cancelled by the kinetic energy), so increasing radial_multiplier is in order.

Since August 2013, FHI-aims verifies the accuracy of its ’radial’ integration grid (the sparse grid of atom-centered radial shells around each atom which is part of the definition of its three-dimensional, overlapping atom-centered integration grids) in comparison to integrals on the dense ’logarithmic’ grid which is used to set up the spherical free atom, all radial functions etc. in one dimension.

These integrals take the form

With our usual definition of basis functions,

we get:

in the non-relativistic case. In the case of scaled ZORA or atomic ZORA scalar relativity, the kinetic energy part is modified and the integrand reads:

These are the integrands to test both the ’logarithmic’ grid and the ’radial’ grid around each atom, where $v(r)$ is set to be the one-dimensional potential of a spherical free atom as calculated at the outset of each run.

If output onsite_integrands is set to be true, the actual integrands $f(r)$ and various of their parts are printed for each radial function. This is mainly useful for debugging purposes, to understand what we are integrating for a given basis function choice. Especially for contracted Gaussian basis functions, $f(r)$ can look quite ugly near the nucleus.

The files that contain the actual integrand $f(r)$ defined above are called
Onsite_r2_phi_h_phi_log.(function).dat and
Onsite_r2_phi_h_phi_rad.(function).dat for the logarithmic and radial grids, respectively, with “(function)” indicating the element and the radial function number in the order used by FHI-aims (for instance, the output basis keyword uses the same order to output the radial functions used in the code). Units are Å for the radial coordinate, but atomic units (Ha/bohr³) for the integrand itself.

For periodic systems, this option allows to write out the k-dependent (Bloch) overlap matrices that correspond to the set of k-points requested with the output band keyword. If output eigenvectors is set in addition, the Kohn-Sham eigenvectors $c_{il}(𝒌)$ will be written into separate files for each spin channel, only for each $k$-point for which band output was requested.

This option can help show if (or if not) the chosen basis set for the full system is close to ill-conditioning (see the comments for keyword basis_threshold).

If this keyword is set, the eigenvalues and occupation numbers for a periodic system are recomputed and written to a file ’Final_KS_eigenvalues.dat’ after the s.c.f. calculation (and, possibly, relaxation, dynamics etc.) is complete. A denser $k$-space grid than during the original s.c.f. calculation can be specified using the dos_kgrid_factors keyword. When combined with the dos_kgrid_factors keyword, the ordering of the k-points may not be in the most logical manner. Similarly, when combined with symmetry_reduced_k_grid .true. there will be a different ordering of the k-points to when symmetry_reduced_k_grid is .false.. If you wish to have a more logical ordering of the k-points, please see the script sort_postscf_eigenvalues.py in the utilities folder.

Note that the resulting output file can become very large. See the header of the ’Final_KS_eigenvalues.dat’ for details and for units used.

Finally, note that an externally generated k-point list k_list.in (see keyword
k_points_external) is not supported by output postscf_eigenvalues and an internally generated, even-spaced k-point grid (also defined in k_list.in) is used instead.

This keyword must be used together with include_spin_orbit and output band or band_mulliken.

FHI-aims normally solves the non-selfconsistent spin-orbit coupled eigenvalue problem (including eigenvectors) in a basis of self-consistent eigenvectors calculated in a preceding, scalar-relativistic calculation (no spin-orbit coupling). This means that the spin-orbit coupled eigenvectors are immediately available and can be interpreted in terms of their origin from an underlying scalar-relativistic basis set.

The spin-orbit coupled eigenvectors in terms of their non-spin-orbit-coupled counterparts will be written for each sampled k point in each specified band in control.in.

As an example, assume that there are $n$ scalar-relativistic (i.e.,pre-SOC) eigenstates in a closed-shell (non-spinpolarized) self-consistent calculation. After the SOC perturbation is applied and the resulting eigenvalue problem solved, this leads to 2$n$ spin-orbit coupled eigenvectors. These are the columns of a (2$n\times$2$n$) matrix, where the first $n$ rows represent $n$ pre-SOC states of the scalar-relativistic spin channel 1 and the following $n$ rows (row numbers ($n$+1) to 2$n$) represent the $n$ pre-SOC states of the scalar-relativistic spin channel 2. Since the original scalar-relativistic eigenvectors form an orthonormal basis set, the absolute value of element ($i$,$j$) represents the contribution from the pre-SOC state corresponding to the $i$th row to the post-SOC state corresponding to the $j$th column.

Note that for parallel calculations, the output is written as binary files, using the ELSI infrastructure (for details, the documentation of ELSI covers this format). You may want to use the script convert_elsi_to_mm.py together with read_elsi.py in the utilities/elsi_matrix directory to create a human-readable format.

This option is based on a Mulliken Analysis and shares its syntax with output dos and output atom_proj_dos. See also section 4.4 for more details.

Different from the atom_proj_dos option, this option adds up the pDOS contributions of all atoms of each species defined in control.in and used in geometry.in. This provides a quick handle to obtain the pDOS contribution of well-defined subgroups of individual atoms, e.g., those of a given layer of a slab, by simply defining a separate species for the desired atoms in control.in.

There are two types of output files for each atom:

• •

  species_l_proj_dos_raw.dat, where species denotes the species name used in geometry.in and control.in. This file contains the total and angular-momentum resolved DOS components as a function of eigenvalue energy (first column) as used internally in FHI-aims. The energy zero is then given by the vacuum level (non-periodic systems) or by the $𝐆$=0 component of the long-range Hartree potential (periodic systems).

• •

  species_l_proj_dos.dat, which gives the same information, except that the energy zero is shifted to the Fermi energy (metallic systems) or valence band maximum (insulators), respectively.

Note that projected densities of states such as given here must be based on some kind of projection orbitals, the choice of which is somewhat arbitrary by necessity. This is thus a tool for qualitative analyses.

In FHI-aims, we project directly on the atom-centered angular-momentum components as defined by the overlapping basis set. This definition becomes the more arbitrary te larger the basis set, just like a mulliken analysis. The closer the full basis comes to completeness, the more ambiguous will a mulliken-like analysis become, since it may not be a priori clear which electrons should be counted towards the basis functions of one atom vs. those of another atom. Thus, do not expect a pDOS to simply converge as the basis set size is increased; use it as a qualitative indicator of trends, but nothing more.

This keyword supports spin-orbit coupling. When spin-orbit coupling is enabled, the file(s) containing the spin-orbit-coupled DOS will have the default filename(s) and the file(s) containing the scalar-relativistic (i.e. no SOC) DOS will have an additional suffix “.no_soc”. Note that if you requested spin collinear in the control.in file, there will be only one file per species (and equivalently for all _raw.dat files) containing the projected DOS of all spin-coupled states:

• •

  species_l_proj_dos.dat

The reason is that the separate spin channels of scalar relativity no longer exist after the SOC operator is applied – the states now form a single set. However, there are two DOS files for the output without spin-orbit coupling; one for each spin channel:

• •

  species_l_proj_dos_spin_up.dat.no_soc

• •

  species_l_proj_dos_spin_dn.dat.no_soc.

As of version 240605 this keyword can now be used in combination with
dos_kgrid_factors. The n_broadening parameter was also introduced in this version.

When this keyword is combined with atom_proj_dos_tetrahedron,
species_proj_dos_tetrahedron, or mulliken, and dos_kgrid_factors is not used, then the n_broadening parameter will be ignored. This is because a Mulliken decomposition over all states is performed for each of these other keywords, which is then re-used for the projected DOS with species_proj_dos.

When this keyword is combined with atom_proj_dos, and different values of Estart, broadening, n_broadening, or Eend are used for each, the code will take the most extreme limit from the values of Estart-broadening*n_broadening and Eend + broadening*n_broadening for the atom and species projected DOS, to ensure that a single Mulliken decomposition is done that utilizes states to encompass both projections, such that two separate decompositions do not need to be carried out.

This option is based on a Mulliken Analysis and shares its syntax with output dos_tetrahedron and output atom_proj_dos_tetrahedron. See also section 4.4 for more details.

Different from the atom_proj_dos_tetrahedron option, this option adds up the pDOS contributions of all atoms of each species defined in control.in and used in geometry.in. This provides a quick handle to obtain the pDOS contribution of well-defined subgroups of individual atoms, e.g., those of a given layer of a slab, by simply defining a separate species for the desired atoms in control.in.

There are two types of output files for each atom:

• •

  species_l_proj_dos_tetrahedron_raw.dat, where species denotes the species name used in geometry.in and control.in. This file contains the total and angular-momentum resolved DOS components as a function of eigenvalue energy (first column) as used internally in FHI-aims. The energy zero is then given by the vacuum level (non-periodic systems) or by the $𝐆$=0 component of the long-range Hartree potential (periodic systems).

• •

  species_l_proj_dos_tetrahedron.dat, which gives the same information, except that the energy zero is shifted to the Fermi energy (metallic systems) or valence band maximum (insulators), respectively.

Note that projected densities of states such as given here must be based on some kind of projection orbitals, the choice of which is somewhat arbitrary by necessity. This is thus a tool for qualitative analyses.

In FHI-aims, we project directly on the atom-centered angular-momentum components as defined by the overlapping basis set. This definition becomes the more arbitrary te larger the basis set, just like a mulliken analysis. The closer the full basis comes to completeness, the more ambiguous will a mulliken-like analysis become, since it may not be a priori clear which electrons should be counted towards the basis functions of one atom vs. those of another atom. Thus, do not expect a pDOS to simply converge as the basis set size is increased; use it as a qualitative indicator of trends, but nothing more.

This keyword supports spin-orbit coupling. When spin-orbit coupling is enabled, the file(s) containing the spin-orbit-coupled DOS will have the default filename(s) and the file(s) containing the scalar-relativistic (i.e. no SOC) DOS will have an additional suffix ”.no_soc”. Note that if you requested spin collinear in the control.in file, there will be only one file per species (and equivalently for all _raw.dat files) containing the projected DOS of all spin-coupled states:

• •

  species_l_proj_dos_tetrahedron.dat

The reason is that the separate spin channels of scalar relativity no longer exist after the SOC operator is applied – the states now form a single set. However, there are two DOS files for the output without spin-orbit coupling; one for each spin channel:

• •

  species_l_proj_dos_tetrahedron_spin_up.dat.no_soc

• •

  species_l_proj_dos_tetrahedron_spin_dn.dat.no_soc.

Note that the meaning of this effective potential depends on the xc option used. For DFT-LDA, this is simply the full local potential. For gradient-corrected (GGA) functionals, the gradient partial derivatives of the exchange-correlation functional are not included in the potential, as they are treated separately by integration by parts (see Ref. [36]). For hybrid functionals or Hartree-Fock, the exchange part of the potential is of course not included.

Note also that this output does not happen on a uniform grid. For further processing, a proper visualization tool is needed, and/or an interpolation onto a uniform grid must be done.

The resulting values are “atoms-in-molecules” like partial charges assigned to each atom, similar to a hirshfeld partitioning (but not identical because a different partitioning function may be used as the hartree_partition_type).