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

$\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 primitive basis function, in general a Discrete Variable Representation (DVR) that is dependent on coordinate $q_{\kappa }$ .^[1] If $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

$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.^[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:

${\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. Its 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

More information Method, Degrees of Freedom Possible ...

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]

Close

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

More information Package Name, Group ...

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]

Close

Example Usage of the Heidelberg Package for NOCl

Input and Operator File

More information nocl0.inp, nocl0.op ...

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

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

Close

Multi-configuration time-dependent Hartree

Methods

Basic algorithm

Wavefunction expansion

Equations of motion

Lagrangian Variational Principle

McLachlan Variational Principle

Dirac-Frenkel Variational Principle

Multi-layer extension

Motivation

Tensor Tree Formalism

Wave function expansion

Equations of motion

Examples of uses in literature

NOCl

Water clusters

Limitations

Software packages implementing the MCTDH method

Example Usage of the Heidelberg Package for NOCl

Input and Operator File

Output absorption spectrum

References

Further reading

Related Articles

Related Articles

"Optimal Mode Combination in the Multiconfiguration Time-Dependent Hartree Method through Multivariate Statistics: Factor Analysis and Hierarchical Clustering"

"The coupling of the hydrated proton to its first solvation shell"

"Quantum Phase Transition in the Spin-Boson Model: A Multilayer Multiconfiguration Time-Dependent Hartree Study"