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.