diff --git a/include/gauxc/c/enums.h b/include/gauxc/c/enums.h index 095b37737..28f1404e9 100644 --- a/include/gauxc/c/enums.h +++ b/include/gauxc/c/enums.h @@ -47,9 +47,10 @@ enum GauXC_AtomicGridSizeDefault { */ enum GauXC_XCWeightAlg { GauXC_XCWeightAlg_NOTPARTITIONED, ///< Not partitioned - GauXC_XCWeightAlg_Becke, ///< The original Becke weighting scheme - GauXC_XCWeightAlg_SSF, ///< The Stratmann-Scuseria-Frisch weighting scheme - GauXC_XCWeightAlg_LKO ///< The Lauqua-Kuessman-Ochsenfeld weighting scheme + GauXC_XCWeightAlg_Becke, ///< The original Becke weighting scheme + GauXC_XCWeightAlg_SSF, ///< The Stratmann-Scuseria-Frisch weighting scheme + GauXC_XCWeightAlg_LKO, ///< The Lauqua-Kuessman-Ochsenfeld weighting scheme + GauXC_XCWeightAlg_Hirshfeld ///< Hirshfeld pro-atom density weighting scheme }; /** diff --git a/include/gauxc/enums.hpp b/include/gauxc/enums.hpp index 9ea604564..f4e7b9b4f 100644 --- a/include/gauxc/enums.hpp +++ b/include/gauxc/enums.hpp @@ -46,7 +46,8 @@ enum class XCWeightAlg { NOTPARTITIONED, ///< Not partitioned Becke, ///< The original Becke weighting scheme SSF, ///< The Stratmann-Scuseria-Frisch weighting scheme - LKO ///< The Lauqua-Kuessman-Ochsenfeld weighting scheme + LKO, ///< The Lauqua-Kuessman-Ochsenfeld weighting scheme + Hirshfeld ///< Hirshfeld pro-atom density weighting scheme }; /** diff --git a/src/xc_integrator/local_work_driver/host/reference/weights.cxx b/src/xc_integrator/local_work_driver/host/reference/weights.cxx index 145bfd1cd..b209234ca 100644 --- a/src/xc_integrator/local_work_driver/host/reference/weights.cxx +++ b/src/xc_integrator/local_work_driver/host/reference/weights.cxx @@ -14,6 +14,8 @@ #include +#include + namespace GauXC { // Reference Becke weights impl @@ -385,6 +387,74 @@ void reference_lko_weights_host( } +void reference_hirshfeld_weights_host( + const Molecule& mol, + const MolMeta& meta, + task_iterator task_begin, + task_iterator task_end +) { + + (void)meta; + + constexpr double pi = 3.141592653589793238462643383279502884; + + const size_t ntasks = std::distance(task_begin,task_end); + const size_t natoms = mol.natoms(); + + struct hirshfeld_atom_data { + double x; + double y; + double z; + double alpha; + double scale; + }; + + std::vector atom_data; + atom_data.reserve(natoms); + + for( const auto& atom : mol ) { + const double radius = default_atomic_radius(atom.Z); + const double alpha = 0.5 / (radius * radius); + const double scale = atom.Z.get() * std::pow(alpha / pi, 1.5); + atom_data.push_back( hirshfeld_atom_data{ atom.x, atom.y, atom.z, alpha, scale } ); + } + + #pragma omp parallel + { + + #pragma omp for + for( size_t iT = 0; iT < ntasks; ++iT ) + for( size_t i = 0; i < (task_begin+iT)->points.size(); ++i ) { + + auto& task = *(task_begin+iT); + auto& weight = task.weights[i]; + const auto& point = task.points[i]; + const auto parent = static_cast(task.iParent); + + double sum = 0.; + double parent_density = 0.; + for(size_t iA = 0; iA < natoms; iA++) { + const auto& atom = atom_data[iA]; + + const double da_x = point[0] - atom.x; + const double da_y = point[1] - atom.y; + const double da_z = point[2] - atom.z; + const double r2 = da_x*da_x + da_y*da_y + da_z*da_z; + const double density = atom.scale * std::exp(-atom.alpha * r2); + + sum += density; + if( iA == parent ) parent_density = density; + + } + + weight *= (sum > 0.) ? parent_density / sum : 0.; + + } // Collapsed loop over tasks and points + + } // OMP context + +} + /** * 1st derivative which expects weight_deri to be preallocated as (ngrid*natoms*3) @@ -682,6 +752,81 @@ void reference_ssf_weights_1st_derivative_host( } } +void reference_hirshfeld_weights_1st_derivative_host( + const Molecule& mol, + const MolMeta& meta, + const XCTask& task, + double* weight_deri +){ + + (void)meta; + + constexpr double pi = 3.141592653589793238462643383279502884; + + const size_t natoms = mol.natoms(); + + struct hirshfeld_atom_data { + double x; + double y; + double z; + double alpha; + double scale; + }; + + std::vector atom_data; + atom_data.reserve(natoms); + for( const auto& atom : mol ) { + const double radius = default_atomic_radius(atom.Z); + const double alpha = 0.5 / (radius * radius); + const double scale = atom.Z.get() * std::pow(alpha / pi, 1.5); + atom_data.push_back( hirshfeld_atom_data{ atom.x, atom.y, atom.z, alpha, scale } ); + } + + std::vector density(natoms); + + for( size_t i = 0; i < task.points.size(); ++i ) { + auto* weight_deri_ith = weight_deri + 3*natoms*i; + std::fill(weight_deri_ith, weight_deri_ith + 3*natoms, 0.); + + const auto& point = task.points[i]; + const double weight = task.weights[i]; + if( std::abs(weight) < 1.e-14 ) continue; + + const size_t iParent = task.iParent; + + double sum = 0.; + for( size_t iA = 0; iA < natoms; ++iA ) { + const auto& atom = atom_data[iA]; + const double da_x = point[0] - atom.x; + const double da_y = point[1] - atom.y; + const double da_z = point[2] - atom.z; + const double r2 = da_x*da_x + da_y*da_y + da_z*da_z; + density[iA] = atom.scale * std::exp(-atom.alpha * r2); + sum += density[iA]; + } + + if( sum <= 0. ) continue; + + auto* weight_deri_iParent = weight_deri_ith + 3*iParent; + for( size_t iB = 0; iB < natoms; ++iB ) { + if( iB == iParent ) continue; + + const auto& atom = atom_data[iB]; + const double scaled_density = density[iB] / sum; + const double coef = -2. * atom.alpha * scaled_density * weight; + + auto* weight_deri_iB = weight_deri_ith + 3*iB; + weight_deri_iB[0] = coef * (point[0] - atom.x); + weight_deri_iB[1] = coef * (point[1] - atom.y); + weight_deri_iB[2] = coef * (point[2] - atom.z); + + weight_deri_iParent[0] -= weight_deri_iB[0]; + weight_deri_iParent[1] -= weight_deri_iB[1]; + weight_deri_iParent[2] -= weight_deri_iB[2]; + } + } +} + /** @@ -983,6 +1128,97 @@ void reference_ssf_weights_1std_contraction_host( } +void reference_hirshfeld_weights_1std_contraction_host( + const Molecule& mol, + const MolMeta& meta, + const XCTask& task, + const double* w_times_f, + double* exc_grad_w +){ + + (void)meta; + + constexpr double pi = 3.141592653589793238462643383279502884; + constexpr double w_times_f_thresh = 1.e-14; + + const size_t natoms = mol.natoms(); + + struct hirshfeld_atom_data { + double x; + double y; + double z; + double alpha; + double scale; + }; + + std::vector atom_data; + atom_data.reserve(natoms); + for( const auto& atom : mol ) { + const double radius = default_atomic_radius(atom.Z); + const double alpha = 0.5 / (radius * radius); + const double scale = atom.Z.get() * std::pow(alpha / pi, 1.5); + atom_data.push_back( hirshfeld_atom_data{ atom.x, atom.y, atom.z, alpha, scale } ); + } + + std::vector density(natoms); + + for( size_t i = 0; i < task.points.size(); ++i ) { + const double w_times_f_i = w_times_f[i]; + if( std::abs(w_times_f_i) < w_times_f_thresh ) continue; + + const auto& point = task.points[i]; + const size_t iParent = task.iParent; + + double sum = 0.; + for( size_t iA = 0; iA < natoms; ++iA ) { + const auto& atom = atom_data[iA]; + const double da_x = point[0] - atom.x; + const double da_y = point[1] - atom.y; + const double da_z = point[2] - atom.z; + const double r2 = da_x*da_x + da_y*da_y + da_z*da_z; + density[iA] = atom.scale * std::exp(-atom.alpha * r2); + sum += density[iA]; + } + + if( sum <= 0. ) continue; + + double parent_x = 0.; + double parent_y = 0.; + double parent_z = 0.; + + for( size_t iB = 0; iB < natoms; ++iB ) { + if( iB == iParent ) continue; + + const auto& atom = atom_data[iB]; + const double scaled_density = density[iB] / sum; + const double coef = -2. * atom.alpha * scaled_density * w_times_f_i; + + const double grad_x = coef * (point[0] - atom.x); + const double grad_y = coef * (point[1] - atom.y); + const double grad_z = coef * (point[2] - atom.z); + + #pragma omp atomic + exc_grad_w[3*iB + 0] += grad_x; + #pragma omp atomic + exc_grad_w[3*iB + 1] += grad_y; + #pragma omp atomic + exc_grad_w[3*iB + 2] += grad_z; + + parent_x -= grad_x; + parent_y -= grad_y; + parent_z -= grad_z; + } + + #pragma omp atomic + exc_grad_w[3*iParent + 0] += parent_x; + #pragma omp atomic + exc_grad_w[3*iParent + 1] += parent_y; + #pragma omp atomic + exc_grad_w[3*iParent + 2] += parent_z; + } + +} + } diff --git a/src/xc_integrator/local_work_driver/host/reference/weights.hpp b/src/xc_integrator/local_work_driver/host/reference/weights.hpp index 7b79a1563..f0295373f 100644 --- a/src/xc_integrator/local_work_driver/host/reference/weights.hpp +++ b/src/xc_integrator/local_work_driver/host/reference/weights.hpp @@ -36,6 +36,13 @@ void reference_lko_weights_host( task_iterator task_end ); +void reference_hirshfeld_weights_host( + const Molecule& mol, + const MolMeta& meta, + task_iterator task_begin, + task_iterator task_end +); + void reference_becke_weights_1st_derivative_host( const Molecule& mol, const MolMeta& meta, @@ -50,6 +57,13 @@ void reference_ssf_weights_1st_derivative_host( double* weight_deri ); +void reference_hirshfeld_weights_1st_derivative_host( + const Molecule& mol, + const MolMeta& meta, + const XCTask& task, + double* weight_deri +); + // Becke weights 1st derivative contracted with integrator void reference_becke_weights_1std_contraction_host( const Molecule& mol, @@ -68,4 +82,12 @@ void reference_ssf_weights_1std_contraction_host( double* exc_grad_w ); +void reference_hirshfeld_weights_1std_contraction_host( + const Molecule& mol, + const MolMeta& meta, + const XCTask& task, + const double* w_times_f, + double* exc_grad_w +); + } diff --git a/src/xc_integrator/local_work_driver/host/reference_local_host_work_driver.cxx b/src/xc_integrator/local_work_driver/host/reference_local_host_work_driver.cxx index 192cfcd33..631507b79 100644 --- a/src/xc_integrator/local_work_driver/host/reference_local_host_work_driver.cxx +++ b/src/xc_integrator/local_work_driver/host/reference_local_host_work_driver.cxx @@ -51,6 +51,9 @@ namespace GauXC { case XCWeightAlg::LKO: reference_lko_weights_host( mol, meta, task_begin, task_end ); break; + case XCWeightAlg::Hirshfeld: + reference_hirshfeld_weights_host( mol, meta, task_begin, task_end ); + break; default: GAUXC_GENERIC_EXCEPTION("Weight Alg Not Supported"); } @@ -66,6 +69,9 @@ namespace GauXC { case XCWeightAlg::SSF: reference_ssf_weights_1std_contraction_host( mol, meta, task, w_times_f, exc_grad_w ); break; + case XCWeightAlg::Hirshfeld: + reference_hirshfeld_weights_1std_contraction_host( mol, meta, task, w_times_f, exc_grad_w ); + break; default: GAUXC_GENERIC_EXCEPTION("Weight Alg Not Supported"); } diff --git a/tests/c_api_test.cxx b/tests/c_api_test.cxx index e91941a42..5a3f0bf9b 100644 --- a/tests/c_api_test.cxx +++ b/tests/c_api_test.cxx @@ -138,6 +138,8 @@ TEST_CASE("C-API Enum Correspondence", "[c-api]") { == static_cast(GauXC_XCWeightAlg_SSF)); CHECK(static_cast(XCWeightAlg::LKO) == static_cast(GauXC_XCWeightAlg_LKO)); + CHECK(static_cast(XCWeightAlg::Hirshfeld) + == static_cast(GauXC_XCWeightAlg_Hirshfeld)); } // --- ExecutionSpace --- @@ -700,7 +702,25 @@ TEST_CASE("C-API MolecularWeights", "[c-api]") { C::gauxc_molecular_weights_factory_delete(&status, &mwf); } - SECTION("Modify weights with LoadBalancer") { + SECTION("Factory creation (Hirshfeld)") { + C::GauXCMolecularWeightsSettings settings; + settings.weight_alg = GauXC_XCWeightAlg_Hirshfeld; + settings.becke_size_adjustment = false; + + C::GauXCMolecularWeightsFactory mwf = C::gauxc_molecular_weights_factory_new( + &status, GauXC_ExecutionSpace_Host, "Default", settings); + CHECK(status.code == 0); + CHECK(mwf.ptr != nullptr); + + C::GauXCMolecularWeights mw = C::gauxc_molecular_weights_factory_get_instance(&status, mwf); + CHECK(status.code == 0); + CHECK(mw.ptr != nullptr); + + C::gauxc_molecular_weights_delete(&status, &mw); + C::gauxc_molecular_weights_factory_delete(&status, &mwf); + } + + auto test_modify_weights_with_load_balancer = [&](GauXC_XCWeightAlg weight_alg) { // Build full pipeline C::GauXCAtom atoms[3]; make_water_atoms(atoms); @@ -743,7 +763,7 @@ TEST_CASE("C-API MolecularWeights", "[c-api]") { REQUIRE(status.code == 0); C::GauXCMolecularWeightsSettings settings; - settings.weight_alg = GauXC_XCWeightAlg_SSF; + settings.weight_alg = weight_alg; settings.becke_size_adjustment = false; C::GauXCMolecularWeightsFactory mwf = C::gauxc_molecular_weights_factory_new( @@ -765,6 +785,14 @@ TEST_CASE("C-API MolecularWeights", "[c-api]") { C::gauxc_molgrid_delete(&status, &mg); C::gauxc_basisset_delete(&status, &basis); C::gauxc_molecule_delete(&status, &mol); + }; + + SECTION("Modify weights with LoadBalancer (SSF)") { + test_modify_weights_with_load_balancer(GauXC_XCWeightAlg_SSF); + } + + SECTION("Modify weights with LoadBalancer (Hirshfeld)") { + test_modify_weights_with_load_balancer(GauXC_XCWeightAlg_Hirshfeld); } } diff --git a/tests/ref_data/benzene_weights_hirshfeld.hdf5 b/tests/ref_data/benzene_weights_hirshfeld.hdf5 new file mode 100644 index 000000000..644919a22 Binary files /dev/null and b/tests/ref_data/benzene_weights_hirshfeld.hdf5 differ diff --git a/tests/weight_derivative_test.cxx b/tests/weight_derivative_test.cxx index ec53daf8b..ac94be40e 100644 --- a/tests/weight_derivative_test.cxx +++ b/tests/weight_derivative_test.cxx @@ -159,6 +159,9 @@ void test_weight_1st_deri_host_fdiff(const std::string& reference_file, XCWeight case XCWeightAlg::SSF: reference_ssf_weights_1st_derivative_host(mol, meta, task, analytical_derivatives.data()); break; + case XCWeightAlg::Hirshfeld: + reference_hirshfeld_weights_1st_derivative_host(mol, meta, task, analytical_derivatives.data()); + break; default: GAUXC_GENERIC_EXCEPTION("Weight Alg Not Supported"); } @@ -363,6 +366,9 @@ TEST_CASE("Weights First Derivative uncontracted HOST fidiff", "[weights_fdiff]" SECTION( "H3 SSF" ) { test_weight_1st_deri_host_fdiff(GAUXC_REF_DATA_PATH "/h3_blyp_cc-pvdz_ssf_gks.bin", XCWeightAlg::SSF, PruningScheme::Unpruned, 1.0e-5, 1.0e-6);} + SECTION( "H3 Hirshfeld" ) { + test_weight_1st_deri_host_fdiff(GAUXC_REF_DATA_PATH "/h3_blyp_cc-pvdz_ssf_gks.bin", XCWeightAlg::Hirshfeld, + PruningScheme::Unpruned, 1.0e-5, 1.0e-6);} } @@ -386,6 +392,9 @@ TEST_CASE("Weights First Derivative contracted HOST fidiff", "[weights_fdiff]") SECTION( "H3 SSF" ) { test_weight_1st_deri_host_fdiff_contracted(GAUXC_REF_DATA_PATH "/h3_blyp_cc-pvdz_ssf_gks.bin", XCWeightAlg::SSF, PruningScheme::Unpruned, 1.0e-5, 1.0e-6);} + SECTION( "H3 Hirshfeld" ) { + test_weight_1st_deri_host_fdiff_contracted(GAUXC_REF_DATA_PATH "/h3_blyp_cc-pvdz_ssf_gks.bin", XCWeightAlg::Hirshfeld, + PruningScheme::Unpruned, 1.0e-5, 1.0e-6);} // SECTION( "Benzene SSF" ) { // test_weight_1st_deri_host_fdiff_contracted(GAUXC_REF_DATA_PATH "/benzene_svwn5_cc-pvdz_ufg_ssf.hdf5", XCWeightAlg::SSF, // PruningScheme::Unpruned, 1.0e-5, 1.0e-6);} diff --git a/tests/weights.cxx b/tests/weights.cxx index 568822686..c779fb059 100644 --- a/tests/weights.cxx +++ b/tests/weights.cxx @@ -11,9 +11,11 @@ */ #include "ut_common.hpp" #include +#include #include #include #include +#include #include #include @@ -39,6 +41,7 @@ TEST_CASE( "Partition Weights", "[weights]" ) { generate_weights_data( mol, basis, "benzene_weights_becke.hdf5", XCWeightAlg::Becke ); generate_weights_data( mol, basis, "benzene_weights_ssf.hdf5", XCWeightAlg::SSF ); generate_weights_data( mol, basis, "benzene_weights_lko.hdf5", XCWeightAlg::LKO ); + generate_weights_data( mol, basis, "benzene_weights_hirshfeld.hdf5", XCWeightAlg::Hirshfeld ); return; #else @@ -52,6 +55,10 @@ TEST_CASE( "Partition Weights", "[weights]" ) { std::string ref_file = GAUXC_REF_DATA_PATH "/benzene_weights_lko.hdf5"; test_host_weights( ref_file, XCWeightAlg::LKO ); } + SECTION("Hirshfeld") { + std::string ref_file = GAUXC_REF_DATA_PATH "/benzene_weights_hirshfeld.hdf5"; + test_host_weights( ref_file, XCWeightAlg::Hirshfeld ); + } #endif @@ -79,4 +86,260 @@ TEST_CASE( "Partition Weights", "[weights]" ) { } } +#ifdef GAUXC_HAS_HOST +namespace { + +double hirshfeld_test_density( + const Atom& atom, + const std::array& point +) { + + constexpr double pi = 3.141592653589793238462643383279502884; + const double radius = default_atomic_radius(atom.Z); + const double alpha = 0.5 / (radius * radius); + + const double da_x = point[0] - atom.x; + const double da_y = point[1] - atom.y; + const double da_z = point[2] - atom.z; + const double r2 = da_x*da_x + da_y*da_y + da_z*da_z; + + return atom.Z.get() * std::pow(alpha / pi, 1.5) * std::exp(-alpha * r2); +} + +double hirshfeld_test_factor( + const Molecule& mol, + size_t parent, + const std::array& point +) { + + double sum = 0.; + for( const auto& atom : mol ) sum += hirshfeld_test_density(atom, point); + + return (sum > 0.) ? hirshfeld_test_density(mol[parent], point) / sum : 0.; +} + +Molecule translated_molecule( + const Molecule& mol, + const std::array& shift +) { + + Molecule translated; + translated.reserve(mol.size()); + for( const auto& atom : mol ) { + translated.emplace_back( + atom.Z, atom.x + shift[0], atom.y + shift[1], atom.z + shift[2] + ); + } + + return translated; +} + +std::vector> translated_points( + const std::vector>& points, + const std::array& shift +) { + + std::vector> translated = points; + for( auto& point : translated ) { + point[0] += shift[0]; + point[1] += shift[1]; + point[2] += shift[2]; + } + + return translated; +} + +} + +TEST_CASE( "Hirshfeld Partition Weights", "[weights]" ) { + + SECTION("H2 symmetric Gaussian pro-atom partition") { + Molecule mol; + mol.emplace_back( AtomicNumber(1), 0.0, 0.0, 0.0 ); + mol.emplace_back( AtomicNumber(1), 1.4, 0.0, 0.0 ); + MolMeta meta(mol); + + std::vector tasks(1); + tasks[0].iParent = 0; + tasks[0].points = {{ + {0.0, 0.0, 0.0}, + {0.7, 0.0, 0.0}, + {1.4, 0.0, 0.0} + }}; + tasks[0].weights = {1.0, 1.0, 1.0}; + + reference_hirshfeld_weights_host( mol, meta, tasks.begin(), tasks.end() ); + + CHECK( tasks[0].weights[0] == Approx(0.987761411136141) ); + CHECK( tasks[0].weights[1] == Approx(0.5) ); + CHECK( tasks[0].weights[2] == Approx(0.0122385888638589) ); + } + + SECTION("CO heteronuclear Gaussian pro-atom partition") { + Molecule mol; + mol.emplace_back( AtomicNumber(6), 0.0, 0.0, 0.0 ); + mol.emplace_back( AtomicNumber(8), 2.2, 0.0, 0.0 ); + MolMeta meta(mol); + + std::vector tasks(2); + tasks[0].iParent = 0; + tasks[0].points = {{ {1.1, 0.0, 0.0} }}; + tasks[0].weights = {1.0}; + tasks[1].iParent = 1; + tasks[1].points = tasks[0].points; + tasks[1].weights = {1.0}; + + reference_hirshfeld_weights_host( mol, meta, tasks.begin(), tasks.end() ); + + CHECK( tasks[0].weights[0] == Approx(0.348581523757377) ); + CHECK( tasks[1].weights[0] == Approx(0.651418476242623) ); + CHECK( tasks[0].weights[0] + tasks[1].weights[0] == Approx(1.0) ); + } + + SECTION("single atom preserves quadrature weights") { + Molecule mol; + mol.emplace_back( AtomicNumber(8), -0.3, 0.4, 1.1 ); + MolMeta meta(mol); + + std::vector tasks(1); + tasks[0].iParent = 0; + tasks[0].points = {{ + {-0.3, 0.4, 1.1}, + { 0.2, 0.9, 1.7}, + {-1.5, 0.2, 0.1} + }}; + tasks[0].weights = {0.25, 2.0, 7.5}; + const auto weights_ref = tasks[0].weights; + + reference_hirshfeld_weights_host( mol, meta, tasks.begin(), tasks.end() ); + + REQUIRE( tasks[0].weights.size() == weights_ref.size() ); + for( size_t i = 0; i < weights_ref.size(); ++i ) { + CHECK( tasks[0].weights[i] == Approx(weights_ref[i]) ); + } + } + + SECTION("partition factors conserve each original point weight") { + Molecule mol; + mol.emplace_back( AtomicNumber(6), -0.7, 0.0, 0.2 ); + mol.emplace_back( AtomicNumber(1), 1.2, -0.4, 0.5 ); + mol.emplace_back( AtomicNumber(8), 0.3, 1.1, -0.3 ); + mol.emplace_back( AtomicNumber(7), -1.4, -0.8, 0.9 ); + MolMeta meta(mol); + + const std::vector> points = {{ + {-0.1, 0.2, 0.0}, + { 0.9, -0.2, 0.4}, + {-1.0, -0.6, 1.3}, + { 1.8, 1.5, -0.9} + }}; + const std::vector point_weights = {0.25, 1.0, 3.5, 10.0}; + + std::vector tasks(mol.natoms()); + for( size_t parent = 0; parent < mol.natoms(); ++parent ) { + tasks[parent].iParent = static_cast(parent); + tasks[parent].points = points; + tasks[parent].weights = point_weights; + } + + reference_hirshfeld_weights_host( mol, meta, tasks.begin(), tasks.end() ); + + for( size_t ipt = 0; ipt < points.size(); ++ipt ) { + double partition_sum = 0.; + for( size_t parent = 0; parent < mol.natoms(); ++parent ) { + const double factor = hirshfeld_test_factor(mol, parent, points[ipt]); + CHECK( tasks[parent].weights[ipt] == Approx(point_weights[ipt] * factor) ); + partition_sum += tasks[parent].weights[ipt]; + } + CHECK( partition_sum == Approx(point_weights[ipt]) ); + } + } + + SECTION("rigid translation leaves partition factors unchanged") { + Molecule mol; + mol.emplace_back( AtomicNumber(6), -0.4, 0.0, 0.1 ); + mol.emplace_back( AtomicNumber(8), 1.6, 0.3, 0.0 ); + mol.emplace_back( AtomicNumber(1), 0.2, 1.5, 0.8 ); + + const std::vector> points = {{ + { 0.1, 0.1, 0.2}, + { 0.8, 0.4, 0.1}, + {-0.9, 1.2, 0.5} + }}; + const std::vector point_weights = {1.0, 2.0, 4.0}; + const std::array shift = {8.0, -3.5, 2.25}; + + auto shifted_mol = translated_molecule(mol, shift); + auto shifted_points = translated_points(points, shift); + + MolMeta meta(mol); + MolMeta shifted_meta(shifted_mol); + + std::vector tasks(1), shifted_tasks(1); + tasks[0].iParent = 1; + tasks[0].points = points; + tasks[0].weights = point_weights; + shifted_tasks[0].iParent = 1; + shifted_tasks[0].points = shifted_points; + shifted_tasks[0].weights = point_weights; + + reference_hirshfeld_weights_host( mol, meta, tasks.begin(), tasks.end() ); + reference_hirshfeld_weights_host( + shifted_mol, shifted_meta, shifted_tasks.begin(), shifted_tasks.end() + ); + + for( size_t ipt = 0; ipt < points.size(); ++ipt ) { + CHECK( shifted_tasks[0].weights[ipt] == Approx(tasks[0].weights[ipt]) ); + } + } + + SECTION("atom order permutation preserves physical atom partitions") { + Molecule mol; + mol.emplace_back( AtomicNumber(6), -0.8, 0.0, 0.0 ); + mol.emplace_back( AtomicNumber(8), 1.3, 0.2, -0.1 ); + mol.emplace_back( AtomicNumber(1), 0.0, -1.1, 0.7 ); + + Molecule permuted; + permuted.push_back(mol[2]); + permuted.push_back(mol[0]); + permuted.push_back(mol[1]); + const std::vector original_to_permuted = {1, 2, 0}; + + const std::vector> points = {{ + {-0.2, -0.1, 0.1}, + { 1.0, 0.4, -0.3}, + { 0.2, -0.9, 0.6} + }}; + + MolMeta meta(mol); + MolMeta permuted_meta(permuted); + + std::vector tasks(mol.natoms()), permuted_tasks(mol.natoms()); + for( size_t parent = 0; parent < mol.natoms(); ++parent ) { + tasks[parent].iParent = static_cast(parent); + tasks[parent].points = points; + tasks[parent].weights = {1.0, 1.0, 1.0}; + } + for( size_t parent = 0; parent < permuted.natoms(); ++parent ) { + permuted_tasks[parent].iParent = static_cast(parent); + permuted_tasks[parent].points = points; + permuted_tasks[parent].weights = {1.0, 1.0, 1.0}; + } + + reference_hirshfeld_weights_host( mol, meta, tasks.begin(), tasks.end() ); + reference_hirshfeld_weights_host( + permuted, permuted_meta, permuted_tasks.begin(), permuted_tasks.end() + ); + + for( size_t original_parent = 0; original_parent < mol.natoms(); ++original_parent ) { + const auto permuted_parent = original_to_permuted[original_parent]; + for( size_t ipt = 0; ipt < points.size(); ++ipt ) { + CHECK( permuted_tasks[permuted_parent].weights[ipt] == + Approx(tasks[original_parent].weights[ipt]) ); + } + } + } +} +#endif + diff --git a/tests/weights_generate.hpp b/tests/weights_generate.hpp index 34788701d..f95276a4e 100644 --- a/tests/weights_generate.hpp +++ b/tests/weights_generate.hpp @@ -66,15 +66,19 @@ void generate_weights_data( const Molecule& mol, const BasisSet& basis, switch( weight_alg ) { case XCWeightAlg::Becke: - reference_becke_weights_host( + reference_becke_weights_host( mol, lb.molmeta(), tasks.begin(), tasks.end() ); break; case XCWeightAlg::SSF: - reference_ssf_weights_host( + reference_ssf_weights_host( mol, lb.molmeta(), tasks.begin(), tasks.end() ); break; case XCWeightAlg::LKO: - reference_lko_weights_host( + reference_lko_weights_host( + mol, lb.molmeta(), tasks.begin(), tasks.end() ); + break; + case XCWeightAlg::Hirshfeld: + reference_hirshfeld_weights_host( mol, lb.molmeta(), tasks.begin(), tasks.end() ); break; } diff --git a/tests/weights_host.hpp b/tests/weights_host.hpp index 821330ae3..0045be2b1 100644 --- a/tests/weights_host.hpp +++ b/tests/weights_host.hpp @@ -26,17 +26,22 @@ void test_host_weights( const std::string& filename, XCWeightAlg weight_alg ) { switch(weight_alg) { case XCWeightAlg::Becke: - reference_becke_weights_host( + reference_becke_weights_host( ref_data.mol, *ref_data.meta, ref_data.tasks_unm.begin(), ref_data.tasks_unm.end() ); break; case XCWeightAlg::SSF: - reference_ssf_weights_host( + reference_ssf_weights_host( ref_data.mol, *ref_data.meta, ref_data.tasks_unm.begin(), ref_data.tasks_unm.end() ); break; case XCWeightAlg::LKO: - reference_lko_weights_host( + reference_lko_weights_host( + ref_data.mol, *ref_data.meta, ref_data.tasks_unm.begin(), + ref_data.tasks_unm.end() ); + break; + case XCWeightAlg::Hirshfeld: + reference_hirshfeld_weights_host( ref_data.mol, *ref_data.meta, ref_data.tasks_unm.begin(), ref_data.tasks_unm.end() ); break;