Generating amorphous silicon from quenching simulations - Revision history

Miguel Caro: /* Analyzing the results */

2022-03-21T08:22:37Z

Analyzing the results

Miguel Caro at 08:13, 21 March 2022

2022-03-21T08:13:27Z

Miguel Caro: /* Analyzing the results */

2022-03-21T08:13:08Z

Analyzing the results

Miguel Caro: /* Analyzing the results */

2022-03-21T06:25:59Z

Analyzing the results

Miguel Caro: /* Fully relax the structures */

2022-03-21T06:21:44Z

Fully relax the structures

Miguel Caro: /* Quenching rates */

2022-03-21T06:20:50Z

Quenching rates

Miguel Caro: /* Fully relax the structure */

2022-03-21T06:17:53Z

Fully relax the structure

Miguel Caro: /* Quenching rates */

2022-03-21T06:14:22Z

Quenching rates

Miguel Caro: /* Quenching rates */

2022-03-20T18:36:24Z

Quenching rates

Miguel Caro: /* Quenching rates */

2022-03-20T18:12:36Z

Quenching rates

← Older revision		Revision as of 08:22, 21 March 2022
Line 288:		Line 288:

	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!).		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.<ref name="bartok_2018" />
		+
		+	{{Reference list}}

← Older revision		Revision as of 08:13, 21 March 2022
Line 1:		Line 1:
−	~~'''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.		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.

← Older revision		Revision as of 08:13, 21 March 2022
Line 274:		Line 274:

	Let us first plot the coordination statistics:		Let us first plot the coordination statistics:
		+
		+	[[File:tutorial3_1.png\|center\|480px\|Fraction of atoms as a function of their number of neighbors.]]
		+
		+	We see that most atoms have 4 neighbors, as expected from the stable sp<sup>3</sup> 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:
		+
		+	[[File:tutorial3_2.png\|center\|480px\|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:
		+
		+	[[File:tutorial3_3.png\|center\|480px\|Energy per atom for different rates with respect to the slowest rate result.]]
		+
		+	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!).

← Older revision		Revision as of 06:25, 21 March 2022
Line 236:		Line 236:

	After that, run the following Python script:		After that, run the following Python script:
		+
		+	<syntaxhighlight lang="python" line="line">
		+	from ase.io import read
		+	from ase.geometry.analysis import Analysis
		+	from ase.neighborlist import NeighborList
		+	import numpy as np
		+
		+	neighbors_cutoff = 2.6
		+	rdf_cutoff = 5.
		+
		+	for nsteps in [500, 1000, 2000, 5000, 10000, 20000]:
		+	atoms = read("relax/relax_%i.xyz" % nsteps)
		+	ana = Analysis(atoms)
		+	cutoffs = np.zeros(len(atoms)) + neighbors_cutoff/2.
		+	nl = NeighborList(cutoffs, skin=0., self_interaction=False, bothways=True)
		+	nl.update(atoms)
		+	cm = nl.get_connectivity_matrix()
		+	coord = [cm[n].count_nonzero() for n in range(0,len(atoms))]
		+	bins = np.bincount(coord)
		+	rdf = ana.get_rdf(rdf_cutoff, 100, return_dists=True)
		+
		+	f = open("neighbors/neighbors_%i.dat" % nsteps, "w+")
		+	for i in range(0, len(bins)):
		+	if bins[i] > 0:
		+	print(i, bins[i]/len(atoms), file=f)
		+
		+	f.close()
		+
		+	f = open("rdf/rdf_%i.dat" % nsteps, "w+")
		+	for i in range(0, len(rdf[0][0])):
		+	print(rdf[0][1][i], rdf[0][0][i], file=f)
		+
		+	f.close()
		+	</syntaxhighlight>
		+
		+	This will generate a series of files. In <code>rdf/</code> you will find the radial distribution functions and in <code>neighbors/</code> 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:

@@ Line 193: / Line 193: @@
 As a final step before analyzing the data, we will perform a geometry optimization of the structures using [[gradient descent]]:
- mkdir -p relax/
+<syntaxhighlight lang="bash" line="line">
 mkdir -p relax/
- for n in 500 1000 2000 5000 10000 20000; do
 for n in 500 1000 2000 5000 10000 20000; do
- echo "Doing $n"
 echo "Doing $n"
- cat>input<<eof
- ! Species-specific info
+cat>input<<eof
- atoms_file = 'quench/quench_$n.xyz'
+! Species-specific info
- pot_file = 'gap_files/silicon.gap'
+atoms_file = 'quench/quench_$n.xyz'
- n_species = 1
+pot_file = 'gap_files/silicon.gap'
- species = Si
+n_species = 1
- masses = 28.09
+species = Si
 masses = 28.09
- ! MD options
- optimize = gd
+! MD options
- md_nsteps = 1000
+optimize = gd
 md_nsteps = 1000
- ! Writeouts
- write_thermo = 1
+! Writeouts
- write_xyz = $n
+write_thermo = 1
- neighbors_buffer = 0.5
+write_xyz = $n
- eof
+neighbors_buffer = 0.5
 eof
- mpirun -np 4 turbogap md > /dev/null
 mpirun -np 4 turbogap md > /dev/null
- mv thermo.log relax/relax_$n.log
- tail -218 trajectory_out.xyz > relax/relax_$n.xyz
+mv thermo.log relax/relax_$n.log
 tail -218 trajectory_out.xyz > relax/relax_$n.xyz
- done
 This script ran in 5.5 minutes on my desktop machine. After it has finished, directory <code>relax/</code> will contain the different structures that we will use for analysis.

@@ Line 133: / Line 133: @@
 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
+<syntaxhighlight lang="bash" line="line">
- for n in 500 1000 2000 5000 10000 20000; do
+mkdir -p quench
-−
- echo "Doing $n"
+for n in 500 1000 2000 5000 10000 20000; do
-−
- cat>input<<eof
+echo "Doing $n"
- ! Species-specific info
- atoms_file = 'cook.xyz'
+cat>input<<eof
- pot_file = 'gap_files/silicon.gap'
+! Species-specific info
- n_species = 1
+atoms_file = 'cook.xyz'
- species = Si
+pot_file = 'gap_files/silicon.gap'
- masses = 28.09
+n_species = 1
 species = Si
- ! MD options
+masses = 28.09
- md_nsteps = $n
- md_step = 2.
+! MD options
 md_nsteps = $n
- ! Temperature
+md_step = 2.
- thermostat = berendsen
- t_beg = 2000
+! Temperature
- t_end = 300
+thermostat = berendsen
- tau_t = 100.
+t_beg = 2000
 t_end = 300
- ! Pressure
+tau_t = 100.
- barostat = berendsen
- barostat_sym = iso
+! Pressure
- p_beg = 1.
+barostat = berendsen
- p_end = 1.
+barostat_sym = iso
- gamma_p = 100.
+p_beg = 1.
- tau_p = 1000.
+p_end = 1.
 gamma_p = 100.
- ! Writeouts
+tau_p = 1000.
- write_thermo = 1
- write_xyz = $n
+! Writeouts
- !write_lv = .true.
+write_thermo = 1
- neighbors_buffer = 0.5
+write_xyz = $n
- eof
+!write_lv = .true.
 neighbors_buffer = 0.5
- mpirun -np 4 turbogap md > /dev/null
+eof
-−
- mv thermo.log quench/quench_$n.log
+mpirun -np 4 turbogap md > /dev/null
- tail -218 trajectory_out.xyz > quench/quench_$n.xyz
 mv thermo.log quench/quench_$n.log
- done
+tail -218 trajectory_out.xyz > quench/quench_$n.xyz
 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:

@@ Line 186: / Line 186: @@
 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 structure ==
+== Fully relax the structures ==
 As a final step before analyzing the data, we will perform a geometry optimization of the structures using [[gradient descent]]:
@@ Line 221: / Line 221: @@
   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.
+This script ran in 5.5 minutes on my desktop machine. After it has finished, directory <code>relax/</code> will contain the different structures that we will use for analysis.

← Older revision		Revision as of 06:14, 21 March 2022
Line 185:		Line 185:

	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.		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 structure ==
		+
		+	As a final step before analyzing the data, we will perform a geometry optimization of the structures using [[gradient descent]]:
		+
		+	mkdir -p relax/
		+
		+	for n in 500 1000 2000 5000 10000 20000; do
		+
		+	echo "Doing $n"
		+
		+	cat>input<<eof
		+	! Species-specific info
		+	atoms_file = 'quench/quench_$n.xyz'
		+	pot_file = 'gap_files/silicon.gap'
		+	n_species = 1
		+	species = Si
		+	masses = 28.09
		+
		+	! MD options
		+	optimize = gd
		+	md_nsteps = 1000
		+
		+	! Writeouts
		+	write_thermo = 1
		+	write_xyz = $n
		+	neighbors_buffer = 0.5
		+	eof
		+
		+	mpirun -np 4 turbogap md > /dev/null
		+
		+	mv thermo.log relax/relax_$n.log
		+	tail -218 trajectory_out.xyz > relax/relax_$n.xyz
		+
		+	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.

@@ Line 178: / Line 178: @@
   done
-This workflow 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:
+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>