diff --git a/Cargo.toml b/Cargo.toml index 62b4df3b10343b48e770cc319e97c504bcc1245d..1516be6e2ebe6d89697860499d256a4a83155b6d 100644 --- a/Cargo.toml +++ b/Cargo.toml @@ -63,7 +63,6 @@ libxc = { version = "0.1", features = ["api-v7_0", "dynamic_loading"] } geometric-pyo3 = { version = "0.1.1", optional = true } approx = "0.5" -assert_float_eq = "1" [build-dependencies] dunce = "1.0" diff --git a/src/dft/mod.rs b/src/dft/mod.rs index ed702c577b716abd56b912fbf61a9930a926da86..e99db17eaab13d8dd5566b319b2a43d5c21698b9 100644 --- a/src/dft/mod.rs +++ b/src/dft/mod.rs @@ -16,12 +16,11 @@ use itertools::{Itertools, izip}; use tensors::{BasicMatrix, MathMatrix, ParMathMatrix}; // use tensors::external_libs::{general_dgemm_f, matr_copy}; use tensors::matrix_blas_lapack::{_dgemm, _dgemm_full, contract_vxc_0_serial}; +use rest_tensors::matrix_blas_lapack::{omp_get_num_threads_wrapper, omp_set_num_threads_wrapper}; //use numgrid::{self, radial_grid_lmg_bse}; // use self::gen_grids::radial_grid_lmg_bse; use rayon::iter::{IntoParallelRefIterator, IndexedParallelIterator, ParallelIterator, IntoParallelRefMutIterator}; use regex::Regex; -#[cfg(test)] -use crate::basis_io::{BasCell, Basis4Elem, cint_norm_factor,}; use crate::basis_io::{gto_1st_value_batch_serial, gto_1st_value_serial, gto_value, gto_value_matrixfull_serial, gto_value_serial, spheric_gto_1st_value_batch, spheric_gto_value_matrixfull}; use crate::molecule_io::Molecule; @@ -38,8 +37,6 @@ use serde::{Deserialize, Serialize}; use libxc::prelude::*; use crate::dft::libxc_helper::{xc_code_fdqc, xc_func_init, lda_exc_vxc, gga_exc_vxc, mgga_exc_vxc, lda_exc, gga_exc, mgga_exc}; -use rest_tensors::matrix_blas_lapack::{omp_get_num_threads_wrapper, omp_set_num_threads_wrapper}; - #[derive(Clone,Debug, PartialEq, Eq, Serialize, Deserialize)] pub enum DFTType { @@ -303,6 +300,16 @@ impl DFA4REST { xc_func_init(*xc_code, self.spin_channel) } + pub fn init_libxc_and_set_param(&self, xc_code: &usize) -> LibXCFunctional { + let mut func = self.init_libxc(xc_code); + if let Some(omega) = self.omega() { + if func.ext_param_names().contains(&"_omega".to_string()) { + func.set_ext_param_by_name("_omega", omega); + } + } + func + } + pub fn get_hybrid_libxc(dfa_compnt_scf: &Vec,spin_channel:usize) -> f64 { let mut hybrid_coeff = None; for xc_func in dfa_compnt_scf { @@ -1190,6 +1197,78 @@ impl DFA4REST { }) } + fn prepare_dft_quantities( + &self, + grids: &Grids, + spin_channel: usize, + mo: &[MatrixFull; 2], + occ: &[Vec; 2], + use_density_gradient: bool, + ) -> (MatrixFull, RIFull, MatrixFull, MatrixFull, MatrixFull) { + let num_grids = grids.coordinates.len(); + let mut rho: MatrixFull = MatrixFull::empty(); + let mut rhop: RIFull = RIFull::empty(); + let mut lapl: MatrixFull = MatrixFull::empty(); + let mut tau: MatrixFull = MatrixFull::empty(); + + if self.use_kinetic_density() { + (rho, rhop, tau) = grids.prepare_tabulated_density_3(mo, occ, spin_channel); + lapl = MatrixFull::new([num_grids, spin_channel], 0.0); + } else { + (rho, rhop) = if !mo[1].data.is_empty() || spin_channel == 1 { + grids.prepare_tabulated_density_2(mo, occ, spin_channel) + } else { + let mut mo_temp = mo.clone(); + mo_temp[1] = mo_temp[0].clone(); + grids.prepare_tabulated_density_2(&mo_temp, occ, spin_channel) + }; + } + + let sigma = if use_density_gradient { + prepare_tabulated_sigma_rayon(&rhop, spin_channel) + } else { + MatrixFull::empty() + }; + + (rho, rhop, sigma, lapl, tau) + } + + fn compute_exc_by_family(&self, xc_func: &LibXCFunctional, spin_channel: usize, rho: &MatrixFull, sigma: &MatrixFull, lapl: &MatrixFull, tau: &MatrixFull) -> MatrixFull { + let num_grids = rho.size()[0]; + match xc_func.family() { + LibXCFamily::LDA => MatrixFull::from_vec([num_grids, 1], if spin_channel == 1 { + lda_exc(xc_func, rho.data_ref().unwrap()) + } else { + lda_exc(xc_func, rho.transpose().data_ref().unwrap()) + }).unwrap(), + LibXCFamily::GGA | LibXCFamily::HybGGA => MatrixFull::from_vec([num_grids, 1], if spin_channel == 1 { + gga_exc(xc_func, rho.data_ref().unwrap(), sigma.data_ref().unwrap()) + } else { + gga_exc(xc_func, rho.transpose().data_ref().unwrap(), sigma.transpose().data_ref().unwrap()) + }).unwrap(), + LibXCFamily::MGGA | LibXCFamily::HybMGGA => MatrixFull::from_vec([num_grids, 1], if spin_channel == 1 { + mgga_exc(xc_func, rho.data_ref().unwrap(), sigma.data_ref().unwrap(), lapl.data_ref().unwrap(), tau.data_ref().unwrap()) + } else { + mgga_exc(xc_func, rho.transpose().data_ref().unwrap(), sigma.transpose().data_ref().unwrap(), lapl.transpose().data_ref().unwrap(), tau.transpose().data_ref().unwrap()) + }).unwrap(), + xc_family => panic!("{xc_family:?} is not yet implemented"), + } + } + + fn integrate_exc(&self, exc: &MatrixFull, rho: &MatrixFull, weights: &[f64], spin_channel: usize) -> (Vec, [f64; 2]) { + let mut exc_total = vec![0.0; spin_channel]; + let mut total_elec = [0.0; 2]; + for i_spin in 0..spin_channel { + let total_elec_s = total_elec.get_mut(i_spin).unwrap(); + exc_total[i_spin] = izip!(exc.data.iter(), rho.iter_column(i_spin), weights.iter()) + .fold(0.0, |acc, (exc, rho, weight)| { + *total_elec_s += rho * weight; + acc + exc * rho * weight + }); + } + (exc_total, total_elec) + } + pub fn use_kinetic_density(&self) -> bool { self .dfa_compnt_pos @@ -1207,21 +1286,15 @@ impl DFA4REST { let mut vxc_ao = vec![MatrixFull::new([num_basis,num_grids],0.0);spin_channel]; let dt0 = utilities::init_timing(); - let (rho,rhop) = if ! mo[1].data.is_empty() || spin_channel == 1 { // RHF or UHF case - grids.prepare_tabulated_density_2(mo, occ, spin_channel) - } else { // ROHF case - let mut mo_temp = mo.clone(); - mo_temp[1] = mo_temp[0].clone(); - grids.prepare_tabulated_density_2(&mo_temp, occ, spin_channel) - }; + let (rho, rhop, sigma, _lapl, _tau) = self.prepare_dft_quantities( + grids, + spin_channel, + mo, + occ, + self.use_density_gradient(), + ); - let dt2 = utilities::timing(&dt0, Some("evaluate rho and rhop")); - let sigma = if self.use_density_gradient() { - prepare_tabulated_sigma_rayon(&rhop, spin_channel) - } else { - MatrixFull::empty() - }; - + let dt2 = utilities::timing(&dt0, Some("evaluate rho/rhop/sigma")); let mut vrho = MatrixFull::new([num_grids,spin_channel],0.0); let mut vsigma=if self.use_density_gradient() && spin_channel==1 { MatrixFull::new([num_grids,1],0.0) @@ -1230,10 +1303,10 @@ impl DFA4REST { } else { MatrixFull::empty() }; - let dt3 = utilities::timing(&dt2, Some("evaluate sigma")); + let dt3 = utilities::timing(&dt2, Some("init vrho/vsigma")); self.dfa_compnt_scf.iter().zip(self.dfa_paramr_scf.iter()).for_each(|(xc_func,xc_para)| { - let xc_func = self.init_libxc(xc_func); + let xc_func = self.init_libxc_and_set_param(xc_func); match xc_func.family() { LibXCFamily::LDA => { if spin_channel==1 { @@ -1390,17 +1463,7 @@ impl DFA4REST { let dt5 = utilities::timing(&dt4, Some("from vrho -> vxc_ao")); - let mut total_elec = [0.0;2]; - for i_spin in 0..spin_channel { - let mut total_elec_s = total_elec.get_mut(i_spin).unwrap(); - exc_total[i_spin] = izip!(exc.data.iter(),rho.iter_column(i_spin),grids.weights.iter()) - .fold(0.0,|acc,(exc,rho,weight)| { - *total_elec_s += rho*weight; - acc + exc * rho * weight - }); - //exc.data.iter_mut().zip(rho.iter_j(i_spin)).for_each(|(exc,rho)| { - // *exc = *exc* rho - } + let (exc_total, total_elec) = self.integrate_exc(&exc, &rho, &grids.weights, spin_channel); if print_level > 0 { if spin_channel==1 { println!("total electron number: {:16.8}", total_elec[0]) @@ -1489,7 +1552,7 @@ impl DFA4REST { }; self.dfa_compnt_scf.iter().zip(self.dfa_paramr_scf.iter()).for_each(|(xc_func,xc_para)| { - let xc_func = self.init_libxc(xc_func); + let xc_func = self.init_libxc_and_set_param(xc_func); match xc_func.family() { LibXCFamily::LDA => { if spin_channel==1 { @@ -1714,17 +1777,7 @@ impl DFA4REST { // println!("{:16.8},{:16.8},{:16.8}", vsigma[[i,0]],vsigma[[i,1]],vsigma[[i,2]]); //}); - let mut loc_total_elec = [0.0;2]; - for i_spin in 0..spin_channel { - let mut loc_total_elec_s = loc_total_elec.get_mut(i_spin).unwrap(); - loc_exc_total[i_spin] = izip!(loc_exc.data.iter(),loc_rho.iter_column(i_spin),loc_weights.iter()) - .fold(0.0,|acc,(exc,rho,weight)| { - *loc_total_elec_s += rho*weight; - acc + exc * rho * weight - }); - //exc.data.iter_mut().zip(rho.iter_j(i_spin)).for_each(|(exc,rho)| { - // *exc = *exc* rho - } + let (loc_exc_total, loc_total_elec) = self.integrate_exc(&loc_exc, &loc_rho, loc_weights, spin_channel); //if let Some(id) = rayon::current_thread_index() { //} @@ -1824,7 +1877,7 @@ impl DFA4REST { //} self.dfa_compnt_scf.iter().zip(self.dfa_paramr_scf.iter()).for_each(|(xc_func,xc_para)| { - let xc_func = self.init_libxc(xc_func); + let xc_func = self.init_libxc_and_set_param(xc_func); match xc_func.family() { LibXCFamily::LDA => { if spin_channel==1 { @@ -2051,17 +2104,7 @@ impl DFA4REST { // println!("{:16.8},{:16.8},{:16.8}", vsigma[[i,0]],vsigma[[i,1]],vsigma[[i,2]]); //}); - let mut loc_total_elec = [0.0;2]; - for i_spin in 0..spin_channel { - let mut loc_total_elec_s = loc_total_elec.get_mut(i_spin).unwrap(); - loc_exc_total[i_spin] = izip!(loc_exc.data.iter(),loc_rho.iter_column(i_spin),loc_weights.iter()) - .fold(0.0,|acc,(exc,rho,weight)| { - *loc_total_elec_s += rho*weight; - acc + exc * rho * weight - }); - //exc.data.iter_mut().zip(rho.iter_j(i_spin)).for_each(|(exc,rho)| { - // *exc = *exc* rho - } + let (loc_exc_total, loc_total_elec) = self.integrate_exc(&loc_exc, &loc_rho, loc_weights, spin_channel); //if let Some(id) = rayon::current_thread_index() { //} @@ -2088,14 +2131,6 @@ impl DFA4REST { let num_grids = grids.coordinates.len(); let num_basis = dm[0].size[0]; let dt0 = utilities::init_timing(); - let (rho,rhop) = if ! mo[1].data.is_empty() || spin_channel == 1 { // RHF or UHF case - grids.prepare_tabulated_density_2(mo, occ, spin_channel) - } else { // ROHF case - let mut mo_temp = mo.clone(); - mo_temp[1] = mo_temp[0].clone(); - grids.prepare_tabulated_density_2(&mo_temp, occ, spin_channel) - }; - //let (rho,rhop) = grids.prepare_tabulated_density(dm, spin_channel); let use_density_gradient = post_xc.iter().fold(false,|flag, x| { let code = DFA4REST::xc_func_init_fdqc(x,spin_channel); let x_flag = code.iter().fold(false, |flag, xc_code| { @@ -2104,26 +2139,29 @@ impl DFA4REST { }); flag || x_flag }); - let sigma = if use_density_gradient { - prepare_tabulated_sigma_rayon(&rhop, spin_channel) - } else { - MatrixFull::empty() - }; + let (rho, rhop, sigma, _lapl, _tau) = self.prepare_dft_quantities( + grids, + spin_channel, + mo, + occ, + use_density_gradient, + ); + let _dt2 = utilities::timing(&dt0, Some("evaluate rho/rhop/sigma")); + //let (rho,rhop) = grids.prepare_tabulated_density(dm, spin_channel); post_xc.iter().for_each(|x| { let mut exc = MatrixFull::new([num_grids,1],0.0); - let mut exc_total =[0.0,0.0]; + let mut exc_total = [0.0, 0.0]; let code = DFA4REST::xc_func_init_fdqc(x,spin_channel); //println!("debug xc_code: {:?}", &code); code.iter().for_each(|xc_code| { exc.par_self_scaled_add(&self.xc_exc_code(xc_code, &rho, &sigma, spin_channel),1.0); }); - for i_spin in 0..spin_channel { - exc_total[i_spin] = izip!(exc.data.iter(),rho.iter_column(i_spin),grids.weights.iter()) - .fold(0.0,|acc,(exc,rho,weight)| { - acc + exc * rho * weight - }); - }; + let (exc_total_vec, _) = self.integrate_exc(&exc, &rho, &grids.weights, spin_channel); + exc_total[0] = exc_total_vec[0]; + if spin_channel == 2 { + exc_total[1] = exc_total_vec[1]; + } //println!("exc_total: {:?}", &exc_total); post_xc_energy.push(exc_total); @@ -2196,13 +2234,13 @@ impl DFA4REST { //} }); - let (rho,rhop) = if ! mo[1].data.is_empty() || spin_channel == 1 { // RHF or UHF case - grids.prepare_tabulated_density_2(mo, occ, spin_channel) - } else { // ROHF case - let mut mo_temp = mo.clone(); - mo_temp[1] = mo_temp[0].clone(); - grids.prepare_tabulated_density_2(&mo_temp, occ, spin_channel) - }; + let (rho, rhop, sigma, _lapl, _tau) = self.prepare_dft_quantities( + grids, + spin_channel, + mo, + occ, + use_density_gradient, + ); str_lines[1..].par_iter_mut().zip(rho.par_iter_column(0)).for_each(|(line,rho)| { line.push_str(&format!("{:20.10} ", rho)); @@ -2213,11 +2251,6 @@ impl DFA4REST { }); } - let sigma = if use_density_gradient { - prepare_tabulated_sigma_rayon(&rhop, spin_channel) - } else { - MatrixFull::empty() - }; if use_density_gradient { str_lines[1..].par_iter_mut().zip(rhop.par_iter_slices_x(0, 0)).for_each(|(line,rhop_i)| { line.push_str(&format!("{:20.10} ", rhop_i)); @@ -2303,38 +2336,30 @@ impl DFA4REST { let num_grids = grids.coordinates.len(); let num_basis = dm[0].size[0]; let dt0 = utilities::init_timing(); - let (rho,rhop) = if ! mo[1].data.is_empty() || spin_channel == 1 { // RHF or UHF case - grids.prepare_tabulated_density_2(mo, occ, spin_channel) - } else { // ROHF case - let mut mo_temp = mo.clone(); - mo_temp[1] = mo_temp[0].clone(); - grids.prepare_tabulated_density_2(&mo_temp, occ, spin_channel) - }; - //let (rho,rhop) = grids.prepare_tabulated_density(dm, spin_channel); let use_density_gradient = xc_code_list.iter().fold(false,|flag, xc_code| { let xc_func = self.init_libxc(xc_code); flag || !matches!(xc_func.family(), LibXCFamily::LDA | LibXCFamily::HybLDA) }); - let sigma = if use_density_gradient { - prepare_tabulated_sigma_rayon(&rhop, spin_channel) - } else { - MatrixFull::empty() - }; + let (rho, rhop, sigma, _lapl, _tau) = self.prepare_dft_quantities( + grids, + spin_channel, + mo, + occ, + use_density_gradient, + ); + let _dt2 = utilities::timing(&dt0, Some("evaluate rho/rhop/sigma")); + //let (rho,rhop) = grids.prepare_tabulated_density(dm, spin_channel); xc_code_list.iter().for_each(|xc_code| { //let mut exc = MatrixFull::new([num_grids,1],0.0); - let mut exc_total =[0.0,0.0]; - //let code = DFA4REST::xc_func_init_fdqc(x,spin_channel); - //println!("debug xc_code: {:?}", &code); - //code.iter().for_each(|xc_code| { + let mut exc_total = [0.0, 0.0]; let exc = self.xc_exc_code(xc_code, &rho, &sigma, spin_channel); //}); - for i_spin in 0..spin_channel { - exc_total[i_spin] = izip!(exc.data.iter(),rho.iter_column(i_spin),grids.weights.iter()) - .fold(0.0,|acc,(exc,rho,weight)| { - acc + exc * rho * weight - }); - }; + let (exc_total_vec, _) = self.integrate_exc(&exc, &rho, &grids.weights, spin_channel); + exc_total[0] = exc_total_vec[0]; + if spin_channel == 2 { + exc_total[1] = exc_total_vec[1]; + } //println!("exc_total: {:?}", &exc_total); xc_energy.push(exc_total); @@ -2350,29 +2375,14 @@ impl DFA4REST { let mut exc = MatrixFull::new([num_grids,1],0.0); let mut exc_total = vec![0.0;spin_channel]; let dt0 = utilities::init_timing(); - let mut rho: MatrixFull = MatrixFull::empty(); - let mut rhop: RIFull = RIFull::empty(); - let mut tau: MatrixFull = MatrixFull::empty(); - let mut lapl: MatrixFull = MatrixFull::empty(); - if self.use_kinetic_density() { - (rho, rhop, tau) = grids.prepare_tabulated_density_3(mo, occ, spin_channel); - lapl = MatrixFull::new([num_grids, spin_channel], 0.0); - } else { - (rho,rhop) = if ! mo[1].data.is_empty() || spin_channel == 1 { // RHF or UHF case - grids.prepare_tabulated_density_2(mo, occ, spin_channel) - } else { // ROHF case - let mut mo_temp = mo.clone(); - mo_temp[1] = mo_temp[0].clone(); - grids.prepare_tabulated_density_2(&mo_temp, occ, spin_channel) - } - }; - // let (rho,rhop) = grids.prepare_tabulated_density_2(mo, occ, spin_channel); - let dt2 = utilities::timing(&dt0, Some("evaluate rho and rhop")); - let sigma = if self.use_density_gradient() { - prepare_tabulated_sigma_rayon(&rhop, spin_channel) - } else { - MatrixFull::empty() - }; + let (rho, rhop, sigma, lapl, tau) = self.prepare_dft_quantities( + grids, + spin_channel, + mo, + occ, + self.use_density_gradient(), + ); + let dt2 = utilities::timing(&dt0, Some("evaluate rho/rhop/sigma")); if iop==0 { // for the SCF energy //let rho_trans = if spin_channel== 1 { @@ -2381,89 +2391,20 @@ impl DFA4REST { // Some(rho.transpose()) //}; self.dfa_compnt_scf.iter().zip(self.dfa_paramr_scf.iter()).for_each(|(xc_func,xc_para)| { - let xc_func = self.init_libxc(xc_func); - match xc_func.family() { - LibXCFamily::LDA => { - let tmp_exc = MatrixFull::from_vec([num_grids,1], - if spin_channel==1 { - lda_exc(&xc_func,rho.data_ref().unwrap()) - } else { - lda_exc(&xc_func,rho.transpose().data_ref().unwrap()) - } - ).unwrap(); - exc.par_self_scaled_add(&tmp_exc,*xc_para); - }, - LibXCFamily::GGA | LibXCFamily::HybGGA => { - let tmp_exc = MatrixFull::from_vec([num_grids,1], - if spin_channel==1 { - gga_exc(&xc_func,rho.data_ref().unwrap(),sigma.data_ref().unwrap()) - } else { - gga_exc(&xc_func,rho.transpose().data_ref().unwrap(),sigma.transpose().data_ref().unwrap()) - } - ).unwrap(); - exc.par_self_scaled_add(&tmp_exc,*xc_para); - }, - LibXCFamily::MGGA | LibXCFamily::HybMGGA => { - let tmp_exc = MatrixFull::from_vec( - [num_grids, 1], - if spin_channel==1 { - mgga_exc(&xc_func,rho.data_ref().unwrap(),sigma.data_ref().unwrap(), lapl.data_ref().unwrap(), tau.data_ref().unwrap()) - } else { - mgga_exc(&xc_func,rho.transpose().data_ref().unwrap(),sigma.transpose().data_ref().unwrap(), lapl.transpose().data_ref().unwrap(), tau.transpose().data_ref().unwrap()) - } - ).unwrap(); - exc.par_self_scaled_add(&tmp_exc,*xc_para); - }, - xc_family => panic!("{xc_family:?} is not yet implemented"), - } + let xc_func = self.init_libxc_and_set_param(xc_func); + let tmp_exc = self.compute_exc_by_family(&xc_func, spin_channel, &rho, &sigma, &lapl, &tau); + exc.par_self_scaled_add(&tmp_exc, *xc_para); }); } else if iop==1 { // for the post-SCF energy calculation if let (Some(dfa_paramr),Some(dfa_compnt)) = (&self.dfa_paramr_pos, &self.dfa_compnt_pos) { dfa_compnt.iter().zip(dfa_paramr.iter()).for_each(|(xc_func,xc_para)| { - let xc_func = self.init_libxc(xc_func); - match xc_func.family() { - LibXCFamily::LDA => { - let tmp_exc = MatrixFull::from_vec([num_grids,1], - if spin_channel==1 { - lda_exc(&xc_func,rho.data_ref().unwrap()) - } else { - lda_exc(&xc_func,rho.transpose().data_ref().unwrap()) - } - ).unwrap(); - exc.par_self_scaled_add(&tmp_exc,*xc_para); - }, - LibXCFamily::GGA | LibXCFamily::HybGGA => { - let tmp_exc = MatrixFull::from_vec([num_grids,1], - if spin_channel==1 { - gga_exc(&xc_func,rho.data_ref().unwrap(),sigma.data_ref().unwrap()) - } else { - gga_exc(&xc_func,rho.transpose().data_ref().unwrap(),sigma.transpose().data_ref().unwrap()) - } - ).unwrap(); - exc.par_self_scaled_add(&tmp_exc,*xc_para); - }, - xc_family => panic!("{xc_family:?} is not yet implemented") - } + let xc_func = self.init_libxc_and_set_param(xc_func); + let tmp_exc = self.compute_exc_by_family(&xc_func, spin_channel, &rho, &sigma, &lapl, &tau); + exc.par_self_scaled_add(&tmp_exc, *xc_para); }); } } - let mut total_elec = [0.0;2]; - for i_spin in 0..spin_channel { - let mut total_elec_s = total_elec.get_mut(i_spin).unwrap(); - exc_total[i_spin] = izip!(exc.data.iter(),rho.iter_column(i_spin),grids.weights.iter()) - .fold(0.0,|acc,(exc,rho,weight)| { - *total_elec_s += rho*weight; - acc + exc * rho * weight - }); - //exc.data.iter_mut().zip(rho.iter_j(i_spin)).for_each(|(exc,rho)| { - // *exc = *exc* rho - } - //if spin_channel==1 { - // println!("total electron number: {:16.8}", total_elec[0]) - //} else { - // println!("electron number in alpha-channel: {:12.8}", total_elec[0]); - // println!("electron number in beta-channel: {:12.8}", total_elec[1]); - //} + let (exc_total, _total_elec) = self.integrate_exc(&exc, &rho, &grids.weights, spin_channel); let global_exc_total = if let Some(mpi_op) = &mpi_operator { let my_rank = mpi_op.rank; let mut global_exc_total = mpi_reduce(&mpi_op.world, &exc_total , 0, &SystemOperation::sum()); @@ -2478,33 +2419,9 @@ impl DFA4REST { } pub fn xc_exc_code(&self, xc_code: &usize, rho: &MatrixFull, sigma:&MatrixFull, spin_channel: usize) -> MatrixFull { - let xc_func = self.init_libxc(xc_code); + let xc_func = self.init_libxc_and_set_param(xc_code); let num_grids = rho.size()[0]; - let tmp_exc = match xc_func.family() { - LibXCFamily::LDA => { - MatrixFull::from_vec([num_grids,1], - if spin_channel==1 { - lda_exc(&xc_func,rho.data_ref().unwrap()) - } else { - lda_exc(&xc_func,rho.transpose().data_ref().unwrap()) - } - ).unwrap() - //exc.par_self_scaled_add(&tmp_exc,*xc_para); - }, - LibXCFamily::GGA | LibXCFamily::HybGGA => { - //println!("debug, {:?}, {:?}, {:?}", xc_code, xc_func, sigma.data.len()); - MatrixFull::from_vec([num_grids,1], - if spin_channel==1 { - gga_exc(&xc_func,rho.data_ref().unwrap(),sigma.data_ref().unwrap()) - } else { - gga_exc(&xc_func,rho.transpose().data_ref().unwrap(),sigma.transpose().data_ref().unwrap()) - } - ).unwrap() - //exc.par_self_scaled_add(&tmp_exc,*xc_para); - }, - xc_family => panic!("{xc_family:?} is not yet implemented"), - }; - tmp_exc + self.compute_exc_by_family(&xc_func, spin_channel, rho, sigma, &MatrixFull::empty(), &MatrixFull::empty()) } } @@ -3863,382 +3780,8 @@ pub fn numerical_density_rayon(grid: &Grids, mol: &Molecule, dm: &Vec = HashMap::new(); - alpha_min_h.insert(angular,0.122); - let mut alpha_max_h: f64 = 0.122; - let mut basis4elem = vec![Basis4Elem { - electron_shells: vec![ - BasCell { - function_type: None, - region: None, - angular_momentum: vec![angular as i32], - exponents: vec![0.122], - coefficients: vec![vec![1.0/cint_norm_factor(angular as i32, 0.122)]], - native_coefficients: vec![vec![1.0/cint_norm_factor(angular as i32, 0.122)]] - } - ], - references: None, - ecp_electrons: None, - ecp_potentials: None, - global_index: (0,0) - }]; - let mut center_coordinates_bohr = vec![(0.0,0.0,0.0)]; - let mut proton_charges = vec![1]; - let grids = Grids::build_nonstd( - center_coordinates_bohr.clone(), - proton_charges.clone(), - vec![alpha_min_h], - vec![alpha_max_h], &mut None); - let mut total_density = 0.0; - let mut count:usize =0; - grids.coordinates.iter().zip(grids.weights.iter()).for_each(|(r,w)| { - let mut density_r_sum = 0.0; - let mut density_r:Vec = vec![]; - basis4elem.iter().zip(center_coordinates_bohr.iter()).for_each(|(elem,geom_nonstd)| { - let geom = [geom_nonstd.0,geom_nonstd.1,geom_nonstd.2]; - let mut tmp_geom = [0.0;3]; - tmp_geom.iter_mut().zip(geom.iter()).for_each(|value| {*value.0 = *value.1}); - //density_r.extend(gto_value(r, &tmp_geom, elem, &mol.ctrl.basis_type)); - //let tmp_vec = gto_value(r, &tmp_geom, elem, &"spheric".to_string()); - let tmp_vec = gto_value(r, &tmp_geom, elem, &"spheric".to_string()); - //println!("debug 0: len {}", &tmp_vec.len()); - density_r.extend(tmp_vec); - //if count<=10 {println!("{:?},{:?},{:?}", density_r,elem.electron_shells[0].exponents, elem.electron_shells[0].coefficients[0])}; - }); - //println!{"debug 1"}; - let mut density_rr = MatrixFull::from_vec([num_basis,1],density_r).unwrap(); - //println!{"debug 2"}; - dm.iter_mut().for_each(|dm_s| { - let mut tmp_mat = MatrixFull::new([num_basis,1],0.0); - tmp_mat.lapack_dgemm(&mut density_rr, dm_s, 'T', 'N', 1.0, 0.0); - if count<=10 {println!("count: {},{:?}",count, &tmp_mat.data)}; - density_r_sum += tmp_mat.data.iter().zip(density_rr.data.iter()).fold(0.0, |acc,(a,b)| {acc + a*b}); - }); - if count<=10 {println!("{:?},{},{}", r,w,density_r_sum)}; - count += 1; - //println!{"debug 3"}; - total_density += density_r_sum * w; - }); - - println!("Total density: {}", total_density); - -} - -#[test] -fn test_libxc() { - let mut rho:Vec = vec![0.1,0.2,0.3,0.4,0.5,0.6,0.8]; - let sigma:Vec = vec![0.2,0.3,0.4,0.5,0.6,0.7]; - //&rho.par_iter().for_each(|c| {println!("{:16.8}",c)}); - //let mut exc:Vec = vec![0.0,0.0,0.0,0.0,0.0]; - //let mut vrho:Vec = vec![0.0,0.0,0.0,0.0,0.0]; - //let func_id: usize = ffi_xc::XC_GGA_X_XPBE as usize; - let spin_channel: usize = 1; - - let mut my_xc = DFA4REST::parse_scf("lda_x_slater", spin_channel); - //let mut my_xc = XcFuncType::xc_func_init_fdqc(&"pw-lda",spin_channel); - - - - let mut exc = MatrixFull::new([rho.len()/spin_channel,1],0.0); - let mut vrho = MatrixFull::new([rho.len()/spin_channel,spin_channel],0.0); - my_xc.dfa_compnt_scf.iter().zip(my_xc.dfa_paramr_scf.iter()).for_each(|(xc_func, xc_para)| { - //let mut new_vec = vec![-0.34280861230056237, -0.43191178672272906, -0.494415573788165, -0.5441747517896713, -0.586194481347579, -0.622924588811561, -0.6856172246011247]; - //let tmp_c = (new_vec.as_mut_ptr(), new_vec.len(),new_vec.capacity()); - //let new_vec = unsafe{Vec::from_raw_parts(tmp_c.0, tmp_c.1, tmp_c.2)}; - //new_vec.par_iter().for_each(|c| {println!("{:16.8e}",c)}); - let xc_func = my_xc.init_libxc(xc_func); - - let (tmp_exc, tmp_vrho) = lda_exc_vxc(&xc_func,&rho); - //let tmp_exc_2 = tmp_exc.clone(); - //println!("WARNNING:: unsolved rayon par_iter problem. It should be relevant to be the fact that tmp_exc is prepared by libxc via ffi"); - //println!("tmp_vec_2 copied from tmp_exc: {:?},{},{}", &tmp_exc_2, tmp_exc_2.len(),tmp_exc_2.capacity()); - //println!("tmp_vec"); - //&tmp_exc.iter().for_each(|c| { - // println!("{:16.8e}",c); - //}); - //println!("tmp_vec_2"); - //&tmp_exc_2.par_iter().for_each(|c| { - // println!("{:16.8e}",c); - //}); - //println!("tmp_vec_2"); - //&tmp_exc.iter().for_each(|c| { - // println!("{:16.8e}",c); - //}); - let mut tmp_exc = MatrixFull::from_vec([rho.len()/spin_channel,1],tmp_exc).unwrap(); - let mut tmp_vrho = MatrixFull::from_vec([rho.len()/spin_channel,spin_channel],tmp_vrho).unwrap(); - //println!("{:?}", &tmp_exc.data); - //println!("{:?}", &tmp_vrho.data); - exc.par_self_scaled_add(&tmp_exc,*xc_para); - vrho.par_self_scaled_add(&tmp_vrho,*xc_para); - //exc.data.par_iter_mut().zip(tmp_exc.data.par_iter()).for_each(|(c,p)| { - // println!("{:16.8e},{:16.8e}",c,p); - //}); - //println!("{:?}", &exc.data); - //println!("{:?}", &vrho.data); - }); - println!("{:?}", exc.data); - println!("{:?}", vrho.data); -} - -#[test] -#[ignore] -fn read_grid() { - let mut grids_file = std::fs::File::open("/home/igor/Documents/Package-Pool/Rust/rest/grids").unwrap(); - let mut content = String::new(); - grids_file.read_to_string(&mut content); - //println!("{}",&content); - let re1 = Regex::new(r"(?x)\s* - (?P[\+-]?\d+.\d+[eE][\+-]?\d+)\s*,# the 'x' position - \s+ - (?P[\+-]?\d+.\d+[eE][\+-]?\d+)\s*,# the 'y' position - \s+ - (?P[\+-]?\d+.\d+[eE][\+-]?\d+)\s*,# the 'z' position - \s+ - (?P[\+-]?\d+.\d+[eE][\+-]?\d+)\s*# the 'w' weight - \s*\n").unwrap(); - //if let Some(cap) = re1.captures(&content) { - // println!("{:?}", &cap) - //} - for cap in re1.captures_iter(&content) { - let x:f64 = cap[1].parse().unwrap(); - let y:f64 = cap[2].parse().unwrap(); - let z:f64 = cap[3].parse().unwrap(); - let w:f64 = cap[4].parse().unwrap(); - println!("{:16.8} {:16.8} {:16.8} {:16.8}", x,y,z,w); - } -} -#[test] -fn debug_transpose() { - let len_a = 111_usize; - let len_b = 40000_usize; - let orig_a:Vec = (0..len_a*len_b).map(|i| {i as f64}).collect(); - let a_mat = MatrixFull::from_vec([len_a,len_b],orig_a).unwrap(); - let dt0 = utilities::init_timing(); - let b_mat = a_mat.transpose_and_drop(); - //b_mat.formated_output(10, "full"); - let dt1 = utilities::timing(&dt0, Some("old transpose")); - let orig_a:Vec = (0..len_a*len_b).map(|i| {i as f64}).collect(); - let a_mat = MatrixFull::from_vec([len_a,len_b],orig_a).unwrap(); - let dt0 = utilities::init_timing(); - let c_mat = a_mat.transpose_and_drop(); - //b_mat.formated_output(10, "full"); - let dt1 = utilities::timing(&dt0, Some("new transpose")); - b_mat.data.iter().zip(c_mat.data.iter()).for_each(|(b,c)| { - assert!(*b==*c); - }); -} - - - #[test] fn test_balancing() { let dd = balancing(550, 23); println!("{:?}",dd); } - -#[test] -fn test_rsh_cam_coeff_raw() { - // Test raw libxc CAM coefficients match pyscf expectations - // Each entry: (libxc_id, expected_omega, expected_alpha, expected_beta, expected_hyb) - // hyb = alpha + beta (pyscf convention for RSH) - let ref_data = vec![ - (464, 0.300, 1.0, -0.842294, 0.157706), // WB97X - (466, 0.300, 1.0, -0.833, 0.167), // WB97X-V - (471, 0.200, 1.0, -0.777964, 0.222036), // WB97X-D - (433, 0.330, 0.65, -0.46, 0.19), // CAM-B3LYP - (400, 0.330, 1.0, -1.0, 0.0), // LC-BLYP - (478, 0.400, 1.0, -1.0, 0.0), // LC-wPBE - (428, 0.110, 0.0, 0.25, 0.25), // HSE06 - (427, 0.1061, 0.0, 0.25, 0.25), // HSE03 - ]; - - for (libxc_id, exp_omega, exp_alpha, exp_beta, exp_hyb) in &ref_data { - let func = xc_func_init(*libxc_id, 1); - assert!(func.is_hyb_cam(), "ID {} should be RSH/CAM", libxc_id); - let (omega, alpha, beta) = func.cam_coef().unwrap_or((0.0, 0.0, 0.0)); - let hyb = alpha + beta; - assert!((omega - exp_omega).abs() < 1e-3, "ID {}: omega mismatch: got {} expected {}", libxc_id, omega, exp_omega); - assert!((alpha - exp_alpha).abs() < 1e-3, "ID {}: alpha mismatch: got {} expected {}", libxc_id, alpha, exp_alpha); - assert!((beta - exp_beta).abs() < 1e-3, "ID {}: beta mismatch: got {} expected {}", libxc_id, beta, exp_beta); - assert!((hyb - exp_hyb).abs() < 1e-4, "ID {}: hyb mismatch: got {} expected {}", libxc_id, hyb, exp_hyb); - println!("ID {}: omega={:.4} alpha={:.4} beta={:+.4} hyb={:.6} OK", libxc_id, omega, alpha, beta, hyb); - } - println!("All {} RSH libxc raw CAM coefficient tests passed.", ref_data.len()); -} - -#[test] -fn test_rsh_parameters_wb97x() { - let dfa = DFA4REST::new("wb97x", 1, 0); - assert!(dfa.is_hybrid()); - assert!(dfa.is_rsh()); - assert!((dfa.omega().unwrap() - 0.3).abs() < 1e-9); - assert!((dfa.rsh_alpha().unwrap() - 1.0).abs() < 1e-9); - assert!((dfa.dfa_hybrid_scf - 0.157706).abs() < 1e-5); - assert_eq!(dfa.dfa_compnt_scf, vec![464]); -} - -#[test] -fn test_rsh_parameters_cam_b3lyp() { - let dfa = DFA4REST::new("cam-b3lyp", 1, 0); - assert!(dfa.is_hybrid()); - assert!(dfa.is_rsh()); - assert!((dfa.omega().unwrap() - 0.33).abs() < 1e-9); - assert!((dfa.rsh_alpha().unwrap() - 0.65).abs() < 1e-5); - assert!((dfa.dfa_hybrid_scf - 0.19).abs() < 1e-5); - assert_eq!(dfa.dfa_compnt_scf, vec![433]); - // also test case variation - let dfa2 = DFA4REST::new("CAM-B3LYP", 1, 0); - assert!((dfa.omega().unwrap() - dfa2.omega().unwrap()).abs() < 1e-9); - assert!((dfa.dfa_hybrid_scf - dfa2.dfa_hybrid_scf).abs() < 1e-9); -} - -#[test] -fn test_rsh_parameters_lc_blyp() { - let dfa = DFA4REST::new("lc-blyp", 1, 0); - assert!(dfa.is_hybrid()); - assert!(dfa.is_rsh()); - assert!((dfa.omega().unwrap() - 0.33).abs() < 1e-9); - assert!((dfa.rsh_alpha().unwrap() - 1.0).abs() < 1e-9); - // LC-BLYP has hyb=0 (alpha=1, beta=-1 → hyb=0, LR-only HF) - assert!((dfa.dfa_hybrid_scf - 0.0).abs() < 1e-8); - assert_eq!(dfa.dfa_compnt_scf, vec![400]); -} - -#[test] -fn test_rsh_parameters_hse06() { - let dfa = DFA4REST::new("hse06", 1, 0); - assert!(dfa.is_hybrid()); - assert!(dfa.is_rsh()); - assert!((dfa.omega().unwrap() - 0.11).abs() < 1e-9); - // HSE06 has alpha=0 (no additional LR-HF), beta=0.25 (SR-HF only) - assert!((dfa.rsh_alpha().unwrap() - 0.0).abs() < 1e-5); - assert!((dfa.dfa_hybrid_scf - 0.25).abs() < 1e-5); - assert_eq!(dfa.dfa_compnt_scf, vec![428]); -} - -#[test] -fn test_non_rsh_b3lyp() { - let dfa = DFA4REST::new("b3lyp", 1, 0); - assert!(dfa.is_hybrid(), "B3LYP should be hybrid"); - assert!(!dfa.is_rsh(), "B3LYP should NOT be range-separated"); - assert!(dfa.omega().is_none()); - assert!(dfa.rsh_alpha().is_none()); -} - -#[test] -fn test_non_rsh_pbe() { - let dfa = DFA4REST::new("pbe", 1, 0); - assert!(!dfa.is_hybrid()); - assert!(!dfa.is_rsh()); - assert!(dfa.omega().is_none()); - assert!(dfa.rsh_alpha().is_none()); - assert_eq!(dfa.dfa_hybrid_scf, 0.0); -} - -#[test] -fn test_cam_b3lyp_case_insensitive() { - let names = ["cam-b3lyp", "CAMB3LYP"]; - let mut dfas = vec![]; - for name in &names { - dfas.push(DFA4REST::new(name, 1, 0)); - } - for i in 1..dfas.len() { - assert!((dfas[0].omega().unwrap() - dfas[i].omega().unwrap()).abs() < 1e-9); - assert!((dfas[0].rsh_alpha().unwrap() - dfas[i].rsh_alpha().unwrap()).abs() < 1e-9); - assert!((dfas[0].dfa_hybrid_scf - dfas[i].dfa_hybrid_scf).abs() < 1e-9); - } -} - -#[test] -fn test_rsh_spin_polarized() { - // Test RSH initialization with spin=2 (open-shell) - let dfa_closed = DFA4REST::new("cam-b3lyp", 1, 0); - let dfa_open = DFA4REST::new("cam-b3lyp", 2, 0); - // RSH parameters should be spin-independent - assert!((dfa_closed.omega().unwrap() - dfa_open.omega().unwrap()).abs() < 1e-9); - assert!((dfa_closed.rsh_alpha().unwrap() - dfa_open.rsh_alpha().unwrap()).abs() < 1e-9); - assert!((dfa_closed.dfa_hybrid_scf - dfa_open.dfa_hybrid_scf).abs() < 1e-9); -} - -#[test] -fn test_rsh_hyb_override() { - // Verify that for RSH functionals, dfa_hybrid_scf is overridden to alpha+beta - // NOT the raw xc_hyb_exx_coeff (which may return 1.0) - let dfa = DFA4REST::new("lc-blyp", 1, 0); - // LC-BLYP: raw xc_hyb_exx_coef = 1.0 (because it's a hybrid), - // but alpha+beta = 0.0 (pure LR-HF, no SR-HF) - // The override should set hyb to 0.0 - assert!((dfa.dfa_hybrid_scf - 0.0).abs() < 1e-8, - "LC-BLYP hyb should be 0 (alpha+beta=0), got {}", dfa.dfa_hybrid_scf); - - let dfa2 = DFA4REST::new("cam-b3lyp", 1, 0); - // CAM-B3LYP: raw xc_hyb_exx_coef = ?, alpha+beta = 0.19 - assert!((dfa2.dfa_hybrid_scf - 0.19).abs() < 1e-3, - "CAM-B3LYP hyb should be ~0.19 (alpha+beta), got {}", dfa2.dfa_hybrid_scf); -} - -#[test] -fn test_rsh_nonstd_parse() { - // Test parse_scf_nonstd with RSH components (e.g. custom mixing) - let dfa = DFA4REST::parse_scf_nonstd( - &vec!["HYB_GGA_XC_WB97X_V".to_string()], - &vec![0.5], // 50% scaling - &0.167, // hyb = alpha+beta - 1, - ); - assert!(dfa.is_rsh(), "Nonstd WB97X-V should be RSH"); - assert!((dfa.omega().unwrap() - 0.3).abs() < 1e-9); - assert!((dfa.rsh_alpha().unwrap() - 1.0).abs() < 1e-9); - assert!((dfa.dfa_hybrid_scf - 0.167).abs() < 1e-5); - - // Test with HSE06 nonstd - let dfa2 = DFA4REST::parse_scf_nonstd( - &vec!["HYB_GGA_XC_HSE06".to_string()], - &vec![1.0], - &0.25, - 1, - ); - assert!(dfa2.is_rsh()); - assert!((dfa2.omega().unwrap() - 0.11).abs() < 1e-9); - assert!((dfa2.rsh_alpha().unwrap() - 0.0).abs() < 1e-5); -} - -#[test] -fn test_rsh_summary_does_not_panic() { - // Ensure summary() works without panicking for RSH and non-RSH - for name in &["cam-b3lyp", "lc-blyp", "hse06", "b3lyp", "pbe"] { - let dfa = DFA4REST::new(name, 1, 0); - dfa.summary(); - } -} - -#[test] -fn test_all_rsh_via_auto_resolver() { - // Test that the auto-resolver (generic libxc name lookup) works for RSH - let dfa = DFA4REST::new("HYB_GGA_XC_WB97X_V", 1, 0); - assert!(dfa.is_rsh()); - assert!((dfa.omega().unwrap() - 0.3).abs() < 1e-9); - assert_eq!(dfa.dfa_compnt_scf, vec![466]); - - let dfa2 = DFA4REST::new("HYB_GGA_XC_HSE06", 1, 0); - assert!(dfa2.is_rsh()); - assert!((dfa2.omega().unwrap() - 0.11).abs() < 1e-9); - assert_eq!(dfa2.dfa_compnt_scf, vec![428]); -} diff --git a/src/dft/parse_xc/family_bdh_omega.json b/src/dft/parse_xc/family_bdh_omega.json new file mode 100644 index 0000000000000000000000000000000000000000..dc6e913d24e747a3f5b2c042bfedbcf1818e0cfd --- /dev/null +++ b/src/dft/parse_xc/family_bdh_omega.json @@ -0,0 +1,26 @@ +{ + "wB2PLYP": { + "code": "0.53*HF + 0.47*LR_HF(0.3) + 0.47*GGA_X_ITYH, 0.73*LYP + 0.27*MP2", + "ref": "10.1021/acs.jctc.9b00013", + "tested": "ORCA 5.0.4", + "alias": ["ωB2PLYP"] + }, + "wB2GP_PLYP": { + "code": "0.65*HF + 0.35*LR_HF(0.27) + 0.35*GGA_X_ITYH, 0.64*LYP + 0.36*MP2", + "ref": "10.1021/acs.jctc.9b00013", + "tested": "ORCA 5.0.4", + "alias": ["ωB2GP_PLYP"] + }, + "wB88PP86": { + "code": "0.65*HF + 0.35*LR_HF(0.2) + 0.35*GGA_X_ITYH, 0.58*GGA_C_P86VWN + 0.42*MP2", + "ref": "10.1021/acs.jctc.1c00535", + "tested": "ORCA 5.0.4", + "alias": ["ωB88PP86"] + }, + "wPBEPP86": { + "code": "0.70*HF + 0.30*LR_HF(0.18) + 0.30*GGA_X_ITYH_PBE, 0.68*GGA_C_P86VWN + 0.48*MP2", + "ref": "10.1021/acs.jctc.1c00535", + "tested": "ORCA 5.0.4", + "alias": ["ωPBEPP86"] + } +} \ No newline at end of file diff --git a/src/dft/parse_xc/parse.rs b/src/dft/parse_xc/parse.rs index 1f6cd7a50a956624f943c1e13c8fad48efb5c8a3..ac096aded0b430ebaeb1b687f147fcff7b985000 100644 --- a/src/dft/parse_xc/parse.rs +++ b/src/dft/parse_xc/parse.rs @@ -184,6 +184,10 @@ impl DFAComponent { // #[cached] pub fn get_hybrid(&self, spin_channel: usize) -> f64 { + if self.is_rsh() { + let (omega, alpha, beta) = self.get_rsh(spin_channel); + return self.factor * (alpha + beta); + } if self.component_type == ComponentType::HF { return self.factor; } else if self.component_type == ComponentType::Libxc { @@ -201,8 +205,48 @@ impl DFAComponent { } } + pub fn get_rsh(&self, spin_channel: usize) -> (Option, f64, f64) { + if !self.is_rsh() { + return (None, self.get_hybrid(spin_channel), 0.0); + } + if self.component_type == ComponentType::Libxc { + let xcfunc = xc_func_init(self.id, spin_channel); + if let Some((omega, alpha, beta)) = xcfunc.cam_coef() { + return (Some(omega), alpha, beta); + } else { + panic!("Error: component {} is detected as RSH but cam_coef is not available in libxc", self.func); + } + } else if self.component_type == ComponentType::RSHF { + let param = &self.param_positional; + match self.func.as_str() { + "LR_HF" | "SR_HF" => { + if param.len() != 1 { + panic!("Error: SR/LR HF component with func {} should have exactly 1 positional parameter, but got {}", self.func, param.len()); + } + }, + "RSH" => { + if param.len() != 3 { + panic!("Error: RSH component should have exactly 3 positional parameters (omega, alpha, beta), but got {}", param.len()); + } + }, + _ => {}, + } + match self.func.as_str() { + "LR_HF" => return (Some(param[0]), self.factor, -1.0*self.factor), + "SR_HF" => return (Some(param[0]), 0.0, self.factor), + "RSH" => return (Some(param[0]), param[1], param[2]), + _ => return (None, 0.0, 0.0), + } + } else { + return (None, 0.0, 0.0); + } + } + pub fn is_rsh(&self) -> bool { match self.component_type { + ComponentType::RSHF => { + true + }, ComponentType::Libxc => { let xcfunc = xc_func_init(self.id, 1); xcfunc.is_hyb_cam() @@ -241,6 +285,9 @@ impl DFAComponent { pub fn canonicalize(&self) -> Self { match self.component_type { + ComponentType::RSHF => { + return self.normalize_rsh(); + } ComponentType::Disp => { // todo: check default let mut comp = self.clone(); @@ -266,6 +313,15 @@ impl DFAComponent { self.clone() } + pub fn normalize_rsh(&self) -> Self { + assert!(self.component_type == ComponentType::RSHF); + let mut new_component = self.clone(); + let (omega, alpha, beta) = self.get_rsh(1); + new_component.func = "RSH".to_string(); + new_component.param_positional = vec![omega.unwrap(), alpha, beta]; + new_component.factor = 1.0; + new_component + } pub fn normalize_pt2_param(&self) -> Self { let mut new_component = DFAComponent::new(1.0, self.func.clone()); @@ -328,10 +384,13 @@ impl DFAComponent { pub fn is_normalized(&self) -> bool { match self.component_type { - ComponentType::PT2 | ComponentType::SCSRPA => { + ComponentType::PT2 | ComponentType::SCSRPA | ComponentType::SBGE2 => { self.factor == 1.0 }, - _ => panic!("Error: only PT2 and SCSRPA components can be checked for normalization, but got {}", self.component_type.as_str()), + ComponentType::RSHF => { + self.factor == 1.0 && self.func == "RSH" && self.param_positional.len() == 3 + }, + _ => panic!("Error: only PT2, SCSRPA, SBGE2, and RSHF components can be checked for normalization, but got {}", self.component_type.as_str()), } } @@ -366,6 +425,13 @@ impl Addable for DFAComponent { return false; } }, + ComponentType::RSHF => { + if self.is_normalized() && other.is_normalized() && self.param_positional[0] == other.param_positional[0] { + return true; + } else { + return false; + } + } _ => { return false; }, @@ -397,6 +463,25 @@ impl std::ops::Add for DFAComponent { reference: self.reference.clone(), // todo: check if reference is the same } }, + ComponentType::RSHF => { + DFAComponent { + factor: 1.0, // should be normalized before add + func: self.func.clone(), + func_full_name: self.func_full_name.clone(), + id: self.id, + // direct add positional parameters + param_positional: { + let mut p = self.param_positional.clone(); + p[1] += other.param_positional[1]; + p[2] += other.param_positional[2]; + p + }, + param_keyword: self.param_keyword.clone(), // should be empty + component_type: self.component_type.clone(), + // xcfunc: None, + reference: self.reference.clone(), // todo: check if reference is the same + } + }, _ => { DFAComponent { factor: self.factor + other.factor, @@ -423,7 +508,7 @@ pub struct DFAdef { pub spin_channel: usize, pub dfa_hybrid_scf: f64, // (omega, alpha, beta) in libxc convention; note beta is usually not used in computation. - pub dfa_rsh_scf: Option<(f64, f64, f64)>, + pub dfa_rsh_scf: (Option, f64, f64), pub dfa_hybrid_nscf: Option, } @@ -436,7 +521,7 @@ impl DFAdef { reference: Vec::new(), spin_channel: 1, dfa_hybrid_scf: 0.0, - dfa_rsh_scf: None, + dfa_rsh_scf: (None, 0.0, 0.0), dfa_hybrid_nscf: None, } } @@ -455,7 +540,13 @@ impl DFAdef { } } xc_data.dfa_hybrid_scf = self.dfa_hybrid_scf; - xc_data.dfa_rsh_scf = self.dfa_rsh_scf; + let alpha = self.dfa_rsh_scf.1; + let beta = self.dfa_rsh_scf.2; + if let Some(omega) = self.dfa_rsh_scf.0 { + xc_data.dfa_rsh_scf = Some((omega, alpha, beta)); + } else { + xc_data.dfa_rsh_scf = None; + } if let Some(nscf_components) = &self.xc_nscf { let mut dfa_compnt_pos = Vec::new(); let mut dfa_paramr_pos = Vec::new(); @@ -497,6 +588,9 @@ impl DFAdef { // println!("Info for SCF functional:"); // let hyb_0 = self.get_hybrid_scf(1); result.push_str(&format!(" Total hybrid: {}\n", self.dfa_hybrid_scf)); + if self.is_rsh() { + result.push_str(&format!(" RSH parameters (omega, alpha, beta): {} {} {}\n", self.dfa_rsh_scf.0.unwrap(), self.dfa_rsh_scf.1, self.dfa_rsh_scf.2)); + } if let Some(nscf_components) = &self.xc_nscf { result.push_str("Final energy components:\n"); for comp in nscf_components { @@ -525,22 +619,33 @@ impl DFAdef { } } - pub fn get_rsh_scf(&self, spin_channel: usize) -> Option<(f64, f64, f64)> { + pub fn get_rsh_scf(&self, spin_channel: usize) -> (Option, f64, f64) { // check if any component in self.xc is RSH, if so, return (omega, alpha, beta) - let mut result = None; + let mut result = (None, 0.0, 0.0); if let Some(components) = &self.xc_scf { for comp in components.iter() { if comp.is_rsh() { - let xcfunc = xc_func_init(comp.id, spin_channel); - let (omega, alpha, beta) = xcfunc.cam_coef().unwrap_or((0.0, 0.0, 0.0)); + let (omega, alpha, beta) = comp.get_rsh(spin_channel); + // println!("Component {}: omega = {:?}, alpha = {}, beta = {}", comp.func, omega, alpha, beta); // only if alpha and beta are both close to zero, we consider it as not range-separated (pure zero). if alpha.abs() < 1e-10 && beta.abs() < 1e-10 { continue; } - if result.is_some() { - panic!("Multiple RSH functionals are specified in the DFA components for SCF. Currently this is not supported."); + if result.0.is_some() { + if result.0.unwrap() != omega.unwrap() { + panic!("Multiple omega detected in the DFA components for SCF. Currently this is not supported."); + } + } else { + result.0 = omega; } - result = Some((omega, alpha, beta)); + // result = Some((omega, result.unwrap().1 + alpha, result.unwrap().2 + beta)); + result.1 += alpha; + result.2 += beta; + + } else { + // if it's not RSH but has hybrid, we consider it as alpha with beta = 0. + let hyb = comp.get_hybrid(spin_channel); + result.1 += hyb; } } } diff --git a/src/dft/parse_xc/xc_helper.rs b/src/dft/parse_xc/xc_helper.rs index 3df6f3896c897785f936473ef8c6b5f1ebea2b75..877bfc8f1b9d250db07db8d467cb940850aa2562 100644 --- a/src/dft/parse_xc/xc_helper.rs +++ b/src/dft/parse_xc/xc_helper.rs @@ -6,6 +6,7 @@ use crate::dft::DFAFamily; #[derive(Clone,PartialEq,Debug)] pub enum ComponentType { HF, + RSHF, PT2, RPA, SCSRPA, @@ -27,6 +28,7 @@ impl ComponentType { pub fn as_str(&self) -> &'static str { match self { ComponentType::HF => "HF", + ComponentType::RSHF => "RSHF", ComponentType::PT2 => "PT2", ComponentType::RPA => "RPA", ComponentType::SCSRPA => "SCSRPA", @@ -40,6 +42,7 @@ impl ComponentType { pub fn to_dfa_family(&self) -> Option { match self { ComponentType::HF => None, + ComponentType::RSHF => None, ComponentType::PT2 => Some(DFAFamily::PT2), ComponentType::RPA => Some(DFAFamily::RPA), ComponentType::SCSRPA => Some(DFAFamily::SCSRPA), @@ -197,7 +200,9 @@ lazy_static! { pub static ref WHITELIST_NONLIBXC:HashMap<&'static str, ComponentType> = HashMap::from([ ("HF", ComponentType::HF), - // todo SR_HF + ("SR_HF", ComponentType::RSHF), + ("LR_HF", ComponentType::RSHF), + ("RSH", ComponentType::RSHF), ("MP2", ComponentType::PT2), ("MP2_OS", ComponentType::PT2), ("MP2_SS", ComponentType::PT2), @@ -276,6 +281,8 @@ pub fn load_json_functionals() -> HashMap { map1.extend(load_json(jsonfile2)); let jsonfile3 = include_str!("./family_bdh.json"); map1.extend(load_json(jsonfile3)); + let jsonfile4 = include_str!("./family_bdh_omega.json"); + map1.extend(load_json(jsonfile4)); map1 } diff --git a/tests/rsh/rsh_energy_ri_direct.rs b/tests/rsh/rsh_energy_ri_direct.rs index 745a3707231b825dc044df1900a59f205b3a5198..068c9fea50cdc264bbf1c3e35f10d30998d0b456 100644 --- a/tests/rsh/rsh_energy_ri_direct.rs +++ b/tests/rsh/rsh_energy_ri_direct.rs @@ -8,6 +8,7 @@ static CTRL_TEMPLATE: &str = r##" [ctrl] print_level = 1 xc = "XC_PLACEHOLDER" + xc_parser = "parse_xc" basis_path = "basis-set-pool/sto-3g" auxbas_path = "basis-set-pool/def2-universal-jkfit" eri_type = "ri-v" @@ -50,6 +51,16 @@ fn test_cam_b3lyp_energy() { "CAM-B3LYP energy mismatch: REST={} pyscf={}", scf.scf_energy, pyscf_ref); } +#[test] +fn test_cam_b3lyp_another_energy() { + let scf = run_scf("RSH(0.33,0.65,-0.46) + 0.46*ITYH + 0.35*B88, 0.19*VWN5 + 0.81*LYP"); + let pyscf_ref = -55.315789868819; + println!("REST CAM-B3LYP: {:.12}", scf.scf_energy); + println!("pyscf: {:.12}", pyscf_ref); + assert!((scf.scf_energy - pyscf_ref).abs() < 1e-5, + "CAM-B3LYP energy mismatch: REST={} pyscf={}", scf.scf_energy, pyscf_ref); +} + #[test] fn test_wb97x_energy() { let scf = run_scf("wb97x"); diff --git a/tests/test_dft.rs b/tests/test_dft.rs new file mode 100644 index 0000000000000000000000000000000000000000..3f79b45668f021cb6d3c4175bc11cafb15410a35 --- /dev/null +++ b/tests/test_dft.rs @@ -0,0 +1,385 @@ +use std::collections::HashMap; +use std::io::Read; +use rest_tensors::{MatrixFull, ParMathMatrix}; +use pyrest::basis_io::{BasCell, Basis4Elem, cint_norm_factor, gto_value, }; +use pyrest::dft::{DFA4REST, Grids}; +use pyrest::dft::libxc_helper::{xc_func_init, lda_exc_vxc}; +use pyrest::utilities; +use regex::Regex; + +#[test] +fn debug_num_density_for_atom() { + let angular = 2; + let num_basis = 2*angular+1; + let mut dm = [ + //MatrixFull::from_vec([1,1], vec![2.0]).unwrap(), + //MatrixFull::from_vec([num_basis,num_basis], vec![ + // 1.0,0.0,0.0, + // 0.0,1.0,0.0, + // 0.0,0.0,0.0]).unwrap(), + MatrixFull::from_vec([5,5], vec![ + 2.0,0.0,0.0,0.0,0.0, + 0.0,0.0,0.0,0.0,0.0, + 0.0,0.0,0.0,0.0,0.0, + 0.0,0.0,0.0,0.0,0.0, + 0.0,0.0,0.0,0.0,0.0]).unwrap(), + MatrixFull::empty()]; + let mut alpha_min_h: HashMap = HashMap::new(); + alpha_min_h.insert(angular,0.122); + let alpha_max_h: f64 = 0.122; + let basis4elem = vec![Basis4Elem { + electron_shells: vec![ + BasCell { + function_type: None, + region: None, + angular_momentum: vec![angular as i32], + exponents: vec![0.122], + coefficients: vec![vec![1.0/cint_norm_factor(angular as i32, 0.122)]], + native_coefficients: vec![vec![1.0/cint_norm_factor(angular as i32, 0.122)]] + } + ], + references: None, + ecp_electrons: None, + ecp_potentials: None, + global_index: (0,0) + }]; + let center_coordinates_bohr = vec![(0.0,0.0,0.0)]; + let proton_charges = vec![1]; + let grids = Grids::build_nonstd( + center_coordinates_bohr.clone(), + proton_charges.clone(), + vec![alpha_min_h], + vec![alpha_max_h], &mut None); + let mut total_density = 0.0; + let mut count:usize =0; + grids.coordinates.iter().zip(grids.weights.iter()).for_each(|(r,w)| { + let mut density_r_sum = 0.0; + let mut density_r:Vec = vec![]; + basis4elem.iter().zip(center_coordinates_bohr.iter()).for_each(|(elem,geom_nonstd)| { + let geom = [geom_nonstd.0,geom_nonstd.1,geom_nonstd.2]; + let mut tmp_geom = [0.0;3]; + tmp_geom.iter_mut().zip(geom.iter()).for_each(|value| {*value.0 = *value.1}); + //density_r.extend(gto_value(r, &tmp_geom, elem, &mol.ctrl.basis_type)); + //let tmp_vec = gto_value(r, &tmp_geom, elem, &"spheric".to_string()); + let tmp_vec = gto_value(r, &tmp_geom, elem, &"spheric".to_string()); + //println!("debug 0: len {}", &tmp_vec.len()); + density_r.extend(tmp_vec); + //if count<=10 {println!("{:?},{:?},{:?}", density_r,elem.electron_shells[0].exponents, elem.electron_shells[0].coefficients[0])}; + }); + //println!{"debug 1"}; + let mut density_rr = MatrixFull::from_vec([num_basis,1],density_r).unwrap(); + //println!{"debug 2"}; + dm.iter_mut().for_each(|dm_s| { + let mut tmp_mat = MatrixFull::new([num_basis,1],0.0); + tmp_mat.lapack_dgemm(&mut density_rr, dm_s, 'T', 'N', 1.0, 0.0); + if count<=10 {println!("count: {},{:?}",count, &tmp_mat.data)}; + density_r_sum += tmp_mat.data.iter().zip(density_rr.data.iter()).fold(0.0, |acc,(a,b)| {acc + a*b}); + }); + if count<=10 {println!("{:?},{},{}", r,w,density_r_sum)}; + count += 1; + //println!{"debug 3"}; + total_density += density_r_sum * w; + }); + + println!("Total density: {}", total_density); + +} + +#[test] +fn test_libxc() { + let rho:Vec = vec![0.1,0.2,0.3,0.4,0.5,0.6,0.8]; + let _sigma:Vec = vec![0.2,0.3,0.4,0.5,0.6,0.7]; + //&rho.par_iter().for_each(|c| {println!("{:16.8}",c)}); + //let mut exc:Vec = vec![0.0,0.0,0.0,0.0,0.0]; + //let mut vrho:Vec = vec![0.0,0.0,0.0,0.0,0.0]; + //let func_id: usize = ffi_xc::XC_GGA_X_XPBE as usize; + let spin_channel: usize = 1; + + let my_xc = DFA4REST::parse_scf("lda_x_slater", spin_channel); + //let mut my_xc = XcFuncType::xc_func_init_fdqc(&"pw-lda",spin_channel); + + + + let mut exc = MatrixFull::new([rho.len()/spin_channel,1],0.0); + let mut vrho = MatrixFull::new([rho.len()/spin_channel,spin_channel],0.0); + my_xc.dfa_compnt_scf.iter().zip(my_xc.dfa_paramr_scf.iter()).for_each(|(xc_func, xc_para)| { + //let mut new_vec = vec![-0.34280861230056237, -0.43191178672272906, -0.494415573788165, -0.5441747517896713, -0.586194481347579, -0.622924588811561, -0.6856172246011247]; + //let tmp_c = (new_vec.as_mut_ptr(), new_vec.len(),new_vec.capacity()); + //let new_vec = unsafe{Vec::from_raw_parts(tmp_c.0, tmp_c.1, tmp_c.2)}; + //new_vec.par_iter().for_each(|c| {println!("{:16.8e}",c)}); + let xc_func = my_xc.init_libxc(xc_func); + + let (tmp_exc, tmp_vrho) = lda_exc_vxc(&xc_func,&rho); + //let tmp_exc_2 = tmp_exc.clone(); + //println!("WARNNING:: unsolved rayon par_iter problem. It should be relevant to be the fact that tmp_exc is prepared by libxc via ffi"); + //println!("tmp_vec_2 copied from tmp_exc: {:?},{},{}", &tmp_exc_2, tmp_exc_2.len(),tmp_exc_2.capacity()); + //println!("tmp_vec"); + //&tmp_exc.iter().for_each(|c| { + // println!("{:16.8e}",c); + //}); + //println!("tmp_vec_2"); + //&tmp_exc_2.par_iter().for_each(|c| { + // println!("{:16.8e}",c); + //}); + //println!("tmp_vec_2"); + //&tmp_exc.iter().for_each(|c| { + // println!("{:16.8e}",c); + //}); + let tmp_exc = MatrixFull::from_vec([rho.len()/spin_channel,1],tmp_exc).unwrap(); + let tmp_vrho = MatrixFull::from_vec([rho.len()/spin_channel,spin_channel],tmp_vrho).unwrap(); + //println!("{:?}", &tmp_exc.data); + //println!("{:?}", &tmp_vrho.data); + exc.par_self_scaled_add(&tmp_exc,*xc_para); + vrho.par_self_scaled_add(&tmp_vrho,*xc_para); + //exc.data.par_iter_mut().zip(tmp_exc.data.par_iter()).for_each(|(c,p)| { + // println!("{:16.8e},{:16.8e}",c,p); + //}); + //println!("{:?}", &exc.data); + //println!("{:?}", &vrho.data); + }); + println!("{:?}", exc.data); + println!("{:?}", vrho.data); +} + +#[test] +#[ignore] +#[allow(unused)] +fn read_grid() { + let mut grids_file = std::fs::File::open("/home/igor/Documents/Package-Pool/Rust/rest/grids").unwrap(); + let mut content = String::new(); + grids_file.read_to_string(&mut content); + //println!("{}",&content); + let re1 = Regex::new(r"(?x)\s* + (?P[\+-]?\d+.\d+[eE][\+-]?\d+)\s*,# the 'x' position + \s+ + (?P[\+-]?\d+.\d+[eE][\+-]?\d+)\s*,# the 'y' position + \s+ + (?P[\+-]?\d+.\d+[eE][\+-]?\d+)\s*,# the 'z' position + \s+ + (?P[\+-]?\d+.\d+[eE][\+-]?\d+)\s*# the 'w' weight + \s*\n").unwrap(); + //if let Some(cap) = re1.captures(&content) { + // println!("{:?}", &cap) + //} + for cap in re1.captures_iter(&content) { + let x:f64 = cap[1].parse().unwrap(); + let y:f64 = cap[2].parse().unwrap(); + let z:f64 = cap[3].parse().unwrap(); + let w:f64 = cap[4].parse().unwrap(); + println!("{:16.8} {:16.8} {:16.8} {:16.8}", x,y,z,w); + } +} +#[test] +fn debug_transpose() { + let len_a = 111_usize; + let len_b = 40000_usize; + let orig_a:Vec = (0..len_a*len_b).map(|i| {i as f64}).collect(); + let a_mat = MatrixFull::from_vec([len_a,len_b],orig_a).unwrap(); + let dt0 = utilities::init_timing(); + let b_mat = a_mat.transpose_and_drop(); + //b_mat.formated_output(10, "full"); + let _dt1 = utilities::timing(&dt0, Some("old transpose")); + let orig_a:Vec = (0..len_a*len_b).map(|i| {i as f64}).collect(); + let a_mat = MatrixFull::from_vec([len_a,len_b],orig_a).unwrap(); + let dt0 = utilities::init_timing(); + let c_mat = a_mat.transpose_and_drop(); + //b_mat.formated_output(10, "full"); + let _dt1 = utilities::timing(&dt0, Some("new transpose")); + b_mat.data.iter().zip(c_mat.data.iter()).for_each(|(b,c)| { + assert!(*b==*c); + }); +} + + + + + +#[test] +fn test_rsh_cam_coeff_raw() { + // Test raw libxc CAM coefficients match pyscf expectations + // Each entry: (libxc_id, expected_omega, expected_alpha, expected_beta, expected_hyb) + // hyb = alpha + beta (pyscf convention for RSH) + let ref_data = vec![ + (464, 0.300, 1.0, -0.842294, 0.157706), // WB97X + (466, 0.300, 1.0, -0.833, 0.167), // WB97X-V + (471, 0.200, 1.0, -0.777964, 0.222036), // WB97X-D + (433, 0.330, 0.65, -0.46, 0.19), // CAM-B3LYP + (400, 0.330, 1.0, -1.0, 0.0), // LC-BLYP + (478, 0.400, 1.0, -1.0, 0.0), // LC-wPBE + (428, 0.110, 0.0, 0.25, 0.25), // HSE06 + (427, 0.1061, 0.0, 0.25, 0.25), // HSE03 + ]; + + for (libxc_id, exp_omega, exp_alpha, exp_beta, exp_hyb) in &ref_data { + let func = xc_func_init(*libxc_id, 1); + assert!(func.is_hyb_cam(), "ID {} should be RSH/CAM", libxc_id); + let (omega, alpha, beta) = func.cam_coef().unwrap_or((0.0, 0.0, 0.0)); + let hyb = alpha + beta; + assert!((omega - exp_omega).abs() < 1e-3, "ID {}: omega mismatch: got {} expected {}", libxc_id, omega, exp_omega); + assert!((alpha - exp_alpha).abs() < 1e-3, "ID {}: alpha mismatch: got {} expected {}", libxc_id, alpha, exp_alpha); + assert!((beta - exp_beta).abs() < 1e-3, "ID {}: beta mismatch: got {} expected {}", libxc_id, beta, exp_beta); + assert!((hyb - exp_hyb).abs() < 1e-4, "ID {}: hyb mismatch: got {} expected {}", libxc_id, hyb, exp_hyb); + println!("ID {}: omega={:.4} alpha={:.4} beta={:+.4} hyb={:.6} OK", libxc_id, omega, alpha, beta, hyb); + } + println!("All {} RSH libxc raw CAM coefficient tests passed.", ref_data.len()); +} + +#[test] +fn test_rsh_parameters_wb97x() { + let dfa = DFA4REST::new("wb97x", 1, 0); + assert!(dfa.is_hybrid()); + assert!(dfa.is_rsh()); + assert!((dfa.omega().unwrap() - 0.3).abs() < 1e-9); + assert!((dfa.rsh_alpha().unwrap() - 1.0).abs() < 1e-9); + assert!((dfa.dfa_hybrid_scf - 0.157706).abs() < 1e-5); + assert_eq!(dfa.dfa_compnt_scf, vec![464]); +} + +#[test] +fn test_rsh_parameters_cam_b3lyp() { + let dfa = DFA4REST::new("cam-b3lyp", 1, 0); + assert!(dfa.is_hybrid()); + assert!(dfa.is_rsh()); + assert!((dfa.omega().unwrap() - 0.33).abs() < 1e-9); + assert!((dfa.rsh_alpha().unwrap() - 0.65).abs() < 1e-5); + assert!((dfa.dfa_hybrid_scf - 0.19).abs() < 1e-5); + assert_eq!(dfa.dfa_compnt_scf, vec![433]); + // also test case variation + let dfa2 = DFA4REST::new("CAM-B3LYP", 1, 0); + assert!((dfa.omega().unwrap() - dfa2.omega().unwrap()).abs() < 1e-9); + assert!((dfa.dfa_hybrid_scf - dfa2.dfa_hybrid_scf).abs() < 1e-9); +} + +#[test] +fn test_rsh_parameters_lc_blyp() { + let dfa = DFA4REST::new("lc-blyp", 1, 0); + assert!(dfa.is_hybrid()); + assert!(dfa.is_rsh()); + assert!((dfa.omega().unwrap() - 0.33).abs() < 1e-9); + assert!((dfa.rsh_alpha().unwrap() - 1.0).abs() < 1e-9); + // LC-BLYP has hyb=0 (alpha=1, beta=-1 → hyb=0, LR-only HF) + assert!((dfa.dfa_hybrid_scf - 0.0).abs() < 1e-8); + assert_eq!(dfa.dfa_compnt_scf, vec![400]); +} + +#[test] +fn test_rsh_parameters_hse06() { + let dfa = DFA4REST::new("hse06", 1, 0); + assert!(dfa.is_hybrid()); + assert!(dfa.is_rsh()); + assert!((dfa.omega().unwrap() - 0.11).abs() < 1e-9); + // HSE06 has alpha=0 (no additional LR-HF), beta=0.25 (SR-HF only) + assert!((dfa.rsh_alpha().unwrap() - 0.0).abs() < 1e-5); + assert!((dfa.dfa_hybrid_scf - 0.25).abs() < 1e-5); + assert_eq!(dfa.dfa_compnt_scf, vec![428]); +} + +#[test] +fn test_non_rsh_b3lyp() { + let dfa = DFA4REST::new("b3lyp", 1, 0); + assert!(dfa.is_hybrid(), "B3LYP should be hybrid"); + assert!(!dfa.is_rsh(), "B3LYP should NOT be range-separated"); + assert!(dfa.omega().is_none()); + assert!(dfa.rsh_alpha().is_none()); +} + +#[test] +fn test_non_rsh_pbe() { + let dfa = DFA4REST::new("pbe", 1, 0); + assert!(!dfa.is_hybrid()); + assert!(!dfa.is_rsh()); + assert!(dfa.omega().is_none()); + assert!(dfa.rsh_alpha().is_none()); + assert_eq!(dfa.dfa_hybrid_scf, 0.0); +} + +#[test] +fn test_cam_b3lyp_case_insensitive() { + let names = ["cam-b3lyp", "CAMB3LYP"]; + let mut dfas = vec![]; + for name in &names { + dfas.push(DFA4REST::new(name, 1, 0)); + } + for i in 1..dfas.len() { + assert!((dfas[0].omega().unwrap() - dfas[i].omega().unwrap()).abs() < 1e-9); + assert!((dfas[0].rsh_alpha().unwrap() - dfas[i].rsh_alpha().unwrap()).abs() < 1e-9); + assert!((dfas[0].dfa_hybrid_scf - dfas[i].dfa_hybrid_scf).abs() < 1e-9); + } +} + +#[test] +fn test_rsh_spin_polarized() { + // Test RSH initialization with spin=2 (open-shell) + let dfa_closed = DFA4REST::new("cam-b3lyp", 1, 0); + let dfa_open = DFA4REST::new("cam-b3lyp", 2, 0); + // RSH parameters should be spin-independent + assert!((dfa_closed.omega().unwrap() - dfa_open.omega().unwrap()).abs() < 1e-9); + assert!((dfa_closed.rsh_alpha().unwrap() - dfa_open.rsh_alpha().unwrap()).abs() < 1e-9); + assert!((dfa_closed.dfa_hybrid_scf - dfa_open.dfa_hybrid_scf).abs() < 1e-9); +} + +#[test] +fn test_rsh_hyb_override() { + // Verify that for RSH functionals, dfa_hybrid_scf is overridden to alpha+beta + // NOT the raw xc_hyb_exx_coeff (which may return 1.0) + let dfa = DFA4REST::new("lc-blyp", 1, 0); + // LC-BLYP: raw xc_hyb_exx_coef = 1.0 (because it's a hybrid), + // but alpha+beta = 0.0 (pure LR-HF, no SR-HF) + // The override should set hyb to 0.0 + assert!((dfa.dfa_hybrid_scf - 0.0).abs() < 1e-8, + "LC-BLYP hyb should be 0 (alpha+beta=0), got {}", dfa.dfa_hybrid_scf); + + let dfa2 = DFA4REST::new("cam-b3lyp", 1, 0); + // CAM-B3LYP: raw xc_hyb_exx_coef = ?, alpha+beta = 0.19 + assert!((dfa2.dfa_hybrid_scf - 0.19).abs() < 1e-3, + "CAM-B3LYP hyb should be ~0.19 (alpha+beta), got {}", dfa2.dfa_hybrid_scf); +} + +#[test] +fn test_rsh_nonstd_parse() { + // Test parse_scf_nonstd with RSH components (e.g. custom mixing) + let dfa = DFA4REST::parse_scf_nonstd( + &vec!["HYB_GGA_XC_WB97X_V".to_string()], + &vec![0.5], // 50% scaling + &0.167, // hyb = alpha+beta + 1, + ); + assert!(dfa.is_rsh(), "Nonstd WB97X-V should be RSH"); + assert!((dfa.omega().unwrap() - 0.3).abs() < 1e-9); + assert!((dfa.rsh_alpha().unwrap() - 1.0).abs() < 1e-9); + assert!((dfa.dfa_hybrid_scf - 0.167).abs() < 1e-5); + + // Test with HSE06 nonstd + let dfa2 = DFA4REST::parse_scf_nonstd( + &vec!["HYB_GGA_XC_HSE06".to_string()], + &vec![1.0], + &0.25, + 1, + ); + assert!(dfa2.is_rsh()); + assert!((dfa2.omega().unwrap() - 0.11).abs() < 1e-9); + assert!((dfa2.rsh_alpha().unwrap() - 0.0).abs() < 1e-5); +} + +#[test] +fn test_rsh_summary_does_not_panic() { + // Ensure summary() works without panicking for RSH and non-RSH + for name in &["cam-b3lyp", "lc-blyp", "hse06", "b3lyp", "pbe"] { + let dfa = DFA4REST::new(name, 1, 0); + dfa.summary(); + } +} + +#[test] +fn test_all_rsh_via_auto_resolver() { + // Test that the auto-resolver (generic libxc name lookup) works for RSH + let dfa = DFA4REST::new("HYB_GGA_XC_WB97X_V", 1, 0); + assert!(dfa.is_rsh()); + assert!((dfa.omega().unwrap() - 0.3).abs() < 1e-9); + assert_eq!(dfa.dfa_compnt_scf, vec![466]); + + let dfa2 = DFA4REST::new("HYB_GGA_XC_HSE06", 1, 0); + assert!(dfa2.is_rsh()); + assert!((dfa2.omega().unwrap() - 0.11).abs() < 1e-9); + assert_eq!(dfa2.dfa_compnt_scf, vec![428]); +} \ No newline at end of file diff --git a/tests/test_parse_xc.rs b/tests/test_parse_xc.rs index 6ca4e90118e048f7a5d14a674bc8d8078bf2b4db..3a9921967771f906a7ec5457c73d55724ee57174 100644 --- a/tests/test_parse_xc.rs +++ b/tests/test_parse_xc.rs @@ -1,18 +1,18 @@ use pyrest::dft::parse_xc::parse::{parse, parse_1step, parse_and_derive, merge_components, parse_pass3}; use pyrest::dft::parse_xc::ComponentType; -use assert_float_eq::{assert_f64_near}; + #[test] fn test_parse_xc_param() { let input1 = "0.5*PBE + 0.5*B88, PBE(_beta=0.1)"; let final_components = parse_1step(input1); assert_eq!(final_components.len(), 3); - assert_f64_near!(final_components[0].factor, 0.5, 9); + assert!((final_components[0].factor - 0.5).abs() <= 1e-9); assert_eq!(final_components[0].id, 101); assert_eq!(final_components[1].id, 106); assert_eq!(final_components[2].factor, 1.0); assert_eq!(final_components[2].id, 130); - assert_f64_near!(final_components[2].param_keyword.get("_beta").unwrap().as_f64().unwrap(), 0.1, 9); + assert!((final_components[2].param_keyword.get("_beta").unwrap().as_f64().unwrap() - 0.1).abs() <= 1e-9); } #[test] @@ -28,10 +28,10 @@ fn test_parse_xc_hybrid() { assert_eq!(final_components[4].id, 8); let input2 = "B3LYP"; let final_dfa = parse(input2); - assert_f64_near!(final_dfa.get_hybrid_scf(1), 0.2, 9); + assert!((final_dfa.get_hybrid_scf(1) - 0.2).abs() <= 1e-9); let input3 = "0.2*HF + 0.5*B3LYP"; let final_dfa3 = parse(input3); - assert_f64_near!(final_dfa3.get_hybrid_scf(1), 0.3, 9); + assert!((final_dfa3.get_hybrid_scf(1) - 0.3).abs() <= 1e-9); } #[test] @@ -72,13 +72,13 @@ fn test_parse_xc_merge() { let final_components = parse_1step(input1); let merged = merge_components(final_components); assert_eq!(merged.len(), 2); - assert_f64_near!(merged[0].factor, 1.1, 9); + assert!((merged[0].factor - 1.1).abs() <= 1e-9); assert_eq!(merged[0].id, 106); let input2 = "BLYP + 0.1*LYP"; let final_components2 = parse_1step(input2); let merged2 = merge_components(final_components2); assert_eq!(merged2.len(), 2); - assert_f64_near!(merged2[1].factor, 1.1, 9); + assert!((merged2[1].factor - 1.1).abs() <= 1e-9); } #[test] @@ -128,4 +128,46 @@ fn test_parse_xc_vv10() { let dfa2 = parse_and_derive(input2, 1, 0); let (san2,_) = dfa2.check_sanity(); assert!(!san2); -} \ No newline at end of file +} + +#[test] +fn test_parse_xc_cam_b3lyp() { + let input1 = "CAM-B3LYP"; + let dfa = parse_and_derive(input1, 1, 0); + let (san,_) = dfa.check_sanity(); + assert!(san); + let (omega, alpha, _beta) = dfa.get_rsh_scf(1); + assert!((omega.unwrap() - 0.33).abs() <= 1e-7); + assert!((alpha - 0.65).abs() <= 1e-7); + assert!((dfa.dfa_hybrid_scf - 0.19).abs() <= 1e-7); +} + +#[test] +fn test_parse_xc_rsh() { + let input1 = "0.2*SR_HF(0.2) + 0.3*LR_HF(0.2) + 0.5*B3LYP"; + let dfa = parse_and_derive(input1, 1, 0); + let (san,_) = dfa.check_sanity(); + assert!(san); + let (omega, alpha, beta) = dfa.get_rsh_scf(1); + assert!((omega.unwrap() - 0.2).abs() <= 1e-7); + assert!((alpha - 0.4).abs() <= 1e-7); + assert!(((alpha + beta) - 0.3).abs() <= 1e-7); + assert!((dfa.dfa_hybrid_scf - 0.3).abs() <= 1e-7); +} + +#[test] +fn test_parse_xc_equivalent_rsh() { + let input1 = "RSH(0.33,0.65,-0.46) + 0.46*ITYH + 0.35*B88, 0.19*VWN5 + 0.81*LYP"; + let input2 = "CAM-B3LYP"; + let dfa1 = parse_and_derive(input1, 1, 0); + let dfa2 = parse_and_derive(input2, 1, 0); + assert!(dfa1.get_rsh_scf(1) == dfa2.get_rsh_scf(1)); + assert!((dfa1.dfa_hybrid_scf - dfa2.dfa_hybrid_scf).abs() <= 1e-7); +} + +#[test] +#[should_panic] +fn test_parse_xc_rsh_fail() { + let input1 = "0.2*RSH(0.1, 0.5, -0.1) + 0.5*CAM-B3LYP"; + let _dfa = parse_and_derive(input1, 1, 0); +}