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.