diff --git a/src/ML-RANN/rann_fingerprint_bond.cpp b/src/ML-RANN/rann_fingerprint_bond.cpp
index acc3318ba0..536859dc8f 100644
--- a/src/ML-RANN/rann_fingerprint_bond.cpp
+++ b/src/ML-RANN/rann_fingerprint_bond.cpp
@@ -36,6 +36,8 @@ DISTRIBUTION A. Approved for public release; distribution unlimited. OPSEC#4918
 
 using namespace LAMMPS_NS::RANN;
 
+static constexpr double SMALL = 1.0e-12;
+
 Fingerprint_bond::Fingerprint_bond(PairRANN *_pair) : Fingerprint(_pair)
 {
   n_body_type = 3;
@@ -317,7 +319,8 @@ void Fingerprint_bond::do3bodyfeatureset_singleneighborloop(double * features,do
   i = ilist[ii];
   itype = pair->map[type[i]];
   int f = pair->net[itype].dimensions[0];
-  double expr[jnum][kmax+12];
+  std::vector<double> row(kmax+12, 0.0);
+  std::vector<std::vector<double>> expr(jnum, row);
   int p = kmax;
   int countmb=((mlength)*(mlength+1))>>1;
   // calculate interpolation expr, rinvs and dfc, for each neighbor
@@ -354,13 +357,13 @@ void Fingerprint_bond::do3bodyfeatureset_singleneighborloop(double * features,do
     expr[jj][p+1]=dely*rinvs;
     expr[jj][p+2]=delz*rinvs;
     //Hack to avoid nan when x y or z component of radial vector is exactly 0. Shouldn't affect accuracy.
-    if (expr[jj][p]*expr[jj][p]<0.000000000001) {
+    if (expr[jj][p]*expr[jj][p] < SMALL) {
       expr[jj][p] = 0.000001;
     }
-    if (expr[jj][p+1]*expr[jj][p+1]<0.000000000001) {
+    if (expr[jj][p+1]*expr[jj][p+1] < SMALL) {
       expr[jj][p+1] = 0.000001;
     }
-    if (expr[jj][p+2]*expr[jj][p+2]<0.000000000001) {
+    if (expr[jj][p+2]*expr[jj][p+2] < SMALL) {
       expr[jj][p+2] = 0.000001;
     }
     expr[jj][p+3] = -dfc*expr[jj][p];
@@ -377,7 +380,7 @@ void Fingerprint_bond::do3bodyfeatureset_singleneighborloop(double * features,do
   int kb = kmax;
   int mb = mlength;
   count = startingneuron;
-  double Bb[mb];
+  std::vector<double> Bb(mb, 0.0);
   double dBbx;
   double dBby;
   double dBbz;
@@ -397,7 +400,7 @@ void Fingerprint_bond::do3bodyfeatureset_singleneighborloop(double * features,do
     }
     for (n=0;n<kb;n++) {
       for (a1=0;a1<mb;a1++) {
-        Bb[a1]=0;
+        Bb[a1]=0.0;
       }
       ai = n;
       double y1 = alpha_k[ai]/re;
@@ -409,7 +412,7 @@ void Fingerprint_bond::do3bodyfeatureset_singleneighborloop(double * features,do
           continue;
         }
         double yprod = expr[jj][ai];
-        double *y4 = &expr[jj][p];
+        double *y4 = expr[jj].data() + p;
         for (a2=0;a2<a;a2++) {
           yprod *= y4[M[a2+1]];
         }
@@ -419,10 +422,7 @@ void Fingerprint_bond::do3bodyfeatureset_singleneighborloop(double * features,do
         }
       }
       if (atomtypes[1]!=atomtypes[2]) {//Bb!=Bg
-        double Bg[mb];
-        for (a1=0;a1<mb;a1++) {
-          Bg[a1]=0;
-        }
+        std::vector<double> Bg(mb, 0.0);
         ai = n;
         double y1 = alpha_k[ai]/re;
         //loop over ktype to get Bg
@@ -433,7 +433,7 @@ void Fingerprint_bond::do3bodyfeatureset_singleneighborloop(double * features,do
             continue;
           }
           double yprod = expr[jj][ai];
-          double *y4 = &expr[jj][p];
+          double *y4 = expr[jj].data() + p;
           for (a2=0;a2<a;a2++) {
             yprod *= y4[M[a2+1]];
           }
@@ -450,8 +450,8 @@ void Fingerprint_bond::do3bodyfeatureset_singleneighborloop(double * features,do
           if (atomtypes[2] != nelements && atomtypes[2] != jtype) {
             continue;
           }
-          double *y3 = &expr[jj][p+3];
-          double *y4 = &expr[jj][p];
+          double *y3 = expr[jj].data() + p + 3;
+          double *y4 = expr[jj].data() + p;
           ai = n;
           yprod = expr[jj][ai];
           for (a2=0;a2<a;a2++) {
@@ -477,8 +477,8 @@ void Fingerprint_bond::do3bodyfeatureset_singleneighborloop(double * features,do
           if (atomtypes[1] != nelements && atomtypes[1] != jtype) {
             continue;
           }
-          double *y3 = &expr[jj][p+3];
-          double *y4 = &expr[jj][p];
+          double *y3 = expr[jj].data() + p + 3;
+          double *y4 = expr[jj].data() + p;
           ai = n;
           yprod = expr[jj][ai];
           for (a2=0;a2<a;a2++) {
@@ -512,8 +512,8 @@ void Fingerprint_bond::do3bodyfeatureset_singleneighborloop(double * features,do
           if (atomtypes[1] != nelements && atomtypes[1] != jtype) {
             continue;
           }
-          double *y3 = &expr[jj][p+3];
-          double *y4 = &expr[jj][p];
+          double *y3 = expr[jj].data() + p + 3;
+          double *y4 = expr[jj].data() + p;
           ai = n;
           yprod = expr[jj][ai];
           for (a2=0;a2<a;a2++) {
@@ -572,16 +572,18 @@ void Fingerprint_bond::do3bodyfeatureset_doubleneighborloop(double * features,do
   i = ilist[ii];
   itype = pair->map[type[i]];
   int f = pair->net[itype].dimensions[0];
-  double expr[jnum][kmax];
-  double y[jnum][3];
-  double ri[jnum];
-  double dfc[jnum];
+  std::vector<double> row(kmax, 0.0);
+  std::vector<std::vector<double>> expr(jnum, row);
+  std::vector<double> yrow(3, 0.0);
+  std::vector<std::vector<double>> y(jnum, yrow);
+  std::vector<double> ri(jnum, 0.0);
+  std::vector<double> dfc(jnum, 0.0);
   int kb = kmax;
   int mb = mlength;
-  double c41[kmax];
-  double c51[kmax];
-  double c61[kmax];
-  double ct[kmax];
+  std::vector<double> c41(kmax, 0.0);
+  std::vector<double> c51(kmax, 0.0);
+  std::vector<double> c61(kmax, 0.0);
+  std::vector<double> ct(kmax, 0.0);
   for (jj = 0; jj < jnum; jj++) {
     jtype = tn[jj];
     if (jtypes != nelements && jtypes != jtype && ktypes != nelements && ktypes != jtype) {
diff --git a/src/ML-RANN/rann_fingerprint_bondscreened.cpp b/src/ML-RANN/rann_fingerprint_bondscreened.cpp
index 8dcc2e4389..018bbc6df1 100644
--- a/src/ML-RANN/rann_fingerprint_bondscreened.cpp
+++ b/src/ML-RANN/rann_fingerprint_bondscreened.cpp
@@ -36,6 +36,8 @@ DISTRIBUTION A. Approved for public release; distribution unlimited. OPSEC#4918
 
 using namespace LAMMPS_NS::RANN;
 
+static constexpr double SMALL = 1.0e-12;
+
 Fingerprint_bondscreened::Fingerprint_bondscreened(PairRANN *_pair) : Fingerprint(_pair)
 {
   n_body_type = 3;
@@ -319,7 +321,8 @@ void Fingerprint_bondscreened::do3bodyfeatureset_singleneighborloop(double * fea
   i = ilist[ii];
   itype = pair->map[type[i]];
   int f = pair->net[itype].dimensions[0];
-  double expr[jnum][kmax+12];
+  std::vector<double> row(kmax+12, 0.0);
+  std::vector<std::vector<double>> expr(jnum, row);
   int p = kmax;
   int countmb=((mlength)*(mlength+1))>>1;
   // calculate interpolation expr, rinvs and dfc, for each neighbor
@@ -358,13 +361,13 @@ void Fingerprint_bondscreened::do3bodyfeatureset_singleneighborloop(double * fea
     expr[jj][p+1]=dely*rinvs;
     expr[jj][p+2]=delz*rinvs;
     //Hack to avoid nan when x y or z component of radial vector is exactly 0. Shouldn't affect accuracy.
-    if (expr[jj][p]*expr[jj][p]<0.000000000001) {
+    if (expr[jj][p]*expr[jj][p] < SMALL) {
       expr[jj][p] = 0.000001;
     }
-    if (expr[jj][p+1]*expr[jj][p+1]<0.000000000001) {
+    if (expr[jj][p+1]*expr[jj][p+1] < SMALL) {
       expr[jj][p+1] = 0.000001;
     }
-    if (expr[jj][p+2]*expr[jj][p+2]<0.000000000001) {
+    if (expr[jj][p+2]*expr[jj][p+2] < SMALL) {
       expr[jj][p+2] = 0.000001;
     }
     expr[jj][p+3] = -dfc*expr[jj][p]-dSikx[jj];
@@ -381,20 +384,20 @@ void Fingerprint_bondscreened::do3bodyfeatureset_singleneighborloop(double * fea
   int kb = kmax;
   int mb = mlength;
   count = startingneuron;
-  double Bb[mb];
+  std::vector<double> Bb(mb, 0.0);
   double dBbx;
   double dBby;
   double dBbz;
-  double dBbx1[mb];
-  double dBby1[mb];
-  double dBbz1[mb];
+  std::vector<double> dBbx1(mb, 0.0);
+  std::vector<double> dBby1(mb, 0.0);
+  std::vector<double> dBbz1(mb, 0.0);
   double yprod;
   for (mcount=0;mcount<countmb;mcount++) {
     int *M = Mf[mcount];
     int *_coeffx = coeffx[mcount];
-    int *_coeffy =  coeffy[mcount];
-    int *_coeffz =  coeffz[mcount];
-    int *_coeff =  coeff[mcount];
+    int *_coeffy = coeffy[mcount];
+    int *_coeffz = coeffz[mcount];
+    int *_coeff = coeff[mcount];
     a = mb+1;
     for (a1=0;a1<mb;a1++) {
       if (Mf[mcount][a1+1]==0) {
@@ -420,7 +423,7 @@ void Fingerprint_bondscreened::do3bodyfeatureset_singleneighborloop(double * fea
           continue;
         }
         double yprod = expr[jj][ai];
-        double *y4 = &expr[jj][p];
+        double *y4 = expr[jj].data() + p;
         for (a2=0;a2<a;a2++) {
           yprod *= y4[M[a2+1]];
         }
@@ -430,10 +433,7 @@ void Fingerprint_bondscreened::do3bodyfeatureset_singleneighborloop(double * fea
         }
       }
       if (atomtypes[1]!=atomtypes[2]) {//Bb!=Bg
-        double Bg[mb];
-        for (a1=0;a1<mb;a1++) {
-          Bg[a1]=0;
-        }
+        std::vector<double> Bg(mb, 0.0);
         ai = n;
         double y1 = alpha_k[ai]/re;
         //loop over ktype to get Bg
@@ -445,7 +445,7 @@ void Fingerprint_bondscreened::do3bodyfeatureset_singleneighborloop(double * fea
             continue;
           }
           double yprod = expr[jj][ai];
-          double *y4 = &expr[jj][p];
+          double *y4 = expr[jj].data() + p;
           for (a2=0;a2<a;a2++) {
             yprod *= y4[M[a2+1]];
           }
@@ -463,8 +463,8 @@ void Fingerprint_bondscreened::do3bodyfeatureset_singleneighborloop(double * fea
           if (atomtypes[2] != nelements && atomtypes[2] != jtype) {
             continue;
           }
-          double *y3 = &expr[jj][p+3];
-          double *y4 = &expr[jj][p];
+          double *y3 = expr[jj].data() + p + 3;
+          double *y4 = expr[jj].data() + p;
           ai = n;
           yprod = expr[jj][ai];
           for (a2=0;a2<a;a2++) {
@@ -500,8 +500,8 @@ void Fingerprint_bondscreened::do3bodyfeatureset_singleneighborloop(double * fea
           if (atomtypes[1] != nelements && atomtypes[1] != jtype) {
             continue;
           }
-          double *y3 = &expr[jj][p+3];
-          double *y4 = &expr[jj][p];
+          double *y3 = expr[jj].data() + p + 3;
+          double *y4 = expr[jj].data() + p;
           ai = n;
           yprod = expr[jj][ai];
           for (a2=0;a2<a;a2++) {
@@ -545,8 +545,8 @@ void Fingerprint_bondscreened::do3bodyfeatureset_singleneighborloop(double * fea
           if (atomtypes[1] != nelements && atomtypes[1] != jtype) {
             continue;
           }
-          double *y3 = &expr[jj][p+3];
-          double *y4 = &expr[jj][p];
+          double *y3 = expr[jj].data() + p + 3;
+          double *y4 = expr[jj].data() + p;
           ai = n;
           yprod = expr[jj][ai];
           for (a2=0;a2<a;a2++) {
@@ -615,17 +615,20 @@ void Fingerprint_bondscreened::do3bodyfeatureset_doubleneighborloop(double * fea
   i = ilist[ii];
   itype = pair->map[type[i]];
   int f = pair->net[itype].dimensions[0];
-  double expr[jnum][kmax];
-  double y[jnum][3];
-  double ri[jnum];
-  double dfc[jnum];
+  std::vector<double> row(kmax, 0.0);
+  std::vector<std::vector<double>> expr(jnum, row);
+  std::vector<double> yrow(3, 0.0);
+  std::vector<std::vector<double>> y(jnum, yrow);
+  std::vector<double> ri(jnum, 0.0);
+  std::vector<double> dfc(jnum, 0.0);
   int kb = kmax;
   int mb = mlength;
-  double Bijk[kb][mb];
-  double c41[kmax];
-  double c51[kmax];
-  double c61[kmax];
-  double ct[kmax];
+  std::vector<double> brow(mb, 0.0);
+  std::vector<std::vector<double>> Bijk(kb, brow);
+  std::vector<double> c41(kmax, 0.0);
+  std::vector<double> c51(kmax, 0.0);
+  std::vector<double> c61(kmax, 0.0);
+  std::vector<double> ct(kmax, 0.0);
   for (jj = 0; jj < jnum; jj++) {
     if (Bij[jj]==false) {continue;}
     jtype = tn[jj];
diff --git a/src/ML-RANN/rann_fingerprint_bondscreenedspin.cpp b/src/ML-RANN/rann_fingerprint_bondscreenedspin.cpp
index cbd3b2a6dc..26fa46b4eb 100644
--- a/src/ML-RANN/rann_fingerprint_bondscreenedspin.cpp
+++ b/src/ML-RANN/rann_fingerprint_bondscreenedspin.cpp
@@ -35,6 +35,8 @@ DISTRIBUTION A. Approved for public release; distribution unlimited. OPSEC#4918
 
 using namespace LAMMPS_NS::RANN;
 
+static constexpr double SMALL = 1.0e-12;
+
 Fingerprint_bondscreenedspin::Fingerprint_bondscreenedspin(PairRANN *_pair) : Fingerprint(_pair)
 {
   n_body_type = 3;
@@ -320,7 +322,8 @@ void Fingerprint_bondscreenedspin::do3bodyfeatureset_singleneighborloop(double *
   i = ilist[ii];
   itype = pair->map[type[i]];
   int f = pair->net[itype].dimensions[0];
-  double expr[jnum][kmax+12];
+  std::vector<double> row(kmax+12, 0.0);
+  std::vector<std::vector<double>> expr(jnum, row);
   int p = kmax;
   int countmb=((mlength)*(mlength+1))>>1;
   double *si = sim->s[i];
@@ -360,13 +363,13 @@ void Fingerprint_bondscreenedspin::do3bodyfeatureset_singleneighborloop(double *
     expr[jj][p+1]=dely*rinvs;
     expr[jj][p+2]=delz*rinvs;
     //Hack to avoid nan when x y or z component of radial vector is exactly 0. Shouldn't affect accuracy.
-    if (expr[jj][p]*expr[jj][p]<0.000000000001) {
+    if (expr[jj][p]*expr[jj][p] < SMALL) {
       expr[jj][p] = 0.000001;
     }
-    if (expr[jj][p+1]*expr[jj][p+1]<0.000000000001) {
+    if (expr[jj][p+1]*expr[jj][p+1] < SMALL) {
       expr[jj][p+1] = 0.000001;
     }
-    if (expr[jj][p+2]*expr[jj][p+2]<0.000000000001) {
+    if (expr[jj][p+2]*expr[jj][p+2] < SMALL) {
       expr[jj][p+2] = 0.000001;
     }
     expr[jj][p+3] = -dfc*expr[jj][p]-dSikx[jj];
@@ -383,8 +386,8 @@ void Fingerprint_bondscreenedspin::do3bodyfeatureset_singleneighborloop(double *
   int kb = kmax;
   int mb = mlength;
   count = startingneuron;
-  double Bb[mb];
-  double Bbs[mb];
+  std::vector<double> Bb(mb, 0.0);
+  std::vector<double> Bbs(mb, 0.0);
   double dBbx;
   double dBby;
   double dBbz;
@@ -421,7 +424,7 @@ void Fingerprint_bondscreenedspin::do3bodyfeatureset_singleneighborloop(double *
         double *sj = sim->s[j];
         double sp = si[0]*sj[0]+si[1]*sj[1]+si[2]*sj[2];
         double yprod = expr[jj][ai];
-        double *y4 = &expr[jj][p];
+        double *y4 = expr[jj].data() + p;
         for (a2=0;a2<a;a2++) {
           yprod *= y4[M[a2+1]];
         }
@@ -432,12 +435,8 @@ void Fingerprint_bondscreenedspin::do3bodyfeatureset_singleneighborloop(double *
         }
       }
       if (atomtypes[1]!=atomtypes[2]) {//Bb!=Bg
-        double Bg[mb];
-        double Bgs[mb];
-        for (a1=0;a1<mb;a1++) {
-          Bg[a1]=0;
-          Bgs[a1]=0;
-        }
+        std::vector<double> Bg(mb, 0.0);
+        std::vector<double> Bgs(mb, 0.0);
         ai = n;
         double y1 = alpha_k[ai]/re;
         //loop over ktype to get Bg
@@ -452,7 +451,7 @@ void Fingerprint_bondscreenedspin::do3bodyfeatureset_singleneighborloop(double *
           double *sj = sim->s[j];
           double sp = si[0]*sj[0]+si[1]*sj[1]+si[2]*sj[2];
           double yprod = expr[jj][ai];
-          double *y4 = &expr[jj][p];
+          double *y4 = expr[jj].data() + p;
           for (a2=0;a2<a;a2++) {
             yprod *= y4[M[a2+1]];
           }
@@ -474,8 +473,8 @@ void Fingerprint_bondscreenedspin::do3bodyfeatureset_singleneighborloop(double *
           j = jl[jj];
           double *sj = sim->s[j];
           double sp = si[0]*sj[0]+si[1]*sj[1]+si[2]*sj[2];
-          double *y3 = &expr[jj][p+3];
-          double *y4 = &expr[jj][p];
+          double *y3 = expr[jj].data() + p + 3;
+          double *y4 = expr[jj].data() + p;
           ai = n;
           yprod = expr[jj][ai];
           for (a2=0;a2<a;a2++) {
@@ -517,8 +516,8 @@ void Fingerprint_bondscreenedspin::do3bodyfeatureset_singleneighborloop(double *
           j = jl[jj];
           double *sj = sim->s[j];
           double sp = si[0]*sj[0]+si[1]*sj[1]+si[2]*sj[2];
-          double *y3 = &expr[jj][p+3];
-          double *y4 = &expr[jj][p];
+          double *y3 = expr[jj].data() + p + 3;
+          double *y4 = expr[jj].data() + p;
           ai = n;
           yprod = expr[jj][ai];
           for (a2=0;a2<a;a2++) {
@@ -568,8 +567,8 @@ void Fingerprint_bondscreenedspin::do3bodyfeatureset_singleneighborloop(double *
           j = jl[jj];
           double *sj = sim->s[j];
           double sp = si[0]*sj[0]+si[1]*sj[1]+si[2]*sj[2];
-          double *y3 = &expr[jj][p+3];
-          double *y4 = &expr[jj][p];
+          double *y3 = expr[jj].data() + p + 3;
+          double *y4 = expr[jj].data() + p;
           ai = n;
           yprod = expr[jj][ai];
           for (a2=0;a2<a;a2++) {
@@ -641,18 +640,21 @@ void Fingerprint_bondscreenedspin::do3bodyfeatureset_doubleneighborloop(double *
   i = ilist[ii];
   itype = pair->map[type[i]];
   int f = pair->net[itype].dimensions[0];
-  double expr[jnum][kmax];
-  double y[jnum][3];
-  double ri[jnum];
-  double dfc[jnum];
+  std::vector<double> row(kmax, 0.0);
+  std::vector<std::vector<double>> expr(jnum, row);
+  std::vector<double> yrow(3, 0.0);
+  std::vector<std::vector<double>> y(jnum, yrow);
+  std::vector<double> ri(jnum, 0.0);
+  std::vector<double> dfc(jnum, 0.0);
   int kb = kmax;
   int mb = mlength;
-  double Bijk[kb][mb];
-  double c41[kmax];
-  double c51[kmax];
-  double c61[kmax];
-  double ct[kmax];
   double *si = sim->s[i];
+  std::vector<double> brow(mb, 0.0);
+  std::vector<std::vector<double>> Bijk(kb, brow);
+  std::vector<double> c41(kmax, 0.0);
+  std::vector<double> c51(kmax, 0.0);
+  std::vector<double> c61(kmax, 0.0);
+  std::vector<double> ct(kmax, 0.0);
   for (jj = 0; jj < jnum; jj++) {
     if (Bij[jj]==false) {continue;}
       jtype = tn[jj];
diff --git a/src/ML-RANN/rann_fingerprint_bondspin.cpp b/src/ML-RANN/rann_fingerprint_bondspin.cpp
index 1ccdf8ff62..5ecd37b89c 100644
--- a/src/ML-RANN/rann_fingerprint_bondspin.cpp
+++ b/src/ML-RANN/rann_fingerprint_bondspin.cpp
@@ -35,6 +35,8 @@ DISTRIBUTION A. Approved for public release; distribution unlimited. OPSEC#4918
 
 using namespace LAMMPS_NS::RANN;
 
+static constexpr double SMALL = 1.0e-12;
+
 Fingerprint_bondspin::Fingerprint_bondspin(PairRANN *_pair) : Fingerprint(_pair)
 {
   n_body_type = 3;
@@ -317,7 +319,8 @@ void Fingerprint_bondspin::do3bodyfeatureset_singleneighborloop(double * feature
   i = ilist[ii];
   itype = pair->map[type[i]];
   int f = pair->net[itype].dimensions[0];
-  double expr[jnum][kmax+12];
+  std::vector<double> row(kmax+12, 0.0);
+  std::vector<std::vector<double>> expr(jnum, row);
   int p = kmax;
   int countmb=((mlength)*(mlength+1))>>1;
   double *si = sim->s[i];
@@ -355,13 +358,13 @@ void Fingerprint_bondspin::do3bodyfeatureset_singleneighborloop(double * feature
     expr[jj][p+1]=dely*rinvs;
     expr[jj][p+2]=delz*rinvs;
     //Hack to avoid nan when x y or z component of radial vector is exactly 0. Shouldn't affect accuracy.
-    if (expr[jj][p]*expr[jj][p]<0.000000000001) {
+    if (expr[jj][p]*expr[jj][p] < SMALL) {
       expr[jj][p] = 0.000001;
     }
-    if (expr[jj][p+1]*expr[jj][p+1]<0.000000000001) {
+    if (expr[jj][p+1]*expr[jj][p+1] < SMALL) {
       expr[jj][p+1] = 0.000001;
     }
-    if (expr[jj][p+2]*expr[jj][p+2]<0.000000000001) {
+    if (expr[jj][p+2]*expr[jj][p+2] < SMALL) {
       expr[jj][p+2] = 0.000001;
     }
     expr[jj][p+3] = -dfc*expr[jj][p];
@@ -378,8 +381,8 @@ void Fingerprint_bondspin::do3bodyfeatureset_singleneighborloop(double * feature
   int kb = kmax;
   int mb = mlength;
   count = startingneuron;
-  double Bb[mb];
-  double Bbs[mb];
+  std::vector<double> Bb(mb, 0.0);
+  std::vector<double> Bbs(mb, 0.0);
   double dBbx;
   double dBby;
   double dBbz;
@@ -415,7 +418,7 @@ void Fingerprint_bondspin::do3bodyfeatureset_singleneighborloop(double * feature
         double *sj = sim->s[j];
         double sp = si[0]*sj[0]+si[1]*sj[1]+si[2]*sj[2];
         double yprod = expr[jj][ai];
-        double *y4 = &expr[jj][p];
+        double *y4 = expr[jj].data() + p;
         for (a2=0;a2<a;a2++) {
           yprod *= y4[M[a2+1]];
         }
@@ -426,12 +429,8 @@ void Fingerprint_bondspin::do3bodyfeatureset_singleneighborloop(double * feature
         }
       }
       if (atomtypes[1]!=atomtypes[2]) {//Bb!=Bg
-        double Bg[mb];
-        double Bgs[mb];
-        for (a1=0;a1<mb;a1++) {
-          Bg[a1]=0;
-          Bgs[a1]=0;
-        }
+        std::vector<double> Bg(mb, 0.0);
+        std::vector<double> Bgs(mb, 0.0);
         ai = n;
         double y1 = alpha_k[ai]/re;
         //loop over ktype to get Bg
@@ -445,7 +444,7 @@ void Fingerprint_bondspin::do3bodyfeatureset_singleneighborloop(double * feature
           double *sj = sim->s[j];
           double sp = si[0]*sj[0]+si[1]*sj[1]+si[2]*sj[2];
           double yprod = expr[jj][ai];
-          double *y4 = &expr[jj][p];
+          double *y4 = expr[jj].data() + p;
           for (a2=0;a2<a;a2++) {
             yprod *= y4[M[a2+1]];
           }
@@ -466,8 +465,8 @@ void Fingerprint_bondspin::do3bodyfeatureset_singleneighborloop(double * feature
           j = jl[jj];
           double *sj = sim->s[j];
           double sp = si[0]*sj[0]+si[1]*sj[1]+si[2]*sj[2];
-          double *y3 = &expr[jj][p+3];
-          double *y4 = &expr[jj][p];
+          double *y3 = expr[jj].data() + p + 3;
+          double *y4 = expr[jj].data() + p;
           ai = n;
           yprod = expr[jj][ai];
           for (a2=0;a2<a;a2++) {
@@ -502,8 +501,8 @@ void Fingerprint_bondspin::do3bodyfeatureset_singleneighborloop(double * feature
           j = jl[jj];
           double *sj = sim->s[j];
           double sp = si[0]*sj[0]+si[1]*sj[1]+si[2]*sj[2];
-          double *y3 = &expr[jj][p+3];
-          double *y4 = &expr[jj][p];
+          double *y3 = expr[jj].data() + p + 3;
+          double *y4 = expr[jj].data() + p;
           ai = n;
           yprod = expr[jj][ai];
           for (a2=0;a2<a;a2++) {
@@ -546,8 +545,8 @@ void Fingerprint_bondspin::do3bodyfeatureset_singleneighborloop(double * feature
           j = jl[jj];
           double *sj = sim->s[j];
           double sp = si[0]*sj[0]+si[1]*sj[1]+si[2]*sj[2];
-          double *y3 = &expr[jj][p+3];
-          double *y4 = &expr[jj][p];
+          double *y3 = expr[jj].data() + p + 3;
+          double *y4 = expr[jj].data() + p;
           ai = n;
           yprod = expr[jj][ai];
           for (a2=0;a2<a;a2++) {
@@ -612,16 +611,18 @@ void Fingerprint_bondspin::do3bodyfeatureset_doubleneighborloop(double * feature
   i = ilist[ii];
   itype = pair->map[type[i]];
   int f = pair->net[itype].dimensions[0];
-  double expr[jnum][kmax];
-  double y[jnum][3];
-  double ri[jnum];
-  double dfc[jnum];
+  std::vector<double> row(kmax, 0.0);
+  std::vector<std::vector<double>> expr(jnum, row);
+  std::vector<double> yrow(3, 0.0);
+  std::vector<std::vector<double>> y(jnum, yrow);
+  std::vector<double> ri(jnum, 0.0);
+  std::vector<double> dfc(jnum, 0.0);
   int kb = kmax;
   int mb = mlength;
-  double c41[kmax];
-  double c51[kmax];
-  double c61[kmax];
-  double ct[kmax];
+  std::vector<double> c41(kmax, 0.0);
+  std::vector<double> c51(kmax, 0.0);
+  std::vector<double> c61(kmax, 0.0);
+  std::vector<double> ct(kmax, 0.0);
   double *si = sim->s[i];
   for (jj = 0; jj < jnum; jj++) {
     jtype = tn[jj];