2.3. A more complex system: H\(_{2}\)O#

We will now move to a bit more complex system, with more than one atom. We will focus on the geometry of the system and the effect of basis set used.

Again, set up the environment and import the required modules:

import psi4
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import sys
sys.path.append("..")
from helpers import *

plt.style.use(['seaborn-poster', 'seaborn-ticks'])

psi4.set_memory('2 GB')
psi4.set_num_threads(2)

  Memory set to   1.863 GiB by Python driver.
  Threads set to 2 by Python driver.

2.3.1. Building the Geometry#

In the previous example the geometry file was quite trivial since we had a single atom. Now we will need to set up the geometry for a slightly more complicated system: a water molecule.

Define the geometry for a water molecule in the next cell, making sure to use the right syntax. Remember that you can use either Z-matrix or Cartesian format (one of the two will make your life easier!) and that information about the total charge of the system and spin multiplicity also need to be provided. It may be useful to refer to the dropdown menu in the previous section about atomic coordinates representation.

h2o = psi4.geometry("""
## FILL HERE
""")

drawXYZ(h2o)

You appear to be running in JupyterLab (or JavaScript failed to load for some other reason). You need to install the 3dmol extension:
jupyter labextension install jupyterlab_3dmol

Exercise 9

Include in your report the coordinates you set for your water molecules and a screenshot of how it looks like in 3D representation.

Exercise 10

How do you calculate the spin multiplicity of a species? Compare the case of the water molecule to the previous example of the hydrogen atom.

Extending what we have seen with the hydrogen atom, we will now investigate the effect of the choice of the basis sets when we have heavy atoms (remember that in the context of computational chemsitry, we consider as heavy each atom other than hydrogen).

Exercise 11

Compute energy calculations for the different basis sets indicated below with a loop similar to the one used before. Include a table of the the calculated energies. Specify what is the difference between the basis sets that we used.

Bonus: which of the additional functions introduced with respect to the 6-31G basis set are more reasonable for the present case, i.e. a water molecule? Why?

basissets = ['6-31G', '6-311G', '6-31+G', '6-31G**', '6-31++G**'] # these are the basis sets we are going to use

psi4.set_options({'reference':'UHF'}) # We are using Unrestriced Hartee Fock

for basis in basissets:
    psi4.core.set_output_file(f'{basis}-h2o-output.log', False) # save in seperate log files
    E = None  # define here the energy calculation with psi4 for single point calculations

    print(basis, E)

Few things important to keep in mind:

  • As you have seen, and will see in the next exercises, the value of total energy is highly dependent on the method and the basis set. In the hydrogen example, we were lucky to have an exact solution to compare with, but most time it is useless to discuss about total energy! What we are usually interested in are energies changes for the same systems with different states, or binding energy.

  • All the calculations performed so far on the water molecule are referred to the geometry that you defined. By running such a single point energy calculation we find the lowest energy solution for the Schrödinger equation at the current geometry. This however does not guarantee it is the lowest-energy conformation, i.e. out calculations may refer to a geometry in practice never observed in nature. What we usually do is performing a geometry optimization, i.e. starting from a given geometry and arriving an optimal geometry. This is extremely important expecially when we are dealing with complex molecules, for which it is possible we don’t know how the molecule looks like, for example not knowing the equilibrium bond length or the angle between the atoms. This concept will be used a lot in this course, and one of the sections of the next exercise session will focus on geometry optimization.

2.4. Output Files#

So far we only focused on the energy given as ouput from the psi4.energy() function. However, for each calculation a .log file is produced, containing detailed information of the procedure.

You can list the files in your directory (by typing ls in your terminal) to see what new files have been generated. Indeed, you should find the .log files created by Psi4. As this files may be very large and we have no intention of editing it.

As example, you can use the command head followed by the name of the last .log file created (6-31++G*-output.log) to see the last lines of that file. By using the option -20, the last 20 lines will be printed (where you can use any number instad of 50). Make sure to specify the file name as 6-31++G\*\*-h2o-output.log, where the *needs to be preceded by a \ symbol. This is because otherwise in Linux syntax the * is a special character that will be interpreted as ‘zero or more characters’, i.e. by typing 6-31++G**-h2o-output.log the system will look for all files which have zero or more characters between 6-31++G and -h2o-output.log. However, in our case we want to use * as a proper character (which is the file name), that’s the reason of using \*.

!tail -20  6-31++G\*\*-h2o-output.log

  Electronic Dipole Moment: [e a0]
     X:     0.0000      Y:     0.0000      Z:    -0.0864

  Dipole Moment: [e a0]
     X:     0.0000      Y:     0.0000      Z:     0.9077     Total:     0.9077

  Dipole Moment: [D]
     X:     0.0000      Y:     0.0000      Z:     2.3071     Total:     2.3071


*** tstop() called on noto.epfl.ch at Wed Oct 12 16:46:05 2022
Module time:
	user time   =       0.65 seconds =       0.01 minutes
	system time =       0.03 seconds =       0.00 minutes
	total time  =          0 seconds =       0.00 minutes
Total time:
	user time   =      18.78 seconds =       0.31 minutes
	system time =       0.91 seconds =       0.02 minutes
	total time  =        493 seconds =       8.22 minutes

If you are interested in extracting the lines of a file containing one particular keyword another useful command to know is grep. For example, considering again the last .log file created, it is possible to extract the final energy with the command:

!grep 'Final Energy' 6-31++G\*\*-h2o-output.log

  @DF-UHF Final Energy:   -76.02966176454109

Exercise 12

Have a look at the log files produced so far and answer the following questions:

  1. What is the significance of the statement Energy and wave function converged?

  2. What is the meaning of the different iter preceding Energy and wave function converged? Compare the number of cycles for the different basis sets.