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>
43 #include <madness/world/MADworld.h>
44 #include <math.h>
45 #include <madness/madness_config.h>
46 
47 typedef double doublereal;
48 typedef MADNESS_FORINT integer;
49 typedef 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 
74 static inline double pow(const double* a, const double* b) {
75  return pow(*a, *b);
76 }
77 
78 static double c_b2 = .333333333333333333333333333333333;
79 static double c_b7 = .333333333333333333333333333333;
80 static double c_b8 = .5;
81 static double c_b14 = 1.333333333333333333333333333333;
82 
83 
84 inline /* 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 /* ----------------------------------------------------------------------- */
125 inline /* 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 
166 inline /* 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 /* ----------------------------------------------------------------------- */
246 inline /* 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 
431 inline /* 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 
712 inline /* 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 
830 const double THRESH_RHO = 1e-8;
831 const double THRESH_GRHO = 1e-20;
832 
833 inline 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 
841 inline 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 //***************************************************************************
848 inline 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 std::vector< Slice > ___
Entire dimension.
Definition: slice.h:128
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