Hi all, I’m trying to generate amorphous silica using the melt / quench method and Vashishta potential. I’m aiming for a final density of ~2.2 g/cc, however regardless of quenching rate, the final density ends up being around 2.6 g/cc. I’ve tried rates of 10, 20, 50, 100, and 1000 K/ps and gotten nearly identical results. The final structure is amorphous though. LAMMPS version 27 Oct 2021.
Here is my process flow:
- Specify positions of 24 basis atoms for structure → replicate to create 5x5x5 supercell
- NPT at 300 K to relax structure (I’ve also tried starting the simulation immediately at 4000 K and achieved the same results)
- Ramp to 4000 K at a rate of 20 K/ps
- Hold at 4000 K for 500 ps (density of ~2.2 g/cc)
- NPT quench at rates specified above (10, 20, 50, 100, 1000 K/ps) (also tried NVT to fix density → relaxed with NPT at 300 K afterwards for 50 ps and final density was 2.4 g/cc… any comment on this could be useful as well)
- NPT at 300 K to relax structure (density of ~2.6 g/cc)
It is quite strange to me that the final densities are so similar even with such drastically different quenching rates. I’m not going to point my finger at the potential or LAMMPS as there must be something I’m doing wrong here… Changing Pdamp from 0.1 to 1 did not change much either (timestep = 0.1 fs). I’ve pasted my files to this thread.
Input script:
# Initial Structure: b-cristobalite
# Initial Lattice Parameter: a = 7.16 Angstroms
# Supercell Size: 5x5x5
# Heating Rate: 20 K/ps
# Melt for: 500 ps
# Quench Rate: 10 K/ps
# Thermalize for: 50 ps
# ------------------------ INITIALIZATION ----------------------------
units metal
dimension 3
boundary p p p
atom_style charge
pair_style vashishta
neighbor 1.0 bin
neigh_modify delay 1
# beta-cristobalite SiO2
variable latparam equal 7.16
# ----------------------- SIMULATION BOX DEFINITION ----------------------------
read_data b-cristobalite.data
# ----------------------- REPLICATE BASIS ATOMS ----------------------------
replicate 5 5 5
# ------------------------ SETTINGS ------------------------------
pair_coeff * * SiO.1990.vashishta.txt O Si
# ----------------------- INITIALIZE STRUCTURE AT T = 300 K ----------------------------
reset_timestep 0
timestep 0.0001
velocity all create 300 123456 mom yes rot yes
fix initialize all npt temp 300 300 0.01 iso 0 0 0.1
# Set thermo output for initializing at T = 300 K
thermo 500
thermo_style custom step temp press pe ke etotal vol density
# Dump images
#removed for simplicity
# 10 ps
run 100000 # energy and density plateau starting at ~3.5 ps
unfix initialize
undump 1
# ----------------------- HEAT FROM T = 300 K TO T = 2150 K ----------------------------
reset_timestep 0
timestep 0.0001
fix heat1 all npt temp 300 2150 0.01 iso 0 0 1
# Set thermo output for heating from T = 300 K to T = 2150 K
thermo 500
thermo_style custom step temp press pe ke etotal vol density
# Dump images
#removed for simplicity
# 92.5 ps
run 925000 # 20 K/ps heating rate
unfix heat1
undump 1
# ----------------------- HEAT FROM T = 2150 K TO T = 4000 K ----------------------------
reset_timestep 0
timestep 0.0001
fix heat2 all npt temp 2150 4000 0.01 iso 0 0 1
# Set thermo output for heating from T = 2150 K to T = 4000 K
thermo 500
thermo_style custom step temp press pe ke etotal vol density
# Dump images
#removed for simplicity
# 92.5 ps
run 925000 # 20 K/ps heating rate
unfix heat2
undump 1
# ----------------------- MELT AT T = 4000 K ----------------------------
reset_timestep 0
timestep 0.0001
fix melt all npt temp 4000 4000 0.01 iso 0 0 1
# Set thermo output for melting at T = 4000 K
thermo 500
thermo_style custom step temp press pe ke etotal vol density
# Dump images
#removed for simplicity
# 500 ps
run 5000000
unfix melt
undump 1
# ----------------------- QUENCH FROM T = 4000 K to T = 300 K ----------------------------
reset_timestep 0
timestep 0.0001
#fix quench all nvt temp 4000 300 0.01
fix quench all npt temp 4000 300 0.01 iso 0 0 1
# Set thermo output for quenching from T = 4000 K to T = 300 K
thermo 500
thermo_style custom step temp press pe ke etotal vol density
# Dump images
#removed for simplicity
# 370 ps
run 3700000 # 10 K/ps quench rate
unfix quench
undump 1
# ----------------------- THERMALIZE AT T = 300 K (NPT TO RELAX STRUCTURE) ----------------------------
reset_timestep 0
timestep 0.0001
fix thermalize all npt temp 300 300 0.01 iso 0 0 1
# Set thermo output for thermalization at T = 300 K
thermo 500
thermo_style custom step temp press pe ke etotal vol density
# Dump images
#removed for simplicity
# 50 ps
run 500000
unfix thermalize
undump 1
b-cristobalite.data file:
LAMMPS Description
24 atoms
2 atom types
0.0 7.16 xlo xhi
0.0 7.16 ylo yhi
0.0 7.16 zlo zhi
Masses
1 15.9994 # O
2 28.0000 # Si
Atoms
1 2 2.4 0.0000 0.0000 0.0000
2 2 2.4 0.0000 3.5800 3.5800
3 2 2.4 3.5800 3.5800 0.0000
4 2 2.4 3.5800 0.0000 3.5800
5 2 2.4 5.3700 1.7900 5.3700
6 2 2.4 1.7900 1.7900 1.7900
7 2 2.4 1.7900 5.3700 5.3700
8 2 2.4 5.3700 5.3700 1.7900
9 1 -1.2 0.8950 0.8950 0.8950
10 1 -1.2 6.2650 2.6850 4.4750
11 1 -1.2 2.6850 4.4750 6.2650
12 1 -1.2 4.4750 6.2650 2.6850
13 1 -1.2 2.6850 6.2650 4.4750
14 1 -1.2 6.2650 4.4750 2.6850
15 1 -1.2 4.4750 2.6850 6.2650
16 1 -1.2 0.8950 4.4750 4.4750
17 1 -1.2 6.2650 6.2650 0.8950
18 1 -1.2 2.6850 0.8950 2.6850
19 1 -1.2 2.6850 2.6850 0.8950
20 1 -1.2 6.2650 0.8950 6.2650
21 1 -1.2 4.4750 0.8950 4.4750
22 1 -1.2 0.8950 6.2650 6.2650
23 1 -1.2 0.8950 2.6850 2.6850
24 1 -1.2 4.4750 4.4750 0.8950
Vasishta potential:
# DATE: 2015-10-14 CONTRIBUTOR: Aidan Thompson, [email protected] CITATION: P. Vashishta, R. K. Kalia, J. P. Rino, and I. Ebbsjo, Phys. Rev. B 41, 12197 (1990).
#
# Vashishta potential file for SiO2, P. Vashishta, R. K. Kalia, J. P. Rino, and I. Ebbsjo,
# Phys. Rev. B 41, 12197 (1990).
#
# These parameters, some inferred indirectly, give a good
# match to the energy-volume curve for alpha-quartz in Fig. 2 of the paper.
#
# These entries are in LAMMPS "metal" units:
# H = eV*Angstroms^eta; Zi, Zj = |e| (e = electronic charge);
# lambda1, lambda4, rc, r0, gamma = Angstroms;
# D = eV*Angstroms^4; W = eV*Angstroms^6; B = eV;
# other quantities are unitless
#
# element1 element2 element3
# H eta Zi Zj lambda1 D lambda4
# W rc B gamma r0 C cos(theta)
Si Si Si 0.82023 11 1.6 1.6 999 0.0 4.43
0.0 10.0 0.0 0.0 0.0 0.0 0.0
O O O 743.848 7 -0.8 -0.8 999 22.1179 4.43
0.0 10.0 0.0 0.0 0.0 0.0 0.0
O Si Si 163.859 9 -0.8 1.6 999 44.2357 4.43
0.0 10.0 20.146 1.0 2.60 0.0 -0.77714596
Si O O 163.859 9 1.6 -0.8 999 44.2357 4.43
0.0 10.0 5.0365 1.0 2.60 0.0 -0.333333333333
Si O Si 0.0 0.0 0.0 0.0 0.0 0.0 0.0
0.0 0.0 0.0 0.0 0.0 0.0 0.0
Si Si O 0.0 0.0 0.0 0.0 0.0 0.0 0.0
0.0 0.0 0.0 0.0 0.0 0.0 0.0
O Si O 0.0 0.0 0.0 0.0 0.0 0.0 0.0
0.0 0.0 0.0 0.0 0.0 0.0 0.0
O O Si 0.0 0.0 0.0 0.0 0.0 0.0 0.0
0.0 0.0 0.0 0.0 0.0 0.0 0.0
Thanks in advance.