From 859aa63db4c1575af27a7ff000c44e3cdd2a4703 Mon Sep 17 00:00:00 2001 From: MRJD Date: Sun, 14 Jan 2024 15:59:29 +0100 Subject: [PATCH] Midst rewriting UHF with struct --- src/calc_type/mod.rs | 119 ++++++- src/calc_type/rhf.rs | 269 ++-------------- src/calc_type/rhf_linscal.rs | 160 ++++++++++ src/calc_type/uhf.rs | 585 ++++++++++++++++++++++------------- src/main.rs | 9 +- 5 files changed, 665 insertions(+), 477 deletions(-) create mode 100644 src/calc_type/rhf_linscal.rs diff --git a/src/calc_type/mod.rs b/src/calc_type/mod.rs index bda2885..2cecc78 100644 --- a/src/calc_type/mod.rs +++ b/src/calc_type/mod.rs @@ -1,8 +1,10 @@ #![allow(non_snake_case)] use crate::basisset::BasisSet; use crate::calc_type::rhf::calc_cmp_idx; -use crate::mol_int_and_deriv::oe_int::{calc_overlap_int_cgto, calc_kinetic_int_cgto, calc_pot_int_cgto}; -use crate::mol_int_and_deriv::te_int::calc_ERI_int_cgto; +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::molecule::Molecule; use ndarray::parallel::prelude::*; use ndarray::{Array1, Array2, Zip}; @@ -12,6 +14,7 @@ use std::ops::{Index, IndexMut}; pub mod guess; pub mod rhf; pub mod uhf; +pub mod rhf_linscal; #[allow(non_camel_case_types)] pub(crate) enum HF_Ref { @@ -21,12 +24,6 @@ pub(crate) enum HF_Ref { } pub(crate) trait HF { - // fn new( - // basis: &BasisSet, - // calc_sett: &CalcSettings, - // ref_type: HF_Ref, - // ) -> Self; - fn run_scf( &mut self, calc_sett: &CalcSettings, @@ -35,6 +32,108 @@ pub(crate) trait HF { mol: &crate::molecule::Molecule, ) -> SCF; + /// ### Description + /// Calculate the 1e integrals for the given basis set and molecule. + /// Returns a tuple of the overlap, kinetic and potential energy matrices. + /// + /// ### Note + /// This is not the non-redudant version of the integrals, i.e. each function for the computation gets called separately. + /// + /// + /// ### Arguments + /// * `basis` - The basis set. + /// * `mol` - The molecule. + /// + fn calc_1e_int_matrs(basis: &BasisSet, mol: &Molecule) -> (Array2, Array2) { + let mut S_matr = Array2::::zeros((basis.no_bf(), basis.no_bf())); + let mut T_matr = Array2::::zeros((basis.no_bf(), basis.no_bf())); + let mut V_matr = Array2::::zeros((basis.no_bf(), basis.no_bf())); + + for (sh_idx1, shell1) in basis.shell_iter().enumerate() { + for sh_idx2 in 0..=sh_idx1 { + let shell2 = basis.shell(sh_idx2); + for (cgto_idx1, cgto1) in shell1.cgto_iter().enumerate() { + let mu = basis.sh_len_offset(sh_idx1) + cgto_idx1; + for (cgto_idx2, cgto2) in shell2.cgto_iter().enumerate() { + let nu = basis.sh_len_offset(sh_idx2) + cgto_idx2; + + // Overlap + S_matr[(mu, nu)] = if mu == nu { + 1.0 + } else { + calc_overlap_int_cgto(cgto1, cgto2) + }; + S_matr[(nu, mu)] = S_matr[(mu, nu)]; + + // Kinetic + T_matr[(mu, nu)] = calc_kinetic_int_cgto(cgto1, cgto2); + T_matr[(nu, mu)] = T_matr[(mu, nu)]; + + // Potential energy + V_matr[(mu, nu)] = calc_pot_int_cgto(cgto1, cgto2, mol); + V_matr[(nu, mu)] = V_matr[(mu, nu)]; + } + } + } + } + + // Return ovelap and core hamiltonian (T + V) + (S_matr, T_matr + V_matr) + } + + fn calc_2e_int_matr(basis: &BasisSet) -> EriArr1 { + let mut eri = EriArr1::new(basis.no_bf()); + + let no_shells = basis.no_shells(); + + for (sh_idx1, shell1) in basis.shell_iter().enumerate() { + for (cgto_idx1, cgto1) in shell1.cgto_iter().enumerate() { + let mu = basis.sh_len_offset(sh_idx1) + cgto_idx1; + + for sh_idx2 in 0..=sh_idx1 { + let shell2 = basis.shell(sh_idx2); + for (cgto_idx2, cgto2) in shell2.cgto_iter().enumerate() { + let nu = basis.sh_len_offset(sh_idx2) + cgto_idx2; + + if mu >= nu { + let munu = calc_cmp_idx(mu, nu); + + for sh_idx3 in 0..no_shells { + let shell3 = basis.shell(sh_idx3); + for (cgto_idx3, cgto3) in shell3.cgto_iter().enumerate() { + let lambda = basis.sh_len_offset(sh_idx3) + cgto_idx3; + + for sh_idx4 in 0..=sh_idx3 { + let shell4 = basis.shell(sh_idx4); + for (cgto_idx4, cgto4) in shell4.cgto_iter().enumerate() { + let sigma = basis.sh_len_offset(sh_idx4) + cgto_idx4; + + if lambda >= sigma { + let lambsig = calc_cmp_idx(lambda, sigma); + if munu >= lambsig { + let cmp_idx = calc_cmp_idx(munu, lambsig); + eri[cmp_idx] = calc_ERI_int_cgto( + cgto1, cgto2, cgto3, cgto4, + ); + // println!("{}: {}", cmp_idx, eri[cmp_idx]); + } else { + continue; + } + } + } + } + } + } + } else { + continue; + } + } + } + } + } + + eri + } } #[derive(Debug)] @@ -117,8 +216,8 @@ pub struct HFMatrices { #[derive(Debug, Default)] pub struct SCF { tot_scf_iter: usize, - pub E_tot_conv: f64, - pub E_scf_conv: f64, + E_tot_conv: f64, + E_scf_conv: f64, C_matr_conv_alph: Array2, P_matr_conv_alph: Array2, // [ ] TODO: pot. change this to sparse matrix C_matr_conv_beta: Option>, diff --git a/src/calc_type/rhf.rs b/src/calc_type/rhf.rs index ec0a703..a4df652 100644 --- a/src/calc_type/rhf.rs +++ b/src/calc_type/rhf.rs @@ -1,10 +1,10 @@ -use super::{CalcSettings, HFMatrices, SCF, HF}; +use super::{CalcSettings, HFMatrices, HF, SCF}; use crate::{ basisset::BasisSet, calc_type::{EriArr1, HF_Ref, 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, calc_schwarz_est_int_inp}, + te_int::{calc_ERI_int_cgto, calc_schwarz_est_int_inp}, }, molecule::Molecule, print_utils::{fmt_f64, print_scf::print_scf_header_and_settings, ExecTimes}, @@ -251,110 +251,6 @@ impl RHF { } } - /// ### Description - /// Calculate the 1e integrals for the given basis set and molecule. - /// Returns a tuple of the overlap, kinetic and potential energy matrices. - /// - /// ### Note - /// This is not the non-redudant version of the integrals, i.e. each function for the computation gets called separately. - /// - /// - /// ### Arguments - /// * `basis` - The basis set. - /// * `mol` - The molecule. - /// - pub fn calc_1e_int_matrs(basis: &BasisSet, mol: &Molecule) -> (Array2, Array2) { - let mut S_matr = Array2::::zeros((basis.no_bf(), basis.no_bf())); - let mut T_matr = Array2::::zeros((basis.no_bf(), basis.no_bf())); - let mut V_matr = Array2::::zeros((basis.no_bf(), basis.no_bf())); - - for (sh_idx1, shell1) in basis.shell_iter().enumerate() { - for sh_idx2 in 0..=sh_idx1 { - let shell2 = basis.shell(sh_idx2); - for (cgto_idx1, cgto1) in shell1.cgto_iter().enumerate() { - let mu = basis.sh_len_offset(sh_idx1) + cgto_idx1; - for (cgto_idx2, cgto2) in shell2.cgto_iter().enumerate() { - let nu = basis.sh_len_offset(sh_idx2) + cgto_idx2; - - // Overlap - S_matr[(mu, nu)] = if mu == nu { - 1.0 - } else { - calc_overlap_int_cgto(cgto1, cgto2) - }; - S_matr[(nu, mu)] = S_matr[(mu, nu)]; - - // Kinetic - T_matr[(mu, nu)] = calc_kinetic_int_cgto(cgto1, cgto2); - T_matr[(nu, mu)] = T_matr[(mu, nu)]; - - // Potential energy - V_matr[(mu, nu)] = calc_pot_int_cgto(cgto1, cgto2, mol); - V_matr[(nu, mu)] = V_matr[(mu, nu)]; - } - } - } - } - - // Return ovelap and core hamiltonian (T + V) - (S_matr, T_matr + V_matr) - } - - pub fn calc_2e_int_matr(basis: &BasisSet) -> EriArr1 { - let mut eri = EriArr1::new(basis.no_bf()); - - let no_shells = basis.no_shells(); - - for (sh_idx1, shell1) in basis.shell_iter().enumerate() { - for (cgto_idx1, cgto1) in shell1.cgto_iter().enumerate() { - let mu = basis.sh_len_offset(sh_idx1) + cgto_idx1; - - for sh_idx2 in 0..=sh_idx1 { - let shell2 = basis.shell(sh_idx2); - for (cgto_idx2, cgto2) in shell2.cgto_iter().enumerate() { - let nu = basis.sh_len_offset(sh_idx2) + cgto_idx2; - - if mu >= nu { - let munu = calc_cmp_idx(mu, nu); - - for sh_idx3 in 0..no_shells { - let shell3 = basis.shell(sh_idx3); - for (cgto_idx3, cgto3) in shell3.cgto_iter().enumerate() { - let lambda = basis.sh_len_offset(sh_idx3) + cgto_idx3; - - for sh_idx4 in 0..=sh_idx3 { - let shell4 = basis.shell(sh_idx4); - for (cgto_idx4, cgto4) in shell4.cgto_iter().enumerate() { - let sigma = basis.sh_len_offset(sh_idx4) + cgto_idx4; - - if lambda >= sigma { - let lambsig = calc_cmp_idx(lambda, sigma); - if munu >= lambsig { - let cmp_idx = calc_cmp_idx(munu, lambsig); - eri[cmp_idx] = calc_ERI_int_cgto( - cgto1, cgto2, cgto3, cgto4, - ); - // println!("{}: {}", cmp_idx, eri[cmp_idx]); - } else { - continue; - } - } - } - } - } - } - } else { - continue; - } - } - } - } - } - - eri - } - - fn calc_1e_int_matrs_inp(&mut self, basis: &BasisSet, mol: &Molecule) { println!("Calculating 1e integrals ..."); @@ -395,9 +291,10 @@ impl RHF { println!("FINSIHED calculating 1e integrals ..."); println!("Starting orthogonalization matrix calculation ..."); - self.hf_matrs.S_matr_inv_sqrt = Self::inv_ssqrt(&self.hf_matrs.S_matr, UPLO::Upper); + self.hf_matrs.S_matr_inv_sqrt = matr_inv_ssqrt(&self.hf_matrs.S_matr, UPLO::Upper); + println!("FINISHED orthogonalization matrix calculation ..."); } - + fn calc_2e_int_matr_inp(eri_arr: &mut EriArr1, basis: &BasisSet) { let no_shells = basis.no_shells(); @@ -474,8 +371,7 @@ impl RHF { } } - - fn calc_P_matr_rhf(P_matr: &mut Array2, C_matr: &Array2, no_occ: usize) { + pub(crate) 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) } @@ -503,7 +399,7 @@ impl RHF { } /// Calc Fock matrix in a direct SCF fashion (i.e. without precomputed eri tensor) - fn calc_new_F_matr_dir_scf_rhf( + pub(crate) fn calc_new_F_matr_dir_scf_rhf( F_matr: &mut Array2, delta_P_matr: &Array2, Schwarz_est_int: &Array2, @@ -585,7 +481,11 @@ impl RHF { .for_each(|(f, gg)| *f += 0.5 * *gg); } - fn calc_E_scf_rhf(P_matr: &Array2, H_core: &Array2, F_matr: &Array2) -> f64 { + pub(crate) fn calc_E_scf_rhf( + P_matr: &Array2, + H_core: &Array2, + F_matr: &Array2, + ) -> f64 { Zip::from(P_matr) .and(H_core) .and(F_matr) @@ -603,14 +503,6 @@ impl RHF { .sum::() / (matr1.len() as f64).sqrt() } - - pub(crate) 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); - let result = v.dot(&e_inv_sqrt_diag).dot(&v.t()); - result - } } #[inline(always)] @@ -621,138 +513,13 @@ pub fn calc_cmp_idx(idx1: usize, idx2: usize) -> usize { idx2 * (idx2 + 1) / 2 + idx1 } } -/// My try of a implementation of the density matrix based approach to solve the RHF SCF equations. -/// This is NOT actually linear scaling, since I am not using any tricks to reduce the computational cost. -/// No sparse matrices or anything like that. -/// -/// Missing: CFMM for Coulomb and LinK for Exchange matrix -#[allow(unused)] -fn rhf_scf_linscal( - calc_sett: &CalcSettings, - exec_times: &mut ExecTimes, - basis: &BasisSet, - mol: &Molecule, -) { - print_scf_header_and_settings(calc_sett, HF_Ref::RHF_ref); - const SHOW_ALL_CONV_CRIT: bool = false; - - let mut is_scf_conv = false; - let mut scf = SCF::default(); - - 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) = RHF::calc_1e_int_matrs(basis, mol); - println!("FINSIHED calculating 1e integrals ..."); - - let mut eri; - let schwarz_est_matr; - if calc_sett.use_direct_scf { - eri = 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 = Some(RHF::calc_2e_int_matr(basis)); - println!("FINSIHED calculating 2e integrals ..."); - } - - // Init matrices for SCF loop - 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 = P_matr.clone(); - let mut diis_str = ""; - - // Initial guess -> H_core - let mut F_matr = H_core.clone(); - - let S_matr_inv = S_matr.invh().unwrap(); - - // Print SCF iteration Header - 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 { - let S_matr_inv_sqrt = RHF::inv_ssqrt(&S_matr, UPLO::Upper); - let F_matr_pr = S_matr_inv_sqrt.dot(&F_matr).dot(&S_matr_inv_sqrt); - - let (orb_ener, C_matr_MO) = F_matr_pr.eigh(UPLO::Upper).unwrap(); - let C_matr_AO = S_matr_inv_sqrt.dot(&C_matr_MO); - - RHF::calc_P_matr_rhf(&mut P_matr, &C_matr_AO, basis.no_occ()); - println!("Guess P_matr:\n {:12.6}", &P_matr); - println!( - "Eigenvalues of guess: {:12.6}", - &P_matr.eigh(UPLO::Upper).unwrap().0 - ); - delta_P_matr = P_matr.clone(); - } else { - // 1. First new Fock matrix - RHF::calc_new_F_matr_dir_scf_rhf( - &mut F_matr, - &delta_P_matr, - schwarz_est_matr.as_ref().unwrap(), - basis, - ); - // 2. Calc E_scf - let E_scf_curr = RHF::calc_E_scf_rhf(&P_matr, &H_core, &F_matr); - scf.E_tot_conv = E_scf_curr + V_nuc; - - let delta_E = E_scf_curr - E_scf_prev; - println!( - "{:>3} {:>20.12} {:>20.12} {}", - scf_iter, - E_scf_curr, - scf.E_tot_conv, - fmt_f64(delta_E, 20, 8, 2), - ); - - // 3. Save old density matrix - P_matr_old = P_matr.clone(); - // 4. Calculate new density matrix - P_matr = calc_new_P_matr_linscal_contravar(&P_matr, &F_matr, &S_matr, &S_matr_inv); - // 5. Calculate calculate ΔP matrix - delta_P_matr = (&P_matr - &P_matr_old).to_owned(); - E_scf_prev = E_scf_curr; - } - } -} - -#[allow(unused)] -/// Result is contravariant -fn calc_new_P_matr_linscal_contravar( - P_matr_curr: &Array2, - F_matr: &Array2, - S_matr: &Array2, - S_matr_inv: &Array2, -) -> Array2 { - const STEP_WIDTH: f64 = 0.000_001; - - let energy_grad = calc_energy_gradient_p_matr(F_matr, P_matr_curr, S_matr); - P_matr_curr - STEP_WIDTH * S_matr_inv.dot(&energy_grad).dot(S_matr_inv) -} -fn calc_energy_gradient_p_matr( - F_matr: &Array2, - P_matr: &Array2, - S_matr: &Array2, -) -> Array2 { - 3.0 * F_matr.dot(P_matr).dot(S_matr) + 3.0 * S_matr.dot(P_matr).dot(F_matr) - - 2.0 * S_matr.dot(P_matr).dot(F_matr).dot(P_matr).dot(S_matr) - - 2.0 * F_matr.dot(P_matr).dot(S_matr).dot(P_matr).dot(S_matr) - - 2.0 * S_matr.dot(P_matr).dot(S_matr).dot(P_matr).dot(F_matr) +pub(crate) fn matr_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); + let result = v.dot(&e_inv_sqrt_diag).dot(&v.t()); + result } #[cfg(test)] diff --git a/src/calc_type/rhf_linscal.rs b/src/calc_type/rhf_linscal.rs new file mode 100644 index 0000000..b62b639 --- /dev/null +++ b/src/calc_type/rhf_linscal.rs @@ -0,0 +1,160 @@ +use ndarray::Array2; +use ndarray_linalg::{InverseH, Eigh, UPLO}; + +use crate::{ + basisset::BasisSet, + calc_type::{HF_Ref, SCF, rhf::{RHF, matr_inv_ssqrt}, HF}, + molecule::Molecule, + print_utils::{print_scf::print_scf_header_and_settings, ExecTimes, fmt_f64}, mol_int_and_deriv::te_int::calc_schwarz_est_int, +}; + +use super::CalcSettings; + +/// My try of a implementation of the density matrix based approach to solve the RHF SCF equations. +/// This is NOT actually linear scaling, since I am not using any tricks to reduce the computational cost. +/// No sparse matrices or anything like that. +/// +/// Missing: CFMM for Coulomb and LinK for Exchange matrix +#[allow(unused)] +fn rhf_scf_linscal( + calc_sett: &CalcSettings, + exec_times: &mut ExecTimes, + basis: &BasisSet, + mol: &Molecule, +) { + print_scf_header_and_settings(calc_sett, HF_Ref::RHF_ref); + const SHOW_ALL_CONV_CRIT: bool = false; + + let mut is_scf_conv = false; + let mut scf = SCF::default(); + + 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) = RHF::calc_1e_int_matrs(basis, mol); + println!("FINSIHED calculating 1e integrals ..."); + + let mut eri; + let schwarz_est_matr; + if calc_sett.use_direct_scf { + eri = 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 = Some(RHF::calc_2e_int_matr(basis)); + println!("FINSIHED calculating 2e integrals ..."); + } + + // Init matrices for SCF loop + 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 = P_matr.clone(); + let mut diis_str = ""; + + // Initial guess -> H_core + let mut F_matr = H_core.clone(); + + let S_matr_inv = S_matr.invh().unwrap(); + + // Print SCF iteration Header + 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 { + let S_matr_inv_sqrt = matr_inv_ssqrt(&S_matr, UPLO::Upper); + let F_matr_pr = S_matr_inv_sqrt.dot(&F_matr).dot(&S_matr_inv_sqrt); + + let (orb_ener, C_matr_MO) = F_matr_pr.eigh(UPLO::Upper).unwrap(); + let C_matr_AO = S_matr_inv_sqrt.dot(&C_matr_MO); + + RHF::calc_P_matr_rhf(&mut P_matr, &C_matr_AO, basis.no_occ()); + println!("Guess P_matr:\n {:12.6}", &P_matr); + println!( + "Eigenvalues of guess: {:12.6}", + &P_matr.eigh(UPLO::Upper).unwrap().0 + ); + delta_P_matr = P_matr.clone(); + } else { + // 1. First new Fock matrix + RHF::calc_new_F_matr_dir_scf_rhf( + &mut F_matr, + &delta_P_matr, + schwarz_est_matr.as_ref().unwrap(), + basis, + ); + // 2. Calc E_scf + let E_scf_curr = RHF::calc_E_scf_rhf(&P_matr, &H_core, &F_matr); + scf.E_tot_conv = E_scf_curr + V_nuc; + + let delta_E = E_scf_curr - E_scf_prev; + println!( + "{:>3} {:>20.12} {:>20.12} {}", + scf_iter, + E_scf_curr, + scf.E_tot_conv, + fmt_f64(delta_E, 20, 8, 2), + ); + + // 3. Save old density matrix + P_matr_old = P_matr.clone(); + // 4. Calculate new density matrix + P_matr = calc_new_P_matr_linscal_contravar(&P_matr, &F_matr, &S_matr, &S_matr_inv); + // 5. Calculate calculate ΔP matrix + delta_P_matr = (&P_matr - &P_matr_old).to_owned(); + E_scf_prev = E_scf_curr; + } + } +} + +#[allow(unused)] +/// Result is contravariant +fn calc_new_P_matr_linscal_contravar( + P_matr_curr: &Array2, + F_matr: &Array2, + S_matr: &Array2, + S_matr_inv: &Array2, +) -> Array2 { + const STEP_WIDTH: f64 = 0.000_001; + + let energy_grad = calc_energy_gradient_p_matr(F_matr, P_matr_curr, S_matr); + P_matr_curr - STEP_WIDTH * S_matr_inv.dot(&energy_grad).dot(S_matr_inv) +} + +fn calc_energy_gradient_p_matr( + F_matr: &Array2, + P_matr: &Array2, + S_matr: &Array2, +) -> Array2 { + 3.0 * F_matr.dot(P_matr).dot(S_matr) + 3.0 * S_matr.dot(P_matr).dot(F_matr) + - 2.0 * S_matr.dot(P_matr).dot(F_matr).dot(P_matr).dot(S_matr) + - 2.0 * F_matr.dot(P_matr).dot(S_matr).dot(P_matr).dot(S_matr) + - 2.0 * S_matr.dot(P_matr).dot(S_matr).dot(P_matr).dot(F_matr) +} + +#[cfg(test)] +mod tests { + use super::*; + + #[test] + fn test_calc_energy_gradient_p_matr() { + let F_matr = Array2::::zeros((3, 3)); + let P_matr = Array2::::zeros((3, 3)); + let S_matr = Array2::::zeros((3, 3)); + + let grad = calc_energy_gradient_p_matr(&F_matr, &P_matr, &S_matr); + assert_eq!(grad, Array2::::zeros((3, 3))); + } +} diff --git a/src/calc_type/uhf.rs b/src/calc_type/uhf.rs index 6f7a76e..85d15e0 100644 --- a/src/calc_type/uhf.rs +++ b/src/calc_type/uhf.rs @@ -1,11 +1,14 @@ -use super::{CalcSettings, EriArr1, SCF, HFMatrices, HF_Ref}; +use super::{rhf::calc_cmp_idx, CalcSettings, EriArr1, HFMatrices, HF_Ref, SCF}; use crate::{ basisset::BasisSet, calc_type::{ - rhf::RHF, + rhf::{matr_inv_ssqrt, RHF}, DIIS, HF, }, - mol_int_and_deriv::te_int::calc_schwarz_est_int, + mol_int_and_deriv::{ + oe_int::{calc_kinetic_int_cgto, calc_overlap_int_cgto, calc_pot_int_cgto}, + te_int::{calc_schwarz_est_int, calc_schwarz_est_int_inp}, + }, molecule::Molecule, print_utils::{fmt_f64, print_scf::print_scf_header_and_settings}, }; @@ -81,6 +84,373 @@ impl UHF { } } + fn calc_1e_int_matrs_inp(&mut self, basis: &BasisSet, mol: &Molecule) { + println!("Calculating 1e integrals ..."); + + for (sh_idx1, shell1) in basis.shell_iter().enumerate() { + for sh_idx2 in 0..=sh_idx1 { + let shell2 = basis.shell(sh_idx2); + for (cgto_idx1, cgto1) in shell1.cgto_iter().enumerate() { + let mu = basis.sh_len_offset(sh_idx1) + cgto_idx1; + for (cgto_idx2, cgto2) in shell2.cgto_iter().enumerate() { + let nu = basis.sh_len_offset(sh_idx2) + cgto_idx2; + + // Overlap + self.hf_matrs.S_matr[(mu, nu)] = if mu == nu { + 1.0 + } else { + calc_overlap_int_cgto(cgto1, cgto2) + }; + self.hf_matrs.S_matr[(nu, mu)] = self.hf_matrs.S_matr[(mu, nu)]; + + // Kinetic + self.hf_matrs.T_matr[(mu, nu)] = calc_kinetic_int_cgto(cgto1, cgto2); + self.hf_matrs.T_matr[(nu, mu)] = self.hf_matrs.T_matr[(mu, nu)]; + + // Potential energy + self.hf_matrs.V_ne_matr[(mu, nu)] = calc_pot_int_cgto(cgto1, cgto2, mol); + self.hf_matrs.V_ne_matr[(nu, mu)] = self.hf_matrs.V_ne_matr[(mu, nu)]; + } + } + } + } + println!("S_matr: \n{:>12.8}", &self.hf_matrs.S_matr); + println!("T_matr: \n{:>12.8}", &self.hf_matrs.T_matr); + println!("V_matr: \n{:>12.8}", &self.hf_matrs.V_ne_matr); + + Zip::from(self.hf_matrs.T_matr.view()) + .and(self.hf_matrs.V_ne_matr.view()) + .par_map_assign_into(&mut self.hf_matrs.H_core_matr, |&t, &v| t + v); + println!("FINSIHED calculating 1e integrals ..."); + + println!("Starting orthogonalization matrix calculation ..."); + self.hf_matrs.S_matr_inv_sqrt = matr_inv_ssqrt(&self.hf_matrs.S_matr, UPLO::Upper); + println!("FINISHED orthogonalization matrix calculation ..."); + } + + fn calc_2e_int_matr_inp(eri_arr: &mut EriArr1, basis: &BasisSet) { + let no_shells = basis.no_shells(); + + for (sh_idx1, shell1) in basis.shell_iter().enumerate() { + for (cgto_idx1, cgto1) in shell1.cgto_iter().enumerate() { + let mu = basis.sh_len_offset(sh_idx1) + cgto_idx1; + + for sh_idx2 in 0..=sh_idx1 { + let shell2 = basis.shell(sh_idx2); + for (cgto_idx2, cgto2) in shell2.cgto_iter().enumerate() { + let nu = basis.sh_len_offset(sh_idx2) + cgto_idx2; + + if mu >= nu { + let munu = calc_cmp_idx(mu, nu); + + for sh_idx3 in 0..no_shells { + let shell3 = basis.shell(sh_idx3); + for (cgto_idx3, cgto3) in shell3.cgto_iter().enumerate() { + let lambda = basis.sh_len_offset(sh_idx3) + cgto_idx3; + + for sh_idx4 in 0..=sh_idx3 { + let shell4 = basis.shell(sh_idx4); + for (cgto_idx4, cgto4) in shell4.cgto_iter().enumerate() { + let sigma = basis.sh_len_offset(sh_idx4) + cgto_idx4; + + if lambda >= sigma { + let lambsig = calc_cmp_idx(lambda, sigma); + if munu >= lambsig { + let cmp_idx = calc_cmp_idx(munu, lambsig); + eri_arr[cmp_idx] = calc_ERI_int_cgto( + cgto1, cgto2, cgto3, cgto4, + ); + // println!("{}: {}", cmp_idx, eri[cmp_idx]); + } else { + continue; + } + } + } + } + } + } + } else { + continue; + } + } + } + } + } + } + + fn dir_indir_scf_2e_matr(&mut self, basis: &BasisSet, calc_sett: &CalcSettings) { + match calc_sett.use_direct_scf { + true => { + println!("Calculating Schwarz int estimates ..."); + calc_schwarz_est_int_inp(self.hf_matrs.schwarz_est.as_mut().unwrap(), basis); + println!("FINISHED Schwarz int estimates ..."); + } + false => { + println!("Calculating 2e integrals ..."); + Self::calc_2e_int_matr_inp(self.hf_matrs.eri_opt.as_mut().unwrap(), basis); + // println!("i j k l; ijkl: val"); + // for i in 0..basis.no_bf() { + // for j in 0..basis.no_bf() { + // for k in 0..basis.no_bf() { + // for l in 0..basis.no_bf() { + // let ijkl = calc_cmp_idx(calc_cmp_idx(i, j), calc_cmp_idx(k, l)); + // println!("{:>3}{:>3}{:>3}{:>3}; {:>4}: {:>12.8}", i,j,k,l,ijkl, self.hf_matrs.eri_opt.as_ref().unwrap()[ijkl]); + // } + // } + // } + // } + println!("FINSIHED calculating 2e integrals ..."); + } + } + } + + pub(crate) fn uhf_scf_normal( + &mut self, + calc_sett: &CalcSettings, + exec_times: &mut crate::print_utils::ExecTimes, + basis: &BasisSet, + mol: &Molecule, + ) -> SCF { + print_scf_header_and_settings(calc_sett, crate::calc_type::HF_Ref::UHF_ref); + + let mut is_scf_conv = false; + let mut scf = SCF::default(); + + // TODO: [ ] account for Multiplicity != 1 -> differnt input handling + let no_elec_half = mol.no_elec() / 2; + let (no_alpha, no_beta) = if mol.no_elec() % 2 == 0 { + (no_elec_half, no_elec_half) + } else { + (no_elec_half + 1, no_elec_half) + }; + + let mut diis_alph = if calc_sett.use_diis { + Some(DIIS::new( + &calc_sett.diis_sett, + [basis.no_bf(), basis.no_bf()], + )) + } else { + None + }; + let mut diis_beta = diis_alph.clone(); + + let V_nuc: f64 = if mol.no_atoms() > 100 { + mol.calc_core_potential_par() + } else { + mol.calc_core_potential_ser() + }; + + // Calculate 1e ints + exec_times.start("1e ints"); + self.calc_1e_int_matrs_inp(basis, mol); + exec_times.stop("1e ints"); + + // Calculate 2e ints / Schwarz estimates + exec_times.start("2e ints / Schwarz esti."); + self.dir_indir_scf_2e_matr(basis, calc_sett); + exec_times.stop("2e ints / Schwarz esti."); + + // Init matrices for SCF loop + // let (mut orb_E_alph, mut C_matr_MO_alph) = init_diag_F_matr(&H_core, &S_matr_inv_sqrt); + // let mut C_matr_AO_alph = S_matr_inv_sqrt.dot(&C_matr_MO_alph); + // + // let mut C_matr_MO_beta = C_matr_MO_alph.clone(); + // let mut C_matr_AO_beta = C_matr_AO_alph.clone(); + // let mut orb_E_beta = orb_E_alph.clone(); + // + // let mut E_scf_prev = 0.0; + // + // let mut P_matr_alph = Array2::::zeros((basis.no_bf(), basis.no_bf())); + // build_P_matr_uhf(&mut P_matr_alph, &C_matr_AO_alph, no_alpha); + // let mut P_matr_beta = Array2::::zeros((basis.no_bf(), basis.no_bf())); + // build_P_matr_uhf(&mut P_matr_beta, &C_matr_AO_beta, no_beta); + // let mut P_matr_old_alph = P_matr_alph.clone(); + // let mut P_matr_old_beta = P_matr_beta.clone(); + // + // // Initial guess -> H_core + // let mut F_matr_alph = Array2::::zeros((basis.no_bf(), basis.no_bf())); + // let mut F_matr_beta = Array2::::zeros((basis.no_bf(), basis.no_bf())); + // + // let mut F_matr_pr_alph; + // let mut F_matr_pr_beta; + + println!( + "{:>3} {:^20} {:^20} {:^20} {:^20}", + "Iter", "E_scf", "E_tot", "ΔE", "RMS(|FPS - SPF|)" + ); + + let mut diis_str = ""; + for scf_iter in 0..=calc_sett.max_scf_iter { + if scf_iter == 0 { + self.hf_matrs.F_matr_alpha = self.hf_matrs.H_core_matr.clone(); + self.hf_matrs.F_matr_beta = Some(self.hf_matrs.H_core_matr.clone()); + + self.hf_matrs.F_matr_pr_alpha = self + .hf_matrs + .S_matr_inv_sqrt + .dot(&self.hf_matrs.F_matr_alpha) + .dot(&self.hf_matrs.S_matr_inv_sqrt); + self.hf_matrs.F_matr_pr_beta = Some(self + .hf_matrs + .S_matr_inv_sqrt + .dot(self.hf_matrs.F_matr_beta.as_ref().unwrap()) + .dot(&self.hf_matrs.S_matr_inv_sqrt)); + + (self.hf_matrs.orb_ener_alpha, self.hf_matrs.C_matr_MO_alpha) = + self.hf_matrs.F_matr_pr_alpha.eigh(UPLO::Upper).unwrap(); + self.hf_matrs.C_matr_AO_alpha = self + .hf_matrs + .S_matr_inv_sqrt + .dot(&self.hf_matrs.C_matr_MO_alpha); + + (self.hf_matrs.orb_ener_beta, self.hf_matrs.C_matr_MO_beta) = + self.hf_matrs.F_matr_pr_beta.eigh(UPLO::Upper).unwrap(); + self.hf_matrs.C_matr_AO_beta = self + .hf_matrs + .S_matr_inv_sqrt + .dot(&self.hf_matrs.C_matr_MO_beta); + + Self::calc_P_matr_uhf( + &mut self.hf_matrs.P_matr_alpha, + &self.hf_matrs.C_matr_AO_alpha, + no_alpha, + ); + Self::calc_P_matr_uhf( + &mut self.hf_matrs.P_matr_beta, + &self.hf_matrs.C_matr_AO_beta, + no_beta, + ); + if calc_sett.use_direct_scf { + self.hf_matrs.delta_P_matr_alpha = Some(self.hf_matrs.P_matr_alpha.clone()); + self.hf_matrs.delta_P_matr_beta = Some(self.hf_matrs.P_matr_beta.clone()); + } + } else { + /// Direct or Indirect SCF + match eri_opt { + // Indirect SCF + Some(ref eri) => { + calc_new_F_matr_ind_scf_uhf( + &mut F_matr_alph, + &H_core, + &P_matr_alph, + &P_matr_beta, + eri, + true, + ); + calc_new_F_matr_ind_scf_uhf( + &mut F_matr_beta, + &H_core, + &P_matr_alph, + &P_matr_beta, + eri, + false, + ); + } + // Direct SCF + None => { + todo!("Direct SCF for UHF not yet implemented!") + } + } + let E_scf_curr = calc_E_scf_uhf( + &P_matr_alph, + &P_matr_beta, + &H_core, + &F_matr_alph, + &F_matr_beta, + ); + scf.E_tot_conv = E_scf_curr + V_nuc; + let fps_comm_alph = DIIS::calc_FPS_comm(&F_matr_alph, &P_matr_alph, &S_matr); + let fps_comm_beta = DIIS::calc_FPS_comm(&F_matr_beta, &P_matr_beta, &S_matr); + + F_matr_pr_alph = S_matr_inv_sqrt.dot(&F_matr_alph).dot(&S_matr_inv_sqrt); + F_matr_pr_beta = S_matr_inv_sqrt.dot(&F_matr_beta).dot(&S_matr_inv_sqrt); + + if calc_sett.use_diis { + let replace_idx = (scf_iter - 1) % calc_sett.diis_sett.diis_max; // always start with 0 + let err_matr_alph = S_matr_inv_sqrt.dot(&fps_comm_alph).dot(&S_matr_inv_sqrt); + let err_matr_beta = S_matr_inv_sqrt.dot(&fps_comm_beta).dot(&S_matr_inv_sqrt); + diis_alph.as_mut().unwrap().push_to_ring_buf( + &F_matr_pr_alph, + &err_matr_alph, + replace_idx, + ); + diis_beta.as_mut().unwrap().push_to_ring_buf( + &F_matr_pr_beta, + &err_matr_beta, + replace_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_alph = diis_alph.as_ref().unwrap().run_DIIS(err_set_len); + F_matr_pr_beta = diis_beta.as_ref().unwrap().run_DIIS(err_set_len); + diis_str = "DIIS"; + } + } + + (orb_E_alph, C_matr_MO_alph) = F_matr_pr_alph.eigh(UPLO::Upper).unwrap(); + C_matr_AO_alph = S_matr_inv_sqrt.dot(&C_matr_MO_alph); + + (orb_E_beta, C_matr_MO_beta) = F_matr_pr_beta.eigh(UPLO::Upper).unwrap(); + C_matr_AO_beta = S_matr_inv_sqrt.dot(&C_matr_MO_beta); + + let delta_E = E_scf_curr - E_scf_prev; + let rms_comm_val = calc_rms_comm_val_uhf(&fps_comm_alph, &fps_comm_beta); + + 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 scf_conv_check(calc_sett, delta_E, rms_comm_val) { + scf.tot_scf_iter = scf_iter; + scf.E_scf_conv = E_scf_curr; + scf.C_matr_conv_alph = C_matr_AO_alph; + scf.C_matr_conv_beta = Some(C_matr_AO_beta); + scf.P_matr_conv_alph = P_matr_alph; + scf.P_matr_conv_beta = Some(P_matr_beta); + scf.orb_E_conv_alph = orb_E_alph; + scf.orb_E_conv_beta = Some(orb_E_beta); + 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_alph = P_matr_alph.clone(); + P_matr_old_beta = P_matr_beta.clone(); + build_P_matr_uhf(&mut P_matr_alph, &C_matr_AO_alph, no_alpha); + build_P_matr_uhf(&mut P_matr_beta, &C_matr_AO_beta, no_beta); + if calc_sett.use_direct_scf { + todo!() + // 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}", "UHF 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 + } + // pub fn run_scf( // &mut self, // calc_sett: &CalcSettings, @@ -285,214 +655,7 @@ impl UHF { // } } -#[allow(unused)] -pub(crate) fn uhf_scf_normal( - calc_sett: &CalcSettings, - exec_times: &mut crate::print_utils::ExecTimes, - basis: &BasisSet, - mol: &Molecule, -) -> SCF { - print_scf_header_and_settings(calc_sett, crate::calc_type::HF_Ref::UHF_ref); - let mut is_scf_conv = false; - let mut scf = SCF::default(); - - // TODO: [ ] account for Multiplicity != 1 -> differnt input handling - let no_elec_half = mol.no_elec() / 2; - let (no_alpha, no_beta) = if mol.no_elec() % 2 == 0 { - (no_elec_half, no_elec_half) - } else { - (no_elec_half + 1, no_elec_half) - }; - - let mut diis_alph = if calc_sett.use_diis { - Some(DIIS::new(&calc_sett.diis_sett, [basis.no_bf(), basis.no_bf()])) - } else { - None - }; - let mut diis_beta = diis_alph.clone(); - - 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) = RHF::calc_1e_int_matrs(basis, mol); - println!("FINSIHED calculating 1e integrals ..."); - - let 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 ...\n"); - } else { - schwarz_est_matr = None; - - println!("Calculating 2e integrals ..."); - eri_opt = Some(RHF::calc_2e_int_matr(basis)); - println!("FINSIHED calculating 2e integrals ...\n"); - } - - let S_matr_inv_sqrt = RHF::inv_ssqrt(&S_matr, UPLO::Upper); - - // Init matrices for SCF loop - let (mut orb_E_alph, mut C_matr_MO_alph) = init_diag_F_matr(&H_core, &S_matr_inv_sqrt); - let mut C_matr_AO_alph = S_matr_inv_sqrt.dot(&C_matr_MO_alph); - - let mut C_matr_MO_beta = C_matr_MO_alph.clone(); - let mut C_matr_AO_beta = C_matr_AO_alph.clone(); - let mut orb_E_beta = orb_E_alph.clone(); - - let mut E_scf_prev = 0.0; - - let mut P_matr_alph = Array2::::zeros((basis.no_bf(), basis.no_bf())); - build_P_matr_uhf(&mut P_matr_alph, &C_matr_AO_alph, no_alpha); - let mut P_matr_beta = Array2::::zeros((basis.no_bf(), basis.no_bf())); - build_P_matr_uhf(&mut P_matr_beta, &C_matr_AO_beta, no_beta); - let mut P_matr_old_alph = P_matr_alph.clone(); - let mut P_matr_old_beta = P_matr_beta.clone(); - - - // Initial guess -> H_core - let mut F_matr_alph = Array2::::zeros((basis.no_bf(), basis.no_bf())); - let mut F_matr_beta = Array2::::zeros((basis.no_bf(), basis.no_bf())); - - let mut F_matr_pr_alph; - let mut F_matr_pr_beta; - - println!( - "{:>3} {:^20} {:^20} {:^20} {:^20}", - "Iter", "E_scf", "E_tot", "ΔE", "RMS(|FPS - SPF|)" - ); - - let mut diis_str = ""; - for scf_iter in 1..=calc_sett.max_scf_iter { - /// Direct or Indirect SCF - match eri_opt { - // Indirect SCF - Some(ref eri) => { - calc_new_F_matr_ind_scf_uhf( - &mut F_matr_alph, - &H_core, - &P_matr_alph, - &P_matr_beta, - eri, - true, - ); - calc_new_F_matr_ind_scf_uhf( - &mut F_matr_beta, - &H_core, - &P_matr_alph, - &P_matr_beta, - eri, - false, - ); - } - // Direct SCF - None => { - todo!("Direct SCF for UHF not yet implemented!") - } - } - let E_scf_curr = calc_E_scf_uhf( - &P_matr_alph, - &P_matr_beta, - &H_core, - &F_matr_alph, - &F_matr_beta, - ); - scf.E_tot_conv = E_scf_curr + V_nuc; - let fps_comm_alph = DIIS::calc_FPS_comm(&F_matr_alph, &P_matr_alph, &S_matr); - let fps_comm_beta = DIIS::calc_FPS_comm(&F_matr_beta, &P_matr_beta, &S_matr); - - F_matr_pr_alph = S_matr_inv_sqrt.dot(&F_matr_alph).dot(&S_matr_inv_sqrt); - F_matr_pr_beta = S_matr_inv_sqrt.dot(&F_matr_beta).dot(&S_matr_inv_sqrt); - - if calc_sett.use_diis { - let replace_idx = (scf_iter - 1) % calc_sett.diis_sett.diis_max; // always start with 0 - let err_matr_alph = S_matr_inv_sqrt.dot(&fps_comm_alph).dot(&S_matr_inv_sqrt); - let err_matr_beta = S_matr_inv_sqrt.dot(&fps_comm_beta).dot(&S_matr_inv_sqrt); - diis_alph - .as_mut() - .unwrap() - .push_to_ring_buf(&F_matr_pr_alph, &err_matr_alph, replace_idx); - diis_beta - .as_mut() - .unwrap() - .push_to_ring_buf(&F_matr_pr_beta, &err_matr_beta, replace_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_alph = diis_alph.as_ref().unwrap().run_DIIS(err_set_len); - F_matr_pr_beta = diis_beta.as_ref().unwrap().run_DIIS(err_set_len); - diis_str = "DIIS"; - } - } - - (orb_E_alph, C_matr_MO_alph) = F_matr_pr_alph.eigh(UPLO::Upper).unwrap(); - C_matr_AO_alph = S_matr_inv_sqrt.dot(&C_matr_MO_alph); - - (orb_E_beta, C_matr_MO_beta) = F_matr_pr_beta.eigh(UPLO::Upper).unwrap(); - C_matr_AO_beta = S_matr_inv_sqrt.dot(&C_matr_MO_beta); - - let delta_E = E_scf_curr - E_scf_prev; - let rms_comm_val = calc_rms_comm_val_uhf(&fps_comm_alph, &fps_comm_beta); - - 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 scf_conv_check(calc_sett, delta_E, rms_comm_val) { - scf.tot_scf_iter = scf_iter; - scf.E_scf_conv = E_scf_curr; - scf.C_matr_conv_alph = C_matr_AO_alph; - scf.C_matr_conv_beta = Some(C_matr_AO_beta); - scf.P_matr_conv_alph = P_matr_alph; - scf.P_matr_conv_beta = Some(P_matr_beta); - scf.orb_E_conv_alph = orb_E_alph; - scf.orb_E_conv_beta = Some(orb_E_beta); - 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_alph = P_matr_alph.clone(); - P_matr_old_beta = P_matr_beta.clone(); - build_P_matr_uhf(&mut P_matr_alph, &C_matr_AO_alph, no_alpha); - build_P_matr_uhf(&mut P_matr_beta, &C_matr_AO_beta, no_beta); - if calc_sett.use_direct_scf { - todo!() - // 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}", "UHF 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 -} +// #[allow(unused)] #[inline(always)] fn scf_conv_check(calc_sett: &CalcSettings, delta_E: f64, rms_comm_val: f64) -> bool { diff --git a/src/main.rs b/src/main.rs index cd021e7..4ae512e 100644 --- a/src/main.rs +++ b/src/main.rs @@ -13,8 +13,7 @@ mod print_utils; use crate::{calc_type::HF_Ref, print_utils::print_header_logo}; use basisset::BasisSet; -use calc_type::{rhf::RHF, uhf::uhf_scf_normal, CalcSettings, DiisSettings}; -use calc_type::HF; +use calc_type::{HF, rhf::RHF, uhf::uhf_scf_normal, CalcSettings, DiisSettings}; use molecule::Molecule; fn main() { @@ -61,9 +60,9 @@ fn main() { let mut rhf = RHF::new(&basis, &calc_sett, HF_Ref::RHF_ref); let _scf = rhf.run_scf(&calc_sett, &mut exec_times, &basis, &mol); exec_times.stop("RHF DIIS indir"); - // exec_times.start("UHF DIIS indir"); - // let _scf = uhf_scf_normal(&calc_sett, &mut exec_times, &basis, &mol); - // exec_times.stop("UHF DIIS indir"); + exec_times.start("UHF DIIS indir"); + let _scf = uhf_scf_normal(&calc_sett, &mut exec_times, &basis, &mol); + exec_times.stop("UHF DIIS indir"); exec_times.stop("Total");