diff --git a/unittest/force-styles/angle_style.cpp b/unittest/force-styles/angle_style.cpp
index dd4588fa74..9e75536b8b 100644
--- a/unittest/force-styles/angle_style.cpp
+++ b/unittest/force-styles/angle_style.cpp
@@ -26,6 +26,7 @@
 #include "angle.h"
 #include "atom.h"
 #include "compute.h"
+#include "fmt/format.h"
 #include "force.h"
 #include "info.h"
 #include "input.h"
@@ -40,7 +41,6 @@
 #include <ctime>
 #include <mpi.h>
 
-#include <iostream>
 #include <map>
 #include <string>
 #include <utility>
@@ -230,9 +230,7 @@ void generate_yaml_file(const char *outfile, const TestConfig &config)
         return;
     }
 
-    const int natoms  = lmp->atom->natoms;
-    const int bufsize = 256;
-    char buf[bufsize];
+    const int natoms = lmp->atom->natoms;
     std::string block("");
 
     YamlWriter writer(outfile);
@@ -284,20 +282,16 @@ void generate_yaml_file(const char *outfile, const TestConfig &config)
     writer.emit_block("angle_coeff", block);
 
     // equilibrium angle
-    std::stringstream eqstr;
-    eqstr << lmp->atom->nangletypes;
-    for (int i = 0; i < lmp->atom->nangletypes; ++i) {
-        eqstr << " " << lmp->force->angle->equilibrium_angle(i + 1);
-    }
-    writer.emit("equilibrium", eqstr.str());
+    block = fmt::format("{}", lmp->atom->nangletypes);
+    for (int i = 0; i < lmp->atom->nangletypes; ++i)
+        block += fmt::format(" {}", lmp->force->angle->equilibrium_angle(i + 1));
+    writer.emit("equilibrium", block);
 
     // extract
     block.clear();
-    std::stringstream outstr;
-    for (auto &data : config.extract) {
-        outstr << data.first << " " << data.second << std::endl;
-    }
-    writer.emit_block("extract", outstr.str());
+    for (auto data : config.extract)
+        block += fmt::format("{} {}\n", data.first, data.second);
+    writer.emit_block("extract", block);
 
     // natoms
     writer.emit("natoms", natoms);
@@ -306,19 +300,18 @@ void generate_yaml_file(const char *outfile, const TestConfig &config)
     writer.emit("init_energy", lmp->force->angle->energy);
 
     // init_stress
-    double *stress = lmp->force->angle->virial;
-    snprintf(buf, bufsize, "% 23.16e % 23.16e % 23.16e % 23.16e % 23.16e % 23.16e", stress[0],
-             stress[1], stress[2], stress[3], stress[4], stress[5]);
-    writer.emit_block("init_stress", buf);
+    auto stress = lmp->force->angle->virial;
+    block = fmt::format("{:23.16e} {:23.16e} {:23.16e} {:23.16e} {:23.16e} {:23.16e}", stress[0],
+                        stress[1], stress[2], stress[3], stress[4], stress[5]);
+    writer.emit_block("init_stress", block);
 
     // init_forces
     block.clear();
-    double **f  = lmp->atom->f;
-    tagint *tag = lmp->atom->tag;
-    for (int i = 0; i < natoms; ++i) {
-        snprintf(buf, bufsize, "% 3d % 23.16e % 23.16e % 23.16e\n", (int)tag[i], f[i][0], f[i][1],
-                 f[i][2]);
-        block += buf;
+    auto f   = lmp->atom->f;
+    auto tag = lmp->atom->tag;
+    for (int i = 1; i <= natoms; ++i) {
+        const int j = lmp->atom->map(i);
+        block += fmt::format("{:3} {:23.16e} {:23.16e} {:23.16e}\n", i, f[j][0], f[j][1], f[j][2]);
     }
     writer.emit_block("init_forces", block);
 
@@ -330,17 +323,16 @@ void generate_yaml_file(const char *outfile, const TestConfig &config)
 
     // run_stress
     stress = lmp->force->angle->virial;
-    snprintf(buf, bufsize, "% 23.16e % 23.16e % 23.16e % 23.16e % 23.16e % 23.16e", stress[0],
-             stress[1], stress[2], stress[3], stress[4], stress[5]);
-    writer.emit_block("run_stress", buf);
+    block  = fmt::format("{:23.16e} {:23.16e} {:23.16e} {:23.16e} {:23.16e} {:23.16e}", stress[0],
+                        stress[1], stress[2], stress[3], stress[4], stress[5]);
+    writer.emit_block("run_stress", block);
 
     block.clear();
     f   = lmp->atom->f;
     tag = lmp->atom->tag;
-    for (int i = 0; i < natoms; ++i) {
-        snprintf(buf, bufsize, "% 3d % 23.16e % 23.16e % 23.16e\n", (int)tag[i], f[i][0], f[i][1],
-                 f[i][2]);
-        block += buf;
+    for (int i = 1; i <= natoms; ++i) {
+        const int j = lmp->atom->map(i);
+        block += fmt::format("{:3} {:23.16e} {:23.16e} {:23.16e}\n", i, f[j][0], f[j][1], f[j][2]);
     }
     writer.emit_block("run_forces", block);
 
@@ -378,8 +370,9 @@ TEST(AngleStyle, plain)
     ASSERT_EQ(lmp->atom->natoms, nlocal);
 
     double epsilon = test_config.epsilon;
-    double **f     = lmp->atom->f;
-    tagint *tag    = lmp->atom->tag;
+
+    auto f   = lmp->atom->f;
+    auto tag = lmp->atom->tag;
     ErrorStats stats;
     stats.reset();
     const std::vector<coord_t> &f_ref = test_config.init_forces;
@@ -391,8 +384,8 @@ TEST(AngleStyle, plain)
     }
     if (print_stats) std::cerr << "init_forces stats, newton on: " << stats << std::endl;
 
-    Angle *angle   = lmp->force->angle;
-    double *stress = angle->virial;
+    auto angle  = lmp->force->angle;
+    auto stress = angle->virial;
     stats.reset();
     EXPECT_FP_LE_WITH_EPS(stress[0], test_config.init_stress.xx, epsilon);
     EXPECT_FP_LE_WITH_EPS(stress[1], test_config.init_stress.yy, epsilon);
@@ -410,8 +403,10 @@ TEST(AngleStyle, plain)
     run_lammps(lmp);
     if (!verbose) ::testing::internal::GetCapturedStdout();
 
-    f                                 = lmp->atom->f;
-    stress                            = angle->virial;
+    f      = lmp->atom->f;
+    tag    = lmp->atom->tag;
+    stress = angle->virial;
+
     const std::vector<coord_t> &f_run = test_config.run_forces;
     ASSERT_EQ(nlocal + 1, f_run.size());
     stats.reset();
@@ -477,6 +472,7 @@ TEST(AngleStyle, plain)
         if (!verbose) ::testing::internal::GetCapturedStdout();
 
         f      = lmp->atom->f;
+        tag    = lmp->atom->tag;
         stress = angle->virial;
         stats.reset();
         for (int i = 0; i < nlocal; ++i) {
@@ -601,9 +597,11 @@ TEST(AngleStyle, omp)
     ASSERT_EQ(lmp->atom->natoms, nlocal);
 
     // relax error a bit for USER-OMP package
-    double epsilon                    = 5.0 * test_config.epsilon;
-    double **f                        = lmp->atom->f;
-    tagint *tag                       = lmp->atom->tag;
+    double epsilon = 5.0 * test_config.epsilon;
+
+    auto f   = lmp->atom->f;
+    auto tag = lmp->atom->tag;
+
     const std::vector<coord_t> &f_ref = test_config.init_forces;
     ErrorStats stats;
     stats.reset();
@@ -614,8 +612,8 @@ TEST(AngleStyle, omp)
     }
     if (print_stats) std::cerr << "init_forces stats, newton on: " << stats << std::endl;
 
-    Angle *angle   = lmp->force->angle;
-    double *stress = angle->virial;
+    auto angle  = lmp->force->angle;
+    auto stress = angle->virial;
 
     stats.reset();
     EXPECT_FP_LE_WITH_EPS(stress[0], test_config.init_stress.xx, 10 * epsilon);
@@ -634,8 +632,10 @@ TEST(AngleStyle, omp)
     run_lammps(lmp);
     if (!verbose) ::testing::internal::GetCapturedStdout();
 
-    f                                 = lmp->atom->f;
-    stress                            = angle->virial;
+    f      = lmp->atom->f;
+    tag    = lmp->atom->tag;
+    stress = angle->virial;
+
     const std::vector<coord_t> &f_run = test_config.run_forces;
     ASSERT_EQ(nlocal + 1, f_run.size());
     stats.reset();
@@ -703,7 +703,8 @@ TEST(AngleStyle, omp)
         run_lammps(lmp);
         if (!verbose) ::testing::internal::GetCapturedStdout();
 
-        f = lmp->atom->f;
+        f   = lmp->atom->f;
+        tag = lmp->atom->tag;
         stats.reset();
         for (int i = 0; i < nlocal; ++i) {
             EXPECT_FP_LE_WITH_EPS(f[i][0], f_run[tag[i]].x, 10 * epsilon);
@@ -797,14 +798,9 @@ TEST(AngleStyle, single)
     command("atom_modify map array");
     command("region box block -10.0 10.0 -10.0 10.0 -10.0 10.0 units box");
 
-    char buf[10];
-    std::string cmd("create_box 1 box");
-    cmd += " angle/types ";
-    snprintf(buf, 10, "%d", nangletypes);
-    cmd += buf;
-    cmd += " extra/angle/per/atom 2";
-    cmd += " extra/special/per/atom 2";
-    command(cmd);
+    command(fmt::format("create_box 1 box angle/types {} "
+                        "extra/angle/per/atom 2 extra/special/per/atom 2",
+                        nangletypes));
 
     command("pair_style zero 8.0");
     command("pair_coeff * *");
diff --git a/unittest/force-styles/bond_style.cpp b/unittest/force-styles/bond_style.cpp
index 7ba4a08705..77a9fa98aa 100644
--- a/unittest/force-styles/bond_style.cpp
+++ b/unittest/force-styles/bond_style.cpp
@@ -26,6 +26,7 @@
 #include "atom.h"
 #include "bond.h"
 #include "compute.h"
+#include "fmt/format.h"
 #include "force.h"
 #include "info.h"
 #include "input.h"
@@ -40,7 +41,6 @@
 #include <ctime>
 #include <mpi.h>
 
-#include <iostream>
 #include <map>
 #include <string>
 #include <utility>
@@ -230,9 +230,7 @@ void generate_yaml_file(const char *outfile, const TestConfig &config)
         return;
     }
 
-    const int natoms  = lmp->atom->natoms;
-    const int bufsize = 256;
-    char buf[bufsize];
+    const int natoms = lmp->atom->natoms;
     std::string block("");
 
     YamlWriter writer(outfile);
@@ -284,20 +282,16 @@ void generate_yaml_file(const char *outfile, const TestConfig &config)
     writer.emit_block("bond_coeff", block);
 
     // equilibrium distance
-    std::stringstream eqstr;
-    eqstr << lmp->atom->nbondtypes;
-    for (int i = 0; i < lmp->atom->nbondtypes; ++i) {
-        eqstr << " " << lmp->force->bond->equilibrium_distance(i + 1);
-    }
-    writer.emit("equilibrium", eqstr.str());
+    block = fmt::format("{}", lmp->atom->nbondtypes);
+    for (int i = 0; i < lmp->atom->nbondtypes; ++i)
+        block += fmt::format(" {}", lmp->force->bond->equilibrium_distance(i + 1));
+    writer.emit("equilibrium", block);
 
     // extract
     block.clear();
-    std::stringstream outstr;
-    for (auto &data : config.extract) {
-        outstr << data.first << " " << data.second << std::endl;
-    }
-    writer.emit_block("extract", outstr.str());
+    for (auto data : config.extract)
+        block += fmt::format("{} {}\n", data.first, data.second);
+    writer.emit_block("extract", block);
 
     // natoms
     writer.emit("natoms", natoms);
@@ -306,19 +300,18 @@ void generate_yaml_file(const char *outfile, const TestConfig &config)
     writer.emit("init_energy", lmp->force->bond->energy);
 
     // init_stress
-    double *stress = lmp->force->bond->virial;
-    snprintf(buf, bufsize, "% 23.16e % 23.16e % 23.16e % 23.16e % 23.16e % 23.16e", stress[0],
-             stress[1], stress[2], stress[3], stress[4], stress[5]);
-    writer.emit_block("init_stress", buf);
+    auto stress = lmp->force->bond->virial;
+    block = fmt::format("{:23.16e} {:23.16e} {:23.16e} {:23.16e} {:23.16e} {:23.16e}", stress[0],
+                        stress[1], stress[2], stress[3], stress[4], stress[5]);
+    writer.emit_block("init_stress", block);
 
     // init_forces
     block.clear();
-    double **f  = lmp->atom->f;
-    tagint *tag = lmp->atom->tag;
-    for (int i = 0; i < natoms; ++i) {
-        snprintf(buf, bufsize, "% 3d % 23.16e % 23.16e % 23.16e\n", (int)tag[i], f[i][0], f[i][1],
-                 f[i][2]);
-        block += buf;
+    auto f   = lmp->atom->f;
+    auto tag = lmp->atom->tag;
+    for (int i = 1; i <= natoms; ++i) {
+        const int j = lmp->atom->map(i);
+        block += fmt::format("{:3} {:23.16e} {:23.16e} {:23.16e}\n", i, f[j][0], f[j][1], f[j][2]);
     }
     writer.emit_block("init_forces", block);
 
@@ -330,17 +323,16 @@ void generate_yaml_file(const char *outfile, const TestConfig &config)
 
     // run_stress
     stress = lmp->force->bond->virial;
-    snprintf(buf, bufsize, "% 23.16e % 23.16e % 23.16e % 23.16e % 23.16e % 23.16e", stress[0],
-             stress[1], stress[2], stress[3], stress[4], stress[5]);
-    writer.emit_block("run_stress", buf);
+    block  = fmt::format("{:23.16e} {:23.16e} {:23.16e} {:23.16e} {:23.16e} {:23.16e}", stress[0],
+                        stress[1], stress[2], stress[3], stress[4], stress[5]);
+    writer.emit_block("run_stress", block);
 
     block.clear();
     f   = lmp->atom->f;
     tag = lmp->atom->tag;
-    for (int i = 0; i < natoms; ++i) {
-        snprintf(buf, bufsize, "% 3d % 23.16e % 23.16e % 23.16e\n", (int)tag[i], f[i][0], f[i][1],
-                 f[i][2]);
-        block += buf;
+    for (int i = 1; i <= natoms; ++i) {
+        const int j = lmp->atom->map(i);
+        block += fmt::format("{:3} {:23.16e} {:23.16e} {:23.16e}\n", i, f[j][0], f[j][1], f[j][2]);
     }
     writer.emit_block("run_forces", block);
 
@@ -378,8 +370,9 @@ TEST(BondStyle, plain)
     ASSERT_EQ(lmp->atom->natoms, nlocal);
 
     double epsilon = test_config.epsilon;
-    double **f     = lmp->atom->f;
-    tagint *tag    = lmp->atom->tag;
+
+    auto f   = lmp->atom->f;
+    auto tag = lmp->atom->tag;
     ErrorStats stats;
     stats.reset();
     const std::vector<coord_t> &f_ref = test_config.init_forces;
@@ -391,8 +384,8 @@ TEST(BondStyle, plain)
     }
     if (print_stats) std::cerr << "init_forces stats, newton on: " << stats << std::endl;
 
-    Bond *bond     = lmp->force->bond;
-    double *stress = bond->virial;
+    auto bond   = lmp->force->bond;
+    auto stress = bond->virial;
     stats.reset();
     EXPECT_FP_LE_WITH_EPS(stress[0], test_config.init_stress.xx, epsilon);
     EXPECT_FP_LE_WITH_EPS(stress[1], test_config.init_stress.yy, epsilon);
@@ -410,8 +403,10 @@ TEST(BondStyle, plain)
     run_lammps(lmp);
     if (!verbose) ::testing::internal::GetCapturedStdout();
 
-    f                                 = lmp->atom->f;
-    stress                            = bond->virial;
+    f      = lmp->atom->f;
+    tag    = lmp->atom->tag;
+    stress = bond->virial;
+
     const std::vector<coord_t> &f_run = test_config.run_forces;
     ASSERT_EQ(nlocal + 1, f_run.size());
     stats.reset();
@@ -477,6 +472,7 @@ TEST(BondStyle, plain)
         if (!verbose) ::testing::internal::GetCapturedStdout();
 
         f      = lmp->atom->f;
+        tag    = lmp->atom->tag;
         stress = bond->virial;
         stats.reset();
         for (int i = 0; i < nlocal; ++i) {
@@ -601,9 +597,11 @@ TEST(BondStyle, omp)
     ASSERT_EQ(lmp->atom->natoms, nlocal);
 
     // relax error a bit for USER-OMP package
-    double epsilon                    = 5.0 * test_config.epsilon;
-    double **f                        = lmp->atom->f;
-    tagint *tag                       = lmp->atom->tag;
+    double epsilon = 5.0 * test_config.epsilon;
+
+    auto f   = lmp->atom->f;
+    auto tag = lmp->atom->tag;
+
     const std::vector<coord_t> &f_ref = test_config.init_forces;
     ErrorStats stats;
     stats.reset();
@@ -614,8 +612,8 @@ TEST(BondStyle, omp)
     }
     if (print_stats) std::cerr << "init_forces stats, newton on: " << stats << std::endl;
 
-    Bond *bond     = lmp->force->bond;
-    double *stress = bond->virial;
+    auto bond   = lmp->force->bond;
+    auto stress = bond->virial;
     stats.reset();
     EXPECT_FP_LE_WITH_EPS(stress[0], test_config.init_stress.xx, 10 * epsilon);
     EXPECT_FP_LE_WITH_EPS(stress[1], test_config.init_stress.yy, 10 * epsilon);
@@ -633,8 +631,10 @@ TEST(BondStyle, omp)
     run_lammps(lmp);
     if (!verbose) ::testing::internal::GetCapturedStdout();
 
-    f                                 = lmp->atom->f;
-    stress                            = bond->virial;
+    f      = lmp->atom->f;
+    tag    = lmp->atom->tag;
+    stress = bond->virial;
+
     const std::vector<coord_t> &f_run = test_config.run_forces;
     ASSERT_EQ(nlocal + 1, f_run.size());
     stats.reset();
@@ -702,7 +702,8 @@ TEST(BondStyle, omp)
         run_lammps(lmp);
         if (!verbose) ::testing::internal::GetCapturedStdout();
 
-        f = lmp->atom->f;
+        f   = lmp->atom->f;
+        tag = lmp->atom->tag;
         stats.reset();
         for (int i = 0; i < nlocal; ++i) {
             EXPECT_FP_LE_WITH_EPS(f[i][0], f_run[tag[i]].x, 10 * epsilon);
@@ -796,20 +797,15 @@ TEST(BondStyle, single)
     command("atom_modify map array");
     command("region box block -10.0 10.0 -10.0 10.0 -10.0 10.0 units box");
 
-    char buf[10];
-    std::string cmd("create_box 1 box");
-    cmd += " bond/types ";
-    snprintf(buf, 10, "%d", nbondtypes);
-    cmd += buf;
-    cmd += " extra/bond/per/atom 2";
-    cmd += " extra/special/per/atom 2";
-    command(cmd);
+    command(fmt::format("create_box 1 box bond/types {} "
+                        "extra/bond/per/atom 2 extra/special/per/atom 2",
+                        nbondtypes));
 
     command("pair_style zero 8.0");
     command("pair_coeff * *");
 
     command("bond_style " + test_config.bond_style);
-    Bond *bond = lmp->force->bond;
+    auto bond = lmp->force->bond;
 
     for (auto &bond_coeff : test_config.bond_coeff) {
         command("bond_coeff " + bond_coeff);
@@ -1017,9 +1013,10 @@ TEST(BondStyle, extract)
         }
         GTEST_SKIP();
     }
-    Bond *bond = lmp->force->bond;
-    void *ptr  = nullptr;
-    int dim    = 0;
+
+    auto bond = lmp->force->bond;
+    void *ptr = nullptr;
+    int dim   = 0;
     for (auto extract : test_config.extract) {
         ptr = bond->extract(extract.first.c_str(), dim);
         EXPECT_NE(ptr, nullptr);
diff --git a/unittest/force-styles/pair_style.cpp b/unittest/force-styles/pair_style.cpp
index 04356bdf27..593d411349 100644
--- a/unittest/force-styles/pair_style.cpp
+++ b/unittest/force-styles/pair_style.cpp
@@ -25,6 +25,7 @@
 
 #include "atom.h"
 #include "compute.h"
+#include "fmt/format.h"
 #include "force.h"
 #include "info.h"
 #include "input.h"
@@ -42,7 +43,6 @@
 #include <ctime>
 #include <mpi.h>
 
-#include <iostream>
 #include <map>
 #include <string>
 #include <utility>
@@ -231,9 +231,7 @@ void generate_yaml_file(const char *outfile, const TestConfig &config)
         return;
     }
 
-    const int natoms  = lmp->atom->natoms;
-    const int bufsize = 256;
-    char buf[bufsize];
+    const int natoms = lmp->atom->natoms;
     std::string block("");
 
     YamlWriter writer(outfile);
@@ -252,21 +250,21 @@ void generate_yaml_file(const char *outfile, const TestConfig &config)
 
     // prerequisites
     block.clear();
-    for (auto prerequisite : config.prerequisites) {
+    for (auto &prerequisite : config.prerequisites) {
         block += prerequisite.first + " " + prerequisite.second + "\n";
     }
     writer.emit_block("prerequisites", block);
 
     // pre_commands
     block.clear();
-    for (auto command : config.pre_commands) {
+    for (auto &command : config.pre_commands) {
         block += command + "\n";
     }
     writer.emit_block("pre_commands", block);
 
     // post_commands
     block.clear();
-    for (auto command : config.post_commands) {
+    for (auto &command : config.post_commands) {
         block += command + "\n";
     }
     writer.emit_block("post_commands", block);
@@ -286,11 +284,9 @@ void generate_yaml_file(const char *outfile, const TestConfig &config)
 
     // extract
     block.clear();
-    std::stringstream outstr;
-    for (auto data : config.extract) {
-        outstr << data.first << " " << data.second << std::endl;
-    }
-    writer.emit_block("extract", outstr.str());
+    for (auto data : config.extract)
+        block += fmt::format("{} {}\n", data.first, data.second);
+    writer.emit_block("extract", block);
 
     // natoms
     writer.emit("natoms", natoms);
@@ -302,19 +298,18 @@ void generate_yaml_file(const char *outfile, const TestConfig &config)
     writer.emit("init_coul", lmp->force->pair->eng_coul);
 
     // init_stress
-    double *stress = lmp->force->pair->virial;
-    snprintf(buf, bufsize, "% 23.16e % 23.16e % 23.16e % 23.16e % 23.16e % 23.16e", stress[0],
-             stress[1], stress[2], stress[3], stress[4], stress[5]);
-    writer.emit_block("init_stress", buf);
+    auto stress = lmp->force->pair->virial;
+    block = fmt::format("{:23.16e} {:23.16e} {:23.16e} {:23.16e} {:23.16e} {:23.16e}", stress[0],
+                        stress[1], stress[2], stress[3], stress[4], stress[5]);
+    writer.emit_block("init_stress", block);
 
     // init_forces
     block.clear();
-    double **f  = lmp->atom->f;
-    tagint *tag = lmp->atom->tag;
-    for (int i = 0; i < natoms; ++i) {
-        snprintf(buf, bufsize, "% 3d % 23.16e % 23.16e % 23.16e\n", (int)tag[i], f[i][0], f[i][1],
-                 f[i][2]);
-        block += buf;
+    auto f   = lmp->atom->f;
+    auto tag = lmp->atom->tag;
+    for (int i = 1; i <= natoms; ++i) {
+        const int j = lmp->atom->map(i);
+        block += fmt::format("{:3} {:23.16e} {:23.16e} {:23.16e}\n", i, f[j][0], f[j][1], f[j][2]);
     }
     writer.emit_block("init_forces", block);
 
@@ -329,17 +324,16 @@ void generate_yaml_file(const char *outfile, const TestConfig &config)
 
     // run_stress
     stress = lmp->force->pair->virial;
-    snprintf(buf, bufsize, "% 23.16e % 23.16e % 23.16e % 23.16e % 23.16e % 23.16e", stress[0],
-             stress[1], stress[2], stress[3], stress[4], stress[5]);
-    writer.emit_block("run_stress", buf);
+    block  = fmt::format("{:23.16e} {:23.16e} {:23.16e} {:23.16e} {:23.16e} {:23.16e}", stress[0],
+                        stress[1], stress[2], stress[3], stress[4], stress[5]);
+    writer.emit_block("run_stress", block);
 
     block.clear();
     f   = lmp->atom->f;
     tag = lmp->atom->tag;
-    for (int i = 0; i < natoms; ++i) {
-        snprintf(buf, bufsize, "% 3d % 23.16e % 23.16e % 23.16e\n", (int)tag[i], f[i][0], f[i][1],
-                 f[i][2]);
-        block += buf;
+    for (int i = 1; i <= natoms; ++i) {
+        const int j = lmp->atom->map(i);
+        block += fmt::format("{:3} {:23.16e} {:23.16e} {:23.16e}\n", i, f[j][0], f[j][1], f[j][2]);
     }
     writer.emit_block("run_forces", block);
 
@@ -363,7 +357,7 @@ TEST(PairStyle, plain)
     if (!lmp) {
         std::cerr << "One or more prerequisite styles are not available "
                      "in this LAMMPS configuration:\n";
-        for (auto prerequisite : test_config.prerequisites) {
+        for (auto &prerequisite : test_config.prerequisites) {
             std::cerr << prerequisite.first << "_style " << prerequisite.second << "\n";
         }
         GTEST_SKIP();
@@ -382,8 +376,8 @@ TEST(PairStyle, plain)
     if (lmp->force->kspace && lmp->force->kspace->compute_flag)
         if (utils::strmatch(lmp->force->kspace_style, "^pppm")) epsilon *= 2.0e8;
 #endif
-    double **f  = lmp->atom->f;
-    tagint *tag = lmp->atom->tag;
+    auto f   = lmp->atom->f;
+    auto tag = lmp->atom->tag;
     ErrorStats stats;
     stats.reset();
     const std::vector<coord_t> &f_ref = test_config.init_forces;
@@ -395,8 +389,8 @@ TEST(PairStyle, plain)
     }
     if (print_stats) std::cerr << "init_forces stats, newton on: " << stats << std::endl;
 
-    Pair *pair     = lmp->force->pair;
-    double *stress = pair->virial;
+    auto pair   = lmp->force->pair;
+    auto stress = pair->virial;
     stats.reset();
     EXPECT_FP_LE_WITH_EPS(stress[0], test_config.init_stress.xx, epsilon);
     EXPECT_FP_LE_WITH_EPS(stress[1], test_config.init_stress.yy, epsilon);
@@ -415,8 +409,10 @@ TEST(PairStyle, plain)
     run_lammps(lmp);
     if (!verbose) ::testing::internal::GetCapturedStdout();
 
-    f                                 = lmp->atom->f;
-    stress                            = pair->virial;
+    f      = lmp->atom->f;
+    tag    = lmp->atom->tag;
+    stress = pair->virial;
+
     const std::vector<coord_t> &f_run = test_config.run_forces;
     ASSERT_EQ(nlocal + 1, f_run.size());
     stats.reset();
@@ -486,6 +482,7 @@ TEST(PairStyle, plain)
         if (!verbose) ::testing::internal::GetCapturedStdout();
 
         f      = lmp->atom->f;
+        tag    = lmp->atom->tag;
         stress = pair->virial;
         stats.reset();
         for (int i = 0; i < nlocal; ++i) {
@@ -643,7 +640,7 @@ TEST(PairStyle, omp)
     if (!lmp) {
         std::cerr << "One or more prerequisite styles with /omp suffix\n"
                      "are not available in this LAMMPS configuration:\n";
-        for (auto prerequisite : test_config.prerequisites) {
+        for (auto &prerequisite : test_config.prerequisites) {
             std::cerr << prerequisite.first << "_style " << prerequisite.second << "\n";
         }
         GTEST_SKIP();
@@ -663,8 +660,9 @@ TEST(PairStyle, omp)
     if (lmp->force->kspace && lmp->force->kspace->compute_flag)
         if (utils::strmatch(lmp->force->kspace_style, "^pppm")) epsilon *= 2.0e8;
 #endif
-    double **f                        = lmp->atom->f;
-    tagint *tag                       = lmp->atom->tag;
+    auto f   = lmp->atom->f;
+    auto tag = lmp->atom->tag;
+
     const std::vector<coord_t> &f_ref = test_config.init_forces;
     ErrorStats stats;
     stats.reset();
@@ -675,8 +673,8 @@ TEST(PairStyle, omp)
     }
     if (print_stats) std::cerr << "init_forces stats, newton on: " << stats << std::endl;
 
-    Pair *pair     = lmp->force->pair;
-    double *stress = pair->virial;
+    auto pair   = lmp->force->pair;
+    auto stress = pair->virial;
     stats.reset();
     EXPECT_FP_LE_WITH_EPS(stress[0], test_config.init_stress.xx, 10 * epsilon);
     EXPECT_FP_LE_WITH_EPS(stress[1], test_config.init_stress.yy, 10 * epsilon);
@@ -695,8 +693,10 @@ TEST(PairStyle, omp)
     run_lammps(lmp);
     if (!verbose) ::testing::internal::GetCapturedStdout();
 
-    f                                 = lmp->atom->f;
-    stress                            = pair->virial;
+    f      = lmp->atom->f;
+    stress = pair->virial;
+    tag    = lmp->atom->tag;
+
     const std::vector<coord_t> &f_run = test_config.run_forces;
     ASSERT_EQ(nlocal + 1, f_run.size());
     stats.reset();
@@ -763,7 +763,8 @@ TEST(PairStyle, omp)
         run_lammps(lmp);
         if (!verbose) ::testing::internal::GetCapturedStdout();
 
-        f = lmp->atom->f;
+        f   = lmp->atom->f;
+        tag = lmp->atom->tag;
         stats.reset();
         for (int i = 0; i < nlocal; ++i) {
             EXPECT_FP_LE_WITH_EPS(f[i][0], f_run[tag[i]].x, 5 * epsilon);
@@ -859,8 +860,9 @@ TEST(PairStyle, intel)
     const int nlocal = lmp->atom->nlocal;
     ASSERT_EQ(lmp->atom->natoms, nlocal);
 
-    double **f                        = lmp->atom->f;
-    tagint *tag                       = lmp->atom->tag;
+    auto f   = lmp->atom->f;
+    auto tag = lmp->atom->tag;
+
     const std::vector<coord_t> &f_ref = test_config.init_forces;
     ErrorStats stats;
     stats.reset();
@@ -891,8 +893,10 @@ TEST(PairStyle, intel)
     run_lammps(lmp);
     if (!verbose) ::testing::internal::GetCapturedStdout();
 
-    f                                 = lmp->atom->f;
-    stress                            = pair->virial;
+    f      = lmp->atom->f;
+    tag    = lmp->atom->tag;
+    stress = pair->virial;
+
     const std::vector<coord_t> &f_run = test_config.run_forces;
     ASSERT_EQ(nlocal + 1, f_run.size());
     stats.reset();
@@ -967,8 +971,9 @@ TEST(PairStyle, opt)
     if (lmp->force->kspace && lmp->force->kspace->compute_flag)
         if (utils::strmatch(lmp->force->kspace_style, "^pppm")) epsilon *= 2.0e8;
 #endif
-    double **f                        = lmp->atom->f;
-    tagint *tag                       = lmp->atom->tag;
+    auto f   = lmp->atom->f;
+    auto tag = lmp->atom->tag;
+
     const std::vector<coord_t> &f_ref = test_config.init_forces;
     ErrorStats stats;
     stats.reset();
@@ -1000,6 +1005,7 @@ TEST(PairStyle, opt)
     if (!verbose) ::testing::internal::GetCapturedStdout();
 
     f                                 = lmp->atom->f;
+    tag                               = lmp->atom->tag;
     stress                            = pair->virial;
     const std::vector<coord_t> &f_run = test_config.run_forces;
     ASSERT_EQ(nlocal + 1, f_run.size());
@@ -1052,7 +1058,7 @@ TEST(PairStyle, single)
     if (!lmp) {
         std::cerr << "One or more prerequisite styles are not available "
                      "in this LAMMPS configuration:\n";
-        for (auto prerequisite : test_config.prerequisites) {
+        for (auto &prerequisite : test_config.prerequisites) {
             std::cerr << prerequisite.first << "_style " << prerequisite.second << "\n";
         }
         test_config.prerequisites.pop_back();
@@ -1134,11 +1140,7 @@ TEST(PairStyle, single)
     command("atom_modify map array");
     command("region box block -10.0 10.0 -10.0 10.0 -10.0 10.0 units box");
 
-    char buf[10];
-    snprintf(buf, 10, "%d", ntypes);
-    std::string cmd("create_box ");
-    cmd += buf;
-    cmd += " box";
+    auto cmd = fmt::format("create_box {} box", ntypes);
     if (molecular) {
         cmd += " bond/types 1"
                " extra/bond/per/atom 1"
@@ -1308,7 +1310,7 @@ TEST(PairStyle, extract)
         GTEST_SKIP();
     }
 
-    Pair *pair = lmp->force->pair;
+    auto pair = lmp->force->pair;
     if (!pair->compute_flag) {
         std::cerr << "Pair style disabled" << std::endl;
         if (!verbose) ::testing::internal::CaptureStdout();