diff --git a/pyscf/pv/energy.py b/pyscf/pv/energy.py old mode 100755 new mode 100644 index 3b32b00f3..5c4282cee --- 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): @@ -72,6 +73,50 @@ 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) + """ + 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,111 @@ 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] - # Get the Fermi Contact integrals for the selected atom - fc_ao = fc_integrals(mol, mf, atom_index) - - # 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)) + if dm is None: + # ------------------------------------------------------------ + # PV operator in AO basis + # ------------------------------------------------------------ + + # 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) + # ------------------------------------------------------------ + + # 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 + + # ------------------------------------------------------------ + # Density matrix contraction + # ------------------------------------------------------------ + + if dm is not None: + + # Positive-energy MO coefficient matrix + C = mf.mo_coeff[:, nNeg:] + + 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: + # + # 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 + # 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 + ) + + + # ------------------------------------------------------------ + # Return + # ------------------------------------------------------------ + return fac * expval_perorb.real -def Epv_molecule(mol, mf): + + + + +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 +340,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 +357,34 @@ 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 or MO density + if dm.shape[0] == n4c or 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..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 @@ -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)