Error in high energy irradiation simulations with electronic stopping

acgc99 · April 4, 2025, 9:42am

I need to simulate a 4 MeV oxygen being irradiated into a Fe thin sample. This is an example with reduced material width:

units metal
atom_style atomic
boundary p p f
atom_modify map yes
lattice bcc 2.8664999

region r_full block 0 30 0 30 -5 35
create_box 2 r_full

region r_injection block 0 30 0 30 30 35

region r_lamella block 0 30 0 30 0 30
region r_inner_lamella block 2 28 2 28 2 28
create_atoms 1 region r_lamella

group g_lamella region r_lamella
group g_inner_lamella region r_inner_lamella
group g_thermoboundary subtract g_lamella g_inner_lamella

fix integrate all nve
fix reset all dt/reset 10 1e-05 0.001 0.028664999 units box

pair_style hybrid eam/alloy zbl 11.4659996 17.198999399999998  # 4a0 6a0
pair_coeff * * eam/alloy mendelev.eam.alloy Fe NULL
pair_coeff 1 2 zbl 26.0 8.0
pair_coeff 2 2 zbl 8.0 8.0

mass 1 55.85
mass 2 16.00

fix esloss all electron/stopping 10.0 FeO.txt

compute ke all ke/atom

thermo_style custom step dt time temp etotal ke f_esloss
thermo 10
thermo_modify lost ignore

compute temp_lamella g_lamella temp
thermo_modify temp temp_lamella

dump out g_lamella custom 100 sim.xyz id element x y z c_ke
dump_modify out element 26 8

velocity all create 200.0 1 loop geom
run 1000

create_atoms 2 single 15 15 33 units lattice
group injected region r_injection

#velocity injected set -100.0 -100.0 -37181.0 units box # Fe 588 = 1kev
velocity injected set -100.0 -100.0 -69457.0 units box # O 4 MeV

fix control g_thermoboundary langevin 100 100 1.0 1
run 1000

And this is the result:

>>> mpiexec -np 10 lmp -in .\surface.in
LAMMPS (19 Nov 2024)
  using 2 OpenMP thread(s) per MPI task
Loaded 1 plugins from C:\Users\acgc99\AppData\Local\LAMMPS 64-bit 19Nov2024-MSMPI\plugins
Lattice spacing in x,y,z = 2.8664999 2.8664999 2.8664999
Created orthogonal box = (0 0 -14.3325) to (85.994997 85.994997 100.3275)
  1 by 2 by 5 MPI processor grid
Created 54900 atoms
  using lattice units in orthogonal box = (0 0 -14.3325) to (85.994997 85.994997 100.3275)
  create_atoms CPU = 0.002 seconds
54900 atoms in group g_lamella
37259 atoms in group g_inner_lamella
17641 atoms in group g_thermoboundary
WARNING: Temperature for thermo pressure is not for group all (src/thermo.cpp:533)
Neighbor list info ...
  update: every = 1 steps, delay = 0 steps, check = yes
  max neighbors/atom: 2000, page size: 100000
  master list distance cutoff = 19.198999
  ghost atom cutoff = 19.198999
  binsize = 9.5994997, bins = 9 9 12
  4 neighbor lists, perpetual/occasional/extra = 3 1 0
  (1) pair eam/alloy, perpetual, skip trim from (4)
      attributes: half, newton on, cut 7.30067539
      pair build: skip/trim
      stencil: none
      bin: none
  (2) pair zbl, perpetual, skip from (4)
      attributes: half, newton on
      pair build: skip
      stencil: none
      bin: none
  (3) fix electron/stopping, occasional
      attributes: full, newton on
      pair build: full/bin/atomonly
      stencil: full/bin/3d
      bin: standard
  (4) neighbor class addition, perpetual
      attributes: half, newton on
      pair build: half/bin/atomonly/newton
      stencil: half/bin/3d
      bin: standard
Setting up Verlet run ...
  Unit style    : metal
  Current step  : 0
  Time step     : 0.001
