Change Molecule.nelec to the total number of electrons

examples/br_example.py

@@ -42,7 +42,7 @@
 hamiltonian = mol.hamiltonian(oei=oei, tei=tei)
 
 # Check that the hamiltonian with a HF reference ansatz doesn't yield the correct HF energy
-circuit = hf_circuit(mol.nso, sum(mol.nelec))
+circuit = hf_circuit(mol.nso, mol.nelec)
 print(
     f"\nElectronic energy: {expectation(circuit, hamiltonian):.8f} (From the H_core guess, should be > actual HF energy)"
 )
@@ -52,14 +52,14 @@ def electronic_energy(parameters):
     """
     Loss function (Electronic energy) for the basis rotation ansatz
     """
-    circuit = hf_circuit(mol.nso, sum(mol.nelec))
-    circuit += br_circuit(mol.nso, parameters, sum(mol.nelec))
+    circuit = hf_circuit(mol.nso, mol.nelec)
+    circuit += br_circuit(mol.nso, parameters, mol.nelec)
 
     return expectation(circuit, hamiltonian)
 
 
 # Build  a basis rotation_circuit
-params = np.random.rand(sum(mol.nelec) * (mol.nso - sum(mol.nelec)))  # n_occ * n_virt
+params = np.random.rand(mol.nelec * (mol.nso - mol.nelec))  # n_occ * n_virt
 
 best, params, extra = optimize(electronic_energy, params)
 
@@ -70,7 +70,7 @@ def electronic_energy(parameters):
 # print("Optimized parameters:", params)
 
 
-params = np.random.rand(sum(mol.nelec) * (mol.nso - sum(mol.nelec)))  # n_occ * n_virt
+params = np.random.rand(mol.nelec * (mol.nso - mol.nelec))  # n_occ * n_virt
 
 result = minimize(electronic_energy, params)
 best, params = result.fun, result.x

src/qibochem/driver/molecule.py

