# Hamiltonian/Liouvillian Documentation

## General Structure

The operator file is divided into sections containing information for the definition of the Hamiltonian. The structure is, analogous to the input file, based on keywords with arguments. See Input Documentation for further details on keyword placement, special characters etc.

The operator file ends with the keyword

Central to the operator definition is the HAMILTONIAN-SECTION. Here, the Hamiltonian is defined, term for term, using labels. These labels may be real numbers or alphanumeric strings. Various simple mathematical operations can be used to manipulate the labels. The program reads this input and translates it into the form needed, replacing the alphanumeric labels with real numbers or operators as required. Real numbers may also be defined as alphanumeric parameters in the PARAMETER-SECTION. Many operators (e.g. the position operator,q, and simple functions of this) are known internally. Other operators are defined in the LABELS-SECTION. This includes cases were the operators are read from a file, or when the operators need to be parametrised using real numbers or parameters.

## Sections

The section XXX starts with the keyword

and ends with

The various sections are listed below. There is no predefined order for the sections.

XXX | Description |
---|---|

OP_DEFINE | Definition of the operator. |

PARAMETER | Definition of any parameters used in the operator. |

HAMILTONIAN | Definition of the Hamiltonian. |

Dissipation | Definition of the dissipative Liouvillian. |

LABELS | Definition of operators needed. |

## Keywords

Below are tables of keywords for each operator section.

The first table in each section describes the keywords. The number and type of arguments is specified. The type is S for character string, R for real number, and I for integer, e.g.

indicates that the keyword takes two arguments, the first is a string, the second a real number. The status column indicates whether the keyword is compulsory (C), or optional (O).

The second table defines the default values for optional keywords and arguments.

## Op_define-Section

Keyword | Description |
---|---|

title
S end-title |
S is a title for the operator. There is no restriction on the format or number of lines used. This title will be printed e.g. in the log file. |

## Parameter-Section

Parameters are defined using the format:

where parameter is any string (case sensitive).

Internally the program uses atomic units. Parameters in other unit systems can be converted by a relevant unit keyword. Various units are recognized by the program. See Units for a full list. For example to define the value 0.1825eV to a parameter with name lambda,

It is often very convenient to define a parameter as an
arithmetic expression which may contain other (previously defined)
parameters. The arithmetic expression **must not contain
spaces or brackets!** (Spaces are now allowed for arithmetic
in a parameter-section, but not for arithmetic anywhere else). Only
the operators +-*/^ are allowed. Note: An exponent acts only on the
label to which it is attached and */ are evaluated before +-.
Otherwise the order of operation is strictly from left to right
(See also *Hamiltonian-Section* below.) Exponents may be
integer or real signed numbers but not parameters, i.e
`a^-0.5` is allowed whereas `a^b` is an invalid
expression. All parameters are real numbers. Numerical values may
be input in fixed format (e.g 1224.0) or in exponential format
(e.g. 1.224d+3). The exponential format is characterized by the
letter d. It must not be replaced by D, e or E. The exponential
format is not allowed for exponents. A arithmetic expression may be
followed by a unit, but individual entries must not carry units.
E.g `c = 2.0,ev/5.3,cm-1` is invalid. One should use:

Parameters should not be reused, i.e. expressions like

must **not** be used.

If a parameter is differently defined in different parts of the
input, the following rule applies: Parameter definitions as options
(-p option) overwrites definitions provided in a PARAMETER-SECTION
of the input file. The latter overwrite an *alter-parameter*
block and finally the definitions provided in a PARAMETER-SECTION
of the operator file are overwritten by all other definitions.

#### Example

PARAMETER-SECTION mass_rd = 2.0/3.0, H-mass # equivalent to : mass_rd = 1224.766667 mass_rv = 0.5, H-mass # equivalent to : mass_rv = 918.575 jtot = 4.0 centri = jtot^2/2.0/mass_rd+jtot/2.0/mass_rd end-parameter-section

Parameter values cannot only determined by simple arithmetic but also by employing functions. The syntax is

par = FUNCTION[expression]

where FUNCTION is one of the keywords:

These are the usual Fortran routines. Note that the last three
functions depend on two variables. The *expression* is an
arithmetic expression containing numbers and/or parameters.

#### Example

par1 = SIN[ degree*PI/180.0 ] par2 = MAX[ alpha, 3.0*beta-4.5 ] enatural = EXP[ 1.0 ]

### Special Parameters

Some parameters have special meanings.

PI | Pre-defined parameter. PI = 3.141592654... |

