4.7 Restarting FHI-aims calculations

This section gives background information on restarting FHI-aims calculations from the electronic structure.

Note that this section covers two versions of the restart process:

•

The elsi_restart keyword, which is simpler and somewhat more flexible in case the objective is just to restart a given electronic structure calculation from the previous one
•

The restart keyword, which produces differently formatted output and which supports a few more sophisticated tasks, for example, restarting using a rotated version of a structure calculated earlier.

If you just want to restart an electronic structure calculation in a flexible way, we suggest use the elsi_restart keyword. elsi_restart produces the same file format, irrespective of whether the linear algebra is performed in serial or in parallel and irrespective of how many MPI tasks were used. See the definition of that keyword for a brief explanation and for a few more bells and whistles (e.g., some form of density matrix extrapolation along the way, if the geometry is slightly adjusted prior to the restart).

Both elsi_restart and restart are regression tested at the time of writing.

4.7.1 Restart procedure based on the “elsi_restart” keyword

This is the most flexible restart procedure in FHI-aims at the time of writing, if the objective is simply to restart a calculation. It writes a compressed version of the density matrix or, in periodic calculations, a compressed version of the density matrix at each $k-$ point. (This means that many separate restart files can emerge for dense $k$ -point grids.)

To read or write a restart file or set of restart files in every scf iteration, simply insert the following line into your control.in file:

  elsi_restart read_and_write 1

This command will prompt FHI-aims to read its density matrix input files if they exist in the directory in which the code is executed. If they do not exist, the code will proceed assuming that this is a calculation from scratch. However, in either case the code will write new restart files in every single s.c.f. iteration and it will overwrite earlier restart files.

To have restart information written in (say) only every tenth s.c.f. iteration, just modify as follows:

  elsi_restart read_and_write 10

To only write restart files (but not read), use

  elsi_restart write 10

(of course, the number “10,” which defines the number of s.c.f. iteration after which a new file will be written, can be customized).

To only write restart information at the end of the s.c.f cycle, you can substitute the integer number as such:

  elsi_restart write scf_converged

Conversely, to only read existing files but not write new ones, use this version:

  elsi_restart read

Finally, it is possible to write restart files for one density functional and then use them to restart the s.c.f. cycle with another density functional. For example, you might wish to preconverge the electronic structure of a given system using a computationally simple functional first (e.g., PBE) and then restart from that density matrix but using a more costly functional, e.g., a hybrid functional.

This process is supported but it is important to note that by default, FHI-aims works on a reduced (time-reversal symmetric) $k$ -space grid for scalar relativistic semilocal DFT. In contrast, hybrid DFT requires all $k$ -points to be present. To include all $k$ -points in the semilocal DFT calculation as well and write all necessary restart files, please add the keyword symmetry_reduced_k_grid .false. to control.in:

  symmetry_reduced_k_grid .false.

Again, for more information please go to the elsi_restart keyword definition.

4.7.2 Restart procedure based on the “restart” keyword

Calculations can be restarted from previous wavefunctions using the keyword restart. In general, restarting the same calculation on the same computer using the same number of MPI-tasks will always work.

Important — Rotated systems:

FHI-aims wavefunctions are not symmetry-aware. While the same restart file can be used for molecules translated by a vector, this does not work for rotated molecules. For a way to restart systems with rotated geometries or superpositions of molecular densities see Section 4.7.5.

FHI-aims has two main means of saving the wavefunction for a later restart, depending on the used keywords.

Kohn–Sham eigenvector ("wavefunction"):

This is the default case always used for calculations using KS_method serial. The saved file contains some information on the system (number of basis functions, states, spin, k-points), the Kohn–Sham eigenvector, eigenvalues and occupations.

Kohn–Sham density matrix ("density matrix"):

This variant is only used for calculations using KS_method parallel. The saved file contains the density matrix ( $\hat{n}_{ij}=\sum_{l}f_{l}c_{il}c_{jl}$ ), either in full or sparse storage.

The corresponding keywords determining which version of storage is used are KS_method and density_update_method. See Fig. 4.10 for possible options and outcomes.

Figure 4.10: Scheme of different restart types depending on used methods.

For technical reasons there is not one single way, but the above mentioned variants with again, variants.

