The Analyse Programs

Generalities

Analyse is a set of programs able to extract data of interest from the results of an QUANTICS calculation. In essence, there is a library of modules providing an interface to the QUANTICS 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.

Some of these programs are seldom used, may not be very general, and should be used with caution. Some do not work as they were written for earlier versions of the package and need repairing.

Documentation of the Individual Programs

While all the analyse programs are stand alone codes, analysis is a basic interface to the main analysis programs. It uses a set of menus to interactively analyse the results of calculations.

The following programs are available for general analysis and plotting:

gwptraj: Plots the centres of GWP basis functions in G-MCTDH and vMCG calculations.
showd1d: Enables plotting of the 1-dimensional densities with time.
showfield: Enables plotting of a time-dependent field.
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.

The following programs are available to calculate and plot properties:

adpop: Calculates the diabatic and adiabatic populations of the electronic states. The calculation of one and two dimensional adiabatic densities is also possible.
autospec: Computation of the spectrum by fourier-transforming the autocorrelation function.
crosscorr: Reads two QUANTICS 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 from a wavefunction.
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.
flux: Computes the quantum flux by analysing the matrixelements of the CAP. Computes reaction probabilities.
gwptransmission: Transforms a GWP wavefunction to the energy representation and computes the energy and transmission spectrum.
jinpol: Interpolates flux-probabilities for a range of total J.
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 a file.
probrw: Computes reaction probability from Flux for real wavepacket propagation
probsq: Reads the auto file and computes the average of |c(t)|², which is the sum of the squared occupation probabilities of the eigenstates.
rdauto: Reads the auto steps from file "auto" and enables plotting of the data.
rdcheck: Outputs various simple expectation values.
rdeigval: Reads the eigenvalues from file "eigval" and enables plotting of the stick spectrum.
statepop: Outputs the diabatic state populations. If adiabatic populations are not needed this program is preferred to adpop as it is quicker and has interactive plotting.
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.

The following programs are available for checking the data used in a simulation:

enqoper: Queries the status and values of parameters and data stored in the oper file.
maxpe: Calculates the minimum and maximum value of a PES.
rdacoeff: Reads the A-vector and prints it to screen.
rdspfs: Reads the SPFs and prints them to a file.
rdgrid: Reads the DVR grid points from dvr file and prints them to screen.
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.

The following programs are available for checking the convergence of MCTDH simulations:

aoverlap: Calculates the difference between two A-vectors.
herma: Computes the hermiticity error of the A-vector.
lownp.py: Prints minimum over time of the populations for all listed nodes an modes of a ML-run.
norm: Calculates the norm and A-norm of a wavefunction or density operator of type II.
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.
rdautoe: Analyses the error within the autocorrelation function.
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.

The following programs are available for analysis of direct dynamics simulations:

ddgeom: Generates the expectation values for the geometries of a direct-dynamics simulation.
ddtraj: Generates the expectation values for the geometries from a direct-dynamics simulation.

The following programs are available for analysis and manipulation of direct dynamics database:

checkdb: Checks and cleans an SQL direct-dynamics database to remove failed and duplicate points.
cleandb: Cleans an old flat file direct-dynamics database to remove failed and duplicate points.
convdb: Converts an old flat file direct-dynamics database from binary to ascii or vice versa.
convdbfromsql: Converts a direct-dynamics database from SQL format to plain files.
convdbtosql: Converts a direct-dynamics database from plain files to an SQL format.
rddddb: Reads a DD DB and produces a file to plot the geometries.
trafodb: Analysis of DD DB using transformations of the energies and coordinates.
analyse_dd_traj: Analysis of DD Trajectories using range of different statistical techniques.

The following programs are available to manipulate and transform the wavefunction :

concat: Reads two restart files and concatenates the wavefunctions, i.e. psi(q,x)=psi1(q)*psi2(x).
convpsi: Converts a multi-packet wavefunction to the corresponding single-packet ones.
di2ad: Transforms a diabatic wavefunction to an adiabatic representation.
joinpsi: Joins psi files that have been calculated with different number of single particle functions.
ml2mctdh: Converts a multilayer wavefunction to MCTDH form.
oppsi: Applies an operator to a wavefunction.
psi2ex: Converts an QUANTICS 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 (see also rst2psi).
rst2psi: Converts a restart-file to psi-file format. (See also psi2rst).
rdrestart: Reads the info part of an restart file.
sumrst: Allows calculation of the superposition of up to 8 wavefunctions with different weight.
wfproj: Computes the projection of a wavefunction onto a (combined) mode.

