I want to reproduce the bond dissociation energy of the BH molecule as stated in the book 'Molecular Electronic Structure Theory' by Helgaker et al, (chapter 8.3.1), using FCI with 1 frozen core orbital and the basis cc-pCVTZ.
Around the equilibrium bond length (at the experimentally determined one, 2.3289a0) I can reproduce the energy for all digits, so the error is below 1 microHartree, see first block of code below.
The dissociation energy given in the book is then computed w.r.t to the minimum of the computed energies, which the book says is at 2.3347a0. I compute the corresponding energy in the second block of code.
In the third and fourth blocks of code I compute the energy of the dissociated molecule, both by placing the atoms very far apart and by adding up the energies of the individual atoms.
The dissociation energies computed both ways differ only by microHartrees, but they are about 2 milliHartree larger than the reference value from the book. This is not a lot, but it hints at an error or misunderstanding in my calculations that I would like to understand.
Any ideas or hints are greatly appreciated!
from pyscf import scf, gto
from pyscf.mcscf import CASCI, UCASCI
basis_BH = {"B": "ccpcvtz", "H": "ccpvtz"}
# at experimental bond length, table 8.6
xyz_exp = "B 0 0 0; H 2.3289 0 0"
mol_exp = gto.M(atom=xyz_exp, basis=basis_BH, spin=0, symmetry=True, unit="B")
hf_exp = scf.RHF(mol_exp).newton()
hf_exp.run()
cascisolver_exp = CASCI(hf_exp, 56, 4)
e_r_exp = cascisolver_exp.kernel()[0]
e_r_exp_ref = -25.130001 - 0.101723 # 2nd column, 1st + 3rd entries
print(f"around equilibrium | computed: {e_r_exp}, reference: {e_r_exp_ref}\n")
# at equilibrium bond length, table 8.7
xyz_eq = "B 0 0 0; H 2.3347 0 0"
mol_eq = gto.M(atom=xyz_eq, basis=basis_BH, spin=0, symmetry=True, unit="B")
hf_eq = scf.RHF(mol_eq).newton()
hf_eq.run()
cascisolver_eq = CASCI(hf_eq, 56, 4)
e_eq = cascisolver_eq.kernel()[0]
# at large separation, dissociated molecule
xyz_sep = "B 0 0 0; H 1000 0 0"
mol_sep = gto.M(atom=xyz_sep, basis=basis_BH, spin=0, symmetry=True, unit="B")
hf_sep = scf.RHF(mol_sep).newton()
hf_sep.run()
cascisolver_sep = CASCI(hf_sep, 56, 4)
e_sep = cascisolver_sep.kernel()[0]
# compute energy of B (frozen core FCI) and H (FCI) separately (supermolecular approach)
molB = gto.M(atom="B 0 0 0", basis=basis_BH["B"], spin=1, verbose=0, symmetry=True)
uhf_B = scf.UHF(molB).newton()
uhf_B.run()
casci_B = UCASCI(uhf_B, nelecas=(2, 1), ncas=42)
e_B_fci_frozencore = casci_B.kernel()[0]
molH = gto.M(atom="H 0 0 0", basis=basis_BH["H"], spin=1, verbose=0, symmetry=True)
uhf_H = scf.UHF(molH).newton()
uhf_H.run()
casci_H = UCASCI(uhf_H, nelecas=(1, 0), ncas=14)
e_H_fci_frozencore = casci_H.kernel()[0]
e_BH_additive = e_H_fci_frozencore + e_B_fci_frozencore
# dissociation energy
e_diss_ref = 0.13267 # table 8.7
e_diss = e_sep - e_eq
e_diss_supermol = e_BH_additive - e_eq
print(
f"Dissociation energy | direct: {e_diss}, supermolecular: {e_diss_supermol}, reference: {e_diss_ref}"
)
the output is
around equilibrium | computed: -25.231724240266033, reference: -25.231724
Dissociation energy | direct: 0.13285929083973969, supermolecular: 0.13285247917572462, reference: 0.13267
```
