Quantum Chemistry - New Features in Maple 2020 - Maplesoft

What's New in Maple 2020

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 

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, 

> with(QuantumChemistry); 1

[AOLabels, ActiveSpaceCI, ActiveSpaceSCF, AtomicData, BondAngles, BondDistances, Charges, ChargesPlot, CorrelationEnergy, CoupledCluster, DensityFunctional, DensityPlot3D, Dipole, DipolePlot, Energy, ...
[AOLabels, ActiveSpaceCI, ActiveSpaceSCF, AtomicData, BondAngles, BondDistances, Charges, ChargesPlot, CorrelationEnergy, CoupledCluster, DensityFunctional, DensityPlot3D, Dipole, DipolePlot, Energy, ...
[AOLabels, ActiveSpaceCI, ActiveSpaceSCF, AtomicData, BondAngles, BondDistances, Charges, ChargesPlot, CorrelationEnergy, CoupledCluster, DensityFunctional, DensityPlot3D, Dipole, DipolePlot, Energy, ...
[AOLabels, ActiveSpaceCI, ActiveSpaceSCF, AtomicData, BondAngles, BondDistances, Charges, ChargesPlot, CorrelationEnergy, CoupledCluster, DensityFunctional, DensityPlot3D, Dipole, DipolePlot, Energy, ...
[AOLabels, ActiveSpaceCI, ActiveSpaceSCF, AtomicData, BondAngles, BondDistances, Charges, ChargesPlot, CorrelationEnergy, CoupledCluster, DensityFunctional, DensityPlot3D, Dipole, DipolePlot, Energy, ...
[AOLabels, ActiveSpaceCI, ActiveSpaceSCF, AtomicData, BondAngles, BondDistances, Charges, ChargesPlot, CorrelationEnergy, CoupledCluster, DensityFunctional, DensityPlot3D, Dipole, DipolePlot, Energy, ...
[AOLabels, ActiveSpaceCI, ActiveSpaceSCF, AtomicData, BondAngles, BondDistances, Charges, ChargesPlot, CorrelationEnergy, CoupledCluster, DensityFunctional, DensityPlot3D, Dipole, DipolePlot, Energy, ...
 

we retrieve the molecular geometry from a web database with the command MolecularGeometry.

> mol := MolecularGeometry(
 

Typesetting:-mprintslash([mol := [[
Typesetting:-mprintslash([mol := [[
Typesetting:-mprintslash([mol := [[
Typesetting:-mprintslash([mol := [[
Typesetting:-mprintslash([mol := [[
 

ExcitationSpectra 

The command ExcitationSpectra computes the molecule's electronic excited-state spectra. 

> spectra := ExcitationSpectra(mol); 1
 

Typesetting:-msub(Typesetting:-mi(
 

We can set the optional keyword showtable to true to display a fancy table of the spectra. 

> spectra := ExcitationSpectra(mol, showtable = true); -1
 

State

Energy Wavelength Spin Oscillator
1; `+`(`*`(1.37044554916321, `*`(Unit('eV')))); `+`(`*`(904.699916405837, `*`(Unit('nm')))); Triplet 0.9898469198503e-1;
2; `+`(`*`(3.15292785322713, `*`(Unit('eV')))); `+`(`*`(393.235123505186, `*`(Unit('nm')))); Triplet .17908969254738;
3; `+`(`*`(4.02064106156715, `*`(Unit('eV')))); `+`(`*`(308.36922639482, `*`(Unit('nm')))); Triplet 0.986203922715837e-3;
4; `+`(`*`(4.98090933345577, `*`(Unit('eV')))); `+`(`*`(248.918799914494, `*`(Unit('nm')))); Singlet 0.345599558886734e-3;
5; `+`(`*`(5.03831243217264, `*`(Unit('eV')))); `+`(`*`(246.082788722981, `*`(Unit('nm')))); Triplet `+`(`*`(4.02843192978929700, `*`(`^`(10, -6))));
6; `+`(`*`(5.9390320096651, `*`(Unit('eV')))); `+`(`*`(208.761625084527, `*`(Unit('nm')))); Singlet `+`(`*`(2.78738632432828597, `*`(`^`(10, -6))));
7; `+`(`*`(7.32989796449744, `*`(Unit('eV')))); `+`(`*`(169.148599308192, `*`(Unit('nm')))); Singlet .480962825722659;
8; `+`(`*`(7.54582704589124, `*`(Unit('eV')))); `+`(`*`(164.308294667556, `*`(Unit('nm')))); Triplet 0.202745802856289e-1;
9; `+`(`*`(7.92221342183102, `*`(Unit('eV')))); `+`(`*`(156.501965770086, `*`(Unit('nm')))); Triplet 0.274888128221947e-1;
10; `+`(`*`(8.5911851257368, `*`(Unit('eV')))); `+`(`*`(144.315592740807, `*`(Unit('nm')))); Singlet 0.122786201668375e-1;
11; `+`(`*`(8.94012832023308, `*`(Unit('eV')))); `+`(`*`(138.682793955063, `*`(Unit('nm')))); Singlet .11361659434923;
12; `+`(`*`(10.4632506405909, `*`(Unit('eV')))); `+`(`*`(118.49491294387, `*`(Unit('nm')))); Singlet .50944587230632;

 

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 := ExcitationSpectra(mol, method = DensityFunctional, nstates = 1, showtable = true); -1
 

State

Energy Wavelength Spin Oscillator
1; `+`(`*`(3.29944502607623, `*`(Unit('eV')))); `+`(`*`(375.772884217792, `*`(Unit('nm')))); Triplet 1.09430285120644;
2; `+`(`*`(4.21812971960686, `*`(Unit('eV')))); `+`(`*`(293.931684462816, `*`(Unit('nm')))); Singlet 0.361277405328289e-3;


ExcitationSpectraPlot 

The command ExcitationSpectraPlot generates a plot of the molecule's excitation spectra.

> spectra := ExcitationSpectraPlot(mol); 1
 

Plot_2d
 

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. 

> TransitionDipolePlot(mol, state = 1); 1
 

Plot_2d
 

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 := TransitionOrbitals(mol, state = 1); -1
 

The singular values are returned as a vector. 

> sv; 1
 

Typesetting:-msub(Typesetting:-mi(
 

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. 

> TransitionOrbitalPlot(mol, state = 1, orbitalindex = 1); 1
 

Plot_2d Plot_2d

 

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. 

> TransitionOrbitalPlot(mol, state = 1, orbitalindex = 2); 1
 

Plot_2d Plot_2d

 

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 := DensityFunctional(mol, excited_states =
 

Typesetting:-mrow(Typesetting:-mi(
Typesetting:-mrow(Typesetting:-mi(
Typesetting:-mrow(Typesetting:-mi(
Typesetting:-mrow(Typesetting:-mi(
Typesetting:-mrow(Typesetting:-mi(
Typesetting:-mrow(Typesetting:-mi(
Typesetting:-mrow(Typesetting:-mi(
Typesetting:-mrow(Typesetting:-mi(
Typesetting:-mrow(Typesetting:-mi(
 

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_dft[excited_state_energies]; 1

Typesetting:-mfenced(Typesetting:-mrow(Typesetting:-mtable(Typesetting:-mtr(Typesetting:-mtd(Typesetting:-mrow(Typesetting:-mo(
 

> data_dft[excited_state_spins]; 1

Typesetting:-mfenced(Typesetting:-mrow(Typesetting:-mtable(Typesetting:-mtr(Typesetting:-mtd(Typesetting:-ms(
 

> data_dft[oscillator_strengths]; 1

Typesetting:-mfenced(Typesetting:-mrow(Typesetting:-mtable(Typesetting:-mtr(Typesetting:-mtd(Typesetting:-mn(
 

Saving and Restoring Computations 

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. 

Using the Package in the Classroom 

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. 

> R[SI] := evalf(ScientificConstants[Constant]('R')); 1
 

Typesetting:-mprintslash([R[SI] := 8.3144598], [8.3144598])
 

We can define the Maxwell-Boltzmann distribution with Maple. 

> f := proc (v, M, T, R) options operator, arrow, function_assign; `*`(Pi, `*`(`^`(2, `/`(1, 2)), `*`(`^`(`/`(`*`(M), `*`(Pi, `*`(R, `*`(T)))), `/`(3, 2)), `*`(`^`(v, 2), `*`(exp(`+`(`-`(`/`(`*`(`/`(1, ...
 

Typesetting:-mprintslash([f := proc (v, M, T, R) options operator, arrow, function_assign; `*`(Pi, `*`(`^`(2, `/`(1, 2)), `*`(`^`(`/`(`*`(M), `*`(Pi, `*`(R, `*`(T)))), `/`(3, 2)), `*`(`^`(v, 2), `*`(e...
 

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 := int(f(v, .1, 273, R[SI]), v = 0 .. infinity); 1
 

Typesetting:-mprintslash([P_total := .9999999999], [.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. 

> Explore(plot(f(v, M, T, R[SI]), v = 0 .. 2000, axes = boxed, labels = ['v(`/`(`*`(m), `*`(s)))', 'f(v)'], color = blue, thickness = 2), T = 100 .. 1000, M = 0.1e-1 .. 1); 1
Explore(plot(f(v, M, T, R[SI]), v = 0 .. 2000, axes = boxed, labels = ['v(`/`(`*`(m), `*`(s)))', 'f(v)'], color = blue, thickness = 2), T = 100 .. 1000, M = 0.1e-1 .. 1); 1
 

Embedded component 

T 

Embedded component 

   100 

M 

Embedded component 

0.01 

 

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 := [[
mol := [[
mol := [[
mol := [[
mol := [[
mol := [[
mol := [[
mol := [[
mol := [[
mol := [[
 

Typesetting:-mprintslash([mol := [[
Typesetting:-mprintslash([mol := [[
Typesetting:-mprintslash([mol := [[
Typesetting:-mprintslash([mol := [[
Typesetting:-mprintslash([mol := [[
Typesetting:-mprintslash([mol := [[
Typesetting:-mprintslash([mol := [[
Typesetting:-mprintslash([mol := [[
Typesetting:-mprintslash([mol := [[
Typesetting:-mprintslash([mol := [[
 


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 := DensityFunctional(mol, basis =
 

Typesetting:-mrow(Typesetting:-mi(
Typesetting:-mrow(Typesetting:-mi(
Typesetting:-mrow(Typesetting:-mi(
Typesetting:-mrow(Typesetting:-mi(
Typesetting:-mrow(Typesetting:-mi(
 

The ground-state energy is given by

> data[e_tot]; 1

Typesetting:-mprintslash([-1259.07521546698], [HFloat(-1259.0752154669813)])

(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. 

> data[mo_occ]; 1
 

Typesetting:-msub(Typesetting:-mi(
 

> DensityPlot3D(mol, data, basis =
 

Plot_2d

> DensityPlot3D(mol, data, basis =

Plot_2d

The energies of these two orbitals are

> data[mo_energy][57 .. 58]; 1

Typesetting:-mfenced(Typesetting:-mrow(Typesetting:-mtable(Typesetting:-mtr(Typesetting:-mtd(Typesetting:-mrow(Typesetting:-mo(

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