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.
The order and placement of the keywords within a section is usually arbitrary. Exceptions are described below.
rd sin 36 3.800 5.60This input format is used in the PRIMITIVEBASISSECTION and, similarly, in the INIT_WFSECTION. In the operatorfile (*.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 OPERATORSECTION of the inputfile when 'alterlabels/endalterlabels' is used. Note that these latter inputs are not freeformat, the order of the arguments matters! Moreover, the fixedformat inputs must have a line for their own; they must not be followed by another keyword.
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 
cm1  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  10^{7}/cm1 
aJ  atto Joule (10^{18}J)  1.602177d1*ev 
invev  (electron volts)^1  0.036749325398 
debye  unit of electrical dipole moment  1 / 0.39343 
AMU  atomic mass unit  1 / 1822.88848325 
pmass  mass of proton  1 / 1836.15267247 
Hmass  mass of Hydrogen atom  1 / 1837.15264409 
Dmass  mass of Deuterium atom  1 / 3671.482934845 
Angst  Angstroem  0.52917720859 
pm  picometer  52.917720859 
nm  nanometer  0.052917720859 
deg  degree  180.0/pi 
Angst1  Angstroem^1  1/Angst 
pm1  picometer^1  1/pm 
nm1  nanometer^1  1/nm 
eV*AMU  sqrt(2mE)  
fs  femtosecond  1/41.34137333656 = 0.02418884326505 
ps  picosecond  1/41341.37333656 = 2.418884326505d5 
invfs  femtosecond^1  41.34137333656 
The input sections reflect the structure of the program. The
RUNSECTION 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 predefined order for the sections.
XXX  Description 

RUN  Whether propagation, relaxation, or diagonalisation, whether to generate or read DVR information, propagation and outputtimes, etc. 
PRIMITIVEBASIS  Definition of primitive basis, e.g. DVR / FFT. Note: PBASIS is a short form for PRIMITIVEBASIS . 
SPFBASIS  Definition of the singleparticle 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 SPFBASIS . 
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.
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 explicitly stated, will be written.  
overwrite  Any files already in namedirectory will be
overwritten. If this keyword is not given, and files already exist in namedirectory, 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 runsection, 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 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  A special operator file, called pes, will be generated, to be used by showsys to plot the PES. (See below). 
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 Gmatrix of the kinetic energy. (See below). Rather than setting the keyword gengmat one may alternatively give the option gmat when running mctdh. 
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=unmodified  4  Only for ML, like relaxation above. By default the ML Avectors of a relaxation starting as Hartree product will be modified by adding small random numbers. This mitigates the initial singularity of the densities and speedsup the relaxation. However, there are cases, where this modification is unwanted. The argument unmodified switches off this modification. 
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 Ith 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, 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 namedirectory will be performed. 
diagonalisation = I  4  The Hamiltonian will be diagonalised using a Lanczos algorithm with I iterations. The SPFBASISSECTION and the INTEGRATORSECTION are ignored for a (Lanczos) diagonalisation run. The WF is expanded in exact format. See remarks below. 
There are 4 levels of calculation types, reflecting the four
stages of a calculation. Only one runtype 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 Avector 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 INTEGRATORSECTION.
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
INTEGRATORSECTION for details. In short,
DAV is for (general) hermitian Hamiltonians, rDAV
and rrDAV are for realsymmetric Hamiltonians, and
cDAV is for complex nonhermitian 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 Avectors 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 Ith 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 Ith state
of the Davidson (or Lanczos) matrix will in general not be the
Ith 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 Avector 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_WFSection. 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 RunSection.
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 restarted
and the WF is read in via the file statement. Then the
Avector 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. (HE)**2 is diagonalised
within the space of Davidsonvectors 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 JacobiDavidson 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 Avector.
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 singleparticle 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
energyunit for the output to the rlx_info file. One also may
apply an energyshift (e.g. subtract the groundstate 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 blockimprovedrelaxation is performed when the packets keyword is given in the SPFBasisSection. The DAV, RDAV or RRDAV keyword must be given in the IntegratorSection and the keywords lock,quad,backrotate,cn are not implemented for blockimprovedrelaxation. For an example input on blockimprovedrelaxation see blkHONO.inp.
Notes on Diagonalisation.
With the keyword diagonalisation the program will perform
a diagonalisation of the Hamiltonian using a straightforward
(non sophisticated) Lanczos algorithm. For this, the wave function
is in exact format. Although the Lanczos algorithm is part
of the MCTDH package, it is not based on the MCTDH method and
hence does not inherit its efficiency. Moreover, diagonalisation
is not parallelized. Diagonalisation has been
implemented to check filterdiagonalisation. We do not recommend
the use of diagonalisation, in particular for larger
systems. It will not work for wave functions with more than
10^{9} grid points.
Keywords augmenting the calculation type  

exact  A numerically exact wavepacket calculation will be made. 
The exact keyword can be added to any runtype
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.
Keywords modifying the basis set onthefly  

spawn = R1(,R2,R3,R4)  The number of SPFs will be increased onthefly according to
the algorithm described in
J. Chem. Phys. 153, 234114 (2020).
Compulsory parameters:

initorb = (S1,R1(,S2,R2))  This keyword controls the generation of SPFs at time t = 0 and
it is automatically invoked by the keyword "spawn" described
above.
All parameters are optional, thus it should be used when one
requires different settings from the recomended defaults.
Optional parameters:

dynspf = (S1,I,R1(S2,R2))  This keyword controls the generation of SPFs at time t > 0
and it is automatically invoked by the keyword "spawn"
described above.
All parameters are optional, thus it should be used when one
requires different settings from the recomended defaults.
Optional parameters:

Keywords to keep particle number of each SPF constant  

numberproj = S  This is applicable for the second quantization calculation
(MCTDHSQR) where each electronic/bosonic DOF consist
of several Hilbert spaces and you want to keep the particle
number of each SPF constant during propagation or relaxation.
This works only, if the initial wavefunction has sharp
particle numbers for both occupied and unoccupied SPFs.
These particle numbers are then forced to be conserved
during propagation/relaxation.
The file should be in ascii format and must contain the following data: nh(1) nh(2) ... nh(d) [Number of Hilbert space in each DOF separated by space] g(1,1) g(2,1) ... g(nh(1),1) [1st grid point of each Hilbert space of the first DOF separated by space] g(1,2) g(2,2) ... g(nh(2),2) [1st grid point of each Hilbert space of the second DOF separated by space] ... ... g(1,d) g(2,d) ... g(nh(d),d) [1st grid point of each Hilbert space of the dth DOF separated by space] Here, nh(i) is the total number of different Hilbert space in the ith DOF and g(i,j)) is the index of the 1st grid point of the ith Hilbert space in the jth DOF. This keyword is subject to certain restrictions: (1) Both the occupied and unoccupied SPFs must have integer particle numbers initially. (2) Mode combination (by MCTDH) cannot be used when using the numberproj keyword. However, one may modecombine several physical number states and take them as one SQRDOF in MCTDH. 
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 Gmatrix 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. Rather than setting the keyword gengmat one may alternatively give the option gmat when running mctdh. 
These keywords are equivalent to a runtype keyword of the level given.
Keywords defining how the readwrite files will be handled.  

readdvr = S  The DVR information will be taken from the DVR file in directory S. Primitivebasissection is ignored.  
readoper=S  The operator information will be taken from the OPER file
in directory S. Operatorsection is ignored. If the operator is read, it is a good habit to read the DVR file as well. 

readalloc=S  An alloc file is read from directory S. If S is not given,
S=name is assumed. If the DVR and/or operator is read, one may read the alloc file as well. 

readinwf=S  The initial wavefunction will be taken from the RESTART
file in directory S. InitWFsection is ignored. Note: additionally the SPFbasissection will be ignored. This information is read from the restart file. If readdvr, readoper, and readinwf are set, the input file consist only of run and integratorsection. 

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 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 CMFstep. This is useful for improved relaxation. If tpsi is set in addition to tout=all, then the output interval is limited to the value given by tpsi. 
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. However, if tout=all is set, then tpsi is used to limit the output interval. 
tstop =S  The stoptime (realtime) 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 realtime. 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 stoptime (cputime) 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 cputime. NB: If stopfile is removed, the run will stop after the next output. 
twall = S  Similar to tcpu, but the elapsed walltime (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 mctdh86 h). 
thermal = R, I  With this keyword one sets tfinal=1/(2*k*T), where
k denotes the Boltzmann constant and T the temperature,
which is given (in Kelvin) by the first argument R.
The second argument I is the seed for the random number generator
used to create a random initial wavefunction. When used together
with the relaxation keyword one can converge an initial state
for a time propagation and the evaluation of the system observables.
Then averaging of many random phase sets leads to the thermal
observables. Note that the keyword tfinal must not be given. If tout and/or tpsi are given, their values will be adjusted to the computed tfinal. This keyword also requires to specify which degrees of freedom are to be thermalized, which must be done though the ran keyword in the INIT_WFSECTION. The electronic degree of freedom cannot be thermalized. See Building the initial wavepacket for more details, Chem. Phys. Lett. 349, 321328 (2001) for a description of the thermalization algorithm implemented and Chem. Phys. 482, 113123 (2017) for its convergence with the number of seeds. Works for both MCTDH and MLMCTDH. AS MCTDH always normalizes the wavefunction, the loss in norm due to relaxation is stored in the variable normthermal. normthermal is written to the log file (near its end), as well as to the restart file. From the latter it may be read with the aid of the analyze routine rdrestart86. See files inputs/thermal4D_s1.inp and inputs/thermal4D_ML_s1.inp for two respective examples. Please see also the MCTDH User's Guide chapter 2.8 and chapter 15. 
usepthreads = I (,S(,S1(,S2,..)))  The program runs in a shared memory parallelised modus,
using I pthreads. The parallelisation of the different parallel
subroutines may be switched off (in order to save memory and
sometimes CPU time) by
setting the keywords nosummf, nofunka,
nomfields, nophihphi, nofunkphi,
nohlochphi, nogetdavmat,
nodsyev, getdiag. If the keywords memcalcha or memmfields are set then more memory is used for a more efficient parallelisation. The keyword getdiag switches on parallel evaluation of the diagonal Hamiltonian matrix elements (only used within Davidson diagonalisation). It is not default, because it requires additional memory. Note: getdiag requires about the same amount of extra memory as memmfields, i.e., it usually makes sense to use both keywords together. Using the nosummf2 (or alternatively summf1) keyword, MCTDH uses a different kind of parallelisation for the summf routine. With the dsyev = I keyword the minimum matrix size for the parallel dsyev routine can be determined. 
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 nosummf, nocalcha,
nofunka2, nophihphi, nohlochphi,
nomfields, nogetdavmat,
nodsyev or nodav, nogetdiag. Other than for usepthreads parallel evaluation of the diagonal Hamiltonian matrix elements is the default (as it does not require additional memory) and may only be switched off with nogetdiag. With the dav = I keyword the number of Davidson vectors per MPI process can be adjusted in order to reduce communication cost. If 'expect = ...' is present in the RUNSECTION one may also set noexpect in which case the expectation value is not evaluated in parallel (useful for small operators). 
energynoteV  The eVconversion factor is set to 1. This is for running models in dimensionless coordinates. 
timenotfs  The fsconversion factor is set to 1. This is for running models in dimensionless coordinates. 
normstop=R  The program will be stopped if norm < R. 
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 criterion is reached for all modes. If I is the string "any", stop when this criterion 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.0d5,eV 
precon=I  An improved relaxation run with a Davidson integrator
(DAV, rDAV, rrDAV, or cDAV) may use a better preconditioner
than just the diagonal. If I.gt.1 a IxI dimensional block of
the Hamiltonian matrix is build, inverted and used as
preconditioner. If a preconditioner is used, the
Olsencorrection should be enabled. A preconditioner should not be used when symcoeff is set, as its use might destroy symmetry. 
splitrst  If a blockDavidson improved relaxation run
is performed, the use of this keyword will finally split
the multipacket restart file into a series of singlepacket
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 mctdh86 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 namedirectory. splitrst works also for geninwf runs. I.e. replace relaxation with geninwf and use the file statement to read a restart file of a blockimprovedrelaxation. (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 energyunit S to the rlx_info file. Additionally an energyshift 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=energyvalue, S=energyunit), 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 runtype 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. 
rlximin=R(,S)  Similar to rlxemin, only for cDAV. With this keyword one can exclude in a cDAVlock run all states with an imaginary part of the energy lower than rlximin. rlxemin and rlxemax refer to the real part of the cDAV complex energy. 
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 Avector accordingly. This saves both, CPUtime 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 importho in the IntegratorSection is recommended. The sizes of the removed imaginary parts is reported in the logfile. If FFT or PLeg are not used, there is no point in setting realphi. 
allcomplete  This keyword inhibits the propagation of the SPFs (it sets the logical complete to .true. for all modes). This estabishes an alternative way to perform 'exact' or 'standard method' calculatios, as the initial SPFs are not altered. Note that allcomplete requires CMF/fix as integration scheme. allcomplete must not be used for MLruns. 
optcntrl(=S(,S1)) 
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 onthefly calculated new field requires the the additional keayword update, i.e., optcntrl=simprop, update. Note: simprop and pc cannot be used together. 
The tfinal keyword must be given. All other keywords are optional. The following default values are used.
Keyword  Default 

tinit  0.0 
tout  tfinaltinit 
tpsi  tout 
The following keywords define the data calculated and saved. If keywords are omitted, data will not be calculated.
Keywords defining outputdata files to be opened.  

Keyword  Description 
all  All the (optional) files discussed below will be opened. More precisely, a propagation/relaxation run will open the files: auto, autoe, gridpop, output, pdensity, psi, speed, steps, stop, timing, ptiming and update. In a diagonalisation run, only the files output, timing, lanczvec and eigvec are opened. 
auto = S (,S1) (,S2)  The autocorrelation function will be written to the file
auto. If S = twice, autocorrelation function is written twice in interval tout (only for CMF) . If S = once, autocorrelation function is written only once in interval tout. If S = no, no autofile will be opened. This only makes sense, if the keyword all was given previously. If S, S1, or S2 = error, the autoe file will be opened, which contains the autocorrelation function computed with the least important natural orbital omitted. This information is useful for estimating the error. 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 autocorrelation function, respectively, needed in the filterdiagonalisation method. Note that in a multipacket run, i.e. when packets > 1, the files auto, auto1, and auto2 contain cross rather than autocorrelation functions. WARNING: The so called t/2 trick is used, which assumes a real initial state and a symmetric Hamiltonian. If these conditions are not met, use cross instead. 
auto  Synonymous for auto = once. 
cross (= S (,I))  A crosscorrelation
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 autocorrelation function will be calculated (but without the t/2trick). This is also the default (i.e. if no argument is given). If an integer I is given after the path S, then the crosscorrelation 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 crossdata. 
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)Spsi(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 SystemHamiltonian is derived, i.e. the
total systemenergy.) The norm (not norm**2) of Psi is
additionally written to the expectation file. 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 realonly , 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 multiset run will be stored separately. 
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 meanfield operators, are calculated and written to the orben file. Note: orben must be set, when the propagation is in energy orbitals (keyword energyorb, IntegratorSection) 
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 nooutput. 
pdensity (=I1,I2,I3,I4)  The oneparticle density will be written to the file pdensity. If the pdensity keyword is followed by an equal sign and up to four integers, the oneparticle 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 or S1 = no, no psifile will be opened. This only makes sense, if the keyword all was given previously. 
psi  Synonymous for psi = single. 
speed  The CPUtime used within an output interval will be written to the file speed. (default). 
nospeed  The speed file is not opened. 
steps  Information on the integrator step sizes will be written to the file steps. For multilayer runs with mlcmf=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 (realtime and/or CPUtime) to the stop file. (default). 
nostop  The stop file is not opened. 
timing  Program timing information will be written to the file timing. (default). 
ptiming (=all)  For a serial MCTDHrun all parallelized routines are included in the timing file. If pthreads are used then the ptiming file is created. This file includes timing information for each thread. The ptiming file is not read in a continuation run but recreated. If the MPI parallel version of MCTDH is used then additionally the mpitiming file is created. This file contains timing information for the different MPI prozesses. If the option =all is set for each MPI prozess a ptiming file will be created. 
notiming  The timing file is not opened. 
noptiming  The ptiming file is not opened. 
update  If the constantmeanfield integrator with adaptive step size is used, the update time for the meanfields is written to the file update. (default). 
noupdate  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. 
graphviz (Can be abbreviated to graph) 
For MultiLayer runs, generate an input file for
graphviz
by which one can visualize the ML tree. The generated file
will be called "mltree.dot". If the graphviz software is
installed (try "dot help"), one can e.g. use the command dot Tx11 mltree.dot to display the tree in an X window on the screen. Circles represent the multilayer 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 singleparticle functions of the submodes, or the numbers of primitive basis functions, respectively. To save the plot to a file, here to tree.pdf, use the command dot Tpdf mltree.dot > mltree.pdf NOTE: For version 8.5.11.1 and higher, the mctdh program automatically creates the file mltree.pdf, when the keyword graphviz is set. 
Keyword  Default 

