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 rescale_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 atoms *= (3,3,3)
18 
19 # Randomize the velocities at the beginning
20 vel = 0.01 * (np.random.sample([len(atoms),3])-0.5)
21 atoms.set_array("vel", vel)
22 
23 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'd ofter 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, 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