Dynamic Addition of Atoms in LAMMPS Beam Model Causes Overlap and Excessive Force

LAMMPS LAMMPS General Discussion
,
1

Hello Everyone ! Firstly, thank you for taking the time to help out the LAMMPS community !

I am currently facing an issue in my simulation, which requires me to dynamically add small atoms (type 2) next to already existing larger atoms (type 1), and have them grow until they obtain the same size as the existing atoms (type 1).

Of course, when atoms are added dynamically and their diameter grows with every run, the bond radius will also adapt dynamically.

The issue I am facing is that, when I create the new type 2 atoms, if I set their diameters initially (set type 2 diameter 0.001), I experience an overlap with the type 1 atoms, even if I had already mentioned a pair style and pair coefficient between the two types. What is even more interesting is, if I do not set the diameter initially, the overlap does not exist.

Below are the input files I use to run this simulation. If I want to initialize diameter or not, I comment out the ‘set type 2 diameter 0.001’ line.

Beam

create_atoms 1 single 0.000 0.0 0.0
create_atoms 1 single 0.005 0.0 0.0
.
.
.
create_atoms 1 single 0.045 0.0 0.0

Bonds

create_bonds single/bond 1 1 2
.
.
.
create_bonds single/bond 1 9 10

Angles

create_bonds single/angle 1 1 2 3
.
.
.
create_bonds single/angle 1 8 9 10

Input File

.
.
.


# create growing atoms

create_atoms 2 single ${new_x1} ${new_y1} 0.0
create_atoms 2 single ${new_x2} ${new_y2} 0.0
set type 2 density 0.5
set type 2 diameter 0.001

Pair Interaction

pair_style hybrid/overlay granular lj/cut 0.1
pair_coeff * * granular hertz/material 1e5 0 0.5 tangential linear_nohistory 0 0
pair_coeff * * lj/cut 0 0 0
special_bonds lj/coul 0 1.0 1.0

Growth Loop

label loop_a
variable a loop 2500


# rebuild topology

run 0 pre yes post yes

# diameter growth rate
variable time_growth equal "0.001 + (0.004) * ((v_a) / (2500))"
variable r0 equal "0.003 + (0.005-0.003)* ((v_a) / (2500))"

bond_coeff ${lb} 2.1618112757330144 ${r0}
bond_coeff ${rb} 2.1618112757330144 ${r0}

group grow_atoms id ${l} ${r}
fix grow grow_atoms adapt 1 atom diameter v_time_growth


run 1
next a
jump SELF loop_a

Now, with gravity off, if the set type 2 diameter code is commented out, the new atoms I create do not overlap with the existing type 1 atoms. But, at timestep 0, I start of with very larger type 2 atoms, because LAMMPS randomly assigns one to it. However, if the ‘set type 2 diameter’ is defined, then there is an overlap. By overlap, I mean that the growing type 2 atoms invade the personal space of the type 1 atoms. That is the basis of my problem.

The end goal is to grow these atoms, without any overlap, while gravity is turned on.

2

Your post is very hard to read, since you don’t use the correct way to quote. To learn how to do this and also how to properly ask questions and what information to provide, please read the post with guidelines and suggestions for this forum and adjust/edit your post accordingly.

Please also note that it will be much easier to get help, if you post an input that has just a minimum set of commands that highlight the specifics of your problem rather than a large and complex input deck.

3

Hi @akohlmey,

Sorry for the messy post - I cleared it up now! I hope it is easier to read now.

4

Not much. It is still impossible to reconstruct what you are doing. Too much information is available only in fragments. It is not clear at all where your problems are.

Furthermore, you don’t report which LAMMPS version you are using.

Please keep in mind that we cannot read your mind and are not familiar with your research.
So please do not make any assumptions about what people would know.
Furthermore, please keep in mind that the chance of getting help becomes larger, the easier it is to reproduce what bothers you and you explain in detail why and how this is.

5

@akohlmey , I understand - sorry for the back and forth. I will try to break down what I am trying to achieve, and where I am facing an issue. I am using LAMMPS (19 Nov 2024).

Essentially, I am trying to create a simulation that models bacteria growth. I start with a rod of atoms (type 1), and have them grow by dynamically adding (type 2) atoms, then I have this rod split (delete bond), and I now have two rods. I repeat this process for a while, and model bacteria repopulation - starting with 1 rod, and ending with, say, 100 rods.

