Input Documentation

General Structure

• Special characters
• Using the C preprocessor (CPP) with MCTDH
• Keyword placement
• Arguments to keywords
• Units

Special characters

Using the C preprocessor (CPP) with MCTDH

As an example, take a look at the input files allene_e.inp and allene_b.inp with common settings in the include files allene_common1.inc and allene_common2.inc.
These examples are also listed in the examples section below.

Keyword placement

The order and placement of the keywords within a section is usually arbitrary. Exceptions are described below.

Arguments to keywords.

    rd    sin     36   3.800    5.60

CAP_rd = CAP [ 5.0  0.357    3 ]

Units.

SECTIONS

The input sections reflect the structure of the program. The RUN-SECTION defines what sort of calculation is desired. Depending on what is requested here the remaining sections provide the required information. The various sections are listed below. There is no pre-defined order for the sections.
 

Below are tables of keywords for each input section. 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. Default values for keywords and arguments are also listed.

RUN-SECTION

There are 4 levels of calculation types, reflecting the four stages of a calculation. Only one run-type keyword of a particular level is allowed. All lower level keywords are automatically included. Thus
propagation
or
gendvr genoper geninwf propagation
are equivalent.

When relaxation=I is given, the I-th eigenstate (counting from 0) is computed by taking the I-th eigenstate (A-vector) of a Lanczos matrix, rather than propagating the A-vector in imaginary time. (The single-particle functions are still propagated in imaginary time.) This requires, that one uses the CMF integrator in the fix or varphi mode and employs the SIL integrator for the A-vector. (NB: The keyword CMF is interpreted as CMF/varphi when the runtype is improved relaxation. Otherwise CMF is always interpreted as CMF/var, see INTEGRATOR-SECTION). It also requires, that the Hamiltonian is hermitian such that the SIL-integrator and not the Arnoldi-Lanczos (CSIL) integrator is taken. If I is replaced by the string follow then the eigenvector with the largest overlap with the restart vector is computed. (The restart WF may be generated by build or taken from a file. See INIT_WF-SECTION).
The size of the Lanczos matrix used depends on the SIL order and accuracy parameters. The Lanczos iteration is stopped when the SIL order is reached or when the estimated error of the Lanczos eigenvalue (in milli Hartree) is smaller than the SIL accuracy parameter. When the keyword full is given, the very first Lanczos space is build up to maximal order, irrespectively of the accuracy. This feature is useful to ensure, that one starts with an appropriate approximation of a desired excited state. When the keyword ortho is given, a full re-orthonormalisation of the Lanczos vectors (or Krylov vectors) is performed. In general it is recommended to use ortho.
The keyword ortho must not be given when the Davidson integrator DAV is used.
When the chosen integrator mode is CMF/varphi (or just CMF ), then the output starts at t=tinit-tpsi. The initial wavefunction is reported at this time. At t=tinit (i.e. usually at t=0) we report the wavefunction obtained after diagonalising the Hamiltonian in the basis of the single-particle functions (without having relaxed these functions). The next steps include first a relaxation of the single-particle functions followed by a (Lanczos) diagonalisation of the Hamiltonian.
The calculation stops, when the final relaxation time is reached or when the A-vector changes less than the SIL-accuracy. The program then stops with the message stop by internal check. Note, that "relaxation" to exited states may require a rather large Lanczos order. See the steps file and in particular the rlx_info file for details on the computation.
When the calculation starts converging to the desired state but after a few iterations jumps to another state, then the numbers of single-particle functions are too small or the chosen accuracies are to low. Increase the accuracy of CMF/varphi (or decrease the step size of CMF/fix). It may also be necessary to increase the Lanczos accuracy and/or order. Inspect the rlx_info file to see what Lanczos orders are actually taken. Also inspect the update file. If a "*" appears at the beginning of a line, then the corresponding update time was too large.
The inputs relaxation and relaxation=0 both generate the ground-state wavefunction. However, different algorithms are used and the latter calculation will in general be considerably faster. NB: Only the simple relaxation works for density operators, relaxation=I will not work for density operators.
For an example see the input files co2_gs.inp and co2_asym.inp which generate the ground state and the asymmetric stretch excited state, respectively. See also the User's Guide, section 3.4 .

