Velocity matrix transport

alm/accelerate_solver.cpp

@@ -1,8 +1,4 @@
-// accelerate_solver.cpp
-//
-// Translation unit isolating <Eigen/AccelerateSupport> (and thus <Accelerate/Accelerate.h>) 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

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 <Eigen/AccelerateSupport> pulls in <Accelerate/Accelerate.h>, 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

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;
 }
 

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}.

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();

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<IntList> in fcs.cpp
-        // and ensures consistent behavior across platforms
-        // std::cout << "== operator called\n";
+        // Vector equality supports std::unordered_set<IntList> in fcs.cpp.
         return iarray == a.iarray;
     }
 };

alm/constraint.cpp

@@ -344,16 +344,9 @@ auto Constraint::update_constraint_matrix(const std::unique_ptr<System> &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> &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> &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> &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<System> &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);

alm/fcs.cpp

@@ -313,10 +313,8 @@ auto Fcs::get_constraint_symmetry(const size_t nat, const std::unique_ptr<Symmet
                                   const size_t nparams, const double tolerance, ConstraintSparseForm &const_out,
                                   const ReductionAlgo algo_in) -> 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;

alm/input_parser.cpp

@@ -178,10 +178,8 @@ auto InputParser::parse_displacement_and_force_files(std::vector<std::vector<dou
 
 auto InputParser::parse_energies(std::vector<double> &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<int>(transmat_to_super.determinant());
     const auto ntoken_per_snapshot = static_cast<size_t>(6 * nat);
 
@@ -260,14 +258,7 @@ auto InputParser::parse_energies(std::vector<double> &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<std::string> &kdname_vec_out,
                                          std::vector<int> &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);

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();