berquist · berquist · Jul 4, 2026 · Jul 4, 2026 · Jul 4, 2026 · Jul 4, 2026
diff --git a/Justfile b/Justfile
@@ -2,7 +2,7 @@ test:
     python -m pytest
 
 test-pixi:
-    pixi run -e test test
+    pixi run --frozen -e test test
 
 pre:
     pre-commit run -a

diff --git a/pymolresponse/data/coords/BC2H4.xyz b/pymolresponse/data/coords/BC2H4.xyz
@@ -1,5 +1,5 @@
 7
-
+in angstrom
 C     -0.0008286620   -0.0004493055   -0.0003241160
 H      0.0005419438   -0.0005510523    1.0881223132
 B      0.7614042004    0.0000912434   -1.2402303969

diff --git a/pymolresponse/data/coords/water_angstrom.xyz b/pymolresponse/data/coords/water_angstrom.xyz
@@ -0,0 +1,5 @@
+3
+
+O         -1.81298        0.53384       -0.01233
+H         -0.82365        0.49649        0.00870
+H         -2.10234       -0.29131        0.45244
diff --git a/pymolresponse/helpers.py b/pymolresponse/helpers.py
@@ -1,12 +1,16 @@
 """Utility functions that are core to calculating physical values."""
 
 from collections.abc import Sequence
-from typing import Any
+from typing import TYPE_CHECKING, Any
 
 import numpy as np
 import periodictable
 
 
+if TYPE_CHECKING:
+    from pymolresponse.indices import Occupations
+
+
 def get_most_abundant_isotope(element: periodictable.core.Element) -> periodictable.core.Isotope:
     most_abundant_isotope = element.isotopes[0]
     abundance = 0
@@ -28,34 +32,41 @@ def get_isotopic_masses(charges: Sequence[int]) -> np.ndarray:
 
 
 def calc_center_of_mass(coords: np.ndarray, masses: np.ndarray) -> np.ndarray:
+    assert len(coords.shape) == 2
+    assert coords.shape[1] == 3
     assert len(masses.shape) == 1
+    assert masses.shape[0] == coords.shape[0]
     denominator = np.sum(masses)
     numerator = np.sum(coords * masses[..., np.newaxis], axis=0)
     return numerator / denominator
 
 
 def calc_center_of_nuclear_charge(coords: np.ndarray, charges: np.ndarray) -> np.ndarray:
+    assert len(coords.shape) == 2
+    assert coords.shape[1] == 3
+    assert len(charges.shape) == 1
+    assert charges.shape[0] == coords.shape[0]
     dummy = np.zeros(3)
     center = nuclear_dipole_contribution(coords, charges, dummy)
     total_charge = np.sum(charges)
     return center / total_charge
 
 
 def nuclear_dipole_contribution(
-    nuccoords: np.ndarray, nuccharges: np.ndarray, origin_in_bohrs: np.ndarray
+    coords: np.ndarray, charges: np.ndarray, origin: np.ndarray
 ) -> np.ndarray:
-    assert len(nuccoords.shape) == 2
-    assert nuccoords.shape[1] == 3
-    assert nuccoords.shape[0] == nuccharges.shape[0]
-    assert origin_in_bohrs.shape == (3,)
-    assert len(nuccharges.shape) in (1, 2)
-    if len(nuccharges.shape) == 1:
-        charges = nuccharges[..., np.newaxis]
+    assert len(coords.shape) == 2
+    assert coords.shape[1] == 3
+    assert coords.shape[0] == charges.shape[0]
+    assert origin.shape == (3,)
+    assert len(charges.shape) in (1, 2)
+    if len(charges.shape) == 1:
+        charges = charges[..., np.newaxis]
     else:
-        assert nuccharges.shape[1] == 1
-        charges = nuccharges
+        assert charges.shape[1] == 1
+        charges = charges
 
