DMET calculation problem

[commco789]

---

I successfully used DMET-VQE to calculate the energy of the H₁₀ system. However, when I tried using the same pipeline to compute the energy of the BF₄⁻, the simulation crashed with a TypeError: '<' not supported between instances of 'complex' and 'float'.
Here's the minimal code I used:

---

from tangelo import SecondQuantizedMolecule
from tangelo.problem_decomposition import DMETProblemDecomposition
from tangelo.algorithms import VQESolver
bf4="""
F -0.67434965686117 -0.13772650470120 -1.22197648248581
F 1.28027646840210 0.52623822634537 -0.22824137456014
F -0.71946117361893 0.86667506771115 0.83534222327770
F 0.11359664918233 -1.25541074050712 0.61497359584763
B -0.00005885050433 0.00022520685180 -0.00009728207940
"""

mol_bf4 = SecondQuantizedMolecule(bf4, q=-1, spin=0, basis="minao")

options_mol_dmet = {"molecule": mol_bf4,
"fragment_atoms": [1]*5,
"fragment_solvers": "vqe",
"fragment_frozen_orbitals": [[0] for _ in range(5)],
"virtual_orbital_threshold": 1e-10,
"verbose": False
}

dmet_mol = DMETProblemDecomposition(options_mol_dmet)
dmet_mol.build()

dmet_mol.verbose = False
energy_mol_dmet = dmet_mol.simulate()

print(f"DMET energy (hartree): \t {energy_mol_dmet}")

---

and the error message:

---

TypeError Traceback (most recent call last)
Cell In[29], line 2
1 dmet_mol.verbose = False
----> 2 energy_mol_dmet = dmet_mol.simulate()
4 print(f"DMET energy (hartree): \t {energy_mol_dmet}")

File ~/miniforge3/envs/tangelo-env/lib/python3.10/site-packages/tangelo/problem_decomposition/dmet/dmet_problem_decomposition.py:280, in DMETProblemDecomposition.simulate(self)
277 if not self.orbitals:
278 raise RuntimeError("No fragment built. Have you called DMET.build ?")
--> 280 self.chemical_potential = self.optimizer(self._oneshot_loop, self.initial_chemical_potential)
281 self.chemical_potential = self.chemical_potential.real
283 # run one more time to save results

File ~/miniforge3/envs/tangelo-env/lib/python3.10/site-packages/tangelo/problem_decomposition/dmet/dmet_problem_decomposition.py:740, in DMETProblemDecomposition._default_optimizer(self, func, var_params)
726 def _default_optimizer(self, func, var_params):
727 """Function used as a default optimizer for DMET when user does not
728 provide one.
729
(...)
737 float: The chemical potential found by the optimizer.
738 """
--> 740 result = scipy.optimize.newton(func, var_params, tol=1e-5)
742 return result.real

File ~/miniforge3/envs/tangelo-env/lib/python3.10/site-packages/scipy/optimize/_zeros_py.py:385, in newton(func, x0, fprime, args, tol, maxiter, fprime2, x1, rtol, full_output, disp)
383 p0, q0 = p1, q1
384 p1 = p
--> 385 q1 = func(p1, *args)
386 funcalls += 1
388 if disp:

File ~/miniforge3/envs/tangelo-env/lib/python3.10/site-packages/tangelo/problem_decomposition/dmet/dmet_problem_decomposition.py:520, in DMETProblemDecomposition._oneshot_loop(self, chemical_potential, save_results, resample, n_shots, purify, rdm_measurements)
518 solver_fragment = VQESolver({**system, **solver_options})
519 solver_fragment.build()
--> 520 solver_fragment.simulate()
522 if purify and solver_fragment.molecule.n_active_electrons == 2 and not self.uhf:
523 onerdm, twordm = solver_fragment.get_rdm(solver_fragment.optimal_var_params, resample=resample, sum_spin=False)

File ~/miniforge3/envs/tangelo-env/lib/python3.10/site-packages/tangelo/algorithms/variational/vqe_solver.py:259, in VQESolver.simulate(self)
256 if len(self.ansatz.circuit._variational_gates) == 0:
257 raise RuntimeError("No variational gate found in the circuit.")
--> 259 optimal_energy, optimal_var_params = self.optimizer(self.energy_estimation, self.initial_var_params)
261 self.optimal_var_params = optimal_var_params
262 self.optimal_energy = optimal_energy

