Inconsistent energy correction between executable, Python API, and VASP

[Andrew-S-Rosen]

I am getting some inconsistent energy corrections depending on how I fetch them. @tgmaxson, this looks far less notable than in #68 but tagging in case you have some wisdom to share.

With s-dftd3 1.2.1:

s-dftd3 CONTCAR --bj pbe

returns -2.0077614201867E-01 Eh or -5.463397149451446 eV (-5.463397104953209 eV with CODATA 2014, the ASE default).

from ase.io import read
from dftd3.ase import DFTD3

atoms = read("CONTCAR")
atoms.calc = DFTD3(method="PBE", damping="d3bj")
e = atoms.get_potential_energy()
print(e)

returns -5.456782359427048. This difference of ~0.0066 eV is not resolved by setting alp=14.0 in light of #150.

In VASP 6.5.1, it is using DFTD3 V3.0 Rev 1 and returning -5.46486 eV, which is closer to s-dtd3 but still different by 0.0015 eV.

CONTCAR.txt

OUTCAR.txt

[tgmaxson]

If I recall correctly with our original testing, VASP and this package may give slightly different values since the PBE parameter values are hard coded into VASP. It might be good to somehow verify they both use the same parameterization beyond just the simple PBE-D3(BJ).

For MLIPs, our solution was to always generate D3 using simple-dftd3 (or another d3 package) and keep the VASP as solely PBE.

[Andrew-S-Rosen]

@tgmaxson, thanks for the comment!

I can confirm that VASP is using the same s6, s8, a1, and a2 parameters as in s-dftd3, which are reproduced below:

simple-dftd3/assets/parameters.toml

Line 48 in 442fbec

d3.bj = {a1=0.4289, s8=0.7875, a2=4.4407, doi="10.1002/jcc.21759"}

There are also some other parameters reported in the OUTCAR, and I don't know what those are, so maybe those are different.

         DFTD3 V3.0 Rev 1
 DF pbe
 parameters
 k1-k3        =   16.0000    1.3333   -4.0000
 Edisp (eV)   -5.46486

 E6    (eV) :    -3.7180
 E8    (eV) :    -1.7468

The difference between the s-dftd3 and the ASE-based API is a fair bit larger, so there must be something else going on there.

For MLIPs, our solution was to always generate D3 using simple-dftd3 (or another d3 package) and keep the VASP as solely PBE.

Good suggestion. I agree.

[awvwgk]

Can you verify that the cutoff values used in VASP are matching the ones used in here? I am adding a feature in #154 to expose the cutoff values in the ASE calculator. In the Python API those can already be changed using set_realspace_cutoff method.

[Andrew-S-Rosen]

@awvwgk, thanks! That looks like the potential culprit.

VDW_RADIUS=50.0
VDW_CNRADIUS=21.167

I believe both of these numbers are different (lower) than s-dftd3.

[ftran1]

The default values for the two-body interaction cutoff are different in VASP (DFTD3 V3.0 Rev 1) and s-dftd3: $\approx 50.202$ Å in VASP and $\approx 31.751$ Å in s-dftd3. If the calculations with the two implementations are done using the same values, then the results are identical.

The default values for the coordination number cutoff are the same in both implementations ($\approx 21.167$ Å).

[thangckt]

The default values for the two-body interaction cutoff are different in VASP (DFTD3 V3.0 Rev 1) and s-dftd3: ≈ 50.202 Å in VASP and ≈ 31.751 Å in s-dftd3. If the calculations with the two implementations are done using the same values, then the results are identical.

The default values for the coordination number cutoff are the same in both implementations ( ≈ 21.167 Å).

@ftran1 are you sure about this?

Based on what I saw in PR #154 and what mentioned by @Andrew-S-Rosen here, the things seem not like that.

[ftran1]

