QUANTICS: Input Documentation

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

END-INPUT

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

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

keyword # comment

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

#include "some-common-settings.inc"

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

/* some comment */

The '/*' is interpreted by QUANTICS like '#', i.e. everything which follows til 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 ivalue 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.

Keyword Description Conversion factor
au atomic units 1 (may be omitted)
mH milli Hartree 1000
ev electron volts 27.21138386
mev milli electron volts 27211.38386
cm-1 wave number cm^-1 2.1947463137d+5
kcal/mol kcal/mol 627.503
kJ/mol kJoule/mol 2.6255d+3
Kelvin Kelvin 3.15777d+5
nmwl wavelength in nanometer 107/cm-1
aJ atto Joule (10-18J) 1.602177d-1*ev
invev (electron volts)^-1 0.036749325398
debye unit of electrical dipole moment 1 / 0.39343
AMU atomic mass unit 1 / 1822.88848325
p-mass mass of proton 1 / 1836.15267247
H-mass mass of Hydrogen atom 1 / 1837.15264409
D-mass mass of Deuterium atom 1 / 3671.482934845
Angst Angstroem 0.52917720859
pm picometer 52.917720859
nm nanometer 0.052917720859
deg degree 180.0/pi
Angst-1 Angstroem^-1 1/Angst
pm-1 picometer^-1 1/pm
nm-1 nanometer^-1 1/nm
eV*AMU sqrt(2mE)
fs femtosecond 1/41.34137333656 = 0.02418884326505
ps picosecond 1/41341.37333656 = 2.418884326505d-5
invfs femtosecond^-1 41.34137333656

SECTIONS

The section XXX starts with the keyword

XXX-SECTION

and ends with

END-XXX-SECTION

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.

XXX Description
RUN Whether propagation, relaxation, or diagonalisation, whether to generate or read DVR information, propagation and output-times, etc.
PRIMITIVE-BASIS Definition of primitive basis, e.g. DVR / FFT. Note: PBASIS is a short form for PRIMITIVE-BASIS .
SPF-BASIS Definition of the single-particle function basis, e.g. whether to combine any degrees of freedom, how to treat an electronic degree of freedom, how many functions etc. Note: SBASIS is a short form for SPF-BASIS .
ML-BASIS Definition of the single-particle function basis and structure of the wavefunction for a multi-layer calculation. How many layers, how to combine the degrees of freedom, how to treat an electronic degree of freedom, how many functions etc.
SPDO-BASIS This is a special feature of density operators of type I. The section corresponds to the SPF-BASIS section. It accounts for the special structure of this type of density operators.
OPERATOR Which operator to be used, any parameter changes to be made etc.
INIT_WF How to generate the initial wavefunction.
INTEGRATOR Which integrator to be used, and with what parameters.

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

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.

Keywords defining how output files will be treated
name = S The name of the calculation. A directory with name S is required in which files related to the calculation, unless otherwise explicitely stated, will be written.
overwrite Any files already in name-directory will be overwritten.
If this keyword is not given, and files already exist in name-directory, calculation will stop without these files being overwritten. It is recommended not to set overwrite but rather to use the option -w .
title If the keyword title appears in one line of the run-section, then the next line is supposed to contain a headline title of the run. This line will be written to output, timing and log files. Note the title line will be read regardless of special characters like '#' or '!' . Alternatively, when an equal sign, = , follows the keyword title, then everything that follows, from the equal sign till the end of the line, is regarded as title.
The following keywords define the method used to describe the system
Keyword Description
wavefunction The system will be described by a wavefunction (default)
density1 The system will be described by a density operator, using the type I formalism
density2 The system will be described by a density operator, using the type II formalism
The following keywords define the calculation to be made
Keyword Level Description
gendvr 1 A DVR file will be generated.
genoper 2 An operator file will be generated.
genpes 2 An operator file, called pes, will be generated with only the pes part of the operator.
gengmat = I1,I2 2 A special operator file, called pes, will be generated, to be used by showsys to plot the (I1,I2) element of the G-matrix of the kinetic energy. (See below).
geninwf 3 An initial wavefunction will be generated.
test 4 All input files will be checked and all other files, necessary for a propagation, will be created, but no propagation step will be performed.
propagation 4 Propagation in real time.
relaxation 4 The propagation will be in imaginary time i.e. the wavepacket will be relaxed to the ground state.
relaxation = I(,S1(,S2(,S3(,S4))))
relaxation = S(,S1(,S2(,S3(,S4))))
4 Improved relaxation. Requires CMF integration scheme.
If an integer I is specified, the I-th eigenstate will be computed (I < 900).
However, only the use I=0 is recommended, see notes below.
If I is replaced by the string follow then the eigenstate closest to the previous WF is computed.
The strings full and ortho may be given as second or third arguments. (ortho only for SIL).
If the Davidson integrator, DAV, is used, the inputs relaxation= I, relaxation=follow or relaxation=lock may be used. The strings full skip1dav, quad , olsen , backrotate , quadphi , and cn (where n is an integer), may be given as additional arguments.
See remarks below.
continuation 4 A continuation of the run in the name-directory will be performed.
diagonalisation = I 4 The Hamiltonian will be diagonalised using a Lanczos algorithm with I iterations. The SPF-BASIS-SECTION and the INTEGRATOR-SECTION are ignored for a (Lanczos) diagonalisation run. The WF is expanded in exact format.

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

propagation

or

gendvr genoper geninwf propagation

are equivalent.

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.

Keywords augmenting the calculation type
exact A numerically exact wavepacket calculation will be made.
g-mctdh If GWP basis functions are present, the G-MCTDH algorithm will be used with the GWPs used directly as SPFs. If a GWP basis is present and this keyword is not present, an MCTDH/G wavepacket calculation will be made in which the SPFs are expanded in GWPs.

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.

The following keywords define calculations that generate files to help the analysis of calculations.
Keyword Level Description
genoper=S 2 An operator file with the name S will be generated from the .op. This can be used together with the EXPECT program to calculate the time dependence of the expectation value of an operator.
genpes 2 A pes file will be generated from the .op file. This file can be used together with the SHOWSYS program to plot the potential energy surface, or together with VMINMAX program to determine minima and maxima of the potential energy surface.
gengmat = I1,I2 2 A pes file will be generated from the .op file, containing the (I1,I2) element of the G-matrix defined by the kinetic energy part of the Hamiltonian. This file can be used together with the SHOWSYS program to plot a energy surface. This keyword, gengmat, is introduced for testing the kinetic energy operator for correctness.

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

