ChemWiki - User contributions [en]

Rep:SMliqsim

2018-01-31T11:30:13Z

Sm6415: /* Methods */

= Liquid Simulations Computational Lab Assignment - Sam Macer =

== Abstract ==
Molecular dynamics simulations of particles interacting according to a Lennard-Jones potential were performed using the velocity-Verlet algorithm. Temperature dependence of density in the NpT ensemble was monitored, and significant deviation from the ideal gas law was observed at low temperatures due to the no collisions assumption breaking down. Heat capacity was monitored as a function of temperature in the NVT ensemble. The radial distribution functions were obtained for solid, liquid and gas systems, as were diffusion coefficients via the mean squared displacement. Furthermore, the diffusion coefficients obtained were realistic; D (solid) = 1.0 x10<sup>5 </sup>m<sup>2</sup>/s), D (liquid = 3.47 x10<sup>9 </sup>m<sup>2</sup>/s), D (gas = 3.47 x10<sup>10 </sup>m<sup>2</sup>/s). Potential industrial applications for such molecular dynamics simulations are also discussed.

== Introduction ==
Diffusion is a fundamental form of transport in fluid and soft-condensed phases. All un-equilibrated fluid systems will move via diffusion toward thermodynamic equilibrium. Many industrial chemical processes could benefit from an optimisation considering diffusion involving the fluid or soft-condensed system being processed. For example, consider piping different coloured icings onto a cupcake to form a pattern. The manufacturer may wish to
prevent the diffusion of the two colouring molecules between the bulk of the icing, as this would ruin the pattern. Studying the temperature dependence of
this diffusion could allow them to find the maximum temperature at which the cupcake
could be safely packaged without significant diffusion occurring.

One study on the extraction of basil essential oil from the
<em>Ocimum basilicum</em> leaf exemplifies this<sup>1</sup>. The diffusion coefficient was
found to change over the course of the extraction; the vapour pressure of the
basil oil in the surrounding air increases, while the homogeneity of oil
distribution in the leaf diminished as it is depleted from different areas at
different rates, due to different tissue structures inside the leaf. The study used models based on Fick’s second law of diffusion, and ultimately elucidated a
quadratic relationship between the concentration and diffusion coefficient for
the system. This knowledge allows manufacturers to fit the model
mathematically using parameters specific to their distillation setup, and use the result to find the best trade-off between extraction time (time and energy are expensive) and extraction yield.

== Aims & Objectives ==
The aim of this experiment is to use computational models of particles interacting according to a Lennard-Jones potential to extract thermodynamic data about the system and ultimately calculate the diffusion coefficient. The experiment uses boxes of between 1000 and
10000 atoms, with periodic boundary conditions applied to simulate a much larger 'effective system'. The velocity-Verlet
algorithm<sup>2</sup> is used, which iterates the velocity and displacement for every particle in the system, from
an initial state, through the desired number of steps. At each of these steps, thermodynamic
state variables of the system can be extracted, for e.g. Temperature can be obtained using the average kinetic energies in conjunction with the equipartition theorem.

Reduced units throughout the experiment to make the magnitudes of the quantities used more managable.
Given the Lennard-Jones potential:

<div align="center">
[[File:Smliq1.png|188x188px]]

r = r*/σ

E = E*/ε

T = kT*/ε
</div>

''r = inter-particle distance, σ = inter-particle distance when the potential = 0, E = potential, ε = depth of the potential well, T = temperature and k = the Boltzmann constant. Starred quantities represent the non-reduced
versions. Note the unites cancel out when forming these quantities and hence they are unitless.

By cancelling through these reduced units, we obtain that time is also in reduced units and can be converted via:

[[File:Smliq122.png]]

== Methods ==
The open source software, LAMMPS, was used to run the molecular dynamics simulations. The experiment was run as 5 separate parts:
# The effect of the time step on the simulation was probed using a box of 1000 atoms, equally spaced with a number density of 0.8 atoms per unit volume. The ensemble was treated as microcanonical. The following time steps were trialed: 0.0150, 0.0100, 0.0075, 0.0025, 0.0010. Energy, temperature and pressure for the system were collected at each iteration. A timestep of 0.0025 was used for the two subsequent simulations - see results and discussion for justification.
# A set of simulations (1000 atoms) were run in the isobaric-isothermal (NpT) ensemble. At P = 2.5 & P = 3.0, the simulation was run over the following temperatures: T = 1.5, 2.0, 2.5, 3.5, 5. Density of the systems was extracted as an average after the simulation, and plotted against temperature to test the agreement with the ideal gas law.
# An experiment (1000 atoms) was run in the NVT ensemble at densities 0.2 & 0.8, over the temperatures T = 2.0, 2.2, 2.4, 2.6, 2.8. The average heat capacity and temperature were collected for each system.
# Simulations (8000 atoms) were performed for a solid, a liquid and a gas system. A timestep of 0.0020 was used. Parameters for this were selected based on a previous probe into the phase transitions of the Lennard-Jones system<sup>3</sup>; solid (density = 1.2, T = 0.5), liquid (density = 0.8, T = 1.2), gas (density = 0.2, T = 1.2). The output trajectories were used to calculate the RDF (radial distribution function), i.e. the radial frequency of meeting another particle averaged over all spherical angles.
# Simulations (8000 atoms) were again performed for a solid, a liquid and a gas phase system. A timestep of 0.0020 was used for this experiment, and the same density/temperature parameters at experiment 4. The MSD (mean squared displacement) was collected as a function of timestep.

== Results & Discussion ==

=== Experiment 1 ===
<div align="center">
[[File:Smliq2.png|404x404px]]

[[File:Smliq3.png|406x406px]]

[[File:Smliq4.png|468x468px]]

These plots show how the thermodynamic state variables of the system vary with time. For all of the variables , we see a very quick convergence from the initial value to the steady state value. A closer inspection of these graphs reveal
equilibration happens after roughly 0.5 reduced time units. This shows how the system which is initially ordered quickly descends into chaos. After this point, the state variables proceed to oscillate randomly as a noise function. This occurs due to fluctuations in the system which is constantly moving. If the system had more particles, the average would be less noisy because there are more opportunities for the random imbalances to cancel each other out.

[[File:Smliq5.png|514x514px]]

The above plot shows how the energy varies with timestep. As the time step is increased the average energy
increases due to the resolution of the numerical solution decreasing. The largest timestep shows a clear divergence. Presumably the timestep is so large that particles are getting unrealistically close to each other. The reaction in the next timestep cascades such that the problem is transferred to neighbouring particles. At this point the simulation breaks down as the energies and particle positions become increasingly unrealistic. 0.0025 was the largest timestep which showed no sign of this increased energy or divergence, and hence the timestep used for all subsequent experiments was smaller than this.

</div>

=== Experiment 2 ===
<div align="center">

[[File:Smliq6.png|501x501px]]

The density decreases with temperature as expected; the
atoms have more kinetic energy, so push each other apart further during their
random motion. The higher pressure plots show a consistently higher density.
This is intuitive as at higher pressures, the atoms are squashed closer
together on average.

At lower temperatures, the ideal gas law result diverges
from the simulation. This is because the ideal gas law assumes no collisions.
At lower temperatures, particles of an ideal gas continue to be compressed into
a smaller space due to their diminishing kinetic energy, such that they are in
fact overlapping. In the simulation the particles collide, meaning a much
greater kinetic energy (temperature) reduction is required to achieve the same
increase in density. Also note that the ideal gas law allows for infinity as the
limiting density, whereas this becomes infinitely unlikely when using Van Der
Waals potential between the atoms.

</div>

=== Experiment 3 ===
<div align="center">
[[File:Smliq7.png|486x486px]]

The heat capacity is smaller for lower densities
as expected: there are less atoms to sink energy into, so the average kinetic energy
increases more quickly. The heat capacity decreases with temperature because at higher temperatures, the energy states are closer together, meaning it is easier to populate new states (higher density of states). This clearly outweighs any increased degeneracy in states, otherwise we would see the heat capacity increase with temperature as is often the case. The plateau could indicate a first order phase transition where the heat capacity tends to infinity. This would be made more clear with a higher resolution of data points. The radial distribution function could be calculated either side as the less ordered phase would have a smoother curve in the case of a transition.

</div>

=== Experiment 4 ===
<div align="center">
[[File:Smliq8.png|478x478px]]

The radial distribution function (RDF) represents the density of particles around a single particle in the system on average as a function of radius, averaged over all possible angles in a sphere. The RDFs for the different phases are distinct. In the solid case, we se sharp peaks corresponding to when we hit a site in the crystalline lattice. In the liquid, the peaks are much less sharp, as there are many more positions the molecules can occupy easily relative to one another. Peaks are still defined as the distance between particles is roughly constant. A gas behaves like a liquid with very little order. It is less peaks if any (occur at higher densities) are much less well defined. There is still a sharp peak close to the particle as neighbouring particles spend the most time when changing direction during a collision as this is when velocity is the smallest. All three phases exhibit a cut-off displacement due to the steep positive region in the Lennard-Jones potential.

[[File:Smliq9.png|533x533px]]

When the RDF is viewed cumulatively, we see the steepest cumulation for the solid, followed by the liquid then gas. This is due to meeting particles more frequently when walking along the radius in the more dense phase.

</div>

=== Experiment 5 ===
<div align="center">

[[File:Smliq10.png|687x687px]]

The mean squared displacement plot used the same colours as before to represent a solid, liquid and gas. Initially the MSD is slow to rise, in particular for the gas. This is due to the particles taking some time to disperse properly and reach equilibrium velocity.

The gas has the most freedom and kinetic energy so the MSD is the greatest. The opposite is true for the solid. The liquid is closer to the solid in terms of MSD because it is still relatively ordered, and there is some kind of activation barrier to molecules swapping between cages. The gradient for a solid is effectively zero as swapping positions in a crystalline lattice is extremely rare.

The second plot shows a method of least squared fit, taking into account points after timestep = 1000 only, to avoid data from the equilibration period being considered. The gradients were as following:

Solid: 2.94 x10<sup>-8</sup>

Liquid: 1.02 x10<sup>-3</sup>

Gas: 1.00 x10<sup>-2</sup>

The MSD is related to the diffusion coefficient by a simple relationship:

[[File:Smliq11.png]]

By simply using this relation, we obtain D as following (in timestep units):

Solid: 4.9 x10<sup>-9</sup>

Liquid: 1.7 x10<sup>-4</sup>

Gas: 1.7 x10<sup>-3</sup>

To convert to the correct units, we must note 1 timestep corresponds to 48.89 fs, therefore we must divide by 49 x10<sup>-15</sup> to get D in SI units:

Solid: 1.0 x10<sup>5 </sup>m<sup>2</sup>/s

Liquid: 3.47 x10<sup>9 </sup>m<sup>2</sup>/s

Gas: 3.47 x10<sup>10 </sup>m<sup>2</sup>/s

Which is reasonable as diffusion coefficients for liquids are normally ~ 1 x10<sup>9 </sup>m<sup>2</sup>/s

Now it has been show we can model
the diffusion coefficient based purely on molecular dynamics simulations, we can see how similar simulations could be applied to real world
applications. In the basil leaf case introduced in the introduction, we could model the leaf
structure to follow the real leaf tissue, by using different densities. We would then be able to model accurately how the diffusion of basil
essential oil molecules changes as a vapour shell forms and molecules in the
leaf are depleted.
</div>

== Conclusion ==
The primary conclusion of this experiment is that molecular dynamics simulations are suitable for physical modelling, and can give more accurate results than ideal solutions such as the ideal gas law. The next step is to begin applying these models to real systems and test their ability to model accurately, such as in the basil leaf case.

The molecular dynamics simulations gave a better estimate of the thermodynamic state variables than the ideal gas law because collisions were included. This is particularly relevant at low temperatures and high pressures. The heat capacity simulations showed that it is possible to observe phase transitions between phases. The RDF is a unique insight that would be much more difficult to measure in reality. It is a useful tool as is gives us detailed information on the structure of the system being studied. In a solid for example, the RDF peaks can be correlated to crystalline lattice sites. Having access to the displacements and velocities of each particle allows us to calculate a wealth of parameters, the diffusion coefficient being just one example. Obtaining a realistic diffusion coefficient implies that the method models the system realistically.

It is clear that the amount of important data that can be quickly obtained for a system simulated using molecular dynamics is immense. It would be recommended that such methods are tested against real world experiments to gain a better understanding of their accuracy.

== References ==
# J. Silveira, A. Costa and E. Costa Junior, ''Engenharia Agrícola'', 2017, 37, 717-726.
# L. Verlet, ''Physical Review'', 1967, 159, 98-103.
# J. Hansen and L. Verlet, ''Physical Review'', 1969, 184, 151-161.

Weblinks:
# <nowiki>http://www.scielo.br/scielo.php?script=sci_arttext&pid=S0100-69162017000400717#B21</nowiki>
# <nowiki>https://journals.aps.org/pr/abstract/10.1103/PhysRev.159.98</nowiki>
# <nowiki>https://journals.aps.org/pr/abstract/10.1103/PhysRev.184.151</nowiki>

== Appendix ==

=== Previously unanswered tasks:[[File:Smliqt1.png|974x974px]] ===

[[File:Smliqt2.png|953x953px]]

[[File:Smliqt3.png|1021x1021px]]

[[File:Smliqt4.png]]

[[File:Smliqt5.png|958x958px]]

[[File:Smliqt6.png|979x979px]]

=== LAMMPS script for heat capacity experiments: ===

<nowiki>###</nowiki> DEFINE SIMULATION BOX GEOMETRY ###

variable dens equal 0.2

lattice sc ${dens}

region box block 0 15 0 15 0 15

create_box 1 box

create_atoms 1 box

<nowiki>###</nowiki> 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

<nowiki>###</nowiki> SPECIFY THE REQUIRED THERMODYNAMIC STATE ###

variable T equal 2.0

variable timestep equal 0.0075

<nowiki>###</nowiki> ASSIGN ATOMIC VELOCITIES ###

velocity all create ${T} 12345 dist gaussian rot yes mom yes

<nowiki>###</nowiki> SPECIFY ENSEMBLE ###

timestep ${timestep}

fix nve all nve

<nowiki>###</nowiki> THERMODYNAMIC OUTPUT CONTROL ###

thermo_style custom time etotal temp press

thermo 10

<nowiki>###</nowiki> RECORD TRAJECTORY ###

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

<nowiki>###</nowiki> SPECIFY TIMESTEP ###

<nowiki>###</nowiki> RUN SIMULATION TO MELT CRYSTAL ###

run 10000

unfix nve

reset_timestep 0

<nowiki>###</nowiki> BRING SYSTEM TO REQUIRED STATE ###

variable tdamp equal ${timestep}*100

fix nvt all nvt temp ${T} ${T} ${tdamp}

run 10000

reset_timestep 0

unfix nvt

fix nve all nve

<nowiki>###</nowiki> MEASURE SYSTEM STATE ###

thermo_style custom etotal temp atoms vol

variable vol equal vol

variable atoms equal atoms

variable temp equal temp

variable temp2 equal temp*temp

variable N2 equal atoms*atoms

variable E2 equal etotal*etotal

variable E equal etotal

fix aves all ave/time 100 1000 100000 v_temp2 v_E v_E2

run 100000

variable heatcapac equal
${N2}*(f_aves[3]-(f_aves[2]*f_aves[2]))/(f_aves[1])

print "Averages"

print "--------"

print "Heat Capacity: ${heatcapac}"

print "Volume: ${vol}"

