Maximum Overlap Method

Excited-state solutions of the SCF equations are obtained for non-Aufbau orbital occupations. 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 facilitating convergence to the target excited state and avoiding variational collapse to lower energy solutions.

The MOM approach implemented in GPAW is the initial maximum overlap method [1], where the orbitals used as initial guess for an excited-state calculation are taken as fixed reference orbitals. The GPAW implementation is presented in [2] (real-space grid and plane-wave approaches) and [3] (LCAO approach), and supports the use of fractional occupation numbers.

Excited-state calculations can be difficult to convergence because excited states typically correspond to saddle points of the electronic energy surface, and not minima. The recommended approach to use together with MOM for calculating excited states is direct optimization (DO) (see Excited-State Calculations with Direct Optimization), which is an alternative to the diagonalization-based eigensolvers and is better suited for converging on saddle points.

Maximizing wavefunction overlaps

Let \(\{|\psi^0_{n}\rangle\}\) be the set of reference orbitals with occupation numbers \(f_n^0\) and \(\{|\psi_{m}^{(k)}\rangle\}\) the orbitals determined at iteration \(k\) of the wave-function optimization. The method aims to find the updated occupation numbers \(f_m^{(k)}\) for the orbitals at iteration \(k\) such to minimize the electronic distance \(\eta\) defined as:

\[\eta = N - \sum_{n m} \bigl|\langle \psi^0_{n} \mid \psi_{m}^{(k)} \rangle\bigr|^2\]

where \(N\) is the number of electrons and the summations run over the occupied orbitals only.

Naively, this can be achieved by finding a mapping between the orbitals using the wavefunction overlap \(O_{nm}^{(k)} = \langle\psi^0_n | \psi_{m}^{(k)}\rangle\), as a measure of their similarity.

Tip

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

