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 /* ----------------------------------------------------------------------
16 Contributing author: Timothy Sirk (ARL)
17 ------------------------------------------------------------------------- */
18
19 #include "read_dump.h"
20
21 #include "atom.h"
22 #include "atom_vec.h"
23 #include "comm.h"
24 #include "domain.h"
25 #include "error.h"
26 #include "irregular.h"
27 #include "memory.h"
28 #include "reader.h"
29 #include "style_reader.h" // IWYU pragma: keep
30 #include "update.h"
31
32 #include <cstring>
33
34 using namespace LAMMPS_NS;
35
36 #define CHUNK 16384
37
38 // also in reader_native.cpp
39
40 enum{ID,TYPE,X,Y,Z,VX,VY,VZ,Q,IX,IY,IZ,FX,FY,FZ};
41 enum{UNSET,NOSCALE_NOWRAP,NOSCALE_WRAP,SCALE_NOWRAP,SCALE_WRAP};
42 enum{NOADD,YESADD,KEEPADD};
43
44 /* ---------------------------------------------------------------------- */
45
ReadDump(LAMMPS * lmp)46 ReadDump::ReadDump(LAMMPS *lmp) : Command(lmp)
47 {
48 MPI_Comm_rank(world,&me);
49 MPI_Comm_size(world,&nprocs);
50
51 dimension = domain->dimension;
52 triclinic = domain->triclinic;
53
54 nfile = 0;
55 files = nullptr;
56
57 nnew = maxnew = 0;
58 nfield = 0;
59 fieldtype = nullptr;
60 fieldlabel = nullptr;
61 fields = nullptr;
62 buf = nullptr;
63
64 readerstyle = utils::strdup("native");
65
66 nreader = 0;
67 readers = nullptr;
68 nsnapatoms = nullptr;
69 clustercomm = MPI_COMM_NULL;
70 filereader = 0;
71 parallel = 0;
72 }
73
74 /* ---------------------------------------------------------------------- */
75
~ReadDump()76 ReadDump::~ReadDump()
77 {
78 for (int i = 0; i < nfile; i++) delete [] files[i];
79 delete [] files;
80 for (int i = 0; i < nfield; i++) delete [] fieldlabel[i];
81 delete [] fieldlabel;
82 delete [] fieldtype;
83 delete [] readerstyle;
84
85 memory->destroy(fields);
86 memory->destroy(buf);
87
88 for (int i = 0; i < nreader; i++) delete readers[i];
89 delete [] readers;
90 delete [] nsnapatoms;
91
92 MPI_Comm_free(&clustercomm);
93 }
94
95 /* ---------------------------------------------------------------------- */
96
command(int narg,char ** arg)97 void ReadDump::command(int narg, char **arg)
98 {
99 if (domain->box_exist == 0)
100 error->all(FLERR,"Read_dump command before simulation box is defined");
101
102 if (narg < 2) error->all(FLERR,"Illegal read_dump command");
103
104 store_files(1,&arg[0]);
105 bigint nstep = utils::bnumeric(FLERR,arg[1],false,lmp);
106
107 int nremain = narg - 2;
108 if (nremain) nremain = fields_and_keywords(nremain,&arg[narg-nremain]);
109 else nremain = fields_and_keywords(0,nullptr);
110 if (nremain) setup_reader(nremain,&arg[narg-nremain]);
111 else setup_reader(0,nullptr);
112
113 // find the snapshot and read/bcast/process header info
114
115 if (me == 0) utils::logmesg(lmp,"Scanning dump file ...\n");
116
117 bigint ntimestep = seek(nstep,1);
118 if (ntimestep < 0)
119 error->all(FLERR,"Dump file does not contain requested snapshot");
120 header(1);
121
122 // reset timestep to nstep
123
124 update->reset_timestep(nstep);
125
126 // counters
127
128 // read in the snapshot and reset system
129
130 if (me == 0) utils::logmesg(lmp,"Reading snapshot from dump file ...\n");
131
132 bigint natoms_prev = atom->natoms;
133 atoms();
134
135 if (filereader)
136 for (int i = 0; i < nreader; i++)
137 readers[i]->close_file();
138
139 // print out stats
140
141 bigint nsnap_all,npurge_all,nreplace_all,ntrim_all,nadd_all;
142
143 bigint tmp = 0;
144 if (filereader)
145 for (int i = 0; i < nreader; i++)
146 tmp += nsnapatoms[i];
147 MPI_Allreduce(&tmp,&nsnap_all,1,MPI_LMP_BIGINT,MPI_SUM,world);
148
149 tmp = npurge;
150 MPI_Allreduce(&tmp,&npurge_all,1,MPI_LMP_BIGINT,MPI_SUM,world);
151 tmp = nreplace;
152 MPI_Allreduce(&tmp,&nreplace_all,1,MPI_LMP_BIGINT,MPI_SUM,world);
153 tmp = ntrim;
154 MPI_Allreduce(&tmp,&ntrim_all,1,MPI_LMP_BIGINT,MPI_SUM,world);
155 tmp = nadd;
156 MPI_Allreduce(&tmp,&nadd_all,1,MPI_LMP_BIGINT,MPI_SUM,world);
157
158 domain->print_box(" ");
159
160 if (me == 0)
161 utils::logmesg(lmp," {} atoms before read\n"
162 " {} atoms in snapshot\n"
163 " {} atoms purged\n"
164 " {} atoms replaced\n"
165 " {} atoms trimmed\n"
166 " {} atoms added\n"
167 " {} atoms after read\n",natoms_prev,nsnap_all,
168 npurge_all,nreplace_all,ntrim_all,nadd_all,atom->natoms);
169 }
170
171 /* ---------------------------------------------------------------------- */
172
store_files(int nstr,char ** str)173 void ReadDump::store_files(int nstr, char **str)
174 {
175 nfile = nstr;
176 files = new char*[nfile];
177
178 // either all or none of files must have '%' wild-card
179
180 for (int i = 0; i < nfile; i++) {
181 files[i] = utils::strdup(str[i]);
182
183 if (i == 0) {
184 if (strchr(files[i],'%')) multiproc = 1;
185 else multiproc = 0;
186 } else {
187 if (multiproc && !strchr(files[i],'%'))
188 error->all(FLERR,"All read_dump files must be serial or parallel");
189 if (!multiproc && strchr(files[i],'%'))
190 error->all(FLERR,"All read_dump files must be serial or parallel");
191 }
192 }
193 }
194
195 /* ---------------------------------------------------------------------- */
196
setup_reader(int narg,char ** arg)197 void ReadDump::setup_reader(int narg, char **arg)
198 {
199 // setup serial or parallel file reading
200 // multiproc = 0: only one file to read from, only proc 0 is a reader
201 // multiproc_nfile >= nprocs: every proc reads one or more files
202 // multiproc_nfile < nprocs: multiproc_nfile readers, create clusters
203 // see read_dump.h for explanation of these variables
204
205 if (multiproc == 0) {
206 nreader = 1;
207 firstfile = -1;
208 MPI_Comm_dup(world,&clustercomm);
209 } else if (multiproc_nfile >= nprocs) {
210 firstfile = static_cast<int> ((bigint) me * multiproc_nfile/nprocs);
211 int lastfile = static_cast<int> ((bigint) (me+1) * multiproc_nfile/nprocs);
212 nreader = lastfile - firstfile;
213 MPI_Comm_split(world,me,0,&clustercomm);
214 } else if (multiproc_nfile < nprocs) {
215 nreader = 1;
216 int icluster = static_cast<int> ((bigint) me * multiproc_nfile/nprocs);
217 firstfile = icluster;
218 MPI_Comm_split(world,icluster,0,&clustercomm);
219 }
220
221 MPI_Comm_rank(clustercomm,&me_cluster);
222 MPI_Comm_size(clustercomm,&nprocs_cluster);
223 if (me_cluster == 0) filereader = 1;
224 else filereader = 0;
225
226 readers = new Reader*[nreader];
227 nsnapatoms = new bigint[nreader];
228 for (int i=0; i < nreader; ++i) {
229 readers[i] = nullptr;
230 nsnapatoms[i] = 0;
231 }
232
233 // create Nreader reader classes per reader
234 // match readerstyle to options in style_reader.h
235
236 if (0) {
237 return; // dummy line to enable else-if macro expansion
238
239 #define READER_CLASS
240 #define ReaderStyle(key,Class) \
241 } else if (strcmp(readerstyle,#key) == 0) { \
242 for (int i = 0; i < nreader; i++) { \
243 readers[i] = new Class(lmp); \
244 }
245 #include "style_reader.h" // IWYU pragma: keep
246 #undef READER_CLASS
247
248 // unrecognized style
249
250 } else error->all(FLERR,utils::check_packages_for_style("reader",readerstyle,lmp));
251
252 if (utils::strmatch(readerstyle,"^adios")) {
253 // everyone is a reader with adios
254 parallel = 1;
255 filereader = 1;
256 }
257
258 // pass any arguments to readers
259
260 if (narg > 0 && filereader)
261 for (int i = 0; i < nreader; i++)
262 readers[i]->settings(narg,arg);
263 }
264
265 /* ----------------------------------------------------------------------
266 seek Nrequest timestep in one or more dump files
267 if exact = 1, must find exactly Nrequest
268 if exact = 0, find first step >= Nrequest
269 return matching ntimestep or -1 if did not find a match
270 ------------------------------------------------------------------------- */
271
seek(bigint nrequest,int exact)272 bigint ReadDump::seek(bigint nrequest, int exact)
273 {
274 int ifile,eofflag;
275 bigint ntimestep;
276
277 // proc 0 finds the timestep in its first reader
278
279 if (me == 0 || parallel) {
280
281 // exit file loop when dump timestep >= nrequest
282 // or files exhausted
283
284 for (ifile = 0; ifile < nfile; ifile++) {
285 ntimestep = -1;
286 if (multiproc) {
287 std::string multiname = files[ifile];
288 multiname.replace(multiname.find('%'),1,"0");
289 readers[0]->open_file(multiname.c_str());
290 } else readers[0]->open_file(files[ifile]);
291
292 while (1) {
293 eofflag = readers[0]->read_time(ntimestep);
294 if (eofflag) break;
295 if (ntimestep >= nrequest) break;
296 readers[0]->skip();
297 }
298
299 if (ntimestep >= nrequest) break;
300 readers[0]->close_file();
301 }
302
303 currentfile = ifile;
304 if (ntimestep < nrequest) ntimestep = -1;
305 if (exact && ntimestep != nrequest) ntimestep = -1;
306 }
307
308 if (!parallel) {
309 // proc 0 broadcasts timestep and currentfile to all procs
310
311 MPI_Bcast(&ntimestep,1,MPI_LMP_BIGINT,0,world);
312 MPI_Bcast(¤tfile,1,MPI_INT,0,world);
313 }
314
315 // if ntimestep < 0:
316 // all filereader procs close all their files and return
317
318 if (ntimestep < 0) {
319 if (filereader)
320 for (int i = 0; i < nreader; i++)
321 readers[i]->close_file();
322 return ntimestep;
323 }
324
325 // for multiproc mode:
326 // all filereader procs search for same ntimestep in currentfile
327
328 if (multiproc && filereader) {
329 for (int i = 0; i < nreader; i++) {
330 if (me == 0 && i == 0) continue; // proc 0, reader 0 already found it
331 std::string multiname = files[currentfile];
332 multiname.replace(multiname.find('%'),1,fmt::format("{}",firstfile+i));
333 readers[i]->open_file(multiname.c_str());
334
335 bigint step;
336 while (1) {
337 eofflag = readers[i]->read_time(step);
338 if (eofflag) break;
339 if (step == ntimestep) break;
340 readers[i]->skip();
341 }
342
343 if (eofflag)
344 error->one(FLERR,"Read dump parallel files "
345 "do not all have same timestep");
346 }
347 }
348
349 return ntimestep;
350 }
351
352 /* ----------------------------------------------------------------------
353 find next matching snapshot in one or more dump files
354 Ncurrent = current timestep from last snapshot
355 Nlast = match no timestep bigger than Nlast
356 Nevery = only match timesteps that are a multiple of Nevery
357 Nskip = skip every this many timesteps
358 return matching ntimestep or -1 if did not find a match
359 ------------------------------------------------------------------------- */
360
next(bigint ncurrent,bigint nlast,int nevery,int nskip)361 bigint ReadDump::next(bigint ncurrent, bigint nlast, int nevery, int nskip)
362 {
363 int ifile,eofflag;
364 bigint ntimestep;
365
366 // proc 0 finds the timestep in its first reader
367
368 if (me == 0 || parallel) {
369
370 // exit file loop when dump timestep matches all criteria
371 // or files exhausted
372
373 int iskip = 0;
374
375 for (ifile = currentfile; ifile < nfile; ifile++) {
376 ntimestep = -1;
377 if (ifile != currentfile) {
378 if (multiproc) {
379 std::string multiname = files[ifile];
380 multiname.replace(multiname.find('%'),1,"0");
381 readers[0]->open_file(multiname.c_str());
382 } else readers[0]->open_file(files[ifile]);
383 }
384
385 while (1) {
386 eofflag = readers[0]->read_time(ntimestep);
387 if (eofflag) break;
388 if (ntimestep > nlast) break;
389 if (ntimestep <= ncurrent) {
390 readers[0]->skip();
391 continue;
392 }
393 if (iskip == nskip) iskip = 0;
394 iskip++;
395 if (nevery && ntimestep % nevery) readers[0]->skip();
396 else if (iskip < nskip) readers[0]->skip();
397 else break;
398 }
399
400 if (eofflag) readers[0]->close_file();
401 else break;
402 }
403
404 currentfile = ifile;
405 if (eofflag) ntimestep = -1;
406 if (ntimestep <= ncurrent) ntimestep = -1;
407 if (ntimestep > nlast) ntimestep = -1;
408 }
409
410 if (!parallel) {
411 // proc 0 broadcasts timestep and currentfile to all procs
412
413 MPI_Bcast(&ntimestep,1,MPI_LMP_BIGINT,0,world);
414 MPI_Bcast(¤tfile,1,MPI_INT,0,world);
415 }
416
417 // if ntimestep < 0:
418 // all filereader procs close all their files and return
419
420 if (ntimestep < 0) {
421 if (filereader)
422 for (int i = 0; i < nreader; i++)
423 readers[i]->close_file();
424 return ntimestep;
425 }
426
427 // for multiproc mode:
428 // all filereader procs search for same ntimestep in currentfile
429
430 if (multiproc && filereader) {
431 for (int i = 0; i < nreader; i++) {
432 if (me == 0 && i == 0) continue;
433 std::string multiname = files[currentfile];
434 multiname.replace(multiname.find('%'),1,fmt::format("{}",firstfile+i));
435 readers[i]->open_file(multiname.c_str());
436
437 bigint step;
438 while (1) {
439 eofflag = readers[i]->read_time(step);
440 if (eofflag) break;
441 if (step == ntimestep) break;
442 readers[i]->skip();
443 }
444
445 if (eofflag)
446 error->one(FLERR,"Read dump parallel files "
447 "do not all have same timestep");
448 }
449 }
450
451 return ntimestep;
452 }
453
454 /* ----------------------------------------------------------------------
455 read and broadcast and store snapshot header info
456 set nsnapatoms = # of atoms in snapshot
457 ------------------------------------------------------------------------- */
458
header(int fieldinfo)459 void ReadDump::header(int fieldinfo)
460 {
461 int boxinfo, triclinic_snap;
462 int fieldflag,xflag,yflag,zflag;
463
464 if (filereader) {
465 for (int i = 0; i < nreader; i++)
466 nsnapatoms[i] = readers[i]->read_header(box,boxinfo,triclinic_snap,fieldinfo,
467 nfield,fieldtype,fieldlabel,
468 scaleflag,wrapflag,fieldflag,
469 xflag,yflag,zflag);
470 }
471
472 if (!parallel) {
473 MPI_Bcast(nsnapatoms,nreader,MPI_LMP_BIGINT,0,clustercomm);
474 MPI_Bcast(&boxinfo,1,MPI_INT,0,clustercomm);
475 MPI_Bcast(&triclinic_snap,1,MPI_INT,0,clustercomm);
476 MPI_Bcast(&box[0][0],9,MPI_DOUBLE,0,clustercomm);
477 }
478
479 // local copy of snapshot box parameters
480 // used in xfield,yfield,zfield when converting dump atom to absolute coords
481
482 if (boxinfo) {
483 xlo = box[0][0];
484 xhi = box[0][1];
485 ylo = box[1][0];
486 yhi = box[1][1];
487 zlo = box[2][0];
488 zhi = box[2][1];
489
490 if (triclinic_snap) {
491 xy = box[0][2];
492 xz = box[1][2];
493 yz = box[2][2];
494 double xdelta = MIN(0.0,xy);
495 xdelta = MIN(xdelta,xz);
496 xdelta = MIN(xdelta,xy+xz);
497 xlo = xlo - xdelta;
498 xdelta = MAX(0.0,xy);
499 xdelta = MAX(xdelta,xz);
500 xdelta = MAX(xdelta,xy+xz);
501 xhi = xhi - xdelta;
502 ylo = ylo - MIN(0.0,yz);
503 yhi = yhi - MAX(0.0,yz);
504 }
505 xprd = xhi - xlo;
506 yprd = yhi - ylo;
507 zprd = zhi - zlo;
508 }
509
510 // done if not checking fields
511
512 if (!fieldinfo) return;
513
514 MPI_Bcast(&fieldflag,1,MPI_INT,0,clustercomm);
515 MPI_Bcast(&xflag,1,MPI_INT,0,clustercomm);
516 MPI_Bcast(&yflag,1,MPI_INT,0,clustercomm);
517 MPI_Bcast(&zflag,1,MPI_INT,0,clustercomm);
518
519 // error check on current vs new box and fields
520 // boxinfo == 0 means no box info in file
521
522 if (boxflag) {
523 if (!boxinfo)
524 error->all(FLERR,"No box information in dump, must use 'box no'");
525 else if ((triclinic_snap && !triclinic) ||
526 (!triclinic_snap && triclinic))
527 error->one(FLERR,"Read_dump triclinic status does not match simulation");
528 }
529
530 // error check on requested fields existing in dump file
531
532 if (fieldflag < 0)
533 error->one(FLERR,"Read_dump field not found in dump file");
534
535 // all explicitly requested x,y,z must have consistent scaling & wrapping
536
537 int value = MAX(xflag,yflag);
538 value = MAX(zflag,value);
539 if ((xflag != UNSET && xflag != value) ||
540 (yflag != UNSET && yflag != value) ||
541 (zflag != UNSET && zflag != value))
542 error->one(FLERR,
543 "Read_dump xyz fields do not have consistent scaling/wrapping");
544
545 // set scaled/wrapped based on xyz flags
546
547 value = UNSET;
548 if (xflag != UNSET) value = xflag;
549 if (yflag != UNSET) value = yflag;
550 if (zflag != UNSET) value = zflag;
551
552 if (value == UNSET) {
553 scaled = wrapped = 0;
554 } else if (value == NOSCALE_NOWRAP) {
555 scaled = wrapped = 0;
556 } else if (value == NOSCALE_WRAP) {
557 scaled = 0;
558 wrapped = 1;
559 } else if (value == SCALE_NOWRAP) {
560 scaled = 1;
561 wrapped = 0;
562 } else if (value == SCALE_WRAP) {
563 scaled = wrapped = 1;
564 }
565
566 // scaled, triclinic coords require all 3 x,y,z fields, to perform unscaling
567 // set yindex,zindex = column index of Y and Z fields in fields array
568 // needed for unscaling to absolute coords in xfield(), yfield(), zfield()
569
570 if (scaled && triclinic == 1) {
571 int flag = 0;
572 if (xflag == UNSET) flag = 1;
573 if (yflag == UNSET) flag = 1;
574 if (dimension == 3 && zflag == UNSET) flag = 1;
575 if (flag)
576 error->one(FLERR,"All read_dump x,y,z fields must be specified for "
577 "scaled, triclinic coords");
578
579 for (int i = 0; i < nfield; i++) {
580 if (fieldtype[i] == Y) yindex = i;
581 if (fieldtype[i] == Z) zindex = i;
582 }
583 }
584 }
585
586 /* ----------------------------------------------------------------------
587 read and process one snapshot of atoms
588 ------------------------------------------------------------------------- */
589
atoms()590 void ReadDump::atoms()
591 {
592 // initialize counters
593
594 npurge = nreplace = ntrim = nadd = 0;
595
596 // if purgeflag set, delete all current atoms
597
598 if (purgeflag) {
599 if (atom->map_style != Atom::MAP_NONE) atom->map_clear();
600 npurge = atom->nlocal;
601 atom->nlocal = atom->nghost = 0;
602 atom->natoms = 0;
603 }
604
605 // read all the snapshot atoms into fields
606 // each proc will own an arbitrary subset of atoms
607
608 read_atoms();
609
610 // migrate old owned atoms to new procs based on atom IDs
611 // not necessary if purged all old atoms or if only 1 proc
612
613 if (!purgeflag && nprocs > 1) migrate_old_atoms();
614
615 // migrate new snapshot atoms to same new procs based on atom IDs
616 // not necessary if purged all old atoms or if only 1 proc
617
618 if (!purgeflag && nprocs > 1) migrate_new_atoms();
619
620 // must build map if not a molecular system
621 // this will be needed to match new atoms to old atoms
622
623 int mapflag = 0;
624 if (atom->map_style == Atom::MAP_NONE) {
625 mapflag = 1;
626 atom->map_init();
627 atom->map_set();
628 }
629
630 // each proc now owns both old and new info for same subset of atoms
631 // update each local atom with new info
632
633 process_atoms();
634
635 // check that atom IDs are valid
636
637 atom->tag_check();
638
639 // delete atom map if created it above
640 // else reinitialize map for current atoms
641 // do this before migrating atoms to new procs via Irregular
642
643 if (mapflag) {
644 atom->map_delete();
645 atom->map_style = Atom::MAP_NONE;
646 } else {
647 atom->nghost = 0;
648 atom->map_init();
649 atom->map_set();
650 }
651
652 // overwrite simulation box with dump snapshot box if requested
653 // reallocate processors to box
654
655 if (boxflag) {
656 domain->boxlo[0] = xlo;
657 domain->boxhi[0] = xhi;
658 domain->boxlo[1] = ylo;
659 domain->boxhi[1] = yhi;
660 if (dimension == 3) {
661 domain->boxlo[2] = zlo;
662 domain->boxhi[2] = zhi;
663 }
664 if (triclinic) {
665 domain->xy = xy;
666 if (dimension == 3) {
667 domain->xz = xz;
668 domain->yz = yz;
669 }
670 }
671
672 domain->set_initial_box();
673 domain->set_global_box();
674 comm->set_proc_grid(0);
675 domain->set_local_box();
676 }
677
678 // migrate atoms to their new owing proc, based on atom coords
679
680 migrate_atoms_by_coords();
681 }
682
683 /* ----------------------------------------------------------------------
684 read all the snapshot atoms into fields
685 done in different ways for multiproc no/yes and # of procs < or >= nprocs
686 nnew = # of snapshot atoms this proc stores
687 ------------------------------------------------------------------------- */
688
read_atoms()689 void ReadDump::read_atoms()
690 {
691 int count,nread,nsend,nrecv,otherproc;
692 bigint nsnap,ntotal,ofirst,olast,rfirst,rlast,lo,hi;
693 MPI_Request request;
694 MPI_Status status;
695
696 // one reader per cluster of procs
697 // each reading proc reads one file and splits data across cluster
698 // cluster can be all procs or a subset
699
700 if (!parallel && (!multiproc || multiproc_nfile < nprocs)) {
701 nsnap = nsnapatoms[0];
702
703 if (filereader) {
704 if (!buf) memory->create(buf,CHUNK,nfield,"read_dump:buf");
705
706 otherproc = 0;
707 ofirst = (bigint) otherproc * nsnap/nprocs_cluster;
708 olast = (bigint) (otherproc+1) * nsnap/nprocs_cluster;
709 if (olast-ofirst > MAXSMALLINT)
710 error->one(FLERR,"Read dump snapshot is too large for a proc");
711 nnew = static_cast<int> (olast - ofirst);
712
713 if (nnew > maxnew || maxnew == 0) {
714 memory->destroy(fields);
715 maxnew = MAX(nnew,1); // avoid null pointer
716 memory->create(fields,maxnew,nfield,"read_dump:fields");
717 }
718
719 ntotal = 0;
720 while (ntotal < nsnap) {
721 nread = MIN(CHUNK,nsnap-ntotal);
722 readers[0]->read_atoms(nread,nfield,buf);
723 rfirst = ntotal;
724 rlast = ntotal + nread;
725
726 nsend = 0;
727 while (nsend < nread) {
728 lo = MAX(ofirst,rfirst);
729 hi = MIN(olast,rlast);
730 if (otherproc) // send to otherproc or copy to self
731 MPI_Send(&buf[nsend][0],(hi-lo)*nfield,MPI_DOUBLE,
732 otherproc,0,clustercomm);
733 else
734 memcpy(&fields[rfirst][0],&buf[nsend][0],
735 (hi-lo)*nfield*sizeof(double));
736 nsend += hi-lo;
737 if (hi == olast) {
738 otherproc++;
739 ofirst = (bigint) otherproc * nsnap/nprocs_cluster;
740 olast = (bigint) (otherproc+1) * nsnap/nprocs_cluster;
741 }
742 }
743
744 ntotal += nread;
745 }
746
747 } else {
748 ofirst = (bigint) me_cluster * nsnap/nprocs_cluster;
749 olast = (bigint) (me_cluster+1) * nsnap/nprocs_cluster;
750 if (olast-ofirst > MAXSMALLINT)
751 error->one(FLERR,"Read dump snapshot is too large for a proc");
752 nnew = static_cast<int> (olast - ofirst);
753 if (nnew > maxnew || maxnew == 0) {
754 memory->destroy(fields);
755 maxnew = MAX(nnew,1); // avoid null pointer
756 memory->create(fields,maxnew,nfield,"read_dump:fields");
757 }
758
759 nrecv = 0;
760 while (nrecv < nnew) {
761 MPI_Irecv(&fields[nrecv][0],(nnew-nrecv)*nfield,MPI_DOUBLE,0,0,
762 clustercomm,&request);
763 MPI_Wait(&request,&status);
764 MPI_Get_count(&status,MPI_DOUBLE,&count);
765 nrecv += count/nfield;
766 }
767 }
768
769 // every proc is a filereader, reads one or more files
770 // each proc keeps all data it reads, no communication required
771
772 } else if (multiproc_nfile >= nprocs || parallel) {
773 bigint sum = 0;
774 for (int i = 0; i < nreader; i++)
775 sum += nsnapatoms[i];
776 if (sum > MAXSMALLINT)
777 error->one(FLERR,"Read dump snapshot is too large for a proc");
778 nnew = static_cast<int> (sum);
779 if (nnew > maxnew || maxnew == 0) {
780 memory->destroy(fields);
781 maxnew = MAX(nnew,1); // avoid null pointer
782 memory->create(fields,maxnew,nfield,"read_dump:fields");
783 }
784
785 nnew = 0;
786 for (int i = 0; i < nreader; i++) {
787 nsnap = nsnapatoms[i];
788 ntotal = 0;
789 while (ntotal < nsnap) {
790 if (parallel) {
791 // read the whole thing at once
792 nread = nsnap-ntotal;
793 } else {
794 nread = MIN(CHUNK,nsnap-ntotal);
795 }
796 readers[i]->read_atoms(nread,nfield,&fields[nnew+ntotal]);
797 ntotal += nread;
798 }
799 nnew += nsnap;
800 }
801 }
802 }
803
804 /* ----------------------------------------------------------------------
805 update info for each old atom I own based on snapshot info
806 if in replace mode and atom ID matches current atom,
807 overwrite atom info with fields from dump file
808 if in add mode and atom ID does not match any old atom,
809 create new atom with dump file field values
810 ------------------------------------------------------------------------- */
811
process_atoms()812 void ReadDump::process_atoms()
813 {
814 int i,m,ifield,itype;
815 int xbox,ybox,zbox;
816 tagint mtag;
817 int *updateflag,*newflag;
818
819 // updateflag[i] = flag for old atoms, 1 if updated, else 0
820 // newflag[i] = flag for new atoms, 0 if used to update old atom, else 1
821
822 int nlocal = atom->nlocal;
823 memory->create(updateflag,nlocal,"read_dump:updateflag");
824 for (i = 0; i < nlocal; i++) updateflag[i] = 0;
825 memory->create(newflag,nnew,"read_dump:newflag");
826 for (i = 0; i < nnew; i++) newflag[i] = 1;
827
828 // loop over new atoms
829
830 double **x = atom->x;
831 double **v = atom->v;
832 double *q = atom->q;
833 double **f = atom->f;
834 tagint *tag = atom->tag;
835 imageint *image = atom->image;
836 tagint map_tag_max = atom->map_tag_max;
837
838 for (i = 0; i < nnew; i++) {
839
840 // check if new atom matches one I own
841 // setting m = -1 forces new atom not to match
842 // NOTE: atom ID in fields is stored as double, not as ubuf
843 // so can only cast it to tagint, thus cannot be full 64-bit ID
844
845 mtag = static_cast<tagint> (fields[i][0]);
846 if (mtag <= map_tag_max) m = atom->map(mtag);
847 else m = -1;
848 if (m < 0 || m >= nlocal) continue;
849
850 updateflag[m] = 1;
851 newflag[i] = 0;
852
853 if (replaceflag) {
854 nreplace++;
855
856 // current image flags
857
858 xbox = (image[m] & IMGMASK) - IMGMAX;
859 ybox = (image[m] >> IMGBITS & IMGMASK) - IMGMAX;
860 zbox = (image[m] >> IMG2BITS) - IMGMAX;
861
862 // overwrite atom attributes with field info
863 // start from field 1 since 0 = id, 1 will be skipped if type
864
865 for (ifield = 1; ifield < nfield; ifield++) {
866 switch (fieldtype[ifield]) {
867 case X:
868 x[m][0] = xfield(i,ifield);
869 break;
870 case Y:
871 x[m][1] = yfield(i,ifield);
872 break;
873 case Z:
874 x[m][2] = zfield(i,ifield);
875 break;
876 case VX:
877 v[m][0] = fields[i][ifield];
878 break;
879 case Q:
880 q[m] = fields[i][ifield];
881 break;
882 case VY:
883 v[m][1] = fields[i][ifield];
884 break;
885 case VZ:
886 v[m][2] = fields[i][ifield];
887 break;
888 case IX:
889 xbox = static_cast<int> (fields[i][ifield]);
890 break;
891 case IY:
892 ybox = static_cast<int> (fields[i][ifield]);
893 break;
894 case IZ:
895 zbox = static_cast<int> (fields[i][ifield]);
896 break;
897 case FX:
898 f[m][0] = fields[i][ifield];
899 break;
900 case FY:
901 f[m][1] = fields[i][ifield];
902 break;
903 case FZ:
904 f[m][2] = fields[i][ifield];
905 break;
906 }
907 }
908
909 // replace image flag in case changed by ix,iy,iz fields or unwrapping
910
911 if (!wrapped) xbox = ybox = zbox = 0;
912
913 image[m] = ((imageint) (xbox + IMGMAX) & IMGMASK) |
914 (((imageint) (ybox + IMGMAX) & IMGMASK) << IMGBITS) |
915 (((imageint) (zbox + IMGMAX) & IMGMASK) << IMG2BITS);
916 }
917 }
918
919 // if trimflag set, delete atoms not updated by snapshot atoms
920
921 if (trimflag) {
922 AtomVec *avec = atom->avec;
923
924 i = 0;
925 while (i < nlocal) {
926 if (!updateflag[i]) {
927 avec->copy(nlocal-1,i,1);
928 updateflag[i] = updateflag[nlocal-1];
929 nlocal--;
930 ntrim++;
931 } else i++;
932 }
933
934 atom->nlocal = nlocal;
935 bigint nblocal = atom->nlocal;
936 MPI_Allreduce(&nblocal,&atom->natoms,1,MPI_LMP_BIGINT,MPI_SUM,world);
937 }
938
939 // done if cannot add new atoms
940
941 if (addflag == NOADD) {
942 memory->destroy(updateflag);
943 memory->destroy(newflag);
944 return;
945 }
946
947 // ----------------------------------------------------
948 // create new atoms for dump file atoms with ID that matches no old atom
949 // ----------------------------------------------------
950
951 // first check that dump file snapshot has atom type field
952
953 int tflag = 0;
954 for (ifield = 0; ifield < nfield; ifield++)
955 if (fieldtype[ifield] == TYPE) tflag = 1;
956 if (!tflag)
957 error->all(FLERR,"Cannot add atoms if dump file does not store atom type");
958
959 int nlocal_previous = atom->nlocal;
960 double one[3];
961
962 for (i = 0; i < nnew; i++) {
963 if (!newflag[i]) continue;
964
965 // create type and coord fields from dump file
966 // coord = 0.0 unless corresponding dump file field was specified
967
968 itype = 0;
969 one[0] = one[1] = one[2] = 0.0;
970 for (ifield = 1; ifield < nfield; ifield++) {
971 switch (fieldtype[ifield]) {
972 case TYPE:
973 itype = static_cast<int> (fields[i][ifield]);
974 break;
975 case X:
976 one[0] = xfield(i,ifield);
977 break;
978 case Y:
979 one[1] = yfield(i,ifield);
980 break;
981 case Z:
982 one[2] = zfield(i,ifield);
983 break;
984 }
985 }
986
987 // create the atom on proc that owns it
988 // reset v,image ptrs in case they are reallocated
989
990 m = atom->nlocal;
991 atom->avec->create_atom(itype,one);
992 nadd++;
993
994 tag = atom->tag;
995 v = atom->v;
996 q = atom->q;
997 image = atom->image;
998
999 // set atom attributes from other dump file fields
1000
1001 xbox = ybox = zbox = 0;
1002
1003 for (ifield = 0; ifield < nfield; ifield++) {
1004 switch (fieldtype[ifield]) {
1005 case ID:
1006 if (addflag == KEEPADD)
1007 tag[m] = static_cast<tagint> (fields[i][ifield]);
1008 break;
1009 case VX:
1010 v[m][0] = fields[i][ifield];
1011 break;
1012 case VY:
1013 v[m][1] = fields[i][ifield];
1014 break;
1015 case VZ:
1016 v[m][2] = fields[i][ifield];
1017 break;
1018 case Q:
1019 q[m] = fields[i][ifield];
1020 break;
1021 case IX:
1022 xbox = static_cast<int> (fields[i][ifield]);
1023 break;
1024 case IY:
1025 ybox = static_cast<int> (fields[i][ifield]);
1026 break;
1027 case IZ:
1028 zbox = static_cast<int> (fields[i][ifield]);
1029 break;
1030 }
1031
1032 // reset image flag in case changed by ix,iy,iz fields
1033
1034 image[m] = ((imageint) (xbox + IMGMAX) & IMGMASK) |
1035 (((imageint) (ybox + IMGMAX) & IMGMASK) << IMGBITS) |
1036 (((imageint) (zbox + IMGMAX) & IMGMASK) << IMG2BITS);
1037 }
1038 }
1039
1040 // if addflag = YESADD or KEEPADD, update total atom count
1041
1042 if (addflag == YESADD || addflag == KEEPADD) {
1043 bigint nblocal = atom->nlocal;
1044 MPI_Allreduce(&nblocal,&atom->natoms,1,MPI_LMP_BIGINT,MPI_SUM,world);
1045 }
1046
1047 // if addflag = YESADD,
1048 // assign consistent IDs to new snapshot atoms across all procs
1049
1050 if (addflag == YESADD) {
1051 if (atom->natoms < 0 || atom->natoms >= MAXBIGINT)
1052 error->all(FLERR,"Too many total atoms");
1053 if (atom->tag_enable) atom->tag_extend();
1054 }
1055
1056 // init per-atom fix/compute/variable values for created atoms
1057
1058 atom->data_fix_compute_variable(nlocal_previous,atom->nlocal);
1059
1060 // free allocated vectors
1061
1062 memory->destroy(updateflag);
1063 memory->destroy(newflag);
1064 }
1065
1066 /* ----------------------------------------------------------------------
1067 migrate old atoms to new procs based on atom IDs
1068 use migrate_atoms() with explicit processor assignments
1069 ------------------------------------------------------------------------- */
1070
migrate_old_atoms()1071 void ReadDump::migrate_old_atoms()
1072 {
1073 tagint *tag = atom->tag;
1074 int nlocal = atom->nlocal;
1075
1076 int *procassign;
1077 memory->create(procassign,nlocal,"read_dump:procassign");
1078 for (int i = 0; i < nlocal; i++)
1079 procassign[i] = tag[i] % nprocs;
1080
1081 Irregular *irregular = new Irregular(lmp);
1082 irregular->migrate_atoms(1,1,procassign);
1083 delete irregular;
1084
1085 memory->destroy(procassign);
1086 }
1087
1088 /* ----------------------------------------------------------------------
1089 migrate new atoms to same new procs based on atom IDs
1090 ------------------------------------------------------------------------- */
1091
migrate_new_atoms()1092 void ReadDump::migrate_new_atoms()
1093 {
1094 tagint mtag;
1095 int *procassign;
1096 double **newfields;
1097
1098 memory->create(procassign,nnew,"read_dump:procassign");
1099 for (int i = 0; i < nnew; i++) {
1100 mtag = static_cast<tagint> (fields[i][0]);
1101 procassign[i] = mtag % nprocs;
1102 }
1103
1104 Irregular *irregular = new Irregular(lmp);
1105 int nrecv = irregular->create_data(nnew,procassign,1);
1106 int newmaxnew = MAX(nrecv,maxnew);
1107 newmaxnew = MAX(newmaxnew,1); // avoid null pointer
1108 memory->create(newfields,newmaxnew,nfield,"read_dump:newfields");
1109 irregular->exchange_data((char *) &fields[0][0],nfield*sizeof(double),
1110 (char *) &newfields[0][0]);
1111 irregular->destroy_data();
1112 delete irregular;
1113
1114 memory->destroy(fields);
1115 memory->destroy(procassign);
1116
1117 // point fields at newfields
1118
1119 fields = newfields;
1120 maxnew = newmaxnew;
1121 nnew = nrecv;
1122 }
1123
1124 /* ----------------------------------------------------------------------
1125 migrate final atoms to new procs based on atom coords
1126 use migrate_atoms() with implicit processor assignments based on atom coords
1127 move atoms back inside simulation box and to new processors
1128 use remap() instead of pbc() in case atoms moved a long distance
1129 adjust image flags of all atoms (old and new) based on current box
1130 ------------------------------------------------------------------------- */
1131
migrate_atoms_by_coords()1132 void ReadDump::migrate_atoms_by_coords()
1133 {
1134 double **x = atom->x;
1135 imageint *image = atom->image;
1136 int nlocal = atom->nlocal;
1137 for (int i = 0; i < nlocal; i++) domain->remap(x[i],image[i]);
1138
1139 if (triclinic) domain->x2lamda(atom->nlocal);
1140 domain->reset_box();
1141 Irregular *irregular = new Irregular(lmp);
1142 irregular->migrate_atoms(1);
1143 delete irregular;
1144 if (triclinic) domain->lamda2x(atom->nlocal);
1145 }
1146
1147 /* ----------------------------------------------------------------------
1148 process arg list for dump file fields and optional keywords
1149 ------------------------------------------------------------------------- */
1150
fields_and_keywords(int narg,char ** arg)1151 int ReadDump::fields_and_keywords(int narg, char **arg)
1152 {
1153 // per-field vectors, leave space for ID and TYPE
1154
1155 fieldtype = new int[narg+2];
1156 fieldlabel = new char*[narg+2];
1157
1158 // add id and type fields as needed
1159 // scan ahead to see if "add yes/keep" keyword/value is used
1160 // requires extra "type" field from from dump file
1161
1162 int iarg;
1163 for (iarg = 0; iarg < narg; iarg++)
1164 if (strcmp(arg[iarg],"add") == 0)
1165 if (iarg < narg-1 && (strcmp(arg[iarg+1],"yes") == 0 ||
1166 strcmp(arg[iarg+1],"keep") == 0)) break;
1167
1168 nfield = 0;
1169 fieldtype[nfield++] = ID;
1170 if (iarg < narg) fieldtype[nfield++] = TYPE;
1171
1172 // parse fields
1173
1174 iarg = 0;
1175 while (iarg < narg) {
1176 int type = whichtype(arg[iarg]);
1177 if (type < 0) break;
1178 if (type == Q && !atom->q_flag)
1179 error->all(FLERR,"Read dump of atom property that isn't allocated");
1180 fieldtype[nfield++] = type;
1181 iarg++;
1182 }
1183
1184 // check for no fields
1185
1186 if (fieldtype[nfield-1] == ID || fieldtype[nfield-1] == TYPE)
1187 error->all(FLERR,"Illegal read_dump command");
1188
1189 if (dimension == 2) {
1190 for (int i = 0; i < nfield; i++)
1191 if (fieldtype[i] == Z || fieldtype[i] == VZ ||
1192 fieldtype[i] == IZ || fieldtype[i] == FZ)
1193 error->all(FLERR,"Illegal read_dump command");
1194 }
1195
1196 for (int i = 0; i < nfield; i++)
1197 for (int j = i+1; j < nfield; j++)
1198 if (fieldtype[i] == fieldtype[j])
1199 error->all(FLERR,"Duplicate fields in read_dump command");
1200
1201 // parse optional args
1202
1203 multiproc_nfile = 0;
1204 boxflag = 1;
1205 replaceflag = 1;
1206 purgeflag = 0;
1207 trimflag = 0;
1208 addflag = NOADD;
1209 for (int i = 0; i < nfield; i++) fieldlabel[i] = nullptr;
1210 scaleflag = 0;
1211 wrapflag = 1;
1212
1213 while (iarg < narg) {
1214 if (strcmp(arg[iarg],"nfile") == 0) {
1215 if (iarg+2 > narg) error->all(FLERR,"Illegal read_dump command");
1216 multiproc_nfile = utils::inumeric(FLERR,arg[iarg+1],false,lmp);
1217 iarg += 2;
1218 } else if (strcmp(arg[iarg],"box") == 0) {
1219 if (iarg+2 > narg) error->all(FLERR,"Illegal read_dump command");
1220 if (strcmp(arg[iarg+1],"yes") == 0) boxflag = 1;
1221 else if (strcmp(arg[iarg+1],"no") == 0) boxflag = 0;
1222 else error->all(FLERR,"Illegal read_dump command");
1223 iarg += 2;
1224 } else if (strcmp(arg[iarg],"replace") == 0) {
1225 if (iarg+2 > narg) error->all(FLERR,"Illegal read_dump command");
1226 if (strcmp(arg[iarg+1],"yes") == 0) replaceflag = 1;
1227 else if (strcmp(arg[iarg+1],"no") == 0) replaceflag = 0;
1228 else error->all(FLERR,"Illegal read_dump command");
1229 iarg += 2;
1230 } else if (strcmp(arg[iarg],"purge") == 0) {
1231 if (iarg+2 > narg) error->all(FLERR,"Illegal read_dump command");
1232 if (strcmp(arg[iarg+1],"yes") == 0) purgeflag = 1;
1233 else if (strcmp(arg[iarg+1],"no") == 0) purgeflag = 0;
1234 else error->all(FLERR,"Illegal read_dump command");
1235 iarg += 2;
1236 } else if (strcmp(arg[iarg],"trim") == 0) {
1237 if (iarg+2 > narg) error->all(FLERR,"Illegal read_dump command");
1238 if (strcmp(arg[iarg+1],"yes") == 0) trimflag = 1;
1239 else if (strcmp(arg[iarg+1],"no") == 0) trimflag = 0;
1240 else error->all(FLERR,"Illegal read_dump command");
1241 iarg += 2;
1242 } else if (strcmp(arg[iarg],"add") == 0) {
1243 if (iarg+2 > narg) error->all(FLERR,"Illegal read_dump command");
1244 if (strcmp(arg[iarg+1],"yes") == 0) addflag = YESADD;
1245 else if (strcmp(arg[iarg+1],"no") == 0) addflag = NOADD;
1246 else if (strcmp(arg[iarg+1],"keep") == 0) addflag = KEEPADD;
1247 else error->all(FLERR,"Illegal read_dump command");
1248 iarg += 2;
1249 } else if (strcmp(arg[iarg],"label") == 0) {
1250 if (iarg+3 > narg) error->all(FLERR,"Illegal read_dump command");
1251 int type = whichtype(arg[iarg+1]);
1252 int i;
1253 for (i = 0; i < nfield; i++)
1254 if (type == fieldtype[i]) break;
1255 if (i == nfield) error->all(FLERR,"Illegal read_dump command");
1256 fieldlabel[i] = utils::strdup(arg[iarg+2]);
1257 iarg += 3;
1258 } else if (strcmp(arg[iarg],"scaled") == 0) {
1259 if (iarg+2 > narg) error->all(FLERR,"Illegal read_dump command");
1260 if (strcmp(arg[iarg+1],"yes") == 0) scaleflag = 1;
1261 else if (strcmp(arg[iarg+1],"no") == 0) scaleflag = 0;
1262 else error->all(FLERR,"Illegal read_dump command");
1263 iarg += 2;
1264 } else if (strcmp(arg[iarg],"wrapped") == 0) {
1265 if (iarg+2 > narg) error->all(FLERR,"Illegal read_dump command");
1266 if (strcmp(arg[iarg+1],"yes") == 0) wrapflag = 1;
1267 else if (strcmp(arg[iarg+1],"no") == 0) wrapflag = 0;
1268 else error->all(FLERR,"Illegal read_dump command");
1269 iarg += 2;
1270 } else if (strcmp(arg[iarg],"format") == 0) {
1271 if (iarg+2 > narg) error->all(FLERR,"Illegal read_dump command");
1272 delete [] readerstyle;
1273 readerstyle = utils::strdup(arg[iarg+1]);
1274 iarg += 2;
1275 break;
1276 } else error->all(FLERR,"Illegal read_dump command");
1277 }
1278
1279 if (multiproc == 0 && multiproc_nfile)
1280 error->all(FLERR,"Dump file is not a multi-proc file");
1281 if (multiproc && multiproc_nfile == 0)
1282 error->all(FLERR,"Dump file is a multi-proc file");
1283
1284 if (purgeflag && (replaceflag || trimflag))
1285 error->all(FLERR,"If read_dump purges it cannot replace or trim");
1286 if (addflag == KEEPADD && atom->tag_enable == 0)
1287 error->all(FLERR,"Read_dump cannot use 'add keep' without atom IDs");
1288
1289 return narg-iarg;
1290 }
1291
1292 /* ----------------------------------------------------------------------
1293 check if str is a field argument
1294 if yes, return index of which
1295 if not, return -1
1296 ------------------------------------------------------------------------- */
1297
whichtype(char * str)1298 int ReadDump::whichtype(char *str)
1299 {
1300 int type = -1;
1301 if (strcmp(str,"id") == 0) type = ID;
1302 else if (strcmp(str,"type") == 0) type = TYPE;
1303 else if (strcmp(str,"x") == 0) type = X;
1304 else if (strcmp(str,"y") == 0) type = Y;
1305 else if (strcmp(str,"z") == 0) type = Z;
1306 else if (strcmp(str,"vx") == 0) type = VX;
1307 else if (strcmp(str,"vy") == 0) type = VY;
1308 else if (strcmp(str,"vz") == 0) type = VZ;
1309 else if (strcmp(str,"q") == 0) type = Q;
1310 else if (strcmp(str,"ix") == 0) type = IX;
1311 else if (strcmp(str,"iy") == 0) type = IY;
1312 else if (strcmp(str,"iz") == 0) type = IZ;
1313 else if (strcmp(str,"fx") == 0) type = FX;
1314 else if (strcmp(str,"fy") == 0) type = FY;
1315 else if (strcmp(str,"fz") == 0) type = FZ;
1316 return type;
1317 }
1318
1319 /* ----------------------------------------------------------------------
1320 convert XYZ fields in dump file into absolute, unscaled coordinates
1321 depends on scaled vs unscaled and triclinic vs orthogonal
1322 does not depend on wrapped vs unwrapped
1323 ------------------------------------------------------------------------- */
1324
xfield(int i,int j)1325 double ReadDump::xfield(int i, int j)
1326 {
1327 if (!scaled) return fields[i][j];
1328 else if (!triclinic) return fields[i][j]*xprd + xlo;
1329 else if (dimension == 2)
1330 return xprd*fields[i][j] + xy*fields[i][yindex] + xlo;
1331 return xprd*fields[i][j] + xy*fields[i][yindex] + xz*fields[i][zindex] + xlo;
1332 }
1333
yfield(int i,int j)1334 double ReadDump::yfield(int i, int j)
1335 {
1336 if (!scaled) return fields[i][j];
1337 else if (!triclinic) return fields[i][j]*yprd + ylo;
1338 else if (dimension == 2) return yprd*fields[i][j] + ylo;
1339 return yprd*fields[i][j] + yz*fields[i][zindex] + ylo;
1340 }
1341
zfield(int i,int j)1342 double ReadDump::zfield(int i, int j)
1343 {
1344 if (!scaled) return fields[i][j];
1345 return fields[i][j]*zprd + zlo;
1346 }
1347