Keywords defining how the read-write files will be handled.
readdvr = S The DVR information will be taken from the DVR file in directory S.
readoper = S The operator information will be taken from the OPER file in directory S.
If the operator is read, it is a good habit to read the DVR file as well.
readinwf = S The initial wavefunction will be taken from the RESTART file in directory S.
deldvr The DVR file will be deleted at the end of the calculation.
deloper The OPER file will be deleted at the end of the calculation.
Keywords defining a direct dynamics calculation.
Options are specified in the DIRDYN-SECTION. For more details see the How To.
Keyword Description
direct = S A direct dynamics calculation will be run. The string S defines the type of coordinates for the dynamics. This defines the KEO. See notes below.
S = Jacobi Cartesian Jacobi's used
S = Nmodes Normal modes used
ngwp = I[, S] The no. of GWPs to be used in a direct dynamics calculation. If S = multi then I GWPs are used for every electronic states. If S not used, an Ehrenfest-type calculation is made.
partition = I The GWP wavepacket in a direct dynamics is partitioned into I modes with the DOFs distributed evenly among them. The no. of GWPs / mode is taken to be even, with the product N^p >= ngwp.
use_dd_sym The direct dynamics run will use and retain the symmetry of the initial wavefunction if run in normal modes and with correct number of GWPs. The correct number is that which fills the GWP shell (the minimal full shell is 2*N+1, where N is the number of modes).
Keywords associated with a propagation or relaxation calculation
Keyword Description
tfinal = R The propagation will run up to a time of R fs. Length of propagation is tfinal - tinit.
tinit = R The propagation will start at time R fs.
tout = R
tout = S
The output will be written every R fs. If this keyword is omitted, output will only be written at the end of the calculation. With S=all, i.e. tout=all, output will be written after each CMF-step. This is useful for improved relaxation. If tpsi is set in addition to tout=all, then output will be produced at multiples of tpsi also.
tpsi = R The wavefunction vector will be written every R fs. If this keyword is omitted, the vector will be written at the same time as the output. Note: tpsi must be an integer multiple of tout.
tstop =S The stop-time (real-time) is written to the stop file. The format is hhh:mm (i.e. 009:25). The job will be stopped after the first output after hhh:mm real-time. Note: if the stop time is next day, add 24 to the hours. E.g. 057:25 will stop the run two days later after 9:25 (24+24+9=57). NB: If the stopfile is removed, the run will stop after the next output.
tcpu = S The stop-time (cpu-time) is written to the stop file. The format is hh:mm:ss (i.e. 00:49:30). Alternative formats are Is and Im, where I denotes an integer. I.e. the inputs 120s, 2m, and 00:02:00 are equivalent. The job will be stopped after the first output after hh:mm:ss cpu-time. NB: If stop-file is removed, the run will stop after the next output.
twall = S Similar to tcpu, but the elapsed wall-time (i.e. real time) is compared with twall. The allowed input formats are similar to tcpu. Note that tstop, tcpu, and twall can alternatively provided as options (see mctdh84 -h).
energy-not-eV The energy in the output file will be given in au. This is for running models in dimensionless energy units.
time-not-fs The input times are in au, and output times are given in au. This is for running models in dimensionless coordinates.
normstop=R The program will be stopped if norm < R. This is the default with R = 10^-16.
no_normstop The program does not stop reagardless of the norm (normstop criteria ignored)
natpopstop=R(,I(,I1)) The program will be stopped when the lowest natural population for the mode number I exceeds the threshold R. If I is omitted or the string "all", stop when this criterium is reached for all modes. If I is the string "any", stop when this criterium is reached for any one mode. If I1 is given, the check applies only to state I1, otherwise to all states.
converged= R(,S) An improved relaxation run with the Davidson integrator will be stopped if the sum of the two last absolute energy changes is < R. The string S may specify an energy unit. A useful choice is converged=1.0d-5,eV
precon=I An improved relaxation run with the real Davidson integrator (RDAV or RRDAV) may use a better pre-conditioner than just the diagonal. If I.gt.1 a IxI dimensional block of the hamiltonian matrix is build, inverted and used as pre-conditioner.
split-rst If a block-Davidson improved relaxation run is performed, the use of this keyword will finally split the multi-packet restart file into a series of single-packet ones. The latter are called rst000, rst001, etc. Note, if the restart file rstxxx is read by a program which in addition reads the dvr and oper files (e. g. showsys) the latter files have to be created with a consistent primitive basis set. To create these dvr and oper files, run mctdh84 using the previous input file, but with a new name, with genoper as task, and with the packets keyword removed. Then move the files rstxxx to the new name-directory.
split-rst works also for geninwf runs. I.e. replace relaxation with geninwf and use the file statement to read a restart file of a block-improved-relaxation. (Here it is recommended to set noorthopsi when the numbers of SPFs are unchanged).
rlxunit=S(,R) The final energy after each iteration step of an improved relaxation calculation using the DAV "integrator" is output in the energy-unit S to the rlx_info file. Additionally an energy-shift may be applied, i.e. the number R is subtracted form the eigenvalue. If rlxunit is not given, S=eV and R=0 is assumed.
rlxemin=R(,S)
rlxemax=R(,S)
These two keywords define an energy window (R=energy-value, S=energy-unit), in which a relaxation=lock run searches for the eigenstate of the maximal overlap with the inital WF. This allows for convergence towards an eigenstate which is not the one with the (total) maximal overlap. Note: The energies given are with respect to a possible energy shift defined by rlxunit. The input is ignored, if the run-type is not relaxation=lock for a normal improved relaxation or relaxation=0 for a block relaxation. In the latter case convergence is towards eigenstates with lowest energies above rlxemin. rlxemax does not make sense and should not be used for relaxation=0 runs.
reflex If the system operator is build on the reflex algorithm (see JCP 134 (2011), 234307) while using the two "electronic" states simultaneously (identical SPFs for both states), then the reflex keyword splits the A-vector accordingly. This saves both, CPU-time and memory. If the plus and minus states are calculated independently, the reflex keyword must not be given.
freeze=I1(,I2,..) The numbers I1, I2,.. define the modes to be frozen, i.e. these modes are not propagated or relaxed. Useful for checking in an improved relaxation run which modes couple strongly. Probably not useful for propagation runs.
realphi This keyword should only be given in a improved relaxation run, real variant (i.e. RDAV or RRDAV). After each orbital relaxation, the real part of the SPFs will be taken. This is useful, if FFT or PLeg are used, because these representations may contaminate the SPFs with small imaginary parts. To compensate for a possible norm loss of the SPFs when taking real parts only, the use of imp-ortho in the Integrator-Section is recommended. The sizes of the removed imaginary parts is reported in the log-file. If FFT or PLeg are not used, there is no point in setting realphi.
optcntrl(=S) This keyword and its arguments are usually set by the script optcntrl.
optcntrl is set if an optimal control problem is to be solved. optcntrl requires propagation. The optcntrl keyword can have several arguments. If optcntrl=pc is given, then a predictor/corrector algorithm for determining the electric field is assumed. If optcntrl=simprop is given, two instances of MCTDH simultaneously propagate the initial and target states. The instance using the on-the-fly calculated new field requires the the additional keayword update, i.e., optcntrl=simprop, update. Note: simprop and pc cannot be used together.
local-control The run is an local control theory (LCT) run. local-control requires propagation. The type of local control and strength parameter for the pulse is defined in the operator. If a seed pulse is required tinit must be negative. This keyword will open the file field.lct
natfast = V1,V2{n},... An more efficient algorithm for H(natpot)*A operation is used. The A-vector is pre-multiplied with several natpot terms and the results are stored for further use. (Hence fast requires slightly more memory). The labels V1, V2, ... are the labels assigned to a natpot in a Labels-Section of the operator file. The (optional) specification of an order n is possible by attaching this number in curly brackets ({}) to the label. n denotes the number of modes used for pre-multiplication. The larger n is, the larger will be the speed-up. However, n is limited to min(4,nmode-2), where nmode denotes the number of (combined) modes (or MCTDH particles). The n orders are optional, i.e. if no order for the current label is given, n=min(4,nmode-2) will be used. If no potential label is given, i.e. only fast (or fast = all), the maximum order for all natpots will be set. Some information about using the "Fast" is written to op.log file. "Fast" works also for a multipackage but not for a multi-set propagation and not for S-MCTDH.
seed=I1(,I2,..) The numbers I1, I2,.. define seed values for the random number generator. Current implementation uses fortran library random generator, thus number of seeds is compiler dependent (ex. gfortran needs 12 values, ifort needs 2). If wrong number of seeds is used the program will issue a warning.
usempi (=S(,S1(,S2,..))) The program runs in a distributed memory parallelised modus using MPI if started with the mpirun command. The parallelisation of the different parallel subroutines may be switched off by setting the keywords no-summf, no-calcha, no-funka2, no-phihphi, no-hlochphi, no-mfields, no-getdavmat, no-dsyev or no-dav, no-getdiag. With the dav = I keyword the number of Davidson vectors per MPI process can be adjusted in order to reduce communication cost.
openmp=I Define the number of OpenMP threads to use. Parts of the code can specifically be turned on by adding the following keywords:
openmp_qc
Perform quantum chemistry calculations for direct dynamics on different threads.
openmp_mfield
Parallelise mean-field calculation (buggy implementation).

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

Keyword Default
tinit 0.0
tout tfinal-tinit
tpsi tout

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

Keywords defining output-data files to be opened.
Keyword Description
all All the standard optional files discussed below will be opened. More precisely, a propagation/relatation run will open the files: auto, autoe, gridpop, output, pdensity, psi, speed, steps, stop, timing, and update. In a diagonalisation run, only the files output, timing, lanczvec, and eigvec are opened. vMCG-specific files are not included in the set of files opened by all.
activegwp The identities of the active degrees of freedom in a dynamic GWP run are written to activegwp. In this file 1 means active, 0 inactive.
auto = S (,S1) (,S2) The auto-correlation function will be written to the file auto.
If S = twice, auto-correlation function is written twice in interval tout (only for CMF) .
If S = once, auto-correlation function is written only once in interval tout.
If S = no, no auto-file will be opened. This only makes sense, if the keyword all was given previously.
If S, S1 or S2 = order1, the file auto1 will be opened.
If S, S1, or S2 = order2, the files auto1 and auto2 will be opened.
These files contain the first and second order auto-correlation function, respectively, needed in the filter-diagonalisation method.
Note that in a multi-packet run, i.e. when packets > 1, the files auto, auto1, and auto2 contain cross- rather than auto-correlation functions.
auto Synonym for auto = once.
cross (= S (,I)) A cross-correlation function will be calculated and written to the file "cross". The reference wavefunction will be taken from the "restart" file residing in the directory named S.
If S=name, the restart file will be taken from the current name directory, which means that the auto-correlation function will be calculated (but without the t/2-trick). This is also the default (i.e. if no argument is given). If an integer I is given after the path S, then the cross-correlation will be evaluated for the electronic state I only. WARNING: If the restart file is taken, i.e. if no path is given, one cannot continue a calculation, because the continuation run will start with a different restart file leading to false cross-data.
ctrace = D The trace tr(D rho) is written to the file ctrace in a density operator propagation. If no operator is given, D = |s><t|, where s=left_state and t=right_state (see Init_wf-Section).
eigvec In a diagonalisation run, the eigenvectors of the tridiagonal Lanczos matrix are written to the eigvec file.
expect = S (,S1, S2, ...) The expectation value of the operator S, <psi(t)|S|psi(t)> / <psi(t)|psi(t)>, is evaluated and written to the file expectation. Up to maxham operators may be specified. (If S=system then the expectation value of the whole System-Hamiltonian is derived, i.e. the total system-energy.) The norm (not norm**2) of Psi is additionally written to the expectation file.
For density operators the expectation value is tr(H rho)/tr(rho).
There may be more than one expect line. I.e
expect = S, S1, S2 is equivalent to
expect = S
expect = S1
expect = S2

When the first argument to expect, S, is real-only , then only the real parts of the expectation values will be output to file expectation.
expect1 = S (,S1, S2, ...) Same as expect, except that the data is written to the file expect1. This second expectation file is useful for a better organisation of the data if several expectation values are computed. One file may store real, the other complex expectation values.
Important Note: The expect1 keyword(s) must appear in the input file after the expect keyword(s). There must not be an expect1 keyword without a previous expect keyword.
gridpop The grid populations will be written to the file gridpop. Note: The grid populations of the different states will be summed.
gridpop=el The grid populations will be written to the file gridpop. The grid populations of the different states of a multi-set run will be stored separately.
gwpcentres/gwpcenters Similar to sheptraj, the GWP centre coordinates are written to the gwpcentres file.
lanczvec In a diagonalisation run, the Lanczos vectors are written to the lanczvec file.
orben The orbital energies, i.e. the eigenvalues of the trace of the mean-field operators, are calculated and written to the orben file. Note: orben must be set, when the propagation is in energy orbitals (keyword energyorb, Integrator-Section)
output The output will be written to the file output rather than to the screen. (default).
screen The output will be written to the screen rather than to the file output. Alternatively to screen one may give the keyword no-output.
pdensity (=I1,I2,I3,I4) The one-particle density will be written to the file pdensity. If the pdensity keyword is followed by an equal sign and up to four integers, the one-particle density will be output only for the specified (contracted) modes.
psi = S (, S1 or R) The wavefunction will be written to a file every tpsi fs. If no arguments are given, it is written single precision to the file psi.
If S or S1 = single, the wavefunction will be written single precision.
If S or S1 = double, the wavefunction will be written double precision.
If S or S1 = natur, the wavefunction will be written as natural orbitals. This option is automatically taken if natural orbitals are propagated.
If S = compact, the wavefunction is written in natural orbital representation and compact form using the cutoff R.
If S or S1 = no, no psi-file will be opened. This only makes sense, if the keyword all was given previously.
psi Synonym for psi = single.
sheptraj Gaussian centres will be written to sheptraj. This is for vMCG quantum Grow.
speed The CPU-time used within an output interval will be written to the file speed. (default).
no-speed The speed file is not opened.
steps Information on the integrator step sizes will be written to the file steps. the file steps. For multi-layer runs with ml-cmf=split (see below), each mode will have its own steps file, named "nxx.steps" where xx is the mode number.
stop The file stop is created. It allows to stop the run in a controlled way by writing 'stop' or the desired stop time (real-time and/or CPU-time) to the stop file. (default).
no-stop The stop file is not opened.
timing Program timing information will be written to the file timing. (default).
no-timing The timing file is not opened.
update If the constant-mean-field integrator with adaptive step size is used, the update time for the mean-fields is written to the file update. (default).
no-update The update file is not opened.
veigen The eigenvectors and eigenvalues of a 1D operator, set up in the the INIT_WF section with the spf type eigenf, are written to the veigen file.
quadpop If set, the grid populations and state populations are calculated using rho^2 rather than rho. This keyword is needed, if a traceless density operator is propagated. Otherwise, all populations are zero. Only for Density Operators of Typ I.
graphviz For Multi-Layer runs, generate an input file for graphviz by which one can visualize the ML tree. Graphviz is, for example, part of the suse linux distribution. The generated file will be called "mltree.dot". If the graphviz software is installed, one can e.g. use the command "dot -Tx11 mltree.dot" to display the tree in an X window on the screen. The command "dot -Tps -o mltree.ps mltree.dot" generates the postscript file mltree.ps. Circles represent the multi-layer modes; the number inside a circle gives the mode number. Boxes represent the primitive degrees of freedom. The numbers on the edges give the numbers of single-particle functions of the submodes, or the numbers of primitive basis functions, respectively.
Keyword Default
auto S = once
cross S = name
psi S = single
psi = compact R = 1.0d-6

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.