mass_label | The value of this parameter is taken as the mass of the
degree of freedom specified by label (as defined in the
PRIMITIVE-BASIS-SECTION of the
input file. This mass is then used in the KE operator etc. By
default all masses are otherwise set to 1.0 au. |

jtot | The total angular momentum. This is used in the adiabatic correction. |

jbf | The body-fixed angular momentum. If this parameter is specified, it is assumed that the coupled-states approximation is being used. The magnetic quantum numbers for the Legendre functions of both the primitive and single-particle bases are then set to this value. See PRIMITIVE-BASIS-SECTION and Building the initial wavefunction |

### Changing Parameter Values

The value of parameters can be defined not only in this section, but alternatively:

- in a separate file by using the the
*parfile = filename*keyword in the OPERATOR-SECTION of the input file. - in the OPERATOR-SECTION of
the input file using the
*alter-parameters*keyword. - in a PARAMETER-SECTION , a section which may also appear in the input file.
- from the command line, by starting the program with
`mctdh<ver> -p apar rpar`

The order of precedence of parameters defined from different sources is command line, input file, parameter file, operator file.

Thus parameters in the operator file are the default values, which can be altered in a run in a variety of ways.

## Hamiltonian-Section

A Hamiltonian-Section is opened by the keyword

where *xxx* stands for a string naming the operator to be
defined. The first variant is for defining the system Hamiltonian,
whereas the others, referred to as *named operators*
are for defining operators used for other
purposes than propagation (e.g. for calculating expectation values, or
generating eigenfunctions which
serve as initial single-particle functions). A Hamiltonian-Section
is closed by the keyword:

end-hamiltonian-section

Note: **Cases do matter within a
Hamiltonian-Section!**

The first lines of this section is reserved for keywords that can define how the operator is set up and used.

Keyword | Description |
---|---|

nodiag /
usediag |
This keyword controls the implementation of unit operators.
If usediag is used matrix elements of unit operators
are not explicitely used. This is faster, and numerically more
accurate. However it does assume that the bra and ket of the
matrix are the same orthonormal basis. If nodiag is
given, matrix elements of unit operators are explicitly
calculated. This is necessary if the bra and the ket
wavefunctions are different. Such a situation occurs when using
the operate keyword (see Init_wf-Section) or when using
correlated operators for flux analysis. By default the system
Hamiltonian is set with usediag, and all other operators
nodiag. Thus these keywords are usually needed only when
wanting to make use of the diag feature using a
non-system operator. |

addmode=S(,S1(,S2..)) | All DOFs not defined in the SPF-BASIS-SECTION are usually
ignored when the operator is build. The addmode
keyword allows to include additional DOFs. Those DOFs, the
modelabel of which is S ( or S1, or..) are included. This is
useful when building operators to be used for eigenfunction
generation in the Init_WF-Section. See the note in the description of the
Initial_WF-Section. |

The following line(s) of this section define the order of the degrees of freedom in the term specifications. The format for this is

where S1, S2 .... label the degrees of freedom treated by the operator. The labels must agree with the mode labels defined for the primitive basis, but the order may differ. As many lines as necessary may be input, as long as the first keyword in each line is "modes". In named operators, only relevant modes need to be included.

After the modes have been defined, each term of the Hamiltonian is expressed in symbolic form by

where A is a coefficient, and *O* a one-dimensional
operator for the n-th degree of freedom (corresponding to the n-th
mode label in the list defined at the beginning of the section).
See the LaTeX note Built-in
Symbols for a list of possible choices of *O* .

If any degrees of freedom for a term have a unit operator (i.e.
O = 1), then they do not need to be explicitly included. The
separator between the degrees of freedom before and after those not
included must then be changed from "|" to "|`n`", where
`n` is the number of the degree of freedom after those
missed out, i.e. the following lines are equivalent. Note that the
number `n` corresponds to the order in the modes definition
line at the start of the operator.

A |1 O1 |3 O3 | .....

NB there must be no blank spaces in the separator keyword
|`n`. Each Hamiltonian line must be either un-numbered
(i.e. "|") or numbered (i.e. "|`n`"), one must not use
numbered and unnumbered input in one line.

This formalism is also used to include multi-dimensional functions. For example, if the function O1 is a function of the first and second degrees of freedom, one could again use an unnumbered or a numbered Hamiltonian line:

A |1&2 O1 |3 O3 | .....

As a line must not have more than 240 characters, one may use
continuation lines. A continuation line starts with a
*&&&* symbol in the coefficient column. The
following examples are equivalent.

A |1 O1 |2 O2 |3 O3 |4 O4 |5 O5 |6 O6 A |1 O1 |2 O2 |3 O3 &&& |4 O4 |5 O5 |6 O6 A |1 O1 &&& |2 O2 |3 O3 &&& |4 O4 |5 O5 &&& |6 O6

An unnumbered input is also possible. Remember, mode lines are continued by setting the keyword "modes" rather than "&&&".

The coefficient A can be a product of real numbers and parameters. The operators O can be a product of real numbers, parameters, and operators. The factors are separated by a '*' or '/' for multiply or divide respectively, i.e. A = a*b/c. A '-' may also be placed in front of the product. Parameters are defined in the PARAMETER-SECTION, operators are defined internally or in the LABELS-SECTION. NB Real numbers must contain a decimal point.

Powers of numbers and parameters (and also of certain operators) are also possible. This is input as a^b, where b is an integer or a real number. Integers may be positive or negative. These act only on the label to which they are attached, e.g. c*d^r = c*(d^r) or c^r*d^r = (c*d)^r.

The program works factor by factor and so c/d/e = c/(d*e),
whereas c/d*e = (c/d)*e. NB There must be **no
spaces** in the product e.g. if a*b * c is input, the
program would not recognize the factor c. Moreover, there must be
**no brackets** and **no sums** within a
HAMILTONIAN-SECTION.

### Rules:

- All numbers, parameters and a possible sign must be taken care of by the coefficient, i.e. they must appear in the first column.
- The operator columns must contain only simple operator symbols (without parameters) and products thereof. No sums, no signs.
- There must be no lines in the Hamiltonian-Section which differ only in the coefficient. If there are such redundant operator lines, one should sum the coefficients (e.g. in the Parameter-Section) and keep only one of those lines.
- The special Symbol
**I**(imaginary unit) stands for an operator (see Built-in Symbols), it is not a parameter! Parameters are always real! The symbol**I**may appear in a HAMILTONIAN-SECTION only. It may appear in the coefficient column or in any DOF column, including the*Time*column. - Natural potentials (natpots) and potentials in CPD/CANDECOMP
(cpd-file) format and Multi-Layer Potfits are special. They must not be
multiplied with another operator, e.g.:
`1.0 | q*V`

is invalid, if V is a natpot/cpdfile/mlpf. However,`-0.1*I | V`

is a valid construct, i.e. a natpot/cpd/mlpf may be multiplied with a (real or imaginary) coefficient. - The assignment of the natpot/cpd/mlpf to the various DOFs is done
through the modelabels. Hence a Hamiltonian line for a natpot/cpd/mlpf
usually simply reads
`1.0 | V`

. In case one wants to let the natpot/cpd operate on DOFs different from the selection or order given by the modelabels, one may make use of the`|`

construct. In this case one must set*j*&*k**ignore*in the natpot/cpdfile statement of the LABELS-SECTION. In contrast, mlpfs don't support the*ignore*keyword; nevertheless one can use the`|`construct if the mlpf only applies to a subset of the system DOFs. If the order of the system DOFs and the mlpf DOFs differ, the mlpf will be automatically adjusted accordingly, where the DOFs are matched by modelabels.*j*&*k*

### LMR operators:

One may use operator products, e.g. `q*dq`

, but this
works only if all operators of the product are DOF operators. But
if there is a mode operator, `v`, then the construct
`v*dq`

does not work, because it is unclear on which
coordinate of the combined mode `dq` should operate. To
solve this problem, LMR operators are introduced. Three operators,
a left one (L), a middle one (M) and a right hand one (R) may build
an operator product. The middle operator may be left out. To build
the operator

where (x,y) are in a combined mode, one writes:

const |1L dq |1&2M v |2R dq |3 ...

Note that the LMR operators of one line must operate on one combined mode exclusively. Currently the LMR operators must be potential or matrix operators. No complex operators (CAPs), no tensor operators (KLeg operators) and no natpots/cpd/mlpf.

Note that LMR operators are also helpful when an FFT is used as
primitive basis. In this case an operator product like
`q*dq`

is not allowed, because `dq`

is not
represented by a matrix if an FFT is used. However, writing this
product as an LMR operator works:

const |1L q |1R dq |2 ...

## Dissipative-Section

The Dissipative-Section applies only to the propagation of density operators. It allows the treatment of open systems in Markov approximation, i.e. the dissipative functional (Liouvillian) is assumed to be local in time. The format of the section is similar to the Hamiltonian section. The first line defines the order in which the degrees of freedom are addressed:

DISSIPATIVE-section modes | X1 | X2 | ... | Xf | X1 | X2 | ... | Xf ... end-DISSIPATIVE-section

Here X1, X2,..., Xf label the f degrees of freedom. Each degree of freedom is defined twice due to the structure of the dissipative Liouvillian. Typically, three different types of expressions can occur.

Type 2: rho A B

Type 3: A rho B

Thus, the first f columns define the operator A, and the second f columns define the operator B. A distinction between the three types is achieved by additional keywords, which should be given in otherwise empty lines:

*pure-left*

Type 2:

*pure-right*

Type 3:

*left-right*

The order of the labels within the line has to be the same for the first and the second part. Furthermore, the labels must agree with the mode labels defined for the primitive basis, and the order within a single part should be chosen as in the Hamiltonian-Section.

A useful example is the Henon-Heiles operator including
dissipation (A. Raab and H.-D. Meyer, J. Chem. Phys
**112** (2000) 10718), where the dissipative
Liouvillian

L[rho] = -c([x,[x,rho]]+[y,[y,rho]]) = -c( x*x*rho -2*x*rho*x + rho*x*x + y*y*rho -2*y*rho*y + rho*y*y)

is implemented:

modes | X | Y | X | Y pure-left -cdiss | 1 | q | 1 | q -cdiss | q | 1 | q | 1 pure-right -cdiss | 1 | q | 1 | q -cdiss | q | 1 | q | 1 left-right 2.0*cdiss | 1 | q | 1 | q 2.0*cdiss | q | 1 | q | 1

The first term in the dissipative Liouvillian is coded in the fourth line of the example, the second term in line ten, the third term in line seven, and so on. Note: The fourth line, e.g., could as well be written:

-cdiss | q^2 | 1 | 1 | 1

This applies to all lines of the *pure-left* and
*pure-right* cases.

Special **Note** for the propagation of density
operators of type I: If the product A*B of two single operators A,
B referring to the **same** degree of freedom appears
in the Liouvillian, then the MCTDH program multiplies them
internally and uses the exact product. To enable the
multiplication, additional labels have to be given in the Labels- Section (see also the example above).

Labels-Section ... A%B = A*B ... end-labels-section

For the propagation of density operators of type II it is
necessary to represent, for example, the operator product A*B as
two single operators A, B, because of the *projection*
option (c.f. the Operator-Section
in the input documentation).

#### Remarks

- In the output-file, both Norm and Trace are given. The Trace is in any case the trace over the density matrix. The Norm is the "Hilbert-Schmidt-Norm", i.e. the square-root of the trace over the square of rho.
- The density1 case does only work if the product operator of A*B exist explicitly. This means that, e.g., "q%q = q^2" and "rai%low = num" are alright, but not "q%dq = q*dq".
- The "auto" keyword in the Run-Section does not give useful results for density matrix calculations. Use "ctrace".

## Labels-Section

The program needs to translate the labels used in the HAMILTONIAN-SECTION into operators. Labels corresponding to certain simple operators do not need to be defined as they are known internally to the program. For internally defined symbols that are in this category, see Built-in Symbols .

Other labels must be defined in this section. This definition may be used to parametrise an operator, or to include difficult operators via a file. The definition of labels defined here can also be changed in the OPERATOR-SECTION of the input file using the alter-labels keyword.

In order to add new operators, see the description of the operator functions library opfuncs.

The format for the definition, one definition per line, is one of the following depending on the operator:

label = operator [par1, par2, ... ]

label = operator {opt1, opt2, ... }

where label is the label used in the HAMILTONIAN-SECTION, and operator is the name of the operator to be attached to the label.

The first form is used for operators which do not require arguments or options. This can be used to redefine an operator in the input file.

The second form is used to define an operator with parameters,
where [ par1, par2, ... ] are the relevant parameters, enclosed in
square brackets and (optionally) separated by commas. Parameters
can be made up of names listed in the PARAMETER-SECTION, using simple arithmetic rules, or
real numbers. E.g.: the input ```
CAP[ 3.0+x0 1.0d-3*strength 3 1
]
```

is valid. The parameters x0 and strength must, of course,
be previously defined. The use of units is not allowed here. For
possible operators, see Built-in
Symbols.

The third line is to define operators with optional arguments, where { opt1, opt2, .. } are the relevant parameters, enclosed in curly brackets and (optionally) separated by commas.

The following operator types relate to operators stored in files.

Operator | Description |
---|---|

natpot {pname,S} | Natural potential fit is read from the file pname/natpot.
If the keyword "name" is used for pname, the natpot file is
expected to reside in the name-directory. A natpot uses the
modelabels to define on which DOFs it will operate. If the
optional string S=ignore is given, the modelabel information is
ignored. This requires that the DOFs on which the natpot shall
operate are defined in the Hamiltonian-Section through the
|i&j&k (i,j,k=1,2,...) construct. |

cpdfile {pname, S} | A fit of a potential in Canonical Polyadic Decomposition
(CPD) form is read from the file pname/cpd.
If the keyword "name" is used for pname, the cpd file is
expected to reside in the name-directory. A cpdfile, like natpots,
uses the modelabels to define on which DOFs it will operate.
If the optional string S=ignore is given, the modelabel information
is ignored. This requires that the DOFs on which the cpdfile
shall operate are defined in the Hamiltonian-Section through
the |i&j&k (i,j,k=1,2,...) construct. |

mlpf {pname} | Multi-Layer potential fit is read from the file pname/mlpf.dat. An mlpf uses the modelabels to define on which DOFs it will operate. See notes on MLPF below. |

read {file,n} | n th Operator is read from file. For the file format see note below. |

readmap{file,n} | n th Operator is read from file. The file name should end with 'readmap'. For example, XXX = readmap{xyz_readmap,n}. For the file format see note below. |

readsrf {file,S} | S=ascii or S=binary. A multi-dimensional potential energy
surface is read from file. Format: One energy per line (i.e per
read statement). For POTFIT, the order is implicitly defined by
the Primitive-Basis-Section. But in mctdh the order is the one
of the WF, i.e. the order of the SPF-Basis-Section. The
|i&j&k (i,j,k=1,2,...) construct of the
Hamiltonian-Section must reflect this order. For readsrf one
always should use a |i&j&k (i,j,k=1,2,...)
construct and not a simple |. Note: The used order
of the DOFs is written to log-file. The first index (DOF) runs
fastest. There is no check for the consistency of the
data. |

read1d {file,S} | S=ascii or S=binary. A one-dimensional potential energy
curve, operating on one DOF, is read from file. Format: One
energy per line (i.e per read statement). read1d is
very similar to readsrf. Note that read1d
must be used when the potential is one-dimensional. |

external1d {file} | A file is read defining a 1D real function (potential). The
data is to be provided in two columns in ASCII format, i.e.
x . The data is
quadratically interpolated. Blank lines and lines beginning
with a hash '#' are ignored. The _{i} V_{i}x-data must be spaced
equidistantly. There should be at least two data points below
the first grid-point and two above the last grid-point,
respectively. external1d is particularly useful when
numerical time data is to be given to the program. As the time
grid points are not known at the beginning, the interpolation
feature is essential. If a 1D potential is to be read and fitting
is not necessary, read1D should be used. |

shepard | Read the full-dimensional potential as a modified Shepard interpolation. This uses the formalism and data structures developed in the Collins group; the read files are POT, IN_SYSTEM, IN_CNT, IN_ATOMPERMS and possibly CON. Second order Taylor series are used around each data point. ("grow" is a synonym.) |

#### Note on read:

The matrix representations of 1D operators is read from file. The representation is with respect to the DVR basis functions (if a DVR and not e.g. sphFBR is used). See review Appendix B Eqs. (B16,B30) for their definitions.

The structure of the matrix is defined by the value of
*type*, type=1 : full complex matrix; type=2 : real diagonal
matrix; type=3 : full real matrix. The file may be in ascii or
binary format. If in ascii format, the first line must contain the
word *ascii* (starting at column 1). Ascii and binary files
must then contain the following data. (The text in [..] is comment.
Do not write it to file.)

```
nop [number of operators on file]
gdim(1), type(1) [grid-dimension and type of first operator on file. Type=1 -> complex matrix, type=2 -> real diagonal, type=3 -> real
matrix .]
...
...
gdim(nop), type(nop)
[ if type=1 then ]
(( dble(oper(g,g1,1)),g=1,gdim(1)),g1=1,gdim(1))
((dimag(oper(g,g1,2)),g=1,gdim(1)),g1=1,gdim(1))
[ if type=2 then ]
(dble(oper(g,g)),g=1,gdim(1))
[ if type=3 then ]
(( dble(oper(g,g1,1)),g=1,gdim(1)),g1=1,gdim(1))
[Similar for n = 2,..,nop]
```

**NOTE 1:**If you use an ASCII format and if the matrix (type 1 or 3) is written in matrix form, then the transposed of that matrix is read, because the first column read is stored into the first row, etc.

**NOTE 2:**One may read a type 3 operator for a combined mode. In this case gdim is to be replaced by subdim. For diagonal mode operators (type 2) please use readsrf.

*Example* (Read file in ascii free-format)

ascii 2 2 3 3 2 1.5 1.9 1.9 2.3 1.3 2.2 3.4

#### Note on muld-potentials:

Multi-dimensional potentials, called *muld* in MCTDH, do
not satisfy the product structure and their application is hence
slow. Furthermore muld potentials are subject to certain
restrictions. They cannot be multiplied with other operators,
except constants and − if multiset is used − diagonal
state operators (e.g.: S1&1). Mulds can, however, be used
together with time-dependent fields. The DOFs the muld is operating
on are defined by the `|i&j&k` (i,j,k=1,2,...)
construct. This means that the first argument of the muld operates
on the i-th DOF, the second on the j-th, etc. Note that `|
V` is equivalent to `|1&2...&f V` where f is
the number of DOFs and V is a muld. I.e. `| V` operates on
*all* DOFs.

Mulds can operate on a subset of DOFs, but this subset must not
contain "holes", i.e. the DOFs the muld is operating on must be
consecutive in the wavefunction. A muld must also not operate on a
subset of a combined mode. For *readsrf* the DOFs cannot be
re-ordered. The order of the DOFs on the readsrf file must be
identical to the order in the wavefunction. The latter order is
defined by the SPF-BASIS-SECTION. (For potfit, however, the order
is defined by the PRIMITIVE-BASIS-SECTION). Note that the ordering
of the *modes* line of the HAMILTONIAN-SECTION may be
different. The `|i&j&k` construct refers to the
ordering defined in the *modes* line.

If a muld operates entirely on one combined mode, it becomes a mode operator (rather than a muld) and all the restrictions discussed above to not apply.

In general the use of muld potentials should be avoided. One should use potfit to transform a muld potential to product form.

#### Note on readmap:

The information about the matrix representation of 1D operators is read from file. This is useful for a very sparse matrix. Such matrices appear in the SQR form of MCTDH when they represent products of combined annihilation/creation operators. The representation is with respect to the DVR basis functions. The file should be in ascii format and must contain the following data.

nop
[number of operators on file]

nel(1)
[number of non-zero elements of the first
operator ]

...

...

nel(nop)

(row_index(g,1),g=1,nel(1))
[row indices of the non-zero elements of the first operator]

(column_index(g,1),g=1,nel(1))
[column indices of the non-zero elements of the first operator]

(oper(row_index(g,1),column_index(g,1),1),g=1,nel(1))
[values of the corresponding matrix element of the first operator]

[ Similar for n = 2,..,nop]

readmap operators are subject to certain restrictions.

**NOTE 1:**A readmap operator must appear as a single term, it cannot be multiplied with or added to other operators.

**NOTE 2:**There must not be two mapping operators with identical non-zero rows and columns. One has to add such operators 'by hand' outside of Quantics.

**Example**(Read file in ascii free-format)

2

3

2

1 5 6

1 4 9

1.0 1.0 1.0

3 5

6 2

1.0 -1.0

#### Notes on MLPF:

A Multi-Layer Potfit (mlpf) is a potential fit in hierarchical tensor format. It can only be used with ML wavefunctions. Its use is very similar to a natpot, and if the potential term covers more than three DOFs (or combined modes), the ML-MCTDH propagation is likely to be (much) more efficient with an mlpf instead of a natpot, at negligible loss of accuracy. The computational gains increase rapidly with the number of DOFs and with the accuracy of the fit.

The tree structure of the ML wavefunction (mlwf) and of the mlpf must be compatible, i.e. if mlpf-node A is a child of mlpf-node B, then mlwf-node A' must be a descendant (but not necessarily a direct child) of mlwf-node B', where A' and B' are the mlwf-nodes corresponding to the mlpf-nodes A and B. In other words, the mlwf-tree must contain all the nodes used in the mlpf-tree, in the same hierarchical relationship, but the mlwf-tree may also contain additional nodes (which will be the case if the mlpf only applies to a subset of the system DOFs). Additionally, if mode combination is used, the combined modes used by the mlpf must also be present in the mlwf, and the DOF order in these combined modes must be identical.

A program for generating an mlpf from an existing vpot file or natpot file is available at https://bitbucket.org/frankotto/mlpf with accompanying documentation.

It is important to note that an mlpf is **not** represented in
sum-of-products form, but as a separate operator structure called a
*Multi-Layer Operator* (mlop). Technically, when an mlpf is read
by mctdh, it will be expanded into an mlop that applies to the complete
wavefunction. In this process, the mlpf can be combined with other
operators acting on the non-mlpf DOFs as follows (where V, V1, V2 are
assigned to mlpfs in the LABELS-SECTION):

- multiplication with separable operators, as in
`1 |1&2&3&4 V |5 q |6 KE` - multiplication with another mlpf, as in
`1 |1&2&3&4 V1 |5&6&7&8 V2` - combinations of the above two

*Caveat:* Because mlops constitute an operator structure separate
from the usual sum-of-products operator, most analysis programs that work
with operators may not (yet) work correctly if mlops are used. The
documentation of the respective analysis programs should mention if
mlops are known to be properly supported.

### Changing Label Assignments

The assignment of labels can not only be defined in this section, but alternatively:

- in the OPERATOR-SECTION of
the input file using the
*alter-labels*keyword.

## Time-dependent Operators

A time-dependent operator can be treated as long as any time-dependent terms are separable, i.e. of the form F(t).g(R). To use such an operator, a "degree of freedom" for the time must be added in the HAMILTONIAN-SECTION. This is done by adding a mode with the label Time (note this is case sensitive) to the modes command. This must the last mode. Operators can then be defined for the time coordinate as for any other.

For example, if a two dimensional system is coupled to a time-dependent operator, the .op file may contain something like the following:

PARAMETER-SECTION omega = 0.1, ev end-parameter-section LABELS-SECTION cosom=cos[omega] end-labels-section HAMILTONIAN-SECTION ------------------------------------- modes | X | Y | Time ------------------------------------- 1.0 | KE | 1 | 1 0.5 | q^2 | 1 | 1 1.0 | 1 | KE | 1 0.5 | 1 | q^2 | 1 0.01 | q | q | 1 1.0 | q | 1 | cosom ------------------------------------- end-hamiltonian-section

where the system is two coupled harmonic oscillators, with one dipole-coupled to a radiation field. Note: Time-dependent operators require the use of VMF propagation scheme, if a correlated operator is time-dependent. If --- as in the above example --- the time-dependent operator is uncorrelated, i.e. operates on one mode only, the more powerful CMF integration scheme may still be used. Note that in the operator file the time is in atomic units (au). This is different from the input file, where the time in in fs.

## Examples

### Simple Hamiltonian: The Henon-Heiles System

The Henon-Heiles system is a pair of coupled oscillators. The Hamiltonian can be written:

The operator file to implement this Hamiltonian can be written:

OP_DEFINE-SECTION title Henon-Heiles PES end-title end-op_define-section PARAMETER-SECTION mass_X = 1.0 mass_Y = 1.0 lambda = 0.2 end-parameter-section HAMILTONIAN-SECTION modes | X | Y 1.0 | KE | 1 0.5 | q^2 | 1 -lambda/3 | q^3 | 1 lambda^2/16 | q^4 | 1 1.0 | 1 | KE 0.5 | 1 | q^2 lambda^2/16 | 1 | q^4 lambda | q | q^2 lambda^2/8 | q^2 | q^2 end-hamiltonian-section end-operator

In the first section the degrees of freedom, X and Y, are defined. In the second section the parameters are defined, and finally the operator is symbolically displayed.

## Available Operators

### Selection of Operators

- hh.op: modified Henon-Heiles, 2D, WF
- dhh.op: modified
Henon-Heiles including dissipation (A. Raab and H.-D.
Meyer, JCP (2000)
**112**: 10718) - nocl0.op: NOCl ground
state. Analytic fit to PE surface in Jacobian coordinates using
1-dimensional operators.

nocl1.op: NOCl excited state. Analytic fit to surface in Jacobian coordinates using 1-dimensional operators. - noclT.op: NOCl with time-dependent dipol-interaction between S0 and S1.
- pyrmod.op: Pyrazine
24-mode model from Krempl et al [JCP (94)
**100**: 926]. See also Worth et.al. [JCP (1996)**105**: 4412]. - pyrmod4.op: Pyrazine
4-mode model from Krempl et al [JCP (94)
**100**: 926]. - pyr24+.op: Pyrazine
24-mode bi-linear model [JCP (1999)
**110**: 936]. - pyr4+.op: Pyrazine
4-mode bi-linear model [JCP (1999)
**110**: 936]. - dpyr4.op: Pyrazine
4-mode model including dissipation (A. Raab and H.-D.
Meyer, JCP (2000)
**112**: 10718) - h3j0.op: H +
H
_{2}3D (J=0) in Jacobi coordinates - h3kleg.op: H +
D
_{2}in Jacobi coordinates. Exact kinetic energy (KLeg, 4D) - ch3i0.op: Methyliodide ground state
- ch3i1.op: Methyliodide excited state
- h2surf.op: H2/rectangular lattice - surface scattering
- licn.op: LiCN ground
state. Analytic fit to a 2D SCF potential energy surface (the CN
bond lenght is frozen). The fit is taken from: R. Essers et
al.,
*Chem. Phys. Lett.***89**, 223 (1982). See licnsrf.f. - co2.op: CO
_{2}vibrational (*J*= 0) Hamiltonian, 3D, valence coordinates, PES taken from J. Zuniga et al.,*J. Mol. Spec.***195**, 137 (1999).

### Complete List of Operators

Here you can find a complete list of operators that are contained in the MCTDH package.

## Available Surfaces

The following multi-dimensional surfaces are also available.
These may be used by selecting the appropriate label. Note that the
order of the degrees of freedom in the input file must be in the
same order as defined for the surface. The log-files of
*potfit* and *mctdh* protocol the order actually
used.

Options may be supplied in curly brackets, e.g. lsth{tfac=0.66666 jacobian}. The parameter
*tfac* defines the center of mass of the diatom as a
fraction of its internuclear distance. tfac = m1/(m1+m2).

Inspect the file $MCTDH_DIR/source/opfuncs/funcsrf.F to see how the potential energy surfaces are implemented in detail.

All potential energy surfaces, except for *lsth*, are no
longer part of the MCTDH-package. The surfaces are collected in the
directory *addsurf*. When one of the surfaces is needed,
copy it from *addsurf* to
*$MCTDH_DIR/source/surfaces*. Alternatively one may set
links in *surfaces* pointing to *addsurf*. This is
most conveniently done by running the script *mklinks*.

A typical input line to incorporate a potential via the operator
file in *mctdh* surface looks like

LABELS-SECTION V = bkmp2{binding} end-labels-section

and in *potfit* (input file)

OPERATOR-SECTION pes = bkmp2{binding} end-operator-section

#### readsrf

Read the potential energies from a file. One energy per line.
The order is the order defined in the PRIMITIVE-BASIS-SECTION, or
(in mctdh) by the `|i&j&k` (i,j,k=1,2,...)
construct of the Hamiltonian-Section. Note: first index runs
fastest. The file may be formated (ascii) or unformated (binary).
The arguments are the file name and one of the keywords, ascii or
binary. Warning: There is no check if the ordering of the data is
consistent with the ordering assumed by mctdh or potfit.

Example: pes = readsrf{ ~/MCTDH/potfile ascii }

If the potential to be read is one-dimensional, one has to use*read1d*

*. See Labels Section*

#### bkmp2

The BKMP2 potential energy surface of H_{3} for
2D-collinear or 3D, both in binding (i.e. valence) or Jacobian
coordinates. For more details on the BKMP2 surface see
Boothroyd, Keogh, Martin, Peterson, J. Chem. Phys.
**104**, 7139 (1996) and J. Chem. Phys.
**95**, 4343 (1991).

The following options are available:

jacobian (default)

binding

collinear

#### c2h

The C_{2}H potential energy surface of Martial in
Jacobian coordinates is provided for the 1A' and 2A'' electronic
symmetries. (Keywords: c2ha1 and c2hasec)

#### h4srf

The H_{2} + H_{2} potential energy surface of
Aguado et al. J. Chem. Phys. **101**, 4004,
(1994). Jacobian coordinates: R, r_{1},
r_{2}, θ_{1}, θ_{2}, φ.

#### h4bmkp

The H_{2} + H_{2} potential energy surface of
Boothroyd et al. J. Chem. Phys. **116**, 666,
(2002).

Jacobian coordinates: R, r_{1}, r_{2},
θ_{1}, θ_{2}, φ.

The following options are available:

- tfac1 = A (A=0.5 is default)
- tfac2 = B (B=0.5 is default)
- jacobian (default)
- jacobian2 (for reordered coordinates: R, r
_{1}, θ_{1}, r_{2}, θ_{2}, φ.) - cartesian
- interaction
This only calculates the interaction potential, i.e. the contribution of two isolated H

_{2}is subtracted:V(R,r_{1},r_{2},θ_{1},θ_{2},φ) = BMKP(R,r_{1},r_{2},θ_{1},θ_{2},φ) − BMKP(∞,r_{1},r_{0},π/2,π/2,0) − BMKP(∞,r_{0},r_{2},π/2,π/2,0) + BMKP(∞,r_{0},r_{0},π/2,π/2,0)where r

_{0}= 1.4483 a.u. and ∞ = 100 a.u. ;-)**Warning!**The implementation for this is currently very naive and will increase the time spent for evaluating the potential by a factor 3.

