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

add more strict checking of data when parsing molecule files to detect format errors

Browse Source
This commit is contained in:
Axel Kohlmeyer
2017-08-10 14:49:51 -04:00
parent f68c6254d4
commit 2a7d2dee36
1 changed files with 213 additions and 220 deletions
Download Patch File Download Diff File Expand all files Collapse all files

433
src/molecule.cpp
View File

@ -427,47 +427,61 @@ void Molecule::read(int flag)
// search line for header keywords and set corresponding variable // search line for header keywords and set corresponding variable
if (strstr(line,"atoms")) sscanf(line,"%d",&natoms); int nmatch = 0;
else if (strstr(line,"bonds")) sscanf(line,"%d",&nbonds); int nwant = 0;
else if (strstr(line,"angles")) sscanf(line,"%d",&nangles); if (strstr(line,"atoms")) {
else if (strstr(line,"dihedrals")) sscanf(line,"%d",&ndihedrals); nmatch = sscanf(line,"%d",&natoms);
else if (strstr(line,"impropers")) sscanf(line,"%d",&nimpropers); nwant = 1;
} else if (strstr(line,"bonds")) {
else if (strstr(line,"mass")) { nmatch = sscanf(line,"%d",&nbonds);
nwant = 1;
} else if (strstr(line,"angles")) {
nmatch = sscanf(line,"%d",&nangles);
nwant = 1;
} else if (strstr(line,"dihedrals")) {
nmatch = sscanf(line,"%d",&ndihedrals);
nwant = 1;
} else if (strstr(line,"impropers")) {
nmatch = sscanf(line,"%d",&nimpropers);
nwant = 1;
} else if (strstr(line,"mass")) {
massflag = 1; massflag = 1;
sscanf(line,"%lg",&masstotal); nmatch = sscanf(line,"%lg",&masstotal);
nwant = 1;
masstotal *= sizescale*sizescale*sizescale; masstotal *= sizescale*sizescale*sizescale;
} } else if (strstr(line,"com")) {
else if (strstr(line,"com")) {
comflag = 1; comflag = 1;
sscanf(line,"%lg %lg %lg",&com[0],&com[1],&com[2]); nmatch = sscanf(line,"%lg %lg %lg",&com[0],&com[1],&com[2]);
nwant = 3;
com[0] *= sizescale; com[0] *= sizescale;
com[1] *= sizescale; com[1] *= sizescale;
com[2] *= sizescale; com[2] *= sizescale;
if (domain->dimension == 2 && com[2] != 0.0) if (domain->dimension == 2 && com[2] != 0.0)
error->all(FLERR,"Molecule file z center-of-mass must be 0.0 for 2d"); error->all(FLERR,"Molecule file z center-of-mass must be 0.0 for 2d");
} } else if (strstr(line,"inertia")) {
else if (strstr(line,"inertia")) {
inertiaflag = 1; inertiaflag = 1;
sscanf(line,"%lg %lg %lg %lg %lg %lg", nmatch = sscanf(line,"%lg %lg %lg %lg %lg %lg",
&itensor[0],&itensor[1],&itensor[2], &itensor[0],&itensor[1],&itensor[2],
&itensor[3],&itensor[4],&itensor[5]); &itensor[3],&itensor[4],&itensor[5]);
itensor[0] *= sizescale*sizescale*sizescale*sizescale*sizescale; nwant = 6;
itensor[1] *= sizescale*sizescale*sizescale*sizescale*sizescale; const double scale5 = sizescale*sizescale*sizescale*sizescale*sizescale;
itensor[2] *= sizescale*sizescale*sizescale*sizescale*sizescale; itensor[0] *= scale5;
itensor[3] *= sizescale*sizescale*sizescale*sizescale*sizescale; itensor[1] *= scale5;
itensor[4] *= sizescale*sizescale*sizescale*sizescale*sizescale; itensor[2] *= scale5;
itensor[5] *= sizescale*sizescale*sizescale*sizescale*sizescale; itensor[3] *= scale5;
} itensor[4] *= scale5;
else if (strstr(line,"body")) { itensor[5] *= scale5;
} else if (strstr(line,"body")) {
bodyflag = 1; bodyflag = 1;
avec_body = (AtomVecBody *) atom->style_match("body"); avec_body = (AtomVecBody *) atom->style_match("body");
if (!avec_body) if (!avec_body)
error->all(FLERR,"Molecule file requires atom style body"); error->all(FLERR,"Molecule file requires atom style body");
sscanf(line,"%d %d",&nibody,&ndbody); nmatch = sscanf(line,"%d %d",&nibody,&ndbody);
} nwant = 2;
} else break;
else break; if (nmatch != nwant)
error->all(FLERR,"Invalid header in molecule file");
} }
// error checks // error checks
@ -493,7 +507,7 @@ void Molecule::read(int flag)
// loop over sections of molecule file // loop over sections of molecule file
while (strlen(keyword)) { while (strlen(keyword) > 0) {
if (strcmp(keyword,"Coords") == 0) { if (strcmp(keyword,"Coords") == 0) {
xflag = 1; xflag = 1;
if (flag) coords(line); if (flag) coords(line);
@ -517,22 +531,22 @@ void Molecule::read(int flag)
} else if (strcmp(keyword,"Bonds") == 0) { } else if (strcmp(keyword,"Bonds") == 0) {
if (nbonds == 0) if (nbonds == 0)
error->all(FLERR,"Molecule file has bonds but no nbonds setting"); error->all(FLERR,"Molecule file has bonds but no nbonds setting");
bondflag = tag_require = 1; bondflag = tag_require = 1;
bonds(flag,line); bonds(flag,line);
} else if (strcmp(keyword,"Angles") == 0) { } else if (strcmp(keyword,"Angles") == 0) {
if (nangles == 0) if (nangles == 0)
error->all(FLERR,"Molecule file has angles but no nangles setting"); error->all(FLERR,"Molecule file has angles but no nangles setting");
angleflag = tag_require = 1; angleflag = tag_require = 1;
angles(flag,line); angles(flag,line);
} else if (strcmp(keyword,"Dihedrals") == 0) { } else if (strcmp(keyword,"Dihedrals") == 0) {
if (ndihedrals == 0) error->all(FLERR,"Molecule file has dihedrals " if (ndihedrals == 0) error->all(FLERR,"Molecule file has dihedrals "
"but no ndihedrals setting"); "but no ndihedrals setting");
dihedralflag = tag_require = 1; dihedralflag = tag_require = 1;
dihedrals(flag,line); dihedrals(flag,line);
} else if (strcmp(keyword,"Impropers") == 0) { } else if (strcmp(keyword,"Impropers") == 0) {
if (nimpropers == 0) error->all(FLERR,"Molecule file has impropers " if (nimpropers == 0) error->all(FLERR,"Molecule file has impropers "
"but no nimpropers setting"); "but no nimpropers setting");
improperflag = tag_require = 1; improperflag = tag_require = 1;
impropers(flag,line); impropers(flag,line);
@ -552,26 +566,26 @@ void Molecule::read(int flag)
shakeatomflag = tag_require = 1; shakeatomflag = tag_require = 1;
if (shaketypeflag) shakeflag = 1; if (shaketypeflag) shakeflag = 1;
if (!shakeflagflag) if (!shakeflagflag)
error->all(FLERR,"Molecule file shake flags not before shake atoms"); error->all(FLERR,"Molecule file shake flags not before shake atoms");
if (flag) shakeatom_read(line); if (flag) shakeatom_read(line);
else skip_lines(natoms,line); else skip_lines(natoms,line);
} else if (strcmp(keyword,"Shake Bond Types") == 0) { } else if (strcmp(keyword,"Shake Bond Types") == 0) {
shaketypeflag = 1; shaketypeflag = 1;
if (shakeatomflag) shakeflag = 1; if (shakeatomflag) shakeflag = 1;
if (!shakeflagflag) if (!shakeflagflag)
error->all(FLERR,"Molecule file shake flags not before shake bonds"); error->all(FLERR,"Molecule file shake flags not before shake bonds");
if (flag) shaketype_read(line); if (flag) shaketype_read(line);
else skip_lines(natoms,line); else skip_lines(natoms,line);
} else if (strcmp(keyword,"Body Integers") == 0) { } else if (strcmp(keyword,"Body Integers") == 0) {
if (bodyflag == 0 || nibody == 0) if (bodyflag == 0 || nibody == 0)
error->all(FLERR,"Molecule file has body params " error->all(FLERR,"Molecule file has body params "
"but no setting for them"); "but no setting for them");
ibodyflag = 1; ibodyflag = 1;
body(flag,0,line); body(flag,0,line);
} else if (strcmp(keyword,"Body Doubles") == 0) { } else if (strcmp(keyword,"Body Doubles") == 0) {
if (bodyflag == 0 || ndbody == 0) if (bodyflag == 0 || ndbody == 0)
error->all(FLERR,"Molecule file has body params " error->all(FLERR,"Molecule file has body params "
"but no setting for them"); "but no setting for them");
dbodyflag = 1; dbodyflag = 1;
body(flag,1,line); body(flag,1,line);
@ -618,7 +632,7 @@ void Molecule::read(int flag)
// body particle must have natom = 1 // body particle must have natom = 1
// set radius by having body class compute its own radius // set radius by having body class compute its own radius
if (bodyflag) { if (bodyflag) {
radiusflag = 1; radiusflag = 1;
if (natoms != 1) if (natoms != 1)
@ -641,12 +655,9 @@ void Molecule::coords(char *line)
int tmp; int tmp;
for (int i = 0; i < natoms; i++) { for (int i = 0; i < natoms; i++) {
readline(line); readline(line);
if (i == 0) { if (4 != sscanf(line,"%d %lg %lg %lg",&tmp,&x[i][0],&x[i][1],&x[i][2]))
int nwords = atom->count_words(line); error->all(FLERR,"Invalid Coords section in molecule file");
if (nwords != 4)
error->all(FLERR,"Invalid Coords section in molecule file");
}
sscanf(line,"%d %lg %lg %lg",&tmp,&x[i][0],&x[i][1],&x[i][2]);
x[i][0] *= sizescale; x[i][0] *= sizescale;
x[i][1] *= sizescale; x[i][1] *= sizescale;
x[i][2] *= sizescale; x[i][2] *= sizescale;
@ -669,12 +680,8 @@ void Molecule::types(char *line)
int tmp; int tmp;
for (int i = 0; i < natoms; i++) { for (int i = 0; i < natoms; i++) {
readline(line); readline(line);
if (i == 0) { if (2 != sscanf(line,"%d %d",&tmp,&type[i]))
int nwords = atom->count_words(line); error->all(FLERR,"Invalid Types section in molecule file");
if (nwords != 2)
error->all(FLERR,"Invalid Types section in molecule file");
}
sscanf(line,"%d %d",&tmp,&type[i]);
type[i] += toffset; type[i] += toffset;
} }
@ -695,12 +702,8 @@ void Molecule::charges(char *line)
int tmp; int tmp;
for (int i = 0; i < natoms; i++) { for (int i = 0; i < natoms; i++) {
readline(line); readline(line);
if (i == 0) { if (2 != sscanf(line,"%d %lg",&tmp,&q[i]))
int nwords = atom->count_words(line); error->all(FLERR,"Invalid Charges section in molecule file");
if (nwords != 2)
error->all(FLERR,"Invalid Charges section in molecule file");
}
sscanf(line,"%d %lg",&tmp,&q[i]);
} }
} }
@ -714,12 +717,8 @@ void Molecule::diameters(char *line)
maxradius = 0.0; maxradius = 0.0;
for (int i = 0; i < natoms; i++) { for (int i = 0; i < natoms; i++) {
readline(line); readline(line);
if (i == 0) { if (2 != sscanf(line,"%d %lg",&tmp,&radius[i]))
int nwords = atom->count_words(line); error->all(FLERR,"Invalid Diameters section in molecule file");
if (nwords != 2)
error->all(FLERR,"Invalid Diameters section in molecule file");
}
sscanf(line,"%d %lg",&tmp,&radius[i]);
radius[i] *= sizescale; radius[i] *= sizescale;
radius[i] *= 0.5; radius[i] *= 0.5;
maxradius = MAX(maxradius,radius[i]); maxradius = MAX(maxradius,radius[i]);
@ -739,12 +738,8 @@ void Molecule::masses(char *line)
int tmp; int tmp;
for (int i = 0; i < natoms; i++) { for (int i = 0; i < natoms; i++) {
readline(line); readline(line);
if (i == 0) { if (2 != sscanf(line,"%d %lg",&tmp,&rmass[i]))
int nwords = atom->count_words(line); error->all(FLERR,"Invalid Masses section in molecule file");
if (nwords != 2)
error->all(FLERR,"Invalid Masses section in molecule file");
}
sscanf(line,"%d %lg",&tmp,&rmass[i]);
rmass[i] *= sizescale*sizescale*sizescale; rmass[i] *= sizescale*sizescale*sizescale;
} }
@ -773,17 +768,13 @@ void Molecule::bonds(int flag, char *line)
for (int i = 0; i < nbonds; i++) { for (int i = 0; i < nbonds; i++) {
readline(line); readline(line);
if (i == 0) { if (4 != sscanf(line,"%d %d " TAGINT_FORMAT " " TAGINT_FORMAT,
int nwords = atom->count_words(line); &tmp,&itype,&atom1,&atom2))
if (nwords != 4) error->all(FLERR,"Invalid Bonds section in molecule file");
error->all(FLERR,"Invalid Bonds section in molecule file");
}
sscanf(line,"%d %d " TAGINT_FORMAT " " TAGINT_FORMAT,
&tmp,&itype,&atom1,&atom2);
itype += boffset; itype += boffset;
if (atom1 <= 0 || atom1 > natoms || if (atom1 <= 0 || atom1 > natoms ||
atom2 <= 0 || atom2 > natoms) atom2 <= 0 || atom2 > natoms)
error->one(FLERR,"Invalid atom ID in Bonds section of molecule file"); error->one(FLERR,"Invalid atom ID in Bonds section of molecule file");
if (itype <= 0) if (itype <= 0)
error->one(FLERR,"Invalid bond type in Bonds section of molecule file"); error->one(FLERR,"Invalid bond type in Bonds section of molecule file");
@ -795,10 +786,10 @@ void Molecule::bonds(int flag, char *line)
bond_atom[m][num_bond[m]] = atom2; bond_atom[m][num_bond[m]] = atom2;
num_bond[m]++; num_bond[m]++;
if (newton_bond == 0) { if (newton_bond == 0) {
m = atom2-1; m = atom2-1;
bond_type[m][num_bond[m]] = itype; bond_type[m][num_bond[m]] = itype;
bond_atom[m][num_bond[m]] = atom1; bond_atom[m][num_bond[m]] = atom1;
num_bond[m]++; num_bond[m]++;
} }
} else { } else {
count[atom1-1]++; count[atom1-1]++;
@ -835,13 +826,9 @@ void Molecule::angles(int flag, char *line)
for (int i = 0; i < nangles; i++) { for (int i = 0; i < nangles; i++) {
readline(line); readline(line);
if (i == 0) { if (5 != sscanf(line,"%d %d " TAGINT_FORMAT " " TAGINT_FORMAT " " TAGINT_FORMAT,
int nwords = atom->count_words(line); &tmp,&itype,&atom1,&atom2,&atom3))
if (nwords != 5) error->all(FLERR,"Invalid Angles section in molecule file");
error->all(FLERR,"Invalid Angles section in molecule file");
}
sscanf(line,"%d %d " TAGINT_FORMAT " " TAGINT_FORMAT " " TAGINT_FORMAT,
&tmp,&itype,&atom1,&atom2,&atom3);
itype += aoffset; itype += aoffset;
if (atom1 <= 0 || atom1 > natoms || if (atom1 <= 0 || atom1 > natoms ||
@ -860,24 +847,24 @@ void Molecule::angles(int flag, char *line)
angle_atom3[m][num_angle[m]] = atom3; angle_atom3[m][num_angle[m]] = atom3;
num_angle[m]++; num_angle[m]++;
if (newton_bond == 0) { if (newton_bond == 0) {
m = atom1-1; m = atom1-1;
angle_type[m][num_angle[m]] = itype; angle_type[m][num_angle[m]] = itype;
angle_atom1[m][num_angle[m]] = atom1; angle_atom1[m][num_angle[m]] = atom1;
angle_atom2[m][num_angle[m]] = atom2; angle_atom2[m][num_angle[m]] = atom2;
angle_atom3[m][num_angle[m]] = atom3; angle_atom3[m][num_angle[m]] = atom3;
num_angle[m]++; num_angle[m]++;
m = atom3-1; m = atom3-1;
angle_type[m][num_angle[m]] = itype; angle_type[m][num_angle[m]] = itype;
angle_atom1[m][num_angle[m]] = atom1; angle_atom1[m][num_angle[m]] = atom1;
angle_atom2[m][num_angle[m]] = atom2; angle_atom2[m][num_angle[m]] = atom2;
angle_atom3[m][num_angle[m]] = atom3; angle_atom3[m][num_angle[m]] = atom3;
num_angle[m]++; num_angle[m]++;
} }
} else { } else {
count[atom2-1]++; count[atom2-1]++;
if (newton_bond == 0) { if (newton_bond == 0) {
count[atom1-1]++; count[atom1-1]++;
count[atom3-1]++; count[atom3-1]++;
} }
} }
} }
@ -911,14 +898,10 @@ void Molecule::dihedrals(int flag, char *line)
for (int i = 0; i < ndihedrals; i++) { for (int i = 0; i < ndihedrals; i++) {
readline(line); readline(line);
if (i == 0) { if (6 != sscanf(line,"%d %d " TAGINT_FORMAT " " TAGINT_FORMAT " "
int nwords = atom->count_words(line);
if (nwords != 6)
error->all(FLERR,"Invalid Dihedrals section in molecule file");
}
sscanf(line,"%d %d " TAGINT_FORMAT " " TAGINT_FORMAT " "
TAGINT_FORMAT " " TAGINT_FORMAT " ", TAGINT_FORMAT " " TAGINT_FORMAT " ",
&tmp,&itype,&atom1,&atom2,&atom3,&atom4); &tmp,&itype,&atom1,&atom2,&atom3,&atom4))
error->all(FLERR,"Invalid Dihedrals section in molecule file");
itype += doffset; itype += doffset;
if (atom1 <= 0 || atom1 > natoms || if (atom1 <= 0 || atom1 > natoms ||
@ -926,10 +909,10 @@ void Molecule::dihedrals(int flag, char *line)
atom3 <= 0 || atom3 > natoms || atom3 <= 0 || atom3 > natoms ||
atom4 <= 0 || atom4 > natoms) atom4 <= 0 || atom4 > natoms)
error->one(FLERR, error->one(FLERR,
"Invalid atom ID in dihedrals section of molecule file"); "Invalid atom ID in dihedrals section of molecule file");
if (itype <= 0) if (itype <= 0)
error->one(FLERR, error->one(FLERR,
"Invalid dihedral type in dihedrals section of molecule file"); "Invalid dihedral type in dihedrals section of molecule file");
if (flag) { if (flag) {
m = atom2-1; m = atom2-1;
@ -941,34 +924,34 @@ void Molecule::dihedrals(int flag, char *line)
dihedral_atom4[m][num_dihedral[m]] = atom4; dihedral_atom4[m][num_dihedral[m]] = atom4;
num_dihedral[m]++; num_dihedral[m]++;
if (newton_bond == 0) { if (newton_bond == 0) {
m = atom1-1; m = atom1-1;
dihedral_type[m][num_dihedral[m]] = itype; dihedral_type[m][num_dihedral[m]] = itype;
dihedral_atom1[m][num_dihedral[m]] = atom1; dihedral_atom1[m][num_dihedral[m]] = atom1;
dihedral_atom2[m][num_dihedral[m]] = atom2; dihedral_atom2[m][num_dihedral[m]] = atom2;
dihedral_atom3[m][num_dihedral[m]] = atom3; dihedral_atom3[m][num_dihedral[m]] = atom3;
dihedral_atom4[m][num_dihedral[m]] = atom4; dihedral_atom4[m][num_dihedral[m]] = atom4;
num_dihedral[m]++; num_dihedral[m]++;
m = atom3-1; m = atom3-1;
dihedral_type[m][num_dihedral[m]] = itype; dihedral_type[m][num_dihedral[m]] = itype;
dihedral_atom1[m][num_dihedral[m]] = atom1; dihedral_atom1[m][num_dihedral[m]] = atom1;
dihedral_atom2[m][num_dihedral[m]] = atom2; dihedral_atom2[m][num_dihedral[m]] = atom2;
dihedral_atom3[m][num_dihedral[m]] = atom3; dihedral_atom3[m][num_dihedral[m]] = atom3;
dihedral_atom4[m][num_dihedral[m]] = atom4; dihedral_atom4[m][num_dihedral[m]] = atom4;
num_dihedral[m]++; num_dihedral[m]++;
m = atom4-1; m = atom4-1;
dihedral_type[m][num_dihedral[m]] = itype; dihedral_type[m][num_dihedral[m]] = itype;
dihedral_atom1[m][num_dihedral[m]] = atom1; dihedral_atom1[m][num_dihedral[m]] = atom1;
dihedral_atom2[m][num_dihedral[m]] = atom2; dihedral_atom2[m][num_dihedral[m]] = atom2;
dihedral_atom3[m][num_dihedral[m]] = atom3; dihedral_atom3[m][num_dihedral[m]] = atom3;
dihedral_atom4[m][num_dihedral[m]] = atom4; dihedral_atom4[m][num_dihedral[m]] = atom4;
num_dihedral[m]++; num_dihedral[m]++;
} }
} else { } else {
count[atom2-1]++; count[atom2-1]++;
if (newton_bond == 0) { if (newton_bond == 0) {
count[atom1-1]++; count[atom1-1]++;
count[atom3-1]++; count[atom3-1]++;
count[atom4-1]++; count[atom4-1]++;
} }
} }
} }
@ -1002,14 +985,10 @@ void Molecule::impropers(int flag, char *line)
for (int i = 0; i < nimpropers; i++) { for (int i = 0; i < nimpropers; i++) {
readline(line); readline(line);
if (i == 0) { if (6 != sscanf(line,"%d %d " TAGINT_FORMAT " " TAGINT_FORMAT " "
int nwords = atom->count_words(line);
if (nwords != 6)
error->all(FLERR,"Invalid Impropers section in molecule file");
}
sscanf(line,"%d %d " TAGINT_FORMAT " " TAGINT_FORMAT " "
TAGINT_FORMAT " " TAGINT_FORMAT " ", TAGINT_FORMAT " " TAGINT_FORMAT " ",
&tmp,&itype,&atom1,&atom2,&atom3,&atom4); &tmp,&itype,&atom1,&atom2,&atom3,&atom4))
error->all(FLERR,"Invalid Impropers section in molecule file");
itype += ioffset; itype += ioffset;
if (atom1 <= 0 || atom1 > natoms || if (atom1 <= 0 || atom1 > natoms ||
@ -1017,10 +996,10 @@ void Molecule::impropers(int flag, char *line)
atom3 <= 0 || atom3 > natoms || atom3 <= 0 || atom3 > natoms ||
atom4 <= 0 || atom4 > natoms) atom4 <= 0 || atom4 > natoms)
error->one(FLERR, error->one(FLERR,
"Invalid atom ID in impropers section of molecule file"); "Invalid atom ID in impropers section of molecule file");
if (itype <= 0) if (itype <= 0)
error->one(FLERR, error->one(FLERR,
"Invalid improper type in impropers section of molecule file"); "Invalid improper type in impropers section of molecule file");
if (flag) { if (flag) {
m = atom2-1; m = atom2-1;
@ -1032,34 +1011,34 @@ void Molecule::impropers(int flag, char *line)
improper_atom4[m][num_improper[m]] = atom4; improper_atom4[m][num_improper[m]] = atom4;
num_improper[m]++; num_improper[m]++;
if (newton_bond == 0) { if (newton_bond == 0) {
m = atom1-1; m = atom1-1;
improper_type[m][num_improper[m]] = itype; improper_type[m][num_improper[m]] = itype;
improper_atom1[m][num_improper[m]] = atom1; improper_atom1[m][num_improper[m]] = atom1;
improper_atom2[m][num_improper[m]] = atom2; improper_atom2[m][num_improper[m]] = atom2;
improper_atom3[m][num_improper[m]] = atom3; improper_atom3[m][num_improper[m]] = atom3;
improper_atom4[m][num_improper[m]] = atom4; improper_atom4[m][num_improper[m]] = atom4;
num_improper[m]++; num_improper[m]++;
m = atom3-1; m = atom3-1;
improper_type[m][num_improper[m]] = itype; improper_type[m][num_improper[m]] = itype;
improper_atom1[m][num_improper[m]] = atom1; improper_atom1[m][num_improper[m]] = atom1;
improper_atom2[m][num_improper[m]] = atom2; improper_atom2[m][num_improper[m]] = atom2;
improper_atom3[m][num_improper[m]] = atom3; improper_atom3[m][num_improper[m]] = atom3;
improper_atom4[m][num_improper[m]] = atom4; improper_atom4[m][num_improper[m]] = atom4;
num_improper[m]++; num_improper[m]++;
m = atom4-1; m = atom4-1;
improper_type[m][num_improper[m]] = itype; improper_type[m][num_improper[m]] = itype;
improper_atom1[m][num_improper[m]] = atom1; improper_atom1[m][num_improper[m]] = atom1;
improper_atom2[m][num_improper[m]] = atom2; improper_atom2[m][num_improper[m]] = atom2;
improper_atom3[m][num_improper[m]] = atom3; improper_atom3[m][num_improper[m]] = atom3;
improper_atom4[m][num_improper[m]] = atom4; improper_atom4[m][num_improper[m]] = atom4;
num_improper[m]++; num_improper[m]++;
} }
} else { } else {
count[atom2-1]++; count[atom2-1]++;
if (newton_bond == 0) { if (newton_bond == 0) {
count[atom1-1]++; count[atom1-1]++;
count[atom3-1]++; count[atom3-1]++;
count[atom4-1]++; count[atom4-1]++;
} }
} }
} }
@ -1087,13 +1066,9 @@ void Molecule::nspecial_read(int flag, char *line)
for (int i = 0; i < natoms; i++) { for (int i = 0; i < natoms; i++) {
readline(line); readline(line);
if (i == 0) { if (4 != sscanf(line,"%d %d %d %d",&tmp,&c1,&c2,&c3))
int nwords = atom->count_words(line); error->all(FLERR,"Invalid Special Bond Counts section in "
if (nwords != 4) "molecule file");
error->all(FLERR,"Invalid Special Bond Counts section in "
"molecule file");
}
sscanf(line,"%d %d %d %d",&tmp,&c1,&c2,&c3);
if (flag) { if (flag) {
nspecial[i][0] = c1; nspecial[i][0] = c1;
@ -1117,13 +1092,13 @@ void Molecule::special_read(char *line)
nwords = parse(line,words,maxspecial+1); nwords = parse(line,words,maxspecial+1);
if (nwords != nspecial[i][2]+1) if (nwords != nspecial[i][2]+1)
error->all(FLERR,"Molecule file special list " error->all(FLERR,"Molecule file special list "
"does not match special count"); "does not match special count");
for (m = 1; m < nwords; m++) { for (m = 1; m < nwords; m++) {
special[i][m-1] = ATOTAGINT(words[m]); special[i][m-1] = ATOTAGINT(words[m]);
if (special[i][m-1] <= 0 || special[i][m-1] > natoms || if (special[i][m-1] <= 0 || special[i][m-1] > natoms ||
special[i][m-1] == i+1) special[i][m-1] == i+1)
error->all(FLERR,"Invalid special atom index in molecule file"); error->all(FLERR,"Invalid special atom index in molecule file");
} }
} }
@ -1229,7 +1204,8 @@ void Molecule::shakeflag_read(char *line)
int tmp; int tmp;
for (int i = 0; i < natoms; i++) { for (int i = 0; i < natoms; i++) {
readline(line); readline(line);
sscanf(line,"%d %d",&tmp,&shake_flag[i]); if (2 != sscanf(line,"%d %d",&tmp,&shake_flag[i]))
error->all(FLERR,"Invalid Shake Flags section in molecule file");
} }
for (int i = 0; i < natoms; i++) for (int i = 0; i < natoms; i++)
@ -1243,23 +1219,32 @@ void Molecule::shakeflag_read(char *line)
void Molecule::shakeatom_read(char *line) void Molecule::shakeatom_read(char *line)
{ {
int tmp; int tmp, nmatch, nwant;
for (int i = 0; i < natoms; i++) { for (int i = 0; i < natoms; i++) {
readline(line); readline(line);
if (shake_flag[i] == 1) if (shake_flag[i] == 1) {
sscanf(line,"%d " TAGINT_FORMAT " " TAGINT_FORMAT " " TAGINT_FORMAT, nmatch = sscanf(line,"%d " TAGINT_FORMAT " " TAGINT_FORMAT
&tmp,&shake_atom[i][0],&shake_atom[i][1],&shake_atom[i][2]); " " TAGINT_FORMAT,&tmp,&shake_atom[i][0],
else if (shake_flag[i] == 2) &shake_atom[i][1],&shake_atom[i][2]);
sscanf(line,"%d " TAGINT_FORMAT " " TAGINT_FORMAT, nwant = 4;
&tmp,&shake_atom[i][0],&shake_atom[i][1]); } else if (shake_flag[i] == 2) {
else if (shake_flag[i] == 3) nmatch = sscanf(line,"%d " TAGINT_FORMAT " " TAGINT_FORMAT,
sscanf(line,"%d " TAGINT_FORMAT " " TAGINT_FORMAT " " TAGINT_FORMAT, &tmp,&shake_atom[i][0],&shake_atom[i][1]);
&tmp,&shake_atom[i][0],&shake_atom[i][1],&shake_atom[i][2]); nwant = 3;
else if (shake_flag[i] == 4) } else if (shake_flag[i] == 3) {
sscanf(line,"%d " TAGINT_FORMAT " " TAGINT_FORMAT " " nmatch = sscanf(line,"%d " TAGINT_FORMAT " " TAGINT_FORMAT
TAGINT_FORMAT " " TAGINT_FORMAT, " " TAGINT_FORMAT,&tmp,&shake_atom[i][0],
&tmp,&shake_atom[i][0],&shake_atom[i][1], &shake_atom[i][1],&shake_atom[i][2]);
&shake_atom[i][2],&shake_atom[i][3]); nwant = 4;
} else if (shake_flag[i] == 4) {
nmatch = sscanf(line,"%d " TAGINT_FORMAT " " TAGINT_FORMAT " "
TAGINT_FORMAT " " TAGINT_FORMAT,
&tmp,&shake_atom[i][0],&shake_atom[i][1],
&shake_atom[i][2],&shake_atom[i][3]);
nwant = 5;
}
if (nmatch != nwant)
error->all(FLERR,"Invalid shake atom in molecule file");
} }
for (int i = 0; i < natoms; i++) { for (int i = 0; i < natoms; i++) {
@ -1277,19 +1262,27 @@ void Molecule::shakeatom_read(char *line)
void Molecule::shaketype_read(char *line) void Molecule::shaketype_read(char *line)
{ {
int tmp; int tmp,nmatch,nwant;
for (int i = 0; i < natoms; i++) { for (int i = 0; i < natoms; i++) {
readline(line); readline(line);
if (shake_flag[i] == 1) if (shake_flag[i] == 1) {
sscanf(line,"%d %d %d %d",&tmp, nmatch = sscanf(line,"%d %d %d %d",&tmp,&shake_type[i][0],
&shake_type[i][0],&shake_type[i][1],&shake_type[i][2]); &shake_type[i][1],&shake_type[i][2]);
else if (shake_flag[i] == 2) nwant = 4;
sscanf(line,"%d %d",&tmp,&shake_type[i][0]); } else if (shake_flag[i] == 2) {
else if (shake_flag[i] == 3) nmatch = sscanf(line,"%d %d",&tmp,&shake_type[i][0]);
sscanf(line,"%d %d %d",&tmp,&shake_type[i][0],&shake_type[i][1]); nwant = 2;
else if (shake_flag[i] == 4) } else if (shake_flag[i] == 3) {
sscanf(line,"%d %d %d %d",&tmp, nmatch = sscanf(line,"%d %d %d",&tmp,&shake_type[i][0],
&shake_type[i][0],&shake_type[i][1],&shake_type[i][2]); &shake_type[i][1]);
nwant = 3;
} else if (shake_flag[i] == 4) {
nmatch = sscanf(line,"%d %d %d %d",&tmp,&shake_type[i][0],
&shake_type[i][1],&shake_type[i][2]);
nwant = 4;
}
if (nmatch != nwant)
error->all(FLERR,"Invalid shake type data in molecule file");
} }
for (int i = 0; i < natoms; i++) { for (int i = 0; i < natoms; i++) {
@ -1501,46 +1494,46 @@ void Molecule::allocate()
if (bondflag) { if (bondflag) {
memory->create(bond_type,natoms,bond_per_atom, memory->create(bond_type,natoms,bond_per_atom,
"molecule:bond_type"); "molecule:bond_type");
memory->create(bond_atom,natoms,bond_per_atom, memory->create(bond_atom,natoms,bond_per_atom,
"molecule:bond_atom"); "molecule:bond_atom");
} }
if (angleflag) { if (angleflag) {
memory->create(angle_type,natoms,angle_per_atom, memory->create(angle_type,natoms,angle_per_atom,
"molecule:angle_type"); "molecule:angle_type");
memory->create(angle_atom1,natoms,angle_per_atom, memory->create(angle_atom1,natoms,angle_per_atom,
"molecule:angle_atom1"); "molecule:angle_atom1");
memory->create(angle_atom2,natoms,angle_per_atom, memory->create(angle_atom2,natoms,angle_per_atom,
"molecule:angle_atom2"); "molecule:angle_atom2");
memory->create(angle_atom3,natoms,angle_per_atom, memory->create(angle_atom3,natoms,angle_per_atom,
"molecule:angle_atom3"); "molecule:angle_atom3");
} }
if (dihedralflag) { if (dihedralflag) {
memory->create(dihedral_type,natoms,dihedral_per_atom, memory->create(dihedral_type,natoms,dihedral_per_atom,
"molecule:dihedral_type"); "molecule:dihedral_type");
memory->create(dihedral_atom1,natoms,dihedral_per_atom, memory->create(dihedral_atom1,natoms,dihedral_per_atom,
"molecule:dihedral_atom1"); "molecule:dihedral_atom1");
memory->create(dihedral_atom2,natoms,dihedral_per_atom, memory->create(dihedral_atom2,natoms,dihedral_per_atom,
"molecule:dihedral_atom2"); "molecule:dihedral_atom2");
memory->create(dihedral_atom3,natoms,dihedral_per_atom, memory->create(dihedral_atom3,natoms,dihedral_per_atom,
"molecule:dihedral_atom3"); "molecule:dihedral_atom3");
memory->create(dihedral_atom4,natoms,dihedral_per_atom, memory->create(dihedral_atom4,natoms,dihedral_per_atom,
"molecule:dihedral_atom4"); "molecule:dihedral_atom4");
} }
if (improperflag) { if (improperflag) {
memory->create(improper_type,natoms,improper_per_atom, memory->create(improper_type,natoms,improper_per_atom,
"molecule:improper_type"); "molecule:improper_type");
memory->create(improper_atom1,natoms,improper_per_atom, memory->create(improper_atom1,natoms,improper_per_atom,
"molecule:improper_atom1"); "molecule:improper_atom1");
memory->create(improper_atom2,natoms,improper_per_atom, memory->create(improper_atom2,natoms,improper_per_atom,
"molecule:improper_atom2"); "molecule:improper_atom2");
memory->create(improper_atom3,natoms,improper_per_atom, memory->create(improper_atom3,natoms,improper_per_atom,
"molecule:improper_atom3"); "molecule:improper_atom3");
memory->create(improper_atom4,natoms,improper_per_atom, memory->create(improper_atom4,natoms,improper_per_atom,
"molecule:improper_atom4"); "molecule:improper_atom4");
} }
if (shakeflag) { if (shakeflag) {
@ -1653,7 +1646,7 @@ void Molecule::parse_keyword(int flag, char *line, char *keyword)
if (me == 0) { if (me == 0) {
if (fgets(line,MAXLINE,fp) == NULL) eof = 1; if (fgets(line,MAXLINE,fp) == NULL) eof = 1;
while (eof == 0 && strspn(line," \t\n\r") == strlen(line)) { while (eof == 0 && strspn(line," \t\n\r") == strlen(line)) {
if (fgets(line,MAXLINE,fp) == NULL) eof = 1; if (fgets(line,MAXLINE,fp) == NULL) eof = 1;
} }
if (fgets(keyword,MAXLINE,fp) == NULL) eof = 1; if (fgets(keyword,MAXLINE,fp) == NULL) eof = 1;
} }