Excited-State Calculations with Direct Optimization

Direct optimization (DO) is an alternative to the diagonalization-based eigensolvers, which does not use density mixing. Since it is designed to converge to saddle point on the electronic energy surface, this approach can be used to perform variational calculations of excited states. The atomic forces can be obtained directly from the method get_forces and therefore, the method can be used to perform geometry optimization and molecular dynamics in the excited state.

The implementation of DO is based on the exponential transformation (see also Direct Minimization Methods) and uses efficient quasi-Newton algorithms for saddle point convergence. The real-space grid and plan-wave implementation is described in [1], while the LCAO inplementation is described in [2] and [3].

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.lcao_etdm import LCAOETDM

calc.set(eigensolver=LCAOETDM(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.

To reduce the risk of variational collapse to lower energy solutions, DO is typically used together with the Maximum Overlap Method. However, if the target excited state shows pronounced charge transfer character, variational collapse can sometimes not be prevented even if DO and MOM are used in conjunction. In such cases, it can be worthwhile to first perform a constrained optimization in which the electron and hole orbitals involved in the target excitation are frozen, and a minimization is done in the remaining subspace, before performing a full unconstrained optimization. The constrained minimization takes care of a large part of the prominent orbital relaxation effect in charge transfer excited states and thereby significantly simplifies the subsequent saddle point search, preventing variational collapse. Constrained optimization can be performed by using the constraints keyword:

calc.set(eigensolver=LCAOETDM(constraints=[[[h11], [h12],..., [p11], [p12],...],
                                           [[h21], [h22],..., [p21], [p22],...],
                                            ...])

Each hij refers to the index of the j-th hole in the i-th K-point, each pij to the index of the j-th excited electron in the i-th K-point. For example, if an excited state calculation is initialize by promoting an electron from the ground state HOMO to the ground state LUMO, one needs to specify the indices of the ground state HOMO (hole) and LUMO (excited electron) in the spin channel where the excitation is performed. All rotations involving these orbitals are frozen during the constrained optimization resulting in these orbitals remaining unaltered after the optimization. It is also possible to constrain selected orbital rotations without completely freezing the involved orbitals by specifying lists of two orbital indices instead of lists of single orbital indices. However, care has to be taken in that case since constraining a single orbital rotation may not fully prevent mixing between those two orbitals during the constrained optimization.

Charge transfer excited state of N-phenylpyrrole

In this example, a calculation of a charge transfer excited state of the N-phenylpyrrole molecule is carried out. After a ground state calculation, a single excitation is performed from the HOMO to the LUMO in one spin channel. No spin purification is used, meaning that only the mixed-spin open-shell determinant is optimized. If an unconstrained optimization is performed from this initial guess, the calculation collapses to a first-order saddle point with pronounced mixing between the HOMO and LUMO and a small dipole moment of -3.396 D, which is not consistent with the wanted charge transfer excited state. Variational collapse is avoided here by performing first a constrained optimization freezing the hole and excited electron of the initial guess. Then the new orbitals are used as the initial guess of an unconstrained optimization, which converges to a higher-energy saddle point with a large dipole moment of -10.227 D consistent with the target charge transfer state.

from ase.io import read
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

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

atoms = read('N-Phenylpyrrole.xyz')
atoms.center(vacuum=5.0)
atoms.set_pbc(False)
atoms.calc = calc

# Ground state calculation
E_GS = atoms.get_potential_energy()

# Ground state LCAO coefficients and occupation numbers
C_GS = [calc.wfs.kpt_u[x].C_nM.copy() for x in range(len(calc.wfs.kpt_u))]
f_gs = [calc.wfs.kpt_u[x].f_n.copy() for x in range(len(calc.wfs.kpt_u))]

# Direct approach using ground state orbitals with changed occupation numbers
calc.set(eigensolver=LCAOETDM(searchdir_algo={'name': 'l-sr1p'},
                              linesearch_algo={'name': 'max-step'},
                              need_init_orbs=False),
         txt='N-Phenylpyrrole_EX_direct.txt')

# Spin-mixed open-shell occupation numbers
f = excite(calc, 0, 0, spin=(0, 0))

# Direct optimization maximum overlap method calculation
prepare_mom_calculation(calc, atoms, f)
E_EX_direct = atoms.get_potential_energy()

# Reset LCAO coefficients and occupation numbers to ground state solution
for k, kpt in enumerate(calc.wfs.kpt_u):
    kpt.C_nM = C_GS[k]
    kpt.f_n = f_gs[k]

h = 26  # Hole
p = 27  # Excited electron

# Constrained optimization freezing hole and excited electron
calc.set(eigensolver=LCAOETDM(constraints=[[[h], [p]], []],
                              need_init_orbs=False),
         txt='N-Phenylpyrrole_EX_constrained.txt')

# Spin-mixed open-shell occupation numbers
f = excite(calc, 0, 0, spin=(0, 0))

# Direct optimization maximum overlap method calculation
prepare_mom_calculation(calc, atoms, f)
E_EX_constrained = atoms.get_potential_energy()

# Unconstrained optimization using constrained solution as initial guess
calc.set(eigensolver=LCAOETDM(searchdir_algo={'name': 'l-sr1p'},
                              linesearch_algo={'name': 'max-step'},
                              need_init_orbs=False),
         txt='N-Phenylpyrrole_EX_from_constrained.txt')

# Direct optimization maximum overlap method calculation
E_EX_from_constrained = atoms.get_potential_energy()

Excited state of silicon using plane wave

In this example, the singlet excited state of a conventional silicon is calculated using plane waves approach.

from gpaw import GPAW
from ase.build import bulk
from gpaw.directmin.etdm_fdpw import FDPWETDM
from gpaw.mom import prepare_mom_calculation
from gpaw.directmin.tools import excite
from ase.parallel import paropen


atoms = bulk('Si', 'diamond', a=5.44, cubic=True)
atoms.cell[0, 0] += 1e-3

# Step: Set up the GPAW calculator
calc = GPAW(mode={'name': 'pw',    # Use plane wave mode
                  'ecut': 340},   # Cutoff energy
            xc='PBE',
            kpts=(1, 1, 1),
            eigensolver=FDPWETDM(converge_unocc=True),
            mixer={'backend': 'no-mixing'},
            occupations={'name': 'fixed-uniform'},
            spinpol=True,
            )

atoms.calc = calc
E_gs = atoms.get_potential_energy()


calc.set(eigensolver=FDPWETDM(excited_state=True,
                              converge_unocc=False,
                              momevery=10,
                              max_step_inner_loop=0.2,
                              maxiter_inner_loop=20))


f_sn = excite(calc, 0, 0, (0, 0))
prepare_mom_calculation(calc, atoms, f_sn)
E_es = atoms.get_potential_energy()

print(f'Excitation energy: {E_es - E_gs}')

with paropen('si_excited.txt', 'w') as fd:
    print(f'Excitation energy Si: {E_es - E_gs} eV',
          file=fd)

References

[1]

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)

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]

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