#### h4bmkprr

Like **h4bmkp**, but employing the rigid rotor
model for the hydrogen molecules.

Jacobian coordinates: R, θ_{1},
θ_{2}, φ.

Options:

- r0 = A (fixed bond length; default 1.449)

#### h4dj

The H_{2} + H_{2} potential energy surface of
Diep and Johnson. See J.Chem.Phys. **112**, 4465
(2000).

This is a 4D surface as the H_{2} are treated as rigid
rotors (with bond length = 1.449 a.u.).

Coordinates: R, θ_{1}, θ_{2},
φ

#### hoosrf

HO_{2} surface. DMBE IV of Varandas (J.Phys.Chem.
**94**, 8073 (1990)). The PES is available in
Jacobian, Cartesian, and binding (i.e. valence) coordinates.

The following options are available:

- tfac = A . (tfac=0.5 is default)
- jacobian (default)
- binding
- cartesian

#### hcn

The HCN potential energy surface of J.M. Bowman, B. Gazdy,
J.A. Bentley, T.J. Lee, and C.E. Dateo, J. Chem. Phys.
**99** (1993) 308, is implemented in Jacobian
coordinates: R=H-(CN) distance, cos(theta), r=CN distance.

#### licnsrf

The LiCN surface (see R. Essers, J. Tennyson, and P.E.S.
Wwormer, Chem. Phys. Lett. **89**, 223 (1992))
is actually a collection of 1D potential curves V_{l}(R),
with V(R,theta) = Sum V_{l}(R) * P_{l}(cos(theta)),
where R=Li-(CN) distance and theta is the Jacobian angle. Keyword:
**licn:l**, where the integer l takes the values 0 to
9.

