diff --git a/unittest/force-styles/CMakeLists.txt b/unittest/force-styles/CMakeLists.txt
index 1439b06572..cc84352cd8 100644
--- a/unittest/force-styles/CMakeLists.txt
+++ b/unittest/force-styles/CMakeLists.txt
@@ -115,17 +115,6 @@ foreach(TEST ${KSPACE_TESTS})
   set_tests_properties(${TNAME} PROPERTIES ENVIRONMENT "LAMMPS_POTENTIALS=${LAMMPS_POTENTIALS_DIR};PYTHONPATH=${TEST_INPUT_FOLDER}:$ENV{PYTHONPATH}")
 endforeach()
 
-# improper style tester
-add_executable(test_improper_style test_improper_style.cpp)
-target_link_libraries(test_improper_style PRIVATE lammps style_tests)
-
-file(GLOB IMPROPER_TESTS LIST_DIRECTORIES false ${TEST_INPUT_FOLDER}/improper-*.yaml)
-foreach(TEST ${IMPROPER_TESTS})
-  string(REGEX REPLACE "^.*improper-(.*)\.yaml" "ImproperStyle:\\1" TNAME ${TEST})
-  add_test(NAME ${TNAME} COMMAND test_improper_style ${TEST} WORKING_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR})
-  set_tests_properties(${TNAME} PROPERTIES ENVIRONMENT "LAMMPS_POTENTIALS=${LAMMPS_POTENTIALS_DIR};PYTHONPATH=${TEST_INPUT_FOLDER}:$ENV{PYTHONPATH}")
-endforeach()
-
 # tester for timestepping fixes
 add_executable(test_fix_timestep test_fix_timestep.cpp)
 target_link_libraries(test_fix_timestep PRIVATE lammps style_tests)
@@ -137,3 +126,25 @@ foreach(TEST ${FIX_TIMESTEP_TESTS})
   add_test(NAME ${TNAME} COMMAND test_fix_timestep ${TEST} WORKING_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR})
   set_tests_properties(${TNAME} PROPERTIES ENVIRONMENT "LAMMPS_POTENTIALS=${LAMMPS_POTENTIALS_DIR};PYTHONPATH=${TEST_INPUT_FOLDER}:$ENV{PYTHONPATH}")
 endforeach()
+
+# dihedral style tester
+add_executable(test_dihedral_style test_dihedral_style.cpp)
+target_link_libraries(test_dihedral_style PRIVATE lammps style_tests)
+
+file(GLOB DIHEDRAL_TESTS LIST_DIRECTORIES false ${TEST_INPUT_FOLDER}/dihedral-*.yaml)
+foreach(TEST ${DIHEDRAL_TESTS})
+  string(REGEX REPLACE "^.*dihedral-(.*)\.yaml" "DihedralStyle:\\1" TNAME ${TEST})
+  add_test(NAME ${TNAME} COMMAND test_dihedral_style ${TEST} WORKING_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR})
+  set_tests_properties(${TNAME} PROPERTIES ENVIRONMENT "LAMMPS_POTENTIALS=${LAMMPS_POTENTIALS_DIR};PYTHONPATH=${TEST_INPUT_FOLDER}:$ENV{PYTHONPATH}")
+endforeach()
+
+# improper style tester
+add_executable(test_improper_style test_improper_style.cpp)
+target_link_libraries(test_improper_style PRIVATE lammps style_tests)
+
+file(GLOB IMPROPER_TESTS LIST_DIRECTORIES false ${TEST_INPUT_FOLDER}/improper-*.yaml)
+foreach(TEST ${IMPROPER_TESTS})
+  string(REGEX REPLACE "^.*improper-(.*)\.yaml" "ImproperStyle:\\1" TNAME ${TEST})
+  add_test(NAME ${TNAME} COMMAND test_improper_style ${TEST} WORKING_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR})
+  set_tests_properties(${TNAME} PROPERTIES ENVIRONMENT "LAMMPS_POTENTIALS=${LAMMPS_POTENTIALS_DIR};PYTHONPATH=${TEST_INPUT_FOLDER}:$ENV{PYTHONPATH}")
+endforeach()
diff --git a/unittest/force-styles/test_config_reader.cpp b/unittest/force-styles/test_config_reader.cpp
index e4d5386b7e..2a639784fa 100644
--- a/unittest/force-styles/test_config_reader.cpp
+++ b/unittest/force-styles/test_config_reader.cpp
@@ -57,6 +57,8 @@ TestConfigReader::TestConfigReader(TestConfig &config) : YamlReader(), config(co
     consumers["bond_coeff"]  = &TestConfigReader::bond_coeff;
     consumers["angle_style"] = &TestConfigReader::angle_style;
     consumers["angle_coeff"] = &TestConfigReader::angle_coeff;
+    consumers["dihedral_style"] = &TestConfigReader::dihedral_style;
+    consumers["dihedral_coeff"] = &TestConfigReader::dihedral_coeff;
     consumers["improper_style"] = &TestConfigReader::improper_style;
     consumers["improper_coeff"] = &TestConfigReader::improper_coeff;
     consumers["init_energy"] = &TestConfigReader::init_energy;
@@ -261,6 +263,22 @@ void TestConfigReader::angle_coeff(const yaml_event_t &event)
     }
 }
 
