Graphitization simulation with van der Waals corrections

In this tutorial we are going to carry out simulations of carbon graphitization similar to those in the simple molecular dynamics tutorial, but we will also add a barostat so that the simulation box size changes to mimick the effect of coupling to an external pressure. We will also look at how to use van der Waals (vdW) corrections on top of a GAP force field. Note that vdW corrections are generally not available for all potentials. Those potentials without "proper" vdW corrections (but again, not all) may include fixed pair-wise long-range interactions that allow you to incorporate dispersion effects in a limited way (but we will not cover those in this tutorial).

It is strongly recommended that you do the simple molecular dynamics tutorial before this one.

Prerequisites for this tutorial

Required

• A TurboGAP installation

Optional

• gnuplot (for plotting)
• An ASE installation
• VMD (for visualization, ASE can handle visualization but it's slow)

Getting ready

We will start our simulation from a melted carbon sample with 216 atoms, generated as in the simple molecular dynamics tutorial. For convenience, you can download the file from here (melt.xyz). The approximate density of this sample is 2.15 g/cm³, its temperature is about 9000 K and the pressure is over 200 kbar (although this number will vary depending on the potential used for the calculation). Obviously, this initial configuration is very far away from the graphitic carbon that we want to generate. To graphitize this sample, we will follow the recipe of the previous tutorial, but will use:

1. barostating, to simulate the graphitization process at atmospheric conditions (1 bar);
2. van der Waals corrections, to accurately capture the interlayer interaction in these materials.

The GAP used in the previous tutorial does not incorporate vdW corrections. We will instead use the C60 GAP:

wget https://zenodo.org/record/4616343/files/gap_files.tar.gz
tar -xvf gap_files.tar.gz

Adding a barostat for pressure coupling

In MD simulation, there are different ways to mimick the effect of an external pressure, or "barostating". The most popular of those involves (slowly) dynamically rescaling the size of the simulation box and the atomic positions according to the instantaneous pressure and the target pressure. One of the simplest and most popular barostating strategies in MD is the Berendsen barostat. This is also the only barostat (barostat = berendsen) currently implemented in TurboGAP (we expect more barostats will be implemented as soon as the developers find some time to do it).

The box rescaling can be carried out in different ways, depending on the available degrees of freedom. The simplest rescaling involves increasing or decreasing the length of the lattice vectors isotropically (barostat_sym = isotropic), that is, making the box bigger or smaller without changing its shape. This scaling mode uses only the "average" pressure, computed using the trace of the stress tensor plus the kinetic energy contribution. The next level of sophistication is to allow the x, y and z components of the lattice vectors to change independently, where three independent pressure values are now used, computed from the diagonal components of the stress tensor (barostat_sym = diagonal). This mode can be useful if the elastic properties of the system under study are anisotropic (as in the present graphitic carbon example). The most general box scaling strategy is to allow all the components of the lattice vectors to change dynamically, including coupling to the shear stress. This strategy is not currently implemented in TurboGAP, but it is expected to become available at some point. We also note that changing the box shape arbitrarily means working with a triclinic unit cell, which slows down the neighbor list builds in TurboGAP (which are optimized for orthorhombic unit cells only) and can make MD for large systems inefficient.

A sample input file with a barostat (and without any vdW corrections yet!) looks like this:

! 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
barostat = berendsen
barostat_sym = iso
p_beg = 1.
p_end = 1.
gamma_p = 10.
tau_p = 1000.
write_thermo = 1
write_xyz = 10
write_lv = .true.
neighbors_buffer = 0.5

Run this simulation with turbogap md. It will take a couple of hours with 4 CPU cores (mpirun -np 4 turbogap md).

For a barostat, the target pressure and time constant are not sufficient, one also needs to specify the characteristic elastic response of the material or liquid under study (a stiff material will deform slower than a flexible one under the same amount of pressure). One usually provides the bulk modulus or inverse compressibility of the material. In TurboGAP, one provides this quantity, gamma_p, in units of the inverse compressibility of liquid water. For example, gamma_p = 10. above uses 10 times the inverse compressibility of water. It is important that stiff materials are not deformed too fast since this means a lot of potential energy will be pumped into or removed from the system, leading to potentially catastrophic consequences (blowing up).

After graphitization at high temperature, which lasts for 100 ps, we will need to anneal the system down to room temperature. Starting from the last snapshot of the graphitization, we can make a starting configuration for annealing:

tail -218 trajectory_out.xyz > graphitize.xyz

And then just change a few lines in the input file (most notably temperature and number of time steps):

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

! MD options
md_nsteps = 10000
md_step = 1.
thermostat = berendsen
t_beg = 300
t_end = 300
tau_t = 100.
write_thermo = 1
write_xyz = 100
barostat = berendsen
barostat_sym = iso
p_beg = 1.
p_end = 1.
gamma_p = 10.
tau_p = 1000.
write_thermo = 1
write_xyz = 10
write_lv = .true.
neighbors_buffer = 0.5

And, again, run turbogap md, possibly adding MPI support if you have a multicore machine.

Adding van der Waals corrections

Van der Waals (vdW) corrections in TurboGAP can currently be incorporated within the Tkatchenko-Scheffler (TS) scheme. We are curently working on a many-body dispersion (MBD) implementation which, when ready, will be available for any potential compatible with TS.

To add TS corrections to this calculation, we add the following options to the input file:

vdw_type = ts
vdw_rcut = 15.
vdw_sr = 0.94
vdw_d = 20.
vdw_c6_ref = 27.8 ! eV * Angstrom^6
vdw_r0_ref = 1.9 ! Angstrom
vdw_alpha0_ref = 1.778 ! Angstrom^3

The reference values *_ref are element specific. Some defaults (including C) are implemented in TurboGAP, and you do not need to specify these tags for those atomic species unless you want to change the default values. Refer to the original TS paper for a discussion about these values.

Adding anisotropic barostating

To let the lattice vectors change differently along the three Cartesian dimensions (but using the same target pressure, coupling constant and compressibility factor), change barostat_sym = iso to barostat_sym = diag.

Comparing the different simulations

By enabling/disabling vdW corrections and the use of an isotropic or an anisotropic barsotat we have four possible combinations. By comparing them we can shed some light on the effect of the different simulation parameters on the results, with the big caveat in this tutorial that for disordered materials it is necessary to perform configurational sampling to obtain meaningful statistics. I.e., one needs to run several simulations per combination of parameters, changing the random seed for generating the initial structure (positions, velocities, or both).

As we can see in the graphs, the anisotropic barostating stabilizes the system with and without vdW corrections, but it stabilizes it more in the presence of vdW corrections. This is because the simulation box size can be optimized along the direction of stacking of the graphitic layers. In the absence of vdW interactions, the high temperature during the graphitization step allows for a significant expansion of the simulation box along the layer-stacking direction when anisotropic barostating is enabled.