Fix multipole order forwarding and ddCOSMO radius derivatives

[robin-dahl]

Summary

This PR fixes two independent problems affecting multipolar electrostatics and combined Cartesian/radius derivatives.

Multipolar electrostatics

multipole_electrostatics previously passed a hard-coded mmax = 0 to the underlying _0, _1, and _2 routines. Consequently, only monopoles were evaluated even when the caller supplied higher-order multipoles.

The wrapper now forwards the requested mmax in all three paths.

COSMO radius derivatives

contract_grad_U previously updated either the Cartesian force or the radius derivative, depending on whether dr was present. Therefore, requesting radius derivatives caused the Cartesian force to omit the entire U contribution. The routine now always accumulates the Cartesian contribution and additionally accumulates dr when requested.

The combined COSMO derivative now also:

• applies the required -1/2 scaling to the radius derivatives.
• includes the missing radial zeta contribution through zeta_grad_dr.

Finally, the gradient tests now use a separate scratch-force array when evaluating radius derivatives. This prevents the corrected combined-output behavior from modifying Cartesian-gradient results that are checked later.

The existing tests did not expose this discrepancy because they validated the Cartesian and radius contributions separately through helper routines that already applied the correct -1/2 scaling and radial zeta term, rather than testing the affected combined cosmo_solvation_force_dr_terms routine.