1 // clang-format off
2 /* ----------------------------------------------------------------------
3 LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
4 https://www.lammps.org/, Sandia National Laboratories
5 Steve Plimpton, sjplimp@sandia.gov
6
7 Copyright (2003) Sandia Corporation. Under the terms of Contract
8 DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
9 certain rights in this software. This software is distributed under
10 the GNU General Public License.
11
12 See the README file in the top-level LAMMPS directory.
13 ------------------------------------------------------------------------- */
14
15 #include "force.h"
16
17 #include "style_angle.h" // IWYU pragma: keep
18 #include "style_bond.h" // IWYU pragma: keep
19 #include "style_dihedral.h" // IWYU pragma: keep
20 #include "style_improper.h" // IWYU pragma: keep
21 #include "style_kspace.h" // IWYU pragma: keep
22 #include "style_pair.h" // IWYU pragma: keep
23
24 #include "angle_hybrid.h"
25 #include "bond_hybrid.h"
26 #include "dihedral_hybrid.h"
27 #include "improper_hybrid.h"
28 #include "pair_hybrid.h"
29 #include "kspace.h"
30
31 #include "atom.h"
32 #include "comm.h"
33 #include "error.h"
34
35 #include <cstring>
36
37 using namespace LAMMPS_NS;
38
39 /* ---------------------------------------------------------------------- */
40
Force(LAMMPS * lmp)41 Force::Force(LAMMPS *lmp) : Pointers(lmp)
42 {
43 newton = newton_pair = newton_bond = 1;
44
45 special_lj[0] = special_coul[0] = 1.0;
46 special_lj[1] = special_lj[2] = special_lj[3] = 0.0;
47 special_coul[1] = special_coul[2] = special_coul[3] = 0.0;
48 special_angle = special_dihedral = 0;
49 special_extra = 0;
50
51 dielectric = 1.0;
52 qqr2e_lammps_real = 332.06371; // these constants are toggled
53 qqr2e_charmm_real = 332.0716; // by new CHARMM pair styles
54
55 pair = nullptr;
56 bond = nullptr;
57 angle = nullptr;
58 dihedral = nullptr;
59 improper = nullptr;
60 kspace = nullptr;
61
62 pair_style = utils::strdup("none");
63 bond_style = utils::strdup("none");
64 angle_style = utils::strdup("none");
65 dihedral_style = utils::strdup("none");
66 improper_style = utils::strdup("none");
67 kspace_style = utils::strdup("none");
68
69 pair_restart = nullptr;
70 create_factories();
71 }
72
create_factories()73 void _noopt Force::create_factories()
74 {
75 // fill pair map with pair styles listed in style_pair.h
76
77 pair_map = new PairCreatorMap();
78
79 #define PAIR_CLASS
80 #define PairStyle(key,Class) \
81 (*pair_map)[#key] = &pair_creator<Class>;
82 #include "style_pair.h" // IWYU pragma: keep
83 #undef PairStyle
84 #undef PAIR_CLASS
85
86 bond_map = new BondCreatorMap();
87
88 #define BOND_CLASS
89 #define BondStyle(key,Class) \
90 (*bond_map)[#key] = &bond_creator<Class>;
91 #include "style_bond.h" // IWYU pragma: keep
92 #undef BondStyle
93 #undef BOND_CLASS
94
95 angle_map = new AngleCreatorMap();
96
97 #define ANGLE_CLASS
98 #define AngleStyle(key,Class) \
99 (*angle_map)[#key] = &angle_creator<Class>;
100 #include "style_angle.h" // IWYU pragma: keep
101 #undef AngleStyle
102 #undef ANGLE_CLASS
103
104 dihedral_map = new DihedralCreatorMap();
105
106 #define DIHEDRAL_CLASS
107 #define DihedralStyle(key,Class) \
108 (*dihedral_map)[#key] = &dihedral_creator<Class>;
109 #include "style_dihedral.h" // IWYU pragma: keep
110 #undef DihedralStyle
111 #undef DIHEDRAL_CLASS
112
113 improper_map = new ImproperCreatorMap();
114
115 #define IMPROPER_CLASS
116 #define ImproperStyle(key,Class) \
117 (*improper_map)[#key] = &improper_creator<Class>;
118 #include "style_improper.h" // IWYU pragma: keep
119 #undef ImproperStyle
120 #undef IMPROPER_CLASS
121
122 kspace_map = new KSpaceCreatorMap();
123
124 #define KSPACE_CLASS
125 #define KSpaceStyle(key,Class) \
126 (*kspace_map)[#key] = &kspace_creator<Class>;
127 #include "style_kspace.h" // IWYU pragma: keep
128 #undef KSpaceStyle
129 #undef KSPACE_CLASS
130 }
131
132 /* ---------------------------------------------------------------------- */
133
~Force()134 Force::~Force()
135 {
136 delete [] pair_style;
137 delete [] bond_style;
138 delete [] angle_style;
139 delete [] dihedral_style;
140 delete [] improper_style;
141 delete [] kspace_style;
142
143 delete [] pair_restart;
144
145 if (pair) delete pair;
146 if (bond) delete bond;
147 if (angle) delete angle;
148 if (dihedral) delete dihedral;
149 if (improper) delete improper;
150 if (kspace) delete kspace;
151
152 pair = nullptr;
153 bond = nullptr;
154 angle = nullptr;
155 dihedral = nullptr;
156 improper = nullptr;
157 kspace = nullptr;
158
159 delete pair_map;
160 delete bond_map;
161 delete angle_map;
162 delete dihedral_map;
163 delete improper_map;
164 delete kspace_map;
165 }
166
167 /* ---------------------------------------------------------------------- */
168
init()169 void Force::init()
170 {
171 qqrd2e = qqr2e/dielectric;
172
173 // check if pair style must be specified after restart
174 if (pair_restart) {
175 if (!pair)
176 error->all(FLERR,"Must re-specify non-restarted pair style "
177 "({}) after read_restart", pair_restart);
178 }
179
180 if (kspace) kspace->init(); // kspace must come before pair
181 if (pair) pair->init(); // so g_ewald is defined
182 if (bond) bond->init();
183 if (angle) angle->init();
184 if (dihedral) dihedral->init();
185 if (improper) improper->init();
186
187 // print warnings if topology and force field are inconsistent
188
189 if (comm->me == 0) {
190 if (!bond && (atom->nbonds > 0)) {
191 error->warning(FLERR,"Bonds are defined but no bond style is set");
192 if ((special_lj[1] != 1.0) || (special_coul[1] != 1.0))
193 error->warning(FLERR,"Likewise 1-2 special neighbor interactions != 1.0");
194 }
195 if (!angle && (atom->nangles > 0)) {
196 error->warning(FLERR,"Angles are defined but no angle style is set");
197 if ((special_lj[2] != 1.0) || (special_coul[2] != 1.0))
198 error->warning(FLERR,"Likewise 1-3 special neighbor interactions != 1.0");
199 }
200 if (!dihedral && (atom->ndihedrals > 0)) {
201 error->warning(FLERR,"Dihedrals are defined but no dihedral style is set");
202 if ((special_lj[3] != 1.0) || (special_coul[3] != 1.0))
203 error->warning(FLERR,"Likewise 1-4 special neighbor interactions != 1.0");
204 }
205 if (!improper && (atom->nimpropers > 0))
206 error->warning(FLERR,"Impropers are defined but no improper style is set");
207 }
208 }
209
210 /* ---------------------------------------------------------------------- */
211
setup()212 void Force::setup()
213 {
214 if (pair) pair->setup();
215 }
216
217 /* ----------------------------------------------------------------------
218 create a pair style, called from input script or restart file
219 ------------------------------------------------------------------------- */
220
create_pair(const std::string & style,int trysuffix)221 void Force::create_pair(const std::string &style, int trysuffix)
222 {
223 delete [] pair_style;
224 if (pair) delete pair;
225 if (pair_restart) delete [] pair_restart;
226 pair_style = nullptr;
227 pair = nullptr;
228 pair_restart = nullptr;
229
230 int sflag;
231 pair = new_pair(style,trysuffix,sflag);
232 store_style(pair_style,style,sflag);
233 }
234
235 /* ----------------------------------------------------------------------
236 generate a pair class
237 if trysuffix = 1, try first with suffix1/2 appended
238 return sflag = 0 for no suffix added, 1 or 2 or 3 for suffix1/2/p added
239 special case: if suffixp exists only try suffixp, not suffix
240 ------------------------------------------------------------------------- */
241
new_pair(const std::string & style,int trysuffix,int & sflag)242 Pair *Force::new_pair(const std::string &style, int trysuffix, int &sflag)
243 {
244 if (trysuffix && lmp->suffix_enable) {
245 if (lmp->suffixp) {
246 sflag = 3;
247 std::string estyle = style + "/" + lmp->suffixp;
248 if (pair_map->find(estyle) != pair_map->end()) {
249 PairCreator &pair_creator = (*pair_map)[estyle];
250 return pair_creator(lmp);
251 }
252 }
253 if (lmp->suffix && !lmp->suffixp) {
254 sflag = 1;
255 std::string estyle = style + "/" + lmp->suffix;
256 if (pair_map->find(estyle) != pair_map->end()) {
257 PairCreator &pair_creator = (*pair_map)[estyle];
258 return pair_creator(lmp);
259 }
260 }
261 if (lmp->suffix2) {
262 sflag = 2;
263 std::string estyle = style + "/" + lmp->suffix2;
264 if (pair_map->find(estyle) != pair_map->end()) {
265 PairCreator &pair_creator = (*pair_map)[estyle];
266 return pair_creator(lmp);
267 }
268 }
269 }
270
271 sflag = 0;
272 if (style == "none") return nullptr;
273 if (pair_map->find(style) != pair_map->end()) {
274 PairCreator &pair_creator = (*pair_map)[style];
275 return pair_creator(lmp);
276 }
277
278 error->all(FLERR,utils::check_packages_for_style("pair",style,lmp));
279
280 return nullptr;
281 }
282
283 /* ----------------------------------------------------------------------
284 one instance per pair style in style_pair.h
285 ------------------------------------------------------------------------- */
286
287 template <typename T>
pair_creator(LAMMPS * lmp)288 Pair *Force::pair_creator(LAMMPS *lmp)
289 {
290 return new T(lmp);
291 }
292
293 /* ----------------------------------------------------------------------
294 return ptr to Pair class if matches word or matches hybrid sub-style
295 if exact, then style name must be exact match to word
296 if not exact, style name must contain word
297 if nsub > 0, match Nth hybrid sub-style
298 return nullptr if no match or if nsub=0 and multiple sub-styles match
299 ------------------------------------------------------------------------- */
300
pair_match(const std::string & word,int exact,int nsub)301 Pair *Force::pair_match(const std::string &word, int exact, int nsub)
302 {
303 int iwhich,count;
304
305 if (exact && (word == pair_style)) return pair;
306 else if (!exact && utils::strmatch(pair_style,word)) return pair;
307 else if (utils::strmatch(pair_style,"^hybrid")) {
308 PairHybrid *hybrid = (PairHybrid *) pair;
309 count = 0;
310 for (int i = 0; i < hybrid->nstyles; i++)
311 if ((exact && (word == hybrid->keywords[i])) ||
312 (!exact && utils::strmatch(hybrid->keywords[i],word))) {
313 iwhich = i;
314 count++;
315 if (nsub == count) return hybrid->styles[iwhich];
316 }
317 if (count == 1) return hybrid->styles[iwhich];
318 }
319
320 return nullptr;
321 }
322
323 /* ----------------------------------------------------------------------
324 return style name of Pair class that matches Pair ptr
325 called by Neighbor::print_neigh_info()
326 return nullptr if no match
327 ------------------------------------------------------------------------- */
328
pair_match_ptr(Pair * ptr)329 char *Force::pair_match_ptr(Pair *ptr)
330 {
331 if (ptr == pair) return pair_style;
332
333 if (utils::strmatch(pair_style,"^hybrid")) {
334 PairHybrid *hybrid = (PairHybrid *) pair;
335 for (int i = 0; i < hybrid->nstyles; i++)
336 if (ptr == hybrid->styles[i]) return hybrid->keywords[i];
337 }
338
339 return nullptr;
340 }
341
342 /* ----------------------------------------------------------------------
343 create a bond style, called from input script or restart file
344 ------------------------------------------------------------------------- */
345
create_bond(const std::string & style,int trysuffix)346 void Force::create_bond(const std::string &style, int trysuffix)
347 {
348 delete [] bond_style;
349 if (bond) delete bond;
350
351 int sflag;
352 bond = new_bond(style,trysuffix,sflag);
353 store_style(bond_style,style,sflag);
354 }
355
356 /* ----------------------------------------------------------------------
357 generate a bond class, fist with suffix appended
358 ------------------------------------------------------------------------- */
359
new_bond(const std::string & style,int trysuffix,int & sflag)360 Bond *Force::new_bond(const std::string &style, int trysuffix, int &sflag)
361 {
362 if (trysuffix && lmp->suffix_enable) {
363 if (lmp->suffix) {
364 sflag = 1;
365 std::string estyle = style + "/" + lmp->suffix;
366 if (bond_map->find(estyle) != bond_map->end()) {
367 BondCreator &bond_creator = (*bond_map)[estyle];
368 return bond_creator(lmp);
369 }
370 }
371
372 if (lmp->suffix2) {
373 sflag = 2;
374 std::string estyle = style + "/" + lmp->suffix2;
375 if (bond_map->find(estyle) != bond_map->end()) {
376 BondCreator &bond_creator = (*bond_map)[estyle];
377 return bond_creator(lmp);
378 }
379 }
380 }
381
382 sflag = 0;
383 if (style == "none") return nullptr;
384 if (bond_map->find(style) != bond_map->end()) {
385 BondCreator &bond_creator = (*bond_map)[style];
386 return bond_creator(lmp);
387 }
388
389 error->all(FLERR,utils::check_packages_for_style("bond",style,lmp));
390
391 return nullptr;
392 }
393
394 /* ----------------------------------------------------------------------
395 one instance per bond style in style_bond.h
396 ------------------------------------------------------------------------- */
397
398 template <typename T>
bond_creator(LAMMPS * lmp)399 Bond *Force::bond_creator(LAMMPS *lmp)
400 {
401 return new T(lmp);
402 }
403
404 /* ----------------------------------------------------------------------
405 return ptr to current bond class or hybrid sub-class if matches style
406 ------------------------------------------------------------------------- */
407
bond_match(const std::string & style)408 Bond *Force::bond_match(const std::string &style)
409 {
410 if (style == bond_style) return bond;
411 else if (strcmp(bond_style,"hybrid") == 0) {
412 BondHybrid *hybrid = (BondHybrid *) bond;
413 for (int i = 0; i < hybrid->nstyles; i++)
414 if (style == hybrid->keywords[i]) return hybrid->styles[i];
415 }
416 return nullptr;
417 }
418
419 /* ----------------------------------------------------------------------
420 create an angle style, called from input script or restart file
421 ------------------------------------------------------------------------- */
422
create_angle(const std::string & style,int trysuffix)423 void Force::create_angle(const std::string &style, int trysuffix)
424 {
425 delete [] angle_style;
426 if (angle) delete angle;
427
428 int sflag;
429 angle = new_angle(style,trysuffix,sflag);
430 store_style(angle_style,style,sflag);
431 }
432
433 /* ----------------------------------------------------------------------
434 generate an angle class
435 ------------------------------------------------------------------------- */
436
new_angle(const std::string & style,int trysuffix,int & sflag)437 Angle *Force::new_angle(const std::string &style, int trysuffix, int &sflag)
438 {
439 if (trysuffix && lmp->suffix_enable) {
440 if (lmp->suffix) {
441 sflag = 1;
442 std::string estyle = style + "/" + lmp->suffix;
443 if (angle_map->find(estyle) != angle_map->end()) {
444 AngleCreator &angle_creator = (*angle_map)[estyle];
445 return angle_creator(lmp);
446 }
447 }
448
449 if (lmp->suffix2) {
450 sflag = 2;
451 std::string estyle = style + "/" + lmp->suffix2;
452 if (angle_map->find(estyle) != angle_map->end()) {
453 AngleCreator &angle_creator = (*angle_map)[estyle];
454 return angle_creator(lmp);
455 }
456 }
457 }
458
459 sflag = 0;
460 if (style == "none") return nullptr;
461 if (angle_map->find(style) != angle_map->end()) {
462 AngleCreator &angle_creator = (*angle_map)[style];
463 return angle_creator(lmp);
464 }
465
466 error->all(FLERR,utils::check_packages_for_style("angle",style,lmp));
467
468 return nullptr;
469 }
470
471 /* ----------------------------------------------------------------------
472 one instance per angle style in style_angle.h
473 ------------------------------------------------------------------------- */
474
475 template <typename T>
angle_creator(LAMMPS * lmp)476 Angle *Force::angle_creator(LAMMPS *lmp)
477 {
478 return new T(lmp);
479 }
480
481 /* ----------------------------------------------------------------------
482 return ptr to current angle class or hybrid sub-class if matches style
483 ------------------------------------------------------------------------- */
484
angle_match(const std::string & style)485 Angle *Force::angle_match(const std::string &style)
486 {
487 if (style == angle_style) return angle;
488 else if (utils::strmatch(angle_style,"^hybrid")) {
489 AngleHybrid *hybrid = (AngleHybrid *) angle;
490 for (int i = 0; i < hybrid->nstyles; i++)
491 if (style == hybrid->keywords[i]) return hybrid->styles[i];
492 }
493 return nullptr;
494 }
495
496 /* ----------------------------------------------------------------------
497 create a dihedral style, called from input script or restart file
498 ------------------------------------------------------------------------- */
499
create_dihedral(const std::string & style,int trysuffix)500 void Force::create_dihedral(const std::string &style, int trysuffix)
501 {
502 delete [] dihedral_style;
503 if (dihedral) delete dihedral;
504
505 int sflag;
506 dihedral = new_dihedral(style,trysuffix,sflag);
507 store_style(dihedral_style,style,sflag);
508 }
509
510 /* ----------------------------------------------------------------------
511 generate a dihedral class
512 ------------------------------------------------------------------------- */
513
new_dihedral(const std::string & style,int trysuffix,int & sflag)514 Dihedral *Force::new_dihedral(const std::string &style, int trysuffix, int &sflag)
515 {
516 if (trysuffix && lmp->suffix_enable) {
517 if (lmp->suffix) {
518 sflag = 1;
519 std::string estyle = style + "/" + lmp->suffix;
520 if (dihedral_map->find(estyle) != dihedral_map->end()) {
521 DihedralCreator &dihedral_creator = (*dihedral_map)[estyle];
522 return dihedral_creator(lmp);
523 }
524 }
525
526 if (lmp->suffix2) {
527 sflag = 2;
528 std::string estyle = style + "/" + lmp->suffix2;
529 if (dihedral_map->find(estyle) != dihedral_map->end()) {
530 DihedralCreator &dihedral_creator = (*dihedral_map)[estyle];
531 return dihedral_creator(lmp);
532 }
533 }
534 }
535
536 sflag = 0;
537 if (style == "none") return nullptr;
538 if (dihedral_map->find(style) != dihedral_map->end()) {
539 DihedralCreator &dihedral_creator = (*dihedral_map)[style];
540 return dihedral_creator(lmp);
541 }
542
543 error->all(FLERR,utils::check_packages_for_style("dihedral",style,lmp));
544
545 return nullptr;
546 }
547
548 /* ----------------------------------------------------------------------
549 one instance per dihedral style in style_dihedral.h
550 ------------------------------------------------------------------------- */
551
552 template <typename T>
dihedral_creator(LAMMPS * lmp)553 Dihedral *Force::dihedral_creator(LAMMPS *lmp)
554 {
555 return new T(lmp);
556 }
557
558 /* ----------------------------------------------------------------------
559 return ptr to current angle class or hybrid sub-class if matches style
560 ------------------------------------------------------------------------- */
561
dihedral_match(const std::string & style)562 Dihedral *Force::dihedral_match(const std::string &style)
563 {
564 if (style == dihedral_style) return dihedral;
565 else if (utils::strmatch(dihedral_style,"^hybrid")) {
566 DihedralHybrid *hybrid = (DihedralHybrid *) dihedral;
567 for (int i = 0; i < hybrid->nstyles; i++)
568 if (style == hybrid->keywords[i]) return hybrid->styles[i];
569 }
570 return nullptr;
571 }
572
573 /* ----------------------------------------------------------------------
574 create an improper style, called from input script or restart file
575 ------------------------------------------------------------------------- */
576
create_improper(const std::string & style,int trysuffix)577 void Force::create_improper(const std::string &style, int trysuffix)
578 {
579 delete [] improper_style;
580 if (improper) delete improper;
581
582 int sflag;
583 improper = new_improper(style,trysuffix,sflag);
584 store_style(improper_style,style,sflag);
585 }
586
587 /* ----------------------------------------------------------------------
588 generate a improper class
589 ------------------------------------------------------------------------- */
590
new_improper(const std::string & style,int trysuffix,int & sflag)591 Improper *Force::new_improper(const std::string &style, int trysuffix, int &sflag)
592 {
593 if (trysuffix && lmp->suffix_enable) {
594 if (lmp->suffix) {
595 sflag = 1;
596 std::string estyle = style + "/" + lmp->suffix;
597 if (improper_map->find(estyle) != improper_map->end()) {
598 ImproperCreator &improper_creator = (*improper_map)[estyle];
599 return improper_creator(lmp);
600 }
601 }
602
603 if (lmp->suffix2) {
604 sflag = 2;
605 std::string estyle = style + "/" + lmp->suffix2;
606 if (improper_map->find(estyle) != improper_map->end()) {
607 ImproperCreator &improper_creator = (*improper_map)[estyle];
608 return improper_creator(lmp);
609 }
610 }
611 }
612
613 sflag = 0;
614 if (style == "none") return nullptr;
615 if (improper_map->find(style) != improper_map->end()) {
616 ImproperCreator &improper_creator = (*improper_map)[style];
617 return improper_creator(lmp);
618 }
619
620 error->all(FLERR,utils::check_packages_for_style("improper",style,lmp));
621
622 return nullptr;
623 }
624
625 /* ----------------------------------------------------------------------
626 one instance per improper style in style_improper.h
627 ------------------------------------------------------------------------- */
628
629 template <typename T>
improper_creator(LAMMPS * lmp)630 Improper *Force::improper_creator(LAMMPS *lmp)
631 {
632 return new T(lmp);
633 }
634
635 /* ----------------------------------------------------------------------
636 return ptr to current improper class or hybrid sub-class if matches style
637 ------------------------------------------------------------------------- */
638
improper_match(const std::string & style)639 Improper *Force::improper_match(const std::string &style)
640 {
641 if (style == improper_style) return improper;
642 else if (utils::strmatch(improper_style,"^hybrid")) {
643 ImproperHybrid *hybrid = (ImproperHybrid *) improper;
644 for (int i = 0; i < hybrid->nstyles; i++)
645 if (style == hybrid->keywords[i]) return hybrid->styles[i];
646 }
647 return nullptr;
648 }
649
650 /* ----------------------------------------------------------------------
651 new kspace style
652 ------------------------------------------------------------------------- */
653
create_kspace(const std::string & style,int trysuffix)654 void Force::create_kspace(const std::string &style, int trysuffix)
655 {
656 delete [] kspace_style;
657 if (kspace) delete kspace;
658
659 int sflag;
660 kspace = new_kspace(style,trysuffix,sflag);
661 store_style(kspace_style,style,sflag);
662 }
663
664 /* ----------------------------------------------------------------------
665 generate a kspace class
666 ------------------------------------------------------------------------- */
667
new_kspace(const std::string & style,int trysuffix,int & sflag)668 KSpace *Force::new_kspace(const std::string &style, int trysuffix, int &sflag)
669 {
670 if (trysuffix && lmp->suffix_enable) {
671 if (lmp->suffix) {
672 sflag = 1;
673 std::string estyle = style + "/" + lmp->suffix;
674 if (kspace_map->find(estyle) != kspace_map->end()) {
675 KSpaceCreator &kspace_creator = (*kspace_map)[estyle];
676 return kspace_creator(lmp);
677 }
678 }
679
680 if (lmp->suffix2) {
681 sflag = 2;
682 std::string estyle = style + "/" + lmp->suffix2;
683 if (kspace_map->find(estyle) != kspace_map->end()) {
684 KSpaceCreator &kspace_creator = (*kspace_map)[estyle];
685 return kspace_creator(lmp);
686 }
687 }
688 }
689
690 sflag = 0;
691 if (style == "none") return nullptr;
692 if (kspace_map->find(style) != kspace_map->end()) {
693 KSpaceCreator &kspace_creator = (*kspace_map)[style];
694 return kspace_creator(lmp);
695 }
696
697 error->all(FLERR,utils::check_packages_for_style("kspace",style,lmp));
698
699 return nullptr;
700 }
701
702 /* ----------------------------------------------------------------------
703 one instance per kspace style in style_kspace.h
704 ------------------------------------------------------------------------- */
705
706 template <typename T>
kspace_creator(LAMMPS * lmp)707 KSpace *Force::kspace_creator(LAMMPS *lmp)
708 {
709 return new T(lmp);
710 }
711
712 /* ----------------------------------------------------------------------
713 return ptr to Kspace class if matches word
714 if exact, then style name must be exact match to word
715 if not exact, style name must contain word
716 return nullptr if no match
717 ------------------------------------------------------------------------- */
718
kspace_match(const std::string & word,int exact)719 KSpace *Force::kspace_match(const std::string &word, int exact)
720 {
721 if (exact && (word == kspace_style)) return kspace;
722 else if (!exact && utils::strmatch(kspace_style,word)) return kspace;
723 return nullptr;
724 }
725
726 /* ----------------------------------------------------------------------
727 store style name in str allocated here
728 if sflag = 0, no suffix
729 if sflag = 1/2/3, append suffix or suffix2 or suffixp to style
730 ------------------------------------------------------------------------- */
731
store_style(char * & str,const std::string & style,int sflag)732 void Force::store_style(char *&str, const std::string &style, int sflag)
733 {
734 std::string estyle = style;
735
736 if (sflag == 1) estyle += std::string("/") + lmp->suffix;
737 else if (sflag == 2) estyle += std::string("/") + lmp->suffix2;
738 else if (sflag == 3) estyle += std::string("/") + lmp->suffixp;
739 str = utils::strdup(estyle);
740 }
741
742 /* ----------------------------------------------------------------------
743 set special bond values
744 ------------------------------------------------------------------------- */
745
set_special(int narg,char ** arg)746 void Force::set_special(int narg, char **arg)
747 {
748 if (narg == 0) error->all(FLERR,"Illegal special_bonds command");
749
750 // defaults, but do not reset special_extra
751
752 special_lj[1] = special_lj[2] = special_lj[3] = 0.0;
753 special_coul[1] = special_coul[2] = special_coul[3] = 0.0;
754 special_angle = special_dihedral = 0;
755
756 int iarg = 0;
757 while (iarg < narg) {
758 if (strcmp(arg[iarg],"amber") == 0) {
759 if (iarg+1 > narg) error->all(FLERR,"Illegal special_bonds command");
760 special_lj[1] = 0.0;
761 special_lj[2] = 0.0;
762 special_lj[3] = 0.5;
763 special_coul[1] = 0.0;
764 special_coul[2] = 0.0;
765 special_coul[3] = 5.0/6.0;
766 iarg += 1;
767 } else if (strcmp(arg[iarg],"charmm") == 0) {
768 if (iarg+1 > narg) error->all(FLERR,"Illegal special_bonds command");
769 special_lj[1] = 0.0;
770 special_lj[2] = 0.0;
771 special_lj[3] = 0.0;
772 special_coul[1] = 0.0;
773 special_coul[2] = 0.0;
774 special_coul[3] = 0.0;
775 iarg += 1;
776 } else if (strcmp(arg[iarg],"dreiding") == 0) {
777 if (iarg+1 > narg) error->all(FLERR,"Illegal special_bonds command");
778 special_lj[1] = 0.0;
779 special_lj[2] = 0.0;
780 special_lj[3] = 1.0;
781 special_coul[1] = 0.0;
782 special_coul[2] = 0.0;
783 special_coul[3] = 1.0;
784 iarg += 1;
785 } else if (strcmp(arg[iarg],"fene") == 0) {
786 if (iarg+1 > narg) error->all(FLERR,"Illegal special_bonds command");
787 special_lj[1] = 0.0;
788 special_lj[2] = 1.0;
789 special_lj[3] = 1.0;
790 special_coul[1] = 0.0;
791 special_coul[2] = 1.0;
792 special_coul[3] = 1.0;
793 iarg += 1;
794 } else if (strcmp(arg[iarg],"lj/coul") == 0) {
795 if (iarg+4 > narg) error->all(FLERR,"Illegal special_bonds command");
796 special_lj[1] = special_coul[1] = utils::numeric(FLERR,arg[iarg+1],false,lmp);
797 special_lj[2] = special_coul[2] = utils::numeric(FLERR,arg[iarg+2],false,lmp);
798 special_lj[3] = special_coul[3] = utils::numeric(FLERR,arg[iarg+3],false,lmp);
799 iarg += 4;
800 } else if (strcmp(arg[iarg],"lj") == 0) {
801 if (iarg+4 > narg) error->all(FLERR,"Illegal special_bonds command");
802 special_lj[1] = utils::numeric(FLERR,arg[iarg+1],false,lmp);
803 special_lj[2] = utils::numeric(FLERR,arg[iarg+2],false,lmp);
804 special_lj[3] = utils::numeric(FLERR,arg[iarg+3],false,lmp);
805 iarg += 4;
806 } else if (strcmp(arg[iarg],"coul") == 0) {
807 if (iarg+4 > narg) error->all(FLERR,"Illegal special_bonds command");
808 special_coul[1] = utils::numeric(FLERR,arg[iarg+1],false,lmp);
809 special_coul[2] = utils::numeric(FLERR,arg[iarg+2],false,lmp);
810 special_coul[3] = utils::numeric(FLERR,arg[iarg+3],false,lmp);
811 iarg += 4;
812 } else if (strcmp(arg[iarg],"angle") == 0) {
813 if (iarg+2 > narg) error->all(FLERR,"Illegal special_bonds command");
814 if (strcmp(arg[iarg+1],"no") == 0) special_angle = 0;
815 else if (strcmp(arg[iarg+1],"yes") == 0) special_angle = 1;
816 else error->all(FLERR,"Illegal special_bonds command");
817 iarg += 2;
818 } else if (strcmp(arg[iarg],"dihedral") == 0) {
819 if (iarg+2 > narg) error->all(FLERR,"Illegal special_bonds command");
820 if (strcmp(arg[iarg+1],"no") == 0) special_dihedral = 0;
821 else if (strcmp(arg[iarg+1],"yes") == 0) special_dihedral = 1;
822 else error->all(FLERR,"Illegal special_bonds command");
823 iarg += 2;
824 } else error->all(FLERR,"Illegal special_bonds command");
825 }
826
827 for (int i = 1; i <= 3; i++)
828 if (special_lj[i] < 0.0 || special_lj[i] > 1.0 ||
829 special_coul[i] < 0.0 || special_coul[i] > 1.0)
830 error->all(FLERR,"Illegal special_bonds command");
831 }
832
833 /* ----------------------------------------------------------------------
834 memory usage of force classes
835 ------------------------------------------------------------------------- */
836
memory_usage()837 double Force::memory_usage()
838 {
839 double bytes = 0;
840 if (pair) bytes += pair->memory_usage();
841 if (bond) bytes += bond->memory_usage();
842 if (angle) bytes += angle->memory_usage();
843 if (dihedral) bytes += dihedral->memory_usage();
844 if (improper) bytes += improper->memory_usage();
845 if (kspace) bytes += kspace->memory_usage();
846 return bytes;
847 }
848