This tutorial covers the basic usage of the tequila
-madness
interface that allows to construct basis-set-free qubit Hamiltonians described in arxiv:2008.02819 using the MP2-PNO surrogate described in doi:10.1063/1.5141880.
In order to use madness as chemistry backend in tequila we need a slightly altered version.
We can conveniently install it over the anaconda cloud with:
conda install madtequila -c kottmann
Unfortunately this is currently only supported for Linux-64. Mac users can however resort to manual compilation: follow the instructions from this fork.
Madness is quite sensitive about changes in it's dependencies (e.g. MKL or MPICH), so we recommend to install this in a fresh environment without many other packages. See the last cell of this notebook for an example.
You should see output like
Starting madness calculation with executable: /PATH/TO/bin/pno_integrals
output redirected to he_pno_integrals.out logfile
and the time should be a few seconds (~10 seconds on an intel i5 processor)
import tequila as tq
import time
geometry= "H 0.0 0.0 0.0\nH 0.0 0.0 1.5"
start=time.time()
h2_2_4 = tq.Molecule(geometry=geometry)
print("took {:4.2f}s".format(time.time()-start))
H = h2_2_4.make_hamiltonian()
print("created a qubit Hamiltonian with {} qubits".format(H.n_qubits))
h2_2_4.print_basis_info()
Starting madness calculation with executable: /home/jsk/anaconda3/envs/tq-1.8.1/bin/pno_integrals output redirected to h2_pno_integrals.out logfile finished after 9.780884265899658s took 9.79s created a qubit Hamiltonian with 4 qubits basis_type : custom basis_name : custom orthogonal : True functions : 2 reference : [0] Current Orbitals {idx_total:0, idx:0, occ:2.0, pair:(0, 0)} coefficients: [1. 0.] {idx_total:1, idx:1, occ:0.0425791, pair:(0, 0)} coefficients: [0. 1.]
In the previous cell we created an H2 Molecule in a system adapted orbital basis. Note that we did not specify how many orbitals shall be computed. The default for this is: N-Orbitals = N-electrons. The H2 molecule has 2 electrons, so 2 orbitals (leading to 4 spin-orbitals/qubits) were computed.
Following the notation of arxiv:2008.02819 we denote this with H2/MRA-PNO(2,4).
The printed information shows us, that we have one Hartree-Fock orbital with occupation number 2.0 and one pair-natural orbital (PNO) with occupation number 0.04. These occupation numbers are from the surrogate model (MP2-PNO) which is based on perturbation theory. From the 0.04 occupation number in the PNO we see that the perturbation in the surrogate is not large meaning that it is justified. We can therefore expect, that these two orbitals are close to optimal.
Through the 'mol.make_hamiltonian()' function we get the qubit Hamiltonian.
We can create expectation values for tequila and optimize them in the same manner as with any other Hamiltonian.
The chemistry module offers some convenience in creating pre-defined circuits. One example is the so-called separable pair approximation (SPA) described in arxiv:2105.03836.
H = h2_2_4.make_hamiltonian()
U = h2_2_4.make_ansatz(name="SPA")
E = tq.ExpectationValue(H=H, U=U)
result = tq.minimize(E, silent=True)
print("SPA Energy of H2/({},{}) is {:+2.5f}".format(mol.n_electrons, H.n_qubits, result.energy))
SPA Energy of H2/(2,4) is -1.05358
Other circuits can be initialized in the same manner and the expectatiovalues can be modified and combined like all other tequila expectation values (see the main tutorials or here).
Here are some examples (all perform well since H2 on 4 qubits is quite easy)
U1 = mol.make_ansatz(name="UpCCGSD")
U2 = tq.gates.X(0) + tq.gates.Ry(angle="a",target=2) + tq.gates.CNOT(2,0) + tq.gates.CNOT(0,1) + tq.gates.CNOT(2,3)
for U in [U1,U2]:
E = tq.ExpectationValue(H=H, U=U)
result = tq.minimize(E, silent=True)
V = E**2 - tq.ExpectationValue(H=H**2, U=U)
# make it executable
V = tq.compile(V)
# evaluate with optimized variables
variance = V(result.variables)
print("energy = ", result.energy)
print("variance = ", variance)
energy = -1.0535782052050917 variance = 1.1920928955078125e-07 energy = -1.0535782052050917 variance = 1.1920928955078125e-07
We can of course also compute more orbitals with the keyword n_pno
that defines the total number of PNOs that shall be computed. Let's compute 2 and 3 PNOs (so in total 3 and 4 orbitals leading to 6 and 8 qubits) and compare the VQE energies with an separable pair ansatz (SPA).
h2_2_6 = tq.Molecule(geometry=geometry, n_pno=2)
h2_2_8 = tq.Molecule(geometry=geometry, n_pno=3)
Starting madness calculation with executable: /home/jsk/anaconda3/envs/tq-1.8.1/bin/pno_integrals output redirected to h2_pno_integrals.out logfile finished after 21.300454139709473s Starting madness calculation with executable: /home/jsk/anaconda3/envs/tq-1.8.1/bin/pno_integrals output redirected to h2_pno_integrals.out logfile finished after 21.94425892829895s
for mol in [h2_2_4, h2_2_6, h2_2_8]:
H = mol.make_hamiltonian()
U = mol.make_ansatz(name="SPA")
E = tq.ExpectationValue(H=H, U=U)
result = tq.minimize(E, silent=True)
fci = mol.compute_energy("fci") # needs pyscf installed
print("{:3} qubits, energy = {:2.5f}, fci = {:2.5f}".format(H.n_qubits, result.energy, fci))
4 qubits, energy = -1.05358, fci = -1.05358 6 qubits, energy = -1.05683, fci = -1.05700 8 qubits, energy = -1.05969, fci = -1.05988
We see that the SPA wavefunction produces identical energie than full-CI (FCI) which is equivalent to the exact groundstate in the given orbital basis. This is not further surprising as it was shown that SPA is exact for two-electron systems in arxiv:2207.12421 provided that the orbitals that form the basis are in an optimal linear combination. In this case it looks like the PNOs are already optimal.
Lets do the same computation with standard basis-sets of similar size. As tequila can not exploit PNO-structure to automatically create the SPA ansatz we need to provide the information which orbitals are assigned to which edge of the molecular graph (see arxiv:2207.12421) - for H2 this is trivial as there is just one edge/bond.
# we need pyscf installed for this cell to execute
for basis_set in ["sto-3g", "6-31G"]:
mol = tq.Molecule(geometry=geometry, basis_set=basis_set)
H = mol.make_hamiltonian()
U = mol.make_ansatz(name="SPA", edges=[[i for i in range(mol.n_orbitals)],])
E = tq.ExpectationValue(H=H, U=U)
result = tq.minimize(E, silent=True)
fci = mol.compute_energy("fci")
print("{:3} qubits, energy = {:2.5f}, fci = {:2.5f}".format(H.n_qubits, result.energy, fci))
converged SCF energy = -0.910873554594386 4 qubits, energy = -0.99815, fci = -0.99815 converged SCF energy = -0.997497294328357 8 qubits, energy = -1.04200, fci = -1.05435
When comparing the SPA and FCI energies we see, that the orbitals of 6-31G (canonical Hartree-Fock) are not yet optimal (this was also observed in arxiv:2207.12421). We can optimize them with tq.chemistry.optimize_orbitals
(see paper).
Comparing the FCI energies of the MRA-PNO molecules above and the standard basis sets used here we see that the MRA-PNOs lead to lower energies for the same qubit sizes. Due to the variational principle, this means that we created a better basis with the system adapted approach.
Let's do a similar calculation for the BeH2 molecule. The default is again to compute one spatial orbital for each electron (or in other words: two spatial orbitals for each electron pair). In arxiv:2207.12421 this was called a minimally correlated basis.
Let's compute BeH2 close to it's equilibrium geometry and run an VQE with the SPA ansatz
beh2_geom = "Be 0.0 0.0 0.0\nH 0.0 0.0 {R}\nH 0.0 0.0 -{R}"
beh2_4_8_15 = tq.Molecule(geometry=beh2_geom.format(R=1.5))
H1 = beh2_4_8_15.make_hamiltonian()
print("{:2} electrons".format(beh2_4_8_15.n_electrons))
print("{:2} orbitals".format(beh2_4_8_15.n_orbitals))
print("{:2} qubits".format(H1.n_qubits))
Starting madness calculation with executable: /home/jsk/anaconda3/envs/tq-1.8.1/bin/pno_integrals output redirected to beh2_pno_integrals.out logfile finished after 45.55256175994873s 4 electrons 4 orbitals 8 qubits
A BeH2 molecule with 4 electrons and 8 qubits (spin-orbitals) was computed by default. But shouldn't it be 6 electrons? This is correct, the so called core-electrons (the electron pair 'sitting' close to the Be atom) is automatically frozen (the so-called 'frozen-core' approximation, standard in most quantum chemistry applications). We can deactivate this with tq.Molecule(...,frozen_core=False), but for most cases it won't lead to much (see arxiv:2207.12421 for a detailed analysis using LiH).
So let's compute the SPA-VQE. We will see that it is almost the same as the exact (FCI) energy (see https://arxiv.org/abs/2105.03836).
U = beh2_4_8_15.make_ansatz(name="SPA")
H = beh2_4_8_15.make_hamiltonian()
E = tq.ExpectationValue(H=H, U=U)
result = tq.minimize(E, silent=True)
fci = beh2_4_8_15.compute_energy("fci")
print("BeH2 with R=1.5")
print("SPA error = {:2.5f}".format(fci-result.energy))
BeH2 with R=1.5 SPA error = -0.00170
Let's do this again with a larger bond distance
beh2_4_8_45 = tq.Molecule(geometry=beh2_geom.format(R=4.5))
H = beh2_4_8_45.make_hamiltonian()
U = beh2_4_8_45.make_ansatz(name="SPA")
E = tq.ExpectationValue(H=H, U=U)
result = tq.minimize(E, silent=True)
fci = beh2_4_8_45.compute_energy("fci") # needs pyscf
print("BeH2 with R=1.5")
print("SPA error = {:2.5f}".format(fci-result.energy))
Starting madness calculation with executable: /home/jsk/anaconda3/envs/tq-1.8.1/bin/pno_integrals output redirected to beh2_pno_integrals.out logfile finished after 82.53703570365906s BeH2 with R=1.5 SPA error = -0.19778
The SPA error is much larger (also shown in arxiv:2105.03836) but we can now try to re-optimize the orbitals to be optimal for an SPA ansatz in the given orbital basis. Note that the orbital basis is not changed but that we just form new linear combinations out of the existing MRA-PNO orbitals.
# needs pyscf
opt = tq.chemistry.optimize_orbitals(circuit=U, molecule=beh2_4_8_45, initial_guess="random", silent=True)
print("opt-SPA error = {:2.5f}".format(fci-opt.energy))
opt-SPA error = -0.00029
The datafiles created by the madness backend can be used in order to read in the orbital-data (in form of the corresponding one- and two-body integrals). This is useful to avoid recomputation of the orbitals which is costly in madness or to read in pre-computed orbitals when madness is not installed on the system.
The data input/output is controlled over two keywords:
name
: the name of the molecule. If not set this is automatically deduced from the geometrydatadir
: the directory where the data should be stored or loaded from (needs tq version 1.8.1 or larger)n_pno
: if not set a minimally correlated set of orbitals is computed, if set to n_pno="read"
data is read inThe three files generated by madness are
in "Mulliken" (or "Chemist") notation (see the tequila paper arxiv:2011.03057).
The gitub repo github.com/kottmanj/moldata contains for example a small collection of precomputed orbitals.
Here is a small example
import tequila as tq
mol1 = tq.Molecule(name="my_molecule", geometry="he 0.0 0.0 0.0", datadir="my_data_dir")
# see in the printout where the two integral tensors came from
# print(mol)
print("contents of directory: my_data_dir")
import os
print(os.listdir("my_data_dir"))
# now read it in again
mol2 = tq.Molecule(name="my_molecule", geometry="he 0.0 0.0 0.0", datadir="my_data_dir", n_pno="read")
print("new molecule with ", mol2.n_orbitals, " orbitals")
Starting madness calculation with executable: /home/jsk/anaconda3/envs/tq-1.8.1/bin/pno_integrals output redirected to my_molecule_pno_integrals.out logfile finished after 13.88163685798645s contents of directory: my_data_dir ['my_molecule_pno_integrals.out', 'my_molecule_gtensor.npy', 'my_molecule_pnoinfo.txt', 'my_molecule_htensor.npy'] new molecule with 2 orbitals
If you know how the mandess input structure works you can modify the keywords like this (HF section dft
and pno section).
Take for example a look into the 'he_pno_integrals.out' file to see all possible keywords.
Consider however, that not all keywords will have an effect on the specialized pno_integrals routine.
mol = tq.Molecule(geometry="he 0.0 0.0 0.0", dft={"k":9, "L":25.0}, pno={"maxrank":4}, frozen_core=True)
Starting madness calculation with executable: /home/jsk/anaconda3/envs/tq-1.8.1/bin/pno_integrals output redirected to he_pno_integrals.out logfile finished after 30.72634220123291s
One of tequila
s primary aims is to simplify usage of many specialized algorithms and programs. In the following you will find a list of articles that describe some of the technology that is applied behind the scenes when you are using tequila with madness as chemistry backend:
madness
overviewtequila
overviewOpenFermion
: Fermion to qubit mappingsqulacs
: Qulacs quantum circuit simulatorUsed anaconda3 with conda v4.12.0
conda create -n test_env python=3.8
conda activate test_env
conda install madtequila -c kottmann
python -m pip install --upgrade pip
python -m pip install "tequila-basic==1.8.1"
python -m pip install "qulacs==0.3.1"
python -m pip install "pyscf==2.0.1"