Reduced Density Matrix Functional Theory (RDMFT) computes the energy of a many-electron atom or molecule as a functional of the one-electron reduced density matrix rather than the many-electron wavefunction ψ(${\mathrm{r}}_1\, {\mathrm{r}}_2$, ..., ${\mathrm{r}}_{\mathrm{N}}$).

The command adds a universal correction, based on the 1-electron reduced density matrix (1-RDM) rather than the density alone, to any density functional theory (DFT) functional.

A well-known limitation of DFT is its difficulty in predicting the energies and properties of molecules with static correlation, which is important to the computation of charges, van der Waals forces, barrier heights, and bi- and multiradicals.

RDMFT's universal correction to DFT allows the orbital occupations to become fractional and thereby, capture strong correlation.

Practically, RDMFT has a computational scaling of O(N³) similar to that of DFT.  In RDMFT the energy is minimized by semidefinite programming.

RDMFunctional employs "B3LYP" as its default exchange-correlation (xc) functional.  Although "B3LYP" is the default xc functional, the RDMFunctional is currently optimized for semi-local (non-hybrid) functionals such as "PBE" and "SCAN".  The full list of available functionals is provided in the help page Functionals.  Furthermore, the available functionals can be searched with the SearchFunctionals command.   

On the Windows operating system the RDMFunctional command requires the installation of Microsoft's Windows Subsystem for Linux (WSL).  For Windows 10 (version 2004 and higher) and Windows 11 you can install the WSL by opening the Command Prompt in administrator mode and entering the command: wsl --install -d Ubuntu  For additional details, please refer to: https://learn.microsoft.com/en-us/windows/wsl/install

basis = string -- name of the basis set.  See Basis for a list of available basis sets.  Default is "sto-3g".

spin = nonnegint -- twice the total spin S (= 2S). Default is 0.

charge = nonnegint -- net charge of the molecule. Default is 0.

unit = "Angstrom" or "Bohr". Default is "Angstrom".

max_memory = posint -- allowed memory in MB.

xc = string -- The full list of available functionals is provided in the help page Functionals.  Default functional is "B3LYP".

solvent = string -- name of the solvent.  See Solvent.  Default is "None".

ghost = list of lists -- each list has the string of an atom's symbol and the atom's x, y, and z coordinates.  See Ghost Atoms.

nuclear_gradient = boolean -- option to return the analytical nuclear gradient if available. Default is false.

active = list -- optional list [N,r] where r is the maximum number of orbitals allowed to be fractional and N is the number of electrons populating those r orbitals.  By default, all orbitals are allowed to become fractional.

populations = string -- atomic-orbital population analysis: "Mulliken" and "Mulliken/meta-Lowdin". Default is "Mulliken".

conv_tol = float -- converge threshold. Default is ${10}^{-10}\.$

diis = boolean -- whether to employ diis. Default is true.

diis_space = posint -- diis's space size. By default, 8 Fock matrices and errors vector are stored.

diis_start_cycle = posint -- the step to start diis. Default is 1.

direct_scf = boolean -- direct SCF is used by default.

direct_scf_tol = float -- direct SCF cutoff threshold. Default is ${10}^{-13}\.$

level_shift_factor = float/int -- level shift (in au) for virtual space. Default is $0.$

max_cycle = posint -- max number of iterations. Default is 50.

max_rho_cutoff = float -- ${10}^{-7}\.$

small_rho_cutoff = float -- $0.$

grids_atomic_radii = string -- "radi.BRAGG_RADII" (default), "radi.COVALENT_RADII", or "None".

grids_becke_scheme = string -- "gen_grid.original_becke" (default) or "gen_grid.stratmann".

grids_level = (0-9) -- big number for large mesh grids. Default is 3.

grids_prune = "gen_grid.nwchem_prune" (default), "gen_grid.sg1_prune", "gen_grid.treutler_prune", or "None".

grids_radi_method = "radi.treutler" (default), "radi.delley", "radi.mura_knowles", or "radi.gauss_chebyshev". 

grids_radii_adjust = "radi.treutler_atomic_radii_adjust" (default), "radi.becke_atomic_radii_adjust", or "None".

grids_symmetry = boolean -- true or false to symmetrize mesh grids.

D. Gibney, J.-N. Boyn, and D. A. Mazziotti, Phys. Rev. Lett. 131, 243003 (2023). "Universal Generalization of Density Functional Theory for Static Correlation"

D. Gibney, J.-N. Boyn, and D. A. Mazziotti, J. Phys. Chem. Lett. 13, 1382–1388 (2022). "Density Functional Theory Transformed into a One-Electron Reduced-Density-Matrix Functional Theory for the Capture of Static Correlation"

D. Gibney, J.-N. Boyn, and D. A. Mazziotti, J. Chem. Theory Comput. 18, 6600–6607 (2022). "Comparison of Density-Matrix Corrections to Density Functional Theory"

