Merge pull request #302 from tequilahub/devel

Update to 1.8.9

requirements.txt

@@ -1,6 +1,6 @@
 # necessary
 numpy
-scipy >= 1.5 # can in principle be smaller, we recommend >= 1.5 for consistency with our tutorials (numerical gradients mostly)
+scipy < 1.11 # for now, until issues with scipy/pyscf are fixed # can in principle be smaller, we recommend >= 1.5 for consistency with our tutorials (numerical gradients mostly)
 sympy
 #jax
 #jaxlib
@@ -9,7 +9,7 @@ setuptools
 pytest
 openfermion ~= 1.0 # can not be smaller than 1.0 
 #cmake # needed by qulacs, can be removed otherwise, now in qulacs requirements
-qulacs < 0.5.0 # default simulator (best integration), remove if the installation gives you trouble and just install one of the other supported backend. Version restriction only for noise models, otherwise the new version is fine
+qulacs # default simulator (best integration), remove if the installation gives you trouble and just install one of the other supported backend. Version restriction only for noise models, otherwise the new version is fine
 
 #optional quantum backends
 #cirq >= 0.9.2 #

src/tequila/apps/robustness/interval.py

@@ -190,7 +190,7 @@ def _compute_interval(self) -> Tuple[float, float]:
                 pauli_lower_bound = pauli_upper_bound = 1.0
             else:
                 expec_normalized = np.clip(2 * (p_expec - min_eigval) / (max_eigval - min_eigval) - 1, -1, 1,
-                                           dtype=np.float64)
+                                           dtype=float)
 
                 pauli_lower_bound = (max_eigval - min_eigval) / 2.0 * (
                         1 + self._calc_lower_bound(expec_normalized, self.fidelity)) + min_eigval
@@ -337,7 +337,7 @@ def _compute_bound_grouped(self) -> float:
         for eigvals, expec, variance in zip(self._pauligroups_eigenvalues, self._pauligroups_expectations,
                                             self._pauligroups_variances):
             min_eigval = min(eigvals)
-            expec_pos = np.clip(expec - min_eigval, 0, None, dtype=np.float)
+            expec_pos = np.clip(expec - min_eigval, 0, None, dtype=float)
             bound += min_eigval + self._calc_lower_bound(expec_pos, variance, self.fidelity)
 
         return bound

src/tequila/circuit/_gates_impl.py

@@ -180,18 +180,18 @@ def parameter(self, other):
 
     def __init__(self, name, parameter: UnionParam, target: UnionList, control: UnionList = None,
                 generator: QubitHamiltonian = None):
-        super().__init__(name=name, target=target, control=control, generator=generator)
         # failsafe
         if hasattr(parameter, "shape") and parameter.shape not in [tuple()]: # take care of new numpy conventions where scalars have shape ()
+            self._parameter=None
             raise TequilaException("parameter has to be a scalar. Received {}\n{}\n{}".format(repr(parameter), type(parameter), str(parameter)))
         self._parameter = assign_variable(variable=parameter)
+        super().__init__(name=name, target=target, control=control, generator=generator)
 
     def __str__(self):
         result = str(self.name) + "(target=" + str(self.target)
         if not self.is_single_qubit_gate():
             result += ", control=" + str(self.control)
-
-        result += ", parameter=" + repr(self._parameter)
+        result += ", parameter=" + repr(self.parameter)
         result += ")"
         return result
 

src/tequila/circuit/qasm.py

@@ -402,16 +402,21 @@ def get_angle(name: str) -> list:
                 raise TequilaException("Invalid specification {}".format(name))
             angle = angle.replace('pi', '')
             try:
+                sign = 1
+                div = 1
+                if angle.find('-') != -1:
+                    angle = angle.replace('-', '')
+                    sign = -1
+                if angle.find('/') != -1:
+                    div = float(angle[angle.index('/')+1:])
+                    angle = angle[:angle.index('/')]
                 if angle.find('*') != -1:
                     angle = angle.replace('*', '')
-                    phase = float(angle) * pi
-                elif angle.find('/') != -1:
-                    angle = angle.replace('/', '')
-                    phase = pi / float(angle)
+                    phase = sign * float(angle) * pi / div
                 elif len(angle) == 0:
-                    phase = pi
+                    phase = sign * pi / div
                 else:
-                    phase = float(angle)
+                    phase = sign * float(angle) / div
             except ValueError:
                 raise TequilaException("Invalid specification {}".format(name))
         angles.append(phase)

