Hamiltonian Documentation

General Structure

The operator file ends with the keyword
end-operator or END-INPUT

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 various sections are listed below. There is no predefined order for the sections.

Keywords

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.
keyword = S , R
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

Parameter-Section

parameter = value , unit

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,

lambda = 0.1825 , ev

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:
a = 2.0,ev     b = 5.3,cm-1     c = a/b
Parameters should not be reused, i.e. expressions like
c = c/b
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 ]

Example

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

Special Parameters

Changing Parameter Values

• 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
  

Hamiltonian-Section

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

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".

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   Build-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

A   |   O1   |   1   |   O3   | .....
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   |&   O1   |   1   |   O3   | .....
A   |1&2   O1   |3   O3   | .....

    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
  

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 Build-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), potentials in CPD/CANDECOMP (cpd-file) format and Multi-Layer Potfits (mlpfs) are special. They must not be multiplied with another operator, e.g.:   1.0 | q*V    is invalid, if V is a natpot/cpd/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.

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

  d/dx * v(x,y) * d/dy

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

  modes |   x   |   y   | ...
  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/cpds/mlpfs.

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 ...

Labels-Section

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
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.

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),g1=1,gdim)
((dimag(oper(g,g1,2)),g=1,gdim),g1=1,gdim)

[ if type=2 then ]
(dble(oper(g,g)),g=1,gdim(1))

[ if type=3 then ]
(( dble(oper(g,g1,1)),g=1,gdim),g1=1,gdim)

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

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.

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]

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

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

Time-dependent Operators

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

Examples

Simple Hamiltonian: The Henon-Heiles System

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

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₂ 3D (J=0) in Jacobi coordinates

h3kleg.op: H + D₂ 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 in the (supplementary) addsrf directory. Supplementary material can be obtained either from our SVN repository or from the packages website.

co2.op: CO₂ vibrational (J = 0) Hamiltonian, 3D, valence coordinates, PES taken from J. Zuniga et al., J. Mol. Spec. 195, 137 (1999).

Complete List of Operators

Available Surfaces

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

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

readsrf