Excited States from Restricted Open Shell Plane-Wave DFT

This paper presents a plane-wave implementation of restricted open-shell Kohn-Sham (ROKS) density functional theory within VASP that enables accurate, spin-pure excited-state calculations with analytical forces for extended systems, demonstrating performance comparable to time-dependent DFT while retaining the favorable scaling of ground-state methods.

The Gist

Imagine you are trying to predict how a material will behave when it gets a "sugar rush" of energy—like when sunlight hits a solar cell or when an LED lights up. In the world of physics, this is called an excited state.

For a long time, scientists had a choice: use a method that was cheap and fast but often got the details wrong (like a blurry photo), or use a method that was incredibly accurate but so slow it could take years to run on a supercomputer for a single molecule.

This paper introduces a new way to get the best of both worlds. The authors have built a tool inside a famous software program called VASP that can calculate these "sugar rush" states quickly and accurately, even for huge materials like crystals.

Here is how they did it, explained through simple analogies:

1. The Problem: The "Spin" Confusion

Think of electrons in an atom like dancers on a floor.

• Ground State: The dancers are all paired up, holding hands, and spinning in perfect harmony. This is stable and easy to calculate.
• Excited State: One dancer jumps up and starts spinning wildly. Now, the group is unbalanced.

The old, fast methods tried to describe this wild dancer using a single, simple rule. But this caused a problem called "spin contamination." It's like trying to describe a chaotic dance party by pretending everyone is still holding hands in a neat circle. The math gets messy, and the prediction of how much energy the dancer needs to jump is often wrong.

2. The Solution: The "Restricted Open-Shell" (ROKS) Trick

To fix this, the authors used a clever trick called Restricted Open-Shell Kohn-Sham (ROKS).

Imagine you want to know the energy of that chaotic dance party. Instead of guessing, the authors say: "Let's look at two different versions of the party at the same time."

1. Version A: The wild dancer spins one way.
2. Version B: The wild dancer spins the opposite way.

They take the average of these two versions and mix it with a third version where the dancer spins in a specific "triplet" pattern. By mathematically blending these three scenarios, they cancel out the messy "spin contamination" errors. The result is a pure, clean picture of the excited state that is just as accurate as the slow, expensive methods but runs at the speed of the fast, cheap ones.

3. The Engine: Finding the Lowest Point

To find the right answer, the computer has to "climb down a hill" to find the lowest energy point (the most stable state).

• The Old Way: Sometimes the computer would slip and fall into the wrong valley (the ground state) instead of the excited state valley.
• The New Way: The authors built a special "preconditioned" engine. Think of this as giving the computer a pair of high-tech boots with springs. These boots help the computer feel the shape of the hill better, so it can slide down to the correct excited valley without slipping back to the ground. They used two different driving styles for this:
  • Conjugate Gradient (CG): A steady, efficient hiker who checks the path ahead.
  • DIIS: A smart navigator that remembers past steps to correct its course quickly.

4. The Proof: Testing the Tool

The team didn't just build the tool; they tested it rigorously.

• The Small Test: They ran the tool on eight small organic molecules (like ingredients in a perfume or plastic). They compared their results to a gold-standard chemistry program called Q-Chem. The results were nearly identical, with differences so small they were like measuring the width of a human hair against the distance from New York to London.
• The Big Test: They applied it to Magnesium Oxide (MgO), a solid crystal with a tiny hole (a vacancy) in it. This is a real-world material used in things like ceramics and electronics. They calculated how this crystal glows when excited.
  • They compared their results to a method called TDDFT (Time-Dependent DFT), which is the current industry standard for accuracy but is very slow.
  • The Result: Their new method gave answers very close to the slow standard (within about 0.2 electron-volts), but it kept the speed advantage of the fast method.

5. Why This Matters

The paper shows that you don't have to sacrifice speed for accuracy anymore.

• For Materials: Scientists can now study huge, complex materials (like defective crystals or surfaces) to see how they absorb light or store energy.
• For Forces: The tool doesn't just calculate energy; it also calculates forces. This is like knowing not just how high the dancer jumped, but also which way they pushed the floor. This allows scientists to simulate how the atoms move and relax after getting excited, which is crucial for designing better solar cells or light-emitting devices.

In summary: The authors have built a "fast lane" for calculating excited states. They fixed the math errors that used to plague fast calculations, allowing researchers to study complex, real-world materials with high accuracy without waiting years for a computer to finish the job.

Technical Summary

Technical Summary: Excited States from Restricted Open Shell Plane-Wave DFT

Problem Statement
Modeling excited-state properties is critical for designing optoelectronic materials, such as LEDs and solar energy converters. However, realistic simulations of defective bulk materials or surfaces often require hundreds or thousands of atoms, necessitating computationally inexpensive methods. While Variational Excited-State Density Functional Theory (DFT), specifically the $\Delta$-Self-Consistent Field ($\Delta$-SCF) approach, offers ground-state scaling costs, it suffers from two primary limitations:

