View source for Resgrp:comp-photo/sh/gdvb1-b6

= GAUSSIAN 99 B1-B6 development versions =

In the fist section we show how to generate the basic route for a trajectory job.  The next section explains all the trajectory options and how these options are to be combined. Finally the last section, explains all the hopping options.

References:

Original method to compute the trajectories without computing the complete potential energy surface:

* Helgaker, T; Uggerud, E.; Aa. Jensen, H.J. Chem. Phys. Letters 1990 173, 145-150
* Uggerud, E.; Helgaker, T. J. Am. Chem. Soc. 1992 114, 4265-4268 

Fifth order method, and the code that we use:

* Chen, W.; Hase, W.L.; Schlegel, H.B. Chem. Phys. Letters 1994 228, 435-442

Surface Hopping with the trust-region trajectories:

* Smith, B.R.; Bearpark, M.J.; Robb, M.A.; Bernardi, F.; Olivucci, M. Chem. Phys. Letters 1995 242, 27-32

== Constructing the Nonstandard Route ==

A non-standard route is generated by using the testrt alias (see your “Gaussian Users Manual” for more detailed instructions).  The route is usually specified on the command line (enclosed in quotation marks):

<pre>
$ testrt "# rhf/sto-3g"
</pre>

If the specified route is valid, testrt will print out the non-standard route corresponding to it.

The full command line required for a G99b1 trajectory:

<pre>
testrt "#P CAS(4,4,NROOT=2)/6-31G* TEST NOSYM POP=FULL SCF=TIGHT TRAJECTORY IOP(11/42=1) IOP(1/5=30,1/6=1,1/7=11,1/8=200,1/9=1,1/10=4,1/42=300,1/44=1)
IOP(5/17=41000200,10/10=700007) IOP(1/80=1,10/80=1) IOP(5/97=11,10/97=11,10/99=500)"
</pre>

<pre>
----------------------------------------------------------------------
#P cas(4,4,nroot=2)/STO-3G nosymm test trajectory pop=full scf=tight g
 uess=read iop(1/5=30,1/6=1,1/7=11,1/8=300,1/9=0,1/10=4,1/42=30,1/43=4,
 1/44=1) iop(11/42=1) iop(5/17=41000200,10/10=700007) iop(1/80=1,10/80=
 1) iop(5/97=22,10/97=22)
 ----------------------------------------------------------------------
 1/5=30,6=1,7=11,8=300,10=4,38=1,42=30,43=4,44=1,80=1/1,18;
 2/15=1,17=6,18=5/2;
 3/25=1/1,2,3;
 4/5=1,7=6,17=4,18=4/1,5;
 5/5=2,17=41000200,28=2,32=2,97=22/10;
 8/6=4,10=1,11=11/1;
 11/6=1,8=1,9=11,15=111,16=11,31=1,42=1/1,2,10;	
 10/6=1,10=700007,28=2,29=1,31=1,80=1,97=22/3;	
 6/7=3,28=1/1;
 7/10=1,25=1,30=1/1,2,3,16;
 1/5=30,6=1,7=11,8=300,10=4,42=30,43=4,44=1,80=1/18(3);
 2/15=1/2; 	
 7/9=1,25=1,30=1,44=-1/16;
 99//99;
 3/25=1/1,2,3;
 4/5=5,7=6,16=2,17=4,18=4/1;
 5/5=2,17=41000200,23=1,28=2,32=2,38=4,97=22/10;
 8/6=4,10=1,11=11/1;
 11/6=1,8=1,9=11,15=111,16=11,31=1,42=1/1,2,10;	
 10/6=1,10=700007,28=2,29=1,31=1,80=1,97=22/3; 	
 7/10=1,25=1,30=1/1,2,3,16;
 1/5=30,6=1,7=11,8=300,10=4,42=30,43=4,44=1,80=1/18(-7);
 6/7=3,19=2,28=1/1;
 7/9=1,25=1,30=1,44=-1/16;
 99/9=1/99;
</pre>

!!WARNING!! the testrt is often wrong, so pay careful attention the the following details and always check the unaltered parts of the route against jobs with routes that you are familar with.

== Required Changes ==

(a) remove all the calls to link 1102 and 1110.

(b) add (-3) to all the link 10 calls because the second call to L202 is jumped over (see 1/18(3)) there is no need to remove it.  The trajectory keyword ensures that the third call to L202 is removed.

