Reax/c and compute rdf command in vitreous SiO2 calculation

LAMMPS LAMMPS Mailing List Mirror
1

Dear all

version:lammps-28Jun14

Recently, I calculated 12960 atoms of vitreous SiO2 with reax/c potential reported in the following paper.

<http://scitation.aip.org/content/aip/journal/jcp/132/17/10.1063/1.3407433>

The procedure was Energy Minimization -> NVT(4000K-300K) -> NVT(4000K-300K, the final output was ‘npt.out’).

However, the coordination number of Si-O became 4.5 when I calculated rdf by the following input files.

units real

atom_style full

boundary p p p

read_data …/npt.out

timestep 0.5 # (fs)

pair_style reax/c NULL

pair_coeff * * …/SiOH.ff Si O

neighbor 0.5 bin

neigh_modify every 10 check yes

thermo_style custom step temp etotal pe ke press density

fix equil all npt temp 300.0 300.0 100.0 iso 1.0 1.0 500.0

fix equilchg all qeq/reax 1 0.0 10.0 1e-6 reax/c

compute Si_O all rdf 1000 1 2

fix Si_O all ave/time 1 1 1 c_Si_O file Si_O.rdf mode vector

compute Si_Si all rdf 1000 1 1

fix Si_Si all ave/time 1 1 1 c_Si_Si file Si_Si.rdf mode vector

compute O_O all rdf 1000 2 2

fix O_O all ave/time 1 1 1 c_O_O file O_O.rdf mode vector

thermo 1

run 0 (Note that the result was the almost same if I set run more than 10000.)

Surprisingly, the coordination number of Si-O became 4.0 when I used tersoff potential with the SAME output ‘npt.out’.

units metal

atom_style full

boundary p p p

read_data …/npt.out

timestep 0.0005 # (ps)

pair_style tersoff

pair_coeff * * …/SiO2.tersoff Si O

neighbor 0.5 bin

neigh_modify every 10 check yes

fix equil all npt temp 300.0 300.0 0.1 iso 1.0 1.0 1.0

thermo_style custom step temp etotal pe ke press density

compute Si_O all rdf 200 1 2

fix Si_O all ave/time 1 1 1 c_Si_O file Si_O.rdf mode vector

compute Si_Si all rdf 200 1 1

fix Si_Si all ave/time 1 1 1 c_Si_Si file Si_Si.rdf mode vector

compute O_O all rdf 200 2 2

fix O_O all ave/time 1 1 1 c_O_O file O_O.rdf mode vector

thermo 1

run 0 (Note that the result was the almost same if I set run more than 10000.)

I found that the similar bug had already been reported.

<https://sourceforge.net/p/lammps/mailman/message/30849776/>

Did anyone find the reasons?

Best regards

2

Hi Nakano-san,

Why is that there are two NVT stages with the same temperature ramping? Usually one would increase the temperature from room temperature to a high temperature, say 4000K, then hold at 4000K for another duration of time, and then decrease the temperature down to room temperature again. Before the RDF compute run (which is usually at NVT), there is usually another NPT stage to equilibrate the pressure. Currently your procedure combines the RDF together with NPT.

For your ReaxFF simulation, the time step size might be a bit too large - usually 0.2 fs is a good number. Do you get the correct density compared to experimental value for the RDF run?

Best regards,
Ray

3

Dear Ray Shan

Thank you for your prompt reply!

Why is that there are two NVT stages with the same temperature ramping?

I’m sorry. That was my mistype…

The correct procedure was Energy Minimization -> NVT(4000K-300K) -> NPT(4000K-300K).

The final output was ‘npt.out’.

usually 0.2 fs is a good number.

Thank you for your advice! I changed the timestep and retried the procedure.

Do you get the correct density compared to experimental value for the RDF run?

Yes. The density is 2.12 during the rdf calculation (at NVT), which is almost consistent with

the experimental value (2.20) and the calculated value (2.14) for SiO2 glass.

http://scitation.aip.org/content/aip/journal/jcp/132/17/10.1063/1.3407433

Input for rdf calculation

units real

atom_style full

boundary p p p

read_data …/npt.out

timestep 0.2 # (fs)

pair_style reax/c NULL

