Wrapper function for loading a molecular Hamiltonian (#432)

* add load molecular hamiltonian function

* doc

src/openfermion/hamiltonians/_molecular_data.py

@@ -933,3 +933,52 @@ def get_molecular_rdm(self, use_fci=False):
         # Cast to InteractionRDM class.
         rdm = InteractionRDM(one_rdm, two_rdm)
         return rdm
+
+
+def load_molecular_hamiltonian(
+        geometry,
+        basis,
+        multiplicity,
+        description,
+        n_active_electrons=None,
+        n_active_orbitals=None):
+    """Attempt to load a molecular Hamiltonian with the given properties.
+
+    Args:
+        geometry: A list of tuples giving the coordinates of each atom.
+            An example is [('H', (0, 0, 0)), ('H', (0, 0, 0.7414))].
+            Distances in angstrom. Use atomic symbols to
+            specify atoms.
+        basis: A string giving the basis set. An example is 'cc-pvtz'.
+            Only optional if loading from file.
+        multiplicity: An integer giving the spin multiplicity.
+        description: A string giving a description.
+        n_active_electrons: An optional integer specifying the number of
+            electrons desired in the active space.
+        n_active_orbitals: An optional integer specifying the number of
+            spatial orbitals desired in the active space.
+
+    Returns:
+        The Hamiltonian as an InteractionOperator.
+    """
+
+    molecule = MolecularData(
+            geometry, basis, multiplicity, description=description)
+    molecule.load()
+
+    if n_active_electrons is None:
+        n_core_orbitals = 0
+        occupied_indices = None
+    else:
+        n_core_orbitals = (molecule.n_electrons - n_active_electrons) // 2
+        occupied_indices = list(range(n_core_orbitals))
+
+    if n_active_orbitals is None:
+        active_indices = None
+    else:
+        active_indices = list(range(n_core_orbitals,
+                                    n_core_orbitals + n_active_orbitals))
+
+    return molecule.get_molecular_hamiltonian(
+            occupied_indices=occupied_indices,
+            active_indices=active_indices)

src/openfermion/hamiltonians/_molecular_data_test.py

@@ -289,3 +289,19 @@ def test_abstract_molecule(self):
         correct_name = "Jellium_PlaneWave22_singlet"
         self.assertEqual(jellium_molecule.name, correct_name)
         os.remove("{}.hdf5".format(jellium_filename))
+
+    def test_load_molecular_hamiltonian(self):
+        bond_length = 1.45
+        geometry = [('Li', (0., 0., 0.)), ('H', (0., 0., bond_length))]
+
+        lih_hamiltonian = load_molecular_hamiltonian(
+                geometry, 'sto-3g', 1, format(bond_length), 2, 2)
+        self.assertEqual(count_qubits(lih_hamiltonian), 4)
+
+        lih_hamiltonian = load_molecular_hamiltonian(
+                geometry, 'sto-3g', 1, format(bond_length), 2, 3)
+        self.assertEqual(count_qubits(lih_hamiltonian), 6)
+
+        lih_hamiltonian = load_molecular_hamiltonian(
+                geometry, 'sto-3g', 1, format(bond_length), None, None)
+        self.assertEqual(count_qubits(lih_hamiltonian), 12)