<pre>
----------------------------------------------------------------------
#P cas(4,4,nroot=2)/STO-3G nosymm test trajectory pop=full scf=tight g
 uess=read iop(1/5=30,1/6=1,1/7=11,1/8=300,1/9=0,1/10=4,1/42=30,1/43=4,
 1/44=1) iop(11/42=1) iop(5/17=41000200,10/10=700007) iop(1/80=1,10/80=
 1) iop(5/97=22,10/97=22)
 ----------------------------------------------------------------------
 1/5=30,6=1,7=11,8=300,10=4,38=1,42=30,43=4,44=1,80=1/1,18;
 2/15=1,17=6,18=5/2;
 3/25=1/1,2,3;
 4/5=1,7=6,17=4,18=4/1,5;
 5/5=2,17=41000200,28=2,32=2,97=22/10;
 8/6=4,10=1,11=11/1;
 11/6=1,8=1,9=11,15=111,16=11,31=1,42=1/1,2,10;	remove the call to 2 and 10
 10/6=1,10=700007,28=2,29=1,31=1,80=1,97=22/3;	add (-3)
 6/7=3,28=1/1;
 7/10=1,25=1,30=1/1,2,3,16;
 1/5=30,6=1,7=11,8=300,10=4,42=30,43=4,44=1,80=1/18(3);
 2/15=1/2; 	
 7/9=1,25=1,30=1,44=-1/16;
 99//99;
 3/25=1/1,2,3;
 4/5=5,7=6,16=2,17=4,18=4/1;
 5/5=2,17=41000200,23=1,28=2,32=2,38=4,97=22/10;
 8/6=4,10=1,11=11/1;
 11/6=1,8=1,9=11,15=111,16=11,31=1,42=1/1,2,10;	remove the call to 2 and 10
 10/6=1,10=700007,28=2,29=1,31=1,80=1,97=22/3; 	add (-3)
 7/10=1,25=1,30=1/1,2,3,16;
 1/5=30,6=1,7=11,8=300,10=4,42=30,43=4,44=1,80=1/18(-7);
 6/7=3,19=2,28=1/1;
 7/9=1,25=1,30=1,44=-1/16;
 99/9=1/99;
</pre>

== Non-standard route explanation ==

=== Keywords ===

* #P: Additional output is generated. This includes messages at the beginning and end of each link giving assorted machine-dependent information (including execution timing data), as well as convergence information in the SCF.
* CAS(4,4,NRoot=2)/6-31G*: 4-electron, 4-orbital CASSCF calculation on the first excited state using a 6-31G* basis set
* NoSymm: Prevents the reorientation and causes all computations to be performed in the input orientation
* Test: This keyword suppresses the automatic creation of an archive entry (formerly intended for the Browse Quantum Chemistry Database System)
* Trajectory: Request a classical trajectory calculation
* Pop=Full: All orbitals are printed, along with the density matrices and a full (orbital by orbital and atom by atom) Mulliken population analysis. Since the size of the output depends on the square of the size of the molecule, it can become quite substantial for larger molecules. 
* SCF=Tight: Use normal, tight convergence in the SCF. This is the default.
* Guess=Read: Requests that the initial guess be read from the checkpoint file

=== IOps ===

* IOP(1/5=30): Read the kinetic energy for the N lowest modes
* IOp(1/6=1): Maximum number of steps 
* IOp(1/7=11): Convergence on the first derivative and estimated displacement for the optimization, RMS first derivative (ConvF = N*10**-6)
* IOp(1/8=200): Maximum step size allowed during Opt (DXMAXT = 0.01 * N)
* IOp(1/9=1): Generate random initial velocity
* IOp(1/10=4): Second derivative matrix calculated analytically
* IOp(1/42=300): N max steps taken by trajectory
* IOp(1/44=1): Set random number seed to N
* IOP(11/42=1): Force expanded form for compressed file formats. G99 versions  require that the gradient file information is not compressed, this is determined by this comand

== Non-standard route diagram ==

<blockquote>
<table border="1">

<tr><td> ''' BOMD Link(s)''' </td><td> '''General explanation''' </td><td> '''Detailed explanation''' </td> <td> '''Jump''' </td></tr>

<tr><td> L101,L118 </td><td> Read and parse route section. Initial entry: IOp(1/38=1), input data for L118 </td><td>Initializes program, controls overlaying and initialize a trajectory calculation </td> <td> </td></tr>

<tr><td> Freq (no L103 and L99 calls) </td><td> Frequency calculation </td><td> Cancellation of last printing of convergence criteria. Read MO from checkpoint file or Cards</td><td> </td></tr>

<tr><td> L101,L118 (3) </td><td> Read and parse route section. Continuation of run: IOp(1/38=0), start new traj. calc. </td><td> Controls the trajectory calculation and terminates the run if something goes wrong </td><td> Go to 1001 if it is an intermediate traj. step</td> </tr>

<tr><td> L202 </td><td> Determine molecular symmetry </td><td> Reorients coordinates, calculates symmetry and checks variables </td><td> </td> </tr>

<tr><td> L716 </td><td> Processes information for optimizations and frequencies </td><td> Converts cartesian forces and second derivatives to internal coordinates and communicates with optimization control programs </td><td> </td> </tr>

<tr><td> L99 </td><td> Save data to checkpoint file and calculation finalization </td><td> Terminates the run. Generates an archive entry and reformats output matrices </td><td> </td> </tr>

