MCTDH: Direct dynamics

Back to Howto Menu

Back to Main Menu

Running Direct Dynamics Calculations

In the following a brief description of how to set up and run direct dynamics calculations is given. For more details and options available 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 vibration of the O-H bond in ethanol . Three sections are needed, a RUN-SECTION, an INITIAL-GEOMETRY-SECTION and a DIRDYN-SECTION A direct dynamics calculation is specified by the keyword direct in addition to the keyword propagation in the RUN-SECTION. Options can be specified such as direct = cartesian or direct = nmodes to specify the coordinates in which the dynamics will be run. The keyword ngwp = N will specify a basis set of N Gaussian wavepackets. Other keywords such as auto and tout specify the what is calculated and how often.

The DIRDYN-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 = etoh_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 MCTDH 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 GAUSSIAN via a script run_g03. The 2 at the end means the command takes 2 arguments - the input file name and directory where the results will go.

Other keywords specify how the DB will be used. The zero point for the energy is specified by ener0. It can be selected whether you use an ascii_db or a binary_db. The latter is smaller, while the former is transferable between platforms and is human readable. It is possible to convert the format after a calculation using the convdb program. (The calculation always uses a binary database, the ascii_db keyword runs convdb at the end of the calculation to convert the database). Full 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 INITIAL-GEOMETRY-SECTION in which the molecule is specified between the keywords cartesian end-cartesian in XYZ format. If the dynamics are run in Cartesians, keywords are then used to specify, e.g. the width of the GWPs.

If the calculation is to be run in normal modes, the normal mode information needs to be specified between nmode end-nmode. This gives the label, frequency of the modes and any initial displacements. 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 GAUSSIAN calculation with the keyword freq=hpmodes has been run, this can be put in the data directory and the keyword transfile = etohdistfr.log in the DIRDYN-SECTION specifies it should be used.

It is also necessary to set up the information required by the quantum chemistry calculation and deposit it in a directory. The following two files are required:

For multi-state calculations the states are diabatised using either a simple regularisation or a global scheme. For the regularized diabatic states, the keyword dd_diab = regularized is specified in the DIRDYN-SECTION and information about the conical intersection point used for the regularisation must be given. This can be read from a file using the coinfile and dreffile keywords. An example based on the butatriene cation using CASSCF is given in Example 3 with the data in Example 3 data files . The template.dat file includes metastrings to generate the 2 surfaces. Note that the number of basis functions also has to be specified in the DIRDYN-SECTION using the nbasis = I keyword.

To use global diabatic states, the keyword dd_diab = global is specified and no information about the conical intersection point is required. An example is given in Example 4 with the data in Example 4 data files .

Analysis programs: gwptraj, showsys

Using Gaussian

The above examples use Gaussian09 to run either HF or CASSCF calculations. There a number of metastrings possible. They are case sensitive. Put the metastrings starting with $ in the keywords and they will be sustituted in the final input file. The exceptions are $chk$ and $cpus$ that go in the usual place for specifying the checkpoint file name and nproc. The metastrings will de substituted during the creation of input files as described. Any other lines will be simply copied. Remember that Gaussian input files are sensitive to blank lines. supplied with the initial MOs. This file is copied automatically to the checkpoint file to be used.
Gaussian Metastrings.
Metastring Description
The following metastrings closed by < > will be substituted by blocks of data and must be put in correct places
<geometry> present geometry will be written here
<mo-coefficients> MO coefficients used as a guess in a CAS calculations written here when guess supplied by $Guess:Guess=Cards$
The following metastrings are replaced by the relevant values.
$chk$ will be replaced by a checkpoint file name to match the generated input file name. Use with %chk=$chk$
$cpus$ will be replaced by the no. of cores available in a parallel job. Use with %nproc=$cpus$.
$root$ The number of roots to be calculated. For example in a CAS calculation use the keyword CAS(M,N,nroot=$root$)
The following metastrings provide keywords depending on the calculation being run.
$Swap:keyword$ Need to swap the order of states to calculate the gradient of the lower state in a CAS calculation. Typically $Swap:IOp(5/97=100,10/97=100)$ is used so that the IOPs are selected when required.
$Freq:keyword$ The keyword to be used to calculate the frequency. For a simple calculation use $Force:Freq$. For a multi-state CAS calculation $Freq:Freq IOp(5/17=41000200,10/10=700007)$ means that the derivative couplings are also calculated
$Force:keyword$ If the frequency is not being calculated (not required after the first step if using Hessian updating) then the keywords specified by the $Force:$ metastring are used. For example in a multi-state CAS calculation use $Force:Freq IOp(5/17=41000200,10/10=700007,7/10=0)$ to calculate the DCP without diagonalising the Hessian.
$Guess0:keyword$ keyword specifies how the initial guess in a CAS calculation will be made if the DB is empty. Typically "Guess0:Guess=Read" is used and a file start.chk
$Guess:keyword$ keyword specifies how the initial guess in a CAS calculation will be made if the DB is not empty. Typically "Guess:Guess=Cards" is used and the MO coefficients are taken from the nearest DB point and written in cards format in place of the metastring <:mo-coefficients>

