Merge pull request #340 from tequilahub/kottmanj-update-1.9.5

* fixing syntax issue in post_init of dataclass (#327)
* phoenics support dropped
* more convenient randomization initialization for OO, avoiding numpy warnings
* Update qasm.py (#334)
Fixes issue #332
* added methods to create annihilation, creation, sz, sp, sm and s2 operators in qubit representation (#336)
* Update qasm.py (#335)
* Update ci_chemistry_psi4.yml (psi4 --> conda-forge channel)
* orbital names are better tracked through transformations (#338)
* keep basis info when re-initializing with mol.from_tequila
* orbital-optimization gives back molecule in the original basis
* more basis info (#339)
---------

Co-authored-by: davibinco <[email protected]>

Author: kottmanj

Diff for: .github/workflows/ci_chemistry_psi4.yml

@@ -16,7 +16,7 @@ jobs:
     runs-on: ubuntu-latest
     strategy:
       matrix:
-        python-version: [3.8]
+        python-version: [3.9]
 
     steps:
     - uses: actions/checkout@v2
@@ -50,7 +50,7 @@ jobs:
         source $HOME/.bashrc
         source $CONDABASE/bin/activate
         conda activate test_psi4
-        conda install psi4 python=3.8 -c psi4 
+        conda install psi4 python=3.9 -c conda-forge
         python -m pip install --upgrade pip
         python -m pip install -r requirements.txt
         python -m pip install -e .

Diff for: README.md

@@ -14,7 +14,7 @@ Tequila can execute the underlying quantum expectation values on state of the ar
 - [talks and slides](https://kottmanj.github.io/talks_and_material/)  
 
 # Installation
-Recommended Python version is 3.8-3.9.   
+Recommended Python version is 3.9 - 3.10.   
 Tequila supports linux, osx and windows. However, not all optional dependencies are supported on windows.  
 
 ## Install from PyPi
@@ -275,7 +275,7 @@ Currently supported
 ### [Psi4](https://github.com/psi4/psi4).
 In a conda environment this can be installed with
 ```bash
-conda install psi4 -c psi4
+conda install psi4 -c conda-forge
 ```
 Here is a small [tutorial](https://nbviewer.org/github/tequilahub/tequila-tutorials/blob/main/chemistry/ChemistryModule.ipynb) that illustrates the usage.
 

Diff for: src/tequila/circuit/qasm.py

@@ -313,19 +313,33 @@ def parse_command(command: str, custom_gates_map: Dict[str, QCircuit], qregister
         return apply_custom_gate(custom_circuit=custom_circuit, qregisters_values=qregisters_values)
 
     if name in ("x", "y", "z", "h", "cx", "cy", "cz", "ch"):
-        return QCircuit.wrap_gate(gates.impl.QGateImpl(name=(name[1] if name[0] == 'c' else name).upper(),
-                                            control=get_qregister(args[0], qregisters) if name[0] == 'c' else None,
-                                            target=get_qregister(args[1 if name[0] == 'c' else 0], qregisters)))
+        target = get_qregister(args[0], qregisters)
+        control = None
+        if name[0].lower() == 'c':
+            control = get_qregister(args[0], qregisters)
+            target = get_qregister(args[1], qregisters)
+            name = name[1]
+        G = getattr(gates, name.upper())
+        return G(control=control, target=target)
+
     if name in ("ccx", "ccy", "ccz"):
-        return QCircuit.wrap_gate(gates.impl.QGateImpl(name=name.upper()[2],
-                                            control=[get_qregister(args[0], qregisters), get_qregister(args[1], qregisters)],
-                                            target=get_qregister(args[2], qregisters)))
+        G = getattr(gates, name[2].upper())
+        control = [get_qregister(args[0], qregisters), get_qregister(args[1], qregisters)]
+        target = get_qregister(args[2], qregisters)
+        return G(control=control, target=target)
+
     if name.startswith("rx(") or name.startswith("ry(") or name.startswith("rz(") or \
         name.startswith("crx(") or name.startswith("cry(") or name.startswith("crz("):
-        return QCircuit.wrap_gate(gates.impl.RotationGateImpl(axis=name[2 if name[0] == 'c' else 1],
-                                                   angle=get_angle(name)[0],
-                                                   control=get_qregister(args[0], qregisters) if name[0] == 'c' else None,
-                                                   target=get_qregister(args[1 if name[0] == 'c' else 0], qregisters)))
+        angle = get_angle(name)[0]
+        i = name.find('(')
+        name = name[0:i]
+        name = name.upper()
+        name = [x for x in name]
+        name[-1] = name[-1].lower()
+        name = "".join(name)
+        G = getattr(gates, name)
+        return G(angle=angle,control=get_qregister(args[0], qregisters) if name[0] == 'C' else None,target=get_qregister(args[1 if name[0] == 'C' else 0], qregisters))
+            
     if name.startswith("U("):
         angles = get_angle(name)
         return gates.U(theta=angles[0], phi=angles[1], lambd=angles[2],
@@ -362,7 +376,7 @@ def parse_command(command: str, custom_gates_map: Dict[str, QCircuit], qregister
                   control=get_qregister(args[0], qregisters),
                   target=get_qregister(args[1], qregisters))
     if name in ("s", "t", "sdg", "tdg"):
-        g = gates.Phase(pi / (2 if name.startswith("s") else 4),
+        g = gates.Phase(angle=pi / (2 if name.startswith("s") else 4),
                      control=None,
                      target=get_qregister(args[0], qregisters))
         if name.find("dg") != -1:

Diff for: src/tequila/quantumchemistry/__init__.py

@@ -95,6 +95,8 @@ def Molecule(geometry: str = None,
     if backend is None:
         if basis_set is None or basis_set.lower() in ["madness", "mra", "pno"]:
             backend = "madness"
+            basis_set = "mra"
+            parameters.basis_set = basis_set
             if orbital_type is not None and orbital_type.lower() not in ["pno", "mra-pno"]:
                 warnings.warn("only PNOs supported as orbital_type without basis set. Setting to pno - You gave={}".format(orbital_type), TequilaWarning)
             orbital_type = "pno"

Diff for: 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)
@@ -1069,14 +1078,16 @@ def __str__(self):
             result += str(x) + "\n"
         return result
 
-    def print_basis_info(self, *args, **kwargs) -> None:
-        print("{:15} : {}".format("basis_type", self._basis_type), *args, **kwargs)
+    def print_basis_info(self, print_coefficients=True, *args, **kwargs) -> None:
         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)
-        print("{:15} : {}".format("functions", self.one_body_integrals.shape[0]), *args, **kwargs)
+        print("{:15} : {}".format("orthogonal basis", self.basis_is_orthogonal()), *args, **kwargs)
+        print("{:15} : {}".format("basis functions", self.one_body_integrals.shape[0]), *args, **kwargs)
+        print("{:15} : {}".format("active orbitals", [o.idx_total for o in self.active_orbitals]), *args, **kwargs)
         print("{:15} : {}".format("reference", [x.idx_total for x in self.reference_orbitals]), *args, **kwargs)
 
+        if not print_coefficients: return
+
         print("Current Orbitals", *args, **kwargs)
         for i,x in enumerate(self.orbitals):
             print(x, *args, **kwargs)

Diff for: 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:

Diff for: 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()