#### nocl1sch

The spline-fit multi-dimensional NOCl S1 potential surface from R. Schinke. The function is in Jacobian coordinates, in the order R, r, theta.

#### lsth

The LSTH potential energy surface of H_{3} for collinear
or general configurations or the 1D potential curve of
H_{2}. For more details on the LSTH surface see
D.G.Truhlar and C.J.Horowitz, JCP **68**, 2466
(1978) and JCP **71**, 1541
(1979). The 1D H_{2} potential is a 87 node spline
fit of Truhlar and Horowitz to the data of W.Kolos and
L.Wolniewicz JCP **43**, 2429 (1965). The
keyword (operator) v:H2 refers to this H_{2} potential
curve.

The following options are available:

- tfac = A . (tfac=0.5 is default)
- jacobian (default) (Coordinates: R, r, θ or R, r)
- binding (Coordinates: r1, r2, φ or r1, r2)
- collinear
- hyper2D (Hyperspherical Coordinates: hyper-radius and hyper-angle)
- hyper3D (Hyperspherical Coordinates: hyper-radius, hyper-angle, and bond angle)

#### tully

CO - Cu(surface) potential of Tully et al., J. Vac. Sci.
Technol. A **11(4)**, (1993), 1914-1920

The follwoing options are available:

- nxy = A (default=1): the number of Cu atoms in the x and y axes
- nz = B (default=1): how many layers of Cu atoms will be taken. Note: The number of Cu-atoms taken into account is nz*(2*nxy-1)^2.
- 2D : a 2 dimensional potential will be applied, using z,r
- 4Da : a 4 dimensional potential will be applied, using x,y,z,r
- 4Db : a 4 dimensional potential will be applied, using z,r,θ,φ
- 6D : a full 6 dimensional potential will be applied, using x,y,z,r,θ,φ (this is the default)

