Modeling Cubes with body/rounded/polyhedron produces Segmentation fault

LAMMPS LAMMPS Beginners
1

Hello, I am new to LAMMPS and as a starting point I have adapted one of the installed examples. I started with the example: mylammps\examples\body\in.cubes

My goal is to create a stack of cubes within a box, and then release one side of the box. I have been able to create the stack of cubes and save a restart file. My second input file runs for a few 2,000 time steps and then fails with the following error:

*** stack smashing detected ***: terminated
*** Process received signal ***
Signal: Aborted (6)
Signal code: (-6)
*** Process received signal ***
Signal: Segmentation fault (11)
Signal code: (128)
Failing at address: (nil)
Segmentation fault

Any help or advice would be greatly appreciated.

To create the stack of cubes, I have used this input file:

# pouring 3d rounded polyhedron bodies

variable    steps index 10000
variable    steps_1 index 50000

units       lj
boundary    fm fm fm
comm_modify vel yes

atom_style  body rounded/polyhedron 1 8
atom_modify map array

region		reg block 0 25 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 0
variable c_t        equal 0
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     object molecule.cube molecule.tetra toff 1 &
#             molecule.rod3d toff 2 molecule.point3d toff 3

region       slab block 5 20 30 45 25 45 units box
fix          ins all pour 20 0 4767548 vol 0.4 10 region slab mol object
#fix          ins all pour 600 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

# Boundary
fix walls all wall/reflect xlo 0 xhi 25 ylo 25 yhi 50 zlo 0 zhi 50

#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

dump	     2 all image 1000 image.*.jpg type type &
	     zoom 1.5 adiam 1.5 body type 0 0 view 75 15
dump_modify  2 pad 6

thermo_style custom step atoms ke pe etotal press

thermo       1000

# Restart
restart 10000 poly.restart

run	     ${steps}

fix          ins all pour 20 0 4767548 vol 0.4 10 region slab mol object

run	     ${steps}

fix          ins all pour 20 0 4767548 vol 0.4 10 region slab mol object

run	     ${steps}

fix          ins all pour 20 0 4767548 vol 0.4 10 region slab mol object

run	     ${steps}

fix          ins all pour 20 0 4767548 vol 0.4 10 region slab mol object

run	     ${steps_1}

fix          ins all pour 20 0 4767548 vol 0.4 10 region slab mol object

run	     ${steps_1}

fix          ins all pour 20 0 4767548 vol 0.4 10 region slab mol object

run	     ${steps_1}

fix          ins all pour 20 0 4767548 vol 0.4 10 region slab mol object

run	     ${steps_1}

fix          ins all pour 20 0 4767548 vol 0.4 10 region slab mol object

run	     ${steps_1}

fix          ins all pour 20 0 4767548 vol 0.4 10 region slab mol object

run	     ${steps_1}

fix          ins all pour 20 0 4767548 vol 0.4 10 region slab mol object

run	     ${steps_1}

To remove one boundary and let the particles flow out, I have used:

# Reading restart file

# Read restart file
read_restart poly.restart.390000


variable    steps index 200000

#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 100 0 100 0 150 units box
#create_box	4 reg
#change_box all x final 0 200 y final 0 200 z final 0 150 boundary p p f remap units box
#change_box all x final 0 100

## Variables
variable cut_inner  equal 0.5
variable k_n        equal 100

variable k_na       equal 5
variable c_n        equal 0
#20
variable c_t        equal 0
#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 10.0 spherical 0.0 -180.0

# Boundary
fix walls all wall/reflect xlo 0 xhi 25 ylo 0 yhi 50 zlo 0 zhi 50

# Bottom Boundary
fix          4 all wall/body/polyhedron 2000 50 50 zplane 0.0 NULL

thermo_style custom step atoms ke pe etotal press

thermo       1000

# Images
dump	     2 all image 1000 image.*.jpg type type &
	     zoom 1.5 adiam 1.5 body type 0 0 view 75 15
dump_modify  2 pad 6

# Restart
restart 10000 poly.restart

run	     ${steps}
2

What version of LAMMPS are you using? and on what platform?

3

Thank you for the quick response. The version is: LAMMPS (29 Sep 2021 - Update 2). I am using Windows Server 2022, but I am running LAMMPS through WSL2 with Ubuntu.

