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.
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. |
| LABELS | Definition of operators needed. |
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.
| 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 = 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 ]where FUNCTION is one of the keywords:
Example
par1 = SIN[ degree*PI/180.0 ] par2 = MAX[ alpha, 3.0*beta-4.5 ] enatural = EXP[ 1.0 ]
| 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 |
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. |
| usepack / nopack |
This keyword controls the expansion of operators into
multiple packages. Default: nopack.
If usepack is set, packages and electronic states are treated with a combined index. If the system contains N electronic states and P packages, a combined index of the N*P internal states is used where the index of the electronic states defined in the PRIMITIVE-BASIS-SECTION runs fastest. (The nth state of the pth packet will have the combined index N*(p-1) + n) If no electronic basis is present, in the PRIMITIVE-BASIS-SECTION the internal states are labeled with an extra operator mode packets. Note: This only works if both, packages and states (if present), are treated with the multi-set formalism. |
| 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
| modes | S1 | S2 | .. |
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
| A | O1 | O2 | O3 | ...... |
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
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 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) and potentials in CPD/CANDECOMP (cpd-file)
format are special. They must not be
multiplied with another operator, e.g.: 1.0 | q*V
is invalid, if V is a natpot or cpdfile. However,
-0.1*I | V is a valid construct, i.e. a
natpot/cpd 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/cpd.
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 ...
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.
| 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. |
| read {file, n} | n th Operator is read from file. 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. |
| srffile{file,path} | A file is read defining a potential energy surface. This is contained in the path/file.srf, and is simply the HAMILTONIAN-SECTION from a .op file. If the keyword "oppath" is given for the path, the file is in the same directory as the .op file. If "default" is given the default operator path is used. |
| 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. |
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.
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-sectionwhere 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.
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.
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 + H2 3D (J=0) in Jacobi coordinates
h3kleg.op: H + D2 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: CO2 vibrational (J = 0) Hamiltonian, 3D, valence coordinates, PES taken from J. Zuniga et al., J. Mol. Spec. 195, 137 (1999).
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
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)
Usage: like zundel above, however, the component 'x','y', or 'z'
must be given as well.