ICODE-ADC/lammps
4
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

changes for body particles in general triclinic

Browse Source
This commit is contained in:
Steve Plimpton
2023-11-14 17:03:56 -07:00
parent 4da49c6d85
commit a0a21fab64
7 changed files with 137 additions and 19 deletions
Download Patch File Download Diff File Expand all files Collapse all files

30
doc/src/Howto_body.rst
View File

@ -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 <Howto_triclininc>` doc page
for more details.
The :doc:`pair_style body/nparticle <pair_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
<Howto_triclininc>` doc page for more details.
The :doc:`pair_style body/rounded/polygon <pair_body_rounded_polygon>`
command can be used with this body style to compute body/body
interactions. The :doc:`fix wall/body/polygon <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
<Howto_triclininc>` doc page for more details.
The :doc:`pair_style body/rounded/polhedron
<pair_body_rounded_polyhedron>` command can be used with this body
style to compute body/body interactions. The :doc:`fix

2
src/BODY/body_nparticle.cpp
View File

@ -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,

2
src/BODY/body_rounded_polygon.cpp
View File

@ -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,

2
src/BODY/body_rounded_polyhedron.cpp
View File

@ -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,

83
src/atom_vec_body.cpp
View File

@ -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;
}
/* ----------------------------------------------------------------------

6
src/atom_vec_body.h
View File

@ -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

31
src/write_data.cpp
View File

@ -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],