From d50c61ac6f9ba180fdff6189514d65aafd71406b Mon Sep 17 00:00:00 2001 From: juanjoaucar Date: Sat, 25 Jul 2026 19:12:47 -0300 Subject: [PATCH 1/4] PV-energy from density matrix --- pyscf/pv/energy.py | 223 +++++++++++++++++++++++++++++--- pyscf/pv/test/test_energy-pv.py | 10 +- 2 files changed, 213 insertions(+), 20 deletions(-) diff --git a/pyscf/pv/energy.py b/pyscf/pv/energy.py index 3b32b00f3..dc61f9207 100755 --- a/pyscf/pv/energy.py +++ b/pyscf/pv/energy.py @@ -72,6 +72,51 @@ def fc_integrals(mol, mf, atom, **kwargs): return fc_integrals + + +def gamma5_fc_integrals(mol, mf, atom, **kwargs): + """ + Compute the integrals for PV contributions to energy for a given atom. + The integrals are calculated in the four-component spinor basis. + Point nuclear charge distribution is assumed. + + :param mol: Molecule object containing information about the molecular system. + :type mol: pyscf.gto.Mole + :param mf: Mean-field object, obtained from a Dirac-Hartree-Fock (DHF) calculation. + :type mf: pyscf.scf.hf.SCF + :param atom: Index of the atom for which the integrals are computed. + :type atom: int + :return: Matrix of integrals (4-component spinor form) for the specified atom. + :rtype: numpy.ndarray (complex128), shape (n4c, n4c) + + Example: + mol = gto.Mole() + # Setup mol with geometry and basis set... + mf = scf.DHF(mol) + mf.kernel() + atom_index = 0 + fc_matrix = PV_integrals(mol, mf, atom_index) + """ + c = lib.param.LIGHT_SPEED + n4c = mf.mo_coeff.shape[0] + n2c = n4c // 2 + + + # Get the Fermi Contact integrals for the selected atom + fc_ao = fc_integrals(mol, mf, atom) + + # Prepare the empty matrix for the FC*Gamma5 integrals (LL, LS, SL, SS blocks) + gamma5_fc_integrals = numpy.zeros((n4c, n4c), dtype=numpy.complex128) + + # LS and SL blocks + gamma5_fc_integrals[:n2c, n2c:] = fc_ao[:n2c, n2c:].copy() + gamma5_fc_integrals[n2c:, :n2c] = fc_ao[n2c:, :n2c:].copy() + + + + + return gamma5_fc_integrals + def fc_expval(mol, mf, atom): """ Calculate the Fermi Contact (FC) expectation values for each occupied orbital in a molecule, @@ -117,18 +162,32 @@ def fc_expval(mol, mf, atom): fc_ao_SS = fc_ao[n2c:, n2c:] for k in range(nocc): - expval_LL_k = numpy.einsum('pi,ij,qj->pq', Lo[:, k].conjugate(), fc_ao_LL, Lo[:, k]) - expval_SS_k = numpy.einsum('pi,ij,qj->pq', So[:, k].conjugate(), fc_ao_SS, Lo[:, k]) + expval_LL_k = numpy.einsum('i,ij,j->', Lo[:,k].conj(), fc_ao_LL, Lo[:,k]) + expval_SS_k = numpy.einsum('i,ij,j->', So[:,k].conj(), fc_ao_SS, So[:,k]) + expval_perorb[k] = expval_LL_k + expval_SS_k + return fac * expval_perorb.real -def Epv_atom(mol, mf, atom_index): +def Epv_atom(mol, mf, atom_index, dm=None): """ Calculate the parity-violating (PV) contribution to energy for a given atom in a molecule. + The function supports: + - orbital expectation values from DHF orbitals + - density matrix contraction in AO or MO basis + + For a one-electron operator: + E_PV = Tr(D * h_PV) + + If dm is provided: + - AO density matrices are transformed to MO basis + - MO density matrices (e.g. MP2/CCSD PySCF densities) are used directly + + It then returns the contributions from each occupied orbital. (First term of Eq. 4 - https://doi.org/10.1002/wcms.1396) @@ -151,13 +210,16 @@ def Epv_atom(mol, mf, atom_index): """ n4c, nmo = mf.mo_coeff.shape n2c = n4c // 2 - nNeg = nmo // 2 # Molecular orbitals of negative energy - nocc = mf.mol.nelectron + # Molecular orbitals of negative energy + nNeg = nmo // 2 - # Get the positive components of the MO spinors - Lo = mf.mo_coeff[:n2c, nNeg:nNeg + nocc] - So = mf.mo_coeff[n2c:, nNeg:nNeg + nocc] + # Number of occupied electrons + nocc = mf.mol.nelectron + + # ------------------------------------------------------------ + # Nuclear weak charge prefactor + # ------------------------------------------------------------ # Atom masses and atomic numbers masses = mf.mol.atom_mass_list(isotope_avg=False) @@ -166,29 +228,115 @@ def Epv_atom(mol, mf, atom_index): # Neutrons per atom neutrons = masses - atomic_numbers - S2THETAW = 0.23122 # AS DIRAC24 (CODATA 2018) + S2THETAW = 0.23122 # (CODATA 2018) # Weak charge of the nucleus of the selected atom QW = (1 - 4 * S2THETAW) * atomic_numbers[atom_index] - neutrons[atom_index] + # Prefactor + fac = (2.2225 * 10 ** (-14) * QW) / (2 * numpy.sqrt(2)) + + # ------------------------------------------------------------ + # PV operator in AO basis + # ------------------------------------------------------------ + # Get the Fermi Contact integrals for the selected atom fc_ao = fc_integrals(mol, mf, atom_index) + # Get the Gamma5 Fermi-Contact integrals for the selected atom (in AO basis) + gamma5_fc_ao = gamma5_fc_integrals(mol, mf, atom_index) + + # ------------------------------------------------------------ + # Orbital expectation values (4c-DHF or 4c-DFT reference) + # ------------------------------------------------------------ + + # Get the positive components of the MO spinors + Lo = mf.mo_coeff[:n2c, nNeg:nNeg + nocc] + So = mf.mo_coeff[n2c:, nNeg:nNeg + nocc] + # Expectation values expval_LS = numpy.einsum('ij,ji->i', Lo.conjugate().T @ fc_ao[:n2c, n2c:], So) expval_SL = numpy.einsum('ij,ji->i', So.conjugate().T @ fc_ao[n2c:, :n2c], Lo) - - # Sum LS and SL terms expval_perorb = expval_LS + expval_SL - # Prefactor - fac = (2.2225 * 10 ** (-14) * QW) / (2 * numpy.sqrt(2)) + # ------------------------------------------------------------ + # Density matrix contraction + # ------------------------------------------------------------ - return fac * expval_perorb.real + expval_dm = None + + if dm is not None: + + # Positive-energy MO coefficient matrix + C = mf.mo_coeff[:, nNeg:] + + nmo_pos = C.shape[1] + + # Transform operator AO -> MO + # Operator transformation: + # + # hᴹᴼ = C† hᴬᴼ C + # + # where: + # C : MO coefficient matrix + # hᴬᴼ : operator matrix in AO basis + # hᴹᴼ : operator matrix in MO basis + gamma5_fc_mo = C.conj().T @ gamma5_fc_ao @ C + + # -------------------------------------------------------- + # Determine density matrix representation + # -------------------------------------------------------- + + if dm.shape[0] == n4c: + # Transform AO density matrix to MO + # 𝑫ᴹᴼ = C† S 𝑫ᴬᴼ S C + # + # where: + # C : MO coefficient matrix + # S : AO overlap matrix + # 𝑫 : one-particle density matrix + S = mf.get_ovlp() + dm_mo = C.conj().T @ S @ dm @ S @ C + + elif dm.shape[0] == nmo_pos: + # Density already in MO basis + # (MP2/CCSD densities) + dm_mo = dm + else: + raise ValueError( + "Density matrix dimension does not match AO or positive-energy MO space" + ) + + # Total expectation value + # E = Tr(𝑫 h) + expval_dm = numpy.trace(dm_mo @ gamma5_fc_mo) + + # Orbital decomposition in MO basis + expval_dm_perorb = numpy.einsum( + 'ij,ji->i', + gamma5_fc_mo, + dm_mo + ) + # ------------------------------------------------------------ + # Return + # ------------------------------------------------------------ -def Epv_molecule(mol, mf): + if expval_dm is None: + return fac * expval_perorb.real + + else: + return fac * expval_dm_perorb.real + + + + + + + + +def Epv_molecule(mol, mf, dm=None): """ Calculate the weak charge parity-violating (PV) contributions for all atoms in the molecule within a punctual nuclear charge distribution model. @@ -196,6 +344,14 @@ def Epv_molecule(mol, mf): This function iterates over all atoms in the molecule and computes the PV contributions for each atom. + If dm is None: + returns orbital contributions from the reference DHF/DFT orbitals. + + If dm is provided: + returns contributions obtained by contracting the PV operator with + the supplied density matrix (AO or MO representation). + + :param mol: Molecule object containing information about the molecular system. :type mol: pyscf.gto.Mole :param mf: Mean-field object, typically obtained from a Dirac-Hartree-Fock (DHF) calculation. @@ -205,9 +361,38 @@ def Epv_molecule(mol, mf): :rtype: numpy.ndarray (real), shape (n_atoms, n_occ) """ - nocc = mf.mol.nelectron - result = numpy.zeros((mf.mol.natm, nocc)) - for i in range(mf.mol.natm): - result[i, :] = Epv_atom(mol, mf, i) + natm = mf.mol.natm + + # Determine output size from the calculation type + if dm is None: + n4c, nmo = mf.mo_coeff.shape + nocc = mf.mol.nelectron + ncol = nocc + + else: + n4c, nmo = mf.mo_coeff.shape + nNeg = nmo // 2 + nmo_pos = nmo - nNeg + + # AO density + if dm.shape[0] == n4c: + ncol = nmo_pos + + # MO density + elif dm.shape[0] == nmo_pos: + ncol = nmo_pos + + else: + raise ValueError( + "Density matrix dimension does not match AO " + "or positive-energy MO space" + ) + + + result = numpy.zeros((natm, ncol)) + + for i in range(natm): + result[i, :] = Epv_atom(mol, mf, i, dm) + return result diff --git a/pyscf/pv/test/test_energy-pv.py b/pyscf/pv/test/test_energy-pv.py index 62aa3fa3f..f026a9e8b 100644 --- a/pyscf/pv/test/test_energy-pv.py +++ b/pyscf/pv/test/test_energy-pv.py @@ -56,7 +56,15 @@ def tearDownModule(): class KnownValues(unittest.TestCase): def test_pv(self): # Reference results (for PV) compared with DIRAC24 - Epv = Epv_molecule(mol, mf) + Epv = Epv_molecule(mol, mf, dm=None) + Epv_Oxigen = numpy.sum(Epv, axis=1)[0] + self.assertAlmostEqual(Epv_Oxigen, -2.745821e-21, 5) + self.assertAlmostEqual(mf.e_tot,-317.831176378,8) + + def test_pv_dm(self): +# Reference results (for PV) compared with DIRAC24, using the density matrix + density_matrix = mf.make_rdm1() + Epv = Epv_molecule(mol, mf, dm=density_matrix) Epv_Oxigen = numpy.sum(Epv, axis=1)[0] self.assertAlmostEqual(Epv_Oxigen, -2.745821e-21, 5) self.assertAlmostEqual(mf.e_tot,-317.831176378,8) From 1325e9c342c69b0d545f1e39561ef1cecc935259 Mon Sep 17 00:00:00 2001 From: juanjoaucar Date: Sat, 25 Jul 2026 19:21:54 -0300 Subject: [PATCH 2/4] Fix issues in pv --- pyscf/pv/energy.py | 12 ++++-------- 1 file changed, 4 insertions(+), 8 deletions(-) mode change 100755 => 100644 pyscf/pv/energy.py diff --git a/pyscf/pv/energy.py b/pyscf/pv/energy.py old mode 100755 new mode 100644 index dc61f9207..20e778daf --- a/pyscf/pv/energy.py +++ b/pyscf/pv/energy.py @@ -19,8 +19,9 @@ ''' import numpy + +from pyscf import gto, lib from pyscf.scf import dhf -from pyscf import lib, gto def fc_integrals(mol, mf, atom, **kwargs): @@ -97,7 +98,6 @@ def gamma5_fc_integrals(mol, mf, atom, **kwargs): atom_index = 0 fc_matrix = PV_integrals(mol, mf, atom_index) """ - c = lib.param.LIGHT_SPEED n4c = mf.mo_coeff.shape[0] n2c = n4c // 2 @@ -374,12 +374,8 @@ def Epv_molecule(mol, mf, dm=None): nNeg = nmo // 2 nmo_pos = nmo - nNeg - # AO density - if dm.shape[0] == n4c: - ncol = nmo_pos - - # MO density - elif dm.shape[0] == nmo_pos: + # AO or MO density + if dm.shape[0] == n4c or dm.shape[0] == nmo_pos: ncol = nmo_pos else: From de1e459ae2e1f2badd137bfd5ff6c2214508917f Mon Sep 17 00:00:00 2001 From: juanjoaucar Date: Sun, 26 Jul 2026 02:29:25 -0300 Subject: [PATCH 3/4] ruff check on test_energy-pv file --- pyscf/pv/test/test_energy-pv.py | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/pyscf/pv/test/test_energy-pv.py b/pyscf/pv/test/test_energy-pv.py index f026a9e8b..de2b59018 100644 --- a/pyscf/pv/test/test_energy-pv.py +++ b/pyscf/pv/test/test_energy-pv.py @@ -1,4 +1,3 @@ -#!/usr/bin/env python # Copyright 2014-2022 The PySCF Developers. All Rights Reserved. # # Licensed under the Apache License, Version 2.0 (the "License"); @@ -14,13 +13,14 @@ # limitations under the License. # -import numpy import unittest -from pyscf import gto -from pyscf import scf -from pyscf import lib + +import numpy + +from pyscf import gto, lib, scf from pyscf.pv.energy import Epv_molecule + def setUpModule(): global mol, mf From 1069d6b78af41a35e2a98b1352329b0dc008df27 Mon Sep 17 00:00:00 2001 From: juanjoaucar Date: Sun, 26 Jul 2026 17:17:53 -0300 Subject: [PATCH 4/4] differentiate according to the existence of dm --- pyscf/pv/energy.py | 52 +++++++++++++++++++++------------------------- 1 file changed, 24 insertions(+), 28 deletions(-) diff --git a/pyscf/pv/energy.py b/pyscf/pv/energy.py index 20e778daf..5c4282cee 100644 --- a/pyscf/pv/energy.py +++ b/pyscf/pv/energy.py @@ -236,35 +236,31 @@ def Epv_atom(mol, mf, atom_index, dm=None): # Prefactor fac = (2.2225 * 10 ** (-14) * QW) / (2 * numpy.sqrt(2)) - # ------------------------------------------------------------ - # PV operator in AO basis - # ------------------------------------------------------------ - - # Get the Fermi Contact integrals for the selected atom - fc_ao = fc_integrals(mol, mf, atom_index) + if dm is None: + # ------------------------------------------------------------ + # PV operator in AO basis + # ------------------------------------------------------------ - # Get the Gamma5 Fermi-Contact integrals for the selected atom (in AO basis) - gamma5_fc_ao = gamma5_fc_integrals(mol, mf, atom_index) + # Get the Fermi Contact integrals for the selected atom + fc_ao = fc_integrals(mol, mf, atom_index) - # ------------------------------------------------------------ - # Orbital expectation values (4c-DHF or 4c-DFT reference) - # ------------------------------------------------------------ + # ------------------------------------------------------------ + # Orbital expectation values (4c-DHF or 4c-DFT reference) + # ------------------------------------------------------------ - # Get the positive components of the MO spinors - Lo = mf.mo_coeff[:n2c, nNeg:nNeg + nocc] - So = mf.mo_coeff[n2c:, nNeg:nNeg + nocc] + # Get the positive components of the MO spinors + Lo = mf.mo_coeff[:n2c, nNeg:nNeg + nocc] + So = mf.mo_coeff[n2c:, nNeg:nNeg + nocc] - # Expectation values - expval_LS = numpy.einsum('ij,ji->i', Lo.conjugate().T @ fc_ao[:n2c, n2c:], So) - expval_SL = numpy.einsum('ij,ji->i', So.conjugate().T @ fc_ao[n2c:, :n2c], Lo) - expval_perorb = expval_LS + expval_SL + # Expectation values + expval_LS = numpy.einsum('ij,ji->i', Lo.conjugate().T @ fc_ao[:n2c, n2c:], So) + expval_SL = numpy.einsum('ij,ji->i', So.conjugate().T @ fc_ao[n2c:, :n2c], Lo) + expval_perorb = expval_LS + expval_SL # ------------------------------------------------------------ # Density matrix contraction # ------------------------------------------------------------ - expval_dm = None - if dm is not None: # Positive-energy MO coefficient matrix @@ -272,6 +268,9 @@ def Epv_atom(mol, mf, atom_index, dm=None): nmo_pos = C.shape[1] + # Get the Gamma5 Fermi-Contact integrals for the selected atom (in AO basis) + gamma5_fc_ao = gamma5_fc_integrals(mol, mf, atom_index) + # Transform operator AO -> MO # Operator transformation: # @@ -309,10 +308,12 @@ def Epv_atom(mol, mf, atom_index, dm=None): # Total expectation value # E = Tr(𝑫 h) - expval_dm = numpy.trace(dm_mo @ gamma5_fc_mo) + #expval_dm = numpy.trace(dm_mo @ gamma5_fc_mo) # Orbital decomposition in MO basis - expval_dm_perorb = numpy.einsum( + # For correlated densities (MP2/CCSD), the orbital decomposition + # is representation-dependent. Only the total trace is invariant. + expval_perorb = numpy.einsum( 'ij,ji->i', gamma5_fc_mo, dm_mo @@ -323,12 +324,7 @@ def Epv_atom(mol, mf, atom_index, dm=None): # Return # ------------------------------------------------------------ - if expval_dm is None: - return fac * expval_perorb.real - - else: - return fac * expval_dm_perorb.real - + return fac * expval_perorb.real