Skip to content

Commit

Permalink
Adding more printing information and options for SCF
Browse files Browse the repository at this point in the history
  • Loading branch information
@MartinRJDagleish
MartinRJDagleish committed Jan 1, 2024
1 parent 3902519 commit ac81cf8
Show file tree
Hide file tree
Showing 4 changed files with 85 additions and 27 deletions.
3 changes: 1 addition & 2 deletions src/calc_type/mod.rs
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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,
Expand Down Expand Up @@ -40,7 +40,6 @@ pub struct SCF {
C_matr_conv: Array2<f64>,
P_matr_conv: Array2<f64>, // [ ] TODO: pot. change this to sparse matrix
orb_energies_conv: Array1<f64>,
diis: DIIS,
}


Expand Down
56 changes: 40 additions & 16 deletions src/calc_type/rhf.rs
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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};

Expand Down Expand Up @@ -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);
Expand All @@ -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;
Expand All @@ -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);
Expand Down Expand Up @@ -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::<f64>() / 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
}

Expand Down Expand Up @@ -316,10 +338,10 @@ fn calc_rms_2_matr(matr1: &Array2<f64>, matr2: &Array2<f64>) -> f64 {
.into_par_iter()
.map(|(val1, val2)| (val1 - val2).powi(2))
.sum::<f64>()
.sqrt()
/ (matr1.len() as f64).sqrt()
}

fn inv_ssqrt(arr2: Array2<f64>, uplo: UPLO) -> Array2<f64> {
fn inv_ssqrt(arr2: &Array2<f64>, uplo: UPLO) -> Array2<f64> {
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);
Expand All @@ -328,13 +350,14 @@ fn inv_ssqrt(arr2: Array2<f64>, uplo: UPLO) -> Array2<f64> {
}

#[inline(always)]
fn calc_DIIS_error_matr(
F_pr_matr: &Array2<f64>,
/// This is also the DIIS error matrix
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_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
Expand Down Expand Up @@ -433,5 +456,6 @@ mod tests {
};

let _scf = rhf_scf_normal(calc_sett, &basis, &mol);
println!("{:?}", _scf);
}
}
33 changes: 26 additions & 7 deletions src/main.rs
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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() {
//##################################
Expand All @@ -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() {
Expand All @@ -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");

Expand Down
20 changes: 18 additions & 2 deletions src/molecule/mod.rs
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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,
Expand Down Expand Up @@ -213,13 +214,28 @@ impl Molecule {
no_atoms,
}
}

fn print_input_file(reader: &mut BufReader<File>) {
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!"));
Expand All @@ -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);
Expand Down

0 comments on commit ac81cf8

Please sign in to comment.