MADNESS 0.10.1
lda.h
Go to the documentation of this file.
1/*
2 This file is part of MADNESS.
3
4 Copyright (C) 2007,2010 Oak Ridge National Laboratory
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 2 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, write to the Free Software
18 Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
19
20 For more information please contact:
21
22 Robert J. Harrison
23 Oak Ridge National Laboratory
24 One Bethel Valley Road
25 P.O. Box 2008, MS-6367
26
27 email: harrisonrj@ornl.gov
28 tel: 865-241-3937
29 fax: 865-572-0680
30*/
31
32/*
33 * lda.h
34 *
35 * Created on: Nov 10, 2008
36 * Author: wsttiger
37 */
38
39#ifndef LDA_H_
40#define LDA_H_
41
42#include <madness/mra/mra.h>
44#include <math.h>
46
47typedef double doublereal;
48typedef MADNESS_FORINT integer;
49typedef int logical;
50
51#define max(a,b) ((a)>=(b)?(a):(b))
52
53/* static double pow_dd(doublereal* a, doublereal* b) { */
54/* return pow(*a,*b); */
55/* } */
56
57/* lda.f -- translated by f2c (version 20050501).
58 You must link the resulting object file with libf2c:
59 on Microsoft Windows system, link with libf2c.lib;
60 on Linux or Unix systems, link with .../path/to/libf2c.a -lm
61 or, if you install libf2c.a in a standard place, with -lf2c -lm
62 -- in that order, at the end of the command line, as in
63 cc *.o -lf2c -lm
64 Source for libf2c is in /netlib/f2c/libf2c.zip, e.g.,
65
66 http://www.netlib.org/f2c/libf2c.zip
67*/
68
69
70/* Table of constant values */
71
72#include <cmath>
73
74static inline double pow(const double* a, const double* b) {
75 return pow(*a, *b);
76}
77
78static double c_b2 = .333333333333333333333333333333333;
79static double c_b7 = .333333333333333333333333333333;
80static double c_b8 = .5;
81static double c_b14 = 1.333333333333333333333333333333;
82
83
84inline /* Subroutine */ int x_rks_s__(const double *r__, double *f, double *
85 dfdra)
86{
87
88 /* Local variables */
89 double ra13;
90
91
92/* This subroutine evaluates the spin polarised exchange functional */
93/* in the Local Density Approximation [1], and the corresponding */
94/* potential. Often this functional is referred to as the Dirac */
95/* functional [2] or Slater functional. */
96
97/* [1] F. Bloch, Zeitschrift fuer Physik, Vol. 57 (1929) 545. */
98
99/* [2] P.A.M. Dirac, Proceedings of the Cambridge Philosophical */
100/* Society, Vol. 26 (1930) 376. */
101
102/* Parameters: */
103
104/* r the total electron density */
105/* f On return the functional value */
106/* dfdra On return the derivative of f with respect to alpha electron */
107/* density */
108
109
110/* Ax = -3/4*(6/pi)**(1/3) */
111/* Bx = -(6/pi)**(1/3) */
112/* C = (1/2)**(1/3) */
113
114
115
116
117 ra13 = pow(r__, &c_b2) * .793700525984099737375852819636154;
118 *f = *r__ * -.930525736349100025002010218071667 * ra13;
119 *dfdra = ra13 * -1.24070098179880003333601362409556;
120
121 return 0;
122} /* x_rks_s__ */
123
124/* ----------------------------------------------------------------------- */
125inline /* Subroutine */ int x_uks_s__(double *ra, double *rb, double *f,
126 double *dfdra, double *dfdrb)
127{
128 /* Local variables */
129 double ra13, rb13;
130
131
132/* This subroutine evaluates the spin polarised exchange functional */
133/* in the Local Density Approximation [1], and the corresponding */
134/* potential. Often this functional is referred to as the Dirac */
135/* functional [2] or Slater functional. */
136
137/* [1] F. Bloch, Zeitschrift fuer Physik, Vol. 57 (1929) 545. */
138
139/* [2] P.A.M. Dirac, Proceedings of the Cambridge Philosophical */
140/* Society, Vol. 26 (1930) 376. */
141
142/* Parameters: */
143
144/* ra the alpha electron density */
145/* rb the beta electron density */
146/* f On return the functional value */
147/* dfdra On return the derivative of f with respect to ra */
148/* dfdrb On return the derivative of f with respect to rb */
149
150
151/* Ax = -3/4*(6/pi)**(1/3) */
152/* Bx = -(6/pi)**(1/3) */
153
154
155
156
157 ra13 = pow(ra, &c_b2);
158 rb13 = pow(rb, &c_b2);
159 *f = (*ra * ra13 + *rb * rb13) * -.930525736349100025002010218071667;
160 *dfdra = ra13 * -1.24070098179880003333601362409556;
161 *dfdrb = rb13 * -1.24070098179880003333601362409556;
162
163 return 0;
164} /* x_uks_s__ */
165
166inline /* Subroutine */ int c_rks_vwn5__(const double *r__, double *f, double *
167 dfdra)
168{
169 /* Local variables */
170 double a2, b2, c2, d2, i1, i2, i3, p1, p2, p3, p4, t4, t5, t6,
171 t7, iv, alpha_rho13__, iv2, pp1, pp2, inv, srho, srho13;
172
173
174/* This subroutine evaluates the Vosko, Wilk and Nusair correlation */
175/* functional number 5 [1] for the closed shell case, with the */
176/* parametrisation as given in table 5. */
177
178/* The original code was obtained from Dr. Phillip Young, */
179/* with corrections from Dr. Paul Sherwood. */
180
181/* [1] S.H. Vosko, L. Wilk, and M. Nusair */
182/* "Accurate spin-dependent electron liquid correlation energies */
183/* for local spin density calculations: a critical analysis", */
184/* Can.J.Phys, Vol. 58 (1980) 1200-1211. */
185
186/* Parameters: */
187
188/* r the total electron density */
189/* f On return the functional value */
190/* dfdra On return the derivative of f with respect to the alpha */
191/* electron density */
192
193
194
195
196/* VWN interpolation parameters */
197
198/* paramagnetic */
199 a2 = .0621814;
200 b2 = 3.72744;
201 c2 = 12.9352;
202 d2 = -.10498;
203
204/* t4 = (1/(4/3)*pi)**(1/3) */
205 t4 = .620350490899399531;
206
207/* t5 = 0.5/(2**(1/3)-1) */
208 t5 = 1.92366105093153617;
209
210/* t6 = 2/(3*(2**(1/3)-1)) */
211 t6 = 2.56488140124204822;
212
213/* t7 = 2.25*(2**(1/3)-1) */
214 t7 = .584822362263464735;
215
216/* Paramagnetic interpolation constants */
217
218 p1 = 6.1519908197590798;
219 p2 = a2 * .5;
220 p3 = 9.6902277115443745e-4;
221 p4 = .038783294878113009;
222
223/* closed shell case */
224 srho = *r__;
225 srho13 = pow(&srho, &c_b7);
226 alpha_rho13__ = pow(&c_b8, &c_b7) * srho;
227 iv2 = t4 / srho13;
228 iv = sqrt(iv2);
229
230/* paramagnetic */
231 inv = 1. / (iv2 + b2 * iv + c2);
232 i1 = log(iv2 * inv);
233 i2 = log((iv - d2) * (iv - d2) * inv);
234/* corrected b1->b2 ps Apr98 */
235 i3 = atan(p1 / (iv * 2. + b2));
236 pp1 = p2 * i1 + p3 * i2 + p4 * i3;
237 pp2 = a2 * (1. / iv - iv * inv * (b2 / (iv - d2) + 1.));
238
239 *f = pp1 * srho;
240 *dfdra = pp1 - iv * .166666666666666666666666666666 * pp2;
241
242 return 0;
243} /* c_rks_vwn5__ */
244
245/* ----------------------------------------------------------------------- */
246inline /* Subroutine */ int c_uks_vwn5__(double *ra, double *rb, double *
247 f, double *dfdra, double *dfdrb)
248{
249 /* System generated locals */
250 double d__1, d__2;
251
252 /* Local variables */
253 double v, beta_rho13__, a1, b1, c1, d1, a2, b2, c2, d2, a3, b3,
254 c3, d3, f1, f2, f3, p1, p2, p3, s1, t1, t2, s2, t4, t5, t6, t7,
255 s3, s4, p4, f4, i1, i2, i3, iv, alpha_rho13__, ff1, ff2, iv2, pp1,
256 pp2, ss1, ss2, tau, inv, vwn1, vwn2, dtau, zeta, srho, zeta3,
257 zeta4, srho13, inter1, inter2;
258
259
260/* This subroutine evaluates the Vosko, Wilk and Nusair correlation */
261/* functional number 5 [1], with the parametrisation as given in */
262/* table 5. */
263
264/* The original code was obtained from Dr. Phillip Young, */
265/* with corrections from Dr. Paul Sherwood. */
266
267/* [1] S.H. Vosko, L. Wilk, and M. Nusair */
268/* "Accurate spin-dependent electron liquid correlation energies */
269/* for local spin density calculations: a critical analysis", */
270/* Can.J.Phys, Vol. 58 (1980) 1200-1211. */
271
272/* Parameters: */
273
274/* ra the alpha-electron density */
275/* rb the beta-electron density */
276/* f On return the functional value */
277/* dfdra On return the derivative of f with respect to ra */
278/* dfdrb On return the derivative of f with respect to rb */
279
280
281
282/* tn13 = 2**(1/3) */
283/* tn43 = 2**(4/3) */
284
285/* VWN interpolation parameters */
286
287/* spin stiffness */
288 a1 = -.0337737278807791058;
289 b1 = 1.13107;
290 c1 = 13.0045;
291 d1 = -.0047584;
292/* paramagnetic */
293 a2 = .0621814;
294 b2 = 3.72744;
295 c2 = 12.9352;
296 d2 = -.10498;
297/* ferromagnetic */
298/* try cadpac/nwchem value (.5*a2) */
299 a3 = .0310907;
300 b3 = 7.06042;
301 c3 = 18.0578;
302 d3 = -.325;
303
304/* t4 = (1/(4/3)*pi)**(1/3) */
305 t4 = .620350490899399531;
306
307/* t5 = 0.5/(2**(1/3)-1) */
308 t5 = 1.92366105093153617;
309
310/* t6 = 2/(3*(2**(1/3)-1)) */
311 t6 = 2.56488140124204822;
312
313/* t7 = 2.25*(2**(1/3)-1) */
314 t7 = .584822362263464735;
315
316/* Spin stiffness interpolation constants */
317
318 s1 = 7.12310891781811772;
319 s2 = a1 * .5;
320 s3 = -6.9917323507644313e-6;
321 s4 = -.0053650918488835769;
322
323/* Paramagnetic interpolation constants */
324
325 p1 = 6.1519908197590798;
326 p2 = a2 * .5;
327 p3 = 9.6902277115443745e-4;
328 p4 = .038783294878113009;
329
330/* Ferromagnetic interpolation constants */
331
332 f1 = 4.73092690956011364;
333 f2 = a3 * .5;
334
335/* F3 = -0.244185082989490298d-02 *0.5d0 */
336/* F4 = -0.570212323620622186d-01 *0.5d0 */
337
338/* try nwchem values */
339
340 f3 = .00224786709554261133;
341 f4 = .0524913931697809227;
342
343/* Interpolation intervals */
344
345 inter1 = .99999999989999999;
346 inter2 = -.99999999989999999;
347
348/* open shell case */
349 alpha_rho13__ = pow(ra, &c_b7);
350 beta_rho13__ = pow(rb, &c_b7);
351 srho = *ra + *rb;
352 srho13 = pow(&srho, &c_b7);
353 iv2 = t4 / srho13;
354 iv = sqrt(iv2);
355
356/* spin-stiffness */
357 inv = 1. / (iv2 + b1 * iv + c1);
358 i1 = log(iv2 * inv);
359 i2 = log((iv - d1) * (iv - d1) * inv);
360 i3 = atan(s1 / (iv * 2. + b1));
361 ss1 = s2 * i1 + s3 * i2 + s4 * i3;
362 ss2 = a1 * (1. / iv - iv * inv * (b1 / (iv - d1) + 1.));
363
364/* paramagnetic */
365 inv = 1. / (iv2 + b2 * iv + c2);
366 i1 = log(iv2 * inv);
367 i2 = log((iv - d2) * (iv - d2) * inv);
368/* corrected b1->b2 ps Apr98 */
369 i3 = atan(p1 / (iv * 2. + b2));
370 pp1 = p2 * i1 + p3 * i2 + p4 * i3;
371 pp2 = a2 * (1. / iv - iv * inv * (b2 / (iv - d2) + 1.));
372
373/* ferromagnetic */
374 inv = 1. / (iv2 + b3 * iv + c3);
375 i1 = log(iv2 * inv);
376 i2 = log((iv - d3) * (iv - d3) * inv);
377 i3 = atan(f1 / (iv * 2. + b3));
378 ff1 = f2 * i1 + f3 * i2 + f4 * i3;
379 ff2 = a3 * (1. / iv - iv * inv * (b3 / (iv - d3) + 1.));
380
381/* polarisation function */
382
383 zeta = (*ra - *rb) / srho;
384 zeta3 = zeta * zeta * zeta;
385 zeta4 = zeta3 * zeta;
386 if (zeta > inter1) {
387 vwn1 = t5 * .51984209978974638;
388 vwn2 = t6 * 1.25992104989487316476721060727823;
389 } else if (zeta < inter2) {
390 vwn1 = t5 * .51984209978974638;
391 vwn2 = t6 * -1.25992104989487316476721060727823;
392 } else {
393 d__1 = zeta + 1.;
394 d__2 = 1. - zeta;
395 vwn1 = (pow(&d__1, &c_b14) + pow(&d__2, &c_b14) - 2.) * t5;
396 d__1 = zeta + 1.;
397 d__2 = 1. - zeta;
398 vwn2 = (pow(&d__1, &c_b7) - pow(&d__2, &c_b7)) * t6;
399 }
400 ss1 *= t7;
401 ss2 *= t7;
402 tau = ff1 - pp1 - ss1;
403 dtau = ff2 - pp2 - ss2;
404
405 v = pp1 + vwn1 * (ss1 + tau * zeta4);
406 *f = v * srho;
407
408 t1 = v - iv * .166666666666666666666666666667 * (pp2 + vwn1 * (ss2 + dtau
409 * zeta4));
410 t2 = vwn2 * (ss1 + tau * zeta4) + vwn1 * 4. * tau * zeta3;
411 *dfdra = t1 + t2 * (1. - zeta);
412 *dfdrb = t1 - t2 * (zeta + 1.);
413
414 return 0;
415} /* c_uks_vwn5__ */
416
417////***************************************************************************
418//int rks_x_lda__ (integer *ideriv, integer *npt, doublereal *rhoa1,
419// doublereal *sigmaaa1, doublereal *zk, doublereal *vrhoa,
420// doublereal *vsigmaaa, doublereal *v2rhoa2, doublereal *v2rhoasigmaaa,
421// doublereal *v2sigmaaa2);
422////***************************************************************************
423//
424////***************************************************************************
425//int rks_c_vwn5__ (integer *ideriv, integer *npt, doublereal *rhoa1,
426// doublereal *sigmaaa1, doublereal *zk, doublereal *vrhoa,
427// doublereal *vsigmaaa, doublereal *v2rhoa2, doublereal *v2rhoasigmaaa,
428// doublereal *v2sigmaaa2);
429////***************************************************************************
430
431inline /* Subroutine */ int rks_c_vwn5__(integer *ideriv, integer *npt, doublereal *
432 rhoa1, doublereal *sigmaaa1, doublereal *zk, doublereal *vrhoa,
433 doublereal *vsigmaaa, doublereal *v2rhoa2, doublereal *v2rhoasigmaaa,
434 doublereal *v2sigmaaa2)
435{
436 // WSTHORNTON
437 static doublereal c_b2 = .16666666666666666;
438 static doublereal c_b3 = .33333333333333331;
439 /* System generated locals */
440 integer i__1;
441 doublereal d__1, d__2;
442
443 /* Builtin functions */
444 double pow_dd(doublereal *, doublereal *), log(doublereal), atan(
445 doublereal);
446
447 /* Local variables */
448 static integer i__;
449 static doublereal s1, t1, t2, t4, t6, t7, t10, t20, t11, t22, t13, t23,
450 t16, t17, t25, t19, t26, t28, t29, t32, t33, t37, t38, t40, t41,
451 t43, t51, t52, t27, t34, t39, t46, t47, t48, t49, t53, t55, t56,
452 t58, t60, t63, t66, t67, t68, t69, t77, t79, t80, t81, t88, t92,
453 t94, t102, t103, t105, t107, t125, t134, t138, rho;
454
455
456/* S.H. Vosko, L. Wilk, and M. Nusair */
457/* Accurate spin-dependent electron liquid correlation energies for */
458/* local spin density calculations: a critical analysis */
459/* Can. J. Phys. 58 (1980) 1200-1211 */
460
461
462/* CITATION: */
463
464/* Functionals were obtained from the Density Functional Repository */
465/* as developed and distributed by the Quantum Chemistry Group, */
466/* CCLRC Daresbury Laboratory, Daresbury, Cheshire, WA4 4AD */
467/* United Kingdom. Contact Huub van Dam (h.j.j.vandam@dl.ac.uk) or */
468/* Paul Sherwood for further information. */
469
470/* COPYRIGHT: */
471
472/* Users may incorporate the source code into software packages and */
473/* redistribute the source code provided the source code is not */
474/* changed in anyway and is properly cited in any documentation or */
475/* publication related to its use. */
476
477/* ACKNOWLEDGEMENT: */
478
479/* The source code was generated using Maple 8 through a modified */
480/* version of the dfauto script published in: */
481
482/* R. Strange, F.R. Manby, P.J. Knowles */
483/* Automatic code generation in density functional theory */
484/* Comp. Phys. Comm. 136 (2001) 310-318. */
485
486 /* Parameter adjustments */
487 --v2sigmaaa2;
488 --v2rhoasigmaaa;
489 --v2rhoa2;
490 --vsigmaaa;
491 --vrhoa;
492 --zk;
493 --sigmaaa1;
494 --rhoa1;
495
496 /* Function Body */
497 if (*ideriv == 0) {
498 i__1 = *npt;
499 for (i__ = 1; i__ <= i__1; ++i__) {
500/* Computing MAX */
501 d__1 = 0., d__2 = rhoa1[i__];
502 rho = max(d__1,d__2);
503 if (rho > 1e-20) {
504 t1 = 1 / rho;
505 t2 = pow_dd(&t1, &c_b3);
506 t4 = pow_dd(&t1, &c_b2);
507 t7 = 1 / (t2 * .6203504908994 + t4 * 2.935818660072219 +
508 12.9352);
509 t10 = log(t2 * .6203504908994 * t7);
510 t16 = atan(6.15199081975908 / (t4 * 1.575246635799487 +
511 3.72744));
512/* Computing 2nd power */
513 d__1 = t4 * .7876233178997433 + .10498;
514 t20 = d__1 * d__1;
515 t22 = log(t20 * t7);
516 zk[i__] = rho * (t10 * .0310907 + t16 * .03878329487811301 +
517 t22 * 9.690227711544374e-4);
518 } else {
519/* rho */
520 zk[i__] = 0.;
521 }
522/* rho */
523 }
524 } else if (*ideriv == 1) {
525 i__1 = *npt;
526 for (i__ = 1; i__ <= i__1; ++i__) {
527/* Computing MAX */
528 d__1 = 0., d__2 = rhoa1[i__];
529 rho = max(d__1,d__2);
530 if (rho > 1e-20) {
531 t1 = 1 / rho;
532 t2 = pow_dd(&t1, &c_b3);
533 t4 = pow_dd(&t1, &c_b2);
534 t6 = t2 * .6203504908994 + t4 * 2.935818660072219 + 12.9352;
535 t7 = 1 / t6;
536 t10 = log(t2 * .6203504908994 * t7);
537 t11 = t10 * .0310907;
538 t13 = t4 * 1.575246635799487 + 3.72744;
539 t16 = atan(6.15199081975908 / t13);
540 t17 = t16 * .03878329487811301;
541 t19 = t4 * .7876233178997433 + .10498;
542/* Computing 2nd power */
543 d__1 = t19;
544 t20 = d__1 * d__1;
545 t22 = log(t20 * t7);
546 t23 = t22 * 9.690227711544374e-4;
547 zk[i__] = rho * (t11 + t17 + t23);
548/* Computing 2nd power */
549 d__1 = t2;
550 t25 = d__1 * d__1;
551 t26 = 1 / t25;
552/* Computing 2nd power */
553 d__1 = rho;
554 t28 = d__1 * d__1;
555 t29 = 1 / t28;
556/* Computing 2nd power */
557 d__1 = t6;
558 t32 = d__1 * d__1;
559 t33 = 1 / t32;
560/* Computing 2nd power */
561 d__1 = t4;
562 t37 = d__1 * d__1;
563/* Computing 2nd power */
564 d__1 = t37;
565 t38 = d__1 * d__1;
566 t40 = 1 / t38 / t4;
567 t41 = t40 * t29;
568 t43 = t26 * -.2067834969664667 * t29 - t41 *
569 .4893031100120365;
570/* Computing 2nd power */
571 d__1 = t13;
572 t51 = d__1 * d__1;
573 t52 = 1 / t51;
574 vrhoa[i__] = t11 + t17 + t23 + rho * ((t26 *
575 -.2067834969664667 * t7 * t29 - t2 * .6203504908994 *
576 t33 * t43) * .05011795824473985 / t2 * t6 + t52 *
577 .0626408570946439 * t40 * t29 / (t52 * 37.8469910464
578 + 1.) + (t19 * -.2625411059665811 * t7 * t41 - t20 *
579 1. * t33 * t43) * 9.690227711544374e-4 / t20 * t6);
580 } else {
581/* rho */
582 zk[i__] = 0.;
583 vrhoa[i__] = 0.;
584 }
585/* rho */
586 }
587 } else if (*ideriv == 2) {
588 i__1 = *npt;
589 for (i__ = 1; i__ <= i__1; ++i__) {
590/* Computing MAX */
591 d__1 = 0., d__2 = rhoa1[i__];
592 rho = max(d__1,d__2);
593 if (rho > 1e-20) {
594 t1 = 1 / rho;
595 t2 = pow_dd(&t1, &c_b3);
596 t4 = pow_dd(&t1, &c_b2);
597 t6 = t2 * .6203504908994 + t4 * 2.935818660072219 + 12.9352;
598 t7 = 1 / t6;
599 t10 = log(t2 * .6203504908994 * t7);
600 t11 = t10 * .0310907;
601 t13 = t4 * 1.575246635799487 + 3.72744;
602 t16 = atan(6.15199081975908 / t13);
603 t17 = t16 * .03878329487811301;
604 t19 = t4 * .7876233178997433 + .10498;
605/* Computing 2nd power */
606 d__1 = t19;
607 t20 = d__1 * d__1;
608 t22 = log(t20 * t7);
609 t23 = t22 * 9.690227711544374e-4;
610 zk[i__] = rho * (t11 + t17 + t23);
611/* Computing 2nd power */
612 d__1 = t2;
613 t25 = d__1 * d__1;
614 t26 = 1 / t25;
615 t27 = t26 * t7;
616/* Computing 2nd power */
617 d__1 = rho;
618 t28 = d__1 * d__1;
619 t29 = 1 / t28;
620/* Computing 2nd power */
621 d__1 = t6;
622 t32 = d__1 * d__1;
623 t33 = 1 / t32;
624 t34 = t2 * t33;
625/* Computing 2nd power */
626 d__1 = t4;
627 t37 = d__1 * d__1;
628/* Computing 2nd power */
629 d__1 = t37;
630 t38 = d__1 * d__1;
631 t39 = t38 * t4;
632 t40 = 1 / t39;
633 t41 = t40 * t29;
634 t43 = t26 * -.2067834969664667 * t29 - t41 *
635 .4893031100120365;
636 t46 = t27 * -.2067834969664667 * t29 - t34 * .6203504908994 *
637 t43;
638 t47 = 1 / t2;
639 t48 = t46 * t47;
640 t49 = t48 * t6;
641/* Computing 2nd power */
642 d__1 = t13;
643 t51 = d__1 * d__1;
644 t52 = 1 / t51;
645 t53 = t52 * t40;
646 t55 = t52 * 37.8469910464 + 1.;
647 t56 = 1 / t55;
648 t58 = t53 * t29 * t56;
649 t60 = t19 * t7;
650 t63 = t20 * t33;
651 t66 = t60 * -.2625411059665811 * t41 - t63 * 1. * t43;
652 t67 = 1 / t20;
653 t68 = t66 * t67;
654 t69 = t68 * t6;
655 vrhoa[i__] = t11 + t17 + t23 + rho * (t49 *
656 .05011795824473985 + t58 * .0626408570946439 + t69 *
657 9.690227711544374e-4);
658 t77 = 1 / t25 / t1;
659/* Computing 2nd power */
660 d__1 = t28;
661 t79 = d__1 * d__1;
662 t80 = 1 / t79;
663 t81 = t77 * t7 * t80;
664 t88 = 1 / t28 / rho;
665 t92 = 1 / t32 / t6;
666/* Computing 2nd power */
667 d__1 = t43;
668 t94 = d__1 * d__1;
669 t102 = 1 / t39 / t1;
670 t103 = t102 * t80;
671 t105 = t40 * t88;
672 t107 = t77 * -.1378556646443111 * t80 + t26 *
673 .4135669939329333 * t88 - t103 * .4077525916766971 +
674 t105 * .978606220024073;
675 t125 = t80 * t56;
676/* Computing 2nd power */
677 d__1 = t51;
678 t134 = d__1 * d__1;
679/* Computing 2nd power */
680 d__1 = t55;
681 t138 = d__1 * d__1;
682 s1 = t49 * .2004718329789594 + t58 * .2505634283785756;
683 v2rhoa2[i__] = s1 + t69 * .00387609108461775 + rho * 2. * ((
684 t81 * -.1378556646443111 + t26 * .4135669939329333 *
685 t33 * t29 * t43 + t27 * .4135669939329333 * t88 + t2 *
686 1.2407009817988 * t92 * t94 - t34 * .6203504908994 *
687 t107) * .05011795824473985 * t47 * t6 + t46 *
688 .01670598608157995 / t2 / t1 * t6 * t29 + t48 *
689 .05011795824473985 * t43 + .03289159980064473 / t51 /
690 t13 * t77 * t125 + t52 * .05220071424553658 * t102 *
691 t125 - t53 * .1252817141892878 * t88 * t56 -
692 1.244848083156773 / t134 / t13 * t77 * t80 / t138 + (
693 t81 * .03446391616107778 + t19 * .5250822119331622 *
694 t33 * t41 * t43 - t60 * .2187842549721509 * t103 +
695 t60 * .5250822119331622 * t105 + t20 * 2. * t92 * t94
696 - t63 * 1. * t107) * 9.690227711544374e-4 * t67 * t6
697 + t66 * 2.544083100456872e-4 / t20 / t19 * t6 * t40 *
698 t29 + t68 * 9.690227711544374e-4 * t43);
699 } else {
700/* rho */
701 zk[i__] = 0.;
702 vrhoa[i__] = 0.;
703 v2rhoa2[i__] = 0.;
704 }
705/* rho */
706 }
707 }
708/* ideriv */
709 return 0;
710} /* rks_c_vwn5__ */
711
712inline /* Subroutine */ int rks_x_lda__(integer *ideriv, integer *npt, doublereal *
713 rhoa1, doublereal *sigmaaa1, doublereal *zk, doublereal *vrhoa,
714 doublereal *vsigmaaa, doublereal *v2rhoa2, doublereal *v2rhoasigmaaa,
715 doublereal *v2sigmaaa2)
716{
717 // WSTHORNTON
718 static doublereal c_b2 = .33333333333333331;
719
720 /* System generated locals */
721 integer i__1;
722 doublereal d__1, d__2;
723
724 /* Builtin functions */
725 double pow_dd(doublereal *, doublereal *);
726
727 /* Local variables */
728 static integer i__;
729 static doublereal t1, t5, rho;
730
731
732/* P.A.M. Dirac */
733/* Proceedings of the Cambridge Philosophical Society, 26 (1930) 376 */
734
735
736/* CITATION: */
737
738/* Functionals were obtained from the Density Functional Repository */
739/* as developed and distributed by the Quantum Chemistry Group, */
740/* CCLRC Daresbury Laboratory, Daresbury, Cheshire, WA4 4AD */
741/* United Kingdom. Contact Huub van Dam (h.j.j.vandam@dl.ac.uk) or */
742/* Paul Sherwood for further information. */
743
744/* COPYRIGHT: */
745
746/* Users may incorporate the source code into software packages and */
747/* redistribute the source code provided the source code is not */
748/* changed in anyway and is properly cited in any documentation or */
749/* publication related to its use. */
750
751/* ACKNOWLEDGEMENT: */
752
753/* The source code was generated using Maple 8 through a modified */
754/* version of the dfauto script published in: */
755
756/* R. Strange, F.R. Manby, P.J. Knowles */
757/* Automatic code generation in density functional theory */
758/* Comp. Phys. Comm. 136 (2001) 310-318. */
759
760 /* Parameter adjustments */
761 --v2sigmaaa2;
762 --v2rhoasigmaaa;
763 --v2rhoa2;
764 --vsigmaaa;
765 --vrhoa;
766 --zk;
767 --sigmaaa1;
768 --rhoa1;
769
770 /* Function Body */
771 if (*ideriv == 0) {
772 i__1 = *npt;
773 for (i__ = 1; i__ <= i__1; ++i__) {
774/* Computing MAX */
775 d__1 = 0., d__2 = rhoa1[i__];
776 rho = max(d__1,d__2);
777 if (rho > 1e-20) {
778 t1 = pow_dd(&rho, &c_b2);
779 zk[i__] = t1 * -.7385587663820224 * rho;
780 } else {
781/* rho */
782 zk[i__] = 0.;
783 }
784/* rho */
785 }
786 } else if (*ideriv == 1) {
787 i__1 = *npt;
788 for (i__ = 1; i__ <= i__1; ++i__) {
789/* Computing MAX */
790 d__1 = 0., d__2 = rhoa1[i__];
791 rho = max(d__1,d__2);
792 if (rho > 1e-20) {
793 t1 = pow_dd(&rho, &c_b2);
794 zk[i__] = t1 * -.7385587663820224 * rho;
795 vrhoa[i__] = t1 * -.9847450218426965;
796 } else {
797/* rho */
798 zk[i__] = 0.;
799 vrhoa[i__] = 0.;
800 }
801/* rho */
802 }
803 } else if (*ideriv == 2) {
804 i__1 = *npt;
805 for (i__ = 1; i__ <= i__1; ++i__) {
806/* Computing MAX */
807 d__1 = 0., d__2 = rhoa1[i__];
808 rho = max(d__1,d__2);
809 if (rho > 1e-20) {
810 t1 = pow_dd(&rho, &c_b2);
811 zk[i__] = t1 * -.7385587663820224 * rho;
812 vrhoa[i__] = t1 * -.9847450218426965;
813/* Computing 2nd power */
814 d__1 = t1;
815 t5 = d__1 * d__1;
816 v2rhoa2[i__] = -.6564966812284644 / t5;
817 } else {
818/* rho */
819 zk[i__] = 0.;
820 vrhoa[i__] = 0.;
821 v2rhoa2[i__] = 0.;
822 }
823/* rho */
824 }
825 }
826/* ideriv */
827 return 0;
828} /* rks_x_lda__ */
829
830const double THRESH_RHO = 1e-8;
831const double THRESH_GRHO = 1e-20;
832
833inline void wst_munge_grho(int npoint, double *rho, double *grho) {
834 for (int i=0; i<npoint; i++) {
835 if (rho[i]<THRESH_RHO) rho[i] = THRESH_RHO;
836 if ((rho[i] <=THRESH_RHO) ||
837 (grho[i] < THRESH_GRHO)) grho[i] = THRESH_GRHO;
838 }
839}
840
841inline void wst_munge_rho(int npoint, double *rho) {
842 for (int i=0; i<npoint; i++) {
843 if (rho[i]<THRESH_RHO) rho[i] = THRESH_RHO;
844 }
845}
846
847//***************************************************************************
848inline void xc_rks_generic_lda(Tensor<double> rho_alpha, ///< Alpha-spin density at each grid point
849 Tensor<double> f, ///< Value of functional at each grid point
850 Tensor<double> df_drho) ///< Derivative of functional w.r.t. rho_alpha
851 {
852 MADNESS_ASSERT(rho_alpha.iscontiguous());
853 MADNESS_ASSERT(f.iscontiguous());
854 MADNESS_ASSERT(df_drho.iscontiguous());
855
856 rho_alpha = rho_alpha.flat();
857 f = f.flat();
858 df_drho = df_drho.flat();
859
860 integer ideriv = 2;
861 integer npt = rho_alpha.dim(0);
862
863 Tensor<double> gamma_alpha(npt);
864 Tensor<double> tf(npt);
865 Tensor<double> tdf_drho(npt);
866 Tensor<double> tdf_dgamma(npt);
867 Tensor<double> td2f_drho2(npt);
868 Tensor<double> td2f_drhodgamma(npt);
869 Tensor<double> td2f_dgamma2(npt);
870
871 wst_munge_rho(npt, rho_alpha.ptr());
872
873 f.fill(0.0);
874 df_drho.fill(0.0);
875
876 int returnvalue = ::rks_x_lda__(&ideriv, &npt, rho_alpha.ptr(), gamma_alpha.ptr(),
877 tf.ptr(),
878 tdf_drho.ptr(), tdf_dgamma.ptr(),
879 td2f_drho2.ptr(), td2f_drhodgamma.ptr(), td2f_dgamma2.ptr());
880
881 f.gaxpy(1.0, tf, 1.0);
882 df_drho.gaxpy(1.0, tdf_drho, 1.0);
883
884 tf.fill(0.0);
885 tdf_drho.fill(0.0);
886
887 returnvalue = ::rks_c_vwn5__(&ideriv, &npt, rho_alpha.ptr(), gamma_alpha.ptr(),
888 tf.ptr(),
889 tdf_drho.ptr(), tdf_dgamma.ptr(),
890 td2f_drho2.ptr(), td2f_drhodgamma.ptr(), td2f_dgamma2.ptr());
891
892 f.gaxpy(1.0, tf, 1.0);
893 df_drho.gaxpy(1.0, tdf_drho, 1.0);
894
895 }
896 //***************************************************************************
897
898 //***************************************************************************
899 template <int NDIM>
900 inline void dft_xc_lda_V(const Key<NDIM>& key, Tensor<double>& t)
901 {
902 Tensor<double> enefunc = copy(t);
903 Tensor<double> V = copy(t);
904 ::xc_rks_generic_lda(t, enefunc, V);
905 t(___) = V(___);
906 }
907 //***************************************************************************
908
909 //***************************************************************************
910 template <int NDIM>
911 inline void dft_xc_lda_ene(const Key<NDIM>& key, Tensor<double>& t)
912 {
913 Tensor<double> V = copy(t);
914 Tensor<double> enefunc = copy(t);
915 ::xc_rks_generic_lda(t, enefunc, V);
916 t(___) = enefunc(___);
917 }
918 //***************************************************************************
919
920// //***************************************************************************
921// static double munge(double r) {
922// if (r < 1e-12) r = 1e-12;
923// return r;
924// }
925// //***************************************************************************
926//
927// //***************************************************************************
928// static void ldaop(const Key<3>& key, Tensor<double>& t) {
929// UNARY_OPTIMIZED_ITERATOR(double, t, double r=munge(2.0* *_p0); double q; double dq1; double dq2;x_rks_s__(&r, &q, &dq1);c_rks_vwn5__(&r, &q, &dq2); *_p0 = dq1+dq2);
930// }
931// //***************************************************************************
932//
933// //***************************************************************************
934// static void ldaeop(const Key<3>& key, Tensor<double>& t) {
935// UNARY_OPTIMIZED_ITERATOR(double, t, double r=munge(2.0* *_p0); double q1; double q2; double dq;x_rks_s__(&r, &q1, &dq);c_rks_vwn5__(&r, &q2, &dq); *_p0 = q1+q2);
936// }
937// //***************************************************************************
938
939
940
941#endif /* LDA_H_ */
This header should include pretty much everything needed for the parallel runtime.
int integer
Definition crayio.c:25
double(* f)(const coord_3d &)
Definition derivatives.cc:54
double(* f1)(const coord_3d &)
Definition derivatives.cc:55
double(* f3)(const coord_3d &)
Definition derivatives.cc:57
double(* f2)(const coord_3d &)
Definition derivatives.cc:56
Fcwf copy(Fcwf psi)
Definition fcwf.cc:338
static const double v
Definition hatom_sf_dirac.cc:20
double ss1
Definition lapack.cc:62
int c_uks_vwn5__(double *ra, double *rb, double *f, double *dfdra, double *dfdrb)
Definition lda.h:246
const double THRESH_GRHO
Definition lda.h:831
void wst_munge_grho(int npoint, double *rho, double *grho)
Definition lda.h:833
double doublereal
Definition lda.h:47
int c_rks_vwn5__(const double *r__, double *f, double *dfdra)
Definition lda.h:166
static double pow(const double *a, const double *b)
Definition lda.h:74
void xc_rks_generic_lda(Tensor< double > rho_alpha, Tensor< double > f, Tensor< double > df_drho)
Definition lda.h:848
int rks_c_vwn5__(integer *ideriv, integer *npt, doublereal *rhoa1, doublereal *sigmaaa1, doublereal *zk, doublereal *vrhoa, doublereal *vsigmaaa, doublereal *v2rhoa2, doublereal *v2rhoasigmaaa, doublereal *v2sigmaaa2)
Definition lda.h:431
int logical
Definition lda.h:49
static double c_b14
Definition lda.h:81
int x_rks_s__(const double *r__, double *f, double *dfdra)
Definition lda.h:84
void dft_xc_lda_V(const Key< NDIM > &key, Tensor< double > &t)
Definition lda.h:900
void dft_xc_lda_ene(const Key< NDIM > &key, Tensor< double > &t)
Definition lda.h:911
static double c_b8
Definition lda.h:80
MADNESS_FORINT integer
Definition lda.h:48
int rks_x_lda__(integer *ideriv, integer *npt, doublereal *rhoa1, doublereal *sigmaaa1, doublereal *zk, doublereal *vrhoa, doublereal *vsigmaaa, doublereal *v2rhoa2, doublereal *v2rhoasigmaaa, doublereal *v2sigmaaa2)
Definition lda.h:712
int x_uks_s__(double *ra, double *rb, double *f, double *dfdra, double *dfdrb)
Definition lda.h:125
const double THRESH_RHO
Definition lda.h:830
static double c_b2
Definition lda.h:78
static double c_b7
Definition lda.h:79
void wst_munge_rho(int npoint, double *rho)
Definition lda.h:841
#define max(a, b)
Definition lda.h:51
Macros and tools pertaining to the configuration of MADNESS.
#define MADNESS_ASSERT(condition)
Assert a condition that should be free of side-effects since in release builds this might be a no-op.
Definition madness_exception.h:134
Main include file for MADNESS and defines Function interface.
static const double b
Definition nonlinschro.cc:119
static const double a
Definition nonlinschro.cc:118
Definition test_dc.cc:47
static double V(const coordT &r)
Definition tdse.cc:288
void e()
Definition test_sig.cc:75
const double zeta
Definition vnucso.cc:82
const double a2
Definition vnucso.cc:86
const double a1
Definition vnucso.cc:85