Trying to get DIIS to work

MartinRJDagleish · Jan 2, 2024 · 38f4079 · 38f4079
1 parent 636bc32
commit 38f4079
Show file tree

Hide file tree

Showing 2 changed files with 10 additions and 23 deletions.
diff --git a/src/calc_type/mod.rs b/src/calc_type/mod.rs
@@ -2,6 +2,7 @@
 use crate::calc_type::rhf::calc_cmp_idx;
 use ndarray::parallel::prelude::*;
 use ndarray::{Array1, Array2, Zip};
+use ndarray_linalg::SolveH;
 use std::ops::{Index, IndexMut};
 
 pub(crate) mod rhf;
@@ -117,26 +118,16 @@ impl DIIS {
         //     let slice = B_matr.slice(s![i + 1..error_set_len, i]).to_shared();
         //     B_matr.slice_mut(s![i, i + 1..error_set_len]).assign(&slice);
         // }
-        //
-        // // * Add langrange multiplier to B_matr_extended
-        // let new_axis_extension_1 = Array2::from_elem((error_set_len, 1), -1.0_f64);
-        // let mut new_axis_extension_2 = Array2::from_elem((1, error_set_len + 1), -1.0_f64);
-        // new_axis_extension_2[[0, error_set_len]] = 0.0_f64;
-        // let mut B_matr_extended = concatenate![Axis(1), B_matr, new_axis_extension_1];
-        // B_matr_extended = concatenate![Axis(0), B_matr_extended, new_axis_extension_2];
-        //
+
         // // * Calculate the coefficients c_vec
-        // let c_vec = B_matr_extended.solveh(&sol_vec).unwrap();
-        // if is_debug {
-        //     println!("c_vec: {:>8.5}", c_vec);
-        // }
+        let c_vec = B_matr.solveh(&sol_vec).unwrap();
         //
         // // * Calculate the new DIIS Fock matrix for new D_matr
         let no_cgtos = self.err_matr_pr_ring_buf[0].shape()[0];
         let mut _F_matr_DIIS = Array2::<f64>::zeros((no_cgtos, no_cgtos));
-        // for i in 0..error_set_len {
-        //     _F_matr_DIIS = _F_matr_DIIS + c_vec[i] * F_matr_set[i];
-        // }
+        for i in 0..error_set_len {
+            _F_matr_DIIS = _F_matr_DIIS + c_vec[i] * &self.F_matr_pr_ring_buf[i];
+        }
 
         _F_matr_DIIS
     }

diff --git a/src/calc_type/rhf.rs b/src/calc_type/rhf.rs
@@ -260,21 +260,17 @@ pub(crate) fn rhf_scf_normal(
             F_matr_pr = S_matr_inv_sqrt.dot(&F_matr).dot(&S_matr_inv_sqrt);
 
             if calc_sett.use_diis {
-                let repl_idx = scf_iter % calc_sett.diis_sett.diis_max;
+                let repl_idx = scf_iter % calc_sett.diis_sett.diis_max - 1; // always start with 0
                 diis.as_mut().unwrap().push_to_ring_buf(
                     &F_matr_pr,
                     &DIIS::calc_FPS_comm(&F_matr_pr, &P_matr, &S_matr),
                     repl_idx,
                 );
 
                 if scf_iter >= calc_sett.diis_sett.diis_min {
-                    let min_val = std::cmp::min(calc_sett.diis_sett.diis_max, scf_iter);
-                    F_matr_pr.assign(&diis.as_ref().unwrap().run_DIIS(min_val));
-                    // F_matr_pr.assign(&run_DIIS(
-                    //     min_val,
-                    //     &diis.take().unwrap().F_matr_pr_ring_buf,
-                    //     &diis.take().unwrap().err_matr_pr_ring_buf,
-                    // ));
+                    println!("*** Using DIIS ***");
+                    let err_set_len = std::cmp::min(calc_sett.diis_sett.diis_max, scf_iter);
+                    F_matr_pr.assign(&diis.as_ref().unwrap().run_DIIS(err_set_len));
                 }
             }