Computational Methods in Chemistry

Exploring Potential Energy Surface of NH₃ and PH₃

In this excercise, we are going to learn how to use GAMESS quantum chemistry code, and to explore the inversion potential energy surface of NH₃ and PH₃. What we are interested in is how the equilibrium geometry and vibrational frequencies change with respect to the methods as well as the basis set.

We are going to crunch a lot of numbers. So, being systematic is important. I'd organize the input and log file (output file) in such a way that the name of an input file tells which molecule, method, basis set, optimization or hessian calculation. By naming, for example, a file as nh3.rhf.31.opt.inp refers to an input file for NH₃ the RHF/6-31G(d,p) optimization.

You need to do the following calculations:

Each calculation generates a .log file, an output file for the run, .irc and .dat files. The last file contains molecular orbitals, geometry, and gradient of the particular run. In this excercise, you need to look at .log file and .dat file.

When the optimization is finished, by searching for a word "equilib" in the log file, you can copy the unique atom geometry immediately below the word "equilib", and paste it into your input file for the Hessian run, as shown below.

After running a Hessian calculation, you should look at the output file and obtain information about vibrational frequencies and zero-point energy. It is shown below the cutout of the output of the Hessian run. The section is taken from the bottom of the output file. You would see 3 X N frequencies are shown, where N is a number of atoms in the molecule. In this example, the modes 1 to 6 are taken as rotational and translational coordinates. (Three coordinates to describe rotations, and three to translations, hence 3N - 6.) That means you read your vibrational frequencies from the 7th column. Here it starts with the 1142.21 cm^-1, and its IR spectral intensity is 5.13881. The rest of the 7th column include the displacement vectors for the vibrational motion associated with the 1142.21 cm^-1 frequency. As you can see, the motion of the 1142.21 cm ^-1 frequency includes movement of the N atom by + 0.10645190 Bohr, a ₀, in the Z direction, and the second H atom by 0.09379029 a₀ in the X direction, and so on so forth. This gives you the motion of atoms associated with the particular vibrational mode. It is much easier to understand how these motions look like by plotting on your monitor.

You have to create an input deck for plotting the molecular geometry and the vibrational frequencies on your screen (you have to have X window terminal or Mac). You have to edit .dat file. You have to cut everything between the "----- START OF NORMAL MODES FOR -MOLPLT- PROGRAM ----- and "----- END OF NORMAL MODES FOR -MOLPLT- PROGRAM -----", then paste this into xxx.mol file, where xxx is your unique file name. You also need to cut everything between the "-------- START OF -MOLPLT- INPUT FILE ----------" and "-------- END OF -MOLPLT- INPUT FILE ----------", and paste this at the top of the xxx.mol file.

The correct MOLPLT file look like this:

Preparing for Report

After you run all the calculations, you can compile the data. Write a report on how energetics change with respect to the basis set change as well as the mothodology change.

You should have gotten a transition state from the planar molecules. What is the criteria for a transition state? What are the barrier heights in different bases and methods?

You should also prepare the pictures of normal modes of vibration. What are the molecular motions found in NH₃ and PH₃? Are they different? If so, why? If not, why not? How does the imaginary frequency associated with the trasition state look like? Can you predict the product of this molecular motion?

Can you think of a reason or two that PH₃ has much higher barrier height than that of NH₃?

You can write a report based on these questions. Good luck!