radiation damage not annealing but exploding

LAMMPS LAMMPS Mailing List Mirror
1

LAMMPS users, I am having trouble with simulating radiation damage. I am modeling quartz. I equilibrated the system at 300K in a previous run. When I run simulation, the temperature slowly fluctuates, then at some point, it randomly spikes and the simulation crashes.

I am using fix nve because I do not want to impose an artificial drag by implementing a thermostat. Is the reason why I am seeing temperature spike after a few iterations due to not defining a region around the simulation box with a fix nvt? I know this is used to absorb the energy so the cascade doesn’t go through periodic boundaries and interact with itself. However, I cannot imagine how this would have affect the temperature in the middle. I even tried fix nvt and got similar results.

Does anyone see an outright flaw in my simulation method? The timestep is not too large as I used variable time step and when the simulation crashes, it is not on the minimum time step possible.

Also, the region for PKA only contains one atom.

units metal
atom_style full
kspace_style pppm 1.0e-6

read_restart data.restart

pair_style buck/coul/long 15
pair_coeff 1 1 1388.7730 0.36231884058 175.0000
pair_coeff 1 2 18003.7572 0.205204814926 133.5381
pair_coeff 2 2 0 .1 0

region rPKA sphere 40 45 50 1.1 units box
group PKA region rPKA

thermo 1

timestep 0.001

fix 1 all dt/reset 1 0.00001 0.01 0.01 units box

fix 2 all nve

velocity PKA set 548.0 0 0 units box

run 30000

2

Did you try using a lower velocity for the PKA to see if the results changed? (I calculated that the velocity that you assigned to PKA was 250keV.) I also notice that your script doesn’t have a dump command. It would probably by useful to actually visualize the simulation prior to its failure. -Vanessa

3

Yes, I tried changing the velocities and the temperature acted very similarly. The script I provided is not the whole script, just the necessary things that affected the physics. I have several compute and dump commands that I did not include just because I did not think it would be useful in troubleshooting what my error might be.

Ben

4

Hello Ben, can you provide a minimal runnable input (i.e. add a
create_atoms command rather than a read data) that allows us to
reproduce the error?
Daniel

5

Perhaps I may have fixed it but I do not completely understand the physics. I was just using the microcanonical ensemble for all the system earlier when the temperature was exploding. I added a few layers on the outer sides of the simulation box and applied NVT ensemble and it seems to be working fine thus far. I replicated the data file below by 10 10 20 any equilibrated at 300K for 30ps using NVT ensemble.

The data file is:

LAMMPS data file. CGCMM style. atom_style full generated by VMD/TopoTools v1.2 on Thu Sep 27 17:25:35 EDT 2012
36 atoms

2 atom types
0.000 8.509 xlo xhi
0.000 9.826 ylo yhi
0.000 5.405 zlo zhi

Masses

2 28.085501 # Si
1 15.999400 # O

Atoms

1 1 1 -1.2 1.760 0.299 1.157 # O
2 1 1 -1.2 3.116 8.744 2.959 # O
3 1 1 -1.2 3.634 0.783 4.761 # O
4 1 1 -1.2 1.139 1.375 4.248 # O
5 1 1 -1.2 0.621 3.240 2.446 # O
6 1 1 -1.2 2.495 2.755 0.644 # O
7 1 2 2.4 2.000 8.671 1.802 # Si
8 1 2 2.4 0.000 2.310 3.603 # Si
9 1 2 2.4 2.255 1.302 0.000 # Si
10 1 1 -1.2 6.015 7.668 1.157 # O
11 1 1 -1.2 7.371 6.288 2.959 # O
12 1 1 -1.2 7.888 8.153 4.761 # O
13 1 1 -1.2 5.393 8.744 4.248 # O
14 1 1 -1.2 4.876 0.783 2.446 # O
15 1 1 -1.2 6.750 0.299 0.644 # O
16 1 2 2.4 6.255 6.215 1.802 # Si
17 1 2 2.4 4.255 9.679 3.603 # Si
18 1 2 2.4 6.509 8.671 0.000 # Si
19 1 1 -1.2 1.760 5.212 1.157 # O
20 1 1 -1.2 3.116 3.831 2.959 # O
21 1 1 -1.2 3.634 5.696 4.761 # O
22 1 1 -1.2 1.139 6.288 4.248 # O
23 1 1 -1.2 0.621 8.153 2.446 # O
24 1 1 -1.2 2.495 7.668 0.644 # O
25 1 2 2.4 2.000 3.758 1.802 # Si
26 1 2 2.4 0.000 7.223 3.603 # Si
27 1 2 2.4 2.255 6.215 0.000 # Si
28 1 1 -1.2 6.015 2.755 1.157 # O
29 1 1 -1.2 7.371 1.375 2.959 # O
30 1 1 -1.2 7.888 3.240 4.761 # O
31 1 1 -1.2 5.393 3.831 4.248 # O
32 1 1 -1.2 4.876 5.696 2.446 # O
33 1 1 -1.2 6.750 5.212 0.644 # O
34 1 2 2.4 6.255 1.302 1.802 # Si
35 1 2 2.4 4.255 4.766 3.603 # Si
36 1 2 2.4 6.509 3.758 0.000 # Si