In reality, this simulation takes places in some sort of chamber, in which fluid flows in and out, and bacteria consume nutrients to help them grow. This is a little more complex, so for now, I plan to model this ‘fluid flow’ as gravity.

So now, I have gravity in my simulation, as well as 2 rates I identify: one for diameter growth (of the type 2 atoms), and one for the bond cutoff distance.

variable time_growth equal "0.001 + (0.004) * ((v_a) / (2500))"
variable r0 equal "0.003 + (0.005-0.003)* ((v_a) / (2500))"

The type 2 diameters start off with a diameter of 0.001, and grow until their diameter is 0.005, same as those of type 1. If the atoms of two different types are exactly touching, their cutoff distance is the sum of their radii (0.0005 + 0.0025 = 0.003), which is what r0 takes care of.

So where does the issue arise? Well first, I had tested this ‘growth toolkit’ with no gravity turned on. When I dynamically added atoms of type 2, I had specified a density, but never a diameter value, sort of like this:

create_atoms 2 single ${new_x1} ${new_y1} 0.0
create_atoms 2 single ${new_x2} ${new_y2} 0.0
set type 2 density 0.5

The positions of the new atoms were founded by extracting the ‘end atoms’ of the original rod, and then adding a net distance of 0.003 to them. This is so that the new created atoms are exactly touching the existing ones:

# endpoints

if "${i} == 1" then "variable x_one equal x[1]"
if "${i} == 1" then "variable y_one equal y[1]"

if "${i} == 1" then "variable x_two equal x[${m}]"
if "${i} == 1" then "variable y_two equal y[${m}]"

if "${i} != 1" then "variable x_one equal x[${lbefore}]"
if "${i} != 1" then "variable y_one equal y[${lbefore}]"

if "${i} != 1" then "variable x_two equal x[${rbefore}]"
if "${i} != 1" then "variable y_two equal y[${rbefore}]"


# direction and unit vector

variable dx equal "v_x_two - v_x_one"
variable dy equal "v_y_two - v_y_one"

# new atom distance

variable d equal 0.003

# new atom positions

variable theta equal "atan2(v_dy,v_dx)"

variable new_x1 equal "v_x_one - v_d * cos(v_theta)"
variable new_y1 equal "v_y_one - v_d * sin(v_theta)"
variable new_x2 equal "v_x_two + v_d * cos(v_theta)"
variable new_y2 equal "v_y_two + v_d * sin(v_theta)"

Now, if I do not specify the type 2 atom diameter, my first timestep creates atoms that nearly cover the whole simulation box (radius = 0.5). From the research I have done, this is simply LAMMPS randomly assigning a value, so that the created atom is not a point particle. And so, as the simulation continues, I see the growing atom (type 2) increase in diameter, next to one of the initial endpoints of the rod, WITHOUT overlapping. The growing atom (Type 2) pushes into the Type 1 atom, understands that it exists, and bounces back (very small displacement), and essentially ‘understands’ that it should stay next to the Type 1 atom, while growing, without invading it’s personal space.

When gravity is turned on, however, the very first timestep of having super larger type 2 atoms instantly breaks the simulation, and so, the rod is disbanded, and atoms are all over the place.

So, when I limit the initial size by setting an initial diameter (0.001), and gravity is turned OFF, you would expect the same result, right? At least I did. However, this time, with only one line of code added (set type 2 diameter 0.001), now the growing atom does not respect the already existing endpoint of the rod (type 1). What happens, as a result, is overlapping, with the growing atom invading the personal space of the existing atom. Now, this still works with gravity turned on, but of course, the reaction forces pile up, and the system gains a lot of energy after a certain amount of time.

The main question I am trying to answer is: why is initially defining a diameter change the entire outcome? Whether it starts of as an atom with diameter of 0.001, or 0.5, it enter a growth loop which dynamically changes the diameter as a function of the timestep. So why, in one case, does it overlap, but in the other, it does not.

I know this is long, but I hope it was clearer in what I was trying to achieve. Again, I do not want to supply all the code, in the event that it is redundant, but I am happy to provide whatever might be useful. Lastly, thank you again, for your time and effort - I sincerely appreciate it.

Best,
aliesk

6

You are trying to explain your research, but I am not a domain expert and thus am not even going to make an attempt to understand that. So you are wasting your and my time.

Instead, all I need are two simple but complete(!) input decks with the absolute minimum amount of commands and the absolute minimum of differences between them that demonstrate what works as expected and what does not. Simple as that.