HF in VeloxChem#
Define the molecule#
We first define the structure of a water molecule and choose a basis set
import veloxchem as vlx
# Atomic coordinates in units of Angstrom
mol_xyz = """3
O 0.0000000000 0.1178336003 0.0000000000
H -0.7595754146 -0.4713344012 -0.0000000000
H 0.7595754146 -0.4713344012 0.0000000000
"""
molecule = vlx.Molecule.read_xyz_string(mol_xyz)
basis = vlx.MolecularBasis.read(molecule, "cc-pVDZ")
You appear to be running in JupyterLab (or JavaScript failed to load for some other reason). You need to install the 3dmol extension:
jupyter labextension install jupyterlab_3dmol
print('Number of atoms:', molecule.number_of_atoms())
print('Number of electrons:', molecule.number_of_electrons())
print('Number of contracted basis functions:', basis.get_dimensions_of_basis(molecule))
Number of atoms: 3
Number of electrons: 10
Number of contracted basis functions: 24
SCF optimization#
Perform a self-consistent field (SCF) optimization to obtain the Hartree–Fock wave function and the associated ground-state energy.
scf_drv = vlx.ScfRestrictedDriver()
scf_results = scf_drv.compute(molecule, basis)
Self Consistent Field Driver Setup
====================================
Wave Function Model : Spin-Restricted Hartree-Fock
Initial Guess Model : Superposition of Atomic Densities
Convergence Accelerator : Two Level Direct Inversion of Iterative Subspace
Max. Number of Iterations : 50
Max. Number of Error Vectors : 10
Convergence Threshold : 1.0e-06
ERI Screening Scheme : Cauchy Schwarz + Density
ERI Screening Mode : Dynamic
ERI Screening Threshold : 1.0e-12
Linear Dependence Threshold : 1.0e-06
* Info * Nuclear repulsion energy: 9.1561447194 a.u.
* Info * Overlap matrix computed in 0.00 sec.
* Info * Kinetic energy matrix computed in 0.00 sec.
* Info * Nuclear potential matrix computed in 0.00 sec.
* Info * Orthogonalization matrix computed in 0.01 sec.
* Info * SAD initial guess computed in 0.00 sec.
* Info * Starting Reduced Basis SCF calculation...
* Info * ...done. SCF energy in reduced basis set: -75.980526010891 a.u. Time: 0.08 sec.
* Info * Overlap matrix computed in 0.00 sec.
* Info * Kinetic energy matrix computed in 0.00 sec.
* Info * Nuclear potential matrix computed in 0.00 sec.
* Info * Orthogonalization matrix computed in 0.00 sec.
Iter. | Hartree-Fock Energy | Energy Change | Gradient Norm | Max. Gradient | Density Change
--------------------------------------------------------------------------------------------
1 -76.025514170400 0.0000000000 0.09661894 0.01034981 0.00000000
2 -76.026530063455 -0.0010158931 0.01856585 0.00282747 0.02538073
3 -76.026575250863 -0.0000451874 0.00393964 0.00064892 0.00485806
4 -76.026577667709 -0.0000024168 0.00240246 0.00040859 0.00161571
5 -76.026578208026 -0.0000005403 0.00021790 0.00003292 0.00050560
6 -76.026578219600 -0.0000000116 0.00003172 0.00000533 0.00012152
7 -76.026578219828 -0.0000000002 0.00000210 0.00000030 0.00001639
8 -76.026578219829 -0.0000000000 0.00000055 0.00000009 0.00000089
*** SCF converged in 8 iterations. Time: 0.60 sec.
Spin-Restricted Hartree-Fock:
-----------------------------
Total Energy : -76.0265782198 a.u.
Electronic Energy : -85.1827229393 a.u.
Nuclear Repulsion Energy : 9.1561447194 a.u.
------------------------------------
Gradient Norm : 0.0000005463 a.u.
Ground State Information
------------------------
Charge of Molecule : 0.0
Multiplicity (2S+1) : 1.0
Magnetic Quantum Number (M_S) : 0.0
Spin Restricted Orbitals
------------------------
Molecular Orbital No. 1:
--------------------------
Occupation: 2.000 Energy: -20.55120 a.u.
( 1 O 1s : 1.00)
Molecular Orbital No. 2:
--------------------------
Occupation: 2.000 Energy: -1.33473 a.u.
( 1 O 2s : -0.44) ( 1 O 3s : -0.38) ( 2 H 1s : -0.19)
( 3 H 1s : -0.19)
Molecular Orbital No. 3:
--------------------------
Occupation: 2.000 Energy: -0.69693 a.u.
( 1 O 1p+1: -0.49) ( 1 O 2p+1: -0.22) ( 2 H 1s : 0.33)
( 3 H 1s : -0.33)
Molecular Orbital No. 4:
--------------------------
Occupation: 2.000 Energy: -0.56605 a.u.
( 1 O 2s : 0.15) ( 1 O 3s : 0.35) ( 1 O 1p-1: 0.54)
( 1 O 2p-1: 0.37) ( 2 H 1s : -0.21) ( 3 H 1s : -0.21)
Molecular Orbital No. 5:
--------------------------
Occupation: 2.000 Energy: -0.49290 a.u.
( 1 O 1p0 : 0.63) ( 1 O 2p0 : 0.50)
Molecular Orbital No. 6:
--------------------------
Occupation: 0.000 Energy: 0.18486 a.u.
( 1 O 3s : 1.00) ( 1 O 1p-1: -0.19) ( 1 O 2p-1: -0.34)
( 2 H 2s : -0.83) ( 3 H 2s : -0.83)
Molecular Orbital No. 7:
--------------------------
Occupation: 0.000 Energy: 0.25567 a.u.
( 1 O 1p+1: -0.28) ( 1 O 2p+1: -0.67) ( 2 H 2s : -1.44)
( 3 H 2s : 1.44)
Molecular Orbital No. 8:
--------------------------
Occupation: 0.000 Energy: 0.78605 a.u.
( 1 O 1p+1: 0.27) ( 1 O 2p+1: 0.45) ( 2 H 1s : 0.94)
( 2 H 2s : -0.69) ( 2 H 1p-1: -0.15) ( 3 H 1s : -0.94)
( 3 H 2s : 0.69) ( 3 H 1p-1: 0.15)
Molecular Orbital No. 9:
--------------------------
Occupation: 0.000 Energy: 0.85103 a.u.
( 1 O 2s : 0.26) ( 1 O 3s : -0.32) ( 1 O 1p-1: -0.33)
( 2 H 1s : -0.79) ( 2 H 2s : 0.55) ( 2 H 1p+1: 0.29)
( 3 H 1s : -0.79) ( 3 H 2s : 0.55) ( 3 H 1p+1: -0.29)
Molecular Orbital No. 10:
--------------------------
Occupation: 0.000 Energy: 1.16387 a.u.
( 1 O 3s : 0.74) ( 1 O 1p-1: 0.76) ( 1 O 2p-1: -1.27)
( 2 H 1s : -0.54) ( 2 H 1p-1: -0.25) ( 3 H 1s : -0.54)
( 3 H 1p-1: -0.25)
SCF information#
The SCF driver object has a method named get_scf_energy() for retrieving the final energy.
print(f'Hartree–Fock energy: {scf_drv.get_scf_energy():14.10f} a.u.')
Hartree–Fock energy: -76.0265782198 a.u.
The return object from the compute() method is a Python dictionary containing several tensors:
C: molecular orbital coefficients as a NumPy arrayE: orbital energies as a NumPy arrayD: \(\alpha\)- and \(\beta\)-spin density matrices as a tuple of NumPy arraysF: \(\alpha\)- and \(\beta\)-spin Fock matrices as a tuple of NumPy arraysS: overlap integrals as a NumPy array
print('Dictionary keys:\n', scf_results.keys())
print()
print('Orbital energies:\n', scf_results['E_alpha'])
Dictionary keys:
dict_keys(['eri_thresh', 'qq_type', 'scf_type', 'scf_energy', 'restart', 'S', 'C_alpha', 'C_beta', 'E_alpha', 'E_beta', 'D_alpha', 'D_beta', 'F_alpha', 'F_beta', 'C', 'E', 'D', 'F'])
Orbital energies:
[-20.55119961 -1.33473347 -0.69692651 -0.56605275 -0.49289827
0.18486487 0.25566978 0.78604709 0.85103488 1.16387193
1.20030818 1.25366153 1.44422793 1.4755769 1.67329727
1.8679223 1.93111284 2.4430781 2.48048696 3.28268216
3.334531 3.50569214 3.86065795 4.14355137]
Visualizing molecular orbitals#
The resulting molecular orbitals (MOs) can be visualized using OrbitalViewer
viewer = vlx.OrbitalViewer()
viewer.plot(molecule, basis, scf_drv.mol_orbs)