enable and apply clang-format
This commit is contained in:
221
src/group.cpp
221
src/group.cpp
@ -1,4 +1,3 @@
|
||||
// clang-format off
|
||||
/* ----------------------------------------------------------------------
|
||||
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
|
||||
https://www.lammps.org/, Sandia National Laboratories
|
||||
@ -20,18 +19,18 @@
|
||||
#include "domain.h"
|
||||
#include "dump.h"
|
||||
#include "error.h"
|
||||
#include "exceptions.h"
|
||||
#include "fix.h"
|
||||
#include "force.h"
|
||||
#include "input.h"
|
||||
#include "math_extra.h"
|
||||
#include "math_eigen.h"
|
||||
#include "math_extra.h"
|
||||
#include "memory.h"
|
||||
#include "modify.h"
|
||||
#include "output.h"
|
||||
#include "region.h"
|
||||
#include "tokenizer.h"
|
||||
#include "variable.h"
|
||||
#include "exceptions.h"
|
||||
|
||||
#include <cmath>
|
||||
#include <cstring>
|
||||
@ -93,8 +92,7 @@ void Group::assign(int narg, char **arg)
|
||||
{
|
||||
int i;
|
||||
|
||||
if (domain->box_exist == 0)
|
||||
error->all(FLERR,"Group command before simulation box is defined");
|
||||
if (domain->box_exist == 0) error->all(FLERR, "Group command before simulation box is defined");
|
||||
if (narg < 2) utils::missing_cmd_args(FLERR, "group", error);
|
||||
|
||||
// delete the group if not being used elsewhere
|
||||
@ -122,8 +120,7 @@ void Group::assign(int narg, char **arg)
|
||||
int bits = inversemask[igroup];
|
||||
for (i = 0; i < nlocal; i++) mask[i] &= bits;
|
||||
|
||||
if (dynamic[igroup])
|
||||
modify->delete_fix(std::string("GROUP_") + names[igroup]);
|
||||
if (dynamic[igroup]) modify->delete_fix(std::string("GROUP_") + names[igroup]);
|
||||
|
||||
delete[] names[igroup];
|
||||
names[igroup] = nullptr;
|
||||
@ -181,8 +178,7 @@ void Group::assign(int narg, char **arg)
|
||||
region->prematch();
|
||||
|
||||
for (i = 0; i < nlocal; i++)
|
||||
if (region->match(x[i][0],x[i][1],x[i][2]))
|
||||
mask[i] |= bit;
|
||||
if (region->match(x[i][0], x[i][1], x[i][2])) mask[i] |= bit;
|
||||
|
||||
// create an empty group
|
||||
|
||||
@ -200,9 +196,12 @@ void Group::assign(int narg, char **arg)
|
||||
if (narg < 3) utils::missing_cmd_args(FLERR, std::string("group ") + arg[1], error);
|
||||
|
||||
int category = NONE;
|
||||
if (strcmp(arg[1],"type") == 0) category = TYPE;
|
||||
else if (strcmp(arg[1],"molecule") == 0) category = MOLECULE;
|
||||
else if (strcmp(arg[1],"id") == 0) category = ID;
|
||||
if (strcmp(arg[1], "type") == 0)
|
||||
category = TYPE;
|
||||
else if (strcmp(arg[1], "molecule") == 0)
|
||||
category = MOLECULE;
|
||||
else if (strcmp(arg[1], "id") == 0)
|
||||
category = ID;
|
||||
|
||||
if ((category == MOLECULE) && (!atom->molecule_flag))
|
||||
error->all(FLERR, "Group molecule command requires atom attribute molecule");
|
||||
@ -213,39 +212,52 @@ void Group::assign(int narg, char **arg)
|
||||
// args = logical condition
|
||||
|
||||
if (narg > 3 &&
|
||||
(strcmp(arg[2],"<") == 0 || strcmp(arg[2],">") == 0 ||
|
||||
strcmp(arg[2],"<=") == 0 || strcmp(arg[2],">=") == 0 ||
|
||||
strcmp(arg[2],"==") == 0 || strcmp(arg[2],"!=") == 0 ||
|
||||
(strcmp(arg[2], "<") == 0 || strcmp(arg[2], ">") == 0 || strcmp(arg[2], "<=") == 0 ||
|
||||
strcmp(arg[2], ">=") == 0 || strcmp(arg[2], "==") == 0 || strcmp(arg[2], "!=") == 0 ||
|
||||
strcmp(arg[2], "<>") == 0)) {
|
||||
|
||||
int condition = -1;
|
||||
if (strcmp(arg[2],"<") == 0) condition = LT;
|
||||
else if (strcmp(arg[2],"<=") == 0) condition = LE;
|
||||
else if (strcmp(arg[2],">") == 0) condition = GT;
|
||||
else if (strcmp(arg[2],">=") == 0) condition = GE;
|
||||
else if (strcmp(arg[2],"==") == 0) condition = EQ;
|
||||
else if (strcmp(arg[2],"!=") == 0) condition = NEQ;
|
||||
else if (strcmp(arg[2],"<>") == 0) condition = BETWEEN;
|
||||
else error->all(FLERR,"Illegal group command");
|
||||
if (strcmp(arg[2], "<") == 0)
|
||||
condition = LT;
|
||||
else if (strcmp(arg[2], "<=") == 0)
|
||||
condition = LE;
|
||||
else if (strcmp(arg[2], ">") == 0)
|
||||
condition = GT;
|
||||
else if (strcmp(arg[2], ">=") == 0)
|
||||
condition = GE;
|
||||
else if (strcmp(arg[2], "==") == 0)
|
||||
condition = EQ;
|
||||
else if (strcmp(arg[2], "!=") == 0)
|
||||
condition = NEQ;
|
||||
else if (strcmp(arg[2], "<>") == 0)
|
||||
condition = BETWEEN;
|
||||
else
|
||||
error->all(FLERR, "Illegal group command");
|
||||
|
||||
tagint bound1, bound2;
|
||||
if (category == TYPE)
|
||||
bound1 = (tagint) utils::expand_type_int(FLERR, arg[3], Atom::ATOM, lmp);
|
||||
else bound1 = utils::tnumeric(FLERR, arg[3], false, lmp);
|
||||
else
|
||||
bound1 = utils::tnumeric(FLERR, arg[3], false, lmp);
|
||||
bound2 = -1;
|
||||
|
||||
if (condition == BETWEEN) {
|
||||
if (narg != 5) error->all(FLERR, "Illegal group command");
|
||||
if (category == TYPE)
|
||||
bound2 = (tagint) utils::expand_type_int(FLERR, arg[4], Atom::ATOM, lmp);
|
||||
else bound2 = utils::tnumeric(FLERR, arg[4], false, lmp);
|
||||
} else if (narg != 4) error->all(FLERR,"Illegal group command");
|
||||
else
|
||||
bound2 = utils::tnumeric(FLERR, arg[4], false, lmp);
|
||||
} else if (narg != 4)
|
||||
error->all(FLERR, "Illegal group command");
|
||||
|
||||
int *attribute = nullptr;
|
||||
tagint *tattribute = nullptr;
|
||||
if (category == TYPE) attribute = atom->type;
|
||||
else if (category == MOLECULE) tattribute = atom->molecule;
|
||||
else if (category == ID) tattribute = atom->tag;
|
||||
if (category == TYPE)
|
||||
attribute = atom->type;
|
||||
else if (category == MOLECULE)
|
||||
tattribute = atom->molecule;
|
||||
else if (category == ID)
|
||||
tattribute = atom->tag;
|
||||
|
||||
// add to group if meets condition
|
||||
|
||||
@ -270,8 +282,7 @@ void Group::assign(int narg, char **arg)
|
||||
if (attribute[i] != bound1) mask[i] |= bit;
|
||||
} else if (condition == BETWEEN) {
|
||||
for (i = 0; i < nlocal; i++)
|
||||
if (attribute[i] >= bound1 && attribute[i] <= bound2)
|
||||
mask[i] |= bit;
|
||||
if (attribute[i] >= bound1 && attribute[i] <= bound2) mask[i] |= bit;
|
||||
}
|
||||
} else {
|
||||
if (condition == LT) {
|
||||
@ -294,8 +305,7 @@ void Group::assign(int narg, char **arg)
|
||||
if (tattribute[i] != bound1) mask[i] |= bit;
|
||||
} else if (condition == BETWEEN) {
|
||||
for (i = 0; i < nlocal; i++)
|
||||
if (tattribute[i] >= bound1 && tattribute[i] <= bound2)
|
||||
mask[i] |= bit;
|
||||
if (tattribute[i] >= bound1 && tattribute[i] <= bound2) mask[i] |= bit;
|
||||
}
|
||||
}
|
||||
|
||||
@ -304,9 +314,12 @@ void Group::assign(int narg, char **arg)
|
||||
} else {
|
||||
int *attribute = nullptr;
|
||||
tagint *tattribute = nullptr;
|
||||
if (category == TYPE) attribute = atom->type;
|
||||
else if (category == MOLECULE) tattribute = atom->molecule;
|
||||
else if (category == ID) tattribute = atom->tag;
|
||||
if (category == TYPE)
|
||||
attribute = atom->type;
|
||||
else if (category == MOLECULE)
|
||||
tattribute = atom->molecule;
|
||||
else if (category == ID)
|
||||
tattribute = atom->tag;
|
||||
|
||||
char *typestr = nullptr;
|
||||
tagint start, stop, delta;
|
||||
@ -329,12 +342,12 @@ void Group::assign(int narg, char **arg)
|
||||
} else if (utils::strmatch(arg[iarg], "^-?\\d+:-?\\d+:\\d+$")) {
|
||||
stop = values.next_tagint();
|
||||
delta = values.next_tagint();
|
||||
} else throw TokenizerException("Syntax error","");
|
||||
} else
|
||||
throw TokenizerException("Syntax error", "");
|
||||
} catch (TokenizerException &e) {
|
||||
error->all(FLERR, "Incorrect range string '{}': {}", arg[iarg], e.what());
|
||||
}
|
||||
if (delta < 1)
|
||||
error->all(FLERR,"Illegal range increment value");
|
||||
if (delta < 1) error->all(FLERR, "Illegal range increment value");
|
||||
}
|
||||
|
||||
// add to group if attribute matches value or sequence
|
||||
@ -342,11 +355,13 @@ void Group::assign(int narg, char **arg)
|
||||
if (attribute) {
|
||||
for (i = 0; i < nlocal; i++)
|
||||
if (attribute[i] >= start && attribute[i] <= stop &&
|
||||
(attribute[i]-start) % delta == 0) mask[i] |= bit;
|
||||
(attribute[i] - start) % delta == 0)
|
||||
mask[i] |= bit;
|
||||
} else {
|
||||
for (i = 0; i < nlocal; i++)
|
||||
if (tattribute[i] >= start && tattribute[i] <= stop &&
|
||||
(tattribute[i]-start) % delta == 0) mask[i] |= bit;
|
||||
(tattribute[i] - start) % delta == 0)
|
||||
mask[i] |= bit;
|
||||
}
|
||||
}
|
||||
|
||||
@ -388,7 +403,8 @@ void Group::assign(int narg, char **arg)
|
||||
|
||||
add_molecules(igroup, bit);
|
||||
|
||||
} else error->all(FLERR,"Unknown group include keyword {}",arg[2]);
|
||||
} else
|
||||
error->all(FLERR, "Unknown group include keyword {}", arg[2]);
|
||||
|
||||
// style = subtract
|
||||
|
||||
@ -465,8 +481,7 @@ void Group::assign(int narg, char **arg)
|
||||
for (int iarg = 2; iarg < narg; iarg++) {
|
||||
jgroup = find(arg[iarg]);
|
||||
if (jgroup == -1) error->all(FLERR, "Group ID {} does not exist", arg[iarg]);
|
||||
if (dynamic[jgroup])
|
||||
error->all(FLERR,"Cannot intersect groups using a dynamic group");
|
||||
if (dynamic[jgroup]) error->all(FLERR, "Cannot intersect groups using a dynamic group");
|
||||
list[iarg - 2] = jgroup;
|
||||
}
|
||||
|
||||
@ -489,16 +504,14 @@ void Group::assign(int narg, char **arg)
|
||||
} else if (strcmp(arg[1], "dynamic") == 0) {
|
||||
|
||||
if (narg < 4) error->all(FLERR, "Illegal group command");
|
||||
if (strcmp(arg[0],arg[2]) == 0)
|
||||
error->all(FLERR,"Group dynamic cannot reference itself");
|
||||
if (strcmp(arg[0], arg[2]) == 0) error->all(FLERR, "Group dynamic cannot reference itself");
|
||||
if (find(arg[2]) < 0)
|
||||
error->all(FLERR, "Group dynamic parent group {} does not exist", arg[2]);
|
||||
if (igroup == 0) error->all(FLERR, "Group all cannot be made dynamic");
|
||||
|
||||
// if group is already dynamic, delete existing FixGroup
|
||||
|
||||
if (dynamic[igroup])
|
||||
modify->delete_fix(std::string("GROUP_") + names[igroup]);
|
||||
if (dynamic[igroup]) modify->delete_fix(std::string("GROUP_") + names[igroup]);
|
||||
|
||||
dynamic[igroup] = 1;
|
||||
|
||||
@ -514,14 +527,14 @@ void Group::assign(int narg, char **arg)
|
||||
|
||||
if (narg != 2) error->all(FLERR, "Illegal group static command");
|
||||
|
||||
if (dynamic[igroup])
|
||||
modify->delete_fix(std::string("GROUP_") + names[igroup]);
|
||||
if (dynamic[igroup]) modify->delete_fix(std::string("GROUP_") + names[igroup]);
|
||||
|
||||
dynamic[igroup] = 0;
|
||||
|
||||
// not a valid group style
|
||||
|
||||
} else error->all(FLERR,"Unknown group command keyword: {}",arg[1]);
|
||||
} else
|
||||
error->all(FLERR, "Unknown group command keyword: {}", arg[1]);
|
||||
|
||||
} catch (LAMMPSException &e) {
|
||||
// undo created group in case of an error
|
||||
@ -537,7 +550,8 @@ void Group::assign(int narg, char **arg)
|
||||
|
||||
int n;
|
||||
n = 0;
|
||||
for (i = 0; i < nlocal; i++) if (mask[i] & bit) n++;
|
||||
for (i = 0; i < nlocal; i++)
|
||||
if (mask[i] & bit) n++;
|
||||
|
||||
double rlocal = n;
|
||||
double all;
|
||||
@ -560,9 +574,7 @@ void Group::assign(const std::string &groupcmd)
|
||||
auto args = utils::split_words(groupcmd);
|
||||
std::vector<char *> newarg(args.size());
|
||||
int i = 0;
|
||||
for (const auto &arg : args) {
|
||||
newarg[i++] = (char *)arg.c_str();
|
||||
}
|
||||
for (const auto &arg : args) { newarg[i++] = (char *) arg.c_str(); }
|
||||
assign(args.size(), newarg.data());
|
||||
}
|
||||
|
||||
@ -715,8 +727,10 @@ void Group::write_restart(FILE *fp)
|
||||
int n;
|
||||
int count = 0;
|
||||
for (int i = 0; i < MAX_GROUP; i++) {
|
||||
if (names[i]) n = strlen(names[i]) + 1;
|
||||
else n = 0;
|
||||
if (names[i])
|
||||
n = strlen(names[i]) + 1;
|
||||
else
|
||||
n = 0;
|
||||
fwrite(&n, sizeof(int), 1, fp);
|
||||
if (n) {
|
||||
fwrite(names[i], sizeof(char), n, fp);
|
||||
@ -759,7 +773,8 @@ void Group::read_restart(FILE *fp)
|
||||
if (me == 0) utils::sfread(FLERR, names[i], sizeof(char), n, fp, nullptr, error);
|
||||
MPI_Bcast(names[i], n, MPI_CHAR, 0, world);
|
||||
count++;
|
||||
} else names[i] = nullptr;
|
||||
} else
|
||||
names[i] = nullptr;
|
||||
}
|
||||
}
|
||||
|
||||
@ -874,12 +889,10 @@ double Group::mass(int igroup, Region *region)
|
||||
|
||||
if (rmass) {
|
||||
for (int i = 0; i < nlocal; i++)
|
||||
if (mask[i] & groupbit && region->match(x[i][0],x[i][1],x[i][2]))
|
||||
one += rmass[i];
|
||||
if (mask[i] & groupbit && region->match(x[i][0], x[i][1], x[i][2])) one += rmass[i];
|
||||
} else {
|
||||
for (int i = 0; i < nlocal; i++)
|
||||
if (mask[i] & groupbit && region->match(x[i][0],x[i][1],x[i][2]))
|
||||
one += mass[type[i]];
|
||||
if (mask[i] & groupbit && region->match(x[i][0], x[i][1], x[i][2])) one += mass[type[i]];
|
||||
}
|
||||
|
||||
double all;
|
||||
@ -924,8 +937,7 @@ double Group::charge(int igroup, Region *region)
|
||||
|
||||
double qone = 0.0;
|
||||
for (int i = 0; i < nlocal; i++)
|
||||
if (mask[i] & groupbit && region->match(x[i][0],x[i][1],x[i][2]))
|
||||
qone += q[i];
|
||||
if (mask[i] & groupbit && region->match(x[i][0], x[i][1], x[i][2])) qone += q[i];
|
||||
|
||||
double qall;
|
||||
MPI_Allreduce(&qone, &qall, 1, MPI_DOUBLE, MPI_SUM, world);
|
||||
@ -1292,13 +1304,11 @@ double Group::ke(int igroup)
|
||||
if (rmass) {
|
||||
for (int i = 0; i < nlocal; i++)
|
||||
if (mask[i] & groupbit)
|
||||
one += (v[i][0]*v[i][0] + v[i][1]*v[i][1] + v[i][2]*v[i][2]) *
|
||||
rmass[i];
|
||||
one += (v[i][0] * v[i][0] + v[i][1] * v[i][1] + v[i][2] * v[i][2]) * rmass[i];
|
||||
} else {
|
||||
for (int i = 0; i < nlocal; i++)
|
||||
if (mask[i] & groupbit)
|
||||
one += (v[i][0]*v[i][0] + v[i][1]*v[i][1] + v[i][2]*v[i][2]) *
|
||||
mass[type[i]];
|
||||
one += (v[i][0] * v[i][0] + v[i][1] * v[i][1] + v[i][2] * v[i][2]) * mass[type[i]];
|
||||
}
|
||||
|
||||
double all;
|
||||
@ -1329,13 +1339,11 @@ double Group::ke(int igroup, Region *region)
|
||||
if (rmass) {
|
||||
for (int i = 0; i < nlocal; i++)
|
||||
if (mask[i] & groupbit && region->match(x[i][0], x[i][1], x[i][2]))
|
||||
one += (v[i][0]*v[i][0] + v[i][1]*v[i][1] + v[i][2]*v[i][2]) *
|
||||
rmass[i];
|
||||
one += (v[i][0] * v[i][0] + v[i][1] * v[i][1] + v[i][2] * v[i][2]) * rmass[i];
|
||||
} else {
|
||||
for (int i = 0; i < nlocal; i++)
|
||||
if (mask[i] & groupbit && region->match(x[i][0], x[i][1], x[i][2]))
|
||||
one += (v[i][0]*v[i][0] + v[i][1]*v[i][1] + v[i][2]*v[i][2]) *
|
||||
mass[type[i]];
|
||||
one += (v[i][0] * v[i][0] + v[i][1] * v[i][1] + v[i][2] * v[i][2]) * mass[type[i]];
|
||||
}
|
||||
|
||||
double all;
|
||||
@ -1372,8 +1380,10 @@ double Group::gyration(int igroup, double masstotal, double *cm)
|
||||
dx = unwrap[0] - cm[0];
|
||||
dy = unwrap[1] - cm[1];
|
||||
dz = unwrap[2] - cm[2];
|
||||
if (rmass) massone = rmass[i];
|
||||
else massone = mass[type[i]];
|
||||
if (rmass)
|
||||
massone = rmass[i];
|
||||
else
|
||||
massone = mass[type[i]];
|
||||
rg += (dx * dx + dy * dy + dz * dz) * massone;
|
||||
}
|
||||
double rg_all;
|
||||
@ -1412,8 +1422,10 @@ double Group::gyration(int igroup, double masstotal, double *cm, Region *region)
|
||||
dx = unwrap[0] - cm[0];
|
||||
dy = unwrap[1] - cm[1];
|
||||
dz = unwrap[2] - cm[2];
|
||||
if (rmass) massone = rmass[i];
|
||||
else massone = mass[type[i]];
|
||||
if (rmass)
|
||||
massone = rmass[i];
|
||||
else
|
||||
massone = mass[type[i]];
|
||||
rg += (dx * dx + dy * dy + dz * dz) * massone;
|
||||
}
|
||||
double rg_all;
|
||||
@ -1454,8 +1466,10 @@ void Group::angmom(int igroup, double *cm, double *lmom)
|
||||
dx = unwrap[0] - cm[0];
|
||||
dy = unwrap[1] - cm[1];
|
||||
dz = unwrap[2] - cm[2];
|
||||
if (rmass) massone = rmass[i];
|
||||
else massone = mass[type[i]];
|
||||
if (rmass)
|
||||
massone = rmass[i];
|
||||
else
|
||||
massone = mass[type[i]];
|
||||
p[0] += massone * (dy * v[i][2] - dz * v[i][1]);
|
||||
p[1] += massone * (dz * v[i][0] - dx * v[i][2]);
|
||||
p[2] += massone * (dx * v[i][1] - dy * v[i][0]);
|
||||
@ -1496,8 +1510,10 @@ void Group::angmom(int igroup, double *cm, double *lmom, Region *region)
|
||||
dx = unwrap[0] - cm[0];
|
||||
dy = unwrap[1] - cm[1];
|
||||
dz = unwrap[2] - cm[2];
|
||||
if (rmass) massone = rmass[i];
|
||||
else massone = mass[type[i]];
|
||||
if (rmass)
|
||||
massone = rmass[i];
|
||||
else
|
||||
massone = mass[type[i]];
|
||||
p[0] += massone * (dy * v[i][2] - dz * v[i][1]);
|
||||
p[1] += massone * (dz * v[i][0] - dx * v[i][2]);
|
||||
p[2] += massone * (dx * v[i][1] - dy * v[i][0]);
|
||||
@ -1603,8 +1619,7 @@ void Group::inertia(int igroup, double *cm, double itensor[3][3])
|
||||
|
||||
double ione[3][3];
|
||||
for (i = 0; i < 3; i++)
|
||||
for (j = 0; j < 3; j++)
|
||||
ione[i][j] = 0.0;
|
||||
for (j = 0; j < 3; j++) ione[i][j] = 0.0;
|
||||
|
||||
for (i = 0; i < nlocal; i++)
|
||||
if (mask[i] & groupbit) {
|
||||
@ -1612,8 +1627,10 @@ void Group::inertia(int igroup, double *cm, double itensor[3][3])
|
||||
dx = unwrap[0] - cm[0];
|
||||
dy = unwrap[1] - cm[1];
|
||||
dz = unwrap[2] - cm[2];
|
||||
if (rmass) massone = rmass[i];
|
||||
else massone = mass[type[i]];
|
||||
if (rmass)
|
||||
massone = rmass[i];
|
||||
else
|
||||
massone = mass[type[i]];
|
||||
ione[0][0] += massone * (dy * dy + dz * dz);
|
||||
ione[1][1] += massone * (dx * dx + dz * dz);
|
||||
ione[2][2] += massone * (dx * dx + dy * dy);
|
||||
@ -1653,8 +1670,7 @@ void Group::inertia(int igroup, double *cm, double itensor[3][3], Region *region
|
||||
|
||||
double ione[3][3];
|
||||
for (i = 0; i < 3; i++)
|
||||
for (j = 0; j < 3; j++)
|
||||
ione[i][j] = 0.0;
|
||||
for (j = 0; j < 3; j++) ione[i][j] = 0.0;
|
||||
|
||||
for (i = 0; i < nlocal; i++)
|
||||
if (mask[i] & groupbit && region->match(x[i][0], x[i][1], x[i][2])) {
|
||||
@ -1662,8 +1678,10 @@ void Group::inertia(int igroup, double *cm, double itensor[3][3], Region *region
|
||||
dx = unwrap[0] - cm[0];
|
||||
dy = unwrap[1] - cm[1];
|
||||
dz = unwrap[2] - cm[2];
|
||||
if (rmass) massone = rmass[i];
|
||||
else massone = mass[type[i]];
|
||||
if (rmass)
|
||||
massone = rmass[i];
|
||||
else
|
||||
massone = mass[type[i]];
|
||||
ione[0][0] += massone * (dy * dy + dz * dz);
|
||||
ione[1][1] += massone * (dx * dx + dz * dz);
|
||||
ione[2][2] += massone * (dx * dx + dy * dy);
|
||||
@ -1693,8 +1711,7 @@ void Group::omega(double *angmom, double inertia[3][3], double *w)
|
||||
inertia[0][1] * inertia[1][2] * inertia[2][0] +
|
||||
inertia[0][2] * inertia[1][0] * inertia[2][1] -
|
||||
inertia[0][0] * inertia[1][2] * inertia[2][1] -
|
||||
inertia[0][1]*inertia[1][0]*inertia[2][2] -
|
||||
inertia[2][0]*inertia[1][1]*inertia[0][2];
|
||||
inertia[0][1] * inertia[1][0] * inertia[2][2] - inertia[2][0] * inertia[1][1] * inertia[0][2];
|
||||
|
||||
// non-singular I matrix
|
||||
// use L = Iw, inverting I to solve for w
|
||||
@ -1703,31 +1720,23 @@ void Group::omega(double *angmom, double inertia[3][3], double *w)
|
||||
if (determinant > EPSILON) {
|
||||
|
||||
inverse[0][0] = inertia[1][1] * inertia[2][2] - inertia[1][2] * inertia[2][1];
|
||||
inverse[0][1] = -(inertia[0][1]*inertia[2][2] -
|
||||
inertia[0][2]*inertia[2][1]);
|
||||
inverse[0][1] = -(inertia[0][1] * inertia[2][2] - inertia[0][2] * inertia[2][1]);
|
||||
inverse[0][2] = inertia[0][1] * inertia[1][2] - inertia[0][2] * inertia[1][1];
|
||||
|
||||
inverse[1][0] = -(inertia[1][0]*inertia[2][2] -
|
||||
inertia[1][2]*inertia[2][0]);
|
||||
inverse[1][0] = -(inertia[1][0] * inertia[2][2] - inertia[1][2] * inertia[2][0]);
|
||||
inverse[1][1] = inertia[0][0] * inertia[2][2] - inertia[0][2] * inertia[2][0];
|
||||
inverse[1][2] = -(inertia[0][0]*inertia[1][2] -
|
||||
inertia[0][2]*inertia[1][0]);
|
||||
inverse[1][2] = -(inertia[0][0] * inertia[1][2] - inertia[0][2] * inertia[1][0]);
|
||||
|
||||
inverse[2][0] = inertia[1][0] * inertia[2][1] - inertia[1][1] * inertia[2][0];
|
||||
inverse[2][1] = -(inertia[0][0]*inertia[2][1] -
|
||||
inertia[0][1]*inertia[2][0]);
|
||||
inverse[2][1] = -(inertia[0][0] * inertia[2][1] - inertia[0][1] * inertia[2][0]);
|
||||
inverse[2][2] = inertia[0][0] * inertia[1][1] - inertia[0][1] * inertia[1][0];
|
||||
|
||||
for (int i = 0; i < 3; i++)
|
||||
for (int j = 0; j < 3; j++)
|
||||
inverse[i][j] /= determinant;
|
||||
for (int j = 0; j < 3; j++) inverse[i][j] /= determinant;
|
||||
|
||||
w[0] = inverse[0][0]*angmom[0] + inverse[0][1]*angmom[1] +
|
||||
inverse[0][2]*angmom[2];
|
||||
w[1] = inverse[1][0]*angmom[0] + inverse[1][1]*angmom[1] +
|
||||
inverse[1][2]*angmom[2];
|
||||
w[2] = inverse[2][0]*angmom[0] + inverse[2][1]*angmom[1] +
|
||||
inverse[2][2]*angmom[2];
|
||||
w[0] = inverse[0][0] * angmom[0] + inverse[0][1] * angmom[1] + inverse[0][2] * angmom[2];
|
||||
w[1] = inverse[1][0] * angmom[0] + inverse[1][1] * angmom[1] + inverse[1][2] * angmom[2];
|
||||
w[2] = inverse[2][0] * angmom[0] + inverse[2][1] * angmom[1] + inverse[2][2] * angmom[2];
|
||||
|
||||
// handle (nearly) singular I matrix
|
||||
// typically due to 2-atom group or linear molecule
|
||||
|
||||
Reference in New Issue
Block a user