diff --git a/src/ML-POD/fitpod_command.cpp b/src/ML-POD/fitpod_command.cpp
index 8e18f67c28..619f82f4a8 100644
--- a/src/ML-POD/fitpod_command.cpp
+++ b/src/ML-POD/fitpod_command.cpp
@@ -69,7 +69,7 @@ void FitPOD::command(int narg, char **arg)
 
   // get POD coefficients from an input file
 
-  if (coeff_file != "") podptr->podArrayCopy(desc.c, podptr->pod.coeff, podptr->pod.nd);
+  if (coeff_file != "") podArrayCopy(desc.c, podptr->pod.coeff, podptr->pod.nd);
 
   // compute POD coefficients using least-squares method
 
@@ -507,13 +507,13 @@ void FitPOD::get_data(datastruct &data, std::vector<std::string> species)
   }
 
   int len = data.num_atom.size();
-  data.num_atom_min = podptr->podArrayMin(&data.num_atom[0], len);
-  data.num_atom_max = podptr->podArrayMax(&data.num_atom[0], len);
+  data.num_atom_min = podArrayMin(&data.num_atom[0], len);
+  data.num_atom_max = podArrayMax(&data.num_atom[0], len);
   data.num_atom_cumsum.resize(len+1);
-  podptr->podCumsum(&data.num_atom_cumsum[0], &data.num_atom[0], len+1);
+  podCumsum(&data.num_atom_cumsum[0], &data.num_atom[0], len+1);
 
   data.num_config_cumsum.resize(nfiles+1);
-  podptr->podCumsum(&data.num_config_cumsum[0], &data.num_config[0], nfiles+1);
+  podCumsum(&data.num_config_cumsum[0], &data.num_config[0], nfiles+1);
 
   // convert all structures to triclinic system
 
@@ -529,11 +529,11 @@ void FitPOD::get_data(datastruct &data, std::vector<std::string> species)
     double *a2 = &lattice[3];
     double *a3 = &lattice[6];
 
-    podptr->matrix33_inverse(Qmat, a1, a2, a3);
-    podptr->triclinic_lattice_conversion(a1, a2, a3, a1, a2, a3);
-    podptr->matrix33_multiplication(Qmat, lattice, Qmat, DIM);
-    podptr->matrix33_multiplication(x, Qmat, x, natom);
-    podptr->matrix33_multiplication(f, Qmat, f, natom);
+    matrix33_inverse(Qmat, a1, a2, a3);
+    triclinic_lattice_conversion(a1, a2, a3, a1, a2, a3);
+    matrix33_multiplication(Qmat, lattice, Qmat, DIM);
+    matrix33_multiplication(x, Qmat, x, natom);
+    matrix33_multiplication(f, Qmat, f, natom);
   }
 
   if (comm->me == 0) {
@@ -640,12 +640,12 @@ void FitPOD::select_data(datastruct &newdata, datastruct data)
     newdata.num_atom_each_file[file] = num_atom_sum;
   }
   int len = newdata.num_atom.size();
-  newdata.num_atom_min = podptr->podArrayMin(&newdata.num_atom[0], len);
-  newdata.num_atom_max = podptr->podArrayMax(&newdata.num_atom[0], len);
+  newdata.num_atom_min = podArrayMin(&newdata.num_atom[0], len);
+  newdata.num_atom_max = podArrayMax(&newdata.num_atom[0], len);
   newdata.num_atom_cumsum.resize(len+1);
-  podptr->podCumsum(&newdata.num_atom_cumsum[0], &newdata.num_atom[0], len+1);
+  podCumsum(&newdata.num_atom_cumsum[0], &newdata.num_atom[0], len+1);
   newdata.num_atom_sum = newdata.num_atom_cumsum[len];
-  podptr->podCumsum(&newdata.num_config_cumsum[0], &newdata.num_config[0], nfiles+1);
+  podCumsum(&newdata.num_config_cumsum[0], &newdata.num_config[0], nfiles+1);
   newdata.num_config_sum = newdata.num_atom.size();
 
   int n = data.num_config_sum;
