Development of ab initio Hubbard parameter calculation schemes in the k-point sampling real-time TDDFT program in CP2K

This paper presents the implementation of ab initio Hubbard parameter calculation schemes, including a novel linear-response method for energy-dependent parameters, within CP2K's k-point sampling real-time TDDFT program, highlighting the distinct theoretical advantages and dynamical applications of this approach compared to the ACBN0 scheme.

The Gist

The Big Picture: Fixing the "Self-Interference" Problem

Imagine you are trying to predict how a crowd of people behaves in a room. If the people are just walking around casually, standard rules of physics (called DFT) work great. But what if the room is packed with a group of very intense, emotional people who constantly bump into each other and react strongly to every move? In the world of atoms, these are electrons in materials like transition metal oxides (think rust, catalysts, or battery materials).

Standard physics rules fail here because they make a mistake: they let an electron "talk to itself." It's like a person in a crowd shouting at their own echo and getting confused about how loud the room actually is. This is called the Self-Interaction Error.

To fix this, scientists invented a patch called DFT+U. Think of U as a "personal space bubble" parameter. It tells the computer, "Hey, these specific electrons are grumpy and need extra space; don't let them crowd each other too much."

The Problem: The tricky part is figuring out exactly how big that personal space bubble (the Hubbard U) should be.

• Old way: Scientists would guess a number, run the simulation, and tweak it until the result matched an experiment. It was like tuning a radio by ear until the static cleared up.
• New way (This Paper): The authors built a machine that calculates the perfect size of that bubble from first principles, without guessing.

The Two New Tools

The authors implemented two different "machines" inside a powerful software called CP2K to calculate this U value automatically.

1. The "Instant Snapshot" Machine (ACBN0)

• How it works: Imagine taking a high-speed photo of the electron crowd. You look at the photo, count how many people are in a specific spot, and instantly calculate how much space they need based only on that single snapshot.
• The Analogy: It's like a traffic cop who looks at a frozen frame of a traffic jam and immediately assigns a speed limit. It's fast and works great for movies (simulations that change over time, like a laser hitting a material).
• The Catch: The math behind this "snapshot" isn't derived from the deepest laws of physics. It works well, but it's a bit of a "black box." We know it works, but it's hard to predict exactly why it behaves the way it does in complex situations.

2. The "Echo Chamber" Machine (Minimum-Tracking Linear Response)

• How it works: This is more like a sound engineer. You tap a specific electron (the "perturbation") and listen to how the rest of the crowd reacts (the "response"). By measuring how the crowd "screens" or dampens that tap, you can calculate the exact personal space needed.
• The Analogy: Imagine shouting in a canyon. The echo tells you about the shape of the canyon. This method shouts at the electrons and listens to the echo to figure out the rules.
• The Innovation: The authors took this method and made it energy-dependent.
  • Static version: "What is the personal space needed for a calm electron?"
  • Dynamic version: "What is the personal space needed for an electron that is vibrating at 100 Hz? What about 1,000 Hz?"
  • Why this matters: Electrons change their behavior depending on how much energy they have. This new tool allows the computer to say, "At low energy, the bubble is size X. At high energy, the bubble shrinks to size Y."

The Showdown: Which is Better?

The authors tested both machines on various materials (like Nickel Oxide and Manganese Oxide).

• For Static Results (The Photo): Both machines did a decent job. They both predicted the "band gap" (the energy needed to turn the material into a conductor) reasonably well compared to real-world experiments. However, they gave different numbers for the U value itself. It's like two chefs making a cake that tastes the same, but using different amounts of sugar.
• For Dynamic Results (The Movie):
  • The Instant Snapshot (ACBN0) is great for simulating what happens when a laser hits a material. It's fast and handles time-varying situations well.
  • The Echo Chamber (Linear Response) is more theoretically rigorous. Because it's based on deep physics, it can tell you how the "personal space" changes as the energy of the electrons changes. This is crucial for understanding high-energy physics, but it is computationally expensive (it takes a long time to run).

The "Energy-Dependent" Breakthrough

The most exciting part of this paper is the new Energy-Dependent calculation.

Think of the Hubbard U not as a fixed number, but as a chameleon.

