3.50 Localisation of Molecular Orbital with Optimised Unitary Transformations

SCF methods utilised in most quantum chemistry codes rely on matrix diagonalisation to produce a set of eigenvectors orthogonalised against the one-electron Hamiltonian (ie., the Kohn-Sham eigenvalue problem). However, the eigenvectors (Kohn-Sham states) produced by this method are often highly delocalised across the molecular system, taking the form of Canonical Molecular Orbitals (CMOs) that reflect the irreducible representation of the system’s underlying symmetry [31]. Although density, and in turn total energies, are invariant under unitary tranformations of the MOs, the delocalised representation is counterintuitive from a chemical perspective [294], which is focussed on bonding interactions between two or more atoms. Furthermore, approaches such as QM-in-QM embedding [192] must solve distinct eigenvalue problems and calculate properties with differing levels of theory for given subregions. To achieve this, the subregions of electron density and their corresponding molecular orbitals must be well-defined and local for the atoms of interest; delocalised representations of electron density for each state are ill-suited for this purpose.

Refer to caption — Figure 3.15: Pipek-Mizey localisation applied to 1,3,6-octatriene using the unitary optimisation of localised molecular orbitals method

However, by performing an appropriate unitary transformation of the KS eigenvectors, one may produce an energetically equivalent set of MOs with greater localisation,

\displaystyle\Psi_{i}\rightarrow\sum_{i}\mathbf{W}_{ji}\Psi_{j}=\Psi_{{i}^{{}^% {\prime}}}

(3.274)

where $\Psi_{i}$ is the ith molecular orbital, $\mathbf{W}_{ji}$ is an orthogonal unitary transformation matrix, and $\Psi_{{i}^{{}^{\prime}}}$ is the ith molecular orbital produced by rotations with respect to the full set of eigenvectors, $\Psi_{j}$ . The calculation of $\mathbf{W}_{ji}$ is based on the unitary optimisation of localised molecular orbitals approach developed by Lehtola et al. [193]. This approach maximises/minimises a given objective (cost) function $\mathcal{J(\mathbf{W})}$ , which quantifies the degree of delocalisation for the rotated set of MOs for a given metric. The current implementation supports the Pipek-Mizey localisation method [252], which calculates the objective function through the expression,

\displaystyle\mathcal{J}^{PM}(\mathbf{W})=\sum^{N_{at}}_{A=1}\sum^{N_{states}}% _{i=1}[Q_{ii}^{A}]^{2}

(3.275)

where $Q^{A}$ is the state-projected population matrix for the states on atom $A$ . The trace of the population matrix for occupied states yields the total charge centered on a given atom, and the exponent biases distributing the charge to a smaller subset of states on a given atom. Minimisation of the cost function is computed via conjugate gradients under constraints which ensure the matrix $\mathbf{W}$ remains unitary [1] [2].

The state population matrix $Q^{A}$ may be calculated using any charge population analysis that projects onto states. The current implementation uses Mulliken and Hirshfeld based analysis, with the latter preferred as it avoids the basis set dependency found in Mulliken analysis.

Localisation optimisation as implemented can produce more chemically intuitive of the KS eigenvectors for visualisation and further analysis (output eigenvectors). However, the real power of this approach is apparent when paired with features of the code which intrinsically require localisation of the Kohn-Sham states or well-defined spatial separation of charge for different states. Such work is underway for QM-in-QM embedding and $\Delta$ SCF method.

Currently only supports the aperiodic systems, but periodic support is a possible future development.

A note on restarting from rotated KS eigenvectors: Restarting from localised KS eigenvectors is supported through the lmo_in_scf keyword in conjunction with restart. It is recommended to set lmo_in_scf with freq equal to the intervals used to write the KS eigenvectors for restarting.

A note on verbose outputting: Setting output_level full will provide additional (arguably too much) information on the cost function, line search and convergence threshold at each iteration in the optimisation process. This can be used to check whether convergence is being achieved, or whether the algorithm is struggling for your system.

Tags for general section of control.in:

Tag: lmo_use(control.in)

Usage: lmo_use .true.
Purpose: Switches on MO localisation during postprocessing after exiting the SCF loop. KS eigenvectors are converted to their localised form. Default: .false.

Tag: lmo_in_scf(control.in)

