Creating oxygenated amorphous carbon

This tutorial will focus on using molecular dynamics and Grand-Canonical Monte-Carlo (GCMC) simulations to determine equilibrium structures of oxygenated amorphous carbon using TurboGAP.

The structure of this tutorial as as follows:

1. Create an amorphous carbon structure using molecular dynamics, via a melt-quench procedure.
2. Perform a standard GCMC calculation to populate the structure with oxygen.
3. Perform a hybrid Monte-Carlo/MD simulation, to relax the volume.

Introduction

What is TurboGAP?

TurboGAP is a code used to simulate Machine-Learned Potentials, specifically, Gaussian Approximation Potentials.

It has numerous selling points:

1. It is fast.
   • It uses soap turbo descriptors, which are both faster and more accurate than your typical SOAP expansion (also found in QUIP). See the original paper paper by Miguel Caro for more details: Optimizing many-body atomic descriptors for enhanced computational performance of machine learning based interatomic potentials
   • MPI parallelised, (overlapping domain decomposition currently being developed with help of CSC).
   • GPU implementation in progress (with CSC support too).
2. It can perform not just typical molecular statics (with/without box relaxation) and dynamics ([math]\displaystyle{ NVT }[/math] / [math]\displaystyle{ NPT }[/math]), it can perform Grand-Canonical Monte-Carlo simulations.
   • Grand-Canonical Monte Carlo ([math]\displaystyle{ \mu VT }[/math]), with ([math]\displaystyle{ NVT }[/math]) / ([math]\displaystyle{ NPT }[/math]) move types available:
3. ML Van der Waals (by prediction of local hirshfeld volumes) using Tkatchenko-Scheffler Machine learning force fields based on local parametrization of dispersion interactions
   • Heikki Muhli has developed Many-Body Polarizability capability, with multiple optimisations.
4. Prediction of an arbitrary number of local properties.
5. Sneak Peek!: we have added the capability to predict/simulate numerous types of experimental data (ML XPS/XRD) and can allow them to influence simulation. (Talk to Tigany Zarrouk/look out for the papers when they come out on arXiv)!

Installing TurboGAP

For the MLIP workshop

This step is not necessary if you are at the MLIP workshop. TurboGAP is installed in the path (on Mahti)

/projappl/project_2008666/turbogap

The tutorial is in

cd /projappl/project_2008666/turbogap/tutorials/creating_a-COx

Copy the directory creating_a-COx to wherever you want to do the simulations and then change the project in creating_a-COx/sample_submit_script.sh from

#SBATCH --account=project_

to

#SBATCH --account=project_2008666

Installation

To install TurboGAP please run

git clone --recursive http://github.com/mcaroba/turbogap.git /your/turbogap/source/directory

This tutorial in its current form is in the creating_a-COx branch.

git switch creating_a-COx
cd turbogap/tutorials/creating_a-COx

Copy the directory creating_a-COx to wherever you want to do the simulations.

Where /your/turbogap/source/directory is the directory where you're putting the TurboGAP source code. To build the TurboGAP binary and library, you need to select the options that best match your architecture, by editing this line in the Makefile with one of the names of the corresponding makefiles in turbogap/makefiles:

include makefiles/Makefile.Ubuntu_gfortran_mpi

Then just run make

make

Make the potential work in TurboGAP

You must convert potentials which are trained from libAtoms (the xml files) to *.gap files. It can be run by

python3 /path/turbogap/tools/quip_to_xml/make_gap_fit.py  your_potential.xml your_potential.gap {your_hirshfeld.xml}

Tutorial

Create Amorphous Carbon

Here, we perform molecular dynamics simulations to form amorphous carbon from diamond. To do this, we use a simple melt-quench procedure, which is modified from the paper of Wang et al. to create amorphous carbon Structure and Pore Size Distribution in Nanoporous Carbon.

This is also similar to what is done in other tutorials Graphitization simulation with van der Waals corrections Generating amorphous silicon from quenching simulations.

The procedure we will follow is:

1. We heat up the diamond to 9000K, thereby randomizing the structure.
2. We quench to 1000K (actual temp used is 3500K in the real paper).
3. We anneal (partially graphitize) the structure at 1000K, to allow the carbon bonds a chance to reform.
4. Relax the structure, allowing both atomic positions and cell vectors to relax.

To run these calculations, do

cd 1.make_amorphous_carbon
bash 1.run_randomise.sh
# Wait for it to finish

bash 2.run_quench.sh
# Wait for it to finish

bash 3.run_anneal.sh
# Wait for it to finish