Keyword Description
opname = S The operator with name S.op will be used.
oppath = S The path S will be used to find the operator file. If oppath is not given, the program will first look in the startup directory and then the default operator path.
closed Density operators are propagated in a closed system, i.e. possible dissipative operators are ignored and the Hamiltonian is used only. This is the default.
open Density operators are propagated in an open system.
projection Modified equations of motion for the coefficients are used if density operators of type II are propagated in an open system. This ensures that the trace is conserved. To switch off this feature use the key word no_projection.
Note: projection is the default.
alter-parameters
.....
end-alter-parameters
The lines between the keywords define parameters to be used in building the Hamiltonian, using the same format as in the PARAMETER-SECTION of the .op file.
parfile = S The parameters to be used in building the Hamiltonian are listed in the file S. The parameters are defined using the same format as in the PARAMETER-SECTION of the .op file. The file must start and end with
parameter-section
end-parameter-section
alter-labels
.....
end-alter-labels
The lines between the keywords redefine labels specified in the LABELS-SECTION of the operator file opname.op
splinepath = S PES data in a format for spline fitting is contained in the directory S. See the HOWTO on spline fits for more information.
v < R Energy cut-off used for potential energy terms. All potential energy values greater than R are set to R.
v > R Energy cut-off used for potential energy terms. All potential energy values less than R are set to R.
analytic_pes If the operator contains a non-separable potential this will not be stored explicitely on the primitive grid points, but in an analytic form which can be used to generate the potential on-the-fly at any point. This should be set if the CDVR method is to be used.
cutoff = R, unit All real diagonal Hamiltonian terms (except natpots) which are smaller than cutoff are removed. (Note, all non-diagonal Hamiltonian terms which are on all grid-points smaller than 1.d-12 au are removed as well). The default is cutoff=tiny (i.e. 1.d-9 au). For an improved relaxation run it may be useful to set cutoff to a lower value. The value of cutoff and the number of Hamiltonian terms removed are protocoled in the op.log file. NB cutoff is not applied to natural potentials. Use natpotcut for those.
natpotcut {V1,V2,...} = R, unit The real constant R sets the threshold for removing natpot terms. This feature may reduce the number of natpot terms while only marginally reducing the quality of the potfit potential. The labels V1, V2, ... are the labels assigned to a natpot in a Labels-Section of the operator file. The keyword natpotcut may appear multiple times in the input file, if different thresholds for different natural potentials are used. If no potential label is given, i.e. natpotcut = R, unit (the equivalent form: natpotcut {all} = R, unit is also possible), the program will use the same threshold for all natural potentials. The keyword unit denotes the MCTDH units. If unit is not given, au is assumed. Information on the removed natpot terms is given in the op.log file.
The default value for the threshold R is R = tiny (i.e, 1.0d-9 au). I.e., even when the natpotcut keyword is not given, all terms which are smaller then 1.0d-9 au will be removed. In fact, this holds for all operators, not only for natural potentials. Using an increased value for the natpotcut threshold (e.g. R=1.d-6) may speed up the calculation, because several natpot terms have been removed. The accuracy of the potentials may decrease, but this effect is negligible, as long as R is sufficiently small.
print-npot Full information on natpots is printed to op.log file. If print-npot is not given, all npot lines except the first and the last one of each natpot are suppessed.
optcontract The Quantics program sums terms of the same type in the operator during the build process. For example real potential terms will be summed separately from matrix operators in a DVR, etc. This keeps PE and KE terms separate. If the optcontract keyword is given, the program will try and contract more terms together, e.g. real matrices and diagonal matrices can be summed. This mixes PE and KE terms.
gwpintorder = I Integrals of GWP matrix elements are made to order I. Default I = 2 (i.e. use the LHA).
gwp_renorma Renormalise wavepacket by changing the A-vector at end of each step. This is the default unless a CAP is present or it is an independent GWP calculation.
nogwp_renorma Do not renormalise wavepacket by changing the A-vector at end of each step.

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.

Keyword Description
qcprogram = S The quantum chemistry program S will be used. S = Gaussian, Molpro, Molcas
subcmd = S [,2] The quantum chemistry calculations are started using the command (or script) S. E.g. if S = qsub then jobs are started using qsub jobname (jobname is the stem of the qc input file). If the option 2 is added the submission command used is qsub jobname dir (dir is the directory containing the input file)
simulate The quantum chemistry program is not run but output files read and processed. This is useful for creating a DB from a set of pre-run calculations.
method = S The ab initio method used. This enables the progarm to know the format of the data it is looking for. S = CAS, RHF
ener0 = R The value of the zero point of energy.
nroot = I The number of roots to be calculated.
fchk If set, Gaussian will create formatted checkpoint files.
data = S The data from a direct dynamics run is stored in directory S. The default is name/dd_data.
db = S The string S controls the use of the database.
S = rdwr write to and read from database (default)
S = read do not write to database. I.e. do not perform further QC calculations and take PES from DB.
S = only synonym for read
S = write QC calculations performed at each step and results stored in database.
S = none QC calculations performed at each step but not saved.
ascii_db At the end of the propagation the DB will be converted to ascii.
update = S The string S controls the updating of the potential surface. Between update times the evolution takes place by extrapolating across the stored PES. Values for S are:
always PES updated every step
never PES never updated
tout PES updated at time tout (set in RUN-SECTION)
error Update time controlled by an error criterion. error,R
dbmin = R If geometry is less than R distance from a point in the DB then the surface is taken from the data at the DB point.
dbpoint = I Specify the distance measure used to calculate a new point.
I=1 (default) RMSD in Cartesian coordinates
I=2 Max atom displacement
I=3 RMSD in normal modes
dbweight = I Specify the distance measure used in the DB Shepard Interpolation.
I=1 (default) RMSD
I=2 Max atom displacement
mindv = R Sets an energy gap threshold in Hartree below which the data calculated at this point will not be used in the dynamics. This parameter is used to exclude points where the diabatization procedure is expected to give wrong results. Points with adiabatic energy gaps bellow R are calculated and stored on the database but are not used in the dynamics calculation.
dd_adiab Dynamics run in the adiabatic picture (default for single surface case, redundant otherwise)
dd_diab = S The string S controls the diabatisation of the PES Values for S are:
regularized local global
ddtrans = I The dynamic coordinates must be transformed to the I cartesian coordinates of the PES. The Transformation matrices may either be read from a DDTRANS-SECTION or from a file specified by transfile
transfile = S [,e0read, S1] The file S is the output from a quantum chemistry calculation containing the normal mode coordinates that transform between system and PES coordinates. If the e0read argument is added the energy from this file is the zero point energy. S1 specifies the type of quantum chemistry used.
coinfile = S , S1 The file S is the output from a quantum chemistry calculation containing the conical intersection point, derivative coupling and gradient difference vectors. These are required for the regularized diabatisation scheme. Alternatively this information can be read from the input file in a DDPOINTS-SECTION.
dreffile = S The file S is the output from a quantum chemistry calculation containing the reference point (usually the Franck-Condon geometry) required for the regularized diabatisation scheme. Alternatively this information can be read from the input file in a DDPOINTS-SECTION.
dipole Dipoles will also be stored in dip.db.
nbasis=I Number of basis functions used in the quantum chemistry calculation. This is used because for some calculations the molecular orbitals' coefficients are stored in file moco.db.
dbsave Save the DB in memory during the propagation rather than reading from file.
hess_upd Use the Powell method to update the Hessian rather than recalculating at each DB point.
rediabatise Used in conjunction with the global diabatisation option to re-diabatise the PES using a pre-computed DB. This diabatisation is done from the central point outwards rather than along the GWP paths. Needs db = rdwr and it is a run-type of genoper.

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,

spf-basis-section
multi-set
    nfuncs = 3 , 2
end-spf-basis-section

This is useful in a vMCG calculation with many degrees of freedom.

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

Keyword Description Default
multi-set If an electronic basis is defined it is treated using the multi-set formalism, i.e. a set of single-particle functions per state. If this keyword is not present, the single-set formalism is used. not set
single-set If an electronic basis is defined it is treated using the single-set formalism, i.e. a common set of single-particle functions for all electronic states. This keyword need not to be given, single-set is default for electronic states. not set
packets = I I independent wavepackets will be simultaneously propagated. Technically, the packets are propagated on auxiliary electronic states. One may choose between S=single-set or S=multi-set. (multi-set is default here). For block improved relaxation, single-set must be given. A block improved relaxation requires the packet keyword and a DAV, RDAV or RRDAV "integrator". I = 1, multi-set
id,I I denotes the mode-number (particle-number) with which the present mode is to be identified. See the note above for the correct use of the id keyword. not set
no-redundancy-check If this keyword is given, the check on redundant SPFs is suppressed. Use with care. not set

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  
indicates that the top layer supports three particles, each of which is expanded into 5 SPFs. Hence this line must be followed by three 1> lines, and so on. The tree is terminated when the primitive basis (grid) is reached. The primitive basis is indicated by [..], e.g.
   2> [q1]
 
indicates that a SPF of a third layer is represented by the primitive basis (grid) with mode-label q1. If mode combination is used, one would e.g. write
   2> [q1,q2]
 

See the paper:
O. Vendrell and H.-D. Meyer
Multilayer multiconfiguration time-dependent Hartree method:
Implementation and applications to a Henon-Heiles Hamiltonian and to pyrazine.

J.Chem.Phys. 134 (2011), 044135.


Example: ML-tree for 6D, 3-layers Henon-Heiles
 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
   
