From fd360c59ba07b154721854cb6445d5332487c4f4 Mon Sep 17 00:00:00 2001 From: James Brown Date: Tue, 24 May 2022 15:48:58 -0400 Subject: [PATCH 1/6] changed default basis, added ecp option --- .../toolboxes/molecular_computation/molecule.py | 16 ++++++++++------ 1 file changed, 10 insertions(+), 6 deletions(-) diff --git a/tangelo/toolboxes/molecular_computation/molecule.py b/tangelo/toolboxes/molecular_computation/molecule.py index 1265f7329..33bc6971f 100644 --- a/tangelo/toolboxes/molecular_computation/molecule.py +++ b/tangelo/toolboxes/molecular_computation/molecule.py @@ -94,7 +94,7 @@ class Molecule: n_min_orbitals: int = field(init=False) def __post_init__(self): - mol = self.to_pyscf(basis="sto-3g", symmetry=False) + mol = self.to_pyscf() self.n_atoms = mol.natm self.n_electrons = mol.nelectron self.n_min_orbitals = mol.nao_nr() @@ -107,7 +107,7 @@ def elements(self): def coords(self): return np.array([self.xyz[i][1] for i in range(self.n_atoms)]) - def to_pyscf(self, basis="sto-3g", symmetry=False): + def to_pyscf(self, basis="CRENBL", symmetry=False, ecp=False): """Method to return a pyscf.gto.Mole object. Args: @@ -122,6 +122,7 @@ def to_pyscf(self, basis="sto-3g", symmetry=False): mol.charge = self.q mol.spin = self.spin mol.symmetry = symmetry + mol.ecp = ecp mol.build() self.xyz = list() @@ -137,7 +138,7 @@ def to_file(self, filename, format=None): filename (str): The name of the file to output the geometry. format (str): The output type of "raw", "xyz", or "zmat". If None, will be inferred by the filename """ - mol = self.to_pyscf(basis="sto-3g") + mol = self.to_pyscf() mol.tofile(filename, format) def to_openfermion(self, basis="sto-3g"): @@ -161,6 +162,8 @@ class SecondQuantizedMolecule(Molecule): Attributes: basis (string): Basis set. + ecp (dict): The effective core potential (ecp) for any atoms in the molecule. + e.g. {"C": "crenbl"} use CRENBL ecp for Carbon atoms. symmetry (bool or str): Whether to use symmetry in RHF or ROHF calculation. Can also specify point group using pyscf allowed string. e.g. "Dooh", "D2h", "C2v", ... @@ -195,6 +198,7 @@ class SecondQuantizedMolecule(Molecule): actives_mos (list): Active MOs indexes. """ basis: str = "sto-3g" + ecp: bool = False symmetry: bool = False frozen_orbitals: list or int = field(default="frozen_core", repr=False) @@ -273,7 +277,7 @@ def _compute_mean_field(self): (mf_energy, mo_energies, mo_occ, n_mos and n_sos). """ - molecule = self.to_pyscf(self.basis, self.symmetry) + molecule = self.to_pyscf(self.basis, self.symmetry, self.ecp) self.mean_field = scf.RHF(molecule) self.mean_field.verbose = 0 @@ -331,7 +335,7 @@ def _convert_frozen_orbitals(self, frozen_orbitals): virtual and frozen virtual orbital indexes. """ - if frozen_orbitals == "frozen_core": + if frozen_orbitals == "frozen_core" and not self.ecp: frozen_orbitals = get_frozen_core(self.to_pyscf(self.basis)) elif frozen_orbitals is None: frozen_orbitals = 0 @@ -456,7 +460,7 @@ def get_integrals(self, mo_coeff=None, consider_frozen=True): """ # Pyscf molecule to get integrals. - pyscf_mol = self.to_pyscf(self.basis, self.symmetry) + pyscf_mol = self.to_pyscf(self.basis, self.symmetry, self.ecp) if mo_coeff is None: mo_coeff = self.mean_field.mo_coeff From 1cc6a7e157acb0b3fa21c4d7c16b54c5715df258 Mon Sep 17 00:00:00 2001 From: James Brown Date: Wed, 25 May 2022 12:50:44 -0400 Subject: [PATCH 2/6] response to PR --- tangelo/toolboxes/molecular_computation/frozen_orbitals.py | 2 +- tangelo/toolboxes/molecular_computation/molecule.py | 6 ++++-- 2 files changed, 5 insertions(+), 3 deletions(-) diff --git a/tangelo/toolboxes/molecular_computation/frozen_orbitals.py b/tangelo/toolboxes/molecular_computation/frozen_orbitals.py index dc81038e9..6bf5cecb9 100644 --- a/tangelo/toolboxes/molecular_computation/frozen_orbitals.py +++ b/tangelo/toolboxes/molecular_computation/frozen_orbitals.py @@ -42,7 +42,7 @@ def get_frozen_core(molecule): # Counting how many of each element is in the molecule. elements = {i: molecule.elements.count(i) for i in molecule.elements} - frozen_core = sum([v * core_orbitals[k] for k, v in elements.items()]) + frozen_core = sum([v * core_orbitals.get(k, 0) for k, v in elements.items()]) return frozen_core diff --git a/tangelo/toolboxes/molecular_computation/molecule.py b/tangelo/toolboxes/molecular_computation/molecule.py index 33bc6971f..1b44abadb 100644 --- a/tangelo/toolboxes/molecular_computation/molecule.py +++ b/tangelo/toolboxes/molecular_computation/molecule.py @@ -107,11 +107,13 @@ def elements(self): def coords(self): return np.array([self.xyz[i][1] for i in range(self.n_atoms)]) - def to_pyscf(self, basis="CRENBL", symmetry=False, ecp=False): + def to_pyscf(self, basis="CRENBL", symmetry=False, ecp=dict()): """Method to return a pyscf.gto.Mole object. Args: basis (string): Basis set. + symmetry (bool): Flag to turn symmetry on + ecp (dict): Dictionary with ecp definition for each atom e.g. {"Cu": "crenbl"} Returns: pyscf.gto.Mole: PySCF compatible object. @@ -198,7 +200,7 @@ class SecondQuantizedMolecule(Molecule): actives_mos (list): Active MOs indexes. """ basis: str = "sto-3g" - ecp: bool = False + ecp: dict = field(default_factory=dict) symmetry: bool = False frozen_orbitals: list or int = field(default="frozen_core", repr=False) From 54e218112bcdcf1b6c164f234848079adeb57c42 Mon Sep 17 00:00:00 2001 From: James Brown Date: Wed, 25 May 2022 12:58:24 -0400 Subject: [PATCH 3/6] change default ecp to immutable None --- tangelo/toolboxes/molecular_computation/molecule.py | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/tangelo/toolboxes/molecular_computation/molecule.py b/tangelo/toolboxes/molecular_computation/molecule.py index 1b44abadb..97979cd4b 100644 --- a/tangelo/toolboxes/molecular_computation/molecule.py +++ b/tangelo/toolboxes/molecular_computation/molecule.py @@ -107,7 +107,7 @@ def elements(self): def coords(self): return np.array([self.xyz[i][1] for i in range(self.n_atoms)]) - def to_pyscf(self, basis="CRENBL", symmetry=False, ecp=dict()): + def to_pyscf(self, basis="CRENBL", symmetry=False, ecp=None): """Method to return a pyscf.gto.Mole object. Args: @@ -118,6 +118,8 @@ def to_pyscf(self, basis="CRENBL", symmetry=False, ecp=dict()): Returns: pyscf.gto.Mole: PySCF compatible object. """ + if ecp is None: + ecp = dict() mol = gto.Mole(atom=self.xyz) mol.basis = basis From ed5226c7f3ad00aa733d3d8db2b77ee7806eba7d Mon Sep 17 00:00:00 2001 From: James Brown Date: Thu, 26 May 2022 09:43:04 -0400 Subject: [PATCH 4/6] added test --- .../toolboxes/molecular_computation/molecule.py | 4 ++-- .../molecular_computation/tests/test_molecule.py | 15 +++++++++++++++ 2 files changed, 17 insertions(+), 2 deletions(-) diff --git a/tangelo/toolboxes/molecular_computation/molecule.py b/tangelo/toolboxes/molecular_computation/molecule.py index 97979cd4b..8cb7629f5 100644 --- a/tangelo/toolboxes/molecular_computation/molecule.py +++ b/tangelo/toolboxes/molecular_computation/molecule.py @@ -339,8 +339,8 @@ def _convert_frozen_orbitals(self, frozen_orbitals): virtual and frozen virtual orbital indexes. """ - if frozen_orbitals == "frozen_core" and not self.ecp: - frozen_orbitals = get_frozen_core(self.to_pyscf(self.basis)) + if frozen_orbitals == "frozen_core": + frozen_orbitals = get_frozen_core(self.to_pyscf(self.basis)) if not self.ecp else 0 elif frozen_orbitals is None: frozen_orbitals = 0 diff --git a/tangelo/toolboxes/molecular_computation/tests/test_molecule.py b/tangelo/toolboxes/molecular_computation/tests/test_molecule.py index 347b55368..a2f99bf3a 100644 --- a/tangelo/toolboxes/molecular_computation/tests/test_molecule.py +++ b/tangelo/toolboxes/molecular_computation/tests/test_molecule.py @@ -184,6 +184,21 @@ def test_symmetry_label(self): assert(mo_symm_labels == molecule.mo_symm_labels) assert(mo_symm_ids == molecule.mo_symm_ids) + def test_ecp(self): + """Verify that the number of electrons is reduced when ecp is called.""" + + molecule = SecondQuantizedMolecule(xyz="Yb", q=0, spin=0, basis="crenbl", ecp="crenbl") + # "Yb" has 70 electrons but the ecp reduces this to 16 + assert(molecule.n_active_electrons == 16) + assert(molecule.n_active_mos == 96) + + molecule = SecondQuantizedMolecule(xyz="Cu", q=0, spin=1, basis="cc-pvdz", ecp="crenbl", + frozen_orbitals=list(range(8))) + # "Cu" has 29 electrons but the ecp reduces this to 19 + assert(molecule.n_electrons == 19) + assert(molecule.n_active_electrons == 3) + assert(molecule.n_active_mos == 35) + if __name__ == "__main__": unittest.main() From 88c5b1d8e874db16a4efc38e6e76881912822951 Mon Sep 17 00:00:00 2001 From: James Brown Date: Thu, 26 May 2022 09:46:30 -0400 Subject: [PATCH 5/6] changed initialization of ecp --- tangelo/toolboxes/molecular_computation/molecule.py | 4 +--- 1 file changed, 1 insertion(+), 3 deletions(-) diff --git a/tangelo/toolboxes/molecular_computation/molecule.py b/tangelo/toolboxes/molecular_computation/molecule.py index 8cb7629f5..7306e1de1 100644 --- a/tangelo/toolboxes/molecular_computation/molecule.py +++ b/tangelo/toolboxes/molecular_computation/molecule.py @@ -118,15 +118,13 @@ def to_pyscf(self, basis="CRENBL", symmetry=False, ecp=None): Returns: pyscf.gto.Mole: PySCF compatible object. """ - if ecp is None: - ecp = dict() mol = gto.Mole(atom=self.xyz) mol.basis = basis mol.charge = self.q mol.spin = self.spin mol.symmetry = symmetry - mol.ecp = ecp + mol.ecp = ecp if ecp else dict() mol.build() self.xyz = list() From 1f0597385b97898b5a8849bfb0b201453d56c994 Mon Sep 17 00:00:00 2001 From: James Brown Date: Thu, 26 May 2022 10:12:24 -0400 Subject: [PATCH 6/6] fixed test --- tangelo/toolboxes/molecular_computation/tests/test_molecule.py | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/tangelo/toolboxes/molecular_computation/tests/test_molecule.py b/tangelo/toolboxes/molecular_computation/tests/test_molecule.py index a2f99bf3a..cf0a8af12 100644 --- a/tangelo/toolboxes/molecular_computation/tests/test_molecule.py +++ b/tangelo/toolboxes/molecular_computation/tests/test_molecule.py @@ -194,8 +194,7 @@ def test_ecp(self): molecule = SecondQuantizedMolecule(xyz="Cu", q=0, spin=1, basis="cc-pvdz", ecp="crenbl", frozen_orbitals=list(range(8))) - # "Cu" has 29 electrons but the ecp reduces this to 19 - assert(molecule.n_electrons == 19) + # "Cu" has 29 electrons but the ecp reduces this to 19. The active electrons are 19 - 8 * 2 = 3 assert(molecule.n_active_electrons == 3) assert(molecule.n_active_mos == 35)