From ac81cf8699fd5ed147ada0b91212a8bff6a79bf6 Mon Sep 17 00:00:00 2001 From: MRJD Date: Mon, 1 Jan 2024 23:27:45 +0100 Subject: [PATCH] Adding more printing information and options for SCF --- src/calc_type/mod.rs | 3 +-- src/calc_type/rhf.rs | 56 +++++++++++++++++++++++++++++++------------- src/main.rs | 33 ++++++++++++++++++++------ src/molecule/mod.rs | 20 ++++++++++++++-- 4 files changed, 85 insertions(+), 27 deletions(-) diff --git a/src/calc_type/mod.rs b/src/calc_type/mod.rs index 6c167e0..529ae1d 100644 --- a/src/calc_type/mod.rs +++ b/src/calc_type/mod.rs @@ -3,7 +3,7 @@ use crate::calc_type::rhf::calc_cmp_idx; use ndarray::{Array1, Array2}; use std::ops::{Index, IndexMut}; -mod rhf; +pub(crate) mod rhf; pub(crate) enum CalcType { RHF, @@ -40,7 +40,6 @@ pub struct SCF { C_matr_conv: Array2, P_matr_conv: Array2, // [ ] TODO: pot. change this to sparse matrix orb_energies_conv: Array1, - diis: DIIS, } diff --git a/src/calc_type/rhf.rs b/src/calc_type/rhf.rs index a0a78ed..aabd7a0 100644 --- a/src/calc_type/rhf.rs +++ b/src/calc_type/rhf.rs @@ -7,7 +7,7 @@ use crate::mol_int_and_deriv::{ use crate::molecule::Molecule; use ndarray::parallel::prelude::*; use ndarray::{s, Array, Array1, Array2, Zip}; -use ndarray_linalg::{Eigh, SymmetricSqrt, UPLO}; +use ndarray_linalg::{Eigh, UPLO}; use super::{CalcSettings, SCF}; @@ -190,6 +190,7 @@ pub(crate) fn rhf_scf_normal(calc_sett: CalcSettings, basis: &BasisSet, mol: &Mo // - [ ] Print header for RHF SCF // - [ ] Print settings + let mut is_conv = false; let mut scf = SCF::default(); let V_nuc = mol.calc_core_potential(); let (S_matr, H_core) = calc_1e_int_matrs(basis, mol); @@ -199,7 +200,7 @@ pub(crate) fn rhf_scf_normal(calc_sett: CalcSettings, basis: &BasisSet, mol: &Mo eri = calc_2e_int_matr(basis); } - let S_matr_inv_sqrt = inv_ssqrt(S_matr, UPLO::Upper); + let S_matr_inv_sqrt = inv_ssqrt(&S_matr, UPLO::Upper); // Init matrices for SCF loop let mut C_matr_AO; @@ -215,9 +216,9 @@ pub(crate) fn rhf_scf_normal(calc_sett: CalcSettings, basis: &BasisSet, mol: &Mo // Print SCF iteration Header println!( "{:>3} {:^20} {:^20} {:^20} {:^20} {:^20}", - "Iter", "E_scf", "E_tot", "RMS(P)", "Delta(E)", "FPS - SPF" + "Iter", "E_scf", "E_tot", "RMS(P)", "Delta(E)", "RMS(|FPS - SPF|)" ); - for scf_iter in 0..calc_sett.max_scf_iter { + 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); // println!("F_matr_pr:\n {}", &F_matr_pr); @@ -246,30 +247,51 @@ pub(crate) fn rhf_scf_normal(calc_sett: CalcSettings, basis: &BasisSet, mol: &Mo let rms_p_val = calc_rms_2_matr(&P_matr, &P_matr_old.clone()); let delta_E = E_scf_curr - E_scf_prev; + let fps_comm = calc_FPS_comm(&F_matr, &P_matr, &S_matr); + let rms_comm_val = + (fps_comm.par_iter().map(|x| x * x).sum::() / fps_comm.len() as f64).sqrt(); // SCF output println!( - "{:>3} {:>20.12} {:>20.12} {:>20.12} {:>20.12}", - scf_iter, E_scf_curr, scf.E_tot_conv, rms_p_val, delta_E + "{:>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 ); - // TODO! do not use rms_p_val for conv, but rather the commutator FPS - SPF =! 0 - if (delta_E.abs() < calc_sett.e_diff_thrsh) && (rms_p_val < calc_sett.commu_conv_thrsh) + if (delta_E.abs() < calc_sett.e_diff_thrsh) + && (rms_p_val < calc_sett.commu_conv_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 = C_matr_AO.clone(); scf.P_matr_conv = P_matr.clone(); scf.orb_energies_conv = orb_ener.clone(); - println!("SCF CONVERGED!"); - println!("Total SCF iterations: {}", scf.tot_scf_iter); + println!("\nSCF CONVERGED!"); + // println!("Total SCF iterations: {}", scf.tot_scf_iter); + is_conv = true; + break; + } else if scf_iter == calc_sett.max_scf_iter { + println!("\nSCF DID NOT CONVERGE!"); + // println!("Total SCF iterations: {}", scf.tot_scf_iter); break; } - E_scf_prev = E_scf_curr; P_matr_old = P_matr.clone(); calc_P_matr(&mut P_matr, &C_matr_AO, basis.no_occ()); } } + + if is_conv { + println!("{:*<55}", ""); + println!("* {:^51} *", "FINAL RESULTS"); + println!("{:*<55}", ""); + println!(" {:^50}", "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 } @@ -316,10 +338,10 @@ fn calc_rms_2_matr(matr1: &Array2, matr2: &Array2) -> f64 { .into_par_iter() .map(|(val1, val2)| (val1 - val2).powi(2)) .sum::() - .sqrt() + / (matr1.len() as f64).sqrt() } -fn inv_ssqrt(arr2: Array2, uplo: UPLO) -> Array2 { +fn inv_ssqrt(arr2: &Array2, uplo: UPLO) -> Array2 { let (e, v) = arr2.eigh(uplo).unwrap(); let e_inv_sqrt = Array1::from_iter(e.iter().map(|x| x.powf(-0.5))); let e_inv_sqrt_diag = Array::from_diag(&e_inv_sqrt); @@ -328,13 +350,14 @@ fn inv_ssqrt(arr2: Array2, uplo: UPLO) -> Array2 { } #[inline(always)] -fn calc_DIIS_error_matr( - F_pr_matr: &Array2, +/// This is also the DIIS error matrix +fn calc_FPS_comm( + F_matr: &Array2, P_matr: &Array2, S_matr: &Array2, // S_matr_inv_sqrt: &Array2, ) -> Array2 { - F_pr_matr.dot(P_matr).dot(S_matr) - S_matr.dot(P_matr).dot(F_pr_matr) + F_matr.dot(P_matr).dot(S_matr) - S_matr.dot(P_matr).dot(F_matr) } // TODO: fix this function @@ -433,5 +456,6 @@ mod tests { }; let _scf = rhf_scf_normal(calc_sett, &basis, &mol); + println!("{:?}", _scf); } } diff --git a/src/main.rs b/src/main.rs index 525d2d4..ad0d641 100644 --- a/src/main.rs +++ b/src/main.rs @@ -4,16 +4,17 @@ extern crate lazy_static; extern crate ndarray; extern crate openblas_src; - mod basisset; mod calc_type; mod mol_int_and_deriv; mod molecule; mod print_utils; -use crate::{print_utils::print_initial_header, calc_type::CalcType}; +use crate::{calc_type::CalcType, print_utils::print_initial_header}; use basisset::BasisSet; +use calc_type::{DiisSettings, CalcSettings}; use molecule::Molecule; +use crate::calc_type::rhf::rhf_scf_normal; fn main() { //################################## @@ -25,13 +26,13 @@ fn main() { print_initial_header(); exec_times.start("Molecule"); - let _mol = Molecule::new("data/xyz/water90.xyz", 0); - println!("Molecule: {:?}", _mol); + let mol = Molecule::new("data/xyz/water90.xyz", 0); + // println!("Molecule: {:?}", _mol); exec_times.stop("Molecule"); exec_times.start("BasisSet"); - let _basis = BasisSet::new("STO-3G", &_mol); - println!("BasisSet: {:?}", _basis); + let basis = BasisSet::new("STO-3G", &mol); + // println!("BasisSet: {:?}", _basis); // println!("Molecule: {:?}", _basis); // println!("\n\n"); // for shell in basis.shell_iter() { @@ -43,8 +44,26 @@ fn main() { //### BODY ### //################################## // Calculation type - let calc_type = CalcType::RHF; + // + exec_times.start("RHF noDIIS"); + let _calc_type = CalcType::RHF; + let calc_sett = CalcSettings { + max_scf_iter: 100, + e_diff_thrsh: 1e-10, + commu_conv_thrsh: 1e-10, + use_diis: false, + use_direct_scf: false, + diis_sett: DiisSettings { + diis_start: 0, + diis_end: 0, + diis_max: 0, + }, + }; + + let _scf = rhf_scf_normal(calc_sett, &basis, &mol); + exec_times.stop("RHF noDIIS"); + exec_times.stop("Total"); diff --git a/src/molecule/mod.rs b/src/molecule/mod.rs index cd32e96..d506532 100644 --- a/src/molecule/mod.rs +++ b/src/molecule/mod.rs @@ -7,6 +7,7 @@ use getset::CopyGetters; use ndarray::{s, Array1, Array2, Axis}; use ndarray_linalg::{Eigh, InverseH, UPLO}; use std::collections::HashMap; +use std::io::Seek; use std::str::FromStr; use std::{ fs::File, @@ -213,13 +214,28 @@ impl Molecule { no_atoms, } } + + fn print_input_file(reader: &mut BufReader) { + println!("{:*<30}", ""); + println!("* {:^26} *", "Inputfile:"); + println!("{:*<30}\n", ""); + + println!("{:-<75}", ""); + for line in reader.lines() { + println!("> {}", line.unwrap()); + } + println!("{:-<75}\n\n", ""); + } fn read_xyz_xmol_inputfile(geom_filename: &str) -> GeometryResult { println!("Inputfile: {geom_filename}"); println!("Reading geometry from input file...\n"); let geom_file = File::open(geom_filename)?; - let reader = BufReader::new(geom_file); + let mut reader = BufReader::new(geom_file); + Self::print_input_file(&mut reader); + reader.seek(std::io::SeekFrom::Start(0))?; // reset reader to start of file + let mut lines = reader .lines() .map(|line| line.expect("Failed to read line!")); @@ -238,7 +254,7 @@ impl Molecule { geom_matr[(at_idx, cc)] = line_parts.next().unwrap().parse().unwrap(); } } - + //* Convert geom_matr from Angstrom to Bohr (atomic units) */ const AA_TO_BOHR: f64 = 1.0e-10 / physical_constants::BOHR_RADIUS; geom_matr.mapv_inplace(|x| x * AA_TO_BOHR);