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, 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 executable $MCTDH_DIR/bin/$MCTDH_PLATFORM/overlapd with debug information, type

compile -d overlap

For more information on the compile script, see Compiling the programs.

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

Documentation of the Individual Programs

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.
convdb: Converts a direct-dynamics database from binary to ascii or vice versa.
convpsi: Converts a multi-packet wavefunction to the corresponding single-packet ones.
crosscorr: Reads two MCTDH wavefunction (or density operators) ψ and ψ_ref and performs the cross-correlation function c(t) = <ψ_ref(0)|ψ(t)>. Matrix-elements <ψ_ref(0)|Ô|ψ(t)> can also be done (for WFs).
crosspec: Fourier-transforms a cross-correlation function generated by crosscorr.
dengen: Generates 2D-densities form a wavefunction.
di2ad: Transforms a diabatic wavefunction to an adiabatic representation.
diag: Diagonalizes the Hamiltonian matrix calculated by fmat.
efield: Calculates the electric field for optimal control runs. efield exchanges data with quantics via fifos. efield is usually called by the script optcntrl.
edengen: Generates (2D) pair-correlation functions, <x,x|ρ₂|y,y>, 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 <ψ(t1)|Ô|ψ(t2)>. The usual hermitian or the symmetric scalar product may be used. The diagonalisation of the matrix may then be performed by diag.
guessfld: Computes the initial guess field for an optimal control run.
herma: Computes the hermiticity error of the A-vector.
jinpol: Interpolates flux-probabilities 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.
norm: Calculates the norm and A-norm of a wavefunction or density operator of type II.
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 or density operators.
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)|², 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.
rddddb: Reads a DD DB and produces a file to plot the geometries.
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.
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.
xfrog: Computes an XFROG trace of the electrical field obtained with OCT-MCTDH
probrw: Computes reaction probability from Flux for real wavepacket propagation
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.

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

adpop


