Feature Request: `gga_grad` Parameter — Three GGA Gradient Algorithms for Noncollinear Magnetism

[dyzheng]

Background

Summary

In NSPIN=4 (noncollinear magnetic) GGA calculations, the 4-component density (ρ, m_x, m_y, m_z) must be mapped to a 2-component (ρ↑, ρ↓) representation with their gradients before calling the exchange-correlation functional. The current develop branch only supports gga_grad=1 and gga_grad=2 modes. This issue describes the complete functional requirements and algorithm specifications for adding gga_grad=3 (Scalmani-Frisch transform) along with the associated dispatch logic and parameter validation.

---

1. Functional Requirements

1.1 Parameter Definition

Add/modify the gga_grad parameter in input_parameter.h:

int gga_grad = 1;  // 1: gradient of |m|; 2: gradient of m * hat{m}; 3: Scalmani-Frisch transform

Extend parameter validation in read_input_item_general.cpp to allow gga_grad = 1 | 2 | 3:

item.annotation = "1: gradient of |m|; 2: gradient of m * hat{m}; 3: Scalmani-Frisch transform";
// check_value: gga_grad should be 1, 2, or 3

The parameter only takes effect when NSPIN=4 and DOMAG || DOMAG_Z is true; it should be ignored or raise an error otherwise.

---

2. Algorithm Specifications

2.1 Common Preprocessing

All three modes share the same density construction:

ρ↑  = (ρ + |m|) / 2
ρ↓  = (ρ - |m|) / 2
|m| = √(m_x² + m_y² + m_z²)
m̂_μ = m_μ / |m|          (when |m| > threshold; otherwise m̂_μ = 0)

2.2 gga_grad = 1: Gradient of |m|

Algorithm:

1. Compute |m| on real-space grid points
2. Construct ρ↑ = (ρ + ρ_core/2 + |m|/2), ρ↓ = (ρ + ρ_core/2 − |m|/2)
3. Transform ρ↑, ρ↓ to reciprocal space via FFT, then take gradients:
   • ∇ρ↑ = FFT_gradient(ρ↑)
   • ∇ρ↓ = FFT_gradient(ρ↓)
4. Call GGA functional to obtain v↑, v↓, h↑, h↓
5. Potential assembly:
   • v_ρ = ½(v↑ + v↓)
   • v_μ = ½(v↑ − v↓) · m̂_μ
6. div(h) correction:
   • v_ρ −= div(½(h↑ + h↓))
   • v_μ −= div(½(h↑ − v↓)) · m̂_μ (with double projection: div is computed per-μ, then multiplied by m̂_μ)

Known issue: |m| is non-differentiable where the magnetization direction changes abruptly (cusp), causing numerical discontinuities.

Code location: xc_functional_gradcorr.cpp, the else branch in the NSPIN=4 section (when gga_grad != 2).

---

2.3 gga_grad = 2: Gradient of m · m̂

Algorithm:

1. Density construction ρ↑, ρ↓ same as above

2. Gradient computation: instead of taking the gradient of |m| directly, use the projection:

   • ∇ρ↑ = ½∇ρ_total + ½ Σ_μ (m̂_μ · ∇m_μ)
   • ∇ρ↓ = ½∇ρ_total − ½ Σ_μ (m̂_μ · ∇m_μ)

   where ∇ρ_total = FFT_gradient(ρ[0] + ρ_core), ∇m_μ = FFT_gradient(ρ[μ])

3. Call GGA functional to obtain v↑, v↓, h↑, h↓

4. Potential assembly and div(h) correction similar to gga_grad=1, but the magnetic component of div(h) is handled differently:

   • v_ρ += ½(v↑ + v↓) − div(½(h↑ + h↓))
   • v_μ += ½(v↑ − v↓) · m̂_μ − div(½(h↑ − h↓) · m̂_μ) ← single projection onto each μ direction

Advantage: Avoids the non-differentiability of |m|; each magnetic component independently undergoes div(h) projection.

Code locations:

