diff --git a/requirements.txt b/requirements.txt index 5663623a..a67423b5 100644 --- a/requirements.txt +++ b/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 # diff --git a/src/tequila/apps/robustness/interval.py b/src/tequila/apps/robustness/interval.py index 5331efc0..9350f2ab 100644 --- a/src/tequila/apps/robustness/interval.py +++ b/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 diff --git a/src/tequila/circuit/_gates_impl.py b/src/tequila/circuit/_gates_impl.py index 7df34874..68f9a870 100644 --- a/src/tequila/circuit/_gates_impl.py +++ b/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 diff --git a/src/tequila/circuit/qasm.py b/src/tequila/circuit/qasm.py index ef519dc0..65878f76 100644 --- a/src/tequila/circuit/qasm.py +++ b/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) diff --git a/src/tequila/hamiltonian/qubit_hamiltonian.py b/src/tequila/hamiltonian/qubit_hamiltonian.py index ccc77b7f..af90991c 100644 --- a/src/tequila/hamiltonian/qubit_hamiltonian.py +++ b/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() diff --git a/src/tequila/ml/interface_torch.py b/src/tequila/ml/interface_torch.py index b191a574..f96d134b 100644 --- a/src/tequila/ml/interface_torch.py +++ b/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): diff --git a/src/tequila/optimizers/optimizer_base.py b/src/tequila/optimizers/optimizer_base.py index 61af4cf2..23cc1b06 100644 --- a/src/tequila/optimizers/optimizer_base.py +++ b/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. diff --git a/src/tequila/optimizers/optimizer_scipy.py b/src/tequila/optimizers/optimizer_scipy.py index 943db94c..4d4d65a2 100644 --- a/src/tequila/optimizers/optimizer_scipy.py +++ b/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: diff --git a/src/tequila/quantumchemistry/encodings.py b/src/tequila/quantumchemistry/encodings.py index cd42b014..12d5ed46 100644 --- a/src/tequila/quantumchemistry/encodings.py +++ b/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 diff --git a/src/tequila/quantumchemistry/madness_interface.py b/src/tequila/quantumchemistry/madness_interface.py index c2c8fb65..e56d58e3 100644 --- a/src/tequila/quantumchemistry/madness_interface.py +++ b/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( diff --git a/src/tequila/quantumchemistry/orbital_optimizer.py b/src/tequila/quantumchemistry/orbital_optimizer.py index 0da97b7e..34b593ae 100644 --- a/src/tequila/quantumchemistry/orbital_optimizer.py +++ b/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, *args, **kwargs): + vqe_solver_arguments=None, initial_guess=None, return_mcscf=False, use_hcb=False, *args, **kwargs): """ Parameters @@ -49,6 +49,7 @@ def optimize_orbitals(molecule, circuit=None, vqe_solver=None, pyscf_arguments=N A customized object can be passed that needs to be callable with the following signature: vqe_solver(H=H, circuit=self.circuit, molecule=molecule, **self.vqe_solver_arguments) pyscf_arguments: Arguments for the MCSCF structure of PySCF, if None, the defaults are {"max_cycle_macro":10, "max_cycle_micro":3} (see here https://pyscf.org/pyscf_api_docs/pyscf.mcscf.html) silent: silence printout + use_hcb: indicate if the circuit is in hardcore Boson encoding vqe_solver_arguments: Optional arguments for a customized vqe_solver or the default solver for the default solver: vqe_solver_arguments={"optimizer_arguments":A, "restrict_to_hcb":False} where A holds the kwargs for tq.minimize restrict_to_hcb keyword controls if the standard (in whatever encoding the molecule structure has) Hamiltonian is constructed or the hardcore_boson hamiltonian @@ -86,6 +87,16 @@ def optimize_orbitals(molecule, circuit=None, vqe_solver=None, pyscf_arguments=N if circuit is None and vqe_solver is None: raise Exception("optimize_orbitals: Either provide a circuit or a callable vqe_solver") + if use_hcb: + if vqe_solver_arguments is None: + vqe_solver_arguments={} + vqe_solver_arguments["restrict_to_hcb"]=True + # consistency check + n_qubits = len(circuit.qubits) + n_orbitals = molecule.n_orbitals + if n_qubits > n_orbitals: + warnings.warn("Potential inconsistency in orbital optimization: use_hcb is switched on but we have\n n_qubits={} in the circuit\n n_orbital={} in the molecule\n".format(n_qubits,n_orbitals), TequilaWarning) + wrapper = PySCFVQEWrapper(molecule_arguments=pyscf_molecule.parameters, n_electrons=pyscf_molecule.n_electrons, const_part=c, circuit=circuit, vqe_solver_arguments=vqe_solver_arguments, silent=silent, vqe_solver=vqe_solver, *args, **kwargs) @@ -203,18 +214,8 @@ def kernel(self, h1, h2, *args, **kwargs): else: # static ansatz U = self.circuit - - if restrict_to_hcb: - # todo: adapt compute_rdms function to operate in HCB encoding -> faster and less measurements - U = molecule.hcb_to_me(U=U) - warnings.warn("optimize_orbitals: restrict_to_hcb=True, will use HCB for VQE but will map back to JW for the RDMs -> not fully optimized, see lines in code below this warning for potential speedup :-)", TequilaWarning) - # should look like this: - #rdm1 = .... compute rdm1 from hcb wavefunction - #rdm2 = .... compute rdm2 from hcb wavefunction - # then wrap the line below to an else statement - rdm1, rdm2 = molecule.compute_rdms(U=U, variables=result.variables, spin_free=True, get_rdm1=True, - get_rdm2=True) + rdm1, rdm2 = molecule.compute_rdms(U=U, variables=result.variables, spin_free=True, get_rdm1=True, get_rdm2=True, use_hcb=restrict_to_hcb) rdm2 = self.reorder(rdm2, 'dirac', 'mulliken') if not self.silent: print("{:20} : {}".format("energy", result.energy)) diff --git a/src/tequila/quantumchemistry/psi4_interface.py b/src/tequila/quantumchemistry/psi4_interface.py index d75d7760..7fc5ae7f 100644 --- a/src/tequila/quantumchemistry/psi4_interface.py +++ b/src/tequila/quantumchemistry/psi4_interface.py @@ -576,7 +576,7 @@ def rdm2(self) -> tuple: def compute_rdms(self, U: QCircuit = None, variables: Variables = None, spin_free: bool = True, get_rdm1: bool = True, get_rdm2: bool = True, psi4_method: str = None, - psi4_options: dict = {}): + psi4_options: dict = {}, *args, **kwargs): """ Same functionality as qc_base.compute_rdms (look there for more information), plus the additional option to compute 1- and 2-RDM using psi4 by the keyword psi4_rdms @@ -605,8 +605,8 @@ def compute_rdms(self, U: QCircuit = None, variables: Variables = None, spin_fre ------- """ if not psi4_method: - super().compute_rdms(U=U, variables=variables, spin_free=spin_free, - get_rdm1=get_rdm1, get_rdm2=get_rdm2) + return super().compute_rdms(U=U, variables=variables, spin_free=spin_free, + get_rdm1=get_rdm1, get_rdm2=get_rdm2, *args, **kwargs) else: # Get 1- and 2-particle reduced density matrix via Psi4 CISD computation # If "cisd" is chosen, change to "detci" (default is excitation level 2 anyhow) to obtain a CIWavefunction diff --git a/src/tequila/quantumchemistry/pyscf_interface.py b/src/tequila/quantumchemistry/pyscf_interface.py index e7f1b25f..d7217540 100644 --- a/src/tequila/quantumchemistry/pyscf_interface.py +++ b/src/tequila/quantumchemistry/pyscf_interface.py @@ -43,7 +43,7 @@ def __init__(self, parameters: ParametersQC, else: mol.symmetry = True - mol.build() + mol.build(parse_arg=False) # solve restricted HF mf = pyscf.scf.RHF(mol) @@ -120,7 +120,7 @@ def _get_hf(self, do_not_solve=True, **kwargs): mo_occ = numpy.zeros(norb) mo_occ[:nelec // 2] = 2 - pyscf_mol = pyscf.gto.M(verbose=0) + pyscf_mol = pyscf.gto.M(verbose=0, parse_arg=False) pyscf_mol.nelectron = nelec pyscf_mol.incore_anyway = True # ensure that custom integrals are used pyscf_mol.energy_nuc = lambda *args: c diff --git a/src/tequila/quantumchemistry/qc_base.py b/src/tequila/quantumchemistry/qc_base.py index e26cf653..385ced81 100644 --- a/src/tequila/quantumchemistry/qc_base.py +++ b/src/tequila/quantumchemistry/qc_base.py @@ -3,7 +3,7 @@ from tequila import TequilaException, BitString, TequilaWarning from tequila.hamiltonian import QubitHamiltonian -from tequila.hamiltonian.paulis import Sp, Sm +from tequila.hamiltonian.paulis import Sp, Sm, Zero from tequila.circuit import QCircuit, gates from tequila.objective.objective import Variable, Variables, ExpectationValue @@ -11,19 +11,20 @@ from tequila.simulators.simulator_api import simulate from tequila.utils import to_float from .chemistry_tools import ActiveSpaceData, FermionicGateImpl, prepare_product_state, ClosedShellAmplitudes, \ - Amplitudes, ParametersQC, NBodyTensor, IntegralManager, OrbitalData + Amplitudes, ParametersQC, NBodyTensor, IntegralManager from .encodings import known_encodings import typing, numpy, numbers from itertools import product -# if you are experiencing import errors you need to update openfermion -# required is version >= 1.0 -# otherwise replace with from openfermion.hamiltonians import MolecularData -import openfermion + try: + # if you are experiencing import errors you need to update openfermion + # required is version >= 1.0 + # otherwise replace with from openfermion.hamiltonians import MolecularData + import openfermion from openfermion.chem import MolecularData except: try: @@ -105,6 +106,13 @@ def __init__(self, parameters: ParametersQC, self._rdm1 = None self._rdm2 = None + def supports_ucc(self): + """ + check if the current molecule supports UCC operations + (e.g. mol.make_excitation_gate) + """ + return self.transformation.supports_ucc + def _initialize_transformation(self, transformation=None, *args, **kwargs): """ Helper Function to initialize the Fermion-to-Qubit Transformation @@ -192,9 +200,8 @@ def make_excitation_generator(self, 1j*Transformed qubit excitation operator, depends on self.transformation """ - if type(self.transformation).__name__ == "BravyiKitaevFast": - raise TequilaException( - "The Bravyi-Kitaev-Superfast transformation does not support general FermionOperators yet") + if not self.supports_ucc(): + raise TequilaException("Molecule with transformation {} does not support general UCC operations".format(self.transformation)) # check indices and convert to list of tuples if necessary if len(indices) == 0: @@ -304,12 +311,15 @@ def make_hardcore_boson_excitation_gate(self, indices, angle, control=None, assu target = [] for pair in indices: assert len(pair) == 2 - target += [pair[0], pair[1]] - consistency = [x < self.n_orbitals for x in target] + target += [self.transformation.up(pair[0]), self.transformation.up(pair[1])] + if self.transformation.up_then_down: + consistency = [x < self.n_orbitals for x in target] + else: + consistency = [x % 2 == 0 for x in target] if not all(consistency): raise TequilaException( - "make_hardcore_boson_excitation_gate: Inconsistencies in indices={}. Should be indexed from 0 ... n_orbitals={}".format( - indices, self.n_orbitals)) + "make_hardcore_boson_excitation_gate: Inconsistencies in indices={} for encoding: {}".format( + indices, self.transformation)) return gates.QubitExcitation(angle=angle, target=target, assume_real=assume_real, control=control, compile_options=compile_options) @@ -404,6 +414,10 @@ def make_excitation_gate(self, indices, angle, control=None, assume_real=True, * Assume that the wavefunction will always stay real. Will reduce potential gradient costs by a factor of 2 """ + + if not self.supports_ucc(): + raise TequilaException("Molecule with transformation {} does not support general UCC operations".format(self.transformation)) + generator = self.make_excitation_generator(indices=indices, remove_constant_term=control is None) p0 = self.make_excitation_generator(indices=indices, form="P0", remove_constant_term=control is None) @@ -556,7 +570,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") + result = QuantumChemistryBase(parameters=self.parameters, integral_manager=integral_manager, orbital_type="native", transformation=self.transformation) return result @@ -645,17 +659,16 @@ def make_hamiltonian(self, *args, **kwargs) -> QubitHamiltonian: qop.is_hermitian() return qop - def make_hardcore_boson_hamiltonian(self): + def make_hardcore_boson_hamiltonian(self, condensed=False): """ Returns ------- Hamiltonian in Hardcore-Boson approximation (electrons are forced into spin-pairs) Indepdent of Fermion-to-Qubit mapping + condensed: always give Hamiltonian back from qubit 0 to N where N is the number of orbitals + if condensed=False then JordanWigner would give back the Hamiltonian defined on even qubits between 0 to 2N """ - if not self.transformation.up_then_down: - warnings.warn( - "Hardcore-Boson Hamiltonian without reordering will result in non-consecutive Hamiltonians that are eventually not be combinable with other features of tequila. Try transformation=\'ReorderedJordanWigner\' or similar for more consistency", - TequilaWarning) + # integrate with QubitEncoding at some point n_orbitals = self.n_orbitals c, obt, tbt = self.get_integrals() @@ -675,6 +688,9 @@ def make_hardcore_boson_hamiltonian(self): uq = q H += h[p, q] * Sm(up) * Sp(uq) + g[p, q] * Sm(up) * Sp(up) * Sm(uq) * Sp(uq) + if not self.transformation.up_then_down and not condensed: + alpha_map = {k.idx:self.transformation.up(k.idx) for k in self.orbitals} + H = H.map_qubits(alpha_map) return H def make_molecular_hamiltonian(self, occupied_indices=None, active_indices=None): @@ -757,11 +773,11 @@ def prepare_hardcore_boson_reference(self): ------- tq.QCircuit that prepares the HCB reference """ - U = gates.X(target=[i.idx for i in self.reference_orbitals]) + U = gates.X(target=[self.transformation.up(i.idx) for i in self.reference_orbitals]) U.n_qubits = self.n_orbitals return U - def hcb_to_me(self, U=None): + def hcb_to_me(self, U=None, condensed=False): """ Transform a circuit in the hardcore-boson encoding (HCB) to the encoding of this molecule @@ -769,6 +785,7 @@ def hcb_to_me(self, U=None): Parameters ---------- U: HCB circuit (using the alpha qubits) + condensed: assume that incoming U is condensed (HCB on the first n_orbitals; and not, as for example in JW on the first n even orbitals) Returns ------- @@ -784,8 +801,12 @@ def hcb_to_me(self, U=None): self.n_orbitals)) # map to alpha qubits - alpha_map = {k: self.transformation.up(k) for k in range(self.n_orbitals)} - alpha_U = U.map_qubits(qubit_map=alpha_map) + if condensed: + alpha_map = {k: self.transformation.up(k) for k in range(self.n_orbitals)} + alpha_U = U.map_qubits(qubit_map=alpha_map) + else: + alpha_U = U + UX = self.transformation.hcb_to_me() if UX is None: raise TequilaException( @@ -931,13 +952,13 @@ def make_spa_ansatz(self, edges, hcb=False, use_units_of_pi=False, label=None, if len(edges) != self.n_electrons//2: raise TequilaException("number of edges need to be equal to number of active electrons//2\n{} edges given\n{} active electrons\nfrozen core is {}".format(len(edges), self.n_electrons, self.parameters.frozen_core)) # making sure that orbitals are uniquely assigned to edges - for edge in edges: - for orbital in edge: + for edge_qubits in edges: + for q1 in edge_qubits: for edge2 in edges: - if edge2==edge: + if edge2==edge_qubits: continue - elif orbital in edge2: - raise TequilaException("make_spa_ansatz: faulty list of edges, orbitals are overlapping e.g. orbital {} is in edge {} and edge {}".format(orbital, edge, edge2)) + elif q1 in edge2: + raise TequilaException("make_spa_ansatz: faulty list of edges, orbitals are overlapping e.g. orbital {} is in edge {} and edge {}".format(q1, edge_qubits, edge2)) # auto assign if the circuit construction is optimized # depending on the current qubit encoding (if hcb_to_me is implemnented we can optimize) @@ -955,55 +976,66 @@ def make_spa_ansatz(self, edges, hcb=False, use_units_of_pi=False, label=None, # construction of the optimized circuit if optimize: - for edge in edges: - U += gates.X(2*edge[0]) - previous = edge[0] - if len(edge)==1: + # circuit in HCB representation + # depends a bit on the ordering of the spin-orbitals in the encoding + # here we transform it to the qubits representing the up-spins + # the hcb_to_me sequence will then transfer to the actual encoding later + for edge_orbitals in edges: + edge_qubits = [self.transformation.up(i) for i in edge_orbitals] + U += gates.X(edge_qubits[0]) + if len(edge_qubits)==1: continue - for orbital in edge[1:]: - c=previous + for i in range(1,len(edge_qubits)): + q1=edge_qubits[i] + c=edge_qubits[i-1] if not ladder: - c=edge[0] - angle=Variable(name=((c, orbital), "D" ,label)) + c=edge_qubits[0] + angle=Variable(name=((edge_orbitals[i-1], edge_orbitals[i]), "D" ,label)) if use_units_of_pi: angle=angle*numpy.pi - if previous == edge[0]: - U += gates.Ry(angle=angle, target=2*orbital, control=None) + if i-1 == 0: + U += gates.Ry(angle=angle, target=q1, control=None) else: - U += gates.Ry(angle=angle, target=2*orbital, control=2*c) - U += gates.CNOT(2*orbital,2*c) - previous = orbital + U += gates.Ry(angle=angle, target=q1, control=c) + U += gates.CNOT(q1,c) + if not hcb: U += self.hcb_to_me() else: - # construction of the non-optimized circuit - U = self.prepare_reference() + # construction of the non-optimized circuit (UpCCD with paired doubles according to edges) + if hcb: + U = self.prepare_hardcore_boson_reference() + else: + U = self.prepare_reference() # will only work if the first orbitals in the edges are the reference orbitals sane = True reference_orbitals = self.reference_orbitals - for edge in edges: - if self.orbitals[edge[0]] not in reference_orbitals: + for edge_qubits in edges: + if self.orbitals[edge_qubits[0]] not in reference_orbitals: sane=False - if len(edge)>1: - for orbital in edge[1:]: - if self.orbitals[orbital] in reference_orbitals: + if len(edge_qubits)>1: + for q1 in edge_qubits[1:]: + if self.orbitals[q1] in reference_orbitals: sane=False if not sane: raise TequilaException("Non-Optimized SPA (e.g. with encodings that are not JW) will only work if the first orbitals of all SPA edges are occupied reference orbitals and all others are not. You gave edges={} and reference_orbitals are {}".format(edges, reference_orbitals)) - for edge in edges: - previous = edge[0] - if len(edge)>1: - for orbital in edge[1:]: + for edge_qubits in edges: + previous = edge_qubits[0] + if len(edge_qubits)>1: + for q1 in edge_qubits[1:]: c = previous if not ladder: - c = edge[0] - angle = Variable(name=((c,orbital), "D" ,label)) + c = edge_qubits[0] + angle = Variable(name=((c,q1), "D" ,label)) if use_units_of_pi: angle=angle*numpy.pi - U += self.make_excitation_gate(indices=[(2*c,2*orbital),(2*c+1,2*orbital+1)], angle=angle) - previous = orbital + if hcb: + U += self.make_hardcore_boson_excitation_gate(indices=[(q1,c)],angle=angle) + else: + U += self.make_excitation_gate(indices=[(2*c,2*q1),(2*c+1,2*q1+1)], angle=angle) + previous = q1 return U def make_ansatz(self, name: str, *args, **kwargs): @@ -1584,7 +1616,7 @@ def rdm2(self): return None def compute_rdms(self, U: QCircuit = None, variables: Variables = None, spin_free: bool = True, - get_rdm1: bool = True, get_rdm2: bool = True, ordering="dirac"): + get_rdm1: bool = True, get_rdm2: bool = True, ordering="dirac", use_hcb: bool = False): """ Computes the one- and two-particle reduced density matrices (rdm1 and rdm2) given a unitary U. This method uses the standard ordering in physics as denoted below. @@ -1619,7 +1651,6 @@ def compute_rdms(self, U: QCircuit = None, variables: Variables = None, spin_fre # Check whether unitary circuit is not 0 if U is None: raise TequilaException('Need to specify a Quantum Circuit.') - # Check whether transformation is BKSF. # Issue here: when a single operator acts only on a subset of qubits, BKSF might not yield the correct # transformation, because it computes the number of qubits incorrectly in this case. @@ -1627,7 +1658,6 @@ def compute_rdms(self, U: QCircuit = None, variables: Variables = None, spin_fre if type(self.transformation).__name__ == "BravyiKitaevFast": raise TequilaException( "The Bravyi-Kitaev-Superfast transformation does not support general FermionOperators yet.") - # Set up number of spin-orbitals and molecular orbitals respectively n_SOs = 2 * self.n_orbitals n_MOs = self.n_orbitals @@ -1636,6 +1666,22 @@ def compute_rdms(self, U: QCircuit = None, variables: Variables = None, spin_fre if U is None: raise TequilaException('Need to specify a Quantum Circuit.') + def _get_hcb_op(op_tuple): + '''Build the hardcore boson operators: b^\dagger_ib_j + h.c. in qubit encoding ''' + if (len(op_tuple) == 2): + return 2 * Sm(op_tuple[0][0]) * Sp(op_tuple[1][0]) + elif (len(op_tuple) == 4): + if ((op_tuple[0][0] == op_tuple[1][0]) and (op_tuple[2][0] == op_tuple[3][0])): # iijj uddu+duud + return Sm(op_tuple[0][0]) * Sp(op_tuple[2][0]) + Sm(op_tuple[2][0]) * Sp(op_tuple[0][0]) + if ((op_tuple[0][0] == op_tuple[2][0]) and (op_tuple[1][0] == op_tuple[3][0]) and ( + op_tuple[0][0] != op_tuple[1][0]) and (op_tuple[2][0] != op_tuple[3][0])): # ijij uuuu+dddd + return 4 * Sm(op_tuple[0][0]) * Sm(op_tuple[1][0]) * Sp(op_tuple[2][0]) * Sp(op_tuple[3][0]) + if ((op_tuple[0][0] == op_tuple[3][0]) and (op_tuple[1][0] == op_tuple[2][0]) and ( + op_tuple[0][0] != op_tuple[1][0]) and (op_tuple[2][0] != op_tuple[3][0])): # ijji abba + return -2 * Sm(op_tuple[0][0]) * Sm(op_tuple[1][0]) * Sp(op_tuple[2][0]) * Sp(op_tuple[3][0]) + else: + return Zero() + def _get_of_op(operator_tuple): """ Returns operator given by a operator tuple as OpenFermion - Fermion operator """ op = openfermion.FermionOperator(operator_tuple) @@ -1718,9 +1764,7 @@ def _build_2bdy_operators_spinfree() -> list: op_tuple = ((2 * p + 1, 1), (2 * q + 1, 1), (2 * s + 1, 0), (2 * r + 1, 0)) if ( p != q and r != s) else '0.0 []' op += _get_of_op(op_tuple) - ops += [op] - return ops def _assemble_rdm1(evals) -> numpy.ndarray: @@ -1783,17 +1827,67 @@ def _assemble_rdm2_spinfree(evals) -> numpy.ndarray: return rdm2 + def _build_1bdy_operators_hcb() -> list: + """ Returns hcb one-body operators as a symmetry-reduced list of QubitHamiltonians """ + # Exploit symmetry pq = qp (not changed by spin-summation) + ops = [] + for p in range(n_MOs): + for q in range(p + 1): + if (p == q): + if (self.transformation.up_then_down): + op_tuple = ((p, 1), (p, 0)) + op = _get_hcb_op(op_tuple) + else: + op_tuple = ((2 * p, 1), (2 * p, 0)) + op = _get_hcb_op(op_tuple) + ops += [op] + else: + ops += [Zero()] + return ops + + def _build_2bdy_operators_hcb() -> list: + """ Returns hcb two-body operators as a symmetry-reduced list of QubitHamiltonians """ + # Exploit symmetries pqrs = qpsr (due to spin summation, '-pqsr = -qprs' drops out) + # and = rspq + ops = [] + scale = 2 + if self.transformation.up_then_down: + scale = 1 + for p, q, r, s in product(range(n_MOs), repeat=4): + if p * n_MOs + q >= r * n_MOs + s and (p >= q or r >= s): + # Spin abba+ baab allow p=q=r=s orb iijj + op_tuple = ((scale * p, 1), (scale * q, 1), (scale * r, 0), (scale * s, 0)) if ( + p == q and s == r) else '0.0 []' + op = _get_hcb_op(op_tuple) + # Spin abba+ baab dont allow p=q=r=s orb ijij + op_tuple = ((scale * p, 1), (scale * q, 1), (scale * r, 0), (scale * s, 0)) if ( + p != q and r != s and p == r and s == q) else '0.0 []' + op += _get_hcb_op(op_tuple) + # Spin aaaa+ bbbb dont allow p=q=r=s orb ijji + op_tuple = ((scale * p, 1), (scale * q, 1), (scale * r, 0), (scale * s, 0)) if ( + p != q and r != s and p == s and q == r) else '0.0 []' + op += _get_hcb_op(op_tuple) + ops += [op] + return ops + # Build operator lists qops = [] - if spin_free: + if spin_free and not use_hcb: qops += _build_1bdy_operators_spinfree() if get_rdm1 else [] qops += _build_2bdy_operators_spinfree() if get_rdm2 else [] + elif use_hcb: + qops += _build_1bdy_operators_hcb() if get_rdm1 else [] + qops += _build_2bdy_operators_hcb() if get_rdm2 else [] else: + if use_hcb: + raise TequilaException( + "compute_rdms: spin_free={} and use_hcb={} are not compatible".format(spin_free, use_hcb)) qops += _build_1bdy_operators_spinful() if get_rdm1 else [] qops += _build_2bdy_operators_spinful() if get_rdm2 else [] # Transform operator lists to QubitHamiltonians - qops = [_get_qop_hermitian(op) for op in qops] + if (not use_hcb): + qops = [_get_qop_hermitian(op) for op in qops] # Compute expected values evals = simulate(ExpectationValue(H=qops, U=U, shape=[len(qops)]), variables=variables) @@ -1801,7 +1895,7 @@ def _assemble_rdm2_spinfree(evals) -> numpy.ndarray: # If self._rdm1, self._rdm2 exist, reset them if they are of the other spin-type def _reset_rdm(rdm): if rdm is not None: - if spin_free and rdm.shape[0] != n_MOs: + if (spin_free or use_hcb) and rdm.shape[0] != n_MOs: return None if not spin_free and rdm.shape[0] != n_SOs: return None @@ -1811,13 +1905,13 @@ def _reset_rdm(rdm): self._rdm2 = _reset_rdm(self._rdm2) # Split expectation values in 1- and 2-particle expectation values if get_rdm1: - len_1 = n_MOs * (n_MOs + 1) // 2 if spin_free else n_SOs * (n_SOs + 1) // 2 + len_1 = n_MOs * (n_MOs + 1) // 2 if (spin_free or use_hcb) else n_SOs * (n_SOs + 1) // 2 else: len_1 = 0 evals_1, evals_2 = evals[:len_1], evals[len_1:] # Build matrices using the expectation values self._rdm1 = _assemble_rdm1(evals_1) if get_rdm1 else self._rdm1 - if spin_free: + if spin_free or use_hcb: self._rdm2 = _assemble_rdm2_spinfree(evals_2) if get_rdm2 else self._rdm2 else: self._rdm2 = _assemble_rdm2_spinful(evals_2) if get_rdm2 else self._rdm2 diff --git a/src/tequila/simulators/simulator_api.py b/src/tequila/simulators/simulator_api.py index 163d3fe1..63506ca6 100755 --- a/src/tequila/simulators/simulator_api.py +++ b/src/tequila/simulators/simulator_api.py @@ -11,7 +11,7 @@ from tequila.circuit.noise import NoiseModel SUPPORTED_BACKENDS = ["qulacs_gpu", "qulacs",'qibo', "qiskit", "cirq", "pyquil", "symbolic", "qlm"] -SUPPORTED_NOISE_BACKENDS = ["qiskit",'qibo', 'cirq', 'pyquil', 'qulacs', "qulacs_gpu"] +SUPPORTED_NOISE_BACKENDS = ["qiskit", 'cirq', 'pyquil'] # qulacs removed in v.1.9 BackendTypes = namedtuple('BackendTypes', 'CircType ExpValueType') INSTALLED_SIMULATORS = {} INSTALLED_SAMPLERS = {} diff --git a/src/tequila/simulators/simulator_qibo.py b/src/tequila/simulators/simulator_qibo.py index 7d7e7e7d..027b22cb 100755 --- a/src/tequila/simulators/simulator_qibo.py +++ b/src/tequila/simulators/simulator_qibo.py @@ -355,7 +355,7 @@ def do_simulate(self, variables, initial_state=None, *args, **kwargs): """ n_qubits = max(self.highest_qubit + 1, self.n_qubits, self.abstract_circuit.max_qubit() + 1) if initial_state is not None: - if isinstance(initial_state, int) or isinstance(initial_state,np.int): + if isinstance(initial_state, (int, np.int64)): wave = QubitWaveFunction.from_int(i=initial_state, n_qubits=n_qubits) elif isinstance(initial_state, str): wave = QubitWaveFunction.from_string(string=initial_state).to_array() diff --git a/src/tequila/simulators/simulator_qulacs.py b/src/tequila/simulators/simulator_qulacs.py index 204d5fca..29cfca95 100755 --- a/src/tequila/simulators/simulator_qulacs.py +++ b/src/tequila/simulators/simulator_qulacs.py @@ -1,6 +1,8 @@ import qulacs import numbers, numpy -from tequila import TequilaException +import warnings + +from tequila import TequilaException, TequilaWarning from tequila.utils.bitstrings import BitNumbering, BitString, BitStringLSB from tequila.wavefunction.qubit_wavefunction import QubitWaveFunction from tequila.simulators.simulator_base import BackendCircuit, BackendExpectationValue, QCircuit, change_basis @@ -91,6 +93,9 @@ def __init__(self, abstract_circuit, noise=None, *args, **kwargs): super().__init__(abstract_circuit=abstract_circuit, noise=noise, *args, **kwargs) self.has_noise=False if noise is not None: + + warnings.warn("Warning: noise in qulacs module will be dropped. Currently only works for qulacs version 0.5 or lower", TequilaWarning) + self.has_noise=True self.noise_lookup = { 'bit flip': [qulacs.gate.BitFlipNoise], diff --git a/src/tequila/version.py b/src/tequila/version.py index d5bb692d..3b2e275a 100644 --- a/src/tequila/version.py +++ b/src/tequila/version.py @@ -1,2 +1,2 @@ -__version__ = "1.8.8" +__version__ = "1.8.9" __author__ = "Jakob Kottmann and Sumner Alperin-Lea and Teresa Tamayo-Mendoza and Alba Cervera-Lierta and Cyrille Lavigne and Tzu-Ching Yen and Vladyslav Verteletskyi and Philipp Schleich and Abhinav Anand and Matthias Degroote and Skylar Chaney and Maha Kesibi and Naomi Grace Curnow and Brandon Solo and Georgios Tsilimigkounakis and Claudia Zendejas-Morales and Artur F Izmaylov and Alan Aspuru-Guzik ... " diff --git a/tests/test_binary_pauli.py b/tests/test_binary_pauli.py index cfacf155..499fe43f 100644 --- a/tests/test_binary_pauli.py +++ b/tests/test_binary_pauli.py @@ -321,7 +321,7 @@ def test_get_qubit_wise(): assert np.isclose(result_ori, result_integrated_si) # Checking the optimized expectation values are the same - initial_values = {k: np.random.uniform(0.0, 6.0, 1) for k in e_ori.extract_variables()} + initial_values = {k: np.random.uniform(0.0, 6.0) for k in e_ori.extract_variables()} sol1 = tq.minimize(method='bfgs', objective=e_ori, initial_values=initial_values) sol2 = tq.minimize(method='bfgs', objective=e_qwc, initial_values=initial_values) sol3 = tq.minimize(method='bfgs', objective=e_integrated, initial_values=initial_values) diff --git a/tests/test_chemistry.py b/tests/test_chemistry.py index 66831af6..86b2dc33 100644 --- a/tests/test_chemistry.py +++ b/tests/test_chemistry.py @@ -30,7 +30,7 @@ def teardown_function(function): [os.remove(x) for x in glob.glob("qvm.log")] [os.remove(x) for x in glob.glob("*.dat")] -@pytest.mark.skipif(condition=not HAS_PSI4 or not HAS_PYSCF, reason="you don't have psi4 or pyscf") +@pytest.mark.skipif(condition=not HAS_PSI4 and not HAS_PYSCF, reason="you don't have psi4 or pyscf") def test_UR_and_UC(): mol = tq.Molecule(geometry="h 0.0 0.0 0.0\nH 0.0 0.0 1.5", basis_set="sto-3g") mol = mol.use_native_orbitals() @@ -64,7 +64,7 @@ def test_base(trafo): assert len(eigvals) == 16 -@pytest.mark.skipif(condition=not HAS_PSI4 or not HAS_PYSCF, reason="you don't have psi4 or pyscf") +@pytest.mark.skipif(condition=not HAS_PSI4 and not HAS_PYSCF, reason="you don't have psi4 or pyscf") @pytest.mark.parametrize("trafo", ["JordanWigner", "BravyiKitaev", "BravyiKitaevTree"]) def test_prepare_reference(trafo): geometry = "Li 0.0 0.0 0.0\nH 0.0 0.0 1.5" @@ -83,7 +83,7 @@ def test_prepare_reference(trafo): energy2 = tq.simulate(E) assert numpy.isclose(energy, energy2, atol=1.e-4) -@pytest.mark.skipif(condition=not HAS_PSI4 or not HAS_PYSCF, reason="you don't have psi4 or pyscf") +@pytest.mark.skipif(condition=not HAS_PSI4 and not HAS_PYSCF, reason="you don't have psi4 or pyscf") def test_orbital_types(): geometry = "H 0.0 0.0 0.0\nH 0.0 0.0 2.0\nH 0.0 0.0 4.0\nH 0.0 0.0 6.0" @@ -124,15 +124,18 @@ def test_dependencies(): assert key in qc.INSTALLED_QCHEMISTRY_BACKENDS.keys() -@pytest.mark.skipif(condition=not HAS_PSI4 or not HAS_PYSCF, reason="no quantum chemistry backends installed") +@pytest.mark.skipif(condition=not HAS_PSI4 and not HAS_PYSCF, reason="no quantum chemistry backends installed") def test_interface(): molecule = tq.chemistry.Molecule(basis_set='sto-3g', geometry="data/h2.xyz", transformation="JordanWigner") -@pytest.mark.skipif(condition=not HAS_PSI4 or not HAS_PYSCF, reason="you don't have psi4") +@pytest.mark.skipif(condition=not HAS_PSI4, reason="you don't have psi4") def test_h2_hamiltonian_psi4(): do_test_h2_hamiltonian(qc_interface=qc.QuantumChemistryPsi4) +@pytest.mark.skipif(condition=not HAS_PYSCF, reason="you don't have pyscf") +def test_h2_hamiltonian_psi4(): + do_test_h2_hamiltonian(qc_interface=qc.QuantumChemistryPySCF) def do_test_h2_hamiltonian(qc_interface): parameters = tequila.quantumchemistry.chemistry_tools.ParametersQC(geometry="data/h2.xyz", basis_set="sto-3g") @@ -277,7 +280,7 @@ def test_restart_psi4(): @pytest.mark.skipif(condition=not HAS_PSI4, reason="psi4 not found") @pytest.mark.parametrize("active", [{"A1": [2, 3]}, {"B2": [0], "B1": [0]}, {"A1": [0, 1, 2, 3]}, {"B1": [0]}]) -def test_active_spaces(active): +def test_psi4_active_spaces(active): mol = tq.chemistry.Molecule(geometry="data/h2o.xyz", basis_set="sto-3g", active_orbitals=active) H = mol.make_hamiltonian() Uhf = mol.prepare_reference() @@ -301,12 +304,13 @@ def test_rdms_psi4(): assert (numpy.allclose(rdm2, rdm2_ref, atol=1e-8)) -@pytest.mark.skipif(condition=not HAS_PSI4 or not HAS_PYSCF, reason="quantum chemistry backend not found") +@pytest.mark.skipif(condition=not HAS_PSI4 and not HAS_PYSCF, reason="psi4/pyscf not found") @pytest.mark.parametrize("geometry", ["H 0.0 0.0 0.0\nH 0.0 0.0 0.7"]) -@pytest.mark.parametrize("trafo", ["JordanWigner", "BravyiKitaev", "BravyiKitaevTree", "ReorderedJordanWigner", - "ReorderedBravyiKitaev"]) +@pytest.mark.parametrize("trafo", tq.quantumchemistry.encodings.known_encodings()) def test_upccgsd(geometry, trafo): molecule = tq.chemistry.Molecule(geometry=geometry, basis_set="sto-3g", transformation=trafo) + if not molecule.supports_ucc(): + return energy = do_test_upccgsd(molecule) fci = molecule.compute_energy("fci") assert numpy.isclose(fci, energy, atol=1.e-3) @@ -314,7 +318,7 @@ def test_upccgsd(geometry, trafo): assert numpy.isclose(fci, energy2, atol=1.e-3) -@pytest.mark.skipif(condition=not HAS_PSI4 or not HAS_PYSCF, reason="psi4 or pyscf not found") +@pytest.mark.skipif(condition=not HAS_PSI4 and not HAS_PYSCF, reason="psi4 or pyscf not found") def test_upccgsd_singles(): molecule = tq.chemistry.Molecule(geometry="H 0.0 0.0 0.0\nH 0.0 0.0 0.7", basis_set="6-31G") H = molecule.make_hamiltonian() @@ -342,16 +346,16 @@ def do_test_upccgsd(molecule, *args, **kwargs): print(E2.extract_variables()) energy1 = tq.simulate(E, variables=variables) energy2 = tq.simulate(E2, variables=variables2) - assert energy1 == energy2 + assert numpy.isclose(energy1,energy2,atol=1.e-4) return result.energy @pytest.mark.parametrize("backend", tq.simulators.simulator_api.INSTALLED_SIMULATORS.keys()) -@pytest.mark.skipif(condition=not HAS_PSI4 or not HAS_PYSCF, reason="psi4/pyscf not found") +@pytest.mark.skipif(condition=not HAS_PSI4 and not HAS_PYSCF, reason="psi4/pyscf not found") def test_hamiltonian_reduction(backend): mol = tq.chemistry.Molecule(geometry="H 0.0 0.0 0.0\nH 0.0 0.0 0.7", basis_set="6-31G") - hf = mol.energies["hf"] + hf = mol.compute_energy("hf") U = mol.prepare_reference() H = mol.make_hamiltonian() E = tq.simulate(tq.ExpectationValue(H=H, U=U), backend=backend) @@ -362,7 +366,7 @@ def test_hamiltonian_reduction(backend): assert numpy.isclose(E, hf, atol=1.e-4) -@pytest.mark.skipif(condition=not HAS_PSI4 or not HAS_PYSCF, reason="psi4/pyscf not found") +@pytest.mark.skipif(condition=not HAS_PSI4 and not HAS_PYSCF, reason="psi4/pyscf not found") @pytest.mark.parametrize("assume_real", [True, False]) @pytest.mark.parametrize("trafo", ["jordan_wigner", "bravyi_kitaev", "tapered_bravyi_kitaev"]) def test_fermionic_gates(assume_real, trafo): @@ -404,7 +408,7 @@ def test_fermionic_gates(assume_real, trafo): assert numpy.isclose(test2, test2x, atol=1.e-6) -@pytest.mark.skipif(condition=not HAS_PSI4 or not HAS_PYSCF, reason="psi4/pyscf not found") +@pytest.mark.skipif(condition=not HAS_PSI4 and not HAS_PYSCF, reason="psi4/pyscf not found") @pytest.mark.parametrize("trafo", ["JordanWigner", "BravyiKitaev", "BravyiKitaevTree", "ReorderedJordanWigner", "ReorderedBravyiKitaev"]) def test_hcb(trafo): @@ -452,7 +456,6 @@ def test_pyscf_methods(method, geometry, basis_set): @pytest.mark.skipif(condition=not HAS_PYSCF, reason="pyscf not found") -@pytest.mark.skipif(condition=not HAS_PSI4, reason="psi4 not found") def test_orbital_optimization(): from tequila.quantumchemistry import optimize_orbitals mol = tq.Molecule(geometry="Li 0.0 0.0 0.0\nH 0.0 0.0 3.0", basis_set="STO-3G") @@ -467,7 +470,6 @@ def test_orbital_optimization(): assert numpy.isclose(-7.79860454, result.energy, atol=1.e-3) -@pytest.mark.skipif(Qcondition=not HAS_PSI4, reason="psi4 not found") @pytest.mark.skipif(condition=not HAS_PYSCF, reason="pyscf not found") def test_orbital_transformation(): mol0 = tq.Molecule(geometry="Li 0.0 0.0 0.0\nH 0.0 0.0 0.75", basis_set="STO-3G", frozen_core=False) @@ -547,8 +549,12 @@ def test_crosscheck_cis_mp2_large(): assert numpy.isclose(e_1, e_2, atol=1.e-4) -@pytest.mark.skipif(condition=not HAS_PSI4 or not HAS_PYSCF, reason="psi4/pyscf not found") +@pytest.mark.skipif(condition=not HAS_PSI4 and not HAS_PYSCF, reason="psi4/pyscf not found") def test_spa_ansatz_be(): + print("\n\n") + print(HAS_PSI4) + print(HAS_PYSCF) + print(not HAS_PSI4 and not HAS_PYSCF) edges = [(0,), (1, 2, 3, 4)] # doing without frozen-core to test explicitly if single-orbital pairs work mol = tq.Molecule(geometry="be 0.0 0.0 0.0", basis_set="sto-3g", transformation="BravyiKitaev", frozen_core=False) @@ -569,3 +575,153 @@ def test_spa_ansatz_be(): E = tq.ExpectationValue(H=H, U=U) result = tq.minimize(E, silent=True) assert numpy.isclose(energy, result.energy) + print("a") + +#@pytest.mark.skipif(condition=not HAS_PSI4 and not HAS_PYSCF, reason="psi4/pyscf not found") +@pytest.mark.parametrize("name", ["SPA", "HCB-SPA"]) +@pytest.mark.parametrize("optimize", [True, False]) +@pytest.mark.parametrize("geometry", ["H 0.0 0.0 0.0\nH 0.0 0.0 4.5", "Li 0.0 0.0 0.0\nH 0.0 0.0 3.0", "Be 0.0 0.0 0.0\nH 0.0 0.0 3.0\nH 0.0 0.0 -3.0"]) +@pytest.mark.parametrize("transformation", tq.quantumchemistry.encodings.known_encodings()) +@pytest.mark.skipif(condition=not HAS_PSI4 and not HAS_PYSCF, reason="psi4/pyscf not found") +def test_spa_consistency(geometry, name, optimize, transformation): + mol = tq.Molecule(geometry=geometry, basis_set="sto-3g", + transformation=transformation).use_native_orbitals() + + if mol.transformation.hcb_to_me() is None: + return + # doesn't need to make physical sense for the test + edges = [] + i = 0 + for i in range(mol.n_electrons//2): + a = i+mol.n_electrons//2 + edge = (i,a) + edges.append(edge) + U = mol.make_ansatz("SPA", edges=edges, optimize=optimize) + variables = {k: 1.0 for k in U.extract_variables()} + wfn = (tq.simulate(U, variables=variables)) + + UX = mol.make_ansatz(name=name, edges=edges, optimize=optimize) + wfnx = (tq.simulate(UX, variables=variables)) + print(wfnx) + if "HCB" in name: + UX += mol.hcb_to_me() + wfnx = (tq.simulate(UX, variables=variables)) + F = abs(wfn.inner(wfnx)) ** 2 + if not numpy.isclose(F, 1.0): + print(wfn) + print(wfnx) + print(UX) + + assert numpy.isclose(F, 1.0) + +@pytest.mark.skipif(condition=not HAS_PSI4 and not HAS_PYSCF, reason="psi4/pyscf not found") +def test_variable_consistency(): + geometry = "H 0.0 0.0 0.0\nH 0.0 0.0 1.0\nH 0.0 0.0 2.0\nH 0.0 0.0 3.0" + mol = tq.Molecule(geometry=geometry, basis_set="sto-3g", transformation="ReorderedJordanWigner").use_native_orbitals() + U = mol.make_ansatz("SPA",edges=[(0,1),(2,3)]) + variables = {k:1.0 for k in U.extract_variables()} + + H = mol.make_hamiltonian() + E = tq.ExpectationValue(H=H, U=U) + E1 = tq.simulate(E, variables=variables) + + U = mol.make_ansatz("HCB-SPA",edges=[(0,1),(2,3)]) + H = mol.make_hardcore_boson_hamiltonian(condensed=False) + E = tq.ExpectationValue(H=H, U=U) + E2 = tq.simulate(E, variables=variables) + + + mol = tq.Molecule(geometry=geometry, basis_set="sto-3g", transformation="JordanWigner").use_native_orbitals() + U = mol.make_ansatz("SPA",edges=[(0,1),(2,3)]) + H = mol.make_hamiltonian(condensed=False) + E = tq.ExpectationValue(H=H, U=U) + E3 = tq.simulate(E, variables=variables) + + U = mol.make_ansatz("HCB-SPA",edges=[(0,1),(2,3)]) + H = mol.make_hardcore_boson_hamiltonian(condensed=False) + E = tq.ExpectationValue(H=H, U=U) + E4 = tq.simulate(E, variables=variables) + + assert numpy.isclose(E1,E2,atol=1.e-4) + assert numpy.isclose(E1,E3,atol=1.e-4) + assert numpy.isclose(E1,E4,atol=1.e-4) + assert numpy.isclose(E2,E3,atol=1.e-4) + assert numpy.isclose(E2,E4,atol=1.e-4) + +@pytest.mark.skipif(condition=not HAS_PSI4 and not HAS_PYSCF, reason="you don't have psi4 or pyscf") +@pytest.mark.parametrize("geometry", ["H 0.0 0.0 0.0\nH 0.0 0.0 4.5", "Li 0.0 0.0 0.0\nH 0.0 0.0 3.0", "Be 0.0 0.0 0.0\nH 0.0 0.0 3.0\nH 0.0 0.0 -3.0"]) +@pytest.mark.parametrize("optimize", [True, False]) +def test_hcb_rdms(geometry, optimize): + + mol1 = tq.Molecule(geometry=geometry, basis_set="sto-3g", transformation="ReorderedJordanWigner") + H = mol1.make_hamiltonian() + HCB = mol1.make_hardcore_boson_hamiltonian() + + if mol1.n_electrons == 4: + edges=[(0, 2, 5), (1, 3, 4)] + elif mol1.n_electrons == 2: + edges=[(0, 1)] + else: + raise Exception("test only created for n_electrons = 2,4 you gave {}".format(mol1.n_electrons)) + + + U1_1 = mol1.make_ansatz("SPA", edges=edges, optimize=optimize) + U1_2 = mol1.make_ansatz("HCB-SPA", edges=edges, optimize=optimize) + E1 = tq.ExpectationValue(H=H, U=U1_1) + E2 = tq.ExpectationValue(H=HCB, U=U1_2) + + variables = {k: 1.0 for k in U1_1.extract_variables()} + energy1 = tq.simulate(E1, variables=variables) + energy2 = tq.simulate(E2, variables=variables) + + xrdm1_1, xrdm1_2 = mol1.compute_rdms(U1_1, variables=variables, use_hcb=False) + yrdm1_1, yrdm1_2 = mol1.compute_rdms(U1_2, variables=variables, use_hcb=True) + + c,h1,h2 = mol1.get_integrals() + # this is the default for RDMs + # integrals are given in openfermion ordering + h2 = h2.reorder(to="phys").elems + + energy3 = numpy.einsum('qp, pq', h1,xrdm1_1, optimize='greedy') + energy3+= 1 / 2 * numpy.einsum('rspq, pqrs', h2,xrdm1_2, optimize='greedy') + energy3+= c + + + d1_1 = tq.numpy.linalg.norm(xrdm1_1 - yrdm1_1) + d1_2 = tq.numpy.linalg.norm(xrdm1_2 - yrdm1_2) + + assert numpy.isclose(d1_1, 0., atol=1.e-4) + assert numpy.isclose(d1_2, 0., atol=1.e-4) + + assert numpy.isclose(energy1, energy2, atol=1.e-4) + assert numpy.isclose(energy1, energy3, atol=1.e-4) + + +@pytest.mark.skipif(condition=not HAS_PYSCF, reason="you don't have pyscf") +@pytest.mark.parametrize("geometry", ["H 0.0 0.0 0.0\nH 0.0 0.0 4.5", "Li 0.0 0.0 0.0\nH 0.0 0.0 3.0", "Be 0.0 0.0 0.0\nH 0.0 0.0 3.0\nH 0.0 0.0 -3.0"]) +def test_orbital_optimization_hcb(geometry): + + mol = tq.Molecule(geometry=geometry, basis_set="sto-3g") + + + if mol.n_electrons == 4: + edges=[(0, 2, 5), (1, 3, 4)] + elif mol.n_electrons == 2: + edges=[(0, 1)] + else: + raise Exception("test only created for n_electrons = 2,4 you gave {}".format(mol.n_electrons)) + + + U1 = mol.make_ansatz(name="HCB-SPA", edges=edges) + import time + start = time.time() + opt1 = tq.chemistry.optimize_orbitals(circuit=U1, molecule=mol, silent=True, use_hcb=True) + time1 = time.time()-start + U2 = mol.make_ansatz(name="SPA", edges=edges) + start = time.time() + opt2 = tq.chemistry.optimize_orbitals(circuit=U2, molecule=mol, silent=True) + time2 = time.time()-start + + assert numpy.isclose(opt1.energy,opt2.energy,atol=1.e-5) + assert time1 < time2 + assert (numpy.isclose(opt1.mo_coeff,opt2.mo_coeff, atol=1.e-5)).all() diff --git a/tests/test_chemistry_f12.py b/tests/test_chemistry_f12.py index c7c9b9b9..cea2d923 100644 --- a/tests/test_chemistry_f12.py +++ b/tests/test_chemistry_f12.py @@ -40,10 +40,13 @@ def test_correction_psi4_cabsplus(): @pytest.mark.skipif(not( os.path.isfile(str(data_path)+'he-f12_gtensor.npy') and os.path.isfile(str(data_path)+'he-f12_htensor.npy') and os.path.isfile(str(data_path)+'he-f12_f12tensor.npy')), reason="data not there") -def test_correction_madness(): +@pytest.mark.parametrize("trafo", tq.quantumchemistry.encodings.known_encodings()) +def test_correction_madness(trafo): geomstring = "He 0.0 0.0 0.0" - mol = tq.Molecule(name='data/f12/he-f12',geometry=geomstring, basis_set='madness', n_pno=None, - active_orbitals=[0, 1], madness_root_dir=None, keep_mad_files=True) + mol = tq.Molecule(name='data/f12/he-f12',geometry=geomstring, basis_set='madness', n_pno="read", + active_orbitals=[0, 1], madness_root_dir=None, keep_mad_files=trafo) + if not mol.supports_ucc(): + return U = mol.make_upccgsd_ansatz(name='UpCCD') angles = {(((0, 1),), 'D', (None, 0)): -0.1278860185100716} rdminfo = {"U": U, "variables": angles}