print "No of atoms: ${atoms}"

Rep:SMliqsim

2018-01-31T11:21:42Z

Sm6415:

= Liquid Simulations Computational Lab Assignment - Sam Macer =

== Abstract ==
Molecular dynamics simulations of particles interacting according to a Lennard-Jones potential were performed using the velocity-Verlet algorithm. Temperature dependence of density in the NpT ensemble was monitored, and significant deviation from the ideal gas law was observed at low temperatures due to the no collisions assumption breaking down. Heat capacity was monitored as a function of temperature in the NVT ensemble. The radial distribution functions were obtained for solid, liquid and gas systems, as were diffusion coefficients via the mean squared displacement. Furthermore, the diffusion coefficients obtained were realistic; D (solid) = 1.0 x10<sup>5 </sup>m<sup>2</sup>/s), D (liquid = 3.47 x10<sup>9 </sup>m<sup>2</sup>/s), D (gas = 3.47 x10<sup>10 </sup>m<sup>2</sup>/s). Potential industrial applications for such molecular dynamics simulations are also discussed.

== Introduction ==
Diffusion is a fundamental form of transport in fluid and soft-condensed phases. All un-equilibrated fluid systems will move via diffusion toward thermodynamic equilibrium. Many industrial chemical processes could benefit from an optimisation considering diffusion involving the fluid or soft-condensed system being processed. For example, consider piping different coloured icings onto a cupcake to form a pattern. The manufacturer may wish to
prevent the diffusion of the two colouring molecules between the bulk of the icing, as this would ruin the pattern. Studying the temperature dependence of
this diffusion could allow them to find the maximum temperature at which the cupcake
could be safely packaged without significant diffusion occurring.

One study on the extraction of basil essential oil from the
<em>Ocimum basilicum</em> leaf exemplifies this<sup>1</sup>. The diffusion coefficient was
found to change over the course of the extraction; the vapour pressure of the
basil oil in the surrounding air increases, while the homogeneity of oil
distribution in the leaf diminished as it is depleted from different areas at
different rates, due to different tissue structures inside the leaf. The study used models based on Fick’s second law of diffusion, and ultimately elucidated a
quadratic relationship between the concentration and diffusion coefficient for
the system. This knowledge allows manufacturers to fit the model
mathematically using parameters specific to their distillation setup, and use the result to find the best trade-off between extraction time (time and energy are expensive) and extraction yield.

== Aims & Objectives ==
The aim of this experiment is to use computational models of particles interacting according to a Lennard-Jones potential to extract thermodynamic data about the system and ultimately calculate the diffusion coefficient. The experiment uses boxes of between 1000 and
10000 atoms, with periodic boundary conditions applied to simulate a much larger 'effective system'. The velocity-Verlet
algorithm<sup>2</sup> is used, which iterates the velocity and displacement for every particle in the system, from
an initial state, through the desired number of steps. At each of these steps, thermodynamic
state variables of the system can be extracted, for e.g. Temperature can be obtained using the average kinetic energies in conjunction with the equipartition theorem.

Reduced units throughout the experiment to make the magnitudes of the quantities used more managable.
Given the Lennard-Jones potential:

<div align="center">
[[File:Smliq1.png|188x188px]]

r = r*/σ

E = E*/ε

T = kT*/ε
</div>

''r = inter-particle distance, σ = inter-particle distance when the potential = 0, E = potential, ε = depth of the potential well, T = temperature and k = the Boltzmann constant. Starred quantities represent the non-reduced
versions. Note the unites cancel out when forming these quantities and hence they are unitless.

By cancelling through these reduced units, we obtain that time is also in reduced units and can be converted via:

[[File:Smliq122.png]]

== Methods ==
The open source software, LAMMPS, was used to run the molecular dynamics simulations. The experiment was run as 5 separate parts:
# The effect of the time step on the simulation was probed using a box of 1000 atoms, equally spaced with a number density of 0.8 atoms per unit volume. The ensemble was treated as microcanonical. The following time steps were trialed: 0.0150, 0.0100, 0.0075, 0.0025, 0.0010. Energy, temperature and pressure for the system were collected at each iteration. A timestep of 0.0025 was used for the two subsequent simulations - see results and discussion for justification.
# A set of simulations were run in the isobaric-isothermal (NpT) ensemble. At P = 2.5 & P = 3.0, the simulation was run over the following temperatures: T = 1.5, 2.0, 2.5, 3.5, 5. Density of the systems was extracted as an average after the simulation, and plotted against temperature to test the agreement with the ideal gas law.
# An experiment was run in the NVT ensemble at densities 0.2 & 0.8, over the temperatures T = 2.0, 2.2, 2.4, 2.6, 2.8. The average heat capacity and temperature were collected for each system.
# Simulations were performed for a solid, a liquid and a gas system. A timestep of 0.0020 was used. Parameters for this were selected based on a previous probe into the phase transitions of the Lennard-Jones system<sup>3</sup>; solid (density = 1.2, T = 0.5), liquid (density = 0.8, T = 1.2), gas (density = 0.2, T = 1.2). The output trajectories were used to calculate the RDF (radial distribution function), i.e. the radial frequency of meeting another particle averaged over all spherical angles.
# Simulations were again performed for a solid, a liquid and a gas phase system. A timestep of 0.0020 was used for this experiment, and the same density/temperature parameters at experiment 4. The MSD (mean squared displacement) was collected as a function of timestep.

== Results & Discussion ==

=== Experiment 1 ===
<div align="center">
[[File:Smliq2.png|404x404px]]

[[File:Smliq3.png|406x406px]]

[[File:Smliq4.png|468x468px]]

These plots show how the thermodynamic state variables of the system vary with time. For all of the variables , we see a very quick convergence from the initial value to the steady state value. A closer inspection of these graphs reveal
equilibration happens after roughly 0.5 reduced time units. This shows how the system which is initially ordered quickly descends into chaos. After this point, the state variables proceed to oscillate randomly as a noise function. This occurs due to fluctuations in the system which is constantly moving. If the system had more particles, the average would be less noisy because there are more opportunities for the random imbalances to cancel each other out.

[[File:Smliq5.png|514x514px]]

The above plot shows how the energy varies with timestep. As the time step is increased the average energy
increases due to the resolution of the numerical solution decreasing. The largest timestep shows a clear divergence. Presumably the timestep is so large that particles are getting unrealistically close to each other. The reaction in the next timestep cascades such that the problem is transferred to neighbouring particles. At this point the simulation breaks down as the energies and particle positions become increasingly unrealistic. 0.0025 was the largest timestep which showed no sign of this increased energy or divergence, and hence the timestep used for all subsequent experiments was smaller than this.

</div>

=== Experiment 2 ===
<div align="center">

[[File:Smliq6.png|501x501px]]

The density decreases with temperature as expected; the
atoms have more kinetic energy, so push each other apart further during their
random motion. The higher pressure plots show a consistently higher density.
This is intuitive as at higher pressures, the atoms are squashed closer
together on average.

At lower temperatures, the ideal gas law result diverges
from the simulation. This is because the ideal gas law assumes no collisions.
At lower temperatures, particles of an ideal gas continue to be compressed into
a smaller space due to their diminishing kinetic energy, such that they are in
fact overlapping. In the simulation the particles collide, meaning a much
greater kinetic energy (temperature) reduction is required to achieve the same
increase in density. Also note that the ideal gas law allows for infinity as the
limiting density, whereas this becomes infinitely unlikely when using Van Der
Waals potential between the atoms.

</div>

=== Experiment 3 ===
<div align="center">
[[File:Smliq7.png|486x486px]]

The heat capacity is smaller for lower densities
as expected: there are less atoms to sink energy into, so the average kinetic energy
increases more quickly. The heat capacity decreases with temperature because at higher temperatures, the energy states are closer together, meaning it is easier to populate new states (higher density of states). This clearly outweighs any increased degeneracy in states, otherwise we would see the heat capacity increase with temperature as is often the case. The plateau could indicate a first order phase transition where the heat capacity tends to infinity. This would be made more clear with a higher resolution of data points. The radial distribution function could be calculated either side as the less ordered phase would have a smoother curve.

</div>

=== Experiment 4 ===
<div align="center">
[[File:Smliq8.png|478x478px]]

The radial distribution function (RDF) represents the density of particles around a single particle in the system on average as a function of radius, averaged over all possible angles in a sphere. The RDFs for the different phases are distinct. In the solid case, we se sharp peaks corresponding to when we hit a site in the crystalline lattice. In the liquid, the peaks are much less sharp, as there are many more positions the molecules can occupy easily relative to one another. Peaks are still defined as the distance between particles is roughly constant. A gas behaves like a liquid with very little order. It is less peaks if any (occur at higher densities) are much less well defined. There is still a sharp peak close to the particle as neighbouring particles spend the most time when changing direction during a collision as this is when velocity is the smallest. All three phases exhibit a cut-off displacement due to the steep positive region in the Lennard-Jones potential.

[[File:Smliq9.png|533x533px]]

When the RDF is viewed cumulatively, we see the steepest cumulation for the solid, followed by the liquid then gas. This is due to meeting particles more frequently when in the more dense phase.

</div>

=== Experiment 5 ===
<div align="center">

[[File:Smliq10.png|687x687px]]

The mean squared displacement plot used the same colours as before to represent a solid, liquid and gas. Initially the MSD is slow to rise, in particular for the gas. This is due to the particles taking some time to disperse properly and reach equilibrium velocity. The gas has the most freedom and kinetic energy so the MSD is the greatest. The opposite is true for the solid. The liquid is closer to the solid in terms of MSD because it is still relatively ordered, and there is some kind of activation barrier to molecules swapping between cages. The gradient for a solid is effectively zero as swapping positions in a crystalline lattice is extremely rare.

The second plot shows a method of least squared fit, taking into account points after timestep = 1000 only, to avoid data from the equilibration period being considered. The gradients were as following:

Solid: 2.94 x10<sup>-8</sup>

Liquid: 1.02 x10<sup>-3</sup>

Gas: 1.00 x10<sup>-2</sup>

The MSD is related to the diffusion coefficient by a simple relationship:

[[File:Smliq11.png]]

By simply using this relation, we obtain D as following (in timestep units):

Solid: 4.9 x10<sup>-9</sup>

Liquid: 1.7 x10<sup>-4</sup>

Gas: 1.7 x10<sup>-3</sup>

To convert to the correct units, we must note 1 timestep corresponds to 48.89 fs, therefore we must divide by 49 x10<sup>-15</sup> to get D in SI units:

Solid: 1.0 x10<sup>5 </sup>m<sup>2</sup>/s

Liquid: 3.47 x10<sup>9 </sup>m<sup>2</sup>/s

Gas: 3.47 x10<sup>10 </sup>m<sup>2</sup>/s

Which is reasonable as diffusion coefficients for liquids are normally ~ 1 x10<sup>9 </sup>m<sup>2</sup>/s

Now it has been show we can model
the diffusion coefficient based purely on molecular dynamics simulations, we can see how similar simulations could be applied to real world
applications. In the basil leaf case introduced in the introduction, we could model the leaf
structure to follow the real leaf tissue, by using different densities. We would then be able to model accurately how the diffusion of basil
essential oil molecules changes as a vapour shell forms and molecules in the
leaf are depleted.
</div>

== Conclusion ==
The primary conclusion of this experiment is that molecular dynamics simulations are suitable for physical modelling, and can give more accurate results than ideal solutions such as the ideal gas law. The next step is to begin applying these models to real systems and test their ability to model accurately.

The molecular dynamics simulations gave a better estimate of the thermodynamic state variables because collisions were included. This is particularly relevant at low temperatures and high pressures. The heat capacity simulations showed that it is possible to observe phase transitions between phases. The RDF is a unique insight that would be much more difficult to measure in reality. It is a useful tool as is gives us detailed information on the structure of the system being studied. In a solid for example, the RDF peaks can be correlated to crystalline lattice sites. Having access to the displacements and velocities of each particle allows us to calculate a wealth of parameters, the diffusion coefficient being just one example. Obtaining a realistic diffusion coefficient implies that the method models the system realistically.

It is clear that the amount of important data that can be quickly obtained for a system simulated using molecular dynamics is immense. It would be recommended that such methods are tested against real world experiments to gain a better understanding of their accuracy.

== References ==
# J. Silveira, A. Costa and E. Costa Junior, ''Engenharia Agrícola'', 2017, 37, 717-726.
# L. Verlet, ''Physical Review'', 1967, 159, 98-103.
# J. Hansen and L. Verlet, ''Physical Review'', 1969, 184, 151-161.

Weblinks:
# <nowiki>http://www.scielo.br/scielo.php?script=sci_arttext&pid=S0100-69162017000400717#B21</nowiki>
# <nowiki>https://journals.aps.org/pr/abstract/10.1103/PhysRev.159.98</nowiki>
# <nowiki>https://journals.aps.org/pr/abstract/10.1103/PhysRev.184.151</nowiki>

== Appendix ==

=== Previously unanswered tasks:[[File:Smliqt1.png|974x974px]] ===

[[File:Smliqt2.png|953x953px]]

[[File:Smliqt3.png|1021x1021px]]

[[File:Smliqt4.png]]

[[File:Smliqt5.png|958x958px]]

[[File:Smliqt6.png|979x979px]]

=== LAMMPS script for heat capacity experiments: ===

<nowiki>###</nowiki> DEFINE SIMULATION BOX GEOMETRY ###

variable dens equal 0.2

lattice sc ${dens}

region box block 0 15 0 15 0 15

create_box 1 box

create_atoms 1 box

<nowiki>###</nowiki> 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

<nowiki>###</nowiki> SPECIFY THE REQUIRED THERMODYNAMIC STATE ###

variable T equal 2.0

variable timestep equal 0.0075

<nowiki>###</nowiki> ASSIGN ATOMIC VELOCITIES ###

velocity all create ${T} 12345 dist gaussian rot yes mom yes

<nowiki>###</nowiki> SPECIFY ENSEMBLE ###

timestep ${timestep}

fix nve all nve

<nowiki>###</nowiki> THERMODYNAMIC OUTPUT CONTROL ###

thermo_style custom time etotal temp press

thermo 10

<nowiki>###</nowiki> RECORD TRAJECTORY ###

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

<nowiki>###</nowiki> SPECIFY TIMESTEP ###

<nowiki>###</nowiki> RUN SIMULATION TO MELT CRYSTAL ###

run 10000

unfix nve

reset_timestep 0

<nowiki>###</nowiki> BRING SYSTEM TO REQUIRED STATE ###

variable tdamp equal ${timestep}*100

fix nvt all nvt temp ${T} ${T} ${tdamp}

run 10000

reset_timestep 0

unfix nvt

fix nve all nve

<nowiki>###</nowiki> MEASURE SYSTEM STATE ###

thermo_style custom etotal temp atoms vol

variable vol equal vol

variable atoms equal atoms

variable temp equal temp

variable temp2 equal temp*temp

variable N2 equal atoms*atoms

variable E2 equal etotal*etotal

variable E equal etotal

fix aves all ave/time 100 1000 100000 v_temp2 v_E v_E2

run 100000

variable heatcapac equal
${N2}*(f_aves[3]-(f_aves[2]*f_aves[2]))/(f_aves[1])

print "Averages"

print "--------"

print "Heat Capacity: ${heatcapac}"

print "Volume: ${vol}"

print "No of atoms: ${atoms}"

File:Smliqt4.png

2018-01-31T10:49:45Z

Sm6415:

Rep:SMliqsim

