Working SCF implemented

Cargo.toml

@@ -6,7 +6,7 @@ edition = "2021"
 # See more keys and their definitions at https://doc.rust-lang.org/cargo/reference/manifest.html
 
 [dependencies]
-ndarray = { version = "0.15.6", features = ["blas", "approx-0_5"] }
+ndarray = { version = "0.15.6", features = ["blas", "approx-0_5", "rayon"] }
 ndarray-linalg = { version = "0.16.0", features = ["openblas-system"] }
 blas-src = { version = "0.9", features = ["openblas"] }
 openblas-src = { version = "0.10.8", features = ["cblas", "system"] }

src/basisset/mod.rs

@@ -33,29 +33,32 @@ enum BasisSetVariants {
 ///
 /// ## Notes
 /// - `no_bf` = `no_ao` * 2 if UHF; `no_bf` = `no_ao` if RHF
-#[derive(Debug, CopyGetters, Getters)]
+#[derive(Debug, CopyGetters)]
 pub(crate) struct BasisSet<'a> {
     name: String,
     use_pure_am: bool,
     #[getset(get_copy = "pub")]
     no_ao: usize,
     #[getset(get_copy = "pub")]
     no_bf: usize,
+    #[getset(get_copy = "pub")]
+    no_occ: usize,
     shells: Vec<Shell<'a>>,
-    shell_lens: Vec<usize>,
+    sh_len_offsets: Vec<usize>,
 }
 
 #[derive(Debug, CopyGetters)]
 pub struct Shell<'a> {
     is_pure_am: bool,
     cgtos: Vec<CGTO<'a>>, // == basis funcs.
+    ang_mom_type: AngMomChar,
     #[getset(get_copy = "pub")]
     shell_len: usize,
     center_pos: &'a Atom,
 }
 
 #[allow(clippy::upper_case_acronyms)]
-#[derive(Debug, Getters, Clone)]
+#[derive(Clone, Debug, Getters)]
 pub struct CGTO<'a> {
     pgto_vec: Vec<PGTO>,
     no_pgtos: usize,
