diff --git a/doc/src/Howto_body.rst b/doc/src/Howto_body.rst index 15d6aee845..41d3cd8efe 100644 --- a/doc/src/Howto_body.rst +++ b/doc/src/Howto_body.rst @@ -154,12 +154,22 @@ values consistent with the current orientation of the rigid body around its center of mass. The values are with respect to the simulation box XYZ axes, not with respect to the principal axes of the rigid body itself. LAMMPS performs the latter calculation internally. + The coordinates of each sub-particle are specified as its x,y,z displacement from the center-of-mass of the body particle. The center-of-mass position of the particle is specified by the x,y,z values in the *Atoms* section of the data file, as is the total mass of the body particle. +If the data file defines a general triclinic box, then the orientation +of the body particle and its corresponding 6 moments of inertia and +sub-particle displacements should reflect the fact the body is defined +withing a general triclinic box with edge vectors **A**,**B**,**C**. +LAMMPS will rotate the box to convert it to a restricted triclinic +box. This operation will also rotate the orientation of the body +particles. See the :doc:`Howto triclinic ` doc page +for more details. + The :doc:`pair_style body/nparticle ` command can be used with this body style to compute body/body and body/non-body interactions. @@ -226,6 +236,7 @@ values consistent with the current orientation of the rigid body around its center of mass. The values are with respect to the simulation box XYZ axes, not with respect to the principal axes of the rigid body itself. LAMMPS performs the latter calculation internally. + The coordinates of each vertex are specified as its x,y,z displacement from the center-of-mass of the body particle. The center-of-mass position of the particle is specified by the x,y,z values in the @@ -270,6 +281,15 @@ A disk, whose diameter is 3.0, mass 1.0, is specified as follows: 0 0 0 3.0 +If the data file defines a general triclinic box, then the orientation +of the body particle and its corresponding 6 moments of inertia and +polygon vertex displacements should reflect the fact the body is +defined withing a general triclinic box with edge vectors +**A**,**B**,**C**. LAMMPS will rotate the box to convert it to a +restricted triclinic box. This operation will also rotate the +orientation of the body particles. See the :doc:`Howto triclinic +` doc page for more details. + The :doc:`pair_style body/rounded/polygon ` command can be used with this body style to compute body/body interactions. The :doc:`fix wall/body/polygon ` @@ -366,6 +386,7 @@ values consistent with the current orientation of the rigid body around its center of mass. The values are with respect to the simulation box XYZ axes, not with respect to the principal axes of the rigid body itself. LAMMPS performs the latter calculation internally. + The coordinates of each vertex are specified as its x,y,z displacement from the center-of-mass of the body particle. The center-of-mass position of the particle is specified by the x,y,z values in the @@ -435,6 +456,15 @@ A sphere whose diameter is 3.0 and mass 1.0, is specified as follows: The number of edges and faces for a rod or sphere must be listed, but is ignored. +If the data file defines a general triclinic box, then the orientation +of the body particle and its corresponding 6 moments of inertia and +polyhedron vertex displacements should reflect the fact the body is +defined withing a general triclinic box with edge vectors +**A**,**B**,**C**. LAMMPS will rotate the box to convert it to a +restricted triclinic box. This operation will also rotate the +orientation of the body particles. See the :doc:`Howto triclinic +` doc page for more details. + The :doc:`pair_style body/rounded/polhedron ` command can be used with this body style to compute body/body interactions. The :doc:`fix diff --git a/src/BODY/body_nparticle.cpp b/src/BODY/body_nparticle.cpp index 62e6ee802a..05df63ff44 100644 --- a/src/BODY/body_nparticle.cpp +++ b/src/BODY/body_nparticle.cpp @@ -99,7 +99,7 @@ int BodyNparticle::unpack_border_body(AtomVecBody::Bonus *bonus, double *buf) } /* ---------------------------------------------------------------------- - populate bonus data structure with data file values + populate bonus data structure with data file values for one body ------------------------------------------------------------------------- */ void BodyNparticle::data_body(int ibonus, int ninteger, int ndouble, diff --git a/src/BODY/body_rounded_polygon.cpp b/src/BODY/body_rounded_polygon.cpp index 2fb2a991f1..b84a271779 100644 --- a/src/BODY/body_rounded_polygon.cpp +++ b/src/BODY/body_rounded_polygon.cpp @@ -155,7 +155,7 @@ int BodyRoundedPolygon::unpack_border_body(AtomVecBody::Bonus *bonus, } /* ---------------------------------------------------------------------- - populate bonus data structure with data file values + populate bonus data structure with data file values for one body ------------------------------------------------------------------------- */ void BodyRoundedPolygon::data_body(int ibonus, int ninteger, int ndouble, diff --git a/src/BODY/body_rounded_polyhedron.cpp b/src/BODY/body_rounded_polyhedron.cpp index 1d11644618..111e41dbd1 100644 --- a/src/BODY/body_rounded_polyhedron.cpp +++ b/src/BODY/body_rounded_polyhedron.cpp @@ -185,7 +185,7 @@ int BodyRoundedPolyhedron::unpack_border_body(AtomVecBody::Bonus *bonus, double } /* ---------------------------------------------------------------------- - populate bonus data structure with data file values + populate bonus data structure with data file values for one body ------------------------------------------------------------------------- */ void BodyRoundedPolyhedron::data_body(int ibonus, int ninteger, int ndouble, diff --git a/src/atom_vec_body.cpp b/src/atom_vec_body.cpp index fd590ad5fe..108d94bbec 100644 --- a/src/atom_vec_body.cpp +++ b/src/atom_vec_body.cpp @@ -597,6 +597,15 @@ void AtomVecBody::pack_data_pre(int ilocal) body[ilocal] = 1; } +/* ---------------------------------------------------------------------- + unmodify values packed by AtomVec::pack_data() +------------------------------------------------------------------------- */ + +void AtomVecBody::pack_data_post(int ilocal) +{ + body[ilocal] = body_flag; +} + /* ---------------------------------------------------------------------- pack bonus body info for writing to data file if buf is nullptr, just return buffer size @@ -631,12 +640,80 @@ void AtomVecBody::write_data_bonus(FILE *fp, int n, double *buf, int /*flag*/) } /* ---------------------------------------------------------------------- - unmodify values packed by AtomVec::pack_data() + convert read_data file info from general to restricted triclinic + parent class operates on data from Velocities section of data file + child class operates on body quaternion ------------------------------------------------------------------------- */ -void AtomVecBody::pack_data_post(int ilocal) +void AtomVecBody::read_data_general_to_restricted(int nlocal_previous, int nlocal) { - body[ilocal] = body_flag; + int j; + + AtomVec::read_data_general_to_restricted(nlocal_previous, nlocal); + + // quat_g2r = quat that rotates from general to restricted triclinic + // quat_new = body quat converted to restricted triclinic + + double quat_g2r[4],quat_new[4]; + MathExtra::mat_to_quat(domain->rotate_g2r,quat_g2r); + + for (int i = nlocal_previous; i < nlocal; i++) { + if (body[i] < 0) continue; + j = body[i]; + MathExtra::quatquat(quat_g2r,bonus[j].quat,quat_new); + bonus[j].quat[0] = quat_new[0]; + bonus[j].quat[1] = quat_new[1]; + bonus[j].quat[2] = quat_new[2]; + bonus[j].quat[3] = quat_new[3]; + } +} + +/* ---------------------------------------------------------------------- + convert info output by write_data from restricted to general triclinic + parent class operates on x and data from Velocities section of data file + child class operates on body quaternion +------------------------------------------------------------------------- */ + +void AtomVecBody::write_data_restricted_to_general() +{ + AtomVec::write_data_restricted_to_general(); + + memory->create(quat_hold,nlocal_bonus,4,"atomvec:quat_hold"); + + for (int i = 0; i < nlocal_bonus; i++) + memcpy(quat_hold[i],bonus[i].quat,4*sizeof(double)); + + // quat_r2g = quat that rotates from restricted to general triclinic + // quat_new = ellipsoid quat converted to general triclinic + + double quat_r2g[4],quat_new[4]; + MathExtra::mat_to_quat(domain->rotate_r2g,quat_r2g); + + for (int i = 0; i < nlocal_bonus; i++) { + MathExtra::quatquat(quat_r2g,bonus[i].quat,quat_new); + bonus[i].quat[0] = quat_new[0]; + bonus[i].quat[1] = quat_new[1]; + bonus[i].quat[2] = quat_new[2]; + bonus[i].quat[3] = quat_new[3]; + } +} + +/* ---------------------------------------------------------------------- + restore info output by write_data to restricted triclinic + original data is in "hold" arrays + parent class operates on x and data from Velocities section of data file + child class operates on body quaternion +------------------------------------------------------------------------- */ + +void AtomVecBody::write_data_restore_restricted() +{ + AtomVec::write_data_restore_restricted(); + + for (int i = 0; i < nlocal_bonus; i++) + memcpy(bonus[i].quat,quat_hold[i],4*sizeof(double)); + + memory->destroy(quat_hold); + quat_hold = nullptr; } /* ---------------------------------------------------------------------- diff --git a/src/atom_vec_body.h b/src/atom_vec_body.h index e02fd3bbb0..5ba90a6dc1 100644 --- a/src/atom_vec_body.h +++ b/src/atom_vec_body.h @@ -66,6 +66,10 @@ class AtomVecBody : public AtomVec { int pack_data_bonus(double *, int) override; void write_data_bonus(FILE *, int, double *, int) override; + void read_data_general_to_restricted(int, int); + void write_data_restricted_to_general(); + void write_data_restore_restricted(); + // methods used by other classes to query/set body info double radius_body(int, int, int *, double *); @@ -77,6 +81,7 @@ class AtomVecBody : public AtomVec { int *body; double *rmass, *radius; double **angmom; + double **quat_hold; int nghost_bonus, nmax_bonus; int intdoubleratio; // sizeof(double) / sizeof(int) @@ -87,7 +92,6 @@ class AtomVecBody : public AtomVec { void grow_bonus(); void copy_bonus_all(int, int); - // check(int); }; } // namespace LAMMPS_NS diff --git a/src/write_data.cpp b/src/write_data.cpp index e945f85c3b..06512a0f20 100644 --- a/src/write_data.cpp +++ b/src/write_data.cpp @@ -76,7 +76,8 @@ void WriteData::command(int narg, char **arg) lmapflag = 1; // store current (default) setting since we may change it - + + int domain_triclinic_general = domain->triclinic_general; int types_style = atom->types_style; int noinit = 0; @@ -98,9 +99,6 @@ void WriteData::command(int narg, char **arg) fixflag = 0; iarg++; } else if (strcmp(arg[iarg],"triclinic/general") == 0) { - if (!domain->triclinic_general) - error->all(FLERR,"Write_data triclinic/general for system " - "that is not general triclinic"); triclinic_general = 1; iarg++; } else if (strcmp(arg[iarg],"nolabelmap") == 0) { @@ -115,6 +113,14 @@ void WriteData::command(int narg, char **arg) } else error->all(FLERR,"Unknown write_data keyword: {}", arg[iarg]); } + // temporarily disable domain->triclinic_general if output not requested + + if (triclinic_general && !domain->triclinic_general) + error->all(FLERR,"Write_data triclinic/general for system " + "that is not general triclinic"); + if (!triclinic_general && domain->triclinic_general) + domain->triclinic_general = 0; + // init entire system since comm->exchange is done // comm::init needs neighbor::init needs pair::init needs kspace::init, etc // exception is when called by -r command-line switch @@ -145,8 +151,9 @@ void WriteData::command(int narg, char **arg) write(file); - // restore saved setting - + // restore saved settings + + domain->triclinic_general = domain_triclinic_general; atom->types_style = types_style; } @@ -216,10 +223,10 @@ void WriteData::write(const std::string &file) if (coeffflag) force_fields(); } - // if general triclinic requested: + // if general triclinic output: // reset internal per-atom data that needs rotation - if (triclinic_general) atom->avec->write_data_restricted_to_general(); + if (domain->triclinic_general) atom->avec->write_data_restricted_to_general(); // per atom info in Atoms and Velocities sections @@ -250,10 +257,10 @@ void WriteData::write(const std::string &file) if (ifix->wd_section) for (int m = 0; m < ifix->wd_section; m++) fix(ifix,m); - // if general triclinic requested: + // if general triclinic output: // restore internal per-atom data that was rotated - if (triclinic_general) atom->avec->write_data_restore_restricted(); + if (domain->triclinic_general) atom->avec->write_data_restore_restricted(); // close data file @@ -312,7 +319,7 @@ void WriteData::header() // box info: orthogonal, restricted triclinic, or general triclinic (if requested) - if (!triclinic_general) { + if (!domain->triclinic_general) { fmt::print(fp,"\n{} {} xlo xhi\n{} {} ylo yhi\n{} {} zlo zhi\n", domain->boxlo[0],domain->boxhi[0], domain->boxlo[1],domain->boxhi[1], @@ -320,7 +327,7 @@ void WriteData::header() if (domain->triclinic) fmt::print(fp,"{} {} {} xy xz yz\n",domain->xy,domain->xz,domain->yz); - } else if (triclinic_general) { + } else if (domain->triclinic_general) { fmt::print(fp,"\n{} {} {} avec\n{} {} {} bvec\n{} {} {} cvec\n", domain->avec[0],domain->avec[1],domain->avec[2], domain->bvec[0],domain->bvec[1],domain->bvec[2],