2018-01-31T10:47:24Z

Sm6415: /* Appendix */

= Liquid Simulations Computational Lab Assignment - Sam Macer =

== Abstract ==
Molecular dynamics simulations of particles interacting according to a Lennard-Jones potential were performed using the velocity-Verlet algorithm. Temperature dependence of density in the NpT ensemble was monitored, and significant deviation from the ideal gas law was observed at low temperatures due to the no collisions assumption breaking down. Heat capacity was monitored as a function of temperature in the NVT ensemble. The radial distribution functions were obtained for solid, liquid and gas systems, as were diffusion coefficients via the mean squared displacement. Potential industrial applications for such molecular dynamics simulations are also discussed.

== Introduction ==
Diffusion is a fundamental form of transport in fluid and soft-condensed phases. All un-equilibrated fluid systems will move via diffusion toward thermodynamic equilibrium. Many industrial chemical processes could benefit from an optimisation considering diffusion involving the fluid or soft-condensed system being processed. For example, consider piping different coloured icings onto a cupcake to form a pattern. The manufacturer may wish to
prevent the diffusion of the two colouring molecules between the bulk of the icing, as this would ruin the pattern. Studying the temperature dependence of
this diffusion could allow them to find the maximum temperature at which the cupcake
could be safely packaged without significant diffusion occurring.

One study on the extraction of basil essential oil from the
<em>Ocimum basilicum</em> leaf exemplifies this<sup>1</sup>. The diffusion coefficient was
found to change over the course of the extraction; the vapour pressure of the
basil oil in the surrounding air increases, while the homogeneity of oil
distribution in the leaf diminished as it is depleted from different areas at
different rates, due to different tissue structures inside the leaf. The study used models based on Fick’s second law of diffusion, and ultimately elucidated a
quadratic relationship between the concentration and diffusion coefficient for
the system. This knowledge allows manufacturers to fit the model
mathematically using parameters specific to their distillation setup, and use the result to find the best trade-off between extraction time (time and energy are expensive) and extraction yield.

== Aims & Objectives ==
The aim of this experiment is to use computational models of particles interacting according to a Lennard-Jones potential to extract thermodynamic data about the system and ultimately calculate the diffusion coefficient. The experiment uses boxes of between 1000 and
10000 atoms, with periodic boundary conditions applied to simulate a much larger 'effective system'. The velocity-Verlet
algorithm<sup>2</sup> is used, which iterates the velocity and displacement for every particle in the system, from
an initial state, through the desired number of steps. At each of these steps, thermodynamic
state variables of the system can be extracted, for e.g. Temperature can be obtained using the average kinetic energies in conjunction with the equipartition theorem.

Reduced units throughout the experiment to make the magnitudes of the quantities used more managable.
Given the Lennard-Jones potential:

<div align="center">
[[File:Smliq1.png|188x188px]]

r = r*/σ

E = E*/ε

T = kT*/ε
</div>

''r = inter-particle distance, σ = inter-particle distance when the potential = 0, E = potential, ε = depth of the potential well, T = temperature and k = the Boltzmann constant. Starred quantities represent the non-reduced
versions. Note the unites cancel out when forming these quantities and hence they are unitless.

By cancelling through these reduced units, we obtain that time is also in reduced units and can be converted via:

[[File:Smliq122.png]]

This is used for calculating the diffusion coefficient, though may be used to obtain the true timescale for any of the simulations.

== Methods ==
The open source software, LAMMPS, was used to run the molecular dynamics simulations. The experiment was run as 5 separate parts:
# The effect of the time step on the simulation was probed using a box of 1000 atoms, equally spaced with a number density of 0.8 atoms per unit volume. The ensemble was treated as microcanonical. The following time steps were trialed: 0.0150, 0.0100, 0.0075, 0.0025, 0.0010. Energy, temperature and pressure for the system were collected at each iteration. A timestep of 0.0025 was used for the two subsequent simulations - see results and discussion for justification.
# A set of simulations were run in the isobaric-isothermal (NpT) ensemble. At P = 2.5 & P = 3.0, the simulation was run over the following temperatures: T = 1.5, 2.0, 2.5, 3.5, 5. Density of the systems was extracted as an average after the simulation, and plotted against temperature to test the agreement with the ideal gas law.
# An experiment was run in the NVT ensemble at densities 0.2 & 0.8, over the temperatures T = 2.0, 2.2, 2.4, 2.6, 2.8. The average heat capacity and temperature were collected for each system.
# Simulations were performed for a solid, a liquid and a gas system. A timestep of 0.0020 was used. Parameters for this were selected based on a previous probe into the phase transitions of the Lennard-Jones system<sup>3</sup>; solid (density = 1.2, T = 0.5), liquid (density = 0.8, T = 1.2), gas (density = 0.2, T = 1.2). The output trajectories were used to calculate the RDF (radial distribution function), i.e. the radial frequency of meeting another particle averaged over all spherical angles.
# Simulations were again performed for a solid, a liquid and a gas phase system. A timestep of 0.0020 was used for this experiment, and the same density/temperature parameters at experiment 4. The MSD (mean squared displacement) was collected as a function of timestep.

== Results & Discussion ==

=== Experiment 1 ===
<div align="center">
[[File:Smliq2.png|404x404px]]

[[File:Smliq3.png|406x406px]]

[[File:Smliq4.png|468x468px]]

These plots show how the thermodynamic state variables of the system vary with time. For all of the variables , we see a very quick convergence from the initial value to the steady state value. A closer inspection of these graphs reveal
equilibration happens after roughly 0.5 reduced time units. This shows how the system which is initially ordered quickly descends into chaos. After this point, the state variables proceed to oscillate randomly as a noise function. This occurs due to fluctuations in the system which is constantly moving. If the system had more particles, the average would be less noisy because there are more opportunities for the random imbalances to cancel each other out.

[[File:Smliq5.png|514x514px]]

The above plot shows how the energy varies with timestep. As the time step is increased the average energy
increases due to the resolution of the numerical solution decreasing. The largest timestep shows a clear divergence. Presumably the timestep is so large that particles are getting unrealistically close to each other. The reaction in the next timestep cascades such that the problem is transferred to neighbouring particles. At this point the simulation breaks down as the energies and particle positions become increasingly unrealistic. 0.0025 was the largest timestep which showed no sign of this increased energy or divergence, and hence the timestep used for all subsequent experiments was smaller than this.

</div>

=== Experiment 2 ===
<div align="center">

[[File:Smliq6.png|501x501px]]

The density decreases with temperature as expected; the
atoms have more kinetic energy, so push each other apart further during their
random motion. The higher pressure plots show a consistently higher density.
This is intuitive as at higher pressures, the atoms are squashed closer
together on average.

At lower temperatures, the ideal gas law result diverges
from the simulation. This is because the ideal gas law assumes no collisions.
At lower temperatures, particles of an ideal gas continue to be compressed into
a smaller space due to their diminishing kinetic energy, such that they are in
fact overlapping. In the simulation the particles collide, meaning a much
greater kinetic energy (temperature) reduction is required to achieve the same
increase in density. Also note that the ideal gas law allows for infinity as the
limiting density, whereas this becomes infinitely unlikely when using Van Der
Waals potential between the atoms.

</div>

=== Experiment 3 ===
<div align="center">
[[File:Smliq7.png|486x486px]]

The heat capacity is smaller for lower densities
as expected: there are less atoms to sink energy into, so the average kinetic energy
increases more quickly. The heat capacity decreases with temperature because at higher temperatures, the energy states are closer together, meaning it is easier to populate new states (higher density of states). This clearly outweighs any increased degeneracy in states, otherwise we would see the heat capacity increase with temperature as is often the case. The plateau could indicate a first order phase transition where the heat capacity tends to infinity. This would be made more clear with a higher resolution of data points. The radial distribution function could be calculated either side as the less ordered phase would have a smoother curve.

</div>

=== Experiment 4 ===
<div align="center">
[[File:Smliq8.png|478x478px]]

The radial distribution function (RDF) represents the density of particles around a single particle in the system on average as a function of radius, averaged over all possible angles in a sphere. The RDFs for the different phases are distinct. In the solid case, we se sharp peaks corresponding to when we hit a site in the crystalline lattice. In the liquid, the peaks are much less sharp, as there are many more positions the molecules can occupy easily relative to one another. Peaks are still defined as the distance between particles is roughly constant. A gas behaves like a liquid with very little order. It is less peaks if any (occur at higher densities) are much less well defined. There is still a sharp peak close to the particle as neighbouring particles spend the most time when changing direction during a collision as this is when velocity is the smallest. All three phases exhibit a cut-off displacement due to the steep positive region in the Lennard-Jones potential.

[[File:Smliq9.png|533x533px]]

When the RDF is viewed cumulatively, we see the steepest cumulation for the solid, followed by the liquid then gas. This is due to meeting particles more frequently when in the more dense phase.

</div>

=== Experiment 5 ===
<div align="center">

[[File:Smliq10.png|687x687px]]

The mean squared displacement plot used the same colours as before to represent a solid, liquid and gas. Initially the MSD is slow to rise, in particular for the gas. This is due to the particles taking some time to disperse properly and reach equilibrium velocity. The gas has the most freedom and kinetic energy so the MSD is the greatest. The opposite is true for the solid. The liquid is closer to the solid in terms of MSD because it is still relatively ordered, and there is some kind of activation barrier to molecules swapping between cages. The gradient for a solid is effectively zero as swapping positions in a crystalline lattice is extremely rare.

The second plot shows a method of least squared fit, taking into account points after timestep = 1000 only, to avoid data from the equilibration period being considered. The gradients were as following:

Solid: 2.94 x10<sup>-8</sup>

Liquid: 1.02 x10<sup>-3</sup>

Gas: 1.00 x10<sup>-2</sup>

The MSD is related to the diffusion coefficient by a simple relationship:

[[File:Smliq11.png]]

By simply using this relation, we obtain D as following (in timestep units):

Solid: 4.9 x10<sup>-9</sup>

Liquid: 1.7 x10<sup>-4</sup>

Gas: 1.7 x10<sup>-3</sup>

To convert to the correct units, we must simply apply:

[[File:Smliq122.png]]

Now it has been show we can model
the diffusion coefficient based purely on molecular dynamics simulations, we can see how similar simulations could be applied to real world
applications. In the basil leaf case introduced in the introduction, we could model the leaf
structure to follow the real leaf tissue, by using different densities. We would then be able to model accurately how the diffusion of basil
essential oil molecules changes as a vapour shell forms and molecules in the
leaf are depleted.
</div>

== Conclusion ==
The primary conclusion of this experiment is that molecular dynamics simulations are suitable for physical modelling, and can give more accurate results than ideal solutions such as the ideal gas law. The next step is to begin applying these models to real systems and test their ability to model accurately.

The molecular dynamics simulations gave a better estimate of the thermodynamic state variables because collisions were included. This is particularly relevant at low temperatures and high pressures. The heat capacity simulations showed that it is possible to observe phase transitions between phases. The RDF is a unique insight that would be much more difficult to measure in reality. It is a useful tool as is gives us detailed information on the structure of the system being studied. In a solid for example, the RDF peaks can be correlated to crystalline lattice sites. Having access to the displacements and velocities of each particle allows us to calculate a wealth of parameters, the diffusion coefficient being just one example.

It is clear that the amount of important data that can be quickly obtained for a system simulated using molecular dynamics is immense. It would be recommended that such methods are tested against real world experiments to gain a better understanding of their accuracy.

== References ==
# J. Silveira, A. Costa and E. Costa Junior, ''Engenharia Agrícola'', 2017, 37, 717-726.
# L. Verlet, ''Physical Review'', 1967, 159, 98-103.
# J. Hansen and L. Verlet, ''Physical Review'', 1969, 184, 151-161.

Weblinks:
# <nowiki>http://www.scielo.br/scielo.php?script=sci_arttext&pid=S0100-69162017000400717#B21</nowiki>
# <nowiki>https://journals.aps.org/pr/abstract/10.1103/PhysRev.159.98</nowiki>
# <nowiki>https://journals.aps.org/pr/abstract/10.1103/PhysRev.184.151</nowiki>

== Appendix ==

=== Previously unanswered tasks:[[File:Smliqt1.png|974x974px]] ===

[[File:Smliqt2.png|953x953px]]

[[File:Smliqt3.png|1021x1021px]]

[[File:Smliqt4.png]]

[[File:Smliqt5.png|958x958px]]

[[File:Smliqt6.png|979x979px]]

=== LAMMPS script for heat capacity experiments: ===

<nowiki>###</nowiki> DEFINE SIMULATION BOX GEOMETRY ###

variable dens equal 0.2

lattice sc ${dens}

region box block 0 15 0 15 0 15

create_box 1 box

create_atoms 1 box

<nowiki>###</nowiki> 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

<nowiki>###</nowiki> SPECIFY THE REQUIRED THERMODYNAMIC STATE ###

variable T equal 2.0

variable timestep equal 0.0075

<nowiki>###</nowiki> ASSIGN ATOMIC VELOCITIES ###

velocity all create ${T} 12345 dist gaussian rot yes mom yes

<nowiki>###</nowiki> SPECIFY ENSEMBLE ###

timestep ${timestep}

fix nve all nve

<nowiki>###</nowiki> THERMODYNAMIC OUTPUT CONTROL ###

thermo_style custom time etotal temp press

thermo 10

<nowiki>###</nowiki> RECORD TRAJECTORY ###

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

<nowiki>###</nowiki> SPECIFY TIMESTEP ###

<nowiki>###</nowiki> RUN SIMULATION TO MELT CRYSTAL ###

run 10000

unfix nve

reset_timestep 0

<nowiki>###</nowiki> BRING SYSTEM TO REQUIRED STATE ###

variable tdamp equal ${timestep}*100

fix nvt all nvt temp ${T} ${T} ${tdamp}

run 10000

reset_timestep 0

unfix nvt

fix nve all nve

<nowiki>###</nowiki> MEASURE SYSTEM STATE ###

thermo_style custom etotal temp atoms vol

variable vol equal vol

variable atoms equal atoms

variable temp equal temp

variable temp2 equal temp*temp

variable N2 equal atoms*atoms

variable E2 equal etotal*etotal

variable E equal etotal

fix aves all ave/time 100 1000 100000 v_temp2 v_E v_E2

run 100000

variable heatcapac equal
${N2}*(f_aves[3]-(f_aves[2]*f_aves[2]))/(f_aves[1])

print "Averages"

print "--------"

print "Heat Capacity: ${heatcapac}"

print "Volume: ${vol}"

print "No of atoms: ${atoms}"

Rep:SMliqsim

2018-01-31T10:43:46Z

Sm6415:

= Liquid Simulations Computational Lab Assignment - Sam Macer =

== Abstract ==
Molecular dynamics simulations of particles interacting according to a Lennard-Jones potential were performed using the velocity-Verlet algorithm. Temperature dependence of density in the NpT ensemble was monitored, and significant deviation from the ideal gas law was observed at low temperatures due to the no collisions assumption breaking down. Heat capacity was monitored as a function of temperature in the NVT ensemble. The radial distribution functions were obtained for solid, liquid and gas systems, as were diffusion coefficients via the mean squared displacement. Potential industrial applications for such molecular dynamics simulations are also discussed.

