1 /* ----------------------------------------------------------------------
2 This is the
3
4 ██╗ ██╗ ██████╗ ██████╗ ██████╗ ██╗ ██╗████████╗███████╗
5 ██║ ██║██╔════╝ ██╔════╝ ██╔════╝ ██║ ██║╚══██╔══╝██╔════╝
6 ██║ ██║██║ ███╗██║ ███╗██║ ███╗███████║ ██║ ███████╗
7 ██║ ██║██║ ██║██║ ██║██║ ██║██╔══██║ ██║ ╚════██║
8 ███████╗██║╚██████╔╝╚██████╔╝╚██████╔╝██║ ██║ ██║ ███████║
9 ╚══════╝╚═╝ ╚═════╝ ╚═════╝ ╚═════╝ ╚═╝ ╚═╝ ╚═╝ ╚══════╝®
10
11 DEM simulation engine, released by
12 DCS Computing Gmbh, Linz, Austria
13 http://www.dcs-computing.com, office@dcs-computing.com
14
15 LIGGGHTS® is part of CFDEM®project:
16 http://www.liggghts.com | http://www.cfdem.com
17
18 Core developer and main author:
19 Christoph Kloss, christoph.kloss@dcs-computing.com
20
21 LIGGGHTS® is open-source, distributed under the terms of the GNU Public
22 License, version 2 or later. It is distributed in the hope that it will
23 be useful, but WITHOUT ANY WARRANTY; without even the implied warranty
24 of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. You should have
25 received a copy of the GNU General Public License along with LIGGGHTS®.
26 If not, see http://www.gnu.org/licenses . See also top-level README
27 and LICENSE files.
28
29 LIGGGHTS® and CFDEM® are registered trade marks of DCS Computing GmbH,
30 the producer of the LIGGGHTS® software and the CFDEM®coupling software
31 See http://www.cfdem.com/terms-trademark-policy for details.
32
33 -------------------------------------------------------------------------
34 Contributing author and copyright for this file:
35 This file is from LAMMPS, but has been modified. Copyright for
36 modification:
37
38 Copyright 2012- DCS Computing GmbH, Linz
39 Copyright 2009-2012 JKU Linz
40
41 Copyright of original file:
42 LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
43 http://lammps.sandia.gov, Sandia National Laboratories
44 Steve Plimpton, sjplimp@sandia.gov
45
46 Copyright (2003) Sandia Corporation. Under the terms of Contract
47 DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
48 certain rights in this software. This software is distributed under
49 the GNU General Public License.
50 ------------------------------------------------------------------------- */
51
52 #include <stdlib.h>
53 #include <string.h>
54 #include "ctype.h"
55 #include "force.h"
56 #include "style_bond.h"
57 #include "style_angle.h"
58 #include "style_dihedral.h"
59 #include "style_improper.h"
60 #include "style_pair.h"
61 #include "style_kspace.h"
62 #include "atom.h"
63 #include "comm.h"
64 #include "pair.h"
65 #include "pair_gran.h"
66 #include "pair_hybrid.h"
67 #include "pair_hybrid_overlay.h"
68 #include "bond.h"
69 #include "bond_hybrid.h"
70 #include "angle.h"
71 #include "dihedral.h"
72 #include "improper.h"
73 #include "kspace.h"
74 #include "group.h"
75 #include "memory.h"
76 #include "error.h"
77
78 using namespace LAMMPS_NS;
79
80 /* ---------------------------------------------------------------------- */
81
Force(LAMMPS * lmp)82 Force::Force(LAMMPS *lmp) : Pointers(lmp), registry(lmp)
83 {
84 newton = newton_pair = newton_bond = 1;
85
86 special_lj[0] = special_coul[0] = 1.0;
87 special_lj[1] = special_lj[2] = special_lj[3] = 0.0;
88 special_coul[1] = special_coul[2] = special_coul[3] = 0.0;
89 special_angle = special_dihedral = 0;
90 special_extra = 0;
91
92 dielectric = 1.0;
93
94 coarsegraining_ = 1.0;
95 error_coarsegraining_ = true;
96 warn_coarsegraining_ = false;
97
98 pair = NULL;
99 bond = NULL;
100 angle = NULL;
101 dihedral = NULL;
102 improper = NULL;
103 kspace = NULL;
104
105 char *str = (char *) "none";
106 int n = strlen(str) + 1;
107 pair_style = new char[n];
108 strcpy(pair_style,str);
109 bond_style = new char[n];
110 strcpy(bond_style,str);
111 angle_style = new char[n];
112 strcpy(angle_style,str);
113 dihedral_style = new char[n];
114 strcpy(dihedral_style,str);
115 improper_style = new char[n];
116 strcpy(improper_style,str);
117 kspace_style = new char[n];
118 strcpy(kspace_style,str);
119
120 // fill pair map with pair styles listed in style_pair.h
121
122 pair_map = new std::map<std::string,PairCreator>();
123
124 #define PAIR_CLASS
125 #define PairStyle(key,Class) \
126 (*pair_map)[#key] = &pair_creator<Class>;
127 #include "style_pair.h"
128 #undef PairStyle
129 #undef PAIR_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 if (pair) delete pair;
144 if (bond) delete bond;
145 if (angle) delete angle;
146 if (dihedral) delete dihedral;
147 if (improper) delete improper;
148 if (kspace) delete kspace;
149
150 delete pair_map;
151 }
152
153 /* ---------------------------------------------------------------------- */
154
init()155 void Force::init()
156 {
157 registry.init();
158 qqrd2e = qqr2e/dielectric;
159
160 if (kspace) kspace->init(); // kspace must come before pair
161 if (pair) pair->init(); // so g_ewald is defined
162 if (bond) bond->init();
163 if (angle) angle->init();
164 if (dihedral) dihedral->init();
165 if (improper) improper->init();
166 if(cg_active() && warn_cg() && atom->ntypes != int(coarsegrainingTypeBased_.size()))
167 error->warningAll(FLERR,"Coarse graining factor not specified for all atom types. will use maximum CG for unspecified atom types.\n\n");
168 }
169
170 /* ----------------------------------------------------------------------
171 create a pair style, called from input script
172 ------------------------------------------------------------------------- */
173
create_pair(const char * style,const char * suffix)174 void Force::create_pair(const char *style, const char *suffix)
175 {
176 delete [] pair_style;
177 if (pair) delete pair;
178
179 int sflag;
180 pair = new_pair(style,suffix,sflag);
181
182 if (sflag) {
183 char estyle[256];
184 sprintf(estyle,"%s/%s",style,suffix);
185 int n = strlen(estyle) + 1;
186 pair_style = new char[n];
187 strcpy(pair_style,estyle);
188 } else {
189 int n = strlen(style) + 1;
190 pair_style = new char[n];
191 strcpy(pair_style,style);
192 }
193 }
194
195 /* ----------------------------------------------------------------------
196 create a pair style, called from restart file
197 ------------------------------------------------------------------------- */
198
create_pair_from_restart(FILE * fp,const char * style,const char * suffix)199 void Force::create_pair_from_restart(FILE * fp, const char *style, const char *suffix)
200 {
201 delete [] pair_style;
202 if (pair) delete pair;
203
204 int sflag;
205 pair = new_pair_from_restart(fp, style,suffix,sflag);
206
207 if (sflag) {
208 char estyle[256];
209 sprintf(estyle,"%s/%s",style,suffix);
210 int n = strlen(estyle) + 1;
211 pair_style = new char[n];
212 strcpy(pair_style,estyle);
213 } else {
214 int n = strlen(style) + 1;
215 pair_style = new char[n];
216 strcpy(pair_style,style);
217 }
218 }
219
220 /* ----------------------------------------------------------------------
221 generate a pair class
222 try first with suffix appended
223 ------------------------------------------------------------------------- */
224
new_pair(const char * style,const char * suffix,int & sflag)225 Pair *Force::new_pair(const char *style, const char *suffix, int &sflag)
226 {
227 if (suffix && lmp->suffix_enable) {
228 sflag = 1;
229 char estyle[256];
230 sprintf(estyle,"%s/%s",style,suffix);
231
232 if (pair_map->find(estyle) != pair_map->end()) {
233 PairCreator pair_creator = (*pair_map)[estyle];
234 return pair_creator(lmp);
235 }
236
237 }
238
239 sflag = 0;
240
241 if (strcmp(style,"none") == 0) return NULL;
242 if (pair_map->find(style) != pair_map->end()) {
243 PairCreator pair_creator = (*pair_map)[style];
244 return pair_creator(lmp);
245 }
246
247 error->all(FLERR,"Invalid pair style");
248
249 return NULL;
250 }
251
252 /* ----------------------------------------------------------------------
253 generate a pair class from restart file, first with suffix appended
254 ------------------------------------------------------------------------- */
255
new_pair_from_restart(FILE * fp,const char * style,const char * suffix,int & sflag)256 Pair *Force::new_pair_from_restart(FILE * fp, const char *style, const char *suffix, int &sflag)
257 {
258 if (suffix && lmp->suffix_enable) {
259 sflag = 1;
260 char estyle[256];
261 sprintf(estyle,"%s/%s",style,suffix);
262
263 if (pair_map->find(estyle) != pair_map->end()) {
264 PairCreator pair_creator = (*pair_map)[estyle];
265 return pair_creator(lmp);
266 }
267 }
268
269 sflag = 0;
270
271 if (strcmp(style,"none") == 0) return NULL;
272 if (pair_map->find(style) != pair_map->end()) {
273 PairCreator pair_creator = (*pair_map)[style];
274 return pair_creator(lmp);
275 }
276
277 error->all(FLERR,"Invalid pair style");
278 return NULL;
279 }
280
281 /* ----------------------------------------------------------------------
282 one instance per pair style in style_pair.h
283 ------------------------------------------------------------------------- */
284
285 template <typename T>
pair_creator(LAMMPS * lmp)286 Pair *Force::pair_creator(LAMMPS *lmp)
287 {
288 return new T(lmp);
289 }
290
291 /* ----------------------------------------------------------------------
292 return ptr to Pair class if matches word or matches hybrid sub-style
293 if exact, then style name must be exact match to word
294 if not exact, style name must contain word
295 return NULL if no match or multiple sub-styles match
296 ------------------------------------------------------------------------- */
297
pair_match(const char * word,int exact)298 Pair *Force::pair_match(const char *word, int exact)
299 {
300 int iwhich,count;
301
302 if (exact && strcmp(pair_style,word) == 0) return pair;
303 else if (!exact && strstr(pair_style,word)) return pair;
304
305 else if (strstr(pair_style,"hybrid/overlay")) {
306 PairHybridOverlay *hybrid = (PairHybridOverlay *) pair;
307 count = 0;
308 for (int i = 0; i < hybrid->nstyles; i++)
309 if ((exact && strcmp(hybrid->keywords[i],word) == 0) ||
310 (!exact && strstr(hybrid->keywords[i],word))) {
311 iwhich = i;
312 count++;
313 }
314 if (count == 1) return hybrid->styles[iwhich];
315
316 } else if (strstr(pair_style,"hybrid")) {
317 PairHybrid *hybrid = (PairHybrid *) pair;
318 count = 0;
319 for (int i = 0; i < hybrid->nstyles; i++)
320 if ((exact && strcmp(hybrid->keywords[i],word) == 0) ||
321 (!exact && strstr(hybrid->keywords[i],word))) {
322 iwhich = i;
323 count++;
324 }
325 if (count == 1) return hybrid->styles[iwhich];
326 }
327
328 return NULL;
329 }
330
331 /* ----------------------------------------------------------------------
332 create a bond style, called from input script or restart file
333 ------------------------------------------------------------------------- */
334
create_bond(const char * style,const char * suffix)335 void Force::create_bond(const char *style, const char *suffix)
336 {
337 delete [] bond_style;
338 if (bond) delete bond;
339
340 int sflag;
341 bond = new_bond(style,suffix,sflag);
342
343 if (sflag) {
344 char estyle[256];
345 sprintf(estyle,"%s/%s",style,suffix);
346 int n = strlen(estyle) + 1;
347 bond_style = new char[n];
348 strcpy(bond_style,estyle);
349 } else {
350 int n = strlen(style) + 1;
351 bond_style = new char[n];
352 strcpy(bond_style,style);
353 }
354 }
355
356 /* ----------------------------------------------------------------------
357 generate a bond class, fist with suffix appended
358 ------------------------------------------------------------------------- */
359
new_bond(const char * style,const char * suffix,int & sflag)360 Bond *Force::new_bond(const char *style, const char *suffix, int &sflag)
361 {
362 if (suffix && lmp->suffix_enable) {
363 sflag = 1;
364 char estyle[256];
365 sprintf(estyle,"%s/%s",style,suffix);
366
367 if (0) return NULL;
368
369 #define BOND_CLASS
370 #define BondStyle(key,Class) \
371 else if (strcmp(estyle,#key) == 0) return new Class(lmp);
372 #include "style_bond.h"
373 #undef BondStyle
374 #undef BOND_CLASS
375 }
376
377 sflag = 0;
378
379 if (strcmp(style,"none") == 0) return NULL;
380
381 #define BOND_CLASS
382 #define BondStyle(key,Class) \
383 else if (strcmp(style,#key) == 0) return new Class(lmp);
384 #include "style_bond.h"
385 #undef BOND_CLASS
386
387 else error->all(FLERR,"Invalid bond style");
388 return NULL;
389 }
390
391 /* ----------------------------------------------------------------------
392 return ptr to current bond class or hybrid sub-class if matches style
393 ------------------------------------------------------------------------- */
394
bond_match(const char * style)395 Bond *Force::bond_match(const char *style)
396 {
397 if (strcmp(bond_style,style) == 0) return bond;
398 else if (strcmp(bond_style,"hybrid") == 0) {
399 BondHybrid *hybrid = (BondHybrid *) bond;
400 for (int i = 0; i < hybrid->nstyles; i++)
401 if (strcmp(hybrid->keywords[i],style) == 0) return hybrid->styles[i];
402 }
403 return NULL;
404 }
405
406 /* ----------------------------------------------------------------------
407 create an angle style, called from input script or restart file
408 ------------------------------------------------------------------------- */
409
create_angle(const char * style,const char * suffix)410 void Force::create_angle(const char *style, const char *suffix)
411 {
412 delete [] angle_style;
413 if (angle) delete angle;
414
415 int sflag;
416 angle = new_angle(style,suffix,sflag);
417
418 if (sflag) {
419 char estyle[256];
420 sprintf(estyle,"%s/%s",style,suffix);
421 int n = strlen(estyle) + 1;
422 angle_style = new char[n];
423 strcpy(angle_style,estyle);
424 } else {
425 int n = strlen(style) + 1;
426 angle_style = new char[n];
427 strcpy(angle_style,style);
428 }
429 }
430
431 /* ----------------------------------------------------------------------
432 generate an angle class
433 ------------------------------------------------------------------------- */
434
new_angle(const char * style,const char * suffix,int & sflag)435 Angle *Force::new_angle(const char *style, const char *suffix, int &sflag)
436 {
437 if (suffix && lmp->suffix_enable) {
438 sflag = 1;
439 char estyle[256];
440 sprintf(estyle,"%s/%s",style,suffix);
441
442 if (0) return NULL;
443
444 #define ANGLE_CLASS
445 #define AngleStyle(key,Class) \
446 else if (strcmp(estyle,#key) == 0) return new Class(lmp);
447 #include "style_angle.h"
448 #undef AngleStyle
449 #undef ANGLE_CLASS
450
451 }
452
453 sflag = 0;
454
455 if (strcmp(style,"none") == 0) return NULL;
456
457 #define ANGLE_CLASS
458 #define AngleStyle(key,Class) \
459 else if (strcmp(style,#key) == 0) return new Class(lmp);
460 #include "style_angle.h"
461 #undef ANGLE_CLASS
462
463 else error->all(FLERR,"Invalid angle style");
464 return NULL;
465 }
466
467 /* ----------------------------------------------------------------------
468 create a dihedral style, called from input script or restart file
469 ------------------------------------------------------------------------- */
470
create_dihedral(const char * style,const char * suffix)471 void Force::create_dihedral(const char *style, const char *suffix)
472 {
473 delete [] dihedral_style;
474 if (dihedral) delete dihedral;
475
476 int sflag;
477 dihedral = new_dihedral(style,suffix,sflag);
478
479 if (sflag) {
480 char estyle[256];
481 sprintf(estyle,"%s/%s",style,suffix);
482 int n = strlen(estyle) + 1;
483 dihedral_style = new char[n];
484 strcpy(dihedral_style,estyle);
485 } else {
486 int n = strlen(style) + 1;
487 dihedral_style = new char[n];
488 strcpy(dihedral_style,style);
489 }
490 }
491
492 /* ----------------------------------------------------------------------
493 generate a dihedral class
494 ------------------------------------------------------------------------- */
495
new_dihedral(const char * style,const char * suffix,int & sflag)496 Dihedral *Force::new_dihedral(const char *style, const char *suffix, int &sflag)
497 {
498 if (suffix && lmp->suffix_enable) {
499 sflag = 1;
500 char estyle[256];
501 sprintf(estyle,"%s/%s",style,suffix);
502
503 if (0) return NULL;
504
505 #define DIHEDRAL_CLASS
506 #define DihedralStyle(key,Class) \
507 else if (strcmp(estyle,#key) == 0) return new Class(lmp);
508 #include "style_dihedral.h"
509 #undef DihedralStyle
510 #undef DIHEDRAL_CLASS
511
512 }
513
514 sflag = 0;
515
516 if (strcmp(style,"none") == 0) return NULL;
517
518 #define DIHEDRAL_CLASS
519 #define DihedralStyle(key,Class) \
520 else if (strcmp(style,#key) == 0) return new Class(lmp);
521 #include "style_dihedral.h"
522 #undef DihedralStyle
523 #undef DIHEDRAL_CLASS
524
525 else error->all(FLERR,"Invalid dihedral style");
526 return NULL;
527 }
528
529 /* ----------------------------------------------------------------------
530 create an improper style, called from input script or restart file
531 ------------------------------------------------------------------------- */
532
create_improper(const char * style,const char * suffix)533 void Force::create_improper(const char *style, const char *suffix)
534 {
535 delete [] improper_style;
536 if (improper) delete improper;
537
538 int sflag;
539 improper = new_improper(style,suffix,sflag);
540
541 if (sflag) {
542 char estyle[256];
543 sprintf(estyle,"%s/%s",style,suffix);
544 int n = strlen(estyle) + 1;
545 improper_style = new char[n];
546 strcpy(improper_style,estyle);
547 } else {
548 int n = strlen(style) + 1;
549 improper_style = new char[n];
550 strcpy(improper_style,style);
551 }
552 }
553
554 /* ----------------------------------------------------------------------
555 generate a improper class
556 ------------------------------------------------------------------------- */
557
new_improper(const char * style,const char * suffix,int & sflag)558 Improper *Force::new_improper(const char *style, const char *suffix, int &sflag)
559 {
560 if (suffix && lmp->suffix_enable) {
561 sflag = 1;
562 char estyle[256];
563 sprintf(estyle,"%s/%s",style,suffix);
564
565 if (0) return NULL;
566
567 #define IMPROPER_CLASS
568 #define ImproperStyle(key,Class) \
569 else if (strcmp(estyle,#key) == 0) return new Class(lmp);
570 #include "style_improper.h"
571 #undef ImproperStyle
572 #undef IMPROPER_CLASS
573
574 }
575
576 sflag = 0;
577
578 if (strcmp(style,"none") == 0) return NULL;
579
580 #define IMPROPER_CLASS
581 #define ImproperStyle(key,Class) \
582 else if (strcmp(style,#key) == 0) return new Class(lmp);
583 #include "style_improper.h"
584 #undef IMPROPER_CLASS
585
586 else error->all(FLERR,"Invalid improper style");
587 return NULL;
588 }
589
590 /* ----------------------------------------------------------------------
591 new kspace style
592 ------------------------------------------------------------------------- */
593
create_kspace(int narg,char ** arg,const char * suffix)594 void Force::create_kspace(int narg, char **arg, const char *suffix)
595 {
596 delete [] kspace_style;
597 if (kspace) delete kspace;
598
599 int sflag;
600 kspace = new_kspace(narg,arg,suffix,sflag);
601
602 if (sflag) {
603 char estyle[256];
604 sprintf(estyle,"%s/%s",arg[0],suffix);
605 int n = strlen(estyle) + 1;
606 kspace_style = new char[n];
607 strcpy(kspace_style,estyle);
608 } else {
609 int n = strlen(arg[0]) + 1;
610 kspace_style = new char[n];
611 strcpy(kspace_style,arg[0]);
612 }
613 }
614
615 /* ----------------------------------------------------------------------
616 generate a kspace class
617 ------------------------------------------------------------------------- */
618
new_kspace(int narg,char ** arg,const char * suffix,int & sflag)619 KSpace *Force::new_kspace(int narg, char **arg, const char *suffix, int &sflag)
620 {
621 if (suffix && lmp->suffix_enable) {
622 sflag = 1;
623 char estyle[256];
624 sprintf(estyle,"%s/%s",arg[0],suffix);
625
626 if (0) return NULL;
627
628 #define KSPACE_CLASS
629 #define KSpaceStyle(key,Class) \
630 else if (strcmp(estyle,#key) == 0) return new Class(lmp,narg-1,&arg[1]);
631 #include "style_kspace.h"
632 #undef KSpaceStyle
633 #undef KSPACE_CLASS
634
635 }
636
637 sflag = 0;
638
639 if (strcmp(arg[0],"none") == 0) return NULL;
640
641 #define KSPACE_CLASS
642 #define KSpaceStyle(key,Class) \
643 else if (strcmp(arg[0],#key) == 0) return new Class(lmp,narg-1,&arg[1]);
644 #include "style_kspace.h"
645 #undef KSPACE_CLASS
646
647 else error->all(FLERR,"Invalid kspace style");
648 return NULL;
649 }
650
651 /* ----------------------------------------------------------------------
652 return ptr to Kspace class if matches word
653 if exact, then style name must be exact match to word
654 if not exact, style name must contain word
655 return NULL if no match
656 ------------------------------------------------------------------------- */
657
kspace_match(const char * word,int exact)658 KSpace *Force::kspace_match(const char *word, int exact)
659 {
660 if (exact && strcmp(kspace_style,word) == 0) return kspace;
661 else if (!exact && strstr(kspace_style,word)) return kspace;
662 return NULL;
663 }
664
665 /* ----------------------------------------------------------------------
666 set special bond values
667 ------------------------------------------------------------------------- */
668
set_special(int narg,char ** arg)669 void Force::set_special(int narg, char **arg)
670 {
671 if (narg == 0) error->all(FLERR,"Illegal special_bonds command");
672
673 int iarg = 0;
674 while (iarg < narg) {
675 if (strcmp(arg[iarg],"amber") == 0) {
676 if (iarg+1 > narg) error->all(FLERR,"Illegal special_bonds command");
677 special_lj[1] = 0.0;
678 special_lj[2] = 0.0;
679 special_lj[3] = 0.5;
680 special_coul[1] = 0.0;
681 special_coul[2] = 0.0;
682 special_coul[3] = 5.0/6.0;
683 iarg += 1;
684 } else if (strcmp(arg[iarg],"charmm") == 0) {
685 if (iarg+1 > narg) error->all(FLERR,"Illegal special_bonds command");
686 special_lj[1] = 0.0;
687 special_lj[2] = 0.0;
688 special_lj[3] = 0.0;
689 special_coul[1] = 0.0;
690 special_coul[2] = 0.0;
691 special_coul[3] = 0.0;
692 iarg += 1;
693 } else if (strcmp(arg[iarg],"dreiding") == 0) {
694 if (iarg+1 > narg) error->all(FLERR,"Illegal special_bonds command");
695 special_lj[1] = 0.0;
696 special_lj[2] = 0.0;
697 special_lj[3] = 1.0;
698 special_coul[1] = 0.0;
699 special_coul[2] = 0.0;
700 special_coul[3] = 1.0;
701 iarg += 1;
702 } else if (strcmp(arg[iarg],"fene") == 0) {
703 if (iarg+1 > narg) error->all(FLERR,"Illegal special_bonds command");
704 special_lj[1] = 0.0;
705 special_lj[2] = 1.0;
706 special_lj[3] = 1.0;
707 special_coul[1] = 0.0;
708 special_coul[2] = 1.0;
709 special_coul[3] = 1.0;
710 iarg += 1;
711 } else if (strcmp(arg[iarg],"lj/coul") == 0) {
712 if (iarg+4 > narg) error->all(FLERR,"Illegal special_bonds command");
713 special_lj[1] = special_coul[1] = atof(arg[iarg+1]);
714 special_lj[2] = special_coul[2] = atof(arg[iarg+2]);
715 special_lj[3] = special_coul[3] = atof(arg[iarg+3]);
716 iarg += 4;
717 } else if (strcmp(arg[iarg],"lj") == 0) {
718 if (iarg+4 > narg) error->all(FLERR,"Illegal special_bonds command");
719 special_lj[1] = atof(arg[iarg+1]);
720 special_lj[2] = atof(arg[iarg+2]);
721 special_lj[3] = atof(arg[iarg+3]);
722 iarg += 4;
723 } else if (strcmp(arg[iarg],"coul") == 0) {
724 if (iarg+4 > narg) error->all(FLERR,"Illegal special_bonds command");
725 special_coul[1] = atof(arg[iarg+1]);
726 special_coul[2] = atof(arg[iarg+2]);
727 special_coul[3] = atof(arg[iarg+3]);
728 iarg += 4;
729 } else if (strcmp(arg[iarg],"angle") == 0) {
730 if (iarg+2 > narg) error->all(FLERR,"Illegal special_bonds command");
731 if (strcmp(arg[iarg+1],"no") == 0) special_angle = 0;
732 else if (strcmp(arg[iarg+1],"yes") == 0) special_angle = 1;
733 else error->all(FLERR,"Illegal special_bonds command");
734 iarg += 2;
735 } else if (strcmp(arg[iarg],"dihedral") == 0) {
736 if (iarg+2 > narg) error->all(FLERR,"Illegal special_bonds command");
737 if (strcmp(arg[iarg+1],"no") == 0) special_dihedral = 0;
738 else if (strcmp(arg[iarg+1],"yes") == 0) special_dihedral = 1;
739 else error->all(FLERR,"Illegal special_bonds command");
740 iarg += 2;
741 } else if (strcmp(arg[iarg],"extra") == 0) {
742 if (iarg+2 > narg) error->all(FLERR,"Illegal special_bonds command");
743 special_extra = atoi(arg[iarg+1]);
744 iarg += 2;
745 } else error->all(FLERR,"Illegal special_bonds command");
746 }
747
748 for (int i = 1; i <= 3; i++)
749 if (special_lj[i] < 0.0 || special_lj[i] > 1.0 ||
750 special_coul[i] < 0.0 || special_coul[i] > 1.0)
751 error->all(FLERR,"Illegal special_bonds command");
752
753 if (special_extra < 0) error->all(FLERR,"Illegal special_bonds command");
754 }
755
756 /* ----------------------------------------------------------------------
757 compute bounds implied by numeric str with a possible wildcard asterik
758 1 = lower bound, nmax = upper bound
759 5 possibilities:
760 (1) i = i to i, (2) * = nmin to nmax,
761 (3) i* = i to nmax, (4) *j = nmin to j, (5) i*j = i to j
762 return nlo,nhi
763 ------------------------------------------------------------------------- */
764
bounds(char * str,int nmax,int & nlo,int & nhi,int nmin)765 void Force::bounds(char *str, int nmax, int &nlo, int &nhi, int nmin)
766 {
767 char *ptr = strchr(str,'*');
768
769 if (ptr == NULL) {
770 nlo = nhi = atoi(str);
771 } else if (strlen(str) == 1) {
772 nlo = nmin;
773 nhi = nmax;
774 } else if (ptr == str) {
775 nlo = nmin;
776 nhi = atoi(ptr+1);
777 } else if (strlen(ptr+1) == 0) {
778 nlo = atoi(str);
779 nhi = nmax;
780 } else {
781 nlo = atoi(str);
782 nhi = atoi(ptr+1);
783 }
784
785 if (nlo < nmin || nhi > nmax)
786 error->all(FLERR,"Numeric index is out of bounds");
787 }
788
789 /* ----------------------------------------------------------------------
790 read a floating point value from a string
791 generate an error if not a legitimate floating point value
792 called by various commands to check validity of their arguments
793 ------------------------------------------------------------------------- */
794
numeric(const char * file,const int line,const char * const str)795 double Force::numeric(const char *file, const int line, const char *const str)
796 {
797 const unsigned int n = strlen(str);
798 char *dstr;
799 dstr = new char[n+1];
800 for (unsigned int i = 0; i < n; i++)
801 {
802 if (isdigit(str[i]) ||
803 str[i] == '-' ||
804 str[i] == '+' ||
805 str[i] == '.' ||
806 str[i] == 'e' ||
807 str[i] == 'E')
808 dstr[i] = str[i];
809 else if (str[i] == '\r' && i == n-1)
810 dstr[i] = '\0';
811 else
812 error->all(file, line, "Expected floating point parameter in input script or data file");
813 }
814 dstr[n] = '\0';
815
816 const double val = atof(dstr);
817 delete [] dstr;
818 return val;
819 }
820
821 /* ----------------------------------------------------------------------
822 read an integer value from a string
823 generate an error if not a legitimate integer value
824 called by various commands to check validity of their arguments
825 ------------------------------------------------------------------------- */
826
inumeric(const char * file,const int line,const char * const str)827 int Force::inumeric(const char *file, const int line, const char *const str)
828 {
829 const unsigned int n = strlen(str);
830 char *istr;
831 istr = new char[n+1];
832 for (unsigned int i = 0; i < n; i++)
833 {
834 if (isdigit(str[i]) || str[i] == '-' || str[i] == '+')
835 istr[i] = str[i];
836 else if (str[i] == '\r' && i == n-1)
837 istr[i] = '\0';
838 else
839 error->all(file, line, "Expected integer parameter in input script or data file");
840 }
841 istr[n] = '\0';
842
843 const int val = atoi(istr);
844 delete [] istr;
845 return val;
846 }
847
848 /* ----------------------------------------------------------------------
849 memory usage of force classes
850 ------------------------------------------------------------------------- */
851
memory_usage()852 bigint Force::memory_usage()
853 {
854 bigint bytes = 0;
855 if (pair) bytes += static_cast<bigint> (pair->memory_usage());
856 if (bond) bytes += static_cast<bigint> (bond->memory_usage());
857 if (angle) bytes += static_cast<bigint> (angle->memory_usage());
858 if (dihedral) bytes += static_cast<bigint> (dihedral->memory_usage());
859 if (improper) bytes += static_cast<bigint> (improper->memory_usage());
860 if (kspace) bytes += static_cast<bigint> (kspace->memory_usage());
861 return bytes;
862 }
863