Usage: lmo_in_scf freq
Purpose: Switches on MO localisation within the SCF loop. Localisation is performed after the KS eigenvalue problem is evaluated. KS eigenvectors are converted to their localised form.

freq is an integer, specifying how frequently localisation is performed in the SCF loop.

Note: Unless otherwise needed for other features, we do not recommend using this keyword - the KS eigensolver produces a new set of delocalised KS eigenvectors at every iteration in the SCF loop. Unless the localised KS eigenvectors are evaluated at each SCF iteration for a specific application, localisation adds unecessary work. Use lmo_use if you only care about outputting the localised form of your converged KS eigenvectors.

Tag: lmo_restart_write(control.in)

Usage: lmo_restart_write .false.
Purpose: Applies localisation to KS eigenvectors output in the keyword restart. restart must be set, otherwise FHI-aims will stop. Default: .false.

Tag: lmo_max_iter(control.in)

Usage: lmo_max_iter lmo_max_iter
Purpose: Limits the maximum number of iterations used in the optimisation process. Default: 1000

lmo_max_iter is an integer.

If the maximum number of iterations is reached without convergence, two outcomes may occur: if the cost function is ’stuck’ at a saddle-point (not significantly changing from one iteration to the next but failing to converge) a warning will display and the calculation will continue as normal; alternatively, if the cost function is still changing from one iteration to the next and convergence is not achieved, FHI-aims will stop. In either scenario, we suggest setting output_level full to analyse the rate of convergence. If the convergence criteria is still reducing, setting lmo_max_iter to a higher value may resolve your issue.

Tag: lmo_unocc(control.in)

Usage: lmo_unocc .true.
Purpose: Limits localisation to occupied KS eigenvectors only. Default: .true.

Inclusion of the virtual KS eigenvectors introduces additional overhead by introducing additional states to evaluate. More importantly, the rate of convergence for the optimisation process slows down drastically. Therefore the default .true. is recommended unless the unoccupied states are desired.

Tag: lmo_objective_func(control.in)

Usage: lmo_objective_func objective_func
Purpose: The localisation metric used to calculate the objective function. Only Pipek-Mizey (pm) is supported. Default: pm

objective_func is a string. Anything other than pm will exit the simulation.

Tag: lmo_pm_charge_metric(control.in)

Usage: lmo_pm_charge_metric pm_charge_metric
Purpose: Sets the charge metric used in the Pipek-Mizey localisation method. Default: hirshfeld

objective_func is a string, which can be either mulliken or hirshfeld.

We recommend sticking with Hirshfeld unless there is a specific reason to believe Mulliken charge analysis would be preferable.

Tag: lmo_pm_charge_metric(control.in)

Usage: lmo_init_guess init_guess
Purpose: Sets the initial guess of the unitary transformation matrix, $\mathbf{W}$ . Takes the form of orthogonal matrix, $\mathbf{W^{{\dagger}}}\mathbf{W}=\mathbf{I}$ . Default: canon

init_guess is a string, which can be either canon or random.

The two available starting points for $\mathbf{W}$ are:

•

canon : Sets $\mathbf{W}=\mathbf{I}$ i.e., the starting cost function is calculated using the input KS eigenvectors.
•

random : Sets $\mathbf{W}$ to a random orthogonal matrix. This starting point breaks the symmetry of the KS eigenvectors, which can avoid saddle points in the optimisation algorithm. However, this can make testing problematic and may produce differing solutions to the optimisation problem.

Tag: lmo_update_factor(control.in)

Usage: lmo_update_factor update_factor
Purpose: Sets the method for calculating the step direction to either steepest-descent or conjugate gradient. Default: cgpr

update_factor is a string, which can be either sd or cgpr.

The two available starting points for $\mathbf{W}$ are:

•

sd : Calculates the step direction as a simple steepest-descent.
•

cgpr : Calculates the step direction as a conjugate gradient in the Polak-Ribiére-Polyak form.

Tag: lmo_debug_matrix_output(control.in)

Usage: lmo_debug_matrix_output .false
Purpose: Outputs the final unitary transformation matrix, $\mathbf{W}$ . Used mostly for CI tests and debugging. Default: .false.

Tag: lmo_force_sign_conv(control.in)

Usage: lmo_force_sign_conv .false.
Purpose: Used to counter a bug in the CI tests, where different compilations of FHI-aims calculated KS eigenvectors with differing phase factors. Do not use for general calculations. Default: .false.