Building superellipse with region command

LAMMPS LAMMPS Development
1

Dear Users,

I want to restrict my particles inside a superelliptic region. More specifically, I want to create a region described by the superellipse equation (2D) with an exponent of 0.5.

To do this, I divide my square simulation box into 4 sub-squares, then use the four corners of my simulation box as a center and draw a circle of radius half of the simulation box. Each of these circle regions is defined with a “side out” command. Now, the “intersect” of each sub-square with the corresponding circles will give me four arc regions. The “union” of this four arc region with the “side in” command should give me the superellipse region of exponent 0.5 that I desire. Then I fix a wall at this newly defined superelliptic region.

For simplicity, I placed two particles at the center of the box and ran the simulation to check the system. The system ran a few timesteps and then I kept getting the “Particle outside surface of region used in fix wall/region” error. Are the walls not being “seen” by the particles or have I overlooked the functionality of some of the commands? Please suggest any changes.

Below are my input files.

###############################################

LAMMPS input script

variable r index 1

units lj
atom_style ellipsoid
dimension 2
boundary f f p
read_data config.txt

variable T index 1.2
variable E index 3

region circle1 ellipsoid 0 0 0 200 200 0.5 units box side out
region box1 block 0 200 0 200 -0.5 0.5 units box side in
region arc1 intersect 2 circle1 box1 side in units box

region circle2 ellipsoid 400 0 0 200 200 0.5 units box side out
region box2 block 200 400 0 200 -0.5 0.5 units box side in
region arc2 intersect 2 circle2 box2 side in units box

region circle3 ellipsoid 400 400 0 200 200 0.5 units box side out
region box3 block 200 400 200 400 -0.5 0.5 units box side in
region arc3 intersect 2 circle3 box3 side in units box

region circle4 ellipsoid 0 400 0 200 200 0.5 units box side out
region box4 block 0 200 200 400 -0.5 0.5 units box side in
region arc4 intersect 2 circle4 box4 side in units box

region superellipse union 4 arc1 arc2 arc3 arc4 side in units box

velocity all create $T 87287 loop geom

pair_style gayberne 1.0 3.0 1.0 4.0
pair_coeff * * 1.0 1.0 $E $E 1 $E $E 1 4.0

neighbor 1.9 bin
neigh_modify delay 0 every 1 check yes

fix wall all wall/region superellipse lj93 1.0 1.0 2.5

fix 1 all nvt/asphere temp $T $T 5
fix 2 all enforce2d

compute q all property/atom quatw quati quatj quatk
compute sh all property/atom shapex shapey shapez

thermo_style custom step temp press pe ke epair etotal density
thermo 1000
log log.rho0.$r

dump 1 all custom 10 dump$r.patch id type x y z c_q[1] c_q[2] c_q[3] c_q[4] c_sh[1] c_sh[2] c_sh[3]
dump_modify 1 sort id

timestep 0.001
restart 1000000 restart$r
run 1000000

###############################################

configuration file
config.txt
###############################################

LAMMPS data file

2 atoms
2 ellipsoids

2 atom types

00 400.0 xlo xhi
00 400.0 ylo yhi
-0.5 0.5 zlo zhi

Atoms

1 1 1 0.4 199.54 200.1581033398804 0.0
2 2 1 0.4 200.965 199.33537920628518 0.0

Velocities

1 0.0 0.0 0.0 0.0 0.0 0.0
2 0.0 0.0 0.0 0.0 0.0 0.0

Ellipsoids

1 3.0 1.0 1.0 1.0 0.0 0.0 0.0
2 3.0 1.0 1.0 1.0 0.0 0.0 0.0

################################################

2

This is too complex a region for the fix to work. There are too many cases, where it is ambiguous which sub-region has to be applied to a particle.

The best approach is likely to build a wall from explicit particles instead.

3

Many thanks for your reply.
I will check that method.

Also, I tried modifying the region ellipsoid source code to accommodate other values of the superellipse exponent. As long as the exponent is an integer and greater than 1, the system/simulation seems to work. I cannot run a denser system when the exponent is 1. The system runs for a few timesteps, and then I get a lost particle error. And, when the exponent is not an integer (0.5, 1.5, 2.5, etc.), I cannot start the simulation and I get the particle out of the fix wall region error even though the particle is in the middle.

Would it be okay if you take a look at this modified source code and possibly identify issues that might arise when integrating this modified source code with the rest of the lammps code?