angular velocities of water molecules

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 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.

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’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).

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.

Efrem Braun

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.

I think the only error here is the re-use of the i index, likely

b/c it was a cut/paste. Compute omega/chunk is a relatively

new compute.

The inertia matrix is symmetric, so the determinant should

be > 0.0. The check for 0 is just in case it is singular.

If that’s the case you have bigger problems - one of
the moments is 0.0, which is like 0 mass, and nothing

is resisting your omega getting arbitrarily large (assuming the chunk

wants to spin).

Steve

I think the only error here is the re-use of the i index, likely

I made a test and with a standard compliant compiler, this code is functioning, as it should.

I think the only error here is the re-use of the i index, likely

I made a test and with a standard compliant compiler, this code is functioning, as it should.

I was talking about the original code.

Steve

I think the only error here is the re-use of the i index, likely

I made a test and with a standard compliant compiler, this code is functioning, as it should.

I was talking about the original code.

Me, too. The value of the variable i remains unchanged before and after the if statement.

I realized what my issue was: since the ‘chunk’ I was running this ‘compute’ on is a diatomic molecule, the determinant is always supposed to be 0. LAMMPS sometimes calculated the determinant as 0, and sometimes as things like 1.2E-6 or -5.4E-8 (I suppose due to floating point / roundoff errors). If it was 0 or -5.4E-8, LAMMPS returned the inverse without dividing by the determinant, and if it was 1.2E-6, LAMMPS returned the inverse after dividing by the determinant. This of course gives very wacky omegas that can vary very much from one timestep to the next.

So perhaps there should be some kind of check that the determinant isn’t 0? I think that running this ‘compute’ on diatomic or linear molecules isn’t such a silly thing that other users won’t make a similar mistake in the future.

More preferably, it’d be nice to be able to compute omega for such a case. I might start to work on such an implementation, though I’m rather tempted to switch my system to a non-linear one since it doesn’t much matter for the particular thing I’m looking at. I might be interested in learning the linear algebra for its own sake, so I’ll let you know if I end up implementing it. See http://physics.stackexchange.com/questions/277255/calculating-rotational-kinetic-energy-of-a-diatomic-molecule if you’re interested in helping me get started.

Thanks very much for your help.

I realized what my issue was: since the 'chunk' I was running this 'compute'
on is a diatomic molecule, the determinant is always supposed to be 0.
LAMMPS sometimes calculated the determinant as 0, and sometimes as things
like 1.2E-6 or -5.4E-8 (I suppose due to floating point / roundoff errors).
If it was 0 or -5.4E-8, LAMMPS returned the inverse without dividing by the
determinant, and if it was 1.2E-6, LAMMPS returned the inverse after
dividing by the determinant. This of course gives very wacky omegas that can
vary very much from one timestep to the next.

So perhaps there should be some kind of check that the determinant isn't 0?

it is probably better to process the moments of inertia rather than
the determinant (it is too late at that step), as that would also
enable you to identify, which rotational axis is not available.

I think that running this 'compute' on diatomic or linear molecules isn't
such a silly thing that other users won't make a similar mistake in the
future.

at least the test should probably not be for exactly zero but with a
tolerance, say 1.0e-15 or even larger.

More preferably, it'd be nice to be able to compute omega for such a case. I
might start to work on such an implementation, though I'm rather tempted to
switch my system to a non-linear one since it doesn't much matter for the
particular thing I'm looking at. I might be interested in learning the
linear algebra for its own sake, so I'll let you know if I end up
implementing it. See
classical mechanics - calculating rotational kinetic energy of a diatomic molecule - Physics Stack Exchange
if you're interested in helping me get started.

you also may find some inspiration in the fix rigid code, which can
handle linear molecules and, for example, has to account for removed
degrees of freedom differently between the two cases.

axel.

p.s.: i've updated group and compute omega/chunk in my lammps-icms
branch. there also was a memory leak in compute omega/chunk, that was
fixed only last week.

I figured out how to do this, and I made the change to compute_bond_local.cpp instead of compute_omega_chunk.cpp since it made more sense to me. I made a pull request on Github with the changes.

I posted a 7Sep16 patch that should fix this for compute omega/chunk

and the variable group omega function. It should be general for any

linear set of atoms (1 or more). I think that’s the only case that triggers

a singular inertia matrix.

Steve