ChemWiki - User contributions [en]

Rep:Mod:sh5214 transition states

2017-03-10T09:11:34Z

Sh5214: /* Exercise 3 */

==Introduction==
In a potential energy surface, the minimum can either represent stable chemical compounds or it can represent the equilibrium positions between 2 or more compounds in a system. The transition state is the point with the maximum amount of energy going from one minimum to the other. The gradient for both of the points mentioned above is zero, however the curvature for a minimum is a minimum, whereas the curvature for a transition state is a saddle point.

==Exercise 1==
'''MO diagram of the TS'''

[[File:Sh5214 1 MO diagram of TS.png | 800px]]

Note: MO 16 & 17 are not degenerate, and MO 18 & 19 are not degenerate.

From the MO diagram above it can be seen that the orbitals only interact with another orbital with the same symmetry, this mean that for a reaction to occur, the orbitals interacting must have the same symmetry. This means that the orbital overlap integral for symmetric-symmetric and antisymmetric-antisymmetric interactions are non-zero, and that the the orbital overlap integral for symmetric-antisymmetric interactions are zero.

{| class="wikitable"
! colspan="2" | MO of Butadiene
|-
|
<jmol>
<jmolApplet>
<color>white</color>
<size>500</size>
<script>frame 16; mo 11; mo nodots nomesh fill translucent; mo titleformat "HOMO of butadiene (anti symmetric)"</script>
<uploadedFileContents>Sh5214 1 BUTADIENE MINIMISE PM6.LOG</uploadedFileContents>
</jmolApplet>
</jmol>
|
<jmol>
<jmolApplet>
<color>white</color>
<size>500</size>
<script>frame 16; mo 12; mo nodots nomesh fill translucent; mo titleformat "LUMO of butadiene (symmetric)"</script>
<uploadedFileContents>Sh5214 1 BUTADIENE MINIMISE PM6.LOG</uploadedFileContents>
</jmolApplet>
</jmol>
|-
! colspan="2" | MO of Ethene
|-
|
<jmol>
<jmolApplet>
<color>white</color>
<size>500</size>
<script>frame 14; mo 6; mo nodots nomesh fill translucent; mo titleformat "HOMO of ethene (symmetric)"</script>
<uploadedFileContents>Sh5214 1 ETHENE MINIMISE PM6.LOG</uploadedFileContents>
</jmolApplet>
</jmol>
|
<jmol>
<jmolApplet>
<color>white</color>
<size>500</size>
<script>frame 14; mo 7; mo nodots nomesh fill translucent; mo titleformat "LUMO of ethene (anti symmetric)"</script>
<uploadedFileContents>Sh5214 1 ETHENE MINIMISE PM6.LOG</uploadedFileContents>
</jmolApplet>
</jmol>
|-
! colspan="2" | MO of TS
|-
|
<jmol>
<jmolApplet>
<color>white</color>
<size>500</size>
<script>frame 16; mo 16; mo nodots nomesh fill translucent; mo titleformat "MO 16 of TS (anti symmetric)"</script>
<uploadedFileContents>Sh5214 1 BUTADIENE ETHENE TS PM6.LOG</uploadedFileContents>
</jmolApplet>
</jmol>
|
<jmol>
<jmolApplet>
<color>white</color>
<size>500</size>
<script>frame 16; mo 17; mo nodots nomesh fill translucent; mo titleformat "MO 17 of TS (symmetric)"</script>
<uploadedFileContents>Sh5214 1 BUTADIENE ETHENE TS PM6.LOG</uploadedFileContents>
</jmolApplet>
</jmol>
|-
|
<jmol>
<jmolApplet>
<color>white</color>
<size>500</size>
<script>frame 16; mo 18; mo nodots nomesh fill translucent; mo titleformat "MO 18 of TS (symmetric)"</script>
<uploadedFileContents>Sh5214 1 BUTADIENE ETHENE TS PM6.LOG</uploadedFileContents>
</jmolApplet>
</jmol>
|
<jmol>
<jmolApplet>
<color>white</color>
<size>500</size>
<script>frame 16; mo 19; mo nodots nomesh fill translucent; mo titleformat "MO 19 of TS (anti symmetric)"</script>
<uploadedFileContents>Sh5214 1 BUTADIENE ETHENE TS PM6.LOG</uploadedFileContents>
</jmolApplet>
</jmol>
|}

'''Bonds Lengths'''

[[File:Sh5214 1 Reaction scheme.png | 800px]]

{| class="wikitable"
|
! Reactant (Å)
! TS (Å)
! Product (Å)
|-
! Bond 1
| 1.33530
| 1.37975
| 1.50034
|-
! Bond 2
| 1.46835
| 1.41108
| 1.33766
|-
! Bond 3
| 1.3530
| 1.37976
| 1.50034
|-
! Bond 4
| 1.32731
| 1.38174
| 1.54076
|-
! Bond 5
| n/a
| 2.11459
| 1.54003
|-
! Bond 6
| n/a
| 2.11490
| 1.54004
|-
|}

As the reaction progresses, the two double bonds on butadiene and the double bond on ethene increases in length, and the single bond on butadiene decreases in length as we form the transition state. This trend continues as we go from the transition state to the product. The typical length of a sp<sup>3</sup> C-C bond is 1.54 Å, the typical length of a sp<sup>2</sup> C-C bond is 1.47 Å and the typical bond length of a alkene double bond is 1.34 Å. The Van der Waals radius of C atom is 1.7 Å, this is shorter than the partially formed bonds in the TS which had a length of around 2.11 Å.

'''Vibration that corresponds to the reaction pathway of the TS (synchronous )'''
<jmol>
<jmolApplet>
<color>white</color>
<size>500</size>
<script>frame 17; vibration 1; rotate x -20;</script>
<uploadedFileContents>Sh5214 1 BUTADIENE ETHENE TS PM6.LOG</uploadedFileContents>
</jmolApplet>
</jmol>

From the vibration above it can be seen that the formation of the two bonds is synchronous.

==Exercise 2==
'''MO diagram of the exo TS'''

[[File:Sh5214 2 MO diagram of TS.png | 800px]]

From my calculations, the results showed that the energy level ordering in the MO diagram for both the exo and the endo diels-alder reaction are the same. My MO diagram also shows that the LUMO of the 1,3-dioxole has a higher energy than the LUMO of cyclohexadiene. This means that in this reaction, the electron rich species is the dienophile and the electron deficient species is the diene, meaning that the reaction is an inverse demand diels-alder reaction. This result is expected due to the fact that the dienophile, 1,3-dioxole, has electron donating O atoms, hence it is electron rich.

{| class="wikitable"
! colspan="2" | MO of 1,3 dioxole
|-
|
<jmol>
<jmolApplet>
<color>white</color>
<size>500</size>
<script>frame 12; mo 19; mo nodots nomesh fill translucent; mo titleformat "HOMO of 1,3 dioxole"</script>
<uploadedFileContents>Sh5214 2 13DIOXOLE MINIMISE B3LYP.LOG</uploadedFileContents>
</jmolApplet>
</jmol>
|
<jmol>
<jmolApplet>
<color>white</color>
<size>500</size>
<script>frame 12; mo 20; mo nodots nomesh fill translucent; mo titleformat "LUMO of 1,3 dioxole"</script>
<uploadedFileContents>Sh5214 2 13DIOXOLE MINIMISE B3LYP.LOG</uploadedFileContents>
</jmolApplet>
</jmol>
|-
! colspan='2' | MO of cyclodexadiene
|-
|
<jmol>
<jmolApplet>
<color>white</color>
<size>500</size>
<script>frame 12; mo 22; mo nodots nomesh fill translucent; mo titleformat "HOMO of cyclohexadiene"</script>
<uploadedFileContents>Sh5214 2 CYCLOHEXADIENE MINIMISE B3LYP.LOG</uploadedFileContents>
</jmolApplet>
</jmol>
|
<jmol>
<jmolApplet>
<color>white</color>
<size>500</size>
<script>frame 12; mo 23; mo nodots nomesh fill translucent; mo titleformat "LUMO of cyclohexadiene"</script>
<uploadedFileContents>Sh5214 2 CYCLOHEXADIENE MINIMISE B3LYP.LOG</uploadedFileContents>
</jmolApplet>
</jmol>
|-
! colspan='2' | MO of exo TS
|-
|
<jmol>
<jmolApplet>
<color>white</color>
<size>500</size>
<script>frame 2; mo 40; mo nodots nomesh fill translucent; mo titleformat "MO 40 of exo TS"</script>
<uploadedFileContents>Sh5214 2 CYCLOHEXADIENE 13DIOXOLE EXO TS B3LYP MO.LOG</uploadedFileContents>
</jmolApplet>
</jmol>
|
<jmol>
<jmolApplet>
<color>white</color>
<size>500</size>
<script>frame 2; mo 41; mo nodots nomesh fill translucent; mo titleformat "MO 41 of exo TS"</script>
<uploadedFileContents>Sh5214 2 CYCLOHEXADIENE 13DIOXOLE EXO TS B3LYP MO.LOG</uploadedFileContents>
</jmolApplet>
</jmol>
|-
|
<jmol>
<jmolApplet>
<color>white</color>
<size>500</size>
<script>frame 2; mo 42; mo nodots nomesh fill translucent; mo titleformat "MO 42 of exo TS"</script>
<uploadedFileContents>Sh5214 2 CYCLOHEXADIENE 13DIOXOLE EXO TS B3LYP MO.LOG</uploadedFileContents>
</jmolApplet>
</jmol>
|
<jmol>
<jmolApplet>
<color>white</color>
<size>500</size>
<script>frame 2; mo 43; mo nodots nomesh fill translucent; mo titleformat "MO 43 of exo TS"</script>
<uploadedFileContents>Sh5214 2 CYCLOHEXADIENE 13DIOXOLE EXO TS B3LYP MO.LOG</uploadedFileContents>
</jmolApplet>
</jmol>
|-
! colspan='2' | MO of endo TS
|-
|
<jmol>
<jmolApplet>
<color>white</color>
<size>500</size>
<script>frame 2; mo 40; mo nodots nomesh fill translucent; mo titleformat "MO 40 of endo TS"</script>
<uploadedFileContents>Sh5214 2 CYCLOHEXADIENE 13DIOXOLE ENDO TS B3LYP MO.LOG</uploadedFileContents>
</jmolApplet>
</jmol>
|
<jmol>
<jmolApplet>
<color>white</color>
<size>500</size>
<script>frame 2; mo 41; mo nodots nomesh fill translucent; mo titleformat "MO 41 of endo TS"</script>
<uploadedFileContents>Sh5214 2 CYCLOHEXADIENE 13DIOXOLE ENDO TS B3LYP MO.LOG</uploadedFileContents>
</jmolApplet>
</jmol>
|-
|
<jmol>
<jmolApplet>
<color>white</color>
<size>500</size>
<script>frame 2; mo 42; mo nodots nomesh fill translucent; mo titleformat "MO 42 of endo TS"</script>
<uploadedFileContents>Sh5214 2 CYCLOHEXADIENE 13DIOXOLE ENDO TS B3LYP MO.LOG</uploadedFileContents>
</jmolApplet>
</jmol>
|
<jmol>
<jmolApplet>
<color>white</color>
<size>500</size>
<script>frame 2; mo 43; mo nodots nomesh fill translucent; mo titleformat "MO 43 of endo TS"</script>
<uploadedFileContents>Sh5214 2 CYCLOHEXADIENE 13DIOXOLE ENDO TS B3LYP MO.LOG</uploadedFileContents>
</jmolApplet>
</jmol>
|}

'''Thermochemistry'''

{| class="wikitable"
!
! Cyclohexadiene (kJ/mol)
! 1,3-dioxole (kJ/mol)
! TS (kJ/mol)
! Product (kJ/mol)
! Reaction barrier (kJ/mol)
! Reaction energy (kJ/mol)
|-
! Exo
| -612616.977634
| -701203.362889
| -1313662.70287
| -1313891.64648
| 157.637658
| -71.30595993
|-
! Endo
| -612616.977634
| -701203.362889
| -1313670.41921
| -1313895.1594
| 149.921312
| -74.8188792
|}

From the table above, it can be seen that the reaction barrier for the endo product is lower than that for the exo product, thus the endo product is the kinetically favoured product. We can also see that the reaction energy of the endo product is more negative than that of the exo product, thus the endo product is also the thermodynamically favoured product.

<jmol>
<jmolApplet>
<color>white</color>
<size>500</size>
<script>frame 2; mo 41; mo cutoff 0.01; mo nodots nomesh fill translucent; mo titleformat "HOMO (MO 41) of exo TS"</script>
<uploadedFileContents>Sh5214 2 CYCLOHEXADIENE 13DIOXOLE EXO TS B3LYP MO.LOG</uploadedFileContents>
</jmolApplet>
</jmol>

<jmol>
<jmolApplet>
<color>white</color>
<size>500</size>
<script>frame 2; mo 41; mo cutoff 0.01; mo nodots nomesh fill translucent; mo titleformat "HOMO (MO 41) of endo TS"</script>
<uploadedFileContents>Sh5214 2 CYCLOHEXADIENE 13DIOXOLE ENDO TS B3LYP MO.LOG</uploadedFileContents>
</jmolApplet>
</jmol>

From the TS MOs above, it can be seen that in the HOMO of the exo TS only one set of orbitals interact, whereas in the HOMO of the endo TS there are two sets of orbitals interacting. This means that the endo TS is more stable than the exo TS, and therefore the reaction barrier for the endo diels-alder reaction is lower than that for the exo diels-alder reaction.

===Exercise 3===

'''Exo diels alder IRC'''

[[File:Sh5214 3 Exo IRC.png]]

[[File:Sh5214 3 DA exo IRC pm6.gjf]]

'''Endo diels alder IRC'''

[[File:Sh5214 3 Endo IRC.png]]

[[File:Sh5214 3 DA endo IRC pm6.gjf]]

'''Cheletropic IRC'''

[[File:Sh5214 3 Cheletropic IRC.png]]

[[File:Sh5214 3 Cheletropic IRC pm6.gjf]]

Thermo energies

{| class="wikitable"
!
! Sulfur dioxide (kJ/mol)
! Xylylene (kJ/mol)
! TS (kJ/mol)
! Product (kJ/mol)
! Reaction barrier (kJ/mol)
! Reaction energy (kJ/mol)
|-
! Exo
| -311.421081
| 467.932399
| 241.748182
| 56.406246
| 85.236864
| -100.105072
|-
! Endo
| -311.421081
| 467.932399
| 237.762673
| 56.960227
| 81.251355
| -99.551092
|-
! Cheletropic
| -311.421081
| 467.932399
| 260.079424
| 0.013128
| 103.568106
| -156.498191
|-
! Alternate DA exo
| -311.421081
| 467.932399
| 275.819298
| 176.706665
| 119.307980
| 20.195347
|-
! Alternate DA endo
| -311.421081
| 467.932399
| 267.984805
| 172.259068
| 111.473487
| 15.747750
|}

Endo product has the lowest reaction barrier, therefore it is the kinetic product, and the cheletropic product has the most negative reaction energy, therefore it is the thermodynamic product. It can also be seen that diels-alder reaction at the alternative diene site is unfavourable due to the fact that both the endo and exo products have a higher reaction barrier and higher reaction energy than the other 3 reactions.

'''Reaction Profile'''

[[File:Sh5214 3 Reaction profile.png]]

The bonding of the 6 membered ring in xylylene goes from a diene to an 6 membered aromatic ring.

Rep:Mod:sh5214 transition states

2017-03-10T01:52:32Z

Sh5214: /* Exercise 3 */

==Introduction==
In a potential energy surface, the minimum can either represent stable chemical compounds or it can represent the equilibrium positions between 2 or more compounds in a system. The transition state is the point with the maximum amount of energy going from one minimum to the other. The gradient for both of the points mentioned above is zero, however the curvature for a minimum is a minimum, whereas the curvature for a transition state is a saddle point.

==Exercise 1==
'''MO diagram of the TS'''

[[File:Sh5214 1 MO diagram of TS.png | 800px]]