1. Spin Contamination: Single-configuration approaches often fail to describe excited singlet states correctly because a single Slater determinant cannot yield the correct total spin expectation value ($\langle \hat{S}^2 \rangle = 0$). This typically leads to underestimated excitation energies.
2. Collapse: During optimization, the desired excited configuration may collapse to the ground state or a lower-lying excited state.
   Existing techniques to mitigate spin contamination, such as using unrestricted Kohn-Sham (UKS) calculations with spin-purification formulas (e.g., $E_S = 2E(\Phi_1) - E(\Phi_3)$), are approximate because they rely on different orbital sets for singlet and triplet configurations, leaving residual spin contamination.

Methodology
The authors implement Restricted Open-Shell Kohn-Sham (ROKS) DFT within the plane-wave projector augmented-wave (PAW) framework of the VASP code. The core methodology involves:

• Variational Minimization: Instead of solving standard Kohn-Sham equations, the method variationally minimizes a weighted energy functional combining mixed-spin and triplet configurations to recover spin-pure singlet energies. The energy functional is:
  $$E = \sum_L c_L E_L - \sum_{m,n} \gamma_{mn} (\langle m|S|n\rangle - \delta_{mn})$$
  where $L$ labels electronic configurations, $c_L$ are coefficients, and $\gamma_{mn}$ are Lagrange multipliers enforcing orbital orthonormality under the PAW overlap operator $S$.
• Optimization Algorithms:
  • Preconditioned Conjugate-Gradient (CG): The primary implementation uses a preconditioned Fletcher-Reeves CG algorithm. The energy gradient is decomposed into "out-of-subspace" and "in-subspace" components. The out-of-subspace component is preconditioned using a kinetic-energy-based operator, while the in-subspace component is optimized via unitary transformations (diagonalization of a generalized operator matrix $R$).
  • Direct Inversion in the Iterative Subspace (DIIS): A DIIS solver is also implemented to handle saddle points, allowing the modeling of excited states beyond the lowest singlet.
• Force Calculation: Analytical atomic forces are derived using a generalized Hellmann-Feynman theorem. The force expression is adapted to the PAW formalism, leveraging existing VASP routines for restricted and unrestricted DFT forces with minor modifications.
• Stabilization Strategies: To prevent collapse to triplet-like states (where the excited singlet becomes nearly degenerate with the triplet), the authors employ strategies such as suppressing coupling between open-shell orbitals in the CG algorithm or adjusting coupling coefficients ($\lambda_{ab}$) in the DIIS solver.

Key Contributions

1. Implementation: The first implementation of ROKS within a plane-wave PAW framework (VASP), enabling excited-state calculations for extended systems with hundreds of atoms.
2. Analytical Forces: Derivation and implementation of analytical forces compatible with the PAW formalism, enabling geometry relaxation and molecular dynamics in excited states.
3. Algorithmic Efficiency: Demonstration that the preconditioned CG approach converges significantly faster (fewer iterations) than DIIS for ground-state-like optimizations due to superior preconditioning of the gradient.

Results

• Molecular Benchmarking: The implementation was validated against the Q-Chem quantum chemistry code for eight organic molecules.
  • Accuracy: Mean absolute deviations (MAD) between VASP and Q-Chem were approximately 30 meV (0.030 eV) for LDA, PBE, and PBE0 functionals. These differences are attributed to basis set and pseudopotential approximations (plane-wave/PAW vs. Gaussian/all-electron).
  • Challenges: For molecules like cytosine, where the excited state shares symmetry with the ground state, avoiding collapse to a triplet-like state required combining CG and DIIS approaches or specific coupling suppression.
  • Forces: Analytical forces showed excellent agreement with finite-difference numerical forces (MAD < 0.3 meV/Å).
• Solid-State Application (MgO with Oxygen Vacancy): The method was applied to the $F_0$ center (neutral oxygen vacancy) in MgO, calculating the three lowest excitations ($^1T_{1u}$, $^3T_{1u}$, and Conduction Band Minimum).
  • Comparison with TDDFT: Vertical excitation energies from ROKS and Time-Dependent DFT (TDDFT) using a dielectric-dependent hybrid (DDH) functional differed by an average of 0.21 eV.
  • Structural Properties: Franck-Condon shifts and mass-weighted displacements showed average deviations of 0.14 eV and 0.12 amu$^{1/2}$Å, respectively, between ROKS and TDDFT.
  • Functional Dependence: A key finding is that ROKS properties (energies, forces, displacements) are less sensitive to the choice of DFT functional (PBE vs. DDH) compared to TDDFT. For instance, the change in vertical excitation energy between PBE and DDH was 0.93 eV for ROKS but 1.81 eV for TDDFT in the CBM excitation.
  • Triplet Collapse: For the CBM excitation, ROKS tended to collapse to a triplet-like state. However, the authors note that the singlet-triplet gap is small for this specific transition, suggesting the error in vertical excitation energy remains minor.

Significance
The paper claims that ROKS provides excitation energies and excited-state forces with accuracy comparable to TDDFT while retaining the favorable computational scaling of ground-state DFT. This makes the approach a promising tool for affordable excited-state simulations in extended systems, particularly where large supercells are required to model defects or surfaces. The reduced sensitivity to the exchange-correlation functional in ROKS compared to TDDFT further suggests its robustness for predicting structural relaxations and potential energy surfaces in solid-state materials.