Gaussian: Start to Practice Computational Chemistry
I am going to start this post with a simple example of computational chemistry, in which the beginner should know about computational chemistry. When computational chemists want to study the electronic structure or chemical properties of molecules, they have to model and optimize the geometry of the molecule first. Then they continue to compute the electronic properties using a certain theoretical technique. In reality, we cannot compute the calculation by hand because the mathematical equation is complicated. Instead, we use the program or software to help us to perform the calculation.
Geometry Optimization
Potential Energy Surface (PES)
Molecular geometry optimization is the minimization of the energy of a molecule to find its stable form. We first consider the potential energy surface of the typical molecule as shown in Figure 1.
So the minimization of energy getting the coordinate for structure is assumed in Cartesian coordinate (XYZ format) as following equation 1.
Figure 1. Potential Energy Surface (PES)
Energy Minimization
The Born-Oppenheimer approximation is used to reduce the complexity cost of computation. The geometry of a molecule at zero absolute temperature (temp = 273.5 K) corresponds to the minimum of the total energy:
coordinate (x,y,z) = min (< \psi | \hat(H) | \psi(*) >) Eq 1.
where the wave function and its twin are conjugated with each other and H is the Hamiltonian operator. Therefore, the process of finding the coordinate that minimizes the energy is known in the trade as molecular geometry optimization, which is usually the first step in the calculation of the flowchart displayed in Figure 2.
We know that the steepest ascent direction of any pathway is called gradient which is 1st derivative of energy. Also, the 1st derivative of the gradient is Hessian we would consider it as a judgment to identify the state of the molecule.
G_i = d(E_i)/d(r_i) Eq 2.
H_ij = d(E_ij)/[d(r_i) d(r_j)] Eq 3.
If the gradient is equal to zero and with another condition that
if H can be evaluated
then All H eigenvalues positive => Minimum point
else All H eigenvalues negative => Maximum point
else H eigenvalues of both signs => Saddle point (transition state)
fi
Minimization method
Nelder-Mead simplex
Gradient descent: steepest ascent
Conjugated gradients
Newton-Raphson method
Figure 2. Geometry optimization flowchart
Gaussian Input Preparation
To do calculations, you have to give data to a program with the molecule and methods that you want to study and use. This initial information must be sufficient for computing. We call the code that can help us to calculate that "program", a file that includes the information of molecule that "input", and a file that include the results that "outout". Most computational chemistry program packages have the same idea and concept. In this post, I introduce how to prepare input for the Gaussian program, invented by John A. Pople. Gaussian program is the role model of many codes at present. Since it has high performance and provides high accuracy results in a range of comparable experimental results, the whole scientists around the world choose Gaussian as the main program to conduct the research.
Structure of Input File
The syntax or structure of Gaussian input is a friendly interface with the user even with who is novice in computational chemistry, and also can be written flexibly and easily adaptable. The way to create an input file is the following.
Syntax of input file: sample_of_input_file.inp
%chk=your_input_name.chk
%mem=M
%nprocshare=N
# Job task, method of calculation, and optional keyword are here
----- blank line -----
Comment is here, you can write anything on this line.
----- blank line -----
Charge and spin multiplicity of molecule
The first line of Cartesian coordinate (XYZ) of molecule
..
..
..
..
The last line of of Cartesian coordinate (XYZ) of molecule
The assignment of resource requirements is also important.
You have to define two things about machine performance.
1. How much memory do you want to use? = M
Example. %mem=800MW
%mem=5000MB
%mem=2GB
2. How many CPUs do you want to use? = N
Example. %nprocs=1 (single cores/serial method)
%nprocs=4 (run in parallel method)
* Note that the memory unit can be KB, MB, GB, TB, KW, MW, GW or TW.
Example of Gaussian Input File
File: test001.inp
Geometry optimization of H2O molecule.
%chk=test001.chk
%mem=1GB
%nprocs=4
# TEST OPT STO-3G OPTCYC=2 scf=conventional
Gaussian Test Job 10, STO-3G FLETCHER-POWELL OPTIMIZATION OF WATER
0 1
O 0.000 0.000 0.000
H 1.000 0.000 0.000
H 0.000 1.000 0.000
File: test002.inp
Single point energy calculation of O2 molecule.
%chk=test002.chk
%mem=500MW
%nprocs=1
#P TEST STO-3G COMPLEX pop=full scf=tight
Gaussian Test Job 01, SINGLET DELTA STO-3G//STO-3G DIOXYGEN
0 1
O
O 1 R
R 1.220
* You should leave a few blank lines at the end of the input file.
Specify different basis set on different atom
The following is an example of Gaussian input that assigns different basis sets on different atoms in the same calculation.
%chk=test1.chk
%mem=1GB
%nprocs=4
#P B3LYP/GEN gfinput IOP(6/7=3)
Title
<charge> <spin mult>
Cu 0.00000000 0.00000000 0.52583600
O 0.00000000 0.00000000 -1.41498800
H 0.00000000 0.80849900 -1.96466800
C 0.00000000 -0.80849900 -1.96466800
...
...
H 0.02000000 0.80849900 -1.96466800
C O H 0
6-311g(d)
****
Cu 0
lanl2dz
****
For calculation using effective core potential (ECP) basis set, GENECP keyword must be appended as well as specifying basis set at the bottom of input.
%chk=test1.chk
%mem=1GB
%nprocs=4
#P B3LYP/GEN pseudo=read gfinout IOP(6/7=3) ! same as #P B3LYP/GENECP gfinout IOP(6/7=3)
Title
<charge> <spin mult>
Cu 0.00000000 0.00000000 0.52583600
O 0.00000000 0.00000000 -1.41498800
H 0.00000000 0.80849900 -1.96466800
C 0.00000000 -0.80849900 -1.96466800
...
...
H 0.02000000 0.80849900 -1.96466800
C O H 0
6-311g(d)
****
Cu 0
lanl2dz
****
Cu 0 ! The ECP will be used for this atom
lanl2dz
Optimization output analysis
In this part, I computed the geometry optimization of formaldehyde using the Hartree-Fock (HF) method with an STO-3G basis set.
Input deck is following
File: ch2o.inp
%chk=ch2o.chk ! set path to save your checkpoint file.
%mem=1GB ! Total memory usage
%nprocs=2 ! Core
#P OPT HF/STO-3G scf=tight pop=regular ! Calculation set up
HF/STO-3G//HF/STO-3G formaldehyde ! comment
0 1 ! charge & multiplicity
C1 ! atom coordinate for whole molecule in z-format
O2 1 r2
H3 1 r3 2 a3
H4 1 r3 2 a3 3 180.0
r2=1.21672286
r3=1.10137241
a3=122.73666566
1. Beginning of output
Coordinate of molecule in Z-matrix format
Symbolic Z-matrix:
Charge = 0 Multiplicity = 1
C1
O2 1 r2
H3 1 r3 2 a3
H4 1 r3 2 a3 3 180. 0
Variables:
r2 1.21672
r3 1.10137
a3 122.73667
Coordinate of molecule in XYZ format
Input orientation:
---------------------------------------------------------------------
Center Atomic Atomic Coordinates (Angstroms)
Number Number Type X Y Z
---------------------------------------------------------------------
1 6 0 0.000000 0.000000 0.000000
2 8 0 0.000000 0.000000 1.216723
3 1 0 0.926436 0.000000 -0.595599
4 1 0 -0.926436 0.000000 -0.595599
---------------------------------------------------------------------
Chemical formula and symmetry
Stoichiometry CH2O
Framework group C2V[C2(CO),SGV(H2)]
Deg. of freedom 3
Full point group C2V NOp 4
Largest Abelian subgroup C2V NOp 4
Largest concise Abelian subgroup C2 NOp 2
The system is oriented such that the principal axis of the systems runs along the z-axis and that all atoms of the system are located in the yz-plane:
Standard orientation:
---------------------------------------------------------------------
Center Atomic Atomic Coordinates (Angstroms)
Number Number Type X Y Z
---------------------------------------------------------------------
1 6 0 0.000000 0.000000 -0.533912
2 8 0 0.000000 0.000000 0.682811
3 1 0 0.000000 0.926436 -1.129510
4 1 0 0.000000 -0.926436 -1.129510
---------------------------------------------------------------------
Detail in basis set. The number of basis functions is 12 functions.
Standard basis: STO-3G (5D, 7F)
Ernie: Thresh= 0.10000D-02 Tol= 0.10000D-05 Strict=F.
There are 7 symmetry adapted cartesian basis functions of A1 symmetry.
There are 0 symmetry adapted cartesian basis functions of A2 symmetry.
There are 2 symmetry adapted cartesian basis functions of B1 symmetry.
There are 3 symmetry adapted cartesian basis functions of B2 symmetry.
There are 7 symmetry adapted basis functions of A1 symmetry.
There are 0 symmetry adapted basis functions of A2 symmetry.
There are 2 symmetry adapted basis functions of B1 symmetry.
There are 3 symmetry adapted basis functions of B2 symmetry.
12 basis functions, 36 primitive gaussians, 12 cartesian basis functions
8 alpha electrons 8 beta electrons
Convergence criteria used for SCF calculation.
Requested convergence on RMS density matrix=1.00D-08 within 128 cycles.
Requested convergence on MAX density matrix=1.00D-06.
Requested convergence on energy=1.00D-06.
2. Middle of output
SCF calculation
Details in SCF cycle step 1. The energy of formaldehyde given in cycle is -112.337159234657 A.U.
Cycle 1 Pass 1 IDiag 1:
E= -112.337159234657
DIIS: error= 3.66D-02 at cycle 1 NSaved= 1.
NSaved= 1 IEnMin= 1 EnMin= -112.337159234657 IErMin= 1 ErrMin= 3.66D-02
ErrMax= 3.66D-02 0.00D+00 EMaxC= 1.00D-01 BMatC= 1.96D-02 BMatP= 1.96D-02
IDIUse=3 WtCom= 6.34D-01 WtEn= 3.66D-01
Coeff-Com: 0.100D+01
Coeff-En: 0.100D+01
Coeff: 0.100D+01
Gap= 0.609 Goal= None Shift= 0.000
GapD= 0.609 DampG=2.000 DampE=0.500 DampFc=1.0000 IDamp=-1.
RMSDP=1.34D-02 MaxDP=6.76D-02 OVMax= 0.00D+00
Lowest energy
When the SCF does converge, Gaussian will print the energy of optimized structure.
You can use 'grep' command in Linux to search and print you out with entire optimized energies.
SCF Done: E(RHF) = -112.354347141 A.U. after 9 cycles
NFock= 9 Conv=0.93D-08 -V/T= 2.0089
KE= 1.113601121083D+02 PE=-3.271252836584D+02 EE= 7.232359354491D+01
Molecular Orbital
In Gaussian and other quantum chemistry program, molecular orbitals are generally printed by default. There are only four virtual orbitals and only those are printed.
The orbital coefficients are given with respect to the molecule in its "Standard orientation" given in "Population analysis" deck in output file.
Mulliken charge analysis
Mulliken charges:
1
1 C 0.074976
2 O -0.187948
3 H 0.056486
4 H 0.056486
Convergence met !
Item Value Threshold Converged?
Maximum Force 0.000008 0.000450 YES
RMS Force 0.000005 0.000300 YES
Maximum Displacement 0.000019 0.001800 YES
RMS Displacement 0.000013 0.001200 YES
Optimized structural parameter*
The deck below is an example of truncated gaussian output that provides you the optimized parameters that respect to your optimized molecule.
The R, A, and D variables in each line are bond distance, atom angle, dihedral parameters, respectively,
-- Stationary point found.
----------------------------
! Optimized Parameters !
! (Angstroms and Degrees) !
-------------------------- --------------------------
! Name Definition Value Derivative Info. !
--------------------------------------------------------------------------------
! R1 R(1,2) 1.2167 -DE/DX = 0.0 !
! R2 R(1,3) 1.1014 -DE/DX = 0.0 !
! R3 R(1,4) 1.1014 -DE/DX = 0.0 !
! A1 A(2,1,3) 122.7367 -DE/DX = 0.0 !
! A2 A(2,1,4) 122.7367 -DE/DX = 0.0 !
! A3 A(3,1,4) 114.5267 -DE/DX = 0.0 !
! D1 D(2,1,4,3) 180.0 -DE/DX = 0.0 !
--------------------------------------------------------------------------------
Final structure in terms of initial Z-matrix:
C
O,1,r2
H,1,r3,2,a3
H,1,r3,2,a3,3,180.,0
Variables:
r2=1.21672286
r3=1.10137241
a3=122.73666566
Enthalpy energy
You can request Gaussian to compute thermochemistry parameter by adding "Freq" keyword.
If Freq was found, Gaussian will compute the vibrational frequency and relevant energies.
Thermochemistry results is at almost the end of output, in "Thermochemistry" block.
Temperature 298.150 Kelvin. Pressure 1.00000 Atm.
Zero-point correction= 0.031184 (Hartree/Particle)
Thermal correction to Energy= 0.034038
Thermal correction to Enthalpy= 0.034982
Thermal correction to Gibbs Free Energy= 0.010176
Sum of electronic and zero-point Energies= -112.323163
Sum of electronic and thermal Energies= -112.320309
Sum of electronic and thermal Enthalpies= -112.319365
Sum of electronic and thermal Free Energies= -112.344171
E (Thermal) CV S
KCal/Mol Cal/Mol-Kelvin Cal/Mol-Kelvin
Total 21.359 6.264 52.209
Electronic 0.000 0.000 0.000
Translational 0.889 2.981 36.130
Rotational 0.889 2.981 16.026
Vibrational 19.582 0.302 0.053
Q Log10(Q) Ln(Q)
Total Bot 0.208549D-04 -4.680791 -10.777919
Total V=0 0.460044D+10 9.662799 22.249418
Vib (Bot) 0.454918D-14 -14.342067 -33.023830
Vib (V=0) 0.100351D+01 0.001523 0.003507
Electronic 0.100000D+01 0.000000 0.000000
Translational 0.646199D+07 6.810366 15.681448
Rotational 0.709431D+03 2.850910 6.564463
3. Ending output file
At the end of an output file, you should see the quotation message and summary of job time without any error warning.
"Normal termination of Gaussian ..." phrase also confirms that your calculation is done successfully.
WAR ES EIN GOTT DER DIESE ZEICHEN SCHRIEB?
- LUDWIG BOLTZMANN, QUOTING GOETHE, ABOUT MAXWELL'S EQUATIONS.
Job cpu time: 0 days 0 hours 0 minutes 1.5 seconds.
File lengths (MBytes): RWF= 5 Int= 0 D2E= 0 Chk= 1 Scr= 1
Normal termination of Gaussian 09 at Thu May 24 20:36:37 2018.
References
Ilya Kuprov, University of Southampton, 2012
J. Nocedal, S.J. Wright, Numerical Optimization, Springer, 2006.
J.B. Foresman, A. Frisch, Exploring Chemistry with Electronic Structure Methods, Gaussian, 1996.
Rangsiman Ketkaew