diff --git a/src/atom.cpp b/src/atom.cpp
index eac098d953034dcf515719ccf3d0dcf25c65c363..6d7bb5b4edd032c989e2a78fa3609d4fa1e2e68c 100644
--- a/src/atom.cpp
+++ b/src/atom.cpp
@@ -950,7 +950,8 @@ void Atom::data_vels(int n, char *buf, tagint id_offset)
check that atom IDs are > 0 and <= map_tag_max
------------------------------------------------------------------------- */
-void Atom::data_bonds(int n, char *buf, int *count, tagint id_offset)
+void Atom::data_bonds(int n, char *buf, int *count, tagint id_offset,
+ int type_offset)
{
int m,tmp,itype;
tagint atom1,atom2;
@@ -966,6 +967,7 @@ void Atom::data_bonds(int n, char *buf, int *count, tagint id_offset)
atom1 += id_offset;
atom2 += id_offset;
}
+ itype += type_offset;
if (atom1 <= 0 || atom1 > map_tag_max ||
atom2 <= 0 || atom2 > map_tag_max)
@@ -1001,7 +1003,8 @@ void Atom::data_bonds(int n, char *buf, int *count, tagint id_offset)
check that atom IDs are > 0 and <= map_tag_max
------------------------------------------------------------------------- */
-void Atom::data_angles(int n, char *buf, int *count, tagint id_offset)
+void Atom::data_angles(int n, char *buf, int *count, tagint id_offset,
+ int type_offset)
{
int m,tmp,itype;
tagint atom1,atom2,atom3;
@@ -1018,6 +1021,7 @@ void Atom::data_angles(int n, char *buf, int *count, tagint id_offset)
atom2 += id_offset;
atom3 += id_offset;
}
+ itype += type_offset;
if (atom1 <= 0 || atom1 > map_tag_max ||
atom2 <= 0 || atom2 > map_tag_max ||
@@ -1068,7 +1072,8 @@ void Atom::data_angles(int n, char *buf, int *count, tagint id_offset)
check that atom IDs are > 0 and <= map_tag_max
------------------------------------------------------------------------- */
-void Atom::data_dihedrals(int n, char *buf, int *count, tagint id_offset)
+void Atom::data_dihedrals(int n, char *buf, int *count, tagint id_offset,
+ int type_offset)
{
int m,tmp,itype;
tagint atom1,atom2,atom3,atom4;
@@ -1087,6 +1092,7 @@ void Atom::data_dihedrals(int n, char *buf, int *count, tagint id_offset)
atom3 += id_offset;
atom4 += id_offset;
}
+ itype += type_offset;
if (atom1 <= 0 || atom1 > map_tag_max ||
atom2 <= 0 || atom2 > map_tag_max ||
@@ -1153,7 +1159,8 @@ void Atom::data_dihedrals(int n, char *buf, int *count, tagint id_offset)
check that atom IDs are > 0 and <= map_tag_max
------------------------------------------------------------------------- */
-void Atom::data_impropers(int n, char *buf, int *count, tagint id_offset)
+void Atom::data_impropers(int n, char *buf, int *count, tagint id_offset,
+ int type_offset)
{
int m,tmp,itype;
tagint atom1,atom2,atom3,atom4;
@@ -1172,6 +1179,7 @@ void Atom::data_impropers(int n, char *buf, int *count, tagint id_offset)
atom3 += id_offset;
atom4 += id_offset;
}
+ itype += type_offset;
if (atom1 <= 0 || atom1 > map_tag_max ||
atom2 <= 0 || atom2 > map_tag_max ||
@@ -1496,17 +1504,17 @@ void Atom::add_molecule(int narg, char **arg)
if (find_molecule(arg[0]) >= 0)
error->all(FLERR,"Reuse of molecule template ID");
- // may over-allocate if not all args are mol files, but OK for srealloc
-
- molecules = (Molecule **)
- memory->srealloc(molecules,(nmolecule+narg-1)*sizeof(Molecule *),
- "atom::molecules");
-
// 1st molecule in set stores nset = # of mols, others store nset = 0
+ // ifile = count of molecules in set
+ // index = argument index where next molecule starts, updated by constructor
int ifile = 1;
+ int index = 1;
while (1) {
- molecules[nmolecule] = new Molecule(lmp,narg,arg,ifile);
+ molecules = (Molecule **)
+ memory->srealloc(molecules,(nmolecule+1)*sizeof(Molecule *),
+ "atom::molecules");
+ molecules[nmolecule] = new Molecule(lmp,narg,arg,index);
molecules[nmolecule]->nset = 0;
molecules[nmolecule-ifile+1]->nset++;
nmolecule++;
diff --git a/src/atom.h b/src/atom.h
index 1583dada7e852af0f3368908f6730a61b436f927..dd0d9821a1079d09ea35e018b8c67681a40fd064 100644
--- a/src/atom.h
+++ b/src/atom.h
@@ -213,10 +213,10 @@ class Atom : protected Pointers {
void data_atoms(int, char *, tagint, int, int, double *);
void data_vels(int, char *, tagint);
- void data_bonds(int, char *, int *, tagint);
- void data_angles(int, char *, int *, tagint);
- void data_dihedrals(int, char *, int *, tagint);
- void data_impropers(int, char *, int *, tagint);
+ void data_bonds(int, char *, int *, tagint, int);
+ void data_angles(int, char *, int *, tagint, int);
+ void data_dihedrals(int, char *, int *, tagint, int);
+ void data_impropers(int, char *, int *, tagint, int);
void data_bonus(int, char *, class AtomVec *, tagint);
void data_bodies(int, char *, class AtomVecBody *, tagint);
diff --git a/src/molecule.cpp b/src/molecule.cpp
index ee8f9ae424b0a9f8ba7e3f3aa18e66d09d0506bd..83157200eba2e9f3939fae24efeb16c9ee2c06ad 100644
--- a/src/molecule.cpp
+++ b/src/molecule.cpp
@@ -35,11 +35,12 @@ using namespace MathConst;
/* ---------------------------------------------------------------------- */
-Molecule::Molecule(LAMMPS *lmp, int narg, char **arg, int ifile) : Pointers(lmp)
+Molecule::Molecule(LAMMPS *lmp, int narg, char **arg, int &index) :
+ Pointers(lmp)
{
me = comm->me;
- if (ifile >= narg) error->all(FLERR,"Illegal molecule command");
+ if (index >= narg) error->all(FLERR,"Illegal molecule command");
int n = strlen(arg[0]) + 1;
id = new char[n];
@@ -50,21 +51,14 @@ Molecule::Molecule(LAMMPS *lmp, int narg, char **arg, int ifile) : Pointers(lmp)
error->all(FLERR,"Molecule template ID must be "
"alphanumeric or underscore characters");
- // scan args past ifile to reach optional args
- // set last = 1 if no more files in list
-
- last = 0;
- int iarg = ifile+1;
- while (iarg < narg) {
- if (strcmp(arg[iarg],"offset") == 0) break;
- iarg++;
- }
- if (iarg == ifile+1) last = 1;
-
- // parse optional args
+ // parse args until reach unknown arg (next file)
toffset = 0;
boffset = aoffset = doffset = ioffset = 0;
+ sizescale = 1.0;
+
+ int ifile = index;
+ int iarg = ifile+1;
while (iarg < narg) {
if (strcmp(arg[iarg],"offset") == 0) {
@@ -78,9 +72,45 @@ Molecule::Molecule(LAMMPS *lmp, int narg, char **arg, int ifile) : Pointers(lmp)
doffset < 0 || ioffset < 0)
error->all(FLERR,"Illegal molecule command");
iarg += 6;
- } else error->all(FLERR,"Illegal molecule command");
+ } else if (strcmp(arg[iarg],"toff") == 0) {
+ if (iarg+2 > narg) error->all(FLERR,"Illegal molecule command");
+ toffset = force->inumeric(FLERR,arg[iarg+1]);
+ if (toffset < 0) error->all(FLERR,"Illegal molecule command");
+ iarg += 2;
+ } else if (strcmp(arg[iarg],"boff") == 0) {
+ if (iarg+2 > narg) error->all(FLERR,"Illegal molecule command");
+ boffset = force->inumeric(FLERR,arg[iarg+1]);
+ if (boffset < 0) error->all(FLERR,"Illegal molecule command");
+ iarg += 2;
+ } else if (strcmp(arg[iarg],"aoff") == 0) {
+ if (iarg+2 > narg) error->all(FLERR,"Illegal molecule command");
+ aoffset = force->inumeric(FLERR,arg[iarg+1]);
+ if (aoffset < 0) error->all(FLERR,"Illegal molecule command");
+ iarg += 2;
+ } else if (strcmp(arg[iarg],"doff") == 0) {
+ if (iarg+2 > narg) error->all(FLERR,"Illegal molecule command");
+ doffset = force->inumeric(FLERR,arg[iarg+1]);
+ if (doffset < 0) error->all(FLERR,"Illegal molecule command");
+ iarg += 2;
+ } else if (strcmp(arg[iarg],"ioff") == 0) {
+ if (iarg+2 > narg) error->all(FLERR,"Illegal molecule command");
+ ioffset = force->inumeric(FLERR,arg[iarg+1]);
+ if (ioffset < 0) error->all(FLERR,"Illegal molecule command");
+ iarg += 2;
+ } else if (strcmp(arg[iarg],"scale") == 0) {
+ if (iarg+2 > narg) error->all(FLERR,"Illegal molecule command");
+ sizescale = force->numeric(FLERR,arg[iarg+1]);
+ if (sizescale <= 0.0) error->all(FLERR,"Illegal molecule command");
+ iarg += 2;
+ } else break;
}
+ index = iarg;
+
+ // last molecule if have scanned all args
+
+ if (iarg == narg) last = 1;
+
// initialize all fields to empty
initialize();
@@ -392,10 +422,14 @@ void Molecule::read(int flag)
else if (strstr(line,"mass")) {
massflag = 1;
sscanf(line,"%lg",&masstotal);
+ masstotal *= sizescale*sizescale*sizescale;
}
else if (strstr(line,"com")) {
comflag = 1;
sscanf(line,"%lg %lg %lg",&com[0],&com[1],&com[2]);
+ com[0] *= sizescale;
+ com[1] *= sizescale;
+ com[2] *= sizescale;
if (domain->dimension == 2 && com[2] != 0.0)
error->all(FLERR,"Molecule file z center-of-mass must be 0.0 for 2d");
}
@@ -404,6 +438,12 @@ void Molecule::read(int flag)
sscanf(line,"%lg %lg %lg %lg %lg %lg",
&itensor[0],&itensor[1],&itensor[2],
&itensor[3],&itensor[4],&itensor[5]);
+ itensor[0] *= sizescale*sizescale*sizescale*sizescale*sizescale;
+ itensor[1] *= sizescale*sizescale*sizescale*sizescale*sizescale;
+ itensor[2] *= sizescale*sizescale*sizescale*sizescale*sizescale;
+ itensor[3] *= sizescale*sizescale*sizescale*sizescale*sizescale;
+ itensor[4] *= sizescale*sizescale*sizescale*sizescale*sizescale;
+ itensor[5] *= sizescale*sizescale*sizescale*sizescale*sizescale;
}
else break;
@@ -414,8 +454,10 @@ void Molecule::read(int flag)
if (natoms < 1) error->all(FLERR,"No or invalid atom count in molecule file");
if (nbonds < 0) error->all(FLERR,"Invalid bond count in molecule file");
if (nangles < 0) error->all(FLERR,"Invalid angle count in molecule file");
- if (ndihedrals < 0) error->all(FLERR,"Invalid dihedral count in molecule file");
- if (nimpropers < 0) error->all(FLERR,"Invalid improper count in molecule file");
+ if (ndihedrals < 0)
+ error->all(FLERR,"Invalid dihedral count in molecule file");
+ if (nimpropers < 0)
+ error->all(FLERR,"Invalid improper count in molecule file");
// count = vector for tallying bonds,angles,etc per atom
@@ -544,6 +586,9 @@ void Molecule::coords(char *line)
error->all(FLERR,"Invalid Coords section in molecule file");
}
sscanf(line,"%d %lg %lg %lg",&tmp,&x[i][0],&x[i][1],&x[i][2]);
+ x[i][0] *= sizescale;
+ x[i][1] *= sizescale;
+ x[i][2] *= sizescale;
}
if (domain->dimension == 2) {
@@ -614,6 +659,7 @@ void Molecule::diameters(char *line)
error->all(FLERR,"Invalid Diameters section in molecule file");
}
sscanf(line,"%d %lg",&tmp,&radius[i]);
+ radius[i] *= sizescale;
radius[i] *= 0.5;
maxradius = MAX(maxradius,radius[i]);
}
@@ -638,6 +684,7 @@ void Molecule::masses(char *line)
error->all(FLERR,"Invalid Masses section in molecule file");
}
sscanf(line,"%d %lg",&tmp,&rmass[i]);
+ rmass[i] *= sizescale*sizescale*sizescale;
}
for (int i = 0; i < natoms; i++)
diff --git a/src/molecule.h b/src/molecule.h
index 81d13d4e260418b957accf980f025de0fe4483cc..9ecb86b2029c737e47d3f43e2b5a90f47fce3c8c 100644
--- a/src/molecule.h
+++ b/src/molecule.h
@@ -102,7 +102,7 @@ class Molecule : protected Pointers {
double **dxbody; // displacement of each atom relative to COM
// in body frame (diagonalized interia tensor)
- Molecule(class LAMMPS *, int, char **, int);
+ Molecule(class LAMMPS *, int, char **, int &);
~Molecule();
void compute_center();
void compute_mass();
@@ -116,6 +116,7 @@ class Molecule : protected Pointers {
int *count;
int toffset,boffset,aoffset,doffset,ioffset;
int autospecial;
+ double sizescale;
void read(int);
void coords(char *);
diff --git a/src/read_data.cpp b/src/read_data.cpp
index 036e333bdb779b294d2c46f00071c3b616a8d759..7e4b70484b43d717c725a33e5da290cb707f9b6f 100644
--- a/src/read_data.cpp
+++ b/src/read_data.cpp
@@ -279,7 +279,7 @@ void ReadData::command(int narg, char **arg)
// -----------------------------------------------------------------
- // perform 1-pass read if no molecular topoogy in file
+ // perform 1-pass read if no molecular topology in file
// perform 2-pass read if molecular topology,
// first pass calculates max topology/atom
@@ -301,7 +301,7 @@ void ReadData::command(int narg, char **arg)
triclinic = 0;
keyword[0] = '\0';
- int nlocal_previous = atom->nlocal;
+ nlocal_previous = atom->nlocal;
int firstpass = 1;
while (1) {
@@ -315,7 +315,7 @@ void ReadData::command(int narg, char **arg)
// read header info
- header();
+ header(firstpass);
// problem setup using info from header
// only done once, if firstpass and first data file
@@ -648,7 +648,7 @@ void ReadData::command(int narg, char **arg)
// will also observe extra settings even if bond/etc topology not in file
// leaves other atom arrays unchanged, since already nmax in length
- atom->deallocate_topology();
+ if (addflag == NONE) atom->deallocate_topology();
atom->avec->grow(atom->nmax);
}
@@ -788,7 +788,7 @@ void ReadData::command(int narg, char **arg)
some logic differs if adding atoms
------------------------------------------------------------------------- */
-void ReadData::header()
+void ReadData::header(int firstpass)
{
int n;
char *ptr;
@@ -856,7 +856,7 @@ void ReadData::header()
if (strstr(line,"atoms")) {
sscanf(line,BIGINT_FORMAT,&natoms);
if (addflag == NONE) atom->natoms = natoms;
- else atom->natoms += natoms;
+ else if (firstpass) atom->natoms += natoms;
// check for these first
// otherwise "triangles" will be matched as "angles"
@@ -1093,6 +1093,8 @@ void ReadData::velocities()
void ReadData::bonds(int firstpass)
{
+ int nchunk,eof;
+
if (me == 0) {
if (firstpass) {
if (screen) fprintf(screen," scanning bonds ...\n");
@@ -1114,32 +1116,38 @@ void ReadData::bonds(int firstpass)
// read and process bonds
- int nchunk,eof;
bigint nread = 0;
- bigint nbonds = atom->nbonds;
while (nread < nbonds) {
nchunk = MIN(nbonds-nread,CHUNK);
eof = comm->read_lines_from_file(fp,nchunk,MAXLINE,buffer);
if (eof) error->all(FLERR,"Unexpected end of data file");
- atom->data_bonds(nchunk,buffer,count,id_offset);
+ atom->data_bonds(nchunk,buffer,count,id_offset,boffset);
nread += nchunk;
}
// if firstpass: tally max bond/atom and return
+ // if addflag = NONE, store max bond/atom with extra
+ // else just check actual max does not exceed existing max
if (firstpass) {
int max = 0;
- for (int i = 0; i < nlocal; i++) max = MAX(max,count[i]);
+ for (int i = nlocal_previous; i < nlocal; i++) max = MAX(max,count[i]);
int maxall;
MPI_Allreduce(&max,&maxall,1,MPI_INT,MPI_MAX,world);
- maxall += atom->extra_bond_per_atom;
+ if (addflag == NONE) maxall += atom->extra_bond_per_atom;
if (me == 0) {
if (screen) fprintf(screen," %d = max bonds/atom\n",maxall);
if (logfile) fprintf(logfile," %d = max bonds/atom\n",maxall);
}
- atom->bond_per_atom = maxall;
+
+ if (addflag != NONE) {
+ if (maxall > atom->bond_per_atom)
+ error->all(FLERR,"Subsequent read data induced "
+ "too many bonds per atom");
+ } else atom->bond_per_atom = maxall;
+
memory->destroy(count);
return;
}
@@ -1147,7 +1155,7 @@ void ReadData::bonds(int firstpass)
// if 2nd pass: check that bonds were assigned correctly
bigint n = 0;
- for (int i = 0; i < nlocal; i++) n += atom->num_bond[i];
+ for (int i = nlocal_previous; i < nlocal; i++) n += atom->num_bond[i];
bigint sum;
MPI_Allreduce(&n,&sum,1,MPI_LMP_BIGINT,MPI_SUM,world);
int factor = 1;
@@ -1158,7 +1166,7 @@ void ReadData::bonds(int firstpass)
if (logfile) fprintf(logfile," " BIGINT_FORMAT " bonds\n",sum/factor);
}
- if (sum != factor*atom->nbonds)
+ if (sum != factor*nbonds)
error->all(FLERR,"Bonds assigned incorrectly");
}
@@ -1168,6 +1176,8 @@ void ReadData::bonds(int firstpass)
void ReadData::angles(int firstpass)
{
+ int nchunk,eof;
+
if (me == 0) {
if (firstpass) {
if (screen) fprintf(screen," scanning angles ...\n");
@@ -1189,32 +1199,38 @@ void ReadData::angles(int firstpass)
// read and process angles
- int nchunk,eof;
bigint nread = 0;
- bigint nangles = atom->nangles;
while (nread < nangles) {
nchunk = MIN(nangles-nread,CHUNK);
eof = comm->read_lines_from_file(fp,nchunk,MAXLINE,buffer);
if (eof) error->all(FLERR,"Unexpected end of data file");
- atom->data_angles(nchunk,buffer,count,id_offset);
+ atom->data_angles(nchunk,buffer,count,id_offset,aoffset);
nread += nchunk;
}
// if firstpass: tally max angle/atom and return
+ // if addflag = NONE, store max angle/atom with extra
+ // else just check actual max does not exceed existing max
if (firstpass) {
int max = 0;
- for (int i = 0; i < nlocal; i++) max = MAX(max,count[i]);
+ for (int i = nlocal_previous; i < nlocal; i++) max = MAX(max,count[i]);
int maxall;
MPI_Allreduce(&max,&maxall,1,MPI_INT,MPI_MAX,world);
- maxall += atom->extra_angle_per_atom;
+ if (addflag == NONE) maxall += atom->extra_angle_per_atom;
if (me == 0) {
if (screen) fprintf(screen," %d = max angles/atom\n",maxall);
if (logfile) fprintf(logfile," %d = max angles/atom\n",maxall);
}
- atom->angle_per_atom = maxall;
+
+ if (addflag != NONE) {
+ if (maxall > atom->angle_per_atom)
+ error->all(FLERR,"Subsequent read data induced "
+ "too many angles per atom");
+ } else atom->angle_per_atom = maxall;
+
memory->destroy(count);
return;
}
@@ -1222,7 +1238,7 @@ void ReadData::angles(int firstpass)
// if 2nd pass: check that angles were assigned correctly
bigint n = 0;
- for (int i = 0; i < nlocal; i++) n += atom->num_angle[i];
+ for (int i = nlocal_previous; i < nlocal; i++) n += atom->num_angle[i];
bigint sum;
MPI_Allreduce(&n,&sum,1,MPI_LMP_BIGINT,MPI_SUM,world);
int factor = 1;
@@ -1233,7 +1249,7 @@ void ReadData::angles(int firstpass)
if (logfile) fprintf(logfile," " BIGINT_FORMAT " angles\n",sum/factor);
}
- if (sum != factor*atom->nangles)
+ if (sum != factor*nangles)
error->all(FLERR,"Angles assigned incorrectly");
}
@@ -1243,6 +1259,8 @@ void ReadData::angles(int firstpass)
void ReadData::dihedrals(int firstpass)
{
+ int nchunk,eof;
+
if (me == 0) {
if (firstpass) {
if (screen) fprintf(screen," scanning dihedrals ...\n");
@@ -1264,31 +1282,38 @@ void ReadData::dihedrals(int firstpass)
// read and process dihedrals
- int nchunk,eof;
bigint nread = 0;
- bigint ndihedrals = atom->ndihedrals;
while (nread < ndihedrals) {
nchunk = MIN(ndihedrals-nread,CHUNK);
eof = comm->read_lines_from_file(fp,nchunk,MAXLINE,buffer);
if (eof) error->all(FLERR,"Unexpected end of data file");
- atom->data_dihedrals(nchunk,buffer,count,id_offset);
+ atom->data_dihedrals(nchunk,buffer,count,id_offset,doffset);
nread += nchunk;
}
// if firstpass: tally max dihedral/atom and return
+ // if addflag = NONE, store max dihedral/atom with extra
+ // else just check actual max does not exceed existing max
if (firstpass) {
int max = 0;
for (int i = 0; i < nlocal; i++) max = MAX(max,count[i]);
int maxall;
MPI_Allreduce(&max,&maxall,1,MPI_INT,MPI_MAX,world);
- maxall += atom->extra_dihedral_per_atom;
+ if (addflag == NONE) maxall += atom->extra_dihedral_per_atom;
if (me == 0) {
if (screen) fprintf(screen," %d = max dihedrals/atom\n",maxall);
if (logfile) fprintf(logfile," %d = max dihedrals/atom\n",maxall);
}
+
+ if (addflag != NONE) {
+ if (maxall > atom->dihedral_per_atom)
+ error->all(FLERR,"Subsequent read data induced "
+ "too many dihedrals per atom");
+ } else atom->dihedral_per_atom = maxall;
+
atom->dihedral_per_atom = maxall;
memory->destroy(count);
return;
@@ -1297,7 +1322,7 @@ void ReadData::dihedrals(int firstpass)
// if 2nd pass: check that dihedrals were assigned correctly
bigint n = 0;
- for (int i = 0; i < nlocal; i++) n += atom->num_dihedral[i];
+ for (int i = nlocal_previous; i < nlocal; i++) n += atom->num_dihedral[i];
bigint sum;
MPI_Allreduce(&n,&sum,1,MPI_LMP_BIGINT,MPI_SUM,world);
int factor = 1;
@@ -1308,7 +1333,7 @@ void ReadData::dihedrals(int firstpass)
if (logfile) fprintf(logfile," " BIGINT_FORMAT " dihedrals\n",sum/factor);
}
- if (sum != factor*atom->ndihedrals)
+ if (sum != factor*ndihedrals)
error->all(FLERR,"Dihedrals assigned incorrectly");
}
@@ -1318,6 +1343,8 @@ void ReadData::dihedrals(int firstpass)
void ReadData::impropers(int firstpass)
{
+ int nchunk,eof;
+
if (me == 0) {
if (firstpass) {
if (screen) fprintf(screen," scanning impropers ...\n");
@@ -1339,32 +1366,38 @@ void ReadData::impropers(int firstpass)
// read and process impropers
- int nchunk,eof;
bigint nread = 0;
- bigint nimpropers = atom->nimpropers;
while (nread < nimpropers) {
nchunk = MIN(nimpropers-nread,CHUNK);
eof = comm->read_lines_from_file(fp,nchunk,MAXLINE,buffer);
if (eof) error->all(FLERR,"Unexpected end of data file");
- atom->data_impropers(nchunk,buffer,count,id_offset);
+ atom->data_impropers(nchunk,buffer,count,id_offset,ioffset);
nread += nchunk;
}
// if firstpass: tally max improper/atom and return
+ // if addflag = NONE, store max improper/atom
+ // else just check it does not exceed existing max
if (firstpass) {
int max = 0;
- for (int i = 0; i < nlocal; i++) max = MAX(max,count[i]);
+ for (int i = nlocal_previous; i < nlocal; i++) max = MAX(max,count[i]);
int maxall;
MPI_Allreduce(&max,&maxall,1,MPI_INT,MPI_MAX,world);
- maxall += atom->extra_improper_per_atom;
+ if (addflag == NONE) maxall += atom->extra_improper_per_atom;
if (me == 0) {
if (screen) fprintf(screen," %d = max impropers/atom\n",maxall);
if (logfile) fprintf(logfile," %d = max impropers/atom\n",maxall);
}
- atom->improper_per_atom = maxall;
+
+ if (addflag != NONE) {
+ if (maxall > atom->improper_per_atom)
+ error->all(FLERR,"Subsequent read data induced "
+ "too many impropers per atom");
+ } else atom->improper_per_atom = maxall;
+
memory->destroy(count);
return;
}
@@ -1372,7 +1405,7 @@ void ReadData::impropers(int firstpass)
// if 2nd pass: check that impropers were assigned correctly
bigint n = 0;
- for (int i = 0; i < nlocal; i++) n += atom->num_improper[i];
+ for (int i = nlocal_previous; i < nlocal; i++) n += atom->num_improper[i];
bigint sum;
MPI_Allreduce(&n,&sum,1,MPI_LMP_BIGINT,MPI_SUM,world);
int factor = 1;
@@ -1383,7 +1416,7 @@ void ReadData::impropers(int firstpass)
if (logfile) fprintf(logfile," " BIGINT_FORMAT " impropers\n",sum/factor);
}
- if (sum != factor*atom->nimpropers)
+ if (sum != factor*nimpropers)
error->all(FLERR,"Impropers assigned incorrectly");
}
diff --git a/src/read_data.h b/src/read_data.h
index 74c3e2c8b9545ca0bac84c9c7ce50835b073e335..925ec1f07eb071d4b1bb7ed3acfc41ca99aa9665 100644
--- a/src/read_data.h
+++ b/src/read_data.h
@@ -41,6 +41,7 @@ class ReadData : protected Pointers {
bigint id_offset;
+ int nlocal_previous;
bigint natoms;
bigint nbonds,nangles,ndihedrals,nimpropers;
int ntypes;
@@ -81,7 +82,7 @@ class ReadData : protected Pointers {
void open(char *);
void scan(int &, int &, int &, int &);
int reallocate(int **, int, int);
- void header();
+ void header(int);
void parse_keyword(int);
void skip_lines(bigint);
void parse_coeffs(char *, const char *, int, int, int);