diff --git a/include/gauxc/c/enums.h b/include/gauxc/c/enums.h index 095b37737..bb5d33d12 100644 --- a/include/gauxc/c/enums.h +++ b/include/gauxc/c/enums.h @@ -38,7 +38,17 @@ enum GauXC_AtomicGridSizeDefault { GauXC_AtomicGridSizeDefault_UltraFineGrid, ///< Ultrafine grid (appropriate accuracy) GauXC_AtomicGridSizeDefault_SuperFineGrid, ///< Superfine grid (most accurate) GauXC_AtomicGridSizeDefault_GM3, ///< Treutler-Ahlrichs GM3 - GauXC_AtomicGridSizeDefault_GM5 ///< Treutler-Ahlrichs GM5 + GauXC_AtomicGridSizeDefault_GM5, ///< Treutler-Ahlrichs GM5 + GauXC_AtomicGridSizeDefault_PySCF0, ///< PySCF default level 0 + GauXC_AtomicGridSizeDefault_PySCF1, ///< PySCF default level 1 + GauXC_AtomicGridSizeDefault_PySCF2, ///< PySCF default level 2 (angular points ~ fine grid) + GauXC_AtomicGridSizeDefault_PySCF3, ///< PySCF default level 3 + GauXC_AtomicGridSizeDefault_PySCF4, ///< PySCF default level 4 (radial points ~ fine grid, angular points ~ ultrafine grid) + GauXC_AtomicGridSizeDefault_PySCF5, ///< PySCF default level 5 + GauXC_AtomicGridSizeDefault_PySCF6, ///< PySCF default level 6 (radial points ~ ultrafine grid, angular points ~ superfine grid) + GauXC_AtomicGridSizeDefault_PySCF7, ///< PySCF default level 7 + GauXC_AtomicGridSizeDefault_PySCF8, ///< PySCF default level 8 + GauXC_AtomicGridSizeDefault_PySCF9 ///< PySCF default level 9 (radial points ~ superfine grid) }; /** @@ -69,9 +79,13 @@ enum GauXC_SupportedAlg { /// High-level specification of pruning schemes for atomic quadratures enum GauXC_PruningScheme { - GauXC_PruningScheme_Unpruned, ///< Unpruned atomic quadrature - GauXC_PruningScheme_Robust, ///< The "Robust" scheme of Psi4 - GauXC_PruningScheme_Treutler ///< The Treutler-Ahlrichs scheme + GauXC_PruningScheme_Unpruned, ///< Unpruned atomic quadrature + GauXC_PruningScheme_Robust, ///< The "Robust" scheme of Psi4 + GauXC_PruningScheme_Treutler, ///< The Treutler-Ahlrichs scheme + GauXC_PruningScheme_PySCF_Treutler, ///< The Treutler scheme of PySCF + GauXC_PruningScheme_PySCF_SG1, ///< The SG1 scheme of PySCF + GauXC_PruningScheme_PySCF_NWChem, ///< The NWChem scheme of PySCF + GauXC_PruningScheme_PySCF_SGX ///< The SGX scheme of PySCF }; diff --git a/include/gauxc/enums.hpp b/include/gauxc/enums.hpp index 9ea604564..eab1983d5 100644 --- a/include/gauxc/enums.hpp +++ b/include/gauxc/enums.hpp @@ -35,7 +35,17 @@ enum class AtomicGridSizeDefault { UltraFineGrid, ///< Ultrafine grid (appropriate accuracy) SuperFineGrid, ///< Superfine grid (most accurate) GM3, ///< Treutler-Ahlrichs GM3 - GM5 ///< Treutler-Ahlrichs GM5 + GM5, ///< Treutler-Ahlrichs GM5 + PySCF0, ///< PySCF default level 0 + PySCF1, ///< PySCF default level 1 + PySCF2, ///< PySCF default level 2 (angular points ~ fine grid) + PySCF3, ///< PySCF default level 3 + PySCF4, ///< PySCF default level 4 (radial points ~ fine grid, angular points ~ ultrafine grid) + PySCF5, ///< PySCF default level 5 + PySCF6, ///< PySCF default level 6 (radial points ~ ultrafine grid, angular points ~ superfine grid) + PySCF7, ///< PySCF default level 7 + PySCF8, ///< PySCF default level 8 + PySCF9 ///< PySCF default level 9 (radial points ~ superfine grid) }; /** @@ -66,9 +76,13 @@ enum class SupportedAlg { /// High-level specification of pruning schemes for atomic quadratures enum class PruningScheme { - Unpruned, ///< Unpruned atomic quadrature - Robust, ///< The "Robust" scheme of Psi4 - Treutler ///< The Treutler-Ahlrichs scheme + Unpruned, ///< Unpruned atomic quadrature + Robust, ///< The "Robust" scheme of Psi4 + Treutler, ///< The Treutler-Ahlrichs scheme + PySCF_Treutler, ///< The Treutler scheme of PySCF + PySCF_SG1, ///< The SG1 scheme of PySCF + PySCF_NWChem, ///< The NWChem scheme of PySCF + PySCF_SGX ///< The SGX scheme of PySCF }; diff --git a/include/gauxc/grid_factory.hpp b/include/gauxc/grid_factory.hpp index b2ed0a8dc..b377334ff 100644 --- a/include/gauxc/grid_factory.hpp +++ b/include/gauxc/grid_factory.hpp @@ -10,6 +10,7 @@ * See LICENSE.txt for details */ #pragma once +#include #include #include #include @@ -62,12 +63,23 @@ PrunedAtomicGridSpecification treutler_pruning_scheme( UnprunedAtomicGridSpecification ); +/// Generate a Pruning specification according to PySCF's Treutler pruning scheme +PrunedAtomicGridSpecification pyscf_treutler_pruning_scheme( + UnprunedAtomicGridSpecification +); + /// Generate a pruning specification from a specificed pruning scheme and /// an unpruned grid specification PrunedAtomicGridSpecification create_pruned_spec( PruningScheme, UnprunedAtomicGridSpecification ); +/// Generate an atom-aware pruning specification from a specified pruning scheme +/// and an unpruned grid specification +PrunedAtomicGridSpecification create_pruned_spec( + PruningScheme, AtomicNumber, UnprunedAtomicGridSpecification +); + using atomic_grid_variant = std::variant; diff --git a/include/gauxc/molgrid/defaults.hpp b/include/gauxc/molgrid/defaults.hpp index 0565647d0..dd9214d77 100644 --- a/include/gauxc/molgrid/defaults.hpp +++ b/include/gauxc/molgrid/defaults.hpp @@ -18,6 +18,8 @@ namespace GauXC { double slater_radius_64(AtomicNumber); double slater_radius_30(AtomicNumber); double clementi_radius_67(AtomicNumber); + double pyscf_slater_radius_64(AtomicNumber); + double pyscf_gill_radius_93(AtomicNumber); double uff_radius_103(AtomicNumber); double default_atomic_radius(AtomicNumber); @@ -41,9 +43,10 @@ namespace GauXC { template inline static atomic_grid_variant - create_default_pruned_grid_spec( PruningScheme scheme, Args&&... args ) { - return create_pruned_spec( scheme, - create_default_unpruned_grid_spec(std::forward(args)...) + create_default_pruned_grid_spec( PruningScheme scheme, AtomicNumber Z, + Args&&... args ) { + return create_pruned_spec( scheme, Z, + create_default_unpruned_grid_spec(Z, std::forward(args)...) ); } diff --git a/src/atomic_radii.cxx b/src/atomic_radii.cxx index 522753092..eff0f224e 100644 --- a/src/atomic_radii.cxx +++ b/src/atomic_radii.cxx @@ -349,4 +349,52 @@ double uff_radius_103(AtomicNumber _Z) { } return radius_uff_list[Z-1] * RADIUS_UFF_SCALING / DDX_BOHR_TO_ANGSTROM; } + +/// PySCF's Bragg-Slater radii from JCP 41, 3199 (1964); DOI:10.1063/1.1725697. +/// Only elements different from Slater-64 are overridden here. +double pyscf_slater_radius_64(AtomicNumber _Z) { + auto Z = _Z.get(); + if( Z >= 96 && Z <= 130 ) return pm_to_bohr(175.); + + switch(Z) { + case 1: /* H */ return pm_to_bohr(35. ); + case 2: /* He */ return pm_to_bohr(140.); + case 10: /* Ne */ return pm_to_bohr(150.); + case 18: /* Ar */ return pm_to_bohr(180.); + case 36: /* Kr */ return pm_to_bohr(190.); + case 54: /* Xe */ return pm_to_bohr(210.); + case 85: /* At */ return pm_to_bohr(145.); + case 86: /* Rn */ return pm_to_bohr(210.); + case 87: /* Fr */ return pm_to_bohr(180.); + default: return slater_radius_64(_Z); + } +} + +/// PySCF's SG1 radii from +/// P.M.W. Gill, B.G. Johnson, J.A. Pople, Chem. Phys. Letters 209 (1993) 506-512 +double pyscf_gill_radius_93(AtomicNumber _Z) { + auto Z = _Z.get(); + switch(Z) { + case 1: /* H */ return 1.0000; + case 2: /* He */ return 0.5882; + case 3: /* Li */ return 3.0769; + case 4: /* Be */ return 2.0513; + case 5: /* B */ return 1.5385; + case 6: /* C */ return 1.2308; + case 7: /* N */ return 1.0256; + case 8: /* O */ return 0.8791; + case 9: /* F */ return 0.7692; + case 10: /* Ne */ return 0.6838; + case 11: /* Na */ return 4.0909; + case 12: /* Mg */ return 3.1579; + case 13: /* Al */ return 2.5714; + case 14: /* Si */ return 2.1687; + case 15: /* P */ return 1.8750; + case 16: /* S */ return 1.6514; + case 17: /* Cl */ return 1.4754; + case 18: /* Ar */ return 1.3333; + default: return -1.; + } +} + } diff --git a/src/grid_factory.cxx b/src/grid_factory.cxx index 1836653e9..5aa94d82e 100644 --- a/src/grid_factory.cxx +++ b/src/grid_factory.cxx @@ -19,8 +19,123 @@ #include #include +#include +#include + namespace GauXC { +double pyscf_slater_radius_64(AtomicNumber); +double pyscf_gill_radius_93(AtomicNumber); + +namespace detail { + +using ll_type = IntegratorXX::LebedevLaikov; +using ll_traits = IntegratorXX::quadrature_traits; + +template +std::vector radial_points_for_spec(UnprunedAtomicGridSpecification unp) { + RadialType rq(unp.radial_size.get(), unp.radial_scale.get()); + return rq.points(); +} + +std::vector radial_points_for_spec(UnprunedAtomicGridSpecification unp) { + switch( unp.radial_quad ) { + case RadialQuad::Becke: + return radial_points_for_spec>(unp); + case RadialQuad::MuraKnowles: + return radial_points_for_spec>(unp); + case RadialQuad::MurrayHandyLaming: + return radial_points_for_spec>(unp); + case RadialQuad::TreutlerAhlrichs: + return radial_points_for_spec>(unp); + default: + GAUXC_GENERIC_EXCEPTION("Unsupported Radial Quadrature"); + abort(); + } +} + +PrunedAtomicGridSpecification make_pruned_spec( + UnprunedAtomicGridSpecification unp, const std::vector& angular_sizes) { + + if( angular_sizes.size() != static_cast(unp.radial_size.get()) ) + GAUXC_GENERIC_EXCEPTION("Invalid PySCF pruning specification"); + + std::vector pruning_regions; + size_t idx_st = 0; + for( size_t i = 1; i <= angular_sizes.size(); ++i ) { + if( i == angular_sizes.size() || angular_sizes[i] != angular_sizes[idx_st] ) { + pruning_regions.push_back({idx_st, i, AngularSize(angular_sizes[idx_st])}); + idx_st = i; + } + } + + return PrunedAtomicGridSpecification{ + unp.radial_quad, unp.radial_size, unp.radial_scale, pruning_regions + }; +} + +const std::array& pyscf_alphas(AtomicNumber Z) { + static const std::array h_he = {{0.25, 0.5, 1.0, 4.5}}; + static const std::array li_ne = {{0.1667, 0.5, 0.9, 3.5}}; + static const std::array rest = {{0.1, 0.4, 0.8, 2.5}}; + + if( Z.get() <= 2 ) return h_he; + if( Z.get() <= 10 ) return li_ne; + return rest; +} + +std::vector pyscf_alpha_regions( + AtomicNumber Z, UnprunedAtomicGridSpecification unp, double atom_radius ) { + + const auto rads = radial_points_for_spec(unp); + const auto& alphas = pyscf_alphas(Z); + std::vector regions(rads.size()); + + for( size_t i = 0; i < rads.size(); ++i ) { + size_t region = 0; + const auto scaled_rad = rads[i] / atom_radius; + for( const auto alpha : alphas ) + if( scaled_rad > alpha ) ++region; + regions[i] = region; + } + + return regions; +} + +int64_t previous_lebedev_npts(int64_t npts) { + const auto order = ll_traits::algebraic_order_by_npts(npts); + if( order < 0 ) + GAUXC_GENERIC_EXCEPTION("Unsupported PySCF Lebedev angular grid"); + + int64_t previous_order = -1; + auto current_order = ll_traits::next_algebraic_order(3); + while( current_order < order ) { + previous_order = current_order; + current_order = ll_traits::next_algebraic_order(current_order + 1); + } + + if( previous_order < 0 ) + GAUXC_GENERIC_EXCEPTION("Unsupported PySCF Lebedev angular grid"); + + return ll_traits::npts_by_algebraic_order(previous_order); +} + +double pyscf_sg1_radius(AtomicNumber Z) { + const auto radius = pyscf_gill_radius_93(Z); + if( radius <= 0. ) + GAUXC_GENERIC_EXCEPTION("Atomic number is outside the PySCF SG1 radius table"); + return radius; +} + +double pyscf_bragg_radius(AtomicNumber Z) { + const auto radius = pyscf_slater_radius_64(Z); + if( radius <= 0. ) + GAUXC_GENERIC_EXCEPTION("Atomic number is outside the PySCF Bragg radius table"); + return radius; +} + +} // namespace detail + /*****************/ /**** Visitor ****/ /*****************/ @@ -221,6 +336,119 @@ PrunedAtomicGridSpecification treutler_pruning_scheme( } +PrunedAtomicGridSpecification pyscf_sg1_pruning_scheme( + UnprunedAtomicGridSpecification unp, AtomicNumber Z ) { + + static const std::array leb_ngrid = {{6, 38, 86, 194, 86}}; + const auto regions = detail::pyscf_alpha_regions( + Z, unp, detail::pyscf_sg1_radius(Z) + ); + + std::vector angular_sizes; + angular_sizes.reserve(regions.size()); + for( const auto region : regions ) angular_sizes.push_back(leb_ngrid[region]); + + return detail::make_pruned_spec(unp, angular_sizes); +} + +PrunedAtomicGridSpecification pyscf_nwchem_pruning_scheme( + UnprunedAtomicGridSpecification unp, AtomicNumber Z ) { + + const auto n_ang = unp.angular_size.get(); + if( n_ang < 50 ) { + return detail::make_pruned_spec( + unp, std::vector(unp.radial_size.get(), n_ang) + ); + } + + std::array leb_npts; + if( n_ang == 50 ) { + leb_npts = {{ + detail::ll_traits::npts_by_algebraic_order(11), + detail::ll_traits::npts_by_algebraic_order(13), + detail::ll_traits::npts_by_algebraic_order(13), + detail::ll_traits::npts_by_algebraic_order(13), + detail::ll_traits::npts_by_algebraic_order(11) + }}; + } else { + const auto previous_n_ang = detail::previous_lebedev_npts(n_ang); + leb_npts = {{ + detail::ll_traits::npts_by_algebraic_order(11), + detail::ll_traits::npts_by_algebraic_order(15), + previous_n_ang, + n_ang, + previous_n_ang + }}; + } + + const auto regions = detail::pyscf_alpha_regions( + Z, unp, detail::pyscf_bragg_radius(Z) + ); + + std::vector angular_sizes; + angular_sizes.reserve(regions.size()); + for( const auto region : regions ) + angular_sizes.push_back(leb_npts[region]); + + return detail::make_pruned_spec(unp, angular_sizes); +} + +PrunedAtomicGridSpecification pyscf_sgx_pruning_scheme( + UnprunedAtomicGridSpecification unp, AtomicNumber Z ) { + + static const std::array, 8> sgx_ang_mapping = {{ + {{ 5, 7, 11, 11, 7}}, + {{ 5, 7, 11, 17, 11}}, + {{ 7, 11, 17, 23, 17}}, + {{ 7, 17, 23, 29, 23}}, + {{ 7, 23, 29, 35, 29}}, + {{11, 29, 35, 41, 35}}, + {{17, 35, 41, 47, 41}}, + {{47, 47, 47, 47, 47}} + }}; + + size_t level = 0; + const auto n_ang = unp.angular_size.get(); + for( const auto& mapping : sgx_ang_mapping ) { + const auto level_npts = detail::ll_traits::npts_by_algebraic_order(mapping[3]); + if( n_ang > level_npts ) ++level; + } + + if( level >= sgx_ang_mapping.size() ) + GAUXC_GENERIC_EXCEPTION("PySCF SGX pruning is not defined for this angular grid"); + + const auto regions = detail::pyscf_alpha_regions( + Z, unp, detail::pyscf_bragg_radius(Z) + ); + + std::vector angular_sizes; + angular_sizes.reserve(regions.size()); + for( const auto region : regions ) + angular_sizes.push_back( + detail::ll_traits::npts_by_algebraic_order(sgx_ang_mapping[level][region]) + ); + + return detail::make_pruned_spec(unp, angular_sizes); +} + +PrunedAtomicGridSpecification pyscf_treutler_pruning_scheme( + UnprunedAtomicGridSpecification unp ) { + + using angular_type = IntegratorXX::LebedevLaikov; + using traits = IntegratorXX::quadrature_traits; + + const size_t rsz = unp.radial_size.get(); + const size_t r_div_3 = rsz / 3ul; + const size_t r_div_2 = rsz / 2ul; + std::vector angular_sizes(rsz, unp.angular_size.get()); + std::fill(angular_sizes.begin(), angular_sizes.begin() + r_div_3, + traits::npts_by_algebraic_order(5)); + std::fill(angular_sizes.begin() + r_div_3, angular_sizes.begin() + r_div_2, + traits::npts_by_algebraic_order(11)); + + return detail::make_pruned_spec(unp, angular_sizes); +} + PrunedAtomicGridSpecification create_pruned_spec( PruningScheme scheme, UnprunedAtomicGridSpecification unp @@ -231,6 +459,14 @@ PrunedAtomicGridSpecification create_pruned_spec( return robust_psi4_pruning_scheme(unp); case PruningScheme::Treutler: return treutler_pruning_scheme(unp); + case PruningScheme::PySCF_Treutler: + return pyscf_treutler_pruning_scheme(unp); + case PruningScheme::PySCF_SG1: + GAUXC_GENERIC_EXCEPTION("PySCF SG1 pruning requires an atomic number"); + case PruningScheme::PySCF_NWChem: + GAUXC_GENERIC_EXCEPTION("PySCF NWChem pruning requires an atomic number"); + case PruningScheme::PySCF_SGX: + GAUXC_GENERIC_EXCEPTION("PySCF SGX pruning requires an atomic number"); // Default to Unpruned Grid case PruningScheme::Unpruned: @@ -245,4 +481,21 @@ PrunedAtomicGridSpecification create_pruned_spec( } +PrunedAtomicGridSpecification create_pruned_spec( + PruningScheme scheme, AtomicNumber Z, UnprunedAtomicGridSpecification unp +) { + + switch(scheme) { + case PruningScheme::PySCF_SG1: + return pyscf_sg1_pruning_scheme(unp, Z); + case PruningScheme::PySCF_NWChem: + return pyscf_nwchem_pruning_scheme(unp, Z); + case PruningScheme::PySCF_SGX: + return pyscf_sgx_pruning_scheme(unp, Z); + default: + return create_pruned_spec(scheme, unp); + } + +} + } diff --git a/src/molgrid_defaults.cxx b/src/molgrid_defaults.cxx index 61b66830f..a40858e06 100644 --- a/src/molgrid_defaults.cxx +++ b/src/molgrid_defaults.cxx @@ -15,6 +15,40 @@ namespace GauXC { +namespace detail { + +inline RadialSize get_pyscf_radial_size( AtomicNumber Z, int level ) { + if (level < 0 || level > 8) + GAUXC_GENERIC_EXCEPTION("Invalid PySCF grid level: " + std::to_string(level) + ". Valid levels are 0-8."); + + if ( Z.get() <= 2 ) { + if (level == 0) return RadialSize( 10 ); + return RadialSize( 20 + 10*level ); // ..., 30, 40, 50, ... + } else if ( Z.get() <= 10 ) { + if (level == 0) return RadialSize( 15 ); + if (level == 1) return RadialSize( 40 ); + return RadialSize( 30 + 15 * level ); // ..., 60, 75, 90, ... + } else if ( Z.get() <= 18 ) { + if (level == 0) return RadialSize( 20 ); + return RadialSize( 35 + 15 * level ); // ..., 50, 65, 80, ... + } else if ( Z.get() <= 36 ) { + if (level == 0) return RadialSize( 30 ); + return RadialSize( 45 + 15 * level ); // ..., 60, 75, 90, ... + } else if ( Z.get() <= 54 ) { + if (level == 0) return RadialSize( 35 ); + return RadialSize( 50 + 15 * level ); // ..., 65, 80, 95, ... + } else if ( Z.get() <= 86 ) { + if (level == 0) return RadialSize( 40 ); + return RadialSize( 55 + 15 * level ); // ..., 70, 85, 100, ... + } else if ( Z.get() <= 118 ) { + if (level == 0) return RadialSize( 50 ); + return RadialSize( 60 + 15 * level ); // ..., 75, 90, 105, ... + } else { + GAUXC_GENERIC_EXCEPTION("Z > 118 Not Supported for PySCF Grid Defaults"); + } +} +} + RadialScale default_mk_radial_scaling_factor( AtomicNumber _Z ) { auto Z = _Z.get(); switch(Z) { @@ -117,6 +151,36 @@ std::tuple case AtomicGridSizeDefault::GM5: return std::make_tuple( RadialSize(50), AngularSize(302) ); + case AtomicGridSizeDefault::PySCF0: + return std::make_tuple( detail::get_pyscf_radial_size(Z, 0), Z.get() <= 2 ? AngularSize(50) : (Z.get() <= 10 ? AngularSize(86) : AngularSize(110))); + + case AtomicGridSizeDefault::PySCF1: + return std::make_tuple( detail::get_pyscf_radial_size(Z, 1), Z.get() <= 2 ? AngularSize(110) : AngularSize(194)); + + case AtomicGridSizeDefault::PySCF2: + return std::make_tuple( detail::get_pyscf_radial_size(Z, 2), Z.get() <= 2 ? AngularSize(194) : AngularSize(302)); + + case AtomicGridSizeDefault::PySCF3: + return std::make_tuple( detail::get_pyscf_radial_size(Z, 3), Z.get() <= 10 ? AngularSize(302) : AngularSize(434)); + + case AtomicGridSizeDefault::PySCF4: + return std::make_tuple( detail::get_pyscf_radial_size(Z, 4), Z.get() <= 2 ? AngularSize(434) : AngularSize(590)); + + case AtomicGridSizeDefault::PySCF5: + return std::make_tuple( detail::get_pyscf_radial_size(Z, 5), Z.get() <= 2 ? AngularSize(590) : AngularSize(770)); + + case AtomicGridSizeDefault::PySCF6: + return std::make_tuple( detail::get_pyscf_radial_size(Z, 6), Z.get() <= 2 ? AngularSize(770) : AngularSize(974)); + + case AtomicGridSizeDefault::PySCF7: + return std::make_tuple( detail::get_pyscf_radial_size(Z, 7), Z.get() <= 2 ? AngularSize(974) : AngularSize(1202)); + + case AtomicGridSizeDefault::PySCF8: + return std::make_tuple( detail::get_pyscf_radial_size(Z, 8), AngularSize(1202)); + + case AtomicGridSizeDefault::PySCF9: + return std::make_tuple( RadialSize(200), AngularSize(1454) ); + case AtomicGridSizeDefault::FineGrid: return std::make_tuple( RadialSize(75), AngularSize(302) ); diff --git a/tests/c_api_test.cxx b/tests/c_api_test.cxx index e91941a42..793c41c22 100644 --- a/tests/c_api_test.cxx +++ b/tests/c_api_test.cxx @@ -126,6 +126,26 @@ TEST_CASE("C-API Enum Correspondence", "[c-api]") { == static_cast(GauXC_AtomicGridSizeDefault_GM3)); CHECK(static_cast(AtomicGridSizeDefault::GM5) == static_cast(GauXC_AtomicGridSizeDefault_GM5)); + CHECK(static_cast(AtomicGridSizeDefault::PySCF0) + == static_cast(GauXC_AtomicGridSizeDefault_PySCF0)); + CHECK(static_cast(AtomicGridSizeDefault::PySCF1) + == static_cast(GauXC_AtomicGridSizeDefault_PySCF1)); + CHECK(static_cast(AtomicGridSizeDefault::PySCF2) + == static_cast(GauXC_AtomicGridSizeDefault_PySCF2)); + CHECK(static_cast(AtomicGridSizeDefault::PySCF3) + == static_cast(GauXC_AtomicGridSizeDefault_PySCF3)); + CHECK(static_cast(AtomicGridSizeDefault::PySCF4) + == static_cast(GauXC_AtomicGridSizeDefault_PySCF4)); + CHECK(static_cast(AtomicGridSizeDefault::PySCF5) + == static_cast(GauXC_AtomicGridSizeDefault_PySCF5)); + CHECK(static_cast(AtomicGridSizeDefault::PySCF6) + == static_cast(GauXC_AtomicGridSizeDefault_PySCF6)); + CHECK(static_cast(AtomicGridSizeDefault::PySCF7) + == static_cast(GauXC_AtomicGridSizeDefault_PySCF7)); + CHECK(static_cast(AtomicGridSizeDefault::PySCF8) + == static_cast(GauXC_AtomicGridSizeDefault_PySCF8)); + CHECK(static_cast(AtomicGridSizeDefault::PySCF9) + == static_cast(GauXC_AtomicGridSizeDefault_PySCF9)); } // --- XCWeightAlg --- @@ -166,6 +186,14 @@ TEST_CASE("C-API Enum Correspondence", "[c-api]") { == static_cast(GauXC_PruningScheme_Robust)); CHECK(static_cast(PruningScheme::Treutler) == static_cast(GauXC_PruningScheme_Treutler)); + CHECK(static_cast(PruningScheme::PySCF_Treutler) + == static_cast(GauXC_PruningScheme_PySCF_Treutler)); + CHECK(static_cast(PruningScheme::PySCF_SG1) + == static_cast(GauXC_PruningScheme_PySCF_SG1)); + CHECK(static_cast(PruningScheme::PySCF_NWChem) + == static_cast(GauXC_PruningScheme_PySCF_NWChem)); + CHECK(static_cast(PruningScheme::PySCF_SGX) + == static_cast(GauXC_PruningScheme_PySCF_SGX)); } } diff --git a/tests/molgrid_test.cxx b/tests/molgrid_test.cxx index 1de1be986..8cd7b78e0 100644 --- a/tests/molgrid_test.cxx +++ b/tests/molgrid_test.cxx @@ -10,6 +10,7 @@ * See LICENSE.txt for details */ #include "catch2/catch.hpp" +#include #include #include @@ -219,6 +220,258 @@ TEST_CASE("MolGrid Defaults", "[molgrid]") { } } + SECTION("PySCF0 Grid Size") { + int64_t refgr[] = { -1, + 10, 10 , + 15, 15, 15, 15, 15, 15, 15, 15 , + 20, 20, 20, 20, 20, 20, 20, 20, + 30, 30, 30, 30, 30, 30, 30, 30, 30, 30, 30, 30, 30, 30, 30, 30, 30, 30, 35, 35, + 35, 35, 35, 35, 35, 35, 35, 35, 35, 35, 35, 35, 35, 35, 35, 35, + 40, 40, 40, 40, 40, 40, 40, 40, 40, 40, 40, 40, 40, 40, 40, 40, 40, 40, 40, 40, 40, 40, 40, 40, 40, 40, 40, 40, 40, 40, 40, 40, + 50, 50, 50, 50, 50, 50, 50, 50, 50, 50, 50, 50, 50, 50, 50, 50, 50, 50, 50, 50, 50, 50, 50, 50, 50, 50, 50, 50, 50, 50, 50, 50 + }; + int64_t refga[] = { -1, + 50, 50 , + 86, 86, 86, 86, 86, 86, 86, 86 , + 110, 110, 110, 110, 110, 110, 110, 110, + 110, 110, 110, 110, 110, 110, 110, 110, 110, 110, 110, 110, 110, 110, 110, 110, 110, 110, 110, 110, + 110, 110, 110, 110, 110, 110, 110, 110, 110, 110, 110, 110, 110, 110, 110, 110, + 110, 110, 110, 110, 110, 110, 110, 110, 110, 110, 110, 110, 110, 110, 110, 110, 110, 110, 110, 110, 110, 110, 110, 110, 110, 110, 110, 110, 110, 110, 110, 110, + 110, 110, 110, 110, 110, 110, 110, 110, 110, 110, 110, 110, 110, 110, 110, 110, 110, 110, 110, 110, 110, 110, 110, 110, 110, 110, 110, 110, 110, 110, 110, 110 + }; + for(auto Z = 1; Z < 119; ++Z) { + auto [gr, ga] = default_grid_size( AtomicNumber(Z), RadialQuad::MuraKnowles, + AtomicGridSizeDefault::PySCF0 ); + REQUIRE( gr.get() == refgr[Z] ); + REQUIRE( ga.get() == refga[Z] ); + } + } + + SECTION("PySCF1 Grid Size") { + int64_t refgr[] = { -1, + 30, 30 , + 40, 40, 40, 40, 40, 40, 40, 40 , + 50, 50, 50, 50, 50, 50, 50, 50, + 60, 60, 60, 60, 60, 60, 60, 60, 60, 60, 60, 60, 60, 60, 60, 60, 60, 60, 65, 65, + 65, 65, 65, 65, 65, 65, 65, 65, 65, 65, 65, 65, 65, 65, 65, 65, + 70, 70, 70, 70, 70, 70, 70, 70, 70, 70, 70, 70, 70, 70, 70, 70, 70, 70, 70, 70, 70, 70, 70, 70, 70, 70, 70, 70, 70, 70, 70, 70, + 75, 75, 75, 75, 75, 75, 75, 75, 75, 75, 75, 75, 75, 75, 75, 75, 75, 75, 75, 75, 75, 75, 75, 75, 75, 75, 75, 75, 75, 75, 75, 75 + }; + int64_t refga[] = { -1, + 110, 110 , + 194, 194, 194, 194, 194, 194, 194, 194 , + 194, 194, 194, 194, 194, 194, 194, 194, + 194, 194, 194, 194, 194, 194, 194, 194, 194, 194, 194, 194, 194, 194, 194, 194, 194, 194, 194, 194, + 194, 194, 194, 194, 194, 194, 194, 194, 194, 194, 194, 194, 194, 194, 194, 194, + 194, 194, 194, 194, 194, 194, 194, 194, 194, 194, 194, 194, 194, 194, 194, 194, 194, 194, 194, 194, 194, 194, 194, 194, 194, 194, 194, 194, 194, 194, 194, 194, + 194, 194, 194, 194, 194, 194, 194, 194, 194, 194, 194, 194, 194, 194, 194, 194, 194, 194, 194, 194, 194, 194, 194, 194, 194, 194, 194, 194, 194, 194, 194, 194 + }; + for(auto Z = 1; Z < 119; ++Z) { + auto [gr, ga] = default_grid_size( AtomicNumber(Z), RadialQuad::MuraKnowles, + AtomicGridSizeDefault::PySCF1 ); + REQUIRE( gr.get() == refgr[Z] ); + REQUIRE( ga.get() == refga[Z] ); + } + } + + SECTION("PySCF2 Grid Size") { + int64_t refgr[] = { -1, + 40, 40 , + 60, 60, 60, 60, 60, 60, 60, 60 , + 65, 65, 65, 65, 65, 65, 65, 65, + 75, 75, 75, 75, 75, 75, 75, 75, 75, 75, 75, 75, 75, 75, 75, 75, 75, 75, 80, 80, + 80, 80, 80, 80, 80, 80, 80, 80, 80, 80, 80, 80, 80, 80, 80, 80, + 85, 85, 85, 85, 85, 85, 85, 85, 85, 85, 85, 85, 85, 85, 85, 85, 85, 85, 85, 85, 85, 85, 85, 85, 85, 85, 85, 85, 85, 85, 85, 85, + 90, 90, 90, 90, 90, 90, 90, 90, 90, 90, 90, 90, 90, 90, 90, 90, 90, 90, 90, 90, 90, 90, 90, 90, 90, 90, 90, 90, 90, 90, 90, 90 + }; + int64_t refga[] = { -1, + 194, 194 , + 302, 302, 302, 302, 302, 302, 302, 302 , + 302, 302, 302, 302, 302, 302, 302, 302, + 302, 302, 302, 302, 302, 302, 302, 302, 302, 302, 302, 302, 302, 302, 302, 302, 302, 302, 302, 302, + 302, 302, 302, 302, 302, 302, 302, 302, 302, 302, 302, 302, 302, 302, 302, 302, + 302, 302, 302, 302, 302, 302, 302, 302, 302, 302, 302, 302, 302, 302, 302, 302, 302, 302, 302, 302, 302, 302, 302, 302, 302, 302, 302, 302, 302, 302, 302, 302, + 302, 302, 302, 302, 302, 302, 302, 302, 302, 302, 302, 302, 302, 302, 302, 302, 302, 302, 302, 302, 302, 302, 302, 302, 302, 302, 302, 302, 302, 302, 302, 302 + }; + for(auto Z = 1; Z < 119; ++Z) { + auto [gr, ga] = default_grid_size( AtomicNumber(Z), RadialQuad::MuraKnowles, + AtomicGridSizeDefault::PySCF2 ); + REQUIRE( gr.get() == refgr[Z] ); + REQUIRE( ga.get() == refga[Z] ); + } + } + + SECTION("PySCF3 Grid Size") { + int64_t refgr[] = { -1, + 50, 50 , + 75, 75, 75, 75, 75, 75, 75, 75 , + 80, 80, 80, 80, 80, 80, 80, 80, + 90, 90, 90, 90, 90, 90, 90, 90, 90, 90, 90, 90, 90, 90, 90, 90, 90, 90, 95, 95, + 95, 95, 95, 95, 95, 95, 95, 95, 95, 95, 95, 95, 95, 95, 95, 95, + 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, + 105, 105, 105, 105, 105, 105, 105, 105, 105, 105, 105, 105, 105, 105, 105, 105, 105, 105, 105, 105, 105, 105, 105, 105, 105, 105, 105, 105, 105, 105, 105, 105 + }; + int64_t refga[] = { -1, + 302, 302 , + 302, 302, 302, 302, 302, 302, 302, 302 , + 434, 434, 434, 434, 434, 434, 434, 434, + 434, 434, 434, 434, 434, 434, 434, 434, 434, 434, 434, 434, 434, 434, 434, 434, 434, 434, 434, 434, + 434, 434, 434, 434, 434, 434, 434, 434, 434, 434, 434, 434, 434, 434, 434, 434, + 434, 434, 434, 434, 434, 434, 434, 434, 434, 434, 434, 434, 434, 434, 434, 434, 434, 434, 434, 434, 434, 434, 434, 434, 434, 434, 434, 434, 434, 434, 434, 434, + 434, 434, 434, 434, 434, 434, 434, 434, 434, 434, 434, 434, 434, 434, 434, 434, 434, 434, 434, 434, 434, 434, 434, 434, 434, 434, 434, 434, 434, 434, 434, 434 + }; + for(auto Z = 1; Z < 119; ++Z) { + auto [gr, ga] = default_grid_size( AtomicNumber(Z), RadialQuad::MuraKnowles, + AtomicGridSizeDefault::PySCF3 ); + REQUIRE( gr.get() == refgr[Z] ); + REQUIRE( ga.get() == refga[Z] ); + } + } + + SECTION("PySCF4 Grid Size") { + int64_t refgr[] = { -1, + 60, 60 , + 90, 90, 90, 90, 90, 90, 90, 90 , + 95, 95, 95, 95, 95, 95, 95, 95, + 105, 105, 105, 105, 105, 105, 105, 105, 105, 105, 105, 105, 105, 105, 105, 105, 105, 105, 110, 110, + 110, 110, 110, 110, 110, 110, 110, 110, 110, 110, 110, 110, 110, 110, 110, 110, + 115, 115, 115, 115, 115, 115, 115, 115, 115, 115, 115, 115, 115, 115, 115, 115, 115, 115, 115, 115, 115, 115, 115, 115, 115, 115, 115, 115, 115, 115, 115, 115, + 120, 120, 120, 120, 120, 120, 120, 120, 120, 120, 120, 120, 120, 120, 120, 120, 120, 120, 120, 120, 120, 120, 120, 120, 120, 120, 120, 120, 120, 120, 120, 120 + }; + int64_t refga[] = { -1, + 434, 434 , + 590, 590, 590, 590, 590, 590, 590, 590 , + 590, 590, 590, 590, 590, 590, 590, 590, + 590, 590, 590, 590, 590, 590, 590, 590, 590, 590, 590, 590, 590, 590, 590, 590, 590, 590, 590, 590, + 590, 590, 590, 590, 590, 590, 590, 590, 590, 590, 590, 590, 590, 590, 590, 590, + 590, 590, 590, 590, 590, 590, 590, 590, 590, 590, 590, 590, 590, 590, 590, 590, 590, 590, 590, 590, 590, 590, 590, 590, 590, 590, 590, 590, 590, 590, 590, 590, + 590, 590, 590, 590, 590, 590, 590, 590, 590, 590, 590, 590, 590, 590, 590, 590, 590, 590, 590, 590, 590, 590, 590, 590, 590, 590, 590, 590, 590, 590, 590, 590 + }; + for(auto Z = 1; Z < 119; ++Z) { + auto [gr, ga] = default_grid_size( AtomicNumber(Z), RadialQuad::MuraKnowles, + AtomicGridSizeDefault::PySCF4 ); + REQUIRE( gr.get() == refgr[Z] ); + REQUIRE( ga.get() == refga[Z] ); + } + } + + SECTION("PySCF5 Grid Size") { + int64_t refgr[] = { -1, + 70, 70 , + 105, 105, 105, 105, 105, 105, 105, 105 , + 110, 110, 110, 110, 110, 110, 110, 110, + 120, 120, 120, 120, 120, 120, 120, 120, 120, 120, 120, 120, 120, 120, 120, 120, 120, 120, 125, 125, + 125, 125, 125, 125, 125, 125, 125, 125, 125, 125, 125, 125, 125, 125, 125, 125, + 130, 130, 130, 130, 130, 130, 130, 130, 130, 130, 130, 130, 130, 130, 130, 130, 130, 130, 130, 130, 130, 130, 130, 130, 130, 130, 130, 130, 130, 130, 130, 130, + 135, 135, 135, 135, 135, 135, 135, 135, 135, 135, 135, 135, 135, 135, 135, 135, 135, 135, 135, 135, 135, 135, 135, 135, 135, 135, 135, 135, 135, 135, 135, 135 + }; + int64_t refga[] = { -1, + 590, 590 , + 770, 770, 770, 770, 770, 770, 770, 770 , + 770, 770, 770, 770, 770, 770, 770, 770, + 770, 770, 770, 770, 770, 770, 770, 770, 770, 770, 770, 770, 770, 770, 770, 770, 770, 770, 770, 770, + 770, 770, 770, 770, 770, 770, 770, 770, 770, 770, 770, 770, 770, 770, 770, 770, + 770, 770, 770, 770, 770, 770, 770, 770, 770, 770, 770, 770, 770, 770, 770, 770, 770, 770, 770, 770, 770, 770, 770, 770, 770, 770, 770, 770, 770, 770, 770, 770, + 770, 770, 770, 770, 770, 770, 770, 770, 770, 770, 770, 770, 770, 770, 770, 770, 770, 770, 770, 770, 770, 770, 770, 770, 770, 770, 770, 770, 770, 770, 770, 770 + }; + for(auto Z = 1; Z < 119; ++Z) { + auto [gr, ga] = default_grid_size( AtomicNumber(Z), RadialQuad::MuraKnowles, + AtomicGridSizeDefault::PySCF5 ); + REQUIRE( gr.get() == refgr[Z] ); + REQUIRE( ga.get() == refga[Z] ); + } + } + + SECTION("PySCF6 Grid Size") { + int64_t refgr[] = { -1, + 80, 80 , + 120, 120, 120, 120, 120, 120, 120, 120 , + 125, 125, 125, 125, 125, 125, 125, 125, + 135, 135, 135, 135, 135, 135, 135, 135, 135, 135, 135, 135, 135, 135, 135, 135, 135, 135, 140, 140, + 140, 140, 140, 140, 140, 140, 140, 140, 140, 140, 140, 140, 140, 140, 140, 140, + 145, 145, 145, 145, 145, 145, 145, 145, 145, 145, 145, 145, 145, 145, 145, 145, 145, 145, 145, 145, 145, 145, 145, 145, 145, 145, 145, 145, 145, 145, 145, 145, + 150, 150, 150, 150, 150, 150, 150, 150, 150, 150, 150, 150, 150, 150, 150, 150, 150, 150, 150, 150, 150, 150, 150, 150, 150, 150, 150, 150, 150, 150, 150, 150 + }; + int64_t refga[] = { -1, + 770, 770 , + 974, 974, 974, 974, 974, 974, 974, 974 , + 974, 974, 974, 974, 974, 974, 974, 974, + 974, 974, 974, 974, 974, 974, 974, 974, 974, 974, 974, 974, 974, 974, 974, 974, 974, 974, 974, 974, + 974, 974, 974, 974, 974, 974, 974, 974, 974, 974, 974, 974, 974, 974, 974, 974, + 974, 974, 974, 974, 974, 974, 974, 974, 974, 974, 974, 974, 974, 974, 974, 974, 974, 974, 974, 974, 974, 974, 974, 974, 974, 974, 974, 974, 974, 974, 974, 974, + 974, 974, 974, 974, 974, 974, 974, 974, 974, 974, 974, 974, 974, 974, 974, 974, 974, 974, 974, 974, 974, 974, 974, 974, 974, 974, 974, 974, 974, 974, 974, 974 + }; + for(auto Z = 1; Z < 119; ++Z) { + auto [gr, ga] = default_grid_size( AtomicNumber(Z), RadialQuad::MuraKnowles, + AtomicGridSizeDefault::PySCF6 ); + REQUIRE( gr.get() == refgr[Z] ); + REQUIRE( ga.get() == refga[Z] ); + } + } + + SECTION("PySCF7 Grid Size") { + int64_t refgr[] = { -1, + 90, 90 , + 135, 135, 135, 135, 135, 135, 135, 135 , + 140, 140, 140, 140, 140, 140, 140, 140, + 150, 150, 150, 150, 150, 150, 150, 150, 150, 150, 150, 150, 150, 150, 150, 150, 150, 150, 155, 155, + 155, 155, 155, 155, 155, 155, 155, 155, 155, 155, 155, 155, 155, 155, 155, 155, + 160, 160, 160, 160, 160, 160, 160, 160, 160, 160, 160, 160, 160, 160, 160, 160, 160, 160, 160, 160, 160, 160, 160, 160, 160, 160, 160, 160, 160, 160, 160, 160, + 165, 165, 165, 165, 165, 165, 165, 165, 165, 165, 165, 165, 165, 165, 165, 165, 165, 165, 165, 165, 165, 165, 165, 165, 165, 165, 165, 165, 165, 165, 165, 165 + }; + int64_t refga[] = { -1, + 974, 974 , + 1202, 1202, 1202, 1202, 1202, 1202, 1202, 1202 , + 1202, 1202, 1202, 1202, 1202, 1202, 1202, 1202, + 1202, 1202, 1202, 1202, 1202, 1202, 1202, 1202, 1202, 1202, 1202, 1202, 1202, 1202, 1202, 1202, 1202, 1202, 1202, 1202, + 1202, 1202, 1202, 1202, 1202, 1202, 1202, 1202, 1202, 1202, 1202, 1202, 1202, 1202, 1202, 1202, + 1202, 1202, 1202, 1202, 1202, 1202, 1202, 1202, 1202, 1202, 1202, 1202, 1202, 1202, 1202, 1202, 1202, 1202, 1202, 1202, 1202, 1202, 1202, 1202, 1202, 1202, 1202, 1202, 1202, 1202, 1202, 1202, + 1202, 1202, 1202, 1202, 1202, 1202, 1202, 1202, 1202, 1202, 1202, 1202, 1202, 1202, 1202, 1202, 1202, 1202, 1202, 1202, 1202, 1202, 1202, 1202, 1202, 1202, 1202, 1202, 1202, 1202, 1202, 1202 + }; + for(auto Z = 1; Z < 119; ++Z) { + auto [gr, ga] = default_grid_size( AtomicNumber(Z), RadialQuad::MuraKnowles, + AtomicGridSizeDefault::PySCF7 ); + REQUIRE( gr.get() == refgr[Z] ); + REQUIRE( ga.get() == refga[Z] ); + } + } + + SECTION("PySCF8 Grid Size") { + int64_t refgr[] = { -1, + 100, 100 , + 150, 150, 150, 150, 150, 150, 150, 150 , + 155, 155, 155, 155, 155, 155, 155, 155, + 165, 165, 165, 165, 165, 165, 165, 165, 165, 165, 165, 165, 165, 165, 165, 165, 165, 165, 170, 170, + 170, 170, 170, 170, 170, 170, 170, 170, 170, 170, 170, 170, 170, 170, 170, 170, + 175, 175, 175, 175, 175, 175, 175, 175, 175, 175, 175, 175, 175, 175, 175, 175, 175, 175, 175, 175, 175, 175, 175, 175, 175, 175, 175, 175, 175, 175, 175, 175, + 180, 180, 180, 180, 180, 180, 180, 180, 180, 180, 180, 180, 180, 180, 180, 180, 180, 180, 180, 180, 180, 180, 180, 180, 180, 180, 180, 180, 180, 180, 180, 180 + }; + int64_t refga[] = { -1, + 1202, 1202 , + 1202, 1202, 1202, 1202, 1202, 1202, 1202, 1202 , + 1202, 1202, 1202, 1202, 1202, 1202, 1202, 1202, + 1202, 1202, 1202, 1202, 1202, 1202, 1202, 1202, 1202, 1202, 1202, 1202, 1202, 1202, 1202, 1202, 1202, 1202, 1202, 1202, + 1202, 1202, 1202, 1202, 1202, 1202, 1202, 1202, 1202, 1202, 1202, 1202, 1202, 1202, 1202, 1202, + 1202, 1202, 1202, 1202, 1202, 1202, 1202, 1202, 1202, 1202, 1202, 1202, 1202, 1202, 1202, 1202, 1202, 1202, 1202, 1202, 1202, 1202, 1202, 1202, 1202, 1202, 1202, 1202, 1202, 1202, 1202, 1202, + 1202, 1202, 1202, 1202, 1202, 1202, 1202, 1202, 1202, 1202, 1202, 1202, 1202, 1202, 1202, 1202, 1202, 1202, 1202, 1202, 1202, 1202, 1202, 1202, 1202, 1202, 1202, 1202, 1202, 1202, 1202, 1202 + }; + for(auto Z = 1; Z < 119; ++Z) { + auto [gr, ga] = default_grid_size( AtomicNumber(Z), RadialQuad::MuraKnowles, + AtomicGridSizeDefault::PySCF8 ); + REQUIRE( gr.get() == refgr[Z] ); + REQUIRE( ga.get() == refga[Z] ); + } + } + + SECTION("PySCF9 Grid Size") { + for(auto Z = 1; Z < 119; ++Z) { + auto [gr, ga] = default_grid_size( AtomicNumber(Z), RadialQuad::MuraKnowles, + AtomicGridSizeDefault::PySCF9 ); + REQUIRE( gr.get() == 200 ); + REQUIRE( ga.get() == 1454 ); + } + } + } @@ -285,6 +538,19 @@ TEST_CASE("Grid Specification", "[molgrid]") { }; } + SECTION("PySCF Treutler") { + gs = MolGridFactory::create_default_pruned_grid_spec( + PruningScheme::PySCF_Treutler,Z,rq,gsz); + size_t rs = rsz.get(); + size_t r3 = rsz.get() / 3ul; + size_t r2 = rsz.get() / 2ul; + ref_pruning_regions = { + {0ul, r3, AngularSize(14)}, + {r3, r2, AngularSize(50)}, + {r2, rs, AngularSize(590)} + }; + } + std::visit([=](auto& g) { REQUIRE( g.radial_quad == rq ); REQUIRE( g.radial_size == rsz ); @@ -298,6 +564,209 @@ TEST_CASE("Grid Specification", "[molgrid]") { } + SECTION("PySCF Slater-Bragg Radii") { + + const std::vector> pyscf_slater_reference = { + { 1, 6.61404143597771554e-01}, + { 2, 2.64561657439108622e+00}, + { 6, 1.32280828719554311e+00}, + { 10, 2.83458918684759276e+00}, + { 18, 3.40150702421711149e+00}, + { 36, 3.59047963667361714e+00}, + { 54, 3.96842486158663021e+00}, + { 86, 3.96842486158663021e+00}, + {103, 3.30702071798885822e+00}, + {130, 3.30702071798885822e+00}, + }; + + for( const auto& [Z, radius] : pyscf_slater_reference ) { + CHECK( pyscf_slater_radius_64(AtomicNumber(Z)) == Approx(radius) ); + } + + } + +} + + +TEST_CASE("PySCF Pruning Specifications", "[molgrid]") { + + struct PySCFPruningCase { + PruningScheme scheme; + AtomicNumber Z; + AtomicGridSizeDefault grid_size; + std::vector reference_regions; + }; + + const std::vector test_cases = { + {PruningScheme::PySCF_Treutler, AtomicNumber(1), AtomicGridSizeDefault::PySCF0, + {{0ul, 3ul, AngularSize(14)}, {3ul, 10ul, AngularSize(50)}}}, + {PruningScheme::PySCF_Treutler, AtomicNumber(11), AtomicGridSizeDefault::PySCF3, + {{0ul, 26ul, AngularSize(14)}, {26ul, 40ul, AngularSize(50)}, + {40ul, 80ul, AngularSize(434)}}}, + {PruningScheme::PySCF_Treutler, AtomicNumber(8), AtomicGridSizeDefault::PySCF2, + {{0ul, 20ul, AngularSize(14)}, {20ul, 30ul, AngularSize(50)}, + {30ul, 60ul, AngularSize(302)}}}, + {PruningScheme::PySCF_Treutler, AtomicNumber(36), AtomicGridSizeDefault::PySCF4, + {{0ul, 35ul, AngularSize(14)}, {35ul, 52ul, AngularSize(50)}, + {52ul, 105ul, AngularSize(590)}}}, + + {PruningScheme::PySCF_SG1, AtomicNumber(1), AtomicGridSizeDefault::PySCF0, + {{0ul, 3ul, AngularSize(6)}, {3ul, 4ul, AngularSize(38)}, + {4ul, 5ul, AngularSize(86)}, {5ul, 9ul, AngularSize(194)}, + {9ul, 10ul, AngularSize(86)}}}, + {PruningScheme::PySCF_SG1, AtomicNumber(2), AtomicGridSizeDefault::PySCF4, + {{0ul, 17ul, AngularSize(6)}, {17ul, 21ul, AngularSize(38)}, + {21ul, 26ul, AngularSize(86)}, {26ul, 42ul, AngularSize(194)}, + {42ul, 60ul, AngularSize(86)}}}, + {PruningScheme::PySCF_SG1, AtomicNumber(6), AtomicGridSizeDefault::PySCF3, + {{0ul, 22ul, AngularSize(6)}, {22ul, 31ul, AngularSize(38)}, + {31ul, 38ul, AngularSize(86)}, {38ul, 57ul, AngularSize(194)}, + {57ul, 75ul, AngularSize(86)}}}, + {PruningScheme::PySCF_SG1, AtomicNumber(8), AtomicGridSizeDefault::PySCF2, + {{0ul, 17ul, AngularSize(6)}, {17ul, 24ul, AngularSize(38)}, + {24ul, 29ul, AngularSize(86)}, {29ul, 44ul, AngularSize(194)}, + {44ul, 60ul, AngularSize(86)}}}, + {PruningScheme::PySCF_SG1, AtomicNumber(11), AtomicGridSizeDefault::PySCF1, + {{0ul, 17ul, AngularSize(6)}, {17ul, 26ul, AngularSize(38)}, + {26ul, 33ul, AngularSize(86)}, {33ul, 45ul, AngularSize(194)}, + {45ul, 50ul, AngularSize(86)}}}, + {PruningScheme::PySCF_SG1, AtomicNumber(18), AtomicGridSizeDefault::PySCF4, + {{0ul, 25ul, AngularSize(6)}, {25ul, 39ul, AngularSize(38)}, + {39ul, 48ul, AngularSize(86)}, {48ul, 69ul, AngularSize(194)}, + {69ul, 95ul, AngularSize(86)}}}, + + {PruningScheme::PySCF_NWChem, AtomicNumber(1), AtomicGridSizeDefault::PySCF0, + {{0ul, 3ul, AngularSize(50)}, {3ul, 8ul, AngularSize(74)}, + {8ul, 10ul, AngularSize(50)}}}, + {PruningScheme::PySCF_NWChem, AtomicNumber(2), AtomicGridSizeDefault::PySCF4, + {{0ul, 27ul, AngularSize(50)}, {27ul, 34ul, AngularSize(86)}, + {34ul, 42ul, AngularSize(350)}, {42ul, 59ul, AngularSize(434)}, + {59ul, 60ul, AngularSize(350)}}}, + {PruningScheme::PySCF_NWChem, AtomicNumber(6), AtomicGridSizeDefault::PySCF3, + {{0ul, 22ul, AngularSize(50)}, {22ul, 32ul, AngularSize(86)}, + {32ul, 38ul, AngularSize(266)}, {38ul, 58ul, AngularSize(302)}, + {58ul, 75ul, AngularSize(266)}}}, + {PruningScheme::PySCF_NWChem, AtomicNumber(8), AtomicGridSizeDefault::PySCF2, + {{0ul, 18ul, AngularSize(50)}, {18ul, 26ul, AngularSize(86)}, + {26ul, 31ul, AngularSize(266)}, {31ul, 47ul, AngularSize(302)}, + {47ul, 60ul, AngularSize(266)}}}, + {PruningScheme::PySCF_NWChem, AtomicNumber(10), AtomicGridSizeDefault::PySCF3, + {{0ul, 31ul, AngularSize(50)}, {31ul, 43ul, AngularSize(86)}, + {43ul, 52ul, AngularSize(266)}, {52ul, 72ul, AngularSize(302)}, + {72ul, 75ul, AngularSize(266)}}}, + {PruningScheme::PySCF_NWChem, AtomicNumber(36), AtomicGridSizeDefault::PySCF4, + {{0ul, 39ul, AngularSize(50)}, {39ul, 61ul, AngularSize(86)}, + {61ul, 75ul, AngularSize(434)}, {75ul, 99ul, AngularSize(590)}, + {99ul, 105ul, AngularSize(434)}}}, + + {PruningScheme::PySCF_SGX, AtomicNumber(1), AtomicGridSizeDefault::PySCF5, + {{0ul, 21ul, AngularSize(50)}, {21ul, 26ul, AngularSize(302)}, + {26ul, 33ul, AngularSize(434)}, {33ul, 52ul, AngularSize(590)}, + {52ul, 70ul, AngularSize(434)}}}, + {PruningScheme::PySCF_SGX, AtomicNumber(2), AtomicGridSizeDefault::PySCF0, + {{0ul, 4ul, AngularSize(14)}, {4ul, 6ul, AngularSize(26)}, + {6ul, 10ul, AngularSize(50)}}}, + {PruningScheme::PySCF_SGX, AtomicNumber(6), AtomicGridSizeDefault::PySCF4, + {{0ul, 27ul, AngularSize(50)}, {27ul, 38ul, AngularSize(302)}, + {38ul, 46ul, AngularSize(434)}, {46ul, 70ul, AngularSize(590)}, + {70ul, 90ul, AngularSize(434)}}}, + {PruningScheme::PySCF_SGX, AtomicNumber(8), AtomicGridSizeDefault::PySCF2, + {{0ul, 18ul, AngularSize(26)}, {18ul, 26ul, AngularSize(110)}, + {26ul, 31ul, AngularSize(194)}, {31ul, 47ul, AngularSize(302)}, + {47ul, 60ul, AngularSize(194)}}}, + {PruningScheme::PySCF_SGX, AtomicNumber(10), AtomicGridSizeDefault::PySCF3, + {{0ul, 31ul, AngularSize(26)}, {31ul, 43ul, AngularSize(110)}, + {43ul, 52ul, AngularSize(194)}, {52ul, 72ul, AngularSize(302)}, + {72ul, 75ul, AngularSize(194)}}}, + {PruningScheme::PySCF_SGX, AtomicNumber(18), AtomicGridSizeDefault::PySCF4, + {{0ul, 34ul, AngularSize(50)}, {34ul, 52ul, AngularSize(302)}, + {52ul, 65ul, AngularSize(434)}, {65ul, 87ul, AngularSize(590)}, + {87ul, 95ul, AngularSize(434)}}} + }; + + for( const auto& test_case : test_cases ) { + auto gs = MolGridFactory::create_default_pruned_grid_spec( + test_case.scheme, test_case.Z, RadialQuad::TreutlerAhlrichs, + test_case.grid_size + ); + + const auto* pruned = std::get_if(&gs); + REQUIRE( pruned ); + CHECK( pruned->pruning_regions == test_case.reference_regions ); + } + +} + +TEST_CASE("PySCF Pruning Known Differences", "[molgrid]") { + + SECTION("PySCF Treutler radial grids extend past GauXC TA defaults") { + REQUIRE_THROWS( MolGridFactory::create_default_pruned_grid_spec( + PruningScheme::PySCF_NWChem, AtomicNumber(54), + RadialQuad::TreutlerAhlrichs, AtomicGridSizeDefault::PySCF4 + ) ); + } + +} + +TEST_CASE("PySCF Pruning Exceptions", "[molgrid]") { + + const UnprunedAtomicGridSpecification unp{ + RadialQuad::TreutlerAhlrichs, RadialSize(10), RadialScale(1.0), + AngularSize(302) + }; + + SECTION("Atom-aware PySCF schemes require an atomic number") { + REQUIRE_THROWS_AS( create_pruned_spec(PruningScheme::PySCF_SG1, unp), + generic_gauxc_exception ); + REQUIRE_THROWS_AS( create_pruned_spec(PruningScheme::PySCF_NWChem, unp), + generic_gauxc_exception ); + REQUIRE_THROWS_AS( create_pruned_spec(PruningScheme::PySCF_SGX, unp), + generic_gauxc_exception ); + } + + SECTION("SG1 radii are limited to H-Ar; PySCF raises IndexError for Z=19") { + REQUIRE_THROWS_AS( + create_pruned_spec(PruningScheme::PySCF_SG1, AtomicNumber(19), unp), + generic_gauxc_exception + ); + } + + SECTION("PySCF accepts Z=0 SG1 via its dummy radius; GauXC rejects it") { + REQUIRE_THROWS_AS( + create_pruned_spec(PruningScheme::PySCF_SG1, AtomicNumber(0), unp), + generic_gauxc_exception + ); + } + + SECTION("Bragg radii stop at Z=130; PySCF raises IndexError for Z=131") { + REQUIRE_THROWS_AS( + create_pruned_spec(PruningScheme::PySCF_NWChem, AtomicNumber(131), unp), + generic_gauxc_exception + ); + REQUIRE_THROWS_AS( + create_pruned_spec(PruningScheme::PySCF_SGX, AtomicNumber(131), unp), + generic_gauxc_exception + ); + } + + SECTION("NWChem requires a supported Lebedev angular size; PySCF raises IndexError") { + auto invalid_angular = unp; + invalid_angular.angular_size = AngularSize(52); + REQUIRE_THROWS_AS( + create_pruned_spec(PruningScheme::PySCF_NWChem, AtomicNumber(8), invalid_angular), + generic_gauxc_exception + ); + } + + SECTION("SGX is bounded by PySCF's mapping table; PySCF raises IndexError") { + auto too_large_sgx = unp; + too_large_sgx.angular_size = AngularSize(6000); + REQUIRE_THROWS_AS( + create_pruned_spec(PruningScheme::PySCF_SGX, AtomicNumber(8), too_large_sgx), + generic_gauxc_exception + ); + } + } diff --git a/tests/standalone_driver.cxx b/tests/standalone_driver.cxx index 68a9c13aa..35f206f60 100644 --- a/tests/standalone_driver.cxx +++ b/tests/standalone_driver.cxx @@ -187,9 +187,13 @@ int main(int argc, char** argv) { }; std::map< std::string, PruningScheme > prune_map = { - {"UNPRUNED", PruningScheme::Unpruned}, - {"ROBUST", PruningScheme::Robust}, - {"TREUTLER", PruningScheme::Treutler} + {"UNPRUNED", PruningScheme::Unpruned}, + {"ROBUST", PruningScheme::Robust}, + {"TREUTLER", PruningScheme::Treutler}, + {"PYSCF_TREUTLER", PruningScheme::PySCF_Treutler}, + {"PYSCF_SG1", PruningScheme::PySCF_SG1}, + {"PYSCF_NWCHEM", PruningScheme::PySCF_NWChem}, + {"PYSCF_SGX", PruningScheme::PySCF_SGX} }; std::map< std::string, RadialQuad > rad_quad_map = {