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
36 LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
37 http://lammps.sandia.gov, Sandia National Laboratories
38 Steve Plimpton, sjplimp@sandia.gov
39
40 Copyright (2003) Sandia Corporation. Under the terms of Contract
41 DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
42 certain rights in this software. This software is distributed under
43 the GNU General Public License.
44 ------------------------------------------------------------------------- */
45
46 /* ----------------------------------------------------------------------
47 Contributing author: Timothy Sirk (ARL)
48 ------------------------------------------------------------------------- */
49
50 #include "lmptype.h"
51 #include <mpi.h>
52 #include <string.h>
53 #include <stdlib.h>
54 #include "dirent.h"
55 #include "read_dump.h"
56 #include "reader.h"
57 #include "style_reader.h"
58 #include "atom.h"
59 #include "atom_vec.h"
60 #include "update.h"
61 #include "modify.h"
62 #include "fix.h"
63 #include "domain.h"
64 #include "comm.h"
65 #include "irregular.h"
66 #include "error.h"
67 #include "memory.h"
68
69 using namespace LAMMPS_NS;
70
71 #define CHUNK 1024
72 #define EPSILON 1.0e-6
73
74 // also in reader_native.cpp
75
76 enum{ID,TYPE,X,Y,Z,VX,VY,VZ,OMEGAX,OMEGAY,OMEGAZ,Q,IX,IY,IZ,RADIUS,MASS,DENSITY,FX,FY,FZ};
77 enum{UNSET,NOSCALE_NOWRAP,NOSCALE_WRAP,SCALE_NOWRAP,SCALE_WRAP};
78
79 /* ---------------------------------------------------------------------- */
80
ReadDump(LAMMPS * lmp)81 ReadDump::ReadDump(LAMMPS *lmp) : Pointers(lmp)
82 {
83 MPI_Comm_rank(world,&me);
84 MPI_Comm_size(world,&nprocs);
85
86 dimension = domain->dimension;
87 triclinic = domain->triclinic;
88
89 nfile = 0;
90 files = NULL;
91
92 nfield = 0;
93 fieldtype = NULL;
94 fieldlabel = NULL;
95 fields = NULL;
96
97 int n = strlen("native") + 1;
98 readerstyle = new char[n];
99 strcpy(readerstyle,"native");
100
101 reader = NULL;
102 fp = NULL;
103 }
104
105 /* ---------------------------------------------------------------------- */
106
~ReadDump()107 ReadDump::~ReadDump()
108 {
109 for (int i = 0; i < nfile; i++) delete [] files[i];
110 delete [] files;
111 for (int i = 0; i < nfield; i++) delete [] fieldlabel[i];
112 delete [] fieldlabel;
113 delete [] fieldtype;
114 delete [] readerstyle;
115
116 memory->destroy(fields);
117 delete reader;
118 }
119
120 /* ---------------------------------------------------------------------- */
121
command(int narg,char ** arg)122 void ReadDump::command(int narg, char **arg)
123 {
124 if (domain->box_exist == 0)
125 error->all(FLERR,"Read_dump command before simulation box is defined");
126
127 if (narg < 2) error->all(FLERR,"Illegal read_dump command");
128
129 store_files(1,&arg[0]);
130 bigint nstep = ATOBIGINT(arg[1]);
131
132 int nremain = narg - 2;
133 if (nremain) nremain = fields_and_keywords(nremain,&arg[narg-nremain]);
134 else nremain = fields_and_keywords(0,NULL);
135 if (nremain) setup_reader(nremain,&arg[narg-nremain]);
136 else setup_reader(0,NULL);
137
138 // find the snapshot and read/bcast/process header info
139
140 if (me == 0 && screen) fprintf(screen,"Scanning dump file ...\n");
141
142 bigint ntimestep = seek(nstep,1);
143 if (ntimestep < 0)
144 error->all(FLERR,"Dump file does not contain requested snapshot");
145 header(1);
146
147 // reset timestep to nstep
148
149 update->reset_timestep(nstep);
150
151 // counters
152
153 // read in the snapshot and reset system
154
155 if (me == 0 && screen)
156 fprintf(screen,"Reading snapshot from dump file ...\n");
157
158 bigint natoms_prev = atom->natoms;
159 atoms();
160
161 if (me == 0) reader->close_file();
162
163 // print out stats
164
165 bigint npurge_all,nreplace_all,ntrim_all,nadd_all;
166
167 bigint tmp;
168 tmp = npurge;
169 MPI_Allreduce(&tmp,&npurge_all,1,MPI_LMP_BIGINT,MPI_SUM,world);
170 tmp = nreplace;
171 MPI_Allreduce(&tmp,&nreplace_all,1,MPI_LMP_BIGINT,MPI_SUM,world);
172 tmp = ntrim;
173 MPI_Allreduce(&tmp,&ntrim_all,1,MPI_LMP_BIGINT,MPI_SUM,world);
174 tmp = nadd;
175 MPI_Allreduce(&tmp,&nadd_all,1,MPI_LMP_BIGINT,MPI_SUM,world);
176
177 domain->print_box(" ");
178
179 if (me == 0) {
180 if (screen) {
181 fprintf(screen," " BIGINT_FORMAT " atoms before read\n",natoms_prev);
182 fprintf(screen," " BIGINT_FORMAT " atoms in snapshot\n",nsnapatoms);
183 fprintf(screen," " BIGINT_FORMAT " atoms purged\n",npurge_all);
184 fprintf(screen," " BIGINT_FORMAT " atoms replaced\n",nreplace_all);
185 fprintf(screen," " BIGINT_FORMAT " atoms trimmed\n",ntrim_all);
186 fprintf(screen," " BIGINT_FORMAT " atoms added\n",nadd_all);
187 fprintf(screen," " BIGINT_FORMAT " atoms after read\n",atom->natoms);
188 }
189 if (logfile) {
190 fprintf(logfile," " BIGINT_FORMAT " atoms before read\n",natoms_prev);
191 fprintf(logfile," " BIGINT_FORMAT " atoms in snapshot\n",nsnapatoms);
192 fprintf(logfile," " BIGINT_FORMAT " atoms purged\n",npurge_all);
193 fprintf(logfile," " BIGINT_FORMAT " atoms replaced\n",nreplace_all);
194 fprintf(logfile," " BIGINT_FORMAT " atoms trimmed\n",ntrim_all);
195 fprintf(logfile," " BIGINT_FORMAT " atoms added\n",nadd_all);
196 fprintf(logfile," " BIGINT_FORMAT " atoms after read\n",atom->natoms);
197 }
198 }
199 }
200
201 /* ---------------------------------------------------------------------- */
202
store_files(int nstr,char ** str)203 void ReadDump::store_files(int nstr, char **str)
204 {
205
206 if(strrchr(str[0],'*'))
207 {
208 file_search(str[0]);
209 return;
210 }
211
212 nfile = nstr;
213 files = new char*[nfile];
214
215 for (int i = 0; i < nfile; i++) {
216 int n = strlen(str[i]) + 1;
217 files[i] = new char[n];
218 strcpy(files[i],str[i]);
219 }
220 }
221
222 /* ----------------------------------------------------------------------
223 infile contains a "*"
224 search for all files which match the infile pattern
225 replace "*" with any timestep value
226 search dir referenced by initial pathname of file
227 ------------------------------------------------------------------------- */
228
file_search(char * infile)229 void ReadDump::file_search(char *infile)
230 {
231 char *ptr;
232
233 files = new char*[10000];
234 int *middle_index = new int[10000];
235
236 // separate infile into dir + filename
237
238 char *dirname = new char[strlen(infile) + 1];
239 char *filename = new char[strlen(infile) + 1];
240
241 if (strchr(infile,'/')) {
242 ptr = strrchr(infile,'/');
243 *ptr = '\0';
244 strcpy(dirname,infile);
245 strcpy(filename,ptr+1);
246 *ptr = '/';
247 } else {
248 strcpy(dirname,"./");
249 strcpy(filename,infile);
250 }
251
252 char *pattern = new char[strlen(filename) + 1];
253 strcpy(pattern,filename);
254
255 // scan all files in directory, searching for files that match pattern
256 // maxnum = largest int that matches "*"
257
258 size_t n = strlen(pattern) + 16;
259 char *begin = new char[n];
260 char *middle = new char[n];
261 char *end = new char[n];
262
263 ptr = strchr(pattern,'*');
264 *ptr = '\0';
265 strcpy(begin,pattern);
266 strcpy(end,ptr+1);
267 int nbegin = strlen(begin);
268 nfile = 0;
269
270 struct dirent *ep;
271 DIR *dp = opendir(dirname);
272 if (dp == NULL)
273 error->one(FLERR,"Cannot open dir to search for dump file");
274
275 while ((ep = readdir(dp))) {
276 if (strstr(ep->d_name,begin) != ep->d_name) continue;
277 if ((ptr = strstr(&ep->d_name[nbegin],end)) == NULL) continue;
278 if (strlen(end) == 0) ptr = ep->d_name + strlen(ep->d_name);
279 *ptr = '\0';
280 if (strlen(&ep->d_name[nbegin]) < n) {
281 strcpy(middle,&ep->d_name[nbegin]);
282
283 nfile++;
284
285 if(nfile >= 10000)
286 error->one(FLERR,"Currently max. 10000 dump files matching pattern can be read");
287
288 files[nfile-1] = new char[strlen(filename) + 16];
289 middle_index[nfile-1] = atoi(middle);
290 sprintf(files[nfile-1],"%s/%s%s%s",dirname,begin,middle,end);
291
292 }
293 }
294 closedir(dp);
295 if (nfile <= 0) error->one(FLERR,"Found no dump file matching pattern");
296
297 bool swaped;
298 int nswaps_left = nfile;
299 do
300 {
301 swaped = false;
302 for(int i = 0; i < nfile-1; i++)
303 {
304 if(middle_index[i] > middle_index[i+1])
305 {
306 //swap
307 char swapper[512];
308 strcpy(swapper,files[i+1]);
309 strcpy(files[i+1],files[i]);
310 strcpy(files[i],swapper);
311
312 int swapper_i = middle_index[i+1];
313 middle_index[i+1] = middle_index[i];
314 middle_index[i] = swapper_i;
315
316 swaped = true;
317 }
318 }
319 nswaps_left--;
320 } while(swaped && nswaps_left > 0);
321
322 // clean up
323
324 delete [] dirname;
325 delete [] filename;
326 delete [] pattern;
327 delete [] begin;
328 delete [] middle;
329 delete [] end;
330 delete [] middle_index;
331 }
332
333 /* ---------------------------------------------------------------------- */
334
setup_reader(int narg,char ** arg)335 void ReadDump::setup_reader(int narg, char **arg)
336 {
337 // allocate snapshot field buffer
338
339 memory->create(fields,CHUNK,nfield,"read_dump:fields");
340
341 // create reader class
342 // match readerstyle to options in style_reader.h
343
344 if (0) return; // dummy line to enable else-if macro expansion
345
346 #define READER_CLASS
347 #define ReaderStyle(key,Class) \
348 else if (strcmp(readerstyle,#key) == 0) reader = new Class(lmp);
349 #include "style_reader.h"
350 #undef READER_CLASS
351
352 // unrecognized style
353
354 else error->all(FLERR,"Invalid dump reader style");
355
356 // pass any arguments to reader
357
358 if (narg > 0) reader->settings(narg,arg);
359 }
360
361 /* ----------------------------------------------------------------------
362 seek Nrequest timestep in one or more dump files
363 if exact = 1, must find exactly Nrequest
364 if exact = 0, find first step >= Nrequest
365 return matching ntimestep or -1 if did not find a match
366 ------------------------------------------------------------------------- */
367
seek(bigint nrequest,int exact)368 bigint ReadDump::seek(bigint nrequest, int exact)
369 {
370 int ifile,eofflag;
371 bigint ntimestep = 0;
372
373 if (me == 0) {
374
375 // exit file loop when dump timestep >= nrequest
376 // or files exhausted
377
378 for (ifile = 0; ifile < nfile; ifile++) {
379 ntimestep = -1;
380 reader->open_file(files[ifile]);
381 while (1) {
382 eofflag = reader->read_time(ntimestep);
383 if (eofflag) break;
384 if (ntimestep >= nrequest) break;
385 reader->skip();
386 }
387 if (ntimestep >= nrequest) break;
388 reader->close_file();
389 }
390
391 currentfile = ifile;
392 if (ntimestep < nrequest) ntimestep = -1;
393 if (exact && ntimestep != nrequest) ntimestep = -1;
394 if (ntimestep < 0) reader->close_file();
395 }
396
397 MPI_Bcast(&ntimestep,1,MPI_LMP_BIGINT,0,world);
398 return ntimestep;
399 }
400
401 /* ----------------------------------------------------------------------
402 find next matching snapshot in one or more dump files
403 Ncurrent = current timestep from last snapshot
404 Nlast = match no timestep bigger than Nlast
405 Nevery = only match timesteps that are a multiple of Nevery
406 Nskip = skip every this many timesteps
407 return matching ntimestep or -1 if did not find a match
408 ------------------------------------------------------------------------- */
409
next(bigint ncurrent,bigint nlast,int nevery,int nskip)410 bigint ReadDump::next(bigint ncurrent, bigint nlast, int nevery, int nskip)
411 {
412 int ifile,eofflag=0;
413 bigint ntimestep = 0;
414
415 if (me == 0) {
416
417 // exit file loop when dump timestep matches all criteria
418 // or files exhausted
419
420 int iskip = 0;
421
422 for (ifile = currentfile; ifile < nfile; ifile++) {
423 ntimestep = -1;
424 if (ifile != currentfile) reader->open_file(files[ifile]);
425 while (1) {
426 eofflag = reader->read_time(ntimestep);
427 if (iskip == nskip) iskip = 0;
428 iskip++;
429 if (eofflag) break;
430 if (ntimestep <= ncurrent) break;
431 if (ntimestep > nlast) break;
432 if (nevery && ntimestep % nevery) reader->skip();
433 else if (iskip < nskip) reader->skip();
434 else break;
435 }
436 if (eofflag) reader->close_file();
437 else break;
438 }
439
440 currentfile = ifile;
441 if (eofflag) ntimestep = -1;
442 if (ntimestep <= ncurrent) ntimestep = -1;
443 if (ntimestep > nlast) ntimestep = -1;
444 if (ntimestep < 0) reader->close_file();
445 }
446
447 MPI_Bcast(&ntimestep,1,MPI_LMP_BIGINT,0,world);
448 return ntimestep;
449 }
450
451 /* ----------------------------------------------------------------------
452 read and broadcast and store snapshot header info
453 set nsnapatoms = # of atoms in snapshot
454 ------------------------------------------------------------------------- */
455
header(int fieldinfo)456 void ReadDump::header(int fieldinfo)
457 {
458 int triclinic_snap;
459 int fieldflag,xflag,yflag,zflag;
460
461 if (me == 0)
462 nsnapatoms = reader->read_header(box,triclinic_snap,
463 fieldinfo,nfield,fieldtype,fieldlabel,
464 scaleflag,wrapflag,fieldflag,
465 xflag,yflag,zflag);
466
467 MPI_Bcast(&nsnapatoms,1,MPI_LMP_BIGINT,0,world);
468 MPI_Bcast(&triclinic_snap,1,MPI_INT,0,world);
469 MPI_Bcast(&box[0][0],9,MPI_DOUBLE,0,world);
470
471 // local copy of snapshot box parameters
472 // used in xfield,yfield,zfield when converting dump atom to absolute coords
473
474 xlo = box[0][0];
475 xhi = box[0][1];
476 ylo = box[1][0];
477 yhi = box[1][1];
478 zlo = box[2][0];
479 zhi = box[2][1];
480 if (triclinic_snap) {
481 xy = box[0][2];
482 xz = box[1][2];
483 yz = box[2][2];
484 double xdelta = MIN(0.0,xy);
485 xdelta = MIN(xdelta,xz);
486 xdelta = MIN(xdelta,xy+xz);
487 xlo = xlo - xdelta;
488 xdelta = MAX(0.0,xy);
489 xdelta = MAX(xdelta,xz);
490 xdelta = MAX(xdelta,xy+xz);
491 xhi = xhi - xdelta;
492 ylo = ylo - MIN(0.0,yz);
493 yhi = yhi - MAX(0.0,yz);
494 }
495 xprd = xhi - xlo;
496 yprd = yhi - ylo;
497 zprd = zhi - zlo;
498
499 // done if not checking fields
500
501 if (!fieldinfo) return;
502
503 MPI_Bcast(&fieldflag,1,MPI_INT,0,world);
504 MPI_Bcast(&xflag,1,MPI_INT,0,world);
505 MPI_Bcast(&yflag,1,MPI_INT,0,world);
506 MPI_Bcast(&zflag,1,MPI_INT,0,world);
507
508 // error check on current vs new box and fields
509 // triclinic_snap < 0 means no box info in file
510
511 if (triclinic_snap < 0 && boxflag > 0)
512 error->all(FLERR,"No box information in dump. You have to use 'box no'");
513 if (triclinic_snap >= 0) {
514 if ((triclinic_snap && !triclinic) ||
515 (!triclinic_snap && triclinic))
516 error->one(FLERR,"Read_dump triclinic status does not match simulation");
517 }
518
519 // error check on requested fields exisiting in dump file
520
521 if (fieldflag < 0)
522 error->one(FLERR,"Read_dump field not found in dump file");
523
524 // all explicitly requested x,y,z must have consistent scaling & wrapping
525
526 int value = MAX(xflag,yflag);
527 value = MAX(zflag,value);
528 if ((xflag != UNSET && xflag != value) ||
529 (yflag != UNSET && yflag != value) ||
530 (zflag != UNSET && zflag != value))
531 error->one(FLERR,
532 "Read_dump xyz fields do not have consistent scaling/wrapping");
533
534 // set scaled/wrapped based on xyz flags
535
536 value = UNSET;
537 if (xflag != UNSET) value = xflag;
538 if (yflag != UNSET) value = yflag;
539 if (zflag != UNSET) value = zflag;
540
541 if (value == UNSET) {
542 scaled = wrapped = 0;
543 } else if (value == NOSCALE_NOWRAP) {
544 scaled = wrapped = 0;
545 } else if (value == NOSCALE_WRAP) {
546 scaled = 0;
547 wrapped = 1;
548 } else if (value == SCALE_NOWRAP) {
549 scaled = 1;
550 wrapped = 0;
551 } else if (value == SCALE_WRAP) {
552 scaled = wrapped = 1;
553 }
554
555 // scaled, triclinic coords require all 3 x,y,z fields, to perform unscaling
556 // set yindex,zindex = column index of Y and Z fields in fields array
557 // needed for unscaling to absolute coords in xfield(), yfield(), zfield()
558
559 if (scaled && triclinic == 1) {
560 int flag = 0;
561 if (xflag == UNSET) flag = 1;
562 if (yflag == UNSET) flag = 1;
563 if (dimension == 3 && zflag == UNSET) flag = 1;
564 if (flag)
565 error->one(FLERR,"All read_dump x,y,z fields must be specified for "
566 "scaled, triclinic coords");
567
568 for (int i = 0; i < nfield; i++) {
569 if (fieldtype[i] == Y) yindex = i;
570 if (fieldtype[i] == Z) zindex = i;
571 }
572 }
573 }
574
575 /* ---------------------------------------------------------------------- */
576
atoms()577 void ReadDump::atoms()
578 {
579 // initialize counters
580
581 npurge = nreplace = ntrim = nadd = 0;
582
583 // if purgeflag set, delete all current atoms
584
585 if (purgeflag) {
586 if (atom->map_style) atom->map_clear();
587 npurge = atom->nlocal;
588 atom->nlocal = atom->nghost = 0;
589 atom->natoms = 0;
590 }
591
592 // to match existing atoms to dump atoms:
593 // must build map if not a molecular system
594
595 int mapflag = 0;
596 if (atom->map_style == 0) {
597 mapflag = 1;
598 atom->map_style = 1;
599 atom->map_init();
600 atom->map_set();
601 }
602
603 // uflag[i] = 1 for each owned atom appearing in dump
604 // ucflag = similar flag for each chunk atom, used in process_atoms()
605
606 int nlocal = atom->nlocal;
607 memory->create(uflag,nlocal,"read_dump:uflag");
608 for (int i = 0; i < nlocal; i++) uflag[i] = 0;
609 memory->create(ucflag,CHUNK,"read_dump:ucflag");
610 memory->create(ucflag_all,CHUNK,"read_dump:ucflag");
611
612 // read, broadcast, and process atoms from snapshot in chunks
613
614 addproc = -1;
615
616 int nchunk;
617 bigint nread = 0;
618 while (nread < nsnapatoms) {
619 nchunk = MIN(nsnapatoms-nread,CHUNK);
620 if (me == 0) reader->read_atoms(nchunk,nfield,fields);
621 MPI_Bcast(&fields[0][0],nchunk*nfield,MPI_DOUBLE,0,world);
622 process_atoms(nchunk);
623 nread += nchunk;
624 }
625
626 // if addflag set, add tags to new atoms if possible
627
628 if (addflag) {
629 bigint nblocal = atom->nlocal;
630 MPI_Allreduce(&nblocal,&atom->natoms,1,MPI_LMP_BIGINT,MPI_SUM,world);
631 if (atom->natoms < 0 || atom->natoms > MAXBIGINT)
632 error->all(FLERR,"Too many total atoms");
633 // change these to MAXTAGINT when allow tagint = bigint
634 if (atom->natoms > MAXSMALLINT) atom->tag_enable = 0;
635 if (atom->natoms <= MAXSMALLINT) atom->tag_extend();
636 }
637
638 // if trimflag set, delete atoms not replaced by snapshot atoms
639
640 if (trimflag) {
641 delete_atoms();
642 bigint nblocal = atom->nlocal;
643 MPI_Allreduce(&nblocal,&atom->natoms,1,MPI_LMP_BIGINT,MPI_SUM,world);
644 }
645
646 // can now delete uflag arrays
647
648 memory->destroy(uflag);
649 memory->destroy(ucflag);
650 memory->destroy(ucflag_all);
651
652 // delete atom map if created it above
653 // else reinitialize map for current atoms
654 // do this before migrating atoms to new procs via Irregular
655
656 if (mapflag) {
657 atom->map_delete();
658 atom->map_style = 0;
659 } else {
660 atom->nghost = 0;
661 atom->map_init();
662 atom->map_set();
663 }
664
665 // overwrite simulation box with dump snapshot box if requested
666 // reallocate processors to box
667
668 if (boxflag) {
669 domain->boxlo[0] = xlo;
670 domain->boxhi[0] = xhi;
671 domain->boxlo[1] = ylo;
672 domain->boxhi[1] = yhi;
673 if (dimension == 3) {
674 domain->boxlo[2] = zlo;
675 domain->boxhi[2] = zhi;
676 }
677 if (triclinic) {
678 domain->xy = xy;
679 if (dimension == 3) {
680 domain->xz = xz;
681 domain->yz = yz;
682 }
683 }
684
685 domain->set_initial_box();
686 domain->set_global_box();
687 comm->set_proc_grid(0);
688 domain->set_local_box();
689 }
690
691 // move atoms back inside simulation box and to new processors
692 // use remap() instead of pbc() in case atoms moved a long distance
693 // adjust image flags of all atoms (old and new) based on current box
694 // use irregular() in case atoms moved a long distance
695
696 double **x = atom->x;
697 tagint *image = atom->image;
698 nlocal = atom->nlocal;
699 for (int i = 0; i < nlocal; i++) domain->remap(x[i],image[i]);
700
701 if (triclinic) domain->x2lamda(atom->nlocal);
702 domain->reset_box();
703 Irregular *irregular = new Irregular(lmp);
704 irregular->migrate_atoms();
705 delete irregular;
706 if (triclinic) domain->lamda2x(atom->nlocal);
707 }
708
709 /* ----------------------------------------------------------------------
710 process arg list for dump file fields and optional keywords
711 ------------------------------------------------------------------------- */
712
fields_and_keywords(int narg,char ** arg)713 int ReadDump::fields_and_keywords(int narg, char **arg)
714 {
715 // per-field vectors, leave space for ID and TYPE
716
717 fieldtype = new int[narg+2];
718 fieldlabel = new char*[narg+2];
719
720 // add id and type fields as needed
721 // scan ahead to see if "add yes" keyword/value is used
722 // requires extra "type" field from from dump file
723
724 int iarg;
725 for (iarg = 0; iarg < narg; iarg++)
726 if (strcmp(arg[iarg],"add") == 0)
727 if (iarg < narg-1 && strcmp(arg[iarg+1],"yes") == 0) break;
728
729 nfield = 0;
730 fieldtype[nfield++] = ID;
731 /*if (iarg < narg)*/ fieldtype[nfield++] = TYPE;
732
733 bool radius_present = false;
734
735 // parse fields
736
737 iarg = 0;
738 while (iarg < narg) {
739 if (strcmp(arg[iarg],"x") == 0) fieldtype[nfield++] = X;
740 else if (strcmp(arg[iarg],"y") == 0) fieldtype[nfield++] = Y;
741 else if (strcmp(arg[iarg],"z") == 0) fieldtype[nfield++] = Z;
742 else if (strcmp(arg[iarg],"vx") == 0) fieldtype[nfield++] = VX;
743 else if (strcmp(arg[iarg],"vy") == 0) fieldtype[nfield++] = VY;
744 else if (strcmp(arg[iarg],"vz") == 0) fieldtype[nfield++] = VZ;
745 else if (strcmp(arg[iarg],"omegax") == 0) {
746 if (!atom->omega_flag)
747 error->all(FLERR,"Read dump of atom property omegax that isn't allocated");
748 fieldtype[nfield++] = OMEGAX;
749 }
750 else if (strcmp(arg[iarg],"omegay") == 0) {
751 if (!atom->omega_flag)
752 error->all(FLERR,"Read dump of atom property omegay that isn't allocated");
753 fieldtype[nfield++] = OMEGAY;
754 }
755 else if (strcmp(arg[iarg],"omegaz") == 0) {
756 if (!atom->omega_flag)
757 error->all(FLERR,"Read dump of atom property omegaz that isn't allocated");
758 fieldtype[nfield++] = OMEGAZ;
759 }
760 else if (strcmp(arg[iarg],"q") == 0) {
761 if (!atom->q_flag)
762 error->all(FLERR,"Read dump of atom property that isn't allocated");
763 fieldtype[nfield++] = Q;
764 }
765 else if (strcmp(arg[iarg],"radius") == 0) {
766 if (!atom->radius_flag)
767 error->all(FLERR,"Read dump of atom property that isn't allocated");
768 fieldtype[nfield++] = RADIUS;
769 radius_present = true;
770 }
771 else if (strcmp(arg[iarg],"mass") == 0) {
772 if (!atom->rmass_flag)
773 error->all(FLERR,"Read dump of atom property that isn't allocated");
774 fieldtype[nfield++] = MASS;
775 }
776 else if (strcmp(arg[iarg],"density") == 0) {
777 if (!atom->density_flag)
778 error->all(FLERR,"Read dump of atom property that isn't allocated");
779 fieldtype[nfield++] = DENSITY;
780 }
781 else if (strcmp(arg[iarg],"ix") == 0) fieldtype[nfield++] = IX;
782 else if (strcmp(arg[iarg],"iy") == 0) fieldtype[nfield++] = IY;
783 else if (strcmp(arg[iarg],"iz") == 0) fieldtype[nfield++] = IZ;
784 else if (strcmp(arg[iarg],"fx") == 0) fieldtype[nfield++] = FX;
785 else if (strcmp(arg[iarg],"fy") == 0) fieldtype[nfield++] = FY;
786 else if (strcmp(arg[iarg],"fz") == 0) fieldtype[nfield++] = FZ;
787 else break;
788 iarg++;
789 }
790
791 if (!radius_present) error->all(FLERR,"read_dump: the radius field is required to be read");
792
793 // check for no fields
794
795 if (fieldtype[nfield-1] == ID || fieldtype[nfield-1] == TYPE)
796 error->all(FLERR,"Illegal read_dump command");
797
798 if (dimension == 2) {
799 for (int i = 0; i < nfield; i++)
800 if (fieldtype[i] == Z || fieldtype[i] == VZ || fieldtype[i] == IZ)
801 error->all(FLERR,"Illegal read_dump command");
802 }
803
804 for (int i = 0; i < nfield; i++)
805 for (int j = i+1; j < nfield; j++)
806 if (fieldtype[i] == fieldtype[j])
807 error->all(FLERR,"Duplicate fields in read_dump command");
808
809 // parse optional args
810
811 boxflag = 1;
812 replaceflag = 1;
813 purgeflag = 0;
814 trimflag = 0;
815 addflag = 0;
816 for (int i = 0; i < nfield; i++) fieldlabel[i] = NULL;
817 scaleflag = 0;
818 wrapflag = 1;
819
820 while (iarg < narg) {
821 if (strcmp(arg[iarg],"box") == 0) {
822 if (iarg+2 > narg) error->all(FLERR,"Illegal read_dump command");
823 if (strcmp(arg[iarg+1],"yes") == 0) boxflag = 1;
824 else if (strcmp(arg[iarg+1],"no") == 0) boxflag = 0;
825 else error->all(FLERR,"Illegal read_dump command");
826 iarg += 2;
827 } else if (strcmp(arg[iarg],"replace") == 0) {
828 if (iarg+2 > narg) error->all(FLERR,"Illegal read_dump command");
829 if (strcmp(arg[iarg+1],"yes") == 0) replaceflag = 1;
830 else if (strcmp(arg[iarg+1],"no") == 0) replaceflag = 0;
831 else error->all(FLERR,"Illegal read_dump command");
832 iarg += 2;
833 } else if (strcmp(arg[iarg],"purge") == 0) {
834 if (iarg+2 > narg) error->all(FLERR,"Illegal read_dump command");
835 if (strcmp(arg[iarg+1],"yes") == 0) purgeflag = 1;
836 else if (strcmp(arg[iarg+1],"no") == 0) purgeflag = 0;
837 else error->all(FLERR,"Illegal read_dump command");
838 iarg += 2;
839 } else if (strcmp(arg[iarg],"trim") == 0) {
840 if (iarg+2 > narg) error->all(FLERR,"Illegal read_dump command");
841 if (strcmp(arg[iarg+1],"yes") == 0) trimflag = 1;
842 else if (strcmp(arg[iarg+1],"no") == 0) trimflag = 0;
843 else error->all(FLERR,"Illegal read_dump command");
844 iarg += 2;
845 } else if (strcmp(arg[iarg],"add") == 0) {
846 if (iarg+2 > narg) error->all(FLERR,"Illegal read_dump command");
847 if (strcmp(arg[iarg+1],"yes") == 0) addflag = 1;
848 else if (strcmp(arg[iarg+1],"no") == 0) addflag = 0;
849 else error->all(FLERR,"Illegal read_dump command");
850 iarg += 2;
851 } else if (strcmp(arg[iarg],"label") == 0) {
852 if (iarg+3 > narg) error->all(FLERR,"Illegal read_dump command");
853 int i;
854 for (i = 0; i < nfield; i++)
855 if (fieldlabel[i] && strcmp(arg[iarg+1],fieldlabel[i]) == 0) break;
856 if (i == nfield) error->all(FLERR,"Illegal read_dump command");
857 int n = strlen(arg[iarg+2]) + 1;
858 fieldlabel[i] = new char[n];
859 strcpy(fieldlabel[i],arg[iarg+2]);
860 iarg += 3;
861 } else if (strcmp(arg[iarg],"scaled") == 0) {
862 if (iarg+2 > narg) error->all(FLERR,"Illegal read_dump command");
863 if (strcmp(arg[iarg+1],"yes") == 0) scaleflag = 1;
864 else if (strcmp(arg[iarg+1],"no") == 0) scaleflag = 0;
865 else error->all(FLERR,"Illegal read_dump command");
866 iarg += 2;
867 } else if (strcmp(arg[iarg],"wrapped") == 0) {
868 if (iarg+2 > narg) error->all(FLERR,"Illegal read_dump command");
869 if (strcmp(arg[iarg+1],"yes") == 0) wrapflag = 1;
870 else if (strcmp(arg[iarg+1],"no") == 0) wrapflag = 0;
871 else error->all(FLERR,"Illegal read_dump command");
872 iarg += 2;
873 } else if (strcmp(arg[iarg],"format") == 0) {
874 if (iarg+2 > narg) error->all(FLERR,"Illegal read_dump command");
875 delete [] readerstyle;
876 int n = strlen(arg[iarg+1]) + 1;
877 readerstyle = new char[n];
878 strcpy(readerstyle,arg[iarg+1]);
879 iarg += 2;
880 break;
881 } else {
882
883 error->all(FLERR,"Illegal read_dump command");
884 }
885 }
886
887 if (purgeflag && (replaceflag || trimflag))
888 error->all(FLERR,"If read_dump purges it cannot replace or trim");
889
890 return narg-iarg;
891 }
892
893 /* ----------------------------------------------------------------------
894 process each of N atoms in chunk read from dump file
895 if in replace mode and atom ID matches current atom,
896 overwrite atom info with fields from dump file
897 if in add mode and atom ID does not match any current atom,
898 create new atom with dump file field values,
899 and assign to a proc in round-robin manner
900 use round-robin method, b/c atom coords may not be inside simulation box
901 ------------------------------------------------------------------------- */
902
process_atoms(int n)903 void ReadDump::process_atoms(int n)
904 {
905 int i,m,ifield,itype=0,itag;;
906 int xbox,ybox,zbox;
907
908 double **x = atom->x;
909 double **v = atom->v;
910 double **f = atom->f;
911 double **omega = atom->omega;
912 double *q = atom->q;
913 double *radius = atom->radius;
914 double *rmass= atom->rmass;
915 double *density= atom->density;
916 tagint *image = atom->image;
917 int nlocal = atom->nlocal;
918 int map_tag_max = atom->map_tag_max;
919
920 for (i = 0; i < n; i++) {
921 ucflag[i] = 0;
922
923 // check if new atom matches one I own
924 // setting m = -1 forces new atom not to match
925
926 itag = static_cast<int> (fields[i][0]);
927 if (itag <= map_tag_max) m = atom->map(static_cast<int> (fields[i][0]));
928 else m = -1;
929 if (m < 0 || m >= nlocal) continue;
930
931 ucflag[i] = 1;
932 uflag[m] = 1;
933
934 if (replaceflag) {
935 nreplace++;
936
937 // current image flags
938
939 xbox = (image[m] & IMGMASK) - IMGMAX;
940 ybox = (image[m] >> IMGBITS & IMGMASK) - IMGMAX;
941 zbox = (image[m] >> IMG2BITS) - IMGMAX;
942
943 // overwrite atom attributes with field info
944 // start from field 1 since 0 = id, 1 will be skipped if type
945
946 for (ifield = 1; ifield < nfield; ifield++) {
947 switch (fieldtype[ifield]) {
948 case X:
949 x[m][0] = xfield(i,ifield);
950 break;
951 case Y:
952 x[m][1] = yfield(i,ifield);
953 break;
954 case Z:
955 x[m][2] = zfield(i,ifield);
956 break;
957 case VX:
958 v[m][0] = fields[i][ifield];
959 break;
960 case VY:
961 v[m][1] = fields[i][ifield];
962 break;
963 case VZ:
964 v[m][2] = fields[i][ifield];
965 break;
966 case OMEGAX:
967 omega[m][0] = fields[i][ifield];
968 break;
969 case OMEGAY:
970 omega[m][1] = fields[i][ifield];
971 break;
972 case OMEGAZ:
973 omega[m][2] = fields[i][ifield];
974 break;
975 case Q:
976 q[m] = fields[i][ifield];
977 break;
978 case RADIUS:
979 radius[m] = fields[i][ifield];
980 break;
981 case MASS:
982 rmass[m] = fields[i][ifield];
983 break;
984 case DENSITY:
985 density[m] = fields[i][ifield];
986 break;
987 case IX:
988 xbox = static_cast<int> (fields[i][ifield]);
989 break;
990 case IY:
991 ybox = static_cast<int> (fields[i][ifield]);
992 break;
993 case IZ:
994 zbox = static_cast<int> (fields[i][ifield]);
995 break;
996 case FX:
997 f[m][0] = fields[i][ifield];
998 break;
999 case FY:
1000 f[m][1] = fields[i][ifield];
1001 break;
1002 case FZ:
1003 f[m][2] = fields[i][ifield];
1004 break;
1005 }
1006 }
1007
1008 // replace image flag in case changed by ix,iy,iz fields or unwrapping
1009
1010 if (!wrapped) xbox = ybox = zbox = 0;
1011
1012 image[m] = ((tagint) (xbox + IMGMAX) & IMGMASK) |
1013 (((tagint) (ybox + IMGMAX) & IMGMASK) << IMGBITS) |
1014 (((tagint) (zbox + IMGMAX) & IMGMASK) << IMG2BITS);
1015 }
1016 }
1017
1018 // create any atoms in chunk that no processor owned
1019 // add atoms in round-robin sequence on processors
1020 // cannot do it geometrically b/c dump coords may not be in simulation box
1021
1022 if (!addflag) return;
1023
1024 MPI_Allreduce(ucflag,ucflag_all,n,MPI_INT,MPI_SUM,world);
1025
1026 int nlocal_previous = atom->nlocal;
1027 double one[3];
1028
1029 for (i = 0; i < n; i++) {
1030 if (ucflag_all[i]) continue;
1031
1032 // each processor adds every Pth atom
1033
1034 addproc++;
1035 if (addproc == nprocs) addproc = 0;
1036 if (addproc != me) continue;
1037
1038 // create type and coord fields from dump file
1039 // coord = 0.0 unless corresponding dump file field was specified
1040
1041 one[0] = one[1] = one[2] = 0.0;
1042 for (ifield = 1; ifield < nfield; ifield++) {
1043 switch (fieldtype[ifield]) {
1044 case TYPE:
1045 itype = static_cast<int> (fields[i][ifield]);
1046 break;
1047 case X:
1048 one[0] = xfield(i,ifield);
1049 break;
1050 case Y:
1051 one[1] = yfield(i,ifield);
1052 break;
1053 case Z:
1054 one[2] = zfield(i,ifield);
1055 break;
1056 }
1057 }
1058
1059 // create the atom on proc that owns it
1060 // reset v,image ptrs in case they are reallocated
1061
1062 m = atom->nlocal;
1063 atom->avec->create_atom(itype,one);
1064 nadd++;
1065
1066 v = atom->v;
1067 f = atom->f;
1068 omega = atom->omega;
1069 q = atom->q;
1070 radius = atom->radius;
1071 rmass = atom->rmass;
1072 density = atom->density;
1073 image = atom->image;
1074
1075 // set atom attributes from other dump file fields
1076
1077 xbox = ybox = zbox = 0;
1078
1079 for (ifield = 1; ifield < nfield; ifield++) {
1080 switch (fieldtype[ifield]) {
1081 case VX:
1082 v[m][0] = fields[i][ifield];
1083 break;
1084 case VY:
1085 v[m][1] = fields[i][ifield];
1086 break;
1087 case VZ:
1088 v[m][2] = fields[i][ifield];
1089 break;
1090 case OMEGAX:
1091 omega[m][0] = fields[i][ifield];
1092 break;
1093 case OMEGAY:
1094 omega[m][1] = fields[i][ifield];
1095 break;
1096 case OMEGAZ:
1097 omega[m][2] = fields[i][ifield];
1098 break;
1099 case Q:
1100 q[m] = fields[i][ifield];
1101 break;
1102 case RADIUS:
1103 radius[m] = fields[i][ifield];
1104 break;
1105 case MASS:
1106 rmass[m] = fields[i][ifield];
1107 break;
1108 case DENSITY:
1109 density[m] = fields[i][ifield];
1110 break;
1111 case IX:
1112 xbox = static_cast<int> (fields[i][ifield]);
1113 break;
1114 case IY:
1115 ybox = static_cast<int> (fields[i][ifield]);
1116 break;
1117 case IZ:
1118 zbox = static_cast<int> (fields[i][ifield]);
1119 break;
1120 case FX:
1121 f[m][0] = fields[i][ifield];
1122 break;
1123 case FY:
1124 f[m][1] = fields[i][ifield];
1125 break;
1126 case FZ:
1127 f[m][2] = fields[i][ifield];
1128 break;
1129 }
1130
1131 // replace image flag in case changed by ix,iy,iz fields
1132
1133 image[m] = ((tagint) (xbox + IMGMAX) & IMGMASK) |
1134 (((tagint) (ybox + IMGMAX) & IMGMASK) << IMGBITS) |
1135 (((tagint) (zbox + IMGMAX) & IMGMASK) << IMG2BITS);
1136 }
1137 }
1138
1139 // invoke set_arrays() for fixes that need initialization of new atoms
1140 // same as in CreateAtoms
1141
1142 nlocal = atom->nlocal;
1143 for (m = 0; m < modify->nfix; m++) {
1144 Fix *fix = modify->fix[m];
1145 if (fix->create_attribute)
1146 {
1147 fix->pre_set_arrays();
1148 for (i = nlocal_previous; i < nlocal; i++)
1149 fix->set_arrays(i);
1150 }
1151 }
1152 }
1153
1154 /* ----------------------------------------------------------------------
1155 delete atoms not flagged as replaced by dump atoms
1156 ------------------------------------------------------------------------- */
1157
delete_atoms()1158 void ReadDump::delete_atoms()
1159 {
1160 AtomVec *avec = atom->avec;
1161 int nlocal = atom->nlocal;
1162
1163 int i = 0;
1164 while (i < nlocal) {
1165 if (uflag[i] == 0) {
1166 avec->copy(nlocal-1,i,1);
1167 uflag[i] = uflag[nlocal-1];
1168 nlocal--;
1169 ntrim++;
1170 } else i++;
1171 }
1172
1173 atom->nlocal = nlocal;
1174 }
1175
1176 /* ----------------------------------------------------------------------
1177 convert XYZ fields in dump file into absolute, unscaled coordinates
1178 depends on scaled vs unscaled and triclinic vs orthogonal
1179 does not depend on wrapped vs unwrapped
1180 ------------------------------------------------------------------------- */
1181
xfield(int i,int j)1182 double ReadDump::xfield(int i, int j)
1183 {
1184 if (!scaled) return fields[i][j];
1185 else if (!triclinic) return fields[i][j]*xprd + xlo;
1186 else if (dimension == 2)
1187 return xprd*fields[i][j] + xy*fields[i][yindex] + xlo;
1188 return xprd*fields[i][j] + xy*fields[i][yindex] + xz*fields[i][zindex] + xlo;
1189 }
1190
yfield(int i,int j)1191 double ReadDump::yfield(int i, int j)
1192 {
1193 if (!scaled) return fields[i][j];
1194 else if (!triclinic) return fields[i][j]*yprd + ylo;
1195 else if (dimension == 2) return yprd*fields[i][j] + ylo;
1196 return yprd*fields[i][j] + yz*fields[i][zindex] + ylo;
1197 }
1198
zfield(int i,int j)1199 double ReadDump::zfield(int i, int j)
1200 {
1201 if (!scaled) return fields[i][j];
1202 return fields[i][j]*zprd + zlo;
1203 }
1204