Accurate PES for H2 Molecule to investigate vibrational states

I used the basic code below to generate a PE curve for a H2 molecule. The results are promising with a De of 0.174 Ha and a far asymptotic value of -1.0 Ha. However, when I compute the number of bound vibrational states, I find 21 rather than the 14 expected from experimental measurements.

To improve the accuracy of the PE curve, should I do a fragment calculation for large values of R? Would the total spin be zero? Do the spin states become uncorrelated at large distances? How would I blend the fragment calculation with a molecule calculation?

thank you very much for your advice,
Don Burgess

import psi4
import numpy as np

mol_tmpl = """H
H 1 **R**"""

# rvals = np.linspace(0.11,0.50,40)
# rvals = np.append(rvals,0.55)

rvals = np.linspace(0.6,5.0,45)

fp = open('h2_pes.dat', 'w')

for r in rvals:
   mol = psi4.geometry(mol_tmpl.replace("**R**", str(r)))
   energy ="ccsd/aug-cc-pVQZ", molecule=mol)
   print(r,energy, file = fp)



The results are promising with a De = 0.174 Ha and a far asymptotic value of -1.0 Ha.

I don’t understand why you’re so concerned about dissociation. Yes, normally you would expect CCSD to not be reliable for dissociated systems. However, this is H2. For two-electron systems, CCSD is equivalent to full configuration interaction, so your description of electron correlation is just fine within the usual approximations (non-relativistic, Born-Oppenheimer)…

Accordingly, I’m not convinced that this is an issue with the potential curve.

I’m not sure how you determined the number of bound states, but are you sure that’s correct?

Because this is H2, transitions to some of those vibrational states might be spectroscopically forbidden by some symmetry selection rule. Unfortunately, I don’t have access to my notes on this to check my intuition.