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