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

Browse Source
This commit is contained in:
sjplimp
2009-06-30 17:18:31 +00:00
parent 9248f42ba1
commit fc89638be0
3 changed files with 527 additions and 91 deletions
Download Patch File Download Diff File Expand all files Collapse all files

564
src/fix_rigid.cpp
View File

@ -16,7 +16,9 @@
#include "stdlib.h"
#include "string.h"
#include "fix_rigid.h"
#include "math_extra.h"
#include "atom.h"
#include "atom_vec.h"
#include "domain.h"
#include "update.h"
#include "respa.h"
@ -50,12 +52,16 @@ FixRigid::FixRigid(LAMMPS *lmp, int narg, char **arg) :
// perform initial allocation of atom-based arrays
// register with Atom class
extended = dorientflag = qorientflag = 0;
body = NULL;
displace = NULL;
eflags = NULL;
dorient = NULL;
qorient = NULL;
grow_arrays(atom->nmax);
atom->add_callback(0);
// parse command-line args
// set nbody and body[i] for each atom
@ -299,6 +305,17 @@ FixRigid::FixRigid(LAMMPS *lmp, int narg, char **arg) :
for (ibody = 0; ibody < nbody; ibody++)
image[ibody] = (512 << 20) | (512 << 10) | 512;
// bitmasks for properties of extended particles
INERTIA_SPHERE_RADIUS = 1;
INERTIA_SPHERE_SHAPE = 2;
INERTIA_ELLIPSOID = 4;
ORIENT_DIPOLE = 8;
ORIENT_QUAT = 16;
OMEGA = 32;
ANGMOM = 64;
TORQUE = 128;
// print statistics
int nsum = 0;
@ -322,7 +339,10 @@ FixRigid::~FixRigid()
memory->sfree(body);
memory->destroy_2d_double_array(displace);
memory->sfree(eflags);
memory->destroy_2d_double_array(dorient);
memory->destroy_2d_double_array(qorient);
// delete nbody-length arrays
memory->sfree(nrigid);
@ -364,7 +384,7 @@ int FixRigid::setmask()
void FixRigid::init()
{
int i,ibody;
int i,itype,ibody;
triclinic = domain->triclinic;
@ -396,14 +416,80 @@ void FixRigid::init()
if (strcmp(update->integrate_style,"respa") == 0)
step_respa = ((Respa *) update->integrate)->step;
// extended = 1 if any particle in a rigid body is finite size
double *radius = atom->radius;
double *rmass = atom->rmass;
double *mass = atom->mass;
double **shape = atom->shape;
double *dipole = atom->dipole;
int *type = atom->type;
int nlocal = atom->nlocal;
if (atom->radius_flag || atom->avec->shape_type) {
int flag = 0;
for (i = 0; i < nlocal; i++) {
if (body[i] < 0) continue;
if (radius && radius[i] > 0.0) flag = 1;
else if (shape && shape[type[i]][0] > 0.0) flag = 1;
}
MPI_Allreduce(&flag,&extended,1,MPI_INT,MPI_MAX,world);
}
// grow extended arrays and set extended flags for each particle
// vorientflag = 1 if any particles store dipole orientation
// qorientflag = 1 if any particles store quat orientation
if (extended) {
if (atom->mu_flag) dorientflag = 1;
if (atom->quat_flag) qorientflag = 1;
grow_arrays(atom->nmax);
for (i = 0; i < nlocal; i++) {
eflags[i] = 0;
if (body[i] < 0) continue;
// set INERTIA if radius or shape > 0.0
if (radius) {
if (radius[i] > 0.0) eflags[i] |= INERTIA_SPHERE_RADIUS;
} else if (shape) {
if (shape[type[i]][0] > 0.0) {
if (shape[type[i]][0] == shape[type[i]][1] &&
shape[type[i]][0] == shape[type[i]][2])
eflags[i] |= INERTIA_SPHERE_SHAPE;
else eflags[i] |= INERTIA_ELLIPSOID;
}
}
// other flags only set if particle is finite size
// set DIPOLE if atom->mu and atom->dipole exist and dipole[itype] > 0.0
// set QUAT if atom->quat exists (could be ellipsoid or sphere)
// set TORQUE if atom->torque exists
// set exactly one of OMEGA and ANGMOM so particle contributes once
// set OMEGA if either radius or rmass exists
// set ANGMOM if shape and mass exist
// set OMEGE if atom->angmom doesn't exist
if (eflags[i] == 0) continue;
if (atom->mu_flag && dipole && dipole[type[i]] > 0.0)
eflags[i] |= ORIENT_DIPOLE;
if (atom->quat_flag) eflags[i] |= ORIENT_QUAT;
if (atom->torque_flag) eflags[i] |= TORQUE;
if ((radius || rmass) && atom->omega_flag) eflags[i] |= OMEGA;
else if (shape && mass && atom->angmom_flag) eflags[i] |= ANGMOM;
else if (atom->omega_flag) eflags[i] != OMEGA;
else error->one("Could not set finite-size particle attribute "
"in fix rigid");
}
}
// compute masstotal & center-of-mass of each rigid body
double **x = atom->x;
double *rmass = atom->rmass;
double *mass = atom->mass;
int *image = atom->image;
int *type = atom->type;
int nlocal = atom->nlocal;
double xprd = domain->xprd;
double yprd = domain->yprd;
@ -426,7 +512,7 @@ void FixRigid::init()
zbox = (image[i] >> 20) - 512;
if (rmass) massone = rmass[i];
else massone = mass[type[i]];
if (triclinic == 0) {
xunwrap = x[i][0] + xbox*xprd;
yunwrap = x[i][1] + ybox*yprd;
@ -460,10 +546,11 @@ void FixRigid::init()
// compute 6 moments of inertia of each body
// dx,dy,dz = coords relative to center-of-mass
double dx,dy,dz,rad;
for (ibody = 0; ibody < nbody; ibody++)
for (i = 0; i < 6; i++) sum[ibody][i] = 0.0;
double dx,dy,dz;
for (i = 0; i < nlocal; i++) {
if (body[i] < 0) continue;
ibody = body[i];
@ -485,9 +572,10 @@ void FixRigid::init()
dx = xunwrap - xcm[ibody][0];
dy = yunwrap - xcm[ibody][1];
dz = zunwrap - xcm[ibody][2];
if (rmass) massone = rmass[i];
else massone = mass[type[i]];
sum[ibody][0] += massone * (dy*dy + dz*dz);
sum[ibody][1] += massone * (dx*dx + dz*dz);
sum[ibody][2] += massone * (dx*dx + dy*dy);
@ -495,6 +583,55 @@ void FixRigid::init()
sum[ibody][4] -= massone * dy*dz;
sum[ibody][5] -= massone * dx*dz;
}
// extended particles 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 *radius = atom->radius;
double **shape = atom->shape;
for (i = 0; i < nlocal; i++) {
if (body[i] < 0) continue;
ibody = body[i];
itype = type[i];
if (rmass) massone = rmass[i];
else massone = mass[itype];
if (eflags[i] & INERTIA_SPHERE_RADIUS) {
sum[ibody][0] += 0.4 * massone * radius[i]*radius[i];
sum[ibody][1] += 0.4 * massone * radius[i]*radius[i];
sum[ibody][2] += 0.4 * massone * radius[i]*radius[i];
}
if (eflags[i] & INERTIA_SPHERE_SHAPE) {
rad = shape[type[i]][0];
sum[ibody][0] += 0.4 * massone * rad*rad;
sum[ibody][1] += 0.4 * massone * rad*rad;
sum[ibody][2] += 0.4 * massone * rad*rad;
}
if (eflags[i] & INERTIA_ELLIPSOID) {
MathExtra::quat_to_mat(atom->quat[i],p);
MathExtra::quat_to_mat_trans(atom->quat[i],ptrans);
idiag[0] = 0.2*massone *
(shape[itype][1]*shape[itype][1]+shape[itype][2]*shape[itype][2]);
idiag[1] = 0.2*massone *
(shape[itype][0]*shape[itype][0]+shape[itype][2]*shape[itype][2]);
idiag[2] = 0.2*massone *
(shape[itype][0]*shape[itype][0]+shape[itype][1]*shape[itype][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];
}
}
}
MPI_Allreduce(sum[0],all[0],6*nbody,MPI_DOUBLE,MPI_SUM,world);
@ -569,36 +706,65 @@ void FixRigid::init()
// 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
double qc[4];
for (i = 0; i < nlocal; i++) {
if (body[i] < 0) displace[i][0] = displace[i][1] = displace[i][2] = 0.0;
else {
ibody = body[i];
if (body[i] < 0) {
displace[i][0] = displace[i][1] = displace[i][2] = 0.0;
continue;
}
xbox = (image[i] & 1023) - 512;
ybox = (image[i] >> 10 & 1023) - 512;
zbox = (image[i] >> 20) - 512;
ibody = body[i];
if (triclinic == 0) {
xunwrap = x[i][0] + xbox*xprd;
yunwrap = x[i][1] + ybox*yprd;
zunwrap = x[i][2] + zbox*zprd;
} else {
xunwrap = x[i][0] + xbox*xprd + ybox*xy + zbox*xz;
yunwrap = x[i][1] + ybox*yprd + zbox*yz;
zunwrap = x[i][2] + zbox*zprd;
}
xbox = (image[i] & 1023) - 512;
ybox = (image[i] >> 10 & 1023) - 512;
zbox = (image[i] >> 20) - 512;
if (triclinic == 0) {
xunwrap = x[i][0] + xbox*xprd;
yunwrap = x[i][1] + ybox*yprd;
zunwrap = x[i][2] + zbox*zprd;
} else {
xunwrap = x[i][0] + xbox*xprd + ybox*xy + zbox*xz;
yunwrap = x[i][1] + ybox*yprd + zbox*yz;
zunwrap = x[i][2] + zbox*zprd;
}
dx = xunwrap - xcm[ibody][0];
dy = yunwrap - xcm[ibody][1];
dz = zunwrap - xcm[ibody][2];
displace[i][0] = dx*ex_space[ibody][0] + dy*ex_space[ibody][1] +
dz*ex_space[ibody][2];
displace[i][1] = dx*ey_space[ibody][0] + dy*ey_space[ibody][1] +
dz*ey_space[ibody][2];
displace[i][2] = dx*ez_space[ibody][0] + dy*ez_space[ibody][1] +
dz*ez_space[ibody][2];
if (extended) {
double **mu = atom->mu;
double *dipole = atom->dipole;
int *type = atom->type;
dx = xunwrap - xcm[ibody][0];
dy = yunwrap - xcm[ibody][1];
dz = zunwrap - xcm[ibody][2];
displace[i][0] = dx*ex_space[ibody][0] + dy*ex_space[ibody][1] +
dz*ex_space[ibody][2];
displace[i][1] = dx*ey_space[ibody][0] + dy*ey_space[ibody][1] +
dz*ey_space[ibody][2];
displace[i][2] = dx*ez_space[ibody][0] + dy*ez_space[ibody][1] +
dz*ez_space[ibody][2];
if (eflags[i] & ORIENT_DIPOLE) {
dorient[i][0] = mu[i][0]*ex_space[ibody][0] +
mu[i][1]*ex_space[ibody][1] + mu[i][2]*ex_space[ibody][2];
dorient[i][1] = mu[i][0]*ey_space[ibody][0] +
mu[i][1]*ey_space[ibody][1] + mu[i][2]*ey_space[ibody][2];
dorient[i][2] = mu[i][0]*ez_space[ibody][0] +
mu[i][1]*ez_space[ibody][1] + mu[i][2]*ez_space[ibody][2];
MathExtra::snormalize3(dipole[type[i]],dorient[i],dorient[i]);
} else if (dorientflag)
dorient[i][0] = dorient[i][1] = dorient[i][2] = 0.0;
if (eflags[i] & ORIENT_QUAT) {
qconjugate(quat[ibody],qc);
quatquat(qc,atom->quat[i],qorient[i]);
qnormalize(qorient[i]);
} else if (qorientflag)
qorient[i][0] = qorient[i][1] = qorient[i][2] = qorient[i][3] = 0.0;
}
}
@ -606,6 +772,7 @@ void FixRigid::init()
// recompute moments of inertia around new axes
// 3 diagonal moments should equal principal moments
// 3 off-diagonal moments should be 0.0
// extended particles contribute extra terms to moments of inertia
for (ibody = 0; ibody < nbody; ibody++)
for (i = 0; i < 6; i++) sum[ibody][i] = 0.0;
@ -626,7 +793,54 @@ void FixRigid::init()
sum[ibody][4] -= massone * displace[i][1]*displace[i][2];
sum[ibody][5] -= massone * displace[i][0]*displace[i][2];
}
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 *radius = atom->radius;
double **shape = atom->shape;
for (i = 0; i < nlocal; i++) {
if (body[i] < 0) continue;
ibody = body[i];
itype = type[i];
if (rmass) massone = rmass[i];
else massone = mass[itype];
if (eflags[i] & INERTIA_SPHERE_RADIUS) {
sum[ibody][0] += 0.4 * massone * radius[i]*radius[i];
sum[ibody][1] += 0.4 * massone * radius[i]*radius[i];
sum[ibody][2] += 0.4 * massone * radius[i]*radius[i];
}
if (eflags[i] & INERTIA_SPHERE_SHAPE) {
rad = shape[type[i]][0];
sum[ibody][0] += 0.4 * massone * rad*rad;
sum[ibody][1] += 0.4 * massone * rad*rad;
sum[ibody][2] += 0.4 * massone * rad*rad;
}
if (eflags[i] & INERTIA_ELLIPSOID) {
MathExtra::quat_to_mat(qorient[i],p);
MathExtra::quat_to_mat_trans(qorient[i],ptrans);
idiag[0] = 0.2*massone *
(shape[itype][1]*shape[itype][1]+shape[itype][2]*shape[itype][2]);
idiag[1] = 0.2*massone *
(shape[itype][0]*shape[itype][0]+shape[itype][2]*shape[itype][2]);
idiag[2] = 0.2*massone *
(shape[itype][0]*shape[itype][0]+shape[itype][1]*shape[itype][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];
}
}
}
MPI_Allreduce(sum[0],all[0],6*nbody,MPI_DOUBLE,MPI_SUM,world);
for (ibody = 0; ibody < nbody; ibody++) {
@ -663,20 +877,20 @@ void FixRigid::init()
void FixRigid::setup(int vflag)
{
int i,n,ibody;
double massone,radone;
// vcm = velocity of center-of-mass of each rigid body
// fcm = force on center-of-mass of each rigid body
int *type = atom->type;
double **v = atom->v;
double **f = atom->f;
double *rmass = atom->rmass;
double *mass = atom->mass;
int *type = atom->type;
int nlocal = atom->nlocal;
for (ibody = 0; ibody < nbody; ibody++)
for (i = 0; i < 6; i++) sum[ibody][i] = 0.0;
double massone;
for (i = 0; i < nlocal; i++) {
if (body[i] < 0) continue;
@ -753,6 +967,42 @@ void FixRigid::setup(int vflag)
sum[ibody][4] += dz * f[i][0] - dx * f[i][2];
sum[ibody][5] += dx * f[i][1] - dy * f[i][0];
}
// extended particles add their rotation/torque to angmom/torque of body
if (extended) {
double **omega_one = atom->omega;
double **angmom_one = atom->angmom;
double **torque_one = atom->torque;
double *radius = atom->radius;
double **shape = atom->shape;
for (i = 0; i < nlocal; i++) {
if (body[i] < 0) continue;
ibody = body[i];
if (eflags[i] & OMEGA) {
if (rmass) massone = rmass[i];
else massone = mass[type[i]];
if (radius) radone = radius[i];
else radone = shape[type[i]][0];
sum[ibody][0] += 0.4 * massone * radone*radone * omega_one[i][0];
sum[ibody][1] += 0.4 * massone * radone*radone * omega_one[i][1];
sum[ibody][2] += 0.4 * massone * radone*radone * omega_one[i][2];
}
if (eflags[i] & ANGMOM) {
sum[ibody][0] += angmom_one[i][0];
sum[ibody][1] += angmom_one[i][1];
sum[ibody][2] += angmom_one[i][2];
}
if (eflags[i] & TORQUE) {
sum[ibody][3] += torque_one[i][0];
sum[ibody][4] += torque_one[i][1];
sum[ibody][5] += torque_one[i][2];
}
}
}
MPI_Allreduce(sum[0],all[0],6*nbody,MPI_DOUBLE,MPI_SUM,world);
@ -777,7 +1027,7 @@ void FixRigid::setup(int vflag)
ez_space[ibody],inertia[ibody],omega[ibody]);
set_v();
// guestimate virial as 2x the set_v contribution
// guesstimate virial as 2x the set_v contribution
if (vflag_global)
for (n = 0; n < 6; n++) virial[n] *= 2.0;
@ -829,7 +1079,7 @@ void FixRigid::initial_integrate(int vflag)
if (vflag) v_setup(vflag);
else evflag = 0;
// set coords and velocities of atoms in rigid bodies
// set coords/orient and velocity/rotation of atoms in rigid bodies
// from quarternion and omega
set_xv();
@ -856,7 +1106,7 @@ void FixRigid::richardson(double *q, double *w,
qfull[2] = q[2] + dtq * wq[2];
qfull[3] = q[3] + dtq * wq[3];
normalize(qfull);
qnormalize(qfull);
// 1st half update from dq/dt = 1/2 w q
@ -866,7 +1116,7 @@ void FixRigid::richardson(double *q, double *w,
qhalf[2] = q[2] + 0.5*dtq * wq[2];
qhalf[3] = q[3] + 0.5*dtq * wq[3];
normalize(qhalf);
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
@ -883,7 +1133,7 @@ void FixRigid::richardson(double *q, double *w,
qhalf[2] += 0.5*dtq * wq[2];
qhalf[3] += 0.5*dtq * wq[3];
normalize(qhalf);
qnormalize(qhalf);
// corrected Richardson update
@ -892,7 +1142,7 @@ void FixRigid::richardson(double *q, double *w,
q[2] = 2.0*qhalf[2] - qfull[2];
q[3] = 2.0*qhalf[3] - qfull[3];
normalize(q);
qnormalize(q);
exyz_from_q(q,ex,ey,ez);
}
@ -904,7 +1154,7 @@ void FixRigid::final_integrate()
double dtfm;
double xy,xz,yz;
// sum forces and torques on atoms in rigid body
// sum over atoms to get force and torque on rigid body
int *image = atom->image;
double **x = atom->x;
@ -955,7 +1205,24 @@ void FixRigid::final_integrate()
sum[ibody][4] += dz*f[i][0] - dx*f[i][2];
sum[ibody][5] += dx*f[i][1] - dy*f[i][0];
}
// extended particles add their torque to torque of body
if (extended) {
double **torque_one = atom->torque;
for (i = 0; i < nlocal; i++) {
if (body[i] < 0) continue;
ibody = body[i];
if (eflags[i] & TORQUE) {
sum[ibody][3] += torque_one[i][0];
sum[ibody][4] += torque_one[i][1];
sum[ibody][5] += torque_one[i][2];
}
}
}
MPI_Allreduce(sum[0],all[0],6*nbody,MPI_DOUBLE,MPI_SUM,world);
for (ibody = 0; ibody < nbody; ibody++) {
@ -980,7 +1247,7 @@ void FixRigid::final_integrate()
angmom[ibody][2] += dtf * torque[ibody][2] * tflag[ibody][2];
}
// set velocities from angmom & omega
// set velocity/rotation of atoms in rigid bodies
// virial is already setup from initial_integrate
for (ibody = 0; ibody < nbody; ibody++)
@ -1238,7 +1505,7 @@ void FixRigid::rotate(double **matrix, int i, int j, int k, int l,
}
/* ----------------------------------------------------------------------
create quaternion from space-frame ex,ey,ez
create unit quaternion from space-frame ex,ey,ez
ex,ey,ez are columns of a rotation matrix
------------------------------------------------------------------------- */
@ -1277,7 +1544,7 @@ void FixRigid::q_from_exyz(double *ex, double *ey, double *ez, double *q)
} else
error->all("Quaternion creation numeric error");
normalize(q);
qnormalize(q);
}
/* ----------------------------------------------------------------------
@ -1331,17 +1598,30 @@ void FixRigid::quatvec(double *a, double *b, double *c)
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[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];
}
/* ----------------------------------------------------------------------
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::normalize(double *q)
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;
@ -1355,23 +1635,23 @@ void FixRigid::normalize(double *q)
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 / Ibody
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 *inertia, double *w)
double *idiag, double *w)
{
double wbody[3];
if (inertia[0] == 0.0) wbody[0] = 0.0;
else wbody[0] = (m[0]*ex[0] + m[1]*ex[1] + m[2]*ex[2]) / inertia[0];
if (inertia[1] == 0.0) wbody[1] = 0.0;
else wbody[1] = (m[0]*ey[0] + m[1]*ey[1] + m[2]*ey[2]) / inertia[1];
if (inertia[2] == 0.0) wbody[2] = 0.0;
else wbody[2] = (m[0]*ez[0] + m[1]*ez[1] + m[2]*ez[2]) / inertia[2];
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];
@ -1383,18 +1663,19 @@ void FixRigid::omega_from_angmom(double *m, double *ex, double *ey, double *ez,
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 = Ibody wbody
Mbody = Idiag wbody
Mspace = P Mbody
------------------------------------------------------------------------- */
void FixRigid::angmom_from_omega(double *w, double *ex, double *ey, double *ez,
double *inertia, double *m)
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]) * inertia[0];
mbody[1] = (w[0]*ey[0] + w[1]*ey[1] + w[2]*ey[2]) * inertia[1];
mbody[2] = (w[0]*ez[0] + w[1]*ez[1] + w[2]*ez[2]) * inertia[2];
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];
@ -1403,17 +1684,18 @@ void FixRigid::angmom_from_omega(double *w, double *ex, double *ey, double *ez,
/* ----------------------------------------------------------------------
set space-frame coords and velocity of each atom in each rigid body
set orientation and rotation of extended particles
x = Q displace + Xcm, mapped back to periodic box
v = Vcm + (W cross (x - Xcm))
------------------------------------------------------------------------- */
void FixRigid::set_xv()
{
int ibody;
int ibody,itype;
int xbox,ybox,zbox;
double x0,x1,x2,v0,v1,v2,fc0,fc1,fc2,massone;
double xy,xz,yz;
double vr[6];
double ione[3],exone[3],eyone[3],ezone[3],vr[6],p[3][3];
int *image = atom->image;
double **x = atom->x;
@ -1504,7 +1786,6 @@ void FixRigid::set_xv()
if (evflag) {
if (rmass) massone = rmass[i];
else massone = mass[type[i]];
fc0 = massone*(v[i][0] - v0)/dtf - f[i][0];
fc1 = massone*(v[i][1] - v1)/dtf - f[i][1];
fc2 = massone*(v[i][2] - v2)/dtf - f[i][2];
@ -1519,25 +1800,66 @@ void FixRigid::set_xv()
v_tally(1,&i,1.0,vr);
}
}
// set orientation, omega, angmom of each extended particle
if (extended) {
double **omega_one = atom->omega;
double **angmom_one = atom->angmom;
double *dipole = atom->dipole;
double **shape = atom->shape;
for (int i = 0; i < nlocal; i++) {
if (body[i] < 0) continue;
ibody = body[i];
if (eflags[i] & ORIENT_DIPOLE) {
MathExtra::quat_to_mat(quat[ibody],p);
MathExtra::times_column3(p,dorient[i],atom->mu[i]);
MathExtra::snormalize3(dipole[type[i]],atom->mu[i],atom->mu[i]);
}
if (eflags[i] & ORIENT_QUAT) {
quatquat(quat[ibody],qorient[i],atom->quat[i]);
qnormalize(atom->quat[i]);
}
if (eflags[i] & OMEGA) {
omega_one[i][0] = omega[ibody][0];
omega_one[i][1] = omega[ibody][1];
omega_one[i][2] = omega[ibody][2];
}
if (eflags[i] & ANGMOM) {
itype = type[i];
ione[0] = 0.2*mass[itype] *
(shape[itype][1]*shape[itype][1] + shape[itype][2]*shape[itype][2]);
ione[1] = 0.2*mass[itype] *
(shape[itype][0]*shape[itype][0] + shape[itype][2]*shape[itype][2]);
ione[2] = 0.2*mass[itype] *
(shape[itype][0]*shape[itype][0] + shape[itype][1]*shape[itype][1]);
exyz_from_q(atom->quat[i],exone,eyone,ezone);
angmom_from_omega(omega[ibody],exone,eyone,ezone,ione,angmom_one[i]);
}
}
}
}
/* ----------------------------------------------------------------------
set space-frame velocity of each atom in a rigid body
set omega and angmom of extended particles
v = Vcm + (W cross (x - Xcm))
------------------------------------------------------------------------- */
void FixRigid::set_v()
{
int ibody;
int ibody,itype;
int xbox,ybox,zbox;
double dx,dy,dz;
double x0,x1,x2,v0,v1,v2,fc0,fc1,fc2,massone;
double xy,xz,yz;
double vr[6];
double ione[3],exone[3],eyone[3],ezone[3],vr[6];
double **f = atom->f;
double **v = atom->v;
double **x = atom->x;
double **v = atom->v;
double **f = atom->f;
double *rmass = atom->rmass;
double *mass = atom->mass;
int *type = atom->type;
@ -1590,7 +1912,6 @@ void FixRigid::set_v()
if (evflag) {
if (rmass) massone = rmass[i];
else massone = mass[type[i]];
fc0 = massone*(v[i][0] - v0)/dtf - f[i][0];
fc1 = massone*(v[i][1] - v1)/dtf - f[i][1];
fc2 = massone*(v[i][2] - v2)/dtf - f[i][2];
@ -1619,6 +1940,36 @@ void FixRigid::set_v()
v_tally(1,&i,1.0,vr);
}
}
// set omega, angmom of each extended particle
if (extended) {
double **omega_one = atom->omega;
double **angmom_one = atom->angmom;
double **shape = atom->shape;
for (int i = 0; i < nlocal; i++) {
if (body[i] < 0) continue;
ibody = body[i];
if (eflags[i] & OMEGA) {
omega_one[i][0] = omega[ibody][0];
omega_one[i][1] = omega[ibody][1];
omega_one[i][2] = omega[ibody][2];
}
if (eflags[i] & ANGMOM) {
itype = type[i];
ione[0] = 0.2*mass[itype] *
(shape[itype][1]*shape[itype][1] + shape[itype][2]*shape[itype][2]);
ione[1] = 0.2*mass[itype] *
(shape[itype][0]*shape[itype][0] + shape[itype][2]*shape[itype][2]);
ione[2] = 0.2*mass[itype] *
(shape[itype][0]*shape[itype][0] + shape[itype][1]*shape[itype][1]);
exyz_from_q(atom->quat[i],exone,eyone,ezone);
angmom_from_omega(omega[ibody],exone,eyone,ezone,ione,angmom_one[i]);
}
}
}
}
/* ----------------------------------------------------------------------
@ -1631,6 +1982,11 @@ double FixRigid::memory_usage()
double bytes = nmax * sizeof(int);
bytes += nmax*3 * sizeof(double);
bytes += maxvatom*6 * sizeof(double);
if (extended) {
bytes += nmax * sizeof(int);
if (dorientflag) bytes = nmax*3 * sizeof(double);
if (qorientflag) bytes = nmax*4 * sizeof(double);
}
return bytes;
}
@ -1642,6 +1998,14 @@ void FixRigid::grow_arrays(int nmax)
{
body = (int *) memory->srealloc(body,nmax*sizeof(int),"rigid:body");
displace = memory->grow_2d_double_array(displace,nmax,3,"rigid:displace");
if (extended) {
eflags = (int *)
memory->srealloc(eflags,nmax*sizeof(int),"rigid:eflags");
if (dorientflag)
dorient = memory->grow_2d_double_array(dorient,nmax,3,"rigid:dorient");
if (qorientflag)
qorient = memory->grow_2d_double_array(qorient,nmax,4,"rigid:qorient");
}
}
/* ----------------------------------------------------------------------
@ -1654,6 +2018,20 @@ void FixRigid::copy_arrays(int i, int j)
displace[j][0] = displace[i][0];
displace[j][1] = displace[i][1];
displace[j][2] = displace[i][2];
if (extended) {
eflags[j] = eflags[i];
if (dorientflag) {
dorient[j][0] = dorient[i][0];
dorient[j][1] = dorient[i][1];
dorient[j][2] = dorient[i][2];
}
if (qorientflag) {
qorient[j][0] = qorient[i][0];
qorient[j][1] = qorient[i][1];
qorient[j][2] = qorient[i][2];
qorient[j][3] = qorient[i][3];
}
}
}
/* ----------------------------------------------------------------------
@ -1666,7 +2044,22 @@ int FixRigid::pack_exchange(int i, double *buf)
buf[1] = displace[i][0];
buf[2] = displace[i][1];
buf[3] = displace[i][2];
return 4;
if (!extended) return 4;
int m = 4;
buf[m++] = eflags[i];
if (dorientflag) {
buf[m++] = dorient[i][0];
buf[m++] = dorient[i][1];
buf[m++] = dorient[i][2];
}
if (qorientflag) {
buf[m++] = qorient[i][0];
buf[m++] = qorient[i][1];
buf[m++] = qorient[i][2];
buf[m++] = qorient[i][3];
}
return m;
}
/* ----------------------------------------------------------------------
@ -1679,7 +2072,22 @@ int FixRigid::unpack_exchange(int nlocal, double *buf)
displace[nlocal][0] = buf[1];
displace[nlocal][1] = buf[2];
displace[nlocal][2] = buf[3];
return 4;
if (!extended) return 4;
int m = 4;
eflags[nlocal] = static_cast<int> (buf[m++]);
if (dorientflag) {
dorient[nlocal][0] = buf[m++];
dorient[nlocal][0] = buf[m++];
dorient[nlocal][0] = buf[m++];
}
if (qorientflag) {
qorient[nlocal][0] = buf[m++];
qorient[nlocal][0] = buf[m++];
qorient[nlocal][0] = buf[m++];
qorient[nlocal][0] = buf[m++];
}
return m;
}
/* ---------------------------------------------------------------------- */

