Gaussian 03 Online Manual
Last update: 2 October 2006

BOMD

This keyword requests a classical trajectory calculation [177,178,179,180] using a Born-Oppenheimer molecular dynamics model (first described in [181,182]; see [403] for an extended review article). The implementation in Gaussian 03 [184,186,187] 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 [404]. 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.

REQUIRED INPUT

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 
Initial velocity for atom 2: x y z                      (ReadVelocity or ReadMWVelocity) 
... 
Initial velocity for atom N: x y z 
...]                                                  Entire section is repeated NTraj times 
[Atom1, Atom2, E0, Len, De, Be                         Optional Morse params. 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 specified. 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, whole 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-weighed 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 Hessians update 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).

StepSize=n
Sets the step size in dynamics to n*0.0001 femtoseconds.

Sample=type
Specifies the type of sampling, where type is one of these keywords: Orthant, 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.

You may also specify alternative isotopes for BOMD jobs using the standard method.

All semi-empirical, SCF, CASSCF, CIS, MP2 and DFT methods.

ADMP

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 short 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 3.0 or other visualization software to display the trajectory path in three dimensions.