Tutorial - Maple Help

Online Help

All Products    Maple    MapleSim


Quantum Chemistry Package Tutorial Worksheet

Copyright (c) RDMCHEM LLC 2021

 

Overview

Load the Package

Definition of a Molecule

Computations

Visualization

Interactive Maplet

Courses

Help

Overview

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'sNew2021 help page for examples of what's new in the QuantumChemistry Package 2021.  

Load the Package

We begin by loading the commands of the Quantum Chemistry package using the Maple with command

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

(2.1)

We recommend setting the global variable Digits to 15 (double precision) rather than 10 (default):

Digits15;

Digits15

(2.2)

 

Definition of a Molecule

A molecule's geometry can be defined in the worksheet interface in three ways:

1. 

entering the geometry as a Maple list of lists

2. 

reading the geometry from an XYZ file

3. 

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:

1. 

a string of the atomic symbol of the atom, i.e. "H" or "F"

2. 

the x Cartesian coordinate of the atom

3. 

the y Cartesian coordinate of the atom

4. 

the z Cartesian coordinate of the atom

The default unit of distance in the package is Angstroms (1 Angstrom= 1010 m), which we employ in this worksheet.  

hfH,0,0,0.48,F,0,0,0.48;

hfH,0,0,−0.48000000,F,0,0,0.48000000

(3.1.1)

 

Web Database

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

hfMolecularGeometryhydrogen fluoride;

hfF,0,0,0,H,0.43580000,−0.14770000,−0.81880000

(3.2.1)

 

XYZ File

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:

2

# 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;

porphinN,−1.62058500,−1.35841100,−0.00004500,N,1.62042400,−1.35852600,−0.00009500,N,−1.57898000,1.33345500,0.00027900,N,1.57912800,1.33329300,−0.00026600,C,−1.27247000,−2.69668000,0.00024100,C,1.27217700,−2.69676300,0.00031100,C,−3.00281900,−1.21613100,−0.00051500,C,3.00267300,−1.21642300,−0.00028500,C,−0.00016600,−3.27300100,0.00055400,C,−2.51334300,−3.43612200,0.00013200,C,2.51294900,−3.43634000,0.00024600,C,−3.55872200,−2.54069000,−0.00034400,C,3.55842900,−2.54101500,0.00016400,C,−3.63768700,0.02511200,−0.00061100,C,3.63773300,0.02474200,−0.00044500,C,−2.95501700,1.24875700,−0.00024700,C,2.95521500,1.24846000,−0.00035900,C,−1.27258700,2.67294600,0.00036400,C,1.27286200,2.67274100,0.00017600,C,0.00017500,3.25342200,0.00036700,C,−3.55030400,2.58382300,−0.00002800,C,3.55059000,2.58346400,−0.00045100,C,−2.50983200,3.46576200,0.00046900,C,2.51016700,3.46547700,0.00019800,H,−0.00017100,−4.35835700,0.00090700,H,−1.06231600,−0.49153600,0.00028100,H,1.06235700,−0.49148500,0.00004000,H,−2.57314000,−4.51439200,0.00038500,H,2.57252600,−4.51461900,0.00038400,H,−4.61490100,−2.76464200,−0.00062200,H,4.61460100,−2.76500800,0.00017400,H,−4.72223800,0.03505700,−0.00090700,H,4.72228200,0.03449200,−0.00050400,H,0.00026800,4.34141600,0.00063100,H,−4.60947300,2.79979000,−0.00015400,H,4.60975700,2.79940900,−0.00066000,H,−2.55314700,4.54611800,0.00079700,H,2.55354400,4.54583900,0.00050700

(3.3.1)

Computations

Computations of molecules can be performed in one of two ways:

1. 

Maple commands based on the energy or property to be computed  

2. 

Maple commands based on the quantum chemistry method to be employed

Energies/Properties

The QuantumChemistry package contains commands to compute directly the following energies and properties:

• 

energy

• 

correlation energy

• 

nuclear energy

• 

molecular orbital coefficients

• 

molecular orbital occupations

• 