• Gradient computation: xc_functional_gradcorr.cpp, if (XC_Functional::gga_grad == 2) branch
• Potential assembly: xc_functional_gradcorr.cpp:588-675
• Libxc path: xc_functional_vxc.cpp:274-334 (gga_grad=2 branch inside v_xc_libxc)

---

2.4 gga_grad = 3: Scalmani-Frisch Transform

Reference: G. Scalmani, M. J. Frisch, J. Chem. Theory Comput. 2012, 8, 2193-2196.

Core idea: Instead of naively splitting the 4-component problem into two 2-component problems, derive the complete GGA potential expression in the (n, m_x, m_y, m_z) basis.

2.4.1 SF Transform for σ (Libxc sigma)

Gradient variable transformation (SF paper Eqs. 3-5):

n_±  = (n ± |m|) / 2                                    (Eq. 3)

A = ∇n · ∇n
B = Σ_μ ∇m_μ · ∇m_μ
t_μ = ∇n · ∇m_μ
C = √(t_x² + t_y² + t_z²)
f_grad = sgn(t_x·m_x + t_y·m_y + t_z·m_z)

γ_a = (A + B)/4 + (f_grad/2)·C     →  σ_uu = γ_a    (Eq. 4 upper)
γ_b = (A + B)/4 − (f_grad/2)·C     →  σ_dd = γ_b    (Eq. 4 lower)
γ_c = (A − B)/4                      →  σ_ud = γ_c    (Eq. 5)

2.4.2 Potential Formula (SF Paper Eqs. 9-10)

Map the vsigma returned by Libxc to SF coefficients:

V_a = vsigma[0]    (corresponding to ρ↑ρ↑)
V_c = vsigma[1]    (corresponding to ρ↑ρ↓)
V_b = vsigma[2]    (corresponding to ρ↓ρ↓)

pref_n_diag = (V_a + V_b + V_c) / 2
pref_m_diag = (V_a + V_b − V_c) / 2
cross       = f_grad · (V_a − V_b) / (2·C)     (when C > threshold; otherwise cross = 0)

The h-vectors for potential:

h_n   = e² · [pref_n_diag · ∇n + cross · (t_x·∇m_x + t_y·∇m_y + t_z·∇m_z)]
h_m_μ = e² · [pref_m_diag · ∇m_μ + cross · t_μ · ∇n]

v(comp, ir) -= div(h_comp)

2.4.3 Current Simplified Implementation (builtin PBE path)

The current gga_grad=3 builtin implementation uses the same gradient computation as gga_grad=2 (m̂_μ·∇m_μ projection), not the full SF σ transform. The div(h) correction removes the double m̂_μ projection from gga_grad=2:

v(0) += ½(v↑_gga + v↓_gga) − div(½(h1+h2))
v(μ) += ½(v↑_gga − v↓_gga)·m̂_μ − div(½(h1−h2)·m̂_μ)    ← independent div per direction

The difference from gga_grad=2 is that gga_grad=2 performs a double projection Σ_m̂ v(···)·m̂_μ, while gga_grad=3 only projects once onto each μ direction.

2.4.4 Dispatch Logic

Potential calculation (v_xc path):

• When gga_grad==3 && NSPIN==4 && (DOMAG || DOMAG_Z), directly call NCGGA_SF_Builtin::v_xc_ncgga_sf_builtin() and return early
• Does not go through the original LDA loop + gradcorr() path

Stress calculation (gradcorr path):

• When gga_grad==3 && NSPIN==4 && (DOMAG || DOMAG_Z), call NCGGA_SF_Builtin::gradcorr_ncgga_sf_builtin() and return early

2.4.5 Libxc Path

The full SF transform requires Libxc path support (when use_libxc=true). A new file xc_functional_gga_noncol_sf.cpp/h is needed, using Libxc's xc_gga_exc_vxc instead of the builtin PBE formulas, with the same SF transform logic otherwise.

Current status: Libxc path is not yet implemented.

---

