Generating amorphous silicon from quenching simulations

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:

 1 mkdir -p quench
 2 
 3 for n in 500 1000 2000 5000 10000 20000; do
 4 
 5 echo "Doing $n"
 6 
 7 cat>input<<eof
 8 ! Species-specific info
 9 atoms_file = 'cook.xyz'
10 pot_file = 'gap_files/silicon.gap'
11 n_species = 1
12 species = Si
13 masses = 28.09
14 
15 ! MD options
16 md_nsteps = $n
17 md_step = 2.
18 
19 ! Temperature
20 thermostat = berendsen
21 t_beg = 2000
22 t_end = 300
23 tau_t = 100.
24 
25 ! Pressure
26 barostat = berendsen
27 barostat_sym = iso
28 p_beg = 1.
29 p_end = 1.
30 gamma_p = 100.
31 tau_p = 1000.
32 
33 ! Writeouts
34 write_thermo = 1
35 write_xyz = $n
36 !write_lv = .true.
37 neighbors_buffer = 0.5
38 eof
39 
40 mpirun -np 4 turbogap md > /dev/null
41 
42 mv thermo.log quench/quench_$n.log
43 tail -218 trajectory_out.xyz > quench/quench_$n.xyz
44 
45 done

This workflow (which took a bit over two hours on my desktop) will run several quenches, each with a different quench rate, which you can easily calculate from the initial and final temperatures, the number of time steps and the time step:

[math]\displaystyle{ \text{quench 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.

Fully relax the structures

As a final step before analyzing the data, we will perform a geometry optimization of the structures using gradient descent:

 1 mkdir -p relax/
 2 
 3 for n in 500 1000 2000 5000 10000 20000; do
 4 
 5 echo "Doing $n"
 6 
 7 cat>input<<eof
 8 ! Species-specific info
 9 atoms_file = 'quench/quench_$n.xyz'
10 pot_file = 'gap_files/silicon.gap'
11 n_species = 1
12 species = Si
13 masses = 28.09
14 
15 ! MD options
16 optimize = gd
17 md_nsteps = 1000
18 
19 ! Writeouts
20 write_thermo = 1
21 write_xyz = $n
22 neighbors_buffer = 0.5
23 eof
24 
25 mpirun -np 4 turbogap md > /dev/null
26 
27 mv thermo.log relax/relax_$n.log
28 tail -218 trajectory_out.xyz > relax/relax_$n.xyz
29 
30 done

This script ran in 5.5 minutes on my desktop machine. After it has finished, directory relax/ will contain the different structures that we will use for analysis.

Analyzing the results

We will carry out basic radial distribution function calculations and coordination analysis of the a-Si structures generated using different quench rates. Let us prepare these two directories:

mkdir -p neighbors
mkdir -p rdf

After that, run the following Python script:

 1 from ase.io import read
 2 from ase.geometry.analysis import Analysis
 3 from ase.neighborlist import NeighborList
 4 import numpy as np
 5 
 6 neighbors_cutoff = 2.6
 7 rdf_cutoff = 5.
 8 
 9 for nsteps in [500, 1000, 2000, 5000, 10000, 20000]:
10     atoms = read("relax/relax_%i.xyz" % nsteps)
11     ana = Analysis(atoms)
12     cutoffs = np.zeros(len(atoms)) + neighbors_cutoff/2.
13     nl = NeighborList(cutoffs, skin=0., self_interaction=False, bothways=True)
14     nl.update(atoms)
15     cm = nl.get_connectivity_matrix()
16     coord = [cm[n].count_nonzero() for n in range(0,len(atoms))]
17     bins = np.bincount(coord)
18     rdf = ana.get_rdf(rdf_cutoff, 100, return_dists=True)
19 
20     f = open("neighbors/neighbors_%i.dat" % nsteps, "w+")
21     for i in range(0, len(bins)):
22         if bins[i] > 0:
23             print(i, bins[i]/len(atoms), file=f)
24 
25     f.close()
26 
27     f = open("rdf/rdf_%i.dat" % nsteps, "w+")
28     for i in range(0, len(rdf[0][0])):
29         print(rdf[0][1][i], rdf[0][0][i], file=f)
30 
31     f.close()

This will generate a series of files. In rdf/ you will find the radial distribution functions and in neighbors/ you will find coordination statistics: the fraction of atoms with 3, 4 and 5 neighbors. In silicon anything other than 4 neighbors is a coordination defect (unlike carbon, where 2, 3 and 4 can be stable motifs).

Let us first plot the coordination statistics:

We see that most atoms have 4 neighbors, as expected from the stable sp³ bonding in silicon. Because this is a-Si, there is a small but non-negligible number of coordination defects, atoms with both 3 and 5 neighbors. The number of defects goes down as we slow down the quench process, also as expected.

Now let us have a look at the radial distribution function:

There is a sharp peak at around 2.35 Angstrom, corresponding to the first-nearest neighbors distance. Then there is a wide gap and a much smoother peak for second-nearest and farther neighbors. In amorphous materials the structure in the RDF fades away relatively quickly (much quicker than in crystals) and the graph converges to one, which indicates an isotropic medium.

Finally, let us look at the energy per atom as a function of the quench rate:

Clearly, the system is allowed to find lower energy configurations as the quench rate is slowed down. From this graph we can tell that there is probably some extra stability to be gained by allowing for even slower quenching than the slowest one we tried (although note that the horizontal axis is logarithmic, so probably not as mush as may seem!).

Check the paper by Albert Bartók et al. on silicon simulation to compare the results of this tutorial.