Melting Quartz to obtain amorphous silica

Hi Lammps-users,

I have been trying to quench b-cristobalite from 5000K to 300K using the BKS potential but I keep getting an error of lost atoms. I obtained the atom coordinates from Wyckoff, and generated a system of 648 atoms. There is no such issues at lower temperatures of 3000K and below. I have tried different damping from 10.0 to 0.5 but the same thing occurred.

I tried to use the same methods as mentioned in the subject thread as well, to no avail. If anyone with experience doing this could kindly advise me on what is wrong with my method, it would be greatly appreciated.

Initialise

dimension 3
boundary p p p

units metal
atom_style charge

Create geometry

read_data silica.dat

neighbor 10.0 bin
neigh_modify every 1 delay 0 check yes

Potentials

kspace_style pppm 0.00001

pair_style buck/coul/long 8.0

Si type 1, Oxy type 2

pair_coeff 1 1 0.0 0.10 0.0
pair_coeff 1 2 18003.7572 0.205205 133.5381
pair_coeff 2 2 1388.7730 0.362319 175.0

Initialise velocities

compute new all temp
compute 1 all pe/atom
compute 2 all ke/atom
compute 3 all stress/atom pair

velocity all create 5000.0 4928459 temp new

Thermostat

fix 1 all nvt temp 5000.0 5000.0 1.0

Run

timestep 0.0005
thermo 100

#Dump format: dump ID group-ID style N file args
dump 1 all custom 1000 *.dump x y z
dump_modify 1 sort id

#Number of steps to run
run 100000

fix 1 all nvt temp 5000.0 300.0 1.0

run 3000000

write_restart restart.*.heat

Regards,
Yeo Jingjie
Ph.D. Student

School of Mechanical and Aerospace Engineering
Nanyang Technological University, Singapore

This could be due to a common problem with BKS style potentials
at high temps, which is that the exp() core is not strong enough
to prevent hot atoms from getting too close to each other. The
common solution is to spline on a stiffer repulsion for small r.

Paul Crozier can give you details - he has done this for silica
with LAMMPS.

Steve

Yeo Jingjie,

Please see the attached example files and this thread: http://lammps.sandia.gov/threads/msg13332.html

Paul

quench.in (664 Bytes)

silica.tabulated (121 KB)

quench.data (109 KB)

Hi Paul,

Thank you very much for the files. In the thread you mentioned that the BKS potential has "a different steep potential to avoid the "fusion" problem". May I know what was the modification that was done to the potential?

Regards,
Yeo Jingjie
Ph.D. Student
School of Mechanical and Aerospace Engineering
Nanyang Technological University, Singapore

Another problem that I ran into, after using the files you provided to quench, I attempted to change the fix nph and fix langevin to fix nve to do a thermal conductivity test, but upon changing to fix nve, the system began to move in one direction, in this case in the negative z-direction. May I know why is this so? Is there some special consideration I should account for?

Regards,
Yeo Jingjie
Ph.D. Student
School of Mechanical and Aerospace Engineering
Nanyang Technological University, Singapore

You must have net momentum in that dir.

Steve

Paul can answer this, but I think he already gave you pointers
to where this is discussed.

Steve

A steep repulsive potential was added such that at small separation distance, atoms repel instead of attract. This gets around the problem with the "fusion problem" with the BKS potential. For the force/energy vs distance curve specific numbers, take a look at the tabulated potential file that I sent. If you plot the data, it should be pretty obvious where the splice occurs.

Paul

Hi Steve and Paul,

Thank you very much for your patient replies. I believe the steep repulsive potential is similar in spirit to overlaying a 24-6 LJ potential? And I had no idea fix langevin could cause momentum drifts, and I believe this is remedied via fix momentum. This discussion has been incredibly helpful, thank you guys once again.

Regards,
Yeo Jingjie
Ph.D. Student
School of Mechanical and Aerospace Engineering
Nanyang Technological University, Singapore