Dispersion corrected total energy of supercell is not scaling correctly for Cu surface (VASP 5.2.12)

[dsawei]

Hello,
I am not sure if anyone has encounter this before, but I ran into a problem using the vdW-DF optB86b functional. Mainly, when I double the supercell in x and y direction, the total energy of the cell is not scaling with the increase in size.

Simulation:
5 layers of Cu atoms with top 3 layers relaxed to simulate Cu(100) surface. This is the 2 different POSCAR
----------------------------------------------------------------------
Cu100 Surface 1x1
3.59822
1.000 0.000 0.000
0.000 1.000 0.000
0.000 0.000 7.500
Cu
10
Selective dynamics
Direc
0.00 0.00 0.0000000000000000 F F F
0.50 0.50 0.0000000000000000 F F F
0.00 0.50 0.0666666666669968 F F F
0.50 0.00 0.0666666666669968 F F F
0.00 0.00 0.1324608411386063 T T T
0.50 0.50 0.1324608411386063 T T T
0.00 0.50 0.1984807308455974 T T T
0.50 0.00 0.1984807308455974 T T T
0.00 0.00 0.2626299848421017 T T T
0.50 0.50 0.2626299848421017 T T T
----------------------------------------------------------------------
Cu100 Surface 2x2
1.000
7.196 0.000 0.000
0.000 7.196 0.000
0.000 0.000 26.987
Cu
40
Selective dynamics
Direct
0.00 0.00 0.2663985264 T T T
0.00 0.50 0.2663985264 T T T
0.25 0.25 0.2663985264 T T T
0.25 0.75 0.2663985264 T T T
0.50 0.00 0.2663985264 T T T
0.50 0.50 0.2663985264 T T T
0.75 0.25 0.2663985264 T T T
0.75 0.75 0.2663985264 T T T
0.00 0.25 0.2004380277 T T T
0.00 0.75 0.2004380277 T T T
0.25 0.00 0.2004380277 T T T
0.25 0.50 0.2004380277 T T T
0.50 0.25 0.2004380277 T T T
0.50 0.75 0.2004380277 T T T
0.75 0.00 0.2004380277 T T T
0.75 0.50 0.2004380277 T T T
0.00 0.00 0.1335271784 T T T
0.00 0.50 0.1335271784 T T T
0.25 0.25 0.1335271784 T T T
0.25 0.75 0.1335271784 T T T
0.50 0.00 0.1335271784 T T T
0.50 0.50 0.1335271784 T T T
0.75 0.25 0.1335271784 T T T
0.75 0.75 0.1335271784 T T T
0.00 0.25 0.0666666667 F F F
0.00 0.75 0.0666666667 F F F
0.25 0.00 0.0666666667 F F F
0.25 0.50 0.0666666667 F F F
0.50 0.25 0.0666666667 F F F
0.50 0.75 0.0666666667 F F F
0.75 0.00 0.0666666667 F F F
0.75 0.50 0.0666666667 F F F
0.00 0.00 0.0000000000 F F F
0.00 0.50 0.0000000000 F F F
0.25 0.25 0.0000000000 F F F
0.25 0.75 0.0000000000 F F F
0.50 0.00 0.0000000000 F F F
0.50 0.50 0.0000000000 F F F
0.75 0.25 0.0000000000 F F F
0.75 0.75 0.0000000000 F F F
----------------------------------------------------------------------

Basically, the 2x2 POSCAR is a 1x1 tesselated twice in x and y direction. In theory, the 2x2 energy should be 4x that of the 1x1 energy, but this is what the OUTCAR told be about the totally relaxed structure.
----------------------------------------------------------------------
For Cu100 Surface 1x1
Free energy of the ion-electron system (eV)
---------------------------------------------------
alpha Z PSCENC = 365.66130814
Ewald energy TEWEN = 75403.50155529
-1/2 Hartree DENC = -88188.13286699
-exchange EXHF = 0.00000000
-V(xc)+E(xc) XCENC = 427.45615609
PAW double counting = 10982.36246392 -12421.07782025
entropy T*S EENTRO = -0.01203361
eigenvalues EBANDS = -524.93104225
atomic energy EATOM = 13942.92420560
---------------------------------------------------
free energy TOTEN = -12.24807407 eV

