51 #define max(a,b) ((a)>=(b)?(a):(b))
74 static inline double pow(
const double*
a,
const double*
b) {
78 static double c_b2 = .333333333333333333333333333333333;
79 static double c_b7 = .333333333333333333333333333333;
81 static double c_b14 = 1.333333333333333333333333333333;
84 inline int x_rks_s__(
const double *r__,
double *
f,
double *
117 ra13 =
pow(r__, &
c_b2) * .793700525984099737375852819636154;
118 *
f = *r__ * -.930525736349100025002010218071667 * ra13;
119 *dfdra = ra13 * -1.24070098179880003333601362409556;
126 double *dfdra,
double *dfdrb)
159 *
f = (*ra * ra13 + *rb * rb13) * -.930525736349100025002010218071667;
160 *dfdra = ra13 * -1.24070098179880003333601362409556;
161 *dfdrb = rb13 * -1.24070098179880003333601362409556;
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;
205 t4 = .620350490899399531;
208 t5 = 1.92366105093153617;
211 t6 = 2.56488140124204822;
214 t7 = .584822362263464735;
218 p1 = 6.1519908197590798;
220 p3 = 9.6902277115443745e-4;
221 p4 = .038783294878113009;
231 inv = 1. / (iv2 + b2 * iv + c2);
233 i2 = log((iv - d2) * (iv - d2) * inv);
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.));
240 *dfdra = pp1 - iv * .166666666666666666666666666666 * pp2;
247 f,
double *dfdra,
double *dfdrb)
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;
288 a1 = -.0337737278807791058;
305 t4 = .620350490899399531;
308 t5 = 1.92366105093153617;
311 t6 = 2.56488140124204822;
314 t7 = .584822362263464735;
318 s1 = 7.12310891781811772;
320 s3 = -6.9917323507644313e-6;
321 s4 = -.0053650918488835769;
325 p1 = 6.1519908197590798;
327 p3 = 9.6902277115443745e-4;
328 p4 = .038783294878113009;
332 f1 = 4.73092690956011364;
340 f3 = .00224786709554261133;
341 f4 = .0524913931697809227;
345 inter1 = .99999999989999999;
346 inter2 = -.99999999989999999;
349 alpha_rho13__ =
pow(ra, &
c_b7);
357 inv = 1. / (iv2 + b1 * iv + c1);
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.));
365 inv = 1. / (iv2 + b2 * iv + c2);
367 i2 = log((iv - d2) * (iv - d2) * inv);
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.));
374 inv = 1. / (iv2 + b3 * iv + c3);
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.));
383 zeta = (*ra - *rb) / srho;
385 zeta4 = zeta3 *
zeta;
387 vwn1 = t5 * .51984209978974638;
388 vwn2 = t6 * 1.25992104989487316476721060727823;
389 }
else if (
zeta < inter2) {
390 vwn1 = t5 * .51984209978974638;
391 vwn2 = t6 * -1.25992104989487316476721060727823;
402 tau = ff1 - pp1 -
ss1;
403 dtau = ff2 - pp2 - ss2;
405 v = pp1 + vwn1 * (
ss1 + tau * zeta4);
408 t1 =
v - iv * .166666666666666666666666666667 * (pp2 + vwn1 * (ss2 + dtau
410 t2 = vwn2 * (
ss1 + tau * zeta4) + vwn1 * 4. * tau * zeta3;
411 *dfdra = t1 + t2 * (1. -
zeta);
412 *dfdrb = t1 - t2 * (
zeta + 1.);
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;
499 for (i__ = 1; i__ <= i__1; ++i__) {
501 d__1 = 0., d__2 = rhoa1[i__];
502 rho =
max(d__1,d__2);
505 t2 = pow_dd(&t1, &c_b3);
506 t4 = pow_dd(&t1, &
c_b2);
507 t7 = 1 / (t2 * .6203504908994 + t4 * 2.935818660072219 +
509 t10 = log(t2 * .6203504908994 * t7);
510 t16 = atan(6.15199081975908 / (t4 * 1.575246635799487 +
513 d__1 = t4 * .7876233178997433 + .10498;
516 zk[i__] = rho * (t10 * .0310907 + t16 * .03878329487811301 +
517 t22 * 9.690227711544374e-4);
524 }
else if (*ideriv == 1) {
526 for (i__ = 1; i__ <= i__1; ++i__) {
528 d__1 = 0., d__2 = rhoa1[i__];
529 rho =
max(d__1,d__2);
532 t2 = pow_dd(&t1, &c_b3);
533 t4 = pow_dd(&t1, &
c_b2);
534 t6 = t2 * .6203504908994 + t4 * 2.935818660072219 + 12.9352;
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;
546 t23 = t22 * 9.690227711544374e-4;
547 zk[i__] = rho * (t11 + t17 + t23);
568 t43 = t26 * -.2067834969664667 * t29 - t41 *
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);
587 }
else if (*ideriv == 2) {
589 for (i__ = 1; i__ <= i__1; ++i__) {
591 d__1 = 0., d__2 = rhoa1[i__];
592 rho =
max(d__1,d__2);
595 t2 = pow_dd(&t1, &c_b3);
596 t4 = pow_dd(&t1, &
c_b2);
597 t6 = t2 * .6203504908994 + t4 * 2.935818660072219 + 12.9352;
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;
609 t23 = t22 * 9.690227711544374e-4;
610 zk[i__] = rho * (t11 + t17 + t23);
634 t43 = t26 * -.2067834969664667 * t29 - t41 *
636 t46 = t27 * -.2067834969664667 * t29 - t34 * .6203504908994 *
646 t55 = t52 * 37.8469910464 + 1.;
648 t58 = t53 * t29 * t56;
651 t66 = t60 * -.2625411059665811 * t41 - t63 * 1. * t43;
655 vrhoa[i__] = t11 + t17 + t23 + rho * (t49 *
656 .05011795824473985 + t58 * .0626408570946439 + t69 *
657 9.690227711544374e-4);
663 t81 = t77 * t7 * t80;
672 t107 = t77 * -.1378556646443111 * t80 + t26 *
673 .4135669939329333 * t88 - t103 * .4077525916766971 +
674 t105 * .978606220024073;
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);
773 for (i__ = 1; i__ <= i__1; ++i__) {
775 d__1 = 0., d__2 = rhoa1[i__];
776 rho =
max(d__1,d__2);
778 t1 = pow_dd(&rho, &
c_b2);
779 zk[i__] = t1 * -.7385587663820224 * rho;
786 }
else if (*ideriv == 1) {
788 for (i__ = 1; i__ <= i__1; ++i__) {
790 d__1 = 0., d__2 = rhoa1[i__];
791 rho =
max(d__1,d__2);
793 t1 = pow_dd(&rho, &
c_b2);
794 zk[i__] = t1 * -.7385587663820224 * rho;
795 vrhoa[i__] = t1 * -.9847450218426965;
803 }
else if (*ideriv == 2) {
805 for (i__ = 1; i__ <= i__1; ++i__) {
807 d__1 = 0., d__2 = rhoa1[i__];
808 rho =
max(d__1,d__2);
810 t1 = pow_dd(&rho, &
c_b2);
811 zk[i__] = t1 * -.7385587663820224 * rho;
812 vrhoa[i__] = t1 * -.9847450218426965;
816 v2rhoa2[i__] = -.6564966812284644 / t5;
834 for (
int i=0; i<npoint; i++) {
842 for (
int i=0; i<npoint; i++) {
850 Tensor<double> df_drho)
856 rho_alpha = rho_alpha.flat();
858 df_drho = df_drho.flat();
861 integer npt = rho_alpha.dim(0);
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);
876 int returnvalue =
::rks_x_lda__(&ideriv, &npt, rho_alpha.ptr(), gamma_alpha.ptr(),
878 tdf_drho.ptr(), tdf_dgamma.ptr(),
879 td2f_drho2.ptr(), td2f_drhodgamma.ptr(), td2f_dgamma2.ptr());
881 f.gaxpy(1.0, tf, 1.0);
882 df_drho.gaxpy(1.0, tdf_drho, 1.0);
887 returnvalue =
::rks_c_vwn5__(&ideriv, &npt, rho_alpha.ptr(), gamma_alpha.ptr(),
889 tdf_drho.ptr(), tdf_dgamma.ptr(),
890 td2f_drho2.ptr(), td2f_drhodgamma.ptr(), td2f_dgamma2.ptr());
892 f.gaxpy(1.0, tf, 1.0);
893 df_drho.gaxpy(1.0, tdf_drho, 1.0);
902 Tensor<double> enefunc =
copy(t);
903 Tensor<double>
V =
copy(t);
913 Tensor<double>
V =
copy(t);
914 Tensor<double> enefunc =
copy(t);
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