Note: MO 16 & 17 are not degenerate, and MO 18 & 19 are not degenerate.

From the MO diagram above it can be seen that the orbitals only interact with another orbital with the same symmetry, this mean that for a reaction to occur, the orbitals interacting must have the same symmetry. This means that the orbital overlap integral for symmetric-symmetric and antisymmetric-antisymmetric interactions are non-zero, and that the the orbital overlap integral for symmetric-antisymmetric interactions are zero.

{| class="wikitable"
! colspan="2" | MO of Butadiene
|-
|
<jmol>
<jmolApplet>
<color>white</color>
<size>500</size>
<script>frame 16; mo 11; mo nodots nomesh fill translucent; mo titleformat "HOMO of butadiene (anti symmetric)"</script>
<uploadedFileContents>Sh5214 1 BUTADIENE MINIMISE PM6.LOG</uploadedFileContents>
</jmolApplet>
</jmol>
|
<jmol>
<jmolApplet>
<color>white</color>
<size>500</size>
<script>frame 16; mo 12; mo nodots nomesh fill translucent; mo titleformat "LUMO of butadiene (symmetric)"</script>
<uploadedFileContents>Sh5214 1 BUTADIENE MINIMISE PM6.LOG</uploadedFileContents>
</jmolApplet>
</jmol>
|-
! colspan="2" | MO of Ethene
|-
|
<jmol>
<jmolApplet>
<color>white</color>
<size>500</size>
<script>frame 14; mo 6; mo nodots nomesh fill translucent; mo titleformat "HOMO of ethene (symmetric)"</script>
<uploadedFileContents>Sh5214 1 ETHENE MINIMISE PM6.LOG</uploadedFileContents>
</jmolApplet>
</jmol>
|
<jmol>
<jmolApplet>
<color>white</color>
<size>500</size>
<script>frame 14; mo 7; mo nodots nomesh fill translucent; mo titleformat "LUMO of ethene (anti symmetric)"</script>
<uploadedFileContents>Sh5214 1 ETHENE MINIMISE PM6.LOG</uploadedFileContents>
</jmolApplet>
</jmol>
|-
! colspan="2" | MO of TS
|-
|
<jmol>
<jmolApplet>
<color>white</color>
<size>500</size>
<script>frame 16; mo 16; mo nodots nomesh fill translucent; mo titleformat "MO 16 of TS (anti symmetric)"</script>
<uploadedFileContents>Sh5214 1 BUTADIENE ETHENE TS PM6.LOG</uploadedFileContents>
</jmolApplet>
</jmol>
|
<jmol>
<jmolApplet>
<color>white</color>
<size>500</size>
<script>frame 16; mo 17; mo nodots nomesh fill translucent; mo titleformat "MO 17 of TS (symmetric)"</script>
<uploadedFileContents>Sh5214 1 BUTADIENE ETHENE TS PM6.LOG</uploadedFileContents>
</jmolApplet>
</jmol>
|-
|
<jmol>
<jmolApplet>
<color>white</color>
<size>500</size>
<script>frame 16; mo 18; mo nodots nomesh fill translucent; mo titleformat "MO 18 of TS (symmetric)"</script>
<uploadedFileContents>Sh5214 1 BUTADIENE ETHENE TS PM6.LOG</uploadedFileContents>
</jmolApplet>
</jmol>
|
<jmol>
<jmolApplet>
<color>white</color>
<size>500</size>
<script>frame 16; mo 19; mo nodots nomesh fill translucent; mo titleformat "MO 19 of TS (anti symmetric)"</script>
<uploadedFileContents>Sh5214 1 BUTADIENE ETHENE TS PM6.LOG</uploadedFileContents>
</jmolApplet>
</jmol>
|}

'''Bonds Lengths'''

[[File:Sh5214 1 Reaction scheme.png | 800px]]

{| class="wikitable"
|
! Reactant (Å)
! TS (Å)
! Product (Å)
|-
! Bond 1
| 1.33530
| 1.37975
| 1.50034
|-
! Bond 2
| 1.46835
| 1.41108
| 1.33766
|-
! Bond 3
| 1.3530
| 1.37976
| 1.50034
|-
! Bond 4
| 1.32731
| 1.38174
| 1.54076
|-
! Bond 5
| n/a
| 2.11459
| 1.54003
|-
! Bond 6
| n/a
| 2.11490
| 1.54004
|-
|}

As the reaction progresses, the two double bonds on butadiene and the double bond on ethene increases in length, and the single bond on butadiene decreases in length as we form the transition state. This trend continues as we go from the transition state to the product. The typical length of a sp<sup>3</sup> C-C bond is 1.54 Å, the typical length of a sp<sup>2</sup> C-C bond is 1.47 Å and the typical bond length of a alkene double bond is 1.34 Å. The Van der Waals radius of C atom is 1.7 Å, this is shorter than the partially formed bonds in the TS which had a length of around 2.11 Å.

'''Vibration that corresponds to the reaction pathway of the TS (synchronous )'''
<jmol>
<jmolApplet>
<color>white</color>
<size>500</size>
<script>frame 17; vibration 1; rotate x -20;</script>
<uploadedFileContents>Sh5214 1 BUTADIENE ETHENE TS PM6.LOG</uploadedFileContents>
</jmolApplet>
</jmol>

From the vibration above it can be seen that the formation of the two bonds is synchronous.

==Exercise 2==
'''MO diagram of the exo TS'''

[[File:Sh5214 2 MO diagram of TS.png | 800px]]

From my calculations, the results showed that the energy level ordering in the MO diagram for both the exo and the endo diels-alder reaction are the same. My MO diagram also shows that the LUMO of the 1,3-dioxole has a higher energy than the LUMO of cyclohexadiene. This means that in this reaction, the electron rich species is the dienophile and the electron deficient species is the diene, meaning that the reaction is an inverse demand diels-alder reaction. This result is expected due to the fact that the dienophile, 1,3-dioxole, has electron donating O atoms, hence it is electron rich.

{| class="wikitable"
! colspan="2" | MO of 1,3 dioxole
|-
|
<jmol>
<jmolApplet>
<color>white</color>
<size>500</size>
<script>frame 12; mo 19; mo nodots nomesh fill translucent; mo titleformat "HOMO of 1,3 dioxole"</script>
<uploadedFileContents>Sh5214 2 13DIOXOLE MINIMISE B3LYP.LOG</uploadedFileContents>
</jmolApplet>
</jmol>
|
<jmol>
<jmolApplet>
<color>white</color>
<size>500</size>
<script>frame 12; mo 20; mo nodots nomesh fill translucent; mo titleformat "LUMO of 1,3 dioxole"</script>
<uploadedFileContents>Sh5214 2 13DIOXOLE MINIMISE B3LYP.LOG</uploadedFileContents>
</jmolApplet>
</jmol>
|-
! colspan='2' | MO of cyclodexadiene
|-
|
<jmol>
<jmolApplet>
<color>white</color>
<size>500</size>
<script>frame 12; mo 22; mo nodots nomesh fill translucent; mo titleformat "HOMO of cyclohexadiene"</script>
<uploadedFileContents>Sh5214 2 CYCLOHEXADIENE MINIMISE B3LYP.LOG</uploadedFileContents>
</jmolApplet>
</jmol>
|
<jmol>
<jmolApplet>
<color>white</color>
<size>500</size>
<script>frame 12; mo 23; mo nodots nomesh fill translucent; mo titleformat "LUMO of cyclohexadiene"</script>
<uploadedFileContents>Sh5214 2 CYCLOHEXADIENE MINIMISE B3LYP.LOG</uploadedFileContents>
</jmolApplet>
</jmol>
|-
! colspan='2' | MO of exo TS
|-
|
<jmol>
<jmolApplet>
<color>white</color>
<size>500</size>
<script>frame 2; mo 40; mo nodots nomesh fill translucent; mo titleformat "MO 40 of exo TS"</script>
<uploadedFileContents>Sh5214 2 CYCLOHEXADIENE 13DIOXOLE EXO TS B3LYP MO.LOG</uploadedFileContents>
</jmolApplet>
</jmol>
|
<jmol>
<jmolApplet>
<color>white</color>
<size>500</size>
<script>frame 2; mo 41; mo nodots nomesh fill translucent; mo titleformat "MO 41 of exo TS"</script>
<uploadedFileContents>Sh5214 2 CYCLOHEXADIENE 13DIOXOLE EXO TS B3LYP MO.LOG</uploadedFileContents>
</jmolApplet>
</jmol>
|-
|
<jmol>
<jmolApplet>
<color>white</color>
<size>500</size>
<script>frame 2; mo 42; mo nodots nomesh fill translucent; mo titleformat "MO 42 of exo TS"</script>
<uploadedFileContents>Sh5214 2 CYCLOHEXADIENE 13DIOXOLE EXO TS B3LYP MO.LOG</uploadedFileContents>
</jmolApplet>
</jmol>
|
<jmol>
<jmolApplet>
<color>white</color>
<size>500</size>
<script>frame 2; mo 43; mo nodots nomesh fill translucent; mo titleformat "MO 43 of exo TS"</script>
<uploadedFileContents>Sh5214 2 CYCLOHEXADIENE 13DIOXOLE EXO TS B3LYP MO.LOG</uploadedFileContents>
</jmolApplet>
</jmol>
|-
! colspan='2' | MO of endo TS
|-
|
<jmol>
<jmolApplet>
<color>white</color>
<size>500</size>
<script>frame 2; mo 40; mo nodots nomesh fill translucent; mo titleformat "MO 40 of endo TS"</script>
<uploadedFileContents>Sh5214 2 CYCLOHEXADIENE 13DIOXOLE ENDO TS B3LYP MO.LOG</uploadedFileContents>
</jmolApplet>
</jmol>
|
<jmol>
<jmolApplet>
<color>white</color>
<size>500</size>
<script>frame 2; mo 41; mo nodots nomesh fill translucent; mo titleformat "MO 41 of endo TS"</script>
<uploadedFileContents>Sh5214 2 CYCLOHEXADIENE 13DIOXOLE ENDO TS B3LYP MO.LOG</uploadedFileContents>
</jmolApplet>
</jmol>
|-
|
<jmol>
<jmolApplet>
<color>white</color>
<size>500</size>
<script>frame 2; mo 42; mo nodots nomesh fill translucent; mo titleformat "MO 42 of endo TS"</script>
<uploadedFileContents>Sh5214 2 CYCLOHEXADIENE 13DIOXOLE ENDO TS B3LYP MO.LOG</uploadedFileContents>
</jmolApplet>
</jmol>
|
<jmol>
<jmolApplet>
<color>white</color>
<size>500</size>
<script>frame 2; mo 43; mo nodots nomesh fill translucent; mo titleformat "MO 43 of endo TS"</script>
<uploadedFileContents>Sh5214 2 CYCLOHEXADIENE 13DIOXOLE ENDO TS B3LYP MO.LOG</uploadedFileContents>
</jmolApplet>
</jmol>
|}

'''Thermochemistry'''

{| class="wikitable"
!
! Cyclohexadiene (kJ/mol)
! 1,3-dioxole (kJ/mol)
! TS (kJ/mol)
! Product (kJ/mol)
! Reaction barrier (kJ/mol)
! Reaction energy (kJ/mol)
|-
! Exo
| -612616.977634
| -701203.362889
| -1313662.70287
| -1313891.64648
| 157.637658
| -71.30595993
|-
! Endo
| -612616.977634
| -701203.362889
| -1313670.41921
| -1313895.1594
| 149.921312
| -74.8188792
|}

From the table above, it can be seen that the reaction barrier for the endo product is lower than that for the exo product, thus the endo product is the kinetically favoured product. We can also see that the reaction energy of the endo product is more negative than that of the exo product, thus the endo product is also the thermodynamically favoured product.

<jmol>
<jmolApplet>
<color>white</color>
<size>500</size>
<script>frame 2; mo 41; mo cutoff 0.01; mo nodots nomesh fill translucent; mo titleformat "HOMO (MO 41) of exo TS"</script>
<uploadedFileContents>Sh5214 2 CYCLOHEXADIENE 13DIOXOLE EXO TS B3LYP MO.LOG</uploadedFileContents>
</jmolApplet>
</jmol>

<jmol>
<jmolApplet>
<color>white</color>
<size>500</size>
<script>frame 2; mo 41; mo cutoff 0.01; mo nodots nomesh fill translucent; mo titleformat "HOMO (MO 41) of endo TS"</script>
<uploadedFileContents>Sh5214 2 CYCLOHEXADIENE 13DIOXOLE ENDO TS B3LYP MO.LOG</uploadedFileContents>
</jmolApplet>
</jmol>

From the TS MOs above, it can be seen that in the HOMO of the exo TS only one set of orbitals interact, whereas in the HOMO of the endo TS there are two sets of orbitals interacting. This means that the endo TS is more stable than the exo TS, and therefore the reaction barrier for the endo diels-alder reaction is lower than that for the exo diels-alder reaction.

===Exercise 3===

'''Exo diels alder IRC'''

[[File:Sh5214 3 Exo IRC.png]]

[[File:Sh5214 3 DA exo IRC pm6.gjf]]

'''Endo diels alder IRC'''

[[File:Sh5214 3 Endo IRC.png]]

[[File:Sh5214 3 DA endo IRC pm6.gjf]]

'''Cheletropic IRC'''

[[File:Sh5214 3 Cheletropic IRC.png]]

[[File:Sh5214 3 Cheletropic IRC pm6.gjf]]

Thermo energies

{| class="wikitable"
!
! Sulfur dioxide (kJ/mol)
! Xylylene (kJ/mol)
! TS (kJ/mol)
! Product (kJ/mol)
! Reaction barrier (kJ/mol)
! Reaction energy (kJ/mol)
|-
! Exo
| -311.4210807
| 467.932399
| 241.748182
| 56.4062463
| 85.236864
| -100.1050716
|-
! Endo
| -311.4210807
| 467.932399
| 237.762673
| 56.9602268
| 2.98593081
| -28.2766372
|-
! Cheletropic
| -311.4210807
| 467.932399
| 260.079424
| 0.013127501
| 103.568106
| -156.5244454
|}

Endo product has the lowest reaction barrier, therefore it is the kinetic product, and the cheletropic product has the most negative reaction energy, therefore it is the thermodynamic product.

'''Reaction Profile'''

[[File:Sh5214 3 Reaction profile.png]]

The bonding of the 6 membered ring in xylylene goes from a diene to an 6 membered aromatic ring.

Rep:Mod:sh5214 transition states

2017-03-10T01:43:32Z

Sh5214: /* Exercise 3 */

==Introduction==
In a potential energy surface, the minimum can either represent stable chemical compounds or it can represent the equilibrium positions between 2 or more compounds in a system. The transition state is the point with the maximum amount of energy going from one minimum to the other. The gradient for both of the points mentioned above is zero, however the curvature for a minimum is a minimum, whereas the curvature for a transition state is a saddle point.

==Exercise 1==
'''MO diagram of the TS'''

[[File:Sh5214 1 MO diagram of TS.png | 800px]]

Note: MO 16 & 17 are not degenerate, and MO 18 & 19 are not degenerate.

From the MO diagram above it can be seen that the orbitals only interact with another orbital with the same symmetry, this mean that for a reaction to occur, the orbitals interacting must have the same symmetry. This means that the orbital overlap integral for symmetric-symmetric and antisymmetric-antisymmetric interactions are non-zero, and that the the orbital overlap integral for symmetric-antisymmetric interactions are zero.

