From 26af4d0209fb150784ccc008d4c848f39c5538cf Mon Sep 17 00:00:00 2001 From: v-boqinzhang Date: Tue, 24 Feb 2026 16:47:32 -0800 Subject: [PATCH 01/29] preparation on gdm.cpp --- .../microsoft/scf/src/scf/scf_impl.cpp | 3 +- .../microsoft/scf/src/scf_algorithm/gdm.cpp | 304 ++++++++++-------- 2 files changed, 169 insertions(+), 138 deletions(-) diff --git a/cpp/src/qdk/chemistry/algorithms/microsoft/scf/src/scf/scf_impl.cpp b/cpp/src/qdk/chemistry/algorithms/microsoft/scf/src/scf/scf_impl.cpp index 7ca99afdc..b5413cfa0 100644 --- a/cpp/src/qdk/chemistry/algorithms/microsoft/scf/src/scf/scf_impl.cpp +++ b/cpp/src/qdk/chemistry/algorithms/microsoft/scf/src/scf/scf_impl.cpp @@ -1060,7 +1060,8 @@ SCFImpl::evaluate_trial_density_energy_and_fock( #endif if (ctx_.cfg->mpi.world_rank == 0) { - if (ctx_.cfg->scf_orbital_type == SCFOrbitalType::Unrestricted) { + if (ctx_.cfg->scf_orbital_type == SCFOrbitalType::Unrestricted || + ctx_.cfg->scf_orbital_type == SCFOrbitalType::RestrictedOpenShell) { F_matrix += (J_matrix.block(0, 0, num_atomic_orbitals_, num_atomic_orbitals_) + J_matrix.block(num_atomic_orbitals_, 0, num_atomic_orbitals_, diff --git a/cpp/src/qdk/chemistry/algorithms/microsoft/scf/src/scf_algorithm/gdm.cpp b/cpp/src/qdk/chemistry/algorithms/microsoft/scf/src/scf_algorithm/gdm.cpp index 74a6ac52b..09aea04ce 100644 --- a/cpp/src/qdk/chemistry/algorithms/microsoft/scf/src/scf_algorithm/gdm.cpp +++ b/cpp/src/qdk/chemistry/algorithms/microsoft/scf/src/scf_algorithm/gdm.cpp @@ -74,6 +74,70 @@ static void apply_orbital_rotation(RowMajorMatrix& C, const int spin_index, num_molecular_orbitals); } +/** + * @brief Compute restricted/unrestricted orbital gradients for all spins + * @param[in] F Fock matrix in AO basis + * @param[in] C Molecular orbital coefficient matrix + * @param[in] num_electrons Occupied orbital counts per spin component + * @param[in] rotation_offset Starting index for each spin's rotation slice + * @param[in] rotation_size Number of rotation parameters per spin + * (n_occ*n_virt) + * @param[in] num_orbital_spin_blocks Number of spin blocks to iterate + * @param[in] num_molecular_orbitals Total molecular orbitals in the system + * @param[in] scf_orbital_type Spin symmetry used across SCF algorithms + * @param[out] gradient Output gradient vector (concatenated across spins) + */ +static void compute_restricted_unrestricted_gradient( + const RowMajorMatrix& F, const RowMajorMatrix& C, + const std::vector& num_electrons, + const std::vector& rotation_offset, + const std::vector& rotation_size, int num_orbital_spin_blocks, + int num_molecular_orbitals, SCFOrbitalType scf_orbital_type, + Eigen::VectorXd& gradient) { + int total_rotation_size = 0; + for (int i = 0; i < num_orbital_spin_blocks; ++i) { + total_rotation_size += rotation_size[i]; + } + gradient.setZero(total_rotation_size); + + for (int spin_index = 0; spin_index < num_orbital_spin_blocks; ++spin_index) { + const int num_occupied_orbitals = num_electrons[spin_index]; + const int num_virtual_orbitals = + num_molecular_orbitals - num_occupied_orbitals; + const int spin_rotation_size = rotation_size[spin_index]; + + if (spin_rotation_size == 0) { + continue; + } + + RowMajorMatrix F_MO = + C.block(num_molecular_orbitals * spin_index, 0, num_molecular_orbitals, + num_molecular_orbitals) + .transpose() * + F.block(num_molecular_orbitals * spin_index, 0, num_molecular_orbitals, + num_molecular_orbitals) * + C.block(num_molecular_orbitals * spin_index, 0, num_molecular_orbitals, + num_molecular_orbitals); + + // Extract occupied-virtual block and compute gradient + // The -4.0 before F_{ia} comes from derivative of energy w.r.t. kappa + // Reference: Helgaker, T., Jørgensen, P., & Olsen, J. (2000). Molecular + // electronic-structure theory, Eq. 10.8.34 (2013 reprint edition) + // -4.0 is for restricted closed-shell system. For unrestricted systems, the + // gradient is computed separately for each spin component, in that case the + // coefficient before F_{ia, spin} is -2.0 + RowMajorMatrix gradient_matrix = + -((scf_orbital_type == SCFOrbitalType::Unrestricted) ? 2.0 : 4.0) * + F_MO.block(0, num_occupied_orbitals, num_occupied_orbitals, + num_virtual_orbitals); + + gradient.segment(rotation_offset[spin_index], spin_rotation_size) = + Eigen::Map(gradient_matrix.data(), + spin_rotation_size) + .eval(); + } +} + /** * @brief Functor for evaluating GDM line search objective */ @@ -90,22 +154,25 @@ class GDMLineFunctor { * @param rotation_offset Starting index for each spin's rotation slice * @param rotation_size Number of rotation parameters per spin (n_occ*n_virt) * @param num_molecular_orbitals Total molecular orbitals in the system - * @param unrestricted Whether alpha/beta densities are treated separately + * @param scf_orbital_type Spin symmetry used across SCF algorithms */ GDMLineFunctor(const SCFImpl& scf_impl, const RowMajorMatrix& C_pseudo_canonical, const std::vector& num_electrons, const std::vector& rotation_offset, const std::vector& rotation_size, - int num_molecular_orbitals, bool unrestricted) + int num_molecular_orbitals, SCFOrbitalType scf_orbital_type) : scf_impl_(scf_impl), C_pseudo_canonical_(C_pseudo_canonical), num_electrons_(num_electrons), rotation_offset_(rotation_offset), rotation_size_(rotation_size), - num_density_matrices_(unrestricted ? 2 : 1), + num_orbital_spin_blocks_( + scf_orbital_type == SCFOrbitalType::Unrestricted ? 2 : 1), + num_density_matrices_( + scf_orbital_type == SCFOrbitalType::Restricted ? 1 : 2), num_molecular_orbitals_(num_molecular_orbitals), - unrestricted_(unrestricted), + scf_orbital_type_(scf_orbital_type), cached_kappa_(Eigen::VectorXd()), cached_energy_(std::numeric_limits::infinity()) {} @@ -157,9 +224,10 @@ class GDMLineFunctor { const std::vector& rotation_size_; // Value parameters + const int num_orbital_spin_blocks_; const int num_density_matrices_; const int num_molecular_orbitals_; - const bool unrestricted_; + const SCFOrbitalType scf_orbital_type_; // Cache for avoiding redundant Fock matrix computation Eigen::VectorXd cached_kappa_; // Cached kappa vector @@ -181,7 +249,7 @@ double GDMLineFunctor::eval(const Eigen::VectorXd& x) { cached_C_ = C_pseudo_canonical_; // Apply rotation for all spins with kappa_trial - for (int i = 0; i < num_density_matrices_; i++) { + for (int i = 0; i < num_orbital_spin_blocks_; i++) { auto kappa_spin = kappa_trial.segment(rotation_offset_[i], rotation_size_[i]); apply_orbital_rotation(cached_C_, i, kappa_spin, num_electrons_[i], @@ -194,7 +262,8 @@ double GDMLineFunctor::eval(const Eigen::VectorXd& x) { for (int i = 0; i < num_density_matrices_; i++) { const int num_occupied_orbitals = num_electrons_[i]; - const double occupation_factor = unrestricted_ ? 1.0 : 2.0; + const double occupation_factor = + (scf_orbital_type_ == SCFOrbitalType::Unrestricted) ? 1.0 : 2.0; cached_P_.block(num_molecular_orbitals_ * i, 0, num_molecular_orbitals_, num_molecular_orbitals_) = @@ -226,47 +295,11 @@ Eigen::VectorXd GDMLineFunctor::grad(const Eigen::VectorXd& x) { eval(x); } - // Initialize the full gradient vector (concatenated for all spins) - int total_rotation_size = 0; - for (int i = 0; i < num_density_matrices_; i++) { - total_rotation_size += rotation_size_[i]; - } - Eigen::VectorXd gradient = Eigen::VectorXd::Zero(total_rotation_size); - - // Compute gradient for each spin component - for (int i = 0; i < num_density_matrices_; i++) { - const int num_occupied_orbitals = num_electrons_[i]; - const int num_virtual_orbitals = - num_molecular_orbitals_ - num_occupied_orbitals; - - // Transform Fock matrix to MO basis: F_MO = C^T * F * C - RowMajorMatrix F_MO = - cached_C_ - .block(num_molecular_orbitals_ * i, 0, num_molecular_orbitals_, - num_molecular_orbitals_) - .transpose() * - cached_F_.block(num_molecular_orbitals_ * i, 0, num_molecular_orbitals_, - num_molecular_orbitals_) * - cached_C_.block(num_molecular_orbitals_ * i, 0, num_molecular_orbitals_, - num_molecular_orbitals_); - - // Extract occupied-virtual block and compute gradient - // The -4.0 before F_{ia} comes from derivative of energy w.r.t. kappa - // Reference: Helgaker, T., Jørgensen, P., & Olsen, J. (2000). Molecular - // electronic-structure theory, Eq. 10.8.34 (2013 reprint edition) - // -4.0 is for restricted closed-shell system. For unrestricted systems, the - // gradient is computed separately for each spin component, in that case the - // coefficient before F_{ia, spin} is -2.0 - RowMajorMatrix gradient_matrix = - -(unrestricted_ ? 2.0 : 4.0) * F_MO.block(0, num_occupied_orbitals, - num_occupied_orbitals, - num_virtual_orbitals); - - // Flatten matrix to vector and store in appropriate segment - gradient.segment(rotation_offset_[i], rotation_size_[i]) = - Eigen::Map(gradient_matrix.data(), - rotation_size_[i]); - } + Eigen::VectorXd gradient; + compute_restricted_unrestricted_gradient( + cached_F_, cached_C_, num_electrons_, rotation_offset_, rotation_size_, + num_orbital_spin_blocks_, num_molecular_orbitals_, scf_orbital_type_, + gradient); return gradient; } @@ -315,13 +348,16 @@ class GDM { * pseudo-canonical orbital basis, K_new = Uoo^T * K_old * Uvv * @param[in,out] history History matrix block to be transformed (either * history_dgrad or history_kappa) + * @param[in] u_left Left rotation matrix (e.g., Uoo or Uaa) + * @param[in] u_right Right rotation matrix (e.g., Uvv or Uaa) * @param[in] history_size Number of history entries * @param[in] num_occupied_orbitals Number of electrons for current spin * @param[in] num_molecular_orbitals Number of molecular orbitals * */ void transform_history_(Eigen::Block& history, - const int history_size, + const RowMajorMatrix& u_left, + const RowMajorMatrix& u_right, const int history_size, const int num_occupied_orbitals, const int num_molecular_orbitals); @@ -344,6 +380,20 @@ class GDM { Eigen::Block history_dgrad_spin, Eigen::VectorBlock current_gradient_spin); + /** + * @brief Build pseudo-canonical orbitals and initial Hessian across spins + * @param[in] F Fock matrix in AO basis + * @param[in,out] C Molecular orbital coefficient matrix + * @param[in] num_density_matrices Number of density matrices + * @param[in] num_molecular_orbitals Total molecular orbitals in the system + * @param[in] scf_orbital_type Spin symmetry used across SCF algorithms + * @param[out] initial_hessian Output concatenated initial Hessian + */ + void build_restricted_unrestricted_pseudo_canonical_orbitals_hessian( + const RowMajorMatrix& F, RowMajorMatrix& C, int num_density_matrices, + int num_molecular_orbitals, SCFOrbitalType scf_orbital_type, + Eigen::VectorXd& initial_hessian); + /// Reference to SCFContext const SCFContext& ctx_; ///< Reference to SCFContext /// Energy change from the last step @@ -467,6 +517,8 @@ GDM::GDM(const SCFContext& ctx, int history_size_limit) } void GDM::transform_history_(Eigen::Block& history, + const RowMajorMatrix& u_left, + const RowMajorMatrix& u_right, const int history_size, const int num_occupied_orbitals, const int num_molecular_orbitals) { @@ -490,15 +542,15 @@ void GDM::transform_history_(Eigen::Block& history, RowMajorMatrix::Zero(num_occupied_orbitals, num_virtual_orbitals); for (int line = 0; line < history_size; line++) { double* history_line_ptr = history.row(line).data(); - // K_ov (new) = Uoo^T * K_ov * Uvv + // K_ov (new) = U_left(oo)^T * K_ov * U_right(vv) blas::gemm(blas::Layout::RowMajor, blas::Op::NoTrans, blas::Op::NoTrans, num_occupied_orbitals, num_virtual_orbitals, num_virtual_orbitals, 1.0, history_line_ptr, - num_virtual_orbitals, Uvv_.data(), num_virtual_orbitals, 0.0, + num_virtual_orbitals, u_right.data(), num_virtual_orbitals, 0.0, temp.data(), num_virtual_orbitals); blas::gemm(blas::Layout::RowMajor, blas::Op::Trans, blas::Op::NoTrans, num_occupied_orbitals, num_virtual_orbitals, - num_occupied_orbitals, 1.0, Uoo_.data(), num_occupied_orbitals, + num_occupied_orbitals, 1.0, u_left.data(), num_occupied_orbitals, temp.data(), num_virtual_orbitals, 0.0, history_line_ptr, num_virtual_orbitals); } @@ -579,10 +631,10 @@ void GDM::generate_pseudo_canonical_orbital_( // Transform the vectors in history_kappa and history_dgrad to // accommodate current pseudo-canonical orbitals - transform_history_(history_kappa_spin, history_size_, num_occupied_orbitals, - num_molecular_orbitals); - transform_history_(history_dgrad_spin, history_size_, num_occupied_orbitals, - num_molecular_orbitals); + transform_history_(history_kappa_spin, Uoo_, Uvv_, history_size_, + num_occupied_orbitals, num_molecular_orbitals); + transform_history_(history_dgrad_spin, Uoo_, Uvv_, history_size_, + num_occupied_orbitals, num_molecular_orbitals); // Transform the gradient to accommodate current pseudo-canonical orbitals RowMajorMatrix current_gradient_matrix = @@ -595,6 +647,51 @@ void GDM::generate_pseudo_canonical_orbital_( current_gradient_spin = current_gradient_transformed; } +void GDM::build_restricted_unrestricted_pseudo_canonical_orbitals_hessian( + const RowMajorMatrix& F, RowMajorMatrix& C, int num_density_matrices, + int num_molecular_orbitals, SCFOrbitalType scf_orbital_type, + Eigen::VectorXd& initial_hessian) { + initial_hessian.setZero(total_rotation_size_); + + for (int i = 0; i < num_density_matrices; ++i) { + const int num_occupied_orbitals = num_electrons_[i]; + const int num_virtual_orbitals = + num_molecular_orbitals - num_occupied_orbitals; + + auto history_kappa_spin = history_kappa_.block( + 0, rotation_offset_[i], history_size_limit_, rotation_size_[i]); + auto history_dgrad_spin = history_dgrad_.block( + 0, rotation_offset_[i], history_size_limit_, rotation_size_[i]); + auto current_gradient_spin = + current_gradient_.segment(rotation_offset_[i], rotation_size_[i]); + + // Generate pseudo-canonical orbitals and transform gradient and history + generate_pseudo_canonical_orbital_( + F, C, i, history_kappa_spin, history_dgrad_spin, current_gradient_spin); + + // Build this spin's segment of initial Hessian + // Reference: Helgaker, T., Jorgensen, P., & Olsen, J. (2000). Molecular + // electronic-structure theory, Eq. 10.8.56 (2013 reprint edition) + // 4.0 is for restricted closed-shell system. For unrestricted systems, the + // gradient is computed separately for each spin component, in that case the + // coefficient should be 2.0 + double initial_hessian_coeff = + (scf_orbital_type == SCFOrbitalType::Unrestricted) ? 2.0 : 4.0; + for (int j = 0; j < num_occupied_orbitals; j++) { + for (int v = 0; v < num_virtual_orbitals; v++) { + int index = rotation_offset_[i] + j * num_virtual_orbitals + v; + double pseudo_canonical_energy_diff = + std::abs(pseudo_canonical_eigenvalues_(num_occupied_orbitals + v) - + pseudo_canonical_eigenvalues_(j)); + initial_hessian(index) = + std::max(initial_hessian_coeff * (std::abs(delta_energy_) + + pseudo_canonical_energy_diff), + nonpositive_threshold_); + } + } + } +} + void GDM::iterate(SCFImpl& scf_impl) { QDK_LOG_TRACE_ENTERING(); auto& C = scf_impl.orbitals_matrix(); @@ -615,44 +712,14 @@ void GDM::iterate(SCFImpl& scf_impl) { return; } - // Compute current gradient and dgrad for each spin - for (int i = 0; i < num_density_matrices; ++i) { - const int num_occupied_orbitals = num_electrons_[i]; - const int num_virtual_orbitals = - num_molecular_orbitals - num_occupied_orbitals; - const int rotation_size = num_occupied_orbitals * num_virtual_orbitals; - RowMajorMatrix F_MO = - C.block(num_molecular_orbitals * i, 0, num_molecular_orbitals, - num_molecular_orbitals) - .transpose() * - F.block(num_molecular_orbitals * i, 0, num_molecular_orbitals, - num_molecular_orbitals) * - C.block(num_molecular_orbitals * i, 0, num_molecular_orbitals, - num_molecular_orbitals); + compute_restricted_unrestricted_gradient( + F, C, num_electrons_, rotation_offset_, rotation_size_, + num_density_matrices, num_molecular_orbitals, cfg->scf_orbital_type, + current_gradient_); - // Extract occupied-virtual block and compute gradient - // The -4.0 before F_{ia} comes from derivative of energy w.r.t. kappa - // Reference: Helgaker, T., Jørgensen, P., & Olsen, J. (2000). Molecular - // electronic-structure theory, Eq. 10.8.34 (2013 reprint edition) - // -4.0 is for restricted closed-shell system. For unrestricted systems, the - // gradient is computed separately for each spin component, in that case the - // coefficient before F_{ia, spin} is -2.0 - RowMajorMatrix current_gradient_matrix = - -((cfg->scf_orbital_type == SCFOrbitalType::Unrestricted) ? 2.0 : 4.0) * - F_MO.block(0, num_occupied_orbitals, num_occupied_orbitals, - num_virtual_orbitals); - current_gradient_.segment(rotation_offset_[i], rotation_size_[i]) = - Eigen::Map(current_gradient_matrix.data(), - rotation_size); - - if (gdm_step_count_ != 0) { - // Add new gradient difference to history for this spin - history_dgrad_ - .block(0, rotation_offset_[i], history_size_limit_, rotation_size_[i]) - .row(history_size_) = - current_gradient_.segment(rotation_offset_[i], rotation_size_[i]) - - previous_gradient_.segment(rotation_offset_[i], rotation_size_[i]); - } + if (gdm_step_count_ != 0) { + // Add new gradient difference to history for all spins + history_dgrad_.row(history_size_) = current_gradient_ - previous_gradient_; } // Update history size and manage history overflow. History for both spins are @@ -674,46 +741,10 @@ void GDM::iterate(SCFImpl& scf_impl) { } } - // Build concatenated initial Hessian for all spins - Eigen::VectorXd initial_hessian = Eigen::VectorXd::Zero(total_rotation_size_); - - for (int i = 0; i < num_density_matrices; ++i) { - const int num_occupied_orbitals = num_electrons_[i]; - const int num_virtual_orbitals = - num_molecular_orbitals - num_occupied_orbitals; - - auto history_kappa_spin = history_kappa_.block( - 0, rotation_offset_[i], history_size_limit_, rotation_size_[i]); - auto history_dgrad_spin = history_dgrad_.block( - 0, rotation_offset_[i], history_size_limit_, rotation_size_[i]); - auto current_gradient_spin = - current_gradient_.segment(rotation_offset_[i], rotation_size_[i]); - - // Generate pseudo-canonical orbitals and transform gradient and history - generate_pseudo_canonical_orbital_( - F, C, i, history_kappa_spin, history_dgrad_spin, current_gradient_spin); - - // Build this spin's segment of initial Hessian - // Reference: Helgaker, T., Jørgensen, P., & Olsen, J. (2000). Molecular - // electronic-structure theory, Eq. 10.8.56 (2013 reprint edition) - // 4.0 is for restricted closed-shell system. For unrestricted systems, the - // gradient is computed separately for each spin component, in that case the - // coefficient should be 2.0 - double initial_hessian_coeff = - (cfg->scf_orbital_type == SCFOrbitalType::Unrestricted) ? 2.0 : 4.0; - for (int j = 0; j < num_occupied_orbitals; j++) { - for (int v = 0; v < num_virtual_orbitals; v++) { - int index = rotation_offset_[i] + j * num_virtual_orbitals + v; - double pseudo_canonical_energy_diff = - std::abs(pseudo_canonical_eigenvalues_(num_occupied_orbitals + v) - - pseudo_canonical_eigenvalues_(j)); - initial_hessian(index) = - std::max(initial_hessian_coeff * (std::abs(delta_energy_) + - pseudo_canonical_energy_diff), - nonpositive_threshold_); - } - } - } + Eigen::VectorXd initial_hessian; + build_restricted_unrestricted_pseudo_canonical_orbitals_hessian( + F, C, num_density_matrices, num_molecular_orbitals, cfg->scf_orbital_type, + initial_hessian); double latest_inverse_rho = 1.0; // BFGS two-loop recursion on concatenated vectors (runs once for all spins) @@ -828,10 +859,9 @@ void GDM::iterate(SCFImpl& scf_impl) { RowMajorMatrix C_pseudo_canonical = C.eval(); // Create line search functor for energy evaluation - GDMLineFunctor line_functor( - scf_impl, C_pseudo_canonical, num_electrons_, rotation_offset_, - rotation_size_, num_molecular_orbitals, - cfg->scf_orbital_type == SCFOrbitalType::Unrestricted); + GDMLineFunctor line_functor(scf_impl, C_pseudo_canonical, num_electrons_, + rotation_offset_, rotation_size_, + num_molecular_orbitals, cfg->scf_orbital_type); Eigen::VectorXd start_kappa = Eigen::VectorXd::Zero(kappa_.size()); Eigen::VectorXd kappa_dir = kappa_; // Search direction From a27effe542835593b15804dfadc6dafc2cfc7d22 Mon Sep 17 00:00:00 2001 From: v-boqinzhang Date: Wed, 25 Feb 2026 17:22:34 -0800 Subject: [PATCH 02/29] compute_restricted_open_shell_gradient --- .../microsoft/scf/src/scf/scf_impl.cpp | 7 + .../microsoft/scf/src/scf/scf_impl.h | 15 ++ .../microsoft/scf/src/scf_algorithm/gdm.cpp | 236 ++++++++++++++---- 3 files changed, 203 insertions(+), 55 deletions(-) diff --git a/cpp/src/qdk/chemistry/algorithms/microsoft/scf/src/scf/scf_impl.cpp b/cpp/src/qdk/chemistry/algorithms/microsoft/scf/src/scf/scf_impl.cpp index b5413cfa0..089b2d7e9 100644 --- a/cpp/src/qdk/chemistry/algorithms/microsoft/scf/src/scf/scf_impl.cpp +++ b/cpp/src/qdk/chemistry/algorithms/microsoft/scf/src/scf/scf_impl.cpp @@ -1097,4 +1097,11 @@ SCFImpl::evaluate_trial_density_energy_and_fock( return {total_energy, F_matrix}; } +void SCFImpl::build_jk_matrices(const RowMajorMatrix& density_matrix, + RowMajorMatrix& J, RowMajorMatrix& K) const { + QDK_LOG_TRACE_ENTERING(); + auto [alpha, beta, omega] = get_hyb_coeff_(); + eri_->build_JK(density_matrix.data(), J.data(), K.data(), alpha, beta, omega); +} + } // namespace qdk::chemistry::scf diff --git a/cpp/src/qdk/chemistry/algorithms/microsoft/scf/src/scf/scf_impl.h b/cpp/src/qdk/chemistry/algorithms/microsoft/scf/src/scf/scf_impl.h index 67f5ade5d..e0c27a66e 100644 --- a/cpp/src/qdk/chemistry/algorithms/microsoft/scf/src/scf/scf_impl.h +++ b/cpp/src/qdk/chemistry/algorithms/microsoft/scf/src/scf/scf_impl.h @@ -148,6 +148,21 @@ class SCFImpl { */ const RowMajorMatrix& get_fock_matrix() const { return F_; } + /** + * @brief Get the core Hamiltonian matrix + * @return Reference to core Hamiltonian H + */ + const RowMajorMatrix& get_core_hamiltonian() const { return H_; } + + /** + * @brief Build Coulomb (J) and exchange (K) matrices for a density matrix + * @param density_matrix Density matrix (num_density_matrices x NAO x NAO) + * @param J Output Coulomb matrix (same size as density_matrix) + * @param K Output exchange matrix (same size as density_matrix) + */ + void build_jk_matrices(const RowMajorMatrix& density_matrix, + RowMajorMatrix& J, RowMajorMatrix& K) const; + /** * @brief Get the molecular orbital coefficient matrix * @see SCF::get_orbitals_matrix() for API details diff --git a/cpp/src/qdk/chemistry/algorithms/microsoft/scf/src/scf_algorithm/gdm.cpp b/cpp/src/qdk/chemistry/algorithms/microsoft/scf/src/scf_algorithm/gdm.cpp index 09aea04ce..9e69c182a 100644 --- a/cpp/src/qdk/chemistry/algorithms/microsoft/scf/src/scf_algorithm/gdm.cpp +++ b/cpp/src/qdk/chemistry/algorithms/microsoft/scf/src/scf_algorithm/gdm.cpp @@ -84,7 +84,6 @@ static void apply_orbital_rotation(RowMajorMatrix& C, const int spin_index, * (n_occ*n_virt) * @param[in] num_orbital_spin_blocks Number of spin blocks to iterate * @param[in] num_molecular_orbitals Total molecular orbitals in the system - * @param[in] scf_orbital_type Spin symmetry used across SCF algorithms * @param[out] gradient Output gradient vector (concatenated across spins) */ static void compute_restricted_unrestricted_gradient( @@ -92,8 +91,7 @@ static void compute_restricted_unrestricted_gradient( const std::vector& num_electrons, const std::vector& rotation_offset, const std::vector& rotation_size, int num_orbital_spin_blocks, - int num_molecular_orbitals, SCFOrbitalType scf_orbital_type, - Eigen::VectorXd& gradient) { + int num_molecular_orbitals, Eigen::VectorXd& gradient) { int total_rotation_size = 0; for (int i = 0; i < num_orbital_spin_blocks; ++i) { total_rotation_size += rotation_size[i]; @@ -127,7 +125,7 @@ static void compute_restricted_unrestricted_gradient( // gradient is computed separately for each spin component, in that case the // coefficient before F_{ia, spin} is -2.0 RowMajorMatrix gradient_matrix = - -((scf_orbital_type == SCFOrbitalType::Unrestricted) ? 2.0 : 4.0) * + -((num_orbital_spin_blocks == 2) ? 2.0 : 4.0) * F_MO.block(0, num_occupied_orbitals, num_occupied_orbitals, num_virtual_orbitals); @@ -138,6 +136,94 @@ static void compute_restricted_unrestricted_gradient( } } +static void compute_restricted_open_shell_gradient( + const SCFImpl& scf_impl, const RowMajorMatrix& C, + const RowMajorMatrix& density_matrix, const std::vector& num_electrons, + const std::vector& rotation_size, int num_molecular_orbitals, + RowMajorMatrix& generalized_fock_mo, Eigen::VectorXd& gradient) { + const int total_rotation_size = rotation_size[0]; + gradient.setZero(total_rotation_size); + + const int num_closed_orbitals = num_electrons[1]; + const int num_open_orbitals = num_electrons[0] - num_closed_orbitals; + const int num_virtual_orbitals = num_molecular_orbitals - num_electrons[0]; + + if (generalized_fock_mo.rows() != num_molecular_orbitals || + generalized_fock_mo.cols() != num_molecular_orbitals) { + throw std::invalid_argument( + "generalized_fock_mo must be preallocated to MO dimensions."); + } + + const int num_atomic_orbitals = scf_impl.get_num_atomic_orbitals(); + const auto& H_ao_full = scf_impl.get_core_hamiltonian(); + const RowMajorMatrix H_ao = + H_ao_full.block(0, 0, num_atomic_orbitals, num_atomic_orbitals); + + RowMajorMatrix J_ao = + RowMajorMatrix::Zero(2 * num_atomic_orbitals, num_atomic_orbitals); + RowMajorMatrix K_ao = + RowMajorMatrix::Zero(2 * num_atomic_orbitals, num_atomic_orbitals); + scf_impl.build_jk_matrices(density_matrix, J_ao, K_ao); + + const auto J_alpha_ao = Eigen::Map( + J_ao.data(), num_atomic_orbitals, num_atomic_orbitals); + const auto J_beta_ao = Eigen::Map( + J_ao.data() + num_atomic_orbitals * num_atomic_orbitals, + num_atomic_orbitals, num_atomic_orbitals); + const auto K_alpha_ao = Eigen::Map( + K_ao.data(), num_atomic_orbitals, num_atomic_orbitals); + const auto K_beta_ao = Eigen::Map( + K_ao.data() + num_atomic_orbitals * num_atomic_orbitals, + num_atomic_orbitals, num_atomic_orbitals); + + const auto C_ao_mo = + C.block(0, 0, num_atomic_orbitals, num_molecular_orbitals); + + RowMajorMatrix H_mo = C_ao_mo.transpose() * H_ao * C_ao_mo; + RowMajorMatrix J_alpha_mo = C_ao_mo.transpose() * J_alpha_ao * C_ao_mo; + RowMajorMatrix J_beta_mo = C_ao_mo.transpose() * J_beta_ao * C_ao_mo; + RowMajorMatrix K_alpha_mo = C_ao_mo.transpose() * K_alpha_ao * C_ao_mo; + RowMajorMatrix K_beta_mo = C_ao_mo.transpose() * K_beta_ao * C_ao_mo; + + RowMajorMatrix F_I = H_mo + 2.0 * J_beta_mo - K_beta_mo; + RowMajorMatrix F_A = J_alpha_mo - J_beta_mo - 0.5 * (K_alpha_mo - K_beta_mo); + RowMajorMatrix Q = J_alpha_mo - J_beta_mo - (K_alpha_mo - K_beta_mo); + + RowMajorMatrix F_sum = F_I + F_A; + generalized_fock_mo.block(0, 0, num_closed_orbitals, num_molecular_orbitals) = + 2.0 * F_sum.block(0, 0, num_closed_orbitals, num_molecular_orbitals); + generalized_fock_mo.block(num_closed_orbitals, 0, num_open_orbitals, + num_molecular_orbitals) = + (F_I + Q).block(num_closed_orbitals, 0, num_open_orbitals, + num_molecular_orbitals); + + int offset = 0; + RowMajorMatrix grad_iw = + 2.0 * (generalized_fock_mo.block(0, num_closed_orbitals, + num_closed_orbitals, num_open_orbitals) - + generalized_fock_mo + .block(num_closed_orbitals, 0, num_open_orbitals, + num_closed_orbitals) + .transpose()); + gradient.segment(offset, num_closed_orbitals * num_open_orbitals) = + Eigen::Map(grad_iw.data(), grad_iw.size()).eval(); + offset += num_closed_orbitals * num_open_orbitals; + + RowMajorMatrix grad_wa = + 2.0 * generalized_fock_mo.block(num_closed_orbitals, + num_closed_orbitals + num_open_orbitals, + num_open_orbitals, num_virtual_orbitals); + gradient.segment(offset, num_open_orbitals * num_virtual_orbitals) = + Eigen::Map(grad_wa.data(), grad_wa.size()).eval(); + offset += num_open_orbitals * num_virtual_orbitals; + + RowMajorMatrix grad_ia = 2.0 * generalized_fock_mo.block( + 0, num_closed_orbitals + num_open_orbitals, + num_closed_orbitals, num_virtual_orbitals); + gradient.segment(offset, num_closed_orbitals * num_virtual_orbitals) = + Eigen::Map(grad_ia.data(), grad_ia.size()).eval(); +} + /** * @brief Functor for evaluating GDM line search objective */ @@ -296,10 +382,17 @@ Eigen::VectorXd GDMLineFunctor::grad(const Eigen::VectorXd& x) { } Eigen::VectorXd gradient; - compute_restricted_unrestricted_gradient( - cached_F_, cached_C_, num_electrons_, rotation_offset_, rotation_size_, - num_orbital_spin_blocks_, num_molecular_orbitals_, scf_orbital_type_, - gradient); + if (scf_orbital_type_ == SCFOrbitalType::RestrictedOpenShell) { + RowMajorMatrix generalized_fock_mo = + RowMajorMatrix::Zero(num_molecular_orbitals_, num_molecular_orbitals_); + compute_restricted_open_shell_gradient( + scf_impl_, cached_C_, cached_P_, num_electrons_, rotation_size_, + num_molecular_orbitals_, generalized_fock_mo, gradient); + } else { + compute_restricted_unrestricted_gradient( + cached_F_, cached_C_, num_electrons_, rotation_offset_, rotation_size_, + num_orbital_spin_blocks_, num_molecular_orbitals_, gradient); + } return gradient; } @@ -384,14 +477,11 @@ class GDM { * @brief Build pseudo-canonical orbitals and initial Hessian across spins * @param[in] F Fock matrix in AO basis * @param[in,out] C Molecular orbital coefficient matrix - * @param[in] num_density_matrices Number of density matrices - * @param[in] num_molecular_orbitals Total molecular orbitals in the system - * @param[in] scf_orbital_type Spin symmetry used across SCF algorithms + * @param[in] num_molecular_orbitals Total number of molecular orbitals * @param[out] initial_hessian Output concatenated initial Hessian */ void build_restricted_unrestricted_pseudo_canonical_orbitals_hessian( - const RowMajorMatrix& F, RowMajorMatrix& C, int num_density_matrices, - int num_molecular_orbitals, SCFOrbitalType scf_orbital_type, + const RowMajorMatrix& F, RowMajorMatrix& C, int num_molecular_orbitals, Eigen::VectorXd& initial_hessian); /// Reference to SCFContext @@ -428,6 +518,9 @@ class GDM { /// Eigenvalues of pseudo-canonical orbitals, used for building Hessian Eigen::VectorXd pseudo_canonical_eigenvalues_; + /// Generalized Fock matrix in MO basis for ROHF + RowMajorMatrix generalized_fock_mo_; + /// Horizontal rotation matrix of occupied orbitals RowMajorMatrix Uoo_; /// Horizontal rotation matrix of virtual orbitals @@ -438,16 +531,20 @@ class GDM { /// Energy of the last accepted step, used to decide if we rescale the kappa /// vector in this step double last_accepted_energy_; - int gdm_step_count_; // GDM iteration counter - int num_density_matrices_; // Number of density matrices (1 for restricted, 2 - // for unrestricted) + int gdm_step_count_; // GDM iteration counter + const int num_orbital_spin_blocks_; + const int num_density_matrices_; }; GDM::GDM(const SCFContext& ctx, int history_size_limit) : ctx_(ctx), history_size_limit_(history_size_limit), last_accepted_energy_(std::numeric_limits::infinity()), - gdm_step_count_(0) { + gdm_step_count_(0), + num_orbital_spin_blocks_( + ctx.cfg->scf_orbital_type == SCFOrbitalType::Unrestricted ? 2 : 1), + num_density_matrices_( + ctx.cfg->scf_orbital_type == SCFOrbitalType::Restricted ? 1 : 2) { QDK_LOG_TRACE_ENTERING(); const auto& cfg = *ctx.cfg; const auto& mol = *ctx.mol; @@ -466,6 +563,12 @@ GDM::GDM(const SCFContext& ctx, int history_size_limit) num_electrons_ = {num_alpha_electrons, num_beta_electrons}; history_size_ = 0; pseudo_canonical_eigenvalues_ = Eigen::VectorXd::Zero(num_molecular_orbitals); + if (cfg.scf_orbital_type == SCFOrbitalType::RestrictedOpenShell) { + generalized_fock_mo_ = + RowMajorMatrix::Zero(num_molecular_orbitals, num_molecular_orbitals); + } else { + generalized_fock_mo_ = RowMajorMatrix(0, 0); + } if (history_size_limit < 1) { throw std::invalid_argument( "GDM history size limit must be at least 1, got: " + @@ -474,36 +577,57 @@ GDM::GDM(const SCFContext& ctx, int history_size_limit) QDK_LOGGER().debug("GDM initialized with history_size_limit = {}", history_size_limit_); - num_density_matrices_ = - (cfg.scf_orbital_type == SCFOrbitalType::Unrestricted) ? 2 : 1; // Calculate rotation sizes for each spin - rotation_size_.resize(num_density_matrices_); - rotation_offset_.resize(num_density_matrices_); - - total_rotation_size_ = 0; - for (int spin_index = 0; spin_index < num_density_matrices_; spin_index++) { - const int num_occupied_orbitals = num_electrons_[spin_index]; - const int num_virtual_orbitals = - num_molecular_orbitals - num_occupied_orbitals; - // Validate dimensions (negative values indicate invalid input) - // Zero occupied or virtual orbitals is valid for unrestricted calculations - // (e.g., H atom has 0 beta electrons) - if (num_occupied_orbitals < 0) { + rotation_size_.resize(num_orbital_spin_blocks_); + rotation_offset_.resize(num_orbital_spin_blocks_); + + if (cfg.scf_orbital_type == SCFOrbitalType::RestrictedOpenShell) { + int num_closed_orbitals = num_electrons_[1]; + int num_open_orbitals = num_electrons_[0] - num_closed_orbitals; + int num_virtual_orbitals = num_molecular_orbitals - num_electrons_[0]; + if (num_closed_orbitals < 0 || num_open_orbitals < 0 || + num_virtual_orbitals < 0) { throw std::invalid_argument( - std::string("GDM: num_occupied_orbitals must be non-negative, got ") + - std::to_string(num_occupied_orbitals) + " for spin " + - std::to_string(spin_index)); + "Invalid ROHF system: " + "num_closed_orbitals=" + + std::to_string(num_closed_orbitals) + + ", num_open_orbitals=" + std::to_string(num_open_orbitals) + + ", num_virtual_orbitals=" + std::to_string(num_virtual_orbitals)); } - if (num_virtual_orbitals < 0) { - throw std::invalid_argument( - std::string("GDM: num_virtual_orbitals must be non-negative, got ") + - std::to_string(num_virtual_orbitals) + " for spin " + - std::to_string(spin_index)); + rotation_size_[0] = num_closed_orbitals * num_open_orbitals + + num_open_orbitals * num_virtual_orbitals + + num_closed_orbitals * num_virtual_orbitals; + rotation_offset_[0] = 0; + total_rotation_size_ = rotation_size_[0]; + } else { + total_rotation_size_ = 0; + for (int spin_index = 0; spin_index < num_orbital_spin_blocks_; + spin_index++) { + const int num_occupied_orbitals = num_electrons_[spin_index]; + const int num_virtual_orbitals = + num_molecular_orbitals - num_occupied_orbitals; + // Validate dimensions (negative values indicate invalid input) + // Zero occupied or virtual orbitals is valid for unrestricted + // calculations (e.g., H atom has 0 beta electrons) + if (num_occupied_orbitals < 0) { + throw std::invalid_argument( + std::string( + "GDM: num_occupied_orbitals must be non-negative, got ") + + std::to_string(num_occupied_orbitals) + " for spin " + + std::to_string(spin_index)); + } + if (num_virtual_orbitals < 0) { + throw std::invalid_argument( + std::string( + "GDM: num_virtual_orbitals must be non-negative, got ") + + std::to_string(num_virtual_orbitals) + " for spin " + + std::to_string(spin_index)); + } + rotation_size_[spin_index] = num_occupied_orbitals * num_virtual_orbitals; + rotation_offset_[spin_index] = total_rotation_size_; + total_rotation_size_ += rotation_size_[spin_index]; } - rotation_size_[spin_index] = num_occupied_orbitals * num_virtual_orbitals; - rotation_offset_[spin_index] = total_rotation_size_; - total_rotation_size_ += rotation_size_[spin_index]; } // Initialize concatenated matrices and vectors @@ -648,12 +772,11 @@ void GDM::generate_pseudo_canonical_orbital_( } void GDM::build_restricted_unrestricted_pseudo_canonical_orbitals_hessian( - const RowMajorMatrix& F, RowMajorMatrix& C, int num_density_matrices, - int num_molecular_orbitals, SCFOrbitalType scf_orbital_type, + const RowMajorMatrix& F, RowMajorMatrix& C, int num_molecular_orbitals, Eigen::VectorXd& initial_hessian) { initial_hessian.setZero(total_rotation_size_); - for (int i = 0; i < num_density_matrices; ++i) { + for (int i = 0; i < num_orbital_spin_blocks_; ++i) { const int num_occupied_orbitals = num_electrons_[i]; const int num_virtual_orbitals = num_molecular_orbitals - num_occupied_orbitals; @@ -675,8 +798,7 @@ void GDM::build_restricted_unrestricted_pseudo_canonical_orbitals_hessian( // 4.0 is for restricted closed-shell system. For unrestricted systems, the // gradient is computed separately for each spin component, in that case the // coefficient should be 2.0 - double initial_hessian_coeff = - (scf_orbital_type == SCFOrbitalType::Unrestricted) ? 2.0 : 4.0; + double initial_hessian_coeff = (num_orbital_spin_blocks_ == 2) ? 2.0 : 4.0; for (int j = 0; j < num_occupied_orbitals; j++) { for (int v = 0; v < num_virtual_orbitals; v++) { int index = rotation_offset_[i] + j * num_virtual_orbitals + v; @@ -700,8 +822,6 @@ void GDM::iterate(SCFImpl& scf_impl) { const auto* cfg = ctx_.cfg; const int num_molecular_orbitals = static_cast(ctx_.num_molecular_orbitals); - const int num_density_matrices = - (cfg->scf_orbital_type == SCFOrbitalType::Unrestricted) ? 2 : 1; // Check if there are any virtual orbitals for any spin component // If not, orbital rotation is not possible and we should skip GDM iteration @@ -712,10 +832,17 @@ void GDM::iterate(SCFImpl& scf_impl) { return; } - compute_restricted_unrestricted_gradient( - F, C, num_electrons_, rotation_offset_, rotation_size_, - num_density_matrices, num_molecular_orbitals, cfg->scf_orbital_type, - current_gradient_); + if (cfg->scf_orbital_type == SCFOrbitalType::RestrictedOpenShell) { + generalized_fock_mo_.setZero(); + compute_restricted_open_shell_gradient( + scf_impl, C, scf_impl.get_density_matrix(), num_electrons_, + rotation_size_, num_molecular_orbitals, generalized_fock_mo_, + current_gradient_); + } else { + compute_restricted_unrestricted_gradient( + F, C, num_electrons_, rotation_offset_, rotation_size_, + num_orbital_spin_blocks_, num_molecular_orbitals, current_gradient_); + } if (gdm_step_count_ != 0) { // Add new gradient difference to history for all spins @@ -743,8 +870,7 @@ void GDM::iterate(SCFImpl& scf_impl) { Eigen::VectorXd initial_hessian; build_restricted_unrestricted_pseudo_canonical_orbitals_hessian( - F, C, num_density_matrices, num_molecular_orbitals, cfg->scf_orbital_type, - initial_hessian); + F, C, num_molecular_orbitals, initial_hessian); double latest_inverse_rho = 1.0; // BFGS two-loop recursion on concatenated vectors (runs once for all spins) From ea5becb9f792ce4e9ddba51f2a56fa2132a09dc0 Mon Sep 17 00:00:00 2001 From: v-boqinzhang Date: Thu, 26 Feb 2026 14:33:43 -0800 Subject: [PATCH 03/29] add compute_restricted_open_shell_gradient function --- .../microsoft/scf/src/scf_algorithm/gdm.cpp | 155 +++++++++++++++--- 1 file changed, 130 insertions(+), 25 deletions(-) diff --git a/cpp/src/qdk/chemistry/algorithms/microsoft/scf/src/scf_algorithm/gdm.cpp b/cpp/src/qdk/chemistry/algorithms/microsoft/scf/src/scf_algorithm/gdm.cpp index 9e69c182a..a644509cf 100644 --- a/cpp/src/qdk/chemistry/algorithms/microsoft/scf/src/scf_algorithm/gdm.cpp +++ b/cpp/src/qdk/chemistry/algorithms/microsoft/scf/src/scf_algorithm/gdm.cpp @@ -36,10 +36,10 @@ namespace impl { * @param[in] num_occupied_orbitals Number of occupied orbitals for this spin * @param[in] num_molecular_orbitals Number of molecular orbitals */ -static void apply_orbital_rotation(RowMajorMatrix& C, const int spin_index, - const Eigen::VectorXd& kappa_vector, - const int num_occupied_orbitals, - const int num_molecular_orbitals) { +static void apply_restricted_unrestricted_orbital_rotation( + RowMajorMatrix& C, const int spin_index, + const Eigen::VectorXd& kappa_vector, const int num_occupied_orbitals, + const int num_molecular_orbitals) { QDK_LOG_TRACE_ENTERING(); const int num_virtual_orbitals = num_molecular_orbitals - num_occupied_orbitals; @@ -74,6 +74,72 @@ static void apply_orbital_rotation(RowMajorMatrix& C, const int spin_index, num_molecular_orbitals); } +/** + * @brief Construct the ROHF kappa matrix and apply rotation C * exp(kappa) + * + * @param[in,out] C Molecular orbital coefficient matrix + * @param[in] num_electrons Occupied orbital counts (alpha, beta) + * @param[in] kappa_vector Concatenated ROHF rotation parameters + * @param[in] num_molecular_orbitals Total molecular orbitals in the system + */ +static void apply_restricted_open_shell_orbital_rotation( + RowMajorMatrix& C, const std::vector& num_electrons, + const Eigen::VectorXd& kappa_vector, const int num_molecular_orbitals) { + QDK_LOG_TRACE_ENTERING(); + const int num_closed_orbitals = num_electrons[1]; + const int num_open_orbitals = num_electrons[0] - num_closed_orbitals; + const int num_virtual_orbitals = num_molecular_orbitals - num_electrons[0]; + + int offset = 0; + const int iw_size = num_closed_orbitals * num_open_orbitals; + const int wa_size = num_open_orbitals * num_virtual_orbitals; + const int ia_size = num_closed_orbitals * num_virtual_orbitals; + + const RowMajorMatrix kappa_iw = Eigen::Map( + kappa_vector.data() + offset, num_closed_orbitals, num_open_orbitals); + offset += iw_size; + + const RowMajorMatrix kappa_wa = Eigen::Map( + kappa_vector.data() + offset, num_open_orbitals, num_virtual_orbitals); + offset += wa_size; + + const RowMajorMatrix kappa_ia = Eigen::Map( + kappa_vector.data() + offset, num_closed_orbitals, num_virtual_orbitals); + offset += ia_size; + + RowMajorMatrix kappa_complete = + RowMajorMatrix::Zero(num_molecular_orbitals, num_molecular_orbitals); + + kappa_complete.block(0, num_closed_orbitals, num_closed_orbitals, + num_open_orbitals) = kappa_iw / 2.0; + kappa_complete.block(num_closed_orbitals, 0, num_open_orbitals, + num_closed_orbitals) = -kappa_iw.transpose() / 2.0; + + const int open_start = num_closed_orbitals; + const int virtual_start = num_closed_orbitals + num_open_orbitals; + + kappa_complete.block(open_start, virtual_start, num_open_orbitals, + num_virtual_orbitals) = kappa_wa / 2.0; + kappa_complete.block(virtual_start, open_start, num_virtual_orbitals, + num_open_orbitals) = -kappa_wa.transpose() / 2.0; + + kappa_complete.block(0, virtual_start, num_closed_orbitals, + num_virtual_orbitals) = kappa_ia / 2.0; + kappa_complete.block(virtual_start, 0, num_virtual_orbitals, + num_closed_orbitals) = -kappa_ia.transpose() / 2.0; + + RowMajorMatrix exp_kappa = + RowMajorMatrix::Zero(num_molecular_orbitals, num_molecular_orbitals); + matrix_exp(kappa_complete.data(), exp_kappa.data(), num_molecular_orbitals); + + RowMajorMatrix C_before_rotate = C; + blas::gemm(blas::Layout::RowMajor, blas::Op::NoTrans, blas::Op::NoTrans, + num_molecular_orbitals, num_molecular_orbitals, + num_molecular_orbitals, 1.0, C_before_rotate.data(), + num_molecular_orbitals, exp_kappa.data(), num_molecular_orbitals, + 0.0, C.data(), num_molecular_orbitals); +} + /** * @brief Compute restricted/unrestricted orbital gradients for all spins * @param[in] F Fock matrix in AO basis @@ -136,6 +202,21 @@ static void compute_restricted_unrestricted_gradient( } } +/** + * @brief Compute ROHF orbital gradients using the generalized Fock matrix + * + * The gradient is packed as (iw, wa, ia) blocks following the same + * segmentation used by the ROHF kappa vector. + * + * @param[in] scf_impl SCF implementation for J/K construction + * @param[in] C Molecular orbital coefficient matrix + * @param[in] density_matrix Density matrix in AO basis + * @param[in] num_electrons Occupied orbital counts (alpha, beta) + * @param[in] rotation_size Rotation size for the ROHF kappa vector + * @param[in] num_molecular_orbitals Total molecular orbitals in the system + * @param[out] generalized_fock_mo Preallocated generalized Fock matrix in MO + * @param[out] gradient Output gradient vector + */ static void compute_restricted_open_shell_gradient( const SCFImpl& scf_impl, const RowMajorMatrix& C, const RowMajorMatrix& density_matrix, const std::vector& num_electrons, @@ -199,27 +280,29 @@ static void compute_restricted_open_shell_gradient( int offset = 0; RowMajorMatrix grad_iw = - 2.0 * (generalized_fock_mo.block(0, num_closed_orbitals, - num_closed_orbitals, num_open_orbitals) - - generalized_fock_mo - .block(num_closed_orbitals, 0, num_open_orbitals, - num_closed_orbitals) - .transpose()); + -2.0 * + (generalized_fock_mo.block(0, num_closed_orbitals, num_closed_orbitals, + num_open_orbitals) - + generalized_fock_mo + .block(num_closed_orbitals, 0, num_open_orbitals, + num_closed_orbitals) + .transpose()); gradient.segment(offset, num_closed_orbitals * num_open_orbitals) = Eigen::Map(grad_iw.data(), grad_iw.size()).eval(); offset += num_closed_orbitals * num_open_orbitals; RowMajorMatrix grad_wa = - 2.0 * generalized_fock_mo.block(num_closed_orbitals, - num_closed_orbitals + num_open_orbitals, - num_open_orbitals, num_virtual_orbitals); + -2.0 * generalized_fock_mo.block(num_closed_orbitals, + num_closed_orbitals + num_open_orbitals, + num_open_orbitals, num_virtual_orbitals); gradient.segment(offset, num_open_orbitals * num_virtual_orbitals) = Eigen::Map(grad_wa.data(), grad_wa.size()).eval(); offset += num_open_orbitals * num_virtual_orbitals; - RowMajorMatrix grad_ia = 2.0 * generalized_fock_mo.block( - 0, num_closed_orbitals + num_open_orbitals, - num_closed_orbitals, num_virtual_orbitals); + RowMajorMatrix grad_ia = + -2.0 * + generalized_fock_mo.block(0, num_closed_orbitals + num_open_orbitals, + num_closed_orbitals, num_virtual_orbitals); gradient.segment(offset, num_closed_orbitals * num_virtual_orbitals) = Eigen::Map(grad_ia.data(), grad_ia.size()).eval(); } @@ -335,11 +418,16 @@ double GDMLineFunctor::eval(const Eigen::VectorXd& x) { cached_C_ = C_pseudo_canonical_; // Apply rotation for all spins with kappa_trial - for (int i = 0; i < num_orbital_spin_blocks_; i++) { - auto kappa_spin = - kappa_trial.segment(rotation_offset_[i], rotation_size_[i]); - apply_orbital_rotation(cached_C_, i, kappa_spin, num_electrons_[i], - num_molecular_orbitals_); + if (scf_orbital_type_ == SCFOrbitalType::RestrictedOpenShell) { + apply_restricted_open_shell_orbital_rotation( + cached_C_, num_electrons_, kappa_trial, num_molecular_orbitals_); + } else { + for (int i = 0; i < num_orbital_spin_blocks_; i++) { + auto kappa_spin = + kappa_trial.segment(rotation_offset_[i], rotation_size_[i]); + apply_restricted_unrestricted_orbital_rotation( + cached_C_, i, kappa_spin, num_electrons_[i], num_molecular_orbitals_); + } } // Compute P_trial from rotated C (for all spins) @@ -480,10 +568,14 @@ class GDM { * @param[in] num_molecular_orbitals Total number of molecular orbitals * @param[out] initial_hessian Output concatenated initial Hessian */ - void build_restricted_unrestricted_pseudo_canonical_orbitals_hessian( + void build_restricted_unrestricted_pseudo_canonical_orbitals_hessian_( const RowMajorMatrix& F, RowMajorMatrix& C, int num_molecular_orbitals, Eigen::VectorXd& initial_hessian); + void build_restricted_open_shell_pseudo_canonical_orbitals_hessian_( + RowMajorMatrix& C, int num_molecular_orbitals, + Eigen::VectorXd& initial_hessian); + /// Reference to SCFContext const SCFContext& ctx_; ///< Reference to SCFContext /// Energy change from the last step @@ -771,7 +863,7 @@ void GDM::generate_pseudo_canonical_orbital_( current_gradient_spin = current_gradient_transformed; } -void GDM::build_restricted_unrestricted_pseudo_canonical_orbitals_hessian( +void GDM::build_restricted_unrestricted_pseudo_canonical_orbitals_hessian_( const RowMajorMatrix& F, RowMajorMatrix& C, int num_molecular_orbitals, Eigen::VectorXd& initial_hessian) { initial_hessian.setZero(total_rotation_size_); @@ -814,6 +906,14 @@ void GDM::build_restricted_unrestricted_pseudo_canonical_orbitals_hessian( } } +void GDM::build_restricted_open_shell_pseudo_canonical_orbitals_hessian_( + RowMajorMatrix& C, int num_molecular_orbitals, + Eigen::VectorXd& initial_hessian) { + const int num_closed_orbitals = num_electrons_[1]; + const int num_open_orbitals = num_electrons_[0] - num_closed_orbitals; + const int num_virtual_orbitals = num_molecular_orbitals - num_electrons_[0]; +} + void GDM::iterate(SCFImpl& scf_impl) { QDK_LOG_TRACE_ENTERING(); auto& C = scf_impl.orbitals_matrix(); @@ -869,8 +969,13 @@ void GDM::iterate(SCFImpl& scf_impl) { } Eigen::VectorXd initial_hessian; - build_restricted_unrestricted_pseudo_canonical_orbitals_hessian( - F, C, num_molecular_orbitals, initial_hessian); + if (cfg->scf_orbital_type == SCFOrbitalType::RestrictedOpenShell) { + build_restricted_open_shell_pseudo_canonical_orbitals_hessian_( + C, num_molecular_orbitals, initial_hessian); + } else { + build_restricted_unrestricted_pseudo_canonical_orbitals_hessian_( + F, C, num_molecular_orbitals, initial_hessian); + } double latest_inverse_rho = 1.0; // BFGS two-loop recursion on concatenated vectors (runs once for all spins) From e6d06a069dea3b6c00eeba570c53c5812bb91d6a Mon Sep 17 00:00:00 2001 From: v-boqinzhang Date: Thu, 26 Feb 2026 16:01:00 -0800 Subject: [PATCH 04/29] modify GDMLinearFunctor::eval, remove Uoo Uvv from class GDM --- .../microsoft/scf/src/scf_algorithm/gdm.cpp | 83 ++++++++++--------- 1 file changed, 46 insertions(+), 37 deletions(-) diff --git a/cpp/src/qdk/chemistry/algorithms/microsoft/scf/src/scf_algorithm/gdm.cpp b/cpp/src/qdk/chemistry/algorithms/microsoft/scf/src/scf_algorithm/gdm.cpp index a644509cf..f52a78a6e 100644 --- a/cpp/src/qdk/chemistry/algorithms/microsoft/scf/src/scf_algorithm/gdm.cpp +++ b/cpp/src/qdk/chemistry/algorithms/microsoft/scf/src/scf_algorithm/gdm.cpp @@ -95,15 +95,15 @@ static void apply_restricted_open_shell_orbital_rotation( const int wa_size = num_open_orbitals * num_virtual_orbitals; const int ia_size = num_closed_orbitals * num_virtual_orbitals; - const RowMajorMatrix kappa_iw = Eigen::Map( + const auto kappa_iw = Eigen::Map( kappa_vector.data() + offset, num_closed_orbitals, num_open_orbitals); offset += iw_size; - const RowMajorMatrix kappa_wa = Eigen::Map( + const auto kappa_wa = Eigen::Map( kappa_vector.data() + offset, num_open_orbitals, num_virtual_orbitals); offset += wa_size; - const RowMajorMatrix kappa_ia = Eigen::Map( + const auto kappa_ia = Eigen::Map( kappa_vector.data() + offset, num_closed_orbitals, num_virtual_orbitals); offset += ia_size; @@ -111,22 +111,25 @@ static void apply_restricted_open_shell_orbital_rotation( RowMajorMatrix::Zero(num_molecular_orbitals, num_molecular_orbitals); kappa_complete.block(0, num_closed_orbitals, num_closed_orbitals, - num_open_orbitals) = kappa_iw / 2.0; + num_open_orbitals) = kappa_iw; kappa_complete.block(num_closed_orbitals, 0, num_open_orbitals, - num_closed_orbitals) = -kappa_iw.transpose() / 2.0; + num_closed_orbitals) = -kappa_iw.transpose(); const int open_start = num_closed_orbitals; const int virtual_start = num_closed_orbitals + num_open_orbitals; kappa_complete.block(open_start, virtual_start, num_open_orbitals, - num_virtual_orbitals) = kappa_wa / 2.0; + num_virtual_orbitals) = kappa_wa; kappa_complete.block(virtual_start, open_start, num_virtual_orbitals, - num_open_orbitals) = -kappa_wa.transpose() / 2.0; + num_open_orbitals) = -kappa_wa.transpose(); kappa_complete.block(0, virtual_start, num_closed_orbitals, - num_virtual_orbitals) = kappa_ia / 2.0; + num_virtual_orbitals) = kappa_ia; kappa_complete.block(virtual_start, 0, num_virtual_orbitals, - num_closed_orbitals) = -kappa_ia.transpose() / 2.0; + num_closed_orbitals) = -kappa_ia.transpose(); + + // 0.5 is for consistency with the gradient definition + kappa_complete *= 0.5; RowMajorMatrix exp_kappa = RowMajorMatrix::Zero(num_molecular_orbitals, num_molecular_orbitals); @@ -437,17 +440,27 @@ double GDMLineFunctor::eval(const Eigen::VectorXd& x) { for (int i = 0; i < num_density_matrices_; i++) { const int num_occupied_orbitals = num_electrons_[i]; const double occupation_factor = - (scf_orbital_type_ == SCFOrbitalType::Unrestricted) ? 1.0 : 2.0; - - cached_P_.block(num_molecular_orbitals_ * i, 0, num_molecular_orbitals_, - num_molecular_orbitals_) = - occupation_factor * - cached_C_.block(num_molecular_orbitals_ * i, 0, num_molecular_orbitals_, - num_occupied_orbitals) * - cached_C_ - .block(num_molecular_orbitals_ * i, 0, num_molecular_orbitals_, - num_occupied_orbitals) - .transpose(); + (scf_orbital_type_ == SCFOrbitalType::Restricted) ? 2.0 : 1.0; + auto P_block = + cached_P_.block(num_molecular_orbitals_ * i, 0, num_molecular_orbitals_, + num_molecular_orbitals_); + + double* C_block_data = nullptr; + if (scf_orbital_type_ == SCFOrbitalType::RestrictedOpenShell) { + auto C_block = + cached_C_.block(0, 0, num_molecular_orbitals_, num_occupied_orbitals); + C_block_data = C_block.data(); + } else { + auto C_block = + cached_C_.block(num_molecular_orbitals_ * i, 0, + num_molecular_orbitals_, num_occupied_orbitals); + C_block_data = C_block.data(); + } + blas::gemm(blas::Layout::RowMajor, blas::Op::NoTrans, blas::Op::Trans, + num_molecular_orbitals_, num_molecular_orbitals_, + num_occupied_orbitals, occupation_factor, C_block_data, + num_molecular_orbitals_, C_block_data, num_molecular_orbitals_, + 0.0, P_block.data(), num_molecular_orbitals_); } // Evaluate energy and Fock matrix using trial density matrix @@ -613,11 +626,6 @@ class GDM { /// Generalized Fock matrix in MO basis for ROHF RowMajorMatrix generalized_fock_mo_; - /// Horizontal rotation matrix of occupied orbitals - RowMajorMatrix Uoo_; - /// Horizontal rotation matrix of virtual orbitals - RowMajorMatrix Uvv_; - Eigen::VectorXd kappa_; // vertical rotation matrix of this step /// Energy of the last accepted step, used to decide if we rescale the kappa @@ -807,22 +815,23 @@ void GDM::generate_pseudo_canonical_orbital_( // Perform pseudo-canonical transformation and BFGS // Obtain pseudo-canonical orbitals. Foo and Fvv are symmetric matrices, but // the output eigenvectors are column-major - Uoo_ = F_MO.block(0, 0, num_occupied_orbitals, num_occupied_orbitals); - Uvv_ = F_MO.block(num_occupied_orbitals, num_occupied_orbitals, - num_virtual_orbitals, num_virtual_orbitals); + RowMajorMatrix Uoo = + F_MO.block(0, 0, num_occupied_orbitals, num_occupied_orbitals); + RowMajorMatrix Uvv = F_MO.block(num_occupied_orbitals, num_occupied_orbitals, + num_virtual_orbitals, num_virtual_orbitals); // Compute eigenvalues/eigenvectors of occupied-occupied and virtual-virtual // blocks for pseudo-canonical orbital transformation lapack::syev(lapack::Job::Vec, lapack::Uplo::Lower, num_occupied_orbitals, - Uoo_.data(), num_occupied_orbitals, + Uoo.data(), num_occupied_orbitals, pseudo_canonical_eigenvalues_.data()); lapack::syev(lapack::Job::Vec, lapack::Uplo::Lower, num_virtual_orbitals, - Uvv_.data(), num_virtual_orbitals, + Uvv.data(), num_virtual_orbitals, pseudo_canonical_eigenvalues_.data() + num_occupied_orbitals); // Transpose to convert column-major eigenvectors to row-major format - Uoo_.transposeInPlace(); - Uvv_.transposeInPlace(); + Uoo.transposeInPlace(); + Uvv.transposeInPlace(); // Transform occupied orbitals auto C_occ_view = C.block(num_molecular_orbitals * spin_index, 0, @@ -831,7 +840,7 @@ void GDM::generate_pseudo_canonical_orbital_( blas::gemm(blas::Layout::RowMajor, blas::Op::NoTrans, blas::Op::NoTrans, num_molecular_orbitals, num_occupied_orbitals, num_occupied_orbitals, 1.0, C_occ.data(), num_occupied_orbitals, - Uoo_.data(), num_occupied_orbitals, 0.0, C_occ_view.data(), + Uoo.data(), num_occupied_orbitals, 0.0, C_occ_view.data(), num_molecular_orbitals); // Transform virtual orbitals @@ -841,15 +850,15 @@ void GDM::generate_pseudo_canonical_orbital_( RowMajorMatrix C_virt = C_virt_view.eval(); blas::gemm(blas::Layout::RowMajor, blas::Op::NoTrans, blas::Op::NoTrans, num_molecular_orbitals, num_virtual_orbitals, num_virtual_orbitals, - 1.0, C_virt.data(), num_virtual_orbitals, Uvv_.data(), + 1.0, C_virt.data(), num_virtual_orbitals, Uvv.data(), num_virtual_orbitals, 0.0, C_virt_view.data(), num_molecular_orbitals); // Transform the vectors in history_kappa and history_dgrad to // accommodate current pseudo-canonical orbitals - transform_history_(history_kappa_spin, Uoo_, Uvv_, history_size_, + transform_history_(history_kappa_spin, Uoo, Uvv, history_size_, num_occupied_orbitals, num_molecular_orbitals); - transform_history_(history_dgrad_spin, Uoo_, Uvv_, history_size_, + transform_history_(history_dgrad_spin, Uoo, Uvv, history_size_, num_occupied_orbitals, num_molecular_orbitals); // Transform the gradient to accommodate current pseudo-canonical orbitals @@ -857,7 +866,7 @@ void GDM::generate_pseudo_canonical_orbital_( Eigen::Map(current_gradient_spin.data(), num_occupied_orbitals, num_virtual_orbitals); RowMajorMatrix current_gradient_transformed_matrix = - Uoo_.transpose() * current_gradient_matrix * Uvv_; + Uoo.transpose() * current_gradient_matrix * Uvv; Eigen::VectorXd current_gradient_transformed = Eigen::Map( current_gradient_transformed_matrix.data(), rotation_size); current_gradient_spin = current_gradient_transformed; From 0237d01939dabeb9127b57ed115707d18fbcbe3b Mon Sep 17 00:00:00 2001 From: v-boqinzhang Date: Tue, 14 Apr 2026 17:31:57 -0700 Subject: [PATCH 05/29] generate_restricted_open_shell_pseudo_canonical_orbital_ --- .../microsoft/scf/src/scf_algorithm/gdm.cpp | 229 ++++++++++++++---- 1 file changed, 179 insertions(+), 50 deletions(-) diff --git a/cpp/src/qdk/chemistry/algorithms/microsoft/scf/src/scf_algorithm/gdm.cpp b/cpp/src/qdk/chemistry/algorithms/microsoft/scf/src/scf_algorithm/gdm.cpp index f52a78a6e..1dd431b8c 100644 --- a/cpp/src/qdk/chemistry/algorithms/microsoft/scf/src/scf_algorithm/gdm.cpp +++ b/cpp/src/qdk/chemistry/algorithms/microsoft/scf/src/scf_algorithm/gdm.cpp @@ -538,22 +538,21 @@ class GDM { private: /** * @brief Transform history matrices (either history_dgrad or history_kappa) - * using current rotation matrices Uoo and Uvv to transform into the - * pseudo-canonical orbital basis, K_new = Uoo^T * K_old * Uvv + * using current rotation matrices to transform into the + * pseudo-canonical orbital basis, K_new = U_left^T * K_old * U_right * @param[in,out] history History matrix block to be transformed (either * history_dgrad or history_kappa) - * @param[in] u_left Left rotation matrix (e.g., Uoo or Uaa) - * @param[in] u_right Right rotation matrix (e.g., Uvv or Uaa) + * @param[in] u_left Left rotation matrix (e.g., Uii or Uaa) + * @param[in] u_right Right rotation matrix (e.g., Uaa or Uww) * @param[in] history_size Number of history entries - * @param[in] num_occupied_orbitals Number of electrons for current spin - * @param[in] num_molecular_orbitals Number of molecular orbitals + * @param[in] num_rows Number of rows in each unpacked history block + * @param[in] num_cols Number of cols in each unpacked history block * */ void transform_history_(Eigen::Block& history, const RowMajorMatrix& u_left, const RowMajorMatrix& u_right, const int history_size, - const int num_occupied_orbitals, - const int num_molecular_orbitals); + const int num_rows, const int num_cols); /** * @brief Generate pseudo-canonical orbitals and apply transformations @@ -568,12 +567,16 @@ class GDM { * for this spin * */ - void generate_pseudo_canonical_orbital_( + void generate_restricted_unrestricted_pseudo_canonical_orbital_( const RowMajorMatrix& F, RowMajorMatrix& C, const int spin_index, Eigen::Block history_kappa_spin, Eigen::Block history_dgrad_spin, Eigen::VectorBlock current_gradient_spin); + void generate_restricted_open_shell_pseudo_canonical_orbital_( + const RowMajorMatrix& F, RowMajorMatrix& C, RowMajorMatrix& history_kappa, + RowMajorMatrix& history_dgrad, Eigen::VectorXd& current_gradient); + /** * @brief Build pseudo-canonical orbitals and initial Hessian across spins * @param[in] F Fock matrix in AO basis @@ -586,7 +589,7 @@ class GDM { Eigen::VectorXd& initial_hessian); void build_restricted_open_shell_pseudo_canonical_orbitals_hessian_( - RowMajorMatrix& C, int num_molecular_orbitals, + const RowMajorMatrix& F, RowMajorMatrix& C, int num_molecular_orbitals, Eigen::VectorXd& initial_hessian); /// Reference to SCFContext @@ -744,43 +747,34 @@ void GDM::transform_history_(Eigen::Block& history, const RowMajorMatrix& u_left, const RowMajorMatrix& u_right, const int history_size, - const int num_occupied_orbitals, - const int num_molecular_orbitals) { + const int num_rows, + const int num_cols) { QDK_LOG_TRACE_ENTERING(); - const int num_virtual_orbitals = - num_molecular_orbitals - num_occupied_orbitals; // Validate dimensions (negative values indicate invalid input) - if (num_occupied_orbitals < 0 || num_virtual_orbitals < 0) { + if (num_rows < 0 || num_cols < 0) { throw std::invalid_argument( - std::string( - "transform_history_: invalid dimensions (num_occupied_orbitals=") + - std::to_string(num_occupied_orbitals) + - ", num_virtual_orbitals=" + std::to_string(num_virtual_orbitals) + ")"); + std::string("transform_history_: invalid dimensions (num_rows=") + + std::to_string(num_rows) + ", num_cols=" + + std::to_string(num_cols) + ")"); } - // Skip transformation if either dimension is zero (no rotations for this - // spin) - if (num_occupied_orbitals == 0 || num_virtual_orbitals == 0) { + // Skip transformation if either dimension is zero + if (num_rows == 0 || num_cols == 0) { return; } - RowMajorMatrix temp = - RowMajorMatrix::Zero(num_occupied_orbitals, num_virtual_orbitals); + RowMajorMatrix temp = RowMajorMatrix::Zero(num_rows, num_cols); for (int line = 0; line < history_size; line++) { double* history_line_ptr = history.row(line).data(); - // K_ov (new) = U_left(oo)^T * K_ov * U_right(vv) + // K_new = U_left^T * K_old * U_right blas::gemm(blas::Layout::RowMajor, blas::Op::NoTrans, blas::Op::NoTrans, - num_occupied_orbitals, num_virtual_orbitals, - num_virtual_orbitals, 1.0, history_line_ptr, - num_virtual_orbitals, u_right.data(), num_virtual_orbitals, 0.0, - temp.data(), num_virtual_orbitals); + num_rows, num_cols, num_cols, 1.0, history_line_ptr, num_cols, + u_right.data(), num_cols, 0.0, temp.data(), num_cols); blas::gemm(blas::Layout::RowMajor, blas::Op::Trans, blas::Op::NoTrans, - num_occupied_orbitals, num_virtual_orbitals, - num_occupied_orbitals, 1.0, u_left.data(), num_occupied_orbitals, - temp.data(), num_virtual_orbitals, 0.0, history_line_ptr, - num_virtual_orbitals); + num_rows, num_cols, num_rows, 1.0, u_left.data(), num_rows, + temp.data(), num_cols, 0.0, history_line_ptr, num_cols); } } -void GDM::generate_pseudo_canonical_orbital_( +void GDM::generate_restricted_unrestricted_pseudo_canonical_orbital_( const RowMajorMatrix& F, RowMajorMatrix& C, const int spin_index, Eigen::Block history_kappa_spin, Eigen::Block history_dgrad_spin, @@ -815,23 +809,23 @@ void GDM::generate_pseudo_canonical_orbital_( // Perform pseudo-canonical transformation and BFGS // Obtain pseudo-canonical orbitals. Foo and Fvv are symmetric matrices, but // the output eigenvectors are column-major - RowMajorMatrix Uoo = + RowMajorMatrix Uii = F_MO.block(0, 0, num_occupied_orbitals, num_occupied_orbitals); - RowMajorMatrix Uvv = F_MO.block(num_occupied_orbitals, num_occupied_orbitals, + RowMajorMatrix Uaa = F_MO.block(num_occupied_orbitals, num_occupied_orbitals, num_virtual_orbitals, num_virtual_orbitals); // Compute eigenvalues/eigenvectors of occupied-occupied and virtual-virtual // blocks for pseudo-canonical orbital transformation lapack::syev(lapack::Job::Vec, lapack::Uplo::Lower, num_occupied_orbitals, - Uoo.data(), num_occupied_orbitals, + Uii.data(), num_occupied_orbitals, pseudo_canonical_eigenvalues_.data()); lapack::syev(lapack::Job::Vec, lapack::Uplo::Lower, num_virtual_orbitals, - Uvv.data(), num_virtual_orbitals, + Uaa.data(), num_virtual_orbitals, pseudo_canonical_eigenvalues_.data() + num_occupied_orbitals); // Transpose to convert column-major eigenvectors to row-major format - Uoo.transposeInPlace(); - Uvv.transposeInPlace(); + Uii.transposeInPlace(); + Uaa.transposeInPlace(); // Transform occupied orbitals auto C_occ_view = C.block(num_molecular_orbitals * spin_index, 0, @@ -840,7 +834,7 @@ void GDM::generate_pseudo_canonical_orbital_( blas::gemm(blas::Layout::RowMajor, blas::Op::NoTrans, blas::Op::NoTrans, num_molecular_orbitals, num_occupied_orbitals, num_occupied_orbitals, 1.0, C_occ.data(), num_occupied_orbitals, - Uoo.data(), num_occupied_orbitals, 0.0, C_occ_view.data(), + Uii.data(), num_occupied_orbitals, 0.0, C_occ_view.data(), num_molecular_orbitals); // Transform virtual orbitals @@ -850,23 +844,23 @@ void GDM::generate_pseudo_canonical_orbital_( RowMajorMatrix C_virt = C_virt_view.eval(); blas::gemm(blas::Layout::RowMajor, blas::Op::NoTrans, blas::Op::NoTrans, num_molecular_orbitals, num_virtual_orbitals, num_virtual_orbitals, - 1.0, C_virt.data(), num_virtual_orbitals, Uvv.data(), + 1.0, C_virt.data(), num_virtual_orbitals, Uaa.data(), num_virtual_orbitals, 0.0, C_virt_view.data(), num_molecular_orbitals); // Transform the vectors in history_kappa and history_dgrad to // accommodate current pseudo-canonical orbitals - transform_history_(history_kappa_spin, Uoo, Uvv, history_size_, - num_occupied_orbitals, num_molecular_orbitals); - transform_history_(history_dgrad_spin, Uoo, Uvv, history_size_, - num_occupied_orbitals, num_molecular_orbitals); + transform_history_(history_kappa_spin, Uii, Uaa, history_size_, + num_occupied_orbitals, num_virtual_orbitals); + transform_history_(history_dgrad_spin, Uii, Uaa, history_size_, + num_occupied_orbitals, num_virtual_orbitals); // Transform the gradient to accommodate current pseudo-canonical orbitals RowMajorMatrix current_gradient_matrix = Eigen::Map(current_gradient_spin.data(), num_occupied_orbitals, num_virtual_orbitals); RowMajorMatrix current_gradient_transformed_matrix = - Uoo.transpose() * current_gradient_matrix * Uvv; + Uii.transpose() * current_gradient_matrix * Uaa; Eigen::VectorXd current_gradient_transformed = Eigen::Map( current_gradient_transformed_matrix.data(), rotation_size); current_gradient_spin = current_gradient_transformed; @@ -890,7 +884,7 @@ void GDM::build_restricted_unrestricted_pseudo_canonical_orbitals_hessian_( current_gradient_.segment(rotation_offset_[i], rotation_size_[i]); // Generate pseudo-canonical orbitals and transform gradient and history - generate_pseudo_canonical_orbital_( + generate_restricted_unrestricted_pseudo_canonical_orbital_( F, C, i, history_kappa_spin, history_dgrad_spin, current_gradient_spin); // Build this spin's segment of initial Hessian @@ -915,8 +909,143 @@ void GDM::build_restricted_unrestricted_pseudo_canonical_orbitals_hessian_( } } +static void diagonalize_block_rotate_orbitals( + RowMajorMatrix& input_block_output_eigenvectors, const int block_size, + Eigen::VectorXd& eigenvalues, const int eigenvalue_start_index, + Eigen::Block transformed_orbitals, + const int num_atomic_orbitals, const int num_molecular_orbitals) { + lapack::syev(lapack::Job::Vec, lapack::Uplo::Lower, block_size, + input_block_output_eigenvectors.data(), block_size, + eigenvalues.data() + eigenvalue_start_index); + input_block_output_eigenvectors.transposeInPlace(); + RowMajorMatrix copied_orbitals = transformed_orbitals; + blas::gemm(blas::Layout::RowMajor, blas::Op::NoTrans, blas::Op::NoTrans, + num_atomic_orbitals, block_size, block_size, 1.0, + copied_orbitals.data(), block_size, + input_block_output_eigenvectors.data(), block_size, 0.0, + transformed_orbitals.data(), num_molecular_orbitals); +} + +static void transform_gradient_vector(const Eigen::VectorXd& current_vector, + const int start_index, const int num_rows, + const int num_cols, const RowMajorMatrix& u_left, + const RowMajorMatrix& u_right, + Eigen::VectorXd& transformed_vector) { + RowMajorMatrix current_matrix = Eigen::Map( + current_vector.data() + start_index, num_rows, num_cols); + RowMajorMatrix transformed_matrix = u_left.transpose() * current_matrix * u_right; + transformed_vector.segment(start_index, num_rows * num_cols) = + Eigen::Map(transformed_matrix.data(), + num_rows * num_cols); +} + +void GDM::generate_restricted_open_shell_pseudo_canonical_orbital_( + const RowMajorMatrix& F, RowMajorMatrix& C, RowMajorMatrix& history_kappa, + RowMajorMatrix& history_dgrad, Eigen::VectorXd& current_gradient) { + const int num_molecular_orbitals = C.cols(); + const int num_atomic_orbitals = C.rows(); + const int num_closed_orbitals = num_electrons_[1]; + const int num_open_orbitals = num_electrons_[0] - num_closed_orbitals; + const int num_occupied_orbitals = num_electrons_[0]; + const int num_virtual_orbitals = num_molecular_orbitals - num_electrons_[0]; + const int size_closed_open = num_closed_orbitals * num_open_orbitals; + const int size_open_virtual = num_open_orbitals * num_virtual_orbitals; + const int size_closed_virtual = num_closed_orbitals * num_virtual_orbitals; + + RowMajorMatrix F_up_mo = + RowMajorMatrix::Zero(num_molecular_orbitals, num_molecular_orbitals); + RowMajorMatrix F_dn_mo = F_up_mo; + + F_up_mo.noalias() = C.transpose() * F.block(0, 0, num_atomic_orbitals, num_molecular_orbitals) * C; + F_dn_mo.noalias() = C.transpose() * F.block(num_atomic_orbitals, 0, num_atomic_orbitals, + num_molecular_orbitals) * C; + RowMajorMatrix Uii = + F_up_mo.block(0, 0, num_closed_orbitals, num_closed_orbitals) + F_dn_mo.block(0, 0, num_closed_orbitals, num_closed_orbitals); + RowMajorMatrix Uww = F_up_mo.block(num_closed_orbitals, num_closed_orbitals, + num_open_orbitals, num_open_orbitals) + F_dn_mo.block(num_closed_orbitals, num_closed_orbitals, + num_open_orbitals, num_open_orbitals); + RowMajorMatrix Uaa = F_up_mo.block(num_occupied_orbitals, num_occupied_orbitals, + num_virtual_orbitals, num_virtual_orbitals) + F_dn_mo.block(num_occupied_orbitals, num_occupied_orbitals, + num_virtual_orbitals, num_virtual_orbitals); + + // Compute eigenvalues/eigenvectors of occupied-occupied and virtual-virtual + // blocks for pseudo-canonical orbital transformation + // Then transform the orbitals to pseudo-canonical basis + auto C_closed_view = + C.block(0, 0, num_atomic_orbitals, num_closed_orbitals); + diagonalize_block_rotate_orbitals( + Uii, num_closed_orbitals, pseudo_canonical_eigenvalues_, 0, + C_closed_view, num_atomic_orbitals, num_molecular_orbitals); + + auto C_open_view = C.block(0, num_closed_orbitals, num_atomic_orbitals, + num_open_orbitals); + diagonalize_block_rotate_orbitals( + Uww, num_open_orbitals, pseudo_canonical_eigenvalues_, + num_closed_orbitals, C_open_view, num_atomic_orbitals, + num_molecular_orbitals); + + auto C_virt_view = C.block(0, num_occupied_orbitals, num_atomic_orbitals, + num_virtual_orbitals); + diagonalize_block_rotate_orbitals( + Uaa, num_virtual_orbitals, pseudo_canonical_eigenvalues_, + num_occupied_orbitals, C_virt_view, num_atomic_orbitals, + num_molecular_orbitals); + + // Transform the vectors in history_kappa and history_dgrad to + // accommodate current pseudo-canonical orbitals + int offset = 0; + auto history_kappa_iw = + history_kappa.block(0, 0, history_size_, size_closed_open); + transform_history_(history_kappa_iw, Uii, Uww, history_size_, + num_closed_orbitals, num_open_orbitals); + auto history_dgrad_iw = + history_dgrad.block(0, 0, history_size_, size_closed_open); + transform_history_(history_dgrad_iw, Uii, Uww, history_size_, + num_closed_orbitals, num_open_orbitals); + offset += size_closed_open; + + auto history_kappa_wa = + history_kappa.block(0, offset, history_size_, size_open_virtual); + transform_history_(history_kappa_wa, Uww, Uaa, history_size_, + num_open_orbitals, num_virtual_orbitals); + auto history_dgrad_wa = + history_dgrad.block(0, offset, history_size_, size_open_virtual); + transform_history_(history_dgrad_wa, Uww, Uaa, history_size_, + num_open_orbitals, num_virtual_orbitals); + offset += size_open_virtual; + + auto history_kappa_ia = + history_kappa.block(0, offset, history_size_, size_closed_virtual); + transform_history_(history_kappa_ia, Uii, Uaa, history_size_, + num_closed_orbitals, num_virtual_orbitals); + auto history_dgrad_ia = + history_dgrad.block(0, offset, history_size_, size_closed_virtual); + transform_history_(history_dgrad_ia, Uii, Uaa, history_size_, + num_closed_orbitals, num_virtual_orbitals); + + // Transform the gradient to accommodate current pseudo-canonical orbitals + Eigen::VectorXd current_gradient_transformed = Eigen::VectorXd(current_gradient.size()); + + offset = 0; + transform_gradient_vector(current_gradient, offset, num_closed_orbitals, + num_open_orbitals, Uii, Uww, + current_gradient_transformed); + offset += size_closed_open; + + transform_gradient_vector(current_gradient, offset, num_open_orbitals, + num_virtual_orbitals, Uww, Uaa, + current_gradient_transformed); + offset += size_open_virtual; + + transform_gradient_vector(current_gradient, offset, num_closed_orbitals, + num_virtual_orbitals, Uii, Uaa, + current_gradient_transformed); + + current_gradient = current_gradient_transformed; +} + void GDM::build_restricted_open_shell_pseudo_canonical_orbitals_hessian_( - RowMajorMatrix& C, int num_molecular_orbitals, + const RowMajorMatrix& F, RowMajorMatrix& C, int num_molecular_orbitals, Eigen::VectorXd& initial_hessian) { const int num_closed_orbitals = num_electrons_[1]; const int num_open_orbitals = num_electrons_[0] - num_closed_orbitals; @@ -980,7 +1109,7 @@ void GDM::iterate(SCFImpl& scf_impl) { Eigen::VectorXd initial_hessian; if (cfg->scf_orbital_type == SCFOrbitalType::RestrictedOpenShell) { build_restricted_open_shell_pseudo_canonical_orbitals_hessian_( - C, num_molecular_orbitals, initial_hessian); + F, C, num_molecular_orbitals, initial_hessian); } else { build_restricted_unrestricted_pseudo_canonical_orbitals_hessian_( F, C, num_molecular_orbitals, initial_hessian); From eea787b0829d5e00af2a3d614e42eb8a8e987362 Mon Sep 17 00:00:00 2001 From: v-boqinzhang Date: Tue, 21 Apr 2026 11:50:52 -0700 Subject: [PATCH 06/29] ROHF initial Hessian, move function build_rohf_f_p_matrix --- .../qdk/chemistry/scf/core/scf_algorithm.h | 54 +++++ .../microsoft/scf/src/scf_algorithm/diis.cpp | 189 +--------------- .../microsoft/scf/src/scf_algorithm/diis.h | 42 +--- .../microsoft/scf/src/scf_algorithm/gdm.cpp | 170 ++++++++++----- .../scf/src/scf_algorithm/scf_algorithm.cpp | 201 ++++++++++++++++-- 5 files changed, 362 insertions(+), 294 deletions(-) diff --git a/cpp/src/qdk/chemistry/algorithms/microsoft/scf/include/qdk/chemistry/scf/core/scf_algorithm.h b/cpp/src/qdk/chemistry/algorithms/microsoft/scf/include/qdk/chemistry/scf/core/scf_algorithm.h index 96f7d7884..120333bb8 100644 --- a/cpp/src/qdk/chemistry/algorithms/microsoft/scf/include/qdk/chemistry/scf/core/scf_algorithm.h +++ b/cpp/src/qdk/chemistry/algorithms/microsoft/scf/include/qdk/chemistry/scf/core/scf_algorithm.h @@ -113,6 +113,23 @@ class SCFAlgorithm { bool unrestricted, int nelec_alpha, int nelec_beta); + /** + * @brief Try to provide ROHF convergence matrices for OG evaluation + * + * For ROHF, some algorithms evaluate convergence using an effective Fock + * matrix and total density matrix rather than the spin-blocked SCFImpl + * matrices. Derived algorithms can override this hook and provide pointers + * to those matrices. + * + * @param[in] scf_impl SCF implementation object + * @param[out] fock_matrix Pointer to convergence Fock matrix + * @param[out] density_matrix Pointer to convergence density matrix + * @return true if matrices were provided by the algorithm, false otherwise + */ + virtual bool try_get_rohf_convergence_matrices( + const SCFImpl& scf_impl, const RowMajorMatrix*& fock_matrix, + const RowMajorMatrix*& density_matrix); + /** * @brief Calculate orbital gradient (OG) error for convergence checking * @@ -137,6 +154,41 @@ class SCFAlgorithm { RowMajorMatrix& error_matrix, int num_orbital_sets); + /** + * @brief Build ROHF convergence matrices from spin-blocked SCF matrices + * + * Converts spin-blocked Fock/density matrices into the effective ROHF Fock + * and total-density representation used for OG evaluation. + */ + static void build_rohf_f_p_matrix(const RowMajorMatrix& F, + const RowMajorMatrix& C, + const RowMajorMatrix& P, int nelec_alpha, + int nelec_beta, + RowMajorMatrix& effective_fock, + RowMajorMatrix& total_density); + + /** + * @brief Refresh cached ROHF convergence matrices from current SCF state + * + * @return true if ROHF matrices were refreshed, false for non-ROHF runs + */ + bool ensure_rohf_convergence_matrices_(const SCFImpl& scf_impl); + + /** + * @brief Access cached ROHF effective Fock matrix + */ + const RowMajorMatrix& get_rohf_convergence_fock_matrix() const; + + /** + * @brief Access cached ROHF total density matrix + */ + const RowMajorMatrix& get_rohf_convergence_density_matrix() const; + + /** + * @brief Mutable access to cached ROHF total density matrix + */ + RowMajorMatrix& rohf_convergence_density_matrix(); + protected: const SCFContext& ctx_; ///< Reference to SCF context double og_error_ = 0.0; ///< Current orbital gradient error @@ -149,5 +201,7 @@ class SCFAlgorithm { double delta_energy_ = std::numeric_limits::infinity(); ///< Energy change double density_rms_ = 0.0; ///< Last calculated density RMS + RowMajorMatrix rohf_effective_fock_; ///< Cached ROHF effective Fock (AO) + RowMajorMatrix rohf_total_density_; ///< Cached ROHF total density (P_a+P_b) }; } // namespace qdk::chemistry::scf diff --git a/cpp/src/qdk/chemistry/algorithms/microsoft/scf/src/scf_algorithm/diis.cpp b/cpp/src/qdk/chemistry/algorithms/microsoft/scf/src/scf_algorithm/diis.cpp index 67cf66f85..8d3555084 100644 --- a/cpp/src/qdk/chemistry/algorithms/microsoft/scf/src/scf_algorithm/diis.cpp +++ b/cpp/src/qdk/chemistry/algorithms/microsoft/scf/src/scf_algorithm/diis.cpp @@ -6,13 +6,10 @@ #include -#include #include -#include #include #include #include -#include #include "../scf/scf_impl.h" #include "util/macros.h" @@ -121,132 +118,6 @@ class DIIS { std::numeric_limits::infinity(); ///< Current DIIS error }; -/** - * @brief Rebuild the effective ROHF Fock and density matrices - * - * Implements the averaged-block construction (Guest & Saunders 1974, - * Plakhutin & Davidson 2014) to obtain the single-density ROHF view and - * caches the total density $P_\alpha + P_\beta$ alongside it. - * - * @param F Spin-blocked Fock matrix with alpha and beta blocks stacked - * @param C Molecular-orbital coefficients used for block transformations - * @param P Spin-blocked density matrix - * @param nelec_alpha Number of alpha electrons - * @param nelec_beta Number of beta electrons - * @param effective_fock Output effective ROHF Fock matrix (AO basis) - * @param total_density Output total density matrix (P_alpha + P_beta) - */ -void build_rohf_f_p_matrix(const RowMajorMatrix& F, const RowMajorMatrix& C, - const RowMajorMatrix& P, int nelec_alpha, - int nelec_beta, RowMajorMatrix& effective_fock, - RowMajorMatrix& total_density) { - QDK_LOG_TRACE_ENTERING(); - const int num_atomic_orbitals = static_cast(C.rows()); - const int num_molecular_orbitals = static_cast(C.cols()); - if (num_atomic_orbitals != num_molecular_orbitals) { - throw std::runtime_error( - "ROHF build requires number of atomic orbitals to equal number of " - "molecular orbitals!"); - } - - total_density = - P.block(0, 0, num_atomic_orbitals, num_atomic_orbitals) + - P.block(num_atomic_orbitals, 0, num_atomic_orbitals, num_atomic_orbitals); - - if (effective_fock.rows() != num_atomic_orbitals || - effective_fock.cols() != num_atomic_orbitals) { - effective_fock = - RowMajorMatrix::Zero(num_atomic_orbitals, num_atomic_orbitals); - } - - if (C.isZero()) { - effective_fock.noalias() = - F.block(0, 0, num_atomic_orbitals, num_atomic_orbitals); - return; - } - - RowMajorMatrix F_up_mo = - RowMajorMatrix::Zero(num_molecular_orbitals, num_molecular_orbitals); - RowMajorMatrix F_dn_mo = F_up_mo; - RowMajorMatrix effective_F_mo = F_up_mo; - - F_up_mo.noalias() = C.transpose() * - F.block(0, 0, num_atomic_orbitals, num_atomic_orbitals) * - C; - F_dn_mo.noalias() = C.transpose() * - F.block(num_atomic_orbitals, 0, num_atomic_orbitals, - num_atomic_orbitals) * - C; - - auto average_block = [&effective_F_mo, &F_up_mo, &F_dn_mo]( - int row, int col, int rows, int cols) { - if (rows <= 0 || cols <= 0) return; - effective_F_mo.block(row, col, rows, cols).noalias() = - 0.5 * (F_up_mo.block(row, col, rows, cols) + - F_dn_mo.block(row, col, rows, cols)); - }; - auto copy_block = [&effective_F_mo](const RowMajorMatrix& src, int row, - int col, int rows, int cols) { - if (rows <= 0 || cols <= 0) return; - effective_F_mo.block(row, col, rows, cols) = - src.block(row, col, rows, cols); - }; - - const int nd = nelec_beta; - const int ns = nelec_alpha - nelec_beta; - const int nv = num_molecular_orbitals - nelec_alpha; - - average_block(0, 0, nd, nd); - average_block(0, nd + ns, nd, nv); - average_block(nd + ns, 0, nv, nd); - average_block(nd + ns, nd + ns, nv, nv); - average_block(nd, nd, ns, ns); - copy_block(F_dn_mo, 0, nd, nd, ns); - copy_block(F_dn_mo, nd, 0, ns, nd); - copy_block(F_up_mo, nd, nd + ns, ns, nv); - copy_block(F_up_mo, nd + ns, nd, nv, ns); - - // Transform the effective Fock matrix back to AO basis by solving - // C^{-T} * F_MO * C^{-1} = F_AO - // We use LAPACK's getrf/getrs to solve the linear systems involving C^T and - // C without explicitly inverting C - const int matrix_dim = num_molecular_orbitals; - using ColMajorMatrix = - Eigen::Matrix; - // LAPACK expects column-major layout, so we copy the row-major data into a - // column-major matrix without transposing the logical layout - ColMajorMatrix Ct = - Eigen::Map(C.data(), matrix_dim, C.rows()); - // F_MO is symmetric, so we can use it directly as the right-hand side - // without transposing - ColMajorMatrix temp_rhs = effective_F_mo; - std::vector ipiv(matrix_dim); - - auto info = - lapack::getrf(matrix_dim, matrix_dim, Ct.data(), matrix_dim, ipiv.data()); - if (info != 0) { - throw std::runtime_error("getrf failed while factorizing C^T"); - } - - info = lapack::getrs(lapack::Op::NoTrans, matrix_dim, matrix_dim, Ct.data(), - matrix_dim, ipiv.data(), temp_rhs.data(), matrix_dim); - if (info != 0) { - throw std::runtime_error("getrs failed while solving C^T X = F_mo"); - } - - temp_rhs.transposeInPlace(); - info = lapack::getrs(lapack::Op::NoTrans, matrix_dim, matrix_dim, Ct.data(), - matrix_dim, ipiv.data(), temp_rhs.data(), matrix_dim); - if (info != 0) { - throw std::runtime_error("getrs failed while solving C^T X = M^T"); - } - - effective_fock = temp_rhs.transpose(); - if (!effective_fock.isApprox(effective_fock.transpose())) { - effective_fock = 0.5 * (effective_fock + effective_fock.transpose().eval()); - } -} - /** * @brief Reconstruct spin-blocked densities from the ROHF MO matrix * @@ -430,14 +301,6 @@ DIIS::DIIS(const SCFContext& ctx, const size_t subspace_size) : SCFAlgorithm(ctx), diis_impl_(std::make_unique(ctx, subspace_size)) { QDK_LOG_TRACE_ENTERING(); - if (ctx.cfg->scf_orbital_type == SCFOrbitalType::RestrictedOpenShell) { - const int num_atomic_orbitals = - static_cast(ctx.basis_set->num_atomic_orbitals); - rohf_effective_fock_ = - RowMajorMatrix::Zero(num_atomic_orbitals, num_atomic_orbitals); - rohf_total_density_ = - RowMajorMatrix::Zero(num_atomic_orbitals, num_atomic_orbitals); - } } DIIS::~DIIS() noexcept = default; @@ -502,42 +365,6 @@ void DIIS::update_density_matrix(RowMajorMatrix& P, const RowMajorMatrix& C, nelec_beta); } -void DIIS::build_rohf_f_p_matrix(const RowMajorMatrix& F, - const RowMajorMatrix& C, - const RowMajorMatrix& P, int nelec_alpha, - int nelec_beta) { - QDK_LOG_TRACE_ENTERING(); - if (ctx_.cfg->scf_orbital_type != SCFOrbitalType::RestrictedOpenShell) { - throw std::logic_error("ROHF matrix build requested for non-ROHF run"); - } - impl::build_rohf_f_p_matrix(F, C, P, nelec_alpha, nelec_beta, - rohf_effective_fock_, rohf_total_density_); -} - -const RowMajorMatrix& DIIS::get_rohf_fock_matrix() const { - QDK_LOG_TRACE_ENTERING(); - if (rohf_effective_fock_.size() == 0) { - throw std::logic_error("ROHF cache not initialized"); - } - return rohf_effective_fock_; -} - -const RowMajorMatrix& DIIS::get_rohf_density_matrix() const { - QDK_LOG_TRACE_ENTERING(); - if (rohf_total_density_.size() == 0) { - throw std::logic_error("ROHF cache not initialized"); - } - return rohf_total_density_; -} - -RowMajorMatrix& DIIS::rohf_density_matrix() { - QDK_LOG_TRACE_ENTERING(); - if (rohf_total_density_.size() == 0) { - throw std::logic_error("ROHF cache not initialized"); - } - return rohf_total_density_; -} - double DIIS::current_diis_error() const { QDK_LOG_TRACE_ENTERING(); return diis_impl_->get_diis_error(); @@ -546,19 +373,21 @@ double DIIS::current_diis_error() const { RowMajorMatrix& DIIS::select_working_density(SCFImpl& scf_impl) { QDK_LOG_TRACE_ENTERING(); if (ctx_.cfg->scf_orbital_type == SCFOrbitalType::RestrictedOpenShell) { - // The ROHF helper is refreshed inside SCFAlgorithm::check_convergence(), - // so by the time iterate() is invoked the total density view is already - // up to date. - (void)scf_impl; // kept for symmetry with the unrestricted branch - return rohf_density_matrix(); + if (!ensure_rohf_convergence_matrices_(scf_impl)) { + throw std::logic_error("ROHF convergence matrices are unavailable"); + } + return rohf_convergence_density_matrix(); } return scf_impl.density_matrix(); } -const RowMajorMatrix& DIIS::select_working_fock(const SCFImpl& scf_impl) { +const RowMajorMatrix& DIIS::select_working_fock(SCFImpl& scf_impl) { QDK_LOG_TRACE_ENTERING(); if (ctx_.cfg->scf_orbital_type == SCFOrbitalType::RestrictedOpenShell) { - return get_rohf_fock_matrix(); + if (!ensure_rohf_convergence_matrices_(scf_impl)) { + throw std::logic_error("ROHF convergence matrices are unavailable"); + } + return get_rohf_convergence_fock_matrix(); } return scf_impl.get_fock_matrix(); } diff --git a/cpp/src/qdk/chemistry/algorithms/microsoft/scf/src/scf_algorithm/diis.h b/cpp/src/qdk/chemistry/algorithms/microsoft/scf/src/scf_algorithm/diis.h index 6b1628d7f..0dc792eea 100644 --- a/cpp/src/qdk/chemistry/algorithms/microsoft/scf/src/scf_algorithm/diis.h +++ b/cpp/src/qdk/chemistry/algorithms/microsoft/scf/src/scf_algorithm/diis.h @@ -69,44 +69,6 @@ class DIIS : public SCFAlgorithm { bool unrestricted, int nelec_alpha, int nelec_beta) override; - /** - * @brief Build cached ROHF Fock/density matrices for convergence checks - * - * Converts spin-blocked Fock/density matrices into the total-density / - * effective-Fock representation, then saved internally for use in DIIS - * iterations and convergence checks - * - * @param[in] F Spin-blocked Fock matrix from `SCFImpl` - * @param[in] C Molecular orbital coefficients - * @param[in] P Spin-blocked density matrix - * @param[in] nelec_alpha Number of alpha electrons - * @param[in] nelec_beta Number of beta electrons - */ - void build_rohf_f_p_matrix(const RowMajorMatrix& F, const RowMajorMatrix& C, - const RowMajorMatrix& P, int nelec_alpha, - int nelec_beta); - - /** - * @brief Access the cached ROHF-effective Fock matrix - * - * @return Const reference to the total-density ROHF Fock matrix - */ - const RowMajorMatrix& get_rohf_fock_matrix() const; - - /** - * @brief Access the cached total (alpha+beta) density matrix - * - * @return Const reference to the cached ROHF density - */ - const RowMajorMatrix& get_rohf_density_matrix() const; - - /** - * @brief Mutable access to the total density matrix (used inside iterate) - * - * @return Reference to the cached ROHF density - */ - RowMajorMatrix& rohf_density_matrix(); - private: /** * @brief Return most recent Pulay error metric for damping logic @@ -120,11 +82,9 @@ class DIIS : public SCFAlgorithm { /** * @brief Select the Fock matrix view used for the current iteration */ - const RowMajorMatrix& select_working_fock(const SCFImpl& scf_impl); + const RowMajorMatrix& select_working_fock(SCFImpl& scf_impl); std::unique_ptr diis_impl_; ///< Pulay DIIS core implementation - RowMajorMatrix rohf_effective_fock_; ///< Cached ROHF effective Fock (AO) - RowMajorMatrix rohf_total_density_; ///< Cached ROHF total density (P_a+P_b) }; } // namespace qdk::chemistry::scf diff --git a/cpp/src/qdk/chemistry/algorithms/microsoft/scf/src/scf_algorithm/gdm.cpp b/cpp/src/qdk/chemistry/algorithms/microsoft/scf/src/scf_algorithm/gdm.cpp index 1dd431b8c..6c0885626 100644 --- a/cpp/src/qdk/chemistry/algorithms/microsoft/scf/src/scf_algorithm/gdm.cpp +++ b/cpp/src/qdk/chemistry/algorithms/microsoft/scf/src/scf_algorithm/gdm.cpp @@ -538,21 +538,21 @@ class GDM { private: /** * @brief Transform history matrices (either history_dgrad or history_kappa) - * using current rotation matrices to transform into the - * pseudo-canonical orbital basis, K_new = U_left^T * K_old * U_right + * using current rotation matrices to transform into the + * pseudo-canonical orbital basis, K_new = U_left^T * K_old * U_right * @param[in,out] history History matrix block to be transformed (either * history_dgrad or history_kappa) - * @param[in] u_left Left rotation matrix (e.g., Uii or Uaa) - * @param[in] u_right Right rotation matrix (e.g., Uaa or Uww) + * @param[in] u_left Left rotation matrix (e.g., Uii or Uaa) + * @param[in] u_right Right rotation matrix (e.g., Uaa or Uww) * @param[in] history_size Number of history entries - * @param[in] num_rows Number of rows in each unpacked history block - * @param[in] num_cols Number of cols in each unpacked history block + * @param[in] num_rows Number of rows in each unpacked history block + * @param[in] num_cols Number of cols in each unpacked history block * */ void transform_history_(Eigen::Block& history, const RowMajorMatrix& u_left, const RowMajorMatrix& u_right, const int history_size, - const int num_rows, const int num_cols); + const int num_rows, const int num_cols); /** * @brief Generate pseudo-canonical orbitals and apply transformations @@ -746,16 +746,15 @@ GDM::GDM(const SCFContext& ctx, int history_size_limit) void GDM::transform_history_(Eigen::Block& history, const RowMajorMatrix& u_left, const RowMajorMatrix& u_right, - const int history_size, - const int num_rows, + const int history_size, const int num_rows, const int num_cols) { QDK_LOG_TRACE_ENTERING(); // Validate dimensions (negative values indicate invalid input) if (num_rows < 0 || num_cols < 0) { throw std::invalid_argument( std::string("transform_history_: invalid dimensions (num_rows=") + - std::to_string(num_rows) + ", num_cols=" + - std::to_string(num_cols) + ")"); + std::to_string(num_rows) + ", num_cols=" + std::to_string(num_cols) + + ")"); } // Skip transformation if either dimension is zero if (num_rows == 0 || num_cols == 0) { @@ -913,7 +912,7 @@ static void diagonalize_block_rotate_orbitals( RowMajorMatrix& input_block_output_eigenvectors, const int block_size, Eigen::VectorXd& eigenvalues, const int eigenvalue_start_index, Eigen::Block transformed_orbitals, - const int num_atomic_orbitals, const int num_molecular_orbitals) { + const int num_atomic_orbitals, const int num_molecular_orbitals) { lapack::syev(lapack::Job::Vec, lapack::Uplo::Lower, block_size, input_block_output_eigenvectors.data(), block_size, eigenvalues.data() + eigenvalue_start_index); @@ -927,21 +926,23 @@ static void diagonalize_block_rotate_orbitals( } static void transform_gradient_vector(const Eigen::VectorXd& current_vector, - const int start_index, const int num_rows, - const int num_cols, const RowMajorMatrix& u_left, - const RowMajorMatrix& u_right, - Eigen::VectorXd& transformed_vector) { + const int start_index, const int num_rows, + const int num_cols, + const RowMajorMatrix& u_left, + const RowMajorMatrix& u_right, + Eigen::VectorXd& transformed_vector) { RowMajorMatrix current_matrix = Eigen::Map( current_vector.data() + start_index, num_rows, num_cols); - RowMajorMatrix transformed_matrix = u_left.transpose() * current_matrix * u_right; + RowMajorMatrix transformed_matrix = + u_left.transpose() * current_matrix * u_right; transformed_vector.segment(start_index, num_rows * num_cols) = Eigen::Map(transformed_matrix.data(), num_rows * num_cols); } void GDM::generate_restricted_open_shell_pseudo_canonical_orbital_( - const RowMajorMatrix& F, RowMajorMatrix& C, RowMajorMatrix& history_kappa, - RowMajorMatrix& history_dgrad, Eigen::VectorXd& current_gradient) { + const RowMajorMatrix& F, RowMajorMatrix& C, RowMajorMatrix& history_kappa, + RowMajorMatrix& history_dgrad, Eigen::VectorXd& current_gradient) { const int num_molecular_orbitals = C.cols(); const int num_atomic_orbitals = C.rows(); const int num_closed_orbitals = num_electrons_[1]; @@ -956,29 +957,37 @@ void GDM::generate_restricted_open_shell_pseudo_canonical_orbital_( RowMajorMatrix::Zero(num_molecular_orbitals, num_molecular_orbitals); RowMajorMatrix F_dn_mo = F_up_mo; - F_up_mo.noalias() = C.transpose() * F.block(0, 0, num_atomic_orbitals, num_molecular_orbitals) * C; - F_dn_mo.noalias() = C.transpose() * F.block(num_atomic_orbitals, 0, num_atomic_orbitals, - num_molecular_orbitals) * C; + F_up_mo.noalias() = + C.transpose() * + F.block(0, 0, num_atomic_orbitals, num_molecular_orbitals) * C; + F_dn_mo.noalias() = C.transpose() * + F.block(num_atomic_orbitals, 0, num_atomic_orbitals, + num_molecular_orbitals) * + C; RowMajorMatrix Uii = - F_up_mo.block(0, 0, num_closed_orbitals, num_closed_orbitals) + F_dn_mo.block(0, 0, num_closed_orbitals, num_closed_orbitals); - RowMajorMatrix Uww = F_up_mo.block(num_closed_orbitals, num_closed_orbitals, - num_open_orbitals, num_open_orbitals) + F_dn_mo.block(num_closed_orbitals, num_closed_orbitals, - num_open_orbitals, num_open_orbitals); - RowMajorMatrix Uaa = F_up_mo.block(num_occupied_orbitals, num_occupied_orbitals, - num_virtual_orbitals, num_virtual_orbitals) + F_dn_mo.block(num_occupied_orbitals, num_occupied_orbitals, - num_virtual_orbitals, num_virtual_orbitals); + 0.5 * (F_up_mo.block(0, 0, num_closed_orbitals, num_closed_orbitals) + + F_dn_mo.block(0, 0, num_closed_orbitals, num_closed_orbitals)); + RowMajorMatrix Uww = + 0.5 * (F_up_mo.block(num_closed_orbitals, num_closed_orbitals, + num_open_orbitals, num_open_orbitals) + + F_dn_mo.block(num_closed_orbitals, num_closed_orbitals, + num_open_orbitals, num_open_orbitals)); + RowMajorMatrix Uaa = + 0.5 * (F_up_mo.block(num_occupied_orbitals, num_occupied_orbitals, + num_virtual_orbitals, num_virtual_orbitals) + + F_dn_mo.block(num_occupied_orbitals, num_occupied_orbitals, + num_virtual_orbitals, num_virtual_orbitals)); // Compute eigenvalues/eigenvectors of occupied-occupied and virtual-virtual // blocks for pseudo-canonical orbital transformation // Then transform the orbitals to pseudo-canonical basis - auto C_closed_view = - C.block(0, 0, num_atomic_orbitals, num_closed_orbitals); + auto C_closed_view = C.block(0, 0, num_atomic_orbitals, num_closed_orbitals); diagonalize_block_rotate_orbitals( - Uii, num_closed_orbitals, pseudo_canonical_eigenvalues_, 0, - C_closed_view, num_atomic_orbitals, num_molecular_orbitals); + Uii, num_closed_orbitals, pseudo_canonical_eigenvalues_, 0, C_closed_view, + num_atomic_orbitals, num_molecular_orbitals); - auto C_open_view = C.block(0, num_closed_orbitals, num_atomic_orbitals, - num_open_orbitals); + auto C_open_view = + C.block(0, num_closed_orbitals, num_atomic_orbitals, num_open_orbitals); diagonalize_block_rotate_orbitals( Uww, num_open_orbitals, pseudo_canonical_eigenvalues_, num_closed_orbitals, C_open_view, num_atomic_orbitals, @@ -995,51 +1004,52 @@ void GDM::generate_restricted_open_shell_pseudo_canonical_orbital_( // accommodate current pseudo-canonical orbitals int offset = 0; auto history_kappa_iw = - history_kappa.block(0, 0, history_size_, size_closed_open); + history_kappa.block(0, 0, history_size_, size_closed_open); transform_history_(history_kappa_iw, Uii, Uww, history_size_, - num_closed_orbitals, num_open_orbitals); + num_closed_orbitals, num_open_orbitals); auto history_dgrad_iw = - history_dgrad.block(0, 0, history_size_, size_closed_open); + history_dgrad.block(0, 0, history_size_, size_closed_open); transform_history_(history_dgrad_iw, Uii, Uww, history_size_, - num_closed_orbitals, num_open_orbitals); + num_closed_orbitals, num_open_orbitals); offset += size_closed_open; auto history_kappa_wa = - history_kappa.block(0, offset, history_size_, size_open_virtual); + history_kappa.block(0, offset, history_size_, size_open_virtual); transform_history_(history_kappa_wa, Uww, Uaa, history_size_, - num_open_orbitals, num_virtual_orbitals); + num_open_orbitals, num_virtual_orbitals); auto history_dgrad_wa = - history_dgrad.block(0, offset, history_size_, size_open_virtual); + history_dgrad.block(0, offset, history_size_, size_open_virtual); transform_history_(history_dgrad_wa, Uww, Uaa, history_size_, - num_open_orbitals, num_virtual_orbitals); + num_open_orbitals, num_virtual_orbitals); offset += size_open_virtual; auto history_kappa_ia = - history_kappa.block(0, offset, history_size_, size_closed_virtual); + history_kappa.block(0, offset, history_size_, size_closed_virtual); transform_history_(history_kappa_ia, Uii, Uaa, history_size_, - num_closed_orbitals, num_virtual_orbitals); + num_closed_orbitals, num_virtual_orbitals); auto history_dgrad_ia = - history_dgrad.block(0, offset, history_size_, size_closed_virtual); + history_dgrad.block(0, offset, history_size_, size_closed_virtual); transform_history_(history_dgrad_ia, Uii, Uaa, history_size_, - num_closed_orbitals, num_virtual_orbitals); + num_closed_orbitals, num_virtual_orbitals); // Transform the gradient to accommodate current pseudo-canonical orbitals - Eigen::VectorXd current_gradient_transformed = Eigen::VectorXd(current_gradient.size()); + Eigen::VectorXd current_gradient_transformed = + Eigen::VectorXd(current_gradient.size()); offset = 0; transform_gradient_vector(current_gradient, offset, num_closed_orbitals, - num_open_orbitals, Uii, Uww, - current_gradient_transformed); + num_open_orbitals, Uii, Uww, + current_gradient_transformed); offset += size_closed_open; transform_gradient_vector(current_gradient, offset, num_open_orbitals, - num_virtual_orbitals, Uww, Uaa, - current_gradient_transformed); + num_virtual_orbitals, Uww, Uaa, + current_gradient_transformed); offset += size_open_virtual; transform_gradient_vector(current_gradient, offset, num_closed_orbitals, - num_virtual_orbitals, Uii, Uaa, - current_gradient_transformed); + num_virtual_orbitals, Uii, Uaa, + current_gradient_transformed); current_gradient = current_gradient_transformed; } @@ -1047,9 +1057,61 @@ void GDM::generate_restricted_open_shell_pseudo_canonical_orbital_( void GDM::build_restricted_open_shell_pseudo_canonical_orbitals_hessian_( const RowMajorMatrix& F, RowMajorMatrix& C, int num_molecular_orbitals, Eigen::VectorXd& initial_hessian) { + initial_hessian.setZero(total_rotation_size_); + + generate_restricted_open_shell_pseudo_canonical_orbital_( + F, C, history_kappa_, history_dgrad_, current_gradient_); + const int num_closed_orbitals = num_electrons_[1]; const int num_open_orbitals = num_electrons_[0] - num_closed_orbitals; const int num_virtual_orbitals = num_molecular_orbitals - num_electrons_[0]; + + double initial_hessian_coeff = + 2.0; // between close and open shells, or between open and virtual + // orbitals, the coefficient should be 2.0 + int offset = 0; + for (int j = 0; j < num_closed_orbitals; j++) { + for (int v = 0; v < num_open_orbitals; v++) { + int index = j * num_open_orbitals + v; + double pseudo_canonical_energy_diff = + std::abs(pseudo_canonical_eigenvalues_(num_closed_orbitals + v) - + pseudo_canonical_eigenvalues_(j)); + initial_hessian(index) = + std::max(initial_hessian_coeff * + (std::abs(delta_energy_) + pseudo_canonical_energy_diff), + nonpositive_threshold_); + } + } + offset += num_closed_orbitals * num_open_orbitals; + + for (int v = 0; v < num_open_orbitals; v++) { + for (int a = 0; a < num_virtual_orbitals; a++) { + int index = offset + v * num_virtual_orbitals + a; + double pseudo_canonical_energy_diff = + std::abs(pseudo_canonical_eigenvalues_(num_electrons_[0] + a) - + pseudo_canonical_eigenvalues_(num_closed_orbitals + v)); + initial_hessian(index) = + std::max(initial_hessian_coeff * + (std::abs(delta_energy_) + pseudo_canonical_energy_diff), + nonpositive_threshold_); + } + } + offset += num_open_orbitals * num_virtual_orbitals; + + initial_hessian_coeff = 4.0; // between closed and virtual orbitals, the + // coefficient should be 4.0 + for (int j = 0; j < num_closed_orbitals; j++) { + for (int a = 0; a < num_virtual_orbitals; a++) { + int index = offset + j * num_virtual_orbitals + a; + double pseudo_canonical_energy_diff = + std::abs(pseudo_canonical_eigenvalues_(num_electrons_[0] + a) - + pseudo_canonical_eigenvalues_(j)); + initial_hessian(index) = + std::max(initial_hessian_coeff * + (std::abs(delta_energy_) + pseudo_canonical_energy_diff), + nonpositive_threshold_); + } + } } void GDM::iterate(SCFImpl& scf_impl) { diff --git a/cpp/src/qdk/chemistry/algorithms/microsoft/scf/src/scf_algorithm/scf_algorithm.cpp b/cpp/src/qdk/chemistry/algorithms/microsoft/scf/src/scf_algorithm/scf_algorithm.cpp index 48c620e79..c5d66c111 100644 --- a/cpp/src/qdk/chemistry/algorithms/microsoft/scf/src/scf_algorithm/scf_algorithm.cpp +++ b/cpp/src/qdk/chemistry/algorithms/microsoft/scf/src/scf_algorithm/scf_algorithm.cpp @@ -7,9 +7,11 @@ #include #include +#include #include #include #include +#include #include "../scf/scf_impl.h" #ifdef QDK_CHEMISTRY_ENABLE_GPU @@ -47,6 +49,13 @@ SCFAlgorithm::SCFAlgorithm(const SCFContext& ctx) : 1; P_last_ = RowMajorMatrix::Zero(num_density_matrices * num_atomic_orbitals, num_atomic_orbitals); + + if (ctx.cfg->scf_orbital_type == SCFOrbitalType::RestrictedOpenShell) { + rohf_effective_fock_ = + RowMajorMatrix::Zero(num_atomic_orbitals, num_atomic_orbitals); + rohf_total_density_ = + RowMajorMatrix::Zero(num_atomic_orbitals, num_atomic_orbitals); + } } SCFAlgorithm::~SCFAlgorithm() noexcept = default; @@ -69,15 +78,15 @@ std::shared_ptr SCFAlgorithm::create(const SCFContext& ctx) { return std::make_shared(ctx, cfg.scf_algorithm.diis_subspace_size); case SCFAlgorithmName::GDM: - if (rohf_enabled) { - throw std::runtime_error("ROHF-enabled GDM is not supported!"); - } + // if (rohf_enabled) { + // throw std::runtime_error("ROHF-enabled GDM is not supported!"); + // } return std::make_shared(ctx, cfg.scf_algorithm.gdm_config); case SCFAlgorithmName::DIIS_GDM: - if (rohf_enabled) { - throw std::runtime_error("ROHF-enabled DIIS_GDM is not supported!"); - } + // if (rohf_enabled) { + // throw std::runtime_error("ROHF-enabled DIIS_GDM is not supported!"); + // } return std::make_shared(ctx, cfg.scf_algorithm.diis_subspace_size, cfg.scf_algorithm.gdm_config); @@ -194,6 +203,170 @@ void SCFAlgorithm::update_density_matrix(RowMajorMatrix& P, } } +bool SCFAlgorithm::try_get_rohf_convergence_matrices( + const SCFImpl& scf_impl, const RowMajorMatrix*& fock_matrix, + const RowMajorMatrix*& density_matrix) { + QDK_LOG_TRACE_ENTERING(); + if (!ensure_rohf_convergence_matrices_(scf_impl)) { + return false; + } + fock_matrix = &get_rohf_convergence_fock_matrix(); + density_matrix = &get_rohf_convergence_density_matrix(); + return true; +} + +void SCFAlgorithm::build_rohf_f_p_matrix(const RowMajorMatrix& F, + const RowMajorMatrix& C, + const RowMajorMatrix& P, + int nelec_alpha, int nelec_beta, + RowMajorMatrix& effective_fock, + RowMajorMatrix& total_density) { + QDK_LOG_TRACE_ENTERING(); + const int num_atomic_orbitals = static_cast(C.rows()); + const int num_molecular_orbitals = static_cast(C.cols()); + if (num_atomic_orbitals != num_molecular_orbitals) { + throw std::runtime_error( + "ROHF build requires number of atomic orbitals to equal number of " + "molecular orbitals!"); + } + + total_density = + P.block(0, 0, num_atomic_orbitals, num_atomic_orbitals) + + P.block(num_atomic_orbitals, 0, num_atomic_orbitals, num_atomic_orbitals); + + if (effective_fock.rows() != num_atomic_orbitals || + effective_fock.cols() != num_atomic_orbitals) { + effective_fock = + RowMajorMatrix::Zero(num_atomic_orbitals, num_atomic_orbitals); + } + + if (C.isZero()) { + effective_fock.noalias() = + F.block(0, 0, num_atomic_orbitals, num_atomic_orbitals); + return; + } + + RowMajorMatrix F_up_mo = + RowMajorMatrix::Zero(num_molecular_orbitals, num_molecular_orbitals); + RowMajorMatrix F_dn_mo = F_up_mo; + RowMajorMatrix effective_F_mo = F_up_mo; + + F_up_mo.noalias() = C.transpose() * + F.block(0, 0, num_atomic_orbitals, num_atomic_orbitals) * + C; + F_dn_mo.noalias() = C.transpose() * + F.block(num_atomic_orbitals, 0, num_atomic_orbitals, + num_atomic_orbitals) * + C; + + auto average_block = [&effective_F_mo, &F_up_mo, &F_dn_mo]( + int row, int col, int rows, int cols) { + if (rows <= 0 || cols <= 0) return; + effective_F_mo.block(row, col, rows, cols).noalias() = + 0.5 * (F_up_mo.block(row, col, rows, cols) + + F_dn_mo.block(row, col, rows, cols)); + }; + auto copy_block = [&effective_F_mo](const RowMajorMatrix& src, int row, + int col, int rows, int cols) { + if (rows <= 0 || cols <= 0) return; + effective_F_mo.block(row, col, rows, cols) = + src.block(row, col, rows, cols); + }; + + const int nd = nelec_beta; + const int ns = nelec_alpha - nelec_beta; + const int nv = num_molecular_orbitals - nelec_alpha; + + average_block(0, 0, nd, nd); + average_block(0, nd + ns, nd, nv); + average_block(nd + ns, 0, nv, nd); + average_block(nd + ns, nd + ns, nv, nv); + average_block(nd, nd, ns, ns); + copy_block(F_dn_mo, 0, nd, nd, ns); + copy_block(F_dn_mo, nd, 0, ns, nd); + copy_block(F_up_mo, nd, nd + ns, ns, nv); + copy_block(F_up_mo, nd + ns, nd, nv, ns); + + // Transform the effective Fock matrix back to AO basis by solving + // C^{-T} * F_MO * C^{-1} = F_AO + // We use LAPACK's getrf/getrs to solve the linear systems involving C^T and + // C without explicitly inverting C + const int matrix_dim = num_molecular_orbitals; + using ColMajorMatrix = + Eigen::Matrix; + // LAPACK expects column-major layout, so we copy the row-major data into a + // column-major matrix without transposing the logical layout + ColMajorMatrix Ct = + Eigen::Map(C.data(), matrix_dim, C.rows()); + // F_MO is symmetric, so we can use it directly as the right-hand side + // without transposing + ColMajorMatrix temp_rhs = effective_F_mo; + std::vector ipiv(matrix_dim); + + auto info = + lapack::getrf(matrix_dim, matrix_dim, Ct.data(), matrix_dim, ipiv.data()); + if (info != 0) { + throw std::runtime_error("getrf failed while factorizing C^T"); + } + + info = lapack::getrs(lapack::Op::NoTrans, matrix_dim, matrix_dim, Ct.data(), + matrix_dim, ipiv.data(), temp_rhs.data(), matrix_dim); + if (info != 0) { + throw std::runtime_error("getrs failed while solving C^T X = F_mo"); + } + + temp_rhs.transposeInPlace(); + info = lapack::getrs(lapack::Op::NoTrans, matrix_dim, matrix_dim, Ct.data(), + matrix_dim, ipiv.data(), temp_rhs.data(), matrix_dim); + if (info != 0) { + throw std::runtime_error("getrs failed while solving C^T X = M^T"); + } + + effective_fock = temp_rhs.transpose(); + if (!effective_fock.isApprox(effective_fock.transpose())) { + effective_fock = 0.5 * (effective_fock + effective_fock.transpose().eval()); + } +} + +bool SCFAlgorithm::ensure_rohf_convergence_matrices_(const SCFImpl& scf_impl) { + QDK_LOG_TRACE_ENTERING(); + if (ctx_.cfg->scf_orbital_type != SCFOrbitalType::RestrictedOpenShell) { + return false; + } + + const auto nelec_vec = scf_impl.get_num_electrons(); + build_rohf_f_p_matrix( + scf_impl.get_fock_matrix(), scf_impl.get_orbitals_matrix(), + scf_impl.get_density_matrix(), nelec_vec[0], nelec_vec[1], + rohf_effective_fock_, rohf_total_density_); + return true; +} + +const RowMajorMatrix& SCFAlgorithm::get_rohf_convergence_fock_matrix() const { + QDK_LOG_TRACE_ENTERING(); + if (rohf_effective_fock_.size() == 0) { + throw std::logic_error("ROHF convergence cache not initialized"); + } + return rohf_effective_fock_; +} + +const RowMajorMatrix& SCFAlgorithm::get_rohf_convergence_density_matrix() + const { + QDK_LOG_TRACE_ENTERING(); + if (rohf_total_density_.size() == 0) { + throw std::logic_error("ROHF convergence cache not initialized"); + } + return rohf_total_density_; +} + +RowMajorMatrix& SCFAlgorithm::rohf_convergence_density_matrix() { + QDK_LOG_TRACE_ENTERING(); + if (rohf_total_density_.size() == 0) { + throw std::logic_error("ROHF convergence cache not initialized"); + } + return rohf_total_density_; +} + double SCFAlgorithm::calculate_og_error_(const RowMajorMatrix& F, const RowMajorMatrix& P, const RowMajorMatrix& S, @@ -252,22 +425,12 @@ bool SCFAlgorithm::check_convergence(const SCFImpl& scf_impl) { const RowMajorMatrix* F_ptr; const RowMajorMatrix* P_ptr; - std::vector nelec_vec = scf_impl.get_num_electrons(); - const int nelec[2] = {nelec_vec[0], nelec_vec[1]}; if (ctx_.cfg->scf_orbital_type == SCFOrbitalType::RestrictedOpenShell) { - // To be modified when ROHFGDM is implemented: in that case, the pointer - // will come from the DIIS instance saved in DIIS_GDM, like the current - // DIIS_GDM implementation - DIIS* rohf_diis = dynamic_cast(this); - if (rohf_diis == nullptr) { - throw std::logic_error("ROHF convergence requires DIIS implementation"); + if (!try_get_rohf_convergence_matrices(scf_impl, F_ptr, P_ptr)) { + throw std::logic_error( + "ROHF convergence matrices are not provided by this SCF algorithm"); } - rohf_diis->build_rohf_f_p_matrix( - scf_impl.get_fock_matrix(), scf_impl.get_orbitals_matrix(), - scf_impl.get_density_matrix(), nelec[0], nelec[1]); - F_ptr = &rohf_diis->get_rohf_fock_matrix(); - P_ptr = &rohf_diis->get_rohf_density_matrix(); } else { F_ptr = &scf_impl.get_fock_matrix(); P_ptr = &scf_impl.get_density_matrix(); From 5ddf2d705da4290b17754926cb795d1f3396d8f8 Mon Sep 17 00:00:00 2001 From: v-boqinzhang Date: Tue, 21 Apr 2026 14:34:50 -0700 Subject: [PATCH 07/29] delete ROHF_invalid_GDM test --- cpp/tests/test_scf.cpp | 10 ---------- 1 file changed, 10 deletions(-) diff --git a/cpp/tests/test_scf.cpp b/cpp/tests/test_scf.cpp index 56a37a3b0..3c762b819 100644 --- a/cpp/tests/test_scf.cpp +++ b/cpp/tests/test_scf.cpp @@ -150,16 +150,6 @@ TEST_F(ScfTest, OH_ROHF_INCORE_DIIS) { EXPECT_TRUE(wfn_doublet->get_orbitals()->is_restricted()); } -TEST_F(ScfTest, OH_ROHF_Invalid_GDM) { - auto oh = testing::create_oh_structure(); - auto scf_solver = ScfSolverFactory::create(); - scf_solver->settings().set("enable_gdm", true); - scf_solver->settings().set("method", "hf"); - scf_solver->settings().set("scf_type", "restricted"); - - EXPECT_THROW(scf_solver->run(oh, 0, 2, "sto-3g"), std::runtime_error); -} - TEST_F(ScfTest, Oxygen_atom_gdm) { auto oxygen = testing::create_oxygen_structure(); auto scf_solver = ScfSolverFactory::create(); From f2b15c14f84fe68fe5684a2ec835cc691a81ed9f Mon Sep 17 00:00:00 2001 From: v-boqinzhang Date: Wed, 22 Apr 2026 10:24:52 -0700 Subject: [PATCH 08/29] add tests --- .../scf/src/scf_algorithm/scf_algorithm.cpp | 6 ----- cpp/tests/test_scf.cpp | 26 +++++++++++++++++++ 2 files changed, 26 insertions(+), 6 deletions(-) diff --git a/cpp/src/qdk/chemistry/algorithms/microsoft/scf/src/scf_algorithm/scf_algorithm.cpp b/cpp/src/qdk/chemistry/algorithms/microsoft/scf/src/scf_algorithm/scf_algorithm.cpp index c5d66c111..06d9bf1a3 100644 --- a/cpp/src/qdk/chemistry/algorithms/microsoft/scf/src/scf_algorithm/scf_algorithm.cpp +++ b/cpp/src/qdk/chemistry/algorithms/microsoft/scf/src/scf_algorithm/scf_algorithm.cpp @@ -78,15 +78,9 @@ std::shared_ptr SCFAlgorithm::create(const SCFContext& ctx) { return std::make_shared(ctx, cfg.scf_algorithm.diis_subspace_size); case SCFAlgorithmName::GDM: - // if (rohf_enabled) { - // throw std::runtime_error("ROHF-enabled GDM is not supported!"); - // } return std::make_shared(ctx, cfg.scf_algorithm.gdm_config); case SCFAlgorithmName::DIIS_GDM: - // if (rohf_enabled) { - // throw std::runtime_error("ROHF-enabled DIIS_GDM is not supported!"); - // } return std::make_shared(ctx, cfg.scf_algorithm.diis_subspace_size, cfg.scf_algorithm.gdm_config); diff --git a/cpp/tests/test_scf.cpp b/cpp/tests/test_scf.cpp index 3c762b819..ed0a42beb 100644 --- a/cpp/tests/test_scf.cpp +++ b/cpp/tests/test_scf.cpp @@ -150,6 +150,18 @@ TEST_F(ScfTest, OH_ROHF_INCORE_DIIS) { EXPECT_TRUE(wfn_doublet->get_orbitals()->is_restricted()); } +TEST_F(ScfTest, OH_ROKS_invalid) { + auto oh = testing::create_oh_structure(); + auto scf_solver = ScfSolverFactory::create(); + scf_solver->settings().set("enable_gdm", true); + scf_solver->settings().set("method", "pbe"); + scf_solver->settings().set("scf_type", "restricted"); + + // Default should be a singlet + EXPECT_THROW(scf_solver->run(oh, 0, 2, "sto-3g"), + std::invalid_argument); // open-shell dublet +} + TEST_F(ScfTest, Oxygen_atom_gdm) { auto oxygen = testing::create_oxygen_structure(); auto scf_solver = ScfSolverFactory::create(); @@ -217,6 +229,20 @@ TEST_F(ScfTest, Oxygen_atom_charged_doublet_gdm) { EXPECT_FALSE(wfn_doublet->get_orbitals()->is_restricted()); } +TEST_F(ScfTest, OH_ROHF_GDM) { + auto oh = testing::create_oh_structure(); + auto scf_solver = ScfSolverFactory::create(); + scf_solver->settings().set("enable_gdm", true); + scf_solver->settings().set("method", "hf"); + scf_solver->settings().set("scf_type", "restricted"); + auto [E_doublet, wfn_doublet] = scf_solver->run(oh, 0, 2, "sto-3g"); + + EXPECT_NEAR(E_doublet, -74.361530753176, testing::scf_energy_tolerance); + + // Check doublet orbitals + EXPECT_TRUE(wfn_doublet->get_orbitals()->is_restricted()); +} + TEST_F(ScfTest, Oxygen_atom_invalid_energy_thresh_diis_switch_gdm) { auto oxygen = testing::create_oxygen_structure(); auto scf_solver = ScfSolverFactory::create(); From afb3617e2e5cc20871f2db7d87021257fcd47d5b Mon Sep 17 00:00:00 2001 From: v-boqinzhang Date: Wed, 22 Apr 2026 11:07:05 -0700 Subject: [PATCH 09/29] clean the code 1 --- .../qdk/chemistry/scf/core/scf_algorithm.h | 18 ++++++++----- .../microsoft/scf/src/scf_algorithm/diis.cpp | 26 +++++++------------ .../microsoft/scf/src/scf_algorithm/diis.h | 10 +++---- .../scf/src/scf_algorithm/scf_algorithm.cpp | 23 ++++++---------- 4 files changed, 32 insertions(+), 45 deletions(-) diff --git a/cpp/src/qdk/chemistry/algorithms/microsoft/scf/include/qdk/chemistry/scf/core/scf_algorithm.h b/cpp/src/qdk/chemistry/algorithms/microsoft/scf/include/qdk/chemistry/scf/core/scf_algorithm.h index 120333bb8..d87651fa6 100644 --- a/cpp/src/qdk/chemistry/algorithms/microsoft/scf/include/qdk/chemistry/scf/core/scf_algorithm.h +++ b/cpp/src/qdk/chemistry/algorithms/microsoft/scf/include/qdk/chemistry/scf/core/scf_algorithm.h @@ -159,6 +159,17 @@ class SCFAlgorithm { * * Converts spin-blocked Fock/density matrices into the effective ROHF Fock * and total-density representation used for OG evaluation. + * + * @param[in] F Spin-blocked Fock matrix in AO basis with alpha and beta + * blocks stacked by row + * @param[in] C Molecular-orbital coefficient matrix used for AO<->MO + * transformations + * @param[in] P Spin-blocked density matrix in AO basis with alpha and beta + * blocks stacked by row + * @param[in] nelec_alpha Number of alpha electrons + * @param[in] nelec_beta Number of beta electrons + * @param[out] effective_fock Effective ROHF Fock matrix in AO basis + * @param[out] total_density Total AO density matrix (P_alpha + P_beta) */ static void build_rohf_f_p_matrix(const RowMajorMatrix& F, const RowMajorMatrix& C, @@ -167,13 +178,6 @@ class SCFAlgorithm { RowMajorMatrix& effective_fock, RowMajorMatrix& total_density); - /** - * @brief Refresh cached ROHF convergence matrices from current SCF state - * - * @return true if ROHF matrices were refreshed, false for non-ROHF runs - */ - bool ensure_rohf_convergence_matrices_(const SCFImpl& scf_impl); - /** * @brief Access cached ROHF effective Fock matrix */ diff --git a/cpp/src/qdk/chemistry/algorithms/microsoft/scf/src/scf_algorithm/diis.cpp b/cpp/src/qdk/chemistry/algorithms/microsoft/scf/src/scf_algorithm/diis.cpp index 8d3555084..cbd575492 100644 --- a/cpp/src/qdk/chemistry/algorithms/microsoft/scf/src/scf_algorithm/diis.cpp +++ b/cpp/src/qdk/chemistry/algorithms/microsoft/scf/src/scf_algorithm/diis.cpp @@ -310,8 +310,7 @@ void DIIS::iterate(SCFImpl& scf_impl) { // Decide which Fock/density view to feed into Pulay: RHF/UHF use the direct // SCFImpl matrices, while ROHF works on the cached total-density view. - RowMajorMatrix& working_density = select_working_density(scf_impl); - const RowMajorMatrix& working_fock = select_working_fock(scf_impl); + auto [working_density, working_fock] = select_working_matrices(scf_impl); auto& C = scf_impl.orbitals_matrix(); const auto& S = scf_impl.overlap(); @@ -370,26 +369,19 @@ double DIIS::current_diis_error() const { return diis_impl_->get_diis_error(); } -RowMajorMatrix& DIIS::select_working_density(SCFImpl& scf_impl) { +std::pair DIIS::select_working_matrices( + SCFImpl& scf_impl) { QDK_LOG_TRACE_ENTERING(); if (ctx_.cfg->scf_orbital_type == SCFOrbitalType::RestrictedOpenShell) { - if (!ensure_rohf_convergence_matrices_(scf_impl)) { + const RowMajorMatrix* rohf_fock = nullptr; + const RowMajorMatrix* rohf_density = nullptr; + if (!try_get_rohf_convergence_matrices(scf_impl, rohf_fock, + rohf_density)) { throw std::logic_error("ROHF convergence matrices are unavailable"); } - return rohf_convergence_density_matrix(); + return {rohf_convergence_density_matrix(), *rohf_fock}; } - return scf_impl.density_matrix(); -} - -const RowMajorMatrix& DIIS::select_working_fock(SCFImpl& scf_impl) { - QDK_LOG_TRACE_ENTERING(); - if (ctx_.cfg->scf_orbital_type == SCFOrbitalType::RestrictedOpenShell) { - if (!ensure_rohf_convergence_matrices_(scf_impl)) { - throw std::logic_error("ROHF convergence matrices are unavailable"); - } - return get_rohf_convergence_fock_matrix(); - } - return scf_impl.get_fock_matrix(); + return {scf_impl.density_matrix(), scf_impl.get_fock_matrix()}; } } // namespace qdk::chemistry::scf diff --git a/cpp/src/qdk/chemistry/algorithms/microsoft/scf/src/scf_algorithm/diis.h b/cpp/src/qdk/chemistry/algorithms/microsoft/scf/src/scf_algorithm/diis.h index 0dc792eea..261d56c56 100644 --- a/cpp/src/qdk/chemistry/algorithms/microsoft/scf/src/scf_algorithm/diis.h +++ b/cpp/src/qdk/chemistry/algorithms/microsoft/scf/src/scf_algorithm/diis.h @@ -8,6 +8,7 @@ #include #include +#include namespace qdk::chemistry::scf { @@ -76,13 +77,10 @@ class DIIS : public SCFAlgorithm { double current_diis_error() const; /** - * @brief Select the density view used for the current iteration + * @brief Select the density/fock views used for the current iteration */ - RowMajorMatrix& select_working_density(SCFImpl& scf_impl); - /** - * @brief Select the Fock matrix view used for the current iteration - */ - const RowMajorMatrix& select_working_fock(SCFImpl& scf_impl); + std::pair select_working_matrices( + SCFImpl& scf_impl); std::unique_ptr diis_impl_; ///< Pulay DIIS core implementation }; diff --git a/cpp/src/qdk/chemistry/algorithms/microsoft/scf/src/scf_algorithm/scf_algorithm.cpp b/cpp/src/qdk/chemistry/algorithms/microsoft/scf/src/scf_algorithm/scf_algorithm.cpp index 06d9bf1a3..eaadc5cc6 100644 --- a/cpp/src/qdk/chemistry/algorithms/microsoft/scf/src/scf_algorithm/scf_algorithm.cpp +++ b/cpp/src/qdk/chemistry/algorithms/microsoft/scf/src/scf_algorithm/scf_algorithm.cpp @@ -201,9 +201,16 @@ bool SCFAlgorithm::try_get_rohf_convergence_matrices( const SCFImpl& scf_impl, const RowMajorMatrix*& fock_matrix, const RowMajorMatrix*& density_matrix) { QDK_LOG_TRACE_ENTERING(); - if (!ensure_rohf_convergence_matrices_(scf_impl)) { + if (ctx_.cfg->scf_orbital_type != SCFOrbitalType::RestrictedOpenShell) { return false; } + + const auto nelec_vec = scf_impl.get_num_electrons(); + build_rohf_f_p_matrix( + scf_impl.get_fock_matrix(), scf_impl.get_orbitals_matrix(), + scf_impl.get_density_matrix(), nelec_vec[0], nelec_vec[1], + rohf_effective_fock_, rohf_total_density_); + fock_matrix = &get_rohf_convergence_fock_matrix(); density_matrix = &get_rohf_convergence_density_matrix(); return true; @@ -322,20 +329,6 @@ void SCFAlgorithm::build_rohf_f_p_matrix(const RowMajorMatrix& F, } } -bool SCFAlgorithm::ensure_rohf_convergence_matrices_(const SCFImpl& scf_impl) { - QDK_LOG_TRACE_ENTERING(); - if (ctx_.cfg->scf_orbital_type != SCFOrbitalType::RestrictedOpenShell) { - return false; - } - - const auto nelec_vec = scf_impl.get_num_electrons(); - build_rohf_f_p_matrix( - scf_impl.get_fock_matrix(), scf_impl.get_orbitals_matrix(), - scf_impl.get_density_matrix(), nelec_vec[0], nelec_vec[1], - rohf_effective_fock_, rohf_total_density_); - return true; -} - const RowMajorMatrix& SCFAlgorithm::get_rohf_convergence_fock_matrix() const { QDK_LOG_TRACE_ENTERING(); if (rohf_effective_fock_.size() == 0) { From 09e9a9ff861a718dbca984446d5275e3471be80d Mon Sep 17 00:00:00 2001 From: v-boqinzhang Date: Wed, 22 Apr 2026 14:54:40 -0700 Subject: [PATCH 10/29] clean the code 2 --- .../qdk/chemistry/scf/core/scf_algorithm.h | 22 +- .../microsoft/scf/src/scf/scf_impl.h | 16 +- .../microsoft/scf/src/scf_algorithm/diis.cpp | 5 +- .../microsoft/scf/src/scf_algorithm/gdm.cpp | 260 ++++++++++-------- 4 files changed, 175 insertions(+), 128 deletions(-) diff --git a/cpp/src/qdk/chemistry/algorithms/microsoft/scf/include/qdk/chemistry/scf/core/scf_algorithm.h b/cpp/src/qdk/chemistry/algorithms/microsoft/scf/include/qdk/chemistry/scf/core/scf_algorithm.h index d87651fa6..774bb4ddc 100644 --- a/cpp/src/qdk/chemistry/algorithms/microsoft/scf/include/qdk/chemistry/scf/core/scf_algorithm.h +++ b/cpp/src/qdk/chemistry/algorithms/microsoft/scf/include/qdk/chemistry/scf/core/scf_algorithm.h @@ -159,17 +159,17 @@ class SCFAlgorithm { * * Converts spin-blocked Fock/density matrices into the effective ROHF Fock * and total-density representation used for OG evaluation. - * - * @param[in] F Spin-blocked Fock matrix in AO basis with alpha and beta - * blocks stacked by row - * @param[in] C Molecular-orbital coefficient matrix used for AO<->MO - * transformations - * @param[in] P Spin-blocked density matrix in AO basis with alpha and beta - * blocks stacked by row - * @param[in] nelec_alpha Number of alpha electrons - * @param[in] nelec_beta Number of beta electrons - * @param[out] effective_fock Effective ROHF Fock matrix in AO basis - * @param[out] total_density Total AO density matrix (P_alpha + P_beta) + * + * @param[in] F Spin-blocked Fock matrix in AO basis with alpha and beta + * blocks stacked by row + * @param[in] C Molecular-orbital coefficient matrix used for AO<->MO + * transformations + * @param[in] P Spin-blocked density matrix in AO basis with alpha and beta + * blocks stacked by row + * @param[in] nelec_alpha Number of alpha electrons + * @param[in] nelec_beta Number of beta electrons + * @param[out] effective_fock Effective ROHF Fock matrix in AO basis + * @param[out] total_density Total AO density matrix (P_alpha + P_beta) */ static void build_rohf_f_p_matrix(const RowMajorMatrix& F, const RowMajorMatrix& C, diff --git a/cpp/src/qdk/chemistry/algorithms/microsoft/scf/src/scf/scf_impl.h b/cpp/src/qdk/chemistry/algorithms/microsoft/scf/src/scf/scf_impl.h index e0c27a66e..6c6fb6650 100644 --- a/cpp/src/qdk/chemistry/algorithms/microsoft/scf/src/scf/scf_impl.h +++ b/cpp/src/qdk/chemistry/algorithms/microsoft/scf/src/scf/scf_impl.h @@ -155,10 +155,18 @@ class SCFImpl { const RowMajorMatrix& get_core_hamiltonian() const { return H_; } /** - * @brief Build Coulomb (J) and exchange (K) matrices for a density matrix - * @param density_matrix Density matrix (num_density_matrices x NAO x NAO) - * @param J Output Coulomb matrix (same size as density_matrix) - * @param K Output exchange matrix (same size as density_matrix) + * @brief Build Coulomb (J) and exchange (K) matrices for a given density + * + * Utility JK build for a caller-supplied density matrix. + * Currently used by the GDM trial-rotation gradient path + * (GDMLineFunctor / ROHF generalized-gradient evaluation). + * It is not part of the main SCFImpl::iterate_ update loop. + * Keeps this trial-path two-electron build logic in one place. + * + * @param[in] density_matrix Density matrix + * (num_density_matrices x NAO x NAO) + * @param[out] J Output Coulomb matrix (same size as density_matrix) + * @param[out] K Output exchange matrix (same size as density_matrix) */ void build_jk_matrices(const RowMajorMatrix& density_matrix, RowMajorMatrix& J, RowMajorMatrix& K) const; diff --git a/cpp/src/qdk/chemistry/algorithms/microsoft/scf/src/scf_algorithm/diis.cpp b/cpp/src/qdk/chemistry/algorithms/microsoft/scf/src/scf_algorithm/diis.cpp index cbd575492..4a35dfb3d 100644 --- a/cpp/src/qdk/chemistry/algorithms/microsoft/scf/src/scf_algorithm/diis.cpp +++ b/cpp/src/qdk/chemistry/algorithms/microsoft/scf/src/scf_algorithm/diis.cpp @@ -370,13 +370,12 @@ double DIIS::current_diis_error() const { } std::pair DIIS::select_working_matrices( - SCFImpl& scf_impl) { + SCFImpl& scf_impl) { QDK_LOG_TRACE_ENTERING(); if (ctx_.cfg->scf_orbital_type == SCFOrbitalType::RestrictedOpenShell) { const RowMajorMatrix* rohf_fock = nullptr; const RowMajorMatrix* rohf_density = nullptr; - if (!try_get_rohf_convergence_matrices(scf_impl, rohf_fock, - rohf_density)) { + if (!try_get_rohf_convergence_matrices(scf_impl, rohf_fock, rohf_density)) { throw std::logic_error("ROHF convergence matrices are unavailable"); } return {rohf_convergence_density_matrix(), *rohf_fock}; diff --git a/cpp/src/qdk/chemistry/algorithms/microsoft/scf/src/scf_algorithm/gdm.cpp b/cpp/src/qdk/chemistry/algorithms/microsoft/scf/src/scf_algorithm/gdm.cpp index 6c0885626..e255b99aa 100644 --- a/cpp/src/qdk/chemistry/algorithms/microsoft/scf/src/scf_algorithm/gdm.cpp +++ b/cpp/src/qdk/chemistry/algorithms/microsoft/scf/src/scf_algorithm/gdm.cpp @@ -263,6 +263,7 @@ static void compute_restricted_open_shell_gradient( const auto C_ao_mo = C.block(0, 0, num_atomic_orbitals, num_molecular_orbitals); + // Calculate Generalized Fock matrix in MO basis RowMajorMatrix H_mo = C_ao_mo.transpose() * H_ao * C_ao_mo; RowMajorMatrix J_alpha_mo = C_ao_mo.transpose() * J_alpha_ao * C_ao_mo; RowMajorMatrix J_beta_mo = C_ao_mo.transpose() * J_beta_ao * C_ao_mo; @@ -573,6 +574,19 @@ class GDM { Eigen::Block history_dgrad_spin, Eigen::VectorBlock current_gradient_spin); + /** + * @brief Build ROHF pseudo-canonical orbitals and transform ROHF history + * + * Diagonalizes closed/open/virtual MO blocks, rotates orbital coefficients + * to the pseudo-canonical basis, and transforms packed ROHF BFGS history and + * current gradient blocks (iw, wa, ia) into the updated basis. + * + * @param[in] F Spin-blocked Fock matrix in AO basis + * @param[in,out] C ROHF orbital coefficient matrix to rotate + * @param[in,out] history_kappa Packed ROHF kappa-history matrix + * @param[in,out] history_dgrad Packed ROHF gradient-difference history + * @param[in,out] current_gradient Packed ROHF current gradient vector + */ void generate_restricted_open_shell_pseudo_canonical_orbital_( const RowMajorMatrix& F, RowMajorMatrix& C, RowMajorMatrix& history_kappa, RowMajorMatrix& history_dgrad, Eigen::VectorXd& current_gradient); @@ -588,6 +602,18 @@ class GDM { const RowMajorMatrix& F, RowMajorMatrix& C, int num_molecular_orbitals, Eigen::VectorXd& initial_hessian); + /** + * @brief Build ROHF initial Hessian in pseudo-canonical basis + * + * Calls ROHF pseudo-canonical transformation, then fills the diagonal + * preconditioner for packed ROHF rotations (iw, wa, ia) using orbital-energy + * differences and the current absolute energy change. + * + * @param[in] F Spin-blocked Fock matrix in AO basis + * @param[in,out] C ROHF orbital coefficient matrix (updated in place) + * @param[in] num_molecular_orbitals Total number of molecular orbitals + * @param[out] initial_hessian Output ROHF initial Hessian vector + */ void build_restricted_open_shell_pseudo_canonical_orbitals_hessian_( const RowMajorMatrix& F, RowMajorMatrix& C, int num_molecular_orbitals, Eigen::VectorXd& initial_hessian); @@ -635,7 +661,10 @@ class GDM { /// vector in this step double last_accepted_energy_; int gdm_step_count_; // GDM iteration counter + + // Number of spin blocks (1 for restricted, 2 for unrestricted) const int num_orbital_spin_blocks_; + // Number of density matrices (1 for closed shell, 2 for open shell) const int num_density_matrices_; }; @@ -773,6 +802,78 @@ void GDM::transform_history_(Eigen::Block& history, } } +/** + * @brief Diagonalize an MO sub-block and rotate the corresponding orbital block + * + * This helper computes eigenpairs of a symmetric MO-space sub-block using + * LAPACK `syev`, stores the resulting eigenvalues into a target slice, and + * applies the eigenvector rotation to the corresponding orbital columns: + * C_block <- C_block * U. + * + * @param[in,out] input_block_output_eigenvectors On input: symmetric block to + * diagonalize. On output: eigenvectors used for rotation. + * @param[in] block_size Dimension of the square sub-block to diagonalize. + * @param[in,out] eigenvalues Vector receiving eigenvalues for this block. + * @param[in] eigenvalue_start_index Start index in `eigenvalues` where this + * block's eigenvalues are written. + * @param[in,out] transformed_orbitals Orbital coefficient sub-block rotated in + * place. + * @param[in] num_atomic_orbitals Row count used for GEMM in the orbital + * rotation. + * @param[in] num_molecular_orbitals Leading dimension of the parent orbital + * coefficient matrix. + */ +static void calculate_pseudo_canonical_orbital_block( + RowMajorMatrix& input_block_output_eigenvectors, const int block_size, + Eigen::VectorXd& eigenvalues, const int eigenvalue_start_index, + Eigen::Block transformed_orbitals, + const int num_atomic_orbitals, const int num_molecular_orbitals) { + lapack::syev(lapack::Job::Vec, lapack::Uplo::Lower, block_size, + input_block_output_eigenvectors.data(), block_size, + eigenvalues.data() + eigenvalue_start_index); + input_block_output_eigenvectors.transposeInPlace(); + RowMajorMatrix copied_orbitals = transformed_orbitals; + blas::gemm(blas::Layout::RowMajor, blas::Op::NoTrans, blas::Op::NoTrans, + num_atomic_orbitals, block_size, block_size, 1.0, + copied_orbitals.data(), block_size, + input_block_output_eigenvectors.data(), block_size, 0.0, + transformed_orbitals.data(), num_molecular_orbitals); +} + +/** + * @brief Rotate a packed gradient sub-block to the pseudo-canonical basis + * + * Interprets a contiguous slice of current_vector (starting at start_index) + * as a num_rows x num_cols matrix G, applies G' = u_left^T * G * u_right, + * and writes the flattened result back into the corresponding slice of + * transformed_vector. + * + * This keeps the gradient representation aligned with the updated + * pseudo-canonical orbitals after block diagonalization and orbital rotation. + * + * @param[in] current_vector Source packed gradient vector. + * @param[in] start_index Start index of the sub-block inside the packed vector. + * @param[in] num_rows Row count of the unpacked gradient block. + * @param[in] num_cols Column count of the unpacked gradient block. + * @param[in] u_left Left rotation matrix. + * @param[in] u_right Right rotation matrix. + * @param[out] transformed_vector Destination packed vector receiving the + * transformed sub-block. + */ +static void rotate_gradient_to_pseudo_canonical_basis( + const Eigen::Ref& current_vector, + const int start_index, const int num_rows, const int num_cols, + const RowMajorMatrix& u_left, const RowMajorMatrix& u_right, + Eigen::Ref transformed_vector) { + RowMajorMatrix current_matrix = Eigen::Map( + current_vector.data() + start_index, num_rows, num_cols); + RowMajorMatrix transformed_matrix = + u_left.transpose() * current_matrix * u_right; + transformed_vector.segment(start_index, num_rows * num_cols) = + Eigen::Map(transformed_matrix.data(), + num_rows * num_cols); +} + void GDM::generate_restricted_unrestricted_pseudo_canonical_orbital_( const RowMajorMatrix& F, RowMajorMatrix& C, const int spin_index, Eigen::Block history_kappa_spin, @@ -794,7 +895,6 @@ void GDM::generate_restricted_unrestricted_pseudo_canonical_orbital_( if (num_occupied_orbitals == 0 || num_virtual_orbitals == 0) { return; } - const int rotation_size = num_occupied_orbitals * num_virtual_orbitals; RowMajorMatrix F_MO = C.block(num_molecular_orbitals * spin_index, 0, num_molecular_orbitals, @@ -805,64 +905,37 @@ void GDM::generate_restricted_unrestricted_pseudo_canonical_orbital_( C.block(num_molecular_orbitals * spin_index, 0, num_molecular_orbitals, num_molecular_orbitals); - // Perform pseudo-canonical transformation and BFGS - // Obtain pseudo-canonical orbitals. Foo and Fvv are symmetric matrices, but - // the output eigenvectors are column-major + // Perform pseudo-canonical transformation + // Diagonalize occupied/virtual blocks and rotate orbitals to the + // pseudo-canonical basis. RowMajorMatrix Uii = F_MO.block(0, 0, num_occupied_orbitals, num_occupied_orbitals); - RowMajorMatrix Uaa = F_MO.block(num_occupied_orbitals, num_occupied_orbitals, - num_virtual_orbitals, num_virtual_orbitals); - - // Compute eigenvalues/eigenvectors of occupied-occupied and virtual-virtual - // blocks for pseudo-canonical orbital transformation - lapack::syev(lapack::Job::Vec, lapack::Uplo::Lower, num_occupied_orbitals, - Uii.data(), num_occupied_orbitals, - pseudo_canonical_eigenvalues_.data()); - lapack::syev(lapack::Job::Vec, lapack::Uplo::Lower, num_virtual_orbitals, - Uaa.data(), num_virtual_orbitals, - pseudo_canonical_eigenvalues_.data() + num_occupied_orbitals); - - // Transpose to convert column-major eigenvectors to row-major format - Uii.transposeInPlace(); - Uaa.transposeInPlace(); - - // Transform occupied orbitals auto C_occ_view = C.block(num_molecular_orbitals * spin_index, 0, num_molecular_orbitals, num_occupied_orbitals); - RowMajorMatrix C_occ = C_occ_view.eval(); - blas::gemm(blas::Layout::RowMajor, blas::Op::NoTrans, blas::Op::NoTrans, - num_molecular_orbitals, num_occupied_orbitals, - num_occupied_orbitals, 1.0, C_occ.data(), num_occupied_orbitals, - Uii.data(), num_occupied_orbitals, 0.0, C_occ_view.data(), - num_molecular_orbitals); + calculate_pseudo_canonical_orbital_block( + Uii, num_occupied_orbitals, pseudo_canonical_eigenvalues_, 0, C_occ_view, + num_molecular_orbitals, num_molecular_orbitals); - // Transform virtual orbitals + RowMajorMatrix Uaa = F_MO.block(num_occupied_orbitals, num_occupied_orbitals, + num_virtual_orbitals, num_virtual_orbitals); auto C_virt_view = C.block(num_molecular_orbitals * spin_index, num_occupied_orbitals, num_molecular_orbitals, num_virtual_orbitals); - RowMajorMatrix C_virt = C_virt_view.eval(); - blas::gemm(blas::Layout::RowMajor, blas::Op::NoTrans, blas::Op::NoTrans, - num_molecular_orbitals, num_virtual_orbitals, num_virtual_orbitals, - 1.0, C_virt.data(), num_virtual_orbitals, Uaa.data(), - num_virtual_orbitals, 0.0, C_virt_view.data(), - num_molecular_orbitals); + calculate_pseudo_canonical_orbital_block( + Uaa, num_virtual_orbitals, pseudo_canonical_eigenvalues_, + num_occupied_orbitals, C_virt_view, num_molecular_orbitals, + num_molecular_orbitals); - // Transform the vectors in history_kappa and history_dgrad to - // accommodate current pseudo-canonical orbitals + // Transform the vectors in history_kappa, history_dgrad, and + // current_gradient_spin to accommodate current pseudo-canonical orbitals transform_history_(history_kappa_spin, Uii, Uaa, history_size_, num_occupied_orbitals, num_virtual_orbitals); transform_history_(history_dgrad_spin, Uii, Uaa, history_size_, num_occupied_orbitals, num_virtual_orbitals); - // Transform the gradient to accommodate current pseudo-canonical orbitals - RowMajorMatrix current_gradient_matrix = - Eigen::Map(current_gradient_spin.data(), - num_occupied_orbitals, num_virtual_orbitals); - RowMajorMatrix current_gradient_transformed_matrix = - Uii.transpose() * current_gradient_matrix * Uaa; - Eigen::VectorXd current_gradient_transformed = Eigen::Map( - current_gradient_transformed_matrix.data(), rotation_size); - current_gradient_spin = current_gradient_transformed; + rotate_gradient_to_pseudo_canonical_basis( + current_gradient_spin, 0, num_occupied_orbitals, num_virtual_orbitals, + Uii, Uaa, current_gradient_spin); } void GDM::build_restricted_unrestricted_pseudo_canonical_orbitals_hessian_( @@ -908,38 +981,6 @@ void GDM::build_restricted_unrestricted_pseudo_canonical_orbitals_hessian_( } } -static void diagonalize_block_rotate_orbitals( - RowMajorMatrix& input_block_output_eigenvectors, const int block_size, - Eigen::VectorXd& eigenvalues, const int eigenvalue_start_index, - Eigen::Block transformed_orbitals, - const int num_atomic_orbitals, const int num_molecular_orbitals) { - lapack::syev(lapack::Job::Vec, lapack::Uplo::Lower, block_size, - input_block_output_eigenvectors.data(), block_size, - eigenvalues.data() + eigenvalue_start_index); - input_block_output_eigenvectors.transposeInPlace(); - RowMajorMatrix copied_orbitals = transformed_orbitals; - blas::gemm(blas::Layout::RowMajor, blas::Op::NoTrans, blas::Op::NoTrans, - num_atomic_orbitals, block_size, block_size, 1.0, - copied_orbitals.data(), block_size, - input_block_output_eigenvectors.data(), block_size, 0.0, - transformed_orbitals.data(), num_molecular_orbitals); -} - -static void transform_gradient_vector(const Eigen::VectorXd& current_vector, - const int start_index, const int num_rows, - const int num_cols, - const RowMajorMatrix& u_left, - const RowMajorMatrix& u_right, - Eigen::VectorXd& transformed_vector) { - RowMajorMatrix current_matrix = Eigen::Map( - current_vector.data() + start_index, num_rows, num_cols); - RowMajorMatrix transformed_matrix = - u_left.transpose() * current_matrix * u_right; - transformed_vector.segment(start_index, num_rows * num_cols) = - Eigen::Map(transformed_matrix.data(), - num_rows * num_cols); -} - void GDM::generate_restricted_open_shell_pseudo_canonical_orbital_( const RowMajorMatrix& F, RowMajorMatrix& C, RowMajorMatrix& history_kappa, RowMajorMatrix& history_dgrad, Eigen::VectorXd& current_gradient) { @@ -964,44 +1005,44 @@ void GDM::generate_restricted_open_shell_pseudo_canonical_orbital_( F.block(num_atomic_orbitals, 0, num_atomic_orbitals, num_molecular_orbitals) * C; + + // Perform pseudo-canonical transformation + // Diagonalize occupied/virtual blocks and rotate orbitals to the + // pseudo-canonical basis. RowMajorMatrix Uii = 0.5 * (F_up_mo.block(0, 0, num_closed_orbitals, num_closed_orbitals) + F_dn_mo.block(0, 0, num_closed_orbitals, num_closed_orbitals)); + auto C_closed_view = C.block(0, 0, num_atomic_orbitals, num_closed_orbitals); + calculate_pseudo_canonical_orbital_block( + Uii, num_closed_orbitals, pseudo_canonical_eigenvalues_, 0, C_closed_view, + num_atomic_orbitals, num_molecular_orbitals); + RowMajorMatrix Uww = 0.5 * (F_up_mo.block(num_closed_orbitals, num_closed_orbitals, num_open_orbitals, num_open_orbitals) + F_dn_mo.block(num_closed_orbitals, num_closed_orbitals, num_open_orbitals, num_open_orbitals)); - RowMajorMatrix Uaa = - 0.5 * (F_up_mo.block(num_occupied_orbitals, num_occupied_orbitals, - num_virtual_orbitals, num_virtual_orbitals) + - F_dn_mo.block(num_occupied_orbitals, num_occupied_orbitals, - num_virtual_orbitals, num_virtual_orbitals)); - - // Compute eigenvalues/eigenvectors of occupied-occupied and virtual-virtual - // blocks for pseudo-canonical orbital transformation - // Then transform the orbitals to pseudo-canonical basis - auto C_closed_view = C.block(0, 0, num_atomic_orbitals, num_closed_orbitals); - diagonalize_block_rotate_orbitals( - Uii, num_closed_orbitals, pseudo_canonical_eigenvalues_, 0, C_closed_view, - num_atomic_orbitals, num_molecular_orbitals); - auto C_open_view = C.block(0, num_closed_orbitals, num_atomic_orbitals, num_open_orbitals); - diagonalize_block_rotate_orbitals( + calculate_pseudo_canonical_orbital_block( Uww, num_open_orbitals, pseudo_canonical_eigenvalues_, num_closed_orbitals, C_open_view, num_atomic_orbitals, num_molecular_orbitals); + RowMajorMatrix Uaa = + 0.5 * (F_up_mo.block(num_occupied_orbitals, num_occupied_orbitals, + num_virtual_orbitals, num_virtual_orbitals) + + F_dn_mo.block(num_occupied_orbitals, num_occupied_orbitals, + num_virtual_orbitals, num_virtual_orbitals)); auto C_virt_view = C.block(0, num_occupied_orbitals, num_atomic_orbitals, num_virtual_orbitals); - diagonalize_block_rotate_orbitals( + calculate_pseudo_canonical_orbital_block( Uaa, num_virtual_orbitals, pseudo_canonical_eigenvalues_, num_occupied_orbitals, C_virt_view, num_atomic_orbitals, num_molecular_orbitals); - // Transform the vectors in history_kappa and history_dgrad to - // accommodate current pseudo-canonical orbitals + // Transform the vectors in history_kappa, history_dgrad, and current_gradient + // to accommodate current pseudo-canonical orbitals int offset = 0; auto history_kappa_iw = history_kappa.block(0, 0, history_size_, size_closed_open); @@ -1032,24 +1073,23 @@ void GDM::generate_restricted_open_shell_pseudo_canonical_orbital_( transform_history_(history_dgrad_ia, Uii, Uaa, history_size_, num_closed_orbitals, num_virtual_orbitals); - // Transform the gradient to accommodate current pseudo-canonical orbitals Eigen::VectorXd current_gradient_transformed = Eigen::VectorXd(current_gradient.size()); offset = 0; - transform_gradient_vector(current_gradient, offset, num_closed_orbitals, - num_open_orbitals, Uii, Uww, - current_gradient_transformed); + rotate_gradient_to_pseudo_canonical_basis( + current_gradient, offset, num_closed_orbitals, num_open_orbitals, Uii, + Uww, current_gradient_transformed); offset += size_closed_open; - transform_gradient_vector(current_gradient, offset, num_open_orbitals, - num_virtual_orbitals, Uww, Uaa, - current_gradient_transformed); + rotate_gradient_to_pseudo_canonical_basis( + current_gradient, offset, num_open_orbitals, num_virtual_orbitals, Uww, + Uaa, current_gradient_transformed); offset += size_open_virtual; - transform_gradient_vector(current_gradient, offset, num_closed_orbitals, - num_virtual_orbitals, Uii, Uaa, - current_gradient_transformed); + rotate_gradient_to_pseudo_canonical_basis( + current_gradient, offset, num_closed_orbitals, num_virtual_orbitals, Uii, + Uaa, current_gradient_transformed); current_gradient = current_gradient_transformed; } @@ -1066,9 +1106,9 @@ void GDM::build_restricted_open_shell_pseudo_canonical_orbitals_hessian_( const int num_open_orbitals = num_electrons_[0] - num_closed_orbitals; const int num_virtual_orbitals = num_molecular_orbitals - num_electrons_[0]; - double initial_hessian_coeff = - 2.0; // between close and open shells, or between open and virtual - // orbitals, the coefficient should be 2.0 + // between close and open shells, or between open and virtual orbitals, the + // coefficient should be 2.0 + double initial_hessian_coeff = 2.0; int offset = 0; for (int j = 0; j < num_closed_orbitals; j++) { for (int v = 0; v < num_open_orbitals; v++) { @@ -1098,8 +1138,8 @@ void GDM::build_restricted_open_shell_pseudo_canonical_orbitals_hessian_( } offset += num_open_orbitals * num_virtual_orbitals; - initial_hessian_coeff = 4.0; // between closed and virtual orbitals, the - // coefficient should be 4.0 + // between closed and virtual orbitals, the coefficient should be 4.0 + initial_hessian_coeff = 4.0; for (int j = 0; j < num_closed_orbitals; j++) { for (int a = 0; a < num_virtual_orbitals; a++) { int index = offset + j * num_virtual_orbitals + a; From e7e5e89c34802accdb7df3d5542bf8b47e865d9f Mon Sep 17 00:00:00 2001 From: v-boqinzhang Date: Wed, 22 Apr 2026 15:29:45 -0700 Subject: [PATCH 11/29] add python tests --- python/tests/test_scf.py | 43 ++++++++++++++++++++++++++++++++++++++++ 1 file changed, 43 insertions(+) diff --git a/python/tests/test_scf.py b/python/tests/test_scf.py index 077cbc7f5..2bb229c17 100644 --- a/python/tests/test_scf.py +++ b/python/tests/test_scf.py @@ -273,6 +273,49 @@ def test_scf_solver_oh_rohf_diis(self): assert abs(energy - (-74.361530753176)) < scf_energy_tolerance assert orbitals.is_restricted() + def test_scf_solver_oh_rohf_incore_diis(self): + """Test SCF solver on OH system with ROHF/sto-3g using incore ERI.""" + oh_structure = create_oh_structure() + scf_solver = algorithms.create("scf_solver") + + scf_solver.settings().set("enable_gdm", False) + scf_solver.settings().set("method", "hf") + scf_solver.settings().set("scf_type", "restricted") + scf_solver.settings().set("eri_method", "incore") + + energy, wavefunction = scf_solver.run(oh_structure, 0, 2, "sto-3g") + orbitals = wavefunction.get_orbitals() + + assert abs(energy - (-74.361530753176)) < scf_energy_tolerance + assert orbitals.is_restricted() + + def test_scf_solver_oh_rohf_gdm(self): + """Test SCF solver on OH system with ROHF/sto-3g and GDM enabled.""" + oh_structure = create_oh_structure() + scf_solver = algorithms.create("scf_solver") + + scf_solver.settings().set("enable_gdm", True) + scf_solver.settings().set("method", "hf") + scf_solver.settings().set("scf_type", "restricted") + + energy, wavefunction = scf_solver.run(oh_structure, 0, 2, "sto-3g") + orbitals = wavefunction.get_orbitals() + + assert abs(energy - (-74.361530753176)) < scf_energy_tolerance + assert orbitals.is_restricted() + + def test_scf_solver_oh_roks_invalid(self): + """Test restricted open-shell KS request on OH doublet raises error.""" + oh_structure = create_oh_structure() + scf_solver = algorithms.create("scf_solver") + + scf_solver.settings().set("enable_gdm", True) + scf_solver.settings().set("method", "pbe") + scf_solver.settings().set("scf_type", "restricted") + + with pytest.raises(ValueError): + scf_solver.run(oh_structure, 0, 2, "sto-3g") + def test_scf_solver_oxygen_atom_gdm(self): """Test SCF solver on oxygen atom with PBE/cc-pvdz.""" oxygen = create_oxygen_structure() From dfc534c0d264ac134c9d203c46f2ea47d5108a4d Mon Sep 17 00:00:00 2001 From: v-boqinzhang Date: Wed, 22 Apr 2026 16:31:20 -0700 Subject: [PATCH 12/29] pre-commit --- python/tests/test_scf.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/python/tests/test_scf.py b/python/tests/test_scf.py index 2bb229c17..a52971309 100644 --- a/python/tests/test_scf.py +++ b/python/tests/test_scf.py @@ -313,7 +313,7 @@ def test_scf_solver_oh_roks_invalid(self): scf_solver.settings().set("method", "pbe") scf_solver.settings().set("scf_type", "restricted") - with pytest.raises(ValueError): + with pytest.raises(ValueError, match="Restricted open-shell calculations are only supported"): scf_solver.run(oh_structure, 0, 2, "sto-3g") def test_scf_solver_oxygen_atom_gdm(self): From bac2355f723dc1920915d3b48d49269103931ae8 Mon Sep 17 00:00:00 2001 From: v-boqinzhang Date: Thu, 23 Apr 2026 09:46:51 -0700 Subject: [PATCH 13/29] resolve comments 1 --- .../microsoft/scf/src/scf_algorithm/gdm.cpp | 2 +- cpp/tests/test_scf.cpp | 12 +++++------- 2 files changed, 6 insertions(+), 8 deletions(-) diff --git a/cpp/src/qdk/chemistry/algorithms/microsoft/scf/src/scf_algorithm/gdm.cpp b/cpp/src/qdk/chemistry/algorithms/microsoft/scf/src/scf_algorithm/gdm.cpp index e255b99aa..323e8b528 100644 --- a/cpp/src/qdk/chemistry/algorithms/microsoft/scf/src/scf_algorithm/gdm.cpp +++ b/cpp/src/qdk/chemistry/algorithms/microsoft/scf/src/scf_algorithm/gdm.cpp @@ -1106,7 +1106,7 @@ void GDM::build_restricted_open_shell_pseudo_canonical_orbitals_hessian_( const int num_open_orbitals = num_electrons_[0] - num_closed_orbitals; const int num_virtual_orbitals = num_molecular_orbitals - num_electrons_[0]; - // between close and open shells, or between open and virtual orbitals, the + // between closed and open shells, or between open and virtual orbitals, the // coefficient should be 2.0 double initial_hessian_coeff = 2.0; int offset = 0; diff --git a/cpp/tests/test_scf.cpp b/cpp/tests/test_scf.cpp index ed0a42beb..a10eb1d5b 100644 --- a/cpp/tests/test_scf.cpp +++ b/cpp/tests/test_scf.cpp @@ -157,9 +157,8 @@ TEST_F(ScfTest, OH_ROKS_invalid) { scf_solver->settings().set("method", "pbe"); scf_solver->settings().set("scf_type", "restricted"); - // Default should be a singlet - EXPECT_THROW(scf_solver->run(oh, 0, 2, "sto-3g"), - std::invalid_argument); // open-shell dublet + // Restricted ROKS should reject this open-shell doublet case. + EXPECT_THROW(scf_solver->run(oh, 0, 2, "sto-3g"), std::invalid_argument); } TEST_F(ScfTest, Oxygen_atom_gdm) { @@ -250,9 +249,8 @@ TEST_F(ScfTest, Oxygen_atom_invalid_energy_thresh_diis_switch_gdm) { scf_solver->settings().set("method", "pbe"); scf_solver->settings().set("enable_gdm", true); scf_solver->settings().set("energy_thresh_diis_switch", -2e-4); - // Default should be a singlet - EXPECT_THROW(scf_solver->run(oxygen, 0, 1, "cc-pvdz"), - std::invalid_argument); // open-shell dublet + + EXPECT_THROW(scf_solver->run(oxygen, 0, 1, "cc-pvdz"), std::invalid_argument); } TEST_F(ScfTest, Oxygen_atom_invalid_bfgs_history_size_limit_gdm) { @@ -262,7 +260,7 @@ TEST_F(ScfTest, Oxygen_atom_invalid_bfgs_history_size_limit_gdm) { scf_solver->settings().set("method", "pbe"); scf_solver->settings().set("enable_gdm", true); scf_solver->settings().set("gdm_bfgs_history_size_limit", 0); - // Default should be a singlet + EXPECT_THROW(scf_solver->run(oxygen, 0, 1, "cc-pvdz"), std::invalid_argument); } From 5c9b7a9f013cd7e821d3945ff0f14838d93c9083 Mon Sep 17 00:00:00 2001 From: v-boqinzhang Date: Mon, 27 Apr 2026 12:51:05 -0700 Subject: [PATCH 14/29] add function compute_atba_gemm --- .../microsoft/scf/src/scf_algorithm/gdm.cpp | 121 ++++++++++++++---- 1 file changed, 95 insertions(+), 26 deletions(-) diff --git a/cpp/src/qdk/chemistry/algorithms/microsoft/scf/src/scf_algorithm/gdm.cpp b/cpp/src/qdk/chemistry/algorithms/microsoft/scf/src/scf_algorithm/gdm.cpp index 323e8b528..1a9132c6c 100644 --- a/cpp/src/qdk/chemistry/algorithms/microsoft/scf/src/scf_algorithm/gdm.cpp +++ b/cpp/src/qdk/chemistry/algorithms/microsoft/scf/src/scf_algorithm/gdm.cpp @@ -28,6 +28,53 @@ namespace qdk::chemistry::scf { namespace impl { +/** + * @brief Compute C = A^T * B * A using BLAS GEMM only + * + * Matrix shapes: + * A: m x n + * B: m x m + * C: n x n + * + * @param[in] m Row count of A and dimension of square B block. + * @param[in] n Column count of A and dimension of square C block. + * + * Assumed leading dimensions for contiguous storage: + * RowMajor: lda = n, ldb = m, ldc = n + * ColMajor: lda = m, ldb = m, ldc = n + * + * Pointers may reference sub-block starts (not necessarily matrix origin). + * The intermediate workspace is allocated internally and does not reuse B. + */ +static void compute_atba_gemm(const double* A, const double* B, double* C, + int m, int n, + blas::Layout layout = blas::Layout::RowMajor) { + if (A == nullptr || B == nullptr || C == nullptr) { + throw std::invalid_argument("compute_atba_gemm: null matrix pointer."); + } + if (m < 0 || n < 0) { + throw std::invalid_argument("compute_atba_gemm: negative dimensions."); + } + if (m == 0 || n == 0) { + return; + } + + const int lda = (layout == blas::Layout::RowMajor) ? n : m; + const int ldb = m; + const int ldc = n; + + std::vector workspace(static_cast(m) * n); + const int ld_workspace = (layout == blas::Layout::RowMajor) ? n : m; + + // workspace = B * A + blas::gemm(layout, blas::Op::NoTrans, blas::Op::NoTrans, m, n, m, 1.0, B, ldb, + A, lda, 0.0, workspace.data(), ld_workspace); + + // C = A^T * workspace + blas::gemm(layout, blas::Op::Trans, blas::Op::NoTrans, n, n, m, 1.0, A, lda, + workspace.data(), ld_workspace, 0.0, C, ldc); +} + /** * @brief Construct the antisymmetric kappa matrix and apply C * exp(kappa) * @param[in,out] C Molecular orbital coefficient matrix @@ -178,13 +225,13 @@ static void compute_restricted_unrestricted_gradient( } RowMajorMatrix F_MO = - C.block(num_molecular_orbitals * spin_index, 0, num_molecular_orbitals, - num_molecular_orbitals) - .transpose() * - F.block(num_molecular_orbitals * spin_index, 0, num_molecular_orbitals, - num_molecular_orbitals) * - C.block(num_molecular_orbitals * spin_index, 0, num_molecular_orbitals, - num_molecular_orbitals); + RowMajorMatrix::Zero(num_molecular_orbitals, num_molecular_orbitals); + const int spin_block_offset = + num_molecular_orbitals * num_molecular_orbitals * spin_index; + const double* C_block_ptr = C.data() + spin_block_offset; + const double* F_block_ptr = F.data() + spin_block_offset; + compute_atba_gemm(C_block_ptr, F_block_ptr, F_MO.data(), + num_molecular_orbitals, num_molecular_orbitals); // Extract occupied-virtual block and compute gradient // The -4.0 before F_{ia} comes from derivative of energy w.r.t. kappa @@ -262,13 +309,33 @@ static void compute_restricted_open_shell_gradient( const auto C_ao_mo = C.block(0, 0, num_atomic_orbitals, num_molecular_orbitals); + const double* C_ao_mo_ptr = C_ao_mo.data(); // Calculate Generalized Fock matrix in MO basis - RowMajorMatrix H_mo = C_ao_mo.transpose() * H_ao * C_ao_mo; - RowMajorMatrix J_alpha_mo = C_ao_mo.transpose() * J_alpha_ao * C_ao_mo; - RowMajorMatrix J_beta_mo = C_ao_mo.transpose() * J_beta_ao * C_ao_mo; - RowMajorMatrix K_alpha_mo = C_ao_mo.transpose() * K_alpha_ao * C_ao_mo; - RowMajorMatrix K_beta_mo = C_ao_mo.transpose() * K_beta_ao * C_ao_mo; + RowMajorMatrix H_mo = + RowMajorMatrix::Zero(num_molecular_orbitals, num_molecular_orbitals); + compute_atba_gemm(C_ao_mo_ptr, H_ao.data(), H_mo.data(), num_atomic_orbitals, + num_molecular_orbitals); + + RowMajorMatrix J_alpha_mo = + RowMajorMatrix::Zero(num_molecular_orbitals, num_molecular_orbitals); + compute_atba_gemm(C_ao_mo_ptr, J_alpha_ao.data(), J_alpha_mo.data(), + num_atomic_orbitals, num_molecular_orbitals); + + RowMajorMatrix J_beta_mo = + RowMajorMatrix::Zero(num_molecular_orbitals, num_molecular_orbitals); + compute_atba_gemm(C_ao_mo_ptr, J_beta_ao.data(), J_beta_mo.data(), + num_atomic_orbitals, num_molecular_orbitals); + + RowMajorMatrix K_alpha_mo = + RowMajorMatrix::Zero(num_molecular_orbitals, num_molecular_orbitals); + compute_atba_gemm(C_ao_mo_ptr, K_alpha_ao.data(), K_alpha_mo.data(), + num_atomic_orbitals, num_molecular_orbitals); + + RowMajorMatrix K_beta_mo = + RowMajorMatrix::Zero(num_molecular_orbitals, num_molecular_orbitals); + compute_atba_gemm(C_ao_mo_ptr, K_beta_ao.data(), K_beta_mo.data(), + num_atomic_orbitals, num_molecular_orbitals); RowMajorMatrix F_I = H_mo + 2.0 * J_beta_mo - K_beta_mo; RowMajorMatrix F_A = J_alpha_mo - J_beta_mo - 0.5 * (K_alpha_mo - K_beta_mo); @@ -897,13 +964,13 @@ void GDM::generate_restricted_unrestricted_pseudo_canonical_orbital_( } RowMajorMatrix F_MO = - C.block(num_molecular_orbitals * spin_index, 0, num_molecular_orbitals, - num_molecular_orbitals) - .transpose() * - F.block(num_molecular_orbitals * spin_index, 0, num_molecular_orbitals, - num_molecular_orbitals) * - C.block(num_molecular_orbitals * spin_index, 0, num_molecular_orbitals, - num_molecular_orbitals); + RowMajorMatrix::Zero(num_molecular_orbitals, num_molecular_orbitals); + const int spin_block_offset = + num_molecular_orbitals * num_molecular_orbitals * spin_index; + const double* C_block_ptr = C.data() + spin_block_offset; + const double* F_block_ptr = F.data() + spin_block_offset; + compute_atba_gemm(C_block_ptr, F_block_ptr, F_MO.data(), + num_molecular_orbitals, num_molecular_orbitals); // Perform pseudo-canonical transformation // Diagonalize occupied/virtual blocks and rotate orbitals to the @@ -998,13 +1065,15 @@ void GDM::generate_restricted_open_shell_pseudo_canonical_orbital_( RowMajorMatrix::Zero(num_molecular_orbitals, num_molecular_orbitals); RowMajorMatrix F_dn_mo = F_up_mo; - F_up_mo.noalias() = - C.transpose() * - F.block(0, 0, num_atomic_orbitals, num_molecular_orbitals) * C; - F_dn_mo.noalias() = C.transpose() * - F.block(num_atomic_orbitals, 0, num_atomic_orbitals, - num_molecular_orbitals) * - C; + const double* C_block_ptr = C.data(); + const double* F_up_block_ptr = F.data(); + const double* F_dn_block_ptr = + F.data() + num_atomic_orbitals * num_molecular_orbitals; + + compute_atba_gemm(C_block_ptr, F_up_block_ptr, F_up_mo.data(), + num_atomic_orbitals, num_molecular_orbitals); + compute_atba_gemm(C_block_ptr, F_dn_block_ptr, F_dn_mo.data(), + num_atomic_orbitals, num_molecular_orbitals); // Perform pseudo-canonical transformation // Diagonalize occupied/virtual blocks and rotate orbitals to the From 9737c3c7ed2414c1ac2cb831af48dc44aa84e09f Mon Sep 17 00:00:00 2001 From: v-boqinzhang Date: Mon, 27 Apr 2026 13:12:44 -0700 Subject: [PATCH 15/29] mod rotate_gradient_to_pseudo_canonical_basis --- .../microsoft/scf/src/scf_algorithm/gdm.cpp | 15 +++++++++++++-- 1 file changed, 13 insertions(+), 2 deletions(-) diff --git a/cpp/src/qdk/chemistry/algorithms/microsoft/scf/src/scf_algorithm/gdm.cpp b/cpp/src/qdk/chemistry/algorithms/microsoft/scf/src/scf_algorithm/gdm.cpp index 1a9132c6c..a2ce0f638 100644 --- a/cpp/src/qdk/chemistry/algorithms/microsoft/scf/src/scf_algorithm/gdm.cpp +++ b/cpp/src/qdk/chemistry/algorithms/microsoft/scf/src/scf_algorithm/gdm.cpp @@ -934,8 +934,19 @@ static void rotate_gradient_to_pseudo_canonical_basis( Eigen::Ref transformed_vector) { RowMajorMatrix current_matrix = Eigen::Map( current_vector.data() + start_index, num_rows, num_cols); - RowMajorMatrix transformed_matrix = - u_left.transpose() * current_matrix * u_right; + RowMajorMatrix temp_matrix = RowMajorMatrix::Zero(num_rows, num_cols); + RowMajorMatrix transformed_matrix = RowMajorMatrix::Zero(num_rows, num_cols); + + // temp_matrix = current_matrix * u_right + blas::gemm(blas::Layout::RowMajor, blas::Op::NoTrans, blas::Op::NoTrans, + num_rows, num_cols, num_cols, 1.0, current_matrix.data(), num_cols, + u_right.data(), num_cols, 0.0, temp_matrix.data(), num_cols); + + // transformed_matrix = u_left^T * temp_matrix + blas::gemm(blas::Layout::RowMajor, blas::Op::Trans, blas::Op::NoTrans, + num_rows, num_cols, num_rows, 1.0, u_left.data(), num_rows, + temp_matrix.data(), num_cols, 0.0, transformed_matrix.data(), + num_cols); transformed_vector.segment(start_index, num_rows * num_cols) = Eigen::Map(transformed_matrix.data(), num_rows * num_cols); From 709344c51db61e14899923dbc5a51795ce6586b4 Mon Sep 17 00:00:00 2001 From: v-boqinzhang Date: Wed, 29 Apr 2026 14:24:13 -0700 Subject: [PATCH 16/29] move compute_atba_gemm --- .../qdk/chemistry/scf/core/scf_algorithm.h | 31 +++++++ .../microsoft/scf/src/scf_algorithm/gdm.cpp | 83 ++++++------------- .../scf/src/scf_algorithm/scf_algorithm.cpp | 55 ++++++++++-- 3 files changed, 106 insertions(+), 63 deletions(-) diff --git a/cpp/src/qdk/chemistry/algorithms/microsoft/scf/include/qdk/chemistry/scf/core/scf_algorithm.h b/cpp/src/qdk/chemistry/algorithms/microsoft/scf/include/qdk/chemistry/scf/core/scf_algorithm.h index 774bb4ddc..fbc3aa223 100644 --- a/cpp/src/qdk/chemistry/algorithms/microsoft/scf/include/qdk/chemistry/scf/core/scf_algorithm.h +++ b/cpp/src/qdk/chemistry/algorithms/microsoft/scf/include/qdk/chemistry/scf/core/scf_algorithm.h @@ -7,11 +7,42 @@ #include #include +#include #include #include +#include namespace qdk::chemistry::scf { +/** + * @brief Compute C = A^T * B * A using BLAS GEMM only + * + * Matrix shapes: + * A: m x n + * B: m x m + * C: n x n + * + * @param[in] m Row count of A and dimension of square B block. + * @param[in] n Column count of A and dimension of square C block. + * + * Assumed leading dimensions for contiguous storage: + * RowMajor: lda = n, ldb = m, ldc = n + * ColMajor: lda = m, ldb = m, ldc = n + * + * Pointers may reference sub-block starts (not necessarily matrix origin). + * The intermediate workspace is provided by caller to avoid repeated + * allocations across consecutive calls. + * + * @param[in] A Pointer to the A matrix block. + * @param[in] B Pointer to the B matrix block. + * @param[out] C Pointer to output C matrix block. + * @param[in,out] workspace Temporary buffer of at least m*n doubles. + * @param[in] layout BLAS matrix layout (RowMajor by default). + */ +void compute_atba_gemm(const double* A, const double* B, double* C, int m, + int n, std::vector& workspace, + blas::Layout layout = blas::Layout::RowMajor); + // Forward declaration class SCFImpl; /** diff --git a/cpp/src/qdk/chemistry/algorithms/microsoft/scf/src/scf_algorithm/gdm.cpp b/cpp/src/qdk/chemistry/algorithms/microsoft/scf/src/scf_algorithm/gdm.cpp index a2ce0f638..b7cefa6f4 100644 --- a/cpp/src/qdk/chemistry/algorithms/microsoft/scf/src/scf_algorithm/gdm.cpp +++ b/cpp/src/qdk/chemistry/algorithms/microsoft/scf/src/scf_algorithm/gdm.cpp @@ -17,6 +17,7 @@ #include "../scf/scf_impl.h" #include "line_search.h" #include "qdk/chemistry/scf/core/scf.h" +#include "qdk/chemistry/scf/core/scf_algorithm.h" #include "qdk/chemistry/scf/core/types.h" #include "util/matrix_exp.h" #ifdef QDK_CHEMISTRY_ENABLE_GPU @@ -28,53 +29,6 @@ namespace qdk::chemistry::scf { namespace impl { -/** - * @brief Compute C = A^T * B * A using BLAS GEMM only - * - * Matrix shapes: - * A: m x n - * B: m x m - * C: n x n - * - * @param[in] m Row count of A and dimension of square B block. - * @param[in] n Column count of A and dimension of square C block. - * - * Assumed leading dimensions for contiguous storage: - * RowMajor: lda = n, ldb = m, ldc = n - * ColMajor: lda = m, ldb = m, ldc = n - * - * Pointers may reference sub-block starts (not necessarily matrix origin). - * The intermediate workspace is allocated internally and does not reuse B. - */ -static void compute_atba_gemm(const double* A, const double* B, double* C, - int m, int n, - blas::Layout layout = blas::Layout::RowMajor) { - if (A == nullptr || B == nullptr || C == nullptr) { - throw std::invalid_argument("compute_atba_gemm: null matrix pointer."); - } - if (m < 0 || n < 0) { - throw std::invalid_argument("compute_atba_gemm: negative dimensions."); - } - if (m == 0 || n == 0) { - return; - } - - const int lda = (layout == blas::Layout::RowMajor) ? n : m; - const int ldb = m; - const int ldc = n; - - std::vector workspace(static_cast(m) * n); - const int ld_workspace = (layout == blas::Layout::RowMajor) ? n : m; - - // workspace = B * A - blas::gemm(layout, blas::Op::NoTrans, blas::Op::NoTrans, m, n, m, 1.0, B, ldb, - A, lda, 0.0, workspace.data(), ld_workspace); - - // C = A^T * workspace - blas::gemm(layout, blas::Op::Trans, blas::Op::NoTrans, n, n, m, 1.0, A, lda, - workspace.data(), ld_workspace, 0.0, C, ldc); -} - /** * @brief Construct the antisymmetric kappa matrix and apply C * exp(kappa) * @param[in,out] C Molecular orbital coefficient matrix @@ -214,6 +168,9 @@ static void compute_restricted_unrestricted_gradient( } gradient.setZero(total_rotation_size); + std::vector atba_workspace( + static_cast(num_molecular_orbitals) * num_molecular_orbitals); + for (int spin_index = 0; spin_index < num_orbital_spin_blocks; ++spin_index) { const int num_occupied_orbitals = num_electrons[spin_index]; const int num_virtual_orbitals = @@ -231,7 +188,8 @@ static void compute_restricted_unrestricted_gradient( const double* C_block_ptr = C.data() + spin_block_offset; const double* F_block_ptr = F.data() + spin_block_offset; compute_atba_gemm(C_block_ptr, F_block_ptr, F_MO.data(), - num_molecular_orbitals, num_molecular_orbitals); + num_molecular_orbitals, num_molecular_orbitals, + atba_workspace); // Extract occupied-virtual block and compute gradient // The -4.0 before F_{ia} comes from derivative of energy w.r.t. kappa @@ -310,32 +268,38 @@ static void compute_restricted_open_shell_gradient( const auto C_ao_mo = C.block(0, 0, num_atomic_orbitals, num_molecular_orbitals); const double* C_ao_mo_ptr = C_ao_mo.data(); + std::vector atba_workspace(static_cast(num_atomic_orbitals) * + num_molecular_orbitals); // Calculate Generalized Fock matrix in MO basis RowMajorMatrix H_mo = RowMajorMatrix::Zero(num_molecular_orbitals, num_molecular_orbitals); compute_atba_gemm(C_ao_mo_ptr, H_ao.data(), H_mo.data(), num_atomic_orbitals, - num_molecular_orbitals); + num_molecular_orbitals, atba_workspace); RowMajorMatrix J_alpha_mo = RowMajorMatrix::Zero(num_molecular_orbitals, num_molecular_orbitals); compute_atba_gemm(C_ao_mo_ptr, J_alpha_ao.data(), J_alpha_mo.data(), - num_atomic_orbitals, num_molecular_orbitals); + num_atomic_orbitals, num_molecular_orbitals, + atba_workspace); RowMajorMatrix J_beta_mo = RowMajorMatrix::Zero(num_molecular_orbitals, num_molecular_orbitals); compute_atba_gemm(C_ao_mo_ptr, J_beta_ao.data(), J_beta_mo.data(), - num_atomic_orbitals, num_molecular_orbitals); + num_atomic_orbitals, num_molecular_orbitals, + atba_workspace); RowMajorMatrix K_alpha_mo = RowMajorMatrix::Zero(num_molecular_orbitals, num_molecular_orbitals); compute_atba_gemm(C_ao_mo_ptr, K_alpha_ao.data(), K_alpha_mo.data(), - num_atomic_orbitals, num_molecular_orbitals); + num_atomic_orbitals, num_molecular_orbitals, + atba_workspace); RowMajorMatrix K_beta_mo = RowMajorMatrix::Zero(num_molecular_orbitals, num_molecular_orbitals); compute_atba_gemm(C_ao_mo_ptr, K_beta_ao.data(), K_beta_mo.data(), - num_atomic_orbitals, num_molecular_orbitals); + num_atomic_orbitals, num_molecular_orbitals, + atba_workspace); RowMajorMatrix F_I = H_mo + 2.0 * J_beta_mo - K_beta_mo; RowMajorMatrix F_A = J_alpha_mo - J_beta_mo - 0.5 * (K_alpha_mo - K_beta_mo); @@ -980,8 +944,11 @@ void GDM::generate_restricted_unrestricted_pseudo_canonical_orbital_( num_molecular_orbitals * num_molecular_orbitals * spin_index; const double* C_block_ptr = C.data() + spin_block_offset; const double* F_block_ptr = F.data() + spin_block_offset; + std::vector atba_workspace( + static_cast(num_molecular_orbitals) * num_molecular_orbitals); compute_atba_gemm(C_block_ptr, F_block_ptr, F_MO.data(), - num_molecular_orbitals, num_molecular_orbitals); + num_molecular_orbitals, num_molecular_orbitals, + atba_workspace); // Perform pseudo-canonical transformation // Diagonalize occupied/virtual blocks and rotate orbitals to the @@ -1080,11 +1047,15 @@ void GDM::generate_restricted_open_shell_pseudo_canonical_orbital_( const double* F_up_block_ptr = F.data(); const double* F_dn_block_ptr = F.data() + num_atomic_orbitals * num_molecular_orbitals; + std::vector atba_workspace(static_cast(num_atomic_orbitals) * + num_molecular_orbitals); compute_atba_gemm(C_block_ptr, F_up_block_ptr, F_up_mo.data(), - num_atomic_orbitals, num_molecular_orbitals); + num_atomic_orbitals, num_molecular_orbitals, + atba_workspace); compute_atba_gemm(C_block_ptr, F_dn_block_ptr, F_dn_mo.data(), - num_atomic_orbitals, num_molecular_orbitals); + num_atomic_orbitals, num_molecular_orbitals, + atba_workspace); // Perform pseudo-canonical transformation // Diagonalize occupied/virtual blocks and rotate orbitals to the diff --git a/cpp/src/qdk/chemistry/algorithms/microsoft/scf/src/scf_algorithm/scf_algorithm.cpp b/cpp/src/qdk/chemistry/algorithms/microsoft/scf/src/scf_algorithm/scf_algorithm.cpp index eaadc5cc6..07152eec3 100644 --- a/cpp/src/qdk/chemistry/algorithms/microsoft/scf/src/scf_algorithm/scf_algorithm.cpp +++ b/cpp/src/qdk/chemistry/algorithms/microsoft/scf/src/scf_algorithm/scf_algorithm.cpp @@ -32,8 +32,44 @@ #include #endif +#include + namespace qdk::chemistry::scf { +void compute_atba_gemm(const double* A, const double* B, double* C, int m, + int n, std::vector& workspace, + blas::Layout layout) { + if (A == nullptr || B == nullptr || C == nullptr) { + throw std::invalid_argument("compute_atba_gemm: null matrix pointer."); + } + if (m < 0 || n < 0) { + throw std::invalid_argument("compute_atba_gemm: negative dimensions."); + } + if (m == 0 || n == 0) { + return; + } + + const size_t required_workspace_size = static_cast(m) * n; + if (workspace.size() < required_workspace_size) { + throw std::invalid_argument( + "compute_atba_gemm: workspace is smaller than m * n."); + } + + const int lda = (layout == blas::Layout::RowMajor) ? n : m; + const int ldb = m; + const int ldc = n; + + const int ld_workspace = (layout == blas::Layout::RowMajor) ? n : m; + + // workspace = B * A + blas::gemm(layout, blas::Op::NoTrans, blas::Op::NoTrans, m, n, m, 1.0, B, ldb, + A, lda, 0.0, workspace.data(), ld_workspace); + + // C = A^T * workspace + blas::gemm(layout, blas::Op::Trans, blas::Op::NoTrans, n, n, m, 1.0, A, lda, + workspace.data(), ld_workspace, 0.0, C, ldc); +} + SCFAlgorithm::SCFAlgorithm(const SCFContext& ctx) : ctx_(ctx), step_count_(0), @@ -252,13 +288,18 @@ void SCFAlgorithm::build_rohf_f_p_matrix(const RowMajorMatrix& F, RowMajorMatrix F_dn_mo = F_up_mo; RowMajorMatrix effective_F_mo = F_up_mo; - F_up_mo.noalias() = C.transpose() * - F.block(0, 0, num_atomic_orbitals, num_atomic_orbitals) * - C; - F_dn_mo.noalias() = C.transpose() * - F.block(num_atomic_orbitals, 0, num_atomic_orbitals, - num_atomic_orbitals) * - C; + const double* C_block_ptr = C.data(); + const double* F_up_block_ptr = F.data(); + const double* F_dn_block_ptr = + F.data() + num_atomic_orbitals * num_atomic_orbitals; + std::vector atba_workspace(static_cast(num_atomic_orbitals) * + num_molecular_orbitals); + compute_atba_gemm(C_block_ptr, F_up_block_ptr, F_up_mo.data(), + num_atomic_orbitals, num_molecular_orbitals, atba_workspace, + blas::Layout::RowMajor); + compute_atba_gemm(C_block_ptr, F_dn_block_ptr, F_dn_mo.data(), + num_atomic_orbitals, num_molecular_orbitals, atba_workspace, + blas::Layout::RowMajor); auto average_block = [&effective_F_mo, &F_up_mo, &F_dn_mo]( int row, int col, int rows, int cols) { From e15cb9547167bfcf4302b872dc9e7fb22930b03d Mon Sep 17 00:00:00 2001 From: v-boqinzhang Date: Wed, 29 Apr 2026 15:26:01 -0700 Subject: [PATCH 17/29] mod update_density_matrix --- .../qdk/chemistry/scf/core/scf_algorithm.h | 10 +++-- .../microsoft/scf/src/scf_algorithm/diis.cpp | 44 ------------------- .../microsoft/scf/src/scf_algorithm/diis.h | 16 ------- .../scf/src/scf_algorithm/scf_algorithm.cpp | 30 +++++++++++++ 4 files changed, 36 insertions(+), 64 deletions(-) diff --git a/cpp/src/qdk/chemistry/algorithms/microsoft/scf/include/qdk/chemistry/scf/core/scf_algorithm.h b/cpp/src/qdk/chemistry/algorithms/microsoft/scf/include/qdk/chemistry/scf/core/scf_algorithm.h index fbc3aa223..a84d854ba 100644 --- a/cpp/src/qdk/chemistry/algorithms/microsoft/scf/include/qdk/chemistry/scf/core/scf_algorithm.h +++ b/cpp/src/qdk/chemistry/algorithms/microsoft/scf/include/qdk/chemistry/scf/core/scf_algorithm.h @@ -129,10 +129,12 @@ class SCFAlgorithm { int num_molecular_orbitals, int idx_spin); /** - * @brief Update the density matrix for restricted or unrestricted - * calculations. For ASAHF and ROHF calculations, this method will be - * overridden to implement the specific density matrix construction for those - * methods. + * @brief Update the density matrix for RHF, UHF, or ROHF calculations. + * + * The base implementation handles standard RHF/UHF construction and also + * reconstructs spin-blocked alpha/beta densities for ROHF from a shared + * molecular-orbital coefficient matrix. Algorithms with custom density + * construction (for example ASAHF) can still override this method. * * @param[in,out] P Reference to the density matrix to be updated * @param[in] C Reference to the molecular orbital coefficients diff --git a/cpp/src/qdk/chemistry/algorithms/microsoft/scf/src/scf_algorithm/diis.cpp b/cpp/src/qdk/chemistry/algorithms/microsoft/scf/src/scf_algorithm/diis.cpp index 4a35dfb3d..e10466388 100644 --- a/cpp/src/qdk/chemistry/algorithms/microsoft/scf/src/scf_algorithm/diis.cpp +++ b/cpp/src/qdk/chemistry/algorithms/microsoft/scf/src/scf_algorithm/diis.cpp @@ -118,38 +118,6 @@ class DIIS { std::numeric_limits::infinity(); ///< Current DIIS error }; -/** - * @brief Reconstruct spin-blocked densities from the ROHF MO matrix - * - * Generates $P_\alpha$ and $P_\beta$ blocks so we can hand the updated - * density back to SCFImpl after diagonalization. - * - * @param P Spin-blocked density matrix to overwrite - * @param C Molecular-orbital coefficients from latest diagonalization - * @param nelec_alpha Number of alpha electrons - * @param nelec_beta Number of beta electrons - */ -void update_rohf_density_matrix(RowMajorMatrix& P, const RowMajorMatrix& C, - int nelec_alpha, int nelec_beta) { - QDK_LOG_TRACE_ENTERING(); - const int num_atomic_orbitals = static_cast(C.rows()); - - auto build_density = [&](auto&& target, int n_occ) { - if (n_occ <= 0) { - target.setZero(); - return; - } - target.noalias() = C.block(0, 0, num_atomic_orbitals, n_occ) * - C.block(0, 0, num_atomic_orbitals, n_occ).transpose(); - }; - - auto P_alpha = P.block(0, 0, num_atomic_orbitals, num_atomic_orbitals); - auto P_beta = - P.block(num_atomic_orbitals, 0, num_atomic_orbitals, num_atomic_orbitals); - build_density(P_alpha, nelec_alpha); - build_density(P_beta, nelec_beta); -} - DIIS::DIIS(const SCFContext& ctx, const size_t subspace_size) : ctx_(ctx), subspace_size_(subspace_size) { QDK_LOG_TRACE_ENTERING(); @@ -352,18 +320,6 @@ void DIIS::iterate(SCFImpl& scf_impl) { } } -void DIIS::update_density_matrix(RowMajorMatrix& P, const RowMajorMatrix& C, - bool unrestricted, int nelec_alpha, - int nelec_beta) { - QDK_LOG_TRACE_ENTERING(); - if (ctx_.cfg->scf_orbital_type == SCFOrbitalType::RestrictedOpenShell) { - impl::update_rohf_density_matrix(P, C, nelec_alpha, nelec_beta); - return; - } - SCFAlgorithm::update_density_matrix(P, C, unrestricted, nelec_alpha, - nelec_beta); -} - double DIIS::current_diis_error() const { QDK_LOG_TRACE_ENTERING(); return diis_impl_->get_diis_error(); diff --git a/cpp/src/qdk/chemistry/algorithms/microsoft/scf/src/scf_algorithm/diis.h b/cpp/src/qdk/chemistry/algorithms/microsoft/scf/src/scf_algorithm/diis.h index 261d56c56..3cbf4781f 100644 --- a/cpp/src/qdk/chemistry/algorithms/microsoft/scf/src/scf_algorithm/diis.h +++ b/cpp/src/qdk/chemistry/algorithms/microsoft/scf/src/scf_algorithm/diis.h @@ -54,22 +54,6 @@ class DIIS : public SCFAlgorithm { */ void iterate(SCFImpl& scf_impl) override; - /** - * @brief Update the density matrix - * - * Default implementation handles restricted and unrestricted cases, while - * subclasses such as ASAHF may override it - * - * @param[in,out] P Density matrix to overwrite - * @param[in] C Molecular orbital coefficient matrix - * @param[in] unrestricted True if two spin blocks are present - * @param[in] nelec_alpha Number of alpha electrons - * @param[in] nelec_beta Number of beta electrons - */ - void update_density_matrix(RowMajorMatrix& P, const RowMajorMatrix& C, - bool unrestricted, int nelec_alpha, - int nelec_beta) override; - private: /** * @brief Return most recent Pulay error metric for damping logic diff --git a/cpp/src/qdk/chemistry/algorithms/microsoft/scf/src/scf_algorithm/scf_algorithm.cpp b/cpp/src/qdk/chemistry/algorithms/microsoft/scf/src/scf_algorithm/scf_algorithm.cpp index 07152eec3..1332e78ca 100644 --- a/cpp/src/qdk/chemistry/algorithms/microsoft/scf/src/scf_algorithm/scf_algorithm.cpp +++ b/cpp/src/qdk/chemistry/algorithms/microsoft/scf/src/scf_algorithm/scf_algorithm.cpp @@ -206,6 +206,36 @@ void SCFAlgorithm::update_density_matrix(RowMajorMatrix& P, bool unrestricted, int nelec_alpha, int nelec_beta) { QDK_LOG_TRACE_ENTERING(); + if (ctx_.cfg->scf_orbital_type == SCFOrbitalType::RestrictedOpenShell) { + const int num_atomic_orbitals = + static_cast(ctx_.basis_set->num_atomic_orbitals); + if (C.rows() != num_atomic_orbitals) { + throw std::invalid_argument( + "ROHF coefficient matrix row count does not match AO dimension"); + } + if (P.rows() != 2 * num_atomic_orbitals || + P.cols() != num_atomic_orbitals) { + throw std::invalid_argument( + "ROHF density matrix must contain alpha and beta AO blocks"); + } + + auto build_density_block = [&](auto&& target, int n_occ) { + if (n_occ <= 0) { + target.setZero(); + return; + } + target.noalias() = C.block(0, 0, num_atomic_orbitals, n_occ) * + C.block(0, 0, num_atomic_orbitals, n_occ).transpose(); + }; + + auto P_alpha = P.block(0, 0, num_atomic_orbitals, num_atomic_orbitals); + auto P_beta = P.block(num_atomic_orbitals, 0, num_atomic_orbitals, + num_atomic_orbitals); + build_density_block(P_alpha, nelec_alpha); + build_density_block(P_beta, nelec_beta); + return; + } + const int num_orbital_sets = unrestricted ? 2 : 1; const int num_atomic_orbitals = static_cast(ctx_.basis_set->num_atomic_orbitals); From caf2634c79615ead1d262abe17f7630b31b242b7 Mon Sep 17 00:00:00 2001 From: v-boqinzhang Date: Thu, 30 Apr 2026 10:31:51 -0700 Subject: [PATCH 18/29] resolve comments 2 --- .../microsoft/scf/src/scf_algorithm/gdm.cpp | 113 +++++++++++------- 1 file changed, 70 insertions(+), 43 deletions(-) diff --git a/cpp/src/qdk/chemistry/algorithms/microsoft/scf/src/scf_algorithm/gdm.cpp b/cpp/src/qdk/chemistry/algorithms/microsoft/scf/src/scf_algorithm/gdm.cpp index b7cefa6f4..7abc906d4 100644 --- a/cpp/src/qdk/chemistry/algorithms/microsoft/scf/src/scf_algorithm/gdm.cpp +++ b/cpp/src/qdk/chemistry/algorithms/microsoft/scf/src/scf_algorithm/gdm.cpp @@ -35,13 +35,22 @@ namespace impl { * @param[in] spin_index Spin index (0 for alpha, 1 for beta) * @param[in] kappa_vector The kappa vector to apply for rotation * @param[in] num_occupied_orbitals Number of occupied orbitals for this spin - * @param[in] num_molecular_orbitals Number of molecular orbitals + * @param[in] num_orbital_spin_blocks Number of spin blocks (1 for RHF, 2 for + * UHF) */ static void apply_restricted_unrestricted_orbital_rotation( RowMajorMatrix& C, const int spin_index, const Eigen::VectorXd& kappa_vector, const int num_occupied_orbitals, - const int num_molecular_orbitals) { + const int num_orbital_spin_blocks) { QDK_LOG_TRACE_ENTERING(); + const int num_molecular_orbitals = C.cols(); + const int num_atomic_orbitals = + (num_orbital_spin_blocks == 2) ? C.rows() / 2 : C.rows(); + + if (num_orbital_spin_blocks == 2 && (C.rows() % 2 != 0)) { + throw std::invalid_argument("In UHF or UKS, C row count must be even!"); + } + const int num_virtual_orbitals = num_molecular_orbitals - num_occupied_orbitals; @@ -61,18 +70,17 @@ static void apply_restricted_unrestricted_orbital_rotation( matrix_exp(kappa_complete.data(), exp_kappa.data(), num_molecular_orbitals); // Rotate C: C' = C * exp(kappa) + const int C_block_offset = + num_atomic_orbitals * num_molecular_orbitals * spin_index; RowMajorMatrix C_before_rotate = - C.block(num_molecular_orbitals * spin_index, 0, num_molecular_orbitals, + C.block(spin_index * num_atomic_orbitals, 0, num_atomic_orbitals, num_molecular_orbitals); + const double* C_before_rotate_ptr = C_before_rotate.data(); blas::gemm(blas::Layout::RowMajor, blas::Op::NoTrans, blas::Op::NoTrans, - num_molecular_orbitals, num_molecular_orbitals, - num_molecular_orbitals, 1.0, C_before_rotate.data(), + num_atomic_orbitals, num_molecular_orbitals, + num_molecular_orbitals, 1.0, C_before_rotate_ptr, num_molecular_orbitals, exp_kappa.data(), num_molecular_orbitals, - 0.0, - C.block(num_molecular_orbitals * spin_index, 0, - num_molecular_orbitals, num_molecular_orbitals) - .data(), - num_molecular_orbitals); + 0.0, C.data() + C_block_offset, num_molecular_orbitals); } /** @@ -81,12 +89,13 @@ static void apply_restricted_unrestricted_orbital_rotation( * @param[in,out] C Molecular orbital coefficient matrix * @param[in] num_electrons Occupied orbital counts (alpha, beta) * @param[in] kappa_vector Concatenated ROHF rotation parameters - * @param[in] num_molecular_orbitals Total molecular orbitals in the system */ static void apply_restricted_open_shell_orbital_rotation( RowMajorMatrix& C, const std::vector& num_electrons, - const Eigen::VectorXd& kappa_vector, const int num_molecular_orbitals) { + const Eigen::VectorXd& kappa_vector) { QDK_LOG_TRACE_ENTERING(); + const int num_atomic_orbitals = C.rows(); + const int num_molecular_orbitals = C.cols(); const int num_closed_orbitals = num_electrons[1]; const int num_open_orbitals = num_electrons[0] - num_closed_orbitals; const int num_virtual_orbitals = num_molecular_orbitals - num_electrons[0]; @@ -138,7 +147,7 @@ static void apply_restricted_open_shell_orbital_rotation( RowMajorMatrix C_before_rotate = C; blas::gemm(blas::Layout::RowMajor, blas::Op::NoTrans, blas::Op::NoTrans, - num_molecular_orbitals, num_molecular_orbitals, + num_atomic_orbitals, num_molecular_orbitals, num_molecular_orbitals, 1.0, C_before_rotate.data(), num_molecular_orbitals, exp_kappa.data(), num_molecular_orbitals, 0.0, C.data(), num_molecular_orbitals); @@ -153,7 +162,6 @@ static void apply_restricted_open_shell_orbital_rotation( * @param[in] rotation_size Number of rotation parameters per spin * (n_occ*n_virt) * @param[in] num_orbital_spin_blocks Number of spin blocks to iterate - * @param[in] num_molecular_orbitals Total molecular orbitals in the system * @param[out] gradient Output gradient vector (concatenated across spins) */ static void compute_restricted_unrestricted_gradient( @@ -161,7 +169,15 @@ static void compute_restricted_unrestricted_gradient( const std::vector& num_electrons, const std::vector& rotation_offset, const std::vector& rotation_size, int num_orbital_spin_blocks, - int num_molecular_orbitals, Eigen::VectorXd& gradient) { + Eigen::VectorXd& gradient) { + const int num_molecular_orbitals = C.cols(); + const int num_atomic_orbitals = + (num_orbital_spin_blocks == 2) ? C.rows() / 2 : C.rows(); + + if (num_orbital_spin_blocks == 2 && (C.rows() % 2 != 0)) { + throw std::invalid_argument("In UHF or UKS, C row count must be even!"); + } + int total_rotation_size = 0; for (int i = 0; i < num_orbital_spin_blocks; ++i) { total_rotation_size += rotation_size[i]; @@ -183,12 +199,14 @@ static void compute_restricted_unrestricted_gradient( RowMajorMatrix F_MO = RowMajorMatrix::Zero(num_molecular_orbitals, num_molecular_orbitals); - const int spin_block_offset = - num_molecular_orbitals * num_molecular_orbitals * spin_index; - const double* C_block_ptr = C.data() + spin_block_offset; - const double* F_block_ptr = F.data() + spin_block_offset; + const int C_block_offset = + num_atomic_orbitals * num_molecular_orbitals * spin_index; + const int F_block_offset = + num_atomic_orbitals * num_atomic_orbitals * spin_index; + const double* C_block_ptr = C.data() + C_block_offset; + const double* F_block_ptr = F.data() + F_block_offset; compute_atba_gemm(C_block_ptr, F_block_ptr, F_MO.data(), - num_molecular_orbitals, num_molecular_orbitals, + num_atomic_orbitals, num_molecular_orbitals, atba_workspace); // Extract occupied-virtual block and compute gradient @@ -353,11 +371,10 @@ class GDMLineFunctor { /** * @brief Bind functor to a specific SCF state for line search evaluations. * @param scf_impl Reference to `SCFImpl` used to evaluate trial densities - * @param C_pseudo_canonical Molecular orbitals in pseudo-canonical basis + * @param C_pseudo_canonical Pseudo-canonical molecular orbital coefficients * @param num_electrons Occupied orbital counts per spin component * @param rotation_offset Starting index for each spin's rotation slice * @param rotation_size Number of rotation parameters per spin (n_occ*n_virt) - * @param num_molecular_orbitals Total molecular orbitals in the system * @param scf_orbital_type Spin symmetry used across SCF algorithms */ GDMLineFunctor(const SCFImpl& scf_impl, @@ -365,7 +382,7 @@ class GDMLineFunctor { const std::vector& num_electrons, const std::vector& rotation_offset, const std::vector& rotation_size, - int num_molecular_orbitals, SCFOrbitalType scf_orbital_type) + SCFOrbitalType scf_orbital_type) : scf_impl_(scf_impl), C_pseudo_canonical_(C_pseudo_canonical), num_electrons_(num_electrons), @@ -375,10 +392,19 @@ class GDMLineFunctor { scf_orbital_type == SCFOrbitalType::Unrestricted ? 2 : 1), num_density_matrices_( scf_orbital_type == SCFOrbitalType::Restricted ? 1 : 2), - num_molecular_orbitals_(num_molecular_orbitals), + num_atomic_orbitals_(scf_orbital_type == SCFOrbitalType::Unrestricted + ? static_cast(C_pseudo_canonical.rows()) / + 2 + : static_cast(C_pseudo_canonical.rows())), + num_molecular_orbitals_(static_cast(C_pseudo_canonical.cols())), scf_orbital_type_(scf_orbital_type), cached_kappa_(Eigen::VectorXd()), - cached_energy_(std::numeric_limits::infinity()) {} + cached_energy_(std::numeric_limits::infinity()) { + if (scf_orbital_type == SCFOrbitalType::Unrestricted && + (C_pseudo_canonical.rows() % 2 != 0)) { + throw std::invalid_argument("In UHF or UKS, C row count must be even!"); + } + } /** * @brief Evaluate energy at given kappa vector x @@ -430,6 +456,7 @@ class GDMLineFunctor { // Value parameters const int num_orbital_spin_blocks_; const int num_density_matrices_; + const int num_atomic_orbitals_; const int num_molecular_orbitals_; const SCFOrbitalType scf_orbital_type_; @@ -454,45 +481,45 @@ double GDMLineFunctor::eval(const Eigen::VectorXd& x) { // Apply rotation for all spins with kappa_trial if (scf_orbital_type_ == SCFOrbitalType::RestrictedOpenShell) { - apply_restricted_open_shell_orbital_rotation( - cached_C_, num_electrons_, kappa_trial, num_molecular_orbitals_); + apply_restricted_open_shell_orbital_rotation(cached_C_, num_electrons_, + kappa_trial); } else { for (int i = 0; i < num_orbital_spin_blocks_; i++) { auto kappa_spin = kappa_trial.segment(rotation_offset_[i], rotation_size_[i]); - apply_restricted_unrestricted_orbital_rotation( - cached_C_, i, kappa_spin, num_electrons_[i], num_molecular_orbitals_); + apply_restricted_unrestricted_orbital_rotation(cached_C_, i, kappa_spin, + num_electrons_[i], + num_orbital_spin_blocks_); } } // Compute P_trial from rotated C (for all spins) - cached_P_ = RowMajorMatrix::Zero( - num_density_matrices_ * num_molecular_orbitals_, num_molecular_orbitals_); + cached_P_ = RowMajorMatrix::Zero(num_density_matrices_ * num_atomic_orbitals_, + num_atomic_orbitals_); for (int i = 0; i < num_density_matrices_; i++) { const int num_occupied_orbitals = num_electrons_[i]; const double occupation_factor = (scf_orbital_type_ == SCFOrbitalType::Restricted) ? 2.0 : 1.0; - auto P_block = - cached_P_.block(num_molecular_orbitals_ * i, 0, num_molecular_orbitals_, - num_molecular_orbitals_); + auto P_block = cached_P_.block(num_atomic_orbitals_ * i, 0, + num_atomic_orbitals_, num_atomic_orbitals_); double* C_block_data = nullptr; if (scf_orbital_type_ == SCFOrbitalType::RestrictedOpenShell) { auto C_block = - cached_C_.block(0, 0, num_molecular_orbitals_, num_occupied_orbitals); + cached_C_.block(0, 0, num_atomic_orbitals_, num_occupied_orbitals); C_block_data = C_block.data(); } else { auto C_block = - cached_C_.block(num_molecular_orbitals_ * i, 0, - num_molecular_orbitals_, num_occupied_orbitals); + cached_C_.block(num_atomic_orbitals_ * i, 0, num_atomic_orbitals_, + num_occupied_orbitals); C_block_data = C_block.data(); } blas::gemm(blas::Layout::RowMajor, blas::Op::NoTrans, blas::Op::Trans, - num_molecular_orbitals_, num_molecular_orbitals_, + num_atomic_orbitals_, num_atomic_orbitals_, num_occupied_orbitals, occupation_factor, C_block_data, num_molecular_orbitals_, C_block_data, num_molecular_orbitals_, - 0.0, P_block.data(), num_molecular_orbitals_); + 0.0, P_block.data(), num_atomic_orbitals_); } // Evaluate energy and Fock matrix using trial density matrix @@ -524,7 +551,7 @@ Eigen::VectorXd GDMLineFunctor::grad(const Eigen::VectorXd& x) { } else { compute_restricted_unrestricted_gradient( cached_F_, cached_C_, num_electrons_, rotation_offset_, rotation_size_, - num_orbital_spin_blocks_, num_molecular_orbitals_, gradient); + num_orbital_spin_blocks_, gradient); } return gradient; @@ -1046,7 +1073,7 @@ void GDM::generate_restricted_open_shell_pseudo_canonical_orbital_( const double* C_block_ptr = C.data(); const double* F_up_block_ptr = F.data(); const double* F_dn_block_ptr = - F.data() + num_atomic_orbitals * num_molecular_orbitals; + F.data() + num_atomic_orbitals * num_atomic_orbitals; std::vector atba_workspace(static_cast(num_atomic_orbitals) * num_molecular_orbitals); @@ -1232,7 +1259,7 @@ void GDM::iterate(SCFImpl& scf_impl) { } else { compute_restricted_unrestricted_gradient( F, C, num_electrons_, rotation_offset_, rotation_size_, - num_orbital_spin_blocks_, num_molecular_orbitals, current_gradient_); + num_orbital_spin_blocks_, current_gradient_); } if (gdm_step_count_ != 0) { @@ -1383,7 +1410,7 @@ void GDM::iterate(SCFImpl& scf_impl) { // Create line search functor for energy evaluation GDMLineFunctor line_functor(scf_impl, C_pseudo_canonical, num_electrons_, rotation_offset_, rotation_size_, - num_molecular_orbitals, cfg->scf_orbital_type); + cfg->scf_orbital_type); Eigen::VectorXd start_kappa = Eigen::VectorXd::Zero(kappa_.size()); Eigen::VectorXd kappa_dir = kappa_; // Search direction From 584c4af310103bd4ac54d8ddd27fd9f28ca606b3 Mon Sep 17 00:00:00 2001 From: v-boqinzhang Date: Thu, 30 Apr 2026 11:36:03 -0700 Subject: [PATCH 19/29] resolve comments 3 --- .../microsoft/scf/src/scf_algorithm/gdm.cpp | 48 ++++++------------- 1 file changed, 15 insertions(+), 33 deletions(-) diff --git a/cpp/src/qdk/chemistry/algorithms/microsoft/scf/src/scf_algorithm/gdm.cpp b/cpp/src/qdk/chemistry/algorithms/microsoft/scf/src/scf_algorithm/gdm.cpp index 7abc906d4..2066f2921 100644 --- a/cpp/src/qdk/chemistry/algorithms/microsoft/scf/src/scf_algorithm/gdm.cpp +++ b/cpp/src/qdk/chemistry/algorithms/microsoft/scf/src/scf_algorithm/gdm.cpp @@ -239,27 +239,21 @@ static void compute_restricted_unrestricted_gradient( * @param[in] density_matrix Density matrix in AO basis * @param[in] num_electrons Occupied orbital counts (alpha, beta) * @param[in] rotation_size Rotation size for the ROHF kappa vector - * @param[in] num_molecular_orbitals Total molecular orbitals in the system - * @param[out] generalized_fock_mo Preallocated generalized Fock matrix in MO * @param[out] gradient Output gradient vector */ static void compute_restricted_open_shell_gradient( const SCFImpl& scf_impl, const RowMajorMatrix& C, const RowMajorMatrix& density_matrix, const std::vector& num_electrons, - const std::vector& rotation_size, int num_molecular_orbitals, - RowMajorMatrix& generalized_fock_mo, Eigen::VectorXd& gradient) { + const std::vector& rotation_size, Eigen::VectorXd& gradient) { const int total_rotation_size = rotation_size[0]; gradient.setZero(total_rotation_size); + const int num_molecular_orbitals = static_cast(C.cols()); const int num_closed_orbitals = num_electrons[1]; const int num_open_orbitals = num_electrons[0] - num_closed_orbitals; const int num_virtual_orbitals = num_molecular_orbitals - num_electrons[0]; - - if (generalized_fock_mo.rows() != num_molecular_orbitals || - generalized_fock_mo.cols() != num_molecular_orbitals) { - throw std::invalid_argument( - "generalized_fock_mo must be preallocated to MO dimensions."); - } + RowMajorMatrix generalized_fock_mo = + RowMajorMatrix::Zero(num_molecular_orbitals, num_molecular_orbitals); const int num_atomic_orbitals = scf_impl.get_num_atomic_orbitals(); const auto& H_ao_full = scf_impl.get_core_hamiltonian(); @@ -543,11 +537,9 @@ Eigen::VectorXd GDMLineFunctor::grad(const Eigen::VectorXd& x) { Eigen::VectorXd gradient; if (scf_orbital_type_ == SCFOrbitalType::RestrictedOpenShell) { - RowMajorMatrix generalized_fock_mo = - RowMajorMatrix::Zero(num_molecular_orbitals_, num_molecular_orbitals_); - compute_restricted_open_shell_gradient( - scf_impl_, cached_C_, cached_P_, num_electrons_, rotation_size_, - num_molecular_orbitals_, generalized_fock_mo, gradient); + compute_restricted_open_shell_gradient(scf_impl_, cached_C_, cached_P_, + num_electrons_, rotation_size_, + gradient); } else { compute_restricted_unrestricted_gradient( cached_F_, cached_C_, num_electrons_, rotation_offset_, rotation_size_, @@ -710,9 +702,6 @@ class GDM { /// Eigenvalues of pseudo-canonical orbitals, used for building Hessian Eigen::VectorXd pseudo_canonical_eigenvalues_; - /// Generalized Fock matrix in MO basis for ROHF - RowMajorMatrix generalized_fock_mo_; - Eigen::VectorXd kappa_; // vertical rotation matrix of this step /// Energy of the last accepted step, used to decide if we rescale the kappa @@ -753,12 +742,6 @@ GDM::GDM(const SCFContext& ctx, int history_size_limit) num_electrons_ = {num_alpha_electrons, num_beta_electrons}; history_size_ = 0; pseudo_canonical_eigenvalues_ = Eigen::VectorXd::Zero(num_molecular_orbitals); - if (cfg.scf_orbital_type == SCFOrbitalType::RestrictedOpenShell) { - generalized_fock_mo_ = - RowMajorMatrix::Zero(num_molecular_orbitals, num_molecular_orbitals); - } else { - generalized_fock_mo_ = RowMajorMatrix(0, 0); - } if (history_size_limit < 1) { throw std::invalid_argument( "GDM history size limit must be at least 1, got: " + @@ -948,6 +931,7 @@ void GDM::generate_restricted_unrestricted_pseudo_canonical_orbital_( Eigen::Block history_kappa_spin, Eigen::Block history_dgrad_spin, Eigen::VectorBlock current_gradient_spin) { + const int num_atomic_orbitals = C.rows(); const int num_molecular_orbitals = C.cols(); const int num_occupied_orbitals = num_electrons_[spin_index]; const int num_virtual_orbitals = @@ -967,12 +951,12 @@ void GDM::generate_restricted_unrestricted_pseudo_canonical_orbital_( RowMajorMatrix F_MO = RowMajorMatrix::Zero(num_molecular_orbitals, num_molecular_orbitals); - const int spin_block_offset = - num_molecular_orbitals * num_molecular_orbitals * spin_index; - const double* C_block_ptr = C.data() + spin_block_offset; - const double* F_block_ptr = F.data() + spin_block_offset; - std::vector atba_workspace( - static_cast(num_molecular_orbitals) * num_molecular_orbitals); + const double* C_block_ptr = + C.data() + num_atomic_orbitals * num_molecular_orbitals * spin_index; + const double* F_block_ptr = + F.data() + num_molecular_orbitals * num_molecular_orbitals * spin_index; + std::vector atba_workspace(static_cast(num_atomic_orbitals) * + num_molecular_orbitals); compute_atba_gemm(C_block_ptr, F_block_ptr, F_MO.data(), num_molecular_orbitals, num_molecular_orbitals, atba_workspace); @@ -1251,11 +1235,9 @@ void GDM::iterate(SCFImpl& scf_impl) { } if (cfg->scf_orbital_type == SCFOrbitalType::RestrictedOpenShell) { - generalized_fock_mo_.setZero(); compute_restricted_open_shell_gradient( scf_impl, C, scf_impl.get_density_matrix(), num_electrons_, - rotation_size_, num_molecular_orbitals, generalized_fock_mo_, - current_gradient_); + rotation_size_, current_gradient_); } else { compute_restricted_unrestricted_gradient( F, C, num_electrons_, rotation_offset_, rotation_size_, From 51494e8b6ab7207362d81ee1ec7e959e54ffd6b4 Mon Sep 17 00:00:00 2001 From: v-boqinzhang Date: Thu, 30 Apr 2026 12:30:34 -0700 Subject: [PATCH 20/29] small debug --- .../algorithms/microsoft/scf/src/scf_algorithm/gdm.cpp | 8 ++++++-- 1 file changed, 6 insertions(+), 2 deletions(-) diff --git a/cpp/src/qdk/chemistry/algorithms/microsoft/scf/src/scf_algorithm/gdm.cpp b/cpp/src/qdk/chemistry/algorithms/microsoft/scf/src/scf_algorithm/gdm.cpp index 2066f2921..acde645fd 100644 --- a/cpp/src/qdk/chemistry/algorithms/microsoft/scf/src/scf_algorithm/gdm.cpp +++ b/cpp/src/qdk/chemistry/algorithms/microsoft/scf/src/scf_algorithm/gdm.cpp @@ -620,6 +620,7 @@ class GDM { */ void generate_restricted_unrestricted_pseudo_canonical_orbital_( const RowMajorMatrix& F, RowMajorMatrix& C, const int spin_index, + const int num_orbital_spin_blocks, Eigen::Block history_kappa_spin, Eigen::Block history_dgrad_spin, Eigen::VectorBlock current_gradient_spin); @@ -928,10 +929,12 @@ static void rotate_gradient_to_pseudo_canonical_basis( void GDM::generate_restricted_unrestricted_pseudo_canonical_orbital_( const RowMajorMatrix& F, RowMajorMatrix& C, const int spin_index, + const int num_orbital_spin_blocks, Eigen::Block history_kappa_spin, Eigen::Block history_dgrad_spin, Eigen::VectorBlock current_gradient_spin) { - const int num_atomic_orbitals = C.rows(); + const int num_atomic_orbitals = + (num_orbital_spin_blocks == 2) ? C.rows() / 2 : C.rows(); const int num_molecular_orbitals = C.cols(); const int num_occupied_orbitals = num_electrons_[spin_index]; const int num_virtual_orbitals = @@ -1013,7 +1016,8 @@ void GDM::build_restricted_unrestricted_pseudo_canonical_orbitals_hessian_( // Generate pseudo-canonical orbitals and transform gradient and history generate_restricted_unrestricted_pseudo_canonical_orbital_( - F, C, i, history_kappa_spin, history_dgrad_spin, current_gradient_spin); + F, C, i, num_orbital_spin_blocks_, history_kappa_spin, + history_dgrad_spin, current_gradient_spin); // Build this spin's segment of initial Hessian // Reference: Helgaker, T., Jorgensen, P., & Olsen, J. (2000). Molecular From 2cd4678188d6ca12af4d0366c00b7c7121015b67 Mon Sep 17 00:00:00 2001 From: v-boqinzhang Date: Thu, 30 Apr 2026 13:59:51 -0700 Subject: [PATCH 21/29] add cached_J_ and cached_K_ for ROHF --- .../microsoft/scf/src/scf/ks_impl.cpp | 6 +- .../microsoft/scf/src/scf/ks_impl.h | 5 +- .../microsoft/scf/src/scf/scf_impl.cpp | 13 ++++- .../microsoft/scf/src/scf/scf_impl.h | 7 ++- .../microsoft/scf/src/scf_algorithm/gdm.cpp | 56 +++++++++++++------ .../scf/src/scf_algorithm/scf_algorithm.cpp | 2 - 6 files changed, 65 insertions(+), 24 deletions(-) diff --git a/cpp/src/qdk/chemistry/algorithms/microsoft/scf/src/scf/ks_impl.cpp b/cpp/src/qdk/chemistry/algorithms/microsoft/scf/src/scf/ks_impl.cpp index f1d43a6d5..adca940df 100644 --- a/cpp/src/qdk/chemistry/algorithms/microsoft/scf/src/scf/ks_impl.cpp +++ b/cpp/src/qdk/chemistry/algorithms/microsoft/scf/src/scf/ks_impl.cpp @@ -128,12 +128,14 @@ double KSImpl::total_energy_() { std::pair KSImpl::evaluate_trial_density_energy_and_fock( - const RowMajorMatrix& P_matrix, const std::source_location& loc) const { + const RowMajorMatrix& P_matrix, double* J_out, double* K_out, + const std::source_location& loc) const { QDK_LOG_TRACE_ENTERING(); // Fock matrix from base class does not include XC contributions; XC terms are // added below auto [total_energy, F_matrix] = - SCFImpl::evaluate_trial_density_energy_and_fock(P_matrix, loc); + SCFImpl::evaluate_trial_density_energy_and_fock(P_matrix, J_out, K_out, + loc); // Do not update XC_ here: XC_ is a member variable and must not be modified // in this const trial evaluation. double scf_xc_energy = 0.0; diff --git a/cpp/src/qdk/chemistry/algorithms/microsoft/scf/src/scf/ks_impl.h b/cpp/src/qdk/chemistry/algorithms/microsoft/scf/src/scf/ks_impl.h index 85264864a..97fca79b1 100644 --- a/cpp/src/qdk/chemistry/algorithms/microsoft/scf/src/scf/ks_impl.h +++ b/cpp/src/qdk/chemistry/algorithms/microsoft/scf/src/scf/ks_impl.h @@ -86,13 +86,16 @@ class KSImpl : public SCFImpl { * must be fully initialized before calling it. * * @param P_matrix Trial density matrix + * @param J_out Optional output buffer for Coulomb matrix J + * @param K_out Optional output buffer for exchange matrix K * @param loc Source location of the caller (automatically captured) * @return std::pair containing: * - first: total energy including XC in Hartree * - second: Fock matrix in AO basis */ std::pair evaluate_trial_density_energy_and_fock( - const RowMajorMatrix& P_matrix, + const RowMajorMatrix& P_matrix, double* J_out = nullptr, + double* K_out = nullptr, const std::source_location& loc = std::source_location::current()) const override; diff --git a/cpp/src/qdk/chemistry/algorithms/microsoft/scf/src/scf/scf_impl.cpp b/cpp/src/qdk/chemistry/algorithms/microsoft/scf/src/scf/scf_impl.cpp index e3d4a3bf3..b264acfc5 100644 --- a/cpp/src/qdk/chemistry/algorithms/microsoft/scf/src/scf/scf_impl.cpp +++ b/cpp/src/qdk/chemistry/algorithms/microsoft/scf/src/scf/scf_impl.cpp @@ -16,6 +16,7 @@ #endif #include +#include #include #include #include @@ -1035,7 +1036,8 @@ void SCFImpl::write_gradients_(const std::vector& gradients, std::pair SCFImpl::evaluate_trial_density_energy_and_fock( - const RowMajorMatrix& P_matrix, const std::source_location& loc) const { + const RowMajorMatrix& P_matrix, double* J_out, double* K_out, + const std::source_location& loc) const { QDK_LOG_TRACE_ENTERING(); QDK_LOGGER().debug( @@ -1056,6 +1058,15 @@ SCFImpl::evaluate_trial_density_energy_and_fock( num_density_matrices_ * num_atomic_orbitals_, num_atomic_orbitals_); eri_->build_JK(P_matrix.data(), J_matrix.data(), K_matrix.data(), alpha, beta, omega); + + if (J_out != nullptr) { + std::memcpy(J_out, J_matrix.data(), + sizeof(double) * static_cast(J_matrix.size())); + } + if (K_out != nullptr) { + std::memcpy(K_out, K_matrix.data(), + sizeof(double) * static_cast(K_matrix.size())); + } #ifdef QDK_CHEMISTRY_ENABLE_PCM throw std::runtime_error("PCM is not supported in trial density evaluation."); #endif diff --git a/cpp/src/qdk/chemistry/algorithms/microsoft/scf/src/scf/scf_impl.h b/cpp/src/qdk/chemistry/algorithms/microsoft/scf/src/scf/scf_impl.h index 6c6fb6650..5aa314cb7 100644 --- a/cpp/src/qdk/chemistry/algorithms/microsoft/scf/src/scf/scf_impl.h +++ b/cpp/src/qdk/chemistry/algorithms/microsoft/scf/src/scf/scf_impl.h @@ -221,6 +221,10 @@ class SCFImpl { * must be fully initialized before calling it. * * @param P_matrix Trial density matrix + * @param J_out Optional output buffer for Coulomb matrix J (same shape as + * P_matrix) + * @param K_out Optional output buffer for exchange matrix K (same shape as + * P_matrix) * @param loc Source location of the caller (automatically captured) * @return std::pair containing: * - first: total energy in Hartree @@ -228,7 +232,8 @@ class SCFImpl { */ virtual std::pair evaluate_trial_density_energy_and_fock( - const RowMajorMatrix& P_matrix, + const RowMajorMatrix& P_matrix, double* J_out = nullptr, + double* K_out = nullptr, const std::source_location& loc = std::source_location::current()) const; protected: diff --git a/cpp/src/qdk/chemistry/algorithms/microsoft/scf/src/scf_algorithm/gdm.cpp b/cpp/src/qdk/chemistry/algorithms/microsoft/scf/src/scf_algorithm/gdm.cpp index acde645fd..af0423f93 100644 --- a/cpp/src/qdk/chemistry/algorithms/microsoft/scf/src/scf_algorithm/gdm.cpp +++ b/cpp/src/qdk/chemistry/algorithms/microsoft/scf/src/scf_algorithm/gdm.cpp @@ -234,16 +234,18 @@ static void compute_restricted_unrestricted_gradient( * The gradient is packed as (iw, wa, ia) blocks following the same * segmentation used by the ROHF kappa vector. * - * @param[in] scf_impl SCF implementation for J/K construction + * @param[in] scf_impl SCF implementation for core Hamiltonian access * @param[in] C Molecular orbital coefficient matrix - * @param[in] density_matrix Density matrix in AO basis + * @param[in] J_ao Coulomb matrix in AO basis for alpha/beta blocks + * @param[in] K_ao Exchange matrix in AO basis for alpha/beta blocks * @param[in] num_electrons Occupied orbital counts (alpha, beta) * @param[in] rotation_size Rotation size for the ROHF kappa vector * @param[out] gradient Output gradient vector */ static void compute_restricted_open_shell_gradient( const SCFImpl& scf_impl, const RowMajorMatrix& C, - const RowMajorMatrix& density_matrix, const std::vector& num_electrons, + const RowMajorMatrix& J_ao, const RowMajorMatrix& K_ao, + const std::vector& num_electrons, const std::vector& rotation_size, Eigen::VectorXd& gradient) { const int total_rotation_size = rotation_size[0]; gradient.setZero(total_rotation_size); @@ -260,12 +262,6 @@ static void compute_restricted_open_shell_gradient( const RowMajorMatrix H_ao = H_ao_full.block(0, 0, num_atomic_orbitals, num_atomic_orbitals); - RowMajorMatrix J_ao = - RowMajorMatrix::Zero(2 * num_atomic_orbitals, num_atomic_orbitals); - RowMajorMatrix K_ao = - RowMajorMatrix::Zero(2 * num_atomic_orbitals, num_atomic_orbitals); - scf_impl.build_jk_matrices(density_matrix, J_ao, K_ao); - const auto J_alpha_ao = Eigen::Map( J_ao.data(), num_atomic_orbitals, num_atomic_orbitals); const auto J_beta_ao = Eigen::Map( @@ -393,6 +389,14 @@ class GDMLineFunctor { num_molecular_orbitals_(static_cast(C_pseudo_canonical.cols())), scf_orbital_type_(scf_orbital_type), cached_kappa_(Eigen::VectorXd()), + cached_J_(scf_orbital_type == SCFOrbitalType::RestrictedOpenShell + ? RowMajorMatrix::Zero(2 * num_atomic_orbitals_, + num_atomic_orbitals_) + : RowMajorMatrix()), + cached_K_(scf_orbital_type == SCFOrbitalType::RestrictedOpenShell + ? RowMajorMatrix::Zero(2 * num_atomic_orbitals_, + num_atomic_orbitals_) + : RowMajorMatrix()), cached_energy_(std::numeric_limits::infinity()) { if (scf_orbital_type == SCFOrbitalType::Unrestricted && (C_pseudo_canonical.rows() % 2 != 0)) { @@ -458,6 +462,8 @@ class GDMLineFunctor { Eigen::VectorXd cached_kappa_; // Cached kappa vector double cached_energy_; RowMajorMatrix cached_F_; // Needed for gradient computation + RowMajorMatrix cached_J_; // Needed for ROHF gradient computation + RowMajorMatrix cached_K_; // Needed for ROHF gradient computation RowMajorMatrix cached_C_; // For writing back to scf_impl RowMajorMatrix cached_P_; // For writing back to scf_impl }; @@ -517,8 +523,16 @@ double GDMLineFunctor::eval(const Eigen::VectorXd& x) { } // Evaluate energy and Fock matrix using trial density matrix - auto [energy, F_trial] = - scf_impl_.evaluate_trial_density_energy_and_fock(cached_P_); + double energy; + RowMajorMatrix F_trial; + if (scf_orbital_type_ == SCFOrbitalType::RestrictedOpenShell) { + std::tie(energy, F_trial) = + scf_impl_.evaluate_trial_density_energy_and_fock( + cached_P_, cached_J_.data(), cached_K_.data()); + } else { + std::tie(energy, F_trial) = + scf_impl_.evaluate_trial_density_energy_and_fock(cached_P_); + } // Cache all results for potential grad() call at same kappa cached_energy_ = energy; @@ -537,9 +551,9 @@ Eigen::VectorXd GDMLineFunctor::grad(const Eigen::VectorXd& x) { Eigen::VectorXd gradient; if (scf_orbital_type_ == SCFOrbitalType::RestrictedOpenShell) { - compute_restricted_open_shell_gradient(scf_impl_, cached_C_, cached_P_, - num_electrons_, rotation_size_, - gradient); + compute_restricted_open_shell_gradient(scf_impl_, cached_C_, cached_J_, + cached_K_, num_electrons_, + rotation_size_, gradient); } else { compute_restricted_unrestricted_gradient( cached_F_, cached_C_, num_electrons_, rotation_offset_, rotation_size_, @@ -610,6 +624,8 @@ class GDM { * @param[in] F Fock matrix in AO basis * @param[in,out] C Molecular orbital coefficient matrix * @param[in] spin_index Spin index (0 for alpha, 1 for beta) + * @param[in] num_orbital_spin_blocks Number of spin blocks (1 for RHF, 2 + * for UHF) * @param[in,out] history_kappa_spin Block reference to history kappa for this * spin * @param[in,out] history_dgrad_spin Block reference to history dgrad for this @@ -1239,9 +1255,15 @@ void GDM::iterate(SCFImpl& scf_impl) { } if (cfg->scf_orbital_type == SCFOrbitalType::RestrictedOpenShell) { - compute_restricted_open_shell_gradient( - scf_impl, C, scf_impl.get_density_matrix(), num_electrons_, - rotation_size_, current_gradient_); + const int num_atomic_orbitals = scf_impl.get_num_atomic_orbitals(); + RowMajorMatrix J_ao = + RowMajorMatrix::Zero(2 * num_atomic_orbitals, num_atomic_orbitals); + RowMajorMatrix K_ao = + RowMajorMatrix::Zero(2 * num_atomic_orbitals, num_atomic_orbitals); + scf_impl.build_jk_matrices(scf_impl.get_density_matrix(), J_ao, K_ao); + compute_restricted_open_shell_gradient(scf_impl, C, J_ao, K_ao, + num_electrons_, rotation_size_, + current_gradient_); } else { compute_restricted_unrestricted_gradient( F, C, num_electrons_, rotation_offset_, rotation_size_, diff --git a/cpp/src/qdk/chemistry/algorithms/microsoft/scf/src/scf_algorithm/scf_algorithm.cpp b/cpp/src/qdk/chemistry/algorithms/microsoft/scf/src/scf_algorithm/scf_algorithm.cpp index 1332e78ca..34caeeb10 100644 --- a/cpp/src/qdk/chemistry/algorithms/microsoft/scf/src/scf_algorithm/scf_algorithm.cpp +++ b/cpp/src/qdk/chemistry/algorithms/microsoft/scf/src/scf_algorithm/scf_algorithm.cpp @@ -245,8 +245,6 @@ void SCFAlgorithm::update_density_matrix(RowMajorMatrix& P, "Coefficient matrix rows do not match orbital set count"); } - // For ASAHF and ROHF, the density matrix construction is different and - // will be handled in the overridden methods const double occupancy_factor = unrestricted ? 1.0 : 2.0; for (int i = 0; i < num_orbital_sets; ++i) { const int n_occ = (i == 0) ? nelec_alpha : nelec_beta; From 71e312493387a57266cb164cc8391ce6a3b7a540 Mon Sep 17 00:00:00 2001 From: v-boqinzhang Date: Thu, 30 Apr 2026 18:40:07 -0700 Subject: [PATCH 22/29] resolve comments 4 --- .../microsoft/scf/src/scf_algorithm/gdm.cpp | 16 +++++++--------- 1 file changed, 7 insertions(+), 9 deletions(-) diff --git a/cpp/src/qdk/chemistry/algorithms/microsoft/scf/src/scf_algorithm/gdm.cpp b/cpp/src/qdk/chemistry/algorithms/microsoft/scf/src/scf_algorithm/gdm.cpp index af0423f93..c32698f46 100644 --- a/cpp/src/qdk/chemistry/algorithms/microsoft/scf/src/scf_algorithm/gdm.cpp +++ b/cpp/src/qdk/chemistry/algorithms/microsoft/scf/src/scf_algorithm/gdm.cpp @@ -184,8 +184,8 @@ static void compute_restricted_unrestricted_gradient( } gradient.setZero(total_rotation_size); - std::vector atba_workspace( - static_cast(num_molecular_orbitals) * num_molecular_orbitals); + std::vector atba_workspace(static_cast(num_atomic_orbitals) * + num_molecular_orbitals); for (int spin_index = 0; spin_index < num_orbital_spin_blocks; ++spin_index) { const int num_occupied_orbitals = num_electrons[spin_index]; @@ -259,8 +259,6 @@ static void compute_restricted_open_shell_gradient( const int num_atomic_orbitals = scf_impl.get_num_atomic_orbitals(); const auto& H_ao_full = scf_impl.get_core_hamiltonian(); - const RowMajorMatrix H_ao = - H_ao_full.block(0, 0, num_atomic_orbitals, num_atomic_orbitals); const auto J_alpha_ao = Eigen::Map( J_ao.data(), num_atomic_orbitals, num_atomic_orbitals); @@ -282,8 +280,9 @@ static void compute_restricted_open_shell_gradient( // Calculate Generalized Fock matrix in MO basis RowMajorMatrix H_mo = RowMajorMatrix::Zero(num_molecular_orbitals, num_molecular_orbitals); - compute_atba_gemm(C_ao_mo_ptr, H_ao.data(), H_mo.data(), num_atomic_orbitals, - num_molecular_orbitals, atba_workspace); + compute_atba_gemm(C_ao_mo_ptr, H_ao_full.data(), H_mo.data(), + num_atomic_orbitals, num_molecular_orbitals, + atba_workspace); RowMajorMatrix J_alpha_mo = RowMajorMatrix::Zero(num_molecular_orbitals, num_molecular_orbitals); @@ -976,9 +975,8 @@ void GDM::generate_restricted_unrestricted_pseudo_canonical_orbital_( F.data() + num_molecular_orbitals * num_molecular_orbitals * spin_index; std::vector atba_workspace(static_cast(num_atomic_orbitals) * num_molecular_orbitals); - compute_atba_gemm(C_block_ptr, F_block_ptr, F_MO.data(), - num_molecular_orbitals, num_molecular_orbitals, - atba_workspace); + compute_atba_gemm(C_block_ptr, F_block_ptr, F_MO.data(), num_atomic_orbitals, + num_molecular_orbitals, atba_workspace); // Perform pseudo-canonical transformation // Diagonalize occupied/virtual blocks and rotate orbitals to the From c88fc5cf6641c774e2e7cc8f86c7f995f5a04a9c Mon Sep 17 00:00:00 2001 From: v-boqinzhang Date: Fri, 1 May 2026 09:30:22 -0700 Subject: [PATCH 23/29] resolve comments 5 --- .../qdk/chemistry/scf/core/scf_algorithm.h | 31 --------------- .../microsoft/scf/src/scf_algorithm/gdm.cpp | 2 +- .../{line_search.h => scf_helper.h} | 39 ++++++++++++++++++- 3 files changed, 38 insertions(+), 34 deletions(-) rename cpp/src/qdk/chemistry/algorithms/microsoft/scf/src/scf_algorithm/{line_search.h => scf_helper.h} (79%) diff --git a/cpp/src/qdk/chemistry/algorithms/microsoft/scf/include/qdk/chemistry/scf/core/scf_algorithm.h b/cpp/src/qdk/chemistry/algorithms/microsoft/scf/include/qdk/chemistry/scf/core/scf_algorithm.h index a84d854ba..f4c053703 100644 --- a/cpp/src/qdk/chemistry/algorithms/microsoft/scf/include/qdk/chemistry/scf/core/scf_algorithm.h +++ b/cpp/src/qdk/chemistry/algorithms/microsoft/scf/include/qdk/chemistry/scf/core/scf_algorithm.h @@ -7,42 +7,11 @@ #include #include -#include #include #include -#include namespace qdk::chemistry::scf { -/** - * @brief Compute C = A^T * B * A using BLAS GEMM only - * - * Matrix shapes: - * A: m x n - * B: m x m - * C: n x n - * - * @param[in] m Row count of A and dimension of square B block. - * @param[in] n Column count of A and dimension of square C block. - * - * Assumed leading dimensions for contiguous storage: - * RowMajor: lda = n, ldb = m, ldc = n - * ColMajor: lda = m, ldb = m, ldc = n - * - * Pointers may reference sub-block starts (not necessarily matrix origin). - * The intermediate workspace is provided by caller to avoid repeated - * allocations across consecutive calls. - * - * @param[in] A Pointer to the A matrix block. - * @param[in] B Pointer to the B matrix block. - * @param[out] C Pointer to output C matrix block. - * @param[in,out] workspace Temporary buffer of at least m*n doubles. - * @param[in] layout BLAS matrix layout (RowMajor by default). - */ -void compute_atba_gemm(const double* A, const double* B, double* C, int m, - int n, std::vector& workspace, - blas::Layout layout = blas::Layout::RowMajor); - // Forward declaration class SCFImpl; /** diff --git a/cpp/src/qdk/chemistry/algorithms/microsoft/scf/src/scf_algorithm/gdm.cpp b/cpp/src/qdk/chemistry/algorithms/microsoft/scf/src/scf_algorithm/gdm.cpp index c32698f46..4be599b22 100644 --- a/cpp/src/qdk/chemistry/algorithms/microsoft/scf/src/scf_algorithm/gdm.cpp +++ b/cpp/src/qdk/chemistry/algorithms/microsoft/scf/src/scf_algorithm/gdm.cpp @@ -15,10 +15,10 @@ #include #include "../scf/scf_impl.h" -#include "line_search.h" #include "qdk/chemistry/scf/core/scf.h" #include "qdk/chemistry/scf/core/scf_algorithm.h" #include "qdk/chemistry/scf/core/types.h" +#include "scf_helper.h" #include "util/matrix_exp.h" #ifdef QDK_CHEMISTRY_ENABLE_GPU #include "util/gpu/cuda_helper.h" diff --git a/cpp/src/qdk/chemistry/algorithms/microsoft/scf/src/scf_algorithm/line_search.h b/cpp/src/qdk/chemistry/algorithms/microsoft/scf/src/scf_algorithm/scf_helper.h similarity index 79% rename from cpp/src/qdk/chemistry/algorithms/microsoft/scf/src/scf_algorithm/line_search.h rename to cpp/src/qdk/chemistry/algorithms/microsoft/scf/src/scf_algorithm/scf_helper.h index 1824b96e7..5bae7a838 100644 --- a/cpp/src/qdk/chemistry/algorithms/microsoft/scf/src/scf_algorithm/line_search.h +++ b/cpp/src/qdk/chemistry/algorithms/microsoft/scf/src/scf_algorithm/scf_helper.h @@ -5,11 +5,44 @@ #pragma once #include +#include #include #include #include +#include -namespace qdk::chemistry::scf::impl { +namespace qdk::chemistry::scf { + +/** + * @brief Compute C = A^T * B * A using BLAS GEMM only + * + * Matrix shapes: + * A: m x n + * B: m x m + * C: n x n + * + * @param[in] m Row count of A and dimension of square B block. + * @param[in] n Column count of A and dimension of square C block. + * + * Assumed leading dimensions for contiguous storage: + * RowMajor: lda = n, ldb = m, ldc = n + * ColMajor: lda = m, ldb = m, ldc = n + * + * Pointers may reference sub-block starts (not necessarily matrix origin). + * The intermediate workspace is provided by caller to avoid repeated + * allocations across consecutive calls. + * + * @param[in] A Pointer to the A matrix block. + * @param[in] B Pointer to the B matrix block. + * @param[out] C Pointer to output C matrix block. + * @param[in,out] workspace Temporary buffer of at least m*n doubles. + * @param[in] layout BLAS matrix layout (RowMajor by default). + */ +void compute_atba_gemm(const double* A, const double* B, double* C, int m, + int n, std::vector& workspace, + blas::Layout layout = blas::Layout::RowMajor); + +namespace impl { /** * @brief Nocedal-Wright line search with strong Wolfe conditions. @@ -153,4 +186,6 @@ void nocedal_wright_line_search(Functor& op, if (!converged) throw std::runtime_error("Line Search Failed"); } -} // namespace qdk::chemistry::scf::impl +} // namespace impl + +} // namespace qdk::chemistry::scf From ebd54dc61a5c6ef0b565895c26b5c8bd2d8ccc44 Mon Sep 17 00:00:00 2001 From: v-boqinzhang Date: Fri, 1 May 2026 10:12:23 -0700 Subject: [PATCH 24/29] resolve comments 6 --- .../qdk/chemistry/scf/core/scf_algorithm.h | 16 ++++++------ .../microsoft/scf/src/scf_algorithm/diis.cpp | 8 ++---- .../microsoft/scf/src/scf_algorithm/gdm.cpp | 10 ++++---- .../scf/src/scf_algorithm/scf_algorithm.cpp | 25 +++++++++---------- 4 files changed, 26 insertions(+), 33 deletions(-) diff --git a/cpp/src/qdk/chemistry/algorithms/microsoft/scf/include/qdk/chemistry/scf/core/scf_algorithm.h b/cpp/src/qdk/chemistry/algorithms/microsoft/scf/include/qdk/chemistry/scf/core/scf_algorithm.h index f4c053703..9465da8bc 100644 --- a/cpp/src/qdk/chemistry/algorithms/microsoft/scf/include/qdk/chemistry/scf/core/scf_algorithm.h +++ b/cpp/src/qdk/chemistry/algorithms/microsoft/scf/include/qdk/chemistry/scf/core/scf_algorithm.h @@ -116,21 +116,19 @@ class SCFAlgorithm { int nelec_beta); /** - * @brief Try to provide ROHF convergence matrices for OG evaluation + * @brief Provide ROHF convergence matrices for OG evaluation * * For ROHF, some algorithms evaluate convergence using an effective Fock * matrix and total density matrix rather than the spin-blocked SCFImpl - * matrices. Derived algorithms can override this hook and provide pointers - * to those matrices. + * matrices. Derived algorithms can override this hook and provide + * references to those matrices. * * @param[in] scf_impl SCF implementation object - * @param[out] fock_matrix Pointer to convergence Fock matrix - * @param[out] density_matrix Pointer to convergence density matrix - * @return true if matrices were provided by the algorithm, false otherwise + * @return Pair of references to convergence Fock and density matrices + * @throws std::logic_error If ROHF convergence matrices are unavailable */ - virtual bool try_get_rohf_convergence_matrices( - const SCFImpl& scf_impl, const RowMajorMatrix*& fock_matrix, - const RowMajorMatrix*& density_matrix); + virtual std::pair + try_get_rohf_convergence_matrices(const SCFImpl& scf_impl); /** * @brief Calculate orbital gradient (OG) error for convergence checking diff --git a/cpp/src/qdk/chemistry/algorithms/microsoft/scf/src/scf_algorithm/diis.cpp b/cpp/src/qdk/chemistry/algorithms/microsoft/scf/src/scf_algorithm/diis.cpp index e10466388..d643fe65f 100644 --- a/cpp/src/qdk/chemistry/algorithms/microsoft/scf/src/scf_algorithm/diis.cpp +++ b/cpp/src/qdk/chemistry/algorithms/microsoft/scf/src/scf_algorithm/diis.cpp @@ -329,12 +329,8 @@ std::pair DIIS::select_working_matrices( SCFImpl& scf_impl) { QDK_LOG_TRACE_ENTERING(); if (ctx_.cfg->scf_orbital_type == SCFOrbitalType::RestrictedOpenShell) { - const RowMajorMatrix* rohf_fock = nullptr; - const RowMajorMatrix* rohf_density = nullptr; - if (!try_get_rohf_convergence_matrices(scf_impl, rohf_fock, rohf_density)) { - throw std::logic_error("ROHF convergence matrices are unavailable"); - } - return {rohf_convergence_density_matrix(), *rohf_fock}; + const auto rohf_matrices = try_get_rohf_convergence_matrices(scf_impl); + return {rohf_convergence_density_matrix(), rohf_matrices.first}; } return {scf_impl.density_matrix(), scf_impl.get_fock_matrix()}; } diff --git a/cpp/src/qdk/chemistry/algorithms/microsoft/scf/src/scf_algorithm/gdm.cpp b/cpp/src/qdk/chemistry/algorithms/microsoft/scf/src/scf_algorithm/gdm.cpp index 4be599b22..465a3e65f 100644 --- a/cpp/src/qdk/chemistry/algorithms/microsoft/scf/src/scf_algorithm/gdm.cpp +++ b/cpp/src/qdk/chemistry/algorithms/microsoft/scf/src/scf_algorithm/gdm.cpp @@ -972,7 +972,7 @@ void GDM::generate_restricted_unrestricted_pseudo_canonical_orbital_( const double* C_block_ptr = C.data() + num_atomic_orbitals * num_molecular_orbitals * spin_index; const double* F_block_ptr = - F.data() + num_molecular_orbitals * num_molecular_orbitals * spin_index; + F.data() + num_atomic_orbitals * num_atomic_orbitals * spin_index; std::vector atba_workspace(static_cast(num_atomic_orbitals) * num_molecular_orbitals); compute_atba_gemm(C_block_ptr, F_block_ptr, F_MO.data(), num_atomic_orbitals, @@ -983,8 +983,8 @@ void GDM::generate_restricted_unrestricted_pseudo_canonical_orbital_( // pseudo-canonical basis. RowMajorMatrix Uii = F_MO.block(0, 0, num_occupied_orbitals, num_occupied_orbitals); - auto C_occ_view = C.block(num_molecular_orbitals * spin_index, 0, - num_molecular_orbitals, num_occupied_orbitals); + auto C_occ_view = C.block(num_atomic_orbitals * spin_index, 0, + num_atomic_orbitals, num_occupied_orbitals); calculate_pseudo_canonical_orbital_block( Uii, num_occupied_orbitals, pseudo_canonical_eigenvalues_, 0, C_occ_view, num_molecular_orbitals, num_molecular_orbitals); @@ -992,8 +992,8 @@ void GDM::generate_restricted_unrestricted_pseudo_canonical_orbital_( RowMajorMatrix Uaa = F_MO.block(num_occupied_orbitals, num_occupied_orbitals, num_virtual_orbitals, num_virtual_orbitals); auto C_virt_view = - C.block(num_molecular_orbitals * spin_index, num_occupied_orbitals, - num_molecular_orbitals, num_virtual_orbitals); + C.block(num_atomic_orbitals * spin_index, num_occupied_orbitals, + num_atomic_orbitals, num_virtual_orbitals); calculate_pseudo_canonical_orbital_block( Uaa, num_virtual_orbitals, pseudo_canonical_eigenvalues_, num_occupied_orbitals, C_virt_view, num_molecular_orbitals, diff --git a/cpp/src/qdk/chemistry/algorithms/microsoft/scf/src/scf_algorithm/scf_algorithm.cpp b/cpp/src/qdk/chemistry/algorithms/microsoft/scf/src/scf_algorithm/scf_algorithm.cpp index 34caeeb10..0ecdfb2c6 100644 --- a/cpp/src/qdk/chemistry/algorithms/microsoft/scf/src/scf_algorithm/scf_algorithm.cpp +++ b/cpp/src/qdk/chemistry/algorithms/microsoft/scf/src/scf_algorithm/scf_algorithm.cpp @@ -261,12 +261,13 @@ void SCFAlgorithm::update_density_matrix(RowMajorMatrix& P, } } -bool SCFAlgorithm::try_get_rohf_convergence_matrices( - const SCFImpl& scf_impl, const RowMajorMatrix*& fock_matrix, - const RowMajorMatrix*& density_matrix) { +std::pair +SCFAlgorithm::try_get_rohf_convergence_matrices(const SCFImpl& scf_impl) { QDK_LOG_TRACE_ENTERING(); if (ctx_.cfg->scf_orbital_type != SCFOrbitalType::RestrictedOpenShell) { - return false; + throw std::logic_error( + "ROHF convergence matrices are only available for " + "RestrictedOpenShell calculations"); } const auto nelec_vec = scf_impl.get_num_electrons(); @@ -275,9 +276,8 @@ bool SCFAlgorithm::try_get_rohf_convergence_matrices( scf_impl.get_density_matrix(), nelec_vec[0], nelec_vec[1], rohf_effective_fock_, rohf_total_density_); - fock_matrix = &get_rohf_convergence_fock_matrix(); - density_matrix = &get_rohf_convergence_density_matrix(); - return true; + return {get_rohf_convergence_fock_matrix(), + get_rohf_convergence_density_matrix()}; } void SCFAlgorithm::build_rohf_f_p_matrix(const RowMajorMatrix& F, @@ -479,14 +479,13 @@ bool SCFAlgorithm::check_convergence(const SCFImpl& scf_impl) { RowMajorMatrix error_matrix; int num_orbital_sets = scf_impl.get_num_orbital_spin_blocks(); - const RowMajorMatrix* F_ptr; - const RowMajorMatrix* P_ptr; + const RowMajorMatrix* F_ptr = nullptr; + const RowMajorMatrix* P_ptr = nullptr; if (ctx_.cfg->scf_orbital_type == SCFOrbitalType::RestrictedOpenShell) { - if (!try_get_rohf_convergence_matrices(scf_impl, F_ptr, P_ptr)) { - throw std::logic_error( - "ROHF convergence matrices are not provided by this SCF algorithm"); - } + const auto rohf_matrices = try_get_rohf_convergence_matrices(scf_impl); + F_ptr = &rohf_matrices.first; + P_ptr = &rohf_matrices.second; } else { F_ptr = &scf_impl.get_fock_matrix(); P_ptr = &scf_impl.get_density_matrix(); From 7ffd6be546e70cdf0bea53e82dc1fa989de42816 Mon Sep 17 00:00:00 2001 From: v-boqinzhang Date: Mon, 18 May 2026 14:19:13 -0700 Subject: [PATCH 25/29] resolve comments 7 --- .../microsoft/scf/src/scf/ks_impl.cpp | 5 +++++ .../microsoft/scf/src/scf/scf_impl.cpp | 17 +++++++++++++---- .../microsoft/scf/src/scf_algorithm/gdm.cpp | 4 ++-- 3 files changed, 20 insertions(+), 6 deletions(-) diff --git a/cpp/src/qdk/chemistry/algorithms/microsoft/scf/src/scf/ks_impl.cpp b/cpp/src/qdk/chemistry/algorithms/microsoft/scf/src/scf/ks_impl.cpp index adca940df..a399a64f3 100644 --- a/cpp/src/qdk/chemistry/algorithms/microsoft/scf/src/scf/ks_impl.cpp +++ b/cpp/src/qdk/chemistry/algorithms/microsoft/scf/src/scf/ks_impl.cpp @@ -36,6 +36,11 @@ KSImpl::KSImpl(std::shared_ptr mol, const SCFConfig& cfg, std::shared_ptr raw_basis_set) : SCFImpl(mol, cfg, basis_set, raw_basis_set, true) { QDK_LOG_TRACE_ENTERING(); + if (cfg.scf_orbital_type == SCFOrbitalType::RestrictedOpenShell) { + throw std::invalid_argument( + "ROKS (Restricted Open-Shell Kohn-Sham) is not supported. " + "Use ROHF (scf_type=\"restricted\") without DFT functionals."); + } #ifdef ENABLE_NVTX3 NVTX3_FUNC_RANGE(); #endif diff --git a/cpp/src/qdk/chemistry/algorithms/microsoft/scf/src/scf/scf_impl.cpp b/cpp/src/qdk/chemistry/algorithms/microsoft/scf/src/scf/scf_impl.cpp index b264acfc5..f0287322e 100644 --- a/cpp/src/qdk/chemistry/algorithms/microsoft/scf/src/scf/scf_impl.cpp +++ b/cpp/src/qdk/chemistry/algorithms/microsoft/scf/src/scf/scf_impl.cpp @@ -619,10 +619,19 @@ void SCFImpl::iterate_() { auto [alpha, beta, omega] = get_hyb_coeff_(); eri_->build_JK(P_.data(), J_.data(), K_.data(), alpha, beta, omega); update_fock_(); - for (int i = 0; i < num_orbital_spin_blocks_; ++i) { - scf_algorithm_->solve_fock_eigenproblem(F_, S_, X_, C_, eigenvalues_, P_, - nelec_, num_atomic_orbitals_, - num_molecular_orbitals_, i); + + if (ctx_.cfg->scf_orbital_type == SCFOrbitalType::RestrictedOpenShell) { + const auto& [effective_fock, total_density] = + scf_algorithm_->try_get_rohf_convergence_matrices(*this); + scf_algorithm_->solve_fock_eigenproblem( + effective_fock, S_, X_, C_, eigenvalues_, P_, nelec_, + num_atomic_orbitals_, num_molecular_orbitals_, 0); + } else { + for (int i = 0; i < num_orbital_spin_blocks_; ++i) { + scf_algorithm_->solve_fock_eigenproblem( + F_, S_, X_, C_, eigenvalues_, P_, nelec_, num_atomic_orbitals_, + num_molecular_orbitals_, i); + } } scf_algorithm_->update_density_matrix( P_, C_, ctx_.cfg->scf_orbital_type == SCFOrbitalType::Unrestricted, diff --git a/cpp/src/qdk/chemistry/algorithms/microsoft/scf/src/scf_algorithm/gdm.cpp b/cpp/src/qdk/chemistry/algorithms/microsoft/scf/src/scf_algorithm/gdm.cpp index 465a3e65f..530766052 100644 --- a/cpp/src/qdk/chemistry/algorithms/microsoft/scf/src/scf_algorithm/gdm.cpp +++ b/cpp/src/qdk/chemistry/algorithms/microsoft/scf/src/scf_algorithm/gdm.cpp @@ -987,7 +987,7 @@ void GDM::generate_restricted_unrestricted_pseudo_canonical_orbital_( num_atomic_orbitals, num_occupied_orbitals); calculate_pseudo_canonical_orbital_block( Uii, num_occupied_orbitals, pseudo_canonical_eigenvalues_, 0, C_occ_view, - num_molecular_orbitals, num_molecular_orbitals); + num_atomic_orbitals, num_molecular_orbitals); RowMajorMatrix Uaa = F_MO.block(num_occupied_orbitals, num_occupied_orbitals, num_virtual_orbitals, num_virtual_orbitals); @@ -996,7 +996,7 @@ void GDM::generate_restricted_unrestricted_pseudo_canonical_orbital_( num_atomic_orbitals, num_virtual_orbitals); calculate_pseudo_canonical_orbital_block( Uaa, num_virtual_orbitals, pseudo_canonical_eigenvalues_, - num_occupied_orbitals, C_virt_view, num_molecular_orbitals, + num_occupied_orbitals, C_virt_view, num_atomic_orbitals, num_molecular_orbitals); // Transform the vectors in history_kappa, history_dgrad, and From 602ec72b4b0747442f977cf0d83c601ad44af827 Mon Sep 17 00:00:00 2001 From: v-boqinzhang Date: Tue, 19 May 2026 09:57:52 -0700 Subject: [PATCH 26/29] forbid trial_density_energy_fock with MPI --- .../chemistry/algorithms/microsoft/scf/src/scf/scf_impl.cpp | 6 ++++++ 1 file changed, 6 insertions(+) diff --git a/cpp/src/qdk/chemistry/algorithms/microsoft/scf/src/scf/scf_impl.cpp b/cpp/src/qdk/chemistry/algorithms/microsoft/scf/src/scf/scf_impl.cpp index f0287322e..8103b91c9 100644 --- a/cpp/src/qdk/chemistry/algorithms/microsoft/scf/src/scf/scf_impl.cpp +++ b/cpp/src/qdk/chemistry/algorithms/microsoft/scf/src/scf/scf_impl.cpp @@ -1049,6 +1049,12 @@ SCFImpl::evaluate_trial_density_energy_and_fock( const std::source_location& loc) const { QDK_LOG_TRACE_ENTERING(); + if (ctx_.cfg->mpi.world_size > 1) { + throw std::runtime_error( + "Temporary limitation: evaluate_trial_density_energy_and_fock is not " + "supported with MPI world_size > 1."); + } + QDK_LOGGER().debug( "Computing energy and Fock matrix by trial density matrix (called from " "function '{}' at line {})", From d004d96f4a4d60f818e32da7b5bdd276fb4c4daf Mon Sep 17 00:00:00 2001 From: v-boqinzhang Date: Tue, 19 May 2026 10:39:11 -0700 Subject: [PATCH 27/29] forbid GDM with MPI --- .../algorithms/microsoft/scf/src/scf_algorithm/gdm.cpp | 5 +++++ 1 file changed, 5 insertions(+) diff --git a/cpp/src/qdk/chemistry/algorithms/microsoft/scf/src/scf_algorithm/gdm.cpp b/cpp/src/qdk/chemistry/algorithms/microsoft/scf/src/scf_algorithm/gdm.cpp index 530766052..7e2cc2de7 100644 --- a/cpp/src/qdk/chemistry/algorithms/microsoft/scf/src/scf_algorithm/gdm.cpp +++ b/cpp/src/qdk/chemistry/algorithms/microsoft/scf/src/scf_algorithm/gdm.cpp @@ -1240,6 +1240,11 @@ void GDM::iterate(SCFImpl& scf_impl) { const auto& F = scf_impl.get_fock_matrix(); const auto* cfg = ctx_.cfg; + if (cfg->mpi.world_size > 1) { + throw std::runtime_error( + "Temporary limitation: evaluate_trial_density_energy_and_fock is not " + "supported with MPI world_size > 1."); + } const int num_molecular_orbitals = static_cast(ctx_.num_molecular_orbitals); From 32b6eee8d700571dbffc96410be41dab2beafc0a Mon Sep 17 00:00:00 2001 From: v-boqinzhang Date: Thu, 21 May 2026 10:34:02 -0700 Subject: [PATCH 28/29] resolve comments 8 --- .../qdk/chemistry/scf/core/scf_algorithm.h | 5 ++--- .../algorithms/microsoft/scf/src/scf/scf_impl.cpp | 15 ++++++++++++++- 2 files changed, 16 insertions(+), 4 deletions(-) diff --git a/cpp/src/qdk/chemistry/algorithms/microsoft/scf/include/qdk/chemistry/scf/core/scf_algorithm.h b/cpp/src/qdk/chemistry/algorithms/microsoft/scf/include/qdk/chemistry/scf/core/scf_algorithm.h index 9465da8bc..4abb888d0 100644 --- a/cpp/src/qdk/chemistry/algorithms/microsoft/scf/include/qdk/chemistry/scf/core/scf_algorithm.h +++ b/cpp/src/qdk/chemistry/algorithms/microsoft/scf/include/qdk/chemistry/scf/core/scf_algorithm.h @@ -120,14 +120,13 @@ class SCFAlgorithm { * * For ROHF, some algorithms evaluate convergence using an effective Fock * matrix and total density matrix rather than the spin-blocked SCFImpl - * matrices. Derived algorithms can override this hook and provide - * references to those matrices. + * matrices. This helper returns references to those matrices. * * @param[in] scf_impl SCF implementation object * @return Pair of references to convergence Fock and density matrices * @throws std::logic_error If ROHF convergence matrices are unavailable */ - virtual std::pair + std::pair try_get_rohf_convergence_matrices(const SCFImpl& scf_impl); /** diff --git a/cpp/src/qdk/chemistry/algorithms/microsoft/scf/src/scf/scf_impl.cpp b/cpp/src/qdk/chemistry/algorithms/microsoft/scf/src/scf/scf_impl.cpp index 8103b91c9..fd23bf3c3 100644 --- a/cpp/src/qdk/chemistry/algorithms/microsoft/scf/src/scf/scf_impl.cpp +++ b/cpp/src/qdk/chemistry/algorithms/microsoft/scf/src/scf/scf_impl.cpp @@ -621,8 +621,9 @@ void SCFImpl::iterate_() { update_fock_(); if (ctx_.cfg->scf_orbital_type == SCFOrbitalType::RestrictedOpenShell) { - const auto& [effective_fock, total_density] = + const auto rohf_convergence_matrices = scf_algorithm_->try_get_rohf_convergence_matrices(*this); + const auto& effective_fock = std::get<0>(rohf_convergence_matrices); scf_algorithm_->solve_fock_eigenproblem( effective_fock, S_, X_, C_, eigenvalues_, P_, nelec_, num_atomic_orbitals_, num_molecular_orbitals_, 0); @@ -1127,6 +1128,18 @@ SCFImpl::evaluate_trial_density_energy_and_fock( void SCFImpl::build_jk_matrices(const RowMajorMatrix& density_matrix, RowMajorMatrix& J, RowMajorMatrix& K) const { QDK_LOG_TRACE_ENTERING(); + + const int expected_rows = num_density_matrices_ * num_atomic_orbitals_; + const int expected_cols = num_atomic_orbitals_; + VERIFY_INPUT(density_matrix.rows() == expected_rows && + density_matrix.cols() == expected_cols, + "density_matrix shape should be (num_density_matrices_ * " + "num_atomic_orbitals_, num_atomic_orbitals_)"); + VERIFY_INPUT(J.rows() == expected_rows && J.cols() == expected_cols, + "J shape should match density_matrix shape"); + VERIFY_INPUT(K.rows() == expected_rows && K.cols() == expected_cols, + "K shape should match density_matrix shape"); + auto [alpha, beta, omega] = get_hyb_coeff_(); eri_->build_JK(density_matrix.data(), J.data(), K.data(), alpha, beta, omega); } From b85588e213e236e45bbab3f88c8dc93b70459107 Mon Sep 17 00:00:00 2001 From: v-boqinzhang Date: Thu, 21 May 2026 11:55:18 -0700 Subject: [PATCH 29/29] CLI review --- .../algorithms/microsoft/scf/src/scf/scf_impl.cpp | 2 ++ .../microsoft/scf/src/scf_algorithm/gdm.cpp | 14 +++++++++++++- .../scf/src/scf_algorithm/scf_algorithm.cpp | 4 ++-- cpp/tests/test_scf.cpp | 14 ++++++++++++++ python/tests/test_scf.py | 15 +++++++++++++++ 5 files changed, 46 insertions(+), 3 deletions(-) diff --git a/cpp/src/qdk/chemistry/algorithms/microsoft/scf/src/scf/scf_impl.cpp b/cpp/src/qdk/chemistry/algorithms/microsoft/scf/src/scf/scf_impl.cpp index fd23bf3c3..2d0b5efda 100644 --- a/cpp/src/qdk/chemistry/algorithms/microsoft/scf/src/scf/scf_impl.cpp +++ b/cpp/src/qdk/chemistry/algorithms/microsoft/scf/src/scf/scf_impl.cpp @@ -1050,6 +1050,8 @@ SCFImpl::evaluate_trial_density_energy_and_fock( const std::source_location& loc) const { QDK_LOG_TRACE_ENTERING(); + // TODO: This function is called within iterate_ which runs only on rank 0, + // but build_jk_matrices requires participation from all MPI ranks. if (ctx_.cfg->mpi.world_size > 1) { throw std::runtime_error( "Temporary limitation: evaluate_trial_density_energy_and_fock is not " diff --git a/cpp/src/qdk/chemistry/algorithms/microsoft/scf/src/scf_algorithm/gdm.cpp b/cpp/src/qdk/chemistry/algorithms/microsoft/scf/src/scf_algorithm/gdm.cpp index 7e2cc2de7..338d46c76 100644 --- a/cpp/src/qdk/chemistry/algorithms/microsoft/scf/src/scf_algorithm/gdm.cpp +++ b/cpp/src/qdk/chemistry/algorithms/microsoft/scf/src/scf_algorithm/gdm.cpp @@ -138,7 +138,9 @@ static void apply_restricted_open_shell_orbital_rotation( kappa_complete.block(virtual_start, 0, num_virtual_orbitals, num_closed_orbitals) = -kappa_ia.transpose(); - // 0.5 is for consistency with the gradient definition + // kappa_vector stores only unique (upper-triangular) rotation parameters. + // kappa_complete is the full antisymmetric matrix (K = -K^T), so each unique + // parameter appears twice. Scale by 0.5 to match the gradient convention. kappa_complete *= 0.5; RowMajorMatrix exp_kappa = @@ -258,6 +260,14 @@ static void compute_restricted_open_shell_gradient( RowMajorMatrix::Zero(num_molecular_orbitals, num_molecular_orbitals); const int num_atomic_orbitals = scf_impl.get_num_atomic_orbitals(); + + if (J_ao.size() < 2 * num_atomic_orbitals * num_atomic_orbitals || + K_ao.size() < 2 * num_atomic_orbitals * num_atomic_orbitals) { + throw std::invalid_argument( + "J_ao and K_ao must each contain alpha and beta blocks " + "(size >= 2 * num_atomic_orbitals^2)."); + } + const auto& H_ao_full = scf_impl.get_core_hamiltonian(); const auto J_alpha_ao = Eigen::Map( @@ -1240,6 +1250,8 @@ void GDM::iterate(SCFImpl& scf_impl) { const auto& F = scf_impl.get_fock_matrix(); const auto* cfg = ctx_.cfg; + // TODO: This function is called within iterate_ which runs only on rank 0, + // but build_jk_matrices requires participation from all MPI ranks. if (cfg->mpi.world_size > 1) { throw std::runtime_error( "Temporary limitation: evaluate_trial_density_energy_and_fock is not " diff --git a/cpp/src/qdk/chemistry/algorithms/microsoft/scf/src/scf_algorithm/scf_algorithm.cpp b/cpp/src/qdk/chemistry/algorithms/microsoft/scf/src/scf_algorithm/scf_algorithm.cpp index 0ecdfb2c6..594ca16e3 100644 --- a/cpp/src/qdk/chemistry/algorithms/microsoft/scf/src/scf_algorithm/scf_algorithm.cpp +++ b/cpp/src/qdk/chemistry/algorithms/microsoft/scf/src/scf_algorithm/scf_algorithm.cpp @@ -105,7 +105,7 @@ std::shared_ptr SCFAlgorithm::create(const SCFContext& ctx) { switch (cfg.scf_algorithm.method) { case SCFAlgorithmName::ASAHF: if (rohf_enabled) { - throw std::runtime_error("ROHF-enabled ASAHF is not supported!"); + throw std::invalid_argument("ROHF-enabled ASAHF is not supported!"); } return std::make_shared( ctx, cfg.scf_algorithm.diis_subspace_size); @@ -290,7 +290,7 @@ void SCFAlgorithm::build_rohf_f_p_matrix(const RowMajorMatrix& F, const int num_atomic_orbitals = static_cast(C.rows()); const int num_molecular_orbitals = static_cast(C.cols()); if (num_atomic_orbitals != num_molecular_orbitals) { - throw std::runtime_error( + throw std::invalid_argument( "ROHF build requires number of atomic orbitals to equal number of " "molecular orbitals!"); } diff --git a/cpp/tests/test_scf.cpp b/cpp/tests/test_scf.cpp index af51fb94c..751db3037 100644 --- a/cpp/tests/test_scf.cpp +++ b/cpp/tests/test_scf.cpp @@ -243,6 +243,20 @@ TEST_F(ScfTest, OH_ROHF_GDM) { EXPECT_TRUE(wfn_doublet->get_orbitals()->is_restricted()); } +TEST_F(ScfTest, Oxygen_atom_ROHF_GDM) { + auto oxygen = testing::create_oxygen_structure(); + auto scf_solver = ScfSolverFactory::create(); + scf_solver->settings().set("enable_gdm", true); + scf_solver->settings().set("method", "hf"); + scf_solver->settings().set("scf_type", "restricted"); + auto [E_triplet, wfn_triplet] = scf_solver->run(oxygen, 0, 3, "cc-pvdz"); + + EXPECT_NEAR(E_triplet, -74.787513074624, testing::scf_energy_tolerance); + + // Check triplet orbitals are restricted (ROHF) + EXPECT_TRUE(wfn_triplet->get_orbitals()->is_restricted()); +} + TEST_F(ScfTest, Oxygen_atom_invalid_energy_thresh_diis_switch_gdm) { auto oxygen = testing::create_oxygen_structure(); auto scf_solver = ScfSolverFactory::create(); diff --git a/python/tests/test_scf.py b/python/tests/test_scf.py index a52971309..3882ef5bb 100644 --- a/python/tests/test_scf.py +++ b/python/tests/test_scf.py @@ -304,6 +304,21 @@ def test_scf_solver_oh_rohf_gdm(self): assert abs(energy - (-74.361530753176)) < scf_energy_tolerance assert orbitals.is_restricted() + def test_scf_solver_oxygen_atom_rohf_gdm(self): + """Test SCF solver on oxygen atom triplet with ROHF/cc-pvdz and GDM.""" + oxygen = create_oxygen_structure() + scf_solver = algorithms.create("scf_solver") + + scf_solver.settings().set("enable_gdm", True) + scf_solver.settings().set("method", "hf") + scf_solver.settings().set("scf_type", "restricted") + + energy, wavefunction = scf_solver.run(oxygen, 0, 3, "cc-pvdz") + orbitals = wavefunction.get_orbitals() + + assert abs(energy - (-74.787513074624)) < scf_energy_tolerance + assert orbitals.is_restricted() + def test_scf_solver_oh_roks_invalid(self): """Test restricted open-shell KS request on OH doublet raises error.""" oh_structure = create_oh_structure()