auto  S = once 
cross  S = name 
psi  S = single 
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. 
alterparameters ..... endalterparameters 
The lines between the keywords define parameters to be used in building the Hamiltonian, using the same format as in the PARAMETERSECTION 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 PARAMETERSECTION of the .op file.
The file must end with the line endparameterfile 
alterlabels ..... endalterlabels 
The lines between the keywords redefine labels specified in the LABELSSECTION of the operator file opname.op 
v < R  Energy cutoff used for potential energy surface in exact calculations. All potential energy values greater than R are set to R. Note, this keyword is ignored in MCTDH calculations. 
v > R  Energy cutoff used for potential energy surface in exact calculations. All potential energy values less than R are set to R. Note, this keyword is ignored in MCTDH calculations. 
analytic_pes  If the operator contains a nonseparable potential this will not be stored explicitly on the primitive grid points, but in an analytic form which can be used to generate the potential onthefly 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 nondiagonal Hamiltonian terms which are on all gridpoints smaller than 1.d12 au are removed as well). The default is cutoff=tiny (i.e. 1.d9 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 protocolled 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 LabelsSection 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.0d9 au). I.e., even when the natpotcut keyword is not given, all terms which are smaller then 1.0d9 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.d6) 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. 
fast = V1,V2{n},...  An more efficient algorithm for H(natpot)*A operation is used. The Avector is premultiplied 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 LabelsSection 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 premultiplication. The larger n is, the larger will be the speedup. However, n is limited to min(4,nmode2), 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,nmode2) 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 multiset propagation. 
reducepf = I, R_1, R_2, ...  I= number of natpot, R_kappa = reduceweight for mode kappa: wpf(kappa). A (large) natpot may be reduced by ignoring all terms for which sum_kappa j_kappa * wpf(kappa) > 1.0, where the contracted mode is ignored in this sum. Note that the mode numbers are the potfit ones which may differ from the MCTDH ones. The weight of the contracted mode can be given any value, it will be set to zero automatically. 
printnpot  Full information on natpots is printed to op.log file.
If printnpot is not given, all npot lines except the
first and the last one of each natpot are suppressed. This keyword is also used enable to the printing of all kterms to the op.log file. (Only the first 40 terms of each operator are printed by default. 
If no OPERATORSECTION 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.
For a MLMCTDH calculation the SPFBASISSECTION must be replaced by a MLBASISSECTION which defines the MLtree. The top layer (upmost Avector) 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 5indicates 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 modelabel q1. If mode combination is used, one would e.g. write
2> [q1,q2]
mlbasissection 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] endmlbasissectionFor a graphical representation of the tree click here .
MLBasisSection 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] endmlbasissectionFor a graphical representation of the tree click here .
The following lines define the singleparticle function basis to be used in a wave function calculation. 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 multiset 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 singleparticle 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 3mode system, with labels X, Y and Z,
to define an spf basis of 3 functions per mode,
spfbasissection X = 3 Y = 3 Z = 3 endspfbasissectionor
spfbasissection X = 3 Y = 3 Z = 3 endspfbasisTo contract the degrees of freedom X and Y into a single mode,
spfbasis X, Y = 3 Z = 3 endspfbasissectionIf a multiset basis is used, then to have 3 functions in the first set and 2 in the second,
spfbasissection multiset X , Y = 3 , 2 Z = 3,2 endspfbasissection
Note: The electronic SPFBasis is not to be specified in the spfbasissection, as it is always complete. The electronic mode will always be the last mode in a singleset run.
If many electronic states are present, the definition of singleparticle functions for a mode can be continued on a second line by using a continuation mark %, e.g.
spfbasissection X = 5 , 5 , 5 , 5 , 6 , 10 , 5 , 2 , 2 , % 5 , 5 , 2 , 2 , 3 endspfbasissection
If for symmetry reasons one set of SPFs is always identical to another one, the latter set need not to be propagated numerically. In such a situation, e.g. H_{2}O in valence coordinates and for a symmetric initial state, one may tell the program not to propagate the second identical set of SPFs. This is done by specifying with the id keyword the mode to which the present mode is identical. E.g.:
spfbasissection R1 = 8 R2 = id,1 Theta = 9 endspfbasissectionMode 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 multistate or mutipacket runs. Note: A SPFBASISSECTION does not need to exist for an exact calculation, except if packets is specified.
Keyword  Description  Default 