== Introduction ==
Diffusion is a fundamental form of transport in fluid and soft-condensed phases. All un-equilibrated fluid systems will move via diffusion toward thermodynamic equilibrium. Many industrial chemical processes could benefit from an optimisation considering diffusion involving the fluid or soft-condensed system being processed. For example, consider piping different coloured icings onto a cupcake to form a pattern. The manufacturer may wish to
prevent the diffusion of the two colouring molecules between the bulk of the icing, as this would ruin the pattern. Studying the temperature dependence of
this diffusion could allow them to find the maximum temperature at which the cupcake
could be safely packaged without significant diffusion occurring.

One study on the extraction of basil essential oil from the
<em>Ocimum basilicum</em> leaf exemplifies this<sup>1</sup>. The diffusion coefficient was
found to change over the course of the extraction; the vapour pressure of the
basil oil in the surrounding air increases, while the homogeneity of oil
distribution in the leaf diminished as it is depleted from different areas at
different rates, due to different tissue structures inside the leaf. The study used models based on Fick’s second law of diffusion, and ultimately elucidated a
quadratic relationship between the concentration and diffusion coefficient for
the system. This knowledge allows manufacturers to fit the model
mathematically using parameters specific to their distillation setup, and use the result to find the best trade-off between extraction time (time and energy are expensive) and extraction yield.

== Aims & Objectives ==
The aim of this experiment is to use computational models of particles interacting according to a Lennard-Jones potential to extract thermodynamic data about the system and ultimately calculate the diffusion coefficient. The experiment uses boxes of between 1000 and
10000 atoms, with periodic boundary conditions applied to simulate a much larger 'effective system'. The velocity-Verlet
algorithm<sup>2</sup> is used, which iterates the velocity and displacement for every particle in the system, from
an initial state, through the desired number of steps. At each of these steps, thermodynamic
state variables of the system can be extracted, for e.g. Temperature can be obtained using the average kinetic energies in conjunction with the equipartition theorem.

Reduced units throughout the experiment to make the magnitudes of the quantities used more managable.
Given the Lennard-Jones potential:

<div align="center">
[[File:Smliq1.png|188x188px]]

r = r*/σ

E = E*/ε

T = kT*/ε
</div>

''r = inter-particle distance, σ = inter-particle distance when the potential = 0, E = potential, ε = depth of the potential well, T = temperature and k = the Boltzmann constant. Starred quantities represent the non-reduced
versions. Note the unites cancel out when forming these quantities and hence they are unitless.

By cancelling through these reduced units, we obtain that time is also in reduced units and can be converted via:

[[File:Smliq122.png]]

This is used for calculating the diffusion coefficient, though may be used to obtain the true timescale for any of the simulations.

== Methods ==
The open source software, LAMMPS, was used to run the molecular dynamics simulations. The experiment was run as 5 separate parts:
# The effect of the time step on the simulation was probed using a box of 1000 atoms, equally spaced with a number density of 0.8 atoms per unit volume. The ensemble was treated as microcanonical. The following time steps were trialed: 0.0150, 0.0100, 0.0075, 0.0025, 0.0010. Energy, temperature and pressure for the system were collected at each iteration. A timestep of 0.0025 was used for the two subsequent simulations - see results and discussion for justification.
# A set of simulations were run in the isobaric-isothermal (NpT) ensemble. At P = 2.5 & P = 3.0, the simulation was run over the following temperatures: T = 1.5, 2.0, 2.5, 3.5, 5. Density of the systems was extracted as an average after the simulation, and plotted against temperature to test the agreement with the ideal gas law.
# An experiment was run in the NVT ensemble at densities 0.2 & 0.8, over the temperatures T = 2.0, 2.2, 2.4, 2.6, 2.8. The average heat capacity and temperature were collected for each system.
# Simulations were performed for a solid, a liquid and a gas system. A timestep of 0.0020 was used. Parameters for this were selected based on a previous probe into the phase transitions of the Lennard-Jones system<sup>3</sup>; solid (density = 1.2, T = 0.5), liquid (density = 0.8, T = 1.2), gas (density = 0.2, T = 1.2). The output trajectories were used to calculate the RDF (radial distribution function), i.e. the radial frequency of meeting another particle averaged over all spherical angles.
# Simulations were again performed for a solid, a liquid and a gas phase system. A timestep of 0.0020 was used for this experiment, and the same density/temperature parameters at experiment 4. The MSD (mean squared displacement) was collected as a function of timestep.

== Results & Discussion ==

=== Experiment 1 ===
<div align="center">
[[File:Smliq2.png|404x404px]]

[[File:Smliq3.png|406x406px]]

[[File:Smliq4.png|468x468px]]

These plots show how the thermodynamic state variables of the system vary with time. For all of the variables , we see a very quick convergence from the initial value to the steady state value. A closer inspection of these graphs reveal
equilibration happens after roughly 0.5 reduced time units. This shows how the system which is initially ordered quickly descends into chaos. After this point, the state variables proceed to oscillate randomly as a noise function. This occurs due to fluctuations in the system which is constantly moving. If the system had more particles, the average would be less noisy because there are more opportunities for the random imbalances to cancel each other out.

[[File:Smliq5.png|514x514px]]

The above plot shows how the energy varies with timestep. As the time step is increased the average energy
increases due to the resolution of the numerical solution decreasing. The largest timestep shows a clear divergence. Presumably the timestep is so large that particles are getting unrealistically close to each other. The reaction in the next timestep cascades such that the problem is transferred to neighbouring particles. At this point the simulation breaks down as the energies and particle positions become increasingly unrealistic. 0.0025 was the largest timestep which showed no sign of this increased energy or divergence, and hence the timestep used for all subsequent experiments was smaller than this.

</div>

=== Experiment 2 ===
<div align="center">

[[File:Smliq6.png|501x501px]]

The density decreases with temperature as expected; the
atoms have more kinetic energy, so push each other apart further during their
random motion. The higher pressure plots show a consistently higher density.
This is intuitive as at higher pressures, the atoms are squashed closer
together on average.

At lower temperatures, the ideal gas law result diverges
from the simulation. This is because the ideal gas law assumes no collisions.
At lower temperatures, particles of an ideal gas continue to be compressed into
a smaller space due to their diminishing kinetic energy, such that they are in
fact overlapping. In the simulation the particles collide, meaning a much
greater kinetic energy (temperature) reduction is required to achieve the same
increase in density. Also note that the ideal gas law allows for infinity as the
limiting density, whereas this becomes infinitely unlikely when using Van Der
Waals potential between the atoms.

</div>

=== Experiment 3 ===
<div align="center">
[[File:Smliq7.png|486x486px]]

The heat capacity is smaller for lower densities
as expected: there are less atoms to sink energy into, so the average kinetic energy
increases more quickly. The heat capacity decreases with temperature because at higher temperatures, the energy states are closer together, meaning it is easier to populate new states (higher density of states). This clearly outweighs any increased degeneracy in states, otherwise we would see the heat capacity increase with temperature as is often the case. The plateau could indicate a first order phase transition where the heat capacity tends to infinity. This would be made more clear with a higher resolution of data points. The radial distribution function could be calculated either side as the less ordered phase would have a smoother curve.

</div>

=== Experiment 4 ===
<div align="center">
[[File:Smliq8.png|478x478px]]

The radial distribution function (RDF) represents the density of particles around a single particle in the system on average as a function of radius, averaged over all possible angles in a sphere. The RDFs for the different phases are distinct. In the solid case, we se sharp peaks corresponding to when we hit a site in the crystalline lattice. In the liquid, the peaks are much less sharp, as there are many more positions the molecules can occupy easily relative to one another. Peaks are still defined as the distance between particles is roughly constant. A gas behaves like a liquid with very little order. It is less peaks if any (occur at higher densities) are much less well defined. There is still a sharp peak close to the particle as neighbouring particles spend the most time when changing direction during a collision as this is when velocity is the smallest. All three phases exhibit a cut-off displacement due to the steep positive region in the Lennard-Jones potential.

[[File:Smliq9.png|533x533px]]

When the RDF is viewed cumulatively, we see the steepest cumulation for the solid, followed by the liquid then gas. This is due to meeting particles more frequently when in the more dense phase.

</div>

=== Experiment 5 ===
<div align="center">

[[File:Smliq10.png|687x687px]]

The mean squared displacement plot used the same colours as before to represent a solid, liquid and gas. Initially the MSD is slow to rise, in particular for the gas. This is due to the particles taking some time to disperse properly and reach equilibrium velocity. The gas has the most freedom and kinetic energy so the MSD is the greatest. The opposite is true for the solid. The liquid is closer to the solid in terms of MSD because it is still relatively ordered, and there is some kind of activation barrier to molecules swapping between cages. The gradient for a solid is effectively zero as swapping positions in a crystalline lattice is extremely rare.

The second plot shows a method of least squared fit, taking into account points after timestep = 1000 only, to avoid data from the equilibration period being considered. The gradients were as following:

Solid: 2.94 x10<sup>-8</sup>

Liquid: 1.02 x10<sup>-3</sup>

Gas: 1.00 x10<sup>-2</sup>

The MSD is related to the diffusion coefficient by a simple relationship:

[[File:Smliq11.png]]

By simply using this relation, we obtain D as following (in timestep units):

Solid: 4.9 x10<sup>-9</sup>

Liquid: 1.7 x10<sup>-4</sup>

Gas: 1.7 x10<sup>-3</sup>

To convert to the correct units, we must simply apply:

[[File:Smliq122.png]]

Now it has been show we can model
the diffusion coefficient based purely on molecular dynamics simulations, we can see how similar simulations could be applied to real world
applications. In the basil leaf case introduced in the introduction, we could model the leaf
structure to follow the real leaf tissue, by using different densities. We would then be able to model accurately how the diffusion of basil
essential oil molecules changes as a vapour shell forms and molecules in the
leaf are depleted.
</div>

== Conclusion ==
The primary conclusion of this experiment is that molecular dynamics simulations are suitable for physical modelling, and can give more accurate results than ideal solutions such as the ideal gas law. The next step is to begin applying these models to real systems and test their ability to model accurately.

The molecular dynamics simulations gave a better estimate of the thermodynamic state variables because collisions were included. This is particularly relevant at low temperatures and high pressures. The heat capacity simulations showed that it is possible to observe phase transitions between phases. The RDF is a unique insight that would be much more difficult to measure in reality. It is a useful tool as is gives us detailed information on the structure of the system being studied. In a solid for example, the RDF peaks can be correlated to crystalline lattice sites. Having access to the displacements and velocities of each particle allows us to calculate a wealth of parameters, the diffusion coefficient being just one example.

It is clear that the amount of important data that can be quickly obtained for a system simulated using molecular dynamics is immense. It would be recommended that such methods are tested against real world experiments to gain a better understanding of their accuracy.

== References ==
# J. Silveira, A. Costa and E. Costa Junior, ''Engenharia Agrícola'', 2017, 37, 717-726.
# L. Verlet, ''Physical Review'', 1967, 159, 98-103.
# J. Hansen and L. Verlet, ''Physical Review'', 1969, 184, 151-161.

Weblinks:
# <nowiki>http://www.scielo.br/scielo.php?script=sci_arttext&pid=S0100-69162017000400717#B21</nowiki>
# <nowiki>https://journals.aps.org/pr/abstract/10.1103/PhysRev.159.98</nowiki>
# <nowiki>https://journals.aps.org/pr/abstract/10.1103/PhysRev.184.151</nowiki>

== Appendix ==
Previously unanswered tasks:
[[File:Smliqt1.png]]

[[File:Smliqt2.png]]

[[File:Smliqt3.png]]

[[File:Smliqt4.png]]

[[File:Smliqt5.png]]

[[File:Smliqt6.png]]

LAMMPS script for heat capacity experiments:

<nowiki>###</nowiki> DEFINE SIMULATION BOX GEOMETRY ###

variable dens equal 0.2

lattice sc ${dens}

region box block 0 15 0 15 0 15

create_box 1 box

create_atoms 1 box

<nowiki>###</nowiki> 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

<nowiki>###</nowiki> SPECIFY THE REQUIRED THERMODYNAMIC STATE ###

variable T equal 2.0

variable timestep equal 0.0075

<nowiki>###</nowiki> ASSIGN ATOMIC VELOCITIES ###

velocity all create ${T} 12345 dist gaussian rot yes mom yes

<nowiki>###</nowiki> SPECIFY ENSEMBLE ###

timestep ${timestep}

fix nve all nve

<nowiki>###</nowiki> THERMODYNAMIC OUTPUT CONTROL ###

thermo_style custom time etotal temp press

thermo 10

<nowiki>###</nowiki> RECORD TRAJECTORY ###

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

<nowiki>###</nowiki> SPECIFY TIMESTEP ###

<nowiki>###</nowiki> RUN SIMULATION TO MELT CRYSTAL ###

run 10000

unfix nve

reset_timestep 0

<nowiki>###</nowiki> BRING SYSTEM TO REQUIRED STATE ###

variable tdamp equal ${timestep}*100

fix nvt all nvt temp ${T} ${T} ${tdamp}

run 10000

reset_timestep 0

unfix nvt

fix nve all nve

<nowiki>###</nowiki> MEASURE SYSTEM STATE ###

thermo_style custom etotal temp atoms vol

variable vol equal vol

variable atoms equal atoms

variable temp equal temp

variable temp2 equal temp*temp

variable N2 equal atoms*atoms

variable E2 equal etotal*etotal

variable E equal etotal

fix aves all ave/time 100 1000 100000 v_temp2 v_E v_E2

run 100000

variable heatcapac equal
${N2}*(f_aves[3]-(f_aves[2]*f_aves[2]))/(f_aves[1])

print "Averages"

print "--------"

print "Heat Capacity: ${heatcapac}"

print "Volume: ${vol}"

print "No of atoms: ${atoms}"

File:Smliqt6.png

2018-01-31T10:42:31Z

Sm6415:

Rep:SMliqsim

2018-01-31T10:42:24Z

2018-01-31T09:32:23Z

Sm6415: Sm6415 uploaded a new version of File:Smliq12.png

Rep:SMliqsim

2018-01-31T09:06:57Z

Sm6415:

= Liquid Simulations Computational Lab Assignment - Sam Macer =

== Abstract ==

== Introduction ==
Diffusion is a fundamental form of transport in fluid and soft-condensed phases. All un-equilibrated fluid systems will move via diffusion toward thermodynamic equilibrium. Many industrial chemical processes could benefit from an optimisation considering diffusion involving the fluid or soft-condensed system being processed. For example, consider piping different coloured icings onto a cupcake to form a pattern. The manufacturer may wish to
prevent the diffusion of the two colouring molecules between the bulk of the icing, as this would ruin the pattern. Studying the temperature dependence of
this diffusion could allow them to find the maximum temperature at which the cupcake
could be safely packaged without significant diffusion occurring.

One study on the extraction of basil essential oil from the
<em>Ocimum basilicum</em> leaf exemplifies this<sup>1</sup>. The diffusion coefficient was
found to change over the course of the extraction; the vapour pressure of the
basil oil in the surrounding air increases, while the homogeneity of oil
distribution in the leaf diminished as it is depleted from different areas at
different rates, due to different tissue structures inside the leaf. The study used models based on Fick’s second law of diffusion, and ultimately elucidated a
quadratic relationship between the concentration and diffusion coefficient for
the system. This knowledge allows manufacturers to fit the model
mathematically using parameters specific to their distillation setup, and use the result to find the best trade-off between extraction time (time and energy are expensive) and extraction yield.

