DFT - Maple Help

Density Functional Theory

 Overview In 1964 Hohenberg and Kohn proved that the ground-state energy of any atom or molecule can be expressed as a functional of its 1-electron density, known as the Hohenberg-Kohn theorem.  The theorem gave birth to density functional theory (DFT) in which the 1-electron density rather than the many-electron wave function becomes the basic variable of the quantum chemistry. The 1-electron density is the probability of finding an electron at a given spatial position of the molecule.  Like Hartree-Fock theory DFT can be solved through an effective 1-electron Schrödinger equation, but unlike Hartree-Fock theory DFT can in principle recover the correlation energy, the energy difference between the energy from full configuration interaction and the energy from the Hartree-Fock method.  The exact energy functional in DFT is not known, and hence, practically, DFT is implemented with approximate functionals. Figure 1: 1-electron density of diphenyl disulfide

DFT

After loading the Quantum Chemistry package, we explore the ground-state energy and orbitals, dipole moment, and ionization energy of diphenyl disulphide with DFT.

Quantum Chemistry

We set the number of Digits to be used in computations to 15 and load the Quantum Chemistry package using Maple's with command.

 >
 ${\mathrm{Digits}}{≔}{15}$ (2.1.1)
 > $\mathrm{with}\left(\mathrm{QuantumChemistry}\right);$
 $\left[{\mathrm{AOLabels}}{,}{\mathrm{ActiveSpaceCI}}{,}{\mathrm{ActiveSpaceSCF}}{,}{\mathrm{AtomicData}}{,}{\mathrm{BondAngles}}{,}{\mathrm{BondDistances}}{,}{\mathrm{Charges}}{,}{\mathrm{ChargesPlot}}{,}{\mathrm{CorrelationEnergy}}{,}{\mathrm{CoupledCluster}}{,}{\mathrm{DensityFunctional}}{,}{\mathrm{DensityPlot3D}}{,}{\mathrm{Dipole}}{,}{\mathrm{DipolePlot}}{,}{\mathrm{Energy}}{,}{\mathrm{ExcitationEnergies}}{,}{\mathrm{ExcitationSpectra}}{,}{\mathrm{ExcitationSpectraPlot}}{,}{\mathrm{ExcitedStateEnergies}}{,}{\mathrm{ExcitedStateSpins}}{,}{\mathrm{FullCI}}{,}{\mathrm{GeometryOptimization}}{,}{\mathrm{HartreeFock}}{,}{\mathrm{Interactive}}{,}{\mathrm{Isotopes}}{,}{\mathrm{MOCoefficients}}{,}{\mathrm{MODiagram}}{,}{\mathrm{MOEnergies}}{,}{\mathrm{MOIntegrals}}{,}{\mathrm{MOOccupations}}{,}{\mathrm{MOOccupationsPlot}}{,}{\mathrm{MOSymmetries}}{,}{\mathrm{MP2}}{,}{\mathrm{MolecularData}}{,}{\mathrm{MolecularGeometry}}{,}{\mathrm{NuclearEnergy}}{,}{\mathrm{NuclearGradient}}{,}{\mathrm{OscillatorStrengths}}{,}{\mathrm{Parametric2RDM}}{,}{\mathrm{PlotMolecule}}{,}{\mathrm{Populations}}{,}{\mathrm{RDM1}}{,}{\mathrm{RDM2}}{,}{\mathrm{RTM1}}{,}{\mathrm{ReadXYZ}}{,}{\mathrm{Restore}}{,}{\mathrm{Save}}{,}{\mathrm{SaveXYZ}}{,}{\mathrm{SearchBasisSets}}{,}{\mathrm{SearchFunctionals}}{,}{\mathrm{SkeletalStructure}}{,}{\mathrm{Thermodynamics}}{,}{\mathrm{TransitionDipolePlot}}{,}{\mathrm{TransitionDipoles}}{,}{\mathrm{TransitionOrbitalPlot}}{,}{\mathrm{TransitionOrbitals}}{,}{\mathrm{Variational2RDM}}{,}{\mathrm{VibrationalModeAnimation}}{,}{\mathrm{VibrationalModes}}{,}{\mathrm{Video}}\right]$ (2.1.2)
 >

Energy

We compute the energy of diphenyl disulphide.

