reax/c gives shorted bond length for Si-O bonds

Hi all,

I am trying to simulate a Na2O-Al2O3-SiO2 glass with reax/c. As shown in the input file below - Equilibrating at 300K, gives me Si-O average bond length of 1.55-1.56 Angstrom instead of 1.61A, which is reported by many experimental studies. I have also run the code for 3000K which gave me shorter bond length than what is expected. I have also tried different timestep and thermostats. At this moment, I have a very bad glass structure. I am hoping for any suggestion that I can get. Thank you very much for your time.

Sincerely,
Ruhil Dongol

REAX potential for Soda-aluminosilicate glass system

units real
dimension 3
boundary p p p
atom_style charge
#read_data data.NaAlSi-1
read_restart restart.min1.NaAlSi1
#read_restart restart.300K.NaAlSi2

pair_style reax/c NULL
pair_coeff * * ffield.duin O Na Al Si
fix 99 all qeq/reax 1 0.0 10.0 1.0e-6 reax/c
neighbor 2.0 bin
neigh_modify delay 0 every 1 check yes
timestep 0.01
#dump 2 all xyz 1000 md.xyz
#dump 1 all xyz 1000 traj.xyz
thermo 100
thermo_style custom step etotal temp press pe vol lx dt time

#------------------------------300K-------------------------------
print “Equilibration at 300K”
reset_timestep 0
velocity all create 300 858939
fix 1 all npt temp 300 300 10 iso 0 0 10
timestep 0.01

dump 2 all xyz 1000 traj.300K.duin
run 400000
unfix 1
print “Equilibration Completed”

write_restart restart.300K.NaAlSi.duin

Hi all,

I am trying to simulate a Na2O-Al2O3-SiO2 glass with reax/c. As shown in the
input file below - Equilibrating at 300K, gives me Si-O average bond length
of 1.55-1.56 Angstrom instead of 1.61A, which is reported by many
experimental studies. I have also run the code for 3000K which gave me
shorter bond length than what is expected. I have also tried different
timestep and thermostats. At this moment, I have a very bad glass structure.
I am hoping for any suggestion that I can get. Thank you very much for your
time.

if a deviation of only 3.5% from experiment with a classical model
bothers you, perhaps you should not do simulations with such models.
please compare your results with what is documented for the force
field parameters for the conditions that you are simulating.

axel.

You will have to first check the Literature to see what Si-O bond length the particular force field you used give. Please realize that ReaxFF is an empirical potential that are fitted to QM/QC and/or experimental values. If it was fitted to a bond length of 1.55, then it will give a bond length of 1.55.

Thank you very much for taking time to reply.

Sincerely,
Ruhil