== Aims & Objectives ==
The aim of this experiment is to use computational models of particles interacting according to a Lennard-Jones potential to extract thermodynamic data about the system and ultimately calculate the diffusion coefficient. The experiment uses boxes of between 1000 and
10000 atoms, with periodic boundary conditions applied to simulate a much larger 'effective system'. The velocity-Verlet
algorithm<sup>2</sup> is used, which iterates the velocity and displacement for every particle in the system, from
an initial state, through the desired number of steps. At each of these steps, thermodynamic
state variables of the system can be extracted, for e.g. Temperature can be obtained using the average kinetic energies in conjunction with the equipartition theorem.

Reduced units throughout the experiment to make the magnitudes of the quantities used more managable.
Given the Lennard-Jones potential:

<div align="center">
[[File:Smliq1.png|188x188px]]

r = r*/σ

E = E*/ε

T = kT*/ε
</div>

''r = inter-particle distance, σ = inter-particle distance when the potential = 0, E = potential, ε = depth of the potential well, T = temperature and k = the Boltzmann constant. Starred quantities represent the non-reduced
versions. Note the unites cancel out when forming these quantities and hence they are unitless.

By cancelling through these reduced units, we obtain that time is also in reduced units and can be converted via:
[[File:Smliq12.png]]
This is used for calculating the diffusion coefficient, though may be used to obtain the true timescale for any of the simulations.

== Methods ==
The open source software, LAMMPS, was used to run the molecular dynamics simulations. The experiment was run as 5 separate parts:
# The effect of the time step on the simulation was probed using a box of 1000 atoms, equally spaced with a number density of 0.8 atoms per unit volume. The ensemble was treated as microcanonical. The following time steps were trialed: 0.0150, 0.0100, 0.0075, 0.0025, 0.0010. Energy, temperature and pressure for the system were collected at each iteration. A timestep of 0.0025 was used for the two subsequent simulations - see results and discussion for justification.
# A set of simulations were run in the isobaric-isothermal (NpT) ensemble. At P = 2.5 & P = 3.0, the simulation was run over the following temperatures: T = 1.5, 2.0, 2.5, 3.5, 5. Density of the systems was extracted as an average after the simulation, and plotted against temperature to test the agreement with the ideal gas law.
# An experiment was run in the NVT ensemble at densities 0.2 & 0.8, over the temperatures T = 2.0, 2.2, 2.4, 2.6, 2.8. The average heat capacity and temperature were collected for each system.
# Simulations were performed for a solid, a liquid and a gas system. A timestep of 0.0020 was used. Parameters for this were selected based on a previous probe into the phase transitions of the Lennard-Jones system<sup>3</sup>; solid (density = 1.2, T = 0.5), liquid (density = 0.8, T = 1.2), gas (density = 0.2, T = 1.2). The output trajectories were used to calculate the RDF (radial distribution function), i.e. the radial frequency of meeting another particle averaged over all spherical angles.
# Simulations were again performed for a solid, a liquid and a gas phase system. A timestep of 0.0020 was used for this experiment, and the same density/temperature parameters at experiment 4. The MSD (mean squared displacement) was collected as a function of timestep.

== Results & Discussion ==

=== Experiment 1 ===
<div align="center">
[[File:Smliq2.png|404x404px]]

[[File:Smliq3.png|406x406px]]

[[File:Smliq4.png|468x468px]]

These plots show how the thermodynamic state variables of the system vary with time. For all of the variables , we see a very quick convergence from the initial value to the steady state value. A closer inspection of these graphs reveal
equilibration happens after roughly 0.5 reduced time units. This shows how the system which is initially ordered quickly descends into chaos. After this point, the state variables proceed to oscillate randomly as a noise function. This occurs due to fluctuations in the system which is constantly moving. If the system had more particles, the average would be less noisy because there are more opportunities for the random imbalances to cancel each other out.

[[File:Smliq5.png|514x514px]]

The above plot shows how the energy varies with timestep. As the time step is increased the average energy
increases due to the resolution of the numerical solution decreasing. The largest timestep shows a clear divergence. Presumably the timestep is so large that particles are getting unrealistically close to each other. The reaction in the next timestep cascades such that the problem is transferred to neighbouring particles. At this point the simulation breaks down as the energies and particle positions become increasingly unrealistic. 0.0025 was the largest timestep which showed no sign of this increased energy or divergence, and hence the timestep used for all subsequent experiments was smaller than this.

</div>

=== Experiment 2 ===
<div align="center">

[[File:Smliq6.png|501x501px]]

The density decreases with temperature as expected; the
atoms have more kinetic energy, so push each other apart further during their
random motion. The higher pressure plots show a consistently higher density.
This is intuitive as at higher pressures, the atoms are squashed closer
together on average.

At lower temperatures, the ideal gas law result diverges
from the simulation. This is because the ideal gas law assumes no collisions.
At lower temperatures, particles of an ideal gas continue to be compressed into
a smaller space due to their diminishing kinetic energy, such that they are in
fact overlapping. In the simulation the particles collide, meaning a much
greater kinetic energy (temperature) reduction is required to achieve the same
increase in density. Also note that the ideal gas law allows for infinity as the
limiting density, whereas this becomes infinitely unlikely when using Van Der
Waals potential between the atoms.

</div>

=== Experiment 3 ===
<div align="center">
[[File:Smliq7.png|486x486px]]

The heat capacity is smaller for lower densities
as expected: there are less atoms to sink energy into, so the average kinetic energy
increases more quickly. The heat capacity decreases with temperature because at higher temperatures, the energy states are closer together, meaning it is easier to populate new states (higher density of states). This clearly outweighs any increased degeneracy in states, otherwise we would see the heat capacity increase with temperature as is often the case. The plateau could indicate a first order phase transition where the heat capacity tends to infinity. This would be made more clear with a higher resolution of data points. The radial distribution function could be calculated either side as the less ordered phase would have a smoother curve.

</div>

=== Experiment 4 ===
<div align="center">
[[File:Smliq8.png|478x478px]]

The radial distribution function (RDF) represents the density of particles around a single particle in the system on average as a function of radius, averaged over all possible angles in a sphere. The RDFs for the different phases are distinct. In the solid case, we se sharp peaks corresponding to when we hit a site in the crystalline lattice. In the liquid, the peaks are much less sharp, as there are many more positions the molecules can occupy easily relative to one another. Peaks are still defined as the distance between particles is roughly constant. A gas behaves like a liquid with very little order. It is less peaks if any (occur at higher densities) are much less well defined. There is still a sharp peak close to the particle as neighbouring particles spend the most time when changing direction during a collision as this is when velocity is the smallest. All three phases exhibit a cut-off displacement due to the steep positive region in the Lennard-Jones potential.

[[File:Smliq9.png|533x533px]]

When the RDF is viewed cumulatively, we see the steepest cumulation for the solid, followed by the liquid then gas. This is due to meeting particles more frequently when in the more dense phase.

</div>

=== Experiment 5 ===
<div align="center">

[[File:Smliq10.png|687x687px]]

The mean squared displacement plot used the same colours as before to represent a solid, liquid and gas. Initially the MSD is slow to rise, in particular for the gas. This is due to the particles taking some time to disperse properly and reach equilibrium velocity. The gas has the most freedom and kinetic energy so the MSD is the greatest. The opposite is true for the solid. The liquid is closer to the solid in terms of MSD because it is still relatively ordered, and there is some kind of activation barrier to molecules swapping between cages. The gradient for a solid is effectively zero as swapping positions in a crystalline lattice is extremely rare.

The second plot shows a method of least squared fit, taking into account points after timestep = 1000 only, to avoid data from the equilibration period being considered. The gradients were as following:

Solid: 2.94 x10<sup>-8</sup>

Liquid: 1.02 x10<sup>-3</sup>

Gas: 1.00 x10<sup>-2</sup>

The MSD is related to the diffusion coefficient by a simple relationship:

[[File:Smliq11.png]]

Now it has been show we can model
the diffusion coefficient, similar simulations could be applied to real world
applications. For e.g. in the basil leaf example, we could model the leaf
structure with different densities, and see how the diffusion of basil
essential oil molecules changes as a vapour shell forms and molecules in the
leaf are depleted.
</div>

== Conclusion ==

== References ==
# J. Silveira, A. Costa and E. Costa Junior, ''Engenharia Agrícola'', 2017, 37, 717-726.
# L. Verlet, ''Physical Review'', 1967, 159, 98-103.
# J. Hansen and L. Verlet, ''Physical Review'', 1969, 184, 151-161.

Weblinks:
# <nowiki>http://www.scielo.br/scielo.php?script=sci_arttext&pid=S0100-69162017000400717#B21</nowiki>
# <nowiki>https://journals.aps.org/pr/abstract/10.1103/PhysRev.159.98</nowiki>
# <nowiki>https://journals.aps.org/pr/abstract/10.1103/PhysRev.184.151</nowiki>

== Appendix ==

LAMMPS script for heat capacity experiments

Heat capacity LAMMPS input:

<nowiki>###</nowiki> DEFINE SIMULATION BOX GEOMETRY ###

variable dens equal 0.2

lattice sc ${dens}

region box block 0 15 0 15 0 15

create_box 1 box

create_atoms 1 box

<nowiki>###</nowiki> 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

<nowiki>###</nowiki> SPECIFY THE REQUIRED THERMODYNAMIC STATE ###

variable T equal 2.0

variable timestep equal 0.0075

<nowiki>###</nowiki> ASSIGN ATOMIC VELOCITIES ###

velocity all create ${T} 12345 dist gaussian rot yes mom yes

<nowiki>###</nowiki> SPECIFY ENSEMBLE ###

timestep ${timestep}

fix nve all nve

<nowiki>###</nowiki> THERMODYNAMIC OUTPUT CONTROL ###

thermo_style custom time etotal temp press

thermo 10

<nowiki>###</nowiki> RECORD TRAJECTORY ###

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

<nowiki>###</nowiki> SPECIFY TIMESTEP ###

<nowiki>###</nowiki> RUN SIMULATION TO MELT CRYSTAL ###

run 10000

unfix nve

reset_timestep 0

<nowiki>###</nowiki> BRING SYSTEM TO REQUIRED STATE ###

variable tdamp equal ${timestep}*100

fix nvt all nvt temp ${T} ${T} ${tdamp}

run 10000

reset_timestep 0

unfix nvt

fix nve all nve

<nowiki>###</nowiki> MEASURE SYSTEM STATE ###

thermo_style custom etotal temp atoms vol

variable vol equal vol

variable atoms equal atoms

variable temp equal temp

variable temp2 equal temp*temp

variable N2 equal atoms*atoms

variable E2 equal etotal*etotal

variable E equal etotal

fix aves all ave/time 100 1000 100000 v_temp2 v_E v_E2

run 100000

variable heatcapac equal
${N2}*(f_aves[3]-(f_aves[2]*f_aves[2]))/(f_aves[1])

print "Averages"

print "--------"

print "Heat Capacity: ${heatcapac}"

print "Volume: ${vol}"

print "No of atoms: ${atoms}"

Rep:SMliqsim

2018-01-31T09:06:31Z

Sm6415:

= Liquid Simulations Computational Lab Assignment - Sam Macer =

== Abstract ==

== Introduction ==
Diffusion is a fundamental form of transport in fluid and soft-condensed phases. All un-equilibrated fluid systems will move via diffusion toward thermodynamic equilibrium. Many industrial chemical processes could benefit from an optimisation considering diffusion involving the fluid or soft-condensed system being processed. For example, consider piping different coloured icings onto a cupcake to form a pattern. The manufacturer may wish to
prevent the diffusion of the two colouring molecules between the bulk of the icing, as this would ruin the pattern. Studying the temperature dependence of
this diffusion could allow them to find the maximum temperature at which the cupcake
could be safely packaged without significant diffusion occurring.

One study on the extraction of basil essential oil from the
<em>Ocimum basilicum</em> leaf exemplifies this<sup>1</sup>. The diffusion coefficient was
found to change over the course of the extraction; the vapour pressure of the
basil oil in the surrounding air increases, while the homogeneity of oil
distribution in the leaf diminished as it is depleted from different areas at
different rates, due to different tissue structures inside the leaf. The study used models based on Fick’s second law of diffusion, and ultimately elucidated a
quadratic relationship between the concentration and diffusion coefficient for
the system. This knowledge allows manufacturers to fit the model
mathematically using parameters specific to their distillation setup, and use the result to find the best trade-off between extraction time (time and energy are expensive) and extraction yield.

== Aims & Objectives ==
The aim of this experiment is to use computational models of particles interacting according to a Lennard-Jones potential to extract thermodynamic data about the system and ultimately calculate the diffusion coefficient. The experiment uses boxes of between 1000 and
10000 atoms, with periodic boundary conditions applied to simulate a much larger 'effective system'. The velocity-Verlet
algorithm<sup>2</sup> is used, which iterates the velocity and displacement for every particle in the system, from
an initial state, through the desired number of steps. At each of these steps, thermodynamic
state variables of the system can be extracted, for e.g. Temperature can be obtained using the average kinetic energies in conjunction with the equipartition theorem.

Reduced units throughout the experiment to make the magnitudes of the quantities used more managable.
Given the Lennard-Jones potential:

<div align="center">
[[File:Smliq1.png|188x188px]]

r = r*/σ

E = E*/ε

T = kT*/ε
</div>

''r = inter-particle distance, σ = inter-particle distance when the potential = 0, E = potential, ε = depth of the potential well, T = temperature and k = the Boltzmann constant. Starred quantities represent the non-reduced
versions. Note the unites cancel out when forming these quantities and hence they are unitless.

By cancelling through these reduced units, we obtain that time is also in reduced units and can be converted via:

This is used for calculating the diffusion coefficient, though may be used to obtain the true timescale for any of the simulations.

== Methods ==
The open source software, LAMMPS, was used to run the molecular dynamics simulations. The experiment was run as 5 separate parts:
# The effect of the time step on the simulation was probed using a box of 1000 atoms, equally spaced with a number density of 0.8 atoms per unit volume. The ensemble was treated as microcanonical. The following time steps were trialed: 0.0150, 0.0100, 0.0075, 0.0025, 0.0010. Energy, temperature and pressure for the system were collected at each iteration. A timestep of 0.0025 was used for the two subsequent simulations - see results and discussion for justification.
# A set of simulations were run in the isobaric-isothermal (NpT) ensemble. At P = 2.5 & P = 3.0, the simulation was run over the following temperatures: T = 1.5, 2.0, 2.5, 3.5, 5. Density of the systems was extracted as an average after the simulation, and plotted against temperature to test the agreement with the ideal gas law.
# An experiment was run in the NVT ensemble at densities 0.2 & 0.8, over the temperatures T = 2.0, 2.2, 2.4, 2.6, 2.8. The average heat capacity and temperature were collected for each system.
# Simulations were performed for a solid, a liquid and a gas system. A timestep of 0.0020 was used. Parameters for this were selected based on a previous probe into the phase transitions of the Lennard-Jones system<sup>3</sup>; solid (density = 1.2, T = 0.5), liquid (density = 0.8, T = 1.2), gas (density = 0.2, T = 1.2). The output trajectories were used to calculate the RDF (radial distribution function), i.e. the radial frequency of meeting another particle averaged over all spherical angles.
# Simulations were again performed for a solid, a liquid and a gas phase system. A timestep of 0.0020 was used for this experiment, and the same density/temperature parameters at experiment 4. The MSD (mean squared displacement) was collected as a function of timestep.

== Results & Discussion ==

=== Experiment 1 ===
<div align="center">
[[File:Smliq2.png|404x404px]]

[[File:Smliq3.png|406x406px]]

[[File:Smliq4.png|468x468px]]