pair_coeff * * …/SiOH.ff Si O

neighbor 0.5 bin

neigh_modify every 10 check yes

thermo_style custom step temp etotal pe ke press density

fix equil all nvt temp 300.0 300.0 100.0

fix equilchg all qeq/reax 1 0.0 10.0 1e-6 reax/c

compute Si_O all rdf 1000 1 2

fix Si_O all ave/time 1000 1 1000 c_Si_O file Si_O.rdf mode vector

compute Si_Si all rdf 1000 1 1

fix Si_Si all ave/time 1000 1 1000 c_Si_Si file Si_Si.rdf mode vector

compute O_O all rdf 1000 2 2

fix O_O all ave/time 1000 1 1000 c_O_O file O_O.rdf mode vector

thermo 1000

run 10000

Output for rdf calculation

Step Temp TotEng PotEng KinEng Press Density

0 302.978 -882747.94 -888165.76 5417.8162 -480.49822 2.1248183

1000 300.88928 -882924.8 -888305.27 5380.4659 1921.8454 2.1248183

2000 299.17253 -882877.27 -888227.04 5349.7671 -1308.8525 2.1248183

3000 304.76592 -882774.3 -888224.08 5449.7874 1356.3489 2.1248183

…

The coordination number of Si, however, is still 4.5e$B!De(B

Time-averaged data for fix Si_O

TimeStep Number-of-rows

Row c_Si_O[1] c_Si_O[2] c_Si_O[3]

0 1000

1 0.005 0 0

…

147 1.465 0 0

148 1.475 0.0429371 0.0005

149 1.485 0.211804 0.003

150 1.495 1.21208 0.0175

151 1.505 3.58809 0.061

152 1.515 7.89575 0.158

153 1.525 13.8177 0.33

154 1.535 18.7129 0.566

155 1.545 25.0461 0.886

156 1.555 30.1722 1.2765

157 1.565 34.4792 1.7285

158 1.575 35.1348 2.195

159 1.585 33.0196 2.639

160 1.595 28.0903 3.0215

161 1.605 23.5711 3.3465

162 1.615 21.2386 3.643

163 1.625 15.6716 3.8645

164 1.635 12.3005 4.0405

165 1.645 7.59465 4.1505

166 1.655 5.52505 4.2315

167 1.665 3.90883 4.2895

168 1.675 2.56376 4.328

169 1.685 1.18446 4.346

170 1.695 0.877893 4.3595

…

Best regards

Kousuke

4

In the paper you referenced, they used a time step of 0.5 fs. That is concerning by at least a factor of two, as Ray has pointed out.

Also, they defined the coordination number as “Coordination is calculated as the average number with the first coordination shell.” And calculated a value of 3.95 for Si-Si. How does your O-O coordination number compare? Do you get their value of 1.99?

How does the method you used to calculate coordination compare to theirs? What was the effect of lowering the time step? Did you EXACTLY follow their procedure in performing your calculation?

Jim

5

Also, your version of LAMMPS is 2 years old. As Axel usually comments, try the latest version and see what you get.

Jim

6

Even more also, try the Kokkos version of ReaxFF.

Jim

7

Also, you should verify that you are able to accurately reproduce simple structural and energetic results from the paper. Without positive validation, you can not be confident that you are using the same potential.

Aidan

8

Dear James Kress

Thank you for your reply!

Also, your version of LAMMPS is 2 years old. As Axel usually comments, try the latest version and see what you get.

I compiled the latest stable lammps and retried calculations.

I got better density 2.14, which is the same as the reported value (2.14).

However, coordination numbers are still strange…

The following results were calculated by the latest lammps.

Also, they defined the coordination number as “Coordination is calculated as the average number with the first coordination shell.”

And calculated a value of 3.95 for Si-Si. How does your O-O coordination number compare? Do you get their value of 1.99?

No. I got ~2.24 for O-Si coordination number.

The output by RDF calculation

Step Temp TotEng PotEng KinEng Press Density

0 307.41008 -880555.27 -886052.34 5497.07 -422.25745 2.1367502

1000 299.94844 -880775.85 -886139.49 5363.6418 401.27741 2.1367502

2000 302.97604 -880851.23 -886269.01 5417.781 773.02615 2.1367502

….

The results of RDF are

Si-O

Time-averaged data for fix Si_O

TimeStep Number-of-rows

Row c_Si_O[1] c_Si_O[2] c_Si_O[3]

0 1000

1 0.005 0 0

2 0.015 0 0

…

229 2.285 0 4.478

230 2.295 0 4.478

O-Si

Time-averaged data for fix O_Si

TimeStep Number-of-rows

Row c_O_Si[1] c_O_Si[2] c_O_Si[3]

0 1000

1 0.005 0 0

2 0.015 0 0

…

229 2.285 0 2.239

230 2.295 0 2.239

How does the method you used to calculate coordination compare to theirs? What was the effect of lowering the time step?

Did you EXACTLY follow their procedure in performing your calculation?

I believe that I exactly followed their procedures. In detail,

NVT

timestep 0.5 # (fs)

velocity all create 4000 1 mom yes rot yes

Annealing

fix anneal all nvt temp 4000.0 300.0 100.0

fix chgeanneal all qeq/reax 1 0.0 10.0 1e-6 reax/c

thermo 1000

run 296000

NPT

timestep 0.5 # (fs)

velocity all create 4000 1 mom yes rot yes

melting

fix melting all npt temp 4000.0 4000.0 100.0 iso 1.0 1.0 500.0

fix chgemelt all qeq/reax 1 0.0 10.0 1e-6 reax/c

thermo 1000

run 150000

Annealing

fix anneal all npt temp 4000.0 300.0 100.0 iso 1.0 1.0 500.0

fix chgeanneal all qeq/reax 1 0.0 10.0 1e-6 reax/c

thermo 1000

run 296000

I confirmed that coordination numbers were not improved when I changed a time step to 0.2 fs …

Even more also, try the Kokkos version of ReaxFF.

Thank you. I will try it soon.

Best regards

Kousuke

9

Dear Kousuke,

Thanks for your reply. The only question you did not answer was:

“Coordination is calculated as the average number with the first coordination shell.”

Is this the exact same method you used to calculate coordination number?

Also I asked about O-O coordination number, not O-Si.

Jim

10

Dear Kress

I am sorry for my late reply.

“Coordination is calculated as the average number with the first coordination shell.”

I am not sure… I am checking now.

Also I asked about O-O coordination number, not O-Si.

I am sorry. O-O coordination number seems to be around 7.5, but it does not show plateau different from Si-O.

(by LAMMPS)

Time-averaged data for fix O_O

TimeStep Number-of-rows

Row c_O_O[1] c_O_O[2] c_O_O[3]

1000 200

1 0.025 0 0

…

49 2.425 4.41075 1.3897

50 2.475 5.6495 2.3175

51 2.525 6.17976 3.3738

52 2.575 5.90208 4.423

53 2.625 5.13924 5.3724

54 2.675 4.02311 6.1442

55 2.725 2.70692 6.6831

56 2.775 1.73165 7.0406

57 2.825 1.07357 7.2703

58 2.875 0.594781 7.4021

59 2.925 0.360986 7.4849

60 2.975 0.251602 7.5446

61 3.025 0.198103 7.5932

62 3.075 0.196437 7.643

63 3.125 0.199002 7.6951

64 3.175 0.231627 7.7577

65 3.225 0.31308 7.845

66 3.275 0.322735 7.9378

67 3.325 0.4025 8.0571

68 3.375 0.491199 8.2071

69 3.425 0.548832 8.3797

70 3.475 0.613137 8.5782

71 3.525 0.681716 8.8053

72 3.575 0.748893 9.0619

73 3.625 0.791395 9.3407

74 3.675 0.799286 9.6301

75 3.725 0.880652 9.9577

…

Surprisingly, I obtained a better Si-O coordination number (3.99) when I used but LAMMPS but VMD to calculate g®.

Of course, I used the same output file… The reason is not clear.

Best regards.

Kousuke

11

Dear Kress

I am sorry. I mistyped.

Surprisingly, I obtained a better Si-O coordination number (3.99) when I used not LAMMPS but VMD to calculate g®.

Of course, I used the same output file… The reason is not clear.

Best regards.

Kousuke