From 016293ab9ec80d46bef1e9e6668b89628d3bf526 Mon Sep 17 00:00:00 2001 From: juanjoaucar Date: Mon, 11 Nov 2024 12:21:54 -0300 Subject: [PATCH 1/4] Parity violating contributions to molecular energy. --- CHANGELOG | 6 + examples/pv/00.energy-PV.py | 32 +++++ pyscf/pv/energy.py | 213 ++++++++++++++++++++++++++++++++ pyscf/pv/test/test_energy-pv.py | 70 +++++++++++ 4 files changed, 321 insertions(+) create mode 100644 examples/pv/00.energy-PV.py create mode 100755 pyscf/pv/energy.py create mode 100644 pyscf/pv/test/test_energy-pv.py diff --git a/CHANGELOG b/CHANGELOG index 2835c98..adb6fcb 100644 --- a/CHANGELOG +++ b/CHANGELOG @@ -1,3 +1,9 @@ +PySCF-Forge 1.0.2 (2024-11-11) +------------------------------ +* New features + - Parity Violating contributions to energy + + PySCF-Forge 1.0.1 (2024-10-31) ------------------------------ * New features diff --git a/examples/pv/00.energy-PV.py b/examples/pv/00.energy-PV.py new file mode 100644 index 0000000..506b41a --- /dev/null +++ b/examples/pv/00.energy-PV.py @@ -0,0 +1,32 @@ +#!/usr/bin/env python +# +# Author: Juan J. Aucar +# +# +''' +Parity-violating contribution to the molecular energy + +Valid for the point-nuclear model within the four-component relativistic framework. +The total contribution correspond to the first term of Eq. (4) in https://doi.org/10.1002/wcms.1396. +The sum over the nuclei of the system should be done as described in Eq. (4) of https://doi.org/10.1021/acs.jpclett.3c03038 +''' + +import numpy +from pyscf import gto, scf +from pyscf.pv.energy import Epv_molecule + +#Alanine molecule +mol = gto.M( + atom = 'O -0.000008 0.000006 0.473161; O -0.498429 1.617953 -0.942950; N -2.916494 2.018558 0.304530; C -2.245961 0.738717 0.446378; C -2.933825 -0.437589 -0.265779; C -0.836260 0.869228 -0.089564; H -2.164332 0.502686 1.502658; H -2.396710 -1.368611 -0.107150; H -3.940684 -0.559206 0.124631; H -3.002345 -0.251817 -1.334665; H -3.903065 1.915512 0.431392; H -2.755346 2.398021 -0.608200; H 0.850292 0.091697 0.064913', + basis = 'sto3g', + symmetry = False, + verbose = 3 +) + +mf = scf.DHF(mol) +mf.conv_tol = 1e-9 +mf.with_ssss=False +mf.kernel() + +Epv = Epv_molecule(mol, mf) +print('PV contribution to energy from the first nucleus (O) of the system: %.15g, ref = -2.745821e-21' % numpy.sum(Epv, axis=1)[0]) diff --git a/pyscf/pv/energy.py b/pyscf/pv/energy.py new file mode 100755 index 0000000..3b32b00 --- /dev/null +++ b/pyscf/pv/energy.py @@ -0,0 +1,213 @@ +# Copyright 2014-2024 The PySCF Developers. All Rights Reserved. +# +# Licensed under the Apache License, Version 2.0 (the "License"); +# you may not use this file except in compliance with the License. +# You may obtain a copy of the License at +# +# http://www.apache.org/licenses/LICENSE-2.0 +# +# Unless required by applicable law or agreed to in writing, software +# distributed under the License is distributed on an "AS IS" BASIS, +# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +# See the License for the specific language governing permissions and +# limitations under the License. + +''' + +@author Juan Jose Aucar + +''' + +import numpy +from pyscf.scf import dhf +from pyscf import lib, gto + + +def fc_integrals(mol, mf, atom, **kwargs): + """ + Compute the Fermi Contact (FC) integrals for a given atom. + + This function calculates the Fermi Contact integrals in + relativistic form for a specified atom in the molecule using + four-component Dirac-Coulomb spinor wavefunctions. + + The FC integrals consist of four distinct blocks: + - LL: Large-Large component + - LS: Large-Small component + - SL: Small-Large component + - SS: Small-Small component + + :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 FC integrals are computed. + :type atom: int + :return: Matrix of Fermi Contact 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 = fc_integrals(mol, mf, atom_index) + """ + c = lib.param.LIGHT_SPEED + n4c = mf.mo_coeff.shape[0] + n2c = n4c // 2 + coordinates = mf.mol.atom_coords()[[atom]] + + # Obtaining AO integrals in spinor form + ao_spinor = gto.eval_gto(mf.mol, "GTOval_spinor", coordinates, comp=1)[:, 0, :] + ao_spinor_S = gto.eval_gto(mf.mol, "GTOval_sp_spinor", coordinates, comp=1)[:, 0, :] + + # Forming the FC integrals matrix (LL, LS, SL, SS blocks) + fc_integrals = numpy.zeros((n4c, n4c), dtype=numpy.complex128) + fc_integrals[:n2c, :n2c] = numpy.einsum("ip,iq->pq", ao_spinor.conjugate(), ao_spinor) + fc_integrals[:n2c, n2c:] = numpy.einsum("ip,iq->pq", ao_spinor.conjugate(), ao_spinor_S) * (0.5 / c) + fc_integrals[n2c:, :n2c] = numpy.einsum("ip,iq->pq", ao_spinor_S.conjugate(), ao_spinor) * (0.5 / c) + fc_integrals[n2c:, n2c:] = numpy.einsum("ip,iq->pq", ao_spinor_S.conjugate(), ao_spinor_S) * ((0.5 / c) ** 2) + + return fc_integrals + +def fc_expval(mol, mf, atom): + """ + Calculate the Fermi Contact (FC) expectation values for each occupied orbital in a molecule, + focusing on a specific atom. + + The expectation values are calculated using the four-component Dirac-Coulomb spinor + wavefunctions. The function computes the contributions from the large-large (LL) and + small-small (SS) components of the spinor for each occupied orbital. + + :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 FC expectation values are computed. + :type atom: int + :return: Array of FC expectation values for each occupied orbital. + :rtype: numpy.ndarray (real), shape (nocc,) + + Example: + mol = gto.Mole() + # Setup mol with geometry and basis set... + mf = scf.DHF(mol) + mf.kernel() + atom_index = 0 + fc_expval_per_orbital = fc_expval(mol, mf, atom_index) + """ + n4c, nmo = mf.mo_coeff.shape + n2c = n4c // 2 + nocc = mf.mol.nelectron + expval_perorb = numpy.zeros(nocc) + + mo_pos_l = mf.mo_coeff[:n2c, nmo//2:] + mo_pos_s = mf.mo_coeff[n2c:, nmo//2:] + + Lo = mo_pos_l[:, :nocc] + So = mo_pos_s[:, :nocc] + + fac = 8 * numpy.pi / 3 + fc_ao = fc_integrals(mol, mf, atom) + + # Split the fc_integrals matrix into LL and SS blocks + fc_ao_LL = fc_ao[:n2c, :n2c] + 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_perorb[k] = expval_LL_k + expval_SS_k + + return fac * expval_perorb.real + + + +def Epv_atom(mol, mf, atom_index): + """ + Calculate the parity-violating (PV) contribution to energy for a given atom in a molecule. + + It then returns the contributions from each occupied orbital. + (First term of Eq. 4 - https://doi.org/10.1002/wcms.1396) + + :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: Index of the atom for which the PV contribution is computed. + :type atom_index: int + :return: Array of PV expectation values for each occupied orbital. + :rtype: numpy.ndarray (real), shape (nocc,) + + Example: + mol = gto.Mole() + # Setup mol with geometry and basis set... + mf = scf.DHF(mol) + mf.kernel() + atom_index = 0 + pv_values = 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 + + + # Get the positive components of the MO spinors + Lo = mf.mo_coeff[:n2c, nNeg:nNeg + nocc] + So = mf.mo_coeff[n2c:, nNeg:nNeg + nocc] + + # Atom masses and atomic numbers + masses = mf.mol.atom_mass_list(isotope_avg=False) + atomic_numbers = mf.mol.atom_charges() + + # Neutrons per atom + neutrons = masses - atomic_numbers + + S2THETAW = 0.23122 # AS DIRAC24 (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)) + + return fac * expval_perorb.real + + + +def Epv_molecule(mol, mf): + """ + Calculate the weak charge parity-violating (PV) contributions for all atoms in the molecule + within a punctual nuclear charge distribution model. + + This function iterates over all atoms in the molecule and computes the PV contributions + for each atom. + + :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. + :type mf: pyscf.scf.hf.SCF + :return: A 2D array where each row corresponds to an atom and each column contains the PV + contributions for the occupied orbitals of that atom. + :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) + return result + diff --git a/pyscf/pv/test/test_energy-pv.py b/pyscf/pv/test/test_energy-pv.py new file mode 100644 index 0000000..6304bb3 --- /dev/null +++ b/pyscf/pv/test/test_energy-pv.py @@ -0,0 +1,70 @@ +#!/usr/bin/env python +# Copyright 2014-2022 The PySCF Developers. All Rights Reserved. +# +# Licensed under the Apache License, Version 2.0 (the "License"); +# you may not use this file except in compliance with the License. +# You may obtain a copy of the License at +# +# http://www.apache.org/licenses/LICENSE-2.0 +# +# Unless required by applicable law or agreed to in writing, software +# distributed under the License is distributed on an "AS IS" BASIS, +# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +# See the License for the specific language governing permissions and +# limitations under the License. +# + +import numpy +import unittest +from pyscf import gto +from pyscf import scf +from pyscf import lib +from pyscf.pv.energy import Epv_molecule + +def setUpModule(): + global mol, mf + +# alanine molecule + mol = gto.M( + atom = ''' + O -0.000008 0.000006 0.473161 + O -0.498429 1.617953 -0.942950 + N -2.916494 2.018558 0.304530 + C -2.245961 0.738717 0.446378 + C -2.933825 -0.437589 -0.265779 + C -0.836260 0.869228 -0.089564 + H -2.164332 0.502686 1.502658 + H -2.396710 -1.368611 -0.107150 + H -3.940684 -0.559206 0.124631 + H -3.002345 -0.251817 -1.334665 + H -3.903065 1.915512 0.431392 + H -2.755346 2.398021 -0.608200 + H 0.850292 0.091697 0.064913 + ''', + basis = 'sto3g', + verbose = 3) + + mf = scf.DHF(mol_alanine) + mf.conv_tol = 1e-9 + mf.with_ssss=False + mf.kernel() + +def tearDownModule(): + global mf + mol.stdout.close() + del mol, mf + +class KnownValues(unittest.TestCase): + def test_pv(self): +# Reference results (for PV) are obtained from DIRAC24 + Epv = Epv_molecule(mol, mf) + Epv_Oxigen = numpy.sum(Epv, axis=1)[0] + self.assertAlmostEqual(Epv_Oxigen, -2.745821e-21, 5) + self.assertAlmostEqual(mf.e_tot,-317.831176378,8) + + +if __name__ == "__main__": + print("Full Tests for Parity violating quantities") + unittest.main() + + From c5d813f936de0ef9d266044fac5ebde24f9b5a47 Mon Sep 17 00:00:00 2001 From: juanjoaucar Date: Mon, 11 Nov 2024 13:19:12 -0300 Subject: [PATCH 2/4] Fix pv test --- pyscf/pv/test/test_energy-pv.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/pyscf/pv/test/test_energy-pv.py b/pyscf/pv/test/test_energy-pv.py index 6304bb3..5889559 100644 --- a/pyscf/pv/test/test_energy-pv.py +++ b/pyscf/pv/test/test_energy-pv.py @@ -44,7 +44,7 @@ def setUpModule(): basis = 'sto3g', verbose = 3) - mf = scf.DHF(mol_alanine) + mf = scf.DHF(mol) mf.conv_tol = 1e-9 mf.with_ssss=False mf.kernel() From 85436f97fddf07097eabdd950eba4eff62bb66d1 Mon Sep 17 00:00:00 2001 From: juanjoaucar Date: Mon, 11 Nov 2024 13:38:57 -0300 Subject: [PATCH 3/4] Small fix. --- pyscf/pv/test/test_energy-pv.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/pyscf/pv/test/test_energy-pv.py b/pyscf/pv/test/test_energy-pv.py index 5889559..437b1c8 100644 --- a/pyscf/pv/test/test_energy-pv.py +++ b/pyscf/pv/test/test_energy-pv.py @@ -50,13 +50,13 @@ def setUpModule(): mf.kernel() def tearDownModule(): - global mf + global mol, mf mol.stdout.close() del mol, mf class KnownValues(unittest.TestCase): def test_pv(self): -# Reference results (for PV) are obtained from DIRAC24 +# Reference results (for PV) compared with DIRAC24 Epv = Epv_molecule(mol, mf) Epv_Oxigen = numpy.sum(Epv, axis=1)[0] self.assertAlmostEqual(Epv_Oxigen, -2.745821e-21, 5) From 96287fd74eec402da9670849a1027a38039ed9c5 Mon Sep 17 00:00:00 2001 From: juanjoaucar Date: Mon, 11 Nov 2024 14:11:40 -0300 Subject: [PATCH 4/4] Fix PV test --- pyscf/pv/test/test_energy-pv.py | 1 - 1 file changed, 1 deletion(-) diff --git a/pyscf/pv/test/test_energy-pv.py b/pyscf/pv/test/test_energy-pv.py index 437b1c8..62aa3fa 100644 --- a/pyscf/pv/test/test_energy-pv.py +++ b/pyscf/pv/test/test_energy-pv.py @@ -51,7 +51,6 @@ def setUpModule(): def tearDownModule(): global mol, mf - mol.stdout.close() del mol, mf class KnownValues(unittest.TestCase):