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$/[math]\displaystyle{ NPT }[/math]), it can perform Monte-Carlo simulations.
   • Full Grand-Canonical Monte Carlo ([math]\displaystyle{ \mu VT }[/math]), with ([math]\displaystyle{ NVT }[/math]) / ([math]\displaystyle{ NPT }[/math]) and more 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-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 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_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).

You can experiment with different chemical potentials. Using mc_mu = -5.16 eV is related to half the binding energy of O₂ at 300K and 1atm. Try it for yourself!

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.

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.

This gives an expansion of the volume.