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

git-svn-id: svn://svn.icms.temple.edu/lammps-ro/trunk@15557 f3b2605a-c512-4ea7-a41b-209d697bcdaa

Browse Source
This commit is contained in:
sjplimp
2016-09-06 23:19:15 +00:00
parent 39be4185c4
commit bf0c18a0f2
1 changed files with 77 additions and 30 deletions
Download Patch File Download Diff File Expand all files Collapse all files

99
src/group.cpp
View File

@ -1671,39 +1671,86 @@ void Group::inertia(int igroup, double *cm, double itensor[3][3], int iregion)
void Group::omega(double *angmom, double inertia[3][3], double *w) void Group::omega(double *angmom, double inertia[3][3], double *w)
{ {
double idiag[3],ex[3],ey[3],ez[3],cross[3]; double idiag[3],ex[3],ey[3],ez[3],cross[3];
double evectors[3][3]; double evectors[3][3],inverse[3][3];
int ierror = MathExtra::jacobi(inertia,idiag,evectors); // determinant = triple product of rows of inertia matrix
if (ierror) error->all(FLERR,
"Insufficient Jacobi rotations for group::omega");
ex[0] = evectors[0][0]; double determinant = inertia[0][0]*inertia[1][1]*inertia[2][2] +
ex[1] = evectors[1][0]; inertia[0][1]*inertia[1][2]*inertia[2][0] +
ex[2] = evectors[2][0]; inertia[0][2]*inertia[1][0]*inertia[2][1] -
ey[0] = evectors[0][1]; inertia[0][0]*inertia[1][2]*inertia[2][1] -
ey[1] = evectors[1][1]; inertia[0][1]*inertia[1][0]*inertia[2][2] -
ey[2] = evectors[2][1]; inertia[2][0]*inertia[1][1]*inertia[0][2];
ez[0] = evectors[0][2];
ez[1] = evectors[1][2];
ez[2] = evectors[2][2];
// enforce 3 evectors as a right-handed coordinate system // non-singular I matrix
// flip 3rd vector if needed // use L = Iw, inverting I to solve for w
MathExtra::cross3(ex,ey,cross); if (determinant > EPSILON) {
if (MathExtra::dot3(cross,ez) < 0.0) MathExtra::negate3(ez);
// if any principal moment < scaled EPSILON, set to 0.0 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];
double max; inverse[1][0] = -(inertia[1][0]*inertia[2][2] -
max = MAX(idiag[0],idiag[1]); inertia[1][2]*inertia[2][0]);
max = MAX(max,idiag[2]); 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]);
if (idiag[0] < EPSILON*max) idiag[0] = 0.0; inverse[2][0] = inertia[1][0]*inertia[2][1] - inertia[1][1]*inertia[2][0];
if (idiag[1] < EPSILON*max) idiag[1] = 0.0; inverse[2][1] = -(inertia[0][0]*inertia[2][1] -
if (idiag[2] < EPSILON*max) idiag[2] = 0.0; inertia[0][1]*inertia[2][0]);
inverse[2][2] = inertia[0][0]*inertia[1][1] - inertia[0][1]*inertia[1][0];
// calculate omega using diagonalized inertia matrix for (int i = 0; i < 3; i++)
for (int j = 0; j < 3; j++)
inverse[i][j] /= determinant;
MathExtra::angmom_to_omega(angmom,ex,ey,ez,idiag,w); 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];
// handle (nearly) singular I matrix
// due to 2-atom group or linear molecule
// use jacobi() and angmom_to_omega() to calculate valid omega
} else {
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);
}
} }