21
src/fix_rigid.h
View File

@ -61,14 +61,28 @@ class FixRigid : public Fix {
double **torque; // torque on each rigid body in space coords
double **quat; // quaternion of each rigid body
int *image; // image flags of xcm of each rigid body
int *body; // which body each atom is part of (-1 if none)
double **displace; // displacement of each atom in body coords
double **fflag; // flag for on/off of center-of-mass force
double **tflag; // flag for on/off of center-of-mass torque
int *body; // which body each atom is part of (-1 if none)
double **displace; // displacement of each atom in body coords
double **sum,**all; // work vectors for each rigid body
int **remapflag; // PBC remap flags for each rigid body
int extended; // 1 if any particles have extended attributes
int dorientflag; // 1 if particles store dipole orientation
int qorientflag; // 1 if particles store quat orientation
int *eflags; // flags for extended particles
double **qorient; // rotation state of ext particle wrt rigid body
double **dorient; // orientation of dipole mu wrt rigid body
// bitmasks for eflags
int INERTIA_SPHERE_RADIUS,INERTIA_SPHERE_SHAPE,INERTIA_ELLIPSOID;
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 *);
@ -76,7 +90,8 @@ class FixRigid : public Fix {
void vecquat(double *, double *, double *);
void quatvec(double *, double *, double *);
void quatquat(double *, double *, double *);
void normalize(double *);
void qconjugate(double *, double *);
void qnormalize(double *);
void richardson(double *, double *, double *, double *,
double *, double *, double *);
void omega_from_angmom(double *, double *, double *,

33
src/math_extra.h
View File

@ -28,6 +28,7 @@ namespace MathExtra {
// 3 vector operations
inline void normalize3(const double *v, double *ans);
inline void snormalize3(const double, const double *v, double *ans);
inline double dot3(const double *v1, const double *v2);
inline void cross3(const double *v1, const double *v2, double *ans);
@ -86,10 +87,22 @@ namespace MathExtra {
void MathExtra::normalize3(const double *v, double *ans)
{
double den = sqrt(v[0]*v[0]+v[1]*v[1]+v[2]*v[2]);
ans[0] = v[0]/den;
ans[1] = v[1]/den;
ans[2] = v[2]/den;
double scale = 1.0/sqrt(v[0]*v[0]+v[1]*v[1]+v[2]*v[2]);
ans[0] = v[0]*scale;
ans[1] = v[1]*scale;
ans[2] = v[2]*scale;
}
/* ----------------------------------------------------------------------
scale a vector to length
------------------------------------------------------------------------- */
void MathExtra::snormalize3(const double length, const double *v, double *ans)
{
double scale = length/sqrt(v[0]*v[0]+v[1]*v[1]+v[2]*v[2]);
ans[0] = v[0]*scale;
ans[1] = v[1]*scale;
ans[2] = v[2]*scale;
}
/* ----------------------------------------------------------------------
@ -347,12 +360,12 @@ void MathExtra::write3(const double mat[3][3])
void MathExtra::normalize4(double *quat)
{
double den = sqrt(quat[0]*quat[0]+quat[1]*quat[1] +
quat[2]*quat[2]+quat[3]*quat[3]);
quat[0] /= den;
quat[1] /= den;
quat[2] /= den;
quat[3] /= den;
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;
}
/* ----------------------------------------------------------------------