Real-Space Hubbard-Corrected DFT+U

• Real-space Hubbard-corrected DFT+U is a framework that augments standard DFT with explicit on-site (and inter-site) corrections to accurately model strong electron correlations.
• It utilizes high-order finite-difference discretization and parallel implementation strategies to efficiently compute energy, forces, and stress in complex materials.
• The methodology enables optimal Hubbard parameter tuning via hybrid functional matching, ensuring precise phase stability and electronic property predictions for correlated systems.

Real-space Hubbard-corrected density functional theory (DFT+U in real space) is a formalism and computational methodology that augments conventional Kohn–Sham DFT with explicit on-site (and, in some cases, inter-site) corrections to address strong electron correlation and self-interaction errors, especially in systems with localized d and f electrons. The "real-space" aspect emphasizes the finite-difference spatial discretization and non-basis-set-specific implementation, enabling large-scale parallelism and improved efficiency for inhomogeneous, complex materials. This framework provides transparent expressions for the Hubbard-corrected energy, forces, and stress tensor, suited for high-throughput and large-scale electronic structure calculations.

1. Formalism and Derivation of Energy, Forces, and Stress in Real-space DFT+U

In real-space DFT+U, the total energy is

$E_\text{DFT+U} = E_\text{DFT} + E_U$

where the Hubbard correction for each correlated atom $I$, angular momentum $\ell$, and spin $\sigma$ is given in the rotationally invariant Dudarev formalism: $$E_U = \sum_{\sigma, I, \ell} \frac{\tilde{U}^\sigma_{I\ell}}{2} \mathrm{Tr} \bigl\{ \mathcal{D}^{\sigma}_{I\ell} - (\mathcal{D}^{\sigma}_{I\ell})^2 \bigr\}$$ The occupation matrix is constructed from Bloch-periodically mapped, atom-centered orbitals $\phi_{I\ell m}$ and the Kohn–Sham orbitals $\psi_j^\sigma$: $$(\mathcal{D}^{\sigma}_{I\ell})_{mm'} = \sum_j f_j^\sigma \left( \int_{\Omega} \psi_j^{\sigma*}(\mathbf{r}) \phi_{I\ell m'}(\mathbf{r} - \mathbf{R}_I) \, d\mathbf{r} \right) \left( \int_{\Omega} \phi_{I\ell m}^*(\mathbf{r} - \mathbf{R}_I) \psi_j^\sigma(\mathbf{r}) \, d\mathbf{r} \right)$$ This real-space discretization is compatible with high-order finite-difference integration and differentiation, naturally fitting grid-based large-scale implementations (Bhowmik et al., 31 Jul 2025).

The variational derivatives yield key ingredients for efficient simulations:

• Hubbard potential: acting on each orbital,

$$V_U^\sigma \psi_n^\sigma(\mathbf{r}) = \sum_{I, \ell, m, m'} \frac{\tilde{U}^{\sigma}_{I\ell}}{2} \left[ \delta_{mm'} - 2 (\mathcal{D}_{I\ell}^{\sigma})_{mm'} \right] \phi_{I\ell m'}(\mathbf{r}) \int_\Omega \phi_{I\ell m}^*(\mathbf{r}') \psi_n^\sigma(\mathbf{r}') d\mathbf{r}'$$

• Forces: (HeLLMann–Feynman term plus Pulay corrections, all on real-space grid)

$$\mathbf{f}^J_U = -\frac{\partial E_U}{\partial \mathbf{R}_J}$$

giving explicit dependence on derivatives of orbitals and occupation matrices with respect to atomic positions.

• Stress tensor: components are given by

$$(\sigma_U)_{\alpha \beta} = \frac{1}{|\Omega|} \frac{\partial E_U}{\partial F_{\alpha \beta}} \Big|_{\mathcal{G}}$$

where $F_{\alpha \beta}$ is the deformation gradient and $\mathcal{G}$ the ground-state configuration. All these expressions are explicitly compatible with finite-difference strategies.

2. Real-space Finite-Difference Discretization and Parallel Implementation

The formalism is discretized on a uniform real-space grid using high-order centered finite-difference approximations for derivatives and the trapezoidal rule for integrals:

• All essential integrals—for density matrix, Hubbard energy, and their derivatives—are replaced with quadratures over the grid.
• Derivatives and gradients needed for forces and stress are computed numerically.
• The atomic orbitals $\phi_{I\ell m}(\mathbf{r})$ are represented on the grid and often have the same spatial extent and behavior as nonlocal pseudopotential projectors, allowing code reuse and efficient evaluation.

Parallel implementation in state-of-the-art codes is based on four layers:

1. Distribution over spin channels,
2. Distribution over Bloch k-points,
3. Distribution over electron bands,
4. Decomposition of the spatial domain (slab or cube partitioning).

