Output from Voronoi varies with number of processors

LAMMPS LAMMPS Mailing List Mirror
1

  1. The defects identified by the Voronoi analysis are just the tip of the iceberg. You have a lot more information at your disposal that you should examine, such as energy per atom.

  2. That said, with the limited information that you have provided, I will conjecture that you are seeing an effect of superposition of disturbances from periodic images of the initial energy deposition. If the energy from this event radiates outward, it will meet its periodic image near the periodic boundaries,* half way between periodic images of the site where energy is initially deposited. This could give rise to additional defects.

  3. Other similar mechanisms might include dislocations that meet up half way between periodic images of their initial sources. The possibilities are endless.

  4. In order to guide yourself through this forest, you should first learn how to reproduce results in prior publications by experienced practitioners.

Aidan

  • The location in the z direction could be at the impact surface, or it could be on the back surface.
2

Hi Aidan,

Thanks for your reply,

  1. The defects identified by the Voronoi analysis are just the tip of the iceberg. You have a lot more information at your disposal that you should examine, such as energy per atom.

  2. That said, with the limited information that you have provided, I will conjecture that you are seeing an effect of superposition of disturbances from periodic images of the initial energy deposition. If the energy from this event radiates outward, it will meet its periodic image near the periodic boundaries,* half way between periodic images of the site where energy is initially deposited. This could give rise to additional defects.

I was thinking about something like that.
Therefore, I tried non-periodic conditions in x an y, but defects still appear near a boundary.
Any idea how to remove this artefact ?

I attach my input in case you have some time to have a look.

  1. Other similar mechanisms might include dislocations that meet up half way between periodic images of their initial sources. The possibilities are endless.

  2. In order to guide yourself through this forest, you should first learn how to reproduce results in prior publications by experienced practitioners.

In principle this case is not very difficult to reproduce. Actually, the results I obtain are consistent with what is obtained in the literature. The only problem is the production of artefact near the boundary. In papers, very few informations are given on the real setup of the simulation. For instance, when it comes to describe the thermostat, it is just said that a thermostat was applied to mimic heat conduction in a semi-infinite medium. Nothing else. Nothing is said about boundary conditions, or which type of thermostat, maximum distance atoms can move during a dt, etc…

Christophe

package gpu 2

dimension 3
units metal
atom_style atomic
atom_modify map hash

boundary p p p
lattice bcc 2.8553

thickness of the thermostat

variable thick_thermo equal 2
variable delta equal 0.2

variable Tfin equal “temp”

Temperature (K) of thermostat

variable temp_thermo equal 300.

Kinetic energy of the projectile in eV

variable Ekproj equal 10000.0

variable seed1 equal 54753
variable seed2 equal 87123
variable seed3 equal 12967
variable seed4 equal 97608

variable vx equal 0.0
variable vy equal 0.0
variable vz equal 1.0

Definition of simulation domain

region simdom block 0. 64 0. 64 -5 64
create_box 2 simdom

region crystal block 0. 64 0. 64 0. 64

create_atoms 1 region crystal
group crystal region crystal

region cent_atoms block 2 62 2 62 0 62
group cent_atoms region cent_atoms

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

create_atoms 1 single 32 32 -4
region proj block 31.80 32.20 31.80 32.20 -4.20 -3.80

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 300.0 54753 dist gaussian loop local

fix equilibrium crystal langevin 300 300 0.1 87123
fix 3 crystal nve

compute msd1 crystal msd com yes

compute tempthermostat thermal_atoms temp
compute tempcore cent_atoms temp

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

thermo 250

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

run 1000

unfix equilibrium
unfix 3
uncompute msd1

Set the initial velocity of the projectile

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: MOLECULAR DYNAMICS

fix 2 all nve/limit 0.1

fix thermostat thermal_atoms langevin 300.0 300 0.1 12967

thermo_style custom step time atoms c_tempcore c_tempthermostat c_kinecproj c_peproj c_r0 c_r1

thermo_modify lost ignore
thermo 100

run 5000

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

3

If you read the right papers (the literature on MD simulation of ion and radiation bombardment is large) you will find more details about best practices for boundary conditions. The basic idea is to remove energy that reaches the boundary. This can never be done perfectly. One possibility is to apply a strong Langevin thermostat to the layer of atoms near the boundary.

4

If you read the right papers (the literature on MD simulation of ion and radiation bombardment is large) you will find more details about best practices for boundary conditions. The basic idea is to remove energy that reaches the boundary. This can never be done perfectly. One possibility is to apply a strong Langevin thermostat to the layer of atoms near the boundary.

That is what I have in my input. I tried two different thermostats:

  1. 2-3 layers at 300 K with strong dissipative Langevin thermostat (except surface that is bombarded). Damping factor 0.03 (metal units). Is it stron enough ?

  2. A more sophisticated thermostat recommended by Axel. On top of the thermostat above, an outer 2-3 layers of immobile atoms (0K) to squash phonons and avoid surface reconstruction.

None of these thermostats helped solving the problem.

Christophe

5

Just because you failed to get it to work does not mean it is the wrong approach. I have no idea whether 0.03 is strong enough. There are many factors to consider and many tradeoffs to be made. If you are serious about learning how to do quality simulations, you will take my original advice and read the best papers in the field, reproducing their work, sharpening your skills and understanding, before embarking on your own research. You won’t get that kind of in-depth knowledge by repeatedly asking questions on this mailing list.

6

Hi Aidan and Axel,

