Generating amorphous silicon from quenching simulations

WARNING: This tutorial is under construction!!!!!!!!!!!

In this tutorial we will run TurboGAP simulations of silicon amorphization using molecular dynamics and geometry optimization. Check the references to know more about the Si potential we are using and atomistic simulation of silicon in general.

Prerequisites for this tutorial

• A TurboGAP installation
• An ASE installation
• Numpy

Optional

• gnuplot (for plotting)
• VMD (for visualization, ASE can handle visualization but it's slow)

Setting up the initial configuration

We are going to create a diamond lattice with 216 atoms and give the atoms random velocities. We will use ASE for this:

 1 from ase.io import write
 2 from ase import Atoms
 3 import numpy as np
 4 
 5 # Simple cubic lattice parameter
 6 a = 5.43/2.
 7 
 8 atoms = Atoms("Si", cell = [a,a,a], positions=[[0,0,0]], pbc=True)
 9 atoms *= (6, 6, 6)
10 
11 vel = 0.01 * (np.random.sample([len(atoms), 3]) - 0.5)
12 atoms.set_array("velocities", vel)
13 
14 write("lattice.xyz", atoms)

lattice.xyz contains the atomic data.

Melting and barostating

Before running TurboGAP, make sure that you have downloaded the GAP files. In your working directory, do:

wget https://zenodo.org/record/5734463/files/Si_PW91.tar.gz
tar -xvf Si_PW91.tar.gz
mv Si_PW91/gap_files .

Now you have a gap_files/ directory with all the necessary files to run a GAP simulation with TurboGAP.

We are going to melt the simple cubic structure at 5000 K for 4 ps (with time step 2 fs). We choose the Berendsen thermostat with time constant 100 fs. We are also going to use the Berendsen barostat with time constant 1000 fs and inverse compressibility 100 times larger than water, with a target pressure of 1 bar. You need the following input file:

! Species-specific info
atoms_file = 'lattice.xyz'
pot_file = 'gap_files/silicon.gap'
n_species = 1
species = Si
masses = 28.09
!
! MD options
md_nsteps = 2000
md_step = 2.
!
! Temperature
thermostat = berendsen
t_beg = 5000
t_end = 5000
tau_t = 100.
!
! Pressure
barostat = berendsen
barostat_sym = iso
p_beg = 1.
p_end = 1.
gamma_p = 100.
tau_p = 1000.
!
! Writeouts
write_thermo = 1
write_xyz = 1000
!write_lv = .true.
neighbors_buffer = 0.5

Run TurboGAP to generate the Si liquid:

mpirun -np 4 turbogap md
tail -218 trajectory_out.xyz > melt.xyz

On my (oldish) desktop machine, this calculation as run above (4 MPI processes) took a bit over six minutes.

Equilibrating at high temperature

We will equilibrate the Si system at 2000 K, which will enable the formation of reasonably stable but non-crystalline motifs:

! Species-specific info
atoms_file = 'melt.xyz'
pot_file = 'gap_files/silicon.gap'
n_species = 1
species = Si
masses = 28.09
!
! MD options
md_nsteps = 50000
md_step = 2.
!
! Temperature
thermostat = berendsen
t_beg = 2000
t_end = 2000
tau_t = 100.
!
! Pressure
barostat = berendsen
barostat_sym = iso
p_beg = 1.
p_end = 1.
gamma_p = 100.
tau_p = 1000.
!
! Writeouts
write_thermo = 1
write_xyz = 50000
!write_lv = .true.
neighbors_buffer = 0.5

And run:

mpirun -np 4 turbogap md
tail -218 trajectory_out.xyz > cook.xyz

This will be the most expensive part of this tutorial, since we are running 50k time steps. On my desktop machine, it took a bit over 2.5 hours with 4 MPI processes. You can speed it up on a computer cluster (it should scale up nicely within a single node; for multinode execution the number of atoms may be a bit too small to yield good speed up).

Quenching rates

To make a-Si we need to quench the high-temperature equilibrated precursor to low temperature. The quench rate (usually given in terms of degrees Kelvin per ps) will determine the quality of the resulting structure. Usually, the slower the better. We can set up a bash workflow like this:

mkdir -p quench
for n in 500 1000 2000 5000 10000 20000; do

echo "Doing $n"

cat>input<<eof
! Species-specific info
atoms_file = 'cook.xyz'
pot_file = 'gap_files/silicon.gap'
n_species = 1
species = Si
masses = 28.09

! MD options
md_nsteps = $n
md_step = 2.

! Temperature
thermostat = berendsen
t_beg = 2000
t_end = 300
tau_t = 100.

! Pressure
barostat = berendsen
barostat_sym = iso
p_beg = 1.
p_end = 1.
gamma_p = 100.
tau_p = 1000.

! Writeouts
write_thermo = 1
write_xyz = $n
!write_lv = .true.
neighbors_buffer = 0.5
eof

mpirun -np 4 turbogap md > /dev/null

mv thermo.log quench/quench_$n.log
tail -218 trajectory_out.xyz > quench/quench_$n.xyz

done

This workflow will run several quenches, each with a different quench rate, which you can easily calculated from the initial and final temperatures, the number of time steps and the time step:

[math]\displaystyle{ \text{rate} = \frac{T_\text{end} - T_\text{beg}}{\Delta t n_\text{steps}}. }[/math]

In this case, we will be varying the quench rate between 1700 K/ps and 42.5 K/ps. We will use each final structure to study the effect of quench rate on the structure of the a-Si.