Systematic differences between VASP and libxc implementations of LDA

[srdjan_stavric]

Dear VASP developers and users,

I am studying a 2D ferromagnet (CrI₃) and encountered a puzzling discrepancy. Using VASP 6.5.1, I perform spin-polarized calculations (vasp_std) with the LDA-PW92 functional
and compute the energy difference ΔE = E(FM) − E(AFM), for a unit cell containing two Cr atoms.

I obtain significantly different ΔE values depending on whether LDA-PW92 is used via VASP’s internal implementation (GGA = PW92) or via libxc. The difference is systematic and on the order of 20–30%,
and it persists across all plane-wave energy cutoffs tested (see attached figure). I observe the same behavior when using the LDA-PZ functional.

Since both approaches nominally use the same xc functional, I would expect identical or at least very similar results. Could someone clarify the origin of this discrepancy?

Thank you for your time,
Srdjan

[andreas.singraber]

Hello Srdjan,

Thank you for the detailed explanation of your issue. We will have a look into it but please attach the relevant input and output files so we can easily reproduce the behavior (as outlined in the forum posting guidelines). Thank you!

All the best,
Andreas Singraber

[srdjan_stavric]

Dear Andreas,

Thank you for the quick reply.

I have uploaded three tar files. One contains the input files for calculations using VASP’s internal (native) implementation of the xc functional, and the other contains the corresponding inputs using the libxc implementation.
Each directory includes a run.sh script that runs all calculations for the tested energy cutoffs (via nohup bash run.sh). These calculations were performed on my laptop and are computationally lightweight.

In addition, the third .tar file contains two representative OUTCAR files for the FM state: one from the native implementation and one from the libxc implementation.
Uploading all OUTCAR files would be rather cumbersome, but I can provide more if needed.

Thank you for your support.

Best regards,
Srdjan

[andreas.singraber]

Hello Srdjan,

thanks for uploading the files! My colleague Fabien Tran had a closer look and was able to identify the origin of the inconsistencies. It turns out that enabling LASPH=.TRUE. when running the Libxc implementation of LDA-PW92 will give you consistent results. Fabien gave me an explanation which I am citing here:

For the GGA functionals implemented in VASP, there is a separation between the LDA and GGA components. LASPH= .FALSE. (the default) is an approximation that is applied to the GGA component. For the VASP implementation of LDA-PW92, the value of LASPH does not matter since LDA-PW92 has no GGA component. For the Libxc implementation it matters since this is the reverse: LDA-PW92 is processed as a GGA with no LDA component (confusing), such that LASPH= .FALSE. has an effect on the whole functional.

Hence, if you would repeat your Libxc-based simulations with LASPH enabled you should see much better agreement. Generally, it is recommended to enable this feature for better accuracy.

Fabien also mentioned another difference in the implementations which may give you still some slight deviation of the results (but on a much smaller scale): the LDA implementation in VASP includes a modification which accounts for relativistic effects while LDA from Libxc carries out the original non-relativistic method.

All the best,
Andreas Singraber

[srdjan_stavric]

Dear Andreas,

Thank you for the clarification, and please pass my thanks to Fabien as well. I included LASPH = .TRUE. and recalculated the energy differences (E(FM) - E(AFM)).
With this setting, the two LDA-PW92 implementations now give very similar results - the difference is on the order of 0.2 meV (compared to ~8 meV previously).

Best regards,
Srdjan