Anisotropy in liggghts with bonds-redux

Submitted by tjleps on Tue, 02/19/2019 - 21:53

Sorry to repeat post, however after considerable time I've been unable to find any reason why this should be happening, and am curious if anyone else who uses LIGGGHTS with Bonds has been able to reproduce the problem.

When forming bonds between particles in a random close packing or a random disperse distribution, I find bonds are preferentially formed in the x and y directions to the z direction. This occurs both in the initial, random, insert/pack and if I allow the particles to settle under gravity into a container and then bond them (regardless of the direction of gravity). It does not seem to be a problem when fully bonding lattice atoms, for example all 12, nearest neighbor, FCC bonds form without preference to direction. This behavior is extremely troubling for me because
a)It makes my own research challenging. but also
b)Because the only place I can imagine the problem is arising from must be in the neighboring process since there is not much in the actual bonding algorithm to break, and that could have implications throughout the simulation physics.

I know that many people on the forums are much more skilled with "looking under the hood" in LIGGGHTS, so maybe someone will have better luck than me figuring out what's going on.

*Edit
https://www.cfdem.com/forums/anisotropy-liggghts-bonds
I've linked to the previous post I made, which shows some of the behavior from a simulated tensile test using a bonded material. When pulled in the Z direction there are striations where it appears no bonds are made across the XY-plane in large sections of the bulk material. When rotated around the Y-axis 90 degrees and pulled in the X direction, there are no corresponding gaps.

Forums:

mschramm | Tue, 02/19/2019 - 22:18

Hello,
Could you provide a simple script for me?
I use this for my own research but I model long fibers (wheat straw).
I have written matlab code to tell me the average orientation of bonded particles.
When dealing with just formed particles (before settling) I observed nearly uniform bonding in the x-y direction while
the z-direction had most bonds near horizontal that declined slowly to 45 degrees then decreased drastically to vertical orientations.
This was in a 4 inch shear cell so i expected to have more horizontal fibers than vertical ones.

tjleps | Wed, 02/20/2019 - 01:50

Sure, I'll write one up that's fairly stand alone, so you don't have to parse through some of the other stuff I've done trying to isolate it. I also have a script that finds the average bond angle from normal to the plane bounding a half-space, and I find that the average bond angle is consistently much larger than 1 radian using the z-hat vector as the angle reference (corresponding to fewer bonds than expected in the z direction) and smaller when I use x-hat or y-hat. For random bond angles, the average bond angle in a half-space should be 1rad.

I am using Richti's bonds implementation, though I did try with your implementation as well and got qualitatively similar results (I haven't verified that they're completely identical since I increased the rigor of my testing, after I verified the problem existed in both implementations, but using only Richti's implementation).

Daniel Queteschiner | Wed, 02/20/2019 - 16:01

What is the maximum number of neighbours in your simulation and what is the maximum number of bonds per atom you allow to be created?
Without having looked into the code too much, I could imagine that
* if your max. bond number is smaller than the number of neighbours and
* if neighbours are sorted from min z to max z (from the binning procedure in the neighbour list generation)
there may be less bonds created in z-direction.
At least that's something I would check ...

tjleps | Wed, 02/20/2019 - 18:55

I used between 3 and 12 bonds per atom during my testing, and had large neighbor lengths ~15-20 neighbors per atom to try to eliminate these possibilities. I suspect however that you're correct about it being an issue with the neighbor sorting. I was under the impression that it was binning from nearest to furthest neighbor, but perhaps something else is going on. I'll try to post a simple diagnostic script to see if other people are able to recreate the conditions I am seeing.

tjleps | Thu, 02/21/2019 - 19:44

In this script I get an average bond angle of:
X 0.8496rad
Y 1.0095rad
Z 1.1497rad

Which is fairly consistent, qualitatively, between every measure I can think of, ie. different seed values, bond numbers, packing densities, simulation space ect.

atom_style hybrid granular bond/gran n_bondtypes 1 bonds_per_atom 6
#atom_style granular
atom_modify map array
boundary f f f
newton off

