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:
- Create an amorphous carbon structure using molecular dynamics, via a melt-quench procedure.
- Perform a standard GCMC calculation to populate the structure with oxygen.
- Perform a hybrid Monte-Carlo/MD simulations, to relax the system.
Contents
Introduction
What is TurboGAP?
TurboGAP is a code used to simulate Machine-Learned Potentials, specifically, Gaussian Approximation Potentials.
It has numerous selling points:
- 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).
- 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).
- 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.
- 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-COxto where you want to run. - Change
sample_submit_script.shto reflect the type of job scheduler you use (here, it's slurm), the modules you've loaded forturbogapand python, and change thePATHenvironment variable to where you've installedturbogap/bin. - Make sure the project account is correct!
- Load the python module you will use, and make sure
aseis installed by runningpip install ase {--user}. - Change the
srun turbogapcommands in thescript_*.shfiles 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-COxto where you want to run. - Make sure the python package
aseis installed by runningpip install ase {--user}. - Edit the convenience script
change_to_local.shand 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:
- We heat up the diamond to 9000K, thereby randomizing the structure.
- We quench to 1000K (actual temp used is 3500K in the real paper).
- We anneal (partially graphitize) the structure at 1000K, to allow the carbon bonds a chance to reform.
- 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:
trajectory_out.xyzis the extended xyz of the trajectory, written every writexyz steps.thermo.loglogs 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 writexyz steps
mc.log: self-explanatory,- format:
mc_step mc_move accepted E_trial E_current N_sites(trial) N_gcmc_species(trial)
- format:
mc_all.xyz: an appended file which contains all accepted movesmc_trial.xyz: a single configuration which contains the trial movemc_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 O2 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 105 - 106 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.
