orbital names are better tracked through transformations

src/tequila/quantumchemistry/chemistry_tools.py

@@ -804,20 +804,18 @@ class IntegralManager:
     _one_body_integrals: numpy.ndarray = None
     _two_body_integrals: NBodyTensor = None
     _constant_term: float = None
-    _basis_type: str = "unknown"
     _basis_name: str = "unknown"
     _orbital_type: str = "unknown" # e.g. "HF", "PNO", "native"
     _orbital_coefficients: numpy.ndarray = None
     _active_space: ActiveSpaceData = None
     _orbitals: typing.List[OrbitalData] = None
 
-    def __init__(self, one_body_integrals, two_body_integrals, basis_type="custom",
+    def __init__(self, one_body_integrals, two_body_integrals,
                  basis_name="unknown", orbital_type="unknown",
                  constant_term=0.0, orbital_coefficients=None, active_space=None, overlap_integrals=None, orbitals=None, *args, **kwargs):
         self._one_body_integrals = one_body_integrals
         self._two_body_integrals = two_body_integrals
         self._constant_term = constant_term
-        self._basis_type = basis_type
         self._basis_name = basis_name
         self._orbital_type = orbital_type
 
@@ -956,9 +954,16 @@ def transform_to_native_orbitals(self):
         """
         c = self.get_orthonormalized_orbital_coefficients()
         self.orbital_coefficients=c
-        self._orbital_type="orthonormalized-{}-basis".format(self._orbital_type)
+        self._orbital_type="orthonormalized-{}-basis".format(self._basis_name)
 
-    def transform_orbitals(self, U):
+    def is_unitary(self, U):
+        if len(U.shape) != 2: return False
+        if U.shape[0] != U.shape[1]: return False
+        test = (U.conj().T).dot(U) - numpy.eye(U.shape[0])
+        if not numpy.isclose(numpy.linalg.norm(test), 0.0): return False
+        return True
+
+    def transform_orbitals(self, U, name=None):
         """
         Transform orbitals
         Parameters
@@ -969,10 +974,12 @@ def transform_orbitals(self, U):
         -------
         updates the structure with new orbitals: c = cU
         """
-        c = self.orbital_coefficients
-        c = numpy.einsum("ix, xj -> ij", c, U, optimize="greedy")
-        self.orbital_coefficients = c
-        self._orbital_type += "-transformed"
+        assert self.is_unitary(U)
+        self.orbital_coefficients = numpy.einsum("ix, xj -> ij", self.orbital_coefficients, U, optimize="greedy")
+        if name is None:
+            self._orbital_type += "-transformed"
+        else:
+            self._orbital_type = name
 
     def get_integrals(self, orbital_coefficients=None, ordering="openfermion", ignore_active_space=False, *args, **kwargs):
         """
@@ -1001,7 +1008,9 @@ def get_integrals(self, orbital_coefficients=None, ordering="openfermion", ignor
             active_integrals = get_active_space_integrals(one_body_integrals=h, two_body_integrals=g,
                                                           occupied_indices=self._active_space.frozen_reference_orbitals,
                                                           active_indices=self._active_space.active_orbitals)
+
             c = active_integrals[0] + c
+
             h = active_integrals[1]
             g = NBodyTensor(elems=active_integrals[2], ordering="openfermion")
         g.reorder(to=ordering)
@@ -1070,7 +1079,6 @@ def __str__(self):
         return result
 
     def print_basis_info(self, *args, **kwargs) -> None:
-        print("{:15} : {}".format("basis_type", self._basis_type), *args, **kwargs)
         print("{:15} : {}".format("basis_name", self._basis_name), *args, **kwargs)
         print("{:15} : {}".format("orbital_type", self._orbital_type), *args, **kwargs)
         print("{:15} : {}".format("orthogonal", self.basis_is_orthogonal()), *args, **kwargs)

src/tequila/quantumchemistry/orbital_optimizer.py

@@ -37,7 +37,7 @@ def __call__(self, local_data, *args, **kwargs):
         self.iterations += 1
 
 def optimize_orbitals(molecule, circuit=None, vqe_solver=None, pyscf_arguments=None, silent=False,
-                      vqe_solver_arguments=None, initial_guess=None, return_mcscf=False, use_hcb=False, molecule_factory=None, molecule_arguments=None, *args, **kwargs):
+                      vqe_solver_arguments=None, initial_guess=None, return_mcscf=False, use_hcb=False, molecule_factory=None, molecule_arguments=None, restrict_to_active_space=True, *args, **kwargs):
     """
 
     Parameters
@@ -78,7 +78,12 @@ def optimize_orbitals(molecule, circuit=None, vqe_solver=None, pyscf_arguments=N
     if pyscf_arguments is None:
         pyscf_arguments = {"max_cycle_macro": 10, "max_cycle_micro": 3}
     no = molecule.n_orbitals
-    pyscf_molecule = QuantumChemistryPySCF.from_tequila(molecule=molecule, transformation=molecule.transformation)
+
+    if not isinstance(molecule, QuantumChemistryPySCF):
+        pyscf_molecule = QuantumChemistryPySCF.from_tequila(molecule=molecule, transformation=molecule.transformation)
+    else:
+        pyscf_molecule = molecule
+
     mf = pyscf_molecule._get_hf()
     result=OptimizeOrbitalsResult()
     mc = mcscf.CASSCF(mf, pyscf_molecule.n_orbitals, pyscf_molecule.n_electrons)
@@ -140,10 +145,11 @@ def optimize_orbitals(molecule, circuit=None, vqe_solver=None, pyscf_arguments=N
         mc.kernel()
     # make new molecule
 
-    transformed_molecule = pyscf_molecule.transform_orbitals(orbital_coefficients=mc.mo_coeff)
+    mo_coeff = mc.mo_coeff
+    transformed_molecule = pyscf_molecule.transform_orbitals(orbital_coefficients=mo_coeff, name="optimized")
     result.molecule=transformed_molecule
     result.old_molecule=molecule
-    result.mo_coeff=mc.mo_coeff
+    result.mo_coeff=mo_coeff
     result.energy=mc.e_tot
 
     if return_mcscf:

src/tequila/quantumchemistry/pyscf_interface.py

@@ -75,6 +75,7 @@ def __init__(self, parameters: ParametersQC,
             kwargs["two_body_integrals"] = g_ao
             kwargs["one_body_integrals"] = h_ao
             kwargs["orbital_coefficients"] = mo_coeff
+            kwargs["orbital_type"] = "hf"
 
             if "nuclear_repulsion" not in kwargs:
                 kwargs["nuclear_repulsion"] = mol.energy_nuc()

src/tequila/quantumchemistry/qc_base.py

@@ -94,7 +94,7 @@ def __init__(self, parameters: ParametersQC,
         else:
             self.integral_manager = self.initialize_integral_manager(active_orbitals=active_orbitals,
                                                                      reference_orbitals=reference_orbitals,
-                                                                     orbitals=orbitals, frozen_orbitals=frozen_orbitals, orbital_type=orbital_type, *args,
+                                                                     orbitals=orbitals, frozen_orbitals=frozen_orbitals, orbital_type=orbital_type, basis_name=self.parameters.basis_set, *args,
                                                                      **kwargs)
 
         if orbital_type is not None and orbital_type.lower() == "native":
@@ -109,15 +109,23 @@ def __init__(self, parameters: ParametersQC,
 
     @classmethod
     def from_tequila(cls, molecule, transformation=None, *args, **kwargs):
-        c, h1, h2 = molecule.get_integrals()
+        c = molecule.integral_manager.constant_term
+        h1 = molecule.integral_manager.one_body_integrals
+        h2 = molecule.integral_manager.two_body_integrals
+        S = molecule.integral_manager.overlap_integrals
+        active_orbitals = [o.idx_total for o in molecule.integral_manager.active_orbitals]
         if transformation is None:
             transformation = molecule.transformation
+        parameters = molecule.parameters
         return cls(nuclear_repulsion=c,
                    one_body_integrals=h1,
                    two_body_integrals=h2,
-                   n_electrons=molecule.n_electrons,
+                   overlap_integrals = S,
+                   orbital_coefficients = molecule.integral_manager.orbital_coefficients,
+                   active_orbitals= active_orbitals,
                    transformation=transformation,
-                   parameters=molecule.parameters, *args, **kwargs)
+                   orbital_type=molecule.integral_manager._orbital_type,
+                   parameters=parameters, *args, **kwargs)
 
     def supports_ucc(self):
         """
@@ -543,11 +551,13 @@ def initialize_integral_manager(self, *args, **kwargs):
 
         return manager
 
-    def transform_orbitals(self, orbital_coefficients, *args, **kwargs):
+    def transform_orbitals(self, orbital_coefficients, ignore_active_space=False, name=None, *args, **kwargs):
         """
         Parameters
         ----------
-        orbital_coefficients: second index is new orbital indes, first is old orbital index (summed over)
+        orbital_coefficients: second index is new orbital indes, first is old orbital index (summed over), indices are assumed to be defined on the active space
+        ignore_active_space: if true orbital_coefficients are not assumed to be given in the active space
+        name: str, name the new orbitals
         args
         kwargs
 
@@ -556,9 +566,20 @@ def transform_orbitals(self, orbital_coefficients, *args, **kwargs):
         New molecule with transformed orbitals
         """
 
+        U = numpy.eye(self.integral_manager.orbital_coefficients.shape[0])
+        # mo_coeff by default only acts on the active space
+        active_indices = [o.idx_total for o in self.integral_manager.active_orbitals]
+
+        if ignore_active_space:
+            U = orbital_coefficients
+        else:
+            for kk,k in enumerate(active_indices):
+                for ll,l in enumerate(active_indices):
+                    U[k][l] = orbital_coefficients[kk][ll]
+
         # can not be an instance of a specific backend (otherwise we get inconsistencies with classical methods in the backend)
         integral_manager = copy.deepcopy(self.integral_manager)
-        integral_manager.transform_orbitals(U=orbital_coefficients)
+        integral_manager.transform_orbitals(U=U, name=name)
         result = QuantumChemistryBase(parameters=self.parameters, integral_manager=integral_manager, transformation=self.transformation)
         return result
 
@@ -583,7 +604,7 @@ def use_native_orbitals(self, inplace=False):
         else:
             integral_manager = copy.deepcopy(self.integral_manager)
             integral_manager.transform_to_native_orbitals()
-            result = QuantumChemistryBase(parameters=self.parameters, integral_manager=integral_manager, orbital_type="native", transformation=self.transformation)
+            result = QuantumChemistryBase(parameters=self.parameters, integral_manager=integral_manager, transformation=self.transformation)
             return result
 
 
@@ -867,13 +888,13 @@ def hcb_to_me(self, U=None, condensed=False):
         """
         if U is None:
             U = QCircuit()
-
-        # consistency
-        consistency = [x < self.n_orbitals for x in U.qubits]
-        if not all(consistency):
-            warnings.warn(
-                "hcb_to_me: given circuit is not defined on the first {} qubits. Is this a HCB circuit?".format(
-                    self.n_orbitals))
+        else:
+            ups = [self.transformation.up(i.idx) for i in self.orbitals]
+            consistency = [x in ups for x in U.qubits]
+            if not all(consistency):
+                warnings.warn(
+                    "hcb_to_me: given circuit is not defined on all first {} qubits. Is this a HCB circuit?".format(
+                        self.n_orbitals))
 
         # map to alpha qubits
         if condensed: