A more complex system: H_{2}O

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 *

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

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 remember 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)

Exercise 10

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

Exercise 11

How do you calculate the spin multiplicity of a species? Compare the case of the water molecule to the previous examples of the hydrogen atom and the helium 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 12

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. our 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. But for now we still can play with different geometries of the water molecule.

2.4. Linear water molecule and BeH\(_{2}\)#

It is assumed that the water molecule you created previously has the usual bent geometry we are familiar with. But would it be possible to compute the energy of a linear water molecule? Well, let’s figure it out. First, try to modify the water geometry so that it has a linear conformation:

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

drawXYZ(h2o_linear)

Funny shape right? Now we can use it as an input for our energy calculation:

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_linear-output.log', False) # save in seperate log files
    E = None  # define here the energy calculation with psi4 for single point calculations

    print(basis, E)

Exercise 13

What can you say about the relative energies and stability of a linear water molecule compared to our original bent case? And why?

And now we can reproduce the same procedure to gaseous beryllium hydride BeH\(_{2}\). Can you predict the results of the calculation?

beh2_bent = psi4.geometry("""
## FILL HERE
""")
beh2_linear = psi4.geometry("""
## FILL HERE
""")

drawXYZ(beh2_bent)
drawXYZ(beh2_linear)

We will compare only the relative energy using one specific basis set but feel free to try different ones if you are curious.

basis = '6-31++G**'  # these are the basis sets we are going to use

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

psi4.core.set_output_file(f'{basis}-beh2_bent-output.log', False) # save in seperate log files for the bent case
E_bent = None  # define here the energy calculation with psi4 for single point calculations for the bent case
psi4.core.set_output_file(f'{basis}-beh2_linear-output.log', False) # save in seperate log files for the linear case
E_linear = None  # define here the energy calculation with psi4 for single point calculations for the linear case

print("The bent molecule has the energy:")
print(basis, E_bent)
print("The linear molecule has the energy:")
print(basis, E_linear)

Exercise 14

What can you say about the stability of the bent and linear beryllium hydride conformations? How does it compare to the water case and why?

2.5. 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 these files may be very large and we have no intention of editing it.

As an example, you can use the command tail 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

     Multipole            Electronic (a.u.)      Nuclear  (a.u.)        Total (a.u.)
 ------------------------------------------------------------------------------------

 L = 1.  Multiply by 2.5417464519 to convert [e a0] to [Debye]
 Dipole X            :          0.0000000            0.0000000            0.0000000
 Dipole Y            :          0.0000000            0.0000000            0.0000000
 Dipole Z            :         -0.0825639            0.9840423            0.9014784
 Magnitude           :                                                    0.9014784

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

*** tstop() called on noto.epfl.ch at Mon Sep 22 19:33:56 2025
Module time:
	user time   =       0.29 seconds =       0.00 minutes
	system time =       0.01 seconds =       0.00 minutes
	total time  =          0 seconds =       0.00 minutes
Total time:
	user time   =       3.68 seconds =       0.06 minutes
	system time =       0.20 seconds =       0.00 minutes
	total time  =         11 seconds =       0.18 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.02974219694876

These commands allow us to extract important information from the output file that can be used for analysis for instance. We will usually have a look at the output files during the next exercises, but for now it is sufficient if you get a feeling of what kind of information there is inside it.

Contents