diff --git a/alm/accelerate_solver.cpp b/alm/accelerate_solver.cpp index 2ec8dd6c..f2a86009 100644 --- a/alm/accelerate_solver.cpp +++ b/alm/accelerate_solver.cpp @@ -1,8 +1,4 @@ -// accelerate_solver.cpp -// -// Translation unit isolating (and thus ) from the -// hand-rolled Fortran BLAS/LAPACK prototypes in include/blas_wrapper.h / include/lapack_wrapper.h. -// See accelerate_solver.h for the rationale. This file intentionally includes NEITHER wrapper. +// Isolate Accelerate headers from the BLAS/LAPACK wrappers; see accelerate_solver.h. #include "accelerate_solver.h" #ifdef USE_ACCEL_BACKEND diff --git a/alm/accelerate_solver.h b/alm/accelerate_solver.h index aec2eb8d..dbddddc5 100644 --- a/alm/accelerate_solver.h +++ b/alm/accelerate_solver.h @@ -1,12 +1,6 @@ -// accelerate_solver.h -// -// Apple Accelerate sparse KKT solver, deliberately isolated in its own translation unit. -// -// Eigen's pulls in , whose Fortran BLAS/LAPACK -// prototypes (dgemm_, dgesdd_, dgeqp3_, ...) collide with the hand-rolled prototypes in -// include/blas_wrapper.h and include/lapack_wrapper.h. Keeping the Accelerate include out of every -// TU that also includes those wrappers (least_squares.cpp in particular) avoids a hard redeclaration -// conflict. Only this header/source pair includes the Accelerate Eigen module. +// Apple Accelerate sparse KKT solver. Keep Eigen/AccelerateSupport in its own +// translation unit to avoid conflicting BLAS/LAPACK declarations from +// blas_wrapper.h and lapack_wrapper.h. #pragma once #ifdef USE_ACCEL_BACKEND diff --git a/alm/alm.cpp b/alm/alm.cpp index 00621933..50c086bb 100644 --- a/alm/alm.cpp +++ b/alm/alm.cpp @@ -474,10 +474,8 @@ auto ALM::get_number_of_fc_elements(const int fc_order) const -> size_t auto ALM::get_number_of_irred_fc_elements(const int fc_order) -> size_t // harmonic=1, ... { - // Returns the number of irreducible force constants for the given order. - // The irreducible force constant means a set of independent force constants - // reduced by using all available symmetry operations and - // constraints for translational invariance. Rotational invariance is not considered. + // Count independent force constants after crystal-symmetry and translational + // constraints, excluding rotational invariance. const auto order = fc_order - 1; if (!initialized_constraint_class) { @@ -867,9 +865,7 @@ auto ALM::init_fc_table() -> void cluster->init(system, symmetry, get_optimizer_control().periodic_image_conv, verbosity, timer); fcs->init(cluster, symmetry, system->get_supercell(), verbosity, timer); - // Switch off the initialized_constraint_class flag - // because the force constants are updated - // but corresponding constranits are not. + // Invalidate constraints after updating the force constants. initialized_constraint_class = false; } diff --git a/alm/alm.h b/alm/alm.h index 0f248e86..1ca08fdc 100644 --- a/alm/alm.h +++ b/alm/alm.h @@ -73,10 +73,8 @@ class ALM auto set_periodicity(const int is_periodic[3]) const -> void; - // Declare the units of the data passed to the unit-sensitive setters - // (set_cell, set_u_train, set_f_train, set_validation_data, define). - // Must be called before any of them; the setters convert the input to the - // internal canonical units (bohr, Ry/bohr). LENGTH_UNIT / FORCE_UNIT. + // Set input units before set_cell, set_u_train, set_f_train, + // set_validation_data, or define. Setters convert to bohr and Ry/bohr. auto set_input_units(const std::string &length_unit, const std::string &force_unit) -> void; // Canonical names of the declared input units: {length, force}. diff --git a/alm/cluster.cpp b/alm/cluster.cpp index 665b524d..0a96489b 100644 --- a/alm/cluster.cpp +++ b/alm/cluster.cpp @@ -1063,10 +1063,8 @@ void Cluster::set_interaction_cluster(const int order, const size_t natmin, cons const auto rc_tmp = cutoff_radii[order][ikd][jkd]; cell_vector.clear(); - // Loop over the cell images of atom 'jat' and add to the list - // as a candidate for the cluster. - // The periodic images whose distance is larger than the minimum value - // of the distance(iat, jat) can be added to the cell_vector list. + // Collect candidate periodic images of jat, including images beyond + // the minimum iat-jat distance. for (auto k = 0; k < distance_table[iat][jat].distances.size(); ++k) { if (rc_tmp < 0.0 || distance_table[iat][jat].distances[k] <= rc_tmp) { cell_vector.emplace_back(distance_table[iat][jat].cells[k]); @@ -1123,11 +1121,8 @@ void Cluster::set_interaction_cluster(const int order, const size_t natmin, cons // that satisfies the condition of the cluster. if (periodic_image_conv == 0) { - // assign IFCs to periodic images in which the center atom and each of the other atoms - // are nearest. - // The distance between non-center atoms are not considered. - // The IFCs in this convention automatically satisfies the ASR without additional constraint, - // but does not satisfy the permutation symmetry. + // Assign IFCs to images nearest to the center atom, ignoring distances + // between other atoms. This satisfies ASR but not permutation symmetry. pairs_icell.clear(); for (const int jat: intpair_uniq) { @@ -1164,11 +1159,8 @@ void Cluster::set_interaction_cluster(const int order, const size_t natmin, cons } else /* if(mirror_image_conv == 1)*/ { - // assign IFCs to periodic images in which the sum of the distances between the atom pairs - // is the smallest. - // The IFCs made in this convention satisfies the permutation symmetry. - // Additional constraints are imposed in constraint.cpp to make the IFCs satisfy ASR - // after assigning IFCs to the periodic images. + // Assign IFCs to images minimizing the sum of pair distances. This preserves + // permutation symmetry; constraint.cpp imposes the additional ASR constraints. std::sort(distance_list.begin(), distance_list.end(), MinDistList::compare_sum_distance); comb_cell_min.clear(); diff --git a/alm/cluster.h b/alm/cluster.h index a11f74ab..1852383d 100644 --- a/alm/cluster.h +++ b/alm/cluster.h @@ -60,10 +60,7 @@ class IntList auto operator==(const IntList &a) const -> bool { - // Use vector's built-in equality operator for consistency - // This is required for std::unordered_set in fcs.cpp - // and ensures consistent behavior across platforms - // std::cout << "== operator called\n"; + // Vector equality supports std::unordered_set in fcs.cpp. return iarray == a.iarray; } }; diff --git a/alm/constraint.cpp b/alm/constraint.cpp index 31e7fb43..0fd20e45 100644 --- a/alm/constraint.cpp +++ b/alm/constraint.cpp @@ -344,16 +344,9 @@ auto Constraint::update_constraint_matrix(const std::unique_ptr &system, { const auto maxorder = cluster->get_maxorder(); - // For the non-algebraic path (ICONST = 1/2/3) the merged constraint matrix is rank-reduced - // and handed to the equality-constrained least-squares solver (GQR / Pardiso-KKT). Use the - // rank-revealing QR (qrd) there regardless of the configured backend: the rref / - // coord_factorization (Gauss-Jordan) backends use a fixed absolute pivot tolerance that - // mis-detects the numerical rank of large constraint systems -- they accept round-off - // residuals as independent pivots and, after normalising by the tiny pivot, inject - // amplified-noise constraint rows that corrupt the fit (broken ASR / rotational invariance, - // imaginary phonons). The configured backend is honoured for the algebraic path - // (ICONST >= 10), where it builds the per-order elimination map consumed by the - // elastic-net / adaptive-lasso solvers. + // Use rank-revealing QR for ICONST = 1/2/3: fixed absolute pivot tolerances + // can mistake round-off for independent constraints and corrupt the fit. + // For ICONST >= 10, use the configured backend for the elimination map. const auto algo = constraint_algebraic ? algo_in : ReductionAlgo::qrd; // const_symmetry is updated. @@ -433,16 +426,9 @@ auto Constraint::update_constraint_matrix(const std::unique_ptr &system, } } - // The merged (full) constraint matrix const_mat_sparse and its rank reduction are only - // consumed by the non-algebraic solvers (build_constraint_matrix_dense and the dense - // least_squares_with_constraints_gqr / sparse solveGQRSparse paths). In the algebraic - // path (ICONST >= 10, used by elastic-net / adaptive-lasso) the solve instead uses the - // per-order elimination map (const_fix / const_relate / index_bimap) built below, and - // get_exist_constraint() inspects const_self / const_fix / const_relate directly without - // reading number_of_constraints. Building/reducing the merged matrix here is therefore - // wasted work in the algebraic path -- and for large maxorder it dominates the reduction - // time, because get_independent_rows_lapack_sparse densifies to a P x N buffer and runs a - // dense LAPACK QR. Skip it entirely; number_of_constraints is set from the mapping below. + // Only non-algebraic solvers need the merged constraint matrix. + // For ICONST >= 10, use the per-order elimination map below and avoid + // the costly merged-matrix rank reduction. if (!constraint_algebraic) { size_t nparams = 0; for (auto order = 0; order < maxorder; ++order) { @@ -467,9 +453,7 @@ auto Constraint::update_constraint_matrix(const std::unique_ptr &system, } if (constraint_algebraic) { - // The mapping requires const_self in reduced row echelon form. rref and - // coord_factorization already produced it in the per-order loop above; only the qrd / - // none paths still need an explicit echelon reduction here. + // The mapping needs reduced row echelon form; only qrd / none still need reduction. if (algo != ReductionAlgo::rref && algo != ReductionAlgo::coord_factorization) { for (auto order = 0; order < maxorder; ++order) { const auto nparam = fcs->get_nequiv()[order].size(); @@ -484,9 +468,7 @@ auto Constraint::update_constraint_matrix(const std::unique_ptr &system, const_relate.data(), index_bimap); - // number_of_constraints is not consumed by the algebraic solve, but keep it - // meaningful for reporting and get_exist_constraint(): each fixed or related - // parameter corresponds to one independent constraint that eliminates a parameter. + // Count one constraint per fixed or related parameter for reporting. number_of_constraints = 0; for (auto order = 0; order < maxorder; ++order) { number_of_constraints += const_fix[order].size() + const_relate[order].size(); @@ -1880,10 +1862,8 @@ auto Constraint::generate_rotational_constraint(const std::unique_ptr &s rank); } } else if (algo_in == ReductionAlgo::coord_factorization) { - // Mirror the rref branch with the stable partial-pivot kernel. const_rotation_cross - // is inter-order and is NOT merged into const_self, so (unlike the other subsets) it - // must be reduced here; otherwise it would reach build_constraint_matrix_sparse raw - // in the non-algebraic path (ICONST = 2/3). + // Reduce inter-order const_rotation_cross here: it is not merged into + // const_self before the non-algebraic solve (ICONST = 2/3). rref_sparse_pivot(nparams[order], const_rotation_self[order], eps6); if (order > 0) { rref_sparse_pivot(nparams[order - 1] + nparams[order], const_rotation_cross[order], eps6); diff --git a/alm/fcs.cpp b/alm/fcs.cpp index 2c3d983b..cebaf2b5 100644 --- a/alm/fcs.cpp +++ b/alm/fcs.cpp @@ -313,10 +313,8 @@ auto Fcs::get_constraint_symmetry(const size_t nat, const std::unique_ptr void { - // Create constraint matrices arising from the crystal symmetry. - // Necessary for hexagonal systems. - // The QR path intentionally preserves the current auto rank tolerance. `tolerance` is retained - // for API compatibility with callers and the non-integer symmetry helper. + // Build crystal-symmetry constraints (needed for hexagonal systems). + // QR uses automatic rank tolerance; retain tolerance for API compatibility. (void)tolerance; int i; @@ -487,11 +485,7 @@ auto Fcs::get_constraint_symmetry_in_integer(const size_t nat, const std::unique const double tolerance, ConstraintSparseForm &const_out, const ReductionAlgo algo_in) -> void { - // Create constraint matrices arising from the crystal symmetry. - // Necessary for hexagonal systems. - // This does exactly the same thing as get_constraint_symmetry but assumes - // all elements of the constraint matrix are integer. - // This version is expected to be more stable (and fast?). + // Integer-matrix version of get_constraint_symmetry for improved stability. int i; unsigned int isym; diff --git a/alm/input_parser.cpp b/alm/input_parser.cpp index 2de49690..a7b2e443 100644 --- a/alm/input_parser.cpp +++ b/alm/input_parser.cpp @@ -178,10 +178,8 @@ auto InputParser::parse_displacement_and_force_files(std::vector &energies, const DispForceFile &datfile_in) const -> void { - // Parse total-supercell reference energies from the "E_pot (unit):" token in each snapshot's - // comment header of the DFSET, convert to Rydberg, and apply the SAME NSTART/NEND/SKIP filtering - // as parse_displacement_and_force_files so that `energies` aligns with u_train/f_train. - // Exactly one parseable E_pot per ORIGINAL snapshot is required (validated before filtering). + // Read one E_pot per DFSET snapshot and convert to Ry. Validate before + // applying the same NSTART/NEND/SKIP filter as displacements and forces. const auto nat = nat_in * static_cast(transmat_to_super.determinant()); const auto ntoken_per_snapshot = static_cast(6 * nat); @@ -260,14 +258,7 @@ auto InputParser::parse_energies(std::vector &energies, const DispForceF auto InputParser::parse_input(ALM *alm) -> void { - // The order of calling methods in this method is important. - // Since following methods rely on variables those already - // parsed. - // Those below are set as the private class variables. See input_parser.h. - // std::string mode; - // int maxorder; - // int nat_base; - // int nkd; + // Parse in this order: later methods depend on values set by earlier ones. // Parse &general field if (!locate_tag("&general")) { @@ -293,11 +284,8 @@ auto InputParser::parse_input(ALM *alm) -> void nat_in = atomic_types_input.size(); nkd_in = kdname_vec.size(); } else { - // If STRUCTURE_FILE is given, use the structure parameters defined in this file. - // In this case, the &cell and &position entries are ignored. - // A POSCAR is always in Angstrom (VASP convention), independent of LENGTH_UNIT. - // ALM::set_cell converts its input from LENGTH_UNIT to bohr, so re-express the - // POSCAR lattice in LENGTH_UNIT here to end up with a single net conversion. + // STRUCTURE_FILE overrides &cell and &position. Convert the POSCAR lattice + // from Angstrom to LENGTH_UNIT before set_cell converts it to bohr. if (length_unit_input == "bohr") { lavec_poscar /= Bohr_in_Angstrom; } @@ -424,10 +412,7 @@ auto InputParser::parse_general_vars(ALM *alm) -> void } if (mode == "opt") mode = "optimize"; - // We first check if STRUCTURE_FILE field is empty or not. - // If not, the structure data is read from the given file (in a POSCAR format) - // and copy the lattice vectors, element types, and coordinates to - // the corresponding private variables of this class. + // Read lattice vectors, species, and coordinates from STRUCTURE_FILE if set. if (!general_var_dict["STRUCTURE_FILE"].empty()) { structure_file = general_var_dict["STRUCTURE_FILE"]; struct stat buffer; @@ -576,11 +561,8 @@ auto InputParser::parse_general_vars(ALM *alm) -> void format_pattern = "yaml"; } - // Units of the input data (&cell lattice, DFSET displacements and forces, - // &cutoff radii) and of the alamode_h5 output. Validated here; the actual - // conversion to the internal Rydberg atomic units happens in the ALM core - // setters. NOTE: STRUCTURE_FILE (POSCAR) is always in Angstrom and E_pot - // energies keep their own per-snapshot (eV/Ry/Ha) header mechanism. + // Validate input and HDF5 output units; core setters perform conversion. + // POSCAR always uses Angstrom; E_pot uses its per-snapshot unit header. std::string length_unit{"bohr"}, force_unit{"Ry/bohr"}, fcs_unit_output{"Ry/bohr"}; if (!general_var_dict["LENGTH_UNIT"].empty()) { length_unit = general_var_dict["LENGTH_UNIT"]; @@ -920,9 +902,7 @@ auto InputParser::parse_structure_poscar(const std::string &fname_poscar, Eigen: Eigen::MatrixXd &coordinates_out, std::vector &kdname_vec_out, std::vector &atomic_types_out) -> void { - // Parse structure data from a file in the POSCAR format and - // copy the read data to the corresponding private variables of - // this class. + // Read POSCAR structure data into the parser fields. std::ifstream ifs; std::string dummy; @@ -1422,9 +1402,8 @@ auto InputParser::parse_optimize_vars(ALM *alm) -> void exit("parse_optimize_vars", "NDATA, NSTART, NEND and SKIP tags are inconsistent."); } - // Energy-difference loss term: read reference energies from the DFSET headers, but only for a - // production (non-CV) fit. During CV the energy pass is skipped entirely so legacy headerless - // DFSETs are unaffected (alpha selection stays force-only). + // Read DFSET energies for production fits only; force-only CV also accepts + // headerless DFSETs. if (optcontrol.efit_weight > 0.0) { if (optcontrol.cross_validation == 0 || optcontrol.efit_cv) { // Production fit (CV=0), or CV with EFIT_CV=1: the training energies are needed. @@ -1515,9 +1494,7 @@ auto InputParser::parse_optimize_vars(ALM *alm) -> void int ialgo_reduce; if (optimize_var_dict["ALGO_REDUCTION"].empty()) { - // Default: 3 = coord_factorization (partial-pivot RREF). Stable replacement for the - // legacy 1 = rref backend; produces identical maps but never divides by a small pivot. - // Set ALGO_REDUCTION = 1 to reproduce the legacy rref behavior. + // Default to partial-pivot RREF (3); ALGO_REDUCTION = 1 selects legacy RREF. ialgo_reduce = 3; } else { assign_val(ialgo_reduce, "ALGO_REDUCTION", optimize_var_dict); diff --git a/alm/input_setter.cpp b/alm/input_setter.cpp index abc46526..b03a22dc 100644 --- a/alm/input_setter.cpp +++ b/alm/input_setter.cpp @@ -192,9 +192,8 @@ auto InputSetter::set_transformation_matrices(const Eigen::Matrix3d &transmat_su const Eigen::Matrix3d &transmat_prim_in, const int autoset_primcell_in, const bool transpose) -> void { - // if the input transformation matrices are defined by (a_s, b_s, c_s)^T = M (a_p, b_p, c_p)^T, - // which is more understandable for human, we need to transpose the matrices to make it consistent - // with the definition of the lattice vectors used in the code. + // Transpose the row-vector input convention (a_s, b_s, c_s)^T = M (a_p, b_p, c_p)^T + // to match the column-vector convention used internally. if (transpose) { transmat_super = transmat_super_in.transpose(); transmat_prim = transmat_prim_in.transpose(); diff --git a/alm/least_squares.cpp b/alm/least_squares.cpp index d245cef1..efec06c0 100644 --- a/alm/least_squares.cpp +++ b/alm/least_squares.cpp @@ -17,23 +17,15 @@ #ifdef USE_MKL_BACKEND #include #elif defined(USE_ACCEL_BACKEND) -// The Accelerate KKT solve lives in accelerate_solver.cpp; its include -// pulls , whose Fortran BLAS/LAPACK prototypes would clash with the -// hand-rolled ones in blas_wrapper.h / lapack_wrapper.h that this TU includes below. Keeping the -// Accelerate header out of this TU is exactly why that solve is isolated. +// Keep Accelerate headers in accelerate_solver.cpp to avoid conflicts +// with the BLAS/LAPACK wrappers below. #include "accelerate_solver.h" #endif #ifdef USE_SUITESPARSE_BACKEND -// SuiteSparse sparse solvers, orthogonal to the MKL/Accelerate KKT backend above and reachable as -// SPARSESOLVER options: -// - CHOLMOD supernodal Cholesky on the PSD normal matrix A^T A, via Eigen's CholmodSupport wrapper. -// - SuiteSparseQR (rank-revealing multifrontal QR) on the rectangular A, via SuiteSparse's own C -// interface. We deliberately avoid Eigen 3.4's Eigen::SPQR wrapper: with SuiteSparse >= 6 it -// instantiates SuiteSparseQR(...) passing an Eigen::Index ('long') where the templated -// index type resolves to int64_t ('long long' on macOS arm64), and the two deduce to conflicting -// types -- a hard compile error. SuiteSparseQR_C_backslash_default() sidesteps it entirely. -// Both headers live under /include/suitesparse, on the include path via the linked targets. +// SuiteSparse backends: CHOLMOD for A^T A, SuiteSparseQR for rectangular A. +// Use the QR C interface to avoid Eigen 3.4 / SuiteSparse >= 6 index-type +// conflicts on macOS arm64. #include #include #include @@ -43,10 +35,8 @@ #include "lapack_wrapper.h" #ifdef USE_SUITESPARSE_BACKEND -// Solve min_x || A x - b ||_2 (least squares for rectangular A; exact solve for square full-rank A) -// with SuiteSparseQR's C interface. The backslash convenience routine is 64-bit-index only, so the -// matrix is shipped as a CHOLMOD_LONG cholmod_sparse: Eigen stores 32-bit CSC indices, widened here. -// Returns 0 on success, 1 on failure. See the include block above for why Eigen::SPQR is bypassed. +// Solve min_x ||A x - b||_2 with SuiteSparseQR. Widen Eigen CSC indices to +// CHOLMOD_LONG for the 64-bit-only backslash interface. Returns 0 on success, 1 on failure. static auto solve_least_squares_spqr(const Eigen::SparseMatrix &A, const Eigen::VectorXd &b, Eigen::VectorXd &x_out) -> int { @@ -108,15 +98,10 @@ static auto solve_least_squares_spqr(const Eigen::SparseMatrix &A, const return status; } -// Rank-revealing reduction of a HOMOGENEOUS constraint matrix C (rows define C x = 0) via SuiteSparseQR -// -- a sparse, multithreaded replacement for the densifying LAPACK dgeqp3 path. We factorize the TALL -// matrix C^T (N x P) with rank detection and column pivoting: the rank-revealing pivot deflates the -// dependent columns of C^T (= dependent rows of C) to the end of the permutation E, so E[0..rank-1] are -// the independent rows of C. C_red is then those ORIGINAL (sparse) rows of C -- NOT the R factor, which -// is densely filled for these constraints. R (only rank x P here) and Q are discarded. -// Only the homogeneous case is handled (every invariance subset, and the merged matrix when no -// FC2FIX/FC3FIX value is imposed); inhomogeneous d != 0 and any SuiteSparseQR failure fall back to -// dgeqp3. Returns 0 on success, 1 on failure. +// Select independent rows of homogeneous C x = 0 using pivoted QR of C^T. +// Keep the original sparse rows indexed by E[0..rank-1], not the dense R factor. +// Returns 0 on success, 1 on failure; callers use dgeqp3 for inhomogeneous +// constraints or QR failure. static auto get_independent_rows_spqr(const Eigen::SparseMatrix &C, const int verbosity, const double tolerance, Eigen::SparseMatrix &C_red, int &r) -> int { @@ -151,11 +136,9 @@ static auto get_independent_rows_spqr(const Eigen::SparseMatrix &C, cons A.sorted = 1; A.packed = 1; - // Map the caller's auto sentinel (tol < 0, i.e. rank_tolerance_auto = -1) to SuiteSparseQR's default - // tolerance SPQR_DEFAULT_TOL (= -2). Passing -1 literally would be SPQR_NO_TOL (rank detection off). - // An explicit positive tolerance is passed through. NOTE: SuiteSparseQR thresholds column 2-norms - // whereas dgeqp3 thresholds R-diagonals, so the numerical rank may differ by a few rows on a - // borderline (badly-scaled) constraint set. + // Map the auto sentinel (-1) to SPQR_DEFAULT_TOL (-2); SPQR uses -1 to + // disable rank detection. SPQR thresholds column norms, while dgeqp3 + // thresholds R diagonals, so borderline numerical ranks may differ. const double spqr_tol = (tolerance < 0.0) ? SPQR_DEFAULT_TOL : tolerance; cholmod_sparse *R = nullptr; @@ -388,10 +371,8 @@ auto get_independent_rows_lapack_sparse(const Eigen::SparseMatrix &C_spa LOG_IF(verbosity, 1, "P = ", P, ", N = ", N, ", nnz = ", C_sparse.nonZeros(), ".\n"); #ifdef USE_SUITESPARSE_BACKEND - // SuiteSparseQR rank-revealing reduction for HOMOGENEOUS constraints (every invariance subset, and - // the merged matrix when no FC2FIX/FC3FIX value is imposed): sparse + multithreaded, avoiding the - // dense P x N densification and dgeqp3 below. Inhomogeneous d (d != 0, from FC2FIX/FC3FIX) and any - // SuiteSparseQR failure fall through to the dgeqp3 path. + // Use sparse QR for homogeneous constraints; fall back to dgeqp3 for + // inhomogeneous constraints (FC2FIX/FC3FIX) or QR failure. if (dvec.isZero(0)) { if (get_independent_rows_spqr(C_sparse, verbosity, tolerance, C_red, r) == 0) { d_red = Eigen::VectorXd::Zero(r); @@ -742,15 +723,11 @@ auto least_squares_with_constraints_svd(const size_t N, const size_t M, const si return 0; } -// 64-bit-index sparse matrix for the normal matrix A^T A. The default SparseMatrix uses a 32-bit -// StorageIndex, whose cumulative outer index overflows once nnz(A^T A) exceeds ~2.1e9 -- which happens -// for large full-range / no-cutoff fits (e.g. ICONST >= 10 with `*-* None`). The overflow corrupts the -// allocation size and aborts with std::bad_alloc even when memory is free. 64-bit indices avoid it. +// Use 64-bit indices for A^T A to avoid overflow above 2^31-1 nonzeros. using SpMat64 = Eigen::SparseMatrix; -// Form A^T A with 64-bit indices. If A^T A is genuinely too large to allocate (it is N x N and can be -// near-dense for no-cutoff problems), abort with guidance to use SuiteSparseQR, which factorizes the -// rectangular A directly and never forms A^T A. (A^T A is computed here, not the dense factor.) +// Form A^T A with 64-bit indices. On allocation failure, recommend +// SuiteSparseQR, which avoids forming the normal matrix. [[nodiscard]] static auto build_normal_matrix_int64(const Eigen::SparseMatrix &A) -> SpMat64 { try { @@ -857,12 +834,8 @@ auto least_squares_eigen_sparse_solver(const Eigen::SparseMatrix &sp_mat } else if (solver_type_lower == "cholmod") { #ifdef USE_SUITESPARSE_BACKEND - // CHOLMOD supernodal Cholesky on the normal matrix A^T A (symmetric positive (semi)definite). - // Fast for well-conditioned problems; for a rank-deficient A^T A the factorization fails and - // info() reports it, in which case SuiteSparseQR is the recommended fallback. Check the - // factorization BEFORE solving -- solving on a failed factor is undefined. A^T A uses 64-bit - // indices (and CHOLMOD's long interface) so large no-cutoff fits do not overflow the 32-bit - // sparse index; for genuinely huge A^T A, build_normal_matrix_int64 aborts pointing to SuiteSparseQR. + // Factor A^T A with CHOLMOD using 64-bit indices. Check factorization before + // solving; rank-deficient matrices may require SuiteSparseQR. SpMat64 AtA = build_normal_matrix_int64(sp_mat); Eigen::VectorXd AtB = sp_mat.transpose() * sp_bvec; Eigen::CholmodSupernodalLLT chol(AtA); @@ -894,16 +867,10 @@ auto least_squares_eigen_sparse_solver(const Eigen::SparseMatrix &sp_mat } -// Block-diagonal SPD preconditioner for the symmetric-indefinite KKT system -// K = [ H C^T ; C 0 ], H = A^T A (N x N, rank-deficient by the gauge modes), -// for use with MINRES. P = diag(G, S) with -// G = diag(H) + sigma (a positive diagonal -- a Jacobi approximation of H) -// S = C G^{-1} C^T (dense, SPD; the preconditioner's exact Schur complement for that G). -// By Murphy-Golub-Wathen this clusters the spectrum of P^{-1}K, so preconditioned MINRES converges in -// far fewer iterations than the (stalling) unpreconditioned run. A *diagonal* G is used deliberately: -// G^{-1} is then trivial, so the setup needs only sparse products to form S -- no N x P solves, which -// for P ~ 3000 dense rotational constraints would dominate the runtime. sigma keeps G strictly -// positive; it affects only the preconditioner's quality (iteration count), never the final solution. +// SPD preconditioner for MINRES on K = [H C^T; C 0], H = A^T A: +// P = diag(G, S), G = diag(H) + sigma, S = C G^{-1} C^T. +// Diagonal G makes setup inexpensive; sigma keeps it positive and affects +// convergence only, not the solution. class KKTBlockDiagPreconditioner { public: @@ -943,13 +910,9 @@ class KKTBlockDiagPreconditioner const SpMat Cs = C * m_Dinv.asDiagonal(); // P x N (scaled columns) Eigen::MatrixXd S = Eigen::MatrixXd(Cs * C.transpose()); // P x P dense - // S can be extremely ill-conditioned -- even when C is full row rank -- because diag(A^T A) - // has near-zero entries (weakly-sampled parameters) that make Dinv, and hence S, span a huge - // dynamic range. (This is NOT evidence that C is rank-deficient: C is full row rank after the - // merge + rank-reduction in Constraint::update_constraint_matrix.) MINRES requires a positive- - // definite preconditioner, so add a relative ridge to S and grow it until the dense factor is - // SPD. The ridge perturbs only the preconditioner -- never the system matrix K, so it does not - // bias the solution; it only affects how fast MINRES converges. + // Near-zero entries of diag(A^T A) can make S ill-conditioned even for + // full-rank C. Increase a relative ridge until S is SPD, as MINRES requires. + // This changes only the preconditioner, not the solution. const double s_scale = S.diagonal().cwiseAbs().mean() + sigma; double ridge = 1.0e-8 * s_scale; m_ready = false; @@ -1066,13 +1029,9 @@ void solveGQRSparse(const Eigen::SparseMatrix &A, const Eigen::VectorXd Eigen::VectorXd sol; bool solved = false; - // Iterative KKT solve (selected via SPARSESOLVER). The KKT matrix K is sparse (it never forms the - // dense factor that a direct method would), so this is the memory-feasible path when the constraint - // matrix C has dense rows -- e.g. the rotational-invariance constraints of ICONST = 2, 3 -- for - // which the direct factorizations below exhaust memory. K is symmetric *indefinite*, so MINRES is - // the appropriate Krylov method, accelerated by the block-diagonal SPD preconditioner - // diag(diag(A^T A) + sigma, C G^{-1} C^T) (KKTBlockDiagPreconditioner above). An unpreconditioned / - // identity-preconditioned MINRES stalls completely on this ill-conditioned saddle-point system. + // Use preconditioned MINRES for the symmetric-indefinite KKT system when + // direct factors are too large, especially with dense rotational constraints. + // KKTBlockDiagPreconditioner improves convergence on this ill-conditioned system. const auto solver_lower = boost::algorithm::to_lower_copy(solver_type); if (solver_lower == "minres") { LOG_IF(verbosity, 1, "Use MINRES (iterative) to solve the KKT problem.\n"); @@ -1093,13 +1052,8 @@ void solveGQRSparse(const Eigen::SparseMatrix &A, const Eigen::VectorXd } minres.compute(K); // sets the matrix; the preconditioner was already configured above - // MINRES monitors the *preconditioned* residual, which with this block-diagonal preconditioner - // can be far smaller than the TRUE residual -- so a single solve at the user's CONV_TOL may stop - // before the constraints (sum rules) are actually satisfied. Re-solve with a progressively - // tighter stopping tolerance until the TRUE relative KKT residual meets the acceptance bar; the - // preconditioner is cheap and each solve takes only a few iterations. Acceptance is judged on the - // true residual (and the explicit constraint residual ||C x - d||), never on minres.info() alone, - // so an under-converged x that would violate the sum rules is never returned. + // Tighten MINRES tolerance until the true KKT and constraint residuals pass. + // Its preconditioned residual alone can hide unsatisfied sum rules. const double rhs_norm = (rhs.norm() > 0.0) ? rhs.norm() : 1.0; const double accept_tol = std::max(tolerance_iteration, 1.0e-10); double kkt_rel_res = std::numeric_limits::infinity(); @@ -1181,12 +1135,8 @@ void solveGQRSparse(const Eigen::SparseMatrix &A, const Eigen::VectorXd #endif #ifdef USE_SUITESPARSE_BACKEND - // Preferred direct solver when ALM is built with SuiteSparse: SuiteSparseQR is multithreaded and - // rank-revealing, so on a well-conditioned KKT it is far faster than Eigen's serial SparseLU / - // SparseQR. It is tried before SparseLU (after any symmetric-indefinite LDLT backend). (Earlier - // this was a post-SparseLU fallback to avoid QR fill-in blowing up memory when the constraint - // matrix was rank-deficient with dense rows; that pathology came from an unreduced constraint - // matrix and is prevented upstream by the rank-revealing reduction in Constraint.) + // Try multithreaded, rank-revealing SuiteSparseQR before Eigen SparseLU. + // Constraint rows have already been rank-reduced upstream. if (!solved) { LOG_IF(verbosity, 1, "Use SuiteSparseQR to solve the KKT problem.\n"); if (solve_least_squares_spqr(K, rhs, sol) == 0) { diff --git a/alm/least_squares.h b/alm/least_squares.h index 61a3d679..e6e74605 100644 --- a/alm/least_squares.h +++ b/alm/least_squares.h @@ -10,77 +10,27 @@ #include "constraint.h" /** - * @brief Solves an unconstrained least squares problem. - * - * This function minimizes ||A x - b||_2, where: - * - A is the data matrix (M×N) - * - b is the observation vector (length M) - * - x is the solution vector (length N) - * - * The method computes the least squares solution using standard numerical techniques. - * - * @param N Number of columns of A (number of unknowns). - * @param M Number of rows of A (number of data points). - * @param amat Pointer to A in column-major order (size M×N). - * @param bvec Pointer to b (size M). - * @param param_out Output array for solution x (size N). - * @param verbosity If > 0, print diagnostic messages. - * @return 0 for success, >0 for failure. + * Minimize ||A x - b||_2 using SVD. amat is column-major M x N; + * bvec has length M and param_out length N. verbosity > 0 prints diagnostics. + * Returns 0 on success, a nonzero solver status on failure. */ auto least_squares_svd(const size_t N, const size_t M, double *amat, const double *bvec, double *param_out, const int verbosity) -> int; /** - * @brief Solves a constrained least squares problem using GQR-based decomposition. - * - * This function minimizes ||A x - b||_2 subject to C x = d, where: - * - A is the data matrix (M×N) - * - b is the observation vector (length M) - * - C is the constraint matrix (P×N) - * - d is the constraint vector (length P) - * - * The method uses QR decomposition with column pivoting to identify linearly independent - * rows of C, extracts a reduced constraint matrix C_red and corresponding d_red, and - * then calls LAPACK's dgglse to compute the constrained least squares solution. - * - * @param N Number of columns of A (number of unknowns). - * @param M Number of rows of A (number of data points). - * @param P Number of rows of C (number of constraints). - * @param amat Pointer to A in column-major order (size M×N). - * @param bvec Pointer to b (size M). - * @param param_out Output array for solution x (size N). - * @param cmat Array of P pointers, each pointing to an array of length N (C in row-major). - * @param dvec Pointer to d (size P). - * @param verbosity If > 0, print diagnostic messages. - * @return INFO from dgglse (0 for success, >0 for failure). + * Minimize ||A x - b||_2 subject to C x = d using rank reduction and dgglse. + * amat is column-major M x N; cmat is a P x N row-pointer array. + * bvec, dvec, and param_out have lengths M, P, and N. + * verbosity > 0 prints diagnostics. Returns the LAPACK status (0 on success). */ auto least_squares_with_constraints_gqr(const size_t N, const size_t M, const size_t P, double *amat, const double *bvec, double *param_out, const double *const *cmat, const double *dvec, const int verbosity) -> int; /** - * @brief Solves a constrained least squares problem using SVD-based decomposition. - * - * This function minimizes ||A x - b||_2 subject to C x = d, where: - * - A is the data matrix (M×N) - * - b is the observation vector (length M) - * - C is the constraint matrix (P×N) - * - d is the constraint vector (length P) - * - * The method computes the pseudoinverse of C to find a particular solution x₀ = C⁺d, - * projects the residual b' = b - A x₀ onto the null space of C, and solves the - * least squares problem in the null space using SVD. - * - * @param N Number of columns of A (number of unknowns). - * @param M Number of rows of A (number of data points). - * @param P Number of rows of C (number of constraints). - * @param amat Pointer to A in column-major order (size M×N). Can be overwritten. - * @param bvec Pointer to b (size M). Can be overwritten. - * @param param_out Output array for solution x (size N). - * @param cmat Array of P pointers, each pointing to an array of length N (C in row-major). - * @param dvec_orig Pointer to d (size P). Will be copied locally. - * @param verbosity If > 0, print diagnostic messages. - * @return 0 for success, >0 for failure. + * Minimize ||A x - b||_2 subject to C x = d using SVD and the + * null-space projector I - C^+ C, starting from x0 = C^+ d. + * verbosity > 0 prints diagnostics. Returns 0 on success, nonzero on failure. */ auto least_squares_with_constraints_svd(const size_t N, const size_t M, const size_t P, double *amat, // A: (M×N) column-major, can be overwritten @@ -92,22 +42,10 @@ auto least_squares_with_constraints_svd(const size_t N, const size_t M, const si /** - * @brief Solves a sparse least squares problem using Eigen's sparse solvers. - * - * This function minimizes ||A x - b||_2, where: - * - A is a sparse matrix represented using Eigen's SparseMatrix. - * - b is the observation vector (Eigen::VectorXd). - * - x is the solution vector (Eigen::VectorXd). - * - * The method uses a specified sparse solver to compute the least squares solution. - * - * @param sp_mat Sparse matrix A (M×N) in Eigen::SparseMatrix format. - * @param sp_bvec Observation vector b (length M) in Eigen::VectorXd format. - * @param x_out Output vector for the solution x (length N). - * @param solver_type String specifying the solver type (e.g., "CG", "BiCGSTAB"). - * @param tolerance_iteration Convergence tolerance for the iterative solver. - * @param maxnum_iteration Maximum number of iterations for the solver. - * @return 0 for success, >0 for failure. + * Minimize ||A x - b||_2 with the selected sparse solver. + * sp_mat is M x N, sp_bvec has length M, and x_out has length N. + * tolerance_iteration and maxnum_iteration control iterative convergence. + * Returns 0 on success, nonzero on failure. */ auto least_squares_eigen_sparse_solver(const Eigen::SparseMatrix &sp_mat, const Eigen::VectorXd &sp_bvec, Eigen::VectorXd &x_out, const std::string &solver_type, @@ -115,49 +53,18 @@ auto least_squares_eigen_sparse_solver(const Eigen::SparseMatrix &sp_mat /** - * @brief Build reduced constraint matrix C_red (r×N) and vector d_red (length r). - * - * This function takes the original constraint matrix C (size P×N) and constraint vector d - * (length P), determines the numerical row rank r of C via QR with column pivoting on C^T, - * and extracts the first r independent rows. The outputs C_red and d_red are stored in - * column-major order. - * - * @param[in] N Number of variables (number of columns of C) - * @param[in] P Number of constraints (number of rows of C) - * @param[in] cmat Pointer to constraint matrix C: cmat[i][j] is row i, column j (0 ≤ i < P, 0 ≤ j < N) - * @param[in] dvec Original constraint vector d (length P) - * @param[in] verbosity Verbosity level (0: silent, >0: print info) - * @param[out] C_red Output reduced constraint matrix (r×N) in column-major layout - * @param[out] d_red Output reduced constraint vector (length r) - * @param[out] r Computed numerical row rank (0 ≤ r ≤ min(P, N)) - * - * @return 0 on success, nonzero LAPACK info code on failure + * Select independent rows of cmat[P][N] using pivoted QR of C^T. + * Return column-major C_red[r][N], matching d_red[r], and numerical rank r. + * verbosity > 0 prints diagnostics. Returns 0 on success or LAPACK INFO. */ auto get_independent_rows(const size_t N, const size_t P, const double *const *cmat, const double *dvec, const int verbosity, std::vector &C_red, std::vector &d_red, int &r) -> int; /** - * @brief Solves a constrained least squares problem using sparse matrices. - * - * This function minimizes ||A x - b||_2 subject to C x = d, where: - * - A is the data matrix (M×N) in sparse format. - * - b is the observation vector (length M). - * - C is the constraint matrix (P×N) in sparse format. - * - d is the constraint vector (length P). - * - x is the solution vector (length N). - * - lambda is the vector of Lagrange multipliers (length P). - * - * The method uses QR decomposition with column pivoting to handle the constraints - * and solve the least squares problem efficiently in sparse form. - * - * @param[in] A Sparse matrix A (M×N) representing the data matrix. - * @param[in] b Observation vector b (length M). - * @param[in] C Sparse matrix C (P×N) representing the constraint matrix. - * @param[in] d Constraint vector d (length P). - * @param[out] x Output solution vector x (length N). - * @param[out] lambda Output vector of Lagrange multipliers (length P). - * @param[in] verbosity Verbosity level (0: silent, >0: print info). + * Solve sparse min ||A x - b||_2 subject to C x = d. + * A is M x N and C is P x N; b and d have lengths M and P. + * Return x[N] and Lagrange multipliers lambda[P]. */ auto solveGQRSparse(const Eigen::SparseMatrix &A, const Eigen::VectorXd &b, const Eigen::SparseMatrix &C, const Eigen::VectorXd &d, Eigen::VectorXd &x, @@ -165,18 +72,10 @@ auto solveGQRSparse(const Eigen::SparseMatrix &A, const Eigen::VectorXd const double tolerance_iteration = 1.0e-8, const int maxnum_iteration = 10000) -> void; /** - * Given a dense column-major matrix A (size M×N) stored in - * `A_data` (length M*N), find which rows are independent. - * - * @param M # rows of A - * @param N # cols of A - * @param A_data pointer to column-major data (size M*N) - * @param tol optional tolerance; if ≤0, compute tol = max(M,N)*|R₀₀|*ε - * @param rank [out] computed numerical rank - * @param pivots [out] size-rank array of 0-based row indices in pivot order - * @param verbosity [in] verbosity level (0: silent, >0: print info) - * - * @returns 0 on success, LAPACK INFO otherwise. + * Find independent rows of column-major A_data[M][N]. Return rank and + * zero-based row pivots in pivot order. For tol <= 0, use + * max(M,N)*|R(0,0)|*epsilon. verbosity > 0 prints diagnostics. + * Returns 0 on success or LAPACK INFO. */ auto find_independent_rows_dense(int M, int N, double *A_data, double tol, int &rank, std::vector &pivots, const int verbosity = 0) -> int; @@ -187,8 +86,7 @@ auto find_independent_rows_dense(int M, int N, double *A_data, double tol, int & constexpr double rank_tolerance_auto = -1.0; -/// Extract the independent rows of a (P×N) sparse matrix C_sparse, -/// using LAPACK QR-with-pivoting on Cᵀ, just like your original code. +/// Extract independent rows of C_sparse (P x N) using pivoted QR of C^T. /// /// @param C_sparse Input (P×N), row-major sparse /// @param dvec Input length-P diff --git a/alm/optimize.cpp b/alm/optimize.cpp index c1c894ae..12e167ff 100644 --- a/alm/optimize.cpp +++ b/alm/optimize.cpp @@ -46,13 +46,9 @@ using namespace ALM_NS; namespace { -// Coordinate descent builds the Gram matrix Prod = A^T A column by column, lazily: each column is a -// separate OpenMP-parallel GEMV inside the sweep. When the full Gram is affordable -- N (columns) -// not much larger than M (rows) -- building it up front with one BLAS-3 GEMM is markedly faster than -// that per-column build, even if the elastic-net solution does not use every column (the GEMM's -// efficiency beats the lazy build's per-column overhead). For strongly underdetermined problems -// (N >> M) only a few columns are ever needed, so the lazy build wins. The active set is bounded by -// min(N, M), so N <= factor * M is the affordability test. Set ALM_GRAM_LAZY to force the lazy path. +// Build the Gram matrix with one GEMM when N <= factor * M; use lazy +// column-wise GEMVs for strongly underdetermined problems (N >> M). +// ALM_GRAM_LAZY forces the lazy path. constexpr Eigen::Index gram_dense_factor = 2; inline auto use_full_gram(const Eigen::MatrixXd &A) -> bool @@ -71,11 +67,9 @@ inline auto matvec_nthreads(const Eigen::Index ncols) -> int return static_cast(n); } -// res = A * x, parallelized over the columns of the (column-major) A with per-thread partial sums, -// because Eigen/Accelerate run dgemv single-threaded here. `scratch` is (A.rows() x nthreads) and is -// reused across calls. Columns with x(j)==0 are skipped (exploits a sparse iterate). The team is -// pinned with num_threads(nthreads) so omp_get_thread_num() < nthreads == scratch.cols() always -- -// no out-of-bounds even if the runtime would otherwise pick a different team size. +// Compute A*x over column-major A with per-thread sums in reusable +// scratch[A.rows()][nthreads], skipping zero x entries. Fix the team size +// to match scratch.cols(). inline void parallel_Ax(const Eigen::MatrixXd &A, const Eigen::VectorXd &x, const int nthreads, Eigen::MatrixXd &scratch, Eigen::VectorXd &res) { @@ -101,15 +95,9 @@ inline void parallel_Atr(const Eigen::MatrixXd &A, const Eigen::VectorXd &r, con } } -// OLS solve min_x ||A x - b||_2 for the adaptive-LASSO penalty weights, returning the coefficients -// and the effective column rank in `rank_out`. -// -// Fast path: normal equations (A^T A) x = A^T b solved by Cholesky. A^T A is a threaded BLAS-3 -// product (the same kernel coordinate descent uses to build its Gram), whereas Eigen's -// rank-revealing ColPivHouseholderQR is single-threaded and dominated the adaptive-LASSO runtime. -// If A^T A is not numerically positive-definite (A rank-deficient or severely ill-conditioned), -// fall back to the robust QR. The weights are |x_OLS|, used only to scale the L1 penalty, so the -// squared condition number of the normal-equations path is harmless whenever Cholesky succeeds. +// Compute OLS coefficients and rank_out for adaptive-LASSO weights. +// Use threaded normal equations with Cholesky, falling back to +// rank-revealing QR for rank-deficient or ill-conditioned matrices. inline auto solve_ols_for_adalasso(const Eigen::MatrixXd &A, const Eigen::VectorXd &b, Eigen::Index &rank_out) -> Eigen::VectorXd { @@ -119,10 +107,8 @@ inline auto solve_ols_for_adalasso(const Eigen::MatrixXd &A, const Eigen::Vector Eigen::LLT llt(gram); if (llt.info() == Eigen::Success) { - // Guard against the cond(A)^2 amplification of the normal-equations path: if the Cholesky - // factor has a tiny pivot relative to the largest (near rank deficiency), the OLS weights from - // the squared system can be inaccurate. Fall back to the rank-revealing QR in that case. The - // diagonal of the Cholesky factor is a cheap proxy for conditioning. + // Small relative Cholesky pivots indicate inaccurate normal-equation + // weights; fall back to rank-revealing QR. const auto ldiag = llt.matrixLLT().diagonal().cwiseAbs(); const double dmin = ldiag.minCoeff(); const double dmax = ldiag.maxCoeff(); @@ -265,10 +251,8 @@ auto Optimize::optimize_main(const std::unique_ptr &symmetry, std::uni } if (optcontrol.linear_model == 3) { - // Adaptive LASSO standardizes the (reweighted) design matrix for conditioning, just like - // elastic net. A per-column L1 penalty (= factor_std, set in the solver dispatch) keeps the - // standardized solve equivalent to penalizing the un-standardized reweighted coefficients, - // so STANDARDIZE = 1 is no longer "meaningless" -- it is the default and is much faster. + // Per-column L1 penalties (factor_std) preserve the adaptive-LASSO + // objective when standardizing the reweighted design matrix. if (std::abs(optcontrol.displacement_normalization_factor - 1.0) > eps) { if (verbosity > 0) { @@ -755,11 +739,9 @@ auto Optimize::run_manual_cv(const std::string &job_prefix, const int maxorder, } fnorm_validation = std::sqrt(fnorm_validation); - // Energy term inside CV (EFIT_CV): append the centered/weighted/w-scaled energy rows to both - // the training and validation systems, so each fold's fit AND its held-out error include the - // energy; alpha is then selected by the combined (force + w·energy) relative residual. The - // combined norm sqrt(fnorm^2 + enorm^2) keeps the reported error dimensionless. The pre-augment - // force-row counts/norms are kept so solution_path can also report the separate force/energy errors. + // EFIT_CV adds centered, weighted energy rows to training and validation. + // Select alpha by the combined residual normalized by sqrt(fnorm^2 + enorm^2); + // retain force-row counts and norms for separate component errors. const bool efit_in_cv = (optcontrol.efit_weight > 0.0 && optcontrol.efit_cv); size_t nrow_force_train = 0, nrow_force_val = 0; double fnorm_force_train = 0.0, fnorm_force_val = 0.0, enorm_train_cv = 0.0, enorm_val_cv = 0.0; @@ -1733,10 +1715,8 @@ auto Optimize::solution_path(const int maxorder, Eigen::MatrixXd &A, Eigen::Vect validation_error.push_back(std::sqrt(res2)); nonzeros.push_back(nzero_lasso); - // Optional split of the (combined) residual into force-only and energy-only relative errors. - // All four outputs are filled together (require all pointers). A zero reference norm (no force - // data, or a degenerate all-equal energy set) yields 0.0 — the relative error is then undefined - // but trivially satisfied; this does not occur for a normal energy-in-CV run. + // Fill separate force and energy relative errors when all four outputs + // are provided. Report zero for a zero reference norm. if (terr_force != nullptr && terr_energy != nullptr && verr_force != nullptr && verr_energy != nullptr) { const double fres_t = fdiff.head(nrow_force_train).squaredNorm(); const double eres_t = fdiff.tail(M - nrow_force_train).squaredNorm(); @@ -2212,11 +2192,8 @@ auto Optimize::get_estimated_max_alpha(const Eigen::MatrixXd &Amat, const Eigen: C = C.transpose() * bvec; auto max_alpha = 0.0; - // Adaptive LASSO solves the standardized system with a per-column penalty p_j = factor_std_j = - // 1/dev_j, so the KKT zero-solution threshold for column j is |Z_j^T b| / (p_j * M * l1_ratio), - // i.e. the standardized gradient must be divided by p_j (equivalently multiplied by dev_j). - // Omitting this made the alpha grid start one step too high (a CV-selected-alpha shift). For - // STANDARDIZE = 0, dev == 1 so this reduces to the previous behavior. + // Adaptive-LASSO zero-solution threshold: |Z_j^T b| / (p_j * M * l1_ratio), + // where p_j = factor_std_j = 1/dev_j. With STANDARDIZE = 0, dev_j = 1. const bool adalasso = (optcontrol.linear_model == 3); for (auto i = 0; i < ncols; ++i) { const auto grad_abs = adalasso ? std::abs(C(i)) * dev(i) : std::abs(C(i)); @@ -2877,9 +2854,7 @@ auto Optimize::get_matrix_elements2(const int maxorder, const size_t ncycle, con double **amat_mod_tmp = nullptr; std::vector nonzero_omp; - // Parameter-major (amat[param][component]): the data block is the same size as the old - // [component][param] layout, only the row-pointer array grows to O(ncols) entries per - // thread (a few MB at most) — negligible next to the ncols*natmin3 doubles it indexes. + // Parameter-major layout: amat[param][component]. allocate(amat_orig_tmp, ncols, natmin3); if (constraint->get_constraint_algebraic()) { @@ -3211,10 +3186,8 @@ auto Optimize::fill_amat(const int maxorder, const size_t natmin, const size_t n const std::unique_ptr &symmetry, const std::unique_ptr &fcs, double **&amat_orig) -> void { - // amat_orig is stored parameter-major: amat_orig[iparam][k], i.e. the natmin3 force - // components of a given free parameter are contiguous in memory. This layout lets the - // constraint projection (project_constraints) and the matrix scans operate on contiguous - // length-natmin3 vectors instead of striding column-wise across a row-major buffer. + // Store amat_orig[param][component] so projection and matrix scans + // access contiguous force-component vectors. const auto natmin3 = natmin * 3; for (size_t i = 0; i < ncols; ++i) { @@ -3260,10 +3233,7 @@ auto Optimize::project_constraints(const int maxorder, const size_t natmin, cons const auto natmin3 = 3 * natmin; const auto idata = natmin3 * irow; - // amat_orig / amat_mod are parameter-major (amat[param][component]); the natmin3 force - // components of each parameter are contiguous, so the copies and AXPYs below stream over - // memory instead of striding column-wise across a row-major buffer (the dominant cost of - // this routine for large algebraic-constraint fits). + // Parameter-major buffers keep copies and AXPYs contiguous. for (int order = 0; order < maxorder; ++order) { for (const auto fix: constraint->get_const_fix(order)) { @@ -3274,10 +3244,8 @@ auto Optimize::project_constraints(const int maxorder, const size_t natmin, cons } } - // index_bimap(order).left ranges over [0, index_bimap(order).size()); added to the running - // iparam offset, inew covers [0, ncols_compact) exactly once across all orders. So this loop - // overwrites every row of amat_mod on every training row, and the per-thread amat_mod buffer - // (reused across rows) never carries stale data into the const_relate -= update below. + // The mapping overwrites every amat_mod row before const_relate updates, + // so the reused per-thread buffer needs no clearing. for (const auto &it: constraint->get_index_bimap(order)) { inew = it.left + iparam; iold = it.right + ishift; @@ -3322,10 +3290,8 @@ auto Optimize::project_energy_row(const int maxorder, const std::unique_ptr const std::unique_ptr &constraint, const std::vector &e_full, std::vector &e_compact, double &e_rhs) -> void { - // Single-row analogue of project_constraints for the energy row. e_full is indexed by the full - // symmetry-irreducible parameter index (one per nequiv group); e_compact by the constraint- - // compacted free parameter index. The FC2FIX-fixed-coefficient energy is moved onto e_rhs, - // exactly mirroring the const_fix RHS handling for forces (project_constraints, optimize.cpp). + // Project e_full to the constraint-compacted e_compact, moving fixed + // coefficient energy to e_rhs as project_constraints does for forces. std::fill(e_compact.begin(), e_compact.end(), 0.0); e_rhs = 0.0; @@ -3499,12 +3465,9 @@ auto Optimize::gamma(const int n, const int *arr) const -> double auto Optimize::gamma_energy(const int n, const int *arr) const -> double { - // Per-entry energy multiplicity factor = gamma(n,arr) / n, matching the ALAMODE energy - // convention (cf. tools/taylor.py: E_order = (1/order) * sum gamma * Phi * prod_u). - // The fc_table enumerates all index orderings (so the force builder fill_amat, which writes - // only the lead atom, is complete); the energy therefore needs the extra 1/n to avoid - // overcounting. With this factor the per-order Euler identity - // E_order = -(1/n) * sum_a u_a F_a holds exactly against the force matrix. + // Energy multiplicity is gamma(n, arr) / n (see tools/taylor.py). + // The 1/n corrects for fc_table index orderings and gives + // E_order = -(1/n) * sum_a u_a F_a. return gamma(n, arr) / static_cast(n); } @@ -3541,10 +3504,9 @@ auto Optimize::run_energy_selftest(const std::unique_ptr &symmetry, co const std::unique_ptr &constraint, const int maxorder, const int verbosity) const -> bool { - // Phase-1 verification: the energy-row builder must satisfy the per-order Euler identity - // A_E[c][p] == -(1/n_p) * sum_{supercell atoms a} u_a * A_F[a][p], n_p = order(p) + 2, - // where A_F is the existing (trusted) force matrix. The supercell sum is realized by summing - // the natmin primitive force rows over the ntran translation images (cf. fill_bvec mapping). + // Check the per-order Euler identity against the force matrix: + // A_E[c][p] == -(1/n_p) * sum_a u_a * A_F[a][p], n_p = order(p) + 2. + // Sum over all supercell atoms using the ntran translation images. std::cout << "\n [ALM_ENERGY_SELFTEST] Verifying energy-row builder vs force builder (Euler's theorem)\n"; if (!e_train.empty()) { @@ -3606,19 +3568,14 @@ auto Optimize::run_energy_selftest(const std::unique_ptr &symmetry, co } } - // Arbitrary, varied, deterministic parameter vector. We contract A_E and A_F with theta: - // the assembled force A_F * theta is the COMPLETE physical force on each atom (it sums every - // FC entry leading on that atom), so the per-order Euler identity - // E_order(theta) == -(1/n) * sum_a u_a * F_a^order(theta) - // holds. (A single A_F column carries only the lead atom's share, so a column-wise check - // would spuriously fail by (n-1)/n; the theta contraction is the correct test.) + // Contract with deterministic theta to check the complete per-order force: + // E_order(theta) == -(1/n) * sum_a u_a * F_a^order(theta). + // Individual force columns contain only the lead atom contribution. std::vector theta(ncols); for (size_t p = 0; p < ncols; ++p) theta[p] = 1.5 + std::sin(0.3 * static_cast(p) + 1.0); - // Phase 2: compact-projection validation of project_energy_row. We verify, per config, - // A_E_full . theta_full == A_E_compact . theta_c - e_rhs, theta_full = recover(theta_c), - // which ties the projected energy row to the Phase-1-validated full builder through ALM's own - // parameter-recovery map and exercises the FC2FIX fixed-coefficient energy subtraction (e_rhs). + // Check projection and fixed-coefficient subtraction for each configuration: + // A_E_full . recover(theta_c) == A_E_compact . theta_c - e_rhs. const bool algebraic = constraint->get_constraint_algebraic(); size_t ncols_c = 0; for (int o = 0; o < maxorder; ++o) ncols_c += constraint->get_index_bimap(o).size(); @@ -3726,10 +3683,9 @@ auto Optimize::build_energy_block(const std::unique_ptr &symmetry, con const std::vector &e_in, const double emin, std::vector &amat_out, std::vector &evec_out, double &enorm_out) const -> void { - // Build the constraint-compacted, weighted-Frisch-Waugh-centered, w-scaled energy block for the - // configuration set (u_in, e_in). Used by the production fit (full training set) and by each CV - // fold (train and held-out subsets). `emin` is the per-config-weight reference (pass the global - // training E_min so the weighting is identical across folds). Requires an algebraic constraint. + // Build the constraint-compacted, weighted-centered, w-scaled energy block. + // Requires algebraic constraints; use global training E_min as emin + // to keep weights consistent across CV folds. const size_t n_ene = u_in.size(); amat_out.clear(); evec_out.clear(); @@ -4020,10 +3976,8 @@ auto Optimize::coordinate_descent(const int M, const int N, const double alpha, } } ++iloop; - // Eigen SIMD reduction instead of an every-sweep OpenMP parallel-for: the work is a - // tiny O(N) memory-bound sum, so fork/join dominated it. Also deterministic, whereas - // the OpenMP reduction's sum order (and thus the converged iteration) depended on the - // thread count. + // Use Eigen SIMD for this small sum to avoid OpenMP overhead and + // thread-count-dependent reduction order. diff = std::sqrt(delta.squaredNorm() / static_cast(N)); if (diff < optcontrol.tolerance_iteration) break; @@ -4085,10 +4039,9 @@ auto Optimize::estimate_lipschitz_l2(const Eigen::MatrixXd &A) -> double v.normalize(); Eigen::VectorXd Av(nrows), w(ncols); - // Guaranteed-positive lower bound on lambda_max(A^T A): the largest column norm squared - // (= max diagonal of A^T A <= lambda_max). Used as a floor so the estimate is never zero -- e.g. - // when the all-ones start lies in the null space of A and power iteration would otherwise return - // 0. The starting L only needs to be reasonable; fista()'s backtracking guarantees correctness. + // Floor the power estimate at the largest squared column norm, a lower + // bound on lambda_max(A^T A), in case the start lies in A's null space. + // FISTA backtracking handles underestimation. double col_floor = 0.0; for (Eigen::Index j = 0; j < ncols; ++j) col_floor = std::max(col_floor, A.col(j).squaredNorm()); @@ -4150,14 +4103,9 @@ auto Optimize::fista(const int M, const int N, const double alpha, const int war const auto lambda1 = alpha * optcontrol.l1_ratio; const auto lambda2 = alpha * (1.0 - optcontrol.l1_ratio); - // FISTA with backtracking line search (Beck & Teboulle 2009). L starts from the power-iteration - // estimate and is grown by bt_eta until the proximal step from y satisfies the sufficient-decrease - // (descent) condition f(x_new) <= Q_L(x_new, y). This makes convergence robust to an - // under-estimated Lipschitz constant -- power iteration only ever gives a lower bound -- without - // having to construct a guaranteed upper bound. Caching A*x keeps the cost at two matvecs per - // iteration: the A*x_new built for the descent test also yields the next A*y via - // A*y = A*x_new + momentum*(A*x_new - A*x_cur), - // so no separate A*y matvec is needed. + // FISTA with backtracking (Beck & Teboulle 2009): grow L by bt_eta until + // f(x_new) <= Q_L(x_new, y). Reuse A*x_new for the next extrapolation, + // A*y = A*x_new + momentum*(A*x_new - A*x_cur), keeping two matvecs per iteration. auto L = std::max(lipschitz_l2 + lambda2, eps); constexpr auto bt_eta = 2.0; // L growth factor when a step is rejected constexpr auto bt_max = 60; // hard cap on backtracks per iteration @@ -4293,9 +4241,7 @@ auto Optimize::admm(const int M, const int N, const double alpha, const int warm const auto tol = optcontrol.tolerance_iteration; const auto sqrtN = std::sqrt(static_cast(N)); - // Scaled-ADMM state. The primal z is warm-started from x; the scaled dual u is reset to zero on - // every call (a per-alpha dual warm start gave only a marginal speedup and is not worth the - // extra bookkeeping, so each alpha starts the dual cold). + // Warm-start the primal z from x and reset the scaled dual u for each alpha. Eigen::VectorXd z(N), xa(N), xhat(N), z_old(N), rhs(N); Eigen::VectorXd u = Eigen::VectorXd::Zero(N); Eigen::VectorXd res(M); // only used for the verbosity>1 diagnostic A*z - b diff --git a/alm/optimize.h b/alm/optimize.h index 74456b0e..482c89a5 100644 --- a/alm/optimize.h +++ b/alm/optimize.h @@ -328,21 +328,18 @@ class Optimize const std::unique_ptr &symmetry, const std::unique_ptr &fcs, double **&amat_orig) -> void; - // Phase 1 (energy term): build a single energy row for one displacement image. - // energy_row[iparam] += gamma_energy_precomputed * prod_{j=0..order+1} u_sub[elems[j]] - // (product over ALL order+2 indices; cf. fill_amat which drops elems[0] for the force). + // Accumulate one energy row using the product over all order+2 indices; + // fill_amat omits the first index for forces. static auto fill_amat_energy(const int maxorder, const size_t ncols, const std::vector &u_sub, const std::vector> &gamma_energy_precomputed, const std::unique_ptr &fcs, std::vector &energy_row) -> void; - // Per-entry energy multiplicity factor = gamma(n,arr)/n (NOT 1/denom; the two coincide only - // for diagonal clusters). Matches tools/taylor.py (E_order *= 1/order). Caller multiplies by - // fc_table.sign, exactly as for the force gamma table. + // Energy multiplicity is gamma(n, arr)/n, matching tools/taylor.py. + // The caller applies fc_table.sign as for the force gamma table. [[nodiscard]] auto gamma_energy(const int n, const int *arr) const -> double; - // Self-contained verification (env ALM_ENERGY_SELFTEST): checks the energy-row builder - // against the trusted force builder via the per-order Euler identity - // E_order == -(1/n) * sum_a u_a F_a(theta). Returns true on PASS. Does not touch the fit path. + // ALM_ENERGY_SELFTEST checks E_order = -(1/n) * sum_a u_a F_a(theta). + // Returns true on success; leaves the fit unchanged. auto run_energy_selftest(const std::unique_ptr &symmetry, const std::unique_ptr &fcs, const std::unique_ptr &constraint, const int maxorder, const int verbosity) const -> bool; @@ -351,9 +348,8 @@ class Optimize const std::unique_ptr &fcs, const std::unique_ptr &constraint, double **amat_orig, double **&amat_mod, std::vector &bvec_mod) -> void; - // Phase 2: project one full-basis energy row (e_full, length ncols) into the constraint-compacted - // basis (e_compact, length ncols_compact), moving the FC2FIX-fixed-coefficient energy onto the - // RHS scalar e_rhs. Single-row analogue of project_constraints' const_fix/index_bimap/const_relate. + // Project e_full[ncols] to e_compact[ncols_compact], moving the fixed + // FC2FIX energy to e_rhs, as project_constraints does for forces. static auto project_energy_row(const int maxorder, const std::unique_ptr &fcs, const std::unique_ptr &constraint, const std::vector &e_full, std::vector &e_compact, double &e_rhs) -> void; diff --git a/alm/rref.cpp b/alm/rref.cpp index a6a5f0d4..487285d5 100644 --- a/alm/rref.cpp +++ b/alm/rref.cpp @@ -238,14 +238,9 @@ auto rref_sparse(const size_t ncols, ConstraintSparseForm &sp_constraint, const auto rref_sparse_pivot(const size_t ncols, ConstraintSparseForm &sp_constraint, const double tolerance) -> void { - // Same Gauss-Jordan elimination as rref_sparse(), but with PARTIAL (maximum-magnitude) row - // pivoting. For each pivot column we move the row with the largest |entry| in that column to - // the pivot position before normalizing, instead of accepting the first entry above the - // tolerance. The set of pivot columns (the left-to-right linearly independent columns) is - // unchanged, so the reduced echelon form -- and therefore the const_fix / const_relate / - // index_bimap map derived from it -- is mathematically identical to rref_sparse() up to - // round-off, but never divides by a small accepted pivot. This is the coordinate-preserving, - // numerically stable replacement (Policy A) for rref_sparse() in the algebraic path. + // Gauss-Jordan elimination with maximum-magnitude row pivoting. + // Preserves the left-to-right pivot-column order used by the constraint map + // while reducing round-off compared with first-acceptable-pivot RREF. const auto nrows = sp_constraint.size(); if (nrows == 0) return; diff --git a/alm/rref.h b/alm/rref.h index e996796c..54cda62d 100644 --- a/alm/rref.h +++ b/alm/rref.h @@ -9,10 +9,8 @@ auto rref(std::vector> &mat, const double tolerance = 1.0e-1 auto rref_sparse(const size_t ncols, ConstraintSparseForm &sp_constraint, const double tolerance = 1.0e-12) -> void; -// Reduced row echelon form with partial (maximum-magnitude) row pivoting. -// Produces the same echelon structure as rref_sparse (identical pivot columns and reduced -// coefficients up to round-off) but is numerically more stable, because it never divides by a -// small accepted pivot. The output is consumed unchanged by Constraint::get_mapping_constraint. +// Partial-pivot RREF preserving the pivot-column order used by +// Constraint::get_mapping_constraint, with improved numerical stability. auto rref_sparse_pivot(const size_t ncols, ConstraintSparseForm &sp_constraint, const double tolerance = 1.0e-12) -> void; diff --git a/alm/symmetry.cpp b/alm/symmetry.cpp index 42c8f62f..5585f130 100644 --- a/alm/symmetry.cpp +++ b/alm/symmetry.cpp @@ -209,10 +209,8 @@ auto Symmetry::setup_symmetry_operation(const Cell &pcell, const Spin &spin_prim symmetry_data_super.clear(); symmetry_data_prim.clear(); - // First, generate space group operations using the primitive cell. - // Please be noted that the input pcell might not be a true primitive cell - // because one can give PRIMCELL value which does not necessary transform the - // input cell into a true primitive cell. + // Generate space-group operations from pcell, which may be non-primitive + // if the user-supplied PRIMCELL does not fully reduce the input cell. if (spin_prim.lspin && spin_prim.noncollinear) { // if (spg_get_major_version() < 2) { // std::cout << " Please update spglib to version 2 or above to .\n"; diff --git a/alm/symmetry.h b/alm/symmetry.h index 93f34c85..7f4fa491 100644 --- a/alm/symmetry.h +++ b/alm/symmetry.h @@ -152,12 +152,8 @@ class Symmetry size_t nsym_super, ntran_super; size_t nsym_prim, ntran_prim; size_t nat_trueprim; - // nat_trueprim is the number of atoms included in a true primitive cell. - // This value can be different from primcell.number_of_atoms because - // the latter value is calculated from the PRIMCELL value given by users, - // which does not necessary reduces the inputcell to a true primitive cell. - // When nat_trueprim != primcell.number_of_atoms, ntran_prim will be larger - // than 1. + // nat_trueprim counts atoms in the true primitive cell. A user-defined + // PRIMCELL may contain more atoms, in which case ntran_prim > 1. std::vector> map_fullsymmetry_super; // [nat_base, nsym_super] std::vector> map_fullsymmetry_prim; // [nat_base, nsym_super] diff --git a/alm/system.cpp b/alm/system.cpp index a39fdb06..d4350b09 100644 --- a/alm/system.cpp +++ b/alm/system.cpp @@ -479,10 +479,8 @@ auto System::set_default_variables() -> void spin_input.lspin = false; spin_input.noncollinear = 0; spin_input.time_reversal_symm = 1; - // Same default as Symmetry::set_default_variables and the CLI TOLERANCE tag. - // Without it, library users (e.g. the Python interface) that never call - // set_tolerance passed an uninitialized value to spglib, whose failure - // path then crashed (NULL free in det_determine_all). + // Default to the same tolerance as Symmetry and the CLI so library calls + // without set_tolerance still pass a valid value to spglib. symmetry_tolerance = 1.0e-3; autoset_primcell = 0; transmat_to_super = Eigen::Matrix3d::Identity(); @@ -543,12 +541,8 @@ auto System::get_atomtype_group(const std::string &cell) const -> const std::vec auto System::set_atomtype_group(const Cell &cell_in, const Spin &spin_in, std::vector> &atomtype_group_out) -> void { - // In the case of collinear calculation, spin moments are considered as scalar - // variables. Therefore, the same elements with different magnetic moments are - // considered as different types. In noncollinear calculations, - // magnetic moments are not considered in this stage. They will be treated - // separately in symmetry.cpp where spin moments will be rotated and flipped - // using time-reversal symmetry. + // Collinear moments distinguish atom types. Noncollinear moments are + // handled in symmetry.cpp with spin rotations and time reversal. unsigned int i; AtomType type_tmp{}; diff --git a/anphon/anharmonic_core.h b/anphon/anharmonic_core.h index 4448d7fd..95292e1b 100644 --- a/anphon/anharmonic_core.h +++ b/anphon/anharmonic_core.h @@ -149,12 +149,8 @@ class AnharmonicCore: protected Pointers const std::complex *const *const *evec_in, const PhaseFactorCache *phase_storage_in); - // Thread-safe serial variants: the per-triplet reciprocal-FC3 cache - // lives in caller-provided storage (phi3_work of size ngroup_v3 and - // kindex_work[2] initialized to -1) and no OpenMP region is entered. - // Parallelism belongs to the caller's triplet loop; entering a parallel - // region per V3 call is far too fine-grained (measured to give negative - // scaling in the SERTA and IBTE setups). + // Serial kernels with caller-owned caches: phi3_work[ngroup_v3] and + // kindex_work[2], initially -1. Parallelize the caller's triplet loop. std::complex V3(const unsigned int ks[3], const double *const *xk_in, const double *const *eval_in, const std::complex *const *const *evec_in, std::complex *phi3_work, int *kindex_work); diff --git a/anphon/collision_operator.cpp b/anphon/collision_operator.cpp index 0fb0984a..a594db4d 100644 --- a/anphon/collision_operator.cpp +++ b/anphon/collision_operator.cpp @@ -386,12 +386,8 @@ void CollisionOperator::setup_L_smear() anharmonic_core_.prepare_fc3_compressed(); - // The loops run over the flattened triplet index (pairs_emitt/absorb, - // built in get_triplets). The factorized V3 kernel gives |V3|^2 of all - // ns^3 band combinations of a triplet at once; the static schedule mostly - // keeps the triplets of one k point on one thread so that its psi_K is - // reused (a group may straddle a chunk boundary, which only costs one - // extra fold). + // Process flattened triplets with all ns^3 band combinations at once. + // Static scheduling improves reuse of psi_K for triplets sharing k. #ifdef _OPENMP #pragma omp parallel diff --git a/anphon/collision_operator.h b/anphon/collision_operator.h index 99d9dbb7..0186c2a6 100644 --- a/anphon/collision_operator.h +++ b/anphon/collision_operator.h @@ -24,17 +24,11 @@ class Integration; class AnharmonicCore; class TetraNodes; class DymatEigenValue; -// Distributed three-phonon collision operator on the irreducible wedge, -// shared by the iterative-family BTE solvers (SOLVER = IBTE today; a -// variational/CG solver can reuse it later). The irreducible k points are -// distributed round-robin over the MPI ranks; L_emitt/L_absorb store the -// transition probabilities of the local rows, and a symmetry expansion -// table maps wedge values of Cartesian vector fields back onto the full -// grid. Diagonal add-ons beyond 3ph (isotope, boundary, 4ph) remain the -// solver's responsibility. -// All dependencies are explicit constructor arguments (no Pointers base): -// the operator is constructed by Iterativebte::setup, after setup_base(), -// when every input already exists. +// Three-phonon collision operator on the irreducible wedge. MPI ranks own +// round-robin k points and store local transition probabilities in +// L_emitt/L_absorb; symmetry expands vector fields to the full grid. +// The solver handles diagonal add-ons (isotope, boundary, 4ph). +// Constructed by Iterativebte::setup after setup_base(). class CollisionOperator { public: @@ -53,12 +47,9 @@ class CollisionOperator // following integration->ismear). void build_L(); - // Include the elastic isotope-disorder channel (Tamura kernel) in the - // operator: its in-scattering enters calc_W_at and its diagonal is the - // row sum, so the channel conserves the constant mode exactly. Must be - // set before build_L(). isotope_factor_in points at the per-species - // mass-variance factors (Isotope::isotope_factor) and must stay alive - // for the lifetime of this object. + // Enable Tamura isotope in-scattering and its row-sum diagonal, preserving + // the constant mode. Call before build_L(); per-species isotope_factor_in + // must remain valid for this object's lifetime. void set_isotope_channel(const bool flag, const double *isotope_factor_in = nullptr) { with_isotope = flag; diff --git a/anphon/conductivity.cpp b/anphon/conductivity.cpp index 1b38a112..943313b9 100644 --- a/anphon/conductivity.cpp +++ b/anphon/conductivity.cpp @@ -14,6 +14,8 @@ #include #include #include +#include +#include #include #include #include "anharmonic_core.h" @@ -40,6 +42,13 @@ using namespace PHON_NS; +// File-static helpers defined further down; declared here because they are used +// earlier in the file (formulation stamp at result-file open, block report). +static std::string active_transport_formulation(bool nonanalytic); +static void build_block_table(const KpointMeshUniform *kmesh_in, const double *const *eval_in, unsigned int ns, + std::vector> &lo_out, std::vector> &hi_out); + + Conductivity::Conductivity(PHON *phon) : Pointers(phon) { set_default_variables(); @@ -102,6 +111,9 @@ void Conductivity::deallocate_variables() if (velmat) { velmat.clear(); } + if (velblock) { + velblock.clear(); + } if (vel_4ph) { vel_4ph.clear(); } @@ -182,14 +194,24 @@ void Conductivity::setup_kappa() Bohr_in_Angstrom * 1.0e-10 / time_ry, false); + // The full velocity matrix is the memory hog (nk ns^2 x 3 complex) and is needed only + // by the coherent term. The default block-trace Peierls term and the boundary speed + // use velblock, the per-branch block-summed diad (nk ns x 9 doubles), which is + // contracted per k point on the fly so the full matrix is never stored unless asked. + const auto corrected = !PhononVelocity::legacy_velocity(); if (calc_coherent) { - if (mympi->my_rank == 0) { - velmat.resize(nk_3ph, ns, ns, 3); - } else { + if (mympi->my_rank == 0) velmat.resize(nk_3ph, ns, ns, 3); + else velmat.resize(1, 1, 1, 3); - } - phonon_velocity->calc_phonon_velmat_mesh(velmat); - check_velocity_matrix_consistency(dos->kmesh_dos.get(), dos->dymat_dos->get_eigenvalues()); + } + if (corrected) { + if (mympi->my_rank == 0) velblock.resize(nk_3ph, ns, 3, 3); + else + velblock.resize(1, 1, 1, 1); + } + if (calc_coherent || corrected) { + phonon_velocity->calc_phonon_velmat_mesh(calc_coherent ? &velmat : nullptr, corrected ? &velblock : nullptr); + if (calc_coherent) check_velocity_matrix_consistency(dos->kmesh_dos.get(), dos->dymat_dos->get_eigenvalues()); if (calc_coherent == 2) { file_coherent_elems = phon->job_title + ".kc_elem"; } @@ -527,6 +549,7 @@ void Conductivity::setup_result_io(const int mode) if (!result_io_h5) { result_io_h5 = std::make_unique(file_kappa_h5); } + result_io_h5->transport_formulation = active_transport_formulation(dynamical->nonanalytic != 0); result_io_h5->open_or_create(build_kappa_file_meta(), build_kappa_channel_meta(mode), !restart_flag_3ph); @@ -536,6 +559,7 @@ void Conductivity::setup_result_io(const int mode) // setup_kappa_4ph() without the 3ph setup path having // created the file first. result_io_h5 = std::make_unique(file_kappa_h5); + result_io_h5->transport_formulation = active_transport_formulation(dynamical->nonanalytic != 0); result_io_h5->open_or_create(build_kappa_file_meta(), build_kappa_channel_meta(mode), !restart_flag_4ph); @@ -708,6 +732,19 @@ KappaChannelMetaH5 Conductivity::build_kappa_channel_meta(const int mode) const } } + if (mode == 1 && !velblock.empty()) { + meta.velocity_diad.reserve(nequiv_total * ns * 9); + for (auto i = 0; i < meta.nk_irred; ++i) { + for (const auto &kp: kmesh_in->kpoint_irred_all[i]) { + for (auto is = 0; is < ns; ++is) { + for (auto a = 0; a < 3; ++a) { + for (auto b = 0; b < 3; ++b) meta.velocity_diad.push_back(velblock[kp.knum][is][a][b]); + } + } + } + } + } + meta.velocities.reserve(nequiv_total * ns * 3); for (auto i = 0; i < meta.nk_irred; ++i) { for (const auto &kp: kmesh_in->kpoint_irred_all[i]) { @@ -1114,6 +1151,93 @@ void Conductivity::write_result_gamma(const unsigned int ik, const unsigned int } +// Names the transport formulation that actually ran, for the result metadata. Reports +// behaviour rather than requested switches: a switch whose prerequisites were not met +// must not be advertised. +static std::string active_transport_formulation(const bool nonanalytic) +{ + if (PhononVelocity::legacy_velocity()) return "legacy"; + std::string out = "velmat_blocktrace_nosym"; + if (nonanalytic) out += ",nonanalytic"; + return out; +} + +// Error policy for the numerical degeneracy tolerance (see degeneracy_utils.h): a merged +// pair with true splitting dw and summed HWHM G has its band-like weight overestimated by +// 1 + (dw/G)^2. The tolerance cannot tell such a pair from a degenerate one, so instead +// of silently picking a limit the worst dw/G among merged blocks is reported. +void Conductivity::report_unresolved_degenerate_blocks(const KpointMeshUniform *kmesh_in, const double *const *eval_in, + const double *const *gamma_in) const +{ + if (mympi->my_rank != 0 || PhononVelocity::legacy_velocity() || writes->getVerbosity() == 0) return; + + std::vector> lo, hi; + build_block_table(kmesh_in, eval_in, ns, lo, hi); + + auto nblocks = 0, nbad = 0, worst_ik = 0, worst_lo = 0, worst_hi = 0; + auto worst = 0.0; + for (auto ik = 0; ik < kmesh_in->nk_irred; ++ik) { + const auto knum = kmesh_in->kpoint_irred_all[ik][0].knum; + for (auto is = 0; is < static_cast(ns);) { + const auto d = hi[ik][is] - lo[ik][is]; + if (d > 1 && eval_in[knum][is] >= eps8) { + ++nblocks; + const auto dw = in_kayser(eval_in[knum][hi[ik][is] - 1]) - in_kayser(eval_in[knum][is]); + // smallest summed HWHM over any pair in the block and any temperature + auto gmin = std::numeric_limits::max(); + for (auto it = 0; it < ntemp; ++it) { + for (auto a = lo[ik][is]; a < hi[ik][is]; ++a) { + for (auto b = a + 1; b < hi[ik][is]; ++b) { + gmin = std::min(gmin, in_kayser(gamma_in[ik * ns + a][it] + gamma_in[ik * ns + b][it])); + } + } + } + if (gmin > 0.0) { + const auto r = dw / gmin; + if (r > worst) { + worst = r; + worst_ik = ik; + worst_lo = lo[ik][is]; + worst_hi = hi[ik][is]; + } + if (r > 0.1) ++nbad; + } + } + is = hi[ik][is]; + } + } + if (nbad > 0) { + const auto flags = std::cout.flags(); + const auto prec = std::cout.precision(); + std::cout << "\n WARNING: " << nbad << " of " << nblocks + << " degenerate transport blocks have a splitting exceeding 0.1 x their summed linewidth\n" + << " (worst dw/Gamma = " << std::scientific << std::setprecision(2) << worst + << " at irreducible k " << worst_ik + 1 << ", branches " << worst_lo + 1 << "-" << worst_hi + << "). Their band-like weight overestimates the coherent limit by up to 1 + (dw/Gamma)^2.\n" + << " These pairs are treated as degenerate by the configured numerical criterion (" + << transport_block_tol_cm() << " cm^-1); treat kappa from these modes with care.\n\n"; + std::cout.flags(flags); + std::cout.precision(prec); + } +} + +// Blocks of every irreducible k, indexed [ik][branch]. Frequency based only, hence +// temperature independent. +static void build_block_table(const KpointMeshUniform *kmesh_in, const double *const *eval_in, const unsigned int ns, + std::vector> &lo_out, std::vector> &hi_out) +{ + const auto nk_irred = kmesh_in->nk_irred; + const auto tol_cm = transport_block_tol_cm(); + + lo_out.resize(nk_irred); + hi_out.resize(nk_irred); + + for (auto ik = 0; ik < nk_irred; ++ik) { + const auto knum = kmesh_in->kpoint_irred_all[ik][0].knum; + transport_block_bounds(ns, eval_in[knum], tol_cm, lo_out[ik], hi_out[ik]); + } +} + void Conductivity::compute_kappa() { unsigned int i; @@ -1143,14 +1267,38 @@ void Conductivity::compute_kappa() double vel_norm; if (len_boundary > eps) { + // Use the basis-invariant block speed for boundary scattering: + // |v|^2 = (1/d) sum_{j,j' in B} sum_mu V^mu_{jj'} V^mu_{j'j}. + // It is constant within each block, as required by the block-trace weights, + // and reduces to the ordinary speed for non-degenerate branches. + std::vector> bnd_lo, bnd_hi; + const auto use_block_speed = !PhononVelocity::legacy_velocity(); + if (use_block_speed) { + build_block_table(dos->kmesh_dos.get(), dos->dymat_dos->get_eigenvalues(), ns, bnd_lo, bnd_hi); + } + for (iks = 0; iks < dos->kmesh_dos->nk_irred * ns; ++iks) { vel_norm = 0.0; auto knum = dos->kmesh_dos->kpoint_irred_all[iks / ns][0].knum; auto snum = iks % ns; - for (auto j = 0; j < 3; ++j) { - vel_norm += vel[knum][snum][j] * vel[knum][snum][j]; + + if (use_block_speed) { + const auto lo = bnd_lo[iks / ns][snum]; + const auto hi = bnd_hi[iks / ns][snum]; + // (1/d) Tr(P V^a P V^a P): velblock is already summed over the second + // block index, so summing its trace over the block's branches gives + // the full double sum. + for (auto is2 = lo; is2 < hi; ++is2) { + for (auto a = 0; a < 3; ++a) vel_norm += velblock[knum][is2][a][a]; + } + vel_norm /= static_cast(hi - lo); + vel_norm = std::sqrt(std::max(vel_norm, 0.0)); + } else { + for (auto j = 0; j < 3; ++j) { + vel_norm += vel[knum][snum][j] * vel[knum][snum][j]; + } + vel_norm = std::sqrt(vel_norm); // legacy: unchanged expression } - vel_norm = std::sqrt(vel_norm); for (i = 0; i < ntemp; ++i) { gamma_total[iks][i] += (vel_norm / len_boundary) * time_ry; // same unit as gamma @@ -1219,6 +1367,7 @@ void Conductivity::compute_kappa() } lifetime_from_gamma(gamma_total, lifetime); + report_unresolved_degenerate_blocks(dos->kmesh_dos.get(), dos->dymat_dos->get_eigenvalues(), gamma_total); kappa.resize(ntemp, 3, 3); @@ -1236,6 +1385,7 @@ void Conductivity::compute_kappa() if (use_h5_io) { + result_io_h5->transport_formulation = active_transport_formulation(dynamical->nonanalytic != 0); result_io_h5->store_kappa(kappa, fph_rta > 0 ? kappa_3only.ptr() : nullptr, calc_coherent ? kappa_coherent.ptr() : nullptr, @@ -1272,6 +1422,8 @@ void Conductivity::compute_kappa_intraband(const KpointMeshUniform *kmesh_in, co const auto nk_irred = kmesh_in->nk_irred; kappa_mode.resize(ntemp, 9, ns, nk_irred); + const auto use_blocktrace = !PhononVelocity::legacy_velocity(); + for (i = 0; i < ntemp; ++i) { for (unsigned int j = 0; j < 3; ++j) { for (unsigned int k = 0; k < 3; ++k) { @@ -1291,10 +1443,22 @@ void Conductivity::compute_kappa_intraband(const KpointMeshUniform *kmesh_in, co auto vv_tmp = 0.0; const auto nk_equiv = kmesh_in->kpoint_irred_all[ik].size(); - // Accumulate group velocity (diad product) for the reducible k points + // Accumulate group velocity (diad product) for the reducible k points. + // Default: degenerate-block trace of the velocity matrix (see below). + // Legacy opt-out: product of finite-difference velocities. for (auto ieq = 0; ieq < nk_equiv; ++ieq) { const auto ktmp = kmesh_in->kpoint_irred_all[ik][ieq].knum; - vv_tmp += vel[ktmp][is][j] * vel[ktmp][is][k]; + + if (use_blocktrace) { + // Block-summed diad. Summed over the branches of one + // degenerate block this is Tr(P_D V^j P_D V^k P_D), + // invariant under any rotation inside the block; the + // matching same-block pairs are removed from the + // coherent term so nothing is double counted. + vv_tmp += velblock[ktmp][is][j][k]; + } else { + vv_tmp += vel[ktmp][is][j] * vel[ktmp][is][k]; + } } if (thermodynamics->classical) { @@ -1353,6 +1517,11 @@ void Conductivity::compute_kappa_coherent(const KpointMeshUniform *kmesh_in, con const auto nk_irred = kmesh_in->nk_irred; + // With the block-trace Peierls term, pairs inside one degenerate block are + // already counted there; only genuine cross-block pairs stay wave-like. + const auto use_blocktrace = !PhononVelocity::legacy_velocity(); + std::vector> block_lo, block_hi; + std::ofstream ofs; if (calc_coherent == 2) { ofs.open(file_coherent_elems.c_str(), std::ios::out); @@ -1362,6 +1531,8 @@ void Conductivity::compute_kappa_coherent(const KpointMeshUniform *kmesh_in, con kappa_save.resize(ns2, nk_irred); } + if (use_blocktrace) build_block_table(kmesh_in, eval_in, ns, block_lo, block_hi); + for (auto i = 0; i < ntemp; ++i) { for (unsigned int j = 0; j < 3; ++j) { for (unsigned int k = 0; k < 3; ++k) { @@ -1383,6 +1554,14 @@ void Conductivity::compute_kappa_coherent(const KpointMeshUniform *kmesh_in, con const auto omega2 = eval_in[knum][js]; if (omega1 < eps8 || omega2 < eps8) continue; + + // same degenerate block -> already in the Peierls trace. + // Zero the element record too, otherwise it keeps whatever a + // previous temperature left there. + if (use_blocktrace && js >= block_lo[ik][is] && js < block_hi[ik][is]) { + if (calc_coherent == 2 && j == k) kappa_save[ib][ik] = czero; + continue; + } auto vv_tmp = czero; const auto nk_equiv = kmesh_in->kpoint_irred_all[ik].size(); @@ -1538,6 +1717,67 @@ void Conductivity::check_velocity_matrix_consistency(const KpointMeshUniform *km << velmat[max_herm_k][max_herm_j][max_herm_i][max_herm_mu].imag() << '\n'; ofs.close(); + // ALAMODE_CHECK_VELMAT=full compares finite differences of sorted + // eigenvalues with analytic velocity-matrix diagonals. Band crossings and + // little-group symmetrization can cause differences. dw_min is the nearest + // branch gap; ALAMODE_VELMAT_GAPMAX filters the dump by that gap. + const auto *const dump_mode = std::getenv("ALAMODE_CHECK_VELMAT"); + + if (dump_mode && std::string(dump_mode) == "full") { + auto gap_max = -1.0; // negative => keep every mode + if (const auto *const gap_env = std::getenv("ALAMODE_VELMAT_GAPMAX")) { + gap_max = std::atof(gap_env); + } + + const auto dumpname = phon->job_title + ".velmat_dump"; + std::ofstream ofs_dump(dumpname.c_str(), std::ios::out); + if (!ofs_dump) exit("check_velocity_matrix_consistency", "Could not open velmat_dump file"); + + ofs_dump << "# Per-mode group velocity: finite difference vs analytic velocity-matrix diagonal\n"; + ofs_dump << "# vFD : central difference of sorted eigenvalues, no symmetrization at k\n"; + ofs_dump << "# reVmat : Re of the analytic velocity-matrix diagonal as used by the transport terms\n" + "# (unsymmetrized by default; little-group symmetrized under ALAMODE_LEGACY_VELOCITY)\n"; + ofs_dump << "# velocities in m/s; omega and dw_min in cm^-1; xk in fractional coordinates\n"; + if (gap_max > 0.0) { + ofs_dump << "# restricted to modes with dw_min < " << gap_max << " cm^-1\n"; + } + ofs_dump << "# ik kx ky kz branch omega dw_min" + " vFD_x vFD_y vFD_z reVmat_x reVmat_y reVmat_z imVmat_x imVmat_y imVmat_z\n"; + ofs_dump << std::scientific << std::setprecision(10); + + auto ndump = 0ULL; + + for (auto ik = 0u; ik < kmesh_in->nk; ++ik) { + for (auto is = 0u; is < ns; ++is) { + if (eval_in[ik][is] < eps8) continue; + + const auto omega = in_kayser(eval_in[ik][is]); + auto dw_min = std::numeric_limits::max(); + + for (auto js = 0u; js < ns; ++js) { + if (js == is || eval_in[ik][js] < eps8) continue; + dw_min = std::min(dw_min, std::abs(in_kayser(eval_in[ik][js]) - omega)); + } + + if (gap_max > 0.0 && dw_min > gap_max) continue; + + ofs_dump << ik + 1; + for (auto mu = 0; mu < 3; ++mu) ofs_dump << ' ' << kmesh_in->xk[ik][mu]; + ofs_dump << ' ' << is + 1 << ' ' << omega << ' ' << dw_min; + for (auto mu = 0u; mu < 3; ++mu) ofs_dump << ' ' << vel[ik][is][mu]; + for (auto mu = 0u; mu < 3; ++mu) ofs_dump << ' ' << velmat[ik][is][is][mu].real(); + for (auto mu = 0u; mu < 3; ++mu) ofs_dump << ' ' << velmat[ik][is][is][mu].imag(); + ofs_dump << '\n'; + ++ndump; + } + } + ofs_dump.close(); + + if (writes->getVerbosity() > 0) { + std::cout << " Per-mode velocity dump (" << ndump << " modes) written to " << dumpname << '\n'; + } + } + if (writes->getVerbosity() > 0) { const auto flags = std::cout.flags(); const auto precision = std::cout.precision(); diff --git a/anphon/conductivity.h b/anphon/conductivity.h index e1c8031a..e7de2040 100644 --- a/anphon/conductivity.h +++ b/anphon/conductivity.h @@ -98,6 +98,8 @@ class Conductivity: protected Pointers NDArray vel, vel_4ph; NDArray, 4> velmat; + // Per-branch block-summed velocity diad [nk][ns][3][3]; see calc_phonon_velmat_mesh. + NDArray velblock; unsigned int nk_3ph, ns; int nshift_restart, nshift_restart4; std::vector vks_l, vks_done, vks_done4; @@ -160,6 +162,9 @@ class Conductivity: protected Pointers void check_velocity_matrix_consistency(const KpointMeshUniform *kmesh_in, const double *const *eval_in) const; + void report_unresolved_degenerate_blocks(const KpointMeshUniform *kmesh_in, const double *const *eval_in, + const double *const *gamma_in) const; + void interpolate_data(const KpointMeshUniform *kmesh_coarse_in, const KpointMeshUniform *kmesh_dense_in, const double *const *val_coarse_in, double **val_dense_out) const; }; diff --git a/anphon/degeneracy_utils.h b/anphon/degeneracy_utils.h index 773416ff..981e94f7 100644 --- a/anphon/degeneracy_utils.h +++ b/anphon/degeneracy_utils.h @@ -12,6 +12,7 @@ #include #include +#include "constants.h" namespace PHON_NS { @@ -68,4 +69,46 @@ inline void average_over_degenerate_modes(const int ns, const double *eval_at_k, is += ideg_now; } } + +// Transport blocks use the Peierls limit Tr(P V^u P V^v P) for nearly +// degenerate modes. TOL_CM is an empirical threshold and can merge +// numerically resolvable splittings. For splitting dw and summed HWHM G, +// this overestimates the pair weight by 1 + (dw/G)^2; compute_kappa +// reports the worst dw/G. Blocks subdivide the damping-averaging groups +// to keep lifetimes constant; damping values do not determine membership. +inline double transport_block_tol_cm() +{ + return 1.0e-6; +} + +// Return block bounds [lo, hi) for each branch. Subdivide the frequency +// groups from find_degenerate_groups (1e-7 Ry) using distance from each +// block's first frequency, preventing chains of close modes from merging. +// Staying within damping-averaging groups ensures constant block lifetimes. +inline void transport_block_bounds(const unsigned int ns, const double *eval_at_k, const double tol_cm, + std::vector &lo_out, std::vector &hi_out) +{ + lo_out.assign(ns, 0); + hi_out.assign(ns, 0); + + std::vector freq_groups; + find_degenerate_groups(ns, eval_at_k, freq_groups); + + auto gbegin = 0u; + for (const auto ndeg: freq_groups) { + const auto gend = gbegin + static_cast(ndeg); + auto is = gbegin; + while (is < gend) { + const auto anchor = in_kayser(eval_at_k[is]); + auto hi = is + 1; + while (hi < gend && std::abs(in_kayser(eval_at_k[hi]) - anchor) < tol_cm) ++hi; + for (auto k = is; k < hi; ++k) { + lo_out[k] = static_cast(is); + hi_out[k] = static_cast(hi); + } + is = hi; + } + gbegin = gend; + } +} } // namespace PHON_NS diff --git a/anphon/dense_hermitian_eigen.h b/anphon/dense_hermitian_eigen.h index 8c597b1c..677094e4 100644 --- a/anphon/dense_hermitian_eigen.h +++ b/anphon/dense_hermitian_eigen.h @@ -15,18 +15,10 @@ namespace PHON_NS { -// Backend seam for dense Hermitian eigenproblems (used by dynamical-matrix -// diagonalization). -// -// v1 backend: LAPACK zheev on the calling rank. Candidate later backends -// behind the same call: ELPA and MAGMA. -// -// mat_in is row-pointer indexed [i][j] and is copied into column-major scratch -// internally. LWORK is fixed at (2n-1)*10 with no workspace query because zheev -// may take a different (blocked vs unblocked) path for different workspace -// sizes and bit-identical results with the historical call are required. -// compute_evec drives JOBZ ('V'/'N'); evec_out (nullable) independently gates -// the eigenvector write-back. +// Dense Hermitian eigensolver using LAPACK zheev on the calling rank. +// Copy row-pointer mat_in[i][j] to column-major scratch. Keep +// LWORK = (2n-1)*10 to preserve the historical LAPACK path and results. +// compute_evec selects JOBZ; nullable evec_out controls write-back. void solve_dense_hermitian(int n, const std::complex *const *mat_in, double *eval_out, std::complex **evec_out, bool compute_evec, char uplo = 'U'); @@ -36,12 +28,8 @@ void solve_dense_hermitian(int n, const std::complex *const *mat_in, dou int solve_dense_hermitian_info(int n, const std::complex *const *mat_in, double *eval_out, std::complex **evec_out, bool compute_evec, char uplo = 'U'); -// Divide-and-conquer variant (LAPACK zheevd, workspace queried) on Eigen -// matrices, for the SCPH solver: eigenvalues ascending in eval_out; the -// eigenvectors, when evec_out is given, in its columns (the convention of -// Eigen::SelfAdjointEigenSolver::eigenvectors()). zheevd is several times -// faster than Eigen's tridiagonal QR at n of a few hundred and uses the -// threaded MKL/OpenBLAS kernels; results agree to roundoff (eigenvector phases -// may differ, which the SCPH solver is invariant to). +// LAPACK zheevd for SCPH Eigen matrices, with queried workspace. +// Return ascending eigenvalues and, when requested, eigenvectors in columns. +// Eigenvector phases may differ from Eigen's solver. void solve_dense_hermitian_dc(const Eigen::MatrixXcd &mat_in, Eigen::VectorXd &eval_out, Eigen::MatrixXcd *evec_out); } // namespace PHON_NS diff --git a/anphon/dense_symmetric_eigen.h b/anphon/dense_symmetric_eigen.h index 319c61c5..5fc8a2fc 100644 --- a/anphon/dense_symmetric_eigen.h +++ b/anphon/dense_symmetric_eigen.h @@ -14,19 +14,9 @@ namespace PHON_NS { -// Backend seam for dense symmetric eigenproblems (used by SOLVER = DBTE; -// candidate later consumer: batched dynamical-matrix diagonalization). -// -// v1 backend: LAPACK dsyev on the calling rank. Planned backends behind the -// same call: ELPA/ScaLAPACK (memory-distributed; assembly is already -// row-distributed, so only a block-cyclic redistribution is needed) and -// MAGMA/cuSOLVER (single-node GPU). Distributed backends will take this -// call collectively. -// -// A is n x n column-major and is overwritten by the eigenvectors (column j -// = eigenvector of w[j]); eigenvalues are returned in ascending order. -// num_lowest >= 0 asks for only the lowest eigenpairs - the v1 backend -// computes the full spectrum and lets the caller truncate, but -// range-capable backends (dsyevr, ELPA partial, ...) may exploit it. +// Dense symmetric eigensolver using LAPACK dsyev on the calling rank. +// Overwrite column-major A[n][n] with eigenvectors in columns; return +// ascending eigenvalues in w. num_lowest requests a subset, but this +// backend computes the full spectrum for the caller to truncate. void solve_dense_symmetric(int n, std::vector &A, std::vector &w, int num_lowest = -1); } // namespace PHON_NS diff --git a/anphon/dielec.cpp b/anphon/dielec.cpp index 036248da..6585e6c6 100644 --- a/anphon/dielec.cpp +++ b/anphon/dielec.cpp @@ -495,12 +495,9 @@ void Dielec::compute_mode_effective_charge(std::vector> &zst void Dielec::compute_mode_effective_charge(std::vector>> &zstar_mode, const std::complex *const *evec_in) const { - // Compute the mode effective charges defined by Eq. (53) or its numerator of - // Gonze & Lee, PRB 55, 10355 (1997), from mass-weighted Gamma-point - // eigenvectors supplied by the caller. The mass division that converts to - // normal coordinates happens during accumulation; evec_in is not modified. - // The full complex amplitudes are kept so that downstream quantities can be - // made invariant under eigenvector phase choices. + // Compute mode effective charges from mass-weighted Gamma eigenvectors + // [Gonze & Lee, PRB 55, 10355 (1997), Eq. (53)]. Apply mass division during + // accumulation without modifying evec_in; retain complex phases. const auto ns = dynamical->neval; const auto &zstar_atom = borncharge; diff --git a/anphon/diis.h b/anphon/diis.h index 15e0f6a3..b0d58982 100644 --- a/anphon/diis.h +++ b/anphon/diis.h @@ -10,83 +10,44 @@ #include #include -/** - * @brief Generalized Direct Inversion in the Iterative Subspace (GDIIS) optimizer. - * - * DIIS accelerates the convergence of iterative methods by extrapolating - * from a linear combination of previous trial vectors to minimize the error. - * - */ +/// Generalized DIIS optimizer: extrapolate trial vectors to minimize the residual. class GDIIS { public: - /** - * @brief Constructor for GDIIS optimizer. - * - * @param[in] max_history Maximum number of vectors to store in history - * @param[in] mixing_beta Mixing parameter (0 < beta <= 1) - * @param[in] verbosity Verbosity level for logging (0: silent, >0: print info) - */ + /// Store up to max_history vectors with mixing_beta in (0, 1]. + /// verbosity > 0 enables logging. GDIIS(int max_history = 10, double mixing_beta = 0.5, int verbosity = 0); ~GDIIS() = default; - /** - * @brief Updates the DIIS history with a new trial vector and error vector. - * - * @param[in] x_trial Current trial vector - * @param[in] error Error/residual vector at current iteration, e.g., error = f(x) - x - */ + /// Add a trial vector and its residual, e.g. f(x) - x, to the history. void push(const Eigen::VectorXd &x_trial, const Eigen::VectorXd &error); - /** - * @brief Computes the extrapolated vector using DIIS. - * - * @param[out] x_new Extrapolated vector - * @return true if successful, false otherwise - */ + /// Write the DIIS extrapolation to x_new; return true on success. bool extrapolate(Eigen::VectorXd &x_new); - /** - * @brief Clears the DIIS history. - */ + /// Clear the DIIS history. void clear(); - /** - * @brief Returns the current history size. - * - * @return Number of vectors currently stored - */ + /// Return the number of stored vectors. [[nodiscard]] int size() const { return static_cast(history_x.size()); } - /** - * @brief Checks if DIIS has enough history for extrapolation. - * - * @return true if at least 2 vectors are stored - */ + /// Return true when at least two vectors are stored. [[nodiscard]] bool is_ready() const { return size() >= 2; } - /** - * @brief Sets the maximum history size. - * - * @param[in] max_hist New maximum history size - */ + /// Set the maximum history size. void set_max_history(int max_hist) { max_history_ = max_hist; } - /** - * @brief Sets the mixing parameter. - * - * @param[in] beta Mixing parameter (0 < beta <= 1) - */ + /// Set the mixing parameter beta in (0, 1]. void set_mixing_beta(double beta) { mixing_beta_ = beta; @@ -99,11 +60,6 @@ class GDIIS std::deque history_x; ///< History of trial vectors std::deque history_error; ///< History of error vectors - /** - * @brief Solves the DIIS linear system to find optimal coefficients. - * - * @param[out] coeffs Optimal linear combination coefficients - * @return true if successful, false otherwise - */ + /// Solve for DIIS coefficients; return true on success. bool solve_diis_equations(Eigen::VectorXd &coeffs); }; diff --git a/anphon/dynamical.cpp b/anphon/dynamical.cpp index a3986309..fac3ad3f 100644 --- a/anphon/dynamical.cpp +++ b/anphon/dynamical.cpp @@ -1266,19 +1266,11 @@ double Dynamical::freq(const double x) const std::vector Dynamical::detect_acoustic_modes_at_gamma(const std::complex *const *evec_gamma, const double projection_threshold, const bool verbose) const { - // Identify the three acoustic (translational) modes at the Gamma point from the - // eigenvectors instead of the frequencies. A mode is acoustic if and only if it lies - // in the subspace spanned by the three rigid translations, whose mass-weighted, - // orthonormal basis vectors are t_alpha[3*j + beta] = delta_{alpha beta} sqrt(m_j / M). - // The projection P(is) = sum_alpha ||^2 summed over the whole subspace is - // invariant under any rotation or mixing of the degenerate translational modes, so the - // detection works for arbitrary crystal systems and arbitrary orientations of the - // degenerate eigenvectors. The three modes with the largest P are returned; a frequency - // threshold is never consulted, so a soft optical mode collapsing to zero frequency can - // no longer be misclassified as acoustic. - // - // evec_gamma[is][3*j + alpha] : component (j, alpha) of the mass-weighted eigenvector of - // mode is at Gamma. + // Select the three Gamma modes with greatest overlap with rigid translations: + // t_alpha[3*j + beta] = delta_{alpha beta} sqrt(m_j / M), + // P(is) = sum_alpha ||^2. + // This avoids classifying soft optical modes by frequency alone. + // evec_gamma[is][3*j + alpha] holds mass-weighted mode components. const auto ns = neval; const auto natmin = system->get_primcell().number_of_atoms; @@ -1727,14 +1719,10 @@ void Dynamical::exec_interpolation(const unsigned int kmesh_orig[3], std::comple const auto nk = static_cast(nk_dense); int nfail = 0; - // One k-point per thread with its own scratch and one single-threaded LAPACK - // call (MKL and the OpenMP build of OpenBLAS run one thread per call inside a - // parallel region; pthreads OpenBLAS / Accelerate must be pinned, see - // v4_index_transform.h). r2q and the Ewald routines have their own parallel - // regions, which run serially when nested. For a single k-point the region is - // inactive and LAPACK keeps its threads. A failed diagonalization is reported - // after the region (exit() aborts through MPI); the remaining fatal paths - // inside (allocation failure, the Ewald geometry check) abort from the worker. + // Parallelize k points with thread-local scratch and serial nested regions. + // Pin pthreads OpenBLAS / Accelerate threads (see v4_index_transform.h). + // A single k point keeps LAPACK threading. Report eigensolver failures after + // the region; allocation and Ewald geometry failures abort from the worker. #pragma omp parallel for schedule(dynamic) reduction(+ : nfail) if (nk > 1) for (int ik = 0; ik < nk; ++ik) { NDArray, 2> mat_tmp(ns, ns); diff --git a/anphon/elastic_tensor.h b/anphon/elastic_tensor.h index e83dee0f..484be696 100644 --- a/anphon/elastic_tensor.h +++ b/anphon/elastic_tensor.h @@ -44,13 +44,9 @@ struct Tensor6 } }; -// Elastic-constant utilities: readers of the user-provided elastic constants -// used by the SCPH/QHA structural relaxation, and the clamped-ion (Born -// long-wave) stress-energy and elastic tensors computed from the harmonic -// IFCs. All dependencies are explicit constructor arguments (no Pointers -// base). The token parsing of the input files lives in strain_file_parsers.h -// so that it can be unit-tested without a System; the readers here add the -// unit conversion, which needs the volume of the current primitive cell. +// Read SCPH/QHA elastic constants and compute clamped-ion stress and +// elastic tensors from harmonic IFCs. strain_file_parsers.h handles +// parsing; this class converts units using the primitive-cell volume. class ElasticTensor { public: @@ -163,18 +159,10 @@ class ElasticTensor void calc_longwave_brackets3(const std::vector &fcs_cubic, const Eigen::MatrixXd &X, Tensor6 &A_hat) const; - // Third-order (finite-strain) elastic tensor C3_{ij kl mn} in GPa via - // Wallace's Eq. (8.14), combining the restricted brackets with the - // second-order elastic tensor of the same (clamped or relaxed) path. - // symmetrize selects the final projection onto the exact elastic index - // symmetries (minor symmetry within each pair and permutations of the - // three pairs; 48 operations). With rotationally invariant IFCs the - // projection is a no-op; for fitted IFCs (which generally violate the - // cubic rotational invariance) it returns the nearest (least-squares) - // tensor with the exact symmetries. Note that the violated invariance - // relations are NOT index permutations of A_hat itself (A_hat comes out - // exactly symmetric in its own index space), so the projection can only - // be applied at the C3 level. + // Compute C3_{ij kl mn} in GPa using Wallace Eq. (8.14) and C2 from the + // same clamped/relaxed path. symmetrize projects C3 onto minor and pair + // permutation symmetries (48 operations). Apply this at C3, not A_hat: + // rotational-invariance violations are not index permutations of A_hat. void calc_elastic_tensor3(const std::vector &fcs_harmonic, const std::vector &fcs_cubic, bool relax_ions, Tensor6 &C3_gpa, bool symmetrize = true) const; diff --git a/anphon/fcs_phonon.cpp b/anphon/fcs_phonon.cpp index a01d6c0e..4eb7dd6b 100644 --- a/anphon/fcs_phonon.cpp +++ b/anphon/fcs_phonon.cpp @@ -161,13 +161,9 @@ void Fcs_phonon::replicate_force_constants(const int maxorder_in) void Fcs_phonon::replicate_force_constant(const System *system_in, std::vector &fcs_inout) const { - // This function does the following tasks: - // 1. Replicates the force constants originally computed for \Phi_{ij}, \Phi_{ijk}, ..., - // where atom i belongs to the (true) primitive cell to all pairs where atom i belongs to - // the user-defined unit cell. - // 2. Relative vector basis is transformed from the Cartesian to the lattice vector basis - // of the user-defined unit cell. - // 3. Relative vector (relvec) member function is computed from relvec_velocity. + // Replicate IFCs from the true primitive cell to the user-defined cell, + // convert relative vectors to its lattice basis, and derive relvec + // from relvec_velocity. std::vector force_constant_replicate; std::vector relvecs, relvecs_vel; @@ -640,12 +636,9 @@ void Fcs_phonon::append_delta_fc2_from_scph(const std::string &fname_dfc2, std:: } } - // The correction rows are indexed in the cell that the SCPH run used as its - // primitive cell. That cell may be the present primitive cell or an integer - // supercell of it (e.g. an SCPH run with &cell = conventional cell so that - // KMESH_INTERPOLATE matches the DFT supercell): every atom of the SCPH cell is - // folded onto a present primitive atom and translationally equivalent rows - // are added once. + // Fold SCPH correction atoms onto the current primitive cell and count + // translationally equivalent rows once. The SCPH cell may be an integer + // supercell of the current cell. std::vector map_dfc2_to_prim; { Eigen::Matrix3d lavec_dfc2; @@ -717,11 +710,9 @@ void Fcs_phonon::append_delta_fc2_from_scph(const std::string &fname_dfc2, std:: const auto &map_alm = system->get_mapping_super_alm(0); const Eigen::Matrix3d lavec_super_inv = scell.lattice_vector.inverse(); - // Like the rows read from the FCS file, the correction rows start at the reference - // atoms of the FCS file's primitive cell, because replicate_force_constant() spreads - // them with the translations of that cell. Starting at every atom of the present - // primitive cell instead would count the correction natmin(present) / natmin(FCS file) - // times when the present cell is the larger one. + // Start corrections at the FCS-file primitive atoms: replicate_force_constant + // applies that cell's translations. Using all current-cell atoms would + // overcount when the current cell is larger. std::vector> atoms1_s(map_p2s.size()); for (const auto &images: map_alm.from_true_primitive) { atoms1_s[map_s2p[images[0]].atom_num].push_back(images[0]); diff --git a/anphon/four_phonon.cpp b/anphon/four_phonon.cpp index f7508b53..6b52b9ee 100644 --- a/anphon/four_phonon.cpp +++ b/anphon/four_phonon.cpp @@ -137,12 +137,10 @@ inline std::complex phase_factor(const int *q, const int *R, const int * return e; } -// Stage 2 of the Fourier sum for a block of NQ quartets sharing k1: -// A_q[b c d] = sum_{dR} chi_k1(b c d, dR) * E_q(dR), dR = R2 - R3, q = 0..NQ-1. -// Each force-constant term is loaded once per block and the phase tables of -// the block stay in cache. Complex products are written out in real -// arithmetic to keep the loop free of the NaN-checking slow paths of -// std::complex operator*. +// Fourier sum for NQ quartets sharing k1: +// A_q[b c d] = sum_dR chi_k1(b c d, dR) * E_q(dR), dR = R2 - R3. +// Load each IFC once per block; expand complex products into real arithmetic +// to avoid std::complex NaN-checking overhead. template void fourier_stage2_block(const int nrow, const int *row_ptr, const long long *row_bcd, const int *q_diff, const std::complex *chi_k1, const std::complex *const *exp_diff, diff --git a/anphon/ifc_derivative.cpp b/anphon/ifc_derivative.cpp index 8dd737b1..70bc1ee4 100644 --- a/anphon/ifc_derivative.cpp +++ b/anphon/ifc_derivative.cpp @@ -185,16 +185,10 @@ void DerivativeIFC::compute_dV1_dumn(MatrixXcdRowMajor &del_v1_del_umn, } } - // Ad hoc remedy for IFCs that do not satisfy the rotational sum rules - // linking successive orders (symmetrization + offset removal): - // (1) retain only the response to symmetric strains by symmetrizing each - // strain index pair, removing the spurious rigid-rotation response; - // (2) the strain-modulated first-order IFCs still violate the total-force - // acoustic sum rule (the reference-IFC ASR is not sufficient), so - // project out the component conjugate to rigid translations: in the - // mode basis with the mass-weighted metric this is the removal of the - // Gamma-point acoustic components (c_kappa = M_kappa / sum M), which - // leaves every optical generalized force unchanged. + // Correct rotational-sum-rule violations by symmetrizing strain index pairs + // and removing Gamma acoustic components with the mass-weighted metric + // (c_kappa = M_kappa / sum M). This removes rigid-rotation and net-force + // responses while preserving optical generalized forces. for (unsigned int is = 0; is < static_cast(ns); ++is) { for (auto mu = 0; mu < 3; ++mu) { for (auto nu = mu + 1; nu < 3; ++nu) { @@ -265,16 +259,10 @@ void DerivativeIFC::compute_d2V1_dumn2(MatrixXcdRowMajor &del2_v1_del_umn2, } } - // Ad hoc remedy for IFCs that do not satisfy the rotational sum rules - // linking successive orders (symmetrization + offset removal): - // (1) retain only the response to symmetric strains by symmetrizing each - // strain index pair, removing the spurious rigid-rotation response; - // (2) the strain-modulated first-order IFCs still violate the total-force - // acoustic sum rule (the reference-IFC ASR is not sufficient), so - // project out the component conjugate to rigid translations: in the - // mode basis with the mass-weighted metric this is the removal of the - // Gamma-point acoustic components (c_kappa = M_kappa / sum M), which - // leaves every optical generalized force unchanged. + // Correct rotational-sum-rule violations by symmetrizing strain index pairs + // and removing Gamma acoustic components with the mass-weighted metric + // (c_kappa = M_kappa / sum M). This removes rigid-rotation and net-force + // responses while preserving optical generalized forces. for (unsigned int is = 0; is < static_cast(ns); ++is) { std::array, 81> sym_tmp{}; for (auto a1 = 0; a1 < 3; ++a1) @@ -349,16 +337,10 @@ void DerivativeIFC::compute_d3V1_dumn3(MatrixXcdRowMajor &del3_v1_del_umn3, } } - // Ad hoc remedy for IFCs that do not satisfy the rotational sum rules - // linking successive orders (symmetrization + offset removal): - // (1) retain only the response to symmetric strains by symmetrizing each - // strain index pair, removing the spurious rigid-rotation response; - // (2) the strain-modulated first-order IFCs still violate the total-force - // acoustic sum rule (the reference-IFC ASR is not sufficient), so - // project out the component conjugate to rigid translations: in the - // mode basis with the mass-weighted metric this is the removal of the - // Gamma-point acoustic components (c_kappa = M_kappa / sum M), which - // leaves every optical generalized force unchanged. + // Correct rotational-sum-rule violations by symmetrizing strain index pairs + // and removing Gamma acoustic components with the mass-weighted metric + // (c_kappa = M_kappa / sum M). This removes rigid-rotation and net-force + // responses while preserving optical generalized forces. for (unsigned int is = 0; is < static_cast(ns); ++is) { std::array, 729> sym_tmp{}; for (auto a1 = 0; a1 < 3; ++a1) @@ -657,12 +639,9 @@ void DerivativeIFC::compute_dV_dumn_all_real_space(const std::vector &groups, const std::size_t m, const Eigen::Matrix3d &convmat) { - // Computes the m-th derivative of the force constants with respect to strain in real space, - // for all 9^m strain-tensor components in a single scan over fcs_aligned. Each FC entry has - // its tail Cartesian indices mu_j fixed, so it contributes fcs_val * prod_j vec_j[nu_j] to - // the 3^m components sharing its mu-combination (one per nu-combination). - // The input array fcs_aligned is assumed to be sorted by the first (n-m) indices of the - // force constant pairs, where n is the order of the force constants. + // Compute all 9^m strain derivatives in one pass. Each IFC contributes + // fcs_val * prod_j vec_j[nu_j] to the 3^m components with its fixed mu indices. + // fcs_aligned must be sorted by its first n-m indices (n is the IFC order). groups.clear(); if (fcs_aligned.empty()) return; @@ -927,11 +906,9 @@ void DerivativeIFC::compute_dV_dstrain_real_space(const std::vector &strain_dirs, const Eigen::Matrix3d &convmat, const double emit_threshold) { - // This is a helper function that computes the directional derivative of the force constants - // with respect to strain in real space, along the strain tensor strain_dirs[j] for the j-th derivative. - // The derivative order is determined by the size of the strain_dirs vector. - // The input array fcs_aligned is assumed to be sorted by the first (n-m) indices of the force constant pairs, - // where m is the derivative order and n is the order of the force constants (n=2: harmonic, n=3: cubic, etc.). + // Compute real-space directional strain derivatives along strain_dirs[j]. + // With m = strain_dirs.size() and n the IFC order, fcs_aligned must be + // sorted by its first n-m indices. if (fcs_aligned.empty()) { delta_fcs.clear(); return; @@ -1803,14 +1780,9 @@ void DerivativeIFC::process_strain_harmonic_set(const std::vector(symmetry_.SymmListWithMap_ref.size()); const auto &map_p2s = system_.get_map_p2s(0); const auto &map_s2p = system_.get_map_s2p(0); diff --git a/anphon/ifc_derivative.h b/anphon/ifc_derivative.h index 96f65940..f1c0543f 100644 --- a/anphon/ifc_derivative.h +++ b/anphon/ifc_derivative.h @@ -58,12 +58,9 @@ class DerivativeIFC verbosity_ = verbosity; } - // Single-pass computation of the m-th strain derivative of the IFCs in real - // space for ALL 9^m strain-tensor components at once. One scan over - // fcs_aligned replaces the 9^m per-component scans; use - // extract_strain_component to materialize one component's delta IFCs. - // fcs_aligned must be sorted by the first (n-m) indices - // (sort_by_heading_indices(m)). + // Compute all 9^m real-space strain derivatives in one pass. Use + // extract_strain_component for one component. fcs_aligned must be sorted + // by its first n-m indices (sort_by_heading_indices(m)). static void compute_dV_dumn_all_real_space(const std::vector &fcs_aligned, std::vector &groups, std::size_t m, const Eigen::Matrix3d &convmat); @@ -103,21 +100,12 @@ class DerivativeIFC const Eigen::VectorXd &sublattice_displacement, double emit_threshold); - // Displacement field of a homogeneous deformation u plus sublattice - // displacements S: atom (l kappa) moves by - // d_lambda = sum_nu u(lambda, nu) R_nu + S(3*kappa + lambda), - // with R the same relative vector used by the strain kernels (valid by - // the acoustic sum rule). S may be empty (purely affine deformation). - // - // Conventions: u is the dimensionless displacement-gradient tensor - // dX_mu/dx_nu - delta_{mu nu} (a homogeneous deformation of all space, - // independent of any cell choice); S is in Cartesian bohr and is indexed - // by the atoms of the USER-defined primitive cell of the run (the &cell - // input), i.e. the same numbering as pairs[].index/3 after - // replicate_force_constant. A non-primitive &cell is fully supported: - // S is then periodic with that larger cell, which is exactly what is - // needed for SCPH/QHA distortion patterns that break the true primitive - // periodicity. + // Homogeneous deformation plus sublattice displacements: + // d_lambda = sum_nu u(lambda, nu) R_nu + S(3*kappa + lambda). + // R is the strain-kernel relative vector (valid by ASR); u is the + // dimensionless displacement gradient. S is in Cartesian bohr, indexed + // by user-cell atoms as in pairs[].index/3 after replicate_force_constant. + // S may be empty for affine deformation or periodic with a non-primitive cell. struct DeformationField { // Dimensionless displacement-gradient tensor (the u_tensor of the diff --git a/anphon/integration.cpp b/anphon/integration.cpp index dd088a80..d532664f 100644 --- a/anphon/integration.cpp +++ b/anphon/integration.cpp @@ -423,6 +423,9 @@ void Integration::insertion_sort(double *a, int *ind, int n) void AdaptiveSmearingSigma::setup(const PhononVelocity *phvel_class, const KpointMeshUniform *kmesh_in, const Eigen::Matrix3d &lavec_p_in, const Eigen::Matrix3d &rlavec_p_in) { + // Use finite-difference velocities for adaptive widths. Degenerate-mode + // velocities are basis dependent, and widths have a fixed 2e-5 Ry floor, + // so cell dependence can persist. Use ISMEAR = 1 or 0 for cell independence. phvel_class->get_phonon_group_velocity_mesh(*kmesh_in, lavec_p_in, false, vel); for (auto u = 0; u < 3; u++) { diff --git a/anphon/interpolation.cpp b/anphon/interpolation.cpp index d16943dc..6954b004 100644 --- a/anphon/interpolation.cpp +++ b/anphon/interpolation.cpp @@ -196,20 +196,10 @@ void PHON_NS::fourier_dymat_k_to_r(const unsigned int nk1, const unsigned int nk const unsigned int ns, const std::complex *const *const *dymat_k, std::complex ***dymat_r) { - // Forward DFT (k -> r), including the 1/N normalization, of all (is, js) - // components of a coarse-mesh dynamical-matrix array at once: + // Forward DFT of all dynamical-matrix components: // dymat_r[is][js][r] = (1/N) sum_k dymat_k[is][js][k] e^{-2 pi i k.r}. - // - // This used to be done with one FFTW plan per (is, js) pair, re-created at - // every call: fftw_execute() always transforms the arrays its plan was - // created with, and fftw_execute_dft() on different arrays carries strict - // alignment requirements, so the plan could not be hoisted safely. The - // coarse mesh is tiny, so plan creation dominated the transform itself. - // The explicit DFT matrix below reproduces the FFTW_FORWARD convention - // (row-major multi-index, negative exponent) exactly, is alignment-free, - // and applies to all ns^2 components as a single matrix product; the arrays - // from allocate() are contiguous, with (is, js) blocks of length - // nk1*nk2*nk3 each. + // A single matrix product avoids repeated FFTW plan creation on small meshes. + // Arrays are contiguous, with nk1*nk2*nk3 entries per (is, js) block. using namespace Eigen; diff --git a/anphon/iterativebte.cpp b/anphon/iterativebte.cpp index ece097c6..2d548aec 100644 --- a/anphon/iterativebte.cpp +++ b/anphon/iterativebte.cpp @@ -102,8 +102,9 @@ void Iterativebte::setup_iterative() ntemp = conductivity->ntemp; Temperature = conductivity->temperature; - // Full-grid velocities in atomic units on every rank (calc_kappa and - // the boundary rate convert units at the point of use). + // Full-grid velocities in atomic units on every rank; convert at use. + // IBTE results can depend strongly on cell choice; the cause remains + // unresolved, and replacing these velocities did not resolve it. phonon_velocity->gather_group_velocities_mesh(*dos->kmesh_dos.get(), system->get_primcell().lattice_vector, vel, @@ -759,15 +760,10 @@ void Iterativebte::project_wedge_vector(std::vector &v, const std::vecto bool Iterativebte::solve_direct_at_temperature(const int itemp, const double beta, double **sqrt_occ, double **Qfin_loc, int &iterations_out, double &residual_out) { - // SOLVER = DBTE: assemble the multiplicity-symmetrized dense operator - // from the stored L entries (one-application cost), transform it to the - // Omega normalization Omega = D^{-1/2} A D^{-1/2} with D = diag(n(n+1)) - // - whose diagonal is 1/tau, so the eigenvalues are scattering rates - - // take its full eigendecomposition, and report the spectrum diagnostics - // that the matrix-free solvers cannot access: discretization asymmetry, - // positive-semidefiniteness, near-null modes with their overlap onto - // the momentum-drift directions, and the sensitivity of kappa to a - // low-eigenvalue cutoff. Intended for small meshes. + // DBTE diagonalizes Omega = D^{-1/2} A D^{-1/2}, D = diag(n(n+1)), + // using the multiplicity-symmetrized operator. Eigenvalues are scattering + // rates. Report asymmetry, negative/near-null modes, momentum-drift overlaps, + // and kappa sensitivity to the low-eigenvalue cutoff. Intended for small meshes. const auto nk_irred = dos->kmesh_dos->nk_irred; const size_t nrows = static_cast(nk_irred) * ns; const size_t nrows3 = nrows * 3; @@ -891,12 +887,9 @@ bool Iterativebte::solve_direct_at_temperature(const int itemp, const double bet if (mympi->my_rank == 0) { - // Degeneracy-reduced basis: one orthonormal (symmetric) combination - // per degenerate block, u_I = d^{-1/2} sum_{s in I} e_s. This - // performs the degeneracy averaging of the iterative solvers - // exactly, without the 3(d-1) artificial null modes a projector - // P S P would inject into the spectrum. Masked rows never share a - // block with active ones (degenerate partners have equal omega). + // Use one normalized combination u_I = d^{-1/2} sum_{s in I} e_s per + // degenerate block, avoiding the artificial null modes of P S P. + // Keep masked rows separate from active rows. const auto tol_omega = 1.0e-7; std::vector> blocks; // [first,last) active scalar rows within one ik for (unsigned int ik = 0; ik < nk_irred; ++ik) { @@ -965,12 +958,9 @@ bool Iterativebte::solve_direct_at_temperature(const int itemp, const double bet } } - // Restrict to the little-group-invariant subspace. The collision - // operator is only defined on fields with dF(Rk) = R dF(k); the - // complementary vector components at high-symmetry k points carry an - // arbitrary (choice-of-operation dependent, non-symmetric) extension - // that must not enter the spectrum. Per block, an orthonormal basis - // of range(P_littlegroup) is taken from the projector eigenvectors. + // Restrict to little-group-invariant fields dF(Rk) = R dF(k). Use projector + // eigenvectors to span range(P_littlegroup), excluding arbitrary + // complementary components from the spectrum. std::vector Ebasis(blocks.size()); std::vector mdim(blocks.size()), offs(blocks.size()); int ndim = 0; @@ -1244,18 +1234,11 @@ bool Iterativebte::solve_direct_at_temperature(const int itemp, const double bet bool Iterativebte::solve_variational_cg(const int itemp, const double beta, double **sqrt_occ, double **Qfin_loc, const double *x0_wedge, int &iterations_out, double &residual_out) { - // Solve (Q_diag + W) dF = b with preconditioned conjugate gradients. - // In the dF variables the off-diagonal couplings are the equilibrium - // rates n1 n2 (n3+1) |V3|^2 delta(...), which detailed balance makes - // symmetric under exchange of the coupled modes on the energy shell, so - // the operator is self-adjoint under the plain star-multiplicity- - // weighted inner product = sum_{wedge k,s} mult_k u.v used by all - // dot products below; the preconditioner is the diagonal. (Smearing - // spreads slightly off shell, which the self-adjointness check below - // quantifies.) - // Because kappa is the value of the variational functional, its error is - // quadratic in the residual, so ITER_THRESHOLD acts on the relative - // residual here. + // Solve (Q_diag + W) dF = b with diagonal-preconditioned CG and the + // star-weighted inner product = sum_{wedge k,s} mult_k u.v. + // Detailed balance gives on-shell self-adjointness; check discretization + // effects below. ITER_THRESHOLD applies to the relative residual; + // the variational kappa error is quadratic in that residual. const auto nk_irred = dos->kmesh_dos->nk_irred; const size_t nrows = static_cast(nk_irred) * ns; const size_t nrows3 = nrows * 3; diff --git a/anphon/kappa_result_io.cpp b/anphon/kappa_result_io.cpp index 79000e9d..3de939f6 100644 --- a/anphon/kappa_result_io.cpp +++ b/anphon/kappa_result_io.cpp @@ -189,10 +189,11 @@ struct KappaResultIOH5::Impl } } - // Create the group, attributes, static datasets, and preallocated - // gamma/flag (and, in tdep mode, frequency/velocity) datasets of a - // channel. In tdep mode the frequency/velocity content is written later - // (write_basis_slices); only the shapes are taken from cmeta here. + // Preallocate channel datasets; tdep frequency/velocity values are filled + // later by write_basis_slices. Stamp the transport formulation before + // publishing the .part file so interrupted runs can restart correctly. + std::string formulation{"legacy"}; + auto create_channel(HighFive::File &fh, const KappaChannelMetaH5 &cmeta) const -> void { using namespace H5Easy; @@ -259,8 +260,35 @@ struct KappaResultIOH5::Impl } } - // Fill the frequency/velocity slices of this run's temperature columns - // (tdep mode only). Pure raw-data writes into preallocated datasets. + // Write this run's temperature slices and velocity diads, creating the + // diad dataset on demand for older restart files. + auto write_velocity_diad(const KappaChannelMetaH5 &cmeta) const -> void + { + if (cmeta.velocity_diad.empty()) return; + const auto path = channel_path(cmeta.tag); + const auto name = path + "/velocity_diad"; + const size_t nequiv_total = cmeta.velocity_diad.size() / (static_cast(cmeta.ns) * 9); + const auto nt = nt_file(); + + if (!file->exist(name)) { + if (tdep) { + h5_create_dataset_prealloc(*file, name, {nt, nequiv_total, cmeta.ns, 3, 3}) + .createAttribute("unit", std::string("(m/s)^2")); + } else { + file->createDataSet(name, HighFive::DataSpace({nequiv_total, cmeta.ns, 3, 3})) + .createAttribute("unit", std::string("(m/s)^2")); + } + } + auto dset = file->getDataSet(name); + if (tdep) { + for (const auto col: run_cols) { + dset.select({col, 0, 0, 0, 0}, {1, nequiv_total, cmeta.ns, 3, 3}).write_raw(cmeta.velocity_diad.data()); + } + } else { + dset.write_raw(cmeta.velocity_diad.data()); + } + } + auto write_basis_slices(const KappaChannelMetaH5 &cmeta) const -> void { if (!tdep) return; @@ -362,6 +390,7 @@ struct KappaResultIOH5::Impl // Validity marker as a raw-data dataset (not an attribute) so setting // it never touches file metadata. Per temperature in tdep mode. h5_create_dataset_prealloc(fh, "/kappa/valid", {tdep ? nt : 1}); + fh.getGroup("/kappa").createAttribute("formulation", formulation); // Machine-readable summary of the same information on the group. auto gk = fh.getGroup("/kappa"); @@ -806,6 +835,10 @@ void KappaResultIOH5::open_or_create(const KappaFileMetaH5 &fmeta, const KappaCh impl->file_temps = fmeta.temperatures; impl->file_fc2temps.assign(impl->file_temps.size(), fmeta.fc2_temperature); + // Transport formulation recorded in the existing file (files from before the stamp + // used the legacy formulation). + std::string existing_formulation; + auto have_existing_formulation = false; if (std::filesystem::exists(impl->filename)) { const HighFive::File oldfile(impl->filename, HighFive::File::ReadOnly); check_h5_schema(oldfile, h5_schema_kappa_result, kappa_version_tdep); @@ -863,6 +896,13 @@ void KappaResultIOH5::open_or_create(const KappaFileMetaH5 &fmeta, const KappaCh } else { need_rebuild = true; } + if (oldfile.exist("/kappa")) { + have_existing_formulation = true; + existing_formulation = "legacy"; + auto grp = oldfile.getGroup("/kappa"); + if (grp.hasAttribute("formulation")) grp.getAttribute("formulation").read(existing_formulation); + if (existing_formulation == "standard") existing_formulation = "legacy"; + } if (!impl->kappa_group_matches(oldfile)) need_rebuild = true; if (!impl->isotope_group_matches(oldfile)) need_rebuild = true; } else { @@ -877,6 +917,40 @@ void KappaResultIOH5::open_or_create(const KappaFileMetaH5 &fmeta, const KappaCh impl->compute_run_cols(); + // Reject formulation changes only if old kappa rows would survive, as with + // tdep temperatures outside this run. Ordinary restarts recompute all kappa + // from formulation-independent 3ph linewidths and can update the stamp. + { + auto retains_foreign_kappa = false; + if (impl->tdep) { + for (size_t j = 0; j < old_temps.size(); ++j) { + if (j < old_valid.size() && !old_valid[j]) continue; // never assembled: nothing to mix + auto covered = false; + for (const auto t: fmeta.temperatures) { + if (std::abs(t - old_temps[j]) < eps6) { + covered = true; + break; + } + } + if (!covered) { + retains_foreign_kappa = true; + break; + } + } + } + if (have_existing_formulation && retains_foreign_kappa && existing_formulation != transport_formulation) { + exit("KappaResultIOH5", + ("The existing temperature-resolved result file holds kappa for temperatures this run does not " + "recompute, assembled with transport formulation '" + + existing_formulation + "', while this run uses '" + transport_formulation + + "'. Use a different PREFIX, include those temperatures in this run, or set " + "ALAMODE_LEGACY_VELOCITY to match the file.") + .c_str()); + } + } + + impl->formulation = transport_formulation; + if (need_rebuild) { if (impl->tdep) { impl->rebuild_tdep({channel}, drops, old_temps, old_valid); @@ -887,6 +961,7 @@ void KappaResultIOH5::open_or_create(const KappaFileMetaH5 &fmeta, const KappaCh impl->file = std::make_unique(impl->filename, HighFive::File::ReadWrite); } impl->write_basis_slices(channel); + impl->write_velocity_diad(channel); impl->reset_kappa_valid(); } @@ -922,10 +997,13 @@ void KappaResultIOH5::ensure_channel(const KappaChannelMetaH5 &channel, const bo } } impl->write_basis_slices(channel); + impl->write_velocity_diad(channel); } bool KappaResultIOH5::open_or_create_for_ibte(const KappaFileMetaH5 &fmeta, const IbteMetaH5 &imeta, const bool reset) { + impl->formulation = "legacy"; // IBTE uses finite-difference velocities throughout + transport_formulation = "legacy"; impl->fmeta = fmeta; impl->ntemp = fmeta.temperatures.size(); impl->tdep = fmeta.temperature_resolved; @@ -1139,6 +1217,22 @@ void KappaResultIOH5::store_kappa(const double *const *const *kappa_peierls, con const double *const *const *kappa_coherent, const double *const *const *kappa_spec, const double *const *gamma_isotope) { + // Stamp how the transport weights were built. This belongs on the kappa result, not on + // the velocities dataset (which always holds finite-difference data), and it is written + // here because store_kappa is reached on every path, including restart and + // temperature-resolved output, so the stamp cannot go stale. + if (impl->file && !transport_formulation.empty()) { + auto &fh = *impl->file; + if (fh.exist("/kappa")) { + auto grp = fh.getGroup("/kappa"); + if (grp.hasAttribute("formulation")) { + grp.getAttribute("formulation").write(transport_formulation); + } else { + grp.createAttribute("formulation", transport_formulation); + } + } + } + impl->write_tensor_txx("/kappa/kappa_peierls", kappa_peierls); if (impl->fmeta.isotope > 0 && gamma_isotope) { diff --git a/anphon/kappa_result_io.h b/anphon/kappa_result_io.h index 7dfb871c..b555b185 100644 --- a/anphon/kappa_result_io.h +++ b/anphon/kappa_result_io.h @@ -87,6 +87,11 @@ struct KappaChannelMetaH5 std::vector> equiv_knum; // full-mesh k indices of each irreducible star Eigen::MatrixXd frequencies; // [nk_irred, ns], cm^-1 std::vector velocities; // [sum(multiplicity)*ns*3] flattened, m/s + // Velocity diad W[k][s][a][b] = sum_{s' in D(s)} Re(V^a_{ss'} V^b_{s's}), + // flattened [sum(multiplicity)*ns*9], in (m/s)^2. Block sums give + // Tr(P V^a P V^b P); cumulative kappa must use W at degeneracies. + // Empty for legacy transport and channels without a velocity matrix (4ph). + std::vector velocity_diad; }; // Crash-safe writer/reader of PREFIX.kappa.h5. All methods must be called @@ -165,6 +170,9 @@ class KappaResultIOH5 // gamma_isotope is the per-mode isotope linewidth [nk_irred][ns] on the // 3ph mesh (internal units), stored when the file was set up with // isotope scattering enabled. + // How the transport weights were built; stamped onto /kappa by store_kappa. + std::string transport_formulation{"standard"}; + void store_kappa(const double *const *const *kappa_peierls, const double *const *const *kappa_3ph_only, const double *const *const *kappa_coherent, const double *const *const *kappa_spec, const double *const *gamma_isotope); diff --git a/anphon/mode_symmetry.cpp b/anphon/mode_symmetry.cpp index 55fe42dd..56b9f6eb 100644 --- a/anphon/mode_symmetry.cpp +++ b/anphon/mode_symmetry.cpp @@ -212,13 +212,8 @@ void ModeSymmetry::analyze_irreps_at_gamma() } } - // ------------------------------------------------------------------ - // Eigenvectors of the analytic dynamical matrix at Gamma. - // With kvec = 0 the directional nonanalytic term vanishes for - // NONANALYTIC = 1/2; NONANALYTIC = 3 must go through the Ewald path - // with the dipole-free force constants (plain eval_k would add an - // uninitialized nonanalytic matrix in that mode). - // ------------------------------------------------------------------ + // Compute analytic Gamma eigenvectors. kvec = 0 removes the directional + // term for NONANALYTIC = 1/2; mode 3 requires Ewald with dipole-free IFCs. std::vector eval_raw(ns), omega(ns); NDArray, 2> evec; evec.resize(ns, ns); @@ -434,12 +429,8 @@ void ModeSymmetry::analyze_irreps_at_gamma() Eigen::Vector3d zhat = Eigen::Vector3d::Zero(); Eigen::Vector3d xhat = Eigen::Vector3d::Zero(); if (symmetry->has_spg_dataset) { - // Authoritative conventional frame from the spglib dataset: - // L_conv = L_input * P^-1 with P the dataset transformation - // matrix (basis change only, expressed directly in the Cartesian - // frame of the calculation). Sanity-checked against the - // operations: when a unique principal axis (order >= 3) exists, - // the conventional c axis must be parallel to it. + // Conventional frame: L_conv = L_input * P^-1 using spglib's transformation. + // When a unique principal axis of order >= 3 exists, check that c is parallel to it. const Eigen::Matrix3d &lavec = system->get_primcell().lattice_vector; const Eigen::Matrix3d &pmat = symmetry->spg_transformation_matrix; if (std::abs(pmat.determinant()) > 1.0e-8) { @@ -554,14 +545,10 @@ void ModeSymmetry::analyze_irreps_at_gamma() } } - // ------------------------------------------------------------------ - // Representation-level validation of the class->column assignment (the - // hybrid scheme of the plan): among all bijections consistent with the - // (kind, nelem) buckets, keep those whose decompositions of Gamma_V, - // Gamma_total, and every symmetry-clean multiplet are non-negative - // integers with consistent dimensions and Gamma_optic >= 0; the - // geometrically chosen assignment must be a survivor. - // ------------------------------------------------------------------ + // Validate class-column bijections within (kind, nelem) buckets using + // non-negative integer decompositions of Gamma_V, Gamma_total, and clean + // multiplets, consistent dimensions, and Gamma_optic >= 0. + // The geometric assignment must pass. if (labels_ok && pg) { auto assignment_ok = [&](const std::vector &a2c) { const auto n_vec = pointgroup::decompose_representation(*pg, a2c, nelem_of_class, trace_class); diff --git a/anphon/optimizers.cpp b/anphon/optimizers.cpp index d9f80d93..4044a5c7 100644 --- a/anphon/optimizers.cpp +++ b/anphon/optimizers.cpp @@ -270,12 +270,9 @@ double FarkasIII_Optimizer::angle_threshold(const int n_vectors) Eigen::MatrixXd FarkasIII_Optimizer::effective_inverse_Hessian() const { - // H stores the inverse Hessian. Build the RFO level-shifted effective inverse Hessian, i.e. - // (H_direct + lambda*I)^{-1}, with lambda >= 0 chosen so the smallest direct-Hessian - // curvature becomes at least curvature_floor. This bounds the step along near-zero-curvature - // (soft) modes and turns negative-curvature directions into controlled descent steps - // (Farkas-Schlegel eqn 7, H_eff = H_direct + lambda*I). H need NOT be positive definite here - // (a strongly soft initial Hessian can be indefinite). + // Build (H_direct + lambda*I)^{-1} from the inverse Hessian H. Choose + // lambda >= 0 to raise the smallest curvature to curvature_floor, limiting + // soft-mode steps and allowing indefinite H (Farkas-Schlegel Eq. 7). Eigen::SelfAdjointEigenSolver es(H); const Eigen::VectorXd &nu = es.eigenvalues(); // eigenvalues of the inverse Hessian const Eigen::MatrixXd &U = es.eigenvectors(); diff --git a/anphon/optimizers.h b/anphon/optimizers.h index 4c747ec0..23ed9315 100644 --- a/anphon/optimizers.h +++ b/anphon/optimizers.h @@ -101,11 +101,9 @@ class FarkasIII_Optimizer: public Optimizer // angle acceptance criterion. Eigen::VectorXd update(const Eigen::VectorXd &point, const Eigen::VectorXd &gradient); - // Farkas-Schlegel "controlled GDIIS" (Phys. Chem. Chem. Phys. 2002, 4, 11): recomputes the - // error vectors with a per-step RFO level-shifted Hessian, grows the DIIS subspace from the - // most recent point keeping the last acceptable step, applies the four acceptance criteria - // (a) angle, (b) step length, (c) coefficient sum, (d) near-singularity, and permanently - // discards points when the angle exceeds 90 degrees. Enabled by default; disabled by GDIIS_PLAIN = 1. + // Controlled GDIIS (Farkas-Schlegel, PCCP 2002, 4, 11): use RFO-shifted + // errors and check angle, step length, coefficient sum, and singularity. + // Discard points with angles > 90 degrees. GDIIS_PLAIN = 1 disables it. Eigen::VectorXd update_controlled(const Eigen::VectorXd &point, const Eigen::VectorXd &gradient); void set_inverse_Hessian(const int dim, const std::vector> &hessian); diff --git a/anphon/phonon_velocity.cpp b/anphon/phonon_velocity.cpp index 66d67813..b83f492b 100644 --- a/anphon/phonon_velocity.cpp +++ b/anphon/phonon_velocity.cpp @@ -9,10 +9,14 @@ or http://opensource.org/licenses/mit-license.php for information. */ #include "phonon_velocity.h" +#include #include +#include #include +#include #include "cell_shift_table.h" #include "constants.h" +#include "degeneracy_utils.h" #include "dense_hermitian_eigen.h" #include "dynamical.h" #include "error.h" @@ -48,11 +52,72 @@ void PhononVelocity::set_default_variables() void PhononVelocity::deallocate_variables() {} +// Default transport uses the unsymmetrized velocity matrix with nonanalytic +// connection, block-trace Peierls weights, cross-block coherent pairs, +// block boundary speeds, and matrix-diagonal PRINTVEL. +// ALAMODE_LEGACY_VELOCITY=1 restores finite-difference velocities and +// elementwise symmetrization without the nonanalytic velocity term. +bool PhononVelocity::legacy_velocity() +{ + return std::getenv("ALAMODE_LEGACY_VELOCITY") != nullptr; +} + void PhononVelocity::setup_velocity() { MPI_Bcast(&print_velocity, 1, MPI_CXX_BOOL, 0, MPI_COMM_WORLD); } +// Project the velocity-matrix diagonal onto the band-path direction, +// avoiding finite differences across sorted-band crossings. Diagonals +// within degenerate multiplets remain basis dependent; only block traces +// are invariant. +void PhononVelocity::get_phonon_group_velocity_bandstructure_velmat(const KpointBandStructure *kpoint_bs_in, + const Eigen::Matrix3d &lavec_p, + const std::vector &fc2_in, + double **phvel_out) const +{ + const auto nk = kpoint_bs_in->nk; + const auto ns = dynamical->neval; + + NDArray, 3> velmat_k; + NDArray, 2> evec_k; + NDArray eval_k; + velmat_k.resize(ns, ns, 3); + evec_k.resize(ns, ns); + eval_k.resize(ns); + + const auto &fc2_vel = (dynamical->nonanalytic == 3) ? ewald->fc2_without_dipole : fc2_in; + + for (auto ik = 0u; ik < nk; ++ik) { + if (dynamical->nonanalytic == 3) { + dynamical->eval_k_ewald(kpoint_bs_in->xk[ik], + kpoint_bs_in->kvec_na[ik], + ewald->fc2_without_dipole, + eval_k, + evec_k, + true); + } else { + dynamical->eval_k(kpoint_bs_in->xk[ik], kpoint_bs_in->kvec_na[ik], fc2_in, eval_k, evec_k, true); + } + for (auto is = 0u; is < ns; ++is) eval_k[is] = dynamical->freq(eval_k[is]); + + velocity_matrix_analytic(kpoint_bs_in->xk[ik], fc2_vel, eval_k, evec_k, velmat_k); + add_nonanalytic_velocity_matrix(kpoint_bs_in->xk[ik], eval_k, evec_k, velmat_k, kpoint_bs_in->kvec_na[ik]); + + for (auto is = 0u; is < ns; ++is) { + double v[3]; + for (auto j = 0; j < 3; ++j) v[j] = velmat_k[is][is][j].real(); + rotvec(v, v, lavec_p); + auto vproj = 0.0; + for (auto j = 0; j < 3; ++j) vproj += (v[j] / (2.0 * pi)) * kpoint_bs_in->kvec_na[ik][j]; + phvel_out[ik][is] = vproj; + } + } + velmat_k.clear(); + evec_k.clear(); + eval_k.clear(); +} + void PhononVelocity::get_phonon_group_velocity_bandstructure(const KpointBandStructure *kpoint_bs_in, const Eigen::Matrix3d &lavec_p, const Eigen::Matrix3d &rlavec_p, @@ -260,13 +325,10 @@ void PhononVelocity::gather_group_velocities_mesh(const KpointMeshUniform &kmesh NDArray &vel_out, const double unit_factor, const bool bcast_full) const { - // Allocate-and-fill wrapper around get_phonon_group_velocity_mesh_mpi - // shared by the RTA and IBTE setup paths. The gathered velocities live - // on rank 0 only unless bcast_full is set; other ranks get a dummy - // allocation. unit_factor is applied before the broadcast, so every - // caller states its unit convention in one place (1.0 keeps atomic - // units; Bohr_in_Angstrom * 1.0e-10 / time_ry converts to m/s). The - // caller owns and deallocates vel_out. + // Allocate and gather velocities on rank 0, or all ranks if bcast_full. + // Apply unit_factor before broadcasting: 1.0 keeps atomic units; + // Bohr_in_Angstrom * 1.0e-10 / time_ry gives m/s. + // Other ranks receive dummy storage; the caller deallocates vel_out. const auto nk = kmesh_in.nk; const auto neval = dynamical->neval; @@ -293,86 +355,206 @@ void PhononVelocity::gather_group_velocities_mesh(const KpointMeshUniform &kmesh } } -void PhononVelocity::calc_phonon_velmat_mesh(std::complex ****velmat_out) const +// Fill PRINTVEL from the analytic velocity-matrix diagonal, without +// elementwise symmetrization. Units match get_phonon_group_velocity_mesh +// (Cartesian, divided by 2 pi, no SI factor). Compute the full matrix +// but retain only its diagonal. Adaptive smearing keeps finite differences. +void PhononVelocity::get_phonon_group_velocity_mesh_velmat(const KpointMeshUniform &kmesh_in, + const Eigen::Matrix3d &lavec_p, double ***phvel3_out) const { - const auto nk = dos->kmesh_dos->nk; + const auto nk = kmesh_in.nk; const auto ns = dynamical->neval; - NDArray, 4> velmat_loc; - NDArray displs; - NDArray sendcount; - NDArray recvcount; - std::vector nk_proc; - std::vector ik_begin_proc, ik_end_proc; - - const auto factor = Bohr_in_Angstrom * 1.0e-10 / (time_ry * 2.0 * pi); - - if (mympi->my_rank == 0 && writes->getVerbosity() > 0) { - std::cout << " Calculating group velocity matrix of phonons on uniform grid ... "; - } + // Self-contained (diagonalizes for itself) so it does not depend on dos->dymat_dos + // having been filled or on the mesh being kmesh_dos. + NDArray, 3> velmat_k; + NDArray, 2> evec_k; + NDArray eval_k; + velmat_k.resize(ns, ns, 3); + evec_k.resize(ns, ns); + eval_k.resize(ns); + + const auto &fc2_vel = + (dynamical->nonanalytic == 3) ? ewald->fc2_without_dipole : fcs_phonon->force_constant_with_cell[0]; + + for (auto ik = 0u; ik < nk; ++ik) { + double kvec[3]; + for (auto j = 0; j < 3; ++j) kvec[j] = kmesh_in.xk[ik][j]; + rotvec(kvec, kvec, system->get_primcell().reciprocal_lattice_vector, 'T'); + const auto norm = std::sqrt(kvec[0] * kvec[0] + kvec[1] * kvec[1] + kvec[2] * kvec[2]); + if (norm > eps) { + for (auto j = 0; j < 3; ++j) kvec[j] /= norm; + } - sendcount.resize(mympi->nprocs); - recvcount.resize(mympi->nprocs); - nk_proc.resize(mympi->nprocs); + if (dynamical->nonanalytic == 3) { + dynamical->eval_k_ewald(kmesh_in.xk[ik], kvec, ewald->fc2_without_dipole, eval_k, evec_k, true); + } else { + dynamical->eval_k(kmesh_in.xk[ik], kvec, fcs_phonon->force_constant_with_cell[0], eval_k, evec_k, true); + } + for (auto is = 0u; is < ns; ++is) eval_k[is] = dynamical->freq(eval_k[is]); - auto nk_loc = nk / mympi->nprocs; - auto nk_res = nk - nk_loc * mympi->nprocs; + velocity_matrix_analytic(kmesh_in.xk[ik], fc2_vel, eval_k, evec_k, velmat_k); + add_nonanalytic_velocity_matrix(kmesh_in.xk[ik], eval_k, evec_k, velmat_k); - for (auto i = 0; i < mympi->nprocs; ++i) { - nk_proc[i] = nk_loc; - if (i < nk_res) ++nk_proc[i]; - sendcount[i] = 3 * ns * ns * nk_proc[i]; - recvcount[i] = sendcount[i]; + for (auto is = 0u; is < ns; ++is) { + double v[3]; + for (auto j = 0; j < 3; ++j) v[j] = velmat_k[is][is][j].real(); + rotvec(v, v, lavec_p); + for (auto j = 0; j < 3; ++j) phvel3_out[ik][is][j] = v[j] / (2.0 * pi); + } } + velmat_k.clear(); + evec_k.clear(); + eval_k.clear(); +} - if (mympi->my_rank == 0) { - displs.resize(mympi->nprocs); - displs[0] = 0; - for (auto i = 1; i < mympi->nprocs; ++i) { - displs[i] = displs[i - 1] + recvcount[i - 1]; +// Gather k-distributed contiguous records (stride elements per k) to rank 0 in chunks +// whose element counts fit an int. A single MPI_Gatherv with count nk*ns*ns*3 overflows +// 32-bit counts for large systems (e.g. 20^3 mesh x 300 branches = 2.16e9 elements). +template +static void gather_k_records(const T *local, const std::vector &nk_proc, const int my_rank, const int nprocs, + const size_t stride, const MPI_Datatype type, T *out) +{ + // k offsets and chunk endpoints in size_t: only the per-chunk MPI counts and + // displacements are ever narrowed to int, and those are bounded below INT_MAX/2. + std::vector kbeg(nprocs + 1, 0); + for (auto r = 0; r < nprocs; ++r) kbeg[r + 1] = kbeg[r] + static_cast(nk_proc[r]); + const auto nk = kbeg[nprocs]; + + const auto max_elems = static_cast(std::numeric_limits::max()) / 2; + if (stride > max_elems) exit("gather_k_records", "A single k-point record exceeds the MPI count limit."); + const auto chunk_k = std::max(1, max_elems / stride); + + std::vector cnt(nprocs), dsp(nprocs); + for (size_t k0 = 0; k0 < nk; k0 += chunk_k) { + const auto k1 = std::min(nk, k0 + chunk_k); + for (auto r = 0; r < nprocs; ++r) { + const auto lo = std::max(k0, kbeg[r]); + const auto hi = std::min(k1, kbeg[r + 1]); + const auto n = hi > lo ? hi - lo : 0; + cnt[r] = static_cast(n * stride); + dsp[r] = n > 0 ? static_cast((lo - k0) * stride) : 0; } + const auto lo_me = std::max(k0, kbeg[my_rank]); + const auto off_me = (lo_me > kbeg[my_rank] ? lo_me - kbeg[my_rank] : 0) * stride; + MPI_Gatherv(cnt[my_rank] > 0 ? local + off_me : nullptr, + cnt[my_rank], + type, + my_rank == 0 ? out + k0 * stride : nullptr, + my_rank == 0 ? cnt.data() : nullptr, + my_rank == 0 ? dsp.data() : nullptr, + type, + 0, + MPI_COMM_WORLD); } +} - ik_begin_proc.resize(mympi->nprocs); - ik_end_proc.resize(mympi->nprocs); - ik_begin_proc[0] = 0; - ik_end_proc[0] = nk_proc[0]; - for (auto i = 1; i < mympi->nprocs; ++i) { - ik_begin_proc[i] = ik_end_proc[i - 1]; - ik_end_proc[i] = ik_begin_proc[i] + nk_proc[i]; +void PhononVelocity::calc_phonon_velmat_mesh(NDArray, 4> *velmat_out, + NDArray *velblock_out) const +{ + // velmat_out : full velocity matrix [nk][ns][ns][3] on rank 0 (needed only by the + // coherent term; nullptr to skip -- it is the memory hog). + // velblock_out: per-branch block-summed diad [nk][ns][3][3] on rank 0, + // velblock[k][is][a][b] = sum_{js in D(is)} Re(V^a_{is js} V^b_{js is}), + // so that summing it over the branches of one degenerate block gives + // the basis-invariant Tr(P V^a P V^b P). This is all the Peierls term + // and the boundary speed need, and it is O(nk ns) in memory, so the + // default RTA run never stores the full matrix. + if (!velmat_out && !velblock_out) return; + + const auto nk = dos->kmesh_dos->nk; + const auto ns = dynamical->neval; + const auto factor = Bohr_in_Angstrom * 1.0e-10 / (time_ry * 2.0 * pi); + const auto legacy = legacy_velocity(); + + if (mympi->my_rank == 0 && writes->getVerbosity() > 0) { + std::cout << " Calculating group velocity matrix of phonons on uniform grid ... "; } - std::vector klist_proc; - for (auto ik = ik_begin_proc[mympi->my_rank]; ik < ik_end_proc[mympi->my_rank]; ++ik) { - klist_proc.push_back(ik); + // k distribution + std::vector nk_proc(mympi->nprocs); + if (nk > static_cast(std::numeric_limits::max())) { + exit("calc_phonon_velmat_mesh", "Number of k points exceeds the supported range."); } + const auto nk_int = static_cast(nk); + auto nk_loc = nk_int / mympi->nprocs; + const auto nk_res = nk_int - nk_loc * mympi->nprocs; + for (auto i = 0; i < mympi->nprocs; ++i) nk_proc[i] = nk_loc + (i < nk_res ? 1 : 0); + auto ik_begin = 0; + for (auto i = 0; i < mympi->my_rank; ++i) ik_begin += nk_proc[i]; + nk_loc = nk_proc[mympi->my_rank]; + + NDArray, 3> vk; // one k point, discarded after use + NDArray, 4> velmat_loc; // only if the full matrix is wanted + NDArray velblock_loc; + vk.resize(ns, ns, 3); + if (velmat_out) velmat_loc.resize(std::max(nk_loc, 1), ns, ns, 3); + if (velblock_out) velblock_loc.resize(std::max(nk_loc, 1), ns, 3, 3); + + const auto &fc2_vel = + (!legacy && dynamical->nonanalytic == 3) ? ewald->fc2_without_dipole : fcs_phonon->force_constant_with_cell[0]; + const auto eval_all = dos->dymat_dos->get_eigenvalues(); + const auto evec_all = dos->dymat_dos->get_eigenvectors(); + const auto tol_cm = transport_block_tol_cm(); + std::vector blk_lo, blk_hi; - nk_loc = klist_proc.size(); + for (auto i = 0; i < nk_loc; ++i) { + const auto knum = ik_begin + i; + + // For NONANALYTIC = 3 the eigenproblem is solved with the dipole-free force + // constants plus an Ewald long-range matrix, so the velocity matrix has to be + // built from the same decomposition. + velocity_matrix_analytic(dos->kmesh_dos->xk[knum], fc2_vel, eval_all[knum], evec_all[knum], vk); + if (!legacy) add_nonanalytic_velocity_matrix(dos->kmesh_dos->xk[knum], eval_all[knum], evec_all[knum], vk); + + if (legacy) { + // Legacy elementwise little-group averaging treats elements as Cartesian + // vectors, which is valid only for non-degenerate diagonals and can + // suppress degenerate velocities. The default skips this averaging. + double symmetrizer_k[3][3]; + std::vector smallgroup_k; + kpoint->get_symmetrization_matrix_at_k(dos->kmesh_dos->xk[knum], smallgroup_k, symmetrizer_k); + for (auto j = 0u; j < ns; ++j) { + for (auto k = 0u; k < ns; ++k) rotvec(vk[j][k], vk[j][k], symmetrizer_k, 'T'); + } + } + for (auto j = 0u; j < ns; ++j) { + for (auto k = 0u; k < ns; ++k) { + rotvec(vk[j][k], vk[j][k], system->get_primcell().lattice_vector); + for (auto mu = 0; mu < 3; ++mu) vk[j][k][mu] *= factor; + } + } - velmat_loc.resize(nk_loc, ns, ns, 3); + if (velmat_out) { + for (auto j = 0u; j < ns; ++j) { + for (auto k = 0u; k < ns; ++k) { + for (auto mu = 0; mu < 3; ++mu) velmat_loc[i][j][k][mu] = vk[j][k][mu]; + } + } + } - for (auto i = 0; i < nk_loc; ++i) { - auto knum = klist_proc[i]; - velocity_matrix_analytic(dos->kmesh_dos->xk[knum], - fcs_phonon->force_constant_with_cell[0], - dos->dymat_dos->get_eigenvalues()[knum], - dos->dymat_dos->get_eigenvectors()[knum], - velmat_loc[i]); - - double symmetrizer_k[3][3]; - std::vector smallgroup_k; - kpoint->get_symmetrization_matrix_at_k(dos->kmesh_dos->xk[knum], smallgroup_k, symmetrizer_k); - - for (auto j = 0; j < ns; ++j) { - for (auto k = 0; k < ns; ++k) { - rotvec(velmat_loc[i][j][k], velmat_loc[i][j][k], symmetrizer_k, 'T'); - rotvec(velmat_loc[i][j][k], velmat_loc[i][j][k], system->get_primcell().lattice_vector); - for (auto mu = 0; mu < 3; ++mu) { - velmat_loc[i][j][k][mu] *= factor; + if (velblock_out) { + // Blocks taken at the REPRESENTATIVE of this k's star, from the same shared + // partition the coherent term uses to exclude same-block pairs, so the two + // partitions are identical by construction. + const auto irr = dos->kmesh_dos->kmap_to_irreducible[knum]; + const auto krep = dos->kmesh_dos->kpoint_irred_all[irr][0].knum; + transport_block_bounds(ns, eval_all[krep], tol_cm, blk_lo, blk_hi); + + for (auto is = 0u; is < ns; ++is) { + for (auto a = 0; a < 3; ++a) { + for (auto b = 0; b < 3; ++b) { + auto acc = 0.0; + for (auto js = blk_lo[is]; js < blk_hi[is]; ++js) { + acc += (vk[is][js][a] * vk[js][is][b]).real(); + } + velblock_loc[i][is][a][b] = acc; + } } } } } + vk.clear(); #ifdef MPI_CXX_DOUBLE_COMPLEX const auto mpi_complex_type = MPI_CXX_DOUBLE_COMPLEX; @@ -380,20 +562,26 @@ void PhononVelocity::calc_phonon_velmat_mesh(std::complex ****velmat_out const auto mpi_complex_type = MPI_COMPLEX16; #endif - MPI_Gatherv(nk_loc > 0 ? &velmat_loc[0][0][0][0] : nullptr, - sendcount[mympi->my_rank], - mpi_complex_type, - mympi->my_rank == 0 ? &velmat_out[0][0][0][0] : nullptr, - mympi->my_rank == 0 ? &recvcount[0] : nullptr, - mympi->my_rank == 0 ? &displs[0] : nullptr, - mpi_complex_type, - 0, - MPI_COMM_WORLD); - - velmat_loc.clear(); - sendcount.clear(); - recvcount.clear(); - displs.clear(); + if (velmat_out) { + gather_k_records>(nk_loc > 0 ? &velmat_loc[0][0][0][0] : nullptr, + nk_proc, + mympi->my_rank, + mympi->nprocs, + static_cast(ns) * ns * 3, + mpi_complex_type, + mympi->my_rank == 0 ? &(*velmat_out)[0][0][0][0] : nullptr); + velmat_loc.clear(); + } + if (velblock_out) { + gather_k_records(nk_loc > 0 ? &velblock_loc[0][0][0][0] : nullptr, + nk_proc, + mympi->my_rank, + mympi->nprocs, + static_cast(ns) * 9, + MPI_DOUBLE, + mympi->my_rank == 0 ? &(*velblock_out)[0][0][0][0] : nullptr); + velblock_loc.clear(); + } if (mympi->my_rank == 0 && writes->getVerbosity() > 0) { std::cout << "done!\n"; @@ -697,6 +885,174 @@ void PhononVelocity::calc_derivative_dynmat_k(const double *xk_in, const std::ve } } +// Central-difference D_na for NONANALYTIC methods 1/2/3. No unique +// direction-independent gradient exists at Gamma; no directional limit +// is implemented. Convert from the eigenproblem's cell-phase convention +// to the displacement-aware velocity convention using +// (grad_a D_na)_ij = (d_a D_na)_ij + i 2pi (t_j - t_i)_a (D_na)_ij, +// where t is the primitive fractional position. The connection term is +// needed even for diagonal velocities: commutator cancellation applies +// to the full dynamical matrix, not D_na alone. +void PhononVelocity::add_nonanalytic_velocity_matrix(const double *xk_in, const double *omega_in, + std::complex **evec_in, + std::complex ***velmat_inout, + const double *kvec_fixed) const +{ + if (dynamical->nonanalytic == 0) return; + + const auto nmode = dynamical->neval; + const auto h = 1.0e-4; + + // Drop the nonanalytic velocity at Gamma, where no direction-independent + // gradient exists. A central difference can be nonzero there; + // a directional limit is not implemented. + if (std::abs(xk_in[0]) < eps && std::abs(xk_in[1]) < eps && std::abs(xk_in[2]) < eps) { + if (mympi->my_rank == 0 && writes->getVerbosity() > 0) { + static auto warned_gamma = false; + if (!warned_gamma) { + warned_gamma = true; + warn("add_nonanalytic_velocity_matrix", + "Velocity at Gamma with NONANALYTIC != 0 is convention dependent; " + "the nonanalytic contribution is set to zero there."); + } + } + return; + } + + NDArray, 2> dna_plus, dna_minus; + NDArray, 3> ddna; + dna_plus.resize(nmode, nmode); + dna_minus.resize(nmode, nmode); + ddna.resize(nmode, nmode, 3); + + for (auto i = 0u; i < nmode; ++i) { + for (auto j = 0u; j < nmode; ++j) { + for (auto k = 0; k < 3; ++k) ddna[i][j][k] = std::complex(0.0, 0.0); + } + } + + double xk_shift[2][3], kvec[2][3]; + + for (auto idir = 0; idir < 3; ++idir) { + for (auto j = 0; j < 3; ++j) { + xk_shift[0][j] = xk_in[j]; + xk_shift[1][j] = xk_in[j]; + } + xk_shift[0][idir] -= h; + xk_shift[1][idir] += h; + + for (auto ishift = 0; ishift < 2; ++ishift) { + if (kvec_fixed) { + // Hold the band-segment direction kvec_na fixed so the derivative uses + // the same nonanalytic matrix as the eigenproblem. + for (auto j = 0; j < 3; ++j) kvec[ishift][j] = kvec_fixed[j]; + continue; + } + for (auto j = 0; j < 3; ++j) kvec[ishift][j] = xk_shift[ishift][j]; + rotvec(kvec[ishift], kvec[ishift], system->get_primcell().reciprocal_lattice_vector, 'T'); + const auto norm = std::sqrt(kvec[ishift][0] * kvec[ishift][0] + kvec[ishift][1] * kvec[ishift][1] + + kvec[ishift][2] * kvec[ishift][2]); + if (norm > eps) { + for (auto j = 0; j < 3; ++j) kvec[ishift][j] /= norm; + } + } + + auto &dst_minus = dna_minus; + auto &dst_plus = dna_plus; + + for (auto i = 0u; i < nmode; ++i) { + for (auto j = 0u; j < nmode; ++j) { + dst_minus[i][j] = std::complex(0.0, 0.0); + dst_plus[i][j] = std::complex(0.0, 0.0); + } + } + + if (dynamical->nonanalytic == 1) { + dynamical->calc_nonanalytic_k_parlinski(xk_shift[0], kvec[0], dst_minus); + dynamical->calc_nonanalytic_k_parlinski(xk_shift[1], kvec[1], dst_plus); + } else if (dynamical->nonanalytic == 2) { + dynamical->calc_nonanalytic_k_mixedspace(xk_shift[0], kvec[0], dst_minus); + dynamical->calc_nonanalytic_k_mixedspace(xk_shift[1], kvec[1], dst_plus); + } else if (dynamical->nonanalytic == 3) { + ewald->add_longrange_matrix(xk_shift[0], kvec[0], dst_minus); + ewald->add_longrange_matrix(xk_shift[1], kvec[1], dst_plus); + } + + for (auto i = 0u; i < nmode; ++i) { + for (auto j = 0u; j < nmode; ++j) { + ddna[i][j][idir] = (dst_plus[i][j] - dst_minus[i][j]) / (2.0 * h); + } + } + } + + // Connection term: derivative of the sublattice phase that calc_nonanalytic_k_* + // has already folded into D_na. Needs D_na at xk_in itself, not at the shifts. + NDArray, 2> dna0; + dna0.resize(nmode, nmode); + for (auto i = 0u; i < nmode; ++i) { + for (auto j = 0u; j < nmode; ++j) dna0[i][j] = std::complex(0.0, 0.0); + } + + double kvec0[3]; + if (kvec_fixed) { + for (auto j = 0; j < 3; ++j) kvec0[j] = kvec_fixed[j]; + } else { + for (auto j = 0; j < 3; ++j) kvec0[j] = xk_in[j]; + rotvec(kvec0, kvec0, system->get_primcell().reciprocal_lattice_vector, 'T'); + const auto norm0 = std::sqrt(kvec0[0] * kvec0[0] + kvec0[1] * kvec0[1] + kvec0[2] * kvec0[2]); + if (norm0 > eps) { + for (auto j = 0; j < 3; ++j) kvec0[j] /= norm0; + } + } + + if (dynamical->nonanalytic == 1) { + dynamical->calc_nonanalytic_k_parlinski(xk_in, kvec0, dna0); + } else if (dynamical->nonanalytic == 2) { + dynamical->calc_nonanalytic_k_mixedspace(xk_in, kvec0, dna0); + } else if (dynamical->nonanalytic == 3) { + ewald->add_longrange_matrix(xk_in, kvec0, dna0); + } + + const auto &xf_prim = system->get_primcell().x_fractional; + + for (auto i = 0u; i < nmode; ++i) { + for (auto j = 0u; j < nmode; ++j) { + for (auto k = 0; k < 3; ++k) { + const auto dt = xf_prim(j / 3, k) - xf_prim(i / 3, k); + ddna[i][j][k] += im * tpi * dt * dna0[i][j]; + } + } + } + dna0.clear(); + + // Project onto the eigenvectors at xk_in and apply Allen's normalisation. + // ddna is d(D_na)/d(q_fractional) already, so no extra factor of i here. + // E^H ddna E per direction: O(ns^3). + { + Eigen::MatrixXcd E(nmode, nmode); + for (auto i = 0u; i < nmode; ++i) { + for (auto j = 0u; j < nmode; ++j) E(j, i) = evec_in[i][j]; + } + Eigen::MatrixXcd Dk(nmode, nmode); + for (auto k = 0; k < 3; ++k) { + for (auto i = 0u; i < nmode; ++i) { + for (auto j = 0u; j < nmode; ++j) Dk(i, j) = ddna[i][j][k]; + } + const Eigen::MatrixXcd Vk = E.adjoint() * (Dk * E); + for (auto i = 0u; i < nmode; ++i) { + for (auto j = 0u; j < nmode; ++j) { + if (omega_in[i] < eps8 || omega_in[j] < eps8) continue; + velmat_inout[i][j][k] += Vk(i, j) * (0.5 / std::sqrt(omega_in[i] * omega_in[j])); + } + } + } + } + + dna_plus.clear(); + dna_minus.clear(); + ddna.clear(); +} + void PhononVelocity::velocity_matrix_analytic(const double *xk_in, const std::vector &fc2_in, const double *omega_in, std::complex **evec_in, std::complex ***velmat_out) const @@ -735,16 +1091,21 @@ void PhononVelocity::velocity_matrix_analytic(const double *xk_in, const std::ve } } - unsigned int ii, jj; - - for (i = 0; i < nmode; ++i) { - for (j = 0; j < nmode; ++j) { - for (ii = 0; ii < nmode; ++ii) { - for (jj = 0; jj < nmode; ++jj) { - for (k = 0; k < 3; ++k) { - velmat_out[i][j][k] += std::conj(evec_in[i][ii]) * ddymat[ii][jj][k] * evec_in[j][jj]; - } - } + // Project onto the eigenvectors as E^H (dD/dq) E with two matrix products per + // direction: O(ns^3), instead of the O(ns^4) explicit four-index sum. + { + Eigen::MatrixXcd E(nmode, nmode); + for (i = 0; i < nmode; ++i) { + for (j = 0; j < nmode; ++j) E(j, i) = evec_in[i][j]; // column i = eigenvector i + } + Eigen::MatrixXcd Dk(nmode, nmode); + for (k = 0; k < 3; ++k) { + for (i = 0; i < nmode; ++i) { + for (j = 0; j < nmode; ++j) Dk(i, j) = ddymat[i][j][k]; + } + const Eigen::MatrixXcd Vk = E.adjoint() * (Dk * E); + for (i = 0; i < nmode; ++i) { + for (j = 0; j < nmode; ++j) velmat_out[i][j][k] = Vk(i, j); } } } diff --git a/anphon/phonon_velocity.h b/anphon/phonon_velocity.h index e4bda6e9..edbb498c 100644 --- a/anphon/phonon_velocity.h +++ b/anphon/phonon_velocity.h @@ -27,6 +27,8 @@ class PhononVelocity: protected Pointers ~PhononVelocity(); + static bool legacy_velocity(); + void setup_velocity(); void phonon_vel_k(const double *, double **) const; @@ -36,6 +38,9 @@ class PhononVelocity: protected Pointers void get_phonon_group_velocity_mesh(const KpointMeshUniform &kmesh_in, const Eigen::Matrix3d &lavec_p, const bool irreducible_only, double ***phvel3_out) const; + void get_phonon_group_velocity_mesh_velmat(const KpointMeshUniform &kmesh_in, const Eigen::Matrix3d &lavec_p, + double ***phvel3_out) const; + void get_phonon_group_velocity_mesh_mpi(const KpointMeshUniform &kmesh_in, const Eigen::Matrix3d &lavec_p, double ***phvel3_out) const; @@ -43,7 +48,14 @@ class PhononVelocity: protected Pointers NDArray &vel_out, const double unit_factor, const bool bcast_full) const; - void calc_phonon_velmat_mesh(std::complex ****velmat_out) const; + // velmat_out (full matrix, coherent term only) and velblock_out (per-branch + // block-summed diad for the Peierls term / boundary speed) may each be nullptr. + void calc_phonon_velmat_mesh(NDArray, 4> *velmat_out, NDArray *velblock_out) const; + + void get_phonon_group_velocity_bandstructure_velmat(const KpointBandStructure *kpoint_bs_in, + const Eigen::Matrix3d &lavec_p, + const std::vector &fc2_in, + double **phvel_out) const; void get_phonon_group_velocity_bandstructure(const KpointBandStructure *kpoint_bs_in, const Eigen::Matrix3d &lavec_p, const Eigen::Matrix3d &rlavec_p, @@ -51,6 +63,12 @@ class PhononVelocity: protected Pointers const std::vector &fc2_without_dipole, double **phvel_out) const; + // kvec_fixed: hold the nonanalytic direction fixed (band paths, where the + // eigenproblem uses the segment direction). nullptr = radial, as on a mesh. + void add_nonanalytic_velocity_matrix(const double *xk_in, const double *omega_in, std::complex **evec_in, + std::complex ***velmat_inout, + const double *kvec_fixed = nullptr) const; + void velocity_matrix_analytic(const double *xk_in, const std::vector &fc2_in, const double *omega_in, std::complex **evec_in, std::complex ***velmat_out) const; diff --git a/anphon/relaxation.cpp b/anphon/relaxation.cpp index cdc430b6..7895bc98 100644 --- a/anphon/relaxation.cpp +++ b/anphon/relaxation.cpp @@ -820,16 +820,9 @@ void Relaxation::rescue_step_after_scp_failure(RelaxationStructureState &structu const std::vector &harm_optical_modes, double **omega2_harmonic, std::complex ***evec_harmonic) const { - // Called instead of update_cell_coordinate when the SCP equation did not - // converge at the current structure. The forces and stress evaluated from an - // unconverged SCP solution are unreliable, so they are not given to the - // optimizer: its history keeps only data from converged SCP solutions. - // The structure is moved back halfway along the last step, so repeated - // failures bisect toward the last structure where the SCP equation was - // solvable. If there is no previous step to undo (failure at the first - // structure iteration), a strongly damped steepest-descent step on the - // unreliable force is taken so that the optimization can leave the initial - // structure; the cell is kept fixed in that case. + // On SCP failure, keep unreliable forces/stress out of optimizer history + // and halve the last structure step. With no previous step, take a strongly + // damped force step at fixed cell to leave the initial structure. auto &q0 = structure_state.q0; auto &u0 = structure_state.u0; diff --git a/anphon/relaxation.h b/anphon/relaxation.h index 6974691f..9937d85f 100644 --- a/anphon/relaxation.h +++ b/anphon/relaxation.h @@ -226,17 +226,12 @@ class Relaxation: protected Pointers double alpha_steepest_decent; double cell_conv_tol; double mixbeta_cell; - // Optional residual-force convergence threshold for the internal coordinates (gradient - // w.r.t. q0). When > 0, structural optimization is declared converged only if the - // coordinate force norm is also below this value, in addition to the step-size criteria. - // Guards against false convergence (small step at a non-stationary point), which can - // occur with the GDIIS optimizer (relax_algo == 3). + // When positive, also require the q0 force norm below this threshold + // to prevent small-step convergence at a non-stationary point. double gradient_conv_tol; - // Optional residual convergence threshold for the cell gradient (the stress-like - // quantity conjugate to the strain tensor, including the applied-pressure term). When > 0 - // and the cell is relaxed (relax_str == 2), convergence also requires the strain-gradient - // norm to be below this value. Its units differ from gradient_conv_tol, hence a separate - // threshold (cf. coord_conv_tol vs cell_conv_tol). + // When positive and relax_str == 2, also require the cell-gradient norm + // (including pressure) below this threshold. Its units differ from the + // coordinate-force threshold. double cell_gradient_conv_tol; // For relax_algo == 3 (GDIIS): if nonzero, apply the Farkas-Schlegel "controlled GDIIS" // step-acceptance criteria (step-length cap, coefficient/extrapolation cap, and diff --git a/anphon/scph.cpp b/anphon/scph.cpp index 3504b9ad..5eb1c542 100644 --- a/anphon/scph.cpp +++ b/anphon/scph.cpp @@ -47,12 +47,9 @@ using namespace PHON_NS; namespace { -// Hermitian eigensolver for the SCPH iteration. Eigen's tridiagonal QR is the -// faster choice for small matrices; LAPACK zheevd (threaded MKL/OpenBLAS) wins -// from a few hundred modes on (n = 432: 0.07 s against 0.10 s single-threaded -// with OpenBLAS, and MKL threads over the 64 cores of a node). Eigenvectors in -// columns, eigenvalues ascending in both cases; the SCPH solver is invariant to -// the eigenvector phases. +// Use Eigen QR for small Hermitian problems and LAPACK zheevd for larger +// ones. Both return ascending eigenvalues and eigenvectors in columns; +// SCPH is invariant to eigenvector phases. constexpr int lapack_eigen_threshold = 256; void hermitian_eigen(const Eigen::MatrixXcd &mat, Eigen::VectorXd &eval, Eigen::MatrixXcd *evec) @@ -1337,12 +1334,9 @@ void Scph::compute_qmat_and_dmat(const Eigen::MatrixXd &omega_now, const double #pragma omp for for (int ik = 0; ik < static_cast(nk); ++ik) { - // At Gamma, the three translational (acoustic) modes must be excluded from Qmat. - // They are identified from the eigenvectors (via the overlap with the harmonic - // acoustic subspace encoded in cmat_convert), NOT from the frequency magnitude: - // a soft optical mode renormalized to nearly zero frequency would otherwise be - // silently treated as acoustic, which removes its own divergent (2n+1)/2omega - // self-interaction and thereby stabilizes a spurious omega = 0 solution. + // Exclude Gamma translations using overlap with the harmonic acoustic + // subspace in cmat_convert. A frequency cutoff would also exclude soft + // optical modes and spuriously stabilize a zero-frequency solution. std::vector is_acoustic_now; if (ik == ik_gamma_dense) { is_acoustic_now = classify_acoustic_modes_from_cmat(cmat_convert[ik]); @@ -1357,12 +1351,8 @@ void Scph::compute_qmat_and_dmat(const Eigen::MatrixXd &omega_now, const double // omega1 for any finite frequency. auto omega_for_q = std::abs(omega_now(ik, is)); if (omega_for_q < eps8) { - // A non-acoustic mode has (transiently) collapsed to zero frequency. - // Evaluating the 1/omega factor there would give a violent restoring - // kick that destabilizes the fixed-point iteration; instead evaluate it - // at a harmonic frequency scale (the value a cold-started iteration - // would use), which pushes the mode back to a finite frequency at a - // physically reasonable rate. + // Use a harmonic frequency scale for collapsed optical modes to avoid + // a destabilizing 1/omega factor and restore a finite frequency. if (ik == ik_gamma_dense) { omega_for_q = omega_floor_gamma; } else { @@ -1459,13 +1449,9 @@ void Scph::diagonalize_and_symmetrize(const Eigen::MatrixXcd &Fmat, const std::v std::cout << " onsite V4 is positive\n\n"; } - // With a warm start, the previous temperature's converged omega2 at the - // same sorted index is used as the reset target -- but only when it is - // meaningfully positive. The sorted index is not a branch label: at Gamma - // an imaginary soft mode sorts below the acoustic zeros, so the stored - // value can be (exactly) zero, and resetting to it would pin the soft mode - // at zero frequency for the rest of the iteration. In that case, and on - // cold starts, flip the sign of the eigenvalue instead. + // Reuse the previous temperature's omega2 only if positive. Sorted indices + // can match acoustic zeros rather than the same branch; otherwise flip + // the eigenvalue sign, as on a cold start. if (flag_converged && omega2_out[knum][is] > eps15) { ++icount; eval_tmp(is) = omega2_out[knum][is] * std::pow(0.99, icount); @@ -1867,20 +1853,10 @@ void Scph::compute_anharmonic_frequency_diis(double **omega2_out, std::complexelapsed(); - // SCPH iteration accelerated by Pulay/Anderson (DIIS) mixing. - // - // The fixed-point variable is the full set of D matrices on the dense k mesh, - // x = {D_k}, g(x) = K(omega(x), C(x)), - // and the residual handed to DIIS is r_n = g(x_n) - x_n of the same iterate. - // The update x_{n+1} = sum_m c_m (x_m + beta * r_m) with beta = mixalpha - // reduces exactly to the simple mixing of compute_anharmonic_frequency when - // the history holds a single pair, so the early iterations are identical to - // the reference implementation. - // - // A single DIIS history is kept for the concatenated state of all k points - // because the SCP equation couples the k points through the inner sum - // over q1; mixing D (rather than the eigenvalues) makes the residual - // basis-free, so no mode tracking across iterations is required. + // Pulay/Anderson mixing of all dense-mesh D matrices: + // r_n = g(x_n) - x_n, x_{n+1} = sum_m c_m (x_m + mixalpha * r_m). + // One history spans the coupled k points; mixing D avoids mode tracking. + // A single history pair reduces to simple mixing. using namespace Eigen; @@ -2098,14 +2074,9 @@ void Scph::compute_anharmonic_frequency_diis(double **omega2_out, std::complex 1) { std::cout << " DIIS: |r|/|x| = " << std::scientific << rnorm_rel; if (eval_repaired) std::cout << " (imaginary-mode repair active)"; diff --git a/anphon/scph_io.cpp b/anphon/scph_io.cpp index eca63480..ae6a6788 100644 --- a/anphon/scph_io.cpp +++ b/anphon/scph_io.cpp @@ -218,16 +218,10 @@ ScphFc2RowsH5 ScphQhaCommon::build_fc2_rows_h5(const std::complex *const const unsigned int NT, const KpointMeshUniform *kmesh_coarse_in, MinimumDistList ***mindist_list_in, const std::string &variant) const { - // Assemble the renormalized FC2 on the virtual supercell in the - // alamode force-constant schema. The row enumeration is identical to - // write_anharmonic_correction_fc2 above (one row per minimum-distance - // image, multiplicity-split); the base harmonic values come from the - // harmonic dynamical matrix sampled on the coarse mesh and transformed - // to real space through the same pathway as the anharmonic correction, - // so base and correction live on identical rows by construction. - // The folding is exact when KMESH_INTERPOLATE matches the supercell - // dimensions of the original FC2 (the standard setup); the imaginary - // part is dropped, as in the legacy .scph_dfc2 file. + // Build renormalized FC2 on the virtual supercell using the same + // minimum-distance, multiplicity-split rows for harmonic and correction IFCs. + // Folding is exact when KMESH_INTERPOLATE matches the original FC2 supercell; + // drop the imaginary part as in .scph_dfc2. const auto ns = dynamical->neval; const auto natmin = system->get_primcell().number_of_atoms; const auto nk1 = kmesh_coarse_in->nk_i[0]; diff --git a/anphon/scph_qha_common.cpp b/anphon/scph_qha_common.cpp index 07323154..53bba6f7 100644 --- a/anphon/scph_qha_common.cpp +++ b/anphon/scph_qha_common.cpp @@ -67,15 +67,9 @@ ScphQhaCommon::classify_acoustic_modes_from_cmat(const std::complex *con } } - // Majority overlap with the harmonic acoustic subspace marks a mode as acoustic. - // A threshold is used instead of picking the three largest overlaps on purpose: - // when a soft optical mode becomes numerically degenerate with the acoustic modes - // during the SCPH iteration, the eigensolver may return arbitrarily mixed columns, - // and a fixed count could then assign a mostly-translational column as "optical" - // (whose then huge 1/omega occupation factor would destabilize the iteration). - // With a threshold, every mostly-translational column is excluded, transiently - // mixed columns resolve themselves once the degeneracy is lifted, and in the - // clean (non-degenerate) case exactly the three translational modes are flagged. + // Flag modes with majority overlap with the harmonic acoustic subspace. + // A threshold handles mixed soft/acoustic eigenvectors without forcing + // exactly three exclusions and leaving a large 1/omega acoustic contribution. std::vector is_acoustic(ns, false); for (auto js = 0; js < ns; ++js) { is_acoustic[js] = overlap[js] > 0.5; @@ -1082,16 +1076,9 @@ void ScphQhaCommon::renormalize_ifcs_at_structure(StructuralOptWorkspace &ws) u_tensor); print_stage_time("strain renormalization v0..v3", time_stage); - // Renormalize the IFCs by the internal displacement q0 (exact Taylor - // recentering of the quartic PES). The strain-renormalized v1..v3 - // (_with_umn) enter here; v4 enters unrenormalized (the row-distributed - // reference V4 of v4_service) because its - // strain renormalization would require d(v4)/du IFC data, which - // del_v_strain does not include (it stops at d(v3)/du) -- within this - // truncation the strain-renormalized v4 equals the reference v4. - // - // v4 is swept once: the sweep writes v3_renorm and the quartic contraction - // q4_q0 that the v1, v2 and v0 renormalizations consume (q0_contraction.h). + // Taylor-recenter the quartic PES at q0 using strain-renormalized v1..v3. + // Use reference v4: del_v_strain stops at d(v3)/du. One v4 sweep builds + // v3_renorm and q4_q0 for the v1, v2, and v0 updates (q0_contraction.h). const auto ik_gamma_irred = static_cast(kmesh_coarse->kpoint_map_symmetry[0].knum_irred_orig); const auto time_sweep_start = timer->elapsed(); v4_service->q0_sweep(q0.data(), ws.v3_with_umn, ws.v3_renorm, ws.q4_q0); @@ -1194,14 +1181,10 @@ void ScphQhaCommon::build_v4_service(const bool full_tensor, const bool offdiag_ const auto nk2_prod = nk_irred * nk; const auto nprocs = static_cast(mympi->nprocs); - // The band-parallel builder computes every element in units of ns rows and - // distributes with the weighted unit partition; the k-point builder computes - // whole slices (two ns^2 x ns^2 scratch arrays) and distributes slices. The - // latter cannot distribute when there are fewer slices than ranks, and its - // scratch dwarfs the local rows when there are few slices per rank. - // With a single process the user's IALGO is kept as is (the k-point builder's - // 2 ns^4 scratch then costs up to 3x the V4 size for a Gamma-only mesh; IALGO = 1 - // avoids it), so that single-process results stay bitwise reproducible. + // The band builder distributes ns-row units; the k-point builder distributes + // whole slices and needs 2 ns^4 scratch. Prefer band distribution when + // slices are too few for MPI. Keep the user's IALGO on one process + // for bitwise reproducibility. auto band = full_tensor && use_band_parallel_v4(); if (full_tensor && !band && nprocs > 1) { const auto ns4 = static_cast(ns) * ns * ns * ns; @@ -1351,15 +1334,10 @@ void ScphQhaCommon::compute_and_print_step_gradients(const StructuralOptWorkspac { const auto ns = dynamical->neval; - // Residual gradient norms over the optimized degrees of freedom (the same - // gradients the optimizer acts on). A small step (du0/du_tensor) does not - // by itself imply a small gradient for the GDIIS optimizer - // (relax_algo == 3), so these are printed for diagnostics and, when the - // corresponding tolerance is > 0, also required for convergence (SCPH) to - // guard against false convergence at a non-stationary point. The - // coordinate force (gradient w.r.t. q0) and the cell gradient (stress - // conjugate to the strain tensor) have different units, so they are - // checked separately (cf. COORD_CONV_TOL vs CELL_CONV_TOL). + // Report residual norms and require them for SCPH convergence when the + // corresponding tolerance is positive. Coordinate and cell gradients + // have different units and use separate tolerances; small GDIIS steps + // alone do not establish stationarity. grad_norm = 0.0; for (auto is = 0; is < ns - 3; is++) { const double f = v1_eff[ws.harm_optical_modes[is]].real(); diff --git a/anphon/scph_qha_common.h b/anphon/scph_qha_common.h index edd5cbca..26c06363 100644 --- a/anphon/scph_qha_common.h +++ b/anphon/scph_qha_common.h @@ -89,11 +89,8 @@ class ScphQhaCommon: protected Pointers std::vector converged_scph_temp; std::vector converged_str_temp; - // Zeroth-order (static) potential energy V0(T) of the relaxed structure, - // recorded per temperature by the structural-optimization drivers and - // stored in the state file. Owned here (not by Relaxation) because it is - // per-run result state of the SCPH/QHA drivers. The drivers size it on - // every rank (exec entry) before any restart loader broadcasts into it. + // Static relaxed-structure energy V0(T), stored in the SCPH/QHA state file. + // Size on every rank before restart broadcasts. std::vector V0; // Legacy-text restart IO for V0 (PREFIX.V0). The text file also serves @@ -261,12 +258,8 @@ class ScphQhaCommon: protected Pointers static void build_cmat_at_k(unsigned int ns, const Eigen::MatrixXcd &evec_ref_mat, const std::complex *const *evec_new_at_k, std::complex **cmat_out); - // Identify which modes of the CURRENT (renormalized) eigenbasis at Gamma are the three - // translational (acoustic) modes, given the unitary C(k=Gamma) connecting the harmonic - // basis to the current one: overlap(js) = sum_{is in acoustic_harm} |C[is][js]|^2, and the - // three modes with the largest overlap are flagged. Robust against the reshuffling of the - // sorted mode indices that occurs when a soft optical mode becomes nearly degenerate with - // the acoustic modes during the SCPH iteration. + // Identify Gamma acoustic modes by majority overlap with the harmonic + // acoustic subspace: overlap(js) = sum_{is in acoustic_harm} |C[is][js]|^2. std::vector classify_acoustic_modes_from_cmat(const std::complex *const *cmat_at_gamma) const; // Occupation-weighted SCP mode matrix at dense k: diff --git a/anphon/scph_v3v4_elements.cpp b/anphon/scph_v3v4_elements.cpp index bb6c3759..4e0061eb 100644 --- a/anphon/scph_v3v4_elements.cpp +++ b/anphon/scph_v3v4_elements.cpp @@ -32,13 +32,9 @@ using PHON_NS::v4_index_transform::transform_index_gemm; namespace { -// CSC-like skeleton of the quartic-IFC scatter pattern phi4[(a1,a2)][(a3,a4)]. -// The row/column indices depend only on evec_index_v4, which is fixed once in -// AnharmonicCore::setup_quartic(), so the skeleton is built once per kernel call -// and only the values are refilled for each (k1,k2) pair. Duplicate (row,col) -// pairs are merged into one slot; refilling in increasing group order keeps the -// last group's value, reproducing the dense scatter's overwrite semantics (the -// quadruplets are in fact unique after grouping, so this is a safe no-op). +// Cache the CSC-like scatter pattern phi4[(a1,a2)][(a3,a4)] from +// evec_index_v4 and refill values per (k1,k2). Merge duplicate slots, +// keeping the last group's value to match dense-scatter semantics. struct SparsePhi4Skeleton { std::vector row; // size nnz: row index a1*ns+a2 @@ -608,11 +604,9 @@ void ScphQhaCommon::compute_V4_elements_mpi_over_kpoint(v4_distributed::V4RowBlo } } - // The remaining transforms are matrix products: each one contracts the - // outermost mode index of the (ns x ns^3) view of the current buffer and - // appends the new index innermost, so the layout rotates as - // [a2 a3 a4 i] -> [a3 a4 i j] -> [a4 i j k] -> [i j k m], and after the - // fourth transform the flat layout coincides with the owned slice of v4. + // Contract outermost indices with GEMMs on ns x ns^3 views: + // [a2 a3 a4 i] -> [a3 a4 i j] -> [a4 i j k] -> [i j k m]. + // The final layout matches the owned v4 slice. constexpr auto complex_one = std::complex(1.0, 0.0); // transform the second index (v4_tmp1 -> v4_tmp2) @@ -828,11 +822,8 @@ void ScphQhaCommon::compute_V4_elements_mpi_over_band(v4_distributed::V4RowBlock } } - // The remaining transforms are matrix products on the (ns x ns^2) view of - // the buffer: each contracts the outermost mode index and appends the new - // index innermost, so the layout rotates [a2 a3 a4] -> [a3 a4 j] -> [a4 j k] - // -> [j k l]. The last one writes, with the prefactor, straight into the ns - // contiguous rows (is_now, j) of the owned unit, columns k*ns + l. + // Contract ns x ns^2 views: [a2 a3 a4] -> [a3 a4 j] -> [a4 j k] -> [j k l]. + // Write the scaled result to owned rows (is_now, j), columns k*ns + l. constexpr auto complex_one = std::complex(1.0, 0.0); transform_index_gemm(&evec_in[knum][0][0], &v4_tmp1[0][0], &v4_tmp2[0][0], ns, ns2, complex_one); transform_index_gemm(&evec_in[jk_now][0][0], &v4_tmp2[0][0], &v4_tmp1[0][0], ns, ns2, complex_one); diff --git a/anphon/system.cpp b/anphon/system.cpp index 3e746d00..1eb35f39 100644 --- a/anphon/system.cpp +++ b/anphon/system.cpp @@ -1211,12 +1211,8 @@ void System::set_mass_elem_from_database(const unsigned int nkd, const std::vect void System::set_atomtype_group(const Cell &cell_in, const Spin &spin_in, std::vector> &atomtype_group_out) { - // In the case of collinear calculation, spin moments are considered as scalar - // variables. Therefore, the same elements with different magnetic moments are - // considered as different types. In noncollinear calculations, - // magnetic moments are not considered in this stage. They will be treated - // separately in symmetry.cpp where spin moments will be rotated and flipped - // using time-reversal symmetry. + // Collinear moments distinguish atom types. Noncollinear moments are + // handled in symmetry.cpp with spin rotations and time reversal. unsigned int i; AtomType type_tmp{}; diff --git a/anphon/v4_distributed.h b/anphon/v4_distributed.h index f27b19db..f859f296 100644 --- a/anphon/v4_distributed.h +++ b/anphon/v4_distributed.h @@ -42,13 +42,9 @@ namespace PHON_NS::v4_distributed { -// Contiguous partition of weighted units over nprocs ranks. Returns the -// nprocs + 1 boundaries; rank r owns [bounds[r], bounds[r+1]). The boundary -// after rank r is the prefix whose cumulative weight is nearest to -// (r + 1) / nprocs of the total (ties go to the earlier prefix), so every -// rank's share differs from the ideal one by at most half a unit's weight on -// each side. Ranks may end up empty (anywhere in the sequence) when there are -// fewer units than ranks. +// Partition weighted units contiguously; rank r owns [bounds[r], bounds[r+1]). +// Choose prefixes nearest each target cumulative share, breaking ties toward +// earlier prefixes. Return nprocs+1 bounds; ranks may be empty. inline std::vector partition_units(const std::vector &weight, const int nprocs) { const std::size_t nunits = weight.size(); diff --git a/anphon/write_phonons.cpp b/anphon/write_phonons.cpp index 1a736171..286b1391 100644 --- a/anphon/write_phonons.cpp +++ b/anphon/write_phonons.cpp @@ -9,6 +9,17 @@ or http://opensource.org/licenses/mit-license.php for information. */ #include "write_phonons.h" +#include "phonon_velocity.h" + +namespace +{ +// PRINTVEL follows the transport velocity formulation: matrix diagonal by default, +// finite differences under the legacy opt-out. +bool use_velmat_velocities() +{ + return !PHON_NS::PhononVelocity::legacy_velocity(); +} +} // namespace #include #include #include "anharmonic_core.h" @@ -674,12 +685,23 @@ void Writes::writePhononVel() const NDArray phvel_bs; phvel_bs.resize(nk, dynamical->neval); - phonon_velocity->get_phonon_group_velocity_bandstructure(kpoint->kpoint_bs.get(), - system->get_primcell().lattice_vector, - system->get_primcell().reciprocal_lattice_vector, - fcs_phonon->force_constant_with_cell[0], - ewald->fc2_without_dipole, - phvel_bs); + // Same velocity machinery as the transport terms. This makes + // the printed velocities come from the same source; it does NOT make them reproduce + // the conductivity, which treats degenerate multiplets as blocks having no per-mode + // velocity. Printed values at a degeneracy remain one admissible basis choice. + if (use_velmat_velocities()) { + phonon_velocity->get_phonon_group_velocity_bandstructure_velmat(kpoint->kpoint_bs.get(), + system->get_primcell().lattice_vector, + fcs_phonon->force_constant_with_cell[0], + phvel_bs); + } else { + phonon_velocity->get_phonon_group_velocity_bandstructure(kpoint->kpoint_bs.get(), + system->get_primcell().lattice_vector, + system->get_primcell().reciprocal_lattice_vector, + fcs_phonon->force_constant_with_cell[0], + ewald->fc2_without_dipole, + phvel_bs); + } ofs_vel << "# k-axis, |Velocity| [m / sec]\n"; ofs_vel.setf(std::ios::fixed); @@ -732,10 +754,16 @@ void Writes::writePhononVelAll() const phvel.resize(nk, ns); phvel_xyz.resize(nk, ns, 3); - phonon_velocity->get_phonon_group_velocity_mesh(*dos->kmesh_dos.get(), - system->get_primcell().lattice_vector, - false, - phvel_xyz); + if (use_velmat_velocities()) { + phonon_velocity->get_phonon_group_velocity_mesh_velmat(*dos->kmesh_dos.get(), + system->get_primcell().lattice_vector, + phvel_xyz); + } else { + phonon_velocity->get_phonon_group_velocity_mesh(*dos->kmesh_dos.get(), + system->get_primcell().lattice_vector, + false, + phvel_xyz); + } unsigned int ik, is; #ifdef _OPENMP #pragma omp parallel for private(is) @@ -2105,14 +2133,9 @@ void Writes::writeNewFcsXml(const std::string &filename_xml, const std::vector 3 * it.atoms_s[k + 1] + it.pairs[k + 1].index % 3) { diff --git a/docs/source/anphondir/formalism_anphon.rst b/docs/source/anphondir/formalism_anphon.rst index ebc00d6c..82ddae00 100644 --- a/docs/source/anphondir/formalism_anphon.rst +++ b/docs/source/anphondir/formalism_anphon.rst @@ -158,14 +158,70 @@ The group velocity of phonon mode :math:`\boldsymbol{q}j` is given by \boldsymbol{v}_{\boldsymbol{q}j} = \frac{\partial \omega_{\boldsymbol{q}j}}{\partial \boldsymbol{q}}. -To evaluate the group velocity numerically, we employ a central difference where -:math:`\boldsymbol{v}` may approximately be given by +*anphon* evaluates it from the **velocity matrix** (default) or by a **finite difference** (legacy). + +Velocity matrix (default) +~~~~~~~~~~~~~~~~~~~~~~~~~ + +Let :math:`\tilde{D}(\boldsymbol{q})` denote the dynamical matrix whose phase factor carries the full interatomic +separation rather than the lattice vector alone, + +.. math:: + :label: dymat_tilde + + \tilde{D}_{\mu\nu}(\kappa\kappa^{\prime};\boldsymbol{q}) = \frac{1}{\sqrt{M_{\kappa}M_{\kappa^{\prime}}}} + \sum_{\ell^{\prime}}\Phi_{\mu\nu}(\ell\kappa;\ell^{\prime}\kappa^{\prime}) + \exp{\left[i\boldsymbol{q}\cdot(\boldsymbol{r}(\ell^{\prime}\kappa^{\prime})-\boldsymbol{r}(\ell\kappa))\right]}. + +It is related to :eq:`dymat` by the unitary transformation +:math:`\tilde{D}(\boldsymbol{q}) = U^{\dagger}(\boldsymbol{q})D(\boldsymbol{q})U(\boldsymbol{q})` with +:math:`U_{\kappa\kappa^{\prime}}(\boldsymbol{q})=\delta_{\kappa\kappa^{\prime}}\,e^{i\boldsymbol{q}\cdot\boldsymbol{r}(\kappa)}`, +so it has the same eigenvalues :math:`\omega_{\boldsymbol{q}j}^{2}`, and its eigenvectors are +:math:`\tilde{\boldsymbol{e}}_{\boldsymbol{q}j}=U^{\dagger}(\boldsymbol{q})\boldsymbol{e}_{\boldsymbol{q}j}`. +The band off-diagonal generalization of the group velocity [9]_ [11]_ is .. math:: + :label: velmat - \boldsymbol{v}_{\boldsymbol{q}j} \approx \frac{\omega_{\boldsymbol{q}+\Delta\boldsymbol{q}j} - \omega_{\boldsymbol{q}-\Delta\boldsymbol{q}j}}{2\Delta\boldsymbol{q}}. + v_{\boldsymbol{q}jj'}^{\mu} = \frac{1}{2\sqrt{\omega_{\boldsymbol{q}j}\omega_{\boldsymbol{q}j'}}} + (\tilde{\boldsymbol{e}}_{\boldsymbol{q}j}^{*})^{\mathrm{T}} + \frac{\partial \tilde{D}(\boldsymbol{q})}{\partial q_{\mu}} \tilde{\boldsymbol{e}}_{\boldsymbol{q}j'}, -If one needs to save the group velocities, please turn on the ``PRINTVEL``-tag. +whose diagonal :math:`v_{\boldsymbol{q}jj}^{\mu}` is the group velocity. It is :math:`\tilde{D}`, not :math:`D`, that must +be differentiated: the two derivatives differ by :math:`i\,[\,D,\,\mathrm{diag}(r_{\mu}(\kappa))\,]`, whose matrix elements +between eigenvectors carry the factor :math:`\omega_{\boldsymbol{q}j}^{2}-\omega_{\boldsymbol{q}j'}^{2}`. The diagonal and any +element inside a degenerate multiplet are therefore the same for both, but genuinely off-diagonal elements --- those entering the +:ref:`coherent term ` --- are not (the *displacement-aware* convention of Ref. [11]_). + +For polar systems the non-analytic correction (``NONANALYTIC = 1, 2, 3``) contributes to +:math:`\partial\tilde{D}/\partial\boldsymbol{q}`. It is obtained by a central difference of the assembled non-analytic +matrix together with the derivative of the sublattice phase factor relating :math:`D` and :math:`\tilde{D}`; omitting the +latter corrupts even the diagonal velocities. At :math:`\Gamma` the non-analytic velocity is not defined without a +directional convention and is set to zero. On a band path the derivative uses the same direction vector as the eigenproblem +for that segment. + +Inside a degenerate multiplet :math:`\mathcal{B}` the eigenvectors are fixed only up to a unitary rotation, so the individual +:math:`v_{\boldsymbol{q}jj}^{\mu}` are not defined; only block traces such as +:math:`\sum_{j,j'\in\mathcal{B}} v_{\boldsymbol{q}jj'}^{\mu}v_{\boldsymbol{q}j'j}^{\nu}` are. How the transport terms use this is described +in :ref:`the Peierls-term section `. + +Finite difference (legacy) +~~~~~~~~~~~~~~~~~~~~~~~~~~ + +Historically the group velocity was obtained from a central difference, + +.. math:: + + \boldsymbol{v}_{\boldsymbol{q}j} \approx \frac{\omega_{\boldsymbol{q}+\Delta\boldsymbol{q}j} - \omega_{\boldsymbol{q}-\Delta\boldsymbol{q}j}}{2\Delta\boldsymbol{q}}, + +where :math:`j` is the index in the *sorted* eigenvalue list at each shifted point. Wherever branches cross or are degenerate the sorted +index exchanges character between :math:`\boldsymbol{q}\pm\Delta\boldsymbol{q}` and the quotient connects two different branches. +This path is retained for comparison (environment variable ``ALAMODE_LEGACY_VELOCITY=1``, which restores the previous +velocity treatment throughout) and is still used by the adaptive smearing widths (``ISMEAR = 2``) and by the iterative +Boltzmann solvers, which have not been reformulated. + +If one needs to save the group velocities, please turn on the ``PRINTVEL``-tag; the printed values follow the same formulation as the +conductivity in that run. At a degeneracy they are one admissible basis choice, not a unique value. Mode effective charge @@ -438,6 +494,8 @@ The average mass :math:`M_{\kappa}` is substituted by the value specified in the .. _kappa: +.. _kappa_peierls: + Lattice thermal conductivity (Peierls term) ------------------------------------------- @@ -445,15 +503,65 @@ The lattice thermal conductivity tensor :math:`\kappa_{\mathrm{ph}}^{\mu\nu}(T)` .. math:: - \kappa_{\mathrm{ph}}^{\mu\nu}(T) = \frac{1}{V N_{q}} \sum_{\boldsymbol{q},j}c_{\boldsymbol{q}j}(T)v_{\boldsymbol{q}j}^{\mu}v_{\boldsymbol{q}j}^{\nu}\tau_{\boldsymbol{q}j}(T), + \kappa_{\mathrm{ph}}^{\mu\nu}(T) = \frac{1}{V N_{q}} \sum_{\boldsymbol{q},j}c_{\boldsymbol{q}j}(T)\,W_{\boldsymbol{q}j}^{\mu\nu}\,\tau_{\boldsymbol{q}j}(T), + \qquad + W_{\boldsymbol{q}j}^{\mu\nu} = \sum_{j'\in\mathcal{B}(j)} v_{\boldsymbol{q}jj'}^{\mu}v_{\boldsymbol{q}j'j}^{\nu}, + +where :math:`V` is the unit cell volume, :math:`c_{\boldsymbol{q}j} = \hbar\omega_{\boldsymbol{q}j}\partial n_{\boldsymbol{q}j}/\partial T`, :math:`\tau_{\boldsymbol{q}j}(T)` is the phonon lifetime, +and :math:`\mathcal{B}(j)` is the set of branches degenerate with :math:`j` at :math:`\boldsymbol{q}` (see below). +For a non-degenerate branch :math:`W_{\boldsymbol{q}j}^{\mu\nu} = v_{\boldsymbol{q}j}^{\mu}v_{\boldsymbol{q}j}^{\nu}` and the familiar expression is recovered. +For a degenerate multiplet the sum of :math:`W` over its members is + +.. math:: + + \sum_{j\in\mathcal{B}} W_{\boldsymbol{q}j}^{\mu\nu} + = \sum_{j,j'\in\mathcal{B}} v_{\boldsymbol{q}jj'}^{\mu}v_{\boldsymbol{q}j'j}^{\nu} + = \mathrm{Tr}\left(P_{\mathcal{B}}V^{\mu}P_{\mathcal{B}}V^{\nu}\right), + +where :math:`V^{\mu}` is the matrix with elements :math:`v_{\boldsymbol{q}jj'}^{\mu}` over all branches and +:math:`P_{\mathcal{B}}` is the projector onto the multiplet, i.e. in the band basis the diagonal matrix equal to 1 for +:math:`j\in\mathcal{B}` and 0 otherwise. A unitary rotation :math:`\mathcal{W}` of the eigenvectors inside the multiplet maps +:math:`V^{\mu}\to\mathcal{W}^{\dagger}V^{\mu}\mathcal{W}` and leaves :math:`P_{\mathcal{B}}` unchanged, so the trace is +invariant, whereas the individual products :math:`v_{\boldsymbol{q}jj}^{\mu}v_{\boldsymbol{q}jj}^{\nu}` are not. The pairs :math:`j\neq j'` inside one multiplet are counted here and are +excluded from the :ref:`coherent term `, so nothing is double counted; the two together form the invariant. + +This assignment is a convention, not a physical distinction. At exact degeneracy the coherent expression of Ref. [11]_ reduces +to the Peierls form for such a pair, so the total conductivity is the same whichever term the pair is assigned to; the assignment +only makes each part separately basis invariant. The *particle-like* part defined here therefore equals the *populations* of +Ref. [11]_ plus the coherences between degenerate branches, and differs from their populations wherever multiplets occur. For +symmetry-enforced degeneracies the multiplet is a well-defined object and the definition is exact; for accidental near-degeneracies +it depends on the numerical tolerance below, whose effect on the total is the factor reported by the warning described there. +Only the total should be compared between calculations or with Ref. [11]_. + +The multiplets are detected from the frequencies only. The frequency-degeneracy groups over which the self-energy is averaged +(tolerance :math:`10^{-7}` Ry) are subdivided with an anchored tolerance of :math:`10^{-6}\ \mathrm{cm}^{-1}`, so that the lifetime is +constant within every multiplet by construction. This is an empirical numerical criterion: a pair whose true splitting +:math:`\Delta\omega` falls below it is treated as degenerate, which overestimates its weight by the factor +:math:`1+[\Delta\omega/(\Gamma_{\boldsymbol{q}j}+\Gamma_{\boldsymbol{q}j'})]^{2}` relative to the coherent expression. *anphon* reports the +largest such ratio among the merged multiplets and warns when it exceeds 0.1. + +Because :math:`W` is contracted for each :math:`\boldsymbol{q}` as soon as the velocity matrix at that point is formed, the full matrix +(of size :math:`N_{q}(3N_{\kappa})^{2}\times3`) is never stored unless the coherent term is requested. -where :math:`V` is the unit cell volume, :math:`c_{\boldsymbol{q}j} = \hbar\omega_{\boldsymbol{q}j}\partial n_{\boldsymbol{q}j}/\partial T`, and :math:`\tau_{\boldsymbol{q}j}(T)` is the phonon lifetime. The phonon lifetime is estimated using the Matthiessen's rule as .. math:: \tau_{\boldsymbol{q}j}^{-1}(T) = 2 (\Gamma_{\boldsymbol{q}j}^{\mathrm{anh}}(T) + \Gamma_{\boldsymbol{q}j}^{\mathrm{iso}}). +When boundary scattering is included (``LEN_BOUNDARY``), the additional rate :math:`|\boldsymbol{v}_{\boldsymbol{q}j}|/L` uses, for a +degenerate multiplet, the block mean speed :math:`|\boldsymbol{v}|^{2}=d_{\mathcal{B}}^{-1}\sum_{j\in\mathcal{B}}\sum_{\mu}W_{\boldsymbol{q}j}^{\mu\mu}`, +which is likewise basis invariant. + +.. note:: + + With this formulation the conductivity does not depend on the choice of unit cell: a primitive cell and a commensurate + supercell sampling the same :math:`\boldsymbol{q}` points give identical :math:`\kappa` (to :math:`10^{-15}` relative in a + silicon test) when a fixed-width smearing (``ISMEAR = 0`` or ``1``) is used. The tetrahedron method (``ISMEAR = -1``) + interpolates frequencies at fixed sorted branch indices and its linewidths remain cell dependent at the few-percent level; the + adaptive smearing (``ISMEAR = 2``) and the iterative solvers (``SOLVER = IBTE`` and variants) still use finite-difference + velocities and are cell dependent as well. + The lattice thermal conductivity is saved in the file ``PREFIX``.kl. The spectra of the lattice thermal conductivity :math:`\kappa_{\mathrm{ph}}^{\mu\mu}(\omega)` can also be calculated by setting ``KAPPA_SPEC = 1`` in the ``&analysis`` field. :math:`\kappa_{\mathrm{ph}}^{\mu\mu}(\omega)` is defined as @@ -475,6 +583,11 @@ The accumulative lattice thermal conductivity :math:`\kappa_{\mathrm{ph,acc}}^{\ \kappa_{\mathrm{ph,acc}}^{\mu\mu}(L) = \frac{1}{V N_{q}} \sum_{\boldsymbol{q},j}c_{\boldsymbol{q}j}v_{\boldsymbol{q}j}^{\mu}v_{\boldsymbol{q}j}^{\mu}\tau_{\boldsymbol{q}j}\Theta (L-|\boldsymbol{v}_{\boldsymbol{q}j}|\tau_{\boldsymbol{q}j}), where :math:`\Theta(x)` is the step function. This quantity can be calculated by using the script ``analyzer.py`` with ``--calc cumulative`` flag. +With the default velocity-matrix formulation ``analyzer.py`` replaces :math:`v^{\mu}_{\boldsymbol{q}j}v^{\nu}_{\boldsymbol{q}j}` +by the degenerate-block diad :math:`W^{\mu\nu}_{\boldsymbol{q}j}` stored in ``PREFIX.kappa.h5`` (and +:math:`|\boldsymbol{v}_{\boldsymbol{q}j}|` by :math:`\sqrt{\mathrm{tr}\,W_{\boldsymbol{q}j}}`), so that the +large-:math:`L` limit coincides with :math:`\kappa_{\mathrm{P}}` of the ``anphon`` run; files written by older +versions fall back to the finite-difference velocities. One can also use another definition for the accumulative thermal conductivity: .. math:: @@ -497,8 +610,17 @@ The coherent components of lattice thermal conductivity (see Ref. [8]_), which a where :math:`c_{\boldsymbol{q}j} = \hbar\omega_{\boldsymbol{q}j}\partial n_{\boldsymbol{q}j}/\partial T` and :math:`\Gamma_{\boldsymbol{q}j}` is the total phonon linewidth (half width) of phonon mode :math:`\boldsymbol{q}j`. -:math:`\boldsymbol{v}_{\boldsymbol{q}jj'}` is a band off-diagonal generalization of the group velocity [9]_. When ``KAPPA_COHERENT = 1 | 2`` the coherent component is calculated and saved in ``PREFIX``.kl_coherent. When ``KAPPA_COHERENT = 2``, all components of the coherent term before summation are saved in ``PREFIX``.kc_elem. - +:math:`\boldsymbol{v}_{\boldsymbol{q}jj'}` is the band off-diagonal velocity matrix :eq:`velmat`, built from +:math:`\tilde{D}` [9]_ [11]_. The sum runs over pairs belonging to *different* degenerate multiplets; +pairs inside one multiplet are already contained in the :ref:`Peierls term `. For such a pair +:math:`\omega_{\boldsymbol{q}j}=\omega_{\boldsymbol{q}j'}`, the prefactor reduces to :math:`c_{\boldsymbol{q}j}` and the Lorentzian +factor to :math:`1/(\Gamma_{\boldsymbol{q}j}+\Gamma_{\boldsymbol{q}j'})`; since the lifetime is constant within a multiplet +this equals :math:`1/(2\Gamma_{\boldsymbol{q}j})=\tau_{\boldsymbol{q}j}`, i.e. exactly the Peierls weight +:math:`c_{\boldsymbol{q}j}v_{\boldsymbol{q}jj'}^{\mu}v_{\boldsymbol{q}j'j}^{\nu}\tau_{\boldsymbol{q}j}`, so the two terms together +form the basis-invariant block trace. The particle-like/wave-like *split* is therefore basis dependent while their sum is not; only the total +should be compared between calculations. When ``KAPPA_COHERENT = 1 | 2`` the coherent component is calculated and saved in ``PREFIX``.kl_coherent. +When ``KAPPA_COHERENT = 2``, all components of the coherent term before summation are saved in ``PREFIX``.kc_elem. Requesting the coherent +term stores the full velocity matrix, :math:`N_{q}(3N_{\kappa})^{2}\times3` complex numbers on the root process. Delta function -------------- @@ -608,5 +730,7 @@ When ``SELF_OFFDIAG = 1``, the off-diagonal elements are also calculated, and th .. [9] P\. B. Allen and J. L. Feldman, Phys. Rev. B **48**, 12581 (1993). .. [10] X\. Gonze and C. Lee, Phys. Rev. B **55**, 10355 (1997). + +.. [11] M\. Simoncelli, N. Marzari, and F. Mauri, Phys. Rev. X **12**, 041011 (2022). diff --git a/docs/source/anphondir/inputanphon.rst b/docs/source/anphondir/inputanphon.rst index 7a6a1101..0a45722b 100644 --- a/docs/source/anphondir/inputanphon.rst +++ b/docs/source/anphondir/inputanphon.rst @@ -1979,6 +1979,14 @@ Please specify the initial atomic displacements :math:`u^{(0)}_{\alpha \mu}` [Bo ``&general`` field for backward compatibility (deprecated); the ``&kappa`` value wins when both are given. + The file records the velocity formulation used to assemble + :math:`\kappa` (HDF5 attribute ``/kappa/formulation``). A file + written by an older version can be restarted under the current + formulation: the stored three-phonon linewidths do not depend on it + and :math:`\kappa` is reassembled from them. Only a + temperature-resolved file holding :math:`\kappa` for temperatures the + new run does not recompute is refused, to avoid mixing formulations. + ```` .. _anphon_restart_4ph: @@ -2014,6 +2022,9 @@ Please specify the initial atomic displacements :math:`u^{(0)}_{\alpha \mu}` [Bo :Default: 0 :Type: Integer :Description: This flag is available when ``MODE = kappa``. For the theoretical details, please see :ref:`this page `. + The Peierls term uses the velocity matrix regardless of this flag; setting it to 1 or 2 additionally + stores the full band off-diagonal matrix, which costs :math:`N_{q}(3N_{\kappa})^{2}\times 3` complex + numbers on the root process. .. caution:: diff --git a/docs/source/anphondir/outputanphon.rst b/docs/source/anphondir/outputanphon.rst index 34521cf2..31d91fdf 100644 --- a/docs/source/anphondir/outputanphon.rst +++ b/docs/source/anphondir/outputanphon.rst @@ -128,7 +128,12 @@ ANPHON: Output files Unified, crash-safe HDF5 result file of ``MODE = kappa`` (schema ``alamode:kappa_result``). It stores the run metadata, phonon frequencies, - group velocities, the three-phonon (and, when ``QUARTIC = 1``, four-phonon) + group velocities (``/scattering/3ph/velocities``; with the default + velocity-matrix formulation also the degenerate-block velocity diad + :math:`W^{\mu\nu}_{\boldsymbol{q}j}` of :ref:`kappa_peierls` as + ``/scattering/3ph/velocity_diad``, which ``analyzer.py`` uses so that its + cumulative and boundary-limited kappa reproduce ``kappa_peierls``), the + three-phonon (and, when ``QUARTIC = 1``, four-phonon) linewidths with per-mode completion flags, the isotope-scattering linewidths and factors (``/scattering/isotope/gamma``, ``/metadata/isotope_factors``, when ``ISOTOPE > 0``), and the final thermal-conductivity diff --git a/tools/GenDisplacement.py b/tools/GenDisplacement.py index b0c509eb..e9d5892a 100644 --- a/tools/GenDisplacement.py +++ b/tools/GenDisplacement.py @@ -859,9 +859,8 @@ def _get_gaussian_sigma(self, temp, ignore_imag, softmode_only=False): else: sigma[iq, imode] = 0.0 - # 5) Option 1: cap σ so typical displacement (σ*C) per mode - # doesn't exceed nn-fraction. Note: truncation at ±k*σ may - # occasionally allow slightly larger displacements. + # Cap typical per-mode displacement sigma*C at nn-fraction; + # truncation at +/-k*sigma can still allow larger displacements. if self._cap_sigma_by_nn: dnn = self._nearest_neighbor_distances() # Å, per atom Cgain = self._per_mode_realspace_gain() # Å per unit Q @@ -1221,10 +1220,7 @@ def _n_bose(self, omega, temperature): else: temperature_au = self._K_BOLTZMANN * temperature / self._RYDBERG_TO_JOULE x = omega / temperature_au - # Prevent overflow when x is too large - # (e.g., when omega=1e6 for acoustic mode suppression) - # For x > ~700, exp(x) would overflow. - # At large x, n_BE ≈ exp(-x) ≈ 0 + # For large x, n_BE approaches zero; avoid exp(x) overflow above ~700. if x > 100.0: return 0.0 return 1.0 / (math.exp(x) - 1.0) @@ -1462,9 +1458,7 @@ def _identify_gamma_acoustic(self, overlap_thresh=0.95, max_modes=3, verbose=Tru if np.linalg.norm(disp_real[jat]) > 1e-10: max_rel_diff = 1.0 # Different from zero reference - # Uniformity score: 1 - max_relative_difference - # Perfect acoustic mode: max_rel_diff = 0 → score = 1.0 - # Non-uniform mode: max_rel_diff → ∞ → score → 0 + # Uniformity score: 1 for equal displacements, approaching 0 as they diverge. uniformity_scores[imode] = 1.0 / (1.0 + max_rel_diff) # Mark modes with high uniformity score as acoustic diff --git a/tools/analyzer/anphonio.py b/tools/analyzer/anphonio.py index 44a61777..52d7b74b 100644 --- a/tools/analyzer/anphonio.py +++ b/tools/analyzer/anphonio.py @@ -19,6 +19,9 @@ def __init__(self, filename): self.read_result(filename) + vel_diad = None # text files carry finite-difference velocities only + formulation = "legacy" + def read_result(self, filename): total_irq = 0 which_phonon = 0 @@ -390,6 +393,34 @@ def __init__(self, filename, channel="3ph", temperature=None, verbose=True): lo, hi = offsets[ik], offsets[ik + 1] self.vel[ik, :, : hi - lo, :] = np.transpose(vel[lo:hi], (1, 0, 2)) + # Velocity diad W[k,s,a,b] = sum_{s' in D(s)} Re(V^a_{ss'} V^b_{s's}), + # in (m/s)^2. Use W to rebuild kappa at degeneracies; its block sum is + # basis invariant. Absent in older files and legacy transport. + self.vel_diad = None + if "velocity_diad" in g: + diad = np.array(g["velocity_diad"][...], dtype=float) + if self.temperature_resolved: + if diad.ndim != 5: + raise RuntimeError( + "{}: unexpected velocity_diad layout for a temperature-resolved file".format( + filename + ) + ) + diad = diad[it] + diad = diad.reshape((knum.size, ns, 3, 3)) + self.vel_diad = np.zeros((nk, ns, nmax, 3, 3), dtype=float) + for ik in range(nk): + lo, hi = offsets[ik], offsets[ik + 1] + self.vel_diad[ik, :, : hi - lo, :, :] = np.transpose( + diad[lo:hi], (1, 0, 2, 3) + ) + + # Transport formulation that produced /kappa in this file ("legacy" when unstamped). + self.formulation = "legacy" + if "kappa" in f and "formulation" in f["kappa"].attrs: + a = f["kappa"].attrs["formulation"] + self.formulation = a.decode() if isinstance(a, bytes) else str(a) + gamma = np.array(g["gamma"][...], dtype=float) if gamma.ndim == 1: gamma = gamma[:, None] diff --git a/tools/analyzer/calculator.py b/tools/analyzer/calculator.py index 3ea250f7..66ddcace 100644 --- a/tools/analyzer/calculator.py +++ b/tools/analyzer/calculator.py @@ -61,6 +61,10 @@ def __init__( self.gamma4_interpolated = None self.gamma_iso = None # linediwth due to isotope scattering self.vel = None # Velocity array + self.vel_diad = ( + None # Degenerate-block velocity diad (None for legacy/text files) + ) + self.formulation = "legacy" self.vel4 = None self.qpoint_weight = None # Weight array self.qpoint_weight4 = None @@ -162,6 +166,15 @@ def set_variables_3ph(self): else: self.gamma3 = result.gamma self.vel = result.vel + self.vel_diad = getattr(result, "vel_diad", None) + self.formulation = getattr(result, "formulation", "legacy") + if self.vel_diad is None and self.formulation != "legacy": + warnings.warn( + "the file's kappa was assembled with the '{}' velocity formulation but it holds no " + "velocity_diad dataset, so kappa recomputed here (cumulative kappa etc.) follows the legacy " + "finite-difference formulation and will not match kappa_total; a restart with the current " + "anphon adds the dataset".format(self.formulation) + ) self.volume = result.volume self.qpoint_weight = result.multiplicity self.qpoints = result.q_coord @@ -197,10 +210,8 @@ def set_variables_4ph(self): if result.lattice_vector is not None: kinds = result.atomic_kinds if kinds is None: - # the text .result files store the fractional coordinates only; without the - # atomic kinds the symmetry is detected treating all atoms as one species, - # which can only over-count operations; the star-size check in the - # interpolator catches an inconsistent result. + # Without atomic kinds, .result symmetry treats all atoms as one species + # and may overcount operations. The interpolator checks star sizes. warnings.warn( "{} does not store the atomic kinds; detecting the symmetry with all " "atoms treated as one species (use the kappa.h5 file for exact kinds)".format( @@ -585,7 +596,7 @@ def get_group_velocity_norm(self): Returns: numpy.ndarray: Group velocity magnitudes with shape (nk, nmode). """ - return np.linalg.norm(self.vel[:, :, 0, :], axis=2) + return self._speed() def print_lifetime(self, temperature, four_phonon=False, isotope=False): """ @@ -672,7 +683,7 @@ def print_lifetime_mode( # Group velocity magnitude |v| (m/s) is temperature-independent; the mean # free path l = |v| * tau (nm) varies with temperature through tau. - vq = np.linalg.norm(self.vel[index_k, index_mode, 0, :]) + vq = self._speed()[index_k, index_mode] omega_q = self.omega[index_k, index_mode] print( @@ -730,6 +741,38 @@ def print_lifetime_mode( ) print("") + def _vv_per_copy(self): + """Velocity diad per symmetry copy, shape (nk, ns, nmax, 3, 3). + + From the stored degenerate-block diad when available (basis invariant once summed + over a block), otherwise the outer product of finite-difference velocities. + """ + if self.vel_diad is not None: + return self.vel_diad + return self.vel[:, :, :, :, None] * self.vel[:, :, :, None, :] + + def _vvprod(self): + """Diad summed over the symmetry copies of each irreducible k, shape (nk, ns, 3, 3).""" + return np.sum(self._vv_per_copy(), axis=2) + + def _speed(self): + """Per-mode speed |v| used for mean-free-paths, shape (nk, ns). + + sqrt(tr W) of the first copy when the diad is available; for a degenerate block this is + an effective speed whose block sum is invariant, whereas individual copies are not. + """ + if self.vel_diad is not None: + return np.sqrt( + np.maximum(np.einsum("ksaa->ks", self.vel_diad[:, :, 0, :, :]), 0.0) + ) + return np.linalg.norm(self.vel[:, :, 0, :], axis=2) + + def _speed_dir(self): + """Per-copy, per-direction speed |v_d|, shape (nk, ns, nmax, 3).""" + if self.vel_diad is not None: + return np.sqrt(np.maximum(np.einsum("kscaa->ksca", self.vel_diad), 0.0)) + return np.abs(self.vel) + def get_thermal_conductivity( self, four_phonon=False, isotope=False, len_boundary=None, gb_shape="sphere" ): @@ -750,14 +793,7 @@ def get_thermal_conductivity( nk, nmode = self.omega.shape - vvprod = np.zeros((nk, nmode, 3, 3), dtype=float) - for i in range(nk): - for j in range(nmode): - for k in range(3): - for m in range(3): - vvprod[i, j, k, m] = np.dot( - self.vel[i, j, :, k], self.vel[i, j, :, m] - ) + vvprod = self._vvprod() if len_boundary is None: for it, temp in enumerate(self.temperatures): @@ -774,7 +810,7 @@ def get_thermal_conductivity( else: if gb_shape == "sphere": assert len_boundary > 0.0, "The boundary length must be positive" - velnorm = np.linalg.norm(self.vel[:, :, 0, :], axis=2) + velnorm = self._speed() for it, temp in enumerate(self.temperatures): tau = self.get_lifetime( @@ -801,7 +837,8 @@ def get_thermal_conductivity( "3 elements when using 'cube' shape" ) - velnorm = np.abs(self.vel) + velnorm = self._speed_dir() + vv_copy = self._vv_per_copy() mfp = np.zeros_like(velnorm) @@ -821,8 +858,7 @@ def get_thermal_conductivity( cv * tau * np.sum( - self.vel[:, :, :, i] - * self.vel[:, :, :, j] + vv_copy[:, :, :, i, j] * len_boundary[i] / (len_boundary[i] + 2.0 * mfp[:, :, :, i]), axis=(2), @@ -909,15 +945,12 @@ def get_cumulative_kappa( nk, nmode = self.omega.shape - velnorm = np.linalg.norm(self.vel[:, :, 0, :], axis=2) + velnorm = self._speed() mfp = velnorm * tau * 0.001 - # Ignore numerically-zero mean-free-paths (e.g. Gamma acoustic modes) - # when determining the sampling range; 1e-6 nm matches the eps6 - # threshold used historically by the analyze_phonons C++ tool. - # modes whose linewidth was not computed (NaN) are excluded from the length grid - # and, through the comparisons below (False for NaN), from the sums + # Exclude near-zero MFPs from the grid (1e-6 nm matches analyze_phonons). + # NaN linewidths are excluded from both the grid and sums. if np.any(np.isnan(mfp)): warnings.warn( "{} modes have no linewidth (incomplete run) and are excluded from the " @@ -943,14 +976,7 @@ def get_cumulative_kappa( kappa = np.zeros((nsamples, 3, 3), dtype=float) if directions is None: - vvprod = np.zeros((nk, nmode, 3, 3), dtype=float) - for i in range(nk): - for j in range(nmode): - for k in range(3): - for m in range(3): - vvprod[i, j, k, m] = np.dot( - self.vel[i, j, :, k], self.vel[i, j, :, m] - ) + vvprod = self._vvprod() for ilen, len_boundary in enumerate(length_vec): tau_mod = np.where(mfp <= len_boundary, tau, 0.0) @@ -965,8 +991,9 @@ def get_cumulative_kappa( # The directional mean-free-path is evaluated per symmetry copy of each # irreducible k point (third axis of self.vel); zero-padded copies have # zero velocity and therefore never contribute. - mfp_dir = np.abs(self.vel) * tau[:, :, None, None] * 0.001 + mfp_dir = self._speed_dir() * tau[:, :, None, None] * 0.001 ctau = (cv * tau)[:, :, None] + vv_copy = self._vv_per_copy() for ilen, len_boundary in enumerate(length_vec): mask = np.ones(self.vel.shape[:3], dtype=bool) @@ -974,9 +1001,7 @@ def get_cumulative_kappa( mask &= mfp_dir[:, :, :, d] <= len_boundary for i in range(3): for j in range(3): - product = ( - ctau * self.vel[:, :, :, i] * self.vel[:, :, :, j] * mask - ) + product = ctau * vv_copy[:, :, :, i, j] * mask kappa[ilen, i, j] = np.sum(product) factor_toSI = ( diff --git a/tools/dfc2.py b/tools/dfc2.py index 4636900a..8dd3ef09 100644 --- a/tools/dfc2.py +++ b/tools/dfc2.py @@ -210,12 +210,9 @@ def update( prim_kind = fc2_data.atomic_kinds inv_prim = np.linalg.inv(prim_lat) # cart @ inv_prim -> fractional - # --- 1. The dfc2 cell should be the same cell as the HDF5 PrimitiveCell - # (ANPHON prints its own primitive cell as the dfc2 header). If - # they differ only in size/shape -- e.g. corrections fitted at a - # different volume -- the atomic correspondence can still be - # established from fractional coordinates, but the transferred - # force constants are only an approximation. + # The dfc2 cell should match HDF5 PrimitiveCell. Fractional coordinates + # can map atoms between different cell shapes/volumes, but transferring + # those force constants is approximate. lattice_match = np.allclose(dfc2_correction.lattice, prim_lat, atol=1.0e-3) if not lattice_match: report = cls._cell_mismatch_report(dfc2_correction.lattice, prim_lat) @@ -231,10 +228,8 @@ def update( print("WARNING: applying corrections across mismatched primitive cells.") print(report) - # Map each dfc2 atom onto a PrimitiveCell atom index by nearest position - # modulo a lattice translation (an identity map when the cells agree). - # Require a one-to-one correspondence so a deformed/mismatched cell that - # cannot be aligned is rejected rather than silently mis-mapped. + # Map dfc2 atoms one-to-one to PrimitiveCell atoms by nearest position + # modulo lattice translations; reject cells that cannot be aligned. dfc2_to_prim = np.full(len(dfc2_correction.positions), -1, dtype=np.int64) used = {} max_frac_res = 0.0 @@ -350,15 +345,9 @@ def update( if abs(val) > tol: n_nonzero_applied += 1 - # --- 5. Diagnostics. - # A physical correction may be listed several times in the dfc2 - # file (once per translational copy of an atom in the conventional - # cell). An unconsumed nonzero row is an expected duplicate only - # if an *equivalent* correction was actually consumed. Equivalence - # is tested with a translation-invariant signature -- the Cartesian - # bond vector together with the two coordinate directions uniquely - # identifies the force-constant component -- so a genuine miss can - # never be silently reclassified as a duplicate. + # Treat an unused nonzero row as a duplicate only if an equivalent row + # was consumed. Match by Cartesian bond vector and coordinate directions + # to distinguish translational copies from missing corrections. def signature(n): a0 = int(dfc2_to_prim[dfc2_correction.atoms[n, 0]]) a1 = int(dfc2_to_prim[dfc2_correction.atoms[n, 1]]) diff --git a/tools/interface/VASP.py b/tools/interface/VASP.py index 7c579d5d..12745ccc 100644 --- a/tools/interface/VASP.py +++ b/tools/interface/VASP.py @@ -631,9 +631,7 @@ def parse_or_repair_xml_file(file_path): except etree.ParseError as e: print(f"Failed to parse XML file with ElementTree due to: {e}") print("Consider installing lxml for better XML parsing support.") - # Handle the error or attempt a manual repair if necessary - # Note: ElementTree doesn't have a built-in 'recover' mode like lxml, - # so you may need to manually fix the XML or use a different strategy + # ElementTree cannot recover malformed XML. def _get_coordinates_and_forces(self, file_to_parse): hdf5_mode = file_to_parse.lower().split(".")[-1] in ["h5", "hdf5"] diff --git a/tools/strainkit/almfit.py b/tools/strainkit/almfit.py index 93a5e4a5..e0b79801 100644 --- a/tools/strainkit/almfit.py +++ b/tools/strainkit/almfit.py @@ -40,9 +40,8 @@ def model_kwargs(nbody, cutoff, nkd): def make_alm(atoms, verbosity=0): ALM = require_alm() - # The alm Python API takes the lattice vectors as rows (row i = i-th - # vector), i.e. the ase convention; the transposition to the column - # convention of the C++ core happens inside the wrapper. + # The Python API uses lattice-vector rows (ASE convention); + # the wrapper transposes them for the C++ core. return ALM( np.asarray(atoms.cell[:], dtype=float), np.asarray(atoms.get_scaled_positions(wrap=False), dtype=float), @@ -125,9 +124,8 @@ def fit_harmonic( raise ValueError(f"training data must have shape (nsnap, {len(atoms)}, 3)") with make_alm(atoms, verbosity) as alm: if transmat_to_prim is not None: - # ``atoms`` is already the (strained) supercell, so the cell given - # to ALM is the supercell itself (SUPERCELL = identity) and only - # PRIMCELL has to be declared. Must precede define(). + # atoms is already the supercell: use SUPERCELL = identity and set + # PRIMCELL before define(). alm.set_supercell(np.eye(3), transmat_to_prim) _define_harmonic(alm, atoms, nbody, cutoff) alm.set_constraint(translation=True) diff --git a/tools/strainkit/elasticfit.py b/tools/strainkit/elasticfit.py index 3ac5c101..c699d054 100644 --- a/tools/strainkit/elasticfit.py +++ b/tools/strainkit/elasticfit.py @@ -371,11 +371,8 @@ def fit_elastic(data, volume, mode="stress", e_ref=None, weight_floor=1.0e-8): wrow = np.ones(len(b)) block_rms = {} if mode == "both": - # One reweighting step: fit all rows unweighted first, take the residual - # RMS of every block from that common solution, and weight the rows by - # its inverse (floored). This is symmetric in the two blocks and does - # not depend on the single-block fits being full rank (an energy-only - # block is rank deficient on the minimal direction set). + # Fit all rows unweighted, then reweight each block by its floored inverse + # residual RMS. This also works when individual blocks are rank deficient. x0, _, _, _ = _lstsq(A, b) res0 = A @ x0 - b for kind in ("energy", "stress"): diff --git a/tools/taylor.py b/tools/taylor.py index 23ab120b..050d7f68 100644 --- a/tools/taylor.py +++ b/tools/taylor.py @@ -85,10 +85,8 @@ def set_forceconstants(self, fcs_dic, primitive_cell_alm): gamma_values = self.gamma(flatten_indices) gamma_scaled_fcs = gamma_values * self.fcs_values[fckey] - # Sort the entries by (first atom, first coordinate) so that the - # force contributions sharing the same target component become - # contiguous segments. The grouping is independent of the - # translation because each translation permutes the atoms. + # Group entries by (first atom, first coordinate) for contiguous force sums. + # Translations permute atoms but preserve this grouping. sort_keys = 3 * atom_indices_taylor[:, 0] + coord_indices_taylor[:, 0] sort_order = np.argsort(sort_keys, kind="stable") keys_sorted = sort_keys[sort_order] @@ -156,9 +154,7 @@ def compute(self, displacements): "ij,ij->i", ff_tmp, displacements_flat[:, flat_indices[:, 0]] ) - # The entries are pre-sorted by (first atom, first coordinate), - # so the contributions to each force component form contiguous - # segments and the scatter-add stays buffered. + # Sum contiguous force-component segments with buffered scatter-add. target_indices = ( 3 * self.map_translation[group_atoms, itran] + group_coords )