@@ -81,7 +84,8 @@ impl<'a> BasisSet<'a> {
         let basis_set_def_total = parser::BasisSetDefTotal::new(basisset_name);
         // Potential speedup: Preallocate vector with correct size
         let mut shells: Vec<Shell<'_>> = Vec::<Shell>::new();
-        let mut shell_lens: Vec<usize> = Vec::<usize>::new();
+        let mut sh_len_offsets: Vec<usize> = vec![0]; // First element is 0
+        let mut sh_offset = 0;
         let mut no_bf = 0;
         for atom in mol.atoms_iter() {
             let basis_set_def_at = basis_set_def_total
@@ -91,7 +95,8 @@ impl<'a> BasisSet<'a> {
             for shell_idx in 0..no_shells {
                 let shell = Shell::new(atom, shell_idx, basis_set_def_at);
                 no_bf += shell.shell_len;
-                shell_lens.push(shell.shell_len);
+                sh_offset += shell.shell_len;
+                sh_len_offsets.push(sh_offset);
                 shells.push(shell);
             }
         }
@@ -101,8 +106,9 @@ impl<'a> BasisSet<'a> {
             name: basisset_name.to_string(),
             no_ao: no_bf,
             no_bf,
+            no_occ: mol.no_elec() / 2,
             shells,
-            shell_lens,
+            sh_len_offsets,
             use_pure_am: false, // hard code for now
         }
     }
@@ -120,10 +126,10 @@ impl<'a> BasisSet<'a> {
     pub fn no_shells(&self) -> usize {
         self.shells.len()
     }
-    
+
     #[inline(always)]
-    pub fn shell_len(&self, sh_idx: usize) -> usize {
-        self.shell_lens[sh_idx]
+    pub fn sh_len_offset(&self, sh_idx: usize) -> usize {
+        self.sh_len_offsets[sh_idx]
     }
 }
 
@@ -170,6 +176,7 @@ impl<'a> Shell<'a> {
         Self {
             is_pure_am: false,
             cgtos,
+            ang_mom_type: curr_ang_mom_char,
             shell_len: no_cgtos,
             center_pos: atom,
         }
@@ -301,10 +308,21 @@ mod tests {
     use super::*;
 
     #[test]
-    fn create_basis_set() {
-        println!("Test create_basis_set");
+    fn test_create_basis_set() {
         let mol = Molecule::new("data/xyz/water90.xyz", 0);
         let test_basis = BasisSet::new("STO-3G", &mol);
         println!("{:?}", test_basis);
     }
+
+    #[test]
+    fn test_sh_len_offsets() {
+        let mol = Molecule::new("data/xyz/water90.xyz", 0);
+        let test_basis = BasisSet::new("STO-3G", &mol);
+
+        println!("sh_len_offsets: {:?}", test_basis.sh_len_offsets);
+        for shell in test_basis.shells.iter() {
+            println!("Shell: {:?}\n", shell);
+            println!("Type:{:?}", shell.ang_mom_type);
+        }
+    }
 }

src/calc_type/mod.rs

@@ -5,11 +5,18 @@ use std::ops::{Index, IndexMut};
 
 mod rhf;
 
+pub(crate) enum CalcType {
+    RHF,
+    UHF,
+    ROHF,
+}
+
 pub struct CalcSettings {
     pub max_scf_iter: usize,
     pub e_diff_thrsh: f64,
     pub commu_conv_thrsh: f64,
     pub use_diis: bool,
+    pub use_direct_scf: bool,
     pub diis_sett: DiisSettings,
 }
 
@@ -19,19 +26,45 @@ pub struct DiisSettings {
     pub diis_max: usize,
 }
 
+
 #[derive(Debug, Clone)]
-pub struct ERI_Arr1 {
+pub struct EriArr1 {
     eri_arr: Array1<f64>,
 }
 
-impl ERI_Arr1 {
+#[derive(Debug, Default)]
+pub struct SCF {
+    tot_scf_iter: usize,
+    pub E_tot_conv: f64,
+    pub E_scf_conv: f64,
+    C_matr_conv: Array2<f64>,
+    P_matr_conv: Array2<f64>, // [ ] TODO: pot. change this to sparse matrix
+    orb_energies_conv: Array1<f64>,
+    diis: DIIS,
+}
+
+
+#[derive(Debug, Default)]
+struct DIIS {
+    // Better approach
+    F_matr_pr_ring_buf: Vec<Array2<f64>>,
+    err_matr_pr_ring_buf: Vec<Array2<f64>>,
+    // Original approach
+    // F_matr_pr_deq: VecDeque<Array2<f64>>,
+    // err_matr_pr_deq: VecDeque<Array2<f64>>,
+}
+
+
+impl EriArr1 {
     pub fn new(no_bf: usize) -> Self {
-        let eri_arr = Array1::<f64>::zeros(no_bf * no_bf * no_bf * no_bf);
+        let idx1 = calc_cmp_idx(no_bf, no_bf);
+        let eri_max_len = calc_cmp_idx(idx1, idx1) + 1;
+        let eri_arr = Array1::<f64>::zeros(eri_max_len);
         Self { eri_arr }
     }
 }
 
-impl Index<(usize, usize, usize, usize)> for ERI_Arr1 {
+impl Index<(usize, usize, usize, usize)> for EriArr1 {
     type Output = f64; 
     fn index(&self, idx_tup: (usize, usize, usize, usize)) -> &f64 {
         let cmp_idx1 = calc_cmp_idx(idx_tup.0, idx_tup.1);
@@ -43,14 +76,14 @@ impl Index<(usize, usize, usize, usize)> for ERI_Arr1 {
 }
 
 /// This is for the case, where the cmp_idx is already computed
-impl Index<usize> for ERI_Arr1 {
+impl Index<usize> for EriArr1 {
     type Output = f64; 
     fn index(&self, cmp_idx: usize) -> &f64 {
         &self.eri_arr[cmp_idx]
     }
 }
 
-impl IndexMut<(usize, usize, usize, usize)> for ERI_Arr1 {
+impl IndexMut<(usize, usize, usize, usize)> for EriArr1 {
     fn index_mut(&mut self, idx_tup: (usize, usize, usize, usize)) -> &mut f64 {
         let cmp_idx1 = calc_cmp_idx(idx_tup.0, idx_tup.1);
         let cmp_idx2 = calc_cmp_idx(idx_tup.2, idx_tup.3);
@@ -61,66 +94,10 @@ impl IndexMut<(usize, usize, usize, usize)> for ERI_Arr1 {
 }
 
 /// This is for the case, where the cmp_idx is already computed
-impl IndexMut<usize> for ERI_Arr1 {
+impl IndexMut<usize> for EriArr1 {
     fn index_mut(&mut self, cmp_idx: usize) -> &mut f64 {
         &mut self.eri_arr[cmp_idx]
     }
 }
 
-// impl IndexMut<usize> for ERI {
-//     fn index_mut(&mut self, i: usize) -> &mut f64 {
-//         match i {
-//             0 => &mut self.x,
-//             1 => &mut self.y,
-//             2 => &mut self.z,
-//             _ => panic!("Index out of bounds for Atom"),
-//         }
-//     }
-// }
 
-#[derive(Debug, Default)]
-struct SCF {
-    tot_scf_iter: usize,
-    E_tot: f64,
-    E_scf: f64,
-    C_matr_final: Array2<f64>,
-    P_matr_final: Array2<f64>, // [ ] TODO: pot. change this to sparse matrix
-    orb_energies_final: Vec<f64>,
-    diis: DIIS,
-}
-
-#[derive(Debug, Default)]
-struct DIIS {
-    // Better approach
-    F_matr_pr_ring_buf: Vec<Array2<f64>>,
-    err_matr_pr_ring_buf: Vec<Array2<f64>>,
-    // Original approach
-    // F_matr_pr_deq: VecDeque<Array2<f64>>,
-    // err_matr_pr_deq: VecDeque<Array2<f64>>,
-}
-
-impl DIIS {
-    fn new() -> DIIS {
-        DIIS {
-            F_matr_pr_ring_buf: Vec::<Array2<f64>>::new(),
-            err_matr_pr_ring_buf: Vec::<Array2<f64>>::new(),
-            // F_matr_pr_deq: VecDeque::<Array2<f64>>::new(),
-            // err_matr_pr_deq: VecDeque::<Array2<f64>>::new(),
-        }
-    }
-}
-
-impl SCF {
-    fn new() -> SCF {
-        let diis = DIIS::new();
-        SCF {
-            tot_scf_iter: 0,
-            E_tot: 0.0_f64,
-            E_scf: 0.0_f64,
-            C_matr_final: Array2::<f64>::zeros((1, 1)),
-            P_matr_final: Array2::<f64>::zeros((1, 1)),
-            orb_energies_final: Vec::<f64>::new(),
-            diis,
-        }
-    }
-}