Example 3: Propagation of H2O vibration under polynomial PES

run type

wavefunction

backend

Basis

steps

propagation

MPS-SM

Numpy

HO-DVR

5000

1. Import modules

  • Required in any calculations

[1]:

import platform
import sys

from pytdscf import BasInfo, Model, Simulator, __version__

print(
    f"pytdscf version: {__version__}, Python version: {sys.version}, platform: {platform.platform()}"
)

pytdscf version: 1.2.0, Python version: 3.12.2 (main, Feb 25 2024, 04:38:01) [Clang 17.0.6 ], platform: Linux-5.15.167.4-microsoft-standard-WSL2-x86_64-with-glibc2.39

2. Set DVR primitive basis

[2]:

from math import sqrt

import numpy as np
from discvar import HarmonicOscillator as HO
from pympo.utils import import_npz

from pytdscf.dvr_operator_cls import TensorOperator
from pytdscf.hamiltonian_cls import TensorHamiltonian
from pytdscf.potentials.h2o_potential import k_orig

backend = "numpy"
ndim = 3
freqs = [sqrt(k_orig[(k, k)]) for k in range(1, ndim + 1)]  # a.u.

nprims = [9] * ndim  # Number of primitive basis

basis = [
    HO(nprim, omega, units="a.u.")
    for nprim, omega in zip(nprims, freqs, strict=True)
]
basinfo = BasInfo([basis])  # Set basis information object

MPS-MCTDH wavefunction

\[\begin{split}|\Psi_{\rm{MPS-MCTDH}}\rangle = \sum_{\mathbf \{j\}}\sum_{\mathbf \{\tau\}} a\substack{j_1 \\ 1\tau_1}a\substack{j_2 \\ \tau_1\tau_2} \cdots a\substack{j_f \\ \tau_{f-1}1} |\varphi_{j_1}^{(1)}(q_1)\rangle|\varphi_{j_2}^{(2)}(q_2)\rangle \cdots|\varphi_{j_f}^{(f)}(q_f)\rangle\end{split}\]

where SPF is

\[\varphi_{j_p}^{(p)}(q_p) = \sum_{i_p=1}^{n_p} c_{i_p}^{j_p}\chi_{i_p}^{(p)}(q_p) \; (j_p = 1,2,\ldots, N_p)\]

Here, select \(\{\chi_{i_p}^{(p)}(q_p)\}\) as Harmonic Oscillator eigenfunction. See detail in documenation. Here one define \(n_p\) = 9, \(N_p\) = 9. (Standard Method)

NOTE

  • In MPS, \(n = N\) (SM) is usually better than \(n < M\) (MCTDH). Only when using a laptop, MCTDH may be better. (RAM required in MCTDH is smaller than SM.)

3. Set Hamiltonian (Taylor expansion polynomial PES in MPO format)

We have already constructed MPO in the previous relaxation step.

[3]:

pot_mpo = import_npz("h2o_pot_mpo.npz")
kin_mpo = import_npz("h2o_kin_mpo.npz")

Setup Hamiltonian

[4]:

potential = [
    [
        {
            (tuple((i,) for i in range(0, ndim))): TensorOperator(mpo=pot_mpo),
            (tuple((i, i) for i in range(0, ndim))): TensorOperator(
                mpo=kin_mpo
            ),
        }
    ]
]

H = TensorHamiltonian(
    ndof=len(basis), potential=potential, kinetic=None, backend=backend
)

operators = {"hamiltonian": H}

4. Set wavefunction (MPS) and All Model

  • m_aux_max is MPS bond dimension (maximum of auxiliary index \(\tau_p\))

[5]:

model = Model(basinfo, operators)
model.m_aux_max = 9

5. Execute Calculation

  • time step width is defined by stepsize=0.2 fs

In this calculation, one runs

  • Real-time propagation

  • Restart from \(\hat{\mu}|\Psi_{\rm GS}\rangle\) wavefunction. (restart file suffix is _dipole)

F.Y.I., See also documentation

NOTE

  • Runtype cannnot rebind. If you change runtype, you should restart the kernel.

  • JAX is better when simulating more large systems. (f>6, m>10)

  • If AVG Krylov iteration in the log file is much larger than 5, you should set smaller timestep.

[6]:

jobname = "h2o_polynomial"
simulator = Simulator(jobname, model, backend="numpy")
simulator.propagate(
    maxstep=1000,  # 200 fs
    stepsize=0.2,
    restart=True,
    savefile_ext="_prop",
    loadfile_ext="_dipole",
)