6

Perhaps I may have fixed it but I do not completely understand the physics. I was just using the microcanonical ensemble for all the system earlier when the temperature was exploding. I added a few layers on the outer sides of the simulation box and applied NVT ensemble and it seems to be working fine thus far. I replicated the data file below by 10 10 20.

The data file is:

LAMMPS data file. CGCMM style. atom_style full generated by VMD/TopoTools v1.2 on Thu Sep 27 17:25:35 EDT 2012
36 atoms

2 atom types
0.000 8.509 xlo xhi
0.000 9.826 ylo yhi
0.000 5.405 zlo zhi

Masses

2 28.085501 # Si
1 15.999400 # O

Atoms

1 1 1 -1.2 1.760 0.299 1.157 # O
2 1 1 -1.2 3.116 8.744 2.959 # O
3 1 1 -1.2 3.634 0.783 4.761 # O
4 1 1 -1.2 1.139 1.375 4.248 # O
5 1 1 -1.2 0.621 3.240 2.446 # O
6 1 1 -1.2 2.495 2.755 0.644 # O
7 1 2 2.4 2.000 8.671 1.802 # Si
8 1 2 2.4 0.000 2.310 3.603 # Si
9 1 2 2.4 2.255 1.302 0.000 # Si
10 1 1 -1.2 6.015 7.668 1.157 # O
11 1 1 -1.2 7.371 6.288 2.959 # O
12 1 1 -1.2 7.888 8.153 4.761 # O
13 1 1 -1.2 5.393 8.744 4.248 # O
14 1 1 -1.2 4.876 0.783 2.446 # O
15 1 1 -1.2 6.750 0.299 0.644 # O
16 1 2 2.4 6.255 6.215 1.802 # Si
17 1 2 2.4 4.255 9.679 3.603 # Si
18 1 2 2.4 6.509 8.671 0.000 # Si
19 1 1 -1.2 1.760 5.212 1.157 # O
20 1 1 -1.2 3.116 3.831 2.959 # O
21 1 1 -1.2 3.634 5.696 4.761 # O
22 1 1 -1.2 1.139 6.288 4.248 # O
23 1 1 -1.2 0.621 8.153 2.446 # O
24 1 1 -1.2 2.495 7.668 0.644 # O
25 1 2 2.4 2.000 3.758 1.802 # Si
26 1 2 2.4 0.000 7.223 3.603 # Si
27 1 2 2.4 2.255 6.215 0.000 # Si
28 1 1 -1.2 6.015 2.755 1.157 # O
29 1 1 -1.2 7.371 1.375 2.959 # O
30 1 1 -1.2 7.888 3.240 4.761 # O
31 1 1 -1.2 5.393 3.831 4.248 # O
32 1 1 -1.2 4.876 5.696 2.446 # O
33 1 1 -1.2 6.750 5.212 0.644 # O
34 1 2 2.4 6.255 1.302 1.802 # Si
35 1 2 2.4 4.255 4.766 3.603 # Si
36 1 2 2.4 6.509 3.758 0.000 # Si

7

Actually, it did not fix it. After letting the simulation run longer, I am having the same issue.

Ben

8

Another issue, if I recall, is that a Buckingham potential
has a known issue for high temperature (or high velocity
in a cascade model) that 2 atoms can get too close

together and “explode”. For this reason, you often
need to spline in a more repulsive core when using
Buckingham for high T. E.g. Paul (CCd) has done this
for silica and has the spline fit tabulated potential.

Steve

9

Yes, I sent that spline fit tabulated potential for SiO2 to Ben quite some time ago and he has already experimented with it quite a bit. But it is probably not adequate for production PKA simulations like I think Ben wants since it includes a fairly arbitrary (made up) repulsive core.

I agree that with just plain Buck like he’s using below, he can get “explosions”, especially since he’s got a high velocity PKA atom in his simulation. So probably neither force field would be adequate for production PKA simulations.

Paul

10

I wonder if using a DPD like thermostat may help a bit with achieving a certain system stabilization, at least for runs with PKA energies that are not too high but currently make the simulations crash. Not sure what has been done in the lit related to radiation damage + DPD, yet, DPD is a thermostat highly employed in the simulation of liquid systems in non-equilibrium conditions due to its momentum conservation property. I’m aware the repulsive core is the main issue here but perhaps a DPD-like band-aid could help stabilizing some of the runs.
Carlos

11

Right, would that be an issue with PKA of only 15eV? I change the energy and I get the same effect. It seems like whatever energy I use, I get the same effect. Perhaps the incoming velocity is just too high even in the low eV range.

Ben

12

So do you get the same response with a PKA of 0 eV? If so, it shows there is some other problem with your setup

Nigel