From b21f8fb8db85b8443c4c6fdfb82322aa39d685e9 Mon Sep 17 00:00:00 2001
From: Liam Marsh <liam.marsh@epfl.ch>
Date: Thu, 11 Jun 2026 16:31:28 +0200
Subject: [PATCH 1/2] add psi4 to orbital ordering styles

---
 qstack/reorder.py | 6 +++++-
 1 file changed, 5 insertions(+), 1 deletion(-)

diff --git a/qstack/reorder.py b/qstack/reorder.py
index da9255f..f387917 100644
--- a/qstack/reorder.py
+++ b/qstack/reorder.py
@@ -68,6 +68,10 @@ def __getitem__(self, l):
         'gpr': MagneticOrder(
             l_checker=lambda _: False,
             ),
+        'psi4': MagneticOrder(
+            l_checker=lambda l: l>0,
+            m_generator=_orca_get_m,
+            ),
         'orca': MagneticOrder(
             l_checker=lambda l: l>0,
             m_generator=_orca_get_m,
@@ -186,7 +190,7 @@ def reorder_ao(mol, vector, src='pyscf', dest='gpr'):
             If None, returns the indices to reorder and sign multipliers for an 1D vector
             to use as `x = x[idx]*sign`.
         src (str): Current convention. Defaults to 'pyscf'.
-        dest (str): Convention to convert to (available: 'pyscf', 'gpr', 'orca', 'gaussian', 'turbomole'). Defaults to 'gpr'.
+        dest (str): Convention to convert to (available: 'pyscf', 'gpr', 'orca', 'gaussian', 'turbomole', 'psi4'). Defaults to 'gpr'.
 
     Returns:
         numpy.ndarray: Reordered vector or matrix, or tuple of (idx (numpy.ndarray), sign (numpy.ndarray)) if vector is None.

From 38715d79fe474c9890726bab2fb51645841e5e67 Mon Sep 17 00:00:00 2001
From: Liam Marsh <liam.marsh@epfl.ch>
Date: Tue, 30 Jun 2026 20:59:05 +0200
Subject: [PATCH 2/2] add psi4 example

---
 examples/example_psi4.psi4-in | 48 +++++++++++++++++++++++++++++++++++
 1 file changed, 48 insertions(+)
 create mode 100644 examples/example_psi4.psi4-in

diff --git a/examples/example_psi4.psi4-in b/examples/example_psi4.psi4-in
new file mode 100644
index 0000000..35fc238
--- /dev/null
+++ b/examples/example_psi4.psi4-in
@@ -0,0 +1,48 @@
+#!/usr/bin/env psi4
+
+# example: compute the density matrix from a psi4 calculation, use qstack to extract the corresponding density function,
+# then export both to disk
+
+# (first, make sure that the version of python that powers psi4 can import pyscf and qstack)
+import numpy as np
+import pyscf
+import qstack.fields
+import qstack.reorder
+
+name_base="compound_path_without_xyz_extension"
+
+set {
+    basis cc-pVQZ
+    scf_type DF
+    freeze_core True
+}
+
+# create molecule-representing objects for psi4 and pyscf
+mol_psi4 = psi4.geometry(open(name_base+".xyz",'r').read())
+mol_pyscf = pyscf.M(atom=pyscf.gto.mole.fromfile(name_base+'.xyz'), basis='cc-pVQZ')
+
+
+# note: this exact level of theory is not something
+energy, wfn = energy('wb97x-d', molecule=mol_psi4, return_wfn=True)
+wfn.to_file(name_base+".psi4-wfn")  # optional
+
+# if you already stored the wavefunction, you can also reload it from disk
+#wfn = psi4.core.Wavefunction.from_file(name_base+".psi4-wfn")
+
+
+# sanity check! is the basis set the same in pyscf and psi4?
+ovlp_pyscf = mol_pyscf.intor('int1e_ovlp')
+ovlp_psi4 = wfn.S().to_array(dense=True)
+ovlp_psi4 = qstack.reorder.reorder_ao(mol_pyscf, ovlp_psi4, src='psi4', dest='pyscf')
+assert np.allclose(ovlp_psi4, ovlp_pyscf)
+
+
+
+# combine the spin-up and spin-down density matrices,
+# and also make them regular numpy arrays in case symmetry-based optimisations exist
+rdm_psi4 = wfn.Da().to_array(dense=True)+wfn.Db().to_array(dense=True)
+
+rdm_pyscf = qstack.reorder.reorder_ao(mol_pyscf, rdm_psi4, src='psi4', dest='pyscf')
+auxmol_pyscf, rho = qstack.fields.decomposition.decompose(mol_pyscf, rdm_pyscf, 'cc-pVQZ-jkfit')
+np.save(name_base+'.psi4-rho.npy', rho)
+np.save(name_base+'.psi4-rdm.npy', rdm_pyscf)