Fix(xc): fix libxc gga_grad branching for nspin=4 GGA

The gga_grad=3 path (default for nspin=4) incorrectly used the
collinear approximation (gga_grad=1) within libxc instead of the
Scalmani-Frisch gradient method (gga_grad=2). This caused inconsistent
results for noncollinear GGA calculations.

Fix: extend libxc_pot.cpp and libxc_tools.cpp gga_grad==2 checks
to also match gga_grad==3, routing both through the SF gradient path.

Also fix input_parameter.h: restore out_hsr/out_hsk/out_hsr_npz_compat
members removed during cherry-pick that were present in newer develop.

Author: dyzheng

docs/advanced/input_files/input-main.md

@@ -28,6 +28,224 @@ If **"ocp=1"** and **"ocp_set"** is set in INPUT file, the occupations of states
 2. **"nupdown"**
 If **"nupdown"** is set to non-zero, number of spin-up and spin-down electrons will be fixed, and Fermi energy level will split to E_Fermi_up and E_Fermi_down. By the way, total magnetization will also be fixed, and will be the value of **"nupdown"**.
 
+## DeltaSpin (Spin-Constrained DFT)
+
+DeltaSpin is a spin-constrained DFT method that allows users to constrain the magnetic moments on individual atoms to target values during self-consistent field (SCF) calculations. This is useful for studying magnetic excitations, non-collinear magnetic structures, and systems where the magnetic ground state is not known a priori.
+
+The theoretical foundation and implementation details can be found in:
+
+- Cai Z, Wang K, Xu Y, et al., "A self-adaptive first-principles approach for magnetic excited states," *Quantum Frontiers* 2.1 (2023): 21. [DOI: 10.1007/s44214-023-00050-z](https://doi.org/10.1007/s44214-023-00050-z)
+- Zheng D, Peng X, Huang Y, et al., "Integrating deep-learning-based magnetic model and non-collinear spin-constrained method: methodology, implementation and application," *npj Computational Materials* (2026).
+
+### Enabling DeltaSpin
+
+Set `sc_mag_switch 1` in the INPUT file. DeltaSpin is supported for both PW (`basis_type = pw`) and LCAO (`basis_type = lcao`) basis sets, with `nspin = 2` (collinear) or `nspin = 4` (non-collinear).
+
+### Specifying Target Magnetic Moments in STRU
+
+Target magnetic moments and constraint flags are specified per atom in the `ATOMIC_POSITIONS` section of the STRU file, using the `mag` (or `magmom`), `sc`, `lambda`, `angle1`, and `angle2` keywords after the atomic coordinates.
+
+#### Collinear (nspin=2)
+
+For collinear spin, only the z-component of the magnetic moment is constrained:
+
+```
+ATOMIC_POSITIONS
+Direct
+
+Fe
+0.0
+2
+0.00  0.00  0.00  mag  2.0   sc 1
+0.51  0.51  0.51  mag  -2.0  sc 1
+```
+
+- `mag 2.0`: target magnetic moment of 2.0 $\mu_B$ along z-axis
+- `sc 1`: constrain the z-component (1 = constrained, 0 = unconstrained)
+
+#### Non-collinear (nspin=4), vector form
+
+For non-collinear spin, specify the magnetic moment as a vector (mx, my, mz):
+
+```
+ATOMIC_POSITIONS
+Direct
+
+Fe
+0.0
+2
+0.00  0.00  0.00  mag  2.0  0.0  0.0  sc 1 1 1
+0.51  0.51  0.51  mag  0.0  0.0  -2.0  sc 1 1 1
+```
+
+- `mag 2.0 0.0 0.0`: target moment vector in Cartesian coordinates ($\mu_B$)
+- `sc 1 1 1`: constrain x, y, z components respectively
+
+#### Non-collinear (nspin=4), angle form
+
+Alternatively, use `angle1` (polar angle $\theta$) and `angle2` (azimuthal angle $\phi$) in degrees to specify the direction:
+
+```
+0.00  0.00  0.00  mag 2.0  angle1 0  angle2 0    sc 1 1 1
+0.51  0.51  0.51  mag 2.0  angle1 180  angle2 0  sc 1 1 1
+```
+
+The Cartesian components are computed as:
+- $m_z = |\mathbf{m}| \cos\theta$
+- $m_x = |\mathbf{m}| \sin\theta \cos\phi$
+- $m_y = |\mathbf{m}| \sin\theta \sin\phi$
+
+#### Providing initial Lagrange multipliers
+
+Initial lambda values (in eV/$\mu_B$) can be provided via the `lambda` keyword to accelerate convergence:
+
+```
+0.00  0.00  0.00  mag 2.0  lambda 0.01 0.0 0.0  sc 1 1 1
+```
+
+A single value sets $\lambda_z$; three values set $\lambda_x$, $\lambda_y$, $\lambda_z$.
+
+#### Partial constraints
+
+Set `sc 0` for unconstrained components. For example, to constrain only the direction but not the magnitude (use with `sc_direction_only`):
+
+```
+0.00  0.00  0.00  mag 2.0  0.0  0.0  sc 1 1 0
+```
+
+### DeltaSpin INPUT Parameters
+
+| Parameter | Type | Default | Description |
+|-----------|------|---------|-------------|
+| `sc_mag_switch` | Boolean | False | Enable DeltaSpin |
+| `sc_thr` | Real | 1.0e-6 | Convergence criterion for lambda loop (RMS, in $\mu_B$) |
+| `nsc` | Integer | 100 | Maximum number of lambda iterations |
+| `nsc_min` | Integer | 2 | Minimum number of lambda iterations |
+| `sc_scf_nmin` | Integer | 2 | Minimum outer SCF iterations before starting lambda loop |
+| `alpha_trial` | Real | 0.01 | Initial trial step size for lambda (eV/$\mu_B^2$) |
+| `sccut` | Real | 3.0 | Maximum step size for lambda (eV/$\mu_B$) |
+| `sc_drop_thr` | Real | 1.0e-2 | Convergence ratio threshold for adaptive lambda loop |
+| `sc_scf_thr` | Real | 1.0e-4 | Density error threshold for entering lambda loop |
+| `sc_direction_only` | Boolean | False | Constrain only the direction, not the magnitude |
+| `sc_lambda_strategy` | String | bfgs | Lambda update strategy (see below) |
+| `decay_grad_switch` | Boolean | False | Enable gradient-based early exit |
+
+For full parameter details, see the [Spin-Constrained DFT](../input_files/input-main.md#spin-constrained-dft) section of the input keyword list.
+
+### Lambda Update Strategies
+
+The `sc_lambda_strategy` parameter controls how the Lagrange multipliers $\lambda$ are updated during the lambda loop:
+
+- **`bfgs`** (default): BFGS quasi-Newton method with line search. Robust and well-tested for both PW and LCAO. Uses `alpha_trial` and `sccut` to control step size.
+
+- **`linear_response`**: Linear response method (Scheme B). Estimates the magnetic susceptibility $\chi$ from the history of $(\lambda, M)$ pairs and performs a one-step Newton-like update: $\Delta\lambda = \beta (M_{\text{target}} - M) / \chi$, where $\beta$ is a mixing parameter.
+
+- **`augmented_lagrangian`**: Augmented Lagrangian method (Scheme C). Uses a penalty parameter $\mu$ that grows over iterations: $\lambda_{\text{new}} = \lambda + \mu (M - M_{\text{target}})$. The penalty increases until convergence is achieved.
+
+- **`hybrid_delayed`**: Hybrid delayed update (Scheme D). Two-phase approach: in the early phase (SCF not yet converged), lambda updates are gentle; in the late phase (SCF nearly converged), augmented Lagrangian updates are applied.
+
+### Direction-Only Mode
+
+When `sc_direction_only 1` is set, only the **direction** of the magnetic moment is constrained to match the target, while the magnitude is allowed to vary freely. This is useful for:
+
+- Studying magnetic anisotropy energy surfaces
+- Cases where the moment magnitude is determined by the electronic structure
+- Converging to the easy-axis direction without fixing the moment size
+
+In this mode, the lambda vector is projected to be perpendicular to the target moment direction at each iteration, ensuring it can only rotate the magnetization, not stretch it.
+
+### Combining DeltaSpin with DFT+U
+
+DeltaSpin can be combined with DFT+U for strongly correlated systems. When both `sc_mag_switch` and `dft_plus_u` are enabled:
+
+1. DFT+U occupation update runs first in each SCF iteration
+2. DeltaSpin lambda loop runs after, constraining the magnetic moments
+3. The DFT+U-corrected Hamiltonian is used by the lambda loop
+
+Example INPUT for PW DFT+U + DeltaSpin:
+
+```
+INPUT_PARAMETERS
+calculation         scf
+basis_type          pw
+ecutwfc             50
+nspin               2
+dft_plus_u          1
+orbital_corr        -1 2
+hubbard_u           0.0 4.0
+sc_mag_switch       1
+sc_thr              1.0e-6
+sc_scf_thr          1.0e-4
+sc_lambda_strategy  bfgs
+```
+
+### Example: Collinear antiferromagnetic Fe
+
+INPUT file:
+
+```
+INPUT_PARAMETERS
+calculation         scf
+basis_type          pw
+ecutwfc             50
+nspin               2
+sc_mag_switch       1
+sc_thr              1.0e-6
+```
+
+STRU file:
+
+```
+ATOMIC_SPECIES
+Fe 55.845 Fe.upf
+
+LATTICE_CONSTANT
+8.190
+
+LATTICE_VECTORS
+ 1.00  0.50  0.50
+ 0.50  1.00  0.50
+ 0.50  0.50  1.00
+
+ATOMIC_POSITIONS
+Direct
+
+Fe
+0.0
+2
+0.00  0.00  0.00  mag  2.0  sc 1
+0.51  0.51  0.51  mag  -2.0  sc 1
+```
+
+### Example: Non-collinear constrained moments
+
+INPUT file:
+
+```
+INPUT_PARAMETERS
+calculation         scf
+basis_type          pw
+ecutwfc             50
+nspin               4
+noncolin            1
+sc_mag_switch       1
+sc_direction_only   1
+sc_lambda_strategy  bfgs
+```
+
+STRU file:
+
+```
+ATOMIC_POSITIONS
+Direct
+
+Fe
+0.0
+2
+0.00  0.00  0.00  mag  2.0  0.0  0.0  sc 1 1 0
+0.51  0.51  0.51  mag  0.0  0.0  2.0  sc 1 1 0
+```
+
 ## Noncollinear Spin Polarized Calculations
 The spin non-collinear polarization calculation corresponds to setting **"noncolin 1"**, in which case the coupling between spin up and spin down will be taken into account.
 In this case, nspin is automatically set to 4, which is usually not required to be specified manually.
@@ -48,6 +266,53 @@ Note: different settings for "noncolin" and "lspinorb" correspond to different c
 
 **Special case**: `noncolin=0, lspinorb=1` is commonly used for non-magnetic materials with SOC effects (e.g., topological insulators, semiconductors with spin-orbit splitting). In this case, the magnetization is NOT automatically set, implying no magnetic moments in the system.
 
+### Choosing gga_grad for Noncollinear Calculations
+
+When performing noncollinear spin calculations (`nspin=4`), the `gga_grad` parameter controls how the gradient of the magnetization is computed in GGA exchange-correlation functionals. This is critical for obtaining correct magnetic anisotropy energies and spin-orbit coupling effects.
+
+#### Available Options
+
+| gga_grad | Gradient Method | Description | When to Use |
+|----------|----------------|-------------|-------------|
+| 0 | Standard | Conventional GGA gradient (collinear approximation) | Collinear calculations only (`nspin<4`) |
+| 1 | $\nabla\|\mathbf{m}\|$ | Gradient of magnetization magnitude | Noncollinear with SOC, when only magnitude variation matters |
+| 2 | $\nabla\hat{\mathbf{m}}$ | Gradient of magnetization direction (default for `nspin=4`) | **Recommended** for most noncollinear calculations with SOC |
+| 3 | Built-in SF | Special functional for noncollinear SOC (requires `v_xc_ncgga_sf_builtin`) | Advanced noncollinear SOC calculations |
+
+#### Recommendations
+
+**For noncollinear calculations with SOC (`lspinorb=1`)**:
+- **Use `gga_grad=2`** (default when `nspin=4`): This computes the gradient of the magnetization direction unit vector $\hat{\mathbf{m}}$, which is physically correct for spin-orbit coupling effects.
+- This is automatically set when you specify `nspin=4` or `lspinorb=1`, but you can explicitly set it to be clear.
+
+**For noncollinear calculations without SOC (`noncolin=1`, `lspinorb=0`)**:
+- **Use `gga_grad=2`** for spin spiral calculations or frustrated magnetic systems
+- The gradient of the magnetization direction is still important even without SOC
+
+**For collinear calculations (`nspin=1` or `nspin=2`)**:
+- **Use `gga_grad=0`** (default): Standard collinear GGA gradient
+- Do NOT use `gga_grad=1`, `2`, or `3` for collinear calculations
+
+#### Example INPUT for Noncollinear SOC
+
+```
+INPUT_PARAMETERS
+calculation         scf
+basis_type          pw
+ecutwfc             50
+nspin               4
+noncolin            1
+lspinorb            1
+gga_grad            2
+```
+
+#### Important Notes
+
+1. **Automatic setting**: When `nspin=4`, ABACUS automatically sets `gga_grad=2` if not explicitly specified
+2. **Physical correctness**: Using the wrong `gga_grad` value can lead to incorrect magnetic anisotropy energies and spin textures
+3. **Compatibility**: `gga_grad=3` requires the built-in special functional and may not be available for all XC functionals
+4. **Performance**: `gga_grad=2` has similar computational cost to `gga_grad=0` for most systems
+
 ## For the continuation job
 - Continuation job for "nspin 1" need file "SPIN1_CHG.cube" which is generated by setting "out_chg=1" in task before. By setting "init_chg file" in new job's INPUT file, charge density will start from file but not atomic. 
 - Continuation job for "nspin 2" need files "SPIN1_CHG.cube" and "SPIN2_CHG.cube" which are generated by "out_chg 1" with "nspin 2", and refer to spin-up and spin-down charge densities respectively. It should be note that reading "SPIN1_CHG.cube" only for the continuation target magnetic moment job is not supported now.