Fixed bug in CPUNEP implementation of qNEP

Fixes CPUNEP forces for qNEP (charge_mode 1/2) models, which did not match a finite-difference derivative of the reported energy.

Root cause: charge neutrality is enforced via

Q_{\text{used},i} = Q_{\text{raw},i} - \overline{Q_{\text{raw}}}

a global correction depending on all atoms, not just local descriptors. The force chain rule through that shift is

\frac{\partial E_{\text{coulomb}}}{\partial R_{n_1}}
= \sum_i D_{\text{real},i} \frac{\partial Q_{\text{used},i}}{\partial R_{n_1}}
= \sum_i \left( D_{\text{real},i} - \overline{D_{\text{real}}} \right) \frac{\partial Q_{\text{raw},i}}{\partial R_{n_1}}

so the charge-response derivative D_real (\partial E_{\text{coulomb}}/\partial Q) must also be shifted by -\overline{D_{\text{real}}} before being used to assemble forces. This term exists in GPUMD's nep_charge.cu (zero_mean_D_real) but was missing from src/nepy/nep.cpp, so the unshifted, purely local D_{\text{real},i} term was used instead. The resulting force error was small on large systems (diluted by 1/N) but sizable on small ones (e.g. ~9% relative error on a 3-atom water molecule).

Adds zero_mean_D_real next to the existing zero_total_charge, and a regression test (test_qnep_forces_agree_with_finite_differences in tests/test_cpunep.py) using the repo's existing nep_cm1.txt/nep_cm2.txt fixtures.