These plots show how the thermodynamic state variables of the system vary with time. For all of the variables , we see a very quick convergence from the initial value to the steady state value. A closer inspection of these graphs reveal
equilibration happens after roughly 0.5 reduced time units. This shows how the system which is initially ordered quickly descends into chaos. After this point, the state variables proceed to oscillate randomly as a noise function. This occurs due to fluctuations in the system which is constantly moving. If the system had more particles, the average would be less noisy because there are more opportunities for the random imbalances to cancel each other out.

[[File:Smliq5.png|514x514px]]

The above plot shows how the energy varies with timestep. As the time step is increased the average energy
increases due to the resolution of the numerical solution decreasing. The largest timestep shows a clear divergence. Presumably the timestep is so large that particles are getting unrealistically close to each other. The reaction in the next timestep cascades such that the problem is transferred to neighbouring particles. At this point the simulation breaks down as the energies and particle positions become increasingly unrealistic. 0.0025 was the largest timestep which showed no sign of this increased energy or divergence, and hence the timestep used for all subsequent experiments was smaller than this.

</div>

=== Experiment 2 ===
<div align="center">

[[File:Smliq6.png|501x501px]]

The density decreases with temperature as expected; the
atoms have more kinetic energy, so push each other apart further during their
random motion. The higher pressure plots show a consistently higher density.
This is intuitive as at higher pressures, the atoms are squashed closer
together on average.

At lower temperatures, the ideal gas law result diverges
from the simulation. This is because the ideal gas law assumes no collisions.
At lower temperatures, particles of an ideal gas continue to be compressed into
a smaller space due to their diminishing kinetic energy, such that they are in
fact overlapping. In the simulation the particles collide, meaning a much
greater kinetic energy (temperature) reduction is required to achieve the same
increase in density. Also note that the ideal gas law allows for infinity as the
limiting density, whereas this becomes infinitely unlikely when using Van Der
Waals potential between the atoms.

</div>

=== Experiment 3 ===
<div align="center">
[[File:Smliq7.png|486x486px]]

The heat capacity is smaller for lower densities
as expected: there are less atoms to sink energy into, so the average kinetic energy
increases more quickly. The heat capacity decreases with temperature because at higher temperatures, the energy states are closer together, meaning it is easier to populate new states (higher density of states). This clearly outweighs any increased degeneracy in states, otherwise we would see the heat capacity increase with temperature as is often the case. The plateau could indicate a first order phase transition where the heat capacity tends to infinity. This would be made more clear with a higher resolution of data points. The radial distribution function could be calculated either side as the less ordered phase would have a smoother curve.

</div>

=== Experiment 4 ===
<div align="center">
[[File:Smliq8.png|478x478px]]

The radial distribution function (RDF) represents the density of particles around a single particle in the system on average as a function of radius, averaged over all possible angles in a sphere. The RDFs for the different phases are distinct. In the solid case, we se sharp peaks corresponding to when we hit a site in the crystalline lattice. In the liquid, the peaks are much less sharp, as there are many more positions the molecules can occupy easily relative to one another. Peaks are still defined as the distance between particles is roughly constant. A gas behaves like a liquid with very little order. It is less peaks if any (occur at higher densities) are much less well defined. There is still a sharp peak close to the particle as neighbouring particles spend the most time when changing direction during a collision as this is when velocity is the smallest. All three phases exhibit a cut-off displacement due to the steep positive region in the Lennard-Jones potential.

[[File:Smliq9.png|533x533px]]

When the RDF is viewed cumulatively, we see the steepest cumulation for the solid, followed by the liquid then gas. This is due to meeting particles more frequently when in the more dense phase.

</div>

=== Experiment 5 ===
<div align="center">

[[File:Smliq10.png|687x687px]]

The mean squared displacement plot used the same colours as before to represent a solid, liquid and gas. Initially the MSD is slow to rise, in particular for the gas. This is due to the particles taking some time to disperse properly and reach equilibrium velocity. The gas has the most freedom and kinetic energy so the MSD is the greatest. The opposite is true for the solid. The liquid is closer to the solid in terms of MSD because it is still relatively ordered, and there is some kind of activation barrier to molecules swapping between cages. The gradient for a solid is effectively zero as swapping positions in a crystalline lattice is extremely rare.

The second plot shows a method of least squared fit, taking into account points after timestep = 1000 only, to avoid data from the equilibration period being considered. The gradients were as following:

Solid: 2.94 x10<sup>-8</sup>

Liquid: 1.02 x10<sup>-3</sup>

Gas: 1.00 x10<sup>-2</sup>

The MSD is related to the diffusion coefficient by a simple relationship:

[[File:Smliq11.png]]

Now it has been show we can model
the diffusion coefficient, similar simulations could be applied to real world
applications. For e.g. in the basil leaf example, we could model the leaf
structure with different densities, and see how the diffusion of basil
essential oil molecules changes as a vapour shell forms and molecules in the
leaf are depleted.
</div>

== Conclusion ==

== References ==
# J. Silveira, A. Costa and E. Costa Junior, ''Engenharia Agrícola'', 2017, 37, 717-726.
# L. Verlet, ''Physical Review'', 1967, 159, 98-103.
# J. Hansen and L. Verlet, ''Physical Review'', 1969, 184, 151-161.

Weblinks:
# <nowiki>http://www.scielo.br/scielo.php?script=sci_arttext&pid=S0100-69162017000400717#B21</nowiki>
# <nowiki>https://journals.aps.org/pr/abstract/10.1103/PhysRev.159.98</nowiki>
# <nowiki>https://journals.aps.org/pr/abstract/10.1103/PhysRev.184.151</nowiki>

== Appendix ==

LAMMPS script for heat capacity experiments

Heat capacity LAMMPS input:

<nowiki>###</nowiki> DEFINE SIMULATION BOX GEOMETRY ###

variable dens equal 0.2

lattice sc ${dens}

region box block 0 15 0 15 0 15

create_box 1 box

create_atoms 1 box

<nowiki>###</nowiki> 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

<nowiki>###</nowiki> SPECIFY THE REQUIRED THERMODYNAMIC STATE ###

variable T equal 2.0

variable timestep equal 0.0075

<nowiki>###</nowiki> ASSIGN ATOMIC VELOCITIES ###

velocity all create ${T} 12345 dist gaussian rot yes mom yes

<nowiki>###</nowiki> SPECIFY ENSEMBLE ###

timestep ${timestep}

fix nve all nve

<nowiki>###</nowiki> THERMODYNAMIC OUTPUT CONTROL ###

thermo_style custom time etotal temp press

thermo 10

<nowiki>###</nowiki> RECORD TRAJECTORY ###

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

<nowiki>###</nowiki> SPECIFY TIMESTEP ###

<nowiki>###</nowiki> RUN SIMULATION TO MELT CRYSTAL ###

run 10000

unfix nve

reset_timestep 0

<nowiki>###</nowiki> BRING SYSTEM TO REQUIRED STATE ###

variable tdamp equal ${timestep}*100

fix nvt all nvt temp ${T} ${T} ${tdamp}

run 10000

reset_timestep 0

unfix nvt

fix nve all nve

<nowiki>###</nowiki> MEASURE SYSTEM STATE ###

thermo_style custom etotal temp atoms vol

variable vol equal vol

variable atoms equal atoms

variable temp equal temp

variable temp2 equal temp*temp

variable N2 equal atoms*atoms

variable E2 equal etotal*etotal

variable E equal etotal

fix aves all ave/time 100 1000 100000 v_temp2 v_E v_E2

run 100000

variable heatcapac equal
${N2}*(f_aves[3]-(f_aves[2]*f_aves[2]))/(f_aves[1])

print "Averages"

print "--------"

print "Heat Capacity: ${heatcapac}"

print "Volume: ${vol}"

print "No of atoms: ${atoms}"

2018-01-31T08:51:47Z

Sm6415:

Rep:SMliqsim

2018-01-30T23:38:41Z

Sm6415:

= Liquid Simulations Computational Lab Assignment - Sam Macer =

== Abstract ==
Time units converstion

== Introduction ==
Diffusion is a fundamental form of transport in fluid and soft-condensed phases. All un-equilibrated fluid systems will move via diffusion toward thermodynamic equilibrium. Many industrial chemical processes could benefit from an optimisation considering diffusion involving the fluid or soft-condensed system being processed. For example, consider piping different coloured icings onto a cupcake to form a pattern. The manufacturer may wish to
prevent the diffusion of the two colouring molecules between the bulk of the icing, as this would ruin the pattern. Studying the temperature dependence of
this diffusion could allow them to find the maximum temperature at which the cupcake
could be safely packaged without significant diffusion occurring.

One study on the extraction of basil essential oil from the
<em>Ocimum basilicum</em> leaf exemplifies this<sup>1</sup>. The diffusion coefficient was
found to change over the course of the extraction; the vapour pressure of the
basil oil in the surrounding air increases, while the homogeneity of oil
distribution in the leaf diminished as it is depleted from different areas at
different rates, due to different tissue structures inside the leaf. The study used models based on Fick’s second law of diffusion, and ultimately elucidated a
quadratic relationship between the concentration and diffusion coefficient for
the system. This knowledge allows manufacturers to fit the model
mathematically using parameters specific to their distillation setup, and use the result to find the best trade-off between extraction time (time and energy are expensive) and extraction yield.

== Aims & Objectives ==
The aim of this experiment is to use computational models of particles interacting according to a Lennard-Jones potential to extract thermodynamic data about the system and ultimately calculate the diffusion coefficient. The experiment uses boxes of between 1000 and
10000 atoms, with periodic boundary conditions applied to simulate a much larger 'effective system'. The velocity-Verlet
algorithm<sup>2</sup> is used, which iterates the velocity and displacement for every particle in the system, from
an initial state, through the desired number of steps. At each of these steps, thermodynamic
state variables of the system can be extracted, for e.g. Temperature can be obtained using the average kinetic energies in conjunction with the equipartition theorem.

Reduced units throughout the experiment to make the magnitudes of the quantities used more managable.
Given the Lennard-Jones potential:

<div align="center">
[[File:Smliq1.png|188x188px]]

r = r*/σ

E = E*/ε

T = kT*/ε
</div>

''r = inter-particle distance, σ = inter-particle distance when the potential = 0, E = potential, ε = depth of the potential well, T = temperature and k = the Boltzmann constant. Starred quantities represent the non-reduced
versions. Note the unites cancel out when forming these quantities and hence they are unitless.

== Methods ==
The open source software, LAMMPS, was used to run the molecular dynamics simulations. The experiment was run as 5 separate parts:
# The effect of the time step on the simulation was probed using a box of 1000 atoms, equally spaced with a number density of 0.8 atoms per unit volume. The ensemble was treated as microcanonical. The following time steps were trialed: 0.0150, 0.0100, 0.0075, 0.0025, 0.0010. Energy, temperature and pressure for the system were collected at each iteration. A timestep of 0.0025 was used for the two subsequent simulations - see results and discussion for justification.
# A set of simulations were run in the isobaric-isothermal (NpT) ensemble. At P = 2.5 & P = 3.0, the simulation was run over the following temperatures: T = 1.5, 2.0, 2.5, 3.5, 5. Density of the systems was extracted as an average after the simulation, and plotted against temperature to test the agreement with the ideal gas law.
# An experiment was run in the NVT ensemble at densities 0.2 & 0.8, over the temperatures T = 2.0, 2.2, 2.4, 2.6, 2.8. The average heat capacity and temperature were collected for each system.
# Simulations were performed for a solid, a liquid and a gas system. A timestep of 0.0020 was used. Parameters for this were selected based on a previous probe into the phase transitions of the Lennard-Jones system<sup>3</sup>; solid (density = 1.2, T = 0.5), liquid (density = 0.8, T = 1.2), gas (density = 0.2, T = 1.2). The output trajectories were used to calculate the RDF (radial distribution function), i.e. the radial frequency of meeting another particle averaged over all spherical angles.
# Simulations were again performed for a solid, a liquid and a gas phase system. A timestep of 0.0020 was used for this experiment, and the same density/temperature parameters at experiment 4. The MSD (mean squared displacement) was collected as a function of timestep.

== Results & Discussion ==

=== Experiment 1 ===
<div align="center">
[[File:Smliq2.png|404x404px]]

[[File:Smliq3.png|406x406px]]

[[File:Smliq4.png|468x468px]]

These plots show how the thermodynamic state variables of the system vary with time. For all of the variables , we see a very quick convergence from the initial value to the steady state value. A closer inspection of these graphs reveal
equilibration happens after roughly 0.5 reduced time units. This shows how the system which is initially ordered quickly descends into chaos. After this point, the state variables proceed to oscillate randomly as a noise function. This occurs due to fluctuations in the system which is constantly moving. If the system had more particles, the average would be less noisy because there are more opportunities for the random imbalances to cancel each other out.

[[File:Smliq5.png|514x514px]]

The above plot shows how the energy varies with timestep. As the time step is increased the average energy
increases due to the resolution of the numerical solution decreasing. The largest timestep shows a clear divergence. Presumably the timestep is so large that particles are getting unrealistically close to each other. The reaction in the next timestep cascades such that the problem is transferred to neighbouring particles. At this point the simulation breaks down as the energies and particle positions become increasingly unrealistic. 0.0025 was the largest timestep which showed no sign of this increased energy or divergence, and hence the timestep used for all subsequent experiments was smaller than this.

</div>

=== Experiment 2 ===
<div align="center">

[[File:Smliq6.png|501x501px]]

The density decreases with temperature as expected; the
atoms have more kinetic energy, so push each other apart further during their
random motion. The higher pressure plots show a consistently higher density.
This is intuitive as at higher pressures, the atoms are squashed closer
together on average.

At lower temperatures, the ideal gas law result diverges
from the simulation. This is because the ideal gas law assumes no collisions.
At lower temperatures, particles of an ideal gas continue to be compressed into
a smaller space due to their diminishing kinetic energy, such that they are in
fact overlapping. In the simulation the particles collide, meaning a much
greater kinetic energy (temperature) reduction is required to achieve the same
increase in density. Also note that the ideal gas law allows for infinity as the
limiting density, whereas this becomes infinitely unlikely when using Van Der
Waals potential between the atoms.

</div>

=== Experiment 3 ===
<div align="center">
[[File:Smliq7.png|486x486px]]

The heat capacity is smaller for lower densities
as expected: there are less atoms to sink energy into, so the average kinetic energy
increases more quickly. The heat capacity decreases with temperature because at higher temperatures, the energy states are closer together, meaning it is easier to populate new states (higher density of states). This clearly outweighs any increased degeneracy in states, otherwise we would see the heat capacity increase with temperature as is often the case. The plateau could indicate a first order phase transition where the heat capacity tends to infinity. This would be made more clear with a higher resolution of data points. The radial distribution function could be calculated either side as the less ordered phase would have a smoother curve.

</div>

=== Experiment 4 ===
<div align="center">
[[File:Smliq8.png|478x478px]]

The radial distribution function (RDF) represents the density of particles around a single particle in the system on average as a function of radius, averaged over all possible angles in a sphere. The RDFs for the different phases are distinct. In the solid case, we se sharp peaks corresponding to when we hit a site in the crystalline lattice. In the liquid, the peaks are much less sharp, as there are many more positions the molecules can occupy easily relative to one another. Peaks are still defined as the distance between particles is roughly constant. A gas behaves like a liquid with very little order. It is less peaks if any (occur at higher densities) are much less well defined. There is still a sharp peak close to the particle as neighbouring particles spend the most time when changing direction during a collision as this is when velocity is the smallest. All three phases exhibit a cut-off displacement due to the steep positive region in the Lennard-Jones potential.