<tr><td> Freq (no L101, L601 and L99 calls) </td><td> Frequency calculation </td><td> Cancellation of population analysis and last printing of convergence criteria. Read MO from read-write file</td><td>1001 </td></tr>

<tr><td> L101,L118 (-7) </td><td> Read and parse route section. Continuation of run: IOp(1/38=0), traj. intermediate step</td><td> Controls the trajectory calculation and terminates the run if it is the last point </td> <td> Go to 1001 if it is an intermediate traj, step</td></tr>

<tr><td> L601 </td><td> Assign orbital and wavefuntion symmetries, print orbitals and perform Mulliken population analysis </td><td> Population and related analysis </td><td> </td> </tr>

<tr><td> L701,L702,L703,L716 </td><td> Processes information for optimizations and frequencies </td><td> 1-electron integral first or second derivatives, 2-electron integral first or second derivatives (sp). 2-electron integral first or second derivatives (spdf), </td><td> </td> </tr>

<tr><td> L99 </td><td> Save data to checkpoint file and calculation finalization </td><td> Terminates the run. Generates an archive entry and reformats output matrices </td><td> </td> </tr>

</table>
</blockquote>

<blockquote>
<table border="1">
<tr><td> '''Freq Link(s)''' </td><td> '''General explanation''' </td><td> '''Detailed explanation''' </td><td> '''Jump''' </td></tr>
<tr><td> L101 </td><td> Read and parse route section </td><td>Initializes program and controls overlaying. Reads Title, Z-matrix, Variables, Constants. </td><td> </td></tr>
<tr><td> L202 </td><td> Determine molecular symmetry </td><td> Reorients coordinates, calculates symmetry and checks variables </td><td> </td></tr>
<tr><td> L301,L302,L303 </td><td> Set up basis set </td><td> Generates bais set information, calculates overlap, kinetic and potential integrals, calculates multipole integrals </td><td> </td></tr>
<tr><td> L401,L405 </td><td> Generate initial orbitals </td><td> Forms the initial MO guess, initializes an MC-SCF calculation</td><td> </td></tr>
<tr><td> L510 </td><td> Mike Robb's MCSCF code </td><td> MC-SCF</td><td> 1002 </td></tr>
<tr><td> L801 </td><td> Produces molecular orbital coefficient matrix and eigenvalues </td><td> Initializes transformation of 2-electron integrals</td><td> </td></tr>
<tr><td> L1101 </td><td>  </td><td> Computes 1-electron integral derivatives</td><td> </td></tr>
<tr><td> L1003 (-3) </td><td> Computes whether a surface hop takes place or not </td><td>Iteratively solves the the CP-MCSCF equations. If a hop has occured the code must return to L510 and evaluate the gradient for the other state, otherwise continue. </td><td> Go to 1002 if a hop has taken place </td></tr>
<tr><td> L601 </td><td> Assign orbital and wavefuntion symmetries, print orbitals and perform Mulliken population analysis </td><td> Population and related analysis </td> <td> </td></tr>
<tr><td> L701,L702,L703,L716 </td><td>  Processes information for optimizations and frequencies </td><td> 1-electron integral first or second derivatives, 2-electron integral first or second derivatives (sp). 2-electron integral first or second derivatives (spdf) </td> <td> </td></tr>
<tr><td> L103 </td><td> Print convergence criteria </td><td> Controls "Berny" optimization </td> <td> </td></tr>
<tr><td> L99 </td><td> Save data to checkpoint file and calculation finalization </td><td> Terminates the run. Generates an archive entry and reformats output matrices </td><td> </td> </tr>
</table>
</blockquote>

== General guide of IOps ==

==== Trajectory Options ====

* IOP(5): Number of normal modes for which the initial kinetic energy is specified.
* IOP(6): Number of trajectories to compute, normally the value is 1.
* IOP(7): Option about angular forces, normally the value is 11
* IOP(8): The thrust radius, changes in this case it is the default of 0.3au (large step)
* IOP(9): Scale the inital sample, you don't want to, normally the value is 1
* IOP(10): Read the Hessian from CALCFC, normally the value is 4
* IOP(20): Use units of bohr, normally the value is 1
* IOP(38): either initial entry (1) or contonuation of run (0)
* IOP(42): Number of points along the reaction path, changes in this case is 20
* IOP(43): Correction scheme, normally the value is 4 (most accurate)
* IOP(44): Initialalize the random number generator, normally the value is 1


----

==== CASSCF Options ====

For the average trajectory on the excited state you will require that:

* IOP(5/17=41000200,10/10=700007): Includes the CP-MC-SCF correction and specifies state averaged orbitals

However, if you are doing a trajectory on the ground state with no state averaging and no hopping you require only that:

* IOP(5/17=200)

or if the job is very large:

* IOP(5/17=400)

