QUANTICS: Input Documentation

• General Structure
• Units
• Sections
• Run- Section
• Operator- Section
• Dirdyn- Section (DD-vMCG)
• Dirdyn-GB- Section (DD-GB)
• SPF-Basis- Section
• ML-Basis- Section
• SPDO-Basis- Section
• Initial- Geometry- Section
• DD-GB- Section
• Primitive- Basis- Section
• Init_wf- Section
• Integrator- Section
• DDpoints- Section
• Example Inputs

Back to Main Menu

Input Documentation

General Structure

The input file is divided into sections containing information for the various parts of the program. The sections of the input file may appear in arbitrary order. All input is based on keywords with arguments. The keywords are in general not case sensitive. However, certain strings such as mode-labels and file-names are case sensitive!

The input file ends with the keyword

Note, that some keywords may be overwritten by options. See Capabilities and usage/Starting the Program.

The following information is general to the input file.

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

Special characters

All blank lines are ignored. Text after the hash '#' symbol is also ignored. This symbol can thus be used to provide comments, e.g.

or to comment out an entire line.

All keywords that start with an exclamation mark '!' are ignored. Note that only the keyword which carries the '!' is ignored. To ignore 'keyword = arg' one must write '!keyword != !arg'.

Using the C preprocessor (CPP) with QUANTICS

If QUANTICS is invoked with the command line option -cpp, the C preprocessor is run over the input file prior to interpretation by QUANTICS. The output of CPP is saved in a temporary file in the directory where the input file resides. Please make sure that this directory has write permission for QUANTICS. Amongst other features, using CPP allows to use a line of the form

to include the file some-common-settings.inc anywhere in the input file. This way, several input files can share common basis sets, initial states or any other settings.

As CPP directives start with the hash character '#', comments have to be hidden from CPP by using the C construct

The '/*' is interpreted by QUANTICS like '#', i.e. everything which follows until the end of the line is ignored. The end-of-comment sign '*/' is only needed for CPP. Thus a comment /*..*/ must not go over several lines, one has to bracket each line individually.

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

Keywords are free format. Blank spaces, ' ', or semi-colons, ';', can be used to divide the keywords. e.g.

    keyword1
    keyword2
    keyword1 keyword2
    keyword1 ; keyword2
    keyword1;keyword2
    

are all possible.

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

Arguments to keywords.

If a keyword requires arguments, these are normally signified by an '=' symbol, and divided by ',' symbols. Brackets '(' ')' as well as spaces can also be used to help readability, e.g.

    keyword=arg1,arg2,...
    keyword = arg1 , arg2,...
    keyword = (arg1 , arg2,...)
    

are all possible formats.

The arguments are optional in most cases. The equal-sign '=' and/or the comma ',' indicates that the keyword following it is an argument. Each of the three lines given below is a correct example input line:

    output psi timing
    output psi = double timing
    output psi = double, natur timing
    output psi = (double, natur) timing
    

since the keyword psi knows the arguments double and natur. However

    output psi = double, timing
    

is an incorrect input line, because timing is not an argument of psi but a keyword of its own.

In some cases there is no '=' following the keyword and there are no commas separating the arguments. A typical example is:

    rd    sin     36   3.800    5.60

This input format is used in the PRIMITIVE-BASIS-SECTION and, similarly, in the INIT_WF-SECTION. In the operator-file (*.op) (see Hamiltonians) there may appear functions with arguments. Such a construct, as e.g.

CAP_rd = CAP [ 5.0  0.357    3 ] 

may appear in the OPERATOR-SECTION of the input-file when 'alter-labels/end-alter-labels' is used. Note that these latter inputs are not free-format, the order of the arguments matters! Moreover, the fixed-format inputs must have a line for their own; they must not be followed by another keyword.

Units.

The unit fs (femto seconds) is assumed, when times are specified. For all other input variables au (atomic units) is assumed as default. However, one may chose other units by separating a unit-keyword by a comma (i.e. treating it as an argument). Possible units and conversion factors are shown in the following table. There are energy-, mass-, length-, and momentum-conversion factors. Note: the input value is divided by the conversion factor to arrive at atomic units. However, for the keywords nmwl and eV*AMU the conversion is not performed by a simple multiplication with a factor. For a quick information on the units used in QUANTICS, run mhelp -s units.

SECTIONS

The section XXX starts with the keyword

and ends with

These keywords must be alone on a line.

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.

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

The RUN-SECTION contains keywords defining the calculation to be made, and general options for the run. There are different types of keywords, listed below.

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 necessary lower level keywords are automatically included. Thus

or

are equivalent.

The exact keyword can be added to any run-type keyword. For example,

will result in a numerically exact wavepacket propagation (i.e. the wavefunction will be represented in the full product primitive basis). Similarly,

will set up a wavepacket for a numerically exact calculation.

Notes on Improved Relaxation.

