From 708d863e8ae3fd0b04cc298c62f45962cc77a134 Mon Sep 17 00:00:00 2001 From: Shirong_Wang Date: Sat, 9 May 2026 22:39:30 +0800 Subject: [PATCH 01/14] save and load basis --- src/basis_io/mod.rs | 4 +- src/fileop/chkfile.rs | 284 ++++++++++++++++++ src/fileop/mod.rs | 1 + src/initial_guess/chkfile.rs | 17 -- src/initial_guess/enxc.rs | 9 +- src/initial_guess/mod.rs | 168 +++++++---- src/initial_guess/proj.rs | 6 + src/initial_guess/sad.rs | 5 +- src/initial_guess/sap.rs | 4 +- src/lib.rs | 1 + src/molecule_io/mod.rs | 8 + src/post_scf_analysis/cube_build.rs | 3 +- src/post_scf_analysis/mod.rs | 181 +---------- src/post_scf_analysis/molden_build.rs | 8 +- src/post_scf_analysis/mulliken.rs | 6 +- src/post_scf_analysis/rand_wf_real_space.rs | 3 +- src/post_scf_analysis/rrs_pbc.rs | 2 - src/post_scf_analysis/spin_correction.rs | 2 +- .../strong_correlation_correction.rs | 4 +- 19 files changed, 432 insertions(+), 284 deletions(-) create mode 100644 src/fileop/chkfile.rs create mode 100644 src/fileop/mod.rs delete mode 100644 src/initial_guess/chkfile.rs create mode 100644 src/initial_guess/proj.rs diff --git a/src/basis_io/mod.rs b/src/basis_io/mod.rs index 0928983..173f138 100644 --- a/src/basis_io/mod.rs +++ b/src/basis_io/mod.rs @@ -80,7 +80,7 @@ fn import_ecp()-> anyhow::Result<()> { /// - BasCell.angular_momentum: The angular momentums of the GTOs in this cell /// - BasCell.exponents: The exponents /// - BasCell.coefficients: The coefficients -#[derive(Debug,Clone)] +#[derive(Debug, Clone, Serialize, Deserialize)] pub struct BasCell { pub function_type: Option, pub region: Option, @@ -137,7 +137,7 @@ impl BasInfo { /// self.electron_shells: the GTO basis functions organized cell by cell [`BasCell`](BasCell) /// self.references: The reference of each basis cell. /// self.global_index: The global index for the given atom -#[derive(Clone, Debug)] +#[derive(Clone, Debug, Serialize, Deserialize)] pub struct Basis4Elem { pub electron_shells: Vec, pub references: Option>, diff --git a/src/fileop/chkfile.rs b/src/fileop/chkfile.rs new file mode 100644 index 0000000..48c0e50 --- /dev/null +++ b/src/fileop/chkfile.rs @@ -0,0 +1,284 @@ +use std::collections::HashMap; +use std::path::Path; +use hdf5; +use hdf5::types::VarLenUnicode; +use hdf5::types::TypeDescriptor; +use crate::scf_io::{SCF, SCFType}; +use crate::geom_io::get_mass_charge; +use crate::basis_io::Basis4Elem; + +pub fn write_scf_attribute(group: &hdf5::Group, dataset_name: &str, value: &[T]) +where + T: Clone + hdf5::H5Type +{ + if let Ok(dataset) = group.dataset(dataset_name) { + match dataset.write_raw(value) { + Ok(_) => (), + Err(e) => println!("Error writing dataset {}: {:?}", dataset_name, e), + } + } else { + let builder = group.new_dataset_builder(); + match builder.with_data(value).create(dataset_name) { + Ok(_) => (), + Err(e) => println!("Error creating dataset {}: {:?}", dataset_name, e), + } + } +} + +pub fn write_string_scalar(file: &hdf5::File, dataset_name: &str, value: &str) { + use std::str::FromStr; + let unicode_string = VarLenUnicode::from_str(value).unwrap(); + + if let Ok(dataset) = file.dataset(dataset_name) { + // Dataset exists - overwrite it + dataset.write_scalar(&unicode_string).unwrap(); + } else { + // Dataset doesn't exist - create it with string type + let ds = file.new_dataset::().shape(()) + .create(dataset_name).unwrap(); + ds.write_scalar(&unicode_string).unwrap(); + } +} + +pub fn save_chkfile(scf_data: &SCF) { + let chkfile= &scf_data.mol.ctrl.chkfile; + let path = Path::new(chkfile); + //if path.exists() {std::fs::remove_file(chkfile).unwrap()}; + //let file = hdf5::File::create(chkfile).unwrap(); + //let scf = file.create_group("scf").unwrap(); + let file = if path.exists() { + hdf5::File::open_rw(chkfile).unwrap() + } else { + hdf5::File::create(chkfile).unwrap() + }; + println!("write chkfile: {}", chkfile); + let is_exist = file.member_names().unwrap().iter().fold(false,|is_exist,x| {is_exist || x.eq("scf")}); + let scf = if is_exist { + file.group("scf").unwrap() + } else { + file.create_group("scf").unwrap() + }; + + write_scf_attribute(&scf, "e_tot", &[scf_data.scf_energy]); + write_scf_attribute(&scf, "num_basis", &[scf_data.mol.num_basis]); + write_scf_attribute(&scf, "spin_channel", &[scf_data.mol.spin_channel]); + write_scf_attribute(&scf, "num_states", &[scf_data.mol.num_state]); + + // let is_exist = scf.member_names().unwrap().iter().fold(false,|is_exist,x| {is_exist || x.eq("mo_coeff")}); + let mut eigenvectors: Vec = vec![]; + for i_spin in 0..scf_data.mol.spin_channel { + let tmp_eigenvectors = scf_data.eigenvectors[i_spin].transpose(); + eigenvectors.extend(tmp_eigenvectors.data.iter()); + if let SCFType::ROHF = scf_data.scftype { // ROHF: only process i_spin = 0, since alpha/beta eigenvectors are the same + break + } + } + write_scf_attribute(&scf, "mo_coeff", &eigenvectors); + + let mut eigenvalues: Vec = vec![]; + for i_spin in 0..scf_data.mol.spin_channel { + eigenvalues.extend(scf_data.eigenvalues[i_spin].iter()); + if let SCFType::ROHF = scf_data.scftype { // ROHF: only process i_spin = 0, since alpha/beta eigenvalues are the same + break + } + } + write_scf_attribute(&scf, "mo_energy", &eigenvalues); + + let mut occ: Vec = vec![]; + for i_spin in 0..scf_data.mol.spin_channel { + occ.extend(scf_data.occupation[i_spin].iter()); + } + // for compatibility with old rest, may be removed in the future + write_scf_attribute(&scf, "mo_occupation", &occ); + // for compatibility with pyscf + write_scf_attribute(&scf, "mo_occ", &occ); + + let mol = &scf_data.mol; + // let (atm, bas, env) = (mol.cint_atm.clone(), mol.cint_bas.clone(), mol.cint_env.clone()); + // convert hashmap of atm, bas, env to one json string + let mol_info = serde_json::to_string(&serde_json::json!({ + "_atm": mol.cint_atm.clone(), + "_bas": mol.cint_bas.clone(), + "_ecpbas": mol.cint_ecpbas.clone(), + "_env": mol.cint_env.clone(), + })).unwrap(); + write_string_scalar(&file, "mol", &mol_info); + let basis4elem = serde_json::to_string(&mol.basis4elem).unwrap(); + write_string_scalar(&file, "molecule/basis4elem", &basis4elem); + + file.close(); +} + +pub fn save_hamiltonian(scf_data: &SCF) { + let chkfile= &scf_data.mol.ctrl.chkfile; + let path = Path::new(chkfile); + let file = if path.exists() { + hdf5::File::open_rw(chkfile).unwrap() + } else { + hdf5::File::create(chkfile).unwrap() + }; + let is_exist = file.member_names().unwrap().iter().fold(false,|is_exist,x| {is_exist || x.eq("scf")}); + let scf = if is_exist { + file.group("scf").unwrap() + } else { + file.create_group("scf").unwrap() + }; + let mut hamiltonians: Vec = vec![]; + for i_spin in 0..scf_data.mol.spin_channel { + let tmp_eigenvectors = scf_data.hamiltonian[i_spin].to_matrixfull().unwrap(); + hamiltonians.extend(scf_data.eigenvalues[i_spin].iter()); + } + let is_hamiltonian = scf.member_names().unwrap().iter().fold(false,|is_exist,x| {is_exist || x.eq("hamiltonian")}); + if is_hamiltonian { + let dataset = scf.dataset("hamiltonian").unwrap(); + dataset.write(&hamiltonians); + } else { + let builder = scf.new_dataset_builder(); + builder.with_data(&hamiltonians).create("hamiltonian"); + }; + file.close(); +} + +pub fn save_overlap(scf_data: &SCF) { + let chkfile= &scf_data.mol.ctrl.chkfile; + let path = Path::new(chkfile); + let file = if path.exists() { + hdf5::File::open_rw(chkfile).unwrap() + } else { + hdf5::File::create(chkfile).unwrap() + }; + let is_exist = file.member_names().unwrap().iter().fold(false,|is_exist,x| {is_exist || x.eq("scf")}); + let scf = if is_exist { + file.group("scf").unwrap() + } else { + file.create_group("scf").unwrap() + }; + let mut overlap: Vec = vec![]; + let overlap = scf_data.ovlp.to_matrixfull().unwrap(); + let is_exist = scf.member_names().unwrap().iter().fold(false,|is_exist,x| {is_exist || x.eq("overlap")}); + if is_exist { + let dataset = scf.dataset("overlap").unwrap(); + dataset.write(&overlap.data); + } else { + let builder = scf.new_dataset_builder(); + builder.with_data(&overlap.data).create("overlap"); + }; + file.close(); +} + +pub fn save_geometry(scf_data: &SCF) { + let ang = crate::constants::ANG; + let chkfile= &scf_data.mol.ctrl.chkfile; + let path = Path::new(chkfile); + let file = if path.exists() { + hdf5::File::open_rw(chkfile).unwrap() + } else { + hdf5::File::create(chkfile).unwrap() + }; + let is_geom = file.member_names().unwrap().iter().fold(false,|is_exist,x| {is_exist || x.eq("geom")}); + let geom = if is_geom { + file.group("geom").unwrap() + } else { + file.create_group("geom").unwrap() + }; + let mass_charge = get_mass_charge(&scf_data.mol.geom.elem); + //let mut geometry: Vec<(f64,f64,f64,f64)> = vec![]; + let mut geometry: Vec<[f64;4]> = vec![]; + mass_charge.iter().zip(scf_data.mol.geom.position.iter_columns_full()).for_each(|(mass_charge, position)| { + geometry.push([mass_charge.1,position[0]*ang,position[1]*ang,position[2]*ang]); + //geometry.push((mass_charge.1,position[0]*ang,position[1]*ang,position[2]*ang)); + }); + + let is_geom = geom.member_names().unwrap().iter().fold(false,|is_exist,x| {is_exist || x.eq("position")}); + if is_geom { + let dataset = geom.dataset("position").unwrap(); + dataset.write(&geometry); + } else { + let builder = geom.new_dataset_builder(); + builder.with_data(&geometry).create("position"); + } + file.close(); +} + +pub fn has_dm(file: &hdf5::File) -> bool { + let mut has_dm = false; + if let Ok(_init_guess) = file.dataset("init_guess") { + has_dm = true; + } + + has_dm +} + +pub fn has_mo_coeff(file: &hdf5::File) -> bool { + let mut has_mo = false; + if let Ok(_mo_coeff) = file.dataset("scf/mo_coeff") { + has_mo = true; + } + + has_mo +} + +pub fn load_cint_data(chkfile: &String) -> (Option<(Vec>, Vec>, Vec)>, Option>>, Option>) { + let file = hdf5::File::open(chkfile).unwrap(); + let mol_info = file.dataset("mol").unwrap().read_scalar::().unwrap(); + let json_string = mol_info.as_str(); + let mol: HashMap = serde_json::from_str(json_string).unwrap(); + let atm = serde_json::from_value(mol.get("_atm").unwrap().clone()).unwrap(); + let bas = serde_json::from_value(mol.get("_bas").unwrap().clone()).unwrap(); + let env = serde_json::from_value(mol.get("_env").unwrap().clone()).unwrap(); + let ecpbas = if let Some(ecpbas_value) = mol.get("_ecpbas") { + Some(serde_json::from_value(ecpbas_value.clone()).unwrap()) + } else { + None + }; + let basis4elem: Option> = if let Some(basis4elem_value) = mol.get("molecule/basis4elem") { + Some(serde_json::from_value(basis4elem_value.clone()).unwrap()) + } else { + None + }; + + (Some((atm, bas, env)), ecpbas, basis4elem) +} + +pub fn load_shapes(chkfile: &String) -> Option<(usize, usize, usize)> { + let file = hdf5::File::open(chkfile).unwrap(); + let scf = file.group("scf").unwrap(); + let mut num_basis = None; + let mut num_states = None; + let mut spin_channel = None; + if let Ok(nb) = scf.dataset("num_basis").unwrap().read_raw::() { + num_basis = Some(nb[0]); + } + if let Ok(nmo) = scf.dataset("num_states").unwrap().read_raw::() { + num_states = Some(nmo[0]); + } + if let Ok(s) = scf.dataset("spin_channel").unwrap().read_raw::() { + spin_channel = Some(s[0]); + } + // println!("Loaded shapes from chkfile: num_basis = {:?}, num_states = {:?}, spin_channel = {:?}", num_basis, num_states, spin_channel); + if num_basis.is_none() || num_states.is_none() || spin_channel.is_none() { + let mo_coeff_data = scf.dataset("mo_coeff").unwrap(); + match mo_coeff_data.ndim() { + 1 => {}, + 2 => { + num_basis = Some(mo_coeff_data.shape()[0]); + num_states = Some(mo_coeff_data.shape()[1]); + spin_channel = Some(1); + }, + 3 => { + spin_channel = Some(mo_coeff_data.shape()[0]); + num_basis = Some(mo_coeff_data.shape()[1]); + num_states = Some(mo_coeff_data.shape()[2]); + }, + _ => { + panic!("Unsupported shape of mo_coeff dataset: {:?}", mo_coeff_data.shape()); + } + } + } + if num_basis.is_some() && num_states.is_some() && spin_channel.is_some() { + Some((num_basis.unwrap(), num_states.unwrap(), spin_channel.unwrap())) + } else { + None + } + +} \ No newline at end of file diff --git a/src/fileop/mod.rs b/src/fileop/mod.rs new file mode 100644 index 0000000..0983024 --- /dev/null +++ b/src/fileop/mod.rs @@ -0,0 +1 @@ +pub mod chkfile; \ No newline at end of file diff --git a/src/initial_guess/chkfile.rs b/src/initial_guess/chkfile.rs deleted file mode 100644 index 329d9a4..0000000 --- a/src/initial_guess/chkfile.rs +++ /dev/null @@ -1,17 +0,0 @@ -pub fn has_dm(file: &hdf5::File) -> bool { - let mut has_dm = false; - if let Ok(_init_guess) = file.dataset("init_guess") { - has_dm = true; - } - - has_dm -} - -pub fn has_mo_coeff(file: &hdf5::File) -> bool { - let mut has_mo = false; - if let Ok(_mo_coeff) = file.dataset("scf/mo_coeff") { - has_mo = true; - } - - has_mo -} \ No newline at end of file diff --git a/src/initial_guess/enxc.rs b/src/initial_guess/enxc.rs index 69100db..05ad52e 100644 --- a/src/initial_guess/enxc.rs +++ b/src/initial_guess/enxc.rs @@ -2,19 +2,16 @@ use core::panic; use std::fs; use rest_libcint::{CintType, CINTR2CDATA}; use serde::{Deserialize, Serialize}; -use serde_json::{Result,Value}; +use serde_json::{Value}; use anyhow; -use tensors::{MatrixFull, MatrixUpper, TensorSlice, TensorSliceMut}; +use tensors::{MatrixFull, MatrixUpper}; use pyo3::pyclass; -use crate::{molecule_io::Molecule, scf_io::SCF}; +use crate::{molecule_io::Molecule}; use crate::constants::{ATM_NUC_MOD_OF, NUC_ECP, SPECIES_INFO}; use rest_libcint::prelude::rest_libcint_wrapper::ECPscalar; use crate::basis_io::ecp::{PotCell, PotCellRaw}; -#[cfg(target_os = "linux")] -use libc::TCA_DUMP_INVISIBLE; - //#[derive(Clone, Debug,Serialize,Deserialize)] //#[pyclass] diff --git a/src/initial_guess/mod.rs b/src/initial_guess/mod.rs index 0590bef..29a6bfc 100644 --- a/src/initial_guess/mod.rs +++ b/src/initial_guess/mod.rs @@ -1,23 +1,20 @@ -use tensors::{MatrixFull, MatrixUpper, BasicMatrix}; - -use crate::constants::E; -use crate::dft::{numerical_density_rayon, numerical_orbital_population}; +#![warn(unused_imports)] +use tensors::{MatrixFull, MatrixUpper}; use crate::initial_guess::enxc::effective_nxc_matrix; use crate::mpi_io::{mpi_broadcast_matrixfull, MPIOperator}; -use crate::scf_io::{scf, SCFType}; -use crate::{molecule_io::Molecule, scf_io::SCF, dft::Grids}; +use crate::scf_io::{SCF, SCFType}; +use crate::{molecule_io::Molecule, dft::Grids}; use crate::initial_guess::sap::get_vsap; use self::sad::initial_guess_from_sad; -//use crate::initial_guess::*; pub mod sap; pub mod sad; pub mod enxc; +pub mod proj; mod pyrest_enxc; -mod chkfile; -use chkfile::{has_dm, has_mo_coeff}; +use crate::fileop::chkfile; enum RESTART { HDF5, @@ -35,7 +32,7 @@ pub fn initial_guess(scf_data: &mut SCF, mpi_operator: &Option) { assert!(std::path::Path::new(&scf_data.mol.ctrl.guessfile).exists(), "The specified guessfile is missing \n({})", &scf_data.mol.ctrl.guessfile); assert!(scf_data.mol.ctrl.guessfile_type.eq(&"hdf5"), "at present only hdf5 type guess file is supported"); let file = hdf5::File::open(&scf_data.mol.ctrl.guessfile).unwrap(); - if has_dm(&file) { + if chkfile::has_dm(&file) { scf_data.density_matrix = initial_guess_from_hdf5guess(&scf_data.mol); // for DFT methods, it needs the eigenvectors to generate the hamiltonian. In consequence, we use the hf method to prepare the eigenvectors from the guess dm scf_data.generate_hf_hamiltonian_for_guess(); @@ -44,7 +41,7 @@ pub fn initial_guess(scf_data: &mut SCF, mpi_operator: &Option) { scf_data.diagonalize_hamiltonian(mpi_operator); scf_data.generate_occupation(); scf_data.generate_density_matrix(); - } else if has_mo_coeff(&file) { + } else if chkfile::has_mo_coeff(&file) { println!("Read MO coefficients from guessfile: {}", &scf_data.mol.ctrl.guessfile); // let (eigenvectors, eigenvalues, is_occupation) = initial_guess_from_hdf5chk( // &scf_data.mol, &scf_data.scftype, &scf_data.mol.ctrl.guessfile); @@ -180,7 +177,7 @@ pub fn initial_guess(scf_data: &mut SCF, mpi_operator: &Option) { pub fn update_scf_from_hdf5chk(scf_data: &mut SCF, chkfile: String) { let (eigenvectors, eigenvalues, is_occupation) = initial_guess_from_hdf5chk( - &scf_data.mol, &scf_data.scftype, &chkfile); + &mut scf_data.mol, &scf_data.scftype, &chkfile); //============================= // for MOM projection @@ -254,7 +251,7 @@ pub fn initial_guess_from_hdf5guess(mol: &Molecule) -> Vec> { dm } -pub fn initial_guess_from_hdf5chk(mol: &Molecule, scftype: &SCFType, chkfile: &String) -> ([MatrixFull;2],[Vec;2],Option<[Vec;2]>) { +pub fn initial_guess_from_hdf5chk(mol: &mut Molecule, scftype: &SCFType, chkfile: &String) -> ([MatrixFull;2],[Vec;2],Option<[Vec;2]>) { let mut spin_channel: usize; if let &SCFType::ROHF = scftype { spin_channel = 1; @@ -262,24 +259,55 @@ pub fn initial_guess_from_hdf5chk(mol: &Molecule, scftype: &SCFType, chkfile: &S } else { spin_channel = mol.spin_channel; } - initial_guess_from_hdf5chkfile(chkfile, + let (loaded_eigenvectors,loaded_eigenvalues,loaded_occupation) = import_guess_from_hdf5chkfile(chkfile, + spin_channel, + // mol.num_state, + // mol.num_basis, + mol.ctrl.print_level + ); + let (loaded_nbasis, loaded_nmo, loaded_spin_channel) = chkfile::load_shapes(chkfile).unwrap(); + let basis_match = (loaded_nbasis == mol.num_basis) && (loaded_nmo == mol.num_state); + let chkbasis = mol.ctrl.basis_path == "chkfile"; + if chkbasis { + println!("taking basis set information from chkfile"); + let (cint_raw_data, ecp_raw, basis4elem) = chkfile::load_cint_data(&chkfile); + if let Some(cint_raw_data) = cint_raw_data { + mol.cint_atm = cint_raw_data.0; + mol.cint_bas = cint_raw_data.1; + mol.cint_env = cint_raw_data.2; + mol.cint_ecpbas = ecp_raw; + mol.update_from_cint(loaded_nbasis, basis4elem); + return initial_guess_from_raw( + loaded_eigenvectors, + loaded_eigenvalues, + loaded_occupation.unwrap(), + spin_channel, + mol.num_state, + mol.num_basis, + mol.ctrl.print_level + ); + } else { + panic!("Failed to load the basis set information from chkfile"); + } + } else if basis_match { + return initial_guess_from_raw( + loaded_eigenvectors, + loaded_eigenvalues, + loaded_occupation.unwrap(), spin_channel, mol.num_state, mol.num_basis, - mol.ctrl.print_level) -} - -pub fn initial_guess_from_hdf5chkfile( - chkname: &str, - spin_channel: usize, - num_state: usize, - num_basis: usize, - print_level: usize -) -> ([MatrixFull;2],[Vec;2], Option<[Vec;2]>) { + mol.ctrl.print_level + ); + } else { + println!("The basis set in chkfile does not match the input basis set, trying to project..."); + return proj::proj_mo( + ); + } - //let mut tmp_scf = SCF::new(&mol); - //tmp_scf.generate_occupation(); +} +pub fn import_guess_from_hdf5chkfile(chkname: &str, spin_channel: usize, print_level: usize) -> (Vec,Vec, Option>) { let file = hdf5::File::open(chkname).unwrap(); let scf = file.group("scf").unwrap(); let member = scf.member_names().unwrap(); @@ -290,63 +318,85 @@ pub fn initial_guess_from_hdf5chkfile( if print_level>0 {println!("E_tot from chkfile: {:18.10}", e_tot)}; // importing MO coefficients let buf01 = scf.dataset("mo_coeff").unwrap().read_raw::().unwrap(); - if buf01.len() != num_state*num_basis*spin_channel { - panic!("Inconsistency happens when importing the molecular coefficients from \'{}\':\n buf01 length: {}, num_state*num_basis*spin_channel: {}", - chkname, buf01.len(), num_state*num_basis*spin_channel); - } // importing MO eigenvalues let buf02 = scf.dataset("mo_energy").unwrap().read_raw::().unwrap(); - if buf02.len() != num_state*spin_channel { - panic!("Inconsistency happens when importing the molecular orbital eigenvalues from \'{}\':\n buf02 length: {}, num_state*spin_channel: {}", - chkname, buf02.len(), num_state*spin_channel); + // importing MO occupation + // let is_exist = scf.member_names().unwrap().iter().fold(false, |is_exist, x| x.eq("mo_occupation")); + let buf03 = Some(scf.dataset("mo_occupation").unwrap().read_raw::().unwrap()); + (buf01, buf02, buf03) +} + +pub fn initial_guess_from_raw( + loaded_eigenvectors: Vec, + loaded_eigenvalues: Vec, + loaded_occupation: Vec, + spin_channel: usize, + num_state: usize, + num_basis: usize, + print_level: usize +) -> ([MatrixFull;2],[Vec;2], Option<[Vec;2]>) { + + let mut mode = ""; + if loaded_eigenvectors.len() == num_state*num_basis*spin_channel { + mode = "normal"; // r2r or u2u + } else if loaded_eigenvectors.len() == num_state*num_basis && spin_channel==2 { + mode = "r2u"; + } else if loaded_eigenvalues.len() == num_state*num_basis*2 && spin_channel==1 { + mode = "u2r"; + panic!("Importing UHF MOs in RHF calculation is not supported at present"); + } else { + panic!("Inconsistency happens when importing the MOs:\n loaded_eigenvectors length: {}, num_state*num_basis*spin_channel: {}, num_state*num_basis: {}", + loaded_eigenvectors.len(), num_state*num_basis*spin_channel, num_state*num_basis); } let mut tmp_eigenvectors: [MatrixFull;2] = [MatrixFull::empty(),MatrixFull::empty()]; let mut tmp_eigenvalues: [Vec;2] = [vec![],vec![]]; + let mut tmp_occupation: [Vec;2] = [vec![],vec![]]; + match mode { + "normal" => { (0..spin_channel).into_iter().for_each(|i| { let start = (0 + i)*num_state*num_basis; let end = (1 + i)*num_state*num_basis; - let tmp_eigen = MatrixFull::from_vec([num_state,num_basis], buf01[start..end].to_vec()).unwrap(); + let tmp_eigen = MatrixFull::from_vec([num_state,num_basis], loaded_eigenvectors[start..end].to_vec()).unwrap(); tmp_eigenvectors[i]=tmp_eigen.transpose_and_drop(); //tmp_scf.eigenvectors[i] = tmp_eigenvectors[i].transpose(); - tmp_eigenvalues[i]=buf02[ (0+i)*num_state..(1+i)*num_state].to_vec(); + tmp_eigenvalues[i]=loaded_eigenvalues[ (0+i)*num_state..(1+i)*num_state].to_vec(); }); if print_level>3 { (0..spin_channel).into_iter().for_each(|i| { tmp_eigenvectors[i].formated_output(5, "full"); println!("eigenval {:?}", &tmp_eigenvalues[i]); }); - } + } + (0..spin_channel).into_iter().for_each(|i_spin| { + tmp_occupation[i_spin]=loaded_occupation[ (0+i_spin)*num_state..(1+i_spin)*num_state].to_vec(); + }); + println!("tmp_eigenvectors {:?}, tmp_eigenvalues {:?}, tmp_occupation {:?}", &tmp_eigenvectors, &tmp_eigenvalues, &tmp_occupation); + }, + "r2u" => { + println!("Importing MO coefficients in RHF format and converting them to UHF format"); + (0..spin_channel).into_iter().for_each(|i| { + let start = 0; + let end = num_state*num_basis; - // importing MO occupation - let is_exist = scf.member_names().unwrap().iter().fold(false, |is_exist, x| x.eq("mo_occupation")); - let buf03 = if is_exist { - Some(scf.dataset("mo_occupation").unwrap().read_raw::().unwrap()) - } else { - None - }; - let tmp_occupation = if let Some(tmp_buf03) = &buf03 { - if tmp_buf03.len() != num_state*spin_channel { - panic!("Inconsistency happens when importing the molecular orbital occupation from \'{}\':\n buf03 length: {}, num_state*spin_channel: {}", - chkname, tmp_buf03.len(), num_state*spin_channel); - let mut tmp_occupation_0 = [vec![],vec![]]; + let tmp_eigen = MatrixFull::from_vec([num_state,num_basis], loaded_eigenvectors[start..end].to_vec()).unwrap(); + tmp_eigenvectors[i]=tmp_eigen.transpose_and_drop(); + + tmp_eigenvalues[i]=loaded_eigenvalues[ 0..num_state].to_vec(); + }); (0..spin_channel).into_iter().for_each(|i_spin| { - tmp_occupation_0[i_spin]=tmp_buf03[ (0+i_spin)*num_state..(1+i_spin)*num_state].to_vec(); + tmp_occupation[i_spin]=loaded_occupation[ 0..num_state].to_vec(); }); - Some(tmp_occupation_0) - } else { - None + }, + "u2r" => {}, + _ => { + panic!("Unknown mode for importing MO coefficients"); } - } else { - None - }; - - //tmp_scf.generate_density_matrix(); + } - //tmp_scf.density_matrix - (tmp_eigenvectors,tmp_eigenvalues,tmp_occupation) + (tmp_eigenvectors,tmp_eigenvalues,Some(tmp_occupation)) } diff --git a/src/initial_guess/proj.rs b/src/initial_guess/proj.rs new file mode 100644 index 0000000..b90dd40 --- /dev/null +++ b/src/initial_guess/proj.rs @@ -0,0 +1,6 @@ +// use crate::scf_io::SCF; +use tensors::matrix::MatrixFull; + +pub fn proj_mo() -> ([MatrixFull;2],[Vec;2],Option<[Vec;2]>) { + unimplemented!() +} \ No newline at end of file diff --git a/src/initial_guess/sad.rs b/src/initial_guess/sad.rs index eb85420..afe9b9a 100644 --- a/src/initial_guess/sad.rs +++ b/src/initial_guess/sad.rs @@ -1,14 +1,11 @@ use tensors::{MatrixFull, BasicMatrix}; use crate::check_norm::OCCType; -use crate::constants::{F_SHELL, KR_SHELL, NELE_IN_SHELLS, SPECIES_INFO, S_SHELL, XE_SHELL}; use crate::molecule_io::Molecule; use crate::ctrl_io::InputKeywords; use crate::geom_io::{GeomCell, formated_element_name}; use crate::mpi_io::MPIOperator; -use crate::scf_io::{initialize_scf, scf}; -use crate::utilities; +use crate::scf_io::{scf}; use std::collections::HashMap; -use std::num; pub fn initial_guess_from_sad(mol: &Molecule, mpi_operator: &Option) -> Vec> { diff --git a/src/initial_guess/sap.rs b/src/initial_guess/sap.rs index 2dbee40..6d06ce1 100644 --- a/src/initial_guess/sap.rs +++ b/src/initial_guess/sap.rs @@ -1,13 +1,11 @@ use crate::molecule_io; use crate::geom_io; -use crate::dft::Grids as dftgrids; use crate::dft; use rayon::prelude::IndexedParallelIterator; use rayon::prelude::IntoParallelRefIterator; use rayon::prelude::IntoParallelRefMutIterator; use rayon::prelude::ParallelIterator; -//use crate::scf::io; -use rest_tensors::{MatrixFull, matrix_blas_lapack::_einsum_01_rayon,matrix_blas_lapack::_dgemm_nn}; +use rest_tensors::{MatrixFull, matrix_blas_lapack::_einsum_01_rayon}; use tensors::matrix_blas_lapack::_dgemm_full; use tensors::BasicMatrix; use crate::constants::ARR; diff --git a/src/lib.rs b/src/lib.rs index f2cdbb1..c31eccb 100644 --- a/src/lib.rs +++ b/src/lib.rs @@ -77,6 +77,7 @@ pub mod ri_gw; pub mod ri_bse; pub mod solvent; pub mod ri_tddft; +pub mod fileop; //extern crate rest; diff --git a/src/molecule_io/mod.rs b/src/molecule_io/mod.rs index 77c8a32..f123178 100644 --- a/src/molecule_io/mod.rs +++ b/src/molecule_io/mod.rs @@ -340,6 +340,14 @@ impl Molecule { Ok(mol) } + + pub fn update_from_cint(&mut self, num_basis: usize, basis4elem: Option>) { + self.num_basis = num_basis; + if let Some(basis4elem_value) = basis4elem { + self.basis4elem = basis4elem_value; + } + } + pub fn initialize_auxbas(&mut self) { let cint_type = if self.ctrl.basis_type.to_lowercase()==String::from("spheric") { CintType::Spheric diff --git a/src/post_scf_analysis/cube_build.rs b/src/post_scf_analysis/cube_build.rs index 0a5704d..7a3da91 100644 --- a/src/post_scf_analysis/cube_build.rs +++ b/src/post_scf_analysis/cube_build.rs @@ -5,9 +5,8 @@ use crate::basis_io; use crate::scf_io::SCF; use rest_tensors::MatrixFull; use tensors::matrix_blas_lapack::_dgemm_full; -use std::path::Iter; use rayon::iter::{IntoParallelRefIterator, IndexedParallelIterator, ParallelIterator, IntoParallelRefMutIterator}; -use std::fs::{self, File}; +use std::fs::{File}; use std::io::prelude::*; use std::io::LineWriter; diff --git a/src/post_scf_analysis/mod.rs b/src/post_scf_analysis/mod.rs index 390394a..9e62801 100644 --- a/src/post_scf_analysis/mod.rs +++ b/src/post_scf_analysis/mod.rs @@ -1,3 +1,4 @@ +#![warn(unused_imports)] pub mod rand_wf_real_space; pub mod cube_build; pub mod molden_build; @@ -6,26 +7,24 @@ pub mod strong_correlation_correction; pub mod rrs_pbc; pub mod spin_correction; -use std::path::Path; use rest_libcint::prelude::rest_libcint_wrapper::int1e_r; -use tensors::{MathMatrix, MatrixFull, RIFull}; +use tensors::{MathMatrix, RIFull}; -use crate::constants::{ANG, AU2DEBYE, SPECIES_INFO}; -use crate::dft::{DFAFamily}; -use crate::geom_io::get_mass_charge; +use crate::constants::{ANG, AU2DEBYE}; use crate::grad::{formated_force, formated_force_ev, numerical_force}; use crate::mpi_io::MPIOperator; -use crate::ri_pt2::sbge2::{close_shell_sbge2_rayon, open_shell_sbge2_rayon, close_shell_sbge2_detailed_rayon, open_shell_sbge2_detailed_rayon}; -use crate::ri_rpa::scsrpa::{evaluate_osrpa_correlation_rayon, evaluate_spin_response_rayon, evaluate_special_radius_only}; -use crate::ri_rpa::{evaluate_rpa_correlation, evaluate_rpa_correlation_rayon}; +use crate::ri_pt2::sbge2::{close_shell_sbge2_rayon, open_shell_sbge2_rayon}; +use crate::ri_rpa::scsrpa::{evaluate_osrpa_correlation_rayon}; +use crate::ri_rpa::{evaluate_rpa_correlation_rayon}; use crate::ri_gw; use crate::ri_bse; use crate::scf_io::{SCF, SCFType, print_force_for_ghost_point_charges}; use crate::ri_pt2::{close_shell_pt2_rayon, open_shell_pt2_rayon}; use crate::utilities::TimeRecords; -use self::molden_build::{gen_header, gen_molden}; +use self::molden_build::{gen_molden}; use self::strong_correlation_correction::scc15_for_rxdh7; +pub use crate::fileop::chkfile::{save_chkfile, save_hamiltonian, save_overlap, save_geometry}; pub fn post_scf_output(scf_data: &SCF, mpi_operator: &Option) { scf_data.mol.ctrl.outputs.iter().for_each(|output_type| { @@ -153,170 +152,6 @@ pub fn post_scf_output(scf_data: &SCF, mpi_operator: &Option) { }); } -pub fn write_scf_attribute(group: &hdf5::Group, dataset_name: &str, value: &[T]) -where - T: Clone + hdf5::H5Type -{ - if let Ok(dataset) = group.dataset(dataset_name) { - match dataset.write_raw(value) { - Ok(_) => (), - Err(e) => println!("Error writing dataset {}: {:?}", dataset_name, e), - } - } else { - let builder = group.new_dataset_builder(); - match builder.with_data(value).create(dataset_name) { - Ok(_) => (), - Err(e) => println!("Error creating dataset {}: {:?}", dataset_name, e), - } - } -} - -pub fn save_chkfile(scf_data: &SCF) { - let chkfile= &scf_data.mol.ctrl.chkfile; - let path = Path::new(chkfile); - //if path.exists() {std::fs::remove_file(chkfile).unwrap()}; - //let file = hdf5::File::create(chkfile).unwrap(); - //let scf = file.create_group("scf").unwrap(); - let file = if path.exists() { - hdf5::File::open_rw(chkfile).unwrap() - } else { - hdf5::File::create(chkfile).unwrap() - }; - println!("write chkfile: {}", chkfile); - let is_exist = file.member_names().unwrap().iter().fold(false,|is_exist,x| {is_exist || x.eq("scf")}); - let scf = if is_exist { - file.group("scf").unwrap() - } else { - file.create_group("scf").unwrap() - }; - - write_scf_attribute(&scf, "e_tot", &[scf_data.scf_energy]); - write_scf_attribute(&scf, "num_basis", &[scf_data.mol.num_basis]); - write_scf_attribute(&scf, "spin_channel", &[scf_data.mol.spin_channel]); - write_scf_attribute(&scf, "num_states", &[scf_data.mol.num_state]); - - // let is_exist = scf.member_names().unwrap().iter().fold(false,|is_exist,x| {is_exist || x.eq("mo_coeff")}); - let mut eigenvectors: Vec = vec![]; - for i_spin in 0..scf_data.mol.spin_channel { - let tmp_eigenvectors = scf_data.eigenvectors[i_spin].transpose(); - eigenvectors.extend(tmp_eigenvectors.data.iter()); - if let SCFType::ROHF = scf_data.scftype { // ROHF: only process i_spin = 0, since alpha/beta eigenvectors are the same - break - } - } - write_scf_attribute(&scf, "mo_coeff", &eigenvectors); - - let mut eigenvalues: Vec = vec![]; - for i_spin in 0..scf_data.mol.spin_channel { - eigenvalues.extend(scf_data.eigenvalues[i_spin].iter()); - if let SCFType::ROHF = scf_data.scftype { // ROHF: only process i_spin = 0, since alpha/beta eigenvalues are the same - break - } - } - write_scf_attribute(&scf, "mo_energy", &eigenvalues); - - let mut occ: Vec = vec![]; - for i_spin in 0..scf_data.mol.spin_channel { - occ.extend(scf_data.occupation[i_spin].iter()); - } - // for compatibility with old rest, may be removed in the future - write_scf_attribute(&scf, "mo_occupation", &occ); - // for compatibility with pyscf - write_scf_attribute(&scf, "mo_occ", &occ); - - file.close(); -} - -pub fn save_hamiltonian(scf_data: &SCF) { - let chkfile= &scf_data.mol.ctrl.chkfile; - let path = Path::new(chkfile); - let file = if path.exists() { - hdf5::File::open_rw(chkfile).unwrap() - } else { - hdf5::File::create(chkfile).unwrap() - }; - let is_exist = file.member_names().unwrap().iter().fold(false,|is_exist,x| {is_exist || x.eq("scf")}); - let scf = if is_exist { - file.group("scf").unwrap() - } else { - file.create_group("scf").unwrap() - }; - let mut hamiltonians: Vec = vec![]; - for i_spin in 0..scf_data.mol.spin_channel { - let tmp_eigenvectors = scf_data.hamiltonian[i_spin].to_matrixfull().unwrap(); - hamiltonians.extend(scf_data.eigenvalues[i_spin].iter()); - } - let is_hamiltonian = scf.member_names().unwrap().iter().fold(false,|is_exist,x| {is_exist || x.eq("hamiltonian")}); - if is_hamiltonian { - let dataset = scf.dataset("hamiltonian").unwrap(); - dataset.write(&hamiltonians); - } else { - let builder = scf.new_dataset_builder(); - builder.with_data(&hamiltonians).create("hamiltonian"); - }; - file.close(); -} - -pub fn save_overlap(scf_data: &SCF) { - let chkfile= &scf_data.mol.ctrl.chkfile; - let path = Path::new(chkfile); - let file = if path.exists() { - hdf5::File::open_rw(chkfile).unwrap() - } else { - hdf5::File::create(chkfile).unwrap() - }; - let is_exist = file.member_names().unwrap().iter().fold(false,|is_exist,x| {is_exist || x.eq("scf")}); - let scf = if is_exist { - file.group("scf").unwrap() - } else { - file.create_group("scf").unwrap() - }; - let mut overlap: Vec = vec![]; - let overlap = scf_data.ovlp.to_matrixfull().unwrap(); - let is_exist = scf.member_names().unwrap().iter().fold(false,|is_exist,x| {is_exist || x.eq("overlap")}); - if is_exist { - let dataset = scf.dataset("overlap").unwrap(); - dataset.write(&overlap.data); - } else { - let builder = scf.new_dataset_builder(); - builder.with_data(&overlap.data).create("overlap"); - }; - file.close(); -} - -pub fn save_geometry(scf_data: &SCF) { - let ang = crate::constants::ANG; - let chkfile= &scf_data.mol.ctrl.chkfile; - let path = Path::new(chkfile); - let file = if path.exists() { - hdf5::File::open_rw(chkfile).unwrap() - } else { - hdf5::File::create(chkfile).unwrap() - }; - let is_geom = file.member_names().unwrap().iter().fold(false,|is_exist,x| {is_exist || x.eq("geom")}); - let geom = if is_geom { - file.group("geom").unwrap() - } else { - file.create_group("geom").unwrap() - }; - let mass_charge = get_mass_charge(&scf_data.mol.geom.elem); - //let mut geometry: Vec<(f64,f64,f64,f64)> = vec![]; - let mut geometry: Vec<[f64;4]> = vec![]; - mass_charge.iter().zip(scf_data.mol.geom.position.iter_columns_full()).for_each(|(mass_charge, position)| { - geometry.push([mass_charge.1,position[0]*ang,position[1]*ang,position[2]*ang]); - //geometry.push((mass_charge.1,position[0]*ang,position[1]*ang,position[2]*ang)); - }); - - let is_geom = geom.member_names().unwrap().iter().fold(false,|is_exist,x| {is_exist || x.eq("position")}); - if is_geom { - let dataset = geom.dataset("position").unwrap(); - dataset.write(&geometry); - } else { - let builder = geom.new_dataset_builder(); - builder.with_data(&geometry).create("position"); - } - file.close(); -} pub fn print_out_dfa(scf_data: &SCF) { let dfa = crate::dft::DFA4REST::new_xc(scf_data.mol.spin_channel, scf_data.mol.ctrl.print_level); diff --git a/src/post_scf_analysis/molden_build.rs b/src/post_scf_analysis/molden_build.rs index bdb96f4..18dd99f 100644 --- a/src/post_scf_analysis/molden_build.rs +++ b/src/post_scf_analysis/molden_build.rs @@ -1,15 +1,9 @@ use crate::molecule_io::Molecule; use crate::geom_io; -use crate::basis_io; use crate::scf_io::SCF; -use rest_tensors::MatrixFull; use rest_libcint::CINTR2CDATA; -use rest_libcint::CintType; -use std::borrow::BorrowMut; use std::collections::HashMap; -use rayon::iter::{IntoParallelRefIterator, IndexedParallelIterator, ParallelIterator, IntoParallelRefMutIterator}; -//mod lib; -use std::fs::{self, File}; +use std::fs::{File}; use std::io::prelude::*; use std::io::LineWriter; diff --git a/src/post_scf_analysis/mulliken.rs b/src/post_scf_analysis/mulliken.rs index 4097811..10eda5f 100644 --- a/src/post_scf_analysis/mulliken.rs +++ b/src/post_scf_analysis/mulliken.rs @@ -1,8 +1,6 @@ -use std::borrow::BorrowMut; - -use crate::{molecule_io,geom_io,scf_io}; +use crate::{geom_io}; use crate::scf_io::SCF; -use rest_tensors::{BasicMatrix, MathMatrix, MatrixFull}; +use rest_tensors::{MatrixFull}; use tensors::matrix_blas_lapack::_dgemm_full; pub fn mulliken_pop(scf_data: &SCF) -> Vec{ diff --git a/src/post_scf_analysis/rand_wf_real_space.rs b/src/post_scf_analysis/rand_wf_real_space.rs index 8e12e5b..157610d 100644 --- a/src/post_scf_analysis/rand_wf_real_space.rs +++ b/src/post_scf_analysis/rand_wf_real_space.rs @@ -1,10 +1,9 @@ -use std::num; use crate::molecule_io::Molecule; use crate::geom_io; use crate::basis_io; use crate::basis_io::basic_math::factorial; use crate::scf_io::SCF; -use rest_tensors::{TensorOpt,MatrixFull}; +use rest_tensors::{MatrixFull}; use rand::distributions::normal::StandardNormal; use itertools::Itertools; use rayon::iter::{IntoParallelRefIterator, IndexedParallelIterator, ParallelIterator, IntoParallelRefMutIterator}; diff --git a/src/post_scf_analysis/rrs_pbc.rs b/src/post_scf_analysis/rrs_pbc.rs index fa6214e..83c8655 100644 --- a/src/post_scf_analysis/rrs_pbc.rs +++ b/src/post_scf_analysis/rrs_pbc.rs @@ -2,8 +2,6 @@ use std::fs::File; use std::io::Write; use crate::scf_io::SCF; use num_complex::Complex; -use rest_tensors::{MatrixFull,MatrixUpper}; -use crate::tensors::BasicMatrix; use std::ffi::c_char; use std::f64::consts::PI; use rayon::prelude::*; diff --git a/src/post_scf_analysis/spin_correction.rs b/src/post_scf_analysis/spin_correction.rs index ef4b8db..2ad2b0d 100644 --- a/src/post_scf_analysis/spin_correction.rs +++ b/src/post_scf_analysis/spin_correction.rs @@ -1,4 +1,4 @@ -use crate::scf_io::{scf, SCFType, SCF}; +use crate::scf_io::{SCFType, SCF}; use crate::scf_io; use crate::utilities::TimeRecords; use crate::mpi_io::MPIOperator; diff --git a/src/post_scf_analysis/strong_correlation_correction.rs b/src/post_scf_analysis/strong_correlation_correction.rs index 1b8996f..c482592 100644 --- a/src/post_scf_analysis/strong_correlation_correction.rs +++ b/src/post_scf_analysis/strong_correlation_correction.rs @@ -1,12 +1,12 @@ use tensors::MatrixFull; -use crate::constants::{E, EV}; +use crate::constants::{EV}; use crate::mpi_io::MPIOperator; use crate::scf_io::{SCF, SCFType}; use crate::ri_rpa::scsrpa::{evaluate_special_radius_only, evaluate_osrpa_correlation_rayon}; use crate::ri_pt2::sbge2::{close_shell_sbge2_detailed_rayon, open_shell_sbge2_detailed_rayon}; -use libm::{erf,erfc,pow}; +use libm::{erf,erfc}; pub fn scc15_for_rxdh7(scf_data: &mut SCF, mpi_operator: &Option) -> f64 { if let (Some(mpi_op), Some(mpi_ix)) = (mpi_operator, &scf_data.mol.mpi_data) { -- Gitee From ac2029daa4014335e81b20e3d666b469b9f5bacb Mon Sep 17 00:00:00 2001 From: Shirong_Wang Date: Sun, 10 May 2026 21:34:54 +0800 Subject: [PATCH 02/14] chkbasis --- src/fileop/chkfile.rs | 2 + src/initial_guess/mod.rs | 29 +++++----- src/molecule_io/mod.rs | 121 ++++++++++++++++++++------------------- src/scf_io/mod.rs | 23 ++++---- src/scf_io/util.rs | 10 ++++ 5 files changed, 99 insertions(+), 86 deletions(-) create mode 100644 src/scf_io/util.rs diff --git a/src/fileop/chkfile.rs b/src/fileop/chkfile.rs index 48c0e50..c1b2eb4 100644 --- a/src/fileop/chkfile.rs +++ b/src/fileop/chkfile.rs @@ -73,6 +73,7 @@ pub fn save_chkfile(scf_data: &SCF) { break } } + // println!("eigenvectors {:?}", scf_data.eigenvectors); write_scf_attribute(&scf, "mo_coeff", &eigenvectors); let mut eigenvalues: Vec = vec![]; @@ -88,6 +89,7 @@ pub fn save_chkfile(scf_data: &SCF) { for i_spin in 0..scf_data.mol.spin_channel { occ.extend(scf_data.occupation[i_spin].iter()); } + // println!("occupation {:?}", scf_data.occupation); // for compatibility with old rest, may be removed in the future write_scf_attribute(&scf, "mo_occupation", &occ); // for compatibility with pyscf diff --git a/src/initial_guess/mod.rs b/src/initial_guess/mod.rs index 29a6bfc..c612889 100644 --- a/src/initial_guess/mod.rs +++ b/src/initial_guess/mod.rs @@ -2,11 +2,9 @@ use tensors::{MatrixFull, MatrixUpper}; use crate::initial_guess::enxc::effective_nxc_matrix; use crate::mpi_io::{mpi_broadcast_matrixfull, MPIOperator}; -use crate::scf_io::{SCF, SCFType}; +use crate::scf_io::{self, SCF, SCFType}; use crate::{molecule_io::Molecule, dft::Grids}; - use crate::initial_guess::sap::get_vsap; - use self::sad::initial_guess_from_sad; pub mod sap; @@ -230,11 +228,16 @@ pub fn update_scf_from_hdf5chk(scf_data: &mut SCF, chkfile: String) { //============================= scf_data.eigenvalues = eigenvalues; scf_data.eigenvectors = eigenvectors; - if let Some(occupation) = is_occupation { - scf_data.occupation = occupation; - } else { + // if let Some(occupation) = is_occupation { + // scf_data.occupation = occupation; + // let (homo, lumo) = scf_io::util::get_homo_lumo_from_occ_integer(&scf_data.occupation, &scf_data.scftype); + // scf_data.homo = homo; + // scf_data.lumo = lumo; + // } else { + // Since homo, lumo is needed, re-generate occupation temporarily scf_data.generate_occupation(); - } + // } + println!("homo {:?}", scf_data.homo); scf_data.generate_density_matrix(); } @@ -272,11 +275,11 @@ pub fn initial_guess_from_hdf5chk(mol: &mut Molecule, scftype: &SCFType, chkfile println!("taking basis set information from chkfile"); let (cint_raw_data, ecp_raw, basis4elem) = chkfile::load_cint_data(&chkfile); if let Some(cint_raw_data) = cint_raw_data { - mol.cint_atm = cint_raw_data.0; - mol.cint_bas = cint_raw_data.1; - mol.cint_env = cint_raw_data.2; - mol.cint_ecpbas = ecp_raw; - mol.update_from_cint(loaded_nbasis, basis4elem); + // mol.cint_atm = cint_raw_data.0; + // mol.cint_bas = cint_raw_data.1; + // mol.cint_env = cint_raw_data.2; + // mol.cint_ecpbas = ecp_raw; + mol.update_from_cint(loaded_nbasis, basis4elem, cint_raw_data, ecp_raw); return initial_guess_from_raw( loaded_eigenvectors, loaded_eigenvalues, @@ -373,7 +376,7 @@ pub fn initial_guess_from_raw( (0..spin_channel).into_iter().for_each(|i_spin| { tmp_occupation[i_spin]=loaded_occupation[ (0+i_spin)*num_state..(1+i_spin)*num_state].to_vec(); }); - println!("tmp_eigenvectors {:?}, tmp_eigenvalues {:?}, tmp_occupation {:?}", &tmp_eigenvectors, &tmp_eigenvalues, &tmp_occupation); + // println!("tmp_eigenvectors {:?}, tmp_eigenvalues {:?}, tmp_occupation {:?}", &tmp_eigenvectors, &tmp_eigenvalues, &tmp_occupation); }, "r2u" => { println!("Importing MO coefficients in RHF format and converting them to UHF format"); diff --git a/src/molecule_io/mod.rs b/src/molecule_io/mod.rs index c486e82..c737df5 100644 --- a/src/molecule_io/mod.rs +++ b/src/molecule_io/mod.rs @@ -186,7 +186,12 @@ impl Molecule { pub fn build_native(mut ctrl: InputKeywords, mut geom: GeomCell, mpi_data: Option) -> anyhow::Result { - let spin_channel = ctrl.spin_channel; + let mut mol = Molecule::init_mol(); + mol.ctrl = ctrl.clone(); + mol.geom = geom.clone(); + mol.mpi_data = mpi_data; + + mol.spin_channel = ctrl.spin_channel; let cint_type = if ctrl.basis_type.to_lowercase()==String::from("spheric") { CintType::Spheric } else if ctrl.basis_type.to_lowercase()==String::from("cartesian") { @@ -195,25 +200,42 @@ impl Molecule { panic!("Error:: Unknown basis type '{}'. Please use either 'spheric' or 'Cartesian'", ctrl.basis_type); }; - - let (mut basis4elem,mut cint_atm,mut cint_bas,cint_env, - fdqc_bas,cint_fdqc,num_elec,num_basis,num_state, cint_ecpbas) - = Molecule::collect_basis(&mut ctrl, &mut geom); - - - - let bas = &basis4elem; - let ecp_electrons = bas.iter().fold(0, |acc, i| { - //println!("debug {:?}", i); - let ecp_electrons = if let Some(num_ecp) = i.ecp_electrons {num_ecp} else {0}; - acc + ecp_electrons - }); - + mol.cint_type = cint_type; + + let chkbasis = ctrl.basis_path == "chkfile"; + if !chkbasis { + let (mut basis4elem,mut cint_atm,mut cint_bas,cint_env, + fdqc_bas,cint_fdqc,num_elec,num_basis,num_state, cint_ecpbas) + = Molecule::collect_basis(&mut ctrl, &mut geom); + + let bas = &basis4elem; + let ecp_electrons = bas.iter().fold(0, |acc, i| { + let ecp_electrons = if let Some(num_ecp) = i.ecp_electrons {num_ecp} else {0}; + acc + ecp_electrons + }); + mol.basis4elem = basis4elem; + mol.cint_atm = cint_atm; + mol.cint_bas = cint_bas; + mol.cint_env = cint_env; + mol.fdqc_bas = fdqc_bas; + mol.cint_fdqc = cint_fdqc; + mol.num_elec = num_elec; + mol.num_basis = num_basis; + mol.num_state = num_state; + mol.cint_ecpbas = cint_ecpbas; + mol.ecp_electrons = ecp_electrons; + } let (mut auxbas , mut cint_aux_atm,mut cint_aux_bas,cint_aux_env, mut fdqc_aux_bas,mut cint_aux_fdqc,num_auxbas) =(Vec::new(),Vec::new(),Vec::new(),Vec::new(),Vec::new(),Vec::new(),0); - + mol.auxbas4elem = auxbas; + mol.cint_aux_atm = cint_aux_atm; + mol.cint_aux_bas = cint_aux_bas; + mol.cint_aux_env = cint_aux_env; + mol.fdqc_aux_bas = fdqc_aux_bas; + mol.cint_aux_fdqc = cint_aux_fdqc; + mol.num_auxbas = num_auxbas; //let mut cint_data = CINTR2CDATA::new(); //cint_data.set_cint_type(&cint_type); @@ -229,12 +251,12 @@ impl Molecule { println!("Using xc_parser {}", ctrl.xc_parser); match ctrl.xc_parser.as_str() { "legacy" => { - let mut cur_xc_data = DFA4REST::new(&ctrl.xc, spin_channel, ctrl.print_level); + let mut cur_xc_data = DFA4REST::new(&ctrl.xc, ctrl.spin_channel, ctrl.print_level); cur_xc_data.update_pt2_params(ctrl.ri_pt2.os_factor, ctrl.ri_pt2.ss_factor); (None, cur_xc_data) }, "parse_xc" => { - let mut dfadef = parse_xc::parse_and_derive(&ctrl.xc, spin_channel, ctrl.print_level); + let mut dfadef = parse_xc::parse_and_derive(&ctrl.xc, ctrl.spin_channel, ctrl.print_level); let (san, err_strings) = dfadef.check_sanity(); if !san { panic!("{}", err_strings.join("\n")); @@ -249,9 +271,9 @@ impl Molecule { }, DFTType::NonStandard => { println!("Warning: DFTType::NonStandard is about to be deprecated. Please use xc_parser = 'parse_xc' instead."); - (None, DFA4REST::new_nonstandard(spin_channel, ctrl.print_level, &ctrl.xc_namelist, &ctrl.xc_paralist, &ctrl.dfa_hybrid_scf)) + (None, DFA4REST::new_nonstandard(ctrl.spin_channel, ctrl.print_level, &ctrl.xc_namelist, &ctrl.xc_paralist, &ctrl.dfa_hybrid_scf)) }, - DFTType::DeepLearning => {(None, DFA4REST::new_deep_learning(spin_channel, ctrl.print_level, &ctrl.xc_model))} + DFTType::DeepLearning => {(None, DFA4REST::new_deep_learning(ctrl.spin_channel, ctrl.print_level, &ctrl.xc_model))} }; xc_data.summary(); @@ -260,12 +282,15 @@ impl Molecule { std::process::exit(0); } } + mol.dfadef = dfadef; + mol.xc_data = xc_data; - let use_eri = xc_data.use_eri(); + mol.use_eri = mol.xc_data.use_eri(); let mut start_mo = count_frozen_core_states(ctrl.frozen_core_postscf, &geom.elem); - let ecp_orbs = ecp_electrons/2; + + let ecp_orbs = mol.ecp_electrons/2; if ecp_orbs == 0 { if ctrl.print_level > 0 { println!("For SCF calculation, no core orbital is frozen by effective core potential (ECP) approximation"); @@ -283,47 +308,17 @@ impl Molecule { } start_mo = start_mo - ecp_orbs } + mol.start_mo = start_mo; if ctrl.print_level>0 { - xc_data.xc_version(); - println!("nbas: {}, natm: {} for standard basis sets", cint_bas.len(),cint_atm.len()); + mol.xc_data.xc_version(); + println!("nbas: {}, natm: {} for standard basis sets", mol.cint_bas.len(), mol.cint_atm.len()); println!("First valence state for the frozen-core algorithm: {:5}", start_mo); }; - let use_solvent = ctrl.solvent_enabled; - let solvent_model=ctrl.solvent_model; - let solv_epsilon = ctrl.solv_epsilon; - let mut mol = Molecule { - ctrl, - mpi_data, - geom, - dfadef, - xc_data, - use_eri, - num_elec, - num_state, - num_basis, - num_auxbas, - start_mo, - ecp_electrons, - spin_channel, - basis4elem, - auxbas4elem: auxbas, - fdqc_bas, - cint_fdqc, - cint_atm, - cint_bas, - cint_ecpbas, - cint_env, - fdqc_aux_bas, - cint_aux_fdqc, - cint_aux_atm, - cint_aux_bas, - cint_aux_env, - cint_type, - use_solvent, - solvent_model, - solv_epsilon, - }; + mol.use_solvent = ctrl.solvent_enabled; + mol.solvent_model = ctrl.solvent_model; + mol.solv_epsilon = ctrl.solv_epsilon; + // check and prepare the auxiliary basis sets if mol.ctrl.use_auxbas {mol.initialize_auxbas()}; @@ -341,11 +336,17 @@ impl Molecule { } - pub fn update_from_cint(&mut self, num_basis: usize, basis4elem: Option>) { + pub fn update_from_cint(&mut self, num_basis: usize, basis4elem: Option>, cint_raw: (Vec>, Vec>, Vec), ecp_raw: Option>>) { self.num_basis = num_basis; if let Some(basis4elem_value) = basis4elem { self.basis4elem = basis4elem_value; } + self.cint_atm = cint_raw.0; + self.cint_bas = cint_raw.1; + self.cint_env = cint_raw.2; + self.cint_ecpbas = ecp_raw; + // update fdqc_bas, cint_fdqc + // update num_elec num_basis num_state cint_ecpbas ecp_electrons; } pub fn initialize_auxbas(&mut self) { diff --git a/src/scf_io/mod.rs b/src/scf_io/mod.rs index 29a116f..6cb2a12 100644 --- a/src/scf_io/mod.rs +++ b/src/scf_io/mod.rs @@ -9,15 +9,12 @@ use crate::mpi_io::{mpi_broadcast, mpi_broadcast_matrixfull, mpi_broadcast_vecto use crate::utilities::{self, TimeRecords}; use crate::utilities::memory_batch::*; use crate::ctrl_io::ri_jk_io::*; - -mod addons; -mod fchk; -mod pyrest_scf_io; - use mpi::collective::SystemOperation; use pyo3::{pyclass}; use tensors::matrix_blas_lapack::{_dgemm, _dgemm_full, _dgemv, _dinverse, _dspgvx, _dsymm, _dsyrk, _hamiltonian_fast_solver, _power_rayon_for_symmetric_matrix, _dsyevd}; +use tensors::matrix_blas_lapack::{omp_get_num_threads_wrapper,omp_set_num_threads_wrapper}; use tensors::{map_upper_to_full, BasicMatUp, BasicMatrix, ERIFold4, MathMatrix, MatrixFull, MatrixFullSlice, MatrixUpper, MatrixUpperSlice, RIFull, TensorSliceMut}; +use tensors::{TensorOpt,TensorOptMut,TensorSlice}; use itertools::{Itertools}; use rayon::prelude::*; use std::collections::HashMap; @@ -25,13 +22,16 @@ use crossbeam::{channel::{unbounded},thread::{scope}}; use std::sync::mpsc::{channel}; use crate::isdf::{prepare_for_ri_isdf, prepare_m_isdf}; use crate::molecule_io::{Molecule}; -use crate::tensors::{TensorOpt,TensorOptMut,TensorSlice}; use crate::initial_guess::initial_guess; use crate::external_libs::dftd; use crate::constants::{SQRT_THRESHOLD}; use crate::solvent::{PcmObject, PcmScf, solvent_prepare}; use crate::ri_jk; -use tensors::matrix_blas_lapack::{omp_get_num_threads_wrapper,omp_set_num_threads_wrapper}; + +mod addons; +mod fchk; +mod pyrest_scf_io; +pub mod util; #[pyclass] #[derive(Clone)] @@ -5287,7 +5287,7 @@ pub fn generate_density_matrix_outside(scf_data: &SCF) -> Vec>{ let num_basis = scf_data.mol.num_basis; let num_state = scf_data.mol.num_state; let spin_channel = scf_data.mol.spin_channel; - let homo = &scf_data.homo; + // let homo = &scf_data.homo; // println!("homo: {:?}", &homo); let mut dm = vec![ MatrixFull::empty(), @@ -5303,7 +5303,7 @@ pub fn generate_density_matrix_outside(scf_data: &SCF) -> Vec>{ }; let occ_s = &scf_data.occupation[i_spin]; - let nw = scf_data.homo[i_spin]+1; + let nw = util::get_nocc_from_occ(&scf_data.occupation)[i_spin]; //println!("number of occupied orbitals from dm generation: {}", nw); let mut weight_eigv = MatrixFull::new([num_basis, num_state],0.0_f64); @@ -5954,7 +5954,4 @@ pub fn print_force_for_ghost_point_charges(scf_data: &SCF) { } } -#[test] -fn test_max() { - println!("{}, {}, {}",1,2,std::cmp::max(1, 2)); -} + diff --git a/src/scf_io/util.rs b/src/scf_io/util.rs new file mode 100644 index 0000000..89dd707 --- /dev/null +++ b/src/scf_io/util.rs @@ -0,0 +1,10 @@ +// use crate::scf_io::{SCF, SCFType, scf} + +pub fn get_nocc_from_occ(occ: &[Vec;2]) -> [usize;2] { + let mut nocc = [0,0]; + for i_spin in 0..2 { + nocc[i_spin] = occ[i_spin].iter().filter(|&&x| x>1e-6).count(); + } + nocc +} + -- Gitee From 7219ad0b3f78050a78ae637886fce6fea5898f38 Mon Sep 17 00:00:00 2001 From: Shirong_Wang Date: Mon, 11 May 2026 21:50:40 +0800 Subject: [PATCH 03/14] update_cint from GPT 5.2 codex --- src/initial_guess/mod.rs | 64 +++++++++++++-------- src/molecule_io/mod.rs | 117 ++++++++++++++++++++++++++++++--------- src/scf_io/mod.rs | 4 +- 3 files changed, 136 insertions(+), 49 deletions(-) diff --git a/src/initial_guess/mod.rs b/src/initial_guess/mod.rs index c612889..32fa9ec 100644 --- a/src/initial_guess/mod.rs +++ b/src/initial_guess/mod.rs @@ -254,6 +254,46 @@ pub fn initial_guess_from_hdf5guess(mol: &Molecule) -> Vec> { dm } +pub fn update_basis_from_hdf5chk(scf_data: &mut SCF) { + let chkbasis = scf_data.mol.ctrl.basis_path == "chkfile"; + if !chkbasis { + return; + } + if scf_data.mol.ctrl.initial_guess.eq(&"inherit") { + return; + } else if scf_data.mol.ctrl.external_init_guess.is_some() { + let chkfile = match scf_data.mol.ctrl.external_init_guess.as_ref().unwrap().as_str() { + "guessfile" => { + assert!(std::path::Path::new(&scf_data.mol.ctrl.guessfile).exists(), "The specified guessfile is missing \n({})", &scf_data.mol.ctrl.guessfile); + assert!(scf_data.mol.ctrl.guessfile_type.eq(&"hdf5"), "at present only hdf5 type guess file is supported"); + scf_data.mol.ctrl.guessfile.clone() + }, + "chkfile" => { + assert!(std::path::Path::new(&scf_data.mol.ctrl.chkfile).exists(), "The specified chkfile is missing \n({})", &scf_data.mol.ctrl.chkfile); + assert!(scf_data.mol.ctrl.chkfile_type.eq(&"hdf5"), "at present only hdf5 type check file is supported"); + scf_data.mol.ctrl.chkfile.clone() + } + _ => { + panic!("Error: unknown external initial guess file type ({})", &scf_data.mol.ctrl.external_init_guess.as_ref().unwrap()); + } + }; + + + println!("taking basis set information from chkfile"); + let (loaded_nbasis, loaded_nmo, loaded_spin_channel) = chkfile::load_shapes(&chkfile).unwrap(); + let (cint_raw_data, ecp_raw, basis4elem) = chkfile::load_cint_data(&chkfile); + if let Some(cint_raw_data) = cint_raw_data { + // mol.cint_atm = cint_raw_data.0; + // mol.cint_bas = cint_raw_data.1; + // mol.cint_env = cint_raw_data.2; + // mol.cint_ecpbas = ecp_raw; + scf_data.mol.update_from_cint(loaded_nbasis, basis4elem, cint_raw_data, ecp_raw); + } else { + panic!("Failed to load the basis set information from chkfile"); + } + } +} + pub fn initial_guess_from_hdf5chk(mol: &mut Molecule, scftype: &SCFType, chkfile: &String) -> ([MatrixFull;2],[Vec;2],Option<[Vec;2]>) { let mut spin_channel: usize; if let &SCFType::ROHF = scftype { @@ -270,29 +310,7 @@ pub fn initial_guess_from_hdf5chk(mol: &mut Molecule, scftype: &SCFType, chkfile ); let (loaded_nbasis, loaded_nmo, loaded_spin_channel) = chkfile::load_shapes(chkfile).unwrap(); let basis_match = (loaded_nbasis == mol.num_basis) && (loaded_nmo == mol.num_state); - let chkbasis = mol.ctrl.basis_path == "chkfile"; - if chkbasis { - println!("taking basis set information from chkfile"); - let (cint_raw_data, ecp_raw, basis4elem) = chkfile::load_cint_data(&chkfile); - if let Some(cint_raw_data) = cint_raw_data { - // mol.cint_atm = cint_raw_data.0; - // mol.cint_bas = cint_raw_data.1; - // mol.cint_env = cint_raw_data.2; - // mol.cint_ecpbas = ecp_raw; - mol.update_from_cint(loaded_nbasis, basis4elem, cint_raw_data, ecp_raw); - return initial_guess_from_raw( - loaded_eigenvectors, - loaded_eigenvalues, - loaded_occupation.unwrap(), - spin_channel, - mol.num_state, - mol.num_basis, - mol.ctrl.print_level - ); - } else { - panic!("Failed to load the basis set information from chkfile"); - } - } else if basis_match { + if basis_match { return initial_guess_from_raw( loaded_eigenvectors, loaded_eigenvalues, diff --git a/src/molecule_io/mod.rs b/src/molecule_io/mod.rs index c737df5..592e064 100644 --- a/src/molecule_io/mod.rs +++ b/src/molecule_io/mod.rs @@ -125,6 +125,28 @@ pub struct Molecule { } impl Molecule { + fn generate_start_mo(&self, ecp_electrons: usize) -> usize { + let mut start_mo = count_frozen_core_states(self.ctrl.frozen_core_postscf, &self.geom.elem); + let ecp_orbs = ecp_electrons/2; + if ecp_orbs == 0 { + if self.ctrl.print_level > 0 { + println!("For SCF calculation, no core orbital is frozen by effective core potential (ECP) approximation"); + } + } else if start_mo < ecp_orbs { + if self.ctrl.print_level > 0 { + println!(" for the frozen-core post-SCF methods {} is smaller than the number of ecp orbitals {}. Set = 0", + start_mo, ecp_orbs); + } + start_mo = 0; + } else { + if self.ctrl.print_level > 0 { + println!(" for the frozen-core post-SCF methods {} is larger than the number of ecp orbitals {}. Take the -= ecp_orbitals", + start_mo, ecp_orbs); + } + start_mo = start_mo - ecp_orbs + } + start_mo + } pub fn init_mol() -> Molecule { Molecule { ctrl:InputKeywords::init_ctrl(), @@ -288,32 +310,12 @@ impl Molecule { mol.use_eri = mol.xc_data.use_eri(); - let mut start_mo = count_frozen_core_states(ctrl.frozen_core_postscf, &geom.elem); - - let ecp_orbs = mol.ecp_electrons/2; - if ecp_orbs == 0 { - if ctrl.print_level > 0 { - println!("For SCF calculation, no core orbital is frozen by effective core potential (ECP) approximation"); - } - } else if start_mo < ecp_orbs { - if ctrl.print_level > 0 { - println!(" for the frozen-core post-SCF methods {} is smaller than the number of ecp orbitals {}. Set = 0", - start_mo, ecp_orbs); - } - start_mo = 0; - } else { - if ctrl.print_level > 0 { - println!(" for the frozen-core post-SCF methods {} is larger than the number of ecp orbitals {}. Take the -= ecp_orbitals", - start_mo, ecp_orbs); - } - start_mo = start_mo - ecp_orbs - } - mol.start_mo = start_mo; + mol.start_mo = mol.generate_start_mo(mol.ecp_electrons); if ctrl.print_level>0 { mol.xc_data.xc_version(); println!("nbas: {}, natm: {} for standard basis sets", mol.cint_bas.len(), mol.cint_atm.len()); - println!("First valence state for the frozen-core algorithm: {:5}", start_mo); + println!("First valence state for the frozen-core algorithm: {:5}", mol.start_mo); }; mol.use_solvent = ctrl.solvent_enabled; mol.solvent_model = ctrl.solvent_model; @@ -346,7 +348,74 @@ impl Molecule { self.cint_env = cint_raw.2; self.cint_ecpbas = ecp_raw; // update fdqc_bas, cint_fdqc - // update num_elec num_basis num_state cint_ecpbas ecp_electrons; + let mut fdqc_bas: Vec = vec![]; + let mut cint_fdqc: Vec> = vec![]; + let mut bas_start = 0_usize; + self.cint_bas.iter().enumerate().for_each(|(bas_index, bas_cell)| { + let atm_index = bas_cell[0] as usize; + let ang = bas_cell[1] as usize; + let num_primitive = bas_cell[2] as usize; + let num_contracted = bas_cell[3] as usize; + let tmp_bas_num = match &self.cint_type { + CintType::Cartesian => (ang+1)*(ang+2)/2, + CintType::Spheric => ang*2+1, + CintType::Spinor => {panic!("Spinor is not yet implemented")}, + }; + let mut tmp_len = 0_usize; + (0..num_contracted).into_iter().for_each(|index0| { + (0..tmp_bas_num).into_iter().for_each(|index1| { + let bas_type = if num_primitive == 1 { + String::from("Primitive") + } else { + String::from("Contracted") + }; + tmp_len += 1; + fdqc_bas.push(BasInfo { + bas_name: get_basis_name(ang, &self.cint_type, index1), + bas_type, + elem_index0: atm_index, + cint_index0: bas_index, + cint_index1: index0*tmp_bas_num+index1, + }); + }); + }); + cint_fdqc.push(vec![bas_start,tmp_len]); + bas_start += tmp_len; + }); + self.fdqc_bas = fdqc_bas; + self.cint_fdqc = cint_fdqc; + + let num_basis_calc = self.fdqc_bas.len(); + if self.num_basis != num_basis_calc { + if self.ctrl.print_level > 0 { + println!("Warning: num_basis {} from chkfile does not match derived value {}. Using derived value.", self.num_basis, num_basis_calc); + } + self.num_basis = num_basis_calc; + } + self.num_state = self.num_basis; + + // update num_elec + let mut num_elec = [0.0;3]; + self.cint_atm.iter().for_each(|atm| { + num_elec[0] += atm[ATM_NUC] as f64; + }); + num_elec[0] -= self.ctrl.charge; + if self.ctrl.use_int_nelec { + sanity_check_nelec(num_elec[0], self.ctrl.spin); + } + let unpair_elec = (self.ctrl.spin-1.0_f64); + num_elec[1] = (num_elec[0]-unpair_elec)/2.0 + unpair_elec; + num_elec[2] = (num_elec[0]-unpair_elec)/2.0; + self.num_elec = num_elec; + + // update ecp_electrons + let ecp_electrons = self.basis4elem.iter().fold(0, |acc, i| { + let ecp_electrons = if let Some(num_ecp) = i.ecp_electrons {num_ecp} else {0}; + acc + ecp_electrons + }); + self.ecp_electrons = ecp_electrons; + + self.start_mo = self.generate_start_mo(self.ecp_electrons); } pub fn initialize_auxbas(&mut self) { @@ -3507,5 +3576,3 @@ fn test_matrixupper() { let matrixupper_index = tensors::map_upper_to_full(10).unwrap(); dd.iter_submatrix(1..3, 0..2, &matrixupper_index).for_each(|x| {println!("{}",x)}); } - - diff --git a/src/scf_io/mod.rs b/src/scf_io/mod.rs index 6cb2a12..3955852 100644 --- a/src/scf_io/mod.rs +++ b/src/scf_io/mod.rs @@ -22,7 +22,7 @@ use crossbeam::{channel::{unbounded},thread::{scope}}; use std::sync::mpsc::{channel}; use crate::isdf::{prepare_for_ri_isdf, prepare_m_isdf}; use crate::molecule_io::{Molecule}; -use crate::initial_guess::initial_guess; +use crate::initial_guess::{initial_guess, update_basis_from_hdf5chk}; use crate::external_libs::dftd; use crate::constants::{SQRT_THRESHOLD}; use crate::solvent::{PcmObject, PcmScf, solvent_prepare}; @@ -5340,6 +5340,8 @@ pub fn initialize_scf(scf_data: &mut SCF, mpi_operator: &Option) { let position = &scf_data.mol.geom.position; scf_data.mol.cint_env = scf_data.mol.update_geom_poisition_in_cint_env(position); + update_basis_from_hdf5chk(scf_data); + // update the RI-JK algorithms if not clearly specified scf_data.update_jk_algorithms(); -- Gitee From aa7a30054e72be891a695aa5cbe9b4b55430081f Mon Sep 17 00:00:00 2001 From: Shirong_Wang Date: Sat, 16 May 2026 14:34:34 +0800 Subject: [PATCH 04/14] clean --- src/scf_io/util.rs | 10 ---------- 1 file changed, 10 deletions(-) diff --git a/src/scf_io/util.rs b/src/scf_io/util.rs index 89dd707..e69de29 100644 --- a/src/scf_io/util.rs +++ b/src/scf_io/util.rs @@ -1,10 +0,0 @@ -// use crate::scf_io::{SCF, SCFType, scf} - -pub fn get_nocc_from_occ(occ: &[Vec;2]) -> [usize;2] { - let mut nocc = [0,0]; - for i_spin in 0..2 { - nocc[i_spin] = occ[i_spin].iter().filter(|&&x| x>1e-6).count(); - } - nocc -} - -- Gitee From 6a2db3b4cc7d1b1784fbfe8928faf4496539b463 Mon Sep 17 00:00:00 2001 From: Shirong_Wang Date: Sat, 16 May 2026 14:45:19 +0800 Subject: [PATCH 05/14] fix mistake --- src/scf_io/mod.rs | 1 + 1 file changed, 1 insertion(+) diff --git a/src/scf_io/mod.rs b/src/scf_io/mod.rs index 9cd6bfa..c90dc22 100644 --- a/src/scf_io/mod.rs +++ b/src/scf_io/mod.rs @@ -38,6 +38,7 @@ use self::util::occupied_orbital_count; #[derive(Clone)] pub struct SCF { #[pyo3(get,set)] + pub mol: Molecule, pub ovlp: MatrixUpper, pub h_core: MatrixUpper, //pub ijkl: Option>, -- Gitee From 7bbbe94f9a80b67077f88ef708a4b7026e6e103c Mon Sep 17 00:00:00 2001 From: Shirong_Wang Date: Sat, 16 May 2026 18:55:55 +0800 Subject: [PATCH 06/14] finish chkbasis --- src/fileop/chkfile.rs | 12 +++++++++--- 1 file changed, 9 insertions(+), 3 deletions(-) diff --git a/src/fileop/chkfile.rs b/src/fileop/chkfile.rs index c1b2eb4..e2e15d3 100644 --- a/src/fileop/chkfile.rs +++ b/src/fileop/chkfile.rs @@ -229,12 +229,18 @@ pub fn load_cint_data(chkfile: &String) -> (Option<(Vec>, Vec> let bas = serde_json::from_value(mol.get("_bas").unwrap().clone()).unwrap(); let env = serde_json::from_value(mol.get("_env").unwrap().clone()).unwrap(); let ecpbas = if let Some(ecpbas_value) = mol.get("_ecpbas") { - Some(serde_json::from_value(ecpbas_value.clone()).unwrap()) + // Some(serde_json::from_value(ecpbas_value.clone()).unwrap()) + if let Some(ecpbas_vec) = serde_json::from_value(ecpbas_value.clone()).ok() { + Some(ecpbas_vec) + } else { + None + } } else { None }; - let basis4elem: Option> = if let Some(basis4elem_value) = mol.get("molecule/basis4elem") { - Some(serde_json::from_value(basis4elem_value.clone()).unwrap()) + // let bs = file.dataset("molecule/basis4elem").ok(); + let basis4elem: Option> = if let Some(basis4elem_value) = file.dataset("molecule/basis4elem").unwrap().read_scalar::().ok() { + Some(serde_json::from_str(basis4elem_value.as_str()).unwrap()) } else { None }; -- Gitee From 47f67d5e6f4bed53f2ec28f8c1e0493ff1d1357e Mon Sep 17 00:00:00 2001 From: Shirong_Wang Date: Sat, 16 May 2026 19:48:39 +0800 Subject: [PATCH 07/14] prepare for proj --- src/initial_guess/mod.rs | 11 ++++++++++- src/initial_guess/proj.rs | 3 ++- 2 files changed, 12 insertions(+), 2 deletions(-) diff --git a/src/initial_guess/mod.rs b/src/initial_guess/mod.rs index 32fa9ec..e13eb60 100644 --- a/src/initial_guess/mod.rs +++ b/src/initial_guess/mod.rs @@ -322,8 +322,17 @@ pub fn initial_guess_from_hdf5chk(mol: &mut Molecule, scftype: &SCFType, chkfile ); } else { println!("The basis set in chkfile does not match the input basis set, trying to project..."); - return proj::proj_mo( + let (cint_raw_data, ecp_raw, basis4elem) = chkfile::load_cint_data(chkfile); + let (mo1, _, _) = initial_guess_from_raw( + loaded_eigenvectors, + loaded_eigenvalues, + loaded_occupation.unwrap(), + loaded_spin_channel, + loaded_nmo, + loaded_nbasis, + mol.ctrl.print_level ); + return proj::proj_mo(mol, mo1, cint_raw_data.unwrap(), ecp_raw); } } diff --git a/src/initial_guess/proj.rs b/src/initial_guess/proj.rs index b90dd40..559dd2f 100644 --- a/src/initial_guess/proj.rs +++ b/src/initial_guess/proj.rs @@ -1,6 +1,7 @@ // use crate::scf_io::SCF; use tensors::matrix::MatrixFull; +use crate::molecule_io::Molecule; -pub fn proj_mo() -> ([MatrixFull;2],[Vec;2],Option<[Vec;2]>) { +pub fn proj_mo(mol1: &Molecule, mo1: [MatrixFull;2], cint2: (Vec>, Vec>, Vec), ecp2: Option>>) -> ([MatrixFull;2],[Vec;2],Option<[Vec;2]>) { unimplemented!() } \ No newline at end of file -- Gitee From e8d6e2f80ef36abbbe3dd7d564a21bfb572e1380 Mon Sep 17 00:00:00 2001 From: Shirong_Wang Date: Sat, 16 May 2026 20:14:38 +0800 Subject: [PATCH 08/14] intor_cross --- src/initial_guess/proj.rs | 27 +++++++++++++++++++++++++++ src/molecule_io/int_cross.rs | 35 +++++++++++++++++++++++++++++++++++ src/molecule_io/mod.rs | 1 + 3 files changed, 63 insertions(+) create mode 100644 src/molecule_io/int_cross.rs diff --git a/src/initial_guess/proj.rs b/src/initial_guess/proj.rs index 559dd2f..0f1baee 100644 --- a/src/initial_guess/proj.rs +++ b/src/initial_guess/proj.rs @@ -2,6 +2,33 @@ use tensors::matrix::MatrixFull; use crate::molecule_io::Molecule; +// reference implementation: +// +// def project_mo_nr2nr(mol1, mo1, mol2): +// r''' Project orbital coefficients from basis set 1 (C1 for mol1) to basis +// set 2 (C2 for mol2). + +// .. math:: + +// |\psi1\rangle = |AO1\rangle C1 + +// |\psi2\rangle = P |\psi1\rangle = |AO2\rangle S^{-1}\langle AO2| AO1\rangle> C1 = |AO2\rangle> C2 + +// C2 = S^{-1}\langle AO2|AO1\rangle C1 + +// There are three relevant functions: +// :func:`project_mo_nr2nr` is the projection for non-relativistic (scalar) basis. +// :func:`project_mo_nr2r` projects from non-relativistic to relativistic basis. +// :func:`project_mo_r2r` is the projection between relativistic (spinor) basis. +// ''' +// s22 = mol2.intor_symmetric('int1e_ovlp') +// s21 = mole.intor_cross('int1e_ovlp', mol2, mol1) +// if isinstance(mo1, numpy.ndarray) and mo1.ndim == 2: +// return lib.cho_solve(s22, numpy.dot(s21, mo1), strict_sym_pos=False) +// else: +// return [lib.cho_solve(s22, numpy.dot(s21, x), strict_sym_pos=False) +// for x in mo1] + pub fn proj_mo(mol1: &Molecule, mo1: [MatrixFull;2], cint2: (Vec>, Vec>, Vec), ecp2: Option>>) -> ([MatrixFull;2],[Vec;2],Option<[Vec;2]>) { unimplemented!() } \ No newline at end of file diff --git a/src/molecule_io/int_cross.rs b/src/molecule_io/int_cross.rs new file mode 100644 index 0000000..060d6ef --- /dev/null +++ b/src/molecule_io/int_cross.rs @@ -0,0 +1,35 @@ +use rest_libcint::prelude::*; +use rest_tensors::MatrixFull; +use crate::molecule_io::Molecule; + +impl Molecule { + /// Compute cross integrals between basis functions of two molecules: + /// ⟨ μ | op | ν ⟩, where μ ∈ self, ν ∈ other + /// + /// # Arguments + /// - `other`: the second molecule + /// - `op_name`: operator name, supports "ovlp", "kinetic", "nuclear" + /// + /// # Returns + /// - `MatrixFull` with shape [self.num_basis, other.num_basis] + pub fn int_cross(&self, other: &Molecule, op_name: String) -> MatrixFull { + let intor = match op_name.as_str() { + "ovlp" => "int1e_ovlp", + "kinetic" => "int1e_kin", + "nuclear" => "int1e_nuc", + _ => panic!("Error:: unsupported op_name for int_cross: {}. Supported: ovlp, kinetic, nuclear", op_name), + }; + + let cint_data1 = self.initialize_cint(false); + let cint_data2 = other.initialize_cint(false); + + let (out, shape) = CInt::integrate_cross_row_major( + intor, + [&cint_data1, &cint_data2], + "s1", + None, + ).into(); + + MatrixFull::from_vec([shape[0], shape[1]], out).unwrap() + } +} diff --git a/src/molecule_io/mod.rs b/src/molecule_io/mod.rs index 269c310..fb12996 100644 --- a/src/molecule_io/mod.rs +++ b/src/molecule_io/mod.rs @@ -3,6 +3,7 @@ extern crate rest_tensors as tensors; mod pyrest_molecule_io; pub mod with_clause; +mod int_cross; use array_tool::vec::Intersect; use pyo3::{pyclass}; -- Gitee From d20146db581d05de99a3d124e1ec59e188530fe7 Mon Sep 17 00:00:00 2001 From: Shirong_Wang Date: Sat, 16 May 2026 21:28:21 +0800 Subject: [PATCH 09/14] cho_solve --- src/initial_guess/mod.rs | 13 +++--- src/initial_guess/proj.rs | 87 ++++++++++++++++++++++++++++++++------- 2 files changed, 80 insertions(+), 20 deletions(-) diff --git a/src/initial_guess/mod.rs b/src/initial_guess/mod.rs index e13eb60..4ec19d8 100644 --- a/src/initial_guess/mod.rs +++ b/src/initial_guess/mod.rs @@ -228,15 +228,15 @@ pub fn update_scf_from_hdf5chk(scf_data: &mut SCF, chkfile: String) { //============================= scf_data.eigenvalues = eigenvalues; scf_data.eigenvectors = eigenvectors; - // if let Some(occupation) = is_occupation { - // scf_data.occupation = occupation; + if let Some(occupation) = is_occupation { + scf_data.occupation = occupation; // let (homo, lumo) = scf_io::util::get_homo_lumo_from_occ_integer(&scf_data.occupation, &scf_data.scftype); // scf_data.homo = homo; // scf_data.lumo = lumo; - // } else { + } else { // Since homo, lumo is needed, re-generate occupation temporarily scf_data.generate_occupation(); - // } + } println!("homo {:?}", scf_data.homo); scf_data.generate_density_matrix(); } @@ -323,7 +323,7 @@ pub fn initial_guess_from_hdf5chk(mol: &mut Molecule, scftype: &SCFType, chkfile } else { println!("The basis set in chkfile does not match the input basis set, trying to project..."); let (cint_raw_data, ecp_raw, basis4elem) = chkfile::load_cint_data(chkfile); - let (mo1, _, _) = initial_guess_from_raw( + let (mo2, _, occ2) = initial_guess_from_raw( loaded_eigenvectors, loaded_eigenvalues, loaded_occupation.unwrap(), @@ -332,7 +332,8 @@ pub fn initial_guess_from_hdf5chk(mol: &mut Molecule, scftype: &SCFType, chkfile loaded_nbasis, mol.ctrl.print_level ); - return proj::proj_mo(mol, mo1, cint_raw_data.unwrap(), ecp_raw); + let mo = proj::proj_mo(mol, mo2, cint_raw_data.unwrap(), ecp_raw); + return (mo, [vec![], vec![]], occ2); } } diff --git a/src/initial_guess/proj.rs b/src/initial_guess/proj.rs index 0f1baee..802c312 100644 --- a/src/initial_guess/proj.rs +++ b/src/initial_guess/proj.rs @@ -1,25 +1,45 @@ -// use crate::scf_io::SCF; +use rstsr::prelude::*; use tensors::matrix::MatrixFull; +use tensors::BasicMatrix; use crate::molecule_io::Molecule; +use crate::utilities::rstsr_util::RestTensorToRstsrTsrAPI; + +/// Solve S22 * X = B using Cholesky + triangular solve, with fallback to general solve. +/// Equivalent to pyscf's lib.cho_solve(a, b, strict_sym_pos=False). +fn cho_solve(s22: &MatrixFull, b: &MatrixFull) -> MatrixFull { + let device = DeviceBLAS::default(); + let s22_tsr = s22.to_rstsr(&device); + let b_tsr = b.to_rstsr(&device); + + match rt::linalg::cholesky_f((&s22_tsr, Lower)) { + Ok(l) => { + let y = rt::linalg::solve_triangular((&l, &b_tsr, Lower)); + let x = rt::linalg::solve_triangular((l.t(), &y, Upper)); + let shape: [usize; 2] = x.shape().to_vec().try_into().unwrap(); + MatrixFull::from_vec(shape, x.into_shape(-1).into_vec()).unwrap() + }, + Err(_) => { + let x = rt::linalg::solve_general((&s22_tsr, &b_tsr)); + let shape: [usize; 2] = x.shape().to_vec().try_into().unwrap(); + MatrixFull::from_vec(shape, x.into_shape(-1).into_vec()).unwrap() + } + } +} // reference implementation: // // def project_mo_nr2nr(mol1, mo1, mol2): // r''' Project orbital coefficients from basis set 1 (C1 for mol1) to basis // set 2 (C2 for mol2). - +// // .. math:: - +// // |\psi1\rangle = |AO1\rangle C1 - -// |\psi2\rangle = P |\psi1\rangle = |AO2\rangle S^{-1}\langle AO2| AO1\rangle> C1 = |AO2\rangle> C2 - +// +// |\psi2\rangle = P |\psi1\rangle = |AO2\rangle S^{-1}\langle AO2| AO1\rangle C1 = |AO2\rangle C2 +// // C2 = S^{-1}\langle AO2|AO1\rangle C1 - -// There are three relevant functions: -// :func:`project_mo_nr2nr` is the projection for non-relativistic (scalar) basis. -// :func:`project_mo_nr2r` projects from non-relativistic to relativistic basis. -// :func:`project_mo_r2r` is the projection between relativistic (spinor) basis. +// // ''' // s22 = mol2.intor_symmetric('int1e_ovlp') // s21 = mole.intor_cross('int1e_ovlp', mol2, mol1) @@ -29,6 +49,45 @@ use crate::molecule_io::Molecule; // return [lib.cho_solve(s22, numpy.dot(s21, x), strict_sym_pos=False) // for x in mo1] -pub fn proj_mo(mol1: &Molecule, mo1: [MatrixFull;2], cint2: (Vec>, Vec>, Vec), ecp2: Option>>) -> ([MatrixFull;2],[Vec;2],Option<[Vec;2]>) { - unimplemented!() -} \ No newline at end of file +/// Project MO coefficients from a source basis (given as raw cint data) to a target molecule's basis. +/// +/// # Arguments +/// - `mol_target`: the target molecule (basis set 2) +/// - `mo_source`: MO coefficients in the source basis (C1 for mol_source) +/// - `cint_source`: raw cint data for the source basis (basis set 1) +/// - `ecp_source`: optional ECP data for the source basis +/// +/// # Returns +/// - MO coefficients projected to the target basis (C2) +pub fn proj_mo(mol_target: &Molecule, mo_source: [MatrixFull;2], cint_source: (Vec>, Vec>, Vec), ecp_source: Option>>) -> [MatrixFull;2] { + let mut mol_source = Molecule::init_mol(); + mol_source.cint_type = mol_target.cint_type.clone(); + mol_source.ctrl.print_level = mol_target.ctrl.print_level; + mol_source.update_from_cint(0, None, cint_source, ecp_source); + + // S22 = target self-overlap + let s22_full = mol_target.int_ij_matrixupper("ovlp".to_string()).to_matrixfull().unwrap(); + + // S21 = + let s21 = mol_target.int_cross(&mol_source, "ovlp".to_string()); + + let device = DeviceBLAS::default(); + let s21_tsr = s21.to_rstsr(&device); + + let mut mo_target: [MatrixFull; 2] = [MatrixFull::empty(), MatrixFull::empty()]; + for spin in 0..2 { + if mo_source[spin].size()[0] == 0 { + continue; + } + let mo_source_tsr = mo_source[spin].to_rstsr(&device); + let temp_tsr = &s21_tsr % &mo_source_tsr; + let temp = MatrixFull::from_vec( + temp_tsr.shape().to_vec().try_into().unwrap(), + temp_tsr.into_shape(-1).into_vec(), + ).unwrap(); + + mo_target[spin] = cho_solve(&s22_full, &temp); + } + + mo_target +} -- Gitee From 6731ad5892b9a3b9f3a58d667669be40b564e91f Mon Sep 17 00:00:00 2001 From: Shirong_Wang Date: Sat, 16 May 2026 21:40:54 +0800 Subject: [PATCH 10/14] normalize mo --- src/initial_guess/proj.rs | 37 +++++++++++++++++++++++++++++-------- 1 file changed, 29 insertions(+), 8 deletions(-) diff --git a/src/initial_guess/proj.rs b/src/initial_guess/proj.rs index 802c312..cdb773b 100644 --- a/src/initial_guess/proj.rs +++ b/src/initial_guess/proj.rs @@ -2,14 +2,13 @@ use rstsr::prelude::*; use tensors::matrix::MatrixFull; use tensors::BasicMatrix; use crate::molecule_io::Molecule; -use crate::utilities::rstsr_util::RestTensorToRstsrTsrAPI; +use crate::utilities::rstsr_util::{RestTensorToRstsrTsrAPI, RestTensorToRstsrViewAPI}; /// Solve S22 * X = B using Cholesky + triangular solve, with fallback to general solve. /// Equivalent to pyscf's lib.cho_solve(a, b, strict_sym_pos=False). -fn cho_solve(s22: &MatrixFull, b: &MatrixFull) -> MatrixFull { - let device = DeviceBLAS::default(); - let s22_tsr = s22.to_rstsr(&device); - let b_tsr = b.to_rstsr(&device); +fn cho_solve(s22: &MatrixFull, b: &MatrixFull, device: &DeviceBLAS) -> MatrixFull { + let s22_tsr = s22.to_rstsr_view(device); + let b_tsr = b.to_rstsr_view(device); match rt::linalg::cholesky_f((&s22_tsr, Lower)) { Ok(l) => { @@ -26,6 +25,27 @@ fn cho_solve(s22: &MatrixFull, b: &MatrixFull) -> MatrixFull { } } +/// Normalize MO columns so that mo^T S mo = I. +/// Equivalent to pyscf: +/// norm = numpy.einsum('pi,pi->i', mo.conj(), s.dot(mo)) +/// mo /= numpy.sqrt(norm) +fn normalize_mo(mo: &mut MatrixFull, s: &MatrixFull, device: &DeviceBLAS) { + let s_tsr = s.to_rstsr_view(device); + let mo_tsr = mo.to_rstsr_view(device); + + let s_dot_mo = &s_tsr % &mo_tsr; + let norm = (&mo_tsr * &s_dot_mo).sum_axes(0); + + let norm_vec: Vec = norm.to_vec(); + let inv_sqrt_norm: Vec = norm_vec.iter().map(|x| 1.0 / x.sqrt()).collect(); + let inv_sqrt_tsr = rt::asarray((inv_sqrt_norm, device)); + + let mo_normalized = &mo_tsr * inv_sqrt_tsr; + + let shape: [usize; 2] = mo_normalized.shape().to_vec().try_into().unwrap(); + *mo = MatrixFull::from_vec(shape, mo_normalized.into_shape(-1).into_vec()).unwrap(); +} + // reference implementation: // // def project_mo_nr2nr(mol1, mo1, mol2): @@ -58,7 +78,7 @@ fn cho_solve(s22: &MatrixFull, b: &MatrixFull) -> MatrixFull { /// - `ecp_source`: optional ECP data for the source basis /// /// # Returns -/// - MO coefficients projected to the target basis (C2) +/// - MO coefficients projected to the target basis (C2), normalized so that C2^T S22 C2 = I pub fn proj_mo(mol_target: &Molecule, mo_source: [MatrixFull;2], cint_source: (Vec>, Vec>, Vec), ecp_source: Option>>) -> [MatrixFull;2] { let mut mol_source = Molecule::init_mol(); mol_source.cint_type = mol_target.cint_type.clone(); @@ -79,14 +99,15 @@ pub fn proj_mo(mol_target: &Molecule, mo_source: [MatrixFull;2], cint_sourc if mo_source[spin].size()[0] == 0 { continue; } - let mo_source_tsr = mo_source[spin].to_rstsr(&device); + let mo_source_tsr = mo_source[spin].to_rstsr_view(&device); let temp_tsr = &s21_tsr % &mo_source_tsr; let temp = MatrixFull::from_vec( temp_tsr.shape().to_vec().try_into().unwrap(), temp_tsr.into_shape(-1).into_vec(), ).unwrap(); - mo_target[spin] = cho_solve(&s22_full, &temp); + mo_target[spin] = cho_solve(&s22_full, &temp, &device); + normalize_mo(&mut mo_target[spin], &s22_full, &device); } mo_target -- Gitee From f7aa96d47994710ef615da873260be77aea97bb6 Mon Sep 17 00:00:00 2001 From: Shirong_Wang Date: Sun, 17 May 2026 15:12:30 +0800 Subject: [PATCH 11/14] debug --- src/initial_guess/mod.rs | 5 +++-- src/initial_guess/proj.rs | 9 +++++++-- src/scf_io/mod.rs | 8 +++++--- 3 files changed, 15 insertions(+), 7 deletions(-) diff --git a/src/initial_guess/mod.rs b/src/initial_guess/mod.rs index 4ec19d8..97de43f 100644 --- a/src/initial_guess/mod.rs +++ b/src/initial_guess/mod.rs @@ -237,7 +237,8 @@ pub fn update_scf_from_hdf5chk(scf_data: &mut SCF, chkfile: String) { // Since homo, lumo is needed, re-generate occupation temporarily scf_data.generate_occupation(); } - println!("homo {:?}", scf_data.homo); + // println!("homo {:?}", scf_data.homo); + println!("occ {:?}", scf_data.occupation); scf_data.generate_density_matrix(); } @@ -321,7 +322,7 @@ pub fn initial_guess_from_hdf5chk(mol: &mut Molecule, scftype: &SCFType, chkfile mol.ctrl.print_level ); } else { - println!("The basis set in chkfile does not match the input basis set, trying to project..."); + println!("The basis set in chkfile does not match the input basis set (loaded: {} basis, {} MOs; target: {} basis, {} MOs), \ntrying to project...", loaded_nbasis, loaded_nmo, mol.num_basis, mol.num_state); let (cint_raw_data, ecp_raw, basis4elem) = chkfile::load_cint_data(chkfile); let (mo2, _, occ2) = initial_guess_from_raw( loaded_eigenvectors, diff --git a/src/initial_guess/proj.rs b/src/initial_guess/proj.rs index cdb773b..aea8eef 100644 --- a/src/initial_guess/proj.rs +++ b/src/initial_guess/proj.rs @@ -12,12 +12,14 @@ fn cho_solve(s22: &MatrixFull, b: &MatrixFull, device: &DeviceBLAS) -> match rt::linalg::cholesky_f((&s22_tsr, Lower)) { Ok(l) => { + println!("Cholesky decomposition succeeded, using triangular solve."); let y = rt::linalg::solve_triangular((&l, &b_tsr, Lower)); let x = rt::linalg::solve_triangular((l.t(), &y, Upper)); let shape: [usize; 2] = x.shape().to_vec().try_into().unwrap(); MatrixFull::from_vec(shape, x.into_shape(-1).into_vec()).unwrap() }, Err(_) => { + println!("Cholesky decomposition failed, falling back to general solve."); let x = rt::linalg::solve_general((&s22_tsr, &b_tsr)); let shape: [usize; 2] = x.shape().to_vec().try_into().unwrap(); MatrixFull::from_vec(shape, x.into_shape(-1).into_vec()).unwrap() @@ -38,9 +40,11 @@ fn normalize_mo(mo: &mut MatrixFull, s: &MatrixFull, device: &DeviceBL let norm_vec: Vec = norm.to_vec(); let inv_sqrt_norm: Vec = norm_vec.iter().map(|x| 1.0 / x.sqrt()).collect(); - let inv_sqrt_tsr = rt::asarray((inv_sqrt_norm, device)); + let inv_sqrt_tsr = rt::asarray((inv_sqrt_norm.clone(), device)); - let mo_normalized = &mo_tsr * inv_sqrt_tsr; + println!("mo shape: {:?}, norm shape: {:?}, inv_sqrt_norm shape: {:?}", mo_tsr.shape(), norm.shape(), inv_sqrt_tsr.shape()); + println!("inv_sqrt_norm: {:?}", inv_sqrt_norm[..5].to_vec()); + let mo_normalized = &mo_tsr * inv_sqrt_tsr.slice((None, ..)); let shape: [usize; 2] = mo_normalized.shape().to_vec().try_into().unwrap(); *mo = MatrixFull::from_vec(shape, mo_normalized.into_shape(-1).into_vec()).unwrap(); @@ -109,6 +113,7 @@ pub fn proj_mo(mol_target: &Molecule, mo_source: [MatrixFull;2], cint_sourc mo_target[spin] = cho_solve(&s22_full, &temp, &device); normalize_mo(&mut mo_target[spin], &s22_full, &device); } + mo_target[0].formated_output(5, "full"); mo_target } diff --git a/src/scf_io/mod.rs b/src/scf_io/mod.rs index c90dc22..30fa1ca 100644 --- a/src/scf_io/mod.rs +++ b/src/scf_io/mod.rs @@ -5337,8 +5337,9 @@ pub fn generate_occupation_outside(scf_data: &SCF) -> ([Vec;2], [usize;2], pub fn generate_density_matrix_outside(scf_data: &SCF) -> Vec>{ - let num_basis = scf_data.mol.num_basis; - let num_state = scf_data.mol.num_state; + // let num_basis = scf_data.mol.num_basis; + // let num_state = scf_data.mol.num_state; + let [num_basis, num_state] = scf_data.eigenvectors[0].size(); let spin_channel = scf_data.mol.spin_channel; // let homo = &scf_data.homo; // println!("homo: {:?}", &homo); @@ -5357,7 +5358,7 @@ pub fn generate_density_matrix_outside(scf_data: &SCF) -> Vec>{ let occ_s = &scf_data.occupation[i_spin]; let nw = occupied_orbital_count(&scf_data.occupation[i_spin]); - //println!("number of occupied orbitals from dm generation: {}", nw); + println!("number of occupied orbitals from dm generation: {}", nw); let mut weight_eigv = MatrixFull::new([num_basis, num_state],0.0_f64); //let mut weight_eigv = eigv_s.clone(); @@ -5454,6 +5455,7 @@ pub fn scf_without_build(scf_data: &mut SCF, mpi_operator: &Option) scf_data.generate_hf_hamiltonian(mpi_operator); let mut scf_records=ScfTraceRecord::initialize(&scf_data); + println!("dm {:?}", scf_data.density_matrix); if scf_data.mol.ctrl.print_level>0 {println!("The total energy: {:20.10} Ha by the initial guess",scf_data.scf_energy)}; //let mut scf_continue = true; -- Gitee From 6ae700527130b8d7412782b3e3dc7feb4548c68c Mon Sep 17 00:00:00 2001 From: Shirong_Wang Date: Sun, 17 May 2026 17:39:29 +0800 Subject: [PATCH 12/14] fix basisproj --- src/initial_guess/mod.rs | 4 ++-- src/initial_guess/proj.rs | 6 +++--- src/molecule_io/int_cross.rs | 2 +- src/scf_io/mod.rs | 3 +-- 4 files changed, 7 insertions(+), 8 deletions(-) diff --git a/src/initial_guess/mod.rs b/src/initial_guess/mod.rs index 97de43f..29ecc79 100644 --- a/src/initial_guess/mod.rs +++ b/src/initial_guess/mod.rs @@ -2,7 +2,7 @@ use tensors::{MatrixFull, MatrixUpper}; use crate::initial_guess::enxc::effective_nxc_matrix; use crate::mpi_io::{mpi_broadcast_matrixfull, MPIOperator}; -use crate::scf_io::{self, SCF, SCFType}; +use crate::scf_io::{SCF, SCFType}; use crate::{molecule_io::Molecule, dft::Grids}; use crate::initial_guess::sap::get_vsap; use self::sad::initial_guess_from_sad; @@ -238,7 +238,7 @@ pub fn update_scf_from_hdf5chk(scf_data: &mut SCF, chkfile: String) { scf_data.generate_occupation(); } // println!("homo {:?}", scf_data.homo); - println!("occ {:?}", scf_data.occupation); + // println!("occ {:?}", scf_data.occupation); scf_data.generate_density_matrix(); } diff --git a/src/initial_guess/proj.rs b/src/initial_guess/proj.rs index aea8eef..2f1d3a1 100644 --- a/src/initial_guess/proj.rs +++ b/src/initial_guess/proj.rs @@ -42,8 +42,8 @@ fn normalize_mo(mo: &mut MatrixFull, s: &MatrixFull, device: &DeviceBL let inv_sqrt_norm: Vec = norm_vec.iter().map(|x| 1.0 / x.sqrt()).collect(); let inv_sqrt_tsr = rt::asarray((inv_sqrt_norm.clone(), device)); - println!("mo shape: {:?}, norm shape: {:?}, inv_sqrt_norm shape: {:?}", mo_tsr.shape(), norm.shape(), inv_sqrt_tsr.shape()); - println!("inv_sqrt_norm: {:?}", inv_sqrt_norm[..5].to_vec()); + // println!("mo shape: {:?}, norm shape: {:?}, inv_sqrt_norm shape: {:?}", mo_tsr.shape(), norm.shape(), inv_sqrt_tsr.shape()); + // println!("inv_sqrt_norm: {:?}", inv_sqrt_norm[..5].to_vec()); let mo_normalized = &mo_tsr * inv_sqrt_tsr.slice((None, ..)); let shape: [usize; 2] = mo_normalized.shape().to_vec().try_into().unwrap(); @@ -113,7 +113,7 @@ pub fn proj_mo(mol_target: &Molecule, mo_source: [MatrixFull;2], cint_sourc mo_target[spin] = cho_solve(&s22_full, &temp, &device); normalize_mo(&mut mo_target[spin], &s22_full, &device); } - mo_target[0].formated_output(5, "full"); + // mo_target[0].formated_output(5, "full"); mo_target } diff --git a/src/molecule_io/int_cross.rs b/src/molecule_io/int_cross.rs index 060d6ef..db7757d 100644 --- a/src/molecule_io/int_cross.rs +++ b/src/molecule_io/int_cross.rs @@ -23,7 +23,7 @@ impl Molecule { let cint_data1 = self.initialize_cint(false); let cint_data2 = other.initialize_cint(false); - let (out, shape) = CInt::integrate_cross_row_major( + let (out, shape) = CInt::integrate_cross( intor, [&cint_data1, &cint_data2], "s1", diff --git a/src/scf_io/mod.rs b/src/scf_io/mod.rs index 30fa1ca..fbccfb0 100644 --- a/src/scf_io/mod.rs +++ b/src/scf_io/mod.rs @@ -5358,7 +5358,7 @@ pub fn generate_density_matrix_outside(scf_data: &SCF) -> Vec>{ let occ_s = &scf_data.occupation[i_spin]; let nw = occupied_orbital_count(&scf_data.occupation[i_spin]); - println!("number of occupied orbitals from dm generation: {}", nw); + // println!("number of occupied orbitals from dm generation: {}", nw); let mut weight_eigv = MatrixFull::new([num_basis, num_state],0.0_f64); //let mut weight_eigv = eigv_s.clone(); @@ -5455,7 +5455,6 @@ pub fn scf_without_build(scf_data: &mut SCF, mpi_operator: &Option) scf_data.generate_hf_hamiltonian(mpi_operator); let mut scf_records=ScfTraceRecord::initialize(&scf_data); - println!("dm {:?}", scf_data.density_matrix); if scf_data.mol.ctrl.print_level>0 {println!("The total energy: {:20.10} Ha by the initial guess",scf_data.scf_energy)}; //let mut scf_continue = true; -- Gitee From 6a8fde439675aea2af7bd283cc05e91c90f5ba96 Mon Sep 17 00:00:00 2001 From: Shirong_Wang Date: Sun, 17 May 2026 18:53:01 +0800 Subject: [PATCH 13/14] adjust stru --- src/fileop/chkfile.rs | 6 +- src/fileop/mod.rs | 1 + src/initial_guess/proj.rs | 23 +----- src/molecule_io/frozen.rs | 150 +++++++++++++++++++++++++++++++++++++ src/molecule_io/mod.rs | 151 ++------------------------------------ 5 files changed, 160 insertions(+), 171 deletions(-) create mode 100644 src/molecule_io/frozen.rs diff --git a/src/fileop/chkfile.rs b/src/fileop/chkfile.rs index e2e15d3..6955179 100644 --- a/src/fileop/chkfile.rs +++ b/src/fileop/chkfile.rs @@ -2,7 +2,7 @@ use std::collections::HashMap; use std::path::Path; use hdf5; use hdf5::types::VarLenUnicode; -use hdf5::types::TypeDescriptor; +// use hdf5::types::TypeDescriptor; use crate::scf_io::{SCF, SCFType}; use crate::geom_io::get_mass_charge; use crate::basis_io::Basis4Elem; @@ -127,7 +127,7 @@ pub fn save_hamiltonian(scf_data: &SCF) { }; let mut hamiltonians: Vec = vec![]; for i_spin in 0..scf_data.mol.spin_channel { - let tmp_eigenvectors = scf_data.hamiltonian[i_spin].to_matrixfull().unwrap(); + // let tmp_eigenvectors = scf_data.hamiltonian[i_spin].to_matrixfull().unwrap(); hamiltonians.extend(scf_data.eigenvalues[i_spin].iter()); } let is_hamiltonian = scf.member_names().unwrap().iter().fold(false,|is_exist,x| {is_exist || x.eq("hamiltonian")}); @@ -155,7 +155,7 @@ pub fn save_overlap(scf_data: &SCF) { } else { file.create_group("scf").unwrap() }; - let mut overlap: Vec = vec![]; + // let mut overlap: Vec = vec![]; let overlap = scf_data.ovlp.to_matrixfull().unwrap(); let is_exist = scf.member_names().unwrap().iter().fold(false,|is_exist,x| {is_exist || x.eq("overlap")}); if is_exist { diff --git a/src/fileop/mod.rs b/src/fileop/mod.rs index 0983024..c5a319f 100644 --- a/src/fileop/mod.rs +++ b/src/fileop/mod.rs @@ -1 +1,2 @@ +#![warn(unused_imports)] pub mod chkfile; \ No newline at end of file diff --git a/src/initial_guess/proj.rs b/src/initial_guess/proj.rs index 2f1d3a1..e846a35 100644 --- a/src/initial_guess/proj.rs +++ b/src/initial_guess/proj.rs @@ -50,28 +50,7 @@ fn normalize_mo(mo: &mut MatrixFull, s: &MatrixFull, device: &DeviceBL *mo = MatrixFull::from_vec(shape, mo_normalized.into_shape(-1).into_vec()).unwrap(); } -// reference implementation: -// -// def project_mo_nr2nr(mol1, mo1, mol2): -// r''' Project orbital coefficients from basis set 1 (C1 for mol1) to basis -// set 2 (C2 for mol2). -// -// .. math:: -// -// |\psi1\rangle = |AO1\rangle C1 -// -// |\psi2\rangle = P |\psi1\rangle = |AO2\rangle S^{-1}\langle AO2| AO1\rangle C1 = |AO2\rangle C2 -// -// C2 = S^{-1}\langle AO2|AO1\rangle C1 -// -// ''' -// s22 = mol2.intor_symmetric('int1e_ovlp') -// s21 = mole.intor_cross('int1e_ovlp', mol2, mol1) -// if isinstance(mo1, numpy.ndarray) and mo1.ndim == 2: -// return lib.cho_solve(s22, numpy.dot(s21, mo1), strict_sym_pos=False) -// else: -// return [lib.cho_solve(s22, numpy.dot(s21, x), strict_sym_pos=False) -// for x in mo1] + /// Project MO coefficients from a source basis (given as raw cint data) to a target molecule's basis. /// diff --git a/src/molecule_io/frozen.rs b/src/molecule_io/frozen.rs new file mode 100644 index 0000000..5278c9c --- /dev/null +++ b/src/molecule_io/frozen.rs @@ -0,0 +1,150 @@ +use crate::molecule_io::Molecule; +use crate::constants::{ELEM1ST, ELEM2ND, ELEM3RD, ELEM4TH, ELEM5TH, ELEM6TH, ELEMTMS}; + +impl Molecule { + pub fn generate_start_mo(&self, ecp_electrons: usize) -> usize { + let mut start_mo = count_frozen_core_states(self.ctrl.frozen_core_postscf, &self.geom.elem); + let ecp_orbs = ecp_electrons/2; + if ecp_orbs == 0 { + if self.ctrl.print_level > 0 { + println!("For SCF calculation, no core orbital is frozen by effective core potential (ECP) approximation"); + } + } else if start_mo < ecp_orbs { + if self.ctrl.print_level > 0 { + println!(" for the frozen-core post-SCF methods {} is smaller than the number of ecp orbitals {}. Set = 0", + start_mo, ecp_orbs); + } + start_mo = 0; + } else { + if self.ctrl.print_level > 0 { + println!(" for the frozen-core post-SCF methods {} is larger than the number of ecp orbitals {}. Take the -= ecp_orbitals", + start_mo, ecp_orbs); + } + start_mo = start_mo - ecp_orbs + } + start_mo + } +} + + +pub fn count_frozen_core_states(n_frozen_shell: i32, elem: &Vec) -> usize { + let mut n_low_state = 0_usize; + let mut n_tm = 0_usize; + + //let n_frozen_shell = self.ctrl.frozen_core_postscf; + let (n_frozen_shell_1, n_frozen_shell_2) = if n_frozen_shell > 10 { + let n_frozen_shell_1 = n_frozen_shell%10; + let n_frozen_shell_2 = n_frozen_shell/10; + (n_frozen_shell_1, n_frozen_shell_2) + } else { + (n_frozen_shell, n_frozen_shell) + }; + + elem.iter().for_each(|sn| { + let formated_elem = crate::geom_io::formated_element_name(&sn); + let flag_first_row = ELEM1ST.iter().fold(false, |acc, elem| acc || elem.eq(&formated_elem)); + let flag_second_row = if flag_first_row { + false + } else { + ELEM2ND.iter().fold(false, |acc, elem| acc || elem.eq(&formated_elem)) + }; + let flag_third_row = if flag_first_row || flag_second_row { + false + } else { + ELEM3RD.iter().fold(false, |acc, elem| acc || elem.eq(&formated_elem)) + }; + let flag_fourth_row = if flag_first_row || flag_second_row || flag_third_row { + false + } else { + ELEM4TH.iter().fold(false, |acc, elem| acc || elem.eq(&formated_elem)) + }; + let flag_fifth_row = if flag_first_row || flag_second_row || flag_third_row || flag_fourth_row { + false + } else { + ELEM5TH.iter().fold(false, |acc, elem| acc || elem.eq(&formated_elem)) + }; + let flag_sixth_row = if flag_first_row || flag_second_row || flag_third_row || flag_fourth_row || flag_fifth_row { + false + } else { + ELEM6TH.iter().fold(false, |acc, elem| acc || elem.eq(&formated_elem)) + }; + let flag_tm = ELEMTMS.iter().fold(false, |acc, elem| acc || elem.eq(&formated_elem)); + + let n_frozen_shell_curr = if flag_tm { + n_tm += 1; + n_frozen_shell_2 + } else { + n_frozen_shell_1 + }; + + if flag_first_row { + n_low_state += 0 + } else if flag_second_row { + if n_frozen_shell_curr==0 { + n_low_state += 0 + } else if (n_frozen_shell_curr==1) { + n_low_state += 1 + } else { + n_low_state += 0 + } + } else if flag_third_row { + if n_frozen_shell_curr==0 { + n_low_state += 0 + } else if n_frozen_shell_curr==1 { + n_low_state += 5 + } else if n_frozen_shell_curr==2 { + n_low_state += 1 + } else { + n_low_state += 0 + } + } else if flag_fourth_row { + if n_frozen_shell_curr==0 { + n_low_state += 0 + } else if n_frozen_shell_curr==1 { + n_low_state += 9 + } else if n_frozen_shell_curr==2 { + n_low_state += 5 + } else if n_frozen_shell_curr==3 { + n_low_state += 1 + } else { + n_low_state += 0 + } + } else if flag_fifth_row { + if n_frozen_shell_curr==0 { + n_low_state += 0 + } else if n_frozen_shell_curr==1 { + n_low_state += 18 + } else if n_frozen_shell_curr==2 { + // NOTE: for 4d-block elements, 3d orbitals are frozen as well for n_frozen_shell = 2 + if flag_tm {n_low_state += 14} else {n_low_state += 9} + } else if n_frozen_shell_curr==3 { + n_low_state += 5 + } else if n_frozen_shell_curr==4 { + n_low_state += 1 + } else { + n_low_state += 0 + } + } else if flag_sixth_row { + if n_frozen_shell_curr==0 { + n_low_state += 0 + } else if n_frozen_shell_curr==1 { + n_low_state += 34 + } else if n_frozen_shell_curr==2 { + // NOTE: for 5d-block elements, 4d and 4f orbitals are frozen as well for n_frozen_shell = 2 + // It leads to a core shell with 60 electrons + if flag_tm {n_low_state += 30} else {n_low_state += 18} + } else if n_frozen_shell_curr==3 { + n_low_state += 9 + } else if n_frozen_shell_curr==4 { + n_low_state += 5 + } else if n_frozen_shell_curr==5 { + n_low_state += 1 + } else { + n_low_state += 0 + } + }; + }); + + n_low_state + +} \ No newline at end of file diff --git a/src/molecule_io/mod.rs b/src/molecule_io/mod.rs index fb12996..537e484 100644 --- a/src/molecule_io/mod.rs +++ b/src/molecule_io/mod.rs @@ -4,22 +4,25 @@ extern crate rest_tensors as tensors; mod pyrest_molecule_io; pub mod with_clause; mod int_cross; +mod frozen; +pub use frozen::count_frozen_core_states; use array_tool::vec::Intersect; use pyo3::{pyclass}; use rayon::prelude::{IntoParallelRefIterator, IndexedParallelIterator, ParallelIterator}; use rest_libcint::prelude::*; +use rest_libcint::{CINTR2CDATA, CintType}; use rest_tensors::{ERIFull,RIFull,ERIFold4,TensorSlice,TensorSliceMut,TensorOpt, MatrixUpper, MatrixFull}; use tensors::{BasicMatrix, SubMatrixUpper}; use tensors::external_libs::{matr_copy_from_ri}; use tensors::matrix_blas_lapack::{_dgemm, _dgemm_full, _power_rayon_for_symmetric_matrix}; +use tensors::matrix_blas_lapack::{omp_set_num_threads_wrapper}; use std::collections::HashMap; use std::ops::Range; use std::sync::mpsc::channel; -use rest_libcint::{CINTR2CDATA, CintType}; use regex::Regex; use crate::basis_io::etb::{get_etb_elem, etb_gen_for_atom_list, InfoV2}; -use crate::constants::{ATM_NUC, ATM_NUC_MOD_OF, AUXBAS_THRESHOLD, ELEM1ST, ELEM2ND, ELEM3RD, ELEM4TH, ELEM5TH, ELEM6TH, ELEMTMS, ENV_PRT_START, NUC_ECP, NUC_STAD_CHARGE}; +use crate::constants::{ATM_NUC, ATM_NUC_MOD_OF, AUXBAS_THRESHOLD, ENV_PRT_START, NUC_ECP, NUC_STAD_CHARGE}; use crate::dft::{DFTType, DFA4REST, parse_xc}; use crate::geom_io::{GeomCell, get_mass_charge, formated_element_name}; use crate::basis_io::{BasInfo, Basis4Elem}; @@ -28,7 +31,6 @@ use crate::mpi_io::{mpi_isend_irecv_wrt_distribution_v03, MPIData, MPIOperator}; use crate::utilities; use crate::basis_io::bse_downloader::{self, ctrl_element_checker, local_element_checker}; use crate::basis_io::basis_list::{basis_fuzzy_matcher, check_basis_name}; -use tensors::matrix_blas_lapack::{omp_set_num_threads_wrapper}; use crate::solvent::PcmMethod; use crate::ri_jk; @@ -126,28 +128,6 @@ pub struct Molecule { } impl Molecule { - fn generate_start_mo(&self, ecp_electrons: usize) -> usize { - let mut start_mo = count_frozen_core_states(self.ctrl.frozen_core_postscf, &self.geom.elem); - let ecp_orbs = ecp_electrons/2; - if ecp_orbs == 0 { - if self.ctrl.print_level > 0 { - println!("For SCF calculation, no core orbital is frozen by effective core potential (ECP) approximation"); - } - } else if start_mo < ecp_orbs { - if self.ctrl.print_level > 0 { - println!(" for the frozen-core post-SCF methods {} is smaller than the number of ecp orbitals {}. Set = 0", - start_mo, ecp_orbs); - } - start_mo = 0; - } else { - if self.ctrl.print_level > 0 { - println!(" for the frozen-core post-SCF methods {} is larger than the number of ecp orbitals {}. Take the -= ecp_orbitals", - start_mo, ecp_orbs); - } - start_mo = start_mo - ecp_orbs - } - start_mo - } pub fn init_mol() -> Molecule { Molecule { ctrl:InputKeywords::init_ctrl(), @@ -3483,127 +3463,6 @@ pub fn generate_ri3fn_from_rimatr(rimatr: &MatrixFull, basbas2baspar: &Matr ri3fn } -pub fn count_frozen_core_states(n_frozen_shell: i32, elem: &Vec) -> usize { - let mut n_low_state = 0_usize; - let mut n_tm = 0_usize; - - //let n_frozen_shell = self.ctrl.frozen_core_postscf; - let (n_frozen_shell_1, n_frozen_shell_2) = if n_frozen_shell > 10 { - let n_frozen_shell_1 = n_frozen_shell%10; - let n_frozen_shell_2 = n_frozen_shell/10; - (n_frozen_shell_1, n_frozen_shell_2) - } else { - (n_frozen_shell, n_frozen_shell) - }; - - elem.iter().for_each(|sn| { - let formated_elem = crate::geom_io::formated_element_name(&sn); - let flag_first_row = ELEM1ST.iter().fold(false, |acc, elem| acc || elem.eq(&formated_elem)); - let flag_second_row = if flag_first_row { - false - } else { - ELEM2ND.iter().fold(false, |acc, elem| acc || elem.eq(&formated_elem)) - }; - let flag_third_row = if flag_first_row || flag_second_row { - false - } else { - ELEM3RD.iter().fold(false, |acc, elem| acc || elem.eq(&formated_elem)) - }; - let flag_fourth_row = if flag_first_row || flag_second_row || flag_third_row { - false - } else { - ELEM4TH.iter().fold(false, |acc, elem| acc || elem.eq(&formated_elem)) - }; - let flag_fifth_row = if flag_first_row || flag_second_row || flag_third_row || flag_fourth_row { - false - } else { - ELEM5TH.iter().fold(false, |acc, elem| acc || elem.eq(&formated_elem)) - }; - let flag_sixth_row = if flag_first_row || flag_second_row || flag_third_row || flag_fourth_row || flag_fifth_row { - false - } else { - ELEM6TH.iter().fold(false, |acc, elem| acc || elem.eq(&formated_elem)) - }; - let flag_tm = ELEMTMS.iter().fold(false, |acc, elem| acc || elem.eq(&formated_elem)); - - let n_frozen_shell_curr = if flag_tm { - n_tm += 1; - n_frozen_shell_2 - } else { - n_frozen_shell_1 - }; - - if flag_first_row { - n_low_state += 0 - } else if flag_second_row { - if n_frozen_shell_curr==0 { - n_low_state += 0 - } else if (n_frozen_shell_curr==1) { - n_low_state += 1 - } else { - n_low_state += 0 - } - } else if flag_third_row { - if n_frozen_shell_curr==0 { - n_low_state += 0 - } else if n_frozen_shell_curr==1 { - n_low_state += 5 - } else if n_frozen_shell_curr==2 { - n_low_state += 1 - } else { - n_low_state += 0 - } - } else if flag_fourth_row { - if n_frozen_shell_curr==0 { - n_low_state += 0 - } else if n_frozen_shell_curr==1 { - n_low_state += 9 - } else if n_frozen_shell_curr==2 { - n_low_state += 5 - } else if n_frozen_shell_curr==3 { - n_low_state += 1 - } else { - n_low_state += 0 - } - } else if flag_fifth_row { - if n_frozen_shell_curr==0 { - n_low_state += 0 - } else if n_frozen_shell_curr==1 { - n_low_state += 18 - } else if n_frozen_shell_curr==2 { - // NOTE: for 4d-block elements, 3d orbitals are frozen as well for n_frozen_shell = 2 - if flag_tm {n_low_state += 14} else {n_low_state += 9} - } else if n_frozen_shell_curr==3 { - n_low_state += 5 - } else if n_frozen_shell_curr==4 { - n_low_state += 1 - } else { - n_low_state += 0 - } - } else if flag_sixth_row { - if n_frozen_shell_curr==0 { - n_low_state += 0 - } else if n_frozen_shell_curr==1 { - n_low_state += 34 - } else if n_frozen_shell_curr==2 { - // NOTE: for 5d-block elements, 4d and 4f orbitals are frozen as well for n_frozen_shell = 2 - // It leads to a core shell with 60 electrons - if flag_tm {n_low_state += 30} else {n_low_state += 18} - } else if n_frozen_shell_curr==3 { - n_low_state += 9 - } else if n_frozen_shell_curr==4 { - n_low_state += 5 - } else if n_frozen_shell_curr==5 { - n_low_state += 1 - } else { - n_low_state += 0 - } - }; - }); - - n_low_state - -} //pub fn basbas2baspar(index: [usize;2]) -> usize { // (index[1]+1)*index[1]/2 + index[0] -- Gitee From 75baa03020c8ac78a24d8f8a1bc4052e903e9666 Mon Sep 17 00:00:00 2001 From: Shirong_Wang Date: Sun, 17 May 2026 19:56:02 +0800 Subject: [PATCH 14/14] check sanity for basis proj --- src/fileop/chkfile.rs | 21 +++++++++++++++++++-- src/initial_guess/mod.rs | 16 +++++++++++++--- src/initial_guess/proj.rs | 35 ++++++++++++++++++++++------------- src/molecule_io/mod.rs | 4 ++++ 4 files changed, 58 insertions(+), 18 deletions(-) diff --git a/src/fileop/chkfile.rs b/src/fileop/chkfile.rs index 6955179..5ef6033 100644 --- a/src/fileop/chkfile.rs +++ b/src/fileop/chkfile.rs @@ -6,6 +6,7 @@ use hdf5::types::VarLenUnicode; use crate::scf_io::{SCF, SCFType}; use crate::geom_io::get_mass_charge; use crate::basis_io::Basis4Elem; +use rest_libcint::CintType; pub fn write_scf_attribute(group: &hdf5::Group, dataset_name: &str, value: &[T]) where @@ -107,10 +108,20 @@ pub fn save_chkfile(scf_data: &SCF) { write_string_scalar(&file, "mol", &mol_info); let basis4elem = serde_json::to_string(&mol.basis4elem).unwrap(); write_string_scalar(&file, "molecule/basis4elem", &basis4elem); + let cinttype = cint_type_as_str(&mol.cint_type); + write_string_scalar(&file, "molecule/cinttype", &cinttype); file.close(); } +pub fn cint_type_as_str(cint_type: &CintType) -> &'static str { + match cint_type { + CintType::Spheric => "spheric", + CintType::Cartesian => "cartesian", + CintType::Spinor => "spinor", + } +} + pub fn save_hamiltonian(scf_data: &SCF) { let chkfile= &scf_data.mol.ctrl.chkfile; let path = Path::new(chkfile); @@ -220,7 +231,7 @@ pub fn has_mo_coeff(file: &hdf5::File) -> bool { has_mo } -pub fn load_cint_data(chkfile: &String) -> (Option<(Vec>, Vec>, Vec)>, Option>>, Option>) { +pub fn load_cint_data(chkfile: &String) -> (Option<(Vec>, Vec>, Vec)>, Option>>, Option>, Option) { let file = hdf5::File::open(chkfile).unwrap(); let mol_info = file.dataset("mol").unwrap().read_scalar::().unwrap(); let json_string = mol_info.as_str(); @@ -245,7 +256,13 @@ pub fn load_cint_data(chkfile: &String) -> (Option<(Vec>, Vec> None }; - (Some((atm, bas, env)), ecpbas, basis4elem) + let cinttype: Option = if let Some(cinttype_value) = file.dataset("molecule/cinttype").unwrap().read_scalar::().ok() { + Some(cinttype_value.as_str().into()) + } else { + None + }; + + (Some((atm, bas, env)), ecpbas, basis4elem, cinttype) } pub fn load_shapes(chkfile: &String) -> Option<(usize, usize, usize)> { diff --git a/src/initial_guess/mod.rs b/src/initial_guess/mod.rs index 29ecc79..58bea70 100644 --- a/src/initial_guess/mod.rs +++ b/src/initial_guess/mod.rs @@ -1,5 +1,7 @@ #![warn(unused_imports)] +use rest_libcint::cint; use tensors::{MatrixFull, MatrixUpper}; +use crate::constants::CINT_CAR_BAS_INFO_L2; use crate::initial_guess::enxc::effective_nxc_matrix; use crate::mpi_io::{mpi_broadcast_matrixfull, MPIOperator}; use crate::scf_io::{SCF, SCFType}; @@ -282,13 +284,16 @@ pub fn update_basis_from_hdf5chk(scf_data: &mut SCF) { println!("taking basis set information from chkfile"); let (loaded_nbasis, loaded_nmo, loaded_spin_channel) = chkfile::load_shapes(&chkfile).unwrap(); - let (cint_raw_data, ecp_raw, basis4elem) = chkfile::load_cint_data(&chkfile); + let (cint_raw_data, ecp_raw, basis4elem, cint_type) = chkfile::load_cint_data(&chkfile); if let Some(cint_raw_data) = cint_raw_data { // mol.cint_atm = cint_raw_data.0; // mol.cint_bas = cint_raw_data.1; // mol.cint_env = cint_raw_data.2; // mol.cint_ecpbas = ecp_raw; scf_data.mol.update_from_cint(loaded_nbasis, basis4elem, cint_raw_data, ecp_raw); + if let Some(cint_type) = cint_type { + scf_data.mol.cint_type = cint_type; + } } else { panic!("Failed to load the basis set information from chkfile"); } @@ -323,7 +328,7 @@ pub fn initial_guess_from_hdf5chk(mol: &mut Molecule, scftype: &SCFType, chkfile ); } else { println!("The basis set in chkfile does not match the input basis set (loaded: {} basis, {} MOs; target: {} basis, {} MOs), \ntrying to project...", loaded_nbasis, loaded_nmo, mol.num_basis, mol.num_state); - let (cint_raw_data, ecp_raw, basis4elem) = chkfile::load_cint_data(chkfile); + let (cint_raw_data, ecp_raw, basis4elem, cint_type) = chkfile::load_cint_data(chkfile); let (mo2, _, occ2) = initial_guess_from_raw( loaded_eigenvectors, loaded_eigenvalues, @@ -333,7 +338,12 @@ pub fn initial_guess_from_hdf5chk(mol: &mut Molecule, scftype: &SCFType, chkfile loaded_nbasis, mol.ctrl.print_level ); - let mo = proj::proj_mo(mol, mo2, cint_raw_data.unwrap(), ecp_raw); + let mut mol_source = Molecule::init_mol(); + mol_source.cint_type = cint_type.unwrap_or(mol.cint_type); + mol_source.ctrl.print_level = mol.ctrl.print_level; + mol_source.update_from_cint(0, None, cint_raw_data.unwrap(), ecp_raw); + assert!(proj::check_proj_sanity(&mol, &mol_source)); + let mo = proj::proj_mo(mol, &mol_source, mo2); return (mo, [vec![], vec![]], occ2); } diff --git a/src/initial_guess/proj.rs b/src/initial_guess/proj.rs index e846a35..b8efeef 100644 --- a/src/initial_guess/proj.rs +++ b/src/initial_guess/proj.rs @@ -50,29 +50,21 @@ fn normalize_mo(mo: &mut MatrixFull, s: &MatrixFull, device: &DeviceBL *mo = MatrixFull::from_vec(shape, mo_normalized.into_shape(-1).into_vec()).unwrap(); } - - -/// Project MO coefficients from a source basis (given as raw cint data) to a target molecule's basis. +/// Project MO coefficients from a source molecule's basis to a target molecule's basis. /// /// # Arguments /// - `mol_target`: the target molecule (basis set 2) -/// - `mo_source`: MO coefficients in the source basis (C1 for mol_source) -/// - `cint_source`: raw cint data for the source basis (basis set 1) -/// - `ecp_source`: optional ECP data for the source basis +/// - `mol_source`: the source molecule (basis set 1) +/// - `mo_source`: MO coefficients in the source basis /// /// # Returns /// - MO coefficients projected to the target basis (C2), normalized so that C2^T S22 C2 = I -pub fn proj_mo(mol_target: &Molecule, mo_source: [MatrixFull;2], cint_source: (Vec>, Vec>, Vec), ecp_source: Option>>) -> [MatrixFull;2] { - let mut mol_source = Molecule::init_mol(); - mol_source.cint_type = mol_target.cint_type.clone(); - mol_source.ctrl.print_level = mol_target.ctrl.print_level; - mol_source.update_from_cint(0, None, cint_source, ecp_source); - +pub fn proj_mo(mol_target: &Molecule, mol_source: &Molecule, mo_source: [MatrixFull;2]) -> [MatrixFull;2] { // S22 = target self-overlap let s22_full = mol_target.int_ij_matrixupper("ovlp".to_string()).to_matrixfull().unwrap(); // S21 = - let s21 = mol_target.int_cross(&mol_source, "ovlp".to_string()); + let s21 = mol_target.int_cross(mol_source, "ovlp".to_string()); let device = DeviceBLAS::default(); let s21_tsr = s21.to_rstsr(&device); @@ -96,3 +88,20 @@ pub fn proj_mo(mol_target: &Molecule, mo_source: [MatrixFull;2], cint_sourc mo_target } + +pub fn check_proj_sanity(mol_target: &Molecule, mol_source: &Molecule) -> bool { + let mut san = true; + if mol_target.cint_type != mol_source.cint_type { + println!("CintType mismatch between target and source molecule (target: {:?}, source: {:?})", mol_target.cint_type, mol_source.cint_type); + san = false; + } + if mol_target.start_mo != mol_source.start_mo { + println!("start_mo mismatch between target and source molecule (target: {}, source: {})", mol_target.start_mo, mol_source.start_mo); + san = false; + } + if mol_target.has_ecp() != mol_source.has_ecp() { + println!("ECP presence mismatch between target and source molecule (target has ECP: {}, source has ECP: {})", mol_target.has_ecp(), mol_source.has_ecp()); + san = false; + } + san +} \ No newline at end of file diff --git a/src/molecule_io/mod.rs b/src/molecule_io/mod.rs index 537e484..64c6982 100644 --- a/src/molecule_io/mod.rs +++ b/src/molecule_io/mod.rs @@ -321,6 +321,10 @@ impl Molecule { } + pub fn has_ecp(&self) -> bool { + self.cint_ecpbas.is_some() + } + pub fn update_from_cint(&mut self, num_basis: usize, basis4elem: Option>, cint_raw: (Vec>, Vec>, Vec), ecp_raw: Option>>) { self.num_basis = num_basis; if let Some(basis4elem_value) = basis4elem { -- Gitee