From 56945a56aabe0d97aa33e9bc87539d1170161fe7 Mon Sep 17 00:00:00 2001
From: sjplimp <sjplimp@f3b2605a-c512-4ea7-a41b-209d697bcdaa>
Date: Tue, 6 Sep 2016 21:55:39 +0000
Subject: [PATCH] git-svn-id: svn://svn.icms.temple.edu/lammps-ro/trunk@15547
f3b2605a-c512-4ea7-a41b-209d697bcdaa
---
src/compute_omega_chunk.cpp | 127 +++++++++++++++++++++++++-----------
src/group.cpp | 73 +++++++++++----------
2 files changed, 129 insertions(+), 71 deletions(-)
diff --git a/src/compute_omega_chunk.cpp b/src/compute_omega_chunk.cpp
index ca4af7931f..f2a8394964 100644
--- a/src/compute_omega_chunk.cpp
+++ b/src/compute_omega_chunk.cpp
@@ -18,11 +18,14 @@
#include "modify.h"
#include "compute_chunk_atom.h"
#include "domain.h"
+#include "math_extra.h"
#include "memory.h"
#include "error.h"
using namespace LAMMPS_NS;
+#define EPSILON 1.0e-6
+
/* ---------------------------------------------------------------------- */
ComputeOmegaChunk::ComputeOmegaChunk(LAMMPS *lmp, int narg, char **arg) :
@@ -192,49 +195,97 @@ void ComputeOmegaChunk::compute_array()
MPI_Allreduce(&angmom[0][0],&angmomall[0][0],3*nchunk,
MPI_DOUBLE,MPI_SUM,world);
- // compute omega for each chunk from L = Iw, inverting I to solve for w
+ // compute omega for each chunk
- double ione[3][3],inverse[3][3];
+ double determinant,invdeterminant;
+ double idiag[3],ex[3],ey[3],ez[3],cross[3];
+ double ione[3][3],inverse[3][3],evectors[3][3];
+ double *iall,*mall;
for (m = 0; m < nchunk; m++) {
- ione[0][0] = inertiaall[m][0];
- ione[1][1] = inertiaall[m][1];
- ione[2][2] = inertiaall[m][2];
- ione[0][1] = inertiaall[m][3];
- ione[1][2] = inertiaall[m][4];
- ione[0][2] = inertiaall[m][5];
- ione[1][0] = ione[0][1];
- 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] +
- 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)
+
+ // determinant = triple product of rows of inertia matrix
+
+ iall = &inertiaall[m][0];
+ determinant = iall[0] * (iall[1]*iall[2] - iall[4]*iall[4]) +
+ iall[3] * (iall[4]*iall[5] - iall[3]*iall[2]) +
+ iall[5] * (iall[3]*iall[4] - iall[1]*iall[5]);
+
+ ione[0][0] = iall[0];
+ ione[1][1] = iall[1];
+ ione[2][2] = iall[2];
+ ione[0][1] = ione[1][0] = iall[3];
+ ione[1][2] = ione[2][1] = iall[4];
+ ione[0][2] = ione[2][0] = iall[5];
+
+ // non-singular I matrix
+ // use L = Iw, inverting I to solve for w
+
+ if (determinant > EPSILON) {
+ 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];
+
+ invdeterminant = 1.0/determinant;
for (i = 0; i < 3; i++)
for (j = 0; j < 3; j++)
- inverse[i][j] /= determinant;
-
- omega[m][0] = inverse[0][0]*angmom[m][0] + inverse[0][1]*angmom[m][1] +
- inverse[0][2]*angmom[m][2];
- omega[m][1] = inverse[1][0]*angmom[m][0] + inverse[1][1]*angmom[m][1] +
- inverse[1][2]*angmom[m][2];
- omega[m][2] = inverse[2][0]*angmom[m][0] + inverse[2][1]*angmom[m][1] +
- inverse[2][2]*angmom[i][2];
+ inverse[i][j] *= invdeterminant;
+
+ mall = &angmomall[m][0];
+ omega[m][0] = inverse[0][0]*mall[0] + inverse[0][1]*mall[1] +
+ inverse[0][2]*mall[2];
+ omega[m][1] = inverse[1][0]*mall[0] + inverse[1][1]*mall[1] +
+ inverse[1][2]*mall[2];
+ omega[m][2] = inverse[2][0]*mall[0] + inverse[2][1]*mall[1] +
+ inverse[2][2]*mall[2];
+
+ // handle each (nearly) singular I matrix
+ // due to 2-atom chunk or linear molecule
+ // use jacobi() and angmom_to_omega() to calculate valid omega
+
+ } else {
+ int ierror = MathExtra::jacobi(ione,idiag,evectors);
+ if (ierror) error->all(FLERR,
+ "Insufficient Jacobi rotations for omega/chunk");
+
+ ex[0] = evectors[0][0];
+ ex[1] = evectors[1][0];
+ ex[2] = evectors[2][0];
+ ey[0] = evectors[0][1];
+ ey[1] = evectors[1][1];
+ ey[2] = evectors[2][1];
+ ez[0] = evectors[0][2];
+ ez[1] = evectors[1][2];
+ ez[2] = evectors[2][2];
+
+ // enforce 3 evectors as a right-handed coordinate system
+ // flip 3rd vector if needed
+
+ MathExtra::cross3(ex,ey,cross);
+ if (MathExtra::dot3(cross,ez) < 0.0) MathExtra::negate3(ez);
+
+ // if any principal moment < scaled EPSILON, set to 0.0
+
+ double max;
+ max = MAX(idiag[0],idiag[1]);
+ max = MAX(max,idiag[2]);
+
+ if (idiag[0] < EPSILON*max) idiag[0] = 0.0;
+ if (idiag[1] < EPSILON*max) idiag[1] = 0.0;
+ if (idiag[2] < EPSILON*max) idiag[2] = 0.0;
+
+ // calculate omega using diagonalized inertia matrix
+
+ MathExtra::angmom_to_omega(&angmomall[m][0],ex,ey,ez,idiag,&omega[m][0]);
+ }
}
}
diff --git a/src/group.cpp b/src/group.cpp
index 5b572b2070..da0e94fc11 100644
--- a/src/group.cpp
+++ b/src/group.cpp
@@ -29,6 +29,7 @@
#include "input.h"
#include "variable.h"
#include "dump.h"
+#include "math_extra.h"
#include "memory.h"
#include "error.h"
@@ -37,6 +38,7 @@
using namespace LAMMPS_NS;
#define MAX_GROUP 32
+#define EPSILON 1.0e-6
enum{TYPE,MOLECULE,ID};
enum{LT,LE,GT,GE,EQ,NEQ,BETWEEN};
@@ -1662,42 +1664,47 @@ void Group::inertia(int igroup, double *cm, double itensor[3][3], int iregion)
}
/* ----------------------------------------------------------------------
- compute angular velocity omega from L = Iw, inverting I to solve for w
+ compute angular velocity omega from L and I
really not a group/region operation, but L,I were computed for a group/region
+ diagonalize I instead of inverting it, to allow for a singular matrix
------------------------------------------------------------------------- */
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] -
- 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;
-
- 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];
+ double idiag[3],ex[3],ey[3],ez[3],cross[3];
+ double evectors[3][3];
+
+ int ierror = MathExtra::jacobi(inertia,idiag,evectors);
+ if (ierror) error->all(FLERR,
+ "Insufficient Jacobi rotations for group::omega");
+
+ ex[0] = evectors[0][0];
+ ex[1] = evectors[1][0];
+ ex[2] = evectors[2][0];
+ ey[0] = evectors[0][1];
+ ey[1] = evectors[1][1];
+ ey[2] = evectors[2][1];
+ ez[0] = evectors[0][2];
+ ez[1] = evectors[1][2];
+ ez[2] = evectors[2][2];
+
+ // enforce 3 evectors as a right-handed coordinate system
+ // flip 3rd vector if needed
+
+ MathExtra::cross3(ex,ey,cross);
+ if (MathExtra::dot3(cross,ez) < 0.0) MathExtra::negate3(ez);
+
+ // if any principal moment < scaled EPSILON, set to 0.0
+
+ double max;
+ max = MAX(idiag[0],idiag[1]);
+ max = MAX(max,idiag[2]);
+
+ if (idiag[0] < EPSILON*max) idiag[0] = 0.0;
+ if (idiag[1] < EPSILON*max) idiag[1] = 0.0;
+ if (idiag[2] < EPSILON*max) idiag[2] = 0.0;
+
+ // calculate omega using diagonalized inertia matrix
+
+ MathExtra::angmom_to_omega(angmom,ex,ey,ez,idiag,w);
}
--
GitLab