Practical calculation of optimized geometries and frequencies using GAUSSIAN98

This is executed in Gaussian98 using the commands below in the input file for H₂O.

%chk=h2o_631g

# HF/6-31G opt

water

0 1

O 0.00000 0.00000 0.00000

H 0.75756 0.00000 -0.58639

H -0.75756 0.00000 -0.58639

--Link1--

%chk=h2o_631g

# HF/6-31G freq geom=checkpoint guess=checkpoint

calculate frequencies of water

0 1

The checkpoint file contains the force constants, polarizabilities, dipole, dipole derivatives, optimized geometry, and coefficients for the Self-Consistent Field calculation. The first phase of the calculation specified is a Hartree-Fock (HF) calculation with a 6-31G basis set. The opt keyword instructs the program to optimize the geometry. After the link (--Link1--) a frequency calculation is carried out (keyword freq). The optimized geometry and SCF coefficients are read from the checkpoint file (keywords geom=checkpoint guess=checkpoint).

To run the input file we use the script g98run.

g98run h2o_631g (generates the file h2o_631g.job)
qsub h2o_631g.job (submits the file h2o_631g.job to the queue)

Inspection of the output file using the vi editor reveals that number of steps required to converge. Since the input file is named h2o_631g.inp, the output file is named h2o_631g.out.

vi h2o_631g.out

In the vi editor type

/step <return>

to search for the word "step". We will jump to following place in the file. Notice that we can see the input geometry. The number of steps in the run is 20 (maximum).

----------------------------

! Initial Parameters !

! (Angstroms and Degrees) !

------------------------ -------------------------

! Name Definition Value Derivative Info. !

-----------------------------------------------------------------------------

! R1 R(1,2) 0.958 estimate D2E/DX2 !

! R2 R(1,3) 0.958 estimate D2E/DX2 !

! A1 A(2,1,3) 104.5166 estimate D2E/DX2 !

-----------------------------------------------------------------------------

Trust Radius=3.00D-01 FncErr=1.00D-07 GrdErr=1.00D-07

Number of steps in this run= 20 maximum allowed number of steps= 100.

GradGradGradGradGradGradGradGradGradGradGradGradGradGradGradGradGradGrad

Input orientation:

---------------------------------------------------------------------

Center Atomic Atomic Coordinates (Angstroms)

Number Number Type X Y Z

---------------------------------------------------------------------

1 8 0 0.000000 0.000000 0.000000

2 1 0 0.757560 0.000000 -0.586390

3 1 0 -0.757560 0.000000 -0.586390

---------------------------------------------------------------------

Distance matrix (angstroms):

1 2 3

1 O 0.000000

2 H 0.957993 0.000000

3 H 0.957993 1.515120 0.000000

It is very important to check the geometry and make sure it makes chemical sense. Here the distance matrix shows that the O-H bond distance is 0.957 Å. This is reasonable.

It is also important to realize that the computation will be done in the Standard orientation and NOT the Input orientation. We can search for the Standard orientation using the command

/Standard <return>

The screen will jump to

Standard orientation:

---------------------------------------------------------------------

Center Atomic Atomic Coordinates (Angstroms)

Number Number Type X Y Z

---------------------------------------------------------------------

1 8 0 0.000000 0.000000 0.117278

2 1 0 0.000000 0.757560 -0.469112

3 1 0 0.000000 -0.757560 -0.469112

---------------------------------------------------------------------

The standard orientation is rotated by 90^o.

To examine the number of steps required to reach convergence we just search for the word "Step". Type

/Step <return>

Note that the vi editor is case sensitive. We will jump to the following place in the output file

Step number 1 out of a maximum of 20

All quantities printed in internal units (Hartrees-Bohrs-Radians)

Second derivative matrix not updated -- first step.

The second derivative matrix:

R1 R2 A1

R1 0.55907

R2 0.00000 0.55907

A1 0.00000 0.00000 0.16000

Eigenvalues --- 0.16000 0.55907 0.55907

RFO step: Lambda=-2.96700106D-03.

Linear search not attempted -- first point.

Iteration 1 RMS(Cart)= 0.07296121 RMS(Int)= 0.00330803

Iteration 2 RMS(Cart)= 0.00299374 RMS(Int)= 0.00000025

Iteration 3 RMS(Cart)= 0.00000023 RMS(Int)= 0.00000000

Variable Old X -DE/DX Delta X Delta X Delta X New X

R1 1.81034 -0.00406 0.00000 -0.00722 -0.00722 1.80312

R2 1.81034 -0.00406 0.00000 -0.00722 -0.00722 1.80312

A1 1.82416 0.02177 0.00000 0.13359 0.13359 1.95775

Item Value Threshold Converged?

Maximum Force 0.021771 0.000450 NO

