# The Analyse Programs

## Generalities

Analyse is a set of programs able to extract data of interest from the results of an MCTDH calculation. In essence, there is a library of modules providing an interface to the MCTDH program output. These modules can be then built into programs as required. The documentation is in two parts: Firstly the various programs for computing observable quantities are described, and secondly the various shell-scripts for plotting the results using the Gnuplot program are detailed.

The programs can be compiled using the shell script \$MCTDH_DIR/bin/compile. Type:

 compile progname
to produce the executable \$MCTDH_DIR/bin/\$MCTDH_PLATFORM/ progname<ver>, where progname is the name of one of the analysis programs (e.g. overlap or autospec). Typing
 compile analyse
will compile all the analysis programs.

The compile script has various options (type compile -h for a list). For example, to produce the exectuable \$MCTDH_DIR/bin/\$MCTDH_PLATFORM/overlap<ver> d (e.g. overlap84d) with debug information, type

 compile -d overlap

NB: All programs are usually compiled when the package is installed.

## Documentation of the Individual Programs

All program names have the version number apended to them, e.g. the overlap program from version 8 release 3 is called overlap84. The following programs exist:

• adpop : Calculates the diabatic and adiabatic populations of the electronic states. The calculation of one and two dimensional adiabatic densities is also possible.
• adproj : Calculates vpot files for the adiabatic potential matrix elements and the adiabatic projector matrix elements for each electronic state.
• analysis : An interface program that uses a set of menus to interactively analyse the results of calculations.
• aoverlap : Calculates the difference between two A-vectors.
• autospec : Computation of the spectrum by fourier-transforming the autocorrelation function.
• chnpot : Interpolation of natpot terms to another grid.
• compare : Used by the elk_test program to compare different sets of results.
• concat : Reads two restart files and concatenates the wavefunctions, i.e. psi(q,x)=psi1(q)*psi2(x).
• convpsi : Converts a multi-packet wavefunction to the corresponding single-packet ones.
• crosscorr : Reads two MCTDH wavefunction psi and psi_ref and performs the cross-correlation function c(t) = <psi_ref(0)|psi(t)>. Matrix-elements <psi_ref(0)|oper|psi(t)> can also be done.
• crosspec : Fourier-transforms a cross-correlation function generated by crosscorr.
• dengen : Generates 2D-densities from a wavefunction.
• diag : Diagonalizes the Hamiltonian matrix calculated by fmat.
• edengen : Generates (2D) pair-correlation functions, , from the wavefunction.
• expect : Calculates the time evolution of the expectation value of an operator. Reads the output file psi.
• fdcheck : Checks whether a state has been found in all filter data sets.
• fdmatch : Finds the optimal matches between sets of eigenvalues, where each set is obtained with the (real) filter program.
• flux : Computes the quantum flux by analysing the matrixelements of the CAP. Computes reaction probabilities.
• fmat : Computes the matrix elements <psi(t1) | Operator | psi(t2)>. The usual hermitian or the symmetric scalar product may be used. The diagonalisation of the matrix may then be performed by diag.
• jinpol : Interpolates flux-probabitities for a range of total J.
• joinpsi : Joins psi files have been calculated with different number of single particle functions.
• maxpe : Calculates the minimum and maximum value of a PES.
• mmemtest : Checks how much memory is available for MCTDH.
• norm : Calculates the norm and A-norm of a wavefunction.
• npotminmax : Like vminmax. Reads a set of natpot files and runs over the total product grid or a selection of grid points to determine the maxima and minima of a potential energy surface.
• ortho : Checks the orthonormality of the spf's and spdo's. Checks also the hermiticity of the spdo's.
• overlap : Calculates the difference between two wavefunctions.
• pexpect : Calculates the time evolution of the expectation value of a one-particle operator. Reads the output file pdensity.
• prhosub : Reads pdensity file and writes reduced density of one selected DoF on primitive grid to files.
• probsq : Reads the auto file and computes the average of |c(t)|2, which is the sum of the squared occupation probabilities of the eigenstates.
• psi2ex : Converts an MCTDH psi file to an exact (full grid) representation. Reads the output file psi and creates psi.ex.
• psi2rst : Writes a WF from an psi file to a restart file.
• rdacoeff : Reads the A-vector and prints it to screen.
• rdautoe : Analyses the error within the autocorrelation function.
• rdcheck : Analyses the state population and natural weights.
• rdgrid : Reads the DVR grid points from dvr file and prints them to screen.
• rdgpop : Analyses the grid population read from file gridpop.
• rdsteps : Reads the step-file and evaluates statistic measures of the integrator steps.
• rdnatpot : Reads a natpot file and prints information on it to rdnatpot.log.
• rdupdate : Reads the update steps from file "update" and enables plotting of the data.
• rdauto : Reads the auto steps from file "auto" and enables plotting of the data.
• rdeigval : Reads the eigenvalues from file "eigval" and enables plotting of the stick spectrum.
• showd1d : Enables plotting of the 1-dimensional densities with time.
• showpot : Enables plotting of a PES fit from a natpot file. Comparison with true potential can be made if vpot file is present.
• showrst : Enables selective plotting of the single-particle functions from the restart file.
• showspf : Enables selective plotting of the single-particle functions from the psi file.
• showsys : Enables 2D plotting of the wavefunction or the PES . Overlay plots, i.e. WF on PES, are possible as well.
• sumrst : Allows calculation of the superposition of up to 8 wavefunctions with different weight.
• sumspec : Sums spectral components (output from autospec).
• twopdens : Computes the 2-particle density in the basis of the SPFs.
• twprob : Computes state-resolved reaction probabilities with the method of Tannor/Weeks.
• vminmax : Reads a pes-file and runs over the total product grid or a selection of grid points to determine the maxima and minima of a potential energy surface.
• wfproj : Computes the projection of a wavefunction onto a (combined) mode.

Below the options of each analyse program are described. This information can also be obtained by using the -h option, e.g.

 aoverlap84 -h
The information given there is a little bit more reliable since the on-line documentation might not be up to date in some cases.

Purpose: Calculates the diabatic and adiabatic population.
And also one or two dimensional adiabatic densities.
Usage:   adpop84 [OPTIONS] [Degrees of Freedom]

Options  :
-i DIR  : data is stored in directory DIR
-f FILE : the PSI is read from FILE rather than from ./psi
The string FILE may be a relative or a full path-name.
-d FILE : the DVR is read from FILE rather than from ./dvr
The string FILE may be a relative or a full path-name.
-p FILE : the PES is read from FILE rather than from ./pes
The string FILE may be a relative or a full path-name.
-o FILE : The ouput is written to file FILE rather than to ./adp
The string FILE may be a relative or a full path-name.
If FILE=no, i.e. -o no, no output file will be written.
-r      : The PSI is read from restart file rather than from ./psi
-n num  : only the first num WFs are read  from psi-file.
-skip m : the first m WFs are skipped.
-step step : Compute the population only every step-th WF.
-tcut tcut : Compute the population only at time tcut.
-w      : Overwrite the output-file adp.
-q      : The quick algorithm is used.
-mc     : Monte-Carlo integration is used.
-noq    : Quick algorithm not used with MC.
-qtol qtol   : The quick algorithm cutoff set to qtol.
-mctol mctol : The Monte-Carlo convergence set to mctol.
Defaults: qtol=1.0d-6, mctol=1.0d-3
-sleazy : qtol and mctol relaxed   to:  qtol=1.0d-4, mctol=1.0d-1
-tight  : qtol and mctol tightened to:  qtol=1.0d-8, mctol=1.0d-3
-ver    : Version information about the program
-h      : Print this help text.
-?      : Print this help text.

A pes file that can be generated using mctdh with the -pes option is needed.
If an exact WF is used then the pes file must be in exact representation too!

The arguments "Degrees of Freedom" have to be replaced by the names of the
degrees of freedom for which densities are to be calculated.
To calculate 2-dim densities separate the two DOF with a comma.
------------------------------------------------------------------------------

The adpop84 program is used to calculate the populations of the electronic states. This is done for the diabatic and for the adiabatic wavefunction. The results are stored in the file adp for each time step. The data can be plotted with plgen. Also one and two dimensional densities can be calculated with adpop84. The densities are stored in the adp_dof* and adp_dof*_dof* files where dof* is replaced by the name of the corresponding degree of freedom. The densities are calculated with and without weights.
NB: adpop needs a pes file that can be generated with mctdh using the -pes option. Further adpop need a dvr and a psi file. NB: adpop works only for multi-set calculations.
NB: As the calculation runs over the primitive grid, the calculation is slow. Analysing one WF takes about [grid-dimension]*[A-vector length]*1.5*10-8 s on a 3 GHz P4. The calculation can be accelerated by setting the quick option, -q, or by using Monte-Carlo integration, -mc. Using the quick option, one ignores points for which the product of one-dimensional grid-populations (1D-densities) is lower than qtol. The Monte-Carlo integration is, of course, less accurate. For small systems (3D, say) it even may be slower than the direct method, but it allows to treat systems which are otherwise too large to be analyzed with adpop.
NB: The adpop population file can be easily plotted with plgen.
NB: For the adpop density files the script pladpop exists. This script plots the adiabatic densities of the different electronic states with or without weights.

projector.

Options  :
-o DIR  : Output data is stored in directory DIR rather than in current one.
-p FILE : the PES is read from FILE rather than from ./pes
The string FILE may be a relative or a full path-name.
-w      : Overwrite the output files.
-nopot  : Adiabatic PES are not calculated.
-noproj : Projectors will not be calculated.
-s s    : Only adiabatic PES and Projector Matrix Elements
for state s will be calculated.
-me i j : Only Projector Matrix Elements P_i_j will be calcualted.
-ver    : Version information about the program
-h      : Print this help text.
-?      : Print this help text.
------------------------------------------------------------------------------

The adproj84 program is used to calculate the adiabatic potential matrix elements and the adiabatic projector matrix elements for each electronic state. The program must be run in a directory containing a dvr-file and a pes-file (usually the directory of the genpes-run). The results are written into several vpot files that can be fitted with the potfit program. To do this the keyword "pes = none" has to be set in the potfit input file. With a MCTDH-run the projection operator file can be created. To receive the adiabatic populations run expect with the projection operator file and the corresponding psi file.
NB: The vpot files for the adiabatic potential are: apr_vs, where s is the electronic state.
NB: The vpot files for the projection matrices are: apr_ps_ij, where s is the adiabatic electronic state and ij denotes the matrix element, i.e. the diabatic states. Please note that the projection matrices are symmetric, this means that each off-diagonal matrix element is only given once.

aoverlapver

Purpose: Calculates the difference between two A-vectors.
See review Eqs.(159-163) for a definition of the output

Usage  : aoverlap84 [-f -i -o -n -w -ver -h -?]  [psi2] .

Options :
-f FILE : The wavefunction is read from file FILE rather than from ./psi
The string FILE may be a relative or a full path-name.
-i FILE : Same as -f FILE.
-o FILE : The difference is written to file FILE.pl rather than to ./aovl_XXX
The string FILE may be a relative or a full path-name.
If FILE=no, i.e. "-o no", no aovl-file will be opened.
-abs    : Take the absolute value of A, when computing the overlap.
-n step : Compute the A-overlap only every step-th output.
-w      : An existing aoverlap file is overwritten.
-ver    : Version information about the program
-h      : Print this help text.
-?      : Print this help text.
The argument psi2 is the name of the file containing the comparison
wavefunction.
If no argument is given, the user is prompted for the missing argument.
----------------------------------------------------------------------------

The aoverlap84 program is very similar to the overlap84 program, except that it considers the A-vectors exclusively, i.e. it implicitly assumes that the spf's of the two calculations are identical. aoverlap84 is useful when the two calculations to be compared use different DVR's but identical numbers of spf's. When the overlap between two calculations is OK but the aoverlap is not then the two calculations use different spf's, e.g. a different position of the projector. See also overlap84 below.

autospecver

Purpose: Computation of the spectrum by Fourier-transforming
the autocorrelation function.
Three spectra are generated acording to the weight functions
cos^n(pi*t/(2*T)) (n=0,1,2). The option -lin allows to use
a second set of filters. In plspec the filters of this
second set are enabled through the option -g3, -g4, or -g5.
Additionally there is the weight function exp(-(t/tau)^iexp).
The conversion factor from damping time to FWHM energy width
is 1.32 eV fs and 2.2 eV fs for iexp=1 or 2, respectively.
Input a zero for tau if no damping is required.

Usage : autospec84 [-f -i -o -e -n -r -ph -p -q -g -t -FT -EP -Mb -ver -h -?]
[emin emax unit tau iexp]

