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