B.1 Flow of the program

The uppermost level of FHI-aims is the subroutine main.f90, the structure of which is shown as a structogram in Fig. B.1. As a subroutine, main.f90 can also be called by external code as a library subroutine, with the restriction that, for parallel execution, this can happen only once [by definition of the message passing interface (MPI) for parallel communication, there can only be one call to mpi_init per program run].

The purpose of main.f90 is to provide a logical separation of groups of computational tasks by way of high-level wrapper subroutines (listed in typewriter font in Fig. B.1). With the exception of global convergence checks, no outright physical quantities are directly manipulated in main.f90. All physically relevant quantities are handled inside the lower-level structure of the code and, if necessary, are passed between them by way of specific modules. For example, the module physics.f90 handles all variables of tangible physical importance (Hamiltonian and overlap matrices, Kohn-Sham wave function, electron densities, potentials, …). The geometry information for a given electronic structure cycle (coordinates) are found in module geometry.f90; etc.

More details regarding these and other modules are included with the code distribution as a separate document. The point here is that main.f90 should never need to use any but the highest-level modules explicitly, i.e.:

• •

  dimensions.f90 : Inclusion of array dimensions for consistent allocations across the code, and wrapper flags to prevent access to unallocated variables which are not needed for a given task.

• •

  localorb_io.f90 : Wrapper module for consistent writing of output in parallel runs

• •

  mpi_utilities.f90 : Wrapper module for MPI initialization, finalization (first and last tasks of a run, respectively), and task distribution

• •

  timing.f90 : Wrapper module including all timing and accounting information, including the count of s.c.f. iterations, relaxation or MD steps.

The first task of main.f90 is to initialize any accounting (timing etc.) information, the infrastructure required for MPI (or, to silently switch off the use of MPI entirely in the case of non-parallel runs), and to record all this information (initial time stamps, code version, number of parallel tasks and computer system layout) in the standard output.

The obvious next task is to read and process all input information given in control.in and geometry.in. Internally, this is handled in three steps (see the wrapper subroutine read_input_data.f90): First, both input files are parsed once, while extracting only the dimension information needed to set up any necessary arrays / array dimensions needed to house the following input data. Organizing this information is the task of module dimensions.f90. Next, the information in control.in is read and checked for consistency, using read_control.f90 for all general information, and repeated calls to read_species_data.f90 for all species-related information. Finally, subroutine read_geo.f90 reads the input data of geometry.in, and verifies its consistency with the data contained in control.in.

At the end of this step (subroutine read_input_data.f90), all input data from all input files should have been read and processed. It is important that any known conflicts, or incomplete settings, should have been verified at this stage, stopping the code with an error message if outrightly conflicting input information is detected. For completeness, we mention that any technical input settings of global interest (e.g., the handling of spin, relativity, or exchange-correlation) are collected and accessible through the top-level module runtime_choices.f90.

The following steps are the “household” steps of electronic structure theory.

Wrapper subroutine prepare_scf.f90 sets up all structure-independent, fixed pieces of the calculation, and stores them for easy access in the actual self-consistency cycle. This includes all free-atom quantities (densities and potentials for the initialization), radial basis functions for all species, one-dimensional logarithmic and three-dimensional radial and angular integration grids, and fixed coefficients for the analytic long-range part of the Hartree potential.

Wrapper subroutine initialize_scf.f90 performs the initial s.c.f. cycle of the electronic structure calculation. In this step, the full three-dimensional integration grid is filled with fixed initial quantities (superposition of free-atom densities, potentials, and partition functions), the overlap and Hamiltonian matrix integrals are performed for the first time, and the initial Hamiltonian and overlap matrices are used to determine the storage requirements in the event of sparse matrix storage. If a two-electron Coulomb operator is needed (hybrid functionals, Hartree-Fock, MP2, $GW$ etc.), the three-center overlap matrix elements $(ij|\mu )$ of Eq. (3.76) (see Sec. 3.25 for details) and the Coulomb matrix of the auxiliary “resolution of the identity” basis set [denoted $V_{\mu \nu }$ in Eq. (3.76)] are precomputed. The most important task in initialize_scf.f90 is the initial solution of the Kohn-Sham equations Eq. (3.38), providing a first solution of the wave function coefficients $c_{jl}$. These are the starting point of every iteration of the s.c.f. cycle in the following step, subroutine scf_solver.f90.

With all preliminary information available, the task of scf_solver.f90 is to produce a self-consistent electron density, wave function, and all associated observables for a given, fixed nuclear geometry. The order of the cycle until convergence is reached is:

1. 1.

   Calculation of the Kohn-Sham electron density associated with the current wave function, $c_{jl}$

2. 2.

   electron density mixing and preconditioning, to produce the input electron density for the next set of Kohn-Sham equations

3. 3.

   decomposition of the electron density into atom-centered multipole fragments $\delta {\tilde{n}}_{\text{at},lm}(r)$ [see Eq. (Eq;mp)], and construction of the multipole components of the Hartree potential, $\delta {\tilde{v}}_{\text{at},lm}(r)$

4. 4.

   construction of the full electrostatic potential $v_{\text{es}}(𝐫)$ on all points $𝐫$ of the three-dimension integration grid

5. 5.

   Integration of the updated Hamiltonian matrix elements, $h_{ij}$

6. 6.

   possible addition of two-electron exchange matrix elements to $h_{ij}$

7. 7.

   solution of the Kohn-Sham equations Eq. (3.38), to produce an updated wave function $c_{jl}$

8. 8.

   computation of updated total energy components, and check of all convergence criteria.

After convergence is reached, scf_solver.f90 also performs some inevitable post-processing steps, including “scaled ZORA” perturbative corrections for the appropriate relativistic treatment, and band structure and density of states data output.

With a converged self-consistent solution at hand, the code can now perform any number of similar calculations for updated geometries, e.g., for a geometry optimization, molecular dynamics, etc. If so, an updated geometry is first produced by subroutine predict_new_geometry.f90. Note that this simple subroutine should also serve as the starting point for any other calculations involving multiple geometries, such as the calculation of “serial” geometries along a given set of coordinates, etc. For the updated geometry, all geometry-related storage arrays in the calculation and the overlap matrix must be recomputed in subroutine reinitialize_scf.f90. Following this, subroutine scf_solver.f90 is invoked again, and a new self-consistent solution is obtained.

The final step of the code is to produce, by post-processing, any information that can be obtained from the converged self-consistent wave function or electron density, including electrostatic moments, charge analyses, etc. Beyond this, only necessary cleanup tasks follow, most notably the deallocation of all storage arrays and the MPI infrastructure, and the final time accounting information.