Multiwfn WFN interoperability validation¶
This notebook is a compact, reader-facing version of the H/O/H₂O validation used to check the atomref-proatoms PROAIM WFN interoperability path. It tests one representative molecular-minus-atomic plane workflow rather than attempting to prove every element, charge state, basis, or angular momentum case.
The distinction between data products is important. Project-native NPZ and structured SCF artifacts are the efficient internal path for PySCF-side work. Radial profiles, and later compact .rad exports, are the preferred density-only representation for stockholder/Hirshfeld-like analyses. A .wfn file is used here as an interoperability container for Multiwfn-like workflows that need Gaussian primitive and orbital coefficient data.
Local execution¶
To run the notebook locally, copy it to local-data/ in a checked-out repository. Place the current Linux Multiwfn executable at local-data/Multiwfn, or inside a local-data/Multiwfn*/ directory, and place settings.ini in local-data/ when a local configuration is needed. The notebook calls Multiwfn with -silent and adds -set local-data/settings.ini when that file is present.
The notebook runs only the fixed H/O/H₂O validation example. It does not regenerate the committed profile or Multiwfn interoperability data products and does not call the public generator.
from pathlib import Path
import numpy as np
import pandas as pd
from atomref_proatoms.dataio.basis import list_basis_bundles, load_basis_nw_text
from atomref_proatoms.dataio.paths import BASIS_ROOT, STATES_FILE
from atomref_proatoms.engines.pyscf_backend import SCFSettings, run_spherical_uks
from atomref_proatoms.exporters.proaim_wfn import (
atom_wfn_filename,
write_atomref_spherical_wfn,
write_mean_field_wfn,
)
from atomref_proatoms.states.state_tables import load_atom_states
from atomref_proatoms.validation.multiwfn_plane import (
parse_multiwfn_point_density_log,
plane_error_metrics,
read_multiwfn_plane,
run_multiwfn_stdin_job,
)
from atomref_proatoms.validation.wfn_density import (
deformation_density_from_templates,
evaluate_density,
parse_wfn_file,
)
PROJECT_ROOT = Path.cwd().parent if Path.cwd().name == "local-data" else Path.cwd()
LOCAL_DATA = PROJECT_ROOT / "local-data"
OUT_DIR = LOCAL_DATA / "multiwfn_plane_validation"
ATOMWFN_DIR = LOCAL_DATA / "atomwfn"
OUT_DIR.mkdir(parents=True, exist_ok=True)
ATOMWFN_DIR.mkdir(parents=True, exist_ok=True)
BASIS_ID = "x2c-QZVPall"
ATOM_SYMBOLS = ["H", "O"]
H2O_WFN = OUT_DIR / "h2o.wfn"
WFN files written from the same density objects¶
The atom files are written from atomref spherical PySCF calculations. Open-shell atom WFNs use spin orbitals directly: alpha orbitals first, beta orbitals second, every occupation at or below one, and an explicit $MOSPIN block containing only labels 1 and 2. No $MOSPIN=3 spatial-orbital reconstruction is used for the fractional/open-shell atom files.
states = load_atom_states(STATES_FILE)
bundles = {bundle.basis_id: bundle for bundle in list_basis_bundles(BASIS_ROOT)}
settings = SCFSettings(
xc="PBE0",
use_x2c=True,
conv_tol=1e-9,
max_cycle=300,
diis_space=12,
diis_start_cycle=1,
grid_level=4,
grid_prune=None,
verbose=3,
)
def neutral_state(symbol: str):
matches = [state for state in states if state.symbol == symbol and state.charge == 0]
if len(matches) != 1:
raise ValueError(f"expected one neutral state for {symbol}, found {len(matches)}")
return matches[0]
atom_export_rows = []
for symbol in ATOM_SYMBOLS:
state = neutral_state(symbol)
run = run_spherical_uks(state, bundles[BASIS_ID], settings=settings)
atom_export_rows.append(
write_atomref_spherical_wfn(
ATOMWFN_DIR / atom_wfn_filename(symbol),
state,
run,
title=f"atomref-proatoms {state.state_id} {BASIS_ID} PBE0 sf-X2C",
)
)
pd.DataFrame(atom_export_rows)
converged SCF energy = -0.501297367976928 <S^2> = 0.75 2S+1 = 2 converged SCF energy = -75.0117962926939 fractional-occupation spherical ensemble: <S^2> is undefined from the 1-RDM; Nalpha = 5 Nbeta = 3 nominal 2S+1 = 3
| file | n_mos | n_primitives | n_centers | occupation_sum | has_fractional_occupations | max_basis_l | writer | has_mospin | max_occupation | ... | strict_atom_mospin_qa_ok | strict_mospin_values_only_1_2 | strict_no_occupation_above_one | strict_alpha_entries_first | state_id | mo_index_gap_ok_for_multiwfn_fallback | mo_index_gap_n_alpha | mo_index_gap_n_beta | mo_index_gap_first_beta_old_label | mo_index_gap_first_beta_new_label | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 0 | local-data/atomwfn/H .wfn | 1 | 38 | 1 | 1.0 | False | 3 | pyscf.tools.wfn_format.write_mo | True | 1.0 | ... | True | True | True | True | H_q0_mult2_nist | NaN | NaN | NaN | NaN | NaN |
| 1 | local-data/atomwfn/O .wfn | 10 | 92 | 1 | 8.0 | True | 4 | pyscf.tools.wfn_format.write_mo | True | 1.0 | ... | True | True | True | True | O_q0_mult3_nist | True | 5.0 | 5.0 | 6.0 | 7.0 |
2 rows × 29 columns
The molecular WFN uses the same basis and electronic-structure convention on a fixed water geometry in the XY plane. The molecule is not part of the proatomic dataset; it is a compact validation system for testing molecular, atomic, and deformation-density paths together.
from pyscf import dft, gto
from pyscf.gto import basis as pyscf_basis
basis_text = load_basis_nw_text(bundles[BASIS_ID])
basis = {symbol: pyscf_basis.parse(basis_text, symb=symbol) for symbol in ATOM_SYMBOLS}
water_geometry_angstrom = [
("O", (0.0, 0.0, 0.0)),
("H", (0.9572, 0.0, 0.0)),
("H", (-0.2399872, 0.92662721, 0.0)),
]
mol_h2o = gto.M(
atom=water_geometry_angstrom,
basis=basis,
charge=0,
spin=0,
unit="Angstrom",
cart=False,
verbose=3,
)
mf_h2o = dft.RKS(mol_h2o)
mf_h2o.xc = "PBE0"
mf_h2o.grids.level = 4
mf_h2o.grids.prune = None
mf_h2o = mf_h2o.sfx2c1e()
mf_h2o.conv_tol = 1e-9
mf_h2o.kernel()
write_mean_field_wfn(
H2O_WFN,
mf_h2o,
title=f"atomref-proatoms H2O validation {BASIS_ID} PBE0 sf-X2C",
)
converged SCF energy = -76.4392193425727
{'file': 'local-data/multiwfn_plane_validation/h2o.wfn',
'n_mos': 5,
'n_primitives': 168,
'n_centers': 3,
'occupation_sum': 10.0,
'has_fractional_occupations': False,
'max_basis_l': 4,
'writer': 'pyscf.tools.wfn_format.write_mo',
'has_mospin': False,
'max_occupation': 2.0,
'mospin_alpha_count': 0,
'mospin_beta_count': 0,
'mospin_spatial_count': 0,
'alpha_occupation_sum': None,
'beta_occupation_sum': None,
'spin_population_alpha_minus_beta': None,
'mo_index_gap_written': False,
'mo_index_gap_reason': 'not_requested'}
Independent WFN read-back¶
The package-side WFN evaluator is intentionally minimal. It reads the saved center, primitive, MO, occupation, coefficient, and $MOSPIN records and evaluates total, alpha, beta, and spin density directly from Cartesian primitives. It is used here to check the WFN boundary; it is not proposed as the primary package storage format.
h2o_wfn = parse_wfn_file(H2O_WFN)
atom_wfns = {
symbol: parse_wfn_file(ATOMWFN_DIR / atom_wfn_filename(symbol))
for symbol in ATOM_SYMBOLS
}
spot_points_bohr = np.array(
[
[0.0, 0.0, 0.0],
[0.5, 0.0, 0.0],
[0.0, 0.5, 0.0],
[1.0, 1.0, 0.0],
[0.3, -0.4, 0.2],
],
dtype=float,
)
rho_wfn = evaluate_density(h2o_wfn, spot_points_bohr)
# A PySCF reference can be evaluated from mf_h2o.make_rdm1() when the SCF object is still in memory.
pd.DataFrame({"x_bohr": spot_points_bohr[:, 0], "rho_wfn": rho_wfn})
| x_bohr | rho_wfn | |
|---|---|---|
| 0 | 0.0 | 334.014203 |
| 1 | 0.5 | 1.029984 |
| 2 | 0.0 | 1.027950 |
| 3 | 1.0 | 0.155616 |
| 4 | 0.3 | 0.968000 |
A representative read-back run gave a maximum absolute H₂O spot-density difference of 2.55e-7 e bohr^-3 between the saved WFN evaluator and direct PySCF density. Away from the oxygen nucleus, the same check was below 7e-9 e bohr^-3.
Multiwfn runs¶
The command streams below reproduce the validation paths used in the local run. They evaluate the O-atom spin diagnostic, the H₂O deformation-density plane, and standalone total-density planes for the molecule and atom WFNs.
DEFORMATION_PLANE_COMMANDS = ["4", "-2", "1", "1", "121,121", "1", "0", "-6", "-5", "q", "q"]
DENSITY_PLANE_COMMANDS = ["4", "1", "1", "121,121", "1", "0", "-6", "-5", "q", "q"]
POINT_CHECK_COMMANDS = ["1", "0.5,0.5,0.5", "1", "q", "q"]
jobs = []
jobs.append(
run_multiwfn_stdin_job(
label="O_point",
input_file_relative_to_cwd=f"atomwfn/{atom_wfn_filename('O')}",
commands=POINT_CHECK_COMMANDS,
cwd=LOCAL_DATA,
log_path=OUT_DIR / "multiwfn_O_point.log",
)
)
jobs.append(
run_multiwfn_stdin_job(
label="delta_xy",
input_file_relative_to_cwd="multiwfn_plane_validation/h2o.wfn",
commands=DEFORMATION_PLANE_COMMANDS,
cwd=LOCAL_DATA,
log_path=OUT_DIR / "multiwfn_delta_xy.log",
plane_output_path=OUT_DIR / "multiwfn_delta_xy_plane.txt",
)
)
pd.DataFrame([job.to_json() for job in jobs])
| label | returncode | executable | cwd | input_file | log | settings_ini | plane_output | commands_passed_via_stdin | command_file_written | |
|---|---|---|---|---|---|---|---|---|---|---|
| 0 | O_point | 0 | /mnt/d/Work/Science/_Codes/DeloneCommons/atomr... | /mnt/d/Work/Science/_Codes/DeloneCommons/atomr... | atomwfn/O .wfn | /mnt/d/Work/Science/_Codes/DeloneCommons/atomr... | /mnt/d/Work/Science/_Codes/DeloneCommons/atomr... | NaN | True | False |
| 1 | delta_xy | 0 | /mnt/d/Work/Science/_Codes/DeloneCommons/atomr... | /mnt/d/Work/Science/_Codes/DeloneCommons/atomr... | multiwfn_plane_validation/h2o.wfn | /mnt/d/Work/Science/_Codes/DeloneCommons/atomr... | /mnt/d/Work/Science/_Codes/DeloneCommons/atomr... | /mnt/d/Work/Science/_Codes/DeloneCommons/atomr... | True | False |
Validation summaries¶
The following compact tables summarize the executed local validation run. The numerical scale is the important result: the WFN read-back and PySCF references agree tightly, Multiwfn sees the same spin-resolved O-atom density, and the Multiwfn deformation plane agrees with the PySCF molecular-minus-atomic reference to a small residual on the 121×121 plane.
| check | points | reference | maximum absolute error | p95 absolute error | RMSE |
|---|---|---|---|---|---|
| H₂O WFN read-back spot density | 5 | direct PySCF density | 2.55e-7 | — | — |
| O atom spin point, total density | 1 | package WFN evaluator | 2.23e-11 | — | — |
| O atom spin point, alpha density | 1 | package WFN evaluator | 3.73e-11 | — | — |
| O atom spin point, beta density | 1 | package WFN evaluator | 1.50e-11 | — | — |
| H₂O deformation plane | 14,641 | direct PySCF molecular-minus-atomic density | 1.51e-4 | 3.08e-7 | 2.0e-6 |
| package WFN deformation plane | 14,641 | direct PySCF molecular-minus-atomic density | 2.55e-7 | 5.03e-10 | 4.24e-9 |
The O-atom point diagnostic at (0.5, 0.5, 0.5) bohr was spin resolved rather than alpha/beta averaged:
| source | total density | alpha density | beta density | spin density |
|---|---|---|---|---|
| package WFN evaluator | 0.502407 | 0.337266 | 0.165141 | 0.172125 |
| Multiwfn point output | 0.502407 | 0.337266 | 0.165141 | 0.172125 |
point = parse_multiwfn_point_density_log(OUT_DIR / "multiwfn_O_point.log")
plane = read_multiwfn_plane(OUT_DIR / "multiwfn_delta_xy_plane.txt")
py_delta = deformation_density_from_templates(h2o_wfn, atom_wfns, plane.points_bohr)
metrics = plane_error_metrics(plane.values, py_delta, prefix="multiwfn_minus_python_wfn")
pd.DataFrame([{**point, **metrics}])
| multiwfn_total_density | multiwfn_alpha_density | multiwfn_beta_density | multiwfn_spin_density | multiwfn_minus_python_wfn_max_abs_error | multiwfn_minus_python_wfn_p95_abs_error | multiwfn_minus_python_wfn_rmse | multiwfn_minus_python_wfn_relative_rmse_vs_max_abs_reference | |
|---|---|---|---|---|---|---|---|---|
| 0 | 0.502407 | 0.337266 | 0.165141 | 0.172125 | 0.000041 | 3.058176e-07 | 7.941030e-07 | 8.317857e-07 |
Interpretation and limits¶
This validation supports the WFN interoperability path used for Multiwfn-style analyses: the saved WFN files preserve the PySCF density, the current Multiwfn build reads the explicit $MOSPIN spin-orbital convention correctly for the O atom, and the deformation-density plane is consistent with the corresponding PySCF reference.
The scope is deliberately narrow. H, O, and H₂O exercise open-shell atom spin typing, molecular WFN export, spherical-AO to Cartesian-primitive WFN evaluation, and g-type primitive handling in the chosen basis. They do not constitute a universal proof for every element, charge, basis family, diffuse formal anion, high angular momentum, or future Multiwfn version. The committed .rad and .wfn interoperability tree should therefore be read under its own manifest and validation contract, while this notebook remains a representative format-boundary check.