For a graphical representation of the tree click here .


Example: ML-tree for 24D, 6-layers pyrazine
 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
   
For a graphical representation of the tree click here .

A graphical representation of the entered tree can be generated by setting the graphviz keyword in the Run-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

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

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

Parameter Description
el Electronic basis.
elcont Electronic basis including continuum states.
HO Harmonic oscillator (Hermite) DVR.
rHO Radial Harmonic oscillator (odd-Hermite) DVR.
Leg Rotator (Legendre) DVR.
Leg/R Rotator (Legendre) DVR for a restricted range on angles.
Lagu1 Laguerre DVR for boundary condition x1/2.
Lagu2 Laguerre DVR for boundary condition x1.
Lagu3 Laguerre DVR for boundary condition x3/2.
Lagu4 Laguerre DVR for boundary condition x2.
sin Sine (Chebyshev) DVR.
FFT Fast Fourier transform collocation.
exp Exponential DVR. Periodic boundary condition.
cos Cos DVR. "gerage" solutions with periodic boundary condition.
sphFBR Spherical harmonics FBR (deprecated).
KLeg Extended Legendre DVR.
K K-quantum number appearing with KLeg-DVR.
PLeg Two-Dimensional Legendre DVR.
Wigner Wigner-d (three-dimensional rotor) DVR.
Extern External DVR.
GWP Gaussian Wavepacket Basis. Basis_size is not required (see Remark below).
GWPCL A classical Gaussian Wavepacket Basis, i.e. the GWPs follow classical trajectories. Basis_size is not required (see Remark below).

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 = 2m3n 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 or GWPCL 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)

Basis Parameters Parameter description
Elcont Nstates: Total no. of electronic states.
Elbnd: No. of bound electronic states.
HO hoxeq: Equilibrium position of harmonic oscillator basis functions.
hofreq: Frequency of harmonic oscillator basis functions.
homass: Mass of harmonic oscillator basis functions. If no mass is given, then the mass is set to 1.
HO S: String S = xi-xf. This string serves as a switch between the two possible input formats.
xi: First grid point.
xf: Last grid point.
rHO hoxeq: Equilibrium position of harmonic oscillator basis functions, which -- because only the positive half-axis is used -- is the left-hand boundary of the wavefunction.
hofreq: Frequency of harmonic oscillator basis functions.
homass: Mass of harmonic oscillator basis functions. If no mass is given, then the mass is set to 1.
rHO S: String S = xi-xf or S = x0-xf. This string serves as a switch between the possible input formats.
xi or x0: First grid point, or left boundary (i.e. hoxeq).
xf: Last grid point.
Leg blz: Magnetic rotational quantum number. Alternatively to a number one may input the string jbfXXX. blz will then be set to the value of the parameter jbfXXX. Here XXX stand for any characters. I.e. one may use jbf or jbf_1 etc. Note: jbfXXX must be defined in the parameter section of the input-file or via an option on the command line. (alter-parameter and operator file definitions come too late).
string: Controls whether symmetry to be used.
all: no symmetry (use all l -values, l=m, m+1,...,m+N-1 ).
odd: odd symmetry (i.e. l=odd ).
even: even symmetry (i.e. l=even ).
Leg/R blz: see Leg
string: see Leg
theta1 Lowest value of theta (in rad) to be included in the restricted grid.
theta2 Largest value of theta (in rad) to be included in the restricted grid.
Lagu1 x0: Sarting point of the radial interval (usually zero).
b: Length parameter. Chin(x) = 1/b * Sqrt((x-x0)/n) * exp(-(x-x0)/(2*b)) * L1n-1 ((x-x0)/b)
icut: Cut parameter to avoid excessively large kinetic energy contributions. icut=0 leaves the second derivative matrix unmodified. See remarks below
Lagu1 S: String S = xi-xf. This string serves as a switch between the three possible input formats.
xi: First grid point.
xf: Last grid point.
icut: See icut above. See remarks below.
Lagu1 S: String S = x0-xf. This string serves as a switch between the three possible input formats.
x0: Sarting point of the radial interval (usually zero).
xf: Last grid point.
icut: See icut above. See remarks below.
Lagu2 --- Input identical to Lagu1
Lagu3 --- Input identical to Lagu1
Lagu4 --- Input identical to Lagu1
sin xi: First grid point.
xf: Last grid point.
string: short , long (short is the default), and/or sdq , or spin . See note below.
sin string: 2pi or 2pi/m , where m is a positive integer (its default is 1) denoting the multiplicity (e.g. of the rotational axis). The wavefunction is assumed to be periodic on the interval -2pi/m to 2pi/m.
Because only ungerade (asymmetric) wavefunctions are computed, the grid is halved and only grid points for positive x appear.
string: sdq . If sdq is set, the symmetrized first derivative, (sin*d/dx+d/dx*sin)/2 is used rather than the simple first derivative.
fft or exp xi: First grid point.
xf: Last grid point.
string: linear, periodic or s-periodic (linear is the default).
fft or exp string: 2pi , s-2pi , c-2pi or 2pi/m , s-2pi/m , c-2pi/m , where m is a positive integer (its default is 1) denoting the multiplicity (e.g. of the rotational axis). The grid is assumed to be periodic, ranging from 0 to 2pi/m. (See note below).
For an exp-DVR which follows PLeg, the input k= kmin,kmax may follow. The default is kmax=-kmin=(N-1)/2. (See note below,PLeg).
exp
(If combined
with PLeg)
k=kmin,kmax
(optional)
The range of magnetic rotational quantum number of the PLeg. If kmax is omited (k=kmin), the range [-kmin,kmin] is chosen. Default is the full range [-(gdim-1)/2, (gdim-1)/2] (or [-jtot,jtot] if jtot is set).
cos xi: First grid point.
xf: Last grid point.
string: short or long (short is default). See note below.
cos string: 2pi or 2pi/m , where m is a positive integer (its default is 1) denoting the multiplicity (e.g. of the rotational axis). The wavefunction is assumed to be periodic on the interval -2pi/m to 2pi/m.
Because only gerade (symmetric) wavefunctions are computed, the grid is halved and only grid points for positive x appear.
sphFBR jmax: Maximum value of quantum number j of the spherical harmonics basis functions. Note: jmax replaces N, the number of basis-functions/grid-points. N is computed from the input. See log file.
string: nosym: no symmetry (uses all values of j below jmax, j=0,1,...jmax).
sym : symmetry (uses values of j of the same parity as jmax, i.e. jmax-j = even).
thrshld: Threshold for convergence when the FBR integrals are performed. This input is optional, default: thrshld=1.d-10.
j_off: Offset value used when the FBR integrals are performed. The first iteration uses j_max + j_off quadrature points. This input is optional, default: j_off=6.
phiFBR
(Must follow
sphFBR)
mmax:
(optional)
Maximum value of quantum number m of the spherical harmonics basis functions. Uses only values of m, such as |m| < = mmax and |m| <= j. Default is mmax=jmax.
mincr:
(optional)
Increment of m's, starting from mmax (must be given: mmax mincr).
Default is 1.
KLeg string: Controls whether symmetry to be used.
all: no symmetry (use all l-values).
odd: odd symmetry (i.e. l = odd).
even: even symmetry (i.e. l = even).
K
(Must follow
KLeg or Wigner)
kmin: Minimum value of body fixed magnetic quantum number of basis functions.
kmax: Maximum value of body fixed magnetic quantum number of basis functions.
dk: (optional) Delta K. Increment in K-value. Default dk=1.
PLeg string: Controls whether symmetry to be used.
all: no symmetry (use all l-values).
odd: odd symmetry (i.e. l = odd).
even: even symmetry (i.e. l = even).
Wigner string: Controls whether symmetry to be used.
all: no symmetry (use all j-values).
odd: odd symmetry (i.e. j = odd (NOT YET IMPLEMENTED)).
even: even symmetry (i.e. j = even (NOT YET IMPLEMENTED)).
Extern string: Name of file containing grid points and DVR derivative matrices. The file may contain only grid points. Then the derivative matrices will be zeroed. (See remark below).
string: ascii: the file is read in ascii format (default).
binary: the file is read in binary format.
unit: The (optional) string unit is a length unit (e.g. angst, nm, pm, or deg) with which the input data is multiplied.
GWP xi: First grid point.
xf: Last grid point.
GWPCL xi: First grid point.
xf: Last grid point.

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

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

For the 2pi and sdq keywords see remarks on cosine DVR. Note that cos-DVR always uses sdq, whereas for sin-DVR one must explicitly give the sdq keyword. Otherwise dq will be assumed. If the sdq keyword is set the operators dq, qdq, and cdq are not allowed for sin-DVR.

The keyword spin is only allowed for gdim=2. If set, it produces 2x2 derivative matrices with zero diagonal and (1,-1) (first derivative) or (1,1) (second derivative) off-diagonal elements (i.e. i*sigma_y and sigma_x where sigma denotes the Pauli matrices).

Remarks on cosine DVR:

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

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.

Example:
PRIMITIVE-BASIS-SECTION
  beta    wigner  20   all
  gamma   exp     15   2pi
  alpha   k       -7   7
end-primitive-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

Keyword Description
nstates = I I electronic states are included in the problem. Default = 1
init_state = I Initial electronic state of a wavefunction propagation.
mom_distrib Initial distribution of GBFs is made in momentum space (default).
coo_distrib Initial distribution of GBFs is made in coordinate space.
init_basis_dist = S Selects different schemes to place basis functions initially. Possible values are:
seqex
Sequentially excite/move each degree of freedom with respect to the centre leaving other degrees of freedom untouched.
simex
Simultaneously excite/move all degrees of freedom with respect to the centre. (All basis functions end up along a diagonal in coordinate space?)
shell (Default)
Form shells around the centre position by initially exciting/moving each degree of freedom sequentially, and then combining motions on several degrees of freedom.
shell-noq0
Same as shell above but do not place any basis function at the centre.
wigner
Place basis functions randomly sampling a Wigner distribution is positions and momenta. (The first basis function is place at the centre.)
random
Same as wigner above but only positions are sampled (each function has no initial momentum).
wigner_sym
Same as wigner above sampling is symmetric in phase space.
shell-diag
Similar to shell above but place the functions along diagonal instead of along the axis of coordinate space, i.e. do at least one excitation/motion along each coordinate. (This is useful to avoid too symmetric initial conditions.)

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, either Nmode and Atom blocks or Cartesian and Freq blocks are required.

  Cartesian = S
  ATNAM     X     Y    Z     [ keywords ]
  End-Cartesian
  
