Multi-configuration time-dependent Hartree

Multi-configuration time-dependent Hartree (MCTDH) is an approach to quantum molecular dynamics, an algorithm to solve the time-dependent Schrödinger equation for multidimensional dynamical systems consisting of distinguishable particles. The nuclei of molecules is one example of such particles and their vibrational motion is a form of time-dependence. The method uses an overall wavefunction composed of products of single-particle wavefunctions as first proposed by Douglas Hartree in 1927. The "multiconfiguration" part of the method refers to combining multiple such products.^[1]: 37

MCTDH can predict the 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.^[2]

MCTDH is suited for multi-dimensional problems, in particular for problems that are difficult or even impossible to solve in conventional ways.^{[citation needed]}

Methods

Basic algorithm

Wavefunction expansion

${\displaystyle \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 ${\displaystyle n_{1}...n_{f}}$. The single particle functions (SPFs), ${\displaystyle \varphi _{j_{\kappa }}^{(\kappa )}(q_{\kappa },t)}$, are expressed in a time-independent basis set:

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

Equations of motion

Lagrangian Variational Principle

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

Where:

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

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

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

McLachlan Variational Principle

${\displaystyle \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:

${\displaystyle {\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 ${\displaystyle \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:

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

Which simplest and thus preferred method of deriving the equations of motion^[1].

Multi-layer extension

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^[4] which was then generalized by Vendrell and Meyer^[5].

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.^[6]

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

Wave function expansion

The generalized ML expansion of Meyer^[5] can be written as follows:

${\displaystyle {\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

${\displaystyle 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:

${\displaystyle {\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:

${\displaystyle 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 ${\displaystyle {\hat {g}}}$ is a Hermitian gauge operator defined as follows:

${\displaystyle \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.^[7]

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^[8] and Eigen^[9] 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[2]
ML-MCTDH	24+[5]
ML-MCTDH with the Spin-Boson Model	1000+[10]

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.^[11] 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.^[12] These issues can also by circumvented by using the spin-boson or other similar bath models that do not pose the same assignment challenges^[10].

Software packages implementing the MCTDH method

Package Name	Group	University	Link
Heidelberg MCTDH	TC Group	Heidelberg University	Link[13]
QUANTICS	Worth	UCL	Link[14]
MCTDH-X	N/A	ETH Zurich	Link[15]

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