Problem with computing defects in irradiated crystal with expansion

Dear all,

I am working on the irradiation of a Fe crystal (100aox100aox100ao) with an open surface. The projectile is a Fe atom coming from outside the crystal with E = 10keV.
There is a thermostat of 2 layers surrounding the box, except for the open surface.
The equilibration stage is made as follows (following the advice of Axel):

  • Initialization of energies of atoms at T
  • Minimization of potential energy
  • Equilibration during 2.5 ps with Langevin with short time constant (0.3 ps)
  • Equilibrationduring 2.5 ps with NVT with longer time constant (2.0 ps)

In order to calculate the defects (interstitials and vacancies), I use the OVITO code with the Wigner-Seitz analyse. Ir works well in general.
However, when there is an open surface, in some cases, I observed that the number of defects is very high, much larger than the average number of defects expected for this energy. After digging a while, I found out that it is because the atoms of the front and of the back surfaces are interpreted as defects by the WS analysis. My guess is that there is a slight expansion of the material during the cascade. Then, when the WS tries to determine which atoms are defects, it sees that all atoms from the surface have moved and then decides they are all defects. As a result, there are many spurious defects in the output (see attachment where all surface atoms are considered as defects by the WS).

Any idea how to solve that?
Is there a way to minimise expansion? Or should I go for a method to “clean” the output before using the WS analysis?

Many thanks in advance and best regards,
Christophe