+void TestConfigReader::dihedral_style(const yaml_event_t &event)
+{
+    config.dihedral_style = (char *)event.data.scalar.value;
+}
+
+void TestConfigReader::dihedral_coeff(const yaml_event_t &event)
+{
+    config.dihedral_coeff.clear();
+    std::stringstream data((char *)event.data.scalar.value);
+    std::string line;
+
+    while (std::getline(data, line, '\n')) {
+        config.dihedral_coeff.push_back(line);
+    }
+}
+
 void TestConfigReader::improper_style(const yaml_event_t &event)
 {
     config.improper_style = (char *)event.data.scalar.value;
diff --git a/unittest/force-styles/test_config_reader.h b/unittest/force-styles/test_config_reader.h
index 55cafbcc63..840228da3a 100644
--- a/unittest/force-styles/test_config_reader.h
+++ b/unittest/force-styles/test_config_reader.h
@@ -44,6 +44,8 @@ public:
     void bond_coeff(const yaml_event_t &event);
     void angle_style(const yaml_event_t &event);
     void angle_coeff(const yaml_event_t &event);
+    void dihedral_style(const yaml_event_t &event);
+    void dihedral_coeff(const yaml_event_t &event);
     void improper_style(const yaml_event_t &event);
     void improper_coeff(const yaml_event_t &event);
     void equilibrium(const yaml_event_t &event);
diff --git a/unittest/force-styles/test_dihedral_style.cpp b/unittest/force-styles/test_dihedral_style.cpp
new file mode 100644
index 0000000000..4c6eabb0a5
--- /dev/null
+++ b/unittest/force-styles/test_dihedral_style.cpp
@@ -0,0 +1,749 @@
+/* ----------------------------------------------------------------------
+   LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
+   http://lammps.sandia.gov, Sandia National Laboratories
+   Steve Plimpton, sjplimp@sandia.gov
+
+   Copyright (2003) Sandia Corporation.  Under the terms of Contract
+   DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
+   certain rights in this software.  This software is distributed under
+   the GNU General Public License.
+
+   See the README file in the top-level LAMMPS directory.
+------------------------------------------------------------------------- */
+
+// unit tests for dihedral styles intended for molecular systems
+
+#include "error_stats.h"
+#include "test_config.h"
+#include "test_config_reader.h"
+#include "test_main.h"
+#include "yaml_reader.h"
+#include "yaml_writer.h"
+
+#include "gmock/gmock.h"
+#include "gtest/gtest.h"
+
+#include "atom.h"
+#include "dihedral.h"
+#include "compute.h"
+#include "fmt/format.h"
+#include "force.h"
+#include "info.h"
+#include "input.h"
+#include "lammps.h"
+#include "modify.h"
+#include "universe.h"
+
+#include <cctype>
+#include <cstdio>
+#include <cstdlib>
+#include <cstring>
+#include <ctime>
+#include <mpi.h>
+
+#include <map>
+#include <string>
+#include <utility>
+#include <vector>
+
+using ::testing::HasSubstr;
+using ::testing::StartsWith;
+
+using namespace LAMMPS_NS;
+
+static void delete_file(const std::string &filename)
+{
+    remove(filename.c_str());
+};
+
+// Clean auxilliary files generated during the test
+// which are also useful for debugging failing tests
+void cleanup_lammps(LAMMPS *lmp, const TestConfig &cfg)
+{
+    delete_file(cfg.basename + ".restart");
+    delete_file(cfg.basename + ".data");
+    delete_file(cfg.basename + "-coeffs.in");
+    delete lmp;
+}
+
+// Initialize LAMMPS
+// with the certain arguments, test configuration and an optional flag for newton
+LAMMPS *init_lammps(int argc, char **argv, const TestConfig &cfg, const bool newton = true)
+{
+    LAMMPS *lmp;
+
+    lmp = new LAMMPS(argc, argv, MPI_COMM_WORLD);
+
+    // check if prerequisite styles are available
+    Info *info = new Info(lmp);
+    int nfail  = 0;
+    for (auto &prerequisite : cfg.prerequisites) {
+        std::string style = prerequisite.second;
+
+        // if the suffixed version of dihedral style is not available, don't test it
+        if (prerequisite.first == "dihedral") {
+            if (lmp->suffix_enable) {
+                style += "/";
+                style += lmp->suffix;
+            }
+        }
+
+        if (!info->has_style(prerequisite.first, style)) ++nfail;
+    }
+    delete info;
+    // abort if prerequisites are not met
+    if (nfail > 0) {
+        cleanup_lammps(lmp, cfg);
+        return nullptr;
+    }
+
+    // utility lambdas to improve readability
+    // execute a single-line command
+    auto command = [&](const std::string &line) {
+        lmp->input->one(line.c_str());
+    };
+    // parse and execute all commands in a file
+    auto parse_input_script = [&](const std::string &filename) {
+        lmp->input->file(filename.c_str());
+    };
+
+    if (newton) {
+        command("variable newton_bond index on");
+    } else {
+        command("variable newton_bond index off");
+    }
+
+    command("variable input_dir index " + INPUT_FOLDER);
+
+    for (auto &pre_command : cfg.pre_commands) {
+        command(pre_command);
+    }
+
+    std::string input_file = INPUT_FOLDER + PATH_SEP + cfg.input_file;
+    parse_input_script(input_file);
+
+    command("dihedral_style " + cfg.dihedral_style);
+
+    for (auto &dihedral_coeff : cfg.dihedral_coeff) {
+        command("dihedral_coeff " + dihedral_coeff);
+    }
+
+    for (auto &post_command : cfg.post_commands) {
+        command(post_command);
+    }
+
+    command("run 0 post no");
+    // auxilliary files for running and debugging tests
+    command("write_restart " + cfg.basename + ".restart");
+    command("write_data " + cfg.basename + ".data");
+    command("write_coeff " + cfg.basename + "-coeffs.in");
+
+    return lmp;
+}
+
+// Run a very short NVE simulation
+void run_lammps(LAMMPS *lmp)
+{
+    // utility lambda to improve readability
+    auto command = [&](const std::string &line) {
+        lmp->input->one(line.c_str());
+    };
+
+    command("fix 1 all nve");
+    // just measure the relevant part of potential energy
+    command("compute pe all pe/atom dihedral");
+    command("compute sum all reduce sum c_pe");
+    command("thermo_style custom step temp pe press c_sum");
+    command("thermo 2");
+    command("run 4 post no");
+}
+
+// Restart LAMMPS simulation
+// to test "write_restart" and "read_restart" functions of dihedral styles
+void restart_lammps(LAMMPS *lmp, const TestConfig &cfg)
+{
+    // utility lambda to improve readability
+    auto command = [&](const std::string &line) {
+        lmp->input->one(line.c_str());
+    };
+
+    command("clear");
+    command("read_restart " + cfg.basename + ".restart");
+
+    // add the dihedral style if it's not defined already in the restart file
+    if (!lmp->force->dihedral) {
+        command("dihedral_style " + cfg.dihedral_style);
+    }
+
+    // add the dihedral coefficients if hybrid style is used
+    // or somehow they aren't defined already in the restart file
+    if ((cfg.dihedral_style.substr(0, 6) == "hybrid") || !lmp->force->dihedral->writedata) {
+        for (auto &dihedral_coeff : cfg.dihedral_coeff) {
+            command("dihedral_coeff " + dihedral_coeff);
+        }
+    }
+
+    for (auto &post_command : cfg.post_commands) {
+        command(post_command);
+    }
+
+    command("run 0 post no");
+}
+
+// What's the purpose?
+// Reads the input structure of atoms
+// sets some essential variables
+void data_lammps(LAMMPS *lmp, const TestConfig &cfg)
+{
+    // utility lambdas to improve readability
+    auto command = [&](const std::string &line) {
+        lmp->input->one(line.c_str());
+    };
+    auto parse_input_script = [&](const std::string &filename) {
+        lmp->input->file(filename.c_str());
+    };
+
+    command("clear"); // clears everything except variables, log, echo
+    command("variable dihedral_style delete");
+    command("variable data_file  delete");
+    command("variable newton_bond delete");
+    command("variable newton_bond index on");
+
+    for (auto &pre_command : cfg.pre_commands) {
+        command(pre_command);
+    }
+
+    command("variable dihedral_style index '" + cfg.dihedral_style + "'");
+    command("variable data_file index " + cfg.basename + ".data");
+
+    std::string input_file = INPUT_FOLDER + PATH_SEP + cfg.input_file;
+    parse_input_script(input_file);
+
+    for (auto &dihedral_coeff : cfg.dihedral_coeff) {
+        command("dihedral_coeff " + dihedral_coeff);
+    }
+    for (auto &post_command : cfg.post_commands) {
+        command(post_command);
+    }
+    command("run 0 post no");
+}
+
+void generate_yaml_file(const char *outfile, const TestConfig &config)
+{
+    // initialize system geometry
+    const char *args[] = {"DihedralStyle", "-log", "none", "-echo", "screen", "-nocite"};
+
+    char **argv = (char **)args;
+    int argc    = sizeof(args) / sizeof(char *);
+    LAMMPS *lmp = init_lammps(argc, argv, config);
+    if (!lmp) {
+        std::cerr << "One or more prerequisite styles are not available "
+                     "in this LAMMPS configuration:\n";
+        for (auto &prerequisite : config.prerequisites) {
+            std::cerr << prerequisite.first << "_style " << prerequisite.second << "\n";
+        }
+        return;
+    }
+
+    const int natoms = lmp->atom->natoms;
+    std::string block("");
+
+    YamlWriter writer(outfile);
+
+    // lammps_version
+    writer.emit("lammps_version", lmp->universe->version);
+
+    // date_generated
+    std::time_t now = time(NULL);
+    block           = ctime(&now);
+    block           = block.substr(0, block.find("\n") - 1);
+    writer.emit("date_generated", block);
+
+    // epsilon
+    writer.emit("epsilon", config.epsilon);
+
+    // prerequisites
+    block.clear();
+    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) {
+        block += command + "\n";
+    }
+    writer.emit_block("pre_commands", block);
+
+    // post_commands
+    block.clear();
+    for (auto &command : config.post_commands) {
+        block += command + "\n";
+    }
+    writer.emit_block("post_commands", block);
+
+    // input_file
+    writer.emit("input_file", config.input_file);
+
+    // dihedral_style
+    writer.emit("dihedral_style", config.dihedral_style);
+
+    // dihedral_coeff
+    block.clear();
+    for (auto &dihedral_coeff : config.dihedral_coeff) {
+        block += dihedral_coeff + "\n";
+    }
+    writer.emit_block("dihedral_coeff", block);
+
+    // extract
+    block.clear();
+    for (auto data : config.extract)
+        block += fmt::format("{} {}\n", data.first, data.second);
+    writer.emit_block("extract", block);
+
+    // natoms
+    writer.emit("natoms", natoms);
+
+    // init_energy
+    writer.emit("init_energy", lmp->force->dihedral->energy);
+
+    // init_stress
+    auto stress = lmp->force->dihedral->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();
+    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);
+
+    // do a few steps of MD
+    run_lammps(lmp);
+
+    // run_energy
+    writer.emit("run_energy", lmp->force->dihedral->energy);
+
+    // run_stress
+    stress = lmp->force->dihedral->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("run_stress", block);
+
+    block.clear();
+    f   = lmp->atom->f;
+    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("run_forces", block);
+
+    cleanup_lammps(lmp, config);
+    return;
+}
+
+TEST(ImproperStyle, plain)
+{
+    const char *args[] = {"ImproperStyle", "-log", "none", "-echo", "screen", "-nocite"};
+
+    char **argv = (char **)args;
+    int argc    = sizeof(args) / sizeof(char *);
+
+    ::testing::internal::CaptureStdout();
+    LAMMPS *lmp = init_lammps(argc, argv, test_config, true);
+
+    std::string output = ::testing::internal::GetCapturedStdout();
+    if (verbose) std::cout << output;
+
+    if (!lmp) {
+        std::cerr << "One or more prerequisite styles are not available "
+                     "in this LAMMPS configuration:\n";
+        for (auto &prerequisite : test_config.prerequisites) {
+            std::cerr << prerequisite.first << "_style " << prerequisite.second << "\n";
+        }
+        GTEST_SKIP();
+    }
+
+    EXPECT_THAT(output, StartsWith("LAMMPS ("));
+    EXPECT_THAT(output, HasSubstr("Loop time"));
+
+    // abort if running in parallel and not all atoms are local
+    const int nlocal = lmp->atom->nlocal;
+    ASSERT_EQ(lmp->atom->natoms, nlocal);
+
+    double epsilon = test_config.epsilon;
+
+    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;
+    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 << "init_forces stats, newton on: " << stats << std::endl;
+
+    auto dihedral  = lmp->force->dihedral;
+    auto stress = dihedral->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 << "init_stress stats, newton on: " << stats << std::endl;
+
+    stats.reset();
+    EXPECT_FP_LE_WITH_EPS(dihedral->energy, test_config.init_energy, epsilon);
+    if (print_stats) std::cerr << "init_energy stats, newton on: " << stats << std::endl;
+
+    if (!verbose) ::testing::internal::CaptureStdout();
+    run_lammps(lmp);
+    if (!verbose) ::testing::internal::GetCapturedStdout();
+
+    f      = lmp->atom->f;
+    tag    = lmp->atom->tag;
+    stress = dihedral->virial;
+
+    const std::vector<coord_t> &f_run = test_config.run_forces;
+    ASSERT_EQ(nlocal + 1, f_run.size());
+    stats.reset();
+    for (int i = 0; i < nlocal; ++i) {
+        EXPECT_FP_LE_WITH_EPS(f[i][0], f_run[tag[i]].x, 10 * epsilon);
+        EXPECT_FP_LE_WITH_EPS(f[i][1], f_run[tag[i]].y, 10 * epsilon);
+        EXPECT_FP_LE_WITH_EPS(f[i][2], f_run[tag[i]].z, 10 * epsilon);
+    }
+    if (print_stats) std::cerr << "run_forces  stats, newton on: " << stats << std::endl;
+
+    stress = dihedral->virial;
+    stats.reset();
+    EXPECT_FP_LE_WITH_EPS(stress[0], test_config.run_stress.xx, epsilon);
+    EXPECT_FP_LE_WITH_EPS(stress[1], test_config.run_stress.yy, epsilon);
+    EXPECT_FP_LE_WITH_EPS(stress[2], test_config.run_stress.zz, epsilon);
+    EXPECT_FP_LE_WITH_EPS(stress[3], test_config.run_stress.xy, epsilon);
+    EXPECT_FP_LE_WITH_EPS(stress[4], test_config.run_stress.xz, epsilon);
+    EXPECT_FP_LE_WITH_EPS(stress[5], test_config.run_stress.yz, epsilon);
+    if (print_stats) std::cerr << "run_stress  stats, newton on: " << stats << std::endl;
+
+    stats.reset();
+    int id        = lmp->modify->find_compute("sum");
+    double energy = lmp->modify->compute[id]->compute_scalar();
+    EXPECT_FP_LE_WITH_EPS(dihedral->energy, test_config.run_energy, epsilon);
+    EXPECT_FP_LE_WITH_EPS(dihedral->energy, energy, epsilon);
+    if (print_stats) std::cerr << "run_energy  stats, newton on: " << stats << std::endl;
+
+    if (!verbose) ::testing::internal::CaptureStdout();
+    cleanup_lammps(lmp, test_config);
+    lmp = init_lammps(argc, argv, test_config, false);
+    if (!verbose) ::testing::internal::GetCapturedStdout();
+
+    // skip over these tests if newton bond is forced to be on
+    if (lmp->force->newton_bond == 0) {
+
+        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_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 << "init_forces stats, newton off:" << stats << std::endl;
+
+        dihedral  = lmp->force->dihedral;
+        stress = dihedral->virial;
+        stats.reset();
+        EXPECT_FP_LE_WITH_EPS(stress[0], test_config.init_stress.xx, 2 * epsilon);
+        EXPECT_FP_LE_WITH_EPS(stress[1], test_config.init_stress.yy, 2 * epsilon);
+        EXPECT_FP_LE_WITH_EPS(stress[2], test_config.init_stress.zz, 2 * epsilon);
+        EXPECT_FP_LE_WITH_EPS(stress[3], test_config.init_stress.xy, 2 * epsilon);
+        EXPECT_FP_LE_WITH_EPS(stress[4], test_config.init_stress.xz, 2 * epsilon);
+        EXPECT_FP_LE_WITH_EPS(stress[5], test_config.init_stress.yz, 2 * epsilon);
+        if (print_stats) std::cerr << "init_stress stats, newton off:" << stats << std::endl;
+
+        stats.reset();
+        EXPECT_FP_LE_WITH_EPS(dihedral->energy, test_config.init_energy, epsilon);
+        if (print_stats) std::cerr << "init_energy stats, newton off:" << stats << std::endl;
+
+        if (!verbose) ::testing::internal::CaptureStdout();
+        run_lammps(lmp);
+        if (!verbose) ::testing::internal::GetCapturedStdout();
+
+        f      = lmp->atom->f;
+        tag    = lmp->atom->tag;
+        stress = dihedral->virial;
+        stats.reset();
+        for (int i = 0; i < nlocal; ++i) {
+            EXPECT_FP_LE_WITH_EPS(f[i][0], f_run[tag[i]].x, 10 * epsilon);
+            EXPECT_FP_LE_WITH_EPS(f[i][1], f_run[tag[i]].y, 10 * epsilon);
+            EXPECT_FP_LE_WITH_EPS(f[i][2], f_run[tag[i]].z, 10 * epsilon);
+        }
+        if (print_stats) std::cerr << "run_forces  stats, newton off:" << stats << std::endl;
+
+        stress = dihedral->virial;
+        stats.reset();
+        EXPECT_FP_LE_WITH_EPS(stress[0], test_config.run_stress.xx, epsilon);
+        EXPECT_FP_LE_WITH_EPS(stress[1], test_config.run_stress.yy, epsilon);
+        EXPECT_FP_LE_WITH_EPS(stress[2], test_config.run_stress.zz, epsilon);
+        EXPECT_FP_LE_WITH_EPS(stress[3], test_config.run_stress.xy, epsilon);
+        EXPECT_FP_LE_WITH_EPS(stress[4], test_config.run_stress.xz, epsilon);
+        EXPECT_FP_LE_WITH_EPS(stress[5], test_config.run_stress.yz, epsilon);
+        if (print_stats) std::cerr << "run_stress  stats, newton off:" << stats << std::endl;
+
+        stats.reset();
+        id     = lmp->modify->find_compute("sum");
+        energy = lmp->modify->compute[id]->compute_scalar();
+        EXPECT_FP_LE_WITH_EPS(dihedral->energy, test_config.run_energy, epsilon);
+        EXPECT_FP_LE_WITH_EPS(dihedral->energy, 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);
+    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 << "restart_forces stats:" << stats << std::endl;
+
+    dihedral  = lmp->force->dihedral;
+    stress = dihedral->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 << "restart_stress stats:" << stats << std::endl;
+
+    stats.reset();
+    EXPECT_FP_LE_WITH_EPS(dihedral->energy, test_config.init_energy, epsilon);
+    if (print_stats) std::cerr << "restart_energy stats:" << stats << std::endl;
+
+    if (!verbose) ::testing::internal::CaptureStdout();
+    data_lammps(lmp, test_config);
+    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 << "data_forces stats:" << stats << std::endl;
+
+    dihedral  = lmp->force->dihedral;
+    stress = dihedral->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 << "data_stress stats:" << stats << std::endl;
+
+    stats.reset();
+    EXPECT_FP_LE_WITH_EPS(dihedral->energy, test_config.init_energy, epsilon);
+    if (print_stats) std::cerr << "data_energy stats:" << stats << std::endl;
+
+    if (!verbose) ::testing::internal::CaptureStdout();
+    cleanup_lammps(lmp, test_config);
+    if (!verbose) ::testing::internal::GetCapturedStdout();
+};
+
+TEST(ImproperStyle, omp)
+{
+    if (!LAMMPS::is_installed_pkg("USER-OMP")) GTEST_SKIP();
+    const char *args[] = {"ImproperStyle", "-log", "none", "-echo", "screen", "-nocite",
+                          "-pk",        "omp",  "4",    "-sf",   "omp"};
+
+    char **argv = (char **)args;
+    int argc    = sizeof(args) / sizeof(char *);
+
+    ::testing::internal::CaptureStdout();
+    LAMMPS *lmp = init_lammps(argc, argv, test_config, true);
+
+    std::string output = ::testing::internal::GetCapturedStdout();
+    if (verbose) std::cout << output;
+
+    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) {
+            std::cerr << prerequisite.first << "_style " << prerequisite.second << "\n";
+        }
+        GTEST_SKIP();
+    }
+
+    EXPECT_THAT(output, StartsWith("LAMMPS ("));
+    EXPECT_THAT(output, HasSubstr("Loop time"));
+
+    // abort if running in parallel and not all atoms are local
+    const int nlocal = lmp->atom->nlocal;
+    ASSERT_EQ(lmp->atom->natoms, nlocal);
+
+    // relax error a bit for USER-OMP package
+    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();
+    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 << "init_forces stats, newton on: " << stats << std::endl;
+
+    auto dihedral  = lmp->force->dihedral;
+    auto stress = dihedral->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 << "init_stress stats, newton on: " << stats << std::endl;
+
+    stats.reset();
+    EXPECT_FP_LE_WITH_EPS(dihedral->energy, test_config.init_energy, epsilon);
+    if (print_stats) std::cerr << "init_energy stats, newton on: " << stats << std::endl;
+
+    if (!verbose) ::testing::internal::CaptureStdout();
+    run_lammps(lmp);
+    if (!verbose) ::testing::internal::GetCapturedStdout();
+
+    f      = lmp->atom->f;
+    tag    = lmp->atom->tag;
+    stress = dihedral->virial;
+
+    const std::vector<coord_t> &f_run = test_config.run_forces;
+    ASSERT_EQ(nlocal + 1, f_run.size());
+    stats.reset();
+    for (int i = 0; i < nlocal; ++i) {
+        EXPECT_FP_LE_WITH_EPS(f[i][0], f_run[tag[i]].x, 10 * epsilon);
+        EXPECT_FP_LE_WITH_EPS(f[i][1], f_run[tag[i]].y, 10 * epsilon);
+        EXPECT_FP_LE_WITH_EPS(f[i][2], f_run[tag[i]].z, 10 * epsilon);
+    }
+    if (print_stats) std::cerr << "run_forces  stats, newton on: " << stats << std::endl;
+
+    stress = dihedral->virial;
+    stats.reset();
+    EXPECT_FP_LE_WITH_EPS(stress[0], test_config.run_stress.xx, 10 * epsilon);
+    EXPECT_FP_LE_WITH_EPS(stress[1], test_config.run_stress.yy, 10 * epsilon);
+    EXPECT_FP_LE_WITH_EPS(stress[2], test_config.run_stress.zz, 10 * epsilon);
+    EXPECT_FP_LE_WITH_EPS(stress[3], test_config.run_stress.xy, 10 * epsilon);
+    EXPECT_FP_LE_WITH_EPS(stress[4], test_config.run_stress.xz, 10 * epsilon);
+    EXPECT_FP_LE_WITH_EPS(stress[5], test_config.run_stress.yz, 10 * epsilon);
+    if (print_stats) std::cerr << "run_stress  stats, newton on: " << stats << std::endl;
+
+    stats.reset();
+    int id        = lmp->modify->find_compute("sum");
+    double energy = lmp->modify->compute[id]->compute_scalar();
+    EXPECT_FP_LE_WITH_EPS(dihedral->energy, test_config.run_energy, epsilon);
+    // TODO: this is currently broken for USER-OMP with dihedral style hybrid
+    // needs to be fixed in the main code somewhere. Not sure where, though.
+    if (test_config.dihedral_style.substr(0, 6) != "hybrid")
+        EXPECT_FP_LE_WITH_EPS(dihedral->energy, energy, epsilon);
+    if (print_stats) std::cerr << "run_energy  stats, newton on: " << stats << std::endl;
+
+    if (!verbose) ::testing::internal::CaptureStdout();
+    cleanup_lammps(lmp, test_config);
+    lmp = init_lammps(argc, argv, test_config, false);
+    if (!verbose) ::testing::internal::GetCapturedStdout();
+
+    // skip over these tests if newton bond is forced to be on
+    if (lmp->force->newton_bond == 0) {
+
+        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_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 << "init_forces stats, newton off:" << stats << std::endl;
+
+        dihedral  = lmp->force->dihedral;
+        stress = dihedral->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 << "init_stress stats, newton off:" << stats << std::endl;
+
+        stats.reset();
+        EXPECT_FP_LE_WITH_EPS(dihedral->energy, test_config.init_energy, epsilon);
+        if (print_stats) std::cerr << "init_energy stats, newton off:" << stats << std::endl;
+
+        if (!verbose) ::testing::internal::CaptureStdout();
+        run_lammps(lmp);
+        if (!verbose) ::testing::internal::GetCapturedStdout();
+
+        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);
+            EXPECT_FP_LE_WITH_EPS(f[i][1], f_run[tag[i]].y, 10 * epsilon);
+            EXPECT_FP_LE_WITH_EPS(f[i][2], f_run[tag[i]].z, 10 * epsilon);
+        }
+        if (print_stats) std::cerr << "run_forces  stats, newton off:" << stats << std::endl;
+
+        stress = dihedral->virial;
+        stats.reset();
+        EXPECT_FP_LE_WITH_EPS(stress[0], test_config.run_stress.xx, 10 * epsilon);
+        EXPECT_FP_LE_WITH_EPS(stress[1], test_config.run_stress.yy, 10 * epsilon);
+        EXPECT_FP_LE_WITH_EPS(stress[2], test_config.run_stress.zz, 10 * epsilon);
+        EXPECT_FP_LE_WITH_EPS(stress[3], test_config.run_stress.xy, 10 * epsilon);
+        EXPECT_FP_LE_WITH_EPS(stress[4], test_config.run_stress.xz, 10 * epsilon);
+        EXPECT_FP_LE_WITH_EPS(stress[5], test_config.run_stress.yz, 10 * epsilon);
+        if (print_stats) std::cerr << "run_stress  stats, newton off:" << stats << std::endl;
+
+        stats.reset();
+        id     = lmp->modify->find_compute("sum");
+        energy = lmp->modify->compute[id]->compute_scalar();
+        EXPECT_FP_LE_WITH_EPS(dihedral->energy, test_config.run_energy, epsilon);
+        // TODO: this is currently broken for USER-OMP with dihedral style hybrid
+        // needs to be fixed in the main code somewhere. Not sure where, though.
+        if (test_config.dihedral_style.substr(0, 6) != "hybrid")
+            EXPECT_FP_LE_WITH_EPS(dihedral->energy, energy, epsilon);
+        if (print_stats) std::cerr << "run_energy  stats, newton off:" << stats << std::endl;
+    }
+
+    if (!verbose) ::testing::internal::CaptureStdout();
+    cleanup_lammps(lmp, test_config);
+    if (!verbose) ::testing::internal::GetCapturedStdout();
+};
diff --git a/unittest/force-styles/test_improper_style.cpp b/unittest/force-styles/test_improper_style.cpp
index 9fc5412546..bbc69af791 100644
--- a/unittest/force-styles/test_improper_style.cpp
+++ b/unittest/force-styles/test_improper_style.cpp
@@ -1,19 +1,3 @@
-/*
-- there are 7 basic functions
-  - delete_file
-  - cleanup_lammps
-  - init_lammps
-  - run_lammps
-  - restart_lammps
-  - data_lammps
-  - generate_yaml_file
-- move delete_file amd cleanup_lammps to a single file
-- I don't understand utility lambda, but they are being reused multiple times
-- add as many comments as possible, to show my understanding and document the code
-- code for matching forces, energy and stress are repeated 3 times
-- run_lammps looks to be same across all tests - it isn't, there's subtle difference
-*/
-
 /* ----------------------------------------------------------------------
    LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
    http://lammps.sandia.gov, Sandia National Laboratories
@@ -312,12 +296,6 @@ void generate_yaml_file(const char *outfile, const TestConfig &config)
     }
     writer.emit_block("improper_coeff", block);
 
-    // equilibrium improper
-    // block = fmt::format("{}", lmp->atom->nimpropertypes);
-    // for (int i = 0; i < lmp->atom->nimpropertypes; ++i)
-    //     block += fmt::format(" {}", lmp->force->improper->equilibrium_improper(i + 1));
-    // writer.emit("equilibrium", block);
-
     // extract
     block.clear();
     for (auto data : config.extract)