Per MPI rank memory allocation (min/avg/max) = 13.21 | 17.77 | 22.64 Mbytes
   Step           Dt            Time           Temp          TotEng         KinEng        f_esloss
         0   0.001          0              200           -216896.09      1419.2505      0
       100   0.001          0.1            103.08995     -216895.76      731.55236      0
       200   0.001          0.2            90.695147     -216895.73      643.59568      0
       300   0.001          0.3            102.95825     -216895.8       730.61775      0
       400   0.001          0.4            102.36501     -216895.79      726.40797      0
       500   0.001          0.5            98.791262     -216895.77      701.04776      0
       600   0.001          0.6            102.18754     -216895.79      725.14857      0
       700   0.001          0.7            103.03707     -216895.79      731.17705      0
       800   0.001          0.8            100.30207     -216895.77      711.76881      0
       900   0.001          0.9            100.72539     -216895.78      714.77283      0
      1000   0.001          1              101.82163     -216895.78      722.552        0
Loop time of 11.8746 on 20 procs for 1000 steps with 54900 atoms

Performance: 7.276 ns/day, 3.298 hours/ns, 84.214 timesteps/s, 4.623 Matom-step/s
90.9% CPU use with 10 MPI tasks x 2 OpenMP threads

MPI task timing breakdown:
Section |  min time  |  avg time  |  max time  |%varavg| %total
---------------------------------------------------------------
Pair    | 6.1165     | 7.7204     | 8.8674     |  44.2 | 65.02
Neigh   | 0          | 0          | 0          |   0.0 |  0.00
Comm    | 2.4274     | 3.3227     | 4.6453     |  51.0 | 27.98
Output  | 0.12483    | 0.15165    | 0.1727     |   3.8 |  1.28
Modify  | 0.15736    | 0.21472    | 0.26699    |   9.6 |  1.81
Other   |            | 0.4651     |            |       |  3.92

Nlocal:           5490 ave        7200 max        2700 min
Histogram: 2 2 0 0 0 0 0 0 0 6
Nghost:        32957.2 ave       46156 max       21650 min
Histogram: 4 0 0 0 0 0 4 0 0 2
Neighs:              0 ave           0 max           0 min
Histogram: 10 0 0 0 0 0 0 0 0 0
FullNghs:       723240 ave      979200 max      308700 min
Histogram: 4 0 0 0 0 0 0 0 0 6

Total # of neighbors = 7232400
Ave neighs/atom = 131.7377
Neighb
job aborted:
[ranks] message

[0-8] terminated

[9] application aborted
aborting MPI_COMM_WORLD (comm=0x44000000), error 1, comm rank 9

---- error analysis -----

[9] on DESKTOP-0PR1AFH
lmp aborted the job. abort code 1

---- error analysis -----

The error message is not very explicative. I’ve found two things that can solve the issue (separately) but modify the simulation:

Remove electronic stopping
Reducing oxygen velocity by a factor of 100.

This made me suspect about the upper limit in the energy of the electronic stoppings files, but I increased it up to 10 GeV and the issue persists.

Here you can find the electronic stoppings (from SRIM) and interatomic potential:
FeO.txt (5.9 KB)
fe_mendelev_eam_alloy.potential (395.7 KB)

Germain · April 4, 2025, 6:56pm

LAMMPS output is buffered in memory by default and only written when reaching a given size to avoid too much reading/writing and saving some time and resources. You can turn the flush option in the thermo_modify command. On some systems however, the cluster can give you hand back before the error output is written. This is apparently your case. However, your simulation crashes early enough to be run on a single core.

When running locally you can directly get the error output in your terminal. In your case, you have a neighbor list overflow. Please see the new guidelines on how to debug your issue yourself. You’ll also find more information in the neigh_modify page.

acgc99 · April 5, 2025, 9:38am

Thank you!

It turns out the ZBL cutoff I was using is too big. It’s hard to find a reference where they explicitely say that value. I’ve only found one and it says for distances < 0.1 nm. Applying that the issue is resolver.

I would like to point out that when I don’t use parallelization I get:

ERROR on proc 0: Neighbor list overflow, boost neigh_modify one (src/npair_bin.cpp:248)

It would be nice to know the value that’s causing the overflow, to understand the issue and apply neigh_modify if necessary.

Germain · April 5, 2025, 10:28am

This issue should appear whether you use parallelization or not if your neighbor list settings are the same. It depends on the neighbor list of individual atoms, which is affected by the cutoff distances of the pair styles. This is why you solved your issue by reducing the cutoff distance.

The way neighbor lists are built is detailed in the neighbor and neigh_modify pages. The need of either modifying your model or the neighbor list parameters is up to the use case and cannot be determined in advance. The default are fine “for most simulations”.

From experience, most of the time, if some issue appears with the neighbor lists memory use, it is a model problem more often than an “unusual simulation”.