Odd behaviour in diatomic molecule (Cl2) potential energy plot

Psi4Education
1

Hello,

I’m a newbie to Psi4, and indeed, to the field of computational chemistry. I believe in free (as in speech, not only as in beer) software, and a unfortunate thing I perceived when starting my graduate studies is that the field is still dominated by proprietary tools. So, first and foremost, I’d like to praise the Psi4 developers for the job you’re doing for science and all of us.

To learn Psi4, I started with the simplest possible system. A diatomic molecule. I wrote a small python script(PsiAPI) trying to plot potential energy versus bond length, first for hydrogen, and then for chlorine, using four different functionals (hf, scf, b3lyp and mp2) and four different basis sets (sto-3g, 3-21G, 6-311G and 6-311G(d,p)) combined:

'''
Plots E vs. d for diatomic chlorine. To change for hydrogen, you can replace 
both Cl's  in geometry section for H, change range(150) to range(50)
and internuclearD = 1.5 + n/100 to internuclearD = 0.5 + n/100.
'''
import psi4
import time

psi4.set_memory('4GB')

for method in ['hf', 'scf', 'b3lyp', 'mp2']:
    for base in ['sto-3g', '3-21G', '6-311G', '6-311G(d,p)']:
        for n in range(150):
            internuclearD = 1.5 + n/100
            #psi4.core.set_output_file('curvaPotencialCl2_' + str(n) + '.dat', False)
            complexo = psi4.geometry('''
            0 1
            Cl         0.00000        0.00000        0.00000
            Cl         0.00000        0.00000        ''' + str(internuclearD))
            try:
                t1 = time.time()
                energia = psi4.energy(method + '/' + base, molecule = complexo)
                t2 = time.time()
                print('resultado;' + str(method) + ';' + str(base) + ';' + str(internuclearD) + ';' + str(energia) + ';' + str(t2 - t1))
            except:
                print('resultado;' + str(method) + ';' + str(base) + ';;')

For hydrogen the results were just as expected. Energy varied a bit depending on method used, but all had the expected shape, a continuous curve with a minimum arond 0.7 angström:
H2All
For chlorine I got a strange result, for all functionals except b3lyp. When ploting the data, after de minimum the curve loses continuity, splitting in two, with large oscilations in energy between consecutive distance values.

Cl2
Cl2.png982×666 78.7 KB

What does cause this odd behaviour?

2

This is an awkward time for me, so I won’t be free to give a more detailed explanation until over a week. Looking at your data, this looks like a classic case of landing on the wrong electronic state. I don’t know your background, so I’m going to start from fundamentals.

We know that electronic states are eigenfunctions of the Hamiltonian, and all the methods approximate those electronic states. The Hamiltonian is a big matrix, so it has multiple eigenfunctions, and thus multiple electronic states for these methods to approximate. Around equilibrium, the methods all approximate the ground state fairly well, but as you leave equilibrium, the equations get more difficult to solve. Sometimes Psi can solve for the solution corresponding to the lowest energy state, and sometimes it solves for the next state up. That’s the cause of the discontinuities - different computations landing on different electronic states.

That also explains the pattern in the points. “SCF” is the name of a broad category of methods, but usually when people use it as a particular method, they mean “HF”, so the orange and green curves are the same. MP2 is an additive correction to HF, but the problem of landing on the wrong state occurs at the earlier SCF computation, so you expect the MP2 to have trouble where the HF does, which is exactly what you see.

As for what can be done to fix this…

  1. Longer term, we have plans to improve our SCF capabilities, which should hopefully avoid discontinuities like you see here.
  2. You can give the orbitals from the previous geometry as a guess for the next step by using set guess read.

By the way, the word “functional” usually means a Density Functional Theory functional. It is correct to call B3LYP a functional, but not MP2.

1 Like
3

An easy fix is to specify the “docc” option.

1 Like
4

Thank you for the explanation, and advice on the fix.

Now I have another doubt. Is there some complete reference on PsiAPI specification? For example, in Psithon input guide I’ve found how to mix basis sets:

basis {
   assign DZ
   assign C 3-21G
   assign H1 sto-3g
   assign C1 sto-3g
}

But I could not find a way to do the same in PsiAPI.

I have some basic python knowledge from before using Psi4, so I feel a bit more comfortable using the syntax I’m already familiar with. I’m trying to model some metal complexes. The common approach in the literature is to give the metal center a pseudopotential, and use different basis sets for lighter elements.

5

The psiAPI is listed here: http://psicode.org/psi4manual/master/psi4api.html
But often with just bare bone description.

Also check the /samples directory for examples. There should be various on basis sets, e.g. https://github.com/psi4/psi4/blob/master/samples/python/mints2/input.py
https://github.com/psi4/psi4/blob/master/samples/python/mints9/input.py

1 Like
6

Sounds similar to the HF/SCF discussion in this thread:

7

I gave it a cursory glance, and that also seems to be a case of SCF converging to a different electronic state. Do you have a question, or did you just want to make the connection?

8

Just making the connection…thanks for following up.