@@ -49,18 +49,19 @@ def __init__(self, geometry=None, charge=0, multiplicity=1, basis=None, xyz_file
         else:
             self.basis = basis
 
+        self.nelec = None  #: Total number of electrons for the molecule
+        self.norb = None  #: Number of molecular orbitals considered for the molecule
+        self.nso = None  #: Number of molecular spin-orbitals considered for the molecule
+        self.e_hf = None  #: Hartree-Fock energy
+        self.oei = None  #: One-electron integrals
+        self.tei = None  #: Two-electron integrals, order follows the second quantization notation
+
         self.ca = None
         self.pa = None
         self.da = None
-        self.nelec = None
         self.nalpha = None
         self.nbeta = None
-        self.oei = None
-        self.tei = None
-        self.e_hf = None
         self.e_nuc = None
-        self.norb = None
-        self.nso = None
         self.overlap = None
         self.eps = None
         self.fa = None
@@ -70,15 +71,15 @@ def __init__(self, geometry=None, charge=0, multiplicity=1, basis=None, xyz_file
         self.aoeri = None
 
         # For HF embedding
-        self.active = active  # List of active MOs included in the active space
+        self.active = active  #: Iterable of molecular orbitals included in the active space
         self.frozen = None
 
         self.inactive_energy = None
         self.embed_oei = None
         self.embed_tei = None
 
-        self.n_active_e = None
-        self.n_active_orbs = None  # Number of spin-orbitals in the active space
+        self.n_active_e = None  #: Number of electrons included in the active space if HF embedding is used
+        self.n_active_orbs = None  #: Number of spin-orbitals in the active space if HF embedding is used
 
     def _process_xyz_file(self, xyz_file, charge, multiplicity):
         """
@@ -136,9 +137,9 @@ def run_pyscf(self, max_scf_cycles=50):
 
         # Save results from HF calculation
         self.ca = np.asarray(pyscf_job.mo_coeff)  # MO coeffcients
-        self.nelec = pyscf_mol.nelec
-        self.nalpha = self.nelec[0]
-        self.nbeta = self.nelec[1]
+        self.nalpha = pyscf_mol.nelec[0]
+        self.nbeta = pyscf_mol.nelec[1]
+        self.nelec = sum(pyscf_mol.nelec)
         self.e_hf = pyscf_job.e_tot  # HF energy
         self.e_nuc = pyscf_mol.energy_nuc()
         self.norb = self.ca.shape[1]
@@ -220,9 +221,9 @@ def run_pyscf(self, max_scf_cycles=50):
     #     tei = np.einsum("pqrs->prsq", tei)
 
     #     # Fill in the class attributes
-    #     self.nelec = (wavefn.nalpha(), wavefn.nbeta())
-    #     self.nalpha = self.nelec[0]
-    #     self.nbeta = self.nelec[1]
+    #     self.nelec = sum(wavefn.nalpha(), wavefn.nbeta())
+    #     self.nalpha = wavefn.nalpha()
+    #     self.nbeta = wavefn.nbeta()
     #     self.e_hf = ehf
     #     self.e_nuc = psi4_mol.nuclear_repulsion_energy()
     #     self.norb = wavefn.nmo()
@@ -288,7 +289,7 @@ def hf_embedding(self, active=None, frozen=None):
         assert max(active) < self.norb and min(active) >= 0, "Active space must be between 0 " "and the number of MOs"
         if frozen:
             assert not (set(active) & set(frozen)), "Active and frozen space cannot overlap"
-            assert max(frozen) + 1 < sum(self.nelec) // 2 and min(frozen) >= 0, (
+            assert max(frozen) + 1 < self.nelec // 2 and min(frozen) >= 0, (
                 "Frozen orbitals must" " be occupied orbitals"
             )
 
@@ -309,7 +310,7 @@ def hf_embedding(self, active=None, frozen=None):
         self.active = active
         self.frozen = frozen
         self.n_active_orbs = 2 * len(active)
-        self.n_active_e = sum(self.nelec) - 2 * len(self.frozen)
+        self.n_active_e = self.nelec - 2 * len(self.frozen)
 
     def hamiltonian(
         self,
@@ -364,7 +365,7 @@ def hamiltonian(
         if ham_type in ("s", "sym"):
             # Qibo SymbolicHamiltonian
             return symbolic_hamiltonian(ham)
-        # raise NameError(f"Unknown {ham_type}!")  # Shouldn't ever reach here
+        raise NameError(f"Unknown {ham_type}!")  # Shouldn't ever reach here
 
     @staticmethod
     def eigenvalues(hamiltonian):

tests/test_basis_rotation.py

@@ -82,12 +82,12 @@ def test_br_ansatz():
 
     # Define quantum circuit
     circuit = Circuit(mol.nso)
-    circuit.add(gates.X(_i) for _i in range(sum(mol.nelec)))
+    circuit.add(gates.X(_i) for _i in range(mol.nelec))
 
     # Add basis rotation ansatz
     # Initialize all at zero
-    parameters = np.zeros(sum(mol.nelec) * (mol.nso - sum(mol.nelec)))  # n_occ * n_virt
-    circuit += basis_rotation.br_circuit(mol.nso, parameters, sum(mol.nelec))
+    parameters = np.zeros(mol.nelec * (mol.nso - mol.nelec))  # n_occ * n_virt
+    circuit += basis_rotation.br_circuit(mol.nso, parameters, mol.nelec)
 
     def electronic_energy(parameters):
         """

tests/test_hardware_efficient.py

@@ -21,7 +21,7 @@ def test_hea_ansatz():
 
     hea_ansatz = hardware_efficient.hea(nlayers, nqubits)
     qc = models.Circuit(nqubits)
-    qc.add(gates.X(_i) for _i in range(sum(mol.nelec)))
+    qc.add(gates.X(_i) for _i in range(mol.nelec))
     qc.add(hea_ansatz)
     qc.set_parameters(theta)
 
@@ -46,7 +46,7 @@ def test_vqe_hea_ansatz_cost(parameters, circuit, hamiltonian):
 
     hea_ansatz = hardware_efficient.hea(nlayers, nqubits, ["RY", "RZ"], "CNOT")
     qc = models.Circuit(nqubits)
-    qc.add(gates.X(_i) for _i in range(sum(mol.nelec)))
+    qc.add(gates.X(_i) for _i in range(mol.nelec))
     qc.add(hea_ansatz)
     qc.set_parameters(theta)
 

tests/test_hf_circuit.py

@@ -19,7 +19,7 @@ def test_jw_circuit():
     h2.run_pyscf()
 
     # JW-HF circuit
-    circuit = hf_circuit(h2.nso, sum(h2.nelec), ferm_qubit_map=None)
+    circuit = hf_circuit(h2.nso, h2.nelec, ferm_qubit_map=None)
 
     # Molecular Hamiltonian and the HF expectation value
     hamiltonian = h2.hamiltonian()
@@ -38,7 +38,7 @@ def test_bk_circuit_1():
     h2.run_pyscf()
 
     # JW-HF circuit
-    circuit = hf_circuit(h2.nso, sum(h2.nelec), ferm_qubit_map="bk")
+    circuit = hf_circuit(h2.nso, h2.nelec, ferm_qubit_map="bk")
 
     # Molecular Hamiltonian and the HF expectation value
     hamiltonian = h2.hamiltonian(ferm_qubit_map="bk")
@@ -55,7 +55,7 @@ def test_bk_circuit_2():
     lih.run_pyscf()
 
     # JW-HF circuit
-    circuit = hf_circuit(lih.nso, sum(lih.nelec), ferm_qubit_map="bk")
+    circuit = hf_circuit(lih.nso, lih.nelec, ferm_qubit_map="bk")
 
     # Molecular Hamiltonian and the HF expectation value
     hamiltonian = lih.hamiltonian(ferm_qubit_map="bk")

tests/test_molecule.py

@@ -213,18 +213,28 @@ def test_symbolic_hamiltonian():
     assert np.allclose(h2_sym_ham.matrix, ref_sym_ham.matrix)
 
 
-@pytest.mark.skip(reason="Psi4 doesn't offer pip install, so needs to be installed through conda or manually.")
-def test_run_psi4():
-    """PSI4 driver"""
-    # Hardcoded benchmark results
-    h2_ref_energy = -1.117349035
-    h2_ref_hcore = np.array([[-1.14765024, -1.00692423], [-1.00692423, -1.14765024]])
-
+def test_hamiltonian_input_error():
     h2 = Molecule([("H", (0.0, 0.0, 0.0)), ("H", (0.0, 0.0, 0.7))])
-    h2.run_psi4()
-
-    assert h2.e_hf == pytest.approx(h2_ref_energy)
-    assert np.allclose(h2.hcore, h2_ref_hcore)
+    h2.e_nuc = 0.0
+    h2.oei = np.random.rand(4, 4)
+    h2.tei = np.random.rand(4, 4, 4, 4)
+    with pytest.raises(NameError):
+        h2.hamiltonian("ihpc")
+
+
+# Commenting out since not actively supporting PSI4 at the moment
+# @pytest.mark.skip(reason="Psi4 doesn't offer pip install, so needs to be installed through conda or manually.")
+# def test_run_psi4():
+#     """PSI4 driver"""
+#     # Hardcoded benchmark results
+#     h2_ref_energy = -1.117349035
+#     h2_ref_hcore = np.array([[-1.14765024, -1.00692423], [-1.00692423, -1.14765024]])
+#
+#     h2 = Molecule([("H", (0.0, 0.0, 0.0)), ("H", (0.0, 0.0, 0.7))])
+#     h2.run_psi4()
+#
+#     assert h2.e_hf == pytest.approx(h2_ref_energy)
+#     assert np.allclose(h2.hcore, h2_ref_hcore)
 
 
 def test_expectation_value():
@@ -237,7 +247,7 @@ def test_expectation_value():
 
     # JW-HF circuit
     circuit = models.Circuit(h2.nso)
-    circuit.add(gates.X(_i) for _i in range(sum(h2.nelec)))
+    circuit.add(gates.X(_i) for _i in range(h2.nelec))
     # Molecular Hamiltonian and the HF expectation value
     hamiltonian = h2.hamiltonian()
     hf_energy = expectation(circuit, hamiltonian)