55 lines
2.1 KiB
Python
55 lines
2.1 KiB
Python
"""
|
|
H2 atomization energy test.
|
|
|
|
E_atom = E(H2) - 2 E(H) [energy to dissociate H2 into 2 H atoms]
|
|
|
|
Experimental value: 4.52 eV (0.166 Ha)
|
|
LDA typically overestimates by ~10-15%: ~4.9-5.1 eV
|
|
"""
|
|
|
|
from deepchem_dft import Molecule, Cell, run_scf
|
|
|
|
ECUT = 12.0 # Ha — fast demo; increase for convergence
|
|
VACUUM = 6.0 # Bohr vacuum on each side
|
|
TOL = 1e-5
|
|
|
|
EV = 27.2114 # Ha → eV conversion
|
|
|
|
# ── H atom ───────────────────────────────────────────────────────────────────
|
|
print("=" * 55)
|
|
print("H atom (1 electron, singly occupied orbital)")
|
|
print("=" * 55)
|
|
mol_H = Molecule.from_xyz("H 0.0 0.0 0.0", unit='angstrom')
|
|
cell_H, off_H = Cell.from_molecule(mol_H, vacuum=VACUUM, ecut=ECUT)
|
|
print(cell_H)
|
|
res_H = run_scf(cell_H, mol_H, offset=off_H, tol=TOL, verbose=True)
|
|
E_H = res_H.e_total
|
|
|
|
# ── H2 molecule ──────────────────────────────────────────────────────────────
|
|
print()
|
|
print("=" * 55)
|
|
print("H2 molecule (2 electrons, equilibrium bond 0.741 Å)")
|
|
print("=" * 55)
|
|
mol_H2 = Molecule.from_xyz("""
|
|
H 0.0 0.0 0.000
|
|
H 0.0 0.0 0.741
|
|
""", unit='angstrom')
|
|
cell_H2, off_H2 = Cell.from_molecule(mol_H2, vacuum=VACUUM, ecut=ECUT)
|
|
print(cell_H2)
|
|
res_H2 = run_scf(cell_H2, mol_H2, offset=off_H2, tol=TOL, verbose=True)
|
|
E_H2 = res_H2.e_total
|
|
|
|
# ── Atomization energy ────────────────────────────────────────────────────────
|
|
E_atom_Ha = E_H2 - 2 * E_H
|
|
E_atom_eV = E_atom_Ha * EV
|
|
|
|
print()
|
|
print("=" * 55)
|
|
print("Atomization energy (H2 → 2 H)")
|
|
print("=" * 55)
|
|
print(f" E(H) = {E_H:+.6f} Ha")
|
|
print(f" E(H2) = {E_H2:+.6f} Ha")
|
|
print(f" E_atom = E(H2) - 2E(H) = {E_atom_Ha:+.6f} Ha = {E_atom_eV:+.3f} eV")
|
|
print(f" Experiment: -0.1745 Ha = -4.748 eV (ZPE-corrected)")
|
|
print(f" LDA expected: ~-0.180 Ha = ~-4.9 eV (slight overbinding)")
|