# Maximum Overlap Method¶

The maximum overlap method (MOM) is used to perform variational calculations of excited states. It is an alternative to the linear expansion Delta Self-Consistent Field for obtaining excited states within a time-independent DFT framework. Since MOM calculations are variational, atomic forces are readily available from routines for ground-state calculations and can be used to perform geometry optimization and molecular dynamics.

Excited-state solutions of the SCF equations are obtained for non-Aufbau orbital occupations and can correspond to saddle points of the energy as a function of the electronic degrees of freedom (the orbital variations) 1 2 3. Hence, excited-state calculations can be affected by variational collapse to lower-energy solutions. MOM is a simple strategy to choose non-Aufbau occupation numbers consistent with the initial guess for an excited state during optimization of the wave function, thereby avoiding variational collapse.

## Implementation¶

The GPAW implementation of MOM is presented in 1 (real space grid and plane waves approaches) and 2 (LCAO approach).

The orbitals $$\{|\psi_{i}\rangle\}$$ used as initial guess for an excited-state calculation are taken as the reference orbitals for MOM (this approach is also known as initial maximum overlap method, see 4). The implementation in GPAW supports the use of fractional occupation numbers. Let $$\{|\psi_{n}\rangle\}_{s}$$ be a subspace of $$N$$ initial guess orbitals with occupation number $$f_{s}$$ and $$\{|\psi_{m}^{(k)}\rangle\}$$ the orbitals determined at iteration $$k$$ of the wave-function optimization. An occupation number of $$f_{s}$$ is given to the first $$N$$ orbitals with the biggest numerical weights, evaluated as 5:

(1)$P_{m}^{(k)} = \max_{n}\left( |O_{nm}^{(k)}| \right)$

where $$O_{nm}^{(k)} = \langle\psi_n | \psi_{m}^{(k)}\rangle$$. Alternatively, the numerical weights can be evaluated as the following projections onto the manifold $$\{|\psi_{n}\rangle\}_{s}$$ 4:

(2)$P_{m}^{(k)} = \left(\sum_{n=1}^{N} |O_{nm}^{(k)}|^{2} \right)^{1/2}$

In plane-waves or finite-difference modes, the elements of the overlap matrix are calculated from:

$O_{nm}^{(k)} = \langle\tilde{\psi}_n | \tilde{\psi}_{m}^{(k)}\rangle + \sum_{a, i_1, i_2} \langle\tilde{\psi}_n | \tilde{p}_{i_1}^{a}\rangle \left( \langle\phi_{i_1}^{a} | \phi_{i_2}^{a}\rangle - \langle\tilde{\phi}_{i_1}^{a} | \tilde{\phi}_{i_2}^{a}\rangle \right) \langle\tilde{p}_{i_2}^{a} | \tilde{\psi}_{m}^{(k)}\rangle$

where $$|\tilde{\psi}_{n}\rangle$$ and $$|\tilde{\psi}_{m}^{(k)}\rangle$$ are the pseudo orbitals, $$|\tilde{p}_{i_1}^{a}\rangle$$, $$|\phi_{i_1}^{a}\rangle$$ and $$|\tilde{\phi}_{i_1}^{a}\rangle$$ are projector functions, partial waves and pseudo partial waves localized on atom $$a$$, respectively. In LCAO, the overlaps $$O_{nm}^{(k)}$$ are calculated as:

$O_{nm}^{(k)} = \sum_{\mu\nu} c^*_{\mu n}S_{\mu\nu}c^{(k)}_{\nu m}, \qquad S_{\mu\nu} = \langle\Phi_{\mu} | \Phi_{\nu}\rangle + \sum_{a, i_1, i_2} \langle\Phi_{\mu} | \tilde{p}_{i_1}^{a}\rangle \left( \langle\phi_{i_1}^{a} | \phi_{i_2}^{a}\rangle - \langle\tilde{\phi}_{i_1}^{a} | \tilde{\phi}_{i_2}^{a}\rangle \right) \langle\tilde{p}_{i_2}^{a} | \Phi_{\nu}\rangle$

where $$c^*_{\mu n}$$ and $$c^{(k)}_{\nu m}$$ are the expansion coefficients for the initial guess orbitals and orbitals at iteration $$k$$, while $$|\Phi_{\nu}\rangle$$ are the basis functions.

