diff --git a/unittest/force-styles/test_pair_style.cpp b/unittest/force-styles/test_pair_style.cpp
index 057f2b5352..fb38de54ad 100644
--- a/unittest/force-styles/test_pair_style.cpp
+++ b/unittest/force-styles/test_pair_style.cpp
@@ -158,7 +158,7 @@ void run_lammps(LAMMPS *lmp)
     command("run 4 post no");
 }
 
-void restart_lammps(LAMMPS *lmp, const TestConfig &cfg)
+void restart_lammps(LAMMPS *lmp, const TestConfig &cfg, bool nofdotr = false)
 {
     // utility lambda to improve readability
     auto command = [&](const std::string &line) {
@@ -180,6 +180,7 @@ void restart_lammps(LAMMPS *lmp, const TestConfig &cfg)
     for (auto &post_command : cfg.post_commands) {
         command(post_command);
     }
+    if (nofdotr) command("pair_modify nofdotr");
 
     command("run 0 post no");
 }
@@ -518,6 +519,37 @@ TEST(PairStyle, plain)
     EXPECT_FP_LE_WITH_EPS(pair->eng_coul, test_config.init_coul, epsilon);
     if (print_stats) std::cerr << "restart_energy stats:" << stats << std::endl;
 
+    if (!verbose) ::testing::internal::CaptureStdout();
+    restart_lammps(lmp, test_config, true);
+    if (!verbose) ::testing::internal::GetCapturedStdout();
+
+    f   = lmp->atom->f;
+    tag = lmp->atom->tag;
+    stats.reset();
+    ASSERT_EQ(nlocal + 1, f_ref.size());
+    for (int i = 0; i < nlocal; ++i) {
+        EXPECT_FP_LE_WITH_EPS(f[i][0], f_ref[tag[i]].x, epsilon);
+        EXPECT_FP_LE_WITH_EPS(f[i][1], f_ref[tag[i]].y, epsilon);
+        EXPECT_FP_LE_WITH_EPS(f[i][2], f_ref[tag[i]].z, epsilon);
+    }
+    if (print_stats) std::cerr << "nofdotr_forces stats:" << stats << std::endl;
+
+    pair   = lmp->force->pair;
+    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);
+    EXPECT_FP_LE_WITH_EPS(stress[2], test_config.init_stress.zz, epsilon);
+    EXPECT_FP_LE_WITH_EPS(stress[3], test_config.init_stress.xy, epsilon);
+    EXPECT_FP_LE_WITH_EPS(stress[4], test_config.init_stress.xz, epsilon);
+    EXPECT_FP_LE_WITH_EPS(stress[5], test_config.init_stress.yz, epsilon);
+    if (print_stats) std::cerr << "nofdotr_stress stats:" << stats << std::endl;
+
+    stats.reset();
+    EXPECT_FP_LE_WITH_EPS(pair->eng_vdwl, test_config.init_vdwl, epsilon);
+    EXPECT_FP_LE_WITH_EPS(pair->eng_coul, test_config.init_coul, epsilon);
+    if (print_stats) std::cerr << "nofdotr_energy stats:" << stats << std::endl;
+
     if (!verbose) ::testing::internal::CaptureStdout();
     data_lammps(lmp, test_config);
     if (!verbose) ::testing::internal::GetCapturedStdout();
@@ -770,6 +802,38 @@ TEST(PairStyle, omp)
         EXPECT_FP_LE_WITH_EPS((pair->eng_vdwl + pair->eng_coul), energy, epsilon);
         if (print_stats) std::cerr << "run_energy  stats, newton off:" << stats << std::endl;
     }
