1 /*
2
3
4 !
5 ! Dalton, a molecular electronic structure program
6 ! Copyright (C) by the authors of Dalton.
7 !
8 ! This program is free software; you can redistribute it and/or
9 ! modify it under the terms of the GNU Lesser General Public
10 ! License version 2.1 as published by the Free Software Foundation.
11 !
12 ! This program is distributed in the hope that it will be useful,
13 ! but WITHOUT ANY WARRANTY; without even the implied warranty of
14 ! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
15 ! Lesser General Public License for more details.
16 !
17 ! If a copy of the GNU LGPL v2.1 was not distributed with this
18 ! code, you can obtain one at https://www.gnu.org/licenses/old-licenses/lgpl-2.1.en.html.
19 !
20
21 !
22
23 */
24 /*-*-mode: C; c-indentation-style: "bsd"; c-basic-offset: 4; -*-*/
25 /* fun-optx.c:
26 implementation of OPTX exchange functional and its derivatives.
27 #### this is just the gradient corrected term for KT3 functional####
28 Reference: N.C. Handy and A.J. Cohen, Mol. Phys., 99, 403 (2001).
29 Keal, Tozer, in press (2004).
30 implemented by Dave Wilson (davidwi@kjemi.uio.no)
31 NOTE:
32 this file may seem unnecessarily complex but the structure does pay off
33 when implementing multiple functionals depending on different parameters.
34 */
35
36 /* strictly conform to XOPEN ANSI C standard */
37 #if !defined(SYS_DEC)
38 /* XOPEN compliance is missing on old Tru64 4.0E Alphas */
39 #define _XOPEN_SOURCE 500
40 #define _XOPEN_SOURCE_EXTENDED 1
41 #endif
42
43 #include <math.h>
44 #include <stddef.h>
45 #include "general.h"
46
47 #define __CVERSION__
48
49 #include "functionals.h"
50
51 /* INTERFACE PART */
optx_isgga(void)52 static integer optx_isgga(void) { return 1; }
53 static integer optxread(const char* conf_line);
54 static real optx_energy(const FunDensProp* dens_prop);
55 static void optx_first(FunFirstFuncDrv *ds, real factor,
56 const FunDensProp* dens_prop);
57 static void optx_second(FunSecondFuncDrv *ds, real factor,
58 const FunDensProp* dens_prop);
59 static void optx_third(FunThirdFuncDrv *ds, real factor,
60 const FunDensProp* dens_prop);
61
62 Functional OPTXFunctional = {
63 "OPTX_", /* name */
64 optx_isgga, /* gga-corrected */
65 1,
66 optxread, /* set bloody common blocks */
67 NULL, /* reporter */
68 optx_energy,
69 optx_first,
70 optx_second,
71 optx_third
72
73 };
74
75 /* IMPLEMENTATION PART */
76
77 static integer
optxread(const char * conf_line)78 optxread(const char* conf_line)
79 {
80 fun_set_hf_weight(0.0);
81 return 1;
82 }
83
84 /* optx_energy:
85 note that in reality E_OPTX = E_OPTX,alpha + E_OPTX,beta
86 i.e the energy is linear in alpha and beta densities.
87
88 OPTX threshold is needed to avoid numerical problems on 0/0
89 divisions. The problems are small but it is better to be on the
90 safe side.
91 */
92 static const real OPTX_THRESHOLD = 1e-14;
93 static const real GAMMA = 0.006;
94 static real
optx_energy(const FunDensProp * dp)95 optx_energy(const FunDensProp* dp)
96 {
97 real ea,eb;
98 if (dp->rhob<OPTX_THRESHOLD)
99 eb = 0.0;
100 else {
101 real rb43 = pow(dp->rhob,4.0/3.0);
102 real grb = dp->gradb;
103 real xb = grb/rb43;
104 real gxb2 = GAMMA*xb*xb;
105 real ub = gxb2/(1.0 + gxb2);
106 eb = rb43*ub*ub;
107 }
108 if (dp->rhoa<OPTX_THRESHOLD)
109 ea=0;
110 else {
111 real ra43 = pow(dp->rhoa,4.0/3.0);
112 real xa = (dp->grada)/ra43;
113 real gxa2 = GAMMA*xa*xa;
114 real ua = gxa2/(1.0 + gxa2);
115 ea = ra43*ua*ua;
116 }
117 return (ea+eb);
118 }
119
120 static void
optx_first(FunFirstFuncDrv * ds,real factor,const FunDensProp * dp)121 optx_first(FunFirstFuncDrv *ds, real factor, const FunDensProp* dp)
122 {
123 real ra, gra, xa, ra43, ra13, ua, ua2, gxa2;
124 real rb, grb, xb, rb43, rb13, ub, ub2, gxb2;
125 real faR, faZ;
126 real fbR, fbZ;
127
128 if (dp->rhoa >OPTX_THRESHOLD) {
129 ra = dp->rhoa;
130 gra = dp->grada;
131 ra43 = pow(dp->rhoa,4.0/3.0);
132 ra13 = pow(dp->rhoa,1.0/3.0);
133 xa = gra/ra43;
134 gxa2 = GAMMA*xa*xa;
135 ua = gxa2/(1.0 + gxa2);
136 ua2 = ua*ua;
137 faZ = 4.0*ua2*(1.0 - ua)/xa;
138 faR = 4.0*ua2*ra13*(4.0/3.0*ua - 1.0);
139 ds->df1000 += factor*faR;
140 ds->df0010 += factor*faZ;
141 }
142 if (dp->rhob >OPTX_THRESHOLD) {
143 rb = dp->rhob;
144 grb = dp->gradb;
145 rb43 = pow(dp->rhob,4.0/3.0);
146 rb13 = pow(dp->rhob,1.0/3.0);
147 xb = grb/rb43;
148 gxb2 = GAMMA*xb*xb;
149 ub = gxb2/(1.0 + gxb2);
150 ub2 = ub*ub;
151 fbZ = 4.0*ub2*(1.0 - ub)/xb;
152 fbR = 4.0*ub2*rb13*(4.0/3.0*ub - 1.0);
153 ds->df0100 += factor*fbR;
154 ds->df0001 += factor*fbZ;
155 }
156 }
157
158
159 static void
optx_second(FunSecondFuncDrv * ds,real factor,const FunDensProp * dp)160 optx_second(FunSecondFuncDrv *ds, real factor, const FunDensProp* dp)
161 {
162 real ra, gra, xa, ra43, ra13, ua, ua2, gxa2;
163 real rb, grb, xb, rb43, rb13, ub, ub2, gxb2;
164 real faR, faZ, faZZ, faRZ, faRR;
165 real fbR, fbZ, fbZZ, fbRZ, fbRR;
166 real faca, fac2a, facb, fac2b;
167
168 if (dp->rhoa >OPTX_THRESHOLD) {
169 ra = dp->rhoa;
170 ra43 = pow(dp->rhoa,4.0/3.0);
171 ra13 = pow(dp->rhoa,1.0/3.0);
172 gra = dp->grada;
173 xa = (dp->grada)/ra43;
174 gxa2 = GAMMA*xa*xa;
175 ua = gxa2/(1.0 + gxa2);
176 ua2 = ua*ua;
177 faZ = 4.0*ua2*(1.0 - ua)/xa;
178 faR = 4.0/3.0*ra13*ua2*(4.0*ua - 3.0);
179 faca = (1.0 - 2.0*ua)*(1.0 - ua);
180 faZZ = 12.0*ua2/(xa*xa*ra43);
181 faRZ = -16.0*ua2/(xa*ra);
182 fac2a = 4.0*ua2*ra13/(3.0*ra);
183 faRR = 4.0/3.0*ua - 1.0 + 16.0*faca;
184 ds->df1000 += factor*faR;
185 ds->df0010 += factor*faZ;
186 ds->df2000 += factor*faRR*fac2a;
187 ds->df0020 += factor*faZZ*faca;
188 ds->df1010 += factor*faRZ*faca;
189 }
190 if (dp->rhob >OPTX_THRESHOLD) {
191 rb = dp->rhob;
192 rb43 = pow(dp->rhob,4.0/3.0);
193 rb13 = pow(dp->rhob,1.0/3.0);
194 grb = dp->gradb;
195 xb = (dp->gradb)/rb43;
196 gxb2 = GAMMA*xb*xb;
197 ub = gxb2/(1.0 + gxb2);
198 ub2 = ub*ub;
199 fbZ = 4.0*ub2*(1.0 - ub)/xb;
200 fbR = 4.0/3.0*rb13*ub2*(4.0*ub - 3.0);
201 facb = (1.0 - 2.0*ub)*(1.0 - ub);
202 fbZZ = 12.0*ub2/(xb*xb*rb43);
203 fbRZ = -16.0*ub2/(xb*rb);
204 fac2b = 4.0*ub2*rb13/(3.0*rb);
205 fbRR = 4.0/3.0*ub - 1.0 + 16.0*facb;
206 ds->df0100 += factor*fbR;
207 ds->df0001 += factor*fbZ;
208 ds->df0200 += factor*fbRR*fac2b;
209 ds->df0002 += factor*fbZZ*facb;
210 ds->df0101 += factor*fbRZ*facb;
211
212 }
213 }
214
215 static void
optx_third(FunThirdFuncDrv * ds,real factor,const FunDensProp * dp)216 optx_third(FunThirdFuncDrv *ds, real factor, const FunDensProp* dp)
217 {
218 real ra, gra, xa, xa2, ra43, ra13, ua, ua2, gxa2;
219 real rb, grb, xb, xb2, rb43, rb13, ub, ub2, gxb2;
220 real faR, faZ, faZZ, faRZ, faRR;
221 real fbR, fbZ, fbZZ, fbRZ, fbRR;
222 real faRRR, faRRZ, faRZZ, faZZZ;
223 real fbRRR, fbRRZ, fbRZZ, fbZZZ;
224 real faca, fac2a, t3a, t4a, t5a;
225 real facb, fac2b, t3b, t4b, t5b;
226
227 if (dp->rhoa >OPTX_THRESHOLD) {
228 ra = dp->rhoa;
229 ra43 = pow(dp->rhoa,4.0/3.0);
230 ra13 = pow(dp->rhoa,1.0/3.0);
231 gra = dp->grada;
232 xa = (dp->grada)/ra43;
233 xa2 = xa*xa;
234 gxa2 = GAMMA*xa2;
235 ua = gxa2/(1.0 + gxa2);
236 ua2 = ua*ua;
237 faZ = 4.0*ua2*(1.0 - ua)/xa;
238 faR = 4.0/3.0*ra13*ua2*(4.0*ua - 3.0);
239 faca = (1.0 - 2.0*ua)*(1.0 - ua);
240 faZZ = 12.0*ua2/(xa*xa*ra43);
241 faRZ = -16.0*ua2/(xa*ra);
242 fac2a = 4.0*ua2*ra13/(3.0*ra);
243 faRR = 4.0/3.0*ua - 1.0 + 16.0*faca;
244 t3a = 4.0*ua2*(1.0 - ua)/(xa2*ra43);
245 t4a = 16.0*ua2*(1.0 - ua)/(3.0*xa*ra*ra);
246 t5a = -8.0*ua2*ra13/(9.0*ra*ra);
247 faZZZ = 6.0/(xa*ra43)*(8.0*ua2 - 7.0*ua + 1.0);
248 faRZZ = -(4.0/ra)*(4.0*ua - 3.0)*(4.0*ua - 1.0);
249 faRRZ = 64.0*ua2 - 70.0*ua + 15.0;
250 faRRR = 4.0*(1.0-ua)*(128.0*ua2-140.0*ua+30.0) +
251 32.0*ua2 - 140.0/3.0*ua + 15;
252 ds->df1000 += factor*faR;
253 ds->df0010 += factor*faZ;
254 ds->df2000 += factor*faRR*fac2a;
255 ds->df0020 += factor*faZZ*faca;
256 ds->df1010 += factor*faRZ*faca;
257 ds->df0030 += factor*faZZZ*t3a;
258 ds->df1020 += factor*faRZZ*t3a;
259 ds->df2010 += factor*faRRZ*t4a;
260 ds->df3000 += factor*faRRR*t5a;
261 }
262 if (dp->rhob >OPTX_THRESHOLD) {
263 rb = dp->rhob;
264 rb43 = pow(dp->rhob,4.0/3.0);
265 rb13 = pow(dp->rhob,1.0/3.0);
266 grb = dp->gradb;
267 xb = (dp->gradb)/rb43;
268 xb2 = xb*xb;
269 gxb2 = GAMMA*xb2;
270 ub = gxb2/(1.0 + gxb2);
271 ub2 = ub*ub;
272 fbZ = 4.0*ub2*(1.0 - ub)/xb;
273 fbR = 4.0/3.0*rb13*ub2*(4.0*ub - 3.0);
274 facb = (1.0 - 2.0*ub)*(1.0 - ub);
275 fbZZ = 12.0*ub2/(xb*xb*rb43);
276 fbRZ = -16.0*ub2/(xb*rb);
277 fac2b = 4.0*ub2*rb13/(3.0*rb);
278 fbRR = 4.0/3.0*ub - 1.0 + 16.0*facb;
279 t3b = 4.0*ub2*(1.0 - ub)/(xb2*rb43);
280 t4b = 16.0*ub2*(1.0 - ub)/(3.0*xb*rb*rb);
281 t5b = -8.0*ub2*rb13/(9.0*rb*rb);
282 fbZZZ = 6.0/(xb*rb43)*(8.0*ub2 - 7.0*ub + 1.0);
283 fbRZZ = -(4.0/rb)*(4.0*ub - 3.0)*(4.0*ub - 1.0);
284 fbRRZ = 64.0*ub2 - 70.0*ub + 15.0;
285 fbRRR = 4.0*(1.0-ub)*(128.0*ub2-140.0*ub+30.0) +
286 32.0*ub2 - 140.0/3.0*ub + 15;
287 ds->df0100 += factor*fbR;
288 ds->df0001 += factor*fbZ;
289 ds->df0200 += factor*fbRR*fac2b;
290 ds->df0002 += factor*fbZZ*facb;
291 ds->df0101 += factor*fbRZ*facb;
292 ds->df0003 += factor*fbZZZ*t3b;
293 ds->df0102 += factor*fbRZZ*t3b;
294 ds->df0201 += factor*fbRRZ*t4b;
295 ds->df0300 += factor*fbRRR*t5b;
296
297 }
298 }
299
300
301