Some options related to IOP(5/17) and IOP(10/10) are briefly outlined below:

* IOP(5/17) has many options, those of interest are ZY000X00

X=1 Use full 2nd order convergence

X=2 2nd order iteration at end, in preparation for CPMCSCF.

X=3 2nd order iteration using RFO type step + level shift.

X=4 Prepare for CPMCSCF(FREQ): Direct method with no storage of Hessian.

X=5 2nd order iteration using RFO type step+ level shiftand prepare for non-direct CPMCSCF

Y=1 determines that a full diagonalization method rather than Lanczos is used.

Z=1 Use State Average density matrices.(the weights 8F10.8)

Z=2 Do SA and prepare for SA-CPMCSCF.

Z=3 Do SA and prepare for Gradient of Energy difference.

Z=4 Do SA and prepare for SA Second Derivative Computation. (terms involving 2nd order orbital rotation derivatives not included)

* IOP(10/10) has many options, those of interest are: Y0000X

X=7 Computes the SA second derivatives. (The only approximation is neglect of the second order orbital rotation derivatives.)

Y=0 (Scaled Gradient Difference or Fxyz)

Y=1 Derivative coupling(without Division by Energy Diff.)

Y=3 Unscaled Gradient Difference * E2-E1.

Y=5 Read Forces from the input stream (Test purposes)

Y=6 Normalized Gradient Difference * E2-E1 + projected ivec

Y=7 Ivec state

Y=8 force (n-1) intersection search (to be used if GD is small)

----

==== L11 Options ====

I have also noticed that most of the other L11 options generated by the testrt command seem to be spurious, I will leave it up to your judgement if you choose to use them or not.  The required route information include only iops 31, 42, and 45:

* IOP(11/31=1) do not use symmetry in Rhys integral derivatives in L1110
* IOP(11/42=1) force exapanded form of file formats
* IOP(11/45=1) force NAt3 storage of matrices.

the remaining options are given below:

* IOP(11/6=1) contract the integral derivatives with the HF density matrix
* IOP(11/8=1) compute F1’s over AOs
* IOP(11/9=11) take the HF contributions from FX1 (A LA IFHFFX) when processing the two electron HF contributions
* IOP(11/15=111) control the output of derivatives to the rw-files, 100=calculate the one-electron contribution, 10=control output of “old” format, 1= force out-of-core algorithm
* IOP(11/16=11) compute the dipole derivative integral contribution to the HF dipole derivatives???

----

==== Turning on and off the SA ====

Linked to what kind of Hessian you calculate are the options for turning on and off the SA.  The entire computation of the trajectory can be carried out with coefficients 0.5.  However, this may not provide a good description of the system early or late in a trajectory, especially if the trajectory calculation is started some distance away from the conical intersection or it terminates some distance away from the conical intersection, ie where the two surfaces involved are no longer nearly degenerate.  In this case the state averaging can be gradually switched on or off.  The options for determining the state average weights are given in overlays 5 and 10, IOP(55) must be set in both links, the remaining iops can set only in link 5.  There are now two ways of doing state averaging, the old way, and a new way implemented in the G99c1 development version.

* IOP(55) the magnitude of change taken each step once SA is switched on/off in which case a negative value is specified.  This is useful when a trajectory is to be started inside a degenerate region.

0 (default) no action

N change by 1/N each step first change will be an increase from 0.0/1.0

-N change by 1/N each step first change will be a decrease from 0.5/0.5

* IOP(56) Switching on state averaging

0 (default) switch on when the energy difference is 0.050

N switch on when the energy difference is less than N*0.001

* IOP(57) Switching off state averaging

0 (default) switch off when the energy difference is 0.075

N switch on when the energy difference is greater than N*0.001

* IOP(58) When to do state averaging

1 (default) frequency calculation SA coefficients are between 0.0-1.0

2 Always SA second derivatives

* IOP(59) New method of changing SA coefficients

0 (default) use old inflexible method

1 use new flexible method

The old method switches SA on when you reach a certain energy gap and increases the SA coefficients each step after this point regardless of wether the energy gap continues to decrease or not.  This can cause a job sensitive to state averaging in a trajectory to fail to converge.  The new method switches on at a given energy and then checks the energy gap each cycle and each cycle decides wether to increase of decrease the SA coefficients.  See the hypothetical example with nroot=2 given below:

<blockquote>
<table border="1">
<tr><td> '''cycle''' </td> <td> '''old SA  ''' </td> <td> '''new SA''' </td><td> '''Egap''' </td></tr>
<tr><td> '''4 ''' </td> <td> 0.2/0.8 </td> <td> 0.2/0.8 </td> <td> <limit </td></tr>
<tr><td> '''5''' </td> <td> 0.3/0.7 </td> <td> 0.3/0.7 </td> <td> <limit </td></tr>
<tr><td> '''6''' </td> <td> 0.4/0.6 </td> <td> 0.2/0.8 </td> <td> >limit </td></tr>
<tr><td> '''7''' </td> <td> 0.5/0.5 </td> <td> 0.1/0.9 </td> <td> >limit </td></tr>
<tr><td> '''8''' </td> <td> crash </td> <td> 0.0/0.1 </td> <td> >limit </td></tr>
</table>
</blockquote>


