Numerical Propagation Algorithms¶

In the previous section, we derived the Equations of Motion (EOMs) for the expansion coefficients and the Single Particle Functions (SPFs). These equations completely describe the time evolution of the wavefunction, but they cannot be solved analytically for realistic molecular systems.

Instead, MCTDH integrates these coupled differential equations numerically. During propagation, the program repeatedly evaluates the Hamiltonian, updates the coefficients and SPFs, checks numerical stability, and advances the wavepacket to the next time step until the requested propagation time is reached.

Why Numerical Integration is Needed¶

The coupled MCTDH equations have the general form

\[ i\hbar\frac{d\mathbf{y}}{dt} = \mathbf{F}(\mathbf{y},t), \]

where

\[ \mathbf y= \left( A_J, \varphi_j^{(\kappa)} \right) \]

contains every time-dependent quantity in the calculation.

Since

the Hamiltonian depends on the SPFs,
the SPFs depend on the coefficients,
and both depend on time,

there is no closed-form solution.

Instead, the equations are integrated step by step.

Propagation Cycle¶

At every time step, MCTDH performs the following operations:

Current wavefunction
        │
        ▼
Evaluate Hamiltonian
        │
        ▼
Compute mean-field operators
        │
        ▼
Integrate coefficient equations
        │
        ▼
Integrate SPF equations
        │
        ▼
Orthogonalize SPFs
        │
        ▼
Normalize wavefunction
        │
        ▼
Compute observables
        │
        ▼
Advance to next time step

This procedure is repeated until the final propagation time is reached.

Variable Time-Step Integration¶

Unlike simple fixed-step numerical methods,

MCTDH automatically adjusts the integration step during propagation.

If the wavefunction changes slowly,

larger time steps are used.

When rapid changes occur,

the step size is reduced automatically.

This adaptive strategy

improves numerical stability,
preserves accuracy,
reduces unnecessary calculations.

The user only specifies the interval at which results are written,

for example

tout = 0.5 fs

while the internal integration step is determined automatically by the propagator.

Predictor–Corrector Scheme¶

The propagation algorithm follows a predictor–corrector strategy.

Step 1 — Predictor¶

An initial estimate of the wavefunction at

\[ t+\Delta t \]

is generated.

Step 2 — Hamiltonian Update¶

Using this predicted wavefunction,

the Hamiltonian,

mean-field operators,

and density matrices are recalculated.

Step 3 — Corrector¶

The equations of motion are solved again using the updated quantities,

producing a more accurate wavefunction.

This process greatly reduces numerical errors during long propagations.

Maintaining Orthogonality¶

During propagation,

the SPFs evolve continuously.

Numerical integration can slowly destroy their orthogonality.

Therefore,

after every propagation step,

the SPFs are orthogonalized,

ensuring

\[ \langle \varphi_i | \varphi_j \rangle = \delta_{ij}. \]

Maintaining an orthonormal basis is essential for stable ML-MCTDH calculations.

Conservation of the Wavefunction Norm¶

One of the most important checks during propagation is conservation of the wavefunction norm.

The output continuously reports

Norm = 0.99999979

or

Norm = 1.00000000

The norm is defined as

\[ \boxed{ \langle\Psi|\Psi\rangle =1. } \]

Values very close to one indicate that the numerical integration is stable.

A noticeable drift in the norm usually indicates

insufficient SPFs,
poor convergence,
or numerical instability.

Energy Conservation¶

The output also reports

E-tot

E-corr

Delta-E

For time-independent Hamiltonians,

the total energy should remain nearly constant.

For example,

E-tot = 5.667911 eV

Delta-E = 0.0003 meV

shows that the energy changes by only a tiny fraction of a millielectronvolt, indicating excellent numerical accuracy.

Recursive Propagation in ML-MCTDH¶

In ML-MCTDH,

every node of the multilayer tree possesses its own

coefficients,
density matrix,
SPFs,
equations of motion.

The propagation therefore proceeds recursively.

Starting from the lowest layer,

the algorithm updates each node before moving upward through the tree.

Information is exchanged continuously between parent and child nodes until the entire multilayer wavefunction has been updated.

Relation to the Output¶

Every block printed in the propagation output is generated during this numerical integration.

For example,

Time = 0.50 fs

Norm = 0.99999979

Natural weights

Mode expectation values

Population

corresponds to one successfully completed propagation cycle.

Once the wavefunction has been advanced,

MCTDH computes

natural populations,
expectation values,
electronic populations,
total energy,
normalization,

before continuing to the next time step.

Propagation Completion¶

When the requested propagation time has been reached, MCTDH terminates the integration and prints

Propagation was successful.

Total time [h:m:s]
Wall time [h:m:s]

This confirms that all equations of motion were integrated successfully over the complete simulation interval.

Connection to the Next Section¶

After the wavepacket has been propagated, the resulting time-dependent wavefunction can be analyzed to obtain physical observables such as

electronic populations,
autocorrelation functions,
absorption spectra,
expectation values,
reduced density matrices.

The next section explains how these quantities are extracted from the propagated wavefunction and how they are related to the files generated during the MCTDH calculation.