Attention: The degrees of freedom have to be used in this given order.

#### wslfh

The WSLFH potential energy surface for H_{3}O by Wu,
Schatz, Lendvay, Fang and Harding. For details see J. Chem.
Phys. **113**, 3150 (2000).

The following options are available:

- tfac1 = A, tfac2 = B . (tfac1=15.9994/16.9994 and tfac2=0.5 are default)
- jacobian (default)
- binding

The keyword (operator) v:HO used by WSLFH surface is the HO
potential curve taken as a Morse function (for details see J.
Chem. Phys. **113**, 3150 (2000)).

#### yzcl1, yzcl2

The YZCL1 and YZCL2 potential energy surfaces for H_{3}O
by Yang, Zhang, Collins and Lee. For details see J. Chem.
Phys. **115**, 174 (2001).

The following options are available:

- tfac1 = A, tfac2 = B . (tfac1=15.9994/16.9994 and tfac2=0.5 are default)
- jacobian (default)
- binding

Note: for adding this surface use the option "-i yzcl".

#### wdse

The WDSE potential energy surface for H_{3}O by
Walch, Dunning, Schatz and Elgersma, J. Chem. Phys.
**72**, 1303 (1980) and J. Chem. Phys.
**73**, 21 (1980), modified by Echave and
Clary, J. Chem. Phys. **100**, 402 (1994).

