From a4c2f79ec50fd052960d3cf13e730682fd72eafc Mon Sep 17 00:00:00 2001 From: Bill March Date: Fri, 1 Feb 2008 22:35:48 +0000 Subject: [PATCH] unit tests for Ryan --- fastlib2/contrib/march/hf/build.py | 7 ++ fastlib2/contrib/march/hf/core_test.csv | 2 + fastlib2/contrib/march/hf/density_test.csv | 2 + fastlib2/contrib/march/hf/hf.cc | 11 +- fastlib2/contrib/march/hf/hf.h | 108 ++++++++++++++++-- fastlib2/contrib/march/hf/hf_test.cc | 98 +++++++++++++++- fastlib2/contrib/march/hf/kinetic_test.csv | 2 + .../contrib/march/hf/orthogonalizing_test.csv | 2 + fastlib2/contrib/march/hf/overlap_test.csv | 2 + fastlib2/contrib/march/hf/potential_test.csv | 2 + 10 files changed, 218 insertions(+), 18 deletions(-) create mode 100644 fastlib2/contrib/march/hf/build.py create mode 100644 fastlib2/contrib/march/hf/core_test.csv create mode 100644 fastlib2/contrib/march/hf/density_test.csv create mode 100644 fastlib2/contrib/march/hf/kinetic_test.csv create mode 100644 fastlib2/contrib/march/hf/orthogonalizing_test.csv create mode 100644 fastlib2/contrib/march/hf/overlap_test.csv create mode 100644 fastlib2/contrib/march/hf/potential_test.csv diff --git a/fastlib2/contrib/march/hf/build.py b/fastlib2/contrib/march/hf/build.py new file mode 100644 index 0000000000..834eed9426 --- /dev/null +++ b/fastlib2/contrib/march/hf/build.py @@ -0,0 +1,7 @@ +librule( + name = "hf", + headers = ["hf.h"], + #sources = ["hf.cc"], + deplibs = ["fastlib:fastlib"], + tests = ["hf_test.cc"] +) \ No newline at end of file diff --git a/fastlib2/contrib/march/hf/core_test.csv b/fastlib2/contrib/march/hf/core_test.csv new file mode 100644 index 0000000000..2da2b58cda --- /dev/null +++ b/fastlib2/contrib/march/hf/core_test.csv @@ -0,0 +1,2 @@ +-2.457, -0.985 +-0.985, -1.493 diff --git a/fastlib2/contrib/march/hf/density_test.csv b/fastlib2/contrib/march/hf/density_test.csv new file mode 100644 index 0000000000..f4416b0a6c --- /dev/null +++ b/fastlib2/contrib/march/hf/density_test.csv @@ -0,0 +1,2 @@ +1, 0 +0, 1 diff --git a/fastlib2/contrib/march/hf/hf.cc b/fastlib2/contrib/march/hf/hf.cc index 5861e3c4b5..6693f1a718 100644 --- a/fastlib2/contrib/march/hf/hf.cc +++ b/fastlib2/contrib/march/hf/hf.cc @@ -17,7 +17,7 @@ int main(int argc, char *argv[]) { // 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 @@ -26,18 +26,21 @@ int main(int argc, char *argv[]) { // Should this be in the same, or a different class from the linear system // solver? - + 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(); diff --git a/fastlib2/contrib/march/hf/hf.h b/fastlib2/contrib/march/hf/hf.h index 2ef3be26ad..062ddcfa41 100644 --- a/fastlib2/contrib/march/hf/hf.h +++ b/fastlib2/contrib/march/hf/hf.h @@ -8,7 +8,7 @@ #ifndef HF_H #define HF_H - +#include /** * A class that stores the information for a contracted Gaussian basis function. @@ -18,7 +18,7 @@ */ class ContractedGaussian { - FORBID_ACCIDENTAL_COPY(ContractedGaussian); + FORBID_ACCIDENTAL_COPIES(ContractedGaussian); private: @@ -51,22 +51,41 @@ class ContractedGaussian { * Algorithm class for the basic 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 HFSolver { + friend class HartreeFockTest; - FORBID_ACCIDENTAL_COPY(HFSolver); + FORBID_ACCIDENTAL_COPIES(HFSolver); + private: - Matrix fock_matrix_; - Matrix overlap_matrix_; + // 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 coefficient_vector_; + // This isn't a matrix, it's a rank-four tensor + // I need to figure out what to do with this + Matrix two_electron_integrals_; + + + Matrix one_electron_integrals_; + Matrix kinetic_energy_integrals_; + Matrix potential_energy_integrals_; + + Matrix coefficient_matrix_; + Matrix overlap_matrix_; + Matrix density_matrix_; + Matrix fock_matrix_; + Vector energy_vector_; index_t number_of_basis_functions_; index_t number_of_electrons_; + double nuclear_repulsion_energy_; public: @@ -75,20 +94,85 @@ class HFSolver { ~HFSolver() {} /** - * Initialize the class with const references to the Fock matrix and the + * Initialize the class with const references to the electron matrices and the * overlap matrix, both of which should have been computed already. */ - void Init(const Matrix& fock_in, const Matrix& overlap_in) { + void Init(double nuclear_energy, const Matrix& overlap_in, + const Matrix& kinetic_in, const Matrix& potential_in, + const Matrix& two_electron_in, index_t num_electrons) { - fock_matrix_.Copy(fock_in); + nuclear_repulsion_energy_ = nuclear_energy; + 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); + two_electron_integrals_.Copy(two_electron_in); - // I think this will be necessary for the lapack routines - coefficient_vector_.Init(number_of_basis_functions_); - energy_vector_.Init(number_of_basis_functions_); + number_of_basis_functions_ = overlap_matrix_.n_cols(); + // Form the core Hamiltonian + la::AddInit(kinetic_energy_integrals_, potential_energy_integrals_, + &one_electron_integrals_); + + } // Init + /** + * Create the matrix S^{-1/2} using the Schur decomposition. Overwrites + * overlap_matrix_ with S^{-1/2}. + * + * TODO: go over the linear algrebra and make sure it is efficient + */ + void FormOrthogonalizingMatrix() { + + // Form the orthogonalizing matrix S^{-1/2} + // Should change this to SchurExpert eventually + Vector real_eigenvalues; + Vector imaginary_eigenvalues; + Matrix schur_form; + Matrix schur_vectors; + + success_t diagonalize = la::SchurInit(overlap_matrix_, &real_eigenvalues, + &imaginary_eigenvalues, &schur_form, &schur_vectors); + + if (diagonalize == SUCCESS_FAIL) { + // Need to handle this better + FATAL("Schur Decomposition Failed\n"); + } + +#ifdef DEBUG + // Check that the eigenvalues are all real + for (index_t i = 0; i < imaginary_eigenvalues.length(); i++) { + DEBUG_ASSERT(imaginary_eigenvalues[i] == 0.0); + } + + // Also check that the Schur form is strictly diagonal + for (index_t i = 0; i < schur_form.n_rows(); i++) { + for (index_t j = (i+1); j < schur_form.n_cols(); j++) { + DEBUG_ASSERT(schur_form.ref(i,j) == 0.0); + } + } +#endif + + // Compute lambda^{-1/2} + for (index_t i = 0; i < real_eigenvalues.length(); i++) { + real_eigenvalues[i] = 1/sqrt(real_eigenvalues[i]); + } + + Matrix sqrt_lambda; + sqrt_lambda.InitDiagonal(real_eigenvalues); + Matrix lambda_times_u_transpose; + la::MulTransBInit(sqrt_lambda, schur_vectors, &lambda_times_u_transpose); + la::MulOverwrite(schur_vectors, lambda_times_u_transpose, &overlap_matrix_); + + } // FormOrthogonalizingMatrix + + + /** + * Compute an initial density matrix + */ diff --git a/fastlib2/contrib/march/hf/hf_test.cc b/fastlib2/contrib/march/hf/hf_test.cc index 591e198d65..d5b846b1ef 100644 --- a/fastlib2/contrib/march/hf/hf_test.cc +++ b/fastlib2/contrib/march/hf/hf_test.cc @@ -6,20 +6,114 @@ * Unit tests for the Hartree-Fock code. */ -#include "base/test.h" +#include "fastlib/base/test.h" #include "hf.h" +/** + * Unit test class for Hartree-Fock code. + */ class HartreeFockTest { + // Use the example from Leach's book + +public: + static const index_t num_electrons = 2; + + void Init() { + + solver_ = new HFSolver(); + + Matrix overlap; + Matrix kinetic; + Matrix potential; + Matrix two_electron; + + Matrix core; + + double nuclear_energy = 0.0; + + data::Load("overlap_test.csv", &overlap); + data::Load("kinetic_test.csv", &kinetic); + data::Load("potential_test.csv", &potential); + //data::Load("two_electron_test.csv", &two_electron); + two_electron.Init(2, 2); + + data::Load("core_test.csv", &core); + + solver_->Init(nuclear_energy, overlap, kinetic, potential, two_electron, + num_electrons); + + for (index_t i = 0; i < core.n_rows(); i++) { + for (index_t j = 0; j < core.n_cols(); j++) { + TEST_DOUBLE_APPROX(core.ref(i, j), + solver_->one_electron_integrals_.ref(i, j), 0.0001); + } + } + + + } + + void Destruct() { + + delete solver_; + + } + + void TestOrthogonalizingMatrix() { + Init(); + + solver_->FormOrthogonalizingMatrix(); + + Matrix true_orthogonal; + data::Load("orthogonalizing_test.csv", &true_orthogonal); + + for (index_t i = 0; i < true_orthogonal.n_rows(); i++) { + for (index_t j = 0; j < true_orthogonal.n_cols(); j++) { + TEST_DOUBLE_APPROX(true_orthogonal.ref(i, j), + solver_->overlap_matrix_.ref(i, j), 0.0001); + } + } + + Destruct(); + + NONFATAL("Orthogonal matrix correct.\n"); + } + + void TestAll() { + + TestOrthogonalizingMatrix(); + + NONFATAL("All tests passed\n"); + + } + +private: + + HFSolver* solver_; }; +#if 0 +class Bob { +public: + static const double x = .3; + double y; + Bob() { + y = .3; + } + double foo() { + return y += x; + } +}; +#endif int main(int argc, char *argv[]) { + - + HartreeFockTest tester; + tester.TestAll(); return 0; diff --git a/fastlib2/contrib/march/hf/kinetic_test.csv b/fastlib2/contrib/march/hf/kinetic_test.csv new file mode 100644 index 0000000000..5942e8693d --- /dev/null +++ b/fastlib2/contrib/march/hf/kinetic_test.csv @@ -0,0 +1,2 @@ +1.412, 0.081 +0.081, 0.760 diff --git a/fastlib2/contrib/march/hf/orthogonalizing_test.csv b/fastlib2/contrib/march/hf/orthogonalizing_test.csv new file mode 100644 index 0000000000..a7b84a93a5 --- /dev/null +++ b/fastlib2/contrib/march/hf/orthogonalizing_test.csv @@ -0,0 +1,2 @@ +-1.065, -0.217 +-0.217, 1.065 diff --git a/fastlib2/contrib/march/hf/overlap_test.csv b/fastlib2/contrib/march/hf/overlap_test.csv new file mode 100644 index 0000000000..377d520708 --- /dev/null +++ b/fastlib2/contrib/march/hf/overlap_test.csv @@ -0,0 +1,2 @@ +1.0, 0.392 +0.392, 1.0 diff --git a/fastlib2/contrib/march/hf/potential_test.csv b/fastlib2/contrib/march/hf/potential_test.csv new file mode 100644 index 0000000000..57a02be577 --- /dev/null +++ b/fastlib2/contrib/march/hf/potential_test.csv @@ -0,0 +1,2 @@ +-3.869, -1.066 +-1.066, -2.303 \ No newline at end of file