## Notes on usage¶

Typically, one first performs a ground-state calculation. To prepare the calculator for an excited-state calculation, the function mom.prepare_mom_calculation can be used:

from gpaw import mom

mom.prepare_mom_calculation(calc, atoms, f)


where f contains the occupation numbers of the excited state (see examples below). Alternatively, the MOM calculation can be initialized by setting calc.set(occupations={'name': 'mom', 'numbers': f_sn}.

The default is to use eq. (1) to compute the numerical weights used to assign the occupation numbers. This was found to be more stable in the presence of diffuse virtual orbitals 5. In order to use eq. (2), instead, corresponding to the original MOM approach 4, one has to specify:

mom.prepare_mom_calculation(..., use_projections=True, ...)


SCF algorithms based on diagonalization of the Hamiltonian matrix tend to fail when degenerate or nearly degenerate orbitals are unequally occupied, a situation that is common in excited-state calculations. For such cases, it is possible to use a Gaussian smearing of the holes and excited electrons in the MOM calculation to improve convergence. This is done by specifying a width in eV (e.g. width=0.01) for the Gaussian smearing function. For difficult cases, the width can be increased at regular intervals to force convergence by specifying a width_increment=.... Note, however, that too extended smearing can lead to discontinuities in the potentials and forces close to crossings between electronic states 2, so this feature should only be used at geometries far from state crossings.

gpaw.mom.prepare_mom_calculation(calc, atoms, numbers, use_projections=False, update_numbers=True, use_fixed_occupations=False, width=0.0, niter_width_update=40, width_increment=0.0)[source]

Helper function to prepare a calculator for a MOM calculation.

calc: GPAW instance

GPAW calculator object.

atoms: ASE instance

ASE atoms object.

numbers: list (len=nspins) of lists (len=nbands)

Occupation numbers (in the range from 0 to 1). Used for the initialization of the MOM reference orbitals.

use_projections: bool

If True, the occupied orbitals at iteration k are chosen as the orbitals |psi^(k)_m> with the biggest weights P_m evaluated as the projections onto the manifold of reference orbitals |psi_n>: P_m = (Sum_n(|O_nm|^2))^0.5 (O_nm = <psi_n|psi^(k)_m>) see https://doi.org/10.1021/acs.jctc.7b00994. If False (default), the weights are evaluated as: P_m = max_n(|O_nm|) see https://doi.org/10.1021/acs.jctc.0c00488.

update_numbers: bool

If True (default), ‘numbers’ gets updated with the calculated occupation numbers, and when changing atomic positions the MOM reference orbitals will be initialized as the occupied orbitals found at convergence for the previous geometry. If False, when changing positions the MOM reference orbitals will be initialized as the orbitals of the previous geometry corresponding to the user-supplied ‘numbers’.

use_fixed_occupations: bool

If True (default), the MOM algorithm is used. If False, fixed occupations will be used.

width: float

Width of Gaussian function in eV for smearing of holes and excited electrons. The holes and excited electrons are found with respect to the zero-width ground-state occupations. See https://doi.org/10.1021/acs.jctc.0c00597.

niter_width_update: int

Number of iterations after which the width of the Gaussian smearing function is increased.

width_increment: float

How much to increase the width of the Gaussian smearing function.

## Example I: Excitation energy molecular Rydberg states¶

In this example, the excitation energies of the singlet and triplet states of water corresponding to excitation from the HOMO-1 non-bonding ($$n$$) to the LUMO $$3s$$ Rydberg orbitals are calculated. The $$n$$ and $$3s$$ orbitals have the same symmetry (a1), therefore, variational collapse can potentially affect a calculation without MOM. In order to calculate the energy of the open-shell singlet state, first a calculation of the mixed-spin state obtained for excitation within the same spin channel is performed, and, then, the spin-purification formula is used: $$E_s=2E_m-E_t$$, where $$E_m$$ and $$E_t$$ are the energies of the mixed-spin and triplet states, respectively.

import copy
from ase.build import molecule
from ase.parallel import paropen
from gpaw import GPAW
from gpaw.mom import prepare_mom_calculation

atoms = molecule('H2O')
atoms.center(vacuum=7)

calc = GPAW(mode='fd',
basis='dzp',
nbands=6,
h=0.2,
xc='PBE',
spinpol=True,
symmetry='off',
convergence={'bands': -1},
txt='h2o.txt')
atoms.calc = calc

# Ground-state calculation
E_gs = atoms.get_potential_energy()

# Ground-state occupation numbers
f_gs = []
for s in range(2):
f_gs.append(calc.get_occupation_numbers(spin=s))

# Triplet n->3s occupation numbers
f_t = copy.deepcopy(f_gs)
f_t -= 1.  # Remove one electron from homo-1 (n) spin up
f_t += 1.  # Add one electron to lumo (3s) spin down

# MOM calculation for triplet n->3s state
prepare_mom_calculation(calc, atoms, f_t)
E_t = atoms.get_potential_energy()

# Mixed-spin n->3s occupation numbers
f_m = copy.deepcopy(f_gs)
f_m -= 1.  # Remove one electron from homo-1 (n) spin up
f_m += 1.  # Add one electron to lumo (3s) spin up

# MOM calculation for mixed-spin n->3s state
prepare_mom_calculation(calc, atoms, f_m)
E_m = atoms.get_potential_energy()
E_s = 2 * E_m - E_t  # Spin purified energy

with paropen('h2o_energies.txt', 'w') as fd:
print(f'Excitation energy triplet n->3s state: {E_t - E_gs:.2f} eV',
file=fd)
print(f'Excitation energy singlet n->3s state: {E_s - E_gs:.2f} eV',
file=fd)
# https://doi.org/10.1021/acs.jctc.8b00406
print('Experimental excitation energy triplet n->3s state: 9.46 eV',
file=fd)
print('Experimental excitation energy singlet n->3s state: 9.67 eV',
file=fd)


## Example II: Excited-state geometry relaxation¶

Here, the bond length of the carbon monoxide molecule is optimized in the singlet excited state obtained by promotion of an electron from the HOMO $$\sigma$$ orbital to the LUMO $$\pi^*_x$$ or $$\pi^*_y$$ orbital. A practical choice for geometry optimization and dynamics in an open-shell singlet excited state is to employ a spin-paired approach where the occupation numbers of the open-shell orbitals are set to 1 6. This approach delivers pure singlet states while avoiding an additional calculation of the corresponding triplet state needed to employ the spin-purification formula (see Example I: Excitation energy molecular Rydberg states). Since the $$\pi^*_x$$ and $$\pi^*_y$$ orbitals of carbon monoxide are degenerate, diagonalization-based SCF algorithms fail to converge to the $$\sigma\rightarrow\pi^*$$ excited state unless symmetry constraints on the density are used. Here, Gaussian smearing of the excited electron is used to force equal fractional occupations of the two $$\pi^*$$ orbitals to avoid convergence issues.

from ase.build import molecule
from ase.optimize import LBFGS
from ase.parallel import paropen
from gpaw import GPAW
from gpaw.mom import prepare_mom_calculation

atoms = molecule('CO')
atoms.center(vacuum=5)

calc = GPAW(mode='lcao',
basis='dzp',
nbands=8,
h=0.2,
xc='PBE',
symmetry='off',
convergence={'density': 1e-7},
txt='co.txt')
atoms.calc = calc

# Ground-state calculation
E_gs = atoms.get_potential_energy()

# Ground-state occupation numbers
f = [calc.get_occupation_numbers(spin=0) / 2.]

# Singlet sigma->pi* occupation numbers
f -= 0.5  # Remove one electron from homo (sigma)
f += 0.5  # Add one electron to lumo (pi*)

# Excited-state MOM calculation with Gaussian
# smearing of the occupation numbers
prepare_mom_calculation(calc, atoms, f, width=0.01)

opt = LBFGS(atoms, logfile='co.log')
opt.run(fmax=0.05)

d = atoms.get_distance(0, 1)

with paropen('co.log', 'a') as fd:
print(f'Optimized C-O bond length of sigma->pi* state: {d:.2f} Å', file=fd)
# https://doi.org/10.1007/978-1-4757-0961-2
print('Experimental C-O bond length of sigma->pi* state: 1.24 Å', file=fd)