[[File:Smliq9.png|533x533px]]

When the RDF is viewed cumulatively, we see the steepest cumulation for the solid, followed by the liquid then gas. This is due to meeting particles more frequently when in the more dense phase.

</div>

=== Experiment 5 ===
<div align="center">

[[File:Smliq10.png|687x687px]]

</div>

== Conclusion ==

== References ==
# J. Silveira, A. Costa and E. Costa Junior, ''Engenharia Agrícola'', 2017, 37, 717-726.
# L. Verlet, ''Physical Review'', 1967, 159, 98-103.
# J. Hansen and L. Verlet, ''Physical Review'', 1969, 184, 151-161.

Weblinks:
# <nowiki>http://www.scielo.br/scielo.php?script=sci_arttext&pid=S0100-69162017000400717#B21</nowiki>
# <nowiki>https://journals.aps.org/pr/abstract/10.1103/PhysRev.159.98</nowiki>
# <nowiki>https://journals.aps.org/pr/abstract/10.1103/PhysRev.184.151</nowiki>

== Appendix ==

LAMMPS script for heat capacity experiments

Heat capacity LAMMPS input:

<nowiki>###</nowiki> DEFINE SIMULATION BOX GEOMETRY ###

variable dens equal 0.2

lattice sc ${dens}

region box block 0 15 0 15 0 15

create_box 1 box

create_atoms 1 box

<nowiki>###</nowiki> 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

<nowiki>###</nowiki> SPECIFY THE REQUIRED THERMODYNAMIC STATE ###

variable T equal 2.0

variable timestep equal 0.0075

<nowiki>###</nowiki> ASSIGN ATOMIC VELOCITIES ###

velocity all create ${T} 12345 dist gaussian rot yes mom yes

<nowiki>###</nowiki> SPECIFY ENSEMBLE ###

timestep ${timestep}

fix nve all nve

<nowiki>###</nowiki> THERMODYNAMIC OUTPUT CONTROL ###

thermo_style custom time etotal temp press

thermo 10

<nowiki>###</nowiki> RECORD TRAJECTORY ###

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

<nowiki>###</nowiki> SPECIFY TIMESTEP ###

<nowiki>###</nowiki> RUN SIMULATION TO MELT CRYSTAL ###

run 10000

unfix nve

reset_timestep 0

<nowiki>###</nowiki> BRING SYSTEM TO REQUIRED STATE ###

variable tdamp equal ${timestep}*100

fix nvt all nvt temp ${T} ${T} ${tdamp}

run 10000

reset_timestep 0

unfix nvt

fix nve all nve

<nowiki>###</nowiki> MEASURE SYSTEM STATE ###

thermo_style custom etotal temp atoms vol

variable vol equal vol

variable atoms equal atoms

variable temp equal temp

variable temp2 equal temp*temp

variable N2 equal atoms*atoms

variable E2 equal etotal*etotal

variable E equal etotal

fix aves all ave/time 100 1000 100000 v_temp2 v_E v_E2

run 100000

variable heatcapac equal
${N2}*(f_aves[3]-(f_aves[2]*f_aves[2]))/(f_aves[1])

print "Averages"

print "--------"

print "Heat Capacity: ${heatcapac}"

print "Volume: ${vol}"

print "No of atoms: ${atoms}"

Rep:SMliqsim

2018-01-30T23:23:44Z

Sm6415:

= Liquid Simulations Computational Lab Assignment - Sam Macer =

== Abstract ==
Time units converstion

== Introduction ==
Diffusion is a fundamental form of transport in fluid and soft-condensed phases. All un-equilibrated fluid systems will move via diffusion toward thermodynamic equilibrium. Many industrial chemical processes could benefit from an optimisation considering diffusion involving the fluid or soft-condensed system being processed. For example, consider piping different coloured icings onto a cupcake to form a pattern. The manufacturer may wish to
prevent the diffusion of the two colouring molecules between the bulk of the icing, as this would ruin the pattern. Studying the temperature dependence of
this diffusion could allow them to find the maximum temperature at which the cupcake
could be safely packaged without significant diffusion occurring.

One study on the extraction of basil essential oil from the
<em>Ocimum basilicum</em> leaf exemplifies this<sup>1</sup>. The diffusion coefficient was
found to change over the course of the extraction; the vapour pressure of the
basil oil in the surrounding air increases, while the homogeneity of oil
distribution in the leaf diminished as it is depleted from different areas at
different rates, due to different tissue structures inside the leaf. The study used models based on Fick’s second law of diffusion, and ultimately elucidated a
quadratic relationship between the concentration and diffusion coefficient for
the system. This knowledge allows manufacturers to fit the model
mathematically using parameters specific to their distillation setup, and use the result to find the best trade-off between extraction time (time and energy are expensive) and extraction yield.

== Aims & Objectives ==
The aim of this experiment is to use computational models of particles interacting according to a Lennard-Jones potential to extract thermodynamic data about the system and ultimately calculate the diffusion coefficient. The experiment uses boxes of between 1000 and
10000 atoms, with periodic boundary conditions applied to simulate a much larger 'effective system'. The velocity-Verlet
algorithm<sup>2</sup> is used, which iterates the velocity and displacement for every particle in the system, from
an initial state, through the desired number of steps. At each of these steps, thermodynamic
state variables of the system can be extracted, for e.g. Temperature can be obtained using the average kinetic energies in conjunction with the equipartition theorem.

Reduced units throughout the experiment to make the magnitudes of the quantities used more managable.
Given the Lennard-Jones potential:

<div align="center">
[[File:Smliq1.png|188x188px]]

r = r*/σ

E = E*/ε

T = kT*/ε
</div>

''r = inter-particle distance, σ = inter-particle distance when the potential = 0, E = potential, ε = depth of the potential well, T = temperature and k = the Boltzmann constant. Starred quantities represent the non-reduced
versions. Note the unites cancel out when forming these quantities and hence they are unitless.

== Methods ==
The open source software, LAMMPS, was used to run the molecular dynamics simulations. The experiment was run as 5 separate parts:
# The effect of the time step on the simulation was probed using a box of 1000 atoms, equally spaced with a number density of 0.8 atoms per unit volume. The ensemble was treated as microcanonical. The following time steps were trialed: 0.0150, 0.0100, 0.0075, 0.0025, 0.0010. Energy, temperature and pressure for the system were collected at each iteration. A timestep of 0.0025 was used for the two subsequent simulations - see results and discussion for justification.
# A set of simulations were run in the isobaric-isothermal (NpT) ensemble. At P = 2.5 & P = 3.0, the simulation was run over the following temperatures: T = 1.5, 2.0, 2.5, 3.5, 5. Density of the systems was extracted as an average after the simulation, and plotted against temperature to test the agreement with the ideal gas law.
# An experiment was run in the NVT ensemble at densities 0.2 & 0.8, over the temperatures T = 2.0, 2.2, 2.4, 2.6, 2.8. The average heat capacity and temperature were collected for each system.
# Simulations were performed for a solid, a liquid and a gas system. A timestep of 0.0020 was used. Parameters for this were selected based on a previous probe into the phase transitions of the Lennard-Jones system<sup>3</sup>; solid (density = 1.2, T = 0.5), liquid (density = 0.8, T = 1.2), gas (density = 0.2, T = 1.2). The output trajectories were used to calculate the RDF (radial distribution function), i.e. the radial frequency of meeting another particle averaged over all spherical angles.
# Simulations were again performed for a solid, a liquid and a gas phase system. A timestep of 0.0020 was used for this experiment, and the same density/temperature parameters at experiment 4. The MSD (mean squared displacement) was collected as a function of timestep.

== Results & Discussion ==

=== Experiment 1 ===
<div align="center">
[[File:Smliq2.png|404x404px]]

[[File:Smliq3.png|406x406px]]

[[File:Smliq4.png|468x468px]]

These plots show how the thermodynamic state variables of the system vary with time. For all of the variables , we see a very quick convergence from the initial value to the steady state value. A closer inspection of these graphs reveal
equilibration happens after roughly 0.5 reduced time units. This shows how the system which is initially ordered quickly descends into chaos. After this point, the state variables proceed to oscillate randomly as a noise function. This occurs due to fluctuations in the system which is constantly moving. If the system had more particles, the average would be less noisy because there are more opportunities for the random imbalances to cancel each other out.

[[File:Smliq5.png|562x562px]]

The above plot shows how the energy varies with timestep. As the time step is increased the average energy
increases due to the resolution of the numerical solution decreasing. The largest timestep shows a clear divergence. Presumably the timestep is so large that particles are getting unrealistically close to each other. The reaction in the next timestep cascades such that the problem is transferred to neighbouring particles. At this point the simulation breaks down as the energies and particle positions become increasingly unrealistic. 0.0025 was the largest timestep which showed no sign of this increased energy or divergence, and hence the timestep used for all subsequent experiments was smaller than this.

</div>

=== Experiment 2 ===
<div align="center">

[[File:Smliq6.png|780x780px]]

The density decreases with temperature as expected; the
atoms have more kinetic energy, so push each other apart further during their
random motion. The higher pressure plots show a consistently higher density.
This is intuitive as at higher pressures, the atoms are squashed closer
together on average.

At lower temperatures, the ideal gas law result diverges
from the simulation. This is because the ideal gas law assumes no collisions.
At lower temperatures, particles of an ideal gas continue to be compressed into
a smaller space due to their diminishing kinetic energy, such that they are in
fact overlapping. In the simulation the particles collide, meaning a much
greater kinetic energy (temperature) reduction is required to achieve the same
increase in density. Also note that the ideal gas law allows for infinity as the
limiting density, whereas this becomes infinitely unlikely when using Van Der
Waals potential between the atoms.

</div>

=== Experiment 3 ===
<div align="center">
[[File:Smliq7.png|695x695px]]

The heat capacity is smaller for lower densities
as expected: there are less atoms to sink energy into, so the average kinetic energy
increases more quickly. The heat capacity decreases with temperature because at higher temperatures, the energy states are closer together, meaning it is easier to populate new states (higher density of states). This clearly outweighs any increased degeneracy in states, otherwise we would see the heat capacity increase with temperature as is often the case. The plateau could indicate a first order phase transition where the heat capacity tends to infinity. This would be made more clear with a higher resolution of data points. The radial distribution function could be calculated either side as the less ordered phase would have a smoother curve.

</div>

=== Experiment 4 ===
<div align="center">
[[File:Smliq8.png|531x531px]]
[[File:Smliq9.png|555x555px]]

</div>

=== Experiment 5 ===
<div align="center">

[[File:Smliq10.png|687x687px]]

</div>
== References ==
# J. Silveira, A. Costa and E. Costa Junior, ''Engenharia Agrícola'', 2017, 37, 717-726.
# L. Verlet, ''Physical Review'', 1967, 159, 98-103.
# J. Hansen and L. Verlet, ''Physical Review'', 1969, 184, 151-161.

Weblinks:
# <nowiki>http://www.scielo.br/scielo.php?script=sci_arttext&pid=S0100-69162017000400717#B21</nowiki>
# <nowiki>https://journals.aps.org/pr/abstract/10.1103/PhysRev.159.98</nowiki>
# <nowiki>https://journals.aps.org/pr/abstract/10.1103/PhysRev.184.151</nowiki>

== Appendix ==

LAMMPS script for heat capacity experiments

Heat capacity LAMMPS input:

<nowiki>###</nowiki> DEFINE SIMULATION BOX GEOMETRY ###

variable dens equal 0.2

lattice sc ${dens}

region box block 0 15 0 15 0 15

create_box 1 box

create_atoms 1 box

<nowiki>###</nowiki> 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

<nowiki>###</nowiki> SPECIFY THE REQUIRED THERMODYNAMIC STATE ###

variable T equal 2.0

variable timestep equal 0.0075

<nowiki>###</nowiki> ASSIGN ATOMIC VELOCITIES ###

velocity all create ${T} 12345 dist gaussian rot yes mom yes

<nowiki>###</nowiki> SPECIFY ENSEMBLE ###

timestep ${timestep}

fix nve all nve

<nowiki>###</nowiki> THERMODYNAMIC OUTPUT CONTROL ###

thermo_style custom time etotal temp press

thermo 10

<nowiki>###</nowiki> RECORD TRAJECTORY ###

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

<nowiki>###</nowiki> SPECIFY TIMESTEP ###

<nowiki>###</nowiki> RUN SIMULATION TO MELT CRYSTAL ###

run 10000

unfix nve

reset_timestep 0

<nowiki>###</nowiki> BRING SYSTEM TO REQUIRED STATE ###

variable tdamp equal ${timestep}*100

fix nvt all nvt temp ${T} ${T} ${tdamp}

run 10000

reset_timestep 0

unfix nvt

fix nve all nve

<nowiki>###</nowiki> MEASURE SYSTEM STATE ###

thermo_style custom etotal temp atoms vol

variable vol equal vol

variable atoms equal atoms

variable temp equal temp

variable temp2 equal temp*temp

variable N2 equal atoms*atoms

variable E2 equal etotal*etotal

variable E equal etotal

fix aves all ave/time 100 1000 100000 v_temp2 v_E v_E2

run 100000

variable heatcapac equal
${N2}*(f_aves[3]-(f_aves[2]*f_aves[2]))/(f_aves[1])

print "Averages"

print "--------"

print "Heat Capacity: ${heatcapac}"

print "Volume: ${vol}"

print "No of atoms: ${atoms}"

Rep:SMliqsim

2017-11-22T11:35:17Z

Sm6415:

Rep:Sm ising

2017-11-22T11:17:18Z

Sm6415:

= Ising Model Report - Sam Macer =

== The Ising model, energy and magnetisation ==
The Ising model is a simple model which examines the magnetic phase behaviour of a system. The system is set up as an n-dimensional lattice of particles which are either spin up (+1) or spin down (-1). The interaction between any two particles that share a face is counted as an additive contribution towards the overall energy of the system. Each interaction is calculated by multiplying the two spins and then multiplying by some constant J to convert units to energy.

'''Task 1:'''

In a lowest energy configuration, all spins will be oriented in the same direction (+1 or -1) to maximise the magnitude of the sum or the interacting spin pairs.

Per dimension, there are two directions in which each spin can interact given that only neighbours sharing a face can interact (modelling the spins as cubes, weather arranged in a line, plane or box).

This means per particle there are 2D interactions, where D = number of dimensions in the system.

The total number of interactions is therefore 2DN/2, where we divide by 2 as each interaction is double counted when we sum interactions over all the particles, where N = number of particles in the system.

Weather each interaction is -1 with -1 or 1 with 1, the magnitude will always be 1, so the total magnitude of the energy of all the interactions will be J * 1.0 * number of interactions where J is some constant converting the units of the energy from arbitrary units into joules.

This simplifies to: Magnitude of the energy of all interactions = DNJ. We know the energy is negative as all spins in alignment correponds to the greatest energy reduction for the system, so we may conclude:

'''''Lowest possible interaction energy for a system (E) = -DNJ'''''

Only two configurations give this state: when either all spins are +1 or all spins are -1.

'''''The multiplicity is therefore given by 2 * N!/N! (= 2) where N is the number of spins in the system.'''''

Entropy (S) = k ln(M) where k is the Boltzmann constant and M is the number of equivalent microstates, so the lowest entropy of the Ising system is given by:

'''''S = k ln(2)'''''

'''Task 2:'''

