This repository contains input files for simulating the crystallization of silicon using enhanced sampling methods. Silicon is described using the Stillinger-Weber potential Stillinger and Weber, Phys. Rev. B, v. 31, p. 5262, (1985) . A bias potential is constructed using well tempered metadynamics as a function of the collective variable introduced in this article.
- Objectives and learning outcomes
- Requirements and installation
- Introduction
- Example
- Visualization
- Determine SIGMA for the EnvironmentSimilarity CV
- Assignment
In this tutorial the student will learn to perform enhanced sampling simulations. We will study in particular the crystallization of silicon and learn to calculate free energy differences, the equilibrium temperature between two phases, and obtain atomistic insight into a important phase transformation.
This tutorial uses the molecular dynamics engine LAMMPS patched with the PLUMED 2 enhanced sampling plugin.
LAMMPS and PLUMED can be installed with the following commands:
wget https://github.com/lammps/lammps/archive/refs/tags/stable_23Jun2022_update3.tar.gz
tar -xf stable_23Jun2022_update3.tar.gz
cd lammps-stable_23Jun2022_update3/
make lib-plumed args="-b"
make yes-plumed
make yes-manybody
make yes-molecule
make yes-extra-fix
make yes-extra-dump
make mpi -j$(nproc) # or make serial if you don't have an MPI library
lammpsexe=$(pwd)/lmp_mpi # or lammpsexe=$(pwd)/lmp_serial if you don't have an MPI library
cd ../lib/plumed/plumed-2.8.1
source sourceme.sh
Now you should be able to use $lammpsexe and plumed to execute LAMMPS and PLUMED, respectively. If you close the terminal these commands will no longer be available. In order to recover them run the last three commands again.
A better option for more experienced users would be to include these lines in your ~/.bashrc:
lammpsexe=${path_to_lammps}/src/lmp_mpi # or lmp_serial if you don't have an MPI library
source ${path_to_lammps}/lib/plumed/plumed-2.8.1/sourceme.sh
Note that you should replace ${path_to_lammps} with the appropriate path to the LAMMPS folder.
Finally, you can retrieve all files in this repository using:
git clone https://github.com/PabloPiaggi/Crystallization-of-Silicon
cd Crystallization-of-Silicon
Compiling LAMMPS and PLUMED on a computer cluster can be non trivial. You can find an example here.
If you would like to try other compilation options you can find further information in LAMMPS and PLUMED's manual.
During crystallization the disordered atoms of a liquid spontaneously organized into periodic patterns with long range order. The time and lengthscales involved are often too short to be studied with experiments. In this tutorial we will see how we can study this fascinating process using enhanced sampling molecular dynamics simulations. We take as example the case of silicon that crystallizes in the cubic diamond crystal structure.
Below we provide a (very short) summary of the methods that will be employed.
In well tempered metadynamics a bias potential is constructed as a function of some collective variable s.
The bias potential is constructed as a sum of repulsive Gaussians that discourage frequently visited configurations.
In this way, the simulation explores different regions of the free energy surface
of the system.
In the long time limit the bias potential converges to,
where is the bias factor.
The collective variable (CV) that we will use is based on comparing the atomic environments in the simulation with those of a reference crystal structure.
The environments and
are compared using the kernel,
where is the atomic density around environment
.
In this way we obtain one value of the kernel per atom in the system.
We will then use as collective variable the number of
that are larger than some threshold.
This is equivalent to counting the number of atoms that have a crystalline environment.
We will also calculate the average of the
.
You can find more details about the CV in this article.
The folder Metadynamics-1700K contains the input files to run the simulation.
You can run the example with the command:
cd Metadynamics-1700K
$lammpsexe < start.lmp > out.lmp
This will run a 5 ns long metadynamics simulation at 1700 K and 1 bar. Several output files will be created. out.lmp file contains LAMMPS' output. PLUMED's output files are log.plumed and COLVAR. Inspect the COLVAR file, this will contain the collective variable, the bias potential, and other interesting quantities as a function of simulation time.
While the simulation is running you can check its progress by plotting the collective variable as a function of time. Furthermore, you can calculate the FES with the command
plumed sum_hills --hills HILLS --mintozero
This will create a new file fes.dat. Plot the contents of this file and track the convergence of the bias potential.
This analysis is performed in the Jupyter Notebook Metadynamics-1700K/Results/Analysis.ipynb
Explore this notebook online!
We recommend using the software Ovito for the visualization of the trajectories. Ovito is very user friendly, and you can find instructions and download it from their website.
Once that you have loaded the dump file into Ovito you can color the atoms according to the degree of order around them.
Apply the Identify diamond structure modifier that can be chosen from the Add modification dropdown menu.
Below we show liquid and solid configurations colored with the modifier Identify diamond structure.
One of the parameters that has to be chosen is in the definition of the collective variable.
One way to choose it is by analyzing the distributions of
for the liquid and the solid.
One can then determine the overlap between these distributions and choose the value of
that minimizes the overlap, therefore maximizing the ability of the CV to discriminate between structures.
The threshold of
that determines whether an atom belongs to one phase or the other has to be chosen based on this distribution.
Here we show a plot of the overlap as a function of SIGMA:
More details can be found in the Metadynamics-1700K/Results/Analysis.ipynb Jupyter Notebook.
The following tasks are proposed:
- Repeat the simulation at different temperatures, e.g. in steps of 50 o 100 K
- Study the convergence of the free energy surface (FES) as a function of simulation time (a minimum of 2 ns per simulation is required and 5 ns are suggested)
- Plot the FES obtained at different temperatures
- Calculate the liquid-solid free energy difference (one should integrate the FES but the difference in free energy between the minima should be a good approximation)
- Calculate the melting temperature and discuss finite size effects
- Run simulations to determine the optimal SIGMA parameter for the EnvironmentSimilarity CV