First, we define the geometry, which we obtained from the MolecularGeometry command

 >
 ${\mathrm{mol}}{≔}\left[\left[{"S"}{,}{-0.55610000}{,}{-1.83210000}{,}{-0.97360000}\right]{,}\left[{"S"}{,}{0.62630000}{,}{-1.92950000}{,}{0.70830000}\right]{,}\left[{"C"}{,}{-1.79690000}{,}{-0.65310000}{,}{-0.47620000}\right]{,}\left[{"C"}{,}{1.82150000}{,}{-0.64580000}{,}{0.38420000}\right]{,}\left[{"C"}{,}{-2.60700000}{,}{-0.91680000}{,}{0.62830000}\right]{,}\left[{"C"}{,}{1.57350000}{,}{0.65950000}{,}{0.80660000}\right]{,}\left[{"C"}{,}{-1.95740000}{,}{0.53360000}{,}{-1.19130000}\right]{,}\left[{"C"}{,}{3.00290000}{,}{-0.94940000}{,}{-0.29110000}\right]{,}\left[{"C"}{,}{-3.57770000}{,}{0.00620000}{,}{1.01760000}\right]{,}\left[{"C"}{,}{2.50850000}{,}{1.66310000}{,}{0.55330000}\right]{,}\left[{"C"}{,}{-2.92800000}{,}{1.45670000}{,}{-0.80200000}\right]{,}\left[{"C"}{,}{3.93800000}{,}{0.05420000}{,}{-0.54430000}\right]{,}\left[{"C"}{,}{-3.73820000}{,}{1.19290000}{,}{0.30240000}\right]{,}\left[{"C"}{,}{3.69070000}{,}{1.36050000}{,}{-0.12220000}\right]{,}\left[{"H"}{,}{-2.50880000}{,}{-1.83880000}{,}{1.19590000}\right]{,}\left[{"H"}{,}{0.65860000}{,}{0.91480000}{,}{1.33560000}\right]{,}\left[{"H"}{,}{-1.33440000}{,}{0.75710000}{,}{-2.05390000}\right]{,}\left[{"H"}{,}{3.21330000}{,}{-1.96160000}{,}{-0.62770000}\right]{,}\left[{"H"}{,}{-4.21010000}{,}{-0.19990000}{,}{1.87620000}\right]{,}\left[{"H"}{,}{2.31610000}{,}{2.68020000}{,}{0.88210000}\right]{,}\left[{"H"}{,}{-3.05270000}{,}{2.38110000}{,}{-1.35840000}\right]{,}\left[{"H"}{,}{4.85860000}{,}{-0.18120000}{,}{-1.07020000}\right]{,}\left[{"H"}{,}{-4.49420000}{,}{1.91150000}{,}{0.60530000}\right]{,}\left[{"H"}{,}{4.41880000}{,}{2.14190000}{,}{-0.31940000}\right]\right]$ (2.2.1)

Then we compute the energy from DFT with the DensityFunctional command.  We select the B3LYP exchange-correlation (xc) functional and the 6-31g basis set

 >
 ${\mathrm{table}}{}\left({\mathrm{%id}}{=}{18446744437503033502}\right)$ (2.2.2)
 > 

The ground-state energy is given by

 > $\mathrm{data}\left[\mathrm{e_tot}\right];$
 ${-1259.07521547}$ (2.2.3)

(a) Compute the energy using the PBE functional, and compare your result with the energy from the B3LYP functional.  Absolute energies are typically not too important in DFT.

The highest occupied molecular orbital (HOMO) can be visualized with the DensityPlot3D command

 > $\mathrm{data}{\left[\mathrm{mo_occ}\right]}^{};$
 ${{\mathrm{_rtable}}}_{{18446744437335435974}}$ (2.2.4)
 >
 >

The energies of these two orbitals are

 > $\mathrm{data}\left[\mathrm{mo_energy}\right]\left[57..58\right];$
 $\left[\begin{array}{c}-0.228840798\\ -0.0480986315\end{array}\right]$ (2.2.5)
 > 

(b) Compare these orbitals with those from the PBE functional.

Dipole Moment

We have the dipole moment from our previous calculation

 > $\mathrm{data}\left[\mathrm{dipole}\right];$
 $\left[\begin{array}{c}0.114940823\\ 2.90462002\\ 0.194775483\end{array}\right]$ (2.3.1)

We can visualize the direction of the dipole moment with the DipolePlot command (click on molecule to rotate the perspective)

 >

(c) Compute the dipole and its plot with the PBE functional, and compare the results with those from B3LYP.  Does the direction of the dipole make sense?

Ionization Energy

We can compute the ionization energy by computing the energy of the ionized molecule

 >
 ${\mathrm{table}}{}\left({\mathrm{%id}}{=}{18446744437330283102}\right)$ (2.4.1)

Therefore, the ionization energy is given by

 >
 ${\mathrm{IE_hartree}}{≔}{0.29260440}{}⟦{\mathrm{E0}}⟧$ (2.4.2)

which we can convert from hartrees to eV

 >
 ${\mathrm{IE_eV}}{≔}{7.96217130}{}⟦{\mathrm{eV}}⟧$ (2.4.3)
 > 

(d) Compare this ionization energy (in eV) with the ionization energy (in eV) predicted using the molecular orbital energies by Koopman's theorem.

References

 1 R. G. Parr and W. Yang, Density Functional Theory of Atoms and Molecules (Oxford University Press, New York, 1989).
 2 P. Hohenberg and W. Kohn, Phys. Rev. 136, B864 (1964). "Inhomogeneous electron gas"
 3 R. O. Jones, Rev. Mod. Phys. 87, 897 (2015). "Density functional theory: Its origins, rise to prominence, and future"