@@ -867,7 +867,7 @@ int FitPOD::podfullneighborlist(double *y, int *alist, int *neighlist, int *numn
 
   nn = podneighborlist(neighlist, numneigh, y, rcutsq, nx, N, dim);
 
-  podptr->podCumsum(numneighsum, numneigh, nx+1);
+  podCumsum(numneighsum, numneigh, nx+1);
 
   return nn;
 }
@@ -879,9 +879,9 @@ void FitPOD::allocate_memory(datastruct data)
   memory->create(desc.A, nd*nd, "fitpod:desc_A");
   memory->create(desc.b, nd, "fitpod:desc_b");
   memory->create(desc.c, nd, "fitpod:desc_c");
-  podptr->podArraySetValue(desc.A, 0.0, nd*nd);
-  podptr->podArraySetValue(desc.b, 0.0, nd);
-  podptr->podArraySetValue(desc.c, 0.0, nd);
+  podArraySetValue(desc.A, 0.0, nd*nd);
+  podArraySetValue(desc.b, 0.0, nd);
+  podArraySetValue(desc.c, 0.0, nd);
 
   int dim = 3;
   int natom_max = data.num_atom_max;
@@ -1192,7 +1192,7 @@ void FitPOD::least_squares_matrix(datastruct data, int ci)
 
   // least-square matrix for all descriptors: A = A + (we*we)*(gd^T * gd)
 
-  podptr->podKron(desc.A, desc.gd, desc.gd, we2, nd, nd);
+  podKron(desc.A, desc.gd, desc.gd, we2, nd, nd);
 
   // least-square matrix for all descriptors derivatives: A =  A + (wf*wf) * (gdd^T * gdd)
 