bash 4.run_relax.sh

1. Randomise

First, we create a diamond structure using ASE, changing the volume to achieve a given density.

sim_name="md_diamond_randomise"

input_atoms="atoms.xyz"
output_atoms="atoms_randomise.xyz"

ln -sf ../gap_files ./

# 1. Create diamond structure (1000 atoms in atoms.xyz file)
echo "> Running: python create_diamond.py"
python3 create_diamond.py

cp diamond.xyz $input_atoms

Then we create the input file for TurboGAP in input.

! Species-specific info
atoms_file = '${input_atoms}'  ! Input file
pot_file = 'gap_files/CO.gap'  ! path to gap_files
n_species = 2                  ! > Actually the number of species in atoms.xyz is 1 (C),
                               !   but we will add oxygen in future simulations
species = C O
masses = 12.01 15.99

! MD options
md_nsteps = 5000               ! Number of MD steps (5ps randomise - actual time in paper is 20ps)
md_step = 1                    ! MD timestep [fs]
thermostat = berendsen         ! Either berendsen / bussi

t_beg = 9000                   ! Initial temperature [K]
t_end = 9000                   ! Final temperature   [K]
tau_t = 100.                   ! Time constant for berendsen [fs]

! Output
write_thermo = 1               ! Write thermodynamic information every step
                               !       (Step, Time, Temp, Kin_E, Pot_E, Pres)

write_xyz = 200                ! Write extended xyz trajectory (trajectory_out.xyz) every 200 steps

                               ! > Predicted local properties are in the xyz, such as the local energy
                               !   and if specified, hirshfeld volumes, core electron binding energies etc

We run MD/relaxation simulations using the md option

srun turbogap md

This simulation outputs a few files:

1. trajectory_out.xyz is the extended xyz of the trajectory, written every write_xyz steps.
2. thermo.log logs the MD simulation

2. Quench

Here, we cool the system.

The only things necessary to change in the input file (if you rebel against using the provided scripts) is the temperature and the number of MD steps.

! MD options
md_nsteps = 5500               ! 5.5ps quench
md_step = 1.
thermostat = berendsen
t_beg = 9000                   ! Quenching from 9000K
t_end = 1000                   !             to 1000K

3. Anneal

We can anneal the structure to graphitize. The method demonstrated here is different from that of the paper (we are hardly graphitising at 1000K and 10ps), but its here to illustrate the use of a barostat.

The necessary changes to the input file are

! MD options                                                                                                                        
md_nsteps = 10000                 ! Graphitization (actual time in the paper is 200ps)
md_step = 1.
barostat = berendsen              ! Using barostat
t_beg = 1000                                                                                                                        
t_end = 1000
p_beg = 1.0                       ! Initial pressure [bar]
p_end = 1.0                       ! Final pressure   [bar]
tau_t = 100.

4. Relax

Now we relax the structure. We can choose multiple options for this, but here we opt for relaxing the box and the positions.

! MD options
optimize = gd-box-ortho       ! optimize option allows us to specify the type of relaxation
                              ! > Can use "gd" (gradient descent)
                              !           "gd-box-ortho" (gradient descent relaxing diagonal cell vector components)
                              !           "gd-box" (gradient-decent relaxing all cell vector components)

! e_tol = 1.d-6               ! Default energy/force tolerances used
! f_tol = 0.010
md_nsteps = 2000

We can now look at the simulation results by running

cat traj_* > all_traj.xyz

and using your atom viewer of choice.

At the end, you should have a structure which is similar to this:

You can look at the local energies predicted by the GAP (here I use ovito).

• Open structure in ovito
• Use modifier "Create Bonds"
• Use modifier "Color Coding"
  • Change "Input Property" to "Potential Energy" (it might also be called "local energy").

Perform Standard Grand-Canonical Monte-Carlo (GCMC)

Theory of GCMC

In a GCMC simulation, a system of interest is at fixed volume, allowed to thermalize by contact with a heat bath, and it can exchange particles with an infinite reservoir, forming a constant ([math]\displaystyle{ \mu,V,T }[/math]) ensemble.

We perform GCMC using a Markov Chain: starting from an initial pure a-C[math]\displaystyle{ _x }[/math] structure, we generate trial configurations by either randomly displacing a particular atom, or inserting/removing oxygen into/from a random position respectively. These trial configurations are either accepted or rejected using the standard acceptance criteria (see Frenkel 2002) for particle displacement/insertion/removal

[math]\displaystyle{ \mathrm{acc}( \mathrm{move}) = \mathrm{min}\biggl[1, \exp\left\{ -\beta ( E(\mathrm{trial}) - E(\mathrm{current}) ) \right\} \biggr] }[/math]