-    return np.sum((nuccoords - origin_in_bohrs) * charges, axis=0)
+    return np.sum((coords - origin) * charges, axis=0)
 
 
 def get_uhf_values(
@@ -107,3 +118,18 @@ def mat_uhf_to_packed_rohf(
     for idx, pair_rohf in enumerate(indices_display_rohf):
         mat_rohf[idx] = sum(get_uhf_values(mat_alpha, mat_beta, pair_rohf))
     return mat_rohf
+
+
+def make_density(
+    C: np.ndarray[tuple[int, int, int], np.dtype[np.floating]], occupations: "Occupations"
+) -> np.ndarray[tuple[int, int, int], np.dtype[np.floating]]:
+    nspin, nbasis, _ = C.shape
+    D = np.empty(shape=(nspin, nbasis, nbasis))
+    C_occ_a = C[0, :, : occupations[0]]
+    D[0] = C_occ_a @ C_occ_a.T
+    if nspin == 2:
+        C_occ_b = C[1, :, : occupations[2]]
+        D[1] = C_occ_b @ C_occ_b.T
+    else:
+        D *= 2.0
+    return D
diff --git a/pymolresponse/integrals.py b/pymolresponse/integrals.py
@@ -9,7 +9,10 @@
 
 STARS = "********"
 
-PropertyIntegrals = np.ndarray[tuple[int, int, int], np.dtype[np.floating]]
+PropertyIntegrals = (
+    np.ndarray[tuple[int, int], np.dtype[np.floating]]
+    | np.ndarray[tuple[int, int, int], np.dtype[np.floating]]
+)
 
 
 @unique

diff --git a/pymolresponse/interfaces/psi4/helpers.py b/pymolresponse/interfaces/psi4/helpers.py
@@ -0,0 +1,156 @@
+import numpy as np
+import numpy.linalg as npl
+
+import psi4
+
+from pymolresponse.constants import convfac_au_to_debye
+from pymolresponse.helpers import (
+    calc_center_of_mass,
+    calc_center_of_nuclear_charge,
+    get_isotopic_masses,
+    nuclear_dipole_contribution,
+)
+
+
+# TODO
+def calc_center_of_mass_psi4(pyscfmol) -> np.ndarray:
+    charges = pyscfmol.atom_charges()
+    masses = get_isotopic_masses(charges)
+    coords = pyscfmol.atom_coords()
+    return calc_center_of_mass(coords, masses)
+
+
+def calc_center_of_electronic_charge_psi4(
+    D: np.ndarray, psi4wfn: psi4.core.Wavefunction
+) -> np.ndarray:
+    assert len(D.shape) == 2
+    # no linear dependencies!
+    assert D.shape[0] == D.shape[1]
+    zerovec = np.zeros(3)
+    dipole_at_zerovec = electronic_dipole_contribution_psi4(D, psi4wfn, zerovec)
+    nelec = psi4wfn.tot_electrons()
+    return -dipole_at_zerovec / nelec
+
+
+def electronic_dipole_contribution_psi4(
+    D: np.ndarray, psi4wfn: psi4.core.Wavefunction, origin: np.ndarray
+) -> np.ndarray:
+    assert isinstance(D, np.ndarray)
+    assert len(D.shape) == 2
+    assert D.shape[0] == D.shape[1]
+    assert isinstance(origin, np.ndarray)
+    assert origin.shape == (3,)
+
+    mints = psi4.core.MintsHelper(psi4wfn)
+    M_AO = np.asarray(mints.ao_dipole())
+    assert isinstance(M_AO, np.ndarray)
+    assert len(M_AO.shape) == 3
+    assert M_AO.shape[1] == M_AO.shape[2]
+    assert M_AO.shape[0] == 3
+    M100_AO = M_AO[0, :, :]
+    M010_AO = M_AO[1, :, :]
+    M001_AO = M_AO[2, :, :]
+
+    M100_MO = D * M100_AO
+    M010_MO = D * M010_AO
+    M001_MO = D * M001_AO
+
+    dipole_electronic_atomic_units = -np.asarray(
+        [np.sum(M100_MO), np.sum(M010_MO), np.sum(M001_MO)]
+    )
+    return dipole_electronic_atomic_units
+
+
+def calculate_dipole(
+    coords: np.ndarray,
+    charges: np.ndarray,
+    origin: np.ndarray,
+    D: np.ndarray,
+    psi4wfn: psi4.core.Wavefunction,
+    do_print: bool = False,
+) -> np.ndarray:
+    assert origin.shape == (3,)
+    nuclear_components_au = nuclear_dipole_contribution(coords, charges, origin)
+    electronic_components_au = electronic_dipole_contribution_psi4(D, psi4wfn, origin)
+    total_components_au = electronic_components_au + nuclear_components_au
+    if do_print:
+        nuclear_components_debye = nuclear_components_au * convfac_au_to_debye
+        electronic_components_debye = electronic_components_au * convfac_au_to_debye
+        total_components_debye = total_components_au * convfac_au_to_debye
+        nuclear_norm_au = npl.norm(nuclear_components_au)
+        electronic_norm_au = npl.norm(electronic_components_au)
+        total_norm_au = npl.norm(total_components_au)
+        nuclear_norm_debye = nuclear_norm_au * convfac_au_to_debye
+        electronic_norm_debye = electronic_norm_au * convfac_au_to_debye
+        total_norm_debye = total_norm_au * convfac_au_to_debye
+        print(" origin                        [a.u.]: {} {} {}".format(*origin))
+        print(" dipole components, electronic [a.u.]: {} {} {}".format(*electronic_components_au))
+        print(" dipole components, nuclear    [a.u.]: {} {} {}".format(*nuclear_components_au))
+        print(" dipole components, total      [a.u.]: {} {} {}".format(*total_components_au))
+        print(f" dipole moment, electronic     [a.u.]: {electronic_norm_au}")
+        print(f" dipole moment, nuclear        [a.u.]: {nuclear_norm_au}")
+        print(f" dipole moment, total          [a.u.]: {total_norm_au}")
+        print(
+            " dipole components, electronic [D]   : {} {} {}".format(*electronic_components_debye)
+        )
+        print(" dipole components, nuclear    [D]   : {} {} {}".format(*nuclear_components_debye))
+        print(" dipole components, total      [D]   : {} {} {}".format(*total_components_debye))
+        print(f" dipole moment, electronic     [D]   : {electronic_norm_debye}")
+        print(f" dipole moment, nuclear        [D]   : {nuclear_norm_debye}")
+        print(f" dipole moment, total          [D]   : {total_norm_debye}")
+    return total_components_au
+
+
+def calculate_origin(
+    origin_string: str,
+    coords: np.ndarray,
+    charges: np.ndarray,
+    D: np.ndarray,
+    psi4wfn: psi4.core.Wavefunction,
+    do_print: bool = False,
+) -> np.ndarray:
+    assert isinstance(origin_string, str)
+    origin_string = origin_string.lower()
+    assert origin_string in (
+        "explicitly-set",
+        "zero",
+        "com",
+        "centerofmass",
+        "ecc",
+        "centerofelcharge",
+        "ncc",
+        "centerofnuccharge",
+    )
+    zerovec = np.zeros(3)
+
+    if origin_string == "explicitly-set":
+        if do_print:
+            print(" --- Origin: explicitly-set ---")
+        # TODO
+        origin = zerovec
+    elif origin_string == "zero":
+        if do_print:
+            print(" --- Origin: zero ---")
+        origin = zerovec
+    elif origin_string in ("com", "centerofmass"):
+        if do_print:
+            print(" --- Origin: center of mass ---")
+        masses = get_isotopic_masses(charges[:, 0])
+        origin = calc_center_of_mass(coords, masses)
+    elif origin_string in ("ecc", "centerofelcharge"):
+        if do_print:
+            print(" --- Origin: center of electronic charge ---")
+        origin = calc_center_of_electronic_charge_psi4(D, psi4wfn)
+    elif origin_string in ("ncc", "centerofnuccharge"):
+        if do_print:
+            print(" --- Origin: center of nuclear charge ---")
+        origin = calc_center_of_nuclear_charge(coords, charges)
+    else:
+        msg = f"Unknown origin: {origin_string}"
+        raise RuntimeError(msg)
+
+    if do_print:
+        print(" Calculating the dipole at the requested origin...")
+        calculate_dipole(coords, charges, origin, D, psi4wfn, do_print)
+
+    return origin
diff --git a/pymolresponse/interfaces/psi4/integrals.py b/pymolresponse/interfaces/psi4/integrals.py
@@ -11,10 +11,18 @@
     from pymolresponse.integrals import PropertyIntegrals
 
 
-_unused = "unused"
-DIPOLE = IntegralLabel(_unused)
-DIPVEL = IntegralLabel(_unused)
-ANGMOM_COMMON_GAUGE = IntegralLabel(_unused)
+OVERLAP = IntegralLabel("OVERLAP")
+DIPOLE = IntegralLabel("DIPOLE")
+DIPVEL = IntegralLabel("DIPVEL")
+ANGMOM_COMMON_GAUGE = IntegralLabel("ANGMOM_COMMON_GAUGE")
+QUADRUPOLE = IntegralLabel("QUADRUPOLE")
+
+# TODO What is the default origin?
+
+# TODO how to handle multipole of arbitrary order?  The call is
+# np.asarray(mints.ao_multipoles(order=2, origin=[0.0, 0.0, 0.0])), where
+# "order" is including all previous orders, so this gives dipole (3) +
+# quadrupole (6) for a total of 9.
 
 
 class IntegralsPsi4(Integrals):
@@ -30,12 +38,16 @@ def __init__(self, wfn_or_mol: psi4.core.Molecule | psi4.core.Wavefunction) -> N
         self._mints = psi4.core.MintsHelper(wfn)
 
     def _compute(self, label: IntegralLabel) -> "PropertyIntegrals":
-        if label == DIPOLE:
+        if label == OVERLAP:
+            return np.asarray(self._mints.ao_overlap())
+        elif label == DIPOLE:
             return np.stack([np.asarray(Mc) for Mc in self._mints.ao_dipole()])
         elif label == DIPVEL:
             return np.stack([np.asarray(Mc) for Mc in self._mints.ao_nabla()])
         elif label == ANGMOM_COMMON_GAUGE:
             return np.stack([np.asarray(Lc) for Lc in self._mints.ao_angular_momentum()])
+        elif label == QUADRUPOLE:
+            return np.stack([np.asarray(Mc) for Mc in self._mints.ao_quadrupole()])
         else:
             raise RuntimeError