diff --git a/src/atom_vec.h b/src/atom_vec.h index 71ffb1dd46..dd2b056ef7 100644 --- a/src/atom_vec.h +++ b/src/atom_vec.h @@ -40,7 +40,8 @@ class AtomVec : protected Pointers { AtomVec(class LAMMPS *, int, char **); virtual ~AtomVec() {} - virtual void init() {} + void init(); + virtual void init_style() {} virtual void grow(int) = 0; virtual void grow_reset() = 0; @@ -82,6 +83,9 @@ class AtomVec : protected Pointers { protected: int nmax; // local copy of atom->nmax + int deform_vremap; // local copy of domain properties + int deform_groupbit; + double *h_rate; }; } diff --git a/src/atom_vec_charge.cpp b/src/atom_vec_charge.cpp index ed2053a52f..3071d5cf92 100644 --- a/src/atom_vec_charge.cpp +++ b/src/atom_vec_charge.cpp @@ -155,7 +155,7 @@ int AtomVecCharge::pack_comm_vel(int n, int *list, double *buf, int pbc_flag, int *pbc) { int i,j,m; - double dx,dy,dz; + double dx,dy,dz,dvx,dvy,dvz; m = 0; if (pbc_flag == 0) { @@ -178,14 +178,35 @@ int AtomVecCharge::pack_comm_vel(int n, int *list, double *buf, dy = pbc[1]*domain->yprd + pbc[3]*domain->yz; dz = pbc[2]*domain->zprd; } - for (i = 0; i < n; i++) { - j = list[i]; - buf[m++] = x[j][0] + dx; - buf[m++] = x[j][1] + dy; - buf[m++] = x[j][2] + dz; - buf[m++] = v[j][0]; - buf[m++] = v[j][1]; - buf[m++] = v[j][2]; + if (!deform_vremap) { + for (i = 0; i < n; i++) { + j = list[i]; + buf[m++] = x[j][0] + dx; + buf[m++] = x[j][1] + dy; + buf[m++] = x[j][2] + dz; + buf[m++] = v[j][0]; + buf[m++] = v[j][1]; + buf[m++] = v[j][2]; + } + } else { + dvx = pbc[0]*h_rate[0] + pbc[5]*h_rate[5] + pbc[4]*h_rate[4]; + dvy = pbc[1]*h_rate[1] + pbc[3]*h_rate[3]; + dvz = pbc[2]*h_rate[2]; + for (i = 0; i < n; i++) { + j = list[i]; + buf[m++] = x[j][0] + dx; + buf[m++] = x[j][1] + dy; + buf[m++] = x[j][2] + dz; + if (mask[i] & deform_groupbit) { + buf[m++] = v[j][0] + dvx; + buf[m++] = v[j][1] + dvy; + buf[m++] = v[j][2] + dvz; + } else { + buf[m++] = v[j][0]; + buf[m++] = v[j][1]; + buf[m++] = v[j][2]; + } + } } } return m; @@ -305,7 +326,7 @@ int AtomVecCharge::pack_border_vel(int n, int *list, double *buf, int pbc_flag, int *pbc) { int i,j,m; - double dx,dy,dz; + double dx,dy,dz,dvx,dvy,dvz; m = 0; if (pbc_flag == 0) { @@ -332,18 +353,43 @@ int AtomVecCharge::pack_border_vel(int n, int *list, double *buf, dy = pbc[1]; dz = pbc[2]; } - for (i = 0; i < n; i++) { - j = list[i]; - buf[m++] = x[j][0] + dx; - buf[m++] = x[j][1] + dy; - buf[m++] = x[j][2] + dz; - buf[m++] = tag[j]; - buf[m++] = type[j]; - buf[m++] = mask[j]; - buf[m++] = q[j]; - buf[m++] = v[j][0]; - buf[m++] = v[j][1]; - buf[m++] = v[j][2]; + if (!deform_vremap) { + for (i = 0; i < n; i++) { + j = list[i]; + buf[m++] = x[j][0] + dx; + buf[m++] = x[j][1] + dy; + buf[m++] = x[j][2] + dz; + buf[m++] = tag[j]; + buf[m++] = type[j]; + buf[m++] = mask[j]; + buf[m++] = q[j]; + buf[m++] = v[j][0]; + buf[m++] = v[j][1]; + buf[m++] = v[j][2]; + } + } else { + dvx = pbc[0]*h_rate[0] + pbc[5]*h_rate[5] + pbc[4]*h_rate[4]; + dvy = pbc[1]*h_rate[1] + pbc[3]*h_rate[3]; + dvz = pbc[2]*h_rate[2]; + for (i = 0; i < n; i++) { + j = list[i]; + buf[m++] = x[j][0] + dx; + buf[m++] = x[j][1] + dy; + buf[m++] = x[j][2] + dz; + buf[m++] = tag[j]; + buf[m++] = type[j]; + buf[m++] = mask[j]; + buf[m++] = q[j]; + if (mask[i] & deform_groupbit) { + buf[m++] = v[j][0] + dvx; + buf[m++] = v[j][1] + dvy; + buf[m++] = v[j][2] + dvz; + } else { + buf[m++] = v[j][0]; + buf[m++] = v[j][1]; + buf[m++] = v[j][2]; + } + } } } return m; diff --git a/src/domain.cpp b/src/domain.cpp index ac967cb211..be3802be86 100644 --- a/src/domain.cpp +++ b/src/domain.cpp @@ -107,14 +107,14 @@ void Domain::init() // check for fix deform - deform_flag = deform_remap = 0; + deform_flag = deform_vremap = deform_groupbit = 0; for (int i = 0; i < modify->nfix; i++) if (strcmp(modify->fix[i]->style,"deform") == 0) { deform_flag = 1; if (((FixDeform *) modify->fix[i])->remapflag == V_REMAP) { - deform_remap = 1; + deform_vremap = 1; deform_groupbit = modify->fix[i]->groupbit; - } else deform_remap = 0; + } } // region inits @@ -384,7 +384,7 @@ void Domain::pbc() if (xperiodic) { if (x[i][0] < lo[0]) { x[i][0] += period[0]; - if (deform_remap && mask[i] & deform_groupbit) v[i][0] += h_rate[0]; + if (deform_vremap && mask[i] & deform_groupbit) v[i][0] += h_rate[0]; idim = image[i] & 1023; otherdims = image[i] ^ idim; idim--; @@ -394,7 +394,7 @@ void Domain::pbc() if (x[i][0] >= hi[0]) { x[i][0] -= period[0]; x[i][0] = MAX(x[i][0],lo[0]); - if (deform_remap && mask[i] & deform_groupbit) v[i][0] -= h_rate[0]; + if (deform_vremap && mask[i] & deform_groupbit) v[i][0] -= h_rate[0]; idim = image[i] & 1023; otherdims = image[i] ^ idim; idim++; @@ -406,7 +406,7 @@ void Domain::pbc() if (yperiodic) { if (x[i][1] < lo[1]) { x[i][1] += period[1]; - if (deform_remap && mask[i] & deform_groupbit) { + if (deform_vremap && mask[i] & deform_groupbit) { v[i][0] += h_rate[5]; v[i][1] += h_rate[1]; } @@ -419,7 +419,7 @@ void Domain::pbc() if (x[i][1] >= hi[1]) { x[i][1] -= period[1]; x[i][1] = MAX(x[i][1],lo[1]); - if (deform_remap && mask[i] & deform_groupbit) { + if (deform_vremap && mask[i] & deform_groupbit) { v[i][0] -= h_rate[5]; v[i][1] -= h_rate[1]; } @@ -434,7 +434,7 @@ void Domain::pbc() if (zperiodic) { if (x[i][2] < lo[2]) { x[i][2] += period[2]; - if (deform_remap && mask[i] & deform_groupbit) { + if (deform_vremap && mask[i] & deform_groupbit) { v[i][0] += h_rate[4]; v[i][1] += h_rate[3]; v[i][2] += h_rate[2]; @@ -448,7 +448,7 @@ void Domain::pbc() if (x[i][2] >= hi[2]) { x[i][2] -= period[2]; x[i][2] = MAX(x[i][2],lo[2]); - if (deform_remap && mask[i] & deform_groupbit) { + if (deform_vremap && mask[i] & deform_groupbit) { v[i][0] -= h_rate[4]; v[i][1] -= h_rate[3]; v[i][2] -= h_rate[2]; diff --git a/src/domain.h b/src/domain.h index a630416708..d7cb88ce5c 100644 --- a/src/domain.h +++ b/src/domain.h @@ -75,7 +75,7 @@ class Domain : protected Pointers { int box_change; // 1 if box bounds ever change, 0 if fixed int deform_flag; // 1 if fix deform exist, else 0 - int deform_remap; // 1 if fix deform remaps v, else 0 + int deform_vremap; // 1 if fix deform remaps v, else 0 int deform_groupbit; // atom group to perform v remap for class Lattice *lattice; // user-defined lattice diff --git a/src/lammps.cpp b/src/lammps.cpp index 6b1d2f704d..4a52d9fa0e 100644 --- a/src/lammps.cpp +++ b/src/lammps.cpp @@ -331,10 +331,11 @@ void LAMMPS::init() update->init(); force->init(); // pair must come after update due to minimizer domain->init(); - atom->init(); // atom must come after force: + atom->init(); // atom must come after force and domain // atom deletes extra array // used by fix shear_history::unpack_restart() // when force->pair->gran_history creates fix ?? + // atom_vec init uses deform_vremap modify->init(); // modify must come after update, force, atom, domain neighbor->init(); // neighbor must come after force, modify comm->init(); // comm must come after force, modify, neighbor