ATNAM a label that needs to start with the atom type in usual chemical notation, examples are O3, H22, Oend, Nterm . Note that first 2 characters mmay form a name accidentally, e.g. Nend would be neon rather than the "end" nitrogen. isotopes D and H(iso = 2) are known.
X,Y,Z Cartesian coordinates in unit specified by unit keyword.
cartesian = S The string S defines the length unit. Default is a.u.
Optional Keywords
freeze freeze atom
classical atom is treated classically
grid N xi xf optional keywords to specifiy the auxiliary grid used for analysis (see primitive-basis-section for more details).
  Nmode
  Label     Q  Freq [, unit] [keywords]
  End-nmode
  
ATNAM a label that needs to start with the atom type in usual chemical notation, examples are O3, H22, Oend, Nterm . Note that first 2 characters mmay form a name accidentally, e.g. Nend would be neon rather than the "end" nitrogen. isotopes D and H(iso = 2) are known.
X,Y,Z Cartesian coordinates in unit sepcified by unit keyword.
N xi xf optional keywords to specifiy the auxilairy grid used for analysis (see primitive-basis-section for more details).
Optional Keywords
freeze freeze normal mode at Q
classical atom is treated classically
grid N xi xf optional keywords to specifiy the auxiliary grid used for analysis (see primitive-basis-section for more details).

INIT_WF-SECTION

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

Keyword Description
file = S1 (,S2,S3) The initial wavefunction will be read from the restart file in directory S1. If S1 is not specified, the restart file will be taken from the name-directory. For a density operator II propagation, the restart file may contain a wavefunction. It will then automatically be transformed to a pure-state density operator of type II.
S2 = orthopsi : The single-particle functions of the initial wavefunction are transformed to natural orbitals, Schmidt-orthogonalised and transformed back after being read from file (this is the default).
S2 = noorthopsi : The initial wavefunction is not Schmidt-orthogonalised after being read.
S2 = realpsi : The real part of the wavefunction will be Schmidt-orthogonalised as in the orthopsi-case and used for the initial wavefunction.
S3 = ignore : Ignore that primitive bases are different. This is a very dangerous option, because when the grids do not match, your results will be wrong. However, it allows you to use a restart file which is, say, generated with sin-DVR, whereas you want to use FFT during the propagation.
Note: The numbers of SPFs of the current WF may differ from those of the WF read in. However, if the current WF has less SPFs than the one read in, some SPfs are removed which may lead to strange results, in particular if noorthopsi is set.
build
.....
end-build
The initial wavefunction will be build using the data specified between these keywords. See Building the initial wavepacket.
read-inwf
.....
end-read-inwf
The initial wavefunction will be read from one or several restart files. The SPFs (and/or the coefficients) for different electronic states can be read from different restart files. This is a more powerful tool than the one inwoked through the file keyword discussed above. See: Reading the initial wavepacket.
block-SPF = S1 (,S2) The SPFs for the initial wavefunction will be read from the restart file in directory S1. (Alternatively, S1 may denote the path of the restart file). The restart file S1 may be a multi-packet or a single-packet one. The SPFs of a packet basis are ignored and the packet basis of the current WF is newly build. Note that this is for multi-packets single-set runs (e.g. block-improved relaxation) only. For other runs one should use read-inwf.
Note that the A-vector is not yet defined, it must be defined via A-coeff or autoblock or block-A.
S2 = noorthopsi : The initial wavefunction is not Schmidt-orthogonalised after being read.
S2 = orthopsi : The initial wavefunction is Schmidt-orthogonalised after being read (default).
S2 = realpsi : The only the real part of the SPFs will taken and the thus modified SPFs are Schmidt-orthogonalised as in the orthopsi-case and used to build the initial WF. Note that if realpsi is set, it will force the program to take the real part of the A-vectors in case they are read in by a following block-A command. In case block-SPF is used in combination with block-A one MUST set noorthopsi because before orthonormalization the SPFs are transformed to natural orbitals and hence are no longer consistent with the A-vectors read in. See: Generating block initial wavepackets.
block-A = S1 (,S2,S3,..) The A-vectors of a multi-packets single-set initial wavefunction will be read from restart files S1, S2, S3, .... The restart files must be of standard non-packet form. The number of restart files must match exactly the block size. One may give more than one block-A keyword.
The SPFs may be generated by build, or -- more likely -- may be read from another restart file with the aid of the block-SPF keyword. In case block-SPF is used in combination with block-A one MUST set noorthopsi in block-SPF. If realpsi is set in block-SPF it will affect the A-vectors read by block-A as well. See: Generating block initial wavepackets.
realpsi Only the real part of the WF will taken and the thus modified SPFs are Schmidt-orthogonalized and the A-vector is re-normalized. Note that the keyword realpsi accomplishes to take the real part of the WF at the end of geninwf, i.e. after operate, orthogonalize, etc., whereas the argument realpsi given for file, read-inwf, or block-SPF does so only for the WF read. The realpsi feature is useful for improved relaxation in real version (RDAV and RRDAV).

After an initial WF 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:

A-coeff, meigenf, operate, orthogonalise, correction, &nbsprealpsi,

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.

WARNING FORTRAN cannot open a file twice. Hence it may become neccessary to copy a file such that the same information can be read via two different file-paths. Such a problem may arrise 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.

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

Keyword Description
A-coeff
.....
end-A-coeff
The lines between the keywords contain in free format one integer (A-index) and one complex number (A-value) per line, thus defining the initial A-vector. The default is A(1) = (1.0,0.0) and zero for all other coefficients.
For a multi-set run, the input line must contain two integers. The first is the A-index and the second the electronic state. A multi-set density-operator run requires three integers before the complex number (A-value), the A-index and the two electronic state indices s and t. E.g.:
518 (0.717,0.0) (single-set, wavefunction and densities)
518 1 (0.717,0.0) (multi-set, wavefunction)
518 1 1 (0.717,0.0) (multi-set, densities)

There is also a long form of the A-coeff input, which however is only allowed for nmode < 33 (nmode < 17 for density type II). Here the number of the particle for each mode (and the electronic state(s) in case of a multi-set run) is input followed by the value of the A-coefficient. E.g.:
1 4 1 3 1 2 (0.5,0.0)
would be a correct input for either a 6-mode single-set, or a 5-mode multi-set WF-calculation. In the latter case the 2 would specify the electronic state. For densities type I this would be a correct input for either a 6-mode single-set, or a or a 4-mode multi-set calculation. In the latter case the last two indices 1 2 would specify the electronic states. For densities type II this would be a correct input for either a 3-mode single-set (j1, j2,j3,k1,k2,k 3), or a or a 2-mode multi-set calculation (j1, j2,k1,k2,s,t).
This re-definition of the A-vector is also possible, when the initial wavefunction is read in (file keyword, restart run). The A-vector is re-normalised before being used. The non-zero entries of the A-vector are protocolled in the log-file.
autoblock autoblock automatically generates a set of A-vectors for a block improved relaxation run. This set is suitable for converging to the nb lowest eigenenergies, where nb denotes the block size, i.e. the number of packets.
select_state = I Read in WF projected on to Ith state.
meigenf = I,S, I1|S1(,S2,S3,S4,I2) Generate the mode-eigenfunctions using operator S, where S = XX is the name of an operator, defined in the HAMILTONIAN-SECTION_XX. Only that uncorrelated part of the operator S, which acts on the mode I, will be used. The operator must be generated with the usediag keyword. I1 denotes the number of the eigenstate (counting from zero), which is taken as the first SPF. The ordering and the eigenenergies are protocoled in the log file. The number I1 may be replaced by the string S1 = follow. In this case the eigenstate with the largest overlap with the starting vector (usually defined in a build block) is taken as first SPF. If the optional number I2 is given, the maximal Lanczos space is restricted to I2. Otherwise the maximal Lanczos space size is equal to the length of the mode vector (i.e. subdim). If the optional string S2=full is given, the number of Lanczos iterations is limited by I2 (or subdim) only. Otherwise the Lanczos iteration is stopped, when the chosen vector (i.e. I1) is converged. If the optional string S3=select or S3=noselect is given, states with zero overlap (<10-13) will be ignored/not ignored when mapping the eigenvectors to the psi function. When S3 is not given, then noselect is chosen by default when full is given and I2=subdim (e.g. I2 not given). Otherwise the default is select. If the optional string S4=write is given, the eigenvectors will be written to the (ASCII) file meigen_mode_ state.
Examples:
meigenf = 3,oper,0
meigenf = 3,oper,follow,full,125
meigenf = 3,oper,follow,full,select,write,125
operate = S(,S1(,S2,...)) Operate on the initial wavefunction using operator S, where S = XX is the name of an operator, defined in the HAMILTONIAN-SECTION_XX. One may apply up to 15 operators to the WF. operate = O1, O2, O3 will produce O3*O2*O1*|psi>. The wavefunction is re-normalised after each application of an operator. The normalisation factors are protocoled in the log-file.
For density operators it calculates the commutator -i [O1,rho0] .
operate = O1,O2, ... ,On will produce (-i)n [On,...[O2,[O1,rho0]]...] .
operate_iter = I Maximum no. of iterations to be used in applying an operator to an MCTDH wavefunction. Default = 10. Only the A-vector but not the SPFs will be modified, if operate_iter=0 is given. This can be useful when generating an initial WF for improved relaxation.
operate_tol = R Convergence tolerance to be used in applying an operator to an MCTDH wavefunction. Default = 1.0d-8
operate_no-norm The normalisation factor is removed from operator*psi, i.e. the initial wavefunction is no longer normalised.
operate_no-direct A non-iterative algorithm may be applied first before the iterative improvement the wavefunction. With the keyword operate_no-direct, the non-iterative part is skipped. (operate_no-direct is default.)
operate_direct With the keyword operate_direct, first the non-iterative algorithm is applied before the iterative procedure starts (only for nstate=1).
orthogonalisation = S (,S1 (,S2 ..)) S = path_to_foreign_restart_file_or_directory
The initial WF is orthogonalised against the WF on the specified file. One may give several files in one orthogonalisation statement, or one may give more than one of such a statement, in order to orthogonalise against several wavefunctions. This feature is useful in particular for improved relaxation. If the path points to a directory, /psi is automatically appended.
symorb = I1,I2 I1 = Number of first set of SPFs (mode number), I2 = Number of second set of SPFs.
With this keyword one mixes two sets of SPFs. The new order (in both sets) is now: 1st SPF of 1st set, 1st SPF of 2nd set, 2nd SPF of 1st set, 2nd SPF of 2nd set, etc. The SPFs were orthonormalized and those which tend to be linearly dependent are removed. As the two sets of SPFs are now identical, one may define symmetric and anti-symmetric linear combinations with the aid of A-coeff.
Note that the id keyword must not be set in the SPF-Basis-Section. One may, however, create an initial restart file with a geninwf run (without id) and then read this restart file (file=..) in a subsequent propagation run using the id keyword.
sym1d = I ((,S),I1(,S),..)

