Simulated density of amorphous silica is too high

LAMMPS LAMMPS General Discussion
1

Hello everyone,

I am a beginner with LAMMPS and am currently learning to create an amorphous silica atomic model using LAMMPS.

I am using the BKS potential because it is widely used in academic literature, and I have applied some corrections using a 24‑6 Lennard‑Jones potential. The simulation runs without any errors in LAMMPS and produces a final configuration.

However, my target density is 2.2 g/cm³, but the final configuration shows a density around 2.5–2.6 g/cm³, which deviates significantly from my actual sample.

So I would like to ask the community for advice: where could the problem lie?

Is it the BKS potential itself?

The short‑range Buckingham cutoff distance?

The long‑range Coulomb cutoff distance?

The temperature damping coefficient or pressure damping coefficient during the melt‑quench process?

I would be very grateful for any valuable suggestions.

Below is my input file:

text
units metal
atom_style charge
boundary p p p

read_data SiO2_random_L21.63_N669.data

pair_style hybrid/overlay buck/coul/long 10.0 5.5 table linear 10000
pair_coeff 1 1 buck/coul/long 1388.7730 0.3623 175.0
pair_coeff 1 2 buck/coul/long 18003.7572 0.2052 133.5381
pair_coeff 2 2 buck/coul/long 0.0 1.0 0.0
pair_coeff 1 1 table table_24-6.txt OO
pair_coeff 1 2 table table_24-6.txt SiO
pair_coeff 2 2 table table_24-6.txt SiSi
#pair_modify shift yes
kspace_style pppm 1.0e-4

neighbor 2.0 bin
neigh_modify every 1 delay 0 check yes
group O type 1
group Si type 2

---------- 1. Energy minimization ----------

min_style cg
minimize 1.0e-6 1.0e-8 1000 10000

---------- 2. Initialize velocities at 300 K ----------

velocity all create 300 4928459 dist gaussian
timestep 0.001 # 1 fs

---------- 3. Heat from 300 K to 5000 K in NPT (100 ps) ----------

fix heat all npt temp 300 5000 0.1 iso 0.0 0.0 1.0
thermo 1000
thermo_style custom step temp press vol density pe ke enthalpy
dump heatdump all custom 5000 dump.heat_669.lammpstrj id type x y z
run 100000 # 100 ps

---------- 4. NPT equilibration at 5000 K for 100 ps ----------

unfix heat
fix melt all npt temp 5000 5000 0.1 iso 0.0 0.0 1.0
dump meltdump all custom 5000 dump.melt_669.lammpstrj id type x y z
run 100000

---------- 5. Cooling from 5000 K to 300 K at 1 K/ps ----------

unfix melt
fix cool all npt temp 5000 300 0.1 iso 0.0 0.0 1.0
dump cooldump all custom 5000 dump.cool_669.lammpstrj id type x y z
dump_modify cooldump first yes
run 4700000

---------- 6. Final relaxation at 300 K for 100 ps ----------

unfix cool
fix relax all npt temp 300 300 0.1 iso 0.0 0.0 1.0
dump relaxdump all custom 5000 dump.relax_669.lammpstrj id type x y z
dump_modify relaxdump first yes
run 100000

---------- Save final structure ----------

unfix relax
write_data final_glass_669.data
write_restart restart.glass_669.restart

print “Simulation completed for 669 atoms”

Thank you in advance for your help!

2

My guess is that the intermolecular interactions are too strong, since pair_style hybrid/overlay sums the two pair potentials. For instance, plot the sum of the Buckingham and table potentials for each pair of atoms and check its asymptotic behaviour. If it is even slightly more negative, here is the explanation for the higher density.
You can then decide what to do next:

  1. Fit the parameters of a function f(x) which is the sum of the Buckingham potential + correction, using the original potential as reference. Repeat the melt-quench simulation with the corrected potential. Or:
  2. Do a second thermal annealing with the original Buckingham potential, perhaps 900 K will heal the structure.

Let us know how it goes, for future reference.