Multi-configuration time-dependent Hartree

Multi-configuration time-dependent Hartree (MCTDH) is a general algorithm to solve the time-dependent Schrödinger equation for multidimensional dynamical systems consisting of distinguishable particles.

MCTDH can thus determine the quantal motion of the nuclei of a molecular system evolving on one or several coupled electronic potential energy surfaces. It is an approximate method whose numerical efficiency decreases with growing accuracy.

MCTDH is suited for multi-dimensional problems, in particular for problems that are difficult or even impossible to solve in conventional ways.

Methods

MCTDH

Ansatz

$\Psi (q_{i},...,q_{f},t)=\sum _{j_{1}}^{n_{1}}...\sum _{j_{f}}^{n_{f}}A_{j_{1}...j_{f}}(t)\prod _{\kappa =1}^{f}\varphi _{j_{\kappa }}^{(\kappa )}(q_{\kappa },t)$

Where the number of configurations is given by the product $n_{1}...n_{f}$. The single particle functions (SPFs), $\varphi _{j_{\kappa }}^{(\kappa )}(q_{\kappa },t)$, are expressed in a time-independent basis set:

$\varphi _{j_{\kappa }}^{(\kappa )}(q_{\kappa },t)=\sum _{i_{1}=1}^{N_{\kappa }}c_{i_{\kappa }}^{(\kappa ,j_{\kappa })}(t)\;\chi _{i_{\kappa }}^{(\kappa )}(q_{\kappa })$

Where $\chi _{i_{\kappa }}^{(\kappa )}(q_{\kappa })$ is a primative basis function, in general a Discrete Variable Representation (DVR) that is dependent on coordinate $q_{\kappa }$. If $n_{1}...n_{f}=1$, one returns to the Time Dependent Hartree (TDH) approach. In MCTDH, both the coefficients and the basis function are time-dependent and optimized using the variational principle.

Equations of motion

Lagrangian Variational Principle

$L=\langle \Psi |i{\frac {\partial }{\partial t}}-H|\Psi \rangle$

Where:

$\delta \int _{t_{1}}^{t_{2}}L{\text{d}}t=0$

Which is subject to the boundary conditions $\delta L(t_{1})=\delta L(t_{2})=0$. After integration, one obtains:

${\text{Re}}\langle \delta \Psi |i{\frac {\partial }{\partial t}}-H|\Psi \rangle =0$

McLachlan Variational Principle

$\delta ||i{\frac {\partial }{\partial t}}-H\Psi ||^{2}=0$

Where only the time derivative is to be varied. We can rewrite this norm squared term as a scalar product, and vary the bra and ket side of the product:

${\begin{aligned}0&=\delta \langle i{\frac {\partial }{\partial t}}\Psi -H\Psi |i{\frac {\partial }{\partial t}}\Psi -H|\Psi \rangle \\&=\langle i\delta {\frac {\partial }{\partial t}}\Psi |i{\frac {\partial }{\partial t}}-H|\Psi \rangle +\langle (i{\frac {\partial }{\partial t}}-H)\Psi |i\delta {\frac {\partial }{\partial t}}\Psi \rangle \\&=-i\langle \delta \Psi |i{\frac {\partial }{\partial t}}-H|\Psi \rangle +i\langle (i{\frac {\partial }{\partial t}}-H)\Psi |\delta \Psi \rangle \\&=2{\text{Im}}\langle \delta \Psi |i{\frac {\partial }{\partial t}}-H|\Psi \rangle \end{aligned}}$

Dirac-Frenkel Variational Principle

If each variation of $\delta \Psi ,i\delta \Psi$ is an allowed variation, then both the Lagrangian and the McLanchlan Variational Principle turn into the Dirac-Frenkel Variational Principle:

$\langle \delta \Psi |i{\frac {\partial }{\partial t}}-H|\Psi \rangle =0$

Which simplest and thus preferred method of deriving the equations of motion.

ML-MCTDH

Motivation

The original ansatz of MCTDH generates a single layer tensor tree; however, there is a limit to the size and complexity this single layer can handle. This prompted the development of a multilayer (ML)-MCTDH ansatz by Manthe which was then generalized by Vendrell and Meyer.

Tensor Tree Formalism

Multiple layers are generated through the creation of a tensor tree of nodes linking the modes (DOFs). Solving the tree layout is an NP-hard problem, but strides have been taken to automate this process through mode correlations by Mendive-Tapia.

Example MCTDH tree with l representing layers and q1-6 being the modes.

Ansatz

The generalized ML expansion of Meyer can be written as follows:

${\begin{aligned}\varphi _{m}^{z-1,\kappa _{l-1}}(q_{\kappa _{l-1}}^{z-1})&=\sum _{j_{1}=1}^{n_{1}^{z}}...\sum _{j_{p^{z}}=1}^{n_{p^{z}}^{z}}A_{m;j_{1},...,j_{p^{z}}}^{z}\prod _{\kappa _{l}=1}^{p^{z}}\varphi _{j_{\kappa _{l}}}^{z,\kappa _{l}}(q_{\kappa _{l}}^{z})\\&=\sum _{J}A_{m;J}^{z}\cdot \Phi _{J}^{z}(q_{\kappa _{l-1}}^{z-1})\end{aligned}}$

Where the coordinates are combined as

$q_{\kappa _{l-1}}^{z-1}=(q_{1}^{z},...,q_{p^{z}}^{z})$

