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

Merge pull request #3258 from akohlmey/create-atoms-mesh

Browse Source 
Create atoms from STL mesh
This commit is contained in:
Axel Kohlmeyer
2022-05-17 18:11:41 -04:00
committed by GitHub
parent 899764c00c 38d7fcee92
commit 72209b0617
22 changed files with 1409 additions and 73 deletions
Download Patch File Download Diff File Expand all files Collapse all files

3
cmake/Modules/Tools.cmake
View File

@ -3,6 +3,9 @@ if(BUILD_TOOLS)
target_compile_definitions(binary2txt PRIVATE -DLAMMPS_${LAMMPS_SIZES})
install(TARGETS binary2txt DESTINATION ${CMAKE_INSTALL_BINDIR})
add_executable(stl_bin2txt ${LAMMPS_TOOLS_DIR}/stl_bin2txt.cpp)
install(TARGETS stl_bin2txt DESTINATION ${CMAKE_INSTALL_BINDIR})
include(CheckGeneratorSupport)
if(CMAKE_GENERATOR_SUPPORT_FORTRAN)
include(CheckLanguage)

23
doc/src/Tools.rst
View File

@ -56,6 +56,7 @@ Pre-processing tools
* :ref:`moltemplate <moltemplate>`
* :ref:`msi2lmp <msi>`
* :ref:`polybond <polybond>`
* :ref:`stl_bin2txt <stlconvert>`
Post-processing tools
@ -1017,6 +1018,28 @@ For more details please see the README.md file in that folder.
----------
.. _stlconvert:
stl_bin2txt tool
----------------
The file stl_bin2txt.cpp converts binary STL files - like they are
frequently offered for download on the web - into ASCII format STL files
that LAMMPS can read with the :doc:`create_atoms mesh <create_atoms>` or
the :doc:`fix smd/wall_surface <fix_smd_wall_surface>` commands. The syntax
for running the tool is
.. code-block:: bash
stl_bin2txt infile.stl outfile.stl
which creates outfile.stl from infile.stl. This tool must be compiled
on a platform compatible with the byte-ordering that was used to create
the binary file. This usually is a so-called little endian hardware
(like x86).
----------
.. _swig:
SWIG interface

117
doc/src/create_atoms.rst
View File