The following programs are available for analysis of Filter Diagonalisation and Obtaining eigenvalues:

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.
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.
diag: Diagonalizes the Hamiltonian matrix calculated by fmat.

The following programs are available for set-up analysis of Optimal Control Calculations :

efield: Calculates the electric field for optimal control runs. efield exchanges data with quantics via fifos. efield is usually called by the script optcntrl.
guessfld: Computes the initial guess field for an optimal control run.
xfrog: Computes an XFROG trace of the electrical field obtained with OCT-MCTDH

The following program is available for wavefunction analysis:

diffmap: Analyses grid-based wavefunctions using diffusion map or isomaps.

Other specialists programs are also available:

adproj: Calculates vpot files for the adiabatic potential matrix elements and the adiabatic projector matrix elements for each electronic state.
compare: Used by the elk_test program to compare different sets of results.
mmemtest: Checks how much memory is available.
reflex: Computes reflection and transmission of a CAP.
tdpespec: Calculates a pump-probe photo-electron spectrum from a set of calculations.

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

Options

The adpop program is used to calculate the populations of the electronic states from the wavefunction stored in the psi file. Both the diabatic and the adiabatic populations are given. The results are stored in the file adp for each time step. Various options are availale to control the output times so that, e.g. -step or -tcut. The adp file can be easily plotted with plgen.

The full calculation runs over the primitive grid and is slow. It is also restricted to systems for which the full primitive grid fits into memory, typically 4 or fewr DOF. The calculation can be accelerated and memory saved by setting the quick option, -q. Using the quick option, one ignores points for which the product of one-dimensional grid-populations (1D-densities) is lower than qtol. This "quick" algorithm is the default and the default for qtol is 10^-5. It is turned off (to provide the exact calculation) using the option -noq.

One and two dimensional densities can also be calculated with adpop when using the full grid. 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.

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

Larger systems can be treated using some approximate schemes. A TDDVR scheme using the DVR from the SPFs is available with the -tddvr option. For calculations with combined mode SPFs the multi-dimensional DVR of Harrevelt and Manthe (JCP (05) 123: 064106) is used. The CDVR correction can also be used with -cdvr to obtain a more accurate result. This brings the memory requirements down to the length of the A-vector, affordable for all systems, but poor results are obtained with small basis sets.

An alternative approach is to use Monte-Carlo integration, -mc. This is, of course, less accurate and for small systems (3D, say) it is much slower than the direct method. The accuracy is controlled by the usual MC error variance measure, and is controlled by the -mctol keyword, which by default has a value of 0.001, i.e. plotting accuracy. MC integration allows the treatment of systems which are otherwise too large to be analyzed with adpop, but converges extremely slowly for systems with more than 10 DOFs due to the large empty space being sampled.

For single-set vMCG calculations, the -mc option leads to MC integration over each matrix element ⟨g_i | P | g_j⟩; rather than integrating the population expectation value directly. This allows the use of Gaussian sampling and is much more efficient, giving good results for large systems. By default matrix elements with a low density A_i^*⟨g_i | g_j⟩A_j are ignored. The tolerance is controlled by the qtol keyword, and -noq means all matrix elements are integrated. For direct dynamics calculations, by default local DBs are used with only points close to the GWP on the RHS included. The -noq option means that Shepard interpolation is done over all points in the DB for all points in the MC integrations.

For vMCG calculations, a quick and approximate calculation of the adiabatic populations can be made using a saddle-point approximation with -saddle. This evaluates the expectation value using the transformation only at the centre of all the GWPs multiplied by the overlaps. An even quicker (and less accurate) version includes only the diagonal matrix elements A_i^*⟨g_i | g_i⟩A_iP_α. This is obtained using the -saddle0 option.

There is an OpenMP version adpop.omp that speeds up CDVR, MC and saddle point calculations. The number of threads, N, is selected with -omp N.

Recommendations . For small grids use the default. For large grids, CDVR is the best option, but run the TDDVR calculation too. If these are similar the CDVR result is definitely to be trusted. For less than 10D systems MC integration can also be used as a check. For large single-set vMCG and DD-vMCG calculations use MC integration, with saddle-point as a quick check. With the MC integration f it runs quickly enough, try to take qtol to 0.0001 for a better result. If it is taking too long, try a qtol of 0.01.

