[lammps-users] miscalculated coordination

LAMMPS LAMMPS Mailing List Mirror
1

Hello,

I’m stretching a sheet of 2d silica until breaking, and I’m counting the broken bonds during stretching through the coord/atom compute.

However, when I write the dump and check the coordination number, sometimes I found that it’s wrongly computed.

In brief, type=1 for Si, type=2 for O, and c_uncoord is a scalar defined by:

compute coord all coord/atom cutoff 2.016
variable uncoordination atom abs(c_coord-4)*(type==1)+abs(c_coord-2)*(type==2)
compute uncoord all reduce sum v_uncoordination
variable uncoord equal c_uncoord

The attached script+input produces a file sio2_3_bondbreak_10.dump, with the following “uncoordinated” atoms:

#id type c_coord
226 1 3
256 1 3
932 2 0
1241 2 4
1666 1 3
1696 1 3
2372 2 0

But looking with ovito at oxygen atom 1241, it shows that it has only two Si neighbors (not 4) at a distance smaller than the cutoff above.
So, the c_coord is somehow miscalculated. I have encountered this problem using many different initial samples.

Also, from the log it appears that c_uncoord, along with fmax and other thermo quantities, change suddenly after a write_dump:

Step Temp PotEng v_pe Fmax Press Lx v_strain c_stressx v_pressx v_uncoord
478000 0.0054869426 -22512.627 -22512.627 0.060812995 -3556.1136 97.878575 0.1255 1.2360716e+10 7093.2647 10
479000 0.0039295438 -22512.631 -22512.631 0.0098276162 -3556.4664 97.878575 0.1255 1.2361198e+10 7093.5409 10
Loop time of 0.625285 on 16 procs for 1000 steps with 3456 atoms
...
write_dump all custom sio2_3_bondbreak_10.dump id type xu yu zu c_coord c_peratom c_stratom[1] q
run 0
...
Step Temp PotEng v_pe Fmax Press Lx v_strain c_stressx v_pressx v_uncoord
479000 0.0039295438 -22512.631 -22512.631 0.0098227024 -3556.4664 97.878575 0.1255 1.2361198e+10 7093.5409 8
Loop time of 2.85619e-06 on 16 procs for 0 steps with 3456 atoms

Timestep is the same, 479000, but values differ.

Is there something that I’m missing?

Thanks for any reply,
Roberto

config.data (176 KB)

silica_2D_stretch.lmp (7.27 KB)

2

Try building the neighbor list more often. One thing that happens after a restart is building a new neighbor list.

Axel

3

Hello Axel,

indeed in the script I have put

neigh_modify every 1 delay 0 check yes

but it does not help.

RG

4

You can try to switch to using “check no” to force rebuilding every step.

5

Dear Axel,

I have tried with ‘check no’, and it seems to help for the coordination (at least for this specific sample; I’ll try with others although it has become painfully slow).

I cannot understand why ‘check yes’ and ‘check no’ should produce different results: Tersoff cutoff is much larger than coordination cutoff… plus we have the skin.

However, even with ‘check no’ some quantities such as fmax still differ for the same timestep:

Step Temp PotEng v_pe Fmax Press Lx v_strain c_stressx v_pressx v_uncoord
1803000 5.8706988e-05 -22303.435 -22303.435 0.00043827858 -4117.8913 100.31358 0.1535 1.4941997e+10 8366.4085
1804000 3.4199686e-05 -22303.435 -22303.435 0.0064266143 -4117.4425 100.31358 0.1535 1.4940411e+10 8365.5203 2
Loop time of 0.929886 on 16 procs for 1000 steps with 3456 atoms
…
run 0
…
Step Temp PotEng v_pe Fmax Press Lx v_strain c_stressx v_pressx v_uncoord
1804000 3.4199686e-05 -22303.435 -22303.435 0.006423401 -4117.4425 100.31358 0.1535 1.4940411e+10 8365.5203 2
Loop time of 2.70893e-06 on 16 procs for 0 steps with 3456 atoms

I’m running on version 3Mar20.

This sounds like a bug to me, but maybe you have an explanation.

Thank you and best regards,
Roberto

6

check yes will not update the neighbor list unless one atom has moved at least half the neighbor skin distance. if you are deforming your system, one could conceive cases where the neighbor list would need updating, but no atom has moved sufficiently far. please note that the communication cutoff does not impact the neighborlist cutoff, it merely makes certain that a sufficient amount of ghost atoms are present (it must be at least the largest cutoff plus skin).

if it is a known bug, you would have to try with the latest LAMMPS release (29 October 2020) to get the bugfix.

whether the difference in Fmax is genuine or not is very difficult to say since it is somewhat small. if there is a significant amount of cancellation happening (i.e. forces between specific pairs of atoms are very large but canceling each other). when restarting or reneighboring the order of atoms changes and any sum (like forces) over numbers with opposite sign and varying magnitude will depend on the order of the summation due to the non-associative behavior of floating point addition. to look into this more closely, a small test example would be required. with potential functions that do not go smoothly to zero at the cutoff or that have other properties that can change the energy/force significantly on small changes, sometimes such changes may occur. but looking at our unit test “epsilon” parameters, it seems the plain tersoff potential is not so no noisy. all of its variants are more noisy by up to two orders of magnitude.

axel.

7

I agree that cell deformation can affect the neighbors. The coordination issue must be due to that. Thank you.

About Fmax I’m more skeptical to address it as roundoff error for two reasons. First, differences are too large, e.g. I see a 0.00026 becoming 0.00037 after a ‘run 0’. It’s a 1e-4 eV/Ang change… very large. And the number of neighbors is quite small in Tersoff and it has a smooth cutoff.

Second, it occurs also in the absence of cell deformation, during a simple damped dynamics; so it should not be related to the other issue above.

Thank you and best regards,
Roberto