From cf674bbcdbcce86ceb767114e62a9dd4470139ed Mon Sep 17 00:00:00 2001 From: Shirong_Wang Date: Wed, 13 May 2026 16:07:19 +0800 Subject: [PATCH 01/11] rsh in parse_xc --- src/dft/parse_xc/parse.rs | 86 +++++++++++++++++++++++++++++------ src/dft/parse_xc/xc_helper.rs | 4 +- tests/test_parse_xc.rs | 25 ++++++++++ 3 files changed, 100 insertions(+), 15 deletions(-) diff --git a/src/dft/parse_xc/parse.rs b/src/dft/parse_xc/parse.rs index 64a7362..d3affb8 100644 --- a/src/dft/parse_xc/parse.rs +++ b/src/dft/parse_xc/parse.rs @@ -184,17 +184,15 @@ 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).unwrap(); + return self.factor * (alpha + beta); + } if self.component_type == ComponentType::HF { return self.factor; } else if self.component_type == ComponentType::Libxc { let mut xcfunc = XcFuncType::xc_func_init(self.id, spin_channel); - let hybrid_coef = match xcfunc.is_rsh() { - false => xcfunc.get_hybrid(), - true => { - let (omega, alpha, beta) = xcfunc.xc_hyb_cam_coef(); - alpha + beta // for RSH, return alpha + beta (matches pyscf) - }, - }; + let hybrid_coef = xcfunc.get_hybrid(); let hyb = self.factor * hybrid_coef; xcfunc.xc_func_end(); return hyb; @@ -203,8 +201,52 @@ impl DFAComponent { } } + pub fn get_rsh(&self, spin_channel: usize) -> Option<(f64, f64, f64)> { + if !self.is_rsh() { + return None; + } + if self.component_type == ComponentType::Libxc { + let mut xcfunc = XcFuncType::xc_func_init(self.id, spin_channel); + // if xcfunc.is_rsh() { + let (omega, alpha, beta) = xcfunc.xc_hyb_cam_coef(); + xcfunc.xc_func_end(); + return Some((omega, alpha, beta)); + // } else { + // xcfunc.xc_func_end(); + // return None; + // } + } else if self.component_type == ComponentType::HF { + 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() { + "HF" => return Some((0.0, self.factor, 0.0)), + "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, + } + } else { + return None; + } + } + pub fn is_rsh(&self) -> bool { match self.component_type { + ComponentType::HF => { + self.func == "LR_HF" || self.func == "SR_HF" || self.func == "RSH" + }, ComponentType::Libxc => { let mut xcfunc = XcFuncType::xc_func_init(self.id, 1); let is_rsh = xcfunc.is_rsh(); @@ -249,6 +291,11 @@ impl DFAComponent { pub fn canonicalize(&self) -> Self { match self.component_type { + ComponentType::HF => { + if self.is_rsh() { + return self.normalize_rsh(); + } + } ComponentType::Disp => { // todo: check default let mut comp = self.clone(); @@ -274,6 +321,16 @@ impl DFAComponent { self.clone() } + pub fn normalize_rsh(&self) -> Self { + assert!(self.component_type == ComponentType::HF); + assert!(self.is_rsh()); + let mut new_component = self.clone(); + let (omega, alpha, beta) = self.get_rsh(1).unwrap(); + new_component.func = "RSH".to_string(); + new_component.param_positional = vec![omega, 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()); @@ -547,21 +604,22 @@ impl DFAdef { pub fn get_rsh_scf(&self, spin_channel: usize) -> Option<(f64, f64, f64)> { // check if any component in self.xc is RSH, if so, return (omega, alpha, beta) - let mut result = None; + let mut result: Option<(f64, f64, f64)> = None; if let Some(components) = &self.xc_scf { for comp in components.iter() { - if comp.is_rsh() { - let mut xcfunc = XcFuncType::xc_func_init(comp.id, spin_channel); - let (omega, alpha, beta) = xcfunc.xc_hyb_cam_coef(); + if let Some((omega, alpha, beta)) = comp.get_rsh(spin_channel) { // 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.unwrap().0 != omega { + panic!("Multiple omega detected in the DFA components for SCF. Currently this is not supported."); + } + result = Some((omega, result.unwrap().1 + alpha, result.unwrap().2 + beta)); + } else { + result = Some((omega, alpha, beta)); } - xcfunc.xc_func_end(); - result = Some((omega, alpha, beta)); } } } diff --git a/src/dft/parse_xc/xc_helper.rs b/src/dft/parse_xc/xc_helper.rs index 634b90f..fb7e24a 100644 --- a/src/dft/parse_xc/xc_helper.rs +++ b/src/dft/parse_xc/xc_helper.rs @@ -199,7 +199,9 @@ lazy_static! { pub static ref WHITELIST_NONLIBXC:HashMap<&'static str, ComponentType> = HashMap::from([ ("HF", ComponentType::HF), - // todo SR_HF + ("SR_HF", ComponentType::HF), + ("LR_HF", ComponentType::HF), + ("RSH", ComponentType::HF), ("MP2", ComponentType::PT2), ("MP2_OS", ComponentType::PT2), ("MP2_SS", ComponentType::PT2), diff --git a/tests/test_parse_xc.rs b/tests/test_parse_xc.rs index 6ca4e90..f024bc3 100644 --- a/tests/test_parse_xc.rs +++ b/tests/test_parse_xc.rs @@ -128,4 +128,29 @@ fn test_parse_xc_vv10() { let dfa2 = parse_and_derive(input2, 1, 0); let (san2,_) = dfa2.check_sanity(); assert!(!san2); +} + +#[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).unwrap(); + assert_f64_near!(omega, 0.33, 7); + assert_f64_near!(alpha, 0.65, 7); + assert_f64_near!(dfa.dfa_hybrid_scf, 0.19, 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).unwrap(); + assert_f64_near!(omega, 0.2, 7); + assert_f64_near!(alpha, 0.3, 7); + assert_f64_near!(alpha+beta, 0.3, 7); + assert_f64_near!(dfa.dfa_hybrid_scf, 0.3, 7); } \ No newline at end of file -- Gitee From 02306518540cce3e948615367f7e166be48df5fa Mon Sep 17 00:00:00 2001 From: Shirong_Wang Date: Wed, 13 May 2026 17:16:23 +0800 Subject: [PATCH 02/11] update rsh logic --- Cargo.toml | 1 - src/dft/parse_xc/parse.rs | 58 ++++++++++++++++++++++++--------------- tests/test_parse_xc.rs | 34 +++++++++++------------ 3 files changed, 53 insertions(+), 40 deletions(-) diff --git a/Cargo.toml b/Cargo.toml index ecbe752..db0da81 100644 --- a/Cargo.toml +++ b/Cargo.toml @@ -62,7 +62,6 @@ dftd4 = { version = "0.2", optional = true } geometric-pyo3 = { version = "0.1.1", optional = true } approx = "0.5" -assert_float_eq = "1" [build-dependencies] autocxx-build = "0.26.0" diff --git a/src/dft/parse_xc/parse.rs b/src/dft/parse_xc/parse.rs index d3affb8..415ceae 100644 --- a/src/dft/parse_xc/parse.rs +++ b/src/dft/parse_xc/parse.rs @@ -185,7 +185,7 @@ 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).unwrap(); + let (omega, alpha, beta) = self.get_rsh(spin_channel); return self.factor * (alpha + beta); } if self.component_type == ComponentType::HF { @@ -201,16 +201,16 @@ impl DFAComponent { } } - pub fn get_rsh(&self, spin_channel: usize) -> Option<(f64, f64, f64)> { + pub fn get_rsh(&self, spin_channel: usize) -> (Option, f64, f64) { if !self.is_rsh() { - return None; + return (None, self.get_hybrid(spin_channel), 0.0); } if self.component_type == ComponentType::Libxc { let mut xcfunc = XcFuncType::xc_func_init(self.id, spin_channel); // if xcfunc.is_rsh() { let (omega, alpha, beta) = xcfunc.xc_hyb_cam_coef(); xcfunc.xc_func_end(); - return Some((omega, alpha, beta)); + return (Some(omega), alpha, beta); // } else { // xcfunc.xc_func_end(); // return None; @@ -231,14 +231,13 @@ impl DFAComponent { _ => {}, } match self.func.as_str() { - "HF" => return Some((0.0, self.factor, 0.0)), - "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, + "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; + return (None, 0.0, 0.0); } } @@ -325,9 +324,9 @@ impl DFAComponent { assert!(self.component_type == ComponentType::HF); assert!(self.is_rsh()); let mut new_component = self.clone(); - let (omega, alpha, beta) = self.get_rsh(1).unwrap(); + let (omega, alpha, beta) = self.get_rsh(1); new_component.func = "RSH".to_string(); - new_component.param_positional = vec![omega, alpha, beta]; + new_component.param_positional = vec![omega.unwrap(), alpha, beta]; new_component.factor = 1.0; new_component } @@ -500,7 +499,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, } @@ -513,7 +512,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, } } @@ -532,7 +531,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(); @@ -602,24 +607,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: Option<(f64, f64, f64)> = None; + let mut result = (None, 0.0, 0.0); if let Some(components) = &self.xc_scf { for comp in components.iter() { - if let Some((omega, alpha, beta)) = comp.get_rsh(spin_channel) { + if comp.is_rsh() { + 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() { - if result.unwrap().0 != omega { + 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."); } - result = Some((omega, result.unwrap().1 + alpha, result.unwrap().2 + beta)); } else { - result = Some((omega, alpha, beta)); + result.0 = omega; } + // 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/tests/test_parse_xc.rs b/tests/test_parse_xc.rs index f024bc3..2e580fe 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] @@ -136,10 +136,10 @@ fn test_parse_xc_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).unwrap(); - assert_f64_near!(omega, 0.33, 7); - assert_f64_near!(alpha, 0.65, 7); - assert_f64_near!(dfa.dfa_hybrid_scf, 0.19, 7); + 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] @@ -148,9 +148,9 @@ fn test_parse_xc_rsh() { let dfa = parse_and_derive(input1, 1, 0); let (san,_) = dfa.check_sanity(); assert!(san); - let (omega, alpha, beta) = dfa.get_rsh_scf(1).unwrap(); - assert_f64_near!(omega, 0.2, 7); - assert_f64_near!(alpha, 0.3, 7); - assert_f64_near!(alpha+beta, 0.3, 7); - assert_f64_near!(dfa.dfa_hybrid_scf, 0.3, 7); -} \ No newline at end of file + 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); +} -- Gitee From 0fd262e8c60be1c5cd4e8eea7f347357300b1d38 Mon Sep 17 00:00:00 2001 From: Shirong_Wang Date: Wed, 13 May 2026 17:46:56 +0800 Subject: [PATCH 03/11] add type RSHF --- src/dft/parse_xc/parse.rs | 48 +++++++++++++++++++++++++++-------- src/dft/parse_xc/xc_helper.rs | 9 ++++--- tests/test_parse_xc.rs | 17 ++++++++++++- 3 files changed, 59 insertions(+), 15 deletions(-) diff --git a/src/dft/parse_xc/parse.rs b/src/dft/parse_xc/parse.rs index 415ceae..6e496ab 100644 --- a/src/dft/parse_xc/parse.rs +++ b/src/dft/parse_xc/parse.rs @@ -215,7 +215,7 @@ impl DFAComponent { // xcfunc.xc_func_end(); // return None; // } - } else if self.component_type == ComponentType::HF { + } else if self.component_type == ComponentType::RSHF { let param = &self.param_positional; match self.func.as_str() { "LR_HF" | "SR_HF" => { @@ -243,8 +243,8 @@ impl DFAComponent { pub fn is_rsh(&self) -> bool { match self.component_type { - ComponentType::HF => { - self.func == "LR_HF" || self.func == "SR_HF" || self.func == "RSH" + ComponentType::RSHF => { + true }, ComponentType::Libxc => { let mut xcfunc = XcFuncType::xc_func_init(self.id, 1); @@ -290,10 +290,8 @@ impl DFAComponent { pub fn canonicalize(&self) -> Self { match self.component_type { - ComponentType::HF => { - if self.is_rsh() { - return self.normalize_rsh(); - } + ComponentType::RSHF => { + return self.normalize_rsh(); } ComponentType::Disp => { // todo: check default @@ -321,8 +319,7 @@ impl DFAComponent { } pub fn normalize_rsh(&self) -> Self { - assert!(self.component_type == ComponentType::HF); - assert!(self.is_rsh()); + 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(); @@ -392,10 +389,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 + }, + _ => panic!("Error: only PT2, SCSRPA, SBGE2, and RSHF components can be checked for normalization, but got {}", self.component_type.as_str()), } } @@ -442,6 +442,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; }, @@ -473,6 +480,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, diff --git a/src/dft/parse_xc/xc_helper.rs b/src/dft/parse_xc/xc_helper.rs index fb7e24a..84ad1a1 100644 --- a/src/dft/parse_xc/xc_helper.rs +++ b/src/dft/parse_xc/xc_helper.rs @@ -8,6 +8,7 @@ use crate::dft::DFAFamily; #[derive(Clone,PartialEq,Debug)] pub enum ComponentType { HF, + RSHF, PT2, RPA, SCSRPA, @@ -29,6 +30,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", @@ -42,6 +44,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), @@ -199,9 +202,9 @@ lazy_static! { pub static ref WHITELIST_NONLIBXC:HashMap<&'static str, ComponentType> = HashMap::from([ ("HF", ComponentType::HF), - ("SR_HF", ComponentType::HF), - ("LR_HF", ComponentType::HF), - ("RSH", ComponentType::HF), + ("SR_HF", ComponentType::RSHF), + ("LR_HF", ComponentType::RSHF), + ("RSH", ComponentType::RSHF), ("MP2", ComponentType::PT2), ("MP2_OS", ComponentType::PT2), ("MP2_SS", ComponentType::PT2), diff --git a/tests/test_parse_xc.rs b/tests/test_parse_xc.rs index 2e580fe..95c280f 100644 --- a/tests/test_parse_xc.rs +++ b/tests/test_parse_xc.rs @@ -136,7 +136,7 @@ fn test_parse_xc_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); + 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); @@ -153,4 +153,19 @@ fn test_parse_xc_rsh() { 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); + let input2 = "0.1*LYP + CAM-B3LYP"; + let dfa2 = parse_and_derive(input2, 1, 0); + let (san2,_) = dfa2.check_sanity(); + assert!(san2); + let (omega2, alpha2, _beta2) = dfa2.get_rsh_scf(1); + assert!((omega2.unwrap() - 0.33).abs() <= 1e-7); + assert!((alpha2 - 0.65).abs() <= 1e-7); + assert!((dfa2.dfa_hybrid_scf - 0.19).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); } -- Gitee From ce590a6b5a38fbece4530d623cb8ac7ea7461c81 Mon Sep 17 00:00:00 2001 From: Shirong_Wang Date: Wed, 13 May 2026 19:31:04 +0800 Subject: [PATCH 04/11] update test --- src/dft/parse_xc/parse.rs | 3 +++ tests/test_parse_xc.rs | 16 +++++++++------- 2 files changed, 12 insertions(+), 7 deletions(-) diff --git a/src/dft/parse_xc/parse.rs b/src/dft/parse_xc/parse.rs index 6e496ab..9a3fd77 100644 --- a/src/dft/parse_xc/parse.rs +++ b/src/dft/parse_xc/parse.rs @@ -605,6 +605,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 { diff --git a/tests/test_parse_xc.rs b/tests/test_parse_xc.rs index 95c280f..3a99219 100644 --- a/tests/test_parse_xc.rs +++ b/tests/test_parse_xc.rs @@ -153,14 +153,16 @@ fn test_parse_xc_rsh() { 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); - let input2 = "0.1*LYP + CAM-B3LYP"; +} + +#[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); - let (san2,_) = dfa2.check_sanity(); - assert!(san2); - let (omega2, alpha2, _beta2) = dfa2.get_rsh_scf(1); - assert!((omega2.unwrap() - 0.33).abs() <= 1e-7); - assert!((alpha2 - 0.65).abs() <= 1e-7); - assert!((dfa2.dfa_hybrid_scf - 0.19).abs() <= 1e-7); + assert!(dfa1.get_rsh_scf(1) == dfa2.get_rsh_scf(1)); + assert!((dfa1.dfa_hybrid_scf - dfa2.dfa_hybrid_scf).abs() <= 1e-7); } #[test] -- Gitee From 7e8122019996d26d67ad7d7216911e5305ccd013 Mon Sep 17 00:00:00 2001 From: Shirong_Wang Date: Mon, 18 May 2026 15:09:51 +0800 Subject: [PATCH 05/11] update to new api --- src/dft/parse_xc/parse.rs | 12 +++++------- 1 file changed, 5 insertions(+), 7 deletions(-) diff --git a/src/dft/parse_xc/parse.rs b/src/dft/parse_xc/parse.rs index d84c6cc..59f9393 100644 --- a/src/dft/parse_xc/parse.rs +++ b/src/dft/parse_xc/parse.rs @@ -210,15 +210,13 @@ impl DFAComponent { return (None, self.get_hybrid(spin_channel), 0.0); } if self.component_type == ComponentType::Libxc { - let mut xcfunc = XcFuncType::xc_func_init(self.id, spin_channel); + let xcfunc = xc_func_init(self.id, spin_channel); // if xcfunc.is_rsh() { - let (omega, alpha, beta) = xcfunc.xc_hyb_cam_coef(); - xcfunc.xc_func_end(); + if let Some((omega, alpha, beta)) = xcfunc.cam_coef() { return (Some(omega), alpha, beta); - // } else { - // xcfunc.xc_func_end(); - // return None; - // } + } else { + return (None, 0.0, 0.0); + } } else if self.component_type == ComponentType::RSHF { let param = &self.param_positional; match self.func.as_str() { -- Gitee From 1959071335d75840f53dc79de815ec51d18067b6 Mon Sep 17 00:00:00 2001 From: Shirong_Wang Date: Mon, 18 May 2026 16:01:25 +0800 Subject: [PATCH 06/11] set omega, and somw wrap --- src/dft/mod.rs | 176 ++++++++++++++++++++++++++----------------------- 1 file changed, 94 insertions(+), 82 deletions(-) diff --git a/src/dft/mod.rs b/src/dft/mod.rs index ed702c5..3807436 100644 --- a/src/dft/mod.rs +++ b/src/dft/mod.rs @@ -303,6 +303,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 +1200,42 @@ 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) + } + pub fn use_kinetic_density(&self) -> bool { self .dfa_compnt_pos @@ -1207,21 +1253,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 +1270,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 { @@ -1489,7 +1529,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 { @@ -1824,7 +1864,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 { @@ -2088,14 +2128,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,11 +2136,15 @@ 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]; @@ -2196,13 +2232,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 +2249,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,23 +2334,19 @@ 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]; @@ -2350,29 +2377,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,7 +2393,7 @@ 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); + let xc_func = self.init_libxc_and_set_param(xc_func); match xc_func.family() { LibXCFamily::LDA => { let tmp_exc = MatrixFull::from_vec([num_grids,1], @@ -2420,7 +2432,7 @@ impl DFA4REST { } 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); + let xc_func = self.init_libxc_and_set_param(xc_func); match xc_func.family() { LibXCFamily::LDA => { let tmp_exc = MatrixFull::from_vec([num_grids,1], @@ -2478,7 +2490,7 @@ 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 => { -- Gitee From 748abf72e824fcfc25064cf2ea790aad024839a6 Mon Sep 17 00:00:00 2001 From: Shirong_Wang Date: Mon, 18 May 2026 16:15:21 +0800 Subject: [PATCH 07/11] more wrap --- src/dft/mod.rs | 253 ++++++++++++++++++++++++------------------------- 1 file changed, 126 insertions(+), 127 deletions(-) diff --git a/src/dft/mod.rs b/src/dft/mod.rs index 3807436..56ce5f2 100644 --- a/src/dft/mod.rs +++ b/src/dft/mod.rs @@ -1236,6 +1236,84 @@ impl DFA4REST { (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_energy( + &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 @@ -1430,17 +1508,8 @@ 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_energy(&exc, &rho, &grids.weights, spin_channel); if print_level > 0 { if spin_channel==1 { println!("total electron number: {:16.8}", total_elec[0]) @@ -1754,17 +1823,8 @@ 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_energy(&loc_exc, &loc_rho, loc_weights, spin_channel); //if let Some(id) = rayon::current_thread_index() { //} @@ -2091,17 +2151,8 @@ 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_energy(&loc_exc, &loc_rho, loc_weights, spin_channel); //if let Some(id) = rayon::current_thread_index() { //} @@ -2147,19 +2198,19 @@ impl DFA4REST { //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_energy(&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); @@ -2349,19 +2400,19 @@ impl DFA4REST { //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 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 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_energy(&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); @@ -2396,35 +2447,19 @@ impl DFA4REST { let xc_func = self.init_libxc_and_set_param(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); + 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); }, 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); + 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); }, 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); + 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); }, xc_family => panic!("{xc_family:?} is not yet implemented"), } @@ -2435,41 +2470,22 @@ impl DFA4REST { let xc_func = self.init_libxc_and_set_param(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); - }, + 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); + }, + LibXCFamily::GGA | LibXCFamily::HybGGA => { + 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); + }, xc_family => panic!("{xc_family:?} is not yet implemented") } }); } } - 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_energy(&exc, &rho, &grids.weights, spin_channel); //if spin_channel==1 { // println!("total electron number: {:16.8}", total_elec[0]) //} else { @@ -2492,31 +2508,14 @@ 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_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(), + ) } } -- Gitee From f25b76546c8d0a0ea4a876dcba6ff20a3838b979 Mon Sep 17 00:00:00 2001 From: Shirong_Wang Date: Mon, 18 May 2026 16:18:07 +0800 Subject: [PATCH 08/11] reformat --- src/dft/mod.rs | 118 ++++++++++++------------------------------------- 1 file changed, 29 insertions(+), 89 deletions(-) diff --git a/src/dft/mod.rs b/src/dft/mod.rs index 56ce5f2..9389ff1 100644 --- a/src/dft/mod.rs +++ b/src/dft/mod.rs @@ -1236,71 +1236,29 @@ impl DFA4REST { (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 { + 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(), + 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_energy( - &self, - exc: &MatrixFull, - rho: &MatrixFull, - weights: &[f64], - spin_channel: usize, - ) -> (Vec, [f64; 2]) { + fn integrate_exc_energy(&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 { @@ -1508,8 +1466,7 @@ impl DFA4REST { let dt5 = utilities::timing(&dt4, Some("from vrho -> vxc_ao")); - let (exc_total, total_elec) = - self.integrate_exc_energy(&exc, &rho, &grids.weights, spin_channel); + let (exc_total, total_elec) = self.integrate_exc_energy(&exc, &rho, &grids.weights, spin_channel); if print_level > 0 { if spin_channel==1 { println!("total electron number: {:16.8}", total_elec[0]) @@ -1823,8 +1780,7 @@ impl DFA4REST { // println!("{:16.8},{:16.8},{:16.8}", vsigma[[i,0]],vsigma[[i,1]],vsigma[[i,2]]); //}); - let (loc_exc_total, loc_total_elec) = - self.integrate_exc_energy(&loc_exc, &loc_rho, loc_weights, spin_channel); + let (loc_exc_total, loc_total_elec) = self.integrate_exc_energy(&loc_exc, &loc_rho, loc_weights, spin_channel); //if let Some(id) = rayon::current_thread_index() { //} @@ -2151,8 +2107,7 @@ impl DFA4REST { // println!("{:16.8},{:16.8},{:16.8}", vsigma[[i,0]],vsigma[[i,1]],vsigma[[i,2]]); //}); - let (loc_exc_total, loc_total_elec) = - self.integrate_exc_energy(&loc_exc, &loc_rho, loc_weights, spin_channel); + let (loc_exc_total, loc_total_elec) = self.integrate_exc_energy(&loc_exc, &loc_rho, loc_weights, spin_channel); //if let Some(id) = rayon::current_thread_index() { //} @@ -2205,8 +2160,7 @@ impl DFA4REST { exc.par_self_scaled_add(&self.xc_exc_code(xc_code, &rho, &sigma, spin_channel),1.0); }); - let (exc_total_vec, _) = - self.integrate_exc_energy(&exc, &rho, &grids.weights, spin_channel); + let (exc_total_vec, _) = self.integrate_exc_energy(&exc, &rho, &grids.weights, spin_channel); exc_total[0] = exc_total_vec[0]; if spin_channel == 2 { exc_total[1] = exc_total_vec[1]; @@ -2407,8 +2361,7 @@ impl DFA4REST { let exc = self.xc_exc_code(xc_code, &rho, &sigma, spin_channel); //}); - let (exc_total_vec, _) = - self.integrate_exc_energy(&exc, &rho, &grids.weights, spin_channel); + let (exc_total_vec, _) = self.integrate_exc_energy(&exc, &rho, &grids.weights, spin_channel); exc_total[0] = exc_total_vec[0]; if spin_channel == 2 { exc_total[1] = exc_total_vec[1]; @@ -2447,18 +2400,15 @@ impl DFA4REST { let xc_func = self.init_libxc_and_set_param(xc_func); match xc_func.family() { LibXCFamily::LDA => { - let tmp_exc = - self.compute_exc_by_family(&xc_func, spin_channel, &rho, &sigma, &lapl, &tau); + 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); }, LibXCFamily::GGA | LibXCFamily::HybGGA => { - let tmp_exc = - self.compute_exc_by_family(&xc_func, spin_channel, &rho, &sigma, &lapl, &tau); + 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); }, LibXCFamily::MGGA | LibXCFamily::HybMGGA => { - let tmp_exc = - self.compute_exc_by_family(&xc_func, spin_channel, &rho, &sigma, &lapl, &tau); + 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); }, xc_family => panic!("{xc_family:?} is not yet implemented"), @@ -2470,13 +2420,11 @@ impl DFA4REST { let xc_func = self.init_libxc_and_set_param(xc_func); match xc_func.family() { LibXCFamily::LDA => { - let tmp_exc = - self.compute_exc_by_family(&xc_func, spin_channel, &rho, &sigma, &lapl, &tau); + 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); }, LibXCFamily::GGA | LibXCFamily::HybGGA => { - let tmp_exc = - self.compute_exc_by_family(&xc_func, spin_channel, &rho, &sigma, &lapl, &tau); + 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); }, xc_family => panic!("{xc_family:?} is not yet implemented") @@ -2484,8 +2432,7 @@ impl DFA4REST { }); } } - let (exc_total, _total_elec) = - self.integrate_exc_energy(&exc, &rho, &grids.weights, spin_channel); + let (exc_total, _total_elec) = self.integrate_exc_energy(&exc, &rho, &grids.weights, spin_channel); //if spin_channel==1 { // println!("total electron number: {:16.8}", total_elec[0]) //} else { @@ -2508,14 +2455,7 @@ 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_and_set_param(xc_code); let num_grids = rho.size()[0]; - self.compute_exc_by_family( - &xc_func, - spin_channel, - rho, - sigma, - &MatrixFull::empty(), - &MatrixFull::empty(), - ) + self.compute_exc_by_family(&xc_func, spin_channel, rho, sigma, &MatrixFull::empty(), &MatrixFull::empty()) } } -- Gitee From 25951198666da67fa263aef99f9d21fc3f71724b Mon Sep 17 00:00:00 2001 From: Shirong_Wang Date: Mon, 18 May 2026 16:29:55 +0800 Subject: [PATCH 09/11] simplify --- src/dft/mod.rs | 44 +++++++++++--------------------------------- 1 file changed, 11 insertions(+), 33 deletions(-) diff --git a/src/dft/mod.rs b/src/dft/mod.rs index 9389ff1..e4646a9 100644 --- a/src/dft/mod.rs +++ b/src/dft/mod.rs @@ -1258,7 +1258,7 @@ impl DFA4REST { } } - fn integrate_exc_energy(&self, exc: &MatrixFull, rho: &MatrixFull, weights: &[f64], spin_channel: usize) -> (Vec, [f64; 2]) { + 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 { @@ -1466,7 +1466,7 @@ impl DFA4REST { let dt5 = utilities::timing(&dt4, Some("from vrho -> vxc_ao")); - let (exc_total, total_elec) = self.integrate_exc_energy(&exc, &rho, &grids.weights, spin_channel); + 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]) @@ -1780,7 +1780,7 @@ impl DFA4REST { // println!("{:16.8},{:16.8},{:16.8}", vsigma[[i,0]],vsigma[[i,1]],vsigma[[i,2]]); //}); - let (loc_exc_total, loc_total_elec) = self.integrate_exc_energy(&loc_exc, &loc_rho, loc_weights, spin_channel); + 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() { //} @@ -2107,7 +2107,7 @@ impl DFA4REST { // println!("{:16.8},{:16.8},{:16.8}", vsigma[[i,0]],vsigma[[i,1]],vsigma[[i,2]]); //}); - let (loc_exc_total, loc_total_elec) = self.integrate_exc_energy(&loc_exc, &loc_rho, loc_weights, spin_channel); + 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() { //} @@ -2160,7 +2160,7 @@ impl DFA4REST { exc.par_self_scaled_add(&self.xc_exc_code(xc_code, &rho, &sigma, spin_channel),1.0); }); - let (exc_total_vec, _) = self.integrate_exc_energy(&exc, &rho, &grids.weights, spin_channel); + 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]; @@ -2361,7 +2361,7 @@ impl DFA4REST { let exc = self.xc_exc_code(xc_code, &rho, &sigma, spin_channel); //}); - let (exc_total_vec, _) = self.integrate_exc_energy(&exc, &rho, &grids.weights, spin_channel); + 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]; @@ -2398,41 +2398,19 @@ 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_and_set_param(xc_func); - match xc_func.family() { - LibXCFamily::LDA => { - 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); - }, - LibXCFamily::GGA | LibXCFamily::HybGGA => { - 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); - }, - LibXCFamily::MGGA | LibXCFamily::HybMGGA => { - 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); - }, - xc_family => panic!("{xc_family:?} is not yet implemented"), - } + 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_and_set_param(xc_func); - match xc_func.family() { - LibXCFamily::LDA => { - 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); - }, - LibXCFamily::GGA | LibXCFamily::HybGGA => { - 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); - }, - xc_family => panic!("{xc_family:?} is not yet implemented") - } + 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 (exc_total, _total_elec) = self.integrate_exc_energy(&exc, &rho, &grids.weights, spin_channel); + let (exc_total, _total_elec) = self.integrate_exc(&exc, &rho, &grids.weights, spin_channel); //if spin_channel==1 { // println!("total electron number: {:16.8}", total_elec[0]) //} else { -- Gitee From 04d92ab82582c0453875105b9ab7bcf1e307a20e Mon Sep 17 00:00:00 2001 From: Shirong_Wang Date: Mon, 18 May 2026 17:13:21 +0800 Subject: [PATCH 10/11] refine, add wbDH --- src/dft/parse_xc/family_bdh_omega.json | 26 ++++++++++++++++++++++++++ src/dft/parse_xc/parse.rs | 5 ++--- src/dft/parse_xc/xc_helper.rs | 2 ++ 3 files changed, 30 insertions(+), 3 deletions(-) create mode 100644 src/dft/parse_xc/family_bdh_omega.json 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 0000000..dc6e913 --- /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 59f9393..ac096ad 100644 --- a/src/dft/parse_xc/parse.rs +++ b/src/dft/parse_xc/parse.rs @@ -211,11 +211,10 @@ impl DFAComponent { } if self.component_type == ComponentType::Libxc { let xcfunc = xc_func_init(self.id, spin_channel); - // if xcfunc.is_rsh() { if let Some((omega, alpha, beta)) = xcfunc.cam_coef() { return (Some(omega), alpha, beta); } else { - return (None, 0.0, 0.0); + 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; @@ -389,7 +388,7 @@ impl DFAComponent { self.factor == 1.0 }, ComponentType::RSHF => { - self.factor == 1.0 + 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()), } diff --git a/src/dft/parse_xc/xc_helper.rs b/src/dft/parse_xc/xc_helper.rs index deb7831..877bfc8 100644 --- a/src/dft/parse_xc/xc_helper.rs +++ b/src/dft/parse_xc/xc_helper.rs @@ -281,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 } -- Gitee From e045d9cc2961a6334ae8e2f4b681357deda0b143 Mon Sep 17 00:00:00 2001 From: Shirong_Wang Date: Mon, 18 May 2026 20:03:26 +0800 Subject: [PATCH 11/11] move tests --- src/dft/mod.rs | 388 +----------------------------- tests/rsh/rsh_energy_ri_direct.rs | 11 + tests/test_dft.rs | 385 +++++++++++++++++++++++++++++ 3 files changed, 397 insertions(+), 387 deletions(-) create mode 100644 tests/test_dft.rs diff --git a/src/dft/mod.rs b/src/dft/mod.rs index e4646a9..e99db17 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 { @@ -2355,9 +2352,6 @@ impl DFA4REST { 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 exc = self.xc_exc_code(xc_code, &rho, &sigma, spin_channel); //}); @@ -2411,12 +2405,6 @@ impl DFA4REST { } } let (exc_total, _total_elec) = self.integrate_exc(&exc, &rho, &grids.weights, spin_channel); - //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 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()); @@ -3792,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/tests/rsh/rsh_energy_ri_direct.rs b/tests/rsh/rsh_energy_ri_direct.rs index 745a370..068c9fe 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 0000000..3f79b45 --- /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 -- Gitee