Using Molpro

Using Molpro as the electronic structure method, it is possible to use the propagation diabatisation method for CASSCF for more than two 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. It is best to do this using symmetry as otherwise accuracy can be lost. 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.

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

An example template file is given in Molpro template.dat.

In the template file, note that the $Orbread0:....$ construct, which reads in guess orbitals for the first calculation, has the MO coefficient file called init-orbital.orb. 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 .

Also note that separate CASSCF calculations are needed to calculate the gradients/Hessians of each state and individual calculations are needed to calculate the non-adiabatic couplings between each pair of states (2-1, 3-1, 3-2 etc.)

The Quantics input file for Molpro is much the same as for Gaussian. It can be viewed at h2o-4s-dd.inp.

Of note in the input file are the keywords:

use_dd_sym: this turns on the use of symmetry when using DD in normal modes. It relies on the term symbols for each mode to tell whether the PES along that direction should be symmetric. In order for this to work properly a full shell of GWPs should be used so that GWPs are evenly distributed along the modes (the log file tells you whether the shell is full). The minimal number of GWPs for an N mode system is 2N+1 i.e. 1 in the centre, 2 along each mode either side of the centre. N.B> This only works when the centre of the initial wavepacket is at 0 along all symmetric potentials. In this case the wavepacket is initially centred at 0.0 along the 3B2 mode, which has a plane of symmetry at 0.0 in the PES.

grid X Y Z: In the initial-geometry-section, this set the length of the grid (from Y to Z) and number of points (X) along each mode for analysis purposes. The default is X=101, Y=-10, Z=10.

default: This forces the program to use the integrator accuracy defined rather than performing the automatic loosening of a factor of 10.

gwp_eoms=precond: This causes the program to use a preconditioning method to generate the EOMs for the GWPs (to solve iCL=Y).

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

dd_diab=global: This command causes the program to use the propagation diabatisation method to generate diabatic PESs from the Molpro data. This must use the hess_upd keyword (see below).

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.

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

hess_upd: The program uses the Powell update method to create approximate Hessians of the PESs at each DB point after an initial Hessian has been created at the first DB geometry. The propagation diabatisation only works with this. It can also be used with other methods, and can save a huge amount of time.

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 at each new point.

rediabatise: This is commented out here as it is a post-processing option. Along with the genoper keyword in the run-section (don't use propagation), this takes a pre-created DB of adiabatic points and performs the propagation diabatisation again, starting a the central point and moving outwards, diabatising in order of distance from the centre (updating the gradients and Hessians as well). This can smooth surfaces out particularly if a GWP has left the central region then come back again, creating DB points as it goes. The diabatisation is approximate, so shorter paths are to be preferred between points, this allows that.

HAMILTONIAN-SECTION_system: See section below as well. This appends a CAP operator to the system Hamiltonian (see also the LABELS-SECTION) allowing the wavepacket to be annihilated as dissociation occurs.

Using Molcas

It is also possible to use Molcas as the quantum chemistry program. As with Molpro, the propagation diabatisation method for CASSCF allows the study of dynamics with more than two electronic states involved.

The Quantics input file for Molcas is the same as for Molpro or Gaussian. An example can be viewed at tol.inp. In this example, the Molcas electronic structure calculations are submitted via the script molcasrun.

A template example file is given in Molcas template.dat, where the properties of all electronic states are computed in a single Molcas calculation. For the first calculation, the guess orbitals are read from the file called start.RasOrb specified in the $Orbread0:....$ command. Subsequent calculations take their guess orbitals from files created by Quantics (MO coefficients from each Molcas calculation are stored in the DB), and this occurs in the $Orbread:....$ command.

Note that in the Molcas calculation, the 'RASSCF' module is called for each electronic state with the corresponding 'rlxroot' value. The following 'ALASKA' and 'MCKINLEY' modules compute then the gradient and Hessian respectively for each electronic state. The non-adiabatic couplings are finally computed between each pair of states.

Setting up Operators for Analysis

As no operator file is used, if an operator is required for analysis, e.g. to calculate and 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. If you want to append an operator e.g. a CAP to the system Hamiltonian then oper is set to be system. Both a system operator and an extra operator use all the MCTDH operator inout and a PARAMETER-SECTION and LABELS-SECTION can also be used.

Adding A Light Field

In direct dynamics the system Hamiltonian is generated automatically. To add extra terms, such as a light field, add the terms in the input file in a section HAMILTONIAN-SECTION_system. This uses all the MCTDH operator inout and a PARAMETER-SECTION and LABELS-SECTION can also be used. The operators qchemdipx, qchemdipy, qchemdipz and qchemdip are available to use the components or absolute value of the dipole moment as calculated from the electronic structure calculations. The keyword dipole must be added in the DIRDYN-SECTION to generate the database of dipole moments. An example input for a local control calculation is here