molecular orbital energies

• 

orbital symmetries

• 

atomic-orbital populations and charges

• 

dipole moments

• 

nuclear gradients

• 

excited-state energies

• 

excited-state spins

• 

excitation spectra

• 

transition dipole moments

• 

oscillator strengths

• 

natural transition orbitals

• 

thermodynamic properties

• 

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

Energyhf,method=HartreeFock;

−98.57249547

(4.1.1)

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   

Energyhf,method=HartreeFock,basis=cc-pvdz;

−100.01804664

(4.1.2)

We can select a different method for computing the energy

Energyhf,method=DensityFunctional,basis=cc-pvdz;

−100.39764610

(4.1.3)

The nuclear energy from the repulsion of the nuclei in the molecule can be computed with the command NuclearEnergy

NuclearEnergyhf;

5.07069443

(4.1.4)

The energies of the molecular orbitals can be computed with the command MOEnergies

MOEnergieshf,method=HartreeFock;

25.902025751.464436470.5785509840.4635098260.4635098260.606112018

(4.1.5)

The occupations of the molecular orbitals can be computed with the command MOOccupations

MOOccupationshf,method=HartreeFock;

2.000000002.000000002.000000002.000000002.000000000.000000000

(4.1.6)

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

MOOccupationshf,method=Variational2RDM, active=6,4;

2.000000002.000000001.999252851.999252841.974011660.0274825764

(4.1.7)

The 1-electron reduced density matrix (1-RDM) can be calculated from the variational 2-RDM method with the command RDM1

dm1RDM1hf,method=Variational2RDM,active=6,4;

1.973724200.0000000000.0000000000.02365336140.0000000001.999252840.0000000000.0000000000.0000000000.0000000001.999252850.0000000000.02365336140.0000000000.0000000000.0277700441

(4.1.8)

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,

dm1Matrixdm1,shape=symmetric;

1.973724200.0000000000.0000000000.02365336140.0000000001.999252840.0000000000.0000000000.0000000000.0000000001.999252850.0000000000.02365336140.0000000000.0000000000.0277700441

(4.1.9)

and then we compute the eigenvalues

eigsLinearAlgebra:-Eigenvaluesdm1;

0.02748257641.974011661.999252841.99925285

(4.1.10)

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

rdataseqr,r=0.8..2.0,0.05:

moleculesseqH,0,0,0,F,0,0,i,i=0.8..2.0,0.05:

edata_hfmapEnergy,molecules,method=HartreeFock:

pes_hfsplinerdata,edata_hf,R,1:

plot_hfplotpes_hf,R=0.8..2.0,axes=BOXED,font=ROMAN,16,labelfont=ROMAN,16,labels=R,E,color=blue;

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

edata_cimapEnergy,molecules,method=ActiveSpaceCI,active=3,4,5,6:

pes_hfseqrdatai,edata_cii,i=1..nopsrdata:

plot_ciplots:-pointplotpes_hf,axes=BOXED,font=ROMAN,16,labelfont=ROMAN,16,symbolsize=16,labels=R,E,color=black;

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

edata_2rdmmapEnergy,molecules,method=Variational2RDM,symmetry=true,conditions=DQ,active=3,4,5,6:

pes_2rdmsplinerdata,edata_2rdm,R,1:

plot_2rdmplotpes_2rdm,R=0.8..2.0,axes=BOXED,font=ROMAN,16,labelfont=ROMAN,16,labels=R,E,color=red;

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

plots:-displayplot_ci,plot_hf,plot_2rdm;

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

hf_hf,data_hfGeometryOptimizationhf,HartreeFock: hf_hf;

F,−0.00376378,0.00127561,0.00707155,H,0.43956378,−0.14897561,−0.82587155

(4.1.11)

Because the fluorine is located at the origin, the bond distance can be computed as follows

bonddistance_hf  BondDistanceshf_hf,1,2;

bonddistance_hf0.95546270

(4.1.12)

Once the optimized geometry is known, the vibrational frequency of the bond can be determined to be 4459 cm1 from