+
+    if (!verbose) ::testing::internal::CaptureStdout();
+    restart_lammps(lmp, test_config, true);
+    if (!verbose) ::testing::internal::GetCapturedStdout();
+
+    f   = lmp->atom->f;
+    tag = lmp->atom->tag;
+    stats.reset();
+    ASSERT_EQ(nlocal + 1, f_ref.size());
+    for (int i = 0; i < nlocal; ++i) {
+        EXPECT_FP_LE_WITH_EPS(f[i][0], f_ref[tag[i]].x, 5 * epsilon);
+        EXPECT_FP_LE_WITH_EPS(f[i][1], f_ref[tag[i]].y, 5 * epsilon);
+        EXPECT_FP_LE_WITH_EPS(f[i][2], f_ref[tag[i]].z, 5 * epsilon);
+    }
+    if (print_stats) std::cerr << "nofdotr_forces stats:" << stats << std::endl;
+
+    pair   = lmp->force->pair;
+    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);
+    EXPECT_FP_LE_WITH_EPS(stress[2], test_config.init_stress.zz, 10 * epsilon);
+    EXPECT_FP_LE_WITH_EPS(stress[3], test_config.init_stress.xy, 10 * epsilon);
+    EXPECT_FP_LE_WITH_EPS(stress[4], test_config.init_stress.xz, 10 * epsilon);
+    EXPECT_FP_LE_WITH_EPS(stress[5], test_config.init_stress.yz, 10 * epsilon);
+    if (print_stats) std::cerr << "nofdotr_stress stats:" << stats << std::endl;
+
+    stats.reset();
+    EXPECT_FP_LE_WITH_EPS(pair->eng_vdwl, test_config.init_vdwl, 5 * epsilon);
+    EXPECT_FP_LE_WITH_EPS(pair->eng_coul, test_config.init_coul, 5 * epsilon);
+    if (print_stats) std::cerr << "nofdotr_energy stats:" << stats << std::endl;
+
     if (!verbose) ::testing::internal::CaptureStdout();
     cleanup_lammps(lmp, test_config);
     if (!verbose) ::testing::internal::GetCapturedStdout();
@@ -780,8 +844,10 @@ TEST(PairStyle, gpu)
     if (!LAMMPS::is_installed_pkg("GPU")) GTEST_SKIP();
     if (test_config.skip_tests.count(test_info_->name())) GTEST_SKIP();
 
-    const char *args_neigh[] = {"PairStyle", "-log", "none", "-echo", "screen", "-nocite", "-sf", "gpu"};
-    const char *args_noneigh[] = {"PairStyle", "-log", "none", "-echo", "screen", "-nocite", "-sf", "gpu", "-pk", "gpu", "0", "neigh", "no"};
+    const char *args_neigh[]   = {"PairStyle", "-log",    "none", "-echo",
+                                "screen",    "-nocite", "-sf",  "gpu"};
+    const char *args_noneigh[] = {"PairStyle", "-log", "none", "-echo", "screen", "-nocite", "-sf",
+                                  "gpu",       "-pk",  "gpu",  "0",     "neigh",  "no"};
 
     char **argv = (char **)args_neigh;
     int argc    = sizeof(args_neigh) / sizeof(char *);
@@ -789,7 +855,7 @@ TEST(PairStyle, gpu)
     // cannot use GPU neighbor list with hybrid pair style (yet)
     if (test_config.pair_style.substr(0, 6) == "hybrid") {
         argv = (char **)args_noneigh;
-        argc    = sizeof(args_noneigh) / sizeof(char *);
+        argc = sizeof(args_noneigh) / sizeof(char *);
     }
 
     ::testing::internal::CaptureStdout();
@@ -816,14 +882,17 @@ TEST(PairStyle, gpu)
 
     // relax error for GPU package depending on precision setting
     double epsilon = test_config.epsilon;
-    if (Info::has_accelerator_feature("GPU","precision","double"))
+    if (Info::has_accelerator_feature("GPU", "precision", "double"))
         epsilon *= 7.5;
-    else if (Info::has_accelerator_feature("GPU","precision","mixed"))
+    else if (Info::has_accelerator_feature("GPU", "precision", "mixed"))
         epsilon *= 5.0e8;
-    else epsilon *= 1.0e10;
-    // relax test precision when using pppm and single precision FFTs, but only when also running with double precision
+    else
+        epsilon *= 1.0e10;
+        // relax test precision when using pppm and single precision FFTs, but only when also
+        // running with double precision
 #if defined(FFT_SINGLE)
-    if (lmp->force->kspace && lmp->force->kspace->compute_flag && Info::has_accelerator_feature("GPU","precision","double"))
+    if (lmp->force->kspace && lmp->force->kspace->compute_flag &&
+        Info::has_accelerator_feature("GPU", "precision", "double"))
         if (utils::strmatch(lmp->force->kspace_style, "^pppm")) epsilon *= 2.0e8;
 #endif
     const std::vector<coord_t> &f_ref = test_config.init_forces;
@@ -1015,6 +1084,37 @@ TEST(PairStyle, intel)
     EXPECT_FP_LE_WITH_EPS((pair->eng_vdwl + pair->eng_coul), energy, epsilon);
     if (print_stats) std::cerr << "run_energy  stats:" << stats << std::endl;
 
