Numerical construction of MPO
[1]:
import numpy as np
import sympy
from pympo import (
AssignManager,
OpSite,
SumOfProducts,
)
[2]:
HenonHeiles = SumOfProducts([])
omega = sympy.Symbol("omega")
lam = sympy.Symbol("lambda")
ndim = 4
ints_p2 = np.ones((3, 3))
ints_q = np.ones((3,))
ints_q2 = np.ones((3,))
ints_q3 = np.ones((3,))
p2_ops = [OpSite(r"\hat{p}^2" + f"_{i}", i, value=ints_p2) for i in range(ndim)]
q2_ops = [OpSite(r"\hat{q}^2" + f"_{i}", i, value=ints_q2) for i in range(ndim)]
q_ops = [OpSite(r"\hat{q}" + f"_{i}", i, value=ints_q) for i in range(ndim)]
q3_ops = [OpSite(r"\hat{q}^3" + f"_{i}", i, value=ints_q3) for i in range(ndim)]
for i in range(ndim):
HenonHeiles += p2_ops[i] * (omega / 2)
HenonHeiles += q2_ops[i] * (omega / 2)
for i in range(ndim - 1):
HenonHeiles += q2_ops[i] * q_ops[i + 1] * lam
HenonHeiles -= q3_ops[i + 1] * (lam / 3)
HenonHeiles = HenonHeiles.simplify()
print(f"{HenonHeiles.ndim=}, {HenonHeiles.nops=}")
display(HenonHeiles.symbol)
HenonHeiles.ndim=4, HenonHeiles.nops=10
$\displaystyle \lambda \hat{q}^2_0 \hat{q}_1 + \lambda \hat{q}^2_1 \hat{q}_2 + \lambda \hat{q}^2_2 \hat{q}_3 - \frac{\lambda \hat{q}^3_1}{3} - \frac{\lambda \hat{q}^3_2}{3} - \frac{\lambda \hat{q}^3_3}{3} + \frac{\omega \left(\hat{p}^2_0 + \hat{q}^2_0\right)}{2} + \frac{\omega \left(\hat{p}^2_1 + \hat{q}^2_1\right)}{2} + \frac{\omega \left(\hat{p}^2_2 + \hat{q}^2_2\right)}{2} + \frac{\omega \left(\hat{p}^2_3 + \hat{q}^2_3\right)}{2}$
[3]:
am = AssignManager(HenonHeiles)
display(*am.Wsym)
None
None
None
None
[4]:
am.assign()
display(*am.Wsym)
2024-12-31 18:07:22.486 | INFO | pympo.bipartite:assign:286 - assigned 1/4
2024-12-31 18:07:22.487 | INFO | pympo.bipartite:assign:286 - assigned 2/4
2024-12-31 18:07:22.488 | INFO | pympo.bipartite:assign:286 - assigned 3/4
2024-12-31 18:07:22.489 | INFO | pympo.bipartite:assign:286 - assigned 4/4
$\displaystyle \left[\begin{matrix}\hat{q}^2_0 & 1 & \hat{p}^2_0 + \hat{q}^2_0\end{matrix}\right]$
$\displaystyle \left[\begin{matrix}0 & 0 & \lambda \hat{q}_1\\\hat{q}^2_1 & 1 & - \frac{\lambda \hat{q}^3_1}{3} + \frac{\omega \left(\hat{p}^2_1 + \hat{q}^2_1\right)}{2}\\0 & 0 & \frac{\omega}{2}\end{matrix}\right]$
$\displaystyle \left[\begin{matrix}0 & 0 & \lambda \hat{q}_2\\\hat{q}^2_2 & 1 & - \frac{\lambda \hat{q}^3_2}{3} + \frac{\omega \left(\hat{p}^2_2 + \hat{q}^2_2\right)}{2}\\0 & 0 & 1\end{matrix}\right]$
$\displaystyle \left[\begin{matrix}\lambda \hat{q}_3\\- \frac{\lambda \hat{q}^3_3}{3} + \frac{\omega \left(\hat{p}^2_3 + \hat{q}^2_3\right)}{2}\\1\end{matrix}\right]$
[5]:
W_prod = sympy.Mul(*am.Wsym)
print(*[f"W{i}" for i in range(am.ndim)], "=")
display(W_prod[0].expand())
W0 W1 W2 W3 =
$\displaystyle \lambda \hat{q}^2_0 \hat{q}_1 + \lambda \hat{q}^2_1 \hat{q}_2 + \lambda \hat{q}^2_2 \hat{q}_3 - \frac{\lambda \hat{q}^3_1}{3} - \frac{\lambda \hat{q}^3_2}{3} - \frac{\lambda \hat{q}^3_3}{3} + \frac{\omega \hat{p}^2_0}{2} + \frac{\omega \hat{p}^2_1}{2} + \frac{\omega \hat{p}^2_2}{2} + \frac{\omega \hat{p}^2_3}{2} + \frac{\omega \hat{q}^2_0}{2} + \frac{\omega \hat{q}^2_1}{2} + \frac{\omega \hat{q}^2_2}{2} + \frac{\omega \hat{q}^2_3}{2}$
[6]:
am.show_graph()
[7]:
mpo = am.numerical_mpo(
subs={omega: 2.0, lam: 1.0}
) # subs is not necessary if you define its value in advance.
mpo
[7]:
[array([[[[1.+0.j, 1.+0.j, 2.+0.j],
[0.+0.j, 0.+0.j, 1.+0.j],
[0.+0.j, 0.+0.j, 1.+0.j]],
[[0.+0.j, 0.+0.j, 1.+0.j],
[1.+0.j, 1.+0.j, 2.+0.j],
[0.+0.j, 0.+0.j, 1.+0.j]],
[[0.+0.j, 0.+0.j, 1.+0.j],
[0.+0.j, 0.+0.j, 1.+0.j],
[1.+0.j, 1.+0.j, 2.+0.j]]]]),
array([[[[0. +0.j, 0. +0.j, 1. +0.j],
[0. +0.j, 0. +0.j, 0. +0.j],
[0. +0.j, 0. +0.j, 0. +0.j]],
[[0. +0.j, 0. +0.j, 0. +0.j],
[0. +0.j, 0. +0.j, 1. +0.j],
[0. +0.j, 0. +0.j, 0. +0.j]],
[[0. +0.j, 0. +0.j, 0. +0.j],
[0. +0.j, 0. +0.j, 0. +0.j],
[0. +0.j, 0. +0.j, 1. +0.j]]],
[[[1. +0.j, 1. +0.j, 1.66666667+0.j],
[0. +0.j, 0. +0.j, 1. +0.j],
[0. +0.j, 0. +0.j, 1. +0.j]],
[[0. +0.j, 0. +0.j, 1. +0.j],
[1. +0.j, 1. +0.j, 1.66666667+0.j],
[0. +0.j, 0. +0.j, 1. +0.j]],
[[0. +0.j, 0. +0.j, 1. +0.j],
[0. +0.j, 0. +0.j, 1. +0.j],
[1. +0.j, 1. +0.j, 1.66666667+0.j]]],
[[[0. +0.j, 0. +0.j, 1. +0.j],
[0. +0.j, 0. +0.j, 0. +0.j],
[0. +0.j, 0. +0.j, 0. +0.j]],
[[0. +0.j, 0. +0.j, 0. +0.j],
[0. +0.j, 0. +0.j, 1. +0.j],
[0. +0.j, 0. +0.j, 0. +0.j]],
[[0. +0.j, 0. +0.j, 0. +0.j],
[0. +0.j, 0. +0.j, 0. +0.j],
[0. +0.j, 0. +0.j, 1. +0.j]]]]),
array([[[[0. +0.j, 0. +0.j, 1. +0.j],
[0. +0.j, 0. +0.j, 0. +0.j],
[0. +0.j, 0. +0.j, 0. +0.j]],
[[0. +0.j, 0. +0.j, 0. +0.j],
[0. +0.j, 0. +0.j, 1. +0.j],
[0. +0.j, 0. +0.j, 0. +0.j]],
[[0. +0.j, 0. +0.j, 0. +0.j],
[0. +0.j, 0. +0.j, 0. +0.j],
[0. +0.j, 0. +0.j, 1. +0.j]]],
[[[1. +0.j, 1. +0.j, 1.66666667+0.j],
[0. +0.j, 0. +0.j, 1. +0.j],
[0. +0.j, 0. +0.j, 1. +0.j]],
[[0. +0.j, 0. +0.j, 1. +0.j],
[1. +0.j, 1. +0.j, 1.66666667+0.j],
[0. +0.j, 0. +0.j, 1. +0.j]],
[[0. +0.j, 0. +0.j, 1. +0.j],
[0. +0.j, 0. +0.j, 1. +0.j],
[1. +0.j, 1. +0.j, 1.66666667+0.j]]],
[[[0. +0.j, 0. +0.j, 1. +0.j],
[0. +0.j, 0. +0.j, 0. +0.j],
[0. +0.j, 0. +0.j, 0. +0.j]],
[[0. +0.j, 0. +0.j, 0. +0.j],
[0. +0.j, 0. +0.j, 1. +0.j],
[0. +0.j, 0. +0.j, 0. +0.j]],
[[0. +0.j, 0. +0.j, 0. +0.j],
[0. +0.j, 0. +0.j, 0. +0.j],
[0. +0.j, 0. +0.j, 1. +0.j]]]]),
array([[[[1. +0.j],
[0. +0.j],
[0. +0.j]],
[[0. +0.j],
[1. +0.j],
[0. +0.j]],
[[0. +0.j],
[0. +0.j],
[1. +0.j]]],
[[[1.66666667+0.j],
[1. +0.j],
[1. +0.j]],
[[1. +0.j],
[1.66666667+0.j],
[1. +0.j]],
[[1. +0.j],
[1. +0.j],
[1.66666667+0.j]]],
[[[1. +0.j],
[0. +0.j],
[0. +0.j]],
[[0. +0.j],
[1. +0.j],
[0. +0.j]],
[[0. +0.j],
[0. +0.j],
[1. +0.j]]]])]
[8]:
dump = None
for core in mpo:
if dump is None:
dump = core[0, 0, 0, :]
else:
dump = np.einsum("i,ij->j", dump, core[:, 0, 0, :])