Execute Quantum Dynamics with PyTDSCF –Taylor expansion PES
[1]:
try:
import pytdscf
except ModuleNotFoundError:
!uv pip install -U git+https://github.com/QCLovers/PyTDSCF --quiet
try:
import discvar
except ModuleNotFoundError:
!uv pip install -U git+https://github.com/QCLovers/Discvar --quiet
/Users/hinom/GitHub/PyMPO/.venv/lib/python3.12/site-packages/tqdm/auto.py:21: TqdmWarning: IProgress not found. Please update jupyter and ipywidgets. See https://ipywidgets.readthedocs.io/en/stable/user_install.html
from .autonotebook import tqdm as notebook_tqdm
[2]:
import itertools
from collections import defaultdict
from math import factorial, sqrt
import numpy as np
import sympy
from discvar import HarmonicOscillator as HO
from pytdscf import Model, Simulator, units
from pytdscf.potentials.ch2o_potential import k_orig
import pympo
[3]:
backend = "jax"
m = 30
N = 9
Δt = 0.1
cutoff = 0.0 # 1.e-09
\[V - V_0 = \frac{k_{11}}{2!} Q_1^2 + \frac{k_{22}}{2!} Q_2^2 + \frac{k_{33}}{2!} Q_3^2 + \frac{k_{111}}{3!} Q_1^3 + \frac{k_{122}}{2!} Q_1Q_2^2 + \cdots\]
[4]:
print(f"{len(k_orig)=}")
active_modes = sorted(list(set(itertools.chain.from_iterable(k_orig.keys()))))
index = {}
for i, mode in enumerate(active_modes):
index[mode] = i
print(f"{index=}")
f = len(active_modes)
k_new = {}
for key, value in k_orig.items():
new_key = tuple([index[k] for k in key])
if abs(value) > cutoff:
k_new[new_key] = value
print(f"{len(k_new)=}")
len(k_orig)=114
index={1: 0, 2: 1, 3: 2, 4: 3, 5: 4, 6: 5}
len(k_new)=114
[5]:
k_new
[5]:
{(0, 0): 2.881141028397287e-05,
(1, 1): 3.2448494255464496e-05,
(2, 2): 4.725466264016419e-05,
(3, 3): 6.805132227202656e-05,
(4, 4): 0.00017348287106252062,
(5, 5): 0.0001808886981537931,
(2, 2, 2): -1.1756574314828089e-07,
(3, 3, 3): 1.956730822300732e-06,
(4, 4, 4): 9.330094468538158e-06,
(5, 5, 5): -9.893510625593033e-09,
(0, 0, 2): -1.238616135463855e-07,
(1, 1, 2): 2.3333266618281752e-07,
(0, 3, 3): 7.709229058903663e-10,
(0, 0, 3): 6.899760007718777e-08,
(1, 1, 3): 4.7154784410294076e-08,
(2, 3, 3): 3.410048987055054e-07,
(2, 2, 3): 2.858839109343442e-07,
(0, 4, 4): 2.955204472579738e-09,
(0, 0, 4): -9.109739004604497e-07,
(1, 1, 4): -6.881771806581336e-07,
(2, 4, 4): -7.439406041842036e-08,
(2, 2, 4): -3.3560843836427284e-07,
(3, 4, 4): -9.31531844617526e-08,
(3, 3, 4): -2.9128037127557674e-07,
(0, 5, 5): 3.2121787745431932e-09,
(2, 5, 5): 6.435921392674742e-07,
(3, 5, 5): -6.933166666974028e-07,
(4, 5, 5): 1.0040242952014167e-05,
(4, 4, 5): 3.4691530765066486e-09,
(0, 1, 4): -5.13948603926911e-10,
(0, 3, 4): -6.424357549086387e-10,
(2, 3, 4): 2.4720927848884413e-07,
(0, 1, 5): -1.798820113744188e-09,
(0, 2, 5): 1.027897207853822e-09,
(1, 2, 5): 6.678762108030208e-07,
(0, 3, 5): -6.424357549086386e-10,
(1, 3, 5): -4.840110977481683e-07,
(1, 4, 5): -2.980901902776083e-08,
(0, 0, 0, 0): 1.7939020286070263e-08,
(1, 1, 1, 1): 1.4535391374218987e-08,
(2, 2, 2, 2): 5.3808032664810655e-09,
(3, 3, 3, 3): 5.192114024418222e-08,
(4, 4, 4, 4): 4.229558138073721e-07,
(5, 5, 5, 5): 5.069812097824381e-07,
(0, 1, 1, 1): 8.727253620131479e-11,
(0, 0, 1, 1): 4.207138124463382e-09,
(0, 0, 0, 1): 1.5046989000226694e-10,
(0, 0, 2, 2): 2.2841329302344115e-09,
(1, 1, 2, 2): 8.314966121525269e-09,
(0, 3, 3, 3): 2.4075182400362706e-11,
(0, 0, 3, 3): -2.317236306034911e-09,
(1, 1, 3, 3): -9.028193400136016e-11,
(2, 3, 3, 3): 9.928003342349572e-09,
(2, 2, 3, 3): 4.613406827469504e-09,
(2, 2, 2, 3): -3.2351026350487387e-09,
(0, 4, 4, 4): 1.594980834024029e-10,
(0, 0, 4, 4): -8.878024449913753e-08,
(0, 0, 0, 4): -7.222554720108814e-11,
(1, 4, 4, 4): -2.7084580200408045e-11,
(1, 1, 4, 4): -7.305012219830055e-08,
(2, 4, 4, 4): -4.9173560052740825e-09,
(2, 2, 4, 4): -6.238180699713982e-08,
(2, 2, 2, 4): -5.537291952083423e-09,
(3, 4, 4, 4): -6.918605542304234e-09,
(3, 3, 4, 4): -1.3683731796806156e-08,
(3, 3, 3, 4): -4.526134291268189e-09,
(0, 5, 5, 5): -3.9122171400589393e-10,
(0, 0, 5, 5): -1.0699612938281198e-07,
(0, 0, 0, 5): -8.727253620131484e-11,
(1, 5, 5, 5): -3.2606825163491245e-08,
(1, 1, 5, 5): -7.859042354818399e-08,
(1, 1, 1, 5): -4.568265860468825e-09,
(2, 2, 5, 5): -8.568357416289087e-08,
(3, 3, 5, 5): -2.4649977380171364e-08,
(4, 5, 5, 5): -2.828833932042618e-10,
(4, 4, 5, 5): 4.6045591979373705e-07,
(4, 4, 4, 5): 3.220055646048512e-10,
(0, 1, 2, 2): 9.329133180140548e-11,
(0, 1, 1, 2): -2.1065784600317367e-11,
(0, 1, 3, 3): 2.4075182400362706e-11,
(0, 1, 1, 3): 1.2037591200181356e-11,
(0, 0, 2, 3): -2.287142328034457e-09,
(1, 1, 2, 3): -4.083752814661525e-09,
(0, 1, 4, 4): -8.426313840126947e-11,
(0, 1, 1, 4): -2.7084580200408045e-11,
(0, 2, 4, 4): -1.2037591200181355e-11,
(0, 2, 2, 4): -2.1065784600317367e-11,
(0, 0, 2, 4): -2.949209844044432e-09,
(1, 2, 4, 4): 4.514096700068007e-11,
(1, 1, 2, 4): 6.503308645897975e-09,
(0, 3, 4, 4): -1.8056386800272028e-11,
(0, 0, 3, 4): 5.29954952587984e-09,
(1, 3, 4, 4): -2.1065784600317367e-11,
(1, 1, 3, 4): 3.737672067656311e-09,
(2, 3, 4, 4): 3.1707015221277685e-08,
(2, 3, 3, 4): -5.359737481880746e-09,
(2, 2, 3, 4): 7.472334737512576e-09,
(0, 1, 5, 5): -1.0833832080163219e-10,
(0, 0, 1, 5): 3.731653272056219e-10,
(0, 2, 5, 5): -1.2037591200181355e-11,
(1, 2, 5, 5): -4.8150364800725425e-11,
(1, 2, 2, 5): 5.109957464476984e-09,
(0, 3, 5, 5): -2.4075182400362706e-11,
(1, 3, 5, 5): 2.407518240036271e-11,
(1, 3, 3, 5): 2.6001196992391727e-09,
(2, 3, 5, 5): 4.3199905419650836e-08,
(0, 4, 5, 5): 1.4746049220222154e-10,
(0, 4, 4, 5): -3.340431558050326e-10,
(1, 4, 5, 5): -3.611277360054406e-11,
(1, 4, 4, 5): 2.0463905040308303e-10,
(2, 4, 5, 5): 2.8249217149025584e-08,
(2, 4, 4, 5): 4.815036480072542e-11,
(3, 4, 5, 5): -2.4045088422362252e-08,
(3, 4, 4, 5): -2.7084580200408048e-11}
[6]:
freqs = [
sqrt(k_new[(k, k)]) * units.au_in_cm1 for k in range(f)
] # a.u. (sqrt{omega^2} = omega)
freqs # in cm^{-1}
[6]:
[1178.0577664327557,
1250.2060681419946,
1508.7131015740217,
1810.5170221727915,
2890.7639406635694,
2951.821181689732]
[7]:
nprims = [N] * f # Number of primitive basis
dvr_prims = [
# HO(ngrid=nprim, omega=omega, units='cm-1')
HO(nprim, omega)
for nprim, omega in zip(nprims, freqs, strict=True)
]
dvr_prims_redundant = [
# HO(ngrid=nprim, omega=omega, units='cm-1')
HO(nprim+4, omega)
for nprim, omega in zip(nprims, freqs, strict=True)
]
ndof = len(dvr_prims) # Number of degree of freedom
dvr_prims[0].plot_fbr()
[8]:
dq2 = [prim.get_2nd_derivative_matrix_fbr() for prim in dvr_prims]
q_scale = 10
q1 = [prim.get_pos_rep_matrix() / q_scale for prim in dvr_prims_redundant]
q2 = [ints @ ints for ints in q1]
q3 = [ints @ ints @ ints for ints in q1]
q4 = [ints @ ints @ ints @ ints for ints in q1]
qn = [[ints[:N, :N] for ints in q] for q in [q1, q2, q3, q4]]
[9]:
kinetic_sop = pympo.SumOfProducts([])
for isite in range(f):
kinetic_sop -= (
1
/ 2
* pympo.OpSite(
r"\frac{\partial}{\partial q^2_{" + f"{isite}" + r"}}",
isite,
value=dq2[isite],
)
)
kinetic_sop.symbol
[9]:
$\displaystyle - 0.5 \frac{\partial}{\partial q^2_{0}} - 0.5 \frac{\partial}{\partial q^2_{1}} - 0.5 \frac{\partial}{\partial q^2_{2}} - 0.5 \frac{\partial}{\partial q^2_{3}} - 0.5 \frac{\partial}{\partial q^2_{4}} - 0.5 \frac{\partial}{\partial q^2_{5}}$
[10]:
potential_sop = pympo.SumOfProducts([])
coef_symbol = {}
for key, coef in k_new.items():
cnt_site = defaultdict(int)
for isite in key:
cnt_site[isite] += 1
op = 1
coef_sym = sympy.Symbol(f"k_{key}")
for isite, order in cnt_site.items():
coef /= factorial(order)
op *= pympo.OpSite(
f"q^{order}" + r"_{" + f"{isite}" + r"}",
isite,
value=qn[order - 1][isite],
)
op *= coef_sym
coef_symbol[coef_sym] = coef * (q_scale ** len(key))
potential_sop += op
potential_sop.symbol
[10]:
$\displaystyle k_{(0, 0)} q^2_{0} + k_{(0, 0, 0, 0)} q^4_{0} + k_{(0, 0, 0, 1)} q^3_{0} q^1_{1} + k_{(0, 0, 0, 4)} q^3_{0} q^1_{4} + k_{(0, 0, 0, 5)} q^3_{0} q^1_{5} + k_{(0, 0, 1, 1)} q^2_{0} q^2_{1} + k_{(0, 0, 1, 5)} q^2_{0} q^1_{1} q^1_{5} + k_{(0, 0, 2)} q^2_{0} q^1_{2} + k_{(0, 0, 2, 2)} q^2_{0} q^2_{2} + k_{(0, 0, 2, 3)} q^2_{0} q^1_{2} q^1_{3} + k_{(0, 0, 2, 4)} q^2_{0} q^1_{2} q^1_{4} + k_{(0, 0, 3)} q^2_{0} q^1_{3} + k_{(0, 0, 3, 3)} q^2_{0} q^2_{3} + k_{(0, 0, 3, 4)} q^2_{0} q^1_{3} q^1_{4} + k_{(0, 0, 4)} q^2_{0} q^1_{4} + k_{(0, 0, 4, 4)} q^2_{0} q^2_{4} + k_{(0, 0, 5, 5)} q^2_{0} q^2_{5} + k_{(0, 1, 1, 1)} q^1_{0} q^3_{1} + k_{(0, 1, 1, 2)} q^1_{0} q^2_{1} q^1_{2} + k_{(0, 1, 1, 3)} q^1_{0} q^2_{1} q^1_{3} + k_{(0, 1, 1, 4)} q^1_{0} q^2_{1} q^1_{4} + k_{(0, 1, 2, 2)} q^1_{0} q^1_{1} q^2_{2} + k_{(0, 1, 3, 3)} q^1_{0} q^1_{1} q^2_{3} + k_{(0, 1, 4)} q^1_{0} q^1_{1} q^1_{4} + k_{(0, 1, 4, 4)} q^1_{0} q^1_{1} q^2_{4} + k_{(0, 1, 5)} q^1_{0} q^1_{1} q^1_{5} + k_{(0, 1, 5, 5)} q^1_{0} q^1_{1} q^2_{5} + k_{(0, 2, 2, 4)} q^1_{0} q^2_{2} q^1_{4} + k_{(0, 2, 4, 4)} q^1_{0} q^1_{2} q^2_{4} + k_{(0, 2, 5)} q^1_{0} q^1_{2} q^1_{5} + k_{(0, 2, 5, 5)} q^1_{0} q^1_{2} q^2_{5} + k_{(0, 3, 3)} q^1_{0} q^2_{3} + k_{(0, 3, 3, 3)} q^1_{0} q^3_{3} + k_{(0, 3, 4)} q^1_{0} q^1_{3} q^1_{4} + k_{(0, 3, 4, 4)} q^1_{0} q^1_{3} q^2_{4} + k_{(0, 3, 5)} q^1_{0} q^1_{3} q^1_{5} + k_{(0, 3, 5, 5)} q^1_{0} q^1_{3} q^2_{5} + k_{(0, 4, 4)} q^1_{0} q^2_{4} + k_{(0, 4, 4, 4)} q^1_{0} q^3_{4} + k_{(0, 4, 4, 5)} q^1_{0} q^2_{4} q^1_{5} + k_{(0, 4, 5, 5)} q^1_{0} q^1_{4} q^2_{5} + k_{(0, 5, 5)} q^1_{0} q^2_{5} + k_{(0, 5, 5, 5)} q^1_{0} q^3_{5} + k_{(1, 1)} q^2_{1} + k_{(1, 1, 1, 1)} q^4_{1} + k_{(1, 1, 1, 5)} q^3_{1} q^1_{5} + k_{(1, 1, 2)} q^2_{1} q^1_{2} + k_{(1, 1, 2, 2)} q^2_{1} q^2_{2} + k_{(1, 1, 2, 3)} q^2_{1} q^1_{2} q^1_{3} + k_{(1, 1, 2, 4)} q^2_{1} q^1_{2} q^1_{4} + k_{(1, 1, 3)} q^2_{1} q^1_{3} + k_{(1, 1, 3, 3)} q^2_{1} q^2_{3} + k_{(1, 1, 3, 4)} q^2_{1} q^1_{3} q^1_{4} + k_{(1, 1, 4)} q^2_{1} q^1_{4} + k_{(1, 1, 4, 4)} q^2_{1} q^2_{4} + k_{(1, 1, 5, 5)} q^2_{1} q^2_{5} + k_{(1, 2, 2, 5)} q^1_{1} q^2_{2} q^1_{5} + k_{(1, 2, 4, 4)} q^1_{1} q^1_{2} q^2_{4} + k_{(1, 2, 5)} q^1_{1} q^1_{2} q^1_{5} + k_{(1, 2, 5, 5)} q^1_{1} q^1_{2} q^2_{5} + k_{(1, 3, 3, 5)} q^1_{1} q^2_{3} q^1_{5} + k_{(1, 3, 4, 4)} q^1_{1} q^1_{3} q^2_{4} + k_{(1, 3, 5)} q^1_{1} q^1_{3} q^1_{5} + k_{(1, 3, 5, 5)} q^1_{1} q^1_{3} q^2_{5} + k_{(1, 4, 4, 4)} q^1_{1} q^3_{4} + k_{(1, 4, 4, 5)} q^1_{1} q^2_{4} q^1_{5} + k_{(1, 4, 5)} q^1_{1} q^1_{4} q^1_{5} + k_{(1, 4, 5, 5)} q^1_{1} q^1_{4} q^2_{5} + k_{(1, 5, 5, 5)} q^1_{1} q^3_{5} + k_{(2, 2)} q^2_{2} + k_{(2, 2, 2)} q^3_{2} + k_{(2, 2, 2, 2)} q^4_{2} + k_{(2, 2, 2, 3)} q^3_{2} q^1_{3} + k_{(2, 2, 2, 4)} q^3_{2} q^1_{4} + k_{(2, 2, 3)} q^2_{2} q^1_{3} + k_{(2, 2, 3, 3)} q^2_{2} q^2_{3} + k_{(2, 2, 3, 4)} q^2_{2} q^1_{3} q^1_{4} + k_{(2, 2, 4)} q^2_{2} q^1_{4} + k_{(2, 2, 4, 4)} q^2_{2} q^2_{4} + k_{(2, 2, 5, 5)} q^2_{2} q^2_{5} + k_{(2, 3, 3)} q^1_{2} q^2_{3} + k_{(2, 3, 3, 3)} q^1_{2} q^3_{3} + k_{(2, 3, 3, 4)} q^1_{2} q^2_{3} q^1_{4} + k_{(2, 3, 4)} q^1_{2} q^1_{3} q^1_{4} + k_{(2, 3, 4, 4)} q^1_{2} q^1_{3} q^2_{4} + k_{(2, 3, 5, 5)} q^1_{2} q^1_{3} q^2_{5} + k_{(2, 4, 4)} q^1_{2} q^2_{4} + k_{(2, 4, 4, 4)} q^1_{2} q^3_{4} + k_{(2, 4, 4, 5)} q^1_{2} q^2_{4} q^1_{5} + k_{(2, 4, 5, 5)} q^1_{2} q^1_{4} q^2_{5} + k_{(2, 5, 5)} q^1_{2} q^2_{5} + k_{(3, 3)} q^2_{3} + k_{(3, 3, 3)} q^3_{3} + k_{(3, 3, 3, 3)} q^4_{3} + k_{(3, 3, 3, 4)} q^3_{3} q^1_{4} + k_{(3, 3, 4)} q^2_{3} q^1_{4} + k_{(3, 3, 4, 4)} q^2_{3} q^2_{4} + k_{(3, 3, 5, 5)} q^2_{3} q^2_{5} + k_{(3, 4, 4)} q^1_{3} q^2_{4} + k_{(3, 4, 4, 4)} q^1_{3} q^3_{4} + k_{(3, 4, 4, 5)} q^1_{3} q^2_{4} q^1_{5} + k_{(3, 4, 5, 5)} q^1_{3} q^1_{4} q^2_{5} + k_{(3, 5, 5)} q^1_{3} q^2_{5} + k_{(4, 4)} q^2_{4} + k_{(4, 4, 4)} q^3_{4} + k_{(4, 4, 4, 4)} q^4_{4} + k_{(4, 4, 4, 5)} q^3_{4} q^1_{5} + k_{(4, 4, 5)} q^2_{4} q^1_{5} + k_{(4, 4, 5, 5)} q^2_{4} q^2_{5} + k_{(4, 5, 5)} q^1_{4} q^2_{5} + k_{(4, 5, 5, 5)} q^1_{4} q^3_{5} + k_{(5, 5)} q^2_{5} + k_{(5, 5, 5)} q^3_{5} + k_{(5, 5, 5, 5)} q^4_{5}$
[11]:
coef_symbol
[11]:
{k_(0, 0): 0.0014405705141986437,
k_(1, 1): 0.001622424712773225,
k_(2, 2): 0.002362733132008209,
k_(3, 3): 0.003402566113601328,
k_(4, 4): 0.008674143553126032,
k_(5, 5): 0.009044434907689655,
k_(2, 2, 2): -1.959429052471348e-05,
k_(3, 3, 3): 0.00032612180371678867,
k_(4, 4, 4): 0.0015550157447563595,
k_(5, 5, 5): -1.6489184375988388e-06,
k_(0, 0, 2): -6.193080677319275e-05,
k_(1, 1, 2): 0.00011666633309140877,
k_(0, 3, 3): 3.8546145294518314e-07,
k_(0, 0, 3): 3.4498800038593886e-05,
k_(1, 1, 3): 2.3577392205147037e-05,
k_(2, 3, 3): 0.0001705024493527527,
k_(2, 2, 3): 0.0001429419554671721,
k_(0, 4, 4): 1.477602236289869e-06,
k_(0, 0, 4): -0.00045548695023022486,
k_(1, 1, 4): -0.0003440885903290668,
k_(2, 4, 4): -3.719703020921018e-05,
k_(2, 2, 4): -0.0001678042191821364,
k_(3, 4, 4): -4.65765922308763e-05,
k_(3, 3, 4): -0.00014564018563778837,
k_(0, 5, 5): 1.6060893872715966e-06,
k_(2, 5, 5): 0.0003217960696337371,
k_(3, 5, 5): -0.0003466583333487014,
k_(4, 5, 5): 0.005020121476007083,
k_(4, 4, 5): 1.7345765382533244e-06,
k_(0, 1, 4): -5.13948603926911e-07,
k_(0, 3, 4): -6.424357549086388e-07,
k_(2, 3, 4): 0.00024720927848884415,
k_(0, 1, 5): -1.798820113744188e-06,
k_(0, 2, 5): 1.027897207853822e-06,
k_(1, 2, 5): 0.0006678762108030207,
k_(0, 3, 5): -6.424357549086387e-07,
k_(1, 3, 5): -0.0004840110977481683,
k_(1, 4, 5): -2.980901902776083e-05,
k_(0, 0, 0, 0): 7.47459178586261e-06,
k_(1, 1, 1, 1): 6.056413072591244e-06,
k_(2, 2, 2, 2): 2.242001361033777e-06,
k_(3, 3, 3, 3): 2.1633808435075923e-05,
k_(4, 4, 4, 4): 0.00017623158908640504,
k_(5, 5, 5, 5): 0.00021124217074268252,
k_(0, 1, 1, 1): 1.4545422700219133e-07,
k_(0, 0, 1, 1): 1.0517845311158455e-05,
k_(0, 0, 0, 1): 2.507831500037782e-07,
k_(0, 0, 2, 2): 5.710332325586029e-06,
k_(1, 1, 2, 2): 2.0787415303813172e-05,
k_(0, 3, 3, 3): 4.0125304000604505e-08,
k_(0, 0, 3, 3): -5.793090765087277e-06,
k_(1, 1, 3, 3): -2.257048350034004e-07,
k_(2, 3, 3, 3): 1.6546672237249285e-05,
k_(2, 2, 3, 3): 1.1533517068673759e-05,
k_(2, 2, 2, 3): -5.3918377250812315e-06,
k_(0, 4, 4, 4): 2.6583013900400484e-07,
k_(0, 0, 4, 4): -0.00022195061124784382,
k_(0, 0, 0, 4): -1.2037591200181357e-07,
k_(1, 4, 4, 4): -4.5140967000680075e-08,
k_(1, 1, 4, 4): -0.00018262530549575137,
k_(2, 4, 4, 4): -8.195593342123471e-06,
k_(2, 2, 4, 4): -0.00015595451749284955,
k_(2, 2, 2, 4): -9.228819920139039e-06,
k_(3, 4, 4, 4): -1.1531009237173723e-05,
k_(3, 3, 4, 4): -3.420932949201539e-05,
k_(3, 3, 3, 4): -7.543557152113648e-06,
k_(0, 5, 5, 5): -6.520361900098232e-07,
k_(0, 0, 5, 5): -0.00026749032345702993,
k_(0, 0, 0, 5): -1.454542270021914e-07,
k_(1, 5, 5, 5): -5.434470860581874e-05,
k_(1, 1, 5, 5): -0.00019647605887045997,
k_(1, 1, 1, 5): -7.6137764341147074e-06,
k_(2, 2, 5, 5): -0.00021420893540722717,
k_(3, 3, 5, 5): -6.162494345042841e-05,
k_(4, 5, 5, 5): -4.71472322007103e-07,
k_(4, 4, 5, 5): 0.0011511397994843427,
k_(4, 4, 4, 5): 5.366759410080853e-07,
k_(0, 1, 2, 2): 4.6645665900702737e-07,
k_(0, 1, 1, 2): -1.0532892300158684e-07,
k_(0, 1, 3, 3): 1.2037591200181354e-07,
k_(0, 1, 1, 3): 6.018795600090678e-08,
k_(0, 0, 2, 3): -1.1435711640172285e-05,
k_(1, 1, 2, 3): -2.0418764073307625e-05,
k_(0, 1, 4, 4): -4.2131569200634735e-07,
k_(0, 1, 1, 4): -1.3542290100204023e-07,
k_(0, 2, 4, 4): -6.018795600090677e-08,
k_(0, 2, 2, 4): -1.0532892300158684e-07,
k_(0, 0, 2, 4): -1.4746049220222159e-05,
k_(1, 2, 4, 4): 2.2570483500340037e-07,
k_(1, 1, 2, 4): 3.251654322948988e-05,
k_(0, 3, 4, 4): -9.028193400136014e-08,
k_(0, 0, 3, 4): 2.64977476293992e-05,
k_(1, 3, 4, 4): -1.0532892300158684e-07,
k_(1, 1, 3, 4): 1.8688360338281553e-05,
k_(2, 3, 4, 4): 0.00015853507610638843,
k_(2, 3, 3, 4): -2.679868740940373e-05,
k_(2, 2, 3, 4): 3.736167368756288e-05,
k_(0, 1, 5, 5): -5.416916040081609e-07,
k_(0, 0, 1, 5): 1.8658266360281095e-06,
k_(0, 2, 5, 5): -6.018795600090677e-08,
k_(1, 2, 5, 5): -2.4075182400362713e-07,
k_(1, 2, 2, 5): 2.554978732238492e-05,
k_(0, 3, 5, 5): -1.2037591200181354e-07,
k_(1, 3, 5, 5): 1.2037591200181354e-07,
k_(1, 3, 3, 5): 1.3000598496195863e-05,
k_(2, 3, 5, 5): 0.00021599952709825418,
k_(0, 4, 5, 5): 7.373024610111077e-07,
k_(0, 4, 4, 5): -1.670215779025163e-06,
k_(1, 4, 5, 5): -1.805638680027203e-07,
k_(1, 4, 4, 5): 1.023195252015415e-06,
k_(2, 4, 5, 5): 0.0001412460857451279,
k_(2, 4, 4, 5): 2.407518240036271e-07,
k_(3, 4, 5, 5): -0.00012022544211181125,
k_(3, 4, 4, 5): -1.3542290100204023e-07}
[12]:
hamiltonian_sop = pympo.SumOfProducts([])
hamiltonian_sop += kinetic_sop
hamiltonian_sop += potential_sop
hamiltonian_sop = hamiltonian_sop.simplify()
hamiltonian_sop.symbol
[12]:
$\displaystyle k_{(0, 0)} q^2_{0} + k_{(0, 0, 0, 0)} q^4_{0} + k_{(0, 0, 0, 1)} q^3_{0} q^1_{1} + k_{(0, 0, 0, 4)} q^3_{0} q^1_{4} + k_{(0, 0, 0, 5)} q^3_{0} q^1_{5} + k_{(0, 0, 1, 1)} q^2_{0} q^2_{1} + k_{(0, 0, 1, 5)} q^2_{0} q^1_{1} q^1_{5} + k_{(0, 0, 2)} q^2_{0} q^1_{2} + k_{(0, 0, 2, 2)} q^2_{0} q^2_{2} + k_{(0, 0, 2, 3)} q^2_{0} q^1_{2} q^1_{3} + k_{(0, 0, 2, 4)} q^2_{0} q^1_{2} q^1_{4} + k_{(0, 0, 3)} q^2_{0} q^1_{3} + k_{(0, 0, 3, 3)} q^2_{0} q^2_{3} + k_{(0, 0, 3, 4)} q^2_{0} q^1_{3} q^1_{4} + k_{(0, 0, 4)} q^2_{0} q^1_{4} + k_{(0, 0, 4, 4)} q^2_{0} q^2_{4} + k_{(0, 0, 5, 5)} q^2_{0} q^2_{5} + k_{(0, 1, 1, 1)} q^1_{0} q^3_{1} + k_{(0, 1, 1, 2)} q^1_{0} q^2_{1} q^1_{2} + k_{(0, 1, 1, 3)} q^1_{0} q^2_{1} q^1_{3} + k_{(0, 1, 1, 4)} q^1_{0} q^2_{1} q^1_{4} + k_{(0, 1, 2, 2)} q^1_{0} q^1_{1} q^2_{2} + k_{(0, 1, 3, 3)} q^1_{0} q^1_{1} q^2_{3} + k_{(0, 1, 4)} q^1_{0} q^1_{1} q^1_{4} + k_{(0, 1, 4, 4)} q^1_{0} q^1_{1} q^2_{4} + k_{(0, 1, 5)} q^1_{0} q^1_{1} q^1_{5} + k_{(0, 1, 5, 5)} q^1_{0} q^1_{1} q^2_{5} + k_{(0, 2, 2, 4)} q^1_{0} q^2_{2} q^1_{4} + k_{(0, 2, 4, 4)} q^1_{0} q^1_{2} q^2_{4} + k_{(0, 2, 5)} q^1_{0} q^1_{2} q^1_{5} + k_{(0, 2, 5, 5)} q^1_{0} q^1_{2} q^2_{5} + k_{(0, 3, 3)} q^1_{0} q^2_{3} + k_{(0, 3, 3, 3)} q^1_{0} q^3_{3} + k_{(0, 3, 4)} q^1_{0} q^1_{3} q^1_{4} + k_{(0, 3, 4, 4)} q^1_{0} q^1_{3} q^2_{4} + k_{(0, 3, 5)} q^1_{0} q^1_{3} q^1_{5} + k_{(0, 3, 5, 5)} q^1_{0} q^1_{3} q^2_{5} + k_{(0, 4, 4)} q^1_{0} q^2_{4} + k_{(0, 4, 4, 4)} q^1_{0} q^3_{4} + k_{(0, 4, 4, 5)} q^1_{0} q^2_{4} q^1_{5} + k_{(0, 4, 5, 5)} q^1_{0} q^1_{4} q^2_{5} + k_{(0, 5, 5)} q^1_{0} q^2_{5} + k_{(0, 5, 5, 5)} q^1_{0} q^3_{5} + k_{(1, 1)} q^2_{1} + k_{(1, 1, 1, 1)} q^4_{1} + k_{(1, 1, 1, 5)} q^3_{1} q^1_{5} + k_{(1, 1, 2)} q^2_{1} q^1_{2} + k_{(1, 1, 2, 2)} q^2_{1} q^2_{2} + k_{(1, 1, 2, 3)} q^2_{1} q^1_{2} q^1_{3} + k_{(1, 1, 2, 4)} q^2_{1} q^1_{2} q^1_{4} + k_{(1, 1, 3)} q^2_{1} q^1_{3} + k_{(1, 1, 3, 3)} q^2_{1} q^2_{3} + k_{(1, 1, 3, 4)} q^2_{1} q^1_{3} q^1_{4} + k_{(1, 1, 4)} q^2_{1} q^1_{4} + k_{(1, 1, 4, 4)} q^2_{1} q^2_{4} + k_{(1, 1, 5, 5)} q^2_{1} q^2_{5} + k_{(1, 2, 2, 5)} q^1_{1} q^2_{2} q^1_{5} + k_{(1, 2, 4, 4)} q^1_{1} q^1_{2} q^2_{4} + k_{(1, 2, 5)} q^1_{1} q^1_{2} q^1_{5} + k_{(1, 2, 5, 5)} q^1_{1} q^1_{2} q^2_{5} + k_{(1, 3, 3, 5)} q^1_{1} q^2_{3} q^1_{5} + k_{(1, 3, 4, 4)} q^1_{1} q^1_{3} q^2_{4} + k_{(1, 3, 5)} q^1_{1} q^1_{3} q^1_{5} + k_{(1, 3, 5, 5)} q^1_{1} q^1_{3} q^2_{5} + k_{(1, 4, 4, 4)} q^1_{1} q^3_{4} + k_{(1, 4, 4, 5)} q^1_{1} q^2_{4} q^1_{5} + k_{(1, 4, 5)} q^1_{1} q^1_{4} q^1_{5} + k_{(1, 4, 5, 5)} q^1_{1} q^1_{4} q^2_{5} + k_{(1, 5, 5, 5)} q^1_{1} q^3_{5} + k_{(2, 2)} q^2_{2} + k_{(2, 2, 2)} q^3_{2} + k_{(2, 2, 2, 2)} q^4_{2} + k_{(2, 2, 2, 3)} q^3_{2} q^1_{3} + k_{(2, 2, 2, 4)} q^3_{2} q^1_{4} + k_{(2, 2, 3)} q^2_{2} q^1_{3} + k_{(2, 2, 3, 3)} q^2_{2} q^2_{3} + k_{(2, 2, 3, 4)} q^2_{2} q^1_{3} q^1_{4} + k_{(2, 2, 4)} q^2_{2} q^1_{4} + k_{(2, 2, 4, 4)} q^2_{2} q^2_{4} + k_{(2, 2, 5, 5)} q^2_{2} q^2_{5} + k_{(2, 3, 3)} q^1_{2} q^2_{3} + k_{(2, 3, 3, 3)} q^1_{2} q^3_{3} + k_{(2, 3, 3, 4)} q^1_{2} q^2_{3} q^1_{4} + k_{(2, 3, 4)} q^1_{2} q^1_{3} q^1_{4} + k_{(2, 3, 4, 4)} q^1_{2} q^1_{3} q^2_{4} + k_{(2, 3, 5, 5)} q^1_{2} q^1_{3} q^2_{5} + k_{(2, 4, 4)} q^1_{2} q^2_{4} + k_{(2, 4, 4, 4)} q^1_{2} q^3_{4} + k_{(2, 4, 4, 5)} q^1_{2} q^2_{4} q^1_{5} + k_{(2, 4, 5, 5)} q^1_{2} q^1_{4} q^2_{5} + k_{(2, 5, 5)} q^1_{2} q^2_{5} + k_{(3, 3)} q^2_{3} + k_{(3, 3, 3)} q^3_{3} + k_{(3, 3, 3, 3)} q^4_{3} + k_{(3, 3, 3, 4)} q^3_{3} q^1_{4} + k_{(3, 3, 4)} q^2_{3} q^1_{4} + k_{(3, 3, 4, 4)} q^2_{3} q^2_{4} + k_{(3, 3, 5, 5)} q^2_{3} q^2_{5} + k_{(3, 4, 4)} q^1_{3} q^2_{4} + k_{(3, 4, 4, 4)} q^1_{3} q^3_{4} + k_{(3, 4, 4, 5)} q^1_{3} q^2_{4} q^1_{5} + k_{(3, 4, 5, 5)} q^1_{3} q^1_{4} q^2_{5} + k_{(3, 5, 5)} q^1_{3} q^2_{5} + k_{(4, 4)} q^2_{4} + k_{(4, 4, 4)} q^3_{4} + k_{(4, 4, 4, 4)} q^4_{4} + k_{(4, 4, 4, 5)} q^3_{4} q^1_{5} + k_{(4, 4, 5)} q^2_{4} q^1_{5} + k_{(4, 4, 5, 5)} q^2_{4} q^2_{5} + k_{(4, 5, 5)} q^1_{4} q^2_{5} + k_{(4, 5, 5, 5)} q^1_{4} q^3_{5} + k_{(5, 5)} q^2_{5} + k_{(5, 5, 5)} q^3_{5} + k_{(5, 5, 5, 5)} q^4_{5} - 0.5 \frac{\partial}{\partial q^2_{0}} - 0.5 \frac{\partial}{\partial q^2_{1}} - 0.5 \frac{\partial}{\partial q^2_{2}} - 0.5 \frac{\partial}{\partial q^2_{3}} - 0.5 \frac{\partial}{\partial q^2_{4}} - 0.5 \frac{\partial}{\partial q^2_{5}}$
[13]:
hamiltonian_mpo = hamiltonian_sop.to_mpo(subs=coef_symbol)
[14]:
operators = {"hamiltonian": hamiltonian_mpo}
[15]:
model = Model(dvr_prims, operators=operators, bond_dim=m)
model.init_HartreeProduct = [
[
[1.0] + [0.0] * (N-1) for _ in range(f)
]
]
[16]:
jobname = "ch2o-mpo"
simulator = Simulator(jobname=jobname, model=model, backend=backend)
import time
start = time.time()
ener, wf = simulator.relax(savefile_ext="_gs", maxstep=10, stepsize=0.1)#.(maxstep=100, stepsize=Δt)
end = time.time()
print(end-start)
ener * units.au_in_eV
19:42:26 | INFO | Log file is ./ch2o-mpo_relax/main.log
19:42:26 | INFO | Wave function is saved in wf_ch2o-mpo_gs.pkl
19:42:26 | INFO | Start initial step 0.000 [fs]
10%|████████▍ | 1/10 [00:01<00:14, 1.62s/it]19:42:28 | INFO | End 0 step; propagated 0.100 [fs]; AVG Krylov iteration: 18.33
80%|███████████████████████████████████████████████████████████████████▏ | 8/10 [00:02<00:00, 6.03it/s]19:42:29 | INFO | Saved wavefunction 0.900 [fs]
100%|███████████████████████████████████████████████████████████████████████████████████| 10/10 [00:02<00:00, 4.06it/s]
19:42:29 | INFO | End 9 step; propagated 0.900 [fs]; AVG Krylov iteration: 3.17
19:42:29 | INFO | End simulation and save wavefunction
19:42:29 | INFO | Wave function is saved in wf_ch2o-mpo_gs.pkl
2.6824569702148438
[16]:
0.7088296546416561
[17]:
!tail ch2o-mpo_relax/main.log
2025-12-04 19:42:28 | DEBUG | pytdscf.properties:_export_properties:410 - | pop 1.0000 | ene[eV]: 0.7088297 | time[fs]: 0.400 | elapsed[sec]: 1.32
2025-12-04 19:42:28 | DEBUG | pytdscf.properties:_export_properties:410 - | pop 1.0000 | ene[eV]: 0.7088297 | time[fs]: 0.500 | elapsed[sec]: 1.39
2025-12-04 19:42:28 | DEBUG | pytdscf.properties:_export_properties:410 - | pop 1.0000 | ene[eV]: 0.7088297 | time[fs]: 0.600 | elapsed[sec]: 1.45
2025-12-04 19:42:28 | DEBUG | pytdscf.properties:_export_properties:410 - | pop 1.0000 | ene[eV]: 0.7088297 | time[fs]: 0.700 | elapsed[sec]: 1.50
2025-12-04 19:42:29 | DEBUG | pytdscf.properties:_export_properties:410 - | pop 1.0000 | ene[eV]: 0.7088297 | time[fs]: 0.800 | elapsed[sec]: 1.56
2025-12-04 19:42:29 | INFO | pytdscf.simulator_cls:_execute:416 - Saved wavefunction 0.900 [fs]
2025-12-04 19:42:29 | DEBUG | pytdscf.properties:_export_properties:410 - | pop 1.0000 | ene[eV]: 0.7088297 | time[fs]: 0.900 | elapsed[sec]: 1.61
2025-12-04 19:42:29 | INFO | pytdscf.simulator_cls:_execute:465 - End 9 step; propagated 0.900 [fs]; AVG Krylov iteration: 3.17
2025-12-04 19:42:29 | INFO | pytdscf.simulator_cls:_execute:466 - End simulation and save wavefunction
2025-12-04 19:42:29 | INFO | pytdscf.simulator_cls:save_wavefunction:590 - Wave function is saved in wf_ch2o-mpo_gs.pkl
/Users/hinom/.local/share/uv/python/cpython-3.12.2-macos-aarch64-none/lib/python3.12/pty.py:95: RuntimeWarning: os.fork() was called. os.fork() is incompatible with multithreaded code, and JAX is multithreaded, so this will likely lead to a deadlock.
pid, fd = os.forkpty()
Confirm with SOP model
[18]:
from discvar import PrimBas_HO
prim_info = [
PrimBas_HO(0.0, omega, nprim)
for nprim, omega in zip(nprims, freqs, strict=True)
]
ndof = len(prim_info) # Number of degree of freedom, H2O has 3 DOFs
from pytdscf.hamiltonian_cls import read_potential_nMR
hamiltonian = read_potential_nMR(k_orig, cut_off=cutoff)
operators = {"hamiltonian": hamiltonian}
model = Model(prim_info, operators, bond_dim=m)
model.init_HartreeProduct = [
[
[1.0] + [0.0] * (N-1) for _ in range(f)
]
]
jobname = "ch2o-sop"
simulator = Simulator(jobname, model, ci_type="MPS", backend="Numpy")
start = time.time()
ener, wf = simulator.relax(savefile_ext="_gs", maxstep=10, stepsize=0.1)
end = time.time()
print(end-start)
ener * units.au_in_eV
19:42:29 | INFO | Construct anharmonic polynomial operator
19:42:29 | INFO | Log file is ./ch2o-sop_relax/main.log
19:42:29 | WARNING | rebind const doDVR
19:42:29 | WARNING | rebind const use_jax
19:42:29 | WARNING | Employing sum of product Hamiltonian is deprecated.
19:42:29 | WARNING | rebind const use_mpo
19:42:29 | WARNING | Sum-of-Products (SOP) based calculation will be executed, but it is inefficient and will be deprecated, use Conventional MPO or Grid-MPO based calculation instead. See example such as https://qclovers.github.io/PyTDSCF/notebook/TD_reduced_density_exciton.html and https://kenhino.github.io/PyMPO/example/pytdscf-taylor.html
19:42:29 | INFO | Wave function is saved in wf_ch2o-sop_gs.pkl
19:42:29 | INFO | Start initial step 0.000 [fs]
10%|████████▍ | 1/10 [00:04<00:36, 4.10s/it]19:42:33 | INFO | End 0 step; propagated 0.100 [fs]; AVG Krylov iteration: 18.33
90%|███████████████████████████████████████████████████████████████████████████▌ | 9/10 [00:11<00:00, 1.27it/s]19:42:40 | INFO | Saved wavefunction 0.900 [fs]
100%|███████████████████████████████████████████████████████████████████████████████████| 10/10 [00:11<00:00, 1.19s/it]
19:42:41 | INFO | End 9 step; propagated 0.900 [fs]; AVG Krylov iteration: 3.17
19:42:41 | INFO | End simulation and save wavefunction
19:42:41 | INFO | Wave function is saved in wf_ch2o-sop_gs.pkl
11.949464082717896
[18]:
0.7088296546416549
[19]:
!tail ch2o-sop_relax/main.log
2025-12-04 19:42:37 | DEBUG | pytdscf.properties:_export_properties:410 - | pop 1.0000 | ene[eV]: 0.7088297 | time[fs]: 0.400 | elapsed[sec]: 9.18
2025-12-04 19:42:37 | DEBUG | pytdscf.properties:_export_properties:410 - | pop 1.0000 | ene[eV]: 0.7088297 | time[fs]: 0.500 | elapsed[sec]: 9.96
2025-12-04 19:42:38 | DEBUG | pytdscf.properties:_export_properties:410 - | pop 1.0000 | ene[eV]: 0.7088297 | time[fs]: 0.600 | elapsed[sec]: 10.69
2025-12-04 19:42:39 | DEBUG | pytdscf.properties:_export_properties:410 - | pop 1.0000 | ene[eV]: 0.7088297 | time[fs]: 0.700 | elapsed[sec]: 11.38
2025-12-04 19:42:40 | DEBUG | pytdscf.properties:_export_properties:410 - | pop 1.0000 | ene[eV]: 0.7088297 | time[fs]: 0.800 | elapsed[sec]: 12.01
2025-12-04 19:42:40 | INFO | pytdscf.simulator_cls:_execute:416 - Saved wavefunction 0.900 [fs]
2025-12-04 19:42:40 | DEBUG | pytdscf.properties:_export_properties:410 - | pop 1.0000 | ene[eV]: 0.7088297 | time[fs]: 0.900 | elapsed[sec]: 12.57
2025-12-04 19:42:41 | INFO | pytdscf.simulator_cls:_execute:465 - End 9 step; propagated 0.900 [fs]; AVG Krylov iteration: 3.17
2025-12-04 19:42:41 | INFO | pytdscf.simulator_cls:_execute:466 - End simulation and save wavefunction
2025-12-04 19:42:41 | INFO | pytdscf.simulator_cls:save_wavefunction:590 - Wave function is saved in wf_ch2o-sop_gs.pkl