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 2021, this toolbox has significant new features and enhancements: (1) a new command that provides a treasure trove of information about molecules−a molecular dictionary, (2) a new method for molecules that further enhances the package's suite of electronic structure solvers, (3) a new command for purifying density matrices with applications to the mitigation of errors in quantum computing, (4) a new optional parameter for the plotting of molecular orbitals that allows for customized colors, (5) new optional parameters for the generation of molecular integrals that support arbitrary molecular orbitals and active spaces, (6) a new addition to the ≈30 builtin lessons for classroom learning in undergraduate-to-graduate chemistry and physics, (7) updates to the interactive Maplet, and (8) additional enhancements to many methods and commands throughout the package.
Note that the Maple Quantum Chemistry Toolbox (QCT) is required in order to execute the examples in this worksheet.
Molecular Dictionary
Computing with Molecules
Error Mitigation in Quantum Computing
3-Dimensional Density Plots
Molecular Integrals
Using the Package in the Classroom
A molecular dictionary appears in QCT 2021 with access to online data for more than 50,000 molecules. The command MolecularDictionary accepts the name of a molecule as a Maple string or CID integer. For example, we retrieve the entries for the following two molecules: nitrogen ("dinitrogen") gas and the Cys-Ser peptide. Before we begin we load the QuantumChemistry package,
withQuantumChemistry;
AOLabels,ActiveSpaceCI,ActiveSpaceSCF,AtomicData,BondAngles,BondDistances,Charges,ChargesPlot,ContractedSchrodinger,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,MolecularDictionary,MolecularGeometry,NuclearEnergy,NuclearGradient,OscillatorStrengths,Parametric2RDM,PlotMolecule,Populations,Purify2RDM,RDM1,RDM2,RTM1,ReadXYZ,Restore,Save,SaveXYZ,SearchBasisSets,SearchFunctionals,SkeletalStructure,Thermodynamics,TransitionDipolePlot,TransitionDipoles,TransitionOrbitalPlot,TransitionOrbitals,Variational2RDM,VibrationalModeAnimation,VibrationalModes,Video
For nitrogen gas we have
MolecularDictionarynitrogen;
Nitrogen: Dinitrogen is an elemental molecule consisting of two trivalently-bonded nitrogen atoms. It has a role as a member of food packaging gas and a food propellant. It is a diatomic nitrogen, a gas molecular entity and an elemental molecule. It is a conjugate base of a diazynium. -ChEBI Nitrogen: Nitrogen is an element with atomic symbol N, atomic number 7, and atomic weight 14.01. -NCI Thesaurus (NCIt)
while for the Cys-Ser peptide we have
MolecularDictionaryCys-Ser;
Cys-Ser: Cys-Ser is a dipeptide composed of L-cysteine and L-serine joined by a peptide linkage. It has a role as a metabolite. It is a dipeptide, a secondary carboxamide, a primary alcohol, a primary amino compound, a thiol and a carboxylic acid. It derives from a L-cysteine and a L-serine. -ChEBI
Similarly, we can explore many, many other molecules. What about the molecule ozone that is sometimes generated before a thunderstorm?
MolecularDictionaryozone;
Ozone: Ozone is an elemental molecule with formula O3. An explosive, pale blue gas (b.p. -112###) that has a characteristic, pleasant odour, it is continuously produced in the upper atmosphere by the action of solar ultraviolet radiation on atmospheric oxygen. It is an antimicrobial agent used in the production of bottled water, as well as in the treatment of meat, poultry and other foodstuffs. It has a role as a member of greenhouse gas, a disinfectant, a tracer, an electrophilic reagent, a mutagen, an oxidising agent and an antiseptic drug. It is a member of reactive oxygen species, an elemental molecule, a triatomic oxygen and a gas molecular entity. -ChEBI
Or the molecule ascorbic acid (Vitamin C)?
MolecularDictionaryascorbic acid;
Ascorbic acid: L-ascorbic acid is the L-enantiomer of ascorbic acid and conjugate acid of L-ascorbate. It has a role as a water-soluble vitamin, a vitamin C, a coenzyme, a flour treatment agent, a food antioxidant, a plant metabolite and a cofactor. It is a conjugate acid of a L-ascorbate. It is an enantiomer of a D-ascorbic acid. -ChEBI Ascorbic acid: Ascorbic Acid is a natural water-soluble vitamin (Vitamin C). Ascorbic acid is a potent reducing and antioxidant agent that functions in fighting bacterial infections, in detoxifying reactions, and in the formation of collagen in fibrous tissue, teeth, bones, connective tissue, skin, and capillaries. Found in citrus and other fruits, and in vegetables, vitamin C cannot be produced or stored by humans and must be obtained in the diet. (NCI04) -NCI Thesaurus (NCIt)
Or that important molecule in coffee?
MolecularDictionarycaffeine;
Caffeine: Caffeine is a trimethylxanthine in which the three methyl groups are located at positions 1, 3, and 7. A purine alkaloid that occurs naturally in tea and coffee. It has a role as a central nervous system stimulant, an EC 3.1.4.* (phosphoric diester hydrolase) inhibitor, an adenosine receptor antagonist, an EC 2.7.11.1 (non-specific serine/threonine protein kinase) inhibitor, a ryanodine receptor agonist, a fungal metabolite, an adenosine A2A receptor antagonist, a psychotropic drug, a diuretic, a food additive, an adjuvant, a plant metabolite, an environmental contaminant, a xenobiotic, a human blood serum metabolite, a mouse metabolite and a mutagen. It is a purine alkaloid and a trimethylxanthine. -ChEBI Caffeine: Caffeine is a methylxanthine alkaloid found in the seeds, nuts, or leaves of a number of plants native to South America and East Asia that is structurally related to adenosine and acts primarily as an adenosine receptor antagonist with psychotropic and anti-inflammatory activities. Upon ingestion, caffeine binds to adenosine receptors in the central nervous system (CNS), which inhibits adenosine binding. This inhibits the adenosine-mediated downregulation of CNS activity; thus, stimulating the activity of the medullary, vagal, vasomotor, and respiratory centers in the brain. This agent also promotes neurotransmitter release that further stimulates the CNS. The anti-inflammatory effects of caffeine are due the nonselective competitive inhibition of phosphodiesterases (PDEs). Inhibition of PDEs raises the intracellular concentration of cyclic AMP (cAMP), activates protein kinase A, and inhibits leukotriene synthesis, which leads to reduced inflammation and innate immunity. -NCI Thesaurus (NCIt)
The database also contains many more technical molecules such as the antiviral drug remdesivir that was recently repurposed to treat Covid-19.
MolecularDictionaryremdesivir;
Remdesivir: Remdesivir is a carboxylic ester resulting from the formal condensation of the carboxy group of N-[(S)-{[(2R,3S,4R,5R)-5-(4-aminopyrrolo[2,1-f][1,2,4]triazin-7-yl)-5-cyano-3,4-dihydroxytetrahydrofuran-2-yl]methoxy}(phenoxy)phosphoryl]-L-alanine with the hydroxy group of 2-ethylbutan-1-ol. A broad-spectrum antiviral prodrug with potent in vitro antiviral activity against a diverse panel of RNA viruses such as Ebola virus, MERS-CoV and SARS-CoV. It is currently in Phase III clinical trials for the treatment of Covid-19 in adults. It has a role as an antiviral drug, a prodrug and an anticoronaviral agent. It is a carboxylic ester, a pyrrolotriazine, a nitrile, a phosphoramidate ester, a C-nucleoside and an aromatic amine. It derives from a GS-441524. -ChEBI Remdesivir: Remdesivir is a prodrug of an adenosine triphosphate (ATP) analog, with potential antiviral activity against a variety of RNA viruses. Upon administration, remdesivir, being a prodrug, is metabolized into its active form GS-441524. As an ATP analog, GS-441524 competes with ATP for incorporation into RNA and inhibits the action of viral RNA-dependent RNA polymerase. This results in the termination of RNA transcription and decreases viral RNA production. -NCI Thesaurus (NCIt)
The new MolecularDictionary command complements the existing commands MolecularGeometry and SkeletalStructure for retrieving molecular geometries and skeletal structures. For example, here is the skeletal structure for remdesivir
SkeletalStructureremdesivir;
and its molecular geometry
remdesivir ≔ MolecularGeometryremdesivir;
remdesivir≔P,0.71860000,0.10970000,−0.76660000,O,−2.83790000,2.01490000,−0.79800000,O,−4.60280000,3.10730000,2.14720000,O,−3.33010000,4.75570000,0.36760000,O,−0.07640000,1.52190000,−0.93430000,O,−0.38890000,−0.77740000,0.01860000,O,1.05700000,−0.46580000,−2.11880000,O,5.56710000,0.78540000,0.04450000,O,4.25640000,−0.39040000,1.53130000,N,−3.90140000,−0.59850000,−0.58280000,N,−3.34010000,−0.39130000,−1.79850000,N,1.97840000,0.37370000,0.27670000,N,−6.06950000,2.93730000,−1.04000000,N,−3.45180000,−2.78050000,−2.09530000,N,−4.30780000,−4.26940000,−0.51260000,C,−3.95980000,1.84290000,0.10590000,C,−3.52950000,2.58210000,1.38960000,C,−2.58500000,3.65710000,0.88790000,C,−1.87900000,2.94460000,−0.25310000,C,−4.21790000,0.35750000,0.34660000,C,−0.65550000,2.15550000,0.19440000,C,−5.14180000,2.45600000,−0.53610000,C,−4.76810000,−0.28070000,1.43770000,C,−4.24220000,−1.83550000,−0.09990000,C,−4.78420000,−1.67020000,1.15620000,C,−3.99470000,−2.99540000,−0.92270000,C,3.22850000,0.90630000,−0.26190000,C,−3.14980000,−1.48650000,−2.47870000,C,7.99230000,0.90260000,0.02990000,C,6.73740000,0.33420000,0.71240000,C,8.01980000,0.61290000,−1.48750000,C,9.24690000,0.42250000,0.78550000,C,4.37400000,0.34540000,0.55990000,C,3.21180000,2.42720000,−0.19330000,C,−0.18090000,−2.02330000,0.54580000,C,8.09250000,−0.86050000,−1.86400000,C,10.54730000,1.00690000,0.25290000,C,0.74520000,−2.88410000,−0.04300000,C,−0.89510000,−2.43210000,1.67190000,C,0.95730000,−4.15390000,0.49430000,C,−0.68330000,−3.70160000,2.20940000,C,0.24290000,−4.56260000,1.62060000,H,−2.96420000,1.90360000,2.04220000,H,−1.91910000,4.04370000,1.66480000,H,−1.60480000,3.62590000,−1.06620000,H,−0.92710000,1.40590000,0.94310000,H,0.09090000,2.82470000,0.63360000,H,−4.22020000,3.62550000,2.87570000,H,−5.12360000,0.19600000,2.34060000,H,−2.69180000,5.39520000,0.00820000,H,−5.15100000,−2.45610000,1.80150000,H,1.77080000,0.79490000,1.18130000,H,3.37230000,0.56210000,−1.29190000,H,7.93990000,1.99660000,0.12380000,H,−2.69460000,−1.38650000,−3.47770000,H,6.70450000,0.67470000,1.75490000,H,6.76100000,−0.76200000,0.71010000,H,8.86860000,1.13460000,−1.94470000,H,7.12730000,1.04590000,−1.95600000,H,9.15420000,0.70950000,1.84080000,H,9.31030000,−0.67170000,0.76940000,H,−4.72700000,−4.45630000,0.39210000,H,−4.12310000,−5.06830000,−1.10970000,H,2.37960000,2.84230000,−0.77100000,H,4.13510000,2.85010000,−0.60420000,H,3.12470000,2.78400000,0.83960000,H,7.22460000,−1.41360000,−1.49380000,H,8.99880000,−1.33170000,−1.47330000,H,8.10960000,−0.96480000,−2.95390000,H,10.77610000,0.63290000,−0.74910000,H,11.37950000,0.72520000,0.90630000,H,10.50160000,2.09980000,0.21570000,H,1.31380000,−2.60450000,−0.92210000,H,−1.61490000,−1.76410000,2.13670000,H,1.67730000,−4.82510000,0.03530000,H,−1.23940000,−4.01950000,3.08640000,H,0.40750000,−5.55130000,2.03880000
which we can plot with the PlotMolecule command
PlotMoleculeremdesivir;
With the dictionary entry, the skeletal structure, and the ball-and-stick plot we have all of the necessary background information to begin molecular calculations like the ones in the next section.
In the QCT 2021 we implement a new method for the ground and excited states of molecules. The Schrödinger equation, developed by Erwin Schrödinger in 1925, is the key equation for solving quantum systems in non-relativistic physics. However, its exact solution scales exponentially with the number of electrons in the quantum system. To deal with this inefficiency, we solve a projection (or contraction) of the Schrödinger equation onto the space of only two electrons, known as a contracted Schrödinger equation. This contraction is especially useful for improving computational efficiency for strongly correlated molecules−molecules that deviate significantly from the mean-field (Hartree-Fock) picture of one electron in the field of the remaining electrons−for which many standard methods are often inaccurate. Let us consider the peptide of the two amino acids L-cysteine and L-serine whose geometry we can import with the MolecularGeometry command
mol ≔ MolecularGeometryCys-Ser;
mol≔S,−2.29090000,2.10530000,−1.07690000,O,2.52680000,−1.19180000,−1.63380000,O,−1.06460000,−0.04400000,1.69560000,O,2.66940000,1.27060000,1.34030000,O,1.68720000,1.60900000,−0.67760000,N,0.14240000,−0.64180000,−0.19310000,N,−3.00870000,−1.66870000,0.29030000,C,1.44960000,−0.55430000,0.41400000,C,−2.27570000,−0.56710000,−0.33790000,C,−1.02220000,−0.37450000,0.51210000,C,2.44060000,−1.50010000,−0.24810000,C,−3.16790000,0.67650000,−0.36510000,C,1.91400000,0.88100000,0.28020000,H,1.35260000,−0.79450000,1.47910000,H,−2.01660000,−0.86110000,−1.36160000,H,0.09180000,−0.84710000,−1.18670000,H,2.11130000,−2.53950000,−0.14810000,H,3.43800000,−1.39630000,0.19140000,H,−3.48660000,0.95790000,0.64500000,H,−4.06810000,0.49580000,−0.96240000,H,−2.46140000,−2.52630000,0.22140000,H,−3.87150000,−1.84260000,−0.22420000,H,3.16370000,−1.81370000,−2.02500000,H,−3.30210000,2.97810000,−0.95790000,H,2.98800000,2.19320000,1.24240000
and plot with the PlotMolecule command
PlotMoleculemol;
In this example we showcase the use of two active spaces, an active space of 6 electrons in 6 orbitals to account for any potential static (strong) electron correlation and a larger active space to account for additional (dynamic) electron correlation. (This second active space is chosen to be [16,16] to make the calculation quick but much larger second active spaces are possible.)
data_peptide ≔ContractedSchrodingermol, active=6,6,active_cse=16,16;
In the table of computed properties we observe that the peptide has a nonzero dipole moment
data_peptidedipole;
0.95773786−2.37658291−1.59751383
The direction of the dipole moment can be readily visualized with the DipolePlot command
DipolePlotmol, method=ContractedSchrodinger, active=6,6, active_cse=16,16;
We can also compute the atomic charges which are closely associated with the dipole moment
seqmoli1,data_peptidechargesi, i=1..nopsmol;
S,0.10504005,O,−0.30624021,O,−0.20852064,O,−0.28630119,O,−0.26480265,N,−0.37476210,N,−0.39527767,C,0.04398837,C,0.03851104,C,0.26632682,C,0.00879806,C,−0.17897052,C,0.29797645,H,0.09534526,H,0.04389074,H,0.20310443,H,0.05844607,H,0.06152166,H,0.07400052,H,0.05447756,H,0.15148599,H,0.15137748,H,0.18850391,H,−0.04234736,H,0.21442792
Only one of the carbons is negatively charged and only two of the carbons are fairly positively charged which we can also see from a ChargesPlot
ChargesPlotmol, method=ContractedSchrodinger, active=6,6, active_cse=16,16;
where the positively charged atoms are shown in red and the negatively charged atoms are shown in blue.
All of the one- and two-electron properties of molecules are computable from the two-electron reduced density matrix (2-RDM). The 2-RDM of a molecule obeys a set of stringent constraints that are necessary for it to represent an N-electron wave function−constraints known as N-representability conditions. A 2-RDM that is computed by an approximate method or measured in the presence of noise as on a quantum computer can violate some of these conditions. In QCT 2021 we introduce a new algorithm through the command Purify2RDM that purifies the 2-RDM to obey an important set of N-representability conditions. Purify2RDM uses an advanced type of optimization, known as semidefinite programming, to project the 2-RDM onto the set of approximately N-representable 2-RDMs. The command is especially useful in molecular simulations on quantum computers where noise generates errors in the measured 2-RDMs. The correction of errors in quantities measured on a quantum computer is known as error mitigation. Here we demonstrate the command by using it to purify an approximate 2-RDM of hydrogen fluoride generated from second-order many-body perturbation (MP2) theory. First, we generate the data from applying the MP2 command to the geometry of hydrogen fluoride where we include the keyword parameter return_rdm="rdm1_and_rdm2" to return the 2-RDM.
hf ≔ MolecularGeometryhydrogen fluoride;
hf≔F,0,0,0,H,0.43580000,−0.14770000,−0.81880000
data_hf ≔ MP2hf, return_rdm=rdm1_and_rdm2:
Second, we reshape the data stored in a 4-index tensor into a 2-index matrix
D2 ≔ Matrix36,36, shape=symmetric:D2.. ≔ ArrayTools:-Permutedata_hfrdm2, 1,3,2,4;
and compute the eigenvalues of the 2-RDM with the Eigenvalues command in the LinearAlgebra package
eigenvals ≔ LinearAlgebra:-EigenvaluesD2;
Notice that the lowest eigenvalue of the approximate 2-RDM is negative (-0.0143), but because each eigenvalue represents the probability of being in a two-electron quantum state, all of the eigenvalues of the 2-RDM should be non-negative! Finally, we use the Purify2RDM command to correct this violation as well as more subtle violations of the physical requirements for the 2-RDM
data_hf_purified ≔ Purify2RDMdata_hfrdm2;
To check the correction of the 2-RDM, as before we rearrange the outputted 4-index tensor into a 2-index matrix
D2p ≔ Matrix36,36, shape=symmetric:D2p.. ≔ ArrayTools:-Permutedata_hf_purifiedrdm2, 1,3,2,4;
which we diagonalize with the LinearAlgebra:-Eigenvalues command
eigenvals ≔ LinearAlgebra:-EigenvaluesD2p;
All of the eigenvalues, we find, are now non-negative.
The command DensityPlot3D receives new optional keywords bonds and colors that control the appearance of chemical bonds and the color of the orbital densities. Consider the diatomic molecule hydrogen fluoride
After performing an electronic structure calculation such as its solution by the new ContractedSchrodinger command
data_hf ≔ ContractedSchrodingerhf, active=3,6;
we can visualize the orbitals with DensityPlot3D. In QCT 2020 the density plot always uses green and purple colors for the orbitals
DensityPlot3Dhf, data_hf, orbitalindex=5, densitycutoff=0.0005, orientation=−97,−4,−75;
Now we can use the keyword color to change the default colors to any color supported in Maple's ColorTools package. For example,
DensityPlot3Dhf, data_hf, orbitalindex=5, densitycutoff=0.0005, colors=Nautical Red,Nautical Light Blue;
The QCT makes advanced features like obtaining the molecular integrals as simple as a single command MOIntegrals. In QCT 2021 we add two new optional keywords initial_mo and active to the MOIntegrals command that makes generating molecular integrals in any MO basis set with or without an active space as easy as possible. Let's consider the molecule citric acid whose skeletal structure, geometry, and dictionary entry are
SkeletalStructurecitric acid;
citricacid ≔ MolecularGeometrycitric acid;
citricacid≔O,0.02960000,0.20950000,1.60690000,O,−0.85580000,1.79790000,−1.49620000,O,−1.35540000,2.26890000,0.66220000,O,−2.39080000,−2.30460000,−0.78710000,O,3.53510000,−0.40330000,−0.92730000,O,−2.55590000,−0.62650000,0.73170000,O,2.54380000,−1.12450000,0.98350000,C,0.10070000,0.36200000,0.19410000,C,−0.49050000,−0.87910000,−0.49710000,C,1.53740000,0.70960000,−0.22660000,C,−0.76500000,1.57730000,−0.15870000,C,−1.90980000,−1.23050000,−0.11540000,C,2.57670000,−0.35680000,0.03000000,H,0.11200000,−1.75550000,−0.22860000,H,−0.45740000,−0.76140000,−1.58660000,H,1.88420000,1.60800000,0.30000000,H,1.55890000,0.93610000,−1.30040000,H,0.46200000,0.97710000,2.01860000,H,−1.42930000,2.56690000,−1.70150000,H,−3.30780000,−2.52520000,−0.51720000,H,4.19140000,−1.11160000,−0.75460000
and
MolecularDictionarycitric acid;
Citric acid: Citric acid is a tricarboxylic acid that is propane-1,2,3-tricarboxylic acid bearing a hydroxy substituent at position 2. It is an important metabolite in the pathway of all aerobic organisms. It has a role as a food acidity regulator, a chelator, an antimicrobial agent and a fundamental metabolite. It is a conjugate acid of a citrate(1-) and a citrate anion. -ChEBI Citric acid: Anhydrous Citric Acid is a tricarboxylic acid found in citrus fruits. Citric acid is used as an excipient in pharmaceutical preparations due to its antioxidant properties. It maintains stability of active ingredients and is used as a preservative. It is also used as an acidulant to control pH and acts as an anticoagulant by chelating calcium in blood. -NCI Thesaurus (NCIt)
We can generate the molecular orbitals (MOs) in a 26 electrons in 26 orbitals ([26,26]) active space−the 13 highest occupied orbitals (26 electrons) and the 13 lowest unoccupied orbitals−with the active keyword
data ≔ MOIntegralscitricacid, active=26,26;
These molecular integrals are in the basis set of the Hartree-Fock MOs. We can easily use any set of MOs with the keyword initial_mo. For example, we can run a density functional theory (DFT) calculation and use the DFT MOs
data_dft ≔ DensityFunctionalcitricacid;
The initial_mo keyword is assigned to a list of the MO coefficients and MO symmetries
data ≔ MOIntegralscitricacid, active=26,26, initial_mo=data_dftmo_coeff,data_dftmo_symmetry;
Voila! We have the MO integrals with respect to the DFT MOs in a [26,26] active space.
The Maple Quantum Chemistry Toolbox includes approximately 30 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. Detailed lesson plans and curricula are provided for Introductory (General) Chemistry, Physical Chemistry (Quantum Mechanics and Thermodynamics), Thermodynamics (Physics), Quantum Mechanics (Physics), Computational Chemistry, and Quantum Chemistry as well as Advanced Placement (AP) and International Baccalaureate (IB) chemistry courses. Topics include atomic structure, chemical bonding, the Maxwell-Boltzmann distribution, heat capacity, enthalpy, entropy, free energy, particle-in-a-box, vibrational normal modes, infrared spectroscopy, as well as advanced electronic structure methods. Use of the QCT in the classroom is described in a recent paper in J. Chem. Ed. QCT 2021 includes a new lesson for Physical Chemistry and Undergraduate Quantum Mechanics on Vibrational Motion and the Harmonic Oscillator. In the following two subsections we show snippets of the lesson
Reduced mass
Here we set the variables A and B to the strings of the atoms in the diatomic molecule (By changing these values, you can use the worksheet to treat other diatomic molecules!)
A ≔ H;B ≔ F;
A≔H
B≔F
We use the AtomicData command in the Quantum Chemistry package to obtain the masses as well as other data
dataA ≔ AtomicDataA;
dataA≔tableionizationenergy=13.59840000eV,name=hydrogen,meltingpoint=13.81000000K,electronegativity=2.10000000,atomicweight=1.00794000amu,atomicnumber=1,boilingpoint=20.28000000K,names=hydrogen,symbol=H
dataB ≔ AtomicDataB;
dataB≔tableionizationenergy=17.42280000eV,name=fluorine,electronaffinity=3.40119000eV,meltingpoint=53.53000000K,electronegativity=3.98000000,atomicweight=18.99840320amu,atomicnumber=9,boilingpoint=85.03000000K,names=fluorine,symbol=F
Set the masses
mA ≔ dataAatomicweight;
mA≔1.00794000amu
mB ≔ dataBatomicweight;
mB≔18.99840320amu
Compute the reduced mass in amu
mu0 ≔ mA⋅mBmA+mB;
μ0≔0.95715895amu
Convert the reduced mass from amu to kg
mu ≔ convertmu0,units,'kg';
μ≔1.5894009210−27kg
Bond length
To compute the equilibrium bond length, we select a set of bond distances from the roots of the sixth-order Chebyshev polynomial that are suitable for interpolation
bonds ≔ mapx → x/5+1.05, fsolveexpandChebyshevT6,x;
bonds≔0.85681483,0.90857864,0.99823619,1.10176381,1.19142136,1.24318516
We define a list of molecular geometries with each geometry corresponding to one of the bond distances
molecules ≔mapR→ A,0,0,0,B,0,0,R, bonds;
molecules≔H,0,0,0,F,0,0,0.85681483,H,0,0,0,F,0,0,0.90857864,H,0,0,0,F,0,0,0.99823619,H,0,0,0,F,0,0,1.10176381,H,0,0,0,F,0,0,1.19142136,H,0,0,0,F,0,0,1.24318516
The energies for each geometry may be then readily computed with the Energy command in the Quantum Chemistry package.
energies ≔ mapEnergy,molecules,basis=cc-pVDZ;
energies≔−100.01686137,−100.01964383,−100.01018182,−99.98693789,−99.96201956,−99.94683591
We use polynomial interpolation to generate a polynomial in the bond distance R
pes ≔ interpbonds,energies,R;
pes≔−2.84449888R5+17.02510690R4−41.28045230R3+50.73972465R2−31.33797004R−92.31177999
The potential energy surface (curve) can be plotted
plotpesR, R=bonds1..bonds−1, axes=boxed, labels='R','E', color=blue, thickness=3;
Finally, we differentiate the potential energy curve with respect to R and set the derivative to zero.
eq ≔ diffpes, R = 0;
eq≔−14.22249441R4+68.10042762R3−123.84135691R2+101.47944931R−31.33797004=0
Solving the resulting equation yields the equilibrium bond length
R_eq ≔ fsolveeq, R=bonds1..bonds−1;
R_eq≔0.90161288
Download Help Document
