diff --git a/src/ASPHERE/atom_vec_ellipsoid.cpp b/src/ASPHERE/atom_vec_ellipsoid.cpp index 4ddbef607c..e6c7a127e4 100755 --- a/src/ASPHERE/atom_vec_ellipsoid.cpp +++ b/src/ASPHERE/atom_vec_ellipsoid.cpp @@ -37,19 +37,21 @@ AtomVecEllipsoid::AtomVecEllipsoid(LAMMPS *lmp, int narg, char **arg) : AtomVec(lmp, narg, arg) { molecular = 0; - mass_type = 1; - shape_type = 1; comm_x_only = comm_f_only = 0; size_forward = 7; size_reverse = 6; - size_border = 10; + size_border = 13; size_velocity = 6; - size_data_atom = 9; + size_data_atom = 13; size_data_vel = 7; - xcol_data = 3; + xcol_data = 7; - atom->angmom_flag = atom->torque_flag = atom->quat_flag = 1; + atom->ellipsoid_flag = 1; + atom->shape_flag = atom->rmass_flag = atom->quat_flag = + atom->angmom_flag = atom->torque_flag = 1; + + PI = 4.0*atan(1.0); } /* ---------------------------------------------------------------------- @@ -74,6 +76,8 @@ 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"); @@ -92,6 +96,7 @@ 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; } @@ -110,6 +115,10 @@ 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]; @@ -392,6 +401,9 @@ 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]; @@ -415,6 +427,9 @@ 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]; @@ -442,6 +457,9 @@ 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]; @@ -472,6 +490,9 @@ 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]; @@ -495,6 +516,9 @@ 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]; @@ -521,11 +545,14 @@ int AtomVecEllipsoid::pack_border_vel(int n, int *list, double *buf, int AtomVecEllipsoid::pack_border_one(int i, double *buf) { - buf[0] = quat[i][0]; - buf[1] = quat[i][1]; - buf[2] = quat[i][2]; - buf[3] = quat[i][3]; - return 4; + buf[0] = shape[i][0]; + buf[1] = shape[i][1]; + buf[2] = shape[i][2]; + buf[3] = quat[i][0]; + buf[4] = quat[i][1]; + buf[5] = quat[i][2]; + buf[6] = quat[i][3]; + return 7; } /* ---------------------------------------------------------------------- */ @@ -544,6 +571,9 @@ void AtomVecEllipsoid::unpack_border(int n, int first, double *buf) tag[i] = static_cast (buf[m++]); type[i] = static_cast (buf[m++]); mask[i] = static_cast (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++]; @@ -567,6 +597,9 @@ void AtomVecEllipsoid::unpack_border_vel(int n, int first, double *buf) tag[i] = static_cast (buf[m++]); type[i] = static_cast (buf[m++]); mask[i] = static_cast (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++]; @@ -584,11 +617,14 @@ void AtomVecEllipsoid::unpack_border_vel(int n, int first, double *buf) int AtomVecEllipsoid::unpack_border_one(int i, double *buf) { - quat[i][0] = buf[0]; - quat[i][1] = buf[1]; - quat[i][2] = buf[2]; - quat[i][3] = buf[3]; - return 4; + shape[i][0] = buf[0]; + shape[i][1] = buf[1]; + shape[i][2] = buf[2]; + quat[i][0] = buf[3]; + quat[i][1] = buf[4]; + quat[i][2] = buf[5]; + quat[i][3] = buf[6]; + return 7; } /* ---------------------------------------------------------------------- @@ -610,6 +646,10 @@ 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]; @@ -645,6 +685,10 @@ int AtomVecEllipsoid::unpack_exchange(double *buf) mask[nlocal] = static_cast (buf[m++]); image[nlocal] = static_cast (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++]; @@ -672,7 +716,7 @@ int AtomVecEllipsoid::size_restart() int i; int nlocal = atom->nlocal; - int n = 18 * nlocal; + int n = 22 * nlocal; if (atom->nextra_restart) for (int iextra = 0; iextra < atom->nextra_restart; iextra++) @@ -702,6 +746,10 @@ 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]; @@ -743,6 +791,10 @@ 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++]; @@ -781,6 +833,12 @@ void AtomVecEllipsoid::create_atom(int itype, double *coord) v[nlocal][0] = 0.0; 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; @@ -811,14 +869,36 @@ 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"); + } + + double density = atof(values[5]); + if (density <= 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[5]); - quat[nlocal][1] = atof(values[6]); - quat[nlocal][2] = atof(values[7]); - quat[nlocal][3] = atof(values[8]); + 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; @@ -841,13 +921,35 @@ void AtomVecEllipsoid::data_atom(double *coord, int imagetmp, char **values) int AtomVecEllipsoid::data_atom_hybrid(int nlocal, char **values) { - quat[nlocal][0] = atof(values[0]); - quat[nlocal][1] = atof(values[1]); - quat[nlocal][2] = atof(values[2]); - quat[nlocal][3] = atof(values[3]); + 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"); + } + + double density = atof(values[3]); + if (density <= 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; + + 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]); - return 4; + return 8; } /* ---------------------------------------------------------------------- @@ -892,6 +994,8 @@ 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); diff --git a/src/ASPHERE/atom_vec_ellipsoid.h b/src/ASPHERE/atom_vec_ellipsoid.h index 830e225003..e3f0938047 100755 --- a/src/ASPHERE/atom_vec_ellipsoid.h +++ b/src/ASPHERE/atom_vec_ellipsoid.h @@ -60,8 +60,10 @@ class AtomVecEllipsoid : public AtomVec { bigint memory_usage(); private: + double PI; int *tag,*type,*mask,*image; double **x,**v,**f; + double **shape,*density,*rmass; double **angmom,**torque,**quat; }; diff --git a/src/ASPHERE/compute_erotate_asphere.cpp b/src/ASPHERE/compute_erotate_asphere.cpp index b7c6302eba..b08a5553df 100644 --- a/src/ASPHERE/compute_erotate_asphere.cpp +++ b/src/ASPHERE/compute_erotate_asphere.cpp @@ -34,24 +34,10 @@ ComputeERotateAsphere(LAMMPS *lmp, int narg, char **arg) : scalar_flag = 1; extscalar = 1; - memory->create(inertia,atom->ntypes+1,3,"fix_temp_sphere:inertia"); + // error check - // error checks - - if (!atom->angmom_flag || !atom->quat_flag || - !atom->avec->shape_type) - error->all("Compute erotate/asphere requires atom attributes " - "angmom, quat, shape"); - if (atom->radius_flag || atom->rmass_flag) - error->all("Compute erotate/asphere cannot be used with atom attributes " - "diameter or rmass"); -} - -/* ---------------------------------------------------------------------- */ - -ComputeERotateAsphere::~ComputeERotateAsphere() -{ - memory->destroy(inertia); + if (!atom->ellipsoid_flag) + error->all("Compute erotate/asphere requires atom style ellipsoid"); } /* ---------------------------------------------------------------------- */ @@ -62,77 +48,62 @@ void ComputeERotateAsphere::init() // no point particles allowed, spherical is OK double **shape = atom->shape; - int *type = atom->type; int *mask = atom->mask; int nlocal = atom->nlocal; for (int i = 0; i < nlocal; i++) if (mask[i] & groupbit) - if (shape[type[i]][0] == 0.0) + if (shape[i][0] == 0.0) error->one("Compute erotate/asphere requires extended particles"); pfactor = 0.5 * force->mvv2e; - calculate_inertia(); } /* ---------------------------------------------------------------------- */ double ComputeERotateAsphere::compute_scalar() { - int i,itype; - invoked_scalar = update->ntimestep; double **quat = atom->quat; double **angmom = atom->angmom; + double **shape = atom->shape; + double *rmass = atom->rmass; int *mask = atom->mask; - int *type = atom->type; int nlocal = atom->nlocal; // sum rotational energy for each particle // no point particles since divide by inertia - double wbody[3]; + double wbody[3],inertia[3]; double rot[3][3]; double erotate = 0.0; - for (i = 0; i < nlocal; i++) + for (int i = 0; i < nlocal; i++) if (mask[i] & groupbit) { - itype = type[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; // wbody = angular velocity in body frame MathExtra::quat_to_mat(quat[i],rot); MathExtra::transpose_times_column3(rot,angmom[i],wbody); - wbody[0] /= inertia[itype][0]; - wbody[1] /= inertia[itype][1]; - wbody[2] /= inertia[itype][2]; + wbody[0] /= inertia[0]; + wbody[1] /= inertia[1]; + wbody[2] /= inertia[2]; - erotate += inertia[itype][0]*wbody[0]*wbody[0]+ - inertia[itype][1]*wbody[1]*wbody[1]+ - inertia[itype][2]*wbody[2]*wbody[2]; + erotate += inertia[0]*wbody[0]*wbody[0] + + inertia[1]*wbody[1]*wbody[1] + inertia[2]*wbody[2]*wbody[2]; } MPI_Allreduce(&erotate,&scalar,1,MPI_DOUBLE,MPI_SUM,world); scalar *= pfactor; return scalar; } - -/* ---------------------------------------------------------------------- - principal moments of inertia for ellipsoids -------------------------------------------------------------------------- */ - -void ComputeERotateAsphere::calculate_inertia() -{ - double *mass = atom->mass; - double **shape = atom->shape; - - for (int i = 1; i <= atom->ntypes; i++) { - inertia[i][0] = mass[i] * - (shape[i][1]*shape[i][1]+shape[i][2]*shape[i][2]) / 5.0; - inertia[i][1] = mass[i] * - (shape[i][0]*shape[i][0]+shape[i][2]*shape[i][2]) / 5.0; - inertia[i][2] = mass[i] * - (shape[i][0]*shape[i][0]+shape[i][1]*shape[i][1]) / 5.0; - } -} diff --git a/src/ASPHERE/compute_erotate_asphere.h b/src/ASPHERE/compute_erotate_asphere.h index fa0da01f6e..f7f31066a8 100644 --- a/src/ASPHERE/compute_erotate_asphere.h +++ b/src/ASPHERE/compute_erotate_asphere.h @@ -27,15 +27,11 @@ namespace LAMMPS_NS { class ComputeERotateAsphere : public Compute { public: ComputeERotateAsphere(class LAMMPS *, int, char **); - ~ComputeERotateAsphere(); void init(); double compute_scalar(); private: double pfactor; - double **inertia; - - void calculate_inertia(); }; } diff --git a/src/ASPHERE/compute_temp_asphere.cpp b/src/ASPHERE/compute_temp_asphere.cpp index 3a0ccbee4a..df29ad7b47 100755 --- a/src/ASPHERE/compute_temp_asphere.cpp +++ b/src/ASPHERE/compute_temp_asphere.cpp @@ -55,17 +55,11 @@ ComputeTempAsphere::ComputeTempAsphere(LAMMPS *lmp, int narg, char **arg) : } vector = new double[6]; - memory->create(inertia,atom->ntypes+1,3,"fix_temp_sphere:inertia"); // error checks - if (!atom->angmom_flag || !atom->quat_flag || - !atom->avec->shape_type) - error->all("Compute temp/asphere requires atom attributes " - "angmom, quat, shape"); - if (atom->radius_flag || atom->rmass_flag) - error->all("Compute temp/asphere cannot be used with atom attributes " - "diameter or rmass"); + if (!atom->ellipsoid_flag) + error->all("Compute temp/asphere requires atom style ellipsoid"); } /* ---------------------------------------------------------------------- */ @@ -74,7 +68,6 @@ ComputeTempAsphere::~ComputeTempAsphere() { delete [] id_bias; delete [] vector; - memory->destroy(inertia); } /* ---------------------------------------------------------------------- */ @@ -85,13 +78,12 @@ void ComputeTempAsphere::init() // no point particles allowed, spherical is OK double **shape = atom->shape; - int *type = atom->type; int *mask = atom->mask; int nlocal = atom->nlocal; for (int i = 0; i < nlocal; i++) if (mask[i] & groupbit) - if (shape[type[i]][0] == 0.0) + if (shape[i][0] == 0.0) error->one("Compute temp/asphere requires extended particles"); if (tempbias) { @@ -113,8 +105,6 @@ void ComputeTempAsphere::init() for (int i = 0; i < modify->nfix; i++) fix_dof += modify->fix[i]->dof(igroup); dof_compute(); - - calculate_inertia(); } /* ---------------------------------------------------------------------- */ @@ -164,13 +154,12 @@ double ComputeTempAsphere::compute_scalar() double **v = atom->v; double **quat = atom->quat; double **angmom = atom->angmom; - double *mass = atom->mass; - int *type = atom->type; + double **shape = atom->shape; + double *rmass = atom->rmass; int *mask = atom->mask; int nlocal = atom->nlocal; - int itype; - double wbody[3]; + double wbody[3],inertia[3]; double rot[3][3]; double t = 0.0; @@ -180,20 +169,27 @@ double ComputeTempAsphere::compute_scalar() for (int i = 0; i < nlocal; i++) if (mask[i] & groupbit) { - itype = type[i]; - t += (v[i][0]*v[i][0] + v[i][1]*v[i][1] + v[i][2]*v[i][2]) * mass[itype]; + 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; // wbody = angular velocity in body frame MathExtra::quat_to_mat(quat[i],rot); MathExtra::transpose_times_column3(rot,angmom[i],wbody); - wbody[0] /= inertia[itype][0]; - wbody[1] /= inertia[itype][1]; - wbody[2] /= inertia[itype][2]; + wbody[0] /= inertia[0]; + wbody[1] /= inertia[1]; + wbody[2] /= inertia[2]; - t += inertia[itype][0]*wbody[0]*wbody[0] + - inertia[itype][1]*wbody[1]*wbody[1] + - inertia[itype][2]*wbody[2]*wbody[2]; + t += inertia[0]*wbody[0]*wbody[0] + + inertia[1]*wbody[1]*wbody[1] + inertia[2]*wbody[2]*wbody[2]; } if (tempbias) tbias->restore_bias_all(); @@ -220,13 +216,12 @@ void ComputeTempAsphere::compute_vector() double **v = atom->v; double **quat = atom->quat; double **angmom = atom->angmom; - double *mass = atom->mass; - int *type = atom->type; + double **shape = atom->shape; + double *rmass = atom->rmass; int *mask = atom->mask; int nlocal = atom->nlocal; - int itype; - double wbody[3]; + double wbody[3],inertia[3]; double rot[3][3]; double massone,t[6]; for (i = 0; i < 6; i++) t[i] = 0.0; @@ -236,8 +231,7 @@ void ComputeTempAsphere::compute_vector() // translational kinetic energy - itype = type[i]; - massone = mass[itype]; + massone = rmass[i]; t[0] += massone * v[i][0]*v[i][0]; t[1] += massone * v[i][1]*v[i][1]; t[2] += massone * v[i][2]*v[i][2]; @@ -245,22 +239,31 @@ void ComputeTempAsphere::compute_vector() t[4] += massone * v[i][0]*v[i][2]; t[5] += massone * v[i][1]*v[i][2]; + // 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; + // wbody = angular velocity in body frame MathExtra::quat_to_mat(quat[i],rot); MathExtra::transpose_times_column3(rot,angmom[i],wbody); - wbody[0] /= inertia[itype][0]; - wbody[1] /= inertia[itype][1]; - wbody[2] /= inertia[itype][2]; + wbody[0] /= inertia[0]; + wbody[1] /= inertia[1]; + wbody[2] /= inertia[2]; // rotational kinetic energy - t[0] += inertia[itype][0]*wbody[0]*wbody[0]; - t[1] += inertia[itype][1]*wbody[1]*wbody[1]; - t[2] += inertia[itype][2]*wbody[2]*wbody[2]; - t[3] += inertia[itype][0]*wbody[0]*wbody[1]; - t[4] += inertia[itype][1]*wbody[0]*wbody[2]; - t[5] += inertia[itype][2]*wbody[1]*wbody[2]; + t[0] += inertia[0]*wbody[0]*wbody[0]; + t[1] += inertia[1]*wbody[1]*wbody[1]; + t[2] += inertia[2]*wbody[2]*wbody[2]; + t[3] += inertia[0]*wbody[0]*wbody[1]; + t[4] += inertia[1]*wbody[0]*wbody[2]; + t[5] += inertia[2]*wbody[1]*wbody[2]; } if (tempbias) tbias->restore_bias_all(); @@ -269,25 +272,6 @@ void ComputeTempAsphere::compute_vector() for (i = 0; i < 6; i++) vector[i] *= force->mvv2e; } -/* ---------------------------------------------------------------------- - principal moments of inertia for ellipsoids -------------------------------------------------------------------------- */ - -void ComputeTempAsphere::calculate_inertia() -{ - double *mass = atom->mass; - double **shape = atom->shape; - - for (int i = 1; i <= atom->ntypes; i++) { - inertia[i][0] = mass[i] * - (shape[i][1]*shape[i][1]+shape[i][2]*shape[i][2]) / 5.0; - inertia[i][1] = mass[i] * - (shape[i][0]*shape[i][0]+shape[i][2]*shape[i][2]) / 5.0; - inertia[i][2] = mass[i] * - (shape[i][0]*shape[i][0]+shape[i][1]*shape[i][1]) / 5.0; - } -} - /* ---------------------------------------------------------------------- remove velocity bias from atom I to leave thermal velocity ------------------------------------------------------------------------- */ diff --git a/src/ASPHERE/compute_temp_asphere.h b/src/ASPHERE/compute_temp_asphere.h index 17ffbd57ad..1504791e1b 100755 --- a/src/ASPHERE/compute_temp_asphere.h +++ b/src/ASPHERE/compute_temp_asphere.h @@ -38,12 +38,10 @@ class ComputeTempAsphere : public Compute { private: int fix_dof; double tfactor; - double **inertia; char *id_bias; Compute *tbias; // ptr to additional bias compute void dof_compute(); - void calculate_inertia(); }; } diff --git a/src/ASPHERE/fix_nh_asphere.cpp b/src/ASPHERE/fix_nh_asphere.cpp index ddf5d1b8db..9a587e6367 100644 --- a/src/ASPHERE/fix_nh_asphere.cpp +++ b/src/ASPHERE/fix_nh_asphere.cpp @@ -33,22 +33,8 @@ using namespace LAMMPS_NS; FixNHAsphere::FixNHAsphere(LAMMPS *lmp, int narg, char **arg) : FixNH(lmp, narg, arg) { - memory->create(inertia,atom->ntypes+1,3,"fix_nvt_asphere:inertia"); - - if (!atom->quat_flag || !atom->angmom_flag || !atom->torque_flag || - !atom->avec->shape_type) - error->all("Fix nvt/nph/npt asphere requires atom attributes " - "quat, angmom, torque, shape"); - if (atom->radius_flag || atom->rmass_flag) - error->all("Fix nvt/nph/npt asphere cannot be used with atom attributes " - "diameter or rmass"); -} - -/* ---------------------------------------------------------------------- */ - -FixNHAsphere::~FixNHAsphere() -{ - memory->destroy(inertia); + if (!atom->ellipsoid_flag) + error->all("Fix nvt/nph/npt asphere requires atom style ellipsoid"); } /* ---------------------------------------------------------------------- */ @@ -59,17 +45,15 @@ void FixNHAsphere::init() // no point particles allowed, spherical is OK double **shape = atom->shape; - int *type = atom->type; int *mask = atom->mask; int nlocal = atom->nlocal; for (int i = 0; i < nlocal; i++) if (mask[i] & groupbit) - if (shape[type[i]][0] == 0.0) + if (shape[i][0] == 0.0) error->one("Fix nvt/nph/npt asphere requires extended particles"); FixNH::init(); - calculate_inertia(); } /* ---------------------------------------------------------------------- @@ -149,25 +133,6 @@ void FixNHAsphere::omega_from_mq(double *q, double *m, double *inertia, MathExtra::times_column3(rot,wbody,w); } -/* ---------------------------------------------------------------------- - principal moments of inertia for ellipsoids -------------------------------------------------------------------------- */ - -void FixNHAsphere::calculate_inertia() -{ - double *mass = atom->mass; - double **shape = atom->shape; - - for (int i = 1; i <= atom->ntypes; i++) { - inertia[i][0] = 0.2*mass[i] * - (shape[i][1]*shape[i][1]+shape[i][2]*shape[i][2]); - inertia[i][1] = 0.2*mass[i] * - (shape[i][0]*shape[i][0]+shape[i][2]*shape[i][2]); - inertia[i][2] = 0.2*mass[i] * - (shape[i][0]*shape[i][0]+shape[i][1]*shape[i][1]); - } -} - /* ---------------------------------------------------------------------- perform half-step update of angular momentum -----------------------------------------------------------------------*/ @@ -207,7 +172,8 @@ void FixNHAsphere::nve_x() double **quat = atom->quat; double **angmom = atom->angmom; - int *type = atom->type; + double **shape = atom->shape; + double *rmass = atom->rmass; int *mask = atom->mask; int nlocal = atom->nlocal; if (igroup == atom->firstgroup) nlocal = atom->nfirst; @@ -218,11 +184,21 @@ void FixNHAsphere::nve_x() // update quaternion a full step via Richardson iteration // returns new normalized quaternion + // principal moments of inertia - for (int i = 0; i < nlocal; i++) { - if (mask[i] & groupbit) - richardson(quat[i],angmom[i],inertia[type[i]]); - } + 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; + + richardson(quat[i],angmom[i],inertia); + } } /* ---------------------------------------------------------------------- diff --git a/src/ASPHERE/fix_nh_asphere.h b/src/ASPHERE/fix_nh_asphere.h index a35d18fdff..2d66264bba 100644 --- a/src/ASPHERE/fix_nh_asphere.h +++ b/src/ASPHERE/fix_nh_asphere.h @@ -21,7 +21,7 @@ namespace LAMMPS_NS { class FixNHAsphere : public FixNH { public: FixNHAsphere(class LAMMPS *, int, char **); - virtual ~FixNHAsphere(); + virtual ~FixNHAsphere() {} void init(); protected: diff --git a/src/ASPHERE/fix_nve_asphere.cpp b/src/ASPHERE/fix_nve_asphere.cpp index 4ddf05081d..fe5d71a475 100755 --- a/src/ASPHERE/fix_nve_asphere.cpp +++ b/src/ASPHERE/fix_nve_asphere.cpp @@ -34,24 +34,10 @@ using namespace LAMMPS_NS; FixNVEAsphere::FixNVEAsphere(LAMMPS *lmp, int narg, char **arg) : FixNVE(lmp, narg, arg) { - memory->create(inertia,atom->ntypes+1,3,"fix_nve_asphere:inertia"); - // error checks - if (!atom->angmom_flag || !atom->quat_flag || !atom->torque_flag || - !atom->avec->shape_type) - error->all("Fix nve/asphere requires atom attributes " - "angmom, quat, torque, shape"); - if (atom->radius_flag || atom->rmass_flag) - error->all("Fix nve/asphere cannot be used with atom attributes " - "diameter or rmass"); -} - -/* ---------------------------------------------------------------------- */ - -FixNVEAsphere::~FixNVEAsphere() -{ - memory->destroy(inertia); + if (!atom->ellipsoid_flag) + error->all("Fix nve/asphere requires atom style ellipsoid"); } /* ---------------------------------------------------------------------- */ @@ -62,17 +48,15 @@ void FixNVEAsphere::init() // no point particles allowed, spherical is OK double **shape = atom->shape; - int *type = atom->type; int *mask = atom->mask; int nlocal = atom->nlocal; for (int i = 0; i < nlocal; i++) if (mask[i] & groupbit) - if (shape[type[i]][0] == 0.0) + if (shape[i][0] == 0.0) error->one("Fix nve/asphere requires extended particles"); FixNVE::init(); - calculate_inertia(); } /* ---------------------------------------------------------------------- */ @@ -80,6 +64,7 @@ void FixNVEAsphere::init() void FixNVEAsphere::initial_integrate(int vflag) { double dtfm; + double inertia[3]; double **x = atom->x; double **v = atom->v; @@ -87,8 +72,8 @@ void FixNVEAsphere::initial_integrate(int vflag) double **quat = atom->quat; double **angmom = atom->angmom; double **torque = atom->torque; - double *mass = atom->mass; - int *type = atom->type; + double **shape = atom->shape; + double *rmass = atom->rmass; int *mask = atom->mask; int nlocal = atom->nlocal; if (igroup == atom->firstgroup) nlocal = atom->nfirst; @@ -99,7 +84,7 @@ void FixNVEAsphere::initial_integrate(int vflag) for (int i = 0; i < nlocal; i++) if (mask[i] & groupbit) { - dtfm = dtf / mass[type[i]]; + dtfm = dtf / rmass[i]; v[i][0] += dtfm * f[i][0]; v[i][1] += dtfm * f[i][1]; v[i][2] += dtfm * f[i][2]; @@ -115,7 +100,16 @@ void FixNVEAsphere::initial_integrate(int vflag) angmom[i][1] += dtf * torque[i][1]; angmom[i][2] += dtf * torque[i][2]; - richardson(quat[i],angmom[i],inertia[type[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; + + richardson(quat[i],angmom[i],inertia); } } @@ -129,15 +123,14 @@ void FixNVEAsphere::final_integrate() double **f = atom->f; double **angmom = atom->angmom; double **torque = atom->torque; - double *mass = atom->mass; - int *type = atom->type; + double *rmass = atom->rmass; int *mask = atom->mask; int nlocal = atom->nlocal; if (igroup == atom->firstgroup) nlocal = atom->nfirst; for (int i = 0; i < nlocal; i++) if (mask[i] & groupbit) { - dtfm = dtf / mass[type[i]]; + dtfm = dtf / rmass[i]; v[i][0] += dtfm * f[i][0]; v[i][1] += dtfm * f[i][1]; v[i][2] += dtfm * f[i][2]; @@ -224,22 +217,3 @@ void FixNVEAsphere::omega_from_mq(double *q, double *m, double *moments, wbody[2] /= moments[2]; MathExtra::times_column3(rot,wbody,w); } - -/* ---------------------------------------------------------------------- - principal moments of inertia for ellipsoids -------------------------------------------------------------------------- */ - -void FixNVEAsphere::calculate_inertia() -{ - double *mass = atom->mass; - double **shape = atom->shape; - - for (int i = 1; i <= atom->ntypes; i++) { - inertia[i][0] = 0.2*mass[i] * - (shape[i][1]*shape[i][1]+shape[i][2]*shape[i][2]); - inertia[i][1] = 0.2*mass[i] * - (shape[i][0]*shape[i][0]+shape[i][2]*shape[i][2]); - inertia[i][2] = 0.2*mass[i] * - (shape[i][0]*shape[i][0]+shape[i][1]*shape[i][1]); - } -} diff --git a/src/ASPHERE/fix_nve_asphere.h b/src/ASPHERE/fix_nve_asphere.h index 0fc29bbe75..b892377331 100755 --- a/src/ASPHERE/fix_nve_asphere.h +++ b/src/ASPHERE/fix_nve_asphere.h @@ -27,18 +27,15 @@ namespace LAMMPS_NS { class FixNVEAsphere : public FixNVE { public: FixNVEAsphere(class LAMMPS *, int, char **); - ~FixNVEAsphere(); void init(); void initial_integrate(int); void final_integrate(); private: double dtq; - double **inertia; void richardson(double *, double *, double *); void omega_from_mq(double *, double *, double *, double *); - void calculate_inertia(); }; } diff --git a/src/ASPHERE/pair_gayberne.cpp b/src/ASPHERE/pair_gayberne.cpp index 62fd2bbeae..5bc3c01199 100755 --- a/src/ASPHERE/pair_gayberne.cpp +++ b/src/ASPHERE/pair_gayberne.cpp @@ -58,7 +58,8 @@ PairGayBerne::~PairGayBerne() memory->destroy(form); memory->destroy(epsilon); memory->destroy(sigma); - memory->destroy(shape); + memory->destroy(shape1); + memory->destroy(shape2); memory->destroy(well); memory->destroy(cut); memory->destroy(lj1); @@ -109,7 +110,7 @@ void PairGayBerne::compute(int eflag, int vflag) MathExtra::quat_to_mat_trans(quat[i],a1); MathExtra::diag_times3(well[itype],a1,temp); MathExtra::transpose_times3(a1,temp,b1); - MathExtra::diag_times3(shape[itype],a1,temp); + MathExtra::diag_times3(shape2[itype],a1,temp); MathExtra::transpose_times3(a1,temp,g1); } @@ -153,7 +154,7 @@ void PairGayBerne::compute(int eflag, int vflag) MathExtra::quat_to_mat_trans(quat[j],a2); MathExtra::diag_times3(well[jtype],a2,temp); MathExtra::transpose_times3(a2,temp,b2); - MathExtra::diag_times3(shape[jtype],a2,temp); + MathExtra::diag_times3(shape2[jtype],a2,temp); MathExtra::transpose_times3(a2,temp,g2); one_eng = gayberne_lj(j,i,a2,b2,g2,r12,rsq,fforce,rtor); ttor[0] = ttor[1] = ttor[2] = 0.0; @@ -168,7 +169,7 @@ void PairGayBerne::compute(int eflag, int vflag) MathExtra::quat_to_mat_trans(quat[j],a2); MathExtra::diag_times3(well[jtype],a2,temp); MathExtra::transpose_times3(a2,temp,b2); - MathExtra::diag_times3(shape[jtype],a2,temp); + MathExtra::diag_times3(shape2[jtype],a2,temp); MathExtra::transpose_times3(a2,temp,g2); one_eng = gayberne_analytic(i,j,a1,a2,b1,b2,g1,g2,r12,rsq, fforce,ttor,rtor); @@ -232,7 +233,8 @@ void PairGayBerne::allocate() memory->create(form,n+1,n+1,"pair:form"); memory->create(epsilon,n+1,n+1,"pair:epsilon"); memory->create(sigma,n+1,n+1,"pair:sigma"); - memory->create(shape,n+1,3,"pair:shape"); + memory->create(shape1,n+1,3,"pair:shape1"); + memory->create(shape2,n+1,3,"pair:shape2"); memory->create(well,n+1,3,"pair:well"); memory->create(cut,n+1,n+1,"pair:cut"); memory->create(lj1,n+1,n+1,"pair:lj1"); @@ -330,22 +332,23 @@ void PairGayBerne::coeff(int narg, char **arg) void PairGayBerne::init_style() { - if (!atom->quat_flag || !atom->torque_flag || !atom->avec->shape_type) - error->all("Pair gayberne requires atom attributes quat, torque, shape"); - if (atom->radius_flag) - error->all("Pair gayberne cannot be used with atom attribute diameter"); + if (!atom->ellipsoid_flag) + error->all("Pair gayberne requires atom style ellipsoid"); int irequest = neighbor->request(this); // per-type shape precalculations + // require that atom shapes are identical within each type for (int i = 1; i <= atom->ntypes; i++) { + if (!atom->shape_consistency(i,shape1[i][0],shape1[i][1],shape1[i][2])) + error->all("Pair gayberne requires atoms with same type have same shape"); if (setwell[i]) { - double *one = atom->shape[i]; - shape[i][0] = one[0]*one[0]; - shape[i][1] = one[1]*one[1]; - shape[i][2] = one[2]*one[2]; - lshape[i] = (one[0]*one[1]+one[2]*one[2])*sqrt(one[0]*one[1]); + shape2[i][0] = shape1[i][0]*shape1[i][0]; + shape2[i][1] = shape1[i][1]*shape1[i][1]; + shape2[i][2] = shape1[i][2]*shape1[i][2]; + lshape[i] = (shape1[i][0]*shape1[i][1]+shape1[i][2]*shape1[i][2]) * + sqrt(shape1[i][0]*shape1[i][1]); } } } @@ -377,14 +380,14 @@ double PairGayBerne::init_one(int i, int j) } else offset[i][j] = 0.0; int ishape = 0; - if (atom->shape[i][0] != atom->shape[i][1] || - atom->shape[i][0] != atom->shape[i][2] || - atom->shape[i][1] != atom->shape[i][2]) ishape = 1; + if (shape1[i][0] != shape1[i][1] || + shape1[i][0] != shape1[i][2] || + shape1[i][1] != shape1[i][2]) ishape = 1; if (setwell[i] == 1) ishape = 1; int jshape = 0; - if (atom->shape[j][0] != atom->shape[j][1] || - atom->shape[j][0] != atom->shape[j][2] || - atom->shape[j][1] != atom->shape[j][2]) jshape = 1; + if (shape1[j][0] != shape1[j][1] || + shape1[j][0] != shape1[j][2] || + shape1[j][1] != shape1[j][2]) jshape = 1; if (setwell[j] == 1) jshape = 1; if (ishape == 0 && jshape == 0) @@ -640,7 +643,7 @@ double PairGayBerne::gayberne_analytic(const int i,const int j,double a1[3][3], double deta[3]; deta[0] = deta[1] = deta[2] = 0.0; - compute_eta_torque(g12,a1,shape[type[i]],temp); + compute_eta_torque(g12,a1,shape2[type[i]],temp); temp1 = -eta*upsilon; for (int m = 0; m < 3; m++) { for (int y = 0; y < 3; y++) tempv[y] = temp1*temp[m][y]; @@ -655,7 +658,7 @@ double PairGayBerne::gayberne_analytic(const int i,const int j,double a1[3][3], double deta2[3]; if (newton_pair || j < nlocal) { deta2[0] = deta2[1] = deta2[2] = 0.0; - compute_eta_torque(g12,a2,shape[type[j]],temp); + compute_eta_torque(g12,a2,shape2[type[j]],temp); for (int m = 0; m < 3; m++) { for (int y = 0; y < 3; y++) tempv[y] = temp1*temp[m][y]; MathExtra::cross3(a2[m],tempv,tempv2); @@ -707,9 +710,9 @@ double PairGayBerne::gayberne_lj(const int i,const int j,double a1[3][3], // compute distance of closest approach double g12[3][3]; - g12[0][0] = g1[0][0]+shape[type[j]][0]; - g12[1][1] = g1[1][1]+shape[type[j]][0]; - g12[2][2] = g1[2][2]+shape[type[j]][0]; + g12[0][0] = g1[0][0]+shape2[type[j]][0]; + g12[1][1] = g1[1][1]+shape2[type[j]][0]; + g12[2][2] = g1[2][2]+shape2[type[j]][0]; g12[0][1] = g1[0][1]; g12[1][0] = g1[1][0]; g12[0][2] = g1[0][2]; g12[2][0] = g1[2][0]; g12[1][2] = g1[1][2]; g12[2][1] = g1[2][1]; @@ -809,7 +812,7 @@ double PairGayBerne::gayberne_lj(const int i,const int j,double a1[3][3], double deta[3]; deta[0] = deta[1] = deta[2] = 0.0; - compute_eta_torque(g12,a1,shape[type[i]],temp); + compute_eta_torque(g12,a1,shape2[type[i]],temp); temp1 = -eta*upsilon; for (int m = 0; m < 3; m++) { for (int y = 0; y < 3; y++) tempv[y] = temp1*temp[m][y]; diff --git a/src/ASPHERE/pair_gayberne.h b/src/ASPHERE/pair_gayberne.h index 03cc0d7e01..3d465be97e 100755 --- a/src/ASPHERE/pair_gayberne.h +++ b/src/ASPHERE/pair_gayberne.h @@ -43,7 +43,8 @@ class PairGayBerne : public Pair { double **cut; double gamma,upsilon,mu; // Gay-Berne parameters - double **shape; // radii in x, y and z SQUARED + double **shape1; // per-type radii in x, y and z + double **shape2; // per-type radii in x, y and z SQUARED double *lshape; // precalculation based on the shape double **well; // well depth scaling along each axis ^ -1.0/mu double **epsilon,**sigma; // epsilon and sigma values for atom-type pairs diff --git a/src/ASPHERE/pair_resquared.cpp b/src/ASPHERE/pair_resquared.cpp index e61f8c1f3e..6a9b1723df 100755 --- a/src/ASPHERE/pair_resquared.cpp +++ b/src/ASPHERE/pair_resquared.cpp @@ -64,6 +64,7 @@ PairRESquared::~PairRESquared() memory->destroy(form); memory->destroy(epsilon); memory->destroy(sigma); + memory->destroy(shape1); memory->destroy(shape2); memory->destroy(well); memory->destroy(cut); @@ -224,6 +225,7 @@ void PairRESquared::allocate() memory->create(form,n+1,n+1,"pair:form"); memory->create(epsilon,n+1,n+1,"pair:epsilon"); memory->create(sigma,n+1,n+1,"pair:sigma"); + memory->create(shape1,n+1,3,"pair:shape1"); memory->create(shape2,n+1,3,"pair:shape2"); memory->create(well,n+1,3,"pair:well"); memory->create(cut,n+1,n+1,"pair:cut"); @@ -319,22 +321,22 @@ void PairRESquared::coeff(int narg, char **arg) void PairRESquared::init_style() { - if (!atom->quat_flag || !atom->torque_flag || !atom->avec->shape_type) - error->all("Pair resquared requires atom attributes quat, torque, shape"); - if (atom->radius_flag) - error->all("Pair resquared cannot be used with atom attribute diameter"); + if (!atom->ellipsoid_flag) + error->all("Pair resquared requires atom style ellipsoid"); int irequest = neighbor->request(this); // per-type shape precalculations + // require that atom shapes are identical within each type for (int i = 1; i <= atom->ntypes; i++) { + if (!atom->shape_consistency(i,shape1[i][0],shape1[i][1],shape1[i][2])) + error->all("Pair gayberne requires atoms with same type have same shape"); if (setwell[i]) { - double *one = atom->shape[i]; - shape2[i][0] = one[0]*one[0]; - shape2[i][1] = one[1]*one[1]; - shape2[i][2] = one[2]*one[2]; - lshape[i] = one[0]*one[1]*one[2]; + shape2[i][0] = shape1[i][0]*shape1[i][0]; + shape2[i][1] = shape1[i][1]*shape1[i][1]; + shape2[i][2] = shape1[i][2]*shape1[i][2]; + lshape[i] = shape1[i][0]*shape1[i][1]*shape1[i][2]; } } } @@ -345,16 +347,14 @@ void PairRESquared::init_style() double PairRESquared::init_one(int i, int j) { - double **shape = atom->shape; - if (setwell[i] == 0 || setwell[j] == 0) error->all("Pair resquared epsilon a,b,c coeffs are not all set"); int ishape = 0; - if (shape[i][0] != 0.0 && shape[i][1] != 0.0 && shape[i][2] != 0.0) + if (shape1[i][0] != 0.0 && shape1[i][1] != 0.0 && shape1[i][2] != 0.0) ishape = 1; int jshape = 0; - if (shape[j][0] != 0.0 && shape[j][1] != 0.0 && shape[j][2] != 0.0) + if (shape1[j][0] != 0.0 && shape1[j][1] != 0.0 && shape1[j][2] != 0.0) jshape = 1; if (ishape == 0 && jshape == 0) { @@ -546,7 +546,6 @@ double PairRESquared::resquared_analytic(const int i, const int j, double *rtor) { int *type = atom->type; - double **shape = atom->shape; // pair computations for energy, force, torque @@ -649,16 +648,16 @@ double PairRESquared::resquared_analytic(const int i, const int j, tprod = eta*chi*sigh; double stemp = h12/2.0; - Ua = (shape[type[i]][0]+stemp)*(shape[type[i]][1]+stemp)* - (shape[type[i]][2]+stemp)*(shape[type[j]][0]+stemp)* - (shape[type[j]][1]+stemp)*(shape[type[j]][2]+stemp); + Ua = (shape1[type[i]][0]+stemp)*(shape1[type[i]][1]+stemp)* + (shape1[type[i]][2]+stemp)*(shape1[type[j]][0]+stemp)* + (shape1[type[j]][1]+stemp)*(shape1[type[j]][2]+stemp); Ua = (1.0+3.0*tprod)*sprod/Ua; Ua = epsilon[type[i]][type[j]]*Ua/-36.0; stemp = h12/cr60; - Ur = (shape[type[i]][0]+stemp)*(shape[type[i]][1]+stemp)* - (shape[type[i]][2]+stemp)*(shape[type[j]][0]+stemp)* - (shape[type[j]][1]+stemp)*(shape[type[j]][2]+stemp); + Ur = (shape1[type[i]][0]+stemp)*(shape1[type[i]][1]+stemp)* + (shape1[type[i]][2]+stemp)*(shape1[type[j]][0]+stemp)* + (shape1[type[j]][1]+stemp)*(shape1[type[j]][2]+stemp); Ur = (1.0+b_alpha*tprod)*sprod/Ur; Ur = epsilon[type[i]][type[j]]*Ur*pow(sigh,6.0)/2025.0; @@ -704,22 +703,22 @@ double PairRESquared::resquared_analytic(const int i, const int j, spr[1] = 0.5*sigma12p3*s[1]; spr[2] = 0.5*sigma12p3*s[2]; - stemp = 1.0/(shape[type[i]][0]*2.0+h12)+ - 1.0/(shape[type[i]][1]*2.0+h12)+ - 1.0/(shape[type[i]][2]*2.0+h12)+ - 1.0/(shape[type[j]][0]*2.0+h12)+ - 1.0/(shape[type[j]][1]*2.0+h12)+ - 1.0/(shape[type[j]][2]*2.0+h12); + stemp = 1.0/(shape1[type[i]][0]*2.0+h12)+ + 1.0/(shape1[type[i]][1]*2.0+h12)+ + 1.0/(shape1[type[i]][2]*2.0+h12)+ + 1.0/(shape1[type[j]][0]*2.0+h12)+ + 1.0/(shape1[type[j]][1]*2.0+h12)+ + 1.0/(shape1[type[j]][2]*2.0+h12); hsec = h12+3.0*sec; dspu = 1.0/h12-1.0/hsec+stemp; pbsu = 3.0*sigma[type[i]][type[j]]/hsec; - stemp = 1.0/(shape[type[i]][0]*cr60+h12)+ - 1.0/(shape[type[i]][1]*cr60+h12)+ - 1.0/(shape[type[i]][2]*cr60+h12)+ - 1.0/(shape[type[j]][0]*cr60+h12)+ - 1.0/(shape[type[j]][1]*cr60+h12)+ - 1.0/(shape[type[j]][2]*cr60+h12); + stemp = 1.0/(shape1[type[i]][0]*cr60+h12)+ + 1.0/(shape1[type[i]][1]*cr60+h12)+ + 1.0/(shape1[type[i]][2]*cr60+h12)+ + 1.0/(shape1[type[j]][0]*cr60+h12)+ + 1.0/(shape1[type[j]][1]*cr60+h12)+ + 1.0/(shape1[type[j]][2]*cr60+h12); hsec = h12+b_alpha*sec; dspr = 7.0/h12-1.0/hsec+stemp; pbsr = b_alpha*sigma[type[i]][type[j]]/hsec; @@ -832,7 +831,6 @@ double PairRESquared::resquared_lj(const int i, const int j, double *ttor, bool calc_torque) { int *type = atom->type; - double **shape = atom->shape; // pair computations for energy, force, torque @@ -875,9 +873,9 @@ double PairRESquared::resquared_lj(const int i, const int j, double lAtwo[3][3][3]; // A1'*S1^2*wi.lA double scorrect[3]; double half_sigma=sigma[type[i]][type[j]] / 2.0; - scorrect[0] = shape[type[i]][0]+half_sigma; - scorrect[1] = shape[type[i]][1]+half_sigma; - scorrect[2] = shape[type[i]][2]+half_sigma; + scorrect[0] = shape1[type[i]][0]+half_sigma; + scorrect[1] = shape1[type[i]][1]+half_sigma; + scorrect[2] = shape1[type[i]][2]+half_sigma; scorrect[0] = scorrect[0] * scorrect[0] / 2.0; scorrect[1] = scorrect[1] * scorrect[1] / 2.0; scorrect[2] = scorrect[2] * scorrect[2] / 2.0; @@ -909,14 +907,14 @@ double PairRESquared::resquared_lj(const int i, const int j, h12p3 = pow(h12,3.0); double sigmap3 = pow(sigma[type[i]][type[j]],3.0); double stemp = h12/2.0; - Ua = (shape[type[i]][0]+stemp)*(shape[type[i]][1]+stemp)* - (shape[type[i]][2]+stemp)*h12p3/8.0; + Ua = (shape1[type[i]][0]+stemp)*(shape1[type[i]][1]+stemp)* + (shape1[type[i]][2]+stemp)*h12p3/8.0; Ua = (1.0+3.0*tprod)*lshape[type[i]]/Ua; Ua = epsilon[type[i]][type[j]]*Ua*sigmap3*solv_f_a; stemp = h12/cr60; - Ur = (shape[type[i]][0]+stemp)*(shape[type[i]][1]+stemp)* - (shape[type[i]][2]+stemp)*h12p3/60.0; + Ur = (shape1[type[i]][0]+stemp)*(shape1[type[i]][1]+stemp)* + (shape1[type[i]][2]+stemp)*h12p3/60.0; Ur = (1.0+b_alpha*tprod)*lshape[type[i]]/Ur; Ur = epsilon[type[i]][type[j]]*Ur*sigmap3*pow(sigh,6.0)*solv_f_r; @@ -931,17 +929,17 @@ double PairRESquared::resquared_lj(const int i, const int j, spr[1] = 0.5*sigma12p3*s[1]; spr[2] = 0.5*sigma12p3*s[2]; - stemp = 1.0/(shape[type[i]][0]*2.0+h12)+ - 1.0/(shape[type[i]][1]*2.0+h12)+ - 1.0/(shape[type[i]][2]*2.0+h12)+ + stemp = 1.0/(shape1[type[i]][0]*2.0+h12)+ + 1.0/(shape1[type[i]][1]*2.0+h12)+ + 1.0/(shape1[type[i]][2]*2.0+h12)+ 3.0/h12; hsec = h12+3.0*sec; dspu = 1.0/h12-1.0/hsec+stemp; pbsu = 3.0*sigma[type[i]][type[j]]/hsec; - stemp = 1.0/(shape[type[i]][0]*cr60+h12)+ - 1.0/(shape[type[i]][1]*cr60+h12)+ - 1.0/(shape[type[i]][2]*cr60+h12)+ + stemp = 1.0/(shape1[type[i]][0]*cr60+h12)+ + 1.0/(shape1[type[i]][1]*cr60+h12)+ + 1.0/(shape1[type[i]][2]*cr60+h12)+ 3.0/h12; hsec = h12+b_alpha*sec; dspr = 7.0/h12-1.0/hsec+stemp; diff --git a/src/ASPHERE/pair_resquared.h b/src/ASPHERE/pair_resquared.h index dae2dd13e2..85605a33f1 100755 --- a/src/ASPHERE/pair_resquared.h +++ b/src/ASPHERE/pair_resquared.h @@ -42,7 +42,8 @@ class PairRESquared : public Pair { double cut_global; double **cut; - double **shape2; // radii in x, y and z SQUARED + double **shape1; // per-type radii in x, y and z + double **shape2; // per-type radii in x, y and z SQUARED double *lshape; // product of the radii double **well; // well depth scaling along each axis double **epsilon,**sigma; // epsilon and sigma values for atom-type pairs diff --git a/src/COLLOID/Install.sh b/src/COLLOID/Install.sh index 36edc72b7e..eb25061d5a 100644 --- a/src/COLLOID/Install.sh +++ b/src/COLLOID/Install.sh @@ -2,13 +2,11 @@ if (test $1 = 1) then - cp atom_vec_colloid.cpp .. cp fix_wall_colloid.cpp .. cp pair_colloid.cpp .. cp pair_lubricate.cpp .. cp pair_yukawa_colloid.cpp .. - cp atom_vec_colloid.h .. cp fix_wall_colloid.h .. cp pair_colloid.h .. cp pair_lubricate.h .. @@ -16,13 +14,11 @@ if (test $1 = 1) then elif (test $1 = 0) then - rm ../atom_vec_colloid.cpp rm ../fix_wall_colloid.cpp rm ../pair_colloid.cpp rm ../pair_lubricate.cpp rm ../pair_yukawa_colloid.cpp - rm ../atom_vec_colloid.h rm ../fix_wall_colloid.h rm ../pair_colloid.h rm ../pair_lubricate.h diff --git a/src/COLLOID/atom_vec_colloid.cpp b/src/COLLOID/atom_vec_colloid.cpp deleted file mode 100644 index 2f71e33472..0000000000 --- a/src/COLLOID/atom_vec_colloid.cpp +++ /dev/null @@ -1,749 +0,0 @@ -/* ---------------------------------------------------------------------- - LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator - http://lammps.sandia.gov, Sandia National Laboratories - Steve Plimpton, sjplimp@sandia.gov - - Copyright (2003) Sandia Corporation. Under the terms of Contract - DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains - certain rights in this software. This software is distributed under - the GNU General Public License. - - See the README file in the top-level LAMMPS directory. -------------------------------------------------------------------------- */ - -#include "lmptype.h" -#include "stdlib.h" -#include "atom_vec_colloid.h" -#include "atom.h" -#include "force.h" -#include "domain.h" -#include "modify.h" -#include "fix.h" -#include "memory.h" -#include "error.h" - -using namespace LAMMPS_NS; - -#define DELTA 10000 - -/* ---------------------------------------------------------------------- */ - -AtomVecColloid::AtomVecColloid(LAMMPS *lmp, int narg, char **arg) : - AtomVec(lmp, narg, arg) -{ - molecular = 0; - mass_type = 1; - shape_type = 1; - - comm_x_only = 1; - comm_f_only = 0; - size_forward = 3; - size_reverse = 6; - size_border = 6; - size_velocity = 6; - size_data_atom = 5; - size_data_vel = 7; - xcol_data = 3; - - atom->omega_flag = atom->torque_flag = 1; -} - -/* ---------------------------------------------------------------------- - grow atom arrays - n = 0 grows arrays by DELTA - n > 0 allocates arrays to size n -------------------------------------------------------------------------- */ - -void AtomVecColloid::grow(int n) -{ - if (n == 0) nmax += DELTA; - else nmax = n; - atom->nmax = nmax; - if (nmax < 0 || nmax > MAXSMALLINT) - error->one("Per-processor system is too big"); - - tag = memory->grow(atom->tag,nmax,"atom:tag"); - type = memory->grow(atom->type,nmax,"atom:type"); - mask = memory->grow(atom->mask,nmax,"atom:mask"); - image = memory->grow(atom->image,nmax,"atom:image"); - x = memory->grow(atom->x,nmax,3,"atom:x"); - v = memory->grow(atom->v,nmax,3,"atom:v"); - f = memory->grow(atom->f,nmax,3,"atom:f"); - - omega = memory->grow(atom->omega,nmax,3,"atom:omega"); - torque = memory->grow(atom->torque,nmax,3,"atom:torque"); - - if (atom->nextra_grow) - for (int iextra = 0; iextra < atom->nextra_grow; iextra++) - modify->fix[atom->extra_grow[iextra]]->grow_arrays(nmax); -} - -/* ---------------------------------------------------------------------- - reset local array ptrs -------------------------------------------------------------------------- */ - -void AtomVecColloid::grow_reset() -{ - tag = atom->tag; type = atom->type; - mask = atom->mask; image = atom->image; - x = atom->x; v = atom->v; f = atom->f; - omega = atom->omega; torque = atom->torque; -} - -/* ---------------------------------------------------------------------- */ - -void AtomVecColloid::copy(int i, int j) -{ - tag[j] = tag[i]; - type[j] = type[i]; - mask[j] = mask[i]; - image[j] = image[i]; - x[j][0] = x[i][0]; - x[j][1] = x[i][1]; - x[j][2] = x[i][2]; - v[j][0] = v[i][0]; - v[j][1] = v[i][1]; - v[j][2] = v[i][2]; - - omega[j][0] = omega[i][0]; - omega[j][1] = omega[i][1]; - omega[j][2] = omega[i][2]; - - if (atom->nextra_grow) - for (int iextra = 0; iextra < atom->nextra_grow; iextra++) - modify->fix[atom->extra_grow[iextra]]->copy_arrays(i,j); -} - -/* ---------------------------------------------------------------------- */ - -int AtomVecColloid::pack_comm(int n, int *list, double *buf, - int pbc_flag, int *pbc) -{ - int i,j,m; - double dx,dy,dz; - - m = 0; - if (pbc_flag == 0) { - for (i = 0; i < n; i++) { - j = list[i]; - buf[m++] = x[j][0]; - buf[m++] = x[j][1]; - buf[m++] = x[j][2]; - } - } else { - if (domain->triclinic == 0) { - dx = pbc[0]*domain->xprd; - dy = pbc[1]*domain->yprd; - dz = pbc[2]*domain->zprd; - } else { - dx = pbc[0]*domain->xprd + pbc[5]*domain->xy + pbc[4]*domain->xz; - 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; - } - } - return m; -} - -/* ---------------------------------------------------------------------- */ - -int AtomVecColloid::pack_comm_vel(int n, int *list, double *buf, - int pbc_flag, int *pbc) -{ - int i,j,m; - double dx,dy,dz,dvx,dvy,dvz; - - m = 0; - if (pbc_flag == 0) { - for (i = 0; i < n; i++) { - j = list[i]; - buf[m++] = x[j][0]; - buf[m++] = x[j][1]; - buf[m++] = x[j][2]; - buf[m++] = v[j][0]; - buf[m++] = v[j][1]; - buf[m++] = v[j][2]; - buf[m++] = omega[j][0]; - buf[m++] = omega[j][1]; - buf[m++] = omega[j][2]; - } - } else { - if (domain->triclinic == 0) { - dx = pbc[0]*domain->xprd; - dy = pbc[1]*domain->yprd; - dz = pbc[2]*domain->zprd; - } else { - dx = pbc[0]*domain->xprd + pbc[5]*domain->xy + pbc[4]*domain->xz; - dy = pbc[1]*domain->yprd + pbc[3]*domain->yz; - dz = pbc[2]*domain->zprd; - } - 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]; - buf[m++] = omega[j][0]; - buf[m++] = omega[j][1]; - buf[m++] = omega[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]; - } - buf[m++] = omega[j][0]; - buf[m++] = omega[j][1]; - buf[m++] = omega[j][2]; - } - } - } - return m; -} - -/* ---------------------------------------------------------------------- */ - -void AtomVecColloid::unpack_comm(int n, int first, double *buf) -{ - int i,m,last; - - m = 0; - last = first + n; - for (i = first; i < last; i++) { - x[i][0] = buf[m++]; - x[i][1] = buf[m++]; - x[i][2] = buf[m++]; - } -} - - -/* ---------------------------------------------------------------------- */ - -void AtomVecColloid::unpack_comm_vel(int n, int first, double *buf) -{ - int i,m,last; - - m = 0; - last = first + n; - for (i = first; i < last; i++) { - x[i][0] = buf[m++]; - x[i][1] = buf[m++]; - x[i][2] = buf[m++]; - v[i][0] = buf[m++]; - v[i][1] = buf[m++]; - v[i][2] = buf[m++]; - omega[i][0] = buf[m++]; - omega[i][1] = buf[m++]; - omega[i][2] = buf[m++]; - } -} - -/* ---------------------------------------------------------------------- */ - -int AtomVecColloid::pack_reverse(int n, int first, double *buf) -{ - int i,m,last; - - m = 0; - last = first + n; - for (i = first; i < last; i++) { - buf[m++] = f[i][0]; - buf[m++] = f[i][1]; - buf[m++] = f[i][2]; - buf[m++] = torque[i][0]; - buf[m++] = torque[i][1]; - buf[m++] = torque[i][2]; - } - return m; -} - -/* ---------------------------------------------------------------------- */ - -int AtomVecColloid::pack_reverse_one(int i, double *buf) -{ - buf[0] = torque[i][0]; - buf[1] = torque[i][1]; - buf[2] = torque[i][2]; - return 3; -} - -/* ---------------------------------------------------------------------- */ - -void AtomVecColloid::unpack_reverse(int n, int *list, double *buf) -{ - int i,j,m; - - m = 0; - for (i = 0; i < n; i++) { - j = list[i]; - f[j][0] += buf[m++]; - f[j][1] += buf[m++]; - f[j][2] += buf[m++]; - torque[j][0] += buf[m++]; - torque[j][1] += buf[m++]; - torque[j][2] += buf[m++]; - } -} - -/* ---------------------------------------------------------------------- */ - -int AtomVecColloid::unpack_reverse_one(int i, double *buf) -{ - torque[i][0] += buf[0]; - torque[i][1] += buf[1]; - torque[i][2] += buf[2]; - return 3; -} - -/* ---------------------------------------------------------------------- */ - -int AtomVecColloid::pack_border(int n, int *list, double *buf, - int pbc_flag, int *pbc) -{ - int i,j,m; - double dx,dy,dz; - - m = 0; - if (pbc_flag == 0) { - for (i = 0; i < n; i++) { - j = list[i]; - buf[m++] = x[j][0]; - buf[m++] = x[j][1]; - buf[m++] = x[j][2]; - buf[m++] = tag[j]; - buf[m++] = type[j]; - buf[m++] = mask[j]; - } - } else { - if (domain->triclinic == 0) { - dx = pbc[0]*domain->xprd; - dy = pbc[1]*domain->yprd; - dz = pbc[2]*domain->zprd; - } else { - dx = pbc[0]; - 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]; - } - } - return m; -} - -/* ---------------------------------------------------------------------- */ - -int AtomVecColloid::pack_border_vel(int n, int *list, double *buf, - int pbc_flag, int *pbc) -{ - int i,j,m; - double dx,dy,dz,dvx,dvy,dvz; - - m = 0; - if (pbc_flag == 0) { - for (i = 0; i < n; i++) { - j = list[i]; - buf[m++] = x[j][0]; - buf[m++] = x[j][1]; - buf[m++] = x[j][2]; - buf[m++] = tag[j]; - buf[m++] = type[j]; - buf[m++] = mask[j]; - buf[m++] = v[j][0]; - buf[m++] = v[j][1]; - buf[m++] = v[j][2]; - buf[m++] = omega[j][0]; - buf[m++] = omega[j][1]; - buf[m++] = omega[j][2]; - } - } else { - if (domain->triclinic == 0) { - dx = pbc[0]*domain->xprd; - dy = pbc[1]*domain->yprd; - dz = pbc[2]*domain->zprd; - } else { - dx = pbc[0]; - dy = pbc[1]; - dz = pbc[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++] = v[j][0]; - buf[m++] = v[j][1]; - buf[m++] = v[j][2]; - buf[m++] = omega[j][0]; - buf[m++] = omega[j][1]; - buf[m++] = omega[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]; - 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]; - } - buf[m++] = omega[j][0]; - buf[m++] = omega[j][1]; - buf[m++] = omega[j][2]; - } - } - } - return m; -} - -/* ---------------------------------------------------------------------- */ - -void AtomVecColloid::unpack_border(int n, int first, double *buf) -{ - int i,m,last; - - m = 0; - last = first + n; - for (i = first; i < last; i++) { - if (i == nmax) grow(0); - x[i][0] = buf[m++]; - x[i][1] = buf[m++]; - x[i][2] = buf[m++]; - tag[i] = static_cast (buf[m++]); - type[i] = static_cast (buf[m++]); - mask[i] = static_cast (buf[m++]); - } -} - -/* ---------------------------------------------------------------------- */ - -void AtomVecColloid::unpack_border_vel(int n, int first, double *buf) -{ - int i,m,last; - - m = 0; - last = first + n; - for (i = first; i < last; i++) { - if (i == nmax) grow(0); - x[i][0] = buf[m++]; - x[i][1] = buf[m++]; - x[i][2] = buf[m++]; - tag[i] = static_cast (buf[m++]); - type[i] = static_cast (buf[m++]); - mask[i] = static_cast (buf[m++]); - v[i][0] = buf[m++]; - v[i][1] = buf[m++]; - v[i][2] = buf[m++]; - omega[i][0] = buf[m++]; - omega[i][1] = buf[m++]; - omega[i][2] = buf[m++]; - } -} - -/* ---------------------------------------------------------------------- - pack data for atom I for sending to another proc - xyz must be 1st 3 values, so comm::exchange() can test on them -------------------------------------------------------------------------- */ - -int AtomVecColloid::pack_exchange(int i, double *buf) -{ - int m = 1; - buf[m++] = x[i][0]; - buf[m++] = x[i][1]; - buf[m++] = x[i][2]; - buf[m++] = v[i][0]; - buf[m++] = v[i][1]; - buf[m++] = v[i][2]; - buf[m++] = tag[i]; - buf[m++] = type[i]; - buf[m++] = mask[i]; - buf[m++] = image[i]; - - buf[m++] = omega[i][0]; - buf[m++] = omega[i][1]; - buf[m++] = omega[i][2]; - - if (atom->nextra_grow) - for (int iextra = 0; iextra < atom->nextra_grow; iextra++) - m += modify->fix[atom->extra_grow[iextra]]->pack_exchange(i,&buf[m]); - - buf[0] = m; - return m; -} - -/* ---------------------------------------------------------------------- */ - -int AtomVecColloid::unpack_exchange(double *buf) -{ - int nlocal = atom->nlocal; - if (nlocal == nmax) grow(0); - - int m = 1; - x[nlocal][0] = buf[m++]; - x[nlocal][1] = buf[m++]; - x[nlocal][2] = buf[m++]; - v[nlocal][0] = buf[m++]; - v[nlocal][1] = buf[m++]; - v[nlocal][2] = buf[m++]; - tag[nlocal] = static_cast (buf[m++]); - type[nlocal] = static_cast (buf[m++]); - mask[nlocal] = static_cast (buf[m++]); - image[nlocal] = static_cast (buf[m++]); - - omega[nlocal][0] = buf[m++]; - omega[nlocal][1] = buf[m++]; - omega[nlocal][2] = buf[m++]; - - if (atom->nextra_grow) - for (int iextra = 0; iextra < atom->nextra_grow; iextra++) - m += modify->fix[atom->extra_grow[iextra]]-> - unpack_exchange(nlocal,&buf[m]); - - atom->nlocal++; - return m; -} - -/* ---------------------------------------------------------------------- - size of restart data for all atoms owned by this proc - include extra data stored by fixes -------------------------------------------------------------------------- */ - -int AtomVecColloid::size_restart() -{ - int i; - - int nlocal = atom->nlocal; - int n = 14 * nlocal; - - if (atom->nextra_restart) - for (int iextra = 0; iextra < atom->nextra_restart; iextra++) - for (i = 0; i < nlocal; i++) - n += modify->fix[atom->extra_restart[iextra]]->size_restart(i); - - return n; -} - -/* ---------------------------------------------------------------------- - pack atom I's data for restart file including extra quantities - xyz must be 1st 3 values, so that read_restart can test on them - molecular types may be negative, but write as positive -------------------------------------------------------------------------- */ - -int AtomVecColloid::pack_restart(int i, double *buf) -{ - int m = 1; - buf[m++] = x[i][0]; - buf[m++] = x[i][1]; - buf[m++] = x[i][2]; - buf[m++] = tag[i]; - buf[m++] = type[i]; - buf[m++] = mask[i]; - buf[m++] = image[i]; - buf[m++] = v[i][0]; - buf[m++] = v[i][1]; - buf[m++] = v[i][2]; - - buf[m++] = omega[i][0]; - buf[m++] = omega[i][1]; - buf[m++] = omega[i][2]; - - 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]); - - buf[0] = m; - return m; -} - -/* ---------------------------------------------------------------------- - unpack data for one atom from restart file including extra quantities -------------------------------------------------------------------------- */ - -int AtomVecColloid::unpack_restart(double *buf) -{ - int nlocal = atom->nlocal; - if (nlocal == nmax) { - grow(0); - if (atom->nextra_store) - memory->grow(atom->extra,nmax,atom->nextra_store,"atom:extra"); - } - - int m = 1; - x[nlocal][0] = buf[m++]; - x[nlocal][1] = buf[m++]; - x[nlocal][2] = buf[m++]; - tag[nlocal] = static_cast (buf[m++]); - type[nlocal] = static_cast (buf[m++]); - mask[nlocal] = static_cast (buf[m++]); - image[nlocal] = static_cast (buf[m++]); - v[nlocal][0] = buf[m++]; - v[nlocal][1] = buf[m++]; - v[nlocal][2] = buf[m++]; - - omega[nlocal][0] = buf[m++]; - omega[nlocal][1] = buf[m++]; - omega[nlocal][2] = buf[m++]; - - double **extra = atom->extra; - if (atom->nextra_store) { - int size = static_cast (buf[0]) - m; - for (int i = 0; i < size; i++) extra[nlocal][i] = buf[m++]; - } - - atom->nlocal++; - return m; -} - -/* ---------------------------------------------------------------------- - create one atom of itype at coord - set other values to defaults -------------------------------------------------------------------------- */ - -void AtomVecColloid::create_atom(int itype, double *coord) -{ - int nlocal = atom->nlocal; - if (nlocal == nmax) grow(0); - - tag[nlocal] = 0; - type[nlocal] = itype; - x[nlocal][0] = coord[0]; - x[nlocal][1] = coord[1]; - x[nlocal][2] = coord[2]; - mask[nlocal] = 1; - image[nlocal] = (512 << 20) | (512 << 10) | 512; - v[nlocal][0] = 0.0; - v[nlocal][1] = 0.0; - v[nlocal][2] = 0.0; - - omega[nlocal][0] = 0.0; - omega[nlocal][1] = 0.0; - omega[nlocal][2] = 0.0; - - atom->nlocal++; -} - -/* ---------------------------------------------------------------------- - unpack one line from Atoms section of data file - initialize other atom quantities -------------------------------------------------------------------------- */ - -void AtomVecColloid::data_atom(double *coord, int imagetmp, char **values) -{ - int nlocal = atom->nlocal; - if (nlocal == nmax) grow(0); - - tag[nlocal] = atoi(values[0]); - if (tag[nlocal] <= 0) - error->one("Invalid atom ID in Atoms section of data file"); - - type[nlocal] = atoi(values[1]); - if (type[nlocal] <= 0 || type[nlocal] > atom->ntypes) - error->one("Invalid atom type in Atoms section of data file"); - - x[nlocal][0] = coord[0]; - x[nlocal][1] = coord[1]; - x[nlocal][2] = coord[2]; - - image[nlocal] = imagetmp; - - mask[nlocal] = 1; - v[nlocal][0] = 0.0; - v[nlocal][1] = 0.0; - v[nlocal][2] = 0.0; - omega[nlocal][0] = 0.0; - omega[nlocal][1] = 0.0; - omega[nlocal][2] = 0.0; - - atom->nlocal++; -} - -/* ---------------------------------------------------------------------- - unpack one line from Velocities section of data file -------------------------------------------------------------------------- */ - -void AtomVecColloid::data_vel(int m, char **values) -{ - v[m][0] = atof(values[0]); - v[m][1] = atof(values[1]); - v[m][2] = atof(values[2]); - omega[m][0] = atof(values[3]); - omega[m][1] = atof(values[4]); - omega[m][2] = atof(values[5]); -} - -/* ---------------------------------------------------------------------- - unpack hybrid quantities from one line in Velocities section of data file -------------------------------------------------------------------------- */ - -int AtomVecColloid::data_vel_hybrid(int m, char **values) -{ - omega[m][0] = atof(values[0]); - omega[m][1] = atof(values[1]); - omega[m][2] = atof(values[2]); - return 3; -} - -/* ---------------------------------------------------------------------- - return # of bytes of allocated memory -------------------------------------------------------------------------- */ - -bigint AtomVecColloid::memory_usage() -{ - bigint bytes = 0; - - if (atom->memcheck("tag")) bytes += memory->usage(tag,nmax); - if (atom->memcheck("type")) bytes += memory->usage(type,nmax); - if (atom->memcheck("mask")) bytes += memory->usage(mask,nmax); - if (atom->memcheck("image")) bytes += memory->usage(image,nmax); - if (atom->memcheck("x")) bytes += memory->usage(x,nmax,3); - if (atom->memcheck("v")) bytes += memory->usage(v,nmax,3); - if (atom->memcheck("f")) bytes += memory->usage(f,nmax,3); - - if (atom->memcheck("omega")) bytes += memory->usage(omega,nmax,3); - if (atom->memcheck("torque")) bytes += memory->usage(torque,nmax,3); - - return bytes; -} diff --git a/src/COLLOID/atom_vec_colloid.h b/src/COLLOID/atom_vec_colloid.h deleted file mode 100644 index 121cbc82d3..0000000000 --- a/src/COLLOID/atom_vec_colloid.h +++ /dev/null @@ -1,66 +0,0 @@ -/* ---------------------------------------------------------------------- - LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator - http://lammps.sandia.gov, Sandia National Laboratories - Steve Plimpton, sjplimp@sandia.gov - - Copyright (2003) Sandia Corporation. Under the terms of Contract - DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains - certain rights in this software. This software is distributed under - the GNU General Public License. - - See the README file in the top-level LAMMPS directory. -------------------------------------------------------------------------- */ - -#ifdef ATOM_CLASS - -AtomStyle(colloid,AtomVecColloid) - -#else - -#ifndef LMP_ATOM_VEC_COLLOID_H -#define LMP_ATOM_VEC_COLLOID_H - -#include "atom_vec.h" - -namespace LAMMPS_NS { - -class AtomVecColloid : public AtomVec { - public: - AtomVecColloid(class LAMMPS *, int, char **); - virtual ~AtomVecColloid() {} - void grow(int); - void grow_reset(); - void copy(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 *); - void unpack_comm_vel(int, int, double *); - int pack_reverse(int, int, double *); - int pack_reverse_one(int, double *); - void unpack_reverse(int, int *, double *); - int unpack_reverse_one(int, double *); - int pack_border(int, int *, double *, int, int *); - int pack_border_vel(int, int *, double *, int, int *); - void unpack_border(int, int, double *); - void unpack_border_vel(int, int, double *); - int pack_exchange(int, double *); - int unpack_exchange(double *); - int size_restart(); - int pack_restart(int, double *); - int unpack_restart(double *); - void create_atom(int, double *); - void data_atom(double *, int, char **); - void data_vel(int, char **); - int data_vel_hybrid(int, char **); - bigint memory_usage(); - - private: - int *tag,*type,*mask,*image; - double **x,**v,**f; - double **omega,**torque; -}; - -} - -#endif -#endif diff --git a/src/COLLOID/fix_wall_colloid.cpp b/src/COLLOID/fix_wall_colloid.cpp index 37a1cda850..0c03ad57a7 100644 --- a/src/COLLOID/fix_wall_colloid.cpp +++ b/src/COLLOID/fix_wall_colloid.cpp @@ -35,31 +35,19 @@ FixWallColloid::FixWallColloid(LAMMPS *lmp, int narg, char **arg) : void FixWallColloid::init() { - if (!atom->avec->shape_type) - error->all("Fix wall/colloid requires atom attribute shape"); - if (atom->radius_flag) - error->all("Fix wall/colloid cannot be used with atom attribute diameter"); - - // insure all particle shapes are spherical - // can be polydisperse - - for (int i = 1; i <= atom->ntypes; i++) - if ((atom->shape[i][0] != atom->shape[i][1]) || - (atom->shape[i][0] != atom->shape[i][2]) || - (atom->shape[i][1] != atom->shape[i][2])) - error->all("Fix wall/colloid requires spherical particles"); + if (!atom->sphere_flag) + error->all("Fix wall/colloid requires atom style sphere"); // insure all particles in group are extended particles - double **shape = atom->shape; - int *type = atom->type; + double *radius = atom->radius; int *mask = atom->mask; int nlocal = atom->nlocal; int flag = 0; for (int i = 0; i < nlocal; i++) if (mask[i] & groupbit) - if (atom->shape[type[i]][0] == 0.0) flag = 1; + if (radius[i] == 0.0) flag = 1; int flagall; MPI_Allreduce(&flag,&flagall,1,MPI_INT,MPI_SUM,world); @@ -95,9 +83,8 @@ void FixWallColloid::wall_particle(int m, int which, double coord) double **x = atom->x; double **f = atom->f; - double **shape = atom->shape; + double *radius = atom->radius; int *mask = atom->mask; - int *type = atom->type; int nlocal = atom->nlocal; int dim = which / 2; @@ -111,7 +98,7 @@ void FixWallColloid::wall_particle(int m, int which, double coord) if (side < 0) delta = x[i][dim] - coord; else delta = coord - x[i][dim]; if (delta >= cutoff[m]) continue; - rad = shape[type[i]][0]; + rad = radius[i]; if (rad >= delta) { onflag = 1; continue; diff --git a/src/COLLOID/pair_colloid.cpp b/src/COLLOID/pair_colloid.cpp index e79891307c..9746752a71 100644 --- a/src/COLLOID/pair_colloid.cpp +++ b/src/COLLOID/pair_colloid.cpp @@ -23,6 +23,7 @@ #include "atom.h" #include "comm.h" #include "force.h" +#include "neighbor.h" #include "neigh_list.h" #include "memory.h" #include "error.h" diff --git a/src/COLLOID/pair_lubricate.cpp b/src/COLLOID/pair_lubricate.cpp index 8d86029221..6d8dec2690 100644 --- a/src/COLLOID/pair_lubricate.cpp +++ b/src/COLLOID/pair_lubricate.cpp @@ -85,9 +85,9 @@ void PairLubricate::compute(int eflag, int vflag) double **v = atom->v; double **f = atom->f; double **omega = atom->omega; - double **angmom = atom->angmom; double **torque = atom->torque; - double **shape = atom->shape; + double *radius = atom->radius; + double *rmass = atom->rmass; int *type = atom->type; int nlocal = atom->nlocal; int newton_pair = force->newton_pair; @@ -111,7 +111,7 @@ void PairLubricate::compute(int eflag, int vflag) ytmp = x[i][1]; ztmp = x[i][2]; itype = type[i]; - radi = shape[itype][0]; + radi = radius[i]; jlist = firstneigh[i]; jnum = numneigh[i]; @@ -151,19 +151,10 @@ void PairLubricate::compute(int eflag, int vflag) vt3 = vr3 - vn3; // additive rotational velocity = omega_1 + omega_2 - // use omega directly if it exists, else angmom - // angular momentum = I*omega = 2/5 * M*R^2 * omega - if (omega_flag) { - w1 = omega[i][0] + omega[j][0]; - w2 = omega[i][1] + omega[j][1]; - w3 = omega[i][2] + omega[j][2]; - } else { - inv_inertia = 1.0 / (0.4*atom->mass[itype]*radi*radi); - w1 = inv_inertia * (angmom[i][0] + angmom[j][0]); - w2 = inv_inertia * (angmom[i][1] + angmom[j][1]); - w3 = inv_inertia * (angmom[i][2] + angmom[j][2]); - } + w1 = omega[i][0] + omega[j][0]; + w2 = omega[i][1] + omega[j][1]; + w3 = omega[i][2] + omega[j][2]; // relative velocities n X P . (v1-v2) = n X (I-nn) . (v1-v2) @@ -185,15 +176,9 @@ void PairLubricate::compute(int eflag, int vflag) // N.(w1-w2) and P.(w1-w2) - if (omega_flag) { - wr1 = omega[i][0] - omega[j][0]; - wr2 = omega[i][1] - omega[j][1]; - wr3 = omega[i][2] - omega[j][2]; - } else { - wr1 = inv_inertia * (angmom[i][0] - angmom[j][0]); - wr2 = inv_inertia * (angmom[i][1] - angmom[j][1]); - wr3 = inv_inertia * (angmom[i][2] - angmom[j][2]); - } + wr1 = omega[i][0] - omega[j][0]; + wr2 = omega[i][1] - omega[j][1]; + wr3 = omega[i][2] - omega[j][2]; wnnr = wr1*delx + wr2*dely + wr3*delz; wn1 = delx*wnnr / rsq; @@ -393,29 +378,24 @@ void PairLubricate::coeff(int narg, char **arg) void PairLubricate::init_style() { + if (!atom->sphere_flag) + error->all("Pair lubricate requires atom style sphere"); if (comm->ghost_velocity == 0) error->all("Pair lubricate requires ghost atoms store velocity"); - if (!atom->torque_flag || !atom->avec->shape_type) - error->all("Pair lubricate requires atom attributes torque and shape"); - if (atom->radius_flag || atom->rmass_flag) - error->all("Pair lubricate cannot be used with atom attributes " - "diameter or rmass"); - if (atom->omega_flag) omega_flag = 1; - else if (atom->angmom_flag) omega_flag = 0; - else error->all("Pair lubricate requires atom attribute omega or angmom"); - - // insure all particle shapes are finite-size, spherical, and monodisperse - // for pair hybrid, should limit test to types using the pair style - - double radi = atom->shape[1][0]; - if (radi == 0.0) error->all("Pair lubricate requires extended particles"); - for (int i = 1; i <= atom->ntypes; i++) - if (atom->shape[i][0] != radi || atom->shape[i][0] != radi || - atom->shape[i][0] != radi) - error->all("Pair lubricate requires spherical, mono-disperse particles"); - int irequest = neighbor->request(this); + + // require that atom radii are identical within each type require + // monodisperse system with same radii for all types + + double rad,radtype; + for (int i = 1; i <= atom->ntypes; i++) { + if (!atom->radius_consistency(i,radtype)) + error->all("Pair colloid requires atoms with same type have same radius"); + if (i > 1 && radtype != rad) + error->all("Pair colloid requires mono-disperse particles"); + rad = radtype; + } } /* ---------------------------------------------------------------------- diff --git a/src/COLLOID/pair_lubricate.h b/src/COLLOID/pair_lubricate.h index 5f6f1ef6c6..30c7bee63c 100644 --- a/src/COLLOID/pair_lubricate.h +++ b/src/COLLOID/pair_lubricate.h @@ -43,7 +43,7 @@ class PairLubricate : public Pair { double cut_inner_global,cut_global; double t_target,mu; int flag1,flag2,flag3,flag4; - int seed,omega_flag; + int seed; double **cut_inner,**cut; class RanMars *random; diff --git a/src/COLLOID/pair_yukawa_colloid.cpp b/src/COLLOID/pair_yukawa_colloid.cpp index 69c974a584..46f66624ee 100644 --- a/src/COLLOID/pair_yukawa_colloid.cpp +++ b/src/COLLOID/pair_yukawa_colloid.cpp @@ -48,7 +48,7 @@ void PairYukawaColloid::compute(int eflag, int vflag) double **x = atom->x; double **f = atom->f; - double **shape = atom->shape; + double *radius = atom->radius; int *type = atom->type; int nlocal = atom->nlocal; double *special_coul = force->special_coul; @@ -67,7 +67,7 @@ void PairYukawaColloid::compute(int eflag, int vflag) ytmp = x[i][1]; ztmp = x[i][2]; itype = type[i]; - radi = shape[itype][0]; + radi = radius[i]; jlist = firstneigh[i]; jnum = numneigh[i]; @@ -81,7 +81,7 @@ void PairYukawaColloid::compute(int eflag, int vflag) delz = ztmp - x[j][2]; rsq = delx*delx + dely*dely + delz*delz; jtype = type[j]; - radj = shape[jtype][0]; + radj = radius[j]; if (rsq < cutsq[itype][jtype]) { r2inv = 1.0/rsq; @@ -121,22 +121,17 @@ void PairYukawaColloid::compute(int eflag, int vflag) void PairYukawaColloid::init_style() { - if (!atom->avec->shape_type) - error->all("Pair yukawa/colloid requires atom attribute shape"); - if (atom->radius_flag) - error->all("Pair yukawa/colloid cannot be used with " - "atom attribute diameter"); + if (!atom->sphere_flag) + error->all("Pair yukawa/colloid requires atom style sphere"); - // insure all particle shapes are spherical - // can be point particles or polydisperse + int irequest = neighbor->request(this); + + // require that atom radii are identical within each type for (int i = 1; i <= atom->ntypes; i++) - if ((atom->shape[i][0] != atom->shape[i][1]) || - (atom->shape[i][0] != atom->shape[i][2]) || - (atom->shape[i][1] != atom->shape[i][2])) - error->all("Pair yukawa/colloid requires spherical particles"); - - int irequest = neighbor->request(this); + if (!atom->radius_consistency(i,rad[i])) + error->all("Pair yukawa/colloid requires atoms with same type " + "have same radius"); } /* ---------------------------------------------------------------------- @@ -151,9 +146,7 @@ double PairYukawaColloid::init_one(int i, int j) } if (offset_flag) { - double radi = atom->shape[i][0]; - double radj = atom->shape[j][0]; - double screening = exp(-kappa * (cut[i][j] - (radi+radj))); + double screening = exp(-kappa * (cut[i][j] - (rad[i]+rad[j]))); offset[i][j] = a[i][j]/kappa * screening; } else offset[i][j] = 0.0; @@ -170,16 +163,12 @@ double PairYukawaColloid::single(int i, int j, int itype, int jtype, double factor_coul, double factor_lj, double &fforce) { - double r2inv,r,rinv,screening,forceyukawa,phi,radi,radj; - - int *type = atom->type; - radi = atom->shape[itype][0]; - radj = atom->shape[jtype][0]; - + double r2inv,r,rinv,screening,forceyukawa,phi; + r2inv = 1.0/rsq; r = sqrt(rsq); rinv = 1.0/r; - screening = exp(-kappa*(r-(radi+radj))); + screening = exp(-kappa*(r-(rad[itype]+rad[jtype]))); forceyukawa = a[itype][jtype] * screening; fforce = factor_coul*forceyukawa * rinv; diff --git a/src/DIPOLE/atom_vec_dipole.cpp b/src/DIPOLE/atom_vec_dipole.cpp index d033175301..7a1bc6be14 100644 --- a/src/DIPOLE/atom_vec_dipole.cpp +++ b/src/DIPOLE/atom_vec_dipole.cpp @@ -12,6 +12,7 @@ ------------------------------------------------------------------------- */ #include "lmptype.h" +#include "math.h" #include "stdlib.h" #include "atom_vec_dipole.h" #include "atom.h" @@ -32,19 +33,19 @@ AtomVecDipole::AtomVecDipole(LAMMPS *lmp, int narg, char **arg) : { molecular = 0; mass_type = 1; - shape_type = 1; - dipole_type = 1; - comm_x_only = comm_f_only = 0; + comm_x_only = 0; + comm_f_only = 1; size_forward = 6; - size_reverse = 6; - size_border = 10; - size_velocity = 6; + size_reverse = 3; + size_border = 11; + size_velocity = 3; size_data_atom = 9; size_data_vel = 7; xcol_data = 4; - atom->q_flag = atom->mu_flag = atom->omega_flag = atom->torque_flag = 1; + atom->dipole_flag = 1; + atom->q_flag = atom->mu_flag = 1; } /* ---------------------------------------------------------------------- @@ -70,9 +71,7 @@ void AtomVecDipole::grow(int n) f = memory->grow(atom->f,nmax,3,"atom:f"); q = memory->grow(atom->q,nmax,"atom:q"); - mu = memory->grow(atom->mu,nmax,3,"atom:mu"); - omega = memory->grow(atom->omega,nmax,3,"atom:omega"); - torque = memory->grow(atom->torque,nmax,3,"atom:torque"); + mu = memory->grow(atom->mu,nmax,4,"atom:mu"); if (atom->nextra_grow) for (int iextra = 0; iextra < atom->nextra_grow; iextra++) @@ -89,7 +88,6 @@ void AtomVecDipole::grow_reset() mask = atom->mask; image = atom->image; x = atom->x; v = atom->v; f = atom->f; q = atom->q; mu = atom->mu; - omega = atom->omega; torque = atom->torque; } /* ---------------------------------------------------------------------- */ @@ -111,9 +109,7 @@ void AtomVecDipole::copy(int i, int j) mu[j][0] = mu[i][0]; mu[j][1] = mu[i][1]; mu[j][2] = mu[i][2]; - omega[j][0] = omega[i][0]; - omega[j][1] = omega[i][1]; - omega[j][2] = omega[i][2]; + mu[j][3] = mu[i][3]; if (atom->nextra_grow) for (int iextra = 0; iextra < atom->nextra_grow; iextra++) @@ -183,9 +179,6 @@ int AtomVecDipole::pack_comm_vel(int n, int *list, double *buf, buf[m++] = v[j][0]; buf[m++] = v[j][1]; buf[m++] = v[j][2]; - buf[m++] = omega[j][0]; - buf[m++] = omega[j][1]; - buf[m++] = omega[j][2]; } } else { if (domain->triclinic == 0) { @@ -209,9 +202,6 @@ int AtomVecDipole::pack_comm_vel(int n, int *list, double *buf, buf[m++] = v[j][0]; buf[m++] = v[j][1]; buf[m++] = v[j][2]; - buf[m++] = omega[j][0]; - buf[m++] = omega[j][1]; - buf[m++] = omega[j][2]; } } else { dvx = pbc[0]*h_rate[0] + pbc[5]*h_rate[5] + pbc[4]*h_rate[4]; @@ -234,9 +224,6 @@ int AtomVecDipole::pack_comm_vel(int n, int *list, double *buf, buf[m++] = v[j][1]; buf[m++] = v[j][2]; } - buf[m++] = omega[j][0]; - buf[m++] = omega[j][1]; - buf[m++] = omega[j][2]; } } } @@ -289,9 +276,6 @@ void AtomVecDipole::unpack_comm_vel(int n, int first, double *buf) v[i][0] = buf[m++]; v[i][1] = buf[m++]; v[i][2] = buf[m++]; - omega[i][0] = buf[m++]; - omega[i][1] = buf[m++]; - omega[i][2] = buf[m++]; } } @@ -317,25 +301,12 @@ int AtomVecDipole::pack_reverse(int n, int first, double *buf) buf[m++] = f[i][0]; buf[m++] = f[i][1]; buf[m++] = f[i][2]; - buf[m++] = torque[i][0]; - buf[m++] = torque[i][1]; - buf[m++] = torque[i][2]; } return m; } /* ---------------------------------------------------------------------- */ -int AtomVecDipole::pack_reverse_one(int i, double *buf) -{ - buf[0] = torque[i][0]; - buf[1] = torque[i][1]; - buf[2] = torque[i][2]; - return 3; -} - -/* ---------------------------------------------------------------------- */ - void AtomVecDipole::unpack_reverse(int n, int *list, double *buf) { int i,j,m; @@ -346,24 +317,11 @@ void AtomVecDipole::unpack_reverse(int n, int *list, double *buf) f[j][0] += buf[m++]; f[j][1] += buf[m++]; f[j][2] += buf[m++]; - torque[j][0] += buf[m++]; - torque[j][1] += buf[m++]; - torque[j][2] += buf[m++]; } } /* ---------------------------------------------------------------------- */ -int AtomVecDipole::unpack_reverse_one(int i, double *buf) -{ - torque[i][0] += buf[0]; - torque[i][1] += buf[1]; - torque[i][2] += buf[2]; - return 3; -} - -/* ---------------------------------------------------------------------- */ - int AtomVecDipole::pack_border(int n, int *list, double *buf, int pbc_flag, int *pbc) { @@ -384,6 +342,7 @@ int AtomVecDipole::pack_border(int n, int *list, double *buf, buf[m++] = mu[j][0]; buf[m++] = mu[j][1]; buf[m++] = mu[j][2]; + buf[m++] = mu[j][3]; } } else { if (domain->triclinic == 0) { @@ -407,6 +366,7 @@ int AtomVecDipole::pack_border(int n, int *list, double *buf, buf[m++] = mu[j][0]; buf[m++] = mu[j][1]; buf[m++] = mu[j][2]; + buf[m++] = mu[j][3]; } } return m; @@ -434,12 +394,10 @@ int AtomVecDipole::pack_border_vel(int n, int *list, double *buf, buf[m++] = mu[j][0]; buf[m++] = mu[j][1]; buf[m++] = mu[j][2]; + buf[m++] = mu[j][3]; buf[m++] = v[j][0]; buf[m++] = v[j][1]; buf[m++] = v[j][2]; - buf[m++] = omega[j][0]; - buf[m++] = omega[j][1]; - buf[m++] = omega[j][2]; } } else { if (domain->triclinic == 0) { @@ -464,12 +422,10 @@ int AtomVecDipole::pack_border_vel(int n, int *list, double *buf, buf[m++] = mu[j][0]; buf[m++] = mu[j][1]; buf[m++] = mu[j][2]; + buf[m++] = mu[j][3]; buf[m++] = v[j][0]; buf[m++] = v[j][1]; buf[m++] = v[j][2]; - buf[m++] = omega[j][0]; - buf[m++] = omega[j][1]; - buf[m++] = omega[j][2]; } } else { dvx = pbc[0]*h_rate[0] + pbc[5]*h_rate[5] + pbc[4]*h_rate[4]; @@ -487,6 +443,7 @@ int AtomVecDipole::pack_border_vel(int n, int *list, double *buf, buf[m++] = mu[j][0]; buf[m++] = mu[j][1]; buf[m++] = mu[j][2]; + buf[m++] = mu[j][3]; if (mask[i] & deform_groupbit) { buf[m++] = v[j][0] + dvx; buf[m++] = v[j][1] + dvy; @@ -496,9 +453,6 @@ int AtomVecDipole::pack_border_vel(int n, int *list, double *buf, buf[m++] = v[j][1]; buf[m++] = v[j][2]; } - buf[m++] = omega[j][0]; - buf[m++] = omega[j][1]; - buf[m++] = omega[j][2]; } } } @@ -513,7 +467,8 @@ int AtomVecDipole::pack_border_one(int i, double *buf) buf[1] = mu[i][0]; buf[2] = mu[i][1]; buf[3] = mu[i][2]; - return 4; + buf[4] = mu[i][3]; + return 5; } /* ---------------------------------------------------------------------- */ @@ -536,6 +491,7 @@ void AtomVecDipole::unpack_border(int n, int first, double *buf) mu[i][0] = buf[m++]; mu[i][1] = buf[m++]; mu[i][2] = buf[m++]; + mu[i][3] = buf[m++]; } } @@ -559,12 +515,10 @@ void AtomVecDipole::unpack_border_vel(int n, int first, double *buf) mu[i][0] = buf[m++]; mu[i][1] = buf[m++]; mu[i][2] = buf[m++]; + mu[i][3] = buf[m++]; v[i][0] = buf[m++]; v[i][1] = buf[m++]; v[i][2] = buf[m++]; - omega[i][0] = buf[m++]; - omega[i][1] = buf[m++]; - omega[i][2] = buf[m++]; } } @@ -576,7 +530,8 @@ int AtomVecDipole::unpack_border_one(int i, double *buf) mu[i][0] = buf[1]; mu[i][1] = buf[2]; mu[i][2] = buf[3]; - return 4; + mu[i][3] = buf[4]; + return 5; } /* ---------------------------------------------------------------------- @@ -602,9 +557,7 @@ int AtomVecDipole::pack_exchange(int i, double *buf) buf[m++] = mu[i][0]; buf[m++] = mu[i][1]; buf[m++] = mu[i][2]; - buf[m++] = omega[i][0]; - buf[m++] = omega[i][1]; - buf[m++] = omega[i][2]; + buf[m++] = mu[i][3]; if (atom->nextra_grow) for (int iextra = 0; iextra < atom->nextra_grow; iextra++) @@ -637,9 +590,7 @@ int AtomVecDipole::unpack_exchange(double *buf) mu[nlocal][0] = buf[m++]; mu[nlocal][1] = buf[m++]; mu[nlocal][2] = buf[m++]; - omega[nlocal][0] = buf[m++]; - omega[nlocal][1] = buf[m++]; - omega[nlocal][2] = buf[m++]; + mu[nlocal][3] = buf[m++]; if (atom->nextra_grow) for (int iextra = 0; iextra < atom->nextra_grow; iextra++) @@ -660,7 +611,7 @@ int AtomVecDipole::size_restart() int i; int nlocal = atom->nlocal; - int n = 18 * nlocal; + int n = 15 * nlocal; if (atom->nextra_restart) for (int iextra = 0; iextra < atom->nextra_restart; iextra++) @@ -694,9 +645,6 @@ int AtomVecDipole::pack_restart(int i, double *buf) buf[m++] = mu[i][0]; buf[m++] = mu[i][1]; buf[m++] = mu[i][2]; - buf[m++] = omega[i][0]; - buf[m++] = omega[i][1]; - buf[m++] = omega[i][2]; if (atom->nextra_restart) for (int iextra = 0; iextra < atom->nextra_restart; iextra++) @@ -735,9 +683,6 @@ int AtomVecDipole::unpack_restart(double *buf) mu[nlocal][0] = buf[m++]; mu[nlocal][1] = buf[m++]; mu[nlocal][2] = buf[m++]; - omega[nlocal][0] = buf[m++]; - omega[nlocal][1] = buf[m++]; - omega[nlocal][2] = buf[m++]; double **extra = atom->extra; if (atom->nextra_store) { @@ -774,9 +719,7 @@ void AtomVecDipole::create_atom(int itype, double *coord) mu[nlocal][0] = 0.0; mu[nlocal][1] = 0.0; mu[nlocal][2] = 0.0; - omega[nlocal][0] = 0.0; - omega[nlocal][1] = 0.0; - omega[nlocal][2] = 0.0; + mu[nlocal][3] = 0.0; atom->nlocal++; } @@ -808,6 +751,9 @@ void AtomVecDipole::data_atom(double *coord, int imagetmp, char **values) mu[nlocal][0] = atof(values[6]); mu[nlocal][1] = atof(values[7]); mu[nlocal][2] = atof(values[8]); + mu[nlocal][3] = sqrt(mu[nlocal][0]*mu[nlocal][0] + + mu[nlocal][1]*mu[nlocal][1] + + mu[nlocal][2]*mu[nlocal][2]); image[nlocal] = imagetmp; @@ -815,9 +761,6 @@ void AtomVecDipole::data_atom(double *coord, int imagetmp, char **values) v[nlocal][0] = 0.0; v[nlocal][1] = 0.0; v[nlocal][2] = 0.0; - omega[nlocal][0] = 0.0; - omega[nlocal][1] = 0.0; - omega[nlocal][2] = 0.0; atom->nlocal++; } @@ -833,36 +776,12 @@ int AtomVecDipole::data_atom_hybrid(int nlocal, char **values) mu[nlocal][0] = atof(values[1]); mu[nlocal][1] = atof(values[2]); mu[nlocal][2] = atof(values[3]); - + mu[nlocal][3] = sqrt(mu[nlocal][0]*mu[nlocal][0] + + mu[nlocal][1]*mu[nlocal][1] + + mu[nlocal][2]*mu[nlocal][2]); return 4; } -/* ---------------------------------------------------------------------- - unpack one line from Velocities section of data file -------------------------------------------------------------------------- */ - -void AtomVecDipole::data_vel(int m, char **values) -{ - v[m][0] = atof(values[0]); - v[m][1] = atof(values[1]); - v[m][2] = atof(values[2]); - omega[m][0] = atof(values[3]); - omega[m][1] = atof(values[4]); - omega[m][2] = atof(values[5]); -} - -/* ---------------------------------------------------------------------- - unpack hybrid quantities from one line in Velocities section of data file -------------------------------------------------------------------------- */ - -int AtomVecDipole::data_vel_hybrid(int m, char **values) -{ - omega[m][0] = atof(values[0]); - omega[m][1] = atof(values[1]); - omega[m][2] = atof(values[2]); - return 3; -} - /* ---------------------------------------------------------------------- return # of bytes of allocated memory ------------------------------------------------------------------------- */ @@ -880,9 +799,7 @@ bigint AtomVecDipole::memory_usage() if (atom->memcheck("f")) bytes += memory->usage(f,nmax,3); if (atom->memcheck("q")) bytes += memory->usage(q,nmax); - if (atom->memcheck("mu")) bytes += memory->usage(mu,nmax,3); - if (atom->memcheck("omega")) bytes += memory->usage(omega,nmax,3); - if (atom->memcheck("torque")) bytes += memory->usage(torque,nmax,3); + if (atom->memcheck("mu")) bytes += memory->usage(mu,nmax,4); return bytes; } diff --git a/src/DIPOLE/atom_vec_dipole.h b/src/DIPOLE/atom_vec_dipole.h index 010b1f5292..5cfd54c3d1 100644 --- a/src/DIPOLE/atom_vec_dipole.h +++ b/src/DIPOLE/atom_vec_dipole.h @@ -37,9 +37,7 @@ class AtomVecDipole : public AtomVec { void unpack_comm_vel(int, int, double *); int unpack_comm_one(int, double *); int pack_reverse(int, int, double *); - int pack_reverse_one(int, double *); void unpack_reverse(int, int *, double *); - int unpack_reverse_one(int, double *); int pack_border(int, int *, double *, int, int *); int pack_border_vel(int, int *, double *, int, int *); int pack_border_one(int, double *); @@ -54,8 +52,6 @@ class AtomVecDipole : public AtomVec { void create_atom(int, double *); void data_atom(double *, int, char **); int data_atom_hybrid(int, char **); - void data_vel(int, char **); - int data_vel_hybrid(int, char **); bigint memory_usage(); private: diff --git a/src/DIPOLE/pair_dipole_cut.cpp b/src/DIPOLE/pair_dipole_cut.cpp index a2e7c17bd4..e1d3b9345f 100644 --- a/src/DIPOLE/pair_dipole_cut.cpp +++ b/src/DIPOLE/pair_dipole_cut.cpp @@ -79,7 +79,6 @@ void PairDipoleCut::compute(int eflag, int vflag) double **mu = atom->mu; double **torque = atom->torque; int *type = atom->type; - double *dipole = atom->dipole; int nlocal = atom->nlocal; double *special_coul = force->special_coul; double *special_lj = force->special_lj; @@ -137,7 +136,7 @@ void PairDipoleCut::compute(int eflag, int vflag) forcecoulz += pre1*delz; } - if (dipole[itype] > 0.0 && dipole[jtype] > 0.0) { + if (mu[i][3] > 0.0 && mu[j][3] > 0.0) { r3inv = r2inv*rinv; r5inv = r3inv*r2inv; r7inv = r5inv*r2inv; @@ -167,7 +166,7 @@ void PairDipoleCut::compute(int eflag, int vflag) tjzcoul += -crossz + pre3 * (mu[j][0]*dely - mu[j][1]*delx); } - if (dipole[itype] > 0.0 && q[j] != 0.0) { + if (mu[i][3] > 0.0 && q[j] != 0.0) { r3inv = r2inv*rinv; r5inv = r3inv*r2inv; pidotr = mu[i][0]*delx + mu[i][1]*dely + mu[i][2]*delz; @@ -182,7 +181,7 @@ void PairDipoleCut::compute(int eflag, int vflag) tizcoul += pre2 * (mu[i][0]*dely - mu[i][1]*delx); } - if (dipole[jtype] > 0.0 && qtmp != 0.0) { + if (mu[j][3] > 0.0 && qtmp != 0.0) { r3inv = r2inv*rinv; r5inv = r3inv*r2inv; pjdotr = mu[j][0]*delx + mu[j][1]*dely + mu[j][2]*delz; @@ -234,11 +233,11 @@ void PairDipoleCut::compute(int eflag, int vflag) if (eflag) { if (rsq < cut_coulsq[itype][jtype]) { ecoul = qtmp*q[j]*rinv; - if (dipole[itype] > 0.0 && dipole[jtype] > 0.0) + if (mu[i][3] > 0.0 && mu[j][3] > 0.0) ecoul += r3inv*pdotp - 3.0*r5inv*pidotr*pjdotr; - if (dipole[itype] > 0.0 && q[j] != 0.0) + if (mu[i][3] > 0.0 && q[j] != 0.0) ecoul += -q[j]*r3inv*pidotr; - if (dipole[jtype] > 0.0 && qtmp != 0.0) + if (mu[j][3] > 0.0 && qtmp != 0.0) ecoul += qtmp*r3inv*pjdotr; ecoul *= factor_coul*qqrd2e; } else ecoul = 0.0; @@ -357,10 +356,8 @@ void PairDipoleCut::coeff(int narg, char **arg) void PairDipoleCut::init_style() { - if (!atom->q_flag || !atom->mu_flag || - !atom->torque_flag || atom->dipole == NULL) - error->all("Pair dipole/cut requires atom attributes " - "q, mu, torque, dipole"); + if (!atom->q_flag || !atom->mu_flag || !atom->torque_flag) + error->all("Pair dipole/cut requires atom attributes q, mu, torque"); int irequest = neighbor->request(this); } diff --git a/src/GPU/pair_gayberne_gpu.cpp b/src/GPU/pair_gayberne_gpu.cpp index c8fc0cf35b..a815530fac 100644 --- a/src/GPU/pair_gayberne_gpu.cpp +++ b/src/GPU/pair_gayberne_gpu.cpp @@ -134,19 +134,21 @@ void PairGayBerneGPU::init_style() { if (force->pair_match("gpu",0) == NULL) error->all("Cannot use pair hybrid with multiple GPU pair styles"); - if (!atom->quat_flag || !atom->torque_flag || !atom->avec->shape_type) - error->all("Pair gayberne requires atom attributes quat, torque, shape"); - if (atom->radius_flag) - error->all("Pair gayberne cannot be used with atom attribute diameter"); + if (!atom->ellipsoid_flag) + error->all("Pair gayberne requires atom style ellipsoid"); // per-type shape precalculations + // require that atom shapes are identical within each type + for (int i = 1; i <= atom->ntypes; i++) { + if (!atom->shape_consistency(i,shape1[i][0],shape1[i][1],shape1[i][2])) + error->all("Pair gayberne requires atoms with same type have same shape"); if (setwell[i]) { - double *one = atom->shape[i]; - shape[i][0] = one[0]*one[0]; - shape[i][1] = one[1]*one[1]; - shape[i][2] = one[2]*one[2]; - lshape[i] = (one[0]*one[1]+one[2]*one[2])*sqrt(one[0]*one[1]); + shape2[i][0] = shape1[i][0]*shape1[i][0]; + shape2[i][1] = shape1[i][1]*shape1[i][1]; + shape2[i][2] = shape1[i][2]*shape1[i][2]; + lshape[i] = (shape1[i][0]*shape1[i][1]+shape1[i][2]*shape1[i][2]) * + sqrt(shape1[i][0]*shape1[i][1]); } } @@ -169,7 +171,7 @@ void PairGayBerneGPU::init_style() double cell_size = sqrt(maxcut) + neighbor->skin; bool init_ok = gb_gpu_init(atom->ntypes+1, gamma, upsilon, mu, - shape, well, cutsq, sigma, epsilon, lshape, form, + shape1, well, cutsq, sigma, epsilon, lshape, form, lj1, lj2, lj3, lj4, offset, force->special_lj, atom->nlocal, atom->nlocal+atom->nghost, 300, cell_size, gpu_mode, screen); @@ -227,7 +229,7 @@ void PairGayBerneGPU::cpu_compute(int start, int eflag, int vflag) MathExtra::quat_to_mat_trans(quat[i],a1); MathExtra::diag_times3(well[itype],a1,temp); MathExtra::transpose_times3(a1,temp,b1); - MathExtra::diag_times3(shape[itype],a1,temp); + MathExtra::diag_times3(shape2[itype],a1,temp); MathExtra::transpose_times3(a1,temp,g1); } @@ -271,7 +273,7 @@ void PairGayBerneGPU::cpu_compute(int start, int eflag, int vflag) MathExtra::quat_to_mat_trans(quat[j],a2); MathExtra::diag_times3(well[jtype],a2,temp); MathExtra::transpose_times3(a2,temp,b2); - MathExtra::diag_times3(shape[jtype],a2,temp); + MathExtra::diag_times3(shape2[jtype],a2,temp); MathExtra::transpose_times3(a2,temp,g2); one_eng = gayberne_lj(j,i,a2,b2,g2,r12,rsq,fforce,rtor); ttor[0] = ttor[1] = ttor[2] = 0.0; @@ -286,7 +288,7 @@ void PairGayBerneGPU::cpu_compute(int start, int eflag, int vflag) MathExtra::quat_to_mat_trans(quat[j],a2); MathExtra::diag_times3(well[jtype],a2,temp); MathExtra::transpose_times3(a2,temp,b2); - MathExtra::diag_times3(shape[jtype],a2,temp); + MathExtra::diag_times3(shape2[jtype],a2,temp); MathExtra::transpose_times3(a2,temp,g2); one_eng = gayberne_analytic(i,j,a1,a2,b1,b2,g1,g2,r12,rsq, fforce,ttor,rtor); @@ -343,7 +345,7 @@ void PairGayBerneGPU::cpu_compute(int *nbors, int start, int eflag, int vflag) MathExtra::quat_to_mat_trans(quat[i],a1); MathExtra::diag_times3(well[itype],a1,temp); MathExtra::transpose_times3(a1,temp,b1); - MathExtra::diag_times3(shape[itype],a1,temp); + MathExtra::diag_times3(shape2[itype],a1,temp); MathExtra::transpose_times3(a1,temp,g1); } @@ -389,7 +391,7 @@ void PairGayBerneGPU::cpu_compute(int *nbors, int start, int eflag, int vflag) MathExtra::quat_to_mat_trans(quat[j],a2); MathExtra::diag_times3(well[jtype],a2,temp); MathExtra::transpose_times3(a2,temp,b2); - MathExtra::diag_times3(shape[jtype],a2,temp); + MathExtra::diag_times3(shape2[jtype],a2,temp); MathExtra::transpose_times3(a2,temp,g2); one_eng = gayberne_lj(j,i,a2,b2,g2,r12,rsq,fforce,rtor); ttor[0] = ttor[1] = ttor[2] = 0.0; @@ -404,7 +406,7 @@ void PairGayBerneGPU::cpu_compute(int *nbors, int start, int eflag, int vflag) MathExtra::quat_to_mat_trans(quat[j],a2); MathExtra::diag_times3(well[jtype],a2,temp); MathExtra::transpose_times3(a2,temp,b2); - MathExtra::diag_times3(shape[jtype],a2,temp); + MathExtra::diag_times3(shape2[jtype],a2,temp); MathExtra::transpose_times3(a2,temp,g2); one_eng = gayberne_analytic(i,j,a1,a2,b1,b2,g1,g2,r12,rsq, fforce,ttor,rtor); diff --git a/src/GRANULAR/Install.sh b/src/GRANULAR/Install.sh index 0b8c3bf00a..7a60972446 100644 --- a/src/GRANULAR/Install.sh +++ b/src/GRANULAR/Install.sh @@ -2,7 +2,6 @@ if (test $1 = 1) then - cp atom_vec_granular.cpp .. cp fix_freeze.cpp .. cp fix_pour.cpp .. cp fix_wall_gran.cpp .. @@ -10,7 +9,6 @@ if (test $1 = 1) then cp pair_gran_hooke.cpp .. cp pair_gran_hooke_history.cpp .. - cp atom_vec_granular.h .. cp fix_freeze.h .. cp fix_pour.h .. cp fix_wall_gran.h .. @@ -20,7 +18,6 @@ if (test $1 = 1) then elif (test $1 = 0) then - rm ../atom_vec_granular.cpp rm ../fix_freeze.cpp rm ../fix_pour.cpp rm ../fix_wall_gran.cpp @@ -28,7 +25,6 @@ elif (test $1 = 0) then rm ../pair_gran_hooke.cpp rm ../pair_gran_hooke_history.cpp - rm ../atom_vec_granular.h rm ../fix_freeze.h rm ../fix_pour.h rm ../fix_wall_gran.h diff --git a/src/GRANULAR/fix_pour.cpp b/src/GRANULAR/fix_pour.cpp index 73b9d49914..1248c05db0 100644 --- a/src/GRANULAR/fix_pour.cpp +++ b/src/GRANULAR/fix_pour.cpp @@ -472,7 +472,6 @@ void FixPour::pre_exchange() m = atom->nlocal - 1; atom->type[m] = ntype; atom->radius[m] = radtmp; - atom->density[m] = denstmp; atom->rmass[m] = 4.0*PI/3.0 * radtmp*radtmp*radtmp * denstmp; atom->mask[m] = 1 | groupbit; atom->v[m][0] = vxtmp; diff --git a/src/GRANULAR/fix_wall_gran.cpp b/src/GRANULAR/fix_wall_gran.cpp index 57f904f99c..d8da32fc1a 100644 --- a/src/GRANULAR/fix_wall_gran.cpp +++ b/src/GRANULAR/fix_wall_gran.cpp @@ -46,8 +46,8 @@ FixWallGran::FixWallGran(LAMMPS *lmp, int narg, char **arg) : { if (narg < 10) error->all("Illegal fix wall/gran command"); - if (!atom->radius_flag || !atom->omega_flag || !atom->torque_flag) - error->all("Fix wall/gran requires atom attributes radius, omega, torque"); + if (!atom->sphere_flag) + error->all("Fix wall/gran requires atom style sphere"); restart_peratom = 1; create_attribute = 1; diff --git a/src/GRANULAR/pair_gran_hooke_history.cpp b/src/GRANULAR/pair_gran_hooke_history.cpp index cd4f8e4eb6..0ce30565dd 100644 --- a/src/GRANULAR/pair_gran_hooke_history.cpp +++ b/src/GRANULAR/pair_gran_hooke_history.cpp @@ -364,8 +364,8 @@ void PairGranHookeHistory::init_style() // error and warning checks - if (!atom->radius_flag || !atom->omega_flag || !atom->torque_flag) - error->all("Pair granular requires atom attributes radius, omega, torque"); + if (!atom->sphere_flag) + error->all("Pair granular requires atom style sphere"); if (comm->ghost_velocity == 0) error->all("Pair granular requires ghost atoms store velocity"); diff --git a/src/PERI/atom_vec_peri.cpp b/src/PERI/atom_vec_peri.cpp index 68110b329c..d23fe810f3 100644 --- a/src/PERI/atom_vec_peri.cpp +++ b/src/PERI/atom_vec_peri.cpp @@ -47,7 +47,8 @@ AtomVecPeri::AtomVecPeri(LAMMPS *lmp, int narg, char **arg) : size_data_vel = 4; xcol_data = 5; - atom->vfrac_flag = atom->density_flag = atom->rmass_flag = 1; + atom->peri_flag = 1; + atom->vfrac_flag = atom->rmass_flag = 1; } /* ---------------------------------------------------------------------- @@ -73,7 +74,6 @@ void AtomVecPeri::grow(int n) f = memory->grow(atom->f,nmax,3,"atom:f"); vfrac = memory->grow(atom->vfrac,nmax,"atom:vfrac"); - density = memory->grow(atom->density,nmax,"atom:density"); rmass = memory->grow(atom->rmass,nmax,"atom:rmass"); s0 = memory->grow(atom->s0,nmax,"atom:s0"); x0 = memory->grow(atom->x0,nmax,3,"atom:x0"); @@ -92,7 +92,7 @@ void AtomVecPeri::grow_reset() tag = atom->tag; type = atom->type; mask = atom->mask; image = atom->image; x = atom->x; v = atom->v; f = atom->f; - vfrac = atom->vfrac; density = atom->density; rmass = atom->rmass; + vfrac = atom->vfrac; rmass = atom->rmass; s0 = atom->s0; x0 = atom->x0; } @@ -112,7 +112,6 @@ void AtomVecPeri::copy(int i, int j) v[j][2] = v[i][2]; vfrac[j] = vfrac[i]; - density[j] = density[i]; rmass[j] = rmass[i]; s0[j] = s0[i]; x0[j][0] = x0[i][0]; @@ -546,7 +545,6 @@ int AtomVecPeri::pack_exchange(int i, double *buf) buf[m++] = image[i]; buf[m++] = vfrac[i]; - buf[m++] = density[i]; buf[m++] = rmass[i]; buf[m++] = s0[i]; buf[m++] = x0[i][0]; @@ -581,7 +579,6 @@ int AtomVecPeri::unpack_exchange(double *buf) image[nlocal] = static_cast (buf[m++]); vfrac[nlocal] = buf[m++]; - density[nlocal] = buf[m++]; rmass[nlocal] = buf[m++]; s0[nlocal] = buf[m++]; x0[nlocal][0] = buf[m++]; @@ -608,7 +605,7 @@ int AtomVecPeri::size_restart() int i; int nlocal = atom->nlocal; - int n = 18 * nlocal; + int n = 17 * nlocal; if (atom->nextra_restart) for (int iextra = 0; iextra < atom->nextra_restart; iextra++) @@ -639,7 +636,6 @@ int AtomVecPeri::pack_restart(int i, double *buf) buf[m++] = v[i][2]; buf[m++] = vfrac[i]; - buf[m++] = density[i]; buf[m++] = rmass[i]; buf[m++] = s0[i]; buf[m++] = x0[i][0]; @@ -680,7 +676,6 @@ int AtomVecPeri::unpack_restart(double *buf) v[nlocal][2] = buf[m++]; vfrac[nlocal] = buf[m++]; - density[nlocal] = buf[m++]; rmass[nlocal] = buf[m++]; s0[nlocal] = buf[m++]; x0[nlocal][0] = buf[m++]; @@ -719,8 +714,7 @@ void AtomVecPeri::create_atom(int itype, double *coord) v[nlocal][2] = 0.0; vfrac[nlocal] = 1.0; - density[nlocal] = 1.0; - rmass[nlocal] = density[nlocal]; + rmass[nlocal] = 1.0; s0[nlocal] = DBL_MAX; x0[nlocal][0] = coord[0]; x0[nlocal][1] = coord[1]; @@ -748,8 +742,7 @@ void AtomVecPeri::data_atom(double *coord, int imagetmp, char **values) error->one("Invalid atom type in Atoms section of data file"); vfrac[nlocal] = atof(values[2]); - density[nlocal] = atof(values[3]); - rmass[nlocal] = density[nlocal]; + rmass[nlocal] = atof(values[3]); if (rmass[nlocal] <= 0.0) error->one("Invalid mass value"); x[nlocal][0] = coord[0]; @@ -780,8 +773,7 @@ void AtomVecPeri::data_atom(double *coord, int imagetmp, char **values) int AtomVecPeri::data_atom_hybrid(int nlocal, char **values) { vfrac[nlocal] = atof(values[0]); - density[nlocal] = atof(values[1]); - rmass[nlocal] = density[nlocal]; + rmass[nlocal] = atof(values[1]); if (rmass[nlocal] <= 0.0) error->one("Invalid mass value"); s0[nlocal] = DBL_MAX; @@ -809,7 +801,6 @@ bigint AtomVecPeri::memory_usage() if (atom->memcheck("f")) bytes += memory->usage(f,nmax,3); if (atom->memcheck("vfrac")) bytes += memory->usage(vfrac,nmax); - if (atom->memcheck("density")) bytes += memory->usage(density,nmax); if (atom->memcheck("rmass")) bytes += memory->usage(rmass,nmax); if (atom->memcheck("s0")) bytes += memory->usage(s0,nmax); if (atom->memcheck("x0")) bytes += memory->usage(x0,nmax,3); diff --git a/src/SRD/fix_srd.cpp b/src/SRD/fix_srd.cpp index ca2aa28b64..41f9040625 100644 --- a/src/SRD/fix_srd.cpp +++ b/src/SRD/fix_srd.cpp @@ -40,8 +40,6 @@ using namespace LAMMPS_NS; enum{SLIP,NOSLIP}; enum{SPHERE,ELLIPSOID,WALL}; -enum{SPHERE_SHAPE,SPHERE_RADIUS}; -enum{ANGULAR_OMEGA,ANGULAR_ANGMOM}; enum{INSIDE_ERROR,INSIDE_WARN,INSIDE_IGNORE}; enum{BIG_MOVE,SRD_MOVE,SRD_ROTATE}; enum{CUBIC_ERROR,CUBIC_WARN}; @@ -291,10 +289,8 @@ void FixSRD::init() if (bigexist && comm->ghost_velocity == 0) error->all("Fix srd requires ghost atoms store velocity"); - if (bigexist && !atom->radius_flag && !atom->avec->shape_type) - error->all("Fix SRD requires atom attribute radius or shape"); - if (bigexist && !atom->angmom_flag && !atom->omega_flag) - error->all("Fix SRD requires atom attribute angmom or omega"); + if (bigexist && !atom->sphere_flag && !atom->ellipsoid_flag) + error->all("Fix SRD requires atom style sphere or ellipsoid"); if (bigexist && atom->angmom_flag && atom->omega_flag) error->all("Fix SRD cannot have both atom attributes angmom and omega"); if (bigexist && collidestyle == NOSLIP && !atom->torque_flag) @@ -897,15 +893,6 @@ void FixSRD::reset_velocities() vold[1] /= n; vold[2] /= n; - // if (triclinic && streamflag) { - // xlamda = vbin[i].xctr; - // vstream[0] = h_rate[0]*xlamda[0] + h_rate[5]*xlamda[1] + - // h_rate[4]*xlamda[2] + h_ratelo[0]; - // vstream[1] = h_rate[1]*xlamda[1] + h_rate[3]*xlamda[2] + h_ratelo[1]; - // vstream[2] = h_rate[2]*xlamda[2] + h_ratelo[2]; - // vnew = vstream; - //} else vnew = vold; - vnew = vold; irandom = static_cast (6.0*vbin[i].random); @@ -2109,7 +2096,6 @@ void FixSRD::parameterize() double *radius = atom->radius; double **shape = atom->shape; - int *type = atom->type; int *mask = atom->mask; int nlocal = atom->nlocal; @@ -2119,24 +2105,20 @@ void FixSRD::parameterize() for (int i = 0; i < nlocal; i++) if (mask[i] & biggroupbit) { - if (radius) { - if (radius[i] == 0.0) - error->one("Big particle in fix srd cannot be point particle"); + if (radius && radius[i] > 0.0) { maxbigdiam = MAX(maxbigdiam,2.0*radius[i]); minbigdiam = MIN(minbigdiam,2.0*radius[i]); - } else { - if (shape[type[i]][0] == 0.0) - error->one("Big particle in fix srd cannot be point particle"); - maxbigdiam = MAX(maxbigdiam,2.0*shape[type[i]][0]); - maxbigdiam = MAX(maxbigdiam,2.0*shape[type[i]][1]); - maxbigdiam = MAX(maxbigdiam,2.0*shape[type[i]][2]); - minbigdiam = MIN(minbigdiam,2.0*shape[type[i]][0]); - minbigdiam = MIN(minbigdiam,2.0*shape[type[i]][1]); - minbigdiam = MIN(minbigdiam,2.0*shape[type[i]][2]); - if (shape[type[i]][0] != shape[type[i]][1] || - shape[type[i]][0] != shape[type[i]][2]) + } 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 + error->one("Big particle in fix srd cannot be point particle"); } double tmp = maxbigdiam; @@ -2159,6 +2141,7 @@ void FixSRD::parameterize() double *rmass = atom->rmass; double *mass = atom->mass; + int *type = atom->type; int flag = 0; mass_srd = 0.0; @@ -2166,8 +2149,10 @@ void FixSRD::parameterize() for (int i = 0; i < nlocal; i++) if (mask[i] & groupbit) { if (rmass) { - if (mass_srd == 0.0) mass_srd = rmass[i]; - else if (rmass[i] != mass_srd) flag = 1; + if (mass_srd == 0.0) { + mass_srd = rmass[i]; + printf("MASS SRD %g\n",rmass[i]); + } else if (rmass[i] != mass_srd) flag = 1; } else { if (mass_srd == 0.0) mass_srd = mass[type[i]]; else if (mass[type[i]] != mass_srd) flag = 1; @@ -2188,6 +2173,8 @@ 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 @@ -2204,19 +2191,18 @@ void FixSRD::parameterize() if (dimension == 3) { for (int i = 0; i < nlocal; i++) if (mask[i] & biggroupbit) { - if (radius) + if (radius && radius[i] > 0.0) volbig += 4.0/3.0*PI*radius[i]*radius[i]*radius[i]; - else - volbig += 4.0/3.0*PI * - shape[type[i]][0]*shape[type[i]][1]*shape[type[i]][2]; + else if (shape && shape[i][0] > 0.0) + volbig += 4.0/3.0*PI * shape[i][0]*shape[i][1]*shape[i][2]; } } else { for (int i = 0; i < nlocal; i++) if (mask[i] & biggroupbit) { - if (radius) + if (radius && radius[i] > 0.0) volbig += PI*radius[i]*radius[i]; - else - volbig += PI*shape[type[i]][0]*shape[type[i]][1]; + else if (shape && shape[i][0] > 0.0) + volbig += PI*shape[i][0]*shape[i][1]; } } @@ -2236,10 +2222,7 @@ void FixSRD::parameterize() mass_big = 0.0; for (int i = 0; i < nlocal; i++) - if (mask[i] & biggroupbit) { - if (rmass) mass_big += rmass[i]; - else mass_big += mass[type[i]]; - } + if (mask[i] & biggroupbit) mass_big += rmass[i]; tmp = mass_big; MPI_Allreduce(&tmp,&mass_big,1,MPI_DOUBLE,MPI_SUM,world); @@ -2425,40 +2408,21 @@ void FixSRD::big_static() double **shape = atom->shape; int *type = atom->type; - int omega_flag = atom->omega_flag; - double skinhalf = 0.5 * neighbor->skin; for (int k = 0; k < nbig; k++) { i = biglist[k].index; - if (radius) { + if (radius && radius[i] > 0.0) { biglist[k].type = SPHERE; - biglist[k].typesphere = SPHERE_RADIUS; - biglist[k].typeangular = ANGULAR_OMEGA; rad = radfactor*radius[i]; biglist[k].radius = rad; biglist[k].radsq = rad*rad; biglist[k].cutbinsq = (rad+skinhalf) * (rad+skinhalf); - } else if (shape[type[i]][0] == shape[type[i]][1] && - shape[type[i]][0] == shape[type[i]][2]) { - biglist[k].type = SPHERE; - biglist[k].typesphere = SPHERE_SHAPE; - - // either atom->omega is defined or atom->angmom, cannot both be defined - - if (omega_flag) biglist[k].typeangular = ANGULAR_OMEGA; - else biglist[k].typeangular = ANGULAR_ANGMOM; - - rad = radfactor*shape[type[i]][0]; - biglist[k].radius = rad; - biglist[k].radsq = rad*rad; - biglist[k].cutbinsq = (rad+skinhalf) * (rad+skinhalf); - } else { + } else if (shape && shape[i][0] > 0.0) { biglist[k].type = ELLIPSOID; - biglist[k].typeangular = ANGULAR_ANGMOM; - arad = radfactor*shape[type[i]][0]; - brad = radfactor*shape[type[i]][1]; - crad = radfactor*shape[type[i]][2]; + arad = radfactor*shape[i][0]; + brad = radfactor*shape[i][1]; + crad = radfactor*shape[i][2]; biglist[k].aradsqinv = 1.0/(arad*arad); biglist[k].bradsqinv = 1.0/(brad*brad); biglist[k].cradsqinv = 1.0/(crad*crad); @@ -2478,44 +2442,32 @@ void FixSRD::big_static() void FixSRD::big_dynamic() { - int i,itype; - double inertia; + int i; double **omega = atom->omega; double **angmom = atom->angmom; double **quat = atom->quat; double **shape = atom->shape; - double *mass = atom->mass; - int *type = atom->type; + double *rmass = atom->rmass; for (int k = 0; k < nbig; k++) { i = biglist[k].index; - // ellipsoid with angmom - // calculate ex,ey,ez and omega from quaternion and angmom - - if (biglist[k].type == ELLIPSOID) { - exyz_from_q(quat[i],biglist[k].ex,biglist[k].ey,biglist[k].ez); - omega_from_mq(angmom[i],biglist[k].ex,biglist[k].ey,biglist[k].ez, - mass[type[i]],shape[type[i]],biglist[k].omega); - - // sphere with omega and shape or radius + // sphere with omega // set omega from atom->omega directly - } else if (biglist[k].typeangular == ANGULAR_OMEGA) { + if (biglist[k].type == SPHERE) { biglist[k].omega[0] = omega[i][0]; biglist[k].omega[1] = omega[i][1]; biglist[k].omega[2] = omega[i][2]; - // sphere with angmom and shape - // calculate omega from angmom + // ellipsoid with angmom + // calculate ex,ey,ez and omega from quaternion and angmom - } else if (biglist[k].typeangular == ANGULAR_ANGMOM) { - itype = type[i]; - inertia = INERTIA * mass[itype]*shape[itype][0]*shape[itype][0]; - biglist[k].omega[0] = angmom[i][0] / inertia; - biglist[k].omega[1] = angmom[i][1] / inertia; - biglist[k].omega[2] = angmom[i][2] / inertia; + } else if (biglist[k].type == ELLIPSOID) { + exyz_from_q(quat[i],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); } } } diff --git a/src/SRD/fix_srd.h b/src/SRD/fix_srd.h index 3c51f9bda7..b7d1fd188b 100644 --- a/src/SRD/fix_srd.h +++ b/src/SRD/fix_srd.h @@ -91,8 +91,6 @@ class FixSRD : public Fix { struct Big { int index; // local index of particle/wall int type; // SPHERE or ELLIPSOID or WALL - int typesphere; // SPHERE_SHAPE or SPHERE_RADIUS - int typeangular; // ANGULAR_OMEGA or ANGULAR_ANGMOM double radius,radsq; // radius of sphere double aradsqinv; // 3 ellipsoid radii double bradsqinv; diff --git a/src/USER-EFF/atom_vec_electron.cpp b/src/USER-EFF/atom_vec_electron.cpp index 4c4e1d1c3a..7ddb37d004 100644 --- a/src/USER-EFF/atom_vec_electron.cpp +++ b/src/USER-EFF/atom_vec_electron.cpp @@ -48,6 +48,7 @@ AtomVecElectron::AtomVecElectron(LAMMPS *lmp, int narg, char **arg) : size_data_vel = 5; xcol_data = 6; + atom->electron_flag = 1; atom->q_flag = atom->spin_flag = atom->eradius_flag = atom->ervel_flag = atom->erforce_flag = 1; } diff --git a/src/USER-EFF/compute_ke_atom_eff.cpp b/src/USER-EFF/compute_ke_atom_eff.cpp index d0c0fed1a5..33c900f4b3 100644 --- a/src/USER-EFF/compute_ke_atom_eff.cpp +++ b/src/USER-EFF/compute_ke_atom_eff.cpp @@ -43,8 +43,8 @@ ComputeKEAtomEff::ComputeKEAtomEff(LAMMPS *lmp, int narg, char **arg) : // error check - if (!atom->ervel_flag) - error->all("Compute ke/atom/eff requires atom attribute ervel"); + if (!atom->electron_flag) + error->all("Compute ke/atom/eff requires atom style electron"); } /* ---------------------------------------------------------------------- */ diff --git a/src/USER-EFF/compute_ke_eff.cpp b/src/USER-EFF/compute_ke_eff.cpp index 15ee955fdb..68aaa2c25b 100644 --- a/src/USER-EFF/compute_ke_eff.cpp +++ b/src/USER-EFF/compute_ke_eff.cpp @@ -39,8 +39,8 @@ ComputeKEEff::ComputeKEEff(LAMMPS *lmp, int narg, char **arg) : // error check - if (!atom->ervel_flag) - error->all("Compute ke/eff requires atom attribute ervel"); + if (!atom->electron_flag) + error->all("Compute ke/eff requires atom style electron"); } /* ---------------------------------------------------------------------- */ diff --git a/src/USER-EFF/compute_temp_deform_eff.cpp b/src/USER-EFF/compute_temp_deform_eff.cpp index 41430e5b4e..7bb53b9085 100644 --- a/src/USER-EFF/compute_temp_deform_eff.cpp +++ b/src/USER-EFF/compute_temp_deform_eff.cpp @@ -43,8 +43,8 @@ ComputeTempDeformEff::ComputeTempDeformEff(LAMMPS *lmp, int narg, char **arg) : { if (narg != 3) error->all("Illegal compute temp/deform/eff command"); - if (!atom->spin_flag || !atom->ervel_flag) - error->all("Compute temp/deform/eff requires atom attributes spin, ervel"); + if (!atom->electron_flag) + error->all("Compute temp/deform/eff requires atom style electron"); scalar_flag = vector_flag = 1; size_vector = 6; diff --git a/src/USER-EFF/compute_temp_eff.cpp b/src/USER-EFF/compute_temp_eff.cpp index 29f9315995..4f7cc40c1b 100644 --- a/src/USER-EFF/compute_temp_eff.cpp +++ b/src/USER-EFF/compute_temp_eff.cpp @@ -35,10 +35,8 @@ using namespace LAMMPS_NS; ComputeTempEff::ComputeTempEff(LAMMPS *lmp, int narg, char **arg) : Compute(lmp, narg, arg) { - // error check - - if (!atom->spin_flag || !atom->ervel_flag) - error->all("Compute temp/eff requires atom attributes spin, ervel"); + if (!atom->electron_flag) + error->all("Compute temp/eff requires atom style electron"); scalar_flag = vector_flag = 1; size_vector = 6; diff --git a/src/USER-EFF/compute_temp_region_eff.cpp b/src/USER-EFF/compute_temp_region_eff.cpp index 8b1ac402c7..9bcabee2a8 100644 --- a/src/USER-EFF/compute_temp_region_eff.cpp +++ b/src/USER-EFF/compute_temp_region_eff.cpp @@ -35,8 +35,8 @@ using namespace LAMMPS_NS; ComputeTempRegionEff::ComputeTempRegionEff(LAMMPS *lmp, int narg, char **arg) : Compute(lmp, narg, arg) { - if (!atom->spin_flag || !atom->ervel_flag) - error->all("Compute temp/region/eff requires atom attributes spin, ervel"); + if (!atom->electron_flag) + error->all("Compute temp/region/eff requires atom style electron"); if (narg != 4) error->all("Illegal compute temp/region/eff command"); diff --git a/src/USER-EFF/fix_nh_eff.cpp b/src/USER-EFF/fix_nh_eff.cpp index 61a6ea506e..70f2544c2d 100644 --- a/src/USER-EFF/fix_nh_eff.cpp +++ b/src/USER-EFF/fix_nh_eff.cpp @@ -30,10 +30,8 @@ enum{NOBIAS,BIAS}; FixNHEff::FixNHEff(LAMMPS *lmp, int narg, char **arg) : FixNH(lmp, narg, arg) { - if (!atom->spin_flag || !atom->eradius_flag || - !atom->ervel_flag || !atom->erforce_flag) - error->all("Fix nvt/nph/npt/eff requires atom attributes " - "spin, eradius, ervel, erforce"); + if (!atom->electron_flag) + error->all("Fix nvt/nph/npt/eff requires atom style electron"); } /* ---------------------------------------------------------------------- diff --git a/src/USER-EFF/fix_nve_eff.cpp b/src/USER-EFF/fix_nve_eff.cpp index 048b251f9c..134317ddeb 100644 --- a/src/USER-EFF/fix_nve_eff.cpp +++ b/src/USER-EFF/fix_nve_eff.cpp @@ -33,10 +33,8 @@ using namespace LAMMPS_NS; FixNVEEff::FixNVEEff(LAMMPS *lmp, int narg, char **arg) : Fix(lmp, narg, arg) { - if (!atom->spin_flag || !atom->eradius_flag || - !atom->ervel_flag || !atom->erforce_flag) - error->all("Fix nve/eff requires atom attributes " - "spin, eradius, ervel, erforce"); + if (!atom->electron_flag) + error->all("Fix nve/eff requires atom style electron"); time_integrate = 1; } diff --git a/src/atom.cpp b/src/atom.cpp index bfb78d6721..d31f20c1e3 100644 --- a/src/atom.cpp +++ b/src/atom.cpp @@ -70,8 +70,8 @@ Atom::Atom(LAMMPS *lmp) : Pointers(lmp) molecule = NULL; q = NULL; mu = NULL; - quat = omega = angmom = torque = NULL; - radius = density = rmass = NULL; + quat = omega = angmom = torque = shape = NULL; + radius = rmass = NULL; vfrac = s0 = NULL; x0 = NULL; @@ -96,23 +96,20 @@ Atom::Atom(LAMMPS *lmp) : Pointers(lmp) improper_type = improper_atom1 = improper_atom2 = NULL; improper_atom3 = improper_atom4 = NULL; - // initialize atom array existence flags + // initialize atom style and array existence flags // customize by adding new flag - molecule_flag = 0; - q_flag = mu_flag = 0; - quat_flag = omega_flag = angmom_flag = torque_flag = 0; - radius_flag = density_flag = rmass_flag = vfrac_flag = 0; - spin_flag = eradius_flag = ervel_flag = erforce_flag = 0; + sphere_flag = ellipsoid_flag = peri_flag = dipole_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; + vfrac_flag = spin_flag = eradius_flag = ervel_flag = erforce_flag = 0; // ntype-length arrays mass = NULL; mass_setflag = NULL; - shape = NULL; - shape_setflag = NULL; - dipole = NULL; - dipole_setflag = NULL; // callback lists & extra restart info @@ -171,12 +168,12 @@ 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(density); memory->destroy(rmass); memory->destroy(vfrac); memory->destroy(s0); @@ -220,10 +217,6 @@ Atom::~Atom() delete [] mass; delete [] mass_setflag; - memory->destroy(shape); - delete [] shape_setflag; - delete [] dipole; - delete [] dipole_setflag; // delete extra arrays @@ -256,14 +249,16 @@ void Atom::create_avec(const char *style, int narg, char **arg) delete [] atom_style; if (avec) delete avec; - // unset atom array existence flags that may have been set by old avec + // unset atom style and array existence flags + // may have been set by old avec // customize by adding new flag - molecule_flag = 0; - q_flag = mu_flag = 0; - quat_flag = omega_flag = angmom_flag = torque_flag = 0; - radius_flag = density_flag = rmass_flag = vfrac_flag = 0; - spin_flag = eradius_flag = ervel_flag = erforce_flag = 0; + sphere_flag = ellipsoid_flag = peri_flag = dipole_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; + vfrac_flag = spin_flag = eradius_flag = ervel_flag = erforce_flag = 0; avec = new_avec(style,narg,arg); int n = strlen(style) + 1; @@ -309,8 +304,6 @@ void Atom::init() // check arrays that are atom type in length check_mass(); - check_shape(); - check_dipole(); // setup of firstgroup @@ -1066,16 +1059,6 @@ void Atom::allocate_type_arrays() mass_setflag = new int[ntypes+1]; for (int itype = 1; itype <= ntypes; itype++) mass_setflag[itype] = 0; } - if (avec->shape_type) { - memory->create(shape,ntypes+1,3,"atom:shape"); - shape_setflag = new int[ntypes+1]; - for (int itype = 1; itype <= ntypes; itype++) shape_setflag[itype] = 0; - } - if (avec->dipole_type) { - dipole = new double[ntypes+1]; - dipole_setflag = new int[ntypes+1]; - for (int itype = 1; itype <= ntypes; itype++) dipole_setflag[itype] = 0; - } } /* ---------------------------------------------------------------------- @@ -1161,161 +1144,61 @@ void Atom::check_mass() } /* ---------------------------------------------------------------------- - set particle shape and flag it as set - called from reading of data file + check that radii of all particles of itype are the same + return 1 if true, else return 0 + also return the radius value for that type ------------------------------------------------------------------------- */ -void Atom::set_shape(const char *str) +int Atom::radius_consistency(int itype, double &rad) { - if (shape == NULL) error->all("Cannot set shape for this atom style"); - - int itype; - double a,b,c; - int n = sscanf(str,"%d %lg %lg %lg",&itype,&a,&b,&c); - if (n != 4) error->all("Invalid shape line in data file"); - - if (itype < 1 || itype > ntypes) error->all("Invalid type for shape set"); - - // store shape as radius, though specified as diameter - - shape[itype][0] = 0.5*a; - shape[itype][1] = 0.5*b; - shape[itype][2] = 0.5*c; - shape_setflag[itype] = 1; - - if (shape[itype][0] < 0.0 || shape[itype][1] < 0.0 || - shape[itype][2] < 0.0) - error->all("Invalid shape value"); - if (shape[itype][0] > 0.0 || shape[itype][1] > 0.0 || - shape[itype][2] > 0.0) { - if (shape[itype][0] == 0.0 || shape[itype][1] == 0.0 || - shape[itype][2] == 0.0) - error->all("Invalid shape value"); + double value = -1.0; + int flag = 0; + for (int i = 0; i < nlocal; i++) { + if (type[i] != itype) continue; + if (value < 0.0) value = radius[i]; + else if (value != radius[i]) flag = 1; } + + int flagall; + MPI_Allreduce(&flag,&flagall,1,MPI_INT,MPI_SUM,world); + if (flagall) return 0; + + MPI_Allreduce(&value,&rad,1,MPI_DOUBLE,MPI_MAX,world); + return 1; } /* ---------------------------------------------------------------------- - set one or more particle shapes and flag them as set - called from reading of input script + 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 ------------------------------------------------------------------------- */ -void Atom::set_shape(int narg, char **arg) +int Atom::shape_consistency(int itype, + double &shapex, double &shapey, double &shapez) { - if (shape == NULL) error->all("Cannot set shape for this atom style"); - - int lo,hi; - force->bounds(arg[0],ntypes,lo,hi); - if (lo < 1 || hi > ntypes) - error->all("Invalid type for shape set"); - - // store shape as radius, though specified as diameter - - for (int itype = lo; itype <= hi; itype++) { - shape[itype][0] = 0.5*atof(arg[1]); - shape[itype][1] = 0.5*atof(arg[2]); - shape[itype][2] = 0.5*atof(arg[3]); - shape_setflag[itype] = 1; - - if (shape[itype][0] < 0.0 || shape[itype][1] < 0.0 || - shape[itype][2] < 0.0) - error->all("Invalid shape value"); - if (shape[itype][0] > 0.0 || shape[itype][1] > 0.0 || - shape[itype][2] > 0.0) { - if (shape[itype][0] == 0.0 || shape[itype][1] == 0.0 || - shape[itype][2] == 0.0) - error->all("Invalid shape value"); - } + double one[3]; + one[0] = one[1] = one[2] = -1.0; + int flag = 0; + for (int i = 0; i < nlocal; i++) { + if (type[i] != itype) continue; + 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; } -} -/* ---------------------------------------------------------------------- - set all particle shapes as read in from restart file -------------------------------------------------------------------------- */ + int flagall; + MPI_Allreduce(&flag,&flagall,1,MPI_INT,MPI_SUM,world); + if (flagall) return 0; -void Atom::set_shape(double **values) -{ - for (int itype = 1; itype <= ntypes; itype++) { - shape[itype][0] = values[itype][0]; - shape[itype][1] = values[itype][1]; - shape[itype][2] = values[itype][2]; - shape_setflag[itype] = 1; - } -} - -/* ---------------------------------------------------------------------- - check that all particle shapes have been set -------------------------------------------------------------------------- */ - -void Atom::check_shape() -{ - if (shape == NULL) return; - for (int itype = 1; itype <= ntypes; itype++) - if (shape_setflag[itype] == 0) error->all("All shapes are not set"); -} - -/* ---------------------------------------------------------------------- - set a dipole moment and flag it as set - called from reading of data file -------------------------------------------------------------------------- */ - -void Atom::set_dipole(const char *str) -{ - if (dipole == NULL) error->all("Cannot set dipole for this atom style"); - - int itype; - double dipole_one; - int n = sscanf(str,"%d %lg",&itype,&dipole_one); - if (n != 2) error->all("Invalid dipole line in data file"); - - dipole[itype] = dipole_one; - dipole_setflag[itype] = 1; - - if (dipole[itype] < 0.0) error->all("Invalid dipole value"); -} - -/* ---------------------------------------------------------------------- - set one or more dipole moments and flag them as set - called from reading of input script -------------------------------------------------------------------------- */ - -void Atom::set_dipole(int narg, char **arg) -{ - if (dipole == NULL) error->all("Cannot set dipole for this atom style"); - - int lo,hi; - force->bounds(arg[0],ntypes,lo,hi); - if (lo < 1 || hi > ntypes) error->all("Invalid type for dipole set"); - - for (int itype = lo; itype <= hi; itype++) { - dipole[itype] = atof(arg[1]); - dipole_setflag[itype] = 1; - - if (dipole[itype] < 0.0) error->all("Invalid dipole value"); - } -} - -/* ---------------------------------------------------------------------- - set all dipole moments as read in from restart file -------------------------------------------------------------------------- */ - -void Atom::set_dipole(double *values) -{ - for (int itype = 1; itype <= ntypes; itype++) { - dipole[itype] = values[itype]; - dipole_setflag[itype] = 1; - } -} - -/* ---------------------------------------------------------------------- - check that all dipole moments have been set -------------------------------------------------------------------------- */ - -void Atom::check_dipole() -{ - if (dipole == NULL) return; - for (int itype = 1; itype <= ntypes; itype++) - if (dipole_setflag[itype] == 0) - error->all("All dipole moments are not set"); + double oneall[3]; + MPI_Allreduce(one,oneall,3,MPI_DOUBLE,MPI_MAX,world); + shapex = oneall[0]; + shapey = oneall[1]; + shapez = oneall[2]; + return 1; } /* ---------------------------------------------------------------------- diff --git a/src/atom.h b/src/atom.h index 9ffe4c1774..e05929eaa5 100644 --- a/src/atom.h +++ b/src/atom.h @@ -48,8 +48,8 @@ class Atom : protected Pointers { int *molecule; double *q,**mu; - double **quat,**omega,**angmom,**torque; - double *radius,*density,*rmass,*vfrac,*s0; + double **quat,**omega,**angmom,**torque,**shape; + double *radius,*rmass,*vfrac,*s0; double **x0; int *spin; @@ -75,15 +75,15 @@ class Atom : protected Pointers { int **improper_type; int **improper_atom1,**improper_atom2,**improper_atom3,**improper_atom4; - // per-atom array existence flags - // these can be checked before array is allocated + // atom style and per-atom array existence flags // customize by adding new flag - int molecule_flag; - int q_flag,mu_flag; - int quat_flag,omega_flag,angmom_flag,torque_flag; - int radius_flag,density_flag,rmass_flag,vfrac_flag; - int spin_flag,eradius_flag,ervel_flag,erforce_flag; + int sphere_flag,ellipsoid_flag,peri_flag,dipole_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 vfrac_flag,spin_flag,eradius_flag,ervel_flag,erforce_flag; // extra peratom info in restart file destined for fix & diag @@ -91,8 +91,8 @@ class Atom : protected Pointers { // per-type arrays - double *mass,**shape,*dipole; - int *mass_setflag,*shape_setflag,*dipole_setflag; + double *mass; + int *mass_setflag; // callback ptrs for atom arrays managed by fix classes @@ -141,14 +141,9 @@ class Atom : protected Pointers { void set_mass(int, char **); void set_mass(double *); void check_mass(); - void set_shape(const char *); - void set_shape(int, char **); - void set_shape(double **); - void check_shape(); - void set_dipole(const char *); - void set_dipole(int, char **); - void set_dipole(double *); - void check_dipole(); + + int radius_consistency(int, double &); + int shape_consistency(int, double &, double &, double &); void first_reorder(); void sort(); diff --git a/src/atom_vec.cpp b/src/atom_vec.cpp index f43e3bcdf8..24037a2c7a 100644 --- a/src/atom_vec.cpp +++ b/src/atom_vec.cpp @@ -24,7 +24,7 @@ AtomVec::AtomVec(LAMMPS *lmp, int narg, char **arg) : Pointers(lmp) { nmax = 0; bonds_allow = angles_allow = dihedrals_allow = impropers_allow = 0; - mass_type = shape_type = dipole_type = 0; + mass_type = dipole_type = 0; } /* ---------------------------------------------------------------------- diff --git a/src/atom_vec.h b/src/atom_vec.h index 9f57ebacf9..7c399bf44b 100644 --- a/src/atom_vec.h +++ b/src/atom_vec.h @@ -24,7 +24,6 @@ class AtomVec : protected Pointers { int bonds_allow,angles_allow; // 1 if bonds, angles are used int dihedrals_allow,impropers_allow; // 1 if dihedrals, impropers used int mass_type; // 1 if per-type masses - int shape_type; // 1 if per-type shape array int dipole_type; // 1 if per-type dipole moments int comm_x_only; // 1 if only exchange x in forward comm diff --git a/src/atom_vec_hybrid.cpp b/src/atom_vec_hybrid.cpp index 0a64c99ba2..f630b67557 100644 --- a/src/atom_vec_hybrid.cpp +++ b/src/atom_vec_hybrid.cpp @@ -75,7 +75,6 @@ AtomVecHybrid::AtomVecHybrid(LAMMPS *lmp, int narg, char **arg) : dihedrals_allow = MAX(dihedrals_allow,styles[k]->dihedrals_allow); impropers_allow = MAX(impropers_allow,styles[k]->impropers_allow); mass_type = MAX(mass_type,styles[k]->mass_type); - shape_type = MAX(shape_type,styles[k]->shape_type); dipole_type = MAX(dipole_type,styles[k]->dipole_type); comm_x_only = MIN(comm_x_only,styles[k]->comm_x_only); diff --git a/src/GRANULAR/atom_vec_granular.cpp b/src/atom_vec_sphere.cpp similarity index 88% rename from src/GRANULAR/atom_vec_granular.cpp rename to src/atom_vec_sphere.cpp index e934796f61..6b7ff6eb18 100644 --- a/src/GRANULAR/atom_vec_granular.cpp +++ b/src/atom_vec_sphere.cpp @@ -15,7 +15,7 @@ #include "math.h" #include "stdlib.h" #include "string.h" -#include "atom_vec_granular.h" +#include "atom_vec_sphere.h" #include "atom.h" #include "domain.h" #include "modify.h" @@ -31,7 +31,7 @@ using namespace LAMMPS_NS; /* ---------------------------------------------------------------------- */ -AtomVecGranular::AtomVecGranular(LAMMPS *lmp, int narg, char **arg) : +AtomVecSphere::AtomVecSphere(LAMMPS *lmp, int narg, char **arg) : AtomVec(lmp, narg, arg) { molecular = 0; @@ -46,26 +46,31 @@ AtomVecGranular::AtomVecGranular(LAMMPS *lmp, int narg, char **arg) : size_data_vel = 7; xcol_data = 5; - atom->radius_flag = atom->density_flag = atom->rmass_flag = 1; - atom->omega_flag = atom->torque_flag = 1; + atom->sphere_flag = 1; + atom->radius_flag = atom->rmass_flag = atom->omega_flag = + atom->torque_flag = 1; PI = 4.0*atan(1.0); } /* ---------------------------------------------------------------------- */ -void AtomVecGranular::init() +void AtomVecSphere::init() { AtomVec::init(); // set radvary if particle diameters are time-varying due to fix adapt radvary = 0; + comm_x_only = 1; + size_forward = 3; + for (int i = 0; i < modify->nfix; i++) if (strcmp(modify->fix[i]->style,"adapt") == 0) { FixAdapt *fix = (FixAdapt *) modify->fix[i]; if (fix->diamflag) { radvary = 1; + comm_x_only = 0; size_forward = 5; } } @@ -77,7 +82,7 @@ void AtomVecGranular::init() n > 0 allocates arrays to size n ------------------------------------------------------------------------- */ -void AtomVecGranular::grow(int n) +void AtomVecSphere::grow(int n) { if (n == 0) nmax += DELTA; else nmax = n; @@ -94,9 +99,7 @@ void AtomVecGranular::grow(int n) f = memory->grow(atom->f,nmax,3,"atom:f"); radius = memory->grow(atom->radius,nmax,"atom:radius"); - density = memory->grow(atom->density,nmax,"atom:density"); rmass = memory->grow(atom->rmass,nmax,"atom:rmass"); - omega = memory->grow(atom->omega,nmax,3,"atom:omega"); torque = memory->grow(atom->torque,nmax,3,"atom:torque"); @@ -109,18 +112,18 @@ void AtomVecGranular::grow(int n) reset local array ptrs ------------------------------------------------------------------------- */ -void AtomVecGranular::grow_reset() +void AtomVecSphere::grow_reset() { tag = atom->tag; type = atom->type; mask = atom->mask; image = atom->image; x = atom->x; v = atom->v; f = atom->f; - radius = atom->radius; density = atom->density; rmass = atom->rmass; + radius = atom->radius; rmass = atom->rmass; omega = atom->omega; torque = atom->torque; } /* ---------------------------------------------------------------------- */ -void AtomVecGranular::copy(int i, int j) +void AtomVecSphere::copy(int i, int j) { tag[j] = tag[i]; type[j] = type[i]; @@ -134,7 +137,6 @@ void AtomVecGranular::copy(int i, int j) v[j][2] = v[i][2]; radius[j] = radius[i]; - density[j] = density[i]; rmass[j] = rmass[i]; omega[j][0] = omega[i][0]; omega[j][1] = omega[i][1]; @@ -147,7 +149,7 @@ void AtomVecGranular::copy(int i, int j) /* ---------------------------------------------------------------------- */ -int AtomVecGranular::pack_comm(int n, int *list, double *buf, +int AtomVecSphere::pack_comm(int n, int *list, double *buf, int pbc_flag, int *pbc) { int i,j,m; @@ -217,7 +219,7 @@ int AtomVecGranular::pack_comm(int n, int *list, double *buf, /* ---------------------------------------------------------------------- */ -int AtomVecGranular::pack_comm_vel(int n, int *list, double *buf, +int AtomVecSphere::pack_comm_vel(int n, int *list, double *buf, int pbc_flag, int *pbc) { int i,j,m; @@ -361,7 +363,7 @@ int AtomVecGranular::pack_comm_vel(int n, int *list, double *buf, /* ---------------------------------------------------------------------- */ -int AtomVecGranular::pack_comm_one(int i, double *buf) +int AtomVecSphere::pack_comm_one(int i, double *buf) { if (radvary == 0) return 0; @@ -372,7 +374,7 @@ int AtomVecGranular::pack_comm_one(int i, double *buf) /* ---------------------------------------------------------------------- */ -void AtomVecGranular::unpack_comm(int n, int first, double *buf) +void AtomVecSphere::unpack_comm(int n, int first, double *buf) { int i,m,last; @@ -399,7 +401,7 @@ void AtomVecGranular::unpack_comm(int n, int first, double *buf) /* ---------------------------------------------------------------------- */ -void AtomVecGranular::unpack_comm_vel(int n, int first, double *buf) +void AtomVecSphere::unpack_comm_vel(int n, int first, double *buf) { int i,m,last; @@ -438,7 +440,7 @@ void AtomVecGranular::unpack_comm_vel(int n, int first, double *buf) /* ---------------------------------------------------------------------- */ -int AtomVecGranular::unpack_comm_one(int i, double *buf) +int AtomVecSphere::unpack_comm_one(int i, double *buf) { if (radvary == 0) return 0; @@ -449,7 +451,7 @@ int AtomVecGranular::unpack_comm_one(int i, double *buf) /* ---------------------------------------------------------------------- */ -int AtomVecGranular::pack_reverse(int n, int first, double *buf) +int AtomVecSphere::pack_reverse(int n, int first, double *buf) { int i,m,last; @@ -468,7 +470,7 @@ int AtomVecGranular::pack_reverse(int n, int first, double *buf) /* ---------------------------------------------------------------------- */ -int AtomVecGranular::pack_reverse_one(int i, double *buf) +int AtomVecSphere::pack_reverse_one(int i, double *buf) { buf[0] = torque[i][0]; buf[1] = torque[i][1]; @@ -478,7 +480,7 @@ int AtomVecGranular::pack_reverse_one(int i, double *buf) /* ---------------------------------------------------------------------- */ -void AtomVecGranular::unpack_reverse(int n, int *list, double *buf) +void AtomVecSphere::unpack_reverse(int n, int *list, double *buf) { int i,j,m; @@ -496,7 +498,7 @@ void AtomVecGranular::unpack_reverse(int n, int *list, double *buf) /* ---------------------------------------------------------------------- */ -int AtomVecGranular::unpack_reverse_one(int i, double *buf) +int AtomVecSphere::unpack_reverse_one(int i, double *buf) { torque[i][0] += buf[0]; torque[i][1] += buf[1]; @@ -506,7 +508,7 @@ int AtomVecGranular::unpack_reverse_one(int i, double *buf) /* ---------------------------------------------------------------------- */ -int AtomVecGranular::pack_border(int n, int *list, double *buf, +int AtomVecSphere::pack_border(int n, int *list, double *buf, int pbc_flag, int *pbc) { int i,j,m; @@ -552,7 +554,7 @@ int AtomVecGranular::pack_border(int n, int *list, double *buf, /* ---------------------------------------------------------------------- */ -int AtomVecGranular::pack_border_vel(int n, int *list, double *buf, +int AtomVecSphere::pack_border_vel(int n, int *list, double *buf, int pbc_flag, int *pbc) { int i,j,m; @@ -639,7 +641,7 @@ int AtomVecGranular::pack_border_vel(int n, int *list, double *buf, /* ---------------------------------------------------------------------- */ -int AtomVecGranular::pack_border_one(int i, double *buf) +int AtomVecSphere::pack_border_one(int i, double *buf) { buf[0] = radius[i]; buf[1] = rmass[i]; @@ -648,7 +650,7 @@ int AtomVecGranular::pack_border_one(int i, double *buf) /* ---------------------------------------------------------------------- */ -void AtomVecGranular::unpack_border(int n, int first, double *buf) +void AtomVecSphere::unpack_border(int n, int first, double *buf) { int i,m,last; @@ -670,7 +672,7 @@ void AtomVecGranular::unpack_border(int n, int first, double *buf) /* ---------------------------------------------------------------------- */ -void AtomVecGranular::unpack_border_vel(int n, int first, double *buf) +void AtomVecSphere::unpack_border_vel(int n, int first, double *buf) { int i,m,last; @@ -697,7 +699,7 @@ void AtomVecGranular::unpack_border_vel(int n, int first, double *buf) /* ---------------------------------------------------------------------- */ -int AtomVecGranular::unpack_border_one(int i, double *buf) +int AtomVecSphere::unpack_border_one(int i, double *buf) { radius[i] = buf[0]; rmass[i] = buf[1]; @@ -709,7 +711,7 @@ int AtomVecGranular::unpack_border_one(int i, double *buf) xyz must be 1st 3 values, so comm::exchange() can test on them ------------------------------------------------------------------------- */ -int AtomVecGranular::pack_exchange(int i, double *buf) +int AtomVecSphere::pack_exchange(int i, double *buf) { int m = 1; buf[m++] = x[i][0]; @@ -724,7 +726,6 @@ int AtomVecGranular::pack_exchange(int i, double *buf) buf[m++] = image[i]; buf[m++] = radius[i]; - buf[m++] = density[i]; buf[m++] = rmass[i]; buf[m++] = omega[i][0]; buf[m++] = omega[i][1]; @@ -740,7 +741,7 @@ int AtomVecGranular::pack_exchange(int i, double *buf) /* ---------------------------------------------------------------------- */ -int AtomVecGranular::unpack_exchange(double *buf) +int AtomVecSphere::unpack_exchange(double *buf) { int nlocal = atom->nlocal; if (nlocal == nmax) grow(0); @@ -758,7 +759,6 @@ int AtomVecGranular::unpack_exchange(double *buf) image[nlocal] = static_cast (buf[m++]); radius[nlocal] = buf[m++]; - density[nlocal] = buf[m++]; rmass[nlocal] = buf[m++]; omega[nlocal][0] = buf[m++]; omega[nlocal][1] = buf[m++]; @@ -778,7 +778,7 @@ int AtomVecGranular::unpack_exchange(double *buf) include extra data stored by fixes ------------------------------------------------------------------------- */ -int AtomVecGranular::size_restart() +int AtomVecSphere::size_restart() { int i; @@ -799,7 +799,7 @@ int AtomVecGranular::size_restart() molecular types may be negative, but write as positive ------------------------------------------------------------------------- */ -int AtomVecGranular::pack_restart(int i, double *buf) +int AtomVecSphere::pack_restart(int i, double *buf) { int m = 1; buf[m++] = x[i][0]; @@ -814,7 +814,7 @@ int AtomVecGranular::pack_restart(int i, double *buf) buf[m++] = v[i][2]; buf[m++] = radius[i]; - buf[m++] = density[i]; + buf[m++] = rmass[i]; buf[m++] = omega[i][0]; buf[m++] = omega[i][1]; buf[m++] = omega[i][2]; @@ -831,7 +831,7 @@ int AtomVecGranular::pack_restart(int i, double *buf) unpack data for one atom from restart file including extra quantities ------------------------------------------------------------------------- */ -int AtomVecGranular::unpack_restart(double *buf) +int AtomVecSphere::unpack_restart(double *buf) { int nlocal = atom->nlocal; if (nlocal == nmax) { @@ -853,13 +853,7 @@ int AtomVecGranular::unpack_restart(double *buf) v[nlocal][2] = buf[m++]; radius[nlocal] = buf[m++]; - density[nlocal] = buf[m++]; - - if (radius[nlocal] == 0.0) rmass[nlocal] = density[nlocal]; - else - rmass[nlocal] = 4.0*PI/3.0 * - radius[nlocal]*radius[nlocal]*radius[nlocal] * density[nlocal]; - + rmass[nlocal] = buf[m++]; omega[nlocal][0] = buf[m++]; omega[nlocal][1] = buf[m++]; omega[nlocal][2] = buf[m++]; @@ -879,7 +873,7 @@ int AtomVecGranular::unpack_restart(double *buf) set other values to defaults ------------------------------------------------------------------------- */ -void AtomVecGranular::create_atom(int itype, double *coord) +void AtomVecSphere::create_atom(int itype, double *coord) { int nlocal = atom->nlocal; if (nlocal == nmax) grow(0); @@ -896,9 +890,7 @@ void AtomVecGranular::create_atom(int itype, double *coord) v[nlocal][2] = 0.0; radius[nlocal] = 0.5; - density[nlocal] = 1.0; - rmass[nlocal] = 4.0*PI/3.0 * - radius[nlocal]*radius[nlocal]*radius[nlocal] * density[nlocal]; + rmass[nlocal] = 4.0*PI/3.0 * radius[nlocal]*radius[nlocal]*radius[nlocal]; omega[nlocal][0] = 0.0; omega[nlocal][1] = 0.0; omega[nlocal][2] = 0.0; @@ -911,7 +903,7 @@ void AtomVecGranular::create_atom(int itype, double *coord) initialize other atom quantities ------------------------------------------------------------------------- */ -void AtomVecGranular::data_atom(double *coord, int imagetmp, char **values) +void AtomVecSphere::data_atom(double *coord, int imagetmp, char **values) { int nlocal = atom->nlocal; if (nlocal == nmax) grow(0); @@ -928,14 +920,14 @@ void AtomVecGranular::data_atom(double *coord, int imagetmp, char **values) if (radius[nlocal] < 0.0) error->one("Invalid radius in Atoms section of data file"); - density[nlocal] = atof(values[3]); - if (density[nlocal] <= 0.0) + double density = atof(values[3]); + if (density <= 0.0) error->one("Invalid density in Atoms section of data file"); - if (radius[nlocal] == 0.0) rmass[nlocal] = density[nlocal]; + if (radius[nlocal] == 0.0) rmass[nlocal] = density; else rmass[nlocal] = 4.0*PI/3.0 * - radius[nlocal]*radius[nlocal]*radius[nlocal] * density[nlocal]; + radius[nlocal]*radius[nlocal]*radius[nlocal] * density; x[nlocal][0] = coord[0]; x[nlocal][1] = coord[1]; @@ -959,20 +951,20 @@ void AtomVecGranular::data_atom(double *coord, int imagetmp, char **values) initialize other atom quantities for this sub-style ------------------------------------------------------------------------- */ -int AtomVecGranular::data_atom_hybrid(int nlocal, char **values) +int AtomVecSphere::data_atom_hybrid(int nlocal, char **values) { radius[nlocal] = 0.5 * atof(values[0]); if (radius[nlocal] < 0.0) error->one("Invalid radius in Atoms section of data file"); - density[nlocal] = atof(values[1]); - if (density[nlocal] <= 0.0) + double density = atof(values[1]); + if (density <= 0.0) error->one("Invalid density in Atoms section of data file"); - if (radius[nlocal] == 0.0) rmass[nlocal] = density[nlocal]; + if (radius[nlocal] == 0.0) rmass[nlocal] = density; else rmass[nlocal] = 4.0*PI/3.0 * - radius[nlocal]*radius[nlocal]*radius[nlocal] * density[nlocal]; + radius[nlocal]*radius[nlocal]*radius[nlocal] * density; return 2; } @@ -981,7 +973,7 @@ int AtomVecGranular::data_atom_hybrid(int nlocal, char **values) unpack one line from Velocities section of data file ------------------------------------------------------------------------- */ -void AtomVecGranular::data_vel(int m, char **values) +void AtomVecSphere::data_vel(int m, char **values) { v[m][0] = atof(values[0]); v[m][1] = atof(values[1]); @@ -995,7 +987,7 @@ void AtomVecGranular::data_vel(int m, char **values) unpack hybrid quantities from one line in Velocities section of data file ------------------------------------------------------------------------- */ -int AtomVecGranular::data_vel_hybrid(int m, char **values) +int AtomVecSphere::data_vel_hybrid(int m, char **values) { omega[m][0] = atof(values[0]); omega[m][1] = atof(values[1]); @@ -1007,7 +999,7 @@ int AtomVecGranular::data_vel_hybrid(int m, char **values) return # of bytes of allocated memory ------------------------------------------------------------------------- */ -bigint AtomVecGranular::memory_usage() +bigint AtomVecSphere::memory_usage() { bigint bytes = 0; @@ -1020,7 +1012,6 @@ bigint AtomVecGranular::memory_usage() if (atom->memcheck("f")) bytes += memory->usage(f,nmax,3); if (atom->memcheck("radius")) bytes += memory->usage(radius,nmax); - if (atom->memcheck("density")) bytes += memory->usage(density,nmax); if (atom->memcheck("rmass")) bytes += memory->usage(rmass,nmax); if (atom->memcheck("omega")) bytes += memory->usage(omega,nmax,3); if (atom->memcheck("torque")) bytes += memory->usage(torque,nmax,3); diff --git a/src/GRANULAR/atom_vec_granular.h b/src/atom_vec_sphere.h similarity index 90% rename from src/GRANULAR/atom_vec_granular.h rename to src/atom_vec_sphere.h index 9e18ac8038..d1fa9ba42b 100644 --- a/src/GRANULAR/atom_vec_granular.h +++ b/src/atom_vec_sphere.h @@ -13,21 +13,21 @@ #ifdef ATOM_CLASS -AtomStyle(granular,AtomVecGranular) +AtomStyle(sphere,AtomVecSphere) #else -#ifndef LMP_ATOM_VEC_GRANULAR_H -#define LMP_ATOM_VEC_GRANULAR_H +#ifndef LMP_ATOM_VEC_SPHERE_H +#define LMP_ATOM_VEC_SPHERE_H #include "atom_vec.h" namespace LAMMPS_NS { -class AtomVecGranular : public AtomVec { +class AtomVecSphere : public AtomVec { public: - AtomVecGranular(class LAMMPS *, int, char **); - ~AtomVecGranular() {} + AtomVecSphere(class LAMMPS *, int, char **); + ~AtomVecSphere() {} void init(); void grow(int); void grow_reset(); diff --git a/src/compute_erotate_sphere.cpp b/src/compute_erotate_sphere.cpp index 1558196cce..b357501a06 100644 --- a/src/compute_erotate_sphere.cpp +++ b/src/compute_erotate_sphere.cpp @@ -35,40 +35,16 @@ ComputeERotateSphere::ComputeERotateSphere(LAMMPS *lmp, int narg, char **arg) : scalar_flag = 1; extscalar = 1; - // error checks + // error check - if (!atom->omega_flag) - error->all("Compute erotate/sphere requires atom attribute omega"); - if (!atom->radius_flag && !atom->avec->shape_type) - error->all("Compute erotate/sphere requires atom attribute " - "radius or shape"); + if (!atom->sphere_flag) + error->all("Compute erotate/sphere requires atom style sphere"); } /* ---------------------------------------------------------------------- */ void ComputeERotateSphere::init() { - int i,itype; - - // if shape used, check that all particles are spherical - // point particles are allowed - - if (atom->radius == NULL) { - double **shape = atom->shape; - int *type = atom->type; - int *mask = atom->mask; - int nlocal = atom->nlocal; - - for (i = 0; i < nlocal; i++) - if (mask[i] & groupbit) { - itype = type[i]; - if (shape[itype][0] != shape[itype][1] || - shape[itype][0] != shape[itype][2]) - error->one("Compute erotate/sphere requires " - "spherical particle shapes"); - } - } - pfactor = 0.5 * force->mvv2e * INERTIA; } @@ -76,59 +52,22 @@ void ComputeERotateSphere::init() double ComputeERotateSphere::compute_scalar() { - int i,itype; - invoked_scalar = update->ntimestep; double **omega = atom->omega; double *radius = atom->radius; double *rmass = atom->rmass; - double *mass = atom->mass; - double **shape = atom->shape; int *mask = atom->mask; - int *type = atom->type; int nlocal = atom->nlocal; // sum rotational energy for each particle - // point particles will not contribute due to radius or shape = 0 + // point particles will not contribute due to radius = 0 double erotate = 0.0; - - if (radius) { - if (rmass) { - for (i = 0; i < nlocal; i++) - if (mask[i] & groupbit) - erotate += (omega[i][0]*omega[i][0] + omega[i][1]*omega[i][1] + - omega[i][2]*omega[i][2]) * radius[i]*radius[i]*rmass[i]; - } else { - for (i = 0; i < nlocal; i++) - if (mask[i] & groupbit) { - itype = type[i]; - erotate += (omega[i][0]*omega[i][0] + omega[i][1]*omega[i][1] + - omega[i][2]*omega[i][2]) * - radius[i]*radius[i]*mass[itype]; - } - } - - } else { - if (rmass) { - for (i = 0; i < nlocal; i++) - if (mask[i] & groupbit) { - itype = type[i]; - erotate += (omega[i][0]*omega[i][0] + omega[i][1]*omega[i][1] + - omega[i][2]*omega[i][2]) * - shape[itype][0]*shape[itype][0]*rmass[i]; - } - } else { - for (i = 0; i < nlocal; i++) - if (mask[i] & groupbit) { - itype = type[i]; - erotate += (omega[i][0]*omega[i][0] + omega[i][1]*omega[i][1] + - omega[i][2]*omega[i][2]) * - shape[itype][0]*shape[itype][0]*mass[itype]; - } - } - } + for (int i = 0; i < nlocal; i++) + if (mask[i] & groupbit) + erotate += (omega[i][0]*omega[i][0] + omega[i][1]*omega[i][1] + + omega[i][2]*omega[i][2]) * radius[i]*radius[i]*rmass[i]; MPI_Allreduce(&erotate,&scalar,1,MPI_DOUBLE,MPI_SUM,world); scalar *= pfactor; diff --git a/src/compute_property_atom.cpp b/src/compute_property_atom.cpp index 1883703c06..2949b9f469 100644 --- a/src/compute_property_atom.cpp +++ b/src/compute_property_atom.cpp @@ -124,6 +124,11 @@ ComputePropertyAtom::ComputePropertyAtom(LAMMPS *lmp, int narg, char **arg) : error->all("Compute property/atom for " "atom property that isn't allocated"); pack_choice[i] = &ComputePropertyAtom::pack_muz; + } else if (strcmp(arg[iarg],"mu") == 0) { + if (!atom->mu_flag) + error->all("Compute property/atom for " + "atom property that isn't allocated"); + pack_choice[i] = &ComputePropertyAtom::pack_mu; } else if (strcmp(arg[iarg],"radius") == 0) { if (!atom->radius_flag) @@ -161,6 +166,21 @@ ComputePropertyAtom::ComputePropertyAtom(LAMMPS *lmp, int narg, char **arg) : "atom property that isn't allocated"); 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"); + 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"); + 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"); + pack_choice[i] = &ComputePropertyAtom::pack_shapez; } else if (strcmp(arg[iarg],"quatw") == 0) { if (!atom->quat_flag) error->all("Compute property/atom for " @@ -831,6 +851,21 @@ void ComputePropertyAtom::pack_muz(int n) /* ---------------------------------------------------------------------- */ +void ComputePropertyAtom::pack_mu(int n) +{ + double **mu = atom->mu; + int *mask = atom->mask; + int nlocal = atom->nlocal; + + for (int i = 0; i < nlocal; i++) { + if (mask[i] & groupbit) buf[n] = mu[i][3]; + else buf[n] = 0.0; + n += nvalues; + } +} + +/* ---------------------------------------------------------------------- */ + void ComputePropertyAtom::pack_radius(int n) { double *radius = atom->radius; @@ -936,6 +971,51 @@ void ComputePropertyAtom::pack_angmomz(int n) /* ---------------------------------------------------------------------- */ +void ComputePropertyAtom::pack_shapex(int n) +{ + double **shape = atom->shape; + int *mask = atom->mask; + int nlocal = atom->nlocal; + + for (int i = 0; i < nlocal; i++) { + if (mask[i] & groupbit) buf[n] = shape[i][0]; + else buf[n] = 0.0; + n += nvalues; + } +} + +/* ---------------------------------------------------------------------- */ + +void ComputePropertyAtom::pack_shapey(int n) +{ + double **shape = atom->shape; + int *mask = atom->mask; + int nlocal = atom->nlocal; + + for (int i = 0; i < nlocal; i++) { + if (mask[i] & groupbit) buf[n] = shape[i][1]; + else buf[n] = 0.0; + n += nvalues; + } +} + +/* ---------------------------------------------------------------------- */ + +void ComputePropertyAtom::pack_shapez(int n) +{ + double **shape = atom->shape; + int *mask = atom->mask; + int nlocal = atom->nlocal; + + for (int i = 0; i < nlocal; i++) { + if (mask[i] & groupbit) buf[n] = shape[i][2]; + else buf[n] = 0.0; + n += nvalues; + } +} + +/* ---------------------------------------------------------------------- */ + void ComputePropertyAtom::pack_quatw(int n) { double **quat = atom->quat; diff --git a/src/compute_property_atom.h b/src/compute_property_atom.h index 1d8cd07f23..0c3bd21b32 100644 --- a/src/compute_property_atom.h +++ b/src/compute_property_atom.h @@ -76,6 +76,7 @@ class ComputePropertyAtom : public Compute { void pack_mux(int); void pack_muy(int); void pack_muz(int); + void pack_mu(int); void pack_radius(int); void pack_omegax(int); void pack_omegay(int); @@ -83,6 +84,9 @@ class ComputePropertyAtom : public Compute { 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); diff --git a/src/compute_temp_sphere.cpp b/src/compute_temp_sphere.cpp index bd8aa03362..bad55efdb8 100644 --- a/src/compute_temp_sphere.cpp +++ b/src/compute_temp_sphere.cpp @@ -55,11 +55,8 @@ ComputeTempSphere::ComputeTempSphere(LAMMPS *lmp, int narg, char **arg) : // error checks - if (!atom->omega_flag) - error->all("Compute temp/sphere requires atom attribute omega"); - if (!atom->radius_flag && !atom->avec->shape_type) - error->all("Compute temp/sphere requires atom attribute " - "radius or shape"); + if (!atom->sphere_flag) + error->all("Compute temp/sphere requires atom style sphere"); } /* ---------------------------------------------------------------------- */ @@ -74,29 +71,8 @@ ComputeTempSphere::~ComputeTempSphere() void ComputeTempSphere::init() { - int i,itype; - - // if shape used, check that all particles are spherical - // point particles are allowed - - if (atom->radius == NULL) { - double **shape = atom->shape; - int *type = atom->type; - int *mask = atom->mask; - int nlocal = atom->nlocal; - - for (i = 0; i < nlocal; i++) - if (mask[i] & groupbit) { - itype = type[i]; - if (shape[itype][0] != shape[itype][1] || - shape[itype][0] != shape[itype][2]) - error->one("Compute temp/sphere requires " - "spherical particle shapes"); - } - } - if (tempbias) { - i = modify->find_compute(id_bias); + int i = modify->find_compute(id_bias); if (i < 0) error->all("Could not find compute ID for temperature bias"); tbias = modify->compute[i]; if (tbias->tempflag == 0) @@ -111,7 +87,7 @@ void ComputeTempSphere::init() } fix_dof = 0; - for (i = 0; i < modify->nfix; i++) + for (int i = 0; i < modify->nfix; i++) fix_dof += modify->fix[i]->dof(igroup); dof_compute(); } @@ -130,44 +106,22 @@ void ComputeTempSphere::dof_compute() int dimension = domain->dimension; double *radius = atom->radius; - double **shape = atom->shape; - int *type = atom->type; int *mask = atom->mask; int nlocal = atom->nlocal; count = 0; if (dimension == 3) { - if (radius) { - for (int i = 0; i < nlocal; i++) { - if (mask[i] & groupbit) { - if (radius[i] == 0.0) count += 3; - else count += 6; - } + for (int i = 0; i < nlocal; i++) + if (mask[i] & groupbit) { + if (radius[i] == 0.0) count += 3; + else count += 6; } - } else { - for (int i = 0; i < nlocal; i++) { - if (mask[i] & groupbit) { - if (shape[type[i]][0] == 0.0) count += 3; - else count += 6; - } - } - } } else { - if (radius) { - for (int i = 0; i < nlocal; i++) { - if (mask[i] & groupbit) { - if (radius[i] == 0.0) count += 2; - else count += 3; - } + for (int i = 0; i < nlocal; i++) + if (mask[i] & groupbit) { + if (radius[i] == 0.0) count += 2; + else count += 3; } - } else { - for (int i = 0; i < nlocal; i++) { - if (mask[i] & groupbit) { - if (shape[type[i]][0] == 0.0) count += 2; - else count += 3; - } - } - } } MPI_Allreduce(&count,&count_all,1,MPI_INT,MPI_SUM,world); @@ -185,41 +139,21 @@ void ComputeTempSphere::dof_compute() count = 0; if (dimension == 3) { - if (radius) { - for (int i = 0; i < nlocal; i++) - if (mask[i] & groupbit) { - if (tbias->dof_remove(i)) { - if (radius[i] == 0.0) count += 3; - else count += 6; - } + for (int i = 0; i < nlocal; i++) + if (mask[i] & groupbit) { + if (tbias->dof_remove(i)) { + if (radius[i] == 0.0) count += 3; + else count += 6; } - } else { - for (int i = 0; i < nlocal; i++) - if (mask[i] & groupbit) { - if (tbias->dof_remove(i)) { - if (shape[type[i]][0] == 0.0) count += 3; - else count += 6; - } - } - } + } } else { - if (radius) { - for (int i = 0; i < nlocal; i++) - if (mask[i] & groupbit) { - if (tbias->dof_remove(i)) { - if (radius[i] == 0.0) count += 2; - else count += 3; - } + for (int i = 0; i < nlocal; i++) + if (mask[i] & groupbit) { + if (tbias->dof_remove(i)) { + if (radius[i] == 0.0) count += 2; + else count += 3; } - } else { - for (int i = 0; i < nlocal; i++) - if (mask[i] & groupbit) { - if (tbias->dof_remove(i)) { - if (shape[type[i]][0] == 0.0) count += 2; - else count += 3; - } - } - } + } } MPI_Allreduce(&count,&count_all,1,MPI_INT,MPI_SUM,world); @@ -235,8 +169,6 @@ void ComputeTempSphere::dof_compute() double ComputeTempSphere::compute_scalar() { - int i,itype; - invoked_scalar = update->ntimestep; if (tempbias) { @@ -248,65 +180,20 @@ double ComputeTempSphere::compute_scalar() double **omega = atom->omega; double *radius = atom->radius; double *rmass = atom->rmass; - double *mass = atom->mass; - double **shape = atom->shape; - int *type = atom->type; int *mask = atom->mask; int nlocal = atom->nlocal; - // 4 cases depending on radius vs shape and rmass vs mass - // point particles will not contribute rotation due to radius or shape = 0 + // point particles will not contribute rotation due to radius = 0 double t = 0.0; - if (radius) { - if (rmass) { - for (i = 0; i < nlocal; i++) - if (mask[i] & groupbit) { - t += (v[i][0]*v[i][0] + v[i][1]*v[i][1] + v[i][2]*v[i][2]) * - rmass[i]; - t += (omega[i][0]*omega[i][0] + omega[i][1]*omega[i][1] + - omega[i][2]*omega[i][2]) * - INERTIA*radius[i]*radius[i]*rmass[i]; - } - - } else { - for (i = 0; i < nlocal; i++) - if (mask[i] & groupbit) { - itype = type[i]; - t += (v[i][0]*v[i][0] + v[i][1]*v[i][1] + v[i][2]*v[i][2]) * - mass[itype]; - t += (omega[i][0]*omega[i][0] + omega[i][1]*omega[i][1] + - omega[i][2]*omega[i][2]) * - INERTIA*radius[i]*radius[i]*mass[itype]; - } + for (int i = 0; i < nlocal; i++) + if (mask[i] & groupbit) { + t += (v[i][0]*v[i][0] + v[i][1]*v[i][1] + v[i][2]*v[i][2]) * rmass[i]; + t += (omega[i][0]*omega[i][0] + omega[i][1]*omega[i][1] + + omega[i][2]*omega[i][2]) * INERTIA*radius[i]*radius[i]*rmass[i]; } - } else { - if (rmass) { - for (i = 0; i < nlocal; i++) - if (mask[i] & groupbit) { - itype = type[i]; - t += (v[i][0]*v[i][0] + v[i][1]*v[i][1] + v[i][2]*v[i][2]) * - rmass[i]; - t += (omega[i][0]*omega[i][0] + omega[i][1]*omega[i][1] + - omega[i][2]*omega[i][2]) * - INERTIA*shape[itype][0]*shape[itype][0]*rmass[i]; - } - - } else { - for (i = 0; i < nlocal; i++) - if (mask[i] & groupbit) { - itype = type[i]; - t += (v[i][0]*v[i][0] + v[i][1]*v[i][1] + v[i][2]*v[i][2]) * - mass[itype]; - t += (omega[i][0]*omega[i][0] + omega[i][1]*omega[i][1] + - omega[i][2]*omega[i][2]) * - INERTIA*shape[itype][0]*shape[itype][0]*mass[itype]; - } - } - } - if (tempbias) tbias->restore_bias_all(); MPI_Allreduce(&t,&scalar,1,MPI_DOUBLE,MPI_SUM,world); @@ -319,8 +206,6 @@ double ComputeTempSphere::compute_scalar() void ComputeTempSphere::compute_vector() { - int i,itype; - invoked_vector = update->ntimestep; if (tempbias) { @@ -330,112 +215,39 @@ void ComputeTempSphere::compute_vector() double **v = atom->v; double **omega = atom->omega; - double *mass = atom->mass; double *rmass = atom->rmass; double *radius = atom->radius; - double **shape = atom->shape; - int *type = atom->type; int *mask = atom->mask; int nlocal = atom->nlocal; - // 4 cases depending on radius vs shape and rmass vs mass - // point particles will not contribute rotation due to radius or shape = 0 + // point particles will not contribute rotation due to radius = 0 double massone,inertiaone,t[6]; - for (i = 0; i < 6; i++) t[i] = 0.0; + for (int i = 0; i < 6; i++) t[i] = 0.0; - if (radius) { - if (rmass) { - for (i = 0; i < nlocal; i++) - if (mask[i] & groupbit) { - massone = rmass[i]; - t[0] += massone * v[i][0]*v[i][0]; - t[1] += massone * v[i][1]*v[i][1]; - t[2] += massone * v[i][2]*v[i][2]; - t[3] += massone * v[i][0]*v[i][1]; - t[4] += massone * v[i][0]*v[i][2]; - t[5] += massone * v[i][1]*v[i][2]; - - inertiaone = INERTIA*radius[i]*radius[i]*rmass[i]; - t[0] += inertiaone * omega[i][0]*omega[i][0]; - t[1] += inertiaone * omega[i][1]*omega[i][1]; - t[2] += inertiaone * omega[i][2]*omega[i][2]; - t[3] += inertiaone * omega[i][0]*omega[i][1]; - t[4] += inertiaone * omega[i][0]*omega[i][2]; - t[5] += inertiaone * omega[i][1]*omega[i][2]; - } - - } else { - for (i = 0; i < nlocal; i++) - if (mask[i] & groupbit) { - itype = type[i]; - massone = mass[itype]; - t[0] += massone * v[i][0]*v[i][0]; - t[1] += massone * v[i][1]*v[i][1]; - t[2] += massone * v[i][2]*v[i][2]; - t[3] += massone * v[i][0]*v[i][1]; - t[4] += massone * v[i][0]*v[i][2]; - t[5] += massone * v[i][1]*v[i][2]; - - inertiaone = INERTIA*radius[i]*radius[i]*mass[itype]; - t[0] += inertiaone * omega[i][0]*omega[i][0]; - t[1] += inertiaone * omega[i][1]*omega[i][1]; - t[2] += inertiaone * omega[i][2]*omega[i][2]; - t[3] += inertiaone * omega[i][0]*omega[i][1]; - t[4] += inertiaone * omega[i][0]*omega[i][2]; - t[5] += inertiaone * omega[i][1]*omega[i][2]; - } + for (int i = 0; i < nlocal; i++) + if (mask[i] & groupbit) { + massone = rmass[i]; + t[0] += massone * v[i][0]*v[i][0]; + t[1] += massone * v[i][1]*v[i][1]; + t[2] += massone * v[i][2]*v[i][2]; + t[3] += massone * v[i][0]*v[i][1]; + t[4] += massone * v[i][0]*v[i][2]; + t[5] += massone * v[i][1]*v[i][2]; + + inertiaone = INERTIA*radius[i]*radius[i]*rmass[i]; + t[0] += inertiaone * omega[i][0]*omega[i][0]; + t[1] += inertiaone * omega[i][1]*omega[i][1]; + t[2] += inertiaone * omega[i][2]*omega[i][2]; + t[3] += inertiaone * omega[i][0]*omega[i][1]; + t[4] += inertiaone * omega[i][0]*omega[i][2]; + t[5] += inertiaone * omega[i][1]*omega[i][2]; } - } else { - if (rmass) { - for (i = 0; i < nlocal; i++) - if (mask[i] & groupbit) { - itype = type[i]; - massone = rmass[i]; - t[0] += massone * v[i][0]*v[i][0]; - t[1] += massone * v[i][1]*v[i][1]; - t[2] += massone * v[i][2]*v[i][2]; - t[3] += massone * v[i][0]*v[i][1]; - t[4] += massone * v[i][0]*v[i][2]; - t[5] += massone * v[i][1]*v[i][2]; - - inertiaone = INERTIA*shape[itype][0]*shape[itype][0]*rmass[i]; - t[0] += inertiaone * omega[i][0]*omega[i][0]; - t[1] += inertiaone * omega[i][1]*omega[i][1]; - t[2] += inertiaone * omega[i][2]*omega[i][2]; - t[3] += inertiaone * omega[i][0]*omega[i][1]; - t[4] += inertiaone * omega[i][0]*omega[i][2]; - t[5] += inertiaone * omega[i][1]*omega[i][2]; - } - - } else { - for (i = 0; i < nlocal; i++) - if (mask[i] & groupbit) { - itype = type[i]; - massone = mass[itype]; - t[0] += massone * v[i][0]*v[i][0]; - t[1] += massone * v[i][1]*v[i][1]; - t[2] += massone * v[i][2]*v[i][2]; - t[3] += massone * v[i][0]*v[i][1]; - t[4] += massone * v[i][0]*v[i][2]; - t[5] += massone * v[i][1]*v[i][2]; - - inertiaone = INERTIA*shape[itype][0]*shape[itype][0]*mass[itype]; - t[0] += inertiaone * omega[i][0]*omega[i][0]; - t[1] += inertiaone * omega[i][1]*omega[i][1]; - t[2] += inertiaone * omega[i][2]*omega[i][2]; - t[3] += inertiaone * omega[i][0]*omega[i][1]; - t[4] += inertiaone * omega[i][0]*omega[i][2]; - t[5] += inertiaone * omega[i][1]*omega[i][2]; - } - } - } - if (tempbias) tbias->restore_bias_all(); MPI_Allreduce(t,vector,6,MPI_DOUBLE,MPI_SUM,world); - for (i = 0; i < 6; i++) vector[i] *= force->mvv2e; + for (int i = 0; i < 6; i++) vector[i] *= force->mvv2e; } /* ---------------------------------------------------------------------- diff --git a/src/dump_custom.cpp b/src/dump_custom.cpp index f87e53edcc..4f6edd1140 100644 --- a/src/dump_custom.cpp +++ b/src/dump_custom.cpp @@ -35,11 +35,12 @@ using namespace LAMMPS_NS; // same list as in compute_property.cpp, also customize that command 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,RADIUS,OMEGAX,OMEGAY,OMEGAZ,ANGMOMX,ANGMOMY,ANGMOMZ, - QUATW,QUATI,QUATJ,QUATK,TQX,TQY,TQZ,SPIN,ERADIUS,ERVEL,ERFORCE, - COMPUTE,FIX,VARIABLE}; + 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, + COMPUTE,FIX,VARIABLE}; enum{LT,LE,GT,GE,EQ,NEQ}; enum{INT,DOUBLE}; @@ -610,17 +611,22 @@ int DumpCustom::count() if (!atom->mu_flag) error->all("Threshhold for an atom property that isn't allocated"); ptr = &atom->mu[0][0]; - nstride = 3; + nstride = 4; } else if (thresh_array[ithresh] == MUY) { if (!atom->mu_flag) error->all("Threshhold for an atom property that isn't allocated"); ptr = &atom->mu[0][1]; - nstride = 3; + nstride = 4; } else if (thresh_array[ithresh] == MUZ) { if (!atom->mu_flag) error->all("Threshhold for an atom property that isn't allocated"); ptr = &atom->mu[0][2]; - nstride = 3; + nstride = 4; + } else if (thresh_array[ithresh] == MU) { + if (!atom->mu_flag) + error->all("Threshhold for an atom property that isn't allocated"); + ptr = &atom->mu[0][3]; + nstride = 4; } else if (thresh_array[ithresh] == RADIUS) { if (!atom->radius_flag) @@ -658,6 +664,21 @@ int DumpCustom::count() 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"); @@ -944,6 +965,11 @@ void DumpCustom::parse_fields(int narg, char **arg) error->all("Dumping an atom property that isn't allocated"); pack_choice[i] = &DumpCustom::pack_muz; vtype[i] = DOUBLE; + } else if (strcmp(arg[iarg],"mu") == 0) { + if (!atom->mu_flag) + error->all("Dumping an atom property that isn't allocated"); + pack_choice[i] = &DumpCustom::pack_mu; + vtype[i] = DOUBLE; } else if (strcmp(arg[iarg],"radius") == 0) { if (!atom->radius_flag) @@ -981,6 +1007,21 @@ void DumpCustom::parse_fields(int narg, char **arg) 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"); @@ -1300,6 +1341,7 @@ int DumpCustom::modify_param(int narg, char **arg) else if (strcmp(arg[1],"mux") == 0) thresh_array[nthresh] = MUX; else if (strcmp(arg[1],"muy") == 0) thresh_array[nthresh] = MUY; else if (strcmp(arg[1],"muz") == 0) thresh_array[nthresh] = MUZ; + else if (strcmp(arg[1],"mu") == 0) thresh_array[nthresh] = MU; else if (strcmp(arg[1],"radius") == 0) thresh_array[nthresh] = RADIUS; else if (strcmp(arg[1],"omegax") == 0) thresh_array[nthresh] = OMEGAX; @@ -1309,6 +1351,9 @@ int DumpCustom::modify_param(int narg, char **arg) 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; @@ -2039,6 +2084,20 @@ void DumpCustom::pack_muz(int n) /* ---------------------------------------------------------------------- */ +void DumpCustom::pack_mu(int n) +{ + double **mu = atom->mu; + int nlocal = atom->nlocal; + + for (int i = 0; i < nlocal; i++) + if (choose[i]) { + buf[n] = mu[i][3]; + n += size_one; + } +} + +/* ---------------------------------------------------------------------- */ + void DumpCustom::pack_radius(int n) { double *radius = atom->radius; @@ -2137,6 +2196,48 @@ 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; diff --git a/src/dump_custom.h b/src/dump_custom.h index 175629560d..3c9bbd2bd8 100644 --- a/src/dump_custom.h +++ b/src/dump_custom.h @@ -134,6 +134,7 @@ class DumpCustom : public Dump { void pack_mux(int); void pack_muy(int); void pack_muz(int); + void pack_mu(int); void pack_radius(int); void pack_omegax(int); void pack_omegay(int); @@ -141,6 +142,9 @@ 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); diff --git a/src/fix_adapt.cpp b/src/fix_adapt.cpp index 61f3f95973..5e79423a2a 100644 --- a/src/fix_adapt.cpp +++ b/src/fix_adapt.cpp @@ -307,25 +307,31 @@ void FixAdapt::change_settings() } else if (ad->which == ATOM) { // set radius from diameter - // also set rmass if both rmass and density are defined + // also scale rmass to new value if (ad->aparam == DIAMETER) { int mflag = 0; - if (atom->rmass_flag && atom->density_flag) mflag = 1; + if (atom->rmass_flag) mflag = 1; double PI = 4.0*atan(1.0); + double density; double *radius = atom->radius; double *rmass = atom->rmass; - double *density = atom->density; int *mask = atom->mask; int nlocal = atom->nlocal; - for (i = 0; i < nlocal; i++) - if (mask[i] & groupbit) { - radius[i] = 0.5*value; - if (mflag) rmass[i] = 4.0*PI/3.0 * - radius[i]*radius[i]*radius[i] * density[i]; - } + if (mflag == 0) { + for (i = 0; i < nlocal; i++) + if (mask[i] & groupbit) + radius[i] = 0.5*value; + } else { + for (i = 0; i < nlocal; i++) + if (mask[i] & groupbit) { + density = rmass[i] / (4.0*PI/3.0 * radius[i]*radius[i]*radius[i]); + radius[i] = 0.5*value; + rmass[i] = 4.0*PI/3.0 * radius[i]*radius[i]*radius[i] * density; + } + } } } } diff --git a/src/fix_nh_sphere.cpp b/src/fix_nh_sphere.cpp index 6dbf9566e5..d1be4184fd 100644 --- a/src/fix_nh_sphere.cpp +++ b/src/fix_nh_sphere.cpp @@ -31,12 +31,8 @@ using namespace LAMMPS_NS; FixNHSphere::FixNHSphere(LAMMPS *lmp, int narg, char **arg) : FixNH(lmp, narg, arg) { - if (!atom->omega_flag || !atom->torque_flag) - error->all("Fix nvt/nph/npt sphere requires " - "atom attributes omega, torque"); - if (!atom->radius_flag && !atom->avec->shape_type) - error->all("Fix nvt/nph/npt sphere requires " - "atom attribute diameter or shape"); + if (!atom->sphere_flag) + error->all("Fix nvt/nph/npt sphere requires atom style sphere"); } /* ---------------------------------------------------------------------- */ @@ -45,36 +41,17 @@ void FixNHSphere::init() { int i,itype; - // check that all particles are finite-size and spherical + // check that all particles are finite-size // no point particles allowed - if (atom->radius_flag) { - double *radius = atom->radius; - int *mask = atom->mask; - int nlocal = atom->nlocal; + double *radius = atom->radius; + int *mask = atom->mask; + int nlocal = atom->nlocal; - for (i = 0; i < nlocal; i++) - if (mask[i] & groupbit) { - if (radius[i] == 0.0) - error->one("Fix nvt/sphere requires extended particles"); - } - - } else { - double **shape = atom->shape; - int *type = atom->type; - int *mask = atom->mask; - int nlocal = atom->nlocal; - - for (i = 0; i < nlocal; i++) - if (mask[i] & groupbit) { - itype = type[i]; - if (shape[itype][0] == 0.0) - error->one("Fix nvt/sphere requires extended particles"); - if (shape[itype][0] != shape[itype][1] || - shape[itype][0] != shape[itype][2]) - error->one("Fix nvt/sphere requires spherical particle shapes"); - } - } + for (i = 0; i < nlocal; i++) + if (mask[i] & groupbit) + if (radius[i] == 0.0) + error->one("Fix nvt/sphere requires extended particles"); FixNH::init(); } @@ -93,9 +70,6 @@ void FixNHSphere::nve_v() double **torque = atom->torque; double *radius = atom->radius; double *rmass = atom->rmass; - double *mass = atom->mass; - double **shape = atom->shape; - int *type = atom->type; int *mask = atom->mask; int nlocal = atom->nlocal; if (igroup == atom->firstgroup) nlocal = atom->nfirst; @@ -104,57 +78,18 @@ void FixNHSphere::nve_v() double dtfrotate = dtf / INERTIA; double dtirotate; - int itype; // update omega for all particles // d_omega/dt = torque / inertia // 4 cases depending on radius vs shape and rmass vs mass - if (radius) { - if (rmass) { - for (int i = 0; i < nlocal; i++) { - if (mask[i] & groupbit) { - dtirotate = dtfrotate / (radius[i]*radius[i]*rmass[i]); - omega[i][0] += dtirotate*torque[i][0]; - omega[i][1] += dtirotate*torque[i][1]; - omega[i][2] += dtirotate*torque[i][2]; - } - } - } else { - for (int i = 0; i < nlocal; i++) { - if (mask[i] & groupbit) { - dtirotate = dtfrotate / (radius[i]*radius[i]*mass[type[i]]); - omega[i][0] += dtirotate*torque[i][0]; - omega[i][1] += dtirotate*torque[i][1]; - omega[i][2] += dtirotate*torque[i][2]; - } - } + for (int i = 0; i < nlocal; i++) + if (mask[i] & groupbit) { + dtirotate = dtfrotate / (radius[i]*radius[i]*rmass[i]); + omega[i][0] += dtirotate*torque[i][0]; + omega[i][1] += dtirotate*torque[i][1]; + omega[i][2] += dtirotate*torque[i][2]; } - - } else { - if (rmass) { - for (int i = 0; i < nlocal; i++) { - if (mask[i] & groupbit) { - itype = type[i]; - dtirotate = dtfrotate / (shape[itype][0]*shape[itype][0]*rmass[i]); - omega[i][0] += dtirotate*torque[i][0]; - omega[i][1] += dtirotate*torque[i][1]; - omega[i][2] += dtirotate*torque[i][2]; - } - } - } else { - for (int i = 0; i < nlocal; i++) { - if (mask[i] & groupbit) { - itype = type[i]; - dtirotate = dtfrotate / - (shape[itype][0]*shape[itype][0]*mass[itype]); - omega[i][0] += dtirotate*torque[i][0]; - omega[i][1] += dtirotate*torque[i][1]; - omega[i][2] += dtirotate*torque[i][2]; - } - } - } - } } /* ---------------------------------------------------------------------- diff --git a/src/fix_nve_sphere.cpp b/src/fix_nve_sphere.cpp index 7c21b63804..02968e1f13 100644 --- a/src/fix_nve_sphere.cpp +++ b/src/fix_nve_sphere.cpp @@ -53,10 +53,8 @@ FixNVESphere::FixNVESphere(LAMMPS *lmp, int narg, char **arg) : // error checks - if (!atom->omega_flag || !atom->torque_flag) - error->all("Fix nve/sphere requires atom attributes omega, torque"); - if (!atom->radius_flag && !atom->avec->shape_type) - error->all("Fix nve/sphere requires atom attribute diameter or shape"); + if (!atom->sphere_flag) + error->all("Fix nve/sphere requires atom style sphere"); if (extra == DIPOLE && !atom->mu_flag) error->all("Fix nve/sphere requires atom attribute mu"); } @@ -77,38 +75,17 @@ int FixNVESphere::setmask() void FixNVESphere::init() { - int i,itype; - - // check that all particles are finite-size and spherical + // check that all particles are finite-size // no point particles allowed - if (atom->radius_flag) { - double *radius = atom->radius; - int *mask = atom->mask; - int nlocal = atom->nlocal; + double *radius = atom->radius; + int *mask = atom->mask; + int nlocal = atom->nlocal; - for (i = 0; i < nlocal; i++) - if (mask[i] & groupbit) { - if (radius[i] == 0.0) - error->one("Fix nve/sphere requires extended particles"); - } - - } else { - double **shape = atom->shape; - int *type = atom->type; - int *mask = atom->mask; - int nlocal = atom->nlocal; - - for (i = 0; i < nlocal; i++) - if (mask[i] & groupbit) { - itype = type[i]; - if (shape[itype][0] == 0.0) - error->one("Fix nve/sphere requires extended particles"); - if (shape[itype][0] != shape[itype][1] || - shape[itype][0] != shape[itype][2]) - error->one("Fix nve/sphere requires spherical particle shapes"); - } - } + for (int i = 0; i < nlocal; i++) + if (mask[i] & groupbit) + if (radius[i] == 0.0) + error->one("Fix nve/sphere requires extended particles"); FixNVE::init(); } @@ -128,8 +105,6 @@ void FixNVESphere::initial_integrate(int vflag) double **torque = atom->torque; double *radius = atom->radius; double *rmass = atom->rmass; - double *mass = atom->mass; - double **shape = atom->shape; int *type = atom->type; int *mask = atom->mask; int nlocal = atom->nlocal; @@ -141,86 +116,21 @@ void FixNVESphere::initial_integrate(int vflag) // update v,x,omega for all particles // d_omega/dt = torque / inertia - // 4 cases depending on radius vs shape and rmass vs mass - if (radius) { - if (rmass) { - for (int i = 0; i < nlocal; i++) { - if (mask[i] & groupbit) { - dtfm = dtf / rmass[i]; - v[i][0] += dtfm * f[i][0]; - v[i][1] += dtfm * f[i][1]; - v[i][2] += dtfm * f[i][2]; - x[i][0] += dtv * v[i][0]; - x[i][1] += dtv * v[i][1]; - x[i][2] += dtv * v[i][2]; - - dtirotate = dtfrotate / (radius[i]*radius[i]*rmass[i]); - omega[i][0] += dtirotate * torque[i][0]; - omega[i][1] += dtirotate * torque[i][1]; - omega[i][2] += dtirotate * torque[i][2]; - } - } - - } else { - for (int i = 0; i < nlocal; i++) { - if (mask[i] & groupbit) { - itype = type[i]; - dtfm = dtf / mass[itype]; - v[i][0] += dtfm * f[i][0]; - v[i][1] += dtfm * f[i][1]; - v[i][2] += dtfm * f[i][2]; - x[i][0] += dtv * v[i][0]; - x[i][1] += dtv * v[i][1]; - x[i][2] += dtv * v[i][2]; - - dtirotate = dtfrotate / (radius[i]*radius[i]*mass[itype]); - omega[i][0] += dtirotate * torque[i][0]; - omega[i][1] += dtirotate * torque[i][1]; - omega[i][2] += dtirotate * torque[i][2]; - } - } - } - - } else { - if (rmass) { - for (int i = 0; i < nlocal; i++) { - if (mask[i] & groupbit) { - itype = type[i]; - dtfm = dtf / rmass[i]; - v[i][0] += dtfm * f[i][0]; - v[i][1] += dtfm * f[i][1]; - v[i][2] += dtfm * f[i][2]; - x[i][0] += dtv * v[i][0]; - x[i][1] += dtv * v[i][1]; - x[i][2] += dtv * v[i][2]; - - dtirotate = dtfrotate / (shape[itype][0]*shape[itype][0]*rmass[i]); - omega[i][0] += dtirotate * torque[i][0]; - omega[i][1] += dtirotate * torque[i][1]; - omega[i][2] += dtirotate * torque[i][2]; - } - } - - } else { - for (int i = 0; i < nlocal; i++) { - if (mask[i] & groupbit) { - itype = type[i]; - dtfm = dtf / mass[itype]; - v[i][0] += dtfm * f[i][0]; - v[i][1] += dtfm * f[i][1]; - v[i][2] += dtfm * f[i][2]; - x[i][0] += dtv * v[i][0]; - x[i][1] += dtv * v[i][1]; - x[i][2] += dtv * v[i][2]; - - dtirotate = dtfrotate / - (shape[itype][0]*shape[itype][0]*mass[itype]); - omega[i][0] += dtirotate * torque[i][0]; - omega[i][1] += dtirotate * torque[i][1]; - omega[i][2] += dtirotate * torque[i][2]; - } - } + for (int i = 0; i < nlocal; i++) { + if (mask[i] & groupbit) { + dtfm = dtf / rmass[i]; + v[i][0] += dtfm * f[i][0]; + v[i][1] += dtfm * f[i][1]; + v[i][2] += dtfm * f[i][2]; + x[i][0] += dtv * v[i][0]; + x[i][1] += dtv * v[i][1]; + x[i][2] += dtv * v[i][2]; + + dtirotate = dtfrotate / (radius[i]*radius[i]*rmass[i]); + omega[i][0] += dtirotate * torque[i][0]; + omega[i][1] += dtirotate * torque[i][1]; + omega[i][2] += dtirotate * torque[i][2]; } } @@ -230,21 +140,18 @@ void FixNVESphere::initial_integrate(int vflag) if (extra == DIPOLE) { double **mu = atom->mu; - double *dipole = atom->dipole; - for (int i = 0; i < nlocal; i++) { - if (mask[i] & groupbit) { - if (dipole[type[i]] > 0.0) { + for (int i = 0; i < nlocal; i++) + if (mask[i] & groupbit) + if (mu[i][3] > 0.0) { g[0] = mu[i][0] + dtv * (omega[i][1]*mu[i][2]-omega[i][2]*mu[i][1]); g[1] = mu[i][1] + dtv * (omega[i][2]*mu[i][0]-omega[i][0]*mu[i][2]); g[2] = mu[i][2] + dtv * (omega[i][0]*mu[i][1]-omega[i][1]*mu[i][0]); msq = g[0]*g[0] + g[1]*g[1] + g[2]*g[2]; - scale = dipole[type[i]]/sqrt(msq); + scale = mu[i][3]/sqrt(msq); mu[i][0] = g[0]*scale; mu[i][1] = g[1]*scale; mu[i][2] = g[2]*scale; } - } - } } } @@ -252,18 +159,14 @@ void FixNVESphere::initial_integrate(int vflag) void FixNVESphere::final_integrate() { - int itype; double dtfm,dtirotate; double **v = atom->v; double **f = atom->f; double **omega = atom->omega; double **torque = atom->torque; - double *mass = atom->mass; double *rmass = atom->rmass; double *radius = atom->radius; - double **shape = atom->shape; - int *type = atom->type; int *mask = atom->mask; int nlocal = atom->nlocal; if (igroup == atom->firstgroup) nlocal = atom->nfirst; @@ -274,74 +177,17 @@ void FixNVESphere::final_integrate() // update v,omega for all particles // d_omega/dt = torque / inertia - // 4 cases depending on radius vs shape and rmass vs mass - if (radius) { - if (rmass) { - for (int i = 0; i < nlocal; i++) { - if (mask[i] & groupbit) { - dtfm = dtf / rmass[i]; - v[i][0] += dtfm * f[i][0]; - v[i][1] += dtfm * f[i][1]; - v[i][2] += dtfm * f[i][2]; - - dtirotate = dtfrotate / (radius[i]*radius[i]*rmass[i]); - omega[i][0] += dtirotate * torque[i][0]; - omega[i][1] += dtirotate * torque[i][1]; - omega[i][2] += dtirotate * torque[i][2]; - } - } - - } else { - for (int i = 0; i < nlocal; i++) { - if (mask[i] & groupbit) { - itype = type[i]; - dtfm = dtf / mass[itype]; - v[i][0] += dtfm * f[i][0]; - v[i][1] += dtfm * f[i][1]; - v[i][2] += dtfm * f[i][2]; - - dtirotate = dtfrotate / (radius[i]*radius[i]*mass[itype]); - omega[i][0] += dtirotate * torque[i][0]; - omega[i][1] += dtirotate * torque[i][1]; - omega[i][2] += dtirotate * torque[i][2]; - } - } + for (int i = 0; i < nlocal; i++) + if (mask[i] & groupbit) { + dtfm = dtf / rmass[i]; + v[i][0] += dtfm * f[i][0]; + v[i][1] += dtfm * f[i][1]; + v[i][2] += dtfm * f[i][2]; + + dtirotate = dtfrotate / (radius[i]*radius[i]*rmass[i]); + omega[i][0] += dtirotate * torque[i][0]; + omega[i][1] += dtirotate * torque[i][1]; + omega[i][2] += dtirotate * torque[i][2]; } - - } else { - if (rmass) { - for (int i = 0; i < nlocal; i++) { - if (mask[i] & groupbit) { - itype = type[i]; - dtfm = dtf / rmass[i]; - v[i][0] += dtfm * f[i][0]; - v[i][1] += dtfm * f[i][1]; - v[i][2] += dtfm * f[i][2]; - - dtirotate = dtfrotate / (shape[itype][0]*shape[itype][0]*rmass[i]); - omega[i][0] += dtirotate * torque[i][0]; - omega[i][1] += dtirotate * torque[i][1]; - omega[i][2] += dtirotate * torque[i][2]; - } - } - - } else { - for (int i = 0; i < nlocal; i++) { - if (mask[i] & groupbit) { - itype = type[i]; - dtfm = dtf / mass[itype]; - v[i][0] += dtfm * f[i][0]; - v[i][1] += dtfm * f[i][1]; - v[i][2] += dtfm * f[i][2]; - - dtirotate = dtfrotate / - (shape[itype][0]*shape[itype][0]*mass[itype]); - omega[i][0] += dtirotate * torque[i][0]; - omega[i][1] += dtirotate * torque[i][1]; - omega[i][2] += dtirotate * torque[i][2]; - } - } - } - } } diff --git a/src/fix_rigid.cpp b/src/fix_rigid.cpp index 4cfc52d271..bfea6f4a49 100644 --- a/src/fix_rigid.cpp +++ b/src/fix_rigid.cpp @@ -352,8 +352,8 @@ FixRigid::FixRigid(LAMMPS *lmp, int narg, char **arg) : // bitmasks for properties of extended particles - INERTIA_SPHERE_RADIUS = 1; - INERTIA_SPHERE_SHAPE = 2; + INERTIA_POINT = 1; + INERTIA_SPHERE = 2; INERTIA_ELLIPSOID = 4; ORIENT_DIPOLE = 8; ORIENT_QUAT = 16; @@ -462,30 +462,32 @@ void FixRigid::init() step_respa = ((Respa *) update->integrate)->step; // extended = 1 if any particle in a rigid body is finite size + // or has a dipole moment extended = dorientflag = qorientflag = 0; + double **shape = atom->shape; + double **mu = atom->mu; double *radius = atom->radius; double *rmass = atom->rmass; double *mass = atom->mass; - double **shape = atom->shape; - double *dipole = atom->dipole; int *type = atom->type; int nlocal = atom->nlocal; - if (atom->radius_flag || atom->avec->shape_type) { + if (atom->radius_flag || atom->shape_flag) { int flag = 0; for (i = 0; i < nlocal; i++) { if (body[i] < 0) continue; if (radius && radius[i] > 0.0) flag = 1; - else if (shape && shape[type[i]][0] > 0.0) flag = 1; + if (shape && shape[i][0] > 0.0) flag = 1; + if (mu && mu[i][3] > 0.0) flag = 1; } MPI_Allreduce(&flag,&extended,1,MPI_INT,MPI_MAX,world); } // grow extended arrays and set extended flags for each particle - // vorientflag = 1 if any particles store dipole orientation + // dorientflag = 1 if any particles store dipole orientation // qorientflag = 1 if any particles store quat orientation if (extended) { @@ -497,39 +499,23 @@ void FixRigid::init() eflags[i] = 0; if (body[i] < 0) continue; - // set INERTIA if radius or shape > 0.0 + // set INERTIA to POINT or SPHERE or ELLIPSOID - if (radius) { - if (radius[i] > 0.0) eflags[i] |= INERTIA_SPHERE_RADIUS; - } else if (shape) { - if (shape[type[i]][0] > 0.0) { - if (shape[type[i]][0] == shape[type[i]][1] && - shape[type[i]][0] == shape[type[i]][2]) - eflags[i] |= INERTIA_SPHERE_SHAPE; - else eflags[i] |= INERTIA_ELLIPSOID; - } - } + if (radius && radius[i] > 0.0) { + eflags[i] |= INERTIA_SPHERE; + eflags[i] |= OMEGA; + eflags[i] |= TORQUE; + } else if (shape && shape[i][0] > 0.0) { + eflags[i] |= INERTIA_ELLIPSOID; + eflags[i] |= ORIENT_QUAT; + eflags[i] |= ANGMOM; + eflags[i] |= TORQUE; + } else eflags[i] |= INERTIA_POINT; - // other flags only set if particle is finite size - // set DIPOLE if atom->mu and atom->dipole exist and dipole[itype] > 0.0 - // set QUAT if atom->quat exists (could be ellipsoid or sphere) - // set TORQUE if atom->torque exists - // set exactly one of OMEGA and ANGMOM so particle contributes once - // set OMEGA if either radius or rmass exists - // set ANGMOM if shape and mass exist - // set OMEGA if atom->angmom doesn't exist + // set DIPOLE if atom->mu and mu[3] > 0.0 - if (eflags[i] == 0) continue; - - if (atom->mu_flag && dipole && dipole[type[i]] > 0.0) + if (atom->mu_flag && mu[i][3] > 0.0) eflags[i] |= ORIENT_DIPOLE; - if (atom->quat_flag) eflags[i] |= ORIENT_QUAT; - if (atom->torque_flag) eflags[i] |= TORQUE; - if ((radius || rmass) && atom->omega_flag) eflags[i] |= OMEGA; - else if (shape && mass && atom->angmom_flag) eflags[i] |= ANGMOM; - else if (atom->omega_flag) eflags[i] |= OMEGA; - else error->one("Could not set finite-size particle attribute " - "in fix rigid"); } } @@ -638,7 +624,7 @@ void FixRigid::init() sum[ibody][5] -= massone * dx*dz; } - // extended particles contribute extra terms to moments of inertia + // extended particles may contribute extra terms to moments of inertia if (extended) { double ex[3],ey[3],ez[3],idiag[3]; @@ -655,26 +641,20 @@ void FixRigid::init() if (rmass) massone = rmass[i]; else massone = mass[itype]; - if (eflags[i] & INERTIA_SPHERE_RADIUS) { + if (eflags[i] & INERTIA_SPHERE) { sum[ibody][0] += 0.4 * massone * radius[i]*radius[i]; sum[ibody][1] += 0.4 * massone * radius[i]*radius[i]; sum[ibody][2] += 0.4 * massone * radius[i]*radius[i]; } - if (eflags[i] & INERTIA_SPHERE_SHAPE) { - rad = shape[type[i]][0]; - sum[ibody][0] += 0.4 * massone * rad*rad; - sum[ibody][1] += 0.4 * massone * rad*rad; - sum[ibody][2] += 0.4 * massone * rad*rad; - } if (eflags[i] & INERTIA_ELLIPSOID) { MathExtra::quat_to_mat(atom->quat[i],p); MathExtra::quat_to_mat_trans(atom->quat[i],ptrans); idiag[0] = 0.2*massone * - (shape[itype][1]*shape[itype][1]+shape[itype][2]*shape[itype][2]); + (shape[i][1]*shape[i][1]+shape[itype][2]*shape[i][2]); idiag[1] = 0.2*massone * - (shape[itype][0]*shape[itype][0]+shape[itype][2]*shape[itype][2]); + (shape[i][0]*shape[i][0]+shape[i][2]*shape[i][2]); idiag[2] = 0.2*massone * - (shape[itype][0]*shape[itype][0]+shape[itype][1]*shape[itype][1]); + (shape[i][0]*shape[i][0]+shape[i][1]*shape[i][1]); MathExtra::diag_times3(idiag,ptrans,itemp); MathExtra::times3(p,itemp,ispace); sum[ibody][0] += ispace[0][0]; @@ -800,7 +780,6 @@ void FixRigid::init() if (extended) { double **mu = atom->mu; - double *dipole = atom->dipole; int *type = atom->type; if (eflags[i] & ORIENT_DIPOLE) { @@ -810,7 +789,7 @@ void FixRigid::init() mu[i][1]*ey_space[ibody][1] + mu[i][2]*ey_space[ibody][2]; dorient[i][2] = mu[i][0]*ez_space[ibody][0] + mu[i][1]*ez_space[ibody][1] + mu[i][2]*ez_space[ibody][2]; - MathExtra::snormalize3(dipole[type[i]],dorient[i],dorient[i]); + MathExtra::snormalize3(mu[i][3],dorient[i],dorient[i]); } else if (dorientflag) dorient[i][0] = dorient[i][1] = dorient[i][2] = 0.0; @@ -827,7 +806,7 @@ void FixRigid::init() // recompute moments of inertia around new axes // 3 diagonal moments should equal principal moments // 3 off-diagonal moments should be 0.0 - // extended particles contribute extra terms to moments of inertia + // extended particles may contribute extra terms to moments of inertia for (ibody = 0; ibody < nbody; ibody++) for (i = 0; i < 6; i++) sum[ibody][i] = 0.0; @@ -864,26 +843,20 @@ void FixRigid::init() if (rmass) massone = rmass[i]; else massone = mass[itype]; - if (eflags[i] & INERTIA_SPHERE_RADIUS) { + if (eflags[i] & INERTIA_SPHERE) { sum[ibody][0] += 0.4 * massone * radius[i]*radius[i]; sum[ibody][1] += 0.4 * massone * radius[i]*radius[i]; sum[ibody][2] += 0.4 * massone * radius[i]*radius[i]; } - if (eflags[i] & INERTIA_SPHERE_SHAPE) { - rad = shape[type[i]][0]; - sum[ibody][0] += 0.4 * massone * rad*rad; - sum[ibody][1] += 0.4 * massone * rad*rad; - sum[ibody][2] += 0.4 * massone * rad*rad; - } if (eflags[i] & INERTIA_ELLIPSOID) { MathExtra::quat_to_mat(qorient[i],p); MathExtra::quat_to_mat_trans(qorient[i],ptrans); idiag[0] = 0.2*massone * - (shape[itype][1]*shape[itype][1]+shape[itype][2]*shape[itype][2]); + (shape[i][1]*shape[i][1]+shape[i][2]*shape[i][2]); idiag[1] = 0.2*massone * - (shape[itype][0]*shape[itype][0]+shape[itype][2]*shape[itype][2]); + (shape[i][0]*shape[i][0]+shape[i][2]*shape[i][2]); idiag[2] = 0.2*massone * - (shape[itype][0]*shape[itype][0]+shape[itype][1]*shape[itype][1]); + (shape[i][0]*shape[i][0]+shape[i][1]*shape[i][1]); MathExtra::diag_times3(idiag,ptrans,itemp); MathExtra::times3(p,itemp,ispace); sum[ibody][0] += ispace[0][0]; @@ -1032,7 +1005,6 @@ void FixRigid::setup(int vflag) double **angmom_one = atom->angmom; double **torque_one = atom->torque; double *radius = atom->radius; - double **shape = atom->shape; for (i = 0; i < nlocal; i++) { if (body[i] < 0) continue; @@ -1041,8 +1013,7 @@ void FixRigid::setup(int vflag) if (eflags[i] & OMEGA) { if (rmass) massone = rmass[i]; else massone = mass[type[i]]; - if (radius) radone = radius[i]; - else radone = shape[type[i]][0]; + radone = radius[i]; sum[ibody][0] += 0.4 * massone * radone*radone * omega_one[i][0]; sum[ibody][1] += 0.4 * massone * radone*radone * omega_one[i][1]; sum[ibody][2] += 0.4 * massone * radone*radone * omega_one[i][2]; @@ -1965,8 +1936,8 @@ void FixRigid::set_xv() if (extended) { double **omega_one = atom->omega; double **angmom_one = atom->angmom; - double *dipole = atom->dipole; double **shape = atom->shape; + double **mu = atom->mu; for (int i = 0; i < nlocal; i++) { if (body[i] < 0) continue; @@ -1974,8 +1945,8 @@ void FixRigid::set_xv() if (eflags[i] & ORIENT_DIPOLE) { MathExtra::quat_to_mat(quat[ibody],p); - MathExtra::times_column3(p,dorient[i],atom->mu[i]); - MathExtra::snormalize3(dipole[type[i]],atom->mu[i],atom->mu[i]); + MathExtra::times_column3(p,dorient[i],mu[i]); + MathExtra::snormalize3(mu[i][3],mu[i],mu[i]); } if (eflags[i] & ORIENT_QUAT) { quatquat(quat[ibody],qorient[i],atom->quat[i]); @@ -1987,13 +1958,12 @@ void FixRigid::set_xv() omega_one[i][2] = omega[ibody][2]; } if (eflags[i] & ANGMOM) { - itype = type[i]; - ione[0] = 0.2*mass[itype] * - (shape[itype][1]*shape[itype][1] + shape[itype][2]*shape[itype][2]); - ione[1] = 0.2*mass[itype] * - (shape[itype][0]*shape[itype][0] + shape[itype][2]*shape[itype][2]); - ione[2] = 0.2*mass[itype] * - (shape[itype][0]*shape[itype][0] + shape[itype][1]*shape[itype][1]); + 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); angmom_from_omega(omega[ibody],exone,eyone,ezone,ione,angmom_one[i]); } @@ -2117,13 +2087,12 @@ void FixRigid::set_v() omega_one[i][2] = omega[ibody][2]; } if (eflags[i] & ANGMOM) { - itype = type[i]; - ione[0] = 0.2*mass[itype] * - (shape[itype][1]*shape[itype][1] + shape[itype][2]*shape[itype][2]); - ione[1] = 0.2*mass[itype] * - (shape[itype][0]*shape[itype][0] + shape[itype][2]*shape[itype][2]); - ione[2] = 0.2*mass[itype] * - (shape[itype][0]*shape[itype][0] + shape[itype][1]*shape[itype][1]); + 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); 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 efe77587b4..962a83d7bb 100644 --- a/src/fix_rigid.h +++ b/src/fix_rigid.h @@ -96,7 +96,7 @@ class FixRigid : public Fix { int p_chain; // bitmasks for eflags - int INERTIA_SPHERE_RADIUS,INERTIA_SPHERE_SHAPE,INERTIA_ELLIPSOID; + int INERTIA_POINT,INERTIA_SPHERE,INERTIA_ELLIPSOID; int ORIENT_DIPOLE,ORIENT_QUAT; int OMEGA,ANGMOM,TORQUE; diff --git a/src/fix_wall_region.cpp b/src/fix_wall_region.cpp index 9c7cb70648..a337d58b12 100644 --- a/src/fix_wall_region.cpp +++ b/src/fix_wall_region.cpp @@ -98,32 +98,20 @@ void FixWallRegion::init() error->all("Region ID for fix wall/region does not exist"); // error checks for style COLLOID - // insure all particle shapes are spherical - // can be polydisperse // insure all particles in group are extended particles if (style == COLLOID) { - if (!atom->avec->shape_type) - error->all("Fix wall/region colloid requires atom attribute shape"); - if (atom->radius_flag) - error->all("Fix wall/region colloid cannot be used with " - "atom attribute diameter"); + if (!atom->sphere_flag) + error->all("Fix wall/region colloid requires atom style sphere"); - for (int i = 1; i <= atom->ntypes; i++) - if ((atom->shape[i][0] != atom->shape[i][1]) || - (atom->shape[i][0] != atom->shape[i][2]) || - (atom->shape[i][1] != atom->shape[i][2])) - error->all("Fix wall/region colloid requires spherical particles"); - - double **shape = atom->shape; - int *type = atom->type; + double *radius = atom->radius; int *mask = atom->mask; int nlocal = atom->nlocal; int flag = 0; for (int i = 0; i < nlocal; i++) if (mask[i] & groupbit) - if (atom->shape[type[i]][0] == 0.0) flag = 1; + if (radius[i] == 0.0) flag = 1; int flagall; MPI_Allreduce(&flag,&flagall,1,MPI_INT,MPI_SUM,world); @@ -197,7 +185,7 @@ void FixWallRegion::post_force(int vflag) double **x = atom->x; double **f = atom->f; - double **shape = atom->shape; + double *radius = atom->radius; int *type = atom->type; int *mask = atom->mask; int nlocal = atom->nlocal; @@ -207,7 +195,7 @@ void FixWallRegion::post_force(int vflag) // region->match() insures particle is in region or on surface, else error // if returned contact dist r = 0, is on surface, also an error - // in COLLOID case, r <= shape radius is an error + // in COLLOID case, r <= radius is an error for (i = 0; i < nlocal; i++) if (mask[i] & groupbit) { @@ -215,7 +203,7 @@ void FixWallRegion::post_force(int vflag) onflag = 1; continue; } - if (style == COLLOID) tooclose = shape[type[i]][0]; + if (style == COLLOID) tooclose = radius[i]; else tooclose = 0.0; n = region->surface(x[i][0],x[i][1],x[i][2],cutoff); @@ -228,8 +216,7 @@ void FixWallRegion::post_force(int vflag) if (style == LJ93) lj93(region->contact[m].r); else if (style == LJ126) lj126(region->contact[m].r); - else if (style == COLLOID) - colloid(region->contact[m].r,shape[type[i]][0]); + else if (style == COLLOID) colloid(region->contact[m].r,radius[i]); else harmonic(region->contact[m].r); ewall[0] += eng; diff --git a/src/input.cpp b/src/input.cpp index c0bbb6ca66..0861c96cc2 100644 --- a/src/input.cpp +++ b/src/input.cpp @@ -431,7 +431,6 @@ int Input::execute_command() else if (!strcmp(command,"dihedral_coeff")) dihedral_coeff(); else if (!strcmp(command,"dihedral_style")) dihedral_style(); else if (!strcmp(command,"dimension")) dimension(); - else if (!strcmp(command,"dipole")) dipole(); else if (!strcmp(command,"dump")) dump(); else if (!strcmp(command,"dump_modify")) dump_modify(); else if (!strcmp(command,"fix")) fix(); @@ -457,7 +456,6 @@ int Input::execute_command() else if (!strcmp(command,"reset_timestep")) reset_timestep(); else if (!strcmp(command,"restart")) restart(); else if (!strcmp(command,"run_style")) run_style(); - else if (!strcmp(command,"shape")) shape(); else if (!strcmp(command,"special_bonds")) special_bonds(); else if (!strcmp(command,"thermo")) thermo(); else if (!strcmp(command,"thermo_modify")) thermo_modify(); @@ -948,16 +946,6 @@ void Input::dimension() /* ---------------------------------------------------------------------- */ -void Input::dipole() -{ - if (narg != 2) error->all("Illegal dipole command"); - if (domain->box_exist == 0) - error->all("Dipole command before simulation box is defined"); - atom->set_dipole(narg,arg); -} - -/* ---------------------------------------------------------------------- */ - void Input::dump() { output->add_dump(narg,arg); @@ -1208,16 +1196,6 @@ void Input::run_style() /* ---------------------------------------------------------------------- */ -void Input::shape() -{ - if (narg != 4) error->all("Illegal shape command"); - if (domain->box_exist == 0) - error->all("Shape command before simulation box is defined"); - atom->set_shape(narg,arg); -} - -/* ---------------------------------------------------------------------- */ - void Input::special_bonds() { // store 1-3,1-4 and dihedral/extra flag values before change diff --git a/src/input.h b/src/input.h index 7d2946d3e9..cbd27881ef 100644 --- a/src/input.h +++ b/src/input.h @@ -75,7 +75,6 @@ class Input : protected Pointers { void dihedral_coeff(); void dihedral_style(); void dimension(); - void dipole(); void dump(); void dump_modify(); void fix(); @@ -101,7 +100,6 @@ class Input : protected Pointers { void reset_timestep(); void restart(); void run_style(); - void shape(); void special_bonds(); void thermo(); void thermo_modify(); diff --git a/src/lammps.cpp b/src/lammps.cpp index 4a52d9fa0e..15b32bac2f 100644 --- a/src/lammps.cpp +++ b/src/lammps.cpp @@ -338,7 +338,7 @@ void LAMMPS::init() // 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 + comm->init(); // comm must come after force, modify, neighbor, atom output->init(); // output must come after domain, force, modify timer->init(); } diff --git a/src/pair_yukawa.cpp b/src/pair_yukawa.cpp index 430a9604bf..65a38096fb 100644 --- a/src/pair_yukawa.cpp +++ b/src/pair_yukawa.cpp @@ -38,6 +38,7 @@ PairYukawa::~PairYukawa() memory->destroy(setflag); memory->destroy(cutsq); + memory->destroy(rad); memory->destroy(cut); memory->destroy(a); memory->destroy(offset); @@ -139,6 +140,7 @@ void PairYukawa::allocate() memory->create(cutsq,n+1,n+1,"pair:cutsq"); + memory->create(rad,n+1,"pair:rad"); memory->create(cut,n+1,n+1,"pair:cut"); memory->create(a,n+1,n+1,"pair:a"); memory->create(offset,n+1,n+1,"pair:offset"); diff --git a/src/pair_yukawa.h b/src/pair_yukawa.h index 9d1671b734..a7dfab9474 100644 --- a/src/pair_yukawa.h +++ b/src/pair_yukawa.h @@ -41,6 +41,7 @@ class PairYukawa : public Pair { protected: double cut_global; double kappa; + double *rad; double **cut,**a,**offset; void allocate(); diff --git a/src/read_data.cpp b/src/read_data.cpp index 85d040a87a..301498cf62 100644 --- a/src/read_data.cpp +++ b/src/read_data.cpp @@ -163,10 +163,6 @@ void ReadData::command(int narg, char **arg) } else if (strcmp(keyword,"Masses") == 0) { mass(); - } else if (strcmp(keyword,"Shapes") == 0) { - shape(); - } else if (strcmp(keyword,"Dipoles") == 0) { - dipole(); } else if (strcmp(keyword,"Pair Coeffs") == 0) { if (force->pair == NULL) error->all("Must define pair_style before Pair Coeffs"); @@ -292,8 +288,7 @@ void ReadData::header(int flag) char *ptr; char *section_keywords[NSECTIONS] = - {"Atoms","Velocities","Bonds","Angles","Dihedrals","Impropers", - "Masses","Shapes","Dipoles", + {"Atoms","Velocities","Bonds","Angles","Dihedrals","Impropers","Masses", "Pair Coeffs","Bond Coeffs","Angle Coeffs", "Dihedral Coeffs","Improper Coeffs", "BondBond Coeffs","BondAngle Coeffs","MiddleBondTorsion Coeffs", @@ -761,64 +756,6 @@ void ReadData::mass() /* ---------------------------------------------------------------------- */ -void ReadData::shape() -{ - int i,m; - char *buf = new char[atom->ntypes*MAXLINE]; - char *original = buf; - - if (me == 0) { - char *eof; - m = 0; - for (i = 0; i < atom->ntypes; i++) { - eof = fgets(&buf[m],MAXLINE,fp); - if (eof == NULL) error->one("Unexpected end of data file"); - m += strlen(&buf[m]); - buf[m-1] = '\0'; - } - } - - MPI_Bcast(&m,1,MPI_INT,0,world); - MPI_Bcast(buf,m,MPI_CHAR,0,world); - - for (i = 0; i < atom->ntypes; i++) { - atom->set_shape(buf); - buf += strlen(buf) + 1; - } - delete [] original; -} - -/* ---------------------------------------------------------------------- */ - -void ReadData::dipole() -{ - int i,m; - char *buf = new char[atom->ntypes*MAXLINE]; - char *original = buf; - - if (me == 0) { - char *eof; - m = 0; - for (i = 0; i < atom->ntypes; i++) { - eof = fgets(&buf[m],MAXLINE,fp); - if (eof == NULL) error->one("Unexpected end of data file"); - m += strlen(&buf[m]); - buf[m-1] = '\0'; - } - } - - MPI_Bcast(&m,1,MPI_INT,0,world); - MPI_Bcast(buf,m,MPI_CHAR,0,world); - - for (i = 0; i < atom->ntypes; i++) { - atom->set_dipole(buf); - buf += strlen(buf) + 1; - } - delete [] original; -} - -/* ---------------------------------------------------------------------- */ - void ReadData::paircoeffs() { int i,m; @@ -1007,7 +944,6 @@ void ReadData::scan(int &bond_per_atom, int &angle_per_atom, while (strlen(keyword)) { if (strcmp(keyword,"Masses") == 0) skip_lines(atom->ntypes); - else if (strcmp(keyword,"Dipoles") == 0) skip_lines(atom->ntypes); else if (strcmp(keyword,"Atoms") == 0) skip_lines(natoms); else if (strcmp(keyword,"Velocities") == 0) skip_lines(natoms); diff --git a/src/read_data.h b/src/read_data.h index 6f0dfacd71..f5c4357d4f 100644 --- a/src/read_data.h +++ b/src/read_data.h @@ -54,7 +54,6 @@ class ReadData : protected Pointers { void impropers(); void mass(); - void shape(); void dipole(); void paircoeffs(); diff --git a/src/read_restart.cpp b/src/read_restart.cpp index f674e18b5b..0d761dd7c8 100644 --- a/src/read_restart.cpp +++ b/src/read_restart.cpp @@ -56,7 +56,7 @@ enum{VERSION,SMALLINT,TAGINT,BIGINT, SPECIAL_LJ_1,SPECIAL_LJ_2,SPECIAL_LJ_3, SPECIAL_COUL_1,SPECIAL_COUL_2,SPECIAL_COUL_3, XY,XZ,YZ}; -enum{MASS,SHAPE,DIPOLE}; +enum{MASS}; enum{PAIR,BOND,ANGLE,DIHEDRAL,IMPROPER}; #define LB_FACTOR 1.1 @@ -700,21 +700,6 @@ void ReadRestart::type_arrays() atom->set_mass(mass); delete [] mass; - } else if (flag == SHAPE) { - double **shape; - memory->create(shape,atom->ntypes+1,3,"restart:shape"); - if (me == 0) fread(&shape[1][0],sizeof(double),atom->ntypes*3,fp); - MPI_Bcast(&shape[1][0],atom->ntypes*3,MPI_DOUBLE,0,world); - atom->set_shape(shape); - memory->destroy(shape); - - } else if (flag == DIPOLE) { - double *dipole = new double[atom->ntypes+1]; - if (me == 0) fread(&dipole[1],sizeof(double),atom->ntypes,fp); - MPI_Bcast(&dipole[1],atom->ntypes,MPI_DOUBLE,0,world); - atom->set_dipole(dipole); - delete [] dipole; - } else error->all("Invalid flag in type arrays section of restart file"); flag = read_int(); diff --git a/src/replicate.cpp b/src/replicate.cpp index 8bc146794b..b9767e424d 100644 --- a/src/replicate.cpp +++ b/src/replicate.cpp @@ -208,25 +208,6 @@ void Replicate::command(int narg, char **arg) } } - if (atom->shape) { - for (int itype = 1; itype <= atom->ntypes; itype++) { - atom->shape_setflag[itype] = old->shape_setflag[itype]; - if (atom->shape_setflag[itype]) { - atom->shape[itype][0] = old->shape[itype][0]; - atom->shape[itype][1] = old->shape[itype][1]; - atom->shape[itype][2] = old->shape[itype][2]; - } - } - } - - if (atom->dipole) { - for (int itype = 1; itype <= atom->ntypes; itype++) { - atom->dipole_setflag[itype] = old->dipole_setflag[itype]; - if (atom->dipole_setflag[itype]) - atom->dipole[itype] = old->dipole[itype]; - } - } - // set bounds for my proc // if periodic and I am lo/hi proc, adjust bounds by EPSILON // insures all replicated atoms will be owned even with round-off diff --git a/src/set.cpp b/src/set.cpp index a75a6a40c7..e1b1a71036 100644 --- a/src/set.cpp +++ b/src/set.cpp @@ -31,12 +31,12 @@ using namespace LAMMPS_NS; -enum{ATOM,GROUP,REGION}; -enum{TYPE,TYPE_FRACTION,MOLECULE, - X,Y,Z,CHARGE, - DIPOLE,DIPOLE_RANDOM,QUAT,QUAT_RANDOM, - DIAMETER,DENSITY,VOLUME,IMAGE, - BOND,ANGLE,DIHEDRAL,IMPROPER}; +enum{ATOM_SELECT,MOL_SELECT,TYPE_SELECT,GROUP_SELECT,REGION_SELECT}; +enum{TYPE,TYPE_FRACTION,MOLECULE,X,Y,Z,CHARGE,MASS,SHAPE, + DIPOLE,DIPOLE_RANDOM,QUAT,QUAT_RANDOM, + DIAMETER,DENSITY,VOLUME,IMAGE,BOND,ANGLE,DIHEDRAL,IMPROPER}; + +#define BIG INT_MAX /* ---------------------------------------------------------------------- */ @@ -55,11 +55,13 @@ void Set::command(int narg, char **arg) error->all("Set command with no atoms existing"); if (narg < 3) error->all("Illegal set command"); - // style and ID + // style and ID info - if (strcmp(arg[0],"atom") == 0) style = ATOM; - else if (strcmp(arg[0],"group") == 0) style = GROUP; - else if (strcmp(arg[0],"region") == 0) style = REGION; + if (strcmp(arg[0],"atom") == 0) style = ATOM_SELECT; + else if (strcmp(arg[0],"mol") == 0) style = MOL_SELECT; + else if (strcmp(arg[0],"type") == 0) style = TYPE_SELECT; + else if (strcmp(arg[0],"group") == 0) style = GROUP_SELECT; + else if (strcmp(arg[0],"region") == 0) style = REGION_SELECT; else error->all("Illegal set command"); int n = strlen(arg[1]) + 1; @@ -91,6 +93,8 @@ void Set::command(int narg, char **arg) newtype = atoi(arg[iarg+1]); fraction = atof(arg[iarg+2]); ivalue = atoi(arg[iarg+3]); + if (newtype <= 0 || newtype > atom->ntypes) + error->all("Invalid value in set command"); if (fraction < 0.0 || fraction > 1.0) error->all("Invalid value in set command"); if (ivalue <= 0) error->all("Invalid random number seed in set command"); @@ -125,23 +129,48 @@ void Set::command(int narg, char **arg) error->all("Cannot set this attribute for this atom style"); set(CHARGE); iarg += 2; + } else if (strcmp(arg[iarg],"mass") == 0) { + if (iarg+2 > narg) error->all("Illegal set command"); + dvalue = atof(arg[iarg+1]); + if (!atom->rmass_flag) + error->all("Cannot set this attribute for this atom style"); + if (dvalue <= 0.0) error->all("Invalid mass in set command"); + set(MASS); + iarg += 2; + } else if (strcmp(arg[iarg],"shape") == 0) { + if (iarg+4 > narg) error->all("Illegal set command"); + xvalue = atof(arg[iarg+1]); + yvalue = atof(arg[iarg+2]); + zvalue = atof(arg[iarg+3]); + if (!atom->shape_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"); + if (xvalue > 0.0 || yvalue > 0.0 || zvalue > 0.0) { + if (xvalue == 0.0 || yvalue == 0.0 || zvalue == 0.0) + error->one("Invalid shape in set command"); + } + set(SHAPE); + iarg += 4; } else if (strcmp(arg[iarg],"dipole") == 0) { if (iarg+4 > narg) error->all("Illegal set command"); xvalue = atof(arg[iarg+1]); yvalue = atof(arg[iarg+2]); zvalue = atof(arg[iarg+3]); - if (!atom->mu_flag || atom->dipole == NULL) + if (!atom->mu_flag) error->all("Cannot set this attribute for this atom style"); set(DIPOLE); iarg += 4; } else if (strcmp(arg[iarg],"dipole/random") == 0) { - if (iarg+2 > narg) error->all("Illegal set command"); + if (iarg+3 > narg) error->all("Illegal set command"); ivalue = atoi(arg[iarg+1]); - if (!atom->mu_flag || atom->shape == NULL) + dvalue = atof(arg[iarg+2]); + if (!atom->mu_flag) error->all("Cannot set this attribute for this atom style"); if (ivalue <= 0) error->all("Invalid random number seed in set command"); + if (dvalue <= 0.0) error->all("Invalid dipole length in set command"); setrandom(DIPOLE_RANDOM); - iarg += 2; + iarg += 3; } else if (strcmp(arg[iarg],"quat") == 0) { if (iarg+5 > narg) error->all("Illegal set command"); xvalue = atof(arg[iarg+1]); @@ -165,12 +194,13 @@ void Set::command(int narg, char **arg) dvalue = atof(arg[iarg+1]); if (!atom->radius_flag) error->all("Cannot set this attribute for this atom style"); + if (dvalue < 0.0) error->all("Invalid diameter in set command"); set(DIAMETER); iarg += 2; } else if (strcmp(arg[iarg],"density") == 0) { if (iarg+2 > narg) error->all("Illegal set command"); dvalue = atof(arg[iarg+1]); - if (!atom->density_flag) + if (!atom->rmass_flag) error->all("Cannot set this attribute for this atom style"); set(DENSITY); iarg += 2; @@ -261,7 +291,7 @@ void Set::command(int narg, char **arg) } /* ---------------------------------------------------------------------- - select atoms according to ATOM, GROUP, REGION style + select atoms according to ATOM, MOLECULE, TYPE, GROUP, REGION style n = nlocal or nlocal+nghost depending on keyword ------------------------------------------------------------------------- */ @@ -269,18 +299,38 @@ void Set::selection(int n) { delete [] select; select = new int[n]; + int nlo,nhi; - if (style == ATOM) { + if (style == ATOM_SELECT) { if (atom->tag_enable == 0) error->all("Cannot use set atom with no atom IDs defined"); - int idatom = atoi(id); + force->bounds(id,BIG,nlo,nhi); int *tag = atom->tag; for (int i = 0; i < n; i++) - if (idatom == tag[i]) select[i] = 1; + if (tag[i] >= nlo && tag[i] <= nhi) select[i] = 1; else select[i] = 0; - } else if (style == GROUP) { + } else if (style == MOL_SELECT) { + if (atom->molecule_flag == 0) + error->all("Cannot use set mol with no molecule IDs defined"); + if (strcmp(id,"0") == 0) nlo = nhi = 0; + else force->bounds(id,BIG,nlo,nhi); + + int *molecule = atom->molecule; + for (int i = 0; i < n; i++) + if (molecule[i] >= nlo && molecule[i] <= nhi) select[i] = 1; + else select[i] = 0; + + } else if (style == TYPE_SELECT) { + force->bounds(id,atom->ntypes,nlo,nhi); + + int *type = atom->type; + for (int i = 0; i < n; i++) + if (type[i] >= nlo && type[i] <= nhi) select[i] = 1; + else select[i] = 0; + + } else if (style == GROUP_SELECT) { int igroup = group->find(id); if (igroup == -1) error->all("Could not find set group ID"); int groupbit = group->bitmask[igroup]; @@ -290,7 +340,7 @@ void Set::selection(int n) if (mask[i] & groupbit) select[i] = 1; else select[i] = 0; - } else { + } else if (style == REGION_SELECT) { int iregion = domain->find_region(id); if (iregion == -1) error->all("Set region ID does not exist"); @@ -308,83 +358,81 @@ void Set::selection(int n) void Set::set(int keyword) { - if (keyword == DIPOLE) atom->check_dipole(); - selection(atom->nlocal); int nlocal = atom->nlocal; - for (int i = 0; i < nlocal; i++) - if (select[i]) { - if (keyword == TYPE) atom->type[i] = ivalue; - else if (keyword == MOLECULE) atom->molecule[i] = ivalue; - else if (keyword == X) atom->x[i][0] = dvalue; - else if (keyword == Y) atom->x[i][1] = dvalue; - else if (keyword == Z) atom->x[i][2] = dvalue; - else if (keyword == CHARGE) atom->q[i] = dvalue; + for (int i = 0; i < nlocal; i++) { + if (!select[i]) continue; - // set radius from diameter - // set rmass if both rmass and density are defined + if (keyword == TYPE) atom->type[i] = ivalue; + else if (keyword == MOLECULE) atom->molecule[i] = ivalue; + else if (keyword == X) atom->x[i][0] = dvalue; + else if (keyword == Y) atom->x[i][1] = dvalue; + else if (keyword == Z) atom->x[i][2] = dvalue; + else if (keyword == CHARGE) atom->q[i] = dvalue; + else if (keyword == MASS) atom->rmass[i] = dvalue; + else if (keyword == DIAMETER) atom->radius[i] = 0.5 * dvalue; + else if (keyword == VOLUME) atom->vfrac[i] = dvalue; + + // set shape - else if (keyword == DIAMETER) { - atom->radius[i] = 0.5 * dvalue; - if (atom->rmass_flag && atom->density_flag) - atom->rmass[i] = 4.0*PI/3.0 * - atom->radius[i]*atom->radius[i]*atom->radius[i] * atom->density[i]; + 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; + + // set rmass via density + // if radius > 0.0, assume sphere + // if shape > 0.0, assume ellipsoid + // else set rmass to density directly - // set density - // set rmass (granular) if both rmass and radius are defined - // set rmass (peri) if both rmass and vfrac are defined - - } else if (keyword == DENSITY) { - atom->density[i] = dvalue; - if (atom->rmass_flag && atom->radius_flag) - atom->rmass[i] = 4.0*PI/3.0 * - atom->radius[i]*atom->radius[i]*atom->radius[i] * - atom->density[i]; - else if (atom->rmass_flag && atom->vfrac_flag) - atom->rmass[i] = dvalue; + } 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 (keyword == VOLUME) atom->vfrac[i] = dvalue; - - // reset any or all of 3 image flags - - else if (keyword == IMAGE) { - int xbox = (atom->image[i] & 1023) - 512; - int ybox = (atom->image[i] >> 10 & 1023) - 512; - int zbox = (atom->image[i] >> 20) - 512; - if (ximageflag) xbox = ximage; - if (yimageflag) ybox = yimage; - if (zimageflag) zbox = zimage; - atom->image[i] = ((zbox + 512 & 1023) << 20) | - ((ybox + 512 & 1023) << 10) | (xbox + 512 & 1023); - - } else if (keyword == DIPOLE) { - if (atom->dipole[atom->type[i]] > 0.0) { - double **mu = atom->mu; - mu[i][0] = xvalue; - mu[i][1] = yvalue; - mu[i][2] = zvalue; - double lensq = - mu[i][0]*mu[i][0] + mu[i][1]*mu[i][1] + mu[i][2]*mu[i][2]; - double scale = atom->dipole[atom->type[i]]/sqrt(lensq); - mu[i][0] *= scale; - mu[i][1] *= scale; - mu[i][2] *= scale; - } - - } else if (keyword == QUAT) { - double PI = 4.0*atan(1.0); - 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]); - } - count++; + // reset any or all of 3 image flags + + } else if (keyword == IMAGE) { + int xbox = (atom->image[i] & 1023) - 512; + int ybox = (atom->image[i] >> 10 & 1023) - 512; + int zbox = (atom->image[i] >> 20) - 512; + if (ximageflag) xbox = ximage; + if (yimageflag) ybox = yimage; + if (zimageflag) zbox = zimage; + atom->image[i] = ((zbox + 512 & 1023) << 20) | + ((ybox + 512 & 1023) << 10) | (xbox + 512 & 1023); + + // set dipole moment + + } else if (keyword == DIPOLE) { + double **mu = atom->mu; + mu[i][0] = xvalue; + mu[i][1] = yvalue; + mu[i][2] = zvalue; + mu[i][3] = sqrt(mu[i][0]*mu[i][0] + mu[i][1]*mu[i][1] + + mu[i][2]*mu[i][2]); + + // set quaternion orientation + + } else if (keyword == QUAT) { + double PI = 4.0*atan(1.0); + 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]); } + count++; + } } /* ---------------------------------------------------------------------- @@ -397,8 +445,6 @@ void Set::setrandom(int keyword) { int i; - if (keyword == DIPOLE_RANDOM) atom->check_dipole(); - selection(atom->nlocal); RanPark *random = new RanPark(lmp,1); double **x = atom->x; @@ -422,7 +468,6 @@ void Set::setrandom(int keyword) } else if (keyword == DIPOLE_RANDOM) { int *type = atom->type; - double *dipole = atom->dipole; double **mu = atom->mu; int nlocal = atom->nlocal; @@ -431,34 +476,32 @@ void Set::setrandom(int keyword) if (domain->dimension == 3) { for (i = 0; i < nlocal; i++) if (select[i]) { - if (dipole[type[i]] > 0.0) { - random->reset(seed,x[i]); - mu[i][0] = random->uniform() - 0.5; - mu[i][1] = random->uniform() - 0.5; - mu[i][2] = random->uniform() - 0.5; - msq = mu[i][0]*mu[i][0] + mu[i][1]*mu[i][1] + mu[i][2]*mu[i][2]; - scale = dipole[type[i]]/sqrt(msq); - mu[i][0] *= scale; - mu[i][1] *= scale; - mu[i][2] *= scale; - count++; - } + random->reset(seed,x[i]); + mu[i][0] = random->uniform() - 0.5; + mu[i][1] = random->uniform() - 0.5; + mu[i][2] = random->uniform() - 0.5; + msq = mu[i][0]*mu[i][0] + mu[i][1]*mu[i][1] + mu[i][2]*mu[i][2]; + scale = dvalue/sqrt(msq); + mu[i][0] *= scale; + mu[i][1] *= scale; + mu[i][2] *= scale; + mu[i][3] = dvalue; + count++; } } else { for (i = 0; i < nlocal; i++) if (select[i]) { - if (dipole[type[i]] > 0.0) { - random->reset(seed,x[i]); - mu[i][0] = random->uniform() - 0.5; - mu[i][1] = random->uniform() - 0.5; - mu[i][2] = 0.0; - msq = mu[i][0]*mu[i][0] + mu[i][1]*mu[i][1]; - scale = dipole[type[i]]/sqrt(msq); - mu[i][0] *= scale; - mu[i][1] *= scale; - count++; - } + random->reset(seed,x[i]); + mu[i][0] = random->uniform() - 0.5; + mu[i][1] = random->uniform() - 0.5; + mu[i][2] = 0.0; + msq = mu[i][0]*mu[i][0] + mu[i][1]*mu[i][1]; + scale = dvalue/sqrt(msq); + mu[i][0] *= scale; + mu[i][1] *= scale; + mu[i][3] = dvalue; + count++; } } diff --git a/src/write_restart.cpp b/src/write_restart.cpp index e341a2c210..260dcbe3da 100644 --- a/src/write_restart.cpp +++ b/src/write_restart.cpp @@ -56,7 +56,7 @@ enum{VERSION,SMALLINT,TAGINT,BIGINT, SPECIAL_LJ_1,SPECIAL_LJ_2,SPECIAL_LJ_3, SPECIAL_COUL_1,SPECIAL_COUL_2,SPECIAL_COUL_3, XY,XZ,YZ}; -enum{MASS,SHAPE,DIPOLE}; +enum{MASS}; enum{PAIR,BOND,ANGLE,DIHEDRAL,IMPROPER}; enum{IGNORE,WARN,ERROR}; // same as thermo.cpp @@ -394,16 +394,6 @@ void WriteRestart::type_arrays() fwrite(&flag,sizeof(int),1,fp); fwrite(&atom->mass[1],sizeof(double),atom->ntypes,fp); } - if (atom->shape) { - int flag = SHAPE; - fwrite(&flag,sizeof(int),1,fp); - fwrite(&atom->shape[1][0],sizeof(double),atom->ntypes*3,fp); - } - if (atom->dipole) { - int flag = DIPOLE; - fwrite(&flag,sizeof(int),1,fp); - fwrite(&atom->dipole[1],sizeof(double),atom->ntypes,fp); - } // -1 flag signals end of type arrays