RMS Force 0.012999 0.000300 NO

Maximum Displacement 0.068035 0.001800 NO

RMS Displacement 0.071496 0.001200 NO

This part of the file tells we how the geometry will be modified and whether convergence has been reached. Notice that the derivatives and second derivatives are calculated in atomic units (1 Hartrees = 13.16 eV, 1 Bohr = 0.529 Å). We can continue to search for the word "Step" or we can jump to the optimized geometry. Search for

/Optimized

----------------------------

! Optimized Parameters !

! (Angstroms and Degrees) !

------------------------ -------------------------

! Name Definition Value Derivative Info. !

-----------------------------------------------------------------------------

! R1 R(1,2) 0.9494 -DE/DX = 0.0002 !

! R2 R(1,3) 0.9494 -DE/DX = 0.0002 !

! A1 A(2,1,3) 111.5599 -DE/DX = 0. !

-----------------------------------------------------------------------------

GradGradGradGradGradGradGradGradGradGradGradGradGradGradGradGradGradGrad

Input orientation:

---------------------------------------------------------------------

Center Atomic Atomic Coordinates (Angstroms)

Number Number Type X Y Z

---------------------------------------------------------------------

1 8 0 0.000000 0.000000 0.106786

2 1 0 0.785059 0.000000 -0.427142

3 1 0 -0.785059 0.000000 -0.427142

---------------------------------------------------------------------

Distance matrix (angstroms):

1 2 3

1 O 0.000000

2 H 0.949420 0.000000

3 H 0.949420 1.570119 0.000000

These are the optimized coordinates. For the purpose of calculating a B-matrix we need the Standard orientation, however, and so we now search on

/Standard

(or just scroll down using the j command in vi).

Standard orientation:

---------------------------------------------------------------------

Center Atomic Atomic Coordinates (Angstroms)

Number Number Type X Y Z

---------------------------------------------------------------------

1 8 0 0.000000 0.000000 0.106786

2 1 0 0.000000 0.785059 -0.427142

3 1 0 0.000000 -0.785059 -0.427142

---------------------------------------------------------------------

These are the coordinates we need. To write these to a separate file use the following commands.

:set number <return>

Numbers will appear in the left margin. In the present case they will have the appearance

516 Standard orientation:

517 ---------------------------------------------------------------------

518 Center Atomic Atomic Coordinates (Angstroms)

519 Number Number Type X Y Z

520 ---------------------------------------------------------------------

521 1 8 0 0.000000 0.000000 0.106786

522 2 1 0 0.000000 0.785059 -0.427142

523 3 1 0 0.000000 -0.785059 -0.427142

524 ---------------------------------------------------------------------

Now we use

:521,523w h2o_631g.coord

This will write the lines from 521 to 523 to an output file called h2o_631g.coord.

To check the calculated values of the frequencies search for the word "frequencies"

/freq

Harmonic frequencies (cm**-1), IR intensities (KM/Mole),

Raman scattering activities (A**4/AMU), Raman depolarization ratios,

reduced masses (AMU), force constants (mDyne/A) and normal coordinates:

1 2 3

A1 A1 B2

Frequencies -- 1736.2565 3991.4190 4148.5119

Red. masses -- 1.0915 1.0370 1.0888

Frc consts -- 1.9387 9.7342 11.0399

IR Inten -- 123.1789 2.9689 54.4346

Raman Activ -- 10.6127 90.0276 40.0303

Depolar -- 0.3960 0.2178 0.7500

Atom AN X Y Z X Y Z X Y Z

1 8 0.00 0.00 0.07 0.00 0.00 0.04 0.00 0.07 0.00

2 1 0.00 -0.38 -0.59 0.00 0.61 -0.35 0.00 -0.58 0.40

3 1 0.00 0.38 -0.59 0.00 -0.61 -0.35 0.00 -0.58 -0.40

The calculated frequencies are not in very good agreement with experiment for this calculation.

We will need to obtain the force constant and dipole derivatives.

These are found in the h2o_631g.chk file, however, this is a binary file. To convert the checkpoint file to ASCII we need to use the following commands

setenv g98root /usr/local/apps/gaussian
source /usr/local/apps/gaussian/g98/bsd/g98.login
formchk h2o_631g.chk

The new file will have the name h2o_631g.fchk.

Using the vi editor we can search for the "Cartesian Force Constants" in this file. We can identify the line numbers using

:set number

and then write these to an output file using the command

:[beg],[end]w h2o_631g.cmat

where [beg] is the number of the beginning line and [end] is the number of the ending line containing the Cartesian force constants.

Likewise, we can find the dipole derivatives (which will be capitalized

"Dipole Derivatives") and write these to an output file called h2o_631g.dd. This file will be useful for the calculation of infrared intensities.