[lammps-users] Melting Quartz to obtain amorphous silica

LAMMPS LAMMPS Mailing List Mirror
1

Dear Paul,

You provided me your input script for melting quartz two weeks ago. Thank you again.
I wanted to ask you please:
After I perform the melt with :

velocity all create 5000.0 4928459 dist gaussian
fix 1 all nph aniso 1.0 1.0 1.0 1.0 1.0 1.0 1000.0
fix 2 all langevin 5000.0 300.0 33.333 48279

Is there a specific way I should make sure to follow to cool it down?

Should I use minimize and/or fix box/relax at the same time or after cooling it down?

Thank you very much in advance for this help

Thomas Coquil
Research Assistant
UCLA Mechanical and Aerospace Engineering Department
420 Westwood Plaza
43-132 ENGINEERING IV
Los Angeles CA, 90095
Tel: 1 626 429 2110
thomas.coquil@…24…

2

Hi Thomas. The script I sent earlier is a recipe for quenching, so it already has the cooling procedure included. Specifically, take a look at the fix langevin command and corresponding documentation, and note that the simulation starts at 5000 K and ends at 300 K.

units real

atom_style full

pair_style table linear 1000

kspace_style pppm 0.0001

kspace_modify gewald 0.29

bond_style none

angle_style none

dihedral_style none

read_data quench.data

pair_coeff 1 1 silica.tabulated BKS_1_1

pair_coeff 1 2 silica.tabulated BKS_1_2

pair_coeff 2 2 silica.tabulated BKS_2_2

special_bonds 0.0 0.0 0.5

velocity all create 5000.0 4928459 dist gaussian

fix 1 all nph aniso 1.0 1.0 1.0 1.0 1.0 1.0 1000.0

fix 2 all langevin 5000.0 300.0 33.333 48279

timestep 1.0

thermo 500

dump 1 all atom 500000 quench.dump

run 500000

So this is the quenching procedure that I’ve used before and am recommending. No need to use minimize or fix box/relax, but you could certainly do a slower quench (by increasing the run time) if you have the compute resource and patience. Slower quenching would probably be more realistic.

Paul

3

Dear Paul,

I have been using the input script you provided to melt and quench quartz. It has been working perfectly when I use the same 5000K temperature (Tstart of fix langevin) as you did.

However when I set the Tstart of fix langevin to 7000 or 10000K for instance (I would like to start at 10000K ultimately), I get the following error:

“ERROR: Out of range atoms - cannot compute PPPM”

It happens after 80 000 steps (0.9 fs timestep) or more (sometimes up to 300 000 timesteps).

This happens although I set the neigh_modify command to delay 0 every 1 check yes:

kspace_style pppm 0.0001
neighbor 10.0 bin
neigh_modify every 1 delay 0 check yes

My timestep (<1fs) should also be fine.

I believe my pair interactions should be fine too. I use the same ones as previously reported byMccaughey and Kaviany :
(A.J.H. McGaughey, M. Kaviany / International Journal of Heat and Mass Transfer 47 (2004) 1799–1816)
And everything seemed to work fine at lower temperatures.

Question:
So is it that I need to make it more viscous and reduce the damp parameter of the fix langevin when working at larger temperature? If so, how should I evaluate the required value to be used?

Is there anything to do about the fix nph and its damp parameter also? How should this one be chosen?

I would truly appreciate your help once again if you had any ideas on this.

Thank you very much in advance

Sincerely,

Thomas Coquil

p.s: here’s my input script for your information:

#initialization############################################################

units metal
dimension 3
boundary p p p
atom_style charge

Atom definition##########################################################

lattice custom 5.4054 &
a1 0.9095 0.0000 0.0000 &
a2 -0.4547 0.7876 0.0000 &
a3 0.0000 0.0000 1.0000 &
basis 0.4697 0.0000 0.0000 basis 0.0000 0.4697 0.6667 basis 0.5303 0.5303 0.3333 basis 0.4135 0.2669 0.1191 basis 0.2669 0.4135 0.5475 basis 0.7331 0.1466 0.7857 basis 0.5865 0.8534 0.2142 basis 0.8534 0.5865 0.4524 basis 0.1466 0.7331 0.8809

