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 = psi4.energy("ccsd/aug-cc-pVQZ", molecule=mol)
print(r,energy, file = fp)
fp.close()
```

===================================================

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