Resgrp:comp-photo-dyn
Benjamin / Charlotte / David / Marta to create this please! Some Documentation and some tutorials... MY COMMENTS IN UPPER CASE MAR
The direct quantum dynamics package
The direct quantum dynamics method we use is based on the propagation of DD-vMCG (direct dynamics variational multi-configuration Gaussian) wavepackets. This has been implemented in a development version of the Heidelberg MCTDH package [see web page] by Graham A. Worth [see web page] and collaborators.
This package consists of a set of programs for multidimensional quantum dynamics using the MCTDH algorithm (multi-configuration time- dependent Hartree). The vMCG method is a particular case of MCTDH, where the nuclear Global wavepacket (GWP) is expanded in a time-dependent basis set of Gaussian functions (gwp). The direct dynamics implementation allows the electronic energy (potential energy for the nuclei) to be calculated on-the-fly by the quantum chemistry package GAUSSIAN around the centre of each Gaussian function in the expansion, for both electronic states involved in the photochemical process under study.
This approach provides a time-resolved simulation of the photochemical mechanism induced by the initial condition chosen by the user. Non-adiabatic effects and tunneling are better described than with semi-classical trajectory methods such as surface hopping [see recent review].
General overview
The study of photochemical mechanisms using the DD-vMCG method can be approached in four stages.
1) The preliminary stage is a static description of the mechanism under study. This corresponds to a description of the main topological features of the potential energy surfaces of the coupled electronic states, as in any other study of the potential energy surfaces using CASSCF [Benzene CASSCF tutorial (G03)]. In practice, this involves:
- Choice of the active space and basis set for the electronic wavefunction
- Characterisation of the relevant critical points (stationary points and conical intersections)
- Characterisations of the main pathways (reaction paths and crossing seams)
- Definition – and possible selection – of the nuclear coordinates
Some comments:
- Choosing the active space may happen to be delicate when a full-valence description is too expensive (or just too large) to be used for frequency calculations, in particular if sigma and pi orbitals can mix. The user should keep in mind that the active space will be propagated along the trajectories followed by the centres of the Gaussian functions and may break in very distorted geometries.
- The main critical points involved in the mechanism must be identified, optimised and characterised. Stationary points are characterised by frequency calculations, and conical intersection points by examining the nature of the electronic states involved and the branching-space vectors that lifts the degeneracy. The important pathways can be characterised in the form of Mimimum Energy Path (MEP) or the Intrinsic Reaction Coordinate (IRC) linking stationary points. In addition, a development version of GAUSSIAN implemented in our group [Seam calculations of CHD] can be used to characterise the curvatures of the intersection space (crossing seam). The same program can calculate minimum energy paths within the intersection space, linking conical intersection points characterised as seam minima or seam saddle points.
- For now, the nuclear coordinates are defined as the normal coordinates of the electronic ground-state minimum (magnitude of nuclear displacements along the directions of the normal modes) for convenience. Their number can be reduced to keep the most relevant ones and thus simplify the study of the system dynamics. A method based on the development version of GAUSSIAN mentioned above may be used in this context [see reference].
2) The second stage is the generation of the input files. This can be performed automatically by a utility program that asks the user questions in order to set the correct options, but two calculations have to be done first:
- Frequency calculation at the ground-state minimum (log and chk or fchk files), used for (i) the initial geometry, (ii) the definition of the normal coordinates, and (iii) the initial active space.
- Optimised conical intersection (log file), used for the generation of diabatic electronic states from the adiabatic ones (when the so-called "regularization" scheme is used).
3) The third stage corresponds to the dynamics calculation itself. This part is managed by the MCTDH (DD-vMCG) program and should finish normally if the input file were written correctly. However, this may fail for various reasons, often because converging the active space becomes difficult. Also, calculations may have to be redone several times, with conditions that can be improved, for example by increasing the number of Gaussian functions (gwp) in the wavepacket expansion, or by changing the conical intersection point chosen for "regularization" scheme if this approach is to be followed.
4) Finally, the last step is the analysis of results. This involves looking at the time evolution of:
- Normal coordinates at the centre of each Gaussian function ("trajectories")
- Averages and widths of the normal coordinates for the global wavepacket (GWP)
- Adiabatic and diabatic energies at the centre of each Gaussian function ("trajectories")
- Electronic populations for the global wavepacket (GWP)
Some comments:
- Normal coordinates can be transformed into Cartesian or internal coordinates. Cartesian coordinates can be used to further represent trajectories as movies using a utility program.
- The quality of the calculation can be checked with various indicators. The simplest is the preservation of the mean total energy. If this increases too fast, the time step might have to be reduced.
- Stages 3 and 4, probably might have to be repeated several times before a satisfactory result is reached. Certainly, the spirit of this part is guided by the trial and error method. Playing with variables such as the number of nuclear Gaussian functions, the final propagation time, the propagation time step, or the integrator may be required.
Current available MCTDH (DD-vMCG) versions
mctdh90dev
This is the "old" and stable version of the code, which has been used by Marta Araujo and David Asturiol. This version relies on the multi-set formulation for DD-vMCG calculations (distinct Gaussian functions on different electronic states), which complicates the interpretation of results and often generates numerical problems. In addition, particular cases require particular versions of the source code to be recompiled.
N.B.: the potential energy surfaces are calculated in the adiabatic picture by GAUSSIAN and must be transformed to the diabatic picture for MCTDH. A simple scheme, based on the "regularized" diabatic states of Köppel et al. [see reference], has been used to date in which the gradient difference and non-adiabatic coupling at a point on the conical intersection seam can be used to estimate the transformation matrix. As a consequence, the user needs to provide the geometry and the branching-space vectors of a chosen conical intersection. This complicates the generation of the input files, as these three vectors are given in a Cartesian frame that must be re-orientated following some constraints discussed further in this manual. This also introduces a possible bias and also makes the interpretation of results quite delicate in terms of the assignment of the electronic populations in both representations. A work is in progress to implement a new scheme beyond these limitations.
Problems:
- Interpretation of results because of:
- Multi-set formulation: every gaussian wavepacket (gwp) has a different set of gaussian functions on each state.
- Lack of diabatic electronic populations for the global wavepacket (GWP)
- Lack of gaussian individual electronic populations
- Numerical problems: The multi-set formulation allows the individual gaussian wavepacket basis set to explore the two electronic states involved independently. Even thought the global wavepacket (GWP) is stable, the individual ones that form the basis set can explore silly regions of the PES and break the Active Space.
mctdh90.31dv
This version is more portable and more flexible in terms of options. The main new feature is the possibility of using the single-set formulation for DD-vMCG calculations (a unique set of Gaussian functions with varying occupancies on each electronic states).
New features:
- Single-set formulation: the gaussian functions on each state are identical. Therefore, the gaussian wavepacket is the same in both states, but it is represented by two weights corresponding to each electronic state.
Problems:
- Parallel calculations
- Lack of diabatic electronic populations for the global wavepacket (GWP)
- Lack of gaussian individual electronic populations
mctdh90.56dv and mctdh90.56dv_patch
This version should be used from now on, but has not been tested extensively yet.