Variational 2-RDM Method
Copyright (c) RDMCHEM LLC 2020
Overview
V2RDM
References
Because electrons are indistinguishable with pairwise interactions, the energy of any atom or molecule can be expressed as a linear functional of the two-electron reduced density matrix (2-RDM). Not every two-electron density matrix, however, corresponds to an N-electron wave function or density matrix. Consequently, we need to constrain the 2-RDM to represent at least one N-electron density matrix through constraints known as N-representability conditions. An approximate set of N-representability conditions is known as the 2-positivity (or DQG) conditions. The 2-positivity conditions correspond to constraining three matrix forms of the 2-RDM, corresponding to the probability distributions of two particles (D), two holes (Q), and a particle-hole pair (G), to be positive semidefinite. (A matrix is positive semidefinite if an only if its eigenvalues are nonnegative). Formally, we can represent the positive semidefiniteness of the D, Q, and G matrices as
2D ≥ 0
2Q ≥ 0
2G ≥ 0 .
Minimization of the ground-state energy with respect to these conditions generates a special type of optimization problem known as a semidefinite program. In the variational 2-RDM method the semidefinite program with N-representability conditions on the 2-RDM is directly solved for the ground-state energy and 2-RDM without computing the many-electron wave function.
Figure 1: 1-electron density of diphenyl disulfide
After loading the Quantum Chemistry package, we explore the ground-state energy and orbitals, dipole moment, and ionization energy of diphenyl disulphide with V2RDM.
Quantum Chemistry
We set the number of Digits to be used in computations to 15 and load the Quantum Chemistry package using Maple's with command.
Digits ≔ 15;
Digits≔15
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
Energy
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 V2RDM with the Variational2RDM command. We select an active space of 6 electrons in 6 orbitals (the active space is the set of electrons and orbitals that will be treated beyond the Hartree-Fock (mean-field) approximation) in the sto-6g basis set
data ≔ Variational2RDMmol, basis=sto-6g, active=6,6;
The ground-state energy is given by
datae_tot;
−1251.46515496
(a) Compute the energy using the [10,10] active space, and compare your result with the energy from the [6,6] active space.
The V2RDM captures electron correlation of the active orbitals which can be observed from the molecular orbital occupations
datamo_occ55..60;
1.999789251.998647551.992309950.004889680.002403100.00196058
The highest occupied molecular orbital (HOMO) can be visualized with the DensityPlot3D command
DensityPlot3Dmol,data,basis=sto-6g, orbitalindex=57; # HOMO
DensityPlot3Dmol,data,basis=sto-6g, orbitalindex=58; # LUMO
(b) Compare these molecular orbital occupations and densities with those from a [10,10] active space.
Dipole Moment
We have the dipole moment from our previous calculation
datadipole;
0.082043001.119573830.06577923
We can visualize the direction of the dipole moment with the DipolePlot command (click on molecule to rotate the perspective)
DipolePlotmol, 'Variational2RDM', basis=sto-6g,active=6,6;
(c) Compute the dipole and its plot with the [10,10] active space, and compare the results with those from the [6,6] active space.
Ionization Energy
We can compute the ionization energy by computing the energy of the ionized molecule
data2 ≔ Variational2RDMmol, basis=sto-6g, active=5,6, charge=1;
Therefore, the ionization energy is given by
IE_hartree ≔ data2e_tot−datae_tot⋅UnitsUnithartree;
IE_hartree≔0.19229399E0
which we can convert from hartrees to eV
IE_eV ≔ convertIE_hartree,units,eV;
IE_eV≔5.23258600eV
(d) Compare this ionization energy (in eV) with the ionization energy (in eV) predicted using the [10,10] active space.
D. A. Mazziotti, Phys. Rev. Lett. 106, 083001 (2011). "Large-scale semidefinite programming for many-electron quantum mechanics"
D. A. Mazziotti, Phys. Rev. Lett. 93, 213001 (2004). "Realization of quantum chemistry without wave functions through first-order semidefinite programming"
A. W. Schlimgen, C. W. Heaps, and D. A. Mazziotti, J. Phys. Chem. Lett. 7, 627-631 (2016). "Entangled electrons foil synthesis of elusive low-valent vanadium oxo complex"
J. M. Montgomery and D. A. Mazziotti, J. Phys. Chem. A 122, 4988-4996 (2018). "Strong electron correlation in nitrogenase cofactor, FeMoco"
Download Help Document
