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@12866 f3b2605a-c512-4ea7-a41b-209d697bcdaa

Browse Source
This commit is contained in:
sjplimp
2015-01-05 17:15:35 +00:00
parent 626ae3f693
commit 7b8f7ed44d
2 changed files with 50 additions and 0 deletions
Download Patch File Download Diff File Expand all files Collapse all files

49
src/math_extra.cpp
View File

@ -226,6 +226,55 @@ void richardson(double *q, double *m, double *w, double *moments, double dtq)
MathExtra::qnormalize(q);
}
/* ----------------------------------------------------------------------
apply evolution operators to quat, quat momentum
Miller et al., J Chem Phys. 116, 8649-8659 (2002)
------------------------------------------------------------------------- */
void no_squish_rotate(int k, double *p, double *q, double *inertia,
double dt)
{
double phi,c_phi,s_phi,kp[4],kq[4];
// apply permuation operator on p and q, get kp and kq
if (k == 1) {
kq[0] = -q[1]; kp[0] = -p[1];
kq[1] = q[0]; kp[1] = p[0];
kq[2] = q[3]; kp[2] = p[3];
kq[3] = -q[2]; kp[3] = -p[2];
} else if (k == 2) {
kq[0] = -q[2]; kp[0] = -p[2];
kq[1] = -q[3]; kp[1] = -p[3];
kq[2] = q[0]; kp[2] = p[0];
kq[3] = q[1]; kp[3] = p[1];
} else if (k == 3) {
kq[0] = -q[3]; kp[0] = -p[3];
kq[1] = q[2]; kp[1] = p[2];
kq[2] = -q[1]; kp[2] = -p[1];
kq[3] = q[0]; kp[3] = p[0];
}
// obtain phi, cosines and sines
phi = p[0]*kq[0] + p[1]*kq[1] + p[2]*kq[2] + p[3]*kq[3];
if (fabs(inertia[k-1]) < 1e-6) phi *= 0.0;
else phi /= 4.0 * inertia[k-1];
c_phi = cos(dt * phi);
s_phi = sin(dt * phi);
// advance p and q
p[0] = c_phi*p[0] + s_phi*kp[0];
p[1] = c_phi*p[1] + s_phi*kp[1];
p[2] = c_phi*p[2] + s_phi*kp[2];
p[3] = c_phi*p[3] + s_phi*kp[3];
q[0] = c_phi*q[0] + s_phi*kq[0];
q[1] = c_phi*q[1] + s_phi*kq[1];
q[2] = c_phi*q[2] + s_phi*kq[2];
q[3] = c_phi*q[3] + s_phi*kq[3];
}
/* ----------------------------------------------------------------------
compute omega from angular momentum, both in space frame
only know Idiag so need to do M = Iw in body frame