1 /*
2 * XCFun, an arbitrary order exchange-correlation library
3 * Copyright (C) 2020 Ulf Ekström and contributors.
4 *
5 * This file is part of XCFun.
6 *
7 * This Source Code Form is subject to the terms of the Mozilla Public
8 * License, v. 2.0. If a copy of the MPL was not distributed with this
9 * file, You can obtain one at http://mozilla.org/MPL/2.0/.
10 *
11 * For information on the complete list of contributors to the
12 * XCFun library, see: <https://xcfun.readthedocs.io/>
13 */
14
15 #pragma once
16
17 #include "XCFunctional.hpp"
18 #include "config.hpp"
19
20 // When regularizing we shouldn't touch the higher order
21 // parts of the density, so we need this.
regularize(ctaylor<T,N> & x)22 template <typename T, int N> void regularize(ctaylor<T, N> & x) {
23 if (x < xcfun::XCFUN_TINY_DENSITY)
24 x.set(0, xcfun::XCFUN_TINY_DENSITY);
25 }
26
regularize(T & x)27 template <typename T> static void regularize(T & x) {
28 if (x < xcfun::XCFUN_TINY_DENSITY)
29 x = xcfun::XCFUN_TINY_DENSITY;
30 }
31
32 // Variables for expressing functionals, these are redundant because
33 // different functionals have different needs.
34 // TODO: Make sure all variables are handled in the switch.
35 template <typename T> struct densvars {
36 // Fills all density variables that can be filled from vars. Length of d
37 // depends on vars.
densvarsdensvars38 densvars(const XCFunctional * parent, const T * d) {
39 this->parent = parent;
40 switch (parent->vars) {
41 case XC_A_GAA:
42 gaa = d[1];
43 gab = 0;
44 gbb = 0;
45 gnn = gaa;
46 gss = gaa;
47 gns = gaa;
48 case XC_A:
49 b = 0;
50 n = a;
51 regularize(n);
52 s = n;
53 break;
54 case XC_A_B_GAA_GAB_GBB_TAUA_TAUB:
55 taua = d[5];
56 taub = d[6];
57 tau = taua + taub;
58 case XC_A_B_GAA_GAB_GBB:
59 gaa = d[2];
60 gab = d[3];
61 gbb = d[4];
62 gnn = gaa + 2 * gab + gbb;
63 gss = gaa - 2 * gab + gbb;
64 gns = gaa - gbb;
65 case XC_A_B:
66 a = d[0];
67 regularize(a);
68 b = d[1];
69 regularize(b);
70 n = a + b;
71 s = a - b;
72 break;
73 case XC_N_S_GNN_GNS_GSS_TAUN_TAUS:
74 taua = d[5] + d[6];
75 taub = d[5] - d[6];
76 tau = taua + taub;
77 case XC_N_S_GNN_GNS_GSS:
78 gnn = d[2];
79 gns = d[3];
80 gss = d[4];
81 gaa = 0.25 * (gnn + 2 * gns + gss);
82 gab = 0.25 * (gnn - gss);
83 gbb = 0.25 * (gnn - 2 * gns + gss);
84 case XC_N_S:
85 n = d[0];
86 regularize(n);
87 s = d[1];
88 a = n + s;
89 regularize(a);
90 b = n - s;
91 regularize(b);
92 break;
93 case XC_N_GNN_TAUN:
94 taua = d[2] / 2;
95 taub = d[2] / 2;
96 tau = d[2];
97 case XC_N_GNN:
98 gnn = d[1];
99 gss = 0;
100 gns = 0;
101 gaa = 0.25 * gnn;
102 gab = gaa;
103 gbb = gaa;
104 case XC_N:
105 n = d[0];
106 regularize(n);
107 s = 0;
108 a = 0.5 * n;
109 b = a;
110 break;
111 case XC_N_2ND_TAYLOR:
112 lapa = 0.5 * (d[4] + d[7] + d[9]);
113 lapb = lapa;
114 case XC_N_NX_NY_NZ_TAUN:
115 taua = d[4] / 2;
116 taub = d[4] / 2;
117 tau = d[4];
118 case XC_N_NX_NY_NZ:
119 gnn = d[1] * d[1] + d[2] * d[2] + d[3] * d[3];
120 gss = 0;
121 gns = 0;
122 gaa = 0.25 * gnn;
123 gab = gaa;
124 gbb = gaa;
125 n = d[0];
126 regularize(n);
127 s = 0;
128 a = 0.5 * n;
129 b = a;
130 break;
131 case XC_A_B_2ND_TAYLOR: // a gax gay gaz haxx haxy haxz hayy hayz hazz
132 // 0 1 2 3 4 5 6 7 8 9
133 lapa = d[4] + d[7] + d[9];
134 lapb = d[14] + d[17] + d[19];
135 a = d[0];
136 regularize(a);
137 b = d[10];
138 regularize(b);
139 gaa = d[1] * d[1] + d[2] * d[2] + d[3] * d[3];
140 gab = d[1] * d[11] + d[2] * d[12] + d[3] * d[13];
141 gbb = d[11] * d[11] + d[12] * d[12] + d[13] * d[13];
142 gnn = gaa + 2 * gab + gbb;
143 gss = gaa - 2 * gab + gbb;
144 gns = gaa - gbb;
145 n = a + b;
146 s = a - b;
147 break;
148 case XC_A_B_AX_AY_AZ_BX_BY_BZ_TAUA_TAUB:
149 taua = d[8];
150 taub = d[9];
151 tau = d[8] + d[9];
152 case XC_A_B_AX_AY_AZ_BX_BY_BZ:
153 // 0 1 2 3 4 5 6 7
154 a = d[0];
155 regularize(a);
156 b = d[1];
157 regularize(b);
158 gaa = d[2] * d[2] + d[3] * d[3] + d[4] * d[4];
159 gab = d[2] * d[5] + d[3] * d[6] + d[4] * d[7];
160 gbb = d[5] * d[5] + d[6] * d[6] + d[7] * d[7];
161 gnn = gaa + 2 * gab + gbb;
162 gss = gaa - 2 * gab + gbb;
163 gns = gaa - gbb;
164 n = a + b;
165 s = a - b;
166 break;
167 case XC_N_S_NX_NY_NZ_SX_SY_SZ_TAUN_TAUS:
168 taua = 0.5 * (d[8] + d[9]);
169 taub = 0.5 * (d[8] - d[9]);
170 tau = d[8];
171 case XC_N_S_NX_NY_NZ_SX_SY_SZ:
172 // 0 1 2 3 4 5 6 7
173 n = d[0];
174 regularize(n);
175 s = d[1];
176 a = 0.5 * (n + s);
177 regularize(a);
178 b = 0.5 * (n - s);
179 regularize(b);
180 gnn = d[2] * d[2] + d[3] * d[3] + d[4] * d[4];
181 gss = d[5] * d[5] + d[6] * d[6] + d[7] * d[7];
182 gns = d[2] * d[5] + d[3] * d[6] + d[4] * d[7];
183 gaa = 0.25 * (gnn + 2 * gns + gss);
184 gab = 0.25 * (gnn - gss);
185 gbb = 0.25 * (gnn - 2 * gns + gss);
186 break;
187 case XC_A_B_GAA_GAB_GBB_LAPA_LAPB_TAUA_TAUB_JPAA_JPBB:
188 jpaa = d[9];
189 jpbb = d[10];
190 case XC_A_B_GAA_GAB_GBB_LAPA_LAPB_TAUA_TAUB:
191 lapa = d[5];
192 lapb = d[6];
193 taua = d[7];
194 taub = d[8];
195 tau = taua + taub;
196 gaa = d[2];
197 gab = d[3];
198 gbb = d[4];
199 gnn = gaa + 2 * gab + gbb;
200 gss = gaa - 2 * gab + gbb;
201 gns = gaa - gbb;
202 a = d[0];
203 regularize(a);
204 b = d[1];
205 regularize(b);
206 n = a + b;
207 s = a - b;
208 break;
209 default:
210 xcfun::die("Illegal/Not yet implemented vars value in densvars()",
211 parent->vars);
212 }
213 zeta = s / n;
214 r_s = pow(3.0 / (n * 4.0 * M_PI), 1.0 / 3.0); // (3/4pi)^1/3*n^(-1/3) !check
215 n_m13 = pow(n, -1.0 / 3.0);
216 a_43 = pow(a, 4.0 / 3.0);
217 b_43 = pow(b, 4.0 / 3.0);
218 }
219
220 const XCFunctional * parent{nullptr};
get_paramdensvars221 double get_param(enum xc_parameter p) const { return parent->settings[p]; }
222
223 T a{static_cast<T>(0)};
224 T b{static_cast<T>(0)};
225 T gaa{static_cast<T>(0)};
226 T gab{static_cast<T>(0)};
227 T gbb{static_cast<T>(0)};
228 T n{static_cast<T>(0)}; /// na+nb
229 T s{static_cast<T>(0)}; /// na - nb
230 T gnn{static_cast<T>(0)}; /// (grad n) ^ 2
231 T gns{static_cast<T>(0)}; /// (grad n).(grad s)
232 T gss{static_cast<T>(0)}; /// (grad s) ^ 2
233 T tau{static_cast<T>(0)}; /// Kinetic energy density.
234 T taua{static_cast<T>(0)}; /// Alpha kinetic energy density.
235 T taub{static_cast<T>(0)}; /// Beta kinetic energy density.
236 T lapa{static_cast<T>(0)}; /// Alpha Laplacian density.
237 T lapb{static_cast<T>(0)}; /// Beta Laplacian density.
238 T zeta{static_cast<T>(0)}; /// s/n
239 T r_s{static_cast<T>(0)}; /// (3/4pi)^1/3*n^(-1/3)
240 T n_m13{static_cast<T>(0)}; /// pow(n,-1.0/3.0)
241 T a_43{static_cast<T>(0)};
242 T b_43{static_cast<T>(0)}; /// pow(a,4.0/3.0), pow(b,4.0/3.0)
243 T jpaa{static_cast<T>(0)}; /// square of the alpha paramagnetic current vector.
244 T jpbb{static_cast<T>(0)}; /// square of the beta paramagnetic current vector.
245 };
246