anphon: convert kpoint/system/fcs_phonon/write_phonons arrays to NDArray (P8)

Final phase of the allocate()/deallocate() retirement. The three
k-point owner classes (KpointGeneral, KpointMeshUniform,
KpointBandStructure) hold xk/kvec_na/kaxis via NDArray, which makes
their tables deep-const through the const KpointMeshUniform* borrows
everywhere; the read-only view parameters that were sloppily declared
mutable become const along the way (exec_interpolation's
xk_dense/kvec_dense, duplicate_xk_boundary, compute_dV2_dumn and
compute_d2V2_dumn2, and the non-template rotvec overloads' vec_in),
and pointer-copy aliases of the k tables become reference bindings.
Fcs_phonon::force_constant_with_cell converts too, surfacing three
methods that mutated it through a shallow const
(replicate_force_constants, load_fcs_from_file, MPI_Bcast_fcs_array) -
they lose the qualifier, and examine_translational_invariance's switch
cases gain scopes for the now non-trivially-constructed accumulators.
System, symmetry_core, and write_phonons locals convert wholesale.

The only allocate() calls left in anphon are the spglib C-ABI
boundaries (double(*)[3] / int(*)[3][3] in symmetry_core and
relaxation), each carrying a comment saying why; no raw new[] remains.

Full 8-suite regression battery + test_units pass, alm still builds
against the shared mathfunctions.h, and the ndarray unit test passes.
A/B timing vs the pre-migration binary: Si RTA -0.14%, BaTiO3 SCPH
+0.32%.

Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>

Author: ttadano

anphon/anharmonic_core.cpp

@@ -912,7 +912,7 @@ void AnharmonicCore::calc_damping_tetrahedron(const unsigned int ntemp, const do
 
     const auto knum = kmesh_in->kpoint_irred_all[ik_in][0].knum;
     const auto knum_minus = kmesh_in->kindex_minus_xk[knum];
-    const auto xk = kmesh_in->xk;
+    const auto &xk = kmesh_in->xk;
 
     kmap_identity.resize(nk);
 

anphon/dynamical.cpp

@@ -1982,7 +1982,7 @@ void Dynamical::replicate_dymat_for_all_kpoints(const KpointMeshUniform *kmesh_c
 
 
 void Dynamical::exec_interpolation(const unsigned int kmesh_orig[3], std::complex<double> ***dymat_r,
-                                   const unsigned int nk_dense, double **xk_dense, double **kvec_dense,
+                                   const unsigned int nk_dense, const double *const *xk_dense, const double *const *kvec_dense,
                                    double **eval_out, std::complex<double> ***evec_out,
                                    const std::vector<Eigen::MatrixXcd> &dymat_short,
                                    const std::vector<Eigen::MatrixXcd> &dymat_long, MinimumDistList ***mindist_list_in,
@@ -2298,7 +2298,7 @@ void Dynamical::calc_new_dymat_with_evec(std::complex<double> ***dymat_out, doub
     }
 }
 
-void Dynamical::duplicate_xk_boundary(double *xk_in, std::vector<std::vector<double>> &vec_xk)
+void Dynamical::duplicate_xk_boundary(const double *xk_in, std::vector<std::vector<double>> &vec_xk)
 {
     int i;
     int n[3];

anphon/dynamical.h

@@ -175,7 +175,7 @@ class Dynamical: protected Pointers
     const double *const *get_xrs_image() const;
 
     void exec_interpolation(const unsigned int kmesh_orig[3], std::complex<double> ***dymat_r,
-                            const unsigned int nk_dense, double **xk_dense, double **kvec_dense, double **eval_out,
+                            const unsigned int nk_dense, const double *const *xk_dense, const double *const *kvec_dense, double **eval_out,
                             std::complex<double> ***evec_out, const std::vector<Eigen::MatrixXcd> &dymat_short,
                             const std::vector<Eigen::MatrixXcd> &dymat_long, MinimumDistList ***mindist_list_in,
                             const bool use_precomputed_dymat = false, const bool return_sqrt = true) const;
@@ -220,7 +220,7 @@ class Dynamical: protected Pointers
                                Eigen::MatrixXcd &evec_sub) const;
 
 
-    void duplicate_xk_boundary(double *, std::vector<std::vector<double>> &);
+    void duplicate_xk_boundary(const double *, std::vector<std::vector<double>> &);
 
 
     NDArray<double, 2> xshift_s;

anphon/fcs_phonon.cpp

@@ -52,14 +52,10 @@ void Fcs_phonon::set_default_variables()
     file_fc4 = "";
 
     update_fc2 = false;
-    force_constant_with_cell = nullptr;
 }
 
 void Fcs_phonon::deallocate_variables()
 {
-    if (force_constant_with_cell) {
-        deallocate(force_constant_with_cell);
-    }
 }
 
 void Fcs_phonon::setup(const std::string &mode)
@@ -111,7 +107,7 @@ void Fcs_phonon::setup(const std::string &mode)
         // parser (parse_analysis_vars).
     }
 
-    allocate(force_constant_with_cell, maxorder);
+    force_constant_with_cell.resize(maxorder);
 
     if (mympi->my_rank == 0) {
 
@@ -138,7 +134,7 @@ void Fcs_phonon::setup(const std::string &mode)
     replicate_force_constants(maxorder);
 }
 
-void Fcs_phonon::replicate_force_constants(const int maxorder_in) const
+void Fcs_phonon::replicate_force_constants(const int maxorder_in)
 {
     for (auto order = 0; order < maxorder_in; ++order) {
         replicate_force_constant(system.get(), force_constant_with_cell[order]);
@@ -245,7 +241,7 @@ void Fcs_phonon::replicate_force_constant(const System *system_in, std::vector<F
 }
 
 
-void Fcs_phonon::load_fcs_from_file(const int maxorder_in) const
+void Fcs_phonon::load_fcs_from_file(const int maxorder_in)
 {
     std::vector filename_list{file_fc2, file_fc3, file_fc4};
 
@@ -725,9 +721,9 @@ double Fcs_phonon::examine_translational_invariance(const int order, const unsig
     const auto natmin3 = 3 * natmin;
 
     switch (order) {
-    case 0:
-        double **sum2;
-        allocate(sum2, natmin3, 3);
+    case 0: {
+        NDArray<double, 2> sum2;
+        sum2.resize(natmin3, 3);
 
         for (j = 0; j < natmin3; ++j) {
             for (k = 0; k < 3; ++k) {
@@ -747,11 +743,12 @@ double Fcs_phonon::examine_translational_invariance(const int order, const unsig
                 ret = std::max(ret, dev);
             }
         }
-        deallocate(sum2);
+        sum2.clear();
         break;
-    case 1:
-        double ***sum3;
-        allocate(sum3, 3 * natmin, 3 * nat, 3);
+    }
+    case 1: {
+        NDArray<double, 3> sum3;
+        sum3.resize(3 * natmin, 3 * nat, 3);
 
         for (j = 0; j < natmin3; ++j) {
             for (k = 0; k < nat3; ++k) {
@@ -775,11 +772,12 @@ double Fcs_phonon::examine_translational_invariance(const int order, const unsig
                 }
             }
         }
-        deallocate(sum3);
+        sum3.clear();
         break;
-    case 2:
-        double ****sum4;
-        allocate(sum4, natmin3, nat3, nat3, 3);
+    }
+    case 2: {
+        NDArray<double, 4> sum4;
+        sum4.resize(natmin3, nat3, nat3, 3);
 
         for (j = 0; j < natmin3; ++j) {
             for (k = 0; k < nat3; ++k) {
@@ -809,21 +807,22 @@ double Fcs_phonon::examine_translational_invariance(const int order, const unsig
                 }
             }
         }
-        deallocate(sum4);
+        sum4.clear();
         break;
+    }
     default:
         break;
     }
 
     return ret;
 }
 
-void Fcs_phonon::MPI_Bcast_fcs_array(const unsigned int N) const
+void Fcs_phonon::MPI_Bcast_fcs_array(const unsigned int N)
 {
     int j, k;
-    double *fcs_tmp;
-    unsigned int ***ind;
-    double ***relative_vector_tmp;
+    NDArray<double, 1> fcs_tmp;
+    NDArray<unsigned int, 3> ind;
+    NDArray<double, 3> relative_vector_tmp;
 
     std::vector<AtomCellSuper> ivec_array;
     std::vector<unsigned int> atoms_s_tmp;
@@ -840,9 +839,9 @@ void Fcs_phonon::MPI_Bcast_fcs_array(const unsigned int N) const
 
         if (len == 0) continue;
 
-        allocate(fcs_tmp, len);
-        allocate(ind, len, nelem, 4);
-        allocate(relative_vector_tmp, len, nelem - 1, 3);
+        fcs_tmp.resize(len);
+        ind.resize(len, nelem, 4);
+        relative_vector_tmp.resize(len, nelem - 1, 3);
 
         if (mympi->my_rank == 0) {
             for (j = 0; j < len; ++j) {
@@ -892,8 +891,8 @@ void Fcs_phonon::MPI_Bcast_fcs_array(const unsigned int N) const
             }
         }
 
-        deallocate(fcs_tmp);
-        deallocate(ind);
-        deallocate(relative_vector_tmp);
+        fcs_tmp.clear();
+        ind.clear();
+        relative_vector_tmp.clear();
     }
 }

anphon/fcs_phonon.h

@@ -13,6 +13,7 @@
 #include <set>
 #include <string>
 #include <vector>
+#include "ndarray.h"
 #include "mathfunctions.h"
 #include "pointers.h"
 
@@ -193,7 +194,7 @@ class Fcs_phonon: protected Pointers
     unsigned int maxorder;
     std::string file_fcs, file_fc2, file_fc3, file_fc4;
 
-    std::vector<FcsArrayWithCell> *force_constant_with_cell;
+    NDArray<std::vector<FcsArrayWithCell>, 1> force_constant_with_cell;
 
     bool update_fc2;
 
@@ -229,17 +230,17 @@ class Fcs_phonon: protected Pointers
 
     void append_delta_fc2_from_scph(const std::string &fname_dfc2, std::vector<FcsArrayWithCell> &fcs_out) const;
 
-    void load_fcs_from_file(const int maxorder_in) const;
+    void load_fcs_from_file(const int maxorder_in);
 
 
     static double examine_translational_invariance(int order, unsigned int nat, unsigned int natmin,
                                                    const std::vector<FcsArrayWithCell> &fc_in);
 
-    void replicate_force_constants(const int maxorder_in) const;
+    void replicate_force_constants(const int maxorder_in);
 
     void MPI_Bcast_fc_class(unsigned int) const;
 
-    void MPI_Bcast_fcs_array(unsigned int) const;
+    void MPI_Bcast_fcs_array(unsigned int);
 
     void MPI_Bcast_fc2_ext();
 };

anphon/gruneisen.cpp

@@ -132,7 +132,7 @@ void Gruneisen::calc_gruneisen()
 
     if (kpoint->kpoint_bs.get()) {
         const auto nk = kpoint->kpoint_bs->nk;
-        const auto xk = kpoint->kpoint_bs->xk;
+        const auto &xk = kpoint->kpoint_bs->xk;
         const auto eval = dynamical->dymat_band->get_eigenvalues();
         const auto evec = dynamical->dymat_band->get_eigenvectors();
 
@@ -165,7 +165,7 @@ void Gruneisen::calc_gruneisen()
 
     if (dos->kmesh_dos.get()) {
         const auto nk = dos->kmesh_dos->nk;
-        const auto xk = dos->kmesh_dos->xk;
+        const auto &xk = dos->kmesh_dos->xk;
         const auto eval = dos->dymat_dos->get_eigenvalues();
         const auto evec = dos->dymat_dos->get_eigenvectors();
 

anphon/ifc_derivative.cpp

@@ -214,7 +214,7 @@ void DerivativeIFC::compute_d3V1_dumn3(MatrixXcdRowMajor &del3_v1_del_umn3,
 
 void DerivativeIFC::compute_dV2_dumn(std::vector<MatrixXcdRowMajor> &del_v2_del_umn,
                                      const std::complex<double> *const *const *const evec_harmonic,
-                                     const unsigned int nk, double **xk_in) const
+                                     const unsigned int nk, const double *const *xk_in) const
 {
     using namespace Eigen;
 
@@ -270,7 +270,7 @@ void DerivativeIFC::compute_dV2_dumn(std::vector<MatrixXcdRowMajor> &del_v2_del_
 
 void DerivativeIFC::compute_d2V2_dumn2(std::vector<MatrixXcdRowMajor> &del2_v2_del_umn2,
                                        const std::complex<double> *const *const *const evec_harmonic,
-                                       const unsigned int nk, double **xk_in) const
+                                       const unsigned int nk, const double *const *xk_in) const
 {
     using namespace Eigen;
 

anphon/ifc_derivative.h

@@ -54,11 +54,11 @@ class DerivativeIFC
 
     void compute_dV2_dumn(std::vector<MatrixXcdRowMajor> &dV2_dumn,
                           const std::complex<double> *const *const *const evec_harmonic, unsigned int nk,
-                          double **xk_in) const;
+                          const double *const *xk_in) const;
 
     void compute_d2V2_dumn2(std::vector<MatrixXcdRowMajor> &d2V2_dumn2,
                             const std::complex<double> *const *const *const evec_harmonic, unsigned int nk,
-                            double **xk_in) const;
+                            const double *const *xk_in) const;
 
     void compute_dV3_dumn(std::vector<std::vector<MatrixXcdRowMajor>> &dV3_dumn, double **omega2_harmonic,
                           const std::complex<double> *const *const *const evec_harmonic,