D. Gibney, J.-N. Boyn, and D. A. Mazziotti, J. Phys. Chem. Lett. 12, 385–391 (2021). "Toward a Resolution of the Static Correlation Problem in Density Functional Theory from Semidefinite Programming"

R. G. Parr and W. Yang, Density Functional Theory of Atoms and Molecules (Oxford University Press, New York, 1989).

$\mathrm{with}\left(\mathrm{QuantumChemistry}\right)\;$

$\left[\mathrm{AOLabels}\,\mathrm{ActiveSpaceCI}\,\mathrm{ActiveSpaceSCF}\,\mathrm{AtomicData}\,\mathrm{BondAngles}\,\mathrm{BondDistances}\,\mathrm{Charges}\,\mathrm{ChargesPlot}\,\mathrm{Chat}\,\mathrm{ContractedSchrodinger}\,\mathrm{CorrelationEnergy}\,\mathrm{CoupledCluster}\,\mathrm{DensityFunctional}\,\mathrm{DensityPlot3D}\,\mathrm{Dipole}\,\mathrm{DipolePlot}\,\mathrm{Energy}\,\mathrm{ExcitationEnergies}\,\mathrm{ExcitationSpectra}\,\mathrm{ExcitationSpectraPlot}\,\mathrm{ExcitedStateEnergies}\,\mathrm{ExcitedStateSpins}\,\mathrm{ExcitonDensityPlot}\,\mathrm{ExcitonPopulations}\,\mathrm{ExcitonPopulationsPlot}\,\mathrm{FullCI}\,\mathrm{GeometryOptimization}\,\mathrm{HartreeFock}\,\mathrm{Interactive}\,\mathrm{Isotopes}\,\mathrm{LiteratureSearch}\,\mathrm{MOCoefficients}\,\mathrm{MODiagram}\,\mathrm{MOEnergies}\,\mathrm{MOIntegrals}\,\mathrm{MOOccupations}\,\mathrm{MOOccupationsPlot}\,\mathrm{MOSymmetries}\,\mathrm{MP2}\,\mathrm{MolecularData}\,\mathrm{MolecularDictionary}\,\mathrm{MolecularGeometry}\,\mathrm{NuclearEnergy}\,\mathrm{NuclearGradient}\,\mathrm{OscillatorStrengths}\,\mathrm{Parametric2RDM}\,\mathrm{PlotMolecule}\,\mathrm{Populations}\,\mathrm{Purify2RDM}\,\mathrm{QuantumComputing}\,\mathrm{RDM1}\,\mathrm{RDM2}\,\mathrm{RDMFunctional}\,\mathrm{RTM1}\,\mathrm{ReadXYZ}\,\mathrm{Restore}\,\mathrm{Save}\,\mathrm{SaveXYZ}\,\mathrm{SearchBasisSets}\,\mathrm{SearchFunctionals}\,\mathrm{SkeletalStructure}\,\mathrm{SolventDatabase}\,\mathrm{Thermodynamics}\,\mathrm{TransitionDipolePlot}\,\mathrm{TransitionDipoles}\,\mathrm{TransitionOrbitalPlot}\,\mathrm{TransitionOrbitals}\,\mathrm{Variational2RDM}\,\mathrm{VibrationalModeAnimation}\,\mathrm{VibrationalModes}\,\mathrm{Video}\right]$

$\mathrm{bond\_distances} ≔ \mathrm{map}\left(x \to x\+2.0\, \left[\mathrm{fsolve}\left(\mathrm{expand}\left(\mathrm{ChebyshevT}\left(8\,x\right)\right)\right)\right]\right)\;$

$\mathrm{bond\_distances}≔\left[1.01921472\,1.16853039\,1.44442977\,1.80490968\,2.19509032\,2.55557023\,2.83146961\,2.98078528\right]$

$\mathrm{molecules} ≔ \left[\mathrm{seq}\left(\left[\left[''N''\,0\,0\,0\right]\,\left[''N''\,0\,0\,R\right]\right]\, R \mathbf{in} \mathrm{bond\_distances}\right)\right]\;$

$\mathrm{molecules}≔\left[\left[\left[''N''\,0\,0\,0\right]\,\left[''N''\,0\,0\,1.01921472\right]\right]\,\left[\left[''N''\,0\,0\,0\right]\,\left[''N''\,0\,0\,1.16853039\right]\right]\,\left[\left[''N''\,0\,0\,0\right]\,\left[''N''\,0\,0\,1.44442977\right]\right]\,\left[\left[''N''\,0\,0\,0\right]\,\left[''N''\,0\,0\,1.80490968\right]\right]\,\left[\left[''N''\,0\,0\,0\right]\,\left[''N''\,0\,0\,2.19509032\right]\right]\,\left[\left[''N''\,0\,0\,0\right]\,\left[''N''\,0\,0\,2.55557023\right]\right]\,\left[\left[''N''\,0\,0\,0\right]\,\left[''N''\,0\,0\,2.83146961\right]\right]\,\left[\left[''N''\,0\,0\,0\right]\,\left[''N''\,0\,0\,2.98078528\right]\right]\right]$