src/tequila/hamiltonian/qubit_hamiltonian.py

@@ -523,7 +523,7 @@ def split(self, *args, **kwargs) -> tuple:
         hermitian = QubitHamiltonian.zero()
         anti_hermitian = QubitHamiltonian.zero()
         for k, v in self.qubit_operator.terms.items():
-            hermitian.qubit_operator.terms[k] = np.float(v.real)
+            hermitian.qubit_operator.terms[k] = v.real
             anti_hermitian.qubit_operator.terms[k] = 1.j * v.imag
 
         return hermitian.simplify(), anti_hermitian.simplify()

src/tequila/ml/interface_torch.py

@@ -128,7 +128,7 @@ def backward(ctx, grad_backward):
                     g_keys = [j for j in grads.keys()]
                     probe = grads[g_keys[0]]  # first entry will tell us number of output
                     dims = len(g_keys), len(probe)
-                    arr = np.empty(dims, dtype=np.float)
+                    arr = np.empty(dims, dtype=float)
                     for j, key in enumerate(g_keys):
                         line = grads[key]
                         for k, ob in enumerate(line):

src/tequila/optimizers/optimizer_base.py

@@ -36,10 +36,18 @@ def iterations(self):
     angles: typing.List[typing.Dict[str, numbers.Number]] = field(default_factory=list)
 
     # history of all function evaluations
-    energies_calls: typing.List[numbers.Real] = field(default_factory=list)
-    gradients_calls: typing.List[typing.Dict[str, numbers.Real]] = field(default_factory=list)
+    energy_calls: typing.List[numbers.Real] = field(default_factory=list)
+    gradient_calls: typing.List[typing.Dict[str, numbers.Real]] = field(default_factory=list)
     angles_calls: typing.List[typing.Dict[str, numbers.Number]] = field(default_factory=list)
-
+
+    # backward comp.
+    @property
+    def energies_calls(self):
+        return self.energy_calls
+    @property
+    def energies_evaluations(self):
+        return self.energy_calls
+
     def __add__(self, other):
         """
         magic method for convenient combination of history objects.

src/tequila/optimizers/optimizer_scipy.py

@@ -275,15 +275,15 @@ def __call__(self, *args, **kwargs):
 
         if self.save_history:
             self.history.energies = callback.energies
-            self.history.energy_evaluations = E.history
+            self.history.energy_calls = E.history
             self.history.angles = callback.angles
-            self.history.angles_evaluations = E.history_angles
+            self.history.angles_calls = E.history_angles
             self.history.gradients = callback.gradients
             self.history.hessians = callback.hessians
             if dE is not None and not isinstance(dE, str):
-                self.history.gradients_evaluations = dE.history
+                self.history.gradient_calls = dE.history
             if ddE is not None and not isinstance(ddE, str):
-                self.history.hessians_evaluations = ddE.history
+                self.history.hessian_calls = ddE.history
 
             # some methods like "cobyla" do not support callback functions
             if len(self.history.energies) == 0:

src/tequila/quantumchemistry/encodings.py

@@ -26,6 +26,15 @@ def known_encodings():
 
 class EncodingBase:
 
+    # true if the encoding is fully integrated
+    # false: can only do special things (like building the Hamiltionian)
+    # but is not consistent with UCC gate generation
+    _ucc_support=False
+
+    @property
+    def supports_ucc(self):
+        return self._ucc_support
+
     @property
     def name(self):
         prefix=""
@@ -129,6 +138,7 @@ class JordanWigner(EncodingBase):
     """
     OpenFermion::jordan_wigner
     """
+    _ucc_support=True
 
     def do_transform(self, fermion_operator:openfermion.FermionOperator, *args, **kwargs) -> openfermion.QubitOperator:
         return openfermion.jordan_wigner(fermion_operator, *args, **kwargs)
@@ -152,6 +162,8 @@ class BravyiKitaev(EncodingBase):
     Uses OpenFermion::bravyi_kitaev
     """
 
+    _ucc_support=True
+
     def do_transform(self, fermion_operator:openfermion.FermionOperator, *args, **kwargs) -> openfermion.QubitOperator:
         return openfermion.bravyi_kitaev(fermion_operator, n_qubits=self.n_orbitals*2)
 
