3.14 Electronic constraints

To divide up space into regions for the purpose of enforcing an electronic constraint, each atom in the structure is included in a constraint_region.

Simple example of an Na-Cl dimer (geometry.in):

  atom 0. 0. 0. Na
    constraint_region 1
  atom 0. 0. 3. Cl
    constraint_region 2

assigns Na to the first region and Cl to the second region of a constrained calculation.

The special case of only one region (e.g., for a fixed spin moment calculation) needs no explicit constraint_region labels. Note that, apart from an explicit setup by keyword constraint_electrons, the case of an integer fixed spin moment for the whole system (all atoms) can also be called by the shortcut multiplicity.

This is the central keyword that can be used to define a strict constraint on the electron numbers associated (i) with the orbitals of a given subset of atoms (“region”) and / or (ii) with a given spin channel. See the multiplicity keyword for a shortcut for fixed spin moment calculations with an integer spin multiplicity.

The method to determine the constraint potentials that enforce the local electron / spin constraint is set by constraint_mix; for more than one active constraint_region, this determination is always iterative. Keyword constraint_it_lim sets the maximum number of iterations before this search is aborted in case of failed convergence (or too ambitious accuracy requirements from keyword constraint_precision).

Only meaningful for non-standard settings of mixer_constraint, irrelevant for the standard bfgs case.

Only meaningful for non-standard settings of mixer_constraint, irrelevant for the standard bfgs case.

This flag has nothing to do with electron density mixing or the electronic self-consistency loop. Instead, this defines the process to determine constraint potentials that enforce the requested electron number constraints, if the keyword constraint_electrons was invoked. This process happens in each s.c.f. iteration after the Hamiltonian matrix is known, in the course of solving the Kohn-Sham eigenvalue problem.

If more than one constraint regions are requested, determining the constraint potentials to enforce the local constraint on the electron numbers is an iterative process. Technically, FHI-aims supports three different algorithms for type :

• •

  bfgs : The default. A BFGS algorithm that optimizes the constraint potentials to their nearest local optimum. Nothing else should be used unless for experimental purposes.

• •

  linear : Linear mixing algorithm to determine the constraint potentials.

• •

  pulay : Pulay-type mixing algorithm to determine the constraint potentials.

Under normal circumstances, the mixer_constraint keyword should not be needed explicitly.

Only meaningful for non-standard settings of mixer_constraint, irrelevant for the standard bfgs case.

Example:
force_occupation_basis 1 1 atomic 2 1 1 1.3333 6

This choice will constrain the overall occupation of a given basis function (not Kohn-Sham state!) in the system to 1.3333 electrons.

The basis function in question resides on the first atom (number 1) as listed in geometry.in. The first spin channel is constrained.

Since we are interested in constraining an atomic-like orbital, we choose one that is part of the minimal basis (type “atomic”).

In fact, we let the constraint act on a 2$p$ level, $m$=1 (defined by the sequence “2 1 1”).

Only Kohn-Sham orbitals up to the 6th state (counted in the overall system) will be included in the constraint. In general, it is a good idea to constrain the orbital in question out of the occupied space, i.e., choose a value for max_KS_state that indicates a state above the Fermi level.

There is a small problem here: We need to define the occupation of a “basis function,” but really, we here have a non-orthogonal basis. Strictly speaking, only the Kohn-Sham states in the system have a well-defined occupation. What to do?

One thing we could use (and we do) are Mulliken occupation numbers of the basis functions. These are formally always well defined. In practice, however, as the overall basis set becomes larger and larger and approaches (over-)completeness, Mulliken occupations become less and less meaningful because other basis functions are not exactly orthogonal to the one used in our projection, and can “restore” the component that was originally constrained away.

Either way, we use Mulliken occupations, assuming that the atomic core basis functions are sufficiently compact and practically orthogonal to everything else.