4

Can you please try with the latest feature release version, 28 March 2023?

I do not see any issues with the current development tree using an abbreviated input.
Also with a more recent version, you don’t need to use fix pour multiple times, but you can create the objects in one go with something like this:

create_atoms 0 random 200 56474 slab overlap 0.5  mol object 32425
minimize 0.0 0.0 1000 10000
reset_timestep 0

If you still see issues, I would suspect that fix wall/reflect may be causing problems. I would suggest to replace it with fix wall/body/polyhedron

5

Thank you for your suggestion. I downloaded the 28 March 2023 version. I also replaced the wall/reflect with wall/body/polygon. LAMMPS is now through a segmentation fault first step: particle stacking.

Any advice would be greatly appreciated.

# pouring 3d rounded polyhedron bodies

variable    steps index 50000

units       lj
boundary    fm fm fm
comm_modify vel yes

atom_style  body rounded/polyhedron 1 8
atom_modify map array

region		reg block 0 25 0 200 0 200 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 0
variable c_t        equal 0
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

region       slab block 5 25 155 195 5 195 units box

# Create boxes
create_atoms 0 random 800 56474 slab overlap 0.5  mol object 32425

# Minimze the box energy
minimize 0.0 0.0 1000 50000
reset_timestep 0

# Boundary
fix          4 all wall/body/polyhedron 2000 50 50 xplane 0 25 
fix          5 all wall/body/polyhedron 2000 50 50 yplane 150 200 
fix          6 all wall/body/polyhedron 2000 50 50 zplane 0 200 

dump	     2 all image 10000 image.*.jpg type type &
	     zoom 1.5 adiam 1.5 body type 0 0 view 75 15 size 640 480
dump_modify  2 pad 6

thermo_style custom step atoms ke pe etotal press

thermo       1000

# Restart
restart 10000 poly.restart

run	     ${steps}
6

That appears to be a genuine bug in the pair style. Perhaps @ndtrung can have a look.

For better performance, you may want to add the line:

fix          8 all balance 500 0.0 shift xyz 20 1.0

to your input to get to the segfault faster… :wink:

1 Like
7

@tshoemaker For debugging purposes, could you try increase the overlap value in the create_atoms command to make sure that there is no overlap between the cubes in the initial configuration? (Note that in molecule.cube the cube length is 2.0.)

8

FWIW, the segfault is here (also with overlap set to 2.0):

Program terminated with signal SIGSEGV, Segmentation fault.
#0  0x00007fcffd6496e5 in LAMMPS_NS::PairBodyRoundedPolyhedron::rescale_cohesive_forces (this=this@entry=0xb6dfb0, x=x@entry=0x7fcfb3c51010, f=f@entry=0x7fcfb3b4d010, 
    torque=torque@entry=0x7fcfb3a07010, contact_list=contact_list@entry=0x7ffd82440430, num_contacts=@0x7ffd8244042c: 33, itype=1, jtype=1, facc=0x7ffd82441130)
    at /home/akohlmey/compile/lammps/src/BODY/pair_body_rounded_polyhedron.cpp:1734
1734	    f[ibody][0] += fx;

and “ibody” has an out of range value like -591321008 or 1818066054

9

Thanks @akohlmey I will take a closer look. I suppose the segfault persist with larger overlap (say 2\sqrt(3)+rounded_diameter ~ 3.5 + 1.0 = 4.5). Apparently, there’s some issue that needs to be addressed with the pair style. The out of range value of ibody could be a consequence of some missing contacts that happened earlier.

10

The simulations appear to be working with an overlap of 5.0. Thank you for your guidance.

On a related but separate note, what is the best way to export these simulations to VTK? I tried adding the following command:

dump dmpvtk all vtk 100 dump*.myforce.vtk

I received the following error:

ERROR: Unrecognized dump style ‘vtk’ is part of the VTK package which is not enabled in this LAMMPS binary. (src/output.cpp:769)

Must I compile LAMMPS from source?

11

Yes.

12

I tried up to 3.0 but not larger.

Perhaps there could be some test, e.g. in init_style() whether bodies overlap?

13

I was able to get the first part of the simulation running by using a larger spacing between particles. I also switched to using a lattice. The first input file is:

pouring 3d rounded polyhedron bodies

variable steps index 100000

units lj
boundary fm fm fm
comm_modify vel yes

atom_style body rounded/polyhedron 1 8
atom_modify map array

region reg block 0 25 0 200 0 200 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 0.1
variable c_t equal 0.1
variable mu equal 0.5
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

region slab block 2 22 2 47 2 198 units box

Create boxes

lattice sc 0.01
create_atoms 0 region slab mol object 32425

Boundary

fix 4 all wall/body/polyhedron 2000 50 50 xplane 0 25
fix 5 all wall/body/polyhedron 2000 50 50 yplane 0 50
fix 6 all wall/body/polyhedron 2000 50 50 zplane 0 200

Load balancing accross processors

fix 8 all balance 500 0.0 shift xyz 20 1.0

dump 2 all image 1000 image.*.jpg type type &
zoom 1.5 adiam 1.5 body type 0 0 view 75 15 size 5000 5000
dump_modify 2 pad 6

thermo_style custom step atoms ke pe etotal press

thermo 1000

Restart

restart 10000 poly.restart

run ${steps}

The second part of the simulation has an input file:

Reading restart file

Read restart file

read_restart poly.restart.100000

variable steps index 20000

Variables

variable cut_inner equal 0.5
variable k_n equal 100

variable k_na equal 5
variable c_n equal 0
#20
variable c_t equal 0
#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

Boundary

fix 4 all wall/body/polyhedron 2000 50 50 xplane 0 25
fix 5 all wall/body/polyhedron 2000 50 50 yplane 0 200
fix 6 all wall/body/polyhedron 2000 50 50 zplane 0 200

thermo_style custom step atoms ke pe etotal press

thermo 1000

Images

dump 2 all image 1000 image.*.jpg type type &
zoom 1.5 adiam 1.5 body type 0 0 view 75 15 size 5000 5000
dump_modify 2 pad 6

Restart

restart 1000 poly.restart

run ${steps}

The second part of the simulation is hitting a segmentation fault at time step 105,000. Do you have any recommendations? Is this the same issue?

14

@akohlmey PairBodyRoundedPolyhedron::init_style() computes the max enclosing radius of each atom type (stored in the merad array). Also, the BodyRoundedPolyhedron::data_body() function that computes the enclosing radius of the newly created atom and stores in atom->radius, which is then used by FixPour to test for atom overlaps.

15

@tshoemaker Thanks for the update. Could you try with a larger value of k_n, such as 200 and 500? When the particles are closely packed under an external field (gravity here), it could be that the softness (or stiffness) of the repulsive forces for a given time step (0.001) leads to numerical instability. Another way is to reduce the time step to 0.0005.

I am not sure if it is a bug in the code for overlap detection, or if it is the combination of the model parameters and time integration that led to the segfault, or both. Please keep me posted with your tests.

16

Thanks for the advice. I have been able to get past the original segmentation fault issues. I will post my running input files below.

However, I tried to transition to using SI units instead of the LJ units and I have been unsuccessful. Using LJ produces reasonable behavior with the LAMMPS tutorial settings. See image dump examples below:
Time Step = 0
image
Time Step = 500,000
image

When I use spring constants and time steps from published DEM sources, the particles overlap and pass through each other. See image dump examples below:
Time Step = 0
image
Time Step = 315,000 (note overlap)
image
Time Step = 500,000
image

Any advice would be greatly appreciated.

Working Input File with LJ units

# pouring 3d rounded polyhedron bodies

variable    steps index 500000

units       lj
boundary    fm fm fm
comm_modify vel yes

atom_style  body rounded/polyhedron 1 8
atom_modify map array

region		reg block 0 10 0 10 0 20 units box
create_box	4 reg

variable cut_inner  equal 0.0
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.1 bin
neigh_modify every 1 delay 0 check yes

# Set integration and gravity
timestep     0.001
fix          1 all nve/body
fix          2 all gravity 1.0 spherical 0.0 -180.0

# Create Particles
molecule     object molecule.cube.1 
region       slab block 1 9 1 9 1 6 units box
lattice 		sc 0.05
create_atoms 0  region slab  mol object 32425

