Bugs in the LDA_XC_GDSMFB functional implementation

Hi,

I found some inconsistencies in the implementation of the
LDA_XC_GDSMFB functional. The Libxc values do not agree
with the reference implementation for some parameters.

My tests indicate that for spin-unpolarized (\zeta=0)
calculations the underlying issues have no consequences.
As far as I know, the Libxc LDA_XC_GDSMFB functional was only
used for spin-unpolarized calculations in literature.
As such, the issues may not affect published results.

I assume most of the issues originate from the fact
that in the final published article some important equations
and parameters were moved to the supplemental material.
On the other hand, the arXiv version has
all the needed equations in the main text.

I converted the Maple LDA_XC_GDSMFB functional definition
to Python (see LDA_XC_GDSMFB(rs,z,params_a_T,params_a_thetaParam) Python function).
Libxc and my Python reimplementation agree (see Tab.1).
I analyzed the LDA_XC_GDSMFB Maple code and found five issues.
I fixed the issues in a modification of my Python implementation
(see LDA_XC_GDSMFB_MOD(rs,z,params_a_T,params_a_thetaParam) Python function).
Each change to the original Maple version is documented in the Python source code.
The modified Python version agrees with the reference
implementation (see Tab.1) for all tested parameters.
Please find the attached Python file LDA_XC_GDSMFB.py containing
LDA_XC_GDSMFB(rs,z,params_a_T,params_a_thetaParam)
and LDA_XC_GDSMFB_MOD(rs,z,params_a_T,params_a_thetaParam) Python functions
for more information.
The script can be executed to produce Tab.1.

Next, I will outline the issues and my proposed fixes.

Issues/Proposed changes

1.) Change g parameters for LDA_XC_GDSMFB
Libxc currently uses g = [2/3, -0.0139261, 0.183208].
These are the KSDT parameters.
For LDA_XC_GDSMFB these parameters should be
g = [2/3,h1,h2] = [2/3,3.18747258,7.74662802].
The parameters h1 and h2 are given in the supplemental material, Tab.1
Note, that the parameters have the wrong sign in the table.
However, they need to have a positive sign, see reference implementation.

2.) Change l parameters LDA_XC_GDSMFB
Libxc currently used l = np.array([1.064009, 0.572565]).
These are the KSDT parameters.
For LDA_XC_GDSMFB these parameters should be
l = [\lambda_1,\lambda_2] = [1.85909536,0].
The parameters \lambda_1 and \lambda_2 are given in the supplemental material, Tab.1
However, they need to have a positive sign, see reference implementation.

3.) Change power of the mtt function
I think, the reduced temperature \theta should be defined as
\theta = 2 (4/(9\pi))^{(2/3)} T r_{s}^2 (1+\zeta)^{(-2/3)}.
Libxc currently uses
mtt = 2(4/(9\pi))^{(2/3)} T r_{s}^2 (1 + \theta_{param} \zeta)^{(2/3)}.
The final exponent has the wrong sign.

4.) Fix \theta_{param} of the mtt function
Libxc currently uses
mtt = 2(4/(9\pi))^{(2/3)} T r_{s}^2 (1 + \theta_{param} \zeta)^{(2/3)}.
In my opinion this expression is only correct for LDA_XC_GDSMFB when
the \theta_{param} is always set to 1 (besides the final exponent switching sign, see 3.)).

5.) Rewrite f function in terms of \theta_{0} and \theta_{1}
Libxc currently uses the f function as defined in the respective Maple file.
I believe that the correct \theta_{0} (used by f_{xc}^{\zeta=0} and \phi)
utilizing the mtt function should be calculated as
\theta_{0} = mtt(r_{s},\zeta) (1+\zeta)^{(2/3)}.
Then correct \theta_{1} (used by f_{xc}^{\zeta=1}) is then defined as \theta_{1} = \theta_{0} 2^{(2/3)}.
Please, inspect the LDA_XC_GDSMFB_MOD(rs,z,params_a_T,params_a_thetaParam)
Python functions in the attached Python file LDA_XC_GDSMFB.py for more details.

Computational details
- tag reference: https://github.com/agbonitz/xc_functional/blob/master/fxc.py
- tag pylibxc pylibxc.version.__version__='6.2.2'
- tag python: Python reimplementation of GDSMFB based on KSDT Libxc Maple code
- tag python-mod: modified Python implementation (including all proposed fixes)

Tab.1: LDA_XC_GDSMFB comparison.

r_{s}	\zeta	T	f_{\mathrm{xc}}^{\mathrm{reference}}	f_{\mathrm{xc}}^{\mathrm{pylibxc}}	f_{\mathrm{xc}}^{\mathrm{python}}	f_{\mathrm{xc}}^{\mathrm{python-mod}}	f_{\mathrm{xc}}^{\mathrm{reference}} - f_{\mathrm{xc}}^{\mathrm{pylibxc}}	f_{\mathrm{xc}}^{\mathrm{pylibxc}} - f_{\mathrm{xc}}^{\mathrm{python}}	f_{\mathrm{xc}}^{\mathrm{reference}} - f_{\mathrm{xc}}^{\mathrm{python-mod}}
0.4	0.0	1e-08	-1.224080	-1.224080	-1.224080	-1.224080	7.25e-08	-4.44e-16	7.25e-08
0.8	0.4	1e-06	-0.653781	-0.653189	-0.653189	-0.653781	-5.92e-04	0.00e+00	2.59e-08
1.0	0.8	0.001	-0.572615	-0.571629	-0.571629	-0.572615	-9.86e-04	0.00e+00	2.30e-08
2.0	0.9	3.0	-0.125520	-0.125520	-0.125520	-0.125520	3.06e-09	0.00e+00	2.46e-09
3.0	1.0	10.0	-0.036241	-0.036241	-0.036241	-0.036241	3.01e-09	-1.95e-15	3.01e-09

I hope this information is sufficient to adjust the Libxc LDA_XC_GDSMFB functional.

Best,
Sebastian

LDA_XC_GDSMFB.py