ChemWiki - User contributions [en]

Mod:Hunt Research Group/Charges

2017-06-06T12:18:11Z

Jc6613:

Work in progress.

This article outlines how to obtain NBO, ESP and RESP-fitted charges.

== Overview ==

There are two main classes of charge assignment methods: one by the partitioning of orbitals and one by statistical fitting of the electrostatic potential (ESP).

The former class includes the popular Mulliken population analysis (MPA), Lowdin orbital population analysis, and natural bond orbital analysis (NBO). There are also variations on MPA. Many of these are supported by <nowiki>[http://www.gaussian.com/g_tech/g_ur/k_population.htm Gaussian]</nowiki>. Comparing across these orbital-based methods can lead to conflicting conclusions and there is no general consensus regarding which is the most accurate, although it is widely agreed that MPA is perhaps the least accurate. Nonetheless MPA is calculated by Gaussian by default after every structural optimisation.

The latter relies upon statistical analysis on the true ESP of the molecule, by minimising the root-mean-square deviation of the reproduced ESP from the point charges. RMS error is a way to quantify statistical quality of such assignment but an alternative one, the relative RMS, divides the RMS error by the overall molecular charge and is a better measure of quality for charged species.

===Averaging vs. Equivalencing===

Many bonds rotate faster than the MD timescale and therefore in this context, some atoms can be equivalent to each other. An example of this is the methyl hydrogens, where removing conformational dependence at the methyl produces chemically relevant charges.

Gaussian population analysis does not account for the equivalence of atoms. Charges on equivalent atoms can be ''a priori'' averaged so that these charges make more sense when analysed.

A better method is to introduce equivalence before the charge assignment process, which is done by the Restrained ESP fitting method. RESP also does one additional step of ESP fitting in which only hydrogen atoms are allowed to optimise and the heavier atoms have fixed charges. By this nature, RESP charges do not need *a priori* averaging.

===Sampling Schemes (ESP and RESP only)===

ESP-fitting algorithms requires sampling of the true ESP around the molecule. There are many sampling schemes, the most popular of which are the Merz-Kollman (MK) and the CHelpG schemes. In some cases MK/ESP-fitted charges are more conformation independent. In general, ESP-fitted charges vary less with different sampling schemes than orbital-based charges do upon changing the orbital partitioning algorithm.

Sampled charges can stored as .esp files in the AnteChamber ESP format via Gaussian 09. However Gaussian 09 itself can do some basic ESP-fitting with the `pop=` keyword ( [(http://www.gaussian.com/g_tech/g_ur/k_population.htm) and the atomic charges can be extracted from the log file.

== NBO ==
* After obtaining an optimised structure, run a Gaussian job with the following keywords:
pop=nbo
* It is recommended that you use the special version of Gaussian, g09/b01-nbo on the HPC for this calculation. This gaussian has been compiled with a newer version of NBO embedded.
*If you are using the special version of gaussian, the b01 version does not support the empirical dispersion keyword. You should remove the keyword from your input file.
*NBO charges can be obtained using gauss view, or read from the resulting log file.
*You can average the charges of equivalent atoms

== ESP ==

== RESP ==
[[File:Resp flowchart.png|none|thumb|1103x1103px]]

Mod:Hunt Research Group/Charges

2017-06-06T12:17:56Z

Jc6613: /* NBO */

This article outlines how to obtain NBO, ESP and RESP-fitted charges.

== Overview ==

There are two main classes of charge assignment methods: one by the partitioning of orbitals and one by statistical fitting of the electrostatic potential (ESP).

The former class includes the popular Mulliken population analysis (MPA), Lowdin orbital population analysis, and natural bond orbital analysis (NBO). There are also variations on MPA. Many of these are supported by <nowiki>[http://www.gaussian.com/g_tech/g_ur/k_population.htm Gaussian]</nowiki>. Comparing across these orbital-based methods can lead to conflicting conclusions and there is no general consensus regarding which is the most accurate, although it is widely agreed that MPA is perhaps the least accurate. Nonetheless MPA is calculated by Gaussian by default after every structural optimisation.

The latter relies upon statistical analysis on the true ESP of the molecule, by minimising the root-mean-square deviation of the reproduced ESP from the point charges. RMS error is a way to quantify statistical quality of such assignment but an alternative one, the relative RMS, divides the RMS error by the overall molecular charge and is a better measure of quality for charged species.

===Averaging vs. Equivalencing===

Many bonds rotate faster than the MD timescale and therefore in this context, some atoms can be equivalent to each other. An example of this is the methyl hydrogens, where removing conformational dependence at the methyl produces chemically relevant charges.

Gaussian population analysis does not account for the equivalence of atoms. Charges on equivalent atoms can be ''a priori'' averaged so that these charges make more sense when analysed.

A better method is to introduce equivalence before the charge assignment process, which is done by the Restrained ESP fitting method. RESP also does one additional step of ESP fitting in which only hydrogen atoms are allowed to optimise and the heavier atoms have fixed charges. By this nature, RESP charges do not need *a priori* averaging.

===Sampling Schemes (ESP and RESP only)===

ESP-fitting algorithms requires sampling of the true ESP around the molecule. There are many sampling schemes, the most popular of which are the Merz-Kollman (MK) and the CHelpG schemes. In some cases MK/ESP-fitted charges are more conformation independent. In general, ESP-fitted charges vary less with different sampling schemes than orbital-based charges do upon changing the orbital partitioning algorithm.

Sampled charges can stored as .esp files in the AnteChamber ESP format via Gaussian 09. However Gaussian 09 itself can do some basic ESP-fitting with the `pop=` keyword ( [(http://www.gaussian.com/g_tech/g_ur/k_population.htm) and the atomic charges can be extracted from the log file.

== NBO ==
* After obtaining an optimised structure, run a Gaussian job with the following keywords:
pop=nbo
* It is recommended that you use the special version of Gaussian, g09/b01-nbo on the HPC for this calculation. This gaussian has been compiled with a newer version of NBO embedded.
*If you are using the special version of gaussian, the b01 version does not support the empirical dispersion keyword. You should remove the keyword from your input file.
*NBO charges can be obtained using gauss view, or read from the resulting log file.
*You can average the charges of equivalent atoms

== ESP ==

== RESP ==
[[File:Resp flowchart.png|none|thumb|1103x1103px]]

Mod:Hunt Research Group/Charges

2017-06-06T11:14:13Z

Jc6613: /* NBO */

This article outlines how to obtain NBO, ESP and RESP-fitted charges.

== Overview ==

There are two main classes of charge assignment methods: one by the partitioning of orbitals and one by statistical fitting of the electrostatic potential (ESP).

The former class includes the popular Mulliken population analysis (MPA), Lowdin orbital population analysis, and natural bond orbital analysis (NBO). There are also variations on MPA. Many of these are supported by <nowiki>[http://www.gaussian.com/g_tech/g_ur/k_population.htm Gaussian]</nowiki>. Comparing across these orbital-based methods can lead to conflicting conclusions and there is no general consensus regarding which is the most accurate, although it is widely agreed that MPA is perhaps the least accurate. Nonetheless MPA is calculated by Gaussian by default after every structural optimisation.

The latter relies upon statistical analysis on the true ESP of the molecule, by minimising the root-mean-square deviation of the reproduced ESP from the point charges. RMS error is a way to quantify statistical quality of such assignment but an alternative one, the relative RMS, divides the RMS error by the overall molecular charge and is a better measure of quality for charged species.

===Averaging vs. Equivalencing===

Many bonds rotate faster than the MD timescale and therefore in this context, some atoms can be equivalent to each other. An example of this is the methyl hydrogens, where removing conformational dependence at the methyl produces chemically relevant charges.

Gaussian population analysis does not account for the equivalence of atoms. Charges on equivalent atoms can be *a priori* averaged so that these charges make more sense when analysed.

A better method is to introduce equivalence before the charge assignment process, which is done by the Restrained ESP fitting method. RESP also does one additional step of ESP fitting in which only hydrogen atoms are allowed to optimise and the heavier atoms have fixed charges. By this nature, RESP charges do not need *a priori* averaging.

===Sampling Schemes (ESP and RESP only)===

ESP-fitting algorithms requires sampling of the true ESP around the molecule. There are many sampling schemes, the most popular of which are the Merz-Kollman (MK) and the CHelpG schemes. In some cases MK/ESP-fitted charges are more conformation independent. In general, ESP-fitted charges vary less with different sampling schemes than orbital-based charges do upon changing the orbital partitioning algorithm.

Sampled charges can stored as .esp files in the AnteChamber ESP format via Gaussian 09. However Gaussian 09 itself can do some basic ESP-fitting with the `pop=` keyword ( [(http://www.gaussian.com/g_tech/g_ur/k_population.htm G09 Keyword: Population] ) and the atomic charges can be extracted from the log file.

== NBO ==
* After obtaining an optimised structure, run a Gaussian job with the following keywords: pop=nbo
`pop=nbo`

* It is recommended that you use the special version of Gaussian, g09/b01-nbo on the HPC for this calculation. This gaussian has been compiled with a newer version of NBO embedded.
*

== ESP ==

== RESP ==
[[File:Resp flowchart.png|none|thumb|1103x1103px]]

Mod:Hunt Research Group/Charges

2017-06-06T11:14:03Z

Mod:Hunt Research Group

2017-06-06T10:42:23Z

Jc6613: /* Research Notes */

==Hunt Group Wiki==

Back to the main [http://www.ch.ic.ac.uk/hunt web-page]
===HPC Resources===
#Computing resources available in the chemistry department [https://www.ch.ic.ac.uk/wiki/index.php/Mod:Hunt_Research_Group/computing_resources link]
#HPC servers and run scripts [https://www.ch.ic.ac.uk/wiki/index.php/Mod:Hunt_Research_Group/hpc link]
#Setting up a connection to HPC if you have a PC [https://www.ch.ic.ac.uk/wiki/index.php/Mod:Hunt_Research_Group/hpc_connections link]
#How to fix Windows files under UNIX [https://www.ch.ic.ac.uk/wiki/index.php/Mod:Hunt_Research_Group/Windowsfiles link]
#How to make ssh more comfortable [https://www.ch.ic.ac.uk/wiki/index.php/Mod:Hunt_Research_Group/pimpSSH link]
#How to make qsub more comfortable (gfunc) [https://www.ch.ic.ac.uk/wiki/index.php/Mod:Hunt_Research_Group/pimpQSUB link]
#How to set up a SSH keypair [https://www.ch.ic.ac.uk/wiki/index.php/Mod:Hunt_Research_Group/SSHkeyfile link]
#How to use gaussview directly on the HPC [https://www.ch.ic.ac.uk/wiki/index.php/Mod:Hunt_Research_Group/gview link]
#How to comfortably search through old BASH commands [https://www.ch.ic.ac.uk/wiki/index.php/Mod:Hunt_Research_Group/searchbash link]
#How to connect to HPC directory on desktop for file transfers - MacFusion [https://www.ch.ic.ac.uk/wiki/index.php/Mod:Hunt_Research_Group/hpc_Directory_on_desktop link]
#How to set up cx2 [https://wiki.ch.ic.ac.uk/wiki/index.php?title=Talk:Mod:Hunt_Research_Group/cx2 link]

===Visualisation===
*'''density based visualisation'''
#download [http://aim.tkgristmill.com AIMALL]

*'''ESPs'''
#Instructions for visualizing electrostatic potentials (Gaussview)[https://www.ch.ic.ac.uk/wiki/index.php/Mod:Hunt_Research_Group/electrostatic_potentials link]
#Electrostatic Potentials II (Molden) [https://www.ch.ic.ac.uk/wiki/index.php/Mod:Hunt_Research_Group/electrostatic_potentials_2 link]
*'''VMD: a molecular dynamics visualisation package'''
#Download VMD [//wiki.ch.ic.ac.uk/wiki/index.php?title=Talk:Mod:Hunt_Research_Group/download_vmd link]
#Quick reminder [https://www.ch.ic.ac.uk/wiki/index.php/Talk:Mod:Hunt_Research_Group/VMDReminder link]
#Tricks and tips [https://www.ch.ic.ac.uk/wiki/index.php/Talk:Mod:Hunt_Research_Group/VMDTips link]
#Changing the graphical representation of your structures [https://www.ch.ic.ac.uk/wiki/index.php/Talk:Mod:Hunt_Research_Group/vmd link]
#Basic visualisation of a trajectory [https://www.ch.ic.ac.uk/wiki/index.php/Talk:Mod:Hunt_Research_Group/VisualisingyourSimulation link]
#How to turn a Gaussian optimization into a VMD movie [https://www.ch.ic.ac.uk/wiki/index.php/Mod:Hunt_Research_Group/VMDmovie link]
#Using scripts in VMD [https://wiki.ch.ic.ac.uk/wiki/index.php?title=Talk:Mod:Hunt_Research_Group/VmdScripts link]
#Dealing with periodic boundaries and bonding (under construction) [https://wiki.ch.ic.ac.uk/wiki/index.php?title=Talk:Mod:Hunt_Research_Group/VmdScriptsPeriodic link]
#Dealing with bonding (under construction) [https://wiki.ch.ic.ac.uk/wiki/index.php?title=Talk:Mod:Hunt_Research_Group/VmdBonding link]
#Overlapping two structures [https://wiki.ch.ic.ac.uk/wiki/index.php?title=Talk:Mod:Hunt_Research_Group/VmdVisual link]
*'''JMol'''
#Visualising MOs using Jmol [https://wiki.ch.ic.ac.uk/wiki/index.php?title=Mod:basic_jmol_instructions link]
#Surfaces (Solvent-Accessible and Connolly) in Jmol [https://www.ch.ic.ac.uk/wiki/index.php/Mod:Hunt_Research_Group/jmolsurfaces link]
*'''EMO Code'''
#How to use Ling's emo plot code[https://wiki.ch.ic.ac.uk/wiki/index.php?title=Talk:Mod:Hunt_Research_Group/emoplot link]
#How to plot EMOs [https://wiki.ch.ic.ac.uk/wiki/index.php?title=Talk:Mod:Hunt_Research_Group/emo link]
*'''PyGauss'''
#Python API for analysis of Gaussian compuations [https://pygauss.readthedocs.org - Documentation]
*'''MD Post processing'''
#Code to Recentre DL_PLOY HISTORY file [https://wiki.ch.ic.ac.uk/wiki/index.php?title=Talk:Mod:Hunt_Research_Group/recentre_xyz.py link]
#Link to the code to convert the DL_POLY HISTORY file to the multi-frame XYZ file[https://wiki.ch.ic.ac.uk/wiki/index.php?title=Talk:Mod:Hunt_Research_Group/his2xyz.py link]

*'''SDFs'''
#How to generate SDFs [https://wiki.ch.ic.ac.uk/wiki/index.php?title=Talk:Mod:Hunt_Research_Group/sdfs_generate link]

===Gaussian General===
#We are starting a database of common errors encountered when running Gaussian jobs [https://www.ch.ic.ac.uk/wiki/index.php/Mod:Hunt_Research_Group/gaussian_errors link]
# Here is an already existing database of common errors [https://www.ace-net.ca/wiki/Gaussian_Error_Messages link]
# [http://www.ch.ic.ac.uk/hunt/g03_man/index.htm G03 Manual]
#How to run NBO5.9 on the HPC [https://www.ch.ic.ac.uk/wiki/index.php/Mod:Hunt_Research_Group/NBO5.9 link]
#How to include dispersion [https://www.ch.ic.ac.uk/wiki/index.php/Mod:Hunt_Research_Group/dispersion link]
#Basic ONIOM (Mechanical Embedding) [https://www.ch.ic.ac.uk/wiki/index.php/Mod:Hunt_Research_Group/basiconiom link]
#M0n and DFT-D [https://www.ch.ic.ac.uk/wiki/index.php/Mod:Hunt_Research_Group/DFTD link]
#IL ONIOM clusters [https://www.ch.ic.ac.uk/wiki/index.php/Mod:Hunt_Research_Group/oniomclusers link]
#Molecular volume calculations [https://wiki.ch.ic.ac.uk/wiki/index.php?title=Mod:Hunt_Research_Group/molecular_volume link]
#problems with scf convergence [https://www.ch.ic.ac.uk/wiki/index.php/Mod:Hunt_Research_Group/scf_convergence link]
#Using a z-matrix [https://www.ch.ic.ac.uk/wiki/index.php/Mod:Hunt_Research_Group/z-matrix link]
#generating natural transition orbitals [https://www.ch.ic.ac.uk/wiki/index.php/Mod:Hunt_Research_Group/nto link]
#Using solvent models [https://www.ch.ic.ac.uk/wiki/index.php/Mod:Hunt_Research_Group/solvent link]
#Using SMD on ILs [https://wiki.ch.ic.ac.uk/wiki/index.php?title=Mod:Hunt_Research_Group:_Using_SMD_on_ILs link]
#computing excited state polarisabilities [https://wiki.ch.ic.ac.uk/wiki/index.php?title=Mod:Hunt_Research_Group:_ES_alpha link]
#computing deuterated and/or anharmonic spectra [https://wiki.ch.ic.ac.uk/wiki/index.php?title=Mod:Hunt_Research_Group:_Danharm link]
#manipulating checkpoint files [https://wiki.ch.ic.ac.uk/wiki/index.php?title=Mod:Hunt_Research_Group:usingchkfiles link]
#AimAll with pseudo potentials [https://wiki.ch.ic.ac.uk/wiki/index.php?title=Mod:Hunt_Research_Group:aim_pseudopotentials link]
#Script to pull thermodynamic data and low frequencies from log files [https://wiki.ch.ic.ac.uk/wiki/index.php?title=Mod:Hunt_Research_Group:freq_script link]
#General procedure for locating transition state structures [[link]]

===ADF General===
#Submission script [http://www.ch.ic.ac.uk/wiki/index.php/Mod:Hunt_Research_Group/ADF_sricpt link]

===Codes to Help Analysis===
# Extract E2 Values (From NBO Calculations) [https://wiki.ch.ic.ac.uk/wiki/index.php?title=Talk:Mod:Hunt_Research_Group/NBO_Matlab_Code link]
# Calculate pDoS/XP spectra code (under construction) [https://wiki.ch.ic.ac.uk/wiki/index.php?title=Talk:Mod:Hunt_Research_Group/Calc_XPS_Code link]
# Codes to extract frequency data from gaussian .log files and generate vibrational spectra [https://wiki.ch.ic.ac.uk/wiki/index.php?title=Mod:Hunt_Research_Group:frequency_spectrum_script link]
# Optimally Tuned Range Seperated Hybrid Functionals [https://wiki.ch.ic.ac.uk/wiki/index.php?title=Mod:Hunt_Research_Group/OTRSH_Funct link]

===Setup and Running Ab-Initio MD Simulations===
#CPMD: Car-Parrinello Molecular Dynamics [https://www.ch.ic.ac.uk/wiki/index.php/Talk:Mod:Hunt_Research_Group/cpmd link]
#How to run CPMD to study aqueous solutions [https://www.ch.ic.ac.uk/wiki/index.php/Talk:Mod:Hunt_Research_Group/cpmd_water link]
#How to run CP2K [https://wiki.ch.ic.ac.uk/wiki/index.php?title=Talk:Mod:Hunt_Research_Group/cp2k_how link]
#[bmim]Cl using CPMD [https://wiki.ch.ic.ac.uk/wiki/index.php?title=Talk:Mod:Hunt_Research_Group/bmimCl_cpmd link]
#[bmim]Cl using CP2K [https://wiki.ch.ic.ac.uk/wiki/index.php?title=Talk:Mod:Hunt_Research_Group/bmimCl_cp2k link]
#mman using CPMD and Gaussian [https://wiki.ch.ic.ac.uk/wiki/index.php?title=Talk:Mod:Hunt_Research_Group/mman link]
#[emim]SCN using CP2K[https://wiki.ch.ic.ac.uk/wiki/index.php?title=Talk:Mod:Hunt_Research_Group/emimscn link]
#CP2K Donts [https://wiki.ch.ic.ac.uk/wiki/index.php?title=Talk:Mod:Hunt_Research_Group/cp2k link]

===Setup and Running Classical MD Simulations===
#DLPOLY a MD simulation package, Installation on an IMac (old needs to be updated) [https://www.ch.ic.ac.uk/wiki/index.php/Mod:Hunt_Research_Group/dlpoly_install link]
#DL_POLY FAQs [http://www.stfc.ac.uk/cse/DL_POLY/ccp1gui/38621.aspx] from DL_POLY webpage.
#Installing Packmol
#Getting started: generating a solvated structure and "relaxing" it [https://www.ch.ic.ac.uk/wiki/index.php/Talk:Mod:Hunt_Research_Group/Starting_MD link]
#Equilibration and production simulations [https://www.ch.ic.ac.uk/wiki/index.php/Talk:Mod:Hunt_Research_Group/EquilibrationandProduction link]
#How to equilibrate an MD run[https://www.ch.ic.ac.uk/wiki/index.php/Mod:Hunt_Research_Group/equilibration link]
#Getting the Force Field [https://www.ch.ic.ac.uk/wiki/index.php/Talk:Mod:Hunt_Research_Group/Wheretostart link]
#Choosing an Ensemble [https://www.ch.ic.ac.uk/wiki/index.php/Talk:Mod:Hunt_Research_Group/Ensembles link]
#Molten Salt Simulations [https://wiki.ch.ic.ac.uk/wiki/index.php?title=Mod:Hunt_Research_Group/MoltenSaltSimulation link]
#Common Errors [https://wiki.ch.ic.ac.uk/wiki/index.php?title=Mod:Hunt_Research_Group/CommonErrors link]
#Voids in ILs[https://www.ch.ic.ac.uk/wiki/index.php/Talk:Mod:Hunt_Research_Group/voids link]
#Equilibration of [bmim][BF4] and [bmim][NO3][https://www.ch.ic.ac.uk/wiki/index.php/Mod:Hunt_Research_Group/BmimBF4_equilibration link]
#Summary of discussions with Ruth[https://www.ch.ic.ac.uk/wiki/index.php/Mod:Hunt_Research_Group/Aug09QtoRuth link]

===Running QM/MM Simulations in ChemShell===
#ChemShell official website which contains the manual and a tutorial [http://www.stfc.ac.uk/CSE/randd/ccg/36254.aspx link]
#Introduction to ChemShell - Copper in water [https://wiki.ch.ic.ac.uk/wiki/index.php?title=Talk:Mod:Hunt_Research_Group/ChemShell_Introduction link]
#Defining the system: Cu2+ and its first 2 solvation shells [https://wiki.ch.ic.ac.uk/wiki/index.php?title=Talk:Mod:Hunt_Research_Group/ChemShell_System_Aqeuous_Cu(II) link]
#Defining the force field parameters [https://wiki.ch.ic.ac.uk/wiki/index.php?title=Talk:Mod:Hunt_Research_Group/ChemShell_Force_Field_Parameters_Aqueous_Cu(II) link]
#Single point QM/MM energy calculation [https://wiki.ch.ic.ac.uk/wiki/index.php?title=Talk:Mod:Hunt_Research_Group/QMMM_SP_Aqeuous_Cu(II) link]
#QM/MM Optimisation [https://wiki.ch.ic.ac.uk/wiki/index.php?title=Talk:Mod:Hunt_Research_Group/QMMM_OPT_Aqeuous_Cu(II) link]
#QM/MM Molecular Dynamics [https://wiki.ch.ic.ac.uk/wiki/index.php?title=Talk:Mod:Hunt_Research_Group/QMMM_MD_Aqeuous_Cu(II) link]
#Using MolCluster [https://wiki.ch.ic.ac.uk/wiki/index.php?title=Talk:Mod:Hunt_Research_Group/MolCluster link]
#Running ChemShell [https://wiki.ch.ic.ac.uk/wiki/index.php?title=Talk:Mod:Hunt_Research_Group/ChemShell link]

===Research Notes===
#Cl- in water [https://www.ch.ic.ac.uk/wiki/index.php/Talk:Mod:Hunt_Research_Group/wannier_centre link]
#The use of Legendre time correlation functions to study reorientational dynamics in liquids[https://www.ch.ic.ac.uk/wiki/index.php/Talk:Mod:Hunt_Research_Group/legendre link]
#Functional for ILs using CPMD [https://wiki.ch.ic.ac.uk/wiki/index.php?title=Talk:Mod:Hunt_Research_Group/IL_cpmd_functional link]
#Solving the angular part of the Schrödinger equation for a hydrogen atom [https://wiki.ch.ic.ac.uk/wiki/index.php?title=Talk:Mod:Hunt_Research_Group/angular_schrodinger link] (notes by Vincent)
#Systematic conformational scan for ion-pair dimers [https://wiki.ch.ic.ac.uk/wiki/index.php?title=Talk:Mod:Hunt_Research_Group/ion_pair_scan link]
#Obtaining NBO, ESP, and RESP charges [https://wiki.ch.ic.ac.uk/wiki/index.php?title=Mod:Hunt_Research_Group/Charges link]

===Installing and using other packages===
#How to install POLYRATE [https://www.ch.ic.ac.uk/wiki/index.php/Mod:Hunt_Research_Group/polyrate link]
#How to install Geomview [https://www.ch.ic.ac.uk/wiki/index.php/Mod:Hunt_Research_Group/geomview link]
#XMGRACE, gfortran, c compilers for Lion [http://hpc.sourceforge.net/]

===Admin Stuff===
#Not used to writing a wiki, make your test runs [https://www.ch.ic.ac.uk/wiki/index.php/Mod:Hunt_Research_Group/testing on this page]
#How to set-up new macs [https://www.ch.ic.ac.uk/wiki/index.php/Mod:Hunt_Research_Group/mac_setup link]
#How to set-up remote desktop [https://www.ch.ic.ac.uk/wiki/index.php/Mod:Hunt_Research_Group/mac_remote link]
#[https://www.ch.ic.ac.uk/wiki/index.php/Mod:Hunt_Research_Group/calendar Calendar]
#How to switch the printer HP CP3525dn duplex on and off [https://www.ch.ic.ac.uk/wiki/index.php/Mod:Hunt_Research_Group/printing link]

Mod:Hunt Research Group

2017-06-06T10:42:02Z

Jc6613: /* Research Notes */

==Hunt Group Wiki==

Back to the main [http://www.ch.ic.ac.uk/hunt web-page]
===HPC Resources===
#Computing resources available in the chemistry department [https://www.ch.ic.ac.uk/wiki/index.php/Mod:Hunt_Research_Group/computing_resources link]
#HPC servers and run scripts [https://www.ch.ic.ac.uk/wiki/index.php/Mod:Hunt_Research_Group/hpc link]
#Setting up a connection to HPC if you have a PC [https://www.ch.ic.ac.uk/wiki/index.php/Mod:Hunt_Research_Group/hpc_connections link]
#How to fix Windows files under UNIX [https://www.ch.ic.ac.uk/wiki/index.php/Mod:Hunt_Research_Group/Windowsfiles link]
#How to make ssh more comfortable [https://www.ch.ic.ac.uk/wiki/index.php/Mod:Hunt_Research_Group/pimpSSH link]
#How to make qsub more comfortable (gfunc) [https://www.ch.ic.ac.uk/wiki/index.php/Mod:Hunt_Research_Group/pimpQSUB link]
#How to set up a SSH keypair [https://www.ch.ic.ac.uk/wiki/index.php/Mod:Hunt_Research_Group/SSHkeyfile link]
#How to use gaussview directly on the HPC [https://www.ch.ic.ac.uk/wiki/index.php/Mod:Hunt_Research_Group/gview link]
#How to comfortably search through old BASH commands [https://www.ch.ic.ac.uk/wiki/index.php/Mod:Hunt_Research_Group/searchbash link]
#How to connect to HPC directory on desktop for file transfers - MacFusion [https://www.ch.ic.ac.uk/wiki/index.php/Mod:Hunt_Research_Group/hpc_Directory_on_desktop link]
#How to set up cx2 [https://wiki.ch.ic.ac.uk/wiki/index.php?title=Talk:Mod:Hunt_Research_Group/cx2 link]

===Visualisation===
*'''density based visualisation'''
#download [http://aim.tkgristmill.com AIMALL]

*'''ESPs'''
#Instructions for visualizing electrostatic potentials (Gaussview)[https://www.ch.ic.ac.uk/wiki/index.php/Mod:Hunt_Research_Group/electrostatic_potentials link]
#Electrostatic Potentials II (Molden) [https://www.ch.ic.ac.uk/wiki/index.php/Mod:Hunt_Research_Group/electrostatic_potentials_2 link]
*'''VMD: a molecular dynamics visualisation package'''
#Download VMD [//wiki.ch.ic.ac.uk/wiki/index.php?title=Talk:Mod:Hunt_Research_Group/download_vmd link]
#Quick reminder [https://www.ch.ic.ac.uk/wiki/index.php/Talk:Mod:Hunt_Research_Group/VMDReminder link]
#Tricks and tips [https://www.ch.ic.ac.uk/wiki/index.php/Talk:Mod:Hunt_Research_Group/VMDTips link]
#Changing the graphical representation of your structures [https://www.ch.ic.ac.uk/wiki/index.php/Talk:Mod:Hunt_Research_Group/vmd link]
#Basic visualisation of a trajectory [https://www.ch.ic.ac.uk/wiki/index.php/Talk:Mod:Hunt_Research_Group/VisualisingyourSimulation link]
#How to turn a Gaussian optimization into a VMD movie [https://www.ch.ic.ac.uk/wiki/index.php/Mod:Hunt_Research_Group/VMDmovie link]
#Using scripts in VMD [https://wiki.ch.ic.ac.uk/wiki/index.php?title=Talk:Mod:Hunt_Research_Group/VmdScripts link]
#Dealing with periodic boundaries and bonding (under construction) [https://wiki.ch.ic.ac.uk/wiki/index.php?title=Talk:Mod:Hunt_Research_Group/VmdScriptsPeriodic link]
#Dealing with bonding (under construction) [https://wiki.ch.ic.ac.uk/wiki/index.php?title=Talk:Mod:Hunt_Research_Group/VmdBonding link]
#Overlapping two structures [https://wiki.ch.ic.ac.uk/wiki/index.php?title=Talk:Mod:Hunt_Research_Group/VmdVisual link]
*'''JMol'''
#Visualising MOs using Jmol [https://wiki.ch.ic.ac.uk/wiki/index.php?title=Mod:basic_jmol_instructions link]
#Surfaces (Solvent-Accessible and Connolly) in Jmol [https://www.ch.ic.ac.uk/wiki/index.php/Mod:Hunt_Research_Group/jmolsurfaces link]
*'''EMO Code'''
#How to use Ling's emo plot code[https://wiki.ch.ic.ac.uk/wiki/index.php?title=Talk:Mod:Hunt_Research_Group/emoplot link]
#How to plot EMOs [https://wiki.ch.ic.ac.uk/wiki/index.php?title=Talk:Mod:Hunt_Research_Group/emo link]
*'''PyGauss'''
#Python API for analysis of Gaussian compuations [https://pygauss.readthedocs.org - Documentation]
*'''MD Post processing'''
#Code to Recentre DL_PLOY HISTORY file [https://wiki.ch.ic.ac.uk/wiki/index.php?title=Talk:Mod:Hunt_Research_Group/recentre_xyz.py link]
#Link to the code to convert the DL_POLY HISTORY file to the multi-frame XYZ file[https://wiki.ch.ic.ac.uk/wiki/index.php?title=Talk:Mod:Hunt_Research_Group/his2xyz.py link]

*'''SDFs'''
#How to generate SDFs [https://wiki.ch.ic.ac.uk/wiki/index.php?title=Talk:Mod:Hunt_Research_Group/sdfs_generate link]

===Gaussian General===
#We are starting a database of common errors encountered when running Gaussian jobs [https://www.ch.ic.ac.uk/wiki/index.php/Mod:Hunt_Research_Group/gaussian_errors link]
# Here is an already existing database of common errors [https://www.ace-net.ca/wiki/Gaussian_Error_Messages link]
# [http://www.ch.ic.ac.uk/hunt/g03_man/index.htm G03 Manual]
#How to run NBO5.9 on the HPC [https://www.ch.ic.ac.uk/wiki/index.php/Mod:Hunt_Research_Group/NBO5.9 link]
#How to include dispersion [https://www.ch.ic.ac.uk/wiki/index.php/Mod:Hunt_Research_Group/dispersion link]
#Basic ONIOM (Mechanical Embedding) [https://www.ch.ic.ac.uk/wiki/index.php/Mod:Hunt_Research_Group/basiconiom link]
#M0n and DFT-D [https://www.ch.ic.ac.uk/wiki/index.php/Mod:Hunt_Research_Group/DFTD link]
#IL ONIOM clusters [https://www.ch.ic.ac.uk/wiki/index.php/Mod:Hunt_Research_Group/oniomclusers link]
#Molecular volume calculations [https://wiki.ch.ic.ac.uk/wiki/index.php?title=Mod:Hunt_Research_Group/molecular_volume link]
#problems with scf convergence [https://www.ch.ic.ac.uk/wiki/index.php/Mod:Hunt_Research_Group/scf_convergence link]
#Using a z-matrix [https://www.ch.ic.ac.uk/wiki/index.php/Mod:Hunt_Research_Group/z-matrix link]
#generating natural transition orbitals [https://www.ch.ic.ac.uk/wiki/index.php/Mod:Hunt_Research_Group/nto link]
#Using solvent models [https://www.ch.ic.ac.uk/wiki/index.php/Mod:Hunt_Research_Group/solvent link]
#Using SMD on ILs [https://wiki.ch.ic.ac.uk/wiki/index.php?title=Mod:Hunt_Research_Group:_Using_SMD_on_ILs link]
#computing excited state polarisabilities [https://wiki.ch.ic.ac.uk/wiki/index.php?title=Mod:Hunt_Research_Group:_ES_alpha link]
#computing deuterated and/or anharmonic spectra [https://wiki.ch.ic.ac.uk/wiki/index.php?title=Mod:Hunt_Research_Group:_Danharm link]
#manipulating checkpoint files [https://wiki.ch.ic.ac.uk/wiki/index.php?title=Mod:Hunt_Research_Group:usingchkfiles link]
#AimAll with pseudo potentials [https://wiki.ch.ic.ac.uk/wiki/index.php?title=Mod:Hunt_Research_Group:aim_pseudopotentials link]
#Script to pull thermodynamic data and low frequencies from log files [https://wiki.ch.ic.ac.uk/wiki/index.php?title=Mod:Hunt_Research_Group:freq_script link]
#General procedure for locating transition state structures [[link]]

===ADF General===
#Submission script [http://www.ch.ic.ac.uk/wiki/index.php/Mod:Hunt_Research_Group/ADF_sricpt link]

===Codes to Help Analysis===
# Extract E2 Values (From NBO Calculations) [https://wiki.ch.ic.ac.uk/wiki/index.php?title=Talk:Mod:Hunt_Research_Group/NBO_Matlab_Code link]
# Calculate pDoS/XP spectra code (under construction) [https://wiki.ch.ic.ac.uk/wiki/index.php?title=Talk:Mod:Hunt_Research_Group/Calc_XPS_Code link]
# Codes to extract frequency data from gaussian .log files and generate vibrational spectra [https://wiki.ch.ic.ac.uk/wiki/index.php?title=Mod:Hunt_Research_Group:frequency_spectrum_script link]
# Optimally Tuned Range Seperated Hybrid Functionals [https://wiki.ch.ic.ac.uk/wiki/index.php?title=Mod:Hunt_Research_Group/OTRSH_Funct link]

===Setup and Running Ab-Initio MD Simulations===
#CPMD: Car-Parrinello Molecular Dynamics [https://www.ch.ic.ac.uk/wiki/index.php/Talk:Mod:Hunt_Research_Group/cpmd link]
#How to run CPMD to study aqueous solutions [https://www.ch.ic.ac.uk/wiki/index.php/Talk:Mod:Hunt_Research_Group/cpmd_water link]
#How to run CP2K [https://wiki.ch.ic.ac.uk/wiki/index.php?title=Talk:Mod:Hunt_Research_Group/cp2k_how link]
#[bmim]Cl using CPMD [https://wiki.ch.ic.ac.uk/wiki/index.php?title=Talk:Mod:Hunt_Research_Group/bmimCl_cpmd link]
#[bmim]Cl using CP2K [https://wiki.ch.ic.ac.uk/wiki/index.php?title=Talk:Mod:Hunt_Research_Group/bmimCl_cp2k link]
#mman using CPMD and Gaussian [https://wiki.ch.ic.ac.uk/wiki/index.php?title=Talk:Mod:Hunt_Research_Group/mman link]
#[emim]SCN using CP2K[https://wiki.ch.ic.ac.uk/wiki/index.php?title=Talk:Mod:Hunt_Research_Group/emimscn link]
#CP2K Donts [https://wiki.ch.ic.ac.uk/wiki/index.php?title=Talk:Mod:Hunt_Research_Group/cp2k link]

===Setup and Running Classical MD Simulations===
#DLPOLY a MD simulation package, Installation on an IMac (old needs to be updated) [https://www.ch.ic.ac.uk/wiki/index.php/Mod:Hunt_Research_Group/dlpoly_install link]
#DL_POLY FAQs [http://www.stfc.ac.uk/cse/DL_POLY/ccp1gui/38621.aspx] from DL_POLY webpage.
#Installing Packmol
#Getting started: generating a solvated structure and "relaxing" it [https://www.ch.ic.ac.uk/wiki/index.php/Talk:Mod:Hunt_Research_Group/Starting_MD link]
#Equilibration and production simulations [https://www.ch.ic.ac.uk/wiki/index.php/Talk:Mod:Hunt_Research_Group/EquilibrationandProduction link]
#How to equilibrate an MD run[https://www.ch.ic.ac.uk/wiki/index.php/Mod:Hunt_Research_Group/equilibration link]
#Getting the Force Field [https://www.ch.ic.ac.uk/wiki/index.php/Talk:Mod:Hunt_Research_Group/Wheretostart link]
#Choosing an Ensemble [https://www.ch.ic.ac.uk/wiki/index.php/Talk:Mod:Hunt_Research_Group/Ensembles link]
#Molten Salt Simulations [https://wiki.ch.ic.ac.uk/wiki/index.php?title=Mod:Hunt_Research_Group/MoltenSaltSimulation link]
#Common Errors [https://wiki.ch.ic.ac.uk/wiki/index.php?title=Mod:Hunt_Research_Group/CommonErrors link]
#Voids in ILs[https://www.ch.ic.ac.uk/wiki/index.php/Talk:Mod:Hunt_Research_Group/voids link]
#Equilibration of [bmim][BF4] and [bmim][NO3][https://www.ch.ic.ac.uk/wiki/index.php/Mod:Hunt_Research_Group/BmimBF4_equilibration link]
#Summary of discussions with Ruth[https://www.ch.ic.ac.uk/wiki/index.php/Mod:Hunt_Research_Group/Aug09QtoRuth link]

===Running QM/MM Simulations in ChemShell===
#ChemShell official website which contains the manual and a tutorial [http://www.stfc.ac.uk/CSE/randd/ccg/36254.aspx link]
#Introduction to ChemShell - Copper in water [https://wiki.ch.ic.ac.uk/wiki/index.php?title=Talk:Mod:Hunt_Research_Group/ChemShell_Introduction link]
#Defining the system: Cu2+ and its first 2 solvation shells [https://wiki.ch.ic.ac.uk/wiki/index.php?title=Talk:Mod:Hunt_Research_Group/ChemShell_System_Aqeuous_Cu(II) link]
#Defining the force field parameters [https://wiki.ch.ic.ac.uk/wiki/index.php?title=Talk:Mod:Hunt_Research_Group/ChemShell_Force_Field_Parameters_Aqueous_Cu(II) link]
#Single point QM/MM energy calculation [https://wiki.ch.ic.ac.uk/wiki/index.php?title=Talk:Mod:Hunt_Research_Group/QMMM_SP_Aqeuous_Cu(II) link]
#QM/MM Optimisation [https://wiki.ch.ic.ac.uk/wiki/index.php?title=Talk:Mod:Hunt_Research_Group/QMMM_OPT_Aqeuous_Cu(II) link]
#QM/MM Molecular Dynamics [https://wiki.ch.ic.ac.uk/wiki/index.php?title=Talk:Mod:Hunt_Research_Group/QMMM_MD_Aqeuous_Cu(II) link]
#Using MolCluster [https://wiki.ch.ic.ac.uk/wiki/index.php?title=Talk:Mod:Hunt_Research_Group/MolCluster link]
#Running ChemShell [https://wiki.ch.ic.ac.uk/wiki/index.php?title=Talk:Mod:Hunt_Research_Group/ChemShell link]

===Research Notes===
#Cl- in water [https://www.ch.ic.ac.uk/wiki/index.php/Talk:Mod:Hunt_Research_Group/wannier_centre link]
#The use of Legendre time correlation functions to study reorientational dynamics in liquids[https://www.ch.ic.ac.uk/wiki/index.php/Talk:Mod:Hunt_Research_Group/legendre link]
#Functional for ILs using CPMD [https://wiki.ch.ic.ac.uk/wiki/index.php?title=Talk:Mod:Hunt_Research_Group/IL_cpmd_functional link]
#Solving the angular part of the Schrödinger equation for a hydrogen atom [https://wiki.ch.ic.ac.uk/wiki/index.php?title=Talk:Mod:Hunt_Research_Group/angular_schrodinger link] (notes by Vincent)
#Systematic conformational scan for ion-pair dimers [https://wiki.ch.ic.ac.uk/wiki/index.php?title=Talk:Mod:Hunt_Research_Group/ion_pair_scan link]
#Obtaining NBO, ESP, and RESP charges [https://wiki.ch.ic.ac.uk/wiki/index.php?title=Mod:Hunt_Research_Group/Charges]

===Installing and using other packages===
#How to install POLYRATE [https://www.ch.ic.ac.uk/wiki/index.php/Mod:Hunt_Research_Group/polyrate link]
#How to install Geomview [https://www.ch.ic.ac.uk/wiki/index.php/Mod:Hunt_Research_Group/geomview link]
#XMGRACE, gfortran, c compilers for Lion [http://hpc.sourceforge.net/]

===Admin Stuff===
#Not used to writing a wiki, make your test runs [https://www.ch.ic.ac.uk/wiki/index.php/Mod:Hunt_Research_Group/testing on this page]
#How to set-up new macs [https://www.ch.ic.ac.uk/wiki/index.php/Mod:Hunt_Research_Group/mac_setup link]
#How to set-up remote desktop [https://www.ch.ic.ac.uk/wiki/index.php/Mod:Hunt_Research_Group/mac_remote link]
#[https://www.ch.ic.ac.uk/wiki/index.php/Mod:Hunt_Research_Group/calendar Calendar]
#How to switch the printer HP CP3525dn duplex on and off [https://www.ch.ic.ac.uk/wiki/index.php/Mod:Hunt_Research_Group/printing link]

Mod:Hunt Research Group/Charges

2017-06-06T10:36:41Z

Jc6613:

This article outlines how to obtain NBO, ESP and RESP-fitted charges.

== NBO ==

== ESP ==
[[File:Resp flowchart.png|none|thumb|1103x1103px]]

Mod:Hunt Research Group/Charges

2017-06-06T10:34:27Z

Jc6613: Created page with "1103x1103px"

[[File:Resp flowchart.png|none|thumb|1103x1103px]]

File:Resp flowchart.png

2017-06-06T10:33:16Z

Jc6613:

2016-02-11T21:06:00Z

Jc6613: /* Analysis */

This computational experiment aims to investigate thermodynamic properties of the MgO crystal lattice, described as either a rhombohedron primitive unit cell with a 60˚ internal angle, or a conventional face-centred cubic cell with a 90˚ internal angle. Different approaches to computational simulation are employed, and a summary of the thermodynamic properties and a comparison across the two different approaches are provided.

== Initial Calculation on MgO ==
An initial calculation on a single unit cell of MgO at absolute zero yields a Helmholtz free energy of –41.07531835 eV, or the equivalent of -3957.689074 kJ mol-1. The lattice parameter, a = b = c, of a primitive rhombohedric cell is 2.9783 Å, or the equivalent of 2.5792 Å in a conventional face-centred cubic cell (NaCl structure). The lattice parameter is in good agreeement with the LCAO-HF calculation (linear combination of atomic orbitals at the Hatree-Fock level) found in literature (2.573 Å <ref name=":0">''II-VI and I-VII Compounds; Semimagnetic Compounds'', Springer Berlin Heidelberg, Berlin, 1999, 41B, pp. 1-6, doi: [http://link.springer.com/chapter/10.1007%2F10681719_206 10.1007/10681719_206]</ref>), suggesting that the GULP shell-model simulation is in good agreement with quantum mechanical calculations.

== The Phonon Dispersion Relation In MgO ==

=== Phonons ===
A phonon is an excitation of a group of atoms in an elastic arrangement in a crystal lattice. While seen as a quasiparticle, a phonon represents an excited state in a quantised quantum mechanical mode of vibration. Understanding the behavioral changes of phonons at different temperatures is crucial in the study of thermal expansion, thermal conductivity and electrical conductivity of crystals.

The displacement of one atom in a crystal lattice leads to a propagating wave of the lattice, where the amplitude of the wave is the displacement of the atom from its equilibrium position, and the wavelength (λ) is the distance between consecutive points of the same amplitude (phase). The wavenumber is more useful in lattice representations, with <math>k = \frac{2 \pi}{\lambda}</math>.

=== Results and Analysis ===
[[File:MgO Phonon Dispersion jc6613.png|frame|351x351px|Figure 1: Phonon Dispersion Relation of MgO.]]
Using GULP, phonon frequencies were calculated at various k-points along the paths Γ→X→W→K and Γ→L→W within the first Brillouin zone. A total of 50 k-points were sampled to give Figure 1.

The shape of the computed phonon dispersion relation is in good agreement with experimental dispersion relations <ref>G. Peckham, Proc. Phys. Soc. 1967, 90, 657, doi: [http://iopscience.iop.org/article/10.1088/0370-1328/90/3/312/pdf 10.1088/0370-1328/90/3/312]</ref>.

By visualising the vibrational modes with the DLVisualize Animate feature, it is understood that the lowest two branches at Γ are overlapping acoustic transverse vibrations, and the second-lowest branch is the longitudinal vibration. Close to the Γ-symmetry k-point, the energy of acoustic phonons is approximately proportional to the wavenumber, <math> E = \frac{h \omega_0 k a}{2 \pi} </math>. The acoustic transverse vibrations in two perpendicular planes overlap to a greater degree in this calculation than that in experimental literature. Degree-of-freedom arguments suggest that 3N-3 optical vibrations exist above the acoustic modes (N = number of atoms in a unit cell), which holds true since there are 3 optical branches in the phonon dispersion relation graph.

== Phonon Density-of-States Functions In MgO ==

=== Phonon Density of States ===
The Phonon density of states function outlines the number of vibrational states available at each frequency. For a unisolated system, the density distribution of states is continuous.

Analytically the density of states can be obtained by integrating the phonon energy dispersion at an infinitesimally small increment of k. When the branch is flat at the point of interest, there is a high density of vibrational state for the particular k-point, or for the particular energy since energy is a function of the wavenumber k.

When full analytical integration is not available, a numerical integration can be performed by GULP by calculating the vibrational modes at each k-point and summing all densities multiplied by the weight of the k-point. The weighting factor, for a regular grid, is the inverse of the number of grid points. As the number of points approach infinity, the summation resembles the true value.

In practical calculations, a finite number of k-points is sampled over the first Brillouin zone. A Monkhorst-Pack mesh is an unbiased method of choosing a set of k-points uniformly spread over the first Brillouin zone. The uniform mesh is suitable for a regular grid such as MgO, where there are points of special interest are spread evenly. The number of points is defined by shrinking factors in three dimensions and the larger the shrinking factors, the finer and more accurate the DOS function.

=== Results ===
{| class="wikitable"
!Shrinking Factor
!MgO Phonon DOS
!Number of Sampled k Points
!Noticeable Maxima (inverse cm)
|-
|1x1x1
|[[File:MgO Phonon DOS (1x1).png|frameless]]
|'''1:'''

''1/2, 1/2, 1/2 (L)''

|256

381

688

805
|-
|2x2x2
|[[File:MgO Phonon DOS (2x2).png|frameless]]
|'''4:'''

''1/4, 1/4, 1/4''

''1/4, 1/4, 3/4;''

''1/4, 3/4, 3/4;''

''3/4, 1/4, 3/4''
|384

and four other individual peaks
|-
|4x4x4
|[[File:MgO Phonon DOS (4x4).png|frameless]]
|'''32'''

''1/8, 1/8, 1/8''

''1/8, 1/8, 3/8''

''1/8, 1/8, 5/8''

''1/8, 1/8, 7/8;''

''then repeat this sequence for all combinations of k in all directions''
|408

and eleven other individual peaks
|-
|8x8x8
|[[File:MgO Phonon DOS (8x8).png|frameless]]
|'''256'''
|432
|-
|16x16x16
|[[File:MgO Phonon DOS (16x16).png|frameless]]
|'''2048'''
|396
|-
|32x32x32
|[[File:MgO Phonon DOS (32x32).png|frameless]]
|'''16,384'''
|396
|}

The numbers of k points given in the table do not account for intrinsic symmetry of phonon modes in the Brillouin zone. Hence, the actual number of sampling points performed by GULP may be lower than those given in the table.

=== Analysis ===
As the shrinking factor increases, the DOS functions become more and more continuous. This is natural since the vibrational states should become more continuous for increasing size of the system. The shape of the DOS function finer than the 16x16x16 system agrees with calculations found in literature and experimental data<ref>A. R. Oganov, M. J. Gillan, G. D. Price, J. Chem. Phys. 2003, 118 (22), pp. 10174, doi: [http://uspex.stonybrook.edu/pdfs/JCP-2003.pdf 10.1063/1.1570394]</ref>. Systems with smaller grids display abrupt changes at peaks and troughs, while experimental studies suggest a smooth transition of densities between different energies. As the grid size increases, the DOS becomes closer to experimental values as the function is smoother and representative of more k-points, not just the ones sampled. Maximum occurs around the 400 cm-1 frequency as expected because the phonon dispersion curves are the most 'crowded' at that particular frequency.

We can essentially observe the process of summation of the numerical integration as we increase the shrinking factor. For the 1x1x1 grid, only the L symmetry point is sampled, and the four peaks correspond to the frequencies at which the phonon dispersion branches intersect the L symmetry label. The two high-frequency peaks at 688 and 805 cm-1 correspond to the two points of intersection with the optical branches, and the lower-frequency peaks correspond to that with the acoustic branches. At L, the acoustic branches are flatter, giving a larger number of vibrational energies. As the number of sampling points increase, more and more points of intersection and frequencies are added to the DOS curve, giving a much smoother and more continuous curve.

The illustration below outlines the relationship between the phonon dispersion curve and the density-of-states curve:

[[File:MgOPhononOverview jc6613.png|frameless|612x612px]]

=== Applicability To Other Crystal Lattices ===
MgO is a II-VI ionic crystal with a face-centred cubic cell structure (a=b=c=4.211 Å). To ensure that we calculate the same number of vibrational states, the reciprocal space separation between sampled k-points must be similar to obtain the same resolution. Theoretically, the same grid size can be applied to lattices of similar cell parameters and the same ionic charges, such as CaO (a=b=c= 4.8Å) and such computational exercise has been carried out in literature with success <ref>B. B. Karki, R.M. Wentzcovitch, Phys. Rev. B, 2003, 68 (22), pp. 224304, doi: [http://journals.aps.org/prb/pdf/10.1103/PhysRevB.68.224304 10.1103/PhysRevB.68.224304]</ref>. For lattices with different ions, however, the grid size must be changed depending on the nature of the ion. In the MgO lattice, the high polarisability of the oxide anion is suitable for the computational model. While, for example in metallic lithium (a=b=c=3.521 Å), the more delocalised nature of the valence electrons imply large energy variation in a given vibrational state, and a larger grid size may be necessary to obtain the same resolution, in addition to the cell parameter effects. The smaller size of the cell parameter also requires a larger grid in the reciprocal space. However, computations have been carried out with limited success <ref>S. Pal, Phys. Rev. B, 1970, 2 (12), pp. 4741, doi: [http://journals.aps.org/prb/pdf/10.1103/PhysRevB.2.4741 10.1103/PhysRevB.2.4741]</ref>. For crystals with a larger, a smaller grid size is required to sample all points of interest within the symmetry of the unit cell. Computational calculations for such complex lattices, such as a zeolite, will be very costly.

== The Helmholtz Free Energy of MgO ==

=== Quasi-Harmonic Approximation ===
Even at low vibrational states, the energy of a vibrational state deviates from perfect harmonic behaviour. To introduce anharmonicity to the simulation, the quasi-harmonic approximation computes the positions and the bond lengths of individual atoms for vibrations corresponding to different wavenumbers. As the lattice compresses or expands, the zero-point energy (electrostatic) and the force constants vary. A new lattice is then computed at various k-points using the new force constant to minimise the overall free energy. The force constant calculation requires an accurate phonon frequency, hence a thorough sampling of k-points is required.

<math> A = U + \sum_i 0.5 \hbar \omega_i + k_B T \sum_i \ln \left( 1 - e^{-\frac{\hbar \omega_i}{k_B T}} \right)</math>, where A is the Helmholtz energy, U is the internal energy, kT is the Boltzmann constant, h-bar is the reduced Plank constant, ω is the vibrational frequency, T is the temperature.

Since the optimisation of the new force constant are independent of other vibrations, the quasi-harmonic approximation breaks down at the high-temperature range because of significant phonon–phonon interactions. Another shortcoming of the quasi-harmonic approximation, due to its intrinsic harmonic nature even with different force constants, is the failure to model dissociation of the bonds at high internuclear distances and strong repulsions at high compressions, further invalidating this approach for high-temperature conditions.

=== Results and Analysis ===
{| class="wikitable"
!Shrinking
Factor
!Number of
Sampled k Points
!Zero-Point Energy
'''(eV)'''
!Helmholtz Free Energy
(eV)
!Heat Capacity At Constant Volume
(eV/K)
|-
|1x1x1
|1
|0.172063
|'''-40.93'''0301
|0.000350
|-
|2x2x2
|4
|0.174209
|'''-40.92'''6609
|0.000347
|-
|4x4x4
|32
|0.174339
|'''-40.926'''450
|0.000347
|-
|8x8x8
|256
|0.174340
|'''-40.9264'''78
|0.000347
|-
|16x16x16
|2048
|0.174340
|'''-40.92648'''2
|0.000347
|-
|32x32x32
|16,384
|0.174340
|'''-40.92648'''4
|0.000347
|}
Overall the Helmholtz free energy is more positive as the shrinking factor is increased, which is also the trend of the zero-point energy. However a noticeable maxima occurs at the 2x2x2 system, where the energy peaks before dropping off to convergence. Because the zero-point energy trend lacks this feature, it is speculated that the peak in Helmholtz free energy is due to entropic contributions. As spatial sampling becomes finer, the TS term in Helmholtz free energy becomes increasingly negative (and more accurate). This is intuitive since a small shrinking factor does not account for the many modes of vibrations that contribute to overall entropy.

The 2x2x2 grid is suitable for calculations accurate to 0.5 meV per cell, which implies that it is also suitable for calculations accurate to 1 meV. For 0.1 meV accurate calculations, a 4x4x4 grid size is required.

=== Applicability To Other Crystal Lattices ===
Because zero-point internal energy calculations will be identical regardless of the sampling, we should focus on the entropic term when discussing the accuracy in calculating the free energy for different crystal lattices. The entropic term is dependent on accurate energy calculations when occupying the vibrational states. The accuracy of the entropic term, and therefore of the whole Helmholtz free energy is reliant on the accuracy of the density of states. The arguments written in the previous section can be applied here. In short, the same grid size can be applied to a CaO calculation, a larger grid size must be used for metallic lithium, and a smaller grid size is required for zeolites to obtain the same resolution.

== Thermal Expansion of MgO: Quasi-Harmonic Approximation ==

The coefficient of thermal expansion is a description of the change of an object's size with a change in temperature. More specifically, the expansion is described as a fractional change in size per increment of temperature. In this computational experiment, we are considered with the volume expansion. For a isotropic material such as the MgO crystal lattice, the thermal expansion coefficient of the object in 3D is three times the linear thermal expansion coefficient in 1D. Mathematically the volume thermal expansion coefficient is defined as <math>\alpha_V = \frac{1}{V} \left( \frac{\partial V}{\partial T} \right)_p</math> , where αV is the volume thermal expansion coefficient, V is the total volume of the system, T is the temperature of the system.

The physical nature of a crystal expanding upon heat is the promotion of electrons onto higher-energy vibrational states with more anti-bonding character. The anti-bonding nature of the state increases the bond length, leading to the expansion of the whole lattice. For a lattice system with large energy variation between Γ and K, a large amount of thermal energy is required for the same change in volume than a system with vibrational states that are closely packed and similar in energy. If a lattice has large separations between its vibrational states then nature of the bonding in the lattice can be understood as structurally rigid, which fits with our chemical intuition that molecules with strong and rigid bonds are less prone to thermal expansion.

=== Computational Methodology ===
At each given temperature, an optimal configuration of the system is simulated by the method described in the previous section. The Helmholtz free energy is the system is minimised to a pre-set tolerance with respect to the cell volume. According to the variational principle, thermodynamic properties of a system is closest to the true value when the derivative of cell parameters with respect to temperature is zero.

=== Results and Analysis ===
The 8x8x8 lattice grid was used for thermal optimisations for free energy values that are accurate to 0.01 meV.
{| class="wikitable"
!Temperature
(K)
!Helmholtz Free Energy
(eV)
!Lattice Constant
(a=b=c, Angstroms)
!Angle
(α=β=γ, Degrees)
!Cell Parameter
Derivatives
|-
|0
|<nowiki>-40.901906</nowiki>
|2.986565
|59.999934
|0
|-
|100
|<nowiki>-40.902418</nowiki>
|2.986564
|59.999937
|0
|-
|200
|<nowiki>-40.909374</nowiki>
|2.987613
|59.999639
|0
|-
|300
|<nowiki>-40.928119</nowiki>
|2.989403
|59.999461
|<nowiki>-0.000002</nowiki>
|-
|400
|<nowiki>-40.958587</nowiki>
|2.989403
|59.992790
|0
|-
|500
|<nowiki>-40.999427</nowiki>
|2.994156
|59.999094
|0
|-
|600
|<nowiki>-41.049305</nowiki>
|2.996846
|59.998906
|0
|-
|700
|<nowiki>-41.049305</nowiki>
|2.996846
|59.998906
|0
|-
|800
|<nowiki>-41.107108</nowiki>
|2.999674
|59.998715
|0
|-
|900
|<nowiki>-41.171879</nowiki>
|3.002623
|59.998520
|0
|-
|1000
|<nowiki>-41.243004</nowiki>
|3.005674
|59.998323
|0
|-
|1200
|<nowiki>-41.488740</nowiki>
|3.015394
|59.997763
|0
|-
|1400
|<nowiki>-41.675515</nowiki>
|3.022433
|59.997451
|0.000001
|-
|1600
|<nowiki>-41.877981</nowiki>
|3.029988
|59.996775
|0
|-
|1800
|<nowiki>-42.094435</nowiki>
|3.038336
|59.995015
|0.000001
|-
|2000
|<nowiki>-42.323668</nowiki>
|3.046923
|59.994762
|<nowiki>-0.000002</nowiki>
|-
|2200
|<nowiki>-42.564747</nowiki>
|3.056841
|59.994783
|0
|-
|2400
|<nowiki>-42.817003</nowiki>
|3.068085
|59.992390
|<nowiki>-0.00002</nowiki>
|-
|2600
|<nowiki>-43.079967</nowiki>
|3.081449
|59.993314
|0.000002
|-
|2752
|<nowiki>-43.28700906</nowiki>
|3.094344
|59.994204
|0.000002
|-
|2852
|<nowiki>-43.42654086</nowiki>
|3.105499
|59.994105
|<nowiki>-0.000001</nowiki>
|-
|2952
|<nowiki>-43.56951767</nowiki>
|3.124113
|59.994437
|<nowiki>-0.000002</nowiki>
|}
The lattice constant increases linearly only in the range 300 - 1000 K.[[File:Screenshot-Thermal Expansion.xls - LibreOffice Calc.png|none|frame]]
The thermal expansion constant of MgO in this temperature range, assuming temperature independence, is equal to <math>8.74 \times 10^{-6} \text{K}^{-1}</math> with <math>R^2 = 0.99</math>.

[[File:Screenshot-Thermal Expansion.xls - LibreOffice Calc ALL.png|none|frame|Figure 3: Thermal expansion of MgO in all simulated temperatures.]]

If the whole temperature range 0K to 2952K was to be included, the thermal expansion of MgO is non-linear with respect to temperature. A polynomial fit to the fourth order has been determined with <math> R^2 = 1.00 </math> such that <math>\frac{a}{a_{0}} = 2.27 \times 10^{-5} T^4 - 1.04 \times 10^{-11} T^3 + 1.64 \times 10^{-8} T^2 - 8.35 \times 10^{-7} T - 2.60 \times 10^{-5}</math>, where <math>\frac{a}{a_{0}}</math> is the fractional change in the cell parameter relative to the cell parameter at T = 0 and <math> T </math> is the temperature in Kelvins. As the coefficient of thermal expansion is defined as the partial derivative of the fractional change with respect to temperature, the coefficient of thermal expansion can be expressed as <math> \alpha (T) = 9.08 \times 10^{-15} T^3 - 4.16 \times 10^{-11} T^2 + 3.28 \times 10^{-8} T - 5.85 \times 10^{-6} </math>.

== Thermal Expansion of MgO: Molecular Dynamics==

=== Computational Methodology ===
The molecular dynamics technique studies the trajectories of individual atoms by numerically solving Newton's equations of motion for an N-body system of interacting particles, where the forces of the particles are calculated by molecular mechanical force fields such as the Lennard-Jones or the Buckingham potential. It is an ''ab initio'' method of simulation applied to a large number of unit cells in order to obtain an accurate picture of the complex system in which analytical solutions are either inhibitively time-consuming or simply unavailable.

Thermodynamic properties of the isothermal-isobaric (NPT) ensemble are time-averaged over the entire space. For ergodic systems, the time averages should be identical to spatial averages provided by the quasi-harmonic approach. Due to the computational expense of this method, the system must be designed to account for available computational power by tailoring the simulation to a suitable timestep, time duration and scale. The scale of such a system must be large enough to eliminate periodic boundary condition artifacts, which in this case is suitable at the femtoseconds timestep with a lattice containing 32 MgO unit cells. This number of cells in the real space is equivalent to the 8x8x8 shrinking factor in the Helmholtz free energy simulation, as both methods give a total of 64 vibrational bands. This shrinking factor was found to give an accurate Helmholtz free energy correct to 0.1 meV.

=== Results and Analysis ===
{| class="wikitable"
!Temperature
(K)
!Averaged Final Volume
(Cubic Angstroms)
!Volume Per Unit Cell
(Cubic Angstroms)
|-
|300
|602.917248
|18.8411640
|-
|400
|604.188880
|18.8809025
|-
|500
|605.968126
|18.9365039
|-
|600
|608.794097
|19.0248155
|-
|700
|608.745756
|19.0233049
|-
|800
|611.617711
|19.1130535
|-
|900
|611.936388
|19.1230121
|-
|1000
|615.044932
|19.2347085
|-
|1100
|619.298570
|19.3530803
|-
|1200
|620.020014
|19.3756254
|-
|1300
|621.486358
|19.4214487
|-
|1400
|622.667288
|19.4583528
|-
|1600
|626.172143
|19.5678795
|-
|1800
|630.981409
|19.7181690
|-
|2000
|632.410325
|19.7628227
|-
|2400
|642.622061
|20.0819394
|-
|2600
|648.508557
|20.2658924
|-
|2800
|655.021305
|20.6318792
|-
|3000
|660.220134
|20.6318792
|}
[[File:Thermal Expansion MD v QHA.png|none|frame|Figure 4: A graph showing expansion of unit cell with increasing temperature, where the blue series is simulated by quasi-harmonic approximation, and the orange series is simulated by molecular dynamics. ]]The molecular dynamics data suggest a linear relationship between unit cell volume and temperature, and the thermal expansion coefficient, assuming its temperature independence is <math>\alpha_V = 3.4599078 \times 10^{-5} \text{K}^{-1}</math>.

Due to time constraints, I was unable to perform a molecular dynamics simulation at 0K. Extrapolating from obtained data, the volume per unit cell at absolute zero is 18.594625 ‎Å3, whereas with quasi-harmonic approximation the 0 K volume is 17.302424 Å3. It is evident that the quasi-harmonic approximation, with its inaccurate parabolic potential curve and boundary condition artifacts, significantly underestimate the magnitude of repulsive forces in a crystal lattice, giving rise to largely inaccurate high-temperature properties.

== Conclusion ==
The experimental thermal coefficient of MgO varies from <math>31.2 \times 10^{-6} \text{K}^{-1}</math> at 300 K to <math>53.3 \times 10^{-6} \text{K}^{-1}</math> at 2000 K, steadily increasing with increasing temperature.<ref>C. J. Engberg, E. H. Zehms, J. Am. Ceram. Soc., 1959, 42 (6), pp. 300, doi: [http://onlinelibrary.wiley.com/store/10.1111/j.1151-2916.1959.tb12958.x/asset/j.1151-2916.1959.tb12958.x.pdf?v=1&t=ikfzf86u&s=3040cd204ec490ffeaa1cabf505c7097de78aed2&systemMessage=Wiley+Online+Library+will+be+unavailable+on+Saturday+27th+February+from+09%3A00-14%3A00+GMT+%2F+04%3A00-09%3A00+EST+%2F+17%3A00-22%3A00+SGT+for+essential+maintenance.++Apologies+for+the+inconvenience. 10.1111/j.1151-2916.1959.tb12958.x]</ref><ref>O.L. Anderson, K. Zou, J. Phys. Chem. Ref. Data, 1990, 19 (1), pp. 70</ref>

Both molecular dynamics and quasi-harmonic simulations display a linear thermal expansion in the temperature range 0 K – 1200 K. The quasi-harmonic method is in better agreement with experimental data in the linear range, after which the thermal expansion diverges. In the high-temperature range, thermal expansion no longer obeys linearity according to quasi-harmonic simulations as molecules are no longer bound to harmonic force fields due to increased energies. The quasi-harmonic calculation does not simulate bond dissociation as a result of large energies in the high-temperature range, and fails to account for anharmonicity experienced at high bond distances.

On the other hand, molecular dynamics indicates linearity over a much larger range of temperature. At lower temperatures where periodic condition artifacts are more prominent, molecular dynamics simulations contain some inaccuracies within its somewhat simplistic classical force models. At high temperatures, molecular dyanmics is more reliable in calculating phonon interactions, and therefore more suited to high-temperature simulations.

For both methods, the calculations break down at and near the melting point of MgO. In an NPT system, the thermal expansion coefficient and the constant-pressure heat capacity are both proportional to the rate of change in vibrational energy with respect to temperature. We expect the heat capacity to instantaneously become infinite at the transition temperature. However it is not the case as both models predict constantly rising thermal expansion coefficients beyond the transition temperature. It is suspected that the long-range atomic forces are not well employed in both methods of simulation, and further studies with more descriptive potential fields are needed to investigate the thermal behavior of the MgO lattice.

== References ==

2016-02-11T14:06:58Z

Jc6613:

Rep:Mod:pitsnake

2016-02-11T13:59:35Z

Jc6613: /* Computational Methodology */

This computational experiment aims to investigate thermodynamic properties of the MgO crystal lattice, described as either a rhombohedron primitive unit cell with a 60˚ internal angle, or a conventional face-centred cubic cell with a 90˚ internal angle. Different approaches to computational simulation are employed, and a summary of the thermodynamic properties and a comparison across the two different approaches are provided.

== Initial Calculation on MgO ==
An initial calculation on a single unit cell of MgO at absolute zero yields a Helmholtz free energy of –41.07531835 eV, or the equivalent of -3957.689074 kJ mol-1. The lattice parameter, a = b = c, of a primitive rhombohedric cell is 2.9783 Å, or the equivalent of 2.5792 Å in a conventional face-centred cubic cell (NaCl structure). The lattice parameter is in good agreeement with the LCAO-HF calculation (linear combination of atomic orbitals at the Hatree-Fock level) found in literature (2.573 Å <ref name=":0">''II-VI and I-VII Compounds; Semimagnetic Compounds'', Springer Berlin Heidelberg, Berlin, 1999, 41B, pp. 1-6, doi: [http://link.springer.com/chapter/10.1007%2F10681719_206 10.1007/10681719_206]</ref>), suggesting that the GULP shell-model simulation is in good agreement with quantum mechanical calculations.

== The Phonon Dispersion Relation In MgO ==

=== Phonons ===
A phonon is an excitation of a group of atoms in an elastic arrangement in a crystal lattice. While seen as a quasiparticle, a phonon represents an excited state in a quantised quantum mechanical mode of vibration. Understanding the behavioral changes of phonons at different temperatures is crucial in the study of thermal expansion, thermal conductivity and electrical conductivity of crystals.

The displacement of one atom in a crystal lattice leads to a propagating wave of the lattice, where the amplitude of the wave is the displacement of the atom from its equilibrium position, and the wavelength (λ) is the distance between consecutive points of the same amplitude (phase). The wavenumber is more useful in lattice representations, with <math>k = \frac{2 \pi}{\lambda}</math>.

=== Results and Analysis ===
Using GULP, phonon frequencies were calculated at various k-points along the paths Γ→X→W→K and Γ→L→W within the first Brillouin zone. A total of 50 k-points were sampled to give the graph shown below.

'''[INSERT Phonon Dispersion]'''

The shape of the computed phonon dispersion relation is in good agreement with experimental dispersion relations <ref>G. Peckham, Proc. Phys. Soc. 1967, 90, 657, doi: [http://iopscience.iop.org/article/10.1088/0370-1328/90/3/312/pdf 10.1088/0370-1328/90/3/312]</ref>.

By visualising the vibrational modes with the DLVisualize Animate feature, it is understood that the lowest two branches at Γ are overlapping acoustic transverse vibrations, and the second-lowest branch is the longitudinal vibration. Close to the Γ-symmetry k-point, the energy of acoustic phonons is approximately proportional to the wavenumber, <math> E = \frac{h \omega_0 k a}{2 \pi} </math>. The acoustic transverse vibrations in two perpendicular planes overlap to a greater degree in this calculation than that in experimental literature. Degree-of-freedom arguments suggest that 3N-3 optical vibrations exist above the acoustic modes (N = number of atoms in a unit cell), which holds true since there are 3 optical branches in the phonon dispersion relation graph.

== Phonon Density-of-States Functions In MgO ==

=== Phonon Density of States ===
The Phonon density of states function outlines the number of vibrational states available at each frequency. For a unisolated system, the density distribution of states is continuous.

Analytically the density of states can be obtained by integrating the phonon energy dispersion at an infinitesimally small increment of k. When the branch is flat at the point of interest, there is a high density of vibrational state for the particular k-point, or for the particular energy since energy is a function of the wavenumber k.

When full analytical integration is not available, a numerical integration can be performed by GULP by calculating the vibrational modes at each k-point and summing all densities multiplied by the weight of the k-point. The weighting factor, for a regular grid, is the inverse of the number of grid points. As the number of points approach infinity, the summation resembles the true value.

In practical calculations, a finite number of k-points is sampled over the first Brillouin zone. A Monkhorst-Pack mesh is an unbiased method of choosing a set of k-points uniformly spread over the first Brillouin zone. The uniform mesh is suitable for a regular grid such as MgO, where there are points of special interest are spread evenly. The number of points is defined by shrinking factors in three dimensions and the larger the shrinking factors, the finer and more accurate the DOS function.

=== Results ===
{| class="wikitable"
!Shrinking Factor
!MgO Phonon DOS
!Number of Sampled k Points
!Noticeable Maxima (inverse cm)
|-
|1x1x1
|[[File:MgO Phonon DOS (1x1).png|frameless]]
|'''1:'''

''1/2, 1/2, 1/2 (L)''

|256

381

688

805
|-
|2x2x2
|[[File:MgO Phonon DOS (2x2).png|frameless]]
|'''4:'''

''1/4, 1/4, 1/4''

''1/4, 1/4, 3/4;''

''1/4, 3/4, 3/4;''

''3/4, 1/4, 3/4''
|384

and four other individual peaks
|-
|4x4x4
|[[File:MgO Phonon DOS (4x4).png|frameless]]
|'''32'''

''1/8, 1/8, 1/8''

''1/8, 1/8, 3/8''

''1/8, 1/8, 5/8''

''1/8, 1/8, 7/8;''

''then repeat this sequence for all combinations of k in all directions''
|408

and eleven other individual peaks
|-
|8x8x8
|[[File:MgO Phonon DOS (8x8).png|frameless]]
|'''256'''
|432
|-
|16x16x16
|[[File:MgO Phonon DOS (16x16).png|frameless]]
|'''2048'''
|396
|-
|32x32x32
|[[File:MgO Phonon DOS (32x32).png|frameless]]
|'''16,384'''
|396
|}

The numbers of k points given in the table do not account for intrinsic symmetry of phonon modes in the Brillouin zone. Hence, the actual number of sampling points performed by GULP may be lower than those given in the table.

=== Analysis ===
As the shrinking factor increases, the DOS functions become more and more continuous. This is natural since the vibrational states should become more continuous for increasing size of the system. The shape of the DOS function finer than the 16x16x16 system agrees with calculations found in literature and experimental data<ref>A. R. Oganov, M. J. Gillan, G. D. Price, J. Chem. Phys. 2003, 118 (22), pp. 10174, doi: [http://uspex.stonybrook.edu/pdfs/JCP-2003.pdf 10.1063/1.1570394]</ref>. Maximum occurs around the 400 cm-1 frequency as expected because the phonon dispersion curves are the most 'crowded' at that particular frequency.

We can essentially observe the process of summation of the numerical integration as we increase the shrinking factor. For the 1x1x1 grid, only the L symmetry point is sampled, and the four peaks correspond to the frequencies at which the phonon dispersion branches intersect the L symmetry label. The two high-frequency peaks at 688 and 805 cm-1 correspond to the two points of intersection with the optical branches, and the lower-frequency peaks correspond to that with the acoustic branches. At L, the acoustic branches are flatter, giving a larger number of vibrational energies. As the number of sampling points increase, more and more points of intersection and frequencies are added to the DOS curve, giving a much smoother and more continuous curve.

=== Applicability To Other Crystal Lattices ===
'''ASK – Polarizability O2-, lattice parameters'''

CaO<ref>B. B. Karki, R.M. Wentzcovitch, Phys. Rev. B, 2003, 68 (22), pp. 224304, doi: [http://journals.aps.org/prb/pdf/10.1103/PhysRevB.68.224304 10.1103/PhysRevB.68.224304]</ref> Li<ref>S. Pal, Phys. Rev. B, 1970, 2 (12), pp. 4741, doi: [http://journals.aps.org/prb/pdf/10.1103/PhysRevB.2.4741 10.1103/PhysRevB.2.4741]</ref>

== The Helmholtz Free Energy of MgO ==

=== Quasi-Harmonic Approximation ===
Even at low vibrational states, the energy of a vibrational state deviates from perfect harmonic behavior. To introduce anharmonicity to the simulation, the quasi-harmonic approximation computes the positions and the bond lengths of individual atoms for vibrations corresponding to different wavenumbers. As the lattice compresses or expands, the zero-point energy (electrostatic) and the force constants vary. A new lattice is then computed at various k-points using the new force constant to minimise the overall free energy. The force constant calculation requires an accurate phonon frequency, hence a thorough sampling of k-points is required.

<math> A = U + \sum_i 0.5 \hbar \omega_i + k_B T \sum_i \ln \left( 1 - e^{-\frac{\hbar \omega_i}{k_B T}} \right)</math>, where A is the Helmholtz energy, U is the internal energy, kT is the Boltmann constant, h-bar is the reduced Plank constant, ω is the vibrational frequency, T is the temperature.

Since the optimisation of the new force constant are independent of other vibrations, the quasi-harmonic approximation breaks down at the high-temperature range because of significant phonon–phonon interactions. Another shortcoming of the quasi-harmonic approximation, due to its intrinsic harmonic nature even with different force constants, is the failure to model dissociation of the bonds at high internuclear distances and strong repulsions at high compressions, further invalidating this approach for high-temperature conditions.

=== Results and Analysis ===
{| class="wikitable"
!Shrinking
Factor
!Number of
Sampled k Points
!Zero-Point Energy
'''(eV)'''
!Helmholtz Free Energy
(eV)
!Heat Capacity At Constant Volume
(eV/K)
|-
|1x1x1
|1
|0.172063
|'''-40.93'''0301
|0.000350
|-
|2x2x2
|4
|0.174209
|'''-40.926'''609
|0.000347
|-
|4x4x4
|32
|0.174339
|'''-40.926'''450
|0.000347
|-
|8x8x8
|256
|0.174340
|'''-40.9264'''78
|0.000347
|-
|16x16x16
|2048
|0.174340
|'''-40.92648'''2
|0.000347
|-
|32x32x32
|16,384
|0.174340
|'''-40.92648'''4
|0.000347
|}
Overall the Helmholtz free energy is more positive as the shrinking factor is increased, which is also the trend of the zero-point energy. However a noticeable maxima occurs at the 2x2x2 system, where the energy peaks before dropping off to convergence. Because the zero-point energy trend lacks this feature, it is speculated that the peak in Helmholtz free energy is due to entropic contributions. As spatial sampling becomes finer, the TS term in Helmholtz free energy becomes increasingly negative (and more accurate). This is intuitive since a small shrinking factor does not account for the many modes of vibrations that contribute to overall entropy.

The 2x2x2 grid is suitable for calculations accurate to 0.5 meV per cell, which implies that it is also suitable for calculations accurate to 1 meV. For 0.1 meV accurate calculations, a 4x4x4 grid size is required.

=== Applicability To Other Crystal Lattices ===

== Thermal Expansion of MgO: Quasi-Harmonic Approximation ==

The coefficient of thermal expansion is a description of the change of an object's size with a change in temperature. More specifically, the expansion is described as a fractional change in size per increment of temperature. In this computational experiment, we are considered with the volume expansion. For a isotropic material such as the MgO crystal lattice, the thermal expansion coefficient of the object in 3D is three times the linear thermal expansion coefficient in 1D.

Mathematically the volume thermal expansion coefficient is defined as <math>\alpha_V = \frac{1}{V} \left( \frac{\partial V}{\partial T} \right)_p</math> , where αV is the volume thermal expansion coefficient at constant volume, V is the total volume of the system, T is the temperature of the system.

=== Computational Methodology ===
At each given temperature, an optimal configuration of the system is simulated by the method described in the previous section. The Helmholtz free energy is the system is minimised to a pre-set tolerance with respect to the cell volume. According to the variational principle, thermodynamic properties of a system is closest to the true value when the derivative of cell parameters with respect to temperature is zero.

=== Results and Analysis ===
The 8x8x8 lattice grid was used for thermal optimisations for free energy values that are accurate to 0.01 meV.
{| class="wikitable"
!Temperature
(K)
!Helmholtz Free Energy
(eV)
!Lattice Constant
(a=b=c, Angstroms)
!Angle
(α=β=γ, Degrees)
!Cell Parameter
Derivatives
|-
|0
|<nowiki>-40.901906</nowiki>
|2.986565
|59.999934
|0
|-
|100
|<nowiki>-40.902418</nowiki>
|2.986564
|59.999937
|0
|-
|200
|<nowiki>-40.909374</nowiki>
|2.987613
|59.999639
|0
|-
|300
|<nowiki>-40.928119</nowiki>
|2.989403
|59.999461
|<nowiki>-0.000002</nowiki>
|-
|400
|<nowiki>-40.958587</nowiki>
|2.989403
|59.992790
|0
|-
|500
|<nowiki>-40.999427</nowiki>
|2.994156
|59.999094
|0
|-
|600
|<nowiki>-41.049305</nowiki>
|2.996846
|59.998906
|0
|-
|700
|<nowiki>-41.049305</nowiki>
|2.996846
|59.998906
|0
|-
|800
|<nowiki>-41.107108</nowiki>
|2.999674
|59.998715
|0
|-
|900
|<nowiki>-41.171879</nowiki>
|3.002623
|59.998520
|0
|-
|1000
|<nowiki>-41.243004</nowiki>
|3.005674
|59.998323
|0
|-
|1200
|<nowiki>-41.488740</nowiki>
|3.015394
|59.997763
|0
|-
|1400
|<nowiki>-41.675515</nowiki>
|3.022433
|59.997451
|0.000001
|-
|1600
|<nowiki>-41.877981</nowiki>
|3.029988
|59.996775
|0
|-
|1800
|<nowiki>-42.094435</nowiki>
|3.038336
|59.995015
|0.000001
|-
|2000
|<nowiki>-42.323668</nowiki>
|3.046923
|59.994762
|<nowiki>-0.000002</nowiki>
|-
|2200
|<nowiki>-42.564747</nowiki>
|3.056841
|59.994783
|0
|-
|2400
|<nowiki>-42.817003</nowiki>
|3.068085
|59.992390
|<nowiki>-0.00002</nowiki>
|-
|2600
|<nowiki>-43.079967</nowiki>
|3.081449
|59.993314
|0.000002
|-
|2752
|<nowiki>-43.28700906</nowiki>
|3.094344
|59.994204
|0.000002
|-
|2852
|<nowiki>-43.42654086</nowiki>
|3.105499
|59.994105
|<nowiki>-0.000001</nowiki>
|-
|2952
|<nowiki>-43.56951767</nowiki>
|3.124113
|59.994437
|<nowiki>-0.000002</nowiki>
|}
The lattice constant increases linearly only in the range 300 - 1000 K.[[File:Screenshot-Thermal Expansion.xls - LibreOffice Calc.png|none|frame]]
The thermal expansion constant of MgO in this temperature range, assuming temperature independence, is equal to <math>8.74 \times 10^{-6} \text{K}^{-1}</math> with <math>R^2 = 0.99</math>.

[[File:Screenshot-Thermal Expansion.xls - LibreOffice Calc ALL.png|none|frame|Figure 3: Thermal expansion of MgO in all simulated temperatures.]]

If the whole temperature range 0K to 2952K was to be included, the thermal expansion of MgO is non-linear with respect to temperature. A polynomial fit to the fourth order has been determined with <math> R^2 = 1.00 </math> such that <math>\frac{a}{a_{0}} = 2.27 \times 10^{-5} T^4 - 1.04 \times 10^{-11} T^3 + 1.64 \times 10^{-8} T^2 - 8.35 \times 10^{-7} T - 2.60 \times 10^{-5}</math>, where <math>\frac{a}{a_{0}}</math> is the fractional change in the cell parameter relative to the cell parameter at T = 0 and <math> T </math> is the temperature in Kelvins. As the coefficient of thermal expansion is defined as the partial derivative of the fractional change with respect to temperature, the coefficient of thermal expansion can be expressed as <math> \alpha (T) = 9.08 \times 10^{-5} T^3 - 4.16 \times 10^{-11} T^2 + 3.28 \times 10^{-8} T - 5.85 \times 10^{-6} </math>.

== Thermal Expansion of MgO: Molecular Dynamics==

=== Computational Methodology ===
The molecular dynamics technique studies the trajectories of individual atoms by numerically solving Newton's equations of motion for an N-body system of interacting particles, where the forces of the particles are calculated by molecular mechanical force fields such as the Lennard-Jones or the Buckingham potential. It is an ''ab initio'' method of simulation applied to a large number of unit cells in order to obtain an accurate picture of the complex system in which analytical solutions are either inhibitively time-consuming or simply unavailable.

Thermodynamic properties of the isothermal-isobaric (NPT) ensemble are time-averaged over the entire space. For ergodic systems, the time averages should be identical to spatial averages provided by the quasi-harmonic approach. Due to the computational expense of this method, the system must be designed to account for available computational power by tailoring the simulation to a suitable timestep, time duration and scale. The scale of such a system must be large enough to eliminate periodic boundary condition artifacts, which in this case is suitable at the femtoseconds timestep with a lattice containing 32 MgO unit cells.

=== Results and Analysis ===
{| class="wikitable"
!Temperature
(K)
!Averaged Final Volume
(Cubic Angstroms)
!Volume Per Unit Cell
(Cubic Angstroms)
|-
|300
|602.917248
|18.8411640
|-
|400
|604.188880
|18.8809025
|-
|500
|605.968126
|18.9365039
|-
|600
|608.794097
|19.0248155
|-
|700
|608.745756
|19.0233049
|-
|800
|611.617711
|19.1130535
|-
|900
|611.936388
|19.1230121
|-
|1000
|615.0449322
|19.2347085
|-
|1100
|619.298570
|19.3530803
|-
|1200
|620.020014
|19.3756254
|-
|1300
|621.486358
|19.4214487
|-
|1400
|622.667288
|19.4583528
|-
|1600
|626.172143
|19.5678795
|-
|1800
|630.981409
|19.7181690
|-
|2000
|632.410325
|19.7628227
|-
|2400
|642.622061
|20.0819394
|-
|2600
|648.508557
|20.2658924
|-
|2800
|655.021305
|20.6318792
|-
|3000
|660.220134
|20.6318792
|}
[[File:Thermal Expansion MD v QHA.png|none|frame|Figure 4: A graph showing expansion of unit cell with increasing temperature, where the blue series is simulated by quasi-harmonic approximation, and the orange series is simulated by molecular dynamics. ]]The molecular dynamics data suggest a linear relationship between unit cell volume and temperature, and the thermal expansion coefficient, assuming its temperature independence is \alpha_V = 3.4599078 \times 10^{-5} \text{K}^{-1}.

Due to time constraints, I was unable to perform a molecular dynamics simulation at 0K. Extrapolating from obtained data, the volume per unit cell at absolute zero is 18.594625 ‎Å3, whereas with quasi-harmonic approximation the 0 K volume is 17.302424 Å3. It is evident that the quasi-harmonic approximation, with its inaccurate parabolic potential curve and boundary condition artifacts, significantly underestimate the magnitude of repulsive forces in a crystal lattice, giving rise to largely inaccurate high-temperature properties.

== Conclusion ==
'''[INSERT COMPARISON MD V QHA V EXPERIMENTAL'''<ref name=":0" /><ref>C. J. Engberg, E. H. Zehms, J. Am. Ceram. Soc., 1959, 42 (6), pp. 300, doi: [http://onlinelibrary.wiley.com/store/10.1111/j.1151-2916.1959.tb12958.x/asset/j.1151-2916.1959.tb12958.x.pdf?v=1&t=ikfzf86u&s=3040cd204ec490ffeaa1cabf505c7097de78aed2&systemMessage=Wiley+Online+Library+will+be+unavailable+on+Saturday+27th+February+from+09%3A00-14%3A00+GMT+%2F+04%3A00-09%3A00+EST+%2F+17%3A00-22%3A00+SGT+for+essential+maintenance.++Apologies+for+the+inconvenience. 10.1111/j.1151-2916.1959.tb12958.x]</ref>''']'''

Both molecular dynamics and quasi-harmonic simulations display a linear thermal expansion in the temperature range 0 K – 1200 K. The quasi-harmonic method is in better agreement with experimental data in the linear range, after which the thermal expansion diverges. In the high-temperature range, thermal expansion no longer obeys linearity according to quasi-harmonic simulations as molecules are no longer bound to harmonic force fields due to increased energies. The quasi-harmonic calculation does not simulate bond dissociation as a result of large energies in the high-temperature range, and fails to account for anharmonicity experienced at high bond distances.

On the other hand, molecular dynamics indicates linearity over a much larger range of temperature. At lower temperatures where periodic condition artifacts are more prominent, molecular dynamics simulations contain some inaccuracies within its somewhat simplistic classical force models. At high temperatures, molecular dyanmics is more reliable in calculating phonon interactions, and therefore more suited to high-temperature simulations.

For both methods, the calculations break down at and near the melting point of MgO. In an NPT system, the thermal expansion coefficient and the constant-pressure heat capacity are both proportional to the rate of change in vibrational energy with respect to temperature. We expect the heat capacity to instantaneously become infinite at the transition temperature. However it is not the case as both models predict constantly rising thermal expansion coefficients beyond the transition temperature.

Rep:Mod:pitsnake

2016-02-11T13:58:59Z

Jc6613: /* Results and Analysis */

This computational experiment aims to investigate thermodynamic properties of the MgO crystal lattice, described as either a rhombohedron primitive unit cell with a 60˚ internal angle, or a conventional face-centred cubic cell with a 90˚ internal angle. Different approaches to computational simulation are employed, and a summary of the thermodynamic properties and a comparison across the two different approaches are provided.

== Initial Calculation on MgO ==
An initial calculation on a single unit cell of MgO at absolute zero yields a Helmholtz free energy of –41.07531835 eV, or the equivalent of -3957.689074 kJ mol-1. The lattice parameter, a = b = c, of a primitive rhombohedric cell is 2.9783 Å, or the equivalent of 2.5792 Å in a conventional face-centred cubic cell (NaCl structure). The lattice parameter is in good agreeement with the LCAO-HF calculation (linear combination of atomic orbitals at the Hatree-Fock level) found in literature (2.573 Å <ref name=":0">''II-VI and I-VII Compounds; Semimagnetic Compounds'', Springer Berlin Heidelberg, Berlin, 1999, 41B, pp. 1-6, doi: [http://link.springer.com/chapter/10.1007%2F10681719_206 10.1007/10681719_206]</ref>), suggesting that the GULP shell-model simulation is in good agreement with quantum mechanical calculations.

== The Phonon Dispersion Relation In MgO ==

=== Phonons ===
A phonon is an excitation of a group of atoms in an elastic arrangement in a crystal lattice. While seen as a quasiparticle, a phonon represents an excited state in a quantised quantum mechanical mode of vibration. Understanding the behavioral changes of phonons at different temperatures is crucial in the study of thermal expansion, thermal conductivity and electrical conductivity of crystals.

The displacement of one atom in a crystal lattice leads to a propagating wave of the lattice, where the amplitude of the wave is the displacement of the atom from its equilibrium position, and the wavelength (λ) is the distance between consecutive points of the same amplitude (phase). The wavenumber is more useful in lattice representations, with <math>k = \frac{2 \pi}{\lambda}</math>.

=== Results and Analysis ===
Using GULP, phonon frequencies were calculated at various k-points along the paths Γ→X→W→K and Γ→L→W within the first Brillouin zone. A total of 50 k-points were sampled to give the graph shown below.

'''[INSERT Phonon Dispersion]'''

The shape of the computed phonon dispersion relation is in good agreement with experimental dispersion relations <ref>G. Peckham, Proc. Phys. Soc. 1967, 90, 657, doi: [http://iopscience.iop.org/article/10.1088/0370-1328/90/3/312/pdf 10.1088/0370-1328/90/3/312]</ref>.

By visualising the vibrational modes with the DLVisualize Animate feature, it is understood that the lowest two branches at Γ are overlapping acoustic transverse vibrations, and the second-lowest branch is the longitudinal vibration. Close to the Γ-symmetry k-point, the energy of acoustic phonons is approximately proportional to the wavenumber, <math> E = \frac{h \omega_0 k a}{2 \pi} </math>. The acoustic transverse vibrations in two perpendicular planes overlap to a greater degree in this calculation than that in experimental literature. Degree-of-freedom arguments suggest that 3N-3 optical vibrations exist above the acoustic modes (N = number of atoms in a unit cell), which holds true since there are 3 optical branches in the phonon dispersion relation graph.

== Phonon Density-of-States Functions In MgO ==

=== Phonon Density of States ===
The Phonon density of states function outlines the number of vibrational states available at each frequency. For a unisolated system, the density distribution of states is continuous.

Analytically the density of states can be obtained by integrating the phonon energy dispersion at an infinitesimally small increment of k. When the branch is flat at the point of interest, there is a high density of vibrational state for the particular k-point, or for the particular energy since energy is a function of the wavenumber k.

When full analytical integration is not available, a numerical integration can be performed by GULP by calculating the vibrational modes at each k-point and summing all densities multiplied by the weight of the k-point. The weighting factor, for a regular grid, is the inverse of the number of grid points. As the number of points approach infinity, the summation resembles the true value.

In practical calculations, a finite number of k-points is sampled over the first Brillouin zone. A Monkhorst-Pack mesh is an unbiased method of choosing a set of k-points uniformly spread over the first Brillouin zone. The uniform mesh is suitable for a regular grid such as MgO, where there are points of special interest are spread evenly. The number of points is defined by shrinking factors in three dimensions and the larger the shrinking factors, the finer and more accurate the DOS function.

=== Results ===
{| class="wikitable"
!Shrinking Factor
!MgO Phonon DOS
!Number of Sampled k Points
!Noticeable Maxima (inverse cm)
|-
|1x1x1
|[[File:MgO Phonon DOS (1x1).png|frameless]]
|'''1:'''

''1/2, 1/2, 1/2 (L)''

|256

381

688

805
|-
|2x2x2
|[[File:MgO Phonon DOS (2x2).png|frameless]]
|'''4:'''

''1/4, 1/4, 1/4''

''1/4, 1/4, 3/4;''

''1/4, 3/4, 3/4;''

''3/4, 1/4, 3/4''
|384

and four other individual peaks
|-
|4x4x4
|[[File:MgO Phonon DOS (4x4).png|frameless]]
|'''32'''

''1/8, 1/8, 1/8''

''1/8, 1/8, 3/8''

''1/8, 1/8, 5/8''

''1/8, 1/8, 7/8;''

''then repeat this sequence for all combinations of k in all directions''
|408

and eleven other individual peaks
|-
|8x8x8
|[[File:MgO Phonon DOS (8x8).png|frameless]]
|'''256'''
|432
|-
|16x16x16
|[[File:MgO Phonon DOS (16x16).png|frameless]]
|'''2048'''
|396
|-
|32x32x32
|[[File:MgO Phonon DOS (32x32).png|frameless]]
|'''16,384'''
|396
|}

The numbers of k points given in the table do not account for intrinsic symmetry of phonon modes in the Brillouin zone. Hence, the actual number of sampling points performed by GULP may be lower than those given in the table.

=== Analysis ===
As the shrinking factor increases, the DOS functions become more and more continuous. This is natural since the vibrational states should become more continuous for increasing size of the system. The shape of the DOS function finer than the 16x16x16 system agrees with calculations found in literature and experimental data<ref>A. R. Oganov, M. J. Gillan, G. D. Price, J. Chem. Phys. 2003, 118 (22), pp. 10174, doi: [http://uspex.stonybrook.edu/pdfs/JCP-2003.pdf 10.1063/1.1570394]</ref>. Maximum occurs around the 400 cm-1 frequency as expected because the phonon dispersion curves are the most 'crowded' at that particular frequency.

We can essentially observe the process of summation of the numerical integration as we increase the shrinking factor. For the 1x1x1 grid, only the L symmetry point is sampled, and the four peaks correspond to the frequencies at which the phonon dispersion branches intersect the L symmetry label. The two high-frequency peaks at 688 and 805 cm-1 correspond to the two points of intersection with the optical branches, and the lower-frequency peaks correspond to that with the acoustic branches. At L, the acoustic branches are flatter, giving a larger number of vibrational energies. As the number of sampling points increase, more and more points of intersection and frequencies are added to the DOS curve, giving a much smoother and more continuous curve.

=== Applicability To Other Crystal Lattices ===
'''ASK – Polarizability O2-, lattice parameters'''

CaO<ref>B. B. Karki, R.M. Wentzcovitch, Phys. Rev. B, 2003, 68 (22), pp. 224304, doi: [http://journals.aps.org/prb/pdf/10.1103/PhysRevB.68.224304 10.1103/PhysRevB.68.224304]</ref> Li<ref>S. Pal, Phys. Rev. B, 1970, 2 (12), pp. 4741, doi: [http://journals.aps.org/prb/pdf/10.1103/PhysRevB.2.4741 10.1103/PhysRevB.2.4741]</ref>

== The Helmholtz Free Energy of MgO ==

=== Quasi-Harmonic Approximation ===
Even at low vibrational states, the energy of a vibrational state deviates from perfect harmonic behavior. To introduce anharmonicity to the simulation, the quasi-harmonic approximation computes the positions and the bond lengths of individual atoms for vibrations corresponding to different wavenumbers. As the lattice compresses or expands, the zero-point energy (electrostatic) and the force constants vary. A new lattice is then computed at various k-points using the new force constant to minimise the overall free energy. The force constant calculation requires an accurate phonon frequency, hence a thorough sampling of k-points is required.

<math> A = U + \sum_i 0.5 \hbar \omega_i + k_B T \sum_i \ln \left( 1 - e^{-\frac{\hbar \omega_i}{k_B T}} \right)</math>, where A is the Helmholtz energy, U is the internal energy, kT is the Boltmann constant, h-bar is the reduced Plank constant, ω is the vibrational frequency, T is the temperature.

Since the optimisation of the new force constant are independent of other vibrations, the quasi-harmonic approximation breaks down at the high-temperature range because of significant phonon–phonon interactions. Another shortcoming of the quasi-harmonic approximation, due to its intrinsic harmonic nature even with different force constants, is the failure to model dissociation of the bonds at high internuclear distances and strong repulsions at high compressions, further invalidating this approach for high-temperature conditions.

=== Results and Analysis ===
{| class="wikitable"
!Shrinking
Factor
!Number of
Sampled k Points
!Zero-Point Energy
'''(eV)'''
!Helmholtz Free Energy
(eV)
!Heat Capacity At Constant Volume
(eV/K)
|-
|1x1x1
|1
|0.172063
|'''-40.93'''0301
|0.000350
|-
|2x2x2
|4
|0.174209
|'''-40.926'''609
|0.000347
|-
|4x4x4
|32
|0.174339
|'''-40.926'''450
|0.000347
|-
|8x8x8
|256
|0.174340
|'''-40.9264'''78
|0.000347
|-
|16x16x16
|2048
|0.174340
|'''-40.92648'''2
|0.000347
|-
|32x32x32
|16,384
|0.174340
|'''-40.92648'''4
|0.000347
|}
Overall the Helmholtz free energy is more positive as the shrinking factor is increased, which is also the trend of the zero-point energy. However a noticeable maxima occurs at the 2x2x2 system, where the energy peaks before dropping off to convergence. Because the zero-point energy trend lacks this feature, it is speculated that the peak in Helmholtz free energy is due to entropic contributions. As spatial sampling becomes finer, the TS term in Helmholtz free energy becomes increasingly negative (and more accurate). This is intuitive since a small shrinking factor does not account for the many modes of vibrations that contribute to overall entropy.

The 2x2x2 grid is suitable for calculations accurate to 0.5 meV per cell, which implies that it is also suitable for calculations accurate to 1 meV. For 0.1 meV accurate calculations, a 4x4x4 grid size is required.

=== Applicability To Other Crystal Lattices ===

== Thermal Expansion of MgO: Quasi-Harmonic Approximation ==

The coefficient of thermal expansion is a description of the change of an object's size with a change in temperature. More specifically, the expansion is described as a fractional change in size per increment of temperature. In this computational experiment, we are considered with the volume expansion. For a isotropic material such as the MgO crystal lattice, the thermal expansion coefficient of the object in 3D is three times the linear thermal expansion coefficient in 1D.

Mathematically the volume thermal expansion coefficient is defined as <math>\alpha_V = \frac{1}{V} \left( \frac{\partial V}{\partial T} \right)_p</math> , where αV is the volume thermal expansion coefficient at constant volume, V is the total volume of the system, T is the temperature of the system.

=== Computational Methodology ===
At each given temperature, an optimal configuration of the system is simulated by the method described in the previous section. The Helmholtz free energy is the system is minimised to a pre-set tolerance with respect to the cell volume. According to the variational principle, thermodynamic properties of a system is closest to the true value when the derivative of cell parameters with respect to temperature is zero.

=== Results and Analysis ===
The 8x8x8 lattice grid was used for thermal optimisations for free energy values that are accurate to 0.01 meV.
{| class="wikitable"
!Temperature
(K)
!Helmholtz Free Energy
(eV)
!Lattice Constant
(a=b=c, Angstroms)
!Angle
(α=β=γ, Degrees)
!Cell Parameter
Derivatives
|-
|0
|<nowiki>-40.901906</nowiki>
|2.986565
|59.999934
|0
|-
|100
|<nowiki>-40.902418</nowiki>
|2.986564
|59.999937
|0
|-
|200
|<nowiki>-40.909374</nowiki>
|2.987613
|59.999639
|0
|-
|300
|<nowiki>-40.928119</nowiki>
|2.989403
|59.999461
|<nowiki>-0.000002</nowiki>
|-
|400
|<nowiki>-40.958587</nowiki>
|2.989403
|59.992790
|0
|-
|500
|<nowiki>-40.999427</nowiki>
|2.994156
|59.999094
|0
|-
|600
|<nowiki>-41.049305</nowiki>
|2.996846
|59.998906
|0
|-
|700
|<nowiki>-41.049305</nowiki>
|2.996846
|59.998906
|0
|-
|800
|<nowiki>-41.107108</nowiki>
|2.999674
|59.998715
|0
|-
|900
|<nowiki>-41.171879</nowiki>
|3.002623
|59.998520
|0
|-
|1000
|<nowiki>-41.243004</nowiki>
|3.005674
|59.998323
|0
|-
|1200
|<nowiki>-41.488740</nowiki>
|3.015394
|59.997763
|0
|-
|1400
|<nowiki>-41.675515</nowiki>
|3.022433
|59.997451
|0.000001
|-
|1600
|<nowiki>-41.877981</nowiki>
|3.029988
|59.996775
|0
|-
|1800
|<nowiki>-42.094435</nowiki>
|3.038336
|59.995015
|0.000001
|-
|2000
|<nowiki>-42.323668</nowiki>
|3.046923
|59.994762
|<nowiki>-0.000002</nowiki>
|-
|2200
|<nowiki>-42.564747</nowiki>
|3.056841
|59.994783
|0
|-
|2400
|<nowiki>-42.817003</nowiki>
|3.068085
|59.992390
|<nowiki>-0.00002</nowiki>
|-
|2600
|<nowiki>-43.079967</nowiki>
|3.081449
|59.993314
|0.000002
|-
|2752
|<nowiki>-43.28700906</nowiki>
|3.094344
|59.994204
|0.000002
|-
|2852
|<nowiki>-43.42654086</nowiki>
|3.105499
|59.994105
|<nowiki>-0.000001</nowiki>
|-
|2952
|<nowiki>-43.56951767</nowiki>
|3.124113
|59.994437
|<nowiki>-0.000002</nowiki>
|}
The lattice constant increases linearly only in the range 300 - 1000 K.[[File:Screenshot-Thermal Expansion.xls - LibreOffice Calc.png|none|frame]]
The thermal expansion constant of MgO in this temperature range, assuming temperature independence, is equal to <math>8.74 \times 10^{-6} \text{K}^{-1}</math> with <math>R^2 = 0.99</math>.

[[File:Screenshot-Thermal Expansion.xls - LibreOffice Calc ALL.png|none|frame|Figure 3: Thermal expansion of MgO in all simulated temperatures.]]

If the whole temperature range 0K to 2952K was to be included, the thermal expansion of MgO is non-linear with respect to temperature. A polynomial fit to the fourth order has been determined with <math> R^2 = 1.00 </math> such that <math>\frac{a}{a_{0}} = 2.27 \times 10^{-5} T^4 - 1.04 \times 10^{-11} T^3 + 1.64 \times 10^{-8} T^2 - 8.35 \times 10^{-7} T - 2.60 \times 10^{-5}</math>, where <math>\frac{a}{a_{0}}</math> is the fractional change in the cell parameter relative to the cell parameter at T = 0 and <math> T </math> is the temperature in Kelvins. As the coefficient of thermal expansion is defined as the partial derivative of the fractional change with respect to temperature, the coefficient of thermal expansion can be expressed as <math> \alpha (T) = 9.08 \times 10^{-5} T^3 - 4.16 \times 10^{-11} T^2 + 3.28 \times 10^{-8} T - 5.85 \times 10^{-6} </math>.

== Thermal Expansion of MgO: Molecular Dynamics==

=== Computational Methodology ===
The molecular dynamics technique studies the trajectories of individual atoms by numerically solving Newton's equations of motion for an N-body system of interacting particles, where the forces of the particles are calculated by molecular mechanical force fields such as the Lennard-Jones or the Buckingham potential. It is an ''ab initio'' method of simulation applied to a large number of unit cells in order to obtain an accurate picture of the complex system in which analytical solutions are either inhibitively time-consuming or simply unavailable.

Thermodynamic properties of the isothermal-isobaric (NPT) ensemble are time-averaged over the entire space. For ergodic systems, the time averages should be identical to spatial averages provided by the quasi-harmonic approach. Due to the computational expense of this method, the system must be designed to account for available computational power by tailoring the simulation to a suitable timestep, time duration and scale. The scale of such a system must be large enough to eliminate periodic boundary condition artifacts, which in this case is suitable at the femtoseconds timestep and the 32x32x32 scale.

=== Results and Analysis ===
{| class="wikitable"
!Temperature
(K)
!Averaged Final Volume
(Cubic Angstroms)
!Volume Per Unit Cell
(Cubic Angstroms)
|-
|300
|602.917248
|18.8411640
|-
|400
|604.188880
|18.8809025
|-
|500
|605.968126
|18.9365039
|-
|600
|608.794097
|19.0248155
|-
|700
|608.745756
|19.0233049
|-
|800
|611.617711
|19.1130535
|-
|900
|611.936388
|19.1230121
|-
|1000
|615.0449322
|19.2347085
|-
|1100
|619.298570
|19.3530803
|-
|1200
|620.020014
|19.3756254
|-
|1300
|621.486358
|19.4214487
|-
|1400
|622.667288
|19.4583528
|-
|1600
|626.172143
|19.5678795
|-
|1800
|630.981409
|19.7181690
|-
|2000
|632.410325
|19.7628227
|-
|2400
|642.622061
|20.0819394
|-
|2600
|648.508557
|20.2658924
|-
|2800
|655.021305
|20.6318792
|-
|3000
|660.220134
|20.6318792
|}
[[File:Thermal Expansion MD v QHA.png|none|frame|Figure 4: A graph showing expansion of unit cell with increasing temperature, where the blue series is simulated by quasi-harmonic approximation, and the orange series is simulated by molecular dynamics. ]]The molecular dynamics data suggest a linear relationship between unit cell volume and temperature, and the thermal expansion coefficient, assuming its temperature independence is \alpha_V = 3.4599078 \times 10^{-5} \text{K}^{-1}.

Due to time constraints, I was unable to perform a molecular dynamics simulation at 0K. Extrapolating from obtained data, the volume per unit cell at absolute zero is 18.594625 ‎Å3, whereas with quasi-harmonic approximation the 0 K volume is 17.302424 Å3. It is evident that the quasi-harmonic approximation, with its inaccurate parabolic potential curve and boundary condition artifacts, significantly underestimate the magnitude of repulsive forces in a crystal lattice, giving rise to largely inaccurate high-temperature properties.

== Conclusion ==
'''[INSERT COMPARISON MD V QHA V EXPERIMENTAL'''<ref name=":0" /><ref>C. J. Engberg, E. H. Zehms, J. Am. Ceram. Soc., 1959, 42 (6), pp. 300, doi: [http://onlinelibrary.wiley.com/store/10.1111/j.1151-2916.1959.tb12958.x/asset/j.1151-2916.1959.tb12958.x.pdf?v=1&t=ikfzf86u&s=3040cd204ec490ffeaa1cabf505c7097de78aed2&systemMessage=Wiley+Online+Library+will+be+unavailable+on+Saturday+27th+February+from+09%3A00-14%3A00+GMT+%2F+04%3A00-09%3A00+EST+%2F+17%3A00-22%3A00+SGT+for+essential+maintenance.++Apologies+for+the+inconvenience. 10.1111/j.1151-2916.1959.tb12958.x]</ref>''']'''

Both molecular dynamics and quasi-harmonic simulations display a linear thermal expansion in the temperature range 0 K – 1200 K. The quasi-harmonic method is in better agreement with experimental data in the linear range, after which the thermal expansion diverges. In the high-temperature range, thermal expansion no longer obeys linearity according to quasi-harmonic simulations as molecules are no longer bound to harmonic force fields due to increased energies. The quasi-harmonic calculation does not simulate bond dissociation as a result of large energies in the high-temperature range, and fails to account for anharmonicity experienced at high bond distances.

On the other hand, molecular dynamics indicates linearity over a much larger range of temperature. At lower temperatures where periodic condition artifacts are more prominent, molecular dynamics simulations contain some inaccuracies within its somewhat simplistic classical force models. At high temperatures, molecular dyanmics is more reliable in calculating phonon interactions, and therefore more suited to high-temperature simulations.

For both methods, the calculations break down at and near the melting point of MgO. In an NPT system, the thermal expansion coefficient and the constant-pressure heat capacity are both proportional to the rate of change in vibrational energy with respect to temperature. We expect the heat capacity to instantaneously become infinite at the transition temperature. However it is not the case as both models predict constantly rising thermal expansion coefficients beyond the transition temperature.

Rep:Mod:pitsnake

2016-02-11T13:53:52Z

Jc6613: /* Results and Analysis */

This computational experiment aims to investigate thermodynamic properties of the MgO crystal lattice, described as either a rhombohedron primitive unit cell with a 60˚ internal angle, or a conventional face-centred cubic cell with a 90˚ internal angle. Different approaches to computational simulation are employed, and a summary of the thermodynamic properties and a comparison across the two different approaches are provided.

== Initial Calculation on MgO ==
An initial calculation on a single unit cell of MgO at absolute zero yields a Helmholtz free energy of –41.07531835 eV, or the equivalent of -3957.689074 kJ mol-1. The lattice parameter, a = b = c, of a primitive rhombohedric cell is 2.9783 Å, or the equivalent of 2.5792 Å in a conventional face-centred cubic cell (NaCl structure). The lattice parameter is in good agreeement with the LCAO-HF calculation (linear combination of atomic orbitals at the Hatree-Fock level) found in literature (2.573 Å <ref name=":0">''II-VI and I-VII Compounds; Semimagnetic Compounds'', Springer Berlin Heidelberg, Berlin, 1999, 41B, pp. 1-6, doi: [http://link.springer.com/chapter/10.1007%2F10681719_206 10.1007/10681719_206]</ref>), suggesting that the GULP shell-model simulation is in good agreement with quantum mechanical calculations.

== The Phonon Dispersion Relation In MgO ==

=== Phonons ===
A phonon is an excitation of a group of atoms in an elastic arrangement in a crystal lattice. While seen as a quasiparticle, a phonon represents an excited state in a quantised quantum mechanical mode of vibration. Understanding the behavioral changes of phonons at different temperatures is crucial in the study of thermal expansion, thermal conductivity and electrical conductivity of crystals.

The displacement of one atom in a crystal lattice leads to a propagating wave of the lattice, where the amplitude of the wave is the displacement of the atom from its equilibrium position, and the wavelength (λ) is the distance between consecutive points of the same amplitude (phase). The wavenumber is more useful in lattice representations, with <math>k = \frac{2 \pi}{\lambda}</math>.

=== Results and Analysis ===
Using GULP, phonon frequencies were calculated at various k-points along the paths Γ→X→W→K and Γ→L→W within the first Brillouin zone. A total of 50 k-points were sampled to give the graph shown below.

'''[INSERT Phonon Dispersion]'''

The shape of the computed phonon dispersion relation is in good agreement with experimental dispersion relations <ref>G. Peckham, Proc. Phys. Soc. 1967, 90, 657, doi: [http://iopscience.iop.org/article/10.1088/0370-1328/90/3/312/pdf 10.1088/0370-1328/90/3/312]</ref>.

By visualising the vibrational modes with the DLVisualize Animate feature, it is understood that the lowest two branches at Γ are overlapping acoustic transverse vibrations, and the second-lowest branch is the longitudinal vibration. Close to the Γ-symmetry k-point, the energy of acoustic phonons is approximately proportional to the wavenumber, <math> E = \frac{h \omega_0 k a}{2 \pi} </math>. The acoustic transverse vibrations in two perpendicular planes overlap to a greater degree in this calculation than that in experimental literature. Degree-of-freedom arguments suggest that 3N-3 optical vibrations exist above the acoustic modes (N = number of atoms in a unit cell), which holds true since there are 3 optical branches in the phonon dispersion relation graph.

== Phonon Density-of-States Functions In MgO ==

=== Phonon Density of States ===
The Phonon density of states function outlines the number of vibrational states available at each frequency. For a unisolated system, the density distribution of states is continuous.

Analytically the density of states can be obtained by integrating the phonon energy dispersion at an infinitesimally small increment of k. When the branch is flat at the point of interest, there is a high density of vibrational state for the particular k-point, or for the particular energy since energy is a function of the wavenumber k.

When full analytical integration is not available, a numerical integration can be performed by GULP by calculating the vibrational modes at each k-point and summing all densities multiplied by the weight of the k-point. The weighting factor, for a regular grid, is the inverse of the number of grid points. As the number of points approach infinity, the summation resembles the true value.

In practical calculations, a finite number of k-points is sampled over the first Brillouin zone. A Monkhorst-Pack mesh is an unbiased method of choosing a set of k-points uniformly spread over the first Brillouin zone. The uniform mesh is suitable for a regular grid such as MgO, where there are points of special interest are spread evenly. The number of points is defined by shrinking factors in three dimensions and the larger the shrinking factors, the finer and more accurate the DOS function.

=== Results ===
{| class="wikitable"
!Shrinking Factor
!MgO Phonon DOS
!Number of Sampled k Points
!Noticeable Maxima (inverse cm)
|-
|1x1x1
|[[File:MgO Phonon DOS (1x1).png|frameless]]
|'''1:'''

''1/2, 1/2, 1/2 (L)''

|256

381

688

805
|-
|2x2x2
|[[File:MgO Phonon DOS (2x2).png|frameless]]
|'''4:'''

''1/4, 1/4, 1/4''

''1/4, 1/4, 3/4;''

''1/4, 3/4, 3/4;''

''3/4, 1/4, 3/4''
|384

and four other individual peaks
|-
|4x4x4
|[[File:MgO Phonon DOS (4x4).png|frameless]]
|'''32'''

''1/8, 1/8, 1/8''

''1/8, 1/8, 3/8''

''1/8, 1/8, 5/8''

''1/8, 1/8, 7/8;''

''then repeat this sequence for all combinations of k in all directions''
|408

and eleven other individual peaks
|-
|8x8x8
|[[File:MgO Phonon DOS (8x8).png|frameless]]
|'''256'''
|432
|-
|16x16x16
|[[File:MgO Phonon DOS (16x16).png|frameless]]
|'''2048'''
|396
|-
|32x32x32
|[[File:MgO Phonon DOS (32x32).png|frameless]]
|'''16,384'''
|396
|}

The numbers of k points given in the table do not account for intrinsic symmetry of phonon modes in the Brillouin zone. Hence, the actual number of sampling points performed by GULP may be lower than those given in the table.

=== Analysis ===
As the shrinking factor increases, the DOS functions become more and more continuous. This is natural since the vibrational states should become more continuous for increasing size of the system. The shape of the DOS function finer than the 16x16x16 system agrees with calculations found in literature and experimental data<ref>A. R. Oganov, M. J. Gillan, G. D. Price, J. Chem. Phys. 2003, 118 (22), pp. 10174, doi: [http://uspex.stonybrook.edu/pdfs/JCP-2003.pdf 10.1063/1.1570394]</ref>. Maximum occurs around the 400 cm-1 frequency as expected because the phonon dispersion curves are the most 'crowded' at that particular frequency.

We can essentially observe the process of summation of the numerical integration as we increase the shrinking factor. For the 1x1x1 grid, only the L symmetry point is sampled, and the four peaks correspond to the frequencies at which the phonon dispersion branches intersect the L symmetry label. The two high-frequency peaks at 688 and 805 cm-1 correspond to the two points of intersection with the optical branches, and the lower-frequency peaks correspond to that with the acoustic branches. At L, the acoustic branches are flatter, giving a larger number of vibrational energies. As the number of sampling points increase, more and more points of intersection and frequencies are added to the DOS curve, giving a much smoother and more continuous curve.

=== Applicability To Other Crystal Lattices ===
'''ASK – Polarizability O2-, lattice parameters'''

CaO<ref>B. B. Karki, R.M. Wentzcovitch, Phys. Rev. B, 2003, 68 (22), pp. 224304, doi: [http://journals.aps.org/prb/pdf/10.1103/PhysRevB.68.224304 10.1103/PhysRevB.68.224304]</ref> Li<ref>S. Pal, Phys. Rev. B, 1970, 2 (12), pp. 4741, doi: [http://journals.aps.org/prb/pdf/10.1103/PhysRevB.2.4741 10.1103/PhysRevB.2.4741]</ref>

== The Helmholtz Free Energy of MgO ==

=== Quasi-Harmonic Approximation ===
Even at low vibrational states, the energy of a vibrational state deviates from perfect harmonic behavior. To introduce anharmonicity to the simulation, the quasi-harmonic approximation computes the positions and the bond lengths of individual atoms for vibrations corresponding to different wavenumbers. As the lattice compresses or expands, the zero-point energy (electrostatic) and the force constants vary. A new lattice is then computed at various k-points using the new force constant to minimise the overall free energy. The force constant calculation requires an accurate phonon frequency, hence a thorough sampling of k-points is required.

<math> A = U + \sum_i 0.5 \hbar \omega_i + k_B T \sum_i \ln \left( 1 - e^{-\frac{\hbar \omega_i}{k_B T}} \right)</math>, where A is the Helmholtz energy, U is the internal energy, kT is the Boltmann constant, h-bar is the reduced Plank constant, ω is the vibrational frequency, T is the temperature.

Since the optimisation of the new force constant are independent of other vibrations, the quasi-harmonic approximation breaks down at the high-temperature range because of significant phonon–phonon interactions. Another shortcoming of the quasi-harmonic approximation, due to its intrinsic harmonic nature even with different force constants, is the failure to model dissociation of the bonds at high internuclear distances and strong repulsions at high compressions, further invalidating this approach for high-temperature conditions.

=== Results and Analysis ===
{| class="wikitable"
!Shrinking
Factor
!Number of
Sampled k Points
!Zero-Point Energy
'''(eV)'''
!Helmholtz Free Energy
(eV)
!Heat Capacity At Constant Volume
(eV/K)
|-
|1x1x1
|1
|0.172063
|'''-40.93'''0301
|0.000350
|-
|2x2x2
|4
|0.174209
|'''-40.926'''609
|0.000347
|-
|4x4x4
|32
|0.174339
|'''-40.926'''450
|0.000347
|-
|8x8x8
|256
|0.174340
|'''-40.9264'''78
|0.000347
|-
|16x16x16
|2048
|0.174340
|'''-40.92648'''2
|0.000347
|-
|32x32x32
|16,384
|0.174340
|'''-40.92648'''4
|0.000347
|}
Overall the Helmholtz free energy is more positive as the shrinking factor is increased, which is also the trend of the zero-point energy. However a noticeable maxima occurs at the 2x2x2 system, where the energy peaks before dropping off to convergence. Because the zero-point energy trend lacks this feature, it is speculated that the peak in Helmholtz free energy is due to entropic contributions. As spatial sampling becomes finer, the TS term in Helmholtz free energy becomes increasingly negative (and more accurate). This is intuitive since a small shrinking factor does not account for the many modes of vibrations that contribute to overall entropy.

The 2x2x2 grid is suitable for calculations accurate to 0.5 meV per cell, which implies that it is also suitable for calculations accurate to 1 meV. For 0.1 meV accurate calculations, a 4x4x4 grid size is required.

=== Applicability To Other Crystal Lattices ===

== Thermal Expansion of MgO: Quasi-Harmonic Approximation ==

The coefficient of thermal expansion is a description of the change of an object's size with a change in temperature. More specifically, the expansion is described as a fractional change in size per increment of temperature. In this computational experiment, we are considered with the volume expansion. For a isotropic material such as the MgO crystal lattice, the thermal expansion coefficient of the object in 3D is three times the linear thermal expansion coefficient in 1D.

Mathematically the volume thermal expansion coefficient is defined as <math>\alpha_V = \frac{1}{V} \left( \frac{\partial V}{\partial T} \right)_p</math> , where αV is the volume thermal expansion coefficient at constant volume, V is the total volume of the system, T is the temperature of the system.

=== Computational Methodology ===
At each given temperature, an optimal configuration of the system is simulated by the method described in the previous section. The Helmholtz free energy is the system is minimised to a pre-set tolerance with respect to the cell volume. According to the variational principle, thermodynamic properties of a system is closest to the true value when the derivative of cell parameters with respect to temperature is zero.

=== Results and Analysis ===
The 8x8x8 lattice grid was used for thermal optimisations for free energy values that are accurate to 0.01 meV.
{| class="wikitable"
!Temperature
(K)
!Helmholtz Free Energy
(eV)
!Lattice Constant
(a=b=c, Angstroms)
!Angle
(α=β=γ, Degrees)
!Cell Parameter
Derivatives
|-
|0
|<nowiki>-40.901906</nowiki>
|2.986565
|59.999934
|0
|-
|100
|<nowiki>-40.902418</nowiki>
|2.986564
|59.999937
|0
|-
|200
|<nowiki>-40.909374</nowiki>
|2.987613
|59.999639
|0
|-
|300
|<nowiki>-40.928119</nowiki>
|2.989403
|59.999461
|<nowiki>-0.000002</nowiki>
|-
|400
|<nowiki>-40.958587</nowiki>
|2.989403
|59.992790
|0
|-
|500
|<nowiki>-40.999427</nowiki>
|2.994156
|59.999094
|0
|-
|600
|<nowiki>-41.049305</nowiki>
|2.996846
|59.998906
|0
|-
|700
|<nowiki>-41.049305</nowiki>
|2.996846
|59.998906
|0
|-
|800
|<nowiki>-41.107108</nowiki>
|2.999674
|59.998715
|0
|-
|900
|<nowiki>-41.171879</nowiki>
|3.002623
|59.998520
|0
|-
|1000
|<nowiki>-41.243004</nowiki>
|3.005674
|59.998323
|0
|-
|1200
|<nowiki>-41.488740</nowiki>
|3.015394
|59.997763
|0
|-
|1400
|<nowiki>-41.675515</nowiki>
|3.022433
|59.997451
|0.000001
|-
|1600
|<nowiki>-41.877981</nowiki>
|3.029988
|59.996775
|0
|-
|1800
|<nowiki>-42.094435</nowiki>
|3.038336
|59.995015
|0.000001
|-
|2000
|<nowiki>-42.323668</nowiki>
|3.046923
|59.994762
|<nowiki>-0.000002</nowiki>
|-
|2200
|<nowiki>-42.564747</nowiki>
|3.056841
|59.994783
|0
|-
|2400
|<nowiki>-42.817003</nowiki>
|3.068085
|59.992390
|<nowiki>-0.00002</nowiki>
|-
|2600
|<nowiki>-43.079967</nowiki>
|3.081449
|59.993314
|0.000002
|-
|2752
|<nowiki>-43.28700906</nowiki>
|3.094344
|59.994204
|0.000002
|-
|2852
|<nowiki>-43.42654086</nowiki>
|3.105499
|59.994105
|<nowiki>-0.000001</nowiki>
|-
|2952
|<nowiki>-43.56951767</nowiki>
|3.124113
|59.994437
|<nowiki>-0.000002</nowiki>
|}
The lattice constant increases linearly only in the range 300 - 1000 K.[[File:Screenshot-Thermal Expansion.xls - LibreOffice Calc.png|none|frame]]
The thermal expansion constant of MgO in this temperature range, assuming temperature independence, is equal to <math>8.74 \times 10^{-6} \text{K}^{-1}</math> with <math>R^2 = 0.99</math>.

[[File:Screenshot-Thermal Expansion.xls - LibreOffice Calc ALL.png|none|frame|Figure 3: Thermal expansion of MgO in all simulated temperatures.]]

If the whole temperature range 0K to 2952K was to be included, the thermal expansion of MgO is non-linear with respect to temperature. A polynomial fit to the fourth order has been determined with <math> R^2 = 1.00 </math> such that <math>\frac{a}{a_{0}} = 2.27 \times 10^{-5} T^4 - 1.04 \times 10^{-11} T^3 + 1.64 \times 10^{-8} T^2 - 8.35 \times 10^{-7} T - 2.60 \times 10^{-5}</math>, where <math>\frac{a}{a_{0}}</math> is the fractional change in the cell parameter relative to the cell parameter at T = 0 and <math> T </math> is the temperature in Kelvins. As the coefficient of thermal expansion is defined as the partial derivative of the fractional change with respect to temperature, the coefficient of thermal expansion can be expressed as <math> \alpha (T) = 9.08 \times 10^{-5} T^3 - 4.16 \times 10^{-11} T^2 + 3.28 \times 10^{-8} T - 5.85 \times 10^{-6} </math>.

== Thermal Expansion of MgO: Molecular Dynamics==

=== Computational Methodology ===
The molecular dynamics technique studies the trajectories of individual atoms by numerically solving Newton's equations of motion for an N-body system of interacting particles, where the forces of the particles are calculated by molecular mechanical force fields such as the Lennard-Jones or the Buckingham potential. It is an ''ab initio'' method of simulation applied to a large number of unit cells in order to obtain an accurate picture of the complex system in which analytical solutions are either inhibitively time-consuming or simply unavailable.

Thermodynamic properties of the isothermal-isobaric (NPT) ensemble are time-averaged over the entire space. For ergodic systems, the time averages should be identical to spatial averages provided by the quasi-harmonic approach. Due to the computational expense of this method, the system must be designed to account for available computational power by tailoring the simulation to a suitable timestep, time duration and scale. The scale of such a system must be large enough to eliminate periodic boundary condition artifacts, which in this case is suitable at the femtoseconds timestep and the 32x32x32 scale.

=== Results and Analysis ===
{| class="wikitable"
!Temperature
(K)
!Averaged Final Volume
(Cubic Angstroms)
!Volume Per Unit Cell
(Cubic Angstroms)
|-
|300
|602.917248
|18.8411640
|-
|400
|604.188880
|18.8809025
|-
|500
|605.968126
|18.9365039
|-
|600
|608.794097
|19.0248155
|-
|700
|608.745756
|19.0233049
|-
|800
|611.617711
|19.1130535
|-
|900
|611.936388
|19.1230121
|-
|1000
|615.0449322
|19.2347085
|-
|1100
|619.298570
|19.3530803
|-
|1200
|620.020014
|19.3756254
|-
|1300
|621.486358
|19.4214487
|-
|1400
|622.667288
|19.4583528
|-
|1600
|626.172143
|19.5678795
|-
|1800
|630.981409
|19.7181690
|-
|2000
|632.410325
|19.7628227
|-
|2400
|642.622061
|20.0819394
|-
|2600
|648.508557
|20.2658924
|-
|2800
|655.021305
|20.6318792
|-
|3000
|660.220134
|20.6318792
|}
[[File:Thermal Expansion MD v QHA.png|none|frame|Figure 4: A graph showing expansion of unit cell with increasing temperature, where the blue series is simulated by quasi-harmonic approximation, and the orange series is simulated by molecular dynamics. ]]The molecular dynamics data suggest a linear relationship between unit cell volume and temperature, and the thermal expansion coefficient, assuming its temperature independence is \alpha_V = 3.4599078 \times 10^{-5} \text{K}^{-1}.

== Conclusion ==
'''[INSERT COMPARISON MD V QHA V EXPERIMENTAL'''<ref name=":0" /><ref>C. J. Engberg, E. H. Zehms, J. Am. Ceram. Soc., 1959, 42 (6), pp. 300, doi: [http://onlinelibrary.wiley.com/store/10.1111/j.1151-2916.1959.tb12958.x/asset/j.1151-2916.1959.tb12958.x.pdf?v=1&t=ikfzf86u&s=3040cd204ec490ffeaa1cabf505c7097de78aed2&systemMessage=Wiley+Online+Library+will+be+unavailable+on+Saturday+27th+February+from+09%3A00-14%3A00+GMT+%2F+04%3A00-09%3A00+EST+%2F+17%3A00-22%3A00+SGT+for+essential+maintenance.++Apologies+for+the+inconvenience. 10.1111/j.1151-2916.1959.tb12958.x]</ref>''']'''

Both molecular dynamics and quasi-harmonic simulations display a linear thermal expansion in the temperature range 0 K – 1200 K. The quasi-harmonic method is in better agreement with experimental data in the linear range, after which the thermal expansion diverges. In the high-temperature range, thermal expansion no longer obeys linearity according to quasi-harmonic simulations as molecules are no longer bound to harmonic force fields due to increased energies. The quasi-harmonic calculation does not simulate bond dissociation as a result of large energies in the high-temperature range, and fails to account for anharmonicity experienced at high bond distances.

On the other hand, molecular dynamics indicates linearity over a much larger range of temperature. At lower temperatures where periodic condition artifacts are more prominent, molecular dynamics simulations contain some inaccuracies within its somewhat simplistic classical force models. At high temperatures, molecular dyanmics is more reliable in calculating phonon interactions, and therefore more suited to high-temperature simulations.

For both methods, the calculations break down at and near the melting point of MgO. In an NPT system, the thermal expansion coefficient and the constant-pressure heat capacity are both proportional to the rate of change in vibrational energy with respect to temperature. We expect the heat capacity to instantaneously become infinite at the transition temperature. However it is not the case as both models predict constantly rising thermal expansion coefficients beyond the transition temperature.

Rep:Mod:pitsnake

2016-02-11T13:47:35Z

Jc6613: /* Analysis */

This computational experiment aims to investigate thermodynamic properties of the MgO crystal lattice, described as either a rhombohedron primitive unit cell with a 60˚ internal angle, or a conventional face-centred cubic cell with a 90˚ internal angle. Different approaches to computational simulation are employed, and a summary of the thermodynamic properties and a comparison across the two different approaches are provided.

== Initial Calculation on MgO ==
An initial calculation on a single unit cell of MgO at absolute zero yields a Helmholtz free energy of –41.07531835 eV, or the equivalent of -3957.689074 kJ mol-1. The lattice parameter, a = b = c, of a primitive rhombohedric cell is 2.9783 Å, or the equivalent of 2.5792 Å in a conventional face-centred cubic cell (NaCl structure). The lattice parameter is in good agreeement with the LCAO-HF calculation (linear combination of atomic orbitals at the Hatree-Fock level) found in literature (2.573 Å <ref name=":0">''II-VI and I-VII Compounds; Semimagnetic Compounds'', Springer Berlin Heidelberg, Berlin, 1999, 41B, pp. 1-6, doi: [http://link.springer.com/chapter/10.1007%2F10681719_206 10.1007/10681719_206]</ref>), suggesting that the GULP shell-model simulation is in good agreement with quantum mechanical calculations.

== The Phonon Dispersion Relation In MgO ==

=== Phonons ===
A phonon is an excitation of a group of atoms in an elastic arrangement in a crystal lattice. While seen as a quasiparticle, a phonon represents an excited state in a quantised quantum mechanical mode of vibration. Understanding the behavioral changes of phonons at different temperatures is crucial in the study of thermal expansion, thermal conductivity and electrical conductivity of crystals.

The displacement of one atom in a crystal lattice leads to a propagating wave of the lattice, where the amplitude of the wave is the displacement of the atom from its equilibrium position, and the wavelength (λ) is the distance between consecutive points of the same amplitude (phase). The wavenumber is more useful in lattice representations, with <math>k = \frac{2 \pi}{\lambda}</math>.

=== Results and Analysis ===
Using GULP, phonon frequencies were calculated at various k-points along the paths Γ→X→W→K and Γ→L→W within the first Brillouin zone. A total of 50 k-points were sampled to give the graph shown below.

'''[INSERT Phonon Dispersion]'''

The shape of the computed phonon dispersion relation is in good agreement with experimental dispersion relations <ref>G. Peckham, Proc. Phys. Soc. 1967, 90, 657, doi: [http://iopscience.iop.org/article/10.1088/0370-1328/90/3/312/pdf 10.1088/0370-1328/90/3/312]</ref>.

By visualising the vibrational modes with the DLVisualize Animate feature, it is understood that the lowest two branches at Γ are overlapping acoustic transverse vibrations, and the second-lowest branch is the longitudinal vibration. Close to the Γ-symmetry k-point, the energy of acoustic phonons is approximately proportional to the wavenumber, <math> E = \frac{h \omega_0 k a}{2 \pi} </math>. The acoustic transverse vibrations in two perpendicular planes overlap to a greater degree in this calculation than that in experimental literature. Degree-of-freedom arguments suggest that 3N-3 optical vibrations exist above the acoustic modes (N = number of atoms in a unit cell), which holds true since there are 3 optical branches in the phonon dispersion relation graph.

== Phonon Density-of-States Functions In MgO ==

=== Phonon Density of States ===
The Phonon density of states function outlines the number of vibrational states available at each frequency. For a unisolated system, the density distribution of states is continuous.

Analytically the density of states can be obtained by integrating the phonon energy dispersion at an infinitesimally small increment of k. When the branch is flat at the point of interest, there is a high density of vibrational state for the particular k-point, or for the particular energy since energy is a function of the wavenumber k.

When full analytical integration is not available, a numerical integration can be performed by GULP by calculating the vibrational modes at each k-point and summing all densities multiplied by the weight of the k-point. The weighting factor, for a regular grid, is the inverse of the number of grid points. As the number of points approach infinity, the summation resembles the true value.

In practical calculations, a finite number of k-points is sampled over the first Brillouin zone. A Monkhorst-Pack mesh is an unbiased method of choosing a set of k-points uniformly spread over the first Brillouin zone. The uniform mesh is suitable for a regular grid such as MgO, where there are points of special interest are spread evenly. The number of points is defined by shrinking factors in three dimensions and the larger the shrinking factors, the finer and more accurate the DOS function.

=== Results ===
{| class="wikitable"
!Shrinking Factor
!MgO Phonon DOS
!Number of Sampled k Points
!Noticeable Maxima (inverse cm)
|-
|1x1x1
|[[File:MgO Phonon DOS (1x1).png|frameless]]
|'''1:'''

''1/2, 1/2, 1/2 (L)''

|256

381

688

805
|-
|2x2x2
|[[File:MgO Phonon DOS (2x2).png|frameless]]
|'''4:'''

''1/4, 1/4, 1/4''

''1/4, 1/4, 3/4;''

''1/4, 3/4, 3/4;''

''3/4, 1/4, 3/4''
|384

and four other individual peaks
|-
|4x4x4
|[[File:MgO Phonon DOS (4x4).png|frameless]]
|'''32'''

''1/8, 1/8, 1/8''

''1/8, 1/8, 3/8''

''1/8, 1/8, 5/8''

''1/8, 1/8, 7/8;''

''then repeat this sequence for all combinations of k in all directions''
|408

and eleven other individual peaks
|-
|8x8x8
|[[File:MgO Phonon DOS (8x8).png|frameless]]
|'''256'''
|432
|-
|16x16x16
|[[File:MgO Phonon DOS (16x16).png|frameless]]
|'''2048'''
|396
|-
|32x32x32
|[[File:MgO Phonon DOS (32x32).png|frameless]]
|'''16,384'''
|396
|}

The numbers of k points given in the table do not account for intrinsic symmetry of phonon modes in the Brillouin zone. Hence, the actual number of sampling points performed by GULP may be lower than those given in the table.

=== Analysis ===
As the shrinking factor increases, the DOS functions become more and more continuous. This is natural since the vibrational states should become more continuous for increasing size of the system. The shape of the DOS function finer than the 16x16x16 system agrees with calculations found in literature and experimental data<ref>A. R. Oganov, M. J. Gillan, G. D. Price, J. Chem. Phys. 2003, 118 (22), pp. 10174, doi: [http://uspex.stonybrook.edu/pdfs/JCP-2003.pdf 10.1063/1.1570394]</ref>. Maximum occurs around the 400 cm-1 frequency as expected because the phonon dispersion curves are the most 'crowded' at that particular frequency.

We can essentially observe the process of summation of the numerical integration as we increase the shrinking factor. For the 1x1x1 grid, only the L symmetry point is sampled, and the four peaks correspond to the frequencies at which the phonon dispersion branches intersect the L symmetry label. The two high-frequency peaks at 688 and 805 cm-1 correspond to the two points of intersection with the optical branches, and the lower-frequency peaks correspond to that with the acoustic branches. At L, the acoustic branches are flatter, giving a larger number of vibrational energies. As the number of sampling points increase, more and more points of intersection and frequencies are added to the DOS curve, giving a much smoother and more continuous curve.

=== Applicability To Other Crystal Lattices ===
'''ASK – Polarizability O2-, lattice parameters'''

CaO<ref>B. B. Karki, R.M. Wentzcovitch, Phys. Rev. B, 2003, 68 (22), pp. 224304, doi: [http://journals.aps.org/prb/pdf/10.1103/PhysRevB.68.224304 10.1103/PhysRevB.68.224304]</ref> Li<ref>S. Pal, Phys. Rev. B, 1970, 2 (12), pp. 4741, doi: [http://journals.aps.org/prb/pdf/10.1103/PhysRevB.2.4741 10.1103/PhysRevB.2.4741]</ref>

== The Helmholtz Free Energy of MgO ==

=== Quasi-Harmonic Approximation ===
Even at low vibrational states, the energy of a vibrational state deviates from perfect harmonic behavior. To introduce anharmonicity to the simulation, the quasi-harmonic approximation computes the positions and the bond lengths of individual atoms for vibrations corresponding to different wavenumbers. As the lattice compresses or expands, the zero-point energy (electrostatic) and the force constants vary. A new lattice is then computed at various k-points using the new force constant to minimise the overall free energy. The force constant calculation requires an accurate phonon frequency, hence a thorough sampling of k-points is required.

<math> A = U + \sum_i 0.5 \hbar \omega_i + k_B T \sum_i \ln \left( 1 - e^{-\frac{\hbar \omega_i}{k_B T}} \right)</math>, where A is the Helmholtz energy, U is the internal energy, kT is the Boltmann constant, h-bar is the reduced Plank constant, ω is the vibrational frequency, T is the temperature.

Since the optimisation of the new force constant are independent of other vibrations, the quasi-harmonic approximation breaks down at the high-temperature range because of significant phonon–phonon interactions. Another shortcoming of the quasi-harmonic approximation, due to its intrinsic harmonic nature even with different force constants, is the failure to model dissociation of the bonds at high internuclear distances and strong repulsions at high compressions, further invalidating this approach for high-temperature conditions.

=== Results and Analysis ===
{| class="wikitable"
!Shrinking
Factor
!Number of
Sampled k Points
!Zero-Point Energy
'''(eV)'''
!Helmholtz Free Energy
(eV)
!Heat Capacity At Constant Volume
(eV/K)
|-
|1x1x1
|1
|0.172063
|'''-40.93'''0301
|0.000350
|-
|2x2x2
|4
|0.174209
|'''-40.926'''609
|0.000347
|-
|4x4x4
|32
|0.174339
|'''-40.926'''450
|0.000347
|-
|8x8x8
|256
|0.174340
|'''-40.9264'''78
|0.000347
|-
|16x16x16
|2048
|0.174340
|'''-40.92648'''2
|0.000347
|-
|32x32x32
|16,384
|0.174340
|'''-40.92648'''4
|0.000347
|}
Overall the Helmholtz free energy is more positive as the shrinking factor is increased, which is also the trend of the zero-point energy. However a noticeable maxima occurs at the 2x2x2 system, where the energy peaks before dropping off to convergence. Because the zero-point energy trend lacks this feature, it is speculated that the peak in Helmholtz free energy is due to entropic contributions. As spatial sampling becomes finer, the TS term in Helmholtz free energy becomes increasingly negative (and more accurate). This is intuitive since a small shrinking factor does not account for the many modes of vibrations that contribute to overall entropy.

The 2x2x2 grid is suitable for calculations accurate to 0.5 meV per cell, which implies that it is also suitable for calculations accurate to 1 meV. For 0.1 meV accurate calculations, a 4x4x4 grid size is required.

=== Applicability To Other Crystal Lattices ===

== Thermal Expansion of MgO: Quasi-Harmonic Approximation ==

The coefficient of thermal expansion is a description of the change of an object's size with a change in temperature. More specifically, the expansion is described as a fractional change in size per increment of temperature. In this computational experiment, we are considered with the volume expansion. For a isotropic material such as the MgO crystal lattice, the thermal expansion coefficient of the object in 3D is three times the linear thermal expansion coefficient in 1D.

Mathematically the volume thermal expansion coefficient is defined as <math>\alpha_V = \frac{1}{V} \left( \frac{\partial V}{\partial T} \right)_p</math> , where αV is the volume thermal expansion coefficient at constant volume, V is the total volume of the system, T is the temperature of the system.

=== Computational Methodology ===
At each given temperature, an optimal configuration of the system is simulated by the method described in the previous section. The Helmholtz free energy is the system is minimised to a pre-set tolerance with respect to the cell volume. According to the variational principle, thermodynamic properties of a system is closest to the true value when the derivative of cell parameters with respect to temperature is zero.

=== Results and Analysis ===
The 8x8x8 lattice grid was used for thermal optimisations for free energy values that are accurate to 0.01 meV.
{| class="wikitable"
!Temperature
(K)
!Helmholtz Free Energy
(eV)
!Lattice Constant
(a=b=c, Angstroms)
!Angle
(α=β=γ, Degrees)
!Cell Parameter
Derivatives
|-
|0
|<nowiki>-40.901906</nowiki>
|2.986565
|59.999934
|0
|-
|100
|<nowiki>-40.902418</nowiki>
|2.986564
|59.999937
|0
|-
|200
|<nowiki>-40.909374</nowiki>
|2.987613
|59.999639
|0
|-
|300
|<nowiki>-40.928119</nowiki>
|2.989403
|59.999461
|<nowiki>-0.000002</nowiki>
|-
|400
|<nowiki>-40.958587</nowiki>
|2.989403
|59.992790
|0
|-
|500
|<nowiki>-40.999427</nowiki>
|2.994156
|59.999094
|0
|-
|600
|<nowiki>-41.049305</nowiki>
|2.996846
|59.998906
|0
|-
|700
|<nowiki>-41.049305</nowiki>
|2.996846
|59.998906
|0
|-
|800
|<nowiki>-41.107108</nowiki>
|2.999674
|59.998715
|0
|-
|900
|<nowiki>-41.171879</nowiki>
|3.002623
|59.998520
|0
|-
|1000
|<nowiki>-41.243004</nowiki>
|3.005674
|59.998323
|0
|-
|1200
|<nowiki>-41.488740</nowiki>
|3.015394
|59.997763
|0
|-
|1400
|<nowiki>-41.675515</nowiki>
|3.022433
|59.997451
|0.000001
|-
|1600
|<nowiki>-41.877981</nowiki>
|3.029988
|59.996775
|0
|-
|1800
|<nowiki>-42.094435</nowiki>
|3.038336
|59.995015
|0.000001
|-
|2000
|<nowiki>-42.323668</nowiki>
|3.046923
|59.994762
|<nowiki>-0.000002</nowiki>
|-
|2200
|<nowiki>-42.564747</nowiki>
|3.056841
|59.994783
|0
|-
|2400
|<nowiki>-42.817003</nowiki>
|3.068085
|59.992390
|<nowiki>-0.00002</nowiki>
|-
|2600
|<nowiki>-43.079967</nowiki>
|3.081449
|59.993314
|0.000002
|-
|2752
|<nowiki>-43.28700906</nowiki>
|3.094344
|59.994204
|0.000002
|-
|2852
|<nowiki>-43.42654086</nowiki>
|3.105499
|59.994105
|<nowiki>-0.000001</nowiki>
|-
|2952
|<nowiki>-43.56951767</nowiki>
|3.124113
|59.994437
|<nowiki>-0.000002</nowiki>
|}
The lattice constant increases linearly only in the range 300 - 1000 K.[[File:Screenshot-Thermal Expansion.xls - LibreOffice Calc.png|none|frame]]
The thermal expansion constant of MgO in this temperature range, assuming temperature independence, is equal to <math>8.74 \times 10^{-6} \text{K}^{-1}</math> with <math>R^2 = 0.99</math>.

[[File:Screenshot-Thermal Expansion.xls - LibreOffice Calc ALL.png|none|frame|Figure 3: Thermal expansion of MgO in all simulated temperatures.]]

If the whole temperature range 0K to 2952K was to be included, the thermal expansion of MgO is non-linear with respect to temperature. A polynomial fit to the fourth order has been determined with <math> R^2 = 1.00 </math> such that <math>\frac{a}{a_{0}} = 2.27 \times 10^{-5} T^4 - 1.04 \times 10^{-11} T^3 + 1.64 \times 10^{-8} T^2 - 8.35 \times 10^{-7} T - 2.60 \times 10^{-5}</math>, where <math>\frac{a}{a_{0}}</math> is the fractional change in the cell parameter relative to the cell parameter at T = 0 and <math> T </math> is the temperature in Kelvins. As the coefficient of thermal expansion is defined as the partial derivative of the fractional change with respect to temperature, the coefficient of thermal expansion can be expressed as <math> \alpha (T) = 9.08 \times 10^{-5} T^3 - 4.16 \times 10^{-11} T^2 + 3.28 \times 10^{-8} T - 5.85 \times 10^{-6} </math>.

== Thermal Expansion of MgO: Molecular Dynamics==

=== Computational Methodology ===
The molecular dynamics technique studies the trajectories of individual atoms by numerically solving Newton's equations of motion for an N-body system of interacting particles, where the forces of the particles are calculated by molecular mechanical force fields such as the Lennard-Jones or the Buckingham potential. It is an ''ab initio'' method of simulation applied to a large number of unit cells in order to obtain an accurate picture of the complex system in which analytical solutions are either inhibitively time-consuming or simply unavailable.

Thermodynamic properties of the isothermal-isobaric (NPT) ensemble are time-averaged over the entire space. For ergodic systems, the time averages should be identical to spatial averages provided by the quasi-harmonic approach. Due to the computational expense of this method, the system must be designed to account for available computational power by tailoring the simulation to a suitable timestep, time duration and scale. The scale of such a system must be large enough to eliminate periodic boundary condition artifacts, which in this case is suitable at the femtoseconds timestep and the 32x32x32 scale.

=== Results and Analysis ===
{| class="wikitable"
!Temperature
(K)
!Averaged Final Volume
(Cubic Angstroms)
!Volume Per Unit Cell
(Cubic Angstroms)
|-
|300
|602.917248
|18.8411640
|-
|400
|604.188880
|18.8809025
|-
|500
|605.968126
|18.9365039
|-
|600
|608.794097
|19.0248155
|-
|700
|608.745756
|19.0233049
|-
|800
|611.617711
|19.1130535
|-
|900
|611.936388
|19.1230121
|-
|1000
|615.0449322
|19.2347085
|-
|1100
|619.298570
|19.3530803
|-
|1200
|620.020014
|19.3756254
|-
|1300
|621.486358
|19.4214487
|-
|1400
|622.667288
|19.4583528
|-
|1600
|626.172143
|19.5678795
|-
|1800
|630.981409
|19.7181690
|-
|2000
|632.410325
|19.7628227
|-
|2400
|642.622061
|20.0819394
|-
|2600
|648.508557
|20.2658924
|-
|2800
|655.021305
|20.6318792
|-
|3000
|660.220134
|20.6318792
|}
[[File:Thermal Expansion MD v QHA.png|none|frame|Figure 4: A graph showing expansion of unit cell with increasing temperature, where the blue series is simulated by quasi-harmonic approximation, and the orange series is simulated by molecular dynamics. ]]

== Conclusion ==
'''[INSERT COMPARISON MD V QHA V EXPERIMENTAL'''<ref name=":0" /><ref>C. J. Engberg, E. H. Zehms, J. Am. Ceram. Soc., 1959, 42 (6), pp. 300, doi: [http://onlinelibrary.wiley.com/store/10.1111/j.1151-2916.1959.tb12958.x/asset/j.1151-2916.1959.tb12958.x.pdf?v=1&t=ikfzf86u&s=3040cd204ec490ffeaa1cabf505c7097de78aed2&systemMessage=Wiley+Online+Library+will+be+unavailable+on+Saturday+27th+February+from+09%3A00-14%3A00+GMT+%2F+04%3A00-09%3A00+EST+%2F+17%3A00-22%3A00+SGT+for+essential+maintenance.++Apologies+for+the+inconvenience. 10.1111/j.1151-2916.1959.tb12958.x]</ref>''']'''

Both molecular dynamics and quasi-harmonic simulations display a linear thermal expansion in the temperature range 0 K – 1200 K. The quasi-harmonic method is in better agreement with experimental data in the linear range, after which the thermal expansion diverges. In the high-temperature range, thermal expansion no longer obeys linearity according to quasi-harmonic simulations as molecules are no longer bound to harmonic force fields due to increased energies. The quasi-harmonic calculation does not simulate bond dissociation as a result of large energies in the high-temperature range, and fails to account for anharmonicity experienced at high bond distances.

On the other hand, molecular dynamics indicates linearity over a much larger range of temperature. At lower temperatures where periodic condition artifacts are more prominent, molecular dynamics simulations contain some inaccuracies within its somewhat simplistic classical force models. At high temperatures, molecular dyanmics is more reliable in calculating phonon interactions, and therefore more suited to high-temperature simulations.

For both methods, the calculations break down at and near the melting point of MgO. In an NPT system, the thermal expansion coefficient and the constant-pressure heat capacity are both proportional to the rate of change in vibrational energy with respect to temperature. We expect the heat capacity to instantaneously become infinite at the transition temperature. However it is not the case as both models predict constantly rising thermal expansion coefficients beyond the transition temperature.

File:Thermal Expansion MD v QHA.png

2016-02-11T13:44:55Z

Jc6613:

2016-02-09T23:26:35Z

Jc6613: /* Thermal Expansion of MgO: Quasi-Harmonic Approximation */

This computational experiment aims to investigate thermodynamic properties of the MgO crystal lattice, described as either a rhombohedron primitive unit cell with a 60˚ internal angle, or a conventional face-centred cubic cell with a 90˚ internal angle. Different approaches to computational simulation are employed, and a summary of the thermodynamic properties and a comparison across the two different approaches are provided.

== Initial Calculation on MgO ==
An initial calculation on a single unit cell of MgO at absolute zero yields a Helmholtz free energy of –41.07531835 eV, or the equivalent of -3957.689074 kJ mol-1. The lattice parameter, a = b = c, of a primitive rhombohedric cell is 2.9783 Å, or the equivalent of 2.5792 Å in a conventional face-centred cubic cell (NaCl structure). The lattice parameter is in good agreeement with the LCAO-HF calculation (linear combination of atomic orbitals at the Hatree-Fock level) found in literature (2.573 Å <ref name=":0">''II-VI and I-VII Compounds; Semimagnetic Compounds'', Springer Berlin Heidelberg, Berlin, 1999, 41B, pp. 1-6, doi: [http://link.springer.com/chapter/10.1007%2F10681719_206 10.1007/10681719_206]</ref>), suggesting that the GULP shell-model simulation is in good agreement with quantum mechanical calculations.

== The Phonon Dispersion Relation In MgO ==

=== Phonons ===
A phonon is an excitation of a group of atoms in an elastic arrangement in a crystal lattice. While seen as a quasiparticle, a phonon represents an excited state in a quantised quantum mechanical mode of vibration. Understanding the behavioral changes of phonons at different temperatures is crucial in the study of thermal expansion, thermal conductivity and electrical conductivity of crystals.

The displacement of one atom in a crystal lattice leads to a propagating wave of the lattice, where the amplitude of the wave is the displacement of the atom from its equilibrium position, and the wavelength (λ) is the distance between consecutive points of the same amplitude (phase). The wavenumber is more useful in lattice representations, with <math>k = \frac{2 \pi}{\lambda}</math>.

=== Results and Analysis ===
Using GULP, phonon frequencies were calculated at various k-points along the paths Γ→X→W→K and Γ→L→W within the first Brillouin zone. A total of 50 k-points were sampled to give the graph shown below.

'''[INSERT Phonon Dispersion]'''

The shape of the computed phonon dispersion relation is in good agreement with experimental dispersion relations <ref>G. Peckham, Proc. Phys. Soc. 1967, 90, 657, doi: [http://iopscience.iop.org/article/10.1088/0370-1328/90/3/312/pdf 10.1088/0370-1328/90/3/312]</ref>.

By visualising the vibrational modes with the DLVisualize Animate feature, it is understood that the lowest two branches at Γ are overlapping acoustic transverse vibrations, and the second-lowest branch is the longitudinal vibration. Close to the Γ-symmetry k-point, the energy of acoustic phonons is approximately proportional to the wavenumber, <math> E = \frac{h \omega_0 k a}{2 \pi} </math>. The acoustic transverse vibrations in two perpendicular planes overlap to a greater degree in this calculation than that in experimental literature. Degree-of-freedom arguments suggest that 3N-3 optical vibrations exist above the acoustic modes (N = number of atoms in a unit cell), which holds true since there are 3 optical branches in the phonon dispersion relation graph.

== Phonon Density-of-States Functions In MgO ==

=== Phonon Density of States ===
The Phonon density of states function outlines the number of vibrational states available at each frequency. For a unisolated system, the density distribution of states is continuous.

Analytically the density of states can be obtained by integrating the phonon energy dispersion at an infinitesimally small increment of k. When the branch is flat at the point of interest, there is a high density of vibrational state for the particular k-point, or for the particular energy since energy is a function of the wavenumber k.

When full analytical integration is not available, a numerical integration can be performed by GULP by calculating the vibrational modes at each k-point and summing all densities multiplied by the weight of the k-point. The weighting factor, for a regular grid, is the inverse of the number of grid points. As the number of points approach infinity, the summation resembles the true value.

In practical calculations, a finite number of k-points is sampled over the first Brillouin zone. A Monkhorst-Pack mesh is an unbiased method of choosing a set of k-points uniformly spread over the first Brillouin zone. The uniform mesh is suitable for a regular grid such as MgO, where there are points of special interest are spread evenly. The number of points is defined by shrinking factors in three dimensions and the larger the shrinking factors, the finer and more accurate the DOS function.

=== Results ===
{| class="wikitable"
!Shrinking Factor
!MgO Phonon DOS
!Number of Sampled k Points
!Noticeable Maxima (inverse cm)
|-
|1x1x1
|[[File:MgO Phonon DOS (1x1).png|frameless]]
|'''1:'''

''1/2, 1/2, 1/2 (L)''

|256

381

688

805
|-
|2x2x2
|[[File:MgO Phonon DOS (2x2).png|frameless]]
|'''4:'''

''1/4, 1/4, 1/4''

''1/4, 1/4, 3/4;''

''1/4, 3/4, 3/4;''

''3/4, 1/4, 3/4''
|384

and four other individual peaks
|-
|4x4x4
|[[File:MgO Phonon DOS (4x4).png|frameless]]
|'''32'''

''1/8, 1/8, 1/8''

''1/8, 1/8, 3/8''

''1/8, 1/8, 5/8''

''1/8, 1/8, 7/8;''

''then repeat this sequence for all combinations of k in all directions''
|408

and eleven other individual peaks
|-
|8x8x8
|[[File:MgO Phonon DOS (8x8).png|frameless]]
|'''256'''
|432
|-
|16x16x16
|[[File:MgO Phonon DOS (16x16).png|frameless]]
|'''2048'''
|396
|-
|32x32x32
|[[File:MgO Phonon DOS (32x32).png|frameless]]
|'''16,384'''
|396
|}

The numbers of k points given in the table do not account for intrinsic symmetry of phonon modes in the Brillouin zone. Hence, the actual number of sampling points performed by GULP may be lower than those given in the table.

=== Analysis ===
As the shrinking factor increases, the DOS functions become more and more continuous. This is natural since the vibrational states should become more continuous for increasing size of the system. The shape of the DOS function finer than the 16x16x16 system agrees with calculations found in literature and experimental data<ref>A. R. Oganov, M. J. Gillan, G. D. Price, J. Chem. Phys. 2003, 118 (22), pp. 10174, doi: [http://uspex.stonybrook.edu/pdfs/JCP-2003.pdf 10.1063/1.1570394]</ref>. Maximum occurs around the 400 cm-1 frequency as expected because the phonon dispersion curves are the most 'crowded' at that particular frequency.

We can essentially observe the process of summation of the numerical integration as we increase the shrinking factor. For the 1x1x1 grid, only the L symmetry point is sampled, and the four peaks correspond to the frequencies at which the phonon dispersion branches intersect the L symmetry label. The two high-frequency peaks at 688 and 805 cm-1 correspond to the two points of intersection with the optical branches, and the lower-frequency peaks correspond to that with the acoustic branches. At L, the acoustic branches are flatter, giving a larger number of vibrational energies. As the number of sampling points increase, more and more points of intersection and frequencies are added to the DOS curve, giving a much smoother and more continuous curve.

=== Applicability To Other Crystal Lattices ===
'''ASK – Polarizability O2-, lattice parameters'''

CaO<ref>B. B. Karki, R.M. Wentzcovitch, Phys. Rev. B, 2003, 68 (22), pp. 224304, doi: [http://journals.aps.org/prb/pdf/10.1103/PhysRevB.68.224304 10.1103/PhysRevB.68.224304]</ref> Li<ref>S. Pal, Phys. Rev. B, 1970, 2 (12), pp. 4741, doi: [http://journals.aps.org/prb/pdf/10.1103/PhysRevB.2.4741 10.1103/PhysRevB.2.4741]</ref>

== The Helmholtz Free Energy of MgO ==

=== Quasi-Harmonic Approximation ===
Even at low vibrational states, the energy of a vibrational state deviates from perfect harmonic behavior. To introduce anharmonicity to the simulation, the quasi-harmonic approximation computes the positions and the bond lengths of individual atoms for vibrations corresponding to different wavenumbers. As the lattice compresses or expands, the zero-point energy (electrostatic) and the force constants vary. A new lattice is then computed at various k-points using the new force constant to minimise the overall free energy. The force constant calculation requires an accurate phonon frequency, hence a thorough sampling of k-points is required.

<math> A = U + \sum_i 0.5 \hbar \omega_i + k_B T \sum_i \ln \left( 1 - e^{-\frac{\hbar \omega_i}{k_B T}} \right)</math>, where A is the Helmholtz energy, U is the internal energy, kT is the Boltmann constant, h-bar is the reduced Plank constant, ω is the vibrational frequency, T is the temperature.

Since the optimisation of the new force constant are independent of other vibrations, the quasi-harmonic approximation breaks down at the high-temperature range because of significant phonon–phonon interactions. Another shortcoming of the quasi-harmonic approximation, due to its intrinsic harmonic nature even with different force constants, is the failure to model dissociation of the bonds at high internuclear distances and strong repulsions at high compressions, further invalidating this approach for high-temperature conditions.

=== Results and Analysis ===
{| class="wikitable"
!Shrinking
Factor
!Number of
Sampled k Points
!Zero-Point Energy
'''(eV)'''
!Helmholtz Free Energy
(eV)
!Heat Capacity At Constant Volume
(eV/K)
|-
|1x1x1
|1
|0.172063
|'''-40.93'''0301
|0.000350
|-
|2x2x2
|4
|0.174209
|'''-40.926'''609
|0.000347
|-
|4x4x4
|32
|0.174339
|'''-40.926'''450
|0.000347
|-
|8x8x8
|256
|0.174340
|'''-40.9264'''78
|0.000347
|-
|16x16x16
|2048
|0.174340
|'''-40.92648'''2
|0.000347
|-
|32x32x32
|16,384
|0.174340
|'''-40.92648'''4
|0.000347
|}
Overall the Helmholtz free energy is more positive as the shrinking factor is increased, which is also the trend of the zero-point energy. However a noticeable maxima occurs at the 2x2x2 system, where the energy peaks before dropping off to convergence. Because the zero-point energy trend lacks this feature, it is speculated that the peak in Helmholtz free energy is due to entropic contributions. As spatial sampling becomes finer, the TS term in Helmholtz free energy becomes increasingly negative (and more accurate). This is intuitive since a small shrinking factor does not account for the many modes of vibrations that contribute to overall entropy.

The 2x2x2 grid is suitable for calculations accurate to 0.5 meV per cell, which implies that it is also suitable for calculations accurate to 1 meV. For 0.1 meV accurate calculations, a 4x4x4 grid size is required.

== Thermal Expansion of MgO: Quasi-Harmonic Approximation ==

The coefficient of thermal expansion is a description of the change of an object's size with a change in temperature. More specifically, the expansion is described as a fractional change in size per increment of temperature. In this computational experiment, we are considered with the volume expansion. For a isotropic material such as the MgO crystal lattice, the thermal expansion coefficient of the object in 3D is three times the linear thermal expansion coefficient in 1D.

Mathematically the volume thermal expansion coefficient is defined as <math>\alpha_V = \frac{1}{V} \left( \frac{\partial V}{\partial T} \right)_p</math> , where αV is the volume thermal expansion coefficient at constant volume, V is the total volume of the system, T is the temperature of the system.

=== Computational Methodology ===
At each given temperature, an optimal configuration of the system is simulated by the method described in the previous section. The Helmholtz free energy is the system is minimised to a pre-set tolerance with respect to the cell volume. According to the variational principle, thermodynamic properties of a system is closest to the true value when the derivative of cell parameters with respect to temperature is zero.

=== Results and Analysis ===
The 8x8x8 lattice grid was used for thermal optimisations for free energy values that are accurate to 0.01 meV.
{| class="wikitable"
!Temperature
(K)
!Helmholtz Free Energy
(eV)
!Lattice Constant
(a=b=c, Angstroms)
!Angle
(α=β=γ, Degrees)
!Cell Parameter
Derivatives
|-
|0
|<nowiki>-40.901906</nowiki>
|2.986565
|59.999934
|0
|-
|100
|<nowiki>-40.902418</nowiki>
|2.986564
|59.999937
|0
|-
|200
|<nowiki>-40.909374</nowiki>
|2.987613
|59.999639
|0
|-
|300
|<nowiki>-40.928119</nowiki>
|2.989403
|59.999461
|<nowiki>-0.000002</nowiki>
|-
|400
|<nowiki>-40.958587</nowiki>
|2.989403
|59.992790
|0
|-
|500
|<nowiki>-40.999427</nowiki>
|2.994156
|59.999094
|0
|-
|600
|<nowiki>-41.049305</nowiki>
|2.996846
|59.998906
|0
|-
|700
|<nowiki>-41.049305</nowiki>
|2.996846
|59.998906
|0
|-
|800
|<nowiki>-41.107108</nowiki>
|2.999674
|59.998715
|0
|-
|900
|<nowiki>-41.171879</nowiki>
|3.002623
|59.998520
|0
|-
|1000
|<nowiki>-41.243004</nowiki>
|3.005674
|59.998323
|0
|-
|'''1200'''
|'''-41.488740'''
|3.015394
|59.997763
|0
|-
|'''1400'''
|'''-41.675515'''
|3.022433
|59.997451
|0.000001
|-
|'''1600'''
|'''-41.877981'''
|3.029988
|59.996775
|0
|-
|'''1800'''
|'''-42.094435'''
|3.038336
|59.995015
|0.000001
|-
|'''2000'''
|'''-42.323668'''
|3.046923
|59.994762
|<nowiki>-0.000002</nowiki>
|-
|'''2200'''
|'''-42.564747'''
|3.056841
|59.994783
|0
|-
|'''2400'''
|'''-42.817003'''
|3.068085
|59.992390
|<nowiki>-0.00002</nowiki>
|-
|'''2600'''
|'''-43.079967'''
|3.081449
|59.993314
|0.000002
|-
|2752
|<nowiki>-43.28700906</nowiki>
|3.094344
|59.994204
|0.000002
|-
|2852
|<nowiki>-43.42654086</nowiki>
|3.105499
|59.994105
|<nowiki>-0.000001</nowiki>
|-
|2952
|<nowiki>-43.56951767</nowiki>
|3.124113
|59.994437
|<nowiki>-0.000002</nowiki>
|}
The lattice constant increases linearly only in the range 300 - 1000 K.[[File:Screenshot-Thermal Expansion.xls - LibreOffice Calc.png|none|frame]]
The thermal expansion constant of MgO in this temperature range, assuming temperature independence, is equal to <math>8.74 \times 10^{-6} \text{K}^{-1}</math> with <math>R^2 = 0.9944</math>.
[[File:Screenshot-Thermal Expansion.xls - LibreOffice Calc ALL.png|none|frame|Figure 2: Thermal expansion of MgO in all simulated temperatures.]]

If the whole temperature range 0K to 2952K was to be included, the thermal expansion of MgO is non-linear with respect to temperature. A polynomial fit to the fourth order has been determined with <math> R^2 = 1.00 </math> such that <math>\frac{a}{a_{0}} = 2.27 \times 10^{-5} T^4 - 1.04 \times 10^{-11} T^3 + 1.64 \times 10^{-8} T^2 - 8.35 \times 10^{-7} T - 2.60 \times 10^{-5}</math>, where <math>\frac{a}{a_{0}}</math> is the fractional change in the cell parameter relative to the cell parameter at T = 0 and <math> T </math> is the temperature in Kelvins. As the coefficient of thermal expansion is defined as the partial derivative of the fractional change with respect to temperature, the coefficient of thermal expansion can be expressed as <math> \alpha (T) = 9.08 \times 10^{-5} T^3 - 4.16 \times 10^{-11} T^2 + 3.28 \times 10^{-8} T - 5.85 \times 10^{-6} </math>.

== Molecular Dynamics Simulation ==

=== Computational Methodology ===
The molecular dynamics technique studies the trajectories of individual atoms by numerically solving Newton's equations of motion for an N-body system of interacting particles, where the forces of the particles are calculated by molecular mechanical force fields such as the Lennard-Jones or the Buckingham potential. It is an ''ab initio'' method of simulation applied to a large number of unit cells in order to obtain an accurate picture of the complex system in which analytical solutions are either inhibitively time-consuming or simply unavailable.

Thermodynamic properties of the isothermal-isobaric (NPT) ensemble are time-averaged over the entire space. For ergodic systems, the time averages should be identical to spatial averages provided by the quasi-harmonic approach. Due to the computational expense of this method, the system must be designed to account for available computational power by tailoring the simulation to a suitable timestep, time duration and scale. The scale of such a system must be large enough to eliminate periodic boundary condition artifacts, which in this case is suitable at the femtoseconds timestep and the 32x32x32 scale.

=== Results ===
{| class="wikitable"
!Temperature
(K)
!Averaged Final Volume
(Cubic Angstroms)
!Volume Per Unit Cell
(Cubic Angstroms)
|-
|100
|
|
|-
|200
|
|
|-
|300
|
|
|-
|400
|
|
|-
|500
|
|
|-
|600
|
|
|-
|700
|
|
|-
|800
|
|
|-
|900
|
|
|-
|1000
|615.0449322
|
|-
|1100
|
|
|-
|1200
|620.020014
|
|-
|1300
|
|
|-
|1400
|622.667288
|
|-
|1600
|626.172143
|
|-
|1800
|630.981409
|
|-
|2000
|632.410325
|
|-
|2400
|642.622061
|
|-
|2600
|648.508557
|
|-
|2800
|655.021305
|
|}

=== Analysis ===

== Conclusion ==
'''[INSERT COMPARISON MD V QHA V EXPERIMENTAL'''<ref name=":0" /><ref>C. J. Engberg, E. H. Zehms, J. Am. Ceram. Soc., 1959, 42 (6), pp. 300, doi: [http://onlinelibrary.wiley.com/store/10.1111/j.1151-2916.1959.tb12958.x/asset/j.1151-2916.1959.tb12958.x.pdf?v=1&t=ikfzf86u&s=3040cd204ec490ffeaa1cabf505c7097de78aed2&systemMessage=Wiley+Online+Library+will+be+unavailable+on+Saturday+27th+February+from+09%3A00-14%3A00+GMT+%2F+04%3A00-09%3A00+EST+%2F+17%3A00-22%3A00+SGT+for+essential+maintenance.++Apologies+for+the+inconvenience. 10.1111/j.1151-2916.1959.tb12958.x]</ref>''']'''

Both molecular dynamics and quasi-harmonic simulations display a linear thermal expansion in the temperature range 0 K – 1200 K. The quasi-harmonic method is in better agreement with experimental data in the linear range, after which the thermal expansion diverges. In the high-temperature range, thermal expansion no longer obeys linearity according to quasi-harmonic simulations as molecules are no longer bound to harmonic force fields due to increased energies. The quasi-harmonic calculation does not simulate bond dissociation as a result of large energies in the high-temperature range, and fails to account for anharmonicity experienced at high bond distances.

On the other hand, molecular dynamics indicates linearity over a much larger range of temperature. At lower temperatures where periodic condition artifacts are more prominent, molecular dynamics simulations contain some inaccuracies within its somewhat simplistic classical force models. At high temperatures, molecular dyanmics is more reliable in calculating phonon interactions, and therefore more suited to high-temperature simulations.

For both methods, the calculations break down at and near the melting point of MgO. In an NPT system, the thermal expansion coefficient and the constant-pressure heat capacity are both proportional to the rate of change in vibrational energy with respect to temperature. We expect the heat capacity to instantaneously become infinite at the transition temperature. However it is not the case as both models predict constantly rising thermal expansion coefficients beyond the transition temperature.