diff --git a/src/ASPHERE/atom_vec_ellipsoid.cpp b/src/ASPHERE/atom_vec_ellipsoid.cpp
index 50820debe8993fcbdb2ef573ae22fe5ba947bbd5..e9a97127fdbc0da7f28a7b2bd720ef361f8a248c 100755
--- a/src/ASPHERE/atom_vec_ellipsoid.cpp
+++ b/src/ASPHERE/atom_vec_ellipsoid.cpp
@@ -30,6 +30,7 @@
using namespace LAMMPS_NS;
#define DELTA 10000
+#define DELTA_BONUS 10000
/* ---------------------------------------------------------------------- */
@@ -41,17 +42,27 @@ AtomVecEllipsoid::AtomVecEllipsoid(LAMMPS *lmp, int narg, char **arg) :
comm_x_only = comm_f_only = 0;
size_forward = 7;
size_reverse = 6;
- size_border = 13;
+ size_border = 14;
size_velocity = 6;
- size_data_atom = 13;
+ size_data_atom = 7;
size_data_vel = 7;
- xcol_data = 7;
+ size_data_bonus = 8;
+ xcol_data = 5;
atom->ellipsoid_flag = 1;
- atom->shape_flag = atom->rmass_flag = atom->quat_flag =
- atom->angmom_flag = atom->torque_flag = 1;
+ atom->rmass_flag = atom->angmom_flag = atom->torque_flag = 1;
PI = 4.0*atan(1.0);
+
+ nlocal_bonus = nghost_bonus = nmax_bonus = 0;
+ bonus = NULL;
+}
+
+/* ---------------------------------------------------------------------- */
+
+AtomVecEllipsoid::~AtomVecEllipsoid()
+{
+ memory->sfree(bonus);
}
/* ----------------------------------------------------------------------
@@ -76,12 +87,11 @@ void AtomVecEllipsoid::grow(int n)
v = memory->grow(atom->v,nmax,3,"atom:v");
f = memory->grow(atom->f,nmax,3,"atom:f");
- shape = memory->grow(atom->shape,nmax,3,"atom:shape");
rmass = memory->grow(atom->rmass,nmax,"atom:rmass");
- quat = memory->grow(atom->quat,nmax,4,"atom:quat");
angmom = memory->grow(atom->angmom,nmax,3,"atom:angmom");
torque = memory->grow(atom->torque,nmax,3,"atom:torque");
-
+ ellipsoid = memory->grow(atom->ellipsoid,nmax,"atom:ellipsoid");
+
if (atom->nextra_grow)
for (int iextra = 0; iextra < atom->nextra_grow; iextra++)
modify->fix[atom->extra_grow[iextra]]->grow_arrays(nmax);
@@ -96,13 +106,29 @@ void AtomVecEllipsoid::grow_reset()
tag = atom->tag; type = atom->type;
mask = atom->mask; image = atom->image;
x = atom->x; v = atom->v; f = atom->f;
- shape = atom->shape; rmass = atom->rmass;
- quat = atom->quat; angmom = atom->angmom; torque = atom->torque;
+ rmass = atom->rmass; angmom = atom->angmom; torque = atom->torque;
}
-/* ---------------------------------------------------------------------- */
+/* ----------------------------------------------------------------------
+ grow Bonus data for style-specific atom data
+------------------------------------------------------------------------- */
-void AtomVecEllipsoid::copy(int i, int j)
+void AtomVecEllipsoid::grow_bonus()
+{
+ nmax_bonus += DELTA_BONUS;
+ if (nmax_bonus < 0 || nmax_bonus > MAXSMALLINT)
+ error->one("Per-processor system is too big");
+
+ bonus = (Bonus *) memory->srealloc(bonus,nmax_bonus*sizeof(Bonus),
+ "atom:bonus");
+}
+
+/* ----------------------------------------------------------------------
+ copy atom I info to atom J
+ if delflag and atom J has style-specific data, then delete it
+------------------------------------------------------------------------- */
+
+void AtomVecEllipsoid::copy(int i, int j, int delflag)
{
tag[j] = tag[i];
type[j] = type[i];
@@ -115,23 +141,90 @@ void AtomVecEllipsoid::copy(int i, int j)
v[j][1] = v[i][1];
v[j][2] = v[i][2];
- shape[j][0] = shape[i][0];
- shape[j][1] = shape[i][1];
- shape[j][2] = shape[i][2];
rmass[j] = rmass[i];
- quat[j][0] = quat[i][0];
- quat[j][1] = quat[i][1];
- quat[j][2] = quat[i][2];
- quat[j][3] = quat[i][3];
angmom[j][0] = angmom[i][0];
angmom[j][1] = angmom[i][1];
angmom[j][2] = angmom[i][2];
+ if (delflag && ellipsoid[j] >= 0) {
+ copy_bonus(nlocal_bonus-1,ellipsoid[j]);
+ nlocal_bonus--;
+ }
+ ellipsoid[j] = ellipsoid[i];
+ if (ellipsoid[j] >= 0) bonus[ellipsoid[j]].ilocal = j;
+
if (atom->nextra_grow)
for (int iextra = 0; iextra < atom->nextra_grow; iextra++)
modify->fix[atom->extra_grow[iextra]]->copy_arrays(i,j);
}
+/* ----------------------------------------------------------------------
+ copy Bonus for I to J, effectively deleting the J entry
+ insure index pointers between per-atom and bonus data are updated
+------------------------------------------------------------------------- */
+
+void AtomVecEllipsoid::copy_bonus(int i, int j)
+{
+ double *ishape = bonus[i].shape;
+ double *iquat = bonus[i].quat;
+ double *jshape = bonus[j].shape;
+ double *jquat = bonus[j].quat;
+ jshape[0] = ishape[0];
+ jshape[1] = ishape[1];
+ jshape[2] = ishape[2];
+ jquat[0] = iquat[0];
+ jquat[1] = iquat[1];
+ jquat[2] = iquat[2];
+ jquat[3] = iquat[3];
+ int m = bonus[i].ilocal;
+ bonus[j].ilocal = m;
+ ellipsoid[m] = j;
+}
+
+/* ----------------------------------------------------------------------
+ set shape values in Bonus data for particle I
+ this may create or delete entry in Bonus data
+------------------------------------------------------------------------- */
+
+void AtomVecEllipsoid::set_bonus(int i,
+ double shapex, double shapey, double shapez)
+{
+ if (ellipsoid[i] < 0) {
+ if (shapex == 0.0 && shapey == 0.0 && shapez == 0.0) return;
+ if (nlocal_bonus == nmax_bonus) grow_bonus();
+ double *shape = bonus[nlocal_bonus].shape;
+ double *quat = bonus[nlocal_bonus].quat;
+ shape[0] = shapex;
+ shape[1] = shapey;
+ shape[2] = shapez;
+ quat[0] = 1.0;
+ quat[1] = 0.0;
+ quat[2] = 0.0;
+ quat[3] = 0.0;
+ bonus[nlocal_bonus].ilocal = i;
+ ellipsoid[i] = nlocal_bonus++;
+ } else if (shapex == 0.0 && shapey == 0.0 && shapez == 0.0) {
+ copy_bonus(nlocal_bonus-1,ellipsoid[i]);
+ nlocal_bonus--;
+ ellipsoid[i] = -1;
+ } else {
+ double *shape = bonus[ellipsoid[i]].shape;
+ shape[0] = shapex;
+ shape[1] = shapey;
+ shape[2] = shapez;
+ }
+}
+
+/* ----------------------------------------------------------------------
+ clear ghost info in Bonus data
+ called before ghosts are recommunicated in comm and irregular
+------------------------------------------------------------------------- */
+
+void AtomVecEllipsoid::clear_bonus()
+{
+ nghost_bonus = 0;
+}
+
/* ---------------------------------------------------------------------- */
int AtomVecEllipsoid::pack_comm(int n, int *list, double *buf,
@@ -139,6 +232,7 @@ int AtomVecEllipsoid::pack_comm(int n, int *list, double *buf,
{
int i,j,m;
double dx,dy,dz;
+ double *quat;
m = 0;
if (pbc_flag == 0) {
@@ -147,10 +241,13 @@ int AtomVecEllipsoid::pack_comm(int n, int *list, double *buf,
buf[m++] = x[j][0];
buf[m++] = x[j][1];
buf[m++] = x[j][2];
- buf[m++] = quat[j][0];
- buf[m++] = quat[j][1];
- buf[m++] = quat[j][2];
- buf[m++] = quat[j][3];
+ if (ellipsoid[j] >= 0) {
+ quat = bonus[ellipsoid[j]].quat;
+ buf[m++] = quat[0];
+ buf[m++] = quat[1];
+ buf[m++] = quat[2];
+ buf[m++] = quat[3];
+ }
}
} else {
if (domain->triclinic == 0) {
@@ -167,10 +264,13 @@ int AtomVecEllipsoid::pack_comm(int n, int *list, double *buf,
buf[m++] = x[j][0] + dx;
buf[m++] = x[j][1] + dy;
buf[m++] = x[j][2] + dz;
- buf[m++] = quat[j][0];
- buf[m++] = quat[j][1];
- buf[m++] = quat[j][2];
- buf[m++] = quat[j][3];
+ if (ellipsoid[j] >= 0) {
+ quat = bonus[ellipsoid[j]].quat;
+ buf[m++] = quat[0];
+ buf[m++] = quat[1];
+ buf[m++] = quat[2];
+ buf[m++] = quat[3];
+ }
}
}
return m;
@@ -183,6 +283,7 @@ int AtomVecEllipsoid::pack_comm_vel(int n, int *list, double *buf,
{
int i,j,m;
double dx,dy,dz,dvx,dvy,dvz;
+ double *quat;
m = 0;
if (pbc_flag == 0) {
@@ -191,10 +292,13 @@ int AtomVecEllipsoid::pack_comm_vel(int n, int *list, double *buf,
buf[m++] = x[j][0];
buf[m++] = x[j][1];
buf[m++] = x[j][2];
- buf[m++] = quat[j][0];
- buf[m++] = quat[j][1];
- buf[m++] = quat[j][2];
- buf[m++] = quat[j][3];
+ if (ellipsoid[j] >= 0) {
+ quat = bonus[ellipsoid[j]].quat;
+ buf[m++] = quat[0];
+ buf[m++] = quat[1];
+ buf[m++] = quat[2];
+ buf[m++] = quat[3];
+ }
buf[m++] = v[j][0];
buf[m++] = v[j][1];
buf[m++] = v[j][2];
@@ -218,10 +322,13 @@ int AtomVecEllipsoid::pack_comm_vel(int n, int *list, double *buf,
buf[m++] = x[j][0] + dx;
buf[m++] = x[j][1] + dy;
buf[m++] = x[j][2] + dz;
- buf[m++] = quat[j][0];
- buf[m++] = quat[j][1];
- buf[m++] = quat[j][2];
- buf[m++] = quat[j][3];
+ if (ellipsoid[j] >= 0) {
+ quat = bonus[ellipsoid[j]].quat;
+ buf[m++] = quat[0];
+ buf[m++] = quat[1];
+ buf[m++] = quat[2];
+ buf[m++] = quat[3];
+ }
buf[m++] = v[j][0];
buf[m++] = v[j][1];
buf[m++] = v[j][2];
@@ -238,10 +345,13 @@ int AtomVecEllipsoid::pack_comm_vel(int n, int *list, double *buf,
buf[m++] = x[j][0] + dx;
buf[m++] = x[j][1] + dy;
buf[m++] = x[j][2] + dz;
- buf[m++] = quat[j][0];
- buf[m++] = quat[j][1];
- buf[m++] = quat[j][2];
- buf[m++] = quat[j][3];
+ if (ellipsoid[j] >= 0) {
+ quat = bonus[ellipsoid[j]].quat;
+ buf[m++] = quat[0];
+ buf[m++] = quat[1];
+ buf[m++] = quat[2];
+ buf[m++] = quat[3];
+ }
if (mask[i] & deform_groupbit) {
buf[m++] = v[j][0] + dvx;
buf[m++] = v[j][1] + dvy;
@@ -265,14 +375,18 @@ int AtomVecEllipsoid::pack_comm_vel(int n, int *list, double *buf,
int AtomVecEllipsoid::pack_comm_hybrid(int n, int *list, double *buf)
{
int i,j,m;
+ double *quat;
m = 0;
for (i = 0; i < n; i++) {
j = list[i];
- buf[m++] = quat[i][0];
- buf[m++] = quat[i][1];
- buf[m++] = quat[i][2];
- buf[m++] = quat[i][3];
+ if (ellipsoid[j] >= 0) {
+ quat = bonus[ellipsoid[j]].quat;
+ buf[m++] = quat[0];
+ buf[m++] = quat[1];
+ buf[m++] = quat[2];
+ buf[m++] = quat[3];
+ }
}
return m;
}
@@ -282,6 +396,7 @@ int AtomVecEllipsoid::pack_comm_hybrid(int n, int *list, double *buf)
void AtomVecEllipsoid::unpack_comm(int n, int first, double *buf)
{
int i,m,last;
+ double *quat;
m = 0;
last = first + n;
@@ -289,10 +404,13 @@ void AtomVecEllipsoid::unpack_comm(int n, int first, double *buf)
x[i][0] = buf[m++];
x[i][1] = buf[m++];
x[i][2] = buf[m++];
- quat[i][0] = buf[m++];
- quat[i][1] = buf[m++];
- quat[i][2] = buf[m++];
- quat[i][3] = buf[m++];
+ if (ellipsoid[i] >= 0) {
+ quat = bonus[ellipsoid[i]].quat;
+ quat[0] = buf[m++];
+ quat[1] = buf[m++];
+ quat[2] = buf[m++];
+ quat[3] = buf[m++];
+ }
}
}
@@ -301,6 +419,7 @@ void AtomVecEllipsoid::unpack_comm(int n, int first, double *buf)
void AtomVecEllipsoid::unpack_comm_vel(int n, int first, double *buf)
{
int i,m,last;
+ double *quat;
m = 0;
last = first + n;
@@ -308,10 +427,13 @@ void AtomVecEllipsoid::unpack_comm_vel(int n, int first, double *buf)
x[i][0] = buf[m++];
x[i][1] = buf[m++];
x[i][2] = buf[m++];
- quat[i][0] = buf[m++];
- quat[i][1] = buf[m++];
- quat[i][2] = buf[m++];
- quat[i][3] = buf[m++];
+ if (ellipsoid[i] >= 0) {
+ quat = bonus[ellipsoid[i]].quat;
+ quat[0] = buf[m++];
+ quat[1] = buf[m++];
+ quat[2] = buf[m++];
+ quat[3] = buf[m++];
+ }
v[i][0] = buf[m++];
v[i][1] = buf[m++];
v[i][2] = buf[m++];
@@ -326,14 +448,18 @@ void AtomVecEllipsoid::unpack_comm_vel(int n, int first, double *buf)
int AtomVecEllipsoid::unpack_comm_hybrid(int n, int first, double *buf)
{
int i,m,last;
+ double *quat;
m = 0;
last = first + n;
for (i = first; i < last; i++) {
- quat[i][0] = buf[m++];
- quat[i][1] = buf[m++];
- quat[i][2] = buf[m++];
- quat[i][3] = buf[m++];
+ if (ellipsoid[i] >= 0) {
+ quat = bonus[ellipsoid[i]].quat;
+ quat[0] = buf[m++];
+ quat[1] = buf[m++];
+ quat[2] = buf[m++];
+ quat[3] = buf[m++];
+ }
}
return m;
}
@@ -414,6 +540,7 @@ int AtomVecEllipsoid::pack_border(int n, int *list, double *buf,
{
int i,j,m;
double dx,dy,dz;
+ double *shape,*quat;
m = 0;
if (pbc_flag == 0) {
@@ -425,13 +552,19 @@ int AtomVecEllipsoid::pack_border(int n, int *list, double *buf,
buf[m++] = tag[j];
buf[m++] = type[j];
buf[m++] = mask[j];
- buf[m++] = shape[j][0];
- buf[m++] = shape[j][1];
- buf[m++] = shape[j][2];
- buf[m++] = quat[j][0];
- buf[m++] = quat[j][1];
- buf[m++] = quat[j][2];
- buf[m++] = quat[j][3];
+ if (ellipsoid[j] < 0) buf[m++] = 0;
+ else {
+ buf[m++] = 1;
+ shape = bonus[ellipsoid[j]].shape;
+ quat = bonus[ellipsoid[j]].quat;
+ buf[m++] = shape[0];
+ buf[m++] = shape[1];
+ buf[m++] = shape[2];
+ buf[m++] = quat[0];
+ buf[m++] = quat[1];
+ buf[m++] = quat[2];
+ buf[m++] = quat[3];
+ }
}
} else {
if (domain->triclinic == 0) {
@@ -451,13 +584,19 @@ int AtomVecEllipsoid::pack_border(int n, int *list, double *buf,
buf[m++] = tag[j];
buf[m++] = type[j];
buf[m++] = mask[j];
- buf[m++] = shape[j][0];
- buf[m++] = shape[j][1];
- buf[m++] = shape[j][2];
- buf[m++] = quat[j][0];
- buf[m++] = quat[j][1];
- buf[m++] = quat[j][2];
- buf[m++] = quat[j][3];
+ if (ellipsoid[j] < 0) buf[m++] = 0;
+ else {
+ buf[m++] = 1;
+ shape = bonus[ellipsoid[j]].shape;
+ quat = bonus[ellipsoid[j]].quat;
+ buf[m++] = shape[0];
+ buf[m++] = shape[1];
+ buf[m++] = shape[2];
+ buf[m++] = quat[0];
+ buf[m++] = quat[1];
+ buf[m++] = quat[2];
+ buf[m++] = quat[3];
+ }
}
}
return m;
@@ -470,6 +609,7 @@ int AtomVecEllipsoid::pack_border_vel(int n, int *list, double *buf,
{
int i,j,m;
double dx,dy,dz,dvx,dvy,dvz;
+ double *shape,*quat;
m = 0;
if (pbc_flag == 0) {
@@ -481,13 +621,19 @@ int AtomVecEllipsoid::pack_border_vel(int n, int *list, double *buf,
buf[m++] = tag[j];
buf[m++] = type[j];
buf[m++] = mask[j];
- buf[m++] = shape[j][0];
- buf[m++] = shape[j][1];
- buf[m++] = shape[j][2];
- buf[m++] = quat[j][0];
- buf[m++] = quat[j][1];
- buf[m++] = quat[j][2];
- buf[m++] = quat[j][3];
+ if (ellipsoid[j] < 0) buf[m++] = 0;
+ else {
+ buf[m++] = 1;
+ shape = bonus[ellipsoid[j]].shape;
+ quat = bonus[ellipsoid[j]].quat;
+ buf[m++] = shape[0];
+ buf[m++] = shape[1];
+ buf[m++] = shape[2];
+ buf[m++] = quat[0];
+ buf[m++] = quat[1];
+ buf[m++] = quat[2];
+ buf[m++] = quat[3];
+ }
buf[m++] = v[j][0];
buf[m++] = v[j][1];
buf[m++] = v[j][2];
@@ -514,13 +660,18 @@ int AtomVecEllipsoid::pack_border_vel(int n, int *list, double *buf,
buf[m++] = tag[j];
buf[m++] = type[j];
buf[m++] = mask[j];
- buf[m++] = shape[j][0];
- buf[m++] = shape[j][1];
- buf[m++] = shape[j][2];
- buf[m++] = quat[j][0];
- buf[m++] = quat[j][1];
- buf[m++] = quat[j][2];
- buf[m++] = quat[j][3];
+ if (ellipsoid[j] < 0) buf[m++] = 0;
+ else {
+ buf[m++] = 1;
+ quat = bonus[ellipsoid[j]].quat;
+ buf[m++] = shape[0];
+ buf[m++] = shape[1];
+ buf[m++] = shape[2];
+ buf[m++] = quat[0];
+ buf[m++] = quat[1];
+ buf[m++] = quat[2];
+ buf[m++] = quat[3];
+ }
buf[m++] = v[j][0];
buf[m++] = v[j][1];
buf[m++] = v[j][2];
@@ -540,13 +691,19 @@ int AtomVecEllipsoid::pack_border_vel(int n, int *list, double *buf,
buf[m++] = tag[j];
buf[m++] = type[j];
buf[m++] = mask[j];
- buf[m++] = shape[j][0];
- buf[m++] = shape[j][1];
- buf[m++] = shape[j][2];
- buf[m++] = quat[j][0];
- buf[m++] = quat[j][1];
- buf[m++] = quat[j][2];
- buf[m++] = quat[j][3];
+ if (ellipsoid[j] < 0) buf[m++] = 0;
+ else {
+ buf[m++] = 1;
+ shape = bonus[ellipsoid[j]].shape;
+ quat = bonus[ellipsoid[j]].quat;
+ buf[m++] = shape[0];
+ buf[m++] = shape[1];
+ buf[m++] = shape[2];
+ buf[m++] = quat[0];
+ buf[m++] = quat[1];
+ buf[m++] = quat[2];
+ buf[m++] = quat[3];
+ }
if (mask[i] & deform_groupbit) {
buf[m++] = v[j][0] + dvx;
buf[m++] = v[j][1] + dvy;
@@ -570,17 +727,24 @@ int AtomVecEllipsoid::pack_border_vel(int n, int *list, double *buf,
int AtomVecEllipsoid::pack_border_hybrid(int n, int *list, double *buf)
{
int i,j,m;
+ double *shape,*quat;
m = 0;
for (i = 0; i < n; i++) {
j = list[i];
- buf[m++] = shape[j][0];
- buf[m++] = shape[j][1];
- buf[m++] = shape[j][2];
- buf[m++] = quat[j][0];
- buf[m++] = quat[j][1];
- buf[m++] = quat[j][2];
- buf[m++] = quat[j][3];
+ if (ellipsoid[j] < 0) buf[m++] = 0;
+ else {
+ buf[m++] = 1;
+ shape = bonus[ellipsoid[j]].shape;
+ quat = bonus[ellipsoid[j]].quat;
+ buf[m++] = shape[0];
+ buf[m++] = shape[1];
+ buf[m++] = shape[2];
+ buf[m++] = quat[0];
+ buf[m++] = quat[1];
+ buf[m++] = quat[2];
+ buf[m++] = quat[3];
+ }
}
return m;
}
@@ -589,7 +753,8 @@ int AtomVecEllipsoid::pack_border_hybrid(int n, int *list, double *buf)
void AtomVecEllipsoid::unpack_border(int n, int first, double *buf)
{
- int i,m,last;
+ int i,j,m,last;
+ double *shape,*quat;
m = 0;
last = first + n;
@@ -601,13 +766,24 @@ void AtomVecEllipsoid::unpack_border(int n, int first, double *buf)
tag[i] = static_cast<int> (buf[m++]);
type[i] = static_cast<int> (buf[m++]);
mask[i] = static_cast<int> (buf[m++]);
- shape[i][0] = buf[m++];
- shape[i][1] = buf[m++];
- shape[i][2] = buf[m++];
- quat[i][0] = buf[m++];
- quat[i][1] = buf[m++];
- quat[i][2] = buf[m++];
- quat[i][3] = buf[m++];
+ ellipsoid[i] = static_cast<int> (buf[m++]);
+ if (ellipsoid[i] < 0) ellipsoid[i] = 0;
+ else {
+ j = nlocal_bonus + nghost_bonus;
+ if (j == nmax_bonus) grow_bonus();
+ shape = bonus[j].shape;
+ quat = bonus[j].quat;
+ shape[0] = buf[m++];
+ shape[1] = buf[m++];
+ shape[2] = buf[m++];
+ quat[0] = buf[m++];
+ quat[1] = buf[m++];
+ quat[2] = buf[m++];
+ quat[3] = buf[m++];
+ bonus[j].ilocal = i;
+ ellipsoid[i] = j;
+ nghost_bonus++;
+ }
}
}
@@ -615,7 +791,8 @@ void AtomVecEllipsoid::unpack_border(int n, int first, double *buf)
void AtomVecEllipsoid::unpack_border_vel(int n, int first, double *buf)
{
- int i,m,last;
+ int i,j,m,last;
+ double *shape,*quat;
m = 0;
last = first + n;
@@ -627,13 +804,24 @@ void AtomVecEllipsoid::unpack_border_vel(int n, int first, double *buf)
tag[i] = static_cast<int> (buf[m++]);
type[i] = static_cast<int> (buf[m++]);
mask[i] = static_cast<int> (buf[m++]);
- shape[i][0] = buf[m++];
- shape[i][1] = buf[m++];
- shape[i][2] = buf[m++];
- quat[i][0] = buf[m++];
- quat[i][1] = buf[m++];
- quat[i][2] = buf[m++];
- quat[i][3] = buf[m++];
+ ellipsoid[i] = static_cast<int> (buf[m++]);
+ if (ellipsoid[i] < 0) ellipsoid[i] = 0;
+ else {
+ j = nlocal_bonus + nghost_bonus;
+ if (j == nmax_bonus) grow_bonus();
+ shape = bonus[j].shape;
+ quat = bonus[j].quat;
+ shape[0] = buf[m++];
+ shape[1] = buf[m++];
+ shape[2] = buf[m++];
+ quat[0] = buf[m++];
+ quat[1] = buf[m++];
+ quat[2] = buf[m++];
+ quat[3] = buf[m++];
+ bonus[j].ilocal = i;
+ ellipsoid[i] = j;
+ nghost_bonus++;
+ }
v[i][0] = buf[m++];
v[i][1] = buf[m++];
v[i][2] = buf[m++];
@@ -647,18 +835,30 @@ void AtomVecEllipsoid::unpack_border_vel(int n, int first, double *buf)
int AtomVecEllipsoid::unpack_border_hybrid(int n, int first, double *buf)
{
- int i,m,last;
+ int i,j,m,last;
+ double *shape,*quat;
m = 0;
last = first + n;
for (i = first; i < last; i++) {
- shape[i][0] = buf[m++];
- shape[i][1] = buf[m++];
- shape[i][2] = buf[m++];
- quat[i][0] = buf[m++];
- quat[i][1] = buf[m++];
- quat[i][2] = buf[m++];
- quat[i][3] = buf[m++];
+ ellipsoid[i] = static_cast<int> (buf[m++]);
+ if (ellipsoid[i] < 0) ellipsoid[i] = 0;
+ else {
+ j = nlocal_bonus + nghost_bonus;
+ if (j == nmax_bonus) grow_bonus();
+ shape = bonus[j].shape;
+ quat = bonus[j].quat;
+ shape[0] = buf[m++];
+ shape[1] = buf[m++];
+ shape[2] = buf[m++];
+ quat[0] = buf[m++];
+ quat[1] = buf[m++];
+ quat[2] = buf[m++];
+ quat[3] = buf[m++];
+ bonus[j].ilocal = i;
+ ellipsoid[i] = j;
+ nghost_bonus++;
+ }
}
return m;
}
@@ -682,17 +882,23 @@ int AtomVecEllipsoid::pack_exchange(int i, double *buf)
buf[m++] = mask[i];
buf[m++] = image[i];
- buf[m++] = shape[i][0];
- buf[m++] = shape[i][1];
- buf[m++] = shape[i][2];
buf[m++] = rmass[i];
- buf[m++] = quat[i][0];
- buf[m++] = quat[i][1];
- buf[m++] = quat[i][2];
- buf[m++] = quat[i][3];
buf[m++] = angmom[i][0];
buf[m++] = angmom[i][1];
buf[m++] = angmom[i][2];
+
+ if (ellipsoid[i] < 0) buf[m++] = 0;
+ else {
+ buf[m++] = 1;
+ int j = ellipsoid[i];
+ buf[m++] = bonus[j].shape[0];
+ buf[m++] = bonus[j].shape[1];
+ buf[m++] = bonus[j].shape[2];
+ buf[m++] = bonus[j].quat[0];
+ buf[m++] = bonus[j].quat[1];
+ buf[m++] = bonus[j].quat[2];
+ buf[m++] = bonus[j].quat[3];
+ }
if (atom->nextra_grow)
for (int iextra = 0; iextra < atom->nextra_grow; iextra++)
@@ -721,18 +927,27 @@ int AtomVecEllipsoid::unpack_exchange(double *buf)
mask[nlocal] = static_cast<int> (buf[m++]);
image[nlocal] = static_cast<int> (buf[m++]);
- shape[nlocal][0] = buf[m++];
- shape[nlocal][1] = buf[m++];
- shape[nlocal][2] = buf[m++];
rmass[nlocal] = buf[m++];
- quat[nlocal][0] = buf[m++];
- quat[nlocal][1] = buf[m++];
- quat[nlocal][2] = buf[m++];
- quat[nlocal][3] = buf[m++];
angmom[nlocal][0] = buf[m++];
angmom[nlocal][1] = buf[m++];
angmom[nlocal][2] = buf[m++];
-
+
+ ellipsoid[nlocal] = static_cast<int> (buf[m++]);
+ if (ellipsoid[nlocal]) {
+ if (nlocal_bonus == nmax_bonus) grow_bonus();
+ double *shape = bonus[nlocal_bonus].shape;
+ double *quat = bonus[nlocal_bonus].quat;
+ shape[0] = buf[m++];
+ shape[1] = buf[m++];
+ shape[2] = buf[m++];
+ quat[0] = buf[m++];
+ quat[1] = buf[m++];
+ quat[2] = buf[m++];
+ quat[3] = buf[m++];
+ bonus[nlocal_bonus].ilocal = nlocal;
+ ellipsoid[nlocal] = nlocal_bonus++;
+ }
+
if (atom->nextra_grow)
for (int iextra = 0; iextra < atom->nextra_grow; iextra++)
m += modify->fix[atom->extra_grow[iextra]]->
@@ -751,8 +966,11 @@ int AtomVecEllipsoid::size_restart()
{
int i;
+ int n = 0;
int nlocal = atom->nlocal;
- int n = 22 * nlocal;
+ for (i = 0; i < nlocal; i++)
+ if (ellipsoid[i] >= 0) n += 23;
+ else n += 16;
if (atom->nextra_restart)
for (int iextra = 0; iextra < atom->nextra_restart; iextra++)
@@ -763,7 +981,7 @@ int AtomVecEllipsoid::size_restart()
}
/* ----------------------------------------------------------------------
- pack atom I's data for restart file including extra quantities
+ pack atom I's data for restart file including Bonus data
xyz must be 1st 3 values, so that read_restart can test on them
molecular types may be negative, but write as positive
------------------------------------------------------------------------- */
@@ -782,18 +1000,24 @@ int AtomVecEllipsoid::pack_restart(int i, double *buf)
buf[m++] = v[i][1];
buf[m++] = v[i][2];
- buf[m++] = shape[i][0];
- buf[m++] = shape[i][1];
- buf[m++] = shape[i][2];
buf[m++] = rmass[i];
- buf[m++] = quat[i][0];
- buf[m++] = quat[i][1];
- buf[m++] = quat[i][2];
- buf[m++] = quat[i][3];
buf[m++] = angmom[i][0];
buf[m++] = angmom[i][1];
buf[m++] = angmom[i][2];
-
+
+ if (ellipsoid[i] < 0) buf[m++] = 0;
+ else {
+ buf[m++] = 1;
+ int j = ellipsoid[i];
+ buf[m++] = bonus[j].shape[0];
+ buf[m++] = bonus[j].shape[1];
+ buf[m++] = bonus[j].shape[2];
+ buf[m++] = bonus[j].quat[0];
+ buf[m++] = bonus[j].quat[1];
+ buf[m++] = bonus[j].quat[2];
+ buf[m++] = bonus[j].quat[3];
+ }
+
if (atom->nextra_restart)
for (int iextra = 0; iextra < atom->nextra_restart; iextra++)
m += modify->fix[atom->extra_restart[iextra]]->pack_restart(i,&buf[m]);
@@ -803,7 +1027,7 @@ int AtomVecEllipsoid::pack_restart(int i, double *buf)
}
/* ----------------------------------------------------------------------
- unpack data for one atom from restart file including extra quantities
+ unpack data for one atom from restart file including Bonus data
------------------------------------------------------------------------- */
int AtomVecEllipsoid::unpack_restart(double *buf)
@@ -827,17 +1051,26 @@ int AtomVecEllipsoid::unpack_restart(double *buf)
v[nlocal][1] = buf[m++];
v[nlocal][2] = buf[m++];
- shape[nlocal][0] = buf[m++];
- shape[nlocal][1] = buf[m++];
- shape[nlocal][2] = buf[m++];
rmass[nlocal] = buf[m++];
- quat[nlocal][0] = buf[m++];
- quat[nlocal][1] = buf[m++];
- quat[nlocal][2] = buf[m++];
- quat[nlocal][3] = buf[m++];
angmom[nlocal][0] = buf[m++];
angmom[nlocal][1] = buf[m++];
angmom[nlocal][2] = buf[m++];
+
+ ellipsoid[nlocal] = static_cast<int> (buf[m++]);
+ if (ellipsoid[nlocal]) {
+ if (nlocal_bonus == nmax_bonus) grow_bonus();
+ double *shape = bonus[nlocal_bonus].shape;
+ double *quat = bonus[nlocal_bonus].quat;
+ shape[0] = buf[m++];
+ shape[1] = buf[m++];
+ shape[2] = buf[m++];
+ quat[0] = buf[m++];
+ quat[1] = buf[m++];
+ quat[2] = buf[m++];
+ quat[3] = buf[m++];
+ bonus[nlocal_bonus].ilocal = nlocal;
+ ellipsoid[nlocal] = nlocal_bonus++;
+ }
double **extra = atom->extra;
if (atom->nextra_store) {
@@ -870,16 +1103,7 @@ void AtomVecEllipsoid::create_atom(int itype, double *coord)
v[nlocal][1] = 0.0;
v[nlocal][2] = 0.0;
- shape[nlocal][0] = 0.5;
- shape[nlocal][1] = 0.5;
- shape[nlocal][2] = 0.5;
- rmass[nlocal] = 4.0*PI/3.0 *
- shape[nlocal][0]*shape[nlocal][1]*shape[nlocal][2];
-
- quat[nlocal][0] = 1.0;
- quat[nlocal][1] = 0.0;
- quat[nlocal][2] = 0.0;
- quat[nlocal][3] = 0.0;
+ ellipsoid[nlocal] = -1;
angmom[nlocal][0] = 0.0;
angmom[nlocal][1] = 0.0;
angmom[nlocal][2] = 0.0;
@@ -905,38 +1129,19 @@ void AtomVecEllipsoid::data_atom(double *coord, int imagetmp, char **values)
if (type[nlocal] <= 0 || type[nlocal] > atom->ntypes)
error->one("Invalid atom type in Atoms section of data file");
- shape[nlocal][0] = 0.5 * atof(values[2]);
- shape[nlocal][1] = 0.5 * atof(values[3]);
- shape[nlocal][2] = 0.5 * atof(values[4]);
- if (shape[nlocal][0] < 0.0 || shape[nlocal][1] < 0.0 ||
- shape[nlocal][2] < 0.0)
- error->one("Invalid shape in Atoms section of data file");
- if (shape[nlocal][0] > 0.0 || shape[nlocal][1] > 0.0 ||
- shape[nlocal][2] > 0.0) {
- if (shape[nlocal][0] == 0.0 || shape[nlocal][1] == 0.0 ||
- shape[nlocal][2] == 0.0)
- error->one("Invalid shape in Atoms section of data file");
- }
+ ellipsoid[nlocal] = atoi(values[2]);
+ if (ellipsoid[nlocal] == 0) ellipsoid[nlocal] = -1;
+ else if (ellipsoid[nlocal] == 1) ellipsoid[nlocal] = 0;
+ else error->one("Invalid atom type in Atoms section of data file");
- double density = atof(values[5]);
- if (density <= 0.0)
+ rmass[nlocal] = atof(values[3]);
+ if (rmass[nlocal] <= 0.0)
error->one("Invalid density in Atoms section of data file");
- if (shape[nlocal][0] == 0.0) rmass[nlocal] = density;
- else
- rmass[nlocal] = 4.0*PI/3.0 *
- shape[nlocal][0]*shape[nlocal][1]*shape[nlocal][2] * density;
-
x[nlocal][0] = coord[0];
x[nlocal][1] = coord[1];
x[nlocal][2] = coord[2];
- quat[nlocal][0] = atof(values[9]);
- quat[nlocal][1] = atof(values[10]);
- quat[nlocal][2] = atof(values[11]);
- quat[nlocal][3] = atof(values[12]);
- MathExtra::normalize4(quat[nlocal]);
-
image[nlocal] = imagetmp;
mask[nlocal] = 1;
@@ -957,35 +1162,50 @@ void AtomVecEllipsoid::data_atom(double *coord, int imagetmp, char **values)
int AtomVecEllipsoid::data_atom_hybrid(int nlocal, char **values)
{
- shape[nlocal][0] = 0.5 * atof(values[0]);
- shape[nlocal][1] = 0.5 * atof(values[1]);
- shape[nlocal][2] = 0.5 * atof(values[2]);
- if (shape[nlocal][0] < 0.0 || shape[nlocal][1] < 0.0 ||
- shape[nlocal][2] < 0.0)
- error->one("Invalid shape in Atoms section of data file");
- if (shape[nlocal][0] > 0.0 || shape[nlocal][1] > 0.0 ||
- shape[nlocal][2] > 0.0) {
- if (shape[nlocal][0] == 0.0 || shape[nlocal][1] == 0.0 ||
- shape[nlocal][2] == 0.0)
- error->one("Invalid shape in Atoms section of data file");
- }
+ ellipsoid[nlocal] = atoi(values[0]);
+ if (ellipsoid[nlocal] == 0) ellipsoid[nlocal] = -1;
+ else if (ellipsoid[nlocal] == 1) ellipsoid[nlocal] = 0;
+ else error->one("Invalid atom type in Atoms section of data file");
- double density = atof(values[3]);
- if (density <= 0.0)
+ rmass[nlocal] = atof(values[1]);
+ if (rmass[nlocal] <= 0.0)
error->one("Invalid density in Atoms section of data file");
- if (shape[nlocal][0] == 0.0) rmass[nlocal] = density;
- else
- rmass[nlocal] = 4.0*PI/3.0 *
- shape[nlocal][0]*shape[nlocal][1]*shape[nlocal][2] * density;
+ return 2;
+}
+
+/* ----------------------------------------------------------------------
+ unpack one line from Ellipsoids section of data file
+------------------------------------------------------------------------- */
+
+void AtomVecEllipsoid::data_atom_bonus(int m, char **values)
+{
+ if (ellipsoid[m])
+ error->one("Assigning ellipsoid parameters to non-ellipsoid atom");
- quat[nlocal][0] = atof(values[4]);
- quat[nlocal][1] = atof(values[5]);
- quat[nlocal][2] = atof(values[6]);
- quat[nlocal][3] = atof(values[7]);
- MathExtra::normalize4(quat[nlocal]);
+ if (nlocal_bonus == nmax_bonus) grow_bonus();
+ double *shape = bonus[nlocal_bonus].shape;
+ double *quat = bonus[nlocal_bonus].quat;
- return 8;
+ shape[0] = 0.5 * atof(values[0]);
+ shape[1] = 0.5 * atof(values[1]);
+ shape[2] = 0.5 * atof(values[2]);
+ if (shape[0] <= 0.0 || shape[1] <= 0.0 || shape[2] <= 0.0)
+ error->one("Invalid shape in Ellipsoids section of data file");
+
+ quat[0] = atof(values[3]);
+ quat[1] = atof(values[4]);
+ quat[2] = atof(values[5]);
+ quat[3] = atof(values[6]);
+ MathExtra::normalize4(quat);
+
+ // reset ellipsoid mass
+ // previously stored density in rmass
+
+ rmass[m] *= 4.0*PI/3.0 * shape[0]*shape[1]*shape[2];
+
+ bonus[nlocal_bonus].ilocal = m;
+ ellipsoid[m] = nlocal_bonus++;
}
/* ----------------------------------------------------------------------
@@ -1030,11 +1250,11 @@ bigint AtomVecEllipsoid::memory_usage()
if (atom->memcheck("v")) bytes += memory->usage(v,nmax,3);
if (atom->memcheck("f")) bytes += memory->usage(f,nmax,3);
- if (atom->memcheck("shape")) bytes += memory->usage(shape,nmax,3);
if (atom->memcheck("rmass")) bytes += memory->usage(rmass,nmax);
- if (atom->memcheck("quat")) bytes += memory->usage(quat,nmax,4);
if (atom->memcheck("angmom")) bytes += memory->usage(angmom,nmax,3);
if (atom->memcheck("torque")) bytes += memory->usage(torque,nmax,3);
+
+ bytes += nmax_bonus*sizeof(Bonus);
return bytes;
}
diff --git a/src/ASPHERE/atom_vec_ellipsoid.h b/src/ASPHERE/atom_vec_ellipsoid.h
index b1a4a81341426f375df76dfbc68d670dfeee1a72..96998b1796d983e84025d92e8f4b7821ef296420 100755
--- a/src/ASPHERE/atom_vec_ellipsoid.h
+++ b/src/ASPHERE/atom_vec_ellipsoid.h
@@ -26,11 +26,18 @@ namespace LAMMPS_NS {
class AtomVecEllipsoid : public AtomVec {
public:
+ struct Bonus {
+ double shape[3];
+ double quat[4];
+ int ilocal;
+ };
+ struct Bonus *bonus;
+
AtomVecEllipsoid(class LAMMPS *, int, char **);
- virtual ~AtomVecEllipsoid() {}
+ virtual ~AtomVecEllipsoid();
void grow(int);
void grow_reset();
- void copy(int, int);
+ void copy(int, int, int);
int pack_comm(int, int *, double *, int, int *);
int pack_comm_vel(int, int *, double *, int, int *);
int pack_comm_hybrid(int, int *, double *);
@@ -59,12 +66,23 @@ class AtomVecEllipsoid : public AtomVec {
int data_vel_hybrid(int, char **);
bigint memory_usage();
+ // manipulate Bonus data structure for extra atom info
+
+ void grow_bonus();
+ void copy_bonus(int, int);
+ void set_bonus(int, double, double, double);
+ void clear_bonus();
+ void data_atom_bonus(int, char **);
+
private:
double PI;
int *tag,*type,*mask,*image;
double **x,**v,**f;
- double **shape,*density,*rmass;
- double **angmom,**torque,**quat;
+ double *density,*rmass;
+ double **angmom,**torque;
+ int *ellipsoid;
+
+ int nlocal_bonus,nghost_bonus,nmax_bonus;
};
}
diff --git a/src/ASPHERE/compute_erotate_asphere.cpp b/src/ASPHERE/compute_erotate_asphere.cpp
index b08a5553df1c7e9b8b7c4c222ecf7088efaf3983..0b63ba0ca7e149c63b813cc7bf9a9525a5171507 100644
--- a/src/ASPHERE/compute_erotate_asphere.cpp
+++ b/src/ASPHERE/compute_erotate_asphere.cpp
@@ -15,7 +15,7 @@
#include "compute_erotate_asphere.h"
#include "math_extra.h"
#include "atom.h"
-#include "atom_vec.h"
+#include "atom_vec_ellipsoid.h"
#include "update.h"
#include "force.h"
#include "memory.h"
@@ -36,7 +36,8 @@ ComputeERotateAsphere(LAMMPS *lmp, int narg, char **arg) :
// error check
- if (!atom->ellipsoid_flag)
+ avec = (AtomVecEllipsoid *) atom->style_match("ellipsoid");
+ if (!avec)
error->all("Compute erotate/asphere requires atom style ellipsoid");
}
@@ -47,13 +48,13 @@ void ComputeERotateAsphere::init()
// check that all particles are finite-size
// no point particles allowed, spherical is OK
- double **shape = atom->shape;
+ int *ellipsoid = atom->ellipsoid;
int *mask = atom->mask;
int nlocal = atom->nlocal;
for (int i = 0; i < nlocal; i++)
if (mask[i] & groupbit)
- if (shape[i][0] == 0.0)
+ if (ellipsoid[i] < 0)
error->one("Compute erotate/asphere requires extended particles");
pfactor = 0.5 * force->mvv2e;
@@ -65,9 +66,9 @@ double ComputeERotateAsphere::compute_scalar()
{
invoked_scalar = update->ntimestep;
- double **quat = atom->quat;
+ AtomVecEllipsoid::Bonus *bonus = avec->bonus;
+ int *ellipsoid = atom->ellipsoid;
double **angmom = atom->angmom;
- double **shape = atom->shape;
double *rmass = atom->rmass;
int *mask = atom->mask;
int nlocal = atom->nlocal;
@@ -75,6 +76,7 @@ double ComputeERotateAsphere::compute_scalar()
// sum rotational energy for each particle
// no point particles since divide by inertia
+ double *shape,*quat;
double wbody[3],inertia[3];
double rot[3][3];
double erotate = 0.0;
@@ -82,18 +84,18 @@ double ComputeERotateAsphere::compute_scalar()
for (int i = 0; i < nlocal; i++)
if (mask[i] & groupbit) {
+ shape = bonus[ellipsoid[i]].shape;
+ quat = bonus[ellipsoid[i]].quat;
+
// principal moments of inertia
- inertia[0] = rmass[i] *
- (shape[i][1]*shape[i][1]+shape[i][2]*shape[i][2]) / 5.0;
- inertia[1] = rmass[i] *
- (shape[i][0]*shape[i][0]+shape[i][2]*shape[i][2]) / 5.0;
- inertia[2] = rmass[i] *
- (shape[i][0]*shape[i][0]+shape[i][1]*shape[i][1]) / 5.0;
+ inertia[0] = rmass[i] * (shape[1]*shape[1]+shape[2]*shape[2]) / 5.0;
+ inertia[1] = rmass[i] * (shape[0]*shape[0]+shape[2]*shape[2]) / 5.0;
+ inertia[2] = rmass[i] * (shape[0]*shape[0]+shape[1]*shape[1]) / 5.0;
// wbody = angular velocity in body frame
- MathExtra::quat_to_mat(quat[i],rot);
+ MathExtra::quat_to_mat(quat,rot);
MathExtra::transpose_times_column3(rot,angmom[i],wbody);
wbody[0] /= inertia[0];
wbody[1] /= inertia[1];
diff --git a/src/ASPHERE/compute_erotate_asphere.h b/src/ASPHERE/compute_erotate_asphere.h
index f7f31066a887bd23ac2679cc8c2e9bcbea8762ed..0e3ae271d770005afbd45691395bda88890c4e4d 100644
--- a/src/ASPHERE/compute_erotate_asphere.h
+++ b/src/ASPHERE/compute_erotate_asphere.h
@@ -32,6 +32,7 @@ class ComputeERotateAsphere : public Compute {
private:
double pfactor;
+ class AtomVecEllipsoid *avec;
};
}
diff --git a/src/ASPHERE/compute_temp_asphere.cpp b/src/ASPHERE/compute_temp_asphere.cpp
index df29ad7b474668a60fced250f0620d2695c77043..ae73d31b15edc605aa1a103531b37b06cb926a10 100755
--- a/src/ASPHERE/compute_temp_asphere.cpp
+++ b/src/ASPHERE/compute_temp_asphere.cpp
@@ -19,7 +19,7 @@
#include "compute_temp_asphere.h"
#include "math_extra.h"
#include "atom.h"
-#include "atom_vec.h"
+#include "atom_vec_ellipsoid.h"
#include "update.h"
#include "force.h"
#include "domain.h"
@@ -56,9 +56,10 @@ ComputeTempAsphere::ComputeTempAsphere(LAMMPS *lmp, int narg, char **arg) :
vector = new double[6];
- // error checks
+ // error check
- if (!atom->ellipsoid_flag)
+ avec = (AtomVecEllipsoid *) atom->style_match("ellipsoid");
+ if (!avec)
error->all("Compute temp/asphere requires atom style ellipsoid");
}
@@ -77,13 +78,13 @@ void ComputeTempAsphere::init()
// check that all particles are finite-size
// no point particles allowed, spherical is OK
- double **shape = atom->shape;
+ int *ellipsoid = atom->ellipsoid;
int *mask = atom->mask;
int nlocal = atom->nlocal;
for (int i = 0; i < nlocal; i++)
if (mask[i] & groupbit)
- if (shape[i][0] == 0.0)
+ if (ellipsoid[i] < 0)
error->one("Compute temp/asphere requires extended particles");
if (tempbias) {
@@ -151,14 +152,15 @@ double ComputeTempAsphere::compute_scalar()
tbias->remove_bias_all();
}
+ AtomVecEllipsoid::Bonus *bonus = avec->bonus;
+ int *ellipsoid = atom->ellipsoid;
double **v = atom->v;
- double **quat = atom->quat;
double **angmom = atom->angmom;
- double **shape = atom->shape;
double *rmass = atom->rmass;
int *mask = atom->mask;
int nlocal = atom->nlocal;
+ double *shape,*quat;
double wbody[3],inertia[3];
double rot[3][3];
double t = 0.0;
@@ -169,20 +171,20 @@ double ComputeTempAsphere::compute_scalar()
for (int i = 0; i < nlocal; i++)
if (mask[i] & groupbit) {
+ shape = bonus[ellipsoid[i]].shape;
+ quat = bonus[ellipsoid[i]].quat;
+
t += (v[i][0]*v[i][0] + v[i][1]*v[i][1] + v[i][2]*v[i][2]) * rmass[i];
// principal moments of inertia
- inertia[0] = rmass[i] *
- (shape[i][1]*shape[i][1]+shape[i][2]*shape[i][2]) / 5.0;
- inertia[1] = rmass[i] *
- (shape[i][0]*shape[i][0]+shape[i][2]*shape[i][2]) / 5.0;
- inertia[2] = rmass[i] *
- (shape[i][0]*shape[i][0]+shape[i][1]*shape[i][1]) / 5.0;
+ inertia[0] = rmass[i] * (shape[1]*shape[1]+shape[2]*shape[2]) / 5.0;
+ inertia[1] = rmass[i] * (shape[0]*shape[0]+shape[2]*shape[2]) / 5.0;
+ inertia[2] = rmass[i] * (shape[0]*shape[0]+shape[1]*shape[1]) / 5.0;
// wbody = angular velocity in body frame
- MathExtra::quat_to_mat(quat[i],rot);
+ MathExtra::quat_to_mat(quat,rot);
MathExtra::transpose_times_column3(rot,angmom[i],wbody);
wbody[0] /= inertia[0];
wbody[1] /= inertia[1];
@@ -213,14 +215,15 @@ void ComputeTempAsphere::compute_vector()
tbias->remove_bias_all();
}
+ AtomVecEllipsoid::Bonus *bonus = avec->bonus;
+ int *ellipsoid = atom->ellipsoid;
double **v = atom->v;
- double **quat = atom->quat;
double **angmom = atom->angmom;
- double **shape = atom->shape;
double *rmass = atom->rmass;
int *mask = atom->mask;
int nlocal = atom->nlocal;
+ double *shape,*quat;
double wbody[3],inertia[3];
double rot[3][3];
double massone,t[6];
@@ -229,6 +232,9 @@ void ComputeTempAsphere::compute_vector()
for (i = 0; i < nlocal; i++)
if (mask[i] & groupbit) {
+ shape = bonus[ellipsoid[i]].shape;
+ quat = bonus[ellipsoid[i]].quat;
+
// translational kinetic energy
massone = rmass[i];
@@ -241,16 +247,13 @@ void ComputeTempAsphere::compute_vector()
// principal moments of inertia
- inertia[0] = rmass[i] *
- (shape[i][1]*shape[i][1]+shape[i][2]*shape[i][2]) / 5.0;
- inertia[1] = rmass[i] *
- (shape[i][0]*shape[i][0]+shape[i][2]*shape[i][2]) / 5.0;
- inertia[2] = rmass[i] *
- (shape[i][0]*shape[i][0]+shape[i][1]*shape[i][1]) / 5.0;
+ inertia[0] = rmass[i] * (shape[1]*shape[1]+shape[2]*shape[2]) / 5.0;
+ inertia[1] = rmass[i] * (shape[0]*shape[0]+shape[2]*shape[2]) / 5.0;
+ inertia[2] = rmass[i] * (shape[0]*shape[0]+shape[1]*shape[1]) / 5.0;
// wbody = angular velocity in body frame
- MathExtra::quat_to_mat(quat[i],rot);
+ MathExtra::quat_to_mat(quat,rot);
MathExtra::transpose_times_column3(rot,angmom[i],wbody);
wbody[0] /= inertia[0];
wbody[1] /= inertia[1];
diff --git a/src/ASPHERE/compute_temp_asphere.h b/src/ASPHERE/compute_temp_asphere.h
index 1504791e1bb79d8b8c97ee7c4adac14a89085090..19e29ebf1b8f3d497f8b08535f74d2065fad08b2 100755
--- a/src/ASPHERE/compute_temp_asphere.h
+++ b/src/ASPHERE/compute_temp_asphere.h
@@ -39,7 +39,8 @@ class ComputeTempAsphere : public Compute {
int fix_dof;
double tfactor;
char *id_bias;
- Compute *tbias; // ptr to additional bias compute
+ class Compute *tbias; // ptr to additional bias compute
+ class AtomVecEllipsoid *avec;
void dof_compute();
};
diff --git a/src/ASPHERE/fix_nh_asphere.cpp b/src/ASPHERE/fix_nh_asphere.cpp
index 9a587e6367958ed45aba3ebecb774949f0489a3f..19c18e402c9047690440f0738db13052028ccde3 100644
--- a/src/ASPHERE/fix_nh_asphere.cpp
+++ b/src/ASPHERE/fix_nh_asphere.cpp
@@ -21,7 +21,7 @@
#include "math_extra.h"
#include "fix_nh_asphere.h"
#include "atom.h"
-#include "atom_vec.h"
+#include "atom_vec_ellipsoid.h"
#include "group.h"
#include "memory.h"
#include "error.h"
@@ -33,8 +33,9 @@ using namespace LAMMPS_NS;
FixNHAsphere::FixNHAsphere(LAMMPS *lmp, int narg, char **arg) :
FixNH(lmp, narg, arg)
{
- if (!atom->ellipsoid_flag)
- error->all("Fix nvt/nph/npt asphere requires atom style ellipsoid");
+ avec = (AtomVecEllipsoid *) atom->style_match("ellipsoid");
+ if (!avec)
+ error->all("Compute nvt/nph/npt asphere requires atom style ellipsoid");
}
/* ---------------------------------------------------------------------- */
@@ -44,13 +45,13 @@ void FixNHAsphere::init()
// check that all particles are finite-size
// no point particles allowed, spherical is OK
- double **shape = atom->shape;
+ int *ellipsoid = atom->ellipsoid;
int *mask = atom->mask;
int nlocal = atom->nlocal;
for (int i = 0; i < nlocal; i++)
if (mask[i] & groupbit)
- if (shape[i][0] == 0.0)
+ if (ellipsoid[i] < 0)
error->one("Fix nvt/nph/npt asphere requires extended particles");
FixNH::init();
@@ -170,9 +171,9 @@ void FixNHAsphere::nve_x()
FixNH::nve_x();
- double **quat = atom->quat;
+ AtomVecEllipsoid::Bonus *bonus = avec->bonus;
+ int *ellipsoid = atom->ellipsoid;
double **angmom = atom->angmom;
- double **shape = atom->shape;
double *rmass = atom->rmass;
int *mask = atom->mask;
int nlocal = atom->nlocal;
@@ -186,18 +187,19 @@ void FixNHAsphere::nve_x()
// returns new normalized quaternion
// principal moments of inertia
+ double *shape,*quat;
double inertia[3];
for (int i = 0; i < nlocal; i++)
if (mask[i] & groupbit) {
- inertia[0] = rmass[i] *
- (shape[i][1]*shape[i][1]+shape[i][2]*shape[i][2]) / 5.0;
- inertia[1] = rmass[i] *
- (shape[i][0]*shape[i][0]+shape[i][2]*shape[i][2]) / 5.0;
- inertia[2] = rmass[i] *
- (shape[i][0]*shape[i][0]+shape[i][1]*shape[i][1]) / 5.0;
+ shape = bonus[ellipsoid[i]].shape;
+ quat = bonus[ellipsoid[i]].quat;
+
+ inertia[0] = rmass[i] * (shape[1]*shape[1]+shape[2]*shape[2]) / 5.0;
+ inertia[1] = rmass[i] * (shape[0]*shape[0]+shape[2]*shape[2]) / 5.0;
+ inertia[2] = rmass[i] * (shape[0]*shape[0]+shape[1]*shape[1]) / 5.0;
- richardson(quat[i],angmom[i],inertia);
+ richardson(quat,angmom[i],inertia);
}
}
diff --git a/src/ASPHERE/fix_nh_asphere.h b/src/ASPHERE/fix_nh_asphere.h
index 2d66264bbab652503df8baaafd66ef42cda6fea5..bfddb479759044eeaeb75f2b42bffd495137c2bf 100644
--- a/src/ASPHERE/fix_nh_asphere.h
+++ b/src/ASPHERE/fix_nh_asphere.h
@@ -27,6 +27,7 @@ class FixNHAsphere : public FixNH {
protected:
double dtq;
double **inertia;
+ class AtomVecEllipsoid *avec;
void richardson(double *, double *, double *);
void omega_from_mq(double *, double *, double *, double *);
diff --git a/src/ASPHERE/fix_nve_asphere.cpp b/src/ASPHERE/fix_nve_asphere.cpp
index fe5d71a4750f845253c77d83a8350a53a69d3c5d..b444cbf1140bc7a24e4498d7f0eba019bd6e2009 100755
--- a/src/ASPHERE/fix_nve_asphere.cpp
+++ b/src/ASPHERE/fix_nve_asphere.cpp
@@ -21,7 +21,7 @@
#include "fix_nve_asphere.h"
#include "math_extra.h"
#include "atom.h"
-#include "atom_vec.h"
+#include "atom_vec_ellipsoid.h"
#include "force.h"
#include "update.h"
#include "memory.h"
@@ -34,10 +34,9 @@ using namespace LAMMPS_NS;
FixNVEAsphere::FixNVEAsphere(LAMMPS *lmp, int narg, char **arg) :
FixNVE(lmp, narg, arg)
{
- // error checks
-
- if (!atom->ellipsoid_flag)
- error->all("Fix nve/asphere requires atom style ellipsoid");
+ avec = (AtomVecEllipsoid *) atom->style_match("ellipsoid");
+ if (!avec)
+ error->all("Compute nve/asphere requires atom style ellipsoid");
}
/* ---------------------------------------------------------------------- */
@@ -47,13 +46,13 @@ void FixNVEAsphere::init()
// check that all particles are finite-size
// no point particles allowed, spherical is OK
- double **shape = atom->shape;
+ int *ellipsoid = atom->ellipsoid;
int *mask = atom->mask;
int nlocal = atom->nlocal;
for (int i = 0; i < nlocal; i++)
if (mask[i] & groupbit)
- if (shape[i][0] == 0.0)
+ if (ellipsoid[i] < 0)
error->one("Fix nve/asphere requires extended particles");
FixNVE::init();
@@ -65,14 +64,15 @@ void FixNVEAsphere::initial_integrate(int vflag)
{
double dtfm;
double inertia[3];
+ double *shape,*quat;
+ AtomVecEllipsoid::Bonus *bonus = avec->bonus;
+ int *ellipsoid = atom->ellipsoid;
double **x = atom->x;
double **v = atom->v;
double **f = atom->f;
- double **quat = atom->quat;
double **angmom = atom->angmom;
double **torque = atom->torque;
- double **shape = atom->shape;
double *rmass = atom->rmass;
int *mask = atom->mask;
int nlocal = atom->nlocal;
@@ -102,14 +102,14 @@ void FixNVEAsphere::initial_integrate(int vflag)
// principal moments of inertia
- inertia[0] = rmass[i] *
- (shape[i][1]*shape[i][1]+shape[i][2]*shape[i][2]) / 5.0;
- inertia[1] = rmass[i] *
- (shape[i][0]*shape[i][0]+shape[i][2]*shape[i][2]) / 5.0;
- inertia[2] = rmass[i] *
- (shape[i][0]*shape[i][0]+shape[i][1]*shape[i][1]) / 5.0;
+ shape = bonus[ellipsoid[i]].shape;
+ quat = bonus[ellipsoid[i]].quat;
+
+ inertia[0] = rmass[i] * (shape[1]*shape[1]+shape[2]*shape[2]) / 5.0;
+ inertia[1] = rmass[i] * (shape[0]*shape[0]+shape[2]*shape[2]) / 5.0;
+ inertia[2] = rmass[i] * (shape[0]*shape[0]+shape[1]*shape[1]) / 5.0;
- richardson(quat[i],angmom[i],inertia);
+ richardson(quat,angmom[i],inertia);
}
}
diff --git a/src/ASPHERE/fix_nve_asphere.h b/src/ASPHERE/fix_nve_asphere.h
index b89237733110ea29255fc32795dd743b245e1c75..b535a47a14176c69fc09600091158e84aec19c8d 100755
--- a/src/ASPHERE/fix_nve_asphere.h
+++ b/src/ASPHERE/fix_nve_asphere.h
@@ -33,6 +33,7 @@ class FixNVEAsphere : public FixNVE {
private:
double dtq;
+ class AtomVecEllipsoid *avec;
void richardson(double *, double *, double *);
void omega_from_mq(double *, double *, double *, double *);
diff --git a/src/ASPHERE/pair_gayberne.cpp b/src/ASPHERE/pair_gayberne.cpp
index 0af5fe7b4af4f5f85efbd8015eff357753ce2721..629141d0c747d821d4c7c0794983bfde94b926c3 100755
--- a/src/ASPHERE/pair_gayberne.cpp
+++ b/src/ASPHERE/pair_gayberne.cpp
@@ -22,7 +22,7 @@
#include "pair_gayberne.h"
#include "math_extra.h"
#include "atom.h"
-#include "atom_vec.h"
+#include "atom_vec_ellipsoid.h"
#include "comm.h"
#include "force.h"
#include "neighbor.h"
@@ -42,6 +42,10 @@ 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");
+
single_enable = 0;
}
@@ -81,14 +85,16 @@ void PairGayBerne::compute(int eflag, int vflag)
double fforce[3],ttor[3],rtor[3],r12[3];
double a1[3][3],b1[3][3],g1[3][3],a2[3][3],b2[3][3],g2[3][3],temp[3][3];
int *ilist,*jlist,*numneigh,**firstneigh;
+ double *iquat,*jquat;
evdwl = 0.0;
if (eflag || vflag) ev_setup(eflag,vflag);
else evflag = vflag_fdotr = 0;
+ AtomVecEllipsoid::Bonus *bonus = avec->bonus;
+ int *ellipsoid = atom->ellipsoid;
double **x = atom->x;
double **f = atom->f;
- double **quat = atom->quat;
double **tor = atom->torque;
int *type = atom->type;
int nlocal = atom->nlocal;
@@ -107,7 +113,8 @@ void PairGayBerne::compute(int eflag, int vflag)
itype = type[i];
if (form[itype][itype] == ELLIPSE_ELLIPSE) {
- MathExtra::quat_to_mat_trans(quat[i],a1);
+ iquat = bonus[ellipsoid[i]].quat;
+ MathExtra::quat_to_mat_trans(iquat,a1);
MathExtra::diag_times3(well[itype],a1,temp);
MathExtra::transpose_times3(a1,temp,b1);
MathExtra::diag_times3(shape2[itype],a1,temp);
@@ -151,7 +158,8 @@ void PairGayBerne::compute(int eflag, int vflag)
break;
case SPHERE_ELLIPSE:
- MathExtra::quat_to_mat_trans(quat[j],a2);
+ jquat = bonus[ellipsoid[j]].quat;
+ MathExtra::quat_to_mat_trans(jquat,a2);
MathExtra::diag_times3(well[jtype],a2,temp);
MathExtra::transpose_times3(a2,temp,b2);
MathExtra::diag_times3(shape2[jtype],a2,temp);
@@ -166,7 +174,8 @@ void PairGayBerne::compute(int eflag, int vflag)
break;
default:
- MathExtra::quat_to_mat_trans(quat[j],a2);
+ jquat = bonus[ellipsoid[j]].quat;
+ MathExtra::quat_to_mat_trans(jquat,a2);
MathExtra::diag_times3(well[jtype],a2,temp);
MathExtra::transpose_times3(a2,temp,b2);
MathExtra::diag_times3(shape2[jtype],a2,temp);
@@ -332,9 +341,6 @@ void PairGayBerne::coeff(int narg, char **arg)
void PairGayBerne::init_style()
{
- if (!atom->ellipsoid_flag)
- error->all("Pair gayberne requires atom style ellipsoid");
-
int irequest = neighbor->request(this);
// per-type shape precalculations
diff --git a/src/ASPHERE/pair_gayberne.h b/src/ASPHERE/pair_gayberne.h
index 3d465be97e801f21b95dd37f3f1bc4d84eabf13b..ea699061cb6da4ba66b04b7922b736e76af62ae0 100755
--- a/src/ASPHERE/pair_gayberne.h
+++ b/src/ASPHERE/pair_gayberne.h
@@ -53,6 +53,7 @@ class PairGayBerne : public Pair {
double **lj1,**lj2,**lj3,**lj4;
double **offset;
int *setwell;
+ class AtomVecEllipsoid *avec;
void allocate();
double gayberne_analytic(const int i, const int j, double a1[3][3],
diff --git a/src/ASPHERE/pair_resquared.cpp b/src/ASPHERE/pair_resquared.cpp
index 6a9b1723dff865e1a4ca8434bf55534ef309e69f..a6b4ed9b371f2b6b18810c1fb153c9971a88884a 100755
--- a/src/ASPHERE/pair_resquared.cpp
+++ b/src/ASPHERE/pair_resquared.cpp
@@ -22,7 +22,7 @@
#include "pair_resquared.h"
#include "math_extra.h"
#include "atom.h"
-#include "atom_vec.h"
+#include "atom_vec_ellipsoid.h"
#include "comm.h"
#include "force.h"
#include "neighbor.h"
@@ -44,7 +44,12 @@ PairRESquared::PairRESquared(LAMMPS *lmp) : Pair(lmp),
b_alpha(45.0/56.0),
cr60(pow(60.0,1.0/3.0))
{
+ avec = (AtomVecEllipsoid *) atom->style_match("ellipsoid");
+ if (!avec)
+ error->all("Pair gayberne requires atom style ellipsoid");
+
single_enable = 0;
+
cr60 = pow(60.0,1.0/3.0);
b_alpha = 45.0/56.0;
solv_f_a = 3.0/(16.0*atan(1.0)*-36.0);
@@ -314,16 +319,12 @@ void PairRESquared::coeff(int narg, char **arg)
if (count == 0) error->all("Incorrect args for pair coefficients");
}
-
/* ----------------------------------------------------------------------
init specific to this pair style
------------------------------------------------------------------------- */
void PairRESquared::init_style()
{
- if (!atom->ellipsoid_flag)
- error->all("Pair resquared requires atom style ellipsoid");
-
int irequest = neighbor->request(this);
// per-type shape precalculations
@@ -497,10 +498,12 @@ void PairRESquared::read_restart_settings(FILE *fp)
Precompute per-particle temporaries for RE-squared calculation
------------------------------------------------------------------------- */
-void PairRESquared::precompute_i(const int i,RE2Vars &ws) {
+void PairRESquared::precompute_i(const int i,RE2Vars &ws)
+{
double aTs[3][3]; // A1'*S1^2
-
- MathExtra::quat_to_mat_trans(atom->quat[i],ws.A);
+ int *ellipsoid = atom->ellipsoid;
+ AtomVecEllipsoid::Bonus *bonus = avec->bonus;
+ MathExtra::quat_to_mat_trans(bonus[ellipsoid[i]].quat,ws.A);
MathExtra::transpose_times_diag3(ws.A,well[atom->type[i]],ws.aTe);
MathExtra::transpose_times_diag3(ws.A,shape2[atom->type[i]],aTs);
MathExtra::diag_times3(shape2[atom->type[i]],ws.A,ws.sa);
diff --git a/src/ASPHERE/pair_resquared.h b/src/ASPHERE/pair_resquared.h
index 85605a33f17e95738141c66ffebd283524ed016e..6aa2273c4520da6841be4de727e174e7ae786521 100755
--- a/src/ASPHERE/pair_resquared.h
+++ b/src/ASPHERE/pair_resquared.h
@@ -52,6 +52,7 @@ class PairRESquared : public Pair {
double **lj1,**lj2,**lj3,**lj4;
double **offset;
int *setwell;
+ class AtomVecEllipsoid *avec;
// per-particle temporaries for RE-squared calculation
diff --git a/src/DIPOLE/atom_vec_dipole.cpp b/src/DIPOLE/atom_vec_dipole.cpp
index 71c9d1407b0c020a707e7fb272413b05cdaabfef..a15d42933c03ed78527b08178d3aea46f00d21c1 100644
--- a/src/DIPOLE/atom_vec_dipole.cpp
+++ b/src/DIPOLE/atom_vec_dipole.cpp
@@ -89,9 +89,11 @@ void AtomVecDipole::grow_reset()
q = atom->q; mu = atom->mu;
}
-/* ---------------------------------------------------------------------- */
+/* ----------------------------------------------------------------------
+ copy atom I info to atom J
+------------------------------------------------------------------------- */
-void AtomVecDipole::copy(int i, int j)
+void AtomVecDipole::copy(int i, int j, int delflag)
{
tag[j] = tag[i];
type[j] = type[i];
diff --git a/src/DIPOLE/atom_vec_dipole.h b/src/DIPOLE/atom_vec_dipole.h
index 8e8c77bdc88d17a9255782a0f70d718734d2144f..c41a9e50808f670c77006600eeb5ab1c54d979df 100644
--- a/src/DIPOLE/atom_vec_dipole.h
+++ b/src/DIPOLE/atom_vec_dipole.h
@@ -29,7 +29,7 @@ class AtomVecDipole : public AtomVec {
AtomVecDipole(class LAMMPS *, int, char **);
void grow(int);
void grow_reset();
- void copy(int, int);
+ void copy(int, int, int);
int pack_comm(int, int *, double *, int, int *);
int pack_comm_vel(int, int *, double *, int, int *);
int pack_comm_hybrid(int, int *, double *);
diff --git a/src/MOLECULE/atom_vec_angle.cpp b/src/MOLECULE/atom_vec_angle.cpp
index 26e3968bcf022d5419f01d68537f6ee6731e491f..facbf8f10bf7871a268c79d36d89dd272b1f0870 100644
--- a/src/MOLECULE/atom_vec_angle.cpp
+++ b/src/MOLECULE/atom_vec_angle.cpp
@@ -114,9 +114,11 @@ void AtomVecAngle::grow_reset()
angle_atom3 = atom->angle_atom3;
}
-/* ---------------------------------------------------------------------- */
+/* ----------------------------------------------------------------------
+ copy atom I info to atom J
+------------------------------------------------------------------------- */
-void AtomVecAngle::copy(int i, int j)
+void AtomVecAngle::copy(int i, int j, int delflag)
{
int k;
diff --git a/src/MOLECULE/atom_vec_angle.h b/src/MOLECULE/atom_vec_angle.h
index 009a38d8ebc06eaec865f3763582f9c7a1cc5f65..daf18f2198aa22256937184ce1398edbff3911d1 100644
--- a/src/MOLECULE/atom_vec_angle.h
+++ b/src/MOLECULE/atom_vec_angle.h
@@ -29,7 +29,7 @@ class AtomVecAngle : public AtomVec {
AtomVecAngle(class LAMMPS *, int, char **);
void grow(int);
void grow_reset();
- void copy(int, int);
+ void copy(int, int, int);
int pack_comm(int, int *, double *, int, int *);
int pack_comm_vel(int, int *, double *, int, int *);
void unpack_comm(int, int, double *);
diff --git a/src/MOLECULE/atom_vec_bond.cpp b/src/MOLECULE/atom_vec_bond.cpp
index 81660e1e2f65d77a9634f921e9b2aa2a177a2d3e..ba433ea52078c799de3ae2bf628ff80e4aa31adc 100644
--- a/src/MOLECULE/atom_vec_bond.cpp
+++ b/src/MOLECULE/atom_vec_bond.cpp
@@ -102,9 +102,11 @@ void AtomVecBond::grow_reset()
bond_atom = atom->bond_atom;
}
-/* ---------------------------------------------------------------------- */
+/* ----------------------------------------------------------------------
+ copy atom I info to atom J
+------------------------------------------------------------------------- */
-void AtomVecBond::copy(int i, int j)
+void AtomVecBond::copy(int i, int j, int delflag)
{
int k;
diff --git a/src/MOLECULE/atom_vec_bond.h b/src/MOLECULE/atom_vec_bond.h
index 3274e3a0cfbad55d4be649c902e75627121ca91f..6f7c793033934be197ae1a27d8fadbd12c542c52 100644
--- a/src/MOLECULE/atom_vec_bond.h
+++ b/src/MOLECULE/atom_vec_bond.h
@@ -29,7 +29,7 @@ class AtomVecBond : public AtomVec {
AtomVecBond(class LAMMPS *, int, char **);
void grow(int);
void grow_reset();
- void copy(int, int);
+ void copy(int, int, int);
int pack_comm(int, int *, double *, int, int *);
int pack_comm_vel(int, int *, double *, int, int *);
void unpack_comm(int, int, double *);
diff --git a/src/MOLECULE/atom_vec_full.cpp b/src/MOLECULE/atom_vec_full.cpp
index 428c8957cacd2ca807661e6ef674ca2b109c4b24..78e99abc119a4fd05ca5e0490f8d51c7d479b596 100644
--- a/src/MOLECULE/atom_vec_full.cpp
+++ b/src/MOLECULE/atom_vec_full.cpp
@@ -155,9 +155,11 @@ void AtomVecFull::grow_reset()
improper_atom3 = atom->improper_atom3; improper_atom4 = atom->improper_atom4;
}
-/* ---------------------------------------------------------------------- */
+/* ----------------------------------------------------------------------
+ copy atom I info to atom J
+------------------------------------------------------------------------- */
-void AtomVecFull::copy(int i, int j)
+void AtomVecFull::copy(int i, int j, int delflag)
{
int k;
diff --git a/src/MOLECULE/atom_vec_full.h b/src/MOLECULE/atom_vec_full.h
index c2c23881d7ddb593d11b5c327599000d895242a1..408cb518ff52888cc30094257bfd0cf7adfdfbc7 100644
--- a/src/MOLECULE/atom_vec_full.h
+++ b/src/MOLECULE/atom_vec_full.h
@@ -29,7 +29,7 @@ class AtomVecFull : public AtomVec {
AtomVecFull(class LAMMPS *, int, char **);
void grow(int);
void grow_reset();
- void copy(int, int);
+ void copy(int, int, int);
int pack_comm(int, int *, double *, int, int *);
int pack_comm_vel(int, int *, double *, int, int *);
void unpack_comm(int, int, double *);
diff --git a/src/MOLECULE/atom_vec_molecular.cpp b/src/MOLECULE/atom_vec_molecular.cpp
index ac26f00c49c921239362d92b4ec6e622b6d19f52..f1a09990ac3d970375978c10b6fa39c3da8c4e34 100644
--- a/src/MOLECULE/atom_vec_molecular.cpp
+++ b/src/MOLECULE/atom_vec_molecular.cpp
@@ -154,9 +154,11 @@ void AtomVecMolecular::grow_reset()
improper_atom3 = atom->improper_atom3; improper_atom4 = atom->improper_atom4;
}
-/* ---------------------------------------------------------------------- */
+/* ----------------------------------------------------------------------
+ copy atom I info to atom J
+------------------------------------------------------------------------- */
-void AtomVecMolecular::copy(int i, int j)
+void AtomVecMolecular::copy(int i, int j, int delflag)
{
int k;
diff --git a/src/MOLECULE/atom_vec_molecular.h b/src/MOLECULE/atom_vec_molecular.h
index 6249b8bf1e9d59f9ae3bda57c0aeb9e8098a3bd7..6b3fa4003823c812c0f6f33b104531570dc9cf23 100644
--- a/src/MOLECULE/atom_vec_molecular.h
+++ b/src/MOLECULE/atom_vec_molecular.h
@@ -29,7 +29,7 @@ class AtomVecMolecular : public AtomVec {
AtomVecMolecular(class LAMMPS *, int, char **);
void grow(int);
void grow_reset();
- void copy(int, int);
+ void copy(int, int, int);
int pack_comm(int, int *, double *, int, int *);
int pack_comm_vel(int, int *, double *, int, int *);
void unpack_comm(int, int, double *);
diff --git a/src/PERI/atom_vec_peri.cpp b/src/PERI/atom_vec_peri.cpp
index b4d94e3fa8a3f76b39977e217d7b1fcd01b3fa68..d85853d114be03479826d46ac90996177bb42bec 100644
--- a/src/PERI/atom_vec_peri.cpp
+++ b/src/PERI/atom_vec_peri.cpp
@@ -96,9 +96,11 @@ void AtomVecPeri::grow_reset()
s0 = atom->s0; x0 = atom->x0;
}
-/* ---------------------------------------------------------------------- */
+/* ----------------------------------------------------------------------
+ copy atom I info to atom J
+------------------------------------------------------------------------- */
-void AtomVecPeri::copy(int i, int j)
+void AtomVecPeri::copy(int i, int j, int delflag)
{
tag[j] = tag[i];
type[j] = type[i];
diff --git a/src/PERI/atom_vec_peri.h b/src/PERI/atom_vec_peri.h
index 4b14c04b42d7a6c9e4bb0015091632bade547baa..ca59123a1cad46b9073895d63381a8369ddc50f0 100755
--- a/src/PERI/atom_vec_peri.h
+++ b/src/PERI/atom_vec_peri.h
@@ -29,7 +29,7 @@ class AtomVecPeri : public AtomVec {
AtomVecPeri(class LAMMPS *, int, char **);
void grow(int);
void grow_reset();
- void copy(int, int);
+ void copy(int, int, int);
int pack_comm(int, int *, double *, int, int *);
int pack_comm_vel(int, int *, double *, int, int *);
int pack_comm_hybrid(int, int *, double *);
diff --git a/src/PERI/pair_peri_lps.cpp b/src/PERI/pair_peri_lps.cpp
index c4a5cdcae3f55d060c7fefa373131c962aca6046..6a0be34ba7a9fe053ffd0600f3c2abcf5d04e041 100644
--- a/src/PERI/pair_peri_lps.cpp
+++ b/src/PERI/pair_peri_lps.cpp
@@ -423,12 +423,10 @@ void PairPeriLPS::init_style()
{
// error checks
+ if (!atom->peri_flag) error->all("Pair style peri requires atom style peri");
if (atom->map_style == 0)
error->all("Pair peri requires an atom map, see atom_modify");
- if (atom->style_match("peri") == 0)
- error->all("Pair style peri_lps requires atom style peri");
-
if (domain->lattice == NULL)
error->all("Pair peri requires a lattice be defined");
if (domain->lattice->xlattice != domain->lattice->ylattice ||
diff --git a/src/PERI/pair_peri_pmb.cpp b/src/PERI/pair_peri_pmb.cpp
index 3ee46cb938463a00e94c5eca7064888ac0d9a5af..6299884c3dd667c0515001c7490003a9c283b907 100644
--- a/src/PERI/pair_peri_pmb.cpp
+++ b/src/PERI/pair_peri_pmb.cpp
@@ -360,12 +360,10 @@ void PairPeriPMB::init_style()
{
// error checks
+ if (!atom->peri_flag) error->all("Pair style peri requires atom style peri");
if (atom->map_style == 0)
error->all("Pair peri requires an atom map, see atom_modify");
- if (atom->style_match("peri") == 0)
- error->all("Pair style peri_pmb requires atom style peri");
-
if (domain->lattice == NULL)
error->all("Pair peri requires a lattice be defined");
if (domain->lattice->xlattice != domain->lattice->ylattice ||
diff --git a/src/SRD/fix_srd.cpp b/src/SRD/fix_srd.cpp
index ec8716f403e3b55c2496da71d16479656f90cfec..576a43e6adf973b639a60fcf89333aea7e264c43 100644
--- a/src/SRD/fix_srd.cpp
+++ b/src/SRD/fix_srd.cpp
@@ -20,7 +20,7 @@
#include "stdlib.h"
#include "fix_srd.h"
#include "atom.h"
-#include "atom_vec.h"
+#include "atom_vec_ellipsoid.h"
#include "group.h"
#include "update.h"
#include "force.h"
@@ -233,6 +233,10 @@ FixSRD::FixSRD(LAMMPS *lmp, int narg, char **arg) : Fix(lmp, narg, arg)
srd_bin_temp = 0.0;
srd_bin_count = 0;
+ // atom style pointers to particles that store extra info
+
+ avec_ellipsoid = (AtomVecEllipsoid *) atom->style_match("ellipsoid");
+
// fix parameters
if (collidestyle == SLIP) comm_reverse = 3;
@@ -288,12 +292,8 @@ void FixSRD::init()
if (force->newton_pair == 0) error->all("Fix srd requires newton pair on");
if (bigexist && comm->ghost_velocity == 0)
error->all("Fix srd requires ghost atoms store velocity");
-
- if (bigexist && !atom->sphere_flag && !atom->ellipsoid_flag)
- error->all("Fix SRD requires atom style sphere or ellipsoid");
if (bigexist && collidestyle == NOSLIP && !atom->torque_flag)
error->all("Fix SRD no-slip requires atom attribute torque");
-
if (initflag && update->dt != dt_big)
error->all("Cannot change timestep once fix srd is setup");
@@ -2092,8 +2092,9 @@ void FixSRD::parameterize()
// big particle must either have radius > 0 or shape > 0 defined
// apply radfactor at end
+ AtomVecEllipsoid::Bonus *ebonus = avec_ellipsoid->bonus;
double *radius = atom->radius;
- double **shape = atom->shape;
+ int *ellipsoid = atom->ellipsoid;
int *mask = atom->mask;
int nlocal = atom->nlocal;
@@ -2106,15 +2107,15 @@ void FixSRD::parameterize()
if (radius && radius[i] > 0.0) {
maxbigdiam = MAX(maxbigdiam,2.0*radius[i]);
minbigdiam = MIN(minbigdiam,2.0*radius[i]);
- } else if (shape && shape[i][0] > 0.0) {
- maxbigdiam = MAX(maxbigdiam,2.0*shape[i][0]);
- maxbigdiam = MAX(maxbigdiam,2.0*shape[i][1]);
- maxbigdiam = MAX(maxbigdiam,2.0*shape[i][2]);
- minbigdiam = MIN(minbigdiam,2.0*shape[i][0]);
- minbigdiam = MIN(minbigdiam,2.0*shape[i][1]);
- minbigdiam = MIN(minbigdiam,2.0*shape[i][2]);
- if (shape[i][0] != shape[i][1] || shape[i][0] != shape[i][2])
- any_ellipsoids = 1;
+ } else if (ellipsoid && ellipsoid[i] >= 0) {
+ any_ellipsoids = 1;
+ double *shape = ebonus[ellipsoid[i]].shape;
+ maxbigdiam = MAX(maxbigdiam,2.0*shape[0]);
+ maxbigdiam = MAX(maxbigdiam,2.0*shape[1]);
+ maxbigdiam = MAX(maxbigdiam,2.0*shape[2]);
+ minbigdiam = MIN(minbigdiam,2.0*shape[0]);
+ minbigdiam = MIN(minbigdiam,2.0*shape[1]);
+ minbigdiam = MIN(minbigdiam,2.0*shape[2]);
} else
error->one("Big particle in fix srd cannot be point particle");
}
@@ -2171,8 +2172,6 @@ void FixSRD::parameterize()
temperature_srd = force->mvv2e *
(lamda/dt_srd)*(lamda/dt_srd) * mass_srd/force->boltz;
- printf("AAA %g %g\n",mass_srd,lamda);
-
// vmax = maximum velocity of an SRD particle
// dmax = maximum distance an SRD can move = 4*lamda = vmax * dt_srd
@@ -2191,16 +2190,20 @@ void FixSRD::parameterize()
if (mask[i] & biggroupbit) {
if (radius && radius[i] > 0.0)
volbig += 4.0/3.0*PI*radius[i]*radius[i]*radius[i];
- else if (shape && shape[i][0] > 0.0)
- volbig += 4.0/3.0*PI * shape[i][0]*shape[i][1]*shape[i][2];
+ else if (ellipsoid && ellipsoid[i] >= 0) {
+ double *shape = ebonus[ellipsoid[i]].shape;
+ volbig += 4.0/3.0*PI * shape[0]*shape[1]*shape[2];
+ }
}
} else {
for (int i = 0; i < nlocal; i++)
if (mask[i] & biggroupbit) {
if (radius && radius[i] > 0.0)
volbig += PI*radius[i]*radius[i];
- else if (shape && shape[i][0] > 0.0)
- volbig += PI*shape[i][0]*shape[i][1];
+ else if (ellipsoid && ellipsoid[i] >= 0) {
+ double *shape = ebonus[ellipsoid[i]].shape;
+ volbig += PI*shape[0]*shape[1];
+ }
}
}
@@ -2401,9 +2404,11 @@ void FixSRD::big_static()
{
int i;
double rad,arad,brad,crad;
+ double *shape;
+ AtomVecEllipsoid::Bonus *ebonus = avec_ellipsoid->bonus;
double *radius = atom->radius;
- double **shape = atom->shape;
+ int *ellipsoid = atom->ellipsoid;
int *type = atom->type;
double skinhalf = 0.5 * neighbor->skin;
@@ -2416,11 +2421,12 @@ void FixSRD::big_static()
biglist[k].radius = rad;
biglist[k].radsq = rad*rad;
biglist[k].cutbinsq = (rad+skinhalf) * (rad+skinhalf);
- } else if (shape && shape[i][0] > 0.0) {
+ } else if (ellipsoid && ellipsoid[i] >= 0) {
+ shape = ebonus[ellipsoid[i]].shape;
biglist[k].type = ELLIPSOID;
- arad = radfactor*shape[i][0];
- brad = radfactor*shape[i][1];
- crad = radfactor*shape[i][2];
+ arad = radfactor*shape[0];
+ brad = radfactor*shape[1];
+ crad = radfactor*shape[2];
biglist[k].aradsqinv = 1.0/(arad*arad);
biglist[k].bradsqinv = 1.0/(brad*brad);
biglist[k].cradsqinv = 1.0/(crad*crad);
@@ -2441,12 +2447,13 @@ void FixSRD::big_static()
void FixSRD::big_dynamic()
{
int i;
+ double *shape,*quat;
+ AtomVecEllipsoid::Bonus *ebonus = avec_ellipsoid->bonus;
double **omega = atom->omega;
double **angmom = atom->angmom;
- double **quat = atom->quat;
- double **shape = atom->shape;
double *rmass = atom->rmass;
+ int *ellipsoid = atom->ellipsoid;
for (int k = 0; k < nbig; k++) {
i = biglist[k].index;
@@ -2463,9 +2470,11 @@ void FixSRD::big_dynamic()
// calculate ex,ey,ez and omega from quaternion and angmom
} else if (biglist[k].type == ELLIPSOID) {
- exyz_from_q(quat[i],biglist[k].ex,biglist[k].ey,biglist[k].ez);
+ shape = ebonus[ellipsoid[i]].shape;
+ quat = ebonus[ellipsoid[i]].quat;
+ exyz_from_q(quat,biglist[k].ex,biglist[k].ey,biglist[k].ez);
omega_from_mq(angmom[i],biglist[k].ex,biglist[k].ey,biglist[k].ez,
- rmass[i],shape[i],biglist[k].omega);
+ rmass[i],shape,biglist[k].omega);
}
}
}
diff --git a/src/SRD/fix_srd.h b/src/SRD/fix_srd.h
index b7d1fd188b5069898e91eededbfb9260ad35c36f..8552078e6d0568e8284e2c428ec8ae579d03caad 100644
--- a/src/SRD/fix_srd.h
+++ b/src/SRD/fix_srd.h
@@ -63,6 +63,8 @@ class FixSRD : public Fix {
double **fwall;
double walltrigger;
+ class AtomVecEllipsoid *avec_ellipsoid;
+
// for orthogonal box, these are in box units
// for triclinic box, these are in lamda units
diff --git a/src/USER-EFF/atom_vec_electron.cpp b/src/USER-EFF/atom_vec_electron.cpp
index 0a13aa1a761157f4d48c9c499d3290687d570398..e5c47745ab58ef2f6711a5d8e72d13ce9b1931c6 100644
--- a/src/USER-EFF/atom_vec_electron.cpp
+++ b/src/USER-EFF/atom_vec_electron.cpp
@@ -97,9 +97,11 @@ void AtomVecElectron::grow_reset()
eradius = atom->eradius; ervel = atom->ervel; erforce = atom->erforce;
}
-/* ---------------------------------------------------------------------- */
+/* ----------------------------------------------------------------------
+ copy atom I info to atom J
+------------------------------------------------------------------------- */
-void AtomVecElectron::copy(int i, int j)
+void AtomVecElectron::copy(int i, int j, int delflag)
{
tag[j] = tag[i];
type[j] = type[i];
diff --git a/src/USER-EFF/atom_vec_electron.h b/src/USER-EFF/atom_vec_electron.h
index 3ae836ae2397a3afdebcd733058052322a7e622c..b9372f62a633683ffb34e32db5385cf6fba00baf 100644
--- a/src/USER-EFF/atom_vec_electron.h
+++ b/src/USER-EFF/atom_vec_electron.h
@@ -30,7 +30,7 @@ class AtomVecElectron : public AtomVec {
~AtomVecElectron() {}
void grow(int);
void grow_reset();
- void copy(int, int);
+ void copy(int, int, int);
int pack_comm(int, int *, double *, int, int *);
int pack_comm_vel(int, int *, double *, int, int *);
int pack_comm_hybrid(int, int *, double *);
diff --git a/src/atom.cpp b/src/atom.cpp
index 81ca000a40c5efaa210f2be0c6ac7fec1f5a853d..abae1aabe7b0bf2a7fd60d993d86421f288945bf 100644
--- a/src/atom.cpp
+++ b/src/atom.cpp
@@ -19,6 +19,7 @@
#include "atom.h"
#include "style_atom.h"
#include "atom_vec.h"
+#include "atom_vec_ellipsoid.h"
#include "comm.h"
#include "neighbor.h"
#include "force.h"
@@ -70,11 +71,11 @@ Atom::Atom(LAMMPS *lmp) : Pointers(lmp)
molecule = NULL;
q = NULL;
mu = NULL;
- quat = omega = angmom = torque = shape = NULL;
+ omega = angmom = torque = NULL;
radius = rmass = NULL;
vfrac = s0 = NULL;
x0 = NULL;
-
+ ellipsoid = NULL;
spin = NULL;
eradius = ervel = erforce = NULL;
@@ -102,8 +103,7 @@ Atom::Atom(LAMMPS *lmp) : Pointers(lmp)
sphere_flag = ellipsoid_flag = peri_flag = electron_flag = 0;
molecule_flag = q_flag = mu_flag = 0;
- rmass_flag = radius_flag = omega_flag = torque_flag = 0;
- quat_flag = shape_flag = angmom_flag = 0;
+ rmass_flag = radius_flag = omega_flag = torque_flag = angmom_flag = 0;
vfrac_flag = spin_flag = eradius_flag = ervel_flag = erforce_flag = 0;
// ntype-length arrays
@@ -167,18 +167,15 @@ Atom::~Atom()
memory->destroy(q);
memory->destroy(mu);
- memory->destroy(quat);
- memory->destroy(shape);
memory->destroy(omega);
memory->destroy(angmom);
memory->destroy(torque);
-
memory->destroy(radius);
memory->destroy(rmass);
memory->destroy(vfrac);
memory->destroy(s0);
memory->destroy(x0);
-
+ memory->destroy(ellipsoid);
memory->destroy(spin);
memory->destroy(eradius);
memory->destroy(ervel);
@@ -256,8 +253,7 @@ void Atom::create_avec(const char *style, int narg, char **arg)
sphere_flag = ellipsoid_flag = peri_flag = electron_flag = 0;
molecule_flag = q_flag = mu_flag = 0;
- rmass_flag = radius_flag = omega_flag = torque_flag = 0;
- quat_flag = shape_flag = angmom_flag = 0;
+ rmass_flag = radius_flag = omega_flag = torque_flag = angmom_flag = 0;
vfrac_flag = spin_flag = eradius_flag = ervel_flag = erforce_flag = 0;
avec = new_avec(style,narg,arg);
@@ -329,19 +325,20 @@ void Atom::setup()
}
/* ----------------------------------------------------------------------
- return 1 if style matches atom style or hybrid sub-style
- else return 0
+ return ptr to AtomVec class if matches style or to matching hybrid sub-class
+ return NULL if no match
------------------------------------------------------------------------- */
-int Atom::style_match(const char *style)
+AtomVec *Atom::style_match(const char *style)
{
- if (strcmp(atom_style,style) == 0) return 1;
+ if (strcmp(atom_style,style) == 0) return avec;
else if (strcmp(atom_style,"hybrid") == 0) {
AtomVecHybrid *avec_hybrid = (AtomVecHybrid *) avec;
for (int i = 0; i < avec_hybrid->nstyles; i++)
- if (strcmp(avec_hybrid->keywords[i],style) == 0) return 1;
+ if (strcmp(avec_hybrid->keywords[i],style) == 0)
+ return avec_hybrid->styles[i];
}
- return 0;
+ return NULL;
}
/* ----------------------------------------------------------------------
@@ -847,6 +844,53 @@ void Atom::data_vels(int n, char *buf)
delete [] values;
}
+/* ----------------------------------------------------------------------
+ unpack n lines from atom-style specific section of data file
+ check that atom IDs are > 0 and <= map_tag_max
+ call style-specific routine to parse line
+------------------------------------------------------------------------- */
+
+void Atom::data_bonus(int n, char *buf, AtomVec *avec_bonus)
+{
+ int j,m,tagdata;
+ char *next;
+
+ next = strchr(buf,'\n');
+ *next = '\0';
+ int nwords = count_words(buf);
+ *next = '\n';
+
+ if (nwords != avec_bonus->size_data_bonus)
+ error->all("Incorrect bonus data format in data file");
+
+ char **values = new char*[nwords];
+
+ // loop over lines of bonus atom data
+ // tokenize the line into values
+ // if I own atom tag, unpack its values
+
+ for (int i = 0; i < n; i++) {
+ next = strchr(buf,'\n');
+
+ values[0] = strtok(buf," \t\n\r\f");
+ for (j = 1; j < nwords; j++)
+ values[j] = strtok(NULL," \t\n\r\f");
+
+ tagdata = atoi(values[0]);
+ if (tagdata <= 0 || tagdata > map_tag_max)
+ error->one("Invalid atom ID in Bonus section of data file");
+
+ // ok to call child's data_atom_bonus() method thru parent avec_bonus,
+ // since data_bonus() was called with child ptr, and method is virtual
+
+ if ((m = map(tagdata)) >= 0) avec_bonus->data_atom_bonus(m,&values[1]);
+
+ buf = next + 1;
+ }
+
+ delete [] values;
+}
+
/* ----------------------------------------------------------------------
check that atom IDs are > 0 and <= map_tag_max
------------------------------------------------------------------------- */
@@ -1170,23 +1214,32 @@ int Atom::radius_consistency(int itype, double &rad)
/* ----------------------------------------------------------------------
check that shape of all particles of itype are the same
return 1 if true, else return 0
- also return the radius value for that type
+ also return the 3 shape params for itype
------------------------------------------------------------------------- */
int Atom::shape_consistency(int itype,
double &shapex, double &shapey, double &shapez)
{
- double one[3];
- one[0] = one[1] = one[2] = -1.0;
+ double zero[3] = {0.0, 0.0, 0.0};
+ double one[3] = {-1.0, -1.0, -1.0};
+ double *shape;
+
+ AtomVecEllipsoid *avec_ellipsoid =
+ (AtomVecEllipsoid *) style_match("ellipsoid");
+ AtomVecEllipsoid::Bonus *bonus = avec_ellipsoid->bonus;
+
int flag = 0;
for (int i = 0; i < nlocal; i++) {
if (type[i] != itype) continue;
+ if (ellipsoid[i] < 0) shape = zero;
+ else shape = bonus[ellipsoid[i]].shape;
+
if (one[0] < 0.0) {
- one[0] = shape[i][0];
- one[1] = shape[i][1];
- one[2] = shape[i][2];
- } else if (one[0] != shape[i][0] || one[1] != shape[i][1] ||
- one[2] != shape[i][2]) flag = 1;
+ one[0] = shape[0];
+ one[1] = shape[1];
+ one[2] = shape[2];
+ } else if (one[0] != shape[0] || one[1] != shape[1] || one[2] != shape[2])
+ flag = 1;
}
int flagall;
@@ -1223,9 +1276,9 @@ void Atom::first_reorder()
for (int i = 0; i < nlocal; i++) {
if (mask[i] & bitmask && i > nfirst) {
- avec->copy(i,nlocal);
- avec->copy(nfirst,i);
- avec->copy(nlocal,nfirst);
+ avec->copy(i,nlocal,0);
+ avec->copy(nfirst,i,0);
+ avec->copy(nlocal,nfirst,0);
while (nfirst < nlocal && mask[nfirst] & bitmask) nfirst++;
}
}
@@ -1310,13 +1363,13 @@ void Atom::sort()
for (i = 0; i < nlocal; i++) {
if (current[i] == permute[i]) continue;
- avec->copy(i,nlocal);
+ avec->copy(i,nlocal,0);
empty = i;
while (permute[empty] != i) {
- avec->copy(permute[empty],empty);
+ avec->copy(permute[empty],empty,0);
empty = current[empty] = permute[empty];
}
- avec->copy(nlocal,empty);
+ avec->copy(nlocal,empty,0);
current[empty] = permute[empty];
}
diff --git a/src/atom.h b/src/atom.h
index 47bad27ffe2e0c36eb1723296408f6d2259b74ed..c0e5dafa603ec2de81c1d80192262da90b12cc79 100644
--- a/src/atom.h
+++ b/src/atom.h
@@ -48,10 +48,10 @@ class Atom : protected Pointers {
int *molecule;
double *q,**mu;
- double **quat,**omega,**angmom,**torque,**shape;
+ double **omega,**angmom,**torque;
double *radius,*rmass,*vfrac,*s0;
double **x0;
-
+ int *ellipsoid;
int *spin;
double *eradius,*ervel,*erforce;
@@ -81,8 +81,7 @@ class Atom : protected Pointers {
int sphere_flag,ellipsoid_flag,peri_flag,electron_flag;
int molecule_flag,q_flag,mu_flag;
- int rmass_flag,radius_flag,omega_flag,torque_flag;
- int quat_flag,shape_flag,angmom_flag;
+ int rmass_flag,radius_flag,omega_flag,torque_flag,angmom_flag;
int vfrac_flag,spin_flag,eradius_flag,ervel_flag,erforce_flag;
// extra peratom info in restart file destined for fix & diag
@@ -120,7 +119,7 @@ class Atom : protected Pointers {
void init();
void setup();
- int style_match(const char *);
+ class AtomVec *style_match(const char *);
void modify_params(int, char **);
void tag_extend();
int tag_consecutive();
@@ -130,6 +129,8 @@ class Atom : protected Pointers {
void data_atoms(int, char *);
void data_vels(int, char *);
+ void data_bonus(int, char *, class AtomVec *);
+
void data_bonds(int, char *);
void data_angles(int, char *);
void data_dihedrals(int, char *);
diff --git a/src/atom_vec.cpp b/src/atom_vec.cpp
index 24037a2c7a86aa382af84caa56ee52bec450d90d..d39a98539ef1ce52df1e8059b6badd8d3a47c21a 100644
--- a/src/atom_vec.cpp
+++ b/src/atom_vec.cpp
@@ -25,6 +25,7 @@ AtomVec::AtomVec(LAMMPS *lmp, int narg, char **arg) : Pointers(lmp)
nmax = 0;
bonds_allow = angles_allow = dihedrals_allow = impropers_allow = 0;
mass_type = dipole_type = 0;
+ size_data_bonus = 0;
}
/* ----------------------------------------------------------------------
diff --git a/src/atom_vec.h b/src/atom_vec.h
index 594931507e6b5d093007e6ed91277ba76189dfdd..a437d14a6b45386f63dcf2e464ea1e8522e4e409 100644
--- a/src/atom_vec.h
+++ b/src/atom_vec.h
@@ -35,6 +35,7 @@ class AtomVec : protected Pointers {
int size_velocity; // # of velocity based quantities
int size_data_atom; // number of values in Atom line
int size_data_vel; // number of values in Velocity line
+ int size_data_bonus; // number of values in Bonus line
int xcol_data; // column (1-N) where x is in Atom line
AtomVec(class LAMMPS *, int, char **);
@@ -43,7 +44,8 @@ class AtomVec : protected Pointers {
virtual void grow(int) = 0;
virtual void grow_reset() = 0;
- virtual void copy(int, int) = 0;
+ virtual void copy(int, int, int) = 0;
+ virtual void clear_bonus() {}
virtual int pack_comm(int, int *, double *, int, int *) = 0;
virtual int pack_comm_vel(int, int *, double *, int, int *) = 0;
@@ -73,6 +75,7 @@ class AtomVec : protected Pointers {
virtual void create_atom(int, double *) = 0;
virtual void data_atom(double *, int, char **) = 0;
+ virtual void data_atom_bonus(int, char **) {}
virtual int data_atom_hybrid(int, char **) {return 0;}
virtual void data_vel(int, char **);
virtual int data_vel_hybrid(int, char **) {return 0;}
diff --git a/src/atom_vec_atomic.cpp b/src/atom_vec_atomic.cpp
index 0b5d36d15c7fecc3f12097fd70279d2265ffa05a..fd33d6a41d13e1c588b927e19839a4bf499643d8 100644
--- a/src/atom_vec_atomic.cpp
+++ b/src/atom_vec_atomic.cpp
@@ -81,9 +81,11 @@ void AtomVecAtomic::grow_reset()
x = atom->x; v = atom->v; f = atom->f;
}
-/* ---------------------------------------------------------------------- */
+/* ----------------------------------------------------------------------
+ copy atom I info to atom J
+------------------------------------------------------------------------- */
-void AtomVecAtomic::copy(int i, int j)
+void AtomVecAtomic::copy(int i, int j, int delflag)
{
tag[j] = tag[i];
type[j] = type[i];
diff --git a/src/atom_vec_atomic.h b/src/atom_vec_atomic.h
index a0700a26d248b9fd8eb32cd76e725a4af687580e..187f4c22f8077b548edb56d55bb08c2173a35cd1 100644
--- a/src/atom_vec_atomic.h
+++ b/src/atom_vec_atomic.h
@@ -30,7 +30,7 @@ class AtomVecAtomic : public AtomVec {
~AtomVecAtomic() {}
void grow(int);
void grow_reset();
- void copy(int, int);
+ void copy(int, int, int);
int pack_comm(int, int *, double *, int, int *);
int pack_comm_vel(int, int *, double *, int, int *);
void unpack_comm(int, int, double *);
diff --git a/src/atom_vec_charge.cpp b/src/atom_vec_charge.cpp
index 5088361cbde461c11b483fb51fa0dc5e1ece9f8e..a2e9720d33c4d8a8b5d708b197f842b697804f78 100644
--- a/src/atom_vec_charge.cpp
+++ b/src/atom_vec_charge.cpp
@@ -86,9 +86,11 @@ void AtomVecCharge::grow_reset()
q = atom->q;
}
-/* ---------------------------------------------------------------------- */
+/* ----------------------------------------------------------------------
+ copy atom I info to atom J
+------------------------------------------------------------------------- */
-void AtomVecCharge::copy(int i, int j)
+void AtomVecCharge::copy(int i, int j, int delflag)
{
tag[j] = tag[i];
type[j] = type[i];
diff --git a/src/atom_vec_charge.h b/src/atom_vec_charge.h
index a25641f76dc54b94b6adafd2ccad5edd57c8aa41..a19d7d6cde09eb966d95fc6db0bf036dc95a17ec 100644
--- a/src/atom_vec_charge.h
+++ b/src/atom_vec_charge.h
@@ -29,7 +29,7 @@ class AtomVecCharge : public AtomVec {
AtomVecCharge(class LAMMPS *, int, char **);
void grow(int);
void grow_reset();
- void copy(int, int);
+ void copy(int, int, int);
int pack_comm(int, int *, double *, int, int *);
int pack_comm_vel(int, int *, double *, int, int *);
void unpack_comm(int, int, double *);
diff --git a/src/atom_vec_hybrid.cpp b/src/atom_vec_hybrid.cpp
index d4e148245c531ba5a948f49bb465f184625f4d79..a7f48066f0d8f402a31fc2cfb2b1d0b27def2b11 100644
--- a/src/atom_vec_hybrid.cpp
+++ b/src/atom_vec_hybrid.cpp
@@ -157,14 +157,14 @@ void AtomVecHybrid::grow_reset()
}
/* ----------------------------------------------------------------------
- copy array values for all sub-styles
+ copy atom I info to atom J for all sub-styles
------------------------------------------------------------------------- */
-void AtomVecHybrid::copy(int i, int j)
+void AtomVecHybrid::copy(int i, int j, int delflag)
{
int tmp = atom->nextra_grow;
atom->nextra_grow = 0;
- for (int k = 0; k < nstyles; k++) styles[k]->copy(i,j);
+ for (int k = 0; k < nstyles; k++) styles[k]->copy(i,j,delflag);
atom->nextra_grow = tmp;
if (atom->nextra_grow)
@@ -174,6 +174,13 @@ void AtomVecHybrid::copy(int i, int j)
/* ---------------------------------------------------------------------- */
+void AtomVecHybrid::clear_bonus()
+{
+ for (int k = 0; k < nstyles; k++) styles[k]->clear_bonus();
+}
+
+/* ---------------------------------------------------------------------- */
+
int AtomVecHybrid::pack_comm(int n, int *list, double *buf,
int pbc_flag, int *pbc)
{
diff --git a/src/atom_vec_hybrid.h b/src/atom_vec_hybrid.h
index 08f03c361264b154aa725948d82a27c0b7bc3653..8d157b8aba740633ccf299e471d60a65cb2c852a 100644
--- a/src/atom_vec_hybrid.h
+++ b/src/atom_vec_hybrid.h
@@ -35,7 +35,8 @@ class AtomVecHybrid : public AtomVec {
void init();
void grow(int);
void grow_reset();
- void copy(int, int);
+ void copy(int, int, int);
+ void clear_bonus();
int pack_comm(int, int *, double *, int, int *);
int pack_comm_vel(int, int *, double *, int, int *);
void unpack_comm(int, int, double *);
diff --git a/src/atom_vec_sphere.cpp b/src/atom_vec_sphere.cpp
index 1aa886c7ed23711e9d2ee9e1a3fb1debe25e8fac..ab53b79c5df7bcd2924587ab7bca7ee34eb2c62f 100644
--- a/src/atom_vec_sphere.cpp
+++ b/src/atom_vec_sphere.cpp
@@ -121,9 +121,11 @@ void AtomVecSphere::grow_reset()
omega = atom->omega; torque = atom->torque;
}
-/* ---------------------------------------------------------------------- */
+/* ----------------------------------------------------------------------
+ copy atom I info to atom J
+------------------------------------------------------------------------- */
-void AtomVecSphere::copy(int i, int j)
+void AtomVecSphere::copy(int i, int j, int delflag)
{
tag[j] = tag[i];
type[j] = type[i];
diff --git a/src/atom_vec_sphere.h b/src/atom_vec_sphere.h
index 566db3ff9681db6b25af768d67d8aacea86578f5..81570296155b5b35ea0822744e041f987d840274 100644
--- a/src/atom_vec_sphere.h
+++ b/src/atom_vec_sphere.h
@@ -31,7 +31,7 @@ class AtomVecSphere : public AtomVec {
void init();
void grow(int);
void grow_reset();
- void copy(int, int);
+ void copy(int, int, int);
int pack_comm(int, int *, double *, int, int *);
int pack_comm_vel(int, int *, double *, int, int *);
int pack_comm_hybrid(int, int *, double *);
diff --git a/src/comm.cpp b/src/comm.cpp
index 4a3c942a45df7f3e5402557d4c5e871a347194e8..e97d3f34741bf325865e48871cfdeb802b9a735f 100644
--- a/src/comm.cpp
+++ b/src/comm.cpp
@@ -566,7 +566,7 @@ void Comm::exchange()
if (x[i][dim] < lo || x[i][dim] >= hi) {
if (nsend > maxsend) grow_send(nsend,1);
nsend += avec->pack_exchange(i,&buf_send[nsend]);
- avec->copy(nlocal-1,i);
+ avec->copy(nlocal-1,i,1);
nlocal--;
} else i++;
}
@@ -643,9 +643,10 @@ void Comm::borders()
MPI_Status status;
AtomVec *avec = atom->avec;
- // clear old ghosts
+ // clear old ghosts and any ghost bonus data internal to AtomVec
atom->nghost = 0;
+ atom->avec->clear_bonus();
// do swaps over all 3 dimensions
diff --git a/src/compute_property_atom.cpp b/src/compute_property_atom.cpp
index 2949b9f469b8324eae296527e7f3b00588aaad00..a47a216e217f960ff399aa52d59582d697042cf7 100644
--- a/src/compute_property_atom.cpp
+++ b/src/compute_property_atom.cpp
@@ -14,6 +14,7 @@
#include "string.h"
#include "compute_property_atom.h"
#include "atom.h"
+#include "atom_vec_ellipsoid.h"
#include "update.h"
#include "domain.h"
#include "memory.h"
@@ -167,39 +168,39 @@ ComputePropertyAtom::ComputePropertyAtom(LAMMPS *lmp, int narg, char **arg) :
pack_choice[i] = &ComputePropertyAtom::pack_angmomz;
} else if (strcmp(arg[iarg],"shapex") == 0) {
- if (!atom->shape_flag)
- error->all("Compute property/atom for "
- "atom property that isn't allocated");
+ avec = (AtomVecEllipsoid *) atom->style_match("ellipsoid");
+ if (!avec) error->all("Compute property/atom for "
+ "atom property that isn't allocated");
pack_choice[i] = &ComputePropertyAtom::pack_shapex;
} else if (strcmp(arg[iarg],"shapey") == 0) {
- if (!atom->shape_flag)
- error->all("Compute property/atom for "
- "atom property that isn't allocated");
+ avec = (AtomVecEllipsoid *) atom->style_match("ellipsoid");
+ if (!avec) error->all("Compute property/atom for "
+ "atom property that isn't allocated");
pack_choice[i] = &ComputePropertyAtom::pack_shapey;
} else if (strcmp(arg[iarg],"shapez") == 0) {
- if (!atom->shape_flag)
- error->all("Compute property/atom for "
- "atom property that isn't allocated");
+ avec = (AtomVecEllipsoid *) atom->style_match("ellipsoid");
+ if (!avec) error->all("Compute property/atom for "
+ "atom property that isn't allocated");
pack_choice[i] = &ComputePropertyAtom::pack_shapez;
} else if (strcmp(arg[iarg],"quatw") == 0) {
- if (!atom->quat_flag)
- error->all("Compute property/atom for "
- "atom property that isn't allocated");
+ avec = (AtomVecEllipsoid *) atom->style_match("ellipsoid");
+ if (!avec) error->all("Compute property/atom for "
+ "atom property that isn't allocated");
pack_choice[i] = &ComputePropertyAtom::pack_quatw;
} else if (strcmp(arg[iarg],"quati") == 0) {
- if (!atom->quat_flag)
- error->all("Compute property/atom for "
- "atom property that isn't allocated");
+ avec = (AtomVecEllipsoid *) atom->style_match("ellipsoid");
+ if (!avec) error->all("Compute property/atom for "
+ "atom property that isn't allocated");
pack_choice[i] = &ComputePropertyAtom::pack_quati;
} else if (strcmp(arg[iarg],"quatj") == 0) {
- if (!atom->quat_flag)
- error->all("Compute property/atom for "
- "atom property that isn't allocated");
+ avec = (AtomVecEllipsoid *) atom->style_match("ellipsoid");
+ if (!avec) error->all("Compute property/atom for "
+ "atom property that isn't allocated");
pack_choice[i] = &ComputePropertyAtom::pack_quatj;
} else if (strcmp(arg[iarg],"quatk") == 0) {
- if (!atom->quat_flag)
- error->all("Compute property/atom for "
- "atom property that isn't allocated");
+ avec = (AtomVecEllipsoid *) atom->style_match("ellipsoid");
+ if (!avec) error->all("Compute property/atom for "
+ "atom property that isn't allocated");
pack_choice[i] = &ComputePropertyAtom::pack_quatk;
} else if (strcmp(arg[iarg],"tqx") == 0) {
if (!atom->torque_flag)
@@ -973,12 +974,14 @@ void ComputePropertyAtom::pack_angmomz(int n)
void ComputePropertyAtom::pack_shapex(int n)
{
- double **shape = atom->shape;
+ AtomVecEllipsoid::Bonus *bonus = avec->bonus;
+ int *ellipsoid = atom->ellipsoid;
int *mask = atom->mask;
int nlocal = atom->nlocal;
for (int i = 0; i < nlocal; i++) {
- if (mask[i] & groupbit) buf[n] = shape[i][0];
+ if ((mask[i] & groupbit) && ellipsoid[i] >= 0)
+ buf[n] = bonus[ellipsoid[i]].shape[0];
else buf[n] = 0.0;
n += nvalues;
}
@@ -988,12 +991,14 @@ void ComputePropertyAtom::pack_shapex(int n)
void ComputePropertyAtom::pack_shapey(int n)
{
- double **shape = atom->shape;
+ AtomVecEllipsoid::Bonus *bonus = avec->bonus;
+ int *ellipsoid = atom->ellipsoid;
int *mask = atom->mask;
int nlocal = atom->nlocal;
for (int i = 0; i < nlocal; i++) {
- if (mask[i] & groupbit) buf[n] = shape[i][1];
+ if ((mask[i] & groupbit) && ellipsoid[i] >= 0)
+ buf[n] = bonus[ellipsoid[i]].shape[1];
else buf[n] = 0.0;
n += nvalues;
}
@@ -1003,12 +1008,14 @@ void ComputePropertyAtom::pack_shapey(int n)
void ComputePropertyAtom::pack_shapez(int n)
{
- double **shape = atom->shape;
+ AtomVecEllipsoid::Bonus *bonus = avec->bonus;
+ int *ellipsoid = atom->ellipsoid;
int *mask = atom->mask;
int nlocal = atom->nlocal;
for (int i = 0; i < nlocal; i++) {
- if (mask[i] & groupbit) buf[n] = shape[i][2];
+ if ((mask[i] & groupbit) && ellipsoid[i] >= 0)
+ buf[n] = bonus[ellipsoid[i]].shape[2];
else buf[n] = 0.0;
n += nvalues;
}
@@ -1018,12 +1025,14 @@ void ComputePropertyAtom::pack_shapez(int n)
void ComputePropertyAtom::pack_quatw(int n)
{
- double **quat = atom->quat;
+ AtomVecEllipsoid::Bonus *bonus = avec->bonus;
+ int *ellipsoid = atom->ellipsoid;
int *mask = atom->mask;
int nlocal = atom->nlocal;
for (int i = 0; i < nlocal; i++) {
- if (mask[i] & groupbit) buf[n] = quat[i][0];
+ if ((mask[i] & groupbit) && ellipsoid[i] >= 0)
+ buf[n] = bonus[ellipsoid[i]].quat[0];
else buf[n] = 0.0;
n += nvalues;
}
@@ -1033,12 +1042,14 @@ void ComputePropertyAtom::pack_quatw(int n)
void ComputePropertyAtom::pack_quati(int n)
{
- double **quat = atom->quat;
+ AtomVecEllipsoid::Bonus *bonus = avec->bonus;
+ int *ellipsoid = atom->ellipsoid;
int *mask = atom->mask;
int nlocal = atom->nlocal;
for (int i = 0; i < nlocal; i++) {
- if (mask[i] & groupbit) buf[n] = quat[i][1];
+ if ((mask[i] & groupbit) && ellipsoid[i] >= 0)
+ buf[n] = bonus[ellipsoid[i]].quat[1];
else buf[n] = 0.0;
n += nvalues;
}
@@ -1048,12 +1059,14 @@ void ComputePropertyAtom::pack_quati(int n)
void ComputePropertyAtom::pack_quatj(int n)
{
- double **quat = atom->quat;
+ AtomVecEllipsoid::Bonus *bonus = avec->bonus;
+ int *ellipsoid = atom->ellipsoid;
int *mask = atom->mask;
int nlocal = atom->nlocal;
for (int i = 0; i < nlocal; i++) {
- if (mask[i] & groupbit) buf[n] = quat[i][2];
+ if ((mask[i] & groupbit) && ellipsoid[i] >= 0)
+ buf[n] = bonus[ellipsoid[i]].quat[2];
else buf[n] = 0.0;
n += nvalues;
}
@@ -1063,12 +1076,14 @@ void ComputePropertyAtom::pack_quatj(int n)
void ComputePropertyAtom::pack_quatk(int n)
{
- double **quat = atom->quat;
+ AtomVecEllipsoid::Bonus *bonus = avec->bonus;
+ int *ellipsoid = atom->ellipsoid;
int *mask = atom->mask;
int nlocal = atom->nlocal;
for (int i = 0; i < nlocal; i++) {
- if (mask[i] & groupbit) buf[n] = quat[i][3];
+ if ((mask[i] & groupbit) && ellipsoid[i] >= 0)
+ buf[n] = bonus[ellipsoid[i]].quat[3];
else buf[n] = 0.0;
n += nvalues;
}
diff --git a/src/compute_property_atom.h b/src/compute_property_atom.h
index 0c3bd21b32c52a3422f5069339feef9ed386eb22..dcdcd5f5a5e6828cf91c7f6576b13686f3a169da 100644
--- a/src/compute_property_atom.h
+++ b/src/compute_property_atom.h
@@ -38,6 +38,7 @@ class ComputePropertyAtom : public Compute {
double *vector;
double **array;
double *buf;
+ class AtomVecEllipsoid *avec;
typedef void (ComputePropertyAtom::*FnPtrPack)(int);
FnPtrPack *pack_choice; // ptrs to pack functions
diff --git a/src/delete_atoms.cpp b/src/delete_atoms.cpp
index 56879053d609b82b426a8b5a137b3e668d6d126b..dba5ea5ce499a520c1a4be9b7f7af956834fb946 100644
--- a/src/delete_atoms.cpp
+++ b/src/delete_atoms.cpp
@@ -69,7 +69,7 @@ void DeleteAtoms::command(int narg, char **arg)
int i = 0;
while (i < nlocal) {
if (dlist[i]) {
- avec->copy(nlocal-1,i);
+ avec->copy(nlocal-1,i,1);
dlist[i] = dlist[nlocal-1];
nlocal--;
} else i++;
diff --git a/src/dump_cfg.cpp b/src/dump_cfg.cpp
index 4580e4ad6da6b571b5be5cecd03df902f2de3fac..471b5d3bb00a517195c9fa319b503078e804cd16 100755
--- a/src/dump_cfg.cpp
+++ b/src/dump_cfg.cpp
@@ -191,7 +191,7 @@ void DumpCFG::write_header(bigint n)
// scale box dimension to sc lattice for C with sigma = 1.44 Angstroms
double scale;
- if (atom->style_match("peri")) {
+ if (atom->peri_flag) {
int nlocal = atom->nlocal;
double vone = 0.0;
for (int i = 0; i < nlocal; i++) vone += atom->vfrac[i];
diff --git a/src/dump_custom.cpp b/src/dump_custom.cpp
index 3699a518056d9fa63f79de25436bf25b0e1bb6e3..823c5c0dbf01a686136f1f7af1b97f9211790f32 100644
--- a/src/dump_custom.cpp
+++ b/src/dump_custom.cpp
@@ -38,8 +38,7 @@ enum{ID,MOL,TYPE,MASS,
X,Y,Z,XS,YS,ZS,XSTRI,YSTRI,ZSTRI,XU,YU,ZU,XUTRI,YUTRI,ZUTRI,IX,IY,IZ,
VX,VY,VZ,FX,FY,FZ,
Q,MUX,MUY,MUZ,MU,RADIUS,OMEGAX,OMEGAY,OMEGAZ,ANGMOMX,ANGMOMY,ANGMOMZ,
- SHAPEX,SHAPEY,SHAPEZ,
- QUATW,QUATI,QUATJ,QUATK,TQX,TQY,TQZ,SPIN,ERADIUS,ERVEL,ERFORCE,
+ TQX,TQY,TQZ,SPIN,ERADIUS,ERVEL,ERFORCE,
COMPUTE,FIX,VARIABLE};
enum{LT,LE,GT,GE,EQ,NEQ};
enum{INT,DOUBLE};
@@ -663,42 +662,6 @@ int DumpCustom::count()
error->all("Threshhold for an atom property that isn't allocated");
ptr = &atom->angmom[0][2];
nstride = 3;
-
- } else if (thresh_array[ithresh] == SHAPEX) {
- if (!atom->shape_flag)
- error->all("Threshhold for an atom property that isn't allocated");
- ptr = &atom->shape[0][0];
- nstride = 3;
- } else if (thresh_array[ithresh] == SHAPEY) {
- if (!atom->shape_flag)
- error->all("Threshhold for an atom property that isn't allocated");
- ptr = &atom->shape[0][1];
- nstride = 3;
- } else if (thresh_array[ithresh] == SHAPEZ) {
- if (!atom->shape_flag)
- error->all("Threshhold for an atom property that isn't allocated");
- ptr = &atom->shape[0][2];
- nstride = 3;
- } else if (thresh_array[ithresh] == QUATW) {
- if (!atom->quat_flag)
- error->all("Threshhold for an atom property that isn't allocated");
- ptr = &atom->quat[0][0];
- nstride = 4;
- } else if (thresh_array[ithresh] == QUATI) {
- if (!atom->quat_flag)
- error->all("Threshhold for an atom property that isn't allocated");
- ptr = &atom->quat[0][1];
- nstride = 4;
- } else if (thresh_array[ithresh] == QUATJ) {
- if (!atom->quat_flag)
- error->all("Threshhold for an atom property that isn't allocated");
- ptr = &atom->quat[0][2];
- nstride = 4;
- } else if (thresh_array[ithresh] == QUATK) {
- if (!atom->quat_flag)
- error->all("Threshhold for an atom property that isn't allocated");
- ptr = &atom->quat[0][3];
- nstride = 4;
} else if (thresh_array[ithresh] == TQX) {
if (!atom->torque_flag)
error->all("Threshhold for an atom property that isn't allocated");
@@ -1006,42 +969,6 @@ void DumpCustom::parse_fields(int narg, char **arg)
error->all("Dumping an atom property that isn't allocated");
pack_choice[i] = &DumpCustom::pack_angmomz;
vtype[i] = DOUBLE;
-
- } else if (strcmp(arg[iarg],"shapex") == 0) {
- if (!atom->shape_flag)
- error->all("Dumping an atom property that isn't allocated");
- pack_choice[i] = &DumpCustom::pack_shapex;
- vtype[i] = DOUBLE;
- } else if (strcmp(arg[iarg],"shapey") == 0) {
- if (!atom->shape_flag)
- error->all("Dumping an atom property that isn't allocated");
- pack_choice[i] = &DumpCustom::pack_shapey;
- vtype[i] = DOUBLE;
- } else if (strcmp(arg[iarg],"shapez") == 0) {
- if (!atom->shape_flag)
- error->all("Dumping an atom property that isn't allocated");
- pack_choice[i] = &DumpCustom::pack_shapez;
- vtype[i] = DOUBLE;
- } else if (strcmp(arg[iarg],"quatw") == 0) {
- if (!atom->quat_flag)
- error->all("Dumping an atom property that isn't allocated");
- pack_choice[i] = &DumpCustom::pack_quatw;
- vtype[i] = DOUBLE;
- } else if (strcmp(arg[iarg],"quati") == 0) {
- if (!atom->quat_flag)
- error->all("Dumping an atom property that isn't allocated");
- pack_choice[i] = &DumpCustom::pack_quati;
- vtype[i] = DOUBLE;
- } else if (strcmp(arg[iarg],"quatj") == 0) {
- if (!atom->quat_flag)
- error->all("Dumping an atom property that isn't allocated");
- pack_choice[i] = &DumpCustom::pack_quatj;
- vtype[i] = DOUBLE;
- } else if (strcmp(arg[iarg],"quatk") == 0) {
- if (!atom->quat_flag)
- error->all("Dumping an atom property that isn't allocated");
- pack_choice[i] = &DumpCustom::pack_quatk;
- vtype[i] = DOUBLE;
} else if (strcmp(arg[iarg],"tqx") == 0) {
if (!atom->torque_flag)
error->all("Dumping an atom property that isn't allocated");
@@ -1350,14 +1277,6 @@ int DumpCustom::modify_param(int narg, char **arg)
else if (strcmp(arg[1],"angmomx") == 0) thresh_array[nthresh] = ANGMOMX;
else if (strcmp(arg[1],"angmomy") == 0) thresh_array[nthresh] = ANGMOMY;
else if (strcmp(arg[1],"angmomz") == 0) thresh_array[nthresh] = ANGMOMZ;
-
- else if (strcmp(arg[1],"shapex") == 0) thresh_array[nthresh] = SHAPEX;
- else if (strcmp(arg[1],"shapey") == 0) thresh_array[nthresh] = SHAPEY;
- else if (strcmp(arg[1],"shapez") == 0) thresh_array[nthresh] = SHAPEZ;
- else if (strcmp(arg[1],"quatw") == 0) thresh_array[nthresh] = QUATW;
- else if (strcmp(arg[1],"quati") == 0) thresh_array[nthresh] = QUATI;
- else if (strcmp(arg[1],"quatj") == 0) thresh_array[nthresh] = QUATJ;
- else if (strcmp(arg[1],"quatk") == 0) thresh_array[nthresh] = QUATK;
else if (strcmp(arg[1],"tqx") == 0) thresh_array[nthresh] = TQX;
else if (strcmp(arg[1],"tqy") == 0) thresh_array[nthresh] = TQY;
else if (strcmp(arg[1],"tqz") == 0) thresh_array[nthresh] = TQZ;
@@ -2196,104 +2115,6 @@ void DumpCustom::pack_angmomz(int n)
/* ---------------------------------------------------------------------- */
-void DumpCustom::pack_shapex(int n)
-{
- double **shape = atom->shape;
- int nlocal = atom->nlocal;
-
- for (int i = 0; i < nlocal; i++)
- if (choose[i]) {
- buf[n] = shape[i][0];
- n += size_one;
- }
-}
-
-/* ---------------------------------------------------------------------- */
-
-void DumpCustom::pack_shapey(int n)
-{
- double **shape = atom->shape;
- int nlocal = atom->nlocal;
-
- for (int i = 0; i < nlocal; i++)
- if (choose[i]) {
- buf[n] = shape[i][1];
- n += size_one;
- }
-}
-
-/* ---------------------------------------------------------------------- */
-
-void DumpCustom::pack_shapez(int n)
-{
- double **shape = atom->shape;
- int nlocal = atom->nlocal;
-
- for (int i = 0; i < nlocal; i++)
- if (choose[i]) {
- buf[n] = shape[i][2];
- n += size_one;
- }
-}
-
-/* ---------------------------------------------------------------------- */
-
-void DumpCustom::pack_quatw(int n)
-{
- double **quat = atom->quat;
- int nlocal = atom->nlocal;
-
- for (int i = 0; i < nlocal; i++)
- if (choose[i]) {
- buf[n] = quat[i][0];
- n += size_one;
- }
-}
-
-/* ---------------------------------------------------------------------- */
-
-void DumpCustom::pack_quati(int n)
-{
- double **quat = atom->quat;
- int nlocal = atom->nlocal;
-
- for (int i = 0; i < nlocal; i++)
- if (choose[i]) {
- buf[n] = quat[i][1];
- n += size_one;
- }
-}
-
-/* ---------------------------------------------------------------------- */
-
-void DumpCustom::pack_quatj(int n)
-{
- double **quat = atom->quat;
- int nlocal = atom->nlocal;
-
- for (int i = 0; i < nlocal; i++)
- if (choose[i]) {
- buf[n] = quat[i][2];
- n += size_one;
- }
-}
-
-/* ---------------------------------------------------------------------- */
-
-void DumpCustom::pack_quatk(int n)
-{
- double **quat = atom->quat;
- int nlocal = atom->nlocal;
-
- for (int i = 0; i < nlocal; i++)
- if (choose[i]) {
- buf[n] = quat[i][3];
- n += size_one;
- }
-}
-
-/* ---------------------------------------------------------------------- */
-
void DumpCustom::pack_tqx(int n)
{
double **torque = atom->torque;
diff --git a/src/dump_custom.h b/src/dump_custom.h
index 3c9bbd2bd8cd065c452de70a1fe42421a057c4d5..c5c867de5a5e5f0cde2963042c521b964308afe5 100644
--- a/src/dump_custom.h
+++ b/src/dump_custom.h
@@ -142,13 +142,6 @@ class DumpCustom : public Dump {
void pack_angmomx(int);
void pack_angmomy(int);
void pack_angmomz(int);
- void pack_shapex(int);
- void pack_shapey(int);
- void pack_shapez(int);
- void pack_quatw(int);
- void pack_quati(int);
- void pack_quatj(int);
- void pack_quatk(int);
void pack_tqx(int);
void pack_tqy(int);
void pack_tqz(int);
diff --git a/src/fix_evaporate.cpp b/src/fix_evaporate.cpp
index d4a963ecf0c8d60c79a03856f0c0be9708269742..c5db3dca4c0199559823d414eb45629e366001e0 100644
--- a/src/fix_evaporate.cpp
+++ b/src/fix_evaporate.cpp
@@ -287,7 +287,7 @@ void FixEvaporate::pre_exchange()
for (i = nlocal-1; i >= 0; i--) {
if (mark[i]) {
- avec->copy(atom->nlocal-1,i);
+ avec->copy(atom->nlocal-1,i,1);
atom->nlocal--;
}
}
diff --git a/src/fix_move.cpp b/src/fix_move.cpp
index 7d18f6145cfd36fdea9de688330b67e58b15523e..1aa844923455d9a98d4ac1540c2ed690a5886238 100644
--- a/src/fix_move.cpp
+++ b/src/fix_move.cpp
@@ -181,7 +181,7 @@ FixMove::FixMove(LAMMPS *lmp, int narg, char **arg) :
if (atom->angmom_flag && comm->me == 0)
error->warning("Fix move does not update angular momentum");
- if (atom->quat_flag && comm->me == 0)
+ if (atom->ellipsoid_flag && comm->me == 0)
error->warning("Fix move does not update quaternions");
// setup scaling and apply scaling factors to velocity & amplitude
diff --git a/src/fix_rigid.cpp b/src/fix_rigid.cpp
index bfea6f4a49b834336c362bc0240a1329dc8f2712..647162b7eb2dfc53d445b6d53220745819f4278f 100644
--- a/src/fix_rigid.cpp
+++ b/src/fix_rigid.cpp
@@ -18,7 +18,7 @@
#include "fix_rigid.h"
#include "math_extra.h"
#include "atom.h"
-#include "atom_vec.h"
+#include "atom_vec_ellipsoid.h"
#include "domain.h"
#include "update.h"
#include "respa.h"
@@ -361,6 +361,10 @@ FixRigid::FixRigid(LAMMPS *lmp, int narg, char **arg) :
ANGMOM = 64;
TORQUE = 128;
+ // atom style pointers to particles that store extra info
+
+ avec_ellipsoid = (AtomVecEllipsoid *) atom->style_match("ellipsoid");
+
// print statistics
int nsum = 0;
@@ -466,20 +470,21 @@ void FixRigid::init()
extended = dorientflag = qorientflag = 0;
- double **shape = atom->shape;
+ AtomVecEllipsoid::Bonus *ebonus = avec_ellipsoid->bonus;
double **mu = atom->mu;
double *radius = atom->radius;
double *rmass = atom->rmass;
double *mass = atom->mass;
+ int *ellipsoid = atom->ellipsoid;
int *type = atom->type;
int nlocal = atom->nlocal;
- if (atom->radius_flag || atom->shape_flag) {
+ if (atom->radius_flag || atom->ellipsoid_flag || atom->mu_flag) {
int flag = 0;
for (i = 0; i < nlocal; i++) {
if (body[i] < 0) continue;
if (radius && radius[i] > 0.0) flag = 1;
- if (shape && shape[i][0] > 0.0) flag = 1;
+ if (ellipsoid && ellipsoid[i] >= 0) flag = 1;
if (mu && mu[i][3] > 0.0) flag = 1;
}
@@ -492,7 +497,7 @@ void FixRigid::init()
if (extended) {
if (atom->mu_flag) dorientflag = 1;
- if (atom->quat_flag) qorientflag = 1;
+ if (atom->ellipsoid_flag) qorientflag = 1;
grow_arrays(atom->nmax);
for (i = 0; i < nlocal; i++) {
@@ -505,7 +510,7 @@ void FixRigid::init()
eflags[i] |= INERTIA_SPHERE;
eflags[i] |= OMEGA;
eflags[i] |= TORQUE;
- } else if (shape && shape[i][0] > 0.0) {
+ } else if (ellipsoid && ellipsoid[i] >= 0) {
eflags[i] |= INERTIA_ELLIPSOID;
eflags[i] |= ORIENT_QUAT;
eflags[i] |= ANGMOM;
@@ -629,9 +634,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 *radius = atom->radius;
- double **shape = atom->shape;
+ double *shape,*quatatom;
for (i = 0; i < nlocal; i++) {
if (body[i] < 0) continue;
@@ -647,14 +650,13 @@ void FixRigid::init()
sum[ibody][2] += 0.4 * massone * radius[i]*radius[i];
}
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[i][1]*shape[i][1]+shape[itype][2]*shape[i][2]);
- idiag[1] = 0.2*massone *
- (shape[i][0]*shape[i][0]+shape[i][2]*shape[i][2]);
- idiag[2] = 0.2*massone *
- (shape[i][0]*shape[i][0]+shape[i][1]*shape[i][1]);
+ 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];
@@ -744,6 +746,7 @@ void FixRigid::init()
// for extended particles, set their orientation wrt to rigid body
double qc[4];
+ double *quatatom;
for (i = 0; i < nlocal; i++) {
if (body[i] < 0) {
@@ -779,9 +782,6 @@ void FixRigid::init()
dz*ez_space[ibody][2];
if (extended) {
- double **mu = atom->mu;
- int *type = atom->type;
-
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];
@@ -794,8 +794,9 @@ void FixRigid::init()
dorient[i][0] = dorient[i][1] = dorient[i][2] = 0.0;
if (eflags[i] & ORIENT_QUAT) {
+ quatatom = ebonus[ellipsoid[i]].quat;
qconjugate(quat[ibody],qc);
- quatquat(qc,atom->quat[i],qorient[i]);
+ quatquat(qc,quatatom,qorient[i]);
qnormalize(qorient[i]);
} else if (qorientflag)
qorient[i][0] = qorient[i][1] = qorient[i][2] = qorient[i][3] = 0.0;
@@ -831,9 +832,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 *radius = atom->radius;
- double **shape = atom->shape;
+ double *shape;
for (i = 0; i < nlocal; i++) {
if (body[i] < 0) continue;
@@ -849,14 +848,12 @@ void FixRigid::init()
sum[ibody][2] += 0.4 * massone * radius[i]*radius[i];
}
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[i][1]*shape[i][1]+shape[i][2]*shape[i][2]);
- idiag[1] = 0.2*massone *
- (shape[i][0]*shape[i][0]+shape[i][2]*shape[i][2]);
- idiag[2] = 0.2*massone *
- (shape[i][0]*shape[i][0]+shape[i][1]*shape[i][1]);
+ 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];
@@ -1934,10 +1931,13 @@ void FixRigid::set_xv()
// set orientation, omega, angmom of each extended particle
if (extended) {
+ double *shape,*quatatom;
+
+ AtomVecEllipsoid::Bonus *ebonus = avec_ellipsoid->bonus;
double **omega_one = atom->omega;
double **angmom_one = atom->angmom;
- double **shape = atom->shape;
double **mu = atom->mu;
+ int *ellipsoid = atom->ellipsoid;
for (int i = 0; i < nlocal; i++) {
if (body[i] < 0) continue;
@@ -1949,8 +1949,9 @@ void FixRigid::set_xv()
MathExtra::snormalize3(mu[i][3],mu[i],mu[i]);
}
if (eflags[i] & ORIENT_QUAT) {
- quatquat(quat[ibody],qorient[i],atom->quat[i]);
- qnormalize(atom->quat[i]);
+ quatatom = ebonus[ellipsoid[i]].quat;
+ quatquat(quat[ibody],qorient[i],quatatom);
+ qnormalize(quatatom);
}
if (eflags[i] & OMEGA) {
omega_one[i][0] = omega[ibody][0];
@@ -1958,13 +1959,12 @@ void FixRigid::set_xv()
omega_one[i][2] = omega[ibody][2];
}
if (eflags[i] & ANGMOM) {
- ione[0] = 0.2*rmass[i] *
- (shape[i][1]*shape[i][1] + shape[i][2]*shape[i][2]);
- ione[1] = 0.2*rmass[i] *
- (shape[i][0]*shape[i][0] + shape[i][2]*shape[i][2]);
- ione[2] = 0.2*rmass[i] *
- (shape[i][0]*shape[i][0] + shape[i][1]*shape[i][1]);
- exyz_from_q(atom->quat[i],exone,eyone,ezone);
+ shape = ebonus[ellipsoid[i]].shape;
+ quatatom = ebonus[ellipsoid[i]].quat;
+ 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]);
}
}
@@ -2073,9 +2073,12 @@ void FixRigid::set_v()
// set omega, angmom of each extended particle
if (extended) {
+ double *shape,*quatatom;
+
+ AtomVecEllipsoid::Bonus *ebonus = avec_ellipsoid->bonus;
double **omega_one = atom->omega;
double **angmom_one = atom->angmom;
- double **shape = atom->shape;
+ int *ellipsoid = atom->ellipsoid;
for (int i = 0; i < nlocal; i++) {
if (body[i] < 0) continue;
@@ -2087,13 +2090,12 @@ void FixRigid::set_v()
omega_one[i][2] = omega[ibody][2];
}
if (eflags[i] & ANGMOM) {
- ione[0] = 0.2*rmass[i] *
- (shape[i][1]*shape[i][1] + shape[i][2]*shape[i][2]);
- ione[1] = 0.2*rmass[i] *
- (shape[i][0]*shape[i][0] + shape[i][2]*shape[i][2]);
- ione[2] = 0.2*rmass[i] *
- (shape[i][0]*shape[i][0] + shape[i][1]*shape[i][1]);
- exyz_from_q(atom->quat[i],exone,eyone,ezone);
+ shape = ebonus[ellipsoid[i]].shape;
+ quatatom = ebonus[ellipsoid[i]].quat;
+ 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]);
}
}
diff --git a/src/fix_rigid.h b/src/fix_rigid.h
index 962a83d7bbc35488db2051a05e52367887491f66..6428d7c32794d401e21d882112328d21f399f58e 100644
--- a/src/fix_rigid.h
+++ b/src/fix_rigid.h
@@ -95,6 +95,8 @@ class FixRigid : public Fix {
double p_period,p_freq;
int p_chain;
+ class AtomVecEllipsoid *avec_ellipsoid;
+
// bitmasks for eflags
int INERTIA_POINT,INERTIA_SPHERE,INERTIA_ELLIPSOID;
int ORIENT_DIPOLE,ORIENT_QUAT;
diff --git a/src/fix_store_state.cpp b/src/fix_store_state.cpp
index 5dc2746839b613ecce0aa7e8f302a5aac1c6131d..bbae452f96abf0066dde93e6e251f19823e9e681 100644
--- a/src/fix_store_state.cpp
+++ b/src/fix_store_state.cpp
@@ -167,23 +167,6 @@ FixStoreState::FixStoreState(LAMMPS *lmp, int narg, char **arg) :
if (!atom->angmom_flag)
error->all("Fix store/state for atom property that isn't allocated");
pack_choice[nvalues++] = &FixStoreState::pack_angmomz;
-
- } else if (strcmp(arg[iarg],"quatw") == 0) {
- if (!atom->quat_flag)
- error->all("Fix store/state for atom property that isn't allocated");
- pack_choice[nvalues++] = &FixStoreState::pack_quatw;
- } else if (strcmp(arg[iarg],"quati") == 0) {
- if (!atom->quat_flag)
- error->all("Fix store/state for atom property that isn't allocated");
- pack_choice[nvalues++] = &FixStoreState::pack_quati;
- } else if (strcmp(arg[iarg],"quatj") == 0) {
- if (!atom->quat_flag)
- error->all("Fix store/state for atom property that isn't allocated");
- pack_choice[nvalues++] = &FixStoreState::pack_quatj;
- } else if (strcmp(arg[iarg],"quatk") == 0) {
- if (!atom->quat_flag)
- error->all("Fix store/state for atom property that isn't allocated");
- pack_choice[nvalues++] = &FixStoreState::pack_quatk;
} else if (strcmp(arg[iarg],"tqx") == 0) {
if (!atom->torque_flag)
error->all("Fix store/state for atom property that isn't allocated");
@@ -1231,66 +1214,6 @@ void FixStoreState::pack_angmomz(int n)
/* ---------------------------------------------------------------------- */
-void FixStoreState::pack_quatw(int n)
-{
- double **quat = atom->quat;
- int *mask = atom->mask;
- int nlocal = atom->nlocal;
-
- for (int i = 0; i < nlocal; i++) {
- if (mask[i] & groupbit) vbuf[n] = quat[i][0];
- else vbuf[n] = 0.0;
- n += nvalues;
- }
-}
-
-/* ---------------------------------------------------------------------- */
-
-void FixStoreState::pack_quati(int n)
-{
- double **quat = atom->quat;
- int *mask = atom->mask;
- int nlocal = atom->nlocal;
-
- for (int i = 0; i < nlocal; i++) {
- if (mask[i] & groupbit) vbuf[n] = quat[i][1];
- else vbuf[n] = 0.0;
- n += nvalues;
- }
-}
-
-/* ---------------------------------------------------------------------- */
-
-void FixStoreState::pack_quatj(int n)
-{
- double **quat = atom->quat;
- int *mask = atom->mask;
- int nlocal = atom->nlocal;
-
- for (int i = 0; i < nlocal; i++) {
- if (mask[i] & groupbit) vbuf[n] = quat[i][2];
- else vbuf[n] = 0.0;
- n += nvalues;
- }
-}
-
-/* ---------------------------------------------------------------------- */
-
-void FixStoreState::pack_quatk(int n)
-{
- double **quat = atom->quat;
- int *mask = atom->mask;
- int nlocal = atom->nlocal;
-
- for (int i = 0; i < nlocal; i++) {
- if (mask[i] & groupbit) vbuf[n] = quat[i][3];
- else vbuf[n] = 0.0;
- n += nvalues;
- }
-}
-
-/* ---------------------------------------------------------------------- */
-
void FixStoreState::pack_tqx(int n)
{
double **torque = atom->torque;
diff --git a/src/fix_store_state.h b/src/fix_store_state.h
index 5ef97caf98f9665a945e48a71993cb0090d64a49..11d83366abfd37a5f8c6eff6c49285f42cd6fa8e 100644
--- a/src/fix_store_state.h
+++ b/src/fix_store_state.h
@@ -100,10 +100,6 @@ class FixStoreState : public Fix {
void pack_angmomx(int);
void pack_angmomy(int);
void pack_angmomz(int);
- void pack_quatw(int);
- void pack_quati(int);
- void pack_quatj(int);
- void pack_quatk(int);
void pack_tqx(int);
void pack_tqy(int);
void pack_tqz(int);
diff --git a/src/force.cpp b/src/force.cpp
index 23a4707d5ae19a2464e94b4a1de8adfd9aa25145..8068717be5077a768c0bfeb1498c872f1107545f 100644
--- a/src/force.cpp
+++ b/src/force.cpp
@@ -147,7 +147,7 @@ Pair *Force::new_pair(const char *style)
}
/* ----------------------------------------------------------------------
- return ptr to current pair class or hybrid sub-class
+ return ptr to Pair class if matches word or to matching hybrid sub-class
if exact, then style name must be exact match to word
if not exact, style name must contain word
return NULL if no match
diff --git a/src/irregular.cpp b/src/irregular.cpp
index 463ee5fc0e415634077f5d0f4dd74fbac732ac65..d9abe11b88704570dd8acc1087f710a74eb4be12 100644
--- a/src/irregular.cpp
+++ b/src/irregular.cpp
@@ -78,12 +78,13 @@ Irregular::~Irregular()
void Irregular::migrate_atoms()
{
// clear global->local map since atoms move to new procs
- // zero out ghosts so map_set() at end will operate only on local atoms
- // exchange() doesn't need to zero ghosts b/c borders()
- // is called right after and it zeroes ghosts and calls map_set()
+ // clear old ghosts so map_set() at end will operate only on local atoms
+ // exchange() doesn't need to clear ghosts b/c borders()
+ // is called right after and it clears ghosts and calls map_set()
if (map_style) atom->map_clear();
atom->nghost = 0;
+ atom->avec->clear_bonus();
// subbox bounds for orthogonal or triclinic
@@ -123,7 +124,7 @@ void Irregular::migrate_atoms()
sizes[nsendatom] = avec->pack_exchange(i,&buf_send[nsend]);
nsend += sizes[nsendatom];
nsendatom++;
- avec->copy(nlocal-1,i);
+ avec->copy(nlocal-1,i,1);
nlocal--;
} else i++;
} else i++;
diff --git a/src/read_data.cpp b/src/read_data.cpp
index 301498cf6289e95fdd9a65c621a83b9540027f4a..69fd1ef77c2b9041b3a3679922e2b330ba97615f 100644
--- a/src/read_data.cpp
+++ b/src/read_data.cpp
@@ -19,6 +19,7 @@
#include "read_data.h"
#include "atom.h"
#include "atom_vec.h"
+#include "atom_vec_ellipsoid.h"
#include "comm.h"
#include "update.h"
#include "force.h"
@@ -39,7 +40,8 @@ using namespace LAMMPS_NS;
#define CHUNK 1024
#define DELTA 4 // must be 2 or larger
-#define NSECTIONS 22 // change when add to header::section_keywords
+ // customize for new sections
+#define NSECTIONS 23 // change when add to header::section_keywords
#define MIN(a,b) ((a) < (b) ? (a) : (b))
#define MAX(a,b) ((a) > (b) ? (a) : (b))
@@ -54,6 +56,12 @@ ReadData::ReadData(LAMMPS *lmp) : Pointers(lmp)
buffer = new char[CHUNK*MAXLINE];
narg = maxarg = 0;
arg = NULL;
+
+ // customize for new sections
+ // pointers to atom styles that store extra info
+
+ nellipsoids = 0;
+ avec_ellipsoid = (AtomVecEllipsoid *) atom->style_match("ellipsoid");
}
/* ---------------------------------------------------------------------- */
@@ -128,8 +136,8 @@ void ReadData::command(int narg, char **arg)
comm->set_procs();
domain->set_local_box();
+ // customize for new sections
// read rest of file in free format
- // if add a section keyword, add to header::section_keywords and NSECTIONS
int atomflag = 0;
@@ -140,6 +148,12 @@ void ReadData::command(int narg, char **arg)
} else if (strcmp(keyword,"Velocities") == 0) {
if (atomflag == 0) error->all("Must read Atoms before Velocities");
velocities();
+ } else if (strcmp(keyword,"Ellipsoids") == 0) {
+ if (!avec_ellipsoid)
+ error->all("Invalid data file section: Ellipsoids");
+ if (atomflag == 0) error->all("Must read Atoms before Ellipsoids");
+ ellipsoids();
+
} else if (strcmp(keyword,"Bonds") == 0) {
if (atom->avec->bonds_allow == 0)
error->all("Invalid data file section: Bonds");
@@ -287,14 +301,17 @@ void ReadData::header(int flag)
int n;
char *ptr;
+ // customize for new sections
+
char *section_keywords[NSECTIONS] =
- {"Atoms","Velocities","Bonds","Angles","Dihedrals","Impropers","Masses",
- "Pair Coeffs","Bond Coeffs","Angle Coeffs",
+ {"Atoms","Velocities","Ellipsoids",
+ "Bonds","Angles","Dihedrals","Impropers",
+ "Masses","Pair Coeffs","Bond Coeffs","Angle Coeffs",
"Dihedral Coeffs","Improper Coeffs",
"BondBond Coeffs","BondAngle Coeffs","MiddleBondTorsion Coeffs",
"EndBondTorsion Coeffs","AngleTorsion Coeffs",
"AngleAngleTorsion Coeffs","BondBond13 Coeffs","AngleAngle Coeffs"};
-
+
// skip 1st line of file
if (me == 0) {
@@ -302,6 +319,8 @@ void ReadData::header(int flag)
if (eof == NULL) error->one("Unexpected end of data file");
}
+ // customize for new header lines
+
while (1) {
// read a line and bcast length if flag is set
@@ -351,6 +370,12 @@ void ReadData::header(int flag)
else if (strstr(line,"extra bond per atom"))
sscanf(line,"%d",&atom->extra_bond_per_atom);
+ else if (strstr(line,"ellipsoids")) {
+ if (!avec_ellipsoid)
+ error->all("No ellipsoids allowed with this atom style");
+ sscanf(line,BIGINT_FORMAT,&nellipsoids);
+ }
+
else if (strstr(line,"xlo xhi"))
sscanf(line,"%lg %lg",&domain->boxlo[0],&domain->boxhi[0]);
else if (strstr(line,"ylo yhi"))
@@ -543,6 +568,57 @@ void ReadData::velocities()
}
}
+/* ----------------------------------------------------------------------
+ read all ellipsoids
+ to find atoms, must build atom map if not a molecular system
+------------------------------------------------------------------------- */
+
+void ReadData::ellipsoids()
+{
+ int i,m,nchunk;
+
+ int mapflag = 0;
+ if (atom->map_style == 0) {
+ mapflag = 1;
+ atom->map_style = 1;
+ atom->map_init();
+ atom->map_set();
+ }
+
+ bigint nread = 0;
+ bigint natoms = nellipsoids;
+
+ while (nread < natoms) {
+ if (natoms-nread > CHUNK) nchunk = CHUNK;
+ else nchunk = natoms-nread;
+ if (me == 0) {
+ char *eof;
+ m = 0;
+ for (i = 0; i < nchunk; i++) {
+ eof = fgets(&buffer[m],MAXLINE,fp);
+ if (eof == NULL) error->one("Unexpected end of data file");
+ m += strlen(&buffer[m]);
+ }
+ buffer[m++] = '\n';
+ }
+ MPI_Bcast(&m,1,MPI_INT,0,world);
+ MPI_Bcast(buffer,m,MPI_CHAR,0,world);
+
+ atom->data_bonus(nchunk,buffer,avec_ellipsoid);
+ nread += nchunk;
+ }
+
+ if (mapflag) {
+ atom->map_delete();
+ atom->map_style = 0;
+ }
+
+ if (me == 0) {
+ if (screen) fprintf(screen," " BIGINT_FORMAT " ellipsoids\n",natoms);
+ if (logfile) fprintf(logfile," " BIGINT_FORMAT " ellipsoids\n",natoms);
+ }
+}
+
/* ---------------------------------------------------------------------- */
void ReadData::bonds()
@@ -933,6 +1009,7 @@ void ReadData::scan(int &bond_per_atom, int &angle_per_atom,
int natoms = static_cast<int> (atom->natoms);
bond_per_atom = angle_per_atom = dihedral_per_atom = improper_per_atom = 0;
+ // customize for new sections
// allocate topology counting vector
// initially, array length = 1 to natoms
// will grow via reallocate() if atom IDs > natoms
@@ -947,32 +1024,33 @@ void ReadData::scan(int &bond_per_atom, int &angle_per_atom,
else if (strcmp(keyword,"Atoms") == 0) skip_lines(natoms);
else if (strcmp(keyword,"Velocities") == 0) skip_lines(natoms);
- else if (strcmp(keyword,"Pair Coeffs") == 0) {
+ else if (strcmp(keyword,"Ellipsoids") == 0) {
+ if (!avec_ellipsoid)
+ error->all("Invalid data file section: Ellipsoids");
+ skip_lines(nellipsoids);
+
+ } else if (strcmp(keyword,"Pair Coeffs") == 0) {
if (force->pair == NULL)
error->all("Must define pair_style before Pair Coeffs");
skip_lines(atom->ntypes);
-
} else if (strcmp(keyword,"Bond Coeffs") == 0) {
if (atom->avec->bonds_allow == 0)
error->all("Invalid data file section: Bond Coeffs");
if (force->bond == NULL)
error->all("Must define bond_style before Bond Coeffs");
skip_lines(atom->nbondtypes);
-
} else if (strcmp(keyword,"Angle Coeffs") == 0) {
if (atom->avec->angles_allow == 0)
error->all("Invalid data file section: Angle Coeffs");
if (force->angle == NULL)
error->all("Must define angle_style before Angle Coeffs");
skip_lines(atom->nangletypes);
-
} else if (strcmp(keyword,"Dihedral Coeffs") == 0) {
skip_lines(atom->ndihedraltypes);
if (atom->avec->dihedrals_allow == 0)
error->all("Invalid data file section: Dihedral Coeffs");
if (force->dihedral == NULL)
error->all("Must define dihedral_style before Dihedral Coeffs");
-
} else if (strcmp(keyword,"Improper Coeffs") == 0) {
if (atom->avec->impropers_allow == 0)
error->all("Invalid data file section: Improper Coeffs");
@@ -986,49 +1064,42 @@ void ReadData::scan(int &bond_per_atom, int &angle_per_atom,
if (force->angle == NULL)
error->all("Must define angle_style before BondBond Coeffs");
skip_lines(atom->nangletypes);
-
} else if (strcmp(keyword,"BondAngle Coeffs") == 0) {
if (atom->avec->angles_allow == 0)
error->all("Invalid data file section: BondAngle Coeffs");
if (force->angle == NULL)
error->all("Must define angle_style before BondAngle Coeffs");
skip_lines(atom->nangletypes);
-
} else if (strcmp(keyword,"MiddleBondTorsion Coeffs") == 0) {
if (atom->avec->dihedrals_allow == 0)
error->all("Invalid data file section: MiddleBondTorsion Coeffs");
if (force->dihedral == NULL)
error->all("Must define dihedral_style before MiddleBondTorsion Coeffs");
skip_lines(atom->ndihedraltypes);
-
} else if (strcmp(keyword,"EndBondTorsion Coeffs") == 0) {
if (atom->avec->dihedrals_allow == 0)
error->all("Invalid data file section: EndBondTorsion Coeffs");
if (force->dihedral == NULL)
error->all("Must define dihedral_style before EndBondTorsion Coeffs");
skip_lines(atom->ndihedraltypes);
-
} else if (strcmp(keyword,"AngleTorsion Coeffs") == 0) {
if (atom->avec->dihedrals_allow == 0)
error->all("Invalid data file section: AngleTorsion Coeffs");
if (force->dihedral == NULL)
error->all("Must define dihedral_style before AngleTorsion Coeffs");
skip_lines(atom->ndihedraltypes);
-
} else if (strcmp(keyword,"AngleAngleTorsion Coeffs") == 0) {
if (atom->avec->dihedrals_allow == 0)
error->all("Invalid data file section: AngleAngleTorsion Coeffs");
if (force->dihedral == NULL)
error->all("Must define dihedral_style before AngleAngleTorsion Coeffs");
skip_lines(atom->ndihedraltypes);
-
} else if (strcmp(keyword,"BondBond13 Coeffs") == 0) {
if (atom->avec->dihedrals_allow == 0)
error->all("Invalid data file section: BondBond13 Coeffs");
if (force->dihedral == NULL)
error->all("Must define dihedral_style before BondBond13 Coeffs");
skip_lines(atom->ndihedraltypes);
-
} else if (strcmp(keyword,"AngleAngle Coeffs") == 0) {
if (atom->avec->impropers_allow == 0)
error->all("Invalid data file section: AngleAngle Coeffs");
diff --git a/src/read_data.h b/src/read_data.h
index f5c4357d4f15e23de053838f42e5872118674c72..80fd441413b97b222fd4f91e2c26d689aec4d5ae 100644
--- a/src/read_data.h
+++ b/src/read_data.h
@@ -38,6 +38,9 @@ class ReadData : protected Pointers {
int narg,maxarg,compressed;
char **arg;
+ bigint nellipsoids;
+ class AtomVecEllipsoid *avec_ellipsoid;
+
void open(char *);
void scan(int &, int &, int &, int &);
int reallocate(int **, int, int);
@@ -48,14 +51,14 @@ class ReadData : protected Pointers {
void atoms();
void velocities();
+ void ellipsoids();
+
void bonds();
void angles();
void dihedrals();
void impropers();
void mass();
- void dipole();
-
void paircoeffs();
void bondcoeffs();
void anglecoeffs(int);
diff --git a/src/set.cpp b/src/set.cpp
index e1b1a71036dcafdd525b3254223225fe44b0aaec..c21ccdbf1936702ee5bae58ab0dc0363fce92be3 100644
--- a/src/set.cpp
+++ b/src/set.cpp
@@ -18,6 +18,7 @@
#include "set.h"
#include "atom.h"
#include "atom_vec.h"
+#include "atom_vec_ellipsoid.h"
#include "domain.h"
#include "region.h"
#include "group.h"
@@ -142,7 +143,7 @@ void Set::command(int narg, char **arg)
xvalue = atof(arg[iarg+1]);
yvalue = atof(arg[iarg+2]);
zvalue = atof(arg[iarg+3]);
- if (!atom->shape_flag)
+ if (!atom->ellipsoid_flag)
error->all("Cannot set this attribute for this atom style");
if (xvalue < 0.0 || yvalue < 0.0 || zvalue < 0.0)
error->all("Invalid shape in set command");
@@ -177,14 +178,14 @@ void Set::command(int narg, char **arg)
yvalue = atof(arg[iarg+2]);
zvalue = atof(arg[iarg+3]);
wvalue = atof(arg[iarg+4]);
- if (!atom->quat_flag)
+ if (!atom->ellipsoid_flag)
error->all("Cannot set this attribute for this atom style");
set(QUAT);
iarg += 5;
} else if (strcmp(arg[iarg],"quat/random") == 0) {
if (iarg+2 > narg) error->all("Illegal set command");
ivalue = atoi(arg[iarg+1]);
- if (!atom->quat_flag)
+ if (!atom->ellipsoid_flag)
error->all("Cannot set this attribute for this atom style");
if (ivalue <= 0) error->all("Invalid random number seed in set command");
setrandom(QUAT_RANDOM);
@@ -358,6 +359,9 @@ void Set::selection(int n)
void Set::set(int keyword)
{
+ AtomVecEllipsoid *avec_ellipsoid =
+ (AtomVecEllipsoid *) atom->style_match("ellipsoid");
+
selection(atom->nlocal);
int nlocal = atom->nlocal;
@@ -376,25 +380,22 @@ void Set::set(int keyword)
// set shape
- else if (keyword == SHAPE) {
- double **shape = atom->shape;
- shape[i][0] = 0.5 * xvalue;
- shape[i][1] = 0.5 * yvalue;
- shape[i][2] = 0.5 * zvalue;
-
+ else if (keyword == SHAPE)
+ avec_ellipsoid->set_bonus(i,0.5*xvalue,0.5*yvalue,0.5*zvalue);
+
// set rmass via density
// if radius > 0.0, assume sphere
// if shape > 0.0, assume ellipsoid
// else set rmass to density directly
- } else if (keyword == DENSITY) {
+ else if (keyword == DENSITY) {
if (atom->radius_flag && atom->radius[i] > 0.0)
atom->rmass[i] = 4.0*PI/3.0 *
atom->radius[i]*atom->radius[i]*atom->radius[i] * dvalue;
- else if (atom->shape_flag && atom->shape[i][0] > 0.0)
- atom->rmass[i] = 4.0*PI/3.0 *
- atom->shape[i][0]*atom->shape[i][1]*atom->shape[i][2] * dvalue;
- else atom->rmass[i] = dvalue;
+ else if (atom->ellipsoid_flag && atom->ellipsoid[i] >= 0) {
+ double *shape = avec_ellipsoid->bonus[atom->ellipsoid[i]].shape;
+ atom->rmass[i] = 4.0*PI/3.0 * shape[0]*shape[1]*shape[2] * dvalue;
+ } else atom->rmass[i] = dvalue;
// reset any or all of 3 image flags
@@ -419,17 +420,18 @@ void Set::set(int keyword)
mu[i][2]*mu[i][2]);
// set quaternion orientation
-
+
} else if (keyword == QUAT) {
- double PI = 4.0*atan(1.0);
+ if (atom->ellipsoid[i] < 0)
+ error->one("Cannot set quaternion for atom that is not an ellipsoid");
+ double *quat = avec_ellipsoid->bonus[atom->ellipsoid[i]].quat;
double theta2 = 0.5 * PI * wvalue/180.0;
double sintheta2 = sin(theta2);
- double **quat = atom->quat;
- quat[i][0] = cos(theta2);
- quat[i][1] = xvalue * sintheta2;
- quat[i][2] = yvalue * sintheta2;
- quat[i][3] = zvalue * sintheta2;
- MathExtra::normalize4(quat[i]);
+ quat[0] = cos(theta2);
+ quat[1] = xvalue * sintheta2;
+ quat[2] = yvalue * sintheta2;
+ quat[3] = zvalue * sintheta2;
+ MathExtra::normalize4(quat);
}
count++;
}
@@ -509,38 +511,48 @@ void Set::setrandom(int keyword)
// no need to normalize quats since creations algorithms already do
} else if (keyword == QUAT_RANDOM) {
- double **quat = atom->quat;
+ AtomVecEllipsoid *avec_ellipsoid =
+ (AtomVecEllipsoid *) atom->style_match("ellipsoid");
+ AtomVecEllipsoid::Bonus *bonus = avec_ellipsoid->bonus;
+ int *ellipsoid = atom->ellipsoid;
int nlocal = atom->nlocal;
+ double *quat;
if (domain->dimension == 3) {
double s,t1,t2,theta1,theta2;
- double PI = 4.0*atan(1.0);
for (i = 0; i < nlocal; i++)
if (select[i]) {
+ if (ellipsoid[i] < 0)
+ error->one("Cannot set quaternion for atom "
+ "that is not an ellipsoid");
+ quat = bonus[ellipsoid[i]].quat;
random->reset(seed,x[i]);
s = random->uniform();
t1 = sqrt(1.0-s);
t2 = sqrt(s);
theta1 = 2.0*PI*random->uniform();
theta2 = 2.0*PI*random->uniform();
- quat[i][0] = cos(theta2)*t2;
- quat[i][1] = sin(theta1)*t1;
- quat[i][2] = cos(theta1)*t1;
- quat[i][3] = sin(theta2)*t2;
+ quat[0] = cos(theta2)*t2;
+ quat[1] = sin(theta1)*t1;
+ quat[2] = cos(theta1)*t1;
+ quat[3] = sin(theta2)*t2;
count++;
}
} else {
double theta2;
- double PI = 4.0*atan(1.0);
for (i = 0; i < nlocal; i++)
if (select[i]) {
+ if (ellipsoid[i] < 0)
+ error->one("Cannot set quaternion for atom "
+ "that is not an ellipsoid");
+ quat = bonus[ellipsoid[i]].quat;
random->reset(seed,x[i]);
theta2 = PI*random->uniform();
- quat[i][0] = cos(theta2);
- quat[i][1] = 0.0;
- quat[i][2] = 0.0;
- quat[i][3] = sin(theta2);
+ quat[0] = cos(theta2);
+ quat[1] = 0.0;
+ quat[2] = 0.0;
+ quat[3] = sin(theta2);
count++;
}
}