Strange behaviour, error and possible bug with fix wall/region with ellipsoid region

Dear All,

I am experiencing strange behaviour when using “fix wall/region” with a region defined using “style ellipsoid”. Particles go through the wall and this results in an error,

“ERROR on proc 0: Particle outside surface of region used in fix wall/region (…/fix_wall_region.cpp:307)”

I made a minimal input file to test this (in.test1 and corresponding log.test1). It consists of a single particle starting in the centre of the box, surrounded by a repulsive (WCA) ellipsoid. Instead of being repelled by the wall, it goes through and generates an error.

As a sanity check, I then defined a spherical region using the ellipsoid region style by setting a=b=c (in.test2).

region elip ellipsoid 0 0 0 $((v_r)) $((v_r)) $((v_r)) side in

The same thing happens and it generates the same error (log.test2_elip). However, when I define the same spherical region using the sphere region style,

region S sphere 0 0 0 $((v_r)) side in

it works as expected with no errors (in.test2 and log.test2_sphere). So it seems that there is something off with the “fix wall/region” and “region style ellipsoid” combination.

Please see the attached input/log files for reference. Note that species 2 is just there to help visualize the regions. I am using LAMMPS version 2Aug2023.

Can anyone explain this behaviour (and ideally offer a solution)? Any help would be much appreciated!

EDIT: I cannot upload files because I’m a new user. Instead I’ve appended in.test2 to the end of this post. Happy to provide files via email.

Kind regards,

James


units           lj
atom_style      atomic
dimension       3

# sphere radius
variable r equal 8

## box length
variable L equal 21

region box0 block $((-v_L/2)) $((v_L/2)) $((-v_L/2)) $((v_L/2)) $((-v_L/2)) $((v_L/2)) units box
create_box 2 box0

## same regions defined via style ellipsoid and sphere
region elip ellipsoid 0 0 0 $((v_r)) $((v_r)) $((v_r))  side in
region S sphere 0 0 0 $((v_r)) side in

create_atoms 2 random 1000 1818 elip ## just to help visualize ellipsoid
create_atoms 1 single 0 0 0

group 1 type 1

# non-interacting WCA potentials for everything
pair_style lj/cut $((2^(1/6)))
pair_coeff * * 0.0 1.0
mass * 1.0

timestep 0.002

fix wall all wall/region elip lj126 1.0 1.0 $((2^(1/6))) ## does not work
#fix wall all wall/region S lj126 1.0 1.0 $((2^(1/6))) ## works

velocity 1 set 1 1 0

fix nve 1 nve

thermo_style custom step temp pe etotal press vol density
thermo 100

dump traj all custom 100 dump.traj id type x y z

run 100000

You need to quote the sections with triple backquotes (```) or else they will not render properly as the forum software recognizes some special characters as markup.

I have edited the original post so that it displays correctly. Thanks!

Thanks. I have looked into the example and did some debugging. There are definitely some bugs in the specific parts of the ellipse region style here. I have identified one, but that does not fix the issue at hand.

Would you mind to create a bug report issue on GitHub for this? That would make it easy to notify the contributor of this piece of code and I can also add my findings there.

Opened a bug report. Thanks!

@JamesS Please see Fixes #3925 in region_ellipsoid.cpp by evoyiatzis · Pull Request #3931 · lammps/lammps · GitHub for a bugfix.

This has fixed the bug for me. Thank you and evoyiatzis for your help!