# Boundaries
fix          4 all wall/body/polyhedron ${k_n} ${c_n} ${c_t} xplane 0 50
fix          5 all wall/body/polyhedron ${k_n} ${c_n} ${c_t} yplane 0 50
fix          6 all wall/body/polyhedron ${k_n} ${c_n} ${c_t} zplane 0 50

# Load balancing accross processors
fix          8 all balance 10000 0.0 shift xyz 20 1.0

#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       ${steps}

dump	     2 all image 10000 image.*.jpg type type &
	     zoom 1.5 adiam 1.5 body type 0 0 view 75 15
#dump_modify  2 pad 6
#dump 3 all custom 1000 dump id x y z fx fy fz

run	     ${steps}

Working cube file for LJ units

# 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

Input file with SI units that is not properly working

# dropping 3d rounded polyhedron bodies

variable    steps index 500000

units       si
boundary    fm fm fm
comm_modify vel yes

atom_style  body rounded/polyhedron 1 50
atom_modify map array

region		reg block 0.0 0.01 0.0 0.01 0.0 0.02 units box
create_box	4 reg

variable cut_inner  equal 0.0
variable k_n        equal 160000
variable k_na       equal 130000
variable c_n        equal 16000
variable c_t        equal 13000
variable mu         equal 0.60
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.00025 bin
neigh_modify every 1 delay 0 check yes

timestep     0.000000403113

fix          1 all nve/body
fix          2 all gravity 9.81 spherical 0.0 -180.0

# Load molecule shape and define region
molecule     object1 molecule.cube.1mm
region       slab1 block 0.001 0.002 0.001 0.007 0.001 0.005 units box

# Create boxes
lattice sc 0.002
create_atoms 0  region slab1  mol object1 32425

# Boundary
fix          4 all wall/body/polyhedron ${k_n} 0 0 xplane 0 0.01 
fix          5 all wall/body/polyhedron ${k_n} 0 0 yplane 0 0.01 
fix          6 all wall/body/polyhedron ${k_n} 0 0 zplane 0 0.02

# Load balancing accross processors
fix          8 all balance 10000 0.0 shift xyz 20 1.0

# Ouput
dump	     2 all image 1000 image.*.jpg type type center s 0.5 0.5 0.2 &
	     zoom 1.0 adiam 1.5 body type 0 0 view 75 15 size 2000 1000
dump_modify  2 pad 6

thermo_style custom step atoms ke pe etotal press
thermo       1000

# Run
run	     ${steps}

Input file with SI units that is not properly working

# 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 0.0000026

Body Integers

8 12 6

Body Doubles

0.000000166667 0.000000166667 0.000000166667 0 0 0
0.0005 0.0005 0.0005
0.0005 -0.0005 0.0005
-0.0005 -0.0005 0.0005
-0.0005 0.0005 0.0005
0.0005 0.0005 -0.0005
0.0005 -0.0005 -0.0005
-0.0005 -0.0005 -0.0005
-0.0005 0.0005 -0.0005
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.00025
17

@tshoemaker Thank you for the update. Could you try a smaller value of k_na, say 13,000 or 1,300? Maybe the attractive forces plus gravity are a bit too strong relative the repulsive forces.

18

@ndtrung, thank you for your response.

I can certainly try a different value of k_na and report my findings. However, to my understanding, the SI example should not have any attractive forces. The cutoff distance is set to zero:
variable cut_inner equal 0.0
pair_style body/rounded/polyhedron ${c_n} ${c_t} ${mu} ${A_ua} ${cut_inner}

To my understanding, the body/rounded/polyhedron pair style has no attractive forces with this setting. From the LAMMPS user manual:

This is the distance at which two particles no longer interact. If rc is specified as 0.0, then it is a contact-only interaction.

Please let me know if I am misinterpreting the user manual. Thank you.

https://docs.lammps.org/pair_body_rounded_polyhedron.html

19

You are correct, sorry I missed that cut_inner value. I would suggest further reducing c_n and c_nt, and setting mu to zero, to see how the repulsive forces are strong enough to remove overlaps when the bodies start to be at contact.

20

@ndtrung, I tried setting the c_n, c_nt, and mu to zero. The particles continue to pass through each other. I also experimented with different values of k.

image