From 202d72a9bebf3f03545b30d84d140b42fbbe626d Mon Sep 17 00:00:00 2001 From: liwt31 Date: Mon, 29 Apr 2024 16:13:28 +0800 Subject: [PATCH] add ttns examples --- example/run.sh | 11 +- example/ttns/junction_ft.py | 237 +++++++++++++++++++++++++++++ example/ttns/junction_zt.py | 180 ++++++++++++++++++++++ example/ttns/sbm_ft.py | 123 +++++++++++++++ example/ttns/sbm_zt.py | 107 +++++++++++++ renormalizer/mps/mps.py | 139 ++++++++++------- renormalizer/mps/tests/test_mpo.py | 1 + renormalizer/sbm/__init__.py | 2 +- renormalizer/sbm/lib.py | 78 +++++++++- renormalizer/sbm/sbm.py | 5 +- renormalizer/tn/hop_expr.py | 8 +- renormalizer/tn/node.py | 15 ++ renormalizer/tn/tree.py | 74 ++++----- renormalizer/utils/quantity.py | 3 + setup.py | 3 +- 15 files changed, 882 insertions(+), 104 deletions(-) create mode 100644 example/ttns/junction_ft.py create mode 100644 example/ttns/junction_zt.py create mode 100644 example/ttns/sbm_ft.py create mode 100644 example/ttns/sbm_zt.py diff --git a/example/run.sh b/example/run.sh index e2b72f6a..86eefd4f 100644 --- a/example/run.sh +++ b/example/run.sh @@ -1,6 +1,15 @@ export PYTHONPATH=../:PYTHONPATH code=0 -for python_args in fmo.py sbm.py h2o_qc.py "dynamics.py std.yaml" "transport_kubo.py std.yaml"; do +for python_args in fmo.py \ + sbm.py \ + h2o_qc.py \ + "dynamics.py std.yaml" \ + "transport_kubo.py std.yaml" \ + ./ttns/junction_zt.py \ + "./ttns/junction_ft.py 32 1 100" \ + "./ttns/sbm_zt.py 050 001 050"\ + ./ttns/sbm_ft.py +do echo ============================$python_args============================= timeout 20s python $python_args exit_code=$? diff --git a/example/ttns/junction_ft.py b/example/ttns/junction_ft.py new file mode 100644 index 00000000..532fa18d --- /dev/null +++ b/example/ttns/junction_ft.py @@ -0,0 +1,237 @@ +import sys + +from renormalizer import BasisHalfSpin, Op, Quantity, BasisSHO, BasisDummy +from renormalizer.mps.mps import expand_bond_dimension_general +from renormalizer.sbm import ColeDavidsonSDF +from renormalizer.utils.configs import ( + EvolveMethod, + EvolveConfig, + CompressConfig, + CompressCriteria, +) +from renormalizer.utils import constant, log +from renormalizer.utils.log import package_logger as logger +from renormalizer.tn import BasisTree, TreeNodeBasis, TTNS, TTNO + +import numpy as np + + +Ms_str = sys.argv[1] # 8 16 24 32 +initial_str = sys.argv[2] # 1 0 +temperature_str = sys.argv[3] # 000 100 200 300 + +dump_dir = "./" +job_name = f"Ms{Ms_str}_i{initial_str}_temperature{temperature_str}" #################### +log.register_file_output(dump_dir + job_name + ".log", mode="w") + +Ms = int(Ms_str) + +n_ph_mode = 1000 +omega_c = Quantity(500, "cm-1").as_au() +ita = Quantity(2000, "cm-1").as_au() / 2 +beta = 0.5 +upper_limit = Quantity(1, "eV").as_au() * 10 +logger.info(("phonon parameters", omega_c, ita, beta, upper_limit)) +sdf = ColeDavidsonSDF(ita, omega_c, beta, upper_limit) +w, c2 = sdf.Wang1(n_ph_mode) +c = np.sqrt(c2) +logger.info(w) +logger.info(c) + +reno = sdf.reno(w[-1]) +logger.info(f"renormalization constant: {reno}") + +temperature = Quantity(int(temperature_str), "K").to_beta() +logger.info(f"temperature beta: {temperature}") + +n_e_mode = 320 + +beta_e = Quantity(1, "eV").as_au() * reno +alpha_e = Quantity(0.2, "eV").as_au() * reno +v = 0.1 * reno +mu_l = Quantity(v / 2, "eV").as_au() +mu_r = Quantity(-v / 2, "eV").as_au() + + +e_k = np.arange(1, n_e_mode + 1) / (n_e_mode + 1) * 4 * beta_e - 2 * beta_e +rho_e = 1 / (e_k[1] - e_k[0]) +e_k_l = e_k - mu_l +e_k_r = e_k - mu_r + +mode_with_e = [(f"L{i}", e) for i, e in enumerate(e_k_l)] + [(f"R{i}", e) for i, e in enumerate(e_k_r)] +mode_with_e.sort(key=lambda x: x[1]) +logger.info(mode_with_e) + +basis = [] +first_positive = True +for mode, e in mode_with_e: + if e > 0 and first_positive: + first_positive = False + basis.append(BasisHalfSpin("s")) + basis.append(BasisHalfSpin((mode, "p"))) + basis.append(BasisHalfSpin((mode, "q"))) + +dofs = [b.dofs[0] for b in basis] +logger.info(dofs) +s_idx = dofs.index("s") +basis_tree_l = BasisTree.binary_mctdh(basis[:s_idx], dummy_label="EL-dummy") +basis_tree_r = BasisTree.binary_mctdh(basis[s_idx + 1 :], dummy_label="ER-dummy") + + +ham_terms = [] +i_l_terms = [] +i_r_terms = [] +for mode, e in mode_with_e: + if mode[0] == "L": + mu = mu_l + i_terms = i_l_terms + else: + assert mode[0] == "R" + mu = mu_r + i_terms = i_r_terms + + ham_terms.append(Op("+ -", (mode, "p"), e + mu)) + ham_terms.append(Op("+ -", (mode, "q"), -(e + mu))) + v2 = alpha_e**2 / beta_e**2 * np.sqrt(4 * beta_e**2 - (e + mu) ** 2) / 2 / np.pi / rho_e + v = np.sqrt(v2) + theta = np.arctan(np.exp(-temperature * e / 2)) + idx = dofs.index((mode, "p")) + if idx < s_idx: + z_idx = list(range(idx + 1, s_idx)) + else: + assert s_idx < idx + z_idx = list(range(s_idx + 1, idx)) + z_dofs = [dofs[i] for i in z_idx] + op1 = Op( + "+ " + "Z " * len(z_idx) + "-", + [(mode, "p")] + z_dofs + ["s"], + v * np.cos(theta), + ) + op2 = Op( + "- " + "Z " * len(z_idx) + "+", + [(mode, "p")] + z_dofs + ["s"], + v * np.cos(theta), + ) + idx = dofs.index((mode, "q")) + if idx < s_idx: + z_idx = list(range(idx + 1, s_idx)) + else: + assert s_idx < idx + z_idx = list(range(s_idx + 1, idx)) + z_dofs = [dofs[i] for i in z_idx] + op3 = Op( + "- " + "Z " * len(z_idx) + "-", + [(mode, "q")] + z_dofs + ["s"], + v * np.sin(theta), + ) + op4 = Op( + "+ " + "Z " * len(z_idx) + "+", + [(mode, "q")] + z_dofs + ["s"], + v * np.sin(theta), + ) + ham_terms.extend([op1, op2, op3, op4]) + # move 1j to expectation + i_terms.extend(op2 - op1 + op4 - op3) + +initial_occupied = initial_str == "1" + +if initial_occupied: + ham_terms.append(Op("+ -", "s", qn=[0, 0], factor=-4 * (c**2 / w**2).sum())) + +# vibrations at last + +# boson energy +for imode in range(n_ph_mode): + op1 = Op(r"p^2", f"v_{imode}_p", factor=0.5, qn=0) + op2 = Op(r"x^2", f"v_{imode}_p", factor=0.5 * w[imode] ** 2, qn=0) + ham_terms.extend([op1, op2]) + op1 = Op(r"p^2", f"v_{imode}_q", factor=-0.5, qn=0) + op2 = Op(r"x^2", f"v_{imode}_q", factor=-0.5 * w[imode] ** 2, qn=0) + ham_terms.extend([op1, op2]) + +theta_array = np.arctanh(np.exp(-w * temperature / 2)) +# system-boson coupling +for imode in range(n_ph_mode): + sys_op = Op("+ -", "s", qn=[0, 0]) + if initial_occupied: + sys_op = sys_op - Op.identity("s") + + theta = theta_array[imode] + op1 = sys_op * Op(r"x", f"v_{imode}_p", factor=2 * c[imode] * np.cosh(theta), qn=[0]) + op2 = sys_op * Op(r"x", f"v_{imode}_q", factor=2 * c[imode] * np.sinh(theta), qn=[0]) + ham_terms.extend(op1 + op2) + +nbas = np.max([16 * c2 / w**3 * np.cosh(theta_array) ** 2, np.ones(n_ph_mode) * 4], axis=0) +nbas = np.round(nbas).astype(int) +nbas = np.min([nbas, np.ones(n_ph_mode) * 512], axis=0) +logger.info(nbas) +basis_list_phonon = [] +for imode in range(n_ph_mode): + basis_list_phonon.append(BasisSHO(f"v_{imode}_p", w[imode], int(nbas[imode]))) + basis_list_phonon.append(BasisSHO(f"v_{imode}_q", w[imode], int(nbas[imode]))) + +labels = np.array([[nbas > Ms], [nbas > Ms]]).T.ravel() +basis_tree_phonon = BasisTree.binary_mctdh( + basis_list_phonon, + contract_primitive=True, + contract_label=labels, + dummy_label="phonon-dummy", +) +node1 = TreeNodeBasis([BasisDummy("dummy")]) +node1.add_child([basis_tree_l.root, basis_tree_r.root]) +node2 = TreeNodeBasis([basis[s_idx]]) +node2.add_child([node1, basis_tree_phonon.root]) +basis_tree = BasisTree(node2) +basis_tree.print(logger.info) + + +# model = Model(basis, ham_terms) +ttno = TTNO(basis_tree, ham_terms) +i_l_mpo = TTNO(basis_tree, i_l_terms) +i_r_mpo = TTNO(basis_tree, i_r_terms) +# n_l_mpo = TTNO(basis_tree, terms=[Op("+ -", f"L{i}") for i in range(n_e_mode)]) +# n_r_mpo = TTNO(basis_tree, terms=[Op("+ -", f"R{i}") for i in range(n_e_mode)]) +n_s_mpo = TTNO(basis_tree, terms=Op("+ -", "s")) +ttno.print_shape(False, logger.info) +i_l_mpo.print_shape(False, logger.info) +i_r_mpo.print_shape(False, logger.info) +# n_r_mpo.print_shape(False, logger.info) +# n_s_mpo.print_shape(False, logger.info) +# 0 - [1, 0] (spin up) means occupation, 1 - [0, 1] (spin down) means unoccupation +# initial condition is taken care of in the Hamiltonian +condition = {dofs[i]: 1 for i in range(len(dofs))} +if initial_occupied: + condition["s"] = 0 +else: + condition["s"] = 1 + +ttns = TTNS(basis_tree, condition=condition) +ttns.compress_config = CompressConfig(CompressCriteria.fixed, max_bonddim=Ms) +ttns = expand_bond_dimension_general(ttns, ttno, ex_mps=None) +ttns.evolve_config = EvolveConfig(EvolveMethod.tdvp_ps) +ttns.print_shape(print_function=logger.info, full=False) + +step = 0.5 * constant.fs2au +# step = 5 +nsteps = 200 +au2muA = 6.623618237510e3 +i = 0 +current_list = [] +while True: + i_l = (1j * ttns.expectation(i_l_mpo)).real + i_r = (1j * ttns.expectation(i_r_mpo)).real + # n_l = mps.expectation(n_l_mpo) + # n_r = mps.expectation(n_r_mpo) + n_s = ttns.expectation(n_s_mpo) + logger.info((i, ttns.bond_dims)) + current = (i_r - i_l) / 2 * au2muA + # logger.info((n_l, n_r, n_s, i_l*au2muA, i_r*au2muA, current)) + logger.info((n_s, i_l * au2muA, i_r * au2muA, current)) + current_list.append(current) + i += 1 + if i == nsteps: + break + if i > 0: + ttns.evolve_config = EvolveConfig(EvolveMethod.tdvp_ps) + ttns = ttns.evolve(ttno, step) +logger.info(current_list) diff --git a/example/ttns/junction_zt.py b/example/ttns/junction_zt.py new file mode 100644 index 00000000..7f3a283a --- /dev/null +++ b/example/ttns/junction_zt.py @@ -0,0 +1,180 @@ +from renormalizer import BasisHalfSpin, Op, Quantity, BasisSHO, BasisDummy +from renormalizer.mps.mps import expand_bond_dimension_general +from renormalizer.sbm import ColeDavidsonSDF +from renormalizer.utils.configs import ( + EvolveMethod, + EvolveConfig, + CompressConfig, + CompressCriteria, +) +from renormalizer.utils import constant +from renormalizer.utils.log import package_logger as logger +from renormalizer.tn import BasisTree, TreeNodeBasis, TTNS, TTNO + +import numpy as np + + +n_ph_mode = 500 +omega_c = Quantity(500, "cm-1").as_au() +ita = Quantity(2000, "cm-1").as_au() / 2 +beta = 0.5 +upper_limit = Quantity(1, "eV").as_au() * 5 +logger.info(("phonon parameters", omega_c, ita, beta, upper_limit)) +sdf = ColeDavidsonSDF(ita, omega_c, beta, upper_limit) +w, c2 = sdf.Wang1(n_ph_mode) +c = np.sqrt(c2) +logger.info(w) +logger.info(c) + +reno = sdf.reno(w[-1]) +logger.info(f"renormalization constant: {reno}") + + +n_e_mode = 160 + +beta_e = Quantity(1, "eV").as_au() * reno +alpha_e = Quantity(0.2, "eV").as_au() * reno +v = 0.1 * reno +mu_l = Quantity(v / 2, "eV").as_au() +mu_r = Quantity(-v / 2, "eV").as_au() + + +e_k = np.arange(1, n_e_mode + 1) / (n_e_mode + 1) * 4 * beta_e - 2 * beta_e +rho_e = 1 / (e_k[1] - e_k[0]) +e_k_l = e_k - mu_l +e_k_r = e_k - mu_r + +mode_with_e = [(f"L{i}", e) for i, e in enumerate(e_k_l)] + [(f"R{i}", e) for i, e in enumerate(e_k_r)] +mode_with_e.sort(key=lambda x: x[1]) +logger.info(mode_with_e) + +basis = [] +first_positive = True +for mode, e in mode_with_e: + if e > 0 and first_positive: + first_positive = False + basis.append(BasisHalfSpin("s")) + basis.append(BasisHalfSpin(mode)) + +dofs = [b.dofs[0] for b in basis] +logger.info(dofs) +s_idx = dofs.index("s") +basis_tree_l = BasisTree.binary_mctdh(basis[:s_idx], dummy_label="EL-dummy") +basis_tree_r = BasisTree.binary_mctdh(basis[s_idx + 1 :], dummy_label="ER-dummy") + + +ham_terms = [] +i_l_terms = [] +i_r_terms = [] +for mode, e in mode_with_e: + if mode[0] == "L": + mu = mu_l + i_terms = i_l_terms + else: + assert mode[0] == "R" + mu = mu_r + i_terms = i_r_terms + + ham_terms.append(Op("+ -", mode, e + mu)) + v2 = alpha_e**2 / beta_e**2 * np.sqrt(4 * beta_e**2 - (e + mu) ** 2) / 2 / np.pi / rho_e + v = np.sqrt(v2) + idx = dofs.index(mode) + if idx < s_idx: + z_idx = list(range(idx + 1, s_idx)) + else: + assert s_idx < idx + z_idx = list(range(s_idx + 1, idx)) + z_dofs = [dofs[i] for i in z_idx] + op1 = Op("+ " + "Z " * len(z_idx) + "-", [mode] + z_dofs + ["s"], v) + op2 = Op("- " + "Z " * len(z_idx) + "+", [mode] + z_dofs + ["s"], v) + ham_terms.extend([op1, op2]) + # move 1j to expectation + i_terms.extend(op2 - op1) + +initial_occupied = True + +if initial_occupied: + ham_terms.append(Op("+ -", "s", qn=[0, 0], factor=-4 * (c**2 / w**2).sum())) + +# vibrations at last + +# boson energy +for imode in range(n_ph_mode): + op1 = Op(r"p^2", f"v_{imode}", factor=0.5, qn=0) + op2 = Op(r"x^2", f"v_{imode}", factor=0.5 * w[imode] ** 2, qn=0) + ham_terms.extend([op1, op2]) + +# system-boson coupling +for imode in range(n_ph_mode): + sys_op = Op("+ -", "s", qn=[0, 0]) + if initial_occupied: + sys_op = sys_op - Op.identity("s") + + op = sys_op * Op(r"x", f"v_{imode}", factor=2 * c[imode], qn=[0]) + ham_terms.extend(op) + +nbas = np.max([16 * c2 / w**3, np.ones(n_ph_mode) * 4], axis=0) +nbas = np.round(nbas).astype(int) +logger.info(nbas) +basis_list_phonon = [] +for imode in range(n_ph_mode): + basis_list_phonon.append(BasisSHO(f"v_{imode}", w[imode], int(nbas[imode]))) + +basis_tree_phonon = BasisTree.binary_mctdh(basis_list_phonon, dummy_label="phonon-dummy") +node1 = TreeNodeBasis([BasisDummy("dummy")]) +node1.add_child([basis_tree_l.root, basis_tree_r.root]) +node2 = TreeNodeBasis([basis[s_idx]]) +node2.add_child([node1, basis_tree_phonon.root]) +basis_tree = BasisTree(node2) +basis_tree.print(logger.info) + + +# model = Model(basis, ham_terms) +ttno = TTNO(basis_tree, ham_terms) +i_l_mpo = TTNO(basis_tree, i_l_terms) +i_r_mpo = TTNO(basis_tree, i_r_terms) +n_l_mpo = TTNO(basis_tree, terms=[Op("+ -", f"L{i}") for i in range(n_e_mode)]) +n_r_mpo = TTNO(basis_tree, terms=[Op("+ -", f"R{i}") for i in range(n_e_mode)]) +n_s_mpo = TTNO(basis_tree, terms=Op("+ -", "s")) +ttno.print_shape(False, logger.info) +i_l_mpo.print_shape(False, logger.info) +i_r_mpo.print_shape(False, logger.info) +n_r_mpo.print_shape(False, logger.info) +n_s_mpo.print_shape(False, logger.info) +# 0 - [1, 0] (spin up) means occupation, 1 - [0, 1] (spin down) means unoccupation +# initial condition is taken care of in the Hamiltonian +condition = {dofs[i]: 1 for i in range(s_idx + 1, len(dofs))} +if initial_occupied: + condition["s"] = 0 +else: + condition["s"] = 1 + +ttns = TTNS(basis_tree, condition=condition) +ttns.compress_config = CompressConfig(CompressCriteria.fixed, max_bonddim=32) +ttns = expand_bond_dimension_general(ttns, ttno, ex_mps=None) +ttns.evolve_config = EvolveConfig(EvolveMethod.tdvp_ps) +ttns.print_shape(print_function=logger.info, full=False) + +step = 0.5 * constant.fs2au +# step = 5 +nsteps = 100 +au2muA = 6.623618237510e3 +i = 0 +current_list = [] +while True: + i_l = (1j * ttns.expectation(i_l_mpo)).real + i_r = (1j * ttns.expectation(i_r_mpo)).real + n_l = ttns.expectation(n_l_mpo) + n_r = ttns.expectation(n_r_mpo) + n_s = ttns.expectation(n_s_mpo) + logger.info((i, ttns.bond_dims)) + current = (i_r - i_l) / 2 * au2muA + logger.info((n_l, n_r, n_s, i_l * au2muA, i_r * au2muA, current)) + current_list.append(current) + i += 1 + if i == nsteps: + break + if i > 0: + ttns.evolve_config = EvolveConfig(EvolveMethod.tdvp_ps) + ttns = ttns.evolve(ttno, step) +logger.info(current_list) diff --git a/example/ttns/sbm_ft.py b/example/ttns/sbm_ft.py new file mode 100644 index 00000000..a30a1784 --- /dev/null +++ b/example/ttns/sbm_ft.py @@ -0,0 +1,123 @@ +import sys + +from renormalizer.model import Op +from renormalizer.model import basis as ba +from renormalizer.mps.mps import expand_bond_dimension_general +from renormalizer.sbm import ColeDavidsonSDF +from renormalizer.utils import EvolveConfig, CompressConfig, CompressCriteria, EvolveMethod +from renormalizer.utils import log +from renormalizer.tn import BasisTree, TTNO, TTNS, TreeNodeBasis + +import numpy as np + + +# ita_str = sys.argv[1] # 010, 025, 100 +# omega_c_str = sys.argv[2] # 001, 100 +# beta_str = sys.argv[3] # 025, 100 +# temperature_str = sys.argv[4] # 05, 10, 15 +ita_str = "010" +omega_c_str = "010" +beta_str = "025" +temperature_str = "20" + + +ita = int(ita_str) / 10 # 1, 2.5, 5, 10 +eps = 0 +Delta = 1 +omega_c = int(omega_c_str) / 10 # 0.1, 1, 10 + +beta = int(beta_str) / 100 # 0.25, 0.5, 0.75, 1 +temperature = int(temperature_str) / 10 # 0.5, 1, 1.5 + +from renormalizer.utils.log import package_logger +logger = package_logger +dump_dir = "./" +job_name = f"ps1_binary_ita{ita_str}_omega{omega_c_str}_beta{beta_str}_temperature{temperature_str}" #################### +log.register_file_output(dump_dir+job_name+".log", mode="w") + + +nmodes = 1000 +Ms = 20 +upper_limit = 30 +sdf = ColeDavidsonSDF(ita, omega_c, beta, upper_limit) + +w, c2 = sdf.Wang1(nmodes) +c = np.sqrt(c2) +logger.info(f"w:{w}") +logger.info(f"c:{c}") + +reno = sdf.reno(w[-1]) +logger.info(f"renormalization constant: {reno}") +Delta *= reno + +ham_terms = [] + +# h_s +ham_terms.extend([Op("sigma_z","spin",factor=eps, qn=0), + Op("sigma_x","spin",factor=Delta, qn=0)]) + + +# boson energy +for imode in range(nmodes): + op1 = Op(r"p^2",f"v_{imode}_p",factor=0.5, qn=0) + op2 = Op(r"x^2",f"v_{imode}_p",factor=0.5*w[imode]**2, qn=0) + ham_terms.extend([op1,op2]) + op1 = Op(r"p^2",f"v_{imode}_q",factor=-0.5, qn=0) + op2 = Op(r"x^2",f"v_{imode}_q",factor=-0.5*w[imode]**2, qn=0) + ham_terms.extend([op1,op2]) + +theta_array = np.arctanh(np.exp(-w/ temperature / 2)) +# system-boson coupling +for imode in range(nmodes): + theta = theta_array[imode] + op = Op(r"sigma_z x", ["spin", f"v_{imode}_p"], factor=np.cosh(theta) * c[imode], qn=[0,0]) + ham_terms.append(op) + op = Op(r"sigma_z x", ["spin", f"v_{imode}_q"], factor=np.sinh(theta) * c[imode], qn=[0,0]) + ham_terms.append(op) + + +nbas_factor = 2 +nbas = np.max([16 * c2/w**3 * np.cosh(theta_array)**2, np.ones(nmodes)*4], axis=0) +nbas = np.min([nbas, np.ones(nmodes)*512], axis=0) +nbas = np.round(nbas).astype(int) +nbas *= nbas_factor +logger.info(nbas) +basis = [ba.BasisHalfSpin("spin",[0,0])] +for imode in range(nmodes): + basis.append(ba.BasisSHO(f"v_{imode}_p", w[imode], int(nbas[imode]))) + basis.append(ba.BasisSHO(f"v_{imode}_q", w[imode], int(nbas[imode]))) + + +tree_order = 2 +basis_vib = basis[1:] +elementary_nodes = [] + +labels = np.array([[nbas>Ms], [nbas>Ms]]).T.ravel() +root = BasisTree.binary_mctdh(basis_vib, contract_primitive=True, contract_label=labels, dummy_label="n").root + +root.add_child(TreeNodeBasis(basis[:1])) + +basis_tree = BasisTree(root) +basis_tree.print(print_function=logger.info) + +# basis_tree = BasisTree.linear(basis) +ttno = TTNO(basis_tree, ham_terms) +exp_z = TTNO(basis_tree, Op("sigma_z", "spin")) +exp_x = TTNO(basis_tree, Op("sigma_x", "spin")) +ttns = TTNS(basis_tree) +ttns.compress_config = CompressConfig(CompressCriteria.fixed, max_bonddim=Ms) +ttns = expand_bond_dimension_general(ttns, ttno, ex_mps=None) +logger.info(ttns.bond_dims) +logger.info(ttno.bond_dims) +logger.info(len(ttns)) +ttns.evolve_config = EvolveConfig(EvolveMethod.tdvp_ps) +nsteps = 400 +dt = 0.1 +expectations = [] +for i in range(nsteps): + ttns = ttns.evolve(ttno, dt) + z = ttns.expectation(exp_z) + x = ttns.expectation(exp_x) + expectations.append((z, x)) + logger.info((z, x)) +logger.info(expectations) \ No newline at end of file diff --git a/example/ttns/sbm_zt.py b/example/ttns/sbm_zt.py new file mode 100644 index 00000000..0d9bd3b8 --- /dev/null +++ b/example/ttns/sbm_zt.py @@ -0,0 +1,107 @@ +import sys + +from renormalizer.model import Op +from renormalizer.model import basis as ba +from renormalizer.mps.mps import expand_bond_dimension_general +from renormalizer.sbm import ColeDavidsonSDF +from renormalizer.utils import EvolveConfig, CompressConfig, CompressCriteria, EvolveMethod +from renormalizer.utils import log +from renormalizer.tn import BasisTree, TTNO, TTNS, TreeNodeBasis + +import numpy as np + +ita_str = sys.argv[1] # 010, 025, 100 +omega_c_str = sys.argv[2] # 001, 100 +beta_str = sys.argv[3] # 025, 100 +# ita_str = "050" +# omega_c_str = "001" +# beta_str = "050" + +ita = int(ita_str) / 10 # 1, 2.5, 5, 10 +eps = 0 +Delta = 1 +omega_c = int(omega_c_str) / 10 # 0.1, 1, 10 + +beta = int(beta_str) / 100 # 0.25, 0.5, 0.75, 1 + +from renormalizer.utils.log import package_logger +logger = package_logger +dump_dir = "./" +job_name = f"ps1_binary_ita{ita_str}_omega{omega_c_str}_beta{beta_str}" #################### +log.register_file_output(dump_dir+job_name+".log", mode="w") + + +nmodes = 1000 +Ms = 20 +upper_limit = 30 +sdf = ColeDavidsonSDF(ita, omega_c, beta, upper_limit) + +w, c2 = sdf.Wang1(nmodes) +c = np.sqrt(c2) +logger.info(f"w:{w}") +logger.info(f"c:{c}") + +reno = sdf.reno(w[-1]) +logger.info(f"renormalization constant: {reno}") +Delta *= reno + +ham_terms = [] + +# h_s +ham_terms.extend([Op("sigma_z","spin",factor=eps, qn=0), + Op("sigma_x","spin",factor=Delta, qn=0)]) + + +# boson energy +for imode in range(nmodes): + op1 = Op(r"p^2",f"v_{imode}",factor=0.5, qn=0) + op2 = Op(r"x^2",f"v_{imode}",factor=0.5*w[imode]**2, qn=0) + ham_terms.extend([op1,op2]) + +# system-boson coupling +for imode in range(nmodes): + op = Op(r"sigma_z x", ["spin", f"v_{imode}"], + factor=c[imode], qn=[0,0]) + ham_terms.append(op) + +nbas = np.max([16 * c2/w**3, np.ones(nmodes)*4], axis=0) +nbas = np.round(nbas).astype(int) +logger.info(nbas) +basis = [ba.BasisHalfSpin("spin",[0,0])] +for imode in range(nmodes): + basis.append(ba.BasisSHO(f"v_{imode}", w[imode], int(nbas[imode]))) + + +tree_order = 2 +basis_vib = basis[1:] +elementary_nodes = [] + + +root = BasisTree.binary_mctdh(basis_vib, contract_primitive=True, contract_label=nbas>Ms, dummy_label="n").root + +root.add_child(TreeNodeBasis(basis[:1])) + +basis_tree = BasisTree(root) +basis_tree.print(print_function=logger.info) + +# basis_tree = BasisTree.linear(basis) +ttno = TTNO(basis_tree, ham_terms) +exp_z = TTNO(basis_tree, Op("sigma_z", "spin")) +exp_x = TTNO(basis_tree, Op("sigma_x", "spin")) +ttns = TTNS(basis_tree) +ttns.compress_config = CompressConfig(CompressCriteria.fixed, max_bonddim=Ms) +ttns = expand_bond_dimension_general(ttns, ttno, ex_mps=None) +logger.info(ttns.bond_dims) +logger.info(ttno.bond_dims) +logger.info(len(ttns)) +ttns.evolve_config = EvolveConfig(EvolveMethod.tdvp_ps) +nsteps = 200 +dt = 0.2 +expectations = [] +for i in range(nsteps): + ttns = ttns.evolve(ttno, dt) + z = ttns.expectation(exp_z) + x = ttns.expectation(exp_x) + expectations.append((z, x)) + logger.info((z, x)) +logger.info(expectations) \ No newline at end of file diff --git a/renormalizer/mps/mps.py b/renormalizer/mps/mps.py index c5a1b8e6..9f492f6a 100644 --- a/renormalizer/mps/mps.py +++ b/renormalizer/mps/mps.py @@ -590,64 +590,8 @@ def expand_bond_dimension(self, hint_mpo=None, coef=1e-10, include_ex=True): """ expand bond dimension as required in compress_config """ - # expander m target - m_target = self.compress_config.bond_dim_max_value - self.bond_dims_mean - # will be restored at exit - self.compress_config.bond_dim_max_value = m_target - if self.compress_config.criteria is not CompressCriteria.fixed: - logger.warning("Setting compress criteria to fixed") - self.compress_config.criteria = CompressCriteria.fixed - logger.debug(f"target for expander: {m_target}") - if hint_mpo is None: - expander = self.__class__.random(self.model, 1, m_target) - else: - # fill states related to `hint_mpo` - logger.debug( - f"average bond dimension of hint mpo: {hint_mpo.bond_dims_mean}" - ) - # in case of localized `self` - if include_ex: - if self.is_mps: - ex_state: MatrixProduct = self.ground_state(self.model, False) - # for self.qntot >= 1 - assert self.model.qn_size == 1 # otherwise not supported - for i in range(self.qntot[0]): - ex_state = Mpo.onsite(self.model, r"a^\dagger") @ ex_state - elif self.is_mpdm: - assert self.qntot == 1 - ex_state: MatrixProduct = self.max_entangled_ex(self.model) - else: - assert False - ex_state.compress_config = self.compress_config - ex_state.move_qnidx(self.qnidx) - ex_state.to_right = self.to_right - lastone = self + ex_state + return expand_bond_dimension(self, hint_mpo, coef, include_ex) - else: - lastone = self - expander_list: List["MatrixProduct"] = [] - cumulated_m = 0 - while True: - lastone.compress_config.criteria = CompressCriteria.fixed - expander_list.append(lastone) - expander = compressed_sum(expander_list) - if cumulated_m == expander.bond_dims_mean: - # probably a small system, the required bond dimension can't be reached - break - cumulated_m = expander.bond_dims_mean - logger.debug( - f"cumulated bond dimension: {cumulated_m}. lastone bond dimension: {lastone.bond_dims}" - ) - if m_target < cumulated_m: - break - if m_target < 0.8 * (lastone.bond_dims_mean * hint_mpo.bond_dims_mean): - lastone = lastone.canonicalise().compress( - m_target // hint_mpo.bond_dims_mean + 1 - ) - lastone = (hint_mpo @ lastone).normalize("mps_and_coeff") - logger.debug(f"expander bond dimension: {expander.bond_dims}") - self.compress_config.bond_dim_max_value += self.bond_dims_mean - return (self + expander.scale(coef*self.norm, inplace=True)).canonicalise().canonicalise().normalize("mps_norm_to_coeff") def evolve(self, mpo, evolve_dt, normalize=True) -> "Mps": @@ -1834,6 +1778,7 @@ def projector( proj = Iden - proj return proj + def integrand_func_factory( shape, hop, @@ -1919,6 +1864,86 @@ def _mu_regularize(s, epsilon=1e-10): return s + epsilon * np.exp(-s / epsilon) +def expand_bond_dimension(mps, hint_mpo=None, coef=1e-10, include_ex=True): + """ + expand bond dimension as required in compress_config + """ + if hint_mpo is not None and include_ex: + # fill states related to `hint_mpo` + logger.debug( + f"average bond dimension of hint mpo: {hint_mpo.bond_dims_mean}" + ) + # in case of localized `self` + if mps.is_mps: + ex_state: MatrixProduct = mps.ground_state(mps.model, False) + # for self.qntot >= 1 + assert mps.model.qn_size == 1 # otherwise not supported + for i in range(mps.qntot[0]): + ex_state = Mpo.onsite(mps.model, r"a^\dagger") @ ex_state + elif mps.is_mpdm: + assert mps.qntot == 1 + ex_state: MatrixProduct = mps.max_entangled_ex(mps.model) + else: + assert False + ex_state.compress_config = mps.compress_config + ex_state.move_qnidx(mps.qnidx) + ex_state.to_right = mps.to_right + else: + ex_state = None + return expand_bond_dimension_general(mps, hint_mpo, coef, ex_state) + + +def expand_bond_dimension_general(mps, hint_mpo=None, coef=1e-10, ex_mps=None): + """ + expand bond dimension as required in compress_config. works for both mps and ttns + """ + # expander m target + m_target = mps.compress_config.bond_dim_max_value - mps.bond_dims_mean + # will be restored at exit + mps.compress_config.bond_dim_max_value = m_target + if mps.compress_config.criteria is not CompressCriteria.fixed: + logger.warning("Setting compress criteria to fixed") + mps.compress_config.criteria = CompressCriteria.fixed + logger.debug(f"target for expander: {m_target}") + if hint_mpo is None: + expander = mps.__class__.random(mps.model, 1, m_target) + else: + # fill states related to `hint_mpo` + logger.debug( + f"average bond dimension of hint mpo: {hint_mpo.bond_dims_mean}" + ) + # in case of localized `self` + if ex_mps: + lastone = mps + ex_mps + + else: + lastone = mps + expander_list: List["MatrixProduct"] = [] + cumulated_m = 0 + while True: + lastone.compress_config.criteria = CompressCriteria.fixed + expander_list.append(lastone) + expander = compressed_sum(expander_list) + if cumulated_m == expander.bond_dims_mean: + # probably a small system, the required bond dimension can't be reached + break + cumulated_m = expander.bond_dims_mean + logger.debug( + f"cumulated bond dimension: {cumulated_m}. lastone bond dimension: {lastone.bond_dims}" + ) + if m_target < cumulated_m: + break + if m_target < 0.8 * (lastone.bond_dims_mean * hint_mpo.bond_dims_mean): + lastone = lastone.canonicalise().compress( + m_target // hint_mpo.bond_dims_mean + 1 + ) + lastone = (hint_mpo @ lastone).normalize("mps_and_coeff") + logger.debug(f"expander bond dimension: {expander.bond_dims}") + mps.compress_config.bond_dim_max_value += mps.bond_dims_mean + return (mps + expander.scale(coef * mps.norm, inplace=True)).canonicalise().canonicalise().normalize( + "mps_norm_to_coeff") + + class BraKetPair: def __init__(self, bra_mps, ket_mps, mpo=None): self.bra_mps = bra_mps diff --git a/renormalizer/mps/tests/test_mpo.py b/renormalizer/mps/tests/test_mpo.py index cd73d3ae..9e6a8774 100644 --- a/renormalizer/mps/tests/test_mpo.py +++ b/renormalizer/mps/tests/test_mpo.py @@ -129,6 +129,7 @@ def test_scheme4(): mps3.append(np.array([0, 1]).reshape((1,2,1))) mps3.append(np.array([0.707, 0.707]).reshape((1,2,1))) mps3.append(np.array([1, 0]).reshape((1,2,1))) + mps3.build_empty_qn() e3 = mps3.expectation(mpo3) assert pytest.approx(e4) == e3 diff --git a/renormalizer/sbm/__init__.py b/renormalizer/sbm/__init__.py index d8eabfa6..3b1bf1e2 100644 --- a/renormalizer/sbm/__init__.py +++ b/renormalizer/sbm/__init__.py @@ -1,4 +1,4 @@ # -*- coding: utf-8 -*- -from renormalizer.sbm.lib import SpectralDensityFunction, param2mollist +from renormalizer.sbm.lib import SpectralDensityFunction, param2mollist, DebyeSDF, OhmicSDF, ColeDavidsonSDF from renormalizer.sbm.sbm import SpinBosonDynamics \ No newline at end of file diff --git a/renormalizer/sbm/lib.py b/renormalizer/sbm/lib.py index 59dd7f9f..f76c8f83 100644 --- a/renormalizer/sbm/lib.py +++ b/renormalizer/sbm/lib.py @@ -1,5 +1,6 @@ # -*- coding: utf-8 -*- +import logging import numpy as np import numpy.polynomial.laguerre as la import numpy.polynomial.legendre as le @@ -7,10 +8,13 @@ import scipy.special import scipy.optimize -from renormalizer.model import Phonon, Mol, SpinBosonModel +from renormalizer.model import Phonon, SpinBosonModel from renormalizer.utils import Quantity +logger = logging.getLogger(__name__) + + class DebyeSpectralDensityFunction: r""" the Debye-type ohmic spectral density function @@ -28,6 +32,9 @@ def func(self, omega_value): return 2. * self.lamb * omega_value * self.omega_c / (omega_value ** 2 + self.omega_c ** 2) +DebyeSDF = DebyeSpectralDensityFunction + + class SpectralDensityFunction: r""" the ohmic spectral density function @@ -154,6 +161,75 @@ def plot_data(self, x0, x1, n, omega_value, c_j2, sigma=0.1): return x, y_c, y_d +OhmicSDF = SpectralDensityFunction + + +class ColeDavidsonSDF: + """ + the ColeDavidson spectral density function + """ + + def __init__(self, ita, omega_c, beta, omega_limit): + self.ita = ita + self.omega_c = omega_c + self.beta = beta + self.omega_limit = omega_limit + + + def reno(self, omega_l): + def integrate_func(x): + return self.func(x) / x**2 + + res = scipy.integrate.quad(integrate_func, a=omega_l, + b=omega_l*1000) + logger.info(f"integrate: {res[0]}, {res[1]}") + re = np.exp(-res[0]*2/np.pi) + + return re + + + def func(self, omega_value): + """ + the function of the spectral density function + """ + theta = np.arctan(omega_value/self.omega_c) + return self.ita * np.sin(self.beta * theta) / (1 + omega_value**2/self.omega_c**2) ** (self.beta / 2) + + + def _dos_Wang1(self, A, omega_value): + """ + Wang's 1st scheme DOS \rho(\omega) + """ + + return A * self.func(omega_value) / omega_value + + def Wang1(self, nb): + """ + Wang's 1st scheme discretization + """ + def integrate_func(x): + return self.func(x) / x + A = (nb + 1 ) / scipy.integrate.quad(integrate_func, a=0, b=self.omega_limit)[0] + logger.info(scipy.integrate.quad(integrate_func, a=0, b=self.omega_limit)[0] * 4 / np.pi) + logger.info(2*self.ita) + nsamples = int(1e7) + delta = self.omega_limit / nsamples + omega_value_big = np.linspace(delta, self.omega_limit, nsamples) + dos = self._dos_Wang1(A, omega_value_big) + rho_cumint = np.cumsum(dos) * delta + diff = (rho_cumint % 1)[1:] - (rho_cumint % 1)[:-1] + idx = np.where(diff < 0)[0] + omega_value = omega_value_big[idx] + logger.info(len(omega_value)) + assert len(omega_value) == nb + + # general form + c_j2 = 2./np.pi * omega_value * self.func(omega_value) / self._dos_Wang1(A, omega_value) + + + return omega_value, c_j2 + + def param2mollist(alpha: float, raw_delta: Quantity, omega_c: Quantity, renormalization_p: float, n_phonons: int): sdf = SpectralDensityFunction(alpha, omega_c) delta, max_omega = sdf.adiabatic_renormalization(raw_delta, renormalization_p) diff --git a/renormalizer/sbm/sbm.py b/renormalizer/sbm/sbm.py index a10e3521..79f783e9 100644 --- a/renormalizer/sbm/sbm.py +++ b/renormalizer/sbm/sbm.py @@ -1,12 +1,11 @@ # -*- coding: utf-8 -*- import logging -from functools import partial from renormalizer.model import Model from renormalizer.mps import Mpo, Mps -from renormalizer.utils import TdMpsJob, Quantity, CompressConfig -import numpy as np +from renormalizer.utils import TdMpsJob, CompressConfig + logger = logging.getLogger(__name__) diff --git a/renormalizer/tn/hop_expr.py b/renormalizer/tn/hop_expr.py index 2895fedd..6ea850ac 100644 --- a/renormalizer/tn/hop_expr.py +++ b/renormalizer/tn/hop_expr.py @@ -1,6 +1,6 @@ import opt_einsum as oe -from renormalizer.mps.backend import np +from renormalizer.mps.backend import np, backend from renormalizer.mps.matrix import asxp from renormalizer.tn.node import TreeNodeTensor from renormalizer.tn.tree import TTNS, TTNO, TTNEnviron @@ -42,8 +42,6 @@ def hop_expr0(snode: TreeNodeTensor, ttns: TTNS, ttno: TTNO, ttne: TTNEnviron): return expr - - def hop_expr1(snode: TreeNodeTensor, ttns: TTNS, ttno: TTNO, ttne: TTNEnviron, return_hdiag=False): # one site enode = ttne.node_list[ttns.node_idx[snode]] @@ -62,7 +60,7 @@ def hop_expr1(snode: TreeNodeTensor, ttns: TTNS, ttno: TTNO, ttne: TTNEnviron, r # input and output input_indices = ttns.get_node_indices(snode) - output_indices = ttns.get_node_indices(snode, True) + output_indices = ttns.get_node_indices(snode, conj=True) shape = snode.shape # cache the contraction path @@ -105,7 +103,7 @@ def hop_expr2(snode: TreeNodeTensor, ttns: TTNS, ttno: TTNO, ttne: TTNEnviron): # input and output input_indices = ttns.get_node_indices(snode, include_parent=True) - output_indices = ttns.get_node_indices(snode, True, include_parent=True) + output_indices = ttns.get_node_indices(snode, conj=True, include_parent=True) # shape shape = list(snode.shape[:-1]) diff --git a/renormalizer/tn/node.py b/renormalizer/tn/node.py index b45175cd..9d40a678 100644 --- a/renormalizer/tn/node.py +++ b/renormalizer/tn/node.py @@ -99,6 +99,21 @@ def __init__(self): def copy_connection(source_node_list: List[NodeUnion], target_node_list: List[NodeUnion]) -> NodeUnion: + """ + Create a new tree with the same connection structure as the source tree in the target tree. + + Parameters + ---------- + source_node_list : List[NodeUnion] + The list of nodes in the source tree. + target_node_list : List[NodeUnion] + The list of nodes in the target tree. + + Returns + ------- + NodeUnion + The root node of the target tree. + """ node2idx: Dict[NodeUnion, int] = {n:i for i, n in enumerate(source_node_list)} root = None for source_node, target_node in zip(source_node_list, target_node_list): diff --git a/renormalizer/tn/tree.py b/renormalizer/tn/tree.py index 8c2a2897..b95dc609 100644 --- a/renormalizer/tn/tree.py +++ b/renormalizer/tn/tree.py @@ -11,6 +11,7 @@ from renormalizer.mps.svd_qn import add_outer, svd_qn, blockrecover, get_qn_mask from renormalizer.mps.lib import select_basis from renormalizer.utils.configs import CompressConfig, OptimizeConfig, EvolveConfig, EvolveMethod +from renormalizer.utils import calc_vn_entropy from renormalizer.tn.node import TreeNodeTensor, TreeNodeBasis, copy_connection, TreeNodeEnviron from renormalizer.tn.treebase import Tree, BasisTree from renormalizer.tn.symbolic_mpo import construct_symbolic_mpo, symbolic_mo_to_numeric_mo_general @@ -19,7 +20,7 @@ class TTNBase(Tree): # A tree whose tree node is TreeNodeTensor def print_shape(self, full=False, print_function=None): - class print_tn_basis(print_tree): + class _print_shape(print_tree): def get_children(self, node): return node.children @@ -30,11 +31,43 @@ def get_node_str(self, node: TreeNodeTensor): else: return str(node.tensor.shape[-1]) - tree = print_tn_basis(self.root) + tree = _print_shape(self.root) if print_function is not None: for row in tree.rows: print_function(row) + def print_vn_entropy(self, print_function=None): + # todo: move ttns.compress to here + _, s_list = self.compress(ret_s=True) + vn_entropy = [calc_vn_entropy(s ** 2) for s in s_list] + nodes = [TreeNodeTensor([entropy]) for entropy in vn_entropy] + copy_connection(self.node_list, nodes) + + class print_data(print_tree): + + def get_children(self, node): + return node.children + + def get_node_str(self, node: TreeNodeTensor): + return str(node.tensor.shape[-1]) + + tree = print_data(nodes[0]) + if print_function is not None: + for row in tree.rows: + print_function(row) + + @property + def bond_dims(self): + return [node.tensor.shape[-1] for node in self] + + @property + def bond_dims_mean(self) -> int: + # duplicate with Mps + return int(round(np.mean(self.bond_dims))) + + @property + def qntot(self): + return self.root.qn[0] class TTNO(TTNBase): @classmethod @@ -155,21 +188,6 @@ def get_node_indices(self, node, prefix_up, prefix_down): assert len(indices) == node.tensor.ndim return indices - @property - def bond_dims(self): - # duplicate with ttns - return [node.tensor.shape[-1] for node in self] - - @property - def bond_dims_mean(self) -> int: - # duplicate with ttns and Mps - return int(round(np.mean(self.bond_dims))) - - @property - def qntot(self): - # duplicate with ttns - return self.root.qn[0] - def __matmul__(self, other): # duplicate with Mpo return self.apply(other) @@ -705,14 +723,9 @@ def norm(self): ''' return np.linalg.norm(self.coeff) * self.ttns_norm - @property - def qntot(self): - # duplicate with ttno - return self.root.qn[0] - @property def ttns_norm(self): - res = self.expectation(TTNO.identity(self.basis)) + res = self.expectation(TTNO.identity(self.basis)).real if res < 0: assert np.abs(res) < 1e-8 @@ -720,15 +733,6 @@ def ttns_norm(self): res = np.sqrt(res) return float(res) - @property - def bond_dims(self): - return [node.tensor.shape[-1] for node in self] - - @property - def bond_dims_mean(self) -> int: - # duplicate with Mps - return int(round(np.mean(self.bond_dims))) - def scale(self, val, inplace=False): # no need to be canonical # self.check_canonical() @@ -755,7 +759,7 @@ def __init__(self, ttns:TTNS, ttno:TTNO, build_environ=True): copy_connection(ttns.node_list, enodes) super().__init__(enodes[0]) assert self.root.parent is None - self.root.environ_parent = np.array([1]).reshape([1, 1, 1]) + self.root.environ_parent = np.array([1], dtype=backend.real_dtype).reshape([1, 1, 1]) # tensor node to basis node. todo: remove duplication? self.tn2bn = {tn: bn for tn, bn in zip(self.node_list, self.basis.node_list)} self.tn2dofs = {tn: bn.dofs for tn, bn in self.tn2bn.items()} @@ -799,7 +803,7 @@ def update_2site(self, snode, ttns, ttno): self.build_parent_environ_node(snode, ichild, ttns, ttno) def build_children_environ_node(self, snode:TreeNodeTensor, ttns: TTNS, ttno: TTNO): - # build the environment from snode to its parent + # build the environment from snode to its parent and store the environment in its parent if snode.parent is None: return enode = self.node_list[ttns.node_idx[snode]] @@ -831,7 +835,7 @@ def build_children_environ_node(self, snode:TreeNodeTensor, ttns: TTNS, ttno: TT enode.parent.environ_children[ichild] = asnumpy(res) def build_parent_environ_node(self, snode:TreeNodeTensor, ichild: int, ttns: TTNS, ttno: TTNO): - # build the environment from snode to the ith child of snode + # build the environment from snode to the ith child of snode and store the environment in the child enode = self.node_list[ttns.node_idx[snode]] onode = ttno.node_list[ttns.node_idx[snode]] args = [] diff --git a/renormalizer/utils/quantity.py b/renormalizer/utils/quantity.py index 8f94082b..fce2e110 100644 --- a/renormalizer/utils/quantity.py +++ b/renormalizer/utils/quantity.py @@ -4,6 +4,7 @@ from __future__ import division +import math import logging from renormalizer.utils import constant @@ -49,6 +50,8 @@ def as_unit(self, unit): # kelvin to beta au^-1 def to_beta(self): + if self.value == 0: + return math.inf return 1.0 / self.as_au() # a simplified and incomplete model for + - * / diff --git a/setup.py b/setup.py index 79c600b1..b6e571da 100644 --- a/setup.py +++ b/setup.py @@ -6,7 +6,8 @@ req = ["numpy", "scipy", "h5py", - "opt_einsum" + "opt_einsum", + "sympy" ] setuptools.setup(