------------------------------------------------------------------------------
 Purpose: Calculates the diabatic and adiabatic population.
          And also one or two dimensional adiabatic densities.
 Usage:   adpop[-i -f -d -p -o -r -n -skip -step -tcut -w 
            -q -mc -noq -qtol -mctol -sleazy -tight
            -ver -h -?] [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. 

 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.
 EXAMPLES: adpop dof1 
           adpop dof3,dof4
           adpop dof1 dof2 dof3,dof4
------------------------------------------------------------------------------

The adpop 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.

Also one and two dimensional densities can be calculated with adpop. 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: the default calculation runs over the primitive grid and 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.

The adpop population file can be easily plotted with plgen.

For the adpop density files the script pladpop exists. This script plots the adiabatic densities of the different electronic states with or without weights.

adproj

 Purpose: Calculates the adiabatic potential and the adiabatic
          projector.
 Usage:   adproj [Options]

Options  :
 -i DIR  : data is stored in directory DIR
 -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.
 -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 adproj 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 electronic state and ij denotes the matrix element. Please note that the projection matrices are symmetric, this means that each off-diagonal matrix element is only given once.

analysis

Follow the menus, and various information will be plotted.

aoverlap

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

 Usage  : aoverlap [-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 aoverlap program is very similar to the overlap program, except that it considers the A-vectors exclusively, i.e. it implicitly assumes that the spf's of the two calculations are identical. aoverlap 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 overlap below.

autospec

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). Additionally there is the
          weight function exp(-(t/tau)^iexp).

Usage : autospec [-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
 -n       : The maxima of the three spectra are normalised 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  500
            Use "-p 1200" for final plots
 -q I     : Same as -p.
 -lin    : Use a linear filter, g(t) = 1 - t/T, for g0 and a
           combination of this and a cos or cos^2 filter for
           g1 and g2 (i.e. columns 2-4), respectively.
 -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.
 -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)×cosⁿ(π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:

σ (ω, n) = prefac \int_{0}^{T} Re [c (t) \exp (i (ω - E_{offset}) t)] \exp (- {(\frac{t}{tau})}^{iexp}) \cos^{n} (π \frac{t}{2 T}) d t

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 ω runs between E_min and E_max. The pre-factor prefac is 1/π if -FT is set (default), or is ω/π if -EP is set, or is 0.4279828×dipolemoment² if -Mb is set. In the latter case, dipolemoment is the norm of the D×Ψ, where D denotes the dipole operator. This norm is printed in the MCTDH log file, when the operation D×Ψ is performed by operate. The -Mb option allows to plot the absorption spectrum on absolute scale in mega barns (Mb). 1 Mb = 10^-18 cm²

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 E_min (and eventually E_max) may be negative. NB: When the pre-factor is omitted (default), the spectrum is normalized to 1 au.

Example:

autospec -e 13.4071 ev 13 17 ev 50 2

chnpot

  Purpose : Interpolation of natpot terms to another grid.
  Usage   : chnpot [-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. At present, the program can interpolate only one or two dimensional natpot terms. 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. 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 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. H₂O in Jacobian coordinates. If this periodicity is not given, e.g. H₂O 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:

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

compare

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

convdb

Interconverts the database produced during a direct-dynamics calculation between binary and ascii, i.e. if binary files are present (type .db) they are re-written as ascii (type .dba) and vice-versa. The old DB files are deleted at the end.

 Usage  : convdb [-ver]  

Options :
 -ver    : Version information about the program

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

crosscorr

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

 Usage  : crosscorr [-f -i -o -n -r -R -t -el -ig -ver -h -?]  [psi_ref] .

Options :
 -f FILE : The wavefunction (density operator)
           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.
 -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 (density  operator).

 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 (density operators) 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.

crosspec

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 : crosspec [-o -e -p -q -a -g -all -reim -re -im -sin -cos -dc  -t
                    -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 1200" 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 sin(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
 -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
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

\frac{1}{π} \int_{0}^{T} e^{i (E - E_{offset}) t} c (t) dt

without invoking damping functions (except possibly for sin(π t/T) and cos(π 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 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 flux. 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.

dengen

 Purpose: Generates 2D-densities from the wavefunction.
 Usage:   dengen [-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: "dengen -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 "dengen -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.
--------------------------------------------------------------------------------

dengen computes 2D-densities similar as does showsys. However, in contrast to showsys, dengen 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 dengen is of advantage, because dengen may run in the background, overnight or on a batch system.

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

di2ad

 Purpose: Converts a diabatic wavefunction to adiabatic.
 Usage:   di2ad [-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 psi2ex 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 quantics -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.

diag

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   : diag [-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  "fmat" program. If "-t" option is used,
 the second set of Overlap- and  Hamiltonian matrices with different
 scalar product is needed.

 Examples  : diag -es 1.d-8 -hi 60  mtt_h mtt_H1_h
           : diag -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. fmat -sym          ! calculates  mtt_s
        fmat -sym -O H1    ! calculates  mtt_H1_s
        diag -u ev -es 1.d-9 mtt_s mtt_H1_s

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

efield

 Purpose: Calculates the electric field for optimal control

 Usage  : efield [-b -t -A -T -O -F -p -i -ver -h -?]

Options: 
 -b       : Compute the electric field during a backward propagation.
 -t       : Compute overlap with target during a forward propagation.
 -A       : Align phases of overlaps of multiple initial and target states
            i.e. optimize |sum_i <Psi_i(T)|Psi_tar,i>|^2
            instead of    sum_i |<Psi_i(T)|Psi_tar,i>|^2
 -T oper  : Maximise matrix element <Psi|oper|Psi> of operator oper.
 -O oper  : Use oper as dipole operator to compute the field.
 -F pnlty : Parameter to penalize the field strength.
 -p       : Do not read old psi file. Instaed propagate 2 WFs simultaneously.
 -i iter  : Number of current iteration.
 -h or -? : Prints this help text.
 -ver     : Prints version information.

 The matrix element of the dipole operator (-O option) formed with the initial
 and target state is calculated and the E-field is evaluated.
 This routine MUST be started in the background before the invocation of MCTDH.
 It exchanges information with the running propagation via FIFOs.
 Note: This routine is usually  started by the script optcntrl.
------------------------------------------------------------------------------

edengen

 Purpose: Generates pair-density matrix <x,x|rho2|y,y> from the wavefunction.
 Usage:   edengen [-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: edengen -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 edengen -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.
--------------------------------------------------------------------------------

edengen computes 2D (off-diagonal) pair-correlation functions (or pair-density matrices) <x,x|ρ₂|y,y> (similar to the one-body density matrix produced by prhosub). It works in analogy to dengen, which however calculates the diagonal density <x,y|ρ₂|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, edengen may run in the background, overnight or on a batch system.

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

expect

 Purpose: Calculates the time evolution of an expectation value.
 Usage  : expect [-f -i -o -nofs -g -a -r -rst -x -t -p -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).
 -p      : Print GNUplot to lpr (sets -g).
 -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 process 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.

fdcheck

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

Usage  : fdcheck 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 fdmatch.
     Example: fdcheck 4 0 < list
     where list contains the output of fdmatch, or:
     fdmatch f1.eig f2.eig f3.eig | sort -n | fdcheck 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.

fdmatch

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

Usage  : fdmatch 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.

flux

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

Usage : flux [-f -i -o -ed -dvr -t -n -lo -hi -p -P -O -e -x -v -w -r -c -u
        -sm -es -wtt -nofs -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 ./adwkb
           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 (autospec)
           or a cspec.pl (crosspec) 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.
 -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: flux -lo 25 -w 0.3 2.5 ev rd
          flux 1000 3000 cm-1  x left
 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 flux 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.

flux 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 flux 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=|ψ><ψ|, which in turn is read from a (foreign) restart file. The arguments needed by flux 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 flux. 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, flux 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² 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×Ψ||, 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 flux 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 incomming 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) = < Ψ | Δ(H-E) | Ψ >

flux(E) = 2 π < Ψ |Δ(H-E) F Δ(H-E) | Ψ >

where F denotes the flux operator. DELTA(E) can be obtained trough several ways. quanticsh 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). flux will automatically look for such a file. Alternatively (by using the -ed option) one may use a spectrum.pl file (generated by autospec), or a cspec.pl file (generated by crosspec), or a flux file (generated by flux without employing projectors) as energy-distribution file.

Example (H+D₂ reactive scattering):

flux -e lsth -lo 7 0.25 2.75 ev rv

flux 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 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} ⟨ ψ (t') | W | ψ (t') ⟩ d t' + ⟨ ψ (t) | Θ (t - s) | ψ (t) ⟩

where Θ 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):

flux -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
smooth	Smooth the Delta(E) data read from enerd file by averaging over 5 energy points.	-sm
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 Eopnorm*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
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
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
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 `flux -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
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,m_j) 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,m_j) 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,m_j) rotational state. A PLeg-DVR combined with an exp-DVR, i.e. a 2D spf, is required.
KLeg2	(none)	j₁ m₁ j₂ m₂ sym	projects on the (j₁,m₁;j₂,m₂) 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. {\|j₁,m₁,j₂,m₂> ± \|j₂,m₂,j₁,m₁>} / √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.

fmat

  Purpose: Computes the overlap- or the Hamiltonian-matrix.
  Usage  : fmat [-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 : fmat -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).

guessfld name

 Purpose: Calculates the initial guess field for optimal control
          The pulse is a Gaussian times sin(omega_field*t)

 Usage  : guessfld tfinal tout omega_rabi omega_f [frequency unit]

 The times tfinal, tout, are in fs. The field-strength,
 omega_rabi, and the  field-frequency, omega_f, 
 are in eV, au or cm-1. The frequency unit can be
 given with an additional optional string. The default unit is au. 
------------------------------------------------------------------------------

The field is defined by:

E(t) = ω_rabi × cos(ω_f(t-t_final/2)) + *sin( πt/t_final )²

The field will be written in ASCII to the file fieldf.0 for every tout over the interval -2×tout ... t_final + 2×tout.

herma name

 Purpose: Checks the hermiticity of the A-vector.
          Only available for density type II.
 Usage:   herma [-f -i -n -r -ver -h -?] .

Options  :
 -f FILE : The density matrix 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 error 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.
------------------------------------------------------------------------------

jinpol

 Purpose : J-interpolation of flux data.
 Usage   : jinpol [-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: jinpol -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.
------------------------------------------------------------------------------

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

joinpsi name1, name2, ...

 Purpose : Join psi files have been calculated with different number of SPFs.
 Usage   : joinpsi [-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     : joinpsi ../dir1/psi ../dir2/psi ../dir3/psi 
                joinpsi -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.

maxpe name

 Purpose : Calculates the minimum and maximum value of a PES
 Usage   : maxpe [-h -ver ] path-to-name-directory
 First generate an oper file for an exact calculation and then
 run maxpe 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 maxpe. Output is to the screen.

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

mmemtest

 Purpose : Checks, how much memory can be allocated.
 Usage   : mmemtest [-h ]
 Options : -h    Print this help text. 
------------------------------------------------------------------

 Run mmemtest on the computer to see how much memory is available.

norm

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

Options : 
 -f FILE : The wavefunction (density operator) 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. 
------------------------------------------------------------------------------

npotminmax

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

Usage    : npotminmax [-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 srf (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
           like readsrf input).
 -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.	-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.	-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:


 NATPOT-SECTION
  ../POTFIT/pf_Q1_Q2   # directory with natpot of Q1 and Q2
  ../POTFIT/pf_Q2_Q3   # directory with natpot of Q2 and Q3
  ../POTFIT/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

ortho

 Purpose: Checks the orthonormality of the spf's (spdo's).
          For density operators it checks also the
          hermiticity of the one-particle density operators.
 Usage:   ortho [-f -i -n -r -ver -h -?] .

Options : 
 -f FILE : The wavefunction (density operator) 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.
 -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 the sum of the squares of the matrix-elements of (ov-1) where ov denotes the overlap-matrix of the spf's.

overlap [-f -i -o -n -w -h -?] comp_file

 Purpose: Calculates the difference between two wavefunctions
          or density operators. See review Eqs.(159-163) for a
          definition of the output.


 Usage  : overlap [-f -i -o -a -b -c -t -n -w -ig -ver -h -?]  [Psi] .

Options : 
 -f FILE : The wavefunction (density operator) 
           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.
 -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 or density operator.
 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.
 If one calculates the overlap between density operators of different types,
 [Psi] has to be the one of type II.
 NB: The overlap is written to the file ./ovl_XXX, the screen-output
 shows error=||Psi-Psi1||.
------------------------------------------------------------------------------

Calculates the error (overlap) between the wavefunctions (density operators) 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.

overlap ../dir/psi

calculates the error between the wavefunctions (density operators) 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

overlap -i psi1 psi2

would calculate the error between the wavefunctions (density operators) 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, overlap will read the file ovl_dir rather than re-do the calculation.

pexpect

 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.

prhosub

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

probsq

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

Usage : probsq [-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

\frac{1}{T} \int_{0}^{T} g (t) {| c (t) |}^{2} d t

where T denotes the final time of the autocorrelation c(t) and g(t) is some weighting function

g (t) = {\begin{matrix} \begin{matrix} \begin{matrix} 1 \\ \frac{π}{2} \sin (π \frac{t}{T}) \end{matrix} \\ \frac{π}{2} \cos (π \frac{t}{2 T}) \end{matrix} \\ 2 \sin^{2} (π \frac{t}{T}) \end{matrix}

For large T this converges to Σ p_i , where p_i denotes the occupation probability of the i-th eigenstate. This sum is related to the so called quadratic entropy. Output are the values Σ p_i, 1/(Σ p_i), and -Log₂(Σ p_i).

psi2ex

 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  : psi2ex [-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.
           Note: -abs sets -plot.
 -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).

psi2rst

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

 Usage  : psi2rst [-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).
------------------------------------------------------------------------------

rdacoeff

 Analyse program "rdacoeff" (Read A-coefficients)
------------------------------------------------------------------------------
 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: rdacoeff | less
or to re-direct it into a file:  rdacoeff > acoeff
------------------------------------------------------------------------------

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

 Analyse program "rdautoe" (Read autoe-file)
 Purpose : Analyses the error measures of the autocorrelation.
 Usage   : rdautoe [-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 autospec. 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.

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

 Analyse program "rdcheck" (Read check file)
 Purpose : Reads the check-file.
           Analyses the state population and natural weights.
 Usage   : rdcheck [-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: rdcheck 0 , or, to see the development of the lowest natural weights in timesteps of 1.2 fs:

rdcheck -t 1.2 0

The entropies of the natural populations are written to screen with the command rdcheck -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 = Σ_i -p_i Log(p_i)

tr2 = Σ_i p_i²

where p_i denotes the normalized (by norm²) i-th natural weight.

rddddb

Reads direct dynamics DB file geo.db (or geo.dba if ascii DB) and writes the geometries and energies into a formatted file plot.db . This is in a pseudo-Gaussian format and the geometries can be plotted by e.g. the Jmol program. There are no arguments

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

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:

rdgpop -s -oc 2 > file plgen file

Or plotting weights vs coordinate:

rdgpop -s -w 2 > file plgen file 2:3

rdgpop nz dof

 Analyse program "rdgpop" (Read grid populations)
 Purpose : Analyses the grid population read from file gridpop.
 Usage   : rdgpop [-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 rdgpop 1 0 while being in the name directory.

rdsteps

 Purpose : Reads the step-file and evaluates statistic measures
           of the integrator steps.
 Usage   : rdsteps [-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. 
------------------------------------------------------------------------------

rdupdate

 Purpose : Reads the step-file and evaluates statistic measures
           of the integrator steps.
 Usage   : rdsteps [-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. 
$ rdupdate -h
ld.so: Incorrectly built binary which accesses errno or h_errno directly.
ld.so: See /usr/share/doc/libc6/FAQ.gz.
 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.

rdauto

 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.

rdeigval

 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.

showd1d [-T|-S|-I|-P|-E|-M|-Ma] [-f -i -o -a -g -p -x -y -fac -n -no -t -nofs -l -w -sm -nw -pop2 -G -inter -ver -h -?] fx sx

 Purpose: To selectively plot 1-dimensional densitites.
 Usage  : showd1d [-T|-S|-I|-P|-M|-Ma] [-f -i -o -a -g -p -x -y -fac -n -no
                   -t -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 wavefunction 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.
 -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.
 -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 abszisse length (default: auto). Ignored for -T, -Ma.
 -y ymax : Set the ordinate length (default: auto). Ignored for -T, -Ma.
 -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).
 If no argument is given, the user is prompted for the missing argument.
 EXAMPLES: showd1d -a f1
           showd1d -a -M -y 0.3  -n 2  f2 s1 
           showd1d -w -p -I 20 f1
           showd1d -a -p -T f3
------------------------------------------------------------------------------

For instance, showd1d -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.

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.

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

showrst [-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  : showrst [-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. showrst -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.

showspf> [-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  : showspf [-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. showspf -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.

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

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

Options :
 -i DIR  : data is stored in 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.
 -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).
 -pop2   : the basis set occupations are plotted rather than grid
           populations. For FFT this implies momentum space representation.
 -ver    : Version information about the program.
 -h      : Print this help text.
 -?      : Print this help text.

  The program is menu driven. Just type "showsys" or "showsys -pop2"
  or "showsys -pes" or "showsys -rst" and then follow the menu.
  If there are more WFs on psi file than you want to inspect, use the
  options -n, -skip, or -step. Example: "showsys -skip 10 -n 2" will plot
  the WFs no. 11 and 12 (i.e. t=10 and 11 if tpsi=1fs),
  and "showsys -step 10 -n 2" will plot the WFs no. 1 and 11.
  If one wants to plot a potential one should use "shwsys -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
       (menu-point 30).
  NB2: The option -pop2 is 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 quantics with the option -pes.
Move to the directory where the pes file is, and start showpes -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.

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. showsys -n 11 will plot only the densities for the first 11 time-steps. If tout=1 this means from 0 to 10 fs. showsys -skip 50 -n 1 will plot the density no. 51 only, and showsys -skip 5 -n 2 -step 10 will plot the densities no. 6, 16 and 26.

Note that there is the routine dengen which generates 2D densities similar as showsys does.

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

Purpose  : plots potential and/or potential fit interactively
Usage    : showpot [-i -d -v -n -u -vfit -vpot -diff -abs -ver -h -?]

Options  :
 -i NAME : The files are stored in the directory 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.
 -vfit   : Set plottask to vfit.
 -vpot   : Set plottask to vpot.
 -diff   : Set plottask to vfit-vpot.
 -abs    : Set plottask to abs(vfit-vpot).
 -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 direct=1_sp mode.

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 b^n-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*bⁱ , where i runs from 0 to n-1. Note that according to this procedure zmax and/or zmin may be negative!

sumrst

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

Scan input files and determine maximum number of SPFs.
Read first restart (psi2) and multiply with coefficient.
Read next restart (psi1) and multiply with coefficient.
Construct a SPF-basis from SPFs of psi1 and psi2 while skipping linear dependent SPFs.
Transform psi1 and psi2 to the new SPF-basis.
Add psi1 and psi2 (spsi).
Transform spsi to natural orbitals.
Reduce SPF basis of spsi if more files are left.
Copy spsi to psi2.
If more files are left then goto 3.
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
  name      = HADplus
  normalize = 1
end-run-section

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

end-input

sumspec

 Purpose : Sums spectral components (output from autospec) 
 Usage   : sumspec [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").
 -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: sumspec -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,
has to follow the file-name.
------------------------------------------------------------------------------

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.

twopdens

 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 ρ_{j,k;j',k'} = Σ_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 ρ is computed are called m1 and m2. ρ 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.

twprob

Purpose: calculate reaction probability with the Tannor-Weeks
         method (see e.g. JCP 110 (1999) 2761, Section IIB.)

Usage: twprob [ -o -F -e -p -t -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
 -p I   : number of energy points is I (default: 500)
 -t tcut: read crosscorrelation function only up to t=tcut
 -[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 crosscorr).
"edstr[12]" are files that contain the energy distributions of
the two WFs. This can be an adwkb file, or any other file, in
which case the distribution is determined by columns 1 and 2
(useful for autospec/cspec/flux files).

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

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

P_{fi} (E) = \frac{1}{4 π^{2} Δ_{i} (E) Δ_{f} (E)} {| \int_{0}^{T} e^{i E t} c (t) dt |}^{2}

from the correlation function c(t)=<φ_f|e^-iĤt|φ_i>=<φ_f|ψ(t)>, where φ_i denotes the initial state of the propagated wavefunction, and φ_f denotes a reference wavefunction. Δ_i and Δ_f are the energy distributions of these two wavefunctions, respectively.

The reference wavefunction φ_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 adwkb-files. However, for the propagated wavefunction it is also possible to calculate Δ_i from the autocorrelation function (using crosspec).

The usage of twprob is usually as follows:

Do a normal propagation with MCTDH. It's advisable to use the option "correction=dia" (INITWF-section) in order to create the adwkb file. The propagation time T must be such that the overlap with the reference WF vanishes for t→T.
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 adwkb file.
2. Run MCTDH with this input file to create the restart file, which will contain the reference wavefunction.
3. Compute the correlation function: crosscorr -o crossfile -f /path/to/propagated/wavefunction/psi -r /path/to/reference/wavefunction/restart
4. Compute the reaction probability for this channel: twprob -o probfile EMIN EMAX UNIT crossfile /path/to/propagated/wavefunction/adwkb /path/to/reference/wavefunction/adwkb

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 quantics, 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 quantics 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)|² in section 8.6.3 of the MCTDH review.

vminmax

 Purpose: Finds the "num" lowest and largest points on a PES.
 Usage:   vminmax [ -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.
------------------------------------------------------------------------------

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.

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.

wfproj

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  : wfproj [-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 flux 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:  wfproj -w -P 2 read ../inwf 1 1 %

This program uses the same projectors as flux. 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.

xfrog

Purpose : Computation of the XFROG trace of a 1-dim function, using a
          rectrangular gate function with Gaussian shoulders.

Usage   : xfrog [-o -n -rf -rt -fd -td -g -ver -h -?] infile

Options :
 -o FILE     : The output is written to file FILE rather than to infile.xfg
               The string FILE can be a relative or a full path-name.
 -n nfield   : Number of points to be read from infile (default: all)
 -rf nres    : Sampling of outpot in frequency domain (default: nfield+1)
 -rt nres    : Sampling of outpot in time domain (default: nfield)
 -fd min max : Frequency range of the outtput (default: 0 to nfield/T)
 -td min max : Time range of the outtput (default: 0 to T)
 -g w1 w2    : Width of the gatefunction in time domain (default: w1=w2=1):
               w1 width of rectangular part, w2 width of Gaussian shoulders.
 -ver        : Version info of this program.
 -h, -?      : Print this help text.

The argument infile is the file containing the electric field. The
gate function is rectangular pulse of unity height and width w1 with Gaussian
shoulders of the form exp(-t**2/(2*w2**2))
The user will  be prompted for missing arguments. All input is in au.
------------------------------------------------------------------------------

probrw

Purpose : Computation of the reaction probability when using
          real wavepackete propagation.

Usage   : probrw [-w -d -ver -h -?] 

Options :
 -w      : The (un)prob_{dof} file may be overwritten.
-d dof   : Probability along the dof.
 -ver    : Version information about the program
 -h      : Print this help text.
 -?      : Print this help text.

 Example: probrw100 -d rd
------------------------------------------------------------------------------

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 analyse routines is called. 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 called 
          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.
pladpop: Reads the adpop density files and plots the one or two dimensional densities of the choosen electronic states.
pladwkb: 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.
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 filter 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, flux (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 rdgpop (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 rdcheck (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.
plspec: Plots the (absorption) spectrum. It calls autospec (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 rdcheck (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, flux (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/quantics one types cdm doc/m. The cdm script is part of .quanticsrc 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 2^x3^y5^z 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 .quanticsrc 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.
mklinks: Creates links from $MCTDH_DIR/source/surfaces/ to addsurf
phelp: Programmers help. Gives a description of MCTDH variables.
vddiff: Version doc diff. Compares the documentation of two MCTDH directories.
vrdiff: Version rest diff. Compares the 'rest' of two MCTDH directories.
vsdiff: Version source diff. Compares the source of two MCTDH directories.

Try program -h to obtain more information (-h option is not available for cdm). A MCTDH user should be familiar with compile, cdm, menv, minstall, mhelp, mklinks, mdistribute, mkGpatch, and mkMpatch. The other scripts are primarily for MCTDH developers. Of particular importance for developers are the scripts elk_test, elk_test_gen, mcb, mcg, and phelp.

MCTDH: Analysis Programs

The Analyse Programs

Generalities

Documentation of the Individual Programs

adpop

adproj

analysis

aoverlap

autospec

chnpot

compare

convdb

convpsi dir p

crosscorr

crosspec

dengen

di2ad

diag

Notes:

efield

edengen

expect

fdcheck

fdmatch

flux

Example (H+D2 reactive scattering):

Example (Investigation of particle loss):

fmat

guessfld name

herma name

jinpol

joinpsi name1, name2, ...

maxpe name

mmemtest

norm

npotminmax

Output files (in name directory)

ortho

overlap [-f -i -o -n -w -h -?] comp_file

pexpect

prhosub

probsq

psi2ex

psi2rst

rdacoeff

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

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

rddddb

rdgrid dof

rdgpop nz dof

rdsteps

rdupdate

rdauto

rdeigval

showd1d [-T|-S|-I|-P|-E|-M|-Ma] [-f -i -o -a -g -p -x -y -fac -n -no -t -nofs -l -w -sm -nw -pop2 -G -inter -ver -h -?] fx sx

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

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

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

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

sumrst

sumspec

twopdens

twprob

Note

vminmax

wfproj

xfrog

probrw

Gnuplot scripts (pl-scripts)

Utility scripts

Example (H+D₂ reactive scattering):