The following options are available:

- tfac1 = A, tfac2 = B . (tfac1=15.9994/16.9994 and tfac2=0.5 are default)
- jacobian (default)
- binding

#### pjt2

PJT2 potential energy surface for H_{2}O by Polyansky,
Jensen and Tennyson. See JCP **105**, 6490-6497
(1996) and JCP **101**, 7651
(1994).

The following options are available:

- tfac = A (default: tfac=1.0/16.9994 for jacobian1 and tfac=0.5 for jacobian2)
- binding (r1,r2,theta; default)
- binding2 (r1,r2,cos(theta))
- jacobian1 (H-OH System)
- jacobian2 (O-HH System)

Note: for adding this surface use the option "-i h2o".

#### hfco

HFCO surface from Kato, J.Chem.Phys. **107**
(1997), 6114-6122.

#### gauss1d

One-dimensional Gauss function in coordinate difference. V(x1,x2) = exp(-0.5*((x1-x2)/width)**2)/(sqrt(2*pi)*width)

Usage: gauss1d{width=*w*}, where *w* is to be
replaced by an appropriate real number.

If the two underlying grids are 2pi periodic, one should
additionally set the option *periodic*, i.e.:
gauss1d{width=*w* periodic}.

#### gauss2d

Two-dimensional Gauss function in coordinate difference. V(x1,y1,x2,y2) = exp(-0.5*[((x1-x2)/width)**2 + ((y1-y2)/width)**2])/(2*pi*width**2)

