Excited-State Calculations with Maximum Overlap Method and Direct Optimization

The maximum overlap method (MOM) can be 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 the method get_forces and can, therefore, be used to perform geometry optimization and molecular dynamics in the excited state.

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.

Even if MOM is used, an excited-state calculation can still be difficult to convergence with the SCF algorithms based on diagonalization of the Hamiltonian matrix that are commonly employed in ground-state calculations. One of the main problems is that excited states often correspond to saddle points of the energy as a function of the electronic degrees of freedom (the orbital variations), but these algorithms perform better for minima (ground states usually correspond to minima). Moreover, standard SCF algorithms tend to fail when degenerate or nearly degenerate orbitals are unequally occupied, a situation that is more common in excited-state rather than ground-state calculations (see Example II: Geometry relaxation excited-state of carbon monoxide below). In GPAW, excited-state calculations can be performed via a direct optimization (DO) of the orbital (implemented for the moment only in LCAO). DO can converge to a generic stationary point, and not only to a minimum and has been shown to be more robust than diagonalization-based SCF algorithms using density mixing in excited-state calculations of molecules [1] [2] [3]; therefore, it is the recommended method for obtaining excited-state solutions with MOM.

Maximum overlap method

Implementation

The MOM approach implemented in GPAW is the initial maximum overlap method [4]. The implementation 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 fixed reference orbitals for MOM. 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.

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 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, ...)

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)

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 doi: 10.1021/acs.jctc.7b00994. If False (default), the weights are evaluated as: P_m = max_n(|O_nm|) see doi: 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 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.

Direct optimization

Direct optimization (DO) can be performed using the implementation of exponential transformation direct minimization (ETDM) [1] [2] [3] described in Direct Minimization Methods. This method uses the exponential transformation and efficient quasi-Newton algorithms to find stationary points of the energy in the space of unitary matrices. Currently, DO can be performed only in LCAO mode.

For excited-state calculations, the recommended quasi-Newton algorithm is a limited-memory symmetric rank-one (L-SR1) method [2] with unit step. In order to use this algorithm, the following eigensolver has to be specified:

from gpaw.directmin.etdm import ETDM

calc.set(eigensolver=ETDM(searchdir_algo={'name': 'l-sr1p'},
                          linesearch_algo={'name': 'max-step',
                                           'max_step': 0.20})

The maximum step length avoids taking too large steps at the beginning of the wave function optimization. The default maximum step length is 0.20, which has been found to provide an adequate balance between stability and speed of convergence for calculations of excited states of molecules [2]. However, a different value might improve the convergence for specific cases.

Example I: 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 [6] 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 [1].

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)

Example II: 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 Example I: 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 [7]. 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 [1].

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 import ETDM


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',
                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=ETDM(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 [7]. 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 [2], so this feature should be used with caution and only at geometries far from state crossings.

References

[1] (1,2,3,4,5)

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).

[2] (1,2,3,4,5,6)

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).

[3] (1,2)

G. Levi, A. V. Ivanov, H. Jónsson Variational Calculations of Excited States Via Direct Optimization of Orbitals in DFT, Faraday Discuss., 224 448-466 (2020).

[4] (1,2,3)

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).

[5] (1,2)

X. Dong, A. D. Mahler, E. M. Kempfer-Robertson, L. M. Thompson Global Elucidation of Self-Consistent Field Solution Space Using Basin Hopping, J. Chem. Theory Comput., 16 5635−5644 (2020).

[6]

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).

[7] (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).