This design enables scalability for both large system sizes and processor counts, demonstrating strong-scaling efficiency exceeding 87% for TiO₂ and yielding speedups of $8\times$–$47\times$ compared to established planewave codes like Quantum ESPRESSO for the same benchmarks (Bhowmik et al., 31 Jul 2025).

3. Exchange–Correlation Consistency in Orbital Projection and Impact on Predictions

Constructing atomic-like orbitals for projection (and thus $\mathcal{D}_{I \ell}^\sigma$) must be done self-consistently with the exchange–correlation (XC) functional used in the DFT calculation:

• XC inconsistency arises if, for example, atom-centered orbitals from LDA are used in an r²SCAN (meta-GGA) calculation.
• Absolute energies may differ by up to $0.0046$ Hartree per TiO₂ formula unit for rutile when using PBE vs. r²SCAN orbitals, but energy differences (e.g., between polymorphs) are preserved, and correct relative phase stability ordering is maintained for Hubbard $U$ in the tested range. This suggests error cancellation is robust for energy differences but not for total absolute energies.
• The conclusion is that, while absolute DFT+U energies with mismatched XC in orbital construction should not be compared, phase stabilities and similar energy differences can be compared safely if the same projection orbitals are used across all systems (Bhowmik et al., 31 Jul 2025).

4. Optimizing the Hubbard Parameter via Hybrid Functional Matching

Optimal choice of the Hubbard parameter $\tilde{U}$ is critical for quantitative predictions in correlated materials:

• The approach introduced matches the r²SCAN+U ground-state density to the energy obtained from a hybrid functional (e.g., HSE) by post-processing: for each candidate $U$, the non-self-consistent hybrid energy is evaluated on the r²SCAN+U density, and the minimum with respect to $U$ determines $U_\text{opt}$.
• For TiO₂ polymorphs, $U_\text{opt}$ values of 5.66–6.06 eV yield band gaps in rutile, anatase, and brookite in good agreement with experiments (2.82, 3.26, 3.29 eV, respectively).
• Because hybrid functionals are considerably more expensive than meta-GGA or GGA+U, this protocol enables systematic calibration of $U$ with minimal computational effort, refining predictive accuracy for gap and phase stability (Bhowmik et al., 31 Jul 2025).

5. Benchmarking and Scalability: Performance Versus Planewave Codes

The implementation was benchmarked against planewave reference codes using TiO₂ polymorphs:

• Real-space DFT+U with finite-difference grids produces total energies, forces, and phase stabilities matching planewave results to high precision.
• Time-to-solution is significantly reduced. As system size and processor count increase, the real-space approach demonstrates better parallel efficiency and scalability. Minimum times to solution improve by more than an order of magnitude, especially for large-scale problems, with planewave codes failing to execute for some large cases due to communication and memory bottlenecks inherent to the reciprocal-space formalism (Bhowmik et al., 31 Jul 2025).

6. Implications for Materials Simulation and Prospects

The methodology enables large-scale, high-throughput, and high-accuracy calculations for systems with strong electron correlation and inhomogeneous environments, including surfaces, interfaces, and disordered materials. Real-space DFT+U, being non-basis-set-specific and highly parallelizable, is highly suited for:

• Studies requiring flexible boundary conditions or localized basis corrections,
• Simulations of extended, defective, or highly anisotropic crystalline and non-crystalline systems,
• Development and testing of new correlation functionals and nonlocal corrections.

A plausible implication is that real-space Hubbard correction frameworks will become essential for next-generation electronic structure codes targeting exascale supercomputing and materials-by-design applications, especially for strongly correlated and functional oxide systems.

7. Summary Table: Core Formulas in Real-space DFT+U

Quantity	Formula (Real-space Discretization)	Key Role
Energy correction $E_U$	$$E_U = \sum_{\sigma, I, \ell} \frac{\tilde{U}^\sigma_{I\ell}}{2} \mathrm{Tr} \left\{ \mathcal{D}_{I\ell}^\sigma - (\mathcal{D}_{I\ell}^\sigma)^2 \right\}$$	Augments DFT for correlation
Potential $V_U^\sigma$	See Section 2.A above	Adds orbital-dependent correction to Hamiltonian
Forces $f_U^J$	See Section 2.B above	Needed for ionic relaxation and MD
Stress $(\sigma_U)_{\alpha \beta}$	See Section 2.C above	Essential for accurate equation of state

All quantities and operations are evaluated on a finite-difference grid, ensuring fidelity, locality, and parallel efficiency.

---

Real-space Hubbard-corrected DFT thus provides a rigorous, highly efficient, and scalable framework that accurately incorporates strong correlation effects within a grid-based electronic structure formalism, supporting advanced simulation and design of complex correlated materials (Bhowmik et al., 31 Jul 2025).