MCTDH: Hamiltonian Documentation

• General Structure
• Sections
• Keywords
• Op_define- Section
• Parameter- Section
• Hamiltonian- Section
• Dissipative- Section
• Labels- Section
• Time-dependent Operators
• Examples
• Available Operators
• Available Surfaces

Back to Main Menu

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

or

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:

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.

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:

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 |j&k construct. In this case one must set ignore in the natpot/cpdfile statement of the LABELS-SECTION. In contrast, mlpfs don't support the ignore keyword; nevertheless one can use the |j&k   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.

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:

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.

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:

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:

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. xi Vi . The data is quadratically interpolated. Blank lines and lines beginning with a hash '#' are ignored. The 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

Of course, multiplication with a coefficient is allowed as well. However, multiplying an mlpf with natpots or muld potentials is not possible, but one can work around this issue by converting the natpot to an mlpf (even if it covers only two DOFs/modes). Also, multiplying an mlpf with operators acting on one of the mlpf DOFS is not supported (yet). The expansion of the mlpfs into mlops is protocolled in the op.log file.

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

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₃ 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:

tfac = A . (tfac=0.5 is default)
jacobian (default)
binding
collinear

c2h

The C₂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₂ + H₂ potential energy surface of Aguado et al. J. Chem. Phys. 101, 4004, (1994). Jacobian coordinates: R, r₁, r₂, θ₁, θ₂, φ.

h4bmkp

The H₂ + H₂ potential energy surface of Boothroyd et al. J. Chem. Phys. 116, 666, (2002).

Jacobian coordinates: R, r₁, r₂, θ₁, θ₂, φ.

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₁, θ₁, r₂, θ₂, φ.)
• cartesian
• interaction

  This only calculates the interaction potential, i.e. the contribution of two isolated H₂ is subtracted:

  V(R,r₁,r₂,θ₁,θ ₂,φ) = BMKP(R,r₁,r₂,θ₁,θ ₂,φ) − BMKP(∞,r₁,r₀,π/2,π/2,0) − BMKP(∞,r₀,r₂,π/2,π/2,0) + BMKP(∞,r₀,r₀,π/2,π/2,0)

  where r₀ = 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, θ₁, θ₂, φ.

Options:

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

h4dj

The H₂ + H₂ potential energy surface of Diep and Johnson. See J.Chem.Phys. 112, 4465 (2000).

This is a 4D surface as the H₂ are treated as rigid rotors (with bond length = 1.449 a.u.).

Coordinates: R, θ₁, θ₂, φ

hoosrf

HO₂ 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₃ for collinear or general configurations or the 1D potential curve of H₂. 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₂ 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₂ 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₃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₃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₃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₂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₅O₂⁺)

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₅O₂⁺)

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