units si
communicate single vel yes
processors 1 1 1

region reg block -.00515 .00515 -.00515 .00515 -0.00515 .00515 units box
create_box 1 reg

neighbor 0.0006 bin
neigh_modify delay 0 check yes

pair_style gran model hertz tangential history
pair_coeff * *
bond_style gran
bond_coeff 1 1 21e8 21e8 1 1e40 1e40

fix m1 all property/global youngsModulus peratomtype 1.e7
fix m2 all property/global poissonsRatio peratomtype 0.3
fix m3 all property/global coefficientRestitution peratomtypepair 1 0.5
fix m4 all property/global coefficientFriction peratomtypepair 1 0.5

group nve_group region reg

fix pts1 all particletemplate/sphere 160481219 atom_type 1 density constant 4228.75 radius constant 0.00015
fix pdd all particledistribution/discrete 141650963 1 pts1 1.0
fix ins nve_group insert/pack seed 160481183 distributiontemplate pdd &
insert_every once maxattempt 300 overlapcheck yes volumefraction_region 0.5 region &
reg ntry_mc 10000

timestep 0.0000015

run 1
write_data RectBoxUnsettled.data
fix bondcr all bond/create/gran 1 1 1 0.0006 1 6 prob 1 6786797
write_data RectBoxUnsettledBonded.data
dump dmp all custom 1 dump*.liggghts id type x y z

tjleps | Thu, 02/21/2019 - 19:56

This script is written for Richti's "LIGGGHTS with bonds", however the same problem occurs in mschramm's implementation as well, and the script should be quickly adaptable between the two. Thanks for any help anyone is able to provide!

mschramm | Fri, 02/22/2019 - 00:33

Looking at what you provided (with some alterations) and 144180 bonds...
I also allowed upto 13 bonds per atom (should only ever be 12 bonds per atom...)
xy_angle was shown to follow a uniform distribution using kstest with a p-value > 0.05 while
the xz_angle was no where near uniform (p-value ~= 0.0) and showed traits of a beta function...

in contrast, I ran my uni-axial compression tutorial (from github) and got a uniform distribution for the x-y and a near uniform distribution for x-z
This is the difference between pre-generated bonded particles getting randomly placed and oriented vs particles being randomly placed and then forming a bond.

I also noticed that not all spheres were getting all neighbors possible

tjleps | Fri, 02/22/2019 - 01:49

It's a very strange result, I couldn't find anything, anywhere, that should lead to any preference in one direction or another during the bonding process. It's also worth noting, though I don't think I mentioned it before, that there is a bias for bonding in the X direction, that's nearly as strong as the bias against bonds forming in the Z direction, though I find that to have less of an impact on bulk parameters.

tjleps | Fri, 02/22/2019 - 05:47

I see, if I'm reading correctly, that you were actually getting an anomaly in the only direction I saw the correct average bond angle, angle from the xz plane. Am I interpreting that correctly?

Daniel Queteschiner | Fri, 02/22/2019 - 13:57