$\mathrm{energies\_dft}≔ \left[\mathrm{seq}\left(\mathrm{Energy}\left(\mathrm{molecule}\,\mathrm{method}\=\mathrm{DensityFunctional}\, \mathrm{basis}\=''cc-pVDZ''\, \mathrm{xc}\=''B3LYP''\right)\, \mathrm{molecule} \mathbf{in} \mathrm{molecules}\right)\right]\;$

$\mathrm{energies\_dft}≔\left[\mathrm{-109.50786339}\,\mathrm{-109.52331040}\,\mathrm{-109.37588586}\,\mathrm{-109.17741963}\,\mathrm{-109.03358249}\,\mathrm{-108.95146495}\,\mathrm{-108.90997138}\,\mathrm{-108.89291096}\right]$

$\mathrm{pes\_dft} ≔ \mathrm{interp}\left(\mathrm{bond\_distances}\,\mathrm{energies\_dft}\,R\right)\;$

$\mathrm{pes\_dft}≔-0.14339035R^7+2.20942302R^6-14.44738600R^5+51.96535307R^4-110.97891955R^3+140.38786698R^2-96.50876625R-81.97856683$

$\mathrm{energies\_rdm}≔ \left[\mathrm{seq}\left(\mathrm{Energy}\left(\mathrm{molecule}\,\mathrm{method}\=\mathrm{RDMFunctional}\, \mathrm{basis}\=''cc-pVDZ''\, \mathrm{xc}\=''B3LYP''\right)\, \mathrm{molecule} \mathbf{in} \mathrm{molecules}\right)\right]\;$

$\mathrm{energies\_rdm}≔\left[\mathrm{-109.50786339}\,\mathrm{-109.52331040}\,\mathrm{-109.37588586}\,\mathrm{-109.18887170}\,\mathrm{-109.10062254}\,\mathrm{-109.06672079}\,\mathrm{-109.05806429}\,\mathrm{-109.05565802}\right]$

$\mathrm{pes\_rdm} ≔ \mathrm{interp}\left(\mathrm{bond\_distances}\,\mathrm{energies\_rdm}\,R\right)\;$

$\mathrm{pes\_rdm}≔0.00403016R^7+0.13029220R^6-2.19462464R^5+12.96198119R^4-38.72108443R^3+62.54064536R^2-51.32917602R-92.88749381$

$$\begin{gathered} \mathrm{p\_dft}≔ \mathrm{plot}\left(\mathrm{pes\_dft}\, R\=1.0..2.9\, \mathrm{axes}\=\mathrm{boxed}\, \mathrm{labels}\=\left[''Bond\; Distance\; (\AA )''\,"Energy\left(\mathrm{hartree}\right)\"\right]\, \mathrm{color}\=\mathrm{red}\, \mathrm{thickness}\=3\, \mathrm{labeldirections}\=\left[\mathrm{horizontal}\,\mathrm{vertical}\right]\right)\: \\ \mathrm{p\_rdm} ≔ \mathrm{plot}\left(\mathrm{pes\_rdm}\, R\=1.0..2.9\, \mathrm{axes}\=\mathrm{boxed}\, \mathrm{labels}\=\left[''Bond\; Distance\; (\AA )''\,"Energy\left(\mathrm{hartree}\right)\"\right]\, \mathrm{color}\=\mathrm{blue}\, \mathrm{thickness}\=3\right)\: \\ \mathrm{plots}:-\mathrm{display}\left(\left\{\mathrm{p\_dft}\,\mathrm{p\_rdm}\right\}\right)\; \end{gathered}$$

$\mathrm{molecule} ≔ \left[\left[''H''\,0\,0\,0\right]\,\left[''F''\,0\,0\,2.0\right]\right]\;$

$\mathrm{molecule}≔\left[\left[''H''\,0\,0\,0\right]\,\left[''F''\,0\,0\,2.00000000\right]\right]$

$\mathrm{data} ≔ \mathrm{RDMFunctional}\left(\mathrm{molecule}\right)\;$

$\mathrm{data}≔table\left(\left[\mathrm{charges}=\left[\right]\,\mathrm{mo\_coeff}=\left[\right]\,\mathrm{mo\_symmetry}=\left[\right]\,\mathrm{populations}=\left[\right]\,\mathrm{mo\_occ}=\left[\right]\,\mathrm{converged}=1\,\mathrm{group}=''C1''\,\mathrm{aolabels}=\left[\right]\,\mathrm{mo\_energy}=\left[\right]\,\mathrm{dipole}=\left[\right]\,\mathrm{e\_tot}=\mathrm{-98.77669656}\,\mathrm{rdm1}=\left[\right]\right]\right)$