13:54:31 | INFO | Log file is ./h2o_polynomial_prop/main.log
13:54:31 | INFO | Wave function is loaded from wf_h2o_polynomial_dipole.pkl
13:54:31 | INFO | Wave function is saved in wf_h2o_polynomial_prop.pkl
13:54:31 | INFO | Start initial step    0.000 [fs]

13:54:31 | INFO | End     0 step; propagated    0.200 [fs]; AVG Krylov iteration: 7.67
13:54:33 | INFO | End   100 step; propagated   20.200 [fs]; AVG Krylov iteration: 8.00
13:54:36 | INFO | End   200 step; propagated   40.200 [fs]; AVG Krylov iteration: 8.00
13:54:38 | INFO | End   300 step; propagated   60.200 [fs]; AVG Krylov iteration: 8.00
13:54:40 | INFO | End   400 step; propagated   80.200 [fs]; AVG Krylov iteration: 8.00
13:54:42 | INFO | End   500 step; propagated  100.200 [fs]; AVG Krylov iteration: 8.33
13:54:45 | INFO | End   600 step; propagated  120.200 [fs]; AVG Krylov iteration: 8.33
13:54:47 | INFO | End   700 step; propagated  140.200 [fs]; AVG Krylov iteration: 8.33
13:54:49 | INFO | End   800 step; propagated  160.200 [fs]; AVG Krylov iteration: 8.33
13:54:51 | INFO | End   900 step; propagated  180.200 [fs]; AVG Krylov iteration: 8.33
13:54:54 | INFO | Saved wavefunction  199.800 [fs]
13:54:54 | INFO | End   999 step; propagated  199.800 [fs]; AVG Krylov iteration: 8.33
13:54:54 | INFO | End simulation and save wavefunction
13:54:54 | INFO | Wave function is saved in wf_h2o_polynomial_prop.pkl

[6]:

(0.020896110045924655, <pytdscf.wavefunction.WFunc at 0x7f10c78f6840>)

6. Check Log file

See h2o_polynomial_prop/main.log, which is defined as jobname.

[7]:

!tail h2o_polynomial_prop/main.log

2025-09-23 13:54:53 | DEBUG | pytdscf.properties:_export_properties:410 - | autocorr: -0.9275+0.3656i| pop 1.0000 | ene[eV]:  0.5686121 | time[fs]:  198.800 | elapsed[sec]:    19.69
2025-09-23 13:54:53 | DEBUG | pytdscf.properties:_export_properties:410 - | autocorr: -0.7493+0.6580i| pop 1.0000 | ene[eV]:  0.5686121 | time[fs]:  199.000 | elapsed[sec]:    19.70
2025-09-23 13:54:54 | DEBUG | pytdscf.properties:_export_properties:410 - | autocorr: -0.4826+0.8729i| pop 1.0000 | ene[eV]:  0.5686121 | time[fs]:  199.200 | elapsed[sec]:    19.72
2025-09-23 13:54:54 | DEBUG | pytdscf.properties:_export_properties:410 - | autocorr: -0.1589+0.9849i| pop 1.0000 | ene[eV]:  0.5686121 | time[fs]:  199.400 | elapsed[sec]:    19.74
2025-09-23 13:54:54 | DEBUG | pytdscf.properties:_export_properties:410 - | autocorr:  0.1837+0.9808i| pop 1.0000 | ene[eV]:  0.5686121 | time[fs]:  199.600 | elapsed[sec]:    19.76
2025-09-23 13:54:54 | INFO | pytdscf.simulator_cls:_execute:415 - Saved wavefunction  199.800 [fs]
2025-09-23 13:54:54 | DEBUG | pytdscf.properties:_export_properties:410 - | autocorr:  0.5048+0.8609i| pop 1.0000 | ene[eV]:  0.5686121 | time[fs]:  199.800 | elapsed[sec]:    19.78
2025-09-23 13:54:54 | INFO | pytdscf.simulator_cls:_execute:464 - End   999 step; propagated  199.800 [fs]; AVG Krylov iteration: 8.33
2025-09-23 13:54:54 | INFO | pytdscf.simulator_cls:_execute:465 - End simulation and save wavefunction
2025-09-23 13:54:54 | INFO | pytdscf.simulator_cls:save_wavefunction:589 - Wave function is saved in wf_h2o_polynomial_prop.pkl