----

==== Events after the hop ====

After a hop the fifth order correction to the PES calculated on the fly cannot be used because the gradient and hessian of the two steps do not belong to the same surface.  Hence to disable this correction for the step after a surface we use the route control options IOP(1/80=1) and IOP(10/80)=1. ALL hopping methods require this option.

Also after a hop, although the geometry of the molecule is the same, but is now on a different surface.  If it has hopped through an energy gap between the surfaces,  its total energy and the gradient must be calculated again because they will have changed.  This is the origing of the additional call to link L510 that is given in the route section, remember we added (-3) after all the calls to link 1003.

A surface hop needs to be performed with state averaging (SA) and at the hop point the coefficients must be 0.5 and 0.5.  The options for turning on and off state averaging over a trajectory have already been discussed

*  IOP(1/80=1,10/80=1): suppress the 5th order correction after surface hop has been made


----

==== Treating the Adiabatic Event ====

There are four ways of treating the adiabatic event:

* (1) Secular Equation Method
* (2) Vector Rotation Method
* (3) Time Dependent Vector Probabilities
* (4) Mixed State Propagation

	The Secular Equation Method deals with the large adiabatic coupling (H12) between the two states of interest that occurs near a crossing.  This method is good for determining when to an “adiabatic hop” ie to jump a small energy gap rather than going through the point of intersection.

	The Vector Rotation Method relys on the fact that the nature of the surface is quite different on either side of a “crossing”, as this epithet implies this method is good for determining a “diabatic hop” when the molecule passes through a point of intersection of the two important states without losing any energy.

	The Time Dependent Vector Probability method propagates a time dependent electronic vector along with the trajectory.  At the begining of the trajectory the probability of being in the upper state is 1.0 and as time goes on some population may leak to the ground state.  When the population change is essentially instantanious then the molecule has passed through a “crossing”.  When the population change is slow then a “hop” is determined when the population in the upper state reaches a certain pre-defined level, day 0.65.

	The Mixed State Propagation Method is not a classical approximation to a quantum mechanical event, as all the other methods are.  In this method a time dependent vector is propagated but this time not on a pure surface but on a mixed state that is determined by the population of the two important levels.  There will be no “hop” only a gradual change of the PES.

	The choice of method is determined by IOP(5/97 and the identical 10/97) However to complicate matters more options specific to each method need to be set in other links.

<pre>
*IOP(5/97): determines the hop method

10 Secular Equation Method only
11 Vector Crossing and Secular Equation Method
21 Testing option
22 Wavefunction Method
23 Time Dependent Gradient
</pre>

<pre>
* IOP(10/97): is the same as IOP(5/97) it determines the hop method

10 Secular Equation Method (Adiabatic)
11 Vector Crossing (Diabatic) and Secular Equation Method
21 Testing option
22 Wavefunction Method
23 Gradient Update with Wavefunction Method
</pre>

* IOp(5/97=11,10/97=11): Vector Crossing and Secular Equation Method
* IOp(10/99=500): threshold for the offdiagonal element in determing a hop, stop at N*0.0001

Options that must be set in multiple links:

IS IT RIGHT?

<blockquote>
<table border="1">
<tr><td> </td><td> '''L1003''' </td> <td> '''L510''' </td> <td> '''L118''' </td></tr>
<tr><td> '''IOp(97)''' </td> <td> yes </td> <td> yes </td> <td> no </td></tr>
<tr><td> '''IOp(55-59)''' </td> <td> yes </td> <td> yes </td> <td> no </td></tr>
<tr><td> '''IOp(80)''' </td> <td> yes </td> <td> yes </td> <td> yes (only if iop(97=23)) </td></tr>
</table>
</blockquote>

== Aditional features==

In the previous section we introduced the trajectory options without giving too much detail.  In this section the options avaliable to you will be fully discussed and an apendix of the various options is provided for quick reference.

=== Starting information ===

This information is most relevent if you wish start by pushing the molecule in a certain direction, for example you may want facilitate movement towards a conical intersection, or you may wish to excite a specific vibrational state leading to dissociation.  This is achieved by adding energy to specific vibrational modes.  You may also want to represent the ZPE already in the molecule when it gets excited to the first excited state.  This is achieved by randomly distributing energy in all the avaliable modes.  Alternatively, you may wish to restart a trajectory from a certain point on the trajectory, for example from a given geometry and velocity.

==== Vibrational excitation ====