@@ -1340,7 +1340,7 @@ double FitPOD::energyforce_calculation(double *force, double *coeff, datastruct
   podptr->podNeighPairs(rij, nb.y, idxi, ai, aj, ti, tj, nb.pairnum_cumsum, atomtype, nb.pairlist, nb.alist, natom);
 
   double *effectivecoeff = &tmpmem[3*Nij]; // 3*Nij
-  podptr->podArraySetValue(effectivecoeff, 0.0, nd1234);
+  podArraySetValue(effectivecoeff, 0.0, nd1234);
 
   double energy = podptr->energyforce_calculation(force, coeff, effectivecoeff, desc.gd, rij,
     &tmpmem[3*Nij+nd1234], nb.pairnum_cumsum, atomtype, idxi, ai, aj, ti, tj, natom, Nij);
@@ -1496,8 +1496,8 @@ void FitPOD::error_analysis(datastruct data, double *coeff)
         outarray[2 + m*ci] = energy;
         outarray[3 + m*ci] = DFTenergy;
         outarray[4 + m*ci] = fabs(DFTenergy-energy)/natom;
-        outarray[5 + m*ci] = podptr->podArrayNorm(force.data(), nforce);
-        outarray[6 + m*ci] = podptr->podArrayNorm(DFTforce, nforce);
+        outarray[5 + m*ci] = podArrayNorm(force.data(), nforce);
+        outarray[6 + m*ci] = podArrayNorm(DFTforce, nforce);
 
         double diff, sum = 0.0, ssr = 0.0;
         for (int j=0; j<dim*natom; j++) {
@@ -1604,3 +1604,216 @@ void FitPOD::energyforce_calculation(datastruct data, double *coeff)
   if (comm->me == 0)
     utils::logmesg(lmp, "**************** End of Energy/Force Calculation ****************\n");
 }
+
+void FitPOD::print_matrix(const char *desc, int m, int n, double **a, int /*lda*/ )
+{
+  int i, j;
+  printf( "\n %s\n", desc );
+
+  for( i = 0; i < m; i++ )
+  {
+    for( j = 0; j < n; j++ ) printf( " %6.12f", a[j][i] );
+    printf( "\n" );
+  }
+}
+
+void FitPOD::print_matrix(const char *desc, int m, int n, double *a, int lda )
+{
+  int i, j;
+  printf( "\n %s\n", desc );
+
+  for( i = 0; i < m; i++ )
+  {
+    for( j = 0; j < n; j++ ) printf( " %6.12f", a[i+j*lda] );
+    printf( "\n" );
+  }
+}
+
+void FitPOD::print_matrix(const char *desc, int m, int n, int *a, int lda)
+{
+  int i, j;
+  printf( "\n %s\n", desc );
+
+  for( i = 0; i < m; i++ )
+  {
+    for( j = 0; j < n; j++ ) printf( " %d", a[i+j*lda] );
+    printf( "\n" );
+  }
+}
+
+void FitPOD::podArrayFill(int* output, int start, int length)
+{
+        for (int j = 0; j < length; ++j)
+                output[j] = start + j;
+}
+
+double FitPOD::podArraySum(double *a, int n)
+{
+  double e = a[0];
+  for (int i=1; i<n; i++)
+    e += a[i];
+  return e;
+}
+
+double FitPOD::podArrayMin(double *a, int n)
+{
+  double b = a[0];
+  for (int i=1; i<n; i++)
+    if (a[i]<b)
+      b = a[i];
+  return b;
+}
+
+double FitPOD::podArrayMax(double *a, int n)
+{
+  double b = a[0];
+  for (int i=1; i<n; i++)
+    if (a[i]>b)
+      b = a[i];
+  return b;
+}
+
+int FitPOD::podArrayMin(int *a, int n)
+{
+  int b = a[0];
+  for (int i=1; i<n; i++)
+    if (a[i]<b)
+      b = a[i];
+  return b;
+}
+
+int FitPOD::podArrayMax(int *a, int n)
+{
+  int b = a[0];
+  for (int i=1; i<n; i++)
+    if (a[i]>b)
+      b = a[i];
+  return b;
+}
+
+void FitPOD::podKron(double *C, double *A, double *B, double alpha, int M1, int M2)
+{
+  int M = M1*M2;
+  for (int idx=0; idx<M; idx++)
+  {
+    int ib = idx%M2;
+    int ia = (idx-ib)/M2;
+    C[idx] += alpha*A[ia]*B[ib];
+  }
+}
+
+void FitPOD::podCumsum(int* output, int* input, int length)
+{
+        output[0] = 0;
+        for (int j = 1; j < length; ++j)
+                output[j] = input[j - 1] + output[j - 1];
+}
+
+double FitPOD::podArrayNorm(double *a, int n)
+{
+  double e = a[0]*a[0];
+  for (int i=1; i<n; i++)
+    e += a[i]*a[i];
+  return sqrt(e);
+}
+
+double FitPOD::podArrayErrorNorm(double *a, double *b, int n)
+{
+  double e = (a[0]-b[0])*(a[0]-b[0]);
+  for (int i=1; i<n; i++)
+    e += (a[i]-b[i])*(a[i]-b[i]);
+  return sqrt(e);
+}
+
+void FitPOD::podArraySetValue(double *y, double a, int n)
+{
+  for (int i=0; i<n; i++)
+    y[i] = a;
+}
+
+void FitPOD::podArrayCopy(double *y, double *x, int n)
+{
+  for (int i=0; i<n; i++)
+    y[i] = x[i];
+}
+
+void FitPOD::rotation_matrix(double *Rmat, double alpha, double beta, double gamma)
+{
+  double ca = cos(alpha);
+  double cb = cos(beta);
+  double cg = cos(gamma);
+  double sa = sin(alpha);
+  double sb = sin(beta);
+  double sg = sin(gamma);
+
+  Rmat[0] = ca*cg*cb - sa*sg;
+  Rmat[3] = -ca*cb*sg - sa*cg;
+  Rmat[6] = ca*sb;
+
+  Rmat[1] = sa*cg*cb + ca*sg;
+  Rmat[4] = -sa*cb*sg + ca*cg;
+  Rmat[7] = sa*sb;
+
+  Rmat[2] = -sb*cg;
+  Rmat[5] = sb*sg;
+  Rmat[8] = cb;
+}
+
+void FitPOD::matrix33_multiplication(double *xrot, double *Rmat, double *x, int natom)
+{
+  double x1, x2, x3;
+  for (int i=0; i < natom; i++) {
+    x1 = x[0 + 3*i];
+    x2 = x[1 + 3*i];
+    x3 = x[2 + 3*i];
+    xrot[0 + 3*i] = Rmat[0]*x1 + Rmat[3]*x2 + Rmat[6]*x3;
+    xrot[1 + 3*i] = Rmat[1]*x1 + Rmat[4]*x2 + Rmat[7]*x3;
+    xrot[2 + 3*i] = Rmat[2]*x1 + Rmat[5]*x2 + Rmat[8]*x3;
+  }
+}
+
+void FitPOD::matrix33_inverse(double *invA, double *A1, double *A2, double *A3)
+{
+  double a11 = A1[0];
+  double a21 = A1[1];
+  double a31 = A1[2];
+  double a12 = A2[0];
+  double a22 = A2[1];
+  double a32 = A2[2];
+  double a13 = A3[0];
+  double a23 = A3[1];
+  double a33 = A3[2];
+  double detA = (a11*a22*a33 - a11*a23*a32 - a12*a21*a33 + a12*a23*a31 + a13*a21*a32 - a13*a22*a31);
+
+  invA[0] = (a22*a33 - a23*a32)/detA;
+  invA[1] = (a23*a31 - a21*a33)/detA;
+  invA[2] = (a21*a32 - a22*a31)/detA;
+  invA[3] = (a13*a32 - a12*a33)/detA;
+  invA[4] = (a11*a33 - a13*a31)/detA;
+  invA[5] = (a12*a31 - a11*a32)/detA;
+  invA[6] = (a12*a23 - a13*a22)/detA;
+  invA[7] = (a13*a21 - a11*a23)/detA;
+  invA[8] = (a11*a22 - a12*a21)/detA;
+}
+
+void FitPOD::triclinic_lattice_conversion(double *a, double *b, double *c, double *A, double *B, double *C)
+{
+  double Anorm = sqrt(A[0]*A[0] + A[1]*A[1] + A[2]*A[2]);
+  double Bnorm = sqrt(B[0]*B[0] + B[1]*B[1] + B[2]*B[2]);
+  double Cnorm = sqrt(C[0]*C[0] + C[1]*C[1] + C[2]*C[2]);
+
+  double Ahat[3];
+  Ahat[0] = A[0]/Anorm; Ahat[1] = A[1]/Anorm; Ahat[2] = A[2]/Anorm;
+
+  double ax = Anorm;
+  double bx = B[0]*Ahat[0] + B[1]*Ahat[1] + B[2]*Ahat[2]; //dot(B,Ahat);
+  double by = sqrt(Bnorm*Bnorm - bx*bx); //sqrt(Bnorm^2 - bx^2);// #norm(cross(Ahat,B));
+  double cx = C[0]*Ahat[0] + C[1]*Ahat[1] + C[2]*Ahat[2]; // dot(C,Ahat);
+  double cy = (B[0]*C[0] + B[1]*C[1] + B[2]*C[2] - bx*cx)/by; // (dot(B, C) - bx*cx)/by;
+  double cz = sqrt(Cnorm*Cnorm - cx*cx - cy*cy); // sqrt(Cnorm^2 - cx^2 - cy^2);
+
+  a[0] = ax; a[1] = 0.0; a[2] = 0.0;
+  b[0] = bx; b[1] = by;  b[2] = 0.0;
+  c[0] = cx; c[1] = cy;  c[2] = cz;
+}
+
diff --git a/src/ML-POD/fitpod_command.h b/src/ML-POD/fitpod_command.h
index 3cbc84fdbd..7001346522 100644
--- a/src/ML-POD/fitpod_command.h
+++ b/src/ML-POD/fitpod_command.h
@@ -122,6 +122,32 @@ public:
 
   FitPOD(LAMMPS *lmp) : Command(lmp) {}
 
+
+  // functions for collecting/collating arrays
+
+  void print_matrix(const char *desc, int m, int n, int *a, int lda);
+  void print_matrix(const char *desc, int m, int n, double *a, int lda);
+  void print_matrix(const char *desc, int m, int n, double **a, int lda);
+  void podCumsum(int *output, int *input, int length);
+  double podArrayNorm(double *a, int n);
+  double podArrayErrorNorm(double *a, double *b, int n);
+  void podArraySetValue(double *y, double a, int n);
+  void podArrayCopy(double *y, double *x, int n);
+  void podArrayFill(int *output, int start, int length);
+  double podArrayMin(double *a, int n);
+  double podArrayMax(double *a, int n);
+  double podArraySum(double *a, int n);
+  int podArrayMin(int *a, int n);
+  int podArrayMax(int *a, int n);
+  void podKron(double *C, double *A, double *B, double alpha, int M1, int M2);
+  void rotation_matrix(double *Rmat, double alpha, double beta, double gamma);
+  void triclinic_lattice_conversion(double *a, double *b, double *c, double *A, double *B,
+                                    double *C);
+  void matrix33_multiplication(double *xrot, double *Rmat, double *x, int natom);
+  void matrix33_inverse(double *invA, double *A1, double *A2, double *A3);
+
+  // functions for reading input files and fitting
+
   void command(int, char **) override;
   int read_data_file(double *fitting_weights, std::string &file_format, std::string &file_extension,
     std::string &test_path, std::string &training_path, std::string &filenametag, const std::string &data_file);
diff --git a/src/ML-POD/mlpod.cpp b/src/ML-POD/mlpod.cpp
index 6506df4528..2712dd6b65 100644
--- a/src/ML-POD/mlpod.cpp
+++ b/src/ML-POD/mlpod.cpp
@@ -108,43 +108,6 @@ MLPOD::~MLPOD()
 
 // clang-format off
 
-
-void MLPOD::print_matrix(const char *desc, int m, int n, double **a, int /*lda*/ )
-{
-  int i, j;
-  printf( "\n %s\n", desc );
-
-  for( i = 0; i < m; i++ )
-  {
-    for( j = 0; j < n; j++ ) printf( " %6.12f", a[j][i] );
-    printf( "\n" );
-  }
-}
-
-void MLPOD::print_matrix(const char *desc, int m, int n, double *a, int lda )
-{
-  int i, j;
-  printf( "\n %s\n", desc );
-
-  for( i = 0; i < m; i++ )
-  {
-    for( j = 0; j < n; j++ ) printf( " %6.12f", a[i+j*lda] );
-    printf( "\n" );
-  }
-}
-
-void MLPOD::print_matrix(const char *desc, int m, int n, int *a, int lda)
-{
-  int i, j;
-  printf( "\n %s\n", desc );
-
-  for( i = 0; i < m; i++ )
-  {
-    for( j = 0; j < n; j++ ) printf( " %d", a[i+j*lda] );
-    printf( "\n" );
-  }
-}
-
 void MLPOD::podMatMul(double *c, double *a, double *b, int r1, int c1, int c2)
 {
   int i, j, k;
@@ -165,84 +128,6 @@ void MLPOD::podArrayFill(int* output, int start, int length)
                 output[j] = start + j;
 }
 
-double MLPOD::podArraySum(double *a, int n)
-{
-  double e = a[0];
-  for (int i=1; i<n; i++)
-    e += a[i];
-  return e;
-}
-
-double MLPOD::podArrayMin(double *a, int n)
-{
-  double b = a[0];
-  for (int i=1; i<n; i++)
-    if (a[i]<b)
-      b = a[i];
-  return b;
-}
-
-double MLPOD::podArrayMax(double *a, int n)
-{
-  double b = a[0];
-  for (int i=1; i<n; i++)
-    if (a[i]>b)
-      b = a[i];
-  return b;
-}
-
-int MLPOD::podArrayMin(int *a, int n)
-{
-  int b = a[0];
-  for (int i=1; i<n; i++)
-    if (a[i]<b)
-      b = a[i];
-  return b;
-}
-
-int MLPOD::podArrayMax(int *a, int n)
-{
-  int b = a[0];
-  for (int i=1; i<n; i++)
-    if (a[i]>b)
-      b = a[i];
-  return b;
-}
-
-void MLPOD::podKron(double *C, double *A, double *B, double alpha, int M1, int M2)
-{
-  int M = M1*M2;
-  for (int idx=0; idx<M; idx++)
-  {
-    int ib = idx%M2;
-    int ia = (idx-ib)/M2;
-    C[idx] += alpha*A[ia]*B[ib];
-  }
-}
-
-void MLPOD::podCumsum(int* output, int* input, int length)
-{
-        output[0] = 0;
-        for (int j = 1; j < length; ++j)
-                output[j] = input[j - 1] + output[j - 1];
-}
-
-double MLPOD::podArrayNorm(double *a, int n)
-{
-  double e = a[0]*a[0];
-  for (int i=1; i<n; i++)
-    e += a[i]*a[i];
-  return sqrt(e);
-}
-
-double MLPOD::podArrayErrorNorm(double *a, double *b, int n)
-{
-  double e = (a[0]-b[0])*(a[0]-b[0]);
-  for (int i=1; i<n; i++)
-    e += (a[i]-b[i])*(a[i]-b[i]);
-  return sqrt(e);
-}
-
 void MLPOD::podArraySetValue(double *y, double a, int n)
 {
   for (int i=0; i<n; i++)
@@ -255,86 +140,6 @@ void MLPOD::podArrayCopy(double *y, double *x, int n)
     y[i] = x[i];
 }
 
-void MLPOD::rotation_matrix(double *Rmat, double alpha, double beta, double gamma)
-{
-  double ca = cos(alpha);
-  double cb = cos(beta);
-  double cg = cos(gamma);
-  double sa = sin(alpha);
-  double sb = sin(beta);
-  double sg = sin(gamma);
-
-  Rmat[0] = ca*cg*cb - sa*sg;
-  Rmat[3] = -ca*cb*sg - sa*cg;
-  Rmat[6] = ca*sb;
-
-  Rmat[1] = sa*cg*cb + ca*sg;
-  Rmat[4] = -sa*cb*sg + ca*cg;
-  Rmat[7] = sa*sb;
-
-  Rmat[2] = -sb*cg;
-  Rmat[5] = sb*sg;
-  Rmat[8] = cb;
-}
-
-void MLPOD::matrix33_multiplication(double *xrot, double *Rmat, double *x, int natom)
-{
-  double x1, x2, x3;
-  for (int i=0; i < natom; i++) {
-    x1 = x[0 + 3*i];
-    x2 = x[1 + 3*i];
-    x3 = x[2 + 3*i];
-    xrot[0 + 3*i] = Rmat[0]*x1 + Rmat[3]*x2 + Rmat[6]*x3;
-    xrot[1 + 3*i] = Rmat[1]*x1 + Rmat[4]*x2 + Rmat[7]*x3;
-    xrot[2 + 3*i] = Rmat[2]*x1 + Rmat[5]*x2 + Rmat[8]*x3;
-  }
-}
-
-void MLPOD::matrix33_inverse(double *invA, double *A1, double *A2, double *A3)
-{
-  double a11 = A1[0];
-  double a21 = A1[1];
-  double a31 = A1[2];
-  double a12 = A2[0];
-  double a22 = A2[1];
-  double a32 = A2[2];
-  double a13 = A3[0];
-  double a23 = A3[1];
-  double a33 = A3[2];
-  double detA = (a11*a22*a33 - a11*a23*a32 - a12*a21*a33 + a12*a23*a31 + a13*a21*a32 - a13*a22*a31);
-
-  invA[0] = (a22*a33 - a23*a32)/detA;
-  invA[1] = (a23*a31 - a21*a33)/detA;
-  invA[2] = (a21*a32 - a22*a31)/detA;
-  invA[3] = (a13*a32 - a12*a33)/detA;
-  invA[4] = (a11*a33 - a13*a31)/detA;
-  invA[5] = (a12*a31 - a11*a32)/detA;
-  invA[6] = (a12*a23 - a13*a22)/detA;
-  invA[7] = (a13*a21 - a11*a23)/detA;
-  invA[8] = (a11*a22 - a12*a21)/detA;
-}
-
-void MLPOD::triclinic_lattice_conversion(double *a, double *b, double *c, double *A, double *B, double *C)
-{
-  double Anorm = sqrt(A[0]*A[0] + A[1]*A[1] + A[2]*A[2]);
-  double Bnorm = sqrt(B[0]*B[0] + B[1]*B[1] + B[2]*B[2]);
-  double Cnorm = sqrt(C[0]*C[0] + C[1]*C[1] + C[2]*C[2]);
-
-  double Ahat[3];
-  Ahat[0] = A[0]/Anorm; Ahat[1] = A[1]/Anorm; Ahat[2] = A[2]/Anorm;
-
-  double ax = Anorm;
-  double bx = B[0]*Ahat[0] + B[1]*Ahat[1] + B[2]*Ahat[2]; //dot(B,Ahat);
-  double by = sqrt(Bnorm*Bnorm - bx*bx); //sqrt(Bnorm^2 - bx^2);// #norm(cross(Ahat,B));
-  double cx = C[0]*Ahat[0] + C[1]*Ahat[1] + C[2]*Ahat[2]; // dot(C,Ahat);
-  double cy = (B[0]*C[0] + B[1]*C[1] + B[2]*C[2] - bx*cx)/by; // (dot(B, C) - bx*cx)/by;
-  double cz = sqrt(Cnorm*Cnorm - cx*cx - cy*cy); // sqrt(Cnorm^2 - cx^2 - cy^2);
-
-  a[0] = ax; a[1] = 0.0; a[2] = 0.0;
-  b[0] = bx; b[1] = by;  b[2] = 0.0;
-  c[0] = cx; c[1] = cy;  c[2] = cz;
-}
-
 void podsnapshots(double *rbf, double *xij, double *besselparams, double rin, double rcut,
     int besseldegree, int inversedegree, int nbesselpars, int N)
 {
diff --git a/src/ML-POD/mlpod.h b/src/ML-POD/mlpod.h
index 40f20de35b..54e75988be 100644
--- a/src/ML-POD/mlpod.h
+++ b/src/ML-POD/mlpod.h
@@ -253,27 +253,10 @@ class MLPOD : protected Pointers {
 
   // functions for collecting/collating arrays
 
-  void print_matrix(const char *desc, int m, int n, int *a, int lda);
-  void print_matrix(const char *desc, int m, int n, double *a, int lda);
-  void print_matrix(const char *desc, int m, int n, double **a, int lda);
   void podMatMul(double *c, double *a, double *b, int r1, int c1, int c2);
-  void podCumsum(int *output, int *input, int length);
-  double podArrayNorm(double *a, int n);
-  double podArrayErrorNorm(double *a, double *b, int n);
   void podArraySetValue(double *y, double a, int n);
   void podArrayCopy(double *y, double *x, int n);
   void podArrayFill(int *output, int start, int length);
-  double podArrayMin(double *a, int n);
-  double podArrayMax(double *a, int n);
-  double podArraySum(double *a, int n);
-  int podArrayMin(int *a, int n);
-  int podArrayMax(int *a, int n);
-  void podKron(double *C, double *A, double *B, double alpha, int M1, int M2);
-  void rotation_matrix(double *Rmat, double alpha, double beta, double gamma);
-  void triclinic_lattice_conversion(double *a, double *b, double *c, double *A, double *B,
-                                    double *C);
-  void matrix33_multiplication(double *xrot, double *Rmat, double *x, int natom);
-  void matrix33_inverse(double *invA, double *A1, double *A2, double *A3);
 
   // functions for calculating energy and force descriptors