I came back to the basics and simulated an external irradiation in a box simply equilibrated at 300K and all at NVE. No thermostat. Well, the problem was still there.

After some headache I fount out that the problem comes from equilibration stage, namely, the Langevin thermostat. To equilibrate I was using (in pseudo LAMMPS language):

fix all langevin 300K
fix all nve

If instead, I initially use a Berendsen thermostat to equilibrate the crystal at 300K, then there is no problem during the cascade. The problem of defects that form on the boundary disappears.

Then, I reintroduced a thermostat (to mimic heat conduction in semi-infinite bulk) and all was ok. No problem of defects on the thermostat neither with the number of processors.

Maybe Axel has an idea why Langevin would cause such problems ?

Christophe

7

Maybe Axel has an idea why Langevin would cause such problems ?

From Aidan’s earlier email:

  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.

Fix langevin will do different (random) things on different #'s of processors.

So it’s the same issue.

Steve

8

There is nothing wrong with Langevin, when used correctly. In general, it is more physically correct than Berendsen, which tends to over-damp kinetic energy fluctuations. The real problem is that your simulation protocol appears to be very sensitive to minor changes in the way the calculations are performed. This type of fragility is usually indicative of more serious underlying problems, such as sensitive dependence on initial conditions.

9

There is nothing wrong with Langevin, when used correctly. In general, it is more physically correct than Berendsen, which tends to over-damp kinetic energy fluctuations.

That is what a colleague of mine with long experience with MD told me. Therefore, I was using Langevin to equilibrate the Fe crystal.

The real problem is that your simulation protocol appears to be very sensitive to minor changes in the way the calculations are performed. This type of fragility is usually indicative of more serious underlying problems, such as sensitive dependence on initial conditions.

What I want to simulate is pretty simple:

  • cubic Fe crystal equilibrated at 300K

  • a projectile (outside the crystal) is given an energy

  • run (with or without thermostat)

There is nothing special. I can understand that, due to the stochasticity of Nature, results (number of defects and clusters) can vary from one run to another. But what results difficult to understand to me is how the type of thermostat used at the beginning to equilibrate the piece of crystal can influence the result in such way that defects form where they should not (on the boundary far from initial impact). For me, clearly this is an artefact.

Or perhaps, if there is nothing wrong with Langevin, it is the way I equilibrate the crystal that is wrong. This is how I do, maybe you could tell me if you see something suspicious:

Definition of crystal region

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

Definition of simulation domain. It is larger than the crystal region to leave

some space for the external projectile (z = -3)

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

Create simulation domain.

create_box 2 simdom

Create atoms in the crystal

create_atoms 1 region crystal

group crystal region crystal

Equilibration (projectile not yet created)

velocity all create 300. {seed1} dist gaussian mom yes loop local fix thermocrystal crystal langevin 300. 300. 1.0 {seed4}

fix 3 crystal nve

timestep 0.01

run 1000

Do you see anything wrong in the equilibration ?

Christophe

10

If you phone up a car mechanic and tell him your car is making a funny noise, he is not going to be able to help you, unless his name is Click or Clack. This is a similar situation. As I suggested before, find a good paper in the vast literature (I have a paper on MD simulation of radiation damage from 1960 by Vineyard et al.), one that you know you should be able to reproduce, and reproduce it. Do not rule out the possibility that your apparent troubles are merely an artifact of the Voronoi analysis, which you seem to be using as your sole diagnostic tool. Also, the best way to find errors in your input is to find errors in your output. Don’t look just at the Voronoi analysis, but other more fundamental properties, such as per-atom energy. Also, don’t look just at the final result, but how things change over time. You can learn a lot by looking “under the hood” of your car/simulation.

Also, make your timestep smaller.

Aidan

11

If you phone up a car mechanic and tell him your car is making a funny noise, he is not going to be able to > help you, unless his name is Click or Clack.

Or possibly Donald. He would definitely tell you what is wrong.

Steve

12

Hi all,

If you phone up a car mechanic and tell him your car is making a funny noise, he is not going to be able to help you, unless his name is Click or Clack. This is a similar situation. As I suggested before, find a good paper in the vast literature (I have a paper on MD simulation of radiation damage from 1960 by Vineyard et al.), one that you know you should be able to reproduce, and reproduce it. Do not rule out the possibility that your apparent troubles are merely an artifact of the Voronoi analysis, which you seem to be using as your sole diagnostic tool. Also, the best way to find errors in your input is to find errors in your output. Don’t look just at the Voronoi analysis, but other more fundamental properties, such as per-atom energy. Also, don’t look just at the final result, but how things change over time. You can learn a lot by looking “under the hood” of your car/simulation.

Also, make your timestep smaller.

I could not talk to Donald but Rajoy gave me some hints…:smiley:

I understood that the equilibration part was not well done. It was too short. I was only running for 1 or 2 ps (at 300K). In the case where one of the surfaces is free (which is the case), if the equilibration run is too short, the surface continues expanding during production. Consequently, when you use Voronoi to determine where defects are after the cascade, which compares initial and final positions of atoms, you get that there are many defects in the bulk. Actually, Voronoi works well. It tells you that many atoms have moved form their initial position. Not its problem if you haven’t done properly equilibration ! They have not moved because of the cascade, but because of expansion that goes on during production run.
So it was not a problem of Langevin, but equilibration that was too short and timestep also too large.
To solve the problem now I work at 5K instead of 300K (expansion is limited) and do longer run for equilibration.

Thank you all for your help.

Christophe