pympfit.mpfit

pympfit.mpfit.core

Core MPFIT math functions: A matrix, b vector, solid harmonics.

pympfit.mpfit.core.build_A_matrix(nsite: int, xyzmult: ndarray[tuple[Any, ...], dtype[float64]], xyzcharge: ndarray[tuple[Any, ...], dtype[float64]], r1: float, r2: float, maxl: int, A: ndarray[tuple[Any, ...], dtype[float64]]) → ndarray[tuple[Any, ...], dtype[float64]]

Construct A matrix as in J. Comp. Chem. Vol. 12, No. 8, 913-917 (1991).

Returns 3D array A(i,j,k) where i stands for the specific multipole, j,k for the charges.

pympfit.mpfit.core.build_A_matrix_torch(nsite: int, xyzmult: torch.Tensor, xyzcharge: torch.Tensor, r1: float, r2: float, maxl: int) → torch.Tensor

Build A matrix using PyTorch + sphericart (fully differentiable).

Equivalent to build_A_matrix but supports autograd for Bayesian inference.

pympfit.mpfit.core.build_b_vector(nsite: int, xyzmult: ndarray[tuple[Any, ...], dtype[float64]], xyzcharge: ndarray[tuple[Any, ...], dtype[float64]], r1: float, r2: float, maxl: int, multipoles: ndarray[tuple[Any, ...], dtype[float64]], b: ndarray[tuple[Any, ...], dtype[float64]]) → ndarray[tuple[Any, ...], dtype[float64]]

Construct b vector as in J. Comp. Chem. Vol. 12, No. 8, 913-917 (1991).

pympfit.mpfit.solvers

Solvers for the MPFIT charge fitting problem.

class pympfit.mpfit.solvers.ConstrainedMPFITSolver(conchg: float = 1.0, method: str = 'L-BFGS-B', maxiter: int = 10000, ftol: float = 1e-08, gtol: float = 1e-06)

Bases: ABC

Base class for constrained MPFIT solvers.

initial_guess(state: ConstrainedMPFITState) → ndarray

Return a zero initial guess with shape=(n_params,).

classmethod jacobian(p0: ndarray, state: ConstrainedMPFITState, conchg: float) → ndarray

Return the gradient of the objective function w.r.t. reduced parameters.

Computes the full-space gradient of the multipole fitting error and charge conservation penalty, then projects to the reduced parameter space via createdkaisq.

This is a pure function — it does not mutate state.

Parameters:

• p0 – The current reduced parameter vector with shape=(n_params,). Mapped to full charges via expandcharge.

• state – The frozen problem state containing coordinates, multipoles, equivalence constraints, and per-molecule target charges.

• conchg – The Lagrange-like penalty weight enforcing per-molecule charge conservation.

Return type:

The gradient in the reduced parameter space with shape=(n_params,).

classmethod loss(p0: ndarray, state: ConstrainedMPFITState, conchg: float) → float

Return the value of the constrained MPFIT objective function.

Parameters:

• p0 – The current reduced parameter vector with shape=(n_params,). Mapped to full charges via expandcharge.

• state – The frozen problem state containing coordinates, multipoles, equivalence constraints, and per-molecule target charges.

• conchg – The Lagrange-like penalty weight enforcing per-molecule charge conservation.

Return type:

The scalar value of the objective function.

solve(state: ConstrainedMPFITState) → ndarray

Solve the constrained MPFIT problem and return total charge per atom.

Parameters:

state – The frozen problem state containing coordinates, multipoles, equivalence constraints, and per-molecule target charges.

Return type:

Total charge per atom with shape=(n_atoms,).

class pympfit.mpfit.solvers.ConstrainedMPFITState(xyzmult: ndarray, xyzcharge: ndarray, multipoles: ndarray, quse: ndarray, atomtype: tuple[str, ...], rvdw: ndarray, lmax: ndarray, r1: float, r2: float, maxl: int, atom_counts: tuple[int, ...], molecule_charges: tuple[float, ...])

Bases: object

Immutable container for a constrained MPFIT problem.

atom_counts: tuple[int, ...]

atomtype: tuple[str, ...]

lmax: ndarray

maxl: int

molecule_charges: tuple[float, ...]

multipoles: ndarray

quse: ndarray

r1: float

r2: float

rvdw: ndarray

xyzcharge: ndarray

xyzmult: ndarray

