From 65b99c2c4201c041ac15ae03bb0d785e9358b069 Mon Sep 17 00:00:00 2001 From: MRJD Date: Sun, 7 Jan 2024 19:56:56 +0100 Subject: [PATCH] Added HFMatrics but before rewrite --- src/calc_type/mod.rs | 68 +++++++- src/calc_type/rhf.rs | 308 +++++++++++++++++++++++++++++++++-- src/calc_type/uhf.rs | 10 +- src/main.rs | 8 +- src/print_utils/mod.rs | 11 +- src/print_utils/print_rhf.rs | 22 --- src/print_utils/print_scf.rs | 28 ++++ 7 files changed, 403 insertions(+), 52 deletions(-) delete mode 100644 src/print_utils/print_rhf.rs create mode 100644 src/print_utils/print_scf.rs diff --git a/src/calc_type/mod.rs b/src/calc_type/mod.rs index 3c968a3..4d4750e 100644 --- a/src/calc_type/mod.rs +++ b/src/calc_type/mod.rs @@ -5,16 +5,27 @@ use ndarray::{Array1, Array2, Zip}; use ndarray_linalg::SolveH; use std::ops::{Index, IndexMut}; -pub(crate) mod rhf; -pub(crate) mod uhf; -pub(crate) mod guess; +pub mod rhf; +pub mod uhf; +pub mod guess; -pub(crate) enum CalcType { +pub(crate) enum Reference { RHF, UHF, ROHF, } +// trait HF { +// fn run_scf( +// calc_sett: &CalcSettings, +// exec_times: &mut crate::print_utils::ExecTimes, +// basis: &crate::basisset::BasisSet, +// mol: &crate::molecule::Molecule, +// ) -> SCF; +// +// } +// +#[derive(Debug)] pub struct CalcSettings { pub max_scf_iter: usize, pub e_diff_thrsh: f64, @@ -35,6 +46,55 @@ pub struct EriArr1 { eri_arr: Array1, } +pub struct HFMatrices { + //--------------------------------- + /// 1. Constant matrices after the first SCF iteration + S_matr: Array2, + S_matr_inv_sqrt: Array2, + T_matr: Array2, + V_matr: Array2, + H_core_matr: Array2, + //--------------------------------- + + //--------------------------------- + /// 2. ERI + //--------------------------------- + // Option -> Direct vs. indirect SCF + eri_arr: Option, + //--------------------------------- + + //--------------------------------- + // 3. Mutable matrices during iterations + //--------------------------------- + // 3.1 Current iteration + //--------------------------------- + C_matr_alpha: Array2, + P_matr_alpha: Array2, // [ ] TODO: pot. change this to sparse matrix + // Options -> RHF, UHF, ROHF + C_matr_beta: Option>, + P_matr_beta: Option>, // [ ] TODO: pot. change this to sparse matrix + + F_matr_alpha: Array2, + F_matr_pr_alpha: Array2, + // Options -> RHF, UHF, ROHF + F_matr_beta: Option>, + F_matr_pr_beta: Option>, + //--------------------------------- + + //--------------------------------- + // 3.2 Previous iteration + //--------------------------------- + P_matr_prev_alpha: Array2, + P_matr_prev_beta: Option>, + + //--------------------------------- + // 3.3 For direct SCF + //--------------------------------- + delta_P_matr_alpha: Option>, + // for UHF + delta_P_matr_beta: Option>, +} + #[derive(Debug, Default)] pub struct SCF { tot_scf_iter: usize, diff --git a/src/calc_type/rhf.rs b/src/calc_type/rhf.rs index c11245c..3600767 100644 --- a/src/calc_type/rhf.rs +++ b/src/calc_type/rhf.rs @@ -1,26 +1,31 @@ use super::{CalcSettings, SCF}; -use crate::basisset::BasisSet; -use crate::calc_type::{EriArr1, DIIS}; -use crate::mol_int_and_deriv::te_int::calc_schwarz_est_int; -use crate::mol_int_and_deriv::{ - oe_int::{calc_kinetic_int_cgto, calc_overlap_int_cgto, calc_pot_int_cgto}, - te_int::calc_ERI_int_cgto, +use crate::{ + basisset::BasisSet, + calc_type::{EriArr1, DIIS}, + mol_int_and_deriv::{ + oe_int::{calc_kinetic_int_cgto, calc_overlap_int_cgto, calc_pot_int_cgto}, + te_int::{calc_ERI_int_cgto, calc_schwarz_est_int}, + }, + molecule::Molecule, + print_utils::{fmt_f64, print_scf::print_scf_header_and_settings, ExecTimes}, }; -use crate::molecule::Molecule; -use crate::print_utils::{fmt_f64, print_rhf::print_scf_header_and_settings, ExecTimes}; -use ndarray::linalg::general_mat_mul; -use ndarray::parallel::prelude::*; -use ndarray::{s, Array, Array1, Array2, Zip}; -use ndarray_linalg::{Eigh, UPLO, InverseH}; + +use ndarray::{linalg::general_mat_mul, parallel::prelude::*, prelude::*, Zip}; +use ndarray_linalg::{Eigh, InverseH, UPLO}; #[derive(Debug, Clone, Default)] pub(crate) struct RHF { // Temp. matrices P_matr_old: Array2, - E_scf_old: f64, + // E_scf_old: f64, } impl RHF { + + // pub init_calc() { + // + // } + /// ### Description /// Calculate the 1e integrals for the given basis set and molecule. /// Returns a tuple of the overlap, kinetic and potential energy matrices. @@ -124,6 +129,7 @@ impl RHF { eri } + /// This is the main RHF SCF function to be called and run the RHF SCF calculation /// ## Options: /// - DIIS @@ -135,7 +141,7 @@ impl RHF { basis: &BasisSet, mol: &Molecule, ) -> SCF { - print_scf_header_and_settings(&calc_sett, crate::calc_type::CalcType::RHF); + print_scf_header_and_settings(&calc_sett, crate::calc_type::Reference::RHF); const SHOW_ALL_CONV_CRIT: bool = false; let mut is_scf_conv = false; @@ -392,6 +398,278 @@ impl RHF { scf } + + + + ////////////////////// BACKUP before rewrite + // /// This is the main RHF SCF function to be called and run the RHF SCF calculation + // /// ## Options: + // /// - DIIS + // /// - direct vs. indirect SCF + // #[allow(unused)] + // pub fn run_scf( + // calc_sett: &CalcSettings, + // exec_times: &mut ExecTimes, + // basis: &BasisSet, + // mol: &Molecule, + // ) -> SCF { + // print_scf_header_and_settings(&calc_sett, crate::calc_type::Reference::RHF); + // const SHOW_ALL_CONV_CRIT: bool = false; + // + // let mut is_scf_conv = false; + // let mut scf = SCF::default(); + // let mut diis: Option = if calc_sett.use_diis { + // Some(DIIS::new( + // &calc_sett.diis_sett, + // [basis.no_bf(), basis.no_bf()], + // )) + // } else { + // None + // }; + // + // let V_nuc: f64 = if mol.no_atoms() > 100 { + // mol.calc_core_potential_par() + // } else { + // mol.calc_core_potential_ser() + // }; + // + // println!("Calculating 1e integrals ..."); + // let (S_matr, H_core) = Self::calc_1e_int_matrs(basis, mol); + // println!("FINSIHED calculating 1e integrals ..."); + // + // let mut eri_opt; + // let schwarz_est_matr; + // if calc_sett.use_direct_scf { + // eri_opt = None; + // + // println!("Calculating Schwarz int estimates ..."); + // schwarz_est_matr = Some(calc_schwarz_est_int(basis)); + // println!("FINISHED Schwarz int estimates ..."); + // } else { + // schwarz_est_matr = None; + // + // println!("Calculating 2e integrals ..."); + // eri_opt = Some(Self::calc_2e_int_matr(basis)); + // println!("FINSIHED calculating 2e integrals ..."); + // } + // + // let S_matr_inv_sqrt = Self::inv_ssqrt(&S_matr, UPLO::Upper); + // + // // Init matrices for SCF loop + // let mut C_matr_AO; + // let mut C_matr_MO; + // let mut orb_ener; + // let mut E_scf_prev = 0.0; + // + // let mut P_matr = Array2::::zeros((basis.no_bf(), basis.no_bf())); + // let mut P_matr_old = P_matr.clone(); + // let mut delta_P_matr = None; + // // if calc_sett.use_direct_scf { + // // delta_P_matr = Some(P_matr.clone()); + // // } else { + // // delta_P_matr = None; + // // } + // // let mut delta_P_matr: Option> = Some(P_matr.clone()); + // + // let mut F_matr_pr; + // let mut diis_str = ""; + // + // // Initial guess -> H_core + // let mut F_matr = H_core.clone(); + // + // // Print SCF iteration Header + // match SHOW_ALL_CONV_CRIT { + // true => { + // println!( + // "{:>3} {:^20} {:^20} {:^20} {:^20} {:^20}", + // "Iter", "E_scf", "E_tot", "RMS(P)", "ΔE", "RMS(|FPS - SPF|)" + // ); + // } + // false => { + // println!( + // "{:>3} {:^20} {:^20} {:^20} {:^20}", + // "Iter", "E_scf", "E_tot", "ΔE", "RMS(|FPS - SPF|)" + // ); + // } + // } + // for scf_iter in 0..=calc_sett.max_scf_iter { + // if scf_iter == 0 { + // F_matr_pr = S_matr_inv_sqrt.dot(&F_matr).dot(&S_matr_inv_sqrt); + // + // (orb_ener, C_matr_MO) = F_matr_pr.eigh(UPLO::Upper).unwrap(); + // C_matr_AO = S_matr_inv_sqrt.dot(&C_matr_MO); + // + // Self::calc_P_matr_rhf(&mut P_matr, &C_matr_AO, basis.no_occ()); + // ///// Quick test if this is the correct SAD P_matr for H2O RHF in STO-3G basis + // // P_matr = arr2(&[ + // // [ + // // 2.0000000000, + // // 0.0000000000, + // // 0.0000000000, + // // 0.0000000000, + // // 0.0000000000, + // // 0.0000000000, + // // 0.0000000000, + // // ], + // // [ + // // 0.0000000000, + // // 2.0000000000, + // // 0.0000000000, + // // 0.0000000000, + // // 0.0000000000, + // // 0.0000000000, + // // 0.0000000000, + // // ], + // // [ + // // 0.0000000000, + // // 0.0000000000, + // // 1.3333333333, + // // 0.0000000000, + // // 0.0000000000, + // // 0.0000000000, + // // 0.0000000000, + // // ], + // // [ + // // 0.0000000000, + // // 0.0000000000, + // // 0.0000000000, + // // 1.3333333333, + // // 0.0000000000, + // // 0.0000000000, + // // 0.0000000000, + // // ], + // // [ + // // 0.0000000000, + // // 0.0000000000, + // // 0.0000000000, + // // 0.0000000000, + // // 1.3333333333, + // // 0.0000000000, + // // 0.0000000000, + // // ], + // // [ + // // 0.0000000000, + // // 0.0000000000, + // // 0.0000000000, + // // 0.0000000000, + // // 0.0000000000, + // // 1.0000000000, + // // 0.0000000000, + // // ], + // // [ + // // 0.0000000000, + // // 0.0000000000, + // // 0.0000000000, + // // 0.0000000000, + // // 0.0000000000, + // // 0.0000000000, + // // 1.0000000000, + // // ], + // // ]); + // if calc_sett.use_direct_scf { + // delta_P_matr = Some(P_matr.clone()); + // } + // } else { + // /// direct or indirect scf + // match eri_opt { + // Some(ref eri) => { + // Self::calc_new_F_matr_ind_scf_rhf(&mut F_matr, &H_core, &P_matr, eri); + // } + // None => { + // Self::calc_new_F_matr_dir_scf_rhf( + // &mut F_matr, + // delta_P_matr.as_ref().unwrap(), + // schwarz_est_matr.as_ref().unwrap(), + // basis, + // ); + // } + // } + // let E_scf_curr = Self::calc_E_scf_rhf(&P_matr, &H_core, &F_matr); + // scf.E_tot_conv = E_scf_curr + V_nuc; + // let fps_comm = DIIS::calc_FPS_comm(&F_matr, &P_matr, &S_matr); + // + // 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 - 1) % calc_sett.diis_sett.diis_max; // always start with 0 + // let err_matr = S_matr_inv_sqrt.dot(&fps_comm).dot(&S_matr_inv_sqrt); + // diis.as_mut() + // .unwrap() + // .push_to_ring_buf(&F_matr_pr, &err_matr, repl_idx); + // + // if scf_iter >= calc_sett.diis_sett.diis_min { + // let err_set_len = std::cmp::min(calc_sett.diis_sett.diis_max, scf_iter); + // F_matr_pr = diis.as_ref().unwrap().run_DIIS(err_set_len); + // diis_str = "DIIS"; + // } + // } + // + // (orb_ener, C_matr_MO) = F_matr_pr.eigh(UPLO::Upper).unwrap(); + // C_matr_AO = S_matr_inv_sqrt.dot(&C_matr_MO); + // + // let delta_E = E_scf_curr - E_scf_prev; + // let rms_comm_val = (fps_comm.par_iter().map(|x| x * x).sum::() + // / fps_comm.len() as f64) + // .sqrt(); + // if SHOW_ALL_CONV_CRIT { + // let rms_p_val = Self::calc_rms_2_matr(&P_matr, &P_matr_old.clone()); + // println!( + // "{:>3} {:>20.12} {:>20.12} {:>20.12} {:>20.12} {:>20.12}", + // scf_iter, E_scf_curr, scf.E_tot_conv, rms_p_val, delta_E, rms_comm_val + // ); + // } else { + // println!( + // "{:>3} {:>20.12} {:>20.12} {} {} {:>10} ", + // scf_iter, + // E_scf_curr, + // scf.E_tot_conv, + // fmt_f64(delta_E, 20, 8, 2), + // fmt_f64(rms_comm_val, 20, 8, 2), + // diis_str + // ); + // diis_str = ""; + // } + // + // if (delta_E.abs() < calc_sett.e_diff_thrsh) + // && (rms_comm_val < calc_sett.commu_conv_thrsh) + // { + // scf.tot_scf_iter = scf_iter; + // scf.E_scf_conv = E_scf_curr; + // scf.C_matr_conv_alph = C_matr_AO; + // scf.P_matr_conv_alph = P_matr; + // scf.C_matr_conv_beta = None; + // scf.P_matr_conv_beta = None; + // scf.orb_E_conv_alph = orb_ener.clone(); + // println!("\nSCF CONVERGED!\n"); + // is_scf_conv = true; + // break; + // } else if scf_iter == calc_sett.max_scf_iter { + // println!("\nSCF DID NOT CONVERGE!\n"); + // break; + // } + // E_scf_prev = E_scf_curr; + // P_matr_old = P_matr.clone(); + // Self::calc_P_matr_rhf(&mut P_matr, &C_matr_AO, basis.no_occ()); + // if calc_sett.use_direct_scf { + // delta_P_matr = Some((&P_matr - &P_matr_old).to_owned()); + // } + // } + // } + // + // if is_scf_conv { + // println!("{:*<55}", ""); + // println!("* {:^51} *", "FINAL RESULTS"); + // println!("{:*<55}", ""); + // println!(" {:^50}", "RHF SCF (in a.u.)"); + // println!(" {:=^50} ", ""); + // println!(" {:<25}{:>25}", "Total SCF iterations:", scf.tot_scf_iter); + // println!(" {:<25}{:>25.18}", "Final SCF energy:", scf.E_scf_conv); + // println!(" {:<25}{:>25.18}", "Final tot. energy:", scf.E_tot_conv); + // println!("{:*<55}", ""); + // } + // scf + // } + // fn calc_P_matr_rhf(P_matr: &mut Array2, C_matr: &Array2, no_occ: usize) { let C_occ = C_matr.slice(s![.., ..no_occ]); general_mat_mul(2.0_f64, &C_occ, &C_occ.t(), 0.0_f64, P_matr) @@ -550,7 +828,7 @@ fn rhf_scf_linscal( basis: &BasisSet, mol: &Molecule, ) { - print_scf_header_and_settings(calc_sett, crate::calc_type::CalcType::RHF); + print_scf_header_and_settings(calc_sett, crate::calc_type::Reference::RHF); const SHOW_ALL_CONV_CRIT: bool = false; let mut is_scf_conv = false; diff --git a/src/calc_type/uhf.rs b/src/calc_type/uhf.rs index cb69f49..6400097 100644 --- a/src/calc_type/uhf.rs +++ b/src/calc_type/uhf.rs @@ -7,12 +7,18 @@ use crate::{ }, mol_int_and_deriv::te_int::calc_schwarz_est_int, molecule::Molecule, - print_utils::{fmt_f64, print_rhf::print_scf_header_and_settings}, + print_utils::{fmt_f64, print_scf::print_scf_header_and_settings}, }; use ndarray::parallel::prelude::*; use ndarray::{linalg::general_mat_mul, s, Array1, Array2, Zip}; use ndarray_linalg::{Eigh, UPLO}; +pub(crate) struct UHF { + // Temp. matrices + P_matr_old: Array2, + E_scf_old: f64, +} + #[allow(unused)] pub(crate) fn uhf_scf_normal( calc_sett: &CalcSettings, @@ -20,7 +26,7 @@ pub(crate) fn uhf_scf_normal( basis: &BasisSet, mol: &Molecule, ) -> SCF { - print_scf_header_and_settings(calc_sett, crate::calc_type::CalcType::UHF); + print_scf_header_and_settings(calc_sett, crate::calc_type::Reference::UHF); let mut is_scf_conv = false; let mut scf = SCF::default(); diff --git a/src/main.rs b/src/main.rs index 4b59ab8..9402802 100644 --- a/src/main.rs +++ b/src/main.rs @@ -11,7 +11,7 @@ mod mol_int_and_deriv; mod molecule; mod print_utils; -use crate::{calc_type::CalcType, print_utils::print_logo}; +use crate::{calc_type::Reference, print_utils::print_header_logo}; use basisset::BasisSet; use calc_type::{rhf::RHF, uhf::uhf_scf_normal, CalcSettings, DiisSettings}; use molecule::Molecule; @@ -23,7 +23,7 @@ fn main() { let mut exec_times = print_utils::ExecTimes::new(); exec_times.start("Total"); - print_logo(); + print_header_logo(); exec_times.start("Molecule"); let mol = Molecule::new("data/xyz/water90.xyz", 0); @@ -42,7 +42,7 @@ fn main() { // Calculation type // exec_times.start("RHF noDIIS"); - let _calc_type = CalcType::RHF; + let _calc_type = Reference::RHF; let calc_sett = CalcSettings { max_scf_iter: 100, @@ -65,5 +65,5 @@ fn main() { //################################## //### FOOTER ### //################################## - exec_times.print(); + exec_times.print_wo_order(); } diff --git a/src/print_utils/mod.rs b/src/print_utils/mod.rs index a59a12e..8a30f30 100644 --- a/src/print_utils/mod.rs +++ b/src/print_utils/mod.rs @@ -1,8 +1,8 @@ use std::{collections::HashMap, time::Instant}; -pub mod print_rhf; +pub mod print_scf; -pub fn print_logo() { +pub fn print_header_logo() { // const HEADER_V1: &str = r#" // ___ ___ _ _____ _ ______ _____ // | \/ | | | | ___| | | ___ \/ ___| @@ -80,9 +80,10 @@ pub(crate) fn fmt_f64(num: f64, width: usize, precision: usize, exp_pad: usize) let exp = num.split_off(num.find('e').unwrap()); let (sign, exp) = if exp.starts_with("e-") { - ('-', &exp[2..]) + // ('-', &exp[2..]) + ('-', exp.strip_prefix("e-").unwrap()) } else { - ('+', &exp[1..]) + ('+', exp.strip_prefix('e').unwrap()) }; num.push_str(&format!("e{}{:0>pad$}", sign, exp, pad = exp_pad)); @@ -112,7 +113,7 @@ impl ExecTimes { } } - pub fn print(&self) { + pub fn print_wo_order(&self) { print_header_for_section("Execution times"); for (name, timings) in &self.timings_map { let exec_time = timings[1].duration_since(timings[0]); diff --git a/src/print_utils/print_rhf.rs b/src/print_utils/print_rhf.rs deleted file mode 100644 index 1f9c957..0000000 --- a/src/print_utils/print_rhf.rs +++ /dev/null @@ -1,22 +0,0 @@ -use crate::calc_type::{CalcSettings, CalcType}; - -pub(crate) fn print_scf_header_and_settings(calc_sett: &CalcSettings, calc_type: CalcType) { - println!("{:=>35}", ""); - match calc_type { - CalcType::RHF => println!("{:^35}", "RHF SCF"), - CalcType::UHF => println!("{:^35}", "UHF SCF"), - CalcType::ROHF => println!("{:^35}", "ROHF SCF") - } - println!("{:=>35}", ""); - - println!("{:-20}", ""); - println!("SCF settings:"); - println!(" {:<20} {:>10e}", "ΔE THRESH", calc_sett.e_diff_thrsh); - println!( " {:<20} {:>10e}", "RMS FPS THRESH", calc_sett.commu_conv_thrsh); - println!(" {:<20} {:>10}", "Direct SCF", calc_sett.use_direct_scf); - println!(" {:<20} {:>10}", "Use DIIS", calc_sett.use_diis); - println!(" {:<20} {:>10}", "DIIS Settings", ""); - println!(" {:<20} {:>8}", "DIIS MIN", calc_sett.diis_sett.diis_min); - println!(" {:<20} {:>8}", "DIIS MAX", calc_sett.diis_sett.diis_max); - println!("{:-20}", ""); -} diff --git a/src/print_utils/print_scf.rs b/src/print_utils/print_scf.rs new file mode 100644 index 0000000..26c0c3b --- /dev/null +++ b/src/print_utils/print_scf.rs @@ -0,0 +1,28 @@ +use crate::calc_type::{CalcSettings, Reference}; + +pub(crate) fn print_scf_header_and_settings(calc_sett: &CalcSettings, calc_type: Reference) { + println!("{:=>35}", ""); + match calc_type { + Reference::RHF => println!("{:^35}", "RHF SCF"), + Reference::UHF => println!("{:^35}", "UHF SCF"), + Reference::ROHF => println!("{:^35}", "ROHF SCF"), + } + println!("{:=>35}", ""); + const SING_IND: &str = " "; + const DBL_IND: &str = " "; + + // ↓ this is used to not format this part + #[allow(clippy::deprecated_cfg_attr)] + #[cfg_attr(rustfmt, rustfmt_skip)] { + println!("{:-20}", ""); + println!("SCF settings:"); + println!("{}{:<20} {:>10e}", SING_IND, "ΔE THRESH", calc_sett.e_diff_thrsh); + println!("{}{:<20} {:>10e}", SING_IND, "RMS FPS THRESH", calc_sett.commu_conv_thrsh); + println!("{}{:<20} {:>10}", SING_IND, "Direct SCF", calc_sett.use_direct_scf); + println!("{}{:<20} {:>10}", SING_IND, "Use DIIS", calc_sett.use_diis); + println!("{}{:<20} {:>10}", SING_IND, "DIIS Settings", ""); + println!("{}{:<20} {:>8}", DBL_IND, "DIIS MIN", calc_sett.diis_sett.diis_min); + println!("{}{:<20} {:>8}", DBL_IND, "DIIS MAX", calc_sett.diis_sett.diis_max); + println!("{:-20}", ""); + } +}