Options :
-f FILE : The autocorrelation is read from file FILE rather than from ./auto
The string FILE may be a relative or full path-name.
-i DIR  : Take DIR as input directory.
-o FILE : The spectrum is written to file FILE.pl
rather than to ./spectrum.pl
The string FILE may be a relative or full path-name.
-e R Unit : the offset energy is set to R in units Unit; e.g. 0.3 ev
Note, a positive R shifts the spectrum to the right.
-n      : The maxima of the three spectra are normalized to 1.
-r      : Plot the resolution rather than the spectrum
-ph  ph : Use phase correction exp(i*ph*1.d-6*t^2)
-p I    : The no. of plot-points is I rather than 1000
Use "-p 2000" for final plots
-q I    : Same as -p.
-lin    : A second set of filter functions is used, namely
g(t) = cos^3(pi*t/(2*T) (column 2),
g(t) = 1 - t/T (column 3),
g(t) = (1-t/T)*cos(pi*t/T) + sin(pi*t/T)/pi (column 4).
The last two filters yield a positive semi-definite
Fourier-transform.
-g ncos : GNUplot command lines are written to the output file.
ncos is the exponent of the cosine damping function.
-t tcut : The autocorrelation function is read only up to t=tcut
-FT     : No multiplication with the energy prefactor. Fourier
transform only. (default!)
-EP     : Multiplication with the energy prefactor.
-Mb norm: Absorption spectrum in Mb. norm is ||D*Psi|| (in au)
where D is the Dipole moment. -Mb sets -EP.
-ctr    : Perform a sin-transform rather than an exp one.
(This is for ctrace files.)
-ver    : Version information about the program
-h      : Print this help text.
-?      : Print this help text.
Following the options one may provide the arguments (separated by blank):
Emin, Emax, Unit, tau, iexp
If not enough arguments are given, the user is prompted for the missing
arguments. Use the Unit 'no' if the propagation is run with the
keywords: 'time-not-fs','energy-not-ev'.
When using the options -EP or -Mb one must make sure that the zero point
of the energy is set appropriately. Use the option -e if the total energy
differs from photon energy. 1 Mb = 10^-18 cm^2. The norm ||D*Psi|| can be
found in the log file of the operate run.
------------------------------------------------------------------------------

The program computes the spectrum by Fourier transform of the autocorrelation function which is multiplied with the weight function exp(-(t/tau)iexp)*cosn(pi*t/2T), (n = 0, 1, 2). (In the output file there is one column for each n). The variables tau (real [fs]) and iexp (integer) are input arguments.
More precisely:
sigma(omega,n) = prefac*Integral (0, T) Re[ c(t) * exp(i*(omega -Eoffset)*t)] * exp(-(t/tau)iexp)* cosn(pi*t/2T) dt
is the formula which is evaluated. Here c(t) is the autocorrelation function and T denotes the largest time for which c(t) is known. (T may be reduced by the option -t). The variable omega runs between Emin and Emax. The pre-factor prefac is 1/pi if -FT is set (default), or is omega/pi if -EP is set, or is 0.4279828*dipolemoment**2 if -Mb is set. In the latter case, dipolemoment is the norm of the D*Psi, where D denotes the dipole operator. This norm is printed in the MCTDH log file, when the operation D*Psi is performed by operate. The -Mb option allows to plot the absorption spectrum on absolute scale in mega barns (Mb). 1 Mb = 10-18 cm2

When the option -lin is set, then the filter function cosn(pi*t/2T) (n=0,1,2) is replaced by cos3(pi*t/2T) (column 2), 1 - t/T (column 3), and (1-t/T)*cos(pi*t/T) + sin(pi*t/T)/pi (column 4). Note that the latter two filters are non-negative. A more comprehensive discussion can be found in the lecture notes Introduction to MCTDH.

When one of the options -EP or -Mb is set, one must ensure, that the energy is the photon energy. It may hence be necessary to shift the energies from total ones to photon energies by using the -e option.

If the option -r is set, the program computes the resolution, i.e. the spectrum of a single line. The time points are taken from the auto file. If there is no auto file and if the option -t is given, a time step of 1 fs is assumed.

The pre-factor omega is omitted by default, or when the option -FT is given. In this case (and only in this case) the energy Emin (and eventually Emax) may be negative. NB: When the pre-factor is omitted (default), the spectrum is normalized to 1 au.
Example: autospec84 -e 13.4071 ev 13 17 ev 50 2

chnpotver

Purpose : Interpolation of natpot terms to another grid.
Usage   : chnpot84 [-D -d -rd -w -mnd -n -nd -ver-h -?] input

Options  :
-D PATH  : The output is directed to PATH.
This option overwrites name in the run-section of the input file.
-d DIR   : The dvr is read from file DIR/dvr. This option overwrites
dvrdir in the run-section of the input file.
-rd FILE : The new dvr is read from file FILE. This option overwrites
readdvr in the run-section of the input file.
-w       : Overwrite an existing  "natpot" file.
-mnd     : Make Name Directory. Creates the name directory, if not existent.
-nf FILE : The natpot is read from file FILE. This option overwrites
natpotfile in the run-section of the input file.
-nd DIR  : The natpot is read from file DIR/natpot. This option overwrites
natpotdir in the run-section of the input file.
-ver     : Version information about the program
-h       : Print this help text.
-?       : Print this help text.

The argument "input" is the file containing the Run-Section, the new
Primitive-Basis-Section written in MCTDH form and the Fit-Section.
----------------------------------------------------------------------------

The program interpolates natural potential terms calculated by program Potfit (written in "natpot" file) from one grid to another. Normally one runs Potfit on an initial coarse grid, then interpolates the natpot terms to a finer grid. The program reads the old grid from "dvr" file. The input file consist of three sections, a Run-Section, a Primitive-Basis-Section and a Fit-Section (see below) and the final keyword "end-input".

In Run-Section following keywords are possible:

Keywords of chnpot RUN-SECTION
Keyword     compulsory / optional Description equivalent option
name = S   C The name of the calculation. A directory with path S is required in which files will be written.   -D
dvrdir = S   C The old DVR file is located in directory S.   -d
natpotdir = S   O The old natpot file is located in directory S. If this keyword is not given, natpotdir=dvrdir is assumed.   -nd
natpotfile = S   O The path of the old natpot file is S. If this keyword is not given, natpotfile=natpotdir/natpot, or natpotfile=dvrdir/natpot is assumed.   -nf
readdvr = S   O The new DVR should be read from directory S. Without argument, S = name is assumed. Without this keyword the dvr file is created in the name directory.   -rd
overwrite   O Files in name directory will be overwritten.   -w

The Primitive-Basis-Section is as in MCTDH.

In the Fit-Section the parameters for interpolation are given. For translational degrees of freedom one can use cubic (bicubic) spline interpolation or ENO (Essentially Non Oscillatory) interpolation. For angular degrees of freedom (in the case of periodicity) one can select between cosine, sine or fourier interpolation. At present, spline interpolation is implemented for combined modes with an arbitrary number of degrees of fredom while cosine, sine and fourier interpolation are restriced to two dimensions. For more the two dimensions, the spline interpolation is applied recursively to the underlying one-dimensional grids. Note, that the interpolation method should be defined for each degree of freedom, and not for each mode. If grid points for some degree of freedom have not changed, it can be omitted.

The keyword "spline" ("spl") is used for spline interpolation. ENO interpolation is currently implemented for non-combined modes. It consists basically on a polynomial interpolation where the points from the original data set used to obtain the interpolant function are automatically selected to fulfill some smooth criterium. The algorithm will provide a continuous non-oscillatory interpolation near regions with sudden changes in slope, as encountered for example when an energy cut is found in the original data. The formulation of the algorithm can be found in J. Comput. Phys. 71, (1987), 231-303. ENO interpolation is requested using the keyword eno, which can optionally be followed by the keywords rmax=integer, maximum degree of interpolation (default 3), rmin=integer, minimum degree of interpolation (default 1) and dfac=real, weight factor controling how preferential is the symmetric interpolation with respect to the asymmetric one (default 50.0). The ENO implementation is a new feature and by now is rather experimental. Use it with care. The defaults should provide a reasonable behaviour of the interpolator, specifically, setting a higher maximum order or lowering dfac too much may result in a noisy potential energy function at the end. Increasing too much dfac (say, above 200) will tight the interpolator always to a symmetric interpolation, thus not being able to avoid discontinuities. A good idea would be to compare it for a given case to the corresponding spline interpolation and pick up the most satisfactory result.

The keywords "cosine" ("cos"), "sine" ("sin") and "fourier" ("ft") are used for the Fourier-interpolation of angular degrees of freedom. (NB: "ft" or (equivalently) "fourier" means that sine and cosine functions are used). The keyword "sym" after "sin", "cos" and "ft" can be used, indicating that only even functions (i.e. Cos2x, Cos4x, ...) should be used. The angular fit procedure is described in JCP 104,(1996) 7974. The parameter M (default is max(50,2*gdim)) can be given as a keyword following the sin, cos, or ft keyword (see Eq.(37) in JCP 104,(1996) 7974). Finally, the keyword "period=xx" defines the period of the potential curve, for which xx="2pi/N" or "P" format can be used, where N is an integer and P is a real number. If the keyword period is missing, period=2pi is taken as default. If the variable under discussion is not an angle, period should be set to the interval of periodicity. If the potential is not periodic, period should be set to twice the grid length.

The order of the arguments which follow the method keyword is arbitrary.
Here is an example of a Fit-Section:

Fit-Section
R spline
r fourier 120 period=7.85
q eno rmax=3 dfac=60
Q eno dfac=70 rmin=1
theta cosine 60 period=2pi/2
phi fourier sym period=2pi
end-fit-section
Here it was assumed that V(theta=0) = V(theta=pi), e.g. H2O in Jacobian coordinates. If this periodicity is not given, e.g. H2O in valence coordinates, period=2pi should be used.

Note, that the interpolation method for the degrees of freedom in combined mode should be the same (for angular modes cosine, sine and fourier can be mixed). Some parameters of the calculation are written to the "chnpot.log" file. The program writes the new natural potentials file "natpot" and the "dvr" file (if the option -rd or keyword readdvr in Run-Section are not used) to the name directory. For over-writing of the natpot file use "-w" option. For an example input file see chnpot.inp

Example: chnpot84 -w -rd input

Note that chnpot may be used to change the modelabels. Chnpot maps the n-th old DOF onto the n-th new DOF, the modelabels are not used. If you just want to change the modelabels without re-fitting the potential, simply set up the PRIMITIVE-BASIS-SECTION of the chnpot input file using the new modelabels but otherwise the old basis definitions. A FIT-SECTION is not needed in this case.

comparever

Is used by the elk_test program to compare different sets of results.

concatver

Purpose: Build the (Hartree) product of two MCTDH wavefunctions.
Usage : concat84 [-ver -D -dvr -w -h -?] psi1 psi2

Options :
-ver    : Version information about the program
-D DIR  : Data are stored in directory DIR,
rather than in the current directory.
The directory DIR will be made, if not existing.
-dvr    : The dvr-file for the new WF will be build.
-w      : Existing files are overwritten.
-h      : Print this help text.
-?      : Print this help text.
The arguments psi1, psi2 are the names of the restart-files to be
to be concatenated, i.e. psi = psi1*psi2 .
If no argument is given, the user is prompted for the missing argument.
The first WF may have electronic states  (in single- or multi-set),
but the second WF must not have electronic states.
Packets are not allowed.
------------------------------------------------------------------------------

The string restart will be appended to the path if necessary, I.e. concat84 dir1 dir2 is equivalent to concat84 dir1/restart dir2/restart.
The log file cc.log contains useful information. The option -dvr is useful if showsys84 -rst is to be run.

convpsiver dir p

Converts a multi-packet wavefunction to the corresponding single-packet ones, i.e. extracts the pth packet from the multi-packet wavefunction dir/psi and writes it to the file dir/psip (or dir/psi0p if p < 10). For p = 0 all packets contained in psi are extracted and written to the files dir/psi01, dir/psi02, etc.

crosscorrver

Purpose: Calculates the cross-correlation function c(t)
c(t) = <psi_ref(0)|psi(t)>

Usage  : crosscorr84 [-f -i -o -n -r -R -O -t -nofs -el -ig
-spf -diag -ver -h -?]  psi_ref .

Options :
-f FILE : The wavefunction is read from file FILE
rather than from ./psi
The string FILE may be a relative or a full path-name.
-o FILE : The difference is written to file FILE rather than to ./cross_XXX
The string FILE may be a relative or a full path-name.
If FILE=no, i.e. "-o no", no cross-file will be opened
and the output will be directed to screen.
-r      : Take the restart rather than the psi file as reference WF.
-R      : Take the restart rather than the psi file as time-dep. WF.
In general, this makes sense only when also -r is set.
-O oper : Apply operator "oper" before evaluating the overlap.
-t tfin : Compute the overlap only up to tfin.
-nofs   : No transformation to fs. Use when 'time-not-fs' was set in mctdh.
-n step : Compute the overlap only every step-th output.
-i init : Compute the overlap using the init-th (rather than the first)
reference WF. (-i is ignored when -r is set).
-el nst : Compute the cross-correlation for state nst only.
-ig     : Ignore "different bases" error-msg.
-spf    : Print information on overlap of the two orbital sets.
dim - trace(dover^+ * dover) for all m,s.
-diag   : Print diagonal elements of overlap matrix of the two orbital sets.
-ver    : Version information about the program
-h      : Print this help text.
-?      : Print this help text.

The argument psi_ref is the name of the file containing the
(time-independent) comparison wavefunction.

With the aid of the option -O, crosscorr can compute matrix elements
<psi_ref|oper|psi>. This, however, requires that psi_ref and psi
are MCTDH-WFs with identical numbers of SPFs.

If no argument is given, the user is prompted for the missing argument.
Note: The option "-o no" directs output to screen.
Note: If the argument ends with a slash "/", the string "psi" or "restart"
is appended. If the argument ends with "restart", -r is set.
-----------------------------------------------------------------------------
The two wavefunctions must, of course, have the same primitive grids, but they may have different combination schemes and may differ in numbers of SPF's. A reference wavefunction can conveniently be generated by running MCTDH with the -t option.

crosspecver

Purpose: Computation of the spectrum by fourier-transforming
the cross-correlation function.
Output is: Energy, Re(FT), Im(FT), Abs(FT),
where FT denotes the Fourier-transform (prefactor Pi^-1)
of the cross-correlation function.

Usage : crosspec84 [-o -e -p -q -a -g -all -reim -re -im -sin -cos -dc -t
-nofs -G -ph -r -old -h ] [emin emax unit file].

Options :
-o FILE : The spectrum is written to file FILE.pl
rather than to ./cspec.pl
The string FILE may be a relative or full path-name.
-e R Unit : the offset energy is set to R in units Unit; e.g. 0.3 ev
-p I    : The no. of plot-points is I rather than  500
Use "-p 2000" for final plots
-q I    : Same as -p I
-g      : GNUplot command lines are written to the output file.
-a      : GNUplot is called automatically. (sets "-g").
-all    : Re, Im and Abs parts (rather than Abs) are plotted.
(ignored if neither -g or -a are set).
-reim   : Real and Imaginary parts (rather than Abs) are plotted.
(ignored if neither -g or -a are set).
-re     : Real part (rather than Abs) is plotted.
(ignored if neither -g or -a are set).
-im     : Imag. part (rather than Abs) is plotted.
(ignored if neither -g or -a are set).
-sin    : The crosscorrelation function is multiplied with sin(pi*t/T).
-cos    : The crosscorrelation function is multiplied with cos(pi*t/2T).
-dc     : A DC component is subtracted from the cross-corr before FT.
-r      : Plot the resolution rather than the spectrum
-ph  ph : Use phase correction exp(i*ph*1.d-6*t^2)
-t tcut : The crosscorrelation function is read only up to t=tcut
-nofs   : No transformation to fs. Use when 'time-not-fs' was set in mctdh.
-G      : Plot gridlines (meaningfull only if -g or -a is set).
-old    : The old output format is used.
(Different order of colums, output times Pi.)
-ver    : Version information about the program
-h      : Print this help text.
-?      : Print this help text.
Following the options one may provide the arguments (separated by blank):
Emin, Emax, Unit file
where file denotes the path of the crosscorr file.
If not enough arguments are given, the user is prompted for the missing
arguments. Use the Unit 'no' if the propagation is run with the
keywords: 'time-not-fs','energy-not-ev'.
------------------------------------------------------------------------------
This routine simply performs the Fourier integral pi^-1 * Int[0,T] exp(i*(E-E_offset)*t)*c(t) dt without invoking damping functions (except possibly for sin(pi*t/T) and cos(pi*t/2T)). The option -dc subtracts the average of the cross-correlation function from this function before FT. E_offset is zero, if not set by option -e. The option -old removes the prefactor pi^-1 and interchanges the columns. (Now: Energy, Re, Im, Abs; with -old: Energy, Abs, Re, Im). The new output format is chosen to make crosspec compatible with autospec. The output file cspec.pl may be read by flux, using the -ed option of flux84. In this case neither -g nor -a must be given, such that no GNUplot command lines appear in the output file.

If the option -r is set, the program computes the resolution, i.e. the spectrum of a single line. The time points are taken from the crosscorr file. If there is no crosscorr file and if the option -t is given, a time step of 1 fs is assumed.

dengenver

Purpose: Generates 2D-densities from the wavefunction.
Usage:   dengen84 [-i -f -n -skip -step -o -rst -pop2 -nw -ver -h -?] f1 f2

Options :
-i DIR  : data is stored in directory DIR
-f FILE : the wavefunction is read from FILE rather than from ./pes
The string FILE may be a relative or a full path-name.
-o FILE : The output is written to file FILE rather than to ./dens2d
The string FILE may be a relative or a full path-name.
-n num  : only num WFs are processed.
-skip m : the first m WFs are skipped.
-step step : only every step's WF will be processed.
-rst    : the restart file is read,  rather than the psi file.
-nw     : No weights are employed. Plot the "naked" DVR populations.
-pop2   : the basis set occupations are plotted rather than grid
populations. For FFT this implies momentum space representation.
-w      : Overwrite the output file.
-ver    : Version information about the program.
-h      : Print this help text.
-?      : Print this help text.

The arguments f1 and f2 denote the DOFs for which a 2D-density is generated.
The arguments f1 and f2 may be integers or modelabels.

If there are more WFs on psi file than you want to inspect, use the
options -n, -skip, or -step. Example: "dengen84 -skip 10 -n 2  rd rv" will
generate densities of the WFs no. 11 and 12 (i.e. t=10 and 11 if tpsi=1fs),
and "dengen84 -step 10 -n 2  rd rv" will generate densities of the
WFs no. 1 and 11 (i.e. t=0 and 10 if tpsi=1fs).
NB: The option -pop2 is ignored for phifbr,sphfbr,k,external bases.
--------------------------------------------------------------------------------
dengen84 computes 2D-densities similar as does showsys84. However, in contrast to showsys84, dengen84 computes only 2D-densities, is not menu driven, is not interactive and does not plot. It merely writes the data to an output file, which may then visualized by GNUPLOT. If the computation of the 2D-densities takes a large amount of computer time, the use of dengen84 is of advantage, because dengen84 may run in the background, overnight or on a batch system.
dengen84 computes the diabatic densities for multistate wavefunctions. The output file dens2d.dat displayes the coordinate values in the first two columns and then the densities, one column for each electronic state.

Purpose: Converts a diabatic wavefunction to adiabatic.
Usage:   di2ad84 [-i -f -p -o -w -ver -h -?]

Options :
-i DIR  : data is stored in directory DIR
-f FILE : the wavefunction is read from FILE rather than from ./pes
The string FILE may be a relative or a full path-name.
-p FILE : the PES is read from FILE rather than from ./pes
The string FILE may be a relative or a full path-name.
-o FILE : The ouput is written to file FILE rather than to ./psi.ad
The string FILE may be a relative or a full path-name.
-w      : An existing output file is overwritten.
-ver    : Version information about the program
-h      : Print this help text.
-?      : Print this help text.

The wavefunction and the pes-file must be in "exact" format.
One may use psi2ex84 to transform the WF.
--------------------------------------------------------------
This program reads the wavefunction in a psi file, and uses the PES information in a pes file to transform it to the adiabatic representation. The waveunction and pes file must be in a full-grid representation. If the psi file was created using an MCTDH wavefunction, this can first be converted using the psi2ex program. The pes file can be created using mctdh84 -pes with an input file for a numerically exact calculation.

Output is in a full grid representation to the file psi.ad, unless otherwise given using -o FILE.

diagver

Program solves the generalized eigenvalue problem H1 B = H0 B E by a singular value decomposition, where H0 and H1 are overlap and Hamiltonian matrices calculated by Fmat program. The output of the program are the eigenvalues and intensities. The matrix elements are calculated using the hermitian or symmetric scalar product. The widths, in the case of symmetric scalar product, are calculated as -2*Im(E) .

Purpose : Diagonalizes the Hamiltonian matrix calculated by fmat.
Usage   : diag84 [-d -o -lo -hi -n -t -es -e -u -w -ver -h -?] mtt_ovl
mtt_ham [[mtt_ovl_2] [mtt_ham_2]].

Options  :
-d PATH  : The output is directed to PATH rather than to the current directory.
The string PATH may be a relative or a full path-name.
If PATH=no, i.e. "-d no", no output files will be opened.
-o NAME  : Change the name of output file from "diag" to NAME.
-lo klo  : Take matrix elements in Mtt files only from the klo-th index.
-hi khi  : Take matrix elements in Mtt files only up to the khi-th index.
-n step  : Take only every n-th step of the Mtt matrix elements.
-t       : Expand the nXn mtt-matrix to (2n-1)X(2n-1) one including negative
times. In this case the second set of overlap and Hamiltonian
matrices are needed.
-es eps  : Set cut-off parameter for eigenvalues. Default is eps=1.e-10
-e E unit: The offset energy is set to E unit.
Note: "-e lsth" is equivalent to "-e 4.74746 ev".
"-e co2" is equivalent to "-e -2534.5298 cm-1".
-u UNIT  : Set the output energy unit to UNIT. Default is au.
-w       : Overwrite an existing "diag" file.
-ver     : Version information about the program
-h       : Print this help text.
-?       : Print this help text.
The arguments "mtt_ovl" and "mtt_ham" are the Overlap- and Hamiltonian
matrix files calculated by  "fmat84" program. If "-t" option is used,
the second set of Overlap- and  Hamiltonian matrices with different
scalar product is needed.

Examples  : diag84 -es 1.d-8 -hi 60  mtt_h mtt_H1_h
: diag84 -hi 40 -n 2 -t mtt_h mtt_H1_h  mtt_s mtt_H1_s
Important note: Both overlap and Hamiltonian-matrices should be calculated
either as hermitian or as symmetric scalar products.
------------------------------------------------------------------------------

Notes:

The arguments with options -lo,-hi and -n i.e. klo, hki and step are related to the calculated "mtt"-files, and not to the MCTDH calculation. For example, if "mtt"-files were calculated with step=2, then using step=3 in Diag program will be equivalent to step=6 with respect to the MCTDH calculation.

If "-t" option is used, the matrix elements are extended to negative times according to the rule: Psi(-t)=conj(Psi(t)). So, here the second set of overlap- and Hamiltonian matrices needed, calculated using the symmetric scalar product, if the first is hermitian, or hermitian scalar product if the first is symmetric. Note, that each set of overlap and Hamiltonian matrices should be calculated using the same scalar product (hermitian or symmetric).

The option "-u UNIT" changes the output energy unit to UNIT. The available units in program can be seen by shell-command "mhelp -s unit".

The following examples show the steps performing the diagonalization after completing the MCTDH calculation. The first concerns the reactive system (to calculate resonances and widths), the second - the bound system.

1. fmat84 -sym          ! calculates  mtt_s
fmat84 -sym -O H1    ! calculates  mtt_H1_s
diag84 -u ev -es 1.d-9 mtt_s mtt_H1_s

2. fmat84               ! calculates  mtt_h
fmat84 -O H1         ! calculates  mtt_H1_h
fmat84 -sym          ! calculates  mtt_s
fmat84 -sym -O H1    ! calculates  mtt_H1_s
diag84 -u cm-1 -es 1.d-8 -n 2 -t mtt_h mtt_H1_h mtt_s mtt_H1_s

edengenver

Purpose: Generates pair-density matrix  from the wavefunction.
Usage:   edengen84 [-i -f -n -skip -step -o -rst -pop2 -nw -w -ver -h -?] f1 f2

Options :
-i DIR  : data is stored in directory DIR
-f FILE : the wavefunction is read from FILE rather than from ./pes
The string FILE may be a relative or a full path-name.
-o FILE : The ouput is written to file FILE rather than to ./pair2d
The string FILE may be a relative or a full path-name.
-n num  : only num WFs are processed.
-skip m : the first m WFs are skipped.
-step step : only every step's WF will be processed.
-rst    : the restart file is read,  rather than the psi file.
-nw     : No weights are employed. Plot the "naked" DVR populations.
-pop2   : the basis set occupations are plotted rather than grid
populations. For FFT this implies momentum space representation.
-w      : Overwrite the output file.
-ver    : Version information about the program.
-h      : Print this help text.
-?      : Print this help text.

The arguments f1 and f2 denote the DOFs for which a pair-density matrix is
generated. The arguments f1 and f2 may be integers or modelabels.

If there are more WFs on psi file than you want to inspect, use the
options -n, -skip, or -step. Example: "edengen84 -skip 10 -n 2 rd rv" will
generate densities of the WFs no. 11 and 12 (i.e. t=10 and 11 if tpsi=1fs),
and "edengen84 -step 10 -n 2 rd rv" will gen. dens. of the WFs no. 1 and 11.
NB: The option -pop2 is ignored for phifbr,sphfbr,k,external bases.
--------------------------------------------------------------------------------
edengen84 computes 2D (off-diagonal) pair-correlation functions (or pair-density matrices) <x,x|ρ2|y,y> (similar to the one-body density matrix produced by prhosub84). It works in analogy to dengen84, which however calculates the diagonal density <x,y|ρ2|x,y>. The data is written to an output file, which may then visualized via GNUPLOT (e.g with plgen). Note that, if the computation of the pair-correlation function is very time consumptive, edengen84 may run in the background, overnight or on a batch system.
edengen84 computes the diabatic densities for multistate wavefunctions. The output file dens2d.dat displayes the coordinate values in the first two columns and then the densities, one column for each electronic state.

expectver

Purpose: Calculates the time evolution of an expectation value.
Usage  : expect [-f -i -o -nofs -g -a -G -r -rst -x -t -p -pt -u -n
-skip -w -ver -h -?] [operator].

Options :
-f FILE : The operator is read from file FILE rather than from ./oper
The string FILE may be a relative or a full path-name.
-i FILE : The wavefunction is read from file FILE rather than from ./psi
The string FILE may be a relative or a full path-name.
-o FILE : The output is written to file FILE.pl rather than to ./expect.pl
The string FILE may be a relative or a full path-name.
If FILE=no, i.e. "-o no", no output file will be opened,
and the output is directed to screen.
-nofs   : No transformation to fs. Use when 'time-not-fs' was set in mctdh.
-g      : GNUplot comand lines are written to the output file(s).
-a      : GNUplot is called automatically. (sets -g and -w).
-G      : Plot gridlines (meaningfull only if -g or -a is set).
-p      : Print GNUplot to lpr (sets -g).
-pt thr : Use thr pthreads for calculation.
-r      : Do not divide the expectation values by norm^2.
-rst    : Use restart rather than psi file.
-x xi xf: Set the abscissa length.
-t tfin : The expectation values are computed only up to time=tfin.
-u UNIT : The unit used is UNITs (Default is a.u.).
-n step : Compute the property only every step-th output.
-skip nskip : Skip the first nskip wavefunctions.
-w      : An existing property file is overwritten.
-ver    : Version information about the program.
-h      : Print this help text.
-?      : Print this help text.

The argument operator is the name of the operator, as specified in an
HAMILTONIAN-SECTION_xxx, for which the expectation value is required.
If no argument is given, the user is prompted for the missing argument.
------------------------------------------------------------------------------
The operator must first have been build by the MCTDH program. This means that it must first be set up in a .op file. If it is known before a propagation is made, this operator can be included with the system operator in the main oper file (by using a named HAMILTONIAN-SECTION_XXX). Alternatively, e.g. if a propagation has already been made, the file oper_name can be generated using the genoper=name keyword.

Start the program in the directory containing the propagated wavefunction in a psi file. If no arguments are given, the program looks in the oper file to see which operators are there and asks which is to be evaluated. The -o option can use a file other than this. The expectation value of this operator is then calculated for each time, and output to the file expect.pl, unless otherwise instructed. The -a option enables the automatic plotting of the results.

fdcheckver

Purpose: Checks whether a state has been found in all comparison
calculations  and computes the internal error.

Usage  : fdcheck84 k exact datasets

k = number of datasets.
exact = 0 if if there is no exact comaprison data set, =1 else.
datasets = data as created by fdmatch84.
Example: fdcheck84 4 0 < list
where list contains the output of fdmatch84, or:
fdmatch84 f1.eig f2.eig f3.eig | sort -n | fdcheck84 3 0
---------------------------------------------------------------------------

fdcheck analyses the output of fdmatch. It computes the mean value and the variance of eigenvalues and intensities which are grouped together by fdmatch. If there is an exact data set, fdcheck computes the errors of eigenvalues and intensities with respect to the exact data.

fdmatchver

Purpose: Finds the optimal matches between a number of lists of
eigenvalues obtained with the (real) filter program.

Usage  : fdmatch84 file_1 ... file_k | sort -n | less

where file_j (j=1...k) is typically a *.eig file generated by the
filter program. There may be up to 10 files.

-----------------------------------------------------------------------

fdmatch orders the eigenvalues computed by different filter runs such that optimal matches are obtained. The output of fdmatch may be piped to fdcheck.

fluxver

Purpose: Computes the quantum flux.
See review Eqs.(197-201) for a definition of the output

Usage : flux84 [-f -i -o -ed -dvr -t -n -lo -hi -p -P -O -e -x -d -v -w -r -c -exp -u
-sm -es -wtt -nofs -noev -ver -h -?] [emin emax unit CAP-dof CAP-type] [.inp-file].

Options :
-f FILE : The wavefunction is read from file FILE rather than from ./psi
The string FILE may be a relative or a full path-name.
-i FILE : The operator is read from file FILE rather than from ./oper
The string FILE may be a relative or a full path-name.
-dvr FL : The DVR information is read from file FL rather than from ./dvr
The string FL may be a relative or a full path-name.
-ed FL  : The energy distribution is read from file FL rather than from ./enerd
The string FL may be a relative or a full path-name.
One may use a flux-file as energy-distribution file.
In this case, the flux-column will be taken as edstr.
If a "d" follows the file-name, the DELTA(E) column will be taken.
If a "+" follows the file-name, the present flux will
be added to the edstr input.
Alternatively, one may use a spectrum.pl (autospec84)
or a cspec.pl (crosspec84) file as energy distribution.
The options "d" and "+" work here as well. With "d"
the third (rather than the second) column of the ed-file
is used, and with "+" the data is added to DELTA(E)
rather than overwriting it.
-o PATH : The output is directed to PATH rather than to the current directory.
The string PATH may be a relative or a full path-name.
If PATH=no, i.e. "-o no", no output-files will be opened.
-n step : Compute Wtt' only every step-th psi-output.
-es nst : Compute the flux for the electronic state nst only.
-lo klo : Compute Wtt' only from the klo-th psi-output onwards.
-hi khi : Compute Wtt' only up to the khi-th psi-output.
-p I    : The no. of energy-points is I rather than  500.
-e R U  : The offset energy is set to R in units U; e.g. "-e 0.3 ev".
Note: "-e lsth" is equivalent to "-e 4.74746 ev".
Give here the negative of "Energy shift (E_total - E_kinetic)",
which was printed to log-file. (correction=dia).
-s strt : Set the starting point of the flux evaluation to strt.
Default is the grid-point where the cap starts.
-t      : Test run.  Check of the data sets. No output files opened.
-P      : Use projectors. -P is follwed by the mode number, the projector
name and up to 8 parameters. The input is closed by a "%".
There may be several projectors. Example:
-P 3 KLeg 2 1 % -P 2 eigenf vib 6 %
-O oper : Apply operator "oper" before measuring the flux.
-u unit : Transform flux to unit "unit".
(useful when a operator is applied).
Note: Rather than a unit one may give a number.
-w      : The gtau file may be overwritten. Without -w option
gtau is read from the gtau file (if it exist).
-r      : Enforce reading the gtau file. Stop if gtau does not exist.
-v      : Re-compute only the theta part of gtau.
(useful in combination with the -s option).
-x      : Evaluate flux only for the theta part of gtau.
(useful for tests in combination with the -s -v options).
-d      : Compute gtau from derivative of wavefunction.
-c icos : Cosine**icos damping of gtau.
-exp n tau: exp(-(t/tau)**n) damping of gtau
-sm     : smooth Delta(E) by averaging the enerd data over 5 points.
-Mb norm: Flux spectrum in Mb. norm is ||D*Psi|| (in au)
where D is the Dipole moment.
-wtt    : Only Wtt is computed. Note: The arguments
"Emin, Emax, unit" must not be given when -wtt is set.
-nofs   : No transformation to fs.
Use when "time-not-fs" was set in mctdh.
-noev   : No transformation to eV.
Use when "energy-not-ev" was set in mctdh.
-ver    : Version information about the program
-h      : Print this help text.
-?      : Print this help text.

The arguments specify the energy interval and
the degree of freedom of the CAP.
The CAP-dof may be a number or a modelabel.
The CAP-type may "right", "left" or "both". "right" is the default.
If the number given as CAP-dof is larger than the number of DoFs,
this number is interpreted as the H-term, h, of the CAP.
Example: flux84 -lo 25 -w 0.3 2.5 ev rd
flux84 1000 3000 cm-1  x left
flux84 -wtt -s 21  y
If no argument is given, the user is prompted for the missing argument(s).
The path of an input file may be given as the very last argument.
Options will overwrite the input-file.
------------------------------------------------------------------------------

The flux84 analyse routine has currently two restrictions:

• When projectors are used, they must project fully on a MCTDH-particle, i.e. projection on a dof within a combined mode is not possible.
• Projectors or operators operating on the CAP-mode are not allowed. When projectors and operators are used simultaneously, then they must operate on different modes.

flux84 can work with operators and/or projectors, however, they must not operate on the CAP-mode. The operators are to be defined in an HAMILTONIAN-SECTION_xxx, where xxx is the name of the operator to which flux84 refers. Presently there are projectors onto eigenfunctions of one-mode operators (also to be defined in an HAMILTONIAN-SECTION_xxx ) and projectors onto diffractive states (plane-waves) and onto rotational states. For a full list of projectors available, see the table below. New projectors may be coded and added to the file analyse/flxprj.F. A very general and quite useful projector is the read projector. This projector is build from an SPF, P=|psi><psi|, which in turn is read from a (foreign) restart file. The arguments needed by flux84 are the energy-range ( including its unit) and the CAP-mode, i.e.: Emin, Emax, unit, CAP-mode . Often options are necessary to steer the performance of flux84. Rather than options and arguments one may provide an input-file (for an example see below). Note: options will overwrite the statements of the input-file.
The option -t (Test run) is very useful for checking how many data sets are on the psi file, which CAP's have been used, what are the mode labels, etc. The most time consuming part is the generation of the function gtau, as it requires the computation of the matrix elements Wtt' of the CAP. If a gtau file exist, flux84 will try to read this file (provided the -w option is not given). Thus the computation of the flux for an other energy interval or a different flux-unit (-u) is fast, as the gtau file already exist. When the gtau file is corrupted, use -w to overwrite it. The options -n, -hi and -lo may be used to reduce the number of time-points when computing Wtt' and gtau. This is useful for convergence checks. The option -lo must be used when the mctdh calculation was run with auto-cap (ACAP).
The option -s start determines the grid-point from where on the flux is evaluated. To be more precise, start is the grid-point at which the g_theta evaluation is switched on. If the option -s is not given, start is set to the grid-point where the CAP is switched on. When start is larger than gdim(fcap) (i.e. grid length of the CAP-DOF), then the evaluation of the g_theta part is suppressed. See Appendix B of J.Chem.Phys 116 (2002), 10641, or the MCTDH-feature article (Theor. Chem. Acc. 109, 251 (2003)) for more details. It often improves the speed of convergence, in particular with respect to the total propagation time, if the start of the flux evaluation is moved closer towards the scattering center. Only for testing purposes one may chose a start value, which is further away from the scattering center as the default value, i.e. a start value which lies inside the CAP region.
The CAP may be left (i.e. located at small coordinate values, or right (i.e. located at large coordinate values, this is the default), or both (i.e. the CAP may be anywhere, no checking of areas where the CAP vanishes).

When the option -Mb is given, the flux spectrum will be multiplied with E*norm**2 where norm is the argument to -Mb. The value of norm is given in the log-file of the mctdh run where operate was performed (Operate-Norm). norm is equal to ||D*Psi||, where D denotes the operator used by operate. Make sure that a correct energy shift (-e option) is given, as the spectrum is multiplied with E.

Lines like the following are output by flux84 to screen (see also the flux.log file)

------------------------------------------------------------------------------
2000*Integral(W_tt)  =   996.1052     (total outgoing  flux)
1000*Integral(flux)  =   995.6810     (total outgoing  flux)
1000*Integral(DELTA) =   998.8516     (total incoming  flux)
------------------------------------------------------------------------------

The first and second lines gives the total flux going into the specified CAP. The total flux is obtained by (first line) analytically integrating the flux over all energies (see Eq.(118) of Theor. Chem. Acc. 109, 251 (2003)), or by numerically integrating the energy resolved flux over the specified energy interval. The third line is the numerical integral of DELTA(E) over the specified energy interval. Here:
DELTA(E) = < Psi | delta(H-E) | Psi >
flux(E) = 2*pi*< Psi |delta(H-E) F delta(H-E) | Psi >
where F denotes the flux operator. DELTA(E) can be obtained trough several ways. mctdh84 may generate the so called enerd file (by giving the keyword correction=dia in the Init_WF-Section). (NB: The file enerd is called adwkb in older versions). flux84 will automatically look for such a file. Alternatively (by using the -ed option) one may use a spectrum.pl file (generated by autospec84), or a cspec.pl file (generated by crosspec84), or a flux file (generated by flux84 without employing projectors) as energy-distribution file.

Example (H+D2 reactive scattering) :
flux84 -e lsth -lo 7 0.25 2.75 ev rv
or
flux84 flx.inp

where flx.inp is an input file. The extension .inp may be dropped. The input-file is organised similar to the MCTDH input-file, however, there are no sections. The last line of the input-file must read: end-input. Below an example input file is shown.

-------------------------------------------------------------
# flux input file. Surface scattering.
# Average rotational energy for a specific diffraction state.

outdir = flux_rot  low = 26   flux-unit = mev
operator = rot
projector 3 diff2  2  1   #  mode ,name, parameters
energy = 0.05, 0.30, ev    cap = z
end-input
-------------------------------------------------------------
Using the -d option (or the keyword use-deriv) the flux is computed via a well known formula employing the derivative of the wavefunction. Matrix elements of the CAP are not evaluated. Presently this feature works only if there is a full projection and if the CAP degree of freedom is uncombined and represented by FFT. When using the derivative formulation it is recommended to move the evaluation point a few grid-points before the start of the CAP. (-s option).

Using the -wtt option one can disable the the Fourier-transform to energy space and investigate the time-dependence of the flux. This information is stored in the wtt and iwtt files. The iwtt file (integral wtt) contains the values Int(0,t) 2*<psi(t')|W|psi(t')> dt' + <psi(t)|Theta|psi(t)> where Theta denote a Heaviside function starting at the grid-point given by the -s option. (If -s is not given, the starting point of the CAP is used). Hence iwtt provides the probability that the WF has passed the point given by the -s option.

Example (Investigation of particle loss):
flux84 -wtt -s 21 rv

Keywords of flux input-file
Keyword   Description equivalent option
overwrite overwrites gtau file. gtau will be re-computed completely. -w
test performs a test run -t
psi = S read psi file form path S -f
outdir = S write output to directory S -o
dvr = S read dvr file form path S -dvr
edstr = S read the energy-distribution from file S -ed
flux-unit= S The computed flux is converted to unit S. If a number is given, the flux will be multiplied with that number. I.e. flux-unit=ev and flux-unit=27.2114 will produce the same output. This option is useful when an operator is applied. -u
operfile = S read oper file form path S -i
high = I Compute the flux (i.e. gtau) only up to the I-th psi-output -hi
low = I Compute the flux (i.e. gtau) only from the I-th psi-output onwards. -lo
el-state = I Compute the flux for the I-th electronic state only. -es
step = I Compute the flux (i.e. gtau) only every I-th psi-output. -n
epoints = I The number of output energy-points is I rather than 500. -p
cos = I Cosine**I damping of gtau -c
exp = I, R exp(-(t/R)**I) damping of gtau -exp
Mb = R The argument is the value of Operate-Norm to be taken from the log-file of the mctdh operate run. The flux-spectrum is multiplied by E*opnorm**2. Make sure that a correct energy shift is given. -Mb
start = I Evaluate the flux at the grid-point start. The default value for start is the start of the CAP (i.e. the last grid-point with vanishing CAP). -s
theta-only Uses only the theta-part of gtau when computing the flux. This is useful for investigating the importance of the theta-region. Note: The W-part only of gtau will be taken when start is set to a number larger than the number of grid-points of the CAP-mode. -x
gtau-only Compute only the theta-part of gtau. This option requires that a gtau file exist. The time-consuming evaluation of the W-part of gtau is suppressed, but gtau_theta is re-evaluated. This option is useful to test different start points (see keyword start). Note: the keywords low and high will be ignored. There values are taken from the previous run. -v
wtt-only Compute only the diagonal terms Wtt.This rather fast calculation is useful to determine the optimal value for low (-lo). -wtt
use-deriv Compute gtau from spatial derivative of wavefunction. -d
eoffset = R, S The offset energy is set to R in units S; e.g. "eoffset = 0.3, ev". Note: "eoffset = lsth" is equivalent to "-e 4.74746 ev". -e
operator = S Apply operator S before evaluating the flux. -O
projector I S (S1) I1 I2 .. Apply the projector S with(label S1 and) parameters I1 I2 ... to mode I before evaluating the flux. projector has to be the first keyword on a line. There are no equal-sign and no commas. There may be several projectors. One line for each projector. No additional keywords in a projector-line. -P
nofs No transformation to fs. Use when "time-not-fs" was set in mctdh -nofs
noev No transformation to eV. Use when "energy-not-ev" was set in mctdh -noev
energy = R1, R2, S Evaluate the flux between Emin=R1 and Emax=R2 where Emin and Emax are given in unit S. (arguments)
cap = S or I The cap is on mode I. Rather than giving the number of the mode, I, one may give the modelabel S. If I is larger than the number of DOFs, this number is interpreted as the H-term, h, of the CAP (run "flux84 -t" and/or see op.log file). This feature becomes important if there are several CAPs. (arguments)
captype = S The cap-type is S=both,right,left. (arguments)

Projectors

Projectors
projector name label parameters Description
eigenf name of an operator   pop An eigenfunction of a 1D operator is taken as projector. The 1D operator is to be defined in a HAMILTONIAN-SECTION_ XX when the operator-file is build. The operator XX must be of usediag type. A 1D spf is assumed, i.e a combined mode is incompatible with eigenf. The parameter pop defines which eigenfunction is taken. pop=1 refers to the ground-state.
read path of restart directory   pop   state The pop-th SPF of a "foreign" restart file is taken as projector. This "foreign" restart file may be build by an geninwf run. The primitive grids and the combination scheme of the "foreign" restart file and the present calculation must be identical, but the numbers of SPF's and the number of electronic states may differ. If the parameters pop and state are not given, they are set to 1.
sph (none)   j    m   projects on the (j,mj) rotational state. A sphFBR (see Primitive Basis Section) is required, the two modes of which must not be further combined. I.e a 2D spf is assumed.
diff (none)   n   projects on the n-th diffraction channel. (i.e. onto a plane wave). A 1D spf is assumed.
diff2 (none)   n   m   projects on the (n,m)-th diffraction channel. (i.e. onto a 2D plane wave). A 2D spf is assumed.
pop2 (none)   n   projects on the n-th basis function of an underlying DVR. This projector thus can only be used if the (uncombined) DOF is described by an simple DVR.
Leg (none)   j   projects on the j-th rotational state. A Leg or Leg/R DVR is required as primitive basis and the DOF should be uncombined.
KLeg (none)   j    m   projects on the (j,mj) rotational state. A KLeg-DVR combined with a K-DVR, i.e. a 2D spf, is required.
PLeg (none)   j    m   projects on the (j,mj) rotational state. A PLeg-DVR combined with an exp-DVR, i.e. a 2D spf, is required.
KLeg2 (none)  j1 m1 j2 m2 sym projects on the (j1,m1;j2,m2) rotational state of two coupled angular DOFs. The mode must be of the type KLeg/K/KLeg/K, i.e. a 4D mode. If sym>0 (sym<0), the projection is on the symmetrized (antisymmetrized) state, i.e.
{|j1,m1,j2,m2> ± |j2,m2,j1,m1>} / √2  .
Wigner (none)   j   k   m   projects on the (j,k,m) rotational state. A Wigner-DVR combined with two DVRs of either type exp or K, i.e. a 3D spf, is required.

fmatver

Purpose: Computes the overlap- or the Hamiltonian-matrix.
Usage  : fmat84 [-f -i -d -o -dvr -n -lo -hi -O -w -sym -es -ver -h -?] .

Options :
-f FILE : The wavefunction is read from file FILE rather than from ./psi
The string FILE may be a relative or a full path-name.
-i FILE : The operator is read from file FILE rather than from ./oper
The string FILE may be a relative or a full path-name.
-dvr FL : The DVR information is read from file FL rather than from ./dvr
The string FL may be a relative or a full path-name.
-d PATH : The output is directed to PATH directory rather than to the
current directory.
The string PATH may be a relative or a full path-name.
If PATH=no, i.e. "-d no", no output-files will be opened.
-o NAME : Change the name of the output file from "mtt_h"("mtt_s") to NAME.
-n step : Compute Mtt' only every step-th tpsi.
-es nst : Compute Mtt' for the electronic state nst only.
-lo klo : Compute Mtt' only from the klo-th tpsi onwards.
-hi khi : Compute Mtt' only up to the khi-th tpsi.
-O oper : Apply operator "oper" before evaluating the overlap.
-w      : Overwrite an existing Mtt file.
-sym    : Use the symmetric scalar product rather than the hermitian one.
-ver    : Version information about the program
-h      : Print this help text.
-?      : Print this help text.

Example : fmat84 -hi 65 -O ham
Important note: The used operator must be build with the "nodiag" flag!
------------------------------------------------------------------------------

N.B. The outputfile is called  "mtt_h" or "mtt_s" if no operator is used (i.e. if the overlap matrix is calculated) and is called   "mtt_h_oper" or "mtt_s_oper" if the operator oper is employed. The flags "h" and "s" are for the case of hermitian and symmetric scalar products. The operator oper must be defined in a HAMILTONIAN-SECTION_oper and must carry the nodiag flag. The generated matrices may be diagonalised (generalised eigenvalue problem) to obtain the eigenvalues of the operator oper. This should be done with the program diag. Resonances may be calculated as well by using the symmetric scalar product (-sym option).

jinpolver

Purpose : J-interpolation of flux data.
Usage   : jinpol84 [-d -f -o -j -j1 -j2 -m -r0 -rf -e1 -e2 -ed -es
-em -mp -A -B -C -p -ver -h].

Options :
-id IDIR: The probabilities files are read from directory IDIR.
-od ODIR: The probabilities files are written to directory ODIR.
-f FILE : The probabilities  are read from file IDIR/FILEext
where ext denote the extension jKJ (e.g. 205).
-o FILE : The interpolated probabilities are written to ODIR/FILEext
where ext denote the extension jKJ (e.g. 205).
-j   j  : Search for files with j-value j. j is the first number in extension.
-j1  J1 : First J value of interpolation range.
-j2  J2 : Last  J value of interpolation range.
-k   K  : Search for files with K-value K. K is the second number in extension.
-r0  R0 : R0 is the estimated value of the transition-state coordinate.
-rf  Rf : Rf is the fixed value of the transition-state coordinate.
(No optimization).
-m mass : Reduced mass in AMU (default mass = 1.0).
-e1  E1 : First energy of fit-interval.
-e2  E2 : Last  energy of fit-interval.
-ed edim: Number of energy points used to compute the integrals. (def.= 500).
-es shft: The energies read from the probability files are shifted by shft.
-em Emin: The input probabilities are set to zero for E<Emin.
-mp maxp: The input probabilities are set to maxp if prob>maxp.
-A      : Minimize F = SQRT(A*B) - C
-B      : Minimize F/C (default)
-C      : Minimize -C
-p      : Print every iteration step.
-ver    : Version information about the program
-h      : Print this help text.

The defaults settings are: IDIR=., ODIR=., j=0, k=0, r0=3.0, Edim=500,
Eshift=0.0, Emin=-1.d9, maxp=1.d9, mass=1.0 (i.e. 1837.15 au).

Example: jinpol84 -f flux -o mflux -j 2 -j1 1 -j2 4 -e1 0.0 -e2 1.2
This will read the files flux201 and flux204 from and write the
files  mflux201, mflux202, mflux203, and mflux204 to current directory.

NB: The energy interval [E1,E2] is used for determining r0 only.
The energy interval of output is the one of the first flux file.
------------------------------------------------------------------------------

jinpol84 interpolates the (reaction) probabilities for different total J. As input it needs two sets of probabilities (e.g. flux files) for J1 and J2. It then creates probabilities for all J with J1 < J < J2. The algorithm has some similarity with J-shifting, but is much more accurate. This J-interpolation algorithm is described in an appendix of J.Chem.Phys. 116 (2002) 10641.
NB: The input files may be in flux-file format or just two-columns ASCII files (energy,probability).
NB: The energy unit is assumed to be eV. Otherwise the R0 value is incorrect, but the interpolation will work. Likewise, if a wrong mass is put in, only the printed value of R0 will be affected.

joinpsiver name1, name2, ...

Purpose : Join psi files have been calculated with different number of SPFs.
Usage   : joinpsi84 [-w -D -o -l -f -p -ver -h -?] file1 file2 ...

Options      :
-w           : Overwrite existing output file.
-D PATH      : The output is directed to PATH.
The string PATH may be a relative or a full path-name.
-o NAME      : Change the name of output file from "psi" to NAME.
-l LIST      : The list if psi files is read from file LIST, rather from
the input line.
-f n1 n2 ... : Force the number of SPFs to be equal to n1,n2,... .
Otherwise the maximum number from psi files will be chosen.
-p           : Print information to the screen.
-ver         : Version information about the program
-h           : Print this help text.
-?           : Print this help text.

Examples     : joinpsi84 ../dir1/psi ../dir2/psi ../dir3/psi
joinpsi84 -w -f 30 40 20 -l list_psi
------------------------------------------------------------------------------
The program is useful, if the maximum number of SPFs is not needed for all propagation times (for example in the scattering problem). One may perform several calculations with varying numbers of SPFs (usually with smaller numbers of SPFs at the beginning and at the end of the propagation, if it is a scattering process), and then joins the psi-files with the aid of this program for the further analysis (e.g. flux analysis, etc.). The number of SPFs for the resultant psi-file equals the maximum from all psi-files. The artificially augmented SPFs are filled with zeros. The program is working only with psi=single or psi=double formats.

maxpever name

Purpose : Calculates the minimum and maximum value of a PES
Usage   : maxpe84 [-h -ver ] path-to-name-directory
First generate an oper file for an exact calculation and then
run maxpe84 with the path to that directory as argument.

-ver    : Version information about the program
-h      : Print this help text
NB: The routine vminmax is a more powerful alternative which
does not require an operator in exact, but in pes format
-------------------------------------------------------------

Calculates the maximum and minimum value of a PES. This is useful to check for the presence of large eigenvalues, which affect the efficiency of a numerically exact propagation.
• Generate an oper file ready for a numerically exact calculation, i.e. the keyword exact should be present in the RUN-SECTION.
• Run maxpe84. Output is to the screen.

Note that vminmax84 is a more powerful alternative which does not require the operator in exact, but in pes format.

mmemtestver

Purpose : Checks, how much memory can be allocated.
Usage   : mmemtest84 [-h ]
Options : -h    Print this help text.
Run mmemtest84 on the computer, on which you want to run a
program of the MCTDH package. Note, that only mctdh84 and potfit84
can use up to 8 GB (on 64 bit architectures with sufficient memory).
All other MCTDH-programs are limited to a allocate at most 2 GB.
------------------------------------------------------------------

normver

Purpose: Calculates the norm and A-norm of a wavefunction.
Usage:   norm [-f -i -n -r -ver -h -?] .

Options :
-f FILE : The wavefunction is read from file FILE
rather than from ./psi
The string FILE may be a relative or a full path-name.
-i DIR  : data is stored in directory DIR
-n step : Compute the norm only every step-th output.
-r      : Take the restart file rather than the psi file.
-ver    : Version information about the program
-h      : Print this help text.
-?      : Print this help text.
------------------------------------------------------------------------------

npotminmaxver

Purpose  : Find lowest and largest points on a surface generated by one or
more natpot files. This is similar to vminmax84 but skipping the
generation of a PES file.

Usage    : npotminmax84 [-D -mnd -t -c -srf -d -ver -w -h -?] <inp>

Options  :
-D DIR  : The output is directed to DIR. This option overwrites "name"
in the run-section of the input file.
-mnd    : Enables creating a name directory.
-t FILE : Do not run over the complete grid but use dvr indices in FILE
(n blank separated integers per line).
-c FILE : Substract potential in FILE (one energy per line, like readsrf
input) before evaluation.
-d DIR  : The new dvr is read from directory DIR.
-srf    : Write potential to ASCII file "srf.asc" (one energy per line
-n num  : Print num maximal and minimal values.
The default is num = 10. Maximum is 100.
-ver    : Version information about the program
-w      : Existing output is overwritten.
-h -?   : Print this help text.
------------------------------------------------------------------------------
The program is similar to vminmax, except that the generation of a PES file is skipped and a set of natpot files is used to generate the potential operator. The program requires an input file which contains at least a RUN-SECTION specifying job parameters and a NATPOT-SECTION listing paths to natpot files. In addition, either the keyword readdvr in the RUN-SECTION needs to be set or a PRIMITIVE-BASIS-SECTION needs to be present in the input file. Note, that the system DVR may contain more degrees of freedom than the natpot files are defined on. Furthermore, the natpot files can have different sets of degrees of fredom compared to each other as long as all degrees of freedom are contained in the system DVR, i.e., it is similar to constructing the operator-file from different natpot files.

In Run-Section following keywords are possible:

Keywords of npotminmax RUN-SECTION
Keyword     compulsory / optional Description equivalent option
name = S   C The name of the calculation. A directory with path S is required in which files will be written.   -D
readdvr = S   C/O The DVR file containing all system DOF is located in directory S. Compulsory if there is no PRIMITIVE-BASIS-SECTION.   -d
trajectory = S   O The scan is not done over all grid points but only over a subset. S is the path to an ASCII file containing the DVR indices of the grid points, where one grid point is given per line. Note, that the indexing is assumed to be in the same order as the DOF appear in the PRIMITIVE-BASIS-SECTION, i.e., the first index belongs the first coordinate in the PRIMITIVE-BASIS-SECTION, etc. The indices are assigned to the natpotfiles by their modelabels as in the Hamiltonian generation.

A trajectory of DVR indices can be either created using the chkpes program or (more conveniently) with mcpotfit.

-t
compare = S   O Compare the sum of all natpot files to a surface file S. That is, after summing the contributions of all natpot files for the n-th grid point, substract the value from the n-th line in S. If the -t option is given, S must contain at least as many lines as the trajectory file where ist is assumed that the n-th line of S corresponds to the n-th line of the trajectory. If the -t option is not given, S must contain as many lines as there are grid points in the system DVR, where the index of the first coordinate runs fastest, the index of the second coordinate second fastest etc.

Note, that the chkpes as well as mcpotfit programm produce surface files along a dvrindex trajectory as used for the trajectory = S keyword.

-c
wrsrf   O Write the sum of all natpot files to a surface file srf.dat in the name directory.   -srf
overwrite   O Files in name directory will be overwritten.   -w

The NATPOT-SECTION only contains a list of paths to directories holding a natpot file.

Example:

RUN-SECTION
name       = npmm_Q1_Q2_Q3
trajectory = Q1_Q2_Q3/dvrindex
compare    = Q1_Q2_Q3/energies
END-RUN-SECTION

NATPOT-SECTION
../POTFITS/pf_Q1_Q2   # directory with natpot of Q1 and Q2
../POTFITS/pf_Q2_Q3   # directory with natpot of Q2 and Q3
../POTFITS/pf_Q1_Q3   # directory with natpot of Q1 and Q3
end-natpot-section

Output files (in name directory)

npotmm.log
The log file.

dvr
The primitive grid definitions. (Only if a PRIMITIVE-BASIS-SECTION is given in the input file.)

largest.dat
The set ot N largest values found, the coordinates as well as contributions of the single natpot files.

smallest.dat
The set ot N smallest values found, the coordinates as well as contributions of the single natpot files.

smallest_abs.dat
The set ot N smallest absolute values found, the coordinates as well as contributions of the single natpot files.

natpot_statistics
For each natpot file its mean and RMS values as well as largest, smallest and smallest absolute values found plus their coordinates

orthover

Purpose: Checks the orthonormality of the SPFs.
Usage:   ortho<ver> [-f -i -n -r -nofs -ver -h -?] .

Options :
-f FILE : The wavefunction is read from file FILE
rather than from ./psi
The string FILE may be a relative or a full path-name.
-i FILE : Same as -f FILE.
-n step : Compute the orthonormality error only every step-th output.
-r      : Take the restart file rather than the psi file.
-nofs   : No transformation to fs. Use when "time-not-fs" was set in mctdh.
-ver    : Version information about the program
-h      : Print this help text.
-?      : Print this help text.
------------------------------------------------------------------------------

The output is to the screen and shows the square root of (1/n) times the double-sum over the squares of the matrix-elements of (ov-1) where ov denotes the overlap-matrix of the SPFs.

overlapver

Purpose: Calculates the difference between two wavefunctions.
See review Eqs.(159-163) for a definition of the output.
Usage  : overlap84 [-f -i -o -a -b -c -t -n -nofs -w -ig -ver -h -?] [Psi] .

Options :
-f FILE : The wavefunction is read from file FILE rather than from ./psi
The string FILE may be a relative or a full path-name.
-i FILE : Same as -f FILE.
-o FILE : The overlaps are written to file FILE.pl rather than to ./ovl_XXX
The string FILE may be a relative or a full path-name.
If FILE=no, i.e. "-o no", no ovl-file will be opened.
-nofs   : No transformation to fs. Use when "time-not-fs" was set in mctdh.
-a      : Compute the overlaps according to their position on psi-file,
i.e. times (tpsi) are ignored.
-b      : Compute the overlaps <Psi(t_final-t)|Psi1(t)>,
i.e. the bra-state runs backwards in time.
-c      : Compute the overlaps <dconjg(Psi(t))|Psi1(t)>,
i.e. the complex-symmetric scalar product is used.
-s s s1 : Compute the overlaps <Psi(t)_s|Psi1(t)_s1>,
s and s1 denote the electronic states of bra- and ket-WF.
-t tend : Compute the overlap only up to t=tend.
-n step : Compute the overlap only every step-th output.
-w      : An existing overlap file is overwritten.
-ig     : Ignore "different bases" error-msg.
-ver    : Version information about the program
-h      : Print this help text.
-?      : Print this help text.

The argument Psi (bra-states) is the name of the file containing
the comparison wavefunction.
The ket-states are read from ./psi or from the file given by the -f option.
If no argument is given, the user is prompted for the missing argument.
NB: The overlap is written to the file ./ovl_XXX, the screen-output
shows error=||Psi-Psi1||. See chapter 8.2 of the MCTDH review.
------------------------------------------------------------------------------
Calculates the error (overlap) between the wavefunctions in ./psi (or FILE) and comp_file . The results are written to the screen and a file, the name of which depends on the input file names, or the argument of the -o option. E.g.

overlap84 ../dir/psi
calculates the error between the wavefunctions in ./psi and ../dir/psi, and the results is stored in ./ovl_dir (unless the -o option was used). The WFs on ../dir/psi are taken as bra-states and the WFs of the psi-file of the current directory are the ket-states.
In contrast
overlap84 -f psi1 psi2
would calculate the error between the wavefunctions in ./psi1 and ./psi2, storing the results in ./ovl_psi2.

The output (on screen) is as described in the review Eqs.(159-163). In particular: Error = ||psi1-psi2|| = Eq.(159); Hilbert = arccos( |<psi1|psi2>|/(||psi1||*||psi2||)) = Eq.(163); Phase = arctan( Im(<psi1|psi2>|)/ Re(<psi1|psi2>|) = Eq.(161); Ph. corr = Eq.(162).
The overlaps <psi1|psi2> <psi1|psi1> and <psi2|psi2> are written to the file ovl_dir, where dir denotes the name-directory of the second calculation. If this file exist and when the option -w is not given, overlap84 will read the file ovl_dir rather than re-do the calculation.

pexpectver

Purpose: Calculates the time evolution of an expectation value.
Usage  : pexpect [-f -i -o -g -a -r -t -p -u -n -w -ver -h -?] [mode operator].

Options :
-f FILE : The operator is read from file FILE rather than from ./oper
The string FILE may be a relative or a full path-name.
-i FILE : The one particle density is read from file FILE rather
than from ./pdensity
The string FILE may be a relative or a full path-name.
-o FILE : The output is written to file FILE.pl rather than to ./pexpect.pl
The string FILE may be a relative or a full path-name.
If FILE=no, i.e. "-o no", no output file will be opened.
-g      : GNUplot comand lines are written to the output file(s).
-a      : GNUplot is called automatically. (sets -g and -w).
-p      : Print GNUplot to lpr (sets -g).
-r      : Do not divide the expectation values by norm^2.
-t tfin : The expectation values are computed only up to time=tfin.
-u UNIT : The unit used is UNITs (Default is a.u.).
-n step : Compute the property only every step-th output.
-w      : An existing property file is overwritten.
-ver    : Version information about the program.
-h      : Print this help text.
-?      : Print this help text.

The first argument "mode" is the (combined) mode for which the
expectation value is calculated.
The second argument "operator" is the name of the operator, as specified in an
HAMILTONIAN-SECTION_xxx, for which the expectation value is required.
If "system" is chosen as operator the uncorrelated energy of the specified
mode is computed.
If no argument is given, the user is prompted for the missing argument(s).
------------------------------------------------------------------------------
The first argument "mode" is the (combined) mode for which the expectation value is calculated. The second argument "operator" is the name of the operator, as specified in an HAMILTONIAN-SECTION_xxx, for which the expectation value is required. If "system" is chosen as operator the uncorrelated energy of the specified mode is computed. If no argument is given, the user is prompted for the missing argument(s). The operator must first have been processed by the MCTDH program. This means that it must first be set up in a .op file. If it is known before a propagation is made, this operator can be included with the system operator in the main oper file (by using a named HAMILTONIAN-SECTION_XXX). Alternatively, e.g. if a propagation has already been made, the file oper_name can be generated using the genoper=name keyword.

Start the program in the directory containing the one-particle density in a pdensity file. If no operator is given, the program looks in the oper file to see which operators are there and asks which is to be evaluated. The -o option can use a file other than this. The expectation value of this operator is then calculated for each time, and output to the file pexpect.pl, unless otherwise instructed. The -a option enables the automatic plotting of the results. Note that pexpect computes expectation values of one-particle operators. (For combined modes this may include more than one degree of freedom). If the specified operator acts on several modes, pexpect uses only the uncorrelated part for the specified mode.

prhosubver

Purpose: Build and output reduced density from pdensity file.
Usage  : prhosub [-i -n -nw -k -q -ver -h -?] [mode [state]]

Options :
-i FILE : The one particle density is read from file FILE rather
than from ./pdensity
The string FILE may be a relative or a full path-name.
-n step : Compute the density only every step-th output.
-nw     : No weights are employed. Write the "naked" DVR populations.
-k min max: Write the momentum spectrum for the interval min<p<max.
If "-k d" is given, default values for min and max will be taken.
-q pts  : Use the number of pts points for representing the momentum
spectrum. Default: pts =   501
This option is ignored, when -k is not given.
-ver    : Version information about the program.
-h      : Print this help text.
-?      : Print this help text.

The argument "mode" is the mode for which the reduced density is
calculated. "state" is the electronic state.
If an argument is not given, the the default value 1 will be taken.
------------------------------------------------------------------------------
The pdensity file must have been produced with the line pdensity=I in the run-section. Here, I is the number of the degree of freedom of interest. pdensity=I will then contain the reduced density in natural orbital representation. prhosubver produces one file of output for each tout, with the reduced density represented on the primitive grid. The columns are then x, x', real part, imag. part, absolute value.

Start the program in the directory containing the one-particle density in a pdensity file.

probsqver

Purpose: Computation of the quadratic entropy. sum p_i^2.

Usage : probsq84 [-i -t -a -ver -h -?]

Options :
-i FILE  : The autocorrelation is read from file FILE
rather than from ./auto
The string FILE may be a relative or full path-name.
-a tstart: The autocorrelation function is read only from t=tstart on.
-t tcut  : The autocorrelation function is read only up to t=tcut
-ver     : Version information about the program
-h       : Print this help text.
-?       : Print this help text.
------------------------------------------------------------------------------
The program computes the integral (1/T)*Int[0,T] g(t)*|c(t)|**2 dt where T denotes the final time of the autocorrelation c(t) and g(t) is some weighting function. g(t) = 1, (pi/2)*sin(pi*t/T), (pi/2)*cos(pi*t/2T), 2*sin(pi*t/T)**2. For large T this converges to Sum pi2 , where pi denotes the occupation probability of the i-th eigenstate. This sum is related to the so called quadratic entropy. Output are the values Sum pi2, 1/(Sum pi2), and -Log2(Sum pi2). If one assumes that all pi = p for i=1 to N, and zero otherwise, then 1/(Sum pi2) equals the number of states N. (Note Sum pi=1, hence p=1/N).

psi2exver

Purpose: Converts the wavefunction in a psi file from MCTDH
to numerically exact form. Default is single precision.
To inspect the WF, the output may also be in ASCII format
with or without printing coordinates (option -ort).
DVR-weights can be removed from the WF by setting the
option -W (requires that -ort is also set).

Usage  : psi2ex84 [-f -i -o -stp -n -skip -thrs -ort -ascii -W -pop2
-plot -abs -abs2 -rst -w -dp -ver -h -?]

Options :
-f FILE : The wavefunction is read from file FILE rather than from ./psi
The string FILE may be a relative or a full path-name.
-i FILE : Same as -f FILE.
-o FILE : The output is written to file FILE rather than to ./psi.ex
The string FILE may be a relative or a full path-name.
-stp step: Compute the exact WF only every step-th output.
-n num  : only the first num WFs are read  from psi-file.
-skip m : the first m WFs are skipped.
-pop2   : exaxt WF is given in second population
-plot   : only meaningful if -pop2 is activated
and some dof is of fft type.
For each dof of fft type, swaps the
x axis in a convenient way for further
plotting. NB: the resulting WF cannot be
used in a subsequent MCTDH calculation.
-ort    : The WF and their coordinates are written to output file.
Note that -ort sets -ascii.
-W      : The DVR-weights are removed from the WF.
Requires that -ort is set, otherwise ignored.
-thrs thrs : The WF and their coordinates are written
only if(abs(WF).GT.thrs). Works only with -ort.
-ascii  : The output file is in ascii rather than binary format.
In this case the (dummy) A-vector and psiinfo
are not written, i.e. the file contains the WF only.
-abs    : The absolute value of the exact WF is written.
In this case the (dummy) A-vector and psiinfo
are not written, i.e. the file contains abs(WF) only.
-abs2   : Similar to -abs, but the absolute value squared
of the exact WF is written.
-dp     : The (binary) output-file is written in double precision.
-w      : An existing output file is overwritten.
-ver    : Version information about the program
-h      : Print this help text.
-?      : Print this help text.
------------------------------------------------------------------------------
The program reads the MCTDH wavefunction in a psi file and converts it to the numerically exact form, represented on the full grid. This is useful for other analyse routines. The output is by default single precision to the file psi.ex. The argument -dp forces double precision output, and -o FILE changes the output file name. The argument -pop2 gives the exact wavefunction in second population. In case some of the degrees of freedom is of fft type, the x axis (which in such a case is the momentum axis) can be swapped using the option -plot, so that it goes from negative to positive values and can nicely be plotted. In case -plot is set, the momentum grid ranges from -(gdim/2)*dp to (gdim/2-1)*dp in steps of dp, where dp=2*pi/(dx*gdim). The values for dx (grid spacing) and gdim (number of grid points) are listed in the log file. In case -plot is not set, the momentum grid runs form zero to (gdim/2-1)*dp and then from -(gdim/2)*dp to -dp.

The format of the psi file is explained in the MCTDH/Output Documentation.
When the option -abs is given, the absolute value of the wavefunction is written. In this case the file contains only the wavefunction (i.e. no header, no dummy A-vector). The order of the elements is like in a Fortran array: Psi(dof_1,dof_2,...,dof_f). To inspect the WF it may be written in ASCII format (-ascii). Furthermore, the coordinates may be printed together with the WF-values to an ASCII output file (-ort) and DVR-weights may be removed from the WF (-W, requires that -ort is set).

psi2rstver

Purpose: Converts the wavefunction of a psi file to restart format.
The argument "step" is the number of the WF to be converted.

Usage  : psi2rst84 [-f -i -o -w  -ver -h -?]  step

Options :
-f FILE: The wavefunction is read from file FILE rather than from ./psi
The string FILE may be a relative or a full path-name.
-i DIR : The wavefunction is read from file DIR/psi rather than from ./psi
-o DIR : The output is written to file DIR/restart rather than to ./restart
-w     : An existing output file is overwritten.
-ver   : Version information about the program
-h     : Print this help text.
-?     : Print this help text.

NB: The output is written to the file restart in the current directory,
or to DIR/restart, if the option -o is given.
To receive some information on the WF file, run "psi2rst nn", where
nn is zero or a large number (larger than the number of WFs on file).
------------------------------------------------------------------------------

rdacoeffver

------------------------------------------------------------------------------
Purpose: Prints  the A-vector of a wavefunction or restart file.
Usage  : rdacoeff [-f -i -n -n1 -n2 -s -max -thr -uplo -r -abs -abs2
-re -im -t -at -count -ver -h -?] .
Options:
-f FILE : The wavefunction is read from file
FILE rather than from ./psi
The string FILE may be a relative or a full path-name.
-i DIR  : data is read from directory DIR
-n step : Print the A-vector only every step-th output.
-thr thr: Print only coefficients whose absolute values
are larger than thr (default: 10^-8).
-max mx : Print only coefficients whose absolute values
are smaller or equal than mx.
-uplo ul: Print only coefficients from the upper (ul="U") or lower (ul="L")
"triangle" of the A-tensor.
If ul="L" is set, only A_j1_j2... with
(j1-1)/(n1-1) + (j2-1)/(n2-1) + ... <= 1 are taken into account,
if ul="U" is set, only A_j1_j2... with
(j1-1)/(n1-1) + (j2-1)/(n2-1) + ... > 1 are taken into account.
-n1  n1 : Print only coefficients with J.ge.n1.
-n2  n2 : Print only coefficients with J.le.n2.
-t   tm : Print only WFs with time .le. tm.
-at atm : Print only WFs with time .ge. atm.
-s  s   : Print only coefficients for state s.
-count  : Count the number of coeffs between thr and max.
-r      : Take the restart file rather than the psi file.
-rst    : Same as -r.
-abs    : Print abs(psi).
-abs2   : Print abs(psi)**2.
-re     : Print Re(psi).
-im     : Print Im(psi).
-ver    : Version information about the program
-h      : Print this help text.
-?      : Print this help text.

It is in general useful to pipe the output into less: rdacoeff84 | less
or to re-direct it into a file:  rdacoeff84 > acoeff
------------------------------------------------------------------------------

rdautoever [-f -i -o -d -t -h-?] ns nmd

Purpose : Analyses the error measures of the autocorrelation.
Usage   : rdautoe84 [-f -i -o -d -t -ver -h -?] [ns  mode].

Options :
-f FILE : the error measures are read from file FILE
rather than from ./autoe
The string FILE may be a relative or full path-name.
-i FILE : Same as -f FILE
-o FILE : The error measures are written to file FILE.pl
rather than to ./autoe.pl
The string FILE may be a relative or full path-name.
If FILE=no, i.e."-o no", no autoe.pl file will be opened.
-t dt   : Print the error measures every dt fs.
-d      : Write the difference of the autocorrelations to output file.
-g      : GNUplot command lines are written to the output file(s).
-ver    : Version information about the program
-h      : Print this help text.
-?      : Print this help text.

There are two arguments:
The first argument is ns, i.e. the number of electronic state to be considered.
The second argument is nmd, i.e. the runing number of the mode to be
considered. (nmd=0 -> All modes are decreased by one spf).
If not enough arguments are given, the user is prompted for the
missing arguments.
------------------------------------------------------------------------------
The routine rdautoe reads the autoe-file and writes the modified autocorrelation function to the output file (e.g. autoe.pl) which has exactly the same format as the auto file and hence may be read be the routine autospec84. By modified autocorrelation function we denote the autocorrelation which is calculated by dropping the most weakly occupied spf of mode nmd and state ns. A spf in each mode (and state) is dropped, if nmd=0. If the option -d is given, the difference of the modified autocorrelation and the (original) autocorrelation is outputed.

rdcheckver [-f -i -oc -on -g -gq -t -h -?] ns nmd

Analyse program "rdcheck" (Read check file)
Analyses the state population and natural weights.
Usage   : rdcheck84 [-f -i -oc -on -e -g -gq -t -n -ver -h -?] ns nmd.

Options :
-f FILE : The data is read from file FILE rather than from ./check
The string FILE may be a relative or full path-name.
-i DIR  : Data is read from file DIR/check
-oc FILE: The state populations are written to file FILE.pl
rather than to ./chk.pl
The string FILE may be a relative or full path-name.
If FILE=no, i.e."-oc no", no chk.pl file will be opened.
-on FILE: The natural populations are written  to file FILE.pl
rather than to ./nat.pl
The string FILE may be a relative or full path-name.
If FILE=no, i.e."-on no", no nat.pl file will be opened.
-oq FILE: The expextation values <q> and <dq> are written to file FILE.pl
rather than to ./qdq.pl
The string FILE may be a relative or full path-name.
If FILE=no, i.e."-oq no", no qdq.pl file will be opened.
-e      : State energies are written to screen.
-tr2    : Trace{rho^2} is written to screen. (MCTDH 1 part. dens.)
Do not use -e, -tr2, or -entropy simultaniously.
-entropy: Entropy of natural populations is written to screen.
Do not use -e, -tr2, or -entropy simultaniously.
-g      : GNUplot command lines are written to the output file(s).
-gq  n  : Selects the n-th dof for th qdq-plot. Sets -g.
-r      : Renormalise the natural weights via division by state population.
-s      : Suppress informative-only part of output (short output).
-t dt   : Print the natural weights every dt fs.
-last   : Take last WF when computing the
maximum over time of lowest nat.-weight.
-skip n : Skip the first n WFs when computing the
maximum over time of lowest nat.-weight.
-cut n  : Consider only the first n WFs when computing the
maximum over time of lowest nat.-weight.
-n      : No units. Sets fs=1. Use, if 'time-not-fs' is set.
-ver    : Version information about the program.
-h      : Print this help text.
-?      : Print this help text.

There are two arguments:
The first argument is ns, i.e. the state for which the natpops are written.
(ns=0 -> no output files written).
The second argument is nmd, i.e. the runing number of the mode
for which the natpops are written.
If not enough arguments are given, the user is prompted for the
missing arguments.
------------------------------------------------------------------------------
If one just wants to check the convergence of the natural weights, type: rdcheck84 0 , or, to see the development of the lowest natural weights in timesteps of 1.2 fs: rdcheck84 -t 1.2 0
The entropies of the natural populations are written to screen with the command rdcheck84 -entropy s 0, where s is to be substituted by the electronic state to be inspected. If there is only one state, s may be dropped. Similar information is obtained with the -tr2 option.
entropy = sum -pi Log(pi)
tr2 = sum pi2
where pi denotes the normalized (by norm2) i-th natural weight.

rdgridver dof

Analyse program "rdgrid" (Read grid points)
Purpose : Read the DVR grid from dvr file and print it.
Usage   : rdgrid [ -f -i -s -oc -w -pop2 -h -? ]  dof .

Options :
-f FILE: the DVR file is read from file FILE rather
than from ./dvr .
The string FILE may be a relative or full path-name.
-i DIR : dvr file is read from directory DIR.
-s     : Additional information is suppressed.
-oc    : One-column output enforeced.
-w     : DVR weights are also printed.
-pop2  : Grids of 2nd population are printed.
-ver   : Version information about the program.
-h     : Print this help text.
-?     : Print this help text.

The argument dof is the number of the degree of freedom for which
the grid is printed.
If no argument is given, the user is prompted for the missing argument.
------------------------------------------------------------------------------
This program is useful to obtain a quick overview of the actual locations of the grid points. Note that the programs showsys and showpot can show cuts of multi-dimensional potential energy surfaces by cutting at grid points only. The input values for the cuts are hence changed by these programs to take the value of the nearest grid point.
With the option -pop2 rdgrid prints the grids of 2nd popuation.
To plot the DVR points or weights use the -s and -oc or -w options and re-direct the output to a file. The data can then be plotted by plgen. E.g. plotting coordinate 2 vs grid:
rdgpop84 -s -oc 2 > file
plgen file
Or plotting weights vs coordinate:
rdgpop84 -s -w 2 > file
plgen file 2:3

rdgpopver nz dof

Analyse program "rdgpop" (Read grid populations)
Purpose : Analyses the grid population read from file gridpop.
Usage   : rdgpop84 [-f -i -g -o -t -skip -cut -n -ver -h -?] [nz  dof].

Options :
-f FILE : the grid populations (densities) are read from file FILE
rather than from ./gridpop
The string FILE may be a relative or full path-name.
-i FILE : Same as -f FILE
-o FILE : The grid populations are written to file FILE.pl
rather than to ./gpop.pl
The string FILE may be a relative or full path-name.
If FILE=no, i.e."-o no", no gpop.pl file will be opened.
-t dt   : Print the maximal occupations every dt fs.
-skip n : Skip the first n WFs when computing the
maximum over time of the grid populations.
-cut n  : Consider only the first n WFs when computing the
maximum over time of the grid populations
-n      : No units. To be used when 'time-not-fs' is set.
-g      : GNUplot command lines are written to the output file(s).
-ver    : Version information about the program
-h      : Print this help text.
-?      : Print this help text.

There are two arguments:
The first argument is nz, i.e. the number of grid points sumed over.
The second argument is dof, i.e. the runing number of the degree
of freedom to be considered. (dof=0 -> no output to file gpop.pl).
If not enough arguments are given, the user is prompted for the
missing arguments.
------------------------------------------------------------------------------
The grid population is summed over the first nz and the last nz points of both the spatial (DVR) grid and the basis (FBR) occupations. By varying nz one may study how much the grid may be shortened. The integer dof specifies the degree of freedom to be investigated. Depending on the MCTDH input, the grid populations are summed over all electronic states or (if gridpop=el is given in the run-section) treated separately for each electronic state. For a simple check on the convergence with respect to the primitive grids simply type rdgpop84 1 0 while being in the name directory.

rdstepsver

Purpose : Reads the step-file and evaluates statistic measures
of the integrator steps.
Usage   : rdsteps84 [-f -i -ver -h -?].

Options :
-f FILE : the data is read from file FILE
rather than from ./steps
The string FILE may be a relative or full path-name.
-i DIR : data is stored in directory DIR
-ver    : Version information about the program
-h      : Print this help text.
-?      : Print this help text.
------------------------------------------------------------------------------

rdnatpotver

Purpose  : Prints information on a natpot to file rdnatpot.log
On demand it prints natural potentials and D-tensor.

Usage    : rdnatpot84 [-i -D -n -evec -b -contr -ver -h -?]

Options  :
-i NAME : The files are read/stored in the directory NAME
-D DIR  : The output is written to directory DIR. This is useful,
if one has no right to write to the name directory.
The directory DIR must exist.
The string DIR may be a relative or a full path-name.
-n FILE : The potential fit is read from file FILE
If "-n FILE" is not given, FILE=natpot is assumed.
-evec M : The eigenvectors (natural potentials)
of mode M are written to ASCII file evec.
M must not equal modc (contracted mode).
If M=99 then natpots of all modes (except modc) are written.
-b      : Eigenvectors (natural potentials) are additionally
written to binary file evecb, (if -evec is set).
-contr  : The contracted natural potentials of the
contracted mode are written to ASCII file contr.
WARNING: This may be a lot of data.
-ver    : Version information about the program
-h      : Print this help text.
-?      : Print this help text.
------------------------------------------------------------------------------
This program is useful to obtain a quick overview on what is on a given natpot file. Information on the grids, mode combinations, and numbers of SPPs are written to rdnatpot.log. On demand the program writes the natural potentials (SPPs) and the D-tensor to ASCII files.

rdupdatever

Purpose : Reads the update steps from file "update"
Usage   : rdupdate [-f -i -inter -ver -h -?].
Options :
-f FILE : the data is read from file FILE
rather than from ./update
The string FILE may be a relative or full path-name.
-i DIR : data is stored in directory DIR
-inter : enables interactive plotting
-ver    : Version information about the program
-h      : Print this help text.
-?      : Print this help text.
------------------------------------------------------------------------------
Reads the CMF step sizes from file, and removes repetition steps. The use of the -inter argument starts an menu to enable interactive plotting of the data. Without this argument, the data is written to the screen.

rdautover

Purpose : Reads the autocorrelation from file "auto"
Usage   : rdauto [-f -i -inter -ver -h -?].
Options :
-f FILE : the data is read from file FILE
rather than from ./auto
The string FILE may be a relative or full path-name.
-i DIR : data is stored in directory DIR
-inter : enables interactive plotting
-ver    : Version information about the program
-h      : Print this help text.
-?      : Print this help text.
----------------------------------------------------------------------------
The use of the -inter argument starts an menu to enable interactive plotting of the data. Without this argument, the data is written to the screen.

rdeigvalver

Purpose : Reads the eigvalues from file "eigval"
Usage   : rdeigval [-f -i -inter -ver -h -?].
Options :
-f FILE : the data is read from file FILE
rather than from ./eigval
The string FILE may be a relative or full path-name.
-i DIR  : data is stored in directory DIR
-inter  : enables interactive plotting
-ver    : Version information about the program
-h      : Print this help text.
-?      : Print this help text.
----------------------------------------------------------------------------
The use of the -inter argument starts an menu to enable interactive plotting of the data. Without this argument, the data is written to the screen.

showd1dver [-T|-S|-I|-P|-M|-Ma] [-f -i -o -D -a -g -p -x -y -fac -n -w -sm -nw -pop2 -G -inter -ver -h -?] fx sx

Purpose: To selectively plot 1-dimensional densitites.
Usage  : showd1d84 [-T|-S|-I|-P|-E|-M|-Ma]
[-f -i -o -D -a -g -p -x -y -fac -n -no -r -t -skip
-nofs -l -sm -w -nw -pop2 -G -inter -ver -h -?] fx sx.

Output Formats :
-S      : "Step through" by pressing RETURN. (default).
-M      : "Movie", one picture every second.
-I id   : "Index file", id=0 shows the first picture, id=1 the second, etc.
-P pnt  : "Point file", shows the population of the pnt-th grid-point over time.
-E      : Similar to -P, but shows the population of all grid-point over time.
May be used to plot state populations of a single-set run.
-T      : "3D time file" (grid file).
-Ma     : Mathematica file (for dynamic plots).

Options  :
-f FILE : The 1D-density is read from file FILE rather than from ./gridpop
The string FILE may be a relative or a full path-name.
-i DIR  : Use DIR as name-directory.
-o FILE : The output is written to file FILE_f* rather than
to ./den1d_f*
The string FILE may be a relative or a full path-name.
If FILE=no, i.e. "-o no", no output file will be opened.
-D DIR  : The output is written to directory DIR. This is useful,
if one has no right to write to the name directory.
The directory DIR must exist.
The string DIR may be a relative or a full path-name.
-pop2   : 2nd population plotted (i.e. basis set occupations).
-g      : GNUplot comand lines are written to the output file(s).
-a      : GNUplot is called automatically. (sets -g and -w).
-p      : Print GNUplot (sets -g). Allowed only with -I or -T.
-sm     : Smooth the plot-data by splines. (Only for GNUplot, sets -g).
-l      : Use logarithmic scale. (Only for GNUplot).
-n step : Show only every step-th picture. (Only useful for -S, -M).
-t tfinal : Show densities only up to t=tfinal.
-skip n : Skip the first n densities.
-r      : Re-normalize densities.
-nw     : No weights are employed. Plot the "naked" DVR populations.
-no     : No output file opened.
-nofs   : No transformation to fs. Use when 'time-not-fs' was set in mctdh.
-x xmin xmax : Set the abscissa length (default: auto).
-y ymax : Set the upper bound of the ordinate (default: auto).
-z ymin : Set the lower bound of the ordinate (default: 0 or 1d-10 (for -l)).
-x, -y, -z options are ignored when -T or -Ma is set.
-fac factor : Multiply the density values by factor.
-G      : Plot gridlines.
-w      : An existing output file is overwritten.
-h      : Print this help text.
-?      : Print this help text.
-ver    : Version information about the program

The arguments fx and sx select the dof and state (x=integer).
With the argument s0 a sum over all states is performed.
If sx is not given, s1 is assumed.
If no argument is given, the user is prompted for the missing argument.
EXAMPLES: showd1d84 -a f1
showd1d84 -a -M -y 0.3  -n 2  f2 s1
showd1d84 -w -p -I 20 f1
showd1d84 -a -p -T f3
------------------------------------------------------------------------------

For instance, showd1d84 -a f2 would plot the density for the 2nd degree of freedom by automatically calling GNUplot. If RETURN is pressed, the density for the next output-time is shown.

The output is to a file den1d_fx_sx or den1d_fx if no state is specified. Hence, the data may be analyzed or plotted by alternative plotting software.

The state is by default taken as 1. Thus if sx is not specified as an argument the 1st state is selected. Note that there is normally only one (effective) electronic state as the densities of the different states are summed before they are written to the file gridpop. If you want to avoid this and investigate the different electronic states separately, you have to provide the keyword "gridpop=el" in the mctdh input file.

For a multi-state single-set run showd1d can be used to plot the state populations. Use showd1d84 -a -E fx, where x is to be replaced by the DOF-number of the electronic degree of freedom.

The options -T, -S, -I, and -M selects the format for the output.
A Gnuplot index file writes a set of data for each time:
coordinate, Abs(spf), Re(spf), Im(spf)
with 2 blank lines between. Gnuplot can then select a set using the "index" option of the "plot" command. If -I id is specified only the plot number id will be plotted. (This is particularly useful when the plot should be printed. See the -p option). If -S is specified (or, since this is the default, if no format option is given) one can scroll through the pictures by pressing RETURN. If -M is specified one sees a "movie", one picture each second.
A Gnuplot grid file (option -T) writes a set of data for each time:
coordinate, time, Abs(spf), Re(spf), Im(spf)
with 1 blank lines between. Gnuplot can then plot a 3D plot of the spf evolution with time.
The mathematica option (-Ma) writes information to the screen to be read by a mathematica script (plspf.m). This enables dynamic pictures.

If one want to use showd1d to investigate the grid populations, one should use the option -nw (no weights) because the use of weights may lead to a wrong conclusion. The option -l (logarithmic plot) is also useful when checking grid populations.

showrstver [-f -i -o -a -g -abs -abs2 -re -im -reim -sm -p -y -nw -w -ver -h -?] fx sx

Purpose: Plot the single-particle functions of a restart file.
Usage  : showrst84 [-f -i -o -a -g -abs -abs2 -re -im -reim
-sm -p -x -y -nw -w -G -ver -h -?]  [fx sx] .

Options :
-f FILE : The wavefunction is read from file FILE rather than from ./restart
The string FILE may be a relative or a full path-name.
-i DIR  : Use DIR as name directory.
-o FILE : The plot data is written to file FILE_f* rather than to ./rst.pl_f*
The string FILE may be a relative or a full path-name.
If FILE=no, i.e. "-o no", no output-file will be opened.
-g      : GNUplot comand lines are written to the output file(s).
-a      : GNUplot is called automatically. (sets -g and -w).
-sm     : Smooth the plot-data by splines. (only for GNUplot, sets -g).
-p  id  : Print GNUplot (sets -g).
The integer id is the number of the spf to be printed.
-abs    : Plot absolute values. (default).
-abs2   : Plot absolute values squared.
-re     : Plot real part (only for GNUplot, sets -g).
-im     : Plot imaginary part (only for GNUplot, sets -g).
-reim   : Plot both imaginary and real parts (only for GNUplot, sets -g).
-nw     : No weights are employed. Plot the "naked" DVR values.
-x xmin xmax : Set the abszisse length.
-y ymax : Set the ordinate length (default: auto).
-G      : Plot gridlines.
-w      : An existing output file is overwritten.
-ver    : Version information about the program
-h      : Print this help text.
-?      : Print this help text.

The arguments select the dof (fx) and state (sx)
Replace x by appropriate number, any order is allowed.
If no argument is given, the user is prompted for the missing argument.
------------------------------------------------------------------------------
The state is by default taken as 1. Thus if sx is not specified as an argument the 1st state is selected.

E.g. showrst84 -a -re f2 would plot the real-part of the single-particle functions for the 2nd degree of freedom. GNUplot will be called automatically (-a option). The output is to the file rst.pl_f2.

showspfver [-spf -nat -f -i -o -a -g -p -x -y -sm -abs -abs2 -re -im -reim -n -w -h -?] fx sx x

Purpose: To selectively plot (uncombined) natorbs or spf's.
Usage  : showspf84 [-T|-S|-I|-M|-Ma] [-f -i -o -a -g -p -x -y
-sm -abs -abs2 -re -im -reim -n -w -nwi -nofs
-spf -nat -G -ver -h -?]  [fx sx x] .

Output Formats :
-S      : "Step through" by pressing RETURN. (default).
-M      : "Movie", one picture every second.
-I id   : "Index file", id=0 shows the first picture, id=1 the second, etc.
-T      : "3D time file" (grid file).
-Ma     : Mathematica file (for dynamic plots).

Options  :
-f FILE : The wavefunction is read from file FILE rather than from ./psi
The string FILE may be a relative or a full path-name.
-i DIR  : Use DIR as name directory.
-o FILE : The output is written to file FILE*_f* rather than
to ./natorb*_f*
The string FILE may be a relative or a full path-name.
If FILE=no, i.e. "-o no", no output-file will be opened.
-g      : GNUplot comand lines are written to the output file(s).
-a      : GNUplot is called automatically. (sets -g and -w).
-p      : Print GNUplot (sets -g). Allowed only with -I or -T.
-sm     : Smooth the plot-data by splines. (only for GNUplot, sets -g).
-abs    : Plot absolute values. (default).
-abs2   : Plot absolute values squared.
-re     : Plot real part (only for GNUplot, sets -g).
-im     : Plot imaginary part (only for GNUplot, sets -g).
-reim   : Plot both imaginary and real parts (only for GNUplot, sets -g).
-spf    : Plot the of the single-particle functions
rather than natorbs.
-nat    : Plot the of the natural orbitals rather
than spf's (default). Note x=1 -> highest occupied natorb.
-nofs   : No transformation to fs. Use when 'time-not-fs' was set in mctdh.
-nw     : No weights are employed. Plot the "naked" DVR values.
-n step : Show only every step-th picture. (Only useful for -S, -M).
-x xmin xmax : Set the abszisse length.
-y ymax : Set the ordinate length (default: auto). Ignored for -T and -Ma.
-G      : Plot gridlines.
-w      : An existing output file is overwritten.
-h      : Print this help text.
-?      : Print this help text.
-ver    : Version information about the program

The arguments select the spf (x), dof (fx) and state (sx)
Replace x by appropriate number, any order is allowed.
If no argument is given, the user is prompted for the missing argument.
------------------------------------------------------------------------------
E.g. showspf84 -a -spf f2 4 would plot the absolute value of the 4th single-particle function for the 2nd degree of freedom.

By default natural orbitals are output to a file natorbx_fx_sx. If the -spf flag is used, the single-particle functions are output to a file spfx_fx_sx.

The state is by default taken as 1. Thus if sx is not specified as an argument the 1st state is selected.

The options -T, -S, -I, and -M selects the format for the output.
A Gnuplot index file writes a set of data for each time:
coordinate, Abs(spf), Re(spf), Im(spf)
with 2 blank lines between. Gnuplot can then select a set using the "index" option of the "plot" command. If -I id is specified only the plot number id will be ploted. (This is particularly useful when the plot should be printed. See the -p option): If -S is specified (or, since this is the default, no format option is given) one can scroll through the pictures by pressing RETURN. If -M is specified one sees a "movie", one picture each second.
A Gnuplot grid file (option -T) writes a set of data for each time:
coordinate, time, Abs(spf), Re(spf), Im(spf)
with 1 blank lines between. Gnuplot can then plot a 3D plot of the spf evolution with time.
The mathematica option (-Ma) writes information to the screen to be read by a mathematica script (plspf.m). This enables dynamic pictures.

showsysver [-i -D -f -n -o -pes -rst -pop2 -ver -h -?]

Purpose: Enables 2D plotting of the wavefunction or PES.
Usage:   showsys84 [-i -D -f -p -o -pes -nopes -n -skip -step
-rst -nw -nofs -u -pop2 -pop2all -ver -h -?]

Options :
-i DIR  : data is read from name-directory DIR.
-f FILE : the wavefunction is read from FILE rather than from ./psi
The string FILE may be a relative or a full path-name.
-p FILE : the PES is read from FILE rather than from ./pes
The string FILE may be a relative or a full path-name.
-o FILE : The ouput is written to file FILE.pl rather than to ./surface.pl
The string FILE may be a relative or a full path-name.
-D DIR  : The output is written to directory DIR. This is useful,
if one has no right to write to the name directory.
The directory DIR must exist.
The string DIR may be a relative or a full path-name.
-pes    : only the PES file is read, i.e. no wavefunction plotting.
-nopes  : the PES file is not read, no potential plotting, no dia->ad trafo.
-rst    : the restart file is read, rather than the psi file.
-n num  : only num WFs are processed.
-skip m : the first m WFs are skipped.
-step step : only every step's WF will be processed.
-nw     : No weights are employed. Plot the "naked" DVR populations. (WF only).
-nofs   : No transformation to fs. Use when 'time-not-fs' was set in mctdh.
-u UNIT : The energy unit UNIT is applied. (default for pes: eV).
-pop2   : the basis set occupations are plotted rather than grid
populations. For FFT this implies momentum space representation.
Note: -pop2 sets -nw.
-pop2all: Similar to -pop2, but now all DOFs are transformed to second population.
This is only useful for cuts. Note: -pop2all sets -pop2 and -nw.
-ver    : Version information about the program.
-h      : Print this help text.
-?      : Print this help text.

The program is menu driven. Just type "showsys84" or "showsys84 -pop2"
If there are more WFs on psi file than you want to inspect, use the
options -n, -skip, or -step. Example: "showsys84 -skip 10 -n 2" will plot
the WFs no. 11 and 12 (i.e. t=10 and 11 if tpsi=1fs),
and "showsys84 -step 10 -n 2" will plot the WFs no. 1 and 11.
If one wants to plot a potential one should use "shwsys84 -pes",
because otherwise the WF is read, which consumes a large amount of memory.
NB1: The option -step may also be set by the menu-point 285.
The WFs to be processed may also be selected by setting time bounds
NB2: The options -pop2 and -pop2all are ignored for phifbr,sphfbr,k,external bases.
--------------------------------------------------------------------------------

Cuts of a potential energy surface (PES) can be generated from an pes file, and cuts and (1D or 2D) densities of the wavefunction can be generated from a psi file.
• Set up the operator in a .op file (see the documentation on Implementing a Hamiltonian).
• Use the MCTDH program with the genpes run-type to generate the dvr and pes files. It is often more convenient not to edit the input file but to run mctdh with the option -pes.
• Move to the directory where the pes file is, and start showpes84 -pes

If the (1D or 2D) densities of the wavefunction are to be displayed in the second representation, i.e. as basis set occupation, momentum representation or angular momentum representation, use the option -pop2 . (This switch is not available through the menu).

After starting the program a menu is presented with the help of which one can select a plot via an interactive process. The menu offers the following choices:

• Select (using "10 change plot task") between cuts of the adiabatic or diabatic potential, when a potential is plotted. When a wavefunction is plotted, one may chose between plotting cuts of the adiabatic wavefunction (abs), cuts of the diabatic wavefunction (abs) or the diabatic reduced density. The latter is the default. The reduced density is defined by integrating |psi|^2 over all variables except the one or two denoted by "x" or "x" and "y".
• Select a cut using "20 change coordinate section". This requires a coordinate value for each degree of freedom typed on a single line. One coordinate must be denoted "x", and one may be denoted "y". If no "y" coordinate is specified, a 1D plot is made.
• Change the bounds for the one or two modes which define the cut.
• For a 2D cut contour plots can generated for which one can specify the number of contours and the corresponding PES levels. Choose between a linear and a logarithmic spacing (see showpot below for details).
• For a 2D plot, 3D surfaces can also be produced. Toggle between surface, contour, or both.
• The plot can be printed to a postscript file, or a printer, or saved as a gnuplot file for further editing.
• If the system has more than one electronic state, choice can be made as to how many states, and which, will be plotted.

If the psi file contains a number of snapshots these will be displayed a shot at a time. In the default "step through" mode, press enter to change the frame. Using option 280 it is possible to change to "movie" mode, when the frames are shown at 1 second intervals. It may be useful to turn of the key (option 240) as the key (labelling contours) can cause the frames to have different sizes.

Using the "overlay" option it is possible to display a PES under the evolving wavepacket. First plot and save the PES to an xyz file (option 5). Change to plotting a wavepacket (using 10). Then turn on the overlay plots option with 400. Now with 410 the file name should be given containing the PES. The selected file name is shown. Unless the replot options are toggled (9 for the plot, 420 for the overlay) "plot to screen" (1) will now prompt for contours for both plots before plotting.

The plot data is written by showsys to the file tmpplot.pl, where it can be analyzed or plotted by alternative plotting software. See also the note on dengen84 below.

Warning: Plotting densities may take a large amount of computer time, if the dimension is large (f>3). Because of this it may be useful to reduce the number of plots by shortening the time-interval. Use menu point 30 to set the time-bounds appropriately. Alternatively (and simpler) one may use the options -n and/or -skip and/or -step to limit the number of wavefunctions read in and processed. E.g. showsys84 -n 11 will plot only the densities for the first 11 time-steps. If tout=1 this means from 0 to 10 fs. showsys84 -skip 50 -n 1 will plot the density no. 51 only, and showsys84 -skip 5 -n 2 -step 10 will plot the densities no. 6, 16 and 26.
Note that there is the routine dengen84 which generates 2D densities similar as showsys does, but, as it does not request interactive input, it can be run in background.

sumrstver

Purpose : Calculate the superposition of up to 8 wavefunctions
(using restart files). The superposition is written to the
restart file in the name directory.

Usage   : sumrst84 [-mnd -D -w -ver -h -?] file.inp  .

Options :
-mnd    : Make name directory.
-D name : 'name' denotes the directory where files are written to,
(name in inpf file ignored).
-w      : Existing data is overwritten.
-ver    : Version information about the program.
-h      : Print this help text.
-?      : Print this help text.

The program calculates the superposition of up to eight wavefunctions from restart files. The input wavefunctions must have the same primitive grid and mode combinations, but may differ in the number of SPF. The algorithm used is outlined as follows:
1. Scan input files and determine maximum number of SPFs.
2. Read first restart (psi2) and multiply with coefficient.
3. Read next restart (psi1) and multiply with coefficient.
4. Construct a SPF-basis from SPFs of psi1 and psi2 while skipping linear dependent SPFs.
5. Transform psi1 and psi2 to the new SPF-basis.
6. Add psi1 and psi2 (spsi).
7. Transform spsi to natural orbitals.
8. Reduce SPF basis of spsi if more files are left.
9. Copy spsi to psi2.
10. If more files are left then goto 3.
11. Check superposition and write new restart. Done.

In step 8, the number of SPF for a given mode m is reduced if the superposition is constructed from more then 2 input wavefunctions. The new number of SPF for a given mode is the maximum number of SPF for this mode found in the input files. That is, the largest numbers of SPF found in the input files determine the precision of the new wavefunction.
For two wavefunctions the superposition is exact.

The program requires an input file with a RUN-SECTION and an INIT_WF-SECTION.

 Keywords for sumrst RUN-SECTION Keyword Description name=S The name-directory to write the new restart file to normalize(=R|S) Norm of the new wavefunction. The new wavefunction will be normalized to unity if no argument is given. If a real number R is given, the wave function will be scaled to norm R. Alternatively the strings S='yes' or S='no' can be given where S='yes' will result in a normalization to unity, S='no' will leave the wavefunction unchanged (same as not giving the keyword normalize at all). overwrite Enable overwrite of output

 Keywords for sumrst INIT_WF-SECTION Keyword Description file=S Path to a directory with a restart file coeff=R(,R1) The coefficient to be multiplied to the wave function specified in the previous 'file' statement where R is the real part of the coefficient and R1 the imaginary part.

All keywords within the INIT_WF-SECTION are mandatory. Furthermore, a file statement must be followed by a coeff statement. An example input file would look like:

RUN-SECTION
normalize = 1
end-run-section

INIT_WF-SECTION
file  = state0
coeff = 1.0
file  = state1
coeff = 1.0
end-init_wf-section

end-input

showpotver [-i -d -v -n -u -vfit -vpot -diff -abs -ver -h -?]

Purpose  : plots potential and/or potential fit interactively
Usage    : showpot84 [-i -D -d -v -n -u -vfit -vpot -diff -abs
-evec -b -contr -ig -ver -h -?]

Options  :
-i NAME : The files are stored in the directory NAME
-D DIR  : The output is written to directory DIR. This is useful,
if one has no right to write to the name directory.
The directory DIR must exist.
The string DIR may be a relative or a full path-name.
-d FILE : The dvr is read from file FILE/dvr
rather than from ./dvr or NAME/dvr
-v FILE : The potential is read from file FILE
If FILE is not given, FILE=vpot is assumed.
If FILE=no, i.e -v no, vpot file is not read.
-n FILE : The potential fit is read from file FILE
If FILE is not given, FILE=natpot is assumed.
-u UNIT : The energy unit UNIT is applied.
-evec M : The eigenvectors (natural potentials)
of mode M are written to file evec.
M must not equal modc (contracted mode).
If M=99 then natpots of all modes (except modc) are written.
-b      : Eigenvectors (natural potentials) are additionally
written to binary file evecb, (if -evec is set).
-contr  : The contracted natural potentials of the
contracted mode are written to file contr.
WARNING: This may be a lot of data.
-vfit   : Set plottask to vfit.
-vpot   : Set plottask to vpot.
-diff   : Set plottask to vfit-vpot.
-abs    : Set plottask to abs(vfit-vpot).
-ig     : ignore certain errors.
-ver    : Version information about the program
-h      : Print this help text.
-?      : Print this help text.
------------------------------------------------------------------------------

NOTE
- If the -n option is not used,
the program assumes that a natpot file exists in the working
directory (or in the NAME directory if specified). If the -v option
is not used, the program looks in the same place for a vpot
file. If one is present, it can be used as a comparison for the fit.
- If you specify both input files, they MUST be based on the same dvr
and they SHOULD be based on the same PES parameters. If the latter
condition is not met a warning is issued. In this case the user
should have a look at the log file showpot.log where the potfit input
files of the vpot and the natpot file are stored.
- If you specify -v notice that showpot loads the complete vpot file into
a single precision vector held in memory. Hence, if one is dealing with
a huge product grid, the loading may take a while and showpot may
consume quite a bit of memory. Actually its memory consumption is of the
same order as potfit's memory consumption in the single-precision mode.
- The "-evec M" option lets showpot write the natural potentials
(eigenvectors of the potential density) of mode M to an ASCII file
evec (e.g. evec02, if M=2). If M=99 then all natural potentials
(except those of the contracted mode) are written to files.
- The "-contr" option lets showpot write the contracted natural
potentials of the contraced mode (i.e. the coefficients) to
file contraced.

Cuts of a potential energy surface (PES) can be generated from a vpot and a natpot file. Both files are generated by the potfit program. The first one contains the exact potential on the complete product grid, the latter a corresponding fit.

After starting the program a menu is presented by the help of which one can select a cut via an interactive process. The menu offers the following choices:

• Select WHAT is plotted: if both files, a vpot file and a natpot file, are read one can plot the PES or a fit of it or the (absolute) difference of both.
• Select between 1D and 2D cuts.
• Select new energy bounds for the plot.
• Change the bounds for the one or two modes which define the cut.
• For a 2D cut contour plots can be generated for which one can specify the number of contours and the corresponding PES levels. Choose between a linear and a logarithmic spacing (see below).

To define a cut, one must attribute x to one coordinate, and one may attribute y to another coordinate (the latter choice generates a 2D plot). The remaining coordinates must be assigned numbers (coordinate values), but at most one of these coordinates may be assigned with the label min or max. If min or max is given, the program sets this coordinate such that the potential is minimal (maximal)with respect to variations of this particular coordinate. max is useful for the plottask abs(Vfit-Vpot).

Description of the logarithmic spacing if Gnuplot commands are included into the output file:

• Firstly, the program calculates the difference d=zmax-zmin, where zmax (zmin) is the upper (lower) bound chosen for the contour levels. In a second step, the program calculates a basis b<1 such that bn-1=r, where n is the number of contour levels and r<1 is a so-called "reduction factor" provided by the user. Contour level i is now fixed at zmin+d*bi , where i runs from 0 to n-1. Note that according to this procedure zmax and/or zmin may be negative!

sumspecver

Purpose : Sums spectral components (output from autospec)
Usage   : sumspec84 [options] file1 fac1 file2 fac2 [file3 fac3 ...]

Options :
-i DIR  : paths are relative to directory DIR
-o FILE : write to file FILE. The default is: sumspec.pl
-w      : overwrite enabled.
-g ncos : GNUplot command lines are written to the output file.
ncos = 0,1,2 is the exponent of the cosine damping function.
NB: This sets "using 1:(ncos+2)"
-a      : GNUplot is called automatically. (sets "-g 1").
-G      : Plot gridlines. Requires -g or -a option.
-x xmin xmax : The energy range is set to [xmin,xmax].
(ignored if neither -g nor -a is set).
-y ymin ymax : The ordinate scale is set to [ymin,ymax].
(ignored if neither -g nor -a is set).
-s shift: Shift the energy scale by -shift
(ignored if neither -g nor -a is set).
-S shift: Shift the energy scale by -shift
Similar to -s, but this time the output file is modified.
-ver    : Version information about the program
-h      : Print this help text.
-?      : Print this help text.

Example: sumspec84 -a -G b2u/spectrum.pl 1.0  b2g/spectrum.pl 1.0

Following the options there follow the names of the files containing
the spectral data. For all components a factor, giving the weight,
------------------------------------------------------------------------------
The spectral component files are assumed to be generated by the autospec program, but sumspec will add the columns 2,3, and 4 of any sets of files. Hence, one may use sumspec to add flux files. In this case use the option "-g 2" to plot the probabilities.
The factors may be negative, and there may be only one input file.

twopdensver

Purpose: Computes the 2-particle density in the basis of the SPFs.
Usage:   twopdens [-f -i -skip -cut -n -s -r -ver -h -?] .

Options :
-f FILE : The wavefunction is read from file
FILE rather than from ./psi
The string FILE may be a relative or a full path-name.
-i DIR  : data is read from directory DIR
-n step : Compute the 2-particle-density only every step-th output.
-skip n : Skip the first n WFs when computing the 2-particle-density.
-cut n  : Consider only the first n WFs when computing the 2-p-density.
-s  s   : Compute the 2-particle-density only for state s.
-m1 mm1 : Set the first mode index to mm1.
-m2 mm2 : Set the second mode index to mm2.
-noev   : Do not print the eigenvalues.
-r      : Take the restart file rather than the psi file.
-ver    : Version information about the program
-h      : Print this help text.
-?      : Print this help text.

It is in general useful to pipe the output into less.
------------------------------------------------------------------------------
The program computes rho_{j,k;j',k'} = Sum_J   A*_{j,k,J} A_{j',k',J},   where J stands for all remaining indices. Of course, j an k do not need to be the first two indices as in this example. The two modes for which rho is computed are called m1 and m2. rho is interpreted as a matrix where (j,k) and (j',k') are taken as super-indices. This matrix is then diagonalized and its eigenvalues are printed to stadard output.

twprobver

Purpose: calculate reaction probability with the Tannor-Weeks
method (see JCP 110 (1999) 2761.)

Usage: twprob84 [ -o -F -e -p -t -cos -exp -exp2 -h -? ]
emin emax unit file [ edstr1 edstr2 ]

Options:
-o FILE : filename for output
-F x    : multiply output with a factor x
-e R U  : shift energy scale in output by value R in units U
-t tcut : read crosscorrelation function only up to T=tcut
-cos n  : multiply crosscorr.function with cos(pi/2*t/T))**n
-exp tau: multiply crosscorr.function with exp(-t/tau)
-exp2 tau:multiply crosscorr.function with exp(-(t/tau)**2)
-[h?]   : print this help text

"emin" and "emax" together with "unit" determine the energy
range of the output.
"file" contains the time-dependent overlap of the propagated
WF with a reference WF (output of crosscorr84).
"edstr[12]" are files that contain the energy distributions of
the two WFs. This can be an enerd file, or any other file, in
which case the distribution is determined by columns 1 and 2
(useful for autospec/cspec/flux files). If an energy
distribution is omitted, it is replaced by a constant of 1.0.

------------------------------------------------------------------------

This utility calculates state-resolved reaction probabilities according to the formula

from the correlation function $c\left(t\right)=⟨{\phi }_{f}|exp\left(-i\stackrel{^}{H}t\right)|{\phi }_{i}⟩=⟨{\phi }_{f}|\psi \left(t\right)⟩$ , where ${\phi }_{i}$ denotes the initial state of the propagated wavefunction, and ${\phi }_{f}$ denotes a reference wavefunction. Δi and Δf are the energy distributions of these two wavefunctions, respectively. Either energy distribution may be omitted, in which case it is formally replaced by a constant of 1.0.

The reference wavefunction ${\phi }_{f}$ must be a direct product of a (narrow) Gaussian in the translational DOF and an internal state of the products. This internal state is chosen according to the channel you are interested in.

The energy distributions are best obtained from enerd-files. However, for the propagated wavefunction it is also possible to calculate Δi from the autocorrelation function (using crosspecver).

The usage of twprobver is usually as follows:

1. Do a normal propagation with MCTDH. It's advisable to use the option "correction=dia" (INITWF-section) in order to create the enerd file. The propagation time T must be such that the overlap with the reference WF vanishes for t→T.
2. For each channel (i.e. each combination of quantum numbers for the final state) of interest, do:
1. Prepare an input file for a GENINWF run. The generated wavefunction must be a direct product of
• a Gaussian in the translational DOF. The Gaussian should be very narrow (this will make the energy distribution Δf broad and also make the correlation function c(t) disappear more quickly), and its center should be outside of the interaction region.
• an internal state of the products that reflects the chosen channel.
Here too, you should give the option "correction=dia" to produce an enerd file.
2. Run MCTDH with this input file to create the restart file, which will contain the reference wavefunction.
3. Compute the correlation function:
crosscorrver -o crossfile -f /path/to/propagated/wavefunction/psi -r /path/to/reference/wavefunction/restart
4. Compute the reaction probability for this channel:
twprobver -o probfile EMIN EMAX UNIT crossfile /path/to/propagated/wavefunction/enerd /path/to/reference/wavefunction/enerd

In practice, you will probably be interested in a lot of channels (note that with twprob there is no implicit summing like in flux), so you will want to write a script that generates the GENINWF input file, runs mctdh, crosscorr and twprob for each channel. It's advisable to make use of the readoper/readdvr options to speed up these calculations (the geninwf part of mctdh as well as crosscorr and twprob take virtually no time to complete).

The output file consists of two columns: energy E and probability P(E). E runs between the given EMIN and EMAX. Note that P(E) is set to zero for values of E that lie outside either of the two energy distributions.

Note: The distribution Δ(E) in the above formulas is called |Δ(E)|2 in section 8.6.3 of the MCTDH review.

vminmaxver

Purpose: Finds the "num" lowest and largest points on a PES.
Usage:   vminmax84 [ -i -p -o -n -s -d -ham -u -w -c -t -ver -h -?]

Options  :
-i DIR  : DIR is taken as  name-directory rather than ./
-p FILE : the PES is read from FILE rather than from ./pes
The string FILE may be a relative or a full path-name.
-o FILE : The ouput is written to file FILE rather than to name/vmm
The string FILE may be a relative or a full path-name.
If FILE=no, i.e. -o no, no output file will be written.
-n num  : Print num maximal and minimal values.
The default is num = 100 . Maximum is 250.
-s s    : Print maximal and minimal values for adiabatic state s. (Def: s=1).
-d s1 s2: Print maximal and minimal values for diabatic state (s1,s2).
-ham nhm: Use the Hamiltonian nhm. Default is nhm=1 (system).
-srf    : write full potential to ASCII file "srf.asc" (one energy
per line, for use as readsrf input).
-u UNIT : Use energy unit UNIT. (Default: eV).
-w      : Overwrite the output-file vmm.
-c srf  : Substract potential in file srf (one energy
per line, like readsrf input) before evaluation.
-t trj  : Do not run over the complete grid but use dvr indices
in file trj (n blank separated integers per line)
-ver    : Version information about the program
-h      : Print this help text.
-?      : Print this help text.

For a multi-state PES file the lowest adiabatic surface is analysed by default.
To analyse diabatic diagonal or off-diagonal surfaces use the option -d.
E.g. use "-d 1 1" to analyse the 1st diabatic state, "-d 1 2" refers to the coupling.
------------------------------------------------------------------------------
This program is very useful for checking if the potential of the implemented operator assumes unexpectedly large positive or negative values. The pes-file is most conveniently generated by using the MCTDH option -pes.

If Hamiltonians other than the system one are to be inspected (option -ham), then these operators should be generated with the usediag keyword.

The option -t followed by a filename (trj) allows to run over a subset of primitive grid points. The file contains the indices of the grid points of interest as bank separated integers in ASCII format.
In addition, the option -c allows comparison of the pes-file to the exact potential stored in a surface file (srf).

These two features are used to test an implemented cluster-expansion against the exact potential within the chkpes Python script. Here the indices represent the trajectory of a random walker and the surface file contains the exact potential energies at those points.

wfprojver

Purpose: Computes the projection of a WF onto one mode, i.e.: <psi|P|psi>,
where P= |phi><phi| denotes a projector onto one mode.

Usage  : wfproj84 [-f -i -o -dvr -n -lo -hi -P -s -w -ver -h -?]

Options :
-f FILE : The wavefunction is read from file FILE rather than from ./psi
The string FILE may be a relative or a full path-name.
-i FILE : The operator is read from file FILE rather than from ./oper
The string FILE may be a relative or a full path-name.
-dvr FL : The DVR information is read from file FL rather than from ./dvr
The string FL may be a relative or a full path-name.
-o FILE : The output is written to FILE rather than to ./wfproj .
The string FILE may be a relative or a full path-name.
-n step : Compute WFProj only every step-th psi-output.
-s nst  : Compute the WFProj for the electronic state nst only.
-lo klo : Compute WFProj only from the klo-th psi-output onwards.
-hi khi : Compute WFProj only up to the khi-th psi-output.
-P      : Use projectors. -P is follwed by the mode number, the projector
name and up to 8 parameters. The input is closed by a "%".
There may be several projectors. Example:
-P 3 KLeg 2 1 % -P 2 eigenf vib 6 %
Most useful is the projector "read".
See the documentation of flux84 for possible projectors.
-w      : A WFProj file may be overwritten.
-ver    : Version information about the program
-h      : Print this help text.
-?      : Print this help text.

Example:  wfproj84 -w -P 2 read ../inwf 1 1 %
------------------------------------------------------------------------------
This program uses the same projectors as flux84. See the documentation of flux to learn more about the projectors implemented. The most useful projector is possibly the read projector. Here the projector is build from an SPF which is read form a foreign restart file. Note that there may be more than one projector. If the MCTDH wavefunction is expanded in p modes (particles) then there may be at most p-1 projectors.

## Gnuplot scripts (pl-scripts)

There are a number of shell scripts located on \$MCTDH_DIR/bin, which enable a convenient plotting of some results of an MCTDH calculation. In general, these scripts call one of the above analyse routines, write the output to a temporary file and plot it using GNUPLOT. (A GNUPLOT of version 3.7 or higher should be used). The default version of the called analyse routines is ver=84. This may be changed by setting the environment variable \$MCTDH_VERSION or by using the option -v. The name directory mustbe the current directory when executing the pl-scripts. The pl-scripts allow for options and some scripts need arguments. A help text is always generated by executing   plxxx -h .

The different scripts support different options. If appropriate, however, the scripts commonly support the following options:

-h     : print a help text.
-a val : set lower range of x to val.
-x val : set upper range of x to val.
-y val : set upper range of y to val.
-z val : set lower range of y to val.
-G     : draw grid lines.
-l     : use logarithmic scale for y.
-s     : suppress messages of the called analyse program.
-v val : set the version number of the analyse program to be call
to val.
-p : prompt for printing the plot at default printer lpr. -P printer : specify alternative printer (e.g. -P "| lpr -Pps2"), or print to a file (e.g. -P plot.ps).
• plall
Prints an overview of all pl-scripts. This is the only pl-script that does not produce a plot. No arguments are required.
Reads the adpop density files and plots the one or two dimensional densities of the choosen electronic states.
Reads the adwkb file and plots the energy distribution DELTA(E) or the translational mean field. No call to any analyse routines. No arguments are required. (Obsolete, use plenerd).
• plenerd
Reads the enerd file and plots the energy distribution DELTA(E) or the translational mean field. No call to any analyse routines. No arguments are required.
• plauto
Plots the autocorrelation function. It reads the auto file and does not call any analyse routines. No arguments are required.
• plcap
Plots the reflection and transmission probabilites of a CAP. No output file of other programs needed.
• pledstr
Plot of the (free) energy distribution of a gaussian WP. No output file of other programs needed.
• pleigval
Plot of the energies computed in a diagonalisation run. It reads the ./eigval file generated by a mctdh-diagonalisation run. No arguments are required.
• plfdfour
Plots Fourier file created by the filter-diagonalisation program. Requires as argument the name of the file containing the Fourier spectrum generated by the filter84 program. No analyse routines are called.
• plfdspec
Plots line spectra created by the filter-diagonalisation program. Requires as argument the name of the file containing the eigen values. If no file name is given, the default filter.eig is taken. No analyse routines are called.
• plflux
Plots the quantum flux, the energy distribution of the initial wave packet and the reaction- (or scattering-) probability. It reads the flux file and does not call any analyse routines. However, to create the flux file, flux84 (see above) must have been run before excecuting plflux. No arguments are required.
• plgen
General GNUPLOT wrapper. Plots the data of ASCII file(s).   plgen requires as arguments a file name and the argument of the GNUPLOT using statement.
• plgpop
Plots the sum of the grid population of the first and last nz grid points. The occupation of the highest nz FBR basis functions is also shown. plgpop calls rdgpop84 (see above) to read the gridpop file. It requires the arguments nz and dof, where dof denotes the degree of freedom of which the grid populations are checked.
• plnat
Plots the natural weights. It calls rdcheck84 (see above) to read the check file. It requires the arguments state and mode, where state denotes the electronic state and mode the (combined or uncombined) mode for which the natural occupations are shown.
• plpit
Reads the file iteration from the potfit name directory and plots the potfit error measures. It requires the argument(s) task which can take the values: wr (wighted relevant), ur (unwighted relevant), ua (unweighted all), mr (maximum relevant), ma (maximum all), f1 (error measure for the first dof, see eq.(118) of the review), f2 (as before for the second dof),...
• plpweight
Reads the file prodwei from the potfit name directory and plots the separable weights. It requires as argument the number of the degree of freedom for which the weight is depicted.
• plqdq
Plots the expectation values <q> and <dq>. It calls rdcheck (see above) to read the check file. It requires the arguments state and dof, where state denotes the number of the electronic state and dof denotes the degree of freedom of which the expectation values are shown. Note: plqdq requires a check file which is generated by mctdh of version is 8.2.2 or higher.
• plrlx
Plots the energy vs time of an "improved relaxation" run. It reads the rlx_info file and calls plgen. The options passed to plgen must be quoted.
• plbrlx
Similar to plrlx, but for block-relaxation.
• plspec
Plots the (absorption) spectrum. It calls autospec84 (see above) to read the auto file and to compute the spectrum by fourier transform of the autocorrelation function. It requires as arguments the arguments of autospec.
• plspeed
Plots the CPU time used between two outputs. It reads the speed file and does not call any analyse routines. No arguments are required.
• plstate
Plots the electronic state populations. (If there is only one state, this is equivalent to the square of the norm). It calls rdcheck84 (see above) to read the check file. No arguments are required.
• plupdate
Plots the CMF update times or the CMF A and Phi errors. It reads the update file and does not call any analyse routines. No arguments are required.
• plwtt
Plots the CAP expectation values Wtt. It reads the wtt file and does not call any analyse routines. However, to create the wtt file, flux84 (see above) must have been run before executing plwtt. No arguments are required.

## Utility scripts

compile : See Installation and Compilation

cdm : A very convenient cd for the MCTDH directory. E.g cdm linear will cd to \$MCTDH_DIR/source/lib/linear. The destination directory may be abbreviated, as long as the abbreviation is unique. I.e. cdm lin is sufficient. The source dirs come before all others. I.e. cdm m cd-es to source/mctdh. To cd to doc/mctdh one types cdm doc/m. The cdm script is part of .mctdhrc and is available only, when working under bash.

elk_test : See Automatic Program Test

elk_test_gen : See Automatic Program Test

ddiff : compares files of current directory with those of another directory. (not recursive).

find235.py: finds integer numbers of the form 2x3y5z close to a given number -- these are valid grid sizes for the FFT in MCTDH.

mkGpatch : makes a GNU-patch. See Applying Patches

mkMpatch : makes a M-patch. See Applying Patches

mdistribute : Copies the files collected by mkMpatch or mfind to the MCTDH directory.

mdircp : Copies a MCTDH directory to a new location, excluding all those files which are not on the original package before installation (e.g. object, binary, *~, and *.dvi files).   mdircp can also create tar files and may send them to a remote host.

mfind : Finds all files of the MCTDH directory that are newer than a specified date (default 1 day).

mcb : MCTDH Code Browser

mcg : MCTDH Code Grep

mcl : MCTDH Code List

menv : MCTDH Environment

mhelp : MCTDH short help for input

minstall : Sets environment variables and makes another MCTDH directory active. This script is part of .mctdhrc and is available only, when working under bash. Typing " minstall path-of-second-MCTDH-directory " makes the second-MCTDH-directory active.

mbackup : A simple but convenient backup script for MCTDH developers. This backup system uses mdircp and mdiff as well. The command lb lists the backup directory. See MCTDH Backup

mdiff : Compares two MCTDH directories and lists files which differ. It can show the differences using diff or mgdiff.