Skip to content

Commit

Permalink
Working on DIIS
Browse files Browse the repository at this point in the history
  • Loading branch information
@MartinRJDagleish
MartinRJDagleish committed Jan 2, 2024
1 parent ac81cf8 commit 636bc32
Show file tree
Hide file tree
Showing 7 changed files with 224 additions and 114 deletions.
1 change: 1 addition & 0 deletions Cargo.lock
View file

Some generated files are not rendered by default. Learn more about how customized files appear on GitHub.

1 change: 1 addition & 0 deletions Cargo.toml
View file
Original file line number Diff line number Diff line change
Expand Up @@ -18,6 +18,7 @@ lazy_static = "1.4.0"
approx = "0.5.1"
boys = "0.1.0"
getset = "0.1.2"
rayon = "1.8.0"

# [package.metadata.docs.rs]
# rustdoc-args = [
Expand Down
110 changes: 98 additions & 12 deletions src/calc_type/mod.rs
View file
Original file line number Diff line number Diff line change
@@ -1,6 +1,7 @@
#![allow(non_snake_case)]
use crate::calc_type::rhf::calc_cmp_idx;
use ndarray::{Array1, Array2};
use ndarray::parallel::prelude::*;
use ndarray::{Array1, Array2, Zip};
use std::ops::{Index, IndexMut};

pub(crate) mod rhf;
Expand All @@ -14,19 +15,19 @@ pub(crate) enum CalcType {
pub struct CalcSettings {
pub max_scf_iter: usize,
pub e_diff_thrsh: f64,
pub rms_p_matr_thrsh: f64,
pub commu_conv_thrsh: f64,
pub use_diis: bool,
pub use_direct_scf: bool,
pub diis_sett: DiisSettings,
}

#[derive(Debug, Clone, Default)]
pub struct DiisSettings {
pub diis_start: usize,
pub diis_end: usize,
pub diis_min: usize,
pub diis_max: usize,
}


#[derive(Debug, Clone)]
pub struct EriArr1 {
eri_arr: Array1<f64>,
Expand All @@ -42,17 +43,104 @@ pub struct SCF {
orb_energies_conv: Array1<f64>,
}


#[derive(Debug, Default)]
struct DIIS {
pub struct DIIS {
// Better approach
F_matr_pr_ring_buf: Vec<Array2<f64>>,
err_matr_pr_ring_buf: Vec<Array2<f64>>,
pub diis_settings: DiisSettings,
pub F_matr_pr_ring_buf: Vec<Array2<f64>>,
pub 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 {
pub fn new(diis_settings: &DiisSettings, matr_size: [usize; 2]) -> Self {
Self {
diis_settings: diis_settings.clone(),
F_matr_pr_ring_buf: vec![Array2::<f64>::zeros(matr_size); diis_settings.diis_max],
err_matr_pr_ring_buf: vec![Array2::<f64>::zeros(matr_size); diis_settings.diis_max],
}
}

#[inline(always)]
/// This is also the DIIS error matrix
pub fn calc_FPS_comm(
F_matr: &Array2<f64>,
P_matr: &Array2<f64>,
S_matr: &Array2<f64>,
// S_matr_inv_sqrt: &Array2<f64>,
) -> Array2<f64> {
F_matr.dot(P_matr).dot(S_matr) - S_matr.dot(P_matr).dot(F_matr)
}

pub fn push_to_ring_buf(&mut self, F_matr: &Array2<f64>, err_matr: &Array2<f64>, idx: usize) {
self.F_matr_pr_ring_buf[idx].assign(F_matr);
self.err_matr_pr_ring_buf[idx].assign(err_matr);
}

// TODO: fix this function
#[inline]
fn run_DIIS(&self, error_set_len: usize) -> Array2<f64> {
let mut B_matr = Array2::<f64>::zeros((error_set_len + 1, error_set_len + 1));
let mut sol_vec = Array1::<f64>::zeros(error_set_len + 1);
sol_vec[error_set_len] = -1.0;

for i in 0..error_set_len {
for j in 0..=i {
B_matr[(i, j)] = Zip::from(&self.err_matr_pr_ring_buf[i])
.and(&self.err_matr_pr_ring_buf[j])
.into_par_iter()
.map(|(err_mat1, err_mat2)| err_mat1 * err_mat2)
.sum();
B_matr[(j, i)] = B_matr[(i, j)];
}
B_matr[(i, error_set_len)] = -1.0;
B_matr[(error_set_len, i)] = -1.0;
}

println!("B_matr: {:>8.5}", B_matr);

// * ACTUALLY: Frobenius inner product of matrices (B_ij = error_matr_i * error_matr_j)
// * OR: flatten error_matr and do dot product
// Zip::indexed(&mut B_matr).par_for_each(|(idx1, idx2), b_val| {
// if idx1 >= idx2 {
// *b_val = Zip::from(self.err_matr_pr_ring_buf[idx1])
// .and(self.err_matr_pr_ring_buf[idx2])
// .into_par_iter()
// .map(|(err_mat1, err_mat2)| err_mat1 * err_mat2)
// .sum();
// }
// });

// for i in 0..error_set_len - 1 {
// 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);
// }
//
// // * 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];
// }

_F_matr_DIIS
}
}

impl EriArr1 {
pub fn new(no_bf: usize) -> Self {
Expand All @@ -64,7 +152,7 @@ impl EriArr1 {
}

impl Index<(usize, usize, usize, usize)> for EriArr1 {
type Output = f64;
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);
let cmp_idx2 = calc_cmp_idx(idx_tup.2, idx_tup.3);
Expand All @@ -76,7 +164,7 @@ impl Index<(usize, usize, usize, usize)> for EriArr1 {

/// This is for the case, where the cmp_idx is already computed
impl Index<usize> for EriArr1 {
type Output = f64;
type Output = f64;
fn index(&self, cmp_idx: usize) -> &f64 {
&self.eri_arr[cmp_idx]
}
Expand All @@ -98,5 +186,3 @@ impl IndexMut<usize> for EriArr1 {
&mut self.eri_arr[cmp_idx]
}
}


Loading

0 comments on commit 636bc32

Please sign in to comment.