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

proper handle partial initialization from infile

Browse Source
This commit is contained in:
Axel Kohlmeyer
2022-03-20 20:23:55 -04:00
parent 89c7d8f707
commit 550ae15dff
20 changed files with 3228 additions and 323 deletions
Download Patch File Download Diff File Expand all files Collapse all files

157
src/RIGID/fix_rigid.cpp
View File

@ -134,11 +134,9 @@ FixRigid::FixRigid(LAMMPS *lmp, int narg, char **arg) :
int is_double,cols;
int custom_index = atom->find_custom(arg[4]+2,is_double,cols);
if (custom_index == -1)
error->all(FLERR,"Fix rigid custom requires "
"previously defined property/atom");
error->all(FLERR,"Fix rigid custom requires previously defined property/atom");
else if (is_double)
error->all(FLERR,"Fix rigid custom requires "
"integer-valued property/atom vector");
error->all(FLERR,"Fix rigid custom requires integer-valued property/atom vector");
int minval = INT_MAX;
int *value = atom->ivector[custom_index];
for (i = 0; i < nlocal; i++)
@ -155,9 +153,9 @@ FixRigid::FixRigid(LAMMPS *lmp, int narg, char **arg) :
} else if (utils::strmatch(arg[4],"^v_")) {
int ivariable = input->variable->find(arg[4]+2);
if (ivariable < 0)
error->all(FLERR,"Variable name for fix rigid custom does not exist");
error->all(FLERR,"Variable {} for fix rigid/small custom does not exist", arg[4]+2);
if (input->variable->atomstyle(ivariable) == 0)
error->all(FLERR,"Fix rigid custom variable is no atom-style variable");
error->all(FLERR,"Fix rigid custom variable {} is not atom-style variable", arg[4]+2);
double *value = new double[nlocal];
input->variable->compute_atom(ivariable,0,value,1,0);
int minval = INT_MAX;
@ -218,7 +216,7 @@ FixRigid::FixRigid(LAMMPS *lmp, int narg, char **arg) :
}
memory->destroy(ncount);
if (custom_flag) delete [] molecule;
if (custom_flag) delete[] molecule;
// each listed group is a rigid body
// check if all listed groups exist
@ -259,7 +257,7 @@ FixRigid::FixRigid(LAMMPS *lmp, int narg, char **arg) :
if (flagall)
error->all(FLERR,"One or more atoms belong to multiple rigid bodies");
delete [] igroups;
delete[] igroups;
} else error->all(FLERR,"Illegal fix rigid command");
@ -502,7 +500,7 @@ FixRigid::FixRigid(LAMMPS *lmp, int narg, char **arg) :
if (strcmp(arg[iarg+1],"all") == 0) allremap = 1;
else {
allremap = 0;
delete [] id_dilate;
delete[] id_dilate;
id_dilate = utils::strdup(arg[iarg+1]);
int idilate = group->find(id_dilate);
if (idilate == -1)
@ -529,7 +527,7 @@ FixRigid::FixRigid(LAMMPS *lmp, int narg, char **arg) :
} else if (strcmp(arg[iarg],"infile") == 0) {
if (iarg+2 > narg) error->all(FLERR,"Illegal fix rigid command");
delete [] inpfile;
delete[] inpfile;
inpfile = utils::strdup(arg[iarg+1]);
restart_file = 1;
reinitflag = 0;
@ -542,7 +540,7 @@ FixRigid::FixRigid(LAMMPS *lmp, int narg, char **arg) :
} else if (strcmp(arg[iarg],"gravity") == 0) {
if (iarg+2 > narg) error->all(FLERR,"Illegal fix rigid command");
delete [] id_gravity;
delete[] id_gravity;
id_gravity = utils::strdup(arg[iarg+1]);
iarg += 2;
@ -584,7 +582,7 @@ FixRigid::FixRigid(LAMMPS *lmp, int narg, char **arg) :
if (body[i] >= 0) ncount[body[i]]++;
MPI_Allreduce(ncount,nrigid,nbody,MPI_INT,MPI_SUM,world);
delete [] ncount;
delete[] ncount;
for (ibody = 0; ibody < nbody; ibody++)
if (nrigid[ibody] <= 1) error->all(FLERR,"One or zero atoms in rigid body");
@ -615,9 +613,9 @@ FixRigid::~FixRigid()
atom->delete_callback(id,Atom::GROW);
delete random;
delete [] inpfile;
delete [] id_dilate;
delete [] id_gravity;
delete[] inpfile;
delete[] id_dilate;
delete[] id_gravity;
memory->destroy(mol2body);
memory->destroy(body2mol);
@ -674,8 +672,6 @@ int FixRigid::setmask()
void FixRigid::init()
{
int i,ibody;
triclinic = domain->triclinic;
// atom style pointers to particles that store extra info
@ -688,18 +684,16 @@ void FixRigid::init()
// if earlyflag, warn if any post-force fixes come after a rigid fix
int count = 0;
for (i = 0; i < modify->nfix; i++)
if (modify->fix[i]->rigid_flag) count++;
for (auto ifix : modify->get_fix_list())
if (ifix->rigid_flag) count++;
if (count > 1 && me == 0) error->warning(FLERR,"More than one fix rigid");
if (earlyflag) {
int rflag = 0;
for (i = 0; i < modify->nfix; i++) {
if (modify->fix[i]->rigid_flag) rflag = 1;
if (rflag && (modify->fmask[i] & POST_FORCE) &&
!modify->fix[i]->rigid_flag)
error->warning(FLERR,"Fix {} alters forces after fix rigid",
modify->fix[i]->id);
bool rflag = false;
for (auto ifix : modify->get_fix_list()) {
if (ifix->rigid_flag) rflag = true;
if ((comm->me == 0) && rflag && (ifix->setmask() & POST_FORCE) && !ifix->rigid_flag)
error->warning(FLERR,"Fix {} alters forces after fix rigid", ifix->id);
}
}
@ -709,36 +703,30 @@ void FixRigid::init()
// and gravity is not applied correctly
if (inpfile && !id_gravity) {
for (i = 0; i < modify->nfix; i++) {
if (strcmp(modify->fix[i]->style,"gravity") == 0) {
if (comm->me == 0)
error->warning(FLERR,"Gravity may not be correctly applied "
"to rigid bodies if they consist of "
"overlapped particles");
break;
}
}
if (modify->get_fix_by_style("^gravity").size() > 0)
if (comm->me == 0)
error->warning(FLERR,"Gravity may not be correctly applied to rigid "
"bodies if they consist of overlapped particles");
}
// error if a fix changing the box comes before rigid fix
for (i = 0; i < modify->nfix; i++)
if (modify->fix[i]->box_change) break;
if (i < modify->nfix) {
for (int j = i+1; j < modify->nfix; j++)
if (utils::strmatch(modify->fix[j]->style,"^rigid"))
bool boxflag = false;
for (auto ifix : modify->get_fix_list()) {
if (ifix->box_change) boxflag = true;
if (boxflag && utils::strmatch(ifix->style,"^rigid"))
error->all(FLERR,"Rigid fixes must come before any box changing fix");
}
// add gravity forces based on gravity vector from fix
if (id_gravity) {
int ifix = modify->find_fix(id_gravity);
if (ifix < 0) error->all(FLERR,"Fix rigid cannot find fix gravity ID");
if (!utils::strmatch(modify->fix[ifix]->style,"^gravity"))
error->all(FLERR,"Fix rigid gravity fix ID is not a gravity fix style");
auto ifix = modify->get_fix_by_id(id_gravity);
if (!ifix) error->all(FLERR,"Fix rigid cannot find fix gravity ID {}", id_gravity);
if (!utils::strmatch(ifix->style,"^gravity"))
error->all(FLERR,"Fix rigid gravity fix ID {} is not a gravity fix style", id_gravity);
int tmp;
gvec = (double *) modify->fix[ifix]->extract("gvec",tmp);
gvec = (double *) ifix->extract("gvec", tmp);
}
// timestep info
@ -767,7 +755,7 @@ void FixRigid::init()
// temperature scale factor
double ndof = 0.0;
for (ibody = 0; ibody < nbody; ibody++) {
for (int ibody = 0; ibody < nbody; ibody++) {
ndof += fflag[ibody][0] + fflag[ibody][1] + fflag[ibody][2];
ndof += tflag[ibody][0] + tflag[ibody][1] + tflag[ibody][2];
}
@ -1237,7 +1225,7 @@ int FixRigid::dof(int tgroup)
if (!setupflag) {
if (comm->me == 0)
error->warning(FLERR,"Cannot count rigid body degrees-of-freedom "
"before bodies are initialized");
"before bodies are fully initialized");
return 0;
}
@ -1307,10 +1295,10 @@ int FixRigid::dof(int tgroup)
n += 2*nall[ibody] + 3*mall[ibody] - 3;
}
delete [] ncount;
delete [] mcount;
delete [] nall;
delete [] mall;
delete[] ncount;
delete[] mcount;
delete[] nall;
delete[] mall;
return n;
}
@ -1397,12 +1385,9 @@ void FixRigid::set_xv()
MathExtra::matvec(ex_space[ibody],ey_space[ibody],
ez_space[ibody],displace[i],x[i]);
v[i][0] = omega[ibody][1]*x[i][2] - omega[ibody][2]*x[i][1] +
vcm[ibody][0];
v[i][1] = omega[ibody][2]*x[i][0] - omega[ibody][0]*x[i][2] +
vcm[ibody][1];
v[i][2] = omega[ibody][0]*x[i][1] - omega[ibody][1]*x[i][0] +
vcm[ibody][2];
v[i][0] = omega[ibody][1]*x[i][2] - omega[ibody][2]*x[i][1] + vcm[ibody][0];
v[i][1] = omega[ibody][2]*x[i][0] - omega[ibody][0]*x[i][2] + vcm[ibody][1];
v[i][2] = omega[ibody][0]*x[i][1] - omega[ibody][1]*x[i][0] + vcm[ibody][2];
// add center of mass to displacement
// map back into periodic box via xbox,ybox,zbox
@ -1555,12 +1540,9 @@ void FixRigid::set_v()
v2 = v[i][2];
}
v[i][0] = omega[ibody][1]*delta[2] - omega[ibody][2]*delta[1] +
vcm[ibody][0];
v[i][1] = omega[ibody][2]*delta[0] - omega[ibody][0]*delta[2] +
vcm[ibody][1];
v[i][2] = omega[ibody][0]*delta[1] - omega[ibody][1]*delta[0] +
vcm[ibody][2];
v[i][0] = omega[ibody][1]*delta[2] - omega[ibody][2]*delta[1] + vcm[ibody][0];
v[i][1] = omega[ibody][2]*delta[0] - omega[ibody][0]*delta[2] + vcm[ibody][1];
v[i][2] = omega[ibody][0]*delta[1] - omega[ibody][1]*delta[0] + vcm[ibody][2];
// virial = unwrapped coords dotted into body constraint force
// body constraint force = implied force due to v change minus f external
@ -1658,8 +1640,7 @@ void FixRigid::setup_bodies_static()
{
int i,ibody;
// extended = 1 if any particle in a rigid body is finite size
// or has a dipole moment
// extended = 1 if any particle in a rigid body is finite size or has a dipole moment
extended = orientflag = dorientflag = 0;
@ -1822,6 +1803,11 @@ void FixRigid::setup_bodies_static()
int *inbody;
if (inpfile) {
// must call it here so it doesn't override read in data but
// initialize bodies whose dynamic settings not set in inpfile
setup_bodies_dynamic();
memory->create(inbody,nbody,"rigid:inbody");
for (ibody = 0; ibody < nbody; ibody++) inbody[ibody] = 0;
readfile(0,masstotal,xcm,vcm,angmom,imagebody,inbody);
@ -2268,8 +2254,7 @@ void FixRigid::setup_bodies_dynamic()
vxcm vycm vzcm lx ly lz ix iy iz
------------------------------------------------------------------------- */
void FixRigid::readfile(int which, double *vec,
double **array1, double **array2, double **array3,
void FixRigid::readfile(int which, double *vec, double **array1, double **array2, double **array3,
imageint *ivec, int *inbody)
{
int nchunk,id,eofflag,xbox,ybox,zbox;
@ -2278,30 +2263,36 @@ void FixRigid::readfile(int which, double *vec,
char *eof,*start,*next,*buf;
char line[MAXLINE];
// open file and read and parse first non-empty, non-comment line containing the number of bodies
if (me == 0) {
fp = fopen(inpfile,"r");
if (fp == nullptr)
error->one(FLERR,"Cannot open fix rigid file {}: {}", inpfile,utils::getsyserror());
error->one(FLERR,"Cannot open fix rigid infile {}: {}", inpfile, utils::getsyserror());
while (true) {
eof = fgets(line,MAXLINE,fp);
if (eof == nullptr) error->one(FLERR,"Unexpected end of fix rigid file");
if (eof == nullptr) error->one(FLERR,"Unexpected end of fix rigid infile");
start = &line[strspn(line," \t\n\v\f\r")];
if (*start != '\0' && *start != '#') break;
}
nlines = utils::inumeric(FLERR, utils::trim(line), true, lmp);
if (which == 0)
utils::logmesg(lmp, "Reading rigid body data for {} bodies from file {}\n", nlines, inpfile);
if (nlines == 0) fclose(fp);
}
MPI_Bcast(&nlines,1,MPI_INT,0,world);
// empty file with 0 lines is needed to trigger initial restart file generation
if (nlines == 0) return; //error->all(FLERR,"Fix rigid file has no lines");
// empty file with 0 lines is needed to trigger initial restart file
// generation when no infile was previously used.
if (nlines == 0) return;
else if (nlines < 0) error->all(FLERR,"Fix rigid infile has incorrect format");
char *buffer = new char[CHUNK*MAXLINE];
int nread = 0;
while (nread < nlines) {
nchunk = MIN(nlines-nread,CHUNK);
eofflag = utils::read_lines_from_file(fp,nchunk,MAXLINE,buffer,me,world);
if (eofflag) error->all(FLERR,"Unexpected end of fix rigid file");
if (eofflag) error->all(FLERR,"Unexpected end of fix rigid infile");
buf = buffer;
next = strchr(buf,'\n');
@ -2328,7 +2319,7 @@ void FixRigid::readfile(int which, double *vec,
id = values.next_int();
if (rstyle == MOLECULE) {
if (id <= 0 || id > maxmol)
throw TokenizerException("invalid rigid body ID ", std::to_string(id));
throw TokenizerException("invalid rigid molecule ID ", std::to_string(id));
id = mol2body[id];
} else id--;
@ -2365,7 +2356,7 @@ void FixRigid::readfile(int which, double *vec,
array1[id][3] = values.next_double();
}
} catch (TokenizerException &e) {
error->all(FLERR, "Invalid fix rigid file: {}", e.what());
error->all(FLERR, "Invalid fix rigid infile: {}", e.what());
}
buf = next + 1;
}
@ -2384,7 +2375,7 @@ void FixRigid::readfile(int which, double *vec,
void FixRigid::write_restart_file(const char *file)
{
if (me) return;
if (comm->me) return;
auto outfile = std::string(file) + ".rigid";
FILE *fp = fopen(outfile.c_str(),"w");
@ -2414,14 +2405,11 @@ void FixRigid::write_restart_file(const char *file)
ybox = (imagebody[i] >> IMGBITS & IMGMASK) - IMGMAX;
zbox = (imagebody[i] >> IMG2BITS) - IMGMAX;
fprintf(fp,"%d %.16g %.16g %.16g %.16g %.16g %.16g %.16g %.16g %.16g %.16g "
"%.16g %.16g %.16g %.16g %.16g %.16g %d %d %d\n",
id,masstotal[i],xcm[i][0],xcm[i][1],xcm[i][2],
ispace[0][0],ispace[1][1],ispace[2][2],
ispace[0][1],ispace[0][2],ispace[1][2],
vcm[i][0],vcm[i][1],vcm[i][2],
angmom[i][0],angmom[i][1],angmom[i][2],
xbox,ybox,zbox);
fprintf(fp,"%d %-1.16e %-1.16e %-1.16e %-1.16e %-1.16e %-1.16e %-1.16e %-1.16e %-1.16e "
"%-1.16e %-1.16e %-1.16e %-1.16e %-1.16e %-1.16e %-1.16e %d %d %d\n",
id,masstotal[i],xcm[i][0],xcm[i][1],xcm[i][2],ispace[0][0],ispace[1][1],ispace[2][2],
ispace[0][1],ispace[0][2],ispace[1][2],vcm[i][0],vcm[i][1],vcm[i][2],
angmom[i][0],angmom[i][1],angmom[i][2],xbox,ybox,zbox);
}
fclose(fp);
@ -2483,6 +2471,7 @@ void FixRigid::copy_arrays(int i, int j, int /*delflag*/)
displace[j][0] = displace[i][0];
displace[j][1] = displace[i][1];
displace[j][2] = displace[i][2];
if (extended) {
eflags[j] = eflags[i];
for (int k = 0; k < orientflag; k++)