From 10dbe87ae8da3b5c9eefd3c0ab24d9097c9bdf65 Mon Sep 17 00:00:00 2001 From: MRJD Date: Fri, 5 Jan 2024 23:32:26 +0100 Subject: [PATCH] Starting with impl of SAD guess --- src/calc_type/guess.rs | 155 +++++++++++++++++++++++++++++++++++++++++ src/calc_type/mod.rs | 1 + 2 files changed, 156 insertions(+) create mode 100644 src/calc_type/guess.rs diff --git a/src/calc_type/guess.rs b/src/calc_type/guess.rs new file mode 100644 index 0000000..1282c01 --- /dev/null +++ b/src/calc_type/guess.rs @@ -0,0 +1,155 @@ +#![allow(clippy::upper_case_acronyms, non_snake_case, dead_code)] + +// Superposition of Atmoic Densities (SAD) guess implementation +// Source of implementation: Libint +// https://github.com/evaleev/libint/blob/4482d01abc37732a24b24799d3efe1d1f89a2891/include/libint2/chemistry/sto3g_atomic_density.h +// + +#[inline(always)] +pub const fn sad_calc_num_ao(z_val: u32) -> u32 { + match z_val { + 1 | 2 => 1, // H & He + 3..=10 => 5, // Li - Ne (2p incl. even for Li and Be) + 11..=18 => 9, // Na - Ar (3p incl. even for Na and Mg) + 19..=20 => 13, // K & Ca (4p incl. even for K and Ca) + 21..=36 => 18, // Sc - Kr + 37..=38 => 22, // Rb & Sr (5p incl. even for Rb and Sr) + 39..=53 => 27, // Y - I + _ => todo!(), + } +} + +// fn calc_occup_vec(z_val: u32) -> Vec { +// let num_ao = sad_calc_num_ao(z_val) as usize; +// let mut occ_vec = vec![0.0; num_ao]; +// +// let mut no_elec = z_val as usize; +// match z_val { +// 1 | 2 => calc_subsh_occup_vec(&mut occ_vec, &mut no_elec), +// 3..=10 => { +// calc_subsh_occup_vec(&mut occ_vec[0..1], &mut no_elec); +// calc_subsh_occup_vec(&mut occ_vec[1..2], &mut no_elec); +// calc_subsh_occup_vec(&mut occ_vec[2..], &mut no_elec); +// } +// 11..=17 => { +// calc_subsh_occup_vec(&mut occ_vec[0..1], &mut no_elec); +// calc_subsh_occup_vec(&mut occ_vec[1..2], &mut no_elec); +// calc_subsh_occup_vec(&mut occ_vec[2..=4], &mut no_elec); +// calc_subsh_occup_vec(&mut occ_vec[5..], &mut no_elec); +// } +// _ => todo!(), +// } +// +// occ_vec +// } + +fn calc_occup_vec(z_val: u32) -> Vec { + let num_ao = sad_calc_num_ao(z_val) as usize; + let mut occ_vec = vec![0.0; num_ao]; + + let mut no_elec = z_val as usize; + match z_val { + 1 | 2 => calc_subsh_occup_vec(&mut occ_vec, &mut no_elec), + 3..=10 => { + calc_subsh_occup_vec(&mut occ_vec[0..1], &mut no_elec); // 1s + calc_subsh_occup_vec(&mut occ_vec[1..2], &mut no_elec); // 2s + calc_subsh_occup_vec(&mut occ_vec[2..], &mut no_elec); // 2p + } + 11..=17 => { + calc_subsh_occup_vec(&mut occ_vec[0..1], &mut no_elec); // 1s + calc_subsh_occup_vec(&mut occ_vec[1..2], &mut no_elec); // 2s + calc_subsh_occup_vec(&mut occ_vec[2..=4], &mut no_elec); // 2p + calc_subsh_occup_vec(&mut occ_vec[5..6], &mut no_elec); // 3s + calc_subsh_occup_vec(&mut occ_vec[6..], &mut no_elec); // 3p + } + 18..=36 => { + calc_subsh_occup_vec(&mut occ_vec[0..1], &mut no_elec); // 1s + calc_subsh_occup_vec(&mut occ_vec[1..2], &mut no_elec); // 2s + calc_subsh_occup_vec(&mut occ_vec[2..=4], &mut no_elec); // 2p + calc_subsh_occup_vec(&mut occ_vec[5..6], &mut no_elec); // 3s + calc_subsh_occup_vec(&mut occ_vec[6..=8], &mut no_elec); // 3p + let mut no_4s_elec = if z_val == 19 || z_val == 24 || z_val == 29 { + 1 + } else { + 2 + }; + no_elec -= no_4s_elec; + calc_subsh_occup_vec(&mut occ_vec[9..10], &mut no_4s_elec); // 4s + + let mut no_4p_elec = 6.min(if z_val > 30 { z_val - 30 } else { 0 }) as usize; + calc_subsh_occup_vec(&mut occ_vec[10..=12], &mut no_4p_elec); // 4p + calc_subsh_occup_vec(&mut occ_vec[13..], &mut no_elec) // 3d + } + _ => todo!(), + } + + occ_vec +} + +fn calc_subsh_occup_vec(occvec: &mut [f64], no_elec: &mut usize) { + let size = occvec.len(); + let no_elec_alloc = (*no_elec).min(2 * size); + *no_elec -= no_elec_alloc; + let ne_per_orb = if no_elec_alloc % size == 0 { + (no_elec_alloc / size) as f64 + } else { + no_elec_alloc as f64 / size as f64 + }; + + for f in occvec.iter_mut() { + *f = ne_per_orb; + } +} + +#[cfg(test)] +mod tests { + use std::iter::zip; + + use super::*; + + #[test] + fn test_sad_calc_num_ao() { + assert_eq!(sad_calc_num_ao(1), 1); + assert_eq!(sad_calc_num_ao(2), 1); + assert_eq!(sad_calc_num_ao(3), 5); + assert_eq!(sad_calc_num_ao(10), 5); + assert_eq!(sad_calc_num_ao(11), 9); + assert_eq!(sad_calc_num_ao(18), 9); + assert_eq!(sad_calc_num_ao(19), 13); + assert_eq!(sad_calc_num_ao(20), 13); + assert_eq!(sad_calc_num_ao(21), 18); + assert_eq!(sad_calc_num_ao(36), 18); + assert_eq!(sad_calc_num_ao(37), 22); + assert_eq!(sad_calc_num_ao(38), 22); + assert_eq!(sad_calc_num_ao(39), 27); + assert_eq!(sad_calc_num_ao(53), 27); + } + + #[test] + fn test_calc_occup_vec() { + let z_val: u32 = 24; + let occup_vec = calc_occup_vec(z_val); + const ORB_NAMES: [&str; 18] = [ + "1s", "2s", "2px", "2py", "2pz", "3s", "3px", "3py", "3pz", "4s", "4px", "4py", "4pz", + "3dxz", "3dyz", "3dz2", "3dx2y2", "3dxy", + ]; + println!("Z_val: {}", z_val); + for (val, orb_name) in zip(occup_vec.iter(), ORB_NAMES.iter()) { + println!("{:>5}: {:>10.5}", orb_name, val); + } + // println!("{:?}", occup_vec); + } + // #[test] + // fn test_calc_subsh_occup_vec() { + // let mut test_vec = vec![1.2, 2.34, 234.1]; + // let size = 25; + // let mut no_elec = 10; + // println!("{:?}", no_elec); + // calc_subsh_occup_vec(&mut test_vec, size, &mut no_elec); + // println!("{:?}", no_elec); + // calc_subsh_occup_vec(&mut test_vec, size, &mut no_elec); + // println!("{:?}", no_elec); + // calc_subsh_occup_vec(&mut test_vec, size, &mut no_elec); + // println!("{:?}", no_elec); + // } +} diff --git a/src/calc_type/mod.rs b/src/calc_type/mod.rs index 8f7ce32..3c968a3 100644 --- a/src/calc_type/mod.rs +++ b/src/calc_type/mod.rs @@ -7,6 +7,7 @@ use std::ops::{Index, IndexMut}; pub(crate) mod rhf; pub(crate) mod uhf; +pub(crate) mod guess; pub(crate) enum CalcType { RHF,