Fix numerical instability in D3 reference-weight derivatives

[Growl1234]

Fixes #8106.

The current D3 implementation calculates the derivatives of normalized reference weights as follows:

double derivative = raw_derivative * normalization
    - raw[reference] * normalization_derivative
        * normalization * normalization;

This expression is mathematically correct. However, when the coordination number (CN) is much higher than the reference values, the raw weights become extremely small. Floating-point underflow and cancellation can then produce incorrect derivatives, leading to unphysically large forces. This issue can be reproduced with an Mg slab (CN ≈ 12).

This PR replaces the calculation with log-sum-exp normalization and a numerically stable analytical derivative. By subtracting the maximum exponent before evaluating the weights, at least one scaled weight is exactly 1, preventing the normalization sum from underflowing. The derivative is then expressed directly in terms of normalized weights and the weighted average reference CN, avoiding the subtraction of large intermediate terms that should mathematically cancel. This ensures numerically stable derivatives even when the actual CN is far beyond the reference range. A regression test for high-CN Mg is also added.

[Copilot]

Copilot was unable to review this pull request because the user who requested the review has reached their quota limit.

[ZhouXY-PKU]

Maybe we can link an issue with wrong results of Mg slabs?

[ZhouXY-PKU]

Nice job!