Usage: gauss2d{width=*w*}, where *w* is to be
replaced by an appropriate real number.

#### coulomb1d

One-dimensional regularized Coulomb potential in (scaled) coordinate difference.

```
V(x1,x2) = 1/SQRT( (a*x1-b*x2+c)**2
+ d)
```

Usage: coulomb1d{a=.. b=.. c=.. d=..}. Defaults: a=1, b=1, c=0, d=1.

#### cpp

This multi-dimensional "potential" is needed for, J > 0, kinetic energy operators in e.g. Jacobian or Radau coordinates. Usage: cpp{jtot=J,dim=d}, where J is to be replaced by the value of the total angular momentum and d denotes the number of conjugate momenta.

Definition: cpp(k1,k2,...,kd) = sqrt(J(J+1)-(k1+k2+...+kd)(k1+k2+...+kd+1))

#### cmm

Similar to cpp. Definition: cmm(k1,k2,...,kd) = sqrt(J(J+1)-(k1+k2+...+kd)(k1+k2+...+kd-1))

#### zundel

15D potential of the Zundel cation
(H_{5}O_{2}^{+})

Usage: `zundel{m=.. d0=.. jacobi|valence [ztrafo]}`

Where` m=.. `specifies the mass of the hydrogen atoms of the
water fragments in au (default: 1.0, only used for 'jacobi'.)

`d0=.. `is the value of d0 in the ztrafo below (default: 1.5)

`jacobi: `use Jacobi coordinates **or**

`valence: `use valence coordinates (default)

`ztrafo: `transform z -> z/(R-2d0) (see d-option above, default: not set).

#### zundel_dip

15D dipole surface of the Zundel cation
(H_{5}O_{2}^{+})

Usage: like `zundel` above, however, the component 'x','y', or 'z'
must be given as well.