asym1d = I ((,S),I1(,S),..)
I = Number of DOF which is to be (a)symmetrized. S = persist
phi(i) = phi(gdim+1-i), i=1,...,gdim, for symmetrization and phi(i) = -phi(gdim+1-i) for asymmetrization.
In order not to annihilate a function by (a)symmetrization one should use in the build block the types: HO odd, HO even, gauss odd, gauss even, or sym as argument to pop when the type eigenf is used.
The (a)symmetrization is applied to the initial wavefunction only. However, if the keyword persist is given and when an improved relaxation run is performed, then the (a)symmetrization is done additionally after each orbital relaxation.
parity = I ((,S),I1(,S),..) I = Number DOF which is to be (a)symmetrized. S = persist
Similar to sym1d and asym1d. Functions which are predominantly symmetric are symmetrized, and those predominantly anti-symmetric are anti-symmetrized.
The (a)symmetrization is applied to the initial wavefunction only. However, if the keyword persist is given and when an improved relaxation run is performed, then the (a)symmetrization is done additionally after each orbital relaxation.
sym2d = I ((,S),I1(,S),..)

asym2d = I ((,S),I1(,S),..)
I = Number 2D-mode which is to be (a)symmetrized. S = persist
phi(x,y) = phi(y,x), for symmetrization and phi(x,y) = -phi(y,x) for asymmetrization. The mode(s) specified must be 2D combined modes and the primitive grids within these modes must be identical.
The (a)symmetrization is applied to the initial wavefunction only. However, if the keyword persist is given and when an improved relaxation run is performed, then the (a)symmetrization is done additionally after each orbital relaxation.
nsym2d = I ((,S),I1(,S),..) I = Number of the 2D-mode for which the symmetric and antisymmetric combinations of the spfs have to be made. S = persist
The mode(s) specified must be 2D combined modes and the primitive grids within these modes must be identical. NOTE: This also changes the A-vector! (see the example input "H2H2_nsym.inp" in the "Further Example inputs" section)
sym2kleg = I (,I1,...)

asym2kleg = I (,I1,...)
I = Number of 4D-mode (2x KLeg/K) which is to be (a)symmetrized.
phi(θ1,k12,k2) = phi(θ2,k21,k1) for symmetrization and phi(θ1,k12,k2) = −phi(θ2,k21,k1) for asymmetrization. The mode(s) specified must be 4D combined modes of the form KLeg/K/KLeg/K, and the primitive grids for the two KLeg and K DOFs must be identical.
sym3d = I ((,S),I1(,S),..) I = Number 3D-mode which is to be symmetrized. S = persist
phi(x,y,z) = phi(y,z,x) = phi(z,x,y) = phi(y,x,z) = phi(x,z,y) = phi(z,y,x). The mode(s) specified must be 3D combined modes and the primitive grids within these modes must be identical. The symmetrization is applied to the initial wavefunction only. However, if the argument persist is given and when an improved relaxation run is performed, then the symmetrization is done additionally after each orbital relaxation.
symcoeff ( = S,S1,I1,I2,...) S, S1 = persist, dav
I1,I2,... = mode numbers for partial symmetrization.
All arguments are optional. The A-vector is symmetrised by summing all permutations of its indices. This makes only sense, if all modes are identical, i.e. if the id keyword (SPF-Basis-Section) is used. The symmetrization is applied to the initial wavefunction only. However, if the argument persist is given and when an improved relaxation run is performed, then the symmetrization is done additionally after each orbital relaxation. When the argument dav is given, a symmetrization of each Davidson vector is performed in addition. The argument dav should only be given for relaxation runs using the RDAV or RRDAV "integrator", otherwise dav is ignored.
If integers I1, I2, .. are given as argument, only a sub-set of modes (particles) will be symmetrized, (useful for mixtures). E.g.: if symcoeff=2,3,4 is given, only particle 2, 3, and 4 are symmetrized. In such a case there may be several (a)symcoeff statements, as long as they refer to different modes (particles).
asymcoeff ( = S,S1,I1,I2,...) S, S1 = persist, dav
I1,I2,... = mode numbers for partial symmetrization.
The A-vector is asymmetrised by summing sig(P)*A(P(j1,..,jp)), where P denotes a permutation. This makes only sense, if all modes are identical, i.e. if the id keyword (SPF-Basis-Section) is used for the modes to be anti-symmetrized.
For a description of the arguments, see symcoeff
correction = S1,S2 S1,S2 = edstr, dia, ad, hh2.
edstr: compute the (uncorrected) energy distribution for the flux analysis.
dia: compute the diabatically corrected energy distribution.
ad: adiabatic correction.
hh2: use the routines specially written for H+H2 reactive scattering using LSTH PES.
hh2-bkmp2: use the routines specially written for H+H2 reactive scattering using BKMP2 PES.
Note that the translational degree of freedom has to be specified by the trans keyword if it is not the first DOF.
trans = I1(,I2) I1,I2 = dof and state of the translational mode. (Needed for correction).
If the keyword trans is not given, dof=1, state=1 is assumed when computing the correction.
tfac = R R = partition factor for Jacobian coordinates. R = m1/(m1+m2) where m1 and m2 are the masses of the two atoms of the diatomic molecule. If not given, tfac=0.5 is assumed.
s-mctdh = R An S-MCTDH wavefunction will be generated with a cutoff for selecting the configurations of R. (see Eq.(21) of JCP 112, 8322 (2000)). See note below.
Special keywords for build sub-section
Keyword Description
init_state = I Initial electronic state of a wavefunction propagation.
left_state = I Initial electronic state for the ket part of a density operator.
right_state = I Initial electronic state for the bra part of a density operator.
temperature=R Initial temperature for a density operator propagation.
print Print some information to the log file on how multi-dimensional mode-functions are build from 1D-functions. NB Not for KLeg.

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

Building the initial wavefunction
Type Description
HO Harmonic oscillator eigenfunction.
HO odd Odd harmonic oscillator eigenfunctions (n=1, 3, 5 ...).
HO even Even harmonic oscillator eigenfunctions (n=0, 2, 4, ....).
gauss Gaussian wavepacket.
gauss odd Gaussian wavepacket. Odd functions only
gauss even Gaussian wavepacket. Even functions only.
Leg Legendre polynomial.
sphFBR Spherical harmonics.
phiFBR Indicates second coordinate of sphFBR.
eigenf Specified potential eigenfunction.
KLeg Associated Legendre polynomial. Requires KLeg or Pleg in Primitive-Basis-Section.
K Body fixed magnetic quantum number for KLeg (or PLeg).
Wigner Wigner-d function. Requires Wigner in Primitive-Basis-Section.
map Use the initial (1D) single-particle functions of another DOF.
tdpoint The DOF is represented by a single point that is driven by a time-dependent function. It must not be part of a combined mode.
gwp A Gaussian wave packet for a DOF with a type gwp (vMCG and G-MCTDH)
gwpp A Gaussian wave packet for a DOF with a type gwp (vMCG and G-MCTDH). The initial distribution of GWPs is made in momentum space in contract to the keyword "gwp" which distributes the GWPs in coordinate space.
ho-gwp Something else that someone should document.

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

Type Parameters Parameter Description
HO

HO odd

HO even
centre Centre of oscillator potential.
momentum Initial momentum of wavepacket.
frequency Frequency of harmonic oscillator. The frequency may be complex.
mass Mass of harmonic oscillator.
pop = p The pth single-particle function will be populated initially. This parameter is optional, the default is p = 1.
pack = p The Build-line is for the pth packet in a multi-packet run. This parameter is optional, the default is p = 1.
periodic The grid is assumed to be periodic. This keyword is meaningful only if the primitive representation is FFT or exp-DVR and if the primitive grid is periodic (one of the keywords: periodic, s-periodic, 2pi or s-2pi must be given in the primitive basis section). If one places a gaussian at the origin of a 2pi grid, then only half of the gaussian is taken as single-particle function. With the aid of the keyword periodic also the region near 2pi gains intensity.
gauss

gauss odd

