From 21833c99de7a09095850635a9936474c4ccdd544 Mon Sep 17 00:00:00 2001 From: sjplimp Date: Wed, 8 Jan 2014 20:38:38 +0000 Subject: [PATCH] git-svn-id: svn://svn.icms.temple.edu/lammps-ro/trunk@11192 f3b2605a-c512-4ea7-a41b-209d697bcdaa --- src/GRANULAR/fix_pour.cpp | 6 +- src/MISC/fix_deposit.cpp | 4 +- src/atom.cpp | 11 +- src/create_atoms.cpp | 287 +++++++++++++++++++++++++++++++------- src/create_atoms.h | 6 +- src/create_box.cpp | 20 +-- src/molecule.cpp | 5 +- src/read_data.cpp | 2 + 8 files changed, 272 insertions(+), 69 deletions(-) diff --git a/src/GRANULAR/fix_pour.cpp b/src/GRANULAR/fix_pour.cpp index c4003f5301..97b6f72427 100644 --- a/src/GRANULAR/fix_pour.cpp +++ b/src/GRANULAR/fix_pour.cpp @@ -118,7 +118,7 @@ FixPour::FixPour(LAMMPS *lmp, int narg, char **arg) : if (mode == MOLECULE) { if (atom->molecule_flag == 0) - error->all(FLERR,"Fix pour requires atom attribute molecule"); + error->all(FLERR,"Fix pour mol requires atom attribute molecule"); if (onemol->xflag == 0) error->all(FLERR,"Fix pour molecule must have coordinates"); if (onemol->typeflag == 0) @@ -587,7 +587,7 @@ void FixPour::pre_exchange() else atom->avec->create_atom(ntype+onemol->type[m],coords[m]); int n = atom->nlocal - 1; atom->tag[n] = maxtag_all + m+1; - if (mode == MOLECULE) atom->molecule[n] = maxmol_all; + if (mode == MOLECULE) atom->molecule[n] = maxmol_all+1; atom->mask[n] = 1 | groupbit; atom->image[n] = imageflags[m]; atom->v[n][0] = vnew[0]; @@ -836,7 +836,7 @@ void FixPour::options(int narg, char **arg) error->all(FLERR,"Molecule ID for fix pour does not exist"); mode = MOLECULE; onemol = atom->molecules[imol]; - iarg += 2; + iarg += 2; } else if (strcmp(arg[iarg],"rigid") == 0) { if (iarg+2 > narg) error->all(FLERR,"Illegal fix pour command"); int n = strlen(arg[iarg+1]) + 1; diff --git a/src/MISC/fix_deposit.cpp b/src/MISC/fix_deposit.cpp index e5e59acc53..78b075200a 100644 --- a/src/MISC/fix_deposit.cpp +++ b/src/MISC/fix_deposit.cpp @@ -100,7 +100,7 @@ FixDeposit::FixDeposit(LAMMPS *lmp, int narg, char **arg) : if (mode == MOLECULE) { if (atom->molecule_flag == 0) - error->all(FLERR,"Fix deposit requires atom attribute molecule"); + error->all(FLERR,"Fix deposit mol requires atom attribute molecule"); if (onemol->xflag == 0) error->all(FLERR,"Fix deposit molecule must have coordinates"); if (onemol->typeflag == 0) @@ -444,7 +444,7 @@ void FixDeposit::pre_exchange() else atom->avec->create_atom(ntype+onemol->type[m],coords[m]); n = atom->nlocal - 1; atom->tag[n] = maxtag_all + m+1; - if (mode == MOLECULE) atom->molecule[n] = maxmol_all; + if (mode == MOLECULE) atom->molecule[n] = maxmol_all+1; atom->mask[n] = 1 | groupbit; atom->image[n] = imageflags[m]; atom->v[n][0] = vnew[0]; diff --git a/src/atom.cpp b/src/atom.cpp index 714449769c..d4127963ac 100644 --- a/src/atom.cpp +++ b/src/atom.cpp @@ -1205,7 +1205,8 @@ int Atom::find_molecule(char *id) } /* ---------------------------------------------------------------------- - add info for iatom from molecule template onemol to current atom ilocal + add info to current atom ilocal from molecule template onemol and its iatom + offset = atom ID preceeding IDs of atoms in this molecule ------------------------------------------------------------------------- */ void Atom::add_molecule_atom(Molecule *onemol, int iatom, @@ -1255,7 +1256,15 @@ void Atom::add_molecule_atom(Molecule *onemol, int iatom, } } + // error check against maxspecial in case user has not done one of these: + // create_box extra/special/per/atom N + // read_data extra special per atom N + // special_bonds extra N + // if explicitly used special_bonds, may not have maintained extra + if (onemol->specialflag) { + if (onemol->maxspecial > maxspecial) + error->one(FLERR,"Molecule file special bond counts are too large"); nspecial[ilocal][0] = onemol->nspecial[iatom][0]; nspecial[ilocal][1] = onemol->nspecial[iatom][1]; int n = nspecial[ilocal][2] = onemol->nspecial[iatom][2]; diff --git a/src/create_atoms.cpp b/src/create_atoms.cpp index 55a0216e3d..2c86525cc7 100644 --- a/src/create_atoms.cpp +++ b/src/create_atoms.cpp @@ -17,22 +17,29 @@ #include "create_atoms.h" #include "atom.h" #include "atom_vec.h" +#include "molecule.h" #include "comm.h" +#include "irregular.h" #include "modify.h" +#include "force.h" #include "fix.h" #include "domain.h" #include "lattice.h" #include "region.h" #include "random_park.h" +#include "random_mars.h" +#include "math_extra.h" +#include "math_const.h" #include "error.h" -#include "force.h" using namespace LAMMPS_NS; +using namespace MathConst; #define BIG 1.0e30 #define EPSILON 1.0e-6 enum{BOX,REGION,SINGLE,RANDOM}; +enum{ATOM,MOLECULE}; /* ---------------------------------------------------------------------- */ @@ -51,9 +58,7 @@ void CreateAtoms::command(int narg, char **arg) // parse arguments if (narg < 2) error->all(FLERR,"Illegal create_atoms command"); - itype = force->inumeric(FLERR,arg[0]); - if (itype <= 0 || itype > atom->ntypes) - error->all(FLERR,"Invalid atom type in create_atoms command"); + ntype = force->inumeric(FLERR,arg[0]); int iarg; if (strcmp(arg[1],"box") == 0) { @@ -93,18 +98,19 @@ void CreateAtoms::command(int narg, char **arg) int scaleflag = 1; remapflag = 0; + mode = ATOM; + int molseed; nbasis = domain->lattice->nbasis; basistype = new int[nbasis]; - for (int i = 0; i < nbasis; i++) basistype[i] = itype; + for (int i = 0; i < nbasis; i++) basistype[i] = ntype; while (iarg < narg) { if (strcmp(arg[iarg],"basis") == 0) { if (iarg+3 > narg) error->all(FLERR,"Illegal create_atoms command"); int ibasis = force->inumeric(FLERR,arg[iarg+1]); - itype = force->inumeric(FLERR,arg[iarg+2]); - if (ibasis <= 0 || ibasis > nbasis || - itype <= 0 || itype > atom->ntypes) + int itype = force->inumeric(FLERR,arg[iarg+2]); + if (ibasis <= 0 || ibasis > nbasis || itype <= 0 || itype > atom->ntypes) error->all(FLERR,"Invalid basis setting in create_atoms command"); basistype[ibasis-1] = itype; iarg += 3; @@ -114,6 +120,15 @@ void CreateAtoms::command(int narg, char **arg) else if (strcmp(arg[iarg+1],"no") == 0) remapflag = 0; else error->all(FLERR,"Illegal create_atoms command"); iarg += 2; + } else if (strcmp(arg[iarg],"mol") == 0) { + if (iarg+3 > narg) error->all(FLERR,"Illegal fix create_atoms command"); + int imol = atom->find_molecule(arg[iarg+1]); + if (imol == -1) + error->all(FLERR,"Molecule ID for create_atoms does not exist"); + mode = MOLECULE; + onemol = atom->molecules[imol]; + molseed = force->inumeric(FLERR,arg[iarg+2]); + iarg += 3; } else if (strcmp(arg[iarg],"units") == 0) { if (iarg+2 > narg) error->all(FLERR,"Illegal create_atoms command"); if (strcmp(arg[iarg+1],"box") == 0) scaleflag = 0; @@ -125,11 +140,36 @@ void CreateAtoms::command(int narg, char **arg) // error checks + if (mode == ATOM && (ntype <= 0 || ntype > atom->ntypes)) + error->all(FLERR,"Invalid atom type in create_atoms command"); + if (style == RANDOM) { if (nrandom < 0) error->all(FLERR,"Illegal create_atoms command"); if (seed <= 0) error->all(FLERR,"Illegal create_atoms command"); } + // error check and further setup for mode = MOLECULE + + ranmol = NULL; + if (mode == MOLECULE) { + if (atom->molecule_flag == 0) + error->all(FLERR,"Create_atoms mol requires atom attribute molecule"); + if (onemol->xflag == 0) + error->all(FLERR,"Create_atoms molecule must have coordinates"); + if (onemol->typeflag == 0) + error->all(FLERR,"Create_atoms molecule must have atom types"); + if (ntype+onemol->maxtype <= 0 || ntype+onemol->maxtype > atom->ntypes) + error->all(FLERR,"Invalid atom type in create_atoms mol command"); + + // create_atoms uses geoemetric center of molecule for insertion + + onemol->compute_center(); + + // molecule random number generator, same for all procs + + ranmol = new RanMars(lmp,molseed+comm->me); + } + // demand non-none lattice be defined for BOX and REGION // else setup scaling for SINGLE and RANDOM // could use domain->lattice->lattice2box() to do conversion of @@ -192,7 +232,7 @@ void CreateAtoms::command(int narg, char **arg) } } - // add atoms in one of 3 ways + // add atoms/molecules in one of 3 ways bigint natoms_previous = atom->natoms; int nlocal_previous = atom->nlocal; @@ -211,16 +251,148 @@ void CreateAtoms::command(int narg, char **arg) fix->set_arrays(i); } - // clean up - - if (domain->lattice) delete [] basistype; - - // new total # of atoms + // new total # of atoms and error check + // for MOLECULE mode, require atom IDs bigint nblocal = atom->nlocal; MPI_Allreduce(&nblocal,&atom->natoms,1,MPI_LMP_BIGINT,MPI_SUM,world); if (atom->natoms < 0 || atom->natoms > MAXBIGINT) error->all(FLERR,"Too many total atoms"); + if (atom->natoms > MAXSMALLINT) { + if (mode == ATOM) { + if (comm->me == 0) + error->warning(FLERR,"Total atom count exceeds ID limit, " + "atoms will not have individual IDs"); + atom->tag_enable = 0; + } else error->all(FLERR,"Total atom count exceeds ID limit"); + } + + // for ATOM mode: + // add IDs for newly created atoms if IDs are still enabled + // if global map exists, reset it + + if (mode == ATOM) { + if (atom->natoms <= MAXSMALLINT) atom->tag_extend(); + if (atom->map_style) { + atom->nghost = 0; + atom->map_init(); + atom->map_set(); + } + } + + // for MOLECULE mode: + // set atom and molecule IDs for created atoms + // send atoms to new owning procs via irregular comm + // since not all created atoms will be within my sub-domain + + if (mode == MOLECULE) { + + // add atom IDs for newly created atoms and reset global map + // global map must exist since atom->molecular is required to be set + + atom->tag_extend(); + atom->nghost = 0; + atom->map_init(); + atom->map_set(); + + // maxmol = max molecule ID across all procs, for previous atoms + + int *molecule = atom->molecule; + + int max = 0; + for (int i = 0; i < nlocal_previous; i++) max = MAX(max,molecule[i]); + int maxmol; + MPI_Allreduce(&max,&maxmol,1,MPI_INT,MPI_MAX,world); + + // molcreate = # of molecules I created + // moloffset = max molecule ID for all molecules owned by previous procs + // including molecules existing before this creation + + int molcreate = (atom->nlocal - nlocal_previous) / onemol->natoms; + int moloffset; + MPI_Scan(&molcreate,&moloffset,1,MPI_INT,MPI_SUM,world); + moloffset = moloffset - molcreate + maxmol; + + // loop over molecules I created + // set their molecule ID + // reset their bond,angle,etc and special values + + int natoms = onemol->natoms; + + int *tag = atom->tag; + int *num_bond = atom->num_bond; + int *num_angle = atom->num_angle; + int *num_dihedral = atom->num_dihedral; + int *num_improper = atom->num_improper; + int **bond_atom = atom->bond_atom; + int **angle_atom1 = atom->angle_atom1; + int **angle_atom2 = atom->angle_atom2; + int **angle_atom3 = atom->angle_atom3; + int **dihedral_atom1 = atom->dihedral_atom1; + int **dihedral_atom2 = atom->dihedral_atom2; + int **dihedral_atom3 = atom->dihedral_atom3; + int **dihedral_atom4 = atom->dihedral_atom4; + int **improper_atom1 = atom->improper_atom1; + int **improper_atom2 = atom->improper_atom2; + int **improper_atom3 = atom->improper_atom3; + int **improper_atom4 = atom->improper_atom4; + int **nspecial = atom->nspecial; + int **special = atom->special; + + int ilocal = nlocal_previous; + for (int i = 0; i < molcreate; i++) { + int offset = tag[ilocal]-1; + for (int m = 0; m < natoms; m++) { + molecule[ilocal] = moloffset + i+1; + if (onemol->bondflag) + for (int j = 0; j < num_bond[ilocal]; j++) + bond_atom[ilocal][j] += offset; + if (onemol->angleflag) + for (int j = 0; j < num_angle[ilocal]; j++) { + angle_atom1[ilocal][j] += offset; + angle_atom2[ilocal][j] += offset; + angle_atom3[ilocal][j] += offset; + } + if (onemol->dihedralflag) + for (int j = 0; j < num_dihedral[ilocal]; j++) { + dihedral_atom1[ilocal][j] += offset; + dihedral_atom2[ilocal][j] += offset; + dihedral_atom3[ilocal][j] += offset; + dihedral_atom4[ilocal][j] += offset; + } + if (onemol->improperflag) + for (int j = 0; j < num_improper[ilocal]; j++) { + improper_atom1[ilocal][j] += offset; + improper_atom2[ilocal][j] += offset; + improper_atom3[ilocal][j] += offset; + improper_atom4[ilocal][j] += offset; + } + if (onemol->specialflag) + for (int j = 0; j < nspecial[ilocal][2]; j++) + special[ilocal][j] += offset; + ilocal++; + } + } + + // perform irregular comm to migrate atoms to new owning procs + + double **x = atom->x; + tagint *image = atom->image; + int nlocal = atom->nlocal; + for (int i = 0; i < nlocal; i++) domain->remap(x[i],image[i]); + + if (domain->triclinic) domain->x2lamda(atom->nlocal); + domain->reset_box(); + Irregular *irregular = new Irregular(lmp); + irregular->migrate_atoms(); + delete irregular; + if (domain->triclinic) domain->lamda2x(atom->nlocal); + } + + // clean up + + delete ranmol; + if (domain->lattice) delete [] basistype; // print status @@ -232,37 +404,6 @@ void CreateAtoms::command(int narg, char **arg) fprintf(logfile,"Created " BIGINT_FORMAT " atoms\n", atom->natoms-natoms_previous); } - - // reset simulation now that more atoms are defined - // add tags for newly created atoms if possible - // if global map exists, reset it - // change these to MAXTAGINT when allow tagint = bigint - - if (atom->natoms > MAXSMALLINT) { - if (comm->me == 0) - error->warning(FLERR,"Total atom count exceeds ID limit, " - "atoms will not have individual IDs"); - atom->tag_enable = 0; - } - if (atom->natoms <= MAXSMALLINT) atom->tag_extend(); - - if (atom->map_style) { - atom->nghost = 0; - atom->map_init(); - atom->map_set(); - } - - // if a molecular system, set nspecial to 0 for new atoms - // NOTE: 31May12, don't think this is needed, avec->create_atom() does it - - //if (atom->molecular) { - // int **nspecial = atom->nspecial; - // for (int i = nlocal_previous; i < atom->nlocal; i++) { - // nspecial[i][0] = 0; - // nspecial[i][1] = 0; - // nspecial[i][2] = 0; - // } - //} } /* ---------------------------------------------------------------------- @@ -292,8 +433,10 @@ void CreateAtoms::add_single() 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]) - atom->avec->create_atom(itype,xone); + coord[2] >= sublo[2] && coord[2] < subhi[2]) { + if (mode == ATOM) atom->avec->create_atom(ntype,xone); + else add_molecule(xone); + } } @@ -372,8 +515,10 @@ void CreateAtoms::add_random() 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]) - atom->avec->create_atom(itype,xone); + coord[2] >= sublo[2] && coord[2] < subhi[2]) { + if (mode == ATOM) atom->avec->create_atom(ntype,xone); + else add_molecule(xone); + } } // clean-up @@ -445,7 +590,7 @@ void CreateAtoms::add_lattice() // iterate on 3d periodic lattice of unit cells using loop bounds // iterate on nbasis atoms in each unit cell // convert lattice coords to box coords - // add atom if it meets all criteria + // add atom or molecule (on each basis point) if it meets all criteria double **basis = domain->lattice->basis; double x[3],lamda[3]; @@ -481,8 +626,46 @@ void CreateAtoms::add_lattice() coord[1] < sublo[1] || coord[1] >= subhi[1] || coord[2] < sublo[2] || coord[2] >= subhi[2]) continue; - // add the atom to my list of atoms + // add the atom or entire molecule to my list of atoms - atom->avec->create_atom(basistype[m],x); + if (mode == ATOM) atom->avec->create_atom(basistype[m],x); + else add_molecule(x); } } + +/* ---------------------------------------------------------------------- + add a randomly rotated molecule with its center at center +------------------------------------------------------------------------- */ + +void CreateAtoms::add_molecule(double *center) +{ + int n; + double r[3],rotmat[3][3],quat[4],xnew[3]; + + if (domain->dimension == 3) { + r[0] = ranmol->uniform() - 0.5; + r[1] = ranmol->uniform() - 0.5; + r[2] = ranmol->uniform() - 0.5; + } else { + r[0] = r[1] = 0.0; + r[2] = 1.0; + } + double theta = ranmol->uniform() * MY_2PI; + MathExtra::norm3(r); + MathExtra::axisangle_to_quat(r,theta,quat); + MathExtra::quat_to_mat(quat,rotmat); + + // create atoms in molecule with atom ID = 0 and mol ID = 0 + // reset in caller after all moleclues created by all procs + // pass add_molecule_atom an offset of 0 since don't know + // max tag of atoms in previous molecules at this point + + int natoms = onemol->natoms; + for (int m = 0; m < natoms; m++) { + MathExtra::matvec(rotmat,onemol->dx[m],xnew); + MathExtra::add3(xnew,center,xnew); + atom->avec->create_atom(ntype+onemol->type[m],xnew); + n = atom->nlocal - 1; + atom->add_molecule_atom(onemol,m,n,0); + } +} diff --git a/src/create_atoms.h b/src/create_atoms.h index cc40f63624..f337ff6f23 100644 --- a/src/create_atoms.h +++ b/src/create_atoms.h @@ -30,17 +30,21 @@ class CreateAtoms : protected Pointers { void command(int, char **); private: - int itype,style,nregion,nbasis,nrandom,seed; + int ntype,style,mode,nregion,nbasis,nrandom,seed; int *basistype; double xone[3]; int remapflag; + class Molecule *onemol; + class RanMars *ranmol; + int triclinic; double sublo[3],subhi[3]; // epsilon-extended proc sub-box for adding atoms void add_single(); void add_random(); void add_lattice(); + void add_molecule(double *); }; } diff --git a/src/create_box.cpp b/src/create_box.cpp index c0ed91150d..1e7ce40122 100644 --- a/src/create_box.cpp +++ b/src/create_box.cpp @@ -108,54 +108,58 @@ void CreateBox::command(int narg, char **arg) int iarg = 2; while (iarg < narg) { - if (strcmp(arg[iarg],"bond types") == 0) { + if (strcmp(arg[iarg],"bond/types") == 0) { if (iarg+2 > narg) error->all(FLERR,"Illegal create_box command"); if (!atom->avec->bonds_allow) error->all(FLERR,"No bonds allowed with this atom style"); atom->nbondtypes = force->inumeric(FLERR,arg[iarg+1]); iarg += 2; - } else if (strcmp(arg[iarg],"angle types") == 0) { + } else if (strcmp(arg[iarg],"angle/types") == 0) { if (iarg+2 > narg) error->all(FLERR,"Illegal create_box command"); if (!atom->avec->angles_allow) error->all(FLERR,"No angles allowed with this atom style"); atom->nangletypes = force->inumeric(FLERR,arg[iarg+1]); iarg += 2; - } else if (strcmp(arg[iarg],"dihedral types") == 0) { + } else if (strcmp(arg[iarg],"dihedral/types") == 0) { if (iarg+2 > narg) error->all(FLERR,"Illegal create_box command"); if (!atom->avec->dihedrals_allow) error->all(FLERR,"No dihedrals allowed with this atom style"); atom->ndihedraltypes = force->inumeric(FLERR,arg[iarg+1]); iarg += 2; - } else if (strcmp(arg[iarg],"improper types") == 0) { + } else if (strcmp(arg[iarg],"improper/types") == 0) { if (iarg+2 > narg) error->all(FLERR,"Illegal create_box command"); if (!atom->avec->impropers_allow) error->all(FLERR,"No impropers allowed with this atom style"); atom->nimpropertypes = force->inumeric(FLERR,arg[iarg+1]); iarg += 2; - } else if (strcmp(arg[iarg],"extra bond per atom") == 0) { + } else if (strcmp(arg[iarg],"extra/bond/per/atom") == 0) { if (iarg+2 > narg) error->all(FLERR,"Illegal create_box command"); if (!atom->avec->bonds_allow) error->all(FLERR,"No bonds allowed with this atom style"); atom->bond_per_atom = force->inumeric(FLERR,arg[iarg+1]); iarg += 2; - } else if (strcmp(arg[iarg],"extra angle per atom") == 0) { + } else if (strcmp(arg[iarg],"extra/angle/per/atom") == 0) { if (iarg+2 > narg) error->all(FLERR,"Illegal create_box command"); if (!atom->avec->angles_allow) error->all(FLERR,"No angles allowed with this atom style"); atom->angle_per_atom = force->inumeric(FLERR,arg[iarg+1]); iarg += 2; - } else if (strcmp(arg[iarg],"extra dihedral per atom") == 0) { + } else if (strcmp(arg[iarg],"extra/dihedral/per/atom") == 0) { if (iarg+2 > narg) error->all(FLERR,"Illegal create_box command"); if (!atom->avec->dihedrals_allow) error->all(FLERR,"No dihedrals allowed with this atom style"); atom->dihedral_per_atom = force->inumeric(FLERR,arg[iarg+1]); iarg += 2; - } else if (strcmp(arg[iarg],"extra improper per atom") == 0) { + } else if (strcmp(arg[iarg],"extra/improper/per/atom") == 0) { if (iarg+2 > narg) error->all(FLERR,"Illegal create_box command"); if (!atom->avec->impropers_allow) error->all(FLERR,"No impropers allowed with this atom style"); atom->improper_per_atom = force->inumeric(FLERR,arg[iarg+1]); iarg += 2; + } else if (strcmp(arg[iarg],"extra/special/per/atom") == 0) { + if (iarg+2 > narg) error->all(FLERR,"Illegal create_box command"); + atom->maxspecial = force->inumeric(FLERR,arg[iarg+1]); + iarg += 2; } else error->all(FLERR,"Illegal create_box command"); } diff --git a/src/molecule.cpp b/src/molecule.cpp index 5ec4b9ca63..669a34ce13 100644 --- a/src/molecule.cpp +++ b/src/molecule.cpp @@ -482,12 +482,13 @@ void Molecule::read(int flag) error->all(FLERR,"Molecule file bond/angle/etc counts " "per atom are too large"); + // test for maxspecial > atom->maxspecial is done when molecules added + // in Atom::add_molecule_atom() + if ((nspecialflag && !specialflag) || (!nspecialflag && specialflag)) error->all(FLERR,"Molecule file needs both Special Bond sections"); if (specialflag && !bondflag) error->all(FLERR,"Molecule file has special flags but no bonds"); - if (maxspecial > atom->maxspecial) - error->all(FLERR,"Molecule file special bond counts are too large"); if ((shakeflagflag || shakeatomflag || shaketypeflag) && !shakeflag) error->all(FLERR,"Molecule file shake info is incomplete"); diff --git a/src/read_data.cpp b/src/read_data.cpp index 741e96557d..52f7f50820 100644 --- a/src/read_data.cpp +++ b/src/read_data.cpp @@ -508,6 +508,8 @@ void ReadData::header(int flag) 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); else if (strstr(line,"xlo xhi")) sscanf(line,"%lg %lg",&domain->boxlo[0],&domain->boxhi[0]);