1 /* ----------------------------------------------------------------------
2 LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
3 https://www.lammps.org/, Sandia National Laboratories
4 Steve Plimpton, sjplimp@sandia.gov
5
6 Copyright (2003) Sandia Corporation. Under the terms of Contract
7 DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
8 certain rights in this software. This software is distributed under
9 the GNU General Public License.
10
11 See the README file in the top-level LAMMPS directory.
12 ------------------------------------------------------------------------- */
13
14 // unit tests for pair styles intended for molecular systems
15
16 #include "error_stats.h"
17 #include "test_config.h"
18 #include "test_config_reader.h"
19 #include "test_main.h"
20 #include "yaml_reader.h"
21 #include "yaml_writer.h"
22
23 #include "gmock/gmock.h"
24 #include "gtest/gtest.h"
25
26 #include "atom.h"
27 #include "compute.h"
28 #include "force.h"
29 #include "info.h"
30 #include "input.h"
31 #include "kspace.h"
32 #include "lammps.h"
33 #include "modify.h"
34 #include "pair.h"
35 #include "universe.h"
36 #include "utils.h"
37
38 #include <cctype>
39 #include <cstdio>
40 #include <cstdlib>
41 #include <cstring>
42 #include <mpi.h>
43
44 #include <map>
45 #include <string>
46 #include <utility>
47 #include <vector>
48
49 using ::testing::HasSubstr;
50 using ::testing::StartsWith;
51
52 using namespace LAMMPS_NS;
53
delete_file(const std::string & filename)54 static void delete_file(const std::string &filename)
55 {
56 remove(filename.c_str());
57 };
58
cleanup_lammps(LAMMPS * lmp,const TestConfig & cfg)59 void cleanup_lammps(LAMMPS *lmp, const TestConfig &cfg)
60 {
61 delete_file(cfg.basename + ".restart");
62 delete_file(cfg.basename + ".data");
63 delete_file(cfg.basename + "-coeffs.in");
64 delete lmp;
65 lmp = nullptr;
66 }
67
init_lammps(int argc,char ** argv,const TestConfig & cfg,const bool newton=true)68 LAMMPS *init_lammps(int argc, char **argv, const TestConfig &cfg, const bool newton = true)
69 {
70 LAMMPS *lmp;
71
72 lmp = new LAMMPS(argc, argv, MPI_COMM_WORLD);
73
74 // check if prerequisite styles are available
75 Info *info = new Info(lmp);
76 int nfail = 0;
77 for (auto &prerequisite : cfg.prerequisites) {
78 std::string style = prerequisite.second;
79
80 // this is a test for pair styles, so if the suffixed
81 // version is not available, there is no reason to test.
82 if (prerequisite.first == "pair") {
83 if (lmp->suffix_enable) {
84 style += "/";
85 style += lmp->suffix;
86 }
87 }
88
89 if (!info->has_style(prerequisite.first, style)) ++nfail;
90 }
91 delete info;
92 if (nfail > 0) {
93 cleanup_lammps(lmp, cfg);
94 return nullptr;
95 }
96
97 // utility lambdas to improve readability
98 auto command = [&](const std::string &line) {
99 lmp->input->one(line.c_str());
100 };
101 auto parse_input_script = [&](const std::string &filename) {
102 lmp->input->file(filename.c_str());
103 };
104
105 if (newton) {
106 command("variable newton_pair index on");
107 } else {
108 command("variable newton_pair index off");
109 }
110
111 command("variable input_dir index " + INPUT_FOLDER);
112
113 for (auto &pre_command : cfg.pre_commands) {
114 command(pre_command);
115 }
116
117 std::string input_file = INPUT_FOLDER + PATH_SEP + cfg.input_file;
118 parse_input_script(input_file);
119
120 command("pair_style " + cfg.pair_style);
121
122 for (auto &pair_coeff : cfg.pair_coeff) {
123 command("pair_coeff " + pair_coeff);
124 }
125
126 // set this here explicitly to a setting different
127 // from the default, so we can spot YAML files for
128 // long-range interactions that do not include these
129 // settings. they will fail after restart or read data.
130 command("pair_modify table 0");
131 command("pair_modify table/disp 0");
132
133 for (auto &post_command : cfg.post_commands) {
134 command(post_command);
135 }
136
137 command("run 0 post no");
138 command("variable write_data_pair index ii");
139 command("write_restart " + cfg.basename + ".restart");
140 command("write_data " + cfg.basename + ".data pair ${write_data_pair}");
141 command("write_coeff " + cfg.basename + "-coeffs.in");
142
143 return lmp;
144 }
145
run_lammps(LAMMPS * lmp)146 void run_lammps(LAMMPS *lmp)
147 {
148 // utility lambda to improve readability
149 auto command = [&](const std::string &line) {
150 lmp->input->one(line.c_str());
151 };
152
153 command("fix 1 all nve");
154 command("compute pe all pe/atom pair");
155 command("compute sum all reduce sum c_pe");
156 command("thermo_style custom step temp pe press c_sum");
157 command("thermo 2");
158 command("run 4 post no");
159 }
160
restart_lammps(LAMMPS * lmp,const TestConfig & cfg,bool nofdotr=false,bool newton=true)161 void restart_lammps(LAMMPS *lmp, const TestConfig &cfg, bool nofdotr = false, bool newton = true)
162 {
163 // utility lambda to improve readability
164 auto command = [&](const std::string &line) {
165 lmp->input->one(line.c_str());
166 };
167
168 command("clear");
169 if (newton)
170 command("newton on");
171 else
172 command("newton off");
173 command("read_restart " + cfg.basename + ".restart");
174
175 if (!lmp->force->pair) {
176 command("pair_style " + cfg.pair_style);
177 }
178 if (!lmp->force->pair->restartinfo || !lmp->force->pair->writedata) {
179 for (auto &pair_coeff : cfg.pair_coeff) {
180 command("pair_coeff " + pair_coeff);
181 }
182 }
183
184 for (auto &post_command : cfg.post_commands) {
185 command(post_command);
186 }
187 if (nofdotr) command("pair_modify nofdotr");
188
189 command("run 0 post no");
190 }
191
data_lammps(LAMMPS * lmp,const TestConfig & cfg)192 void data_lammps(LAMMPS *lmp, const TestConfig &cfg)
193 {
194 // utility lambdas to improve readability
195 auto command = [&](const std::string &line) {
196 lmp->input->one(line.c_str());
197 };
198 auto parse_input_script = [&](const std::string &filename) {
199 lmp->input->file(filename.c_str());
200 };
201
202 command("clear");
203 command("variable pair_style delete");
204 command("variable data_file delete");
205 command("variable newton_pair delete");
206 command("variable newton_pair index on");
207
208 for (auto &pre_command : cfg.pre_commands) {
209 command(pre_command);
210 }
211
212 command("variable pair_style index '" + cfg.pair_style + "'");
213 command("variable data_file index " + cfg.basename + ".data");
214
215 std::string input_file = INPUT_FOLDER + PATH_SEP + cfg.input_file;
216 parse_input_script(input_file);
217
218 for (auto &pair_coeff : cfg.pair_coeff) {
219 command("pair_coeff " + pair_coeff);
220 }
221 for (auto &post_command : cfg.post_commands) {
222 command(post_command);
223 }
224 command("run 0 post no");
225 }
226
227 // re-generate yaml file with current settings.
228
generate_yaml_file(const char * outfile,const TestConfig & config)229 void generate_yaml_file(const char *outfile, const TestConfig &config)
230 {
231 // initialize system geometry
232 const char *args[] = {"PairStyle", "-log", "none", "-echo", "screen", "-nocite"};
233
234 char **argv = (char **)args;
235 int argc = sizeof(args) / sizeof(char *);
236 LAMMPS *lmp = init_lammps(argc, argv, config);
237 if (!lmp) {
238 std::cerr << "One or more prerequisite styles are not available "
239 "in this LAMMPS configuration:\n";
240 for (auto prerequisite : config.prerequisites) {
241 std::cerr << prerequisite.first << "_style " << prerequisite.second << "\n";
242 }
243 return;
244 }
245
246 const int natoms = lmp->atom->natoms;
247 std::string block("");
248
249 YamlWriter writer(outfile);
250
251 // write yaml header
252 write_yaml_header(&writer, &test_config, lmp->version);
253
254 // pair_style
255 writer.emit("pair_style", config.pair_style);
256
257 // pair_coeff
258 block.clear();
259 for (auto pair_coeff : config.pair_coeff) {
260 block += pair_coeff + "\n";
261 }
262 writer.emit_block("pair_coeff", block);
263
264 // extract
265 block.clear();
266 for (auto data : config.extract)
267 block += fmt::format("{} {}\n", data.first, data.second);
268 writer.emit_block("extract", block);
269
270 // natoms
271 writer.emit("natoms", natoms);
272
273 // init_vdwl
274 writer.emit("init_vdwl", lmp->force->pair->eng_vdwl);
275
276 // init_coul
277 writer.emit("init_coul", lmp->force->pair->eng_coul);
278
279 // init_stress
280 auto stress = lmp->force->pair->virial;
281 // avoid false positives on tiny stresses. force to zero instead.
282 for (int i = 0; i < 6; ++i)
283 if (fabs(stress[i]) < 1.0e-13) stress[i] = 0.0;
284 block = fmt::format("{:23.16e} {:23.16e} {:23.16e} {:23.16e} {:23.16e} {:23.16e}", stress[0],
285 stress[1], stress[2], stress[3], stress[4], stress[5]);
286 writer.emit_block("init_stress", block);
287
288 // init_forces
289 block.clear();
290 auto f = lmp->atom->f;
291 for (int i = 1; i <= natoms; ++i) {
292 const int j = lmp->atom->map(i);
293 block += fmt::format("{:3} {:23.16e} {:23.16e} {:23.16e}\n", i, f[j][0], f[j][1], f[j][2]);
294 }
295 writer.emit_block("init_forces", block);
296
297 // do a few steps of MD
298 run_lammps(lmp);
299
300 // run_vdwl
301 writer.emit("run_vdwl", lmp->force->pair->eng_vdwl);
302
303 // run_coul
304 writer.emit("run_coul", lmp->force->pair->eng_coul);
305
306 // run_stress
307 stress = lmp->force->pair->virial;
308 // avoid false positives on tiny stresses. force to zero instead.
309 for (int i = 0; i < 6; ++i)
310 if (fabs(stress[i]) < 1.0e-13) stress[i] = 0.0;
311 block = fmt::format("{:23.16e} {:23.16e} {:23.16e} {:23.16e} {:23.16e} {:23.16e}", stress[0],
312 stress[1], stress[2], stress[3], stress[4], stress[5]);
313 writer.emit_block("run_stress", block);
314
315 block.clear();
316 f = lmp->atom->f;
317 for (int i = 1; i <= natoms; ++i) {
318 const int j = lmp->atom->map(i);
319 block += fmt::format("{:3} {:23.16e} {:23.16e} {:23.16e}\n", i, f[j][0], f[j][1], f[j][2]);
320 }
321 writer.emit_block("run_forces", block);
322
323 cleanup_lammps(lmp, config);
324 return;
325 }
326
TEST(PairStyle,plain)327 TEST(PairStyle, plain)
328 {
329 if (test_config.skip_tests.count(test_info_->name())) GTEST_SKIP();
330
331 const char *args[] = {"PairStyle", "-log", "none", "-echo", "screen", "-nocite"};
332
333 char **argv = (char **)args;
334 int argc = sizeof(args) / sizeof(char *);
335
336 ::testing::internal::CaptureStdout();
337 LAMMPS *lmp = init_lammps(argc, argv, test_config, true);
338
339 std::string output = ::testing::internal::GetCapturedStdout();
340 if (verbose) std::cout << output;
341
342 if (!lmp) {
343 std::cerr << "One or more prerequisite styles are not available "
344 "in this LAMMPS configuration:\n";
345 for (auto &prerequisite : test_config.prerequisites) {
346 std::cerr << prerequisite.first << "_style " << prerequisite.second << "\n";
347 }
348 GTEST_SKIP();
349 }
350
351 EXPECT_THAT(output, StartsWith("LAMMPS ("));
352 EXPECT_THAT(output, HasSubstr("Loop time"));
353
354 // abort if running in parallel and not all atoms are local
355 const int nlocal = lmp->atom->nlocal;
356 ASSERT_EQ(lmp->atom->natoms, nlocal);
357
358 double epsilon = test_config.epsilon;
359 // relax test precision when using pppm and single precision FFTs
360 #if defined(FFT_SINGLE)
361 if (lmp->force->kspace && lmp->force->kspace->compute_flag)
362 if (utils::strmatch(lmp->force->kspace_style, "^pppm")) epsilon *= 2.0e8;
363 #endif
364 auto f = lmp->atom->f;
365 auto tag = lmp->atom->tag;
366 ErrorStats stats;
367 stats.reset();
368 const std::vector<coord_t> &f_ref = test_config.init_forces;
369 ASSERT_EQ(nlocal + 1, f_ref.size());
370 for (int i = 0; i < nlocal; ++i) {
371 EXPECT_FP_LE_WITH_EPS(f[i][0], f_ref[tag[i]].x, epsilon);
372 EXPECT_FP_LE_WITH_EPS(f[i][1], f_ref[tag[i]].y, epsilon);
373 EXPECT_FP_LE_WITH_EPS(f[i][2], f_ref[tag[i]].z, epsilon);
374 }
375 if (print_stats) std::cerr << "init_forces stats, newton on: " << stats << std::endl;
376
377 auto pair = lmp->force->pair;
378 auto stress = pair->virial;
379 stats.reset();
380 EXPECT_FP_LE_WITH_EPS(stress[0], test_config.init_stress.xx, epsilon);
381 EXPECT_FP_LE_WITH_EPS(stress[1], test_config.init_stress.yy, epsilon);
382 EXPECT_FP_LE_WITH_EPS(stress[2], test_config.init_stress.zz, epsilon);
383 EXPECT_FP_LE_WITH_EPS(stress[3], test_config.init_stress.xy, epsilon);
384 EXPECT_FP_LE_WITH_EPS(stress[4], test_config.init_stress.xz, epsilon);
385 EXPECT_FP_LE_WITH_EPS(stress[5], test_config.init_stress.yz, epsilon);
386 if (print_stats) std::cerr << "init_stress stats, newton on: " << stats << std::endl;
387
388 stats.reset();
389 EXPECT_FP_LE_WITH_EPS(pair->eng_vdwl, test_config.init_vdwl, epsilon);
390 EXPECT_FP_LE_WITH_EPS(pair->eng_coul, test_config.init_coul, epsilon);
391 if (print_stats) std::cerr << "init_energy stats, newton on: " << stats << std::endl;
392
393 if (!verbose) ::testing::internal::CaptureStdout();
394 run_lammps(lmp);
395 if (!verbose) ::testing::internal::GetCapturedStdout();
396
397 f = lmp->atom->f;
398 tag = lmp->atom->tag;
399 stress = pair->virial;
400
401 const std::vector<coord_t> &f_run = test_config.run_forces;
402 ASSERT_EQ(nlocal + 1, f_run.size());
403 stats.reset();
404 for (int i = 0; i < nlocal; ++i) {
405 EXPECT_FP_LE_WITH_EPS(f[i][0], f_run[tag[i]].x, 5 * epsilon);
406 EXPECT_FP_LE_WITH_EPS(f[i][1], f_run[tag[i]].y, 5 * epsilon);
407 EXPECT_FP_LE_WITH_EPS(f[i][2], f_run[tag[i]].z, 5 * epsilon);
408 }
409 if (print_stats) std::cerr << "run_forces stats, newton on: " << stats << std::endl;
410
411 stress = pair->virial;
412 stats.reset();
413 EXPECT_FP_LE_WITH_EPS(stress[0], test_config.run_stress.xx, epsilon);
414 EXPECT_FP_LE_WITH_EPS(stress[1], test_config.run_stress.yy, epsilon);
415 EXPECT_FP_LE_WITH_EPS(stress[2], test_config.run_stress.zz, epsilon);
416 EXPECT_FP_LE_WITH_EPS(stress[3], test_config.run_stress.xy, epsilon);
417 EXPECT_FP_LE_WITH_EPS(stress[4], test_config.run_stress.xz, epsilon);
418 EXPECT_FP_LE_WITH_EPS(stress[5], test_config.run_stress.yz, epsilon);
419 if (print_stats) std::cerr << "run_stress stats, newton on: " << stats << std::endl;
420
421 stats.reset();
422 int id = lmp->modify->find_compute("sum");
423 double energy = lmp->modify->compute[id]->compute_scalar();
424 EXPECT_FP_LE_WITH_EPS(pair->eng_vdwl, test_config.run_vdwl, epsilon);
425 EXPECT_FP_LE_WITH_EPS(pair->eng_coul, test_config.run_coul, epsilon);
426 // skip comparing per-atom energy with total energy for "kim"
427 if (std::string("kim") != lmp->force->pair_style)
428 EXPECT_FP_LE_WITH_EPS((pair->eng_vdwl + pair->eng_coul), energy, epsilon);
429 if (print_stats) std::cerr << "run_energy stats, newton on: " << stats << std::endl;
430
431 if (!verbose) ::testing::internal::CaptureStdout();
432 cleanup_lammps(lmp, test_config);
433 lmp = init_lammps(argc, argv, test_config, false);
434 if (!verbose) ::testing::internal::GetCapturedStdout();
435
436 // skip over these tests if newton pair is forced to be on
437 if (lmp->force->newton_pair == 0) {
438
439 f = lmp->atom->f;
440 tag = lmp->atom->tag;
441 stats.reset();
442 for (int i = 0; i < nlocal; ++i) {
443 EXPECT_FP_LE_WITH_EPS(f[i][0], f_ref[tag[i]].x, epsilon);
444 EXPECT_FP_LE_WITH_EPS(f[i][1], f_ref[tag[i]].y, epsilon);
445 EXPECT_FP_LE_WITH_EPS(f[i][2], f_ref[tag[i]].z, epsilon);
446 }
447 if (print_stats) std::cerr << "init_forces stats, newton off:" << stats << std::endl;
448
449 pair = lmp->force->pair;
450 stress = pair->virial;
451 stats.reset();
452 EXPECT_FP_LE_WITH_EPS(stress[0], test_config.init_stress.xx, 3 * epsilon);
453 EXPECT_FP_LE_WITH_EPS(stress[1], test_config.init_stress.yy, 3 * epsilon);
454 EXPECT_FP_LE_WITH_EPS(stress[2], test_config.init_stress.zz, 3 * epsilon);
455 EXPECT_FP_LE_WITH_EPS(stress[3], test_config.init_stress.xy, 3 * epsilon);
456 EXPECT_FP_LE_WITH_EPS(stress[4], test_config.init_stress.xz, 3 * epsilon);
457 EXPECT_FP_LE_WITH_EPS(stress[5], test_config.init_stress.yz, 3 * epsilon);
458 if (print_stats) std::cerr << "init_stress stats, newton off:" << stats << std::endl;
459
460 stats.reset();
461 EXPECT_FP_LE_WITH_EPS(pair->eng_vdwl, test_config.init_vdwl, epsilon);
462 EXPECT_FP_LE_WITH_EPS(pair->eng_coul, test_config.init_coul, epsilon);
463 if (print_stats) std::cerr << "init_energy stats, newton off:" << stats << std::endl;
464
465 if (!verbose) ::testing::internal::CaptureStdout();
466 run_lammps(lmp);
467 if (!verbose) ::testing::internal::GetCapturedStdout();
468
469 f = lmp->atom->f;
470 tag = lmp->atom->tag;
471 stress = pair->virial;
472 stats.reset();
473 for (int i = 0; i < nlocal; ++i) {
474 EXPECT_FP_LE_WITH_EPS(f[i][0], f_run[tag[i]].x, 5 * epsilon);
475 EXPECT_FP_LE_WITH_EPS(f[i][1], f_run[tag[i]].y, 5 * epsilon);
476 EXPECT_FP_LE_WITH_EPS(f[i][2], f_run[tag[i]].z, 5 * epsilon);
477 }
478 if (print_stats) std::cerr << "run_forces stats, newton off:" << stats << std::endl;
479
480 stress = pair->virial;
481 stats.reset();
482 EXPECT_FP_LE_WITH_EPS(stress[0], test_config.run_stress.xx, epsilon);
483 EXPECT_FP_LE_WITH_EPS(stress[1], test_config.run_stress.yy, epsilon);
484 EXPECT_FP_LE_WITH_EPS(stress[2], test_config.run_stress.zz, epsilon);
485 EXPECT_FP_LE_WITH_EPS(stress[3], test_config.run_stress.xy, epsilon);
486 EXPECT_FP_LE_WITH_EPS(stress[4], test_config.run_stress.xz, epsilon);
487 EXPECT_FP_LE_WITH_EPS(stress[5], test_config.run_stress.yz, epsilon);
488 if (print_stats) std::cerr << "run_stress stats, newton off:" << stats << std::endl;
489
490 stats.reset();
491 id = lmp->modify->find_compute("sum");
492 energy = lmp->modify->compute[id]->compute_scalar();
493 EXPECT_FP_LE_WITH_EPS(pair->eng_vdwl, test_config.run_vdwl, epsilon);
494 EXPECT_FP_LE_WITH_EPS(pair->eng_coul, test_config.run_coul, epsilon);
495 // skip comparing per-atom energy with total energy for "kim"
496 if (std::string("kim") != lmp->force->pair_style)
497 EXPECT_FP_LE_WITH_EPS((pair->eng_vdwl + pair->eng_coul), energy, epsilon);
498 if (print_stats) std::cerr << "run_energy stats, newton off:" << stats << std::endl;
499 }
500
501 if (!verbose) ::testing::internal::CaptureStdout();
502 restart_lammps(lmp, test_config);
503 if (!verbose) ::testing::internal::GetCapturedStdout();
504
505 f = lmp->atom->f;
506 tag = lmp->atom->tag;
507 stats.reset();
508 ASSERT_EQ(nlocal + 1, f_ref.size());
509 for (int i = 0; i < nlocal; ++i) {
510 EXPECT_FP_LE_WITH_EPS(f[i][0], f_ref[tag[i]].x, epsilon);
511 EXPECT_FP_LE_WITH_EPS(f[i][1], f_ref[tag[i]].y, epsilon);
512 EXPECT_FP_LE_WITH_EPS(f[i][2], f_ref[tag[i]].z, epsilon);
513 }
514 if (print_stats) std::cerr << "restart_forces stats:" << stats << std::endl;
515
516 pair = lmp->force->pair;
517 stress = pair->virial;
518 stats.reset();
519 EXPECT_FP_LE_WITH_EPS(stress[0], test_config.init_stress.xx, epsilon);
520 EXPECT_FP_LE_WITH_EPS(stress[1], test_config.init_stress.yy, epsilon);
521 EXPECT_FP_LE_WITH_EPS(stress[2], test_config.init_stress.zz, epsilon);
522 EXPECT_FP_LE_WITH_EPS(stress[3], test_config.init_stress.xy, epsilon);
523 EXPECT_FP_LE_WITH_EPS(stress[4], test_config.init_stress.xz, epsilon);
524 EXPECT_FP_LE_WITH_EPS(stress[5], test_config.init_stress.yz, epsilon);
525 if (print_stats) std::cerr << "restart_stress stats:" << stats << std::endl;
526
527 stats.reset();
528 EXPECT_FP_LE_WITH_EPS(pair->eng_vdwl, test_config.init_vdwl, epsilon);
529 EXPECT_FP_LE_WITH_EPS(pair->eng_coul, test_config.init_coul, epsilon);
530 if (print_stats) std::cerr << "restart_energy stats:" << stats << std::endl;
531
532 // pair style rann does not support pair_modify nofdotr
533 if (test_config.pair_style != "rann") {
534 if (!verbose) ::testing::internal::CaptureStdout();
535 restart_lammps(lmp, test_config, true);
536 if (!verbose) ::testing::internal::GetCapturedStdout();
537
538 f = lmp->atom->f;
539 tag = lmp->atom->tag;
540 stats.reset();
541 ASSERT_EQ(nlocal + 1, f_ref.size());
542 for (int i = 0; i < nlocal; ++i) {
543 EXPECT_FP_LE_WITH_EPS(f[i][0], f_ref[tag[i]].x, epsilon);
544 EXPECT_FP_LE_WITH_EPS(f[i][1], f_ref[tag[i]].y, epsilon);
545 EXPECT_FP_LE_WITH_EPS(f[i][2], f_ref[tag[i]].z, epsilon);
546 }
547 if (print_stats) std::cerr << "nofdotr_forces stats:" << stats << std::endl;
548
549 pair = lmp->force->pair;
550 stress = pair->virial;
551 stats.reset();
552 EXPECT_FP_LE_WITH_EPS(stress[0], test_config.init_stress.xx, epsilon);
553 EXPECT_FP_LE_WITH_EPS(stress[1], test_config.init_stress.yy, epsilon);
554 EXPECT_FP_LE_WITH_EPS(stress[2], test_config.init_stress.zz, epsilon);
555 EXPECT_FP_LE_WITH_EPS(stress[3], test_config.init_stress.xy, epsilon);
556 EXPECT_FP_LE_WITH_EPS(stress[4], test_config.init_stress.xz, epsilon);
557 EXPECT_FP_LE_WITH_EPS(stress[5], test_config.init_stress.yz, epsilon);
558 if (print_stats) std::cerr << "nofdotr_stress stats:" << stats << std::endl;
559
560 stats.reset();
561 EXPECT_FP_LE_WITH_EPS(pair->eng_vdwl, test_config.init_vdwl, epsilon);
562 EXPECT_FP_LE_WITH_EPS(pair->eng_coul, test_config.init_coul, epsilon);
563 if (print_stats) std::cerr << "nofdotr_energy stats:" << stats << std::endl;
564 }
565
566 if (!verbose) ::testing::internal::CaptureStdout();
567 data_lammps(lmp, test_config);
568 if (!verbose) ::testing::internal::GetCapturedStdout();
569
570 f = lmp->atom->f;
571 tag = lmp->atom->tag;
572 stats.reset();
573 ASSERT_EQ(nlocal + 1, f_ref.size());
574 for (int i = 0; i < nlocal; ++i) {
575 EXPECT_FP_LE_WITH_EPS(f[i][0], f_ref[tag[i]].x, epsilon);
576 EXPECT_FP_LE_WITH_EPS(f[i][1], f_ref[tag[i]].y, epsilon);
577 EXPECT_FP_LE_WITH_EPS(f[i][2], f_ref[tag[i]].z, epsilon);
578 }
579 if (print_stats) std::cerr << "data_forces stats:" << stats << std::endl;
580
581 pair = lmp->force->pair;
582 stress = pair->virial;
583 stats.reset();
584 EXPECT_FP_LE_WITH_EPS(stress[0], test_config.init_stress.xx, epsilon);
585 EXPECT_FP_LE_WITH_EPS(stress[1], test_config.init_stress.yy, epsilon);
586 EXPECT_FP_LE_WITH_EPS(stress[2], test_config.init_stress.zz, epsilon);
587 EXPECT_FP_LE_WITH_EPS(stress[3], test_config.init_stress.xy, epsilon);
588 EXPECT_FP_LE_WITH_EPS(stress[4], test_config.init_stress.xz, epsilon);
589 EXPECT_FP_LE_WITH_EPS(stress[5], test_config.init_stress.yz, epsilon);
590 if (print_stats) std::cerr << "data_stress stats:" << stats << std::endl;
591
592 stats.reset();
593 EXPECT_FP_LE_WITH_EPS(pair->eng_vdwl, test_config.init_vdwl, epsilon);
594 EXPECT_FP_LE_WITH_EPS(pair->eng_coul, test_config.init_coul, epsilon);
595 if (print_stats) std::cerr << "data_energy stats:" << stats << std::endl;
596
597 if (pair->respa_enable) {
598 if (!verbose) ::testing::internal::CaptureStdout();
599 cleanup_lammps(lmp, test_config);
600 lmp = init_lammps(argc, argv, test_config, false);
601 lmp->input->one("run_style respa 2 1 inner 1 4.8 5.5 outer 2");
602 run_lammps(lmp);
603 if (!verbose) ::testing::internal::GetCapturedStdout();
604
605 // need to relax error by a large amount with tabulation, since
606 // coul/long styles do not use tabulation in compute_inner()
607 // and compute_middle() so we get a significant deviation.
608 pair = lmp->force->pair;
609 if (pair->ncoultablebits) epsilon *= 5.0e6;
610
611 f = lmp->atom->f;
612 tag = lmp->atom->tag;
613 stats.reset();
614 for (int i = 0; i < nlocal; ++i) {
615 EXPECT_FP_LE_WITH_EPS(f[i][0], f_run[tag[i]].x, 5 * epsilon);
616 EXPECT_FP_LE_WITH_EPS(f[i][1], f_run[tag[i]].y, 5 * epsilon);
617 EXPECT_FP_LE_WITH_EPS(f[i][2], f_run[tag[i]].z, 5 * epsilon);
618 }
619 if (print_stats) std::cerr << "run_forces stats, r-RESPA:" << stats << std::endl;
620
621 stress = pair->virial;
622 stats.reset();
623 EXPECT_FP_LE_WITH_EPS(stress[0], test_config.run_stress.xx, epsilon);
624 EXPECT_FP_LE_WITH_EPS(stress[1], test_config.run_stress.yy, epsilon);
625 EXPECT_FP_LE_WITH_EPS(stress[2], test_config.run_stress.zz, epsilon);
626 EXPECT_FP_LE_WITH_EPS(stress[3], test_config.run_stress.xy, epsilon);
627 EXPECT_FP_LE_WITH_EPS(stress[4], test_config.run_stress.xz, epsilon);
628 EXPECT_FP_LE_WITH_EPS(stress[5], test_config.run_stress.yz, epsilon);
629 if (print_stats) std::cerr << "run_stress stats, r-RESPA:" << stats << std::endl;
630
631 stats.reset();
632 id = lmp->modify->find_compute("sum");
633 energy = lmp->modify->compute[id]->compute_scalar();
634 EXPECT_FP_LE_WITH_EPS(pair->eng_vdwl, test_config.run_vdwl, epsilon);
635 EXPECT_FP_LE_WITH_EPS(pair->eng_coul, test_config.run_coul, epsilon);
636 EXPECT_FP_LE_WITH_EPS((pair->eng_vdwl + pair->eng_coul), energy, epsilon);
637 if (print_stats) std::cerr << "run_energy stats, r-RESPA:" << stats << std::endl;
638 }
639 if (!verbose) ::testing::internal::CaptureStdout();
640 cleanup_lammps(lmp, test_config);
641 if (!verbose) ::testing::internal::GetCapturedStdout();
642 };
643
TEST(PairStyle,omp)644 TEST(PairStyle, omp)
645 {
646 if (!LAMMPS::is_installed_pkg("OPENMP")) GTEST_SKIP();
647 if (test_config.skip_tests.count(test_info_->name())) GTEST_SKIP();
648
649 const char *args[] = {"PairStyle", "-log", "none", "-echo", "screen", "-nocite",
650 "-pk", "omp", "4", "-sf", "omp"};
651
652 // cannot run dpd styles with more than 1 thread due to using multiple pRNGs
653 if (utils::strmatch(test_config.pair_style, "^dpd")) args[8] = "1";
654
655 char **argv = (char **)args;
656 int argc = sizeof(args) / sizeof(char *);
657
658 ::testing::internal::CaptureStdout();
659 LAMMPS *lmp = init_lammps(argc, argv, test_config, true);
660
661 std::string output = ::testing::internal::GetCapturedStdout();
662 if (verbose) std::cout << output;
663
664 if (!lmp) {
665 std::cerr << "One or more prerequisite styles with /omp suffix\n"
666 "are not available in this LAMMPS configuration:\n";
667 for (auto &prerequisite : test_config.prerequisites) {
668 std::cerr << prerequisite.first << "_style " << prerequisite.second << "\n";
669 }
670 GTEST_SKIP();
671 }
672
673 EXPECT_THAT(output, StartsWith("LAMMPS ("));
674 EXPECT_THAT(output, HasSubstr("Loop time"));
675
676 // abort if running in parallel and not all atoms are local
677 const int nlocal = lmp->atom->nlocal;
678 ASSERT_EQ(lmp->atom->natoms, nlocal);
679
680 // relax error a bit for OPENMP package
681 double epsilon = 5.0 * test_config.epsilon;
682 // relax test precision when using pppm and single precision FFTs
683 #if defined(FFT_SINGLE)
684 if (lmp->force->kspace && lmp->force->kspace->compute_flag)
685 if (utils::strmatch(lmp->force->kspace_style, "^pppm")) epsilon *= 2.0e8;
686 #endif
687 auto f = lmp->atom->f;
688 auto tag = lmp->atom->tag;
689
690 const std::vector<coord_t> &f_ref = test_config.init_forces;
691 ErrorStats stats;
692 stats.reset();
693 for (int i = 0; i < nlocal; ++i) {
694 EXPECT_FP_LE_WITH_EPS(f[i][0], f_ref[tag[i]].x, epsilon);
695 EXPECT_FP_LE_WITH_EPS(f[i][1], f_ref[tag[i]].y, epsilon);
696 EXPECT_FP_LE_WITH_EPS(f[i][2], f_ref[tag[i]].z, epsilon);
697 }
698 if (print_stats) std::cerr << "init_forces stats, newton on: " << stats << std::endl;
699
700 auto pair = lmp->force->pair;
701 auto stress = pair->virial;
702 stats.reset();
703 EXPECT_FP_LE_WITH_EPS(stress[0], test_config.init_stress.xx, 10 * epsilon);
704 EXPECT_FP_LE_WITH_EPS(stress[1], test_config.init_stress.yy, 10 * epsilon);
705 EXPECT_FP_LE_WITH_EPS(stress[2], test_config.init_stress.zz, 10 * epsilon);
706 EXPECT_FP_LE_WITH_EPS(stress[3], test_config.init_stress.xy, 10 * epsilon);
707 EXPECT_FP_LE_WITH_EPS(stress[4], test_config.init_stress.xz, 10 * epsilon);
708 EXPECT_FP_LE_WITH_EPS(stress[5], test_config.init_stress.yz, 10 * epsilon);
709 if (print_stats) std::cerr << "init_stress stats, newton on: " << stats << std::endl;
710
711 stats.reset();
712 EXPECT_FP_LE_WITH_EPS(pair->eng_vdwl, test_config.init_vdwl, epsilon);
713 EXPECT_FP_LE_WITH_EPS(pair->eng_coul, test_config.init_coul, epsilon);
714 if (print_stats) std::cerr << "init_energy stats, newton on: " << stats << std::endl;
715
716 if (!verbose) ::testing::internal::CaptureStdout();
717 run_lammps(lmp);
718 if (!verbose) ::testing::internal::GetCapturedStdout();
719
720 f = lmp->atom->f;
721 stress = pair->virial;
722 tag = lmp->atom->tag;
723
724 const std::vector<coord_t> &f_run = test_config.run_forces;
725 ASSERT_EQ(nlocal + 1, f_run.size());
726 stats.reset();
727 for (int i = 0; i < nlocal; ++i) {
728 EXPECT_FP_LE_WITH_EPS(f[i][0], f_run[tag[i]].x, 5 * epsilon);
729 EXPECT_FP_LE_WITH_EPS(f[i][1], f_run[tag[i]].y, 5 * epsilon);
730 EXPECT_FP_LE_WITH_EPS(f[i][2], f_run[tag[i]].z, 5 * epsilon);
731 }
732 if (print_stats) std::cerr << "run_forces stats, newton on: " << stats << std::endl;
733
734 stress = pair->virial;
735 stats.reset();
736 EXPECT_FP_LE_WITH_EPS(stress[0], test_config.run_stress.xx, 10 * epsilon);
737 EXPECT_FP_LE_WITH_EPS(stress[1], test_config.run_stress.yy, 10 * epsilon);
738 EXPECT_FP_LE_WITH_EPS(stress[2], test_config.run_stress.zz, 10 * epsilon);
739 EXPECT_FP_LE_WITH_EPS(stress[3], test_config.run_stress.xy, 10 * epsilon);
740 EXPECT_FP_LE_WITH_EPS(stress[4], test_config.run_stress.xz, 10 * epsilon);
741 EXPECT_FP_LE_WITH_EPS(stress[5], test_config.run_stress.yz, 10 * epsilon);
742 if (print_stats) std::cerr << "run_stress stats, newton on: " << stats << std::endl;
743
744 stats.reset();
745 int id = lmp->modify->find_compute("sum");
746 double energy = lmp->modify->compute[id]->compute_scalar();
747 EXPECT_FP_LE_WITH_EPS(pair->eng_vdwl, test_config.run_vdwl, epsilon);
748 EXPECT_FP_LE_WITH_EPS(pair->eng_coul, test_config.run_coul, epsilon);
749 EXPECT_FP_LE_WITH_EPS((pair->eng_vdwl + pair->eng_coul), energy, epsilon);
750 if (print_stats) std::cerr << "run_energy stats, newton on: " << stats << std::endl;
751
752 // skip over these tests if newton pair is forced to be on
753 if (lmp->force->newton_pair == 0) {
754
755 if (!verbose) ::testing::internal::CaptureStdout();
756 cleanup_lammps(lmp, test_config);
757 lmp = init_lammps(argc, argv, test_config, false);
758 if (!verbose) ::testing::internal::GetCapturedStdout();
759
760 f = lmp->atom->f;
761 tag = lmp->atom->tag;
762 stats.reset();
763 for (int i = 0; i < nlocal; ++i) {
764 EXPECT_FP_LE_WITH_EPS(f[i][0], f_ref[tag[i]].x, epsilon);
765 EXPECT_FP_LE_WITH_EPS(f[i][1], f_ref[tag[i]].y, epsilon);
766 EXPECT_FP_LE_WITH_EPS(f[i][2], f_ref[tag[i]].z, epsilon);
767 }
768 if (print_stats) std::cerr << "init_forces stats, newton off:" << stats << std::endl;
769
770 pair = lmp->force->pair;
771 stress = pair->virial;
772 stats.reset();
773 EXPECT_FP_LE_WITH_EPS(stress[0], test_config.init_stress.xx, 10 * epsilon);
774 EXPECT_FP_LE_WITH_EPS(stress[1], test_config.init_stress.yy, 10 * epsilon);
775 EXPECT_FP_LE_WITH_EPS(stress[2], test_config.init_stress.zz, 10 * epsilon);
776 EXPECT_FP_LE_WITH_EPS(stress[3], test_config.init_stress.xy, 10 * epsilon);
777 EXPECT_FP_LE_WITH_EPS(stress[4], test_config.init_stress.xz, 10 * epsilon);
778 EXPECT_FP_LE_WITH_EPS(stress[5], test_config.init_stress.yz, 10 * epsilon);
779 if (print_stats) std::cerr << "init_stress stats, newton off:" << stats << std::endl;
780
781 stats.reset();
782 EXPECT_FP_LE_WITH_EPS(pair->eng_vdwl, test_config.init_vdwl, epsilon);
783 EXPECT_FP_LE_WITH_EPS(pair->eng_coul, test_config.init_coul, epsilon);
784 if (print_stats) std::cerr << "init_energy stats, newton off:" << stats << std::endl;
785
786 if (!verbose) ::testing::internal::CaptureStdout();
787 run_lammps(lmp);
788 if (!verbose) ::testing::internal::GetCapturedStdout();
789
790 f = lmp->atom->f;
791 tag = lmp->atom->tag;
792 stats.reset();
793 for (int i = 0; i < nlocal; ++i) {
794 EXPECT_FP_LE_WITH_EPS(f[i][0], f_run[tag[i]].x, 5 * epsilon);
795 EXPECT_FP_LE_WITH_EPS(f[i][1], f_run[tag[i]].y, 5 * epsilon);
796 EXPECT_FP_LE_WITH_EPS(f[i][2], f_run[tag[i]].z, 5 * epsilon);
797 }
798 if (print_stats) std::cerr << "run_forces stats, newton off:" << stats << std::endl;
799
800 stress = pair->virial;
801 stats.reset();
802 EXPECT_FP_LE_WITH_EPS(stress[0], test_config.run_stress.xx, 10 * epsilon);
803 EXPECT_FP_LE_WITH_EPS(stress[1], test_config.run_stress.yy, 10 * epsilon);
804 EXPECT_FP_LE_WITH_EPS(stress[2], test_config.run_stress.zz, 10 * epsilon);
805 EXPECT_FP_LE_WITH_EPS(stress[3], test_config.run_stress.xy, 10 * epsilon);
806 EXPECT_FP_LE_WITH_EPS(stress[4], test_config.run_stress.xz, 10 * epsilon);
807 EXPECT_FP_LE_WITH_EPS(stress[5], test_config.run_stress.yz, 10 * epsilon);
808 if (print_stats) std::cerr << "run_stress stats, newton off:" << stats << std::endl;
809
810 stats.reset();
811 id = lmp->modify->find_compute("sum");
812 energy = lmp->modify->compute[id]->compute_scalar();
813 EXPECT_FP_LE_WITH_EPS(pair->eng_vdwl, test_config.run_vdwl, epsilon);
814 EXPECT_FP_LE_WITH_EPS(pair->eng_coul, test_config.run_coul, epsilon);
815 EXPECT_FP_LE_WITH_EPS((pair->eng_vdwl + pair->eng_coul), energy, epsilon);
816 if (print_stats) std::cerr << "run_energy stats, newton off:" << stats << std::endl;
817 }
818
819 if (!verbose) ::testing::internal::CaptureStdout();
820 restart_lammps(lmp, test_config, true);
821 if (!verbose) ::testing::internal::GetCapturedStdout();
822
823 f = lmp->atom->f;
824 tag = lmp->atom->tag;
825 stats.reset();
826 ASSERT_EQ(nlocal + 1, f_ref.size());
827 for (int i = 0; i < nlocal; ++i) {
828 EXPECT_FP_LE_WITH_EPS(f[i][0], f_ref[tag[i]].x, 5 * epsilon);
829 EXPECT_FP_LE_WITH_EPS(f[i][1], f_ref[tag[i]].y, 5 * epsilon);
830 EXPECT_FP_LE_WITH_EPS(f[i][2], f_ref[tag[i]].z, 5 * epsilon);
831 }
832 if (print_stats) std::cerr << "nofdotr_forces stats:" << stats << std::endl;
833
834 pair = lmp->force->pair;
835 stress = pair->virial;
836 stats.reset();
837 EXPECT_FP_LE_WITH_EPS(stress[0], test_config.init_stress.xx, 10 * epsilon);
838 EXPECT_FP_LE_WITH_EPS(stress[1], test_config.init_stress.yy, 10 * epsilon);
839 EXPECT_FP_LE_WITH_EPS(stress[2], test_config.init_stress.zz, 10 * epsilon);
840 EXPECT_FP_LE_WITH_EPS(stress[3], test_config.init_stress.xy, 10 * epsilon);
841 EXPECT_FP_LE_WITH_EPS(stress[4], test_config.init_stress.xz, 10 * epsilon);
842 EXPECT_FP_LE_WITH_EPS(stress[5], test_config.init_stress.yz, 10 * epsilon);
843 if (print_stats) std::cerr << "nofdotr_stress stats:" << stats << std::endl;
844
845 stats.reset();
846 EXPECT_FP_LE_WITH_EPS(pair->eng_vdwl, test_config.init_vdwl, 5 * epsilon);
847 EXPECT_FP_LE_WITH_EPS(pair->eng_coul, test_config.init_coul, 5 * epsilon);
848 if (print_stats) std::cerr << "nofdotr_energy stats:" << stats << std::endl;
849
850 if (!verbose) ::testing::internal::CaptureStdout();
851 cleanup_lammps(lmp, test_config);
852 if (!verbose) ::testing::internal::GetCapturedStdout();
853 };
854
TEST(PairStyle,gpu)855 TEST(PairStyle, gpu)
856 {
857 if (!LAMMPS::is_installed_pkg("GPU")) GTEST_SKIP();
858 if (!Info::has_gpu_device()) GTEST_SKIP();
859 if (test_config.skip_tests.count(test_info_->name())) GTEST_SKIP();
860
861 const char *args_neigh[] = {"PairStyle", "-log", "none", "-echo",
862 "screen", "-nocite", "-sf", "gpu"};
863 const char *args_noneigh[] = {"PairStyle", "-log", "none", "-echo", "screen", "-nocite", "-sf",
864 "gpu", "-pk", "gpu", "0", "neigh", "no"};
865
866 char **argv = (char **)args_neigh;
867 int argc = sizeof(args_neigh) / sizeof(char *);
868
869 // cannot use GPU neighbor list with hybrid pair style (yet)
870 if (test_config.pair_style.substr(0, 6) == "hybrid") {
871 argv = (char **)args_noneigh;
872 argc = sizeof(args_noneigh) / sizeof(char *);
873 }
874
875 ::testing::internal::CaptureStdout();
876 LAMMPS *lmp = init_lammps(argc, argv, test_config, false);
877
878 std::string output = ::testing::internal::GetCapturedStdout();
879 if (verbose) std::cout << output;
880
881 if (!lmp) {
882 std::cerr << "One or more prerequisite styles with /gpu suffix\n"
883 "are not available in this LAMMPS configuration:\n";
884 for (auto &prerequisite : test_config.prerequisites) {
885 std::cerr << prerequisite.first << "_style " << prerequisite.second << "\n";
886 }
887 GTEST_SKIP();
888 }
889
890 EXPECT_THAT(output, StartsWith("LAMMPS ("));
891 EXPECT_THAT(output, HasSubstr("Loop time"));
892
893 // abort if running in parallel and not all atoms are local
894 const int nlocal = lmp->atom->nlocal;
895 ASSERT_EQ(lmp->atom->natoms, nlocal);
896
897 // relax error for GPU package depending on precision setting
898 double epsilon = test_config.epsilon;
899 if (Info::has_accelerator_feature("GPU", "precision", "double"))
900 epsilon *= 7.5;
901 else if (Info::has_accelerator_feature("GPU", "precision", "mixed"))
902 epsilon *= 5.0e8;
903 else
904 epsilon *= 1.0e10;
905 // relax test precision when using pppm and single precision FFTs, but only when also
906 // running with double precision
907 #if defined(FFT_SINGLE)
908 if (lmp->force->kspace && lmp->force->kspace->compute_flag &&
909 Info::has_accelerator_feature("GPU", "precision", "double"))
910 if (utils::strmatch(lmp->force->kspace_style, "^pppm")) epsilon *= 2.0e8;
911 #endif
912 const std::vector<coord_t> &f_ref = test_config.init_forces;
913 const std::vector<coord_t> &f_run = test_config.run_forces;
914 ErrorStats stats;
915
916 auto f = lmp->atom->f;
917 auto tag = lmp->atom->tag;
918 stats.reset();
919 for (int i = 0; i < nlocal; ++i) {
920 EXPECT_FP_LE_WITH_EPS(f[i][0], f_ref[tag[i]].x, epsilon);
921 EXPECT_FP_LE_WITH_EPS(f[i][1], f_ref[tag[i]].y, epsilon);
922 EXPECT_FP_LE_WITH_EPS(f[i][2], f_ref[tag[i]].z, epsilon);
923 }
924 if (print_stats) std::cerr << "init_forces stats, newton off:" << stats << std::endl;
925
926 auto pair = lmp->force->pair;
927 auto stress = pair->virial;
928 stats.reset();
929 EXPECT_FP_LE_WITH_EPS(stress[0], test_config.init_stress.xx, 10 * epsilon);
930 EXPECT_FP_LE_WITH_EPS(stress[1], test_config.init_stress.yy, 10 * epsilon);
931 EXPECT_FP_LE_WITH_EPS(stress[2], test_config.init_stress.zz, 10 * epsilon);
932 EXPECT_FP_LE_WITH_EPS(stress[3], test_config.init_stress.xy, 10 * epsilon);
933 EXPECT_FP_LE_WITH_EPS(stress[4], test_config.init_stress.xz, 10 * epsilon);
934 EXPECT_FP_LE_WITH_EPS(stress[5], test_config.init_stress.yz, 10 * epsilon);
935 if (print_stats) std::cerr << "init_stress stats, newton off:" << stats << std::endl;
936
937 stats.reset();
938 EXPECT_FP_LE_WITH_EPS(pair->eng_vdwl, test_config.init_vdwl, epsilon);
939 EXPECT_FP_LE_WITH_EPS(pair->eng_coul, test_config.init_coul, epsilon);
940 if (print_stats) std::cerr << "init_energy stats, newton off:" << stats << std::endl;
941
942 if (!verbose) ::testing::internal::CaptureStdout();
943 run_lammps(lmp);
944 if (!verbose) ::testing::internal::GetCapturedStdout();
945
946 f = lmp->atom->f;
947 tag = lmp->atom->tag;
948 stats.reset();
949 for (int i = 0; i < nlocal; ++i) {
950 EXPECT_FP_LE_WITH_EPS(f[i][0], f_run[tag[i]].x, 5 * epsilon);
951 EXPECT_FP_LE_WITH_EPS(f[i][1], f_run[tag[i]].y, 5 * epsilon);
952 EXPECT_FP_LE_WITH_EPS(f[i][2], f_run[tag[i]].z, 5 * epsilon);
953 }
954 if (print_stats) std::cerr << "run_forces stats, newton off:" << stats << std::endl;
955
956 stress = pair->virial;
957 stats.reset();
958 EXPECT_FP_LE_WITH_EPS(stress[0], test_config.run_stress.xx, 10 * epsilon);
959 EXPECT_FP_LE_WITH_EPS(stress[1], test_config.run_stress.yy, 10 * epsilon);
960 EXPECT_FP_LE_WITH_EPS(stress[2], test_config.run_stress.zz, 10 * epsilon);
961 EXPECT_FP_LE_WITH_EPS(stress[3], test_config.run_stress.xy, 10 * epsilon);
962 EXPECT_FP_LE_WITH_EPS(stress[4], test_config.run_stress.xz, 10 * epsilon);
963 EXPECT_FP_LE_WITH_EPS(stress[5], test_config.run_stress.yz, 10 * epsilon);
964 if (print_stats) std::cerr << "run_stress stats, newton off:" << stats << std::endl;
965
966 stats.reset();
967 auto id = lmp->modify->find_compute("sum");
968 auto energy = lmp->modify->compute[id]->compute_scalar();
969 EXPECT_FP_LE_WITH_EPS(pair->eng_vdwl, test_config.run_vdwl, epsilon);
970 EXPECT_FP_LE_WITH_EPS(pair->eng_coul, test_config.run_coul, epsilon);
971 EXPECT_FP_LE_WITH_EPS((pair->eng_vdwl + pair->eng_coul), energy, epsilon);
972 if (print_stats) std::cerr << "run_energy stats, newton off:" << stats << std::endl;
973
974 if (!verbose) ::testing::internal::CaptureStdout();
975 cleanup_lammps(lmp, test_config);
976 if (!verbose) ::testing::internal::GetCapturedStdout();
977 };
978
TEST(PairStyle,intel)979 TEST(PairStyle, intel)
980 {
981 if (!LAMMPS::is_installed_pkg("INTEL")) GTEST_SKIP();
982 if (test_config.skip_tests.count(test_info_->name())) GTEST_SKIP();
983
984 const char *args[] = {"PairStyle", "-log", "none", "-echo", "screen", "-nocite",
985 "-pk", "intel", "0", "mode", "double", "omp",
986 "4", "lrt", "no", "-sf", "intel"};
987
988 // cannot use more than 1 thread for dpd styles due to pRNG
989 if (utils::strmatch(test_config.pair_style, "^dpd")) args[12] = "1";
990
991 char **argv = (char **)args;
992 int argc = sizeof(args) / sizeof(char *);
993
994 ::testing::internal::CaptureStdout();
995 LAMMPS *lmp = init_lammps(argc, argv, test_config, true);
996
997 std::string output = ::testing::internal::GetCapturedStdout();
998 if (verbose) std::cout << output;
999
1000 if (!lmp) {
1001 std::cerr << "One or more prerequisite styles with /intel suffix\n"
1002 "are not available in this LAMMPS configuration:\n";
1003 for (auto prerequisite : test_config.prerequisites) {
1004 std::cerr << prerequisite.first << "_style " << prerequisite.second << "\n";
1005 }
1006 GTEST_SKIP();
1007 }
1008
1009 // relax error a bit for INTEL package
1010 double epsilon = 7.5 * test_config.epsilon;
1011 // relax test precision when using pppm and single precision FFTs
1012 #if defined(FFT_SINGLE)
1013 if (lmp->force->kspace && lmp->force->kspace->compute_flag)
1014 if (utils::strmatch(lmp->force->kspace_style, "^pppm")) epsilon *= 2.0e8;
1015 #endif
1016
1017 // we need to relax the epsilon a LOT for tests using long-range
1018 // coulomb with tabulation. seems more like mixed precision or a bug
1019 for (auto post_cmd : test_config.post_commands) {
1020 if (post_cmd.find("pair_modify table") != std::string::npos) {
1021 if (post_cmd.find("pair_modify table 0") == std::string::npos) epsilon *= 1000000.0;
1022 }
1023 }
1024
1025 EXPECT_THAT(output, StartsWith("LAMMPS ("));
1026 EXPECT_THAT(output, HasSubstr("Loop time"));
1027
1028 // abort if running in parallel and not all atoms are local
1029 const int nlocal = lmp->atom->nlocal;
1030 ASSERT_EQ(lmp->atom->natoms, nlocal);
1031
1032 auto f = lmp->atom->f;
1033 auto tag = lmp->atom->tag;
1034
1035 const std::vector<coord_t> &f_ref = test_config.init_forces;
1036 ErrorStats stats;
1037 stats.reset();
1038 for (int i = 0; i < nlocal; ++i) {
1039 EXPECT_FP_LE_WITH_EPS(f[i][0], f_ref[tag[i]].x, epsilon);
1040 EXPECT_FP_LE_WITH_EPS(f[i][1], f_ref[tag[i]].y, epsilon);
1041 EXPECT_FP_LE_WITH_EPS(f[i][2], f_ref[tag[i]].z, epsilon);
1042 }
1043 if (print_stats) std::cerr << "init_forces stats:" << stats << std::endl;
1044
1045 Pair *pair = lmp->force->pair;
1046 double *stress = pair->virial;
1047 stats.reset();
1048 EXPECT_FP_LE_WITH_EPS(stress[0], test_config.init_stress.xx, 10 * epsilon);
1049 EXPECT_FP_LE_WITH_EPS(stress[1], test_config.init_stress.yy, 10 * epsilon);
1050 EXPECT_FP_LE_WITH_EPS(stress[2], test_config.init_stress.zz, 10 * epsilon);
1051 EXPECT_FP_LE_WITH_EPS(stress[3], test_config.init_stress.xy, 10 * epsilon);
1052 EXPECT_FP_LE_WITH_EPS(stress[4], test_config.init_stress.xz, 10 * epsilon);
1053 EXPECT_FP_LE_WITH_EPS(stress[5], test_config.init_stress.yz, 10 * epsilon);
1054 if (print_stats) std::cerr << "init_stress stats:" << stats << std::endl;
1055
1056 stats.reset();
1057 EXPECT_FP_LE_WITH_EPS(pair->eng_vdwl, test_config.init_vdwl, epsilon);
1058 EXPECT_FP_LE_WITH_EPS(pair->eng_coul, test_config.init_coul, epsilon);
1059 if (print_stats) std::cerr << "init_energy stats:" << stats << std::endl;
1060
1061 if (!verbose) ::testing::internal::CaptureStdout();
1062 run_lammps(lmp);
1063 if (!verbose) ::testing::internal::GetCapturedStdout();
1064
1065 f = lmp->atom->f;
1066 tag = lmp->atom->tag;
1067 stress = pair->virial;
1068
1069 const std::vector<coord_t> &f_run = test_config.run_forces;
1070 ASSERT_EQ(nlocal + 1, f_run.size());
1071 stats.reset();
1072 for (int i = 0; i < nlocal; ++i) {
1073 EXPECT_FP_LE_WITH_EPS(f[i][0], f_run[tag[i]].x, 5 * epsilon);
1074 EXPECT_FP_LE_WITH_EPS(f[i][1], f_run[tag[i]].y, 5 * epsilon);
1075 EXPECT_FP_LE_WITH_EPS(f[i][2], f_run[tag[i]].z, 5 * epsilon);
1076 }
1077 if (print_stats) std::cerr << "run_forces stats:" << stats << std::endl;
1078
1079 stress = pair->virial;
1080 stats.reset();
1081 EXPECT_FP_LE_WITH_EPS(stress[0], test_config.run_stress.xx, 10 * epsilon);
1082 EXPECT_FP_LE_WITH_EPS(stress[1], test_config.run_stress.yy, 10 * epsilon);
1083 EXPECT_FP_LE_WITH_EPS(stress[2], test_config.run_stress.zz, 10 * epsilon);
1084 EXPECT_FP_LE_WITH_EPS(stress[3], test_config.run_stress.xy, 10 * epsilon);
1085 EXPECT_FP_LE_WITH_EPS(stress[4], test_config.run_stress.xz, 10 * epsilon);
1086 EXPECT_FP_LE_WITH_EPS(stress[5], test_config.run_stress.yz, 10 * epsilon);
1087 if (print_stats) std::cerr << "run_stress stats:" << stats << std::endl;
1088
1089 stats.reset();
1090 int id = lmp->modify->find_compute("sum");
1091 double energy = lmp->modify->compute[id]->compute_scalar();
1092 EXPECT_FP_LE_WITH_EPS(pair->eng_vdwl, test_config.run_vdwl, epsilon);
1093 EXPECT_FP_LE_WITH_EPS(pair->eng_coul, test_config.run_coul, epsilon);
1094
1095 // rebo family of pair styles will have a large error in per-atom energy for INTEL
1096 if (test_config.pair_style.find("rebo") != std::string::npos) epsilon *= 100000.0;
1097
1098 EXPECT_FP_LE_WITH_EPS((pair->eng_vdwl + pair->eng_coul), energy, epsilon);
1099 if (print_stats) std::cerr << "run_energy stats:" << stats << std::endl;
1100
1101 if (!verbose) ::testing::internal::CaptureStdout();
1102 cleanup_lammps(lmp, test_config);
1103 if (!verbose) ::testing::internal::GetCapturedStdout();
1104 };
1105
TEST(PairStyle,opt)1106 TEST(PairStyle, opt)
1107 {
1108 if (!LAMMPS::is_installed_pkg("OPT")) GTEST_SKIP();
1109 if (test_config.skip_tests.count(test_info_->name())) GTEST_SKIP();
1110
1111 const char *args[] = {"PairStyle", "-log", "none", "-echo", "screen", "-nocite", "-sf", "opt"};
1112
1113 char **argv = (char **)args;
1114 int argc = sizeof(args) / sizeof(char *);
1115
1116 ::testing::internal::CaptureStdout();
1117 LAMMPS *lmp = init_lammps(argc, argv, test_config);
1118
1119 std::string output = ::testing::internal::GetCapturedStdout();
1120 if (verbose) std::cout << output;
1121
1122 if (!lmp) {
1123 std::cerr << "One or more prerequisite styles with /opt suffix\n"
1124 "are not available in this LAMMPS configuration:\n";
1125 for (auto prerequisite : test_config.prerequisites) {
1126 std::cerr << prerequisite.first << "_style " << prerequisite.second << "\n";
1127 }
1128 GTEST_SKIP();
1129 }
1130
1131 EXPECT_THAT(output, StartsWith("LAMMPS ("));
1132 EXPECT_THAT(output, HasSubstr("Loop time"));
1133
1134 // abort if running in parallel and not all atoms are local
1135 const int nlocal = lmp->atom->nlocal;
1136 ASSERT_EQ(lmp->atom->natoms, nlocal);
1137
1138 // relax error a bit for OPT package
1139 double epsilon = 2.0 * test_config.epsilon;
1140 // relax test precision when using pppm and single precision FFTs
1141 #if defined(FFT_SINGLE)
1142 if (lmp->force->kspace && lmp->force->kspace->compute_flag)
1143 if (utils::strmatch(lmp->force->kspace_style, "^pppm")) epsilon *= 2.0e8;
1144 #endif
1145 auto f = lmp->atom->f;
1146 auto tag = lmp->atom->tag;
1147
1148 const std::vector<coord_t> &f_ref = test_config.init_forces;
1149 ErrorStats stats;
1150 stats.reset();
1151 for (int i = 0; i < nlocal; ++i) {
1152 EXPECT_FP_LE_WITH_EPS(f[i][0], f_ref[tag[i]].x, epsilon);
1153 EXPECT_FP_LE_WITH_EPS(f[i][1], f_ref[tag[i]].y, epsilon);
1154 EXPECT_FP_LE_WITH_EPS(f[i][2], f_ref[tag[i]].z, epsilon);
1155 }
1156 if (print_stats) std::cerr << "init_forces stats:" << stats << std::endl;
1157
1158 Pair *pair = lmp->force->pair;
1159 double *stress = pair->virial;
1160 stats.reset();
1161 EXPECT_FP_LE_WITH_EPS(stress[0], test_config.init_stress.xx, 10 * epsilon);
1162 EXPECT_FP_LE_WITH_EPS(stress[1], test_config.init_stress.yy, 10 * epsilon);
1163 EXPECT_FP_LE_WITH_EPS(stress[2], test_config.init_stress.zz, 10 * epsilon);
1164 EXPECT_FP_LE_WITH_EPS(stress[3], test_config.init_stress.xy, 10 * epsilon);
1165 EXPECT_FP_LE_WITH_EPS(stress[4], test_config.init_stress.xz, 10 * epsilon);
1166 EXPECT_FP_LE_WITH_EPS(stress[5], test_config.init_stress.yz, 10 * epsilon);
1167 if (print_stats) std::cerr << "init_stress stats:" << stats << std::endl;
1168
1169 stats.reset();
1170 EXPECT_FP_LE_WITH_EPS(pair->eng_vdwl, test_config.init_vdwl, epsilon);
1171 EXPECT_FP_LE_WITH_EPS(pair->eng_coul, test_config.init_coul, epsilon);
1172 if (print_stats) std::cerr << "init_energy stats:" << stats << std::endl;
1173
1174 if (!verbose) ::testing::internal::CaptureStdout();
1175 run_lammps(lmp);
1176 if (!verbose) ::testing::internal::GetCapturedStdout();
1177
1178 f = lmp->atom->f;
1179 tag = lmp->atom->tag;
1180 stress = pair->virial;
1181 const std::vector<coord_t> &f_run = test_config.run_forces;
1182 ASSERT_EQ(nlocal + 1, f_run.size());
1183 stats.reset();
1184 for (int i = 0; i < nlocal; ++i) {
1185 EXPECT_FP_LE_WITH_EPS(f[i][0], f_run[tag[i]].x, 5 * epsilon);
1186 EXPECT_FP_LE_WITH_EPS(f[i][1], f_run[tag[i]].y, 5 * epsilon);
1187 EXPECT_FP_LE_WITH_EPS(f[i][2], f_run[tag[i]].z, 5 * epsilon);
1188 }
1189 if (print_stats) std::cerr << "run_forces stats:" << stats << std::endl;
1190
1191 stress = pair->virial;
1192 stats.reset();
1193 EXPECT_FP_LE_WITH_EPS(stress[0], test_config.run_stress.xx, 10 * epsilon);
1194 EXPECT_FP_LE_WITH_EPS(stress[1], test_config.run_stress.yy, 10 * epsilon);
1195 EXPECT_FP_LE_WITH_EPS(stress[2], test_config.run_stress.zz, 10 * epsilon);
1196 EXPECT_FP_LE_WITH_EPS(stress[3], test_config.run_stress.xy, 10 * epsilon);
1197 EXPECT_FP_LE_WITH_EPS(stress[4], test_config.run_stress.xz, 10 * epsilon);
1198 EXPECT_FP_LE_WITH_EPS(stress[5], test_config.run_stress.yz, 10 * epsilon);
1199 if (print_stats) std::cerr << "run_stress stats:" << stats << std::endl;
1200
1201 stats.reset();
1202 int id = lmp->modify->find_compute("sum");
1203 double energy = lmp->modify->compute[id]->compute_scalar();
1204 EXPECT_FP_LE_WITH_EPS(pair->eng_vdwl, test_config.run_vdwl, epsilon);
1205 EXPECT_FP_LE_WITH_EPS(pair->eng_coul, test_config.run_coul, epsilon);
1206 EXPECT_FP_LE_WITH_EPS((pair->eng_vdwl + pair->eng_coul), energy, epsilon);
1207 if (print_stats) std::cerr << "run_energy stats:" << stats << std::endl;
1208
1209 if (!verbose) ::testing::internal::CaptureStdout();
1210 restart_lammps(lmp, test_config, true);
1211 if (!verbose) ::testing::internal::GetCapturedStdout();
1212
1213 f = lmp->atom->f;
1214 tag = lmp->atom->tag;
1215 stats.reset();
1216 ASSERT_EQ(nlocal + 1, f_ref.size());
1217 for (int i = 0; i < nlocal; ++i) {
1218 EXPECT_FP_LE_WITH_EPS(f[i][0], f_ref[tag[i]].x, 5 * epsilon);
1219 EXPECT_FP_LE_WITH_EPS(f[i][1], f_ref[tag[i]].y, 5 * epsilon);
1220 EXPECT_FP_LE_WITH_EPS(f[i][2], f_ref[tag[i]].z, 5 * epsilon);
1221 }
1222 if (print_stats) std::cerr << "nofdotr_forces stats:" << stats << std::endl;
1223
1224 pair = lmp->force->pair;
1225 stress = pair->virial;
1226 stats.reset();
1227 EXPECT_FP_LE_WITH_EPS(stress[0], test_config.init_stress.xx, 10 * epsilon);
1228 EXPECT_FP_LE_WITH_EPS(stress[1], test_config.init_stress.yy, 10 * epsilon);
1229 EXPECT_FP_LE_WITH_EPS(stress[2], test_config.init_stress.zz, 10 * epsilon);
1230 EXPECT_FP_LE_WITH_EPS(stress[3], test_config.init_stress.xy, 10 * epsilon);
1231 EXPECT_FP_LE_WITH_EPS(stress[4], test_config.init_stress.xz, 10 * epsilon);
1232 EXPECT_FP_LE_WITH_EPS(stress[5], test_config.init_stress.yz, 10 * epsilon);
1233 if (print_stats) std::cerr << "nofdotr_stress stats:" << stats << std::endl;
1234
1235 stats.reset();
1236 EXPECT_FP_LE_WITH_EPS(pair->eng_vdwl, test_config.init_vdwl, 5 * epsilon);
1237 EXPECT_FP_LE_WITH_EPS(pair->eng_coul, test_config.init_coul, 5 * epsilon);
1238 if (print_stats) std::cerr << "nofdotr_energy stats:" << stats << std::endl;
1239
1240 if (!verbose) ::testing::internal::CaptureStdout();
1241 cleanup_lammps(lmp, test_config);
1242 if (!verbose) ::testing::internal::GetCapturedStdout();
1243 };
1244
TEST(PairStyle,single)1245 TEST(PairStyle, single)
1246 {
1247 if (test_config.skip_tests.count(test_info_->name())) GTEST_SKIP();
1248
1249 const char *args[] = {"PairStyle", "-log", "none", "-echo", "screen", "-nocite"};
1250
1251 char **argv = (char **)args;
1252 int argc = sizeof(args) / sizeof(char *);
1253
1254 // need to add this dependency
1255 test_config.prerequisites.push_back(std::make_pair("atom", "full"));
1256
1257 // create a LAMMPS instance with standard settings to detect the number of atom types
1258 if (!verbose) ::testing::internal::CaptureStdout();
1259 LAMMPS *lmp = init_lammps(argc, argv, test_config);
1260 if (!verbose) ::testing::internal::GetCapturedStdout();
1261
1262 if (!lmp) {
1263 std::cerr << "One or more prerequisite styles are not available "
1264 "in this LAMMPS configuration:\n";
1265 for (auto &prerequisite : test_config.prerequisites) {
1266 std::cerr << prerequisite.first << "_style " << prerequisite.second << "\n";
1267 }
1268 test_config.prerequisites.pop_back();
1269 if (!verbose) ::testing::internal::CaptureStdout();
1270 cleanup_lammps(lmp, test_config);
1271 if (!verbose) ::testing::internal::GetCapturedStdout();
1272 GTEST_SKIP();
1273 }
1274 test_config.prerequisites.pop_back();
1275
1276 // gather some information and skip if unsupported
1277 int ntypes = lmp->atom->ntypes;
1278 int molecular = lmp->atom->molecular;
1279 if (molecular > Atom::MOLECULAR) {
1280 std::cerr << "Only atomic and simple molecular atom styles are supported\n";
1281 if (!verbose) ::testing::internal::CaptureStdout();
1282 cleanup_lammps(lmp, test_config);
1283 if (!verbose) ::testing::internal::GetCapturedStdout();
1284 GTEST_SKIP();
1285 }
1286
1287 Pair *pair = lmp->force->pair;
1288 if (!pair->single_enable) {
1289 std::cerr << "Single method not available for pair style " << test_config.pair_style
1290 << std::endl;
1291 if (!verbose) ::testing::internal::CaptureStdout();
1292 cleanup_lammps(lmp, test_config);
1293 if (!verbose) ::testing::internal::GetCapturedStdout();
1294 GTEST_SKIP();
1295 }
1296
1297 if (!pair->compute_flag) {
1298 std::cerr << "Pair style disabled" << std::endl;
1299 if (!verbose) ::testing::internal::CaptureStdout();
1300 cleanup_lammps(lmp, test_config);
1301 if (!verbose) ::testing::internal::GetCapturedStdout();
1302 GTEST_SKIP();
1303 }
1304
1305 // now start over
1306 if (!verbose) ::testing::internal::CaptureStdout();
1307
1308 // utility lambda to improve readability
1309 auto command = [&](const std::string &line) {
1310 lmp->input->one(line.c_str());
1311 };
1312
1313 command("clear");
1314 command("variable newton_pair delete");
1315 command("variable newton_pair index on");
1316
1317 command("variable input_dir index " + INPUT_FOLDER);
1318
1319 for (auto &pre_command : test_config.pre_commands) {
1320 command(pre_command);
1321 }
1322
1323 command("atom_style full");
1324 command("units ${units}");
1325 command("boundary p p p");
1326 command("newton ${newton_pair} ${newton_bond}");
1327
1328 if (molecular == Atom::MOLECULAR) {
1329 command("special_bonds lj/coul "
1330 "${bond_factor} ${angle_factor} ${dihedral_factor}");
1331 }
1332
1333 command("atom_modify map array");
1334 command("region box block -10.0 10.0 -10.0 10.0 -10.0 10.0 units box");
1335
1336 auto cmd = fmt::format("create_box {} box", ntypes);
1337 if (molecular == Atom::MOLECULAR) {
1338 cmd += " bond/types 1"
1339 " extra/bond/per/atom 1"
1340 " extra/special/per/atom 1";
1341 }
1342 command(cmd);
1343
1344 command("pair_style " + test_config.pair_style);
1345
1346 pair = lmp->force->pair;
1347
1348 for (auto &pair_coeff : test_config.pair_coeff) {
1349 command("pair_coeff " + pair_coeff);
1350 }
1351
1352 // create (only) two atoms
1353
1354 command("mass * 1.0");
1355 command("create_atoms 1 single 0.0 -0.75 0.4 units box");
1356 command("create_atoms 2 single 1.5 0.25 -0.1 units box");
1357 command("set atom 1 charge -0.5");
1358 command("set atom 2 charge 0.5");
1359 command("set atom 1 mol 1");
1360 command("set atom 2 mol 2");
1361 command("special_bonds lj/coul 1.0 1.0 1.0");
1362
1363 if (molecular == Atom::MOLECULAR) {
1364 command("create_bonds single/bond 1 1 2");
1365 command("bond_style zero");
1366 command("bond_coeff 1 2.0");
1367 }
1368
1369 for (auto &post_command : test_config.post_commands) {
1370 command(post_command);
1371 }
1372
1373 command("run 0 post no");
1374 if (!verbose) ::testing::internal::GetCapturedStdout();
1375
1376 int idx1 = lmp->atom->map(1);
1377 int idx2 = lmp->atom->map(2);
1378 double epsilon = test_config.epsilon;
1379 double **f = lmp->atom->f;
1380 double **x = lmp->atom->x;
1381 double delx = x[idx2][0] - x[idx1][0];
1382 double dely = x[idx2][1] - x[idx1][1];
1383 double delz = x[idx2][2] - x[idx1][2];
1384 double rsq = delx * delx + dely * dely + delz * delz;
1385 double fsingle = 0.0;
1386 double epair[4], esngl[4];
1387 double splj = lmp->force->special_lj[1];
1388 double spcl = lmp->force->special_coul[1];
1389 ErrorStats stats;
1390
1391 epair[0] = pair->eng_vdwl + pair->eng_coul;
1392 esngl[0] = pair->single(idx1, idx2, 1, 2, rsq, spcl, splj, fsingle);
1393 EXPECT_FP_LE_WITH_EPS(f[idx1][0], -fsingle * delx, epsilon);
1394 EXPECT_FP_LE_WITH_EPS(f[idx1][1], -fsingle * dely, epsilon);
1395 EXPECT_FP_LE_WITH_EPS(f[idx1][2], -fsingle * delz, epsilon);
1396 EXPECT_FP_LE_WITH_EPS(f[idx2][0], fsingle * delx, epsilon);
1397 EXPECT_FP_LE_WITH_EPS(f[idx2][1], fsingle * dely, epsilon);
1398 EXPECT_FP_LE_WITH_EPS(f[idx2][2], fsingle * delz, epsilon);
1399
1400 if (!verbose) ::testing::internal::CaptureStdout();
1401 command("displace_atoms all random 0.5 0.5 0.5 723456");
1402 command("run 0 post no");
1403 if (!verbose) ::testing::internal::GetCapturedStdout();
1404
1405 f = lmp->atom->f;
1406 x = lmp->atom->x;
1407 idx1 = lmp->atom->map(1);
1408 idx2 = lmp->atom->map(2);
1409 delx = x[idx2][0] - x[idx1][0];
1410 dely = x[idx2][1] - x[idx1][1];
1411 delz = x[idx2][2] - x[idx1][2];
1412 rsq = delx * delx + dely * dely + delz * delz;
1413 fsingle = 0.0;
1414
1415 epair[1] = pair->eng_vdwl + pair->eng_coul;
1416 esngl[1] = pair->single(idx1, idx2, 1, 2, rsq, spcl, splj, fsingle);
1417 EXPECT_FP_LE_WITH_EPS(f[idx1][0], -fsingle * delx, epsilon);
1418 EXPECT_FP_LE_WITH_EPS(f[idx1][1], -fsingle * dely, epsilon);
1419 EXPECT_FP_LE_WITH_EPS(f[idx1][2], -fsingle * delz, epsilon);
1420 EXPECT_FP_LE_WITH_EPS(f[idx2][0], fsingle * delx, epsilon);
1421 EXPECT_FP_LE_WITH_EPS(f[idx2][1], fsingle * dely, epsilon);
1422 EXPECT_FP_LE_WITH_EPS(f[idx2][2], fsingle * delz, epsilon);
1423
1424 if (!verbose) ::testing::internal::CaptureStdout();
1425 command("displace_atoms all random 0.5 0.5 0.5 3456963");
1426 command("run 0 post no");
1427 if (!verbose) ::testing::internal::GetCapturedStdout();
1428
1429 f = lmp->atom->f;
1430 x = lmp->atom->x;
1431 idx1 = lmp->atom->map(1);
1432 idx2 = lmp->atom->map(2);
1433 delx = x[idx2][0] - x[idx1][0];
1434 dely = x[idx2][1] - x[idx1][1];
1435 delz = x[idx2][2] - x[idx1][2];
1436 rsq = delx * delx + dely * dely + delz * delz;
1437 fsingle = 0.0;
1438
1439 epair[2] = pair->eng_vdwl + pair->eng_coul;
1440 esngl[2] = pair->single(idx1, idx2, 1, 2, rsq, spcl, splj, fsingle);
1441 EXPECT_FP_LE_WITH_EPS(f[idx1][0], -fsingle * delx, epsilon);
1442 EXPECT_FP_LE_WITH_EPS(f[idx1][1], -fsingle * dely, epsilon);
1443 EXPECT_FP_LE_WITH_EPS(f[idx1][2], -fsingle * delz, epsilon);
1444 EXPECT_FP_LE_WITH_EPS(f[idx2][0], fsingle * delx, epsilon);
1445 EXPECT_FP_LE_WITH_EPS(f[idx2][1], fsingle * dely, epsilon);
1446 EXPECT_FP_LE_WITH_EPS(f[idx2][2], fsingle * delz, epsilon);
1447
1448 if (!verbose) ::testing::internal::CaptureStdout();
1449 command("displace_atoms all random 0.5 0.5 0.5 9726532");
1450 command("run 0 post no");
1451 if (!verbose) ::testing::internal::GetCapturedStdout();
1452
1453 f = lmp->atom->f;
1454 x = lmp->atom->x;
1455 idx1 = lmp->atom->map(1);
1456 idx2 = lmp->atom->map(2);
1457 delx = x[idx2][0] - x[idx1][0];
1458 dely = x[idx2][1] - x[idx1][1];
1459 delz = x[idx2][2] - x[idx1][2];
1460 rsq = delx * delx + dely * dely + delz * delz;
1461 fsingle = 0.0;
1462
1463 epair[3] = pair->eng_vdwl + pair->eng_coul;
1464 esngl[3] = pair->single(idx1, idx2, 1, 2, rsq, spcl, splj, fsingle);
1465 EXPECT_FP_LE_WITH_EPS(f[idx1][0], -fsingle * delx, epsilon);
1466 EXPECT_FP_LE_WITH_EPS(f[idx1][1], -fsingle * dely, epsilon);
1467 EXPECT_FP_LE_WITH_EPS(f[idx1][2], -fsingle * delz, epsilon);
1468 EXPECT_FP_LE_WITH_EPS(f[idx2][0], fsingle * delx, epsilon);
1469 EXPECT_FP_LE_WITH_EPS(f[idx2][1], fsingle * dely, epsilon);
1470 EXPECT_FP_LE_WITH_EPS(f[idx2][2], fsingle * delz, epsilon);
1471 if (print_stats) std::cerr << "single_force stats:" << stats << std::endl;
1472
1473 if ((test_config.pair_style.find("coul/dsf") != std::string::npos) &&
1474 (test_config.pair_style.find("coul/wolf") != std::string::npos)) {
1475 stats.reset();
1476 EXPECT_FP_LE_WITH_EPS(epair[0], esngl[0], epsilon);
1477 EXPECT_FP_LE_WITH_EPS(epair[1], esngl[1], epsilon);
1478 EXPECT_FP_LE_WITH_EPS(epair[2], esngl[2], epsilon);
1479 EXPECT_FP_LE_WITH_EPS(epair[3], esngl[3], epsilon);
1480 if (print_stats) std::cerr << "single_energy stats:" << stats << std::endl;
1481 } else if (print_stats)
1482 std::cerr << "skipping single_energy test due to self energy\n";
1483
1484 if (!verbose) ::testing::internal::CaptureStdout();
1485 cleanup_lammps(lmp, test_config);
1486 if (!verbose) ::testing::internal::GetCapturedStdout();
1487 }
1488
TEST(PairStyle,extract)1489 TEST(PairStyle, extract)
1490 {
1491 if (test_config.skip_tests.count(test_info_->name())) GTEST_SKIP();
1492
1493 const char *args[] = {"PairStyle", "-log", "none", "-echo", "screen", "-nocite"};
1494
1495 char **argv = (char **)args;
1496 int argc = sizeof(args) / sizeof(char *);
1497
1498 if (!verbose) ::testing::internal::CaptureStdout();
1499 LAMMPS *lmp = init_lammps(argc, argv, test_config, true);
1500 if (!verbose) ::testing::internal::GetCapturedStdout();
1501
1502 if (!lmp) {
1503 std::cerr << "One or more prerequisite styles are not available "
1504 "in this LAMMPS configuration:\n";
1505 for (const auto &prerequisite : test_config.prerequisites) {
1506 std::cerr << prerequisite.first << "_style " << prerequisite.second << "\n";
1507 }
1508 GTEST_SKIP();
1509 }
1510
1511 auto pair = lmp->force->pair;
1512 if (!pair->compute_flag) {
1513 std::cerr << "Pair style disabled" << std::endl;
1514 if (!verbose) ::testing::internal::CaptureStdout();
1515 cleanup_lammps(lmp, test_config);
1516 if (!verbose) ::testing::internal::GetCapturedStdout();
1517 GTEST_SKIP();
1518 }
1519
1520 void *ptr = nullptr;
1521 int dim = 0;
1522 for (auto &extract : test_config.extract) {
1523 ptr = pair->extract(extract.first.c_str(), dim);
1524 EXPECT_NE(ptr, nullptr);
1525 EXPECT_EQ(dim, extract.second);
1526 }
1527 ptr = pair->extract("does_not_exist", dim);
1528 EXPECT_EQ(ptr, nullptr);
1529
1530 // replace pair style with the same.
1531 // should just update setting, but not create new style.
1532
1533 int ntypes = lmp->atom->ntypes;
1534 for (int i = 1; i <= ntypes; ++i) {
1535 for (int j = 1; j <= ntypes; ++j) {
1536 pair->cutsq[i][j] = -1.0;
1537 }
1538 }
1539
1540 // utility lambda to improve readability
1541 auto command = [&](const std::string &line) {
1542 lmp->input->one(line.c_str());
1543 };
1544
1545 if (!verbose) ::testing::internal::CaptureStdout();
1546 command("pair_style " + test_config.pair_style);
1547 EXPECT_EQ(pair, lmp->force->pair);
1548
1549 for (auto &pair_coeff : test_config.pair_coeff) {
1550 command("pair_coeff " + pair_coeff);
1551 }
1552 pair->init();
1553 if (!verbose) ::testing::internal::GetCapturedStdout();
1554
1555 for (int i = 1; i <= ntypes; ++i) {
1556 for (int j = 1; j <= ntypes; ++j) {
1557 EXPECT_GE(pair->cutsq[i][j], 0.0);
1558 }
1559 }
1560
1561 if (!verbose) ::testing::internal::CaptureStdout();
1562 cleanup_lammps(lmp, test_config);
1563 if (!verbose) ::testing::internal::GetCapturedStdout();
1564 }
1565