add single test to angle_style, add support for equilibrium data

unittest/force-styles/angle_style.cpp

@@ -766,3 +766,152 @@ TEST(AngleStyle, omp) {
     if (!verbose) ::testing::internal::GetCapturedStdout();
 };
 
+TEST(AngleStyle, single) {
+    const char *args[] = {"AngleStyle", "-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 nangletypes = lmp->atom->nangletypes;
+    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");
+
+    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("pair_style zero 8.0");
+    command("pair_coeff * *");
+
+    command("angle_style " + test_config.angle_style);
+    Angle *angle = lmp->force->angle;
+
+    for (auto& angle_coeff : test_config.angle_coeff) {
+        command("angle_coeff " + angle_coeff);
+    }
+
+    // create (only) three atoms and one angle
+    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_bonds single/angle 1 1 2 3");
+
+    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];
+
+    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();
+}