Your packing produced by fix insert/pack has a rather low volume fraction (~0.33). That means that there is quite a bit of spacing between the particles which makes it more difficult to guess which particle pairs actually require bonding. (If they where close packed you can go with the actual neighbours and the kissing number is your maximum bond number).
You specify a very large skin distance of 0.0006 (twice the particle diameter!) which gets added to the particle diameter to form the radius of Verlet list.
As you can see from the simulations stats at the end of the run, this leads to an average of about 33 neighbours per particle. With your bond creation settings you can set the skin distance as low as 0.0003, which leads to an average of ~9.5 neighbours per particle. (That's still quite large/a lot of neighbours). Since this is only half of the neighbours, you would need to allow for ~20 bonds per atom to potentially allow bonding to all of them. If you do so, you will end up with about 250000 bonds for ~25000 particles (probably also forming bonds with next nearest neighbours). I would assume that any anisotropy is negligible then.

Just to illustrate the creation of bonds for atom 15870 in this case (max 25 bonds)

creating bond btw atoms 15870 and 15099 (z=0.000753, bin=(20|32|22)) (i has now 1 bonds)
creating bond btw atoms 15870 and 15100 (z=0.000831, bin=(21|32|22)) (i has now 2 bonds)
creating bond btw atoms 15870 and 15101 (z=0.000711, bin=(22|32|22)) (i has now 3 bonds)
creating bond btw atoms 15870 and 15123 (z=0.000837, bin=(23|33|22)) (i has now 4 bonds)
creating bond btw atoms 15870 and 15828 (z=0.001041, bin=(20|30|23)) (i has now 5 bonds)
creating bond btw atoms 15870 and 15830 (z=0.001147, bin=(22|30|23)) (i has now 6 bonds)
creating bond btw atoms 15870 and 15852 (z=0.000918, bin=(21|31|23)) (i has now 7 bonds)
creating bond btw atoms 15870 and 15853 (z=0.000990, bin=(22|31|23)) (i has now 8 bonds)
creating bond btw atoms 15870 and 15869 (z=0.001059, bin=(20|32|23)) (i has now 9 bonds)
creating bond btw atoms 15870 and 15892 (z=0.000996, bin=(22|33|23)) (i has now 10 bonds)
creating bond btw atoms 15870 and 15913 (z=0.001079, bin=(21|34|23)) (i has now 11 bonds)
creating bond btw atoms 15870 and 16592 (z=0.001417, bin=(21|31|24)) (i has now 12 bonds)
creating bond btw atoms 15870 and 16593 (z=0.001507, bin=(22|31|24)) (i has now 13 bonds)
creating bond btw atoms 15870 and 16616 (z=0.001346, bin=(20|32|24)) (i has now 14 bonds)
creating bond btw atoms 15870 and 16617 (z=0.001445, bin=(21|32|24)) (i has now 15 bonds)
creating bond btw atoms 15870 and 16618 (z=0.001259, bin=(22|32|24)) (i has now 16 bonds)
creating bond btw atoms 15870 and 16641 (z=0.001276, bin=(20|33|24)) (i has now 17 bonds)
creating bond btw atoms 15870 and 16642 (z=0.001284, bin=(22|33|24)) (i has now 18 bonds)
creating bond btw atoms 15870 and 16664 (z=0.001462, bin=(21|34|24)) (i has now 19 bonds)
creating bond btw atoms 15870 and 17390 (z=0.001588, bin=(21|32|25)) (i has now 20 bonds)
creating bond btw atoms 15870 and 17410 (z=0.001557, bin=(20|33|25)) (i has now 21 bonds)
creating bond btw atoms 15870 and 17411 (z=0.001623, bin=(22|33|25)) (i has now 22 bonds)

You see that the neighbours are indeed handled according to their bins with the top bins (zbin=25) coming last.

With a max of 12 bonds:

creating bond btw atoms 15870 and 15123 (z=0.000837, bin=(23|33|22)) (i has now 1 bonds)
creating bond btw atoms 15870 and 15828 (z=0.001041, bin=(20|30|23)) (i has now 2 bonds)
creating bond btw atoms 15870 and 15830 (z=0.001147, bin=(22|30|23)) (i has now 3 bonds)
creating bond btw atoms 15870 and 15852 (z=0.000918, bin=(21|31|23)) (i has now 4 bonds)
creating bond btw atoms 15870 and 15853 (z=0.000990, bin=(22|31|23)) (i has now 5 bonds)
creating bond btw atoms 15870 and 15869 (z=0.001059, bin=(20|32|23)) (i has now 6 bonds)
creating bond btw atoms 15870 and 15892 (z=0.000996, bin=(22|33|23)) (i has now 7 bonds)
creating bond btw atoms 15870 and 15913 (z=0.001079, bin=(21|34|23)) (i has now 8 bonds)
creating bond btw atoms 15870 and 16592 (z=0.001417, bin=(21|31|24)) (i has now 9 bonds)
creating bond btw atoms 15870 and 16593 (z=0.001507, bin=(22|31|24)) (i has now 10 bonds)
creating bond btw atoms 15870 and 16616 (z=0.001346, bin=(20|32|24)) (i has now 11 bonds)
creating bond btw atoms 15870 and 16617 (z=0.001445, bin=(21|32|24)) (i has now 12 bonds)

You see that it cannot form bonds with most of the particles in the lowest zbins since those have no free bonds anymore and skips any potential bonds in the top zbins since the max number of bonds have already been created

tjleps | Fri, 02/22/2019 - 17:36

I guess I misunderstood how neighbor lists were being formed, I didn't realize they form from the top down as opposed to in to out. I've had the same problem with a wide range of neighbor lengths, however. I used a large length in this script based on your earlier comment about it potentially being a problem with insufficient neighbor lengths.

Is there a way to randomize the neighbor lists to eliminate the binning bias? For my application I don't want a dense packing of particles before bonding, and would like to be able to vary the number of bonds independently of other variables. Also if I drop the neighbor length so that the average number of neighbors is close to the number of bonds I'm attempting to create, I'd be concerned that I'd end up with a large number of particles that would not have enough neighbors.

Thanks

Daniel Queteschiner | Fri, 02/22/2019 - 18:15

Neighbor lists in LIGGGHTS are a combination of cell lists and Verlet lists:
First the domain is subdivided into equal bins (based on the particle size). It is straight forward to map the particle coordinates to these bins and at this point we only need to search for potential neighbors in adjacent bins (typically running over x-, then y- and finally z direction). To further reduce the number of neighbors (as we want to avoid unnecessary intersection tests and potential force calculations) any particles that are further away than the particle diameter+skin are excluded from the neighbor list. (The 'skin' addition is added to avoid the costly rebuild of the neighbor list each time the particles move; only when a particle has moved further than this distance a rebuild is triggered).

>>I used a large length in this script based on your earlier comment about it potentially being a problem with insufficient neighbor lengths.
What I was actually referring to was the average number of neighbors per atoms and the number of allowed bonds. Of course you want to keep the former low to only have the nearest neighbors while the latter needs to be adjusted to the number of neighbors to establish all bonds. You can tune this interdependence to some extent by restricting the length at which bonds are formed (i.e. the parameter in fix bond/create/gran)

>>Is there a way to randomize the neighbor lists to eliminate the binning bias?
Somewhat difficult without coding, you may adjust your settings so that all potential bonds can be formed and then use the probability option of fix bond/create/gran to reduce the number of bonds to what you actually want ...
A proper solution would be to change the creation code such that it sorts the neighbors by distance (just for the bonding operation of course) and forms bonds between the nearest ones.

tjleps | Fri, 02/22/2019 - 19:00

I'll see about sorting the neighbor list during the bonding operation. I've programmed a few fixes, so I think I should be able to get it figured out.

Thanks for your help!

mschramm | Fri, 02/22/2019 - 18:24

Hello,
Ran some more tests with 25 possible bonds per atom.
Max bond length was 3 times sphere radius.

I also made it so the particles entered a box. I allowed the particles to settle while continuing to cram as many spheres into the box as possible (40000 steps with 27914 atoms and 172216 bonds)
The maximum number of bonds was found to be 13 bonds per atom.
Again the angle of the bonds in the xy plane appeared to be uniform while
the angle between the xy plane and the z direction appeared to be distributed "in a triangle shape..."

tjleps | Fri, 02/22/2019 - 19:05

MSchramm
It seems the problem is indeed in how the neighbor lists are built, they are ordered from the negative z bins to the positive z bins, and so the bonds are only able to form horizontally if there are many more neighbors than bonds created.

I'll have to take a look at your algorithm for forming bonds then placing particles, perhaps that will be a good work around until I have a chance to resort the neighbor list by distance during the bonding process.

Thank you for your help!