The recently introduced Davidson "integrator" DAV is much more powerful for improved relaxation than the SIL. For DAV one of the three inputs: relaxation= I or relaxation=follow or relaxation=lock may be used. For the first input, the program seeks the ground-state (I=0) or for the I-th excited state. However, the program takes the I-th state of the first diagonalisation, which may not be the I-th state of the Hamiltonian. With the option full the program tries to converge as many states as it can within the limited space supplied by the "integrator". full thus makes the counting safer, but the computation more costly. With follow the Davidson diagonalisation takes that eigenvector which is closest to the one of the previous diagonalisation. With lock the Davidson diagonalisation takes that eigenvector which is closest to the initial WF as defined in the Init_WF-Section. lock is hence the safer choice for finding excited states, but it is numerically more demanding than follow. lock requires that the keyword cross is given in the Run-Section.

There are three Davidson-routines implemented: DAV, RDAV, and CDAV. (See INTEGRATOR-SECTION). The RDAV routine allows to additionally use the keywords quad, olsen, backrotate, quadphi, and cn, (where n denotes an integer, e.g. c6) as arguments to relaxation .
The keyword quad lets the program switch to use a quadratic variational principle (i.e. (H-E)**2 is diagonalised within the space of Davidson-vectors rather than H).
With the keyword olsen the program applies the Olsen correction to the residual Davidson vector. This is only useful when a preconditioner is used. (Keyword precon).
The keyword backrotate lets the program rotate the SPFs after relaxation, such that they have maximal overlap with the orbitals prior to the orbital relaxation step. This sometimes improves the overlap with the previous A-vector.
The keyword quadphi lets the program use a different propagation algorithm for orbital relaxation (experimental).
The keyword cn lets the program perform n orbital relaxation steps, before a new diagonalisation step is done.

The file rlx_info contains a lot of information on the relaxation process. If DAV / RDAV / CDAV is employed, use the script rdrlx to read it. With the aid of the keyword rlxunit one may set the energy-unit for the output to the rlx_info file. One also may apply an energy-shift (e.g. subtract the ground-state energy). See the table Keywords associated with a propagation or relaxation calculation below.

The exact keyword can be added to any run-type keyword. For example,
propagation exact
will result in a numerically exact wavepacket propagation (i.e. the wavefunction will be represented in the full product primitive basis). Similarly,
geninwf exact
will set up a wavepacket for a numerically exact calculation.

These keywords are equivalent to a run-type keyword of the level given.

The tfinal keyword must be given. All other keywords are optional. The following default values are used.

The following keywords define the data calculated and saved. If keywords are omitted, data will not be calculated.

OPERATOR-SECTION

If no OPERATOR-SECTION is included, the program looks for the operator information in the input file. This can be useful if a small model operator is studied. A full log of the operator is then automatically output.

SPF-BASIS-SECTION

The following lines define the single-particle function basis to be used in a wave function calculation or a calculation employing density operators of type II. The input defines firstly how many degrees of freedom are contained within a mode. Secondly, the number of spfs are given; a list being needed for a multi-set basis. The format is:

mode_label1 , mode_label2 , ... = nspf1 , nspf2 , ...

For example for a 3-mode system, with labels X, Y and Z,

to define an spf basis of 3 functions per mode,

spf-basis-section
    X = 3
    Y = 3
    Z = 3
end-spf-basis-section

spf-basis-section
    X = 3    Y = 3    Z = 3
end-spf-basis

spf-basis
    X, Y = 3    Z = 3
end-spf-basis-section

spf-basis-section
multi-set
    X , Y = 3 , 2
    Z = 3,2
end-spf-basis-section

Note: The electronic SPF-Basis is not to be specified in the spf-basis-section, as it is always complete. The electronic mode will always be the last mode in a single-set run.

If many electronic states are present, the definition of single-particle functions for a mode can be continued on a second line by using a continuation mark %, e.g.

spf-basis-section
    X    =   5 , 5 , 5 , 5 ,  6 , 10 , 5 , 2 , 2 , %
             5 , 5
end-spf-basis-section

If for symmetry reasons one set of SPFs is always identical to another one, the latter set need not to be propagated numerically. In such a situation, e.g. H₂O in valence coordinates and for a symmetric initial state, one may tell the program not to propagate the second identical set of SPFs. This is done by specifying with the id keyword the mode to which the present mode is identical. E.g.:

spf-basis-section
    R1    = 8
    R2    = id,1
    Theta = 9
end-spf-basis-section

The following keywords, if given, define multi-state or muti-packet runs. Note: A SPF-BASIS-SECTION does not need to exist for an exact calculation, except if packets is specified.

SPDO-BASIS-SECTION