{| class="wikitable"
! colspan="2" | MO of Butadiene
|-
|
<jmol>
<jmolApplet>
<color>white</color>
<size>500</size>
<script>frame 16; mo 11; mo nodots nomesh fill translucent; mo titleformat "HOMO of butadiene (anti symmetric)"</script>
<uploadedFileContents>Sh5214 1 BUTADIENE MINIMISE PM6.LOG</uploadedFileContents>
</jmolApplet>
</jmol>
|
<jmol>
<jmolApplet>
<color>white</color>
<size>500</size>
<script>frame 16; mo 12; mo nodots nomesh fill translucent; mo titleformat "LUMO of butadiene (symmetric)"</script>
<uploadedFileContents>Sh5214 1 BUTADIENE MINIMISE PM6.LOG</uploadedFileContents>
</jmolApplet>
</jmol>
|-
! colspan="2" | MO of Ethene
|-
|
<jmol>
<jmolApplet>
<color>white</color>
<size>500</size>
<script>frame 14; mo 6; mo nodots nomesh fill translucent; mo titleformat "HOMO of ethene (symmetric)"</script>
<uploadedFileContents>Sh5214 1 ETHENE MINIMISE PM6.LOG</uploadedFileContents>
</jmolApplet>
</jmol>
|
<jmol>
<jmolApplet>
<color>white</color>
<size>500</size>
<script>frame 14; mo 7; mo nodots nomesh fill translucent; mo titleformat "LUMO of ethene (anti symmetric)"</script>
<uploadedFileContents>Sh5214 1 ETHENE MINIMISE PM6.LOG</uploadedFileContents>
</jmolApplet>
</jmol>
|-
! colspan="2" | MO of TS
|-
|
<jmol>
<jmolApplet>
<color>white</color>
<size>500</size>
<script>frame 16; mo 16; mo nodots nomesh fill translucent; mo titleformat "MO 16 of TS (anti symmetric)"</script>
<uploadedFileContents>Sh5214 1 BUTADIENE ETHENE TS PM6.LOG</uploadedFileContents>
</jmolApplet>
</jmol>
|
<jmol>
<jmolApplet>
<color>white</color>
<size>500</size>
<script>frame 16; mo 17; mo nodots nomesh fill translucent; mo titleformat "MO 17 of TS (symmetric)"</script>
<uploadedFileContents>Sh5214 1 BUTADIENE ETHENE TS PM6.LOG</uploadedFileContents>
</jmolApplet>
</jmol>
|-
|
<jmol>
<jmolApplet>
<color>white</color>
<size>500</size>
<script>frame 16; mo 18; mo nodots nomesh fill translucent; mo titleformat "MO 18 of TS (symmetric)"</script>
<uploadedFileContents>Sh5214 1 BUTADIENE ETHENE TS PM6.LOG</uploadedFileContents>
</jmolApplet>
</jmol>
|
<jmol>
<jmolApplet>
<color>white</color>
<size>500</size>
<script>frame 16; mo 19; mo nodots nomesh fill translucent; mo titleformat "MO 19 of TS (anti symmetric)"</script>
<uploadedFileContents>Sh5214 1 BUTADIENE ETHENE TS PM6.LOG</uploadedFileContents>
</jmolApplet>
</jmol>
|}

'''Bonds Lengths'''

[[File:Sh5214 1 Reaction scheme.png | 800px]]

{| class="wikitable"
|
! Reactant (Å)
! TS (Å)
! Product (Å)
|-
! Bond 1
| 1.33530
| 1.37975
| 1.50034
|-
! Bond 2
| 1.46835
| 1.41108
| 1.33766
|-
! Bond 3
| 1.3530
| 1.37976
| 1.50034
|-
! Bond 4
| 1.32731
| 1.38174
| 1.54076
|-
! Bond 5
| n/a
| 2.11459
| 1.54003
|-
! Bond 6
| n/a
| 2.11490
| 1.54004
|-
|}

As the reaction progresses, the two double bonds on butadiene and the double bond on ethene increases in length, and the single bond on butadiene decreases in length as we form the transition state. This trend continues as we go from the transition state to the product. The typical length of a sp<sup>3</sup> C-C bond is 1.54 Å, the typical length of a sp<sup>2</sup> C-C bond is 1.47 Å and the typical bond length of a alkene double bond is 1.34 Å. The Van der Waals radius of C atom is 1.7 Å, this is shorter than the partially formed bonds in the TS which had a length of around 2.11 Å.

'''Vibration that corresponds to the reaction pathway of the TS (synchronous )'''
<jmol>
<jmolApplet>
<color>white</color>
<size>500</size>
<script>frame 17; vibration 1; rotate x -20;</script>
<uploadedFileContents>Sh5214 1 BUTADIENE ETHENE TS PM6.LOG</uploadedFileContents>
</jmolApplet>
</jmol>

From the vibration above it can be seen that the formation of the two bonds is synchronous.

==Exercise 2==
'''MO diagram of the exo TS'''

[[File:Sh5214 2 MO diagram of TS.png | 800px]]

From my calculations, the results showed that the energy level ordering in the MO diagram for both the exo and the endo diels-alder reaction are the same. My MO diagram also shows that the LUMO of the 1,3-dioxole has a higher energy than the LUMO of cyclohexadiene. This means that in this reaction, the electron rich species is the dienophile and the electron deficient species is the diene, meaning that the reaction is an inverse demand diels-alder reaction. This result is expected due to the fact that the dienophile, 1,3-dioxole, has electron donating O atoms, hence it is electron rich.

{| class="wikitable"
! colspan="2" | MO of 1,3 dioxole
|-
|
<jmol>
<jmolApplet>
<color>white</color>
<size>500</size>
<script>frame 12; mo 19; mo nodots nomesh fill translucent; mo titleformat "HOMO of 1,3 dioxole"</script>
<uploadedFileContents>Sh5214 2 13DIOXOLE MINIMISE B3LYP.LOG</uploadedFileContents>
</jmolApplet>
</jmol>
|
<jmol>
<jmolApplet>
<color>white</color>
<size>500</size>
<script>frame 12; mo 20; mo nodots nomesh fill translucent; mo titleformat "LUMO of 1,3 dioxole"</script>
<uploadedFileContents>Sh5214 2 13DIOXOLE MINIMISE B3LYP.LOG</uploadedFileContents>
</jmolApplet>
</jmol>
|-
! colspan='2' | MO of cyclodexadiene
|-
|
<jmol>
<jmolApplet>
<color>white</color>
<size>500</size>
<script>frame 12; mo 22; mo nodots nomesh fill translucent; mo titleformat "HOMO of cyclohexadiene"</script>
<uploadedFileContents>Sh5214 2 CYCLOHEXADIENE MINIMISE B3LYP.LOG</uploadedFileContents>
</jmolApplet>
</jmol>
|
<jmol>
<jmolApplet>
<color>white</color>
<size>500</size>
<script>frame 12; mo 23; mo nodots nomesh fill translucent; mo titleformat "LUMO of cyclohexadiene"</script>
<uploadedFileContents>Sh5214 2 CYCLOHEXADIENE MINIMISE B3LYP.LOG</uploadedFileContents>
</jmolApplet>
</jmol>
|-
! colspan='2' | MO of exo TS
|-
|
<jmol>
<jmolApplet>
<color>white</color>
<size>500</size>
<script>frame 2; mo 40; mo nodots nomesh fill translucent; mo titleformat "MO 40 of exo TS"</script>
<uploadedFileContents>Sh5214 2 CYCLOHEXADIENE 13DIOXOLE EXO TS B3LYP MO.LOG</uploadedFileContents>
</jmolApplet>
</jmol>
|
<jmol>
<jmolApplet>
<color>white</color>
<size>500</size>
<script>frame 2; mo 41; mo nodots nomesh fill translucent; mo titleformat "MO 41 of exo TS"</script>
<uploadedFileContents>Sh5214 2 CYCLOHEXADIENE 13DIOXOLE EXO TS B3LYP MO.LOG</uploadedFileContents>
</jmolApplet>
</jmol>
|-
|
<jmol>
<jmolApplet>
<color>white</color>
<size>500</size>
<script>frame 2; mo 42; mo nodots nomesh fill translucent; mo titleformat "MO 42 of exo TS"</script>
<uploadedFileContents>Sh5214 2 CYCLOHEXADIENE 13DIOXOLE EXO TS B3LYP MO.LOG</uploadedFileContents>
</jmolApplet>
</jmol>
|
<jmol>
<jmolApplet>
<color>white</color>
<size>500</size>
<script>frame 2; mo 43; mo nodots nomesh fill translucent; mo titleformat "MO 43 of exo TS"</script>
<uploadedFileContents>Sh5214 2 CYCLOHEXADIENE 13DIOXOLE EXO TS B3LYP MO.LOG</uploadedFileContents>
</jmolApplet>
</jmol>
|-
! colspan='2' | MO of endo TS
|-
|
<jmol>
<jmolApplet>
<color>white</color>
<size>500</size>
<script>frame 2; mo 40; mo nodots nomesh fill translucent; mo titleformat "MO 40 of endo TS"</script>
<uploadedFileContents>Sh5214 2 CYCLOHEXADIENE 13DIOXOLE ENDO TS B3LYP MO.LOG</uploadedFileContents>
</jmolApplet>
</jmol>
|
<jmol>
<jmolApplet>
<color>white</color>
<size>500</size>
<script>frame 2; mo 41; mo nodots nomesh fill translucent; mo titleformat "MO 41 of endo TS"</script>
<uploadedFileContents>Sh5214 2 CYCLOHEXADIENE 13DIOXOLE ENDO TS B3LYP MO.LOG</uploadedFileContents>
</jmolApplet>
</jmol>
|-
|
<jmol>
<jmolApplet>
<color>white</color>
<size>500</size>
<script>frame 2; mo 42; mo nodots nomesh fill translucent; mo titleformat "MO 42 of endo TS"</script>
<uploadedFileContents>Sh5214 2 CYCLOHEXADIENE 13DIOXOLE ENDO TS B3LYP MO.LOG</uploadedFileContents>
</jmolApplet>
</jmol>
|
<jmol>
<jmolApplet>
<color>white</color>
<size>500</size>
<script>frame 2; mo 43; mo nodots nomesh fill translucent; mo titleformat "MO 43 of endo TS"</script>
<uploadedFileContents>Sh5214 2 CYCLOHEXADIENE 13DIOXOLE ENDO TS B3LYP MO.LOG</uploadedFileContents>
</jmolApplet>
</jmol>
|}

'''Thermochemistry'''

{| class="wikitable"
!
! Cyclohexadiene (kJ/mol)
! 1,3-dioxole (kJ/mol)
! TS (kJ/mol)
! Product (kJ/mol)
! Reaction barrier (kJ/mol)
! Reaction energy (kJ/mol)
|-
! Exo
| -612616.977634
| -701203.362889
| -1313662.70287
| -1313891.64648
| 157.637658
| -71.30595993
|-
! Endo
| -612616.977634
| -701203.362889
| -1313670.41921
| -1313895.1594
| 149.921312
| -74.8188792
|}

From the table above, it can be seen that the reaction barrier for the endo product is lower than that for the exo product, thus the endo product is the kinetically favoured product. We can also see that the reaction energy of the endo product is more negative than that of the exo product, thus the endo product is also the thermodynamically favoured product.

<jmol>
<jmolApplet>
<color>white</color>
<size>500</size>
<script>frame 2; mo 41; mo cutoff 0.01; mo nodots nomesh fill translucent; mo titleformat "HOMO (MO 41) of exo TS"</script>
<uploadedFileContents>Sh5214 2 CYCLOHEXADIENE 13DIOXOLE EXO TS B3LYP MO.LOG</uploadedFileContents>
</jmolApplet>
</jmol>

<jmol>
<jmolApplet>
<color>white</color>
<size>500</size>
<script>frame 2; mo 41; mo cutoff 0.01; mo nodots nomesh fill translucent; mo titleformat "HOMO (MO 41) of endo TS"</script>
<uploadedFileContents>Sh5214 2 CYCLOHEXADIENE 13DIOXOLE ENDO TS B3LYP MO.LOG</uploadedFileContents>
</jmolApplet>
</jmol>

From the TS MOs above, it can be seen that in the HOMO of the exo TS only one set of orbitals interact, whereas in the HOMO of the endo TS there are two sets of orbitals interacting. This means that the endo TS is more stable than the exo TS, and therefore the reaction barrier for the endo diels-alder reaction is lower than that for the exo diels-alder reaction.

===Exercise 3===

'''Exo diels alder IRC'''

[[File:Sh5214 3 Exo IRC.png]]

[[File:Sh5214 3 DA exo IRC pm6.gjf]]

'''Endo diels alder IRC'''

[[File:Sh5214 3 Endo IRC.png]]

[[File:Sh5214 3 DA endo IRC pm6.gjf]]

'''Cheletropic IRC'''

[[File:Sh5214 3 Cheletropic IRC.png]]

[[File:Sh5214 3 Cheletropic IRC pm6.gjf]]

Thermo energies

{| class="wikitable"
!
! Sulfur dioxide (kJ/mol)
! Xylylene (kJ/mol)
! TS (kJ/mol)
! Product (kJ/mol)
! Reaction barrier (kJ/mol)
! Reaction energy (kJ/mol)
|-
! Exo
| -311.4210807
| 467.932399
| 241.748182
| 56.4062463
| 85.236864
| -100.1050716
|-
! Endo
| -311.4210807
| 467.932399
| 237.762673
| 56.9602268
| 2.98593081
| -28.2766372
|-
! Cheletropic
| -311.4210807
| 467.932399
| 260.079424
| 0.013127501
| 103.568106
| -156.5244454
|}

Endo product has the lowest reaction barrier, therefore it is the kinetic product, and the cheletropic product has the most negative reaction energy, therefore it is the thermodynamic product.

'''Reaction Profile'''

[[File:Sh5214 3 Reaction profile.png]]

File:Sh5214 3 Reaction profile.png

2017-03-10T01:42:23Z

Sh5214:

Rep:Mod:sh5214 transition states

2017-03-10T01:22:48Z

2017-03-09T14:51:55Z

Sh5214:

File:Sh5214 2 MO diagram of TS.png

2017-03-09T14:40:28Z

Sh5214:

Rep:Mod:sh5214 transition states

==Introduction
In a potential energy surface, the minimas can either represent stable chemical compounds or it can represent the equilibrium positions between 2 or more compounds in a system.

==Exercise 1==
<jmol>
<jmolApplet>
<color>white</color>
<size>500</size>
<script>frame 16; mo 11; mo nodots nomesh fill translucent; mo titleformat "HOMO of butadiene"</script>
<uploadedFileContents>Sh5214 1 BUTADIENE MINIMISE PM6.LOG</uploadedFileContents>
</jmolApplet>
</jmol>

<jmol>
<jmolApplet>
<color>white</color>
<size>500</size>
<script>frame 16; mo 12; mo nodots nomesh fill translucent; mo titleformat "LUMO of butadiene"</script>
<uploadedFileContents>Sh5214 1 BUTADIENE MINIMISE PM6.LOG</uploadedFileContents>
</jmolApplet>
</jmol>

<jmol>
<jmolApplet>
<color>white</color>
<size>500</size>
<script>frame 14; mo 6; mo nodots nomesh fill translucent; mo titleformat "HOMO of ethene"</script>
<uploadedFileContents>Sh5214 1 ETHENE MINIMISE PM6.LOG</uploadedFileContents>
</jmolApplet>
</jmol>

<jmol>
<jmolApplet>
<color>white</color>
<size>500</size>
<script>frame 14; mo 7; mo nodots nomesh fill translucent; mo titleformat "LUMO of ethene"</script>
<uploadedFileContents>Sh5214 1 ETHENE MINIMISE PM6.LOG</uploadedFileContents>
</jmolApplet>
</jmol>

<jmol>
<jmolApplet>
<color>white</color>
<size>500</size>
<script>frame 16; mo 16; mo nodots nomesh fill translucent; mo titleformat "MO 16 of TS"</script>
<uploadedFileContents>Sh5214 1 BUTADIENE ETHENE TS PM6.LOG</uploadedFileContents>
</jmolApplet>
</jmol>

<jmol>
<jmolApplet>
<color>white</color>
<size>500</size>
<script>frame 16; mo 17; mo nodots nomesh fill translucent; mo titleformat "MO 17 of TS"</script>
<uploadedFileContents>Sh5214 1 BUTADIENE ETHENE TS PM6.LOG</uploadedFileContents>
</jmolApplet>
</jmol>

<jmol>
<jmolApplet>
<color>white</color>
<size>500</size>
<script>frame 16; mo 18; mo nodots nomesh fill translucent; mo titleformat "MO 18 of TS"</script>
<uploadedFileContents>Sh5214 1 BUTADIENE ETHENE TS PM6.LOG</uploadedFileContents>
</jmolApplet>
</jmol>

<jmol>
<jmolApplet>
<color>white</color>
<size>500</size>
<script>frame 16; mo 19; mo nodots nomesh fill translucent; mo titleformat "MO 19 of TS"</script>
<uploadedFileContents>Sh5214 1 BUTADIENE ETHENE TS PM6.LOG</uploadedFileContents>
</jmolApplet>
</jmol>

<jmol>
<jmolApplet>
<color>white</color>
<size>500</size>
<script>frame 17; vibration 1; rotate x -20;</script>
<uploadedFileContents>Sh5214 1 BUTADIENE ETHENE TS PM6.LOG</uploadedFileContents>
</jmolApplet>
</jmol>

Rep:Mod:sh5214 transition states

2017-03-02T13:33:34Z

Sh5214:

<jmol>
<jmolApplet>
<color>white</color>
<size>500</size>
<script>frame 16; mo 11; mo nodots nomesh fill translucent; mo titleformat "HOMO of butadiene"</script>
<uploadedFileContents>Sh5214 1 BUTADIENE MINIMISE PM6.LOG</uploadedFileContents>
</jmolApplet>
</jmol>

<jmol>
<jmolApplet>
<color>white</color>
<size>500</size>
<script>frame 16; mo 12; mo nodots nomesh fill translucent; mo titleformat "LUMO of butadiene"</script>
<uploadedFileContents>Sh5214 1 BUTADIENE MINIMISE PM6.LOG</uploadedFileContents>
</jmolApplet>
</jmol>

<jmol>
<jmolApplet>
<color>white</color>
<size>500</size>
<script>frame 14; mo 6; mo nodots nomesh fill translucent; mo titleformat "HOMO of ethene"</script>
<uploadedFileContents>Sh5214 1 ETHENE MINIMISE PM6.LOG</uploadedFileContents>
</jmolApplet>
</jmol>

<jmol>
<jmolApplet>
<color>white</color>
<size>500</size>
<script>frame 14; mo 7; mo nodots nomesh fill translucent; mo titleformat "LUMO of ethene"</script>
<uploadedFileContents>Sh5214 1 ETHENE MINIMISE PM6.LOG</uploadedFileContents>
</jmolApplet>
</jmol>

<jmol>
<jmolApplet>
<color>white</color>
<size>500</size>
<script>frame 16; mo 16; mo nodots nomesh fill translucent; mo titleformat "MO 16 of TS"</script>
<uploadedFileContents>Sh5214 1 BUTADIENE ETHENE TS PM6.LOG</uploadedFileContents>
</jmolApplet>
</jmol>

<jmol>
<jmolApplet>
<color>white</color>
<size>500</size>
<script>frame 16; mo 17; mo nodots nomesh fill translucent; mo titleformat "MO 17 of TS"</script>
<uploadedFileContents>Sh5214 1 BUTADIENE ETHENE TS PM6.LOG</uploadedFileContents>
</jmolApplet>
</jmol>

<jmol>
<jmolApplet>
<color>white</color>
<size>500</size>
<script>frame 16; mo 18; mo nodots nomesh fill translucent; mo titleformat "MO 18 of TS"</script>
<uploadedFileContents>Sh5214 1 BUTADIENE ETHENE TS PM6.LOG</uploadedFileContents>
</jmolApplet>
</jmol>

<jmol>
<jmolApplet>
<color>white</color>
<size>500</size>
<script>frame 16; mo 19; mo nodots nomesh fill translucent; mo titleformat "MO 19 of TS"</script>
<uploadedFileContents>Sh5214 1 BUTADIENE ETHENE TS PM6.LOG</uploadedFileContents>
</jmolApplet>
</jmol>

<jmol>
<jmolApplet>
<color>white</color>
<size>500</size>
<script>frame 17; vibration 1; rotate x -20;</script>
<uploadedFileContents>Sh5214 1 BUTADIENE ETHENE TS PM6.LOG</uploadedFileContents>
</jmolApplet>
</jmol>

Rep:Mod:sh5214 transition states

2017-03-02T11:34:30Z

Sh5214: Created page with "<jmol> <jmolApplet> <color>white</color> <size>500</size> <script>frame 2; mo 11; mo nodots nomesh fill translucent; mo titleformat "HOMO of butadiene"</script> <uploaded..."

<jmol>
<jmolApplet>
<color>white</color>
<size>500</size>
<script>frame 2; mo 11; mo nodots nomesh fill translucent; mo titleformat "HOMO of butadiene"</script>
<uploadedFileContents>Sh5214 1 BUTADIENE MINIMISE PM6.LOG</uploadedFileContents>
</jmolApplet>
</jmol>

<jmol>
<jmolApplet>
<color>white</color>
<size>500</size>
<script>frame 2; mo 12; mo nodots nomesh fill translucent; mo titleformat "LUMO of butadiene"</script>
<uploadedFileContents>Sh5214 1 BUTADIENE MINIMISE PM6.LOG</uploadedFileContents>
</jmolApplet>
</jmol>

<jmol>
<jmolApplet>
<color>white</color>
<size>500</size>
<script>frame 2; mo 6; mo nodots nomesh fill translucent; mo titleformat "HOMO of ethene"</script>
<uploadedFileContents>Sh5214 1 ETHENE MINIMISE PM6.LOG</uploadedFileContents>
</jmolApplet>
</jmol>

<jmol>
<jmolApplet>
<color>white</color>
<size>500</size>
<script>frame 2; mo 7; mo nodots nomesh fill translucent; mo titleformat "LUMO of ethene"</script>
<uploadedFileContents>Sh5214 1 ETHENE MINIMISE PM6.LOG</uploadedFileContents>
</jmolApplet>
</jmol>

<jmol>
<jmolApplet>
<color>white</color>
<size>500</size>
<script>frame 2; mo 16; mo nodots nomesh fill translucent; mo titleformat "MO 16 of TS"</script>
<uploadedFileContents>Sh5214 1 BUTADIENE ETHENE TS PM6.LOG</uploadedFileContents>
</jmolApplet>
</jmol>

<jmol>
<jmolApplet>
<color>white</color>
<size>500</size>
<script>frame 2; mo 17; mo nodots nomesh fill translucent; mo titleformat "MO 17 of TS"</script>
<uploadedFileContents>Sh5214 1 BUTADIENE ETHENE TS PM6.LOG</uploadedFileContents>
</jmolApplet>
</jmol>

<jmol>
<jmolApplet>
<color>white</color>
<size>500</size>
<script>frame 2; mo 18; mo nodots nomesh fill translucent; mo titleformat "MO 18 of TS"</script>
<uploadedFileContents>Sh5214 1 BUTADIENE ETHENE TS PM6.LOG</uploadedFileContents>
</jmolApplet>
</jmol>

<jmol>
<jmolApplet>
<color>white</color>
<size>500</size>
<script>frame 2; mo 19; mo nodots nomesh fill translucent; mo titleformat "MO 19 of TS"</script>
<uploadedFileContents>Sh5214 1 BUTADIENE ETHENE TS PM6.LOG</uploadedFileContents>
</jmolApplet>
</jmol>

File:Sh5214 1 BUTADIENE ETHENE TS PM6.LOG

2017-03-02T11:21:39Z

Sh5214:

File:Sh5214 1 ETHENE MINIMISE PM6.LOG

2017-03-02T11:12:25Z

Sh5214:

File:Sh5214 1 BUTADIENE MINIMISE PM6.LOG

2017-03-02T10:53:52Z

Sh5214:

Rep:Mod:sh5214 liquid simulation

2017-02-24T11:03:50Z

Sh5214: /* */

== ==
'''<big>TASK</big>: Open the file HO.xls. In it, the velocity-Verlet algorithm is used to model the behaviour of a classical harmonic oscillator. Complete the three columns "ANALYTICAL", "ERROR", and "ENERGY": "ANALYTICAL" should contain the value of the classical solution for the position at time <math>t</math>, "ERROR" should contain the ''absolute'' difference between "ANALYTICAL" and the velocity-Verlet solution (i.e. ERROR should always be positive -- make sure you leave the half step rows blank!), and "ENERGY" should contain the total energy of the oscillator for the velocity-Verlet solution. Remember that the position of a classical harmonic oscillator is given by <math> x\left(t\right) = A\cos\left(\omega t + \phi\right)</math> (the values of <math>A</math>, <math>\omega</math>, and <math>\phi</math> are worked out for you in the sheet).'''

[[File:Sh5214 liquid sim theory analytical.png]] [[File:Sh5214 liquid sim theory energy.png]] [[File:Sh5214 liquid sim theory error.png ‎]]

== ==
'''<big>TASK</big>: For the default timestep value, 0.1, estimate the positions of the maxima in the ERROR column as a function of time. Make a plot showing these values as a function of time, and fit an appropriate function to the data.'''

[[File:Sh5214 liquid sim theory error maxima.png ]]

== ==
'''<big>TASK</big>: Experiment with different values of the timestep. What sort of a timestep do you need to use to ensure that the total energy does not change by more than 1% over the course of your "simulation"? Why do you think it is important to monitor the total energy of a physical system when modelling its behaviour numerically?'''

The maximum timestep that doesn't result in more than an 1% change in energy is around 0.45. The reason why the energy is monitored is to make sure the energy change between each timestep is not too large. If there is a large energy difference, it would mean that the overall change in velocity of the particles in the system is too large, which would lead to not obtaining the full amount of information from the simulation.

== ==
'''<big>TASK:</big> For a single Lennard-Jones interaction, <math>\phi\left(r\right) = 4\epsilon \left( \frac{\sigma^{12}}{r^{12}} - \frac{\sigma^6}{r^6} \right)</math>, find the separation, <math>r_0</math>, at which the potential energy is zero. What is the force at this separation? Find the equilibrium separation, <math>r_{eq}</math>, and work out the well depth (<math>\phi\left(r_{eq}\right)</math>). Evaluate the integrals <math>\int_{2\sigma}^\infty \phi\left(r\right)\mathrm{d}r</math>, <math>\int_{2.5\sigma}^\infty \phi\left(r\right)\mathrm{d}r</math>, and <math>\int_{3\sigma}^\infty \phi\left(r\right)\mathrm{d}r</math> when <math>\sigma = \epsilon = 1.0</math>.

The value of <math>r_0</math>, for which potential energy is 0, is <math>r_0=\sigma</math>.

The force at this seperation is <math>F=\frac{24\epsilon}{\sigma}</math>.

The equilibrium separation is <math>r_{eq}=\sqrt[6]{2}\sigma</math>.

The well depth is <math>\phi\left(r_{eq}\right)=-2\epsilon</math>

<math>\int_{2\sigma}^\infty \phi\left(r\right)\mathrm{d}r=-0.0248(3sf)</math>

<math>\int_{2.5\sigma}^\infty \phi\left(r\right)\mathrm{d}r=-0.00818(3sf)</math>

<math>\int_{3\sigma}^\infty \phi\left(r\right)\mathrm{d}r=-0.00329(3sf)</math>

== ==
'''<big>TASK</big>: Estimate the number of water molecules in 1ml of water under standard conditions. Estimate the volume of <math>10000</math> water molecules under standard conditions.'''

Number of water molecules in 1mL water = 0.0556 moles (3sf) = 3.35e22 molecules (3sf)

Volume of 10000 water molecules = 2.99e-19 mL (3sf) = cube with side length of ~ 6.69 nm (3sf)

== ==
'''<big>TASK</big>: Consider an atom at position <math>\left(0.5, 0.5, 0.5\right)</math> in a cubic simulation box which runs from <math>\left(0, 0, 0\right)</math> to <math>\left(1, 1, 1\right)</math>. In a single timestep, it moves along the vector <math>\left(0.7, 0.6, 0.2\right)</math>. At what point does it end up, ''after the periodic boundary conditions have been applied''?'''.

After PBC has been applied the atom will be at (0.2,0.1,0.7)

== ==
'''<big>TASK</big>: The Lennard-Jones parameters for argon are <math>\sigma = 0.34\mathrm{nm}, \epsilon\ /\ k_B= 120 \mathrm{K}</math>. If the LJ cutoff is <math>r^* = 3.2</math>, what is it in real units? What is the well depth in <math>\mathrm{kJ\ mol}^{-1}</math>? What is the reduced temperature <math>T^* = 1.5</math> in real units?

<math>r=r^*\times\sigma=3.2\times0.34=1.088nm</math>

<math>T=T^*\frac{\epsilon}{k_B}=1.5\times120=180K</math>

<math>\phi\left(r\right)=-6.16\times10^{-24}Jmol^{-1}</math>

== ==
'''<big>TASK</big>: Why do you think giving atoms random starting coordinates causes problems in simulations? Hint: what happens if two atoms happen to be generated close together?'''

If atoms are created too close, then we will have a very large potential, meaning that if the timestep is not reduced, then we will have a large change in energy between timesteps which could mean that we do not obtain all the information from the simulation system. However if we decrease the timestep so that the change in energy is reduced to a suitable size, then it would require much more calculations to run the simulation for the same of time.

== ==
'''<big>TASK</big>: Satisfy yourself that this lattice spacing corresponds to a number density of lattice points of <math>0.8</math>. Consider instead a face-centred cubic lattice with a lattice point number density of 1.2. What is the side length of the cubic unit cell?'''

Simple cubic: 0.8 density = 1.25 cell volume = <math>\sqrt[3]{1.25}</math> = 1.07722 cell length.

Face centre cubic: 1.2 density = 0.833 cell volume = <math>\sqrt[3]{0.833}</math> = 0.941 cell length.

== ==
'''<big>TASK</big>: Consider again the face-centred cubic lattice from the previous task. How many atoms would be created by the create_atoms command if you had defined that lattice instead?'''

FCC contains 4 atoms in each lattice cell, SC contains 1 atom in each lattice cell, thus 4000 atoms would be created if FCC was used.

== ==
'''<big>TASK</big>: Using the [http://lammps.sandia.gov/doc/Section_commands.html#cmd_5 LAMMPS manual], find the purpose of the following commands in the input script:'''

<pre>
mass 1 1.0
pair_style lj/cut 3.0
pair_coeff * * 1.0 1.0
</pre>

mass 1 1.0: mass of atom type 1 is 1.0

pair_style lj/cut 3.0: cut-off for lennard jones interaction is 3 unit length

pair_coeff * * 1.0 1.0: pair_coeff refers to the force field coeffecients for pairs of atoms, * * means between any two atom types, 1.0 is the coefficient

== ==
'''<big>TASK</big>: Given that we are specifying <math>\mathbf{x}_i\left(0\right)</math> and <math>\mathbf{v}_i\left(0\right)</math>, which integration algorithm are we going to use?'''

Velocity velvet integration

== ==
'''<big>TASK</big>: Look at the lines below.'''
<pre>
### SPECIFY TIMESTEP ###
variable timestep equal 0.001
variable n_steps equal floor(100/${timestep})
variable n_steps equal floor(100/0.001)
timestep ${timestep}
timestep 0.001

### RUN SIMULATION ###
run ${n_steps}
run 100000
</pre>
'''The second line (starting "variable timestep...") tells LAMMPS that if it encounters the text ${timestep} on a subsequent line, it should replace it by the value given. In this case, the value ${timestep} is always replaced by 0.001. In light of this, what do you think the purpose of these lines is? Why not just write:'''
<pre>
timestep 0.001
run 100000
</pre>

The reason why a variable is used to hold the value of timestep is so that it is more convenient to change the value of timestep, as you would only have to change the line that defines the timestep instead of going through the whole code and changing all instances where timestep is used.

== ==
'''<big>TASK</big>: make plots of the energy, temperature, and pressure, against time for the 0.001 timestep experiment (attach a picture to your report). Does the simulation reach equilibrium? How long does this take? When you have done this, make a single plot which shows the energy versus time for all of the timesteps (again, attach a picture to your report). Choosing a timestep is a balancing act: the shorter the timestep, the more accurately the results of your simulation will reflect the physical reality; short timesteps, however, mean that the same number of simulation steps cover a shorter amount of actual time, and this is very unhelpful if the process you want to study requires observation over a long time. Of the five timesteps that you used, which is the largest to give acceptable results? Which one of the five is a ''particularly'' bad choice? Why?'''

