From a33c0081cd6c64b39e2a037ba5590949f1c6c9b0 Mon Sep 17 00:00:00 2001 From: abhishandy Date: Sat, 8 Aug 2020 15:58:41 -0400 Subject: [PATCH] setup a simple improper-harmonic test --- unittest/force-styles/CMakeLists.txt | 2 +- unittest/force-styles/test_improper_style.cpp | 827 +++++++++++++++++- .../force-styles/tests/improper-harmonic.yaml | 83 ++ 3 files changed, 904 insertions(+), 8 deletions(-) create mode 100644 unittest/force-styles/tests/improper-harmonic.yaml diff --git a/unittest/force-styles/CMakeLists.txt b/unittest/force-styles/CMakeLists.txt index b8628f8c4d..9c2b93c6a9 100644 --- a/unittest/force-styles/CMakeLists.txt +++ b/unittest/force-styles/CMakeLists.txt @@ -96,7 +96,7 @@ endforeach() add_executable(test_angle_style test_angle_style.cpp) target_link_libraries(test_angle_style PRIVATE lammps style_tests) -file(GLOB IMPROPER_TESTS LIST_DIRECTORIES false ${TEST_INPUT_FOLDER}/angle-*.yaml) +file(GLOB ANGLE_TESTS LIST_DIRECTORIES false ${TEST_INPUT_FOLDER}/angle-*.yaml) foreach(TEST ${ANGLE_TESTS}) string(REGEX REPLACE "^.*angle-(.*)\.yaml" "AngleStyle:\\1" TNAME ${TEST}) add_test(NAME ${TNAME} COMMAND test_angle_style ${TEST} WORKING_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR}) diff --git a/unittest/force-styles/test_improper_style.cpp b/unittest/force-styles/test_improper_style.cpp index d40c2b2cfd..a82fe5305e 100644 --- a/unittest/force-styles/test_improper_style.cpp +++ b/unittest/force-styles/test_improper_style.cpp @@ -72,6 +72,8 @@ 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"); @@ -80,6 +82,82 @@ void cleanup_lammps(LAMMPS *lmp, const TestConfig &cfg) 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 improper style is not available, don't test it + if (prerequisite.first == "improper") { + 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("improper_style " + cfg.improper_style); + + for (auto &improper_coeff : cfg.improper_coeff) { + command("improper_coeff " + improper_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 @@ -88,6 +166,7 @@ void run_lammps(LAMMPS *lmp) }; command("fix 1 all nve"); + // just measure the relevant part of potential energy command("compute pe all pe/atom improper"); command("compute sum all reduce sum c_pe"); command("thermo_style custom step temp pe press c_sum"); @@ -95,14 +174,748 @@ void run_lammps(LAMMPS *lmp) command("run 4 post no"); } -LAMMPS *init_lammps(int argc, char **argv, const TestConfig &cfg, const bool newton = true){} +// Restart LAMMPS simulation +// to test "write_restart" and "read_restart" functions of improper 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()); + }; -void restart_lammps(LAMMPS *lmp, const TestConfig &cfg){} + command("clear"); + command("read_restart " + cfg.basename + ".restart"); -void data_lammps(LAMMPS *lmp, const TestConfig &cfg){} + // add the improper style if it's not defined already in the restart file + if (!lmp->force->improper) { + command("improper_style " + cfg.improper_style); + } -void generate_yaml_file(const char *outfile, const TestConfig &config){} + // add the improper coefficients if hybrid style is used + // or somehow they aren't defined already in the restart file + if ((cfg.improper_style.substr(0, 6) == "hybrid") || !lmp->force->improper->writedata) { + for (auto &improper_coeff : cfg.improper_coeff) { + command("improper_coeff " + improper_coeff); + } + } -TEST(ImproperStyle, plain){} -TEST(ImproperStyle, omp){} -TEST(ImproperStyle, single){} + 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 improper_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 improper_style index '" + cfg.improper_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 &improper_coeff : cfg.improper_coeff) { + command("improper_coeff " + improper_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[] = {"ImproperStyle", "-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); + + // improper_style + writer.emit("improper_style", config.improper_style); + + // improper_coeff + block.clear(); + for (auto &improper_coeff : config.improper_coeff) { + block += improper_coeff + "\n"; + } + 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) + 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->improper->energy); + + // init_stress + auto stress = lmp->force->improper->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->improper->energy); + + // run_stress + stress = lmp->force->improper->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 &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 improper = lmp->force->improper; + auto stress = improper->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(improper->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 = improper->virial; + + const std::vector &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 = improper->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(improper->energy, test_config.run_energy, epsilon); + EXPECT_FP_LE_WITH_EPS(improper->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; + + improper = lmp->force->improper; + stress = improper->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(improper->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 = improper->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 = improper->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(improper->energy, test_config.run_energy, epsilon); + EXPECT_FP_LE_WITH_EPS(improper->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; + + improper = lmp->force->improper; + stress = improper->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(improper->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; + + improper = lmp->force->improper; + stress = improper->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(improper->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 &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 improper = lmp->force->improper; + auto stress = improper->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(improper->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 = improper->virial; + + const std::vector &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 = improper->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(improper->energy, test_config.run_energy, epsilon); + // TODO: this is currently broken for USER-OMP with improper style hybrid + // needs to be fixed in the main code somewhere. Not sure where, though. + if (test_config.improper_style.substr(0, 6) != "hybrid") + EXPECT_FP_LE_WITH_EPS(improper->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; + + improper = lmp->force->improper; + stress = improper->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(improper->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 = improper->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(improper->energy, test_config.run_energy, epsilon); + // TODO: this is currently broken for USER-OMP with improper style hybrid + // needs to be fixed in the main code somewhere. Not sure where, though. + if (test_config.improper_style.substr(0, 6) != "hybrid") + EXPECT_FP_LE_WITH_EPS(improper->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(); +}; + +// TODO: Improper styles do not have single() routines +/* +TEST(ImproperStyle, single) +{ + const char *args[] = {"ImproperStyle", "-log", "none", "-echo", "screen", "-nocite"}; + + char **argv = (char **)args; + int argc = sizeof(args) / sizeof(char *); + + // create a LAMMPS instance with standard settings to detect the number of atom types + if (!verbose) ::testing::internal::CaptureStdout(); + LAMMPS *lmp = init_lammps(argc, argv, test_config); + if (!verbose) ::testing::internal::GetCapturedStdout(); + + 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(); + } + + // gather some information and skip if unsupported + int nimpropertypes = lmp->atom->nimpropertypes; + int molecular = lmp->atom->molecular; + if (molecular != 1) { + std::cerr << "Only simple molecular atom styles are supported\n"; + if (!verbose) ::testing::internal::CaptureStdout(); + cleanup_lammps(lmp, test_config); + if (!verbose) ::testing::internal::GetCapturedStdout(); + GTEST_SKIP(); + } + + // utility lambda to improve readability + auto command = [&](const std::string &line) { + lmp->input->one(line.c_str()); + }; + + // now start over + if (!verbose) ::testing::internal::CaptureStdout(); + command("clear"); + command("variable newton_bond delete"); + command("variable newton_bond index on"); + + command("variable input_dir index " + INPUT_FOLDER); + + for (auto &pre_command : test_config.pre_commands) { + command(pre_command); + } + + command("atom_style molecular"); + command("units ${units}"); + command("boundary p p p"); + command("newton ${newton_pair} ${newton_bond}"); + command("special_bonds lj/coul " + "${bond_factor} ${angle_factor} ${dihedral_factor}"); + + command("atom_modify map array"); + command("region box block -10.0 10.0 -10.0 10.0 -10.0 10.0 units box"); + + command(fmt::format("create_box 1 box improper/types {} " + "extra/improper/per/atom 2 extra/special/per/atom 2", + nimpropertypes)); + + command("pair_style zero 8.0"); + command("pair_coeff * *"); + + command("improper_style " + test_config.improper_style); + Improper *improper = lmp->force->improper; + + for (auto &improper_coeff : test_config.improper_coeff) { + command("improper_coeff " + improper_coeff); + } + + // create (only) four atoms and one improper + command("mass * 1.0"); + command("create_atoms 1 single 5.0 -0.75 0.4 units box"); + command("create_atoms 1 single 5.5 0.25 -0.1 units box"); + command("create_atoms 1 single 5.0 0.75 0.4 units box"); + command("create_atoms 1 single 5.5 1.25 -0.1 units box"); + command("create_bonds single/improper 1 1 2 3 4"); + + for (auto &post_command : test_config.post_commands) { + command(post_command); + } + + command("run 0 post no"); + if (!verbose) ::testing::internal::GetCapturedStdout(); + + int idx1 = lmp->atom->map(1); + int idx2 = lmp->atom->map(2); + int idx3 = lmp->atom->map(3); + double epsilon = test_config.epsilon; + double eangle[4], esingle[4]; + // TODO: why are atoms displaced randomly? + + eangle[0] = angle->energy; + esingle[0] = angle->single(1, idx1, idx2, idx3); + + if (!verbose) ::testing::internal::CaptureStdout(); + command("displace_atoms all random 0.5 0.5 0.5 23456"); + command("run 0 post no"); + if (!verbose) ::testing::internal::GetCapturedStdout(); + + idx1 = lmp->atom->map(1); + idx2 = lmp->atom->map(2); + idx3 = lmp->atom->map(3); + eangle[1] = angle->energy; + esingle[1] = angle->single(1, idx1, idx2, idx3); + + if (!verbose) ::testing::internal::CaptureStdout(); + command("displace_atoms all random 0.5 0.5 0.5 456963"); + command("run 0 post no"); + if (!verbose) ::testing::internal::GetCapturedStdout(); + + idx1 = lmp->atom->map(1); + idx2 = lmp->atom->map(2); + idx3 = lmp->atom->map(3); + eangle[2] = angle->energy; + esingle[2] = angle->single(1, idx1, idx2, idx3); + + if (!verbose) ::testing::internal::CaptureStdout(); + command("displace_atoms all random 0.5 0.5 0.5 9726532"); + command("run 0 post no"); + if (!verbose) ::testing::internal::GetCapturedStdout(); + + idx1 = lmp->atom->map(1); + idx2 = lmp->atom->map(2); + idx3 = lmp->atom->map(3); + eangle[3] = angle->energy; + esingle[3] = angle->single(1, idx1, idx2, idx3); + + ErrorStats stats; + EXPECT_FP_LE_WITH_EPS(eangle[0], esingle[0], epsilon); + EXPECT_FP_LE_WITH_EPS(eangle[1], esingle[1], epsilon); + EXPECT_FP_LE_WITH_EPS(eangle[2], esingle[2], epsilon); + EXPECT_FP_LE_WITH_EPS(eangle[3], esingle[3], epsilon); + if (print_stats) std::cerr << "single_energy stats:" << stats << std::endl; + + // int i = 0; + // for (auto &dist : test_config.equilibrium) + // EXPECT_NEAR(dist, angle->equilibrium_angle(++i), 0.00001); + + if (!verbose) ::testing::internal::CaptureStdout(); + cleanup_lammps(lmp, test_config); + if (!verbose) ::testing::internal::GetCapturedStdout(); +} +*/ diff --git a/unittest/force-styles/tests/improper-harmonic.yaml b/unittest/force-styles/tests/improper-harmonic.yaml new file mode 100644 index 0000000000..001216ed17 --- /dev/null +++ b/unittest/force-styles/tests/improper-harmonic.yaml @@ -0,0 +1,83 @@ +--- +lammps_version: 30 Jun 2020 +date_generated: Sun Jul 12 19:14:28 202 +epsilon: 2.5e-13 +prerequisites: ! | + atom full + improper harmonic +pre_commands: ! "" +post_commands: ! "" +input_file: in.fourmol +improper_style: harmonic +improper_coeff: ! | + 1 75.0 0.1 + 2 45.0 1.0 +extract: ! "" +natoms: 29 +init_energy: 41.530817896491 +init_stress: ! |2- + 8.9723357320869255e+01 -8.7188643750026529e+01 -2.5347135708427588e+00 9.2043419883119697e+01 -2.8187238090404989e+01 -1.5291148024927028e+00 +init_forces: ! |2 + 1 4.7865489310693519e+01 7.8760925902181782e+00 -3.2694525514709809e+01 + 2 -1.1124882516177386e+00 -9.0075464203887741e+00 -7.2431691227364725e+00 + 3 -5.9057050592858884e+00 5.3263619873546183e+01 5.2353380124691370e+01 + 4 -1.6032230038990651e+01 -2.4560529343731371e+01 1.2891625920422349e+01 + 5 -4.4802331573497668e+01 -4.8300919461089343e+01 -2.3310767889219321e+01 + 6 4.7083124388174838e+01 -9.5212933434476135e+00 -3.2526392870546800e+01 + 7 -1.6208182775476303e+01 1.4458587960739102e+01 -3.5314745459502710e+00 + 8 -6.5664612141880827e+00 -2.5126850154274180e+01 8.2187944731423329e+01 + 9 -1.5504395262358301e+01 1.6121044185227817e+01 -4.2007069622477866e-01 + 10 9.9863759179364777e+00 4.1873540105704535e+01 -6.6085640966037388e+01 + 11 -2.0441876158908627e+01 -6.5186824168985984e+00 9.0023620309811072e+00 + 12 -1.0772126658369636e+01 -1.0807367300158273e+01 -9.6049647456797036e+00 + 13 2.8847886813946415e+00 7.2973241014859491e+00 -1.0414233993845645e-01 + 14 1.5267407478336423e+01 -9.4754911480231527e+00 -6.6307012925544306e+00 + 15 1.2402914209534794e+01 -6.2644630791613860e+00 1.8484576795819905e+01 + 16 3.8927757686513598e-01 1.0690061587911142e+01 6.1542759189377589e+00 + 17 1.4664194297570587e+00 -1.9971277376602155e+00 1.0776844613215855e+00 + 18 1.5785371874873189e-01 1.6495665212200037e+00 -6.6944747776989910e+00 + 19 -1.9328033033421517e+00 -2.4078805870919515e+00 2.8669575541313312e+00 + 20 1.7749495845934198e+00 7.5831406587194772e-01 3.8275172235676602e+00 + 21 3.4186149299343622e+00 4.2795410364249307e+00 -1.2789555411020601e+01 + 22 -6.0875600315279446e+00 -4.1504951869796436e+00 4.5212856070195588e+00 + 23 2.6689451015935823e+00 -1.2904584944528708e-01 8.2682698040010418e+00 + 24 -1.3053945393770472e+00 5.0741459325182818e+00 -3.0209518576072751e+00 + 25 -1.0471133765834191e+00 -3.5082261409793545e+00 5.7374874908500728e-01 + 26 2.3525079159604663e+00 -1.5659197915389276e+00 2.4472031085222676e+00 + 27 -2.8720725187343144e-01 2.3577465459556626e+00 -8.0312673032167137e-01 + 28 -6.2799575211499037e-01 -1.4097313073755557e+00 3.2747938980615732e-02 + 29 9.1520300398842180e-01 -9.4801523858010661e-01 7.7037879134105569e-01 +run_energy: 29.2286477697792 +run_stress: ! |2- + 6.7161703985479804e+01 -7.4680138065367487e+01 7.5184340798876628e+00 5.7521437901240859e+01 -2.7304190748521741e+01 -1.4932945649428730e+01 +run_forces: ! |2 + 1 3.6220193547187421e+01 1.1585587142479543e+01 -1.7238241972840832e+01 + 2 -1.7637583558698005e+00 -1.3758538851839576e+01 -1.1469071109412811e+01 + 3 -7.1131175159873123e+00 3.1523130842685575e+01 3.9133327910920059e+01 + 4 -1.1598618479874565e+01 -1.4069946914275834e+01 1.1631964860700649e+01 + 5 -3.5092143924598361e+01 -3.4761325046913356e+01 -2.1062123798094685e+01 + 6 4.3880849952881221e+01 -6.1732167988806808e+00 -2.6034788339533922e+01 + 7 -1.4366916728367741e+01 1.3135040103127027e+01 -3.0387136809768127e+00 + 8 -8.5333281419768383e+00 -2.4737111100998256e+01 6.5176870591189868e+01 + 9 -1.2996571868203590e+01 1.3674206710496604e+01 -6.7871105914534047e-01 + 10 1.1736432849972278e+01 3.5147252452549246e+01 -4.9691358934493337e+01 + 11 -1.9930599656448706e+01 -3.2836744898198571e+00 7.6150969595859577e+00 + 12 -5.8293065548538978e+00 -1.3423749355667645e+01 -5.2738511429383701e+00 + 13 7.7658418286980746e-01 2.0512357329017221e+00 1.8932242747136039e+00 + 14 1.2984672111772401e+01 -7.2763363322049823e-01 -4.8588140465034959e+00 + 15 7.4743142834562555e+00 -3.4640965582760748e+00 9.4855052919847029e+00 + 16 3.6043562512330887e+00 8.0382102532623243e+00 4.0048096667552189e+00 + 17 5.4695804680833793e-01 -7.5537048761027403e-01 4.0487452808954988e-01 + 18 -1.1709942604030477e-01 -8.9468235295761567e-01 3.5008479949198765e+00 + 19 9.8256914435044085e-01 1.2515894863018922e+00 -1.5372413162209382e+00 + 20 -8.6546971831013608e-01 -3.5690713334427648e-01 -1.9636066786989381e+00 + 21 6.0524643645662579e-01 7.0314728523699110e-01 -2.2349906198624576e+00 + 22 -1.0299517357238845e+00 -6.7850914711871291e-01 8.1029011311054200e-01 + 23 4.2470529926725875e-01 -2.4638138118278141e-02 1.4247005067519156e+00 + 24 4.1516966455944304e-01 -1.5912776575035221e+00 9.3767616296745859e-01 + 25 3.0070697634261156e-01 1.0957067953103508e+00 -1.8209981159009775e-01 + 26 -7.1587664090205461e-01 4.9557086219317126e-01 -7.5557635137736079e-01 + 27 -1.4554421797154906e-02 1.1056084388051551e-01 -3.5467198058544314e-02 + 28 -2.8257906393508312e-02 -6.6089589395261994e-02 5.7412954785297787e-04 + 29 4.2812328190663218e-02 -4.4471254485253520e-02 3.4893068510691336e-02 +...