Equations of motion

Where the equations of motion are now represented as follows:

${\begin{aligned}i{\frac {\partial A_{J}^{1}}{\partial t}}&=\sum _{K}\langle \Phi _{J}^{1}|{\hat {H}}-\sum _{\kappa _{1}=1}^{p^{1}}{\hat {g}}^{1,\kappa _{1}}|\Phi _{K}^{1}\rangle A_{K}^{1}\\&=\sum _{K}\langle \Phi _{J}^{1}|{\hat {H}}|\Phi _{K}^{1}\rangle A_{K}^{1}-\sum _{\kappa _{1}=1}^{p^{1}}\sum _{i=1}^{n_{\kappa _{1}}^{1}}{\hat {g}}_{j_{\kappa _{1}}i}^{1,\kappa _{1}}A_{j_{1}...i...j_{p^{1}}}\end{aligned}}$

The SPF EOMs are formally defined the same for all layers:

$i{\frac {\partial \varphi _{n}^{z,\kappa _{l}}}{\partial t}}=(1-P^{z,\kappa _{l}})\sum _{j,m=1}^{n_{\kappa _{l}}^{z}}(p^{z,\kappa _{l}})_{nj}^{-1}\cdot \langle {\hat {H}}\rangle _{jm}^{z,\kappa _{l}}\varphi _{m}^{z,\kappa _{l}}+\sum _{j=1}^{n_{\kappa _{l}}^{z}}{g_{jn}^{z,\kappa _{l}}\varphi _{j}^{z,\kappa _{l}}}$

Where ${\hat {g}}$ is a Hermitian gauge operator defined as follows:

$\langle \varphi _{j}^{z,\kappa _{l}}|i{\frac {\partial }{\partial t}}\varphi _{k}^{z,\kappa _{l}}\rangle =\langle \varphi _{j}^{z,\kappa _{l}}|{\hat {g}}^{z,\kappa _{l}}|\varphi _{k}^{z,\kappa _{l}}\rangle =g_{jk}^{z,\kappa _{l}}$

Examples of uses in literature

NOCl

The first verification of the MCTDH method was with the NOCl molecule. It's size and asymmetry makes it a perfect test bed for MCTDH: it is small and simple enough for its numerics to be manually verified, yet complicated enough for it to already squeeze advantages against conventional product-basis methods.

Water clusters

The solvation of the hydronium ion is a topic of continued research. Researchers have been able to successfully use MCTDH to model the Zundel and Eigen ions in close agreement with experiment.

Limitations

Approximate Degree of Freedom Allowance for Each Computational Method

Method	Degrees of Freedom Possible
Conventional Methods (e.g. TDH)	6
MCTDH	12
ML-MCTDH	24+
ML-MCTDH with the Spin-Boson Model	1000+

For a typical input in ML-MCTDH to be run, a node tree, potential energy surface, and equations of motion must be generated by the user. These prerequisites—along with total compute time—soft-cap the size of systems able to be studied with ML-MCTDH; however, advances in neural networks have been shown to address the difficulty of the generation of potential energy surfaces. These issues can also by circumvented by using the spin-boson or other similar bath models that do not pose the same assignment challenges.

Software packages implementing the MCTDH method

Package Name	Group	University	Link
Heidelberg MCTDH	TC Group	Heidelberg University	Link
QUANTICS	Worth	UCL	Link
MCTDH-X	N/A	ETH Zurich	Link

Example Usage of the Heidelberg Package for NOCl

Input and Operator File

nocl0.inp

nocl0.op

RUN-SECTION relaxation tfinal= 50.0 tout= 10.0 name = nocl0 overwrite output psi=double timing end-run-section OPERATOR-SECTION opname = nocl0 end-operator-section SBASIS-SECTION rd = 5 rv = 5 theta = 5 end-sbasis-section pbasis-section #Label DVR N Parameter rd sin 36 3.800 5.600 rv HO 24 2.136 0.272,ev 13615.5 theta Leg 60 0 0 end-pbasis-section INTEGRATOR-SECTION CMF/var = 0.50 , 1.0d-5 BS/spf = 10 , 1.0d-7 SIL/A = 12 , 1.0d-7 end-integrator-section INIT_WF-SECTION build rd gauss 4.315 0.0 0.0794 rv HO 2.151 0.0 0.218,eV 13615.5 theta gauss 2.22 0.0 0.0745 end-build end-init_wf-section ALLOC-SECTION maxkoe=160 maxhtm=220 maxhop=220 maxsub=60 maxLMR=1 maxdef=85 maxedim=1 maxfac=25 maxmuld=1 maxnhtmshift=1 end-alloc-section end-input

OP_DEFINE-SECTION title NOCl S0 surface end-title end-op_define-section PARAMETER-SECTION mass_rd = 16.1538, AMU mass_rv = 7.4667, AMU end-parameter-section HAMILTONIAN-SECTION --------------------------------------------------------- modes | rd | rv | theta --------------------------------------------------------- 0.5/mass_rd | q^-2 | 1 | j^2 0.5/mass_rv | 1 | q^-2 | j^2 1.0 | KE | 1 | 1 1.0 | 1 | KE | 1 1.0 |1&2&3 V --------------------------------------------------------- end-hamiltonian-section LABELS-SECTION V = srffile {nocl0um, default} end-labels-section end-operator

Output absorption spectrum

The absorption spectrum for the NOCl molecule on excitation to the S1 state