[math]\displaystyle{ \mathrm{acc}(N \rightarrow N+1) = \mathrm{min}\biggl[ 1, \frac{V}{\lambda^3 (N+1)} \exp\left\{ - \beta ( E(N+1) - E(N) - \mu ) \right\} \biggr] }[/math]

[math]\displaystyle{ \mathrm{acc}(N \rightarrow N-1) = \mathrm{min}\biggl[1, \frac{\lambda^3 N}{V} \exp\left\{ -\beta ( E(N-1) - E(N) + \mu ) \right\} \biggr] }[/math]

where [math]\displaystyle{ \lambda }[/math] is the thermal de-Broglie wavelength which is given by [math]\displaystyle{ \lambda = \sqrt{\frac{2\pi \hslash^2}{mk_BT}} }[/math]. We then repeat the procedure with the last accepted configuration until the the maximum number of iterations has been reached.

GCMC in TurboGAP

Now we perform Grand-Canonical Monte-Carlo to obtain an oxygenated amorphous carbon structure. This is a [math]\displaystyle{ \mu VT }[/math] ensemble, hence we must specify the chemical potential, [math]\displaystyle{ \mu }[/math], and the species which we want to insert: here, just oxygen.

The format is similar to above, but there are more options:

mc_nsteps = 10000                         ! Number of mc trial steps to be performed

n_mc_types = 3                           ! Number of mc trial types
mc_types = 'move' 'insertion' 'removal'  ! MC types can be: 'insertion' 'removal' 'md' 'swap' 'move'
                                         !                  'volume'
mc_acceptance = 1 1 1                    ! Ratios for choosing the respective trial moves (all equally likely here)

mc_move_max = 0.2                        ! Maximum distance for the move of a particular particle

mc_mu = 0.0                              ! gcmc: Chemical potential [eV], using a large one here
mc_species = 'O'                         ! gcmc: species to insert / remove
mc_min_dist = 0.1                        ! gcmc: minimum distance between particles for insertion

To run it, we run TurboGAP in mc mode:

srun turbogap mc

Which can be done by To run these calculations, do

cd ../2.standard_gcmc
bash 1.run_gcmc.sh

The following files are written to every write_xyz steps

• mc.log: self-explanatory,
  • format: mc_step mc_move accepted E_trial E_current N_sites(trial) N_gcmc_species(trial)
• mc_all.xyz: an appended file which contains all accepted moves
• mc_trial.xyz: a single configuration which contains the trial move
• mc_current.xyz: a single configuration which contains the current accepted

You should find a structure which looks like this (the last configuration in mc_all.xyz).

Here, we use a chemical potential of $μ = 0.0$eV. You can experiment with different chemical potentials. Using $μ = -5.16$eV is related to half the binding energy of O₂ at 300K and 1atm. Try it for yourself (in your own time)!

For this short simulation and this size of box, we will not reach convergence of the oxygen content (MC simulation steps are on the order of 10⁵ - 10⁶ steps). We can run

module load python-data/3.9
python analyse_O_content.py

to see the oxygen content.

To see all the options for Monte-Carlo, visit Monte-Carlo.

Here, we can see the local energy for the GCMC configuration.

[math]\displaystyle{ NPT }[/math] using volume moves

We can do constant pressure Monte-Carlo by doing volume moves. In fact, we can mix volume moves and MD [math]\displaystyle{ NPT }[/math] moves to create trial configurations, if we so desire.

Looking at the input file, we see that we've specified a volume type MC move in which the acceptance criterion is given by [math]\displaystyle{ \mathrm{acc}(V \rightarrow V') = \mathrm{min}\biggl[1, \exp\left\{ -\beta ( E(V') - E(V) + P(V-V') - (N+1)\ln (V/V')/\beta ) \right\} \biggr] }[/math]

mc_nsteps = 200

n_mc_types = 3
mc_types = 'move' 'volume' 'md'                  ! We now specify an MD move type for doing NPT

mc_acceptance = 1 1 1
mc_move_max = 0.2

mc_lnvol_max = 0.01                              ! Maximum lnvol to modify the volume

! Specify MD configuration
t_beg = 300
t_end = 300
p_beg = 1.0
p_end = 1.0

md_nsteps = 30
md_step = 0.1                                    ! Reducing the timestep
barostat = berendsen
tau_t = 100.

To run it we do

cd ../3.volume_mc
bash 1.run_volume.sh

This gives an expansion of the volume.