Optimizing electronic structure#
Here we consider the optimization of the electronic structure using Hartree–Fock (HF) and Kohn-Sham density functional theory (KS-DFT), and how to plot the resulting molecular orbitals. Our example is methanol, using the xTB structure from the previous section.
import py3Dmol as p3d
import veloxchem as vlx
Specify structure#
The molecular geometry is written in xyz-coordinates, and visualized using py3Dmol.
methanol_xyz = """6
C -0.362051846542 0.008519935441 0.024906559609
O 0.933468973663 0.509856441784 -0.193822841148
H -0.483698601596 -0.368885697455 1.046180408743
H -0.612499226929 -0.790059463155 -0.682044832355
H -1.045997965440 0.841433017766 -0.127103979461
H 1.570778666846 -0.200864234382 -0.068115315388
"""
viewer = p3d.view(width=300, height=200)
viewer.addModel(methanol_xyz, 'xyz')
viewer.setViewStyle({"style": "outline", "width": 0.05})
viewer.setStyle({"stick":{},"sphere": {"scale":0.25}})
viewer.show()
3Dmol.js failed to load for some reason. Please check your browser console for error messages.
SCF optimization#
At the HF/6-311G level of theory, the electronic structure is optimized by creating the molecule
and basis
objects, setting up an unrestricted SCF driver, and performing the SCF calculation:
molecule = vlx.Molecule.read_xyz_string(methanol_xyz)
basis = vlx.MolecularBasis.read(molecule, "6-311G")
scf_hf = vlx.ScfRestrictedDriver()
hf_results = scf_hf.compute(molecule, basis)
print(f'The total (HF) energy is {scf_hf.get_scf_energy(): .5f} Hartree')
The total (HF) energy is -115.01894 Hartree
Visualizing molecular orbitals#
The MOs can be visualized using the OrbitalViewer
function:
viewer = vlx.OrbitalViewer()
viewer.plot(molecule, basis, scf_hf.mol_orbs)
Here you can only see the default MO selection (the HOMO), but if you run this in a notebook you can visualize any MO, change isovalues, and more.
KS-DFT optimization#
With the B3LYP exchange-correlation functional, the SCF is optimized using:
molecule = vlx.Molecule.read_xyz_string(methanol_xyz)
basis = vlx.MolecularBasis.read(molecule, "6-311G")
scf_b3lyp = vlx.ScfRestrictedDriver()
scf_b3lyp.xcfun = 'b3lyp'
b3lyp_results = scf_b3lyp.compute(molecule, basis)
print(f'The total (B3LYP) energy is {scf_b3lyp.get_scf_energy(): .5f} Hartree')
The total (B3LYP) energy is -115.71437 Hartree
And the MOs are visualized as:
viewer = vlx.OrbitalViewer()
viewer.plot(molecule, basis, scf_b3lyp.mol_orbs)