My Community

Hello everyone,

I have just noticed that EXX energies for isolated molecules, i.e., energies from a calculation with the following INCAR, starting from the WAVECAR of a PBE run on the same geometry:

5 GGA = PE
6
7 LHFCALC= .TRUE.
8 AEXX = 1.0
9 ALDAC = 0.0
10 AGGAC = 0.0
11
12 Electronic Relaxation
13
14 ENCUT = 520
15 EDIFF = 1E-08
16 ALGO = EIGENVAL
17 NELM = 1
18
19
20 Ionic relaxation
21
22 EDIFFG = -0.01
23 IBRION = -1
24 NSW = 0
25
26 Box relaxation
27
28 ISIF = 2
29
30 DOS related values
31
32 ISMEAR = 0
33 SIGMA = 0.05

are heavily dependent on the chosen box. This is, increasing the vacuum layer (beyond a vacuum level where no effects of mirror images are observed in the PBE run) does change the obtained EXX energy. The obtained changes are quite large (up to 0.06 eV for a H20 molecule), see data below:

l / A E0(DFT) / eV E0(EXX) / eV
10.000 -14.265 -29.212
15.000 -14.264 -29.171
20.000 -14.262 -29.109
25.000 -14.261 -29.102
30.000 -14.261 -29.099
35.000 -14.261 -29.097
40.000 -14.261 -29.096

(The chosen box is cubic with side length l). This shows that the PBE energy is converged to 0.01 eV for 15 A, but the EXX energy would require a huge box of 25 A for the same accuracy. It get's even worse when calculating the binding energies of H20 dimers. The same effect can be seen by employing more k-points then the Gamma point. For a box with l=10 A, I obtain the following values by employing an increasing amount of k-points (MP-grid: n x n x n):

#kpt (n) E0(DFT) E0(EXX)
1.000 -14.265 -29.212
2.000 -14.265 -29.114
3.000 -14.265 -29.104
4.000 -14.265 -29.102
5.000 -14.265 -29.101
6.000 -14.265 -29.100

I think this problem is linked to the handling of the discontinuity in the HF exchange, i.e., this dependence on k-points/box size is an error, which becomes smaller by increasing either parameter. In the case of k-points one could argue that this is the case because the problematic contributions <phi_nk phi_nk|phi_nk phi_nk> become smaller with increasing k-point grid.

So this raises two questions:

1) Is this a bug or known behavior?
2) How to deal with it, i.e., increasing the k-point grid or the box size (I think the latter is computationally more demanding)?

The increase of the box volume leads to the increase of the number of plane waves.
Although the molecules remains the same, the size of the basis set increases.
Therefore the calculated energies (both E0 and EXX) differ.
The best way how to deal with it is to use the same box for all calculations.

While this is in principle true, don't you think that this effect should play a minor role for l > 10-15 A (this certainly shows in the DFT energies) and that the effect of increasing the k-mesh hints to another problem?