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

change angular velocity computation for groups/chunks to handle (near) zero determinants more consistently and be more efficient.

Browse Source
This commit is contained in:
Axel Kohlmeyer
2016-08-30 16:47:10 -04:00
parent 7875009218
commit 2af907a10d
2 changed files with 41 additions and 38 deletions
Download Patch File Download Diff File Expand all files Collapse all files

37
src/compute_omega_chunk.cpp
View File

@ -23,6 +23,8 @@
using namespace LAMMPS_NS;
#define SMALL 1.0e-15
/* ---------------------------------------------------------------------- */
ComputeOmegaChunk::ComputeOmegaChunk(LAMMPS *lmp, int narg, char **arg) :
@ -207,27 +209,28 @@ void ComputeOmegaChunk::compute_array()
ione[2][1] = ione[1][2];
ione[2][0] = ione[0][2];
inverse[0][0] = ione[1][1]*ione[2][2] - ione[1][2]*ione[2][1];
inverse[0][1] = -(ione[0][1]*ione[2][2] - ione[0][2]*ione[2][1]);
inverse[0][2] = ione[0][1]*ione[1][2] - ione[0][2]*ione[1][1];
inverse[1][0] = -(ione[1][0]*ione[2][2] - ione[1][2]*ione[2][0]);
inverse[1][1] = ione[0][0]*ione[2][2] - ione[0][2]*ione[2][0];
inverse[1][2] = -(ione[0][0]*ione[1][2] - ione[0][2]*ione[1][0]);
inverse[2][0] = ione[1][0]*ione[2][1] - ione[1][1]*ione[2][0];
inverse[2][1] = -(ione[0][0]*ione[2][1] - ione[0][1]*ione[2][0]);
inverse[2][2] = ione[0][0]*ione[1][1] - ione[0][1]*ione[1][0];
double determinant = ione[0][0]*ione[1][1]*ione[2][2] +
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];
if (determinant > 0.0)
for (int i = 0; i < 3; i++)
for (int j = 0; j < 3; j++)
inverse[i][j] /= determinant;
// avoid division by (near) zero for (near) singular matrix. inverse will be set to zero matrix instead.
if (fabs(invdet) < SMALL)
invdet = 1.0/invdet;
else
invdet = 0.0;
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]);
omega[i][0] = inverse[0][0]*angmomall[i][0] +
inverse[0][1]*angmomall[i][1] + inverse[0][2]*angmomall[i][2];

42
src/group.cpp
View File

@ -42,6 +42,7 @@ enum{TYPE,MOLECULE,ID};
enum{LT,LE,GT,GE,EQ,NEQ,BETWEEN};
#define BIG 1.0e20
#define SMALL 1.0e-15
// allocate space for static class variable
@ -1670,29 +1671,28 @@ void Group::omega(double *angmom, double inertia[3][3], double *w)
{
double inverse[3][3];
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][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][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[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][2] = inertia[0][0]*inertia[1][1] - inertia[0][1]*inertia[1][0];
double determinant = inertia[0][0]*inertia[1][1]*inertia[2][2] +
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] -
double invdet = inertia[0][0]*inertia[1][1]*inertia[2][2] +
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];
if (determinant > 0.0)
for (int i = 0; i < 3; i++)
for (int j = 0; j < 3; j++)
inverse[i][j] /= determinant;
// avoid division by (near) zero for (near) singular matrix. inverse will be set to zero matrix instead.
if (fabs(invdet) < SMALL)
invdet = 1.0/invdet;
else
invdet = 0.0;
inverse[0][0] = invdet*(inertia[1][1]*inertia[2][2] - inertia[1][2]*inertia[2][1]);
inverse[0][1] = -invdet*(inertia[0][1]*inertia[2][2] - inertia[0][2]*inertia[2][1]);
inverse[0][2] = invdet*(inertia[0][1]*inertia[1][2] - inertia[0][2]*inertia[1][1]);
inverse[1][0] = -invdet*(inertia[1][0]*inertia[2][2] - inertia[1][2]*inertia[2][0]);
inverse[1][1] = invdet*(inertia[0][0]*inertia[2][2] - inertia[0][2]*inertia[2][0]);
inverse[1][2] = -invdet*(inertia[0][0]*inertia[1][2] - inertia[0][2]*inertia[1][0]);
inverse[2][0] = invdet*(inertia[1][0]*inertia[2][1] - inertia[1][1]*inertia[2][0]);
inverse[2][1] = -invdet*(inertia[0][0]*inertia[2][1] - inertia[0][1]*inertia[2][0]);
inverse[2][2] = invdet*(inertia[0][0]*inertia[1][1] - inertia[0][1]*inertia[1][0]);
w[0] = inverse[0][0]*angmom[0] + inverse[0][1]*angmom[1] +
inverse[0][2]*angmom[2];