gauss even
center Center of initial Gaussian wavepacket.
momentum Initial momentum of wavepacket.
width Denotes width of initial Gaussian wavepacket. The width here is defined as the standard deviation of the initial Gaussian, i.e. Sqrt(<x2>-<x>2). The width parameter may be complex.
pop = p The pth single-particle function will be populated initially. This parameter is optional, the default is p = 1.
pack = p The Build-line is for the pth packet in a multi-packet run. This parameter is optional, the default is p = 1.
periodic Same as for HO.
Leg m Magnetic rotational quantum number. Alternatively to a number one may input the string jbfXXX or slzXXX or blz. m will then be set to the value of the parameter jbfXXX or slzXXX, or to the value of the magnetic rotational quantum number, blz, as defined in the primitive basis set. Here XXX stand for any characters. I.e. one may use slz or slz_1 etc. Note: if jbfXXX or slzXXX is used, they must be defined in the parameter section of the input-file or via an option on the command line. (operator file definitions come too late).
l The l-quantum number of the first single-particle function. Alternatively to a number one may input the string sl0XXX. l will then be set to the value of the parameter sl0XXX. Here XXX stand for any characters. I.e. one may use sl0 or sl0_1 etc. Note: if sl0XXX is used, it must be defined in the parameter section of the input-file or via an option on the command line. (operator file definitions come too late). Note: l >= m.
symmetry Indicates symmetry.
nosym: no symmetry (all l values are used).
sym: symmetry used (gerade if l is gerade and vice versa).
sphFBR j Quantum number j of the initial spherical harmonic.
phiFBR m Quantum number m of the initial spherical harmonic.
eigenf potential Name of one-dimensional potential curve, XX, defined in the HAMILTONIAN-SECTION_XX. (See note below).
pop = p (,S(,S1)) The pth eigenfunction will be populated initially. This parameter is optional, the default is the lowest, p = 1. Note : p=1 -> ground state, p=2 -> first excited state, etc. If additionally S=sym is given, only every second eigenstate is taken. I.e. for a symmetric potential pop=1,sym selects even, and pop=2,sym odd states as initial single particle functions. If additionally S1=check is given, the program tests for the correct symmetry. This is useful if the symmetric potential is a double well, such that the odd/even character of the eigenfunctions do not follow a simple alternating pattern. The selection of the eigenstates is protocoled in the log-file.
NB: eigenf works with real and complex potentials. The Hamilton to be diagonalized, however, must be represented by a matrix. Hence eigenf does not work for FFT-DVR. If you want the eigenvectors printed, specify "veigen" in the RUN-SECTION.
KLeg l Initial rotational quantum number. Alternatively to a number one may input the string sl0XXX. l will then be set to the value of the parameter sl0XXX. Here XXX stand for any characters. I.e. one may use sl0 or sl0_1 etc. Note: if sl0XXX is used, it must be defined in the parameter section of the input-file or via an option on the command line. (operator file definitions come too late).
symmetry Indicates symmetry.
nosym: no symmetry (all l values are used.)
sym: symmetry used (gerade if l is gerade and vice versa).
excite=s Excitation algorithm for the generation of (unoccupied) spf's.
s=j: Excitation of j-states is preferred (default).
s=m: Excitation of m-states is preferred (j goes down).
s=m-old: Excitation of m-states is preferred (j goes up; old behaviour, deprecated).
print Optional. Prints (j,m) quantum numbers of the generated spf's to the 'log'-file. (This option may tell you the difference between the excite=s options.)
nspf=I Optional. Only used for multi-KLeg modes (i.e. modes with more than one KLeg). This sets the number of 1D SPFs which are initially generated; the real (multi-D) mode SPFs are later built from these 1D SPFs. If omitted, or set to zero, an automatic value is chosen. Use the print option to see how many (and which) 1D SPFs are generated.
K k Body fixed magnetic quantum number of initial wavefunction (note: |k| <= l). Alternatively to a number one may input the string slz or jbf. k will then be set to the value of the parameter slz or jbf, respectively. Note: if slz or jbf is used, it must be defined in the parameter section of the input-file or via an option on the command line. (operator file definitions come too late).
kmin Minimum value of the body fixed magnetic quantum number for the initial wavefunction.
Default is the corresponding value given in the primitive-basis-section.
kmax Maximum value of the body fixed magnetic quantum number for the initial wavefunction.
Default is the corresponding value given in the primitive-basis-section.
dk Step of the body fixed magnetic quantum number k.
Default is 1.
Wigner j Initial rotational quantum number.
symmetry Indicates symmetry.
nosym: all j values are used.
sym: either odd or even j values are is used, depending on whether j is odd or even (NOT YET IMPLEMENTED).
excite=s Excitation algorithm for the generation of (unoccupied) spfs.
s=mkj: Excitation of states is in order of preference: m (third DOF), then k (second DOF), then j.
s=kmj: Excitation of states is in order of preference: k (second DOF), then m (third DOF), then j.
print Optional. Prints (j,k,m) quantum numbers of the generated spfs to the 'log'-file. (This option may tell you the difference between the excite=s options.)
map label Modelabel of the DOF from where to take the initial orbitals. (See note below).
tdpoint Function Name of time-dependent function defined in HAMILTONIAN-SECTION_Function.
gwp type Distribution type: seqex, simex, shell or shell-noq0.
centre The initial wave packet centre position.
momentum The initial wave packet momentum.
width The initial wave packet width parameter. Can be real,imag.
(R|S [= R]) The distribution option.
Either just space GWPs apart equidistantly by R
or if S is one of:
  • random: random placement. Default parameter R=0.8. A larger R results in wider spacing.
  • random-pq: random placement and momentum. Default parameter R=0.9. A larger R results in wider spacing.
  • ovl: Fixed overlap. Default parameter R=0.5. A larger R results in closer spacing. Beware: If R gets close to 1, very narrow spacing results in almost singular overlap-matrices and thus numerical difficulties.
gwpwidth=gwpwidth Width of the GWP basis functions. Optional. Default is width. If gwpwidth=width, only the first GWP will be populated. Beware: This will go all wonky if shell-noq0 is chosen for type!!! No checking for this is done in the code.
scw = R Scaling of GWP widths. Optional. Default is R=0.
widthi=gwpwidth*(1.0+R)i.
ho-gwp

The initial single-particle functions are formed as follows:

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

modelabel type parameters state = s

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

modelabel type parameters pack = p

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:

modelabel type parameters pop = i

which populates the i th function.

Reading the initial wavefunction

Reading the initial wavefunction
Keyword Description
file = S The string S is a path of a diretory which contains a restart file. Up to 8 restart files may be read. Each file statement must be followed by SPF and/or A line(s).
SPF SPF i -> k : Write the SPFs of electronic state i of the restart file just read to the electronic state k of the system.
SPF i -> k1,k2,k3,... : Write the SPFs to states k1,k2,k3,..., i.e. those states will have identical sets of SPFs.
A A i -> k : Write the A-vector block of electronic state i of the restart file just being read to the electronic state k of the system.
init_state = I If an A-vector is not specified, a Hartree product is assumed and placed on the electronic state of number I. If init_state is not given, init_state=1 is assumed.
orthopsi The SPFs are transformed to natural orbitals and then orthonormalized. orthopsi is the default.
noorthopsi The SPFs remain untouched, they are not re-orthonormalized
realpsi Only the real part of the wavefunction is used and Schmidt-orthonormalised as in the orthopsi-case.
ignore The different primitive bases error will be ignored.

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

If the initial state should be a simple Hartree product one may drop the A line and use init_state to specify the initial electronic state.

