3.10 SCF Cycle: Initialization, density mixing, preconditioning, convergence

There is a simple tool that can be used to visualize the s.c.f. convergence behavior of FHI-aims graphically for a given run. Preferably do this analysis on your desktop computer (i.e., copy over the necessary files). This is what you need to do:

• •

  Install the Grace 2D visualization program (http://plasma-gate.weizmann.ac.il/Grace/) on your computer.

• •

  Go to the directory with the FHI-aims output file you wish to analyze.

• •

  From the FHI-aims utilities directory, copy over the file scf_convergence_template.agr .

• •

  At the commandline, call the FHI-aims utility
  plot_scf_convergence.pl [FHI-aims output file]
  to visualize the s.c.f. convergence behavior of the FHI-aims output file. plot_scf_convergence.pl can be found in the FHI-aims utilities directory and must (of course) be called with the correct directory path preceding the file name.

If successful, this procedure will assemble and open a graphical representation of the s.c.f convergence of FHI-aims.

FHI-aims includes an on-the-fly convergence analysis that monitors the behavior of the s.c.f cycle and can make decisions to stop the s.c.f cycle early if convergence is unlikely to be achieved within a reasonable number of iterations. This feature is intended to avoid wasting computational resources on calculations that are converging too slowly or not at all, since the default in FHI-aims is to run until convergence is achieved or 1000 s.c.f. iterations are completed, see sc_iter_limit.

The convergence analysis tracks the density changes over multiple s.c.f iterations, fits linear trends to the natural logarithm of these changes, and predicts how many additional iterations would be needed to reach the target accuracy determined by the keyword sc_accuracy_rho. Based on these predictions, the code can decide whether to continue or stop the calculation.

The convergence analysis begins after a minimum number of s.c.f iterations (controlled by scf_analysis_start). At this point and in each subsequent iteration, it:

1. 1.

   Analyzes the trend of density changes over the most recent iterations

2. 2.

   Calculates the convergence rate (slope) of the natural logarithm of the density changes for both charge and spin densities over the last 15 s.c.f. iterations

3. 3.

   Predicts how many additional s.c.f iterations are needed to reach convergence based on the slope

4. 4.

   Checks for stagnation (when the density changes stop decreasing)

The predicted steps to convergence are evaluated as:

$$N_{\text{predicted}}=\lceil \frac{\ln (\mathrm{\Delta }{n}_{\text{target}})-\ln (\mathrm{\Delta }{n}_{\text{current}})}{\text{slope}}\rceil$$

, where $\mathrm{\Delta }{n}_{\text{target}}$ is the target density change defined by sc_accuracy_rho,

If you run a calculation and reach scf_analysis_start s.c.f. iterations, the code outputs the estimated number of steps to convergence in this manner:

This allows you to monitor the convergence behavior of a running calculation and may give you an early indication if something is going wrong or simply an estimate of how long the calculation will take to converge.

At a predefined decision point (controlled by scf_analysis_decision_point), the code evaluates whether convergence is proceeding acceptably. If the predicted number of steps to convergence exceeds a threshold (controlled by scf_max_allowed_steps_at_decision), the calculation may be stopped.

In case the calculation is allowed to continue, multiple evaluations at regular intervals (every scf_max_allowed_steps_at_decision s.c.f iterations after the first decision point) are performed. The number of these evaluations is controlled by scf_max_decisions_allowed. If convergence remains too slow after all allowed decision points, the calculation is stopped automatically.

The automated convergence analysis will not stop the s.c.f cycle in the following cases:

• •

  Single free atom calculations: These calculations are considered trivial in cost.

• •

  Very fast s.c.f iterations: When individual s.c.f iterations are extremely cheap, namely being less than 1 second per s.c.f. iteration, the analysis will not stop the s.c.f cycle.

• •

  When sc_iter_limit is explicitly set: If you manually specify a maximum number of s.c.f iterations using sc_iter_limit, the automated convergence analysis is bypassed entirely, and the calculation will run up to your specified limit.

The initial_charge keyword always applies to the last atom previously specified in input file geometry.in. The charge is introduced by using an ionic instead of neutral spherical free-atom density on that site in the initial superposition-of-free-atoms density. Note that initial charge densities are generated by the functional specified with xc for DFT-LDA/GGA, but refer to pw-lda densities for all other functionals (hybrid functionals, Hartree-Fock, …). Note also that the sum of initial_charge values should be equal to the value of the charge specified in control.in (or equal to 0 if charge is not specified). If this is not the case FHI-aims will stop, unless the keyword override_initial_charge_check is used.

The initial_moment keyword always applies to the last atom previously specified in input file geometry.in. The moment is introduced by using a spin-polarized instead of an unpolarized spherical free-atom density on that site in the initial superposition-of-free-atoms density. Note that initial charge densities are generated by the functional specified with xc for DFT-LDA/GGA, but refer to pw-lda densities for all other functionals (hybrid functionals, Hartree-Fock, …).

This keyword decides whether key s.c.f. convergence parameters will be adjusted automatically during the s.c.f. cycle, based on a simple estimate of the system character according to its approximate HOMO-LUMO gap or band gap.

The number keyword determines in which iteration of the s.c.f. cycle the adjustment will be attempted. A zero value corresponds to the initial s.c.f. cycle; values of 1, 2, etc. correspond to an update after the first, second, etc. s.c.f. iterations are almost complete and their eigenvalue spectra known.

The frequency keyword determines for which full s.c.f. cycle (i.e., for which geometry step) an adjustment will be made:

• •

  once means that an adjustment of the s.c.f. parameters will only be made in the first geometry in a geometry relaxation or MD run.

• •

  always indicates that an adjustment will be made for every new s.c.f. cycle, i.e., for every new geometry in a run.

• •

  never indicates that no adjustment will be attempted.

Parameters will only be adjusted if they are not explicitly set by a keyword in control.in. Any parameter that is included in control.in will not be modified by the adjust_scf keyword. For frequency once or always, the following parameters may be adjusted:

• •

  The initial default value of charge_mix_param for mixer pulay will be set to 0.05 (i.e., an overall cautious value). For frequency never, the default value of the charge_mix_param keyword remains at its usual default value 0.2 (for many metallic or spin-polarized systems, this is a fairly aggressive value).

• •

  In s.c.f. iteration number, the current estimated value of the HOMO-LUMO gap (for solids, the band gap) is checked. If the system shows fractional occupation numbers or if the estimated gap has a value of less than 0.2 eV at this point, the system is likely near-degenerate or metallic and the s.c.f. cycle could be diffcult to converge. In this case, the charge_mix_param is set to 0.02 – a cautious value but, in conjunction with the Pulay mixer, still surprisingly effective. The broadening of the occupation numbers near the Fermi level is increased to occupation_type option 0.05 [eV]. These default values of 0.02 for the charge_mix_param and the width of 0.05 eV for occupation_type can be modified via the keyword adjust_scf_param.

• •

  If, instead, the gap is found to be equal or greater than 0.2 eV in s.c.f. iteration number, the rather aggressive default charge_mix_param 0.2 is kept for the Pulay mixer, and the default broadening value (suitable for non-metallic systems) occupation_type option 0.01 [eV] is also kept. These default values of 0.2 for the charge_mix_param and the width of 0.01 eV for occupation_type can be modified via the keyword adjust_scf_param.

Example:

adjust_scf_param lowgap charge_mix_param 0.05
adjust_scf_param lowgap occupation_width 0.1
adjust_scf_param largegap charge_mix_param 0.5

In this example:

• •

  In the case of a low gap (i.e. less than 0.2 eV), the charge_mix_param will be adjusted to 0.05, after s.c.f iteration number as specified with adjust_scf. At the same time, the width of occupation_type will be adjusted to 0.1 eV.

• •

  In the case of a large gap (i.e. greater than 0.2 eV), the charge_mix_param will be adjusted to 0.5, after s.c.f iteration number as specified with adjust_scf.

If adjust_scf is not used, adjust_scf_param will do nothing.

The keywords xc_pre and elsi_restart both allow one to specify a non-standard (and usually faster) initialization of the electronic structure that is later calculated self-consistently for the exchange-correlation method specified by the usual xc keyword.

One frequent problem is that they can be accidentally set togethe. In that case, the initialization using xc_pre can undo the often highly valuable initialization already present from the elsi_restart information. For example, xc_pre can replace the highly expensive information from a previous hybrid DFT run and replace it with the much less suited information from, say, a PBE pre-initialization.

This is, ostensibly, a user error, but can cost immense amounts of CPU time especially for the highest-value calculations. Therefore, we safeguard against this possibility.

When the xc_pre keyword is set, elsi_restart can still be used at the same time – but only if no elsi_restart restart files are actually available to read (i.e., if all that elsi_restart would do is write). If elsi_restart actually does have information to initialize the calculation, keyword xc_pre will only be allowed to be used if allow_restart_xc_pre is explicitly set to be true. If it isn’t, then the code will stop and explain to the user that the two should not be used together.

The allow_restart_xc_pre keyword has no effect when used together with the older restart keyword.

See Ref. [36] for details regarding the available density mixing algorithms. In the simplest case of a linear mixer, value specifies a constant value ${\widehat{G}}^1$ to mix the output density of the Kohn-Sham Equations in iteration number $\mu$, $n_{\text{KS}}^{(\mu )}$, with the (already mixed) input density that defined those equations, $n^{(\mu -1)}$:

$$n_{\text{dmp}}^{(\mu )}=n^{(\mu -1)}+{\widehat{G}}^1(n_{\text{KS}}^{(\mu )}-n^{(\mu -1)}).$$

(3.39)

If a preconditioner is specified, charge_mix_param defines an additional linear factor to that preconditioner. In case of a pulay mixer, all density residuals are mixed with this factor.

In general, the best choice for value is system-dependent, and also depends on the chosen mixer algorithm. In general, please also see Appendix A.8, since several keywords can be used to alleviate s.c.f. mixing instabilities.

• •

  In principle, a linear mixer will always converge with a sufficiently small value. In easy cases, we recommend value=0.1-0.2, but in difficult cases, “sufficiently small” can mean one to three orders of magnitude(!) lower, i.e., the process can be apallingly slow.

• •

  For a straight pulay mixer, our default value is adjusted according to the adjust_scf keyword, depending on the estimated band gap / HOMO-LUMO gap. For non-metallic systems, we choose a conservative value of 0.2. In metallic systems or systems that are otherwise problematic, the default value set by adjust_scf is value=0.02. Note that this small value does not necessarily correspond to slow mixing since the Pulay mixer will learn over time and accelerate the mixing process. These defaults can be modified via the adjust_scf_param keyword.

• •

  In metallic systems, density oscillations can occur from one iteration to the next (charge sloshing). This can be alleviated by a preconditioner. With a preconditioner and pulay mixer specified together, value is still important and may be chosen around 0.05 .

See also the mixer and preconditioner keywords.

relative_fp_charge_mix determines the under-relaxation of the fixed point part of the Broyden mixer s.c.f. cycle together with charge_mix_param. relative_fp_charge_mix and charge_mix_param are multiplicative, and if no history is included the effective under-relaxation is relative_fp_charge_mix times charge_mix_param.

Sets the default initial spin moment per atom for all atoms whose initial_moments are not specified explicitly in geometry.in.

However, do not use the default_initial_moment keyword unless you know exactly what you are doing. default_initial_moment is a very dangerous keyword if used unwisely. It sets a per-atom moment. Remember that, if a structure contains (say) ten atoms, the initial moment per structure will be ten times the default_initial_moment per atom - possibly, assigning entirely unphysical initial spin moments to atoms that should never have seen a local spin polarization in the first place. The result will be a terrible s.c.f initialization, leading to unphysical results or to no s.c.f. convergence at all.

Instead, think about which exact atoms should carry an initial spin moment in your structure, and what that moment should be. Then, use initial_moment keywords for each atom in your geometry.in file. Doing so will usually allow you to specify a much better initial spin density (and, thus, better s.c.f. convergence) than any blanket default_initial_moment keyword can ever do.

We keep the default_initial_moment keyword for convenience when dealing with some magnetic materials. However, with a single, unfortunate choice of default_initial_moment, one can initialize a complex structure with a completely unphysical spin density. This will lead to severe issues with the s.c.f. cycle.

If there is at least a single initial_moment keyword specified in geometry.in and if the keyword default_initial_moment is not set explicitly, then the default_initial_moment will be zero for all atoms for which no initial_moment is specified explicitly.

Note that at least one moment in the system must be set to a non-zero value in order to reach any spin-polarized state at all. If the entire initial spin polarization is zero, the final s.c.f. result will also not be spin-polarized, no matter how magnetic the system is in reality.

Only if no initial_moment is included in geometry.in at all in a spin-polarized calculation, the default_initial_moment must be specified explicitly by the user for the code to run at all.

For isolated free-atom calculations, default_initial_moment hund can be used. This will result in the usual high-spin atom initialization characteristic of free atoms. However, it is much better (and often the only way) to use the force_occupation_basis keyword to specify the expected orbital occupation of a free atom explicitly. Many free atom calculations will not converge in a s.c.f. cycle if the expected symmetry breaking (and, thus, orbital occupation) is not specified.

Warning: It is not advisable to set a default_initial_moment other than zero in structures in which only a few atoms actually carry spin (such as a large molecule with a few transition metal atoms). Rather, a good choice would be to set a zero default_initial_moment and then set finite initial_moment values for individual atoms in geometry.in.

And again: Setting default_initial_moment to zero and not specifying any initial_moment values in geometry.in will lead to a zero spin state in the final result, simply because the system is never given any indication which way to break its symmetry. So one does need to set a finite moment somewhere if a finite-spin converged solution is expected in the calculation.

The upshot is: It pays to think about the right spin initialization. Do not use the default_initial_moment keyword unless you know exactly what you are doing. Instead, use individual initial_moment keywords in geometry.in. There may be multiple different self-consistent solutions, and in spin-polarized systems, this can happen for very natural reasons (e.g., ferromagnetic vs. antiferromagnetic states). Similar to geometry optimization, starting from a very bad initial guess can cause terrible problems, particularly for s.c.f. convergence; but physically unreasonable solutions may also result. Conversely, a good starting point may greatly simplify a calculation. It really does pay to think about a detailed, physically reasonable spin initialization.

The force_complex_eigenvectors keyword can be useful e.g. when combined with elsi_restart, to ensure you have the same type of density matrix written as that needed by a restarted calculation. It can be the case that the density matrix is different, e.g. if you perform an SCF calculation with less than 3 k-points in each direction, then the calculation uses real eigenvectors and hence a real density matrix would be saved with elsi_restart write scf_converged. However, if restarting with elsi_restart read and including either output band, output band_mulliken, output polarization, output Z2_invariant, or use_spin_texture_scf, then complex eigenvectors and a complex density matrix are required in this calculation. elsi_restart read can handle for this, reading in the real density matrix from file, and setting the complex density matrix array in the calculation based on this with complex elements to be 0.0. However this could require a few more SCF steps than if you could restart with a complex density matrix from file - which could be achieved by combining force_complex_eigenvectors .true. with elsi_restart.

The opposite case could also happen - where a complex density matrix is saved and then a restarted calculation would ordinarily use a real density matrix. In this case, elsi_restart read would stop with an error, to avoid data loss from not being able to read the complex part of the density matrix from file. force_complex_eigenvectors .true. could also be used here, in the restarted calculation, to ensure consistent density matrix types.

This option is not required under normal circumstances. It is mainly useful to produce / test a fast, non-self-consistent solution for a superposition-of-atoms potential that yields only the sum of eigenvalues as a result. If a non-self-consistent total energy is needed (correct only for the non-spinpolarized free-atom density!), running a normal calculation with sc_iter_limit=0 is the better way.

Options for type are:

• •

  sc: Self-consistent Kohn-Sham potential in each s.c.f. iteration

• •

  superpos_pot: Superposition of free-atom potentials, evaluation only once (no self-consistency cycle). Restriction: This method works only for self-adapting angular_grids (i.e., angular_acc not equal zero for at least one species). This also means that the option will not perform well in periodic boundary conditions. A fix is simple, contact us if needed.

• •

  superpos_rho: Superposition potential created from sum of free-atom densities; evaluation only once (no self-consistency cycle)

• •

  non-self-consistent: Same as superpos_rho.

A description of the different versions of the algorithm used to set the occupation numbers of the free atom minimal basis in the initial density superposition for a charged or spin polarized calculation is as follows:

init_atom_occs_version 1 is the original implementation, and proceeds as:

1. 1.

   The total number of electrons for the atom/ion to be initialized is passed in, which is equal to the nuclear charge (specified by the keyword nucleus) minus the initial_charge.

2. 2.

   If an initial_moment is specified, the number of electrons in each spin channel is corrected to achieve this value.

3. 3.

   Hund’s rule/the Aufbau principle is applied for the electrons in each spin channel.

4. 4.

   A selection of non-Aufbau occupations is corrected for (species with 29, 47, 58, 64 and 79 electrons, i.e. Cu, Ag, Ce, Gd and Au -like species).

5. 5.

   The occupation numbers for the minimal basis functions are output and set.

From a chemically intuitive perspective, this algorithm has a number of deficiencies:

• •

  The algorithm has no knowledge of whether a species is a neutral atom, an anion or a cation. This is consequential because, e.g. for transition metal/lanthanide/actinide cations, electrons are removed in a different order to which they are added. As an example, for first row transition metal cations, 4s electrons should be removed before 3d. This algorithm does not do this, and would predict both neutral V and Fe³⁺ to have the same electron configuration - 4s² 3d³, whereas Fe³⁺ should be 4s⁰ 3d⁵.

• •

  The algorithm does not take into account all non-Aufbau occupations, in particular, it is missing Cr, Nb, Mo, Ru, Rh, Pd, La, Pt, Ac, Th, Pa, U, Np, and Cm.

init_atom_occs_version 2 was introduced in version stamp 260316 to try to correct for these deficiencies, with a more sophisticated algorithm, which is as follows:

1. 1.

   Both the nuclear charge and the initial_charge are passed as arguments, rather than simply the total number of electrons, to be able to determine whether a species is a neutral atom, cation, or anion and fill appropriately.

2. 2.

   Based on the elemental species (specified by either the closest integer number of electrons to nucleus, an unambiguous 1 or 2 character definition of species, or specification of the element keyword otherwise) the neutral elemental species is filled with electrons according to Hund’s rule/the Aufbau principle.

3. 3.

   A correction is made for a complete list of non-Aufbau occupations.

4. 4.

   A determination is made based on the nuclear charge and intitial_charge as to whether the species is an anion or cation.

   • •

     If it is cation, then electrons are removed from the highest energy orbital (e.g. 4s before 3d) and paired before unpaired, until the desired charge is achieved.

   • •

     If it is an anion, then electrons are added starting from the lowest energy unoccupied/partially occupied orbital, until the desired charge is achieved.

5. 5.

   If an initial_moment has been requested, then a check is made as to whether the current moment is already equal to the initial_moment, or whether the current moment needs to be modified.

   • •

     If the magnitude of the current moment is greater than or equal to the initial_moment, then electrons in partially occupied sub-shells, starting from the highest energy one, are flipped until the correct moment is achieved.

   • •

     If the magnitude of the current moment is less than the desired initial_moment then electrons are removed from one spin channel and added to another until the desired initial_moment is achieved. The rules for adding and removing electrons in each spin channel follow those above for the anion and cation, respectively. There is a preference to only touch partially occupied valence electrons/orbitals if the desired spin moment is possible by only using these.

6. 6.

   The occupation numbers for the minimal basis functions are output and set.

After some initial testing, init_atom_occs_version 2, has been shown to lead to faster convergence for some systems, in particular spin polarized systems with transition metals. As it is a new feature, more testing will be done before it is potentially set as the default option.

Try only if the standard / preconditioned pulay mixer definitely fails. Keywords ini_linear_mix_param, ini_spin_mix_param can be used to specify separate mixing parameters for the initial linear mixing.

ini_linear_mixing should only be tried if the standard algorithms provably fail. In that case, value should be relatively small.

ini_spin_mix_param should only be tried if the standard algorithms provably fail. In that case, value should be relatively small.

FHI-aims provides three mainstream density mixing algorithms across the s.c.f. cycle, type linear, pulay, and broyden. We here only give a brief summary of options, please see Ref. [36] for further references and for the exact mathematical details.

For most practical purposes (non-pathological systems), Pulay’s DIIS mixing algorithm [254] is robust and fast, and should be the algorithm of choice. For this algorithm, n_max_pulay $n$ determines the number of past iterations $\mu -k$ ($k$=1,…,$n$) to be mixed with the Kohn-Sham output density of iteration $\mu$. charge_mix_param determines an additional (system-dependent) linear factor that is multiplied with the output density change of the Pulay mixer. Normally, this (and perhaps a preconditioner) is all you need to do to ensure convergence.

In some pathological cases, reaching self-consistency is a more tricky problem. Broadly speaking, these are systems with a small HOMO-LUMO gap (band gap) and/or several competing possibilities for a self-consistent solution. Specifically, these difficult cases include:

• •

  Large metallic systems (e.g., slabs), where charge may “slosh” from one end of the system to another before reaching self-consistency. In that case, the pulay mixer may be used together with a large charge_mix_param and a preconditioner (see that keyword) to dampen the resulting oscillations. Also make sure that occupation_type is set to a sufficiently large broadening of occupation numbers near the Fermi level in metallic systems.

• •

  Spin-polarized systems with competing spin states. A classic. If problems arise, playing with ini_linear_mixing, the charge_mix_param and spin_mix_param and further options listed in this section may help. Likewise, setting a specific fixed_spin_moment may be helpful. Finally, different initial_moment settings may easily switch between different metastable self-consistent spin states (e.g., ferromagnetic vs. antiferromagnetic), and should be tested separately if different competing spin states are suspected.

• •

  Systems near a level crossing (even dimers, if two or more Kohn-Sham levels of different symmetry come close for a given binding distance). Apart from the usual mixing mechanisms, keyword occupation_type with a larger broadening near the Fermi level may help alleviate this situation.

• •

  Spin-polarized free atoms. The simplest conceivable systems may exhibit unexpected problems towards self-consistency, likely because the electron density can rotate between several fully equivalent spin states. Here, demanding a specific orbital occupation using the force_occupation_basis keyword may be useful.

In principle, even in the toughest cases a linear mixer will always converge with a sufficiently small charge_mix_param. Unfortunately, “sufficiently small” can mean a charge_mix_param of 10^-2-10^-4, i.e., the process can be apallingly slow. Playing with the pulay mixer settings is usually the better strategy, unless a proof of principle is sought.

The broyden mixer is an improvement on the linear mixer. The broyden mixer works by effectively dividing the s.c.f. cycle into two separte parts: the space in which we have local information gained by previous evaluations of the s.c.f. cycle, and the remaining space in which we do not have information. charge_mix_param determines the linear factor which under-relaxes the density predicted by the broyden mixer, n_max_broyden n controls the number of past iterations used to construct the next estimate, and relative_fp_charge_mix is the additional (multiplicative) under-relaxation affecting only the step length in the space where we have no further information.

Note that a modification is needed when going beyond DFT-LDA/GGA (Hartree-Fock, hybrid functionals, …). In that case, the density implicitly enters the two-electron exchange operator (via the density matrix, ${\widehat{n}}_{ij}={\sum }_l{f}_l{c}_{il}{c}_{jl}$, where $i$ and $j$ label basis functions, and $l$ label the Kohn-Sham states), and should also be mixed.

By default, for linear mixing, we do not mix the exchange operator, unless keyword use_density_matrix_hf is enabled. The latter is the default if the pulay mixer is selected. Then, the density matrix is submitted to the same Pulay mixing factors as the normal charge density $n(r)$ before constructing the exchange operator. Note that we do not have a formal density matrix available that corresponds to the initial superposition of free-atom densities, making this form of mixing slightly less efficient than for normal Kohn-Sham DFT-LDA/GGA.

This option is perhaps useful when there are definite convergence problems with the standard mixing algorithms, but can otherwise safely be ignored.

A larger number of stored iterations can sometimes lead to a stabilization of the mixing process. Choosing number too large (e.g., 20 and above), though, may destabilize the Pulay matrix, which can become near-singular.

Note that the storage effort associated with Pulay mixing is significant on systems with few CPUs / low memory. Each additional stored iteration requires the storage of two charge density residuals, and two times three charge density gradient residual components. For large systems, low memory, and overall stable mixing, reducing n_max_pulay may be a way to get a given calculation below the most difficult memory barriers.

A larger number of stored iterations can sometimes lead to a stabilization of the mixing process. Choosing number too large (e.g., 20 and above), though, may destabilize the Broyden matrix, which can become near-singular.

Note that the storage effort associated with Broyden mixing is significant on systems with few CPUs / low memory. Each additional stored iteration requires the storage of two charge density residuals, and two times three charge density gradient residual components. For large systems, low memory, and overall stable mixing, reducing n_max_broyden may be a way to get a given calculation below the most difficult memory barriers.

It’s our tentative observation that a larger mixing parameter (0.5-0.8) is sometimes helpful with the preconditioner, but after the preconditioner is switched off, a smaller mixing parameter (as set by charge_mix_param) may be desirable. prec_mix_param can provide the needed separate setting, if desirable.

Default values:

• •

  Non-periodic systems: preconditioner kerker off (no preconditioner).

• •

  Periodic systems: preconditioner kerker 2.0 (Preconditioner with a momentum of $q_0$=2.0 bohr^-1).

See Ref. [36] regarding the mathematical definition of the Kerker-type [208, 229, 164] preconditioner, which is the currently implemented form.

See keyword precondition_max_l for the angular momentum cutoff specified for the kerker preconditioner.

keyword can have the following forms, controlling various aspects of the preconditioner:

• •

  kerker :

  • –

    if followed by off : No preconditioner used.

  • –

    if followed by value : value is a real positive number, indicating the characteristic momentum $q_0$ associated with the Thomas-Fermi type screening assumed in the preconditioner (in bohr^-1). Typical values in the literature range around 1.0-2.0 bohr^-1, but larger values may be useful in small clusters.

• •

  turnoff : To avoid any residual influence on the s.c.f. cycle, the preconditioner can be switched off at a given level of s.c.f. convergence, leaving the remaining convergence to a pure pulay mixer. Possible types of turnoff criteria are:

  • –

    charge : value refers to the root-mean-square deviation between $n^{(\mu -1)}$ (the mixed and preconditioned input density to the Kohn-Sham Equations) and $n_{\text{KS}}^{(\mu )}$ (the unmixed output density from the Kohn-Sham Equations). Default: sc_accuracy_rho.

  • –

    energy : value refers to the total energy difference between two successive iterations [in eV].

  • –

    sum_ev : value refers to the difference in the eigenvalue sums between two successive iterations [in eV].

  All requested convergence criteria for the preconditioner must be fulfilled. If no explicit turnoff criterion is set, the preconditioner turnoff charge, energy and sum_ev defaults to the same values as sc_accuracy_rho, sc_accuracy_etot, and sc_accuracy_eev, respectively.

• •

  dielectric : The inverse of the microscopic dielectric function, ${ϵ}^{-1}(𝐫,{𝐫}^{\prime };0)$ is used to precondition the density. This option is experimental and computationally expensive but the the dielectric function can be theoretically justified to be a good preconditioner. Works only for non-periodic systems.

• •

  none or off : No preconditioner used.

We use a partitioned atom-centered multipole rewrite of the Kerker preconditioner in angular momentum space, similar in spirit to the Hartree potential (see Ref. [36] for details). In principle, precondition_max_l is thus the equivalent of the l_hartree angular momentum cutoff for the expansion. However, since we here precondition a density difference (which reduces to zero as we approach self-consistency), and since we are combatting charge sloshing across potentially faraway parts of the systems, preconditioning the atom-centered monopole component of the density difference (value=0) is often all that is needed for the preconditioner to work. This is also the numerically most efficient way of running the preconditioner.

This keyword is experimental and derived fully empirically. For notorious cases the keyword may help to achieve convergence (faster). For $q_0=1$ the method is identical to the standard kerker preconditioner.

Note: A more generic restart functionality is provided by the elsi_restart – please see that keyword for a description. The elsi_restart keyword is the supported restart functionality going forward unless you have special reasons to prefer the simple restart keyword.

The restart keyword summarizes a set of different cases for which the electronic structure information in FHI-aims can be saved and then later used to restart a calculation from that point. restart comes in many internal variants and, for example, restarting a calculation with more MPI tasks than a previous run is not possible. The alternative elsi_restart keyword (which uses a different storage format than the restart keyword) is much more flexible and may be the better choice.

If file is not yet present, the calculation simply writes that file during the run. If file is already present, it is read and the wave function contained therein is used to restart the calculation, instead of a fresh superposition of free atoms initialization. Mind that restart is currently not supported for keywords load_balancing and use_local_index so that file will not be written or read when either keyword is set.

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-self-consistent starting point, this can lead to unexpected jumps in the calculated non-self-consistent total energy. Only the self-consistent total energy is truly meaningful and this (the self-consistent) total energy should be the same for the same stationary density. (Note also that some systems may exhibit several different self-consistent stationary densities – 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 one with the lowest energy and must be found by varying the initial density guess.)

In parallel runs, there is one file for each process, numbered as fileXXX. See also restart_read_only and restart_write_only. There are limited checks on whether or not the restart file provided is actually from the same system, but ensuring that a given restart file works is mainly the user’s responsibility. See also MD_restart for more information on the separate restarting process for molecular dynamics.

A much more complete overview of the restart infrastructure in FHI-aims can be found in the dedicated Section 4.7.

Similar to keyword restart, but does not overwrite file at any time (this may facilitate another restart from the same file later on).

Note: A more generic restart functionality is provided by the elsi_restart, described further below. We recommend to use elsi_restart unless you have special reasons to prefer the simple restart keyword.

Similar to keyword restart, but does not read the restart file in case it already exists.

Note: A more generic restart functionality is provided by the elsi_restart, described further below. We recommend to use elsi_restart unless you have special reasons to prefer the simple restart keyword.

See also restart.

Only relevant when using keyword restart. For technical and efficiency reasons FHI-aims uses different types of restartfiles depending on the eigenproblem solver (see Section 4.7 for details). In cases where the wavefunction is needed (e.g. for external post-processing, …), the default restart handling might make it difficult to do so.

This option only works for cluster (non-periodic) and periodic $\mathrm{\Gamma }$-only calculations.

Available options for task are:

• •

  write : Writes restart info (in ELSI format) to files. The info will be written to files in every freq s.c.f. iterations, and will always be written out after a converged cycle. The number of files depends on the choice of elsi_restart_use_overlap.

• •

  read : Reads restart info from files and uses it (instead of free atom superpositions) as the initial guess of an s.c.f. cycle. The code will search for files written by the write option. If relevant files do not exist, the code will abort. freq is not used for this option. The restarted run can use any number of MPI tasks, i.e., not necessarily the same number as used in the original run.

• •

  read_and_write : Performs what the write and read options do. Note that if the restart files do not exist, the code will still proceed normally.

Restarting hybrid functional calculations from the density matrix obtained using a different (i.e. less expensive) functional may reduce the number of SCF steps taken by the hybrid calculation.

Note that if a calculation is restarted using a different functional than that used to obtain the density matrix, the atomic_solver should be set to that used to obtain the density matrix.

Note also that the elsi_restart keyword – when it has actual past restart information available to read – cannot be used together with the xc_pre keyword, unless explicitly allowed using the allow_restart_xc_pre (see there for more information). We do this to prevent unintended duelling inializations. The two keywords can still be used together if all that elsi_restart would do is write information, not read it.

If you have a calculation with less than 3 k-points in each direction, and plan to restart a calculation with additional keywords, for example, output band, output band_mulliken, output polarization, output Z2_invariant, or use_spin_texture_scf then it can be useful to combine elsi_restart with the keyword elsi_restart with force_complex_eigenvectors .true.

The reason for this is that a calculation with less than 3 k-points in each direction will natively use real eigenvectors, and hence a real density matrix. However, the keywords output band, output band_mulliken, output polarization, output Z2_invariant, or use_spin_texture_scf requires the use of complex eigenvectors and hence a complex density matrix. Therefore, by either adding or removing these keywords to calculation with less than 3 k-points in each direction that you wish to restart, you could end up with a mismatch between the type of density matrix (real or complex) that is written, and the type that is required by the restarted calculation.

elsi_restart read can handle the case where it reads a real density matrix from file but a complex density matrix is required by the calculation, by simply setting the complex parts of the calculation’s density matrix to 0.0. However, elsi_restart read would throw an error in the opposite case, where a complex density matrix from file is attempted to be read by a calculation that needs a real density matrix, to avoid data loss from not being able to read the complex part of the density matrix from file. In this latter case, using force_complex_eigenvectors .true. in combination with elsi_restart read can ensure the calculation is possible. In the former case, force_complex_eigenvectors .true. could also be useful in the elsi_restart write stage, as sometimes reading a real density matrix into a calculation that needs a complex density matrix can result in a few extra SCF steps than if a complex density matrix was able to be read. See also the definition of force_complex_eigenvectors for further details.

Available options for task are:

• •

  write : Writes restart coefficients to a file ’ri_restart_coeffs.out’. The information will be only be written to files after a converged cycle. value defines a cutoff radius for each product basis function, such that basis functions centred on an atom further than this distance from the current integration point are neglected when calculating the overlap matrix. This radius is given by value multiplied by the charge radius of the product basis function in question. The default value should always be a safe choice.

• •

  read : Reads restart info from file ’ri_restart_coeffs.out’ and uses it (instead of free atom superpositions) to build an approximate density as the initial guess of an s.c.f. cycle. If relevant files do not exist, the code will abort. value sets sc_iter_limit; if set to 0, the calculation will be initialised with the RI density but no SCF steps will be completed; the density will be assumed to be converged.

• •

  read_and_write : Performs what the write and read options do. Note that if the restart files do not exist, the code will still proceed normally.

The input file "eigenstate_occs.in" must be constructed and specified by the user, and contain the custom occupation of each eigenstate in the set of (if applicable, k-point symmetrised) eigenstates. The keyword eigenstate_info can be used in the first SCF step of a two step SCF + RI procedure (see ri_skip_scf) to output information on the eigenstates. "eigenstate_occs.in" must contain one occupation value per line, and contain the same number of lines as the file "eigenstate_info.out" produced by the keyword eigenstate_info.

RI decomposition is then performed on the electronic density constructed from these custom occupations, using the standard RI keywords ri_density_restart write and ri_full_output.

This keyword allows to catch obviously ludicrous runs. If it triggers, something went seriously wrong during the mixing procedure and the settings for mixers, occupation broadening, preconditioner etc. should all be revisited very carefully.

The alternative setting sc_abandon_etot never switches the abort off, restoring the previous behaviour.

Very sensitive criterion for s.c.f. convergence. Usually, value=10^-3 eV is enough to indicate a reliable total-energy and force convergence. If value is set to zero or not given, the sum of eigenvalues will not be used as a convergence criterion.

The Harris-Foulkes form of the functional is used as the total energy in FHI-aims (see Ref. [36] for a brief discussion). A typical tight convergence criterion is value=10^-6 eV. If value is set to zero or not given, the total energy will not be used as a convergence criterion.

Attention: If keyword sc_accuracy_forces is set in control.in, forces are by default computed, regardless of whether or not they are later needed. The rationale is that the only way to check a requested force convergence criterion is to compute the necessary forces, despite the added numerical effort. For single-point calculations (no relaxation required), sc_accuracy_forces should therefore not be set unless explicitly needed for some reason.

One can explicitly set the keyword value to the string ’not_checked’. In this case, the forces are computed in only a single shot, their s.c.f. convergence will not be checked.

In this case, the code now relies on the default convergence criterion for the electron density itself (sc_accuracy_rho) to be sufficiently tight to guarantee good forces. This avoids excessive numbers of force evaluations in production runs.

Calculating forces is relatively expensive, so this convergence criterion is either not checked (default behavior), or it is checked after the purely electronic / energetic criteria, sc_accuracy_eev, sc_accuracy_etot, sc_accuracy_rho, are all fulfilled. To avoid too many iterations with force computations before the forces are converged, it is important to set the other criteria (especially sc_accuracy_eev, which checks a non-variational quantity) to tight convergence, as indicated above. For sc_accuracy_forces itself, e.g., value=10^-4 eV/Å is a reliable and robust criterion to avoid noise in geometry relaxations.

For simple structure relaxation, not_checked will often be good enough if the electronic criteria (especially sc_accuracy_rho) are set reasonably tightly.

For molecular dynamics, however, even slightly underconverged forces may lead to long-term energy drifts in (nominally) constant-energy runs. Here, it is a good idea to try out in explicit tests how tightly sc_accuracy_eev must be set to guarantee drift-free trajectories.

By default, FHI-aims checks separately the convergence of the charge density and the spin density, using the same criterion. Specifically, the unmixed output density of the Kohn-Sham Equations, $n_{\text{KS}}^{(\mu )}$, is checked against the input density $n^{(\mu -1)}$ to the same equations. A typical tight convergence criterion is value=10^-5.

sc_accuracy_rho is the most important s.c.f. convergence criterion. If the keyword is not set in control.in, the code chooses its initial default as follows:

1. 1.

   In general, but especially if the stress tensor is needed and/or if the output keyword is used, use fairly tight criteria:

   • •

     Up to 6 atoms in geometry.in: value=10^-6

   • •

     Between 6 and 6000 atoms in geometry.in: value=10${}_{}{}^{-6}\cdot \sqrt{\text{n\_atoms}/6}$

   • •

     Above 6000 atoms in geometry.in: value=10${}_{}{}^{-6}\sqrt{1000}$

2. 2.

   If forces are corrected for s.c.f. convergence incompleteness by the force_correction keyword or if forces are not used at all:

   • •

     Multiply the defaults of item 1 by a factor of eight, since the convergence of total energies or of corrected forces is already very good.

Warning. Setting sc_accuracy_stress to a finite non-negative value will result in a disproportionately large computational cost. In this case, the number of stress evaluations per relaxation step is at least doubled. This is normally by far the most expensive part of the calculation. Only set this to a finite value if you have a good reason to do so.

The default for value is not_checked, i.e. the convergence of the analytical stress will not be checked. However, you have to ensure that other convergence criteria (especially sc_accuracy_eev, which checks a non-variational quantity) are set to tight values, as indicated above.

Calculating the analytical stress is relatively expensive, so this convergence criterion is only checked after the purely electronic / energetic criteria, sc_accuracy_eev, sc_accuracy_etot, sc_accuracy_rho, are all fulfilled. To avoid too many iterations with analytical stress computations before the analytical stress is converged, it is important to set the energetic criteria to tight convergence.

This keyword only makes sense (and is only accepted) for periodic slab calculations with the option use_dipole_correction set. If you are interested in the work function or vacuum level shifts explicitly, it is recommended to use this flag. A typical tight convergence criterion is value=10^-4.

See the sc_init_iter keyword for a more complete description of this behavior.

sc_init_iter ends only the s.c.f. cycle for the first geometry of a run. The idea is to do a step that looks exactly like a new geometry step, but to not alter the geometry. Rather, reinitializing from the exact orbitals reached after (number) iterations ensures that the mixer and other parts of the calculation will not drag along some misguided information of the original superposition-of-free-atoms initialization. In some cases, such a clean start can help converge the s.c.f. cycle of a calculation that otherwise has difficulties to converge at all. The default is to perform this procedure only once but the user can specify how many times they want it to be performed with the option (times).

See also the sc_init_factor keyword, which can be used to tune this behavior further.

sc_iter_limit is a keyword that should be set in every run. Note: You must ensure for every run that the self-consistency cycle was actually converged. If this is not the case, a loud warning is issued in the standard output of FHI-aims at the end of the s.c.f. cycle, and relaxations, molecular dynamics, and postprocessing may continue anyway depending on the postprocess_anyway setting.

If the end of the s.c.f. cycle is reached in this way, forces are computed regardless of whether the electronic convergence was reached.

spin_mix_param may be different from charge_mix_param, but there is not usually a clear recipe how it should be different. This option is thus only needed if the standard algorithms provably fail.

This is an experimental keyword that allows to switch on an initial homogeneous_field, e.g., to lock in the symmetry of a free atom in order to enforce smoother s.c.f. convergence. The field is switched off after number iterations, before self-consistency is reached.

If type is not safe, the actual value of homogeneous_field and occupation_type should be set explicitly in geometry.in and control.in, respectively. The default homogeneous field is 10^-3 eV/Å.

If type is set to safe, a homogeneous field of 10^-3 eV/Å and a Gaussian occupation with very small (10^-5 eV) broadening are automatically enforced.

This flag is mainly used to artificially break the symmetry of spin-polarized free atoms with open $d$ and $f$ shells, which are sometimes very hard to converge otherwise (see special cases listed for mixer). Remember that FHI-aims does not allow to enforce a given symmetry automatically.

This is an experimental keyword that allows to perform a Boys localization on one or more subspaces of the KS eigenvector. Boys localization is only applicable in non-periodic calculations. Since each localization requires calculation of the transition dipole matrix, this adds a considerable overhead in computation time if it is performed in each SCF cycle. The currently recommended setting is either 0 or 2.

One can specify any xc-type from the LDA, GGA and meta-GGA options available as listed for xc. One can also use the dfauto formalism, i.e. xc_pre dfauto xc-type steps for non-hybrid XC functionals. See the dfauto section for tag xc for more detail.

Note that the elsi_restart keyword, when it has actual past restart information available to read, cannot be used together with the xc_pre keyword, unless explicitly allowed using the allow_restart_xc_pre (see there for more information). We do this to prevent unintended duelling inializations. The two keywords can still be used together if all that elsi_restart would do is write information, not read it.

The convergence analysis of the s.c.f. cycle needs sufficient data to make reliable predictions about convergence behavior. Setting scf_analysis_start adjusts the starting point of the analysis, but values below 15 s.c.f. iterations are not allowed.

This keyword controls when the convergence analysis first evaluates whether to continue or stop the s.c.f. cycle. The decision point should be set late enough that the calculation has had sufficient time to demonstrate its convergence behavior, but early enough to avoid wasting resources on poorly converging calculations.

After the initial decision point, and if the calculation is allowed to continue, additional evaluations occur at regular intervals (every scf_max_allowed_steps_at_decision iterations by default).

At each decision point, the convergence analysis predicts how many additional s.c.f. iterations are needed to reach convergence. If this prediction exceeds number, the s.c.f. cycle will be stopped (also depending on other factors such as whether convergence has stalled).

Larger values are more permissive and allow slower-converging calculations to continue. Smaller values are more aggressive in stopping calculations that are converging slowly.

The convergence analysis of the s.c.f. cycle performss number evaluations at regular intervals (every scf_max_allowed_steps_at_decision iterations) starting from scf_analysis_decision_point. If convergence remains unachievable in a reasonable amount of time through all these evaluation points, the calculation is stopped automatically at the final check point.

For example, with the default settings:

• •

  First evaluation at iteration 100 (scf_analysis_decision_point)

• •

  Second evaluation at iteration 125

• •

  Third evaluation at iteration 150

• •

  Forced stop at iteration 175 if still not converging well