The SPDO-BASIS-SECTION is a special feature of a calculation using density operators of type I. It is organised almost as the SPF-BASIS-SECTION. For a single-set calculation these sections are in fact identical. In a multi-set calculation differences occur due to the special structure of density operators of type I.

The number of sets in a multi-set calculation using a density operator of type I is the squared number of sets used in the corresponding calculation emplyoing a wave function or a density operator of type II. This can be seen in the following example:

If the SPF-BASIS-SECTION reads (two sets for each DOF)

spf-basis-section
multi-set
    X , Y = 3, 2
    Z     = 4, 2
end-spf-basis-section

the SPDO-BASIS-SECTION may be chosen as (four sets for each DOF)

spdo-basis-section
multi-set
    X , Y = 2, 2, 2, 2
    Z     = 2, 2, 2, 2
end-spdo-basis-section

The ordering is here as
n(1,1), n(1,2), n(1,3), ..., n(2,1), n(2,2), n(3,2), ...
where n(s,t) is the number of SPDOs for the pair of states (s,t).

PRIMITIVE-BASIS-SECTION

mode_label   basis_type   basis_size    parameters

mode_label is an alphanumeric string (case sensitive) labelling the degree of freedom.

basis_type must be one of the following:
 

basis_size is an integer specifying the primitive basis size, e.g. grid points or vector elements etc. Note that for an FFT-representation basis_size (i.e. gdim in the program) must have a prime factor decomposition with only 2's, 3's and 5's but should have a decomposition with only 2's and 3's for optimal performance, i. e. basis_size = 2^m * 3^n where m and n are positive integers. Note also that basis_size must be odd for the exp-DVR.
For a sphFBR-representation, basis_size is not required in input. The number of basis functions is calculated by the program itself, according to the type of basis (see parameter description).
For a K-DVR basis_size is not required in input. It will be calculated from the kmin,kmax parameters given in this section. NB: The basis_size is called gdim.

The parameters to be input depend on the basis type as follows:

Remarks on continuum electronic basis:

    el     elcont      4    2
    Elcont     sin  101  0.0, ev  3.75, ev
  

   el        gauss   1.0, ev  0.0  0.3, ev

Remarks on harmonic oscillator DVR:

The HO-DVR depends only on the product hofreq*homass. If the homass entry is missing, the program sets homass to 1. Alternatively, one may specify the first and last grid-point. The program then computes the corresponding product hofreq*homass.
Example: The following lines are equivalent.

    Y    HO    32    0.00       0.10        1822.89
    Y    HO    32    0.00       2.721,eV    1822.89,au
    Y    HO    32    0.00       0.1,au      1.0,AMU
    Y    HO    32    0.00   21947.46,cm-1   1.0,AMU
    Y    HO    32    xi-xf     -0.528       0.528

Remarks on Laguerre DVR:

The Lagua-DVR is build from the basis functions:
phi(n,x) = Sqrt((n-1)!/(n+a-1)!) * x^(a/2) * exp(-x/2) * L(n-1,a,x) ,
where a = 1,2,3 or 4 (Lagu1 - Lagu4). Hence, the boundary condition for x -> 0 is phi(x) ~ x^(a/2). (With the aid of the parameters x0 and b the coordinates may be shifted and scaled, i.e. phi(x) <- phi((x-x0)/b)) ). The distribution of the grid points is very uneven, being very dense for small x and widely spaced for large x. The matrix of second derivatives may have very large negative eigenvalues. These will slow down the integrator. The integer parameter icut helps to fix this problem. The matrix of second derivatives is diagonalized and the first icut eigenvalues (these are the largest negative eigenvalues) are replaced by the icut+1st one. The thus modified eigenvalues and the eigenvectors are then used to build the working matrix of second derivatives. Note that the FBR matrix representation of { d^2/dx^2 - c/x^2) } is analytically evaluated ( c = -1/4, 0, 3/4, 2 for a = 1, 2, 3, 4). After this matrix is transformed to the DVR representation, the centrifugal term c/x^2 is removed by substraction. The width parameter b has to be chosen carefully. Its optimal numerical value will depend on N, the number of grid points. Alternatively, one may use the input formats xi-xf (not recommended in general) or x0-xf. These formats compute b.

Remarks on sine DVR:

short: xi and xf denote, as usual, the first and last grid point. short is default and need not to be given.
long : xi and xf denote the boundaries of the sine basis functions (i.e. boundaries of the 'box') and not the first and last grid point.
Example: The following lines are equivalent.

  x   sin     19    1.00    19.0     short
  x   sin     19    0.00    20.0     long