@@ -161,6 +173,8 @@ class BravyiKitaevTree(EncodingBase):
     Uses OpenFermion::bravyi_kitaev_tree
     """
 
+    _ucc_support=True
+
     def do_transform(self, fermion_operator:openfermion.FermionOperator, *args, **kwargs) -> openfermion.QubitOperator:
         return openfermion.bravyi_kitaev_tree(fermion_operator, n_qubits=self.n_orbitals*2)
 
@@ -169,6 +183,8 @@ class BravyiKitaevFast(EncodingBase):
     Uses OpenFermion::bravyi_kitaev_tree
     """
 
+    _ucc_support=False
+
     def do_transform(self, fermion_operator:openfermion.FermionOperator, *args, **kwargs) -> openfermion.QubitOperator:
         n_qubits = openfermion.count_qubits(fermion_operator)
         if n_qubits != self.n_orbitals*2:
@@ -177,6 +193,9 @@ def do_transform(self, fermion_operator:openfermion.FermionOperator, *args, **kw
         return openfermion.bravyi_kitaev_fast(op)
 
 class TaperedBravyKitaev(EncodingBase):
+
+    _ucc_support=False
+
     """
     Uses OpenFermion::symmetry_conserving_bravyi_kitaev (tapered bravyi_kitaev_tree arxiv:1701.07072)
     Reduces Hamiltonian by 2 qubits

src/tequila/quantumchemistry/madness_interface.py

@@ -37,6 +37,47 @@ def find_executable(madness_root_dir=None):
             executable = shutil.which("{}/src/apps/pno/pno_integrals".format(madness_root_dir))
         return executable
 
+    def plot2cube(self, orbital, filename=None, *args, **kwargs):
+        """
+        plot orbitals to cube file (needs madtequila backend installed)
+        Parameters
+        ----------
+        method: orbital, the orbital index (starting from 0 on the active orbitals)
+                if you want to plot frozen orbitals you can hand in a Tequila Orbital structure with idx_total defined
+        filename: name of the cubefile (default: mra_orbital_X.cube where X is the total index of the active orbital)
+        args: further arguments for plot2cube
+        kwargs further keyword arguments for plot2cube
+
+        see here for more https://github.com/kottmanj/madness/tree/tequila/src/apps/plot
+        """
+
+        plot2cube = shutil.which("plot2cube")
+        if plot2cube is None:
+            raise TequilaMadnessException(
+                "can't plot to cube file. Couldn't find plot2cube executable.\n\nTry installing\n\t conda install madtequila -c kottmann\nand assure the version is >2.3")
+
+        if hasattr(orbital,"idx"):
+            idx = orbital.idx
+        else:
+            idx = self.orbitals[orbital].idx_total
+
+        callist = [plot2cube, "file=mra_orbital_{}".format(idx)]
+
+        if filename is not None:
+            callist.append("outfile={}".format(filename))
+        for k, v in kwargs.items():
+            callist.append("{}={}".format(k, v))
+        for k in args:
+            callist.append("{}".format(k))
+
+        import subprocess
+        try:
+            with open("plot2cube_{}.log".format(orbital), "w") as logfile:
+                subprocess.call(callist, stdout=logfile)
+        except:
+            print("plotting failed ....")
+            print("see plot2cube_{}.log".format(orbital))
+
     def __init__(self, parameters: ParametersQC,
                  transformation: typing.Union[str, typing.Callable] = None,
                  active_orbitals: list = "auto",
@@ -93,10 +134,10 @@ def __init__(self, parameters: ParametersQC,
                     warnings.warn("MADNESS did not terminate as expected! status = {}".format(status), TequilaWarning)
                 status += str(madness_status) + "\n"
             except Exception as E:
-                status += "madness_run={}\n".format(str(E))
+                status += str(E) + "\n"
 
             # will read the binary files, convert them and save them with the right name
-            h, g, pinfo= self.convert_madness_output_from_bin_to_npy(name=name, datadir=datadir)
+            h, g, pinfo = self.convert_madness_output_from_bin_to_npy(name=name, datadir=datadir)
             status += "found {}_htensor.npy={}\n".format(name, "failed" not in h)
             status += "found {}_gtensor.npy={}\n".format(name, "failed" not in g)
             status += "found {}_pnoinfo.txt={}\n".format(name, "failed" not in pinfo)
@@ -106,15 +147,17 @@ def __init__(self, parameters: ParametersQC,
             status += str(g)
             status += "pnoinfo report:\n"
             status += str(pinfo)
+
+            solution = "Solution 1: Assuming precomputed files are available:\n    provide {name}_gtensor.npy, {name}_htensor.npy and {name}_pnoinfo.txt\n    and call the Molecule constructor with n_pno='read' keyword \n\nSolution 2: Try installing with conda\n    conda install madtequila -c kottmann\n\nSolution 3: Install from source\n    follow instructions on github.com/kottmanj/madness".format(
+                name=name)
+            if self.executable is not None:
+                solution = "madness executable was found, but calculation did not succeed, check {name}_pno_integrals.out for clues".format(
+                    name=name)
+
             if "failed" in h or "failed" in g:
                 raise TequilaMadnessException("Could not initialize the madness interface\n"
                                               "Status report is\n"
-                                              "{status}\n"
-                                              "either provide {name}_gtensor.npy and {name}_htensor.npy files\n"
-                                              "or provide the number of pnos over by giving the n_pnos keyword to run madness\n"
-                                              "in order for madness to run you need to make sure that the pno_integrals executable can be found in your environment\n"
-                                              "alternatively you can provide the path to the madness_root_dir: the directory where you compiled madness\n".format(
-                    name=name, status=status))
+                                              "{status}\n\n".format(status=status) + solution)
         # get additional information from madness file
         nuclear_repulsion = 0.0
         pairinfo = None
@@ -142,12 +185,12 @@ def __init__(self, parameters: ParametersQC,
 
         if pairinfo is None:
             raise TequilaMadnessException("Pairinfo from madness calculation not found\nPlease provide pnoinfo.txt")
-        
+
         n_orbitals_total = h.shape[0]
         if "n_orbitals" in kwargs:
             # this would be the active orbitals
             kwargs.pop("n_orbitals")
-            
+
         assert h.shape[1] == n_orbitals_total
         assert sum(g.shape) == 4 * n_orbitals_total
         assert len(g.shape) == 4
@@ -219,8 +262,7 @@ def cleanup(self, warn=False, delete_all_files=False):
 
     def run_madness(self, *args, **kwargs):
         if self.executable is None:
-            return "pno_integrals executable not found\n" \
-                   "pass over executable keyword or export MAD_ROOT_DIR to system environment"
+            return "\n\n----> pno_integrals executable not found <----\n\n"
         self.write_madness_input(n_pno=self.n_pno, n_virt=self.n_virt, *args, **kwargs)
 
         # prevent reading in old files
@@ -492,15 +534,19 @@ def make_hardcore_boson_pno_upccd_ansatz(self, pairs=None, label=None, include_r
                         c = [a.idx, a.idx]
                     else:
                         c = [x.idx, x.idx]
+
+            alpha_map = {k.idx: self.transformation.up(k.idx) for k in self.orbitals}
+            U = U.map_qubits(alpha_map)
         else:
             for x in pairs:
                 if include_reference:
-                    U += gates.X(x.idx)
+                    U += gates.X(self.transformation.up(x.idx))
                 for a in self.get_pair_orbitals(i=x, j=x):
                     if x == a:
                         continue
                     idx = self.format_excitation_indices([(x.idx, a.idx)])
                     U += self.make_hardcore_boson_excitation_gate(indices=idx, angle=(idx, "D", label))
+
         return U
 
     def make_upccgsd_indices(self, label=None, name="UpCCGD", exclude: list = None, *args, **kwargs):
@@ -612,18 +658,20 @@ def write_madness_input(self, n_pno=None, n_virt=0, filename="input", maxrank=No
         n_pairs = n_electrons // 2
         if n_orbitals is None:
             n_orbitals = n_electrons  # minimal correlated (each active pair will have one virtual)
-        
+
         if n_pno is None:
             n_pno = n_orbitals - n_pairs
 
         if maxrank is None:
             # need at least maxrank=1, otherwise no PNOs are computed
             # this was a bug in <=v1.8.5 
-            maxrank = max(1,int(numpy.ceil(n_pno // n_pairs)))
-
-        if maxrank<=0:
-            warnings.warn("maxrank={} in tequila madness backend! No PNOs will be computed. Set the value when initializing the Molecule as tq.Molecule(..., pno={\"maxrank\":1, ...})".format(maxrank), TequilaWarning)
-
+            maxrank = max(1, int(numpy.ceil(n_pno // n_pairs)))
+
+        if maxrank <= 0:
+            warnings.warn(
+                "maxrank={} in tequila madness backend! No PNOs will be computed. Set the value when initializing the Molecule as tq.Molecule(..., pno={\"maxrank\":1, ...})".format(
+                    maxrank), TequilaWarning)
+
         data = {}
         if self.parameters.multiplicity != 1:
             raise TequilaMadnessException(