Merge pull request #962 from ndtrung81/body-dem

Discrete element models for the BODY package
This commit is contained in:
Steve Plimpton 2018-07-20 14:37:41 -06:00 committed by GitHub
commit 5c21d2aff9
No known key found for this signature in database
GPG Key ID: 4AEE18F83AFDEB23
43 changed files with 9777 additions and 67 deletions

Binary file not shown.

After

Width:  |  Height:  |  Size: 144 KiB

View File

@ -0,0 +1,13 @@
\documentstyle[12pt]{article}
\begin{document}
\begin{eqnarray*}
F_n &=& k_n \delta_n - c_n v_n, \qquad \delta_n \le 0 \\
&=& -k_{na} \delta_n - c_n v_n, \qquad 0 < \delta_n \le r_c \\
&=& 0 \qquad \qquad \qquad \qquad \delta_n > r_c \\
F_t &=& \mu k_n \delta_n - c_t v_t, \qquad \delta_n \le 0 \\
&=& 0 \qquad \qquad \qquad \qquad \delta_n > 0
\end{eqnarray*}
\end{document}

Binary file not shown.

After

Width:  |  Height:  |  Size: 35 KiB

View File

@ -678,6 +678,8 @@ USER-INTEL, k = KOKKOS, o = USER-OMP, t = OPT.
"vector"_fix_vector.html,
"viscosity"_fix_viscosity.html,
"viscous"_fix_viscous.html,
"wall/body/polygon"_fix_wall_body_polygon.html,
"wall/body/polyhedron"_fix_wall_body_polyhedron.html,
"wall/colloid"_fix_wall.html,
"wall/gran"_fix_wall_gran.html,
"wall/gran/region"_fix_wall_gran_region.html,
@ -930,7 +932,9 @@ KOKKOS, o = USER-OMP, t = OPT.
"airebo (oi)"_pair_airebo.html,
"airebo/morse (oi)"_pair_airebo.html,
"beck (go)"_pair_beck.html,
"body"_pair_body.html,
"body/nparticle"_pair_body_nparticle.html,
"body/rounded/polygon"_pair_body_rounded/polygon.html,
"body/rounded/polyhedron"_pair_body_rounded/polyhedron.html,
"bop"_pair_bop.html,
"born (go)"_pair_born.html,
"born/coul/dsf"_pair_born.html,

View File