adproj

Options

The adproj86 program is used to calculate the adiabatic potential matrix elements and the adiabatic projector matrix elements for each electronic state. The program must be run in a directory containing a dvr-file and a pes-file (usually the directory of the genpes-run). The results are written into several vpot files that can be fitted with the potfit program. To do this the keyword "pes = none" has to be set in the potfit input file. With a MCTDH-run the projection operator file can be created. To receive the adiabatic populations run expect with the projection operator file and the corresponding psi file.
NB: The vpot files for the adiabatic potential are: apr_vs, where s is the electronic state.
NB: The vpot files for the projection matrices are: apr_ps_ij, where s is the adiabatic electronic state and ij denotes the matrix element, i.e. the diabatic states. Please note that the projection matrices are symmetric, this means that each off-diagonal matrix element is only given once.
NB: When the option -od is set, the program computes the ij-th off-diagonal element of the electronic density matrix and nothing else. The vpot files produced are called apr_odij_kl, where k and l run from 1 to nstate. Note that these projection matrices are not symmetric, i.e. apr_odij_kl .ne. apr_odij_lk .
For more details see the MCTDH-guide chapter 12.8.3

analysis

Options

Follow the menus, and various information will be plotted.

analyse_dd_traj

Options

Use ddtraj to prepare the list of wavepacket trajectories

aoverlap

Options

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

Options

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 QUANTICS 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 the option -lin is set, then the filter function cosⁿ(pi*t/2T) (n=0,1,2) is replaced by cos³(pi*t/2T) (column 2), 1 - t/T (column 3), and (1-t/T)*cos(pi*t/T) + sin(pi*t/T)/pi (column 4) and cos²(pi*t/2T) (column 5). Note that the filters in colums 3 and 4 are non-negative. A more comprehensive discussion can be found in the lecture notes Introduction to MCTDH. Note, the filter of column 4 with -lin appears also on column 5 of the output created without -lin.

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. The energy given as argument to -e is added (the spectrum is shifted to the right). To substract a ground-state energy E_GS one must set -e -E_GS.

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

checkdb

Options

Check the properties of an SQL DB. It can also be used to remove points from failed calculations and duplicate points. These can be removed using a tolerance to define how far apart "duplicate" points can be. As well as the calculated data, essential files such as template.dat and the file containing the frequencies if normal modes are used may also be stored. These can be extracted using the -cre (create) option.

compare

Options

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

concat

Options

Concatenates the wavefunctions in two restart files.

cleandb

Options

Cleans an old style flat file DB to remove points from failed calculations. Duplicate points can also be removed using a tolerance to define how far apart "duplicate" points can be.

convdb

Options

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.

convdbfromsql

Options

Converts the database produced during a direct-dynamics calculation from an SQL database to a set of plain files. By default ascii files are produced (type .dba). The old DB files are deleted at the end.

convdbtosql

Options

Converts the database produced during a direct-dynamics calculation from an old set of plain files, either ascii or binary, to an SQL format. The old DB files are deleted at the end.

convpsi dir p

Not Available

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

Options

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 QUANTICS with the -t option.

crosspec

Options

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.

ddgeom

Options

ddgeom computes the expectation values for the geometries of a direct-dynamics calculation.

ddgeom can analyse distances, angles and dihedrals for the expectation values of the geometries of a direct-dynamics calculation.

ddtraj

Options

ddtraj computes the expectation values for the coordinates from a direct-dynamics calculation to generate trajectories. On start up a menu will be presented with options.

By default the molecular geometries are output in Angstrom. It is possible to select the trajectory of the centre of a single GWP by entering the GWP number. The average structure can also be generated using the "av" option, or the average structure for a single electronic state using "avS" and then entering the selected state. If Jmol is installed it will be launched automatically to animate the molecular structure along the "trajectory". Additionally the geometries are written to ddtraj_X.pl files, with the tag X being the GWP number, or "av" for the average etc.

If the option -notr is given and the propagation is in normal modes, then the values of the normal modes will be used and plotted using gnuplot. Again the trajectory for a single GWP or the averaged structures can be selected. The data is saved in ddtrajq_X.pl.

