Flying Ice cube

LAMMPS LAMMPS General Discussion
1

Hi…I perform the standard melt quench technique for obtaining amorphous silica. But for some reason I see the flying ice cube effect right after my first equilibration step at 300K. I have tried reducing the timestep, modifying the initial structure but I still get the same flying ice cube effect. I don’t understand where I am going wrong.

Any help is much appreciated. Thanks in advance.

This is my code

# ---------------INITIALIZATION------------------------------
units metal
dimension 3
boundary p p p
atom_style charge

# ---------------STRUCTURE DEFINITION------------------------
lattice custom 5.067922 a1 1.000000 0.000000 0.000000 a2 0.000000 1.000000 0.000000 a3 0.000000 0.000000 1.393663 &
                        basis 0.000000 0.500000 0.250000 basis 0.000000 0.000000 0.000000 basis 0.500000 0.000000 0.750000 basis 0.500000 0.500000 0.500000 &
                        basis 0.250000 0.412034 0.375000 basis 0.087966 0.750000 0.125000 basis 0.750000 0.587966 0.375000 basis 0.912034 0.250000 0.125000 &
                        basis 0.750000 0.912034 0.875000 basis 0.587966 0.250000 0.625000 basis 0.250000 0.087966 0.875000 basis 0.412034 0.750000 0.625000

region simbox block 0.0 7.0 0.0 7.0 0.0 7.0 

create_box 2 simbox
create_atoms 1 box basis 1 1 basis 2 1 basis 3 1 basis 4 1 &
                   basis 5 2 basis 6 2 basis 7 2 basis 8 2 &
                   basis 9 2 basis 10 2 basis 11 2 basis 12 2
mass 1 28.0855
mass 2 15.9994

group siliconatoms type 1
group oxygenatoms type 2

set group siliconatoms charge 2.4
set group oxygenatoms charge -1.2


# ---------------TEMPERATURE AND TIMESTEP--------------------

velocity all create 300.0 12345 dist gaussian
timestep 0.002

thermo 10000
thermo_style custom step etotal pe ke temp density vol lx ly lz


# ---------------INTERATOMIC POTENTIAL-----------------------
pair_style comb
pair_coeff * * ffield.comb Si O
neighbor 2.0 bin
neigh_modify every 10 delay 0 check yes
fix charge all qeq/comb 100 0.05 file qeq.comb


# ---------------EQUILIBRATION-------------------------
dump equib all atom 10000 equib.atom
fix equilibration all nvt temp 300.0 300.0 $(100.0*dt)
run 100000
unfix equilibration

minimize 0.0 0.0 100000 100000
undump equib

# ---------------MELTING-------------------------------------

dump melt all atom 10000 melt.atom
fix melting1 all nvt temp 300.0 1000.0 $(100.0*dt) 
run 87500
unfix melting1
fix melting2 all nvt temp 1000.0 2000.0 $(100.0*dt)
run 125000 
unfix melting2
fix melting3 all nvt temp 2000.0 3000.0 $(100.0*dt) 
run 125000
unfix melting3



fix equilibration all nvt temp 3000.0 3000.0 $(100.0*dt)
run 100000
unfix equilibration

minimize 0.0 0.0 100000 100000
undump melt

# ---------------QUENCHING-----------------------------------

dump quench all atom 10000 quench.atom


fix quenching4 all nvt temp 3000.0 2000.0 $(100.0*dt) 
run 500000
unfix quenching4
fix quenching5 all nvt temp 2000.0 1000.0 $(100.0*dt) 
run 500000
unfix quenching5
fix quenching6 all nvt temp 1000.0 300.0 $(100.0*dt) 
run 350000
unfix quenching6


minimize 0.0 0.0 100000 100000

fix equilibration all nvt temp 300.0 300.0 $(100.0*dt)
run 100000
unfix equilibration
undump quench

write_data final_amorphous_bulk_SiO2.dat
write_restart final_amorphous_bulk_SiO2.restart
2

You can use the velocity command to more rapidly cool down the system by using:
velocity all set 0.0 0.0 0.0
The best place for that would be after each minimize command.
Or after reequilibration after the quench you can use the velocity zero linear command to remove any center of mass momentum that has developed. Typically that would only be needed once after the system has been quenched and reequilibrated. If this keeps happening, you may need to try shortening the timestep, and if that doesn’t help either, you may consider using fix momentum.

3

ok …Thank you so much. I will try these out.