1 /* Ergo, version 3.8, a program for linear scaling electronic structure
2 * calculations.
3 * Copyright (C) 2019 Elias Rudberg, Emanuel H. Rubensson, Pawel Salek,
4 * and Anastasia Kruchinina.
5 *
6 * This program is free software: you can redistribute it and/or modify
7 * it under the terms of the GNU General Public License as published by
8 * the Free Software Foundation, either version 3 of the License, or
9 * (at your option) any later version.
10 *
11 * This program is distributed in the hope that it will be useful,
12 * but WITHOUT ANY WARRANTY; without even the implied warranty of
13 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14 * GNU General Public License for more details.
15 *
16 * You should have received a copy of the GNU General Public License
17 * along with this program. If not, see <http://www.gnu.org/licenses/>.
18 *
19 * Primary academic reference:
20 * Ergo: An open-source program for linear-scaling electronic structure
21 * calculations,
22 * Elias Rudberg, Emanuel H. Rubensson, Pawel Salek, and Anastasia
23 * Kruchinina,
24 * SoftwareX 7, 107 (2018),
25 * <http://dx.doi.org/10.1016/j.softx.2018.03.005>
26 *
27 * For further information about Ergo, see <http://www.ergoscf.org>.
28 */
29
30 /*-*-mode: C; c-indentation-style: "bsd"; c-basic-offset: 4; -*-*/
31 /** @file fun-kt.c
32 Implementation of KT functional and its derivatives.
33 Or exactly: KT GGA correction to the functional for KT1,KT2 total functional
34 energy is E_LDA+E_KT).
35 Reference: Keal, Tozer, J. Chem. Phys., 119, 3015 (2003).
36 GAMMA is included in the KTx definition in fun-gga.c
37 implemented by Dave Wilson (davidwi@kjemi.uio.no)
38 NOTE:
39 this file may seem unnecessarily complex but the structure does pay off
40 when implementing multiple functionals depending on different parameters.
41 */
42
43 #define FOURTH_ORDER_DERIVATIVES
44
45 #include <math.h>
46 #include <stddef.h>
47
48 #define __CVERSION__
49
50 #include "functionals.h"
51
52 /* INTERFACE PART */
kt_isgga(void)53 static int kt_isgga(void) { return 1; }
54 static int kt_read(const char* conf_line);
55 static real kt_energy(const FunDensProp* dens_prop);
56 static void kt_first(FunFirstFuncDrv *ds, real factor,
57 const FunDensProp* dens_prop);
58 static void kt_second(FunSecondFuncDrv *ds, real factor,
59 const FunDensProp* dens_prop);
60 static void kt_third(FunThirdFuncDrv *ds, real factor,
61 const FunDensProp* dens_prop);
62 #ifdef FOURTH_ORDER_DERIVATIVES
63 static void kt_fourth(FunFourthFuncDrv *ds, real factor,
64 const FunDensProp* dens_prop);
65 #endif
66
67
68 Functional KTFunctional = {
69 "KT", /* name */
70 kt_isgga, /* gga-corrected */
71 kt_read, /* set bloody common blocks */
72 NULL, /* reporter */
73 kt_energy,
74 kt_first,
75 kt_second,
76 kt_third
77 #ifdef FOURTH_ORDER_DERIVATIVES
78 ,kt_fourth
79 #endif
80 };
81
82 /* IMPLEMENTATION PART */
83
84 static int
kt_read(const char * conf_line)85 kt_read(const char* conf_line)
86 {
87 fun_set_hf_weight(0.0);
88 return 1;
89 }
90
91 /* kt_energy:
92 note that in reality E_KT = E_KT,alpha + E_KT,beta
93 i.e the energy is linear in alpha and beta densities.
94
95 KT threshold is needed to avoid numerical problems on 0/0
96 divisions. The problems are small but it is better to be on the
97 safe side.
98 */
99 static const real KT_THRESHOLD = 1e-14;
100 static const real DELTA = 0.1;
101 static real
kt_energy(const FunDensProp * dp)102 kt_energy(const FunDensProp* dp)
103 {
104 real ea,eb;
105 if (dp->rhob<KT_THRESHOLD)
106 eb = 0.0;
107 else {
108 real xb = dp->gradb;
109 real rb = POW(dp->rhob,4.0/3.0);
110 real denomb = rb + DELTA;
111 eb = xb*xb/denomb;
112 }
113 if (dp->rhoa<KT_THRESHOLD)
114 ea=0;
115 else {
116 real xa = dp->grada;
117 real ra = POW(dp->rhoa,4.0/3.0);
118 real denoma = ra + DELTA;
119 ea = xa*xa/denoma;
120 }
121 return (ea+eb);
122 }
123
124 static void
kt_first(FunFirstFuncDrv * ds,real factor,const FunDensProp * dp)125 kt_first(FunFirstFuncDrv *ds, real factor, const FunDensProp* dp)
126 {
127 real xa, ra43, ra13;
128 real xb, rb43, rb13;
129 real denoma, denoma2;
130 real denomb, denomb2;
131 real faZ, faR;
132 real fbZ, fbR;
133
134 if (dp->rhoa >KT_THRESHOLD) {
135 xa = dp->grada;
136 ra43 = POW(dp->rhoa,4.0/3.0);
137 ra13 = POW(dp->rhoa,1.0/3.0);
138 denoma = ra43 + DELTA;
139 denoma2= denoma*denoma;
140 faR = -4.0/3.0*xa*xa*ra13/denoma2;
141 faZ = 2.0*xa/denoma;
142 ds->df1000 += factor*faR;
143 ds->df0010 += factor*faZ;
144 }
145 if (dp->rhob >KT_THRESHOLD) {
146 xb = dp->gradb;
147 rb43 = POW(dp->rhob,4.0/3.0);
148 rb13 = POW(dp->rhob,1.0/3.0);
149 denomb = rb43 + DELTA;
150 denomb2= denomb*denomb;
151 fbR = -4.0/3.0*xb*xb*rb13/denomb2;
152 fbZ = 2.0*xb/denomb;
153 ds->df0100 += factor*fbR;
154 ds->df0001 += factor*fbZ;
155 }
156 }
157
158
159 static void
kt_second(FunSecondFuncDrv * ds,real factor,const FunDensProp * dp)160 kt_second(FunSecondFuncDrv *ds, real factor, const FunDensProp* dp)
161 {
162 real xa, ra43, ra13, ram2;
163 real xb, rb43, rb13, rbm2;
164 real denoma, denoma2, t1a;
165 real denomb, denomb2, t1b;
166 real faR, faZ, faRZ, faZZ, faRR;
167 real fbR, fbZ, fbRZ, fbZZ, fbRR;
168
169 if (dp->rhoa >KT_THRESHOLD) {
170 xa = dp->grada;
171 ra43 = POW(dp->rhoa,4.0/3.0);
172 ra13 = POW(dp->rhoa,1.0/3.0);
173 ram2 = POW(dp->rhoa,-2.0/3.0);
174 denoma = ra43 + DELTA;
175 denoma2= denoma*denoma;
176 faR = -(4.0/3.0)*xa*xa*ra13/denoma2;
177 faZ = 2.0*xa/denoma;
178 t1a = 4.0/9.0*xa*xa/denoma2;
179 faRR = -ram2 + 8.0*ra13*ra13/denoma;
180 faRZ = -8.0/3.0*xa*ra13/denoma2;
181 faZZ = 2.0/denoma;
182 ds->df1000 += factor*faR;
183 ds->df0010 += factor*faZ;
184 ds->df1010 += factor*faRZ;
185 ds->df2000 += factor*faRR*t1a;
186 ds->df0020 += factor*faZZ;
187 }
188 if (dp->rhob >KT_THRESHOLD) {
189 xb = dp->gradb;
190 rb43 = POW(dp->rhob,4.0/3.0);
191 rb13 = POW(dp->rhob,1.0/3.0);
192 rbm2 = POW(dp->rhob,-2.0/3.0);
193 denomb = rb43 + DELTA;
194 denomb2= denomb*denomb;
195 fbR = -4.0/3.0*xb*xb*rb13/denomb2;
196 fbZ = 2.0*xb/denomb;
197 t1b = 4.0/9.0*xb*xb/denomb2;
198 fbRR = -rbm2 + 8.0*rb13*rb13/denomb;
199 fbRZ = -8.0/3.0*xb*rb13/denomb2;
200 fbZZ = 2.0/denomb;
201 ds->df0100 += factor*fbR;
202 ds->df0001 += factor*fbZ;
203 ds->df0101 += factor*fbRZ;
204 ds->df0200 += factor*fbRR*t1b;
205 ds->df0002 += factor*fbZZ;
206 }
207 }
208
209
210 static void
kt_third(FunThirdFuncDrv * ds,real factor,const FunDensProp * dp)211 kt_third(FunThirdFuncDrv *ds, real factor, const FunDensProp* dp)
212 {
213 real ra, xa, ra43, ra13, ram1, ram2;
214 real rb, xb, rb43, rb13, rbm1, rbm2;
215 real xad2, t1a, t2a;
216 real xbd2, t1b, t2b;
217 real denoma, denoma2;
218 real denomb, denomb2;
219 real faR, faZ, faRR, faZZ, faRZ;
220 real fbR, fbZ, fbRR, fbZZ, fbRZ;
221 real faRRR, faRRZ, faRZZ, faZZZ;
222 real fbRRR, fbRRZ, fbRZZ, fbZZZ;
223
224 if (dp->rhoa >KT_THRESHOLD) {
225 xa = dp->grada;
226 ra = dp->rhoa;
227 ra13 = POW(dp->rhoa,1.0/3.0);
228 ra43 = ra13*ra;
229 ram2 = ra13/ra;
230 ram1 = ram2*ra13;
231 denoma = ra43 + DELTA;
232 denoma2= denoma*denoma;
233 xad2 = xa/denoma2;
234 t1a = 8.0/9.0*xad2*xa*ram1;
235 t2a = 1.0 - 8.0*ra43/denoma;
236 faR = -4.0/3.0*xad2*xa*ra13;
237 faZ = 2.0*xa/denoma;
238 faRR = -4.0/9.0*xad2*xa*ram2;
239 faRZ = -8.0/3.0*xad2*ra13;
240 faZZ = 2.0/denoma;
241 faZZZ = 0.0;
242 faRZZ = -8.0/3.0*ra13/denoma2;
243 faRRZ = -8.0/9.0*xad2*ram2;
244 faRZZ = -8.0/3.0*ra13/denoma2;
245 faRRR = 1.0/(3.0*ra43) + 4.0/denoma -
246 16.0*ra43/denoma2;
247 ds->df1000 += factor*faR;
248 ds->df0010 += factor*faZ;
249 ds->df1010 += factor*faRZ;
250 ds->df2000 += factor*faRR*t2a;
251 ds->df0020 += factor*faZZ;
252 ds->df3000 += factor*faRRR*t1a;
253 ds->df0030 += factor*faZZZ;
254 ds->df2010 += factor*faRRZ*t2a;
255 ds->df1020 += factor*faRZZ;
256 }
257 if (dp->rhob >KT_THRESHOLD) {
258 xb = dp->gradb;
259 rb = dp->rhob;
260 rb13 = POW(dp->rhob,1.0/3.0);
261 rb43 = rb13*rb;
262 rbm2 = rb13/rb;
263 rbm1 = rbm2*rb13;
264 denomb = rb43 + DELTA;
265 denomb2= denomb*denomb;
266 xbd2 = xb/denomb2;
267 t1b = 8.0/9.0*xbd2*xb*rbm1;
268 t2b = 1.0 - 8.0*rb43/denomb;
269 fbR = -4.0/3.0*xbd2*xb*rb13;
270 fbZ = 2.0*xb/denomb;
271 fbRR = -4.0/9.0*xbd2*xb*rbm2;
272 fbRZ = -8.0/3.0*xbd2*rb13;
273 fbZZ = 2.0/denomb;
274 fbZZZ = 0.0;
275 fbRRZ = -8.0/9.0*xbd2*rbm2;
276 fbRZZ = -8.0/3.0*rb13/denomb2;
277 fbRRR = 1.0/(3.0*rb43) + 4.0/denomb -
278 16.0*rb43/denomb2;
279 ds->df0100 += factor*fbR;
280 ds->df0001 += factor*fbZ;
281 ds->df0101 += factor*fbRZ;
282 ds->df0200 += factor*fbRR*t2b;
283 ds->df0002 += factor*fbZZ;
284 ds->df0300 += factor*fbRRR*t1b;
285 ds->df0003 += factor*fbZZZ;
286 ds->df0201 += factor*fbRRZ*t2b;
287 ds->df0102 += factor*fbRZZ;
288 }
289 }
290
291
292 #ifdef FOURTH_ORDER_DERIVATIVES
293 static void
kt_fourth(FunFourthFuncDrv * ds,real factor,const FunDensProp * dp)294 kt_fourth(FunFourthFuncDrv *ds, real factor, const FunDensProp* dp)
295 {
296 real ra, xa, ra43, ra13, ra23, ram1, ram2;
297 real rb, xb, rb43, rb13, rb23, rbm1, rbm2;
298 real xad2, xa2d2, t1a, t2a, t3a;
299 real xbd2, xb2d2, t1b, t2b, t3b;
300 real denoma, denoma2, denoma3;
301 real denomb, denomb2, denomb3;
302 real faR, faZ, faRR, faZZ, faRZ;
303 real fbR, fbZ, fbRR, fbZZ, fbRZ;
304 real faRRR, faRRZ, faRZZ, faZZZ;
305 real fbRRR, fbRRZ, fbRZZ, fbZZZ;
306 real faRRRR, faRRRZ, faRRZZ, faRZZZ, faZZZZ;
307 real fbRRRR, fbRRRZ, fbRRZZ, fbRZZZ, fbZZZZ;
308
309 if (dp->rhoa >KT_THRESHOLD) {
310 xa = dp->grada;
311 ra = dp->rhoa;
312 ra13 = POW(dp->rhoa,1.0/3.0);
313 ra43 = ra13*ra;
314 ra23 = ra13*ra13;
315 ram2 = ra13/ra;
316 ram1 = ram2*ra13;
317 denoma = ra43 + DELTA;
318 denoma2= denoma*denoma;
319 denoma3= denoma2*denoma;
320 xad2 = xa/denoma2;
321 xa2d2 = xa*xa/denoma2;
322 t1a = 1.0/(3.0*ra43) + 4.0/denoma
323 - 16.0*ra43/denoma2;
324 t2a = 1.0 - 8.0*ra43/denoma;
325 t3a = -8.0/27.0*xa2d2*ram1;
326 faR = -4.0/3.0*xad2*xa*ra13;
327 faZ = 2.0*xa/denoma;
328 faRR = -4.0/9.0*xad2*xa*ram2;
329 faRZ = -8.0/3.0*xad2*ra13;
330 faZZ = 2.0/denoma;
331 faZZZ = 0.0;
332 faRZZ = -8.0/3.0*ra13/denoma2;
333 faRRZ = -8.0/9.0*xad2*ram2;
334 faRZZ = -8.0/3.0*ra13/denoma2;
335 faRRR = 8.0/9.0*xad2*xa*ram1;
336 faZZZZ = 0.0;
337 faRZZZ = 0.0;
338 faRRZZ = -8.0/(9.0*ra23*denoma2);
339 faRRRZ = 16.0/9.0*xad2*ram1;
340 faRRRR = 5.0/(3.0*ra*ra43)
341 + 20.0/(3.0*ra*denoma)
342 + 96.0*ra13/denoma2
343 - 256.0*ra*ra23/denoma3;
344 ds->df1000 += factor*faR;
345 ds->df0010 += factor*faZ;
346 ds->df1010 += factor*faRZ;
347 ds->df2000 += factor*faRR*t2a;
348 ds->df0020 += factor*faZZ;
349 ds->df3000 += factor*faRRR*t1a;
350 ds->df0030 += factor*faZZZ;
351 ds->df2010 += factor*faRRZ*t2a;
352 ds->df1020 += factor*faRZZ;
353 ds->df4000 += factor*faRRRR*t3a;
354 ds->df0040 += factor*faZZZZ;
355 ds->df3010 += factor*faRRRZ*t1a;
356 ds->df1030 += factor*faRZZZ;
357 ds->df2020 += factor*faRRZZ*t2a;
358 }
359 if (dp->rhob >KT_THRESHOLD) {
360 xb = dp->gradb;
361 rb = dp->rhob;
362 rb13 = POW(dp->rhob,1.0/3.0);
363 rb43 = rb13*rb;
364 rb23 = rb13*rb13;
365 rbm2 = rb13/rb;
366 rbm1 = rbm2*rb13;
367 denomb = rb43 + DELTA;
368 denomb2= denomb*denomb;
369 denomb3= denomb2*denomb;
370 xbd2 = xb/denomb2;
371 xb2d2 = xb*xb/denomb2;
372 t1b = 1.0/(3.0*rb43) + 4.0/denomb
373 - 16.0*rb43/denomb2;
374 t2b = 1.0 - 8.0*rb43/denomb;
375 t3b = -8.0/27.0*xb2d2/rb13;
376 fbR = -4.0/3.0*xbd2*xb*rb13;
377 fbZ = 2.0*xb/denomb;
378 fbRR = -4.0/9.0*xbd2*xb*rbm2;
379 fbRZ = -8.0/3.0*xbd2*rb13;
380 fbZZ = 2.0/denomb;
381 fbZZZ = 0.0;
382 fbRRZ = -8.0/9.0*xbd2*rbm2;
383 fbRZZ = -8.0/3.0*rb13/denomb2;
384 fbRRR = 8.0/9.0*xbd2*xb*rbm1;
385 fbZZZZ = 0.0;
386 fbRZZZ = 0.0;
387 fbRRZZ = -8.0/(9.0*rb23*denomb2);
388 fbRRRZ = 16.0/9.0*xbd2*rbm1;
389 fbRRRR = +5.0/(3.0*rb*rb43)
390 + 20.0/(3.0*rb*denomb)
391 + 96.0*rb13/denomb2
392 - 256.0*rb*rb23/denomb3;
393 ds->df0100 += factor*fbR;
394 ds->df0001 += factor*fbZ;
395 ds->df0101 += factor*fbRZ;
396 ds->df0200 += factor*fbRR*t2b;
397 ds->df0002 += factor*fbZZ;
398 ds->df0300 += factor*fbRRR*t1b;
399 ds->df0003 += factor*fbZZZ;
400 ds->df0201 += factor*fbRRZ*t2b;
401 ds->df0102 += factor*fbRZZ;
402 ds->df0400 += factor*fbRRRR*t3b;
403 ds->df0004 += factor*fbZZZZ;
404 ds->df0301 += factor*fbRRRZ*t1b;
405 ds->df0103 += factor*fbRZZZ;
406 ds->df0202 += factor*fbRRZZ*t2b;
407 }
408 }
409 #endif
410
411
412