Read-Inwf
  file = <restart-directory>
  SPF 1 -> 1,2
  init_state = 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 only for MCTDH wavefunctions, not for density operators or selected CI wavefunctions. It does, however, works for exact wavefunctions, if both, the wavefunction read and the final wavefunction of exact format. 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 CO2 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 direcory which contains a restart file (e.g. name-direcory 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 transformend 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 neccessary to copy a file to have the same information on two different file-paths.

INTEGRATOR-SECTION

To control how often the mean field matrices are updated, one of the following keywords may be chosen. Default is VMF.
Keyword Description
VMF Variable mean fields: The mean fields are calculated at each integration step.
CMF = R, R1 Constant mean fields: The mean fields are kept constant over a variable time interval.
Equivalent to CMF/var for propagation and relaxation but equivalent to CMF/varphi for improved relaxation.
R is the initial time interval (in fs).
R1 is an accuracy parameter for the time-step control.
CMF/fix = R The mean fields are kept constant over a fixed time interval, specified by R (in fs).
CMF/var = R, R1 The mean fields are kept constant over a variable time interval, determined by the errors of both the MCTDH coefficients and the single-particle functions.
var is default; CMF/var is equivalent to CMF, except for improved relaxation.
R is the initial time interval (in fs).
R1 is an accuracy parameter for the time-step control.
CMF/var2 = R, R1 The mean fields are kept constant over a variable time interval, determined by the errors of both the MCTDH coefficients and the single-particle functions Using a mid-point scheme. Coefficients and SPFs propagated by the same integrator.
var2 is default for ML-MCTDH;
R is the initial time interval (in fs).
R1 is an accuracy parameter for the time-step control.
CMF/varphi = R, R1 The mean fields are kept constant over a variable time interval, determined by the error of the single-particle functions only.
In case of improved relaxation varphi is default making CMF/varphi equivalent to CMF in this case.
R is the initial time interval (in fs).
R1 is an accuracy parameter for the time-step control.
CMF/vara = R, R1 The mean fields are kept constant over a variable time interval, determined by the error of the MCTDH coefficients (A-vector) only.
R is the initial time interval (in fs).
R1 is an accuracy parameter for the time-step control.
ML-CMF = S In a multi-layer calculation, the CMF scheme is used as specified by the following options (see note below).
S = unified
S = split
The following keywords define the integrator to be used.
Keyword Description
ABM/S = I, R, R1 Adams-Bashforth-Moulton predictor-corrector integrator used for S = all, spf, A.
I = order
R = accuracy
R1 = initial stepsize
BS/S = I, R, R1 Bulirsch-Stoer extrapolation integrator used for S = all, spf, A.
I = maximal order
R = accuracy
R1 = initial stepsize
RK5/S = R, R1 (,S) Runge-Kutta integrator of order 5, used for S = all, spf, A.
R = accuracy
R1 = initial stepsize (can be omitted or set to zero, then the integrator guesses a suitable value)
S=imp-ortho (improved orthogonality of SPFs, see below)
RK8/S = R, R1 (,S) Runge-Kutta integrator of order 8, used for S = all, spf, A.
R = accuracy
R1 = initial stepsize (can be omitted or set to zero, then the integrator guesses a suitable value)
S=imp-ortho (improved orthogonality of SPFs, see below)
SIL/S = I, R, S1 SIL integrator used for S = all, spf, A.
I = maximal order
R = accuracy
S1 = standard: The standard error estimate is used (default).
S1 = novel: The improved error criterion is taken.
See note below!
CSIL/S = I, R, S1 complex-SIL integrator used for S = all, spf, A.
Same as SIL/S (see above), but the use of the (complex) Lanczos-Arnoldi integrator is enforced.
(SIL tries to make the choice Lanczos/Lanczos-Arnoldi automatically).
DAV = I, R
rDAV = I, R
rrDAV = I, R
cDAV = I, R
Davidson "integrator" allowed for improved relaxation only.
The keywords DAV and DAV/A are equivalent. See note below.
I = maximal order, R = accuracy
There are three routines: DAV is for hermitian Hamiltonians,
rDAV and rrDav are for real Hamiltonians and real wavefunctions. Memory is saved as only real data is stored. rrDAV uses more real arithmetic and is hence faster, but not general.(See note below).
cDAV is for non-hermitian Hamiltonians (resonances).
LSODA = R Lsoda integrator from ODEpack. Experimental. Currently implemented for VMF only.
R = accuracy
ZVODE/S = [R] {,Pn} ZVODE integrator used for S = all, spf, A. Experimental.
R = accuracy
Pn = A comma seperated list of other options:
Initial stepsize if real (can be omitted or set to zero, then the integrator determines a suitable value),
or 'imp-ortho' (improved orthogonality of SPFs, see below),
or 'stiff' (default is the non-stiff version of the integrator),
or 'no_overshoot' (integrates exactly to desired time, not past; avoids problems with the energy calculation).
or 'minstp=FLOAT' (minimum integration step size, if smaller steps are necessary to meet the accuracy requirements, fail; default=0.0)
or 'maxstp=FLOAT' (maximum integration step size, no default)
or 'first=FLOAT' (initial integration step size, no default)
or 'maxit=INT' (maximum number of integration steps per t_out period, if more steps are required, fail; default=10000)
if maxit<0, both of the above failures will be caught, the GWP with the highest contribution to overall integration error is removed from the dynamics and the integration is restarted at the last (t_out) integration step

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, RK5/8, LSODA or ZVODE 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 (and some yet-to-be-determined combinations for ZVODE, likely all). Default is BS/spf + SIL/A.

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:

@n for mode number n
@n@m for modes number n and m; and so on
def for all modes not specified otherwise (default integrator/settings; can also be chosen with S=all or S=spf)

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

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.

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 following keywords may be used to further define the method of propagation.
Keyword Description
eps_inv = R R is the value used to regularise the inverse of the reduced density matrices. (Default: 10-8. See eq.(82) review.)
proj-h One dimensional Hamiltonians are not extracted to be in front of the projector. (See eq.(45) review)
h-proj One dimensional Hamiltonians are extracted to be in front of the projector. (See eq.(44) review)
energyorb Energy orbitals are propagated in place of spfs. Only for CMF. The spfs are propagated normally, but the whole wavefunction is transformed to energy orbital picture after each output (out1). energyorb is default when the Davidson "integrator" is employed (improved relaxation). energyorb requires that the orben file is set (Run-Section).
stdorb Standard orbitals are propagated. This is default, except for improved relaxation using DAV.
simple-proj The simple projector is used rather than the improved one. The inversion of the spf-overlap-matrix is thus avoided.
nohsym The symmetry of the operators determined by the program is not used to calculate the operator matirx elements.
CDVR Multi-dimensional potential terms are evaluated using CDVR . The analytic_pes keyword needs to be set in the OPERATOR-SECTION
TDDVR Multi-dimensional potential terms are evaluated using TDDVR. The analytic_pes keyword needs to be set in the OPERATOR-SECTION
direct-multi-d Multi-dimensional potential terms are evaluated using a direct algorithm, i.e. the potential is not stored on the full grid but re-calculated each time it is needed. This is much slower, but needs minimal memory. The analytic_pes keyword needs to be set in the OPERATOR-SECTION
update = R For numerically exact calculations, if the SIL integrator has been chosen the update time (i.e. step size) is set to R. The default is R = tout. Note that R should be an integer fraction of tout and should be chosen such that the SIL executes between 4 and 12 iterations (inspect the steps file).
lmf The equations of motion of the Linear mean-field approach are used to propagate density operators of type II (no effect for type I or wave functions).
df The equations of motion derived from the Dirac-Frenkel/McLachlan variational principle are used to propagate density operators of type II (no effect for type I or wave functions). This is the default.
imp-ortho Improve the orthonormality of the SPFs after each integration step. The long form improved-orthonormality is also possible. The keyword may also be given as a third argument to the RK5/8 integrators.
Note: currently works only with the RK5/8 integrators.
default Forces the integrator to use the accuracy parameters specified by the user rather than lossening them as occurs with GWP runs.
The following keywords may be used to further define the method of propagation for GWPs.
gwp_model = S1 Choose the model for the GWP propagation. S1 may be:
vmcg: GWPs follow full variational trajectories
classical: GWPs follow classical trajectories
independent: GWPs follow classical trajectories and wavepacket expansion coefficients are time-independent (independent GWP approximation)
static: GWPs are static (time-independent basis set) time-independent: as static
bel: GWPs are propagated using Basis Expansion Leaping algorithm (see below)
gwp_eoms = S1 [,S2] Choose details of how the GWP equations of motion are integrated. S1 or S2 may be:
inverse:
The eoms are solved by inversion of the C-matrix using its eigenvalues.
gauss or lu_decomp:
The eoms are solved by LU decomposition (same as Gauss elimination method).
conj_grad:
The eoms are solved using a conjugate gradient iterative method
direct:
The eoms are solved by direct inversion of the C-matrix
tikhonov:
Use Tikhonov regularization to solve eoms: ∂tΛ=(C2+diag(tikh_fac))-1 C Y, where tikh_fac is a small scalar factor that can be input after a comma (default value is 1×10-12). This is the Default option.
svd:
The eoms are solved using a singular value decomposition of the C-matrix
cholesky:
Solve the eoms via a Cholesky decomposition of the C-matrix, after equilibration/scalling of the system of equations
clsq:
Undocumented!
In addition, S1 or S2 may be:
ovr_eigen:
The inverse of the overlap matrix is formed using the eigenvalues (Default)
ovr_direct:
The inverse of the overlap matrix is formed directly
cx = S Controls the use of the Y = CX decomposition. If used, ∂tΛ = X + C-1 YR rather than ∂tΛ=C-1 Y. S can be:
dynamic :
decomposition used.
static :
decomposition used, but X ignored (propagate only using C-1 YR).
none :
no decomposition used.
no_cx The GWP EOMS are integrated using the full C-1 Y rather than making use of the the Y = CX decomposition. Same as cx = none
eps_gwp = R The eps value for inverting the C-matrix is set to R.
Default R = 1.0d-8.
lha Hamiltonian matrix elements are evaluated by operating on the GWP on the RHS and using an LHA of the PES.
lhasymm Hamiltonian matrix elements are evaluated by operating on the mid-point between 2 GWPs and using an LHA of the PES. This produces Hermitian matrices.
mcintegration Hamiltonian matrix elements are evaluated using Mont-Carlo integration.
lhastep = R The step used for numerical derivatives in the LHA is R. Default R = 0.01
gwp_coupling = S1 [,S2[,S3]]
or
dynamic_gwps = S1 [,S2[,S3]]
This keyword controls which GWPs are coupled together by including in the C matrix. This controls numerical problems due to the projector. By default, the GWP parameters are only included the diagonal value of the C-matrix is above the parameter defined by eps_dyn (see below). This behaviour can be changed by choosing S1, S2 or S3 to be:
all: all parameters are included all of the time.
none: no parameters are included all of the time.
min: parameters are made time-dependent using eps_gwp = 1.0d-4
normal: parameters are made time-dependent using eps_gwp = 1.0d-8 (default).
In addition, S1, S2 or S3 may be:
whole_g: The parameters of Gaussian functions will be included or excluded for all degrees of freedom, not just those that meet the eps_gwp criterion.
restart_int: When the identity of the included dynamic parameters changes the integrator is restart (i.e. the history is forgotten). Currently only implemented for ...
tind: parameters no coupled are kept stationary.

NB: This keyword is affected by the CX flag. If the CX decomposition is not used (no_cx), then uncoupled GWP parameters are stationary.
eps_dyn = R The value of epsdyn used to select the parameters included in the C-matrix. Default is 1.0d-8.
cysinv Use the inverse of the overlap matrix in the propagation of the SPFs and in the projector used to form the C-matrix and Y-vector, rather than forming an orthonormal projector.
smat_ld = S1 Specifies the method by which linear dependencies in the overlap matrix are determined. S1 may be either:
overlaps: Set GWPs to be linearly dependent if their overlaps exceed the value set by ovl_thresh (default is 0.95).
eigen: Calculates the eigenvalues and eigenvectors of the overlap matrix and sets GWPs to be dependent if they make a large contribution to the eigenvectors with eigenvalue less than eps_ld (1.0d-6 by default).
eps_ld = R Parameter to determine linear dependence when used with smat_ld=eigen (default is 1.0d-6).

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 that 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. Thus the functions should be time-independent and their EOMs undefined.

Whether a GWP parameter should be time-depedent can be determined by the relevant diagonal C-matrix element. If this is above a certain parameter 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.

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:

BEL = R Time propagation is only performed in the expansion coefficients with re-expansions of the wave packet into a new GWP set at time intervals R. BEL requires
  • VMF
  • no_cx
  • dynamic_gwps = none
and in the INIT_WF-SECTION
  • gwp_plain
If this value is negative, its absolute will be used as the interval and no re-expansion will be performed before the start of the computation.
BELflags = R[,R[,R...]] Flags for BEL. All are optional but order is significant. If a later one is to specified but earlier ones are to be left at default values, an appropriate number of commas without values has to be supplied. Default values are given in parenthesis. +/- after parenthesis indicates that higher accuracy results from larger/smaller values.
  1. maximum acceptable deviation of the reassignment (1e-2)-
    If this value is negative, it's absolute will be used and the interval between re-expansions will not be adjusted if the maximum number of GWPs is too low to keep the deviation low.
  2. minimum significant density value (1e-2)-
  3. LHA compliance limit (1e-2)-
  4. number of additional GWPs that are to be added after deviation limit is reached (40)+
  5. absolute minimum width of GWPs (1e-3)-
  6. absolute maximum width of GWPs (3e-0)+
  7. 1: recompute func in integrator only once
    0: every time
    (1)
    no influence on accuracy, but a lot faster if ==1
  8. 1: minimize deviation by optimization of GWP momentum parameter only
    2: optimize GWP momentum and width
    (2)
    more accurate if ==2 but minimum width for LHA critical potentials has to be chosen judiciously
  9. maximum number of simplex minimization steps (4000)+
  10. width of GWPs = width of wave packet * BELflags[10] (1)
    subject to subsequent optimization if BELflags[8]==2
  11. minimum overlap matrix eigenvalue accepted before switching off non-orthogonal projection (1e-5)
  12. scale factor for the wave function residual density contribution to the position selection cost function (1)
BELprint = I [: R1[,R1..] : R2[,R2..] : R3[,R3..]] Debug output for BEL. At every reassignment, the wavefunction and new and old values for GWP parameters are logged to BELout file. If the Rs are given, plot from R1 to R2 with step size R3 for every dimension. If no Rs are given, I is the total number of points distributed equidistantly in a N dimensional rectangular region spanned by the maximum extension of the GWPs and their widths. A negative sign for I may be specified to have only the reduced densities along each dimension plotted with |I| points.

Omitted integrator parameters default to:

Keyword Default
CMF/fix R = min(1.0,tout)
CMF/var R = min(1.0,tout) R1 = 1.0e-6
CMF/varphi R = min(1.0,tout) R1 = 1.0e-6
ABM I = 6 R = 1.0e-5 R1 = 1.0e-4
BS I = 8 R = 1.0e-6 R1 = update-time/2 for CMF or
R1 = 0.2 for VMF
SIL I = 30 R = 1.0e-6 S1 = standard
RK5 R = 1.0e-6 R1 = 0 (i.e. auto-guess)
RK8 R = 1.0e-6 R1 = 0 (i.e. auto-guess)
proj-h / h-proj in VMF mode: h-proj
in CMF mode: proj-h
for relaxation: proj-h
eps_inv R = 1.0e-8
eps_no R = 1.0e-8
eps_gwp R = 1.0e-8

Example Inputs for wave functions

Example Inputs for use of CPP with QUANTICS

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

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

Example Inputs for density operators of type I

Example Inputs for density operators of type II