A certan amount of initial kinetic energy can be placed directly into a mode this can be specified in kcal/mol or alternatively as the number of vibrational quanta to be excited.  The remaining modes are sampled and the displacement and velocity vectors associated with these modes randomly determined.  The total number of modes must, of course, be 3N-6.

* IOP(1/5) gives the first M lowest energy modes for which the starting energy is to be specified
* IOP(1/70) determines they type of sampling that is carried out on the remaining (3N-6)-M modes.

IOP(5) therefore moves the boundary between the number of modes for which a kinetic energy is specified and those that are sampled.  For example, lets assume you have a molecule with 6 vibrational mode, there are this various options for distributing energy in the molecule.  You can sample all of the modes (see option1 in the table below), or you can specify the energy of the first mode and sample all the other modes (option2), or you can specify the energy of all of the modes up to mode 4 and then sample the remaining two modes (option3), or you can specify the energy of all of the modes (option4).  Unfortunately you CANNOT sample the first mode and then specify the energy of any higher mode like mode 4 (option 5) [HOWEVER there is a fix for this problem].

<blockquote>
<table border="1">
<tr><td> '''Mode''' </td><td> '''Option 1''' </td> <td> '''Option 2'''  </td> <td> '''Option 3'''  </td><td> '''Option 4'''  </td><td> '''Option 5'''  </td></tr>
<tr><td> 1 </td> <td> Sample </td> <td> 1.0 kcal/mol </td> <td> 0.0 kcal/mol </td><td> 0.0 kcal/mol </td><td> Sample </td></tr>
<tr><td> 2 </td> <td> Sample </td> <td> Sample </td> <td> 0.0 kcal/mol </td><td> 0.0 kcal/mol </td><td> Sample </td></tr>
<tr><td> 3 </td> <td> Sample </td> <td> Sample </td> <td> 0.0 kcal/mol </td><td> 0.0 kcal/mol </td><td> Sample </td></tr>
<tr><td> 4 </td> <td> Sample </td> <td> Sample </td> <td> 1.0 kcal/mol </td><td> 0.0 kcal/mol </td><td> 1.0 kcal/mol </td></tr>
<tr><td> 5 </td> <td> Sample </td> <td> Sample </td> <td> Sample </td><td> 0.0 kcal/mol </td><td> Sample </td></tr>
<tr><td> 6 </td> <td> Sample </td> <td> Sample</td> <td> Sample</td><td> 0.0 kcal/mol </td><td> Sample </td></tr>
</table>
</blockquote>

<blockquote>
<table border="1">
<tr><td> '''IOp(5)''' </td><td> '''-1''' </td> <td> '''0'''  </td> <td> '''4'''  </td><td> '''6'''  </td><td> '''Impossible'''  </td></tr>
<tr><td> '''Mode''' </td><td> '''Option 1''' </td> <td> '''Option 2'''  </td> <td> '''Option 3'''  </td><td> '''Option 4'''  </td><td> '''Option 5'''  </td></tr>
<tr><td> 1 </td> <td>  </td> <td> 1 1.0 </td> <td> 1 1.0 </td><td> 1 0.0 </td><td> </td></tr>
<tr><td> 2 </td> <td>  </td> <td>  </td> <td> 2 0.0 </td><td> 2 0.0 </td><td> </td></tr>
<tr><td> 3 </td> <td>  </td> <td>  </td> <td> 3 0.0 </td><td> 3 0.0 </td><td> </td></tr>
<tr><td> 4 </td> <td>  </td> <td> </td> <td> 4 1.0 </td><td> 4 0.0 </td><td> </td></tr>
<tr><td> 5 </td> <td>  </td> <td>  </td> <td>  </td><td> 5 0.0 </td><td> </td></tr>
<tr><td> 6 </td> <td>  </td> <td> </td> <td> </td><td> 6 0.0  </td><td> </td></tr>
</table>
</blockquote>

	If you are going to use IOP(5) then the same geometry as that used to calculate the frequencies must be used.  If the molecule has symmetry the orientation of the molecule will be altered within the calculation, this can be avoided by using the NOSYMM keyword.  Always use the default option for IOP(1/70) I know nothing about this option!!

	Once you have decided how you are going to treat each vibation, for those that you are specifying the energy you must add this information into the com file and the format for this is as given below.  The vibration (IMode) and then the associated energy (VibEng) specified on the same line with at least one space between the numbers. Values are entered consecutively in the same order as the modes appear in a frequency calculation.  Multiple entries can run over more than one line.

Specification Format:

* Coordiantes
* (blank line)
* NPath
* IMode(1) VibEng(1) IMode(2) VibEng(2) IMode(3) VibEng(3) IMode(4) VibEng(4) IMode(5) VibEng(5) ...etc
* (blank line)