multiset  If an electronic basis is defined it is treated using the multiset formalism, i.e. a set of singleparticle functions per state. If this keyword is not present, the singleset formalism is used.  not set 
singleset  If an electronic basis is defined it is treated using the singleset formalism, i.e. a common set of singleparticle functions for all electronic states. This keyword need not to be given, singleset is default for electronic states.  not set 
packets = I (,S)  I independent wavepackets will be simultaneously propagated. Technically, the packets are propagated on auxiliary electronic states. One may choose between S=singleset or S=multiset. (multiset is default here). For block improved relaxation, singleset must be given. A block improved relaxation requires the packet keyword and a DAV, RDAV or RRDAV "integrator".  I = 1, multiset 
id,I  I denotes the modenumber (particlenumber) with which the present mode is to be identified. See the note above for the correct use of the id keyword.  not set 
noredundancycheck  If this keyword is given, the check on redundant SPFs is suppressed. Use with care.  not set 
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. 
HO  Harmonic oscillator (Hermite) DVR. 
rHO  Radial Harmonic oscillator (oddHermite) DVR. 
Leg  Rotator (Legendre) DVR. 
Leg/R  Rotator (Legendre) DVR for a restricted range on angles. 
Lagu1  Laguerre DVR for boundary condition x^{1/2}. 
Lagu2  Laguerre DVR for boundary condition x^{1}. 
Lagu3  Laguerre DVR for boundary condition x^{3/2}. 
Lagu4  Laguerre DVR for boundary condition x^{2}. 
sin  Sine (Chebyshev) DVR. 
FFT  Fast Fourier transform collocation. 
exp  Exponential DVR. Periodic boundary condition. 
cos  Cos DVR. "gerade" solutions with periodic boundary condition. 
sphFBR  Spherical harmonics FBR (deprecated). 
KLeg  Extended Legendre DVR. 
K  Kquantum number appearing with KLegDVR. 
PLeg  TwoDimensional Legendre DVR. 
Wigner  Wignerd (threedimensional rotor) DVR. 
Extern  External DVR. 
basis_size is an integer specifying the primitive basis
size, e.g. grid points or vector elements etc. Note that for an
FFTrepresentation basis_size (i.e. gdim in the
program) must have a prime factor decomposition with only 2's,
3's and 5's but should have a decomposition with only 2's and 3's
for optimal performance, i. e. basis_size =
2^{m}3^{n} where m
and n are positive integers. 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
expDVR.
For a sphFBRrepresentation, 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 KDVR basis_size is not required in input. It will
be calculated from the kmin,kmax parameters given in this
section. NB: The basis_size is called gdim.
The parameters to be input depend on the basis type as follows:
Basis  Parameters  Parameter description 

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 = xixf. 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 halfaxis is used  is the lefthand 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 = xixf or S = x0xf. 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 inputfile or via an option on the command line. (alterparameter 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+N1 ). 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:  Starting point of the radial interval (usually zero). 
b:  Length parameter. Chi_{n}(x) = 1/b * Sqrt((xx_{0})/n) * exp((xx_{0})/(2*b)) * L^{1}_{n1} ((xx_{0})/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 = xixf. 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 = x0xf. This string serves as a switch between the three possible input formats. 
x0:  Starting 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 speriodic (linear is the default).  
fft or exp  string:  2pi , s2pi , c2pi or
2pi/m , s2pi/m , c2pi/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 expDVR which follows PLeg, the input k= kmin,kmax may follow. The default is kmax=kmin=(N1)/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 omitted (k=kmin), the range [kmin,kmin] is chosen. Default is the full range [(gdim1)/2, (gdim1)/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 (See Remarks on sphFBR) 
jmax:  Maximum value of quantum number j of the spherical harmonics basis functions. Note: jmax replaces N, the number of basisfunctions/gridpoints. 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. jmaxj = even). 

thrshld:  Threshold for convergence when the FBR integrals are performed. This input is optional, default: thrshld=1.d10.  
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 lvalues). 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 Kvalue. Default dk=1.  
PLeg  string:  Controls whether symmetry to be used. all: no symmetry (use all lvalues). odd: odd symmetry (i.e. l = odd). even: even symmetry (i.e. l = even). 
Wigner (See Remarks on Wigner DVR 
string:  Controls whether symmetry to be used. all: no symmetry (use all jvalues, default). 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. 
Remarks on harmonic oscillator DVR:
The HODVR 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 gridpoint. 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,cm1 1.0,AMU Y HO 32 xixf 0.528 0.528
Remarks on Laguerre DVR:
The LaguaDVR is build from the basis functions:
phi(n,x) = Sqrt((n1)!/(n+a1)!) * x^(a/2) *
exp(x/2) * L(n1,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((xx0)/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 subtraction. 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 xixf (not recommended in general) or
x0xf. 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 longFor the 2pi and sdq keywords see remarks on cosine DVR.
s1 sin 2 0.5 0.5 spinIf spin is set, 2x2 derivative matrices with zero diagonal and (1/2,1/2) (first derivative) or (1/2,1/2) (second derivative) offdiagonal elements are generated. I.e. (i/2)*sigma_y and (1/2)*sigma_x replace the derivative matrices, where sigma denotes the Pauli matrices. In other words, the Pauli matrices sigma_x, sigma_y, and sigma_z are given by the operators 2*dq^2, 2*I*dq, and 2*q, respectively. (Of course, factors should be moved to the coefficient column of a HamiltonianSection, but it is preferable to let the imaginary unit "I" stay with the operator).
Remarks on cosine DVR:
The basis functions underlying the DVR are 1/sqrt(L),
(2/sqrt(L))*cos[(j*pi/L)*(xx_{0})]. The symmetrized
derivative, sdq = 0.5*( sin[(pi/L)*(xx_{0})] * d/dx +
d/dx * sin[(pi/L)*(xx_{0})] ) is used as first
derivative. Because only gerade (symmetric) wavefunctions are
computed, we consider only the interval
[x_{0},x_{0}+L], although the wavefunction is
periodic on the interval [x_{0}L,x_{0}+L]. xi
and xf are the first and last gridpoint, but when long is
given, the input is interpreted as x_{0} and
x_{0}+L.
The use of the operators dq, qdq, and cdq is not allowed
for cosDVR.
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 expDVR have an identical set of input
parameters (if expDVR is not combined with
PLegDVR). These two methods are in fact largely
equivalent and produce identical results (for same parameters).
The numeric, however, is different and the expDVR will be
faster for small grids whereas FFT is faster for long
grids. Around 30 grid points both methods are of similar speed.
The use of the interactionpicture is possible for
expDVR.
FFT and expDVR 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 gridspacing is dx = (xfxi)/(N1). 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 gridspacing is dx = (xfxi)/N. The routine eingabe.f rescales xf > xfdx.
speriodic: xi and xf are considered as identical due to the periodic boundary conditions. The gridspacing is dx = (xfxi)/N. The grid points, however, are now placed symmetrically on the interval (xi,xf) and eingabe.f rescales xi > xi+dx/2, xf > xfdx/2.
2pi/mult or s2pi/mult or
c2pi/mult: The routine eingabe.f sets xi=0 and
xf=2*pi/mult and then performs according to the periodic
or speriodic keyword. For c2pi/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 speriodic x fft 32 s2pi or x fft 32 3.0434179 3.0434179 linear x fft 32 3.1415927 3.1415927 speriodic x fft 32 c2pi
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 elkeyword we assume that there are three electronic states.
primitivebasissection el el 3 X FFT 32 2.00 2.00 linear Y HO 32 0.00 5.2,eV 1.00 endprimitivebasissection
Remarks on External DVR:
ExternDVR is an external DVR, where grid points and
DVR derivative matrices are read from a file. The file can be
read in ascii or binary format:
in ascii format (default):
do i=1,gdim read(unit,*) ort(g) enddo do i=1,gdim read(unit,*) (dif2mat(j,i),j=1,gdim) enddo do i=1,gdim read(unit,*) (dif1mat(j,i),j=1,gdim) enddo
do i=1,gdim read(unit) ort(g) enddo do i=1,gdim read(unit) (dif2mat(j,i),j=1,gdim) enddo do i=1,gdim read(unit) (dif1mat(j,i),j=1,gdim) enddo
where gdim is the basis_size. The file can have absolute or relative path. If file contains only grid points (for example in POTFIT program, when only grid points are used), the DVR matrices will be zeroed. If only second derivative matrix dif2mat is given, dif1mat will be zeroed.
Example:
x extern 30 x_data y extern 50 y_data binary z extern 42 z_data binary angst theta extern 23 theta_data deg
Remarks on KLeg, PLeg and sphFBR:
KLeg, PLeg and sphFBR all define 2D singleparticle 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 Fouriertransform the potential (with projection86) 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), 619632. See also H2H2.inp and H2H2.op on the inputs or operatos directory, respectively
Example:
PRIMITIVEBASISSECTION 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 endprimitivebasissectionThe 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 Krange is mainly for tests.
The KLeg/K, PLeg/exp and sphFBR/phiFB combinations generate modeoperators. Hence the DOFs (theta2, K_th), (theta1,phi), and (alpha,beta) must be combined with each other in the SPFBasisSection.
Remarks on Wigner DVR:
Wigner DVR defines 3D singleparticle 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 SPFBasisSection.
Important Note: The order of the three WignerDOFs must be beta, gamma, alpha or beta, k, m when k denotes the BF and m the SFcomponent of the angular momentum. This ordering must hold in the PRIMITIVEBASISSECTION, in the combination scheme defined in the SPFBASISSECTION, and in the modes line of the HAMILTONIANSECTION. Note that the angles (beta,gamma,alpha) correspondent to the angular momenta (j,k,m), respectively.
Example:
PRIMITIVEBASISSECTION ........ beta wigner 20 all gamma exp 15 2pi alpha k 7 7 endprimitivebasissection SPFBASISSECTION ........ beta,gamma,alpha = 20 endspfbasissection
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 namedirectory. S2 = orthopsi : The wavefunction is transformed to natural orbital representation and the singleparticle functions are Schmidtorthogonalised (this is the default). S2 = noorthopsi : The initial wavefunction is not transformed and the SPFs are not Schmidtorthogonalised after being read. S2 = realpsi : Only the real part of the wavefunction will be taken, the imaginary part is ignored. The SPFs are Schmidtorthogonalised as in the orthopsicase. Note: Not implemented for multilayer runs. 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 sinDVR, 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. For multilayer runs, the read wavefunction must either be a multilayer wavefunction with the same tree structure (though the number of SPFs may differ), or it must be a standard MCTDH wavefunction (using the same primitive modes in the same order and with identical modecombinations). In the latter case, the SPFs for the bottom layer are taken directly from the MCTDH wavefunction, while the SPFs for upper layers are determined via a truncated Hierarchical SVD. The estimated truncation error will be printed to the log file. The numbers of bottom layer SPFs of the MLwavefunction may differ from mumbers of SPFs of the MCTDHwavefunction. 
select (S1 S2) I1 I2 ....  With the help of the keyword select one may read
a selection of the SPFs from a restart file given by the
file keyword. In input select must come after
file and each select statement must be a line of
its own. The first symbol S1 reads sk, where k is
an integer giving the electronic state (of a multiset run).
Similarly, the second symbol S2 reads mk,
wherek is an integer specifying the mode, as defined in
the SPFBASISSECTION. The following integers give the
numbers of the SPFs to be included. There must be exactly as
many numbers as there are SPFs for that mode (and state).
If the state symbol sk is not given, the first state,
s1, is assumed. (Similarly for mk).
Example equivalent input file=geninwfrun file=geninwfrun select s1 m1 1 3 4 select 1 3 4 select s1 m2 1 4 6 3 select m2 1 4 6 3 select s2 m2 1 3 5 select s2 m2 1 3 5 
build ..... endbuild 
The initial wavefunction will be build using the data specified between these keywords. See Building the initial wavepacket. 
readinwf ..... endreadinwf 
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. The wavefunctions must be in multiset format. This is a more powerful tool than the one invoked through the file keyword discussed above. See: Reading the initial wavepacket. 
blockSPF = 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
multipacket or a singlepacket 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 multipackets singleset
runs (e.g. blockimproved relaxation) only.
For other runs one should use readinwf. Note that the Avector is not yet defined, it must be defined via Acoeff or autoblock or blockA. S2 = noorthopsi : The initial wavefunction is not Schmidtorthogonalised after being read. S2 = orthopsi : The initial wavefunction is Schmidtorthogonalised after being read (default). S2 = realpsi : The only the real part of the SPFs will taken and the thus modified SPFs are Schmidtorthogonalised as in the orthopsicase 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 Avectors in case they are read in by a following blockA command. In case blockSPF is used in combination with blockA one MUST set noorthopsi because before orthonormalization the SPFs are transformed to natural orbitals and hence are no longer consistent with the Avectors read in. See: Generating block initial wavepackets. 
blockA = S1 (,S2,S3,..)  The Avectors of a multipackets singleset
initial wavefunction will be read from restart files
S1, S2, S3, .... The restart files must be of standard
nonpacket form. The number of restart files must match
exactly the block size. One may give more than one
blockA keyword. The SPFs may be generated by build, or  more likely  may be read from another restart file with the aid of the blockSPF keyword. In case blockSPF is used in combination with blockA one MUST set noorthopsi in blockSPF. If realpsi is set in blockSPF it will affect the Avectors read by blockA as well. See: Generating block initial wavepackets. 
realpsi  Only the real part of the WF will taken and the
thus modified SPFs are Schmidtorthogonalized and the
Avector is renormalized. 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,
readinwf, or blockSPF 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, readinwf, or blockSPF / blockA /
autoblock this WF
may be modified by one or several of the following procedures.
The program will take actions in the order:
Acoeff, meigenf,
operate, orthogonalise, correction,
realpsi,
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 necessary to copy a file such that the same information can be read via two different filepaths. Such a problem may arise when using blockSPF and blockA, or when using orthogonalise.
The information needed to generate the initial wavefunction is written one line per degree of freedom between the keywords build and endbuild. 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 primitivebasis 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.
Rather than init_state one may use euclid (see below)
to specify the initial state. This becomes important, if there are
several (singleset) electronic states.
To summarize:
Keyword  Description 

Acoeff ..... endAcoeff 
The lines between the keywords contain in free format one
integer (Aindex) and one complex number (Avalue) per line,
thus defining the initial Avector. The default is A(1) =
(1.0,0.0) and zero for all other coefficients. For a multiset run, the input line must contain two integers. The first is the Aindex and the second the electronic state. E.g.: 518 (0.717,0.0) (singleset) 518 1 (0.717,0.0) (multiset) There is also a long form of the Acoeff input. Here the number of the particle for each mode (and the electronic state(s) in case of a multiset run) is input, followed by the value of the Acoefficient. E.g.: 1 4 1 3 1 2 (0.5,0.0) would be a correct input for either a 6mode singleset, or a 5mode multiset WFcalculation. In the latter case the last integer, "2", would specify the electronic state. This redefinition of the Avector is also possible, when the initial wavefunction is read in (file keyword, restart run). The Avector is renormalised before being used. The nonzero entries of the Avector are protocolled in the logfile. 
autoblock  autoblock automatically generates a set of Avectors 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. 
meigenf = I,S,I1,S1 (,S2,S3,S4,I2,S5(=I3))  Generate the modeeigenfunctions using operator S, where
S = XX is the name of an operator, defined in the
HAMILTONIANSECTION_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 protocolled 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. If the initial wavepacket
is thermalized (see thermal keyword) the optional string
S5=ran needs to be given if this mode should be randomized.
By default the randomization will be averaged up to the
{basis_size1}th state, unless an alternative optional value
I3 is declared. Examples: meigenf = 3,oper,0 meigenf = 3,oper,follow,full,125 meigenf = 3,oper,follow,full,select,write,125 meigenf = 3,oper,0,full,noselect,write,ran=30 
operate = S(,S1(,S2,...))  Operate on the initial wavefunction using operator S,
where S = XX is the name of an operator, defined in the
HAMILTONIANSECTION_XX. One may apply up to 15 operators to
the WF. operate = O1, O2, O3 will produce
O3*O2*O1*psi>. The wavefunction is renormalised after
each application of an operator. The normalisation factors
are protocolled in the logfile. Note: MLwavefunctions and operators containing multilayer operators (such as mlpf) are currently not supported for this operation. 
operate_iter = I  Maximum no. of iterations to be used in applying an operator to an MCTDH wavefunction. Default = 10. Only the Avector 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.0d8 
operate_nonorm  Depreciated! The normalisation factor is removed from operator*psi, i.e. the initial wavefunction is no longer normalised. 
operate_ortho  Operate on the initial wavefunction using the operator (Dshift), where D is the operator defined through the operate keyword and where the number shift is automatically adjusted such, that the operated wavefunction is orthogonal to the initial one. For the usage of operate_ortho it is requested, that in the HamiltonianSection of the operator under discussion appears a line with a unit operator and coefficient 1.234d9. Please add the line 1.234d9  1 to the operator tableau. 
orthogonalise = 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 antisymmetric linear combinations with the aid of Acoeff. Note that the id keyword must not be set in the SPFBasisSection. 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+1i), i=1,...,gdim, for symmetrization and phi(i) = phi(gdim+1i) 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 antisymmetric are antisymmetrized. 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 2Dmode 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 2Dmode 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 Avector! (see the example input "H2H2_nsym.inp" in the "Further Example inputs" section) 
sym2kleg = I (,I1,...) asym2kleg = I (,I1,...) 
I = Number of 4Dmode (2x KLeg/K) which is to be
(a)symmetrized. phi(θ_{1},k_{1},θ_{2},k_{2}) = phi(θ_{2},k_{2},θ_{1},k_{1}) for symmetrization and phi(θ_{1},k_{1},θ_{2},k_{2}) = −phi(θ_{2},k_{2},θ_{1},k_{1}) 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 3Dmode 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 Avector is symmetrised by summing all permutations of its indices. This makes only sense, if all modes are identical, i.e. if the id keyword (SPFBasisSection) 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. It is recommended to set dav for block relaxations with RDAV or RRDAV. If integers I1, I2, .. are given as argument, only a subset 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 Avector 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 (SPFBasisSection) is used for the modes to be antisymmetrized. 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. This should be used as default. ad: adiabatic correction. hh2: use the routines specially written for H+H_{2} reactive scattering using LSTH PES. hh2bkmp2: use the routines specially written for H+H_{2} 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. Note also that the mass for the translational degree of freedom, mass_<modelabel of translationalDOF>, has to be defined in a ParameterSection. 
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. 
Special keywords for build subsection  
Keyword  Description 
init_state = I  Initial electronic state of a wavefunction propagation. 
Print some information to the log file on how multidimensional modefunctions are build from 1Dfunctions. 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 PrimitiveBasisSection. 
K  Body fixed magnetic quantum number for KLeg (or PLeg). 
Wigner  Wignerd function. Requires Wigner in PrimitiveBasisSection. 
map  Use the initial (1D) singleparticle functions of another DOF. May be used to add up to three functions. 
readspf  Read the initial 1D SPF from an ASCII file. 
euclid  The initial 1D SPF is set to unity for one grid point, and zero everywhere else. 
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 singleparticle function will be populated initially. This parameter is optional, the default is p = 1.  
pack = p  The Buildline is for the pth packet in a multipacket 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 expDVR and if the primitive grid is periodic (one of the keywords: periodic, speriodic, 2pi or s2pi 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 singleparticle function. With the aid of the keyword periodic also the region near 2pi gains intensity.  
ran = I  This parameter is necesary only when thermalizing the initial wavepacket. In such case the corresponding degree of freedom will be randomized up to the Ith state of the function. I cannot exceed the primitive basis size and when omitted a default value of ran = basis_size1 will be used. (See Remarks on thermal relaxation calculations).  
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(<x^{2}><x>^{2}). The width parameter may be complex.  
pop = p  The pth singleparticle function will be populated initially. This parameter is optional, the default is p = 1.  
pack = p  The Buildline is for the pth packet in a multipacket run. This parameter is optional, the default is p = 1.  
periodic  Same as for HO.  
ran = I  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 inputfile or via an option on the command line. (operator file definitions come too late). 
l  The lquantum number of the first singleparticle 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 inputfile 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). 

ran = I  Same as for HO.  
sphFBR  j  Quantum number j of the initial spherical harmonic. 
phiFBR  m  Quantum number m of the initial spherical harmonic. 
eigenf  potential  Name of onedimensional potential curve, XX, defined in
the HAMILTONIANSECTION_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 protocolled in
the logfile. NB: eigenf works with real and complex matrices, but the matrix to be diagonalized has to be hermitian. If you want the eigenvectors printed, specify "veigen" in the RUNSECTION. 

ran = I  Same as for HO.  
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 inputfile 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 jstates is preferred (default). s=m: Excitation of mstates is preferred (j goes down). s=mold: Excitation of mstates is preferred (j goes up; old behavior, deprecated). 

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 multiKLeg modes (i.e. modes with more than one KLeg). This sets the number of 1D SPFs which are initially generated; the real (multiD) 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 inputfile 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 primitivebasissection. 

kmax  Maximum value of the body fixed magnetic quantum number
for the initial wavefunction. Default is the corresponding value given in the primitivebasissection. 

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. Note that the angles (beta,gamma,alpha) correspondent to the angular momenta (j,k,m), respectively.  
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). 
factor  Optional. The mapped function is multiplied with the factor factor. If no factor is given, 1 is assumed. If one number is given, R1, a real factor is assumed. If two numbers are given, R1 R2, they are interpreted as real and imaginary part of a complex number. A complex factor may also be given as (R1,R2). A factor is useful if map is used to add functions. (See note below).  
readspf  filename  Name of the ASCII file from which to read the initial 1D SPF(s) (which later may be combined to combined modes). The number of lines in the file must be equal to the size of the primitive basis for this DOF, and each line must contain two real numbers: the real and imaginary part of the SPF. There must be no blank or comment line. If one wants to input the next higher SPF one should write next in the following line and then the data of the next SPF. And so on. Missing higher order SPFs will be generated automatically. For a singleset calculation it is possible to read the SPFs of the electronic degree of fredom. 
addweight  If the optional keyword addweight is given, the function values at the grid points are multiplied with SQRT(weight), i.e. the inputed values are assumed not to be DVR amplitudes, but simply the value of the function at the grid point. (See Eq.(B.25) of the MCTDH review (2000)).  
euclid  point  An integer which spcifies the gridpoint at which the initial 1D SPF ist set to unity. On all other gridpoints the SPF is set to zero. One may use euclid also for an electronic DOF (replacing init_state). 
ran = I  Same as for HO.  
The initial singleparticle functions are formed as follows:
y eigenf Hspf pop=2Then each of the following three Hamiltonians can be used equally well to generate the desired initial singleparticle functions.
HAMILTONIANSECTION_Hspf modes  y 1.0  KE omy  q^2 endhamiltoniansection
HAMILTONIANSECTION_Hspf  modes  x  y  z  el  1.0  1  KE  1  1 omy  1  q^2  1  1  endhamiltoniansection
HAMILTONIANSECTION_Hspf  modes  x  y  z  1.0  KE  1  1 omx  q^2  1  1 1.0  1  KE  1 omy  1  q^2  1 1.0  1  1  KE omz  1  1  q^2 const  q^3  sin  q  endhamiltoniansectionThe lines 1,2 and 5,6 of the tableau of the third Hamiltonian are ignored, as they refer to other degrees of freedom than the y one. The last line is ignored, as it represents a correlated Hamiltonian term.
PRIMITIVEBASISSECTION  # Mode DVR N xi xf  z fft 96 0.75 4.5000 linear x fft 24 0.00 5.3669 periodic y fft 24 0.00 5.3669 periodic theta PLeg 36 even phi exp 9 2pi  ENDPRIMITIVEBASISSECTION INIT_WFSECTION build  # mode type center moment. freq. mass  z HO 1.80 24.d0 694.1683673,eV 1.d0 x HO 0.00 0.0 0.0 1.d0 # flat function y HO 0.00 0.0 0.0 1.d0 # flat function theta KLeg 2 sym # j of initial spherical harmonic phi K 1 # m of initial spherical harmonic  endbuild endinit_wfsection
For Wigner initialwf:  DOF init_wf allowed primitivebasis types     1 wigner  wigner only 2,3 K only  K or exp  For Wigner primitivebasis:  DOF pbasis  allowed initialwf types     1 wigner  wigner  gauss or leg  any 1D type (HO,gauss,leg,eigenf,map) 2,3 K or exp  K  K (for pbasis=K only)  HO,gauss,eigenf,map (pbasis=exp only) An example involving the use of Wignerinwf:
PRIMITIVEBASISSECTION  # Mode DVR N center freq mass  r HO 8 1.626 0.01837 1739.96 beta wigner 16 all gamma exp 15 2pi m k 7 7  ENDPRIMITIVEBASISSECTION INIT_WFSECTION build  # mode type q.n. center mom. freq. mass  r HO 1.426 0.0 0.01837 1739.96 beta wigner 10 nosym excite=mkj print # initialj, use odd and even j gamma k 1 2 2 1 # initialk, kmin, kmax, dk m k 2 5 5 1 # initialm, mmin, mmax, dm  endbuild ENDINIT_WFSECTION
PRIMITIVEBASISSECTION . . R FFT 192 0.6 3.0 R_dum sin 60 0.5057591623 0.6062827225 endprimitivebasissection INIT_WFSECTION build . . R map R_dum R_dum eigenf operator pop=1 endbuild endinit_wfsectionSince the DOF R_dum (in our example) is not defined in the SPFBASISSECTION it will be ignored, when the operator is build. To avoid this, one has to use the addmode keyword.
HAMILTONIANSECTION_operator addmode = R_dum  modes ......  R_dum  ....  . .The map keyword can be used to sum up to three functions. A corresponding INIT_WFSECTION may read:
INIT_WFSECTION build rd map dum1 1.0 rd map dum2 1.5 rd map dum3 (0.5,0.3) rv HO 2.151 0.0 0.218,eV 13615.5 theta gauss 2.22 0.0 0.0745 dum1 gauss 4.310 0.0 0.075 pop=1 dum2 gauss 4.315 0.0 0.080 pop=2 dum3 gauss 4.320 0.0 0.085 pop=3The initial SPF for rd is a weighted sum of three harmonic oscillator functions. The first two coefficients are real, the third is complex. After summation, the function is normalized. The DOFs rd, dum1, dum2, and dum3 must be defined in the PRIMITIVEBASISSECTION with identical grids.
Remarks on multiset calculations:
For multiset calculations, the initial functions can be defined differently for each state by using the state keyword:
modelabel type parameters state = sThe 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 multipacket calculations, i.e. when packets > 1, the initial function definition must be given for each packet as
modelabel type parameters pack = pwhere 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 singleparticle function can be selected using the pop keyword:
modelabel type parameters pop = iwhich populates the i th function.
Remarks on thermal relaxation calculations:
Which degrees of freedom are to be thermalized must be specified through the ran keyword in the INIT_WFSECTION according to the order of parameters given in the table above labelled "parameter description". In the following we show four possible illustrative cases with an example. For the first mode the Ith state up to which the randomization will be averaged is not specified, thus the default value basis_size1 will be employed. Meanwhile for the second mode, I is declared with a value of 30. Following this in the case of the third mode, the spf type eigenf is used together with the ran keyword. Finally, when modeeigenfunctions are generated, the ran argument has to be included in the respective parameter's list line (meigenf keyword) as it is the case for the fourth mode. The meigenf procedure needs an initial SPF generated in the buildblock. This initial SPF should not be an eigenfunction of the operator used in meigenf, but should have an overlap with several eigenfunctions. For this reason we have dispaced the origion of the HOfunction by 0.25 au. Independently of the type of mode, the ran keyword will always be the last argument. Note that an electronic degree of freedom cannot be randomized. Electronic states are always considered as part of the system, not of a thermalized bath.
Example:
INIT_WFSECTION build init_state = 1 q0001 HO 0.00 0.0 1.0 ran q0002 HO 0.00 0.0 1.0 ran=30 q0003 eigenf oper ran=30 q0004 HO 0.25 0.0 1.0 endbuild meigenf = 4,oper,0,full,noselect,write,ran=30 endinit_wfsection
Reading the initial wavefunction ReadInwf ... endreadinwf  