• In the past, we treated the chameleon as a static color (e.g., "It's always green").
• This new method realizes the chameleon changes color based on the lighting (energy).
• The authors showed that by using the "Echo Chamber" method, they can calculate exactly how the "color" (the U value) shifts as the energy changes. This is a huge step forward because it allows for much more accurate simulations of materials under extreme conditions, like those hit by powerful lasers.

The Bottom Line

The authors have successfully built two new tools inside the CP2K software:

1. ACBN0: A fast, practical tool for simulating how materials react to light and time.
2. Linear Response: A rigorous, deep-physics tool that can calculate how electron interactions change with energy.

They proved that while both tools are good at predicting static properties, they serve different masters. If you need to simulate a movie of a material reacting to a laser, use the fast snapshot tool. If you need to understand the deep, energy-dependent physics of how electrons interact, use the rigorous echo-chamber tool.

This work bridges the gap between simple, fast calculations and complex, high-precision physics, giving scientists better ways to design new materials for batteries, solar cells, and quantum computers.

Technical Summary

1. Problem Statement

Density Functional Theory (DFT) with local or semilocal exchange-correlation (xc) functionals fails to accurately describe materials with strong electron-electron correlations, such as transition metal (TM) oxides, due to self-interaction errors. The DFT+U method addresses this by introducing explicit on-site interaction parameters (Hubbard $U$ and Hund's $J$).

However, a critical challenge remains:

• Parameter Determination: Empirical fitting of $U$ and $J$ is inconsistent and environment-dependent. Ab initio methods are required.
• Static vs. Dynamic: Most existing ab initio schemes calculate static (energy-independent) parameters. However, in non-equilibrium dynamics (e.g., laser-induced processes) or for excited states far from the Fermi level, these parameters should be energy-dependent ($U(\omega)$) or time-dependent ($U(t)$) to reflect dynamic screening.
• Limitations of Current Methods:
  • Constrained RPA (cRPA): Can calculate energy-dependent parameters but often overestimates screening due to missing higher-order diagrams and struggles with strong orbital hybridization.
  • ACBN0 (Mean-field-like): Applicable to real-time dynamics but lacks a rigorous derivation from many-body theory, making the physical interpretation of its time-dependent behavior difficult.
  • Linear Response (Static): Historically limited to static calculations and often restricted to $\Gamma$-point sampling, limiting applicability to periodic solids.

2. Methodology

The authors implemented and extended two distinct ab initio schemes within the CP2K software suite, specifically integrating them into a k-point sampling Real-Time Time-Dependent DFT (RT-TDDFT) framework.

A. Implementation of ACBN0

• Approach: A mean-field-like approach where Hubbard parameters are derived from the static correlations embedded in the density matrix and Kohn-Sham orbitals.
• Formulation: The authors adapted the ACBN0 formalism to their k-point sampling projection scheme. They defined a "renormalized" occupation number based on the projected density matrix of localized orbitals (atomic orbitals of isolated atoms).
• Key Feature: $U$ and $J$ are calculated algebraically from the density matrix, allowing for immediate application to time-dependent simulations without solving complex response equations at every time step.

B. Implementation of Minimum-Tracking Linear Response (Static)

• Approach: Extends the method proposed by Moynihan et al. to k-point sampling.
• Formulation: A small perturbation is applied to a specific set of correlated orbitals. The Hubbard parameters are derived from the linear response of the Hartree-xc potential ($V_{Hxc}$) to changes in orbital occupation.
• Advantage: Unlike the standard Cococcioni linear response (which requires perturbations on all atoms), this "minimum-tracking" method calculates parameters based on the response of the perturbed site, reducing computational overhead in specific contexts, though it still requires supercell calculations for convergence.

C. New Proposal: Energy-Dependent Linear Response

• Innovation: The authors extended the minimum-tracking linear-response theory to the frequency domain ($\omega$) to calculate energy-dependent Hubbard parameters $U(\omega)$ and $J(\omega)$.
• Derivation:
  • They reformulated the static linear response expression to include the time-dependent Hartree-xc kernel ($f_{Hxc}$) and the density response function ($\chi$).
  • By applying an impulsive perturbation ($\delta(t)$) in RT-TDDFT and collecting the time-series response of occupation numbers and potentials, they Fourier-transformed the data to obtain the frequency-dependent response kernels.
  • Formula: $U(\omega)$ is derived from the ratio of the Fourier-transformed potential response to the occupation response, incorporating the retarded screening function $\epsilon^{-1(R)} = 1 + \chi^{(R)} f_{Hxc}$.
• Significance: This approach naturally includes the xc-functional dependence of the screening, unlike cRPA which relies solely on Coulombic screening.

3. Key Results

Static Calculations

• Materials Tested: ScN, Rutile TiO$_2$, MnO, NiO, Wurtzite ZnO, and CuO.
• ACBN0 Performance:
  • Reproduced band gaps and magnetic moments reasonably well compared to experiments (e.g., TiO$_2$ gap ~3.1 eV vs exp 3.3 eV).
  • Calculated $U_{eff}$ values showed discrepancies with previous literature (Agapito et al., Tancogne-Dejean et al.), attributed to differences in pseudopotentials and basis sets, but physical properties remained consistent.
• Linear Response Performance:
  • Successfully converged $U$ and $J$ values for the same materials.
  • Discrepancy: The $U$ values obtained via linear response differed significantly from ACBN0 (e.g., for NiO, ACBN0 $U \approx 3.6$ eV vs Linear Response $U \approx 5.4$ eV).
  • Accuracy: Neither method was universally superior; ACBN0 better reproduced the ScN band gap, while Linear Response better reproduced MnO.
  • Cost: Linear response was computationally much more expensive due to the need for large supercells and iterative self-consistency.

Dynamical Calculations

• ACBN0 in Strong Fields:
  • Simulated NiO under a strong laser pulse (760 nm, 1.1 TW/cm$^2$).
  • Result: Observed a time-dependent reduction in $U_{eff}$ (by ~140 meV) during and after the pulse, attributed to electronic excitation and enhanced screening. This aligns with previous experimental and theoretical findings, validating ACBN0 for non-equilibrium dynamics.
• Energy-Dependent Linear Response:
  • Calculated $U(\omega)$ and $J(\omega)$ for NiO and MnO.
  • Consistency: At $\omega = 0$, the dynamic results matched the static linear response values, confirming the theoretical extension.
  • Asymptotic Behavior: For NiO and O 2p orbitals, $U(\omega)$ correctly asymptoted to the unscreened Coulomb interaction at high frequencies ($\hbar\omega \approx 200$ eV).
  • Issue: For Mn 3d orbitals, the expected asymptotic behavior was not achieved, likely due to numerical instabilities where both numerator and denominator in the response ratio approach zero at high frequencies.

4. Significance and Contributions

1. First k-point RT-TDDFT Implementation: The work successfully integrates ab initio Hubbard parameter calculations into a k-point sampling RT-TDDFT framework in CP2K, enabling the study of periodic correlated materials under dynamic conditions.
2. New Energy-Dependent Scheme: The proposal of a linear-response-based method for calculating $U(\omega)$ provides a rigorous alternative to cRPA. It explicitly includes xc-functional effects and avoids the strict separation of "correlated" vs. "screening" orbitals, which is beneficial for materials with strong orbital hybridization.
3. Comparative Analysis: The paper provides a critical comparison between the mean-field-like ACBN0 and the rigorous linear-response methods. It highlights that while ACBN0 is computationally efficient and effective for dynamical simulations, its theoretical foundation is less transparent. Conversely, the linear-response method is theoretically robust but computationally demanding.
4. Foundation for Future DFT+U(ω): The developed framework lays the groundwork for advanced DFT+U(ω) theories (e.g., by Vanzini and Marzari), where energy-dependent parameters can be used to improve the description of low-energy effective models and high-energy excitations.

5. Conclusion

The authors have successfully implemented and tested two distinct ab initio schemes for Hubbard parameters in CP2K. While ACBN0 proves to be a practical tool for real-time dynamics, the newly proposed energy-dependent linear-response method offers a theoretically rigorous path forward for calculating frequency-dependent interactions. The study underscores the trade-off between computational cost and theoretical rigor, suggesting that the choice of method depends on the specific application (static properties vs. non-equilibrium dynamics) and the material system. Future work will focus on resolving numerical instabilities in the high-frequency limit of the linear-response scheme.