+    if (!verbose) ::testing::internal::CaptureStdout();
+    restart_lammps(lmp, test_config, true);
+    if (!verbose) ::testing::internal::GetCapturedStdout();
+
+    f   = lmp->atom->f;
+    tag = lmp->atom->tag;
+    stats.reset();
+    ASSERT_EQ(nlocal + 1, f_ref.size());
+    for (int i = 0; i < nlocal; ++i) {
+        EXPECT_FP_LE_WITH_EPS(f[i][0], f_ref[tag[i]].x, 5 * epsilon);
+        EXPECT_FP_LE_WITH_EPS(f[i][1], f_ref[tag[i]].y, 5 * epsilon);
+        EXPECT_FP_LE_WITH_EPS(f[i][2], f_ref[tag[i]].z, 5 * epsilon);
+    }
+    if (print_stats) std::cerr << "nofdotr_forces stats:" << stats << std::endl;
+
+    pair   = lmp->force->pair;
+    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);
+    EXPECT_FP_LE_WITH_EPS(stress[2], test_config.init_stress.zz, 10 * epsilon);
+    EXPECT_FP_LE_WITH_EPS(stress[3], test_config.init_stress.xy, 10 * epsilon);
+    EXPECT_FP_LE_WITH_EPS(stress[4], test_config.init_stress.xz, 10 * epsilon);
+    EXPECT_FP_LE_WITH_EPS(stress[5], test_config.init_stress.yz, 10 * epsilon);
+    if (print_stats) std::cerr << "nofdotr_stress stats:" << stats << std::endl;
+
+    stats.reset();
+    EXPECT_FP_LE_WITH_EPS(pair->eng_vdwl, test_config.init_vdwl, 5 * epsilon);
+    EXPECT_FP_LE_WITH_EPS(pair->eng_coul, test_config.init_coul, 5 * epsilon);
+    if (print_stats) std::cerr << "nofdotr_energy stats:" << stats << std::endl;
+
     if (!verbose) ::testing::internal::CaptureStdout();
     cleanup_lammps(lmp, test_config);
     if (!verbose) ::testing::internal::GetCapturedStdout();
@@ -1123,6 +1223,37 @@ TEST(PairStyle, opt)
     EXPECT_FP_LE_WITH_EPS((pair->eng_vdwl + pair->eng_coul), energy, epsilon);
     if (print_stats) std::cerr << "run_energy  stats:" << stats << std::endl;
 
+    if (!verbose) ::testing::internal::CaptureStdout();
+    restart_lammps(lmp, test_config, true);
+    if (!verbose) ::testing::internal::GetCapturedStdout();
+
+    f   = lmp->atom->f;
+    tag = lmp->atom->tag;
+    stats.reset();
+    ASSERT_EQ(nlocal + 1, f_ref.size());
+    for (int i = 0; i < nlocal; ++i) {
+        EXPECT_FP_LE_WITH_EPS(f[i][0], f_ref[tag[i]].x, 5 * epsilon);
+        EXPECT_FP_LE_WITH_EPS(f[i][1], f_ref[tag[i]].y, 5 * epsilon);
+        EXPECT_FP_LE_WITH_EPS(f[i][2], f_ref[tag[i]].z, 5 * epsilon);
+    }
+    if (print_stats) std::cerr << "nofdotr_forces stats:" << stats << std::endl;
+
+    pair   = lmp->force->pair;
+    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);
+    EXPECT_FP_LE_WITH_EPS(stress[2], test_config.init_stress.zz, 10 * epsilon);
+    EXPECT_FP_LE_WITH_EPS(stress[3], test_config.init_stress.xy, 10 * epsilon);
+    EXPECT_FP_LE_WITH_EPS(stress[4], test_config.init_stress.xz, 10 * epsilon);
+    EXPECT_FP_LE_WITH_EPS(stress[5], test_config.init_stress.yz, 10 * epsilon);
+    if (print_stats) std::cerr << "nofdotr_stress stats:" << stats << std::endl;
+
+    stats.reset();
+    EXPECT_FP_LE_WITH_EPS(pair->eng_vdwl, test_config.init_vdwl, 5 * epsilon);
+    EXPECT_FP_LE_WITH_EPS(pair->eng_coul, test_config.init_coul, 5 * epsilon);
+    if (print_stats) std::cerr << "nofdotr_energy stats:" << stats << std::endl;
+
     if (!verbose) ::testing::internal::CaptureStdout();
     cleanup_lammps(lmp, test_config);
     if (!verbose) ::testing::internal::GetCapturedStdout();