This keyword requests a classical trajectory calculation [Bunker71, Raff85, Hase91, Thompson98] using a Born-Oppenheimer molecular dynamics model (first described in [Helgaker90, Uggerud92]; see [Bolton98] for an extended review article). The implementation in Gaussian 16 [Chen94, Millam99, Li00] extends the usual methodology by using a very accurate Hessian-based algorithm that incorporates a predictor step on the local quadratic surface followed by a corrector step. The latter uses a fifth-order polynomial or a rational function fitted to the energy, gradient, and Hessian at the beginning and end points of each step. This method for generating the correction step enables an increase in the step size of a factor of 10 or more over previous implementations.
The selection of the initial conditions using quasi-classical fixed normal mode sampling and the final product analysis are carried out in the same manner as in the classical trajectory program VENUS [Hase96]. Alternatively, initial Cartesian coordinates and velocities may be read in.
Note that the ADMP method provides equivalent functionality at substantially lower computational cost at the Hartree-Fock and DFT levels.
All BOMD jobs must specify the number of dissociation paths; for many jobs, this value will be zero (a blank line is also allowed), and no other BOMD input will be used. In this case, the trajectory is integrated for a fixed number of steps, as specified by the program default of 100 or the value of the MaxPoints option.
If NPath is set to –1, the dissociation pathways will be detected automatically and a gradient criteria (Hartree/Bohr) will be used in place of the usual fragmentation pathway and stopping criteria.
When the number of dissociation paths is greater than zero, the full BOMD job input has the following general structure:
| NPath |
Number of dissociation paths (maximum=20) |
| IFrag1, …, IFragNAtoms |
Fragmentation information |
| … |
…Repeated NPath times |
| [R1, R2, R3, R4, G5, ITest, IAtom, JAtom, R6 |
Optional stopping criteria (ReadStop option) |
| …] |
…Repeated NPath times |
| [Estart,DelE,SBeta,Ef,DPert,IFlag] |
Optional simulated annealing params. (SimAnneal) |
| [Mode-num, VibEng(Mode-num), …] |
Optional initial normal mode energies (NSample) |
| [Initial velocity for atom 1: x y z |
Optional initial velocities (ReadVelocity, ReadMWVelocity) |
| Initial velocity for atom 2: x y z
|
|
| …
|
|
| Initial velocity for atom N: x y z
|
|
| …] |
Entire section repeated NTraj times |
| [Atom1, Atom2, E0, Len, De, Be |
Optional Morse parameters for each diatomic product |
| …]
|
|
| |
Terminate subsection with a blank line |
The input line(s) following NPath define the fragmentation information for each path. The value in each position specifies that the corresponding atom belongs to the specified fragment number (i.e., atom i belongs to fragment number IFragi). Note that fragment information for each path must begin on a new line, but the ones for any individual path may be continued over as many lines as necessary.
Stopping criteria are specified next when the ReadStop option is included. Up to six stopping criteria may be specified for each path. A trajectory is terminated when all of the active criteria are satisfied. However, a value of zero for any parameter turns off testing for the corresponding stopping criteria. The stopping criteria tests are defined as follows (default parameter values are in parentheses):
-
Minimum distance between the centers of mass for any pair of fragments > R1 (18)
-
Minimum distance between atoms located in different fragments > R2 (20)
-
Maximum distance between the center of mass and any atom in the same fragment < R3 (0)
-
The maximum distance between any pair of atoms in the same fragment < R4 (0)
-
Interfragment gradient < G5 (10-6)
-
If ITest=1, distance between atoms IAtom and JAtom > R6 (0). Otherwise, distance between atoms IAtom and JAtom < R6 (0)
All distances are specified in Bohr, and the units of the gradient G5 are Hartrees/Bohr.
Parameters for simulated annealing/fragmentation follow the stopping criteria in the input stream when the SimAnneal option is specified:
-
Estart is the desired initial kinetic energy (Hartrees).
-
DelE is the energy gain/loss in Hartrees.
-
SBeta is the Fermi-Dirac inverse temperature (1/Hartrees).
-
Ef is the Fermi energy (wavenumbers): all modes corresponding to a frequency in wavenumbers below Ef will be enhanced, while those above Ef will be reduced. The reverse will happen if SBeta is negative.
-
DPert is the size of the random perturbation.
-
IFlag determines the algorithm for applying an energy perturbation for simulated annealing (i.e., adding/removing energy from the eigenmodes). It has the following possible values: 0 (weigh each eigencomponent according to its frequency), 1 (add DelE in a random fashion), 2 (combination of 0 and 1), 00 (if near a transition state, add all energy along that mode), and 10 (ignore any nearby transition state).
The next part of the input specifies how much energy is in each normal mode when the NSample option is used. For each mode, VibEng is the translational energy in kcal/mol in the forward direction along the transition vector. If VibEng < 0, then the initial velocity is in the reverse direction. (You can explicitly specify the forward direction using the Phase option.)
Next, the initial velocity for each atom is read if the ReadVelocity or ReadMWVelocity option is included. Each initial velocity is specified as a Cartesian velocity in atomic units (Bohr/sec) or as a mass-weighted Cartesian velocity (in amu1/2*Bohr/sec), respectively. One complete set of velocities is read for each requested trajectory computation.
Finally, Morse parameter data can be specified for each diatomic product. The Morse parameter data is used to determine the vibrational excitation of diatomic fragments using the EBK quantization rules. It consists of the atomic symbols for the two atoms, the bond length between them (Len, in Angstroms), the energy at that distance (E0 in Hartrees), and the Morse curve parameters De (Hartrees) and Be (Angstroms-1). This input subsection is terminated by a blank line.
MaxPoints=n
Specifies the maximum number of steps that may be taken in each trajectory (the default is 100). If a trajectory job is restarted, the maximum number of steps will default to the number specified in the original calculation.
Phase=(N1, N2[,N3[,N4]])
Defines the phase for the transition vector such that forward motion along the transition vector corresponds to an increase in the specified internal coordinate. If two atom numbers are given, the coordinate is a bond stretch between the two atoms; three atom numbers specify an angle bend and four atoms define a dihedral angle.
ReadVelocity
Read initial Cartesian velocities from the input stream. Note that the velocities must have the same symmetry orientation as the molecule. This option suppresses the fifth-order anharmonicity correction.
ReadMWVelocity
Read initial mass-weighted Cartesian velocities from the input stream. Note that the velocities must have the same symmetry orientation as the molecule. This option suppresses the fifth-order anharmonicity correction.
SimAnneal
Use simulated annealing (the initial velocity is randomly generated). Additional parameters are read in, as described above.
Only one of ReadVelocity, ReadMWVelocity and SimAnneal can be specified.
ReadStop
Read in alternative stopping criteria.
RTemp=N
Specifies the rotational temperature. The default is to choose the initial rotational energy from a thermal distribution assuming a symmetric top (the temperature defaults to 0 K).
NSample=N
Read in initial kinetic energies for the first N normal modes (the default is 0). The energies for the remaining modes are determined by thermal sampling by default.
NTraj=N
Compute N trajectories.
Update=n
By default BOMD does second derivatives at every point. Using the Update keyword causes the program to perform Hessian updates for n gradient points before doing a new analytic Hessian. GradOnly requests that only gradients be done and that the Hessian be updated all the time (full second derivatives are not computed). ReCalcFC is a synonym for this option.
RandomVelocity
Randomly generate initial velocities for dynamics without using any second derivative information. This is the default for GradientOnly dynamics.
StepSize=n
Sets the dynamic step size to n*0.0001, in the appropriate units. The default step size is 0.25 amu1/2*Bohr except for GradientOnly calculations where it is 0.025 femtoseconds.
Sample=type
Specifies the type of sampling, where type is one of these keywords: Microcanonical, Fixed, and Local. The default is Fixed normal mode energy unless RTemp was specified, in which case Local mode sampling is implied.
Restart
Restart a trajectory calculation from the checkpoint file. Note that options set in the original job will continue to be in effect and cannot be modified.
ReadIsotopes
This option allows you to specify alternatives to the default temperature, pressure, frequency scale factor and/or isotopes—298.15 K, 1 atmosphere, no scaling, and the most abundant isotopes (respectively). It is useful when you want to rerun an analysis using different parameters from the data in a checkpoint file.
Be aware, however, that all of these can be specified in the route section (Temperature, Pressure and Scale keywords) and molecule specification (the Iso parameter), as in this example:
#T Method/6-31G(d) JobType Temperature=300.0 …
…
0 1
C(Iso=13)
…
ReadIsotopes input has the following format:
| temp pressure [scale] |
Values must be real numbers. |
| isotope mass for atom 1
|
|
| isotope mass for atom 2
|
|
| …
|
|
| isotope mass for atom n
|
|
Where temp, pressure, and scale are the desired temperature, pressure, and an optional scale factor for frequency data when used for thermochemical analysis (the default is unscaled). The remaining lines hold the isotope masses for the various atoms in the molecule, arranged in the same order as they appeared in the molecule specification section. If integers are used to specify the atomic masses, the program will automatically use the corresponding actual exact isotopic mass (e.g., 18 specifies 18O, and Gaussian uses the value 17.99916).
All semi-empirical, HF, CASSCF, CIS, MP2 and DFT methods.
The following sample BOMD input file illustrates many of the available options. It will calculate a trajectory for H2CO dissociating to H2 + CO, starting at the transition state. There is one fragmentation pathway: C and O belong to fragment 1, and the two hydrogens belong to fragment 2.
Stopping criteria are also specified in this example job. The trajectory will be stopped if the distance between the centers of mass of H2 and CO exceed 13 Bohr, the closest distance between H2 and CO exceeds 11 Bohr, all atoms in a fragment are less than 1.3 Bohr from the center of mass of the fragment, any atom in the fragment is less than 2.5 Bohr from all other atoms in the fragment, the gradient for the separation of the fragments is less than 0.0000005 Hartree/Bohr, and the distance between atoms 1 and 3 is greater than 12.8 Bohr.
The initial kinetic energy along the transition vector is 5.145 kcal/mol, in the direction of the products (the forward direction is characterized by an increase in the larger C-H distance). The Morse parameters for H2 and CO are specified to determine the vibrational excitation of the product diatomics; they were computed in a previous calculation. The calculation will be carried out at 300 K.
| # HF/3-21G BOMD(Phase=(1,3),RTemp=300,NSample=1,ReadStop) Geom=Crowd |
|
| |
|
| HF/3-21G dissociation of H2CO –> H2 + CO |
|
| |
|
| 0 1 |
|
| C |
|
| O 1 r1 |
|
| H 1 r2 2 a |
|
| H 1 r3 3 b 2 180. |
|
| |
|
| r1 1.15275608 |
|
| r2 1.74415774 |
|
| r3 1.09413376 |
|
| a 114.81897892 |
|
| b 49.08562961 |
|
| |
|
| 1 |
|
| 1 1 2 2 |
|
| 13.0 11.0 1.3 2.5 0.0000005 1 1 3 12.8 |
|
| 1 5.145 |
|
| C O -112.09329898 1.12895435 0.49458169 2.24078955 |
|
| H H -1.12295984 0.73482237 0.19500473 1.94603924 |
|
| |
Final blank line |
Note that all six stopping criteria are used here only for illustrative purposes. In most cases, one or two stopping criteria are sufficient.
At the beginning of a BOMD calculation, the parameters used for the job are displayed in the output:
TRJ-TRJ-TRJ-TRJ-TRJ-TRJ-TRJ-TRJ-TRJ-TRJ-TRJ-TRJ-TRJ-TRJ-TRJ-TRJ-TRJ
-------------------------------------------------------------------
INPUT DATA FOR L118
-------------------------------------------------------------------
General parameters:
Max. points for each Traj. = 100
Total Number of Trajectories = 1
Random Number Generator Seed = 398465
Trajectory Step Size = 0.250 sqrt(amu)*bohr
Sampling parameters:
Vib Energy Sampling Option = Thermal sampling
Vib Sampling Temperature = 300.0 K
Sampling direction = Forward
Rot Energy Sampling Option = Thermal distribution (symmetric top)
Rot Sampling Temperature = 300.0 K
Start point scaling criteria = 1.000D-05 Hartree
…
Reaction Path 1
****************
Fragment 1 center 1 ( C ) 2 ( O )
Fragment 2 center 3 ( H ) 4 ( H )
Termination criteria:
The CM distances are larger than 13.000 bohr
The min atomic distances among fragments are larger than 11.0 bohr
The max atomic and CM distances in frags are shorter than 1.3 bohr
The max atomic distances in fragments are shorter than 2.500 bohr
The change of gradient along CM is less than 5.00D-07 Hartree/bohr
Distance between atom center 1 ( C ) and 3 ( H ) is GE 12.800 bohr
Morse parameters for diatomic fragments:
E0 Re De Be
C O -112.0932990 1.1289544 0.4945817 2.2407896
H H -1.1229598 0.7348224 0.1950047 1.9460392
---------------------------------------------------------------------
The initial kinetic energies for the normal modes appear at the beginning of each trajectory step:
-------------------------------------------------------
Thermal Sampling of Vibrational Modes
Mode Wavenumber Vib. quant.# Energy (kcal/mol)
-------------------------------------------------------
1 -2212.761 5.14500
2 837.330 0 1.19702
3 1113.182 0 1.59137
4 1392.476 0 1.99064
5 2026.859 0 2.89754
6 3168.689 0 4.52987
-------------------------------------------------------
After the trajectory computation is complete, summary information is displayed in the output:
Trajectory summary for trajectory 1
Energy/gradient evaluations 76
Hessian evaluations 76
Trajectory summary
Time (fs) Kinetic (au) Potent (au) Delta E (au) Delta A (h-bar)
0.000000 0.0214192 -113.0388912 0.0000000 0.0000000000000000
1.169296 0.0293490 -113.0468302 -0.0000091 0.0000000000053006
2.161873 0.0407383 -113.0582248 -0.0000144 0.0000000000045404
…
The information is given for each time step in the trajectory. In addition, the output includes geometrical parameters for the atoms in each fragment, the distances between fragments, and the relative mass-weighted velocities for each fragments and between fragments, all reported at each time step. You can also use GaussView or other visualization software to display the trajectory path as an animation.
