Quantum Chemistry Package Tutorial Worksheet
Copyright (c) RDMCHEM LLC 2023
Load the Package
Definition of a Molecule
The QuantumChemistry package is an environment for computing and visualizing the quantum energies and properties of many-electron atoms and molecules.
In this tutorial we will teach the user to use the following basic features of the package:
Definition of a molecule's geometry from an interface, an XYZ file, or a web database containing more than 96 million molecules (see section below)
Computation of a molecule's energies and properties from advanced wavefunction, density functional, and reduced density matrix methods (see section below)
Visualization of molecular geometries, densities, and vibrations through 3D interactive plots and animations (see section below)
Launch of an
interactive Maplet (see section below) for a graphical user interface to the package commands for defining, computing, and visualizing molecular data
Integration of lessons and curricula in chemistry and physics (see section below) which can be used for formal courses as well as self study
Launch of the Maple help pages for the QuantumChemistry package for detailed information on each command (see section below)
See the What'sNew2023 help page for examples of what's new in the QuantumChemistry Package 2023.
We begin by loading the commands of the Quantum Chemistry package using the Maple with command
We recommend setting the global variable Digits to 15 (double precision) rather than 10 (default):
A molecule's geometry can be defined in the worksheet interface in three ways:
entering the geometry as a Maple list of lists
reading the geometry from an XYZ file
retrieving the geometry from a web database
Maple List of Lists
Let us first enter the geometry of the molecule hydrogen fluoride HF by using a Maple list of lists. Each list has four elements:
a string of the atomic symbol of the atom, i.e. "H" or "F"
the x Cartesian coordinate of the atom
the y Cartesian coordinate of the atom
the z Cartesian coordinate of the atom
The default unit of distance in the package is Angstroms (1 Angstrom= 10−10 m), which we employ in this worksheet.
Second, we can retrieve the geometry of a molecule from a web database containing more than 96 million molecules. The database is especially useful for molecules with more than a few atoms. To perform the search, we use the command MolecularGeometry (the computer must be connected to the web).
Second, we can read the geometry of a molecule from an XYZ file, that is a file containing the Cartesian coordinates of each atom. For example, the XYZ file of hydrogen fluoride has the following format:
# Comment Line: Note that we give the total number of atoms in the first line, a comment in the second line, and then the Cartesian coordinates of the atoms in the remaining lines
H 0 0 0
F 0 0 0.95
To save an XYZ file, we use the command SaveXYZ to make an XYZ file from the Maple list of lists. To read the resulting XYZ file, we use the command ReadXYZ.
Invoking ReadXYZ without a parameter ReadXYZ() opens a file dialogue for selecting an XYZ file. Invoking ReadXYZ with a string parameter opens the file associated with the string name. If the file is not present in Maple's current working directory, the string must contain the full path to the file or Maple's current working directory must modified by using the currentdir command. For example, the XYZ file for porphin "porphin.xyz" can be read from the data folder of the package. First we construct a full path to the file
file ≔ FileTools:-JoinPathkerneloptstoolboxdir=QuantumChemistry,data,porphin.xyz;
Then we call ReadXYZ to read the XYZ file into the QuantumChemistry package's molecular geometry format
porphin ≔ ReadXYZfile;
Computations of molecules can be performed in one of two ways:
Maple commands based on the energy or property to be computed
Maple commands based on the quantum chemistry method to be employed
The QuantumChemistry package contains commands to compute directly the following energies and properties:
molecular orbital coefficients
molecular orbital occupations
molecular orbital energies
atomic-orbital populations and charges
transition dipole moments
natural transition orbitals
1-electron reduced density matrix (1-RDM)
2-electron reduced density matrix (2-RDM)
1-electron reduced transition matrices (1-RTMs)
Unless otherwise noted, all energies are in atomic units (a.u.).
The energy of hydrogen fluoride HF, for example, can be directly computed from the Hartree-Fock method from the command Energy
Each energy/property command can take all of the options, specified by optional keywords, available for the given method. For instance, we can specify a larger basis set with the basis keyword
We can select a different method for computing the energy
The nuclear energy from the repulsion of the nuclei in the molecule can be computed with the command NuclearEnergy
The energies of the molecular orbitals can be computed with the command MOEnergies
The occupations of the molecular orbitals can be computed with the command MOOccupations
The occupations being 2 (filled) or 0 (unfilled) reflect that the Hartree-Fock method does not correlate the orbitals describing the electrons. Using the command MOOccupations with a correlated method like the variational 2-RDM method generates fractional occupations, that is occupations between 0 and 2
The 1-electron reduced density matrix (1-RDM) can be calculated from the variational 2-RDM method with the command RDM1
The eigenvalues of the 1-RDM are the occupations of its eigenvectors, known as natural orbitals. It is the occupations of the natural orbitals that are returned by the command command MOOccupations. To confirm that this is true for the variational 2-RDM method, we can calculate the eigenvalues of the 1-RDM with the Eigenvalues command in the LinearAlgebra package. First, we ensure that dm1 is a symmetric matrix which will have real eigenvalues,
and then we compute the eigenvalues
The eigenvalues reveal the number of electrons in each natural orbital. The deviation of the eigenvalues from 0 and 2 reveals the presence of electron correlation. With a knowledge of the 1-RDM other important one-electron properties can be computed.
The Energy command is especially useful for generating potential energy curves (and surfaces) with just a few lines of code. For example, the following five lines generate the potential energy of hydrogen fluoride from the Hartree-Fock method with the internuclear distance ranging from 0.8 Angstroms to 2.0 Angstroms in 0.05 increments
A similar point plot can be generated from the full solution of the Schrodinger equation in the given basis set, known as full configuration interaction
Finally, a plot can be generated from the variational 2-RDM method in which the correlated 2-RDM is directly computed without computing the many-electron wavefunction
The three plots can be readily display together to reveal that the energies from the variational 2-RDM calculations exactly match those from the full configuration interaction. .
In contrast to the 2-RDM and full configuration methods, the Hartree-Fock method yields a potential energy curve whose minimum energy occurs at an internuclear distance R that is too short and whose curvature at the minimum energy is too large.
Optimization of a molecule's geometry to reach a first-order critical point such as a minimum or a transition state can be performed with the command GeometryOptimization. Furthermore, the vibrational frequencies of a molecule in the normal-mode approximation can be computed with the command VibrationalModes. For example, we can optimize the internuclear distance of the hydrogen fluoride molecule by the Hartree-Fock method
Because the fluorine is located at the origin, the bond distance can be computed as follows
bonddistance_hf ≔ BondDistanceshf_hf,1,2;
Once the optimized geometry is known, the vibrational frequency of the bond can be determined to be 4459 cm−1 from
Similarly, we can find optimize the internuclear distance by the variational 2-RDM method, using the D and Q N-representability conditions and freezing the core orbitals,
In this case the equilibrium bond distance, the distance at which the minimum energy occurs, is
The vibrational frequency of the bond can be determined in cm−1 from
As in the plots above, we observe that the treatment of electron correlation by the variational 2-RDM method causes the equilibrium bond length to increase to 0.993 from 0.955 Angstroms and the vibrational frequency in wavenumbers to decrease significantly from 4459 cm−1.
The QuantumChemistry package can solve the Schrodinger equation of atoms and molecules by several different methods. The methods can be employed in the Energy and other commands of the previous section, and they can also be called directly. Supported methods include:
Hartree-Fock (HF) methods
Density Functional Theory (DFT) methods
Variational Two-electron Reduced Density Matrix (2-RDM) method
Parametric Two-electron Reduced Density Matrix (2-RDM) method
Anti-Hermitian Contracted Schrödinger Equation (CSE) method
Full Configuration Interaction (FCI) method
Second-order Many-body Perturbation Theory (MP2) method
Coupled Cluster Singles-Doubles (CCSD) method
Complete-Active-Space Configuration Interaction (CAS-CI) method
Complete-Active-Space Self-Consistent-Field (CAS-SCF) method
Time-dependent Hartree-Fock (TDHF) method (via the HartreeFock command)
Configuration-Interaction Singles (CIS) method (via the HartreeFock command)
Time-dependent Density Functional Theory (TDDFT) method (via the DensityFunctional command)
Tamm-Dancoff Approximation (TDA) method (via the DensityFunctional command)
The Hartree-Fock method, for example, can be invoked by the HartreeFock command
As described in the help page for the HartreeFock command, which can be opened by clicking on the hyperlink, each method has the molecule's geometry as a Maple list of lists as its first argument. Additional keyword arguments can be provided to control the charge, spin, symmetry, and basis set as well as technical details of the solution method. Complete information about these options is provided on the help page. Running a method's command returns a Maple table containing entries that correspond to generated data such as the total electronic energy, the coefficients of the molecular orbitals in terms of the atomic orbitals, the electron occupations of the molecular orbitals, the energies of the orbitals, the symmetries of the orbitals, the 1- and 2-electron reduced density matrices, the name of the point-group symmetry, and a parameter indicating convergence (set to one) or non-convergence (set to zero). The total energy from the data table can be accessed as follows:
Similarly, the energies of the molecular orbitals can be retrieved as follows
Rerunning the command takes a negligible amount of CPU time because the result is stored in a cache table.
The variational 2-RDM method, as another example, can be invoked by the Variational2RDM command
The basis set of the molecule can be changed using the basis keyword. The basis keyword supports a wide range of basis sets including mixed basis sets in which each atom type has its own basis set as well as effective core potential (ECP) basis sets for heavy atoms (see Basis for details).
The charge of the molecule can be changed using the charge and spin keywords
From the two energies we can estimate the ionization energy, the energy to remove an electron as 0.531 a.u.
which (at least at this level of approximation) is slightly greater than the energy to ionize the hydrogen atom.
Density functional theory (DFT) can be invoked by the DensityFunctional command
As with the Hartree-Fock method the charge of the molecule can be changed using the charge and spin keywords
From DFT with the default B3LYP functional the energy to remove an electron is 0.596 a.u.
Second-order many-body perturbation theory (MP2) can be run with the MP2 command
As with the Hartree-Fock and the DFT method the molecule can be treated with the +1 charge
From MP2 the energy to remove an electron is 0.552 atomic units of energy
MP2 in the aug-cc-pVDZ basis set predicts a lower ionization energy for hydrogen fluoride than DFT with the B3LYP functional.
The parametric 2-RDM method can be run with the Parametric2RDM command
From parametric 2-RDM the energy to remove an electron is 0.586 atomic units of energy
Parametric 2-RDM in the aug-cc-pVDZ basis set predicts a slightly lower ionization energy for hydrogen fluoride than DFT with the B3LYP functional.
The anti-Hermitian contracted Schrödinger equation (ACSE) method can be run with the ContractedSchrodinger command
data1≔ContractedSchrodingerhf,basis=aug-cc-pvdz, active=6,4, active_cse=10,10:
As with Hartree-Fock, DFT, and the Parametric2RDM method the molecule can be treated with the +1 charge
data2≔ContractedSchrodingerhf,basis=aug-cc-pvdz,charge=1,spin=1, active=5,4, active_cse=9,10:
From the ACSE the energy to remove an electron is 0.539 atomic units of energy
The ACSE ionization energy is higher than that from the Hartree-Fock method but lower than that from the parametric 2-RDM method.
The coupled cluster method with single and double excitations can be run with the CoupledCluster command
From coupled cluster the energy to remove an electron is 0.585 atomic units of energy
The coupled cluster singles-doubles result is close to that from the parametric 2-RDM method.
The excited spectra of hydrogen fluoride can be computed from time-dependent density functional theory (TDDFT) with the ExcitationSpectra command
ex_table ≔ ExcitationSpectrahf, method=DensityFunctional, basis=aug-cc-pvdz,showtable=true: