From 1c46ce9dad08b4c4aa36e820ab3d5ff439585618 Mon Sep 17 00:00:00 2001 From: "Christian B. Mendl" Date: Thu, 25 Jan 2024 20:17:01 +0100 Subject: [PATCH] adding 1D Fermi-Hubbard model example (with modified wrap-around hopping term to preserve translational invariance) --- .../fermi_hubbard1d_blockenc_opt.py} | 0 ...bard1d_dynamics_approx_larger_systems.hdf5 | Bin 0 -> 2112 bytes ...ubbard1d_dynamics_approx_larger_systems.py | 132 +++++++++++ .../fermi_hubbard1d_dynamics_circuit.py | 209 ++++++++++++++++++ .../fermi_hubbard1d_dynamics_opt.py | 178 +++++++++++++++ .../fermi_hubbard1d_dynamics_opt_n13.hdf5 | Bin 0 -> 7640 bytes .../fermi_hubbard1d_dynamics_opt_n21.hdf5 | Bin 0 -> 9688 bytes .../fermi_hubbard1d_dynamics_opt_n5.hdf5 | Bin 0 -> 4312 bytes .../fermi_hubbard1d_dynamics_opt_n9.hdf5 | Bin 0 -> 6616 bytes 9 files changed, 519 insertions(+) rename examples/{fermihubbard1d/fermihubbard_blockenc_opt.py => fermi_hubbard1d/fermi_hubbard1d_blockenc_opt.py} (100%) create mode 100644 examples/fermi_hubbard1d/fermi_hubbard1d_dynamics_approx_larger_systems.hdf5 create mode 100644 examples/fermi_hubbard1d/fermi_hubbard1d_dynamics_approx_larger_systems.py create mode 100644 examples/fermi_hubbard1d/fermi_hubbard1d_dynamics_circuit.py create mode 100644 examples/fermi_hubbard1d/fermi_hubbard1d_dynamics_opt.py create mode 100644 examples/fermi_hubbard1d/fermi_hubbard1d_dynamics_opt_n13.hdf5 create mode 100644 examples/fermi_hubbard1d/fermi_hubbard1d_dynamics_opt_n21.hdf5 create mode 100644 examples/fermi_hubbard1d/fermi_hubbard1d_dynamics_opt_n5.hdf5 create mode 100644 examples/fermi_hubbard1d/fermi_hubbard1d_dynamics_opt_n9.hdf5 diff --git a/examples/fermihubbard1d/fermihubbard_blockenc_opt.py b/examples/fermi_hubbard1d/fermi_hubbard1d_blockenc_opt.py similarity index 100% rename from examples/fermihubbard1d/fermihubbard_blockenc_opt.py rename to examples/fermi_hubbard1d/fermi_hubbard1d_blockenc_opt.py diff --git a/examples/fermi_hubbard1d/fermi_hubbard1d_dynamics_approx_larger_systems.hdf5 b/examples/fermi_hubbard1d/fermi_hubbard1d_dynamics_approx_larger_systems.hdf5 new file mode 100644 index 0000000000000000000000000000000000000000..973b56e13274b5edbd2fc1178d81d473b567df4c GIT binary patch literal 2112 zcmeD5aB<`1lHy_j0S*oZ76t(@6Gr@p0tXI=2#gPtPk=HQp>zk7Ucm%mFfy<+faD~g z;sQ|fD_9`{1yGG4L9VV0K$S4_VKh`5g8>VK1$72Yc$6L?A>iTa7y$BW1k}H<^mGE6 z&S1fioLQ6{pITG|%GMYz#cWwCt0USqw5m0GMDPc|-vs&xo3Tkn;^2Pz4Wz2Sgn(3_#}O zIg&@jH){^wh-&s++b v_8(Wwh+w*GXrFxeK%w}|J@!nd_vISScG&YxRKFk^UTD8VJ$AlislGh`+O%Mz literal 0 HcmV?d00001 diff --git a/examples/fermi_hubbard1d/fermi_hubbard1d_dynamics_approx_larger_systems.py b/examples/fermi_hubbard1d/fermi_hubbard1d_dynamics_approx_larger_systems.py new file mode 100644 index 0000000..ea6721e --- /dev/null +++ b/examples/fermi_hubbard1d/fermi_hubbard1d_dynamics_approx_larger_systems.py @@ -0,0 +1,132 @@ +import numpy as np +import scipy +import matplotlib.pyplot as plt +import matplotlib.colors as mc +import h5py +import qib +import rqcopt as oc + + +def compute_circuit_errors(J, u, Llist, t, nlayers): + """ + Compute circuit approximation errors for various system sizes. + """ + expiH = {} + for L in Llist: + # construct Hamiltonian + # using open boundary conditions and manually adding modified wrap-around hopping term + latt = qib.lattice.IntegerLattice((L,), pbc=False) + field = qib.field.Field(qib.field.ParticleType.FERMION, qib.lattice.LayeredLattice(latt, 2)) + H = qib.FermiHubbardHamiltonian(field, float(J), float(u), spin=True).as_matrix().todense() + # construct hopping term between first and second site and then shift circularly + adj_single = np.zeros((L, L), dtype=int) + adj_single[0, 1] = 1 + adj_single[1, 0] = 1 + hop_single = qib.operator.FieldOperatorTerm([qib.operator.IFODesc(field, qib.operator.IFOType.FERMI_CREATE), + qib.operator.IFODesc(field, qib.operator.IFOType.FERMI_ANNIHIL)], + -J * np.kron(np.identity(2), adj_single)) + T_wrap = np.asarray(qib.FieldOperator([hop_single]).as_matrix().todense()) + # circular shift + T_wrap = oc.permute_operation(T_wrap, list(np.roll(range(L), -1)) + list(np.roll(range(L, 2*L), -1))) + H += T_wrap + # reference time evolution operator + expiH[L] = scipy.linalg.expm(-1j*H*t) + + perm_set = {} + for L in Llist: + # site permutations for correct mapping of two-qubit gates + hop_even_sites = None # equivalent to list(range(2*L)) + hop_odd_sites = list(np.roll(range(L), 1)) + list(np.roll(range(L, 2*L), 1)) + intpot_sites = [j + L*s for j in range(L) for s in range(2)] + perm_set[L] = [hop_even_sites, # even hopping + hop_odd_sites, # odd hopping + intpot_sites] # interaction term + + # load optimized unitaries from disk + Vlist = len(nlayers)*[None] + indices = len(nlayers)*[None] + err_opt = len(nlayers)*[None] + for j, n in enumerate(nlayers): + with h5py.File(f"fermi_hubbard1d_dynamics_opt_n{n}.hdf5", "r") as f: + # parameters must agree + assert f.attrs["L"] == Llist[0] + assert f.attrs["J"] == J + assert f.attrs["u"] == u + assert f.attrs["t"] == t + Vlist[j] = f["Vlist"][:] + assert Vlist[j].shape[0] == n + err_opt[j] = f["err_iter"][-1] + assert (n - 1) % 4 == 0 + indices[j] = [0] + ((n - 1) // 4)*[1, 2, 1, 0] + + # approximation error of optimized circuits for larger system sizes + circ_err = np.zeros((len(Llist), len(nlayers))) + for i, L in enumerate(Llist): + for j in range(len(nlayers)): + perms = [perm_set[L][k] for k in indices[j]] + circ_err[i, j] = np.linalg.norm(oc.brickwall_unitary(Vlist[j], 2*L, perms) - expiH[L], ord=2) + + print("error computation consistency check:", np.linalg.norm(err_opt - circ_err[0], np.inf)) + + return circ_err + + +def main(recompute=True): + + # Hamiltonian parameters + J = 1 + u = 4 + + # various system sizes + Llist = [4, 6] + + # time + t = 0.25 + + # number of circuit layers + nlayers = [5, 9, 13, 21] + + if recompute: + circ_err = compute_circuit_errors(J, u, Llist, t, nlayers) + # save errors to disk + with h5py.File("fermi_hubbard1d_dynamics_approx_larger_systems.hdf5", "w") as f: + f.create_dataset("circ_err", data=circ_err) + # store parameters + f.attrs["J"] = float(J) + f.attrs["u"] = float(u) + f.attrs["t"] = float(t) + f.attrs["Llist"] = Llist + f.attrs["nlayers"] = nlayers + else: + # load errors from disk + with h5py.File("fermi_hubbard1d_dynamics_approx_larger_systems.hdf5", "r") as f: + # parameters must agree + assert f.attrs["J"] == J + assert f.attrs["u"] == u + assert f.attrs["t"] == t + assert np.array_equal(f.attrs["Llist"], Llist) + assert np.array_equal(f.attrs["nlayers"], nlayers) + circ_err = f["circ_err"][:] + print(circ_err) + + # define plot colors + clr_base = mc.to_rgb(plt.rcParams['axes.prop_cycle'].by_key()['color'][0]) + clrs = len(Llist)*[None] + for i in range(len(Llist)): + s = i / (len(Llist) - 1) + clrs[i] = ((1 - s)*clr_base[0], (1 - s)*clr_base[1], (1 - s)*clr_base[2]) + + for i, L in enumerate(Llist): + plt.loglog(nlayers, circ_err[i], '.-', color=clrs[i], label=f"L = {L}") + xt = [5, 6, 9, 13, 17] + plt.xticks(xt, [rf"$\mathdefault{{{l}}}$" if l % 2 == 1 else "" for l in xt]) + plt.xlabel("number of layers") + plt.ylabel("error") + plt.legend(loc="upper right") + plt.title(rf"$\mathrm{{approximating }}\ e^{{-i H^{{\mathrm{{FH}}}} t}} \ \mathrm{{for}} \ J = {J}, u = {u}, t = {t}$") + plt.savefig("fermi_hubbard1d_dynamics_approx_larger_systems.pdf") + plt.show() + + +if __name__ == "__main__": + main() diff --git a/examples/fermi_hubbard1d/fermi_hubbard1d_dynamics_circuit.py b/examples/fermi_hubbard1d/fermi_hubbard1d_dynamics_circuit.py new file mode 100644 index 0000000..5c17dcc --- /dev/null +++ b/examples/fermi_hubbard1d/fermi_hubbard1d_dynamics_circuit.py @@ -0,0 +1,209 @@ +import numpy as np +import scipy +import matplotlib.pyplot as plt +import h5py +import qib +import rqcopt as oc + + +def construct_kinetic_term(J: float): + """ + Construct the kinetic hopping term of the Fermi-Hubbard Hamiltonian + on a one-dimensional lattice, based on Jordan-Wigner encoding. + """ + return -J * np.array([ + [0., 0., 0., 0.], + [0., 0., 1., 0.], + [0., 1., 0., 0.], + [0., 0., 0., 0.]]) + + +def construct_interaction_term(u: float): + """ + Construct the local interaction term of the Fermi-Hubbard Hamiltonian + on a one-dimensional lattice, based on Jordan-Wigner encoding. + """ + return u * np.diag([0., 0., 0., 1.]) + + +def trotterized_time_evolution(L: int, hloc, perm_set, method: oc.SplittingMethod, dt: float, nsteps: int): + """ + Compute the numeric ODE flow operator of the quantum time evolution + based on the provided splitting method. + """ + Vlist = [] + perms = [] + for i, c in zip(method.indices, method.coeffs): + Vlist.append(scipy.linalg.expm(-1j*c*dt*hloc[i])) + perms.append(perm_set[i]) + V = oc.brickwall_unitary(Vlist, L, perms) + return np.linalg.matrix_power(V, nsteps) + + +def main(): + + # side length of lattice (spinful sites) + L = 4 + + # Hamiltonian parameters + J = 1 + u = 4 + + # construct Hamiltonian + # using open boundary conditions and manually adding modified wrap-around hopping term + latt = qib.lattice.IntegerLattice((L,), pbc=False) + field = qib.field.Field(qib.field.ParticleType.FERMION, qib.lattice.LayeredLattice(latt, 2)) + H = qib.FermiHubbardHamiltonian(field, float(J), float(u), spin=True).as_matrix().todense() + # construct hopping term between first and second site and then shift circularly + adj_single = np.zeros((L, L), dtype=int) + adj_single[0, 1] = 1 + adj_single[1, 0] = 1 + hop_single = qib.operator.FieldOperatorTerm([qib.operator.IFODesc(field, qib.operator.IFOType.FERMI_CREATE), + qib.operator.IFODesc(field, qib.operator.IFOType.FERMI_ANNIHIL)], + -J * np.kron(np.identity(2), adj_single)) + T_wrap = np.asarray(qib.FieldOperator([hop_single]).as_matrix().todense()) + # circular shift + T_wrap = oc.permute_operation(T_wrap, list(np.roll(range(L), -1)) + list(np.roll(range(L, 2*L), -1))) + H += T_wrap + # visualize spectrum + λ = np.linalg.eigvalsh(H) + plt.plot(λ, '.') + plt.xlabel(r"$j$") + plt.ylabel(r"$\lambda_j$") + plt.title(f"modified Fermi-Hubbard Hamiltonian for J = {J}, u = {u} on a 1D lattice with {L} sites") + plt.show() + + # reference global unitary + t = 0.25 + expiH = scipy.linalg.expm(-1j*H*t) + + # local Hamiltonian terms + tkin = construct_kinetic_term(J) + vpot = construct_interaction_term(u) + hloc = [tkin, tkin, vpot] + hop_even_sites = list(range(2*L)) + hop_odd_sites = list(np.roll(range(L), 1)) + list(np.roll(range(L, 2*L), 1)) + intpot_sites = [j + L*s for j in range(L) for s in range(2)] + print("hop_even_sites:", hop_even_sites) + print("hop_odd_sites: ", hop_odd_sites) + print("intpot_sites: ", intpot_sites) + perm_set = [None, # even hopping + hop_odd_sites, # odd hopping + intpot_sites] # interaction term + + # real-time evolution via Strang splitting for different time steps + nsteps_stra = np.array([2**i for i in range(4)]) + err_stra = np.zeros(len(nsteps_stra)) + stra = oc.SplittingMethod.suzuki(3, 1) + print("stra.coeffs:", stra.coeffs) + for i, nsteps in enumerate(nsteps_stra): + print("nsteps:", nsteps) + dt = t / nsteps + print("dt:", dt) + W = trotterized_time_evolution(2*L, hloc, perm_set, stra, dt, nsteps) + err_stra[i] = np.linalg.norm(W - expiH, ord=2) + print(f"err_stra[{i}]: {err_stra[i]}") + # convergence plot + dt_list = t / nsteps_stra + plt.loglog(dt_list, err_stra, '.-', label="Strang") + plt.loglog(dt_list, 2*np.array(dt_list)**2, '--', label="~Δt^2") + plt.xlabel("Δt") + plt.ylabel("error") + plt.legend() + plt.title(f"real-time evolution up to t = {t} using Strang splitting") + plt.show() + + # real-time evolution via Suzuki's order-4 splitting for different time steps + nsteps_suz4 = np.array(sorted([2**i for i in range(3)])) + err_suz4 = np.zeros(len(nsteps_suz4)) + suz4 = oc.SplittingMethod.suzuki(3, 2) + for i, nsteps in enumerate(nsteps_suz4): + print("nsteps:", nsteps) + dt = t / nsteps + print("dt:", dt) + W = trotterized_time_evolution(2*L, hloc, perm_set, suz4, dt, nsteps) + err_suz4[i] = np.linalg.norm(W - expiH, ord=2) + print(f"err_suz4[{i}]: {err_suz4[i]}") + # convergence plot + dt_list = t / nsteps_suz4 + plt.loglog(dt_list, err_suz4, '.-', label="Suzuki ") + plt.loglog(dt_list, 0.5*np.array(dt_list)**4, '--', label="~Δt^4") + plt.xlabel("Δt") + plt.ylabel("error") + plt.legend() + plt.title(f"real-time evolution up to t = {t} using Suzuki's order-4 splitting") + plt.show() + + # real-time evolution using Yoshida's minimal order-4 method for different time steps + nsteps_yosh = np.array([2**i for i in range(4)]) + err_yosh = np.zeros(len(nsteps_yosh)) + yosh = oc.SplittingMethod.yoshida4(3) + for i, nsteps in enumerate(nsteps_yosh): + print("nsteps:", nsteps) + dt = t / nsteps + print("dt:", dt) + W = trotterized_time_evolution(2*L, hloc, perm_set, yosh, dt, nsteps) + err_yosh[i] = np.linalg.norm(W - expiH, ord=2) + print(f"err_yosh[{i}]: {err_yosh[i]}") + # convergence plot + dt_list = t / nsteps_yosh + plt.loglog(dt_list, err_yosh, '.-', label="Yoshida order 4") + plt.loglog(dt_list, 10*np.array(dt_list)**4, '--', label="~Δt^4") + plt.xlabel("Δt") + plt.ylabel("error") + plt.legend() + plt.title(f"real-time evolution up to t = {t} using Yoshida's method of order-4") + plt.show() + + # real-time evolution using the order-6 method AY 15-6 by Auzinger et al. for different time steps + nsteps_auzi = np.array([2**i for i in range(3)]) + err_auzi = np.zeros(len(nsteps_auzi)) + # recommended method is m = 5; we use the m = 4 method anyway since it requires two substeps less + auzi = oc.SplittingMethod.auzinger15_6() + for i, nsteps in enumerate(nsteps_auzi): + print("nsteps:", nsteps) + dt = t / nsteps + print("dt:", dt) + W = trotterized_time_evolution(2*L, hloc, perm_set, auzi, dt, nsteps) + err_auzi[i] = np.linalg.norm(W - expiH, ord=2) + print(f"err_auzi[{i}]: {err_auzi[i]}") + # convergence plot + dt_list = t / nsteps_auzi + plt.loglog(dt_list, err_auzi, '.-', label="Auzinger AY 15-6") + plt.loglog(dt_list, 6*np.array(dt_list)**6, '--', label="~Δt^6") + plt.xlabel("Δt") + plt.ylabel("error") + plt.legend() + plt.title(f"real-time evolution up to t = {t} using the AY 15-6 method of order 6") + plt.show() + + # optimized circuits + err_iter_opt = {} + for nlayers in [5, 9, 13, 21]: + with h5py.File(f"fermi_hubbard1d_dynamics_opt_n{nlayers}.hdf5", "r") as f: + # parameters must agree + assert f.attrs["L"] == L + assert f.attrs["J"] == J + assert f.attrs["u"] == u + assert f.attrs["t"] == t + err_iter_opt[nlayers] = f["err_iter"][:] + + # compare in terms of number of layers + plt.loglog([5, 9, 13, 21], + [err_iter_opt[n][-1] for n in [5, 9, 13, 21]], 'o-', linewidth=2, label="opt. circuit") + plt.loglog((stra.num_layers-1)*nsteps_stra + 1, err_stra, '.-', label="Strang") + plt.loglog((suz4.num_layers-1)*nsteps_suz4 + 1, err_suz4, 'v-', label="Suzuki order 4") + plt.loglog((yosh.num_layers-1)*nsteps_yosh + 1, err_yosh, '^-', label="Yoshida order 4") + plt.loglog((auzi.num_layers-1)*nsteps_auzi + 1, err_auzi, 'D-', label="Auzinger AY 15-6") + xt = [3, 5, 7, 9, 13, 25, 49] + plt.xticks(xt, [rf"$\mathdefault{{{l}}}$" for l in xt]) + plt.xlabel("number of layers") + plt.ylabel("error") + plt.legend(loc="upper right") + plt.title(rf"$\mathrm{{approximating }}\ e^{{-i H^{{\mathrm{{FH}}}} t}} \ \mathrm{{on}} \ \mathrm{{1D}} \ \mathrm{{lattice}} \ \mathrm{{with}} \ {L} \ \mathrm{{sites}}, J = {J}, u = {u}, t = {t}$") + plt.savefig("fermi_hubbard1d_dynamics_circuit_error.pdf") + plt.show() + + +if __name__ == "__main__": + main() diff --git a/examples/fermi_hubbard1d/fermi_hubbard1d_dynamics_opt.py b/examples/fermi_hubbard1d/fermi_hubbard1d_dynamics_opt.py new file mode 100644 index 0000000..089cf05 --- /dev/null +++ b/examples/fermi_hubbard1d/fermi_hubbard1d_dynamics_opt.py @@ -0,0 +1,178 @@ +import numpy as np +import scipy +import qib +import h5py +import rqcopt as oc +import matplotlib.pyplot as plt + + +def construct_kinetic_term(J: float): + """ + Construct the kinetic hopping term of the Fermi-Hubbard Hamiltonian + on a one-dimensional lattice, based on Jordan-Wigner encoding. + """ + return -J * np.array([ + [0., 0., 0., 0.], + [0., 0., 1., 0.], + [0., 1., 0., 0.], + [0., 0., 0., 0.]]) + + +def construct_interaction_term(u: float): + """ + Construct the local interaction term of the Fermi-Hubbard Hamiltonian + on a one-dimensional lattice, based on Jordan-Wigner encoding. + """ + return u * np.diag([0., 0., 0., 1.]) + + +def fermi_hubbard1d_dynamics_opt(nlayers: int, bootstrap: bool, method_start: oc.SplittingMethod=None, **kwargs): + """ + Optimize the quantum gates in a brickwall layout to approximate + the time evolution governed by the Fermi-Hubbard Hamiltonian. + """ + # side length of lattice (spinful sites) + L = 4 + # Hamiltonian parameters + J = 1 + u = 4 + + # construct Hamiltonian + # using open boundary conditions and manually adding modified wrap-around hopping term + latt = qib.lattice.IntegerLattice((L,), pbc=False) + field = qib.field.Field(qib.field.ParticleType.FERMION, qib.lattice.LayeredLattice(latt, 2)) + H = qib.FermiHubbardHamiltonian(field, float(J), float(u), spin=True).as_matrix().todense() + # construct hopping term between first and second site and then shift circularly + adj_single = np.zeros((L, L), dtype=int) + adj_single[0, 1] = 1 + adj_single[1, 0] = 1 + hop_single = qib.operator.FieldOperatorTerm([qib.operator.IFODesc(field, qib.operator.IFOType.FERMI_CREATE), + qib.operator.IFODesc(field, qib.operator.IFOType.FERMI_ANNIHIL)], + -J * np.kron(np.identity(2), adj_single)) + T_wrap = np.asarray(qib.FieldOperator([hop_single]).as_matrix().todense()) + # circular shift + T_wrap = oc.permute_operation(T_wrap, list(np.roll(range(L), -1)) + list(np.roll(range(L, 2*L), -1))) + H += T_wrap + + # time + t = 0.25 + + # reference global unitary + expiH = scipy.linalg.expm(-1j*H*t) + + # site permutations for correct mapping of two-qubit gates + hop_even_sites = list(range(2*L)) + hop_odd_sites = list(np.roll(range(L), 1)) + list(np.roll(range(L, 2*L), 1)) + intpot_sites = [j + L*s for j in range(L) for s in range(2)] + print("hop_even_sites:", hop_even_sites) + print("hop_odd_sites: ", hop_odd_sites) + print("intpot_sites: ", intpot_sites) + perm_set = [None, # even hopping + hop_odd_sites, # odd hopping + intpot_sites] # interaction term + + # unitaries used as starting point for optimization + if bootstrap: + # load optimized unitaries for nlayers - 4 from disk + with h5py.File(f"fermi_hubbard1d_dynamics_opt_n{nlayers-4}.hdf5", "r") as f: + # parameters must agree + assert f.attrs["L"] == L + assert f.attrs["J"] == J + assert f.attrs["u"] == u + assert f.attrs["t"] == t + Vlist_start = f["Vlist"][:] + assert Vlist_start.shape[0] == nlayers - 4 + # pad identity matrices + id4 = np.identity(4).reshape((1, 4, 4)) + Vlist_start = np.concatenate((Vlist_start, id4, id4, id4, id4), axis=0) + assert Vlist_start.shape[0] == nlayers + assert (nlayers - 1) % 4 == 0 + indices = [0] + ((nlayers - 1) // 4)*[1, 2, 1, 0] + perms = [perm_set[i] for i in indices] + else: + assert method_start.nterms == 3 + assert len(method_start.coeffs) == nlayers + # local Hamiltonian terms + tkin = construct_kinetic_term(J) + vpot = construct_interaction_term(u) + Vlist_start = [] + perms = [] + for i, c in zip(method_start.indices, method_start.coeffs): + Vlist_start.append(scipy.linalg.expm(-1j*c*t*(tkin if i <= 1 else vpot))) + perms.append(perm_set[i]) + # perform optimization + Vlist, f_iter, err_iter = oc.optimize_brickwall_circuit(2*L, expiH, Vlist_start, perms, **kwargs) + + # visualize optimization progress + print(f"err_iter before: {err_iter[0]}") + print(f"err_iter after {len(err_iter)-1} iterations: {err_iter[-1]}") + plt.semilogy(range(len(err_iter)), err_iter) + plt.xlabel("iteration") + plt.ylabel("spectral norm error") + plt.title(f"optimization progress for a quantum circuit with {len(Vlist)} layers") + plt.show() + # rescaled and shifted target function + plt.semilogy(range(len(f_iter)), 1 + np.array(f_iter) / 2**(2*L)) + plt.xlabel("iteration") + plt.ylabel(r"$1 + f(\mathrm{Vlist})/2^{2L}$") + plt.title(f"optimization target function for a quantum circuit with {len(Vlist)} layers") + plt.show() + + # save results to disk + f_iter = np.array(f_iter) + err_iter = np.array(err_iter) + with h5py.File(f"fermi_hubbard1d_dynamics_opt_n{nlayers}.hdf5", "w") as f: + f.create_dataset("Vlist", data=Vlist) + f.create_dataset("f_iter", data=f_iter) + f.create_dataset("err_iter", data=err_iter) + # store parameters + f.attrs["L"] = L + f.attrs["J"] = float(J) + f.attrs["u"] = float(u) + f.attrs["t"] = float(t) + + +def main(): + + # Note: running this code might take several hours + + # 5 layers + # use a single Strang splitting step as starting point for optimization + strang = oc.SplittingMethod.suzuki(3, 1) + print("strang.coeffs:", strang.coeffs) + print("strang.indices:", strang.indices) + fermi_hubbard1d_dynamics_opt(5, False, strang, niter=10) + + # 9 layers + # use two Strang splitting steps as starting point for optimization + strang = oc.SplittingMethod.suzuki(3, 1) + indices_start_n9, coeffs_start_n9 = oc.merge_layers(2*strang.indices, 2*strang.coeffs) + # divide by 2 since we are taking two steps + coeffs_start_n9 = [0.5*c for c in coeffs_start_n9] + print("coeffs_start_n9:", coeffs_start_n9) + print("indices_start_n9:", indices_start_n9) + strang2 = oc.SplittingMethod(3, indices_start_n9, coeffs_start_n9, 2) + fermi_hubbard1d_dynamics_opt(9, False, strang2, niter=16) + + # 13 layers + # use three Strang splitting steps as starting point for optimization + strang = oc.SplittingMethod.suzuki(3, 1) + indices_start_n13, coeffs_start_n13 = oc.merge_layers(3*strang.indices, 3*strang.coeffs) + # divide by 3 since we are taking three steps + coeffs_start_n13 = [c/3 for c in coeffs_start_n13] + print("coeffs_start_n13:", coeffs_start_n13) + print("indices_start_n13:", indices_start_n13) + strang3 = oc.SplittingMethod(3, indices_start_n13, coeffs_start_n13, 2) + fermi_hubbard1d_dynamics_opt(13, False, strang3, niter=30) + + # 21 layers + # use the order-4 Suzuki method as starting point for optimization + suz4 = oc.SplittingMethod.suzuki(3, 2) + print("suz4.coeffs:", suz4.coeffs) + print("suz4.indices:", suz4.indices) + # note: not fully converged yet, can reach even higher precision with more iterations + fermi_hubbard1d_dynamics_opt(21, False, suz4, niter=80) + + +if __name__ == "__main__": + main() diff --git a/examples/fermi_hubbard1d/fermi_hubbard1d_dynamics_opt_n13.hdf5 b/examples/fermi_hubbard1d/fermi_hubbard1d_dynamics_opt_n13.hdf5 new file mode 100644 index 0000000000000000000000000000000000000000..669a60ba56e08a88f442db8c3abee814090cfc00 GIT binary patch literal 7640 zcmeHMc~BEq7~c?(D=25hAcvq7O*x8+D9A%g5EV+nqg09)f&@{NisC^_H7FvEXssuL zP$f>Q7Lkffl?B9BDRu%XAeUCHfS{~{jDU>k%f7c3)R1VH)^XTB@_Wbc`{nz-_q}~P zyW6}xJax3VT7>k{&>&O^UAoEqNbpr6Iy}e^*2k?E*fERxXT`B=MP#cVwah);o&y5g#wHkIp z7YJ5|L=i+tR1|I@22-aFjDZ|v5o&b(q9196I@1;`Hd4~T_P9LZ{y~?laW)Z0OML+S zmew+KQK3uE0uUNB)14TnO^lPU2||D^(YS1x*V(E;HE9GYDk^j!8hD4ouoyWegkHk4 zKK1or5$qQv3+P8*7qjmZ2KB{4xsWnq)cWSVFSD=K70MKgh3d)-2g)EDu8(x)b!#c5 z_OYHuNQy%6;y2+%>s@Hb%ZE93lEJsn2Zi&|@lx3DRz8?>lECL~750`xpgR|fP$4*Q z7kcAyUASw}y~lzcK~=&eoD2Ki^l6Fx2Fe zktJf9|6ue*TF*i#@q>DD^C#Q0o--7PXE&LNtcSRH44b@IB54J?+g2vsHk+=reugR1 zi|g7o7V0UzUanS1kE_14@K!(QiA&A?v;vQSrHKdc?!D{SKMEMPA5xnWqO0`!+kMX| zOL^Q10t?q#ZEFc7HLvkAPSyv3GJT!FYrIa$&2Oj(B!RjR8czqFC)6d>GiFkY2S&T%a&~0(=IJ+SuA?i3TBR=F3%hv4~^~{qJJf@u6*Qpwjv=6Oi6_ z!?G4ZJ-P9|>Ra!Tc`iuLFX)nQf9Qozw7i(;{0JP_`z)K%x`E@Rz~drgj!WYRrPqJG z$oPkoaQ%6|WIpvv5d!@ozaKJL??lEY9SajEu>PCk#*H1m-U-;3r_y}ujrz3Sxw7(P z{}iFp>w~gO0lPc7j_GTddZ)nMa);>iyN4&T+W^tz+`e+&jFUOpd*qlozrD zG}~1L_$#pgw-)mELOvPEt)TgaO+s3aaLj9fdP>(56aTRd^BcWVnRqF1B+j(GaGM^! zAM~!lH|brs86J8N-v5h>&Pieo!pW@auKb5a^GPEcx3Dk;?*BzK-~i+U?$&8E|KUgz zS}(5J%Z@Faex3*dn@?qS`>z8k z>v+}rvCF`z{H3i1A^`|{UYS+<*T1-L%W~C>Q3AkVt+8v+=?MbBEol!H!EpxvNUHK% zYbpQ?{umWJ{j#|LSlsz~mxhG^h;A%t(}LLfaK`_FZ*mVG;0|P9AOiy#7|6i;p8-Am zJ(79*%Y3)wFyDFLdsajEo^=LJg{=4Q%==bEqrl7d(Qjn0>;AI6^qZFSQtbc5;2+|0 BX6FC^ literal 0 HcmV?d00001 diff --git a/examples/fermi_hubbard1d/fermi_hubbard1d_dynamics_opt_n21.hdf5 b/examples/fermi_hubbard1d/fermi_hubbard1d_dynamics_opt_n21.hdf5 new file mode 100644 index 0000000000000000000000000000000000000000..c36bad4ca1829cd78cc5ed081ef384a2655cba5c GIT binary patch literal 9688 zcmeHNdpuQH8{Upvj9gMi>Ee_~RPH2McztDjuZJX(X1Bbw=rT4+_|{85({o z12+4HVwi7eD0wFEpz2hI8bE_0m^i9m!3R;HCTN6+sDNn3Q9SaHZkC}1!cNXQ6sg+f^r9VKB1 z@gnC1)*_hbppz<$<4Xs%?&BoBi+s|}$3ikr3jhE4U{4ZI#r>t{gSt*))Gvbx zb}kW&fyuyC=tQ>5!&NKTdsI=2po*wA(&;mhT;T;TwLB15W-r_HofP`zIKg?wz~ zbL^M8-n8TCf=wg!#>%E^Jb`-GJtgimpnlg9^B6`!G2M8&e%}z2B@b(k8s zf7g#2szR7})b0~^6{pebpPe{xDJhrtAD;j`;c}TQ$M>gyxf< z{atuD_Z4_!;B6+orxvI3zc$vPD{>g}-%%OABfggh3THnYHFCE!;L2{0IvFql@99YO z^Ie=tJHHz5H;&=^OaP82!~zKag;lAD57XJjm~1X3w0Fb(Q`c_ZuU%_lHkr6Y5eb{7cE@3XZS@%f+U>O($Sy>$6v)P9fn zM*6?C8V?*K=v(ue9q*duYSLr@R_CM)+Gr zFf9qbno*(kQ_O>f&VMVZH?Jae#~P?NFX4Sm_I-H&^XRaK!>>g`^SdiIBEG!QTC*S@ ztzNo(C9hD}NyNw8erS zogb>+Ul~=28$i9BX2JSTzylvAH`JKZ<%6!jR6jlSJX}v|h3lJ_ODk6iT~CzWe14=-64ZM#y%DW{ zpVApBGuFNnTEDU9iETW%UyFyWc(lcH5$@f2Lcbu^8BBYd`8G43e*aHPn{V#A%o$JT zsHGi){wkMj_Y-{Fyy(>QELWlP&mr|{{it?G3+m>wkRRkB8bn2!H)gOYTQ} zO6&RiAfM3jmR+eveEi#ojt&Q3w^$_ZMO^=xF`W?wkMJyVK5}p4&t?uo{%aTagyy{B zf$oM7nFrV0fw$gSt?N(C!G?;F$MRHiY3F0s`X;+4XUy@hes=03{Q0ByyApiZi;JsH zY+WLB{b!PTZwKFuo(1*Fr@TS;|6Ym-EvY;e+W6?!|H|;N#T%EK4ns`;kfb= zLhqd}xcuS%GN4cddZluNUjLDiqAhIr{C8wYuN=$F6vx_^%jO#}fP$&6c&<+-ZTt=r zo^v*7G4SW}3Xs1~j1S_&Y91NdPco$)PnRD`uTpZ@L^-Vw+i0^80M%{4mX|Fh?)&QUXQ}5k+!P#AlQ}eNRsh?`k zFyue%NuS(g_Pn%}8%j)!q4f#Hk|aHsf$jYbmPgzdQ+M_75M|1Dut_-@J5 zcZ3}=W$$|eGqWEdJgG3Jhp=VFWFkD5cn@JXJHD4_KObL=u*Bmd!u*rAABgsMQ85Vp zBJLv$G8@xJv^UypLRjQ>9$|-4Jwl^y8@R8-?|I;Z&_(=PgfRoz2>TD?2wk_z_Y;1% zEHn{z__K!KG}8^p{^s?*A^1Y~6taJ=FCBsoUpFIr>_8_%r{}Wp1~Y_9;rZgmI3?jzj@vrC9q(ji<5HG?IQEw)s82yUsAM|MfIBEtD;$_TGoO&o%5q3Xz9 zV=x_|m&&nA6}qdyjod=s!7A55yc3nZWTWarpjK^qXWUE|pgdmxan5MgAYKk?Gb#vR z0iI*P46NA+xF+%X8}mQmqQtWw6TWr;`Et=ma#I|D^@}E#5=(n}UtB-u-C0w89pQ825H$prkEW_Azqm#4;POC2oO!TaXRl ziI&dVbIKZU3cG#B&$b5jwX4dm9<%~WUmG8vplk&sH*bIDTQLs=7X)FO<>mqBa8A@| zsd>OjB_+0bhb0(Y>Eg(Cvjna!X^k_^S^$-6`i5JkS^(F_M$Abf79f}t!^j-<35yO( zdY9*zgYPb?hO*Ag9mEL#c4L{9b3w>U;{mHWvq99Xy73{^#|-Rp9%a|BV+K@bHTttG z<_u!7rQ*GsLQ~K$s@3vv<|mA1xhWjt&jz4gU9TZ!HaNVhr}QG#2NrD`{ZV$E9_Y+E z;K8`9EA+lfcOCzy|DPH7E6#ud`F|wpZ8-J6B^l~}9`HTuX!xGhfeill-oI1tTLo^f a!EeWh{i*Tt7QY=&{HKLD{@C}|fd2q7o;7L! literal 0 HcmV?d00001 diff --git a/examples/fermi_hubbard1d/fermi_hubbard1d_dynamics_opt_n5.hdf5 b/examples/fermi_hubbard1d/fermi_hubbard1d_dynamics_opt_n5.hdf5 new file mode 100644 index 0000000000000000000000000000000000000000..92f15fbc7a9603d98e7f97d20b72f11941de3177 GIT binary patch literal 4312 zcmeD5aB<`1lHy_j0S*oZ76t(@6Gr@pf*Vj_7#}8|0A(;j=?*Bpf(gQ4WMF3i$w@%P z1)%27V220N=S0C_UUmz{AxMoKzy99@T)-C!pyJ z=Krvq%;FLtCoMj+B(;cvA+@Lo%7UnXrBef#1_CrA11mWFLI_X_wE*j2WQLf;3e^uy zp1}MA0l37O85+Pwi2>aWO3Oee5(NYm>=`(~DjgV1xEW0TGczz`LUlsKaT|k8FDuYA zkPn!cn7|(4faXx3w_wJ>ynteGu%EvRBLg#Bkpj>i_y8>5cd+481qlm~iVCc@;@3MA zWiitgB%~M_A?XTIATuz)Xl7vifMl@6CCLAH;*%E=KH%&FVS>pHB=*2F355qY)QK>f z{O|~Xh6k)9f%$!Kg$JyjgjL5kpw8rBfRtQ}3``6>z+C796&C>JT1eTW02YCyBA6Id z5l-rZz2O1z{hgby?7^voNQpmMH!^v5y|*_#6eO`-Zwb(|C{%0v13ks7SnIT<33SHrRU~SQ%@+2d3T`YF^L;)4&|5#Z>eEji<%$ z_HRF7Q_!}1^S%G??S1#~J%1)rPB4B5#UIQ(mqTFl_7hT1NFHXM0hP>~lyoiQUBN^9 zW!3>v?pIzx-S@$ML792V^HK8292){8(t#g|0LsO%@<4R;pxQ{8vo#M z>R3;;_`k*fqT^-W33&RSA+<7Juk zOn)blE_?oc_8F7?_7leo9*|Vzqi%_NNeu8N)}BZ*hZyCLx8Lh z5QH``V67imznurx&jGeGM1buKFR1Cb+MlrY21MBhdvyJ!BVKh3Yz(>?l3lL%_q;5u8*apdQtL(kGzl z4Ceo^oXp}9ASW$8vm~{Mfg!c12+D$}fTdFdm<9qgBLgcq{Xz&(3bg?1U}T1v#0k|8 zO`gE~0|B_inHd_uMu`F44NA*EClUn&73>)}z$zUWOt=|L{xdT$WI}a9#Bm#gO)o3Z zG>{LNn3%vG;eh5)ptoSg!MuQCaIl}h3nK$FT#*9M9ryq&-*>R#Qw0eNkctXc{OX2_ zEM~fbgcKtqBwfMEIT+0hj31B;wzvd^37+^AgoIfGgo2YDaMoa85DE{dCm3Kf`Qgz3 zH3HrL==^~V6<9sV0Ihd#K%K+G04cc`8JHM&fVt2IDlP!bwUDw$0W1PZMKCd_BAnC* z`?d{cM&%D@?8T{sNQo_prkr=azPDe!mA}R=bEYlNmX*5~-fytAj*eThqI{)2A!GjZ zGE4MDx9)|h-wZYH&pQv-^4XJZ3F#+J9_U{&PmT4DVgCJK|KRPbRjwv&6r2BLo!-M! zU*6l_uGE}9mnj=;-oNK9dv8^y2HiWlilG0`nz3grVra3w4ORacY97d64a}t!yAKjR zAMDkp9AZ^5YO|F(-}LqW|M&JkgK`qiztG%IT=?NuB?L83PslS9XdXfT;x-8)BL+3^ z)}+=HYM3W;SjjF?~xy56@8BniIsue}06S`XR}by)gC6Q1d{(Nl2}y*nN;t`CxzLT>UkJgK@S3 z(ERmi>CY@`ew7hIuMQko>=%PY5OdJ4|YHn!yrJark!a zhm`-PyUxVV{MiXMFDSR!R%LY{QRRP^^LYVLhqk>+Q1!h~^FZ;|@Ufd>_rcZwssNw+^X>y%6_PN`A06$m%gYB$#Ep&i+@+rvLc|7<_xy-3HQ5o6;w~0Mb*o zb58pPq{WyHF9Xpo2DG7bFWg&hcz>?FdGf?B&(=(^-*Y`cdG+Exd&z6%Yl{Ci*)LyI z)Y%xp>M2~~Z?N$eh_Vm%==w`Z(l6kEt{*g}4Fbg3f59FAHlj(t literal 0 HcmV?d00001