ICODE-ADC/lammps
4
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

git-svn-id: svn://svn.icms.temple.edu/lammps-ro/trunk@13643 f3b2605a-c512-4ea7-a41b-209d697bcdaa

Browse Source
This commit is contained in:
sjplimp
2015-07-17 23:20:48 +00:00
parent 6d665f1860
commit 1ebb654b4d
6 changed files with 475 additions and 180 deletions
Download Patch File Download Diff File Expand all files Collapse all files

94
src/atom.cpp
View File

@ -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");

18
src/atom.h
View File

@ -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 *);

20
src/group.cpp
View File

@ -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

1
src/group.h
View File

@ -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 *);

475
src/read_data.cpp
View File

@ -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;
// 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;
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);
// Atom class type settings are only set by first data file
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);
} 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;
// 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;
}
}
}
/* ----------------------------------------------------------------------

41
src/read_data.h
View File

@ -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();