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 now 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 the second most important pair of orbitals associated with the second singular value
TransitionOrbitalPlotmol, state = 1, orbitalindex=2;
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), i.e.
data_dftexcited_state_energies;
−2954.22143496−2954.19773755
data_dftexcited_state_spins;
TripletSinglet
data_dftoscillator_strengths;
0.772668880.00048702
ActiveSpaceCI and ActiveSpaceSCF
Excited states are now easily computed with the ActiveSpaceCI and ActiveSpaceSCF commands with the new keyword state. For example, we can readily compute the first excited-state energy and properties of the molecule 5-bromouracil
data_casci ≔ ActiveSpaceCImol,active=6,6, state=1;
The new entry spin_squared in the output table gives us the expectation value of the total electron spin squared, equal to S(S+1). Here we observe that this expectation value is 2, which indicates that the first excited state is a triplet (S =1 or 2S+1=3).
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 (i.e. between the Maple worksheet and the Maple command line or between the Maple worksheet and the interactive Maplet), between operating systems (i.e. Windows, MacOS, 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 and then perform your analysis and visualization in the Maple worksheet (or the interactive Maplet). On 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