From the lowest energy configuration, if one of the spins were to instantaneously flip, the new configuration would be all spins aligned but one; i.e all +1 with one -1 or all -1 with one +1.

Multiplicity therefore = 2 * N!/(N - 1)!

The entropy for the D =3, N =1000 system is therefore k ln(2 * (1000!/999!))

'''''S = k ln(2 * (1000!/999!))'''''

'''''S = k ln(2000)'''''

In terms of the energy, one spin interacts with 6 neighbours in a 3D lattice. One positive interaction changing into one negative interaction increases the energy by 2J, so overall, flipping one spin leads to a 12J increase in energy to the system.

=== Magnetisation ===
The magnetisation of the system is defined as the sum of all the spins in the system.

'''Task 3:'''

System 1:

+ + + - -

For system 1, the total magnetisation is 3 - 2:

'''''Magnetisation = + 1'''''

System 1:

+ + + + -

+ + + - -

+ + + - -

+ + - - -

+ - - - -

For system 1, the total magnetisation is 13 - 12:

'''''Magnetisation = + 1'''''

At absolute zero the entropic contribution would be minimised, so one would expect the energetic contribution to dominate and align all the spins.

'''''At absolute zero therefore, for a 3D system of 1000 spins, magnetisation = 1000'''''

== Python (Energy and magnetisation) ==
'''Task 4'''

Functions to extract the energy and magnetisation of a system were then written in python;

def energy(self):
"Return the total energy of the current lattice configuration."
# Function works by cycling through each element and then interacting with each neighbour on a grid of coordinates. The energy is halved to prevent double counting and made negative.
energy = 0.0
net_interaction = 0
y = 0
x = 0
for row in self.lattice:
for element in row:
# interact right
if x != self.n_cols - 1:
net_interaction += (element * self.lattice[y, x + 1])
else:
# if we are at the end element, interact with the adjacent system
net_interaction += (element * self.lattice[y, 0])
# interact left
if x != 0:
net_interaction += (element * self.lattice[y, x - 1])
else:
# if we are at the start element, interact with the adjacent system
net_interaction += (element * self.lattice[y, self.n_cols - 1])
# interact bottom
if y != self.n_rows - 1:
net_interaction += (element * self.lattice[y + 1, x])
else:
# if we are at the end element, interact with the adjacent system
net_interaction += (element * self.lattice[0, x])
# interact top
if y != 0:
net_interaction += (element * self.lattice[y - 1, x])
else:
# if we are at the start element, interact with the adjacent system
net_interaction += (element * self.lattice[self.n_rows - 1, x])
# iterate selected spin
if x == self.n_cols - 1:
x = 0
y += 1
else:
x += 1
interaction_energy = - 0.5 * net_interaction
return interaction_energy

def magnetisation(self):
"Return the total magnetisation of the current lattice configuration."
magnetisation = 0
for row in self.lattice:
for element in row:
magnetisation += element
return magnetisation
'''Task 5'''

The 'ILcheck.py' script was then ran multiple times to check the energy and magnetisation functions gave consistently correct results:

[[File:Ising001.png|400px]][[File:Ising002.png|400px]][[File:Ising003.png|400px]]

== Python (Monte Carlo Function) ==
'''Task 6'''

Let us consider a system with 100 spins. Considering each spin has two possible states, there are <math>2^{100} = 1.27 *10^{30}</math> potential configurations of the system. Even if we could analyse <math>10^9</math> configurations per second, it would still take us <math>1.27 *10^{21} s</math> to calculate the magnetisations of all of the configurations to take an average. Clearly this would take too long so we must look for an alternative method of gathering data on the behaviour of the energy. One way of doing this is using a Monte Carlo simulation where we start with a random configuration and randomly flip spins. We accept the new configuration if there is an energy reduction. Otherwise we accept the new configuration at a rate proportional to the Boltzmann factor as this reflects how often a particle would theoretically have enough energy to flip to a higher energy state.

'''Task 7'''

A Monte-Carlo simulation was coded in python as follows:

def montecarlostep(self, T):
"preforms a single Monte Carlo step"
energy = self.energy()
current_lattice = self.lattice
# create a new lattice with on spin flipped
#the following two lines select the coordinates of the random spin
random_i = np.random.choice(range(self.n_rows))
random_j = np.random.choice(range(self.n_cols))
self.lattice[random_i][random_j] = -self.lattice[random_i][random_j]
new_E = self.energy()
# compare energies and decided what configuration to adopt, then update parameters
dE = new_E - energy
if dE < 0:
# print('new energy lower: lattice updated')
pass
else:
# the following line will choose a random number in the range (0,1)
randNumb = np.random.random()
if randNumb <= math.exp(-dE/T):
# print('lattice was updated')
pass
else:
# print('lattice was not updated')
self.lattice[random_i][random_j] = -self.lattice[random_i][random_j]
self.E += self.energy()
self.E2 += self.energy()**2
self.M += self.magnetisation()
self.M2 += self.magnetisation()**2
# converts cycle index x to cycle number x + 1 (the next machine cycle step number)
# In other words increments the counter of number of cycles
self.n_cycles += 1
return self.energy(), self.magnetisation()

def statistics(self):
"returns statistical data collected during a Monte Carlo simulation"
avE = self.E / self.n_cycles
avE2 = self.E2 / self.n_cycles
avM = self.M / self.n_cycles
avM2 = self.M2 / self.n_cycles
N = self.n_cycles
return avE, avE2, avM, avM2, N
'''Task 8'''

When the temperature of the system is above the Curie Temperature (Tc), the system is paramagnetic at thermodynamic equilibrium, i.e. the spins are oriented randomly in the absence of an external magnetic field. When the temperature of the system is below the Curie Temperature, the system is ferromagnetic at thermodynamic equilibrium (spins align). In the case T < Tc, one would therefore expect the system to spontaneously magnetise (spins align) as it attains the most thermodynamically stable state.

Indeed, when we run our Monte Carlo simulation for an 8x8 lattice at T = 0.5 K, it can be seen that the spins tend toward alignment in on direction as the energy minimises;

[[File:Ising004.png|center]]
Averaged quantities:
E = -1.96841896186
E*E = 3.88120282707
M = -0.991326800847
M*M = 0.983257746292

== Accelerating the code ==
'''Task 9'''

In an attempt to accelerate the code and save time, the energy() and magnetisation() functions were made more efficient using functions within the numpy python library:

def energy(self):
"Improved energy function which implements numpy's roll and multiply functions"
# multiply by a shifted grid to apply vertical energy contribution
vertical_energy = -np.sum(np.multiply(self.lattice, np.roll(self.lattice, 1, axis = 0)))
# multiply by a shifted grid to apply the horizontal energy contribution
total_energy = vertical_energy - np.sum(np.multiply(self.lattice, np.roll(self.lattice, 1, axis = 1)))
return total_energy

def magnetisation(self):
"Improved magnetisation function that implements numpy's sum function:"
return np.sum(self.lattice)

'''Task 10'''

Time trials were used to compare the speed of the old and updated energy functions:

New function:
2.4918089999999893 s
2.4761589999999956 s
2.540256999999997 s
2.4700420000000065 s
2.5968489999999917 s

Average = 2.52 s
Standard error of average = 0.024 s

Old function:
10.844294 s
10.496529000000002 s
10.151754999999998 s
10.489855999999996 s
9.988599999999998 s

Average = 10.39 s
Standard error of average = 0.149 s

The result is the improved function only takes 24.3 % of the time took by the old function to do the same thing.
==The effect of temperature==
The following plots were produced by running 150000 steps of the Monte Carlo simulation at T = 0.5 K and T = 2.0 K respectively:

[[File:Ising005.png|500px]][[File:Ising006.png|500px]]

[[File:Ising007.png|500px]][[File:Ising008.png|500px]]

Note, in the second system, we see that there must be a realistic chance for the whole system to flip the net spin direction as this is what can be seen from the graph to occur. This happens because in this system there are two global minimum points on the potential energy surface separated by a barrier. At this temperature there is enough energy to only occasionally overcome this barrier.

'''Task 11'''

Looking at the dynamics of the systems, we can see in the temperature region in which we are working, the region after 3000 steps should be a sufficient representation of the system at equilibrium. The Monte Carlo step and statistics functions were updated to only begin the statistics functionality after 3000 steps:

# Defines first cycle in which statistical data is collected
statsStartStep = 3000
# Converts cycle nummber into the machine read index of the cycle (humans count from 1, machine indexes from 0)
statsStartStep -= 1

def montecarlostep(self, T):
"preforms a single Monte Carlo step"
energy = self.energy()
current_lattice = self.lattice
# create a new lattice with on spin flipped
#the following two lines select the coordinates of the random spin
random_i = np.random.choice(range(self.n_rows))
random_j = np.random.choice(range(self.n_cols))
self.lattice[random_i][random_j] = -self.lattice[random_i][random_j]
new_E = self.energy()
# compare energies and decided what configuration to adopt, then update parameters
dE = new_E - energy
if dE < 0:
# print('new energy lower: lattice updated')
pass
else:
# the following line will choose a random number in the range (0,1)
randNumb = np.random.random()
if randNumb <= math.exp(-dE/T):
# print('lattice was updated')
pass
else:
# print('lattice was not updated')
self.lattice[random_i][random_j] = -self.lattice[random_i][random_j]
# Update the running total IF we are in the region we are gathering stats for
if self.n_cycles >= self.statsStartStep:
self.E += self.energy()
self.E2 += self.energy()**2
self.M += self.magnetisation()
self.M2 += self.magnetisation()**2
# converts cycle index x to cycle number x + 1 (the next machine cycle step number)
self.n_cycles += 1
# print(self.n_cycles)
return self.energy(), self.magnetisation()

def statistics(self):
# Updated statistics function
if self.n_cycles >= self.statsStartStep:
avE = self.E / (self.n_cycles - self.statsStartStep)
avE2 = self.E2 / (self.n_cycles - self.statsStartStep)
avM = self.M / (self.n_cycles - self.statsStartStep)
avM2 = self.M2 / (self.n_cycles - self.statsStartStep)
N = self.n_cycles
return avE, avE2, avM, avM2, N

== The effect of system size ==
'''Task 12'''

A function was then used to plot average energy per spin and average magnetisation per spin against temperature from T = 0.25 K to T = 5 K in 0.25 K intervals. Error bars represent one standard deviation's worth of uncertainty for the calculation. This was done for 2x2, 4x4, 8x8, 16x16 and 32x32 size lattices, graphed in the same order:

[[File:Ising009.png|500px]]

2x2

[[File:Ising010.png|500px]]

4x4

[[File:Ising011.png|500px]]

8x8

[[File:Ising012.png|500px]]

16x16

[[File:Ising013.png|500px]]

32x32

All of the graphs show a convergence to a limiting value at high temperatures showing the Temperature a) ceases to affect magnetisation (the paramagnetic state) and b) subsequently ceases to affect energy stabilisation due to spin alignment (there are no more static alignment effects). The magnetisation graphs show an oscillating convergence because the spin of the hole system can flip near the Curie temperature as explained earlier. Both +1 and -1 spin systems provide the same energy stabilisation to the system, so energy does not oscillate in it's convergence.

The error bars for both quantities peak around the Curie temperate presumably as at this point the system is most often undergoing full system spin flips, causing a large variance. At the beginning of the plot the variance is low as expected because there is little dynamic activity. At the end of the plot, even though there is a lot of noise and hence wild oscillation of the spins in the system, there is so much energy that almost any change is accepted, and so the extreme configurations (i.e. all aligned) rarely occur because statistically they are so unlikely to occur. This means the variation in energy and magnetisation for the most part varies within a smaller range and so the errorbars of one standard deviation (95% confidence) are smaller.

The variation between the plots themselves can also be elucidated:

'''Task 13'''

The average energy per spin plots were superimposed. What we see is the convergence of the energy to a certain curve for large lattices:

[[File:Ising014.png|500px]]

This shows the lattice size is large enough to properly model long range effects as increasing the lattice size has little further effect.

[[File:Pl2mag.png|500px]]

The plot superimposing the magnetisations is also interesting. Firstly we see the magnitude of magnetisation per spin is pretty constant which is exactly what we would expect. Secondly we can see the oscillating convergence much more clearly. Finally we see the convergence is quickest for the smaller lattices, presumably because it's much easier kinetically to have an even distribution of spins.

== Determining the heat capacity ==
'''Task 14'''

We know the average energy is equal to the probability of each energy times the energy. In terms of the partition function (Z, Statistical thermodynamics), we can write this and rearrange to get: <math>\left\langle E\right\rangle = -\frac{\partial}{\partial \beta}logZ</math>.

Considering <math>\mathrm{Var}[E] = \left\langle E^2\right\rangle - \left\langle E\right\rangle^2</math>, we can rearrange (by equating the squared average energy in both equations) to:

<math>\mathrm{Var}[E] = -\frac{\partial \left\langle E\right\rangle}{\partial \beta}</math>

Using (The defenition of heat capacity):

<math>C = \frac{\partial \left\langle E\right\rangle}{\partial T}</math>

We obtain via integration:

<math>C = \frac{\mathrm{Var}[E]}{k_B T^2}</math>

'''Task 15'''

[[File:Ising015.png|500px]]

'''Task 16'''

Heat capacity can therefore be worked out using:
x2 = np.array(np.loadtxt('2x2.dat'))
varEx2 = x2[:,2]-x2[:,1]**2
cx2 = varEx2 / x2[:,0]**2
# Plot and make the heat capacity per spin by dividing by lattice size
pl.plot(x2[:,0], cx2/4, color = 'green')
for each lattice size where 2x2.dat is a table with columns [Temperature, Average energy, Average squared energy, Average Magnetisation, Average Squared Magnetisation]

The Heat capacity per spin was subsequently plotted as a function of temperature and lattice size:

[[File:Ising016.png|500px]]

== Finding the Curie temperature ==
'''Task 17'''
[[File:Ising017.png|500px]]

'''Task 18'''
import numpy as np
import matplotlib.pyplot as pl
data = np.loadtxt('32x32C.dat') #assume data is now a 2D array containing two columns, T and C
T = data[:,0] #get the first column
C = data[:,5] # get the second column
#first we fit the polynomial to the data
fit = np.polyfit(T, C, 3) # fit a third order polynomial
#now we generate interpolated values of the fitted polynomial over the range of our function
T_min = np.min(T)
T_max = np.max(T)
T_range = np.linspace(T_min, T_max, 1000) #generate 1000 evenly spaced points between T_min and T_max
fitted_C_values = np.polyval(fit, T_range) # use the fit object to generate the corresponding values of C
pl.plot(T_range, fitted_C_values, color = 'green')
fit50 = np.polyfit(T, C, 50)
fitted_C_values50 = np.polyval(fit50, T_range) # use the fit object to generate the corresponding values of C
pl.plot(T_range, fitted_C_values50, color = 'red')
x32C = np.array(np.loadtxt('32x32C.dat'))
pl.plot(x32C[:,0], x32C[:,5], color = 'black')
pl.ylabel('Heat capacity per spin')
pl.xlabel('Temperature')
pl.legend(['Third Order Fit', '50 Order fit', 'Raw C++ Data'])
pl.show()

[[File:Ising018.png|500px]]

'''Task 19'''

[[File:Ising019.png|500px]]

The fit was repeated for a smaller range. The peak was found to be: 1.85691773891 [ 2.29369369]

'''Task 20'''

[[File:Ising020.png|500px]]

Tc for an infinite lattice: 2.27916772 K