Creating oxygenated amorphous carbon

This tutorial will focus on using Grand-Canonical Monte-Carlo (GCMC) to determine equilibrium structures and properties 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, using the Hamiltonian MD approach, for an increased acceptance rate.

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.
   • 2b/3b/core potentials (of course).
   • 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 molecular statics (with/without box relaxation) and dynamics (NVT/NPT), it can perform Monte-Carlo simulations.
   • Full Grand-Canonical Monte Carlo (NVT) / (NPT) with multiple move types available:
   • Insertion / removal / displacement / swap / volume / hybrid MD / Hamiltonian.
3. Prediction of an arbitrary number of local properties.
4. ML Van der Waals (by prediction of local hirshfeld volumes) using Tkachenko-Scheffler
   • Heikki Muhli has developed full Many-Body Polarizability capability, with multiple optimisations.
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

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

/projappl/project_2008666/turbogap

To install TurboGAP please run

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

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 nanoporous carbon Structure and Pore Size Distribution in Nanoporous Carbon.

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" (

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

Perform Standard GCMC

Now we perform Grand-Canonical Monte-Carlo to obtain an oxygenated amorphous carbon structure. This is a mu,V,T ensemble, hence we must specify the chemical potential, mu, and the species which we want to insert: here, just oxygen.

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

mc_nsteps = 5000                         ! 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

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_gcmcspecies(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

Monte-Carlo options

Keyword	Definition	Optional	Type	Default	Used when	Example
mc_nsteps	Number of MC steps	N	Int	0	turbogap mc	mc_nsteps = 1000
n_mc_types	Number of MC types	N	Int	0	turbogap mc	n_mc_types = 2
mc_types	Types of MC trials	N	Str(s)	None	n_mc_types > 0	mc_types = 'volume' 'move'
mc_acceptance	Ratios of MC trial moves	N	Int(s)	None	n_mc_types > 0	mc_acceptance = 2 1
n_mc_swaps	Number of swaps	Y	Int	0	'swap' in mc_types	n_mc_swaps = 2
mc_swaps	Swap species	Y	Strs	None	'swap' in mc_types	mc_swaps = 'C' 'O' 'N' 'C'
mc_move_max	Maximum displacement [A]	Y	Float	1.0	'move' in mc_types	mc_move_max = 0.5
mc_lnvol_max	Log volume max for volume moves.	Y	Float	0.01	'volume' in mc_types	mc_lnvol_max = 0.02
mc_min_dist	Minimum distance for insertion [A]	Y	Float	0.2	'insertion in mc_types	mc_min_dist = 0.1
mc_mu	Chemical potential [eV]	Y	Float	0.0	'insertion'/'removal' in mc_types	mc_mu = -5.16
mc_species	GCMC species	Y	Str	None	'insertion'/'removal' in mc_types	mc_species = 'O'
mc_relax	Relax MC trials prior to acc. evaluation	Y	Bool	.false.	turbogap mc	mc_relax = .true.
n_mc_relax_after	Number of specific trials to relax	Y	Int(s)	None	mc_relax = .true.	n_mc_relax_after = 1
mc_relax_after	Relax specific MC trial types	Y	Str(s)	None	n_mc_relax_after > 0	mc_relax_after = 'volume'
mc_nrelax	Number of relaxation steps after trial	Y	Int	0	mc_relax = .true.	mc_nrelax = 50
mc_relax_opt	Optimisation for relaxing after steps	Y	Str	'gd'	mc_relax = .true.	mc_relax_opt = 'gd-box-ortho'
mc_hamiltonian	Use NVE for 'md' trials	Y	Bool	.false.	'md' in mc_types	mc_hamiltonian = .true.
mc_hybrid_opt	Optimisation for 'md' steps	Y	Str	'vv'	'md' in mc_types	mc_hybrid_opt = 'vv'
mc_write_xyz	Debug: read and write to files every step	Y	Bool	.false.	turbogap mc	mc_write_xyz = .true.