- Fix libxc_tools.cpp vtxc: use chr->rho[is][ir] (nspin=4 density)
  instead of libxc-format rho (up/dn spin density)
- Revert libxc_setup.cpp to develop (the accel branch changes
  had map key collision, PBEH param count error, and weakened
  param validation)
- Revert input_parameter.h: restore develop defaults for device,
  sc_* params, xc_exch_ext/xc_corr_ext; keep only gga_grad=3
- Revert docs to develop version

- libxc_pot.cpp: comments for SF gradient decomposition, gga_grad branching
- libxc_tools.cpp: comment for nspin=4 SF dh path with 4-component potential
- xc_pot.cpp: comment for non-libxc SF builtin dispatch
- xc_functional_gga_noncol_sf_builtin.cpp: header comment block with
  Scalmani-Frisch reference, chain-rule formulas, and step-by-step labels

…thm (gga_grad=0)

- Add gga_grad parameter to XC_Functional::v_xc and
  XC_Functional_Libxc::v_xc_libxc instead of reading PARAM.inp.gga_grad
  inside the xc module; update all call sites.
- Pass use_sf/has_mag down to convert_vtxc_v, convert_v_nspin4 and
  convert_v_nspin4_sf; drop PARAM reads from libxc_pot.cpp,
  libxc_tools.cpp, xc_pot.cpp and the SF builtin (net global-dependency
  delta of the PR is now zero).
- gga_grad=0 (new default) selects the original algorithm so existing
  results are unchanged; gga_grad=1 keeps the collinear approximation,
  2/3 select the SF-based paths.
- Revert unrelated default arguments added to libxc_abacus.h interfaces
  and restore the init_func doc comment.
- Forward-declare Charge in xc_functional_gga_noncol_sf_builtin.h.

…_vtxc_v

The gga_grad=2/3 branch of convert_vtxc_v wrote the 4-component dh into
v(2,:) and v(3,:) of a 2-row matrix (out-of-bounds). Map dh back to the
spin-up/down representation with the inverse of the SF decomposition:
v_up -= dh0 + sum_mu dh_mu * m_hat_mu, v_dn -= dh0 - sum_mu dh_mu * m_hat_mu.
The vtxc accumulation over the 4 nspin=4 channels is unchanged.

Methods 2 and 3 share the same chain-rule spin-up/down gradients
grad(rho_up/dn) = (grad(rho) +/- m_hat . grad(m))/2, but differ in the
divergence of h = df/d(grad rho) in the potential:
- gga_grad=2 (projected): v_mu -= m_hat_mu * div((h_up - h_dn)/2),
  dropping the (h_up - h_dn) . grad(m_hat_mu) cross terms;
- gga_grad=3 (full SF): v_mu -= div((h_up - h_dn)/2 * m_hat_mu),
  retaining all cross terms.
Implement the projected variant in cal_dh_sf (LIBXC path) and in
v_xc_ncgga_sf_builtin (built-in path, now also dispatched for
gga_grad=2), and document the formula-level differences in the code
comments and in the gga_grad INPUT documentation.

cal_gdr_sf computed the 'total density' gradient from rho[ir*nspin+0],
which is rho_up in the LIBXC layout, not the total density. Use
chr->rho[0] + chr->rho_core instead, matching the up/down decomposition
in convert_rho_amag_nspin4 and the SF builtin.

Add gga_grad unit tests (test_xc5.cpp): exact values for
compute_mag_part_nspin4, convert_v_nspin4 and convert_v_nspin4_sf; the
projection identity of the gga_grad=2 divergence in cal_dh_sf;
dispatch/crash-free checks and regression anchors for the built-in and
LIBXC SF paths of v_xc. Note that the mocked pointwise FFT cannot
distinguish methods 2 and 3 numerically (verified on real grids).

Also sync the input-main.md index/blank lines with
docs/generate_input_main.py (CI documentation consistency check).

The original algorithm (gga_grad=0) uses the global magnetization
direction inferred from the STRU initial moments (lsign_/ux_): when all
initial moments are collinear, their common direction serves as a global
quantization axis and up/down densities are defined w.r.t. it via
sign(m . u). This matches nspin=2 collinear results exactly but makes
the magnetic potential energy surface discontinuous when moments are
tilted slightly away from collinearity.

gga_grad=1 now forces neg=1 in noncolin_rho (local |m| decomposition),
keeping the magnetic potential energy surface continuous at the price of
exact nspin=2 consistency. gradcorr takes gga_grad explicitly; the
stress path is updated accordingly. Unit tests cover the noncolin_rho
sign mechanism and that gga_grad=1 reproduces the lsign_=false result.