1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
|
//
// Created by Eduard Valeyev on 8/11/18.
//
#include "catch.hpp"
#include <libint2.hpp>
#if defined(LIBINT2_SUPPORT_ERI) && LIBINT2_MAX_AM_eri >= 1
extern "C" {
void init_c_api(unsigned int max_am);
double *
compute_eri(unsigned int am1, double alpha1, double *A,
unsigned int am2, double alpha2, double *B,
unsigned int am3, double alpha3, double *C,
unsigned int am4, double alpha4, double *D);
void finalize_c_api();
}
TEST_CASE("C API", "[c-api]") {
int am1, am2, am3, am4;
double alpha1, alpha2, alpha3, alpha4;
double A[3], B[3], C[3], D[3];
am1 = am2 = am3 = am4 = 1;
alpha1 = 1.1;
alpha2 = 2.3;
alpha3 = 3.4;
alpha4 = 4.8;
A[0] = 0.0; A[1] = 1.0; A[2] = 2.0;
B[0] = 1.0; B[1] = 2.0; B[2] = 0.0;
C[0] = 2.0; C[1] = 0.0; C[2] = 1.0;
D[0] = 0.0; D[1] = 1.0; D[2] = 2.0;
using std::max;
auto max_am = max(max(am1,am2),max(am3,am4));
init_c_api(max_am);
auto* c_result = compute_eri(am1, alpha1, A, am2, alpha2, B, am3, alpha3, C, am4, alpha4, D);
const double* cpp_result;
using libint2::Shell;
Shell sh1{{alpha1}, {{am1, false, {1.0}}}, {A[0], A[1], A[2]}};
Shell sh2{{alpha2}, {{am2, false, {1.0}}}, {B[0], B[1], B[2]}};
Shell sh3{{alpha3}, {{am3, false, {1.0}}}, {C[0], C[1], C[2]}};
Shell sh4{{alpha4}, {{am4, false, {1.0}}}, {D[0], D[1], D[2]}};
libint2::Engine engine(libint2::Operator::coulomb, 1, max_am);
engine.compute(sh1, sh2, sh3, sh4);
cpp_result = engine.results()[0];
unsigned int n1, n2, n3, n4;
int a, b, c, d, abcd;
n1 = (am1 + 1) * (am1 + 2)/2;
n2 = (am2 + 1) * (am2 + 2)/2;
n3 = (am3 + 1) * (am3 + 2)/2;
n4 = (am4 + 1) * (am4 + 2)/2;
const auto norm_factor = sh1.contr[0].coeff[0] * sh2.contr[0].coeff[0] * sh3.contr[0].coeff[0] * sh4.contr[0].coeff[0];
for(a=0, abcd=0; a<n1; a++) {
for(b=0; b<n2; b++) {
for(c=0; c<n3; c++) {
for(d=0; d<n4; d++, ++abcd) {
//printf("a = %d b = %d c = %d d = %d (ab|cd) = %20.15lf , ref (ab|cd) = %20.15lf\n", a, b, c, d, c_result[abcd]*norm_factor, cpp_result[abcd]);
REQUIRE(c_result[abcd]*norm_factor == Approx(cpp_result[abcd]));
}
}
}
}
finalize_c_api();
}
#endif
|
