unit tests for Ryan

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"]
+)

fastlib2/contrib/march/hf/core_test.csv

@@ -0,0 +1,2 @@
+-2.457, -0.985
+-0.985, -1.493

fastlib2/contrib/march/hf/density_test.csv

@@ -0,0 +1,2 @@
+1, 0
+0, 1

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();

fastlib2/contrib/march/hf/hf.h

@@ -8,7 +8,7 @@
 
 #ifndef HF_H
 #define HF_H
-
+#include <fastlib/fastlib.h>
 
 /**
  * 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
+   */
 
   
   

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;
   

fastlib2/contrib/march/hf/kinetic_test.csv

@@ -0,0 +1,2 @@
+1.412, 0.081
+0.081, 0.760

fastlib2/contrib/march/hf/orthogonalizing_test.csv

@@ -0,0 +1,2 @@
+-1.065, -0.217
+-0.217, 1.065

fastlib2/contrib/march/hf/overlap_test.csv

@@ -0,0 +1,2 @@
+1.0, 0.392
+0.392, 1.0

fastlib2/contrib/march/hf/potential_test.csv

@@ -0,0 +1,2 @@
+-3.869, -1.066
+-1.066, -2.303