3. Development Tasks for the develop Branch

3.1 Input Parameter Extension

File	Change
source/module_parameter/input_parameter.h	Add option 3 to gga_grad annotation
source/module_io/read_input_item_general.cpp	Allow gga_grad=3; update error message to "should be 1, 2, or 3"

3.2 New Files

File	Purpose
source/module_hamilt_general/module_xc/xc_functional_gga_noncol_sf_builtin.cpp	SF potential and stress implementation using builtin PBE
source/module_hamilt_general/module_xc/xc_functional_gga_noncol_sf_builtin.h	Header declaring functions in the NCGGA_SF_Builtin namespace
source/module_hamilt_general/module_xc/xc_functional_gga_noncol_sf.cpp	(future, optional) SF implementation using Libxc
source/module_hamilt_general/module_xc/xc_functional_gga_noncol_sf.h	(future, optional) Libxc SF header

3.3 Modified Files — Dispatch Logic

File	Change
source/module_hamilt_general/module_xc/xc_functional_vxc.cpp	Add gga_grad==3 check before LDA/v_xc computation; call SF builtin and return
source/module_hamilt_general/module_xc/xc_functional_gradcorr.cpp	Add gga_grad==3 check in the stress branch; call SF builtin stress function and return
source/module_hamilt_general/module_xc/CMakeLists.txt	Add xc_functional_gga_noncol_sf_builtin.cpp to the build

3.4 Core Algorithm Implementation Specification

3.4.1 v_xc_ncgga_sf_builtin Function

Inputs: nrxx, omega, tpiba, chr (pointer to Charge object)

Steps:

1. Precondition check: NSPIN==4 && (DOMAG || DOMAG_Z), throw exception otherwise
2. Density construction: Compute ρ↑, ρ↓, |m|, m̂_μ
3. Gradient computation: Compute ∇ρ↑, ∇ρ↓ using gga_grad=2 method (m̂_μ·∇m_μ projection)
4. LDA potential: Pointwise loop, call xc_spin for LDA part (e_xc, v↑_lda, v↓_lda), assemble 4-component LDA potential v(0..3)
5. GGA gradient correction (only when is_gga == true):
   • Pointwise loop, call gcx_spin and gcc_spin to obtain v1xup, v1xdw, v2xup, v2xdw, v2c, etc.
   • Assemble v↑_gga, v↓_gga, h1, h2
   • Density-derivative part: v(0) += ½(v↑_gga + v↓_gga), v(μ) += ½(v↑_gga − v↓_gga)·m̂_μ
   • div(h) correction:
     • v(0) -= div(½(h1+h2))
     • v(μ) -= div(½(h1−h2)·m̂_μ) ← independent div per μ direction
6. Statistics: Accumulate etxc, vtxc; MPI reduce then multiply by omega/nxyz

3.4.2 gradcorr_ncgga_sf_builtin Function

Inputs: chr, rhopw, ucell; Output: stress_gga

Steps:

1. Density and gradient construction same as v_xc_ncgga_sf_builtin
2. Pointwise loop, call gcx_spin / gcc_spin to obtain v2xup, v2xdw, v2c, etc.
3. Stress tensor accumulation:

   σ_lm += e²(∇ρ↑_l·∇ρ↑_m·v2xup + ∇ρ↓_l·∇ρ↓_m·v2xdw)
   σ_lm += e²(∇ρ↑_l·∇ρ↑_m·v2cup + ∇ρ↓_l·∇ρ↓_m·v2cdw + (∇ρ↑_l·∇ρ↓_m + ∇ρ↓_l·∇ρ↑_m)·v2cud)
   

3.5 Important Notes