class pympfit.mpfit.solvers.ConstrainedSciPySolver(conchg: float = 1.0, method: str = 'L-BFGS-B', maxiter: int = 10000, ftol: float = 1e-08, gtol: float = 1e-06)

Bases: ConstrainedMPFITSolver

Constrained MPFIT solver using scipy.optimize.minimize.

class pympfit.mpfit.solvers.MPFITSVDSolver(svd_threshold: float = 0.0001)

Bases: MPFITSolver

Solver that uses SVD to find charges that reproduce multipole moments.

This solver processes each multipole site independently using SVD decomposition, accumulating charge contributions via quse_masks.

solve(design_matrix: ndarray, reference_values: ndarray, ancillary_arrays: dict | None = None) → ndarray

Solve for charges using SVD for each multipole site.

Parameters:

• design_matrix – An object array where each element is a site-specific A matrix.

• reference_values – An object array where each element is a site-specific b vector.

• ancillary_arrays – Dictionary containing ‘quse_masks’ - list of boolean masks indicating which atoms affect each multipole site.

Returns:

• Charge values with shape=(n_atoms + n_vsites, 1). When no vsites

• are present, this is simply (n_atoms, 1).

class pympfit.mpfit.solvers.MPFITSolver

Bases: ABC

Base class for MPFIT solvers.

MPFIT uses per-site fitting where each atom has its own design matrix (A) and reference vector (b). The solver receives object arrays containing these per-site matrices and must solve each site independently.

abstractmethod solve(design_matrix: ndarray, reference_values: ndarray, ancillary_arrays: dict | None = None) → ndarray

Solve the MPFIT fitting problem.

Parameters:

• design_matrix – An object array where each element is a site-specific A matrix.

• reference_values – An object array where each element is a site-specific b vector.

• ancillary_arrays – Dictionary containing additional arrays needed for the solver. For SVD solver, must include ‘quse_masks’ - list of boolean masks indicating which atoms affect each multipole site.

Return type:

The set of charge values with shape=(n_atoms + n_vsites, 1)

exception pympfit.mpfit.solvers.MPFITSolverError

Bases: Exception

An exception raised when an MPFIT solver fails.

pympfit.mpfit.solvers.build_quse_matrix(xyzmult: ndarray, xyzcharge: ndarray, rvdw: ndarray) → ndarray

Build binary mask: quse[s, i] = 1 if atom i affects site s.

pympfit.mpfit.solvers.count_parameters(state: ConstrainedMPFITState) → int

Count free parameters after applying atom-type equivalence constraints.

pympfit.mpfit.solvers.createdkaisq(dparam1: ndarray, state: ConstrainedMPFITState) → ndarray

Project full-space gradient to reduced parameter space.

pympfit.mpfit.solvers.expandcharge(p0: ndarray, state: ConstrainedMPFITState) → tuple[ndarray, ndarray]

Map reduced parameters to full charges with atom-type constraints.

Returns:

• allcharge – Per-site charge contributions, shape (n_sites, n_atoms).

• qstore – Total charge per atom, shape (n_atoms,).

pympfit.mpfit._mpfit

pympfit.mpfit._mpfit.generate_constrained_mpfit_charge_parameter(multipole_records: list[MoleculeGDMARecord | MoleculeMBISRecord], molecules: list[Molecule], solver: ConstrainedMPFITSolver | None = None, atom_type_labels: list[list[str]] | None = None, radius: int = 2, equivalize_between_methyl_carbons: bool = True, equivalize_between_methyl_hydrogens: bool = True, equivalize_between_other_heavy_atoms: bool = True, equivalize_between_other_hydrogen_atoms: bool = True) → list[LibraryChargeParameter]

Fit constrained charges for one or more molecules.

When multiple molecules are provided, atoms in different molecules that share the same atom type label are constrained to carry the same charge. For a single molecule, within-molecule symmetry constraints are still applied. Labels can be generated automatically via Morgan fingerprints and within-molecule symmetry, or supplied directly.

Parameters:

• multipole_records – One multipole record (GDMA or MBIS) per molecule (one conformer each).

• molecules – The molecules corresponding to each multipole record, in the same order. The formal charge of each molecule is read from molecule.total_charge.

• solver – The constrained solver to use. Defaults to ConstrainedSciPySolver.

• atom_type_labels – Per-molecule atom type labels. If None, labels are generated automatically via generate_global_atom_type_labels.

• radius – Morgan fingerprint radius passed to generate_global_atom_type_labels when atom_type_labels is None.

