diff --git a/src/atom.cpp b/src/atom.cpp index 5b0908221f..b2abc07511 100644 --- a/src/atom.cpp +++ b/src/atom.cpp @@ -516,7 +516,7 @@ void Atom::modify_params(int narg, char **arg) error->all(FLERR, "Atom_modify id command after simulation box is defined"); if (strcmp(arg[iarg+1],"yes") == 0) tag_enable = 1; - else if (strcmp(arg[iarg+1],"no") == 0) tag_enable = 2; + else if (strcmp(arg[iarg+1],"no") == 0) tag_enable = 0; else error->all(FLERR,"Illegal atom_modify command"); iarg += 2; } else if (strcmp(arg[iarg],"map") == 0) { @@ -559,9 +559,11 @@ void Atom::modify_params(int narg, char **arg) check that atom IDs are valid error if any atom ID < 0 or atom ID = MAXTAGINT if any atom ID > 0, error if any atom ID == 0 + if any atom ID > 0, error if tag_enable = 0 if all atom IDs = 0, tag_enable must be 0 - OK if atom IDs > natoms - NOTE: not checking that atom IDs are unique + if max atom IDs < natoms, must be duplicates + OK if max atom IDs > natoms + NOTE: not fully checking that atom IDs are unique ------------------------------------------------------------------------- */ void Atom::tag_check() @@ -578,12 +580,16 @@ void Atom::tag_check() MPI_Allreduce(&min,&minall,1,MPI_LMP_TAGINT,MPI_MIN,world); MPI_Allreduce(&max,&maxall,1,MPI_LMP_TAGINT,MPI_MAX,world); - if (minall < 0) error->all(FLERR,"Atom ID is negative"); - if (maxall >= MAXTAGINT) error->all(FLERR,"Atom ID is too big"); - if (maxall > 0 && minall == 0) error->all(FLERR,"Atom ID is zero"); - // this last message is wrong - if (maxall == 0 && tag_enable && natoms) - error->all(FLERR,"Not all atom IDs are 0"); + if (minall < 0) error->all(FLERR,"One or more Atom IDs is negative"); + if (maxall >= MAXTAGINT) error->all(FLERR,"One or more atom IDs is too big"); + if (maxall > 0 && minall == 0) + error->all(FLERR,"One or more atom IDs is zero"); + if (maxall > 0 && tag_enable == 0) + error->all(FLERR,"Non-zero atom IDs with atom_modify id = no"); + if (maxall == 0 && natoms && tag_enable) + error->all(FLERR,"All atom IDs = 0 but atom_modify id = yes"); + if (tag_enable && maxall < natoms) + error->all(FLERR,"Duplicate atom IDs exist"); } /* ---------------------------------------------------------------------- @@ -725,7 +731,8 @@ void Atom::deallocate_topology() call style-specific routine to parse line ------------------------------------------------------------------------- */ -void Atom::data_atoms(int n, char *buf) +void Atom::data_atoms(int n, char *buf, tagint id_offset, int type_offset, + int shiftflag, double *shift) { int m,xptr,iptr; imageint imagedata; @@ -833,6 +840,12 @@ void Atom::data_atoms(int n, char *buf) xdata[0] = atof(values[xptr]); xdata[1] = atof(values[xptr+1]); xdata[2] = atof(values[xptr+2]); + if (shiftflag) { + xdata[0] += shift[0]; + xdata[1] += shift[1]; + xdata[2] += shift[2]; + } + domain->remap(xdata,imagedata); if (triclinic) { domain->x2lamda(xdata,lamda); @@ -841,8 +854,15 @@ void Atom::data_atoms(int n, char *buf) if (coord[0] >= sublo[0] && coord[0] < subhi[0] && coord[1] >= sublo[1] && coord[1] < subhi[1] && - coord[2] >= sublo[2] && coord[2] < subhi[2]) + coord[2] >= sublo[2] && coord[2] < subhi[2]) { avec->data_atom(xdata,imagedata,values); + if (id_offset) tag[nlocal-1] += id_offset; + if (type_offset) { + type[nlocal-1] += type_offset; + if (type[nlocal-1] > ntypes) + error->one(FLERR,"Invalid atom type in Atoms section of data file"); + } + } buf = next + 1; } @@ -856,7 +876,7 @@ void Atom::data_atoms(int n, char *buf) call style-specific routine to parse line ------------------------------------------------------------------------- */ -void Atom::data_vels(int n, char *buf) +void Atom::data_vels(int n, char *buf, tagint id_offset) { int j,m; tagint tagdata; @@ -883,7 +903,7 @@ void Atom::data_vels(int n, char *buf) for (j = 1; j < nwords; j++) values[j] = strtok(NULL," \t\n\r\f"); - tagdata = ATOTAGINT(values[0]); + tagdata = ATOTAGINT(values[0]) + id_offset; if (tagdata <= 0 || tagdata > map_tag_max) error->one(FLERR,"Invalid atom ID in Velocities section of data file"); if ((m = map(tagdata)) >= 0) avec->data_vel(m,&values[1]); @@ -901,7 +921,7 @@ void Atom::data_vels(int n, char *buf) check that atom IDs are > 0 and <= map_tag_max ------------------------------------------------------------------------- */ -void Atom::data_bonds(int n, char *buf, int *count) +void Atom::data_bonds(int n, char *buf, int *count, tagint id_offset) { int m,tmp,itype; tagint atom1,atom2; @@ -913,6 +933,11 @@ void Atom::data_bonds(int n, char *buf, int *count) *next = '\0'; sscanf(buf,"%d %d " TAGINT_FORMAT " " TAGINT_FORMAT, &tmp,&itype,&atom1,&atom2); + if (id_offset) { + atom1 += id_offset; + atom2 += id_offset; + } + if (atom1 <= 0 || atom1 > map_tag_max || atom2 <= 0 || atom2 > map_tag_max) error->one(FLERR,"Invalid atom ID in Bonds section of data file"); @@ -947,7 +972,7 @@ void Atom::data_bonds(int n, char *buf, int *count) check that atom IDs are > 0 and <= map_tag_max ------------------------------------------------------------------------- */ -void Atom::data_angles(int n, char *buf, int *count) +void Atom::data_angles(int n, char *buf, int *count, tagint id_offset) { int m,tmp,itype; tagint atom1,atom2,atom3; @@ -959,6 +984,12 @@ void Atom::data_angles(int n, char *buf, int *count) *next = '\0'; sscanf(buf,"%d %d " TAGINT_FORMAT " " TAGINT_FORMAT " " TAGINT_FORMAT, &tmp,&itype,&atom1,&atom2,&atom3); + if (id_offset) { + atom1 += id_offset; + atom2 += id_offset; + atom3 += id_offset; + } + if (atom1 <= 0 || atom1 > map_tag_max || atom2 <= 0 || atom2 > map_tag_max || atom3 <= 0 || atom3 > map_tag_max) @@ -1008,7 +1039,7 @@ void Atom::data_angles(int n, char *buf, int *count) check that atom IDs are > 0 and <= map_tag_max ------------------------------------------------------------------------- */ -void Atom::data_dihedrals(int n, char *buf, int *count) +void Atom::data_dihedrals(int n, char *buf, int *count, tagint id_offset) { int m,tmp,itype; tagint atom1,atom2,atom3,atom4; @@ -1021,6 +1052,13 @@ void Atom::data_dihedrals(int n, char *buf, int *count) sscanf(buf,"%d %d " TAGINT_FORMAT " " TAGINT_FORMAT " " TAGINT_FORMAT " " TAGINT_FORMAT, &tmp,&itype,&atom1,&atom2,&atom3,&atom4); + if (id_offset) { + atom1 += id_offset; + atom2 += id_offset; + atom3 += id_offset; + atom4 += id_offset; + } + if (atom1 <= 0 || atom1 > map_tag_max || atom2 <= 0 || atom2 > map_tag_max || atom3 <= 0 || atom3 > map_tag_max || @@ -1086,7 +1124,7 @@ void Atom::data_dihedrals(int n, char *buf, int *count) check that atom IDs are > 0 and <= map_tag_max ------------------------------------------------------------------------- */ -void Atom::data_impropers(int n, char *buf, int *count) +void Atom::data_impropers(int n, char *buf, int *count, tagint id_offset) { int m,tmp,itype; tagint atom1,atom2,atom3,atom4; @@ -1099,6 +1137,13 @@ void Atom::data_impropers(int n, char *buf, int *count) sscanf(buf,"%d %d " TAGINT_FORMAT " " TAGINT_FORMAT " " TAGINT_FORMAT " " TAGINT_FORMAT, &tmp,&itype,&atom1,&atom2,&atom3,&atom4); + if (id_offset) { + atom1 += id_offset; + atom2 += id_offset; + atom3 += id_offset; + atom4 += id_offset; + } + if (atom1 <= 0 || atom1 > map_tag_max || atom2 <= 0 || atom2 > map_tag_max || atom3 <= 0 || atom3 > map_tag_max || @@ -1163,7 +1208,7 @@ void Atom::data_impropers(int n, char *buf, int *count) call style-specific routine to parse line ------------------------------------------------------------------------- */ -void Atom::data_bonus(int n, char *buf, AtomVec *avec_bonus) +void Atom::data_bonus(int n, char *buf, AtomVec *avec_bonus, tagint id_offset) { int j,m,tagdata; char *next; @@ -1189,7 +1234,7 @@ void Atom::data_bonus(int n, char *buf, AtomVec *avec_bonus) for (j = 1; j < nwords; j++) values[j] = strtok(NULL," \t\n\r\f"); - tagdata = ATOTAGINT(values[0]); + tagdata = ATOTAGINT(values[0]) + id_offset; if (tagdata <= 0 || tagdata > map_tag_max) error->one(FLERR,"Invalid atom ID in Bonus section of data file"); @@ -1210,7 +1255,8 @@ void Atom::data_bonus(int n, char *buf, AtomVec *avec_bonus) call style-specific routine to parse line ------------------------------------------------------------------------- */ -void Atom::data_bodies(int n, char *buf, AtomVecBody *avec_body) +void Atom::data_bodies(int n, char *buf, AtomVecBody *avec_body, + tagint id_offset) { int j,m,tagdata,ninteger,ndouble; @@ -1222,8 +1268,8 @@ void Atom::data_bodies(int n, char *buf, AtomVecBody *avec_body) // if I own atom tag, unpack its values for (int i = 0; i < n; i++) { - if (i == 0) tagdata = ATOTAGINT(strtok(buf," \t\n\r\f")); - else tagdata = ATOTAGINT(strtok(NULL," \t\n\r\f")); + if (i == 0) tagdata = ATOTAGINT(strtok(buf," \t\n\r\f")) + id_offset; + else tagdata = ATOTAGINT(strtok(NULL," \t\n\r\f")) + id_offset; ninteger = atoi(strtok(NULL," \t\n\r\f")); ndouble = atoi(strtok(NULL," \t\n\r\f")); @@ -1260,9 +1306,10 @@ void Atom::allocate_type_arrays() /* ---------------------------------------------------------------------- set a mass and flag it as set called from reading of data file + type_offset may be used when reading multiple data files ------------------------------------------------------------------------- */ -void Atom::set_mass(const char *str) +void Atom::set_mass(const char *str, int type_offset) { if (mass == NULL) error->all(FLERR,"Cannot set mass for this atom style"); @@ -1270,6 +1317,7 @@ void Atom::set_mass(const char *str) double mass_one; int n = sscanf(str,"%d %lg",&itype,&mass_one); if (n != 2) error->all(FLERR,"Invalid mass line in data file"); + itype += type_offset; if (itype < 1 || itype > ntypes) error->all(FLERR,"Invalid type for mass set"); diff --git a/src/atom.h b/src/atom.h index 29da7d8797..afd07b1bd2 100644 --- a/src/atom.h +++ b/src/atom.h @@ -188,19 +188,19 @@ class Atom : protected Pointers { void deallocate_topology(); - void data_atoms(int, char *); - void data_vels(int, char *); + void data_atoms(int, char *, tagint, int, int, double *); + void data_vels(int, char *, tagint); - void data_bonds(int, char *, int *); - void data_angles(int, char *, int *); - void data_dihedrals(int, char *, int *); - void data_impropers(int, char *, int *); + 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_bonus(int, char *, class AtomVec *); - void data_bodies(int, char *, class AtomVecBody *); + void data_bonus(int, char *, class AtomVec *, tagint); + void data_bodies(int, char *, class AtomVecBody *, tagint); virtual void allocate_type_arrays(); - void set_mass(const char *); + void set_mass(const char *, int); void set_mass(int, double); void set_mass(int, char **); void set_mass(double *); diff --git a/src/group.cpp b/src/group.cpp index 738c9dde9d..f0d16c6df6 100644 --- a/src/group.cpp +++ b/src/group.cpp @@ -587,6 +587,26 @@ int Group::find(const char *name) return -1; } +/* ---------------------------------------------------------------------- + find group with name or create group if it doesn't exist + return group index +------------------------------------------------------------------------- */ + +int Group::find_or_create(const char *name) +{ + int igroup = find(name); + if (igroup >= 0) return igroup; + + if (ngroup == MAX_GROUP) error->all(FLERR,"Too many groups"); + igroup = find_unused(); + int n = strlen(name) + 1; + names[igroup] = new char[n]; + strcpy(names[igroup],name); + ngroup++; + + return igroup; +} + /* ---------------------------------------------------------------------- return index of first available group should never be called when group limit has been reached diff --git a/src/group.h b/src/group.h index 6def01ce55..cef184102e 100644 --- a/src/group.h +++ b/src/group.h @@ -33,6 +33,7 @@ class Group : protected Pointers { void assign(int, char **); // assign atoms to a group void create(char *, int *); // add flagged atoms to a group int find(const char *); // lookup name in list of groups + int find_or_create(const char *); // lookup name or create new group void write_restart(FILE *); void read_restart(FILE *); diff --git a/src/read_data.cpp b/src/read_data.cpp index 65b363b167..3a90b4c066 100644 --- a/src/read_data.cpp +++ b/src/read_data.cpp @@ -29,6 +29,7 @@ #include "atom_vec_tri.h" #include "force.h" #include "molecule.h" +#include "group.h" #include "comm.h" #include "update.h" #include "modify.h" @@ -41,6 +42,7 @@ #include "dihedral.h" #include "improper.h" #include "special.h" +#include "irregular.h" #include "error.h" #include "memory.h" @@ -55,6 +57,8 @@ using namespace LAMMPS_NS; // customize for new sections #define NSECTIONS 25 // change when add to header::section_keywords +enum{NONE,APPEND,VALUE,MERGE}; + // pair style suffixes to ignore // when matching Pair Coeffs comment to currently-defined pair style @@ -115,8 +119,16 @@ void ReadData::command(int narg, char **arg) // optional args - addflag = mergeflag = 0; - offset[0] = offset[1] = offset[2] = 0.0; + addflag = NONE; + id_offset = 0; + offsetflag = shiftflag = 0; + toffset = boffset = aoffset = doffset = ioffset = 0; + shift[0] = shift[1] = shift[2] = 0.0; + extra_atom_types = extra_bond_types = extra_angle_types = + extra_dihedral_types = extra_improper_types = 0; + + groupbit = 0; + nfix = 0; fix_index = NULL; fix_header = NULL; @@ -125,18 +137,81 @@ void ReadData::command(int narg, char **arg) int iarg = 1; while (iarg < narg) { if (strcmp(arg[iarg],"add") == 0) { - addflag = 1; - iarg++; - } else if (strcmp(arg[iarg],"merge") == 0) { - mergeflag = 1; - iarg++; + if (iarg+2 > narg) error->all(FLERR,"Illegal read_data command"); + if (strcmp(arg[iarg+1],"append") == 0) addflag = APPEND; + else if (strcmp(arg[iarg+1],"merge") == 0) addflag = MERGE; + else { + addflag = VALUE; + bigint offset = force->bnumeric(FLERR,arg[iarg+1]); + if (offset > MAXTAGINT) + error->all(FLERR,"Read data add offset is too big"); + id_offset = offset; + } + iarg += 2; } else if (strcmp(arg[iarg],"offset") == 0) { - if (iarg+4 > narg) + if (iarg+6 > narg) error->all(FLERR,"Illegal read_data command"); + offsetflag = 1; + toffset = force->inumeric(FLERR,arg[iarg+1]); + boffset = force->inumeric(FLERR,arg[iarg+2]); + aoffset = force->inumeric(FLERR,arg[iarg+3]); + doffset = force->inumeric(FLERR,arg[iarg+4]); + ioffset = force->inumeric(FLERR,arg[iarg+5]); + if (toffset < 0 || boffset < 0 || aoffset < 0 || + doffset < 0 || ioffset < 0) error->all(FLERR,"Illegal read_data command"); - offset[0] = force->numeric(FLERR,arg[iarg+1]); - offset[1] = force->numeric(FLERR,arg[iarg+2]); - offset[2] = force->numeric(FLERR,arg[iarg+3]); + iarg += 6; + } else if (strcmp(arg[iarg],"shift") == 0) { + if (iarg+4 > narg) error->all(FLERR,"Illegal read_data command"); + shiftflag = 1; + shift[0] = force->numeric(FLERR,arg[iarg+1]); + shift[1] = force->numeric(FLERR,arg[iarg+2]); + shift[2] = force->numeric(FLERR,arg[iarg+3]); + if (domain->dimension == 2 && shift[2] != 0.0) + error->all(FLERR,"Non-zero read_data shift z value for 2d simulation"); iarg += 4; + + } else if (strcmp(arg[iarg],"extra/atom/types") == 0) { + if (iarg+2 > narg) error->all(FLERR,"Illegal read_data command"); + extra_atom_types = force->inumeric(FLERR,arg[iarg+1]); + if (extra_atom_types < 0) error->all(FLERR,"Illegal read_data command"); + iarg += 2; + } else if (strcmp(arg[iarg],"extra/bond/types") == 0) { + if (iarg+2 > narg) error->all(FLERR,"Illegal read_data command"); + if (!atom->avec->bonds_allow) + error->all(FLERR,"No bonds allowed with this atom style"); + extra_bond_types = force->inumeric(FLERR,arg[iarg+1]); + if (extra_bond_types < 0) error->all(FLERR,"Illegal read_data command"); + iarg += 2; + } else if (strcmp(arg[iarg],"extra/angle/types") == 0) { + if (iarg+2 > narg) error->all(FLERR,"Illegal read_data command"); + if (!atom->avec->angles_allow) + error->all(FLERR,"No angles allowed with this atom style"); + extra_angle_types = force->inumeric(FLERR,arg[iarg+1]); + if (extra_angle_types < 0) error->all(FLERR,"Illegal read_data command"); + iarg += 2; + } else if (strcmp(arg[iarg],"extra/dihedral/types") == 0) { + if (iarg+2 > narg) error->all(FLERR,"Illegal read_data command"); + if (!atom->avec->dihedrals_allow) + error->all(FLERR,"No dihedrals allowed with this atom style"); + extra_dihedral_types = force->inumeric(FLERR,arg[iarg+1]); + if (extra_dihedral_types < 0) + error->all(FLERR,"Illegal read_data command"); + iarg += 2; + } else if (strcmp(arg[iarg],"extra/improper/types") == 0) { + if (iarg+2 > narg) error->all(FLERR,"Illegal read_data command"); + if (!atom->avec->impropers_allow) + error->all(FLERR,"No impropers allowed with this atom style"); + extra_improper_types = force->inumeric(FLERR,arg[iarg+1]); + if (extra_improper_types < 0) + error->all(FLERR,"Illegal read_data command"); + iarg += 2; + + } else if (strcmp(arg[iarg],"group") == 0) { + if (iarg+2 > narg) error->all(FLERR,"Illegal read_data command"); + int igroup = group->find_or_create(arg[iarg+1]); + groupbit = group->bitmask[igroup]; + iarg += 2; + } else if (strcmp(arg[iarg],"fix") == 0) { if (iarg+4 > narg) error->all(FLERR,"Illegal read_data command"); @@ -161,25 +236,54 @@ void ReadData::command(int narg, char **arg) strcpy(fix_section[nfix],arg[iarg+3]); nfix++; iarg += 4; + } else error->all(FLERR,"Illegal read_data command"); } // error checks - if (domain->box_exist && !addflag && !mergeflag) - error->all(FLERR,"Cannot read_data after simulation box is defined"); - if (addflag && mergeflag) error->all(FLERR,"Cannot read_data add and merge"); - if (domain->dimension == 2 && offset[2] != 0.0) - error->all(FLERR,"Cannot use non-zero z offset in read_data " - "for 2d simulation"); - if (domain->dimension == 2 && domain->zperiodic == 0) error->all(FLERR,"Cannot run 2d simulation with nonperiodic Z dimension"); + if (domain->box_exist && !addflag) + error->all(FLERR,"Cannot read_data without add keyword " + "after simulation box is defined"); + if (!domain->box_exist && addflag) + error->all(FLERR,"Cannot use read_data add before " + "simulation box is defined"); + if (offsetflag && addflag == NONE) + error->all(FLERR,"Cannot use read_data offset without add flag"); + if (shiftflag && addflag == NONE) + error->all(FLERR,"Cannot use read_data shift without add flag"); + if (addflag != NONE && + (extra_atom_types || extra_bond_types || extra_angle_types || + extra_dihedral_types || extra_improper_types)) + error->all(FLERR,"Cannot use read_data extra with add flag"); + + // first time system initialization + + if (addflag == NONE) { + domain->box_exist = 1; + update->ntimestep = 0; + } + + // compute atomID offset for addflag = MERGE + + if (addflag == APPEND) { + tagint *tag = atom->tag; + int nlocal = atom->nlocal; + tagint max = 0; + for (int i = 0; i < nlocal; i++) max = MAX(max,tag[i]); + MPI_Allreduce(&max,&id_offset,1,MPI_LMP_TAGINT,MPI_MAX,world); + } + + // ----------------------------------------------------------------- // perform 1-pass read if no molecular topoogy in file // perform 2-pass read if molecular topology, // first pass calculates max topology/atom + // flags for this data file + int atomflag,topoflag; int bondflag,angleflag,dihedralflag,improperflag; int ellipsoidflag,lineflag,triflag,bodyflag; @@ -188,6 +292,14 @@ void ReadData::command(int narg, char **arg) bondflag = angleflag = dihedralflag = improperflag = 0; ellipsoidflag = lineflag = triflag = bodyflag = 0; + // values in this data file + + natoms = ntypes = 0; + nbonds = nangles = ndihedrals = nimpropers = 0; + nbondtypes = nangletypes = ndihedraltypes = nimpropertypes = 0; + triclinic = 0; + + int nlocal_previous = atom->nlocal; int firstpass = 1; while (1) { @@ -204,16 +316,12 @@ void ReadData::command(int narg, char **arg) header(); // problem setup using info from header - // 1st pass only - - if (firstpass) { - domain->box_exist = 1; - update->ntimestep = 0; - - // apply extra settings before grow(), even if no topology in file - // deallocate() insures new settings are used for topology arrays - // if per-atom topology is in file, another grow() is done below + // only done once, if firstpass and first data file + // apply extra settings before grow(), even if no topology in file + // deallocate() insures new settings are used for topology arrays + // if per-atom topology is in file, another grow() is done below + if (firstpass && addflag == NONE) { atom->bond_per_atom = atom->extra_bond_per_atom; atom->angle_per_atom = atom->extra_angle_per_atom; atom->dihedral_per_atom = atom->extra_dihedral_per_atom; @@ -227,6 +335,37 @@ void ReadData::command(int narg, char **arg) atom->deallocate_topology(); atom->avec->grow(n); + domain->boxlo[0] = boxlo[0]; domain->boxhi[0] = boxhi[0]; + domain->boxlo[1] = boxlo[1]; domain->boxhi[1] = boxhi[1]; + domain->boxlo[2] = boxlo[2]; domain->boxhi[2] = boxhi[2]; + + if (triclinic) { + domain->xy = xy; domain->xz = xz; domain->yz = yz; + } + + domain->print_box(" "); + domain->set_initial_box(); + domain->set_global_box(); + comm->set_proc_grid(); + domain->set_local_box(); + } + + // change simulation box to be union of existing box and new box + shift + // only done if firstpass and not first data file + + if (firstpass && addflag != NONE) { + domain->boxlo[0] = MIN(domain->boxlo[0],boxlo[0]+shift[0]); + domain->boxhi[0] = MAX(domain->boxhi[0],boxhi[0]+shift[0]); + domain->boxlo[1] = MIN(domain->boxlo[1],boxlo[1]+shift[1]); + domain->boxhi[1] = MAX(domain->boxhi[1],boxhi[1]+shift[1]); + domain->boxlo[2] = MIN(domain->boxlo[2],boxlo[2]+shift[2]); + domain->boxhi[2] = MAX(domain->boxhi[2],boxhi[2]+shift[2]); + + // NOTE: not sure what to do about this: + //if (triclinic) { + // domain->xy = xy; domain->xz = xz; domain->yz = yz; + // } + domain->print_box(" "); domain->set_initial_box(); domain->set_global_box(); @@ -261,34 +400,34 @@ void ReadData::command(int narg, char **arg) error->warning(FLERR,"Atom style in data file differs " "from currently defined atom style"); atoms(); - } else skip_lines(atom->natoms); + } else skip_lines(natoms); } else if (strcmp(keyword,"Velocities") == 0) { if (atomflag == 0) error->all(FLERR,"Must read Atoms before Velocities"); if (firstpass) velocities(); - else skip_lines(atom->natoms); + else skip_lines(natoms); } else if (strcmp(keyword,"Bonds") == 0) { topoflag = bondflag = 1; - if (atom->nbonds == 0) + if (nbonds == 0) error->all(FLERR,"Invalid data file section: Bonds"); if (atomflag == 0) error->all(FLERR,"Must read Atoms before Bonds"); bonds(firstpass); } else if (strcmp(keyword,"Angles") == 0) { topoflag = angleflag = 1; - if (atom->nangles == 0) + if (nangles == 0) error->all(FLERR,"Invalid data file section: Angles"); if (atomflag == 0) error->all(FLERR,"Must read Atoms before Angles"); angles(firstpass); } else if (strcmp(keyword,"Dihedrals") == 0) { topoflag = dihedralflag = 1; - if (atom->ndihedrals == 0) + if (ndihedrals == 0) error->all(FLERR,"Invalid data file section: Dihedrals"); if (atomflag == 0) error->all(FLERR,"Must read Atoms before Dihedrals"); dihedrals(firstpass); } else if (strcmp(keyword,"Impropers") == 0) { topoflag = improperflag = 1; - if (atom->nimpropers == 0) + if (nimpropers == 0) error->all(FLERR,"Invalid data file section: Impropers"); if (atomflag == 0) error->all(FLERR,"Must read Atoms before Impropers"); impropers(firstpass); @@ -325,7 +464,7 @@ void ReadData::command(int narg, char **arg) } else if (strcmp(keyword,"Masses") == 0) { if (firstpass) mass(); - else skip_lines(atom->ntypes); + else skip_lines(ntypes); } else if (strcmp(keyword,"Pair Coeffs") == 0) { if (force->pair == NULL) error->all(FLERR,"Must define pair_style before Pair Coeffs"); @@ -334,7 +473,7 @@ void ReadData::command(int narg, char **arg) error->warning(FLERR,"Pair style in data file differs " "from currently defined pair style"); paircoeffs(); - } else skip_lines(atom->ntypes); + } else skip_lines(ntypes); } else if (strcmp(keyword,"PairIJ Coeffs") == 0) { if (force->pair == NULL) error->all(FLERR,"Must define pair_style before PairIJ Coeffs"); @@ -343,7 +482,7 @@ void ReadData::command(int narg, char **arg) error->warning(FLERR,"Pair style in data file differs " "from currently defined pair style"); pairIJcoeffs(); - } else skip_lines(atom->ntypes*(atom->ntypes+1)/2); + } else skip_lines(ntypes*(ntypes+1)/2); } else if (strcmp(keyword,"Bond Coeffs") == 0) { if (atom->avec->bonds_allow == 0) error->all(FLERR,"Invalid data file section: Bond Coeffs"); @@ -354,7 +493,7 @@ void ReadData::command(int narg, char **arg) error->warning(FLERR,"Bond style in data file differs " "from currently defined bond style"); bondcoeffs(); - } else skip_lines(atom->nbondtypes); + } else skip_lines(nbondtypes); } else if (strcmp(keyword,"Angle Coeffs") == 0) { if (atom->avec->angles_allow == 0) error->all(FLERR,"Invalid data file section: Angle Coeffs"); @@ -365,7 +504,7 @@ void ReadData::command(int narg, char **arg) error->warning(FLERR,"Angle style in data file differs " "from currently defined angle style"); anglecoeffs(0); - } else skip_lines(atom->nangletypes); + } else skip_lines(nangletypes); } else if (strcmp(keyword,"Dihedral Coeffs") == 0) { if (atom->avec->dihedrals_allow == 0) error->all(FLERR,"Invalid data file section: Dihedral Coeffs"); @@ -376,7 +515,7 @@ void ReadData::command(int narg, char **arg) error->warning(FLERR,"Dihedral style in data file differs " "from currently defined dihedral style"); dihedralcoeffs(0); - } else skip_lines(atom->ndihedraltypes); + } else skip_lines(ndihedraltypes); } else if (strcmp(keyword,"Improper Coeffs") == 0) { if (atom->avec->impropers_allow == 0) error->all(FLERR,"Invalid data file section: Improper Coeffs"); @@ -387,7 +526,7 @@ void ReadData::command(int narg, char **arg) error->warning(FLERR,"Improper style in data file differs " "from currently defined improper style"); impropercoeffs(0); - } else skip_lines(atom->nimpropertypes); + } else skip_lines(nimpropertypes); } else if (strcmp(keyword,"BondBond Coeffs") == 0) { if (atom->avec->angles_allow == 0) @@ -395,14 +534,14 @@ void ReadData::command(int narg, char **arg) if (force->angle == NULL) error->all(FLERR,"Must define angle_style before BondBond Coeffs"); if (firstpass) anglecoeffs(1); - else skip_lines(atom->nangletypes); + else skip_lines(nangletypes); } else if (strcmp(keyword,"BondAngle Coeffs") == 0) { if (atom->avec->angles_allow == 0) error->all(FLERR,"Invalid data file section: BondAngle Coeffs"); if (force->angle == NULL) error->all(FLERR,"Must define angle_style before BondAngle Coeffs"); if (firstpass) anglecoeffs(2); - else skip_lines(atom->nangletypes); + else skip_lines(nangletypes); } else if (strcmp(keyword,"MiddleBondTorsion Coeffs") == 0) { if (atom->avec->dihedrals_allow == 0) @@ -413,7 +552,7 @@ void ReadData::command(int narg, char **arg) "Must define dihedral_style before " "MiddleBondTorsion Coeffs"); if (firstpass) dihedralcoeffs(1); - else skip_lines(atom->ndihedraltypes); + else skip_lines(ndihedraltypes); } else if (strcmp(keyword,"EndBondTorsion Coeffs") == 0) { if (atom->avec->dihedrals_allow == 0) error->all(FLERR,"Invalid data file section: EndBondTorsion Coeffs"); @@ -421,7 +560,7 @@ void ReadData::command(int narg, char **arg) error->all(FLERR, "Must define dihedral_style before EndBondTorsion Coeffs"); if (firstpass) dihedralcoeffs(2); - else skip_lines(atom->ndihedraltypes); + else skip_lines(ndihedraltypes); } else if (strcmp(keyword,"AngleTorsion Coeffs") == 0) { if (atom->avec->dihedrals_allow == 0) error->all(FLERR,"Invalid data file section: AngleTorsion Coeffs"); @@ -429,7 +568,7 @@ void ReadData::command(int narg, char **arg) error->all(FLERR, "Must define dihedral_style before AngleTorsion Coeffs"); if (firstpass) dihedralcoeffs(3); - else skip_lines(atom->ndihedraltypes); + else skip_lines(ndihedraltypes); } else if (strcmp(keyword,"AngleAngleTorsion Coeffs") == 0) { if (atom->avec->dihedrals_allow == 0) error->all(FLERR, @@ -439,7 +578,7 @@ void ReadData::command(int narg, char **arg) "Must define dihedral_style before " "AngleAngleTorsion Coeffs"); if (firstpass) dihedralcoeffs(4); - else skip_lines(atom->ndihedraltypes); + else skip_lines(ndihedraltypes); } else if (strcmp(keyword,"BondBond13 Coeffs") == 0) { if (atom->avec->dihedrals_allow == 0) error->all(FLERR,"Invalid data file section: BondBond13 Coeffs"); @@ -447,7 +586,7 @@ void ReadData::command(int narg, char **arg) error->all(FLERR, "Must define dihedral_style before BondBond13 Coeffs"); if (firstpass) dihedralcoeffs(5); - else skip_lines(atom->ndihedraltypes); + else skip_lines(ndihedraltypes); } else if (strcmp(keyword,"AngleAngle Coeffs") == 0) { if (atom->avec->impropers_allow == 0) @@ -456,7 +595,7 @@ void ReadData::command(int narg, char **arg) error->all(FLERR, "Must define improper_style before AngleAngle Coeffs"); if (firstpass) impropercoeffs(1); - else skip_lines(atom->nimpropertypes); + else skip_lines(nimpropertypes); } else { char str[128]; @@ -469,7 +608,7 @@ void ReadData::command(int narg, char **arg) // error if natoms > 0 yet no atoms were read - if (atom->natoms > 0 && atomflag == 0) + if (natoms > 0 && atomflag == 0) error->all(FLERR,"No atoms in data file"); // close file @@ -486,9 +625,8 @@ void ReadData::command(int narg, char **arg) // at end of 1st pass, error check for required sections // customize for new sections - if ((atom->nbonds && !bondflag) || (atom->nangles && !angleflag) || - (atom->ndihedrals && !dihedralflag) || - (atom->nimpropers && !improperflag)) + if ((nbonds && !bondflag) || (nangles && !angleflag) || + (ndihedrals && !dihedralflag) || (nimpropers && !improperflag)) error->one(FLERR,"Needed molecular topology not in data file"); if ((nellipsoids && !ellipsoidflag) || (nlines && !lineflag) || @@ -510,6 +648,15 @@ void ReadData::command(int narg, char **arg) atom->avec->grow(atom->nmax); } + // assign atoms added by this data file to specified group + + if (groupbit) { + int *mask = atom->mask; + int nlocal = atom->nlocal; + for (int i = nlocal_previous; i < nlocal; i++) + mask[i] |= groupbit; + } + // create special bond lists for molecular systems if (atom->molecular == 1) { @@ -590,6 +737,19 @@ void ReadData::command(int narg, char **arg) if (atom->molecular == 2) atom->avec->onemols[0]->check_attributes(1); + // if adding atoms, migrate atoms to new processors + // use irregular() b/c box size could have changed dramaticaly + // resulting in procs now owning very different subboxes + // with their previously owned atoms now far outside the subbox + + if (addflag != NONE) { + if (domain->triclinic) domain->x2lamda(atom->nlocal); + Irregular *irregular = new Irregular(lmp); + irregular->migrate_atoms(1); + delete irregular; + if (domain->triclinic) domain->lamda2x(atom->nlocal); + } + // shrink-wrap the box if necessary and move atoms to new procs // if atoms are lost is b/c data file box was far from shrink-wrapped // do not use irregular() comm, which would not lose atoms, @@ -621,6 +781,7 @@ void ReadData::command(int narg, char **arg) header ends with EOF or non-blank line containing no header keyword if EOF, line is set to blank line else line has first keyword line for rest of file + some logic differs if adding atoms ------------------------------------------------------------------------- */ void ReadData::header() @@ -648,7 +809,7 @@ void ReadData::header() while (1) { - // read a line and bcast length if flag is set + // read a line and bcast length if (me == 0) { if (fgets(line,MAXLINE,fp) == NULL) n = 0; @@ -689,7 +850,9 @@ void ReadData::header() // customize for new header lines if (strstr(line,"atoms")) { - sscanf(line,BIGINT_FORMAT,&atom->natoms); + sscanf(line,BIGINT_FORMAT,&natoms); + if (addflag == NONE) atom->natoms = natoms; + else atom->natoms += natoms; // check for these first // otherwise "triangles" will be matched as "angles" @@ -710,43 +873,70 @@ void ReadData::header() if (!avec_body) error->all(FLERR,"No bodies allowed with this atom style"); sscanf(line,BIGINT_FORMAT,&nbodies); - } - else if (strstr(line,"bonds")) sscanf(line,BIGINT_FORMAT,&atom->nbonds); - else if (strstr(line,"angles")) sscanf(line,BIGINT_FORMAT,&atom->nangles); - else if (strstr(line,"dihedrals")) sscanf(line,BIGINT_FORMAT, - &atom->ndihedrals); - else if (strstr(line,"impropers")) sscanf(line,BIGINT_FORMAT, - &atom->nimpropers); + } else if (strstr(line,"bonds")) { + sscanf(line,BIGINT_FORMAT,&nbonds); + if (addflag == NONE) atom->nbonds = nbonds; + else atom->nbonds += nbonds; + } else if (strstr(line,"angles")) { + sscanf(line,BIGINT_FORMAT,&nangles); + if (addflag == NONE) atom->nangles = nangles; + else atom->nangles += nangles; + } else if (strstr(line,"dihedrals")) { + sscanf(line,BIGINT_FORMAT,&ndihedrals); + if (addflag == NONE) atom->ndihedrals = ndihedrals; + else atom->ndihedrals += natoms; + } else if (strstr(line,"impropers")) { + sscanf(line,BIGINT_FORMAT,&nimpropers); + if (addflag == NONE) atom->nimpropers = nimpropers; + else atom->nimpropers += nimpropers; + + // Atom class type settings are only set by first data file - else if (strstr(line,"atom types")) sscanf(line,"%d",&atom->ntypes); - else if (strstr(line,"bond types")) sscanf(line,"%d",&atom->nbondtypes); - else if (strstr(line,"angle types")) sscanf(line,"%d",&atom->nangletypes); - else if (strstr(line,"dihedral types")) - sscanf(line,"%d",&atom->ndihedraltypes); - else if (strstr(line,"improper types")) - sscanf(line,"%d",&atom->nimpropertypes); + } else if (strstr(line,"atom types")) { + sscanf(line,"%d",&ntypes); + if (addflag == NONE) atom->ntypes = ntypes + extra_atom_types; + } else if (strstr(line,"bond types")) { + sscanf(line,"%d",&nbondtypes); + if (addflag == NONE) atom->nbondtypes = nbondtypes + extra_bond_types; + } else if (strstr(line,"angle types")) { + sscanf(line,"%d",&nangletypes); + if (addflag == NONE) atom->nangletypes = nangletypes + extra_angle_types; + } else if (strstr(line,"dihedral types")) { + sscanf(line,"%d",&ndihedraltypes); + if (addflag == NONE) + atom->ndihedraltypes = ndihedraltypes + extra_dihedral_types; + } else if (strstr(line,"improper types")) { + sscanf(line,"%d",&nimpropertypes); + if (addflag == NONE) + atom->nimpropertypes = nimpropertypes + extra_improper_types; - else if (strstr(line,"extra bond per atom")) - sscanf(line,"%d",&atom->extra_bond_per_atom); - else if (strstr(line,"extra angle per atom")) - sscanf(line,"%d",&atom->extra_angle_per_atom); - else if (strstr(line,"extra dihedral per atom")) - sscanf(line,"%d",&atom->extra_dihedral_per_atom); - else if (strstr(line,"extra improper per atom")) - sscanf(line,"%d",&atom->extra_improper_per_atom); - else if (strstr(line,"extra special per atom")) - sscanf(line,"%d",&force->special_extra); + // these settings only used by first data file + + } else if (strstr(line,"extra bond per atom")) { + if (addflag == NONE) sscanf(line,"%d",&atom->extra_bond_per_atom); + } else if (strstr(line,"extra angle per atom")) { + if (addflag == NONE) sscanf(line,"%d",&atom->extra_angle_per_atom); + } else if (strstr(line,"extra dihedral per atom")) { + if (addflag == NONE) sscanf(line,"%d",&atom->extra_dihedral_per_atom); + } else if (strstr(line,"extra improper per atom")) { + if (addflag == NONE) sscanf(line,"%d",&atom->extra_improper_per_atom); + } else if (strstr(line,"extra special per atom")) { + if (addflag == NONE) sscanf(line,"%d",&force->special_extra); + + // local copy of box info + // so can treat differently for first vs subsequent data files + + } else if (strstr(line,"xlo xhi")) { + sscanf(line,"%lg %lg",&boxlo[0],&boxhi[0]); + } else if (strstr(line,"ylo yhi")) { + sscanf(line,"%lg %lg",&boxlo[1],&boxhi[1]); + } else if (strstr(line,"zlo zhi")) { + sscanf(line,"%lg %lg",&boxlo[2],&boxhi[2]); + } else if (strstr(line,"xy xz yz")) { + triclinic = 1; + sscanf(line,"%lg %lg %lg",&xy,&xz,&yz); - else if (strstr(line,"xlo xhi")) - sscanf(line,"%lg %lg",&domain->boxlo[0],&domain->boxhi[0]); - else if (strstr(line,"ylo yhi")) - sscanf(line,"%lg %lg",&domain->boxlo[1],&domain->boxhi[1]); - else if (strstr(line,"zlo zhi")) - sscanf(line,"%lg %lg",&domain->boxlo[2],&domain->boxhi[2]); - else if (strstr(line,"xy xz yz")) { - domain->triclinic = 1; - sscanf(line,"%lg %lg %lg",&domain->xy,&domain->xz,&domain->yz); } else break; } @@ -815,27 +1005,28 @@ void ReadData::atoms() } bigint nread = 0; - bigint natoms = atom->natoms; while (nread < natoms) { nchunk = MIN(natoms-nread,CHUNK); eof = comm->read_lines_from_file(fp,nchunk,MAXLINE,buffer); if (eof) error->all(FLERR,"Unexpected end of data file"); - atom->data_atoms(nchunk,buffer); + atom->data_atoms(nchunk,buffer,id_offset,toffset,shiftflag,shift); nread += nchunk; } // check that all atoms were assigned correctly - bigint tmp = atom->nlocal; - MPI_Allreduce(&tmp,&natoms,1,MPI_LMP_BIGINT,MPI_SUM,world); + bigint n = atom->nlocal; + bigint sum; + MPI_Allreduce(&n,&sum,1,MPI_LMP_BIGINT,MPI_SUM,world); + bigint nassign = sum - (atom->natoms - natoms); if (me == 0) { - if (screen) fprintf(screen," " BIGINT_FORMAT " atoms\n",natoms); - if (logfile) fprintf(logfile," " BIGINT_FORMAT " atoms\n",natoms); + if (screen) fprintf(screen," " BIGINT_FORMAT " atoms\n",nassign); + if (logfile) fprintf(logfile," " BIGINT_FORMAT " atoms\n",nassign); } - if (natoms != atom->natoms) + if (sum != atom->natoms) error->all(FLERR,"Did not assign all atoms correctly"); // check that atom IDs are valid @@ -872,13 +1063,12 @@ void ReadData::velocities() } bigint nread = 0; - bigint natoms = atom->natoms; while (nread < natoms) { nchunk = MIN(natoms-nread,CHUNK); eof = comm->read_lines_from_file(fp,nchunk,MAXLINE,buffer); if (eof) error->all(FLERR,"Unexpected end of data file"); - atom->data_vels(nchunk,buffer); + atom->data_vels(nchunk,buffer,id_offset); nread += nchunk; } @@ -928,7 +1118,7 @@ void ReadData::bonds(int firstpass) 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); + atom->data_bonds(nchunk,buffer,count,id_offset); nread += nchunk; } @@ -1003,7 +1193,7 @@ void ReadData::angles(int firstpass) 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); + atom->data_angles(nchunk,buffer,count,id_offset); nread += nchunk; } @@ -1078,7 +1268,7 @@ void ReadData::dihedrals(int firstpass) 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); + atom->data_dihedrals(nchunk,buffer,count,id_offset); nread += nchunk; } @@ -1153,7 +1343,7 @@ void ReadData::impropers(int firstpass) 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); + atom->data_impropers(nchunk,buffer,count,id_offset); nread += nchunk; } @@ -1216,7 +1406,7 @@ void ReadData::bonus(bigint nbonus, AtomVec *ptr, const char *type) nchunk = MIN(natoms-nread,CHUNK); eof = comm->read_lines_from_file(fp,nchunk,MAXLINE,buffer); if (eof) error->all(FLERR,"Unexpected end of data file"); - atom->data_bonus(nchunk,buffer,ptr); + atom->data_bonus(nchunk,buffer,ptr,id_offset); nread += nchunk; } @@ -1296,7 +1486,7 @@ void ReadData::bodies(int firstpass) MPI_Bcast(&m,1,MPI_INT,0,world); MPI_Bcast(buffer,m,MPI_CHAR,0,world); - if (firstpass) atom->data_bodies(nchunk,buffer,avec_body); + if (firstpass) atom->data_bodies(nchunk,buffer,avec_body,id_offset); nread += nchunk; } @@ -1316,16 +1506,16 @@ void ReadData::bodies(int firstpass) void ReadData::mass() { char *next; - char *buf = new char[atom->ntypes*MAXLINE]; + char *buf = new char[ntypes*MAXLINE]; - int eof = comm->read_lines_from_file(fp,atom->ntypes,MAXLINE,buf); + int eof = comm->read_lines_from_file(fp,ntypes,MAXLINE,buf); if (eof) error->all(FLERR,"Unexpected end of data file"); char *original = buf; - for (int i = 0; i < atom->ntypes; i++) { + for (int i = 0; i < ntypes; i++) { next = strchr(buf,'\n'); *next = '\0'; - atom->set_mass(buf); + atom->set_mass(buf,toffset); buf = next + 1; } delete [] original; @@ -1336,16 +1526,16 @@ void ReadData::mass() void ReadData::paircoeffs() { char *next; - char *buf = new char[atom->ntypes*MAXLINE]; + char *buf = new char[ntypes*MAXLINE]; - int eof = comm->read_lines_from_file(fp,atom->ntypes,MAXLINE,buf); + int eof = comm->read_lines_from_file(fp,ntypes,MAXLINE,buf); if (eof) error->all(FLERR,"Unexpected end of data file"); char *original = buf; - for (int i = 0; i < atom->ntypes; i++) { + for (int i = 0; i < ntypes; i++) { next = strchr(buf,'\n'); *next = '\0'; - parse_coeffs(buf,NULL,1); + parse_coeffs(buf,NULL,1,2,toffset); if (narg == 0) error->all(FLERR,"Unexpected end of PairCoeffs section"); force->pair->coeff(narg,arg); buf = next + 1; @@ -1360,18 +1550,18 @@ void ReadData::pairIJcoeffs() int i,j; char *next; - int nsq = atom->ntypes* (atom->ntypes+1) / 2; + int nsq = ntypes * (ntypes+1) / 2; char *buf = new char[nsq * MAXLINE]; int eof = comm->read_lines_from_file(fp,nsq,MAXLINE,buf); if (eof) error->all(FLERR,"Unexpected end of data file"); char *original = buf; - for (i = 0; i < atom->ntypes; i++) - for (j = i; j < atom->ntypes; j++) { + for (i = 0; i < ntypes; i++) + for (j = i; j < ntypes; j++) { next = strchr(buf,'\n'); *next = '\0'; - parse_coeffs(buf,NULL,0); + parse_coeffs(buf,NULL,0,2,toffset); if (narg == 0) error->all(FLERR,"Unexpected end of PairCoeffs section"); force->pair->coeff(narg,arg); buf = next + 1; @@ -1384,16 +1574,16 @@ void ReadData::pairIJcoeffs() void ReadData::bondcoeffs() { char *next; - char *buf = new char[atom->nbondtypes*MAXLINE]; + char *buf = new char[nbondtypes*MAXLINE]; - int eof = comm->read_lines_from_file(fp,atom->nbondtypes,MAXLINE,buf); + int eof = comm->read_lines_from_file(fp,nbondtypes,MAXLINE,buf); if (eof) error->all(FLERR,"Unexpected end of data file"); char *original = buf; - for (int i = 0; i < atom->nbondtypes; i++) { + for (int i = 0; i < nbondtypes; i++) { next = strchr(buf,'\n'); *next = '\0'; - parse_coeffs(buf,NULL,0); + parse_coeffs(buf,NULL,0,1,boffset); if (narg == 0) error->all(FLERR,"Unexpected end of BondCoeffs section"); force->bond->coeff(narg,arg); buf = next + 1; @@ -1406,18 +1596,18 @@ void ReadData::bondcoeffs() void ReadData::anglecoeffs(int which) { char *next; - char *buf = new char[atom->nangletypes*MAXLINE]; + char *buf = new char[nangletypes*MAXLINE]; - int eof = comm->read_lines_from_file(fp,atom->nangletypes,MAXLINE,buf); + int eof = comm->read_lines_from_file(fp,nangletypes,MAXLINE,buf); if (eof) error->all(FLERR,"Unexpected end of data file"); char *original = buf; - for (int i = 0; i < atom->nangletypes; i++) { + for (int i = 0; i < nangletypes; i++) { next = strchr(buf,'\n'); *next = '\0'; - if (which == 0) parse_coeffs(buf,NULL,0); - else if (which == 1) parse_coeffs(buf,"bb",0); - else if (which == 2) parse_coeffs(buf,"ba",0); + if (which == 0) parse_coeffs(buf,NULL,0,1,aoffset); + else if (which == 1) parse_coeffs(buf,"bb",0,1,aoffset); + else if (which == 2) parse_coeffs(buf,"ba",0,1,aoffset); if (narg == 0) error->all(FLERR,"Unexpected end of AngleCoeffs section"); force->angle->coeff(narg,arg); buf = next + 1; @@ -1430,21 +1620,21 @@ void ReadData::anglecoeffs(int which) void ReadData::dihedralcoeffs(int which) { char *next; - char *buf = new char[atom->ndihedraltypes*MAXLINE]; + char *buf = new char[ndihedraltypes*MAXLINE]; - int eof = comm->read_lines_from_file(fp,atom->ndihedraltypes,MAXLINE,buf); + int eof = comm->read_lines_from_file(fp,ndihedraltypes,MAXLINE,buf); if (eof) error->all(FLERR,"Unexpected end of data file"); char *original = buf; - for (int i = 0; i < atom->ndihedraltypes; i++) { + for (int i = 0; i < ndihedraltypes; i++) { next = strchr(buf,'\n'); *next = '\0'; - if (which == 0) parse_coeffs(buf,NULL,0); - else if (which == 1) parse_coeffs(buf,"mbt",0); - else if (which == 2) parse_coeffs(buf,"ebt",0); - else if (which == 3) parse_coeffs(buf,"at",0); - else if (which == 4) parse_coeffs(buf,"aat",0); - else if (which == 5) parse_coeffs(buf,"bb13",0); + if (which == 0) parse_coeffs(buf,NULL,0,1,doffset); + else if (which == 1) parse_coeffs(buf,"mbt",0,1,doffset); + else if (which == 2) parse_coeffs(buf,"ebt",0,1,doffset); + else if (which == 3) parse_coeffs(buf,"at",0,1,doffset); + else if (which == 4) parse_coeffs(buf,"aat",0,1,doffset); + else if (which == 5) parse_coeffs(buf,"bb13",0,1,doffset); if (narg == 0) error->all(FLERR,"Unexpected end of DihedralCoeffs section"); force->dihedral->coeff(narg,arg); buf = next + 1; @@ -1457,17 +1647,17 @@ void ReadData::dihedralcoeffs(int which) void ReadData::impropercoeffs(int which) { char *next; - char *buf = new char[atom->nimpropertypes*MAXLINE]; + char *buf = new char[nimpropertypes*MAXLINE]; - int eof = comm->read_lines_from_file(fp,atom->nimpropertypes,MAXLINE,buf); + int eof = comm->read_lines_from_file(fp,nimpropertypes,MAXLINE,buf); if (eof) error->all(FLERR,"Unexpected end of data file"); char *original = buf; - for (int i = 0; i < atom->nimpropertypes; i++) { + for (int i = 0; i < nimpropertypes; i++) { next = strchr(buf,'\n'); *next = '\0'; - if (which == 0) parse_coeffs(buf,NULL,0); - else if (which == 1) parse_coeffs(buf,"aa",0); + if (which == 0) parse_coeffs(buf,NULL,0,1,ioffset); + else if (which == 1) parse_coeffs(buf,"aa",0,1,ioffset); if (narg == 0) error->all(FLERR,"Unexpected end of ImproperCoeffs section"); force->improper->coeff(narg,arg); buf = next + 1; @@ -1634,9 +1824,11 @@ void ReadData::skip_lines(bigint n) if 2nd word starts with letter, then is hybrid style, add addstr after it else add addstr before 2nd word if dupflag, duplicate 1st word, so pair_coeff "2" becomes "2 2" + if noffset, add offset to first noffset args, which are atom/bond/etc types ------------------------------------------------------------------------- */ -void ReadData::parse_coeffs(char *line, const char *addstr, int dupflag) +void ReadData::parse_coeffs(char *line, const char *addstr, + int dupflag, int noffset, int offset) { char *ptr; if ((ptr = strchr(line,'#'))) *ptr = '\0'; @@ -1655,6 +1847,17 @@ void ReadData::parse_coeffs(char *line, const char *addstr, int dupflag) if (dupflag && narg == 1) arg[narg++] = word; word = strtok(NULL," \t\n\r\f"); } + + if (noffset) { + int value = force->inumeric(FLERR,arg[0]); + sprintf(argoffset1,"%d",value+offset); + arg[0] = argoffset1; + if (noffset == 2) { + value = force->inumeric(FLERR,arg[1]); + sprintf(argoffset2,"%d",value+offset); + arg[1] = argoffset2; + } + } } /* ---------------------------------------------------------------------- diff --git a/src/read_data.h b/src/read_data.h index 938e9da4ba..74c3e2c8b9 100644 --- a/src/read_data.h +++ b/src/read_data.h @@ -32,19 +32,19 @@ class ReadData : protected Pointers { void command(int, char **); private: + int me,compressed; char *line,*keyword,*buffer,*style; FILE *fp; char **arg; - int me,narg,maxarg,compressed; + int narg,maxarg; + char argoffset1[8],argoffset2[8]; - // optional args + bigint id_offset; - int addflag,mergeflag; - double offset[3]; - int nfix; - int *fix_index; - char **fix_header; - char **fix_section; + bigint natoms; + bigint nbonds,nangles,ndihedrals,nimpropers; + int ntypes; + int nbondtypes,nangletypes,ndihedraltypes,nimpropertypes; bigint nellipsoids; class AtomVecEllipsoid *avec_ellipsoid; @@ -55,13 +55,36 @@ class ReadData : protected Pointers { bigint nbodies; class AtomVecBody *avec_body; + // box info + + double boxlo[3],boxhi[3]; + double xy,xz,yz; + int triclinic; + + // optional args + + int addflag,offsetflag,shiftflag; + tagint addvalue; + int toffset,boffset,aoffset,doffset,ioffset; + double shift[3]; + int extra_atom_types,extra_bond_types,extra_angle_types; + int extra_dihedral_types,extra_improper_types; + int groupbit; + + int nfix; + int *fix_index; + char **fix_header; + char **fix_section; + + // methods + void open(char *); void scan(int &, int &, int &, int &); int reallocate(int **, int, int); void header(); void parse_keyword(int); void skip_lines(bigint); - void parse_coeffs(char *, const char *, int); + void parse_coeffs(char *, const char *, int, int, int); int style_match(const char *, const char *); void atoms();