How to calculate DFT energy with density from another level of theory?

Question

I need to calculate the DFT energy given the density of another level of theory, for example Hartree Fock.

$E[\rho_{HF}]$

This is sometimes used to fix various electron delocalization errors in conjugated molecules, such as what is done here to fix torsional barriers in isomerization reactions https://arxiv.org/pdf/2102.06842.pdf

This of course means that there is no SCF procedure done.

How can I do this in various codes? I am specifically interested in input examples for Q-Chem, but if you know of a way to do it for other codes that could be useful for others. The linked paper uses ORCA.

Answer 1

So this question is related to the recent discussion in Procedure to classify errors in Kohn-Sham DFT where I already gave the general procedure: get the density from some level of theory, and use it to run a single-point DFT calculation.

This question was about how to do this in Q-Chem; it's actually quite simple. Here's an example input for a calculation on the water molecule, where I first run Hartree-Fock to get the reference density, and then read in the wave function to a calculation with the PBE functional

$molecule
0 1
O        0.000000    0.000000    0.117790 
H        0.000000    0.755453   -0.471161 
H        0.000000   -0.755453   -0.471161 
$end

$rem
basis pcseg-2
method hf
$end

@@@

$molecule
read
$end

$rem
basis pcseg-2
scf_guess read
method pbe
max_scf_cycles 0
$end

Setting the maximum number of SCF cycles to 0 in Q-Chem means that the program will just evaluate a single-point energy without changing the orbitals. The output for the second part of the job is

 The restricted Kohn-Sham energy will be evaluated
 The orbitals will not be altered
 Exchange:  PBE      Correlation:  PBE
 Using SG-1 standard quadrature grid
 Using Q-Chem read-in guess as SCF_GUESS READ specified.

Avoiding writing the new MOs to disk !!! 
 ---------------------------------------
  Cycle       Energy         DIIS Error
 ---------------------------------------
    1     -76.3736280198      6.70E-03 Convergence criterion met
 ---------------------------------------
 SCF time:  CPU 0.17 s  wall 0.17 s
 SCF   energy in the final basis set = -76.37362802
 Total energy in the final basis set = -76.37362802

If you also need the converged PBE energy, removing the line that sets max_scf_cycles to zero leads to the following output

 Exchange:  PBE      Correlation:  PBE
 Using SG-1 standard quadrature grid
 A restricted SCF calculation will be
 performed using DIIS
 SCF converges when DIIS error is below 1.0e-05
 ---------------------------------------
  Cycle       Energy         DIIS error
 ---------------------------------------
    1     -76.3736280198      6.70e-03  
    2     -76.3810410874      4.41e-04  
    3     -76.3810699365      3.92e-04  
    4     -76.3811039575      2.38e-04  
    5     -76.3811190178      2.47e-05  
    6     -76.3811192100      3.84e-07  Convergence criterion met
 ---------------------------------------

Since the DIIS error is quite small in the first iteration, and the convergence is very rapid after that, you see that the HF density is a pretty good guess for the PBE functional; this was one of the findings in my assessment of initial guesses, J. Chem. Theory Comput. 15, 1593 (2019). You also see that optimizing the electron density lowers the energy by 7.5 $mE_h$, but this is how you would also reintroduce the density errors.

Answer 2

Incidentally, I initiated some discussions along these lines in these questions for HF-DFT (also called 'density-corrected' DFT) here, and here. I'm familiar with Quantum ESPRESSO, so I can detail a procedure with QE. If you want to generate a Hartree-Fock density, you can just turn on a hybrid functional (say HSE06) and turn up the fraction of exact exchange to 1. This can be done by setting exx_fraction = 1. This way, the generated Hartree-Fock density albeit in a Kohn-Sham scheme (i.e. KS eigenvalues and orbitals), should coincide with the energy and density of a calculation involving self-consistent hartree-fock orbitals.

You can then use this density as the starting point, and turn on your regular functionals (say PBE, LDA etc).