• equivalize_between_methyl_carbons – Whether topologically symmetric methyl(ene) carbons (matched by [#6X4H3,#6H4,#6X4H2:1]) in different molecules should be assigned an equivalent charge.

• equivalize_between_methyl_hydrogens – Whether topologically symmetric methyl(ene) hydrogens (attached to a methyl(ene) carbon) in different molecules should be assigned an equivalent charge.

• equivalize_between_other_heavy_atoms – Whether topologically symmetric heavy atoms that are not methyl(ene) carbons in different molecules should be assigned an equivalent charge.

• equivalize_between_other_hydrogen_atoms – Whether topologically symmetric hydrogens that are not methyl(ene) hydrogens in different molecules should be assigned an equivalent charge.

Return type:

One LibraryChargeParameter per molecule.

pympfit.mpfit._mpfit.generate_global_atom_type_labels(molecules: list[Molecule], radius: int = 2, equivalize_between_methyl_carbons: bool = True, equivalize_between_methyl_hydrogens: bool = True, equivalize_between_other_heavy_atoms: bool = True, equivalize_between_other_hydrogen_atoms: bool = True) → list[list[str]]

Generate consistent atom type labels across multiple molecules.

This function implements a hybrid, hierarchical atom typing approach, in which:

1. Within-molecule equivalence is set via get_atom_symmetries for topologically equivalent atoms for the same molecule

2. Cross-molecule equivalence is set via Morgan fingerprints, where local chemical environments establish equivalence.

Two atoms receive the same label only if they share both the same Morgan environment hash and the same within-molecule symmetry group (relative to their Morgan group).

Parameters:

• molecules – The molecules to generate labels for.

• radius – The Morgan fingerprint radius (in bonds) controlling cross-molecule equivalence scope. Smaller values (2-3) match shared substructures (e.g., imidazolium ring across EMIM/BMIM). Larger values require more of the surrounding structure to match.

• equivalize_between_methyl_carbons – Whether topologically symmetric methyl(ene) carbons (matched by [#6X4H3,#6H4,#6X4H2:1]) in different molecules should be assigned an equivalent charge.

• equivalize_between_methyl_hydrogens – Whether topologically symmetric methyl(ene) hydrogens (attached to a methyl(ene) carbon) in different molecules should be assigned an equivalent charge.

• equivalize_between_other_heavy_atoms – Whether topologically symmetric heavy atoms that are not methyl(ene) carbons in different molecules should be assigned an equivalent charge.

• equivalize_between_other_hydrogen_atoms – Whether topologically symmetric hydrogens that are not methyl(ene) hydrogens in different molecules should be assigned an equivalent charge.

Returns:

• A list of label lists, one per molecule. Each inner list contains

• one string label per atom, ordered by atom index.

pympfit.mpfit._mpfit.generate_mpfit_charge_parameter(multipole_records: list[MoleculeGDMARecord | MoleculeMBISRecord], solver: MPFITSolver | None = None, vsite_collection: VirtualSiteCollection | None = None, fit_limit: int | None = None) → LibraryChargeParameter | tuple[LibraryChargeParameter, ndarray]

Generate point charges that reproduce the distributed multipole analysis data.

For multiple conformers, charges are fit independently for each conformer and then averaged.

Parameters:

• multipole_records – The records containing the distributed multipole data (GDMA or MBIS). If multiple records are provided, charges are fit independently for each and averaged.

• solver – The solver to use when finding the charges that minimize the MPFIT loss function. By default, the SVD solver is used.

• vsite_collection – Optional virtual site collection defining extra charge sites beyond atoms. When provided, charges are fit at both atom and vsite positions.

• fit_limit – Optional maximum multipole rank (0-indexed) for fitting. When provided and less than the record’s available rank, the multipole tensor is truncated so only terms up to this rank are used. Allows running the QM/multipole step once at a high rank (e.g., GDMA limit=8 or MBIS max_radial_moment=4) and fitting charges at multiple lower ranks without rerunning the QM calculation.

Returns:

• When vsite_collection is None (LibraryChargeParameter with atom charges.)

• When vsite_collection is provided (Tuple of (LibraryChargeParameter,)

• vsite_charges) where vsite_charges is a numpy array of shape (n_vsites,).

pympfit.mpfit._mpfit.molecule_to_mpfit_library_charge(molecule: Molecule) → LibraryChargeParameter

Create a library charge parameter from a molecule.

Parameters:

molecule – The molecule to create the SMILES pattern from.

Return type:

The library charge parameter with one charge value per atom.