Remarks on cosine DVR:

The basis functions underlying the DVR are 1/sqrt(L), (2/sqrt(L))*cos[(j*pi/L)*(x-x₀)]. The symmetrized derivative, sdq = 0.5*( sin[(pi/L)*(x-x₀)] * d/dx + d/dx * sin[(pi/L)*(x-x₀)] ) is used as first derivative. Because only gerade (symmetric) wavefunctions are computed, we consider only the interval [x₀,x₀+L], although the wavefunction is periodic on the interval [x₀-L,x₀+L]. xi and xf are the first and last grid-point, but when long is given, the input is interpreted as x₀ and x₀+L.
Example: The following lines are equivalent.

  x   cos     36    2pi/2
  x   cos     36    0.00      3.1415926535897  long
  x   sin     36    0.0436332313  3.097959422  short
  x   sin     36    0.0436332313  3.097959422

Remarks on FFT and exponential DVR:

FFT and exp-DVR have an identical set of input parameters (if exp-DVR is not combined with PLeg-DVR). These two methods are in fact largely equivalent and produce identical results (for same parameters). The numeric, however, is different and the exp-DVR will be faster for small grids whereas FFT is faster for long grids. Around 30 grid points both methods are of similar speed. The use of the interaction-picture is possible for exp-DVR.
FFT and exp-DVR enforce periodic boundary conditions, but they are often used for ordinary coordinates. A set of keywords adapts the input to the various situations:

linear: xi and xf are the coordinates of first and the last point of the grid. The grid-spacing is dx = (xf-xi)/(N-1). Due to the periodic boundary conditions the first grid point and the one following the last grid point are to be identified.

periodic: xi and xf are considered as identical due to the periodic boundary conditions. The grid-spacing is dx = (xf-xi)/N. The routine eingabe.f re-scales xf -> xf-dx.

s-periodic: xi and xf are considered as identical due to the periodic boundary conditions. The grid-spacing is dx = (xf-xi)/N. The grid points, however, are now placed symmetrically on the interval (xi,xf) and eingabe.f re-scales xi -> xi+dx/2, xf -> xf-dx/2.

2pi/mult or s-2pi/mult or c-2pi/mult: The routine eingabe.f sets xi=0 and xf=2*pi/mult and then performs according to the periodic or s-periodic keyword. For c-2pi/mult the grid is shifted such that xi=-xf.
Example: The following sets of lines are equivalent.

  x   fft     32    0.00       3.0434179    linear
  x   fft     32    0.00       3.1415927    periodic
  x   fft     32    2pi/2

or

  x   fft     32    0.0981748  6.1850105    linear
  x   fft     32    0.00       6.2831853    s-periodic
  x   fft     32    s-2pi

or

  x   fft     32   -3.0434179  3.0434179    linear
  x   fft     32   -3.1415927  3.1415927    s-periodic
  x   fft     32    c-2pi

Example:

For a system with two degrees of freedom, labeled X and Y, the following would define for X an FFT grid of 32 points from -2 to 2, and for Y a DVR basis of 32 harmonic oscillator functions generated with the given parameters. To show the use of the el-keyword we assume that there are three electronic states.

primitive-basis-section
  el     el     3
   X    FFT    32    -2.00    2.00    linear
   Y     HO    32     0.00    5.2,eV   1.00
end-primitive-basis-section

Remarks on External DVR:

Extern-DVR is an external DVR, where grid points and DVR derivative matrices are read from a file. The file can be read in ascii or binary format:

in ascii format (default):

do i=1,gdim
   read(unit,*) ort(g)
enddo
do i=1,gdim
   read(unit,*) (dif2mat(j,i),j=1,gdim)
enddo
do i=1,gdim
   read(unit,*) (dif1mat(j,i),j=1,gdim)
enddo

do i=1,gdim
   read(unit) ort(g)
enddo
do i=1,gdim
   read(unit) (dif2mat(j,i),j=1,gdim)
enddo
do i=1,gdim
   read(unit) (dif1mat(j,i),j=1,gdim)
enddo

where gdim is the basis_size. The file can have absolute or relative path. If file containes only grid points (for example in POTFIT program, when only grid points are used), the DVR matrices will be zeroed. If only second derivative matrix dif2mat is given, dif1mat will be zeroed.

Example:

     x   extern  30   x_data
     y   extern  50   y_data     binary
     z   extern  42   z_data     binary  angst
 theta   extern  23   theta_data deg

