Maple Quantum Chemistry Toolbox
The Maple Quantum Chemistry Toolbox from RDMChem, a separate add-on product to Maple, is a powerful environment for the computation and visualization of the electronic structure of molecules. In Maple 2020, this toolbox provides more capabilities for computing and visualizing excited states, new commands for saving and restoring your computations, additional lessons for classroom learning from introductory chemistry to advanced computational chemistry, updates to the interactive Maplet, and additional enhancements to many methods and commands throughout the package.
Note that the Maple Quantum Chemistry Toolbox is required in order to execute the examples in this worksheet.
Computing and Visualizing Excited States
Saving and Restoring Computations
Using the Package in the Classroom
More than ten new commands have been introduced for the computation and visualization of excited states. Here we examine five of these commands. Consider the molecule 5-bromouracil which can substitute for thymine in DNA. After loading the QuantumChemistry package,
withQuantumChemistry;
AOLabels,ActiveSpaceCI,ActiveSpaceSCF,AtomicData,BondAngles,BondDistances,Charges,ChargesPlot,CorrelationEnergy,CoupledCluster,DensityFunctional,DensityPlot3D,Dipole,DipolePlot,Energy,ExcitationEnergies,ExcitationSpectra,ExcitationSpectraPlot,ExcitedStateEnergies,ExcitedStateSpins,FullCI,GeometryOptimization,HartreeFock,Interactive,Isotopes,MOCoefficients,MODiagram,MOEnergies,MOIntegrals,MOOccupations,MOOccupationsPlot,MOSymmetries,MP2,MolecularData,MolecularGeometry,NuclearEnergy,NuclearGradient,OscillatorStrengths,Parametric2RDM,PlotMolecule,Populations,RDM1,RDM2,RTM1,ReadXYZ,Restore,Save,SaveXYZ,SearchBasisSets,SearchFunctionals,SkeletalStructure,Thermodynamics,TransitionDipolePlot,TransitionDipoles,TransitionOrbitalPlot,TransitionOrbitals,Variational2RDM,VibrationalModeAnimation,VibrationalModes,Video
we retrieve the molecular geometry from a web database with the command MolecularGeometry.
mol ≔ MolecularGeometry5-bromouracil;
mol≔Br,−2.97410000,0.79650000,0.00030000,O,−1.40900000,−1.87630000,−0.00010000,O,2.88650000,−0.19640000,0.00040000,N,0.74060000,−1.04550000,−0.00020000,N,1.10820000,1.26360000,−0.00010000,C,−0.64370000,−0.91570000,0.00010000,C,1.67030000,−0.00310000,−0.00010000,C,−1.13730000,0.49180000,−0.00010000,C,−0.24160000,1.48520000,−0.00010000,H,1.11300000,−1.99090000,−0.00010000,H,1.72120000,2.07350000,0,H,−0.52880000,2.53050000,−0.00010000
ExcitationSpectra
The command ExcitationSpectra computes the molecule's electronic excited-state spectra.
spectra ≔ ExcitationSpectramol;
We can set the optional keyword showtable to true to display a fancy table of the spectra.
spectra ≔ ExcitationSpectramol, showtable=true:
State
Energy
Wavelength
Spin
Oscillator
1
1.37044555eV
904.69991641nm
Triplet
0.09898469
2
3.15292785eV
393.23512351nm
0.17908969
3
4.02064106eV
308.36922639nm
0.00098620
4
4.98090933eV
248.91879991nm
Singlet
0.00034560
5
5.03831243eV
246.08278872nm
4.0284319310−6
6
5.93903201eV
208.76162508nm
2.7873863210−6
7
7.32989796eV
169.14859931nm
0.48096283
8
7.54582705eV
164.30829467nm
0.02027458
9
7.92221342eV
156.50196577nm
0.02748881
10
8.59118513eV
144.31559274nm
0.01227862
11
8.94012832eV
138.68279396nm
0.11361659
12
10.46325064eV
118.49491294nm
0.50944587
The table contains information such as the excitation energy, wavelength, spin, and oscillator strength for the first 12 states. Additional or fewer states may be requested with the keyword nstates. The default method is time-dependent Hartree-Fock, but other supported methods like time-dependent density functional theory can be selected by changing the base method.
spectra ≔ ExcitationSpectramol, method=DensityFunctional,nstates=1, showtable=true:
3.29944503eV
375.77288422nm
1.09430285
4.21812972eV
293.93168446nm
0.00036128
ExcitationSpectraPlot
The command ExcitationSpectraPlot generates a plot of the molecule's excitation spectra.
spectra ≔ ExcitationSpectraPlotmol;
We can use the keyword onlysinglets to show only the excitations from the ground state to singlet excited states.
TransitionDipolePlot
The command TransitionDipolePlot visualizes the transition dipole moment of a specified excited state.
TransitionDipolePlotmol, state=1;
TransitionOrbitals
The command TransitionOrbitals computes the most important orbitals for the ground-to-excited-state transition using the singular value decomposition.
sv,orbital_ground,orbital_excited ≔ TransitionOrbitalsmol, state=1:
The singular values are returned as a vector.
sv;
We observe that the first three singular values are "large" with the remaining values being negligible. The nth singular value corresponds to a transition from the molecular orbital in the corresponding column of orbital_ground to the molecular orbital in the corresponding column of orbital_excited.
TransitionOrbitalPlot
The command TransitionOrbitalPlot visualizes the pair of orbitals involved in a ground-to-excited-state transition for the nth singular value. For example, for the transition from the ground state to the first excited state, we can visualize the most important pair of orbitals associated with the first singular value.
TransitionOrbitalPlotmol, state = 1, orbitalindex=1;
Note that the change in electron density from the left (ground-state occupied orbital) to the right (excited-state occupied orbital) is consistent with the direction of the transition dipole moment (computed above).
Similarly, we can visualize the second most important pair of orbitals associated with the second singular value.
TransitionOrbitalPlotmol, state = 1, orbitalindex=2;
Again we observe that the change in electron density agrees with the direction of the transition dipole moment.
Optional Keywords
The excited-state data can also be computed directly in a supporting method. For example, we can perform a DensityFunctional computation that automatically computed the excited-state data by setting the optional keyword excited_states to true, "TDDFT" (Time-dependent Density Functional Theory) , or "TDA" (Tamm-Dancoff Approximation). Furthermore, we can control the number of computed singlet and triplet states with the optional keyword nstates whose default is 6 for 6 singlet states and 6 triplet states.
data_dft ≔ DensityFunctionalmol, excited_states=TDA,nstates=1;
The table contains entries for the excited_state_energies, excited_state_spins, transition_dipoles, oscillator_strengths, and the 1-electron reduced transition matrices (1-RTMs), for example:
data_dftexcited_state_energies;
−2954.22143496−2954.19773755
data_dftexcited_state_spins;
TripletSinglet
data_dftoscillator_strengths;
0.772668880.00048702
Two new commands Save and Restore have been added to the package to make it easier than ever to save your data for future analysis. The Save command saves the cache tables associated with your calculations while the Restore command restores these cache tables in a new Maple session. This pair of commands makes it easy to switch between Maple interfaces (such as between the Maple worksheet and the Maple command-line version or between the Maple worksheet and the interactive Maplet), between operating systems (Windows, Mac OS X, or Linux), or even between different computers without redoing your computations. For advanced computations it is convenient to run computations with the Maple command-line version and then perform your analysis and visualization in the Maple worksheet (or the interactive Maplet). In the Maple worksheet or Maple document mode the commands Save() and Restore() open a directory/file browser to facilitate file selection.
The Maple Quantum Chemistry Toolbox includes many lessons that can be used in chemistry and physics courses from advanced high school courses through the graduate level. These lessons and associated curricula provide instructors and students with real-time quantum chemistry computations and visualizations that quickly deepen understanding of molecular concepts. With Maple 2020, the support for Introductory (General) Chemistry, Computational Chemistry, and Quantum Chemistry has been extended even further, with the addition of over 13 new topics, including atomic structure, chemical bonding, valence-shell electron pair repulsion (VSEPR), the Maxwell-Boltzmann distribution, heat capacity, enthalpy, entropy, free energy, as well as advanced electronic structure methods. The general chemistry lessons are ideal not only for college-level introductory chemistry courses but also secondary-school Advanced Placement (AP) and International Baccalaureate (IB) courses. Computational chemistry lessons cover electronic structure theory and methods including Hartree-Fock method, Moller-Plessett perturbation theory, density functional theory, and 2-RDM theory for use in undergraduate and graduate computational chemistry courses. In the following two subsections we show snippets of the lessons on the Maxwell-Boltzmann distribution for Introductory Chemistry and density functional theory for Computational Chemistry.
Maxwell-Boltzmann Distribution
Before we begin, we define the ideal gas constant in SI units.
RSI ≔ evalfScientificConstantsConstant'R';
RSI≔8.3144598000
We can define the Maxwell-Boltzmann distribution with Maple.
fv,M,T,R ≔ 4⋅Pi⋅M2⋅Pi⋅R⋅T32v2exp−M⋅v22⋅R⋅T;
f≔v,M,T,R↦π⋅2⋅Mπ⋅R⋅T32⋅v2⋅ⅇ−M⋅v22⋅R⋅T
Let us check that the total probability of the particle having some speed is unity for a gas with a molar mass of 0.1 kg/mol at a temperature of 273 K.
P_total ≔ intfv,0.1,273,RSI, v=0..infinity;
P_total≔0.9999999999
(a) Does the total probability equal the expected value of unity?
We can explore the Maxwell-Boltzmann distribution as a function of the temperature and molar mass of the gas.
Exploreplotfv,M,T,RSI, v=0..2000, axes=boxed,labels='vms','fv', color=blue, thickness=2, T=100..1000,M=0.01..1;
T
M
Using the slide rules to change temperature T in K and molar mass M in kg/mol (SI units), answer the following questions:
(b) Approximately from the graph, what is the most probably velocity for water H2O at 275 K?
(c) Approximately from the graph, what is the most probably velocity for methane CH4 at 550 K?
Density Functional Theory
We compute the energy of diphenyl disulphide. First, we define the geometry, which we obtained from the MolecularGeometry command.
mol ≔ S,−0.55610000,−1.83210000,−0.97360000,S,0.62630000,−1.92950000,0.70830000,C,−1.79690000,−0.65310000,−0.47620000,C,1.82150000,−0.64580000,0.38420000,C,−2.60700000,−0.91680000,0.62830000,C,1.57350000,0.65950000,0.80660000,C,−1.95740000,0.53360000,−1.19130000,C,3.00290000,−0.94940000,−0.29110000,C,−3.57770000,0.00620000,1.01760000,C,2.50850000,1.66310000,0.55330000,C,−2.92800000,1.45670000,−0.80200000,C,3.93800000,0.05420000,−0.54430000,C,−3.73820000,1.19290000,0.30240000,C,3.69070000,1.36050000,−0.12220000,H,−2.50880000,−1.83880000,1.19590000,H,0.65860000,0.91480000,1.33560000,H,−1.33440000,0.75710000,−2.05390000,H,3.21330000,−1.96160000,−0.62770000,H,−4.21010000,−0.19990000,1.87620000,H,2.31610000,2.68020000,0.88210000,H,−3.05270000,2.38110000,−1.35840000,H,4.85860000,−0.18120000,−1.07020000,H,−4.49420000,1.91150000,0.60530000,H,4.41880000,2.14190000,−0.31940000;
mol≔S,−0.55610000,−1.83210000,−0.97360000,S,0.62630000,−1.92950000,0.70830000,C,−1.79690000,−0.65310000,−0.47620000,C,1.82150000,−0.64580000,0.38420000,C,−2.60700000,−0.91680000,0.62830000,C,1.57350000,0.65950000,0.80660000,C,−1.95740000,0.53360000,−1.19130000,C,3.00290000,−0.94940000,−0.29110000,C,−3.57770000,0.00620000,1.01760000,C,2.50850000,1.66310000,0.55330000,C,−2.92800000,1.45670000,−0.80200000,C,3.93800000,0.05420000,−0.54430000,C,−3.73820000,1.19290000,0.30240000,C,3.69070000,1.36050000,−0.12220000,H,−2.50880000,−1.83880000,1.19590000,H,0.65860000,0.91480000,1.33560000,H,−1.33440000,0.75710000,−2.05390000,H,3.21330000,−1.96160000,−0.62770000,H,−4.21010000,−0.19990000,1.87620000,H,2.31610000,2.68020000,0.88210000,H,−3.05270000,2.38110000,−1.35840000,H,4.85860000,−0.18120000,−1.07020000,H,−4.49420000,1.91150000,0.60530000,H,4.41880000,2.14190000,−0.31940000
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.
data ≔ DensityFunctionalmol, basis=6-31g, xc =B3LYP;
The ground-state energy is given by
datae_tot;
−1259.07521547
(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.
datamo_occ;
DensityPlot3Dmol,data,basis=6-31g, orbitalindex=57; # HOMO
DensityPlot3Dmol,data,basis=6-31g, orbitalindex=58; # LUMO
The energies of these two orbitals are
datamo_energy57..58;
−0.22884080−0.04809863
(b) Compare these orbitals with those from the PBE functional.
Download Help Document
