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 "read_restart.h"
16 
17 #include "angle.h"
18 #include "atom.h"
19 #include "atom_vec.h"
20 #include "bond.h"
21 #include "comm.h"
22 #include "dihedral.h"
23 #include "domain.h"
24 #include "error.h"
25 #include "fix_read_restart.h"
26 #include "force.h"
27 #include "group.h"
28 #include "improper.h"
29 #include "irregular.h"
30 #include "memory.h"
31 #include "modify.h"
32 #include "mpiio.h"
33 #include "pair.h"
34 #include "special.h"
35 #include "update.h"
36 
37 #include <cstring>
38 #include <dirent.h>
39 
40 #include "lmprestart.h"
41 
42 using namespace LAMMPS_NS;
43 
44 /* ---------------------------------------------------------------------- */
45 
ReadRestart(LAMMPS * lmp)46 ReadRestart::ReadRestart(LAMMPS *lmp) : Command(lmp), mpiio(nullptr) {}
47 
48 /* ---------------------------------------------------------------------- */
49 
command(int narg,char ** arg)50 void ReadRestart::command(int narg, char **arg)
51 {
52   if (narg != 1 && narg != 2) error->all(FLERR,"Illegal read_restart command");
53 
54   if (domain->box_exist)
55     error->all(FLERR,"Cannot read_restart after simulation box is defined");
56 
57   MPI_Barrier(world);
58   double time1 = MPI_Wtime();
59 
60   MPI_Comm_rank(world,&me);
61   MPI_Comm_size(world,&nprocs);
62 
63   // check for remap option
64 
65   int remapflag = 0;
66   if (narg == 2) {
67     if (strcmp(arg[1],"remap") == 0) remapflag = 1;
68     else error->all(FLERR,"Illegal read_restart command");
69   }
70 
71   // if filename contains "*", search dir for latest restart file
72 
73   char *file;
74   if (strchr(arg[0],'*')) {
75     int n=0;
76     if (me == 0) {
77       auto fn = file_search(arg[0]);
78       n = fn.size()+1;
79       file = utils::strdup(fn);
80     }
81     MPI_Bcast(&n,1,MPI_INT,0,world);
82     if (me != 0) file = new char[n];
83     MPI_Bcast(file,n,MPI_CHAR,0,world);
84   } else file = utils::strdup(arg[0]);
85 
86   // check for multiproc files and an MPI-IO filename
87 
88   if (strchr(arg[0],'%')) multiproc = 1;
89   else multiproc = 0;
90   if (strstr(arg[0],".mpiio")) mpiioflag = 1;
91   else mpiioflag = 0;
92 
93   if (multiproc && mpiioflag)
94     error->all(FLERR,"Read restart MPI-IO input not allowed with % in filename");
95 
96   if (mpiioflag) {
97     mpiio = new RestartMPIIO(lmp);
98     if (!mpiio->mpiio_exists)
99       error->all(FLERR,"Reading from MPI-IO filename when MPIIO package is not installed");
100   }
101 
102   // open single restart file or base file for multiproc case
103 
104   if (me == 0) {
105     utils::logmesg(lmp,"Reading restart file ...\n");
106     std::string hfile = file;
107     if (multiproc) {
108       hfile.replace(hfile.find('%'),1,"base");
109     }
110     fp = fopen(hfile.c_str(),"rb");
111     if (fp == nullptr)
112       error->one(FLERR,"Cannot open restart file {}: {}", hfile, utils::getsyserror());
113   }
114 
115   // read magic string, endian flag, format revision
116 
117   magic_string();
118   endian();
119   format_revision();
120 
121   // read header info which creates simulation box
122 
123   header();
124   domain->box_exist = 1;
125 
126   // problem setup using info from header
127 
128   int n;
129   if (nprocs == 1) n = static_cast<int> (atom->natoms);
130   else n = static_cast<int> (LB_FACTOR * atom->natoms / nprocs);
131 
132   atom->allocate_type_arrays();
133   atom->deallocate_topology();
134 
135   // allocate atom arrays to size N, rounded up by AtomVec->DELTA
136 
137   bigint nbig = n;
138   nbig = atom->avec->roundup(nbig);
139   n = static_cast<int> (nbig);
140   atom->avec->grow(n);
141   n = atom->nmax;
142 
143   domain->print_box("  ");
144   domain->set_initial_box(0);
145   domain->set_global_box();
146   comm->set_proc_grid();
147   domain->set_local_box();
148 
149   // read groups, ntype-length arrays, force field, fix info from file
150   // nextra = max # of extra quantities stored with each atom
151 
152   group->read_restart(fp);
153   type_arrays();
154   force_fields();
155 
156   int nextra = modify->read_restart(fp);
157   atom->nextra_store = nextra;
158   memory->create(atom->extra,n,nextra,"atom:extra");
159 
160   // read file layout info
161 
162   file_layout();
163 
164   // close header file if in multiproc mode
165 
166   if (multiproc && me == 0) {
167     fclose(fp);
168     fp = nullptr;
169   }
170 
171   // read per-proc info
172 
173   AtomVec *avec = atom->avec;
174 
175   int maxbuf = 0;
176   double *buf = nullptr;
177   int m,flag;
178 
179   // MPI-IO input from single file
180 
181   if (mpiioflag) {
182     mpiio->openForRead(file);
183     memory->create(buf,assignedChunkSize,"read_restart:buf");
184     mpiio->read((headerOffset+assignedChunkOffset),assignedChunkSize,buf);
185     mpiio->close();
186 
187     // can calculate number of atoms from assignedChunkSize
188 
189     if (!nextra) {
190       atom->nlocal = 1; // temporarily claim there is one atom...
191       int perAtomSize = avec->size_restart(); // ...so we can get its size
192       atom->nlocal = 0; // restore nlocal to zero atoms
193       int atomCt = (int) (assignedChunkSize / perAtomSize);
194       if (atomCt > atom->nmax) avec->grow(atomCt);
195     }
196 
197     m = 0;
198     while (m < assignedChunkSize) m += avec->unpack_restart(&buf[m]);
199   }
200 
201   // input of single native file
202   // nprocs_file = # of chunks in file
203   // proc 0 reads a chunk and bcasts it to other procs
204   // each proc unpacks the atoms, saving ones in it's sub-domain
205   // if remapflag set, remap the atom to box before checking sub-domain
206   // check for atom in sub-domain differs for orthogonal vs triclinic box
207 
208   else if (multiproc == 0) {
209 
210     int triclinic = domain->triclinic;
211     imageint *iptr;
212     double *x,lamda[3];
213     double *coord,*sublo,*subhi;
214     if (triclinic == 0) {
215       sublo = domain->sublo;
216       subhi = domain->subhi;
217     } else {
218       sublo = domain->sublo_lamda;
219       subhi = domain->subhi_lamda;
220     }
221 
222     for (int iproc = 0; iproc < nprocs_file; iproc++) {
223       if (read_int() != PERPROC)
224         error->all(FLERR,"Invalid flag in peratom section of restart file");
225 
226       n = read_int();
227       if (n > maxbuf) {
228         maxbuf = n;
229         memory->destroy(buf);
230         memory->create(buf,maxbuf,"read_restart:buf");
231       }
232       read_double_vec(n,buf);
233 
234       m = 0;
235       while (m < n) {
236         x = &buf[m+1];
237         if (remapflag) {
238           iptr = (imageint *) &buf[m+7];
239           domain->remap(x,*iptr);
240         }
241 
242         if (triclinic) {
243           domain->x2lamda(x,lamda);
244           coord = lamda;
245         } else coord = x;
246 
247         if (coord[0] >= sublo[0] && coord[0] < subhi[0] &&
248             coord[1] >= sublo[1] && coord[1] < subhi[1] &&
249             coord[2] >= sublo[2] && coord[2] < subhi[2]) {
250           m += avec->unpack_restart(&buf[m]);
251         } else m += static_cast<int> (buf[m]);
252       }
253     }
254 
255     if (me == 0) {
256       fclose(fp);
257       fp = nullptr;
258     }
259   }
260 
261   // input of multiple native files with procs <= files
262   // # of files = multiproc_file
263   // each proc reads a subset of files, striding by nprocs
264   // each proc keeps all atoms in all perproc chunks in its files
265 
266   else if (nprocs <= multiproc_file) {
267 
268     for (int iproc = me; iproc < multiproc_file; iproc += nprocs) {
269       std::string procfile = file;
270       procfile.replace(procfile.find('%'),1,fmt::format("{}",iproc));
271       fp = fopen(procfile.c_str(),"rb");
272       if (fp == nullptr)
273         error->one(FLERR,"Cannot open restart file {}: {}",
274                                      procfile, utils::getsyserror());
275       utils::sfread(FLERR,&flag,sizeof(int),1,fp,nullptr,error);
276       if (flag != PROCSPERFILE)
277         error->one(FLERR,"Invalid flag in peratom section of restart file");
278       int procsperfile;
279       utils::sfread(FLERR,&procsperfile,sizeof(int),1,fp,nullptr,error);
280 
281       for (int i = 0; i < procsperfile; i++) {
282         utils::sfread(FLERR,&flag,sizeof(int),1,fp,nullptr,error);
283         if (flag != PERPROC)
284           error->one(FLERR,"Invalid flag in peratom section of restart file");
285 
286         utils::sfread(FLERR,&n,sizeof(int),1,fp,nullptr,error);
287         if (n > maxbuf) {
288           maxbuf = n;
289           memory->destroy(buf);
290           memory->create(buf,maxbuf,"read_restart:buf");
291         }
292         utils::sfread(FLERR,buf,sizeof(double),n,fp,nullptr,error);
293 
294         m = 0;
295         while (m < n) m += avec->unpack_restart(&buf[m]);
296       }
297 
298       fclose(fp);
299       fp = nullptr;
300     }
301   }
302 
303   // input of multiple native files with procs > files
304   // # of files = multiproc_file
305   // cluster procs based on # of files
306   // 1st proc in each cluster reads per-proc chunks from file
307   // sends chunks round-robin to other procs in its cluster
308   // each proc keeps all atoms in its perproc chunks in file
309 
310   else {
311 
312     // nclusterprocs = # of procs in my cluster that read from one file
313     // filewriter = 1 if this proc reads file, else 0
314     // fileproc = ID of proc in my cluster who reads from file
315     // clustercomm = MPI communicator within my cluster of procs
316 
317     int nfile = multiproc_file;
318     int icluster = static_cast<int> ((bigint) me * nfile/nprocs);
319     int fileproc = static_cast<int> ((bigint) icluster * nprocs/nfile);
320     int fcluster = static_cast<int> ((bigint) fileproc * nfile/nprocs);
321     if (fcluster < icluster) fileproc++;
322     int fileprocnext =
323       static_cast<int> ((bigint) (icluster+1) * nprocs/nfile);
324     fcluster = static_cast<int> ((bigint) fileprocnext * nfile/nprocs);
325     if (fcluster < icluster+1) fileprocnext++;
326     int nclusterprocs = fileprocnext - fileproc;
327     int filereader = 0;
328     if (me == fileproc) filereader = 1;
329     MPI_Comm clustercomm;
330     MPI_Comm_split(world,icluster,0,&clustercomm);
331 
332     if (filereader) {
333       std::string procfile = file;
334       procfile.replace(procfile.find('%'),1,fmt::format("{}",icluster));
335       fp = fopen(procfile.c_str(),"rb");
336       if (fp == nullptr)
337         error->one(FLERR,"Cannot open restart file {}: {}", procfile, utils::getsyserror());
338     }
339 
340     int procsperfile;
341 
342     if (filereader) {
343       utils::sfread(FLERR,&flag,sizeof(int),1,fp,nullptr,error);
344       if (flag != PROCSPERFILE)
345         error->one(FLERR,"Invalid flag in peratom section of restart file");
346       utils::sfread(FLERR,&procsperfile,sizeof(int),1,fp,nullptr,error);
347     }
348     MPI_Bcast(&procsperfile,1,MPI_INT,0,clustercomm);
349 
350     int tmp,iproc;
351     MPI_Request request;
352 
353     for (int i = 0; i < procsperfile; i++) {
354       if (filereader) {
355         utils::sfread(FLERR,&flag,sizeof(int),1,fp,nullptr,error);
356         if (flag != PERPROC)
357           error->one(FLERR,"Invalid flag in peratom section of restart file");
358 
359         utils::sfread(FLERR,&n,sizeof(int),1,fp,nullptr,error);
360         if (n > maxbuf) {
361           maxbuf = n;
362           memory->destroy(buf);
363           memory->create(buf,maxbuf,"read_restart:buf");
364         }
365         utils::sfread(FLERR,buf,sizeof(double),n,fp,nullptr,error);
366 
367         if (i % nclusterprocs) {
368           iproc = me + (i % nclusterprocs);
369           MPI_Send(&n,1,MPI_INT,iproc,0,world);
370           MPI_Recv(&tmp,0,MPI_INT,iproc,0,world,MPI_STATUS_IGNORE);
371           MPI_Rsend(buf,n,MPI_DOUBLE,iproc,0,world);
372         }
373 
374       } else if (i % nclusterprocs == me - fileproc) {
375         MPI_Recv(&n,1,MPI_INT,fileproc,0,world,MPI_STATUS_IGNORE);
376         if (n > maxbuf) {
377           maxbuf = n;
378           memory->destroy(buf);
379           memory->create(buf,maxbuf,"read_restart:buf");
380         }
381         MPI_Irecv(buf,n,MPI_DOUBLE,fileproc,0,world,&request);
382         MPI_Send(&tmp,0,MPI_INT,fileproc,0,world);
383         MPI_Wait(&request,MPI_STATUS_IGNORE);
384       }
385 
386       if (i % nclusterprocs == me - fileproc) {
387         m = 0;
388         while (m < n) m += avec->unpack_restart(&buf[m]);
389       }
390     }
391 
392     if (filereader && fp != nullptr) {
393       fclose(fp);
394       fp = nullptr;
395     }
396     MPI_Comm_free(&clustercomm);
397   }
398 
399   // clean-up memory
400 
401   delete[] file;
402   memory->destroy(buf);
403 
404   // for multiproc or MPI-IO files:
405   // perform irregular comm to migrate atoms to correct procs
406 
407   if (multiproc || mpiioflag) {
408 
409     // if remapflag set, remap all atoms I read back to box before migrating
410 
411     if (remapflag) {
412       double **x = atom->x;
413       imageint *image = atom->image;
414       int nlocal = atom->nlocal;
415 
416       for (int i = 0; i < nlocal; i++)
417         domain->remap(x[i],image[i]);
418     }
419 
420     // create a temporary fix to hold and migrate extra atom info
421     // necessary b/c irregular will migrate atoms
422 
423     if (nextra)
424       modify->add_fix(fmt::format("_read_restart all READ_RESTART {} {}",
425                                   nextra,modify->nfix_restart_peratom));
426 
427     // move atoms to new processors via irregular()
428     // turn sorting on in migrate_atoms() to avoid non-reproducible restarts
429     // in case read by different proc than wrote restart file
430     // first do map_init() since irregular->migrate_atoms() will do map_clear()
431 
432     if (atom->map_style != Atom::MAP_NONE) {
433       atom->map_init();
434       atom->map_set();
435     }
436     if (domain->triclinic) domain->x2lamda(atom->nlocal);
437     Irregular *irregular = new Irregular(lmp);
438     irregular->migrate_atoms(1);
439     delete irregular;
440     if (domain->triclinic) domain->lamda2x(atom->nlocal);
441 
442     // put extra atom info held by fix back into atom->extra
443     // destroy temporary fix
444 
445     if (nextra) {
446       memory->destroy(atom->extra);
447       memory->create(atom->extra,atom->nmax,nextra,"atom:extra");
448       int ifix = modify->find_fix("_read_restart");
449       FixReadRestart *fix = (FixReadRestart *) modify->fix[ifix];
450       int *count = fix->count;
451       double **extra = fix->extra;
452       double **atom_extra = atom->extra;
453       int nlocal = atom->nlocal;
454       for (int i = 0; i < nlocal; i++)
455         for (int j = 0; j < count[i]; j++)
456           atom_extra[i][j] = extra[i][j];
457       modify->delete_fix("_read_restart");
458     }
459   }
460 
461   // check that all atoms were assigned to procs
462 
463   bigint natoms;
464   bigint nblocal = atom->nlocal;
465   MPI_Allreduce(&nblocal,&natoms,1,MPI_LMP_BIGINT,MPI_SUM,world);
466 
467   if (me == 0)
468     utils::logmesg(lmp,"  {} atoms\n",natoms);
469 
470   if (natoms != atom->natoms)
471     error->all(FLERR,"Did not assign all restart atoms correctly");
472 
473   if ((atom->molecular == Atom::TEMPLATE) && (me == 0)) {
474     std::string mesg;
475 
476     if (atom->nbonds)
477       mesg += fmt::format("  {} template bonds\n",atom->nbonds);
478     if (atom->nangles)
479       mesg += fmt::format("  {} template angles\n",atom->nangles);
480     if (atom->ndihedrals)
481       mesg += fmt::format("  {} template dihedrals\n",atom->ndihedrals);
482     if (atom->nimpropers)
483       mesg += fmt::format("  {} template impropers\n",atom->nimpropers);
484 
485     utils::logmesg(lmp,mesg);
486   }
487 
488   if ((atom->molecular == Atom::MOLECULAR) && (me == 0)) {
489     std::string mesg;
490     if (atom->nbonds)
491       mesg += fmt::format("  {} bonds\n",atom->nbonds);
492     if (atom->nangles)
493       mesg += fmt::format("  {} angles\n",atom->nangles);
494     if (atom->ndihedrals)
495       mesg += fmt::format("  {} dihedrals\n",atom->ndihedrals);
496     if (atom->nimpropers)
497       mesg += fmt::format("  {} impropers\n",atom->nimpropers);
498 
499     utils::logmesg(lmp,mesg);
500   }
501 
502   // check that atom IDs are valid
503 
504   atom->tag_check();
505 
506   // create global mapping of atoms
507 
508   if (atom->map_style != Atom::MAP_NONE) {
509     atom->map_init();
510     atom->map_set();
511   }
512 
513   // create special bond lists for molecular systems
514 
515   if (atom->molecular == Atom::MOLECULAR) {
516     Special special(lmp);
517     special.build();
518   }
519 
520   // total time
521 
522   MPI_Barrier(world);
523 
524   if (comm->me == 0)
525     utils::logmesg(lmp,"  read_restart CPU = {:.3f} seconds\n",MPI_Wtime()-time1);
526 
527   delete mpiio;
528 }
529 
530 /* ----------------------------------------------------------------------
531    inpfile contains a "*"
532    search for all files which match the inpfile pattern
533    replace "*" with latest timestep value to create outfile name
534    search dir referenced by initial pathname of file
535    if inpfile also contains "%", use "base" when searching directory
536    only called by proc 0
537 ------------------------------------------------------------------------- */
538 
file_search(const std::string & inpfile)539 std::string ReadRestart::file_search(const std::string &inpfile)
540 {
541   // separate inpfile into dir + filename
542 
543   auto dirname = utils::path_dirname(inpfile);
544   auto filename = utils::path_basename(inpfile);
545 
546   // if filename contains "%" replace "%" with "base"
547 
548   auto pattern = filename;
549   auto loc = pattern.find('%');
550   if (loc != std::string::npos) pattern.replace(loc,1,"base");
551 
552   // scan all files in directory, searching for files that match regexp pattern
553   // maxnum = largest integer that matches "*"
554 
555   bigint maxnum = -1;
556   loc = pattern.find('*');
557   if (loc != std::string::npos) {
558     // convert pattern to equivalent regexp
559     pattern.replace(loc,1,"\\d+");
560     struct dirent *ep;
561     DIR *dp = opendir(dirname.c_str());
562     if (dp == nullptr)
563       error->one(FLERR,"Cannot open directory {} to search for restart file: {}",
564                  dirname, utils::getsyserror());
565 
566     while ((ep = readdir(dp))) {
567       std::string candidate(ep->d_name);
568       if (utils::strmatch(candidate,pattern)) {
569         bigint num = ATOBIGINT(utils::strfind(candidate.substr(loc),"\\d+").c_str());
570         if (num > maxnum) maxnum = num;
571       }
572     }
573     closedir(dp);
574     if (maxnum < 0) error->one(FLERR,"Found no restart file matching pattern");
575     filename.replace(filename.find('*'),1,std::to_string(maxnum));
576   }
577   return utils::path_join(dirname,filename);
578 }
579 
580 /* ----------------------------------------------------------------------
581    read header of restart file
582 ------------------------------------------------------------------------- */
583 
header()584 void ReadRestart::header()
585 {
586   int xperiodic(-1),yperiodic(-1),zperiodic(-1);
587 
588   // read flags and fields until flag = -1
589 
590   int flag = read_int();
591   while (flag >= 0) {
592 
593     // check restart file version, warn if different
594 
595     if (flag == VERSION) {
596       char *version = read_string();
597       if (me == 0)
598         utils::logmesg(lmp,"  restart file = {}, LAMMPS = {}\n",
599                        version,lmp->version);
600       delete[] version;
601 
602       // we have no forward compatibility, thus exit with error
603 
604       if (revision > FORMAT_REVISION)
605         error->all(FLERR,"Restart file format revision incompatible "
606                    "with current LAMMPS version");
607 
608       // warn when attempting to read older format revision
609 
610       if ((me == 0) && (revision < FORMAT_REVISION))
611         error->warning(FLERR,"Old restart file format revision. "
612                        "Switching to compatibility mode.");
613 
614     // check lmptype.h sizes, error if different
615 
616     } else if (flag == SMALLINT) {
617       int size = read_int();
618       if (size != sizeof(smallint))
619         error->all(FLERR,"Smallint setting in lmptype.h is not compatible");
620     } else if (flag == IMAGEINT) {
621       int size = read_int();
622       if (size != sizeof(imageint))
623         error->all(FLERR,"Imageint setting in lmptype.h is not compatible");
624     } else if (flag == TAGINT) {
625       int size = read_int();
626       if (size != sizeof(tagint))
627         error->all(FLERR,"Tagint setting in lmptype.h is not compatible");
628     } else if (flag == BIGINT) {
629       int size = read_int();
630       if (size != sizeof(bigint))
631         error->all(FLERR,"Bigint setting in lmptype.h is not compatible");
632 
633     // reset unit_style only if different
634     // so that timestep,neighbor-skin are not changed
635 
636     } else if (flag == UNITS) {
637       char *style = read_string();
638       if (strcmp(style,update->unit_style) != 0) update->set_units(style);
639       delete[] style;
640 
641     } else if (flag == NTIMESTEP) {
642       update->ntimestep = read_bigint();
643 
644     // set dimension from restart file
645 
646     } else if (flag == DIMENSION) {
647       int dimension = read_int();
648       domain->dimension = dimension;
649       if (domain->dimension == 2 && domain->zperiodic == 0)
650         error->all(FLERR,
651                    "Cannot run 2d simulation with non-periodic Z dimension");
652 
653     // read nprocs from restart file, warn if different
654 
655     } else if (flag == NPROCS) {
656       nprocs_file = read_int();
657       if (nprocs_file != comm->nprocs && me == 0)
658         error->warning(FLERR,"Restart file used different # of processors: "
659                        "{} vs. {}",nprocs_file,comm->nprocs);
660 
661     // don't set procgrid, warn if different
662 
663     } else if (flag == PROCGRID) {
664       int procgrid[3];
665       read_int();
666       read_int_vec(3,procgrid);
667       flag = 0;
668       if (comm->user_procgrid[0] != 0 &&
669           procgrid[0] != comm->user_procgrid[0]) flag = 1;
670       if (comm->user_procgrid[1] != 0 &&
671           procgrid[1] != comm->user_procgrid[1]) flag = 1;
672       if (comm->user_procgrid[2] != 0 &&
673           procgrid[2] != comm->user_procgrid[2]) flag = 1;
674       if (flag && me == 0)
675         error->warning(FLERR,"Restart file used different 3d processor grid");
676 
677     // don't set newton_pair, leave input script value unchanged
678     // set newton_bond from restart file
679     // warn if different and input script settings are not default
680 
681     } else if (flag == NEWTON_PAIR) {
682       int newton_pair_file = read_int();
683       if (force->newton_pair != 1) {
684         if (newton_pair_file != force->newton_pair && me == 0)
685           error->warning(FLERR,
686                          "Restart file used different newton pair setting, "
687                          "using input script value");
688       }
689     } else if (flag == NEWTON_BOND) {
690       int newton_bond_file = read_int();
691       if (force->newton_bond != 1) {
692         if (newton_bond_file != force->newton_bond && me == 0)
693           error->warning(FLERR,
694                          "Restart file used different newton bond setting, "
695                          "using restart file value");
696       }
697       force->newton_bond = newton_bond_file;
698       if (force->newton_pair || force->newton_bond) force->newton = 1;
699       else force->newton = 0;
700 
701     // set boundary settings from restart file
702     // warn if different and input script settings are not default
703 
704     } else if (flag == XPERIODIC) {
705       xperiodic = read_int();
706     } else if (flag == YPERIODIC) {
707       yperiodic = read_int();
708     } else if (flag == ZPERIODIC) {
709       zperiodic = read_int();
710     } else if (flag == BOUNDARY) {
711       int boundary[3][2];
712       read_int();
713       read_int_vec(6,&boundary[0][0]);
714 
715       if (domain->boundary[0][0] || domain->boundary[0][1] ||
716           domain->boundary[1][0] || domain->boundary[1][1] ||
717           domain->boundary[2][0] || domain->boundary[2][1]) {
718         if (boundary[0][0] != domain->boundary[0][0] ||
719             boundary[0][1] != domain->boundary[0][1] ||
720             boundary[1][0] != domain->boundary[1][0] ||
721             boundary[1][1] != domain->boundary[1][1] ||
722             boundary[2][0] != domain->boundary[2][0] ||
723             boundary[2][1] != domain->boundary[2][1]) {
724           if (me == 0)
725             error->warning(FLERR,
726                            "Restart file used different boundary settings, "
727                            "using restart file values");
728         }
729       }
730 
731       domain->boundary[0][0] = boundary[0][0];
732       domain->boundary[0][1] = boundary[0][1];
733       domain->boundary[1][0] = boundary[1][0];
734       domain->boundary[1][1] = boundary[1][1];
735       domain->boundary[2][0] = boundary[2][0];
736       domain->boundary[2][1] = boundary[2][1];
737 
738       if (xperiodic < 0 || yperiodic < 0 || zperiodic < 0)
739         error->all(FLERR,"Illegal or unset periodicity in restart");
740 
741       domain->periodicity[0] = domain->xperiodic = xperiodic;
742       domain->periodicity[1] = domain->yperiodic = yperiodic;
743       domain->periodicity[2] = domain->zperiodic = zperiodic;
744 
745       domain->nonperiodic = 0;
746       if (xperiodic == 0 || yperiodic == 0 || zperiodic == 0) {
747         domain->nonperiodic = 1;
748         if (boundary[0][0] >= 2 || boundary[0][1] >= 2 ||
749             boundary[1][0] >= 2 || boundary[1][1] >= 2 ||
750             boundary[2][0] >= 2 || boundary[2][1] >= 2)
751           domain->nonperiodic = 2;
752       }
753 
754     } else if (flag == BOUNDMIN) {
755       double minbound[6];
756       read_int();
757       read_double_vec(6,minbound);
758       domain->minxlo = minbound[0]; domain->minxhi = minbound[1];
759       domain->minylo = minbound[2]; domain->minyhi = minbound[3];
760       domain->minzlo = minbound[4]; domain->minzhi = minbound[5];
761 
762     // create new AtomVec class using any stored args
763 
764     } else if (flag == ATOM_STYLE) {
765       char *style = read_string();
766       int nargcopy = read_int();
767       char **argcopy = new char*[nargcopy];
768       for (int i = 0; i < nargcopy; i++)
769         argcopy[i] = read_string();
770       atom->create_avec(style,nargcopy,argcopy,1);
771       if (comm->me ==0)
772         utils::logmesg(lmp,"  restoring atom style {} from restart\n",style);
773       for (int i = 0; i < nargcopy; i++) delete[] argcopy[i];
774       delete[] argcopy;
775       delete[] style;
776 
777     } else if (flag == NATOMS) {
778       atom->natoms = read_bigint();
779     } else if (flag == NTYPES) {
780       atom->ntypes = read_int();
781     } else if (flag == NBONDS) {
782       atom->nbonds = read_bigint();
783     } else if (flag == NBONDTYPES) {
784       atom->nbondtypes = read_int();
785     } else if (flag == BOND_PER_ATOM) {
786       atom->bond_per_atom = read_int();
787     } else if (flag == NANGLES) {
788       atom->nangles = read_bigint();
789     } else if (flag == NANGLETYPES) {
790       atom->nangletypes = read_int();
791     } else if (flag == ANGLE_PER_ATOM) {
792       atom->angle_per_atom = read_int();
793     } else if (flag == NDIHEDRALS) {
794       atom->ndihedrals = read_bigint();
795     } else if (flag == NDIHEDRALTYPES) {
796       atom->ndihedraltypes = read_int();
797     } else if (flag == DIHEDRAL_PER_ATOM) {
798       atom->dihedral_per_atom = read_int();
799     } else if (flag == NIMPROPERS) {
800       atom->nimpropers = read_bigint();
801     } else if (flag == NIMPROPERTYPES) {
802       atom->nimpropertypes = read_int();
803     } else if (flag == IMPROPER_PER_ATOM) {
804       atom->improper_per_atom = read_int();
805 
806     } else if (flag == TRICLINIC) {
807       domain->triclinic = read_int();
808     } else if (flag == BOXLO) {
809       read_int();
810       read_double_vec(3,domain->boxlo);
811     } else if (flag == BOXHI) {
812       read_int();
813       read_double_vec(3,domain->boxhi);
814     } else if (flag == XY) {
815       domain->xy = read_double();
816     } else if (flag == XZ) {
817       domain->xz = read_double();
818     } else if (flag == YZ) {
819       domain->yz = read_double();
820 
821     } else if (flag == SPECIAL_LJ) {
822       read_int();
823       read_double_vec(3,&force->special_lj[1]);
824     } else if (flag == SPECIAL_COUL) {
825       read_int();
826       read_double_vec(3,&force->special_coul[1]);
827 
828     } else if (flag == TIMESTEP) {
829       update->dt = read_double();
830 
831     } else if (flag == ATOM_ID) {
832       atom->tag_enable = read_int();
833     } else if (flag == ATOM_MAP_STYLE) {
834       atom->map_style = read_int();
835     } else if (flag == ATOM_MAP_USER) {
836       atom->map_user  = read_int();
837     } else if (flag == ATOM_SORTFREQ) {
838       atom->sortfreq = read_int();
839     } else if (flag == ATOM_SORTBIN) {
840       atom->userbinsize = read_double();
841 
842     } else if (flag == COMM_MODE) {
843       comm->mode = read_int();
844     } else if (flag == COMM_CUTOFF) {
845       comm->cutghostuser = read_double();
846     } else if (flag == COMM_VEL) {
847       comm->ghost_velocity = read_int();
848 
849     } else if (flag == EXTRA_BOND_PER_ATOM) {
850       atom->extra_bond_per_atom = read_int();
851     } else if (flag == EXTRA_ANGLE_PER_ATOM) {
852       atom->extra_angle_per_atom = read_int();
853     } else if (flag == EXTRA_DIHEDRAL_PER_ATOM) {
854       atom->extra_dihedral_per_atom = read_int();
855     } else if (flag == EXTRA_IMPROPER_PER_ATOM) {
856       atom->extra_improper_per_atom = read_int();
857     } else if (flag == ATOM_MAXSPECIAL) {
858       atom->maxspecial = read_int();
859     } else if (flag == NELLIPSOIDS) {
860       atom->nellipsoids = read_bigint();
861     } else if (flag == NLINES) {
862       atom->nlines = read_bigint();
863     } else if (flag == NTRIS) {
864       atom->ntris = read_bigint();
865     } else if (flag == NBODIES) {
866       atom->nbodies = read_bigint();
867 
868       // for backward compatibility
869     } else if (flag == EXTRA_SPECIAL_PER_ATOM) {
870       force->special_extra = read_int();
871 
872     } else error->all(FLERR,"Invalid flag in header section of restart file");
873 
874     flag = read_int();
875   }
876 }
877 
878 /* ---------------------------------------------------------------------- */
879 
type_arrays()880 void ReadRestart::type_arrays()
881 {
882   int flag = read_int();
883   while (flag >= 0) {
884 
885     if (flag == MASS) {
886       read_int();
887       double *mass = new double[atom->ntypes+1];
888       read_double_vec(atom->ntypes,&mass[1]);
889       atom->set_mass(mass);
890       delete[] mass;
891 
892     } else error->all(FLERR,
893                       "Invalid flag in type arrays section of restart file");
894 
895     flag = read_int();
896   }
897 }
898 
899 /* ---------------------------------------------------------------------- */
900 
force_fields()901 void ReadRestart::force_fields()
902 {
903   char *style;
904 
905   int flag = read_int();
906   while (flag >= 0) {
907 
908     if (flag == PAIR) {
909       style = read_string();
910       force->create_pair(style,1);
911       delete[] style;
912       if (comm->me ==0)
913         utils::logmesg(lmp,"  restoring pair style {} from restart\n",
914                        force->pair_style);
915       force->pair->read_restart(fp);
916 
917     } else if (flag == NO_PAIR) {
918       style = read_string();
919       if (comm->me ==0)
920         utils::logmesg(lmp,"  pair style {} stores no restart info\n", style);
921       force->create_pair("none",0);
922       force->pair_restart = style;
923 
924     } else if (flag == BOND) {
925       style = read_string();
926       force->create_bond(style,1);
927       delete[] style;
928       if (comm->me ==0)
929         utils::logmesg(lmp,"  restoring bond style {} from restart\n",
930                        force->bond_style);
931       force->bond->read_restart(fp);
932 
933     } else if (flag == ANGLE) {
934       style = read_string();
935       force->create_angle(style,1);
936       delete[] style;
937       if (comm->me ==0)
938         utils::logmesg(lmp,"  restoring angle style {} from restart\n",
939                        force->angle_style);
940       force->angle->read_restart(fp);
941 
942     } else if (flag == DIHEDRAL) {
943       style = read_string();
944       force->create_dihedral(style,1);
945       delete[] style;
946       if (comm->me ==0)
947         utils::logmesg(lmp,"  restoring dihedral style {} from restart\n",
948                        force->dihedral_style);
949       force->dihedral->read_restart(fp);
950 
951     } else if (flag == IMPROPER) {
952       style = read_string();
953       force->create_improper(style,1);
954       delete[] style;
955       if (comm->me ==0)
956         utils::logmesg(lmp,"  restoring improper style {} from restart\n",
957                        force->improper_style);
958       force->improper->read_restart(fp);
959 
960     } else error->all(FLERR,
961                       "Invalid flag in force field section of restart file");
962 
963     flag = read_int();
964   }
965 }
966 
967 /* ---------------------------------------------------------------------- */
968 
file_layout()969 void ReadRestart::file_layout()
970 {
971   int flag = read_int();
972   while (flag >= 0) {
973 
974     if (flag == MULTIPROC) {
975       multiproc_file = read_int();
976       if (multiproc == 0 && multiproc_file)
977         error->all(FLERR,"Restart file is not a multi-proc file");
978       if (multiproc && multiproc_file == 0)
979         error->all(FLERR,"Restart file is a multi-proc file");
980 
981     } else if (flag == MPIIO) {
982       int mpiioflag_file = read_int();
983       if (mpiioflag == 0 && mpiioflag_file)
984         error->all(FLERR,"Restart file is a MPI-IO file");
985       if (mpiioflag && mpiioflag_file == 0)
986         error->all(FLERR,"Restart file is not a MPI-IO file");
987 
988       if (mpiioflag) {
989         bigint *nproc_chunk_offsets;
990         memory->create(nproc_chunk_offsets,nprocs,
991                        "write_restart:nproc_chunk_offsets");
992         bigint *nproc_chunk_sizes;
993         memory->create(nproc_chunk_sizes,nprocs,
994                        "write_restart:nproc_chunk_sizes");
995 
996         // on rank 0 read in the chunk sizes that were written out
997         // then consolidate them and compute offsets relative to the
998         // end of the header info to fit the current partition size
999         // if the number of ranks that did the writing is different
1000 
1001         if (me == 0) {
1002           int ndx;
1003           int *all_written_send_sizes;
1004           memory->create(all_written_send_sizes,nprocs_file,
1005                          "write_restart:all_written_send_sizes");
1006           int *nproc_chunk_number;
1007           memory->create(nproc_chunk_number,nprocs,
1008                          "write_restart:nproc_chunk_number");
1009 
1010           utils::sfread(FLERR,all_written_send_sizes,sizeof(int),nprocs_file,fp,nullptr,error);
1011 
1012           if ((nprocs != nprocs_file) && !(atom->nextra_store)) {
1013             // nprocs differ, but atom sizes are fixed length, yeah!
1014             atom->nlocal = 1; // temporarily claim there is one atom...
1015             int perAtomSize = atom->avec->size_restart(); // ...so we can get its size
1016             atom->nlocal = 0; // restore nlocal to zero atoms
1017 
1018             bigint total_size = 0;
1019             for (int i = 0; i < nprocs_file; ++i) {
1020               total_size += all_written_send_sizes[i];
1021             }
1022             bigint total_ct = total_size / perAtomSize;
1023 
1024             bigint base_ct = total_ct / nprocs;
1025             bigint leftover_ct = total_ct  - (base_ct * nprocs);
1026             bigint current_ByteOffset = 0;
1027             base_ct += 1;
1028             bigint base_ByteOffset = base_ct * (perAtomSize * sizeof(double));
1029             for (ndx = 0; ndx < leftover_ct; ++ndx) {
1030               nproc_chunk_offsets[ndx] = current_ByteOffset;
1031               nproc_chunk_sizes[ndx] = base_ct * perAtomSize;
1032               current_ByteOffset += base_ByteOffset;
1033             }
1034             base_ct -= 1;
1035             base_ByteOffset -= (perAtomSize * sizeof(double));
1036             for (; ndx < nprocs; ++ndx) {
1037               nproc_chunk_offsets[ndx] = current_ByteOffset;
1038               nproc_chunk_sizes[ndx] = base_ct * perAtomSize;
1039               current_ByteOffset += base_ByteOffset;
1040             }
1041           } else { // we have to read in based on how it was written
1042             int init_chunk_number = nprocs_file/nprocs;
1043             int num_extra_chunks = nprocs_file - (nprocs*init_chunk_number);
1044 
1045             for (int i = 0; i < nprocs; i++) {
1046               if (i < num_extra_chunks)
1047                 nproc_chunk_number[i] = init_chunk_number+1;
1048               else
1049                 nproc_chunk_number[i] = init_chunk_number;
1050             }
1051 
1052             int all_written_send_sizes_index = 0;
1053             bigint current_offset = 0;
1054             for (int i=0;i<nprocs;i++) {
1055               nproc_chunk_offsets[i] = current_offset;
1056               nproc_chunk_sizes[i] = 0;
1057               for (int j=0;j<nproc_chunk_number[i];j++) {
1058                 nproc_chunk_sizes[i] +=
1059                   all_written_send_sizes[all_written_send_sizes_index];
1060                 current_offset +=
1061                   (all_written_send_sizes[all_written_send_sizes_index] *
1062                    sizeof(double));
1063                 all_written_send_sizes_index++;
1064               }
1065 
1066             }
1067           }
1068           memory->destroy(all_written_send_sizes);
1069           memory->destroy(nproc_chunk_number);
1070         }
1071 
1072         // scatter chunk sizes and offsets to all procs
1073 
1074         MPI_Scatter(nproc_chunk_sizes, 1, MPI_LMP_BIGINT,
1075                     &assignedChunkSize , 1, MPI_LMP_BIGINT, 0,world);
1076         MPI_Scatter(nproc_chunk_offsets, 1, MPI_LMP_BIGINT,
1077                     &assignedChunkOffset , 1, MPI_LMP_BIGINT, 0,world);
1078 
1079         memory->destroy(nproc_chunk_sizes);
1080         memory->destroy(nproc_chunk_offsets);
1081       }
1082     }
1083 
1084     flag = read_int();
1085   }
1086 
1087   // if MPI-IO file, broadcast the end of the header offste
1088   // this allows all ranks to compute offset to their data
1089 
1090   if (mpiioflag) {
1091     if (me == 0) headerOffset = ftell(fp);
1092     MPI_Bcast(&headerOffset,1,MPI_LMP_BIGINT,0,world);
1093   }
1094 }
1095 
1096 // ----------------------------------------------------------------------
1097 // ----------------------------------------------------------------------
1098 // low-level fread methods
1099 // ----------------------------------------------------------------------
1100 // ----------------------------------------------------------------------
1101 
1102 /* ----------------------------------------------------------------------
1103 ------------------------------------------------------------------------- */
1104 
magic_string()1105 void ReadRestart::magic_string()
1106 {
1107   int n = strlen(MAGIC_STRING) + 1;
1108   char *str = new char[n];
1109 
1110   int count;
1111   if (me == 0) count = fread(str,sizeof(char),n,fp);
1112   MPI_Bcast(&count,1,MPI_INT,0,world);
1113   if (count < n)
1114     error->all(FLERR,"Invalid LAMMPS restart file");
1115   MPI_Bcast(str,n,MPI_CHAR,0,world);
1116   if (strcmp(str,MAGIC_STRING) != 0)
1117     error->all(FLERR,"Invalid LAMMPS restart file");
1118   delete[] str;
1119 }
1120 
1121 /* ----------------------------------------------------------------------
1122 ------------------------------------------------------------------------- */
1123 
endian()1124 void ReadRestart::endian()
1125 {
1126   int endian = read_int();
1127   if (endian == ENDIAN) return;
1128   if (endian == ENDIANSWAP)
1129     error->all(FLERR,"Restart file byte ordering is swapped");
1130   else error->all(FLERR,"Restart file byte ordering is not recognized");
1131 }
1132 
1133 /* ----------------------------------------------------------------------
1134 ------------------------------------------------------------------------- */
1135 
format_revision()1136 void ReadRestart::format_revision()
1137 {
1138   revision = read_int();
1139 }
1140 
1141 /* ----------------------------------------------------------------------
1142 ------------------------------------------------------------------------- */
1143 
check_eof_magic()1144 void ReadRestart::check_eof_magic()
1145 {
1146   // no check for revision 0 restart files
1147   if (revision < 1) return;
1148 
1149   int n = strlen(MAGIC_STRING) + 1;
1150   char *str = new char[n];
1151 
1152   // read magic string at end of file and restore file pointer
1153 
1154   if (me == 0) {
1155     long curpos = ftell(fp);
1156     fseek(fp,(long)-n,SEEK_END);
1157     utils::sfread(FLERR,str,sizeof(char),n,fp,nullptr,error);
1158     fseek(fp,curpos,SEEK_SET);
1159   }
1160 
1161   MPI_Bcast(str,n,MPI_CHAR,0,world);
1162   if (strcmp(str,MAGIC_STRING) != 0)
1163     error->all(FLERR,"Incomplete or corrupted LAMMPS restart file");
1164 
1165   delete[] str;
1166 }
1167 
1168 /* ----------------------------------------------------------------------
1169    read an int from restart file and bcast it
1170 ------------------------------------------------------------------------- */
1171 
read_int()1172 int ReadRestart::read_int()
1173 {
1174   int value;
1175   if ((me == 0) && (fread(&value,sizeof(int),1,fp) < 1))
1176     value = -1;
1177   MPI_Bcast(&value,1,MPI_INT,0,world);
1178   return value;
1179 }
1180 
1181 /* ----------------------------------------------------------------------
1182    read a bigint from restart file and bcast it
1183 ------------------------------------------------------------------------- */
1184 
read_bigint()1185 bigint ReadRestart::read_bigint()
1186 {
1187   bigint value;
1188   if ((me == 0) && (fread(&value,sizeof(bigint),1,fp) < 1))
1189     value = -1;
1190   MPI_Bcast(&value,1,MPI_LMP_BIGINT,0,world);
1191   return value;
1192 }
1193 
1194 /* ----------------------------------------------------------------------
1195    read a double from restart file and bcast it
1196 ------------------------------------------------------------------------- */
1197 
read_double()1198 double ReadRestart::read_double()
1199 {
1200   double value;
1201   if ((me == 0) && (fread(&value,sizeof(double),1,fp) < 1))
1202     value = 0.0;
1203   MPI_Bcast(&value,1,MPI_DOUBLE,0,world);
1204   return value;
1205 }
1206 
1207 /* ----------------------------------------------------------------------
1208    read a char string (including nullptr) and bcast it
1209    str is allocated here, ptr is returned, caller must deallocate
1210 ------------------------------------------------------------------------- */
1211 
read_string()1212 char *ReadRestart::read_string()
1213 {
1214   int n = read_int();
1215   if (n < 0) error->all(FLERR,"Illegal size string or corrupt restart");
1216   char *value = new char[n];
1217   if (me == 0) utils::sfread(FLERR,value,sizeof(char),n,fp,nullptr,error);
1218   MPI_Bcast(value,n,MPI_CHAR,0,world);
1219   return value;
1220 }
1221 
1222 /* ----------------------------------------------------------------------
1223    read vector of N ints from restart file and bcast them
1224 ------------------------------------------------------------------------- */
1225 
read_int_vec(int n,int * vec)1226 void ReadRestart::read_int_vec(int n, int *vec)
1227 {
1228   if (n < 0) error->all(FLERR,"Illegal size integer vector read requested");
1229   if (me == 0) utils::sfread(FLERR,vec,sizeof(int),n,fp,nullptr,error);
1230   MPI_Bcast(vec,n,MPI_INT,0,world);
1231 }
1232 
1233 /* ----------------------------------------------------------------------
1234    read vector of N doubles from restart file and bcast them
1235 ------------------------------------------------------------------------- */
1236 
read_double_vec(int n,double * vec)1237 void ReadRestart::read_double_vec(int n, double *vec)
1238 {
1239   if (n < 0) error->all(FLERR,"Illegal size double vector read requested");
1240   if (me == 0) utils::sfread(FLERR,vec,sizeof(double),n,fp,nullptr,error);
1241   MPI_Bcast(vec,n,MPI_DOUBLE,0,world);
1242 }
1243