@ -27,18 +27,16 @@ styles supported by LAMMPS are as follows. The name in the first
column is used as the {bstyle} argument for the "atom_style
body"_atom_style.html command.
{nparticle} | rigid body with N sub-particles |
{rounded/polygon} | 2d convex polygon with N vertices :tb(c=2,s=|)
{nparticle} : rigid body with N sub-particles
{rounded/polygon} : 2d polygons with N vertices
{rounded/polyhedron} : 3d polyhedra with N vertices, E edges and F faces :tb(s=:)
The body style determines what attributes are stored for each body and
thus how they can be used to compute pairwise body/body or
bond/non-body (point particle) interactions. More details of each
style are described below.
NOTE: The rounded/polygon style listed in the table above and
described below has not yet been relesed in LAMMPS. It will be soon.
We hope to add more styles in the future. See "Section
More styles may be added in the future. See "Section
10.12"_Section_modify.html#mod_12 for details on how to add a new body
style to the code.
@ -61,7 +59,7 @@ the simple particles.
By contrast, when body particles are used, LAMMPS treats an entire
body as a single particle for purposes of computing pairwise
interactions, building neighbor lists, migrating particles between
processors, outputting particles to a dump file, etc. This means that
processors, output of particles to a dump file, etc. This means that
interactions between pairs of bodies or between a body and non-body
(point) particle need to be encoded in an appropriate pair style. If
such a pair style were to mimic the "fix rigid"_fix_rigid.html model,
@ -72,17 +70,20 @@ single body/body interaction was computed.
Thus it only makes sense to use body particles and develop such a pair
style, when particle/particle interactions are more complex than what
the "fix rigid"_fix_rigid.html command can already calculate. For
example, if particles have one or more of the following attributes:
example, consider particles with one or more of the following
attributes:
represented by a surface mesh
represented by a collection of geometric entities (e.g. planes + spheres)
deformable
internal stress that induces fragmentation :ul
then the interaction between pairs of particles is likely to be more
complex than the summation of simple sub-particle interactions. An
example is contact or frictional forces between particles with planar
surfaces that inter-penetrate.
For these models, the interaction between pairs of particles is likely
to be more complex than the summation of simple pairwise interactions.
An example is contact or frictional forces between particles with
planar surfaces that inter-penetrate. Likewise, the body particle may
store internal state, such as a stress tensor used to compute a
fracture criterion.
These are additional LAMMPS commands that can be used with body
particles of different styles
@ -130,7 +131,9 @@ x1 y1 z1
...
xN yN zN :pre
N is the number of sub-particles in the body particle. M = 6 + 3*N.
where M = 6 + 3*N, and N is the number of sub-particles in the body
particle.
The integer line has a single value N. The floating point line(s)
list 6 moments of inertia followed by the coordinates of the N
sub-particles (x1 to zN) as 3N values. These values can be listed on
@ -175,15 +178,18 @@ The {bflag2} argument is ignored.
[Specifics of body style rounded/polygon:]
NOTE: Aug 2016 - This body style has not yet been added to LAMMPS.
The info below is a placeholder.
The {rounded/polygon} body style represents body particles as a 2d
polygon with a variable number of N vertices. This style can only be
used for 2d models; see the "boundary"_boundary.html command. See the
"pair_style body/rounded/polygon" doc page for a diagram of two
squares with rounded circles at the vertices. Special cases for N = 1
(circle) and N = 2 (rod with rounded ends) can also be specified.
The {rounded/polygon} body style represents body particles as a convex
polygon with a variable number N > 2 of vertices, which can only be
used for 2d models. One example use of this body style is for 2d
discrete element models, as described in "Fraige"_#Fraige. Similar to
body style {nparticle}, the atom_style body command for this body
style takes two additional arguments:
One use of this body style is for 2d discrete element models, as
described in "Fraige"_#body-Fraige.
Similar to body style {nparticle}, the atom_style body command for
this body style takes two additional arguments:
atom_style body rounded/polygon Nmin Nmax
Nmin = minimum # of vertices in any body in the system
@ -203,17 +209,20 @@ x1 y1 z1
...
xN yN zN
i j j k k ...
radius :pre
diameter :pre
N is the number of vertices in the body particle. M = 6 + 3*N + 2*N +
1. The integer line has a single value N. The floating point line(s)
where M = 6 + 3*N + 2*N + 1, and N is the number of vertices in the
body particle.
The integer line has a single value N. The floating point line(s)
list 6 moments of inertia followed by the coordinates of the N
vertices (x1 to zN) as 3N values, followed by 2N vertex indices
corresponding to the end points of the N edges, followed by a single
radius value = the smallest circle encompassing the polygon. That
last value is used to facilitate the body/body contact detection.
These floating-point values can be listed on as many lines as you
wish; see the "read_data"_read_data.html command for more details.
vertices (x1 to zN) as 3N values (with z = 0.0 for each), followed by
2N vertex indices corresponding to the end points of the N edges,
followed by a single diameter value = the rounded diameter of the
circle that surrounds each vertex. The diameter value can be different
for each body particle. These floating-point values can be listed on
as many lines as you wish; see the "read_data"_read_data.html command
for more details.
The 6 moments of inertia (ixx,iyy,izz,ixy,ixz,iyz) should be the
values consistent with the current orientation of the rigid body
@ -225,8 +234,11 @@ 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.
For example, the following information would specify a square
particles whose edge length is sqrt(2):
For example, the following information would specify a square particle
whose edge length is sqrt(2) and rounded diameter is 1.0. The
orientation of the square is aligned with the xy coordinate axes which
is consistent with the 6 moments of inertia: ixx iyy izz ixy ixz iyz =
1 1 4 0 0 0. Note that only Izz matters in 2D simulations.
3 1 27
4
@ -235,12 +247,178 @@ particles whose edge length is sqrt(2):
-0.7071 0.7071 0
0.7071 0.7071 0
0.7071 -0.7071 0
0 1 1 2 2 3 3 0
0 1
1 2
2 3
3 0
1.0 :pre
A rod in 2D, whose length is 4.0, mass 1.0, rounded at two ends
by circles of diameter 0.5, is specified as follows:
1 1 13
2
1 1 1.33333 0 0 0
-2 0 0
2 0 0
0.5 :pre
A disk, whose diameter is 3.0, mass 1.0, is specified as follows:
1 1 10
1
1 1 4.5 0 0 0
0 0 0
3.0 :pre
The "pair_style body/rounded/polygon"_pair_body_rounded_polygon.html
command can be used with this body style to compute body/body
interactions.
interactions. The "fix wall/body/polygon"_fix_wall_body_polygon.html
command can be used with this body style to compute the interaction of
body particles with a wall.
:line
[Specifics of body style rounded/polyhedron:]
The {rounded/polyhedron} body style represents body particles as a 3d
polyhedron with a variable number of N vertices, E edges and F faces.
This style can only be used for 3d models; see the
"boundary"_boundary.html command. See the "pair_style
body/rounded/polygon" doc page for a diagram of a two 2d squares with
rounded circles at the vertices. A 3d cube with rounded spheres at
the 8 vertices and 12 rounded edges would be similar. Special cases
for N = 1 (sphere) and N = 2 (rod with rounded ends) can also be
specified.
This body style is for 3d discrete element models, as described in
"Wang"_#body-Wang.
Similar to body style {rounded/polygon}, the atom_style body command
for this body style takes two additional arguments:
atom_style body rounded/polyhedron Nmin Nmax
Nmin = minimum # of vertices in any body in the system
Nmax = maximum # of vertices in any body in the system :pre
The Nmin and Nmax arguments are used to bound the size of data
structures used internally by each particle.
When the "read_data"_read_data.html command reads a data file for this
body style, the following information must be provided for each entry
in the {Bodies} section of the data file:
atom-ID 3 M
N E F
ixx iyy izz ixy ixz iyz
x1 y1 z1
...
xN yN zN
0 1
1 2
2 3
...
0 1 2 -1
0 2 3 -1
...
1 2 3 4
diameter :pre
where M = 6 + 3*N + 2*E + 4*F + 1, and N is the number of vertices in
the body particle, E = number of edges, F = number of faces.
The integer line has three values: number of vertices (N), number of
edges (E) and number of faces (F). The floating point line(s) list 6
moments of inertia followed by the coordinates of the N vertices (x1
to zN) as 3N values, followed by 2N vertex indices corresponding to
the end points of the E edges, then 4*F vertex indices defining F
faces. The last value is the diameter value = the rounded diameter of
the sphere that surrounds each vertex. The diameter value can be
different for each body particle. These floating-point values can be
listed on as many lines as you wish; see the
"read_data"_read_data.html command for more details. Because the
maxmimum vertices per face is hard-coded to be 4
(i.e. quadrilaterals), faces with more than 4 vertices need to be
split into triangles or quadrilaterals. For triangular faces, the
last vertex index should be set to -1.
The ordering of the 4 vertices within a face should follow
the right-hand rule so that the normal vector of the face points
outwards from the center of mass.
The 6 moments of inertia (ixx,iyy,izz,ixy,ixz,iyz) should be the
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
{Atoms} section of the data file.
For example, the following information would specify a cubic particle
whose edge length is 2.0 and rounded diameter is 0.5.
The orientation of the cube is aligned with the xyz coordinate axes
which is consistent with the 6 moments of inertia: ixx iyy izz ixy ixz
iyz = 0.667 0.667 0.667 0 0 0.
1 3 79
8 12 6
0.667 0.667 0.667 0 0 0
1 1 1
1 -1 1
-1 -1 1
-1 1 1
1 1 -1
1 -1 -1
-1 -1 -1
-1 1 -1
0 1
1 2
2 3
3 0
4 5
5 6
6 7
7 4
0 4
1 5
2 6
3 7
0 1 2 3
4 5 6 7
0 1 5 4
1 2 6 5
2 3 7 6
3 0 4 7
0.5 :pre
A rod in 3D, whose length is 4.0, mass 1.0 and rounded at two ends
by circles of diameter 0.5, is specified as follows:
1 1 13
2
0 1.33333 1.33333 0 0 0
-2 0 0
2 0 0
0.5 :pre
A sphere whose diameter is 3.0 and mass 1.0, is specified as follows:
1 1 10
1
0.9 0.9 0.9 0 0 0
0 0 0
3.0 :pre
The "pair_style
body/rounded/polhedron"_pair_body_rounded_polyhedron.html command can
be used with this body style to compute body/body interactions. The
"fix wall/body/polyhedron"_fix_wall_body_polygon.html command can be
used with this body style to compute the interaction of body particles
with a wall.
:line
For output purposes via the "compute
body/local"_compute_body_local.html and "dump local"_dump.html
@ -257,10 +435,10 @@ the body particle itself. These values are calculated using the
current COM and orientation of the body particle.
For images created by the "dump image"_dump_image.html command, if the
{body} keyword is set, then each body particle is drawn as a convex
polygon consisting of N line segments. Note that the line segments
are drawn between the N vertices, which does not correspond exactly to
the physical extent of the body (because the "pair_style
{body} keyword is set, then each body particle is drawn as a polygon
consisting of N line segments. Note that the line segments are drawn
between the N vertices, which does not correspond exactly to the
physical extent of the body (because the "pair_style
rounded/polygon"_pair_body_rounded_polygon.html defines finite-size
spheres at those point and the line segments between the spheres are
tangent to the spheres). The drawn diameter of each line segment is
@ -269,6 +447,10 @@ determined by the {bflag1} parameter for the {body} keyword. The
:line
:link(Fraige)
:link(body-Fraige)
[(Fraige)] F. Y. Fraige, P. A. Langston, A. J. Matchett, J. Dodds,
Particuology, 6, 455 (2008).
:link(body-Wang)
[(Wang)] J. Wang, H. S. Yu, P. A. Langston, F. Y. Fraige, Granular
Matter, 13, 1 (2011).

View File

@ -0,0 +1,104 @@
"LAMMPS WWW Site"_lws - "LAMMPS Documentation"_ld - "LAMMPS Commands"_lc :c
:link(lws,http://lammps.sandia.gov)
:link(ld,Manual.html)
:link(lc,Section_commands.html#comm)
:line
fix wall/body/polygon command :h3
[Syntax:]
fix ID group-ID wall/body/polygon k_n c_n c_t wallstyle args keyword values ... :pre
ID, group-ID are documented in "fix"_fix.html command :ulb,l
wall/body/polygon = style name of this fix command :l
k_n = normal repulsion strength (force/distance or pressure units) :l
c_n = normal damping coefficient (force/distance or pressure units) :l
c_t = tangential damping coefficient (force/distance or pressure units) :l
wallstyle = {xplane} or {yplane} or {zplane} or {zcylinder} :l
args = list of arguments for a particular style :l
{xplane} or {yplane} args = lo hi
lo,hi = position of lower and upper plane (distance units), either can be NULL)
{zcylinder} args = radius
radius = cylinder radius (distance units) :pre
zero or more keyword/value pairs may be appended to args :l
keyword = {wiggle} :l
{wiggle} values = dim amplitude period
dim = {x} or {y} or {z}
amplitude = size of oscillation (distance units)
period = time of oscillation (time units) :pre
:ule
[Examples:]
fix 1 all wall/body/polygon 1000.0 20.0 5.0 xplane -10.0 10.0
[Description:]
This fix is for use with 2d models of body particles of style
{rounded/polygon}. It bounds the simulation domain with wall(s). All
particles in the group interact with the wall when they are close
enough to touch it. The nature of the interaction between the wall
and the polygon particles is the same as that between the polygon
particles themselves, which is similar to a Hookean potential. See
"Section 6.14"_Section_howto.html#howto_14 of the manual and the
"body"_body.html doc page for more details on using body particles.
The parameters {k_n}, {c_n}, {c_t} have the same meaning and units as
those specified with the "pair_style
body/rounded/polygon"_pair_body_rounded_polygon.html command.
The {wallstyle} can be planar or cylindrical. The 2 planar options
specify a pair of walls in a dimension. Wall positions are given by
{lo} and {hi}. Either of the values can be specified as NULL if a
single wall is desired. For a {zcylinder} wallstyle, the cylinder's
axis is at x = y = 0.0, and the radius of the cylinder is specified.
Optionally, the wall can be moving, if the {wiggle} keyword is
appended.
For the {wiggle} keyword, the wall oscillates sinusoidally, similar to
the oscillations of particles which can be specified by the "fix
move"_fix_move.html command. This is useful in packing simulations of
particles. The arguments to the {wiggle} keyword specify a dimension
for the motion, as well as it's {amplitude} and {period}. Note that
if the dimension is in the plane of the wall, this is effectively a
shearing motion. If the dimension is perpendicular to the wall, it is
more of a shaking motion. A {zcylinder} wall can only be wiggled in
the z dimension.
Each timestep, the position of a wiggled wall in the appropriate {dim}
is set according to this equation:
position = coord + A - A cos (omega * delta) :pre
where {coord} is the specified initial position of the wall, {A} is
the {amplitude}, {omega} is 2 PI / {period}, and {delta} is the time
elapsed since the fix was specified. The velocity of the wall is set
to the derivative of this expression.
[Restart, fix_modify, output, run start/stop, minimize info:]
None of the "fix_modify"_fix_modify.html options are relevant to this
fix. No global or per-atom quantities are stored by this fix for
access by various "output commands"_Section_howto.html#howto_15. No
parameter of this fix can be used with the {start/stop} keywords of
the "run"_run.html command. This fix is not invoked during "energy
minimization"_minimize.html.
[Restrictions:]
This fix is part of the BODY package. It is only enabled if LAMMPS
was built with that package. See the "Making
LAMMPS"_Section_start.html#start_3 section for more info.
Any dimension (xy) that has a wall must be non-periodic.
[Related commands:]
"atom_style body"_atom_style.html, "pair_style
body/rounded/polygon"_pair_body_rounded_polygon.html
[Default:] none

View File

@ -0,0 +1,103 @@
"LAMMPS WWW Site"_lws - "LAMMPS Documentation"_ld - "LAMMPS Commands"_lc :c
:link(lws,http://lammps.sandia.gov)
:link(ld,Manual.html)
:link(lc,Section_commands.html#comm)
:line
fix wall/body/polyhedron command :h3
[Syntax:]
fix ID group-ID wall/body/polyhedron k_n c_n c_t wallstyle args keyword values ... :pre
ID, group-ID are documented in "fix"_fix.html command :ulb,l
wall/body/polyhedron = style name of this fix command :l
k_n = normal repulsion strength (force/distance units or pressure units - see discussion below) :l
c_n = normal damping coefficient (force/distance units or pressure units - see discussion below) :l
c_t = tangential damping coefficient (force/distance units or pressure units - see discussion below) :l
wallstyle = {xplane} or {yplane} or {zplane} or {zcylinder} :l
args = list of arguments for a particular style :l
{xplane} or {yplane} args = lo hi
lo,hi = position of lower and upper plane (distance units), either can be NULL)
{zcylinder} args = radius
radius = cylinder radius (distance units) :pre
zero or more keyword/value pairs may be appended to args :l
keyword = {wiggle} :l
{wiggle} values = dim amplitude period
dim = {x} or {y} or {z}
amplitude = size of oscillation (distance units)
period = time of oscillation (time units) :pre
:ule
[Examples:]
fix 1 all wall/body/polyhedron 1000.0 20.0 5.0 xplane -10.0 10.0
[Description:]
This fix is for use with 3d models of body particles of style
{rounded/polyhedron}. It bounds the simulation domain with wall(s).
All particles in the group interact with the wall when they are close
enough to touch it. The nature of the interaction between the wall
and the polygon particles is the same as that between the polygon
particles themselves, which is similar to a Hookean potential. See
"Section 6.14"_Section_howto.html#howto_14 of the manual and the
"body"_body.html doc page for more details on using body particles.
The parameters {k_n}, {c_n}, {c_t} have the same meaning and units as
those specified with the "pair_style
body/rounded/polyhedron"_pair_body_rounded_polyhedron.html command.
The {wallstyle} can be planar or cylindrical. The 3 planar options
specify a pair of walls in a dimension. Wall positions are given by
{lo} and {hi}. Either of the values can be specified as NULL if a
single wall is desired. For a {zcylinder} wallstyle, the cylinder's
axis is at x = y = 0.0, and the radius of the cylinder is specified.
Optionally, the wall can be moving, if the {wiggle} keyword is appended.
For the {wiggle} keyword, the wall oscillates sinusoidally, similar to
the oscillations of particles which can be specified by the "fix
move"_fix_move.html command. This is useful in packing simulations of
particles. The arguments to the {wiggle} keyword specify a dimension
for the motion, as well as it's {amplitude} and {period}. Note that
if the dimension is in the plane of the wall, this is effectively a
shearing motion. If the dimension is perpendicular to the wall, it is
more of a shaking motion. A {zcylinder} wall can only be wiggled in
the z dimension.
Each timestep, the position of a wiggled wall in the appropriate {dim}
is set according to this equation:
position = coord + A - A cos (omega * delta) :pre
where {coord} is the specified initial position of the wall, {A} is
the {amplitude}, {omega} is 2 PI / {period}, and {delta} is the time
elapsed since the fix was specified. The velocity of the wall is set
to the derivative of this expression.
[Restart, fix_modify, output, run start/stop, minimize info:]
None of the "fix_modify"_fix_modify.html options are relevant to this
fix. No global or per-atom quantities are stored by this fix for
access by various "output commands"_Section_howto.html#howto_15. No
parameter of this fix can be used with the {start/stop} keywords of
the "run"_run.html command. This fix is not invoked during "energy
minimization"_minimize.html.
[Restrictions:]
This fix is part of the BODY package. It is only enabled if
LAMMPS was built with that package. See the "Making
LAMMPS"_Section_start.html#start_3 section for more info.
Any dimension (xyz) that has a wall must be non-periodic.
[Related commands:]
"atom_style body"_atom_style.html, "pair_style
body/rounded/polyhedron"_pair_body_rounded_polyhedron.html
[Default:] none

View File

@ -283,6 +283,8 @@ fix_vector.html
fix_viscosity.html
fix_viscous.html
fix_wall.html
fix_wall_body_polygon.html
fix_wall_body_polyhedron.html
fix_wall_ees.html
fix_wall_gran.html
fix_wall_gran_region.html
@ -424,8 +426,9 @@ pair_agni.html
pair_airebo.html
pair_awpmd.html
pair_beck.html
pair_body.html
pair_body_nparticle.html
pair_body_rounded_polygon.html
pair_body_rounded_polyhedron.html
pair_bop.html
pair_born.html
pair_brownian.html

View File

@ -10,21 +10,21 @@ pair_style body command :h3
[Syntax:]
pair_style body cutoff :pre
pair_style body/nparticle cutoff :pre
cutoff = global cutoff for interactions (distance units)
[Examples:]
pair_style body 3.0
pair_style body/nparticle 3.0
pair_coeff * * 1.0 1.0
pair_coeff 1 1 1.0 1.5 2.5 :pre
[Description:]
Style {body} is for use with body particles and calculates pairwise
body/body interactions as well as interactions between body and
point-particles. See "Section 6.14"_Section_howto.html#howto_14
Style {body/nparticle} is for use with body particles and calculates
pairwise body/body interactions as well as interactions between body
and point-particles. See "Section 6.14"_Section_howto.html#howto_14
of the manual and the "body"_body.html doc page for more details on
using body particles.

View File

@ -8,12 +8,127 @@
pair_style body/rounded/polygon command :h3
[Syntax:]
pair_style body/rounded/polygon c_n c_t mu delta_ua cutoff :pre
c_n = normal damping coefficient
c_t = tangential damping coefficient
mu = normal friction coefficient during gross sliding
delta_ua = multiple contact scaling factor
cutoff = global separation cutoff for interactions (distance units), see below for definition :pre
[Examples:]
pair_style body/rounded/polygon 20.0 5.0 0.0 1.0 0.5
pair_coeff * * 100.0 1.0
pair_coeff 1 1 100.0 1.0 :pre
[Description:]
Note: This feature is not yet implemented.
Style {body/rounded/polygon} is for use with 2d models of body
particles of style {rounded/polygon}. It calculates pairwise
body/body interactions which can include body particles modeled as
1-vertex circular disks with a specified diameter. See "Section
6.14"_Section_howto.html#howto_14 of the manual and the
"body"_body.html doc page for more details on using body
rounded/polygon particles.
This pairwise interaction between rounded polygons is described in
"Fraige"_#pair-Fraige, where a polygon does not have sharp corners,
but is rounded at its vertices by circles centered on each vertex with
a specified diameter. The edges of the polygon are defined between
pairs of adjacent vertices. The circle diameter for each polygon is
specified in the data file read by the "read data"_read_data.html
command. This is a 2d discrete element model (DEM) which allows for
multiple contact points.
Note that when two particles interact, the effective surface of each
polygon particle is displaced outward from each of its vertices and
edges by half its circle diameter (as in the diagram below of a gray
and yellow square particle). The interaction forces and energies
between two particles are defined with respect to the separation of
their respective rounded surfaces, not by the separation of the
vertices and edges themselves.
This means that the specified cutoff in the pair_style command is the
cutoff distance, r_c, for the surface separation, \delta_n (see figure
below). This is the distance at which two particles no longer
interact. If r_c is specified as 0.0, then it is a contact-only
interaction. I.e. the two particles must overlap in order to exert a
repulsive force on each other. If r_c > 0.0, then the force between
two particles will be attractive for surface separations from 0 to
r_c, and repulsive once the particles overlap.
Note that unlike for other pair styles, the specified cutoff is not
the distance between the centers of two particles at which they stop
interacting. This center-to-center distance depends on the shape and
size of the two particles and their relative orientation. LAMMPS
takes that into account when computing the surface separation distance
and applying the r_c cutoff.
The forces between vertex-vertex, vertex-edge, and edge-edge overlaps
are given by:
:c,image(Eqs/pair_body_rounded.jpg)
:c,image(JPG/pair_body_rounded.jpg)
Note that F_n and F_t are functions of the surface separation \delta_n
= d - (R_i + R_j). In this model, when (R_i + R_j) < d < (R_i + R_j)
+ r_c, that is, 0 < \delta_n < r_c, the cohesive region of the two
surfaces overlap and the two surfaces are attractive to each other.
In "Fraige"_#pair-Fraige, the tangential friction force between two
particles that are in contact is modeled differently prior to gross
sliding (i.e. static friction) and during gross-sliding (kinetic
friction). The latter takes place when the tangential deformation
exceeds the Coulomb frictional limit. In the current implementation,
however, we do not take into account frictional history, i.e. we do
not keep track of how many time steps the two particles have been in
contact nor calculate the tangential deformation. Instead, we assume
that gross sliding takes place as soon as two particles are in
contact.
The following coefficients must be defined for each pair of atom types
via the "pair_coeff"_pair_coeff.html command as in the examples above,
or in the data file read by the "read_data"_read_data.html command:
k_n (energy/distance^2 units)
k_na (energy/distance^2 units) :ul
Effectively, k_n and k_na are the slopes of the red lines in the plot
above for force versus surface separation, for \delta_n < 0 and 0 <
\delta_n < r_c respectively.
[Mixing, shift, table, tail correction, restart, rRESPA info]:
This pair style does not support the "pair_modify"_pair_modify.html
mix, shift, table, and tail options.
This pair style does not write its information to "binary restart
files"_restart.html. Thus, you need to re-specify the pair_style and
pair_coeff commands in an input script that reads a restart file.
This pair style can only be used via the {pair} keyword of the
"run_style respa"_run_style.html command. It does not support the
{inner}, {middle}, {outer} keywords.
[Restrictions:]
These pair styles are part of the BODY package. They are only enabled
if LAMMPS was built with that package. See the "Making
LAMMPS"_Section_start.html#start_3 section for more info.
This pair style requires the "newton"_newton.html setting to be "on"
for pair interactions.
[Related commands:]
"pair_style body"_pair_body.html
"pair_coeff"_pair_coeff.html
[Default:] none
:link(pair-Fraige)
[(Fraige)] F. Y. Fraige, P. A. Langston, A. J. Matchett, J. Dodds,
Particuology, 6, 455 (2008).

View File

@ -0,0 +1,130 @@
"LAMMPS WWW Site"_lws - "LAMMPS Documentation"_ld - "LAMMPS Commands"_lc :c
:link(lws,http://lammps.sandia.gov)
:link(ld,Manual.html)
:link(lc,Section_commands.html#comm)
:line
pair_style body/rounded/polyhedron command :h3
[Syntax:]
pair_style body/rounded/polyhedron c_n c_t mu delta_ua cutoff :pre
c_n = normal damping coefficient
c_t = tangential damping coefficient
mu = normal friction coefficient during gross sliding
delta_ua = multiple contact scaling factor
cutoff = global separation cutoff for interactions (distance units), see below for definition :pre
[Examples:]
pair_style body/rounded/polyhedron 20.0 5.0 0.0 1.0 0.5
pair_coeff * * 100.0 1.0
pair_coeff 1 1 100.0 1.0 :pre
[Description:]
Style {body/rounded/polygon} is for use with 3d models of body
particles of style {rounded/polyhedron}. It calculates pairwise
body/body interactions which can include body particles modeled as
1-vertex spheres with a specified diameter. See "Section
6.14"_Section_howto.html#howto_14 of the manual and the
"body"_body.html doc page for more details on using body
rounded/polyhedron particles.
This pairwise interaction between the rounded polyhedra is described
in "Wang"_#pair-Wang, where a polyhedron does not have sharp corners
and edges, but is rounded at its vertices and edges by spheres
centered on each vertex with a specified diameter. The edges if the
polyhedron are defined between pairs of adjacent vertices. Its faces
are defined by a loop of edges. The sphere diameter for each polygon
is specified in the data file read by the "read data"_read_data.html
command. This is a discrete element model (DEM) which allows for
multiple contact points.
Note that when two particles interact, the effective surface of each
polyhedron particle is displaced outward from each of its vertices,
edges, and faces by half its sphere diameter. The interaction forces
and energies between two particles are defined with respect to the
separation of their respective rounded surfaces, not by the separation
of the vertices, edges, and faces themselves.
This means that the specified cutoff in the pair_style command is the
cutoff distance, r_c, for the surface separation, \delta_n (see figure
below). This is the distance at which two particles no longer
interact. If r_c is specified as 0.0, then it is a contact-only
interaction. I.e. the two particles must overlap in order to exert a
repulsive force on each other. If r_c > 0.0, then the force between
two particles will be attractive for surface separations from 0 to
r_c, and repulsive once the particles overlap.
Note that unlike for other pair styles, the specified cutoff is not
the distance between the centers of two particles at which they stop
interacting. This center-to-center distance depends on the shape and
size of the two particles and their relative orientation. LAMMPS
takes that into account when computing the surface separation distance
and applying the r_c cutoff.
The forces between vertex-vertex, vertex-edge, vertex-face, edge-edge,
and edge-face overlaps are given by:
:c,image(Eqs/pair_body_rounded.jpg)
:c,image(JPG/pair_body_rounded.jpg)
In "Wang"_#pair-Wang, the tangential friction force between two
particles that are in contact is modeled differently prior to gross
sliding (i.e. static friction) and during gross-sliding (kinetic
friction). The latter takes place when the tangential deformation
exceeds the Coulomb frictional limit. In the current implementation,
however, we do not take into account frictional history, i.e. we do
not keep track of how many time steps the two particles have been in
contact nor calculate the tangential deformation. Instead, we assume
that gross sliding takes place as soon as two particles are in
contact.
The following coefficients must be defined for each pair of atom types
via the "pair_coeff"_pair_coeff.html command as in the examples above,
or in the data file read by the "read_data"_read_data.html command:
k_n (energy/distance^2 units)
k_na (energy/distance^2 units) :ul
Effectively, k_n and k_na are the slopes of the red lines in the plot
above for force versus surface separation, for \delta_n < 0 and 0 <
\delta_n < r_c respectively.
[Mixing, shift, table, tail correction, restart, rRESPA info]:
This pair style does not support the "pair_modify"_pair_modify.html
mix, shift, table, and tail options.
This pair style does not write its information to "binary restart
files"_restart.html. Thus, you need to re-specify the pair_style and
pair_coeff commands in an input script that reads a restart file.
This pair style can only be used via the {pair} keyword of the
"run_style respa"_run_style.html command. It does not support the
{inner}, {middle}, {outer} keywords.
[Restrictions:]
These pair styles are part of the BODY package. They are only enabled
if LAMMPS was built with that package. See the "Making
LAMMPS"_Section_start.html#start_3 section for more info.
This pair style requires the "newton"_newton.html setting to be "on"
for pair interactions.
[Related commands:]
"pair_coeff"_pair_coeff.html
[Default:] none
:link(pair-Wang)
[(Wang)] J. Wang, H. S. Yu, P. A. Langston, F. Y. Fraige, Granular
Matter, 13, 1 (2011).

76
examples/body/data.cubes Normal file
View File

@ -0,0 +1,76 @@
LAMMPS data file for polygons: cubes, moment of inertia I = m edge^2/ 6
2 atoms
2 bodies
1 atom types
0 6 xlo xhi
0 6 ylo yhi
0 6 zlo zhi
Atoms
1 1 1 1 1.5 1.5 1.5
2 1 1 1 4.0 4.0 4.0
Bodies
1 3 79
8 12 6
0.667 0.667 0.667 0 0 0
1 1 1
1 -1 1
-1 -1 1
-1 1 1
1 1 -1
1 -1 -1
-1 -1 -1
-1 1 -1
0 1
1 2
2 3
3 0
4 5
5 6
6 7
7 4
0 4
1 5
2 6
3 7
0 1 2 3
4 5 6 7
0 1 5 4
1 2 6 5
2 3 7 6
3 0 4 7
0.5
2 3 79
8 12 6
0.667 0.667 0.667 0 0 0
1 1 1
1 -1 1
-1 -1 1
-1 1 1
1 1 -1
1 -1 -1
-1 -1 -1
-1 1 -1
0 1
1 2
2 3
3 0
4 5
5 6
6 7
7 4
0 4
1 5
2 6
3 7
0 1 2 3
4 5 6 7
0 1 5 4
1 2 6 5
2 3 7 6
3 0 4 7
0.5

32
examples/body/data.squares Executable file
View File

@ -0,0 +1,32 @@
LAMMPS data file for polygons: squares of edge length L: Izz = 1/6mL^2
2 atoms
2 bodies
1 atom types
0 12 xlo xhi
0 12 ylo yhi
-0.5 0.5 zlo zhi
Atoms
1 1 1 1 4 5 0
2 1 1 1 9 6 0
Bodies
1 1 19
4
1 1 2.67 0 0 0
-2 -2 0
-2 2 0
2 2 0
2 -2 0
0.5
2 1 19
4
1 1 2.67 0 0 0
-2 -2 0
-2 2 0
2 2 0
2 -2 0
0.5

View File

@ -8,7 +8,7 @@ read_data data.body
velocity all create 1.44 87287 loop geom
pair_style body 5.0
pair_style body/nparticle 5.0
pair_coeff * * 1.0 1.0
neighbor 0.5 bin

53
examples/body/in.cubes Normal file
View File

@ -0,0 +1,53 @@
# 3d rounded cubes
variable r index 3
variable steps index 10000
units lj
dimension 3
atom_style body rounded/polyhedron 1 10
read_data data.cubes
replicate $r $r $r
velocity all create 1.2 187287 dist gaussian mom yes rot yes
variable cut_inner equal 0.5
variable k_n equal 100
variable k_na equal 1
variable c_n equal 20
variable c_t equal 5
variable mu equal 0
variable A_ua equal 1
pair_style body/rounded/polyhedron ${c_n} ${c_t} ${mu} ${A_ua} ${cut_inner}
pair_coeff * * ${k_n} ${k_na}
comm_modify vel yes
neighbor 0.5 bin
neigh_modify every 1 delay 0 check yes
timestep 0.001
#fix 1 all nve/body
fix 1 all nvt/body temp 1.2 1.2 0.1
#fix 1 all npt/body temp 1.2 1.2 0.1 iso 0.002 0.02 1.0
compute p2 all pressure 1_temp
#compute 1 all body/local id 1 2 3
#dump 1 all local 1000 dump.* index c_1[1] c_1[2] c_1[3] c_1[4]
#dump 2 all image 1000 image.*.jpg type type &
# zoom 1.5 adiam 1.5 body type 0 0 view 60 15
#dump_modify 2 pad 6
thermo_style custom step ke pe etotal c_p2 c_1_temp
thermo 1000
run ${steps}

57
examples/body/in.pour3d Normal file
View File

@ -0,0 +1,57 @@
# pouring 3d rounded polyhedron bodies
variable steps index 6000
units lj
boundary p p fm
comm_modify vel yes
atom_style body rounded/polyhedron 1 8
atom_modify map array
region reg block 0 50 0 50 0 50 units box
create_box 4 reg
variable cut_inner equal 0.5
variable k_n equal 100
variable k_na equal 5
variable c_n equal 20
variable c_t equal 5
variable mu equal 0
variable A_ua equal 1
pair_style body/rounded/polyhedron ${c_n} ${c_t} ${mu} ${A_ua} ${cut_inner}
pair_coeff * * ${k_n} ${k_na}
neighbor 0.5 bin
neigh_modify every 1 delay 0 check yes
timestep 0.001
fix 1 all nve/body
fix 2 all gravity 1.0 spherical 0.0 -180.0
molecule object molecule.cube molecule.tetra toff 1 &
molecule.rod3d toff 2 molecule.point3d toff 3
region slab block 5 45 5 45 25 35 units box
fix ins all pour 500 0 4767548 vol 0.4 10 region slab mol object &
molfrac 0.25 0.25 0.25 0.25
fix 4 all wall/body/polyhedron 2000 50 50 zplane 0.0 NULL
#compute 1 all body/local type 1 2 3
#dump 1 all local 1000 dump.polyhedron index c_1[1] c_1[2] c_1[3] c_1[4]
#dump 10 all custom 1000 tmp.dump id type x y z radius
thermo_style custom step atoms ke pe etotal press
thermo 1000
#dump 2 all image 500 image.*.jpg type type &
# zoom 1.5 adiam 1.5 body type 0 0 view 75 15
#dump_modify 2 pad 6
run ${steps}

55
examples/body/in.squares Executable file
View File

@ -0,0 +1,55 @@
# 2d rounded polygon bodies
variable r index 4
variable steps index 100000
variable T index 0.5
variable P index 0.1
variable seed index 980411
units lj
dimension 2
atom_style body rounded/polygon 1 6
atom_modify map array
read_data data.squares
replicate $r $r 1
velocity all create $T ${seed} dist gaussian mom yes rot yes
variable cut_inner equal 0.5
variable k_n equal 100
variable k_na equal 2
variable c_n equal 1
variable c_t equal 1
variable mu equal 0.1
variable delta_ua equal 0.5
pair_style body/rounded/polygon ${c_n} ${c_t} ${mu} ${delta_ua} ${cut_inner}
pair_coeff * * ${k_n} ${k_na}
comm_modify vel yes
neighbor 0.5 bin
neigh_modify every 1 delay 0 check yes
timestep 0.001
#fix 1 all nve/body
#fix 1 all nvt/body temp $T $T 1.0
fix 1 all npt/body temp $T $T 1.0 x 0.001 $P 1.0 &
y 0.001 $P 1.0 couple xy fixedpoint 0 0 0
fix 2 all enforce2d
#compute 1 all body/local id 1 2 3
#dump 1 all local 100000 dump.polygon.* index c_1[1] c_1[2] c_1[3] c_1[4]
thermo_style custom step ke pe etotal press
thermo 1000
#dump 2 all image 10000 image.*.jpg type type zoom 2.0 &
# adiam 1.5 body type 0 0
#dump_modify 2 pad 6
run ${steps}

57
examples/body/in.wall2d Executable file
View File

@ -0,0 +1,57 @@
# 2d rounded polygon bodies
variable r index 4
variable steps index 100000
variable T index 0.5
variable P index 0.1
variable seed index 980411
units lj
dimension 2
atom_style body rounded/polygon 1 6
atom_modify map array
read_data data.squares
replicate $r $r 1
velocity all create $T ${seed} dist gaussian mom yes rot yes
change_box all boundary p f p
variable cut_inner equal 0.5
variable k_n equal 100
variable k_na equal 2
variable c_n equal 0.1
variable c_t equal 0.1
variable mu equal 0.1
variable delta_ua equal 0.5
pair_style body/rounded/polygon ${c_n} ${c_t} ${mu} ${delta_ua} ${cut_inner}
pair_coeff * * ${k_n} ${k_na}
comm_modify vel yes
neighbor 0.5 bin
neigh_modify every 1 delay 0 check yes
timestep 0.001
#fix 1 all nve/body
#fix 1 all nvt/body temp $T $T 1.0
fix 1 all npt/body temp $T $T 1.0 x 0.001 $P 1.0 fixedpoint 0 0 0
fix 2 all enforce2d
fix 3 all wall/body/polygon 2000 50 50 yplane 0.0 48.0
#compute 1 all body/local id 1 2 3
#dump 1 all local 100000 dump.polygon.* index c_1[1] c_1[2] c_1[3] c_1[4]
thermo_style custom step ke pe etotal press
thermo 1000
#dump 2 all image 10000 image.*.jpg type type zoom 2.0 &
# adiam 1.5 body type 0 0
#dump_modify 2 pad 6
run ${steps}

View File

@ -0,0 +1,125 @@
LAMMPS (29 Jun 2018)
# 3d rounded cubes
variable r index 3
variable steps index 10000
units lj
dimension 3
atom_style body rounded/polyhedron 1 10
read_data data.cubes
orthogonal box = (0 0 0) to (6 6 6)
1 by 1 by 1 MPI processor grid
reading atoms ...
2 atoms
2 bodies
replicate $r $r $r
replicate 3 $r $r
replicate 3 3 $r
replicate 3 3 3
orthogonal box = (0 0 0) to (18 18 18)
1 by 1 by 1 MPI processor grid
54 atoms
Time spent = 0.000217915 secs
velocity all create 1.2 187287 dist gaussian mom yes rot yes
variable cut_inner equal 0.5
variable k_n equal 100
variable k_na equal 1
variable c_n equal 20
variable c_t equal 5
variable mu equal 0
variable A_ua equal 1
pair_style body/rounded/polyhedron ${c_n} ${c_t} ${mu} ${A_ua} ${cut_inner}
pair_style body/rounded/polyhedron 20 ${c_t} ${mu} ${A_ua} ${cut_inner}
pair_style body/rounded/polyhedron 20 5 ${mu} ${A_ua} ${cut_inner}
pair_style body/rounded/polyhedron 20 5 0 ${A_ua} ${cut_inner}
pair_style body/rounded/polyhedron 20 5 0 1 ${cut_inner}
pair_style body/rounded/polyhedron 20 5 0 1 0.5
pair_coeff * * ${k_n} ${k_na}
pair_coeff * * 100 ${k_na}
pair_coeff * * 100 1
comm_modify vel yes
neighbor 0.5 bin
neigh_modify every 1 delay 0 check yes
timestep 0.001
#fix 1 all nve/body
fix 1 all nvt/body temp 1.2 1.2 0.1
#fix 1 all npt/body temp 1.2 1.2 0.1 iso 0.002 0.02 1.0
compute p2 all pressure 1_temp
#compute 1 all body/local id 1 2 3
#dump 1 all local 1000 dump.* index c_1[1] c_1[2] c_1[3] c_1[4]
#dump 2 all image 1000 image.*.jpg type type # zoom 1.5 adiam 1.5 body type 0 0 view 60 15
#dump_modify 2 pad 6
thermo_style custom step ke pe etotal c_p2 c_1_temp
thermo 1000
run ${steps}
run 10000
Neighbor list info ...
update every 1 steps, delay 0 steps, check yes
max neighbors/atom: 2000, page size: 100000
master list distance cutoff = 3.9641
ghost atom cutoff = 3.9641
binsize = 1.98205, bins = 10 10 10
1 neighbor lists, perpetual/occasional/extra = 1 0 0
(1) pair body/rounded/polyhedron, perpetual
attributes: half, newton on
pair build: half/bin/atomonly/newton
stencil: half/bin/3d/newton
bin: standard
Per MPI rank memory allocation (min/avg/max) = 4.952 | 4.952 | 4.952 Mbytes
Step KinEng PotEng TotEng c_p2 c_1_temp
0 1.7666667 0 1.7666667 0.01090535 0.59439252
1000 3.1462962 0.17392649 3.3202227 0.02361912 1.1654694
2000 2.9311648 0.13836102 3.0695258 0.021748224 1.1950624
3000 3.090491 0.16511199 3.255603 0.018691142 1.23672
4000 2.7401565 0.17792155 2.9180781 0.015093853 1.1180839
5000 3.0880849 0.17587085 3.2639557 0.030563042 1.2831154
6000 3.2180776 0.19732251 3.4154001 0.028338151 1.258839
7000 2.9514882 0.25088882 3.202377 0.025296925 1.1746326
8000 3.0101226 0.28825968 3.2983823 0.027273454 1.2138056
9000 3.0164253 0.1901733 3.2065986 0.033228915 1.3095914
10000 2.3780401 0.34082434 2.7188644 0.031838531 1.0208679
Loop time of 38.5686 on 1 procs for 10000 steps with 54 atoms
Performance: 22401.653 tau/day, 259.278 timesteps/s
100.0% CPU use with 1 MPI tasks x no OpenMP threads
MPI task timing breakdown:
Section | min time | avg time | max time |%varavg| %total
---------------------------------------------------------------
Pair | 38.426 | 38.426 | 38.426 | 0.0 | 99.63
Neigh | 0.0043154 | 0.0043154 | 0.0043154 | 0.0 | 0.01
Comm | 0.047616 | 0.047616 | 0.047616 | 0.0 | 0.12
Output | 0.00017595 | 0.00017595 | 0.00017595 | 0.0 | 0.00
Modify | 0.082948 | 0.082948 | 0.082948 | 0.0 | 0.22
Other | | 0.007761 | | | 0.02
Nlocal: 54 ave 54 max 54 min
Histogram: 1 0 0 0 0 0 0 0 0 0
Nghost: 96 ave 96 max 96 min
Histogram: 1 0 0 0 0 0 0 0 0 0
Neighs: 100 ave 100 max 100 min
Histogram: 1 0 0 0 0 0 0 0 0 0
Total # of neighbors = 100
Ave neighs/atom = 1.85185
Neighbor list builds = 268
Dangerous builds = 0
Total wall time: 0:00:38

View File

@ -0,0 +1,125 @@
LAMMPS (29 Jun 2018)
# 3d rounded cubes
variable r index 3
variable steps index 10000
units lj
dimension 3
atom_style body rounded/polyhedron 1 10
read_data data.cubes
orthogonal box = (0 0 0) to (6 6 6)
1 by 2 by 2 MPI processor grid
reading atoms ...
2 atoms
2 bodies
replicate $r $r $r
replicate 3 $r $r
replicate 3 3 $r
replicate 3 3 3
orthogonal box = (0 0 0) to (18 18 18)
1 by 2 by 2 MPI processor grid
54 atoms
Time spent = 0.00103807 secs
velocity all create 1.2 187287 dist gaussian mom yes rot yes
variable cut_inner equal 0.5
variable k_n equal 100
variable k_na equal 1
variable c_n equal 20
variable c_t equal 5
variable mu equal 0
variable A_ua equal 1
pair_style body/rounded/polyhedron ${c_n} ${c_t} ${mu} ${A_ua} ${cut_inner}
pair_style body/rounded/polyhedron 20 ${c_t} ${mu} ${A_ua} ${cut_inner}
pair_style body/rounded/polyhedron 20 5 ${mu} ${A_ua} ${cut_inner}
pair_style body/rounded/polyhedron 20 5 0 ${A_ua} ${cut_inner}
pair_style body/rounded/polyhedron 20 5 0 1 ${cut_inner}
pair_style body/rounded/polyhedron 20 5 0 1 0.5
pair_coeff * * ${k_n} ${k_na}
pair_coeff * * 100 ${k_na}
pair_coeff * * 100 1
comm_modify vel yes
neighbor 0.5 bin
neigh_modify every 1 delay 0 check yes
timestep 0.001
#fix 1 all nve/body
fix 1 all nvt/body temp 1.2 1.2 0.1
#fix 1 all npt/body temp 1.2 1.2 0.1 iso 0.002 0.02 1.0
compute p2 all pressure 1_temp
#compute 1 all body/local id 1 2 3
#dump 1 all local 1000 dump.* index c_1[1] c_1[2] c_1[3] c_1[4]
#dump 2 all image 1000 image.*.jpg type type # zoom 1.5 adiam 1.5 body type 0 0 view 60 15
#dump_modify 2 pad 6
thermo_style custom step ke pe etotal c_p2 c_1_temp
thermo 1000
run ${steps}
run 10000
Neighbor list info ...
update every 1 steps, delay 0 steps, check yes
max neighbors/atom: 2000, page size: 100000
master list distance cutoff = 3.9641
ghost atom cutoff = 3.9641
binsize = 1.98205, bins = 10 10 10
1 neighbor lists, perpetual/occasional/extra = 1 0 0
(1) pair body/rounded/polyhedron, perpetual
attributes: half, newton on
pair build: half/bin/atomonly/newton
stencil: half/bin/3d/newton
bin: standard
Per MPI rank memory allocation (min/avg/max) = 4.879 | 5.068 | 5.256 Mbytes
Step KinEng PotEng TotEng c_p2 c_1_temp
0 1.7666667 0 1.7666667 0.01090535 0.59439252
1000 3.1462962 0.17392649 3.3202227 0.02361912 1.1654694
2000 2.9311648 0.13836102 3.0695258 0.021748224 1.1950624
3000 3.090491 0.16511199 3.255603 0.018691142 1.23672
4000 2.7401565 0.17792155 2.9180781 0.015093853 1.1180839
5000 3.0880849 0.17587085 3.2639557 0.030563042 1.2831154
6000 3.2180776 0.19732251 3.4154001 0.028338151 1.258839
7000 2.9514882 0.25088882 3.202377 0.025296925 1.1746326
8000 3.0101226 0.28825968 3.2983823 0.027273454 1.2138056
9000 3.0164253 0.1901733 3.2065986 0.033228915 1.3095914
10000 2.3780401 0.34082434 2.7188644 0.031838531 1.0208679
Loop time of 20.5306 on 4 procs for 10000 steps with 54 atoms
Performance: 42083.509 tau/day, 487.078 timesteps/s
100.0% CPU use with 4 MPI tasks x no OpenMP threads
MPI task timing breakdown:
Section | min time | avg time | max time |%varavg| %total
---------------------------------------------------------------
Pair | 7.5288 | 10.878 | 19.952 | 159.0 | 52.98
Neigh | 0.0014424 | 0.0016552 | 0.0021195 | 0.7 | 0.01
Comm | 0.50623 | 9.5805 | 12.93 | 169.4 | 46.66
Output | 0.00011921 | 0.00014341 | 0.00021386 | 0.0 | 0.00
Modify | 0.044663 | 0.047684 | 0.05382 | 1.6 | 0.23
Other | | 0.023 | | | 0.11
Nlocal: 13.5 ave 17 max 9 min
Histogram: 1 0 0 1 0 0 0 0 1 1
Nghost: 63.5 ave 68 max 58 min
Histogram: 1 0 0 1 0 0 0 0 0 2
Neighs: 25 ave 38 max 6 min
Histogram: 1 0 0 0 0 1 0 0 1 1
Total # of neighbors = 100
Ave neighs/atom = 1.85185
Neighbor list builds = 268
Dangerous builds = 0
Total wall time: 0:00:20

View File

@ -0,0 +1,138 @@
LAMMPS (29 Jun 2018)
# pouring 3d rounded polyhedron bodies
variable steps index 6000
units lj
boundary p p fm
comm_modify vel yes
atom_style body rounded/polyhedron 1 8
atom_modify map array
region reg block 0 50 0 50 0 50 units box
create_box 4 reg
Created orthogonal box = (0 0 0) to (50 50 50)
1 by 1 by 1 MPI processor grid
variable cut_inner equal 0.5
variable k_n equal 100
variable k_na equal 5
variable c_n equal 20
variable c_t equal 5
variable mu equal 0
variable A_ua equal 1
pair_style body/rounded/polyhedron ${c_n} ${c_t} ${mu} ${A_ua} ${cut_inner}
pair_style body/rounded/polyhedron 20 ${c_t} ${mu} ${A_ua} ${cut_inner}
pair_style body/rounded/polyhedron 20 5 ${mu} ${A_ua} ${cut_inner}
pair_style body/rounded/polyhedron 20 5 0 ${A_ua} ${cut_inner}
pair_style body/rounded/polyhedron 20 5 0 1 ${cut_inner}
pair_style body/rounded/polyhedron 20 5 0 1 0.5
pair_coeff * * ${k_n} ${k_na}
pair_coeff * * 100 ${k_na}
pair_coeff * * 100 5
neighbor 0.5 bin
neigh_modify every 1 delay 0 check yes
timestep 0.001
fix 1 all nve/body
fix 2 all gravity 1.0 spherical 0.0 -180.0
molecule object molecule.cube molecule.tetra toff 1 molecule.rod3d toff 2 molecule.point3d toff 3
Read molecule object:
1 atoms with max type 1
0 bonds with max type 0
0 angles with max type 0
0 dihedrals with max type 0
0 impropers with max type 0
Read molecule object:
1 atoms with max type 2
0 bonds with max type 0
0 angles with max type 0
0 dihedrals with max type 0
0 impropers with max type 0
Read molecule object:
1 atoms with max type 3
0 bonds with max type 0
0 angles with max type 0
0 dihedrals with max type 0
0 impropers with max type 0
Read molecule object:
1 atoms with max type 4
0 bonds with max type 0
0 angles with max type 0
0 dihedrals with max type 0
0 impropers with max type 0
region slab block 5 45 5 45 25 35 units box
fix ins all pour 500 0 4767548 vol 0.4 10 region slab mol object molfrac 0.25 0.25 0.25 0.25
Particle insertion: 134 every 4472 steps, 500 by step 13417
fix 4 all wall/body/polyhedron 2000 50 50 zplane 0.0 NULL
#compute 1 all body/local type 1 2 3
#dump 1 all local 1000 dump.polyhedron index c_1[1] c_1[2] c_1[3] c_1[4]
#dump 10 all custom 1000 tmp.dump id type x y z radius
thermo_style custom step atoms ke pe etotal press
thermo 1000
#dump 2 all image 500 image.*.jpg type type # zoom 1.5 adiam 1.5 body type 0 0 view 75 15
#dump_modify 2 pad 6
run ${steps}
run 6000
Neighbor list info ...
update every 1 steps, delay 0 steps, check yes
max neighbors/atom: 2000, page size: 100000
master list distance cutoff = 5
ghost atom cutoff = 5
binsize = 2.5, bins = 20 20 20
1 neighbor lists, perpetual/occasional/extra = 1 0 0
(1) pair body/rounded/polyhedron, perpetual
attributes: half, newton on
pair build: half/bin/atomonly/newton
stencil: half/bin/3d/newton
bin: standard
Per MPI rank memory allocation (min/avg/max) = 0.5065 | 0.5065 | 0.5065 Mbytes
Step Atoms KinEng PotEng TotEng Press
0 0 -0 0 0 0
1000 134 -0 0.00083010524 0.00083010524 -2.1515152e-06
2000 134 -0 -0.00069962476 -0.00069962476 -1.4170663e-08
3000 134 -0 -0.00069962687 -0.00069962687 -4.1478181e-11
4000 134 -0 -0.00069962687 -0.00069962687 -1.2141026e-13
5000 268 -0 0.014969705 0.014969705 3.0797164e-05
6000 268 -0 0.042467887 0.042467887 0.00056148005
Loop time of 0.634737 on 1 procs for 6000 steps with 268 atoms
Performance: 816716.196 tau/day, 9452.734 timesteps/s
100.0% CPU use with 1 MPI tasks x no OpenMP threads
MPI task timing breakdown:
Section | min time | avg time | max time |%varavg| %total
---------------------------------------------------------------
Pair | 0.41391 | 0.41391 | 0.41391 | 0.0 | 65.21
Neigh | 0.010547 | 0.010547 | 0.010547 | 0.0 | 1.66
Comm | 0.0030921 | 0.0030921 | 0.0030921 | 0.0 | 0.49
Output | 0.00011492 | 0.00011492 | 0.00011492 | 0.0 | 0.02
Modify | 0.19736 | 0.19736 | 0.19736 | 0.0 | 31.09
Other | | 0.009719 | | | 1.53
Nlocal: 268 ave 268 max 268 min
Histogram: 1 0 0 0 0 0 0 0 0 0
Nghost: 3 ave 3 max 3 min
Histogram: 1 0 0 0 0 0 0 0 0 0
Neighs: 68 ave 68 max 68 min
Histogram: 1 0 0 0 0 0 0 0 0 0
Total # of neighbors = 68
Ave neighs/atom = 0.253731
Neighbor list builds = 168
Dangerous builds = 0
Total wall time: 0:00:00

View File

@ -0,0 +1,221 @@
LAMMPS (29 Jun 2018)
# 2d rounded polygon bodies
variable r index 4
variable steps index 100000
variable T index 0.5
variable P index 0.1
variable seed index 980411
units lj
dimension 2
atom_style body rounded/polygon 1 6
atom_modify map array
read_data data.squares
orthogonal box = (0 0 -0.5) to (12 12 0.5)
1 by 1 by 1 MPI processor grid
reading atoms ...
2 atoms
2 bodies
replicate $r $r 1
replicate 4 $r 1
replicate 4 4 1
orthogonal box = (0 0 -0.5) to (48 48 0.5)
1 by 1 by 1 MPI processor grid
32 atoms
Time spent = 0.00020504 secs
velocity all create $T ${seed} dist gaussian mom yes rot yes
velocity all create 0.5 ${seed} dist gaussian mom yes rot yes
velocity all create 0.5 980411 dist gaussian mom yes rot yes
variable cut_inner equal 0.5
variable k_n equal 100
variable k_na equal 2
variable c_n equal 1
variable c_t equal 1
variable mu equal 0.1
variable delta_ua equal 0.5
pair_style body/rounded/polygon ${c_n} ${c_t} ${mu} ${delta_ua} ${cut_inner}
pair_style body/rounded/polygon 1 ${c_t} ${mu} ${delta_ua} ${cut_inner}
pair_style body/rounded/polygon 1 1 ${mu} ${delta_ua} ${cut_inner}
pair_style body/rounded/polygon 1 1 0.1 ${delta_ua} ${cut_inner}
pair_style body/rounded/polygon 1 1 0.1 0.5 ${cut_inner}
pair_style body/rounded/polygon 1 1 0.1 0.5 0.5
pair_coeff * * ${k_n} ${k_na}
pair_coeff * * 100 ${k_na}
pair_coeff * * 100 2
comm_modify vel yes
neighbor 0.5 bin
neigh_modify every 1 delay 0 check yes
timestep 0.001
#fix 1 all nve/body
#fix 1 all nvt/body temp $T $T 1.0
fix 1 all npt/body temp $T $T 1.0 x 0.001 $P 1.0 y 0.001 $P 1.0 couple xy fixedpoint 0 0 0
fix 1 all npt/body temp 0.5 $T 1.0 x 0.001 $P 1.0 y 0.001 $P 1.0 couple xy fixedpoint 0 0 0
fix 1 all npt/body temp 0.5 0.5 1.0 x 0.001 $P 1.0 y 0.001 $P 1.0 couple xy fixedpoint 0 0 0
fix 1 all npt/body temp 0.5 0.5 1.0 x 0.001 0.1 1.0 y 0.001 $P 1.0 couple xy fixedpoint 0 0 0
fix 1 all npt/body temp 0.5 0.5 1.0 x 0.001 0.1 1.0 y 0.001 0.1 1.0 couple xy fixedpoint 0 0 0
fix 2 all enforce2d
#compute 1 all body/local id 1 2 3
#dump 1 all local 100000 dump.polygon.* index c_1[1] c_1[2] c_1[3] c_1[4]
thermo_style custom step ke pe etotal press
thermo 1000
#dump 2 all image 10000 image.*.jpg type type zoom 2.0 # adiam 1.5 body type 0 0
#dump_modify 2 pad 6
run ${steps}
run 100000
Neighbor list info ...
update every 1 steps, delay 0 steps, check yes
max neighbors/atom: 2000, page size: 100000
master list distance cutoff = 6.15685
ghost atom cutoff = 6.15685
binsize = 3.07843, bins = 16 16 1
1 neighbor lists, perpetual/occasional/extra = 1 0 0
(1) pair body/rounded/polygon, perpetual
attributes: half, newton on
pair build: half/bin/atomonly/newton
stencil: half/bin/2d/newton
bin: standard
Per MPI rank memory allocation (min/avg/max) = 4.781 | 4.781 | 4.781 Mbytes
Step KinEng PotEng TotEng Press
0 0.484375 0.25 0.734375 0.0067274306
1000 0.39423376 0.0017918048 0.39602557 0.0021941612
2000 0.42284177 0.01346585 0.43630762 0.0029377883
3000 0.58154405 0.011321689 0.59286574 0.003667871
4000 0.73518304 0.034603175 0.76978621 0.0018689207
5000 0.84367476 0.025292163 0.86896692 0.0089161373
6000 0.70803236 0.0085631016 0.71659546 0.0045552895
7000 0.56206452 0.10453031 0.66659483 0.010255161
8000 0.64538994 0.088817673 0.73420761 0.0037633655
9000 0.90540819 0.063696004 0.96910419 0.0077673359
10000 0.68632042 0.093265016 0.77958544 0.0057864838
11000 0.59118074 0.025654748 0.61683549 0.012518759
12000 0.67522767 0.038176401 0.71340407 0.01741153
13000 0.7644843 0.10429844 0.86878274 0.013161339
14000 0.56152694 0.067836655 0.62936359 0.016852121
15000 0.41895506 0.019513348 0.43846841 0.015225695
16000 0.55799421 0.1564559 0.71445011 0.011703561
17000 0.59391964 0.034450221 0.62836986 0.026215002
18000 0.75911858 0.030885726 0.7900043 0.018396366
19000 0.64417995 0.12110912 0.76528907 0.010247952
20000 0.57751435 0.16965651 0.74717086 0.023392323
21000 0.7613368 0.13405354 0.89539034 0.021498982
22000 0.57676692 0.18011879 0.75688571 0.024469161
23000 0.54043723 0.11842026 0.65885749 0.019799067
24000 0.62276061 0.038967924 0.66172853 0.019080086
25000 0.53157536 0.11651937 0.64809473 0.017019298
26000 0.72213293 0.039012448 0.76114538 0.015434904
27000 0.62157832 0.13697494 0.75855326 0.028711011
28000 0.41323738 0.16301101 0.57624839 0.041792632
29000 0.45774328 0.17569066 0.63343394 0.019975231
30000 0.78901796 0.099791386 0.88880934 0.024116947
31000 0.85205397 0.11977547 0.97182945 0.026667489
32000 0.37137095 0.1232622 0.49463315 0.00087637364
33000 0.26860871 0.26056381 0.52917252 0.036110517
34000 0.3018636 0.21336905 0.51523265 0.040315549
35000 0.39915129 0.28245957 0.68161085 0.034876856
36000 0.25761236 0.2352705 0.49288286 0.022772767
37000 0.1071233 0.31692858 0.42405188 0.017994666
38000 0.083729577 0.28473145 0.36846103 -0.0045370431
39000 0.070355565 0.26682083 0.33717639 0.017921556
40000 0.075894079 0.20077896 0.27667304 0.014873186
41000 0.05891028 0.15989064 0.21880092 0.025547873
42000 0.1225107 0.16583605 0.28834675 0.038842785
43000 0.17049189 0.14323991 0.3137318 0.029550161
44000 0.26823939 0.15208257 0.42032196 0.028113612
45000 0.10172203 0.1729706 0.27469264 -0.013769913
46000 0.14841355 0.19085074 0.33926429 -0.00073741985
47000 0.27654927 0.19097937 0.46752864 0.04021431
48000 0.53432331 0.080769923 0.61509323 0.029932845
49000 0.69111634 0.13064951 0.82176585 0.028985406
50000 0.24520806 0.18317453 0.42838258 0.05179746
51000 0.23541368 0.14281364 0.37822732 0.071884238
52000 0.25464996 0.095730242 0.3503802 0.034488204
53000 0.53677633 0.1058745 0.64265084 0.059932498
54000 0.32970921 0.27979128 0.60950049 0.062869716
55000 0.49094054 0.096735015 0.58767556 0.04728005
56000 0.54398249 0.2216472 0.76562969 0.056712022
57000 0.60869068 0.2338422 0.84253288 0.077143302
58000 0.72175509 0.18687368 0.90862877 0.019357656
59000 0.79442757 0.092502981 0.88693055 0.066882632
60000 0.6810555 0.077699385 0.75875488 0.095975173
61000 0.63178834 0.05071143 0.68249977 0.043586668
62000 0.76589344 0.044615704 0.81050914 0.085718411
63000 0.84815889 0.030527848 0.87868674 0.053072795
64000 0.7309043 0.051938637 0.78284294 0.058887766
65000 0.62498816 0.034474465 0.65946262 0.068446407
66000 0.69817494 0.068546004 0.76672094 0.062634433
67000 0.86444275 0.010184259 0.87462701 0.073635055
68000 0.77820319 0.0079319524 0.78613515 0.090330925
69000 0.56938919 0.0092629332 0.57865213 0.061838729
70000 0.61870712 0.010047381 0.6287545 0.066501338
71000 0.71651803 0.0088366199 0.72535465 0.079136316
72000 0.76278925 0.008828151 0.77161741 0.063672771
73000 0.75447428 0.0083985526 0.76287283 0.078256913
74000 0.66185251 0.0091910052 0.67104351 0.069840511
75000 0.58458829 0.0097671568 0.59435544 0.076123422
76000 0.7487564 0.0100022 0.7587586 0.076171741
77000 0.89505465 0.009250681 0.90430533 0.074921699
78000 0.73738164 0.0092029279 0.74658457 0.078835344
79000 0.65735281 0.010099528 0.66745233 0.077940627
80000 0.70247542 0.010306464 0.71278189 0.079560093
81000 0.74839505 0.010199092 0.75859415 0.080835104
82000 0.75193767 0.010274058 0.76221173 0.081086684
83000 0.71392598 0.010495573 0.72442156 0.082746145
84000 0.58498928 0.011027388 0.59601667 0.08356465
85000 0.59022869 0.011729474 0.60195817 0.084519397
86000 0.81753578 0.011208964 0.82874475 0.085490261
87000 0.83480682 0.010542579 0.8453494 0.086268527
88000 0.67322538 0.011170734 0.68439611 0.08751623
89000 0.62637389 0.012033316 0.6384072 0.088548094
90000 0.92828557 0.011750388 0.94003596 0.089199823
91000 0.96072564 0.010324509 0.97105015 0.090204803
92000 0.72105071 0.011484152 0.73253486 0.09140819
93000 0.65762527 0.012558219 0.67018349 0.092453474
94000 0.73991591 0.01261909 0.752535 0.093373477
95000 0.91791653 0.011980455 0.92989699 0.094182136
96000 0.76562561 0.011807085 0.7774327 0.095323684
97000 0.57292104 0.013610205 0.58653124 0.096505977
98000 0.68141076 0.013863204 0.69527396 0.097380069
99000 0.82390969 0.013002341 0.83691203 0.098235926
100000 0.77639728 0.012989342 0.78938662 0.099274147
Loop time of 3.88899 on 1 procs for 100000 steps with 32 atoms
Performance: 2221655.884 tau/day, 25713.610 timesteps/s
99.9% CPU use with 1 MPI tasks x no OpenMP threads
MPI task timing breakdown:
Section | min time | avg time | max time |%varavg| %total
---------------------------------------------------------------
Pair | 3.056 | 3.056 | 3.056 | 0.0 | 78.58
Neigh | 0.0051048 | 0.0051048 | 0.0051048 | 0.0 | 0.13
Comm | 0.091444 | 0.091444 | 0.091444 | 0.0 | 2.35
Output | 0.0011995 | 0.0011995 | 0.0011995 | 0.0 | 0.03
Modify | 0.69909 | 0.69909 | 0.69909 | 0.0 | 17.98
Other | | 0.03616 | | | 0.93
Nlocal: 32 ave 32 max 32 min
Histogram: 1 0 0 0 0 0 0 0 0 0
Nghost: 21 ave 21 max 21 min
Histogram: 1 0 0 0 0 0 0 0 0 0
Neighs: 57 ave 57 max 57 min
Histogram: 1 0 0 0 0 0 0 0 0 0
Total # of neighbors = 57
Ave neighs/atom = 1.78125
Neighbor list builds = 1445
Dangerous builds = 0
Total wall time: 0:00:03

View File

@ -0,0 +1,221 @@
LAMMPS (29 Jun 2018)
# 2d rounded polygon bodies
variable r index 4
variable steps index 100000
variable T index 0.5
variable P index 0.1
variable seed index 980411
units lj
dimension 2
atom_style body rounded/polygon 1 6
atom_modify map array
read_data data.squares
orthogonal box = (0 0 -0.5) to (12 12 0.5)
2 by 2 by 1 MPI processor grid
reading atoms ...
2 atoms
2 bodies
replicate $r $r 1
replicate 4 $r 1
replicate 4 4 1
orthogonal box = (0 0 -0.5) to (48 48 0.5)
2 by 2 by 1 MPI processor grid
32 atoms
Time spent = 0.000324011 secs
velocity all create $T ${seed} dist gaussian mom yes rot yes
velocity all create 0.5 ${seed} dist gaussian mom yes rot yes
velocity all create 0.5 980411 dist gaussian mom yes rot yes
variable cut_inner equal 0.5
variable k_n equal 100
variable k_na equal 2
variable c_n equal 1
variable c_t equal 1
variable mu equal 0.1
variable delta_ua equal 0.5
pair_style body/rounded/polygon ${c_n} ${c_t} ${mu} ${delta_ua} ${cut_inner}
pair_style body/rounded/polygon 1 ${c_t} ${mu} ${delta_ua} ${cut_inner}
pair_style body/rounded/polygon 1 1 ${mu} ${delta_ua} ${cut_inner}
pair_style body/rounded/polygon 1 1 0.1 ${delta_ua} ${cut_inner}
pair_style body/rounded/polygon 1 1 0.1 0.5 ${cut_inner}
pair_style body/rounded/polygon 1 1 0.1 0.5 0.5
pair_coeff * * ${k_n} ${k_na}
pair_coeff * * 100 ${k_na}
pair_coeff * * 100 2
comm_modify vel yes
neighbor 0.5 bin
neigh_modify every 1 delay 0 check yes
timestep 0.001
#fix 1 all nve/body
#fix 1 all nvt/body temp $T $T 1.0
fix 1 all npt/body temp $T $T 1.0 x 0.001 $P 1.0 y 0.001 $P 1.0 couple xy fixedpoint 0 0 0
fix 1 all npt/body temp 0.5 $T 1.0 x 0.001 $P 1.0 y 0.001 $P 1.0 couple xy fixedpoint 0 0 0
fix 1 all npt/body temp 0.5 0.5 1.0 x 0.001 $P 1.0 y 0.001 $P 1.0 couple xy fixedpoint 0 0 0
fix 1 all npt/body temp 0.5 0.5 1.0 x 0.001 0.1 1.0 y 0.001 $P 1.0 couple xy fixedpoint 0 0 0
fix 1 all npt/body temp 0.5 0.5 1.0 x 0.001 0.1 1.0 y 0.001 0.1 1.0 couple xy fixedpoint 0 0 0
fix 2 all enforce2d
#compute 1 all body/local id 1 2 3
#dump 1 all local 100000 dump.polygon.* index c_1[1] c_1[2] c_1[3] c_1[4]
thermo_style custom step ke pe etotal press
thermo 1000
#dump 2 all image 10000 image.*.jpg type type zoom 2.0 # adiam 1.5 body type 0 0
#dump_modify 2 pad 6
run ${steps}
run 100000
Neighbor list info ...
update every 1 steps, delay 0 steps, check yes
max neighbors/atom: 2000, page size: 100000
master list distance cutoff = 6.15685
ghost atom cutoff = 6.15685
binsize = 3.07843, bins = 16 16 1
1 neighbor lists, perpetual/occasional/extra = 1 0 0
(1) pair body/rounded/polygon, perpetual
attributes: half, newton on
pair build: half/bin/atomonly/newton
stencil: half/bin/2d/newton
bin: standard
Per MPI rank memory allocation (min/avg/max) = 4.774 | 4.774 | 4.774 Mbytes
Step KinEng PotEng TotEng Press
0 0.484375 0.25 0.734375 0.0067274306
1000 0.39423376 0.0017918048 0.39602557 0.0021941612
2000 0.42284177 0.01346585 0.43630762 0.0029377883
3000 0.58154405 0.011321689 0.59286574 0.003667871
4000 0.73518304 0.034603175 0.76978621 0.0018689207
5000 0.84367476 0.025292163 0.86896692 0.0089161373
6000 0.70803236 0.0085631016 0.71659546 0.0045552895
7000 0.56206452 0.10453031 0.66659483 0.010255161
8000 0.64538994 0.088817673 0.73420761 0.0037633655
9000 0.90540819 0.063696004 0.96910419 0.0077673359
10000 0.68632042 0.093265016 0.77958544 0.0057864837
11000 0.59118074 0.025654748 0.61683549 0.012518759
12000 0.67522767 0.038176401 0.71340407 0.01741153
13000 0.7644843 0.10429844 0.86878274 0.013161339
14000 0.56152694 0.067836656 0.6293636 0.016852113
15000 0.41895505 0.019513353 0.43846841 0.015225696
16000 0.55799443 0.15645637 0.7144508 0.011703646
17000 0.59385248 0.03451986 0.62837234 0.025482966
18000 0.75902169 0.031103586 0.79012527 0.018263354
19000 0.64266826 0.12535314 0.76802141 0.014884119
20000 0.57836261 0.16581188 0.74417449 0.024667165
21000 0.78281936 0.11877527 0.90159464 -0.0090089213
22000 0.5312006 0.13300874 0.66420934 0.025797278
23000 0.56458861 0.084369128 0.64895774 0.024630917
24000 0.65126875 0.06122992 0.71249867 0.034377198
25000 0.55173441 0.15694886 0.70868327 0.021634086
26000 0.59121615 0.17071182 0.76192797 0.024758366
27000 0.6394843 0.17442949 0.81391378 0.034919937
28000 0.31144221 0.41243036 0.72387256 0.074115225
29000 0.13516917 0.3075419 0.44271107 0.023861298
30000 0.14094934 0.24407203 0.38502137 0.037030438
31000 0.26313749 0.087395422 0.35053291 0.042347005
32000 0.51602457 0.063012079 0.57903664 0.018550299
33000 0.55628829 0.200213 0.75650129 0.026507686
34000 0.97399408 0.082504517 1.0564986 0.037889878
35000 0.64710533 0.17662002 0.82372535 0.058295508
36000 0.45769083 0.08241194 0.54010277 0.014957415
37000 0.72850105 0.053874061 0.78237512 0.037194593
38000 0.44177995 0.28939498 0.73117493 0.045194029
39000 0.46828451 0.077630686 0.54591519 0.089849009
40000 0.46786451 0.092828423 0.56069294 0.028042052
41000 0.71861856 0.097085715 0.81570427 0.036473296
42000 0.74121021 0.10553127 0.84674148 0.054058843
43000 0.62945489 0.12770673 0.75716161 0.047267994
44000 0.49900638 0.085150056 0.58415644 0.054798793
45000 0.70199572 0.063415877 0.7654116 0.038363546
46000 0.49513142 0.10649384 0.60162526 0.059392561
47000 0.3858898 0.079458749 0.46534855 0.051825764
48000 0.62585854 0.028585902 0.65444444 0.054074424
49000 0.65934482 0.51865062 1.1779954 -0.035272836
50000 0.5420438 0.082056756 0.62410056 0.031187494
51000 0.36685223 0.14224019 0.50909241 0.073790397
52000 0.19044627 0.15368389 0.34413016 0.059034266
53000 0.26847678 0.075693324 0.3441701 0.032276915
54000 0.3593711 0.19034549 0.54971659 0.070827883
55000 0.21659198 0.1929074 0.40949939 0.035916364
56000 0.28242715 0.12313241 0.40555956 0.062083926
57000 0.34067475 0.14711992 0.48779467 0.059321458
58000 0.4842796 0.16143425 0.64571385 0.059048247
59000 0.84438871 0.076546849 0.92093556 0.048046901
60000 0.92794849 0.054331626 0.98228012 0.058392272
61000 0.6916736 0.076168342 0.76784194 0.058654987
62000 0.63317965 0.094506389 0.72768604 0.061044719
63000 0.63317266 0.038785593 0.67195825 0.097236147
64000 0.81696668 0.121811 0.93877769 0.064935373
65000 0.82644758 0.25188344 1.078331 0.093352359
66000 0.64975019 0.17930857 0.82905876 0.058805254
67000 0.63487678 0.16877059 0.80364737 0.070254696
68000 0.79140717 0.11631004 0.9077172 0.064646394
69000 0.85687272 0.057835331 0.91470805 0.071057291
70000 0.67785976 0.040686768 0.71854653 0.074687222
71000 0.60594577 0.032193155 0.63813893 0.069349268
72000 0.77586745 0.024068533 0.79993598 0.083394193
73000 0.88877625 0.025746326 0.91452258 0.081511105
74000 0.73507888 0.036574786 0.77165367 0.075360233
75000 0.68787782 0.042098622 0.72997644 0.068651098
76000 0.72515745 0.04360868 0.76876613 0.069594624
77000 0.77580944 0.041826702 0.81763614 0.071937144
78000 0.76640394 0.039285046 0.80568899 0.074274921
79000 0.62504309 0.039593585 0.66463667 0.076443295
80000 0.60001642 0.043468215 0.64348464 0.094547719
81000 0.82175037 0.045608873 0.86735924 0.080186295
82000 0.85783276 0.042692576 0.90052534 0.081576548
83000 0.71367707 0.042172193 0.75584926 0.08256625
84000 0.68532406 0.044724759 0.73004882 0.083672013
85000 0.72576789 0.046982462 0.77275035 0.084789331
86000 0.75597701 0.04765086 0.80362787 0.085758056
87000 0.74190598 0.047629096 0.78953507 0.086679976
88000 0.60967704 0.049906172 0.65958321 0.085526191
89000 0.54490288 0.054768238 0.59967112 0.090604027
90000 0.75398341 0.057153453 0.81113686 0.091900858
91000 0.84577472 0.052753512 0.89852823 0.091913909
92000 0.7176235 0.050677427 0.76830093 0.092032507
93000 0.61699446 0.054097013 0.67109147 0.092071275
94000 0.76330752 0.057398618 0.82070614 0.092435043
95000 0.98754458 0.053801311 1.0413459 0.093526707
96000 0.7405897 0.052135628 0.79272533 0.095011929
97000 0.65587599 0.057011962 0.71288795 0.096692123
98000 0.72345634 0.060700171 0.78415651 0.097510345
99000 0.88283624 0.061795247 0.94463149 0.09799633
100000 0.86303812 0.058912988 0.92195111 0.09892993
Loop time of 2.80074 on 4 procs for 100000 steps with 32 atoms
Performance: 3084895.573 tau/day, 35704.810 timesteps/s
99.9% CPU use with 4 MPI tasks x no OpenMP threads
MPI task timing breakdown:
Section | min time | avg time | max time |%varavg| %total
---------------------------------------------------------------
Pair | 0.81169 | 0.89466 | 0.97669 | 8.4 | 31.94
Neigh | 0.0017524 | 0.0018129 | 0.0018773 | 0.1 | 0.06
Comm | 0.91307 | 0.99193 | 1.0691 | 7.3 | 35.42
Output | 0.00076914 | 0.00093722 | 0.0013936 | 0.0 | 0.03
Modify | 0.75335 | 0.75779 | 0.76346 | 0.4 | 27.06
Other | | 0.1536 | | | 5.48
Nlocal: 8 ave 10 max 4 min
Histogram: 1 0 0 0 0 0 1 0 0 2
Nghost: 17.25 ave 19 max 15 min
Histogram: 1 0 1 0 0 0 0 0 0 2
Neighs: 13.5 ave 21 max 5 min
Histogram: 1 0 0 0 1 0 1 0 0 1
Total # of neighbors = 54
Ave neighs/atom = 1.6875
Neighbor list builds = 1443
Dangerous builds = 0
Total wall time: 0:00:02

View File

@ -0,0 +1,223 @@
LAMMPS (29 Jun 2018)
# 2d rounded polygon bodies
variable r index 4
variable steps index 100000
variable T index 0.5
variable P index 0.1
variable seed index 980411
units lj
dimension 2
atom_style body rounded/polygon 1 6
atom_modify map array
read_data data.squares
orthogonal box = (0 0 -0.5) to (12 12 0.5)
1 by 1 by 1 MPI processor grid
reading atoms ...
2 atoms
2 bodies
replicate $r $r 1
replicate 4 $r 1
replicate 4 4 1
orthogonal box = (0 0 -0.5) to (48 48 0.5)
1 by 1 by 1 MPI processor grid
32 atoms
Time spent = 0.00029707 secs
velocity all create $T ${seed} dist gaussian mom yes rot yes
velocity all create 0.5 ${seed} dist gaussian mom yes rot yes
velocity all create 0.5 980411 dist gaussian mom yes rot yes
change_box all boundary p f p
variable cut_inner equal 0.5
variable k_n equal 100
variable k_na equal 2
variable c_n equal 0.1
variable c_t equal 0.1
variable mu equal 0.1
variable delta_ua equal 0.5
pair_style body/rounded/polygon ${c_n} ${c_t} ${mu} ${delta_ua} ${cut_inner}
pair_style body/rounded/polygon 0.1 ${c_t} ${mu} ${delta_ua} ${cut_inner}
pair_style body/rounded/polygon 0.1 0.1 ${mu} ${delta_ua} ${cut_inner}
pair_style body/rounded/polygon 0.1 0.1 0.1 ${delta_ua} ${cut_inner}
pair_style body/rounded/polygon 0.1 0.1 0.1 0.5 ${cut_inner}
pair_style body/rounded/polygon 0.1 0.1 0.1 0.5 0.5
pair_coeff * * ${k_n} ${k_na}
pair_coeff * * 100 ${k_na}
pair_coeff * * 100 2
comm_modify vel yes
neighbor 0.5 bin
neigh_modify every 1 delay 0 check yes
timestep 0.001
#fix 1 all nve/body
#fix 1 all nvt/body temp $T $T 1.0
fix 1 all npt/body temp $T $T 1.0 x 0.001 $P 1.0 fixedpoint 0 0 0
fix 1 all npt/body temp 0.5 $T 1.0 x 0.001 $P 1.0 fixedpoint 0 0 0
fix 1 all npt/body temp 0.5 0.5 1.0 x 0.001 $P 1.0 fixedpoint 0 0 0
fix 1 all npt/body temp 0.5 0.5 1.0 x 0.001 0.1 1.0 fixedpoint 0 0 0
fix 2 all enforce2d
fix 3 all wall/body/polygon 2000 50 50 yplane 0.0 48.0
#compute 1 all body/local id 1 2 3
#dump 1 all local 100000 dump.polygon.* index c_1[1] c_1[2] c_1[3] c_1[4]
thermo_style custom step ke pe etotal press
thermo 1000
#dump 2 all image 10000 image.*.jpg type type zoom 2.0 # adiam 1.5 body type 0 0
#dump_modify 2 pad 6
run ${steps}
run 100000
Neighbor list info ...
update every 1 steps, delay 0 steps, check yes
max neighbors/atom: 2000, page size: 100000
master list distance cutoff = 6.15685
ghost atom cutoff = 6.15685
binsize = 3.07843, bins = 16 16 1
1 neighbor lists, perpetual/occasional/extra = 1 0 0
(1) pair body/rounded/polygon, perpetual
attributes: half, newton on
pair build: half/bin/atomonly/newton
stencil: half/bin/2d/newton
bin: standard
Per MPI rank memory allocation (min/avg/max) = 4.771 | 4.771 | 4.771 Mbytes
Step KinEng PotEng TotEng Press
0 0.484375 0.25 0.734375 0.0067274306
1000 0.49241101 0.0031318767 0.49554289 0.017768281
2000 0.56118632 0.0026068888 0.56379321 0.003410416
3000 0.75565115 0.025578366 0.78122951 0.0071862988
4000 0.72298647 0.093150646 0.81613712 0.003190158
5000 0.51684166 0.049164868 0.56600653 0.0096960168
6000 0.56627905 0.048132853 0.6144119 0.020733586
7000 0.58122129 0.018223718 0.59944501 0.0038160759
8000 0.64297977 0.025934821 0.66891459 0.0041091784
9000 0.41748404 0.0077890042 0.42527305 0.0039270065
10000 0.35738377 0.078487805 0.43587158 3.9079782e-05
11000 0.41529308 0.13619284 0.55148592 -0.0067482285
12000 0.43274718 0.071315497 0.50406268 0.007006378
13000 0.4748331 0.069904647 0.54473775 0.0010384372
14000 0.6287791 0.12721033 0.75598943 0.0047792448
15000 0.4692413 0.12344005 0.59268136 0.018033616
16000 0.43157074 0.14306789 0.57463862 0.042356676
17000 0.53085999 0.22126296 0.75212294 0.027509646
18000 0.52688968 0.13225282 0.6591425 0.0021558013
19000 0.55032328 0.12513047 0.67545375 0.025036251
20000 0.48465097 0.1431055 0.62775647 0.017193781
21000 0.53166734 0.21928574 0.75095307 0.011564317
22000 0.62177353 0.09296159 0.71473512 0.017660922
23000 0.6972939 0.12434123 0.82163514 0.024432327
24000 0.42767372 0.22152311 0.64919684 -0.013712449
25000 0.4816037 0.19272865 0.67433236 0.052386055
26000 0.72642579 0.19697046 0.92339625 0.020407694
27000 0.39649144 0.15058326 0.5470747 0.023705766
28000 0.44896324 0.18500106 0.6339643 -0.0089410286
29000 0.5565759 0.11085772 0.66743362 0.048437166
30000 0.58173584 0.21773281 0.79946865 0.0057357773
31000 0.49199415 0.23601982 0.72801397 0.046744152
32000 0.55665496 0.20542161 0.76207658 -0.0038756805
33000 0.62730739 0.24460524 0.87191263 0.045330682
34000 0.58107044 0.16395278 0.74502322 -0.0049496051
35000 0.56838849 0.21842922 0.78681771 0.0062086036
36000 0.45910273 0.28464172 0.74374445 -0.011700747
37000 0.37092037 0.27646862 0.647389 0.022305679
38000 0.7278047 0.30674438 1.0345491 0.07698342
39000 0.5132923 0.27395066 0.78724295 0.026898634
40000 0.62348649 0.24424644 0.86773293 0.039403899
41000 0.3658401 0.15512326 0.52096337 0.022559003
42000 0.4912253 0.35712978 0.84835508 -0.010336341
43000 0.70225957 0.36314638 1.0654059 0.004148866
44000 0.56958157 0.25488927 0.82447084 0.067537066
45000 0.45854352 0.30149439 0.76003791 -0.017002401
46000 0.62787247 0.34567995 0.97355242 0.11894801
47000 0.61348914 0.29378625 0.90727539 0.067873976
48000 0.71301829 0.34135284 1.0543711 0.021077736
49000 0.53520804 0.30593196 0.84113999 0.0059257647
50000 0.44966403 0.35370793 0.80337195 0.0020395669
51000 0.5236113 0.32296924 0.84658054 -0.051011506
52000 0.53905573 0.351771 0.89082672 0.013720106
53000 0.55978158 0.41293947 0.97272106 0.068558589
54000 0.52170459 0.2718066 0.7935112 0.0093138985
55000 0.61078876 0.43353897 1.0443277 0.045377392
56000 0.51300655 0.33182278 0.84482933 -0.018418487
57000 0.54882822 0.38380093 0.93262915 0.10249946
58000 0.72106212 0.45361279 1.1746749 0.030313481
59000 0.55871447 0.63823029 1.1969448 0.019079703
60000 0.49395192 0.58283102 1.0767829 0.0179349
61000 0.45991079 0.62540573 1.0853165 0.074398804
62000 0.4655788 0.60862262 1.0742014 0.11472976
63000 0.55634524 0.63069255 1.1870378 -0.0025676135
64000 0.57688903 0.45435264 1.0312417 0.0083813852
65000 0.57168922 0.42217005 0.99385927 0.044931269
66000 0.6206044 0.46727538 1.0878798 0.019686229
67000 0.61037155 0.41840109 1.0287726 0.0195109
68000 0.63848598 0.41305347 1.0515395 0.072940144
69000 0.49244916 0.3834095 0.87585866 0.07963677
70000 0.41847062 0.51907975 0.93755037 0.18447904
71000 0.45198986 0.52973709 0.98172695 0.078419371
72000 0.47064262 0.37808165 0.84872427 -0.00046308054
73000 0.6690143 0.37549359 1.0445079 0.061208432
74000 0.60444955 0.33779636 0.94224592 -0.068840321
75000 0.61762382 0.3916421 1.0092659 0.16253292
76000 0.63657961 0.50277989 1.1393595 0.013857508
77000 0.52524028 0.43597896 0.96121924 -0.03296482
78000 0.43803533 0.33172284 0.76975817 0.078763029
79000 0.67156089 0.55272177 1.2242827 0.080822223
80000 0.68678238 0.46061627 1.1473987 0.0027036992
81000 0.64956678 0.44959229 1.0991591 0.11201483
82000 0.51060477 0.43508342 0.9456882 0.028000608
83000 0.59550548 0.69026083 1.2857663 -0.0015809004
84000 0.64222145 0.38768816 1.0299096 0.014153173
85000 0.7661229 0.43445261 1.2005755 0.048034534
86000 0.60025257 0.53027929 1.1305319 0.0056865157
87000 0.46220939 0.47470035 0.93690974 0.075311946
88000 0.54123847 0.62899839 1.1702369 0.13260162
89000 0.61212272 0.6114241 1.2235468 0.033284822
90000 0.63924773 0.6916249 1.3308726 0.045088296
91000 0.49316865 0.51037033 1.003539 0.023203598
92000 0.57572123 0.43496319 1.0106844 0.297092
93000 0.65187559 0.56815972 1.2200353 0.1538215
94000 0.64107331 0.58948521 1.2305585 0.031117778
95000 0.64584158 0.6364688 1.2823104 0.096154676
96000 0.60509093 0.601487 1.2065779 0.03457172
97000 0.68837218 0.77974186 1.468114 0.17801164
98000 0.62725266 0.64137144 1.2686241 0.17449001
99000 0.46861221 0.67000291 1.1386151 0.2429588
100000 0.5879119 0.7140612 1.3019731 0.064634257
Loop time of 2.50594 on 1 procs for 100000 steps with 32 atoms
Performance: 3447804.126 tau/day, 39905.140 timesteps/s
100.0% CPU use with 1 MPI tasks x no OpenMP threads
MPI task timing breakdown:
Section | min time | avg time | max time |%varavg| %total
---------------------------------------------------------------
Pair | 1.5639 | 1.5639 | 1.5639 | 0.0 | 62.41
Neigh | 0.0086911 | 0.0086911 | 0.0086911 | 0.0 | 0.35
Comm | 0.058926 | 0.058926 | 0.058926 | 0.0 | 2.35
Output | 0.0012379 | 0.0012379 | 0.0012379 | 0.0 | 0.05
Modify | 0.83537 | 0.83537 | 0.83537 | 0.0 | 33.34
Other | | 0.03781 | | | 1.51
Nlocal: 32 ave 32 max 32 min
Histogram: 1 0 0 0 0 0 0 0 0 0
Nghost: 20 ave 20 max 20 min
Histogram: 1 0 0 0 0 0 0 0 0 0
Neighs: 57 ave 57 max 57 min
Histogram: 1 0 0 0 0 0 0 0 0 0
Total # of neighbors = 57
Ave neighs/atom = 1.78125
Neighbor list builds = 2705
Dangerous builds = 0
Total wall time: 0:00:02

View File

@ -0,0 +1,223 @@
LAMMPS (29 Jun 2018)
# 2d rounded polygon bodies
variable r index 4
variable steps index 100000
variable T index 0.5
variable P index 0.1
variable seed index 980411
units lj
dimension 2
atom_style body rounded/polygon 1 6
atom_modify map array
read_data data.squares
orthogonal box = (0 0 -0.5) to (12 12 0.5)
2 by 2 by 1 MPI processor grid
reading atoms ...
2 atoms
2 bodies
replicate $r $r 1
replicate 4 $r 1
replicate 4 4 1
orthogonal box = (0 0 -0.5) to (48 48 0.5)
2 by 2 by 1 MPI processor grid
32 atoms
Time spent = 0.000386 secs
velocity all create $T ${seed} dist gaussian mom yes rot yes
velocity all create 0.5 ${seed} dist gaussian mom yes rot yes
velocity all create 0.5 980411 dist gaussian mom yes rot yes
change_box all boundary p f p
variable cut_inner equal 0.5
variable k_n equal 100
variable k_na equal 2
variable c_n equal 0.1
variable c_t equal 0.1
variable mu equal 0.1
variable delta_ua equal 0.5
pair_style body/rounded/polygon ${c_n} ${c_t} ${mu} ${delta_ua} ${cut_inner}
pair_style body/rounded/polygon 0.1 ${c_t} ${mu} ${delta_ua} ${cut_inner}
pair_style body/rounded/polygon 0.1 0.1 ${mu} ${delta_ua} ${cut_inner}
pair_style body/rounded/polygon 0.1 0.1 0.1 ${delta_ua} ${cut_inner}
pair_style body/rounded/polygon 0.1 0.1 0.1 0.5 ${cut_inner}
pair_style body/rounded/polygon 0.1 0.1 0.1 0.5 0.5
pair_coeff * * ${k_n} ${k_na}
pair_coeff * * 100 ${k_na}
pair_coeff * * 100 2
comm_modify vel yes
neighbor 0.5 bin
neigh_modify every 1 delay 0 check yes
timestep 0.001
#fix 1 all nve/body
#fix 1 all nvt/body temp $T $T 1.0
fix 1 all npt/body temp $T $T 1.0 x 0.001 $P 1.0 fixedpoint 0 0 0
fix 1 all npt/body temp 0.5 $T 1.0 x 0.001 $P 1.0 fixedpoint 0 0 0
fix 1 all npt/body temp 0.5 0.5 1.0 x 0.001 $P 1.0 fixedpoint 0 0 0
fix 1 all npt/body temp 0.5 0.5 1.0 x 0.001 0.1 1.0 fixedpoint 0 0 0
fix 2 all enforce2d
fix 3 all wall/body/polygon 2000 50 50 yplane 0.0 48.0
#compute 1 all body/local id 1 2 3
#dump 1 all local 100000 dump.polygon.* index c_1[1] c_1[2] c_1[3] c_1[4]
thermo_style custom step ke pe etotal press
thermo 1000
#dump 2 all image 10000 image.*.jpg type type zoom 2.0 # adiam 1.5 body type 0 0
#dump_modify 2 pad 6
run ${steps}
run 100000
Neighbor list info ...
update every 1 steps, delay 0 steps, check yes
max neighbors/atom: 2000, page size: 100000
master list distance cutoff = 6.15685
ghost atom cutoff = 6.15685
binsize = 3.07843, bins = 16 16 1
1 neighbor lists, perpetual/occasional/extra = 1 0 0
(1) pair body/rounded/polygon, perpetual
attributes: half, newton on
pair build: half/bin/atomonly/newton
stencil: half/bin/2d/newton
bin: standard
Per MPI rank memory allocation (min/avg/max) = 4.773 | 4.773 | 4.773 Mbytes
Step KinEng PotEng TotEng Press
0 0.484375 0.25 0.734375 0.0067274306
1000 0.49241101 0.0031318767 0.49554289 0.017768281
2000 0.56118632 0.0026068888 0.56379321 0.003410416
3000 0.75565115 0.025578366 0.78122951 0.0071862988
4000 0.72298647 0.093150646 0.81613712 0.003190158
5000 0.51684166 0.049164868 0.56600653 0.0096960168
6000 0.56627905 0.048132853 0.6144119 0.020733586
7000 0.58122129 0.018223718 0.59944501 0.0038160759
8000 0.64297977 0.025934821 0.66891459 0.0041091784
9000 0.41748404 0.0077890042 0.42527305 0.0039270065
10000 0.35738377 0.078487805 0.43587158 3.9079865e-05
11000 0.41529307 0.13619284 0.55148591 -0.0067482285
12000 0.43274718 0.071315527 0.50406271 0.007006369
13000 0.4748324 0.069905666 0.54473807 0.0010385254
14000 0.62603727 0.098905625 0.7249429 0.0048876764
15000 0.44512086 0.10415235 0.54927321 0.01902062
16000 0.47460177 0.18053316 0.65513493 0.045013976
17000 0.52742676 0.10110706 0.62853382 0.013615471
18000 0.46111734 0.096118795 0.55723613 0.0073676834
19000 0.59668439 0.13652292 0.73320731 0.029403553
20000 0.46840192 0.11611719 0.58451911 -0.00034412499
21000 0.53550533 0.096457461 0.6319628 0.0019785732
22000 0.46599715 0.13206373 0.59806087 0.031970672
23000 0.49280776 0.20404726 0.69685501 0.03657433
24000 0.60901688 0.18255214 0.79156902 0.044955017
25000 0.47345185 0.13671357 0.61016542 0.020313539
26000 0.47653832 0.12448225 0.60102057 0.01878099
27000 0.50008212 0.24740634 0.74748845 0.021862639
28000 0.41627204 0.2519463 0.66821834 0.054683701
29000 0.55608273 0.23100212 0.78708485 -0.0043318497
30000 0.53884537 0.3001584 0.83900377 -0.012838186
31000 0.53036238 0.2300328 0.76039518 -0.0061688449
32000 0.42666792 0.20536256 0.63203048 0.045305282
33000 0.62908185 0.1652033 0.79428515 0.0072777588
34000 0.47028154 0.388736 0.85901754 0.04332288
35000 0.54602322 0.2775624 0.82358562 0.02898206
36000 0.59860544 0.21824655 0.81685199 0.0025936194
37000 0.62467827 0.11983499 0.74451326 0.050052743
38000 0.72594229 0.36584781 1.0917901 0.04280621
39000 0.51129656 0.23859043 0.74988699 0.050817447
40000 0.53263836 0.24212889 0.77476725 0.036245922
41000 0.50288088 0.36668283 0.86956371 0.018381415
42000 0.46653688 0.21974887 0.68628574 0.012661062
43000 0.61738785 0.32131037 0.93869821 0.012709433
44000 0.56603903 0.26515554 0.83119457 0.03315102
45000 0.56231638 0.32111693 0.88343331 0.06079756
46000 0.7096208 0.2570131 0.96663391 0.048770468
47000 0.588755 0.1880748 0.7768298 0.035962604
48000 0.56296339 0.25783519 0.82079858 0.053019928
49000 0.419885 0.42328618 0.84317118 0.038105269
50000 0.63073351 0.41426285 1.0449964 0.0015271048
51000 0.59357935 0.184222 0.77780136 0.015996218
52000 0.60608471 0.36247533 0.96856003 0.10984665
53000 0.5227842 0.27686739 0.79965159 0.02761699
54000 0.39435923 0.34197355 0.73633278 0.061183263
55000 0.46748455 0.34230903 0.80979358 0.077441382
56000 0.59819827 0.29212061 0.89031889 0.043772353
57000 0.61682559 0.32788566 0.94471124 0.03992069
58000 0.52702478 0.24891506 0.77593984 0.058480883
59000 0.66925719 0.4109031 1.0801603 0.072434423
60000 0.66807714 0.39233068 1.0604078 0.082370324
61000 0.5724275 0.43308567 1.0055132 0.0072945426
62000 0.49433556 0.38453743 0.87887299 0.0036097443
63000 0.57575143 0.54067119 1.1164226 0.073339638
64000 0.68045383 0.38246533 1.0629192 0.025314593
65000 0.59843527 0.42928622 1.0277215 -0.030096445
66000 0.60274797 0.50186417 1.1046121 0.069797184
67000 0.47450407 0.52689807 1.0014021 0.008758012
68000 0.5514135 0.64113187 1.1925454 0.093863314
69000 0.52008074 0.45749565 0.97757639 -0.066061381
70000 0.69042662 0.50416006 1.1945867 0.014128617
71000 0.63925854 0.35153425 0.9907928 -0.01134957
72000 0.52088835 0.47626986 0.99715821 0.10198133
73000 0.46333852 0.5515537 1.0148922 0.00060582772
74000 0.53481418 0.50409531 1.0389095 0.00919451
75000 0.67182749 0.50380162 1.1756291 0.043301985
76000 0.70492289 0.4112122 1.1161351 0.14880484
77000 0.59781817 0.50197661 1.0997948 -0.057111711
78000 0.51677429 0.4348232 0.95159749 -0.0074619446
79000 0.50663297 0.55000424 1.0566372 0.0052071216
80000 0.59392006 0.48394003 1.0778601 -0.018990234
81000 0.66323593 0.40358336 1.0668193 -0.02961345
82000 0.61596979 0.49177944 1.1077492 0.1314853
83000 0.63917554 0.61656584 1.2557414 0.11908351
84000 0.49305291 0.46161646 0.95466937 0.033558488
85000 0.52552044 0.54250555 1.068026 0.13015174
86000 0.55140914 0.38924725 0.94065638 0.047412499
87000 0.60952504 0.52603688 1.1355619 0.039230066
88000 0.50119735 0.547539 1.0487364 0.019659933
89000 0.40331401 0.50331134 0.90662535 -0.056906034
90000 0.47067839 0.51306911 0.9837475 0.11918166
91000 0.45564995 0.38693455 0.8425845 0.12040045
92000 0.64163032 0.34232532 0.98395564 0.0057051641
93000 0.70375593 0.53646186 1.2402178 0.16044241
94000 0.53378112 0.51971406 1.0534952 0.11389004
95000 0.47055342 0.50396004 0.97451346 0.079424215
96000 0.59543473 0.40204536 0.99748009 0.096813093
97000 0.64821917 0.50051728 1.1487365 0.054071312
98000 0.55723937 0.4945909 1.0518303 0.047316424
99000 0.56044424 0.50773312 1.0681774 0.0149959
100000 0.68254229 0.32704484 1.0095871 0.0069212661
Loop time of 2.20043 on 4 procs for 100000 steps with 32 atoms
Performance: 3926501.701 tau/day, 45445.622 timesteps/s
100.0% CPU use with 4 MPI tasks x no OpenMP threads
MPI task timing breakdown:
Section | min time | avg time | max time |%varavg| %total
---------------------------------------------------------------
Pair | 0.41008 | 0.41366 | 0.41719 | 0.4 | 18.80
Neigh | 0.0027823 | 0.0030481 | 0.0034747 | 0.5 | 0.14
Comm | 0.74581 | 0.7675 | 0.78684 | 2.0 | 34.88
Output | 0.00082111 | 0.0010884 | 0.0016899 | 1.1 | 0.05
Modify | 0.83828 | 0.85329 | 0.86656 | 1.4 | 38.78
Other | | 0.1618 | | | 7.36
Nlocal: 8 ave 9 max 7 min
Histogram: 1 0 0 0 0 2 0 0 0 1
Nghost: 12.75 ave 14 max 12 min
Histogram: 2 0 0 0 0 1 0 0 0 1
Neighs: 11 ave 19 max 5 min
Histogram: 1 0 0 2 0 0 0 0 0 1
Total # of neighbors = 44
Ave neighs/atom = 1.375
Neighbor list builds = 2663
Dangerous builds = 0
Total wall time: 0:00:02

View File

@ -0,0 +1,52 @@
# 3d polygon body: cubes, moment of inertia I = m edge^2/ 6
1 atoms
3 79 body
Coords
1 0 0 0
Types
1 1
Masses
1 1.0
Body Integers
8 12 6
Body Doubles
0.667 0.667 0.667 0 0 0
1 1 1
1 -1 1
-1 -1 1
-1 1 1
1 1 -1
1 -1 -1
-1 -1 -1
-1 1 -1
0 1
1 2
2 3
3 0
4 5
5 6
6 7
7 4
0 4
1 5
2 6
3 7
0 1 2 3
4 5 6 7
0 1 5 4
1 2 6 5
2 3 7 6
3 0 4 7
0.5

View File

@ -0,0 +1,26 @@
# 2d polygon body: disks Izz = 1/2 m radius^2
1 atoms
3 10 body
Coords
1 0 0 0
Types
1 1
Masses
1 1.0
Body Integers
1 0 0
Body Doubles
1 1 1.125 0 0 0
0 0 0
3.0

View File

@ -0,0 +1,27 @@
# 2d polygon body: rods Izz = 1/12 m length^2
1 atoms
3 13 body
Coords
1 0 0 0
Types
1 1
Masses
1 1.0
Body Integers
2 1 0
Body Doubles
1 1 1.333 0 0 0
-2 0 0
2 0 0
0.5

View File

@ -0,0 +1,39 @@
# 3d polygon body: regular tetrahedra, moment of inertia I = m edge^2/ 20
1 atoms
3 47 body
Coords
1 0 0 0
Types
1 1
Masses
1 1.0
Body Integers
4 6 4
Body Doubles
0.4 0.4 0.4 0 0 0
1 1 1
1 -1 -1
-1 1 -1
-1 -1 1
0 1
0 2
0 3
1 2
2 3
3 1
0 1 2 -1
0 1 3 -1
0 2 3 -1
1 2 3 -1
0.5

View File

@ -0,0 +1,452 @@
/* ----------------------------------------------------------------------
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.
------------------------------------------------------------------------- */
/* ----------------------------------------------------------------------
Contributing author: Trung Dac Nguyen (ndactrung@gmail.com)
------------------------------------------------------------------------- */
#include <cstdlib>
#include "body_rounded_polygon.h"
#include "atom_vec_body.h"
#include "atom.h"
#include "force.h"
#include "domain.h"
#include "math_extra.h"
#include "memory.h"
#include "error.h"
using namespace LAMMPS_NS;
#define EPSILON 1.0e-7
enum{SPHERE,LINE}; // also in DumpImage
/* ---------------------------------------------------------------------- */
BodyRoundedPolygon::BodyRoundedPolygon(LAMMPS *lmp, int narg, char **arg) :
Body(lmp, narg, arg)
{
if (narg != 3) error->all(FLERR,"Invalid body rounded/polygon command");
if (domain->dimension != 2)
error->all(FLERR,"Atom_style body rounded/polygon "
"can only be used in 2d simulations");
// nmin and nmax are minimum and maximum number of vertices
int nmin = force->inumeric(FLERR,arg[1]);
int nmax = force->inumeric(FLERR,arg[2]);
if (nmin <= 0 || nmin > nmax)
error->all(FLERR,"Invalid body rounded/polygon command");
size_forward = 0;
// 1 integer for number of vertices,
// 3*nmax doubles for vertex coordinates + 2*nmax doubles for edge ends
// 1 double for the enclosing radius
// 1 double for the rounded radius
size_border = 1 + 3*nmax + 2*nmax + 1 + 1;
// NOTE: need to set appropriate nnbin param for dcp
icp = new MyPoolChunk<int>(1,1);
dcp = new MyPoolChunk<double>(3*nmin+2*nmin+1+1,3*nmax+2*nmax+1+1);
memory->create(imflag,nmax,"body/rounded/polygon:imflag");
memory->create(imdata,nmax,7,"body/nparticle:imdata");
}
/* ---------------------------------------------------------------------- */
BodyRoundedPolygon::~BodyRoundedPolygon()
{
delete icp;
delete dcp;
memory->destroy(imflag);
memory->destroy(imdata);
}
/* ---------------------------------------------------------------------- */
int BodyRoundedPolygon::nsub(AtomVecBody::Bonus *bonus)
{
return bonus->ivalue[0];
}
/* ---------------------------------------------------------------------- */
double *BodyRoundedPolygon::coords(AtomVecBody::Bonus *bonus)
{
return bonus->dvalue;
}
/* ---------------------------------------------------------------------- */
int BodyRoundedPolygon::nedges(AtomVecBody::Bonus *bonus)
{
int nvertices = bonus->ivalue[0];
if (nvertices == 1) return 0;
else if (nvertices == 2) return 1;
return nvertices;
}
/* ---------------------------------------------------------------------- */
double *BodyRoundedPolygon::edges(AtomVecBody::Bonus *bonus)
{
return bonus->dvalue+3*nsub(bonus);
}
/* ---------------------------------------------------------------------- */
double BodyRoundedPolygon::enclosing_radius(struct AtomVecBody::Bonus *bonus)
{
int nvertices = bonus->ivalue[0];
if (nvertices == 1 || nvertices == 2)
return *(bonus->dvalue+3*nsub(bonus)+2);
return *(bonus->dvalue + 3*nsub(bonus) + 2*nsub(bonus));
}
/* ---------------------------------------------------------------------- */
double BodyRoundedPolygon::rounded_radius(struct AtomVecBody::Bonus *bonus)
{
int nvertices = bonus->ivalue[0];
if (nvertices == 1 || nvertices == 2)
return *(bonus->dvalue+3*nsub(bonus)+2+1);
return *(bonus->dvalue + 3*nsub(bonus) + 2*nsub(bonus)+1);
}
/* ---------------------------------------------------------------------- */
int BodyRoundedPolygon::pack_border_body(AtomVecBody::Bonus *bonus, double *buf)
{
int nsub = bonus->ivalue[0];
buf[0] = nsub;
memcpy(&buf[1],bonus->dvalue,(3*nsub+2*nsub+1+1)*sizeof(double));
return 1+(3*nsub+2*nsub+1+1);
}
/* ---------------------------------------------------------------------- */
int BodyRoundedPolygon::unpack_border_body(AtomVecBody::Bonus *bonus,
double *buf)
{
int nsub = static_cast<int> (buf[0]);
bonus->ivalue[0] = nsub;
memcpy(bonus->dvalue,&buf[1],(3*nsub+2*nsub+1+1)*sizeof(double));
return 1+(3*nsub+2*nsub+1+1);
}
/* ----------------------------------------------------------------------
populate bonus data structure with data file values
------------------------------------------------------------------------- */
void BodyRoundedPolygon::data_body(int ibonus, int ninteger, int ndouble,
int *ifile, double *dfile)
{
AtomVecBody::Bonus *bonus = &avec->bonus[ibonus];
// set ninteger, ndouble in bonus and allocate 2 vectors of ints, doubles
if (ninteger != 1)
error->one(FLERR,"Incorrect # of integer values in "
"Bodies section of data file");
int nsub = ifile[0];
if (nsub < 1)
error->one(FLERR,"Incorrect integer value in "
"Bodies section of data file");
// nentries = number of double entries to be read from Body section:
// 6 for inertia + 3*nsub for vertex coords + 1 for rounded radius
int nentries = 6 + 3*nsub + 1;
if (ndouble != nentries)
error->one(FLERR,"Incorrect # of floating-point values in "
"Bodies section of data file");
bonus->ninteger = 1;
bonus->ivalue = icp->get(bonus->iindex);
bonus->ivalue[0] = nsub;
bonus->ndouble = 3*nsub + 2*nsub + 1 + 1;
bonus->dvalue = dcp->get(bonus->ndouble,bonus->dindex);
// diagonalize inertia tensor
double tensor[3][3];
tensor[0][0] = dfile[0];
tensor[1][1] = dfile[1];
tensor[2][2] = dfile[2];
tensor[0][1] = tensor[1][0] = dfile[3];
tensor[0][2] = tensor[2][0] = dfile[4];
tensor[1][2] = tensor[2][1] = dfile[5];
double *inertia = bonus->inertia;
double evectors[3][3];
int ierror = MathExtra::jacobi(tensor,inertia,evectors);
if (ierror) error->one(FLERR,
"Insufficient Jacobi rotations for body nparticle");
// if any principal moment < scaled EPSILON, set to 0.0
double max;
max = MAX(inertia[0],inertia[1]);
max = MAX(max,inertia[2]);
if (inertia[0] < EPSILON*max) inertia[0] = 0.0;
if (inertia[1] < EPSILON*max) inertia[1] = 0.0;
if (inertia[2] < EPSILON*max) inertia[2] = 0.0;
// exyz_space = principal axes in space frame
double ex_space[3],ey_space[3],ez_space[3];
ex_space[0] = evectors[0][0];
ex_space[1] = evectors[1][0];
ex_space[2] = evectors[2][0];
ey_space[0] = evectors[0][1];
ey_space[1] = evectors[1][1];
ey_space[2] = evectors[2][1];
ez_space[0] = evectors[0][2];
ez_space[1] = evectors[1][2];
ez_space[2] = evectors[2][2];
// enforce 3 evectors as a right-handed coordinate system
// flip 3rd vector if needed
double cross[3];
MathExtra::cross3(ex_space,ey_space,cross);
if (MathExtra::dot3(cross,ez_space) < 0.0) MathExtra::negate3(ez_space);
// create initial quaternion
MathExtra::exyz_to_q(ex_space,ey_space,ez_space,bonus->quat);
// bonus->dvalue = the first 3*nsub elements are sub-particle displacements
// find the enclosing radius of the body from the maximum displacement
int i,m;
double delta[3], rsq, erad, rrad;
double erad2 = 0;
int j = 6;
int k = 0;
for (i = 0; i < nsub; i++) {
delta[0] = dfile[j];
delta[1] = dfile[j+1];
delta[2] = dfile[j+2];
MathExtra::transpose_matvec(ex_space,ey_space,ez_space,
delta,&bonus->dvalue[k]);
rsq = delta[0] * delta[0] + delta[1] * delta[1] +
delta[2] * delta[2];
if (rsq > erad2) erad2 = rsq;
j += 3;
k += 3;
}
// .. the next 2*nsub elements are edge ends
int nedges;
if (nsub == 1) { // spheres
nedges = 0;
bonus->dvalue[k] = 0;
*(&bonus->dvalue[k]+1) = 0;
k += 2;
// the last element of bonus->dvalue is the rounded & enclosing radius
rrad = 0.5 * dfile[j];
bonus->dvalue[k] = rrad;
erad = rrad;
k++;
bonus->dvalue[k] = rrad;
atom->radius[bonus->ilocal] = erad;
} else if (nsub == 2) { // rods
nedges = 1;
for (i = 0; i < nedges; i++) {
bonus->dvalue[k] = 0;
*(&bonus->dvalue[k]+1) = 1;
k += 2;
}
erad = sqrt(erad2);
bonus->dvalue[k] = erad;
// the last element of bonus->dvalue is the rounded radius
rrad = 0.5 * dfile[j];
k++;
bonus->dvalue[k] = rrad;
atom->radius[bonus->ilocal] = erad + rrad;
} else { // polygons
nedges = nsub;
for (i = 0; i < nedges; i++) {
bonus->dvalue[k] = i;
m = i+1;
if (m == nedges) m = 0;
*(&bonus->dvalue[k]+1) = m;
k += 2;
}
// the next to last element is the enclosing radius
erad = sqrt(erad2);
bonus->dvalue[k] = erad;
// the last element of bonus->dvalue is the rounded radius
rrad = 0.5 * dfile[j];
k++;
bonus->dvalue[k] = rrad;
atom->radius[bonus->ilocal] = erad + rrad;
}
}
/* ----------------------------------------------------------------------
return radius of body particle defined by ifile/dfile params
params are ordered as in data file
called by Molecule class which needs single body size
------------------------------------------------------------------------- */
double BodyRoundedPolygon::radius_body(int ninteger, int ndouble,
int *ifile, double *dfile)
{
int nsub = ifile[0];
if (nsub < 1)
error->one(FLERR,"Incorrect integer value in "
"Bodies section of data file");
if (ndouble != 6 + 3*nsub + 1)
error->one(FLERR,"Incorrect # of floating-point values in "
"Bodies section of data file");
// sub-particle coords are relative to body center at (0,0,0)
// offset = 6 for sub-particle coords
double onerad;
double maxrad = 0.0;
double delta[3];
int offset = 6;
for (int i = 0; i < nsub; i++) {
delta[0] = dfile[offset];
delta[1] = dfile[offset+1];
delta[2] = dfile[offset+2];
offset += 3;
onerad = MathExtra::len3(delta);
maxrad = MAX(maxrad,onerad);
}
// add in radius of rounded corners
return maxrad + 0.5*dfile[offset];
}
/* ---------------------------------------------------------------------- */
int BodyRoundedPolygon::noutcol()
{
// the number of columns for the vertex coordinates
return 3;
}
/* ---------------------------------------------------------------------- */
int BodyRoundedPolygon::noutrow(int ibonus)
{
// only return the first nsub rows for the vertex coordinates
return avec->bonus[ibonus].ivalue[0];
}
/* ---------------------------------------------------------------------- */
void BodyRoundedPolygon::output(int ibonus, int m, double *values)
{
AtomVecBody::Bonus *bonus = &avec->bonus[ibonus];
double p[3][3];
MathExtra::quat_to_mat(bonus->quat,p);
MathExtra::matvec(p,&bonus->dvalue[3*m],values);
double *x = atom->x[bonus->ilocal];
values[0] += x[0];
values[1] += x[1];
values[2] += x[2];
}
/* ---------------------------------------------------------------------- */
int BodyRoundedPolygon::image(int ibonus, double flag1, double flag2,
int *&ivec, double **&darray)
{
int j;
double p[3][3];
double *x, rrad;
AtomVecBody::Bonus *bonus = &avec->bonus[ibonus];
int n = bonus->ivalue[0];
if (n == 1) {
for (int i = 0; i < n; i++) {
imflag[i] = SPHERE;
MathExtra::quat_to_mat(bonus->quat,p);
MathExtra::matvec(p,&bonus->dvalue[3*i],imdata[i]);
rrad = enclosing_radius(bonus);
x = atom->x[bonus->ilocal];
imdata[i][0] += x[0];
imdata[i][1] += x[1];
imdata[i][2] += x[2];
if (flag1 <= 0) imdata[i][3] = 2*rrad;
else imdata[i][3] = flag1;
}
} else {
// first end pt of each line
for (int i = 0; i < n; i++) {
imflag[i] = LINE;
MathExtra::quat_to_mat(bonus->quat,p);
MathExtra::matvec(p,&bonus->dvalue[3*i],imdata[i]);
rrad = rounded_radius(bonus);
x = atom->x[bonus->ilocal];
imdata[i][0] += x[0];
imdata[i][1] += x[1];
imdata[i][2] += x[2];
if (flag1 <= 0) imdata[i][6] = 2*rrad;
else imdata[i][6] = flag1;
}
// second end pt of each line
for (int i = 0; i < n; i++) {
j = i+1;
if (j == n) j = 0;
imdata[i][3] = imdata[j][0];
imdata[i][4] = imdata[j][1];
imdata[i][5] = imdata[j][2];
}
}
ivec = imflag;
darray = imdata;
return n;
}

View File

@ -0,0 +1,86 @@
/* -*- c++ -*- ----------------------------------------------------------
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 BODY_CLASS
BodyStyle(rounded/polygon,BodyRoundedPolygon)
#else
#ifndef LMP_BODY_ROUNDED_POLYGON_H
#define LMP_BODY_ROUNDED_POLYGON_H
#include "body.h"
#include "atom_vec_body.h"
namespace LAMMPS_NS {
class BodyRoundedPolygon : public Body {
public:
BodyRoundedPolygon(class LAMMPS *, int, char **);
~BodyRoundedPolygon();
int nsub(struct AtomVecBody::Bonus *);
double *coords(struct AtomVecBody::Bonus *);
int nedges(struct AtomVecBody::Bonus *);
double *edges(struct AtomVecBody::Bonus *);
double enclosing_radius(struct AtomVecBody::Bonus *);
double rounded_radius(struct AtomVecBody::Bonus *);
int pack_border_body(struct AtomVecBody::Bonus *, double *);
int unpack_border_body(struct AtomVecBody::Bonus *, double *);
void data_body(int, int, int, int *, double *);
double radius_body(int, int, int *, double *);
int noutrow(int);
int noutcol();
void output(int, int, double *);
int image(int, double, double, int *&, double **&);
private:
int *imflag;
double **imdata;
};
}
#endif
#endif
/* ERROR/WARNING messages:
E: Invalid body rounded/polygon command
Arguments in atom-style command are not correct.
E: Invalid format in Bodies section of data file
The specified number of integer or floating point values does not
appear.
E: Incorrect # of integer values in Bodies section of data file
See doc page for body style.
E: Incorrect integer value in Bodies section of data file
See doc page for body style.
E: Incorrect # of floating-point values in Bodies section of data file
See doc page for body style.
E: Insufficient Jacobi rotations for body nparticle
Eigensolve for rigid body was not sufficiently accurate.
*/

View File

@ -0,0 +1,526 @@
/* ----------------------------------------------------------------------
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.
------------------------------------------------------------------------- */
/* ----------------------------------------------------------------------
Contributing author: Trung Dac Nguyen (ndactrung@gmail.com)
------------------------------------------------------------------------- */
#include <cstdlib>
#include "body_rounded_polyhedron.h"
#include "atom_vec_body.h"
#include "atom.h"
#include "force.h"
#include "domain.h"
#include "math_extra.h"
#include "memory.h"
#include "error.h"
using namespace LAMMPS_NS;
#define EPSILON 1.0e-7
#define MAX_FACE_SIZE 4 // maximum number of vertices per face (for now)
enum{SPHERE,LINE}; // also in DumpImage
/* ---------------------------------------------------------------------- */
BodyRoundedPolyhedron::BodyRoundedPolyhedron(LAMMPS *lmp, int narg, char **arg) :
Body(lmp, narg, arg)
{
if (narg != 3) error->all(FLERR,"Invalid body rounded/polygon command");
// nmin and nmax are minimum and maximum number of vertices
int nmin = force->inumeric(FLERR,arg[1]);
int nmax = force->inumeric(FLERR,arg[2]);
if (nmin <= 0 || nmin > nmax)
error->all(FLERR,"Invalid body rounded/polyhedron command");
size_forward = 0;
// 3 integers: 1 for no. of vertices, 1 for no. of edges, 1 for no. of faces
// 3*nmax doubles for vertex coordinates + 2*nmax doubles for edge ends +
// (MAX_FACE_SIZE+1)*nmax for faces
// 1 double for the enclosing radius
// 1 double for the rounded radius
size_border = 3 + 3*nmax + 2*nmax + MAX_FACE_SIZE*nmax + 1 + 1;
// NOTE: need to set appropriate nnbin param for dcp
icp = new MyPoolChunk<int>(1,3);
dcp = new MyPoolChunk<double>(3*nmin+2+1+1,
3*nmax+2*nmax+MAX_FACE_SIZE*nmax+1+1);
memory->create(imflag,2*nmax,"body/rounded/polyhedron:imflag");
memory->create(imdata,2*nmax,7,"body/polyhedron:imdata");
}
/* ---------------------------------------------------------------------- */
BodyRoundedPolyhedron::~BodyRoundedPolyhedron()
{
delete icp;
delete dcp;
memory->destroy(imflag);
memory->destroy(imdata);
}
/* ---------------------------------------------------------------------- */
int BodyRoundedPolyhedron::nsub(AtomVecBody::Bonus *bonus)
{
return bonus->ivalue[0];
}
/* ---------------------------------------------------------------------- */
double *BodyRoundedPolyhedron::coords(AtomVecBody::Bonus *bonus)
{
return bonus->dvalue;
}
/* ---------------------------------------------------------------------- */
int BodyRoundedPolyhedron::nedges(AtomVecBody::Bonus *bonus)
{
int nvertices = bonus->ivalue[0];
int nedges = bonus->ivalue[1];
int nfaces = bonus->ivalue[2];
if (nvertices == 1) return 0;
else if (nvertices == 2) return 1;
return nedges; //(nvertices+nfaces-2); // Euler's polyon formula: V-E+F=2
}
/* ---------------------------------------------------------------------- */
double *BodyRoundedPolyhedron::edges(AtomVecBody::Bonus *bonus)
{
return bonus->dvalue+3*nsub(bonus);
}
/* ---------------------------------------------------------------------- */
int BodyRoundedPolyhedron::nfaces(AtomVecBody::Bonus *bonus)
{
return bonus->ivalue[2];
}
/* ---------------------------------------------------------------------- */
double *BodyRoundedPolyhedron::faces(AtomVecBody::Bonus *bonus)
{
int nvertices = bonus->ivalue[0];
if (nvertices == 1 || nvertices == 2) return NULL;
return bonus->dvalue+3*nsub(bonus)+2*nedges(bonus);
}
/* ---------------------------------------------------------------------- */
double BodyRoundedPolyhedron::enclosing_radius(struct AtomVecBody::Bonus *bonus)
{
int nvertices = bonus->ivalue[0];
if (nvertices == 1 || nvertices == 2)
return *(bonus->dvalue+3*nsub(bonus)+2);
return *(bonus->dvalue+3*nsub(bonus) + 2*nedges(bonus) +
MAX_FACE_SIZE*nfaces(bonus));
}
/* ---------------------------------------------------------------------- */
double BodyRoundedPolyhedron::rounded_radius(struct AtomVecBody::Bonus *bonus)
{
int nvertices = bonus->ivalue[0];
if (nvertices == 1 || nvertices == 2)
return *(bonus->dvalue+3*nsub(bonus)+2+1);
return *(bonus->dvalue+3*nsub(bonus) + 2*nedges(bonus) +
MAX_FACE_SIZE*nfaces(bonus)+1);
}
/* ---------------------------------------------------------------------- */
int BodyRoundedPolyhedron::pack_border_body(AtomVecBody::Bonus *bonus, double *buf)
{
int nsub = bonus->ivalue[0];
int ned = bonus->ivalue[1];
int nfac = bonus->ivalue[2];
buf[0] = nsub;
buf[1] = ned;
buf[2] = nfac;
int ndouble;
if (nsub == 1 || nsub == 2) ndouble = 3*nsub+2+MAX_FACE_SIZE*nfac+1+1;
else ndouble = 3*nsub+2*ned+MAX_FACE_SIZE*nfac+1+1;
memcpy(&buf[3],bonus->dvalue,ndouble*sizeof(double));
return 3+ndouble;
}
/* ---------------------------------------------------------------------- */
int BodyRoundedPolyhedron::unpack_border_body(AtomVecBody::Bonus *bonus,
double *buf)
{
int nsub = static_cast<int> (buf[0]);
int ned = static_cast<int> (buf[1]);
int nfac = static_cast<int> (buf[2]);
bonus->ivalue[0] = nsub;
bonus->ivalue[1] = ned;
bonus->ivalue[2] = nfac;
int ndouble;
if (nsub == 1 || nsub == 2) ndouble = 3*nsub+2+MAX_FACE_SIZE*nfac+1+1;
else ndouble = 3*nsub+2*ned+MAX_FACE_SIZE*nfac+1+1;
memcpy(bonus->dvalue,&buf[3],ndouble*sizeof(double));
return 3+ndouble;
}
/* ----------------------------------------------------------------------
populate bonus data structure with data file values
------------------------------------------------------------------------- */
void BodyRoundedPolyhedron::data_body(int ibonus, int ninteger, int ndouble,
int *ifile, double *dfile)
{
AtomVecBody::Bonus *bonus = &avec->bonus[ibonus];
// set ninteger, ndouble in bonus and allocate 2 vectors of ints, doubles
if (ninteger != 3)
error->one(FLERR,"Incorrect # of integer values in "
"Bodies section of data file");
int nsub = ifile[0];
int ned = ifile[1];
int nfac = ifile[2];
if (nsub < 1)
error->one(FLERR,"Incorrect integer value in "
"Bodies section of data file");
// nentries = number of double entries to be read from Body section:
// nsub == 1 || nsub == 2 || nsub == 3:
// 6 for inertia + 3*nsub for vertex coords + 1 for rounded radius
// nsub > 3:
// 6 for inertia + 3*nsub for vertex coords + 2*nsub for edges +
// 3*nfaces + 1 for rounded radius
int nedges,nentries;
if (nsub == 1 || nsub == 2) {
nentries = 6 + 3*nsub + 1;
} else {
nedges = ned; //nsub + nfac - 2;
nentries = 6 + 3*nsub + 2*nedges + MAX_FACE_SIZE*nfac + 1;
}
if (ndouble != nentries)
error->one(FLERR,"Incorrect # of floating-point values in "
"Bodies section of data file");
bonus->ninteger = 3;
bonus->ivalue = icp->get(bonus->iindex);
bonus->ivalue[0] = nsub;
bonus->ivalue[1] = ned;
bonus->ivalue[2] = nfac;
if (nsub == 1 || nsub == 2) bonus->ndouble = 3*nsub + 2*nsub + 1 + 1;
else bonus->ndouble = 3*nsub + 2*nedges + MAX_FACE_SIZE*nfac + 1 + 1;
bonus->dvalue = dcp->get(bonus->ndouble,bonus->dindex);
// diagonalize inertia tensor
double tensor[3][3];
tensor[0][0] = dfile[0];
tensor[1][1] = dfile[1];
tensor[2][2] = dfile[2];
tensor[0][1] = tensor[1][0] = dfile[3];
tensor[0][2] = tensor[2][0] = dfile[4];
tensor[1][2] = tensor[2][1] = dfile[5];
double *inertia = bonus->inertia;
double evectors[3][3];
int ierror = MathExtra::jacobi(tensor,inertia,evectors);
if (ierror) error->one(FLERR,
"Insufficient Jacobi rotations for body nparticle");
// if any principal moment < scaled EPSILON, set to 0.0
double max;
max = MAX(inertia[0],inertia[1]);
max = MAX(max,inertia[2]);
if (inertia[0] < EPSILON*max) inertia[0] = 0.0;
if (inertia[1] < EPSILON*max) inertia[1] = 0.0;
if (inertia[2] < EPSILON*max) inertia[2] = 0.0;
// exyz_space = principal axes in space frame
double ex_space[3],ey_space[3],ez_space[3];
ex_space[0] = evectors[0][0];
ex_space[1] = evectors[1][0];
ex_space[2] = evectors[2][0];
ey_space[0] = evectors[0][1];
ey_space[1] = evectors[1][1];
ey_space[2] = evectors[2][1];
ez_space[0] = evectors[0][2];
ez_space[1] = evectors[1][2];
ez_space[2] = evectors[2][2];
// enforce 3 evectors as a right-handed coordinate system
// flip 3rd vector if needed
double cross[3];
MathExtra::cross3(ex_space,ey_space,cross);
if (MathExtra::dot3(cross,ez_space) < 0.0) MathExtra::negate3(ez_space);
// create initial quaternion
MathExtra::exyz_to_q(ex_space,ey_space,ez_space,bonus->quat);
// bonus->dvalue = the first 3*nsub elements are sub-particle displacements
// find the enclosing radius of the body from the maximum displacement
int i,m;
double delta[3], rsq, erad, rrad;
double erad2 = 0;
int j = 6;
int k = 0;
for (i = 0; i < nsub; i++) {
delta[0] = dfile[j];
delta[1] = dfile[j+1];
delta[2] = dfile[j+2];
MathExtra::transpose_matvec(ex_space,ey_space,ez_space,
delta,&bonus->dvalue[k]);
rsq = delta[0] * delta[0] + delta[1] * delta[1] +
delta[2] * delta[2];
if (rsq > erad2) erad2 = rsq;
j += 3;
k += 3;
}
// .. the next 2*nsub elements are edge ends
if (nsub == 1) { // spheres
nedges = 0;
bonus->dvalue[k] = 0;
*(&bonus->dvalue[k]+1) = 0;
k += 2;
rrad = 0.5 * dfile[j];
bonus->dvalue[k] = rrad;
erad = rrad; // enclosing radius = rounded_radius
// the last element of bonus->dvalue is the rounded radius
k++;
bonus->dvalue[k] = rrad;
atom->radius[bonus->ilocal] = erad;
} else if (nsub == 2) { // rods
nedges = 1;
for (i = 0; i < nedges; i++) {
bonus->dvalue[k] = 0;
*(&bonus->dvalue[k]+1) = 1;
k += 2;
}
erad = sqrt(erad2);
bonus->dvalue[k] = erad;
// the last element of bonus->dvalue is the rounded radius
rrad = 0.5 * dfile[j];
k++;
bonus->dvalue[k] = rrad;
atom->radius[bonus->ilocal] = erad + rrad;
} else { // polyhedra
// edges
for (i = 0; i < nedges; i++) {
bonus->dvalue[k] = dfile[j];
*(&bonus->dvalue[k]+1) = dfile[j+1];
k += 2;
j += 2;
}
// faces
for (i = 0; i < nfac; i++) {
for (m = 0; m < MAX_FACE_SIZE; m++)
*(&bonus->dvalue[k]+m) = dfile[j+m];
k += MAX_FACE_SIZE;
j += MAX_FACE_SIZE;
}
// the next to last element is the enclosing radius
erad = sqrt(erad2);
bonus->dvalue[k] = erad;
// the last element bonus-> dvalue is the rounded radius
rrad = 0.5 * dfile[j];
k++;
bonus->dvalue[k] = rrad;
atom->radius[bonus->ilocal] = erad + rrad;
}
}
/* ----------------------------------------------------------------------
return radius of body particle defined by ifile/dfile params
params are ordered as in data file
called by Molecule class which needs single body size
------------------------------------------------------------------------- */
double BodyRoundedPolyhedron::radius_body(int ninteger, int ndouble,
int *ifile, double *dfile)
{
int nsub = ifile[0];
int ned = ifile[1];
int nfac = ifile[2];
int nedges = ned; //nsub + nfac - 2;
int nentries;
if (nsub == 1 || nsub == 2) nentries = 6 + 3*nsub + 1;
else nentries = 6 + 3*nsub + 2*nedges + MAX_FACE_SIZE*nfac + 1;
if (nsub < 1)
error->one(FLERR,"Incorrect integer value in "
"Bodies section of data file");
if (ndouble != nentries)
error->one(FLERR,"Incorrect # of floating-point values in "
"Bodies section of data file");
// sub-particle coords are relative to body center at (0,0,0)
// offset = 6 for sub-particle coords
double onerad;
double maxrad = 0.0;
double delta[3];
int offset = 6;
for (int i = 0; i < nsub; i++) {
delta[0] = dfile[offset];
delta[1] = dfile[offset+1];
delta[2] = dfile[offset+2];
offset += 3;
onerad = MathExtra::len3(delta);
maxrad = MAX(maxrad,onerad);
}
if (nsub > 2) offset += (2*nedges+MAX_FACE_SIZE*nfac);
// add in radius of rounded corners
return maxrad + 0.5*dfile[offset];
}
/* ---------------------------------------------------------------------- */
int BodyRoundedPolyhedron::noutcol()
{
// the number of columns for the vertex coordinates
return 3;
}
/* ---------------------------------------------------------------------- */
int BodyRoundedPolyhedron::noutrow(int ibonus)
{
// only return the first nsub rows for the vertex coordinates
return avec->bonus[ibonus].ivalue[0];
}
/* ---------------------------------------------------------------------- */
void BodyRoundedPolyhedron::output(int ibonus, int m, double *values)
{
AtomVecBody::Bonus *bonus = &avec->bonus[ibonus];
double p[3][3];
MathExtra::quat_to_mat(bonus->quat,p);
MathExtra::matvec(p,&bonus->dvalue[3*m],values);
double *x = atom->x[bonus->ilocal];
values[0] += x[0];
values[1] += x[1];
values[2] += x[2];
}
/* ---------------------------------------------------------------------- */
int BodyRoundedPolyhedron::image(int ibonus, double flag1, double flag2,
int *&ivec, double **&darray)
{
int j, nelements;
double p[3][3];
double *x, rrad;
AtomVecBody::Bonus *bonus = &avec->bonus[ibonus];
int nvertices = bonus->ivalue[0];
if (nvertices == 1) { // spheres
for (int i = 0; i < nvertices; i++) {
imflag[i] = SPHERE;
MathExtra::quat_to_mat(bonus->quat,p);
MathExtra::matvec(p,&bonus->dvalue[3*i],imdata[i]);
rrad = enclosing_radius(bonus);
x = atom->x[bonus->ilocal];
imdata[i][0] += x[0];
imdata[i][1] += x[1];
imdata[i][2] += x[2];
if (flag1 <= 0) imdata[i][3] = 2*rrad;
else imdata[i][3] = flag1;
}
nelements = nvertices;
} else {
int nfaces = bonus->ivalue[2];
int nedges = bonus->ivalue[1]; //nvertices + nfaces - 2;
if (nvertices == 2) nedges = 1; // special case: rods
double* edge_ends = &bonus->dvalue[3*nvertices];
int pt1, pt2;
for (int i = 0; i < nedges; i++) {
imflag[i] = LINE;
pt1 = static_cast<int>(edge_ends[2*i]);
pt2 = static_cast<int>(edge_ends[2*i+1]);
MathExtra::quat_to_mat(bonus->quat,p);
MathExtra::matvec(p,&bonus->dvalue[3*pt1],imdata[i]);
MathExtra::matvec(p,&bonus->dvalue[3*pt2],&imdata[i][3]);
rrad = rounded_radius(bonus);
x = atom->x[bonus->ilocal];
imdata[i][0] += x[0];
imdata[i][1] += x[1];
imdata[i][2] += x[2];
imdata[i][3] += x[0];
imdata[i][4] += x[1];
imdata[i][5] += x[2];
if (flag1 <= 0) imdata[i][6] = 2*rrad;
else imdata[i][6] = flag1;
}
nelements = nedges;
}
ivec = imflag;
darray = imdata;
return nelements;
}

View File

@ -0,0 +1,88 @@
/* -*- c++ -*- ----------------------------------------------------------
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 BODY_CLASS
BodyStyle(rounded/polyhedron,BodyRoundedPolyhedron)
#else
#ifndef LMP_BODY_ROUNDED_POLYHEDRON_H
#define LMP_BODY_ROUNDED_POLYHEDRON_H
#include "body.h"
#include "atom_vec_body.h"
namespace LAMMPS_NS {
class BodyRoundedPolyhedron : public Body {
public:
BodyRoundedPolyhedron(class LAMMPS *, int, char **);
~BodyRoundedPolyhedron();
int nsub(struct AtomVecBody::Bonus *);
double *coords(struct AtomVecBody::Bonus *);
int nedges(struct AtomVecBody::Bonus *);
double *edges(struct AtomVecBody::Bonus *);
int nfaces(struct AtomVecBody::Bonus *);
double *faces(struct AtomVecBody::Bonus *);
double enclosing_radius(struct AtomVecBody::Bonus *);
double rounded_radius(struct AtomVecBody::Bonus *);
int pack_border_body(struct AtomVecBody::Bonus *, double *);
int unpack_border_body(struct AtomVecBody::Bonus *, double *);
void data_body(int, int, int, int *, double *);
double radius_body(int, int, int *, double *);
int noutrow(int);
int noutcol();
void output(int, int, double *);
int image(int, double, double, int *&, double **&);
private:
int *imflag;
double **imdata;
};
}
#endif
#endif
/* ERROR/WARNING messages:
E: Invalid body rounded/polyhedron command
Arguments in atom-style command are not correct.
E: Invalid format in Bodies section of data file
The specified number of integer or floating point values does not
appear.
E: Incorrect # of integer values in Bodies section of data file
See doc page for body style.
E: Incorrect integer value in Bodies section of data file
See doc page for body style.
E: Incorrect # of floating-point values in Bodies section of data file
See doc page for body style.
E: Insufficient Jacobi rotations for body rounded/polyhedron
Eigensolve for rigid body was not sufficiently accurate.
*/

View File

@ -0,0 +1,832 @@
/* ----------------------------------------------------------------------
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.
------------------------------------------------------------------------- */
/* ----------------------------------------------------------------------
Contributing author: Trung Dac Nguyen (ndactrung@gmail.com)
------------------------------------------------------------------------- */
#include <cmath>
#include <cstdlib>
#include <cstring>
#include "fix_wall_body_polygon.h"
#include "atom.h"
#include "atom_vec_body.h"
#include "body_rounded_polygon.h"
#include "domain.h"
#include "update.h"
#include "force.h"
#include "pair.h"
#include "modify.h"
#include "respa.h"
#include "math_const.h"
#include "math_extra.h"
#include "memory.h"
#include "error.h"
using namespace LAMMPS_NS;
using namespace FixConst;
using namespace MathConst;
enum{XPLANE=0,YPLANE=1,ZCYLINDER}; // XYZ PLANE need to be 0,1,2
enum{HOOKE,HOOKE_HISTORY};
enum {INVALID=0,NONE=1,VERTEX=2};
enum {FAR=0,XLO,XHI,YLO,YHI};
//#define _POLYGON_DEBUG
#define DELTA 10000
#define EPSILON 1e-2
#define BIG 1.0e20
#define MAX_CONTACTS 4 // maximum number of contacts for 2D models
#define EFF_CONTACTS 2 // effective contacts for 2D models
/* ---------------------------------------------------------------------- */
FixWallBodyPolygon::FixWallBodyPolygon(LAMMPS *lmp, int narg, char **arg) :
Fix(lmp, narg, arg)
{
if (narg < 7) error->all(FLERR,"Illegal fix wall/body/polygon command");
if (!atom->body_flag)
error->all(FLERR,"Fix wall/body/polygon requires "
"atom style body/rounded/polygon");
restart_peratom = 1;
create_attribute = 1;
// wall/particle coefficients
kn = force->numeric(FLERR,arg[3]);
c_n = force->numeric(FLERR,arg[4]);
if (strcmp(arg[5],"NULL") == 0) c_t = 0.5 * c_n;
else c_t = force->numeric(FLERR,arg[5]);
if (kn < 0.0 || c_n < 0.0 || c_t < 0.0)
error->all(FLERR,"Illegal fix wall/body/polygon command");
// wallstyle args
int iarg = 6;
if (strcmp(arg[iarg],"xplane") == 0) {
if (narg < iarg+3) error->all(FLERR,"Illegal fix wall/body/polygon command");
wallstyle = XPLANE;
if (strcmp(arg[iarg+1],"NULL") == 0) lo = -BIG;
else lo = force->numeric(FLERR,arg[iarg+1]);
if (strcmp(arg[iarg+2],"NULL") == 0) hi = BIG;
else hi = force->numeric(FLERR,arg[iarg+2]);
iarg += 3;
} else if (strcmp(arg[iarg],"yplane") == 0) {
if (narg < iarg+3) error->all(FLERR,"Illegal fix wall/body/polygon command");
wallstyle = YPLANE;
if (strcmp(arg[iarg+1],"NULL") == 0) lo = -BIG;
else lo = force->numeric(FLERR,arg[iarg+1]);
if (strcmp(arg[iarg+2],"NULL") == 0) hi = BIG;
else hi = force->numeric(FLERR,arg[iarg+2]);
iarg += 3;
} else if (strcmp(arg[iarg],"zcylinder") == 0) {
if (narg < iarg+2) error->all(FLERR,"Illegal fix wall/body/polygon command");
wallstyle = ZCYLINDER;
lo = hi = 0.0;
cylradius = force->numeric(FLERR,arg[iarg+1]);
iarg += 2;
}
// check for trailing keyword/values
wiggle = 0;
while (iarg < narg) {
if (strcmp(arg[iarg],"wiggle") == 0) {
if (iarg+4 > narg) error->all(FLERR,"Illegal fix wall/body/polygon command");
if (strcmp(arg[iarg+1],"x") == 0) axis = 0;
else if (strcmp(arg[iarg+1],"y") == 0) axis = 1;
else if (strcmp(arg[iarg+1],"z") == 0) axis = 2;
else error->all(FLERR,"Illegal fix wall/body/polygon command");
amplitude = force->numeric(FLERR,arg[iarg+2]);
period = force->numeric(FLERR,arg[iarg+3]);
wiggle = 1;
iarg += 4;
} else error->all(FLERR,"Illegal fix wall/body/polygon command");
}
if (wallstyle == XPLANE && domain->xperiodic)
error->all(FLERR,"Cannot use wall in periodic dimension");
if (wallstyle == YPLANE && domain->yperiodic)
error->all(FLERR,"Cannot use wall in periodic dimension");
if (wallstyle == ZCYLINDER && (domain->xperiodic || domain->yperiodic))
error->all(FLERR,"Cannot use wall in periodic dimension");
if (wiggle && wallstyle == ZCYLINDER && axis != 2)
error->all(FLERR,"Invalid wiggle direction for fix wall/body/polygon");
// setup oscillations
if (wiggle) omega = 2.0*MY_PI / period;
time_origin = update->ntimestep;
dmax = nmax = 0;
discrete = NULL;
dnum = dfirst = NULL;
edmax = ednummax = 0;
edge = NULL;
ednum = edfirst = NULL;
enclosing_radius = NULL;
rounded_radius = NULL;
}
/* ---------------------------------------------------------------------- */
FixWallBodyPolygon::~FixWallBodyPolygon()
{
memory->destroy(discrete);
memory->destroy(dnum);
memory->destroy(dfirst);
memory->destroy(edge);
memory->destroy(ednum);
memory->destroy(edfirst);
memory->destroy(enclosing_radius);
memory->destroy(rounded_radius);
}
/* ---------------------------------------------------------------------- */
int FixWallBodyPolygon::setmask()
{
int mask = 0;
mask |= POST_FORCE;
return mask;
}
/* ---------------------------------------------------------------------- */
void FixWallBodyPolygon::init()
{
dt = update->dt;
avec = (AtomVecBody *) atom->style_match("body");
if (!avec)
error->all(FLERR,"Pair body/rounded/polygon requires atom style body");
if (strcmp(avec->bptr->style,"rounded/polygon") != 0)
error->all(FLERR,"Pair body/rounded/polygon requires "
"body style rounded/polygon");
bptr = (BodyRoundedPolygon *) avec->bptr;
// set pairstyle from body/polygonular pair style
if (force->pair_match("body/rounded/polygon",1))
pairstyle = HOOKE;
else error->all(FLERR,"Fix wall/body/polygon is incompatible with Pair style");
}
/* ---------------------------------------------------------------------- */
void FixWallBodyPolygon::setup(int vflag)
{
if (strstr(update->integrate_style,"verlet"))
post_force(vflag);
}
/* ---------------------------------------------------------------------- */
void FixWallBodyPolygon::post_force(int vflag)
{
double vwall[3],dx,dy,dz,del1,del2,delxy,delr,rsq,eradi,rradi,wall_pos;
int i,ni,npi,ifirst,nei,iefirst,side;
double facc[3];
// set position of wall to initial settings and velocity to 0.0
// if wiggle, set wall position and velocity accordingly
double wlo = lo;
double whi = hi;
vwall[0] = vwall[1] = vwall[2] = 0.0;
if (wiggle) {
double arg = omega * (update->ntimestep - time_origin) * dt;
if (wallstyle == axis) {
wlo = lo + amplitude - amplitude*cos(arg);
whi = hi + amplitude - amplitude*cos(arg);
}
vwall[axis] = amplitude*omega*sin(arg);
}
// loop over all my atoms
// rsq = distance from wall
// dx,dy,dz = signed distance from wall
// for rotating cylinder, reset vwall based on particle position
// skip atom if not close enough to wall
// if wall was set to NULL, it's skipped since lo/hi are infinity
// compute force and torque on atom if close enough to wall
// via wall potential matched to pair potential
double **x = atom->x;
double **v = atom->v;
double **f = atom->f;
int *body = atom->body;
double *radius = atom->radius;
double **torque = atom->torque;
double **angmom = atom->angmom;
int *mask = atom->mask;
int nlocal = atom->nlocal;
// grow the per-atom lists if necessary and initialize
if (atom->nmax > nmax) {
memory->destroy(dnum);
memory->destroy(dfirst);
memory->destroy(ednum);
memory->destroy(edfirst);
memory->destroy(enclosing_radius);
memory->destroy(rounded_radius);
nmax = atom->nmax;
memory->create(dnum,nmax,"fix:dnum");
memory->create(dfirst,nmax,"fix:dfirst");
memory->create(ednum,nmax,"fix:ednum");
memory->create(edfirst,nmax,"fix:edfirst");
memory->create(enclosing_radius,nmax,"fix:enclosing_radius");
memory->create(rounded_radius,nmax,"fix:rounded_radius");
}
ndiscrete = nedge = 0;
for (i = 0; i < nlocal; i++)
dnum[i] = ednum[i] = 0;
for (i = 0; i < nlocal; i++) {
if (mask[i] & groupbit) {
if (body[i] < 0) continue;
dx = dy = dz = 0.0;
side = FAR;
if (wallstyle == XPLANE) {
del1 = x[i][0] - wlo;
del2 = whi - x[i][0];
if (del1 < del2) {
dx = del1;
wall_pos = wlo;
side = XLO;
} else {
dx = -del2;
wall_pos = whi;
side = XHI;
}
} else if (wallstyle == YPLANE) {
del1 = x[i][1] - wlo;
del2 = whi - x[i][1];
if (del1 < del2) {
dy = del1;
wall_pos = wlo;
side = YLO;
} else {
dy = -del2;
wall_pos = whi;
side = YHI;
}
} else if (wallstyle == ZCYLINDER) {
delxy = sqrt(x[i][0]*x[i][0] + x[i][1]*x[i][1]);
delr = cylradius - delxy;
if (delr > eradi) dz = cylradius;
else {
dx = -delr/delxy * x[i][0];
dy = -delr/delxy * x[i][1];
}
}
rsq = dx*dx + dy*dy + dz*dz;
if (rsq > radius[i]*radius[i]) continue;
double r = sqrt(rsq);
double rsqinv = 1.0 / rsq;
if (dnum[i] == 0) body2space(i);
npi = dnum[i];
ifirst = dfirst[i];
nei = ednum[i];
iefirst = edfirst[i];
eradi = enclosing_radius[i];
rradi = rounded_radius[i];
// reset vertex and edge forces
for (ni = 0; ni < npi; ni++) {
discrete[ifirst+ni][3] = 0;
discrete[ifirst+ni][4] = 0;
discrete[ifirst+ni][5] = 0;
}
for (ni = 0; ni < nei; ni++) {
edge[iefirst+ni][2] = 0;
edge[iefirst+ni][3] = 0;
edge[iefirst+ni][4] = 0;
}
int interact, num_contacts, done;
double delta_a, delta_ua, j_a;
Contact contact_list[MAX_CONTACTS];
num_contacts = 0;
facc[0] = facc[1] = facc[2] = 0;
interact = vertex_against_wall(i, wall_pos, x, f, torque, side,
contact_list, num_contacts, facc);
if (num_contacts >= 2) {
// find the first two distinct contacts
done = 0;
for (int m = 0; m < num_contacts-1; m++) {
for (int n = m+1; n < num_contacts; n++) {
delta_a = contact_separation(contact_list[m], contact_list[n]);
if (delta_a > 0) {
delta_ua = 1.0;
j_a = delta_a / (EFF_CONTACTS * delta_ua);
if (j_a < 1.0) j_a = 1.0;
// scale the force at both contacts
contact_forces(contact_list[m], j_a, x, v, angmom, f, torque,
vwall, facc);
contact_forces(contact_list[n], j_a, x, v, angmom, f, torque,
vwall, facc);
done = 1;
break;
}
}
if (done == 1) break;
}
} else if (num_contacts == 1) {
// if there's only one contact, it should be handled here
// since forces/torques have not been accumulated from vertex2wall()
contact_forces(contact_list[0], 1.0, x, v, angmom, f, torque,
vwall, facc);
}
} // group bit
}
}
/* ---------------------------------------------------------------------- */
void FixWallBodyPolygon::reset_dt()
{
dt = update->dt;
}
/* ----------------------------------------------------------------------
convert N sub-particles in body I to space frame using current quaternion
store sub-particle space-frame displacements from COM in discrete list
------------------------------------------------------------------------- */
void FixWallBodyPolygon::body2space(int i)
{
int ibonus = atom->body[i];
AtomVecBody::Bonus *bonus = &avec->bonus[ibonus];
int nsub = bptr->nsub(bonus);
double *coords = bptr->coords(bonus);
int body_num_edges = bptr->nedges(bonus);
double* vertices = bptr->edges(bonus);
double eradius = bptr->enclosing_radius(bonus);
double rradius = bptr->rounded_radius(bonus);
// get the number of sub-particles (vertices)
// and the index of the first vertex of my body in the list
dnum[i] = nsub;
dfirst[i] = ndiscrete;
// grow the vertex list if necessary
// the first 3 columns are for coords, the last 3 for forces
if (ndiscrete + nsub > dmax) {
dmax += DELTA;
memory->grow(discrete,dmax,6,"fix:discrete");
}
double p[3][3];
MathExtra::quat_to_mat(bonus->quat,p);
for (int m = 0; m < nsub; m++) {
MathExtra::matvec(p,&coords[3*m],discrete[ndiscrete]);
discrete[ndiscrete][3] = 0;
discrete[ndiscrete][4] = 0;
discrete[ndiscrete][5] = 0;
ndiscrete++;
}
// get the number of edges (vertices)
// and the index of the first edge of my body in the list
ednum[i] = body_num_edges;
edfirst[i] = nedge;
// grow the edge list if necessary
// the first 2 columns are for vertex indices within body,
// the last 3 for forces
if (nedge + body_num_edges > edmax) {
edmax += DELTA;
memory->grow(edge,edmax,5,"fix:edge");
}
for (int m = 0; m < body_num_edges; m++) {
edge[nedge][0] = static_cast<int>(vertices[2*m+0]);
edge[nedge][1] = static_cast<int>(vertices[2*m+1]);
edge[nedge][2] = 0;
edge[nedge][3] = 0;
edge[nedge][4] = 0;
nedge++;
}
enclosing_radius[i] = eradius;
rounded_radius[i] = rradius;
}
/* ----------------------------------------------------------------------
Determine the interaction mode between i's vertices against the wall
i = atom i (body i)
x = atoms' coordinates
f = atoms' forces
torque = atoms' torques
Return:
contact_list = list of contacts between i and the wall
num_contacts = number of contacts between i's vertices and the wall
interact = 0 no interaction with the wall
1 there's at least one vertex of i interacts
with the wall
---------------------------------------------------------------------- */
int FixWallBodyPolygon::vertex_against_wall(int i, double wall_pos,
double** x, double** f, double** torque, int side,
Contact* contact_list, int &num_contacts, double* facc)
{
int ni, npi, ifirst, interact;
double xpi[3], xpj[3], dist, eradi, rradi;
double fx, fy, fz, rx, ry, rz;
int nlocal = atom->nlocal;
npi = dnum[i];
ifirst = dfirst[i];
eradi = enclosing_radius[i];
rradi = rounded_radius[i];
interact = 0;
// loop through body i's vertices
for (ni = 0; ni < npi; ni++) {
// convert body-fixed coordinates to space-fixed, xi
xpi[0] = x[i][0] + discrete[ifirst+ni][0];
xpi[1] = x[i][1] + discrete[ifirst+ni][1];
xpi[2] = x[i][2] + discrete[ifirst+ni][2];
int mode, contact, p2vertex;
double d, R, hi[3], t, delx, dely, delz, fpair, shift;
double xj[3], rij;
// compute the distance from the vertex xpi to the wall
mode = compute_distance_to_wall(xpi, rradi, wall_pos, side,
d, hi, contact);
if (mode == INVALID || mode == NONE) continue;
if (mode == VERTEX) {
interact = 1;
// vertex i interacts with the wall
delx = xpi[0] - hi[0];
dely = xpi[1] - hi[1];
delz = xpi[2] - hi[2];
// R = surface separation = d shifted by the rounded radius
// R = d - p1.rounded_radius;
// note: the force is defined for R, not for d
// R > 0: no interaction
// R <= 0: deformation between vertex i and the wall
rij = sqrt(delx*delx + dely*dely + delz*delz);
R = rij - rradi;
// the normal frictional term -c_n * vn will be added later
if (R <= 0) { // deformation occurs
fpair = -kn * R;
} else fpair = 0.0;
fx = delx*fpair/rij;
fy = dely*fpair/rij;
fz = delz*fpair/rij;
#ifdef _POLYGON_DEBUG
printf(" Interaction between vertex %d of %d and wall:", ni);
printf(" mode = %d; contact = %d; d = %f; rij = %f\n",
mode, contact, d, rij);
printf(" R = %f\n", R);
printf(" fpair = %f\n", fpair);
#endif
if (contact == 1) {
// vertex ni of body i contacts with edge nj of body j
contact_list[num_contacts].ibody = i;
contact_list[num_contacts].jbody = -1;
contact_list[num_contacts].vertex = ni;
contact_list[num_contacts].edge = -1;
contact_list[num_contacts].xv[0] = xpi[0];
contact_list[num_contacts].xv[1] = xpi[1];
contact_list[num_contacts].xv[2] = xpi[2];
contact_list[num_contacts].xe[0] = hi[0];
contact_list[num_contacts].xe[1] = hi[1];
contact_list[num_contacts].xe[2] = hi[2];
contact_list[num_contacts].separation = R;
num_contacts++;
// store forces to vertex ni to be rescaled later,
// if there are 2 contacts
discrete[ifirst+ni][3] = fx;
discrete[ifirst+ni][4] = fy;
discrete[ifirst+ni][5] = fz;
#ifdef _POLYGON_DEBUG
printf(" Stored forces at vertex and edge for accumulating later.\n");
#endif
} else { // no contact
// accumulate force and torque to the body directly
f[i][0] += fx;
f[i][1] += fy;
f[i][2] += fz;
sum_torque(x[i], xpi, fx, fy, fz, torque[i]);
} // end if contact
} // end if mode
} // end for looping through the vertices of body i
return interact;
}
/* -------------------------------------------------------------------------
Compute the distance between a vertex to the wall
another body
Input:
x0 = coordinate of the tested vertex
rradi = rounded radius of the vertex
wall_pos = position of the wall
Output:
d = Distance from a point x0 to an wall
hi = coordinates of the projection of x0 on the wall
contact = 0 no contact between the queried vertex and the wall
1 contact detected
return NONE if there is no interaction
EDGE if the tested vertex interacts with the wall
------------------------------------------------------------------------- */
int FixWallBodyPolygon::compute_distance_to_wall(double* x0, double rradi,
double wall_pos, int side, double &d, double hi[3], int &contact)
{
int mode;
double delxy;
// h0 = position of the projection of x0 on the wall
if (wallstyle == XPLANE) {
hi[0] = wall_pos;
hi[1] = x0[1];
hi[2] = x0[2];
} else if (wallstyle == YPLANE) {
hi[0] = x0[0];
hi[1] = wall_pos;
hi[2] = x0[2];
} else if (wallstyle == ZCYLINDER) {
delxy = sqrt(x0[0]*x0[0] + x0[1]*x0[1]);
hi[0] = x0[0]*cylradius/delxy;
hi[1] = x0[1]*cylradius/delxy;
hi[2] = x0[2];
}
// distance from x0 to the wall = distance from x0 to hi
distance(hi, x0, d);
// determine the interaction mode
if (d < rradi) {
mode = VERTEX;
contact = 1;
} else {
if (side == XLO) {
if (x0[0] < wall_pos) mode = VERTEX;
else mode = NONE;
} else if (side == XHI) {
if (x0[0] > wall_pos) mode = VERTEX;
else mode = NONE;
} else if (side == YLO) {
if (x0[1] < wall_pos) mode = VERTEX;
else mode = NONE;
} else if (side == YHI) {
if (x0[1] > wall_pos) mode = VERTEX;
else mode = NONE;
}
}
if (mode == NONE) contact = 0;
else contact = 1;
return mode;
}
/* ----------------------------------------------------------------------
Compute the contact forces between two bodies
modify the force stored at the vertex and edge in contact by j_a
sum forces and torque to the corresponding bodies
fn = normal friction component
ft = tangential friction component (-c_t * vrt)
------------------------------------------------------------------------- */
void FixWallBodyPolygon::contact_forces(Contact& contact, double j_a,
double** x, double** v, double** angmom, double** f,
double** torque, double* vwall, double* facc)
{
int ibody,ibonus,ifirst, jefirst, ni;
double fx,fy,fz,delx,dely,delz,rsq,rsqinv;
double vr1,vr2,vr3,vnnr,vn1,vn2,vn3,vt1,vt2,vt3;
double fn[3],ft[3],vi[3];
double *quat, *inertia;
AtomVecBody::Bonus *bonus;
ibody = contact.ibody;
// compute the velocity of the vertex in the space-fixed frame
ibonus = atom->body[ibody];
bonus = &avec->bonus[ibonus];
quat = bonus->quat;
inertia = bonus->inertia;
total_velocity(contact.xv, x[ibody], v[ibody], angmom[ibody],
inertia, quat, vi);
// vector pointing from the vertex to the point on the wall
delx = contact.xv[0] - contact.xe[0];
dely = contact.xv[1] - contact.xe[1];
delz = contact.xv[2] - contact.xe[2];
rsq = delx*delx + dely*dely + delz*delz;
rsqinv = 1.0/rsq;
// relative translational velocity
vr1 = vi[0] - vwall[0];
vr2 = vi[1] - vwall[1];
vr3 = vi[2] - vwall[2];
// normal component
vnnr = vr1*delx + vr2*dely + vr3*delz;
vn1 = delx*vnnr * rsqinv;
vn2 = dely*vnnr * rsqinv;
vn3 = delz*vnnr * rsqinv;
// tangential component
vt1 = vr1 - vn1;
vt2 = vr2 - vn2;
vt3 = vr3 - vn3;
// normal friction term at contact
fn[0] = -c_n * vn1;
fn[1] = -c_n * vn2;
fn[2] = -c_n * vn3;
// tangential friction term at contact
// excluding the tangential deformation term for now
ft[0] = -c_t * vt1;
ft[1] = -c_t * vt2;
ft[2] = -c_t * vt3;
// only the cohesive force is scaled by j_a
ifirst = dfirst[ibody];
ni = contact.vertex;
fx = discrete[ifirst+ni][3] * j_a + fn[0] + ft[0];
fy = discrete[ifirst+ni][4] * j_a + fn[1] + ft[1];
fz = discrete[ifirst+ni][5] * j_a + fn[2] + ft[2];
f[ibody][0] += fx;
f[ibody][1] += fy;
f[ibody][2] += fz;
sum_torque(x[ibody], contact.xv, fx, fy, fz, torque[ibody]);
// accumulate forces to the vertex only
facc[0] += fx; facc[1] += fy; facc[2] += fz;
#ifdef _POLYGON_DEBUG
printf("From contact forces: vertex fx %f fy %f fz %f\n"
" torque body %d: %f %f %f\n",
discrete[ifirst+ni][3], discrete[ifirst+ni][4], discrete[ifirst+ni][5],
atom->tag[ibody],torque[ibody][0],torque[ibody][1],torque[ibody][2]);
#endif
}
/* ----------------------------------------------------------------------
Determine the length of the contact segment, i.e. the separation between
2 contacts
------------------------------------------------------------------------- */
double FixWallBodyPolygon::contact_separation(const Contact& c1,
const Contact& c2)
{
double x1 = c1.xv[0];
double y1 = c1.xv[1];
double x2 = c1.xe[0];
double y2 = c1.xe[1];
double x3 = c2.xv[0];
double y3 = c2.xv[1];
double delta_a = 0.0;
if (fabs(x2 - x1) > EPSILON) {
double A = (y2 - y1) / (x2 - x1);
delta_a = fabs(y1 - A * x1 - y3 + A * x3) / sqrt(1 + A * A);
} else {
delta_a = fabs(x1 - x3);
}
return delta_a;
}
/* ----------------------------------------------------------------------
Accumulate torque to body from the force f=(fx,fy,fz) acting at point x
------------------------------------------------------------------------- */
void FixWallBodyPolygon::sum_torque(double* xm, double *x, double fx,
double fy, double fz, double* torque)
{
double rx = x[0] - xm[0];
double ry = x[1] - xm[1];
double rz = x[2] - xm[2];
double tx = ry * fz - rz * fy;
double ty = rz * fx - rx * fz;
double tz = rx * fy - ry * fx;
torque[0] += tx;
torque[1] += ty;
torque[2] += tz;
}
/* ----------------------------------------------------------------------
Calculate the total velocity of a point (vertex, a point on an edge):
vi = vcm + omega ^ (p - xcm)
------------------------------------------------------------------------- */
void FixWallBodyPolygon::total_velocity(double* p, double *xcm,
double* vcm, double *angmom, double *inertia,
double *quat, double* vi)
{
double r[3],omega[3],ex_space[3],ey_space[3],ez_space[3];
r[0] = p[0] - xcm[0];
r[1] = p[1] - xcm[1];
r[2] = p[2] - xcm[2];
MathExtra::q_to_exyz(quat,ex_space,ey_space,ez_space);
MathExtra::angmom_to_omega(angmom,ex_space,ey_space,ez_space,
inertia,omega);
vi[0] = omega[1]*r[2] - omega[2]*r[1] + vcm[0];
vi[1] = omega[2]*r[0] - omega[0]*r[2] + vcm[1];
vi[2] = omega[0]*r[1] - omega[1]*r[0] + vcm[2];
}
/* ---------------------------------------------------------------------- */
void FixWallBodyPolygon::distance(const double* x2, const double* x1,
double& r) {
r = sqrt((x2[0] - x1[0]) * (x2[0] - x1[0])
+ (x2[1] - x1[1]) * (x2[1] - x1[1])
+ (x2[2] - x1[2]) * (x2[2] - x1[2]));
}

View File

@ -0,0 +1,131 @@
/* -*- c++ -*- ----------------------------------------------------------
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 FIX_CLASS
FixStyle(wall/body/polygon,FixWallBodyPolygon)
#else
#ifndef LMP_FIX_WALL_BODY_POLYGON_H
#define LMP_FIX_WALL_BODY_POLYGON_H
#include "fix.h"
namespace LAMMPS_NS {
class FixWallBodyPolygon : public Fix {
public:
FixWallBodyPolygon(class LAMMPS *, int, char **);
virtual ~FixWallBodyPolygon();
int setmask();
void init();
void setup(int);
virtual void post_force(int);
void reset_dt();
struct Contact {
int ibody, jbody; // body (i.e. atom) indices (not tags)
int vertex; // vertex of the first polygon
int edge; // edge of the second polygon
double xv[3]; // coordinates of the vertex
double xe[3]; // coordinates of the projection of the vertex on the edge
double separation;// separation at contact
};
protected:
int wallstyle,pairstyle,wiggle,axis;
double kn; // normal repulsion strength
double c_n; // normal damping coefficient
double c_t; // tangential damping coefficient
double lo,hi,cylradius;
double amplitude,period,omega;
double dt;
int time_origin;
class AtomVecBody *avec;
class BodyRoundedPolygon *bptr;
double **discrete; // list of all sub-particles for all bodies
int ndiscrete; // number of discretes in list
int dmax; // allocated size of discrete list
int *dnum; // number of discretes per line, 0 if uninit
int *dfirst; // index of first discrete per each line
int nmax; // allocated size of dnum,dfirst vectors
double **edge; // list of all edge for all bodies
int nedge; // number of edge in list
int edmax; // allocated size of edge list
int *ednum; // number of edges per line, 0 if uninit
int *edfirst; // index of first edge per each line
int ednummax; // allocated size of ednum,edfirst vectors
double *enclosing_radius; // enclosing radii for all bodies
double *rounded_radius; // rounded radii for all bodies
void body2space(int);
int vertex_against_wall(int ibody, double wall_pos, double** x,
double** f, double** torque, int side,
Contact* contact_list, int &num_contacts,
double* facc);
int compute_distance_to_wall(double* x0, double rradi, double wall_pos,
int side, double &d, double hi[3], int &contact);
double contact_separation(const Contact& c1, const Contact& c2);
void contact_forces(Contact& contact, double j_a, double** x,
double** v, double** angmom, double** f, double** torque,
double* vwall, double* facc);
void sum_torque(double* xm, double *x, double fx,
double fy, double fz, double* torque);
void total_velocity(double* p, double *xcm, double* vcm, double *angmom,
double *inertia, double *quat, double* vi);
void distance(const double* x2, const double* x1, double& r);
};
}
#endif
#endif
/* ERROR/WARNING messages:
E: Illegal ... command
Self-explanatory. Check the input script syntax and compare to the
documentation for the command. You can use -echo screen as a
command-line option when running LAMMPS to see the offending line.
E: Fix wall/body/polygon requires atom style body rounded/polygon
Self-explanatory.
E: Cannot use wall in periodic dimension
Self-explanatory.
E: Cannot wiggle and shear fix wall/body/polygon
Cannot specify both options at the same time.
E: Invalid wiggle direction for fix wall/body/polygon
Self-explanatory.
E: Fix wall/body/polygon is incompatible with Pair style
Must use a body pair style to define the parameters needed for
this fix.
*/

View File

@ -0,0 +1,945 @@
/* ----------------------------------------------------------------------
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.
------------------------------------------------------------------------- */
/* ----------------------------------------------------------------------
Contributing author: Trung Dac Nguyen (ndactrung@gmail.com)
------------------------------------------------------------------------- */
#include <cmath>
#include <cstdlib>
#include <cstring>
#include "fix_wall_body_polyhedron.h"
#include "atom.h"
#include "atom_vec_body.h"
#include "body_rounded_polyhedron.h"
#include "domain.h"
#include "update.h"
#include "force.h"
#include "pair.h"
#include "modify.h"
#include "respa.h"
#include "math_const.h"
#include "math_extra.h"
#include "memory.h"
#include "error.h"
using namespace LAMMPS_NS;
using namespace FixConst;
using namespace MathConst;
enum{XPLANE=0,YPLANE=1,ZPLANE}; // XYZ PLANE need to be 0,1,2
enum{HOOKE,HOOKE_HISTORY};
enum {INVALID=0,NONE=1,VERTEX=2};
enum {FAR=0,XLO,XHI,YLO,YHI,ZLO,ZHI};
//#define _POLYHEDRON_DEBUG
#define DELTA 10000
#define EPSILON 1e-2
#define BIG 1.0e20
#define MAX_CONTACTS 4 // maximum number of contacts for 2D models
#define EFF_CONTACTS 2 // effective contacts for 2D models
/* ---------------------------------------------------------------------- */
FixWallBodyPolyhedron::FixWallBodyPolyhedron(LAMMPS *lmp, int narg, char **arg) :
Fix(lmp, narg, arg)
{
if (narg < 7) error->all(FLERR,"Illegal fix wall/body/polyhedron command");
if (!atom->body_flag)
error->all(FLERR,"Fix wall/body/polyhedron requires "
"atom style body/rounded/polyhedron");
restart_peratom = 1;
create_attribute = 1;
// wall/particle coefficients
kn = force->numeric(FLERR,arg[3]);
c_n = force->numeric(FLERR,arg[4]);
if (strcmp(arg[5],"NULL") == 0) c_t = 0.5 * c_n;
else c_t = force->numeric(FLERR,arg[5]);
if (kn < 0.0 || c_n < 0.0 || c_t < 0.0)
error->all(FLERR,"Illegal fix wall/body/polyhedron command");
// wallstyle args
int iarg = 6;
if (strcmp(arg[iarg],"xplane") == 0) {
if (narg < iarg+3) error->all(FLERR,"Illegal fix wall/body/polyhedron command");
wallstyle = XPLANE;
if (strcmp(arg[iarg+1],"NULL") == 0) lo = -BIG;
else lo = force->numeric(FLERR,arg[iarg+1]);
if (strcmp(arg[iarg+2],"NULL") == 0) hi = BIG;
else hi = force->numeric(FLERR,arg[iarg+2]);
iarg += 3;
} else if (strcmp(arg[iarg],"yplane") == 0) {
if (narg < iarg+3) error->all(FLERR,"Illegal fix wall/body/polyhedron command");
wallstyle = YPLANE;
if (strcmp(arg[iarg+1],"NULL") == 0) lo = -BIG;
else lo = force->numeric(FLERR,arg[iarg+1]);
if (strcmp(arg[iarg+2],"NULL") == 0) hi = BIG;
else hi = force->numeric(FLERR,arg[iarg+2]);
iarg += 3;
} else if (strcmp(arg[iarg],"zplane") == 0) {
if (narg < iarg+3) error->all(FLERR,"Illegal fix wall/body/polyhedron command");
wallstyle = ZPLANE;
if (strcmp(arg[iarg+1],"NULL") == 0) lo = -BIG;
else lo = force->numeric(FLERR,arg[iarg+1]);
if (strcmp(arg[iarg+2],"NULL") == 0) hi = BIG;
else hi = force->numeric(FLERR,arg[iarg+2]);
iarg += 3;
}
// check for trailing keyword/values
wiggle = 0;
while (iarg < narg) {
if (strcmp(arg[iarg],"wiggle") == 0) {
if (iarg+4 > narg) error->all(FLERR,"Illegal fix wall/body/polyhedron command");
if (strcmp(arg[iarg+1],"x") == 0) axis = 0;
else if (strcmp(arg[iarg+1],"y") == 0) axis = 1;
else if (strcmp(arg[iarg+1],"z") == 0) axis = 2;
else error->all(FLERR,"Illegal fix wall/body/polyhedron command");
amplitude = force->numeric(FLERR,arg[iarg+2]);
period = force->numeric(FLERR,arg[iarg+3]);
wiggle = 1;
iarg += 4;
} else error->all(FLERR,"Illegal fix wall/body/polyhedron command");
}
if (wallstyle == XPLANE && domain->xperiodic)
error->all(FLERR,"Cannot use wall in periodic dimension");
if (wallstyle == YPLANE && domain->yperiodic)
error->all(FLERR,"Cannot use wall in periodic dimension");
if (wallstyle == ZPLANE && domain->zperiodic)
error->all(FLERR,"Cannot use wall in periodic dimension");
// setup oscillations
if (wiggle) omega = 2.0*MY_PI / period;
time_origin = update->ntimestep;
dmax = nmax = 0;
discrete = NULL;
dnum = dfirst = NULL;
edmax = ednummax = 0;
edge = NULL;
ednum = edfirst = NULL;
facmax = facnummax = 0;
face = NULL;
facnum = facfirst = NULL;
enclosing_radius = NULL;
rounded_radius = NULL;
}
/* ---------------------------------------------------------------------- */
FixWallBodyPolyhedron::~FixWallBodyPolyhedron()
{
memory->destroy(discrete);
memory->destroy(dnum);
memory->destroy(dfirst);
memory->destroy(edge);
memory->destroy(ednum);
memory->destroy(edfirst);
memory->destroy(face);
memory->destroy(facnum);
memory->destroy(facfirst);
memory->destroy(enclosing_radius);
memory->destroy(rounded_radius);
}
/* ---------------------------------------------------------------------- */
int FixWallBodyPolyhedron::setmask()
{
int mask = 0;
mask |= POST_FORCE;
return mask;
}
/* ---------------------------------------------------------------------- */
void FixWallBodyPolyhedron::init()
{
dt = update->dt;
avec = (AtomVecBody *) atom->style_match("body");
if (!avec)
error->all(FLERR,"Pair body/rounded/polyhedron requires atom style body");
if (strcmp(avec->bptr->style,"rounded/polyhedron") != 0)
error->all(FLERR,"Pair body/rounded/polyhedron requires "
"body style rounded/polyhedron");
bptr = (BodyRoundedPolyhedron *) avec->bptr;
// set pairstyle from body/polyhedronular pair style
if (force->pair_match("body/rounded/polyhedron",1))
pairstyle = HOOKE;
else error->all(FLERR,"Fix wall/body/polyhedron is incompatible with Pair style");
}
/* ---------------------------------------------------------------------- */
void FixWallBodyPolyhedron::setup(int vflag)
{
if (strstr(update->integrate_style,"verlet"))
post_force(vflag);
}
/* ---------------------------------------------------------------------- */
void FixWallBodyPolyhedron::post_force(int vflag)
{
double vwall[3],dx,dy,dz,del1,del2,delxy,delr,rsq,eradi,rradi,wall_pos;
int i,ni,npi,ifirst,nei,iefirst,nfi,iffirst,side;
double facc[3];
// set position of wall to initial settings and velocity to 0.0
// if wiggle, set wall position and velocity accordingly
double wlo = lo;
double whi = hi;
vwall[0] = vwall[1] = vwall[2] = 0.0;
if (wiggle) {
double arg = omega * (update->ntimestep - time_origin) * dt;
if (wallstyle == axis) {
wlo = lo + amplitude - amplitude*cos(arg);
whi = hi + amplitude - amplitude*cos(arg);
}
vwall[axis] = amplitude*omega*sin(arg);
}
// loop over all my atoms
// rsq = distance from wall
// dx,dy,dz = signed distance from wall
// for rotating cylinder, reset vwall based on particle position
// skip atom if not close enough to wall
// if wall was set to NULL, it's skipped since lo/hi are infinity
// compute force and torque on atom if close enough to wall
// via wall potential matched to pair potential
double **x = atom->x;
double **v = atom->v;
double **f = atom->f;
int *body = atom->body;
double *radius = atom->radius;
double **torque = atom->torque;
double **angmom = atom->angmom;
int *mask = atom->mask;
int nlocal = atom->nlocal;
// grow the per-atom lists if necessary and initialize
if (atom->nmax > nmax) {
memory->destroy(dnum);
memory->destroy(dfirst);
memory->destroy(ednum);
memory->destroy(edfirst);
memory->destroy(facnum);
memory->destroy(facfirst);
memory->destroy(enclosing_radius);
memory->destroy(rounded_radius);
nmax = atom->nmax;
memory->create(dnum,nmax,"fix:dnum");
memory->create(dfirst,nmax,"fix:dfirst");
memory->create(ednum,nmax,"fix:ednum");
memory->create(edfirst,nmax,"fix:edfirst");
memory->create(facnum,nmax,"fix:facnum");
memory->create(facfirst,nmax,"fix:facfirst");
memory->create(enclosing_radius,nmax,"fix:enclosing_radius");
memory->create(rounded_radius,nmax,"fix:rounded_radius");
}
ndiscrete = nedge = nface = 0;
for (i = 0; i < nlocal; i++)
dnum[i] = ednum[i] = facnum[i] = 0;
for (i = 0; i < nlocal; i++) {
if (mask[i] & groupbit) {
if (body[i] < 0) continue;
dx = dy = dz = 0.0;
side = FAR;
if (wallstyle == XPLANE) {
del1 = x[i][0] - wlo;
del2 = whi - x[i][0];
if (del1 < del2) {
dx = del1;
wall_pos = wlo;
side = XLO;
} else {
dx = -del2;
wall_pos = whi;
side = XHI;
}
} else if (wallstyle == YPLANE) {
del1 = x[i][1] - wlo;
del2 = whi - x[i][1];
if (del1 < del2) {
dy = del1;
wall_pos = wlo;
side = YLO;
} else {
dy = -del2;
wall_pos = whi;
side = YHI;
}
} else if (wallstyle == ZPLANE) {
del1 = x[i][2] - wlo;
del2 = whi - x[i][2];
if (del1 < del2) {
dy = del1;
wall_pos = wlo;
side = ZLO;
} else {
dy = -del2;
wall_pos = whi;
side = ZHI;
}
}
rsq = dx*dx + dy*dy + dz*dz;
if (rsq > radius[i]*radius[i]) continue;
double r = sqrt(rsq);
double rsqinv = 1.0 / rsq;
if (dnum[i] == 0) body2space(i);
npi = dnum[i];
ifirst = dfirst[i];
nei = ednum[i];
iefirst = edfirst[i];
nfi = facnum[i];
iffirst = facfirst[i];
eradi = enclosing_radius[i];
rradi = rounded_radius[i];
if (npi == 1) {
sphere_against_wall(i, wall_pos, side, vwall, x, v, f, angmom, torque);
continue;
}
// reset vertex and edge forces
for (ni = 0; ni < npi; ni++) {
discrete[ifirst+ni][3] = 0;
discrete[ifirst+ni][4] = 0;
discrete[ifirst+ni][5] = 0;
discrete[ifirst+ni][6] = 0;
}
for (ni = 0; ni < nei; ni++) {
edge[iefirst+ni][2] = 0;
edge[iefirst+ni][3] = 0;
edge[iefirst+ni][4] = 0;
edge[iefirst+ni][5] = 0;
}
int interact, num_contacts, done;
double delta_a, delta_ua, j_a;
Contact contact_list[MAX_CONTACTS];
num_contacts = 0;
facc[0] = facc[1] = facc[2] = 0;
interact = edge_against_wall(i, wall_pos, side, vwall, x, f, torque,
contact_list, num_contacts, facc);
} // group bit
}
}
/* ---------------------------------------------------------------------- */
void FixWallBodyPolyhedron::reset_dt()
{
dt = update->dt;
}
/* ----------------------------------------------------------------------
convert N sub-particles in body I to space frame using current quaternion
store sub-particle space-frame displacements from COM in discrete list
------------------------------------------------------------------------- */
void FixWallBodyPolyhedron::body2space(int i)
{
int ibonus = atom->body[i];
AtomVecBody::Bonus *bonus = &avec->bonus[ibonus];
int nsub = bptr->nsub(bonus);
double *coords = bptr->coords(bonus);
int body_num_edges = bptr->nedges(bonus);
double* edge_ends = bptr->edges(bonus);
int body_num_faces = bptr->nfaces(bonus);
double* face_pts = bptr->faces(bonus);
double eradius = bptr->enclosing_radius(bonus);
double rradius = bptr->rounded_radius(bonus);
// get the number of sub-particles (vertices)
// and the index of the first vertex of my body in the list
dnum[i] = nsub;
dfirst[i] = ndiscrete;
// grow the vertex list if necessary
// the first 3 columns are for coords, the last 3 for forces
if (ndiscrete + nsub > dmax) {
dmax += DELTA;
memory->grow(discrete,dmax,7,"fix:discrete");
}
double p[3][3];
MathExtra::quat_to_mat(bonus->quat,p);
for (int m = 0; m < nsub; m++) {
MathExtra::matvec(p,&coords[3*m],discrete[ndiscrete]);
discrete[ndiscrete][3] = 0;
discrete[ndiscrete][4] = 0;
discrete[ndiscrete][5] = 0;
discrete[ndiscrete][6] = 0;
ndiscrete++;
}
// get the number of edges (vertices)
// and the index of the first edge of my body in the list
ednum[i] = body_num_edges;
edfirst[i] = nedge;
// grow the edge list if necessary
// the first 2 columns are for vertex indices within body,
// the last 3 for forces
if (nedge + body_num_edges > edmax) {
edmax += DELTA;
memory->grow(edge,edmax,6,"fix:edge");
}
for (int m = 0; m < body_num_edges; m++) {
edge[nedge][0] = static_cast<int>(edge_ends[2*m+0]);
edge[nedge][1] = static_cast<int>(edge_ends[2*m+1]);
edge[nedge][2] = 0;
edge[nedge][3] = 0;
edge[nedge][4] = 0;
edge[nedge][5] = 0;
nedge++;
}
// get the number of faces and the index of the first face
facnum[i] = body_num_faces;
facfirst[i] = nface;
// grow the face list if necessary
// the first 3 columns are for vertex indices within body, the last 3 for forces
if (nface + body_num_faces > facmax) {
facmax += DELTA;
memory->grow(face,facmax,6,"pair:face");
}
for (int m = 0; m < body_num_faces; m++) {
face[nface][0] = static_cast<int>(face_pts[3*m+0]);
face[nface][1] = static_cast<int>(face_pts[3*m+1]);
face[nface][2] = static_cast<int>(face_pts[3*m+2]);
face[nface][3] = 0;
face[nface][4] = 0;
face[nface][5] = 0;
nface++;
}
enclosing_radius[i] = eradius;
rounded_radius[i] = rradius;
}
/* ----------------------------------------------------------------------
Determine the interaction mode between a sphere against the wall
i = atom i (body i)
x = atoms' coordinates
f = atoms' forces
torque = atoms' torques
---------------------------------------------------------------------- */
int FixWallBodyPolyhedron::sphere_against_wall(int i, double wall_pos,
int side, double* vwall, double** x, double** v, double** f,
double** angmom, double** torque)
{
int mode;
double rradi,hi[3],d,delx,dely,delz,R,fpair,fx,fy,fz;
rradi = rounded_radius[i];
mode = NONE;
if (wallstyle == XPLANE) {
hi[0] = wall_pos;
hi[1] = x[i][1];
hi[2] = x[i][2];
} else if (wallstyle == YPLANE) {
hi[0] = x[i][0];
hi[1] = wall_pos;
hi[2] = x[i][2];
} else if (wallstyle == ZPLANE) {
hi[0] = x[i][0];
hi[1] = x[i][1];
hi[2] = wall_pos;
}
distance(hi, x[i], d);
if (d <= rradi) {
delx = x[i][0] - hi[0];
dely = x[i][1] - hi[1];
delz = x[i][2] - hi[2];
R = d - rradi;
fpair = -kn * R;
fx = delx*fpair/d;
fy = dely*fpair/d;
fz = delz*fpair/d;
contact_forces(i, 1.0, x[i], hi, delx, dely, delz,
fx, fy, fz, x, v, angmom, f, torque, vwall);
mode = VERTEX;
}
return mode;
}
/* ----------------------------------------------------------------------
Determine the interaction mode between i's vertices against the wall
i = atom i (body i)
x = atoms' coordinates
f = atoms' forces
torque = atoms' torques
Output:
contact_list = list of contacts between i and the wall
num_contacts = number of contacts between i's vertices and the wall
Return:
number of contacts of the edge to the wall (0, 1 or 2)
---------------------------------------------------------------------- */
int FixWallBodyPolyhedron::edge_against_wall(int i, double wall_pos,
int side, double* vwall, double** x, double** f, double** torque,
Contact* contact_list, int &num_contacts, double* facc)
{
int ni, nei, mode, contact;
double rradi;
int nlocal = atom->nlocal;
nei = ednum[i];
rradi = rounded_radius[i];
contact = 0;
// loop through body i's edges
for (ni = 0; ni < nei; ni++)
mode = compute_distance_to_wall(i, ni, x[i], rradi, wall_pos, side, vwall,
contact);
return contact;
}
/* -------------------------------------------------------------------------
Compute the distance between a vertex to the wall
another body
Input:
x0 = coordinate of the tested vertex
rradi = rounded radius of the vertex
wall_pos = position of the wall
Output:
d = Distance from a point x0 to an wall
hi = coordinates of the projection of x0 on the wall
contact = 0 no contact between the queried vertex and the wall
1 contact detected
return NONE if there is no interaction
VERTEX if the tested vertex interacts with the wall
------------------------------------------------------------------------- */
int FixWallBodyPolyhedron::compute_distance_to_wall(int ibody, int edge_index,
double *xmi, double rounded_radius_i, double wall_pos,
int side, double* vwall, int &contact)
{
int mode,ifirst,iefirst,npi1,npi2;
double d1,d2,xpi1[3],xpi2[3],hi[3];
double fx,fy,fz,fpair,delx,dely,delz,R;
double** x = atom->x;
double** v = atom->v;
double** f = atom->f;
double** torque = atom->torque;
double** angmom = atom->angmom;
// two ends of the edge from body i
ifirst = dfirst[ibody];
iefirst = edfirst[ibody];
npi1 = static_cast<int>(edge[iefirst+edge_index][0]);
npi2 = static_cast<int>(edge[iefirst+edge_index][1]);
xpi1[0] = xmi[0] + discrete[ifirst+npi1][0];
xpi1[1] = xmi[1] + discrete[ifirst+npi1][1];
xpi1[2] = xmi[2] + discrete[ifirst+npi1][2];
xpi2[0] = xmi[0] + discrete[ifirst+npi2][0];
xpi2[1] = xmi[1] + discrete[ifirst+npi2][1];
xpi2[2] = xmi[2] + discrete[ifirst+npi2][2];
// determine the intersection of the edge to the wall
mode = NONE;
double j_a = 1.0;
if (wallstyle == XPLANE) {
hi[0] = wall_pos;
hi[1] = xpi1[1];
hi[2] = xpi1[2];
} else if (wallstyle == YPLANE) {
hi[0] = xpi1[0];
hi[1] = wall_pos;
hi[2] = xpi1[2];
} else if (wallstyle == ZPLANE) {
hi[0] = xpi1[0];
hi[1] = xpi1[1];
hi[2] = wall_pos;
}
distance(hi, xpi1, d1);
if (d1 <= rounded_radius_i && static_cast<int>(discrete[ifirst+npi1][6]) == 0) {
delx = xpi1[0] - hi[0];
dely = xpi1[1] - hi[1];
delz = xpi1[2] - hi[2];
R = d1 - rounded_radius_i;
fpair = -kn * R;
fx = delx*fpair/d1;
fy = dely*fpair/d1;
fz = delz*fpair/d1;
contact_forces(ibody, j_a, xpi1, hi, delx, dely, delz,
fx, fy, fz, x, v, angmom, f, torque, vwall);
discrete[ifirst+npi1][6] = 1;
contact++;
mode = VERTEX;
}
if (wallstyle == XPLANE) {
hi[0] = wall_pos;
hi[1] = xpi2[1];
hi[2] = xpi2[2];
} else if (wallstyle == YPLANE) {
hi[0] = xpi2[0];
hi[1] = wall_pos;
hi[2] = xpi2[2];
} else if (wallstyle == ZPLANE) {
hi[0] = xpi2[0];
hi[1] = xpi2[1];
hi[2] = wall_pos;
}
distance(hi, xpi2, d2);
if (d2 <= rounded_radius_i && static_cast<int>(discrete[ifirst+npi2][6]) == 0) {
delx = xpi2[0] - hi[0];
dely = xpi2[1] - hi[1];
delz = xpi2[2] - hi[2];
R = d2 - rounded_radius_i;
fpair = -kn * R;
fx = delx*fpair/d2;
fy = dely*fpair/d2;
fz = delz*fpair/d2;
contact_forces(ibody, j_a, xpi2, hi, delx, dely, delz,
fx, fy, fz, x, v, angmom, f, torque, vwall);
discrete[ifirst+npi2][6] = 1;
contact++;
mode = VERTEX;
}
return mode;
}
/* ----------------------------------------------------------------------
Compute contact forces between two bodies
modify the force stored at the vertex and edge in contact by j_a
sum forces and torque to the corresponding bodies
fn = normal friction component
ft = tangential friction component (-c_t * v_t)
------------------------------------------------------------------------- */
void FixWallBodyPolyhedron::contact_forces(int ibody,
double j_a, double *xi, double *xj, double delx, double dely, double delz,
double fx, double fy, double fz, double** x, double** v, double** angmom,
double** f, double** torque, double* vwall)
{
int ibonus,jbonus;
double fxt,fyt,fzt,rsq,rsqinv;
double vr1,vr2,vr3,vnnr,vn1,vn2,vn3,vt1,vt2,vt3;
double fn[3],ft[3],vi[3],vj[3];
double *quat, *inertia;
AtomVecBody::Bonus *bonus;
// compute the velocity of the vertex in the space-fixed frame
ibonus = atom->body[ibody];
bonus = &avec->bonus[ibonus];
quat = bonus->quat;
inertia = bonus->inertia;
total_velocity(xi, x[ibody], v[ibody], angmom[ibody],
inertia, quat, vi);
// vector pointing from the contact point on ibody to the wall
rsq = delx*delx + dely*dely + delz*delz;
rsqinv = 1.0/rsq;
// relative translational velocity
vr1 = vi[0] - vwall[0];
vr2 = vi[1] - vwall[1];
vr3 = vi[2] - vwall[2];
// normal component
vnnr = vr1*delx + vr2*dely + vr3*delz;
vn1 = delx*vnnr * rsqinv;
vn2 = dely*vnnr * rsqinv;
vn3 = delz*vnnr * rsqinv;
// tangential component
vt1 = vr1 - vn1;
vt2 = vr2 - vn2;
vt3 = vr3 - vn3;
// normal friction term at contact
fn[0] = -c_n * vn1;
fn[1] = -c_n * vn2;
fn[2] = -c_n * vn3;
// tangential friction term at contact
// excluding the tangential deformation term for now
ft[0] = -c_t * vt1;
ft[1] = -c_t * vt2;
ft[2] = -c_t * vt3;
fxt = fx; fyt = fy; fzt = fz;
fx = fxt * j_a + fn[0] + ft[0];
fy = fyt * j_a + fn[1] + ft[1];
fz = fzt * j_a + fn[2] + ft[2];
f[ibody][0] += fx;
f[ibody][1] += fy;
f[ibody][2] += fz;
sum_torque(x[ibody], xi, fx, fy, fz, torque[ibody]);
#ifdef _POLYHEDRON_DEBUG
printf("From contact forces: vertex fx %f fy %f fz %f\n"
" torque body %d: %f %f %f\n"
" torque body %d: %f %f %f\n",
fxt, fyt, fzt,
atom->tag[ibody],torque[ibody][0],torque[ibody][1],torque[ibody][2],
atom->tag[jbody],torque[jbody][0],torque[jbody][1],torque[jbody][2]);
#endif
}
/* ----------------------------------------------------------------------
Compute the contact forces between two bodies
modify the force stored at the vertex and edge in contact by j_a
sum forces and torque to the corresponding bodies
fn = normal friction component
ft = tangential friction component (-c_t * vrt)
------------------------------------------------------------------------- */
void FixWallBodyPolyhedron::contact_forces(Contact& contact, double j_a,
double** x, double** v, double** angmom, double** f,
double** torque, double* vwall, double* facc)
{
int ibody,ibonus,ifirst, jefirst, ni;
double fx,fy,fz,delx,dely,delz,rsq,rsqinv;
double vr1,vr2,vr3,vnnr,vn1,vn2,vn3,vt1,vt2,vt3;
double fn[3],ft[3],vi[3];
double *quat, *inertia;
AtomVecBody::Bonus *bonus;
ibody = contact.ibody;
// compute the velocity of the vertex in the space-fixed frame
ibonus = atom->body[ibody];
bonus = &avec->bonus[ibonus];
quat = bonus->quat;
inertia = bonus->inertia;
total_velocity(contact.xv, x[ibody], v[ibody], angmom[ibody],
inertia, quat, vi);
// vector pointing from the vertex to the point on the wall
delx = contact.xv[0] - contact.xe[0];
dely = contact.xv[1] - contact.xe[1];
delz = contact.xv[2] - contact.xe[2];
rsq = delx*delx + dely*dely + delz*delz;
rsqinv = 1.0/rsq;
// relative translational velocity
vr1 = vi[0] - vwall[0];
vr2 = vi[1] - vwall[1];
vr3 = vi[2] - vwall[2];
// normal component
vnnr = vr1*delx + vr2*dely + vr3*delz;
vn1 = delx*vnnr * rsqinv;
vn2 = dely*vnnr * rsqinv;
vn3 = delz*vnnr * rsqinv;
// tangential component
vt1 = vr1 - vn1;
vt2 = vr2 - vn2;
vt3 = vr3 - vn3;
// normal friction term at contact
fn[0] = -c_n * vn1;
fn[1] = -c_n * vn2;
fn[2] = -c_n * vn3;
// tangential friction term at contact
// excluding the tangential deformation term for now
ft[0] = -c_t * vt1;
ft[1] = -c_t * vt2;
ft[2] = -c_t * vt3;
// only the cohesive force is scaled by j_a
ifirst = dfirst[ibody];
ni = contact.vertex;
fx = discrete[ifirst+ni][3] * j_a + fn[0] + ft[0];
fy = discrete[ifirst+ni][4] * j_a + fn[1] + ft[1];
fz = discrete[ifirst+ni][5] * j_a + fn[2] + ft[2];
f[ibody][0] += fx;
f[ibody][1] += fy;
f[ibody][2] += fz;
sum_torque(x[ibody], contact.xv, fx, fy, fz, torque[ibody]);
// accumulate forces to the vertex only
facc[0] += fx; facc[1] += fy; facc[2] += fz;
#ifdef _POLYHEDRON_DEBUG
printf("From contact forces: vertex fx %f fy %f fz %f\n"
" torque body %d: %f %f %f\n",
discrete[ifirst+ni][3], discrete[ifirst+ni][4], discrete[ifirst+ni][5],
atom->tag[ibody],torque[ibody][0],torque[ibody][1],torque[ibody][2]);
#endif
}
/* ----------------------------------------------------------------------
Determine the length of the contact segment, i.e. the separation between
2 contacts
------------------------------------------------------------------------- */
double FixWallBodyPolyhedron::contact_separation(const Contact& c1,
const Contact& c2)
{
double x1 = c1.xv[0];
double y1 = c1.xv[1];
double x2 = c1.xe[0];
double y2 = c1.xe[1];
double x3 = c2.xv[0];
double y3 = c2.xv[1];
double delta_a = 0.0;
if (fabs(x2 - x1) > EPSILON) {
double A = (y2 - y1) / (x2 - x1);
delta_a = fabs(y1 - A * x1 - y3 + A * x3) / sqrt(1 + A * A);
} else {
delta_a = fabs(x1 - x3);
}
return delta_a;
}
/* ----------------------------------------------------------------------
Accumulate torque to body from the force f=(fx,fy,fz) acting at point x
------------------------------------------------------------------------- */
void FixWallBodyPolyhedron::sum_torque(double* xm, double *x, double fx,
double fy, double fz, double* torque)
{
double rx = x[0] - xm[0];
double ry = x[1] - xm[1];
double rz = x[2] - xm[2];
double tx = ry * fz - rz * fy;
double ty = rz * fx - rx * fz;
double tz = rx * fy - ry * fx;
torque[0] += tx;
torque[1] += ty;
torque[2] += tz;
}
/* ----------------------------------------------------------------------
Calculate the total velocity of a point (vertex, a point on an edge):
vi = vcm + omega ^ (p - xcm)
------------------------------------------------------------------------- */
void FixWallBodyPolyhedron::total_velocity(double* p, double *xcm,
double* vcm, double *angmom, double *inertia,
double *quat, double* vi)
{
double r[3],omega[3],ex_space[3],ey_space[3],ez_space[3];
r[0] = p[0] - xcm[0];
r[1] = p[1] - xcm[1];
r[2] = p[2] - xcm[2];
MathExtra::q_to_exyz(quat,ex_space,ey_space,ez_space);
MathExtra::angmom_to_omega(angmom,ex_space,ey_space,ez_space,
inertia,omega);
vi[0] = omega[1]*r[2] - omega[2]*r[1] + vcm[0];
vi[1] = omega[2]*r[0] - omega[0]*r[2] + vcm[1];
vi[2] = omega[0]*r[1] - omega[1]*r[0] + vcm[2];
}
/* ---------------------------------------------------------------------- */
void FixWallBodyPolyhedron::distance(const double* x2, const double* x1,
double& r) {
r = sqrt((x2[0] - x1[0]) * (x2[0] - x1[0])
+ (x2[1] - x1[1]) * (x2[1] - x1[1])
+ (x2[2] - x1[2]) * (x2[2] - x1[2]));
}

View File

@ -0,0 +1,143 @@
/* -*- c++ -*- ----------------------------------------------------------
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 FIX_CLASS
FixStyle(wall/body/polyhedron,FixWallBodyPolyhedron)
#else
#ifndef LMP_FIX_WALL_BODY_POLYHERON_H
#define LMP_FIX_WALL_BODY_POLYHERON_H
#include "fix.h"
namespace LAMMPS_NS {
class FixWallBodyPolyhedron : public Fix {
public:
FixWallBodyPolyhedron(class LAMMPS *, int, char **);
virtual ~FixWallBodyPolyhedron();
int setmask();
void init();
void setup(int);
virtual void post_force(int);
void reset_dt();
struct Contact {
int ibody, jbody; // body (i.e. atom) indices (not tags)
int vertex; // vertex of the first polygon
int edge; // edge of the second polygon
double xv[3]; // coordinates of the vertex
double xe[3]; // coordinates of the projection of the vertex on the edge
double separation;// separation at contact
};
protected:
int wallstyle,pairstyle,wiggle,axis;
double kn,c_n,c_t;
double lo,hi,cylradius;
double amplitude,period,omega;
double dt;
int time_origin;
class AtomVecBody *avec;
class BodyRoundedPolyhedron *bptr;
double **discrete; // list of all sub-particles for all bodies
int ndiscrete; // number of discretes in list
int dmax; // allocated size of discrete list
int *dnum; // number of discretes per line, 0 if uninit
int *dfirst; // index of first discrete per each line
int nmax; // allocated size of dnum,dfirst vectors
double **edge; // list of all edge for all bodies
int nedge; // number of edge in list
int edmax; // allocated size of edge list
int *ednum; // number of edges per line, 0 if uninit
int *edfirst; // index of first edge per each line
int ednummax; // allocated size of ednum,edfirst vectors
double **face; // list of all edge for all bodies
int nface; // number of faces in list
int facmax; // allocated size of face list
int *facnum; // number of faces per line, 0 if uninit
int *facfirst; // index of first face per each line
int facnummax; // allocated size of facnum,facfirst vectors
double *enclosing_radius; // enclosing radii for all bodies
double *rounded_radius; // rounded radii for all bodies
void body2space(int);
int edge_against_wall(int ibody, double wall_pos, int side, double* vwall,
double** x, double** f, double** torque, Contact* contact_list,
int &num_contacts, double* facc);
int sphere_against_wall(int i, double wall_pos, int side, double* vwall,
double** x, double** v, double** f, double** angmom, double** torque);
int compute_distance_to_wall(int ibody, int edge_index, double *xmi,
double rounded_radius_i, double wall_pos, int side,
double* vwall, int &contact);
double contact_separation(const Contact& c1, const Contact& c2);
void contact_forces(int ibody, double j_a, double *xi, double *xj,
double delx, double dely, double delz,
double fx, double fy, double fz, double** x, double** v,
double** angmom, double** f, double** torque, double* vwall);
void contact_forces(Contact& contact, double j_a, double** x,
double** v, double** angmom, double** f, double** torque,
double* vwall, double* facc);
void sum_torque(double* xm, double *x, double fx,
double fy, double fz, double* torque);
void total_velocity(double* p, double *xcm, double* vcm, double *angmom,
double *inertia, double *quat, double* vi);
void distance(const double* x2, const double* x1, double& r);
};
}
#endif
#endif
/* ERROR/WARNING messages:
E: Illegal ... command
Self-explanatory. Check the input script syntax and compare to the
documentation for the command. You can use -echo screen as a
command-line option when running LAMMPS to see the offending line.
E: Fix wall/body/polyhedron requires atom style body rounded/polyhedron
Self-explanatory.
E: Cannot use wall in periodic dimension
Self-explanatory.
E: Cannot wiggle and shear fix wall/body/polygon
Cannot specify both options at the same time.
E: Invalid wiggle direction for fix wall/body/polygon
Self-explanatory.
E: Fix wall/body/polygon is incompatible with Pair style
Must use a body pair style to define the parameters needed for
this fix.
*/

View File

@ -15,7 +15,7 @@
#include <cstdio>
#include <cstdlib>
#include <cstring>
#include "pair_body.h"
#include "pair_body_nparticle.h"
#include "math_extra.h"
#include "atom.h"
#include "atom_vec_body.h"
@ -32,7 +32,7 @@ using namespace LAMMPS_NS;
/* ---------------------------------------------------------------------- */
PairBody::PairBody(LAMMPS *lmp) : Pair(lmp)
PairBodyNparticle::PairBodyNparticle(LAMMPS *lmp) : Pair(lmp)
{
dmax = nmax = 0;
discrete = NULL;
@ -44,7 +44,7 @@ PairBody::PairBody(LAMMPS *lmp) : Pair(lmp)
/* ---------------------------------------------------------------------- */
PairBody::~PairBody()
PairBodyNparticle::~PairBodyNparticle()
{
memory->destroy(discrete);
memory->destroy(dnum);
@ -66,7 +66,7 @@ PairBody::~PairBody()
/* ---------------------------------------------------------------------- */
void PairBody::compute(int eflag, int vflag)
void PairBodyNparticle::compute(int eflag, int vflag)
{
int i,j,ii,jj,inum,jnum,itype,jtype;
int ni,nj,npi,npj,ifirst,jfirst;
@ -336,7 +336,7 @@ void PairBody::compute(int eflag, int vflag)
allocate all arrays
------------------------------------------------------------------------- */
void PairBody::allocate()
void PairBodyNparticle::allocate()
{
allocated = 1;
int n = atom->ntypes;
@ -361,7 +361,7 @@ void PairBody::allocate()
global settings
------------------------------------------------------------------------- */
void PairBody::settings(int narg, char **arg)
void PairBodyNparticle::settings(int narg, char **arg)
{
if (narg != 1) error->all(FLERR,"Illegal pair_style command");
@ -381,7 +381,7 @@ void PairBody::settings(int narg, char **arg)
set coeffs for one or more type pairs
------------------------------------------------------------------------- */
void PairBody::coeff(int narg, char **arg)
void PairBodyNparticle::coeff(int narg, char **arg)
{
if (narg < 4 || narg > 5)
error->all(FLERR,"Incorrect args for pair coefficients");
@ -415,12 +415,12 @@ void PairBody::coeff(int narg, char **arg)
init specific to this pair style
------------------------------------------------------------------------- */
void PairBody::init_style()
void PairBodyNparticle::init_style()
{
avec = (AtomVecBody *) atom->style_match("body");
if (!avec) error->all(FLERR,"Pair body requires atom style body");
if (!avec) error->all(FLERR,"Pair body/nparticle requires atom style body");
if (strcmp(avec->bptr->style,"nparticle") != 0)
error->all(FLERR,"Pair body requires body style nparticle");
error->all(FLERR,"Pair body/nparticle requires body style nparticle");
bptr = (BodyNparticle *) avec->bptr;
neighbor->request(this,instance_me);
@ -430,7 +430,7 @@ void PairBody::init_style()
init for one type pair i,j and corresponding j,i
------------------------------------------------------------------------- */
double PairBody::init_one(int i, int j)
double PairBodyNparticle::init_one(int i, int j)
{
if (setflag[i][j] == 0) {
epsilon[i][j] = mix_energy(epsilon[i][i],epsilon[j][j],
@ -459,7 +459,7 @@ double PairBody::init_one(int i, int j)
store sub-particle space-frame displacements from COM in discrete list
------------------------------------------------------------------------- */
void PairBody::body2space(int i)
void PairBodyNparticle::body2space(int i)
{
int ibonus = atom->body[i];
AtomVecBody::Bonus *bonus = &avec->bonus[ibonus];

View File

@ -13,21 +13,21 @@
#ifdef PAIR_CLASS
PairStyle(body,PairBody)
PairStyle(body/nparticle,PairBodyNparticle)
#else
#ifndef LMP_PAIR_BODY_H
#define LMP_PAIR_BODY_H
#ifndef LMP_PAIR_BODY_NPARTICLE_H
#define LMP_PAIR_BODY_NPARTICLE_H
#include "pair.h"
namespace LAMMPS_NS {
class PairBody : public Pair {
class PairBodyNparticle : public Pair {
public:
PairBody(class LAMMPS *);
~PairBody();
PairBodyNparticle(class LAMMPS *);
~PairBodyNparticle();
void compute(int, int);
void settings(int, char **);
void coeff(int, char **);

File diff suppressed because it is too large Load Diff

View File

@ -0,0 +1,137 @@
/* -*- c++ -*- ----------------------------------------------------------
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 PAIR_CLASS
PairStyle(body/rounded/polygon,PairBodyRoundedPolygon)
#else
#ifndef LMP_PAIR_BODY_ROUNDED_POLYGON_H
#define LMP_PAIR_BODY_ROUNDED_POLYGON_H
#include "pair.h"
namespace LAMMPS_NS {
class PairBodyRoundedPolygon : public Pair {
public:
PairBodyRoundedPolygon(class LAMMPS *);
~PairBodyRoundedPolygon();
void compute(int, int);
void settings(int, char **);
void coeff(int, char **);
void init_style();
double init_one(int, int);
struct Contact {
int ibody, jbody; // body (i.e. atom) indices (not tags)
int vertex; // vertex of the first polygon
int edge; // edge of the second polygon
double xv[3]; // coordinates of the vertex
double xe[3]; // coordinates of the projection of the vertex on the edge
double separation;// separation at contact
};
protected:
double **k_n; // normal repulsion strength
double **k_na; // normal attraction strength
double c_n; // normal damping coefficient
double c_t; // tangential damping coefficient
double mu; // normal friction coefficient during gross sliding
double delta_ua; // contact line (area for 3D models) modification factor
double cut_inner; // cutoff for interaction between vertex-edge surfaces
class AtomVecBody *avec;
class BodyRoundedPolygon *bptr;
double **discrete; // list of all sub-particles for all bodies
int ndiscrete; // number of discretes in list
int dmax; // allocated size of discrete list
int *dnum; // number of discretes per line, 0 if uninit
int *dfirst; // index of first discrete per each line
int nmax; // allocated size of dnum,dfirst vectors
double **edge; // list of all edge for all bodies
int nedge; // number of edge in list
int edmax; // allocated size of edge list
int *ednum; // number of edges per line, 0 if uninit
int *edfirst; // index of first edge per each line
int ednummax; // allocated size of ednum,edfirst vectors
double *enclosing_radius; // enclosing radii for all bodies
double *rounded_radius; // rounded radii for all bodies
double *maxerad; // per-type maximum enclosing radius
void allocate();
void body2space(int);
// sphere-sphere interaction
void sphere_against_sphere(int i, int j, double delx, double dely, double delz,
double rsq, double k_n, double k_na,
double** x, double** v, double** f, int evflag);
// vertex-edge interaction
int vertex_against_edge(int i, int j, double k_n, double k_na,
double** x, double** f, double** torque,
tagint* tag, Contact* contact_list,
int &num_contacts, double &evdwl, double* facc);
// compute distance between a point and an edge from another body
int compute_distance_to_vertex(int ibody, int edge_index, double* xmi,
double rounded_radius, double* x0,
double x0_rounded_radius, double cut_inner,
double &d, double hi[3], double &t,
int &contact);
// compute contact forces if contact points are detected
void contact_forces(Contact& contact, double j_a,
double** x, double** v, double** angmom, double** f,
double** torque, double &evdwl, double* facc);
// compute the separation between two contacts
double contact_separation(const Contact& c1, const Contact& c2);
// accumulate torque to a body given a force at a given point
void sum_torque(double* xm, double *x, double fx,
double fy, double fz, double* torque);
// helper functions
int opposite_sides(double* x1, double* x2, double* a, double* b);
void total_velocity(double* p, double *xcm, double* vcm, double *angmom,
double *inertia, double *quat, double* vi);
inline void distance(const double* x2, const double* x1, double& r);
};
}
#endif
#endif
/* ERROR/WARNING messages:
E: Illegal ... command
Self-explanatory. Check the input script syntax and compare to the
documentation for the command. You can use -echo screen as a
command-line option when running LAMMPS to see the offending line.
E: Incorrect args for pair coefficients
Self-explanatory. Check the input script or data file.
E: Pair body/rounded/polygon requires atom style body rounded/polygon
Self-explanatory.
E: Pair body requires body style rounded/polygon
This pair style is specific to the rounded/polygon body style.
*/

File diff suppressed because it is too large Load Diff

View File

@ -0,0 +1,200 @@
/* -*- c++ -*- ----------------------------------------------------------
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 PAIR_CLASS
PairStyle(body/rounded/polyhedron,PairBodyRoundedPolyhedron)
#else
#ifndef LMP_PAIR_BODY_ROUNDED_POLYHEDRON_H
#define LMP_PAIR_BODY_ROUNDED_POLYHEDRON_H
#include "pair.h"
namespace LAMMPS_NS {
class PairBodyRoundedPolyhedron : public Pair {
public:
PairBodyRoundedPolyhedron(class LAMMPS *);
~PairBodyRoundedPolyhedron();
void compute(int, int);
void settings(int, char **);
void coeff(int, char **);
void init_style();
double init_one(int, int);
virtual void kernel_force(double R, int itype, int jtype,
double& energy, double& fpair);
struct Contact {
int ibody, jbody; // body (i.e. atom) indices (not tags)
int type; // 0 = VERTEX-FACE; 1 = EDGE-EDGE
double fx,fy,fz; // unscaled cohesive forces at contact
double xi[3]; // coordinates of the contact point on ibody
double xj[3]; // coordinates of the contact point on jbody
double separation; // contact surface separation
int unique;
};
protected:
double **k_n; // normal repulsion strength
double **k_na; // normal attraction strength
double c_n; // normal damping coefficient
double c_t; // tangential damping coefficient
double mu; // normal friction coefficient during gross sliding
double A_ua; // characteristic contact area
double cut_inner; // cutoff for interaction between vertex-edge surfaces
class AtomVecBody *avec;
class BodyRoundedPolyhedron *bptr;
double **discrete; // list of all sub-particles for all bodies
int ndiscrete; // number of discretes in list
int dmax; // allocated size of discrete list
int *dnum; // number of discretes per line, 0 if uninit
int *dfirst; // index of first discrete per each line
int nmax; // allocated size of dnum,dfirst vectors
double **edge; // list of all edge for all bodies
int nedge; // number of edge in list
int edmax; // allocated size of edge list
int *ednum; // number of edges per line, 0 if uninit
int *edfirst; // index of first edge per each line
int ednummax; // allocated size of ednum,edfirst vectors
double **face; // list of all edge for all bodies
int nface; // number of faces in list
int facmax; // allocated size of face list
int *facnum; // number of faces per line, 0 if uninit
int *facfirst; // index of first face per each line
int facnummax; // allocated size of facnum,facfirst vectors
double *enclosing_radius; // enclosing radii for all bodies
double *rounded_radius; // rounded radii for all bodies
double *maxerad; // per-type maximum enclosing radius
void allocate();
void body2space(int);
// sphere-sphere interaction
void sphere_against_sphere(int ibody, int jbody, int itype, int jtype,
double delx, double dely, double delz, double rsq,
double** v, double** f, int evflag);
// sphere-edge interaction
void sphere_against_edge(int ibody, int jbody, int itype, int jtype,
double** x, double** v, double** f, double** torque,
double** angmom, int evflag);
// sphere-face interaction
void sphere_against_face(int ibody, int jbody, int itype, int jtype,
double** x, double** v, double** f, double** torque,
double** angmom, int evflag);
// edge-edge interactions
int edge_against_edge(int ibody, int jbody, int itype, int jtype,
double** x,Contact* contact_list, int &num_contacts,
double &evdwl, double* facc);
// edge-face interactions
int edge_against_face(int ibody, int jbody, int itype, int jtype,
double** x, Contact* contact_list, int &num_contacts,
double &evdwl, double* facc);
// a face vs. a single edge
int interaction_face_to_edge(int ibody, int face_index, double* xmi,
double rounded_radius_i, int jbody, int edge_index,
double* xmj, double rounded_radius_j,
int itype, int jtype, double cut_inner,
Contact* contact_list, int &num_contacts,
double& energy, double* facc);
// an edge vs. an edge from another body
int interaction_edge_to_edge(int ibody, int edge_index_i, double* xmi,
double rounded_radius_i, int jbody, int edge_index_j,
double* xmj, double rounded_radius_j,
int itype, int jtype, double cut_inner,
Contact* contact_list, int &num_contacts,
double& energy, double* facc);
// compute contact forces if contact points are detected
void contact_forces(int ibody, int jbody, double *xi, double *xj,
double delx, double dely, double delz, double fx, double fy, double fz,
double** x, double** v, double** angmom, double** f, double** torque,
double* facc);
// compute force and torque between two bodies given a pair of interacting points
void pair_force_and_torque(int ibody, int jbody, double* pi, double* pj,
double r, double contact_dist, int itype, int jtype,
double** x, double** v, double** f, double** torque,
double** angmom, int jflag, double& energy, double* facc);
// rescale the cohesive forces if a contact area is detected
void rescale_cohesive_forces(double** x, double** f, double** torque,
Contact* contact_list, int &num_contacts,
int itype, int jtype, double* facc);
// compute the separation between two contacts
double contact_separation(const Contact& c1, const Contact& c2);
// detect the unique contact points (as there may be double counts)
void find_unique_contacts(Contact* contact_list, int& num_contacts);
// accumulate torque to a body given a force at a given point
void sum_torque(double* xm, double *x, double fx, double fy, double fz, double* torque);
// find the intersection point (if any) between an edge and a face
int edge_face_intersect(double* x1, double* x2, double* x3, double* a, double* b,
double* hi1, double* hi2, double& d1, double& d2,
int& inside_a, int& inside_b);
// helper functions
int opposite_sides(double* n, double* x0, double* a, double* b);
void project_pt_plane(const double* q, const double* p,
const double* n, double* q_proj, double &d);
void project_pt_plane(const double* q, const double* x1, const double* x2,
const double* x3, double* q_proj, double &d, int& inside);
void project_pt_line(const double* q, const double* xi1, const double* xi2,
double* h, double& d, double& t);
void inside_polygon(int ibody, int face_index, double* xmi,
const double* q1, const double* q2, int& inside1, int& inside2);
void distance_bt_edges(const double* x1, const double* x2,
const double* x3, const double* x4,
double* h1, double* h2, double& t1, double& t2, double& r);
void total_velocity(double* p, double *xcm, double* vcm, double *angmom,
double *inertia, double *quat, double* vi);
void sanity_check();
};
}
#endif
#endif
/* ERROR/WARNING messages:
E: Illegal ... command
Self-explanatory. Check the input script syntax and compare to the
documentation for the command. You can use -echo screen as a
command-line option when running LAMMPS to see the offending line.
E: Incorrect args for pair coefficients
Self-explanatory. Check the input script or data file.
E: Pair body/rounded/polyhedron requires atom style body rounded/polyhedron
Self-explanatory.
E: Pair body requires body style rounded/polyhedron
This pair style is specific to the rounded/polyhedron body style.
*/