Hi Steve,
I've been working with 'compute omega/chunk,' and I was getting some rather
bizarre values. I checked compute_omega_chunk.cpp, and I think that there
are some issues there. The code I'm referring to reads:
if (determinant > 0.0)
for (int i = 0; i < 3; i++)
for (int j = 0; j < 3; j++)
inverse[i][j] /= determinant;
This code segment is missing braces, so the angular momentum matrix's
i don't think that braces are required here. this code block in itself
is syntactically and logically correct. while many people favor using
braces in such code constructs for additional clarity (i am one of
them), they are not part of the LAMMPS coding style.
inverse does not end up having its values divided by the matrix's
determinant. The 'i' and 'j' indices should also perhaps be changed to 'j'
and 'k' since this code segment is within a larger 'for' loop in which the
index 'i' is used.
that should also not be a problem (again, i support making this change
for code clarity). the "for scope'd" int i declaration will hide the
outer loop i variable and make inaccessible, but after this block goes
out of scope, the original i is visible and accessible again. this is
standard C++ scoping behavior (AFAIK, only some older visual c++
compilers do not follow the common interpretation of the standard
here).
I'm also not sure why the 'if' statement is there in the first place. I
think the inverse matrix should get divided by the determinant no matter the
sign of the determinant.
i think so, too. looking at the code, however, it is a bit surprising,
that this went unnoticed for at least 7 years.
I'm guessing the 'if' statement was supposed to
check if the determinant is equal to 0, and if so, spit out an error (this
error checking would be good for folks like me who were trying to calculate
the omega of a diatomic molecule without knowing that the moment of inertia
tensor is singular for a group of 2 particles).
you usually don't put if statements with error exits inside of compute
kernels, but rather return "empty" data, if this is a non-critical
component (i.e. you would exit, if your force computation has gone
ballistic, but not for analysis). in this case, like when your matrix
is singular, the inverse is set to zero.
The same code is in group.cpp, so I'm guessing that the error that was
originally there got lifted into compute_omega_chunk.cpp.
possibly. i would program this kind of segment a bit differently and
avoid the extra nested loop with possibly costly divisions, and have
something like this:
double invdet = ione[0][0]*ione[1][1]*ione[2][2] +
ione[0][1]*ione[1][2]*ione[2][0] + ione[0][2]*ione[1][0]*ione[2][1] -
ione[0][0]*ione[1][2]*ione[2][1] - ione[0][1]*ione[1][0]*ione[2][2] -
ione[2][0]*ione[1][1]*ione[0][2];
// avoid division by zero for singular matrix. inverse will be set
to zero matrix,
// instead and thus omega will become zero for such cases as well.
if (invdet != 0.0) invdet = 1.0/invdet;
inverse[0][0] = invdet*(ione[1][1]*ione[2][2] - ione[1][2]*ione[2][1]);
inverse[0][1] = -invdet*(ione[0][1]*ione[2][2] - ione[0][2]*ione[2][1]);
inverse[0][2] = invdet*(ione[0][1]*ione[1][2] - ione[0][2]*ione[1][1]);
inverse[1][0] = -invdet*(ione[1][0]*ione[2][2] - ione[1][2]*ione[2][0]);
inverse[1][1] = invdet*(ione[0][0]*ione[2][2] - ione[0][2]*ione[2][0]);
inverse[1][2] = -invdet*(ione[0][0]*ione[1][2] - ione[0][2]*ione[1][0]);
inverse[2][0] = invdet*(ione[1][0]*ione[2][1] - ione[1][1]*ione[2][0]);
inverse[2][1] = -invdet*(ione[0][0]*ione[2][1] - ione[0][1]*ione[2][0]);
inverse[2][2] = invdet*(ione[0][0]*ione[1][1] - ione[0][1]*ione[1][0]);
same for the code in group and then double check with the
documentation, and make sure that it mentions the restriction.
axel.