Doing a dynamical walk

The theory behind the ``direct dynamics'' as implemented in DALTON is described in Ref. [28]. The main idea behind this approach is that Newton's equations of motion for the nuclei are integrated in the presence of the quantum mechanical potential set up by the electrons. Thus one may follow a molecular reaction from a given starting point (usually a transition state) as it would behave if the nuclei could be treated exactly as classical particles. One should also keep in mind that the Hamiltonian used is constructed within the framework of the Born-Oppenheimer approximation, which may turn out not to be a good approximation at given points during the reaction. Furthermore, the calculation describes the way molecules with a predefined orientation and momentum will react. Thus the trajectory obtained is only one of a large number of possible trajectories depending on the initial state of the molecule.

The necessary input in order to do a dynamical walk of for instance protonated formaldehyde would look like:

**DALTON INPUT
.WALK
.MAX IT
 200
*WALK
.DYNAMI
.FRAGME
 5
 1 1 1 2 2
.MOMENT
 1
1 -.00001
.MODE
 1
**WAVE FUNCTIONS
.HF
**END OF DALTON INPUT

The walk is specified to be a dynamic walk through the keyword .DYNAMI. The starting trust radius will in dynamical calculations be changed to a new default value of 0.005.

The keyword .FRAGME dictates which atoms belong to which molecular fragment. In this particular case, we assume that two protons leave the protonated formaldehyde as a hydrogen molecule, and that the leaving hydrogen atoms are the last atoms of the MOLECULE input file. This partitioning is mainly needed in order to get proper values for the relative translational energy between the two fragments, as well as for deciding how much of the energy has been distributed into internal degrees of freedom.

The default isotopic substitution is that the most abundant isotopes are to be used in the calculation. Isotopic substitution is important as the masses of the nuclei enters when Newton's equations of motion are integrated. The specification of the isotopic constitution of the molecule is given in the MOLECULE.INP file, as described in Chapter 23.

In this calculation we start the calculation at a transition state, and in order to get the reaction started we need to give the molecule a slight push. This is achieved by the keyword .MOMENT. In the next line the user then specifies the number of modes in which there is an initial momentum, followed by lines containing pairs of numbers, of which the first specifies the mode, and the second the momentum in this mode. There must be as many pairs of modes and momenta as specified in the line after the .MOMENT keyword. It is impossible to predict in advance which way the reaction will proceed, and the calculation should be checked after a few iterations, in order to ensure that it proceeds in the right direction. If not, the calculation should be started from the transition state again with a different sign on the initial momentum. The DALTON.TRJ must also be removed as discussed below.

It is in principle possible to start a calculation from any point on a molecular potential energy surface, and in cases where these starting points do not correspond to a stationary point, .MOMENT may be skipped, as there exist a downward slope (in other words, an attractive force) driving the molecule(s) in a specific direction. One may of course also start the molecule with a given initial momentum in different energy modes.

During the dynamical calculation, care has also to be taken in order to ensure that the steps taken are not too long. If this occurs, the initial trust radius and/or the trust radius increment should be reduced by the keyword .TRUST. In the DALTON output one will find ``Accumulated kinetic energy since start'', and this property will be calculated in two ways: From conservation of the total energy, and from integrated momenta. If the difference between these numbers is larger than approximately 1% , the calculation should be stopped and the starting trust radius be decreased and the calculation restarted from the starting point again after removal of the DALTON.TRJ-file.

The calculation of dynamical walks may take from about 70 to 1200 iterations (as a general rule) and one must therefore adjust the maximum number of iterations allowed. This is done by the .MAX IT keyword. In the present example the maximum number of iterations have been reset to 200. If the calculation cannot be closely monitored, it is recommended not to set the maximum number of iterations too high, and rather restart the calculation if this turns out to be necessary. This can be accomplished by specifying the iteration at which the calculation will restarted by the keyword .ITERAT in the **DALTON input module.

The calculation should be stopped (at least for ordinary hydrogen elimination reactions) when the ``Relative velocity'' starts to decrease, as this indicates that the molecules are so far apart that basis set superposition errors become apparent. We will below return to how one then calculates translational energy release of the reaction.

During the whole calculation, a file DALTON.TRJ is updated. This file contains information from the entire dynamic walk. If a walk is restarted from a given point, the new information will be appended to the old DALTON.TRJ-file. Note that this also implies that if you need to restart the calculation from the beginning (because the reaction went the wrong way or because of a too large trust radius), the DALTON.TRJ-file must be removed. Thus, it may often be advisable to take a backup of this file in certain parts of the calculation. As this file contains all the information about the dynamical walk, this file can be used to generate a video-sequence of the molecular reaction along this specific trajectory with the correct time-scaling [38].