If -wr is selected, the files are written but no animation shown. The option "all" writes all of the trajectory files withpout animation. The files can all be animated using jmol (or gnuplot for ddtrajq files).

dengen

Options

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

Options

This program reads the wavefunction in a psi file, and uses the PES information in an oper file to transform it to the adiabatic representation. The waveunction and oper 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 oper file can be created in a genoper calculation 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

Not Available

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

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

diffmap

Options

Analyses an MCTDH or exact wavefunction in a psi file using diffusion map or isomap methods. Works by sampling gridpoints of the wavefunction at each selected time, locating points where the density is some proportion of the density at the centre of the wavefunction. Selected points are then analysed using eith diffusion map (default) or isomap in order to provide insight into the actual modes which most well describe the dynamics. For more details see Richings and Habershon, Molecules, 2021, 26, 7621.

Start the program in the directory containing the propagated wavefunction in a psi file.

Recommendation: The sampling can take long time as the program looks for gridpoints where the density is sufficiently large (as determined by the thresholds given in the options), so this is probably best for use when analysing relatively low dimensional problems.

efield

Not Available

edengen

Options

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

Options

The operator must first have been built by the QUANTICS 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

Options

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

Options

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

Options

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 must be of the form P=|psi><psi|, where psi is a mode-wavefunction.
Projectors operating on the CAP-mode are not allowed. When a operator operates on the CAP-mode, a symmetrization procedure has to be specified (-Os, see below). When projectors and operators are used simultaneously, then they must operate on different modes.

flux can work with operators. The operators are to be defined in a 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 number of 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.

Operating on the CAP-mode
Operating on the CAP-mode requires that the product F*O has to be symmetrized. F denotes the flux operator. The option -Os (which replaces -O) symmetrizes the product FO as 0.5*(OF+FO). Note that the computation time is almost doubled.

Finally, some remarks on setting the usediag or nodiag flag when generation the operators. Correlated operators must be generated with nodiag, when the are to be used with flux. Uncorrelated operators may use either flag, but usediag is to be preferred, in particular when the wavefunction is large. However, if the -Oq option is used, the operator (which actually is S) must be uncorrelated and generated with the usediag flag.

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 QUANTICS 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 (rather than from curent dir)	-f
outdir = S	write output to directory S	-o
make-outdir	Create the output directory.	-om
dvr = S	read dvr file form path S (rather than from curent dir)	-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 (rather than from curent dir)	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
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
symmetrized-operator	Apply the operator symmetrically, i.e. use 0.5(OF+F*O).	-Os
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.
p-1d	(none)	i1 i2	projects on a momentum state. (i.e. onto a plane wave) with momentum p = i1/i2. The momentum is in atomic units and is given as fraction because only integer parameters are allowed. A 1D SPF is assumed.
p-2d	(none)	i1 i2 i3	projects on a 2D momentum state. (i.e. onto a plane wave) with momenta p1 = i1/i3 and p2 = i2/i3. The momenta are in atomic units and are given as fractions because only integer parameters are allowed. A 2D SPF is assumed.
pop1	(none)	n	projects on the n-th grid-point of an underlying DVR. This projector can be used for combined modes as well, but then the user has to be careful to define n correctly. (e.g. n=n1+(n2-1)*gdim(f1), if m=(f1,f2))
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

Not Available

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

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.

QUANTICS: Analysis Programs

The Analyse Programs

Generalities

Documentation of the Individual Programs

adpop

adproj

analysis

analyse_dd_traj

aoverlap

autospec

checkdb

compare

concat

cleandb

convdb

convdbfromsql

convdbtosql

convpsi dir p

crosscorr

crosspec

ddgeom

ddtraj

dengen

di2ad

diag

Notes:

diffmap

efield

edengen

expect

fdcheck

fdmatch

flux

Example (H+D2 reactive scattering):

Example (Investigation of particle loss):

fmat

guessfld name

gwptraj

herma

gwptransmission

jinpol

joinpsi name1, name2, ...

enqoper

lownp.py

maxpe name

ml2mctdh name

ml2mctdh name

mmemtest

norm

oppsi operator

ortho

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

pexpect

prhosub

probsq

psi2ex

psi2rst

rst2psi

rdrestart

rdacoeff

rdspfs

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

reflex

rdspf

showd1d

showfield

showrst

showspf>

showsys

statepop

sumrst

sumspec

Example (H+D₂ reactive scattering):