diff --git a/src/QMCWaveFunctions/BsplineFactory/HybridRepCenterOrbitals.h b/src/QMCWaveFunctions/BsplineFactory/HybridRepCenterOrbitals.h index a0c75fb8a8..462496a73b 100644 --- a/src/QMCWaveFunctions/BsplineFactory/HybridRepCenterOrbitals.h +++ b/src/QMCWaveFunctions/BsplineFactory/HybridRepCenterOrbitals.h @@ -409,6 +409,31 @@ class HybridRepCenterOrbitals using RealType = typename DistanceTable::RealType; using PosType = typename DistanceTable::PosType; + enum class Region + { + INSIDE, // within the buffer shell + BUFFER, // in the buffer region + INTER // interstitial area + }; + + struct LocationSmoothingInfo + { + ///r from distance table + RealType dist_r; + ///dr from distance table + PosType dist_dr; + ///for APBC + PointType r_image; + /// region of the location + Region region; + ///smooth function value + RealType f; + ///smooth function first derivative + RealType df_dr; + ///smooth function second derivative + RealType d2f_dr2; + }; + private: ///atomic centers std::vector> AtomicCenters; @@ -416,18 +441,6 @@ class HybridRepCenterOrbitals int myTableID; ///mapping supercell to primitive cell std::vector Super2Prim; - ///r from distance table - RealType dist_r; - ///dr from distance table - PosType dist_dr; - ///for APBC - PointType r_image; - ///smooth function value - RealType f; - ///smooth function first derivative - RealType df_dr; - ///smooth function second derivative - RealType d2f_dr2; ///smoothing schemes enum class smoothing_schemes { @@ -438,6 +451,28 @@ class HybridRepCenterOrbitals /// smoothing function smoothing_functions smooth_func_id; + /// select a region (within the buffer shell, in the buffer, interstitial region) and compute the smoothing function if in the buffer. + inline void selectRegionAndComputeSmoothing(const ST& cutoff_buffer, + const ST& cutoff, + LocationSmoothingInfo& info) const + { + const RealType r = info.dist_r; + if (r < cutoff_buffer) + info.region = Region::INSIDE; + else if (r < cutoff) + { + constexpr RealType cone(1); + const RealType scale = cone / (cutoff - cutoff_buffer); + const RealType x = (r - cutoff_buffer) * scale; + info.f = smoothing(smooth_func_id, x, info.df_dr, info.d2f_dr2); + info.df_dr *= scale; + info.d2f_dr2 *= scale * scale; + info.region = Region::BUFFER; + } + else + info.region = Region::INTER; + } + public: HybridRepCenterOrbitals() {} @@ -553,7 +588,10 @@ class HybridRepCenterOrbitals } template - inline int get_bc_sign(const PointType& r, const Cell& PrimLattice, TinyVector& HalfG) + inline int get_bc_sign(const PointType& r, + const PointType& r_image, + const Cell& PrimLattice, + TinyVector& HalfG) const { int bc_sign = 0; PointType shift_unit = PrimLattice.toUnit(r - r_image); @@ -567,21 +605,18 @@ class HybridRepCenterOrbitals //evaluate only V template - inline RealType evaluate_v(const ParticleSet& P, const int iat, VV& myV) + inline void evaluate_v(const ParticleSet& P, const int iat, VV& myV, LocationSmoothingInfo& info) { const auto& ei_dist = P.getDistTableAB(myTableID); - const int center_idx = ei_dist.get_first_neighbor(iat, dist_r, dist_dr, P.getActivePtcl() == iat); - if (center_idx < 0) - abort(); - auto& myCenter = AtomicCenters[Super2Prim[center_idx]]; - if (dist_r < myCenter.getCutoff()) + const int center_idx = ei_dist.get_first_neighbor(iat, info.dist_r, info.dist_dr, P.getActivePtcl() == iat); + auto& myCenter = AtomicCenters[Super2Prim[center_idx]]; + selectRegionAndComputeSmoothing(myCenter.getCutoffBuffer(), myCenter.getCutoff(), info); + if (info.region != Region::INTER) { - PointType dr(-dist_dr[0], -dist_dr[1], -dist_dr[2]); - r_image = myCenter.getCenterPos() + dr; - myCenter.evaluate_v(dist_r, dr, myV); - return smooth_function(myCenter.getCutoffBuffer(), myCenter.getCutoff(), dist_r); + PointType dr(-info.dist_dr[0], -info.dist_dr[1], -info.dist_dr[2]); + info.r_image = myCenter.getCenterPos() + dr; + myCenter.evaluate_v(info.dist_r, dr, myV); } - return RealType(-1); } /* check if the batched algorithm is safe to operate @@ -593,7 +628,7 @@ class HybridRepCenterOrbitals * The batched algorthm forces the evaluation on the reference center and introduce some error. * In this case, the non-batched algorithm should be used. */ - bool is_batched_safe(const VirtualParticleSet& VP) + bool is_batched_safe(const VirtualParticleSet& VP) const { const int center_idx = VP.refSourcePtcl; auto& myCenter = AtomicCenters[Super2Prim[center_idx]]; @@ -603,88 +638,75 @@ class HybridRepCenterOrbitals // C2C, C2R cases template - inline RealType evaluateValuesC2X(const VirtualParticleSet& VP, VM& multi_myV) + inline void evaluateValuesC2X(const VirtualParticleSet& VP, VM& multi_myV, LocationSmoothingInfo& info) { const int center_idx = VP.refSourcePtcl; - dist_r = VP.getRefPS().getDistTableAB(myTableID).getDistRow(VP.refPtcl)[center_idx]; + info.dist_r = VP.getRefPS().getDistTableAB(myTableID).getDistRow(VP.refPtcl)[center_idx]; auto& myCenter = AtomicCenters[Super2Prim[center_idx]]; - if (dist_r < myCenter.getCutoff()) - { - myCenter.evaluateValues(VP.getDistTableAB(myTableID).getDisplacements(), center_idx, dist_r, multi_myV); - return smooth_function(myCenter.getCutoffBuffer(), myCenter.getCutoff(), dist_r); - } - return RealType(-1); + selectRegionAndComputeSmoothing(myCenter.getCutoffBuffer(), myCenter.getCutoff(), info); + if (info.region != Region::INTER) + myCenter.evaluateValues(VP.getDistTableAB(myTableID).getDisplacements(), center_idx, info.dist_r, multi_myV); } // R2R case template - inline RealType evaluateValuesR2R(const VirtualParticleSet& VP, - const Cell& PrimLattice, - TinyVector& HalfG, - VM& multi_myV, - SV& bc_signs) + inline void evaluateValuesR2R(const VirtualParticleSet& VP, + const Cell& PrimLattice, + TinyVector& HalfG, + VM& multi_myV, + SV& bc_signs, + LocationSmoothingInfo& info) { const int center_idx = VP.refSourcePtcl; - dist_r = VP.getRefPS().getDistTableAB(myTableID).getDistRow(VP.refPtcl)[center_idx]; + info.dist_r = VP.getRefPS().getDistTableAB(myTableID).getDistRow(VP.refPtcl)[center_idx]; auto& myCenter = AtomicCenters[Super2Prim[center_idx]]; - if (dist_r < myCenter.getCutoff()) + selectRegionAndComputeSmoothing(myCenter.getCutoffBuffer(), myCenter.getCutoff(), info); + if (info.region != Region::INTER) { const auto& displ = VP.getDistTableAB(myTableID).getDisplacements(); for (int ivp = 0; ivp < VP.getTotalNum(); ivp++) - { - r_image = myCenter.getCenterPos() - displ[ivp][center_idx]; - bc_signs[ivp] = get_bc_sign(VP.R[ivp], PrimLattice, HalfG); - ; - } - myCenter.evaluateValues(displ, center_idx, dist_r, multi_myV); - return smooth_function(myCenter.getCutoffBuffer(), myCenter.getCutoff(), dist_r); + bc_signs[ivp] = get_bc_sign(VP.R[ivp], myCenter.getCenterPos() - displ[ivp][center_idx], PrimLattice, HalfG); + myCenter.evaluateValues(displ, center_idx, info.dist_r, multi_myV); } - return RealType(-1); } //evaluate only VGL template - inline RealType evaluate_vgl(const ParticleSet& P, const int iat, VV& myV, GV& myG, VV& myL) + inline void evaluate_vgl(const ParticleSet& P, const int iat, VV& myV, GV& myG, VV& myL, LocationSmoothingInfo& info) { const auto& ei_dist = P.getDistTableAB(myTableID); - const int center_idx = ei_dist.get_first_neighbor(iat, dist_r, dist_dr, P.getActivePtcl() == iat); - if (center_idx < 0) - abort(); - auto& myCenter = AtomicCenters[Super2Prim[center_idx]]; - if (dist_r < myCenter.getCutoff()) + const int center_idx = ei_dist.get_first_neighbor(iat, info.dist_r, info.dist_dr, P.getActivePtcl() == iat); + auto& myCenter = AtomicCenters[Super2Prim[center_idx]]; + selectRegionAndComputeSmoothing(myCenter.getCutoffBuffer(), myCenter.getCutoff(), info); + if (info.region != Region::INTER) { - PointType dr(-dist_dr[0], -dist_dr[1], -dist_dr[2]); - r_image = myCenter.getCenterPos() + dr; - myCenter.evaluate_vgl(dist_r, dr, myV, myG, myL); - return smooth_function(myCenter.getCutoffBuffer(), myCenter.getCutoff(), dist_r); + const PointType dr(-info.dist_dr[0], -info.dist_dr[1], -info.dist_dr[2]); + info.r_image = myCenter.getCenterPos() + dr; + myCenter.evaluate_vgl(info.dist_r, dr, myV, myG, myL); } - return RealType(-1); } //evaluate only VGH template - inline RealType evaluate_vgh(const ParticleSet& P, const int iat, VV& myV, GV& myG, HT& myH) + inline void evaluate_vgh(const ParticleSet& P, const int iat, VV& myV, GV& myG, HT& myH, LocationSmoothingInfo& info) { const auto& ei_dist = P.getDistTableAB(myTableID); - const int center_idx = ei_dist.get_first_neighbor(iat, dist_r, dist_dr, P.getActivePtcl() == iat); - if (center_idx < 0) - abort(); - auto& myCenter = AtomicCenters[Super2Prim[center_idx]]; - if (dist_r < myCenter.getCutoff()) + const int center_idx = ei_dist.get_first_neighbor(iat, info.dist_r, info.dist_dr, P.getActivePtcl() == iat); + auto& myCenter = AtomicCenters[Super2Prim[center_idx]]; + selectRegionAndComputeSmoothing(myCenter.getCutoffBuffer(), myCenter.getCutoff(), info); + if (info.region != Region::INTER) { - PointType dr(-dist_dr[0], -dist_dr[1], -dist_dr[2]); - r_image = myCenter.getCenterPos() + dr; - myCenter.evaluate_vgh(dist_r, dr, myV, myG, myH); - return smooth_function(myCenter.getCutoffBuffer(), myCenter.getCutoff(), dist_r); + const PointType dr(-info.dist_dr[0], -info.dist_dr[1], -info.dist_dr[2]); + info.r_image = myCenter.getCenterPos() + dr; + myCenter.evaluate_vgh(info.dist_r, dr, myV, myG, myH); } - return RealType(-1); } // interpolate buffer region, value only template - inline void interpolate_buffer_v(VV& psi, const VV& psi_AO) const + inline void interpolate_buffer_v(VV& psi, const VV& psi_AO, const RealType f) const { - const RealType cone(1); + constexpr RealType cone(1); for (size_t i = 0; i < psi.size(); i++) psi[i] = psi_AO[i] * f + psi[i] * (cone - f); } @@ -696,10 +718,15 @@ class HybridRepCenterOrbitals VV& d2psi, const VV& psi_AO, const GV& dpsi_AO, - const VV& d2psi_AO) const + const VV& d2psi_AO, + const LocationSmoothingInfo& info) const { - const RealType cone(1), ctwo(2); - const RealType rinv(1.0 / dist_r); + constexpr RealType cone(1), ctwo(2); + const RealType rinv(1.0 / info.dist_r); + auto& dist_dr = info.dist_dr; + auto& f = info.f; + auto& df_dr = info.df_dr; + auto& d2f_dr2 = info.d2f_dr2; if (smooth_scheme == smoothing_schemes::CONSISTENT) for (size_t i = 0; i < psi.size(); i++) { // psi, dpsi, d2psi are all consistent @@ -726,19 +753,6 @@ class HybridRepCenterOrbitals throw std::runtime_error("Unknown smooth scheme!"); } - inline RealType smooth_function(const ST& cutoff_buffer, const ST& cutoff, const RealType r) - { - const RealType cone(1); - if (r < cutoff_buffer) - return cone; - const RealType scale = cone / (cutoff - cutoff_buffer); - const RealType x = (r - cutoff_buffer) * scale; - f = smoothing(smooth_func_id, x, df_dr, d2f_dr2); - df_dr *= scale; - d2f_dr2 *= scale * scale; - return f; - } - template friend class qmcplusplus::HybridRepSetReader; }; diff --git a/src/QMCWaveFunctions/BsplineFactory/HybridRepCplx.h b/src/QMCWaveFunctions/BsplineFactory/HybridRepCplx.h index 9dec1398a0..dcc6aeb87d 100644 --- a/src/QMCWaveFunctions/BsplineFactory/HybridRepCplx.h +++ b/src/QMCWaveFunctions/BsplineFactory/HybridRepCplx.h @@ -47,9 +47,12 @@ class HybridRepCplx : public SPLINEBASE, private HybridRepCenterOrbitals> multi_myV; + typename HYBRIDBASE::LocationSmoothingInfo info; using SPLINEBASE::myG; using SPLINEBASE::myH; @@ -95,24 +98,17 @@ class HybridRepCplx : public SPLINEBASE, private HybridRepCenterOrbitals> myV_one(multi_myV[iat], myV.size()); - SPLINEBASE::assign_v(r, myV_one, psi, 0, myV.size() / 2); + SPLINEBASE::assign_v(VP.R[iat], myV_one, psi, 0, myV.size() / 2); } else { - const PointType& r = VP.R[iat]; Vector> myV_one(multi_myV[iat], myV.size()); - SPLINEBASE::assign_v(r, myV_one, psi_AO, 0, myV.size() / 2); + SPLINEBASE::assign_v(VP.R[iat], myV_one, psi_AO, 0, myV.size() / 2); SPLINEBASE::evaluateValue(VP, iat, psi); - HYBRIDBASE::interpolate_buffer_v(psi, psi_AO); + HYBRIDBASE::interpolate_buffer_v(psi, psi_AO, info.f); } ratios[iat] = simd::dot(psi.data(), psiinv.data(), psi.size()); } @@ -172,26 +165,19 @@ class HybridRepCplx : public SPLINEBASE, private HybridRepCenterOrbitals> multi_myV; + typename HYBRIDBASE::LocationSmoothingInfo info; using SPLINEBASE::HalfG; using SPLINEBASE::myG; @@ -97,26 +100,21 @@ class HybridRepReal : public SPLINEBASE, private HybridRepCenterOrbitals bc_signs(VP.getTotalNum()); - const RealType smooth_factor = HYBRIDBASE::evaluateValuesR2R(VP, PrimLattice, HalfG, multi_myV, bc_signs); - const RealType cone(1); + HYBRIDBASE::evaluateValuesR2R(VP, PrimLattice, HalfG, multi_myV, bc_signs, info); for (int iat = 0; iat < VP.getTotalNum(); ++iat) { - if (smooth_factor < 0) + if (info.region == Region::INTER) SPLINEBASE::evaluateValue(VP, iat, psi); - else if (smooth_factor == cone) + else if (info.region == Region::INSIDE) { Vector> myV_one(multi_myV[iat], myV.size()); SPLINEBASE::assign_v(bc_signs[iat], myV_one, psi, 0, myV.size()); @@ -148,19 +145,17 @@ class HybridRepReal : public SPLINEBASE, private HybridRepCenterOrbitals> myV_one(multi_myV[iat], myV.size()); SPLINEBASE::assign_v(bc_signs[iat], myV_one, psi_AO, 0, myV.size()); SPLINEBASE::evaluateValue(VP, iat, psi); - HYBRIDBASE::interpolate_buffer_v(psi, psi_AO); + HYBRIDBASE::interpolate_buffer_v(psi, psi_AO, info.f); } ratios[iat] = simd::dot(psi.data(), psiinv.data(), psi.size()); } } else - { for (int iat = 0; iat < VP.getTotalNum(); ++iat) { evaluateValue(VP, iat, psi); ratios[iat] = simd::dot(psi.data(), psiinv.data(), psi.size()); } - } } void mw_evaluateDetRatios(const RefVectorWithLeader& spo_list, @@ -174,28 +169,21 @@ class HybridRepReal : public SPLINEBASE, private HybridRepCenterOrbitalsevaluateVGL(elec, 0, psiref_0, dpsiref_0, d2psiref_0); @@ -545,15 +543,10 @@ void test_Ne(bool transform) sposet->evaluateValue(elec, 0, values); - std::cout << "values = " << values << std::endl; - // Generated from gen_mo.py for position [1e-05, 0.0, 0.0] CHECK(values[0] == Approx(-16.11819042)); sposet->evaluateVGL(elec, 0, values, dpsi, d2psi); - std::cout << "values = " << values << std::endl; - std::cout << "dpsi = " << dpsi << std::endl; - std::cout << "d2psi = " << d2psi << std::endl; // Generated from gen_mo.py for position [1e-05, 0.0, 0.0] CHECK(values[0] == Approx(-16.11819042)); CHECK(dpsi[0][0] == Approx(0.1747261458)); diff --git a/src/QMCWaveFunctions/tests/test_hybridrep.cpp b/src/QMCWaveFunctions/tests/test_hybridrep.cpp index 2ab073fdbf..b56b37d096 100644 --- a/src/QMCWaveFunctions/tests/test_hybridrep.cpp +++ b/src/QMCWaveFunctions/tests/test_hybridrep.cpp @@ -62,9 +62,10 @@ TEST_CASE("Hybridrep SPO from HDF diamond_1x1x1", "[wavefunction]") elec_.setName("elec"); ptcl.addParticleSet(std::move(elec_uptr)); - elec_.create({2}); + elec_.create({3}); elec_.R[0] = {0.4, 0.0, 0.0}; elec_.R[1] = {0.0, 1.0, 0.0}; + elec_.R[2] = {0.0, 0.5, 0.5}; SpeciesSet& tspecies = elec_.getSpeciesSet(); int upIdx = tspecies.addSpecies("u"); int chargeIdx = tspecies.addAttribute("charge"); @@ -97,7 +98,6 @@ TEST_CASE("Hybridrep SPO from HDF diamond_1x1x1", "[wavefunction]") SPOSet::ValueMatrix d2psiM(elec_.R.size(), spo->getOrbitalSetSize()); spo->evaluate_notranspose(elec_, 0, elec_.R.size(), psiM, dpsiM, d2psiM); -#if defined(QMC_COMPLEX) // due to the different ordering of bands skip the tests on CUDA+Real builds // checking evaluations, reference values are not independently generated. // electron 0 @@ -113,11 +113,7 @@ TEST_CASE("Hybridrep SPO from HDF diamond_1x1x1", "[wavefunction]") CHECK(std::real(dpsiM[0][1][2]) == Approx(-0.6386129856)); // lapl CHECK(std::real(d2psiM[0][0]) == Approx(-4.1090884209)); -#if defined(MIXED_PRECISION) CHECK(std::real(d2psiM[0][1]) == Approx(22.3851032257).epsilon(3e-5)); -#else - CHECK(std::real(d2psiM[0][1]) == Approx(22.3851032257)); -#endif // electron 1 // value @@ -133,7 +129,21 @@ TEST_CASE("Hybridrep SPO from HDF diamond_1x1x1", "[wavefunction]") // lapl CHECK(std::real(d2psiM[1][0]) == Approx(1.3313053846)); CHECK(std::real(d2psiM[1][1]) == Approx(-4.712583065)); -#endif + + // electron 2 + // value + CHECK(std::real(psiM[2][0]) == Approx(-0.8694150782)); + CHECK(std::real(psiM[2][1]) == Approx(0.780789871)); + // grad + CHECK(std::real(dpsiM[2][0][0]) == Approx(-0.0602699547)); + CHECK(std::real(dpsiM[2][0][1]) == Approx(-0.2770632381)); + CHECK(std::real(dpsiM[2][0][2]) == Approx(-0.2770632488)); + CHECK(std::real(dpsiM[2][1][0]) == Approx(-1.5982062406)); + CHECK(std::real(dpsiM[2][1][1]) == Approx(1.0020672904)); + CHECK(std::real(dpsiM[2][1][2]) == Approx(-1.9794520201)); + // lapl + CHECK(std::real(d2psiM[2][0]) == Approx(1.1232769428).epsilon(3e-5)); + CHECK(std::real(d2psiM[2][1]) == Approx(-4.9779265738).epsilon(3e-5)); } TEST_CASE("Hybridrep SPO from HDF diamond_2x1x1", "[wavefunction]") @@ -170,9 +180,10 @@ TEST_CASE("Hybridrep SPO from HDF diamond_2x1x1", "[wavefunction]") elec_.setName("elec"); ptcl.addParticleSet(std::move(elec_uptr)); - elec_.create({2}); + elec_.create({3}); elec_.R[0] = {0.4, 0.0, 0.0}; elec_.R[1] = {0.0, 1.0, 0.0}; + elec_.R[2] = {0.0, 0.5, 0.5}; SpeciesSet& tspecies = elec_.getSpeciesSet(); int upIdx = tspecies.addSpecies("u"); int chargeIdx = tspecies.addAttribute("charge"); @@ -202,9 +213,9 @@ TEST_CASE("Hybridrep SPO from HDF diamond_2x1x1", "[wavefunction]") ParticleSet::RealType r; ParticleSet::PosType dr; elec_.getDistTable(0).get_first_neighbor(0, r, dr, false); - std::cout << std::setprecision(14) << "check r^2 against dr^2. " + app_log() << std::setprecision(14) << "check r^2 against dr^2. " << "r = " << r << " dr = " << dr << std::endl; - std::cout << "abs(r^2 - dr^2) = " << std::abs(r * r - dot(dr, dr)) + app_log() << "abs(r^2 - dr^2) = " << std::abs(r * r - dot(dr, dr)) << " epsilon = " << std::numeric_limits::epsilon() << std::endl; #if defined(MIXED_PRECISION) REQUIRE(std::abs(r * r - dot(dr, dr)) < std::numeric_limits::epsilon() * 1e8); @@ -220,9 +231,6 @@ TEST_CASE("Hybridrep SPO from HDF diamond_2x1x1", "[wavefunction]") #if defined(QMC_COMPLEX) // real part - // due to the different ordering of bands skip the tests on CUDA+Real builds - // checking evaluations, reference values are not independently generated. - // electron 0 // value CHECK(std::real(psiM[0][0]) == Approx(0.6776432991)); @@ -252,9 +260,22 @@ TEST_CASE("Hybridrep SPO from HDF diamond_2x1x1", "[wavefunction]") // lapl CHECK(std::real(d2psiM[1][0]) == Approx(-1.3757134676)); CHECK(std::real(d2psiM[1][1]) == Approx(-2.4803137779)); -#endif -#if defined(QMC_COMPLEX) + // electron 2 + // value + CHECK(std::real(psiM[2][0]) == Approx(0.8841851282)); + CHECK(std::real(psiM[2][1]) == Approx(1.0331713017)); + // grad + CHECK(std::real(dpsiM[2][0][0]) == Approx(0.0688574613)); + CHECK(std::real(dpsiM[2][0][1]) == Approx(0.2735091889)); + CHECK(std::real(dpsiM[2][0][2]) == Approx(0.2666900514)); + CHECK(std::real(dpsiM[2][1][0]) == Approx(0.5398793935)); + CHECK(std::real(dpsiM[2][1][1]) == Approx(0.7525391523)); + CHECK(std::real(dpsiM[2][1][2]) == Approx(-0.1224437827)); + // lapl + CHECK(std::real(d2psiM[2][0]) == Approx(-1.1799273657).epsilon(1e-4)); + CHECK(std::real(d2psiM[2][1]) == Approx(-2.0339757673).epsilon(1e-4)); + // imaginary part // electron 0 // value @@ -270,6 +291,7 @@ TEST_CASE("Hybridrep SPO from HDF diamond_2x1x1", "[wavefunction]") // lapl CHECK(std::imag(d2psiM[0][0]) == Approx(4.0779790878).epsilon(1e-4)); CHECK(std::imag(d2psiM[0][1]) == Approx(-0.7897151113).epsilon(1e-4)); + // electron 1 // value CHECK(std::imag(psiM[1][0]) == Approx(0.9008999467)); @@ -284,6 +306,21 @@ TEST_CASE("Hybridrep SPO from HDF diamond_2x1x1", "[wavefunction]") // lapl CHECK(std::imag(d2psiM[1][0]) == Approx(-1.3757134676)); CHECK(std::imag(d2psiM[1][1]) == Approx(-2.4919104576)); + + // electron 2 + // value + CHECK(std::imag(psiM[2][0]) == Approx(0.8841851282)); + CHECK(std::imag(psiM[2][1]) == Approx(1.0291547321)); + // grad + CHECK(std::imag(dpsiM[2][0][0]) == Approx(0.0688574613)); + CHECK(std::imag(dpsiM[2][0][1]) == Approx(0.2735091889)); + CHECK(std::imag(dpsiM[2][0][2]) == Approx(0.2666900514)); + CHECK(std::imag(dpsiM[2][1][0]) == Approx(0.5376840012)); + CHECK(std::imag(dpsiM[2][1][1]) == Approx(0.7550587239)); + CHECK(std::imag(dpsiM[2][1][2]) == Approx(-0.1228616516)); + // lapl + CHECK(std::imag(d2psiM[2][0]) == Approx(-1.2044699918).epsilon(1e-4)); + CHECK(std::imag(d2psiM[2][1]) == Approx(-1.885562226).epsilon(1e-4)); #endif // test batched interfaces