Coverage for eminus/xc/lda_c_chachiyo.py: 100.00%
36 statements
« prev ^ index » next coverage.py v7.6.9, created at 2024-12-18 08:43 +0000
« prev ^ index » next coverage.py v7.6.9, created at 2024-12-18 08:43 +0000
1# SPDX-FileCopyrightText: 2021 The eminus developers
2# SPDX-License-Identifier: Apache-2.0
3"""Chachiyo LDA correlation.
5Reference: J. Chem. Phys. 145, 021101.
6"""
8import numpy as np
11def lda_c_chachiyo(n, **kwargs):
12 """Chachiyo parametrization of the correlation functional (spin-paired).
14 Corresponds to the functional with the label LDA_C_CHACHIYO and ID 287 in Libxc.
16 Reference: J. Chem. Phys. 145, 021101.
18 Args:
19 n: Real-space electronic density.
21 Keyword Args:
22 **kwargs: Throwaway arguments.
24 Returns:
25 Chachiyo correlation energy density and potential.
26 """
27 a = -0.01554535 # (np.log(2) - 1) / (2 * np.pi**2)
28 b = 20.4562557
30 rs = (3 / (4 * np.pi * n)) ** (1 / 3)
31 rs2 = rs**2
32 ecinner = 1 + b / rs + b / rs2
34 ec = a * np.log(ecinner)
36 vc = ec + a * b * (2 + rs) / (3 * (b + b * rs + rs2))
37 return ec, np.array([vc]), None
40def chachiyo_scaling(zeta):
41 """Weighting factor between the paramagnetic and the ferromagnetic case.
43 Reference: J. Chem. Phys. 145, 021101.
45 Args:
46 zeta: Relative spin polarization.
48 Returns:
49 Weighting factor and its derivative.
50 """
51 fzeta = ((1 + zeta) ** (4 / 3) + (1 - zeta) ** (4 / 3) - 2) / (2 * (2 ** (1 / 3) - 1))
53 dfdzeta = (2 * (1 - zeta) ** (1 / 3) - 2 * (1 + zeta) ** (1 / 3)) / (3 - 3 * 2 ** (1 / 3))
54 return fzeta, dfdzeta
57def lda_c_chachiyo_spin(n, zeta, weight_function=chachiyo_scaling, **kwargs):
58 """Chachiyo parametrization of the correlation functional (spin-polarized).
60 Corresponds to the functional with the label LDA_C_CHACHIYO and ID 287 in Libxc.
62 Reference: J. Chem. Phys. 145, 021101.
64 Args:
65 n: Real-space electronic density.
66 zeta: Relative spin polarization.
68 Keyword Args:
69 weight_function: Functional function.
70 **kwargs: Throwaway arguments.
72 Returns:
73 Chachiyo correlation energy density and potential.
74 """
75 a0 = -0.01554535 # (np.log(2) - 1) / (2 * np.pi**2)
76 a1 = -0.007772675 # (np.log(2) - 1) / (4 * np.pi**2)
77 b0 = 20.4562557
78 b1 = 27.4203609
80 rs = (3 / (4 * np.pi * n)) ** (1 / 3)
81 rs2 = rs**2
83 fzeta, dfdzeta = weight_function(zeta)
85 ec0inner = 1 + b0 / rs + b0 / rs2
86 ec1inner = 1 + b1 / rs + b1 / rs2
87 ec0 = a0 * np.log(ec0inner)
88 ec1 = a1 * np.log(ec1inner)
90 ec = ec0 + (ec1 - ec0) * fzeta
92 factor = -1 / rs2 - 2 / rs**3
93 dec0drs = a0 / ec0inner * b0 * factor
94 dec1drs = a1 / ec1inner * b1 * factor
95 decdrs = dec0drs + (dec1drs - dec0drs) * fzeta
96 prefactor = ec - rs / 3 * decdrs
97 decdf = (ec1 - ec0) * dfdzeta
99 vc_up = prefactor + decdf * (1 - zeta)
100 vc_dw = prefactor - decdf * (1 + zeta)
101 return ec, np.array([vc_up, vc_dw]), None