The density matrix based restart always writes one file per spin-channel, no matter how many parallel MPI tasks are used. The wavefunction based restart generally writes a single restart file for all cluster calculations and one file per MPI-task for periodic calculations. Due to this periodic calculations using the wavefunction restart need to be restarted with the same number of MPI-tasks.

4.7.3 Mixing variants - the "force_single_restartfile" option

For a subset of possible calculations the above described differences are not valid and it is possible to restart different calculations using one type of restart file. This is the case for all possible cluster calculations (using either KS_method serial or KS_method parallel) and for periodic calculations with only 1 k-point, $\Gamma$ -only (using either KS_method serial or KS_method parallel.

This is implemented via the keyword force_single_restartfile. Using this functionality, it is for example possible to profit from the parallel ELPA solver for huge systems and still get a file containing the KS-wavefunction for post-processing.

4.7.4 Comments on the ’restart’ starting point and on self-consistency

As noted in the description of the restart keyword, it is important to note that the restart infrastructure corresponds to a restart from the last Kohn-Sham orbitals, not from the last density.

In practice, this means that the code will restart from the last unmixed Kohn-Sham density, not from the last mixed density. When restarting from a non-selfconsistent starting point, this can lead to unexpected jumps in the calculated non-selfconsistent total energy between the “old” and the “new” (restarted) FHI-aims run.

Only the self-consistent total energy is truly meaningful. This quantity (the self-consistent) total energy should be the same for the same stationary density, when approached from different starting densities. However, note additionally that some systems may exhibit several different self-consistent stationary densities even for the exact same atomic positions and for the exact same density functional. A simple example are antiferromagnetic vs. ferromagnetic spin states in some systems. In such cases, the true ground state in a DFT sense is the stationary density that yields the lowest energy. It can be found by way of a global search for different stationary densities, usually by varying the initial density guess.

4.7.5 Rotating the FHI-aims wavefunction

This feature is new and tested for the case of weakly interacting organic molecules as can be found in organic crystals. If you encounter any problems or strange behaviour, please let us know!

In many molecular systems a rotation does not change the electron density (or wavefunction) of a system. Therefore, the same wavefunction can be used again to restart such a calculaction.

This is especially true for weakly interacting systems such as van-der-Waals bound organic crystals. In these cases, the wavefunction of each single molecule in the crystals is not too different from the isolated wavefunction of a single molecule. This can be exploited to start a FHI-aims calculation of huge molecular systems from a superposition of molecular densities instead of the usual superposition of atomic densities.

This functionality is implemented as an external tool using Python and can be found in the FHI-aims repository itself (Python-package aimsutils) or at https://gitlab.lrz.de/theochem/aimsutils.

If you use this functionality, please cite C. Schober, K. Reuter, H. Oberhofer, J. Chem. Phys., 144:054103, 2016[272].

Important pre-requisites

This functionality only works for restart files containing the wavefunction of the system, not the density matrix. Therefore, please make sure you have wavefunction restart files or use the option force_single_restartfile. The FHI-aims output file of the calculation is also necessary.

Rotations can be defined via Euler angles in zyz convention or quaternions. To avoid any confusions with rotational conventions (especially with Euler angles), the input geometry will be rotated and saved using the same rotation matrix. Please be sure to check if your expected and the actual rotation are the same (or directly use the generated geometry.in).

The final rotated restart file is again in the wavefunction format.

- example_calc
|
+-- control.in
+-- geometry.in             # necessary
+-- aims.out                # necessary
+-- restart.wavefunction    # name can be anything, will be parsed from outfile

Examples

The two examples are shown using pre-defined scripts, but it is totally possibly to use the Python package directly. The full Python API is available via the Sphinx-documentation of the aimsutils package.

Rotation of a single system:

The script RotateSingle.py can be used to create the rotated geometry and restart for a single calculation.

        $: # RotateSingle alpha beta gamma -t x y z
                          ’-----.-----’    ’---.---’
                          Euler angles   translation vector
        $: RotateSingle 45 0 123 5.4 2.3 6.3

        $: # RotateSingle x_i x_j x_k x_l -t x y z
                          ’------.------’ ’---.---’
                            quaternion   translation vector
        $: RotateSingle 0.33 -0.12 0.44 0.84 5.4 2.3 6.3

This will create a subfolder rotated with the new restart.combined and geometry.in.

Rotation and combination of wavefunctions

The script RotateMany.py reads a instruction file (rotations.in) with the following format:

            # Comments or empty lines are ignored

            # title will also be foldername
            title project_x
            # lattice vector for new cell if any
            lattice_vector 11.2 0.00 0.00
            lattice_vector 0.00 10.5 0.00
            lattice_vector 0.00 0.00 10.7
            # now all rotated molecules
            # Euler angles (zyz)
            # alpha beta gamma x y z parent_molecule
            33.4 12.4 167.0   10.  0.  0.  calc_2
             0.  45.    0.     3.5 2.1 5.4 calc_1
           130.  90.   23.     0.  2.6 12. calc_3
            # OR Quaternion
            # x_i x_j x_k  x_l   x y z
            -0.4 0.15 0.33 -0.84 4 0 0 calc_2
            -0.2 0.45 0.13  0.84 0 2 0 calc_1
             0.0 0.00 0.00 -0.32 4 4 5 calc_3

The necessary converged FHI-aims calculations must be in subfolders with the appropriate name:

- calc_something/ (main calculation folder)
|
+-- rotations.in
+-- calc1/
    |
    +-- ...
+-- calc2/
    |
    +-- ...
+-- calc3/
    |
    +-- ...

This will create a subfolder project_x with the rotated wavefunction assembled from the individual calculations defined in the script (restart.combined) and geometry (geometry.in).

Theory

The following should give a summary of how the rotation of wavefunctions can be done with FHI-aims. For details on the used methods and mathematics, please have a look at the references.

In FHIaims, a basis function is defined as

\displaystyle\Phi_{i,lm}=\frac{u_{i}(r)}{r}\cdot S_{l,m}(\theta,\phi),

(4.8)

with $S_{l,m}(\theta,\phi)$ being real valued spherical harmonics. These are obtained from the complex spherical harmonics $Y_{lm}$ via

\displaystyle S_{l,m}(\theta,\phi)=\begin{cases}\frac{(-1)^{m}}{\sqrt{(}2)}(Y_% {lm}+Y_{lm}^{*})&m>0\\ Y_{l0}&m=0\\ \frac{(-1)^{m}}{i\sqrt{(}2)}(Y_{l|m|}-Y_{l|m|}^{*})&m<0\\ \end{cases}

(4.9)

A wavefunction is then given by

\displaystyle\Psi_{k}(r)=\sum_{i=1}^{\text{n\_basis}}c_{i}^{k}\Phi_{i}(r)

(4.10)

with the coefficients $c_{i}^{k}$ .

A rotation of the molecule (rotation matrix $\mathbf{R}$ leads to the same set of YLMs $\mathcal{Y}$ (as they are fixed with respect to the xyz-coordinate system), but with different coefficients $\mathbf{c}$ (new linear combination of basis functions):

\displaystyle\mathbf{c^{\prime}}=\mathbf{R}\mathbf{c}

(4.11)

The rotation of complex YLMs can be done with Wigner D matrices ( $D^{l}_{mm^{\prime}}$ ). For real YLMs different schemes are available[196, 16, 34]. Due to the non-default YLM convention in FHIaims we construct our rotation matrices for different $l$ from the complex Wigner D matrices via the transformation matrix $\mathbf{C}^{l}$ (reference: [34]):

\displaystyle S_{l,m}=\mathbf{C^{l}}Y_{l,m}

(4.12)

The matrix $\mathbf{C}$ is constructed according to [34] with the constraint of the different sign convention in FHIaims. The real rotation matrix $\Delta^{l}(R)$ ,

\displaystyle\Delta^{l}(R)=\left(\mathbf{C}^{l}\right)^{*}\mathbf{D}^{l}(R)% \left(\mathbf{C}^{l}\right)^{t},

(4.13)

is then used to obtain the rotated coefficients for each $l, m$ for each basis function in the system.

\displaystyle\mathbf{c}^{\prime l}=\Delta^{l}(R)\mathbf{c}^{l}

(4.14)