Output from Voronoi varies with number of processors

Dear LAMMPs users,

I am working on external irradiation of Fe. The projectile has an energy of 10 keV and a velocity vector along the z axis.

Since there can be atoms that are sputtered from the surface, I use f boundary conditions on z (advice from Axel). Box size = 64x64*96

In order to calculate the number of vacancies and interstitials that are formed during the cascade, I use the Voronoi package with voronoi/atom occupation command. The output is then postprocessed with a Python script that I wrote (I tested it in other cases and I believe it is ok).

I observed that depending on the number of processors, the number of vacancies and interstitials can strongly vary.

For 24 procs I obtain around 200 V and 100 SIAs while for 12 processors, I obtain around 1500 V and 800 SIAs. Very different.

  • Is there something wrong in my input or is it an issue with Voronoi package ?

  • Since some atoms can be sputtered out of the crystal, is it correct to use voronoi/atom occupation on group ‘crystal’ or should I use ‘all’ ?

Many thanks in advance for any help you could bring me.
Christophe

PS: I also observed that when I use m boundary condition on z axis, Voronoi can give very large number (eg 1e8), which I believe is a bug.

dimension 3
units metal
atom_style atomic
atom_modify map hash

boundary p p f
lattice bcc 2.8553

include in.setup_variables

Definition of simulation domain

region simdom block 0. {ncellx} 0. {ncelly} -5. ${ncellz}

Definition of crystal region

region crystal block 0. {ncellx} 0. {ncelly} 0. ${ncellz}

Definition of central atom region

region cent_atoms block {xminthermo} {xmaxthermo} {yminthermo} {ymaxthermo} 0. ${zmaxthermo}

create_box 2 simdom

create_atoms 1 region crystal

group crystal type 1
group cent_atoms region cent_atoms

Definition of thermostat

group thermal_atoms subtract crystal cent_atoms

define the mass of all atoms of Fe (*)

mass * 55.845

#Potential

pair_style eam/fs
pair_coeff * * PotentialBFe-H.fs Fe Fe
neighbor 2.0 bin
neigh_modify every 1

run_style verlet

CREATION OF PROJECTILE

Projectile is created outside the crystal (negative z)

create_atoms 1 single (v_xproj) (v_yproj) $(v_zproj)

region proj block (v_xproj-v_delta) (v_xproj+v_delta) (v_yproj-v_delta) (v_yproj+v_delta) (v_zproj-v_delta) (v_zproj+v_delta)
group proj region proj

INITIAL VORONOI TESSELATION

compute v1 crystal voronoi/atom occupation
compute r0 crystal reduce sum c_v1[1]
compute r1 crystal reduce sum c_v1[2]

thermo_style custom c_r0 c_r1
run 0

write_dump crystal custom coord_atoms_init.xyz id x y z

EQUILIBRIUM OF CRYSTAL

velocity crystal create {temp_thermo} {seed1} dist gaussian
compute msd1 crystal msd com yes

fix equilibrium crystal langevin {temp_thermo} {temp_thermo} 0.1 ${seed2}
fix 3 crystal nve

thermo_style custom step atoms time temp press lx ly lz c_msd1[4]

thermo 250

timestep 0.001

Equilibration

run 2000

unfix equilibrium
unfix 3
uncompute msd1
reset_timestep 0

Set the initial velocity of the projectile in A/ps

variable velproj equal 0.01sqrt(2${Ekproj}1.602e-196.022e23/55.845e-3)

variable velprojx equal {velproj}*{vx}
variable velprojy equal {velproj}*{vy}
variable velprojz equal {velproj}*{vz}

velocity proj set {velprojx} {velprojy} ${velprojz} units box

CASCADE

Thermostat

fix 2 all nve/limit 0.1
fix thermostat thermal_atoms langevin {temp_thermo} {temp_thermo} 0.1 ${seed3}

thermo_style custom step time atoms temp_r0 c_r1

thermo_modify lost ignore
thermo 100

fix adaptdt all dt/reset 1 NULL NULL 0.1 units box

run 2000

Store coordinates of atoms with Voronoi output to determine vacancies and SIAs

write_dump crystal custom coord_atoms_final.xyz id x y z c_v1[1] c_v1[2]

Dear LAMMPs users,

I am working on external irradiation of Fe. The projectile has an energy of
10 keV and a velocity vector along the z axis.

Since there can be atoms that are sputtered from the surface, I use f
boundary conditions on z (advice from Axel). Box size = 64x64*96

In order to calculate the number of vacancies and interstitials that are
formed during the cascade, I use the Voronoi package with voronoi/atom
occupation command. The output is then postprocessed with a Python script
that I wrote (I tested it in other cases and I believe it is ok).

I observed that depending on the number of processors, the number of
vacancies and interstitials can strongly vary.
For 24 procs I obtain around 200 V and 100 SIAs while for 12 processors, I
obtain around 1500 V and 800 SIAs. Very different.

are you sure, this is due to the voronoi compute giving different
results with a different number of processors, and not due to simply
getting different trajectories?
do you get the same processor dependence of your results, if you
record a trajectory without any analysis and then do the analysis via
rerun using different numbers of processors.

my suspicion is, that your kind of simulation is highly sensitive to
even small changes of how atoms are moving and thus you might need to
do multiple simulations artificially inducing such small changes and
your result is not just the numbers you get from a single individual
run/sample, but the distribution you get from a sufficiently large
number of runs.

axel.

  1. Changing the number of processors changes the order of operations.
  2. Floating point arithmetic is almost but not perfectly commutative.
  3. The chaotic nature of MD trajectories means that very small numerical differences eventually grow into large differences.

The logical conclusion is that you should expect results to vary “to some extent” when you change the number of processors. If this is happening, it will show up not just in the Voronoi analysis, but also in other properties such as temperature, energy, etc. One additional thing to consider is that Voronoi tesselations of points arranged on perfect lattices are often degenerate. This means that changing the order of operations can result in a different voronoi tesselation, before any atoms have been moved.

Aidan