Keyword  Description  
file = S  The string S is a path of a directory 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 Avector 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 Avector 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 transformed to natural orbitals and not reorthonormalized  
realpsi  Only the real part of the wavefunction is used and Schmidtorthonormalised as in the orthopsicase.  
ignore  The different primitive bases error will be ignored. 
Note: Here all wavefunctions are assumed to be of multiset format
Examples:
For a two state WF the following input is equivalent to a simple
file = <restartdirectory> statement.
ReadInwf file = <restartdirectory> orthopsi # This is default, opposite: noorthopsi SPF 1 > 1 SPF 2 > 2 A 1 > 1 A 2 > 2 endreadinwfThe following two inputs are equivalent. In both cases the SPFs of the first state of the gscalculation (which may have only one state) is put on state 1, 2 and 3 of the initial WF. The Avector of the first state of the gscalculation is put on the second state of the initial WF. The Avector for states 1 and 3 is zero.
ReadInwf file = gs SPF 1 > 1 SPF 1 > 2 SPF 1 > 3 A 1 > 2 endreadinwf ReadInwf file = gs SPF 1 > 1,2,3 A 1 > 2 endreadinwfIf 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.
ReadInwf file = <restartdirectory> SPF 1 > 1,2 init_state = 2 endreadinwfOne may use more than one (up to eight) restart files to build the initial wavefunction.
ReadInwf file = run1 SPF 1 > 1 file = run2 SPF 1 > 2 A 1 > 1 file = run3 A 2 > 2 endreadinwfIt is possible to multiply the Avector 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.
ReadInwf file = <restartdirectory> SPF 1 > 1 SPF 2 > 2 A 1 > 1 weight = 0.65 A 2 > 2 weight = (2.3, 0.96) endreadinwfIn the first Aline, weight is assumed to be real, i.e. the input is equivalent to weight=(0.65,0.0).
ReadInwf file = run1 SPF 1 > 2 endreadinwf
INIT_WFSECTION 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 endbuild Acoeff 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) endAcoeff endinit_wfsectionIf the block size grows, the definition of the Avector with Acoeff becomes rather cumbersome. Here autoblock is very useful, as it generates zero order Avectors for the lowest states. E.g.:
INIT_WFSECTION 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 endbuild autoblock endinit_wfsectionNote that the autoblock command is executed in the propwf step. Hence a geninwf run alone will yield a zero Avector.
INIT_WFSECTION blockSPF = GS autoblock endinit_wfsectionHere GS denotes a directory which contains a restart file (e.g. namedirectory 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 blockSPF = 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. blockSPF=GS,realpsi. By default the SPFs are reorthonormalized. To switch this feature off, set the option noorthopsi.
INIT_WFSECTION blockSPF = GS/rst000,noorthopsi blockA = GS/rst001,GS/rst003,GS/rst005 blockA = GS/rst007,GS/rst009,GS/rst011 endinit_wfsectionThere may be more than one blockA statement. The number of Avectors read must match the block size. It is important to set the option noorthopsi here, because otherwise the SPFs will be transformed and will no longer be consistent with the Avectors. One usually wishes to set relaxation=follow when Avectors were read in (relaxation=lock is not implemented for block improved relaxation).
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 timestep 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 singleparticle 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 timestep control. 
CMF/varphi = R, R1  The mean fields are kept constant over a variable time
interval, determined by the error of the singleparticle
functions only. (Not available for multilayer runs, as ML
doesn't distinguish between A and SPFs.) 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 timestep control. 
CMF/vara = R, R1  The mean fields are kept constant over a variable time
interval, determined by the error of the MCTDH coefficients
(Avector) only. (Not available for multilayer runs, as ML
doesn't distinguish between A and SPFs.) R is the initial time interval (in fs). R1 is an accuracy parameter for the timestep control. 
MLCMF = S  When CMF is set for a MLcalculation, then the keyword MLCMF must follow. The argument S can take the values unified or split. See note below. 
The following keywords define the integrator to be used.  

Keyword  Description 
ABM/S = I, R, R1  AdamsBashforthMoulton predictorcorrector integrator
used for S = all, spf, A. I = order R = accuracy R1 = initial stepsize 
BS/S = I, R, R1  BulirschStoer extrapolation integrator used for S = all,
spf, A. I = maximal order R = accuracy R1 = initial stepsize 
RKn/S = R, R1 (,S1,S2)  RungeKutta integrator of fixed order n = 5 or 8,
used for S = all, spf, A, or modespec (see below). R = accuracy R1 = initial stepsize (can be omitted or set to zero, then the integrator tries to guess a suitable value) S1 = importho : Orthonormality of the SPFs is improved after each RungeKutta step. For multilayer runs, this does not work in VMF mode. S2 = rhonorm : The error is determined as a densityweighted L2norm of a difference SFPvector. Unoccupied SPFs then do not limit the step size. For MLruns setting rhonorm is allowed only for CMF and mlcmf=split integrator setting. If rhonorm is set it is advisable to set importho as well. 
SIL/S = I, R, S1  SIL integrator used for S = all, spf, A, @1. 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  complexSIL integrator used for S = all, spf, A, @1. Same as SIL/S (see above), but the use of the (complex) LanczosArnoldi integrator is enforced. (SIL tries to make the choice Lanczos/LanczosArnoldi 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 nonhermitian Hamiltonians (resonances). 
The string S=all after the integrator name may be omitted, i.e. ABM is equivalent to ABM/ALL.
For VMF calculations, only the BS, ABM or RK5/8 integrator may be used for the differential equations. Default is ABM. For VMF calculations the integrator must carry the extension S=all (or no extension at all), i.e. there is only one integrator within the VMF scheme.
For CMF calculations with standard (not multilayer) MCTDH, the following combinations of integrators are possible: ABM/spf + SIL/A, BS/spf + SIL/A, RKx/spf + SIL/A, BS/all, ABM/all, RKx/all. Default is BS/spf + SIL/A.
For CMF calculations with multilayer MCTDH, there are two different schemes of integrating the equations of motion. The scheme can be selected by the keyword mlcmf (see below). 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 MLBASISSECTION. One can doublecheck the mode numbers by inspecting the visualization of the ML tree generated by the graphviz keyword [RUNSECTION]. For example inputs see test38.inp and test39.inp which are on $MCTDH_DIR/elk_inputs.
For numerically exact calculations (i.e. the keyword exact has been specified in the RUNSECTION), 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 "Apropagation" 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 (only real operators, e.g. no p).
NB: The propagation of the SPFs is always done
in complex arithmetic. cDAV is for nonhermitian 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 blockrelaxations. To this end set the packets
keyword in the SPFBasisSection. Note that cDAV cannot be used in
block form. All Davidson routines allow for parallelization
(pthreads and/or mpi).
NOTE: There are two Lanczos integrators, a real version for hermitian Hamiltonians and a complex version (LanczosArnoldi) for nonhermitian 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 nonhermitian, one should use CSIL. The use of SIL in case of a nonhermitian Hamiltonian may lead to wrong results.
The following keywords may be used to further define the method of propagation.  

Keyword  Description 
eps_svd = R  R is the value used to regularize the EOM via an SVD as described in Meyer and Wang, JCP 148 124105 (2018). If eps_svd=0 (which is the default), then no SVD, but normal regularization via regularizing the density matrices is performed. 
eps_inv = R  R is the value used to regularise the inverse of the reduced density matrices. (Default: 10^{8}. See eq.(82) review.) When SVD regularization is used and the lowest natural population increases eps_inv, then SVD regularization is switched off and the program returns to normal regularization. 
eps_no = R  R is the value used to regularise the "natural orbital Hamiltonians". (Default: 10^{8}) 
projh  One dimensional Hamiltonians are not extracted to be in front of the projector. (See eq.(45) review). This is default for CMF, relaxation, and MLruns. 
hproj  One dimensional Hamiltonians are extracted to be in front of the projector. (See eq.(44) review). This is not implemented for MLruns. 
natorb  Natural orbitals are propagated in place of spfs. When the CMFintegrator is employed, the spfs are propagated normally, but the whole wavefunction is transformed to natural orbital picture after each CMF step. In case of a relaxation run the orbitals are reorthonormalised after each CMF step. (Not supported for multilayer runs.) 
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 (RunSection). (Not supported for multilayer runs.) 
stdorb  Standard orbitals are propagated. This is default, except for improved relaxation using DAV. 
interpic  Spfs are propagated using the interaction picture. interpic is not allowed for CMFintegration. 
simpleproj  The simple projector is used rather than the improved one. The inversion of the spfoverlapmatrix is thus avoided. 
nohsym  The symmetry of the operators determined by the program is not used to calculate the operator matrix elements. 
CDVR  Multidimensional potential terms are evaluated using CDVR . The analytic_pes keyword needs to be set in the OPERATORSECTION 
TDDVR  Multidimensional potential terms are evaluated using TDDVR. The analytic_pes keyword needs to be set in the OPERATORSECTION 
directmultid  Multidimensional potential terms are evaluated using a direct algorithm, i.e. the potential is not stored on the full grid but recalculated each time it is needed. This is much slower, but needs minimal memory. The analytic_pes keyword needs to be set in the OPERATORSECTION 
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). 
importho  Reorthonormalize all SPFs after each step (VMF: output step; CMF: update step). The long form improvedorthonormality is also possible. Note: This currently only works for multilayer runs. For standard MCTDH, the RungeKuttaspecific importho is tried instead (see above). 
mlcmf = S  For multilayer runs in CMF mode, the SPFs of each mode can be propagated together (S = unified) or separately (S = split). In the latter case, the choice of integrator is currently restricted to RK5/8 or (C)SIL. Split propagation is usually faster, because each mode can adjust its stepsize independently. The drawback is a somewhat increased usage of memory. (Default: S = unified) 
Omitted integrator parameters default to:
Keyword  Default  

CMF/fix  R = min(1.0,tout)  
CMF/var  R = min(1.0,tout)  R1 = 1.0e6  
CMF/varphi  R = min(1.0,tout)  R1 = 1.0e6  
ABM  I = 6  R = 1.0e5  R1 = 1.0e4 
BS  I = 8  R = 1.0e6  R1 = updatetime/2 for CMF or R1 = 0.2 for VMF 
SIL  I = 30  R = 1.0e6  S1 = standard 
RK5  R = 1.0e6  R1 = 0 (i.e. autoguess)  
RK8  R = 1.0e6  R1 = 0 (i.e. autoguess)  
projh / hproj  in VMF mode: hproj in CMF mode: projh multilayer: projh for relaxation: projh 

eps_inv  R = 1.0e8  
eps_no  R = 1.0e8 
List of parameters of ALLOCSECTION.  

Parameter  Example (NOCl)  Default value  Description 
maxdim  5  autom. detec.  Number of DOFs+2 
maxsta  1  autom. detec.  Number of multiset electronic states 
maxtape  0  autom. detec.  Length of MLtree 
maxdepth  0  autom. detec.  Depth of MLtree 
maxham  3  autom. detec.  Maximal number of of Hamiltonians 
maxpar  203  autom. detec.  Maximal number of parameters 
maxpes  1  autom. detec.  Maximal number of potential energy surfaces 
maxrdf  0  autom. detec.  Maximal number of rdfiles (e.g. natpots) 
maxkoe  160  4000  Maximal number of Hamiltonian terms 
maxhtm  220  4000  Maximal number of operators for building the Hamiltonian(s) 
maxhop  220  4000  Maximal number of labels used in operators 
maxfac  25  200  Maximal number of factors/terms in an operator 
maxdef  85  240  Maximal number of operators that can be defined 
maxmuld  1  10  Maximal number of multidimensional operators 
maxsub  40  100  Maximal number of subroutines used in timing analysis 
maxLMR  0  20  Maximal allowed number of LMRterms 
maxedim  0  50  Maximal number of eigenvalues in in eigenf (veigen) and in adiabatic correction 
maxreadspf  0  10  Maximal number of readspf keywords (buildsection) 
maxnhtmshift  1  10  Maximal number of different htmshift values 