File ~/miniforge3/envs/tangelo-env/lib/python3.10/site-packages/tangelo/algorithms/variational/vqe_solver.py:709, in VQESolver._default_optimizer(self, func, var_params)
706 from scipy.optimize import minimize
708 with HiddenPrints() if not self.verbose else nullcontext():
--> 709 result = minimize(func, var_params, method="SLSQP",
710 options={"disp": True, "maxiter": 2000, "eps": 1e-5, "ftol": 1e-5})
712 if self.verbose:
713 print(f"VQESolver optimization results:")

File ~/miniforge3/envs/tangelo-env/lib/python3.10/site-packages/scipy/optimize/_minimize.py:750, in minimize(fun, x0, args, method, jac, hess, hessp, bounds, constraints, tol, callback, options)
747 res = _minimize_cobyqa(fun, x0, args, bounds, constraints, callback,
748 **options)
749 elif meth == 'slsqp':
--> 750 res = _minimize_slsqp(fun, x0, args, jac, bounds,
751 constraints, callback=callback, **options)
752 elif meth == 'trust-constr':
753 res = _minimize_trustregion_constr(fun, x0, args, jac, hess, hessp,
754 bounds, constraints,
755 callback=callback, **options)

File ~/miniforge3/envs/tangelo-env/lib/python3.10/site-packages/scipy/optimize/_slsqp_py.py:381, in _minimize_slsqp(func, x0, args, jac, bounds, constraints, maxiter, ftol, iprint, disp, eps, callback, finite_diff_rel_step, **unknown_options)
378 xu[infbnd[:, 1]] = np.nan
380 # ScalarFunction provides function and gradient evaluation
--> 381 sf = _prepare_scalar_function(func, x, jac=jac, args=args, epsilon=eps,
382 finite_diff_rel_step=finite_diff_rel_step,
383 bounds=new_bounds)
384 # gh11403 SLSQP sometimes exceeds bounds by 1 or 2 ULP, make sure this
385 # doesn't get sent to the func/grad evaluator.
386 wrapped_fun = _clip_x_for_func(sf.fun, new_bounds)

File ~/miniforge3/envs/tangelo-env/lib/python3.10/site-packages/scipy/optimize/_optimize.py:291, in _prepare_scalar_function(fun, x0, jac, args, bounds, epsilon, finite_diff_rel_step, hess)
287 bounds = (-np.inf, np.inf)
289 # ScalarFunction caches. Reuse of fun(x) during grad
290 # calculation reduces overall function evaluations.
--> 291 sf = ScalarFunction(fun, x0, args, grad, hess,
292 finite_diff_rel_step, bounds, epsilon=epsilon)
294 return sf

File ~/miniforge3/envs/tangelo-env/lib/python3.10/site-packages/scipy/optimize/_differentiable_functions.py:223, in ScalarFunction.init(self, fun, x0, args, grad, hess, finite_diff_rel_step, finite_diff_bounds, epsilon)
220 finite_diff_options["as_linear_operator"] = True
222 # Initial function evaluation
--> 223 self._update_fun()
225 # Initial gradient evaluation
226 self._wrapped_grad, self._ngev = _wrapper_grad(
227 grad,
228 fun=self._wrapped_fun,
229 args=args,
230 finite_diff_options=finite_diff_options
231 )

File ~/miniforge3/envs/tangelo-env/lib/python3.10/site-packages/scipy/optimize/_differentiable_functions.py:296, in ScalarFunction._update_fun(self)
294 if not self.f_updated:
295 fx = self._wrapped_fun(self.x)
--> 296 if fx < self._lowest_f:
297 self._lowest_x = self.x
298 self._lowest_f = fx

TypeError: '<' not supported between instances of 'complex' and 'float'

---

I don't know how to deal with the situation. It seems to trigger a complex-valued energy or objective, which causes a failure in scipy.optimize.minimize. Any suggestions on how to fix or work around this would be greatly appreciated. Thanks for your excellent work on this project!

Run Environment:

Machine: Mac, M2

Python version: 3.10.18

Tangelo version: latest (installed via pip)

[alexfleury-sb]

Hello @commco789, thanks for raising this issue!

I was able to replicate the error. This is due to the circuit backend producing complex values when computing the Hamiltonian eigenvalues, which ideally should be real. This is likely due to small numerical errors, where a tiny imaginary part is inadvertently added. Cirq is the default backend for Tangelo in a vanilla environment; however, it tends to produce imaginary values. There are 2 ways to circumvent this situation.

1. Switch the backend. The easiest way is to install Qulacs in your environment. If Qulacs is installed, Tangelo switches its default backend to it, as we developers prefer it for its speed.

This can be done with this command in your environment:

pip install qulacs