emode,vmodeVibrationalModeshf_hf,HartreeFock;

emode,vmode4458.50363636,0.10414645−0.03529699−0.19567489−0.452153250.153242390.84952527

(4.1.13)

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,

hf_2rdm,data_2rdmGeometryOptimizationhf,Variational2RDM,conditions=DQ,active=3,4,5,6:
hf_2rdm;

F,−0.01249964,0.00423634,0.02348487,H,0.44829964,−0.15193634,−0.84228487

(4.1.14)

In this case the equilibrium bond distance, the distance at which the minimum energy occurs, is

bonddistance_2rdmBondDistanceshf_2rdm,1,2;

bonddistance_2rdm0.99311789

(4.1.15)

The vibrational frequency of the bond can be determined in cm1 from

emode,vmodeVibrationalModeshf_2rdm,Variational2RDM,conditions=DQ,active=3,4,5,6;

emode,vmode3780.92884067,0.10408256−0.03519865−0.19572664−0.451876320.152815440.84974949

(4.1.16)

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

Methods

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

dataHartreeFockhf, basis=dz;

table%id=18446744143386586334

(4.2.1)

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:

datae_tot

−100.02155522

(4.2.2)

Similarly, the energies of the molecular orbitals can be retrieved as follows

datamo_energy

_rtable18446744078522350878

(4.2.3)

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

dataVariational2RDMhf,basis=cc-pvdz, active=3,6,conditions=DQ:

 

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

data1HartreeFockhf,basis=aug-cc-pvdz:

The charge of the molecule can be changed using the charge and spin keywords

data2HartreeFockhf,basis=aug-cc-pvdz,charge=1,spin=1:

From the two energies we can estimate the ionization energy, the energy to remove an electron as 0.531 a.u.

iedata2e_totdata1e_tot

ie0.53141295

(4.2.4)

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

data1DensityFunctionalhf,basis=aug-cc-pvdz:

As with the Hartree-Fock method the charge of the molecule can be changed using the charge and spin keywords

data2DensityFunctionalhf,basis=aug-cc-pvdz,charge=1,spin=1:

From DFT with the default B3LYP functional the energy to remove an electron is 0.596 a.u.

iedata2e_totdata1e_tot

ie0.59619956

(4.2.5)

 

Second-order many-body perturbation theory (MP2) can be run with the MP2 command

data1MP2hf,basis=aug-cc-pvdz:

As with the Hartree-Fock and the DFT method the molecule can be treated with the +1 charge

data2MP2hf,basis=aug-cc-pvdz,charge=1,spin=1:

From MP2 the energy to remove an electron is 0.552 atomic units of energy

iedata2e_totdata1e_tot

ie0.55160015

(4.2.6)

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

data1Parametric2RDMhf,basis=aug-cc-pvdz:

As with the Hartree-Fock and the DFT method the molecule can be treated with the +1 charge

data2Parametric2RDMhf,basis=aug-cc-pvdz,charge=1,spin=1:

From parametric 2-RDM the energy to remove an electron is 0.586 atomic units of energy

iedata2e_totdata1e_tot

ie0.58562861

(4.2.7)

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

data1ContractedSchrodingerhf,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

data2ContractedSchrodingerhf,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

iedata2e_totdata1e_tot

ie0.53884898

(4.2.8)

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

data1CoupledClusterhf,basis=aug-cc-pvdz:

As with the Hartree-Fock and the DFT method the molecule can be treated with the +1 charge

data2CoupledClusterhf,basis=aug-cc-pvdz,charge=1,spin=1:

From coupled cluster the energy to remove an electron is 0.585 atomic units of energy

iedata2e_totdata1e_tot

ie0.58533448

(4.2.9)

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:

State

Energy

Wavelength

Spin

Oscillator

1

8.83272265eV

140.36917303nm

Triplet

0.03753967

2

8.83272575eV

140.36912368nm

Triplet

0.03753965

3

9.25846857eV

133.91436875nm

Singlet

0.03516694

4

9.25847182eV

133.91432172nm

Singlet

0.03516693

5