# Running Direct Dynamics Calculations using Grid-Based Dynamics

This section describes the procedure used to carry out direct-dynamics
using grid-based dynamics methods i.e. the standard method or MCTDH. This
method uses expansions of the PES in terms of Gaussian functions and is
based on the Gausssian Process Regression (GPR) procedure. There are
different,
but related, versions of this fitting procedure which have been
implemented. In the following description we assume that the dynamics is
proceeding in `F` degrees-of-freedom (DOFs) (in the current
implementation, these DOFs are mass-frequency-scaled normal modes).

All of the fitting procedures use a database of electronic energies
calculated at a set of randomly selected points in the F-dimensional
space defined by the choice of DOFs. The energies are calculated
on-the-fly by converting each point in the space to molecular Cartesian
coordinates, then passing that geometry to an electronic structure
package. For single-state calculations, only the electronic energies are
required, but for multi-state calculations, it may also be necessary to
calculate non-adiabatic coupling terms between all states. The coupling
terms are used in the propagation diabatisation procedure to diabatise
the energies, giving smooth PESs on which the dynamics proceeds. NB The NACTs are not
needed for the Procrustes diabatisation. This
sampling is carried out at pre-selected times in the dynamics; the
default is for the update to occur at the time of an output
(i.e. at `tout`, but this can be altered to integer multiples of
the output time to reduce frequency if desired.

The sampled points are chosen randomly within a range centred on the
wavepacket centre (given by the position expectation in each DOF). The
coordinates along each DOF are randomly chosen within a given number of
standard deviations of the wavepacket centre (the SD is given by the
wavepacket width along each DOF). The default is +/-3SD, but this can be
altered in the `dirdyn-gb-section` of the input. Not all points
are accepted, because, for each choice of geometry, the GPR can measure
the variance of the fit there. If the variance is below some threshold
(default is 1.0d-6, but again this can be altered), then the fit of the
PES is considered adequate based on the points already in the database, and
the point is rejected. If the variance is above the threshold, the point
is accepted, the energy calculated and then added to the database. It is also
possible to set an energy threshold above which points are rejected (ground-state
energy), which can be particularly useful with Morse type potentials with a step
potential wall. Once
the selection of points is complete, there are two possibilities: first,
if no extra points have been added, the PES fit from the previous step is
used and the dynamics proceeds. The second possibility is that some extra
points have been added to the database, in which case a new GPR fit of
the potential needs to be created. The following paragraphs describe the
different ways in which this can be done.

The basic implementation of the fitting is to expand the PES in terms of a set of N F-dimensional Gaussian functions centred at the randomly selected points in the F-dimensional space (and hence corresponding to the N different molecular geometries stored in the database). The weights of these Gaussians are determined by the procedure outlined in the papers referenced in the Citations. This method is preferred for exact dynamics using 2 or 3 DOFs.

An alternative expansion is through the use of an additive kernel of Gaussian functions. This kernel is simply a sum of F 1-D Gaussians with an optional (and strongly recommended) sum of F(F-1)/2 2-D Gaussians (i.e. all possible combinations of the F DOFs), and even a further optional (and less necessary except in cases with strong 3-mode coupling as the extra computational cost is significant), F(F-1)(F-2)/6 3-D Gaussians. Each of these sums is assigned a weight by the same procedure as outlined in the previous paragraph. For high-dimensional problems, this form of the expansion should result in fewer points needing to be added to the database of energies, than when F-D Gaussians are used, but there is an extra cost in that each of the N points now contributes F+F(F-1)/2 terms (plus F(F-1)(F-2)/6 if 3-D terms are included) to the Hamiltonian rather than just 1.

In order to overcome the slow propagation, with MCTDH, when using either of the above expansions of the PES for more than a few DOFs, a method for approximating the Gaussian expansion of the PES has been implemented, which has the goal of significantly reducing the number of potential terms in the Hamiltonian by fitting a series of 1-, 2- and 3-D terms to an underlying GPR expansion (of either form outlined in the above two paragraphs).

The secondary fitting procedure attempts to represent the GPR PES in terms of a small number of 1-, 2- and, optionally, 3-D functions, by disposing of the use of Gaussians, instead simply assigning values of the terms at the DVR gridpoints. The 1-D terms are simply the values of the GPR potential at the DVR gridpoints along the relevant DOFs (the first DOF includes the value of the PES at the origin). The fit of the 2-D terms relies on minimising the difference between the fitted values and the residual of the PES i.e. the GPR PES minus the relevant 1-D terms, by using Singular Value Decomposition to get the fit. By evaluating the residual at the points in configuration space, (qf,qg) where {qf} and {qg} are the positions of all DVR gridpoints along the DOFs, f and g, respectively, we have a 2-D array of values. This array is decomposed into a weighted sum of vector outer products by the SVD procedure. Each vector represents the values of the potential term at the DVR gridpoints along each DOF. This can be taken further by using 3-D tensor decomposition to get vectors along 3 DOFs to represent 3-mode couplings. Including 3-D terms should done with caution as the extra accuracy obtained is likely offset by the much increased effort required to propagate the wavefunction. However, it may be necessary in some cases e.g. malonaldehyde.

As the method is very much a combination of the concepts of the original MCTDH and standard methods, and the direct-dynamics techniques using the vMCG method, the structure of the input file needed to run such a calculation is essentially a mixture of the elements in the input files used in those methods. Therefore it is recommended that one reads the sections describing MCTDH and DD-vMCG inputs prior to reading this section. In this section, the similarities will only be summarised, so, for more detail, please refer to those sections.

The DD-GB method relies on the construction and propagation of a wavefunction
by expansion upon a "grid" of DVR (or similar) basis functions in an
identical way to the standard and MCTDH methods originally implemented in
the Quantics package. The choice of numbers of SPFs, mode combinations,
DVR-types, initial wavefunction, and integrator is the same as for the
standards grid-based methods. As such the input file should contain a
`run-section`, `pbasis-section`, `sbasis-section`,
`init_wf-section`
and optional `integrator-section`. Please note that only the VMF
integrator has been implemented so far for MCTDH calculations, so the
program will give an error if CMF is selected.

In order to indicate to the program that a DD-GB calculation is
requires, the keyword `dd-gb=nmodes` should be added to
the `run-section` (currently only normal mode coordinates have
been implemented).

In common with the DD-vMCG method, we need to provide certain files in order for the necessary electronic structure calculations to be performed so as to be able to generate the PES.

As with the DD-vMCG method it is necessary to give the program an
initial, reference geometry, in Cartesian coordinates, for the molecule
under consideration. If using normal modes as the coordinate system, it
is also necessary to list the modes along with their frequencies
(generated by a previous calculation using an electronic structure
program). Both the Cartesian coordinates and normal mode list are
included in the `DD-GB-section`, the latter between
keywords `nmode` and `end-nmode`. This section is thus
analogous to the `initial-geometry-section` used in DD-vMCG
calculations.

The `DD-GB-section` differs from
the `initial-geoemtry-section` in what is omitted: as the initial
wavefunction is defined in the `init_wf-section`, neither the
initial state (`init_state` keyword) nor the initial coordinates
of the wavefunction are defined (listed in the `nmode` subsection
when using vMCG) in this section. Similarly, because the grid-basis is
defined in
the `pbasis-section`, the `nstates` keyword is omitted as
are the keywords describing widths and positions of the GWPs and length
of analysis grid.

Please note that, when using normal modes, the names of the modes in
the `DD-GB-section` must match those used to describe the modes in
the `pbasis`, `sbasis` and `init_wf` sections.

The `Dirdyn-GB-section` is also needed for DD-GB calculations and
follows the form of the similarly named section when using
DD-vMCG. Crucially this section is used to define the directory location
of the database of electronic structures energies, the electronic
structure method to be used and the command to use to actually run the
electronic structure program..

In the following a brief description of how to set up and run direct dynamics calculations is given. For more details and the options available, please look in the main manual. An input file will be required, and a directory set up to contain the quantum chemistry files and DB. What is described below is not the only way to set up these calculations, and more control over the integration and other advanced parts of the program can be obtained if desired.

A simple example of a direct dynamics input is shown for the
internal conversion in water.
Six sections are needed, a `RUN-SECTION`, a `PBASIS-SECTION`,
an `SPF-BASIS-SECTION`, a `DDGB-SECTION`,
a `DIRDYN-GB-SECTION` and a `INIT_WF-SECTION`. A
direct dynamics calculation is specified by the keyword
`dd-gb=nmodes` in addition to the keyword `propagation` in
the RUN-SECTION. Other keywords such as `auto` and `tout`
are as for standard dynamics calculations.

The `PBASIS-SECTION`, `SPF-BASIS-SECTION`
and `INIT_WF-SECTION` have the same form as for a standard
grid-based calculation. Details can be found
in the Input Documentation

The `DIRDYN-GB-SECTION` contains keywords relevant for the
direct dynamics. The most basic is the name of the directory
containing the quantum chemistry files defined by `data =
h2o-4s-ddgb_dddata` . Also the `qcprogram` and `method`
keywords must be given to specify the quantum chemistry program to
be used and the electronic structure method. This allows the input
files to be correctly set up and output files read by the interface
between Quantics and the quantum chemistry program. The final
important keyword is `subcmd`. This gives the command that
would be used on the command line to submit the quantum chemistry
job. In the example the computer being
used submits to Molpro via a script `run_molpr`. The 2 at the
end means the command takes 2 arguments - the input file name and optional
directory.

Other keywords specify how the DB will be used and control the method of
Gaussian fit used. The zero point
for the energy is specified by `ener0`. The `ddmorder` keyword
specifies which kernel to use for the GPR expansion (0 indicates F-D
Gaussian, 1 a sum of 1-D Gaussians, 2 a sum of 1- and 2-D
Gaussians and 3 a sum of 1-, 2- and 3-D
Gaussians). In addition the `ddgbsvd` keyword turns on the
SVD based fitting algorithm, whilst the `ddgbtendec` keyword can be used to
switch on the SVd fitting as well as the further 3-D tensor fitting. For a full list
of options
see the main manual.

The easiest way to set up the initial geometry, as done in the
examples here, is using an
`DD-GB-SECTION` in which the molecule is
specified between the keywords `cartesian end-cartesian` in
XYZ format. In addition, the calculation should be run in normal modes,
with the normal mode information being specified between `nmode
end-nmode`. This gives the labels and frequencies of the modes.
The Cartesian coordinates information is then ignored except the
atom types. To set up the transformation matrices between cartesian
and normal mode coordinates a file can be read containing the
results of a quantum chemistry vibrational frequency calculation.
For example if
a Molpro
calculation with the keyword `frequencies` has been run,
this can be put in the data directory and the keyword `transfile =
water-opt-nosym.out` in the `DIRDYN-GB-SECTION` specifies it
should be used.

It is also necessary to place certain extra files in the database directory, particularly the template.dat file which gives a skeleton input file, which is used to create input files (by adding e.g. the current geometry) for the electronic structure program. More details on these files will be given below.

## Using Molpro

Using Molpro as the electronic structure method, it is possible to use the propagation diabatisation method for CASSCF for two or more electronic states.

A complete example for 4 state water is included in the `inputs`
directory, the input file being `h2o-4s-dd.inp` and the associated
directory for the database etc. being `h2o-4s-dd_dddata`. The
following explanation uses that example.

In order to perform DD with Molpro using normal modes, it is first necessary to optimise the molecular geometry using CAS with Molpro, and also to run a frequency calculation. It is best to do this using symmetry as otherwise accuracy can be lost, however, when using state-averaged CAS in earlier versions of Molpro the frequencies are not available when using symmetry, so an initial optimisation can be done with symmetry, before running a frequency calculation at the optimised geometry, without symmetry. The file produced contains the normal modes which are used in transforming between Cartesians and normal modes in Quantics. The next calculation to perform is another Molpro calculation using the same active space and the optimised geometry, but this time without symmetry, in order to generate a file containing the MO coefficients which will be used as an initial guess for the calculations performed during the DD (this can also be combined with the no-symmetry frequency calculation).

Both the frequency output and orbital files should be copied to the
auxiliary directory where the `template.dat` and `template_ao.dat`
files, and the database are kept.

An example template file is given in the Molpro template.dat file which provides the commands necessary to run the correct Molpro calculations as the dynamics proceeds.

We also need a template_ao file, an example of which is given in template_ao.dat. This is needed to generate the overlap matrix between basis functions at neighbouring geometries, which is then used to overlap the CAS wavefunctions at those geometries. Doing this allows us to track the signs of the wavefunctions of each state (which can flip randomly), hence correcting the signs of the non-adiabatic coupling terms, and thus allowing the propagation/integration of the NACTs to proceed smoothly.

Please note that individual calculations are needed to calculate the non-adiabatic couplings between each pair of states (2-1, 3-1, 3-2 etc.). Also note that we have to produce a punch file, containing the coefficients of the CI vectors, for use in the wavefunction overlaps described above.

Note that a file
called init-orbital.orb is
required to provide the initial guess orbitals for the CASSCF
calculation. This file must have been produced using a Molpro calculation
with symmetry turned off; the Molpro calculations performed during the
dynamics are done without symmetry (because most geometries sampled will
be of C1 symmetry anyway), therefore the guess orbitals must have been
produced in this way, otherwise Molpro produces a list of symmetrised
orbitals which will fail to be read in by the new calculations.
Subsequent calculations take their orbitals from files created by Quantics
(MO coefficients from each Molpro calculation are stored in the DB), and
this occurs in the $Orbread:....$ construct. The name of the orbital file
is filled in at `orbfile`.

Of note in the input file are the keywords in the `DIRDYN-GB-SECTION`:

`qcprogram=molpro`: Alerts Quantics to look for Molpro
inputs/outputs.

`method=cas`: Alerts Quantics to look for Molpro (SA)-CASSCF
inputs/outputs.

`dd_diab=global`: This command causes the program to use the
propagation diabatisation method to generate diabatic PESs from the
Molpro data.

`nbasis=X`: This gives the total number of basis functions
(contracted) used by the electronic structure calculations. Needed to set
up array sizes properly in the program.

`nelectrons=X`: This gives the total number of electrons. Needed to overlap
the wavefunctions at adjacent geometries.

`nocc=X`: This gives the total number of occupied orbitals i.e. the closed
orbitals plus those in the active space. Needed for the overlap of
the wavefunctions at adjacent geometries.

`nclosed=X`: This gives the total number of closed (doubly-occupied)
orbitals. Needed for the overlap of the wavefunctions at adjacent geometries.

`subcmd=...`: Gives the location (on your machine) of the script
to run Molpro.

`dbsave`: Causes the DB to be saved in RAM during the course of
the dynamics. This saves a lot of time. The DB is still written to file after each
period of sampling.

`gapvartol=X`: Sets the GPR variance parameter used to determine whether to
calculate energies at a sampled geometry. Larger values give a coarser sampling.

`gapsample=X`: Gives the number of geometries, around the centre of the
wavepacket, to sample on each state.

Alternatively to the propagation diabatisation method we can use the Procrustes method, which relies on overlapping CAS wavefunctions at neighbouring geometries, then rotating one of them to best match the other. The same setup is used as for the propagation diabatisation, the main differences being the use of a template.dat file without the need to calculate NACTs and the use of the keyword:

`dd_diab=procrustes`: This command causes the program to use the
Procrustes diabatisation method to generate diabatic PESs from the
Molpro data.

Example files for a Procrustes diabatisation with Molpro are given: but.inp, template.dat, template_ao.dat, init-orbital.orb, and but-opt.out.

## Studying states of different spin multiplicities: including SOCs

With DD-GB methods we can also study systems with states of different spin
multiplicities, e.g. singlets with triplets, when using Molpro. This can be done both
with the propagation and Procrustes diabatisation methods, where the states of each
individual multiplicity are diabatised (e.g. all the singlets are diabatised,
followed by the triplets), after which the spin-orbit couplings (SOCs) between the
two manifolds are transformed into the diabatic representation; we thus end up with
diabatised, pure spin states. Alternatively, with the propagation diabatisation
method only, we can diagonalise the spin-orbit Hamiltonian matrix, which produces a
set of adiabatic states, which are not pure spin states due to the mixing of the
SOCs. The matrix which diagonalises the SOC Hamiltonian is then used to transform the
NACTs into that basis, after which the propagation diabatisation procedure is applied
to all states at once. We thus finish with a set of diabatic, non-pure spin states.
The keywords `spinbas` (default) and `socbas` switch between these
methods.

Example files for a Procrustes diabatisation with SOCs (spinbas) using Molpro are given: tform.inp, template.dat, template_ao.dat, init-orbital.orb, and tform-sacas-freq.out.

Note the more complicated `template.dat` file, which uses MRCI (with no
excitations, so actually CASSCF) to calculate the SOCs.

In the `tform.inp` file we have the following extra keywords:

`multmap=a,b,c,d,e,...`: A list of the multiplicities of each state.

`statemap=a,b,c,d,e,...`: A list of the states showing which ones are
degenerate, so in the example the first 2 states are singlets, so they have no
degenerate complements, whilst the final 6 states are triplets, but ordered so that
the degenerate triplets are states (3,5,7) and (4,6,8). The individual states differ
in the SOCs with other states.

`socdegen`: This keyword, which does not work with the SVD fitting yet, can
be set so as to use a single state to represent each degenerate set of states. This
is achieved by taking the magnitude of the SOCs between each set of states. The
`multmap` and `statemap` keywords must be altered when doing this, e.g.
in the above example `multmap=1,1,3,3` and `statemap=1,2,3,4`.

## Setting up Operators for Analysis or to Augment the System Hamiltonian

As no `operator` file is used, if an operator is required
for analysis, e.g. to calculate an expectation value, this can be
added to the `input` file. Simply add the operator in a
section `HAMILTONIAN-SECTION_oper`, where *oper* is
the name of the operator. The expecation value of the operator can be
obtained during the course of the propagation by
adding `expect=oper` to the `RUN-SECTION`

By adding a section called `HAMILTONIAN-SECTION_system`
(i.e. `oper` in the above is set to `system`, extra
terms can be added to the System Hamiltonian to modify the wavepacket
propagation. This can be used, for example, to
append a CAP operator to the system Hamiltonian (see also
the `LABELS-SECTION`) allowing the wavepacket to be annihilated
as dissociation occurs.

Both a system operator and an extra operator use all the MCTDH operator
input and a `PARAMETER-SECTION` and `LABELS-SECTION` can
also be used. Note that multiple `analysis` operators can be
added (add `expect=oper1,oper2,...` to the `RUN-SECTION` to
output some or all expectation values) as
well as a `system` operator.