@ -11,7 +11,7 @@ Syntax
create_atoms type style args keyword values ...
* type = atom type (1-Ntypes) of atoms to create (offset for molecule creation)
* style = *box* or *region* or *single* or *random*
* style = *box* or *region* or *single* or *mesh* or *random*
.. parsed-literal::
@ -20,6 +20,8 @@ Syntax
region-ID = particles will only be created if contained in the region
*single* args = x y z
x,y,z = coordinates of a single particle (distance units)
*mesh* args = STL-file
STL-file = file with triangle mesh in STL format
*random* args = N seed region-ID
N = number of particles to create
seed = random # seed (positive integer)
@ -47,6 +49,14 @@ Syntax
*set* values = dim name
dim = *x* or *y* or *z*
name = name of variable to set with x, y, or z atom position
*radscale* value = factor
factor = scale factor for setting atom radius
*meshmode* values = mode arg
mode = *bisect* or *qrand*
*bisect* arg = radthresh
radthresh = threshold value for *mesh* to determine when to split triangles (distance units)
*qrand* arg = density
density = minimum number density for atoms place on *mesh* triangles (inverse distance squared units)
*rotate* values = theta Rx Ry Rz
theta = rotation angle for single molecule (degrees)
Rx,Ry,Rz = rotation vector for single molecule
@ -69,21 +79,23 @@ Examples
create_atoms 3 single 0 0 5
create_atoms 1 box var v set x xpos set y ypos
create_atoms 2 random 50 12345 NULL overlap 2.0 maxtry 50
create_atoms 1 mesh open_box.stl meshmode qrand 0.1 units box
create_atoms 1 mesh funnel.stl meshmode bisect 4.0 units box radscale 0.9
Description
"""""""""""
This command creates atoms (or molecules) within the simulation box,
either on a lattice, or a single atom (or molecule), or a random
collection of atoms (or molecules). It is an alternative to reading
in atom coordinates explicitly via a :doc:`read_data <read_data>` or
:doc:`read_restart <read_restart>` command. A simulation box must
already exist, which is typically created via the :doc:`create_box
<create_box>` command. Before using this command, a lattice must also
be defined using the :doc:`lattice <lattice>` command, unless you
specify the *single* style with units = box or the *random* style.
For the remainder of this doc page, a created atom or molecule is
referred to as a "particle".
either on a lattice, or a single atom (or molecule), or on a surface
defined by a triangulated mesh, or a random collection of atoms (or
molecules). It is an alternative to reading in atom coordinates
explicitly via a :doc:`read_data <read_data>` or :doc:`read_restart
<read_restart>` command. A simulation box must already exist, which is
typically created via the :doc:`create_box <create_box>` command.
Before using this command, a lattice must also be defined using the
:doc:`lattice <lattice>` command, unless you specify the *single* style
with units = box or the *random* style. For the remainder of this doc
page, a created atom or molecule is referred to as a "particle".
If created particles are individual atoms, they are assigned the
specified atom *type*, though this can be altered via the *basis*
@ -119,6 +131,62 @@ the specified coordinates. This can be useful for debugging purposes
or to create a tiny system with a handful of particles at specified
positions.
.. figure:: img/marble_race.jpg
:figwidth: 33%
:align: right
:target: _images/marble_race.jpg
For the *mesh* style, a file with a triangle mesh in `STL format
<https://en.wikipedia.org/wiki/STL_(file_format)>`_ is read and one or
more particles are placed into the area of each triangle. The reader
supports both ASCII and binary files conforming to the format on the
Wikipedia page. Binary STL files (e.g. as frequently offered for
3d-printing) can also be first converted to ASCII for editing with the
:ref:`stl_bin2txt tool <stlconvert>`. The use of the *units box* option
is required. There are two algorithms available for placing atoms:
*bisect* and *qrand*. They can be selected via the *meshmode* option;
*bisect* is the default. If the atom style allows it, the radius will
be set to a value depending on the algorithm and the value of the
*radscale* parameter (see below), and the atoms created from the mesh
are assigned a new molecule ID.
In *bisect* mode a particle is created at the center of each triangle
unless the average distance of the triangle vertices from its center is
larger than the *radthresh* value (default is lattice spacing in
x-direction). In case the average distance is over the threshold, the
triangle is recursively split into two halves along the the longest side
until the threshold is reached. There will be at least one sphere per
triangle. The value of *radthresh* is set as an argument to *meshmode
bisect*. The average distance of the vertices from the center is also
used to set the radius.
In *qrand* mode a quasi-random sequence is used to distribute particles
on mesh triangles using an approach by :ref:`(Roberts) <Roberts2019>`.
Particles are added to the triangle until the minimum number density is
met or exceeded such that every triangle will have at least one
particle. The minimum number density is set as an argument to the
*qrand* option. The radius will be set so that the sum of the area of
the radius of the particles created in place of a triangle will be equal
to the area of that triangle.
.. note::
The atom placement algorithms in the *mesh* style benefit from meshes
where triangles are close to equilateral. It is therefore
recommended to pre-process STL files to optimize the mesh
accordingly. There are multiple open source and commercial software
tools available with the capability to generate optimized meshes.
.. note::
In most cases the atoms created in *mesh* style will become an
immobile or rigid object that would not be time integrated or moved
by :doc:`fix move <fix_move>` or :doc:`fix rigid <fix_rigid>`. For
computational efficiency *and* to avoid undesired contributions to
pressure and potential energy due to close contacts, it is usually
beneficial to exclude computing interactions between the created
particles using :doc:`neigh_modify exclude <neigh_modify>`.
For the *random* style, *N* particles are added to the system at
randomly generated coordinates, which can be useful for generating an
amorphous system. The particles are created one by one using the
@ -316,6 +384,13 @@ the atoms around the rotation axis is consistent with the right-hand
rule: if your right-hand's thumb points along *R*, then your fingers
wrap around the axis in the direction of rotation.
The *radscale* keyword only applies to the *mesh* style and adjusts the
radius of created particles (see above), provided this is supported by
the atom style. Its value is a prefactor (must be > 0.0, default is
1.0) that is applied to the atom radius inferred from the size of the
individual triangles in the triangle mesh that the particle corresponds
to.
The *overlap* keyword only applies to the *random* style. It prevents
newly created particles from being created closer than the specified
*Doverlap* distance from any other particle. When the particles being
@ -424,9 +499,11 @@ values specified in the file read by the :doc:`molecule <molecule>`
command. E.g. the file typically defines bonds (angles,etc) between
atoms in the molecule, and can optionally define charges on each atom.
Note that the *sphere* atom style sets the default particle diameter
to 1.0 as well as the density. This means the mass for the particle
is not 1.0, but is PI/6 \* diameter\^3 = 0.5236.
Note that the *sphere* atom style sets the default particle diameter to
1.0 as well as the density. This means the mass for the particle is not
1.0, but is PI/6 \* diameter\^3 = 0.5236. When using the *mesh* style,
the particle diameter is adjusted from the size of the individual
triangles in the triangle mesh.
Note that the *ellipsoid* atom style sets the default particle shape
to (0.0 0.0 0.0) and the density to 1.0 which means it is a point
@ -460,5 +537,13 @@ Default
The default for the *basis* keyword is that all created atoms are
assigned the argument *type* as their atom type (when single atoms are
being created). The other defaults are *remap* = no, *rotate* =
random, *overlap* not checked, *maxtry* = 10, and *units* = lattice.
being created). The other defaults are *remap* = no, *rotate* = random,
*radscale* = 1.0, *radthresh* = x-lattice spacing, *overlap* not
checked, *maxtry* = 10, and *units* = lattice.
----------
.. _Roberts2019:
**(Roberts)** R. Roberts (2019) "Evenly Distributing Points in a Triangle." Extreme Learning.
`<http://extremelearning.com.au/evenly-distributing-points-in-a-triangle/>`_

2
doc/src/dump_image.rst
View File

@ -44,7 +44,7 @@ Syntax
color = *type*
bflag1,bflag2 = 2 numeric flags to affect how bodies are drawn
*fix* = fixID color fflag1 fflag2
fixID = ID of fix that generates objects to dray
fixID = ID of fix that generates objects to draw
color = *type*
fflag1,fflag2 = 2 numeric flags to affect how fix objects are drawn
*size* values = width height = size of images

BIN
doc/src/img/marble_race.jpg Normal file
View File

Binary file not shown.
Side by Side

After

Width:  |  Height:  |  Size: 39 KiB

1
doc/utils/sphinx-config/false_positives.txt
View File

@ -3216,6 +3216,7 @@ Stesmans
stiffnesses
Stillinger
stk
stl
stochastically
stochasticity
Stockmayer

1
examples/README
View File

@ -86,6 +86,7 @@ mc: MC package models: GCMC, Widom, fix mol/swap
mdi: use of the MDI package and MolSSI MDI code coupling library
meam: MEAM test for SiC and shear (same as shear examples)
melt: rapid melt of 3d LJ system
mesh: create_atoms mesh command
micelle: self-assembly of small lipid-like molecules into 2d bilayers
min: energy minimization of 2d LJ melt
mliap: examples for using several bundled MLIAP potentials

43
examples/mesh/in.marble_race Normal file
View File

@ -0,0 +1,43 @@
units real
lattice sc 5.0
region box block -110 60 -30 220 -90 130 units box
create_box 2 box
region particles cylinder y 0 -30 47 130 190 units box
create_atoms 1 region particles
region lid cylinder y 0 -30 47 190 200 units box
group mobile type 1
create_atoms 2 mesh race_track.stl units box
group mesh type 2
mass * 39.95
pair_style lj/cut 8.76
pair_coeff 1 1 0.2339 3.504
pair_coeff 1 2 0.2339 7.008 $(7.008*2^(1.0/6.0))
pair_coeff 2 2 0.0 1.0
balance 1.1 shift xyz 10 1.01
neigh_modify exclude type 2 2
timestep 1.0
fix track mesh setforce 0.0 0.0 0.0
fix pull mobile addforce 0.0 -0.05 0.0 region particles
fix dir mobile oneway 10 lid -y
fix move mobile nve
fix load all balance 1000 1.1 shift xyz 10 1.01 weight neigh 0.5 weight group 2 mesh 0.1 mobile 1.0
minimize 0.0 0.0 1000 1000
reset_timestep 0 time 0.0
velocity mobile create 150.0 54634234
compute ptemp mobile temp
thermo_modify temp ptemp
thermo 1000
# dump 1 all atom 1000 race.lammpstrj
run 10000

49
examples/mesh/in.mesh_box Normal file
View File

@ -0,0 +1,49 @@
units real
atom_style hybrid sphere bond
lattice sc 5.0
region box block 50 250 50 250 50 250 units box
create_box 2 box
region particles block 110 190 110 190 110 190 units box
create_atoms 1 region particles
region lid block 100 110 50 250 50 250 units box
group mobile type 1
set type 1 diameter 7.0
# create_atoms 2 mesh open_box.stl meshmode bisect 4.0 units box
create_atoms 2 mesh open_box.stl meshmode qrand 0.1 units box
group mesh type 2
pair_style lj/cut 8.76
pair_coeff 1 1 0.2339 3.504
pair_coeff 1 2 0.2339 7.008 $(7.008*2^(1.0/6.0))
pair_coeff 2 2 0.0 1.0
mass * 39.95
neigh_modify exclude type 2 2
timestep 1.0
run 0 post no
fix dir mobile oneway 10 lid x
fix move mobile nve
fix load all balance 1000 1.1 shift xyz 10 1.01 weight neigh 0.5 weight group 2 mesh 0.1 mobile 1.0
fix rot mesh move rotate 150.0 150.0 150.0 1.0 0.0 0.0 500000.0 units box
reset_timestep 0 time 0.0
velocity mobile create 150.0 54634234
compute ptemp mobile temp
thermo_modify temp ptemp
thermo 200
compute ke all ke/atom
#dump 2 all movie 200 mesh.mkv c_ke radius size 960 1440 zoom 1.5 box no 0.0 view 120 180
#dump_modify 2 bitrate 4000 framerate 12 color orange 1.0 0.5 0.0 amap min max cf 0.0 6 min blue 0.1 fuchsia 0.2 red 0.4 orange 0.6 yellow max white
#dump 1 all custom 500 open_box.lammpstrj id type mol x y z vx vy vz
run 5000

156
examples/mesh/log.11May22.marble_race.g++.1 Normal file
View File

@ -0,0 +1,156 @@
LAMMPS (4 May 2022)
OMP_NUM_THREADS environment is not set. Defaulting to 1 thread. (src/comm.cpp:98)
using 1 OpenMP thread(s) per MPI task
units real
lattice sc 5.0
Lattice spacing in x,y,z = 5 5 5
region box block -110 60 -30 220 -90 130 units box
create_box 2 box
Created orthogonal box = (-110 -30 -90) to (60 220 130)
1 by 1 by 1 MPI processor grid
region particles cylinder y 0 -30 47 130 190 units box
create_atoms 1 region particles
Created 3601 atoms
using lattice units in orthogonal box = (-110 -30 -90) to (60 220 130)
create_atoms CPU = 0.001 seconds
region lid cylinder y 0 -30 47 190 200 units box
group mobile type 1
3601 atoms in group mobile
create_atoms 2 mesh race_track.stl units box
Reading STL object Georgs Murmelbahn from file race_track.stl
read 9472 triangles with 1.00 atoms per triangle added in recursive bisection mode
Created 9472 atoms
using box units in orthogonal box = (-110 -30 -90) to (60 220 130)
create_atoms CPU = 0.040 seconds
group mesh type 2
9472 atoms in group mesh
mass * 39.95
pair_style lj/cut 8.76
pair_coeff 1 1 0.2339 3.504
pair_coeff 1 2 0.2339 7.008 $(7.008*2^(1.0/6.0))
pair_coeff 1 2 0.2339 7.008 7.8662140345520858986
pair_coeff 2 2 0.0 1.0
balance 1.1 shift xyz 10 1.01
Balancing ...
Generated 0 of 1 mixed pair_coeff terms from geometric mixing rule
Neighbor list info ...
update every 1 steps, delay 10 steps, check yes
max neighbors/atom: 2000, page size: 100000
master list distance cutoff = 10.76
ghost atom cutoff = 10.76
binsize = 5.38, bins = 32 47 41
1 neighbor lists, perpetual/occasional/extra = 1 0 0
(1) pair lj/cut, perpetual
attributes: half, newton on
pair build: half/bin/atomonly/newton
stencil: half/bin/3d
bin: standard
neigh_modify exclude type 2 2
timestep 1.0
fix track mesh setforce 0.0 0.0 0.0
fix pull mobile addforce 0.0 -0.05 0.0 region particles
fix dir mobile oneway 10 lid -y
fix move mobile nve
fix load all balance 1000 1.1 shift xyz 10 1.01 weight neigh 0.5 weight group 2 mesh 0.1 mobile 1.0
minimize 0.0 0.0 1000 1000
WARNING: Using 'neigh_modify every 1 delay 0 check yes' setting during minimization (src/min.cpp:187)
Generated 0 of 1 mixed pair_coeff terms from geometric mixing rule
Per MPI rank memory allocation (min/avg/max) = 4.803 | 4.803 | 4.803 Mbytes
Step Temp E_pair E_mol TotEng Press Volume
0 0 737062.81 0 737062.81 21986.781 9350000
67 0 -2063.91 0 -2063.91 -5.0227698 9350000
Loop time of 0.517938 on 1 procs for 67 steps with 13073 atoms
99.8% CPU use with 1 MPI tasks x 1 OpenMP threads
Minimization stats:
Stopping criterion = linesearch alpha is zero
Energy initial, next-to-last, final =
737062.806250078 -2063.90998808136 -2063.90998808136
Force two-norm initial, final = 689296.27 22.226599
Force max component initial, final = 336546.89 0.90593277
Final line search alpha, max atom move = 1.2850327e-11 1.1641532e-11
Iterations, force evaluations = 67 393
MPI task timing breakdown:
Section | min time | avg time | max time |%varavg| %total
---------------------------------------------------------------
Pair | 0.38475 | 0.38475 | 0.38475 | 0.0 | 74.28
Neigh | 0.037593 | 0.037593 | 0.037593 | 0.0 | 7.26
Comm | 0.0015299 | 0.0015299 | 0.0015299 | 0.0 | 0.30
Output | 0 | 0 | 0 | 0.0 | 0.00
Modify | 0.032572 | 0.032572 | 0.032572 | 0.0 | 6.29
Other | | 0.0615 | | | 11.87
Nlocal: 13073 ave 13073 max 13073 min
Histogram: 1 0 0 0 0 0 0 0 0 0
Nghost: 590 ave 590 max 590 min
Histogram: 1 0 0 0 0 0 0 0 0 0
Neighs: 66778 ave 66778 max 66778 min
Histogram: 1 0 0 0 0 0 0 0 0 0
Total # of neighbors = 66778
Ave neighs/atom = 5.1080854
Neighbor list builds = 5
Dangerous builds = 0
reset_timestep 0 time 0.0
velocity mobile create 150.0 54634234
compute ptemp mobile temp
thermo_modify temp ptemp
WARNING: Temperature for thermo pressure is not for group all (src/thermo.cpp:520)
thermo 1000
# dump 1 all atom 1000 race.lammpstrj
run 10000
Generated 0 of 1 mixed pair_coeff terms from geometric mixing rule
Per MPI rank memory allocation (min/avg/max) = 3.678 | 3.678 | 3.678 Mbytes
Step Temp E_pair E_mol TotEng Press Volume
0 150 -2063.91 0 -454.27256 2.8467694 9350000
1000 197.01825 -1960.8576 0 153.32871 13.788868 9350000
2000 216.32291 -2037.8959 0 283.44712 13.915645 9350000
3000 239.06547 -2059.8437 0 505.54782 14.550975 9350000
4000 266.60476 -2076.3375 0 784.57583 17.457495 9350000
5000 299.6816 -2109.6562 0 1106.202 19.375766 9350000
6000 335.17037 -2129.5487 0 1467.1364 24.636144 9350000
7000 367.9265 -2101.7855 0 1846.4029 33.591291 9350000
8000 404.7304 -2144.8541 0 2198.2739 43.134333 9350000
9000 435.75368 -2180.3183 0 2495.7179 53.466409 9350000
10000 457.96804 -2194.6681 0 2719.7486 64.521177 9350000
Loop time of 12.1711 on 1 procs for 10000 steps with 13073 atoms
Performance: 70.988 ns/day, 0.338 hours/ns, 821.617 timesteps/s
99.9% CPU use with 1 MPI tasks x 1 OpenMP threads
MPI task timing breakdown:
Section | min time | avg time | max time |%varavg| %total
---------------------------------------------------------------
Pair | 9.4572 | 9.4572 | 9.4572 | 0.0 | 77.70
Neigh | 1.1196 | 1.1196 | 1.1196 | 0.0 | 9.20
Comm | 0.04293 | 0.04293 | 0.04293 | 0.0 | 0.35
Output | 0.00042467 | 0.00042467 | 0.00042467 | 0.0 | 0.00
Modify | 1.241 | 1.241 | 1.241 | 0.0 | 10.20
Other | | 0.31 | | | 2.55
Nlocal: 13073 ave 13073 max 13073 min
Histogram: 1 0 0 0 0 0 0 0 0 0
Nghost: 590 ave 590 max 590 min
Histogram: 1 0 0 0 0 0 0 0 0 0
Neighs: 112709 ave 112709 max 112709 min
Histogram: 1 0 0 0 0 0 0 0 0 0
Total # of neighbors = 112709
Ave neighs/atom = 8.62151
Neighbor list builds = 129
Dangerous builds = 0
Total wall time: 0:00:12

163
examples/mesh/log.11May22.marble_race.g++.4 Normal file
View File

@ -0,0 +1,163 @@
LAMMPS (4 May 2022)
OMP_NUM_THREADS environment is not set. Defaulting to 1 thread. (src/comm.cpp:98)
using 1 OpenMP thread(s) per MPI task
units real
lattice sc 5.0
Lattice spacing in x,y,z = 5 5 5
region box block -110 60 -30 220 -90 130 units box
create_box 2 box
Created orthogonal box = (-110 -30 -90) to (60 220 130)
1 by 2 by 2 MPI processor grid
region particles cylinder y 0 -30 47 130 190 units box
create_atoms 1 region particles
Created 3601 atoms
using lattice units in orthogonal box = (-110 -30 -90) to (60 220 130)
create_atoms CPU = 0.001 seconds
region lid cylinder y 0 -30 47 190 200 units box
group mobile type 1
3601 atoms in group mobile
create_atoms 2 mesh race_track.stl units box
Reading STL object Georgs Murmelbahn from file race_track.stl
read 9472 triangles with 1.00 atoms per triangle added in recursive bisection mode
Created 9472 atoms
using box units in orthogonal box = (-110 -30 -90) to (60 220 130)
create_atoms CPU = 0.036 seconds
group mesh type 2
9472 atoms in group mesh
mass * 39.95
pair_style lj/cut 8.76
pair_coeff 1 1 0.2339 3.504
pair_coeff 1 2 0.2339 7.008 $(7.008*2^(1.0/6.0))
pair_coeff 1 2 0.2339 7.008 7.8662140345520858986
pair_coeff 2 2 0.0 1.0
balance 1.1 shift xyz 10 1.01
Balancing ...
Generated 0 of 1 mixed pair_coeff terms from geometric mixing rule
Neighbor list info ...
update every 1 steps, delay 10 steps, check yes
max neighbors/atom: 2000, page size: 100000
master list distance cutoff = 10.76
ghost atom cutoff = 10.76
binsize = 5.38, bins = 32 47 41
1 neighbor lists, perpetual/occasional/extra = 1 0 0
(1) pair lj/cut, perpetual
attributes: half, newton on
pair build: half/bin/atomonly/newton
stencil: half/bin/3d
bin: standard
rebalancing time: 0.001 seconds
iteration count = 17
initial/final maximal load/proc = 6175 4316
initial/final imbalance factor = 1.8893903 1.3205844
x cuts: 0 1
y cuts: 0 0.5859375 1
z cuts: 0 0.36376953 1
neigh_modify exclude type 2 2
timestep 1.0
fix track mesh setforce 0.0 0.0 0.0
fix pull mobile addforce 0.0 -0.05 0.0 region particles
fix dir mobile oneway 10 lid -y
fix move mobile nve
fix load all balance 1000 1.1 shift xyz 10 1.01 weight neigh 0.5 weight group 2 mesh 0.1 mobile 1.0
minimize 0.0 0.0 1000 1000
WARNING: Using 'neigh_modify every 1 delay 0 check yes' setting during minimization (src/min.cpp:187)
Generated 0 of 1 mixed pair_coeff terms from geometric mixing rule
Per MPI rank memory allocation (min/avg/max) = 4.506 | 4.54 | 4.57 Mbytes
Step Temp E_pair E_mol TotEng Press Volume
0 0 737062.81 0 737062.81 21986.781 9350000
67 0 -2063.91 0 -2063.91 -5.0227698 9350000
Loop time of 0.373091 on 4 procs for 67 steps with 13073 atoms
98.7% CPU use with 4 MPI tasks x 1 OpenMP threads
Minimization stats:
Stopping criterion = linesearch alpha is zero
Energy initial, next-to-last, final =
737062.806250145 -2063.90998808134 -2063.90998808134
Force two-norm initial, final = 689296.27 22.226599
Force max component initial, final = 336546.89 0.90593277
Final line search alpha, max atom move = 1.2850327e-11 1.1641532e-11
Iterations, force evaluations = 67 393
MPI task timing breakdown:
Section | min time | avg time | max time |%varavg| %total
---------------------------------------------------------------
Pair | 0.0023911 | 0.098676 | 0.3067 | 39.6 | 26.45
Neigh | 0.0054604 | 0.0097001 | 0.017713 | 4.9 | 2.60
Comm | 0.0072159 | 0.22319 | 0.32344 | 27.2 | 59.82
Output | 0 | 0 | 0 | 0.0 | 0.00
Modify | 0.003044 | 0.0084861 | 0.018696 | 6.5 | 2.27
Other | | 0.03304 | | | 8.86
Nlocal: 3268.25 ave 4314 max 2222 min
Histogram: 2 0 0 0 0 0 0 0 0 2
Nghost: 1010.75 ave 1101 max 792 min
Histogram: 1 0 0 0 0 0 0 0 1 2
Neighs: 16694.5 ave 52643 max 0 min
Histogram: 2 0 1 0 0 0 0 0 0 1
Total # of neighbors = 66778
Ave neighs/atom = 5.1080854
Neighbor list builds = 5
Dangerous builds = 0
reset_timestep 0 time 0.0
velocity mobile create 150.0 54634234
compute ptemp mobile temp
thermo_modify temp ptemp
WARNING: Temperature for thermo pressure is not for group all (src/thermo.cpp:520)
thermo 1000
# dump 1 all atom 1000 race.lammpstrj
run 10000
Generated 0 of 1 mixed pair_coeff terms from geometric mixing rule
Per MPI rank memory allocation (min/avg/max) = 3.481 | 3.598 | 3.711 Mbytes
Step Temp E_pair E_mol TotEng Press Volume
0 150 -2063.91 0 -454.27256 2.8467694 9350000
1000 197.01825 -1960.8576 0 153.32871 13.788868 9350000
2000 216.32291 -2037.8959 0 283.44712 13.915645 9350000
3000 239.06547 -2059.8437 0 505.54782 14.550975 9350000
4000 266.60476 -2076.3375 0 784.57583 17.457495 9350000
5000 299.6816 -2109.6562 0 1106.202 19.375766 9350000
6000 335.17037 -2129.5487 0 1467.1364 24.636144 9350000
7000 367.9265 -2101.7855 0 1846.4029 33.591291 9350000
8000 404.7304 -2144.8541 0 2198.2739 43.134332 9350000
9000 435.7537 -2180.3187 0 2495.7178 53.466393 9350000
10000 457.96586 -2194.6411 0 2719.7522 64.522003 9350000
Loop time of 4.60636 on 4 procs for 10000 steps with 13073 atoms
Performance: 187.567 ns/day, 0.128 hours/ns, 2170.909 timesteps/s
99.2% CPU use with 4 MPI tasks x 1 OpenMP threads
MPI task timing breakdown:
Section | min time | avg time | max time |%varavg| %total
---------------------------------------------------------------
Pair | 2.0632 | 2.4782 | 2.963 | 25.3 | 53.80
Neigh | 0.18255 | 0.30042 | 0.52984 | 25.5 | 6.52
Comm | 0.40668 | 1.1595 | 1.7224 | 51.6 | 25.17
Output | 0.00032705 | 0.00053493 | 0.00067573 | 0.0 | 0.01
Modify | 0.22563 | 0.32471 | 0.53003 | 21.1 | 7.05
Other | | 0.3429 | | | 7.44
Nlocal: 3268.25 ave 6890 max 1643 min
Histogram: 1 2 0 0 0 0 0 0 0 1
Nghost: 1701 ave 2074 max 1456 min
Histogram: 2 0 0 0 0 1 0 0 0 1
Neighs: 28177.5 ave 34170 max 21435 min
Histogram: 1 1 0 0 0 0 0 0 0 2
Total # of neighbors = 112710
Ave neighs/atom = 8.6215865
Neighbor list builds = 129
Dangerous builds = 0
Total wall time: 0:00:05

147
examples/mesh/log.11May22.mesh_box.g++.1 Normal file
View File

@ -0,0 +1,147 @@
LAMMPS (4 May 2022)
OMP_NUM_THREADS environment is not set. Defaulting to 1 thread. (src/comm.cpp:98)
using 1 OpenMP thread(s) per MPI task
units real
atom_style hybrid sphere bond
WARNING: Atom style hybrid defines both, per-type and per-atom masses; both must be set, but only per-atom masses will be used (src/atom_vec_hybrid.cpp:130)
lattice sc 5.0
Lattice spacing in x,y,z = 5 5 5
region box block 50 250 50 250 50 250 units box
create_box 2 box
Created orthogonal box = (50 50 50) to (250 250 250)
1 by 1 by 1 MPI processor grid
region particles block 110 190 110 190 110 190 units box
create_atoms 1 region particles
Created 4913 atoms
using lattice units in orthogonal box = (50 50 50) to (250 250 250)
create_atoms CPU = 0.001 seconds
region lid block 100 110 50 250 50 250 units box
group mobile type 1
4913 atoms in group mobile
set type 1 diameter 7.0
Setting atom values ...
4913 settings made for diameter
# create_atoms 2 mesh open_box.stl meshmode bisect 4.0 units box
create_atoms 2 mesh open_box.stl meshmode qrand 0.1 units box
Reading STL object Open Box from file open_box.stl
read 10 triangles with 500.00 atoms per triangle added in quasi-random mode
Created 5000 atoms
using box units in orthogonal box = (50 50 50) to (250 250 250)
create_atoms CPU = 0.001 seconds
group mesh type 2
5000 atoms in group mesh
pair_style lj/cut 8.76
pair_coeff 1 1 0.2339 3.504
pair_coeff 1 2 0.2339 7.008 $(7.008*2^(1.0/6.0))
pair_coeff 1 2 0.2339 7.008 7.8662140345520858986
pair_coeff 2 2 0.0 1.0
mass * 39.95
neigh_modify exclude type 2 2
timestep 1.0
run 0 post no
WARNING: No fixes with time integration, atoms won't move (src/verlet.cpp:60)
Generated 0 of 1 mixed pair_coeff terms from geometric mixing rule
Neighbor list info ...
update every 1 steps, delay 10 steps, check yes
max neighbors/atom: 2000, page size: 100000
master list distance cutoff = 10.76
ghost atom cutoff = 10.76
binsize = 5.38, bins = 38 38 38
1 neighbor lists, perpetual/occasional/extra = 1 0 0
(1) pair lj/cut, perpetual
attributes: half, newton on
pair build: half/bin/newton
stencil: half/bin/3d
bin: standard
Per MPI rank memory allocation (min/avg/max) = 5.64 | 5.64 | 5.64 Mbytes
Step Temp E_pair E_mol TotEng Press
0 0 -1778.6527 0 -1778.6527 -27.271044
Loop time of 1.514e-06 on 1 procs for 0 steps with 9913 atoms
fix dir mobile oneway 10 lid x
fix move mobile nve
fix load all balance 1000 1.1 shift xyz 10 1.01 weight neigh 0.5 weight group 2 mesh 0.1 mobile 1.0
fix rot mesh move rotate 150.0 150.0 150.0 1.0 0.0 0.0 500000.0 units box
reset_timestep 0 time 0.0
velocity mobile create 150.0 54634234
compute ptemp mobile temp
thermo_modify temp ptemp
WARNING: Temperature for thermo pressure is not for group all (src/thermo.cpp:520)
thermo 200
compute ke all ke/atom
#dump 2 all movie 200 mesh.mkv c_ke radius size 960 1440 zoom 1.5 box no 0.0 view 120 180
#dump_modify 2 bitrate 4000 framerate 12 color orange 1.0 0.5 0.0 amap min max cf 0.0 6 min blue 0.1 fuchsia 0.2 red 0.4 orange 0.6 yellow max white
#dump 1 all custom 500 open_box.lammpstrj id type mol x y z vx vy vz
run 5000
Generated 0 of 1 mixed pair_coeff terms from geometric mixing rule
Per MPI rank memory allocation (min/avg/max) = 6.155 | 6.155 | 6.155 Mbytes
Step Temp E_pair E_mol TotEng Press Volume
0 150 -1778.6527 0 417.60814 -14.721534 8000000
200 179.67539 -2344.2923 0 286.46789 14.546173 8000000
400 181.6261 -2331.1731 0 328.14895 15.522647 8000000
600 183.6116 -2362.5284 0 325.86464 14.828678 8000000
800 187.69451 -2422.3021 0 325.87189 14.99914 8000000
1000 186.32256 -2393.187 0 334.89931 14.49224 8000000
1200 186.57797 -2406.6141 0 325.21184 15.514285 8000000
1400 188.39717 -2403.2075 0 355.25478 14.282381 8000000
1600 185.774 -2391.15 0 328.90446 16.025507 8000000
1800 186.95292 -2423.2897 0 314.02613 15.51055 8000000
2000 186.0637 -2414.6095 0 309.68668 14.894421 8000000
2200 187.03262 -2409.2253 0 329.25754 14.90094 8000000
2400 186.66862 -2409.987 0 323.16626 15.497536 8000000
2600 184.98646 -2386.3118 0 322.21161 14.935837 8000000
2800 185.83051 -2392.5258 0 328.35604 14.863905 8000000
3000 184.68382 -2397.8429 0 306.24943 15.696458 8000000
3200 187.56974 -2409.8941 0 336.45315 14.352166 8000000
3400 187.5721 -2425.0188 0 321.36292 14.7297 8000000
3600 185.97304 -2391.4399 0 331.52886 15.586758 8000000
3800 185.40034 -2401.6336 0 312.94973 15.742308 8000000
4000 187.71377 -2409.3588 0 339.09729 15.102297 8000000
4200 186.4676 -2394.1921 0 336.01789 15.312368 8000000
4400 186.98262 -2396.3842 0 341.36649 14.764489 8000000
4600 185.98808 -2397.7342 0 325.45468 15.379472 8000000
4800 187.6927 -2422.0727 0 326.07474 14.252141 8000000
5000 188.21075 -2428.1325 0 327.60023 14.694135 8000000
Loop time of 9.49245 on 1 procs for 5000 steps with 9913 atoms
Performance: 45.510 ns/day, 0.527 hours/ns, 526.734 timesteps/s
99.8% CPU use with 1 MPI tasks x 1 OpenMP threads
MPI task timing breakdown:
Section | min time | avg time | max time |%varavg| %total
---------------------------------------------------------------
Pair | 5.3698 | 5.3698 | 5.3698 | 0.0 | 56.57
Bond | 0.00036337 | 0.00036337 | 0.00036337 | 0.0 | 0.00
Neigh | 3.0911 | 3.0911 | 3.0911 | 0.0 | 32.56
Comm | 0.039329 | 0.039329 | 0.039329 | 0.0 | 0.41
Output | 0.0010197 | 0.0010197 | 0.0010197 | 0.0 | 0.01
Modify | 0.86293 | 0.86293 | 0.86293 | 0.0 | 9.09
Other | | 0.1279 | | | 1.35
Nlocal: 9913 ave 9913 max 9913 min
Histogram: 1 0 0 0 0 0 0 0 0 0
Nghost: 0 ave 0 max 0 min
Histogram: 1 0 0 0 0 0 0 0 0 0
Neighs: 97756 ave 97756 max 97756 min
Histogram: 1 0 0 0 0 0 0 0 0 0
Total # of neighbors = 97756
Ave neighs/atom = 9.8613941
Ave special neighs/atom = 0
Neighbor list builds = 375
Dangerous builds = 1
Total wall time: 0:00:09

147
examples/mesh/log.11May22.mesh_box.g++.4 Normal file
View File

@ -0,0 +1,147 @@
LAMMPS (4 May 2022)
OMP_NUM_THREADS environment is not set. Defaulting to 1 thread. (src/comm.cpp:98)
using 1 OpenMP thread(s) per MPI task
units real
atom_style hybrid sphere bond
WARNING: Atom style hybrid defines both, per-type and per-atom masses; both must be set, but only per-atom masses will be used (src/atom_vec_hybrid.cpp:130)
lattice sc 5.0
Lattice spacing in x,y,z = 5 5 5
region box block 50 250 50 250 50 250 units box
create_box 2 box
Created orthogonal box = (50 50 50) to (250 250 250)
1 by 2 by 2 MPI processor grid
region particles block 110 190 110 190 110 190 units box
create_atoms 1 region particles
Created 4913 atoms
using lattice units in orthogonal box = (50 50 50) to (250 250 250)
create_atoms CPU = 0.001 seconds
region lid block 100 110 50 250 50 250 units box
group mobile type 1
4913 atoms in group mobile
set type 1 diameter 7.0
Setting atom values ...
4913 settings made for diameter
# create_atoms 2 mesh open_box.stl meshmode bisect 4.0 units box
create_atoms 2 mesh open_box.stl meshmode qrand 0.1 units box
Reading STL object Open Box from file open_box.stl
read 10 triangles with 2000.00 atoms per triangle added in quasi-random mode
Created 5000 atoms
using box units in orthogonal box = (50 50 50) to (250 250 250)
create_atoms CPU = 0.001 seconds
group mesh type 2
5000 atoms in group mesh
pair_style lj/cut 8.76
pair_coeff 1 1 0.2339 3.504
pair_coeff 1 2 0.2339 7.008 $(7.008*2^(1.0/6.0))
pair_coeff 1 2 0.2339 7.008 7.8662140345520858986
pair_coeff 2 2 0.0 1.0
mass * 39.95
neigh_modify exclude type 2 2
timestep 1.0
run 0 post no
WARNING: No fixes with time integration, atoms won't move (src/verlet.cpp:60)
Generated 0 of 1 mixed pair_coeff terms from geometric mixing rule
Neighbor list info ...
update every 1 steps, delay 10 steps, check yes
max neighbors/atom: 2000, page size: 100000
master list distance cutoff = 10.76
ghost atom cutoff = 10.76
binsize = 5.38, bins = 38 38 38
1 neighbor lists, perpetual/occasional/extra = 1 0 0
(1) pair lj/cut, perpetual
attributes: half, newton on
pair build: half/bin/newton
stencil: half/bin/3d
bin: standard
Per MPI rank memory allocation (min/avg/max) = 5.455 | 5.474 | 5.491 Mbytes
Step Temp E_pair E_mol TotEng Press
0 0 -1778.6527 0 -1778.6527 -27.271044
Loop time of 2.8725e-06 on 4 procs for 0 steps with 9913 atoms
fix dir mobile oneway 10 lid x
fix move mobile nve
fix load all balance 1000 1.1 shift xyz 10 1.01 weight neigh 0.5 weight group 2 mesh 0.1 mobile 1.0
fix rot mesh move rotate 150.0 150.0 150.0 1.0 0.0 0.0 500000.0 units box
reset_timestep 0 time 0.0
velocity mobile create 150.0 54634234
compute ptemp mobile temp
thermo_modify temp ptemp
WARNING: Temperature for thermo pressure is not for group all (src/thermo.cpp:520)
thermo 200
compute ke all ke/atom
#dump 2 all movie 200 mesh.mkv c_ke radius size 960 1440 zoom 1.5 box no 0.0 view 120 180
#dump_modify 2 bitrate 4000 framerate 12 color orange 1.0 0.5 0.0 amap min max cf 0.0 6 min blue 0.1 fuchsia 0.2 red 0.4 orange 0.6 yellow max white
#dump 1 all custom 500 open_box.lammpstrj id type mol x y z vx vy vz
run 5000
Generated 0 of 1 mixed pair_coeff terms from geometric mixing rule
Per MPI rank memory allocation (min/avg/max) = 6.011 | 6.046 | 6.072 Mbytes
Step Temp E_pair E_mol TotEng Press Volume
0 150 -1778.6527 0 417.60814 -14.721534 8000000
200 178.90737 -2298.4303 0 321.08483 15.681113 8000000
400 182.18748 -2339.8994 0 327.64208 13.983929 8000000
600 182.78421 -2363.0316 0 313.24702 15.533691 8000000
800 185.27077 -2397.561 0 315.12532 14.162041 8000000
1000 185.59249 -2394.3379 0 323.05898 14.45888 8000000
1200 184.9041 -2365.2474 0 342.07022 14.281891 8000000
1400 186.59058 -2397.5471 0 334.46342 15.095686 8000000
1600 184.44411 -2386.6701 0 313.91242 16.825132 8000000
1800 188.29985 -2432.366 0 324.6713 15.109591 8000000
2000 187.05565 -2423.8405 0 314.97946 15.67693 8000000
2200 187.02332 -2411.5557 0 326.79085 14.702305 8000000
2400 185.67262 -2386.3932 0 332.17677 15.796702 8000000
2600 185.97931 -2411.0423 0 312.01823 16.15444 8000000
2800 187.2186 -2403.7584 0 337.44748 15.197009 8000000
3000 185.18778 -2388.9601 0 322.51107 16.289848 8000000
3200 186.38078 -2396.9152 0 332.02344 15.860688 8000000
3400 184.95984 -2386.3632 0 321.77041 16.513329 8000000
3600 187.94957 -2408.3205 0 343.58806 13.931044 8000000
3800 187.1454 -2413.0585 0 327.07559 14.801401 8000000
4000 186.12517 -2394.4218 0 330.77436 14.257144 8000000
4200 186.91955 -2407.2113 0 329.61591 14.869086 8000000
4400 184.72142 -2382.9411 0 321.70176 15.898793 8000000
4600 186.8458 -2408.8059 0 326.94154 14.766765 8000000
4800 187.28977 -2397.6811 0 344.56676 15.301534 8000000
5000 185.17839 -2378.709 0 332.62466 16.043707 8000000
Loop time of 3.41925 on 4 procs for 5000 steps with 9913 atoms
Performance: 126.344 ns/day, 0.190 hours/ns, 1462.309 timesteps/s
99.3% CPU use with 4 MPI tasks x 1 OpenMP threads
MPI task timing breakdown:
Section | min time | avg time | max time |%varavg| %total
---------------------------------------------------------------
Pair | 1.011 | 1.4158 | 1.8652 | 31.9 | 41.41
Bond | 0.00026998 | 0.00029559 | 0.00031513 | 0.0 | 0.01
Neigh | 0.60535 | 0.82114 | 1.0415 | 21.6 | 24.02
Comm | 0.19046 | 0.88593 | 1.5316 | 63.7 | 25.91
Output | 0.00067893 | 0.00085803 | 0.001033 | 0.0 | 0.03
Modify | 0.20164 | 0.23804 | 0.27649 | 6.5 | 6.96
Other | | 0.05722 | | | 1.67
Nlocal: 2478.25 ave 2954 max 2037 min
Histogram: 1 1 0 0 0 0 0 0 1 1
Nghost: 1043.25 ave 1146 max 940 min
Histogram: 2 0 0 0 0 0 0 0 0 2
Neighs: 24605.8 ave 32065 max 17721 min
Histogram: 2 0 0 0 0 0 0 0 1 1
Total # of neighbors = 98423
Ave neighs/atom = 9.9286795
Ave special neighs/atom = 0
Neighbor list builds = 371
Dangerous builds = 0
Total wall time: 0:00:03

72
examples/mesh/open_box.stl Normal file
View File

@ -0,0 +1,72 @@
solid Open Box
facet normal 0.0 0.0 1.0
outer loop
vertex 100.0 100.0 100.0
vertex 200.0 100.0 100.0
vertex 200.0 200.0 100.0
endloop
endfacet
facet normal 0.0 0.0 1.0
outer loop
vertex 100.0 100.0 100.0
vertex 100.0 200.0 100.0
vertex 200.0 200.0 100.0
endloop
endfacet
facet normal 0.0 0.0 -1.0
outer loop
vertex 100.0 100.0 200.0
vertex 200.0 100.0 200.0
vertex 200.0 200.0 200.0
endloop
endfacet
facet normal 0.0 0.0 -1.0
outer loop
vertex 100.0 100.0 200.0
vertex 100.0 200.0 200.0
vertex 200.0 200.0 200.0
endloop
endfacet
facet normal 0.0 1.0 0.0
outer loop
vertex 100.0 100.0 100.0
vertex 200.0 100.0 100.0
vertex 200.0 100.0 200.0
endloop
endfacet
facet normal 0.0 1.0 0.0
outer loop
vertex 100.0 100.0 100.0
vertex 100.0 100.0 200.0
vertex 200.0 100.0 200.0
endloop
endfacet
facet normal 0.0 -1.0 0.0
outer loop
vertex 100.0 200.0 100.0
vertex 200.0 200.0 100.0
vertex 200.0 200.0 200.0
endloop
endfacet
facet normal 0.0 -1.0 0.0
outer loop
vertex 100.0 200.0 100.0
vertex 100.0 200.0 200.0
vertex 200.0 200.0 200.0
endloop
endfacet
facet normal -1.0 0.0 0.0
outer loop
vertex 200.0 100.0 100.0
vertex 200.0 200.0 100.0
vertex 200.0 200.0 200.0
endloop
endfacet
facet normal -1.0 0.0 0.0
outer loop
vertex 200.0 100.0 100.0
vertex 200.0 100.0 200.0
vertex 200.0 200.0 200.0
endloop
endfacet
endsolid

1
examples/mesh/race_track.stl Symbolic link
View File

@ -0,0 +1 @@
../PACKAGES/machdyn/funnel_flow/boundary.stl

3
src/MACHDYN/fix_smd_wall_surface.cpp
View File

@ -260,8 +260,7 @@ void FixSMDWallSurface::read_triangles(int pass) {
double r1 = (center - vert[0]).norm();
double r2 = (center - vert[1]).norm();
double r3 = (center - vert[2]).norm();
double r = MAX(r1, r2);
r = MAX(r, r3);
double r = MAX(MAX(r1, r2), r3);
/*
* if atom/molecule is in my subbox, create it

4
src/atom.cpp
View File

@ -1722,7 +1722,7 @@ void Atom::set_mass(const char *file, int line, const char *str, int type_offset
void Atom::set_mass(const char *file, int line, int itype, double value)
{
if (mass == nullptr) error->all(file,line,"Cannot set mass for atom style {}", atom_style);
if (mass == nullptr) error->all(file,line, "Cannot set mass for atom style {}", atom_style);
if (itype < 1 || itype > ntypes)
error->all(file,line,"Invalid type {} for atom mass {}", itype, value);
if (value <= 0.0) error->all(file,line,"Invalid atom mass value {}", value);
@ -1738,7 +1738,7 @@ void Atom::set_mass(const char *file, int line, int itype, double value)
void Atom::set_mass(const char *file, int line, int /*narg*/, char **arg)
{
if (mass == nullptr) error->all(file,line,"Cannot set atom mass for atom style", atom_style);
if (mass == nullptr) error->all(file,line, "Cannot set atom mass for atom style {}", atom_style);
int lo,hi;
utils::bounds(file,line,arg[0],1,ntypes,lo,hi,error);

440
src/create_atoms.cpp
View File

@ -35,23 +35,30 @@
#include "random_park.h"
#include "region.h"
#include "special.h"
#include "text_file_reader.h"
#include "variable.h"
#include <cstring>
using namespace LAMMPS_NS;
using namespace MathConst;
using MathConst::MY_2PI;
using MathConst::MY_PI;
using MathConst::THIRD;
#define BIG 1.0e30
#define EPSILON 1.0e-6
#define LB_FACTOR 1.1
#define DEFAULT_MAXTRY 1000
static constexpr double BIG = 1.0e30;
static constexpr double EPSILON = 1.0e-6;
static constexpr double LB_FACTOR = 1.1;
static constexpr double INV_P_CONST = 0.7548777;
static constexpr double INV_SQ_P_CONST = 0.5698403;
static constexpr int DEFAULT_MAXTRY = 1000;
enum { BOX, REGION, SINGLE, RANDOM };
enum { BOX, REGION, SINGLE, RANDOM, MESH };
enum { ATOM, MOLECULE };
enum { COUNT, INSERT, INSERT_SELECTED };
enum { NONE, RATIO, SUBSET };
enum { BISECTION, QUASIRANDOM };
static constexpr const char *mesh_name[] = {"recursive bisection", "quasi-random"};
/* ---------------------------------------------------------------------- */
CreateAtoms::CreateAtoms(LAMMPS *lmp) : Command(lmp), basistype(nullptr) {}
@ -63,9 +70,7 @@ void CreateAtoms::command(int narg, char **arg)
if (domain->box_exist == 0)
error->all(FLERR, "Create_atoms command before simulation box is defined");
if (modify->nfix_restart_peratom)
error->all(FLERR,
"Cannot create_atoms after "
"reading restart file with per-atom info");
error->all(FLERR, "Cannot create_atoms after reading restart file with per-atom info");
// check for compatible lattice
@ -81,9 +86,10 @@ void CreateAtoms::command(int narg, char **arg)
// parse arguments
if (narg < 2) error->all(FLERR, "Illegal create_atoms command");
if (narg < 2) utils::missing_cmd_args(FLERR, "create_atoms", error);
ntype = utils::inumeric(FLERR, arg[0], false, lmp);
const char *meshfile;
int iarg;
if (strcmp(arg[1], "box") == 0) {
style = BOX;
@ -91,7 +97,7 @@ void CreateAtoms::command(int narg, char **arg)
region = nullptr;
} else if (strcmp(arg[1], "region") == 0) {
style = REGION;
if (narg < 3) error->all(FLERR, "Illegal create_atoms command");
if (narg < 3) utils::missing_cmd_args(FLERR, "create_atoms region", error);
region = domain->get_region_by_id(arg[2]);
if (!region) error->all(FLERR, "Create_atoms region {} does not exist", arg[2]);
region->init();
@ -99,18 +105,18 @@ void CreateAtoms::command(int narg, char **arg)
iarg = 3;
} else if (strcmp(arg[1], "single") == 0) {
style = SINGLE;
if (narg < 5) error->all(FLERR, "Illegal create_atoms command");
if (narg < 5) utils::missing_cmd_args(FLERR, "create_atoms single", error);
xone[0] = utils::numeric(FLERR, arg[2], false, lmp);
xone[1] = utils::numeric(FLERR, arg[3], false, lmp);
xone[2] = utils::numeric(FLERR, arg[4], false, lmp);
iarg = 5;
} else if (strcmp(arg[1], "random") == 0) {
style = RANDOM;
if (narg < 5) error->all(FLERR, "Illegal create_atoms command");
if (narg < 5) utils::missing_cmd_args(FLERR, "create_atoms random", error);
nrandom = utils::inumeric(FLERR, arg[2], false, lmp);
if (nrandom < 0) error->all(FLERR, "Illegal create_atoms command");
if (nrandom < 0) error->all(FLERR, "Illegal create_atoms number of random atoms {}", nrandom);
seed = utils::inumeric(FLERR, arg[3], false, lmp);
if (seed <= 0) error->all(FLERR, "Illegal create_atoms command");
if (seed <= 0) error->all(FLERR, "Illegal create_atoms random seed {}", seed);
if (strcmp(arg[4], "NULL") == 0)
region = nullptr;
else {
@ -120,8 +126,13 @@ void CreateAtoms::command(int narg, char **arg)
region->prematch();
}
iarg = 5;
} else if (strcmp(arg[1], "mesh") == 0) {
style = MESH;
if (narg < 3) utils::missing_cmd_args(FLERR, "create_atoms mesh", error);
meshfile = arg[2];
iarg = 3;
} else
error->all(FLERR, "Illegal create_atoms command");
error->all(FLERR, "Unknown create_atoms command option {}", arg[1]);
// process optional keywords
@ -138,31 +149,33 @@ void CreateAtoms::command(int narg, char **arg)
int subsetseed;
overlapflag = 0;
maxtry = DEFAULT_MAXTRY;
radscale = 1.0;
radthresh = domain->lattice->xlattice;
mesh_style = BISECTION;
mesh_density = 1.0;
nbasis = domain->lattice->nbasis;
basistype = new int[nbasis];
for (int i = 0; i < nbasis; i++) basistype[i] = ntype;
while (iarg < narg) {
if (strcmp(arg[iarg], "basis") == 0) {
if (iarg + 3 > narg) error->all(FLERR, "Illegal create_atoms command");
if (iarg + 3 > narg) utils::missing_cmd_args(FLERR, "create_atoms basis", error);
int ibasis = utils::inumeric(FLERR, arg[iarg + 1], false, lmp);
int itype = utils::inumeric(FLERR, arg[iarg + 2], false, lmp);
if (ibasis <= 0 || ibasis > nbasis || itype <= 0 || itype > atom->ntypes)
error->all(FLERR, "Invalid basis setting in create_atoms command");
error->all(FLERR, "Out of range basis setting '{} {}' in create_atoms command", ibasis,
itype);
basistype[ibasis - 1] = itype;
iarg += 3;
} else if (strcmp(arg[iarg], "remap") == 0) {
if (iarg + 2 > narg) error->all(FLERR, "Illegal create_atoms command");
if (iarg + 2 > narg) utils::missing_cmd_args(FLERR, "create_atoms remap", error);
remapflag = utils::logical(FLERR, arg[iarg + 1], false, lmp);
iarg += 2;
} else if (strcmp(arg[iarg], "mol") == 0) {
if (iarg + 3 > narg) error->all(FLERR, "Illegal create_atoms command");
if (iarg + 3 > narg) utils::missing_cmd_args(FLERR, "create_atoms mol", error);
int imol = atom->find_molecule(arg[iarg + 1]);
if (imol == -1)
error->all(FLERR,
"Molecule template ID for "
"create_atoms does not exist");
error->all(FLERR, "Molecule template ID {} for create_atoms does not exist", arg[iarg + 1]);
if ((atom->molecules[imol]->nset > 1) && (comm->me == 0))
error->warning(FLERR, "Molecule template for create_atoms has multiple molecules");
mode = MOLECULE;
@ -170,24 +183,25 @@ void CreateAtoms::command(int narg, char **arg)
molseed = utils::inumeric(FLERR, arg[iarg + 2], false, lmp);
iarg += 3;
} else if (strcmp(arg[iarg], "units") == 0) {
if (iarg + 2 > narg) error->all(FLERR, "Illegal create_atoms command");
if (iarg + 2 > narg) utils::missing_cmd_args(FLERR, "create_atoms units", error);
if (strcmp(arg[iarg + 1], "box") == 0)
scaleflag = 0;
else if (strcmp(arg[iarg + 1], "lattice") == 0)
scaleflag = 1;
else
error->all(FLERR, "Illegal create_atoms command");
error->all(FLERR, "Unknown create_atoms units option {}", arg[iarg + 1]);
iarg += 2;
} else if (strcmp(arg[iarg], "var") == 0) {
if (style == SINGLE) error->all(FLERR, "Illegal create_atoms command: "
"can't combine 'var' keyword with 'single' style!");
if (iarg + 2 > narg) error->all(FLERR, "Illegal create_atoms command");
if (style == SINGLE)
error->all(FLERR,
"Illegal create_atoms command: cannot use 'var' keyword with 'single' style!");
if (iarg + 2 > narg) utils::missing_cmd_args(FLERR, "create_atoms var", error);
delete[] vstr;
vstr = utils::strdup(arg[iarg + 1]);
varflag = 1;
iarg += 2;
} else if (strcmp(arg[iarg], "set") == 0) {
if (iarg + 3 > narg) error->all(FLERR, "Illegal create_atoms command");
if (iarg + 3 > narg) utils::missing_cmd_args(FLERR, "create_atoms set", error);
if (strcmp(arg[iarg + 1], "x") == 0) {
delete[] xstr;
xstr = utils::strdup(arg[iarg + 2]);
@ -198,10 +212,10 @@ void CreateAtoms::command(int narg, char **arg)
delete[] zstr;
zstr = utils::strdup(arg[iarg + 2]);
} else
error->all(FLERR, "Illegal create_atoms command");
error->all(FLERR, "Unknown create_atoms set option {}", arg[iarg + 2]);
iarg += 3;
} else if (strcmp(arg[iarg], "rotate") == 0) {
if (iarg + 5 > narg) error->all(FLERR, "Illegal create_atoms command");
if (iarg + 5 > narg) utils::missing_cmd_args(FLERR, "create_atoms rotate", error);
quat_user = 1;
double thetaone;
double axisone[3];
@ -210,33 +224,34 @@ void CreateAtoms::command(int narg, char **arg)
axisone[1] = utils::numeric(FLERR, arg[iarg + 3], false, lmp);
axisone[2] = utils::numeric(FLERR, arg[iarg + 4], false, lmp);
if (axisone[0] == 0.0 && axisone[1] == 0.0 && axisone[2] == 0.0)
error->all(FLERR, "Illegal create_atoms command");
error->all(FLERR, "Illegal create_atoms rotate arguments");
if (domain->dimension == 2 && (axisone[0] != 0.0 || axisone[1] != 0.0))
error->all(FLERR, "Invalid create_atoms rotation vector for 2d model");
MathExtra::norm3(axisone);
MathExtra::axisangle_to_quat(axisone, thetaone, quatone);
iarg += 5;
} else if (strcmp(arg[iarg], "ratio") == 0) {
if (iarg + 3 > narg) error->all(FLERR, "Illegal create_atoms command");
if (iarg + 3 > narg) utils::missing_cmd_args(FLERR, "create_atoms ratio", error);
subsetflag = RATIO;
subsetfrac = utils::numeric(FLERR, arg[iarg + 1], false, lmp);
subsetseed = utils::inumeric(FLERR, arg[iarg + 2], false, lmp);
if (subsetfrac <= 0.0 || subsetfrac > 1.0 || subsetseed <= 0)
error->all(FLERR, "Illegal create_atoms command");
error->all(FLERR, "Illegal create_atoms ratio settings");
iarg += 3;
} else if (strcmp(arg[iarg], "subset") == 0) {
if (iarg + 3 > narg) error->all(FLERR, "Illegal create_atoms command");
if (iarg + 3 > narg) utils::missing_cmd_args(FLERR, "create_atoms subset", error);
subsetflag = SUBSET;
nsubset = utils::bnumeric(FLERR, arg[iarg + 1], false, lmp);
subsetseed = utils::inumeric(FLERR, arg[iarg + 2], false, lmp);
if (nsubset <= 0 || subsetseed <= 0) error->all(FLERR, "Illegal create_atoms command");
if ((nsubset <= 0) || (subsetseed <= 0))
error->all(FLERR, "Illegal create_atoms subset settings");
iarg += 3;
} else if (strcmp(arg[iarg], "overlap") == 0) {
if (style != RANDOM)
error->all(FLERR, "Create_atoms overlap can only be used with random style");
if (iarg + 2 > narg) error->all(FLERR, "Illegal create_atoms command");
overlap = utils::numeric(FLERR, arg[iarg + 1], false, lmp);
if (overlap <= 0) error->all(FLERR, "Illegal create_atoms command");
if (overlap <= 0) error->all(FLERR, "Illegal create_atoms overlap value: {}", overlap);
overlapflag = 1;
iarg += 2;
} else if (strcmp(arg[iarg], "maxtry") == 0) {
@ -244,22 +259,42 @@ void CreateAtoms::command(int narg, char **arg)
error->all(FLERR, "Create_atoms maxtry can only be used with random style");
if (iarg + 2 > narg) error->all(FLERR, "Illegal create_atoms command");
maxtry = utils::inumeric(FLERR, arg[iarg + 1], false, lmp);
if (maxtry <= 0) error->all(FLERR,"Illegal create_atoms command");
if (maxtry <= 0) error->all(FLERR, "Illegal create_atoms command");
iarg += 2;
} else if (strcmp(arg[iarg], "meshmode") == 0) {
if (style != MESH)
error->all(FLERR, "Create_atoms meshmode can only be used with mesh style");
if (iarg + 3 > narg) utils::missing_cmd_args(FLERR, "create_atoms meshmode", error);
if (strcmp(arg[iarg + 1], "bisect") == 0) {
mesh_style = BISECTION;
radthresh = utils::numeric(FLERR, arg[iarg + 2], false, lmp);
} else if (strcmp(arg[iarg + 1], "qrand") == 0) {
mesh_style = QUASIRANDOM;
mesh_density = utils::numeric(FLERR, arg[iarg + 2], false, lmp);
} else
error->all(FLERR, "Unknown create_atoms meshmode {}", arg[iarg + 2]);
iarg += 3;
} else if (strcmp(arg[iarg], "radscale") == 0) {
if (iarg + 2 > narg) utils::missing_cmd_args(FLERR, "create_atoms radscale", error);
if (style != MESH)
error->all(FLERR, "Create_atoms radscale can only be used with mesh style");
if (!atom->radius_flag)
error->all(FLERR, "Must have atom attribute radius to set radscale factor");
radscale = utils::numeric(FLERR, arg[iarg + 1], false, lmp);
if (radscale <= 0.0) error->all(FLERR, "Illegal create_atoms radscale value: {}", radscale);
iarg += 2;
} else
error->all(FLERR, "Illegal create_atoms command");
error->all(FLERR, "Illegal create_atoms command option {}", arg[iarg]);
}
// error checks and further setup for mode = MOLECULE
if (mode == ATOM) {
if (ntype <= 0 || ntype > atom->ntypes)
if ((ntype <= 0) || (ntype > atom->ntypes))
error->all(FLERR, "Invalid atom type in create_atoms command");
} else if (mode == MOLECULE) {
if (onemol->xflag == 0)
error->all(FLERR, "Create_atoms molecule must have coordinates");
if (onemol->typeflag == 0)
error->all(FLERR, "Create_atoms molecule must have atom types");
if (onemol->xflag == 0) error->all(FLERR, "Create_atoms molecule must have coordinates");
if (onemol->typeflag == 0) error->all(FLERR, "Create_atoms molecule must have atom types");
if (ntype + onemol->ntypes <= 0 || ntype + onemol->ntypes > atom->ntypes)
error->all(FLERR, "Invalid atom type in create_atoms mol command");
if (onemol->tag_require && !atom->tag_enable)
@ -274,6 +309,12 @@ void CreateAtoms::command(int narg, char **arg)
ranmol = new RanMars(lmp, molseed + comm->me);
}
if (style == MESH) {
if (mode == MOLECULE)
error->all(FLERR, "Create_atoms mesh is not compatible with the 'mol' option");
if (scaleflag) error->all(FLERR, "Create_atoms mesh must use 'units box' option");
}
ranlatt = nullptr;
if (subsetflag != NONE) ranlatt = new RanMars(lmp, subsetseed + comm->me);
@ -316,7 +357,7 @@ void CreateAtoms::command(int narg, char **arg)
// lattice to box, but not consistent with other uses of units=lattice
// triclinic remapping occurs in add_single()
if (style == BOX || style == REGION) {
if ((style == BOX) || (style == REGION) || (style == MESH)) {
if (nbasis == 0) error->all(FLERR, "Cannot create atoms with undefined lattice");
} else if (scaleflag == 1) {
xone[0] *= domain->lattice->xlattice;
@ -415,6 +456,8 @@ void CreateAtoms::command(int narg, char **arg)
add_single();
else if (style == RANDOM)
add_random();
else if (style == MESH)
add_mesh(meshfile);
else
add_lattice();
@ -670,7 +713,7 @@ void CreateAtoms::add_random()
if (overlapflag) {
double odist = overlap;
if (mode == MOLECULE) odist += onemol->molradius;
odistsq = odist*odist;
odistsq = odist * odist;
}
// random number generator, same for all procs
@ -718,7 +761,7 @@ void CreateAtoms::add_random()
// insert Nrandom new atom/molecule into simulation box
int ntry,success;
int ntry, success;
int ninsert = 0;
for (int i = 0; i < nrandom; i++) {
@ -766,7 +809,7 @@ void CreateAtoms::add_random()
dely = xone[1] - x[i][1];
delz = xone[2] - x[i][2];
domain->minimum_image(delx, dely, delz);
distsq = delx*delx + dely*dely + delz*delz;
distsq = delx * delx + dely * dely + delz * delz;
if (distsq < odistsq) {
reject = 1;
break;
@ -813,6 +856,302 @@ void CreateAtoms::add_random()
delete random;
}
/* ----------------------------------------------------------------------
add atom at center of triangle
or split into two halves along the longest side and recurse
------------------------------------------------------------------------- */
int CreateAtoms::add_bisection(const double vert[3][3], tagint molid)
{
double center[3], temp[3];
MathExtra::add3(vert[0], vert[1], center);
MathExtra::add3(center, vert[2], temp);
MathExtra::scale3(THIRD, temp, center);
MathExtra::sub3(center, vert[0], temp);
double ravg = MathExtra::len3(temp);
MathExtra::sub3(center, vert[1], temp);
ravg += MathExtra::len3(temp);
MathExtra::sub3(center, vert[2], temp);
ravg += MathExtra::len3(temp);
ravg *= THIRD;
// if the average distance of the vertices from the center is larger than the
// lattice parameter, the triangle is split it in half along its longest side
int ilocal = 0;
if (ravg > radthresh) {
double vert1[3][3], vert2[3][3], side[3][3];
// determine side vectors
MathExtra::sub3(vert[0], vert[1], side[0]);
MathExtra::sub3(vert[1], vert[2], side[1]);
MathExtra::sub3(vert[2], vert[0], side[2]);
// determine longest side
const double l1 = MathExtra::len3(side[0]);
const double l2 = MathExtra::len3(side[1]);
const double l3 = MathExtra::len3(side[2]);
int isplit = 0;
if (l1 < l2) {
isplit = 1;
if (l2 < l3) isplit = 2;
} else {
if (l1 < l3) isplit = 2;
}
MathExtra::scaleadd3(-0.5, side[isplit], vert[isplit], temp);
int inext = (isplit + 1) % 3;
// create two new triangles
for (int i = 0; i < 3; ++i) {
for (int j = 0; j < 3; ++j) {
vert1[i][j] = vert2[i][j] = vert[i][j];
vert1[isplit][j] = temp[j];
vert2[inext][j] = temp[j];
}
}
ilocal = add_bisection(vert1, molid);
ilocal += add_bisection(vert2, molid);
} else {
/*
* if atom/molecule is in my subbox, create it
* ... use x to hold triangle center
* ... radius is the mmaximal distance from triangle center to all vertices
*/
if ((center[0] >= sublo[0]) && (center[0] < subhi[0]) && (center[1] >= sublo[1]) &&
(center[1] < subhi[1]) && (center[2] >= sublo[2]) && (center[2] < subhi[2])) {
atom->avec->create_atom(ntype, center);
int idx = atom->nlocal - 1;
if (atom->radius_flag) atom->radius[idx] = ravg * radscale;
if (atom->molecule_flag) atom->molecule[idx] = molid;
++ilocal;
}
}
return ilocal;
}
/* ----------------------------------------------------------------------
add monodisperse atoms by mapping quasirandom sequence onto a triangle
method by Martin Roberts 2019 http://extremelearning.com.au/evenly-distributing-points-in-a-triangle/
------------------------------------------------------------------------- */
int CreateAtoms::add_quasirandom(const double vert[3][3], tagint molid)
{
double ab[3], ac[3], bc[3], temp[3], point[3], ref[3];
double lab, lac, lbc, area, xi, yi;
double seed = 0.5;
// Order vertices such that a has the largest angle
MathExtra::sub3(vert[1], vert[0], ab);
MathExtra::sub3(vert[2], vert[0], ac);
MathExtra::sub3(vert[2], vert[1], bc);
lab = MathExtra::len3(ab);
lac = MathExtra::len3(ac);
lbc = MathExtra::len3(bc);
if (lac > lab && lac > lbc) {
// B has the largest angle, relabel as A
MathExtra::scale3(-1.0, ab);
MathExtra::copy3(bc, ac);
MathExtra::copy3(vert[1], ref);
} else if (lab > lac && lab > lbc) {
// C has the largest angle, relabel as A
MathExtra::scale3(-1.0, ac);
MathExtra::copy3(bc, ab);
MathExtra::scale3(-1.0, ab);
MathExtra::copy3(vert[2], ref);
} else {
MathExtra::copy3(vert[0], ref);
}
// Estimate number of particles from area
MathExtra::cross3(ab, ac, temp);
area = 0.5 * MathExtra::len3(temp);
int nparticles = ceil(mesh_density * area);
// estimate radius from number of particles and area
double rad = sqrt(area/MY_PI/nparticles);
for (int i = 0; i < nparticles; i++) {
// Define point in unit square
xi = (i + 1) * INV_P_CONST;
yi = (i + 1) * INV_SQ_P_CONST;
xi += seed;
yi += seed;
xi = std::fmod(xi, 1.0);
yi = std::fmod(yi, 1.0);
// Map to triangle using parallelogram method
if ((xi + yi) < 1) {
MathExtra::scale3(xi, ac, point);
MathExtra::scale3(yi, ab, temp);
MathExtra::add3(point, temp, point);
} else {
xi = 1.0 - xi;
yi = 1.0 - yi;
MathExtra::scale3(xi, ac, point);
MathExtra::scale3(yi, ab, temp);
MathExtra::add3(point, temp, point);
}
MathExtra::add3(point, ref, point);
// if atom/molecule is in my subbox, create it
if ((point[0] >= sublo[0]) && (point[0] < subhi[0]) && (point[1] >= sublo[1]) &&
(point[1] < subhi[1]) && (point[2] >= sublo[2]) && (point[2] < subhi[2])) {
atom->avec->create_atom(ntype, point);
int idx = atom->nlocal - 1;
if (atom->molecule_flag) atom->molecule[idx] = molid;
if (atom->radius_flag) atom->radius[idx] = rad * radscale;
}
}
return nparticles;
}
/* ----------------------------------------------------------------------
add atoms at center of triangulated mesh
------------------------------------------------------------------------- */
void CreateAtoms::add_mesh(const char *filename)
{
double vert[3][3];
tagint *mol = atom->molecule;
int nlocal_previous = atom->nlocal;
bigint atomlocal = 0, atomall = 0, ntriangle = 0;
// find next molecule ID, if available
tagint molid = 0;
if (atom->molecule_flag) {
tagint max = 0;
for (int i = 0; i < nlocal_previous; i++) max = MAX(max, mol[i]);
tagint maxmol;
MPI_Allreduce(&max, &maxmol, 1, MPI_LMP_TAGINT, MPI_MAX, world);
molid = maxmol + 1;
}
FILE *fp = fopen(filename, "rb");
if (fp == nullptr) error->one(FLERR, "Cannot open file {}: {}", filename, utils::getsyserror());
// first try reading the file in ASCII format
TextFileReader reader(fp, "STL mesh");
try {
char *line = reader.next_line();
if (!line || !utils::strmatch(line, "^solid"))
throw TokenizerException("Invalid STL mesh file format", "");
line += 6;
if (comm->me == 0)
utils::logmesg(lmp, "Reading STL object {} from text file {}\n", utils::trim(line), filename);
while ((line = reader.next_line())) {
// next line is facet line with 5 words
auto values = utils::split_words(line);
// otherwise stop reading
if ((values.size() != 5) || !utils::strmatch(values[0], "^facet")) break;
// ignore normal
line = reader.next_line(2);
if (!line || !utils::strmatch(line, "^ *outer *loop"))
throw TokenizerException("Error reading outer loop", "");
for (int k = 0; k < 3; ++k) {
line = reader.next_line(4);
values = utils::split_words(line);
if ((values.size() != 4) || !utils::strmatch(values[0], "^vertex"))
throw TokenizerException("Error reading vertex", "");
vert[k][0] = utils::numeric(FLERR, values[1], false, lmp);
vert[k][1] = utils::numeric(FLERR, values[2], false, lmp);
vert[k][2] = utils::numeric(FLERR, values[3], false, lmp);
}
line = reader.next_line(1);
if (!line || !utils::strmatch(line, "^ *endloop"))
throw TokenizerException("Error reading endloop", "");
line = reader.next_line(1);
if (!line || !utils::strmatch(line, "^ *endfacet"))
throw TokenizerException("Error reading endfacet", "");
// now we have the three vertices ... proceed with adding atoms
++ntriangle;
if (mesh_style == BISECTION) {
// add in center of triangle or bisecting recursively along longest edge
// as needed to get the desired atom density/radii
atomlocal += add_bisection(vert, molid);
} else if (mesh_style == QUASIRANDOM) {
// add quasi-random distribution of atoms
atomlocal += add_quasirandom(vert, molid);
}
}
} catch (std::exception &e) {
// if ASCII failed for the first line, try reading as binary
if (utils::strmatch(e.what(), "^Invalid STL mesh file format")) {
char title[80];
float triangle[12];
uint32_t ntri;
uint16_t attr;
size_t count;
rewind(fp);
count = fread(title, sizeof(char), 80, fp);
title[79] = '\0';
count = fread(&ntri, sizeof(ntri), 1, fp);
if (count <= 0) {
error->all(FLERR, "Error reading STL file {}: {}", filename, utils::getsyserror());
} else {
if (comm->me == 0)
utils::logmesg(lmp, "Reading STL object {} from binary file {}\n", utils::trim(title),
filename);
}
for (uint32_t i = 0U; i < ntri; ++i) {
count = fread(triangle, sizeof(float), 12, fp);
if (count != 12)
error->all(FLERR, "Error reading STL file {}: {}", filename, utils::getsyserror());
count = fread(&attr, sizeof(attr), 1, fp);
for (int j = 0; j < 3; ++j)
for (int k = 0; k < 3; ++k) vert[j][k] = triangle[3 * j + 3 + k];
++ntriangle;
if (mesh_style == BISECTION) {
// add in center of triangle or bisecting recursively along longest edge
// as needed to get the desired atom density/radii
atomlocal += add_bisection(vert, molid);
} else if (mesh_style == QUASIRANDOM) {
// add quasi-random distribution of atoms
atomlocal += add_quasirandom(vert, molid);
}
}
} else {
error->all(FLERR, "Error reading triangles from file {}: {}", filename, e.what());
}
}
if (ntriangle > 0) {
MPI_Allreduce(&atomlocal, &atomall, 1, MPI_LMP_BIGINT, MPI_SUM, world);
double ratio = (double) atomall / (double) ntriangle;
if (comm->me == 0)
utils::logmesg(lmp, " read {} triangles with {:.2f} atoms per triangle added in {} mode\n",
ntriangle, ratio, mesh_name[mesh_style]);
}
}
/* ----------------------------------------------------------------------
add many atoms by looping over lattice
------------------------------------------------------------------------- */
@ -982,8 +1321,7 @@ void CreateAtoms::loop_lattice(int action)
// if variable test specified, eval variable
if (varflag && vartest(x) == 0)
continue;
if (varflag && vartest(x) == 0) continue;
// test if atom/molecule position is in my subbox

5
src/create_atoms.h
View File

@ -41,6 +41,7 @@ class CreateAtoms : public Command {
double subsetfrac;
int *basistype;
double xone[3], quatone[4];
double radthresh, radscale, mesh_density;
int varflag, vvar, xvar, yvar, zvar;
char *vstr, *xstr, *ystr, *zstr;
@ -53,6 +54,7 @@ class CreateAtoms : public Command {
int *flag; // flag subset of particles to insert on lattice
int *next;
int mesh_style;
class Region *region;
class Molecule *onemol;
@ -64,6 +66,9 @@ class CreateAtoms : public Command {
void add_single();
void add_random();
void add_mesh(const char *);
int add_bisection(const double [3][3], tagint);
int add_quasirandom(const double [3][3], tagint);
void add_lattice();
void loop_lattice(int);
void add_molecule(double *);

5
tools/Makefile
View File

@ -5,7 +5,7 @@
#
all:
$(MAKE) binary2txt chain micelle2d
$(MAKE) binary2txt chain micelle2d stl_bin2txt
binary2txt: binary2txt.o
g++ -g binary2txt.o -o binary2txt
@ -19,6 +19,9 @@ micelle2d: micelle2d.o
thermo_extract: thermo_extract.o
gcc -g thermo_extract.o -o thermo_extract
stl_bin2txt: stl_bin2txt.o
g++ -g stl_bin2txt.o -o stl_bin2txt
clean:
rm binary2txt chain micelle2d
rm thermo_extract

1
tools/README
View File

@ -46,6 +46,7 @@ replica tool to reorder LAMMPS replica trajectories according to
singularity Singularity container descriptions suitable for LAMMPS development
smd convert Smooth Mach Dynamics triangles to VTK
spin perform a cubic polynomial interpolation of a GNEB MEP
stl_bin2txt convert binary STL files to ASCII
swig Interface file and demo scripts for SWIG wrappers for the LAMMPS C library interface
valgrind suppression files for use with valgrind's memcheck tool
vim add-ons to VIM editor for editing LAMMPS input scripts

99
tools/stl_bin2txt.cpp Normal file
View File

@ -0,0 +1,99 @@
/* -----------------------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
https://www.lammps.org/
Steve Plimpton, sjplimp@sandia.gov, Sandia National Laboratories
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.
------------------------------------------------------------------------ */
// Convert a binary STL file to ASCII format
// Contributing author: Axel Kohlmeyer, Temple U, akohlmey at gmail.com
//
// Specs of the format taken from: https://en.wikipedia.org/wiki/STL_(file_format)
#include <cerrno>
#include <cstdint>
#include <cstdio>
#include <cstring>
int main(int argc, char **argv)
{
FILE *in, *out;
char title[80];
float normal[3], vert1[3], vert2[3], vert3[3];
uint32_t ntriangles;
size_t count;
uint16_t attributes;
if (argc != 3) {
printf("Usage: %s <input file> <output file>\n", argv[0]);
return 1;
}
in = fopen(argv[1], "rb");
if (!in) {
printf("Error opening input file %s: %s\n", argv[1], strerror(errno));
return 2;
}
out = fopen(argv[2], "w");
if (!out) {
printf("Error opening output file %s: %s\n", argv[1], strerror(errno));
return 3;
}
/* read header */
count = fread(title, sizeof(char), sizeof(title), in);
if (count != sizeof(title)) {
printf("Error reading binary STL header: %s\n", strerror(errno));
return 4;
}
count = strlen(title);
if (count == 0) snprintf(title, 80, "STL object from file %s", argv[1]);
/* read triangle count */
count = fread(&ntriangles, sizeof(uint32_t), 1, in);
if (count != 1) {
printf("Error reading binary STL triangle count: %s\n", strerror(errno));
return 5;
}
/* write header */
printf("Converting: %s with %u triangles\n", title, ntriangles);
fprintf(out, "solid %s\n", title);
/* loop over triangles */
for (uint32_t i = 0; i < ntriangles; ++i) {
count = fread(normal, sizeof(float), 3, in);
count += fread(vert1, sizeof(float), 3, in);
count += fread(vert2, sizeof(float), 3, in);
count += fread(vert3, sizeof(float), 3, in);
if (count != 12) {
printf("Error reading binary STL vertices: %s\n", strerror(errno));
return 6;
}
count = fread(&attributes, sizeof(uint16_t), 1, in);
if (count != 1) {
printf("Error reading binary STL facet attributes: %s\n", strerror(errno));
return 7;
}
fprintf(out, " facet normal %e %e %e\n", normal[0], normal[1], normal[2]);
fputs(" outer loop\n", out);
fprintf(out, " vertex %e %e %e\n", vert1[0], vert1[1], vert1[2]);
fprintf(out, " vertex %e %e %e\n", vert2[0], vert2[1], vert2[2]);
fprintf(out, " vertex %e %e %e\n", vert3[0], vert3[1], vert3[2]);
fputs(" endloop\n endfacet\n", out);
if (ferror(out)) {
printf("Error writing text STL facet: %s\n", strerror(errno));
return 7;
}
}
fprintf(out, "endsolid %s\n", title);
fclose(in);
fclose(out);
return 0;
}