diff --git a/fastlib2/contrib/march/hf/build.py b/fastlib2/contrib/march/hf/build.py index 7f446d738d..d5f05a7d59 100644 --- a/fastlib2/contrib/march/hf/build.py +++ b/fastlib2/contrib/march/hf/build.py @@ -1,10 +1,3 @@ -librule( - name = "scf_solver", - headers = ["scf_solver.h"], - deplibs = ["fastlib:fastlib"], - tests = ["scf_solver_test.cc"] -) - librule( name = "dual_tree_integrals", headers = ["dual_tree_integrals.h", "square_tree.h"], @@ -14,10 +7,11 @@ librule( ) librule( - name = "hf", - headers = ["hf.h"], - deplibs = [":scf_solver", ":dual_tree_integrals", "fastlib:fastlib"], - tests = ["hf_test.cc"] + name = "scf_solver", + headers = ["scf_solver.h"], + sources = ["scf_solver.cc"], + deplibs = ["fastlib:fastlib", ":dual_tree_integrals"], + tests = ["scf_solver_test.cc"] ) binrule( @@ -32,4 +26,10 @@ binrule( headers = ["dual_tree_integrals.h", "naive_fock_matrix.h"], deplibs = ["fastlib:fastlib", ":dual_tree_integrals"], sources = ["fock_matrix_test.cc"] -) \ No newline at end of file +) + +binrule( + name = "hf", + deplibs = [":scf_solver", ":dual_tree_integrals", "fastlib:fastlib"], + sources = ["hf.cc"] +) diff --git a/fastlib2/contrib/march/hf/dual_tree_integrals.cc b/fastlib2/contrib/march/hf/dual_tree_integrals.cc index 50adec95e2..2ae7aac691 100644 --- a/fastlib2/contrib/march/hf/dual_tree_integrals.cc +++ b/fastlib2/contrib/march/hf/dual_tree_integrals.cc @@ -11,7 +11,7 @@ * Also, the integral project notes have a slightly different definition of * this function. I should make sure they're compatible. */ -double DualTreeIntegrals::ErfLikeFunction_(double z) { +double DualTreeIntegrals::ErfLikeFunction(double z) { if (z == 0) { return 1.0; @@ -39,7 +39,7 @@ double DualTreeIntegrals::ComputeSingleIntegral_(double mu_nu_dist, normalization_constant_fourth_; // the 0.25 comes from the four center distance identity - return_value = return_value * ErfLikeFunction_(bandwidth_ * four_way_dist); + return_value = return_value * ErfLikeFunction(bandwidth_ * four_way_dist); // I added a factor of 1/2, another mistake I think return_value = return_value * @@ -89,7 +89,7 @@ double DualTreeIntegrals::ComputeSingleIntegral_(const Vector& mu_center, // F(\alpha d^2(x_{i j} x_{k l})) // equivalent to F(\alpha return_value = return_value * - ErfLikeFunction_(bandwidth_ * four_centers_dists); + ErfLikeFunction(bandwidth_ * four_centers_dists); // exp(-\alpha (d^2(x_i, x_j) + d^2(x_k, x_l)) @@ -149,137 +149,3 @@ index_t DualTreeIntegrals::CountOnDiagonal_(SquareIntegralTree* rho_sigma) { } // CountOnDiagonal_() -/* -int main(int argc, char* argv[]) { - - - fx_init(argc, argv, NULL); - - DualTreeIntegrals integrals; - - const char* centers_file = fx_param_str(NULL, "centers", "test_centers.csv"); - Matrix centers_in; - data::Load(centers_file, ¢ers_in); - - index_t num_funs = centers_in.n_cols(); - - la::Scale((double)num_funs, ¢ers_in); - - data::Save("gaussian_test300.csv", centers_in); - - const char* density_file = fx_param_str(NULL, "density", "test_density.csv"); - Matrix density_in; - data::Load(density_file, &density_in); - - - Matrix core_in; - core_in.Init(num_funs, num_funs); - core_in.SetZero(); - - double bandwidth = fx_param_double(NULL, "bandwidth", 100.0); - - struct datanode* dual_mod = fx_submodule(NULL, "multi", "multi_tree"); - - integrals.Init(centers_in, dual_mod, density_in, core_in, bandwidth); - - fx_timer_start(dual_mod, "multi_tree"); - integrals.ComputeFockMatrix(); - fx_timer_stop(dual_mod, "multi_tree"); - - Matrix tree_coulomb; - Matrix tree_exchange; - ArrayList old_from_new; - - printf("MULTI-TREE:\n"); - integrals.OutputFockMatrix(&tree_coulomb, &tree_exchange, &old_from_new); - - Matrix naive_coulomb; - Matrix naive_exchange; - if (data::Load("naive_exchange.csv", &naive_exchange) == SUCCESS_FAIL) { - // printf("failed to load"); - naive_exchange.Destruct(); - NaiveFockMatrix naive; - - struct datanode* naive_mod = fx_submodule(NULL, "naive", "naive"); - naive.Init(centers_in, naive_mod, density_in, core_in, bandwidth); - - fx_timer_start(naive_mod, "naive"); - naive.ComputeFockMatrix(); - fx_timer_stop(naive_mod, "naive"); - - - - printf("\n\nNAIVE:\n"); - naive.PrintFockMatrix(&naive_coulomb, &naive_exchange); - - data::Save("naive_coulomb.csv", naive_coulomb); - data::Save("naive_exchange.csv", naive_exchange); - - printf("\n\n"); - } - else { - data::Load("naive_coulomb.csv", &naive_coulomb); - } - - - //printf("\n\n DIFFERENCE:\n"); - Matrix difference_mat; - - la::Scale(-1.0, &naive_coulomb); - difference_mat.Copy(tree_coulomb); - - Vector diff; - for (index_t i = 0; i < difference_mat.n_cols(); i++) { - - - - } // i - - double coulomb_error = 0.0; - double max_coulomb_error = 0.0; - double exchange_error = 0.0; - double max_exchange_error = 0.0; - - for (index_t i = 0; i < tree_exchange.n_cols(); i++) { - - for (index_t j = 0; j < tree_exchange.n_rows(); j++) { - - double this_val = fabs(naive_exchange.get(old_from_new[i], old_from_new[j]) - - tree_exchange.get(i, j)); - //printf("i:%d, j:%d, this_val:%g\n", i, j, this_val); - exchange_error = exchange_error + this_val; - - if (this_val > max_exchange_error) { - max_exchange_error = this_val; - //printf("max_error: (%d, %d)\n", i, j); - } - - this_val = fabs(naive_coulomb.get(old_from_new[i], old_from_new[j]) - - tree_coulomb.get(i, j)); - //printf("i:%d, j:%d, this_val:%g\n", i, j, this_val); - coulomb_error = coulomb_error + this_val; - - if (this_val > max_coulomb_error) { - max_coulomb_error = this_val; - //printf("max_error: (%d, %d)\n", i, j); - } - - } // j - - } // i - - exchange_error = exchange_error/(tree_exchange.n_cols() * tree_exchange.n_rows()); - coulomb_error = coulomb_error/(tree_exchange.n_cols() * tree_exchange.n_rows()); - - fx_format_result(NULL, "ave_abs_coulomb_error", "%g", coulomb_error); - fx_format_result(NULL, "max_abs_coulomb_error", "%g", max_coulomb_error); - - fx_format_result(NULL, "ave_abs_exchange_error", "%g", exchange_error); - fx_format_result(NULL, "max_abs_exchange_error", "%g", max_exchange_error); - - fx_done(NULL); - - return 0; - -} -*/ \ No newline at end of file diff --git a/fastlib2/contrib/march/hf/dual_tree_integrals.h b/fastlib2/contrib/march/hf/dual_tree_integrals.h index 7499050487..5e1164e9fb 100644 --- a/fastlib2/contrib/march/hf/dual_tree_integrals.h +++ b/fastlib2/contrib/march/hf/dual_tree_integrals.h @@ -148,7 +148,6 @@ class DualTreeIntegrals { // the square tree SquareIntegralTree* square_tree_; - // The centers of the identical width, spherical Gaussian basis functions Matrix centers_; @@ -180,6 +179,8 @@ class DualTreeIntegrals { double epsilon_coulomb_absolute_; double epsilon_exchange_absolute_; + // Values that can be guaranteed to be less than this value are pruned by the + // absolute criterion double hybrid_cutoff_; // The return values are stored here @@ -209,8 +210,6 @@ class DualTreeIntegrals { // The normalization constant to the fourth power double normalization_constant_fourth_; - index_t num_mu_nu_base_cases_; - index_t num_absolute_prunes_; index_t num_relative_prunes_; @@ -246,9 +245,6 @@ class DualTreeIntegrals { } // PreOrderTraversal - /** - * Returns the global upper bound for the given - */ double ComputeUpperBound_() { return 0.0; @@ -588,18 +584,16 @@ class DualTreeIntegrals { - - double ErfLikeFunction_(double z); - - double ComputeSingleIntegral_(double mu_nu_dist, double rho_sigma_dist, double four_way_dist); - + + double ComputeSingleIntegral_(const Vector& mu_center, const Vector& nu_center, const Vector& rho_center, const Vector& sigma_center); - + + bool RectangleOnDiagonal_(IntegralTree* mu, IntegralTree* nu); @@ -654,9 +648,12 @@ class DualTreeIntegrals { // Multiply by normalization to the fourth, since it appears // once in each of the four integrals - double this_integral = density_matrix_.ref(rho_index, sigma_index) * + + double this_integral = density_matrix_.get(rho_index, sigma_index) * ComputeSingleIntegral_(mu_vec, nu_vec, rho_vec, sigma_vec); + // this line gets it right + // double this_integral = ComputeSingleIntegral_(mu_vec, nu_vec, rho_vec, sigma_vec); if (rho != sigma) { this_integral = this_integral * 2; } @@ -1246,34 +1243,92 @@ class DualTreeIntegrals { } // SetEntryBounds_ + void ResetTreeForExchange_(SquareIntegralTree* root) { + + if (root != NULL) { + + root->stat().set_remaining_references(number_of_basis_functions_ * + number_of_basis_functions_); + + root->stat().set_approximation_val(0.0); + + ResetTreeForExchange_(root->left()); + ResetTreeForExchange_(root->right()); + + } + + } // ResetTreeForExchange_() + /** - * Resets the tree after the Coulomb computation and before the exchange + * Resets the tree after the density matrix changes. + * + * This is also used between Coulomb and Exchange computations, which is + * probably wrong. */ void ResetTree_(SquareIntegralTree* root) { - if (root != NULL) { - root->stat().set_remaining_references(number_of_basis_functions_ * - number_of_basis_functions_); - - root->stat().set_approximation_val(0.0); + double max_density; + double min_density; + + if (root->is_leaf()) { + + max_density = -DBL_INF; + min_density = DBL_INF; + + for (index_t i = root->query1()->begin(); i < root->query1()->end(); + i++) { + + for (index_t j = root->query2()->begin(); j < root->query2()->end(); + j++) { + + double this_density = density_matrix_.ref(i, j); + if (this_density > max_density) { + max_density = this_density; + } + if (this_density < min_density) { + min_density = this_density; + } + + } // j + + } // i + + } // leaf + else { ResetTree_(root->left()); ResetTree_(root->right()); + + max_density = max(root->left()->stat().density_upper_bound(), + root->right()->stat().density_upper_bound()); + min_density = min(root->left()->stat().density_lower_bound(), + root->right()->stat().density_lower_bound()); + + } // non-leaf - } + root->stat().set_density_upper_bound(max_density); + root->stat().set_density_lower_bound(min_density); + + root->stat().set_remaining_references(number_of_basis_functions_ * + number_of_basis_functions_); + + root->stat().set_approximation_val(0.0); } // ResetTree_() public: + double ErfLikeFunction(double z); + + /** * Initialize the class with the centers of the data points, the fx module, * bandwidth */ - void Init(const Matrix& centers_in, struct datanode* mod, - const Matrix& density_in, const Matrix& core_in, double band) { + void Init(const Matrix& centers_in, struct datanode* mod, double band) { + // Needs to be copied because it will be permuted centers_.Copy(centers_in); module_ = mod; @@ -1306,12 +1361,6 @@ public: // The common normalization constant of all the Gaussians normalization_constant_fourth_ = pow((2 * bandwidth_ / math::PI), 3); - // Not sure this is right, might want to consider starting with the previous - // iteration's version - fock_matrix_.Copy(core_in); - - density_matrix_.Copy(density_in); - coulomb_matrix_.Init(number_of_basis_functions_, number_of_basis_functions_); coulomb_matrix_.SetZero(); @@ -1319,6 +1368,9 @@ public: exchange_matrix_.Init(number_of_basis_functions_, number_of_basis_functions_); exchange_matrix_.SetZero(); + + fock_matrix_.Init(number_of_basis_functions_, number_of_basis_functions_); + fock_matrix_.SetZero(); leaf_size_ = fx_param_int(module_, "leaf_size", 20); @@ -1327,19 +1379,58 @@ public: // Set up the indices of the nodes for symmetry pruning traversal_index_ = 0; - PreOrderTraversal_(tree_); + //PreOrderTraversal_(tree_); //tree_->Print(); square_tree_ = new SquareIntegralTree(); - square_tree_->Init(tree_, tree_, density_matrix_); + square_tree_->Init(tree_, tree_, number_of_basis_functions_); //square_tree_->Print(); - SetEntryBounds_(); - - - num_mu_nu_base_cases_ = 0; + //SetEntryBounds_(); } // Init + /** + * Returns the permutation of the basis centers + */ + void GetPermutation(ArrayList* perm) { + + perm->Copy(old_from_new_centers_); + + } // GetPermutation() + + /** + * Call this after the density matrix is permuted in the SCF solver + */ + void GetDensity(const Matrix& updated_density) { + + density_matrix_.Copy(updated_density); + + ResetTree_(square_tree_); + + + SetEntryBounds_(); + + } + + /** + * Updates the density matrix and clears the Fock matrix between iterations + * Should also reset the tree for the next iteration. + */ + void UpdateMatrices(const Matrix& new_density) { + + //This isn't necessary since it's already an alias + density_matrix_.CopyValues(new_density); + + // Reset tree density bounds + ResetTree_(square_tree_); + + SetEntryBounds_(); + + coulomb_matrix_.SetZero(); + exchange_matrix_.SetZero(); + + } // UpdateMatrices() + /** * Drives the computation, assuming that all the parameters are correct */ @@ -1353,7 +1444,7 @@ public: // matrix // I think this is the only resetting the tree will need SetEntryBounds_(); - ResetTree_(square_tree_); + ResetTreeForExchange_(square_tree_); ComputeExchangeRecursion_(square_tree_, square_tree_); // Then by adding both into the Fock matrix la::AddTo(coulomb_matrix_, &fock_matrix_); @@ -1362,11 +1453,16 @@ public: } // ComputeTwoElectronIntegrals + const Matrix& FockMatrix() const { + return fock_matrix_; + } + /** * Returns the computed Fock matrix. For now, it just prints it, but should * eventually return it in a useable form for the SCF procedure. */ - void OutputFockMatrix(Matrix* coulomb_out, Matrix* exchange_out, + void OutputFockMatrix(Matrix* fock_out, Matrix* coulomb_out, + Matrix* exchange_out, ArrayList* old_from_new) { //printf("number_of_approximations_ = %d\n", number_of_approximations_); @@ -1386,20 +1482,33 @@ public: fx_format_result(module_, "abs_prunes", "%d", num_absolute_prunes_); fx_format_result(module_, "rel_prunes", "%d", num_relative_prunes_); - //fx_format_result(module_, "num_0_3_bases", "%d", num_mu_nu_base_cases_); - /* printf("Multi-tree Coulomb:\n"); coulomb_matrix_.PrintDebug(); printf("Multi-tree Exchange:\n"); exchange_matrix_.PrintDebug(); */ + if (fock_out) { + fock_out->Copy(fock_matrix_); + } + if (coulomb_out) { + coulomb_out->Copy(coulomb_matrix_); + } + if (exchange_out) { + exchange_out->Copy(exchange_matrix_); + } + + if (old_from_new) { + old_from_new->Copy(old_from_new_centers_); + } + + // Need to output the Fock matrix + // should I unpermute here? + // Maybe keep it permuted in the other code, and unpermute at the end? + + // For now, unpermute it here - coulomb_out->Copy(coulomb_matrix_); - exchange_out->Copy(exchange_matrix_); - old_from_new->Copy(old_from_new_centers_); - } // OutputFockMatrix() diff --git a/fastlib2/contrib/march/hf/fock_matrix_test.cc b/fastlib2/contrib/march/hf/fock_matrix_test.cc new file mode 100644 index 0000000000..9b5b30f3df --- /dev/null +++ b/fastlib2/contrib/march/hf/fock_matrix_test.cc @@ -0,0 +1,348 @@ +#include "dual_tree_integrals.h" +#include "naive_fock_matrix.h" +//#include + +class FockMatrixTest { + +public: + + void Setup(const char* centers_name, const char* density_name, + struct datanode* multi_mod, struct datanode* naive_mod, + double band) { + + Matrix test_centers; + data::Load(centers_name, &test_centers); + + Matrix test_density; + data::Load(density_name, &test_density); + + Matrix test_core; + test_core.Init(test_centers.n_cols(), test_centers.n_cols()); + test_core.SetZero(); + + multi_ = new DualTreeIntegrals(); + multi_->Init(test_centers, multi_mod, band); + multi_->GetDensity(test_density); + + naive_ = new NaiveFockMatrix(); + naive_->Init(test_centers, naive_mod, test_density, test_core, band); + + } // Setup + + + void Destruct() { + + old_from_new.Destruct(); + delete multi_; + delete naive_; + + } // Destruct + + void CompareMatrices(struct datanode* mod) { + + double coulomb_error = 0.0; + double max_coulomb_error = 0.0; + double exchange_error = 0.0; + double max_exchange_error = 0.0; + + double rel_coulomb_error = 0.0; + double max_rel_coulomb_error = 0.0; + double rel_exchange_error = 0.0; + double max_rel_exchange_error = 0.0; + + double min_coulomb = DBL_INF; + double max_coulomb = -DBL_INF; + double min_exchange = DBL_INF; + double max_exchange = -DBL_INF; + + Matrix multi_exchange; + multi_exchange.Alias(multi_->exchange_matrix_); + Matrix multi_coulomb; + multi_coulomb.Alias(multi_->coulomb_matrix_); + Matrix naive_exchange; + naive_exchange.Alias(naive_->exchange_matrix_); + Matrix naive_coulomb; + naive_coulomb.Alias(naive_->coulomb_matrix_); + + /*printf("MULTI:\n"); + multi_exchange.PrintDebug(); + printf("NAIVE:\n"); + naive_exchange.PrintDebug(); + printf("\n\n"); + */ + + index_t num_rows = multi_exchange.n_rows(); + //printf("num_rows: %d\n", num_rows); + + for (index_t i = 0; i < num_rows; i++) { + + for (index_t j = 0; j < num_rows; j++) { + + double this_val; + double this_naive; + this_naive = fabs(naive_exchange.get(old_from_new[i], old_from_new[j])); + //printf("i:%d, j:%d, this_val:%g\n", i, j, this_val); + if (this_naive > max_exchange) { + max_exchange = this_naive; + } + if (this_naive < min_exchange) { + min_exchange = this_naive; + } + + this_val = fabs(multi_exchange.get(i, j) - + naive_exchange.get(old_from_new[i], old_from_new[j])); + exchange_error = exchange_error + this_val; + + if (this_val > max_exchange_error) { + max_exchange_error = this_val; + //printf("max_error: (%d, %d)\n", i, j); + } + + this_val = this_val/this_naive; + rel_exchange_error = rel_exchange_error + this_val; + + if (this_val > max_rel_exchange_error) { + max_rel_exchange_error = this_val; + } + + //////////// Coulomb ////////////// + + this_naive = fabs(naive_coulomb.get(old_from_new[i], old_from_new[j])); + this_val = fabs(naive_coulomb.get(old_from_new[i], old_from_new[j]) + - multi_coulomb.get(i, j)); + + if (this_naive > max_coulomb) { + max_coulomb = this_naive; + } + if (this_naive < min_coulomb) { + min_coulomb = this_naive; + } + + coulomb_error = coulomb_error + this_val; + + if (this_val > max_coulomb_error) { + max_coulomb_error = this_val; + } + + this_val = this_val/this_naive; + rel_coulomb_error = rel_coulomb_error + this_val; + + if (this_val > max_rel_coulomb_error) { + max_rel_coulomb_error = this_val; + } + + } // j + + } // i + + exchange_error = exchange_error/(num_rows * num_rows); + coulomb_error = coulomb_error/(num_rows * num_rows); + + rel_exchange_error = rel_exchange_error/(num_rows * num_rows); + rel_coulomb_error = rel_coulomb_error/(num_rows * num_rows); + + fx_format_result(mod, "max_coulomb", "%g", max_coulomb); + fx_format_result(mod, "min_coulomb", "%g", min_coulomb); + fx_format_result(mod, "max_exchange", "%g", max_exchange); + fx_format_result(mod, "min_exchange", "%g", min_exchange); + + + fx_format_result(mod, "ave_abs_coulomb_error", "%g", coulomb_error); + fx_format_result(mod, "max_abs_coulomb_error", "%g", max_coulomb_error); + + fx_format_result(mod, "ave_abs_exchange_error", "%g", exchange_error); + fx_format_result(mod, "max_abs_exchange_error", "%g", max_exchange_error); + + fx_format_result(mod, "ave_rel_coulomb_error", "%g", rel_coulomb_error); + fx_format_result(mod, "max_rel_coulomb_error", "%g", max_rel_coulomb_error); + + fx_format_result(mod, "ave_rel_exchange_error", "%g", rel_exchange_error); + fx_format_result(mod, "max_rel_exchange_error", "%g", + max_rel_exchange_error); + + } // CompareMatrices() + + void TestMatricesSmall() { + + + struct datanode* multi_mod_small = fx_submodule(NULL, "multi_small", + "multi_small"); + + fx_set_param(multi_mod_small, "epsilon", "0.0"); + + struct datanode* naive_mod_small = fx_submodule(NULL, "naive_small", + "naive_small"); + + fx_set_param(naive_mod_small, "coulomb_output", "naive_coulomb_4_1.0.csv"); + fx_set_param(naive_mod_small, "exchange_output", + "naive_exchange_4_1.0.csv"); + + Setup("test_centers.csv", "test_density.csv", + multi_mod_small, naive_mod_small, 1.0); + + multi_->ComputeFockMatrix(); + Matrix cou; + Matrix exc; + multi_->OutputFockMatrix(NULL, &cou, &exc, &old_from_new); + + naive_->ComputeFockMatrix(); + + CompareMatrices(multi_mod_small); + + + Destruct(); + + } // TestMatricesSmall() + + + void TestMatricesMidNoPrune() { + + struct datanode* multi_mod_noprune = fx_submodule(NULL, "multi_noprune", + "multi_noprune"); + + fx_set_param(multi_mod_noprune, "epsilon", "0.0"); + fx_set_param(multi_mod_noprune, "hybrid_cutoff", "0.0"); + fx_set_param(multi_mod_noprune, "epsilon_absolute", "0.0"); + fx_set_param(multi_mod_noprune, "leaf_size", "1"); + + struct datanode* naive_mod_noprune = fx_submodule(NULL, "naive_noprune", + "naive_noprune"); + + fx_set_param(naive_mod_noprune, "coulomb_output", "naive_coulomb_10_1.0.csv"); + fx_set_param(naive_mod_noprune, "exchange_output", + "naive_exchange_10_1.0.csv"); + + Setup("test_centers_10.csv", "test_density_10.csv", + multi_mod_noprune, naive_mod_noprune, 1.0); + + multi_->ComputeFockMatrix(); + Matrix cou; + Matrix exc; + multi_->OutputFockMatrix(NULL, &cou, &exc, &old_from_new); + + + naive_->ComputeFockMatrix(); + + CompareMatrices(multi_mod_noprune); + + + Destruct(); + + + + } // TestMatricesMidNoPrune + + void TestMatricesMidPrune() { + + struct datanode* multi_mod_prune = fx_submodule(NULL, "multi_prune", + "multi_prune"); + + fx_set_param(multi_mod_prune, "epsilon", "1.0"); + fx_set_param(multi_mod_prune, "leaf_size", "1"); + + + struct datanode* naive_mod_prune = fx_submodule(NULL, "naive_prune", + "naive_prune"); + + fx_set_param(naive_mod_prune, "coulomb_output", "naive_coulomb_10_1.0.csv"); + fx_set_param(naive_mod_prune, "exchange_output", + "naive_exchange_10_1.0.csv"); + + Setup("test_centers_10.csv", "test_density_10.csv", + multi_mod_prune, naive_mod_prune, 1.0); + + multi_->ComputeFockMatrix(); + Matrix cou; + Matrix exc; + multi_->OutputFockMatrix(NULL, &cou, &exc, &old_from_new); + + naive_->ComputeFockMatrix(); + + CompareMatrices(multi_mod_prune); + + + Destruct(); + + } // TestMatricesMidPrune + + void TestMatricesLarge() { + + struct datanode* multi_mod_large = fx_submodule(NULL, "multi_large", + "multi_large"); + + fx_set_param(multi_mod_large, "epsilon", "0.05"); + fx_set_param(multi_mod_large, "leaf_size", "5"); + fx_set_param(multi_mod_large, "hybrid_cutoff", "50"); + fx_set_param(multi_mod_large, "epsilon_absolute", "25"); + + + struct datanode* naive_mod_large = fx_submodule(NULL, "naive_large", + "naive_large"); + + fx_set_param(naive_mod_large, "coulomb_output", + "naive_coulomb_100_0.01.csv"); + fx_set_param(naive_mod_large, "exchange_output", + "naive_exchange_100_0.01.csv"); + + Setup("test_centers_100.csv", "test_density_100.csv", + multi_mod_large, naive_mod_large, 0.01); + + fx_timer_start(multi_mod_large, "multi"); + multi_->ComputeFockMatrix(); + fx_timer_stop(multi_mod_large, "multi"); + Matrix cou; + Matrix exc; + multi_->OutputFockMatrix(NULL, &cou, &exc, &old_from_new); + + fx_timer_start(naive_mod_large, "naive"); + naive_->ComputeFockMatrix(); + fx_timer_stop(naive_mod_large, "naive"); + + CompareMatrices(multi_mod_large); + + + Destruct(); + + } // TestMatricesLarge + + void TestMatrices() { + + NONFATAL("SMALL Test\n"); + TestMatricesSmall(); + + NONFATAL("NOPRUNE Test\n"); + TestMatricesMidNoPrune(); + + NONFATAL("PRUNE Test\n"); + TestMatricesMidPrune(); + + + if (fx_param_exists(NULL, "large")) { + NONFATAL("LARGE Test\n"); + TestMatricesLarge(); + } + } // TestMatrices + +private: + + DualTreeIntegrals* multi_; + + NaiveFockMatrix* naive_; + + ArrayList old_from_new; + +}; //class FockMatrixTest + + +int main(int argc, char* argv[]) { + + fx_init(argc, argv); + + FockMatrixTest tester; + tester.TestMatrices(); + + fx_done(); + + return 0; + +} // main \ No newline at end of file diff --git a/fastlib2/contrib/march/hf/hf.cc b/fastlib2/contrib/march/hf/hf.cc index 6e6ed96db1..99281a2701 100644 --- a/fastlib2/contrib/march/hf/hf.cc +++ b/fastlib2/contrib/march/hf/hf.cc @@ -10,40 +10,55 @@ int main(int argc, char *argv[]) { - fx_init(argc, argv, NULL); + fx_init(argc, argv); - ////////////// Read in data ////////////// + struct datanode* mod = fx_submodule(NULL, "hf", "hf"); - // How will the data be organized? - // What is the best format to read in basis functions? - // I will likely need my own function to parse basis functions - // Check out the PSI3 code + int num_electrons = fx_param_int_req(NULL, "num_electrons"); + + const char* centers_file = fx_param_str(NULL, "basis_centers", + "test_centers.csv"); + Matrix centers; + data::Load(centers_file, ¢ers); + + const char* nuclear_file = fx_param_str(NULL, "nuclear_centers", + "test_nuclear_centers.csv"); + + Matrix nuclear; + data::Load(nuclear_file, &nuclear); + + const char* nuclear_mass_file = fx_param_str(NULL, "nuclear_masses", + "test_nuclear_masses.csv"); + + Matrix nuclear_masses; + data::Load(nuclear_mass_file, &nuclear_masses); + + // Need to double check if this is right + if (nuclear.n_cols() != nuclear_masses.n_rows()) { + FATAL("Number of masses must equal number of nuclear coordinates!\n"); + } + + Vector nuclear_mass; + nuclear_masses.MakeColumnVector(0, &nuclear_mass); + + Matrix density; + if (fx_param_exists(NULL, "initial_density")) { + const char* density_file = fx_param_str_req(NULL, "initial_density"); + data::Load(density_file, &density); + } + else { + density.Init(centers.n_cols(), centers.n_cols()); + density.SetZero(); + } + SCFSolver solver; - ////////////// Compute the integrals /////////// + solver.Init(mod, num_electrons, centers, density, nuclear, nuclear_mass); - // Should this be in the same, or a different class from the linear system - // solver? + solver.ComputeWavefunction(); - Matrix fock_matrix; - Matrix overlap_matrix; - - ////////////// Solve the linear system ///////////// - - //HFSolver solver; - //solver.Init(fock_matrix, overlap_matrix); - - - - - //////////// Output the results /////////////////// - - // Total energy - // Spin orbitals: both filled and virtual - - - fx_done(NULL); + fx_done(); return 0; diff --git a/fastlib2/contrib/march/hf/hf.h b/fastlib2/contrib/march/hf/hf.h index d18ad3f867..70d445975c 100644 --- a/fastlib2/contrib/march/hf/hf.h +++ b/fastlib2/contrib/march/hf/hf.h @@ -3,48 +3,15 @@ * * @author Bill March (march@gatech.edu) * - * Contains classes for the Hartree-Fock implementation. + * Header file for the HF implementation as a whole. */ #ifndef HF_H #define HF_H +#include "dual_tree_integrals.h" +#include "scf_solver.h" #include -/** - * A class that stores the information for a contracted Gaussian basis function. - * - * TODO: Should this class also have functions for computations among contracted - * functions? - */ -class ContractedGaussian { - - FORBID_ACCIDENTAL_COPIES(ContractedGaussian); - - private: - - // The number of primitive Gaussians that make up this function - index_t number_of_primitives_; - ArrayList bandwidths_; - ArrayList coefficients_; - - - public: - - ContractedGaussian() {} - - ~ContractedGaussian() {} - - void Init(index_t num, const ArrayList& band, - const ArrayList& coeff) { - - number_of_primitives_ = num; - bandwidths_.Copy(band); - coefficients_.Copy(coeff); - - } - - -}; diff --git a/fastlib2/contrib/march/hf/scf_solver.h b/fastlib2/contrib/march/hf/scf_solver.h index 750310496c..56e3e31a0f 100644 --- a/fastlib2/contrib/march/hf/scf_solver.h +++ b/fastlib2/contrib/march/hf/scf_solver.h @@ -10,14 +10,12 @@ #define SCF_SOLVER_H #include +#include "dual_tree_integrals.h" /** * Algorithm class for the SCF part of the HF computation. This class assumes * the integrals have been computed and does the SVD-like part of the * computation. - * - * For now, this is simply an implementation of the basic algorithm. In the - * future, I should examine how this could be done better. */ class SCFSolver { @@ -26,20 +24,21 @@ class SCFSolver { FORBID_ACCIDENTAL_COPIES(SCFSolver); private: - - // I can probably be more efficient in terms of storing these matrices - // I don't want to store many matrices of this size in the final code - - Vector two_electron_integrals_; - Matrix one_electron_integrals_; // T + V + Matrix basis_centers_; + Matrix nuclear_centers_; + + Vector nuclear_masses_; + + Matrix core_matrix_; // T + V Matrix kinetic_energy_integrals_; // T Matrix potential_energy_integrals_; // V Matrix coefficient_matrix_; // C or C' // Consider changing this name to reflect that it's the change of basis matrix - Matrix overlap_matrix_; // S or S^{-1/2} + Matrix overlap_matrix_; // S + Matrix change_of_basis_matrix_; // S^{-1/2} Matrix density_matrix_; // D Matrix fock_matrix_; // F or F', depending on the basis @@ -47,7 +46,10 @@ class SCFSolver { index_t number_of_basis_functions_; // N index_t number_of_electrons_; // K + index_t number_of_nuclei_; + index_t number_to_fill_; + // I think I'll have to compute this in the beginning double nuclear_repulsion_energy_; ArrayList total_energy_; @@ -62,105 +64,134 @@ class SCFSolver { struct datanode* module_; + DualTreeIntegrals integrals_; + + ArrayList occupied_indices_; + + ArrayList old_from_new_centers_; + + double bandwidth_; + public: SCFSolver() {} ~SCFSolver() {} - /** - * Initialize the class with const references to the electron matrices and - * the overlap matrix, both of which should have been computed already. - */ - void Init(double nuclear_energy, const Matrix& overlap_in, - const Matrix& kinetic_in, const Matrix& potential_in, - const Vector& two_electron_in, index_t num_electrons, - double converged, struct datanode* mod) { + + void Init(struct datanode* mod, index_t num_electrons, + const Matrix& basis_centers, const Matrix& density, + const Matrix& nuclear, const Vector& nuclear_mass) { - nuclear_repulsion_energy_ = nuclear_energy; + module_ = mod; number_of_electrons_ = num_electrons; - // Read in integrals - overlap_matrix_.Copy(overlap_in); - kinetic_energy_integrals_.Copy(kinetic_in); - potential_energy_integrals_.Copy(potential_in); + struct datanode* integral_mod = fx_submodule(module_, "integrals", + "integrals"); - // Is this really the best thing to do? - // Copying this whole thing might be too expensive - two_electron_integrals_.Copy(two_electron_in); + bandwidth_ = fx_param_double(module_, "bandwidth", 0.1); - number_of_basis_functions_ = overlap_matrix_.n_cols(); + integrals_.Init(basis_centers, integral_mod, bandwidth_); - DEBUG_ASSERT(number_of_basis_functions_ >= number_of_electrons_); + // Need to get out the permutation from the integrals_, then use it to + // permute the basis centers - // Form the core Hamiltonian - la::AddInit(kinetic_energy_integrals_, potential_energy_integrals_, - &one_electron_integrals_); + integrals_.GetPermutation(&old_from_new_centers_); + + PermuteMatrix_(basis_centers, &basis_centers_, old_from_new_centers_); + PermuteMatrix_(density, &density_matrix_, old_from_new_centers_); + + integrals_.GetDensity(density_matrix_); + + nuclear_centers_.Copy(nuclear); + + nuclear_masses_.Copy(nuclear_mass); + + number_of_nuclei_ = nuclear_centers_.n_cols(); + + number_to_fill_ = (index_t)ceil((double)number_of_electrons_/2); + occupied_indices_.Init(number_to_fill_); + + DEBUG_ASSERT(number_of_nuclei_ == nuclear_masses_.length()); + + + number_of_basis_functions_ = basis_centers_.n_cols(); + + DEBUG_ASSERT(number_of_basis_functions_ >= number_to_fill_); // Empty inits to prevent errors on closing + overlap_matrix_.Init(number_of_basis_functions_, + number_of_basis_functions_); + kinetic_energy_integrals_.Init(number_of_basis_functions_, + number_of_basis_functions_); + potential_energy_integrals_.Init(number_of_basis_functions_, + number_of_basis_functions_); + coefficient_matrix_.Init(number_of_basis_functions_, number_of_basis_functions_); - density_matrix_.Init(number_of_basis_functions_, - number_of_basis_functions_); - fock_matrix_.Init(number_of_basis_functions_, number_of_basis_functions_); energy_vector_.Init(number_of_basis_functions_); total_energy_.Init(expected_number_of_iterations_); - convergence_tolerance_ = converged; + convergence_tolerance_ = fx_param_double(module_, "convergence_tolerance", + 0.1); + // Need to double check that this is right density_matrix_frobenius_norm_ = DBL_MAX; - module_ = mod; + current_iteration_ = 0; - } // Init + } // Init() private: + + double ComputeOverlapIntegral_(double dist); + double ComputeKineticIntegral_(double dist); + + double ComputeNuclearIntegral_(const Vector& nuclear_position, + const Vector& mu, const Vector& nu); + /** - * Finds a pairwise index for use in finding the entire index of an integral. - */ - index_t FindIntegralIndexHelper_(index_t mu, index_t nu) { + * Permutes the matrix mat according to the permutation given. The permuted + * matrix is written to new_mat, overwriting whatever was there before + */ + void PermuteMatrix_(const Matrix& old_mat, Matrix* new_mat, + const ArrayList& perm) { + + index_t num_cols = old_mat.n_cols(); + DEBUG_ASSERT(num_cols == perm.size()); + + new_mat->Init(old_mat.n_rows(), num_cols); + + for (index_t i = 0; i < num_cols; i++) { - //DEBUG_ASSERT(mu >= nu); - if (mu >= nu) { - return (((mu * (mu + 1))/2) + nu); - } - else { - return (((nu * (nu + 1))/2) + mu); - } - } // FindIntegralIndexHelper_ - - /** - * Find the index of the two electron integral (mu nu | rho sigma) in the - * two-electron integrals array - */ - index_t FindIntegralIndex_(index_t mu, index_t nu, index_t rho, index_t sig) { - - //DEBUG_ASSERT(mu >= nu); - //DEBUG_ASSERT(rho >= sig); - - index_t mu_nu = FindIntegralIndexHelper_(mu, nu); - index_t rho_sig = FindIntegralIndexHelper_(rho, sig); - - //DEBUG_ASSERT(mu_nu >= rho_sig); - if (mu_nu >= rho_sig) { - return (FindIntegralIndexHelper_(mu_nu, rho_sig)); - } - else { - return (FindIntegralIndexHelper_(rho_sig, mu_nu)); - } + Vector old_vec; + old_mat.MakeColumnVector(i, &old_vec); + Vector new_vec; + new_mat->MakeColumnVector(perm[i], &new_vec); - }// FindIntegralIndex_ + new_vec.CopyValues(old_vec); + + } + } // PermuteMatrix_() - /** - * Create the matrix S^{-1/2} using the eigenvector decomposition. Overwrites - * overlap_matrix_ with S^{-1/2}. + /** + * Given the basis set and nuclear coordinates, compute and store the one + * electron matrices. * + * For now, just using loops. In the future, it's an N-body problem but + * probably a very small fraction of the total running time. */ - void FormOrthogonalizingMatrix_() { + void ComputeOneElectronMatrices_(); + + + /** + * Create the matrix S^{-1/2} using the eigenvector decomposition. * + */ + void FormChangeOfBasisMatrix_() { Vector eigenvalues; Matrix eigenvectors; @@ -185,22 +216,21 @@ class SCFSolver { Matrix lambda_times_u_transpose; la::MulTransBInit(sqrt_lambda, eigenvectors, &lambda_times_u_transpose); - la::MulOverwrite(eigenvectors, lambda_times_u_transpose, &overlap_matrix_); + la::MulInit(eigenvectors, lambda_times_u_transpose, + &change_of_basis_matrix_); - } // FormOrthogonalizingMatrix_ + } // FormChangeOfBasisMatrix_() /** * Compute the density matrix. * * TODO: Consider an SVD or some eigenvalue solver that will find the - * eigenvalues in decending order. + * eigenvalues in ascending order. */ void ComputeDensityMatrix_() { - ArrayList occupied_orbitals; - - FillOrbitals_(&occupied_orbitals); + FillOrbitals_(); // I MUST find a smarter way to do this for large problems // This could also probably be a separate function @@ -210,23 +240,28 @@ class SCFSolver { density_row++) { // Columns of density matrix - for (index_t density_column = 0; + for (index_t density_column = density_row; density_column < number_of_basis_functions_; density_column++) { // Occupied orbitals double this_sum = 0.0; for (index_t occupied_index = 0; - occupied_index < occupied_orbitals.size(); occupied_index++) { + occupied_index < occupied_indices_.size(); occupied_index++) { // By multiplying by 2, I'm assuming all the orbitals are full this_sum = this_sum + ( 2 * coefficient_matrix_.ref( - density_row, occupied_orbitals[occupied_index]) * + density_row, occupied_indices_[occupied_index]) * coefficient_matrix_.ref( - density_column, occupied_orbitals[occupied_index])); + density_column, occupied_indices_[occupied_index])); } // occupied_index + // I think this is necessary, but not sure + if (likely(density_row != density_column)) { + this_sum = 2 * this_sum; + } + // Computing the frobenius norm of the difference between this // iteration's density matrix and the previous one for testing // convergence @@ -237,11 +272,18 @@ class SCFSolver { } density_matrix_.set(density_row, density_column, this_sum); + if (density_row != density_column) { + density_matrix_.set(density_column, density_row, this_sum); + } } // density_column } //density_row + printf("density_matrix_norm: %g\n", density_matrix_frobenius_norm_); + + density_matrix_.PrintDebug(); + } // ComputeDensityMatrix_ /** @@ -266,7 +308,7 @@ class SCFSolver { #endif // 3. Find the untransformed eigenvector matrix - la::MulOverwrite(overlap_matrix_, coefficients_prime, + la::MulOverwrite(change_of_basis_matrix_, coefficients_prime, &coefficient_matrix_); } // DiagonalizeFockMatrix_ @@ -276,52 +318,48 @@ class SCFSolver { * Determine the K/2 lowest energy orbitals. * * TODO: If K is odd, then the last entry here is the orbital that should - * have one electron. + * have one electron. I think the closed-shell RHF formulation I'm using + * forbids an odd number of electons. * * I should use the built in C++ iterator-driven routines. Ryan suggested * that I write iterators for ArrayLists, since that would open up a lot of * functionality, including sorting. */ - void FillOrbitals_(ArrayList* indices) { - - index_t number_to_fill = (index_t)ceil((double)number_of_electrons_/2); - - indices->Init(number_to_fill); + void FillOrbitals_() { double max_energy_kept = -DBL_INF; index_t next_to_go = 0; - for (index_t i = 0; i < number_of_basis_functions_; i++) { + for (index_t i = 0; i < number_to_fill_; i++) { - if (unlikely(i < number_to_fill)) { - - (*indices)[i] = i; - if (energy_vector_[i] > max_energy_kept) { - max_energy_kept = energy_vector_[i]; - next_to_go = i; - } - + occupied_indices_[i] = i; + if (energy_vector_[i] > max_energy_kept) { + max_energy_kept = energy_vector_[i]; + next_to_go = i; } - else { + + } + + + for (index_t i = number_to_fill_; i < number_of_basis_functions_; i++) { + + double this_energy = energy_vector_[i]; + if (this_energy < max_energy_kept) { + occupied_indices_[next_to_go] = i; - double this_energy = energy_vector_[i]; - if (this_energy < max_energy_kept) { - (*indices)[next_to_go] = i; - - // Find the new index to throw out - double new_max = -DBL_INF; - next_to_go = -1; - for (index_t j = 0; j < number_to_fill; j++) { - if (energy_vector_[(*indices)[j]] > new_max) { - new_max = energy_vector_[(*indices)[j]]; - next_to_go = j; - } + // Find the new index to throw out + double new_max = -DBL_INF; + next_to_go = -1; + for (index_t j = 0; j < number_to_fill_; j++) { + if (energy_vector_[occupied_indices_[j]] > new_max) { + new_max = energy_vector_[occupied_indices_[j]]; + next_to_go = j; } - max_energy_kept = new_max; - DEBUG_ASSERT(!isinf(max_energy_kept)); - DEBUG_ASSERT(next_to_go >= 0); } - + max_energy_kept = new_max; + DEBUG_ASSERT(!isinf(max_energy_kept)); + DEBUG_ASSERT(next_to_go >= 0); + DEBUG_ASSERT(next_to_go < number_of_basis_functions_); } } @@ -329,48 +367,7 @@ class SCFSolver { } //FillOrbitals_ - /** - * Do step 4a. in Sherrill's notes. This is a key step that could be turned - * into an N-body computation. - * - * TODO: Figure out how to make a matrix that is the density weighted two - * electron integrals. That would eliminate at least two of the for loops. - * - * BUG: I'm counting each integral more than once. - */ - void UpdateFockMatrix_() { - for (index_t mu = 0; mu < number_of_basis_functions_; mu++) { - - for (index_t nu = 0; nu < number_of_basis_functions_; nu++) { - - double new_value = one_electron_integrals_.ref(mu, nu); - - for (index_t rho = 0; rho < number_of_basis_functions_; rho++) { - - for (index_t sigma = 0; sigma <= rho; sigma++) { - - printf("%d,%d,%d,%d\n", mu, nu, rho, sigma); - - index_t first_index = FindIntegralIndex_(mu, nu, rho, sigma); - index_t second_index = FindIntegralIndex_(mu, rho, nu, sigma); - - new_value = new_value + density_matrix_.ref(rho, sigma) * - (2 * two_electron_integrals_[first_index] - - two_electron_integrals_[second_index]); - - } - - } - - fock_matrix_.set(mu, nu, new_value); - - } - - } - - } // UpdateFockMatrix_ - /** * Find the energy of the electrons in the ground state of the current * wavefunction. @@ -383,8 +380,10 @@ class SCFSolver { for (index_t nu = 0; nu < number_of_basis_functions_; nu++) { + // I don't think this is right + // I guess the density matrix needs to multiply the one electron too? total_energy = total_energy + density_matrix_.ref(mu, nu) - * (one_electron_integrals_.ref(mu, nu) + fock_matrix_.ref(mu, nu)); + * (core_matrix_.ref(mu, nu) + fock_matrix_.ref(mu, nu)); } @@ -402,7 +401,6 @@ class SCFSolver { bool is_converged = true; if (unlikely(current_iteration_ == 0)) { - density_matrix_frobenius_norm_ = 0.0; return false; } @@ -429,13 +427,25 @@ class SCFSolver { void TransformFockBasis_() { Matrix orthogonal_transpose_times_fock; - la::MulTransAInit(overlap_matrix_, fock_matrix_, + la::MulTransAInit(change_of_basis_matrix_, fock_matrix_, &orthogonal_transpose_times_fock); la::MulOverwrite(orthogonal_transpose_times_fock, - overlap_matrix_, &fock_matrix_); + change_of_basis_matrix_, &fock_matrix_); } // TransformFockBasis_ + void UpdateFockMatrix_() { + + // Needs to call something from the object, preferably updating it first? + + integrals_.UpdateMatrices(density_matrix_); + + integrals_.ComputeFockMatrix(); + + la::AddOverwrite(core_matrix_, integrals_.FockMatrix(), &fock_matrix_); + + } + /** * Does the SCF iterations to find the HF wavefunction */ @@ -473,30 +483,33 @@ class SCFSolver { } // FindSCFSolution_ + /** + * Returns the nuclear repulsion energy for the nuclei given in + * nuclear_centers_ and nuclear_masses_ + * + * I'm only counting each pair once, which I think is correct. + */ + double ComputeNuclearRepulsion_(); + /** * Sets up the matrices for the SCF iterations */ void Setup_() { - FormOrthogonalizingMatrix_(); + nuclear_repulsion_energy_ = ComputeNuclearRepulsion_(); - // 1. Form the core Fock matrix in transformed basis - // F' = S^{-1/2} H S^{-1/2} - // For now, we assume the initial density matrix is zero - // In the future, I should support non-zero initialization - - fock_matrix_.CopyValues(one_electron_integrals_); + ComputeOneElectronMatrices_(); - TransformFockBasis_(); // fock_matrix_ should now be F' - - // 2. Solve the transformed Fock matrix eigenvalue problem + FormChangeOfBasisMatrix_(); + + fock_matrix_.Alias(core_matrix_); + + TransformFockBasis_(); DiagonalizeFockMatrix_(); ComputeDensityMatrix_(); - density_matrix_frobenius_norm_ = 0.0; - } //Setup_ /** @@ -525,6 +538,16 @@ class SCFSolver { energy_vector_matrix.AliasColVector(energy_vector_); data::Save(energy_vector_file, energy_vector_matrix); + fx_format_result(module_, "density_matrix_norm", "%g", + density_matrix_frobenius_norm_); + + fx_format_result(module_, "num_iterations", "%d", current_iteration_); + + fx_format_result(module_, "total_energy", "%g", + total_energy_[current_iteration_-1]); + + integrals_.OutputFockMatrix(NULL, NULL, NULL, NULL); + } public: @@ -535,9 +558,13 @@ class SCFSolver { */ void ComputeWavefunction() { + fx_timer_start(module_, "SCF_Setup"); Setup_(); + fx_timer_stop(module_, "SCF_Setup"); + fx_timer_start(module_, "SCF_Iterations"); FindSCFSolution_(); + fx_timer_stop(module_, "SCF_Iterations"); OutputResults_(); @@ -545,8 +572,10 @@ class SCFSolver { void PrintMatrices() { - printf("One electron integrals:\n"); - ot::Print(one_electron_integrals_); + // These should be changed to print debug or something + + printf("Core Matrix:\n"); + ot::Print(core_matrix_); printf("Coefficient matrix:\n"); ot::Print(coefficient_matrix_); diff --git a/fastlib2/contrib/march/hf/scf_solver_test.cc b/fastlib2/contrib/march/hf/scf_solver_test.cc index 66354ba15f..62bb3c9ca6 100644 --- a/fastlib2/contrib/march/hf/scf_solver_test.cc +++ b/fastlib2/contrib/march/hf/scf_solver_test.cc @@ -151,89 +151,14 @@ public: } // TestFillOrbitals - void TestFindIntegralIndex() { - index_t test1 = solver_->FindIntegralIndexHelper_(1, 0); - - index_t test2 = solver_->FindIntegralIndexHelper_(5, 4); - - index_t test3 = solver_->FindIntegralIndex_(0, 0, 0, 0); - - index_t test4 = solver_->FindIntegralIndex_(1, 1, 0, 0); - - index_t test5 = solver_->FindIntegralIndex_(1, 0, 1, 0); - - index_t test6 = solver_->FindIntegralIndex_(0, 0, 1, 1); - - index_t test7 = solver_->FindIntegralIndexHelper_(4, 5); - - TEST_ASSERT(test1 == 1); - TEST_ASSERT(test2 == 19); - TEST_ASSERT(test3 == 0); - TEST_ASSERT(test4 == 3); - TEST_ASSERT(test5 == 2); - TEST_ASSERT(test6 == test4); - TEST_ASSERT(test7 == test2); - - NONFATAL("FindIntegralIndex correct.\n"); - - } // TestFindIntegralIndex - void TestDiagonalizeFockMatrix() { NONFATAL("TestDiagonalizeFockMatrix not implemented!\n"); } // TestDiagonalizeFockMatrix - void TestUpdateFockMatrix() { - - Setup(); - - /*index_t test1 = solver_->FindIntegralIndex_(0, 0, 0, 0); - printf("0,0,0,0 = %d\n", test1); - index_t test2 = solver_->FindIntegralIndex_(1, 0, 0, 0); - printf("1,0,0,0 = %d\n", test2); - index_t test3 = solver_->FindIntegralIndex_(1, 0, 1, 0); - printf("1,0,1,0 = %d\n", test3); - index_t test4 = solver_->FindIntegralIndex_(1, 1, 0, 0); - printf("1,1,0,0 = %d\n", test4); - index_t test5 = solver_->FindIntegralIndex_(1, 1, 1, 0); - printf("1,1,1,0 = %d\n", test5); - index_t test6 = solver_->FindIntegralIndex_(1, 1, 1, 1); - printf("1,1,1,1 = %d\n", test6); - */ - - Matrix true_density; - data::Load("density_test.csv", &true_density); - solver_->density_matrix_.CopyValues(true_density); - - solver_->Setup_(); - - solver_->UpdateFockMatrix_(); - - Matrix true_updated_fock; - data::Load("updated_fock_test.csv", &true_updated_fock); - - printf("true_updated_fock\n"); - ot::Print(true_updated_fock); - printf("fock_matrix_\n"); - ot::Print(solver_->fock_matrix_); - - //solver_->PrintMatrices(); - - for (index_t i = 0; i < true_updated_fock.n_rows(); i++) { - for (index_t j = 0; j < true_updated_fock.n_cols(); j++) { - TEST_DOUBLE_APPROX(true_updated_fock.ref(i, j), - solver_->fock_matrix_.ref(i, j), eps); - } - } - - Destruct(); - - NONFATAL("TestUpdateFockMatrix correct.\n"); - - } // TestUpdateFockMatrix - + void TestTestConvergence() { // Not quite sure how to do this one @@ -243,6 +168,8 @@ public: NONFATAL("TestTestConvergence not implemented!\n"); } // TestTestConvergence + + void TestComputeElectronicEnergy() { @@ -268,6 +195,46 @@ public: } // TestComputeElectronicEnergy + void TestComputeOverlapIntegral() { + + Setup(); + + double dist1 = 0.5; + + double test_integral = solver_->ComputeOverlapIntegral_(dist1); + + double correct_integral = 0; + + Destruct(); + + NONFATAL("TestComputeOverlapIntegral NOT IMPLEMENTED.\n"); + + } // TestComputeOverlapIntegral() + + void TestComputeKineticIntegral() { + + NONFATAL("TestComputeKineticIntegral NOT IMPLEMENTED.\n"); + + } // TestComputeKineticIntegral() + + void TestComputeNuclearIntegral() { + + NONFATAL("TestComputeNuclearIntegral NOT IMPLEMENTED.\n"); + + } // TestComputeNuclearIntegral() + + void TestComputeOneElectronMatrices() { + + NONFATAL("TestComputeOneElectronMatrices NOT IMPLEMENTED.\n"); + + } // TestComputeOneElectronMatrices() + + void TestComputeNuclearRepulsion() { + + NONFATAL("TestComputeNuclearRepulsion NOT IMPLEMENTED.\n"); + + } // TestComputeNuclearRepulsion() + void TestAll() { @@ -277,12 +244,8 @@ public: TestFillOrbitals(); - TestFindIntegralIndex(); - TestDiagonalizeFockMatrix(); - TestUpdateFockMatrix(); - TestTestConvergence(); TestComputeElectronicEnergy(); diff --git a/fastlib2/contrib/march/hf/square_tree.h b/fastlib2/contrib/march/hf/square_tree.h index 9b990c687c..188c4d3f65 100644 --- a/fastlib2/contrib/march/hf/square_tree.h +++ b/fastlib2/contrib/march/hf/square_tree.h @@ -38,36 +38,6 @@ private: public: - - // Needs null checks - /*void Init(const SquareIntegralStat& left_left, - const SquareIntegralStat& left_right, - const SquareIntegralStat& right_left, - const SquareIntegralStat& right_right) { - - density_upper_bound_ = max(left_left.density_upper_bound(), - left_right.density_upper_bound()); - density_upper_bound_ = max(density_upper_bound_, - right_left.density_upper_bound()); - density_upper_bound_ = max(density_upper_bound_, - right_right.density_upper_bound()); - - density_lower_bound_ = min(left_left.density_lower_bound(), - left_right.density_lower_bound()); - density_lower_bound_ = min(density_lower_bound_, - right_left.density_lower_bound()); - density_lower_bound_ = min(density_lower_bound_, - right_right.density_lower_bound()); - - entry_lower_bound_ = 0.0; - entry_upper_bound_ = 0.0; - - approximation_val_ = 0.0; - - remaining_references_ = left_left.remaining_references(); - - } // void Init(children) - */ void Init(const SquareIntegralStat& left, const SquareIntegralStat& right) { @@ -88,8 +58,9 @@ public: } // void Init (2 children) void Init(index_t start1, index_t end1, index_t start2, index_t end2, - const Matrix& density) { + index_t num_funs) { + /* double min_density = DBL_MAX; double max_density = -DBL_MAX; @@ -113,13 +84,18 @@ public: density_lower_bound_ = min_density; DEBUG_ASSERT(density_upper_bound_ > -DBL_MAX); DEBUG_ASSERT(density_lower_bound_ < DBL_MAX); + */ + + density_upper_bound_ = DBL_MAX; + density_lower_bound_ = -DBL_MAX; + entry_upper_bound_ = 0.0; entry_lower_bound_ = 0.0; approximation_val_ = 0.0; - remaining_references_ = density.n_cols() * density.n_cols(); + remaining_references_ = num_funs * num_funs; } // void Init(leaf) @@ -220,7 +196,7 @@ class SquareTree { * greater than that of q2 */ void Init(QueryTree1* query1_root, QueryTree2* query2_root, - const Matrix& density) { + index_t num_funs) { query1_ = query1_root; query2_ = query2_root; @@ -236,7 +212,7 @@ class SquareTree { right_child_ = NULL; stat_.Init(query1_->begin(), query1_->end(), query2_->begin(), - query2_->end(), density); + query2_->end(), num_funs); } // I'm assuming that query1_ will always have two significant children @@ -249,8 +225,8 @@ class SquareTree { DEBUG_ASSERT(query1_->right()->end() > query2_->begin()); DEBUG_ASSERT(query1_->left()->end() > query2_->begin()); - left_child_->Init(query1_->left(), query2_, density); - right_child_->Init(query1_->right(), query2_, density); + left_child_->Init(query1_->left(), query2_, num_funs); + right_child_->Init(query1_->right(), query2_, num_funs); stat_.Init(left_child_->stat(), right_child_->stat()); @@ -266,8 +242,8 @@ class SquareTree { left_child_ = new SquareTree(); right_child_ = new SquareTree(); - left_child_->Init(query1_, query2_->left(), density); - right_child_->Init(query1_, query2_->right(), density); + left_child_->Init(query1_, query2_->left(), num_funs); + right_child_->Init(query1_, query2_->right(), num_funs); stat_.Init(left_child_->stat(), right_child_->stat()); @@ -285,7 +261,7 @@ class SquareTree { DEBUG_ASSERT(query1_->end() > query2_->begin()); stat_.Init(query1_->begin(), query1_->end(), query2_->begin(), - query2_->end(), density); + query2_->end(), num_funs); } // Idea: since query2 isn't necessary, go farther down query2 @@ -299,8 +275,8 @@ class SquareTree { DEBUG_ASSERT(query1_->end() > query2_->left()->begin()); DEBUG_ASSERT(query1_->end() > query2_->right()->begin()); - left_child_->Init(query1_, query2_->left(), density); - right_child_->Init(query1_, query2_->right(), density); + left_child_->Init(query1_, query2_->left(), num_funs); + right_child_->Init(query1_, query2_->right(), num_funs); stat_.Init(left_child_->stat(), right_child_->stat()); @@ -315,8 +291,8 @@ class SquareTree { left_child_ = new SquareTree(); right_child_ = new SquareTree(); - left_child_->Init(query1_->left(), query2_, density); - right_child_->Init(query1_->right(), query2_, density); + left_child_->Init(query1_->left(), query2_, num_funs); + right_child_->Init(query1_->right(), query2_, num_funs); stat_.Init(left_child_->stat(), right_child_->stat()); @@ -332,7 +308,7 @@ class SquareTree { DEBUG_ASSERT(query1_->end() > query2_->begin()); stat_.Init(query1_->begin(), query1_->end(), query2_->begin(), - query2_->end(), density); + query2_->end(), num_funs); } // q1 is split twice @@ -346,17 +322,19 @@ class SquareTree { DEBUG_ASSERT(query1_->left()->end() > query2_->begin()); DEBUG_ASSERT(query1_->right()->end() > query2_->begin()); - left_child_->Init(query1_->left(), query2_, density); - right_child_->Init(query1_->right(), query2_, density); + left_child_->Init(query1_->left(), query2_, num_funs); + right_child_->Init(query1_->right(), query2_, num_funs); stat_.Init(left_child_->stat(), right_child_->stat()); } } // q1 higher - + + /* DEBUG_ASSERT(stat_.density_upper_bound() < DBL_MAX); DEBUG_ASSERT(stat_.density_lower_bound() > -DBL_MAX); + */ } // Init() (two-children) diff --git a/fastlib2/contrib/march/hf/square_tree_test.cc b/fastlib2/contrib/march/hf/square_tree_test.cc new file mode 100644 index 0000000000..28b59bb010 --- /dev/null +++ b/fastlib2/contrib/march/hf/square_tree_test.cc @@ -0,0 +1,44 @@ +#include "square_tree.h" +#include "dual_tree_integrals.h" + +class SquareTreeTester { + +private: + + typedef SquareTree + SqrIntegralTree; + + SqrIntegralTree* tree_; + + void Setup_() { + + tree_ = new SqrIntegralTree(); + + + } + + void Destruct_() { + + } + +public: + + void TestAll() { + + Setup_(); + + Destruct_(); + + } // TestAll() + +}; // class SquareTreeTester + + +int main(int argc, char* argv[]) { + + SquareTreeTester tester; + tester.TestAll(); + + return 0; + +} // main() \ No newline at end of file