Note: If there is a blank line where the mode energies are to be specified, by default GAUSSIAN puts zero energy into those modes

	The default option, ie if you specify no IOP(5) option, is to read in a kinetic energy for the transition state mode and to sample all other modes, ie it is identical to IOP(1/5=1).  This is necessary because the original code was writen for trajectories always starting from a transition state, hence the energy of the single imaginary mode needed to be specified, all the other vibrational modes are sampled.  Whatever value you have for IOP(1/5=N) you must specify a value for every mode so defined, ie N modes (unless the energy is to be zero).  The maximium input is for IOP(5) is to specify energy for all 3N-6 modes ie IOP(1/5=3N-6).  Any modes you wish to give zero energy need not be specified.

	For example you may want to place the molecule on an excited state with no extra kinetic or vibrational energy in this case you set IOP(1/5=3N-6) and leave a single blank line at the end of the .com file (which will activate a default option to put zero energy into all the modes).

	Alternatively you may wish to determine an initital excited state structure and velocity for a molecule on the ground state subject to zero point energy vibrations.  In this case you set IOP(1/5=-1) and give no addition information at the end of the .com file.  This option can be used to find a statistical distribution of inital geometries and velocites for trajectories starting on the excited state that represent various starting points for molecules that have been instantaneously excited from a ground state subject to the effects of a zero point vibrational energy.

	If you give VibEng as a real (ie with a decimal place) the energy is assumed to be in kcal/mol.  If you specify an integer value (ie with no decimal place) then mode is excited by a given quanta (option4).  If VibEng is positive energy goes in the forward direction of the transition vector (+VibEng), if VibEng is negative energy goes in the reverse direction (-VibEng).  Alternatively the PHASE keyword can be used, but I know nothing about this option.  Input is terminated by a blank line.

	You can have some modes with a specific energy in kcal/mol and other modes with a specified quanta (option6), however you must specify the energy kcal/mol prior to the quanta (option7).  You can specify some modes in kcal/mol, others in quanta and require others to be sampled (option 8).  However, due to the original programming you cannot, for example, excite mode 3 by 2 vibrational quanta and then put 1.0 kcal/mol into mode 6 (option 9).

<blockquote>
<table border="1">
<tr><td> '''Mode''' </td><td> '''Option 6''' </td> <td> '''Option 7'''  </td> <td> '''Option 8'''  </td><td> '''Option 9'''  </td></tr>
<tr><td> 1 </td> <td> 2 quanta </td> <td> 1.0 kcal/mol </td> <td> 1.0 kcal/mol </td><td> 2 quanta </td></tr>
<tr><td> 2 </td> <td> 2 quanta </td> <td> 0.0 kcal/mol </td> <td> 0.0 kcal/mol </td><td> 2 quanta </td></tr>
<tr><td> 3 </td> <td> 2 quanta </td> <td> 0.0 kcal/mol </td> <td> 2 quanta </td><td> 2 quanta </td></tr>
<tr><td> 4 </td> <td> 2 quanta </td> <td> 2 quanta </td> <td> 2 quanta </td><td> 1.0 kcal/mol </td></tr>
<tr><td> 5 </td> <td> 2 quanta </td> <td> 2 quanta </td> <td> Sample </td><td> 0.0 kcal/mol </td></tr>
<tr><td> 6 </td> <td> 2 quanta </td> <td> 2 quanta </td> <td> Sample</td><td> 0.0 kcal/mol </td></tr>
</table>
</blockquote>

<blockquote>
<table border="1">
<tr><td> '''IOp(5)''' </td><td> '''6''' </td> <td> '''6'''  </td> <td> '''4'''  </td><td> '''Impossible'''  </td></tr>
<tr><td> '''Mode''' </td><td> '''Option 6''' </td> <td> '''Option 7'''  </td> <td> '''Option 8'''  </td><td> '''Option 9'''  </td></tr>
<tr><td> 1 </td> <td> 1 2  </td> <td> 1 1.0 </td> <td> 1 1.0 </td><td> </td></tr>
<tr><td> 2 </td> <td> 2 2  </td> <td> 2 0.0  </td> <td> 2 0.0 </td><td>  </td></tr>
<tr><td> 3 </td> <td> 3 2  </td> <td> 3 0.0  </td> <td> 3 2  </td><td>  </td></tr>
<tr><td> 4 </td> <td> 4 2  </td> <td> 4 2  </td> <td> 4 2  </td><td>  </td></tr>
<tr><td> 5 </td> <td> 5 2  </td> <td> 5 2  </td> <td>  </td><td>  </td></tr>
<tr><td> 6 </td> <td> 6 2  </td> <td> 6 2  </td> <td> </td><td>  </td></tr>
</table>
</blockquote>

==== Random number generation ====