Remarks on KLeg, PLeg and sphFBR:

KLeg, PLeg and sphFBR all define 2D single-particle functions, although the basis definition is for each degree of freedom individually. KLeg must hence be followed by K and PLeg by exp and sphFBR by phiFBR.

Example:

PRIMITIVE-BASIS-SECTION
  alpha   sphFBR   30   sym
  beta    phiFBR    5
  theta1  PLeg     31   even
  phi     exp      15   2pi   k=-6,6
  theta2  KLeg     31   even
  K_th     K       -5    5
end-primitive-basis-section

The KLeg/K, PLeg/exp and sphFBR/phiFB combinations generate mode-operators. Hence the DOFs (theta2, K_th), (theta1,phi), and (alpha,beta) must be combined with each other and must not be combined with any other mode. This excludes the use of these DVR's when the WF is expanded in exact format, because the exact format combines all modes into one particle. Note that the exact format is used by a diagonalisation run.

The above restriction has been relaxed somewhat in recent versions of MCTDH. Since version 8.3.13, you can combine several KLeg/K into one mode, even along with other DOFs. This makes it also possible to use KLegs in exact calculations. However, since this feature is rather new, you are advised to check your results carefully if you use it.

INIT_WF-SECTION

After some initial WF is generated, either by file,or build, or read-inwf, this WF may be modified by one or several of the following procedures. The program will take actions in the order:
     A-coeff,  meigenf,   operate,    orthogonalisation,    correction,
irrespectively of the order in the input file.

The INIT_WF-SECTION is used also for generating an initial density operator (of type I and II). For the generation of a pure state the initial wave function ket is multiplied with the corresponding bra, i.e.

rho = | Psi > < Psi |

Hence, in this case it is necessary only to specify an initial wave function. To obtain a thermalised initial state the temperature key word (see below) has to be used. The initial wave function is then assumed to represent the ground state, and the parameter values are taken to generate the appropriate thermalised state.

S-MCTDH, selected CI

The selection of configurations is controlled by the cut-off parameter which is given as an argument to the keyword s-mctdh. The method is described in G. A. Worth, J. Chem. Phys. 112, 8322 (2000). The default value for the selection parameter is 2.0. Smaller values may lead to rather inaccurate results and larger values may not lead to a speed-up. In general, S-MCTDH will only be useful, if there are several particles (5 or more) and if the propagation of the A-vector is much more costly than the propagation of the SPFs. S-MCTDH is incompatible with operate and cross, requires that the wavefunction is build and not read in from a file, and requires CMF propagation with natorb. Finally, several of the analyse routines are not able to handle S-MCTDH wavefunctions.

Building the initial wavefunction

The information needed to generate the initial wavefunction is written one line per degree of freedom between the keywords build and end-build. The input is free format, with blanks dividing the various parameters. The format for each line is

modelabel    type    parameters

The modelabel is an alphanumeric string attached to the degree of freedom. This must match a label specified in the primitive-basis section.

If one uses an electronic basis (i.e. one degree of freedom has the primitive basis type el), then the electronic initial state can be specified by the init_state keyword,

init_state = s

where the integer s specifies the initial state. This statement must be the only one on a line. If the lowest electronic state is selected, this keyword is not required.

If one uses electronic states in a density operator propagation the keywords corresponding to init_state are left_state and right_state (for both types). The usage is

left_state = s
right_state = t

where the integers s,t specify the initial states of the density operator. In a ket-bra notation this means |s><<I>t|. Each statement must be the only one on a line. If s or t denotes the lowest electronic state, the corresponding keyword is not required.

To obtain a thermalised initial state for a density-operator propagation the keyword

temperature = T

has to be used where T is the temperature. For density operators of type I an analytical formula is applied to generate the thermalised initial state. Here the size of the primitive grid has to be chosen large enough to obtain an accurate representation. For density operators of type II the number of SPFs is crucial for an accurate representation.

To summarize:

When building the initial wavefunction, the type of the 1D functions can be chosen from one of the following:

The parameters needed for each type of function are as follows:

The initial single-particle functions are formed as follows:

• If there is an electronic degree of freedom treated in the single-set formulation, the single-particle functions are vectors representing the population of the states. The first vector is populated in the state specified by init_state, while the remaining vectors are generated orthonormal to this vector and fully populated in the remaining states. For example, if there are three states and the second state is initially populated, three vectors are produced, (0 1 0), (0 0 1) and (1 0 0), respectively.
• For the type gauss the first single particle function is a normalised Gaussian given an initial momentum momentum of the form f(x) = N exp(-1/4 ((x-centre)/width)²) exp(I momentum (x-centre)). Note: width may be chosen to be complex! (In this case the input reads "Re,Im,unit", e.g. "1.2,0.5,eV"). Higher single-particle functions are produced by multiplying the preceding function by x (so producing a series of powers of x), followed by Schmidt-orthogonalisation onto the lower functions.
• For the type HO the first single-particle function is a normalised harmonic oscillator of the form f(x) = N exp(-1/2 frequency mass (x-centre)²) exp(I momentum (x-centre)). Note: frequency may be chosen to be complex! (In this case the input reads "Re,Im,unit", e.g. "1.2,0.5,eV"). Higher single-particle functions are produced by multiplying the preceding function by x (so producing a series of powers of x), followed by Schmidt-orthogonalisation onto the lower functions.
• For the type HO odd the first single-particle function is a normalised harmonic oscillator of the form f(x) = N x exp(-1/2 frequency mass (x-centre)²) exp(I momentum (x-centre)). Note: frequency may be chosen to be complex! (In this case the input reads "Re,Im,unit", e.g. "1.2,0.5,eV"). Higher single-particle functions are produced by multiplying the preceding function by x^2, followed by Schmidt-orthogonalisation onto the lower functions.
• For the type HO even the first single-particle function is a normalised harmonic oscillator of the form f(x) = N exp(-1/2 frequency mass (x-centre)²) exp(I momentum (x-centre)). Note: frequency may be chosen to be complex! (In this case the input reads "Re,Im,unit", e.g. "1.2,0.5,eV"). Higher single-particle functions are produced by multiplying the preceding function by x^2, followed by Schmidt-orthogonalisation onto the lower functions.
• If the type Leg is chosen the single-particle functions are normalised (associated) Legendre polynomials. If the corresponding primitive basis is not Leg, then m must be zero. The parameter l defines which Legendre polynomial is taken as the first single-particle function. If the parameter symmetry equals nosym then all (even and odd l) Legendre polynomials are taken, and if sym, then only even (if l is even) or odd(if l is odd) Legendre polynomials are taken. Note the following restrictions: l >= m; if the corresponding primitive basis is Leg, then blz = l and ibleg = symmetry as defined in primitive-basis-section.
• In the sphFBR case the first single-particle function is a normalised spherical harmonic |j,m> = |j,m>. Higher single-particle functions are other spherical harmonics whose energies are as close as possible from the energy of the first one. Of course, the initial spherical harmonic must be part of the basis set.
• If eigenf is given, the single-particle functions are eigenfunctions of a one-dimensional Hamiltonian defined in the HAMILTONIAN-SECTION_XX of the operator file. The keyword pop selects which eigenfunction becomes the initial single-particle function. The ground state corresponds to pop=1, the first excited state to pop=2, and so on. If pop is set to a negative integer, only every second eigenfunction will be taken as initial single-particle function. This is useful, if the one-dimensional Hamiltonian to be diagonalised is symmetric. Thus, pop=-1 selects gerade eigenstates, and pop=-2 the ungerade ones. The selection (together with the computed eigenvalues) is protocoled in the log-file. If you want the eigenvectors too, specify "veigen" in the RUN-SECTION.
  Operators defined in a HAMILTONIAN-SECTION_XX are usually multi-dimensional operators. As a one-dimensional operator is required for diagonalisation, the program takes the uncorrelated part of the desired degree of freedom from the specified Hamiltonian. This is most conveniently described by an example. Assume that there are the DOFs x,y,z, and el. The eigenf-input reads:

   y  eigenf Hspf pop=2
  

  Then each of the following three Hamiltonians can be used equally well to generate the desired initial single-particle functions.

  HAMILTONIAN-SECTION_Hspf
  modes    |  y
  1.0      |  KE
  omy      |  q^2
  end-hamiltonian-section
  

  HAMILTONIAN-SECTION_Hspf
  ------------------------------------
  modes    |  x   |  y   |  z  |  el
  ------------------------------------
  1.0      |  1   |  KE  |  1  |  1
  omy      |  1   |  q^2 |  1  |  1
  ------------------------------------
  end-hamiltonian-section
  

  HAMILTONIAN-SECTION_Hspf
  ------------------------------
  modes    |  x   |  y   |  z
  ------------------------------
  1.0      |  KE  |  1   |  1
  omx      |  q^2 |  1   |  1
  1.0      |  1   |  KE  |  1
  omy      |  1   |  q^2 |  1
  1.0      |  1   |  1   |  KE
  omz      |  1   |  1   |  q^2
  const    | q^3  | sin  |  q
  ------------------------------
  end-hamiltonian-section
  

  The lines 1,2 and 5,6 of the tableau of the third Hamiltonian are ignored, as they refer to other degrees of freedom than the y one. The last line is ignored, as it represents a correlated Hamiltonian term.
  NB: There is no longer a need to set S1&1 for the electronic DOF, as it was the case in older versions. The old operator files, however, still work.
