Fix Molden GTO normalization and coordinate conversion

[Stardust0831]

Title: Improve ABACUS Molden output for wavefunction analysis

Summary:
This PR fixes several Molden conversion issues in tools/molden/molden.py while keeping the default workflow unchanged as much as possible.

Changes:

• Correct the primitive Gaussian coefficient convention when writing Molden GTO data. The NAO-to-GTO fit uses unnormalized radial primitives, while Molden readers usually normalize primitive Gaussian functions internally.
• Fix Cartesian_angstrom coordinate conversion. Coordinates in Angstrom are now converted to Bohr for the Molden [Atoms] AU section by dividing by 0.529177210903.
• Add optional multi-start NAO-to-GTO fitting. A single -r value keeps the old single-start behavior; comma-separated -r values enable multi-start fitting and keep the fit with the lowest nonlinear error.
• Add optional Molden [Nval] output via --write-nval. The values are read from UPF z_valence. This option is disabled by default.

Notes:

• The changes are limited to the Molden converter.

Validation:

• Ran the existing molden.py unit tests successfully.
• Checked that default output does not contain [Nval].
• Checked that --write-nval writes C/O/H valence charges for the PhenolDimer test case.
• Checked that Cartesian_angstrom coordinates are written at the correct Bohr scale.

[QuantumMisaka]

@kirk0830 have a look

[kirk0830]

Thanks for your contribution!

This PR LGTM mostly, except the GTO-normalization procedure.

Can you provide some references supporting your idea about normalizing the GTO in molden file? Because as far as I remember, it is also the point I tried to ensure for times when I implemented such a post-processing script, and that is also the reason why I left a warning message informing users of the potential inconsistent number of electrons. But if I was wrong, that would be a good news for users who expect analyzing wavefunction with Multiwfn.

Thanks ;)

[Stardust0831]

@kirk0830
Thanks for the careful review.

I checked this point against the actual Molden handling in Multiwfn and CP2K, and I think the normalization correction is needed.

What the PR Changes

In the current ABACUS fitting code, the radial primitive is built without a primitive normalization factor:

g = c * np.exp(-a * r**2) * r**l

In this PR, before writing the fitted coefficient to the Molden [GTO] section, the coefficient is converted by

def _gto_radial_norm(a, l):
    return math.sqrt(2.0 * (2.0 * a)**(l + 1.5) / math.gamma(l + 1.5))

c /= GTORadials._gto_radial_norm(a, l)

This is not adding an extra normalization to the fitted ABACUS radial function. It converts the fitted coefficient to the convention expected by Molden readers such as Multiwfn.

Multiwfn Behavior

In Multiwfn fileIO.f90, readmolden declares normgau:

real*8,external :: normgau

Source:
https://github.com/digital-chemistry-laboratory/multiwfn-mac-build/blob/source_dist/fileIO.f90#L4078-L4093

It reads the Molden [GTO] exponent and contraction coefficient into primexp and concoeff:

read(10,*) exptmp,concoefftmp
primexp(iprimshell)=exptmp
concoeff(iprimshell)=concoefftmp

Source:
https://github.com/digital-chemistry-laboratory/multiwfn-mac-build/blob/source_dist/fileIO.f90#L4384-L4396

Later, when building internal primitive GTF and MO coefficients, Multiwfn multiplies the Molden coefficient by normgau:

tnormgau=normgau(b(k)%type,b(k)%exp)
temp=concoeff(iexp+l-1)
primconnorm(k)=temp*tnormgau
CO(imo,k)=CObasa(ibasis,imo)*temp*tnormgau

Source:
https://github.com/digital-chemistry-laboratory/multiwfn-mac-build/blob/source_dist/fileIO.f90#L4788-L4796

Therefore, if ABACUS writes the fitted coefficient c directly, Multiwfn effectively reconstructs

c * N(a,l) * exp(-a*r^2) * r^l

instead of the fitted primitive

c * exp(-a*r^2) * r^l

Writing c / N(a,l) avoids this extra factor.

CP2K Behavior

CP2K's Molden writer uses the same idea. In molden_utils.F, CP2K states that it undoes normalization factors included in its internal gccs, and writes the Molden coefficient as

gcc(ipgf, ishell, iset)/(prefac*zet(ipgf, iset)**expzet)

Source:
https://doxygen.cp2k.org/dc/dff/molden__utils_8F_source.html

Relevant source lines are around file lines 228-234 in molden_utils.F; the same location can also be checked in the CP2K repository:
https://github.com/cp2k/cp2k/blob/master/src/molden_utils.F#L228-L234

This supports the same convention conversion: internal GTO coefficients and Molden [GTO] coefficients are not necessarily stored in the same normalization convention, so the writer may need to remove the primitive normalization factor before writing Molden.

Numerical Check

This correction also improves the electron count significantly.

For H2O, the ideal valence electron number is 8. With a 0.1 Bohr density grid:

old script: 7.35492379066527
new script: 8.00488755774764

This does not mean the NAO-to-GTO fitting becomes exact, so I think the original warning should remain. The fitting itself can still introduce residual errors. But the normalization correction avoids the extra primitive normalization factor applied by Molden readers such as Multiwfn, and the H2O density integral is consistent with that interpretation.

[mohanchen]

Thanks for your contribution! LGTM.