energy without entropy = -12.23604046 energy(sigma->0) = -12.24406287
----------------------------------------------------------------------------------------------
For Cu100 Surface 2x2
Free energy of the ion-electron system (eV)
---------------------------------------------------
alpha Z PSCENC = 1462.64252260
Ewald energy TEWEN = 288329.27520804
-1/2 Hartree DENC = -339585.80386113
-exchange EXHF = 0.00000000
-V(xc)+E(xc) XCENC = 1676.83433140
PAW double counting = 43836.61277353 -49601.89921897
entropy T*S EENTRO = 0.02312542
eigenvalues EBANDS = -1901.64506847
atomic energy EATOM = 55771.69682238
---------------------------------------------------
free energy TOTEN = -12.26336519 eV

energy without entropy = -12.28649061 energy(sigma->0) = -12.27107366
----------------------------------------------------------------------
As far as VASP is concerned, there are almost no energy differences between the two structures even thought the 2x2 is 4x the size of 1x1. Can someone tell me what is going on?

Some Support Info.
----------------------------------------------------------------------
KPOINTS 8x8x1 (same for both)

POTCAR PAW_PBE (same for both)

INCAR (same for both)

Startparameter for this Run:
ISTART = 0
LCHARG = .False.
LWAVE = .False.

Electronic Relaxation
ENMAX = 500 eV !specify different value than given in POTCAR
PREC = Normal !Medium=default, Low, High; affects ENMAX, mesh, pspot
EDIFF = 1E-04 !energy stopping-criterion for electr. iterations
EDIFFG = -.03 !force stopping-criterion for geometry steps
LREAL = .FALSE !real-space projection

Ionic Relaxation
NSW = 500 !max number of geometry steps
IBRION = 2 !ionic relax: 0-MD 1-quasi-New (RMM-DIIS for NEB)
ISIF = 2 !(0:force=y stress=n ions=y shape=n volume=n,
ISYM = 1 !1=use symmetry, 0 = no symmetry
POTIM = .5 !initial time step for geo-opt (increase for soft sys)

Dispersion Correction
LUSE_VDW = .TRUE.
AGGAC = 0.0000
GGA = MK
PARAM1 = 0.1234
PARAM2 = 1.0000

DOS related values:
ISMEAR = 1 ! -4-tet -1-fermi 1=Methfessel/Paxton 1.order
SIGMA = .2 !broadening in eV
----------------------------------------------------------------------
<span class='smallblacktext'>[ Edited ]</span>

[kelum]

Hi,
I tried to reproduce this problem but without success, I used your INCAR, POSCAR, the Cu
POTCAR, and only changed to cheaper 4x4x1 k-point mesh and did only 5 ionic
relaxation steps, both changes shouldn't make a significant difference.
For the 1x1 cell I get (from OSZICAR) for your geometry
1 F= -.12108180E+02
which agrees with your number and for the 2x2 cell I get:
1 F= -.48992073E+02
which is close to four times the 1x1 number. That's reasonable given that the k-points
density is not the same for the two calculations, if the k-point meshes correspond to each
other the 2x2 energy is four times the 1x1 energy +-1 meV. So I'm not sure what went
wrong.
Could you please show what are the values of the vdW corrections for the two runs,
to see if there is some problem there or not? I get about 21.5 eV for the small surface
and about 86.4 eV for the 2x2 surface.
I assume you get the right numbers using just PBE?

[dsawei]

Hi,

If you are getting reasonable result, then I am really confused.

For 1x1
Total vdW correction in eV: 21.5894753738756

Fpr 2x2
Total vdW correction in eV: 123.876143821605

There appears to be a significant over correction for the 2x2, but I am not really sure what could be causing that problem.

[dsawei]

Hi,

After doubling checking all my calculations, the PBE and the Grimme's D2 dispersion correction all gets the right number, but the vDW-DF correction seems to have this problem.

Is it possible for you to email me your input files (including VASP version number etc.)? I want to compare it to my own input. I might've missed something that is obvious to other people.

[dsawei]

The problem was from an incorrect (or corrupt/something wrong with it) vdw_kernel.bindat file. After obtaining a new vdw_kernel.bindat file, the energy is scaled correctly.

Just fyi if anyone else ran into similar problem.