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

Browse Source
This commit is contained in:
sjplimp
2011-04-19 00:09:41 +00:00
parent 6284969800
commit fb1e5c2389
12 changed files with 132 additions and 509 deletions
Download Patch File Download Diff File Expand all files Collapse all files

12
src/ASPHERE/fix_nh_asphere.cpp
View File

@ -71,14 +71,14 @@ void FixNHAsphere::richardson(double *q, double *m, double *moments)
// full update from dq/dt = 1/2 w q
double wq[4];
MathExtra::multiply_vec_quat(w,q,wq);
MathExtra::vecquat(w,q,wq);
double qfull[4];
qfull[0] = q[0] + dtq * wq[0];
qfull[1] = q[1] + dtq * wq[1];
qfull[2] = q[2] + dtq * wq[2];
qfull[3] = q[3] + dtq * wq[3];
MathExtra::normalize4(qfull);
MathExtra::qnormalize(qfull);
// 1st half of update from dq/dt = 1/2 w q
@ -87,13 +87,13 @@ void FixNHAsphere::richardson(double *q, double *m, double *moments)
qhalf[1] = q[1] + 0.5*dtq * wq[1];
qhalf[2] = q[2] + 0.5*dtq * wq[2];
qhalf[3] = q[3] + 0.5*dtq * wq[3];
MathExtra::normalize4(qhalf);
MathExtra::qnormalize(qhalf);
// re-compute omega at 1/2 step from m at 1/2 step and q at 1/2 step
// recompute wq
omega_from_mq(qhalf,m,moments,w);
MathExtra::multiply_vec_quat(w,qhalf,wq);
MathExtra::vecquat(w,qhalf,wq);
// 2nd half of update from dq/dt = 1/2 w q
@ -101,7 +101,7 @@ void FixNHAsphere::richardson(double *q, double *m, double *moments)
qhalf[1] += 0.5*dtq * wq[1];
qhalf[2] += 0.5*dtq * wq[2];
qhalf[3] += 0.5*dtq * wq[3];
MathExtra::normalize4(qhalf);
MathExtra::qnormalize(qhalf);
// corrected Richardson update
@ -109,7 +109,7 @@ void FixNHAsphere::richardson(double *q, double *m, double *moments)
q[1] = 2.0*qhalf[1] - qfull[1];
q[2] = 2.0*qhalf[2] - qfull[2];
q[3] = 2.0*qhalf[3] - qfull[3];
MathExtra::normalize4(q);
MathExtra::qnormalize(q);
}
/* ----------------------------------------------------------------------

12
src/ASPHERE/fix_nve_asphere.cpp
View File

@ -155,14 +155,14 @@ void FixNVEAsphere::richardson(double *q, double *m, double *moments)
// full update from dq/dt = 1/2 w q
double wq[4];
MathExtra::multiply_vec_quat(w,q,wq);
MathExtra::vecquat(w,q,wq);
double qfull[4];
qfull[0] = q[0] + dtq * wq[0];
qfull[1] = q[1] + dtq * wq[1];
qfull[2] = q[2] + dtq * wq[2];
qfull[3] = q[3] + dtq * wq[3];
MathExtra::normalize4(qfull);
MathExtra::qnormalize(qfull);
// 1st half of update from dq/dt = 1/2 w q
@ -171,13 +171,13 @@ void FixNVEAsphere::richardson(double *q, double *m, double *moments)
qhalf[1] = q[1] + 0.5*dtq * wq[1];
qhalf[2] = q[2] + 0.5*dtq * wq[2];
qhalf[3] = q[3] + 0.5*dtq * wq[3];
MathExtra::normalize4(qhalf);
MathExtra::qnormalize(qhalf);
// re-compute omega at 1/2 step from m at 1/2 step and q at 1/2 step
// recompute wq
omega_from_mq(qhalf,m,moments,w);
MathExtra::multiply_vec_quat(w,qhalf,wq);
MathExtra::vecquat(w,qhalf,wq);
// 2nd half of update from dq/dt = 1/2 w q
@ -185,7 +185,7 @@ void FixNVEAsphere::richardson(double *q, double *m, double *moments)
qhalf[1] += 0.5*dtq * wq[1];
qhalf[2] += 0.5*dtq * wq[2];
qhalf[3] += 0.5*dtq * wq[3];
MathExtra::normalize4(qhalf);
MathExtra::qnormalize(qhalf);
// corrected Richardson update
@ -193,7 +193,7 @@ void FixNVEAsphere::richardson(double *q, double *m, double *moments)
q[1] = 2.0*qhalf[1] - qfull[1];
q[2] = 2.0*qhalf[2] - qfull[2];
q[3] = 2.0*qhalf[3] - qfull[3];
MathExtra::normalize4(q);
MathExtra::qnormalize(q);
}
/* ----------------------------------------------------------------------

3
src/ASPHERE/pair_gayberne.cpp
View File

@ -43,8 +43,7 @@ enum{SPHERE_SPHERE,SPHERE_ELLIPSE,ELLIPSE_SPHERE,ELLIPSE_ELLIPSE};
PairGayBerne::PairGayBerne(LAMMPS *lmp) : Pair(lmp)
{
avec = (AtomVecEllipsoid *) atom->style_match("ellipsoid");
if (!avec)
error->all("Pair gayberne requires atom style ellipsoid");
if (!avec) error->all("Pair gayberne requires atom style ellipsoid");
single_enable = 0;
}

3
src/ASPHERE/pair_resquared.cpp
View File

@ -45,8 +45,7 @@ PairRESquared::PairRESquared(LAMMPS *lmp) : Pair(lmp),
cr60(pow(60.0,1.0/3.0))
{
avec = (AtomVecEllipsoid *) atom->style_match("ellipsoid");
if (!avec)
error->all("Pair gayberne requires atom style ellipsoid");
if (!avec) error->all("Pair resquared requires atom style ellipsoid");
single_enable = 0;

9
src/SRD/fix_srd.cpp
View File

@ -2092,7 +2092,8 @@ void FixSRD::parameterize()
// big particle must either have radius > 0 or shape > 0 defined
// apply radfactor at end
AtomVecEllipsoid::Bonus *ebonus = avec_ellipsoid->bonus;
AtomVecEllipsoid::Bonus *ebonus;
if (avec_ellipsoid) ebonus = avec_ellipsoid->bonus;
double *radius = atom->radius;
int *ellipsoid = atom->ellipsoid;
int *mask = atom->mask;
@ -2406,7 +2407,8 @@ void FixSRD::big_static()
double rad,arad,brad,crad;
double *shape;
AtomVecEllipsoid::Bonus *ebonus = avec_ellipsoid->bonus;
AtomVecEllipsoid::Bonus *ebonus;
if (avec_ellipsoid) ebonus = avec_ellipsoid->bonus;
double *radius = atom->radius;
int *ellipsoid = atom->ellipsoid;
int *type = atom->type;
@ -2449,7 +2451,8 @@ void FixSRD::big_dynamic()
int i;
double *shape,*quat;
AtomVecEllipsoid::Bonus *ebonus = avec_ellipsoid->bonus;
AtomVecEllipsoid::Bonus *ebonus;
if (avec_ellipsoid) ebonus = avec_ellipsoid->bonus;
double **omega = atom->omega;
double **angmom = atom->angmom;
double *rmass = atom->rmass;

2
src/atom_vec_ellipsoid.cpp
View File

@ -1198,7 +1198,7 @@ void AtomVecEllipsoid::data_atom_bonus(int m, char **values)
quat[1] = atof(values[4]);
quat[2] = atof(values[5]);
quat[3] = atof(values[6]);
MathExtra::normalize4(quat);
MathExtra::qnormalize(quat);
// reset ellipsoid mass
// previously stored density in rmass

420
src/fix_rigid.cpp
View File

@ -34,7 +34,6 @@ using namespace LAMMPS_NS;
#define TOLERANCE 1.0e-6
#define EPSILON 1.0e-7
#define MAXJACOBI 50
#define MIN(A,B) ((A) < (B)) ? (A) : (B)
#define MAX(A,B) ((A) > (B)) ? (A) : (B)
@ -470,7 +469,8 @@ void FixRigid::init()
extended = dorientflag = qorientflag = 0;
AtomVecEllipsoid::Bonus *ebonus = avec_ellipsoid->bonus;
AtomVecEllipsoid::Bonus *ebonus;
if (avec_ellipsoid) ebonus = avec_ellipsoid->bonus;
double **mu = atom->mu;
double *radius = atom->radius;
double *rmass = atom->rmass;
@ -589,6 +589,7 @@ void FixRigid::init()
pre_neighbor();
// compute 6 moments of inertia of each body
// stored as 6-vector in Voigt ordering
// dx,dy,dz = coords relative to center-of-mass
double dx,dy,dz,rad;
@ -632,8 +633,7 @@ void FixRigid::init()
// extended particles may contribute extra terms to moments of inertia
if (extended) {
double ex[3],ey[3],ez[3],idiag[3];
double p[3][3],ptrans[3][3],ispace[3][3],itemp[3][3];
double ivec[6];
double *shape,*quatatom;
for (i = 0; i < nlocal; i++) {
@ -652,19 +652,13 @@ void FixRigid::init()
if (eflags[i] & INERTIA_ELLIPSOID) {
shape = ebonus[ellipsoid[i]].shape;
quatatom = ebonus[ellipsoid[i]].quat;
MathExtra::quat_to_mat(quatatom,p);
MathExtra::quat_to_mat_trans(quatatom,ptrans);
idiag[0] = 0.2*massone * (shape[1]*shape[1]+shape[2]*shape[2]);
idiag[1] = 0.2*massone * (shape[0]*shape[0]+shape[2]*shape[2]);
idiag[2] = 0.2*massone * (shape[0]*shape[0]+shape[1]*shape[1]);
MathExtra::diag_times3(idiag,ptrans,itemp);
MathExtra::times3(p,itemp,ispace);
sum[ibody][0] += ispace[0][0];
sum[ibody][1] += ispace[1][1];
sum[ibody][2] += ispace[2][2];
sum[ibody][3] += ispace[0][1];
sum[ibody][4] += ispace[1][2];
sum[ibody][5] += ispace[0][2];
MathExtra::inertia_ellipsoid(shape,quatatom,massone,ivec);
sum[ibody][0] += ivec[0];
sum[ibody][1] += ivec[1];
sum[ibody][2] += ivec[2];
sum[ibody][3] += ivec[3];
sum[ibody][4] += ivec[4];
sum[ibody][5] += ivec[5];
}
}
}
@ -674,12 +668,9 @@ void FixRigid::init()
// inertia = 3 eigenvalues = principal moments of inertia
// ex_space,ey_space,ez_space = 3 eigenvectors = principal axes of rigid body
double **tensor,**evectors;
memory->create(tensor,3,3,"fix_rigid:tensor");
memory->create(evectors,3,3,"fix_rigid:evectors");
int ierror;
double tensor[3][3],evectors[3][3];
double ez0,ez1,ez2;
int ierror;
for (ibody = 0; ibody < nbody; ibody++) {
tensor[0][0] = all[ibody][0];
@ -689,7 +680,7 @@ void FixRigid::init()
tensor[1][2] = tensor[2][1] = all[ibody][4];
tensor[0][2] = tensor[2][0] = all[ibody][5];
ierror = jacobi(tensor,inertia[ibody],evectors);
ierror = MathExtra::jacobi(tensor,inertia[ibody],evectors);
if (ierror) error->all("Insufficient Jacobi rotations for rigid body");
ex_space[ibody][0] = evectors[0][0];
@ -733,14 +724,10 @@ void FixRigid::init()
// create initial quaternion
q_from_exyz(ex_space[ibody],ey_space[ibody],ez_space[ibody],quat[ibody]);
MathExtra::exyz_to_q(ex_space[ibody],ey_space[ibody],ez_space[ibody],
quat[ibody]);
}
// free temporary memory
memory->destroy(tensor);
memory->destroy(evectors);
// displace = initial atom coords in basis of principal axes
// set displace = 0.0 for atoms not in any rigid body
// for extended particles, set their orientation wrt to rigid body
@ -795,9 +782,9 @@ void FixRigid::init()
if (eflags[i] & ORIENT_QUAT) {
quatatom = ebonus[ellipsoid[i]].quat;
qconjugate(quat[ibody],qc);
quatquat(qc,quatatom,qorient[i]);
qnormalize(qorient[i]);
MathExtra::qconjugate(quat[ibody],qc);
MathExtra::quatquat(qc,quatatom,qorient[i]);
MathExtra::qnormalize(qorient[i]);
} else if (qorientflag)
qorient[i][0] = qorient[i][1] = qorient[i][2] = qorient[i][3] = 0.0;
}
@ -805,6 +792,7 @@ void FixRigid::init()
// test for valid principal moments & axes
// recompute moments of inertia around new axes
// stored as 6-vector in Voigt ordering
// 3 diagonal moments should equal principal moments
// 3 off-diagonal moments should be 0.0
// extended particles may contribute extra terms to moments of inertia
@ -830,8 +818,7 @@ void FixRigid::init()
}
if (extended) {
double ex[3],ey[3],ez[3],idiag[3];
double p[3][3],ptrans[3][3],ispace[3][3],itemp[3][3];
double ivec[6];
double *shape;
for (i = 0; i < nlocal; i++) {
@ -849,19 +836,13 @@ void FixRigid::init()
}
if (eflags[i] & INERTIA_ELLIPSOID) {
shape = ebonus[ellipsoid[i]].shape;
MathExtra::quat_to_mat(qorient[i],p);
MathExtra::quat_to_mat_trans(qorient[i],ptrans);
idiag[0] = 0.2*massone * (shape[1]*shape[1]+shape[2]*shape[2]);
idiag[1] = 0.2*massone * (shape[0]*shape[0]+shape[2]*shape[2]);
idiag[2] = 0.2*massone * (shape[0]*shape[0]+shape[1]*shape[1]);
MathExtra::diag_times3(idiag,ptrans,itemp);
MathExtra::times3(p,itemp,ispace);
sum[ibody][0] += ispace[0][0];
sum[ibody][1] += ispace[1][1];
sum[ibody][2] += ispace[2][2];
sum[ibody][3] += ispace[0][1];
sum[ibody][4] += ispace[1][2];
sum[ibody][5] += ispace[0][2];
MathExtra::inertia_ellipsoid(shape,qorient[i],massone,ivec);
sum[ibody][0] += ivec[0];
sum[ibody][1] += ivec[1];
sum[ibody][2] += ivec[2];
sum[ibody][3] += ivec[3];
sum[ibody][4] += ivec[4];
sum[ibody][5] += ivec[5];
}
}
}
@ -1048,7 +1029,7 @@ void FixRigid::setup(int vflag)
// set velocities from angmom & omega
for (ibody = 0; ibody < nbody; ibody++)
omega_from_angmom(angmom[ibody],ex_space[ibody],ey_space[ibody],
MathExtra::angmom_to_omega(angmom[ibody],ex_space[ibody],ey_space[ibody],
ez_space[ibody],inertia[ibody],omega[ibody]);
set_v();
@ -1209,7 +1190,7 @@ void FixRigid::final_integrate()
angmom[ibody][1] += dtf * torque[ibody][1] * tflag[ibody][1];
angmom[ibody][2] += dtf * torque[ibody][2] * tflag[ibody][2];
omega_from_angmom(angmom[ibody],ex_space[ibody],ey_space[ibody],
MathExtra::angmom_to_omega(angmom[ibody],ex_space[ibody],ey_space[ibody],
ez_space[ibody],inertia[ibody],omega[ibody]);
}
@ -1227,20 +1208,19 @@ void FixRigid::richardson(double *q, double *w,
{
// compute omega at 1/2 step from m at 1/2 step and q at 0
omega_from_angmom(m,ex,ey,ez,moments,w);
MathExtra::angmom_to_omega(m,ex,ey,ez,moments,w);
// full update from dq/dt = 1/2 w q
double wq[4];
vecquat(w,q,wq);
MathExtra::vecquat(w,q,wq);
double qfull[4];
qfull[0] = q[0] + dtq * wq[0];
qfull[1] = q[1] + dtq * wq[1];
qfull[2] = q[2] + dtq * wq[2];
qfull[3] = q[3] + dtq * wq[3];
qnormalize(qfull);
MathExtra::qnormalize(qfull);
// 1st half update from dq/dt = 1/2 w q
@ -1249,16 +1229,15 @@ void FixRigid::richardson(double *q, double *w,
qhalf[1] = q[1] + 0.5*dtq * wq[1];
qhalf[2] = q[2] + 0.5*dtq * wq[2];
qhalf[3] = q[3] + 0.5*dtq * wq[3];
qnormalize(qhalf);
MathExtra::qnormalize(qhalf);
// udpate ex,ey,ez from qhalf
// re-compute omega at 1/2 step from m at 1/2 step and q at 1/2 step
// recompute wq
exyz_from_q(qhalf,ex,ey,ez);
omega_from_angmom(m,ex,ey,ez,moments,w);
vecquat(w,qhalf,wq);
MathExtra::q_to_exyz(qhalf,ex,ey,ez);
MathExtra::angmom_to_omega(m,ex,ey,ez,moments,w);
MathExtra::vecquat(w,qhalf,wq);
// 2nd half update from dq/dt = 1/2 w q
@ -1266,8 +1245,7 @@ void FixRigid::richardson(double *q, double *w,
qhalf[1] += 0.5*dtq * wq[1];
qhalf[2] += 0.5*dtq * wq[2];
qhalf[3] += 0.5*dtq * wq[3];
qnormalize(qhalf);
MathExtra::qnormalize(qhalf);
// corrected Richardson update
@ -1275,9 +1253,9 @@ void FixRigid::richardson(double *q, double *w,
q[1] = 2.0*qhalf[1] - qfull[1];
q[2] = 2.0*qhalf[2] - qfull[2];
q[3] = 2.0*qhalf[3] - qfull[3];
MathExtra::qnormalize(q);
qnormalize(q);
exyz_from_q(q,ex,ey,ez);
MathExtra::q_to_exyz(q,ex,ey,ez);
}
/* ----------------------------------------------------------------------
@ -1511,304 +1489,6 @@ void FixRigid::deform(int flag)
domain->lamda2x(xcm[ibody],xcm[ibody]);
}
/* ----------------------------------------------------------------------
compute evalues and evectors of 3x3 real symmetric matrix
based on Jacobi rotations
adapted from Numerical Recipes jacobi() function
------------------------------------------------------------------------- */
int FixRigid::jacobi(double **matrix, double *evalues, double **evectors)
{
int i,j,k;
double tresh,theta,tau,t,sm,s,h,g,c,b[3],z[3];
for (i = 0; i < 3; i++) {
for (j = 0; j < 3; j++) evectors[i][j] = 0.0;
evectors[i][i] = 1.0;
}
for (i = 0; i < 3; i++) {
b[i] = evalues[i] = matrix[i][i];
z[i] = 0.0;
}
for (int iter = 1; iter <= MAXJACOBI; iter++) {
sm = 0.0;
for (i = 0; i < 2; i++)
for (j = i+1; j < 3; j++)
sm += fabs(matrix[i][j]);
if (sm == 0.0) return 0;
if (iter < 4) tresh = 0.2*sm/(3*3);
else tresh = 0.0;
for (i = 0; i < 2; i++) {
for (j = i+1; j < 3; j++) {
g = 100.0*fabs(matrix[i][j]);
if (iter > 4 && fabs(evalues[i])+g == fabs(evalues[i])
&& fabs(evalues[j])+g == fabs(evalues[j]))
matrix[i][j] = 0.0;
else if (fabs(matrix[i][j]) > tresh) {
h = evalues[j]-evalues[i];
if (fabs(h)+g == fabs(h)) t = (matrix[i][j])/h;
else {
theta = 0.5*h/(matrix[i][j]);
t = 1.0/(fabs(theta)+sqrt(1.0+theta*theta));
if (theta < 0.0) t = -t;
}
c = 1.0/sqrt(1.0+t*t);
s = t*c;
tau = s/(1.0+c);
h = t*matrix[i][j];
z[i] -= h;
z[j] += h;
evalues[i] -= h;
evalues[j] += h;
matrix[i][j] = 0.0;
for (k = 0; k < i; k++) rotate(matrix,k,i,k,j,s,tau);
for (k = i+1; k < j; k++) rotate(matrix,i,k,k,j,s,tau);
for (k = j+1; k < 3; k++) rotate(matrix,i,k,j,k,s,tau);
for (k = 0; k < 3; k++) rotate(evectors,k,i,k,j,s,tau);
}
}
}
for (i = 0; i < 3; i++) {
evalues[i] = b[i] += z[i];
z[i] = 0.0;
}
}
return 1;
}
/* ----------------------------------------------------------------------
perform a single Jacobi rotation
------------------------------------------------------------------------- */
void FixRigid::rotate(double **matrix, int i, int j, int k, int l,
double s, double tau)
{
double g = matrix[i][j];
double h = matrix[k][l];
matrix[i][j] = g-s*(h+g*tau);
matrix[k][l] = h+s*(g-h*tau);
}
/* ----------------------------------------------------------------------
create unit quaternion from space-frame ex,ey,ez
ex,ey,ez are columns of a rotation matrix
------------------------------------------------------------------------- */
void FixRigid::q_from_exyz(double *ex, double *ey, double *ez, double *q)
{
// squares of quaternion components
double q0sq = 0.25 * (ex[0] + ey[1] + ez[2] + 1.0);
double q1sq = q0sq - 0.5 * (ey[1] + ez[2]);
double q2sq = q0sq - 0.5 * (ex[0] + ez[2]);
double q3sq = q0sq - 0.5 * (ex[0] + ey[1]);
// some component must be greater than 1/4 since they sum to 1
// compute other components from it
if (q0sq >= 0.25) {
q[0] = sqrt(q0sq);
q[1] = (ey[2] - ez[1]) / (4.0*q[0]);
q[2] = (ez[0] - ex[2]) / (4.0*q[0]);
q[3] = (ex[1] - ey[0]) / (4.0*q[0]);
} else if (q1sq >= 0.25) {
q[1] = sqrt(q1sq);
q[0] = (ey[2] - ez[1]) / (4.0*q[1]);
q[2] = (ey[0] + ex[1]) / (4.0*q[1]);
q[3] = (ex[2] + ez[0]) / (4.0*q[1]);
} else if (q2sq >= 0.25) {
q[2] = sqrt(q2sq);
q[0] = (ez[0] - ex[2]) / (4.0*q[2]);
q[1] = (ey[0] + ex[1]) / (4.0*q[2]);
q[3] = (ez[1] + ey[2]) / (4.0*q[2]);
} else if (q3sq >= 0.25) {
q[3] = sqrt(q3sq);
q[0] = (ex[1] - ey[0]) / (4.0*q[3]);
q[1] = (ez[0] + ex[2]) / (4.0*q[3]);
q[2] = (ez[1] + ey[2]) / (4.0*q[3]);
} else
error->all("Quaternion creation numeric error");
qnormalize(q);
}
/* ----------------------------------------------------------------------
compute space-frame ex,ey,ez from current quaternion q
ex,ey,ez = space-frame coords of 1st,2nd,3rd principal axis
operation is ex = q' d q = Q d, where d is (1,0,0) = 1st axis in body frame
------------------------------------------------------------------------- */
void FixRigid::exyz_from_q(double *q, double *ex, double *ey, double *ez)
{
ex[0] = q[0]*q[0] + q[1]*q[1] - q[2]*q[2] - q[3]*q[3];
ex[1] = 2.0 * (q[1]*q[2] + q[0]*q[3]);
ex[2] = 2.0 * (q[1]*q[3] - q[0]*q[2]);
ey[0] = 2.0 * (q[1]*q[2] - q[0]*q[3]);
ey[1] = q[0]*q[0] - q[1]*q[1] + q[2]*q[2] - q[3]*q[3];
ey[2] = 2.0 * (q[2]*q[3] + q[0]*q[1]);
ez[0] = 2.0 * (q[1]*q[3] + q[0]*q[2]);
ez[1] = 2.0 * (q[2]*q[3] - q[0]*q[1]);
ez[2] = q[0]*q[0] - q[1]*q[1] - q[2]*q[2] + q[3]*q[3];
}
/* ----------------------------------------------------------------------
vector-quaternion multiply: c = a*b, where a = (0,a)
------------------------------------------------------------------------- */
void FixRigid::vecquat(double *a, double *b, double *c)
{
c[0] = -a[0]*b[1] - a[1]*b[2] - a[2]*b[3];
c[1] = b[0]*a[0] + a[1]*b[3] - a[2]*b[2];
c[2] = b[0]*a[1] + a[2]*b[1] - a[0]*b[3];
c[3] = b[0]*a[2] + a[0]*b[2] - a[1]*b[1];
}
/* ----------------------------------------------------------------------
quaternion-vector multiply: c = a*b, where b = (0,b)
------------------------------------------------------------------------- */
void FixRigid::quatvec(double *a, double *b, double *c)
{
c[0] = -a[1]*b[0] - a[2]*b[1] - a[3]*b[2];
c[1] = a[0]*b[0] + a[2]*b[2] - a[3]*b[1];
c[2] = a[0]*b[1] + a[3]*b[0] - a[1]*b[2];
c[3] = a[0]*b[2] + a[1]*b[1] - a[2]*b[0];
}
/* ----------------------------------------------------------------------
quaternion-quaternion multiply: c = a*b
------------------------------------------------------------------------- */
void FixRigid::quatquat(double *a, double *b, double *c)
{
c[0] = a[0]*b[0] - a[1]*b[1] - a[2]*b[2] - a[3]*b[3];
c[1] = a[0]*b[1] + b[0]*a[1] + a[2]*b[3] - a[3]*b[2];
c[2] = a[0]*b[2] + b[0]*a[2] + a[3]*b[1] - a[1]*b[3];
c[3] = a[0]*b[3] + b[0]*a[3] + a[1]*b[2] - a[2]*b[1];
}
/* ----------------------------------------------------------------------
quaternion multiply: c = inv(a)*b
a is a quaternion
b is a four component vector
c is a three component vector
------------------------------------------------------------------------- */
void FixRigid::invquatvec(double *a, double *b, double *c)
{
c[0] = -a[1]*b[0] + a[0]*b[1] + a[3]*b[2] - a[2]*b[3];
c[1] = -a[2]*b[0] - a[3]*b[1] + a[0]*b[2] + a[1]*b[3];
c[2] = -a[3]*b[0] + a[2]*b[1] - a[1]*b[2] + a[0]*b[3];
}
/* ----------------------------------------------------------------------
conjugate of a quaternion: qc = conjugate of q
assume q is of unit length
------------------------------------------------------------------------- */
void FixRigid::qconjugate(double *q, double *qc)
{
qc[0] = q[0];
qc[1] = -q[1];
qc[2] = -q[2];
qc[3] = -q[3];
}
/* ----------------------------------------------------------------------
normalize a quaternion
------------------------------------------------------------------------- */
void FixRigid::qnormalize(double *q)
{
double norm = 1.0 / sqrt(q[0]*q[0] + q[1]*q[1] + q[2]*q[2] + q[3]*q[3]);
q[0] *= norm;
q[1] *= norm;
q[2] *= norm;
q[3] *= norm;
}
/* ----------------------------------------------------------------------
matvec_rows: c = Ab, where rows of A are x, y, z
------------------------------------------------------------------------- */
void FixRigid::matvec_rows(double *x, double *y, double *z,
double *b, double *c)
{
c[0] = x[0]*b[0] + x[1]*b[1] + x[2]*b[2];
c[1] = y[0]*b[0] + y[1]*b[1] + y[2]*b[2];
c[2] = z[0]*b[0] + z[1]*b[1] + z[2]*b[2];
}
/* ----------------------------------------------------------------------
matvec_cols: c = Ab, where columns of A are x, y, z
------------------------------------------------------------------------- */
void FixRigid::matvec_cols(double *x, double *y, double *z,
double *b, double *c)
{
c[0] = x[0]*b[0] + y[0]*b[1] + z[0]*b[2];
c[1] = x[1]*b[0] + y[1]*b[1] + z[1]*b[2];
c[2] = x[2]*b[0] + y[2]*b[1] + z[2]*b[2];
}
/* ----------------------------------------------------------------------
compute omega from angular momentum, both in space frame
only know Idiag so need to do M = Iw in body frame
ex,ey,ez are column vectors of rotation matrix P
Mbody = P_transpose Mspace
wbody = Mbody / Idiag
wspace = P wbody
set wbody component to 0.0 if inertia component is 0.0
otherwise body can spin easily around that axis
------------------------------------------------------------------------- */
void FixRigid::omega_from_angmom(double *m, double *ex, double *ey, double *ez,
double *idiag, double *w)
{
double wbody[3];
if (idiag[0] == 0.0) wbody[0] = 0.0;
else wbody[0] = (m[0]*ex[0] + m[1]*ex[1] + m[2]*ex[2]) / idiag[0];
if (idiag[1] == 0.0) wbody[1] = 0.0;
else wbody[1] = (m[0]*ey[0] + m[1]*ey[1] + m[2]*ey[2]) / idiag[1];
if (idiag[2] == 0.0) wbody[2] = 0.0;
else wbody[2] = (m[0]*ez[0] + m[1]*ez[1] + m[2]*ez[2]) / idiag[2];
w[0] = wbody[0]*ex[0] + wbody[1]*ey[0] + wbody[2]*ez[0];
w[1] = wbody[0]*ex[1] + wbody[1]*ey[1] + wbody[2]*ez[1];
w[2] = wbody[0]*ex[2] + wbody[1]*ey[2] + wbody[2]*ez[2];
}
/* ----------------------------------------------------------------------
compute angular momentum from omega, both in space frame
only know Idiag so need to do M = Iw in body frame
ex,ey,ez are column vectors of rotation matrix P
wbody = P_transpose wspace
Mbody = Idiag wbody
Mspace = P Mbody
------------------------------------------------------------------------- */
void FixRigid::angmom_from_omega(double *w,
double *ex, double *ey, double *ez,
double *idiag, double *m)
{
double mbody[3];
mbody[0] = (w[0]*ex[0] + w[1]*ex[1] + w[2]*ex[2]) * idiag[0];
mbody[1] = (w[0]*ey[0] + w[1]*ey[1] + w[2]*ey[2]) * idiag[1];
mbody[2] = (w[0]*ez[0] + w[1]*ez[1] + w[2]*ez[2]) * idiag[2];
m[0] = mbody[0]*ex[0] + mbody[1]*ey[0] + mbody[2]*ez[0];
m[1] = mbody[0]*ex[1] + mbody[1]*ey[1] + mbody[2]*ez[1];
m[2] = mbody[0]*ex[2] + mbody[1]*ey[2] + mbody[2]*ez[2];
}
/* ----------------------------------------------------------------------
set space-frame coords and velocity of each atom in each rigid body
set orientation and rotation of extended particles
@ -1933,7 +1613,8 @@ void FixRigid::set_xv()
if (extended) {
double *shape,*quatatom;
AtomVecEllipsoid::Bonus *ebonus = avec_ellipsoid->bonus;
AtomVecEllipsoid::Bonus *ebonus;
if (avec_ellipsoid) ebonus = avec_ellipsoid->bonus;
double **omega_one = atom->omega;
double **angmom_one = atom->angmom;
double **mu = atom->mu;
@ -1950,8 +1631,8 @@ void FixRigid::set_xv()
}
if (eflags[i] & ORIENT_QUAT) {
quatatom = ebonus[ellipsoid[i]].quat;
quatquat(quat[ibody],qorient[i],quatatom);
qnormalize(quatatom);
MathExtra::quatquat(quat[ibody],qorient[i],quatatom);
MathExtra::qnormalize(quatatom);
}
if (eflags[i] & OMEGA) {
omega_one[i][0] = omega[ibody][0];
@ -1964,8 +1645,9 @@ void FixRigid::set_xv()
ione[0] = 0.2*rmass[i] * (shape[1]*shape[1] + shape[2]*shape[2]);
ione[1] = 0.2*rmass[i] * (shape[0]*shape[0] + shape[2]*shape[2]);
ione[2] = 0.2*rmass[i] * (shape[0]*shape[0] + shape[1]*shape[1]);
exyz_from_q(quatatom,exone,eyone,ezone);
angmom_from_omega(omega[ibody],exone,eyone,ezone,ione,angmom_one[i]);
MathExtra::q_to_exyz(quatatom,exone,eyone,ezone);
MathExtra::omega_to_angmom(omega[ibody],exone,eyone,ezone,ione,
angmom_one[i]);
}
}
}
@ -2075,7 +1757,8 @@ void FixRigid::set_v()
if (extended) {
double *shape,*quatatom;
AtomVecEllipsoid::Bonus *ebonus = avec_ellipsoid->bonus;
AtomVecEllipsoid::Bonus *ebonus;
if (avec_ellipsoid) ebonus = avec_ellipsoid->bonus;
double **omega_one = atom->omega;
double **angmom_one = atom->angmom;
int *ellipsoid = atom->ellipsoid;
@ -2095,8 +1778,9 @@ void FixRigid::set_v()
ione[0] = 0.2*rmass[i] * (shape[1]*shape[1] + shape[2]*shape[2]);
ione[1] = 0.2*rmass[i] * (shape[0]*shape[0] + shape[2]*shape[2]);
ione[2] = 0.2*rmass[i] * (shape[0]*shape[0] + shape[1]*shape[1]);
exyz_from_q(quatatom,exone,eyone,ezone);
angmom_from_omega(omega[ibody],exone,eyone,ezone,ione,angmom_one[i]);
MathExtra::q_to_exyz(quatatom,exone,eyone,ezone);
MathExtra::omega_to_angmom(omega[ibody],exone,eyone,ezone,ione,
angmom_one[i]);
}
}
}

18
src/fix_rigid.h
View File

@ -102,27 +102,9 @@ class FixRigid : public Fix {
int ORIENT_DIPOLE,ORIENT_QUAT;
int OMEGA,ANGMOM,TORQUE;
int jacobi(double **, double *, double **);
void rotate(double **, int, int, int, int, double, double);
void q_from_exyz(double *, double *, double *, double *);
void exyz_from_q(double *, double *, double *, double *);
void vecquat(double *, double *, double *);
void quatvec(double *, double *, double *);
void quatquat(double *, double *, double *);
void invquatvec(double *, double *, double *);
void qconjugate(double *, double *);
void qnormalize(double *);
void matvec_rows(double *, double *, double *, double *, double *);
void matvec_cols(double *, double *, double *, double *, double *);
void richardson(double *, double *, double *, double *,
double *, double *, double *);
void no_squish_rotate(int, double *, double *, double *, double);
void omega_from_angmom(double *, double *, double *,
double *, double *, double *);
void angmom_from_omega(double *, double *, double *,
double *, double *, double *);
void set_xv();
void set_v();
};

27
src/fix_rigid_nve.cpp
View File

@ -21,6 +21,7 @@
#include "stdio.h"
#include "string.h"
#include "fix_rigid_nve.h"
#include "math_extra.h"
#include "atom.h"
#include "domain.h"
#include "update.h"
@ -57,9 +58,9 @@ void FixRigidNVE::setup(int vflag)
double mbody[3];
for (int ibody = 0; ibody < nbody; ibody++) {
matvec_rows(ex_space[ibody],ey_space[ibody],ez_space[ibody],
MathExtra::matvec_rows(ex_space[ibody],ey_space[ibody],ez_space[ibody],
angmom[ibody],mbody);
quatvec(quat[ibody],mbody,conjqm[ibody]);
MathExtra::quatvec(quat[ibody],mbody,conjqm[ibody]);
conjqm[ibody][0] *= 2.0;
conjqm[ibody][1] *= 2.0;
conjqm[ibody][2] *= 2.0;
@ -98,9 +99,9 @@ void FixRigidNVE::initial_integrate(int vflag)
torque[ibody][1] *= tflag[ibody][1];
torque[ibody][2] *= tflag[ibody][2];
matvec_rows(ex_space[ibody],ey_space[ibody],ez_space[ibody],
MathExtra::matvec_rows(ex_space[ibody],ey_space[ibody],ez_space[ibody],
torque[ibody],tbody);
quatvec(quat[ibody],tbody,fquat);
MathExtra::quatvec(quat[ibody],tbody,fquat);
conjqm[ibody][0] += dtf2 * fquat[0];
conjqm[ibody][1] += dtf2 * fquat[1];
@ -119,17 +120,17 @@ void FixRigidNVE::initial_integrate(int vflag)
// transform p back to angmom
// update angular velocity
exyz_from_q(quat[ibody],ex_space[ibody],ey_space[ibody],
MathExtra::q_to_exyz(quat[ibody],ex_space[ibody],ey_space[ibody],
ez_space[ibody]);
invquatvec(quat[ibody],conjqm[ibody],mbody);
matvec_cols(ex_space[ibody],ey_space[ibody],ez_space[ibody],
MathExtra::invquatvec(quat[ibody],conjqm[ibody],mbody);
MathExtra::matvec_cols(ex_space[ibody],ey_space[ibody],ez_space[ibody],
mbody,angmom[ibody]);
angmom[ibody][0] *= 0.5;
angmom[ibody][1] *= 0.5;
angmom[ibody][2] *= 0.5;
omega_from_angmom(angmom[ibody],ex_space[ibody],ey_space[ibody],
MathExtra::angmom_to_omega(angmom[ibody],ex_space[ibody],ey_space[ibody],
ez_space[ibody],inertia[ibody],omega[ibody]);
}
@ -247,24 +248,24 @@ void FixRigidNVE::final_integrate()
torque[ibody][1] *= tflag[ibody][1];
torque[ibody][2] *= tflag[ibody][2];
matvec_rows(ex_space[ibody],ey_space[ibody],ez_space[ibody],
MathExtra::matvec_rows(ex_space[ibody],ey_space[ibody],ez_space[ibody],
torque[ibody],tbody);
quatvec(quat[ibody],tbody,fquat);
MathExtra::quatvec(quat[ibody],tbody,fquat);
conjqm[ibody][0] += dtf2 * fquat[0];
conjqm[ibody][1] += dtf2 * fquat[1];
conjqm[ibody][2] += dtf2 * fquat[2];
conjqm[ibody][3] += dtf2 * fquat[3];
invquatvec(quat[ibody],conjqm[ibody],mbody);
matvec_cols(ex_space[ibody],ey_space[ibody],ez_space[ibody],
MathExtra::invquatvec(quat[ibody],conjqm[ibody],mbody);
MathExtra::matvec_cols(ex_space[ibody],ey_space[ibody],ez_space[ibody],
mbody,angmom[ibody]);
angmom[ibody][0] *= 0.5;
angmom[ibody][1] *= 0.5;
angmom[ibody][2] *= 0.5;
omega_from_angmom(angmom[ibody],ex_space[ibody],ey_space[ibody],
MathExtra::angmom_to_omega(angmom[ibody],ex_space[ibody],ey_space[ibody],
ez_space[ibody],inertia[ibody],omega[ibody]);
}

30
src/fix_rigid_nvt.cpp
View File

@ -22,6 +22,7 @@
#include "stdlib.h"
#include "string.h"
#include "fix_rigid_nvt.h"
#include "math_extra.h"
#include "atom.h"
#include "domain.h"
#include "update.h"
@ -163,9 +164,9 @@ void FixRigidNVT::setup(int vflag)
double mbody[3];
for (int ibody = 0; ibody < nbody; ibody++) {
matvec_rows(ex_space[ibody],ey_space[ibody],ez_space[ibody],
MathExtra::matvec_rows(ex_space[ibody],ey_space[ibody],ez_space[ibody],
angmom[ibody],mbody);
quatvec(quat[ibody],mbody,conjqm[ibody]);
MathExtra::quatvec(quat[ibody],mbody,conjqm[ibody]);
conjqm[ibody][0] *= 2.0;
conjqm[ibody][1] *= 2.0;
conjqm[ibody][2] *= 2.0;
@ -224,9 +225,9 @@ void FixRigidNVT::initial_integrate(int vflag)
// step 1.3 - apply torque (body coords) to quaternion momentum
matvec_rows(ex_space[ibody],ey_space[ibody],ez_space[ibody],
MathExtra::matvec_rows(ex_space[ibody],ey_space[ibody],ez_space[ibody],
torque[ibody],tbody);
quatvec(quat[ibody],tbody,fquat);
MathExtra::quatvec(quat[ibody],tbody,fquat);
conjqm[ibody][0] += dtf2 * fquat[0];
conjqm[ibody][1] += dtf2 * fquat[1];
@ -249,17 +250,17 @@ void FixRigidNVT::initial_integrate(int vflag)
// transform p back to angmom
// update angular velocity
exyz_from_q(quat[ibody],ex_space[ibody],ey_space[ibody],
MathExtra::q_to_exyz(quat[ibody],ex_space[ibody],ey_space[ibody],
ez_space[ibody]);
invquatvec(quat[ibody],conjqm[ibody],mbody);
matvec_cols(ex_space[ibody],ey_space[ibody],ez_space[ibody],
MathExtra::invquatvec(quat[ibody],conjqm[ibody],mbody);
MathExtra::matvec_cols(ex_space[ibody],ey_space[ibody],ez_space[ibody],
mbody,angmom[ibody]);
angmom[ibody][0] *= 0.5;
angmom[ibody][1] *= 0.5;
angmom[ibody][2] *= 0.5;
omega_from_angmom(angmom[ibody],ex_space[ibody],ey_space[ibody],
MathExtra::angmom_to_omega(angmom[ibody],ex_space[ibody],ey_space[ibody],
ez_space[ibody],inertia[ibody],omega[ibody]);
akin_r += angmom[ibody][0]*omega[ibody][0] +
@ -397,12 +398,12 @@ void FixRigidNVT::final_integrate()
// convert torque to the body frame
matvec_rows(ex_space[ibody],ey_space[ibody],ez_space[ibody],
MathExtra::matvec_rows(ex_space[ibody],ey_space[ibody],ez_space[ibody],
torque[ibody],tbody);
// compute "force" for quaternion
quatvec(quat[ibody],tbody,fquat);
MathExtra::quatvec(quat[ibody],tbody,fquat);
// update the conjugate quaternion momentum (conjqm)
@ -411,10 +412,11 @@ void FixRigidNVT::final_integrate()
conjqm[ibody][2] = scale_r * conjqm[ibody][2] + dtf2 * fquat[2];
conjqm[ibody][3] = scale_r * conjqm[ibody][3] + dtf2 * fquat[3];
// compute angular momentum in the body frame then convert to the space-fixed frame
// compute angular momentum in the body frame
// then convert to the space-fixed frame
invquatvec(quat[ibody],conjqm[ibody],mbody);
matvec_cols(ex_space[ibody],ey_space[ibody],ez_space[ibody],
MathExtra::invquatvec(quat[ibody],conjqm[ibody],mbody);
MathExtra::matvec_cols(ex_space[ibody],ey_space[ibody],ez_space[ibody],
mbody,angmom[ibody]);
angmom[ibody][0] *= 0.5;
@ -423,7 +425,7 @@ void FixRigidNVT::final_integrate()
// compute new angular velocity
omega_from_angmom(angmom[ibody],ex_space[ibody],ey_space[ibody],
MathExtra::angmom_to_omega(angmom[ibody],ex_space[ibody],ey_space[ibody],
ez_space[ibody],inertia[ibody],omega[ibody]);
}

57
src/math_extra.h
View File

@ -64,8 +64,8 @@ namespace MathExtra {
void write3(const double mat[3][3]);
int mldivide3(const double mat[3][3], const double *vec, double *ans);
int jacobi(double **matrix, double *evalues, double **evectors);
void rotate(double **matrix, int i, int j, int k, int l,
int jacobi(double matrix[3][3], double *evalues, double evectors[3][3]);
void rotate(double matrix[3][3], int i, int j, int k, int l,
double s, double tau);
// shape matrix operations
@ -75,13 +75,8 @@ namespace MathExtra {
double *ans);
// quaternion operations
inline void normalize4(double *quat);
inline void axisangle_to_quat(const double *v, const double angle,
double *quat);
inline void multiply_quat_quat(const double *one, const double *two,
double *ans);
inline void multiply_vec_quat(const double *one, const double *two,
double *ans);
inline void vecquat(double *a, double *b, double *c);
inline void quatvec(double *a, double *b, double *c);
inline void quatquat(double *a, double *b, double *c);
@ -93,9 +88,9 @@ namespace MathExtra {
inline void matvec_cols(double *x, double *y, double *z,
double *b, double *c);
void omega_from_angmom(double *m, double *ex, double *ey, double *ez,
void angmom_to_omega(double *m, double *ex, double *ey, double *ez,
double *idiag, double *w);
void angmom_from_omega(double *w, double *ex, double *ey, double *ez,
void omega_to_angmom(double *w, double *ex, double *ey, double *ez,
double *idiag, double *m);
void exyz_to_q(double *ex, double *ey, double *ez, double *q);
void q_to_exyz(double *q, double *ex, double *ey, double *ez);
@ -424,20 +419,6 @@ void MathExtra::multiply_shape_shape(const double *one, const double *two,
ans[5] = one[0]*two[5] + one[5]*two[1];
}
/* ----------------------------------------------------------------------
normalize a quaternion
------------------------------------------------------------------------- */
void MathExtra::normalize4(double *quat)
{
double scale = 1.0/sqrt(quat[0]*quat[0]+quat[1]*quat[1] +
quat[2]*quat[2]+quat[3]*quat[3]);
quat[0] *= scale;
quat[1] *= scale;
quat[2] *= scale;
quat[3] *= scale;
}
/* ----------------------------------------------------------------------
compute quaternion from axis-angle rotation
v MUST be a unit vector
@ -454,35 +435,6 @@ void MathExtra::axisangle_to_quat(const double *v, const double angle,
quat[3] = v[2]*sina;
}
/* ----------------------------------------------------------------------
multiply 2 quaternions
effectively concatenates rotations
NOT a commutative operation
------------------------------------------------------------------------- */
void MathExtra::multiply_quat_quat(const double *one, const double *two,
double *ans)
{
ans[0] = one[0]*two[0]-one[1]*two[1]-one[2]*two[2]-one[3]*two[3];
ans[1] = one[0]*two[1]+one[1]*two[0]+one[2]*two[3]-one[3]*two[2];
ans[2] = one[0]*two[2]-one[1]*two[3]+one[2]*two[0]+one[3]*two[1];
ans[3] = one[0]*two[3]+one[1]*two[2]-one[2]*two[1]+one[3]*two[0];
}
/* ----------------------------------------------------------------------
multiply 3 vector times quaternion
3 vector one is treated as quaternion (0,one)
------------------------------------------------------------------------- */
void MathExtra::multiply_vec_quat(const double *one, const double *two,
double *ans)
{
ans[0] = -one[0]*two[1]-one[1]*two[2]-one[2]*two[3];
ans[1] = one[0]*two[0]+one[1]*two[3]-one[2]*two[2];
ans[2] = -one[0]*two[3]+one[1]*two[0]+one[2]*two[1];
ans[3] = one[0]*two[2]-one[1]*two[1]+one[2]*two[0];
}
/* ----------------------------------------------------------------------
vector-quaternion multiply: c = a*b, where a = (0,a)
------------------------------------------------------------------------- */
@ -509,6 +461,7 @@ void MathExtra::quatvec(double *a, double *b, double *c)
/* ----------------------------------------------------------------------
quaternion-quaternion multiply: c = a*b
NOT a commutative operation
------------------------------------------------------------------------- */
void MathExtra::quatquat(double *a, double *b, double *c)

2
src/set.cpp
View File

@ -431,7 +431,7 @@ void Set::set(int keyword)
quat[1] = xvalue * sintheta2;
quat[2] = yvalue * sintheta2;
quat[3] = zvalue * sintheta2;
MathExtra::normalize4(quat);
MathExtra::qnormalize(quat);
}
count++;
}