\[O_{nm}^{(k)} = \langle\tilde{\psi}^0_n | \tilde{\psi}_{m}^{(k)}\rangle + \sum_{a, i_1, i_2} \langle\tilde{\psi}^0_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}^0_{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^{*0}_{\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^{*0}_{\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.

Effectively, we want to find the permutation \(\mathcal P\) of the updated occupation numbers such that the sum of the absolute values of the diagonal elements of the overlap matrix is maximized:

\[\max_{\mathcal P} \sum_{nm} {\mathcal P} _{nm}\,\bigl|O^{(k)}_{nm}\bigr| = \max_{\mathcal P} \text{Tr} \left( \mathcal P |O^{(k)}| \right) \rightarrow \mathcal P^\max\]

where \({\mathcal P}_{nm} \in \{0,1\}\). With the matrix representation of the permutation \(\mathcal P^\max_{nm}\) the occupation numbers are updated according to:

\[f_m^{(k)} = \sum_n \mathcal P^\max_{nm} f_n^0\]

Given the wavefunction overlaps \(|O_{nm}^{(k)}|\), the optimal permutation can be found using scipy.optimize.linear_sum_assignment. The figure shows the absolute values of the overlap matrix for a fictional system with 8 bands and the initial occupations \(f_n^0\).

From the visual representation of the overlap matrix in the left panel of the figure above it is immediately clear how the updated occupations \(f_m^{(k)}\) should look like.

Maximizing subspace projections

An alternative approach consists in maximizing the sum of the projections of the orbitals onto subspaces of equally occupied reference orbitals, \({s(f_s)}\):

\[P_{m}^{(k)}(f_s) = \left(\sum_{n \in s} |O_{nm}^{(k)}|^{2} \right)^{1/2}\]

where \(n \in s\) denotes that only orbitals from the subspace \({s(f_s)}\) are taken into account. Thus, one maximizes the sum:

\[\sum_{m} P_{m}^{(k)}(f_m^{(k)})\]

If only one subspace of equally occupied orbitals is present, the method is equivalent to the initial maximum overlap method presented in [1].

The right panel in the figure above shows the weight matrix calculated from the overlaps of the above example. Again, the assignment, which for more than one subspace is done using scipy.optimize.linear_sum_assignment numerically, can be directly seen in this example case.

How to use MOM

Initial guess orbitals for the excited-state calculation are first needed. Typically, they are obtained from a ground-state calculation. Then, to prepare the calculator for a MOM 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}. A helper function can be used to create the list of excited-state occupation numbers:

from gpaw.directmin.tools import excite
f = excite(calc, i, a, spin=(si, sa))

which will promote an electron from occupied orbital i in spin channel si to unoccupied orbital a in spin channel sa (the index of HOMO and LUMO is 0). For example, excite(calc, -1, 2, spin=(0, 1)) will remove an electron from the HOMO-1 in spin channel 0 and add an electron to LUMO+2 in spin channel 1.

The default is to use the subspace projections to assign the occupation numbers. In order to use the approach of matching the orbitals directly based on their overlaps, one has to specify:

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

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

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 (default), 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 doi: 10.1021/acs.jctc.7b00994. If False, try to match the orbitals directly based on their overlaps.

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 False (default), the MOM algorithm is used. If True, 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 doi: 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.

Excitation energy Rydberg state of water

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. 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 [4] 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. The calculations use the Finite Difference mode to obtain an accurate representation of the diffuse Rydberg orbital [2].

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[0][2] -= 1.  # Remove one electron from homo-1 (n) spin up
f_t[1][4] += 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[0][2] -= 1.  # Remove one electron from homo-1 (n) spin up
f_m[0][4] += 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)

Geometry relaxation excited-state of carbon monoxide

In this example, the bond length of the carbon monoxide molecule in the lowest singlet \(\Pi(\sigma\rightarrow \pi^*)\) excited state is optimized using two types of calculations, each based on a different approximation to the potential energy curve of an open-shell excited singlet state. The first is a spin-polarized calculation of the mixed-spin state as defined in Excitation energy Rydberg state of water. The second is a spin-paired calculation where the occupation numbers of the open-shell orbitals are set to 1 [5]. Both calculations use LCAO basis and the direct optimization (DO) method.

In order to obtain the correct angular momentum of the excited state, the electron is excited into a complex \(\pi^*_{+1}\) or \(\pi^*_{-1}\) orbital, where +1 or −1 is the eigenvalue of the z-component angular momentum operator. The use of complex orbitals provides an excited-state density with the uniaxial symmetry consistent with the symmetry of the molecule [6].

from ase.build import molecule
from ase.optimize import BFGS
from ase.parallel import paropen
from gpaw import GPAW, LCAO
from gpaw.mom import prepare_mom_calculation
from gpaw.directmin.tools import excite
from gpaw.directmin.etdm_lcao import LCAOETDM


for spinpol in [True, False]:
    if spinpol:
        tag = 'spinpol'
    else:
        tag = 'spinpaired'

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

    calc = GPAW(xc='PBE',
                mode=LCAO(force_complex_dtype=True),
                h=0.2,
                basis='dzp',
                spinpol=spinpol,
                eigensolver='etdm-lcao',
                occupations={'name': 'fixed-uniform'},
                mixer={'backend': 'no-mixing'},
                nbands='nao',
                symmetry='off',
                txt='co_' + tag + '.txt')
    atoms.calc = calc

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

    # Prepare initial guess for complex pi* orbitals by taking
    # linear combination of real pi*x and pi*y orbitals
    lumo = 5  # lumo is pi*x or pi*y orbital
    for kpt in calc.wfs.kpt_u:
        pp = kpt.C_nM[lumo] + 1.0j * kpt.C_nM[lumo + 1]
        pm = kpt.C_nM[lumo] - 1.0j * kpt.C_nM[lumo + 1]
        kpt.C_nM[lumo][:] = pm
        kpt.C_nM[lumo + 1][:] = pp

    calc.set(eigensolver=LCAOETDM(searchdir_algo={'name': 'l-sr1p'},
                                  linesearch_algo={'name': 'max-step'},
                                  need_init_orbs=False))

    # Occupation numbers for sigma->pi* excited state:
    # Remove one electron from homo (sigma) and add one electron to lumo (pi*)
    f = excite(calc, 0, 0, spin=(0, 0))
    if not spinpol:
        f[0] /= 2

    # Prepare excited-state DO-MOM calculation
    prepare_mom_calculation(calc, atoms, f)

    opt = BFGS(atoms, logfile='co_' + tag + '.log', maxstep=0.05)
    opt.run(fmax=0.05)

    d = atoms.get_distance(0, 1)

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

The electronic configuration of the \(\Pi(\sigma\rightarrow \pi^*)\) state includes two unequally occupied, degenerate \(\pi^*\) orbitals. Because of this, convergence to this excited state is more difficult when using SCF eigensolvers with density mixing instead of DO, unless symmetry constraints on the density are enforced during the calculation. Convergence of such excited-state calculations with an SCF eigensolver can be improved by using a Gaussian smearing of the holes and excited electrons [5]. Gaussian smearing is implemented in MOM and can be used by specifying a width in eV for the Gaussian smearing function:

mom.prepare_mom_calculation(..., width=0.01, ...)

For difficult cases, the width can be increased at regular intervals 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 [7], so this feature should be used with caution and only at geometries far from state crossings.

References

[1] (1,2)

G. M. J. Barca, A. T. B. Gilbert, P. M. W. Gill Simple Models for Difficult Electronic Excitations, J. Chem. Theory Comput., 14 1501-1509 (2018).

[2] (1,2)

A. V. Ivanov, G. Levi, H. Jónsson Method for Calculating Excited Electronic States Using Density Functionals and Direct Orbital Optimization with Real Space Grid or Plane-Wave Basis Set, J. Chem. Theory Comput., (2021).

[3]

G. Levi, A. V. Ivanov, H. Jónsson Variational Density Functional Calculations of Excited States via Direct Optimization, J. Chem. Theory Comput., 16 6968–6982 (2020).

[4]

T. Ziegler, A. Rauk, E. J. Baerends On the calculation of multiplet energies by the hartree-fock-slater method Theoret. Chim. Acta, 43 261–271 (1977).

[5] (1,2)

G. Levi, M. Pápai, N. E. Henriksen, A. O. Dohn, K. B. Møller Solution structure and ultrafast vibrational relaxation of the PtPOP complex revealed by ∆SCF-QM/MM Direct Dynamics simulations, J. Phys. Chem. C, 122 7100-7119 (2018).

[6]

A. V. Ivanov, G. Levi, H. Jónsson Method for Calculating Excited Electronic States Using Density Functionals and Direct Orbital Optimization with Real Space Grid or Plane-Wave Basis Set, J. Chem. Theory Comput., (2021).

[7]

G. Levi, A. V. Ivanov, H. Jónsson Variational Density Functional Calculations of Excited States via Direct Optimization, J. Chem. Theory Comput., 16 6968–6982 (2020).