Resgrp:comp-photo/sh

Introduction

This wiki deals with construction of Gaussian routes for trajectory calculations, either with or without a surface hop. Trajectories are obtained by integrating the Newtonian equations of motion for nuclei over a sequence of steps on a potential energy surface (PES). This method does not require complete computation of the potential energy surface, but rather uses the local gradient and hessian to integrate over a small region. Furthermore a fifth-order fit refinement of the PES allows larger step sizes than have previously been possible, and hence speeds up total calculation time considerably.

If we want to investigate photochemical reactions, the trajectory must at some point hop from an excited state, higher energy surface to the ground state (and lower energy) surface. Three types of 'surface hopping algorithm' have been implemented, and instructions for using these will be discussed.

The Born Oppenheimer Molecular Dynamics (BOMD) algorithm

The BOMD keyword in GAUSSIAN requests a classical trajectory calculation [ref] using a Born-Oppenheimer molecular dynamics model [ref]. Some of the current features are:

• Predictor step on the local quadratic surface
• Corrector step using a fifth-order polynomial or a rational function fitted to the energy
• Gradient and Hessian evaluations at the beginning and end points of each step
• Selection of the initial conditions using quasi-classical fixed normal mode sampling
• Final product analysis carried out in the same manner as in the classical trajectory program VENUS
• Initial Cartesian coordinates and velocities may be read in as well

Preparing a Trajectory Job

Example 1: test379

Example of BOMD job using the surface hopping algorithm:

#P CAS(3,3,NRoot=2)/6-31G* test
# NoSymm Pop=Full scf=tight
# BOMD=(NTraj=1,StepSize=300,ReadVelocity,MaxPoints=6) units=au
# iop(1/44=1,1/7=11) iop(1/80=1,10/80=1) iop(5/97=11,10/97=11,10/99=300)
# iop(5/17=41000200,10/10=700007)

• IOp(1/44=1): set seed for random number generator to 1.
• IOp(1/7=11): Project out angular force option. No projection in correction step, No projection in prediction step.
• IOp(5/17=41000200,10/10=700007): includes the CP-MC-SCF correction and specifies state averaged orbitals.
• IOp(5/97=11,10/97=11): Options for Trajectory Surface Hopping calculations

Method 2:

     97=11...Make a hop based on the secular equation (adiabatic hop) this option also includes a hop decision based on the vector following method (diabatic hop) 
    (Needs option 99 and 80, see below)

             These options must be set in multiple links:
                          L1003   L510  L118
             iop(97)      yes     yes    no
             iop(55-59)   yes     yes    no
             iop(80)      yes     no     yes

             These options must be set for the following links:
                     55-58 80   97   98   99
             l118    no    yes  no   no   no
             l510    yes   no   yes  no   no
             l1003   yes  yes  yes  yes  yes

• IOp(10/99=300):
• IOp(1/80=1,10/80=1): surpress the 5th order correction after surface hop has been made in Trajectory Surface Hopping calculations.

Sampling

Ground state sampling

In the first part of the job, the sampling is conducted on the S0 PES after Force-Constant (FC) being computed. This provides the starting coordinates and velocities to be used in the second part of the job.

In the second part of the job, the dynamics in the S1 PES is carried out. The coordinates are read from the checkpoint file. However, the velocities can not be read from the checkpoint file (for technical reasons), which forces the user to copy-paste by hand the starting velocities. Therefore, the first part of the job has to be executed once. Then, after the user has to copy the velocity and run the two parts together. Note that the first part has to be run again anyway in order to get the initial structure after sampling.

%Chk={path}/{name}.chk
#P CAS(3,3,NRoot=1)/6-31G* test
# NoSymm Pop=Full scf=tight
# BOMD=(StepSize=300,maxpoint=1,sample=fixed) units=au
# iop(1/44=1,1/7=11) iop(1/80=1,10/80=1) iop(5/97=11,10/97=11,10/99=300)
# iop(5/17=41000200,10/10=700007)

[CHARGE, SPIN]
[STRUCTURE]

0

--Link1--
%Chk={path}/{name}.chk
#P CAS(3,3,NRoot=2)/6-31G* test
# NoSymm Pop=Full scf=tight
# BOMD=(StepSize=300,ReadMWVelocity,MaxPoints=6) units=au geom=check
# iop(1/44=1,1/7=11) iop(1/80=1,10/80=1) iop(5/97=11,10/97=11,10/99=300)
# iop(5/17=41000200,10/10=700007)

[CHARGE, SPIN]

0

[VELOCITY]