• If KLeg is chosen, the single-particle functions are normalised associated Legendre polynomials. The corresponding primitive basis must be KLeg or PLeg. The input line that follows must be for the type K. The parameter l defines first index (orbital angular momentum) of associated Legendre polynomial taken as the first single-particle function. The second (magnetic) index is given in the next K initial-wavefunction input-line. If the parameter symmetry equals nosym then all (even and odd l) Legendre polynomials are taken, and if sym, then only even (if l is even) or odd (if l is odd) Legendre polynomials are taken.
• The type K may appear after Leg, KLeg or gauss only. If the previous initial wavefunction is of type Leg or KLeg, the parameters here define the second (magnetic) index of the associated Legendre polynomials defined in the Leg or KLeg initial wavefunction. Then k is the initial value and kmin,kmax,dk are the bounds and the step of k. If the previous initial wavefunction is of type gauss then k indicates the initially occupied value of k-quantum number. In this case the last three parameters must be omitted. Note: The type K requires that the primitive basis must be KLeg and K or or PLeg and exp.
  The following example (inelastic N₂/LiF surface scattering) exemplifies the use of KLeg and K.

  PRIMITIVE-BASIS-SECTION
  ---------------------------------------------------
  # Mode    DVR      N      xi       xf
  ---------------------------------------------------
     z      fft     96   -0.75   4.5000  linear
     x      fft     24    0.00   5.3669  periodic
     y      fft     24    0.00   5.3669  periodic
     theta  PLeg    36    even
     phi    exp      9    2pi
  ---------------------------------------------------
  END-PRIMITIVE-BASIS-SECTION
  
  INIT_WF-SECTION
  build
  -----------------------------------------------------------
  # mode    type center  moment.  freq.         mass
  -----------------------------------------------------------
      z      HO   1.80   -24.d0 694.1683673,eV  1.d0
      x      HO   0.00     0.0    0.0           1.d0  # flat function
      y      HO   0.00     0.0    0.0           1.d0  # flat function
   theta     KLeg  2       sym      # j of initial spherical harmonic
     phi      K    1                # m of initial spherical harmonic
  -----------------------------------------------------------
  end-build
  end-init_wf-section
  

• Mapping of initial (1D) single-particle functions is particularly useful when the initial SPF shall be generated by eigenf but an FFT is desired for the propagation step. The dummy grid may have fewer grid-points than the propagation grid, i.e the grid to which the SPF is mapped. However, the grid points of the dummy grid must coincide with those of the propagation grid. This, in practice, limits mapping to equidistant grids, i.e. FFT, exp-DVR and sin-DVR. The following example may be useful. Here the dummy grid coincides with the grid-points 6-65 of the propagation grid.

  PRIMITIVE-BASIS-SECTION
             .
             .
    R       FFT    192   -0.6          3.0
    R_dum   sin     60   -0.5057591623 0.6062827225
  end-primitive-basis-section
  
  INIT_WF-SECTION
  build
             .
             .
    R      map     R_dum
    R_dum  eigenf  operator    pop=1
  end-build
  end-init_wf-section
  

  Since the DOF   R_dum (in our example) is not defined in the SPF-BASIS-SECTION it will be ignored, when the operator is build. To avoid this, one has to use the addmode keyword.

  HAMILTONIAN-SECTION_operator
  addmode = R_dum
  ----------------------------------------
      modes     ......     | R_dum  | ....
  ----------------------------------------
             .
             .
  

For multi-set calculations, the initial functions can be defined differently for each state by using the state keyword:

modelabel    type    parameters    state = s

For multi-packet calculations, i.e. when packets > 1, the initial function definition must be given for each packet as

modelabel    type    parameters    pack = p

By default the 1st single particle function of each degree of freedom is used to form the initial Hartree product wavefunction. The initially populated single-particle function can be selected using the pop keyword:

modelabel    type    parameters    pop = i

Reading the initial wavefunction