[[File:Sh5214 liquid sim intro 001 energy.png]] [[File:Sh5214 liquid sim intro 001 temp.png]] [[File:Sh5214 liquid sim intro 001 pressure.png]] [[File:Sh5214 liquid sim intro all energy.png]] [[File:Sh5214 liquid sim intro all energy zoom.png]]

From the plot of energy vs time for timestep = 0.001 we can see that equilibrium is reached, and from the last plot it can be seen that it takes about 0.4 units of time to reach equilibrium. The largest timestep that results in the equilibration of energy is timestep = 0.01. Timestep = 0.015 is a bad choice because the system does not equilibrate and instead diverges.

== ==
'''<big>TASK</big>: We need to choose <math>\gamma</math> so that the temperature is correct <math>T = \mathfrak{T}</math> if we multiply every velocity <math>\gamma</math>. We can write two equations:'''

<math>\frac{1}{2}\sum_i m_i v_i^2 = \frac{3}{2} N k_B T</math>

<math>\frac{1}{2}\sum_i m_i \left(\gamma v_i\right)^2 = \frac{3}{2} N k_B \mathfrak{T}</math>

'''Solve these to determine <math>\gamma</math>.'''

<math>\frac{1}{2}\sum_i m_i \left(\gamma v_i\right)^2 = \frac{3}{2} N k_B </math>

<math>\implies \frac{1}{2}\sum_i m_i v_i^2=\frac{3Nk_B\mathfrak{T}}{2\gamma^2}=\frac{3}{2}Nk_BT</math>

<math>\implies \gamma^2=\frac{\mathfrak{T}}{T} \implies\gamma=\sqrt{\frac{\mathfrak{T}}{T}}</math>

== ==
'''<big>TASK</big>: Use the [http://lammps.sandia.gov/doc/fix_ave_time.html manual page] to find out the importance of the three numbers ''100 1000 100000''. How often will values of the temperature, etc., be sampled for the average? How many measurements contribute to the average? Looking to the following line, how much time will you simulate?'''

100: number of timesteps between each sample

1000: number of samples for each average

100000: number of timesteps required to obtain one average

== ==
'''<big>TASK</big>: When your simulations have finished, download the log files as before. At the end of the log file, LAMMPS will output the values and errors for the pressure, temperature, and density <math>\left(\frac{N}{V}\right)</math>. Use software of your choice to plot the density as a function of temperature for both of the pressures that you simulated. Your graph(s) should include error bars in both the x and y directions. You should also include a line corresponding to the density predicted by the ideal gas law at that pressure. Is your simulated density lower or higher? Justify this. Does the discrepancy increase or decrease with pressure?'''

[[File:Sh5214 liquid sim npt density temperature.png]]

My simulated density is lower than the density predicted by the idea gas law, this discrepancy increases with pressure. The reason why the deal gas law gives a higher density is due to the fact that it doesn't take into account the interactions between atoms, hence repulsions between atoms are not accounted for in the ideal gas law, this means that the atoms would be closer together in the idea gas law and so a higher density.

== ==
'''<big>TASK</big>: As in the last section, you need to run simulations at ten phase points. In this section, we will be in density-temperature <math>\left(\rho^*, T^*\right)</math> phase space, rather than pressure-temperature phase space. The two densities required at <math>0.2</math> and <math>0.8</math>, and the temperature range is <math>2.0, 2.2, 2.4, 2.6, 2.8</math>. Plot <math>C_V/V</math> as a function of temperature, where <math>V</math> is the volume of the simulation cell, for both of your densities (on the same graph). Is the trend the one you would expect? Attach an example of one of your input scripts to your report.'''

[[File:Sh5214 nvt heat capacity.png]]

From the graph it can be seen that specific heat capacity decreases with temperature, this is expected because there is a limit to how much energy a molecule can hold, as the temperature of the molecule increases, the amount of available rotational/vibrational energy level that the molecule can be excited into decreases, meaning that the molecule will be able to store less energy as its temperature increases.

We can also see that the specific heat capacity for the system with 0.8 density is higher than the system with 0.2 density, this is expected because a higher density means that there are more molecules per unit volume, and therefore more energy could be stored per unit volume.

<pre>
### DEFINE SIMULATION BOX GEOMETRY ###
variable D equal 0.2
lattice sc ${D}
region box block 0 15 0 15 0 15
create_box 1 box
create_atoms 1 box

### DEFINE PHYSICAL PROPERTIES OF ATOMS ###
mass 1 1.0
pair_style lj/cut/opt 3.0
pair_coeff 1 1 1.0 1.0
neighbor 2.0 bin

### SPECIFY THE REQUIRED THERMODYNAMIC STATE ###
variable T equal 2.0
variable timestep equal 0.0025

### ASSIGN ATOMIC VELOCITIES ###
velocity all create ${T} 12345 dist gaussian rot yes mom yes

### SPECIFY ENSEMBLE ###
timestep ${timestep}
fix nve all nve

### THERMODYNAMIC OUTPUT CONTROL ###
thermo_style custom time etotal temp press
thermo 10

### RECORD TRAJECTORY ###
dump traj all custom 1000 output-1 id x y z

### SPECIFY TIMESTEP ###

### RUN SIMULATION TO MELT CRYSTAL ###
run 10000
unfix nve
reset_timestep 0

### BRING SYSTEM TO REQUIRED STATE ###
variable tdamp equal ${timestep}*100
variable vdamp equal ${timestep}*1000
fix nvt all nvt temp ${T} ${T} ${tdamp}
run 10000
reset_timestep 0

## SWITCH OFF THERMOSTAT ##
unfix nvt
fix nve all nve

### MEASURE SYSTEM STATE ###
thermo_style custom step etotal temp press density atoms vol
variable vol equal vol
variable N2 equal atoms*atoms
variable dens equal density
variable dens2 equal density*density
variable temp equal temp
variable temp2 equal temp*temp
variable press equal press
variable press2 equal press*press
variable etotal equal etotal
variable etotal2 equal etotal*etotal
fix aves all ave/time 100 1000 100000 v_dens v_temp v_press v_etotal v_dens2 v_temp2 v_press2 v_etotal2 v_vol
run 100000

variable avetemp equal f_aves[2]
variable errtemp equal sqrt(f_aves[6]-f_aves[2]*f_aves[2])
variable cv equal (${N2}*(f_aves[8]-f_aves[4]*f_aves[4])/(${T}*${T}))
variable V equal f_aves[9]

print "Averages"
print "--------"
print "Temperature: ${avetemp}"
print "Stderr: ${errtemp}"
print "cv: ${cv}"
print "volume: ${V}"
</pre>

== ==
'''<big>TASK</big>: perform simulations of the Lennard-Jones system in the three phases. When each is complete, download the trajectory and calculate <math>g(r)</math> and <math>\int g(r)\mathrm{d}r</math>. Plot the RDFs for the three systems on the same axes, and attach a copy to your report. Discuss qualitatively the differences between the three RDFs, and what this tells you about the structure of the system in each phase. In the solid case, illustrate which lattice sites the first three peaks correspond to. What is the lattice spacing? What is the coordination number for each of the first three peaks?'''

[[File:Sh5214 liquid sim rdf plot.png]]

Firstly the range of peaks for solids is longer than liquids which is longer than gas, this means that in the system of solids, the atoms are more ordered than the atoms in the system of liquids, which are more ordered than the atoms in the system of gases. Another observation that could be made is that the magnitude of the peaks decreases as r increases, this is because g(r) is calculated by number of atoms in the shell divided by volume of shell and the density of the system, as r increases, the volume of the shell will increase at a faster rate than the number of atoms in the shell, and thus g(r) decreases. Also it can be seen that magnitude of the peaks of the RDF of the solid is greater than both liquid and gas, this is because solids are more closely packed and more ordered than both gas and liquid, hence at certain r values there will be a high number of atoms and at other r values there will be a very low number of atoms, and that's why solids have higher peaks and lower troughs.

The first peak for the solid RDF is at r=1.025 and the second peak is at r=1.425, thus the lattice spacing is 0.4 times the interatomic separation. There are three different planes in a fcc crystal, they have miller indices of (1,1,1) (1,1,0) and (1,0,0). The (1,1,1) plane corresponds to the first peak and has a coordination number of 6 atoms, the (1,1,0) plane corresponds to the second peak and has a coordination number of 2, the (1,0,0) plane corresponds to the third peak and has a coordination number of 4. The coordination number of lattice sites will affect the magnitude of the peak, and that's why magnitude of first peak > third peak > second peak

== ==
'''<big>TASK</big>: make a plot for each of your simulations (solid, liquid, and gas), showing the mean squared displacement (the "total" MSD) as a function of timestep. Are these as you would expect? Estimate <math>D</math> in each case. Be careful with the units! Repeat this procedure for the MSD data that you were given from the one million atom simulations.'''

[[File:Sh5214 liquid sim d gas msd.png]] [[File:Sh5214 liquid sim d liquid msd.png]] [[File:Sh5214 liquid sim d solid msd.png]]

Gas MSD (my data):
<pre>
t=2000, MSD=29.2
t=5000, MSD=104.7
gradient = (104.7-29.2)/3000 = 0.0252 (3sf)
D = gradient/6 = 0.00419 (3sf)
</pre>

Gas MSD (1e6 atoms):
<pre>
t=2000, MSD=36.4
t=5000, MSD=144.4
gradient = (144.4-36.4)/3000 = 0.036
D = gradient/6 = 0.006 (3sf)
</pre>

Liquid MSD (my data):
<pre>
t=0, MSD=0
t=5000, MSD=11
gradient = 11/5000 = 0.0022
D = gradient/6 = 0.000367 (3sf)
</pre>

Liquid MSD (1e6 atoms):
<pre>
t=0, MSD=0
t=5000, MSD=5.19
gradient = 5.19/5000 = 0.00104 (3sf)
D = gradient/6 = 0.000173
</pre>

Solid MSD (my data): D=0

Solid MSD (1e6 atoms): D=0

MSD is the average distance traveled by a particle in a system, it can be seen that the MSD of the solid plateaus very early on and at a very low value this is due to the fact that atoms in a solid are rigid and fixed in place therefore average distance traveled by a atom in a solid is very small. We see that the MSD plot of the gas starts off as a parabola before turning linear, the parabola part corresponds to when gas atom doesn't collide with anything else meaning that during this time the velocity of the atom is constant and so distance increases proportional to time squared, when more and more collisions occur the curve goes from parabolic to linear. In the liquid MSD plot, we can see that at the start of the curve for a short amount of time the curve is not linear, this is because little to no collisions occur right at the start of the simulation, for most of the plot the line is linear due to the fact that there are collisions between liquid atoms.

The magnitude of MSD for gas > liquid > solid, this is due to the fact that density of gas < liquid < solid, a lower density means that less collisions occur and thus the average translational velocity of atoms is higher, resulting in greater distance traveled.

== ==
'''<big>TASK</big>: In the theoretical section at the beginning, the equation for the evolution of the position of a 1D harmonic oscillator as a function of time was given. Using this, evaluate the normalised velocity autocorrelation function for a 1D harmonic oscillator (it is analytic!):

<math>C\left(\tau\right) = \frac{\int_{-\infty}^{\infty} v\left(t\right)v\left(t + \tau\right)\mathrm{d}t}{\int_{-\infty}^{\infty} v^2\left(t\right)\mathrm{d}t}</math>

'''Be sure to show your working in your writeup. On the same graph, with x range 0 to 500, plot <math>C\left(\tau\right)</math> with <math>\omega = 1/2\pi</math> and the VACFs from your liquid and solid simulations. What do the minima in the VACFs for the liquid and solid system represent? Discuss the origin of the differences between the liquid and solid VACFs. The harmonic oscillator VACF is very different to the Lennard Jones solid and liquid. Why is this? Attach a copy of your plot to your writeup.'''

<math>x(t)=Acos(\omega t+\phi)</math>

<math>v(t)=\frac{dx}{dt}=-A\omega sin(\omega t+\phi)</math>

<math>C(\tau)=\frac{\int_{-\infty}^{\infty}-A\omega sin(\omega t+\phi).-A\omega sin(\omega t+\omega\tau+\phi)\mathrm{d}t}{\int_{-\infty}^{\infty}A^2\omega^2sin^2(\omega t+\phi)\mathrm{d}t}</math>

<math>=\frac{\int_{-\infty}^{\infty}0.5A^2\omega^2[cos(\omega\tau)-cos(2\omega t+\omega\tau+2\phi)]dt}{\int_{-\infty}^{\infty}0.5A^2\omega^2-0.5A^2\omega^2cos(2\omega t+2\phi)}</math>

<math>=\frac{[0.5A^2\omega^2cos(\omega\tau)t+0.25A^2\omega sin(2\omega t+\omega\tau+2\phi)]_{-\infty}^\infty}{[0.5A^2\omega^2t+0.25A^2\omega cos(2\omega t+2\phi)]_{-\infty}^\infty}</math>

We know that sine and cosine has a range of -1 to 1, hense the above equation can be simplified to:

<math>\approx\frac{[0.5A^2\omega^2cos(\omega\tau)t]_{-\infty}^\infty}{[0.5A^2\omega^2t]_{-\infty}^\infty}</math>

<math>\approx cos(\omega\tau)</math>

[[File:Sh5214 liquide sim d solid liquid cos vacf normalised.png]]

The minima in the VACFs for the liquid and solid system represents the time at which the velocity of the atoms in the system is at a maximum. We can see that the VACF for solids oscillates about zero, with the amplitude of the oscillation will decrease over time, this is due to the fact that the atoms in the solids are locked up in a lattice, these atoms in the lattice will want to be at the optimal position where the repulsive forces and attractive forces balance out, they will therefore oscillate about this position which will also result in an oscillation in velocity, the reason why the magnitude of these oscillation decreases over time is because of other forces that perturb the oscillation. We can see that there is also a very slight minima on the liquid VACF curve, which shows that liquids also oscillate but the oscillation disappears very quickly, the reason why atoms in liquids don't oscillate for as long as solids is because liquid atoms are not locked in a lattice like solids and so are free to move, this means that any oscillation is very quickly dampened out by diffusion.

The reason why the harmonic oscillator VACF is so different is due to the fact that it doesn't take into account dampening of the oscillation, hence the magnitude of the oscillation never decreases.

== ==
'''<big>TASK</big>: Use the trapezium rule to approximate the integral under the velocity autocorrelation function for the solid, liquid, and gas, and use these values to estimate <math>D</math> in each case. You should make a plot of the running integral in each case. Are they as you expect? Repeat this procedure for the VACF data that you were given from the one million atom simulations. What do you think is the largest source of error in your estimates of D from the VACF?'''

[[File:Sh5214 liquid sim d gas integral vacf.png]] [[File:Sh5214 liquid sim d liquid integral vacf.png]] [[File:Sh5214 liquid sim d solid integral vacf.png]]

Gas VACF (my data):
<pre>
t=5000, Running Integral=6.54 (3sf)
D = Running Integral/3 = 2.18 (3sf)
</pre>

Gas VACF (1e6 atoms):
<pre>
t=5000, Running Integral=9.81 (3sf)
D = Running Integral/3 = 3.27 (3sf)
</pre>

Liquid VACF (my data):
<pre>
t=5000, Running Integral=0.569 (3sf)
D = Running Integral/3 = 0.190 (3sf)
</pre>

Liquid VACF (1e6 atoms):
<pre>
t=5000, Running Integral=0.270 (3sf)
D = Running Integral/3 = 0.0901 (3sf)
</pre>

Solid VACF (my data): D = 0

Solid VACF (1e6 atoms): D = 0

Rep:Mod:sh5214 liquid simulation

2017-02-24T07:57:45Z

Sh5214: /* ? */

== ==
'''<big>TASK</big>: Open the file HO.xls. In it, the velocity-Verlet algorithm is used to model the behaviour of a classical harmonic oscillator. Complete the three columns "ANALYTICAL", "ERROR", and "ENERGY": "ANALYTICAL" should contain the value of the classical solution for the position at time <math>t</math>, "ERROR" should contain the ''absolute'' difference between "ANALYTICAL" and the velocity-Verlet solution (i.e. ERROR should always be positive -- make sure you leave the half step rows blank!), and "ENERGY" should contain the total energy of the oscillator for the velocity-Verlet solution. Remember that the position of a classical harmonic oscillator is given by <math> x\left(t\right) = A\cos\left(\omega t + \phi\right)</math> (the values of <math>A</math>, <math>\omega</math>, and <math>\phi</math> are worked out for you in the sheet).'''

[[File:Sh5214 liquid sim theory analytical.png]] [[File:Sh5214 liquid sim theory energy.png]] [[File:Sh5214 liquid sim theory error.png ‎]]

== ==
'''<big>TASK</big>: For the default timestep value, 0.1, estimate the positions of the maxima in the ERROR column as a function of time. Make a plot showing these values as a function of time, and fit an appropriate function to the data.'''

[[File:Sh5214 liquid sim theory error maxima.png ]]

== ==
'''<big>TASK</big>: Experiment with different values of the timestep. What sort of a timestep do you need to use to ensure that the total energy does not change by more than 1% over the course of your "simulation"? Why do you think it is important to monitor the total energy of a physical system when modelling its behaviour numerically?'''

The maximum timestep that doesn't result in more than an 1% change in energy is around 0.45. The reason why the energy is monitored is to make sure the energy change between each timestep is not too large. If there is a large energy difference, it would mean that the overall change in velocity of the particles in the system is too large, which would lead to not obtaining the full amount of information from the simulation.

== ==
'''<big>TASK:</big> For a single Lennard-Jones interaction, <math>\phi\left(r\right) = 4\epsilon \left( \frac{\sigma^{12}}{r^{12}} - \frac{\sigma^6}{r^6} \right)</math>, find the separation, <math>r_0</math>, at which the potential energy is zero. What is the force at this separation? Find the equilibrium separation, <math>r_{eq}</math>, and work out the well depth (<math>\phi\left(r_{eq}\right)</math>). Evaluate the integrals <math>\int_{2\sigma}^\infty \phi\left(r\right)\mathrm{d}r</math>, <math>\int_{2.5\sigma}^\infty \phi\left(r\right)\mathrm{d}r</math>, and <math>\int_{3\sigma}^\infty \phi\left(r\right)\mathrm{d}r</math> when <math>\sigma = \epsilon = 1.0</math>.

