Simple molecular dynamics

In this tutorial we will run a very simple MD simulation with TurboGAP. We will melt a diamond lattice and then quench the liquid down to 3500K and bake it at that temperature for some time to make graphitic carbon. In this tutorial we will not worry about pressure, but will use the scale_box functionality.

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 a = 3.5
 6 cell = [a, a, a]
 7 positions = [[0, 0, 0],
 8              [a/4, a/4, a/4],
 9              [0, a/2, a/2],
10              [a/4, 3*a/4, 3*a/4],
11              [a/2, 0, a/2],
12              [3*a/4, a/4, 3*a/4],
13              [a/2, a/2, 0],
14              [3*a/4, 3*a/4, a/4]]
15 
16 atoms = Atoms("C8", positions = positions, cell = cell, pbc=True)
17 # If you want a bigger or smaller system, adjust these numbers below
18 atoms *= (3,3,3)
19 
20 # Randomize the velocities at the beginning
21 vel = 0.01 * (np.random.sample([len(atoms),3])-0.5)
22 atoms.set_array("vel", vel)
23 
24 write("diamond.xyz", atoms)

diamond.xyz contains the atomic data.

Melting and expanding

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

wget https://zenodo.org/record/4000211/files/carbon.gap.tar.gz
tar -xvf carbon.gap.tar.gz

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 diamond structure at 9000 K for 3 ps. We choose the Berendsen thermostat with time step 1 fs and time constant 100 fs. At the same time, we are going to expand each of the lattice vectors by a factor of 1.2, since the density of diamond is too high to generate graphitic carbon. You need the following input file:

! Species-specific info
atoms_file = 'diamond.xyz'
pot_file = 'gap_files/carbon.gap'
n_species = 1
species = C
masses = 12.01

! MD options
md_nsteps = 3000
md_step = 1.
thermostat = berendsen
t_beg = 9000
t_end = 9000
tau_t = 100.
write_thermo = 1
write_xyz = 10
scale_box = .true.
box_scaling_factor = 1.2

Then, run TurboGAP. This MD took 6 minutes on my oldish desktop machine, running on 8 cores, so take this into account. You can either run the code serially:

turbogap md

or with parallel support:

mpirun -np 8 turbogap md

Adjust the number of processors and the MPI command as appropriate for your system (e.g., on modern HPC environments, you will often use srun instead of mpirun). For a system of this size you can expect good scaling within a compute node, but not beyond that.

Thermodynamic information will be printed to thermo.log. The first column gives time step, the second gives time, then temperature, kinetic energy, potential energy and pressure. One can get more info with specific keywords in the input file, but we will not cover that here. We can check with gnuplot the evolution of temperature and potential energy:

Temperature and potential energy during the melting of the initial diamond structure.

set multiplot layout 1,2
set ylabel 'Temperature (K)'
set xlabel 'MD time (ps)'
plot "thermo.log" u ($2/1000.):($3) w l not
set ylabel "Potential energy (eV)"
plot "thermo.log" u ($2/1000.):($5) w l not

You can save the info in the trajectory_out.xyz and thermo.log by renaming the file, since those will get overwritten in subsequent steps of this tutorial. You also need to make sure to extract the last snapshot of the melt to use it as starting point for the graphitization. This will do:

1 tail -218 trajectory_out.xyz > melt.xyz

If you have adjusted the number of atoms in your diamond lattice (in the Python script above), then adjust the tail command argument accordingly (use the number of atoms plus 2).

Graphitization

Graphitization is the process whereby a carbon material turns into graphitic carbon at the correct conditions of temperature and pressure. Using MD, we can get graphitic carbon reasonably fast for low pressures and temperatures around 3500 K. "Low" compared to the pressures at which diamond forms, not the pressures at which your ears start to hurt when you are diving. With the GAP we are using here, you will start to see some graphitization around 20 ps into the baking stage, and things are usually quite settled by 100 ps (although it is safest to keep going until the potential energy does not change much, 200 ps should be a good rule of thumb).

Use this input file:

! Species-specific info
atoms_file = 'melt.xyz'
pot_file = 'gap_files/carbon.gap'
n_species = 1
species = C
masses = 12.01

! MD options
md_nsteps = 100000
md_step = 1.
thermostat = berendsen
t_beg = 3500
t_end = 3500
tau_t = 100.
write_thermo = 1
write_xyz = 100

but be mindful about how much time this simulation will take. For my good old desktop, it took about 2h 30min to run all 100 ps of MD on 8 cores. If you just want to get a taste of it, run 20 ps to start seeing graphitic layers forming. If you can run this on a 48 core compute note you will probably get the results quite fast. You can also choose to run the simulation in several steps (e.g., 10 ps of MD each time). In any case, just run TurboGAP as before:

mpirun -np 8 turbogap md

We can check again the evolution of temperature and potential energy:

Additionally, you can choose to quench the structure to, e.g., room temperature. In that case, just extract the last snapshot from the graphitization run, and use it as starting point for another MD run with both t_beg and t_end set to 300, and reduce the number of steps to 10000 or so.

Here you can see the graphitization process and resulting structure: