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.

This tutorial is found on the TurboGAP wiki: turbogap.fi (-> tutorials -> Creating oxygenated amorphous carbon) https://turbogap.fi/wiki/index.php/Creating_oxygenated_amorphous_carbon.

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 simulations, to relax the system.

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.
   • Adaptive time-scale MD (by Uttiyoarnab Saha).
3. Prediction of an arbitrary number of local properties.
   • 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 Dispersion capability, with multiple optimisations.
   • Max Veit is developing electrostatics.
4. 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

If you have a CSC account and can ssh into Mahti/Puhti there is no need to install anything. TurboGAP is installed in the path (on Mahti / Puhti)

/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 check the project in creating_a-COx/sample_submit_script.sh, change it from

#SBATCH --account=project_

to

#SBATCH --account=project_2008666

This tutorial depends on ase, so install it as so (after loading the python module)

module load python-data/3.9-22.04
pip install ase --user

Each of the simulations should take ~5-8 minutes.

Installation

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

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

Running this tutorial on a cluster

• Copy the directory turbogap/tutorials/creating_a-COx to where you want to run.
• Change sample_submit_script.sh to reflect the type of job scheduler you use (here, it's slurm), the modules you've loaded for turbogap and python, and change the PATH environment variable to where you've installed turbogap/bin.
• Make sure the project account is correct!
• Load the python module you will use, and make sure ase is installed by running pip install ase {--user}.
• Change the srun turbogap commands in the script_*.sh files to the standard for your cluster (e.g. mpirun -np $N turbogap).
• In each of the directories enumerated with "1.,2., etc", run the scripts in order after the preceding job has finished. They are enumerated with "1.,2., etc" with bash, e.g. bash 1.run_randomise.sh.

Running this tutorial locally

• Copy the directory turbogap/tutorials/creating_a-COx to where you want to run.
• Make sure the python package ase is installed by running pip install ase {--user}.
• Edit the convenience script change_to_local.sh and run with bash.
• Change the path environment variable in sample_submit_script.sh.

Note: how to make a 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_files.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
# Wait for it to finish

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

Note: In create_diamond.py we make 1000 atoms. You can change this to a smaller number, if you want things to run faster.

atoms *= (3,3,3)

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 = bussi             ! Either bussi / berendsen

t_beg = 9000                   ! Initial temperature [K]
t_end = 9000                   ! Final temperature   [K]
tau_t = 100.                   ! Time constant [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 = bussi
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                 ! 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

n_mc_mu = 1
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 see all the options for Monte-Carlo, visit Monte-Carlo.

To run it, we run TurboGAP in mc mode:

srun turbogap mc

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 [math]\displaystyle{ \mu = 0.0 }[/math] eV. You can experiment with different chemical potentials. Using [math]\displaystyle{ \mu = -5.16 }[/math] 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.

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

Hamiltonian Monte-Carlo

Observing the local energies above, we notice that there are some rather high values. To relax them using MC, we can use Hamiltonian MC, which uses the results of NVE molecular dynamics as trial configurations for MC displacements. This gives very high acceptance rates in comparison to other move types due to the energy being approximately conserved from the symplectic (velocity verlet) integrator (hence your acceptance criterion will be approximately 1).

md_nsteps = 20                           ! specifying the number of steps for velocity verlet
md_step = 0.1                            ! 0.1fs timestep

mc_nsteps = 50                           ! Number of mc trial steps to be performed
n_mc_types = 1
mc_types = 'md'
mc_acceptance = 1                        ! Relative rate of choosing the respective trial moves
mc_hamiltonian = .true.                  ! For Hamiltonian Monte-Carlo
                                         !   (NVE ensemble used to increase trial acceptance) for md
                                         !   move type

cd ../3.hamiltonian_mc
bash 1.run_hamiltonian_mc.sh

[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 ../4.volume_mc
bash 1.run_volume.sh

This gives an expansion of the volume.