apaam/mlpack
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

scf working

Browse Source
This commit is contained in:
Bill March
2008-06-02 14:17:56 +00:00
parent 4ab707e53d
commit d17d101eaa
10 changed files with 875 additions and 556 deletions
Download Patch File Download Diff File Expand all files Collapse all files
fastlib2/contrib/march/hf/build.py
+12 -12
View File
@@ -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"]
)
)
binrule(
name = "hf",
deplibs = [":scf_solver", ":dual_tree_integrals", "fastlib:fastlib"],
sources = ["hf.cc"]
)
fastlib2/contrib/march/hf/dual_tree_integrals.cc
+3 -137
View File
@@ -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, &centers_in);
index_t num_funs = centers_in.n_cols();
la::Scale((double)num_funs, &centers_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<index_t> 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;
}
*/
fastlib2/contrib/march/hf/dual_tree_integrals.h
+151 -42
View File
@@ -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<index_t>* 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<index_t>* 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()
fastlib2/contrib/march/hf/fock_matrix_test.cc
+348
View File
@@ -0,0 +1,348 @@
#include "dual_tree_integrals.h"
#include "naive_fock_matrix.h"
//#include <fastlib/base/test.h>
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<index_t> old_from_new;
}; //class FockMatrixTest
int main(int argc, char* argv[]) {
fx_init(argc, argv);
FockMatrixTest tester;
tester.TestMatrices();
fx_done();
return 0;
} // main
fastlib2/contrib/march/hf/hf.cc
+42 -27
View File
@@ -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, &centers);
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;
fastlib2/contrib/march/hf/hf.h
+3 -36
View File
@@ -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 <fastlib/fastlib.h>
/**
* 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<double> bandwidths_;
ArrayList<double> coefficients_;
public:
ContractedGaussian() {}
~ContractedGaussian() {}
void Init(index_t num, const ArrayList<double>& band,
const ArrayList<double>& coeff) {
number_of_primitives_ = num;
bandwidths_.Copy(band);
coefficients_.Copy(coeff);
}
};
fastlib2/contrib/march/hf/scf_solver.h
+204 -175
View File
@@ -10,14 +10,12 @@
#define SCF_SOLVER_H
#include <fastlib/fastlib.h>
#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<double> total_energy_;
@@ -62,105 +64,134 @@ class SCFSolver {
struct datanode* module_;
DualTreeIntegrals integrals_;
ArrayList<index_t> occupied_indices_;
ArrayList<index_t> 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<index_t>& 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<index_t> 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<index_t>* 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_);
fastlib2/contrib/march/hf/scf_solver_test.cc
+43 -80
View File
@@ -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();
fastlib2/contrib/march/hf/square_tree.h
+25 -47
View File
@@ -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)
fastlib2/contrib/march/hf/square_tree_test.cc
+44
View File
@@ -0,0 +1,44 @@
#include "square_tree.h"
#include "dual_tree_integrals.h"
class SquareTreeTester {
private:
typedef SquareTree<IntegralTree, IntegralTree, SquareIntegralStat>
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()