When the keyword relaxation is given, a normal relaxation, i.e. a propagation in negative imaginary time, is performed. However, relaxation=<...> will enforce an improved relaxation run, where the A-vector is not determined by relaxation but by diagonalisation of the Hamiltonian in the basis of the SPFs. Improved relaxation requires that one uses a CMF integrator scheme in the fix or varphi mode. (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 is recommended to use CMF without any extension.

For diagonalising the Hamiltonian in the basis of the SPFs one may use the SIL "integrator" (now, of coarse, a diagonaliser), or, which is usually the better choice, a Davidson routine. There are actually several Davidson routines implemented, called DAV, rDAV, rrDAV, and cDAV. See INTEGRATOR-SECTION for details. In short, DAV is for (general) hermitian Hamiltonians, rDAV and rrDAV are for real-symmetric Hamiltonians, and cDAV is for complex non-hermitian Hamiltonians, i.e. for computing resonances. Actually, the rrDAV keyword calls the rDAV routine but real arithmetic is then used for the H*A operation. rDAV alone stores the A-vectors as real. All DAV "integrators", except cDAV, can also be used in block form, i.e. for converging a set of eigenstates simultaneously. If the keyword relaxation=0 is given, the algorithm will simply converge to the ground state, or, in a block improved relaxation run, to a set of states of lowest energies. If relaxation=I is used, with 0 < I < 900, then the I-th state (counting from 0) of the very first diagonalisation is taken, and for the following diagonalisations the state selection is done with follow (see below). However, the I-th state of the Davidson (or Lanczos) matrix will in general not be the I-th state of the Hamiltonian, but a higher one. To minimize this effect there is the keyword follow which forces the routine to take as many Davidson/Lanczos iterations as allowed by input. In any case, relaxation=0 with I > 0 is only for testing, to compute excited states one should use relaxation=follow or relaxation=lock. With follow the Davidson diagonalisation takes that eigenvector which is closest to its start vector, i.e. to the A-vector of the initial WF in case of the very first diagonalisation, or the selected eigenvector 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. I.e. it accounts for changes in the SPFs. 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.

The rDAV routine allows to additionally use the keywords skip1dav, quad, olsen, backrotate, and cn, (where n denotes an integer, e.g. c6) as arguments to relaxation.

The keyword skip1dav lets the program skip the first Davidson diagonalisation, i.e. the calculation starts with orbital relaxation. This is only useful, if a relaxation is re-started and the WF is read in via the file statement. Then the A-vector was already obtained by diagonalisation and it does not make sense to diagonalise again. However, when build or autoblock is used, skip1dav should not be set!

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 useful particularly when a preconditioner is used. (Keyword precon). The Olsen correction in essence turns the Davidson algorithm into a Jacobi-Davidson one.

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 cn lets the program perform n orbital relaxation steps, before a new diagonalisation step is done.

When a calculation starts converging to a desired state but after a few iterations jumps to another state, then the numbers of single-particle functions are too small. It may also be necessary to increase the Davidson order. Inspect the rlx_info file (by running the script rdrlx) to see what Davidson 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 file rlx_info contains a lot of information on the relaxation process. Use the script rdrlx to read it. The first line of the output of the script rdrlx is labeled with a negative time. This line shows the expectation value of the Hamiltonian with respect to the provided initial state. The second line, labeled with time=0, shows the energy after the first diagonalisation, but without SPF relaxation. The following lines refer to SPF relaxation and subsequent diagonalisation.

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.

For an example input see the 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 .

A block-improved-relaxation is performed when the packets keyword is given in the SPF-Basis-Section. The DAV, RDAV or RRDAV keyword must be given in the Integrator-Section and the keywords lock,quad,backrotate,cn are not implemented for block-improved-relaxation. For an example input on block-improved-relaxation see blkHONO.inp.

Notes on Diagonalisation.

Notes on Thermal Relaxation.

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

In the OPERATOR-SECTION the operator to be used is specified. Parameters and labels defined in the operator file may also be altered here. The opname keyword is compulsory, all others are optional.

The default value for oppath is the path of the operator directory created during the installation of the QUANTICS package.

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.

DIRDYN-SECTION

In the DIRDYN-SECTION keywords required to run a direct dynamics calculation can be given. The keyword direct in the run-section must be given to activate this.

DDPOINTS-SECTION

    DDPOINTS-SECTION   
       .......  
    END-DDPOINTS-SECTION
    

    atnam   x  y   z
    

If non-cart is given the input is in maas-frequency scaled normal mode coordinates

    label   Q
    

In addition to the DDPOINTS-SECTON, if a transfile is not provided, the transformation matrices for the coordinate transformation can be given in a DDTRANS-SECTION. This uses the format in VCHAM, which can be used to generate the matrices. The structure is the usual

    DDPOINTS-SECTION   
       .......  
    END-DDPOINTS-SECTION
    

    reference
    atnam   x  y   z
    .....
    

    masses
    atnam   mass [, unit]
    .....
    

    forward
    f11  f21  f31 ....

    f12  f22  f32 ....

    .....
    

    backward
    b11  b12  b13 ....

    b21  b22  b23 ....

    .....
    

DIABREF-SECTION

If a different initial point for the DB generation is needed this can be provided in a DIABREF-SECTION. This is necessary if the initial geometry for the propagation (centre of the WP at t=0) is at a point with degeneracies between states where the QC cannot be trusted and the derivative couplings are undefined. It is best to choose a point stepped away from the degeneracy along a totally symmetric mode so as to provide a good definition of the diabatic states, reatining their symmetries. If using symmetry, this point must be far enough away from a high symmetry point that any symmetry replicas are at least mindb away or else an averaged structure will be produced, which could be back at the high symmetry intersection point.

The coordinates of the point can be given in either normal modes (5 to a line)

    nmodes
    q   q1   q2   q3   q4   q5
    q6  q7   ....
    

    cartesian
    ATNAM  X   Y   Z
    ....
    

DIRDYN-GB-SECTION

In the DIRDYN-GB-SECTION keywords required to run a direct dynamics calculation using a grid-based dynamics method (standard or MCTDH methods) can be given. The keyword dd-gb in the run-section must be given to activate this.

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

where the degrees of freedom labelled mode_label1, mode_label2 etc. are contracted together in a single mode, and the number of single-particle functions for this mode are nspf1 in the first set, nspf2 in the second set etc. More than one mode definition can be written on a line.

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

or

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

To contract the degrees of freedom X and Y into a single mode,

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

If a multi-set basis is used, then to have 3 functions in the first set and 2 in the second,

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

If all degrees of freedom are combined into a single mode then the keyword Nfuncs can be used. This is useful in a vMCG calculation with many degrees of freedom.

spf-basis-section
multi-set
    nfuncs = 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

Mode 2 is now identified with mode 1 and the calculation is faster, because the propagation of mode 2 is skipped. The id keyword is, of course, very important when identical particles, e.g. bosons, are treated.

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.

The following keywords, if given, define the type of GWPs to be used in a G-MCTDH or vMCG calculation. The default type is separable frozen, normalised and phaseless functions.

ML-BASIS-SECTION

For a ML-MCTDH calculation the SPF-BASIS-SECTION must be replaced by a ML-BASIS-SECTION which defines the ML-tree. The top layer (up-most A-vector) is indicated by 0> and the following ones by 1>, 2>, etc. The numbers following these symbols gives the numbers of SPFs. E.g.

  0> 5 5 5  

   2> [q1]
 

   2> [q1,q2]
 

 mlbasis-section
 0> 5 5 5
    1> 5 5
        2> [q1]
        2> [q2]
    1> 5 5
        2> [q3]
        2> [q4]
    1> 5 5
        2> [q5]
        2> [q6]
 end-mlbasis-section
   

 ML-Basis-Section
 0> 2 2
   # Electronic
   1> [el]
   # Vibrations
   1> 4 4
     # Main system
     2> 4 4
       3> [v10a v6a]
       3> [v1 v9a v8a]
     # Bath
     2> 3 3
       3> 2 2 2
         4> [v2 v6b v8b]
         4> [v4 v5 v3]
         4> [v16a v12 v13]
       3> 3 2 2
         4> 2 2
           5> [v19b]
           5> [v18b]
         4> 2 2
           5> [v18a v14]
           5> [v19a v17a]
         4> 2 2
           5> [v20b v11 v7b]
           5> [v16b]
 end-mlbasis-section
   

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

where n(s,t) is the number of SPDOs for the pair of states (s,t).

PRIMITIVE-BASIS-SECTION

The definition of the primitive basis used to describe the system being studied is written on one line per degree of freedom. The input is free format, with blanks dividing the various parameters. The format for each line is

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^m3ⁿ where m and n are positive integers. One can use the utility script find235.py to find numbers that have this required form.

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.

For a GWP basis basis_size is not required in the input.

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

(See Remarks on sphFBR)

Remarks on electronic basis:

A discrete electronic basis can be defined with the basis type el, with the basis size representing the number of electronic states. By default, the el basis is taken as a complete, time-independent basis and nothing further is required, i.e. SPFs do not need to be defined in the SBASIS-SECTION. It can be used in either single-set or multi-set representations. SPFs can be defined for an electronic basis to provide an MCTDH-like contraction. But this cannot be used in the multi-set formalism.

Remarks on continuum electronic basis:

To define continuum electronic states two consecutive lines are required in the section. The first uses the elcont basis type to define the number of electronic states in the problem, while the second line defines the DVR used to discretise the continuum. Thus the lines

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

define an electronic basis with 4 states, the first 2 of which are bound states while states 3 and 4 are coupled to a continuum. The continuum is then discretised using a sin DVR which sets 101 points between 0 and 3.75 eV. The single-particle functions for the continuum are then set up as for a normal mode using, e.g. a set of gaussian functions for the el mode in the SBASIS-SECTION, e.g.

   el        gauss   1.0, ev  0.0  0.3, ev

In a ML calculation, the elcont DOF needs to be explictely added to the tree. E.g. to add 8 spfs on the continumm grid of a system with 3 electronic states:

    0>  2 8 3
      1> [el] 
      1> [Elcont]
      1>  2 2 
         2> [x]
         2> [y]

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:

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
    

    s1   sin  2  -0.5  0.5  spin
    

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.

The use of the operators dq, qdq, and cdq is not allowed for cos-DVR.

Example: The following lines are equivalent.

  x   cos     36    2pi/2
  x   cos     36    0.00      3.1415926535897  long
  x   cos     36    0.0436332313  3.097959422  short
  x   cos     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.

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

in binary format:

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

The use if sphFBR is deprecated, PLeg or KLeg are preferred. WARNING: sphFBR does not work correctly if the potential is a potfit. Use PLeg instead. One may also consider to Fourier-transform the potential (with projection84) and then use KLeg. This approach is described e.g. in J.Chem.Phys. 123 (2005), 174311; J.Chem.Phys. 128 (2008), 064305; Mol.Phys. 110 (2012), 619-632. See also H2H2.inp and H2H2.op on the inputs or operatos directory, respectively

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 exp line end with the keywords k=-6,6. These are optional. Without this statement, K would span the full range from -7 to 7 (yielding 15 points). The restriction of the K-range is mainly for tests.

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.

Remarks on Wigner DVR:

Wigner DVR defines 3D single-particle functions, although the basis definition is for each degree of freedom individually. The TWO degrees of freedom following Wigner are part of the combined mode and must be either K or exp, or any combination; i.e. Wigner/K/K, Wigner/K/exp, Wigner/exp/K, and Wigner/exp/exp are all valid choices. As these combinations generate mode operators, all three DOFs must be combined in the SPF-Basis-Section.

Important Note: The order of the three Wigner-DOFs must be   beta, gamma, alpha   or   beta, k, m   when k denotes the BF- and m the SF-component of the angular momentum. This ordering must hold in the PRIMITIVE-BASIS-SECTION, in the combination scheme defined in the SPF-BASIS-SECTION, and in the modes line of the HAMILTONIAN-SECTION. Note that the angles (beta,gamma,alpha) correspondent to the angular momenta (j,k,m), respectively.

Example:

PRIMITIVE-BASIS-SECTION
  ........
  beta    wigner  20   all
  gamma   exp     15   2pi
  alpha   k       -7   7
end-primitive-basis-section

SPF-BASIS-SECTION
  ........
  beta,gamma,alpha = 20
end-spf-basis-section
   

Remarks on Gaussian Wavepacket basis:

A Gaussian Wavepacket basis does not require a grid as the SPFs are represented by parametrised Gaussian functions. A grid is, however, set up and used to help the analysis. For example the showd1d program uses the one-dimensional densities saved on grid points in the gridpop file. By default this auxiliary grid has 101 points and runs from -10.0 to 10.0. This is the case if no grid is specified. Specifying a grid tailors the calculation to the problem at hand, for example

    X   GWP
    Y   GWP    51  -2  2
  

defines a default grid for X and a shorter grid for Y.

INITIAL-GEOMETRY-SECTION

For a direct dynamics calculation (activated by the keyword direct in the run-section) the initial geometry does not need to be specified via a DVR. Instead it can be input in an INITIAL-GEOMETRY-SECTION. The first lines are for optional keywords, formatted one keyword to a line

After the keywords the geometry is specified in one or two blocks depending on the type of coordinates specified in the coords keyword. For Jacobi coordiates only a Cartesian block is required. For normal modes, both Nmode and Cartesian are required.

  Cartesian = S
  ATNAM     X     Y    Z     [ keywords ]
  End-Cartesian
  

  Nmode
  Label     Q  Freq [, unit] [keywords]
  End-nmode
  

DD-GB-SECTION

For a direct dynamics calculation using a grid-based dynamics method (activated by the keyword dd-gb in the run-section) the initial geometry is not specified via the DVR, but in the DD-GB-SECTION.

The geometry is specified in two blocks. For normal modes both Cartesian and nmode blocks are required.

	    Cartesian = S
	    ATNAM     X     Y    Z     [ keywords ]
	    End-Cartesian
	    
   

	    Nmode
	    Label     Freq [, unit]
	    End-nmode
  
   

INIT_WF-SECTION

In order to specify the initial wavefunction, one of the following options must be given.

      file=geninwfrun                  file=geninwfrun
      select s1 m1 1 3 4               select  1 3 4
      select s1 m2 1 4 6 3             select m2 1 4 6 3
      select s2 m2 1 3 5               select s2 m2 1 3 5 

WARNING FORTRAN cannot open a file twice. Hence it may become necessary to copy a file such that the same information can be read via two different file-paths. Such a problem may arise when using block-SPF and block-A, or when using orthogonalise.

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.

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.

After an initial WF (or density matrix) is generated, either by file, build, read-inwf, or block-SPF / block-A / autoblock this WF may be modified by one or several of the following procedures. The program will take actions in the order:

irrespectively of the order in the input file. Note that autoblock is executed in the propagation/relaxation step. Hence it is not meaningful to use operate or orthogonalise in connection with autoblock.

Full list of the possible options:

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

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,

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

where the integers s,t specify the initial states of the density operator. In a ket-bra notation this means |s><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

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.

Same as temperature above, but the frequencies defining the oscillators to be thermalised are taken from the operator parameters. This is required for calculations using mass-frequency scaled coordinates. The frequencies must be defined as parameters in the oper file, i.e. listed in the PARAMETER-SECTION of the .op file, or listed in the OPERATOR-SECTION of the .inp file in an ALTER-PARAMETERS block. The frequencies must be given the names freq0_I = freq, where I is the DOF number and freq the frequency.

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 protocolled in the log-file. If you want the eigenvectors too, specify "veigen" in the RUN-SECTION.
  The operator defined in a HAMILTONIAN-SECTION_XX must be a one-dimensional operator, i.e. only have terms that act on the approriate degree of freedom. For 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      |2  KE  
  omy      |2  q^2 
  ------------------------------
  end-hamiltonian-section
  

• 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, gauss, or Wigner 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. If the previous initial wavefunction is of type Wigner, the parameters here define the projections of angular momenta on the body-fixed and space-fixed z-axes for the first and second DOFs after Wigner, respectively, in the 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 PLeg and exp or Wigner and exp/K and exp/K.
  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  
  
  

• If Wigner is chosen, the single-particle functions are normalized Wigner (small)-d functions. The corresponding primitive basis must be Wigner. The initial wave function must be K for the following TWO degrees of freedom. The j parameter defines the first index (total angular momentum) of the Wigner function taken as the first single particle function. The indices of the angular momenta about the body-fixed and space-fixed z-axes are given in the following two K-inwf lines. If the parameter symmetry equals nosym then all (even and odd j) Wigner functions are used; if sym, then only even or odd Wigner functions are used, depending on whether j is odd or even.
  Note that the angles (beta,gamma,alpha) correspondent to the angular momenta (j,k,m), respectively.
  The following table shows which combinations of primitive-basis and initial-wf types involving Wigner are allowed:

  For Wigner initial-wf:
  ---------------------------------------------
  DOF  init_wf    allowed primitive-basis types
  ---  ------- || -----------------------------
   1   wigner  || wigner only
  2,3  K only  || K or exp
  ---------------------------------------------
  
  For Wigner primitive-basis:
  ----------------------------------------------------------------------------------------
  DOF  p-basis  || allowed initial-wf types
  ---  -------  || -----------------------------------------------------------------------
   1   wigner   || wigner | gauss or leg           | any 1D type (HO,gauss,leg,eigenf,map)
  2,3  K or exp ||   K    | K (for p-basis=K only) | HO,gauss,eigenf,map (p-basis=exp only)
  ----------------------------------------------------------------------------------------
      
  

  An example involving the use of Wigner-inwf:

  PRIMITIVE-BASIS-SECTION
  ---------------------------------------------------
  # Mode    DVR      N  center  freq     mass
  ---------------------------------------------------
    r       HO       8  1.626   0.01837  1739.96
    beta    wigner  16    all
    gamma   exp     15    2pi
    m       k     -7  7
  ---------------------------------------------------
  END-PRIMITIVE-BASIS-SECTION
  
  INIT_WF-SECTION
  build
  -----------------------------------------------------------
  # mode    type  q.n. center  mom.  freq.    mass
  -----------------------------------------------------------
    r       HO         1.426   0.0   0.01837  1739.96
    beta   wigner 10   nosym excite=mkj  print  # initial-j, use odd and even j
    gamma    k     1  -2  2  1                  # initial-k, kmin, kmax, dk
    m        k    -2  -5  5  1                  # initial-m, mmin, mmax, dm
  -----------------------------------------------------------
  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  | ....
  ----------------------------------------
             . 
             . 
  

  The map keyword can be used to sum up to three functions. A corresponding INIT_WF-SECTION may read:

  INIT_WF-SECTION
  build
   rd    map dum1  1.0
   rd    map dum2 -1.5
   rd    map dum3  (0.5,0.3)
   rv    HO     2.151  0.0  0.218,eV  13615.5
   theta gauss  2.22   0.0  0.0745
   dum1  gauss  4.310  0.0  0.075 pop=1
   dum2  gauss  4.315  0.0  0.080 pop=2
   dum3  gauss  4.320  0.0  0.085 pop=3
  

  The initial SPF for rd is a weighted sum of three harmonic oscillator functions. The first two coefficients are real, the third is complex. After summation, the function is normalized. The DOFs rd, dum1, dum2, and dum3 must be defined in the PRIMITIVE-BASIS-SECTION with identical grids.

Remarks on multi-set calculations:

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

The same function type must be used for each state, but the parameters can be different. It is however possible to choose between the different harmonic oscillator function sets HO, HO odd and HO even. Thus even functions can be used on one state and odd functions on the other. If the state keyword is not used, the same function parameters are used to generate the spfs for each state.

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

where the integer number p defines to which initial packet the corresponding input line belongs. It is possible to specify more initial packets in this section than given by the packets argument. All data with p > packets is then ignored.

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:

which populates the i th function.

Remarks on thermal relaxation calculations:

Which degrees of freedom are to be thermalized must be specified through the ran keyword in the INIT_WF-SECTION according to the order of parameters given in the table above labelled "parameter description". In the following we show four possible illustrative cases with an example. For the first mode the Ith state up to which the randomization will be averaged is not specified, thus the default value basis_size-1 will be employed. Meanwhile for the second mode, I is declared with a value of 30. Following this in the case of the third mode, the spf type eigenf is used together with the ran keyword. Finally, when mode-eigenfunctions are generated, the ran argument has to be included in the respective parameter's list line (meigenf keyword) as it is the case for the fourth mode. The meigenf procedure needs an initial SPF generated in the build-block. This initial SPF should not be an eigenfunction of the operator used in meigenf, but should have an overlap with several eigenfunctions. For this reason we have dispaced the origion of the HO-function by 0.25 au. Independently of the type of mode, the ran keyword will always be the last argument. Note that an electronic degree of freedom cannot be randomized. Electronic states are always considered as part of the system, not of a thermalized bath.

Example:

INIT_WF-SECTION
build
  init_state = 1
  q0001 HO     0.00  0.0  1.0 ran
  q0002 HO     0.00  0.0  1.0 ran=30
  q0003 eigenf oper    ran=30
  q0004 HO     0.25  0.0  1.0
end-build
meigenf = 4,oper,0,full,noselect,write,ran=30
end-init_wf-section
  

Initial Distribution of GWPs

1. Shell
Initially a set of 1D GWPs are built for each DOF with the first GWP at the q0,p0 specified by the input parameters. Subsequent GWPs are generated by stepping away from this point by a distance s specified either explicitely using step = R or to give the overlap specified by ovl = R. The steps are to negative and positive displacements in the order

 0  +s  -s  +2s  -2s .... 

  
DOF     1      2        
GWP
 1      0      0    
 2     +s      0    
 3      0     +s    
 4     -s      0    
 5      0     -s    
 6      s      s    
 7     -s      s    
 8     +s     -s    
 9     -s     -s    

  
DOF     1      2        
GWP
10    +2s      0
11      0    +2s
12    -2s      0
13      0    -2s
14    +2s     +s
15     +s    +2s
16    -2s     +s
17    +2s     -s
18     -s    +2s
19     +s    -2s
20    -2s     -s
21    +2s    +2s
22     -s    -2s
23    -2s    +2s
24    +2s    -2s
25    -2s    -2s

  
DOF     1      2      3  
GWP
 1      0      0      0
 2     +s      0      0
 3      0     +s      0
 4      0      0     +s
 5     -s      0      0
 6      0     -s      0
 7      0      0     -s
 8      s      s      0
 9      s      0      s
10      0      s      s
11     -s      s      0
12     -s      0      s
13      s     -s      0
14      s      0     -s
15      0     -s      s
       .....

  
 1     +s      0     0
 2     -s      0     0
       .....

 0  +s  -s  +2s  -2s .... 

  
DOF     1      2     3    
GWP
 1      0      0     0
 2     +s      0     0
 3     -s      0     0
 4      0     +s     0
 5      0     -s     0
 6      0      0    +s
 7      0      0    -s
 8    +2s      0     0
 9    -2s      0     0
10      0    +2s     0
11      0    -2s     0
      .......

  
 1     +s      0     0
 2     -s      0     0
       .....

if wigner_sym is used, the points are symmetrically distributed. I.e. if the random point q,p is chosen (where q, p are displacements from q0, p0), then 4 GWPs are created with (+q,+p), (-q,+p), (+q,-p), (-q,-p) . Hence the number of GWPs has to be 4n+1 .

4. Random
As Wigner, but the distribution is only made in coordinate space.

5. Read
The parameters for the full multi-dimensional GWPs are explicitely defined in a special section as follows

  
GWP-DEFINE-SECTION
define Ngwp [Nmode]  [Nstate]
sigma1  q1   p1
sigma2  q2   p2
....

END-GWP-DEFINE-SECTION

sigmaX qX pX are the standard deviation defining the GWP width, the centre coordinate and centre momentum respectively. One line is required per DOF and must be given in the order of the DOFs in the mode.

This option should be used with care as few checks are made.

6. Pqgrid, nq, np
The initial 1D GWPs are set up on a PQ grid, i.e. displacements are made in both P and Q, with dimensions nq, np as given by the input parameters. The grid is formed by distributing the GWPs in coordinate space as for the Shell option, and then adding displacements in momentum space. This sets up, e.g. for a 5,5 grid (nq = 5, np = 5) with displacements q and p, 25 GWPs are generated

  
GWP     Q      P
 1      0      0    
 2     +q      0    
 3     -q      0    
 4    +2q      0    
 5    -2q      0    
 6      0     +p    
 7     +q     +p    
 8     -q     +p    
 9    +2q     +p    
10    -2q     +p    
11      0     -p    
12     +q     -p    
13     -q     -p    
14    +2q     -p    
15    -2q     -p    
16      0    +2p    
17     +q    +2p    
18     -q    +2p    
19    +2q    +2p    
20    -2q    +2p    
21      0    -2p    
22     +q    -2p    
23     -q    -2p    
24    +2q    -2p    
25    -2q    -2p    

7. Simex
As for Shell, initially a set of 1D GWPs are built for each DOF with the first GWP at the q0,p0 specified by the input parameters. Subsequent GWPs are generated by stepping away from this point by a distance s specified either explicitely using step = R or to give the overlap specified by ovl = R. The steps are to negative and positive displacements in the order

 0  +s  -s  +2s  -2s .... 

DOF     1      2     3
GWP
 1      0      0     0
 2     +s     +s    +s
 3     -s     -s    -s
 4    +2s    +2s   +2s
 5    -2s    -2s   -2s
      .......

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

The following two inputs are equivalent. In both cases the SPFs of the first state of the gs-calculation (which may have only one state) is put on state 1, 2 and 3 of the initial WF. The A-vector of the first state of the gs-calculation is put on the second state of the initial WF. The A-vector for states 1 and 3 is zero.

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

One may use more than one (up to eight) restart files to build the initial wavefunction.

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

It is possible to multiply the A-vector of a specific electronic state by a complex factor. For this the keyword weight is used. If weight is not given, it is set to 1 by default.

Read-Inwf
  file = <restart-directory>
  SPF 1 -> 1
  SPF 2 -> 2
  A   1 -> 1   weight = 0.65
  A   2 -> 2   weight = (2.3, -0.96)
end-read-inwf

In the first A-line, weight is assumed to be real, i.e. the input is equivalent to weight=(0.65,0.0).

In a multi-packet run the different packets are treated as different electronic states. The Read-Inwf statement does not know packets and one has to do the assignment by "hand" using the formula s1 = (p-1)*nstate + s, where p and s denote the actual packet and state and s1 is the effective state number, which one must enter in the Read-Inwf block.

Read-Inwf works for MCTDH wavefunctions and density operators, but not selected CI wavefunctions. It also works for exact wavefunctions. In this case the syntax is slightly different, as there is no A-vector.

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

Generating block initial wavepackets

For a multi-packets single-set run, e.g. for block improved relaxation, there are several ways to generate an initial wavepacket. The obvious one is to use build and A-coeff. For example, for a CO₂ run with 4 packets (i.e. packets = 4,single-set in SPF-BASIS-SECTION) one may use:

INIT_WF-SECTION
 BUILD
  r1     gauss  1.20,Angst  0.0  0.030,Angst
  theta  gauss  3.14159     0.0  0.0524
  r2     gauss  1.25,Angst  0.0  0.025,Angst
 end-build

 A-coeff
  1 1 1 1 (1.d0,0.d0)
  2 1 1 2 (1.d0,0.d0)
  1 2 1 3 (1.d0,0.d0)
  1 1 2 4 (1.d0,0.d0)
 end-A-coeff
end-init_wf-section

If the block size grows, the definition of the A-vector with A-coeff becomes rather cumbersome. Here autoblock is very useful, as it generates zero order A-vectors for the lowest states. E.g.:

INIT_WF-SECTION
 BUILD
  r1     gauss  1.20,Angst  0.0  0.030,Angst
  theta  gauss  3.14159     0.0  0.0524
  r2     gauss  1.25,Angst  0.0  0.025,Angst
 end-build
 autoblock
end-init_wf-section

Note that the autoblock command is executed in the propwf step. Hence a geninwf run alone will yield a zero A-vector.

The SPFs may be read from a "foreign" restart file rather than be build. This foreign restart file may contain a standard single-packet WF (generated by relaxation or improved relaxation) or a multi-packets single-set WF (generated by block improved relaxation).

INIT_WF-SECTION
 block-SPF = GS
 autoblock
end-init_wf-section

Here GS denotes a directory which contains a restart file (e.g. name-directory of a relaxation to the ground state). Rather then giving a directory one may give the path of the restart file. Hence the above statement is equivalent to block-SPF = GS/restart. This feature is helpful in case the restart file is not simply called restart. In case the previous calculation had used an FFT it will be useful to take the real part of the SPFs only. This can be accomplished by setting the option realpsi, i.e. block-SPF=GS,realpsi. By default the SPFs are re-orthonormalized. To switch this feature off, set the option noorthopsi.

Initial A-vectors can be read with the aid of the block-A keyword.

INIT_WF-SECTION
 block-SPF = GS/rst000,noorthopsi
 block-A   = GS/rst001,GS/rst003,GS/rst005
 block-A   = GS/rst007,GS/rst009,GS/rst011
end-init_wf-section

There may be more than one block-A statement. The number of A-vectors read must match the block size. It is important to set the option noorthopsi here, because otherwise the SPFs will be transformed and will no longer be consistent with the A-vectors. One usually wishes to set relaxation=follow when A-vectos were read in (relaxation=lock is not implemented for block improved relaxation).

Important: Some compilers, in particular gfortran, do not allow to attach two read/write channels to the same file. Hence the file-path appearing in the block-SPF statement must not be identical to any of the file-paths appearing in the block-A statement(s). In rare cases it may become necessary to copy a file to have the same information on two different file-paths.

INTEGRATOR-SECTION

The string S=all after the integrator name may be omitted, i.e. ABM is equivalent to ABM/ALL.

For CMF calculations with multi-layer MCTDH, there are two different schemes of integrating the equations of motion. The scheme can be selected by the keyword ml-cmf. For unified propagations, the complete EoM is integrated as one, hence the choice of integrator is restricted to ABM/all, BS/all, or RKx/all. (For legacy reasons, here /spf has the same effect as /all.) For split propagations, the EoM for each mode is integrated by itself, hence for each mode one can choose its individual integrator/parameters. (Note: Currently the choice of integration algorithm is restricted to RK5/8 or (C)SIL; (C)SIL only makes sense for mode 1.) This is done by using S=modespec, where the mode specification modespec is one of the following:

The modes are numbered consecutively according to the ML-BASIS-SECTION. One can double-check the mode numbers by inspecting the visualization of the ML tree generated by the graphviz keyword [RUN-SECTION]. An example input using SIL for the top layer and a more accurate RK8 for node 3 is:

INTEGRATOR-SECTION
  ml-cmf  = split
  CMF/var = 0.03 , 1.d-6
  SIL/@1  = 20   , 1.d-6
  RK8/def = 1.d-6, 0.03, imp-ortho
  RK8/@3  = 1.d-8, 0.03, imp-ortho # for R
end-integrator-section
    

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.

The Davidson routines DAV, rDAV, and rrDAV can be used to perform block-relaxations. To this end set the packets keyword in the SPF-Basis-Section. Note that cDAV cannot be used in block form. All Davidson routines allow for parallelization (mpi).

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.

The equations for propagating vMCG GWPs have numerical problems due to (1) inversion of the C-matrix that may be singular (when GWPs are unpopulated) (2) inversion of the overlap matrix that may also be singular (when GWPs are linearly dependent). There are a number of options that control these problems. In general the default parameters work OK, but if problems during propagation are seen, it may be necessary to try other ways.

The inverse of the overlap matrix by default uses SVD to remove singular values below the threshold of epsgwp. This effectively projects out any linearly dependent functions. It is noted in the log file if the SVD removes any functions. If these events are causing obvious spikes in the dynamics, or else the integrator crashes with tine step sizes the eps_gsmat parameter can be changed, or the alternative of using regularisation for the inversion can be tried.

For stability of propagation it is useful to keep GWPs stationary that are not required to describe the evolving wavepacket. To understand the problem consider the limit when ony 1 GWP is populated. The trajectories of the unpopulated GWPs are undefined and random, introducing noise into the propagation. Also in the limit of a complete basis set, The C-matrix will be singular (0 everywhere) as it contains the projector out of the space covered by the functions. Again the trajetcories of the GWPs are not defined. In these limits, the undefined functions should be time-independent. For this reason, the CX decomposition that enables the GWP trajectory to be written as a classical trajectory + coupling is not used by default.

The inversion of the singular C-matrix is controlled in 2 ways. The first is how the coupling between GWPs is treated, controlled by the gwp_coupling keyword. The default is to use dynamic coupling, where only GWPs with significant values of C-matrix elements are included in the propagation. The others, by default, are stationary unless the CX option is selected in which case they follow un-coupled classical trajectories.

In dynamical coupling, whether a GWP parameter should be time-depedent is determined by the relevant diagonal C-matrix element. If this is above a certain parameter, controlled by eps_dyn, they are allowed to move. At the end of the log file there is some information that helps see this behaviour. For example the following output is from a 2-mode 2-state calculation.

**** GWP basis dynamic allocation ****
 mode  state      Ndyn /   N      Tmax      Cmin
    1    1          4 /   10      0.500  0.110122E-29
    2    1          6 /    8     28.500  0.303003E-08
    1    2         10 /   10      8.250  0.239499E-06
    2    2          5 /    8     26.000  0.128838E-08

N is the number of GWP parameters (this is a calculation using 2D frozen GWPs so there are 2 parameters per function ignoring the phase as this is defined by the time-evolution of the other paramaters). Ndyn is the maximum number of these parameters that become dynamical. Tmax is the time at which this state is reached, i.e. when the latest parameter is added to the dynamical set. Cmin is the smallest value of the diagonal C-matrix elements. This helps to show the quality of the basis as if not all the functions are required it should be reasonable.

In addition to the GWP coupling, the inversion is controlled by using a reularisation scheme. The default is Tikonov regularisation, which has been found to be the most rebust, but other schemes can be chosen, such as regularisation of the eigenvalues, or SVD. It is also possible to use full, direct inversion of both the C-matrix and the overlap matrix if greater accuracy is required. These may fail though if the EOMs are too badly conditioned.

Finally, the integrator also plays a role. The default 5th order Runge-Kutta parameters are selected to balance speed with accuracy. If greater accuracy is required, the tolerance should be decreased, or else 8th order RK tried. If there is severe numerical instability, e.g. only a narrow range of inversion and integration parameters work, this probably points to a problem with the Hamiltonian. Note also that if an LHA is being used, the energye will not be conserved no matter how accurate the integration.

If the numerical stability is an issue, Basis Expansion Leaping may be used for propagation. This keeps all GWP parameters fixed and the propagation is performed only in the expansion coefficients. The basis will be adapted repeatedly after a short time interval by re-expanding the wave packet into a new set. Keeping the parameters constant eliminates stability issues and drastically lowers the computational cost such that large sets of houndreds of GWPs can be used. See (WK,TJF PRL 2012) for details.

Tests have so far only be performed for exclusively GWP basis sets in joined dimensions with one state.

Three keywords in the integrator section control BEL:

Omitted integrator parameters default to:

DDPOINTS-SECTION

  symmetry_reference [ = unit ]
  ATNAM     X     Y    Z     
  end-symmetry_reference

        

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 explicitly included)
• Further example inputs

Example Inputs for use of CPP with QUANTICS

Input containing #include directives, invoke QUANTICS with -cpp:

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

Include files for use with both of the above input files:

• 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)