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 "molecule.h"
16
17 #include "atom.h"
18 #include "atom_vec.h"
19 #include "atom_vec_body.h"
20 #include "comm.h"
21 #include "domain.h"
22 #include "error.h"
23 #include "force.h"
24 #include "math_extra.h"
25 #include "math_eigen.h"
26 #include "memory.h"
27 #include "tokenizer.h"
28
29 #include <cmath>
30 #include <cstring>
31
32 using namespace LAMMPS_NS;
33
34 #define MAXLINE 256
35 #define EPSILON 1.0e-7
36 #define BIG 1.0e20
37
38 #define SINERTIA 0.4 // moment of inertia prefactor for sphere
39
40 /* ---------------------------------------------------------------------- */
41
Molecule(LAMMPS * lmp,int narg,char ** arg,int & index)42 Molecule::Molecule(LAMMPS *lmp, int narg, char **arg, int &index) :
43 Pointers(lmp), id(nullptr), x(nullptr), type(nullptr), molecule(nullptr), q(nullptr), radius(nullptr),
44 rmass(nullptr), num_bond(nullptr), bond_type(nullptr), bond_atom(nullptr),
45 num_angle(nullptr), angle_type(nullptr), angle_atom1(nullptr), angle_atom2(nullptr),
46 angle_atom3(nullptr), num_dihedral(nullptr), dihedral_type(nullptr), dihedral_atom1(nullptr),
47 dihedral_atom2(nullptr), dihedral_atom3(nullptr), dihedral_atom4(nullptr), num_improper(nullptr),
48 improper_type(nullptr), improper_atom1(nullptr), improper_atom2(nullptr),
49 improper_atom3(nullptr), improper_atom4(nullptr), nspecial(nullptr), special(nullptr),
50 shake_flag(nullptr), shake_atom(nullptr), shake_type(nullptr), avec_body(nullptr), ibodyparams(nullptr),
51 dbodyparams(nullptr), fragmentmask(nullptr),
52 dx(nullptr), dxcom(nullptr), dxbody(nullptr), quat_external(nullptr), fp(nullptr), count(nullptr)
53 {
54 me = comm->me;
55
56 if (index >= narg) error->all(FLERR,"Illegal molecule command");
57
58 id = utils::strdup(arg[0]);
59 if (!utils::is_id(id))
60 error->all(FLERR,"Molecule template ID must have only "
61 "alphanumeric or underscore characters");
62
63 // parse args until reach unknown arg (next file)
64
65 toffset = 0;
66 boffset = aoffset = doffset = ioffset = 0;
67 sizescale = 1.0;
68
69 int ifile = index;
70 int iarg = ifile+1;
71
72 while (iarg < narg) {
73 if (strcmp(arg[iarg],"offset") == 0) {
74 if (iarg+6 > narg) error->all(FLERR,"Illegal molecule command");
75 toffset = utils::inumeric(FLERR,arg[iarg+1],false,lmp);
76 boffset = utils::inumeric(FLERR,arg[iarg+2],false,lmp);
77 aoffset = utils::inumeric(FLERR,arg[iarg+3],false,lmp);
78 doffset = utils::inumeric(FLERR,arg[iarg+4],false,lmp);
79 ioffset = utils::inumeric(FLERR,arg[iarg+5],false,lmp);
80 if (toffset < 0 || boffset < 0 || aoffset < 0 ||
81 doffset < 0 || ioffset < 0)
82 error->all(FLERR,"Illegal molecule command");
83 iarg += 6;
84 } else if (strcmp(arg[iarg],"toff") == 0) {
85 if (iarg+2 > narg) error->all(FLERR,"Illegal molecule command");
86 toffset = utils::inumeric(FLERR,arg[iarg+1],false,lmp);
87 if (toffset < 0) error->all(FLERR,"Illegal molecule command");
88 iarg += 2;
89 } else if (strcmp(arg[iarg],"boff") == 0) {
90 if (iarg+2 > narg) error->all(FLERR,"Illegal molecule command");
91 boffset = utils::inumeric(FLERR,arg[iarg+1],false,lmp);
92 if (boffset < 0) error->all(FLERR,"Illegal molecule command");
93 iarg += 2;
94 } else if (strcmp(arg[iarg],"aoff") == 0) {
95 if (iarg+2 > narg) error->all(FLERR,"Illegal molecule command");
96 aoffset = utils::inumeric(FLERR,arg[iarg+1],false,lmp);
97 if (aoffset < 0) error->all(FLERR,"Illegal molecule command");
98 iarg += 2;
99 } else if (strcmp(arg[iarg],"doff") == 0) {
100 if (iarg+2 > narg) error->all(FLERR,"Illegal molecule command");
101 doffset = utils::inumeric(FLERR,arg[iarg+1],false,lmp);
102 if (doffset < 0) error->all(FLERR,"Illegal molecule command");
103 iarg += 2;
104 } else if (strcmp(arg[iarg],"ioff") == 0) {
105 if (iarg+2 > narg) error->all(FLERR,"Illegal molecule command");
106 ioffset = utils::inumeric(FLERR,arg[iarg+1],false,lmp);
107 if (ioffset < 0) error->all(FLERR,"Illegal molecule command");
108 iarg += 2;
109 } else if (strcmp(arg[iarg],"scale") == 0) {
110 if (iarg+2 > narg) error->all(FLERR,"Illegal molecule command");
111 sizescale = utils::numeric(FLERR,arg[iarg+1],false,lmp);
112 if (sizescale <= 0.0) error->all(FLERR,"Illegal molecule command");
113 iarg += 2;
114 } else break;
115 }
116
117 index = iarg;
118
119 // last molecule if have scanned all args
120
121 if (iarg == narg) last = 1;
122 else last = 0;
123
124 // initialize all fields to empty
125
126 initialize();
127
128 // scan file for sizes of all fields and allocate storage for them
129
130 if (me == 0) {
131 fp = fopen(arg[ifile],"r");
132 if (fp == nullptr)
133 error->one(FLERR,"Cannot open molecule file {}: {}",
134 arg[ifile], utils::getsyserror());
135 }
136 read(0);
137 if (me == 0) fclose(fp);
138 allocate();
139
140 // read file again to populate all fields
141
142 if (me == 0) fp = fopen(arg[ifile],"r");
143 read(1);
144 if (me == 0) fclose(fp);
145
146 // stats
147
148 if (me == 0)
149 utils::logmesg(lmp,"Read molecule template {}:\n {} molecules\n"
150 " {} fragments\n"
151 " {} atoms with max type {}\n"
152 " {} bonds with max type {}\n"
153 " {} angles with max type {}\n"
154 " {} dihedrals with max type {}\n"
155 " {} impropers with max type {}\n", id,nmolecules,
156 nfragments,natoms,ntypes,nbonds,nbondtypes,nangles,
157 nangletypes,ndihedrals,ndihedraltypes,nimpropers,
158 nimpropertypes);
159 }
160
161 /* ---------------------------------------------------------------------- */
162
~Molecule()163 Molecule::~Molecule()
164 {
165 delete [] id;
166 deallocate();
167 }
168
169 /* ----------------------------------------------------------------------
170 compute center = geometric center of molecule
171 also compute:
172 dx = displacement of each atom from center
173 molradius = radius of molecule from center
174 including finite-size particles or body particles
175 ------------------------------------------------------------------------- */
176
compute_center()177 void Molecule::compute_center()
178 {
179 if (centerflag) return;
180 centerflag = 1;
181
182 center[0] = center[1] = center[2] = 0.0;
183 for (int i = 0; i < natoms; i++) {
184 center[0] += x[i][0];
185 center[1] += x[i][1];
186 center[2] += x[i][2];
187 }
188 center[0] /= natoms;
189 center[1] /= natoms;
190 center[2] /= natoms;
191
192 memory->destroy(dx);
193 memory->create(dx,natoms,3,"molecule:dx");
194
195 for (int i = 0; i < natoms; i++) {
196 dx[i][0] = x[i][0] - center[0];
197 dx[i][1] = x[i][1] - center[1];
198 dx[i][2] = x[i][2] - center[2];
199 }
200
201 molradius = 0.0;
202 for (int i = 0; i < natoms; i++) {
203 double rad = MathExtra::len3(dx[i]);
204 if (radiusflag) rad += radius[i];
205 molradius = MAX(molradius,rad);
206 }
207 }
208
209 /* ----------------------------------------------------------------------
210 compute masstotal = total mass of molecule
211 could have been set by user, otherwise calculate it
212 ------------------------------------------------------------------------- */
213
compute_mass()214 void Molecule::compute_mass()
215 {
216 if (massflag) return;
217 massflag = 1;
218
219 atom->check_mass(FLERR);
220
221 masstotal = 0.0;
222 for (int i = 0; i < natoms; i++) {
223 if (rmassflag) masstotal += rmass[i];
224 else masstotal += atom->mass[type[i]];
225 }
226 }
227
228 /* ----------------------------------------------------------------------
229 compute com = center of mass of molecule
230 could have been set by user, otherwise calculate it
231 works for finite size particles assuming no overlap
232 also compute:
233 dxcom = displacement of each atom from COM
234 comatom = which atom (1-Natom) is nearest the COM
235 maxextent = furthest any atom in molecule is from comatom (not COM)
236 ------------------------------------------------------------------------- */
237
compute_com()238 void Molecule::compute_com()
239 {
240 if (!comflag) {
241 comflag = 1;
242
243 atom->check_mass(FLERR);
244
245 double onemass;
246 com[0] = com[1] = com[2] = 0.0;
247 for (int i = 0; i < natoms; i++) {
248 if (rmassflag) onemass = rmass[i];
249 else onemass = atom->mass[type[i]];
250 com[0] += x[i][0]*onemass;
251 com[1] += x[i][1]*onemass;
252 com[2] += x[i][2]*onemass;
253 }
254 if (masstotal > 0.0) {
255 com[0] /= masstotal;
256 com[1] /= masstotal;
257 com[2] /= masstotal;
258 }
259 }
260
261 memory->destroy(dxcom);
262 memory->create(dxcom,natoms,3,"molecule:dxcom");
263
264 for (int i = 0; i < natoms; i++) {
265 dxcom[i][0] = x[i][0] - com[0];
266 dxcom[i][1] = x[i][1] - com[1];
267 dxcom[i][2] = x[i][2] - com[2];
268 }
269
270 double rsqmin = BIG;
271 for (int i = 0; i < natoms; i++) {
272 double rsq = MathExtra::lensq3(dxcom[i]);
273 if (rsq < rsqmin) {
274 comatom = i;
275 rsqmin = rsq;
276 }
277 }
278
279 double rsqmax = 0.0;
280 for (int i = 0; i < natoms; i++) {
281 double dx = x[comatom][0] - x[i][0];
282 double dy = x[comatom][1] - x[i][1];
283 double dz = x[comatom][2] - x[i][2];
284 double rsq = dx*dx + dy*dy + dz*dz;
285 rsqmax = MAX(rsqmax,rsq);
286 }
287
288 comatom++;
289 maxextent = sqrt(rsqmax);
290 }
291
292 /* ----------------------------------------------------------------------
293 compute itensor = 6 moments of inertia of molecule around xyz axes
294 could have been set by user, otherwise calculate it
295 accounts for finite size spheres, assuming no overlap
296 also compute:
297 inertia = 3 principal components of inertia
298 ex,ey,ez = principal axes in space coords
299 quat = quaternion for orientation of molecule
300 dxbody = displacement of each atom from COM in body frame
301 ------------------------------------------------------------------------- */
302
compute_inertia()303 void Molecule::compute_inertia()
304 {
305 if (!inertiaflag) {
306 inertiaflag = 1;
307
308 atom->check_mass(FLERR);
309
310 double onemass,dx,dy,dz;
311 for (int i = 0; i < 6; i++) itensor[i] = 0.0;
312 for (int i = 0; i < natoms; i++) {
313 if (rmassflag) onemass = rmass[i];
314 else onemass = atom->mass[type[i]];
315 dx = dxcom[i][0];
316 dy = dxcom[i][1];
317 dz = dxcom[i][2];
318 itensor[0] += onemass * (dy*dy + dz*dz);
319 itensor[1] += onemass * (dx*dx + dz*dz);
320 itensor[2] += onemass * (dx*dx + dy*dy);
321 itensor[3] -= onemass * dy*dz;
322 itensor[4] -= onemass * dx*dz;
323 itensor[5] -= onemass * dx*dy;
324 }
325
326 if (radiusflag) {
327 for (int i = 0; i < natoms; i++) {
328 if (rmassflag) onemass = rmass[i];
329 else onemass = atom->mass[type[i]];
330 itensor[0] += SINERTIA*onemass * radius[i]*radius[i];
331 itensor[1] += SINERTIA*onemass * radius[i]*radius[i];
332 itensor[2] += SINERTIA*onemass * radius[i]*radius[i];
333 }
334 }
335 }
336
337 // diagonalize inertia tensor for each body via Jacobi rotations
338 // inertia = 3 eigenvalues = principal moments of inertia
339 // evectors and exzy = 3 evectors = principal axes of rigid body
340
341 double cross[3];
342 double tensor[3][3],evectors[3][3];
343
344 tensor[0][0] = itensor[0];
345 tensor[1][1] = itensor[1];
346 tensor[2][2] = itensor[2];
347 tensor[1][2] = tensor[2][1] = itensor[3];
348 tensor[0][2] = tensor[2][0] = itensor[4];
349 tensor[0][1] = tensor[1][0] = itensor[5];
350
351 if (MathEigen::jacobi3(tensor,inertia,evectors))
352 error->all(FLERR,"Insufficient Jacobi rotations for rigid molecule");
353
354 ex[0] = evectors[0][0];
355 ex[1] = evectors[1][0];
356 ex[2] = evectors[2][0];
357 ey[0] = evectors[0][1];
358 ey[1] = evectors[1][1];
359 ey[2] = evectors[2][1];
360 ez[0] = evectors[0][2];
361 ez[1] = evectors[1][2];
362 ez[2] = evectors[2][2];
363
364 // if any principal moment < scaled EPSILON, set to 0.0
365
366 double max;
367 max = MAX(inertia[0],inertia[1]);
368 max = MAX(max,inertia[2]);
369
370 if (inertia[0] < EPSILON*max) inertia[0] = 0.0;
371 if (inertia[1] < EPSILON*max) inertia[1] = 0.0;
372 if (inertia[2] < EPSILON*max) inertia[2] = 0.0;
373
374 // enforce 3 evectors as a right-handed coordinate system
375 // flip 3rd vector if needed
376
377 MathExtra::cross3(ex,ey,cross);
378 if (MathExtra::dot3(cross,ez) < 0.0) MathExtra::negate3(ez);
379
380 // create quaternion
381
382 MathExtra::exyz_to_q(ex,ey,ez,quat);
383
384 // compute displacements in body frame defined by quat
385
386 memory->destroy(dxbody);
387 memory->create(dxbody,natoms,3,"molecule:dxbody");
388
389 for (int i = 0; i < natoms; i++)
390 MathExtra::transpose_matvec(ex,ey,ez,dxcom[i],dxbody[i]);
391 }
392
393 /* ----------------------------------------------------------------------
394 read molecule info from file
395 flag = 0, just scan for sizes of fields
396 flag = 1, read and store fields
397 ------------------------------------------------------------------------- */
398
read(int flag)399 void Molecule::read(int flag)
400 {
401 char line[MAXLINE];
402 char *eof;
403
404 // skip 1st line of file
405
406 if (me == 0) {
407 eof = fgets(line,MAXLINE,fp);
408 if (eof == nullptr) error->one(FLERR,"Unexpected end of molecule file");
409 }
410
411 // read header lines
412 // skip blank lines or lines that start with "#"
413 // stop when read an unrecognized line
414
415 while (1) {
416
417 readline(line);
418
419 // trim comments. if line is blank, continue
420
421 auto text = utils::trim(utils::trim_comment(line));
422 if (text.empty()) continue;
423
424 // search line for header keywords and set corresponding variable
425 try {
426 ValueTokenizer values(text);
427
428 int nmatch = values.count();
429 int nwant = 0;
430 if (values.contains("atoms")) {
431 natoms = values.next_int();
432 nwant = 2;
433 } else if (values.contains("bonds")) {
434 nbonds = values.next_int();
435 nwant = 2;
436 } else if (values.contains("angles")) {
437 nangles = values.next_int();
438 nwant = 2;
439 } else if (values.contains("dihedrals")) {
440 ndihedrals = values.next_int();
441 nwant = 2;
442 } else if (values.contains("impropers")) {
443 nimpropers = values.next_int();
444 nwant = 2;
445 } else if (values.contains("fragments")) {
446 nfragments = values.next_int();
447 nwant = 2;
448 } else if (values.contains("mass")) {
449 massflag = 1;
450 masstotal = values.next_double();
451 nwant = 2;
452 masstotal *= sizescale*sizescale*sizescale;
453 } else if (values.contains("com")) {
454 comflag = 1;
455 com[0] = values.next_double();
456 com[1] = values.next_double();
457 com[2] = values.next_double();
458 nwant = 4;
459 com[0] *= sizescale;
460 com[1] *= sizescale;
461 com[2] *= sizescale;
462 if (domain->dimension == 2 && com[2] != 0.0)
463 error->all(FLERR,"Molecule file z center-of-mass must be 0.0 for 2d");
464 } else if (values.contains("inertia")) {
465 inertiaflag = 1;
466 itensor[0] = values.next_double();
467 itensor[1] = values.next_double();
468 itensor[2] = values.next_double();
469 itensor[3] = values.next_double();
470 itensor[4] = values.next_double();
471 itensor[5] = values.next_double();
472 nwant = 7;
473 const double scale5 = sizescale*sizescale*sizescale*sizescale*sizescale;
474 itensor[0] *= scale5;
475 itensor[1] *= scale5;
476 itensor[2] *= scale5;
477 itensor[3] *= scale5;
478 itensor[4] *= scale5;
479 itensor[5] *= scale5;
480 } else if (values.contains("body")) {
481 bodyflag = 1;
482 avec_body = (AtomVecBody *) atom->style_match("body");
483 if (!avec_body)
484 error->all(FLERR,"Molecule file requires atom style body");
485 nibody = values.next_int();
486 ndbody = values.next_int();
487 nwant = 3;
488 } else {
489 // unknown header keyword
490 if (utils::strmatch(text,"^\\d+\\s+\\S+")) {
491 values.next_int();
492 auto keyword = values.next_string();
493 error->one(FLERR,"Invalid header keyword: {}",keyword);
494 } else break;
495 }
496 if (nmatch != nwant)
497 error->one(FLERR,"Invalid header line format in molecule file");
498 } catch (TokenizerException &e) {
499 error->one(FLERR, "Invalid header in molecule file\n"
500 "{}", e.what());
501 }
502 }
503
504 // error checks
505
506 if (natoms < 1)
507 error->all(FLERR,"No atoms or invalid atom count in molecule file");
508 if (nbonds < 0) error->all(FLERR,"Invalid bond count in molecule file");
509 if (nangles < 0) error->all(FLERR,"Invalid angle count in molecule file");
510 if (ndihedrals < 0)
511 error->all(FLERR,"Invalid dihedral count in molecule file");
512 if (nimpropers < 0)
513 error->all(FLERR,"Invalid improper count in molecule file");
514
515 // count = vector for tallying bonds,angles,etc per atom
516
517 if (flag == 0) memory->create(count,natoms,"molecule:count");
518
519 // grab keyword and skip next line
520
521 std::string keyword = parse_keyword(0,line);
522 readline(line);
523
524 // loop over sections of molecule file
525
526 while (!keyword.empty()) {
527 if (keyword == "Coords") {
528 xflag = 1;
529 if (flag) coords(line);
530 else skip_lines(natoms,line,keyword);
531 } else if (keyword == "Types") {
532 typeflag = 1;
533 if (flag) types(line);
534 else skip_lines(natoms,line,keyword);
535 } else if (keyword == "Molecules") {
536 moleculeflag = 1;
537 if (flag) molecules(line);
538 else skip_lines(natoms,line,keyword);
539 } else if (keyword == "Fragments") {
540 if (nfragments == 0)
541 error->all(FLERR,"Found Fragments section but no nfragments setting in header");
542 fragmentflag = 1;
543 if (flag) fragments(line);
544 else skip_lines(nfragments,line,keyword);
545 } else if (keyword == "Charges") {
546 qflag = 1;
547 if (flag) charges(line);
548 else skip_lines(natoms,line,keyword);
549 } else if (keyword == "Diameters") {
550 radiusflag = 1;
551 if (flag) diameters(line);
552 else skip_lines(natoms,line,keyword);
553 } else if (keyword == "Masses") {
554 rmassflag = 1;
555 if (flag) masses(line);
556 else skip_lines(natoms,line,keyword);
557
558 } else if (keyword == "Bonds") {
559 if (nbonds == 0)
560 error->all(FLERR,"Found Bonds section but no nbonds setting in header");
561 bondflag = tag_require = 1;
562 bonds(flag,line);
563 } else if (keyword == "Angles") {
564 if (nangles == 0)
565 error->all(FLERR,"Found Angles section but no nangles setting in header");
566 angleflag = tag_require = 1;
567 angles(flag,line);
568 } else if (keyword == "Dihedrals") {
569 if (ndihedrals == 0) error->all(FLERR,"Found Dihedrals section"
570 "but no ndihedrals setting in header");
571 dihedralflag = tag_require = 1;
572 dihedrals(flag,line);
573 } else if (keyword == "Impropers") {
574 if (nimpropers == 0) error->all(FLERR,"Found Impropers section"
575 "but no nimpropers setting in header");
576 improperflag = tag_require = 1;
577 impropers(flag,line);
578
579 } else if (keyword == "Special Bond Counts") {
580 nspecialflag = 1;
581 nspecial_read(flag,line);
582 } else if (keyword == "Special Bonds") {
583 specialflag = tag_require = 1;
584 if (flag) special_read(line);
585 else skip_lines(natoms,line,keyword);
586
587 } else if (keyword == "Shake Flags") {
588 shakeflagflag = 1;
589 if (flag) shakeflag_read(line);
590 else skip_lines(natoms,line,keyword);
591 } else if (keyword == "Shake Atoms") {
592 shakeatomflag = tag_require = 1;
593 if (shaketypeflag) shakeflag = 1;
594 if (!shakeflagflag)
595 error->all(FLERR,"Shake Flags section must come before Shake Atoms section");
596 if (flag) shakeatom_read(line);
597 else skip_lines(natoms,line,keyword);
598 } else if (keyword == "Shake Bond Types") {
599 shaketypeflag = 1;
600 if (shakeatomflag) shakeflag = 1;
601 if (!shakeflagflag)
602 error->all(FLERR,"Shake Flags section must come before Shake Bonds section");
603 if (flag) shaketype_read(line);
604 else skip_lines(natoms,line,keyword);
605
606 } else if (keyword == "Body Integers") {
607 if (bodyflag == 0 || nibody == 0)
608 error->all(FLERR,"Found Body Integers section but no setting in header");
609 ibodyflag = 1;
610 body(flag,0,line);
611 } else if (keyword == "Body Doubles") {
612 if (bodyflag == 0 || ndbody == 0)
613 error->all(FLERR,"Found Body Doubles section but no setting in header");
614 dbodyflag = 1;
615 body(flag,1,line);
616 } else {
617
618 // Error: Either a too long/short section or a typo in the keyword
619
620 if (utils::strmatch(keyword,"^[A-Za-z ]+$"))
621 error->one(FLERR,"Unknown section '{}' in molecule "
622 "file\n",keyword);
623 else error->one(FLERR,"Unexpected line in molecule file "
624 "while looking for the next "
625 "section:\n{}",line);
626 }
627 keyword = parse_keyword(1,line);
628 }
629
630 // error check
631
632 if (flag == 0) {
633 if ((nspecialflag && !specialflag) || (!nspecialflag && specialflag))
634 error->all(FLERR,"Molecule file needs both Special Bond sections");
635 if (specialflag && !bondflag)
636 error->all(FLERR,"Molecule file has special flags but no bonds");
637 if ((shakeflagflag || shakeatomflag || shaketypeflag) && !shakeflag)
638 error->all(FLERR,"Molecule file shake info is incomplete");
639 if (bodyflag && nibody && ibodyflag == 0)
640 error->all(FLERR,"Molecule file has no Body Integers section");
641 if (bodyflag && ndbody && dbodyflag == 0)
642 error->all(FLERR,"Molecule file has no Body Doubles section");
643 if (nfragments > 0 && !fragmentflag)
644 error->all(FLERR,"Molecule file has no Fragments section");
645 }
646
647 // auto-generate special bonds if needed and not in file
648
649 if (bondflag && specialflag == 0) {
650 if (domain->box_exist == 0)
651 error->all(FLERR,"Cannot auto-generate special bonds before "
652 "simulation box is defined");
653
654 if (flag) {
655 special_generate();
656 specialflag = 1;
657 nspecialflag = 1;
658 }
659 }
660
661 // body particle must have natom = 1
662 // set radius by having body class compute its own radius
663
664 if (bodyflag) {
665 radiusflag = 1;
666 if (natoms != 1)
667 error->all(FLERR,"Molecule natoms must be 1 for body particle");
668 if (sizescale != 1.0)
669 error->all(FLERR,"Molecule sizescale must be 1.0 for body particle");
670 if (flag) {
671 radius[0] = avec_body->radius_body(nibody,ndbody,ibodyparams,dbodyparams);
672 maxradius = radius[0];
673 }
674 }
675
676 // clean up
677
678 if (flag) memory->destroy(count);
679 }
680
681 /* ----------------------------------------------------------------------
682 read coords from file
683 ------------------------------------------------------------------------- */
684
coords(char * line)685 void Molecule::coords(char *line)
686 {
687 for (int i = 0; i < natoms; i++) count[i] = 0;
688 try {
689 for (int i = 0; i < natoms; i++) {
690 readline(line);
691
692 ValueTokenizer values(utils::trim_comment(line));
693 if (values.count() != 4)
694 error->all(FLERR,"Invalid line in Coords section of "
695 "molecule file: {}",line);
696
697 int iatom = values.next_int() - 1;
698 if (iatom < 0 || iatom >= natoms)
699 error->all(FLERR,"Invalid atom index in Coords section of molecule file");
700 count[iatom]++;
701 x[iatom][0] = values.next_double();
702 x[iatom][1] = values.next_double();
703 x[iatom][2] = values.next_double();
704
705 x[iatom][0] *= sizescale;
706 x[iatom][1] *= sizescale;
707 x[iatom][2] *= sizescale;
708 }
709 } catch (TokenizerException &e) {
710 error->all(FLERR,"Invalid line in Coords section of "
711 "molecule file: {}\n{}",e.what(),line);
712 }
713
714 for (int i = 0; i < natoms; i++)
715 if (count[i] == 0) error->all(FLERR,"Atom {} missing in Coords "
716 "section of molecule file",i+1);
717
718 if (domain->dimension == 2) {
719 for (int i = 0; i < natoms; i++)
720 if (x[i][2] != 0.0)
721 error->all(FLERR,"Z coord in molecule file for atom {} "
722 "must be 0.0 for 2d-simulation.",i+1);
723 }
724 }
725
726 /* ----------------------------------------------------------------------
727 read types from file
728 set ntypes = max of any atom type
729 ------------------------------------------------------------------------- */
730
types(char * line)731 void Molecule::types(char *line)
732 {
733 for (int i = 0; i < natoms; i++) count[i] = 0;
734 try {
735 for (int i = 0; i < natoms; i++) {
736 readline(line);
737
738 ValueTokenizer values(utils::trim_comment(line));
739 if (values.count() != 2)
740 error->all(FLERR,"Invalid line in Types section of "
741 "molecule file: {}",line);
742
743 int iatom = values.next_int() - 1;
744 if (iatom < 0 || iatom >= natoms)
745 error->all(FLERR,"Invalid atom index in Types section of molecule file");
746 count[iatom]++;
747 type[iatom] = values.next_int();
748 type[iatom] += toffset;
749 }
750 } catch (TokenizerException &e) {
751 error->all(FLERR, "Invalid line in Types section of "
752 "molecule file: {}\n{}", e.what(),line);
753 }
754
755 for (int i = 0; i < natoms; i++) {
756 if (count[i] == 0) error->all(FLERR,"Atom {} missing in Types "
757 "section of molecule file",i+1);
758
759 if ((type[i] <= 0) || (domain->box_exist && (type[i] > atom->ntypes)))
760 error->all(FLERR,"Invalid atom type {} for atom {} "
761 "in molecule file",type[i],i+1);
762
763 ntypes = MAX(ntypes,type[i]);
764 }
765 }
766
767 /* ----------------------------------------------------------------------
768 read molecules from file
769 set nmolecules = max of any molecule type
770 ------------------------------------------------------------------------- */
771
molecules(char * line)772 void Molecule::molecules(char *line)
773 {
774 for (int i = 0; i < natoms; i++) count[i] = 0;
775 try {
776 for (int i = 0; i < natoms; i++) {
777 readline(line);
778 ValueTokenizer values(utils::trim_comment(line));
779 if (values.count() != 2)
780 error->all(FLERR,"Invalid line in Molecules section of "
781 "molecule file: {}",line);
782
783 int iatom = values.next_int() - 1;
784 if (iatom < 0 || iatom >= natoms)
785 error->all(FLERR,"Invalid atom index in Molecules section of molecule file");
786 count[iatom]++;
787 molecule[iatom] = values.next_tagint();
788 // molecule[iatom] += moffset; // placeholder for possible molecule offset
789 }
790 } catch (TokenizerException &e) {
791 error->all(FLERR, "Invalid line in Molecules section of "
792 "molecule file: {}\n{}",e.what(),line);
793 }
794
795 for (int i = 0; i < natoms; i++)
796 if (count[i] == 0) error->all(FLERR,"Atom {} missing in Molecules "
797 "section of molecule file",i+1);
798
799 for (int i = 0; i < natoms; i++)
800 if (molecule[i] < 0)
801 error->all(FLERR,"Invalid molecule ID {} for atom {} "
802 "in molecule file",molecule[i],i+1);
803
804 for (int i = 0; i < natoms; i++)
805 nmolecules = MAX(nmolecules,molecule[i]);
806 }
807
808 /* ----------------------------------------------------------------------
809 read fragments from file
810 ------------------------------------------------------------------------- */
811
fragments(char * line)812 void Molecule::fragments(char *line)
813 {
814 try {
815 for (int i = 0; i < nfragments; i++) {
816 readline(line);
817
818 ValueTokenizer values(utils::trim_comment(line));
819
820 if ((int)values.count() > natoms+1)
821 error->all(FLERR,"Too many atoms per fragment in Fragments "
822 "section of molecule file");
823
824 fragmentnames[i] = values.next_string();
825
826 while (values.has_next()) {
827 int iatom = values.next_int()-1;
828 if (iatom < 0 || iatom >= natoms)
829 error->all(FLERR,"Invalid atom ID {} for fragment {} in "
830 "Fragments section of molecule file",
831 iatom+1, fragmentnames[i]);
832 fragmentmask[i][iatom] = 1;
833 }
834 }
835 } catch (TokenizerException &e) {
836 error->all(FLERR, "Invalid atom ID in Fragments section of "
837 "molecule file: {}\n{}", e.what(),line);
838 }
839 }
840
841 /* ----------------------------------------------------------------------
842 read charges from file
843 ------------------------------------------------------------------------- */
844
charges(char * line)845 void Molecule::charges(char *line)
846 {
847 for (int i = 0; i < natoms; i++) count[i] = 0;
848 try {
849 for (int i = 0; i < natoms; i++) {
850 readline(line);
851
852 ValueTokenizer values(utils::trim_comment(line));
853 if ((int)values.count() != 2)
854 error->all(FLERR,"Invalid line in Charges section of "
855 "molecule file: {}",line);
856
857 int iatom = values.next_int() - 1;
858 if (iatom < 0 || iatom >= natoms)
859 error->all(FLERR,"Invalid atom index in Charges section of molecule file");
860
861 count[iatom]++;
862 q[iatom] = values.next_double();
863 }
864 } catch (TokenizerException &e) {
865 error->all(FLERR, "Invalid line in Charges section of "
866 "molecule file: {}.\n{}",e.what(),line);
867 }
868
869 for (int i = 0; i < natoms; i++)
870 if (count[i] == 0) error->all(FLERR,"Atom {} missing in Charges "
871 "section of molecule file",i+1);
872 }
873
874 /* ----------------------------------------------------------------------
875 read diameters from file and set radii
876 ------------------------------------------------------------------------- */
877
diameters(char * line)878 void Molecule::diameters(char *line)
879 {
880 for (int i = 0; i < natoms; i++) count[i] = 0;
881 try {
882 maxradius = 0.0;
883 for (int i = 0; i < natoms; i++) {
884 readline(line);
885
886 ValueTokenizer values(utils::trim_comment(line));
887 if (values.count() != 2)
888 error->all(FLERR,"Invalid line in Diameters section of "
889 "molecule file: {}",line);
890 int iatom = values.next_int() - 1;
891 if (iatom < 0 || iatom >= natoms)
892 error->all(FLERR,"Invalid atom index in Diameters section of molecule file");
893 count[iatom]++;
894 radius[iatom] = values.next_double();
895 radius[iatom] *= sizescale;
896 radius[iatom] *= 0.5;
897 maxradius = MAX(maxradius,radius[iatom]);
898 }
899 } catch (TokenizerException &e) {
900 error->all(FLERR, "Invalid line in Diameters section of "
901 "molecule file: {}\n{}",e.what(),line);
902 }
903
904 for (int i = 0; i < natoms; i++) {
905 if (count[i] == 0) error->all(FLERR,"Atom {} missing in Diameters "
906 "section of molecule file",i+1);
907 if (radius[i] < 0.0)
908 error->all(FLERR,"Invalid atom diameter {} for atom {} "
909 "in molecule file", radius[i], i+1);
910 }
911 }
912
913 /* ----------------------------------------------------------------------
914 read masses from file
915 ------------------------------------------------------------------------- */
916
masses(char * line)917 void Molecule::masses(char *line)
918 {
919 for (int i = 0; i < natoms; i++) count[i] = 0;
920 try {
921 for (int i = 0; i < natoms; i++) {
922 readline(line);
923
924 ValueTokenizer values(utils::trim_comment(line));
925 if (values.count() != 2)
926 error->all(FLERR,"Invalid line in Masses section of "
927 "molecule file: {}",line);
928
929 int iatom = values.next_int() - 1;
930 if (iatom < 0 || iatom >= natoms)
931 error->all(FLERR,"Invalid atom index in Masses section of molecule file");
932 count[iatom]++;
933 rmass[iatom] = values.next_double();
934 rmass[iatom] *= sizescale*sizescale*sizescale;
935 }
936 } catch (TokenizerException &e) {
937 error->all(FLERR, "Invalid line in Masses section of "
938 "molecule file: {}\n{}",e.what(),line);
939 }
940
941 for (int i = 0; i < natoms; i++) {
942 if (count[i] == 0) error->all(FLERR,"Atom {} missing in Masses "
943 "section of molecule file",i+1);
944 if (rmass[i] <= 0.0)
945 error->all(FLERR,"Invalid atom mass {} for atom {} "
946 "in molecule file", radius[i], i+1);
947 }
948 }
949
950 /* ----------------------------------------------------------------------
951 read bonds from file
952 set nbondtypes = max type of any bond
953 store each with both atoms if newton_bond = 0
954 if flag = 0, just count bonds/atom
955 if flag = 1, store them with atoms
956 ------------------------------------------------------------------------- */
957
bonds(int flag,char * line)958 void Molecule::bonds(int flag, char *line)
959 {
960 int itype;
961 tagint m,atom1,atom2;
962 int newton_bond = force->newton_bond;
963
964 if (flag == 0)
965 for (int i = 0; i < natoms; i++) count[i] = 0;
966 else
967 for (int i = 0; i < natoms; i++) num_bond[i] = 0;
968
969 for (int i = 0; i < nbonds; i++) {
970 readline(line);
971
972 try {
973 ValueTokenizer values(utils::trim_comment(line));
974 if (values.count() != 4)
975 error->all(FLERR,"Invalid line in Bonds section of "
976 "molecule file: {}",line);
977 values.next_int();
978 itype = values.next_int();
979 atom1 = values.next_tagint();
980 atom2 = values.next_tagint();
981 } catch (TokenizerException &e) {
982 error->all(FLERR, "Invalid line in Bonds section of "
983 "molecule file: {}\n{}",e.what(),line);
984 }
985
986 itype += boffset;
987
988 if ((atom1 <= 0) || (atom1 > natoms) ||
989 (atom2 <= 0) || (atom2 > natoms) || (atom1 == atom2))
990 error->all(FLERR,"Invalid atom ID in Bonds section of molecule file");
991 if ((itype <= 0) || (domain->box_exist && (itype > atom->nbondtypes)))
992 error->all(FLERR,"Invalid bond type in Bonds section of molecule file");
993
994 if (flag) {
995 m = atom1-1;
996 nbondtypes = MAX(nbondtypes,itype);
997 bond_type[m][num_bond[m]] = itype;
998 bond_atom[m][num_bond[m]] = atom2;
999 num_bond[m]++;
1000 if (newton_bond == 0) {
1001 m = atom2-1;
1002 bond_type[m][num_bond[m]] = itype;
1003 bond_atom[m][num_bond[m]] = atom1;
1004 num_bond[m]++;
1005 }
1006 } else {
1007 count[atom1-1]++;
1008 if (newton_bond == 0) count[atom2-1]++;
1009 }
1010 }
1011
1012 // bond_per_atom = max of count vector
1013
1014 if (flag == 0) {
1015 bond_per_atom = 0;
1016 for (int i = 0; i < natoms; i++)
1017 bond_per_atom = MAX(bond_per_atom,count[i]);
1018 }
1019 }
1020
1021 /* ----------------------------------------------------------------------
1022 read angles from file
1023 store each with all 3 atoms if newton_bond = 0
1024 if flag = 0, just count angles/atom
1025 if flag = 1, store them with atoms
1026 ------------------------------------------------------------------------- */
1027
angles(int flag,char * line)1028 void Molecule::angles(int flag, char *line)
1029 {
1030 int itype;
1031 tagint m,atom1,atom2,atom3;
1032 int newton_bond = force->newton_bond;
1033
1034 if (flag == 0)
1035 for (int i = 0; i < natoms; i++) count[i] = 0;
1036 else
1037 for (int i = 0; i < natoms; i++) num_angle[i] = 0;
1038
1039 for (int i = 0; i < nangles; i++) {
1040 readline(line);
1041
1042 try {
1043 ValueTokenizer values(utils::trim_comment(line));
1044 if (values.count() != 5)
1045 error->all(FLERR,"Invalid line in Angles section of "
1046 "molecule file: {}",line);
1047 values.next_int();
1048 itype = values.next_int();
1049 atom1 = values.next_tagint();
1050 atom2 = values.next_tagint();
1051 atom3 = values.next_tagint();
1052 } catch (TokenizerException &e) {
1053 error->all(FLERR, "Invalid line in Angles section of "
1054 "molecule file: {}\n{}",e.what(),line);
1055 }
1056
1057 itype += aoffset;
1058
1059 if ((atom1 <= 0) || (atom1 > natoms) ||
1060 (atom2 <= 0) || (atom2 > natoms) ||
1061 (atom3 <= 0) || (atom3 > natoms) ||
1062 (atom1 == atom2) || (atom1 == atom3) || (atom2 == atom3))
1063 error->all(FLERR,"Invalid atom ID in Angles section of molecule file");
1064 if ((itype <= 0) || (domain->box_exist && (itype > atom->nangletypes)))
1065 error->all(FLERR,"Invalid angle type in Angles section of molecule file");
1066
1067 if (flag) {
1068 m = atom2-1;
1069 nangletypes = MAX(nangletypes,itype);
1070 angle_type[m][num_angle[m]] = itype;
1071 angle_atom1[m][num_angle[m]] = atom1;
1072 angle_atom2[m][num_angle[m]] = atom2;
1073 angle_atom3[m][num_angle[m]] = atom3;
1074 num_angle[m]++;
1075 if (newton_bond == 0) {
1076 m = atom1-1;
1077 angle_type[m][num_angle[m]] = itype;
1078 angle_atom1[m][num_angle[m]] = atom1;
1079 angle_atom2[m][num_angle[m]] = atom2;
1080 angle_atom3[m][num_angle[m]] = atom3;
1081 num_angle[m]++;
1082 m = atom3-1;
1083 angle_type[m][num_angle[m]] = itype;
1084 angle_atom1[m][num_angle[m]] = atom1;
1085 angle_atom2[m][num_angle[m]] = atom2;
1086 angle_atom3[m][num_angle[m]] = atom3;
1087 num_angle[m]++;
1088 }
1089 } else {
1090 count[atom2-1]++;
1091 if (newton_bond == 0) {
1092 count[atom1-1]++;
1093 count[atom3-1]++;
1094 }
1095 }
1096 }
1097
1098 // angle_per_atom = max of count vector
1099
1100 if (flag == 0) {
1101 angle_per_atom = 0;
1102 for (int i = 0; i < natoms; i++)
1103 angle_per_atom = MAX(angle_per_atom,count[i]);
1104 }
1105 }
1106
1107 /* ----------------------------------------------------------------------
1108 read dihedrals from file
1109 store each with all 4 atoms if newton_bond = 0
1110 if flag = 0, just count dihedrals/atom
1111 if flag = 1, store them with atoms
1112 ------------------------------------------------------------------------- */
1113
dihedrals(int flag,char * line)1114 void Molecule::dihedrals(int flag, char *line)
1115 {
1116 int itype;
1117 tagint m,atom1,atom2,atom3,atom4;
1118 int newton_bond = force->newton_bond;
1119
1120 if (flag == 0)
1121 for (int i = 0; i < natoms; i++) count[i] = 0;
1122 else
1123 for (int i = 0; i < natoms; i++) num_dihedral[i] = 0;
1124
1125 for (int i = 0; i < ndihedrals; i++) {
1126 readline(line);
1127
1128 try {
1129 ValueTokenizer values(utils::trim_comment(line));
1130 if (values.count() != 6)
1131 error->all(FLERR,"Invalid line in Dihedrals section of "
1132 "molecule file: {}",line);
1133
1134 values.next_int();
1135 itype = values.next_int();
1136 atom1 = values.next_tagint();
1137 atom2 = values.next_tagint();
1138 atom3 = values.next_tagint();
1139 atom4 = values.next_tagint();
1140 } catch (TokenizerException &e) {
1141 error->all(FLERR, "Invalid line in Dihedrals section of "
1142 "molecule file: {}\n{}",e.what(),line);
1143 }
1144
1145 itype += doffset;
1146
1147 if ((atom1 <= 0) || (atom1 > natoms) ||
1148 (atom2 <= 0) || (atom2 > natoms) ||
1149 (atom3 <= 0) || (atom3 > natoms) ||
1150 (atom4 <= 0) || (atom4 > natoms) ||
1151 (atom1 == atom2) || (atom1 == atom3) || (atom1 == atom4) ||
1152 (atom2 == atom3) || (atom2 == atom4) || (atom3 == atom4))
1153 error->all(FLERR,
1154 "Invalid atom ID in dihedrals section of molecule file");
1155 if ((itype <= 0) || (domain->box_exist && (itype > atom->ndihedraltypes)))
1156 error->all(FLERR,"Invalid dihedral type in Dihedrals section of molecule file");
1157
1158 if (flag) {
1159 m = atom2-1;
1160 ndihedraltypes = MAX(ndihedraltypes,itype);
1161 dihedral_type[m][num_dihedral[m]] = itype;
1162 dihedral_atom1[m][num_dihedral[m]] = atom1;
1163 dihedral_atom2[m][num_dihedral[m]] = atom2;
1164 dihedral_atom3[m][num_dihedral[m]] = atom3;
1165 dihedral_atom4[m][num_dihedral[m]] = atom4;
1166 num_dihedral[m]++;
1167 if (newton_bond == 0) {
1168 m = atom1-1;
1169 dihedral_type[m][num_dihedral[m]] = itype;
1170 dihedral_atom1[m][num_dihedral[m]] = atom1;
1171 dihedral_atom2[m][num_dihedral[m]] = atom2;
1172 dihedral_atom3[m][num_dihedral[m]] = atom3;
1173 dihedral_atom4[m][num_dihedral[m]] = atom4;
1174 num_dihedral[m]++;
1175 m = atom3-1;
1176 dihedral_type[m][num_dihedral[m]] = itype;
1177 dihedral_atom1[m][num_dihedral[m]] = atom1;
1178 dihedral_atom2[m][num_dihedral[m]] = atom2;
1179 dihedral_atom3[m][num_dihedral[m]] = atom3;
1180 dihedral_atom4[m][num_dihedral[m]] = atom4;
1181 num_dihedral[m]++;
1182 m = atom4-1;
1183 dihedral_type[m][num_dihedral[m]] = itype;
1184 dihedral_atom1[m][num_dihedral[m]] = atom1;
1185 dihedral_atom2[m][num_dihedral[m]] = atom2;
1186 dihedral_atom3[m][num_dihedral[m]] = atom3;
1187 dihedral_atom4[m][num_dihedral[m]] = atom4;
1188 num_dihedral[m]++;
1189 }
1190 } else {
1191 count[atom2-1]++;
1192 if (newton_bond == 0) {
1193 count[atom1-1]++;
1194 count[atom3-1]++;
1195 count[atom4-1]++;
1196 }
1197 }
1198 }
1199
1200 // dihedral_per_atom = max of count vector
1201
1202 if (flag == 0) {
1203 dihedral_per_atom = 0;
1204 for (int i = 0; i < natoms; i++)
1205 dihedral_per_atom = MAX(dihedral_per_atom,count[i]);
1206 }
1207 }
1208
1209 /* ----------------------------------------------------------------------
1210 read impropers from file
1211 store each with all 4 atoms if newton_bond = 0
1212 if flag = 0, just count impropers/atom
1213 if flag = 1, store them with atoms
1214 ------------------------------------------------------------------------- */
1215
impropers(int flag,char * line)1216 void Molecule::impropers(int flag, char *line)
1217 {
1218 int itype;
1219 tagint m,atom1,atom2,atom3,atom4;
1220 int newton_bond = force->newton_bond;
1221
1222 if (flag == 0)
1223 for (int i = 0; i < natoms; i++) count[i] = 0;
1224 else
1225 for (int i = 0; i < natoms; i++) num_improper[i] = 0;
1226
1227 for (int i = 0; i < nimpropers; i++) {
1228 readline(line);
1229
1230 try {
1231 ValueTokenizer values(utils::trim_comment(line));
1232 if (values.count() != 6)
1233 error->all(FLERR,"Invalid line in Impropers section of "
1234 "molecule file: {}",line);
1235 values.next_int();
1236 itype = values.next_int();
1237 atom1 = values.next_tagint();
1238 atom2 = values.next_tagint();
1239 atom3 = values.next_tagint();
1240 atom4 = values.next_tagint();
1241 } catch (TokenizerException &e) {
1242 error->all(FLERR, "Invalid line in Impropers section of "
1243 "molecule file: {}\n{}",e.what(),line);
1244 }
1245
1246 itype += ioffset;
1247
1248 if ((atom1 <= 0) || (atom1 > natoms) ||
1249 (atom2 <= 0) || (atom2 > natoms) ||
1250 (atom3 <= 0) || (atom3 > natoms) ||
1251 (atom4 <= 0) || (atom4 > natoms) ||
1252 (atom1 == atom2) || (atom1 == atom3) || (atom1 == atom4) ||
1253 (atom2 == atom3) || (atom2 == atom4) || (atom3 == atom4))
1254 error->all(FLERR,
1255 "Invalid atom ID in impropers section of molecule file");
1256 if ((itype <= 0) || (domain->box_exist && (itype > atom->nimpropertypes)))
1257 error->all(FLERR,"Invalid improper type in Impropers section of molecule file");
1258
1259 if (flag) {
1260 m = atom2-1;
1261 nimpropertypes = MAX(nimpropertypes,itype);
1262 improper_type[m][num_improper[m]] = itype;
1263 improper_atom1[m][num_improper[m]] = atom1;
1264 improper_atom2[m][num_improper[m]] = atom2;
1265 improper_atom3[m][num_improper[m]] = atom3;
1266 improper_atom4[m][num_improper[m]] = atom4;
1267 num_improper[m]++;
1268 if (newton_bond == 0) {
1269 m = atom1-1;
1270 improper_type[m][num_improper[m]] = itype;
1271 improper_atom1[m][num_improper[m]] = atom1;
1272 improper_atom2[m][num_improper[m]] = atom2;
1273 improper_atom3[m][num_improper[m]] = atom3;
1274 improper_atom4[m][num_improper[m]] = atom4;
1275 num_improper[m]++;
1276 m = atom3-1;
1277 improper_type[m][num_improper[m]] = itype;
1278 improper_atom1[m][num_improper[m]] = atom1;
1279 improper_atom2[m][num_improper[m]] = atom2;
1280 improper_atom3[m][num_improper[m]] = atom3;
1281 improper_atom4[m][num_improper[m]] = atom4;
1282 num_improper[m]++;
1283 m = atom4-1;
1284 improper_type[m][num_improper[m]] = itype;
1285 improper_atom1[m][num_improper[m]] = atom1;
1286 improper_atom2[m][num_improper[m]] = atom2;
1287 improper_atom3[m][num_improper[m]] = atom3;
1288 improper_atom4[m][num_improper[m]] = atom4;
1289 num_improper[m]++;
1290 }
1291 } else {
1292 count[atom2-1]++;
1293 if (newton_bond == 0) {
1294 count[atom1-1]++;
1295 count[atom3-1]++;
1296 count[atom4-1]++;
1297 }
1298 }
1299 }
1300
1301 // improper_per_atom = max of count vector
1302
1303 if (flag == 0) {
1304 improper_per_atom = 0;
1305 for (int i = 0; i < natoms; i++)
1306 improper_per_atom = MAX(improper_per_atom,count[i]);
1307 }
1308 }
1309
1310 /* ----------------------------------------------------------------------
1311 read 3 special bonds counts from file
1312 if flag = 0, just tally maxspecial
1313 if flag = 1, store them with atoms
1314 ------------------------------------------------------------------------- */
1315
nspecial_read(int flag,char * line)1316 void Molecule::nspecial_read(int flag, char *line)
1317 {
1318 if (flag == 0) maxspecial = 0;
1319
1320 for (int i = 0; i < natoms; i++) {
1321 readline(line);
1322
1323 int c1, c2, c3;
1324
1325 try {
1326 ValueTokenizer values(utils::trim_comment(line));
1327 if (values.count() != 4)
1328 error->all(FLERR,"Invalid line in Special Bond Counts section of "
1329 "molecule file: {}",line);
1330 values.next_int();
1331 c1 = values.next_tagint();
1332 c2 = values.next_tagint();
1333 c3 = values.next_tagint();
1334 } catch (TokenizerException &e) {
1335 error->all(FLERR, "Invalid line in Special Bond Counts section of "
1336 "molecule file: {}\n{}",e.what(),line);
1337 }
1338
1339 if (flag) {
1340 nspecial[i][0] = c1;
1341 nspecial[i][1] = c1+c2;
1342 nspecial[i][2] = c1+c2+c3;
1343 } else maxspecial = MAX(maxspecial,c1+c2+c3);
1344 }
1345 }
1346
1347 /* ----------------------------------------------------------------------
1348 read special bond indices from file
1349 ------------------------------------------------------------------------- */
1350
special_read(char * line)1351 void Molecule::special_read(char *line)
1352 {
1353 try {
1354 for (int i = 0; i < natoms; i++) {
1355 readline(line);
1356
1357 ValueTokenizer values(utils::trim_comment(line));
1358 int nwords = values.count();
1359
1360 if (nwords != nspecial[i][2]+1)
1361 error->all(FLERR,"Molecule file special list "
1362 "does not match special count");
1363
1364 values.next_int(); // ignore
1365
1366 for (int m = 1; m < nwords; m++) {
1367 special[i][m-1] = values.next_tagint();
1368 if (special[i][m-1] <= 0 || special[i][m-1] > natoms ||
1369 special[i][m-1] == i+1)
1370 error->all(FLERR,"Invalid atom index in Special Bonds "
1371 "section of molecule file");
1372 }
1373 }
1374 } catch (TokenizerException &e) {
1375 error->all(FLERR, "Invalid line in Special Bonds section of "
1376 "molecule file: {}\n{}",e.what(),line);
1377 }
1378 }
1379
1380 /* ----------------------------------------------------------------------
1381 auto generate special bond info
1382 ------------------------------------------------------------------------- */
1383
special_generate()1384 void Molecule::special_generate()
1385 {
1386 int newton_bond = force->newton_bond;
1387 tagint atom1,atom2;
1388
1389 // temporary array for special atoms
1390
1391 tagint **tmpspecial;
1392 memory->create(tmpspecial,natoms,atom->maxspecial,"molecule:tmpspecial");
1393 memset(&tmpspecial[0][0],0,sizeof(tagint)*natoms*atom->maxspecial);
1394
1395 for (int i = 0; i < natoms; i++) count[i] = 0;
1396
1397 // 1-2 neighbors
1398
1399 if (newton_bond) {
1400 for (int i = 0; i < natoms; i++) {
1401 for (int j = 0; j < num_bond[i]; j++) {
1402 atom1 = i;
1403 atom2 = bond_atom[i][j]-1;
1404 nspecial[i][0]++;
1405 nspecial[atom2][0]++;
1406 if (count[i] >= atom->maxspecial || count[atom2] >= atom->maxspecial)
1407 error->all(FLERR,"Molecule auto special bond generation overflow");
1408 tmpspecial[i][count[i]++] = atom2 + 1;
1409 tmpspecial[atom2][count[atom2]++] = i + 1;
1410 }
1411 }
1412 } else {
1413 for (int i = 0; i < natoms; i++) {
1414 nspecial[i][0] = num_bond[i];
1415 for (int j = 0; j < num_bond[i]; j++) {
1416 atom1 = i;
1417 atom2 = bond_atom[i][j];
1418 if (count[atom1] >= atom->maxspecial)
1419 error->all(FLERR,"Molecule auto special bond generation overflow");
1420 tmpspecial[i][count[atom1]++] = atom2;
1421 }
1422 }
1423 }
1424
1425 // 1-3 neighbors with no duplicates
1426
1427 for (int i = 0; i < natoms; i++) nspecial[i][1] = nspecial[i][0];
1428
1429 int dedup;
1430 for (int i = 0; i < natoms; i++) {
1431 for (int m = 0; m < nspecial[i][0]; m++) {
1432 for (int j = 0; j < nspecial[tmpspecial[i][m]-1][0]; j++) {
1433 dedup = 0;
1434 for (int k =0; k < count[i]; k++) {
1435 if (tmpspecial[tmpspecial[i][m]-1][j] == tmpspecial[i][k] ||
1436 tmpspecial[tmpspecial[i][m]-1][j] == i+1) {
1437 dedup = 1;
1438 }
1439 }
1440 if (!dedup) {
1441 if (count[i] >= atom->maxspecial)
1442 error->all(FLERR,"Molecule auto special bond generation overflow");
1443 tmpspecial[i][count[i]++] = tmpspecial[tmpspecial[i][m]-1][j];
1444 nspecial[i][1]++;
1445 }
1446 }
1447 }
1448 }
1449
1450 // 1-4 neighbors with no duplicates
1451
1452 for (int i = 0; i < natoms; i++) nspecial[i][2] = nspecial[i][1];
1453
1454 for (int i = 0; i < natoms; i++) {
1455 for (int m = nspecial[i][0]; m < nspecial[i][1]; m++) {
1456 for (int j = 0; j < nspecial[tmpspecial[i][m]-1][0]; j++) {
1457 dedup = 0;
1458 for (int k =0; k < count[i]; k++) {
1459 if (tmpspecial[tmpspecial[i][m]-1][j] == tmpspecial[i][k] ||
1460 tmpspecial[tmpspecial[i][m]-1][j] == i+1) {
1461 dedup = 1;
1462 }
1463 }
1464 if (!dedup) {
1465 if (count[i] >= atom->maxspecial)
1466 error->all(FLERR,"Molecule auto special bond generation overflow");
1467 tmpspecial[i][count[i]++] = tmpspecial[tmpspecial[i][m]-1][j];
1468 nspecial[i][2]++;
1469 }
1470 }
1471 }
1472 }
1473
1474 maxspecial = 0;
1475 for (int i = 0; i < natoms; i++)
1476 maxspecial = MAX(maxspecial,nspecial[i][2]);
1477
1478 memory->create(special,natoms,maxspecial,"molecule:special");
1479 for (int i = 0; i < natoms; i++)
1480 for (int j = 0; j < nspecial[i][2]; j++)
1481 special[i][j] = tmpspecial[i][j];
1482
1483 memory->destroy(tmpspecial);
1484 }
1485
1486 /* ----------------------------------------------------------------------
1487 read SHAKE flags from file
1488 ------------------------------------------------------------------------- */
1489
shakeflag_read(char * line)1490 void Molecule::shakeflag_read(char *line)
1491 {
1492 try {
1493 for (int i = 0; i < natoms; i++) {
1494 readline(line);
1495
1496 ValueTokenizer values(utils::trim_comment(line));
1497
1498 if (values.count() != 2)
1499 error->all(FLERR,"Invalid Shake Flags section in molecule file");
1500
1501 values.next_int();
1502 shake_flag[i] = values.next_int();
1503 }
1504 } catch (TokenizerException &e) {
1505 error->all(FLERR, "Invalid Shake Flags section in molecule file\n"
1506 "{}", e.what());
1507 }
1508
1509 for (int i = 0; i < natoms; i++)
1510 if (shake_flag[i] < 0 || shake_flag[i] > 4)
1511 error->all(FLERR,"Invalid shake flag in molecule file");
1512 }
1513
1514 /* ----------------------------------------------------------------------
1515 read SHAKE atom info from file
1516 ------------------------------------------------------------------------- */
1517
shakeatom_read(char * line)1518 void Molecule::shakeatom_read(char *line)
1519 {
1520 int nmatch=0, nwant=0;
1521 try {
1522 for (int i = 0; i < natoms; i++) {
1523 readline(line);
1524
1525 ValueTokenizer values(utils::trim_comment(line));
1526 nmatch = values.count();
1527
1528 switch (shake_flag[i]) {
1529 case 1:
1530 values.next_int();
1531 shake_atom[i][0] = values.next_tagint();
1532 shake_atom[i][1] = values.next_tagint();
1533 shake_atom[i][2] = values.next_tagint();
1534 nwant = 4;
1535 break;
1536
1537 case 2:
1538 values.next_int();
1539 shake_atom[i][0] = values.next_tagint();
1540 shake_atom[i][1] = values.next_tagint();
1541 nwant = 3;
1542 break;
1543
1544 case 3:
1545 values.next_int();
1546 shake_atom[i][0] = values.next_tagint();
1547 shake_atom[i][1] = values.next_tagint();
1548 shake_atom[i][2] = values.next_tagint();
1549 nwant = 4;
1550 break;
1551
1552 case 4:
1553 values.next_int();
1554 shake_atom[i][0] = values.next_tagint();
1555 shake_atom[i][1] = values.next_tagint();
1556 shake_atom[i][2] = values.next_tagint();
1557 shake_atom[i][3] = values.next_tagint();
1558 nwant = 5;
1559 break;
1560
1561 case 0:
1562 values.next_int();
1563 nwant = 1;
1564 break;
1565
1566 default:
1567 error->all(FLERR,"Invalid shake atom in molecule file");
1568 }
1569
1570 if (nmatch != nwant)
1571 error->all(FLERR,"Invalid shake atom in molecule file");
1572 }
1573
1574 } catch (TokenizerException &e) {
1575 error->all(FLERR,"Invalid shake atom in molecule file\n"
1576 "{}", e.what());
1577 }
1578
1579 for (int i = 0; i < natoms; i++) {
1580 int m = shake_flag[i];
1581 if (m == 1) m = 3;
1582 for (int j = 0; j < m; j++)
1583 if (shake_atom[i][j] <= 0 || shake_atom[i][j] > natoms)
1584 error->all(FLERR,"Invalid shake atom in molecule file");
1585 }
1586 }
1587
1588 /* ----------------------------------------------------------------------
1589 read SHAKE bond type info from file
1590 ------------------------------------------------------------------------- */
1591
shaketype_read(char * line)1592 void Molecule::shaketype_read(char *line)
1593 {
1594 try {
1595 int nmatch=0, nwant=0;
1596 for (int i = 0; i < natoms; i++) {
1597 readline(line);
1598
1599 ValueTokenizer values(utils::trim_comment(line));
1600 nmatch = values.count();
1601
1602 switch (shake_flag[i]) {
1603 case 1:
1604 values.next_int();
1605 shake_type[i][0] = values.next_int();
1606 shake_type[i][1] = values.next_int();
1607 shake_type[i][2] = values.next_int();
1608 nwant = 4;
1609 break;
1610
1611 case 2:
1612 values.next_int();
1613 shake_type[i][0] = values.next_int();
1614 nwant = 2;
1615 break;
1616
1617 case 3:
1618 values.next_int();
1619 shake_type[i][0] = values.next_int();
1620 shake_type[i][1] = values.next_int();
1621 nwant = 3;
1622 break;
1623
1624 case 4:
1625 values.next_int();
1626 shake_type[i][0] = values.next_int();
1627 shake_type[i][1] = values.next_int();
1628 shake_type[i][2] = values.next_int();
1629 nwant = 4;
1630 break;
1631
1632 case 0:
1633 values.next_int();
1634 nwant = 1;
1635 break;
1636
1637 default:
1638 error->all(FLERR,"Invalid shake type data in molecule file");
1639 }
1640
1641 if (nmatch != nwant)
1642 error->all(FLERR,"Invalid shake type data in molecule file");
1643 }
1644 } catch (TokenizerException &e) {
1645 error->all(FLERR, "Invalid shake type data in molecule file\n",
1646 "{}", e.what());
1647 }
1648
1649 for (int i = 0; i < natoms; i++) {
1650 int m = shake_flag[i];
1651 if (m == 1) m = 3;
1652 for (int j = 0; j < m-1; j++)
1653 if (shake_type[i][j] <= 0)
1654 error->all(FLERR,"Invalid shake bond type in molecule file");
1655 if (shake_flag[i] == 1)
1656 if (shake_type[i][2] <= 0)
1657 error->all(FLERR,"Invalid shake angle type in molecule file");
1658 }
1659 }
1660
1661 /* ----------------------------------------------------------------------
1662 read body params from file
1663 pflag = 0/1 for integer/double params
1664 ------------------------------------------------------------------------- */
1665
body(int flag,int pflag,char * line)1666 void Molecule::body(int flag, int pflag, char *line)
1667 {
1668 int nparam = nibody;
1669 if (pflag) nparam = ndbody;
1670
1671 int nword = 0;
1672
1673 try {
1674 while (nword < nparam) {
1675 readline(line);
1676
1677 ValueTokenizer values(utils::trim_comment(line));
1678 int ncount = values.count();
1679
1680 if (ncount == 0)
1681 error->all(FLERR,"Too few values in body section of molecule file");
1682 if (nword+ncount > nparam)
1683 error->all(FLERR,"Too many values in body section of molecule file");
1684
1685 if (flag) {
1686 if (pflag == 0) {
1687 while (values.has_next()) {
1688 ibodyparams[nword++] = values.next_int();
1689 }
1690 } else {
1691 while (values.has_next()) {
1692 dbodyparams[nword++] = values.next_double();
1693 }
1694 }
1695 } else nword += ncount;
1696 }
1697 } catch (TokenizerException &e) {
1698 error->all(FLERR, "Invalid body params in molecule file\n",
1699 "{}", e.what());
1700 }
1701 }
1702
1703 /* ----------------------------------------------------------------------
1704 return fragment index if name matches existing fragment, -1 if no such fragment
1705 ------------------------------------------------------------------------- */
1706
findfragment(const char * name)1707 int Molecule::findfragment(const char *name)
1708 {
1709 for (int i = 0; i < nfragments; i++)
1710 if (fragmentnames[i] == name) return i;
1711 return -1;
1712 }
1713
1714 /* ----------------------------------------------------------------------
1715 error check molecule attributes and topology against system settings
1716 flag = 0, just check this molecule
1717 flag = 1, check all molecules in set, this is 1st molecule in set
1718 ------------------------------------------------------------------------- */
1719
check_attributes(int flag)1720 void Molecule::check_attributes(int flag)
1721 {
1722 int n = 1;
1723 if (flag) n = nset;
1724 int imol = atom->find_molecule(id);
1725
1726 for (int i = imol; i < imol+n; i++) {
1727 Molecule *onemol = atom->molecules[imol];
1728
1729 // check per-atom attributes of molecule
1730 // warn if not a match
1731
1732 int mismatch = 0;
1733 if (onemol->qflag && !atom->q_flag) mismatch = 1;
1734 if (onemol->radiusflag && !atom->radius_flag) mismatch = 1;
1735 if (onemol->rmassflag && !atom->rmass_flag) mismatch = 1;
1736
1737 if (mismatch && me == 0)
1738 error->warning(FLERR,
1739 "Molecule attributes do not match system attributes");
1740
1741 // for all atom styles, check nbondtype,etc
1742
1743 mismatch = 0;
1744 if (atom->nbondtypes < onemol->nbondtypes) mismatch = 1;
1745 if (atom->nangletypes < onemol->nangletypes) mismatch = 1;
1746 if (atom->ndihedraltypes < onemol->ndihedraltypes) mismatch = 1;
1747 if (atom->nimpropertypes < onemol->nimpropertypes) mismatch = 1;
1748
1749 if (mismatch)
1750 error->all(FLERR,"Molecule topology type exceeds system topology type");
1751
1752 // for molecular atom styles, check bond_per_atom,etc + maxspecial
1753 // do not check for atom style template, since nothing stored per atom
1754
1755 if (atom->molecular == Atom::MOLECULAR) {
1756 if (atom->avec->bonds_allow &&
1757 atom->bond_per_atom < onemol->bond_per_atom) mismatch = 1;
1758 if (atom->avec->angles_allow &&
1759 atom->angle_per_atom < onemol->angle_per_atom) mismatch = 1;
1760 if (atom->avec->dihedrals_allow &&
1761 atom->dihedral_per_atom < onemol->dihedral_per_atom) mismatch = 1;
1762 if (atom->avec->impropers_allow &&
1763 atom->improper_per_atom < onemol->improper_per_atom) mismatch = 1;
1764 if (atom->maxspecial < onemol->maxspecial) mismatch = 1;
1765
1766 if (mismatch)
1767 error->all(FLERR,"Molecule topology/atom exceeds system topology/atom");
1768 }
1769
1770 // warn if molecule topology defined but no special settings
1771
1772 if (onemol->bondflag && !onemol->specialflag)
1773 if (me == 0) error->warning(FLERR,"Molecule has bond topology "
1774 "but no special bond settings");
1775 }
1776 }
1777
1778 /* ----------------------------------------------------------------------
1779 init all data structures to empty
1780 ------------------------------------------------------------------------- */
1781
initialize()1782 void Molecule::initialize()
1783 {
1784 natoms = 0;
1785 nbonds = nangles = ndihedrals = nimpropers = 0;
1786 ntypes = 0;
1787 nmolecules = 1;
1788 nbondtypes = nangletypes = ndihedraltypes = nimpropertypes = 0;
1789 nibody = ndbody = 0;
1790 nfragments = 0;
1791
1792 bond_per_atom = angle_per_atom = dihedral_per_atom = improper_per_atom = 0;
1793 maxspecial = 0;
1794
1795 xflag = typeflag = moleculeflag = fragmentflag = qflag = radiusflag = rmassflag = 0;
1796 bondflag = angleflag = dihedralflag = improperflag = 0;
1797 nspecialflag = specialflag = 0;
1798 shakeflag = shakeflagflag = shakeatomflag = shaketypeflag = 0;
1799 bodyflag = ibodyflag = dbodyflag = 0;
1800
1801 centerflag = massflag = comflag = inertiaflag = 0;
1802 tag_require = 0;
1803
1804 x = nullptr;
1805 type = nullptr;
1806 q = nullptr;
1807 radius = nullptr;
1808 rmass = nullptr;
1809
1810 num_bond = nullptr;
1811 bond_type = nullptr;
1812 bond_atom = nullptr;
1813
1814 num_angle = nullptr;
1815 angle_type = nullptr;
1816 angle_atom1 = angle_atom2 = angle_atom3 = nullptr;
1817
1818 num_dihedral = nullptr;
1819 dihedral_type = nullptr;
1820 dihedral_atom1 = dihedral_atom2 = dihedral_atom3 = dihedral_atom4 = nullptr;
1821
1822 num_improper = nullptr;
1823 improper_type = nullptr;
1824 improper_atom1 = improper_atom2 = improper_atom3 = improper_atom4 = nullptr;
1825
1826 nspecial = nullptr;
1827 special = nullptr;
1828
1829 shake_flag = nullptr;
1830 shake_atom = nullptr;
1831 shake_type = nullptr;
1832
1833 ibodyparams = nullptr;
1834 dbodyparams = nullptr;
1835
1836 dx = nullptr;
1837 dxcom = nullptr;
1838 dxbody = nullptr;
1839 }
1840
1841 /* ----------------------------------------------------------------------
1842 allocate all data structures
1843 also initialize values for data structures that are always allocated
1844 ------------------------------------------------------------------------- */
1845
allocate()1846 void Molecule::allocate()
1847 {
1848 if (xflag) memory->create(x,natoms,3,"molecule:x");
1849 if (typeflag) memory->create(type,natoms,"molecule:type");
1850 if (moleculeflag) memory->create(molecule,natoms,"molecule:molecule");
1851 if (fragmentflag) {
1852 fragmentnames.resize(nfragments);
1853 memory->create(fragmentmask,nfragments,natoms,"molecule:fragmentmask");
1854 for (int i = 0; i < nfragments; i++)
1855 for (int j = 0; j < natoms; j++) fragmentmask[i][j] = 0;
1856 }
1857 if (qflag) memory->create(q,natoms,"molecule:q");
1858 if (radiusflag) memory->create(radius,natoms,"molecule:radius");
1859 if (rmassflag) memory->create(rmass,natoms,"molecule:rmass");
1860
1861 // always allocate num_bond,num_angle,etc and nspecial
1862 // even if not in molecule file, initialize to 0
1863 // this is so methods that use these arrays don't have to check they exist
1864
1865 memory->create(num_bond,natoms,"molecule:num_bond");
1866 for (int i = 0; i < natoms; i++) num_bond[i] = 0;
1867 memory->create(num_angle,natoms,"molecule:num_angle");
1868 for (int i = 0; i < natoms; i++) num_angle[i] = 0;
1869 memory->create(num_dihedral,natoms,"molecule:num_dihedral");
1870 for (int i = 0; i < natoms; i++) num_dihedral[i] = 0;
1871 memory->create(num_improper,natoms,"molecule:num_improper");
1872 for (int i = 0; i < natoms; i++) num_improper[i] = 0;
1873 memory->create(nspecial,natoms,3,"molecule:nspecial");
1874 for (int i = 0; i < natoms; i++)
1875 nspecial[i][0] = nspecial[i][1] = nspecial[i][2] = 0;
1876
1877 if (specialflag)
1878 memory->create(special,natoms,maxspecial,"molecule:special");
1879
1880 if (bondflag) {
1881 memory->create(bond_type,natoms,bond_per_atom,
1882 "molecule:bond_type");
1883 memory->create(bond_atom,natoms,bond_per_atom,
1884 "molecule:bond_atom");
1885 }
1886
1887 if (angleflag) {
1888 memory->create(angle_type,natoms,angle_per_atom,
1889 "molecule:angle_type");
1890 memory->create(angle_atom1,natoms,angle_per_atom,
1891 "molecule:angle_atom1");
1892 memory->create(angle_atom2,natoms,angle_per_atom,
1893 "molecule:angle_atom2");
1894 memory->create(angle_atom3,natoms,angle_per_atom,
1895 "molecule:angle_atom3");
1896 }
1897
1898 if (dihedralflag) {
1899 memory->create(dihedral_type,natoms,dihedral_per_atom,
1900 "molecule:dihedral_type");
1901 memory->create(dihedral_atom1,natoms,dihedral_per_atom,
1902 "molecule:dihedral_atom1");
1903 memory->create(dihedral_atom2,natoms,dihedral_per_atom,
1904 "molecule:dihedral_atom2");
1905 memory->create(dihedral_atom3,natoms,dihedral_per_atom,
1906 "molecule:dihedral_atom3");
1907 memory->create(dihedral_atom4,natoms,dihedral_per_atom,
1908 "molecule:dihedral_atom4");
1909 }
1910
1911 if (improperflag) {
1912 memory->create(improper_type,natoms,improper_per_atom,
1913 "molecule:improper_type");
1914 memory->create(improper_atom1,natoms,improper_per_atom,
1915 "molecule:improper_atom1");
1916 memory->create(improper_atom2,natoms,improper_per_atom,
1917 "molecule:improper_atom2");
1918 memory->create(improper_atom3,natoms,improper_per_atom,
1919 "molecule:improper_atom3");
1920 memory->create(improper_atom4,natoms,improper_per_atom,
1921 "molecule:improper_atom4");
1922 }
1923
1924 if (shakeflag) {
1925 memory->create(shake_flag,natoms,"molecule:shake_flag");
1926 memory->create(shake_atom,natoms,4,"molecule:shake_flag");
1927 memory->create(shake_type,natoms,3,"molecule:shake_flag");
1928 }
1929
1930 if (bodyflag) {
1931 if (nibody) memory->create(ibodyparams,nibody,"molecule:ibodyparams");
1932 if (ndbody) memory->create(dbodyparams,ndbody,"molecule:dbodyparams");
1933 }
1934 }
1935
1936 /* ----------------------------------------------------------------------
1937 deallocate all data structures
1938 ------------------------------------------------------------------------- */
1939
deallocate()1940 void Molecule::deallocate()
1941 {
1942 memory->destroy(x);
1943 memory->destroy(type);
1944 memory->destroy(molecule);
1945 memory->destroy(q);
1946 memory->destroy(radius);
1947 memory->destroy(rmass);
1948
1949 memory->destroy(molecule);
1950 memory->destroy(fragmentmask);
1951
1952 if (fragmentflag) {
1953 fragmentnames.clear();
1954 }
1955
1956 memory->destroy(num_bond);
1957 memory->destroy(bond_type);
1958 memory->destroy(bond_atom);
1959
1960 memory->destroy(num_angle);
1961 memory->destroy(angle_type);
1962 memory->destroy(angle_atom1);
1963 memory->destroy(angle_atom2);
1964 memory->destroy(angle_atom3);
1965
1966 memory->destroy(num_dihedral);
1967 memory->destroy(dihedral_type);
1968 memory->destroy(dihedral_atom1);
1969 memory->destroy(dihedral_atom2);
1970 memory->destroy(dihedral_atom3);
1971 memory->destroy(dihedral_atom4);
1972
1973 memory->destroy(num_improper);
1974 memory->destroy(improper_type);
1975 memory->destroy(improper_atom1);
1976 memory->destroy(improper_atom2);
1977 memory->destroy(improper_atom3);
1978 memory->destroy(improper_atom4);
1979
1980 memory->destroy(nspecial);
1981 memory->destroy(special);
1982
1983 memory->destroy(shake_flag);
1984 memory->destroy(shake_atom);
1985 memory->destroy(shake_type);
1986
1987 memory->destroy(dx);
1988 memory->destroy(dxcom);
1989 memory->destroy(dxbody);
1990
1991 memory->destroy(ibodyparams);
1992 memory->destroy(dbodyparams);
1993 }
1994
1995 /* ----------------------------------------------------------------------
1996 read and bcast a line
1997 ------------------------------------------------------------------------- */
1998
readline(char * line)1999 void Molecule::readline(char *line)
2000 {
2001 int n;
2002 if (me == 0) {
2003 if (fgets(line,MAXLINE,fp) == nullptr) n = 0;
2004 else n = strlen(line) + 1;
2005 }
2006 MPI_Bcast(&n,1,MPI_INT,0,world);
2007 if (n == 0) error->all(FLERR,"Unexpected end of molecule file");
2008 MPI_Bcast(line,n,MPI_CHAR,0,world);
2009 }
2010
2011 /* ----------------------------------------------------------------------
2012 extract keyword from line
2013 flag = 0, read and bcast line
2014 flag = 1, line has already been read
2015 ------------------------------------------------------------------------- */
2016
parse_keyword(int flag,char * line)2017 std::string Molecule::parse_keyword(int flag, char *line)
2018 {
2019 char line2[MAXLINE];
2020 if (flag) {
2021
2022 // read upto non-blank line plus 1 following line
2023 // eof is set to 1 if any read hits end-of-file
2024
2025 int eof = 0;
2026 if (me == 0) {
2027 if (fgets(line,MAXLINE,fp) == nullptr) eof = 1;
2028 while (eof == 0 && strspn(line," \t\n\r") == strlen(line)) {
2029 if (fgets(line,MAXLINE,fp) == nullptr) eof = 1;
2030 }
2031 if (fgets(line2,MAXLINE,fp) == nullptr) eof = 1;
2032 }
2033
2034 // if eof, set keyword empty and return
2035
2036 MPI_Bcast(&eof,1,MPI_INT,0,world);
2037 if (eof) {
2038 return std::string("");
2039 }
2040
2041 // bcast keyword line to all procs
2042
2043 MPI_Bcast(line,MAXLINE,MPI_CHAR,0,world);
2044 }
2045
2046 // return non-whitespace and non-comment portion of line
2047
2048 return utils::trim(utils::trim_comment(line));
2049 }
2050
2051 /* ----------------------------------------------------------------------
2052 skip N lines of file. Check if non-numeric content (e.g. keyword).
2053 ------------------------------------------------------------------------- */
2054
skip_lines(int n,char * line,const std::string & section)2055 void Molecule::skip_lines(int n, char *line, const std::string §ion)
2056 {
2057 for (int i = 0; i < n; i++) {
2058 readline(line);
2059 if (utils::strmatch(utils::trim(utils::trim_comment(line)),"^[A-Za-z ]+$"))
2060 error->one(FLERR,"Unexpected line in molecule file while "
2061 "skipping {} section:\n{}",section,line);
2062 }
2063 }
2064
2065 /* ----------------------------------------------------------------------
2066 proc 0 prints molecule params
2067 ------------------------------------------------------------------------- */
2068
2069 /*
2070
2071 void Molecule::print()
2072 {
2073 printf("MOLECULE %s\n",id);
2074 printf(" %d natoms\n",natoms);
2075 if (nbonds) printf(" %d nbonds\n",nbonds);
2076 if (nangles) printf(" %d nangles\n",nangles);
2077 if (ndihedrals) printf(" %d ndihedrals\n",ndihedrals);
2078 if (nimpropers) printf(" %d nimpropers\n",nimpropers);
2079
2080 if (xflag) {
2081 printf( "Coords:\n");
2082 for (int i = 0; i < natoms; i++)
2083 printf(" %d %g %g %g\n",i+1,x[i][0],x[i][1],x[i][2]);
2084 }
2085 if (typeflag) {
2086 printf( "Types:\n");
2087 for (int i = 0; i < natoms; i++)
2088 printf(" %d %d\n",i+1,type[i]);
2089 }
2090 if (qflag) {
2091 printf( "Charges:\n");
2092 for (int i = 0; i < natoms; i++)
2093 printf(" %d %g\n",i+1,q[i]);
2094 }
2095 if (radiusflag) {
2096 printf( "Radii:\n");
2097 for (int i = 0; i < natoms; i++)
2098 printf(" %d %g\n",i+1,radius[i]);
2099 }
2100 if (rmassflag) {
2101 printf( "Masses:\n");
2102 for (int i = 0; i < natoms; i++)
2103 printf(" %d %g\n",i+1,rmass[i]);
2104 }
2105
2106 if (bondflag) {
2107 printf( "Bonds:\n");
2108 for (int i = 0; i < natoms; i++) {
2109 printf(" %d %d\n",i+1,num_bond[i]);
2110 for (int j = 0; j < num_bond[i]; j++)
2111 printf(" %d %d %d %d\n",j+1,bond_type[i][j],i+1,bond_atom[i][j]);
2112 }
2113 }
2114 if (angleflag) {
2115 printf( "Angles:\n");
2116 for (int i = 0; i < natoms; i++) {
2117 printf(" %d %d\n",i+1,num_angle[i]);
2118 for (int j = 0; j < num_angle[i]; j++)
2119 printf(" %d %d %d %d %d\n",
2120 j+1,angle_type[i][j],
2121 angle_atom1[i][j],angle_atom2[i][j],angle_atom3[i][j]);
2122 }
2123 }
2124 if (dihedralflag) {
2125 printf( "Dihedrals:\n");
2126 for (int i = 0; i < natoms; i++) {
2127 printf(" %d %d\n",i+1,num_dihedral[i]);
2128 for (int j = 0; j < num_dihedral[i]; j++)
2129 printf(" %d %d %d %d %d %d\n",
2130 j+1,dihedral_type[i][j],
2131 dihedral_atom1[i][j],dihedral_atom2[i][j],
2132 dihedral_atom3[i][j],dihedral_atom4[i][j]);
2133 }
2134 }
2135 if (improperflag) {
2136 printf( "Impropers:\n");
2137 for (int i = 0; i < natoms; i++) {
2138 printf(" %d %d\n",i+1,num_improper[i]);
2139 for (int j = 0; j < num_improper[i]; j++)
2140 printf(" %d %d %d %d %d %d\n",
2141 j+1,improper_type[i][j],
2142 improper_atom1[i][j],improper_atom2[i][j],
2143 improper_atom3[i][j],improper_atom4[i][j]);
2144 }
2145 }
2146
2147 if (specialflag) {
2148 printf( "Special neighs:\n");
2149 for (int i = 0; i < natoms; i++) {
2150 printf(" %d %d %d %d\n",i+1,
2151 nspecial[i][0],nspecial[i][1]-nspecial[i][0],
2152 nspecial[i][2]-nspecial[i][1]);
2153 printf(" ");
2154 for (int j = 0; j < nspecial[i][2]; j++)
2155 printf(" %d",special[i][j]);
2156 printf("\n");
2157 }
2158 }
2159 }
2160
2161 */
2162