The value of <math>r_0</math>, for which potential energy is 0, is <math>r_0=\sigma</math>.

The force at this seperation is <math>F=\frac{24\epsilon}{\sigma}</math>.

The equilibrium separation is <math>r_{eq}=\sqrt[6]{2}\sigma</math>.

The well depth is <math>\phi\left(r_{eq}\right)=-2\epsilon</math>

<math>\int_{2\sigma}^\infty \phi\left(r\right)\mathrm{d}r=-0.0248(3sf)</math>

<math>\int_{2.5\sigma}^\infty \phi\left(r\right)\mathrm{d}r=-0.00818(3sf)</math>

<math>\int_{3\sigma}^\infty \phi\left(r\right)\mathrm{d}r=-0.00329(3sf)</math>

== ==
'''<big>TASK</big>: Estimate the number of water molecules in 1ml of water under standard conditions. Estimate the volume of <math>10000</math> water molecules under standard conditions.'''

Number of water molecules in 1mL water = 0.0556 moles (3sf) = 3.35e22 molecules (3sf)

Volume of 10000 water molecules = 2.99e-19 mL (3sf) = cube with side length of ~ 6.69 nm (3sf)

== ==
'''<big>TASK</big>: Consider an atom at position <math>\left(0.5, 0.5, 0.5\right)</math> in a cubic simulation box which runs from <math>\left(0, 0, 0\right)</math> to <math>\left(1, 1, 1\right)</math>. In a single timestep, it moves along the vector <math>\left(0.7, 0.6, 0.2\right)</math>. At what point does it end up, ''after the periodic boundary conditions have been applied''?'''.

After PBC has been applied the atom will be at (0.2,0.1,0.7)