I found that the error does not appear when using your example. However, I did not test to completion, so I cannot guarantee the error won't emerge later in the calculation.

2. A more robust way would be to patch the package to force real energies in VQESolver. This can be done by changing this part of the code

Tangelo/tangelo/algorithms/variational/vqe_solver.py

Line 322 in 480a8b1

return energy

to something like

return energy.real

This would be a more robust solution. I recommend this as a primary option if you are comfortable patching code. If not, please let us know, and we can update the official GitHub package.

[ValentinS4t1qbit]

Thank you both ! I recall enforcing this to be real in very early versions of Tangelo, but relaxed this.

=> Do we have any reason to believe that there are cases where end users may run VQE while playing with numbers with a significant imaginary part (e.g clearly non-zero, not because of small rounding errors) ? For instance, if an end user specifies an input Hamiltonian that they found / computed somewhere. They would then likely get "weird results".

What about:

• check whether energy.imag > 0.1, print a warning if it is
• return energy.real regardless

Or do we feel strongly that VQE should enforce real values and everything else is the user's responsibility ?

[alexfleury-sb]

Printing a warning could be a good approach. The function is called energy_estimation, so we should expect reals for Hamiltonians.

[commco789]

@alexfleury-sb @ValentinS4t1qbit Thank you both for your insightful replies! I’ve tried using Qulacs as the backend, and there are no more errors! However, I’ve noticed that the computation time is still quite long for such a simple system(BF₄⁻). Even with a fragment containing just a single atom, the calculation has been running for several days using 8 cores, and it still hasn't finished. Is there any way to improve the computation speed?

[alexfleury-sb]

Glad to hear you've unblocked your work for now! Regarding the speed, it looks like your first fragment is spanning 16 qubits. For a VQE on a system this large, you'll have a significant number of parameters to optimize, requiring many energy evaluations. These factors likely contribute to the intensive calculation.

Before diving into quantum optimizations, I recommend verifying your DMET fragmentation by using a classical solver like FCI or CCSD on the fragment. This will confirm the correctness of your setup.

Once that's confirmed, here are some strategies to explore for improving performance:

• Change the qubit mapping: Switch to scbk. This can immediately reduce the number of qubits by two without losing information.
• Try a different ansatz: The default UCCSD ansatz can be computationally demanding. Consider using a Hardware-Efficient ansatz instead. While it might come with an accuracy penalty, it's generally much easier to optimize.
• Reduce the active space size: You can further reduce the active space using the fragment_frozen_orbitals argument.

You can modify all these parameters using the solver_options in your DMETProblemDecomposition object.

[commco789]

Thank you so much for your helpful response! As I'm still new to Tangelo, I really appreciate your helpful suggestions! The suggestions you've provided will be extremely valuable as I continue my research, and I will certainly explore them further.
I had one quick question out of curiosity. I noticed that methods in the "translator" (tangelo/linq/translator) module seem more abundant than those in the "simulator" (tangelo/linq/target) module. For instance, qml has a translator but no simulator. I'm wondering if there are reasons why certain methods don't fit well into the "simulator" category for integration with Tangelo?
Thanks again for your support, and I look forward to diving deeper into Tangelo!

[ValentinS4t1qbit]

Hey @commco789. I can take that last question:

A LOT of simulators exist (wanna write your own? 😉 ), but when we looked at the landscape (usage in the community, scalability, performance, ...) we realized only 3-4 of them were covering 99% of the people and supported our design in Tangelo. In our early days (2019) we supported Rigetti, Microsoft Q# ... but realized those were never used, were not performing better than others we supported, were too cumbersome to support, or their API wasn't working out for our algorithm library. That is: they did not really add much, except that they were yet another dependency to maintain.

With translator and simulator we basically decoupled the representation of a quantum object from simulating it. Instead of supporting a big cartesian product of whatever format or backend you may have designed in your garage, we simply provided conversion functions so that people can easily switch between framework, and run their circuits on the main backends available today.

[commco789]

Thank you very much for your detailed explanation!❤️❤️❤️ I now have a better understanding of why you chose to focus on supporting a few simulators, especially considering the trade-offs in community usage, scalability, and performance. Indeed, supporting too many backends not only increases maintenance costs, but can also become a burden if they don’t offer clear practical advantages.
Your design of decoupling the representation of quantum objects (translator) from their execution (simulator) makes conversions between different frameworks much clearer and more flexible. This is a really clever design and shows thoughtful consideration for both user experience and long-term maintainability. I plan to try writing a simulator myself to gain a deeper understanding of the underlying principles. Thank you again for sharing your insights and inspiration!