This assumption will work well for localized basis functions such as the 1$s$ levels of most elements. As a rule, the constraint is expected to be less and less unique if applied to more delocalized basis functions – there can even be multiple different self-consistent constrained solutions for the same formal constraint. For instance, Si 2$p$ basis functions can exhibit this problem if the basis sets on the surrounding atoms – which overlap with the Si atom – become too large. Here, a smaller basis set (tier 1) can indeed be the way to keep a qualitatively meaningful Mulliken-type constraint.

Example:
force_occupation_projector 8 1 0.0 6 10
force_occupation_projector 9 2 1.0 6 10

This enforces 0.0 occupation in state 8 of spin channel 1 and 1.0 occupation for state 9 of spin channel 2. KS states between 6 and 10 are checked for overlap with state 8 and 9 of previous SCF steps. If 8/1 was occupied and 9/2 was unoccupied in the ground state, this corresponds to a triplet excitation 8$\to$9.

To simulate XPS energies with $n$ inequivalent atoms of the same species (called excitation centre in the following) a total of $n+1$ single runs is required: One ground state calculation and one force_occupation_basis/force_occupation_projector calculation for each of the excitation centres.

The ground state calculation should use the restart_write_only or the restart flag to create a restart file that is needed for the force_occupation_basis or alternatively the force_occupation_projector run. Therefore the geometry.in and all other parameters (except the charge) such as basis sets have to match. In practise it is often beneficial for the interpretation of the XP spectra to use output cube eigenstate to have an idea of the localization of the different core levels.

For the simulation of the XPS spectra typically a full core hole is introduced at the excitation centre (for example in ref. [75]). Relying on initial state effects alone (i.e., defining the ionization energies using ground state eigenvalues) neglects the screening of the core hole by valence electrons [203] and is therefore not a good approximation for XPS. For example, to force the occupation of eigenstate 3 to 1.0, where states 1-4 are (near-)degenerate or at least very similar in energy and type:

force_occupation_projector 3 1 1.0 1 4

Results in the following occupation (the introduction of the core hole leads to a re-ordering):

Note that charge was set to 1 to take into account the reduced electron number and a restart from the ground state run was made using restart.

XPS energies can then be calculated as the difference of the total energy obtained in the ground state calculation and the total energy of the core-hole excited simulation (corresponding to the definition of the ionization energy). That means that for each excitation center an ionization energy is calculated. For the core level shifts only relative energy differences are relevant, which are already directly reflected in the differences of total energies of the core-hole excited states. If, however, absolute energies are of interest, note that experiments are referenced to either the vacuum level or the Fermi level, and that simulations including an extended surface might differ by the workfunction from those for isolated molecules. The ionization energies can then by broadened with Gaussian functions of same amplitude (assuming no preferential direction, especially valid for $1s$ spectra) and summed up to obtain the total XPS spectrum.

Simulating NEXAFS spectra can be less straightforward as there are different approximations to account for the core hole and the excited electron. One possibility is to use the transition potential approach [301], where instead of a full core only half a core hole is used, i.e., $n=0.5$ in one spin channel. Independently from the approximation used for the core hole: To obtain the dipole matrix elements that give information about the transition probability the flag compute_dipolematrix needs to be used. Note that to use this option the FHI-aims binary has to be compiled enabling hdf5, as the output is a hdf5 container containing eigenvalues and matrix elements. It is recommended to include additional empty_states, depending on the amount of unoccupied states you want to probe. In this case the ground state calculation has already to include the same number of empty states, otherwise a restart is not possible.

Example:
force_occupation_smearing 0.05

This keyword helps to converge systems with state degeneracies, which are constrained by force_occupation_projector. Specifically when calculating electronic excited states in the frontier orbital regime, many state degeneracies can occur. If one of two degenerate states is constrained to a different than the ground state occupation, convergence can be hindered. If this happens, this keyword can enable convergence for the price of a minimally different final occupation and therefore also a small error in excitation energy. One should be very careful with this keyword and only employ it if convergence can not be reached without.

WARNING: It is very easy to generate reandom numbers when using this keyword. The smearing value should never exceed 0.15