The Mathematica Notebook He_Atom_PT1.nb (PDF He_Atom_PT1.pdf) illustrates how to obtain the first-order energy correction for the helium atom after analytically solving the Schrodinger equation for a zeroth-order Hamiltonian which neglects the electron-electron repulsion. It also demonstrates the capabilities of Mathematica to calculate fairly difficult integrals. Examine the file to solidify your understanding of ideas behind perturbation theory; notice that the first-order energy is rather far from the true energy. In the next part, you will use more elaborate calculations to obtain better energy for the helium atom 1S ground state.
You will be using a computational chemistry program Gaussian03 (or Gaussian 09) to perform some calculations. Gaussian is a widely used commercial computational chemistry program, and it is important that you learn well how to use it. One valuable resource is their on-line technical documentation. Later in the course you can decide when to use Gaussian (e.g. geometry optimization with analytical derivatives) and when to use Firefly/PC GAMESS (e.g. single point calculations with s,p,d,f,g functions).
Similar to most quantum chemistry programs, Gaussian input consists of a text file that contains the instructions for performing the calculation, and the description of molecular geometry. A typical Gaussian input file will have a section that allocates computational resources (Link 0 commands that start with %) followed by the route section (lines that start with #), the title section (free format text), and molecule specification. There are three blank lines: one separating the route section from the title section, the other separating the title section from the molecule specifications, and the last one at the end of the molecule specification. The section that you added tells the Gaussian program to allocate 16 megawords (128 MB) of memory and perform optimization while evaluating the Hessian at each point. It also specifies the charge (0) and spin multiplicity (1) for the molecule. Both the charge and the spin multiplicity are integers. Charge is the total charge of the system (e.g. zero for a zwitterionic amino acid) and the multiplicity is 1 for closed-shell molecules in the ground electronic state. Multiplicity is 2 for single radicals and 3 for diradicals, such as molecular oxygen.
Using text editor of your choice (e.g. bluefish create an input file (He_SCF_STO3G.dat) for performing SCF calculation for the helium atom in the STO-3G basis. The minimal content of the file would be:
# HF/STO-3G He atom ground state (comment line) 0 1 He
Inspect your file. There should be an empty line at the end. Save the file as He_SCF_STO3G.dat in your directory. Execute the calculation by entering into the Unix shell the following command g09 < He_SCF_STO3G.dat > He_SCF_STO3G.out & . The calculation will perform a single-point Hartree-Fock energy calculation (there is no nuclear positions to potimize in atomic calculations) and because of small basis set, will finish the calculation quickly. You can easily check for the energy value by using the grep command: grep 'SCF Done' He_SCF_STO3G.out. Compare this energy value with the Hylleraas and first-order perturbation theory result you obtained earlier in the class.
When performing calculations on molecules, you need to input the initial geometry either as a set of XYZ coordinates of atoms, or as a connectivity matrix, known as the Z-matrix. For simple systems, such as diatomic or triatomic molecules, you can build Z-matrices by hand. For more complex molecules, it is easier to let molecular editors create Z-matrices. For example, you can use MOLDEN to build and save Gaussian Z-matrices. Open MOLDEN, build formaldehyde, and save it as aGaussian Z-matrix. Open the saved file, and insert the following block of lines just before the structure definition.
%Mem=48MW # HF/6-31+G(d,p) Opt=CalcAll MaxDisk=1GW Formaldehyde geometry optimization by calculating Hessian at every step 0 1
Inspect your file. There should be an empty line before and after the comment line, and at the end of the file.
Execute the calculation. You can call the program directly as in g09 < input.dat > output.out & .
The output of a typical ab initio calculation is quite daunting to novices. So many numbers, so many foreign-sounding terms, and then they repeat over nineteen times before the end of the file is reached ... and there is no answer in big bold red letters at the last line of the output. Gaussian at least emits reassuring messages from the past intellectuals, like this one:
"THE TEST OF A FIRST RATE INTELLIGENCE IS THE ABILITY TO HOLD TWO OPPOSED
IDEAS IN THE MIND AT THE SAME TIME, AND STILL RETAIN THE ABILITY TO FUNCTION.
ONE SHOULD, FOR EXAMPLE, BE ABLE TO SEE THAT THINGS ARE HOPELESS AND YET BE
DETERMINED TO MAKE THEM OTHERWISE."
-- F. SCOTT FITZGERALD
Most quantum chemistry calculations that you will be doing fall into one of the two categories: geometry optimization, and single point calculations. The calculation that you did with formaldehyde in the previous tutorial using Opt=CalcAll was an example of geometry optimization involving calculation of the Hessian at every step. In this calculation, a line such as "SCF Done: E(RHF) = -113.873389817 A.U. after 11 cycles" reports on the molecular energy in Hartree-Fock approximation at each step in optmization. You could use Unix command grep to see how the energy decreases during optimization:
grep 'SCF Done:' formaldehyde.out
SCF Done: E(RHF) = -113.873389817 A.U. after 11 cycles SCF Done: E(RHF) = -113.874554607 A.U. after 12 cycles SCF Done: E(RHF) = -113.874558296 A.U. after 12 cycles
To confirm that the calculation really converged, you would look also at the forces, which should become close to zero, and the size of structural changes from step to step (displacements), which should also be close to zero. These are reported at the end of each geometry optimization ctycle in a block like this:
Item Value Threshold Converged?
Maximum Force 0.033774 0.000450 NO
RMS Force 0.013621 0.000300 NO
Maximum Displacement 0.044624 0.001800 NO
RMS Displacement 0.022543 0.001200 NO
You can use grep 'Maximum Force ' formaldehyde.out to monitor the convergence. The convergence of energy, forces, and displacements is graphically shown by MOLDEN and other programs for the visualization of results from electronic structure calculations.