1. No changes needed for gga_grad=1/2: These two modes already have complete implementations on the develop branch; this task only adds gga_grad=3 and its dispatch
2. gga_grad=3 only affects NSPIN=4: When NSPIN≠4 or there is no magnetic constraint, the original code path should be followed
3. Current builtin SF implementation is simplified: The gradient part reuses the gga_grad=2 m̂·∇m projection rather than the full SF σ transform. Implementing the full SF transform can be considered as future work
4. Libxc path for gga_grad=3 is not implemented: Requires new xc_functional_gga_noncol_sf.cpp/h files with gga_grad==3 dispatch in v_xc_libxc. When use_libxc=true, v_xc does not take the builtin path, so a separate Libxc SF implementation is needed
5. Stress calculation needs verification: The current builtin SF stress implementation has not been thoroughly tested; test cases need to be written on the develop branch

---

4. Testing Requirements

4.1 Basic Functionality Tests

Test	Description
Parameter validation	gga_grad=0 or gga_grad=4 should cause an error exit
Default value	When gga_grad is not set, the default should be 1
NSPIN≠4	When NSPIN=1 or NSPIN=2, gga_grad should not affect results
NSPIN=4	Setting gga_grad=1/2/3 should all complete SCF normally

4.2 Physical Correctness Tests

Test	Description
gga_grad=1 vs 2	The two modes should give different but physically reasonable results for the same system; gga_grad=2 should be smoother at magnetization direction discontinuities
gga_grad=3 results	Should differ from gga_grad=2 (due to different div(h) projection), but should not show AMAG deviation >5% or SCF convergence failure
Energy comparison	Total energies from all three modes should differ by <0.1 eV/atom; gga_grad=3 results should match expectations from the SF paper

Recommended test system: Fe bcc, 2 atoms, NSPIN=4, noncolin=true, PBE functional.

---

5. Key File Index

File	Path	Role
Input parameter	source/module_parameter/input_parameter.h	gga_grad parameter definition
Input parsing	source/module_io/read_input_item_general.cpp	gga_grad parameter reading & validation
v_xc entry point	source/module_hamilt_general/module_xc/xc_functional_vxc.cpp	LDA potential + dispatch to gradcorr/SF
gradcorr	source/module_hamilt_general/module_xc/xc_functional_gradcorr.cpp	GGA gradient correction (gga_grad=1/2 path + gga_grad=3 dispatch)
SF Builtin impl	source/module_hamilt_general/module_xc/xc_functional_gga_noncol_sf_builtin.cpp	gga_grad=3 builtin PBE implementation
SF Builtin header	source/module_hamilt_general/module_xc/xc_functional_gga_noncol_sf_builtin.h	SF builtin declarations
CMake	source/module_hamilt_general/module_xc/CMakeLists.txt	Build configuration

---

6. References

1. G. Scalmani, M. J. Frisch, J. Chem. Theory Comput. 2012, 8, 2193-2196. — Original SF transform paper
2. ABACUS existing gga_grad=1 implementation: xc_functional_gradcorr.cpp NSPIN=4 else branch
3. ABACUS existing gga_grad=2 implementation: xc_functional_gradcorr.cpp NSPIN=4 if branch + xc_functional_vxc.cpp Libxc path

Describe the solution you'd like

raise PR to merge gga_grad3 from dyzheng repo.

Task list only for developers

• Notice possible changes of behavior
• Explain the changes of codes in core modules of ESolver, HSolver, ElecState, Hamilt, Operator or Psi

Notice Possible Changes of Behavior (Reminder only for developers)

No response

Notice any changes of core modules (Reminder only for developers)

No response

Notice Possible Changes of Core Modules (Reminder only for developers)

No response

Additional Context

No response

Task list for Issue attackers (only for developers)

• Review and understand the proposed feature and its importance.
• Research on the existing solutions and relevant research articles/resources.
• Discuss with the team to evaluate the feasibility of implementing the feature.
• Create a design document outlining the proposed solution and implementation details.
• Get feedback from the team on the design document.
• Develop the feature following the agreed design.
• Write unit tests and integration tests for the feature.
• Update the documentation to include the new feature.
• Perform code review and address any issues.
• Merge the feature into the main branch.
• Monitor for any issues or bugs reported by users after the feature is released.
• Address any issues or bugs reported by users and continuously improve the feature.