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

Merge pull request #2450 from jrgissing/molecule-custom_id_order

Browse Source 
molecule: use user-provided IDs in molecule files
This commit is contained in:
Richard Berger
2020-11-18 17:52:16 -05:00
committed by GitHub
parent 61ccccf908 145d688fa4
commit 96db39f08b
2 changed files with 69 additions and 35 deletions
Download Patch File Download Diff File Expand all files Collapse all files

11
doc/src/molecule.rst
View File

@ -201,11 +201,12 @@ bonds between nuclear cores and Drude electrons in a different manner.
Each section is listed below in alphabetic order. The format of each
section is described including the number of lines it must contain and
rules (if any) for whether it can appear in the data file. In each
case the ID is ignored; it is simply included for readability, and
should be a number from 1 to Nlines for the section, indicating which
atom (or bond, etc) the entry applies to. The lines are assumed to be
listed in order from 1 to Nlines, but LAMMPS does not check for this.
rules (if any) for whether it can appear in the data file. For per-
atom sections, entries should be numbered from 1 to Natoms (where
Natoms is the number of atoms in the template), indicating which atom
(or bond, etc) the entry applies to. Per-atom sections need to
include a setting for every atom, but the atoms can be listed in any
order.
----------

93
src/molecule.cpp
View File

@ -509,7 +509,6 @@ void Molecule::read(int flag)
// count = vector for tallying bonds,angles,etc per atom
if (flag == 0) memory->create(count,natoms,"molecule:count");
else count = nullptr;
// grab keyword and skip next line
@ -616,10 +615,6 @@ void Molecule::read(int flag)
parse_keyword(1,line,keyword);
}
// clean up
memory->destroy(count);
// error check
if (flag == 0) {
@ -665,6 +660,10 @@ void Molecule::read(int flag)
maxradius = radius[0];
}
}
// clean up
if (flag) memory->destroy(count);
}
/* ----------------------------------------------------------------------
@ -673,6 +672,7 @@ void Molecule::read(int flag)
void Molecule::coords(char *line)
{
for (int i = 0; i < natoms; i++) count[i] = 0;
try {
for (int i = 0; i < natoms; i++) {
readline(line);
@ -680,20 +680,25 @@ void Molecule::coords(char *line)
ValueTokenizer values(line);
if (values.count() != 4) error->one(FLERR,"Invalid Coords section in molecule file");
values.next_int();
x[i][0] = values.next_double();
x[i][1] = values.next_double();
x[i][2] = values.next_double();
int iatom = values.next_int() - 1;
if (iatom < 0 || iatom >= natoms) error->one(FLERR,"Invalid Coords section in molecule file");
count[iatom]++;
x[iatom][0] = values.next_double();
x[iatom][1] = values.next_double();
x[iatom][2] = values.next_double();
x[i][0] *= sizescale;
x[i][1] *= sizescale;
x[i][2] *= sizescale;
x[iatom][0] *= sizescale;
x[iatom][1] *= sizescale;
x[iatom][2] *= sizescale;
}
} catch (TokenizerException &e) {
error->one(FLERR, fmt::format("Invalid Coords section in molecule file\n"
"{}", e.what()));
}
for (int i = 0; i < natoms; i++)
if (count[i] == 0) error->all(FLERR,"Invalid Coords section in molecule file");
if (domain->dimension == 2) {
for (int i = 0; i < natoms; i++)
if (x[i][2] != 0.0)
@ -708,6 +713,7 @@ void Molecule::coords(char *line)
void Molecule::types(char *line)
{
for (int i = 0; i < natoms; i++) count[i] = 0;
try {
for (int i = 0; i < natoms; i++) {
readline(line);
@ -715,15 +721,20 @@ void Molecule::types(char *line)
ValueTokenizer values(line);
if (values.count() != 2) error->one(FLERR,"Invalid Types section in molecule file");
values.next_int();
type[i] = values.next_int();
type[i] += toffset;
int iatom = values.next_int() - 1;
if (iatom < 0 || iatom >= natoms) error->one(FLERR,"Invalid Types section in molecule file");
count[iatom]++;
type[iatom] = values.next_int();
type[iatom] += toffset;
}
} catch (TokenizerException &e) {
error->one(FLERR, fmt::format("Invalid Types section in molecule file\n"
"{}", e.what()));
}
for (int i = 0; i < natoms; i++)
if (count[i] == 0) error->all(FLERR,"Invalid Types section in molecule file");
for (int i = 0; i < natoms; i++)
if ((type[i] <= 0) || (domain->box_exist && (type[i] > atom->ntypes)))
error->all(FLERR,"Invalid atom type in molecule file");
@ -739,21 +750,27 @@ void Molecule::types(char *line)
void Molecule::molecules(char *line)
{
for (int i = 0; i < natoms; i++) count[i] = 0;
try {
for (int i = 0; i < natoms; i++) {
readline(line);
ValueTokenizer values(line);
if (values.count() != 2) error->one(FLERR,"Invalid Molecules section in molecule file");
values.next_int();
molecule[i] = values.next_int();
// molecule[i] += moffset; // placeholder for possible molecule offset
int iatom = values.next_int() - 1;
if (iatom < 0 || iatom >= natoms) error->one(FLERR,"Invalid Molecules section in molecule file");
count[iatom]++;
molecule[iatom] = values.next_int();
// molecule[iatom] += moffset; // placeholder for possible molecule offset
}
} catch (TokenizerException &e) {
error->one(FLERR, fmt::format("Invalid Molecules section in molecule file\n"
"{}", e.what()));
}
for (int i = 0; i < natoms; i++)
if (count[i] == 0) error->all(FLERR,"Invalid Molecules section in molecule file");
for (int i = 0; i < natoms; i++)
if (molecule[i] <= 0)
error->all(FLERR,"Invalid molecule ID in molecule file");
@ -798,6 +815,7 @@ void Molecule::fragments(char *line)
void Molecule::charges(char *line)
{
for (int i = 0; i < natoms; i++) count[i] = 0;
try {
for (int i = 0; i < natoms; i++) {
readline(line);
@ -805,13 +823,18 @@ void Molecule::charges(char *line)
ValueTokenizer values(line);
if ((int)values.count() != 2) error->one(FLERR,"Invalid Charges section in molecule file");
values.next_int();
q[i] = values.next_double();
int iatom = values.next_int() - 1;
if (iatom < 0 || iatom >= natoms) error->one(FLERR,"Invalid Charges section in molecule file");
count[iatom]++;
q[iatom] = values.next_double();
}
} catch (TokenizerException &e) {
error->one(FLERR, fmt::format("Invalid Charges section in molecule file\n"
"{}", e.what()));
}
for (int i = 0; i < natoms; i++)
if (count[i] == 0) error->all(FLERR,"Invalid Charges section in molecule file");
}
/* ----------------------------------------------------------------------
@ -820,6 +843,7 @@ void Molecule::charges(char *line)
void Molecule::diameters(char *line)
{
for (int i = 0; i < natoms; i++) count[i] = 0;
try {
maxradius = 0.0;
for (int i = 0; i < natoms; i++) {
@ -828,17 +852,22 @@ void Molecule::diameters(char *line)
ValueTokenizer values(line);
if (values.count() != 2) error->one(FLERR,"Invalid Diameters section in molecule file");
values.next_int();
radius[i] = values.next_double();
radius[i] *= sizescale;
radius[i] *= 0.5;
maxradius = MAX(maxradius,radius[i]);
int iatom = values.next_int() - 1;
if (iatom < 0 || iatom >= natoms) error->one(FLERR,"Invalid Diameters section in molecule file");
count[iatom]++;
radius[iatom] = values.next_double();
radius[iatom] *= sizescale;
radius[iatom] *= 0.5;
maxradius = MAX(maxradius,radius[iatom]);
}
} catch (TokenizerException &e) {
error->one(FLERR, fmt::format("Invalid Diameters section in molecule file\n"
"{}", e.what()));
}
for (int i = 0; i < natoms; i++)
if (count[i] == 0) error->all(FLERR,"Invalid Diameters section in molecule file");
for (int i = 0; i < natoms; i++)
if (radius[i] < 0.0)
error->all(FLERR,"Invalid atom diameter in molecule file");
@ -850,6 +879,7 @@ void Molecule::diameters(char *line)
void Molecule::masses(char *line)
{
for (int i = 0; i < natoms; i++) count[i] = 0;
try {
for (int i = 0; i < natoms; i++) {
readline(line);
@ -857,15 +887,20 @@ void Molecule::masses(char *line)
ValueTokenizer values(line);
if (values.count() != 2) error->one(FLERR,"Invalid Masses section in molecule file");
values.next_int();
rmass[i] = values.next_double();
rmass[i] *= sizescale*sizescale*sizescale;
int iatom = values.next_int() - 1;
if (iatom < 0 || iatom >= natoms) error->one(FLERR,"Invalid Masses section in molecule file");
count[iatom]++;
rmass[iatom] = values.next_double();
rmass[iatom] *= sizescale*sizescale*sizescale;
}
} catch (TokenizerException &e) {
error->one(FLERR, fmt::format("Invalid Masses section in molecule file\n"
"{}", e.what()));
}
for (int i = 0; i < natoms; i++)
if (count[i] == 0) error->all(FLERR,"Invalid Masses section in molecule file");
for (int i = 0; i < natoms; i++)
if (rmass[i] <= 0.0) error->all(FLERR,"Invalid atom mass in molecule file");
}
@ -1296,7 +1331,6 @@ void Molecule::special_generate()
{
int newton_bond = force->newton_bond;
tagint atom1,atom2;
int *count = new int[natoms];
// temporary array for special atoms
@ -1382,7 +1416,6 @@ void Molecule::special_generate()
}
}
}
delete[] count;
maxspecial = 0;
for (int i = 0; i < natoms; i++)