== ==
'''<big>TASK</big>: The Lennard-Jones parameters for argon are <math>\sigma = 0.34\mathrm{nm}, \epsilon\ /\ k_B= 120 \mathrm{K}</math>. If the LJ cutoff is <math>r^* = 3.2</math>, what is it in real units? What is the well depth in <math>\mathrm{kJ\ mol}^{-1}</math>? What is the reduced temperature <math>T^* = 1.5</math> in real units?

<math>r=r^*\times\sigma=3.2\times0.34=1.088nm</math>

<math>T=T^*\frac{\epsilon}{k_B}=1.5\times120=180K</math>

<math>\phi\left(r\right)=-6.16\times10^{-24}Jmol^{-1}</math>

== ==
'''<big>TASK</big>: Why do you think giving atoms random starting coordinates causes problems in simulations? Hint: what happens if two atoms happen to be generated close together?'''

If atoms are created too close, then we will have a very large potential, meaning that if the timestep is not reduced, then we will have a large change in energy between timesteps which could mean that we do not obtain all the information from the simulation system. However if we decrease the timestep so that the change in energy is reduced to a suitable size, then it would require much more calculations to run the simulation for the same of time.

== ==
'''<big>TASK</big>: Satisfy yourself that this lattice spacing corresponds to a number density of lattice points of <math>0.8</math>. Consider instead a face-centred cubic lattice with a lattice point number density of 1.2. What is the side length of the cubic unit cell?'''

Simple cubic: 0.8 density = 1.25 cell volume = <math>\sqrt[3]{1.25}</math> = 1.07722 cell length.

Face centre cubic: 1.2 density = 0.833 cell volume = <math>\sqrt[3]{0.833}</math> = 0.941 cell length.

== ==
'''<big>TASK</big>: Consider again the face-centred cubic lattice from the previous task. How many atoms would be created by the create_atoms command if you had defined that lattice instead?'''

FCC contains 4 atoms in each lattice cell, SC contains 1 atom in each lattice cell, thus 4000 atoms would be created if FCC was used.

== ==
'''<big>TASK</big>: Using the [http://lammps.sandia.gov/doc/Section_commands.html#cmd_5 LAMMPS manual], find the purpose of the following commands in the input script:'''

<pre>
mass 1 1.0
pair_style lj/cut 3.0
pair_coeff * * 1.0 1.0
</pre>

mass 1 1.0: mass of atom type 1 is 1.0

pair_style lj/cut 3.0: cut-off for lennard jones interaction is 3 unit length

pair_coeff * * 1.0 1.0: pair_coeff refers to the force field coeffecients for pairs of atoms, * * means between any two atom types, 1.0 is the coefficient

== ==
'''<big>TASK</big>: Given that we are specifying <math>\mathbf{x}_i\left(0\right)</math> and <math>\mathbf{v}_i\left(0\right)</math>, which integration algorithm are we going to use?'''

Velocity velvet integration

== ==
'''<big>TASK</big>: Look at the lines below.'''
<pre>
### SPECIFY TIMESTEP ###
variable timestep equal 0.001
variable n_steps equal floor(100/${timestep})
variable n_steps equal floor(100/0.001)
timestep ${timestep}
timestep 0.001

### RUN SIMULATION ###
run ${n_steps}
run 100000
</pre>
'''The second line (starting "variable timestep...") tells LAMMPS that if it encounters the text ${timestep} on a subsequent line, it should replace it by the value given. In this case, the value ${timestep} is always replaced by 0.001. In light of this, what do you think the purpose of these lines is? Why not just write:'''
<pre>
timestep 0.001
run 100000
</pre>

The reason why a variable is used to hold the value of timestep is so that it is more convenient to change the value of timestep, as you would only have to change the line that defines the timestep instead of going through the whole code and changing all instances where timestep is used.

== ==
'''<big>TASK</big>: make plots of the energy, temperature, and pressure, against time for the 0.001 timestep experiment (attach a picture to your report). Does the simulation reach equilibrium? How long does this take? When you have done this, make a single plot which shows the energy versus time for all of the timesteps (again, attach a picture to your report). Choosing a timestep is a balancing act: the shorter the timestep, the more accurately the results of your simulation will reflect the physical reality; short timesteps, however, mean that the same number of simulation steps cover a shorter amount of actual time, and this is very unhelpful if the process you want to study requires observation over a long time. Of the five timesteps that you used, which is the largest to give acceptable results? Which one of the five is a ''particularly'' bad choice? Why?'''

[[File:Sh5214 liquid sim intro 001 energy.png]] [[File:Sh5214 liquid sim intro 001 temp.png]] [[File:Sh5214 liquid sim intro 001 pressure.png]] [[File:Sh5214 liquid sim intro all energy.png]] [[File:Sh5214 liquid sim intro all energy zoom.png]]

From the plot of energy vs time for timestep = 0.001 we can see that equilibrium is reached, and from the last plot it can be seen that it takes about 0.4 units of time to reach equilibrium. The largest timestep that results in the equilibration of energy is timestep = 0.01. Timestep = 0.015 is a bad choice because the system does not equilibrate and instead diverges.

== ==
'''<big>TASK</big>: We need to choose <math>\gamma</math> so that the temperature is correct <math>T = \mathfrak{T}</math> if we multiply every velocity <math>\gamma</math>. We can write two equations:'''

<math>\frac{1}{2}\sum_i m_i v_i^2 = \frac{3}{2} N k_B T</math>

<math>\frac{1}{2}\sum_i m_i \left(\gamma v_i\right)^2 = \frac{3}{2} N k_B \mathfrak{T}</math>

'''Solve these to determine <math>\gamma</math>.'''

<math>\frac{1}{2}\sum_i m_i \left(\gamma v_i\right)^2 = \frac{3}{2} N k_B </math>

<math>\implies \frac{1}{2}\sum_i m_i v_i^2=\frac{3Nk_B\mathfrak{T}}{2\gamma^2}=\frac{3}{2}Nk_BT</math>

<math>\implies \gamma^2=\frac{\mathfrak{T}}{T} \implies\gamma=\sqrt{\frac{\mathfrak{T}}{T}}</math>

== ==
'''<big>TASK</big>: Use the [http://lammps.sandia.gov/doc/fix_ave_time.html manual page] to find out the importance of the three numbers ''100 1000 100000''. How often will values of the temperature, etc., be sampled for the average? How many measurements contribute to the average? Looking to the following line, how much time will you simulate?'''

100: number of timesteps between each sample

1000: number of samples for each average

100000: number of timesteps required to obtain one average

== ==
'''<big>TASK</big>: When your simulations have finished, download the log files as before. At the end of the log file, LAMMPS will output the values and errors for the pressure, temperature, and density <math>\left(\frac{N}{V}\right)</math>. Use software of your choice to plot the density as a function of temperature for both of the pressures that you simulated. Your graph(s) should include error bars in both the x and y directions. You should also include a line corresponding to the density predicted by the ideal gas law at that pressure. Is your simulated density lower or higher? Justify this. Does the discrepancy increase or decrease with pressure?'''

[[File:Sh5214 liquid sim npt density temperature.png]]

My simulated density is lower than the density predicted by the idea gas law, this discrepancy increases with pressure. The reason why the deal gas law gives a higher density is due to the fact that it doesn't take into account the interactions between atoms, hence repulsions between atoms are not accounted for in the ideal gas law, this means that the atoms would be closer together in the idea gas law and so a higher density.

== ==
'''<big>TASK</big>: As in the last section, you need to run simulations at ten phase points. In this section, we will be in density-temperature <math>\left(\rho^*, T^*\right)</math> phase space, rather than pressure-temperature phase space. The two densities required at <math>0.2</math> and <math>0.8</math>, and the temperature range is <math>2.0, 2.2, 2.4, 2.6, 2.8</math>. Plot <math>C_V/V</math> as a function of temperature, where <math>V</math> is the volume of the simulation cell, for both of your densities (on the same graph). Is the trend the one you would expect? Attach an example of one of your input scripts to your report.'''

[[File:Sh5214 nvt heat capacity.png]]

From the graph it can be seen that specific heat capacity decreases with temperature, this is expected because there is a limit to how much energy a molecule can hold, as the temperature of the molecule increases, the amount of available rotational/vibrational energy level that the molecule can be excited into decreases, meaning that the molecule will be able to store less energy as its temperature increases.

We can also see that the specific heat capacity for the system with 0.8 density is higher than the system with 0.2 density, this is expected because a higher density means that there are more molecules per unit volume, and therefore more energy could be stored per unit volume.

<pre>
### DEFINE SIMULATION BOX GEOMETRY ###
lattice sc 0.2
region box block 0 15 0 15 0 15
create_box 1 box
create_atoms 1 box

### DEFINE PHYSICAL PROPERTIES OF ATOMS ###
mass 1 1.0
pair_style lj/cut/opt 3.0
pair_coeff 1 1 1.0 1.0
neighbor 2.0 bin

### SPECIFY THE REQUIRED THERMODYNAMIC STATE ###
variable T equal 2.0
variable timestep equal 0.0025

### ASSIGN ATOMIC VELOCITIES ###
velocity all create ${T} 12345 dist gaussian rot yes mom yes

### SPECIFY ENSEMBLE ###
timestep ${timestep}
fix nve all nve

### THERMODYNAMIC OUTPUT CONTROL ###
thermo_style custom time etotal temp press
thermo 10

### RECORD TRAJECTORY ###
dump traj all custom 1000 output-1 id x y z

### SPECIFY TIMESTEP ###

### RUN SIMULATION TO MELT CRYSTAL ###
run 10000
unfix nve
reset_timestep 0

### BRING SYSTEM TO REQUIRED STATE ###
variable tdamp equal ${timestep}*100
variable vdamp equal ${timestep}*1000
fix nvt all nvt temp ${T} ${T} ${tdamp}
run 10000
reset_timestep 0

## SWITCH OFF THERMOSTAT ##
unfix nvt
fix nve all nve

### MEASURE SYSTEM STATE ###
thermo_style custom step etotal temp press density atoms vol
variable vol equal vol
variable N2 equal atoms*atoms
variable dens equal density
variable dens2 equal density*density
variable temp equal temp
variable temp2 equal temp*temp
variable press equal press
variable press2 equal press*press
variable etotal equal etotal
variable etotal2 equal etotal*etotal
fix aves all ave/time 100 1000 100000 v_dens v_temp v_press v_etotal v_dens2 v_temp2 v_press2 v_etotal2 v_vol
run 100000

variable avetemp equal f_aves[2]
variable errtemp equal sqrt(f_aves[6]-f_aves[2]*f_aves[2])
variable cv equal (${N2}*(f_aves[8]-f_aves[4]*f_aves[4])/(${T}*${T}))
variable V equal f_aves[9]

print "Averages"
print "--------"
print "Temperature: ${avetemp}"
print "Stderr: ${errtemp}"
print "cv: ${cv}"
print "volume: ${V}"
</pre>

== ==
'''<big>TASK</big>: perform simulations of the Lennard-Jones system in the three phases. When each is complete, download the trajectory and calculate <math>g(r)</math> and <math>\int g(r)\mathrm{d}r</math>. Plot the RDFs for the three systems on the same axes, and attach a copy to your report. Discuss qualitatively the differences between the three RDFs, and what this tells you about the structure of the system in each phase. In the solid case, illustrate which lattice sites the first three peaks correspond to. What is the lattice spacing? What is the coordination number for each of the first three peaks?'''

[[File:Sh5214 liquid sim rdf plot.png]]

Firstly the range of peaks for solids is longer than liquids which is longer than gas, this means that in the system of solids, the atoms are more ordered than the atoms in the system of liquids, which are more ordered than the atoms in the system of gases. Another observation that could be made is that the magnitude of the peaks decreases as r increases, this is because g(r) is calculated by number of atoms in the shell divided by volume of shell and the density of the system, as r increases, the volume of the shell will increase at a faster rate than the number of atoms in the shell, and thus g(r) decreases. Also it can be seen that magnitude of the peaks of the RDF of the solid is greater than both liquid and gas, this is because solids are more closely packed and more ordered than both gas and liquid, hence at certain r values there will be a high number of atoms and at other r values there will be a very low number of atoms, and that's why solids have higher peaks and lower troughs.

The first peak for the solid RDF is at r=1.025 and the second peak is at r=1.425, thus the lattice spacing is 0.4 times the interatomic separation. There are three different planes in a fcc crystal, they have miller indices of (1,1,1) (1,1,0) and (1,0,0). The (1,1,1) plane corresponds to the first peak and has a coordination number of 6 atoms, the (1,1,0) plane corresponds to the second peak and has a coordination number of 2, the (1,0,0) plane corresponds to the third peak and has a coordination number of 4. The coordination number of lattice sites will affect the magnitude of the peak, and that's why magnitude of first peak > third peak > second peak

== ==
'''<big>TASK</big>: make a plot for each of your simulations (solid, liquid, and gas), showing the mean squared displacement (the "total" MSD) as a function of timestep. Are these as you would expect? Estimate <math>D</math> in each case. Be careful with the units! Repeat this procedure for the MSD data that you were given from the one million atom simulations.'''

[[File:Sh5214 liquid sim d gas msd.png]] [[File:Sh5214 liquid sim d liquid msd.png]] [[File:Sh5214 liquid sim d solid msd.png]]

Gas MSD (my data):
<pre>
t=2000, MSD=29.2
t=5000, MSD=104.7
gradient = (104.7-29.2)/3000 = 0.0252 (3sf)
D = gradient/6 = 0.00419 (3sf)
</pre>

Gas MSD (1e6 atoms):
<pre>
t=2000, MSD=36.4
t=5000, MSD=144.4
gradient = (144.4-36.4)/3000 = 0.036
D = gradient/6 = 0.006 (3sf)
</pre>

Liquid MSD (my data):
<pre>
t=0, MSD=0
t=5000, MSD=11
gradient = 11/5000 = 0.0022
D = gradient/6 = 0.000367 (3sf)
</pre>

Liquid MSD (1e6 atoms):
<pre>
t=0, MSD=0
t=5000, MSD=5.19
gradient = 5.19/5000 = 0.00104 (3sf)
D = gradient/6 = 0.000173
</pre>

Solid MSD (my data): D=0

Solid MSD (1e6 atoms): D=0

MSD is the average distance traveled by a particle in a system, it can be seen that the MSD of the solid plateaus very early on and at a very low value this is due to the fact that atoms in a solid are rigid and fixed in place therefore average distance traveled by a atom in a solid is very small. We see that the MSD plot of the gas starts off as a parabola before turning linear, the parabola part corresponds to when gas atom doesn't collide with anything else meaning that during this time the velocity of the atom is constant and so distance increases proportional to time squared, when more and more collisions occur the curve goes from parabolic to linear. In the liquid MSD plot, we can see that at the start of the curve for a short amount of time the curve is not linear, this is because little to no collisions occur right at the start of the simulation, for most of the plot the line is linear due to the fact that there are collisions between liquid atoms.

The magnitude of MSD for gas > liquid > solid, this is due to the fact that density of gas < liquid < solid, a lower density means that less collisions occur and thus the average translational velocity of atoms is higher, resulting in greater distance traveled.

== ==
'''<big>TASK</big>: In the theoretical section at the beginning, the equation for the evolution of the position of a 1D harmonic oscillator as a function of time was given. Using this, evaluate the normalised velocity autocorrelation function for a 1D harmonic oscillator (it is analytic!):

<math>C\left(\tau\right) = \frac{\int_{-\infty}^{\infty} v\left(t\right)v\left(t + \tau\right)\mathrm{d}t}{\int_{-\infty}^{\infty} v^2\left(t\right)\mathrm{d}t}</math>

'''Be sure to show your working in your writeup. On the same graph, with x range 0 to 500, plot <math>C\left(\tau\right)</math> with <math>\omega = 1/2\pi</math> and the VACFs from your liquid and solid simulations. What do the minima in the VACFs for the liquid and solid system represent? Discuss the origin of the differences between the liquid and solid VACFs. The harmonic oscillator VACF is very different to the Lennard Jones solid and liquid. Why is this? Attach a copy of your plot to your writeup.'''

<math>x(t)=Acos(\omega t+\phi)</math>

<math>v(t)=\frac{dx}{dt}=-A\omega sin(\omega t+\phi)</math>

<math>C(\tau)=\frac{\int_{-\infty}^{\infty}-A\omega sin(\omega t+\phi).-A\omega sin(\omega t+\omega\tau+\phi)\mathrm{d}t}{\int_{-\infty}^{\infty}A^2\omega^2sin^2(\omega t+\phi)\mathrm{d}t}</math>

<math>=\frac{\int_{-\infty}^{\infty}0.5A^2\omega^2[cos(\omega\tau)-cos(2\omega t+\omega\tau+2\phi)]dt}{\int_{-\infty}^{\infty}0.5A^2\omega^2-0.5A^2\omega^2cos(2\omega t+2\phi)}</math>

<math>=\frac{[0.5A^2\omega^2cos(\omega\tau)t+0.25A^2\omega sin(2\omega t+\omega\tau+2\phi)]_{-\infty}^\infty}{[0.5A^2\omega^2t+0.25A^2\omega cos(2\omega t+2\phi)]_{-\infty}^\infty}</math>

We know that sine and cosine has a range of -1 to 1, hense the above equation can be simplified to:

<math>\approx\frac{[0.5A^2\omega^2cos(\omega\tau)t]_{-\infty}^\infty}{[0.5A^2\omega^2t]_{-\infty}^\infty}</math>

<math>\approx cos(\omega\tau)</math>

[[File:Sh5214 liquide sim d solid liquid cos vacf normalised.png]]

The minima in the VACFs for the liquid and solid system represents the time at which the velocity of the atoms in the system is at a maximum. We can see that the VACF for solids oscillates about zero, with the amplitude of the oscillation will decrease over time, this is due to the fact that the atoms in the solids are locked up in a lattice, these atoms in the lattice will want to be at the optimal position where the repulsive forces and attractive forces balance out, they will therefore oscillate about this position which will also result in an oscillation in velocity, the reason why the magnitude of these oscillation decreases over time is because of other forces that perturb the oscillation. We can see that there is also a very slight minima on the liquid VACF curve, which shows that liquids also oscillate but the oscillation disappears very quickly, the reason why atoms in liquids don't oscillate for as long as solids is because liquid atoms are not locked in a lattice like solids and so are free to move, this means that any oscillation is very quickly dampened out by diffusion.

The reason why the harmonic oscillator VACF is so different is due to the fact that it doesn't take into account dampening of the oscillation, hence the magnitude of the oscillation never decreases.

== ==
'''<big>TASK</big>: Use the trapezium rule to approximate the integral under the velocity autocorrelation function for the solid, liquid, and gas, and use these values to estimate <math>D</math> in each case. You should make a plot of the running integral in each case. Are they as you expect? Repeat this procedure for the VACF data that you were given from the one million atom simulations. What do you think is the largest source of error in your estimates of D from the VACF?'''

[[File:Sh5214 liquid sim d gas integral vacf.png]] [[File:Sh5214 liquid sim d liquid integral vacf.png]] [[File:Sh5214 liquid sim d solid integral vacf.png]]

Gas VACF (my data):
<pre>
t=5000, Running Integral=6.54 (3sf)
D = Running Integral/3 = 2.18 (3sf)
</pre>

Gas VACF (1e6 atoms):
<pre>
t=5000, Running Integral=9.81 (3sf)
D = Running Integral/3 = 3.27 (3sf)
</pre>

Liquid VACF (my data):
<pre>
t=5000, Running Integral=0.569 (3sf)
D = Running Integral/3 = 0.190 (3sf)
</pre>

Liquid VACF (1e6 atoms):
<pre>
t=5000, Running Integral=0.270 (3sf)
D = Running Integral/3 = 0.0901 (3sf)
</pre>

Solid VACF (my data): D = 0

Solid VACF (1e6 atoms): D = 0

Rep:Mod:sh5214 liquid simulation

2017-02-24T07:57:12Z

Sh5214: /* ? */

== ==
'''<big>TASK</big>: Open the file HO.xls. In it, the velocity-Verlet algorithm is used to model the behaviour of a classical harmonic oscillator. Complete the three columns "ANALYTICAL", "ERROR", and "ENERGY": "ANALYTICAL" should contain the value of the classical solution for the position at time <math>t</math>, "ERROR" should contain the ''absolute'' difference between "ANALYTICAL" and the velocity-Verlet solution (i.e. ERROR should always be positive -- make sure you leave the half step rows blank!), and "ENERGY" should contain the total energy of the oscillator for the velocity-Verlet solution. Remember that the position of a classical harmonic oscillator is given by <math> x\left(t\right) = A\cos\left(\omega t + \phi\right)</math> (the values of <math>A</math>, <math>\omega</math>, and <math>\phi</math> are worked out for you in the sheet).'''

[[File:Sh5214 liquid sim theory analytical.png]] [[File:Sh5214 liquid sim theory energy.png]] [[File:Sh5214 liquid sim theory error.png ‎]]

== ==
'''<big>TASK</big>: For the default timestep value, 0.1, estimate the positions of the maxima in the ERROR column as a function of time. Make a plot showing these values as a function of time, and fit an appropriate function to the data.'''

[[File:Sh5214 liquid sim theory error maxima.png ]]

== ==
'''<big>TASK</big>: Experiment with different values of the timestep. What sort of a timestep do you need to use to ensure that the total energy does not change by more than 1% over the course of your "simulation"? Why do you think it is important to monitor the total energy of a physical system when modelling its behaviour numerically?'''

The maximum timestep that doesn't result in more than an 1% change in energy is around 0.45. The reason why the energy is monitored is to make sure the energy change between each timestep is not too large. If there is a large energy difference, it would mean that the overall change in velocity of the particles in the system is too large, which would lead to not obtaining the full amount of information from the simulation.

== ==
'''<big>TASK:</big> For a single Lennard-Jones interaction, <math>\phi\left(r\right) = 4\epsilon \left( \frac{\sigma^{12}}{r^{12}} - \frac{\sigma^6}{r^6} \right)</math>, find the separation, <math>r_0</math>, at which the potential energy is zero. What is the force at this separation? Find the equilibrium separation, <math>r_{eq}</math>, and work out the well depth (<math>\phi\left(r_{eq}\right)</math>). Evaluate the integrals <math>\int_{2\sigma}^\infty \phi\left(r\right)\mathrm{d}r</math>, <math>\int_{2.5\sigma}^\infty \phi\left(r\right)\mathrm{d}r</math>, and <math>\int_{3\sigma}^\infty \phi\left(r\right)\mathrm{d}r</math> when <math>\sigma = \epsilon = 1.0</math>.

The value of <math>r_0</math>, for which potential energy is 0, is <math>r_0=\sigma</math>.

The force at this seperation is <math>F=\frac{24\epsilon}{\sigma}</math>.

The equilibrium separation is <math>r_{eq}=\sqrt[6]{2}\sigma</math>.

The well depth is <math>\phi\left(r_{eq}\right)=-2\epsilon</math>

<math>\int_{2\sigma}^\infty \phi\left(r\right)\mathrm{d}r=-0.0248(3sf)</math>

<math>\int_{2.5\sigma}^\infty \phi\left(r\right)\mathrm{d}r=-0.00818(3sf)</math>

<math>\int_{3\sigma}^\infty \phi\left(r\right)\mathrm{d}r=-0.00329(3sf)</math>

== ==
'''<big>TASK</big>: Estimate the number of water molecules in 1ml of water under standard conditions. Estimate the volume of <math>10000</math> water molecules under standard conditions.'''

Number of water molecules in 1mL water = 0.0556 moles (3sf) = 3.35e22 molecules (3sf)

Volume of 10000 water molecules = 2.99e-19 mL (3sf) = cube with side length of ~ 6.69 nm (3sf)

== ==
'''<big>TASK</big>: Consider an atom at position <math>\left(0.5, 0.5, 0.5\right)</math> in a cubic simulation box which runs from <math>\left(0, 0, 0\right)</math> to <math>\left(1, 1, 1\right)</math>. In a single timestep, it moves along the vector <math>\left(0.7, 0.6, 0.2\right)</math>. At what point does it end up, ''after the periodic boundary conditions have been applied''?'''.

After PBC has been applied the atom will be at (0.2,0.1,0.7)

== ==
'''<big>TASK</big>: The Lennard-Jones parameters for argon are <math>\sigma = 0.34\mathrm{nm}, \epsilon\ /\ k_B= 120 \mathrm{K}</math>. If the LJ cutoff is <math>r^* = 3.2</math>, what is it in real units? What is the well depth in <math>\mathrm{kJ\ mol}^{-1}</math>? What is the reduced temperature <math>T^* = 1.5</math> in real units?

<math>r=r^*\times\sigma=3.2\times0.34=1.088nm</math>

<math>T=T^*\frac{\epsilon}{k_B}=1.5\times120=180K</math>

<math>\phi\left(r\right)=-6.16\times10^{-24}Jmol^{-1}</math>

== ==
'''<big>TASK</big>: Why do you think giving atoms random starting coordinates causes problems in simulations? Hint: what happens if two atoms happen to be generated close together?'''

If atoms are created too close, then we will have a very large potential, meaning that if the timestep is not reduced, then we will have a large change in energy between timesteps which could mean that we do not obtain all the information from the simulation system. However if we decrease the timestep so that the change in energy is reduced to a suitable size, then it would require much more calculations to run the simulation for the same of time.

== ==
'''<big>TASK</big>: Satisfy yourself that this lattice spacing corresponds to a number density of lattice points of <math>0.8</math>. Consider instead a face-centred cubic lattice with a lattice point number density of 1.2. What is the side length of the cubic unit cell?'''

Simple cubic: 0.8 density = 1.25 cell volume = <math>\sqrt[3]{1.25}</math> = 1.07722 cell length.

Face centre cubic: 1.2 density = 0.833 cell volume = <math>\sqrt[3]{0.833}</math> = 0.941 cell length.

== ==
'''<big>TASK</big>: Consider again the face-centred cubic lattice from the previous task. How many atoms would be created by the create_atoms command if you had defined that lattice instead?'''

FCC contains 4 atoms in each lattice cell, SC contains 1 atom in each lattice cell, thus 4000 atoms would be created if FCC was used.

== ==
'''<big>TASK</big>: Using the [http://lammps.sandia.gov/doc/Section_commands.html#cmd_5 LAMMPS manual], find the purpose of the following commands in the input script:'''

<pre>
mass 1 1.0
pair_style lj/cut 3.0
pair_coeff * * 1.0 1.0
</pre>

mass 1 1.0: mass of atom type 1 is 1.0

pair_style lj/cut 3.0: cut-off for lennard jones interaction is 3 unit length

pair_coeff * * 1.0 1.0: pair_coeff refers to the force field coeffecients for pairs of atoms, * * means between any two atom types, 1.0 is the coefficient

== ==
'''<big>TASK</big>: Given that we are specifying <math>\mathbf{x}_i\left(0\right)</math> and <math>\mathbf{v}_i\left(0\right)</math>, which integration algorithm are we going to use?'''

Velocity velvet integration

== ==
'''<big>TASK</big>: Look at the lines below.'''
<pre>
### SPECIFY TIMESTEP ###
variable timestep equal 0.001
variable n_steps equal floor(100/${timestep})
variable n_steps equal floor(100/0.001)
timestep ${timestep}
timestep 0.001

### RUN SIMULATION ###
run ${n_steps}
run 100000
</pre>
'''The second line (starting "variable timestep...") tells LAMMPS that if it encounters the text ${timestep} on a subsequent line, it should replace it by the value given. In this case, the value ${timestep} is always replaced by 0.001. In light of this, what do you think the purpose of these lines is? Why not just write:'''
<pre>
timestep 0.001
run 100000
</pre>

The reason why a variable is used to hold the value of timestep is so that it is more convenient to change the value of timestep, as you would only have to change the line that defines the timestep instead of going through the whole code and changing all instances where timestep is used.

== ==
'''<big>TASK</big>: make plots of the energy, temperature, and pressure, against time for the 0.001 timestep experiment (attach a picture to your report). Does the simulation reach equilibrium? How long does this take? When you have done this, make a single plot which shows the energy versus time for all of the timesteps (again, attach a picture to your report). Choosing a timestep is a balancing act: the shorter the timestep, the more accurately the results of your simulation will reflect the physical reality; short timesteps, however, mean that the same number of simulation steps cover a shorter amount of actual time, and this is very unhelpful if the process you want to study requires observation over a long time. Of the five timesteps that you used, which is the largest to give acceptable results? Which one of the five is a ''particularly'' bad choice? Why?'''

[[File:Sh5214 liquid sim intro 001 energy.png]] [[File:Sh5214 liquid sim intro 001 temp.png]] [[File:Sh5214 liquid sim intro 001 pressure.png]] [[File:Sh5214 liquid sim intro all energy.png]] [[File:Sh5214 liquid sim intro all energy zoom.png]]

From the plot of energy vs time for timestep = 0.001 we can see that equilibrium is reached, and from the last plot it can be seen that it takes about 0.4 units of time to reach equilibrium. The largest timestep that results in the equilibration of energy is timestep = 0.01. Timestep = 0.015 is a bad choice because the system does not equilibrate and instead diverges.

== ==
'''<big>TASK</big>: We need to choose <math>\gamma</math> so that the temperature is correct <math>T = \mathfrak{T}</math> if we multiply every velocity <math>\gamma</math>. We can write two equations:'''

<math>\frac{1}{2}\sum_i m_i v_i^2 = \frac{3}{2} N k_B T</math>

<math>\frac{1}{2}\sum_i m_i \left(\gamma v_i\right)^2 = \frac{3}{2} N k_B \mathfrak{T}</math>

'''Solve these to determine <math>\gamma</math>.'''

<math>\frac{1}{2}\sum_i m_i \left(\gamma v_i\right)^2 = \frac{3}{2} N k_B </math>

<math>\implies \frac{1}{2}\sum_i m_i v_i^2=\frac{3Nk_B\mathfrak{T}}{2\gamma^2}=\frac{3}{2}Nk_BT</math>

<math>\implies \gamma^2=\frac{\mathfrak{T}}{T} \implies\gamma=\sqrt{\frac{\mathfrak{T}}{T}}</math>

== ==
'''<big>TASK</big>: Use the [http://lammps.sandia.gov/doc/fix_ave_time.html manual page] to find out the importance of the three numbers ''100 1000 100000''. How often will values of the temperature, etc., be sampled for the average? How many measurements contribute to the average? Looking to the following line, how much time will you simulate?'''

100: number of timesteps between each sample

1000: number of samples for each average

100000: number of timesteps required to obtain one average

== ==
'''<big>TASK</big>: When your simulations have finished, download the log files as before. At the end of the log file, LAMMPS will output the values and errors for the pressure, temperature, and density <math>\left(\frac{N}{V}\right)</math>. Use software of your choice to plot the density as a function of temperature for both of the pressures that you simulated. Your graph(s) should include error bars in both the x and y directions. You should also include a line corresponding to the density predicted by the ideal gas law at that pressure. Is your simulated density lower or higher? Justify this. Does the discrepancy increase or decrease with pressure?'''

[[File:Sh5214 liquid sim npt density temperature.png]]

My simulated density is lower than the density predicted by the idea gas law, this discrepancy increases with pressure. The reason why the deal gas law gives a higher density is due to the fact that it doesn't take into account the interactions between atoms, hence repulsions between atoms are not accounted for in the ideal gas law, this means that the atoms would be closer together in the idea gas law and so a higher density.

== ==
'''<big>TASK</big>: As in the last section, you need to run simulations at ten phase points. In this section, we will be in density-temperature <math>\left(\rho^*, T^*\right)</math> phase space, rather than pressure-temperature phase space. The two densities required at <math>0.2</math> and <math>0.8</math>, and the temperature range is <math>2.0, 2.2, 2.4, 2.6, 2.8</math>. Plot <math>C_V/V</math> as a function of temperature, where <math>V</math> is the volume of the simulation cell, for both of your densities (on the same graph). Is the trend the one you would expect? Attach an example of one of your input scripts to your report.'''

[[File:Sh5214 nvt heat capacity.png]]

From the graph it can be seen that specific heat capacity decreases with temperature, this is expected because there is a limit to how much energy a molecule can hold, as the temperature of the molecule increases, the amount of available rotational/vibrational energy level that the molecule can be excited into decreases, meaning that the molecule will be able to store less energy as its temperature increases.

We can also see that the specific heat capacity for the system with 0.8 density is higher than the system with 0.2 density, this is expected because a higher density means that there are more molecules per unit volume, and therefore more energy could be stored per unit volume.

<pre>
### DEFINE SIMULATION BOX GEOMETRY ###
lattice sc 0.2
region box block 0 15 0 15 0 15
create_box 1 box
create_atoms 1 box

### DEFINE PHYSICAL PROPERTIES OF ATOMS ###
mass 1 1.0
pair_style lj/cut/opt 3.0
pair_coeff 1 1 1.0 1.0
neighbor 2.0 bin

### SPECIFY THE REQUIRED THERMODYNAMIC STATE ###
variable T equal 2.0
variable timestep equal 0.0025

### ASSIGN ATOMIC VELOCITIES ###
velocity all create ${T} 12345 dist gaussian rot yes mom yes

### SPECIFY ENSEMBLE ###
timestep ${timestep}
fix nve all nve

### THERMODYNAMIC OUTPUT CONTROL ###
thermo_style custom time etotal temp press
thermo 10

### RECORD TRAJECTORY ###
dump traj all custom 1000 output-1 id x y z

### SPECIFY TIMESTEP ###

### RUN SIMULATION TO MELT CRYSTAL ###
run 10000
unfix nve
reset_timestep 0

### BRING SYSTEM TO REQUIRED STATE ###
variable tdamp equal ${timestep}*100
variable vdamp equal ${timestep}*1000
fix nvt all nvt temp ${T} ${T} ${tdamp}
run 10000
reset_timestep 0

## SWITCH OFF THERMOSTAT ##
unfix nvt
fix nve all nve

### MEASURE SYSTEM STATE ###
thermo_style custom step etotal temp press density atoms vol
variable vol equal vol
variable N2 equal atoms*atoms
variable dens equal density
variable dens2 equal density*density
variable temp equal temp
variable temp2 equal temp*temp
variable press equal press
variable press2 equal press*press
variable etotal equal etotal
variable etotal2 equal etotal*etotal
fix aves all ave/time 100 1000 100000 v_dens v_temp v_press v_etotal v_dens2 v_temp2 v_press2 v_etotal2 v_vol
run 100000

variable avetemp equal f_aves[2]
variable errtemp equal sqrt(f_aves[6]-f_aves[2]*f_aves[2])
variable cv equal (${N2}*(f_aves[8]-f_aves[4]*f_aves[4])/(${T}*${T}))
variable V equal f_aves[9]

print "Averages"
print "--------"
print "Temperature: ${avetemp}"
print "Stderr: ${errtemp}"
print "cv: ${cv}"
print "volume: ${V}"
</pre>

== ==
'''<big>TASK</big>: perform simulations of the Lennard-Jones system in the three phases. When each is complete, download the trajectory and calculate <math>g(r)</math> and <math>\int g(r)\mathrm{d}r</math>. Plot the RDFs for the three systems on the same axes, and attach a copy to your report. Discuss qualitatively the differences between the three RDFs, and what this tells you about the structure of the system in each phase. In the solid case, illustrate which lattice sites the first three peaks correspond to. What is the lattice spacing? What is the coordination number for each of the first three peaks?'''

[[File:Sh5214 liquid sim rdf plot.png]]

Firstly the range of peaks for solids is longer than liquids which is longer than gas, this means that in the system of solids, the atoms are more ordered than the atoms in the system of liquids, which are more ordered than the atoms in the system of gases. Another observation that could be made is that the magnitude of the peaks decreases as r increases, this is because g(r) is calculated by number of atoms in the shell divided by volume of shell and the density of the system, as r increases, the volume of the shell will increase at a faster rate than the number of atoms in the shell, and thus g(r) decreases. Also it can be seen that magnitude of the peaks of the RDF of the solid is greater than both liquid and gas, this is because solids are more closely packed and more ordered than both gas and liquid, hence at certain r values there will be a high number of atoms and at other r values there will be a very low number of atoms, and that's why solids have higher peaks and lower troughs.

The first peak for the solid RDF is at r=1.025 and the second peak is at r=1.425, thus the lattice spacing is 0.4 times the interatomic separation. There are three different planes in a fcc crystal, they have miller indices of (1,1,1) (1,1,0) and (1,0,0). The (1,1,1) plane corresponds to the first peak and has a coordination number of 6 atoms, the (1,1,0) plane corresponds to the second peak and has a coordination number of 2, the (1,0,0) plane corresponds to the third peak and has a coordination number of 4. The coordination number of lattice sites will affect the magnitude of the peak, and that's why magnitude of first peak > third peak > second peak

== ==
'''<big>TASK</big>: make a plot for each of your simulations (solid, liquid, and gas), showing the mean squared displacement (the "total" MSD) as a function of timestep. Are these as you would expect? Estimate <math>D</math> in each case. Be careful with the units! Repeat this procedure for the MSD data that you were given from the one million atom simulations.'''

[[File:Sh5214 liquid sim d gas msd.png]] [[File:Sh5214 liquid sim d liquid msd.png]] [[File:Sh5214 liquid sim d solid msd.png]]

Gas MSD (my data):
<pre>
t=2000, MSD=29.2
t=5000, MSD=104.7
gradient = (104.7-29.2)/3000 = 0.0252 (3sf)
D = gradient/6 = 0.00419 (3sf)
</pre>

Gas MSD (1e6 atoms):
<pre>
t=2000, MSD=36.4
t=5000, MSD=144.4
gradient = (144.4-36.4)/3000 = 0.036
D = gradient/6 = 0.006 (3sf)
</pre>

Liquid MSD (my data):
<pre>
t=0, MSD=0
t=5000, MSD=11
gradient = 11/5000 = 0.0022
D = gradient/6 = 0.000367 (3sf)
</pre>

Liquid MSD (1e6 atoms):
<pre>
t=0, MSD=0
t=5000, MSD=5.19
gradient = 5.19/5000 = 0.00104 (3sf)
D = gradient/6 = 0.000173
</pre>

Solid MSD (my data): D=0

Solid MSD (1e6 atoms): D=0

MSD is the average distance traveled by a particle in a system, it can be seen that the MSD of the solid plateaus very early on and at a very low value this is due to the fact that atoms in a solid are rigid and fixed in place therefore average distance traveled by a atom in a solid is very small. We see that the MSD plot of the gas starts off as a parabola before turning linear, the parabola part corresponds to when gas atom doesn't collide with anything else meaning that during this time the velocity of the atom is constant and so distance increases proportional to time squared, when more and more collisions occur the curve goes from parabolic to linear. In the liquid MSD plot, we can see that at the start of the curve for a short amount of time the curve is not linear, this is because little to no collisions occur right at the start of the simulation, for most of the plot the line is linear due to the fact that there are collisions between liquid atoms.

The magnitude of MSD for gas > liquid > solid, this is due to the fact that density of gas < liquid < solid, a lower density means that less collisions occur and thus the average translational velocity of atoms is higher, resulting in greater distance traveled.

== ==
'''<big>TASK</big>: In the theoretical section at the beginning, the equation for the evolution of the position of a 1D harmonic oscillator as a function of time was given. Using this, evaluate the normalised velocity autocorrelation function for a 1D harmonic oscillator (it is analytic!):

<math>C\left(\tau\right) = \frac{\int_{-\infty}^{\infty} v\left(t\right)v\left(t + \tau\right)\mathrm{d}t}{\int_{-\infty}^{\infty} v^2\left(t\right)\mathrm{d}t}</math>

'''Be sure to show your working in your writeup. On the same graph, with x range 0 to 500, plot <math>C\left(\tau\right)</math> with <math>\omega = 1/2\pi</math> and the VACFs from your liquid and solid simulations. What do the minima in the VACFs for the liquid and solid system represent? Discuss the origin of the differences between the liquid and solid VACFs. The harmonic oscillator VACF is very different to the Lennard Jones solid and liquid. Why is this? Attach a copy of your plot to your writeup.'''

<math>x(t)=Acos(\omega t+\phi)</math>

<math>v(t)=\frac{dx}{dt}=-A\omega sin(\omega t+\phi)</math>

<math>C(\tau)=\frac{\int_{-\infty}^{\infty}-A\omega sin(\omega t+\phi).-A\omega sin(\omega t+\omega\tau+\phi)\mathrm{d}t}{\int_{-\infty}^{\infty}A^2\omega^2sin^2(\omega t+\phi)\mathrm{d}t}</math>

<math>=\frac{\int_{-\infty}^{\infty}0.5A^2\omega^2[cos(\omega\tau)-cos(2\omega t+\omega\tau+2\phi)]dt}{\int_{-\infty}^{\infty}0.5A^2\omega^2-0.5A^2\omega^2cos(2\omega t+2\phi)}</math>

<math>=\frac{[0.5A^2\omega^2cos(\omega\tau)t+0.25A^2\omega sin(2\omega t+\omega\tau+2\phi)]_{-\infty}^\infty}{[0.5A^2\omega^2t+0.25A^2\omega cos(2\omega t+2\phi)]_{-\infty}^\infty}</math>

We know that sine and cosine has a range of -1 to 1, hense the above equation can be simplified to:

<math>\approx\frac{[0.5A^2\omega^2cos(\omega\tau)t]_{-\infty}^\infty}{[0.5A^2\omega^2t]_{-\infty}^\infty}</math>

<math>\approx cos(\omega\tau)</math>

[[File:Sh5214 liquide sim d solid liquid cos vacf normalised.png]]

The minima in the VACFs for the liquid and solid system represents the time at which the velocity of the atoms in the system is at a maximum. We can see that the VACF for solids oscillates about zero, with the amplitude of the oscillation will decrease over time, this is due to the fact that the atoms in the solids are locked up in a lattice, these atoms in the lattice will want to be at the optimal position where the repulsive forces and attractive forces balance out, they will therefore oscillate about this position which will also result in an oscillation in velocity, the reason why the magnitude of these oscillation decreases over time is because of other forces that perturb the oscillation. We can see that there is also a very slight minima on the liquid VACF curve, which shows that liquids also oscillate but the oscillation disappears very quickly, the reason why atoms in liquids don't oscillate for as long as solids is because liquid atoms are not locked in a lattice like solids and so are free to move, this means that any oscillation is very quickly dampened out by diffusion.

The reason why the harmonic oscillator VACF is so different is due to the fact that it doesn't take into account dampening of the oscillation, hence the magnitude of the oscillation never decreases.

==?==
'''<big>TASK</big>: Use the trapezium rule to approximate the integral under the velocity autocorrelation function for the solid, liquid, and gas, and use these values to estimate <math>D</math> in each case. You should make a plot of the running integral in each case. Are they as you expect? Repeat this procedure for the VACF data that you were given from the one million atom simulations. What do you think is the largest source of error in your estimates of D from the VACF?'''

[[File:Sh5214 liquid sim d gas integral vacf.png]] [[File:Sh5214 liquid sim d liquid integral vacf.png]] [[File:Sh5214 liquid sim d solid integral vacf.png]]

Gas VACF (my data):
<pre>
t=5000, Running Integral=6.54 (3sf)
D = Running Integral/3 = 2.18 (3sf)
</pre>

Gas VACF (1e6 atoms):
<pre>
t=5000, Running Integral=9.81 (3sf)
D = Running Integral/3 = 3.27 (3sf)
</pre>

Liquid VACF (my data):
<pre>
t=5000, Running Integral=0.569 (3sf)
D = Running Integral/3 = 0.190 (3sf)
</pre>

Liquid VACF (1e6 atoms):
<pre>
t=5000, Running Integral=0.270 (3sf)
D = Running Integral/3 = 0.0901 (3sf)
</pre>

Solid VACF (my data): D = 0

Solid VACF (1e6 atoms): D = 0