The default values for the two-body interaction cutoff are different in VASP (DFTD3 V3.0 Rev 1) and s-dftd3: ≈ 50.202 Å in VASP and ≈ 31.751 Å in s-dftd3. If the calculations with the two implementations are done using the same values, then the results are identical.
The default values for the coordination number cutoff are the same in both implementations ( ≈ 21.167 Å).

@ftran1 are you sure about this?

Based on what I saw in PR #154 and what mentioned by @Andrew-S-Rosen here, the things seem not like that.

Yes. See file cutoff.f90 of s-dftd3 (value of disp2_default is in Bohr) and VDW_RADIUS for VASP (IVDW=11 or 12). Importantly however, by default the value passed by VASP to s-dftd3 (IVDW=15) is $\approx50.20$ Å, so that when using VASP, the default value is the same ($\approx50.20$ Å) whatever is the implementation of DFT-D3 that is used.

[thangckt]

Hi @ftran1

I see these lines in Python API

    disp2 = self.parameters.realspace_cutoff.get("disp2", 60.0 * Bohr) / Bohr 
    disp3 = self.parameters.realspace_cutoff.get("disp3", 40.0 * Bohr) / Bohr 
    cn = self.parameters.realspace_cutoff.get("cn", 40.0 * Bohr) / Bohr

I think default values for 2 body interaction and coord cutoffs are 60 and 40 Angstroms? Am I wrong?

[awvwgk]

The ASE unit conversion can be somewhat confusing. 60 * Bohr converts 60 Bohr to ASE's default unit which is Angstrom. Since I am expecting the user input in the ASE calculator in Angstrom the value gets converted to Bohr before passing it to dftd3 (via / Bohr).

[thangckt]

@awvwgk Thanks for guiding.

I am still a little confused:

1. By default the cutoff values in s-dftd3 Python is 60 * Bohr ≈31.7 Å, or does it use ≈50.20 Å?
2. Are these default values in s-dftd3 Python the same as in versions before PR Allow modifying cutoff values from ASE calculator #154 ?

Thanks.

[awvwgk]

The default cutoff is 60 Bohr both in Python and Fortran

simple-dftd3/src/dftd3/cutoff.f90

Lines 27 to 28 in 5391f03

	!> Two-body interaction cutoff
	real(wp), parameter :: disp2_default = 60.0_wp

[thangckt]

@awvwgk Thanks you for your guide.

In other D3 package, triple interaction is not available. When one want to compare results from such packages, is it safety to set something like this?

calc_custom.set(realspace_cutoff={"disp2": 31, "disp3": 0, "cn": 21})

is this equivalent to only account 2 body interaction?

[awvwgk]

The ATM term can be disabled by setting s9=0, you can explicitly set this one in the input to the ASE calculator.

simple-dftd3/python/dftd3/ase.py

Lines 47 to 59 in 5391f03

	The params_tweaks dict contains the damping parameters, at least s8, a1 and a2
	must be provided
	======================== =========== ============================================
	Tweakable parameter Default Description
	======================== =========== ============================================
	s6 1.0 Scaling of the dipole-dipole dispersion
	s8 None Scaling of the dipole-quadrupole dispersion
	s9 1.0 Scaling of the three-body dispersion energy
	a1 None Scaling of the critical radii
	a2 None Offset of the critical radii
	alp 14.0 Exponent of the zero damping (ATM only)
	======================== =========== ============================================

Something along the lines of

from dftd3.ase import DFTD3
from dftd3.parameters import get_damping_param

param = get_damping_param("pbe", defaults=["bj"])
param["s9"] = 0.0

calc = DFTD3(damping="d3bj")
calc.set(params_tweaks=param)

Setting the disp3 threshold to zero would work, but is a quite expensive way to remove the ATM since the iterations over all triples will be done at least for L=0 in the current implementation to discard them.

[awvwgk]

In other D3 package, triple interaction is not available. When one want to compare results from such packages, is it safety to set something like this?

I am keeping a comparison between different D3 packages in https://dftd3.readthedocs.io/en/latest/comparison.html, if there is one I missed or the comparison could be made more complete, let me know please.