When a trajectory job is started with energy in the normal modes this must be distributed into the potential energy, ie a geometry distortion, and the kinetic energy, ie an initial velocity, of the molecule.  When any kind of sampling is taking place the ZPE is randomly distributed equally between the PE and KE for a mode.  There are two random number generator IOPs for L118 IOP(1/38) and IOP(1/44).  IOP(1/44) determines the seed for the random number generator, then an array of (random) numbers is generated by a routine where IOP(1/38) determines the point in the array from which a (random) number is picked.  For two consecutive trajectory jobs to give identical results you require  both IOP(1/38) and IOP(44) to be identical.  If running a series of jobs from one starting point it is normally the case that IOP(1/44)=1 is set to a constant seed value and that IOP(1/38)=1,2,3...is varied, ie the position in the array from which the random number is taken is varied.

==== Initial velocity ====

It is possible to start a trajectory with a specific geometry and velocity, this option is most often used to restart a trajectory.  This option works in conjunction with IOP(1/20) or the UNITS keyword, and IOP(1/9), which determines the units of the geometric and velocities vectors.  Either you use IOP(1/20)=0 or UNITS=ANG, which determines the geometric units as Angstom and degrees and IOP(1/9)=3, which determines the velocity units as bohr/sec.  Or you use IOP(1/20)=1 or UNITS=AU, which determines the geometric units as Bohr and degrees IOP(1/9)=4, which determines the velocity units as amu**(1/2)*bohr/sec, ie the initial velocity of the molecule is specified in terms of mass weighted cartesian velocities for each atom.

Along with specifying the units of the input geometry, the velocity vectors themselves must be given.  These are entered after the geometry and stopping critera.  The x,y and z vector components of each velocity vector for each atom are then given.

Along with specifying the units of the input geometry, the velocity vectors themselves must be given.  These are entered after the geometry and stopping critera.  The x,y and z vector components of each velocity vector for each atom are then given.  For G98 input they are eneterd on the same line, while for a G99 input each component is on a new line, atom by atom, in both cases the velocities for each atom must be enetered in the SAME ORDER as in the geometry specification.  A starting velocity is required for each trajectory.  Hence, assuming you are running a single G98 trajectory, there should be atoms new lines each holding three real number seperated by a blank space.  While for a single G99 trajectory, there should be 3*atoms new lines each holding a real number. Input is terminated by an empty line.

Specification Format:

* Coordiantes
* (blank line)
* Npath
* (blank line)
* Atom Vx(1) Vy(1) Vz(1)
* Atom Vx(2) Vy(2) Vz(2)
* ...etc
* (blank line)

There has been some confusion about the units of the velocities given by the G98 output, this is because the units of the cartesian velocities are not indicated.  To ensure you have made the right choice check the first step velocities if you are restarting a job from velocites to ensure they are the same.  This especially applies to using G98 velocities to start a G99 job!

To the best of my knowledge the velocity units in the G98 output are mass weighted atomic units.  So in restarting either a G98 or G99 output the cartesian velocities you will be required to restart with IOP(1/20)=1 (bohr) and IOP(1/9)=4 (mass weighted bohr/sec).  A common mistake is to forget to change the geometric coordinates into bohr in the new input file!

=== Stopping information ===

There is no need to specify exact terminating criteria, the trajectory can just keep going until the maximum number of steps is reached.  In this case a blank card is entered.  Alternatively you can specify critera for when the trajectory is to stop, for example when it dissociates, or when it reaches the products or returns to the reactants.

* Maximium Number of Steps

	The maximium number of steps is given by IOP(1/42), the default for this is 50 steps.  Note that this is not a time critera, so if you choose to set a large step size 50 the trajectory will finish much later than if you set a small step size.  A typical number of steps for the default step size option is 50-300.  If you are not using any other stopping critera you need to put a zero on the line after the blank line after the coords, also finish with a blank line.

=== Restarting a trajectory ===

In general there are two ways to restart a trajectory:
(1) using the geometry and velocity from a particular step in the trajectory
(2) starting from a checkpoint file

	When restarting using a geometry and velocity always use the PREDICTED coords and velocity from a specific step, use this data as already outlined in the initial velocity section.

	If you wish to negate the input velocity vector then you use IOP(1/86)=1.  In the G99b1-b6 versions only you may not wish to take the large first step that normally initiates a trajectory and to avoid this you set IOP(1/89)=1, but be awear that this option may well comprimise energy conservation.  Using this method there is no restarting of the time dependent code, so you must restart outside of the hopping region.  You may also have already hopped in which case you need to swap the roots right from the begining, to do this set IOP(5/92)=xx1.

	Alternatively can restart directly from a checkpoint file.  First you add the IOP(1/35=1) then you remove the first call to link2 and then you alter or add any other options as you require.  For example if you have reached the maximium step allowance but want to continue you will need to increase IOP(1/42), otherwiswe the job will halt on the first step.  You can also alter other parameters such at the trajectory step size, timedep step size and hopping region limits.  Only avaliabe in G99c1, if the time dependent code (vector option only) was active at the finish then you need to set IOP(1/88)=1.  Note you cannot restart the time-dependent gradient option.

= Testjobs =