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

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.

• No time right now but check out Eunji Sim's answer to the milestones question! – Nike Dattani Mar 5 at 1:36

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

• Ahh brilliant, I thought It would be like two calculations one with SCF and the other without SCF. But this looks like what I'm looking for. – Cody Aldaz Mar 5 at 7:52
• I'm actually curious to receive more answers to your question. Because the method I described is a bit expensive - The computation of the Fock operator on the wavefunction is atleast one order of magnitude more expensive than PBE or LDA calculation. Maybe there is a way to generate Hartree-Fock density without going through this whole rigmarole. – livars98 Mar 5 at 7:59
• I think the procedure described in this answer is problematic due to several reasons. 1. The resulting energy depends on when in the SCF iteration the ingredients for the total energy are collected. Is QE using the input density or the output density for this? Do they use the Kohn-Sham eigenvalues? Maybe not all ingredients are from the starting point you want. 2. HSE has a screening built in. That may give you a different starting point than you think. Probably QE can do HF calculations directly. – Gregor Michalicek Mar 5 at 9:02
• @GregorMichalicek The corresponding energy from the SCF calculation should be taken from the first iteration, I believe. I'm not sure if there are other ways of doing HF calculations on QE. I found this discussion on a QE forum: mail-archive.com/users@lists.quantum-espresso.org/msg28524.html. – Xivi76 Mar 5 at 9:21
• Of course, you take it from the first iteration. I don't doubt that. But there are also different places within an iteration. If you ask me: If QE doesn't write that out by default take a look at the code. It is open source. Maybe you can modify the code in a simple way to get what you want in a consistent way. – Gregor Michalicek Mar 5 at 9:43

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

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