Interaction energy of O2 dimers - script doubt

Dear colleagues,

I am calculating the interaction energy of O2 dimers and the result values ​​are “sounding strange” to me (based on energy values ​​I got for other dimers).

Searching the forum, I found this post

From what I noticed the problem would be related to the representation of oxygen (ground state perhaps). I confess that I did not understand it very well, because I am not from the area.

Could someone more experienced tell me if my script is correct or how to fix it?


molecule oxygen_dimer{
   O 0.000000000 0.000000000 0.621515000
   O 0.000000000 0.000000000 -0.621515000
   O 0.310757500 2.379657500 0.439477471
   O -0.310757500 1.758142500 -0.439477471

set basis aug-cc-pVTZ

E = energy('MP2',molecule=oxygen_dimer, bsse_type='cp')
Efinal = E* psi_hartree2kcalmol
psi4.print_out("MP2/aug-cc-pVTZ = ")
psi4.print_out("%10.6f" % (Efinal))

This is the output I got for the input example above:

Thanks in advance for your time.

Your script is indeed incorrect.

For this to make any sense, you need to know a bit about spin in quantum mechanics. Molecules have a electron spin state (quantum number S) and a projection of the electron spin state on the z-axis (quantum number M_S). Psi assumes that M_S is as high as possible for a given S.

Normally, all electrons are paired up. That locks S = M_S = 0, and there’s no need to think about spin at all. This is not the case for molecular oxygen. In its ground electronic state, you have S = 1. Depending on how the electrons pairs up, you have M_S = 1, 0, or -1. That, we can handle. Psi determines how many unpaired electrons a molecule has by looking at the user-specified spin multiplicity. For a neutral molecule with two unpaired electrons, so you need to tell Psi 0 3.

The dimer is much trickier. The rules for how spins add tell us that we could have 0, 1, or 2, and those three different choices have different energies. Choosing S = 2 is the one possibility that can be well-represented with single-reference methods, and sadly, that is not the ground state.

You can either target the S = 2 excited state, drop O2 dimer from your study, or pursue multireference methods.

Thank you very much for the reply and the time dedicated to the explanation.

I graduated from a course where quantum mechanics didn’t go in depth and found myself working with some stuff in the field at the graduate level.

I will read the indicated paper.

Best regards.

Explaining what I’m trying to do: my idea is to obtain Lennard-Jones potential parameters from ab initio calculations. That’s why I thought of molecular O2.

So I would have to calculate the interaction energy of various molecular O2 conformers to fit the parameters of the interatomic potential curve.

To conclude, a final question arose: so I have to think about how many unpaired electrons the molecule has and pass this information on to psi4, correct?

Many thanks for your time!

It’s more complicated than that.

Let’s start with a simple MO-type model of O2. You have two electrons in the 2 pi antibonding orbitals. You can imagine up-up, up-down, down-up, and down-down. All of these states have two unpaired electrons. However, there is an S = 0 state (with the M_S = 0 substate up-down - down-up) and the S = 1 state (with the M_S = 1 substate up-up, the M_S = -1 substate down-down, and the M_S = 0 substate up-down + down-up). These all have the same number of unpaired electrons, but different S and M_S. We know from both experiment and theory that the three wavefunctions belonging to S = 1 are identical in energy, whereas the S = 0 state is considerably higher.[1] S = 0 and S = 1 represent entirely different potential energy surfaces, even though they have the same number of unpaired electrons.

If you assume that all unpaired electrons are the same spin, do you get the potential energy surface that you care about? If the answer is yes, then “think about how many unpaired electrons the molecule has and pass this information on to psi4”. In particular, take the number of unpaired electrons and add 1. For O2, this gives you 3, which is how I get to 0 3. If the answer is no, use a multireference method.

[1] = I’m ignoring relativistic effects. If relativity matters, the degeneracy may break. For O2, this should only matter for extremely high precision spectroscopy.