region simbox block 0 4 0 4 0 16 units lattice

create_box 2 simbox
create_atoms 1 box basis 1 1 basis 2 1 basis 3 1 &
basis 4 2 basis 5 2 basis 6 2 &
basis 7 2 basis 8 2 basis 9 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

Atoms interactions settings##################################

Si type 1, O type 2

pair_style hybrid/overlay table linear 3000 coul/long 8.0
pair_coeff 1 1 table BKSLJ.table SiSi 8.0
pair_coeff 1 2 table BKSLJ.table SiO 8.0
pair_coeff 2 2 table BKSLJ.table OO 8.0
pair_coeff * * coul/long

kspace_style pppm 0.0001
neighbor 10.0 bin
neigh_modify every 1 delay 0 check yes

pair_write 1 2 1000 rsq 0.01 10.0 table.txt BKS_Si_O 2.4 -1.2
pair_write 2 2 1000 rsq 0.01 10.0 table.txt BKS_O_O -1.2 -1.2
pair_write 1 1 1000 rsq 0.01 10.0 table.txt BKS_Si_Si 2.4 2.4

velocity all create 10000.0 4928459 dist gaussian
timestep 0.000905
thermo_style custom step temp press &
ke pe etotal vol lx ly lz epair

# Quench in NPH with Langevin ##############################

fix 2 all nph aniso 1.01325 1.01325 1.01325 1.01325 1.01325 1.01325 1.0
fix 3 all langevin 10000.0 300.0 0.033333 48279

thermo 2000
dump 2 all atom 10000 quenchPROCESS.dump
dump 3 all atom 3000000 quenchFINAL.dump
run 3000000
write_restart restart.amorF4x4x12

4

Thomas,

Are you sure your initial configuration does not have any overlapping atoms? The error message you are seeing typically means that either atoms are overlapping or huge forces are being computed somewhere. Check this by visualizing your configuration, dumping forces, dumping pressures, etc.

I’d recommend using larger damp values in your fix nph and fix langevin commands. The values you have there look rather small. Note they are in units of time, and that longer time means smoother action. Please carefully read the documentation here:

http://lammps.sandia.gov/doc/fix_langevin.html

http://lammps.sandia.gov/doc/fix_nh.html

There’s no exact “correct” value to use for these damping parameters, but you do need to be in the right ballpark. I don’t know of a magic formula for choosing a good value, but as a rule of thumb you might start with around 10000 timesteps’ worth (in time units).

A tighter kspace tolerance might help:

kspace_style pppm 0.00001

One more idea: you might consider trying fix dt/reset so that LAMMPS automatically chooses an appropriate timestep size for you. At such high and rapidly changing temperatures, it might make sense to let LAMMPS figure this out for you automatically and on the fly. Something like this would probably work for you:

fix 9 all dt/reset 1 0.001 2.0 0.1 units box

See: http://lammps.sandia.gov/doc/fix_dt_reset.html

Paul

5

Dear Paul,

Thank you very much, it seems like it is working with larger damping coefficients indeed.

I am just wondering. By doing this… Should the density of the final amorphous phase automatically adjust itself to that of amorphous silica?

I mean I am starting from a quartz phase with a density of about 2.6. Should that density naturally drop to about 2.2 for amorphous silica after the quench in NPH or is it not that straightforward at all?

When I quench it with your input script at 5000K for 3M steps of 0.9fs, I get a final phase of density 2.5… not quite 2.2 yet.

Should I rather start with a crystalline phase of density already modified to 2.2 to get closer to the real amorphous phase density?

what would be more realistic you think?

Thank you very much for all your insight on this.

Best regards,

Thomas Coquil
Research Assistant
UCLA Mechanical and Aerospace Engineering Department
420 Westwood Plaza
43-132 ENGINEERING IV
Los Angeles CA, 90095
Tel: 1 626 429 2110
thomas.coquil@…24…