Examples:
For a two state WF the following input is equivalent to a simple file = <restart-directory> statement.

Read-Inwf
  orthopsi  # This is default, opposite: noorthopsi
  file = <restart-directory>
  SPF 1 -> 1
  SPF 2 -> 2
  A   1 -> 1
  A   2 -> 2
end-read-inwf

Read-Inwf
  file = gs
  SPF 1 -> 1
  SPF 1 -> 2
  SPF 1 -> 3
  A   1 -> 2
end-read-inwf

Read-Inwf
  file = gs
  SPF 1 -> 1,2,3
  A   1 -> 2
end-read-inwf

Read-Inwf
  file = run1
  SPF 1 -> 1
  file = run2
  SPF 1 -> 2
  A   1 -> 1
  file = run3
  A   2 -> 2
end-read-inwf

Read-Inwf
  file = run1
  SPF 1 -> 2
end-read-inwf

INTEGRATOR-SECTION

The string S=all after the integrator name may be omitted, i.e. ABM is equivalent to ABM/ALL. For VMF calculations, only the BS, ABM or RK5/8 integrator may be used for the differential equations. Default is ABM. For VMF calculations the integrator must carry the extension S=all (or no extension at all), i.e. there is only one integrator within the VMF scheme. For CMF calculations, the following combinations of integrators are possible: ABM/spf + SIL/A, BS/spf + SIL/A, RKx/spf + SIL/A, BS/all, ABM/all, RKx/all. Default is BS/spf + SIL/A. For numerically exact calculations (i.e. the keyword exact has been specified in the RUN-SECTION), any of the integrators may be chosen, with S = all, or S = spf. Default is ABM, but more efficient is usually SIL (or CSIL).

The Davidson "integrator", DAV, (actually a diagonaliser, of course) is the optimal choice for the "A-propagation" of improved relaxation. The accuracy parameter is an upper limit (in au) for the error of the eigenvalue. Improved relaxation requires CMF/varphi or CMF/fix. (Simply use CMF). The routine DAV is for hermitian Hamiltonians and general wavefunctions. rDAV is for real Hamiltonians (i.e. H*psi = real for all real psi) and real (initial) WF. Memory is saved as the Davidson vectors are stored as real. Most of the arithmetic (in particular H*psi), however, is still complex. This problem is partly solved by rrDAV. rrDAV uses the same Davidson routine, but a simplified routine is employed to perform the matrix product H*A. This routine is faster, because it uses more real arithmetic, but works only for simple Hamiltonians (no electronic states, only real operators, e.g. no p. NB: The propagation of the SPFs is always done in complex arithmetic. cDAV is for non-hermitian Hamiltonians. This allows to compute resonances of CAP augmented Hamiltonians. The convergence, however, is slower. cDAV is chosen automatically, if CAPs are present.

NOTE: There are two Lanczos integrators, a real version for hermitian Hamiltonians and a complex version (Lanczos-Arnoldi) for non-hermitian ones. If one gives the keyword SIL, then the program tries to detect whether the Hamiltonian is hermitian or not and it makes the appropriate choice. This, however, works safely only for CAPs. When one knows that the Hamiltonian is non-hermitian, one should use CSIL. The use of SIL in case of a non-hermitian Hamiltonian may lead to wrong results.

Omitted integrator parameters default to:

Example Inputs for wave functions

• NOCl S0 (long)
• NOCl S0 (short)
• Relaxation on NOCl S0 surface
• Propagation on NOCl S1 surface
• CH₃I (relaxation, single-set)
• CH₃I (propagation, single-set)
• CH₃I (relaxation, multi-set)
• CH₃I (propagation, multi-set)
• H+D₂ reactive scattering (Coupled States Approximation)
• H+D₂ reactive scattering
• Pyrazine (6 modes)
• N₂-LiF surface scattering
• 1D harmonic oscillator (diagonalisation)
• 2D LiCN model (diagonalisation)
• CO₂ vibrational spectrum
• Benzene photoelectron spectrum (continuum explicitely included)
• Further example inputs

Example Inputs for use of CPP with MCTDH

• Allene E state vibrational spectrum
• Allene B state vibrational spectrum

• Include file #1 with common settings
• Include file #2 with common settings

Example Inputs for density operators of type I

• Open Henon-Heiles system
• Open Pyrazine system (four modes)

Example Inputs for density operators of type II

• Open Henon-Heiles system
• Open Pyrazine system (four modes)
• Closed Pyrazine system (four modes, 300 K)