1 #include "catch.hpp"
2 #include "fixture.h"
3 
4 TEST_CASE_METHOD(libint2::unit::DefaultFixture, "electrostatic potential", "[engine][1-body]") {
5 #if defined(LIBINT2_SUPPORT_ONEBODY)
6   if (LIBINT_SHGSHELL_ORDERING != LIBINT_SHGSHELL_ORDERING_STANDARD)
7     return;
8 
9   std::vector<Shell> obs{
10       Shell{{1.0, 3.0}, {{2, true, {1.0, 0.3}}}, {{0.0, 0.0, 0.0}}},
11       Shell{{2.0, 5.0}, {{2, true, {1.0, 0.2}}}, {{1.0, 1.0, 1.0}}}};
12   {
13     const auto lmax = std::min(3, LIBINT2_MAX_AM_elecpot);
14     if (lmax >= 2) {
15       auto engine = Engine(Operator::nuclear, 2, lmax);
16       engine.set_params(make_point_charges(atoms));
17 
18       const auto scale = 2.3;
19       engine.prescale_by(scale);
20       engine.compute(obs[0], obs[0]);
21       {
22         std::vector<double> shellset_ref = {
23             -1.238239259091998e+01, 0.000000000000000e+00,
24             0.000000000000000e+00,  -5.775996163160049e-02,
25             0.000000000000000e+00,  0.000000000000000e+00,
26             -1.301230978657952e+01, -6.796143730068988e-02,
27             0.000000000000000e+00,  1.139389632827834e-01,
28             0.000000000000000e+00,  -6.796143730068988e-02,
29             -1.343732979153083e+01, 0.000000000000000e+00,
30             -1.478824785355970e-02, -5.775996163160049e-02,
31             0.000000000000000e+00,  0.000000000000000e+00,
32             -1.284475452992947e+01, 0.000000000000000e+00,
33             0.000000000000000e+00,  1.139389632827834e-01,
34             -1.478824785355970e-02, 0.000000000000000e+00,
35             -1.241040347301479e+01};
36         for (int i = 0; i != 25; ++i) {
37           REQUIRE(engine.results()[0][i]/scale == Approx(shellset_ref[i]));
38         }
39       }
40 
41       engine.prescale_by(1);
42       engine.compute(obs[0], obs[1]);
43       {
44         std::vector<double> shellset_ref = {
45             -4.769186621041819e-01, -9.303619356400431e-01,
46             -1.559058302243514e+00, -9.290824121864600e-01,
47             -5.835786921473129e-04, -1.159266418436018e+00,
48             -3.770080831197964e-01, 9.572841308198474e-01,
49             -8.291498398421207e-01, -1.663667687168316e+00,
50             -2.171951144148577e+00, 1.074249956874296e+00,
51             2.128355904665372e+00,  1.074590109905394e+00,
52             -3.485163651594458e-03, -1.160865205880651e+00,
53             -8.344173649626901e-01, 9.566621490332916e-01,
54             -3.760919234260182e-01, 1.660514988916377e+00,
55             -1.120272634615116e-03, -1.385603731947886e+00,
56             -2.105750177166632e-03, 1.380654897976564e+00,
57             2.115041199099945e+00};
58         for (int i = 0; i != 25; ++i) {
59           REQUIRE(engine.results()[0][i] == Approx(shellset_ref[i]));
60         }
61       }
62     }
63   }
64 
65   // see https://github.com/evaleev/libint/issues/199
66   {
67     const auto lmax = std::min(3, LIBINT2_MAX_AM_elecpot);
68     if (lmax >= 2) {
69       const auto deriv_order = INCLUDE_ONEBODY;
70       auto engine = Engine(Operator::nuclear, 2, lmax, deriv_order);
71       engine.set_params(make_point_charges(atoms));
72       const auto& buf = engine.results();
73       REQUIRE(libint2::num_geometrical_derivatives(atoms.size() + 2, deriv_order) == buf.size());
74     }
75   }
76 
77 #endif // LIBINT2_SUPPORT_ONEBODY
78 }
79