ChemWiki - User contributions [en]

Rep:Liquids GE715

2018-01-26T19:41:58Z

Ge715:

== Abstract ==

In this experiment properties of liquids were modelled using molecular dynamics experiments. The effect of the temperature on the density as well as the effect of temperature on the heat capacity, and the radial distribution functions were explored.

== Introduction ==

Molecular dynamics is a modern computational technique concerned with calculating properties of substances by modelling them and their behaviour on atomic level with powerful computers.

Lennard Jones potential equals <math>\phi\left(r\right) = 4\epsilon \left( \frac{\sigma^{12}}{r^{12}} - \frac{\sigma^6}{r^6} \right)</math>. Potential energy is zero at the point when <math>r=\sigma</math>, thus <math>r_0=\sigma</math>.

Force is simply the derivative of the potential, thus <math>F = 4\epsilon \left( -12 \frac{\sigma^{12}}{r^{13}} + 6 \frac{\sigma^6}{r^7} \right)</math>. At <math>r=\sigma</math> this equals <math>F = - \frac{24\epsilon}{\sigma}</math>. Well depth is <math>\phi\left(r_{eq}\right) = - \epsilon</math>, where <math>r_{eq}=\sigma \sqrt[6]{2}</math>.

The following integrals are equal to: <math>\int_{2\sigma}^\infty \phi\left(r\right)\mathrm{d}r=-0.0248</math>, <math>\int_{2.5\sigma}^\infty \phi\left(r\right)\mathrm{d}r=-0.00818</math>, and <math>\int_{3\sigma}^\infty \phi\left(r\right)\mathrm{d}r=-0.00329</math>.

=== Introductory tasks ===

==== Water ====

There are roughly <math>3.3*10^{22}</math> molecules in 1 mL of water. However, 10000 atoms of water take up about <math>3*10^{-19} mL</math> or <math>300 nm^{3}</math>. This is a good illustration of how different atomic scale is from the scale on which the properties of substances tend to be observed.

==== Atoms ====

Atom at position <math>(0.5,0.5,0.5)</math> moved for <math>(0.7,0.6,0.2)</math> would be found at <math>(0.2,0.1,0.7)</math>, since it would be copied in both x and y direction (assuming x,y,z).

==== Reduced units ====

Lennard-Jones parameters for argon were given, the macroscopic properties are calculated here.

<math>LJ_{cut-off}=1.088 nm</math>

<math>T=180 K</math>

<math>r=0.382 nm</math>

<math>E=1 \frac{kJ}{mol}</math>

== Aims and Objectives ==

The aim of this experiment was to explore the modern use of molecular dynamics for simulating simple liquids.

The objectives were to look at how timesteps used in experiment affect the quality of the results, how temperature affects density and find radial distribution functions and their integrals for different phases of matter.

== Methods ==

The calculations were performed on a high performance computer using LAMMPS software package, while the calculations of radial distribution functions were done using VMD (delta r=0.05).

== Equilibration ==

If two atoms were generated very close together that would mean the force between them would be enormous, which would have a potential to heavily disturb the simulation (very high potential energy would increase the temperature).

=== Commands ====

mass 1 1.0 signifies that the mass of the first kind of particles (1) is 1.0.

pair_style lj/cut 3.0 tells us where the interactions are cut off, as seen in the integral above, <math>3.0 \sigma</math> is a very good cut-off, since there is no significant perturbation.

The algorithm used for the calculations is velocity-Verlet, because the classical one does not use velocities.

=== Simulations ===

In the gallery below graphs are shown. It can be seen that as timesteps are increasing, the quality of results is decreasing. Timpestep of 0.015 s is wholy inappropriate, but I have chosen to run the experiments at the timestep of 0.0025 s. Results would probably also be almost as good using 0.01 s. We can see this by the fact that equilibrium is reached very quickly, for the first for timesteps in the space of 0.1 s -0.3 s. The highest timestep will mimic the reality better, but it might be too demanding computationally.

<gallery mode="packed" heights=200 px>
Image:GE715L1.PNG|Figure 1: Energy at timestep 0.001 s
Image:GE715L2.PNG|Figure 2: Temperature at timestep 0.001 s
Image:GE715L3.PNG|Figure 3: Pressure at timestep 0.001 s
Image:GE715L4.PNG|Figure 4: Energy at timestep 0.0025 s
Image:GE715L5.PNG|Figure 5: Energy at timestep 0.0075 s
Image:GE715L6.PNG|Figure 6: Energy at timestep 0.01 s
Image:GE715L7.PNG|Figure 7: Energy at timestep 0.015 s
</gallery>

== Specific conditions ==

The next task was to run simulations varying temperatures at constant pressure. This shows the dependence of density on temperature at a constant pressure. Considering the ideal gas law the density will be inversely proportional with temperature.

If we run simulations under different conditions, we get the results shown in the gallery below. Bear in mind that some errors are almost negligible.

<gallery mode="packed" heights=400 px>
Image:GE715L8.PNG|Figure 8: Density vs temperature at p=2
Image:GE715L9.PNG|Figure 9: Density vs temperature at p=3
Image:GE715L16.PNG|Figure 10: Inverse of density vs temperature at p=2
Image:GE715L17.PNG|Figure 11: Inverse of density vs temperature at p=3
</gallery>

The inverses are of densities are shown to prove that this simulation obeys the proportionality proposed ideal gas law, where the density and temperature are inversely proportional.

== Heat capacity calculation ==

Heat capacity at constant volume may be calculated using equation <math>C_V = N^2\frac{\left\langle E^2\right\rangle - \left\langle E\right\rangle^2}{k_B T^2}</math>, which LAMMPS software enables us to do. This was done at two different densities (0.2 and 0.8) and is shown in the gallery below.

<gallery mode="packed" heights=250 px>
Image:GE715L18.PNG|Figure 12: Heat capacity per volume vs temperature at d=0.2
Image:GE715L19.PNG|Figure 13: Heat capacity per volume vs temperature at d=0.8

</gallery>

Any trends in heat capacity over a range of temperatures at a fixed volume is not entirely clear even though we would normally expect a gradual decrease in heat capacity with higher temperatures.

An example of the code used to calculate this is shown below (d=0.2, T=2.0). Heat capacity was finally divided by a volume.

<pre>
### 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

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

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

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

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

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

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

### SPECIFY TIMESTEP ###

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

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

### MEASURE SYSTEM STATE ###
thermo_style custom atoms etotal temp vol
variable N2 equal atoms*atoms
variable E equal etotal
variable E2 equal etotal*etotal
variable temp equal temp
variable temp2 equal temp*temp
variable vol equal vol
fix aves all ave/time 100 1000 100000 v_E v_E2 v_temp2
run 100000

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

print "Averages"
print "--------"
print "Heat Capacity: ${heatcapac}"
print "Volume: ${vol}"

</pre>

== Radial distribution functions ==

The radial distribution functions from the calculations are given in the gallery below.

<gallery mode="packed" heights=300 px>
Image:GE715L10.PNG|Figure 14: RDF for a gas
Image:GE715L11.PNG|Figure 15: RDF for a liquid
Image:GE715L12.PNG|Figure 16: RDF for a solid
Image:GE715L13.PNG|Figure 17: Integral of RDF for a gas
Image:GE715L14.PNG|Figure 18: Integral of RDF for a liquid
Image:GE715L15.PNG|Figure 19: Integral of RDF for a solid
</gallery>

We can see from the graphs the differences between different phases and their order. It is clear that in a solid lattice structure is fairly rigid (peaks are the neighbouring atoms), liquid is slightly less ordered, while gas has almost. no real order at all. Lattice spacing in a solid is <math>\sigma</math>, while the first peak represents the nearest neighbouring atoms, while the next two represent the atom that are not nearest neighbours, but are less than two atoms away in a straight line from the central atom. The interactions are dropping of with distance, as expected in Lennard-Jones model.

== Conclusion ==

It has been determined that appropriate timestep for dynamics calculations is 0.0025 s.

Simulations were performed to determine that temperature and density of a liquid are inversely proportional.

Another finding was that heat capacity does not necessarily have a clear strong dependence on temperature at constant volume.

It was also found that radial distribution function depends heavily on the phase of matter, with order being clearly seen in solid lattice.

File:GE715L19.PNG

2018-01-26T19:29:36Z

Ge715:

File:GE715L18.PNG

2018-01-26T19:29:21Z

Ge715:

Rep:Liquids GE715

2018-01-26T18:16:20Z

Ge715:

== Abstract ==

In this experiment properties of liquids were modelled using molecular dynamics experiments. The effect of the temperature on the density as well as the radial distribution functions were explored.

== Introduction ==

Molecular dynamics is a modern computational technique concerned with calculating properties of substances by modelling them and their behaviour on atomic level with powerful computers.

Lennard Jones potential equals <math>\phi\left(r\right) = 4\epsilon \left( \frac{\sigma^{12}}{r^{12}} - \frac{\sigma^6}{r^6} \right)</math>. Potential energy is zero at the point when <math>r=\sigma</math>, thus <math>r_0=\sigma</math>.

Force is simply the derivative of the potential, thus <math>F = 4\epsilon \left( -12 \frac{\sigma^{12}}{r^{13}} + 6 \frac{\sigma^6}{r^7} \right)</math>. At <math>r=\sigma</math> this equals <math>F = - \frac{24\epsilon}{\sigma}</math>. Well depth is <math>\phi\left(r_{eq}\right) = - \epsilon</math>, where <math>r_{eq}=\sigma \sqrt[6]{2}</math>.

The following integrals are equal to: <math>\int_{2\sigma}^\infty \phi\left(r\right)\mathrm{d}r=-0.0248</math>, <math>\int_{2.5\sigma}^\infty \phi\left(r\right)\mathrm{d}r=-0.00818</math>, and <math>\int_{3\sigma}^\infty \phi\left(r\right)\mathrm{d}r=-0.00329</math>.

=== Introductory tasks ===

==== Water ====

There are roughly <math>3.3*10^{22}</math> molecules in 1 mL of water. However, 10000 atoms of water take up about <math>3*10^{-19} mL</math> or <math>300 nm^{3}</math>. This is a good illustration of how different atomic scale is from the scale on which the properties of substances tend to be observed.

==== Atoms ====

Atom at position <math>(0.5,0.5,0.5)</math> moved for <math>(0.7,0.6,0.2)</math> would be found at <math>(0.2,0.1,0.7)</math>, since it would be copied in both x and y direction (assuming x,y,z).

==== Reduced units ====

Lennard-Jones parameters for argon were given, the macroscopic properties are calculated here.

<math>LJ_{cut-off}=1.088 nm</math>

<math>T=180 K</math>

<math>r=0.382 nm</math>

<math>E=1 \frac{kJ}{mol}</math>

== Aims and Objectives ==

The aim of this experiment was to explore the modern use of molecular dynamics for simulating simple liquids.

The objectives were to look at how timesteps used in experiment affect the quality of the results, how temperature affects density and find radial distribution functions and their integrals for different phases of matter.

== Methods ==

The calculations were performed on a high performance computer using LAMMPS software package, while the calculations of radial distribution functions were done using VMD (delta r=0.05).

== Equilibration ==

If two atoms were generated very close together that would mean the force between them would be enormous, which would have a potential to heavily disturb the simulation (very high potential energy would increase the temperature).

=== Commands ====

mass 1 1.0 signifies that the mass of the first kind of particles (1) is 1.0.

pair_style lj/cut 3.0 tells us where the interactions are cut off, as seen in the integral above, <math>3.0 \sigma</math> is a very good cut-off, since there is no significant perturbation.

The algorithm used for the calculations is velocity-Verlet, because the classical one does not use velocities.

=== Simulations ===

In the gallery below graphs are shown. It can be seen that as timesteps are increasing, the quality of results is decreasing. Timpestep of 0.015 s is wholy inappropriate, but I have chosen to run the experiments at the timestep of 0.0025 s. Results would probably also be almost as good using 0.01 s. We can see this by the fact that equilibrium is reached very quickly, for the first for timesteps in the space of 0.1 s -0.3 s. The highest timestep will mimic the reality better, but it might be too demanding computationally.

<gallery mode="packed" heights=200 px>
Image:GE715L1.PNG|Figure 1: Energy at timestep 0.001 s
Image:GE715L2.PNG|Figure 2: Temperature at timestep 0.001 s
Image:GE715L3.PNG|Figure 3: Pressure at timestep 0.001 s
Image:GE715L4.PNG|Figure 4: Energy at timestep 0.0025 s
Image:GE715L5.PNG|Figure 5: Energy at timestep 0.0075 s
Image:GE715L6.PNG|Figure 6: Energy at timestep 0.01 s
Image:GE715L7.PNG|Figure 7: Energy at timestep 0.015 s
</gallery>

== Specific conditions ==

The next task was to run simulations varying temperatures at constant pressure. This shows the dependence of density on temperature at a constant pressure. Considering the ideal gas law the density will be inversely proportional with temperature.

If we run simulations under different conditions, we get the results shown in the gallery below. Bear in mind that some errors are almost negligible.

<gallery mode="packed" heights=400 px>
Image:GE715L8.PNG|Figure 8: Density vs temperature at p=2
Image:GE715L9.PNG|Figure 9: Density vs temperature at p=3
Image:GE715L16.PNG|Figure 10: Inverse of density vs temperature at p=2
Image:GE715L17.PNG|Figure 11: Inverse of density vs temperature at p=3
</gallery>

The inverses are of densities are shown to prove that this simulation obeys the proportionality proposed ideal gas law, where the density and temperature are inversely proportional.

== Radial distribution functions ==

The radial distribution functions from the calculations are given in the gallery below.

<gallery mode="packed" heights=400 px>
Image:GE715L10.PNG|Figure 12: RDF for gas
Image:GE715L11.PNG|Figure 13: RDF for liquid
Image:GE715L12.PNG|Figure 14: RDF for solid
Image:GE715L13.PNG|Figure 15: Integral of RDF for gas
Image:GE715L14.PNG|Figure 16: Integral of RDF for liquid
Image:GE715L15.PNG|Figure 17: Integral of RDF for solid
</gallery>

We can see from the graphs the differences between different phases and their order. It is clear that in a solid lattice structure is fairly rigid (peaks are the neighbouring atoms), liquid is slightly less ordered, while gas has almost. no real order at all. Lattice spacing in a solid is <math>\sigma</math>, while the first peak represents the nearest neighbouring atoms, while the next two represent the atom that are not nearest neighbours, but are less than two atoms away in a straight line from the central atom. The interactions are dropping of with distance, as expected in Lennard-Jones model.

== Conclusion ==

It has been determined that appropriate timestep for dynamics calculations is 0.0025 s. Simulations were performed to determine that temperature and density of a liquid are inversely proportional. It was also found that radial distribution function depends heavily on the phase of matter, with order being clearly seen in solid lattice.

Rep:Liquids GE715

2018-01-26T18:11:14Z

Ge715:

File:GE715L17.PNG

2018-01-26T18:11:04Z

Ge715:

File:GE715L16.PNG

2018-01-26T18:10:43Z

Ge715:

Rep:Liquids GE715

2018-01-26T17:42:34Z

Ge715:

== Abstract ==

In this experiment properties of liquids were modelled using molecular dynamics experiments. The effect of the temperature on the density as well as the radial distribution functions were explored.

== Introduction ==

Molecular dynamics is a modern computational technique concerned with calculating properties of substances by modelling them and their behaviour on atomic level with powerful computers.

Lennard Jones potential equals <math>\phi\left(r\right) = 4\epsilon \left( \frac{\sigma^{12}}{r^{12}} - \frac{\sigma^6}{r^6} \right)</math>. Potential energy is zero at the point when <math>r=\sigma</math>, thus <math>r_0=\sigma</math>.

Force is simply the derivative of the potential, thus <math>F = 4\epsilon \left( -12 \frac{\sigma^{12}}{r^{13}} + 6 \frac{\sigma^6}{r^7} \right)</math>. At <math>r=\sigma</math> this equals <math>F = - \frac{24\epsilon}{\sigma}</math>. Well depth is <math>\phi\left(r_{eq}\right) = - \epsilon</math>, where <math>r_{eq}=\sigma \sqrt[6]{2}</math>.

The following integrals are equal to: <math>\int_{2\sigma}^\infty \phi\left(r\right)\mathrm{d}r=-0.0248</math>, <math>\int_{2.5\sigma}^\infty \phi\left(r\right)\mathrm{d}r=-0.00818</math>, and <math>\int_{3\sigma}^\infty \phi\left(r\right)\mathrm{d}r=-0.00329</math>.

=== Water ===

There are roughly <math>3.3*10^{22}</math> molecules in 1 mL of water. However, 10000 atoms of water take up about <math>3*10^{-19} mL</math> or <math>300 nm^{3}</math>.

=== Atoms ===

Atom at position <math>(0.5,0.5,0.5)</math> moved for <math>(0.7,0.6,0.2)</math> would be found at <math>(0.2,0.1,0.7)</math>.

=== Reduced units ===

<math>LJ_{cut-off}=1.088 nm</math>

<math>T=180 K</math>

<math>r=0.382 nm</math>

<math>E=1 \frac{kJ}{mol}</math>

== Aims and Objectives ==

The aim of this experiment was to explore the modern use of molecular dynamics for simulating simple liquids.

The objectives were to look at how timesteps used in experiment affect the quality of the results, how temperature affects density and find radial distribution functions and their integrals for different phases of matter.

== Methods ==

The calculations were performed on a high performance computer using LAMMPS software package, while the calculations of radial distribution functions were done using VMD (delta r=0.05).

== Equilibration ==

If two atoms were generated very close together that would mean the force between them would be huge, which would have a potential to heavily disturb the simulation.

=== Commands ====

pair_style lj/cut 3.0 tells us where the interactions are cut off, as seen in the integral above, <math>3.0 \sigma</math> is a very good cut-off, since there is no significant perturbation.

=== Simulations ===

In the gallery below graphs are shown. It can be seen that as timesteps are increasing, the quality of results is decreasing. Timpestep of 0.015 s is wholy inappropriate, but I have chosen to run the experiments at the timestep of 0.0025 s. Results would probably also be almost as good using 0.01 s.

<gallery mode="packed" heights=200 px>
Image:GE715L1.PNG|Figure 1: Energy at timestep 0.001 s
Image:GE715L2.PNG|Figure 2: Temperature at timestep 0.001 s
Image:GE715L3.PNG|Figure 3: Pressure at timestep 0.001 s
Image:GE715L4.PNG|Figure 4: Energy at timestep 0.0025 s
Image:GE715L5.PNG|Figure 5: Energy at timestep 0.0075 s
Image:GE715L6.PNG|Figure 6: Energy at timestep 0.01 s
Image:GE715L7.PNG|Figure 7: Energy at timestep 0.015 s
</gallery>

== Specific conditions ==

If we run simulations under different conditions, we get the results shown in the gallery below. Bear in mind that some errors are almost negligible.

<gallery mode="packed" heights=400 px>
Image:GE715L8.PNG|Figure 8: Density vs temperature at p=2
Image:GE715L9.PNG|Figure 9: Density vs temperature at p=3
</gallery>

== Radial distribution functions ==

The radial distribution functions from the calculations are given in the gallery below.

<gallery mode="packed" heights=400 px>
Image:GE715L10.PNG|Figure 10: RDF for gas
Image:GE715L11.PNG|Figure 11: RDF for liquid
Image:GE715L12.PNG|Figure 12: RDF for solid
Image:GE715L13.PNG|Figure 13: Integral of RDF for gas
Image:GE715L14.PNG|Figure 14: Integral of RDF for liquid
Image:GE715L15.PNG|Figure 15: Integral of RDF for solid
</gallery>

We can see from the graphs the differences between different phases and their order. It is clear that in a solid lattice structure is fairly rigid (peaks are the neighbouring atoms), liquid is slightly less ordered, while gas has almost. no real order at all. Lattice spacing in a solid is <math>\sigma</math>, while the first peak represents the nearest neighbouring atoms, while the next two represent the atom that are not nearest neighbours, but are less than two atoms away in a straight line from the central atom.

Rep:Liquids GE715

2018-01-26T17:37:00Z

Ge715:

== Abstract ==

In this experiment properties of liquids were modelled using molecular dynamics experiments. The effect of the temperature on the density as well as the radial distribution functions were explored.

== Introduction ==

Molecular dynamics is a modern computational technique concerned with calculating properties of substances by modelling them and their behaviour on atomic level with powerful computers.

Lennard Jones potential equals <math>\phi\left(r\right) = 4\epsilon \left( \frac{\sigma^{12}}{r^{12}} - \frac{\sigma^6}{r^6} \right)</math>. Potential energy is zero at the point when <math>r=\sigma</math>, thus <math>r_0=\sigma</math>.

Force is simply the derivative of the potential, thus <math>F = 4\epsilon \left( -12 \frac{\sigma^{12}}{r^{13}} + 6 \frac{\sigma^6}{r^7} \right)</math>. At <math>r=\sigma</math> this equals <math>F = - \frac{24\epsilon}{\sigma}</math>. Well depth is <math>\phi\left(r_{eq}\right) = 1 \frac{kJ}{mol}</math>.

The following integrals are equal to: <math>\int_{2\sigma}^\infty \phi\left(r\right)\mathrm{d}r=-0.0248</math>, <math>\int_{2.5\sigma}^\infty \phi\left(r\right)\mathrm{d}r=-0.00818</math>, and <math>\int_{3\sigma}^\infty \phi\left(r\right)\mathrm{d}r=-0.00329</math>.

=== Water ===

There are roughly <math>3.3*10^{22}</math> molecules in 1 mL of water. However, 10000 atoms of water take up about <math>3*10^{-19} mL</math> or <math>300 nm^{3}</math>.

=== Atoms ===

Atom at position <math>(0.5,0.5,0.5)</math> moved for <math>(0.7,0.6,0.2)</math> would be found at <math>(0.2,0.1,0.7)</math>.

=== Reduced units ===

<math>LJ_{cut-off}=1.088 nm</math>

<math>T=180 K</math>

<math>r=0.382 nm</math>

<math>E=1 \frac{kJ}{mol}</math>

== Aims and Objectives ==

The aim of this experiment was to explore the modern use of molecular dynamics for simulating simple liquids.

The objectives were to look at how timesteps used in experiment affect the quality of the results, how temperature affects density and find radial distribution functions and their integrals for different phases of matter.

== Methods ==

The calculations were performed on a high performance computer using LAMMPS software package, while the calculations of radial distribution functions were done using VMD (delta r=0.05).

== Equilibration ==

If two atoms were generated very close together that would mean the force between them would be huge, which would have a potential to heavily disturb the simulation.

=== Commands ====

pair_style lj/cut 3.0 tells us where the interactions are cut off, as seen in the integral above, <math>3.0 \sigma</math> is a very good cut-off, since there is no significant perturbation.

=== Simulations ===

In the gallery below graphs are shown. It can be seen that as timesteps are increasing, the quality of results is decreasing. Timpestep of 0.015 s is wholy inappropriate, but I have chosen to run the experiments at the timestep of 0.0025 s. Results would probably also be almost as good using 0.01 s.

<gallery mode="packed" heights=200 px>
Image:GE715L1.PNG|Figure 1: Energy at timestep 0.001 s
Image:GE715L2.PNG|Figure 2: Temperature at timestep 0.001 s
Image:GE715L3.PNG|Figure 3: Pressure at timestep 0.001 s
Image:GE715L4.PNG|Figure 4: Energy at timestep 0.0025 s
Image:GE715L5.PNG|Figure 5: Energy at timestep 0.0075 s
Image:GE715L6.PNG|Figure 6: Energy at timestep 0.01 s
Image:GE715L7.PNG|Figure 7: Energy at timestep 0.015 s
</gallery>

== Specific conditions ==

If we run simulations under different conditions, we get the results shown in the gallery below. Bear in mind that some errors are almost negligible.

<gallery mode="packed" heights=400 px>
Image:GE715L8.PNG|Figure 8: Density vs temperature at p=2
Image:GE715L9.PNG|Figure 9: Density vs temperature at p=3
</gallery>

== Radial distribution functions ==

The radial distribution functions from the calculations are given in the gallery below.

<gallery mode="packed" heights=400 px>
Image:GE715L10.PNG|Figure 10: RDF for gas
Image:GE715L11.PNG|Figure 11: RDF for liquid
Image:GE715L12.PNG|Figure 12: RDF for solid
Image:GE715L13.PNG|Figure 13: Integral of RDF for gas
Image:GE715L14.PNG|Figure 14: Integral of RDF for liquid
Image:GE715L15.PNG|Figure 15: Integral of RDF for solid
</gallery>

We can see from the graphs the differences between different phases and their order. It is clear that in a solid lattice structure is fairly rigid (peaks are the neighbouring atoms), liquid is slightly less ordered, while gas has almost. no real order at all. Lattice spacing in a solid is <math>\sigma</math>, while the first peak represents the nearest neighbouring atoms, while the next two represent the atom that are not nearest neighbours, but are less than two atoms away in a straight line from the central atom.

Rep:Liquids GE715

2018-01-26T17:19:21Z

Ge715:

== Introduction ==

Lennard Jones potential equals <math>\phi\left(r\right) = 4\epsilon \left( \frac{\sigma^{12}}{r^{12}} - \frac{\sigma^6}{r^6} \right)</math>. Potential energy is zero at the point when <math>r=\sigma</math>, thus <math>r_0=\sigma</math>.

Force is simply the derivative of the potential, thus <math>F = 4\epsilon \left( -12 \frac{\sigma^{12}}{r^{13}} + 6 \frac{\sigma^6}{r^7} \right)</math>. At <math>r=\sigma</math> this equals <math>F = - \frac{24\epsilon}{\sigma}</math>. Well depth is <math>\phi\left(r_{eq}\right) = 1 \frac{kJ}{mol}</math>.

The following integrals are equal to: <math>\int_{2\sigma}^\infty \phi\left(r\right)\mathrm{d}r=-0.0248</math>, <math>\int_{2.5\sigma}^\infty \phi\left(r\right)\mathrm{d}r=-0.00818</math>, and <math>\int_{3\sigma}^\infty \phi\left(r\right)\mathrm{d}r=-0.00329</math>.

=== Water ===

There are roughly <math>3.3*10^{22}</math> molecules in 1 mL of water. However, 10000 atoms of water take up about <math>3*10^{-19} mL</math> or <math>300 nm^{3}</math>.

=== Atoms ===

Atom at position <math>(0.5,0.5,0.5)</math> moved for <math>(0.7,0.6,0.2)</math> would be found at <math>(0.2,0.1,0.7)</math>.

=== Reduced units ===

<math>LJ cut-off=1.088 nm</math>
<math>T=180 K</math>
<math>r=0.382 nm</math>
<math>E=1 \frac{kJ}{mol}</math>

== Equilibration ==

If two atoms were generated very close together that would mean the force between them would be huge, which would have a potential to heavily disturb the simulation.

=== Commands ====

pair_style lj/cut 3.0 tells us where the interactions are cut off, as seen in the integral above, <math>3.0 \sigma</math> is a very good cut-off, since there is no significant perturbation.

=== Simulations ===

In the gallery below graphs are shown. It can be seen that as timesteps are increasing, the quality of results is decreasing. Timpestep of 0.015 s is wholy inappropriate, but I have chosen to run the experiments at the timestep of 0.0025 s. Results would probably also be almost as good using 0.01 s.

<gallery mode="packed" heights=200 px>
Image:GE715L1.PNG|Figure 1: Energy at timestep 0.001 s
Image:GE715L2.PNG|Figure 2: Temperature at timestep 0.001 s
Image:GE715L3.PNG|Figure 3: Pressure at timestep 0.001 s
Image:GE715L4.PNG|Figure 4: Energy at timestep 0.0025 s
Image:GE715L5.PNG|Figure 5: Energy at timestep 0.0075 s
Image:GE715L6.PNG|Figure 6: Energy at timestep 0.01 s
Image:GE715L7.PNG|Figure 7: Energy at timestep 0.015 s
</gallery>

== Specific conditions ==

If we run simulations under different conditions, we get the results shown in the gallery below. Bear in mind that some errors are almost negligible.

<gallery mode="packed" heights=400 px>
Image:GE715L8.PNG|Figure 8: Density vs temperature at p=2
Image:GE715L9.PNG|Figure 9: Density vs temperature at p=3
</gallery>

== Radial distribution functions ==

The radial distribution functions from the calculations are given in the gallery below.

<gallery mode="packed" heights=400 px>
Image:GE715L10.PNG|Figure 10: RDF for gas
Image:GE715L11.PNG|Figure 11: RDF for liquid
Image:GE715L12.PNG|Figure 12: RDF for solid
Image:GE715L13.PNG|Figure 13: Integral of RDF for gas
Image:GE715L14.PNG|Figure 14: Integral of RDF for liquid
Image:GE715L15.PNG|Figure 15: Integral of RDF for solid
</gallery>

We can see from the graphs the differences between different phases and their order. It is clear that in a solid lattice structure is fairly rigid (peaks are the neighbouring atoms), liquid is slightly less ordered, while gas has almost. no real order at all. Lattice spacing in a solid is <math>\sigma</math>, while the first peak represents the nearest neighbouring atoms, while the next two represent the atom that are not nearest neighbours, but are less than two atoms away in a straight line from the central atom.

File:GE715L15.PNG

2018-01-26T17:12:45Z

Ge715:

File:GE715L14.PNG

2018-01-26T17:12:30Z

Ge715:

File:GE715L13.PNG

2018-01-26T17:12:00Z

Ge715:

Rep:Liquids GE715

2018-01-26T17:11:22Z

Ge715:

File:GE715L12.PNG

2018-01-26T17:06:35Z

Ge715:

Rep:Liquids GE715

2018-01-26T16:58:16Z

Ge715:

File:GE715L11.PNG

2018-01-26T16:57:50Z

Ge715:

Rep:Liquids GE715

2018-01-26T16:55:56Z

Ge715:

File:GE715L10.PNG

2018-01-26T16:55:43Z

Ge715:

File:GE715L9.PNG

2018-01-26T16:51:07Z

Ge715:

Rep:Liquids GE715

2018-01-26T16:45:00Z

Ge715:

Rep:Liquids GE715

2018-01-26T16:44:07Z

Ge715:

== Introduction ==

Lennard Jones potential equals <math>\phi\left(r\right) = 4\epsilon \left( \frac{\sigma^{12}}{r^{12}} - \frac{\sigma^6}{r^6} \right)</math>. Potential energy is zero at the point when <math>r=\sigma</math>, thus <math>r_0=\sigma</math>.

Force is simply the derivative of the potential, thus <math>F = 4\epsilon \left( -12 \frac{\sigma^{12}}{r^{13}} + 6 \frac{\sigma^6}{r^7} \right)</math>. At <math>r=\sigma</math> this equals <math>F = - \frac{24\epsilon}{\sigma}</math>. Well depth is <math>\phi\left(r_{eq}\right) = 1 \frac{kJ}{mol}</math>.

The following integrals are equal to: <math>\int_{2\sigma}^\infty \phi\left(r\right)\mathrm{d}r=-0.0248</math>, <math>\int_{2.5\sigma}^\infty \phi\left(r\right)\mathrm{d}r=-0.00818</math>, and <math>\int_{3\sigma}^\infty \phi\left(r\right)\mathrm{d}r=-0.00329</math>.

=== Water ===

There are roughly <math>3.3*10^{22}</math> molecules in 1 mL of water. However, 10000 atoms of water take up about <math>3*10^{-19} mL</math> or <math>300 nm^{3}</math>.

=== Atoms ===

Atom at position <math>(0.5,0.5,0.5)</math> moved for <math>(0.7,0.6,0.2)</math> would be found at <math>(0.2,0.1,0.7)</math>.

=== Reduced units ===

<math>LJ cut-off=1.088 nm</math>
<math>T=180 K</math>
<math>r=0.382 nm</math>
<math>E=1 \frac{kJ}{mol}</math>

== Equilibration ==

If two atoms were generated very close together that would mean the force between them would be huge, which would have a potential to heavily disturb the simulation.

=== Commands ====

pair_style lj/cut 3.0 tells us where the interactions are cut off, as seen in the integral above, <math>3.0 \sigma</math> is a very good cut-off, since there is no significant perturbation.

=== Simulations ===

In the gallery below graphs are shown. It can be seen that as timesteps are increasing, the quality of results is decreasing. Timpestep of 0.015 s is wholy inappropriate, but I have chosen to run the experiments at the timestep of 0.0025 s. Results would probably also be almost as good using 0.01 s.

<gallery mode="packed" heights=200 px>
Image:GE715L1.PNG|Figure 1: Energy at timestep 0.001 s
Image:GE715L2.PNG|Figure 2: Temperature at timestep 0.001 s
Image:GE715L3.PNG|Figure 3: Pressure at timestep 0.001 s
Image:GE715L4.PNG|Figure 4: Energy at timestep 0.0025 s
Image:GE715L5.PNG|Figure 5: Energy at timestep 0.0075 s
Image:GE715L6.PNG|Figure 6: Energy at timestep 0.01 s
Image:GE715L7.PNG|Figure 7: Energy at timestep 0.015 s
</gallery>

== Specific conditions ==

If we run simulations under different conditions, we get the results shown in the gallery below.

<gallery mode="packed" heights=800 px>
Image:GE715L8.PNG|Figure 8: Density vs temperature at p=2
Image:GE715L9.PNG|Figure 9: Density vs temperature at p=3
</gallery>

Rep:Liquids GE715

2018-01-26T16:43:52Z

Ge715:

== Introduction ==

Lennard Jones potential equals <math>\phi\left(r\right) = 4\epsilon \left( \frac{\sigma^{12}}{r^{12}} - \frac{\sigma^6}{r^6} \right)</math>. Potential energy is zero at the point when <math>r=\sigma</math>, thus <math>r_0=\sigma</math>.

Force is simply the derivative of the potential, thus <math>F = 4\epsilon \left( -12 \frac{\sigma^{12}}{r^{13}} + 6 \frac{\sigma^6}{r^7} \right)</math>. At <math>r=\sigma</math> this equals <math>F = - \frac{24\epsilon}{\sigma}</math>. Well depth is <math>\phi\left(r_{eq}\right) = 1 \frac{kJ}{mol}</math>.

The following integrals are equal to: <math>\int_{2\sigma}^\infty \phi\left(r\right)\mathrm{d}r=-0.0248</math>, <math>\int_{2.5\sigma}^\infty \phi\left(r\right)\mathrm{d}r=-0.00818</math>, and <math>\int_{3\sigma}^\infty \phi\left(r\right)\mathrm{d}r=-0.00329</math>.

=== Water ===

There are roughly <math>3.3*10^{22}</math> molecules in 1 mL of water. However, 10000 atoms of water take up about <math>3*10^{-19} mL</math> or <math>300 nm^{3}</math>.

=== Atoms ===

Atom at position <math>(0.5,0.5,0.5)</math> moved for <math>(0.7,0.6,0.2)</math> would be found at <math>(0.2,0.1,0.7)</math>.

=== Reduced units ===

<math>LJ cut-off=1.088 nm</math>
<math>T=180 K</math>
<math>r=0.382 nm</math>
<math>E=1 \frac{kJ}{mol}</math>

== Equilibration ==

If two atoms were generated very close together that would mean the force between them would be huge, which would have a potential to heavily disturb the simulation.

=== Commands ====

pair_style lj/cut 3.0 tells us where the interactions are cut off, as seen in the integral above, <math>3.0 \sigma</math> is a very good cut-off, since there is no significant perturbation.

=== Simulations ===

In the gallery below graphs are shown. It can be seen that as timesteps are increasing, the quality of results is decreasing. Timpestep of 0.015 s is wholy inappropriate, but I have chosen to run the experiments at the timestep of 0.0025 s. Results would probably also be almost as good using 0.01 s.

<gallery mode="packed" heights=200 px>
Image:GE715L1.PNG|Figure 1: Energy at timestep 0.001 s
Image:GE715L2.PNG|Figure 2: Temperature at timestep 0.001 s
Image:GE715L3.PNG|Figure 3: Pressure at timestep 0.001 s
Image:GE715L4.PNG|Figure 4: Energy at timestep 0.0025 s
Image:GE715L5.PNG|Figure 5: Energy at timestep 0.0075 s
Image:GE715L6.PNG|Figure 6: Energy at timestep 0.01 s
Image:GE715L7.PNG|Figure 7: Energy at timestep 0.015 s
</gallery>

== Specific conditions ==

If we run simulations under different conditions, we get the results shown in the gallery below.

<gallery mode="packed" heights=200 px>
Image:GE715L8.PNG|Figure 8: Density vs temperature at p=2
Image:GE715L9.PNG|Figure 9: Density vs temperature at p=3
</gallery>

File:GE715L8.PNG

2018-01-26T16:41:54Z

Ge715:

File:GE715L7.PNG

2018-01-26T16:27:01Z

Ge715:

File:GE715L6.PNG

2018-01-26T16:26:46Z

Ge715:

File:GE715L5.PNG

2018-01-26T16:26:32Z

Ge715:

File:GE715L4.PNG

2018-01-26T16:26:19Z

Ge715:

File:GE715L3.PNG

2018-01-26T16:25:59Z

Ge715:

File:GE715L2.PNG

2018-01-26T16:25:45Z

Ge715:

File:GE715L1.PNG

2018-01-26T16:25:29Z

Ge715:

Rep:Liquids GE715

2018-01-26T16:24:56Z

Ge715: Created page with "== Introduction == Lennard Jones potential equals <math>\phi\left(r\right) = 4\epsilon \left( \frac{\sigma^{12}}{r^{12}} - \frac{\sigma^6}{r^6} \right)</math>. Potential ener..."

Rep:CMPY3 GE715

2017-11-22T11:57:41Z

Ge715: /* Monte Carlo simulation */ Note

== Ising Model ==

Ising model is built on the assumption that a lattice that can be of any dimensionality contains spins which can either be up (+1 or down (-1) and interactions can be studied using this model. The model is extremely simplified, since in real systems the spins can point in any direction, but the basic properties of magnetic systems, such as magnetisation, energy levels and Curie temperature, may be studied using this simplification.

The interaction energy in Ising model is defined as <math>E=- \frac{1}{2} J \sum_i^N \sum_{j\ \in\ \mathrm{neighbours}\left(i\right)} s_i s_j</math>. The energy will be the lowest possible when all the spins point in the same direction. In that case <math>s_i s_j=1</math> for all possible <math>i,j</math> pairs. The number of pairs for each <math>i</math> depends on the number of dimensions.

<math>P=2D</math>, where <math>P</math> represents the number of pairs and <math>D</math> the number of dimensions. The lowest possible energy will thus be <math>E=- \frac{1}{2} JNP</math>, where <math>N</math> represents the total number of spins. Thus the lowest possible energy may also be simplified as <math>E=-DNJ</math>, as requested.

The multiplicity of this particular state is 2 (all spins up or all spins down), thus <math>\Omega=2</math>, entropy of the system is <math>S=k_b ln \Omega=9.6*10^{-24} JK^{-1}</math>.

=== Flip of a spin ===

If one of the spins flips, the energy increases. In the given case (<math>D=3, N=1000</math>), the <math>E_i=-3000 J</math>, while the flip of one spin increases the energy by <math>12 J</math> (6 pairs, <math>2 J</math> per pair), so that <math>E_f=-2988 J</math>. Note that <math>J</math> does not represent joules in this case.

Multiplicity has increased to <math>\Omega=2000</math>, thus <math>S=1.05*10^{-22} JK^{-1}</math>.

<math>\Delta S=1.04*10^{-22} JK^{-1}</math>.

This means that entropy has increased by more than 100 times.

=== Magnetisation in the sketch ===

The 1D system shown has the magnetisation of <math>M=1</math>. This is also true for the 2D system shown. For the model system with <math>N=1000, D=3</math> the expected magnetisation at absolute zero would be <math>M=\pm 1000</math>.

== Code for magnetisation and energy ==

Python code was written to calculate the energy and magnetisation of a configuration and tested a few times. The test results are shown in figure 1 below.

<gallery mode="packed" heights=800 px>
Image:GE715_CMP_0.png|Figure 1: Test results
</gallery>

== Monte Carlo simulation ==

If no weighing is used, the calculation would be far too demanding even for current supercomputers, because 100 spins can take about <math>1.27*10^{30}</math> arrangements, which means calculating the average for such system would take more than <math>1.3*10^{21} s</math>, which is 40 trillion years, or about 3 times the age of the Universe, thus a better system must be used to do this quicker. One option is using Monte Carlo simulation, which mimics real transitions. It is essentially optimisation algorithm that is looking for stable state, but with a possibility of going into higher energy state should it have enough energy.

Below Curie temperature it is to be expected, that absolute value of magnetisation will be close to LxL, where L is the size of the lattice (works also in more than 2 dimensions).

The old code was accidentally deleted before it could be added to the report.

== Accelerating the code ==

The code at the beginning was written using a loop for magnetisation and a double loop for energy. 2000 Monte Carlo steps needed was <math>M=9.32 s \pm 0.10 s</math>. After optimizing magnetisation only this has dropped to <math>M=8.88 s \pm 0.04 s</math>, while after optimizing both magnetisation and energy with using Numpy's sum, roll and multiply functions, the time was further reduced to <math>M=0.0796 s \pm 0.018 s</math>. It was later realised that only one roll in each direction is needed to produce the correct energy, so the final time for 2000 steps was successfully reduced to <math>M=0.418 s \pm 0.005 s</math>.

The whole of the code used is given below.

<pre>
import numpy as np

class IsingLattice:

E = 0.0
E2 = 0.0
M = 0.0
M2 = 0.0
energy_x=0
energy2_x=0
magnetisation_x=0
magnetisation2_x=0

n_cycles = 0

def __init__(self, n_rows, n_cols):
self.n_rows = n_rows
self.n_cols = n_cols
self.lattice = np.random.choice([-1,1], size=(n_rows, n_cols))

def energy(self):
a=np.roll(self.lattice,1,axis=0)
a1=np.multiply(a,self.lattice)
x1=np.sum(a1)
c=np.roll(self.lattice,1,axis=1)
c1=np.multiply(c,self.lattice)
x2=np.sum(c1)
x=x1+x2
energy = -x
return energy

def magnetisation(self):
"Return the total magnetisation of the current lattice configuration."
return np.sum(self.lattice)

def energy_flip(self,random_i,random_j):
"Return the total energy of the current lattice after flip configuration."
self.lattice[random_i][random_j]=-self.lattice[random_i][random_j]
a=np.roll(self.lattice,1,axis=0)
a1=np.multiply(a,self.lattice)
x1=np.sum(a1)
c=np.roll(self.lattice,1,axis=1)
c1=np.multiply(c,self.lattice)
x2=np.sum(c1)
x=x1+x2
energy = -x
return energy

def montecarlostep(self, T):
# complete this function so that it performs a single Monte Carlo step
energy = self.energy()
#the following two lines will select the coordinates of the random spin for you
random_i = np.random.choice(range(0, self.n_rows))
random_j = np.random.choice(range(0, self.n_cols))
#the following line will choose a random number in the range [0,1) for you
random_number = np.random.random()
self.lattice[random_i][random_j] = -self.lattice[random_i][random_j]
energy2 = self.energy()
dE = energy2 - energy

if energy2<energy or random_number<=np.e**(-(energy2-energy)/T):
yyy=True
else:
self.lattice[random_i][random_j] = -self.lattice[random_i][random_j]

self.n_cycles+=1

self.E += self.energy()
self.E2 += self.energy()**2
self.M += self.magnetisation()
self.M2 += self.magnetisation()**2
return self.energy(), self.magnetisation()

def statistics(self):
# complete this function so that it calculates the correct values for the averages of E, E*E (E2), M, M*M (M2), and returns them
avE = self.E / self.n_cycles
avE2 = self.E2 / self.n_cycles
avM = self.M / self.n_cycles
avM2 = self.M / self.n_cycles
return avE, avE2, avM, avM2, self.n_cycles
</pre>

== Temperature effects ==

=== Correction of averaging code ===

A few experiments were performed at different temperatures using ILfinalframe.py, to see how long it usually takes for the system to equilibrate. At low temperatures it usually happened before 1000 steps were performed, thus 1000 Monte Carlo steps was decided to be appropriate cut-off period for calculation of average values of magnetisation and energy. The figures 2-4, showing the experiments performed are shown in the gallery below.

<gallery mode="packed" heights=250 px>
Image:GE715_CMP_1.png|Figure 2: Test result @ 0.5 K
Image:GE715_CMP_2.png|Figure 3: Test result @ 1 K
Image:GE715_CMP_3.png|Figure 4: Test result @ 1.5 K
</gallery>

The definitions of two class functions were changed as shown below to accommodate this change.

<pre>
def montecarlostep(self, T):
# complete this function so that it performs a single Monte Carlo step
energy = self.energy()
#the following two lines will select the coordinates of the random spin for you
random_i = np.random.choice(range(0, self.n_rows))
random_j = np.random.choice(range(0, self.n_cols))
#the following line will choose a random number in the range [0,1) for you
random_number = np.random.random()
self.lattice[random_i][random_j] = -self.lattice[random_i][random_j]
energy2 = self.energy()
dE = energy2 - energy

if energy2<energy or random_number<=np.e**(-(energy2-energy)/T):
yyy=True
else:
self.lattice[random_i][random_j] = -self.lattice[random_i][random_j]

self.n_cycles+=1
if self.n_cycles==1000:
self.E=0
self.E2=0
self.M=0
self.M2=0
self.E += self.energy()
self.E2 += self.energy()**2
self.M += self.magnetisation()
self.M2 += self.magnetisation()**2
return self.energy(), self.magnetisation()

def statistics(self):
# complete this function so that it calculates the correct values for the averages of E, E*E (E2), M, M*M (M2), and returns them
avE = self.E / self.n_cycles
avE2 = self.E2 / self.n_cycles
avM = self.M / self.n_cycles
avM2 = self.M / self.n_cycles
if self.n_cycles>1000:
avE = self.E / (self.n_cycles-1000)
avE2 = self.E2 / (self.n_cycles-1000)
avM = self.M / (self.n_cycles-1000)
avM2 = self.M / (self.n_cycles-1000)
return avE, avE2, avM, avM2, self.n_cycles
</pre>

=== Range of temperatures ===

The experiment was run for 8x8 lattice at a temperature range of 0.25 K-5.00 K, with a point every 0.05 K. The experiments using a different number of total steps were performed and the comparison is shown below in figures 5-7.

<gallery mode="packed" heights=250 px>
Image:GE715_CMP_4.png|Figure 5: Test result for 8x8 lattice with 11000 steps
Image:GE715_CMP_5.png|Figure 6: Test result for 8x8 lattice with 21000 steps
Image:GE715_CMP_6.png|Figure 7: Test result for 8x8 lattice with 101000 steps
</gallery>

There is a small, but noticeable improvement in the smoothness of the curves obtained for energy with using more steps, thus the highest number of steps (101000 steps) was used in the further experiments regardless of the higher need for computational power.

== Lattice size effects ==

With increasing lattice size the energy and magnetisation curves become more realistic, since the long range effects are better accounted for. 7 different lattice sizes were tested and energy and magnetisation for each are shown below in figure 8-14 and 15-21 respectively.

<gallery mode="packed" heights=250 px>
Image:GE715_CMP_7.png|Figure 8: Energy result for 2x2 lattice
Image:GE715_CMP_8.png|Figure 9: Energy result for 4x4 lattice
Image:GE715_CMP_9.png|Figure 10: Energy result for 8x8 lattice
Image:GE715_CMP_10.png|Figure 11: Energy result for 16x16 lattice
Image:GE715_CMP_11.png|Figure 12: Energy result for 32x32 lattice
Image:GE715_CMP_12.png|Figure 13: Energy result for 64x64 lattice
Image:GE715_CMP_13.png|Figure 14: Energy result for 128x128 lattice
Image:GE715_CMP_14.png|Figure 15: Magnetisation result for 2x2 lattice
Image:GE715_CMP_15.png|Figure 16: Magnetisation result for 4x4 lattice
Image:GE715_CMP_16.png|Figure 17: Magnetisation result for 8x8 lattice
Image:GE715_CMP_17.png|Figure 18: Magnetisation result for 16x16 lattice
Image:GE715_CMP_18.png|Figure 19: Magnetisation result for 32x32 lattice
Image:GE715_CMP_19.png|Figure 20: Magnetisation result for 64x64 lattice
Image:GE715_CMP_20.png|Figure 21: Magnetisation result for 128x128 lattice
</gallery>

Overall comparison of energy and magnetisation per spin in each lattice was made and can be seen in figures 22 and 23 below.

<gallery mode="packed" heights=500 px>
Image:GE715_CMP_21.png|Figure 22: Energy comparison
Image:GE715_CMP_22.png|Figure 23: Magnetisation comparison
</gallery>

The code used for producing figure 23 above is shown below. Other codes are analagous. Example for defining ls2 is also shown Others are similar..

<pre>
ls2=loadtxt("2x2(100).dat")

figsize(12,8)
plot(ls2[:,0],ls2[:,3]/4,color="red",label="2x2")
plot(ls4[:,0],ls4[:,3]/16,color="orange",label="4x4")
plot(ls8[:,0],ls8[:,3]/64,color="yellow",label="8x8")
plot(ls16[:,0],ls16[:,3]/256,color="lime",label="16x16")
plot(ls32[:,0],ls32[:,3]/1024,color="blue",label="32x32")
plot(ls64[:,0],ls64[:,3]/4096,color="indigo",label="64x64")
plot(ls128[:,0],ls128[:,3]/16384,color="violet",label="128x128")
title("Magnetisation per spin")
xlabel("Temperature (K)")
ylabel("Magnetisation per spin")
legend()
savefig("Magnetisation.png")
show()
</pre>

Significant changes may be observed from 2x2 lattice up to 16x16 lattice, later the additional long range effects become less noticeable, but new effects are still present. It seems that most of the long range effects are captured by 16x16 lattice, but higher dimensions bring additional improvement.

== Heat Capacity ==

This formula holds true by definition: <math>C = \frac{\partial \left\langle E\right\rangle}{\partial T}</math>.

Variance of energy can be calculated in the following way: <math>\mathrm{Var}[X] = \left\langle E^2\right\rangle - \left\langle E\right\rangle^2</math>.

Since expected energy can also be expressed as partition function: <math>\left\langle E\right\rangle = \frac{\partial}{\partial \beta}lnZ</math>, variance may be expressed as <math>\mathrm{Var}[X] = \frac{\partial ^2}{\partial \beta ^2}lnZ = - \frac{\partial \left\langle E\right\rangle}{\partial \beta}</math>.

We can rearrange the equations given, because <math>T = \frac{1}{k_BT}</math>, so that we obtain <math>C = \frac{\mathrm{Var}[E]}{k_B T^2}</math>.

=== Heat capacity of the model ===

Heat capacities of the models were calculated using the script given below (example for 2x2 lattice).

<pre>
figsize(6,4)
plot(ls2[:,0],(ls2[:,2]-ls2[:,1]**2)/ls2[:,0],color="red",label="2x2")
xlabel("Temperature (K)")
ylabel("C")
savefig("C2.png")
show()
</pre>

Graphs of heat capacities are shown in the gallery below (figures 24-30). It can be seen that the model breaks down at low temperatures for bigger lattices, but it still gives a peak heat capacity.

<gallery mode="packed" heights=250 px>
Image:GE715_CMP_23.png|Figure 24: C result for 2x2 lattice
Image:GE715_CMP_24.png|Figure 25: C result for 4x4 lattice
Image:GE715_CMP_25.png|Figure 26: C result for 8x8 lattice
Image:GE715_CMP_26.png|Figure 27: C result for 16x16 lattice
Image:GE715_CMP_27.png|Figure 28: C result for 32x32 lattice
Image:GE715_CMP_28.png|Figure 29: C result for 64x64 lattice
Image:GE715_CMP_29.png|Figure 30: C result for 128x128 lattice
</gallery>

Comparison of heat capacity per spin is given in the figure 31 below.

<gallery mode="packed" heights=500 px>
Image:GE715_CMP_30.png|Figure 31: Heat capacity comparison
</gallery>

== Curie temperature ==

Curie temperature can be deduced from where there is a peak in heat capacity per spin. The graphs of heat capacity vs temperature is shown below (figures 32 and 33) for two different lattice sizes, once using Python script and once using C++ data provided.

<gallery mode="packed" heights=500 px>
Image:GE715_CMP_31.png|Figure 32: 16x16 lattice
Image:GE715_CMP_32.png|Figure 33: 64x64 lattice
</gallery>

=== Polynomial fit ===

The C++ data was also fitted with a polynomial, firstly over the whole range and afterwards just in the peak range. Fitting a polynomial proved to be hard as expected (since data is close to delta-like function), but it became more feasible when smaller range was used. Examples of both fits are shown in figures 34 and 35 below.

<gallery mode="packed" heights=500 px>
Image:GE715_CMP_33.png|Figure 34: Whole region fit
Image:GE715_CMP_34.png|Figure 35: Peak region fit
</gallery>

=== Finding Curie temperature ===

The last step was to find peak heat capacity for different lattice sizes in C++ data. The results are shown in the table below.

{| class="wikitable"
|-
! Lattice size
! Curie temperature
|-
|2x2
|2.48 K
|-
|4x4
|2.40 K
|-
|8x8
|2.32 K
|-
|16x16
|2.25 K
|-
|32x32
|2.30 K
|-
|64x64
|2.25 K
|}

The curie temperature for infinite lattice is known to be calculated from that of a finite lattice using equation <math>T_{C, L} = \frac{A}{L} + T_{C,\infty}</math>, where <math>T_{C, L}</math> is the Curie temperature of an <math>L\times L</math>lattice, <math>T_{C,\infty}</math> is the Curie temperature of an infinite lattice, and <math>A</math> is a constant. This was used to draw a graph shown in figure 36 below.

<gallery mode="packed" heights=500 px>
Image:GE715_CMP_35.png|Figure 36: Curie temperature for infinite lattice
</gallery>

The Curie temperature for infinite lattice was found to be <math>T=2.26 K</math>.

The code used for drawing the graph in figure 36 is shown below.

<pre>
L=array([1/2,1/4,1/8,1/16,1/32,1/64])
T=array([2.48,2.4,2.32,2.25,2.3,2.25])
X=array([1/2,0])
figsize(9,6)
plot(L,T,marker="D",linestyle=" ")
xlabel("1/L")
ylabel("T")
a=polyfit(L,T,1)
ylim(2.23,2.5)
plot(X,a[0]*X+a[1],color="lime")
title("Curie temperature")
savefig("Curie.png")
show()
</pre>

Rep:CMPY3 GE715

2017-11-22T11:54:20Z

Ge715: Monte Carlo

== Ising Model ==

Ising model is built on the assumption that a lattice that can be of any dimensionality contains spins which can either be up (+1 or down (-1) and interactions can be studied using this model. The model is extremely simplified, since in real systems the spins can point in any direction, but the basic properties of magnetic systems, such as magnetisation, energy levels and Curie temperature, may be studied using this simplification.

The interaction energy in Ising model is defined as <math>E=- \frac{1}{2} J \sum_i^N \sum_{j\ \in\ \mathrm{neighbours}\left(i\right)} s_i s_j</math>. The energy will be the lowest possible when all the spins point in the same direction. In that case <math>s_i s_j=1</math> for all possible <math>i,j</math> pairs. The number of pairs for each <math>i</math> depends on the number of dimensions.

<math>P=2D</math>, where <math>P</math> represents the number of pairs and <math>D</math> the number of dimensions. The lowest possible energy will thus be <math>E=- \frac{1}{2} JNP</math>, where <math>N</math> represents the total number of spins. Thus the lowest possible energy may also be simplified as <math>E=-DNJ</math>, as requested.

The multiplicity of this particular state is 2 (all spins up or all spins down), thus <math>\Omega=2</math>, entropy of the system is <math>S=k_b ln \Omega=9.6*10^{-24} JK^{-1}</math>.

=== Flip of a spin ===

If one of the spins flips, the energy increases. In the given case (<math>D=3, N=1000</math>), the <math>E_i=-3000 J</math>, while the flip of one spin increases the energy by <math>12 J</math> (6 pairs, <math>2 J</math> per pair), so that <math>E_f=-2988 J</math>. Note that <math>J</math> does not represent joules in this case.

Multiplicity has increased to <math>\Omega=2000</math>, thus <math>S=1.05*10^{-22} JK^{-1}</math>.

<math>\Delta S=1.04*10^{-22} JK^{-1}</math>.

This means that entropy has increased by more than 100 times.

=== Magnetisation in the sketch ===

The 1D system shown has the magnetisation of <math>M=1</math>. This is also true for the 2D system shown. For the model system with <math>N=1000, D=3</math> the expected magnetisation at absolute zero would be <math>M=\pm 1000</math>.

== Code for magnetisation and energy ==

Python code was written to calculate the energy and magnetisation of a configuration and tested a few times. The test results are shown in figure 1 below.

<gallery mode="packed" heights=800 px>
Image:GE715_CMP_0.png|Figure 1: Test results
</gallery>

== Monte Carlo simulation ==

If no weighing is used, the calculation would be far too demanding even for current supercomputers, because 100 spins can take about <math>1.27*10^{30}</math> arrangements, which means calculating the average for such system would take more than <math>1.3*10^{21} s</math>, which is 40 trillion years, or about 3 times the age of the Universe, thus a better system must be used to do this quicker. One option is using Monte Carlo simulation, which mimics real transitions. It is essentially optimisation algorithm that is looking for stable state, but with a possibility of going into higher energy state should it have enough energy.

Below Curie temperature it is to be expected, that absolute value of magnetisation will be close to LxL, where L is the size of the lattice (works also in more than 2 dimensions).

== Accelerating the code ==

The code at the beginning was written using a loop for magnetisation and a double loop for energy. 2000 Monte Carlo steps needed was <math>M=9.32 s \pm 0.10 s</math>. After optimizing magnetisation only this has dropped to <math>M=8.88 s \pm 0.04 s</math>, while after optimizing both magnetisation and energy with using Numpy's sum, roll and multiply functions, the time was further reduced to <math>M=0.0796 s \pm 0.018 s</math>. It was later realised that only one roll in each direction is needed to produce the correct energy, so the final time for 2000 steps was successfully reduced to <math>M=0.418 s \pm 0.005 s</math>.

The whole of the code used is given below.

<pre>
import numpy as np

class IsingLattice:

E = 0.0
E2 = 0.0
M = 0.0
M2 = 0.0
energy_x=0
energy2_x=0
magnetisation_x=0
magnetisation2_x=0

n_cycles = 0

def __init__(self, n_rows, n_cols):
self.n_rows = n_rows
self.n_cols = n_cols
self.lattice = np.random.choice([-1,1], size=(n_rows, n_cols))

def energy(self):
a=np.roll(self.lattice,1,axis=0)
a1=np.multiply(a,self.lattice)
x1=np.sum(a1)
c=np.roll(self.lattice,1,axis=1)
c1=np.multiply(c,self.lattice)
x2=np.sum(c1)
x=x1+x2
energy = -x
return energy

def magnetisation(self):
"Return the total magnetisation of the current lattice configuration."
return np.sum(self.lattice)

def energy_flip(self,random_i,random_j):
"Return the total energy of the current lattice after flip configuration."
self.lattice[random_i][random_j]=-self.lattice[random_i][random_j]
a=np.roll(self.lattice,1,axis=0)
a1=np.multiply(a,self.lattice)
x1=np.sum(a1)
c=np.roll(self.lattice,1,axis=1)
c1=np.multiply(c,self.lattice)
x2=np.sum(c1)
x=x1+x2
energy = -x
return energy

def montecarlostep(self, T):
# complete this function so that it performs a single Monte Carlo step
energy = self.energy()
#the following two lines will select the coordinates of the random spin for you
random_i = np.random.choice(range(0, self.n_rows))
random_j = np.random.choice(range(0, self.n_cols))
#the following line will choose a random number in the range [0,1) for you
random_number = np.random.random()
self.lattice[random_i][random_j] = -self.lattice[random_i][random_j]
energy2 = self.energy()
dE = energy2 - energy

if energy2<energy or random_number<=np.e**(-(energy2-energy)/T):
yyy=True
else:
self.lattice[random_i][random_j] = -self.lattice[random_i][random_j]

self.n_cycles+=1

self.E += self.energy()
self.E2 += self.energy()**2
self.M += self.magnetisation()
self.M2 += self.magnetisation()**2
return self.energy(), self.magnetisation()

def statistics(self):
# complete this function so that it calculates the correct values for the averages of E, E*E (E2), M, M*M (M2), and returns them
avE = self.E / self.n_cycles
avE2 = self.E2 / self.n_cycles
avM = self.M / self.n_cycles
avM2 = self.M / self.n_cycles
return avE, avE2, avM, avM2, self.n_cycles
</pre>

== Temperature effects ==

=== Correction of averaging code ===

A few experiments were performed at different temperatures using ILfinalframe.py, to see how long it usually takes for the system to equilibrate. At low temperatures it usually happened before 1000 steps were performed, thus 1000 Monte Carlo steps was decided to be appropriate cut-off period for calculation of average values of magnetisation and energy. The figures 2-4, showing the experiments performed are shown in the gallery below.

<gallery mode="packed" heights=250 px>
Image:GE715_CMP_1.png|Figure 2: Test result @ 0.5 K
Image:GE715_CMP_2.png|Figure 3: Test result @ 1 K
Image:GE715_CMP_3.png|Figure 4: Test result @ 1.5 K
</gallery>

The definitions of two class functions were changed as shown below to accommodate this change.

<pre>
def montecarlostep(self, T):
# complete this function so that it performs a single Monte Carlo step
energy = self.energy()
#the following two lines will select the coordinates of the random spin for you
random_i = np.random.choice(range(0, self.n_rows))
random_j = np.random.choice(range(0, self.n_cols))
#the following line will choose a random number in the range [0,1) for you
random_number = np.random.random()
self.lattice[random_i][random_j] = -self.lattice[random_i][random_j]
energy2 = self.energy()
dE = energy2 - energy

if energy2<energy or random_number<=np.e**(-(energy2-energy)/T):
yyy=True
else:
self.lattice[random_i][random_j] = -self.lattice[random_i][random_j]

self.n_cycles+=1
if self.n_cycles==1000:
self.E=0
self.E2=0
self.M=0
self.M2=0
self.E += self.energy()
self.E2 += self.energy()**2
self.M += self.magnetisation()
self.M2 += self.magnetisation()**2
return self.energy(), self.magnetisation()

def statistics(self):
# complete this function so that it calculates the correct values for the averages of E, E*E (E2), M, M*M (M2), and returns them
avE = self.E / self.n_cycles
avE2 = self.E2 / self.n_cycles
avM = self.M / self.n_cycles
avM2 = self.M / self.n_cycles
if self.n_cycles>1000:
avE = self.E / (self.n_cycles-1000)
avE2 = self.E2 / (self.n_cycles-1000)
avM = self.M / (self.n_cycles-1000)
avM2 = self.M / (self.n_cycles-1000)
return avE, avE2, avM, avM2, self.n_cycles
</pre>

=== Range of temperatures ===

The experiment was run for 8x8 lattice at a temperature range of 0.25 K-5.00 K, with a point every 0.05 K. The experiments using a different number of total steps were performed and the comparison is shown below in figures 5-7.

<gallery mode="packed" heights=250 px>
Image:GE715_CMP_4.png|Figure 5: Test result for 8x8 lattice with 11000 steps
Image:GE715_CMP_5.png|Figure 6: Test result for 8x8 lattice with 21000 steps
Image:GE715_CMP_6.png|Figure 7: Test result for 8x8 lattice with 101000 steps
</gallery>

There is a small, but noticeable improvement in the smoothness of the curves obtained for energy with using more steps, thus the highest number of steps (101000 steps) was used in the further experiments regardless of the higher need for computational power.

== Lattice size effects ==

With increasing lattice size the energy and magnetisation curves become more realistic, since the long range effects are better accounted for. 7 different lattice sizes were tested and energy and magnetisation for each are shown below in figure 8-14 and 15-21 respectively.

<gallery mode="packed" heights=250 px>
Image:GE715_CMP_7.png|Figure 8: Energy result for 2x2 lattice
Image:GE715_CMP_8.png|Figure 9: Energy result for 4x4 lattice
Image:GE715_CMP_9.png|Figure 10: Energy result for 8x8 lattice
Image:GE715_CMP_10.png|Figure 11: Energy result for 16x16 lattice
Image:GE715_CMP_11.png|Figure 12: Energy result for 32x32 lattice
Image:GE715_CMP_12.png|Figure 13: Energy result for 64x64 lattice
Image:GE715_CMP_13.png|Figure 14: Energy result for 128x128 lattice
Image:GE715_CMP_14.png|Figure 15: Magnetisation result for 2x2 lattice
Image:GE715_CMP_15.png|Figure 16: Magnetisation result for 4x4 lattice
Image:GE715_CMP_16.png|Figure 17: Magnetisation result for 8x8 lattice
Image:GE715_CMP_17.png|Figure 18: Magnetisation result for 16x16 lattice
Image:GE715_CMP_18.png|Figure 19: Magnetisation result for 32x32 lattice
Image:GE715_CMP_19.png|Figure 20: Magnetisation result for 64x64 lattice
Image:GE715_CMP_20.png|Figure 21: Magnetisation result for 128x128 lattice
</gallery>

Overall comparison of energy and magnetisation per spin in each lattice was made and can be seen in figures 22 and 23 below.

<gallery mode="packed" heights=500 px>
Image:GE715_CMP_21.png|Figure 22: Energy comparison
Image:GE715_CMP_22.png|Figure 23: Magnetisation comparison
</gallery>

The code used for producing figure 23 above is shown below. Other codes are analagous. Example for defining ls2 is also shown Others are similar..

<pre>
ls2=loadtxt("2x2(100).dat")

figsize(12,8)
plot(ls2[:,0],ls2[:,3]/4,color="red",label="2x2")
plot(ls4[:,0],ls4[:,3]/16,color="orange",label="4x4")
plot(ls8[:,0],ls8[:,3]/64,color="yellow",label="8x8")
plot(ls16[:,0],ls16[:,3]/256,color="lime",label="16x16")
plot(ls32[:,0],ls32[:,3]/1024,color="blue",label="32x32")
plot(ls64[:,0],ls64[:,3]/4096,color="indigo",label="64x64")
plot(ls128[:,0],ls128[:,3]/16384,color="violet",label="128x128")
title("Magnetisation per spin")
xlabel("Temperature (K)")
ylabel("Magnetisation per spin")
legend()
savefig("Magnetisation.png")
show()
</pre>

Significant changes may be observed from 2x2 lattice up to 16x16 lattice, later the additional long range effects become less noticeable, but new effects are still present. It seems that most of the long range effects are captured by 16x16 lattice, but higher dimensions bring additional improvement.

== Heat Capacity ==

This formula holds true by definition: <math>C = \frac{\partial \left\langle E\right\rangle}{\partial T}</math>.

Variance of energy can be calculated in the following way: <math>\mathrm{Var}[X] = \left\langle E^2\right\rangle - \left\langle E\right\rangle^2</math>.

Since expected energy can also be expressed as partition function: <math>\left\langle E\right\rangle = \frac{\partial}{\partial \beta}lnZ</math>, variance may be expressed as <math>\mathrm{Var}[X] = \frac{\partial ^2}{\partial \beta ^2}lnZ = - \frac{\partial \left\langle E\right\rangle}{\partial \beta}</math>.

We can rearrange the equations given, because <math>T = \frac{1}{k_BT}</math>, so that we obtain <math>C = \frac{\mathrm{Var}[E]}{k_B T^2}</math>.

=== Heat capacity of the model ===

Heat capacities of the models were calculated using the script given below (example for 2x2 lattice).

<pre>
figsize(6,4)
plot(ls2[:,0],(ls2[:,2]-ls2[:,1]**2)/ls2[:,0],color="red",label="2x2")
xlabel("Temperature (K)")
ylabel("C")
savefig("C2.png")
show()
</pre>

Graphs of heat capacities are shown in the gallery below (figures 24-30). It can be seen that the model breaks down at low temperatures for bigger lattices, but it still gives a peak heat capacity.

<gallery mode="packed" heights=250 px>
Image:GE715_CMP_23.png|Figure 24: C result for 2x2 lattice
Image:GE715_CMP_24.png|Figure 25: C result for 4x4 lattice
Image:GE715_CMP_25.png|Figure 26: C result for 8x8 lattice
Image:GE715_CMP_26.png|Figure 27: C result for 16x16 lattice
Image:GE715_CMP_27.png|Figure 28: C result for 32x32 lattice
Image:GE715_CMP_28.png|Figure 29: C result for 64x64 lattice
Image:GE715_CMP_29.png|Figure 30: C result for 128x128 lattice
</gallery>

Comparison of heat capacity per spin is given in the figure 31 below.

<gallery mode="packed" heights=500 px>
Image:GE715_CMP_30.png|Figure 31: Heat capacity comparison
</gallery>

== Curie temperature ==

Curie temperature can be deduced from where there is a peak in heat capacity per spin. The graphs of heat capacity vs temperature is shown below (figures 32 and 33) for two different lattice sizes, once using Python script and once using C++ data provided.

<gallery mode="packed" heights=500 px>
Image:GE715_CMP_31.png|Figure 32: 16x16 lattice
Image:GE715_CMP_32.png|Figure 33: 64x64 lattice
</gallery>

=== Polynomial fit ===

The C++ data was also fitted with a polynomial, firstly over the whole range and afterwards just in the peak range. Fitting a polynomial proved to be hard as expected (since data is close to delta-like function), but it became more feasible when smaller range was used. Examples of both fits are shown in figures 34 and 35 below.

<gallery mode="packed" heights=500 px>
Image:GE715_CMP_33.png|Figure 34: Whole region fit
Image:GE715_CMP_34.png|Figure 35: Peak region fit
</gallery>

=== Finding Curie temperature ===

The last step was to find peak heat capacity for different lattice sizes in C++ data. The results are shown in the table below.

{| class="wikitable"
|-
! Lattice size
! Curie temperature
|-
|2x2
|2.48 K
|-
|4x4
|2.40 K
|-
|8x8
|2.32 K
|-
|16x16
|2.25 K
|-
|32x32
|2.30 K
|-
|64x64
|2.25 K
|}

The curie temperature for infinite lattice is known to be calculated from that of a finite lattice using equation <math>T_{C, L} = \frac{A}{L} + T_{C,\infty}</math>, where <math>T_{C, L}</math> is the Curie temperature of an <math>L\times L</math>lattice, <math>T_{C,\infty}</math> is the Curie temperature of an infinite lattice, and <math>A</math> is a constant. This was used to draw a graph shown in figure 36 below.

<gallery mode="packed" heights=500 px>
Image:GE715_CMP_35.png|Figure 36: Curie temperature for infinite lattice
</gallery>

The Curie temperature for infinite lattice was found to be <math>T=2.26 K</math>.

The code used for drawing the graph in figure 36 is shown below.

<pre>
L=array([1/2,1/4,1/8,1/16,1/32,1/64])
T=array([2.48,2.4,2.32,2.25,2.3,2.25])
X=array([1/2,0])
figsize(9,6)
plot(L,T,marker="D",linestyle=" ")
xlabel("1/L")
ylabel("T")
a=polyfit(L,T,1)
ylim(2.23,2.5)
plot(X,a[0]*X+a[1],color="lime")
title("Curie temperature")
savefig("Curie.png")
show()
</pre>

Rep:CMPY3 GE715

2017-11-22T11:42:07Z

Ge715: /* Finding Curie temperature */ Code

== Ising Model ==

Ising model is built on the assumption that a lattice that can be of any dimensionality contains spins which can either be up (+1 or down (-1) and interactions can be studied using this model. The model is extremely simplified, since in real systems the spins can point in any direction, but the basic properties of magnetic systems, such as magnetisation, energy levels and Curie temperature, may be studied using this simplification.

The interaction energy in Ising model is defined as <math>E=- \frac{1}{2} J \sum_i^N \sum_{j\ \in\ \mathrm{neighbours}\left(i\right)} s_i s_j</math>. The energy will be the lowest possible when all the spins point in the same direction. In that case <math>s_i s_j=1</math> for all possible <math>i,j</math> pairs. The number of pairs for each <math>i</math> depends on the number of dimensions.

<math>P=2D</math>, where <math>P</math> represents the number of pairs and <math>D</math> the number of dimensions. The lowest possible energy will thus be <math>E=- \frac{1}{2} JNP</math>, where <math>N</math> represents the total number of spins. Thus the lowest possible energy may also be simplified as <math>E=-DNJ</math>, as requested.

The multiplicity of this particular state is 2 (all spins up or all spins down), thus <math>\Omega=2</math>, entropy of the system is <math>S=k_b ln \Omega=9.6*10^{-24} JK^{-1}</math>.

=== Flip of a spin ===

If one of the spins flips, the energy increases. In the given case (<math>D=3, N=1000</math>), the <math>E_i=-3000 J</math>, while the flip of one spin increases the energy by <math>12 J</math> (6 pairs, <math>2 J</math> per pair), so that <math>E_f=-2988 J</math>. Note that <math>J</math> does not represent joules in this case.

Multiplicity has increased to <math>\Omega=2000</math>, thus <math>S=1.05*10^{-22} JK^{-1}</math>.

<math>\Delta S=1.04*10^{-22} JK^{-1}</math>.

This means that entropy has increased by more than 100 times.

=== Magnetisation in the sketch ===

The 1D system shown has the magnetisation of <math>M=1</math>. This is also true for the 2D system shown. For the model system with <math>N=1000, D=3</math> the expected magnetisation at absolute zero would be <math>M=\pm 1000</math>.

== Code for magnetisation and energy ==

Python code was written to calculate the energy and magnetisation of a configuration and tested a few times. The test results are shown in figure 1 below.

<gallery mode="packed" heights=800 px>
Image:GE715_CMP_0.png|Figure 1: Test results
</gallery>

== Accelerating the code ==

The code at the beginning was written using a loop for magnetisation and a double loop for energy. 2000 Monte Carlo steps needed was <math>M=9.32 s \pm 0.10 s</math>. After optimizing magnetisation only this has dropped to <math>M=8.88 s \pm 0.04 s</math>, while after optimizing both magnetisation and energy with using Numpy's sum, roll and multiply functions, the time was further reduced to <math>M=0.0796 s \pm 0.018 s</math>. It was later realised that only one roll in each direction is needed to produce the correct energy, so the final time for 2000 steps was successfully reduced to <math>M=0.418 s \pm 0.005 s</math>.

The whole of the code used is given below.

<pre>
import numpy as np

class IsingLattice:

E = 0.0
E2 = 0.0
M = 0.0
M2 = 0.0
energy_x=0
energy2_x=0
magnetisation_x=0
magnetisation2_x=0

n_cycles = 0

def __init__(self, n_rows, n_cols):
self.n_rows = n_rows
self.n_cols = n_cols
self.lattice = np.random.choice([-1,1], size=(n_rows, n_cols))

def energy(self):
a=np.roll(self.lattice,1,axis=0)
a1=np.multiply(a,self.lattice)
x1=np.sum(a1)
c=np.roll(self.lattice,1,axis=1)
c1=np.multiply(c,self.lattice)
x2=np.sum(c1)
x=x1+x2
energy = -x
return energy

def magnetisation(self):
"Return the total magnetisation of the current lattice configuration."
return np.sum(self.lattice)

def energy_flip(self,random_i,random_j):
"Return the total energy of the current lattice after flip configuration."
self.lattice[random_i][random_j]=-self.lattice[random_i][random_j]
a=np.roll(self.lattice,1,axis=0)
a1=np.multiply(a,self.lattice)
x1=np.sum(a1)
c=np.roll(self.lattice,1,axis=1)
c1=np.multiply(c,self.lattice)
x2=np.sum(c1)
x=x1+x2
energy = -x
return energy

def montecarlostep(self, T):
# complete this function so that it performs a single Monte Carlo step
energy = self.energy()
#the following two lines will select the coordinates of the random spin for you
random_i = np.random.choice(range(0, self.n_rows))
random_j = np.random.choice(range(0, self.n_cols))
#the following line will choose a random number in the range [0,1) for you
random_number = np.random.random()
self.lattice[random_i][random_j] = -self.lattice[random_i][random_j]
energy2 = self.energy()
dE = energy2 - energy

if energy2<energy or random_number<=np.e**(-(energy2-energy)/T):
yyy=True
else:
self.lattice[random_i][random_j] = -self.lattice[random_i][random_j]

self.n_cycles+=1

self.E += self.energy()
self.E2 += self.energy()**2
self.M += self.magnetisation()
self.M2 += self.magnetisation()**2
return self.energy(), self.magnetisation()

def statistics(self):
# complete this function so that it calculates the correct values for the averages of E, E*E (E2), M, M*M (M2), and returns them
avE = self.E / self.n_cycles
avE2 = self.E2 / self.n_cycles
avM = self.M / self.n_cycles
avM2 = self.M / self.n_cycles
return avE, avE2, avM, avM2, self.n_cycles
</pre>

== Temperature effects ==

=== Correction of averaging code ===

A few experiments were performed at different temperatures using ILfinalframe.py, to see how long it usually takes for the system to equilibrate. At low temperatures it usually happened before 1000 steps were performed, thus 1000 Monte Carlo steps was decided to be appropriate cut-off period for calculation of average values of magnetisation and energy. The figures 2-4, showing the experiments performed are shown in the gallery below.

<gallery mode="packed" heights=250 px>
Image:GE715_CMP_1.png|Figure 2: Test result @ 0.5 K
Image:GE715_CMP_2.png|Figure 3: Test result @ 1 K
Image:GE715_CMP_3.png|Figure 4: Test result @ 1.5 K
</gallery>

The definitions of two class functions were changed as shown below to accommodate this change.

<pre>
def montecarlostep(self, T):
# complete this function so that it performs a single Monte Carlo step
energy = self.energy()
#the following two lines will select the coordinates of the random spin for you
random_i = np.random.choice(range(0, self.n_rows))
random_j = np.random.choice(range(0, self.n_cols))
#the following line will choose a random number in the range [0,1) for you
random_number = np.random.random()
self.lattice[random_i][random_j] = -self.lattice[random_i][random_j]
energy2 = self.energy()
dE = energy2 - energy

if energy2<energy or random_number<=np.e**(-(energy2-energy)/T):
yyy=True
else:
self.lattice[random_i][random_j] = -self.lattice[random_i][random_j]

self.n_cycles+=1
if self.n_cycles==1000:
self.E=0
self.E2=0
self.M=0
self.M2=0
self.E += self.energy()
self.E2 += self.energy()**2
self.M += self.magnetisation()
self.M2 += self.magnetisation()**2
return self.energy(), self.magnetisation()

def statistics(self):
# complete this function so that it calculates the correct values for the averages of E, E*E (E2), M, M*M (M2), and returns them
avE = self.E / self.n_cycles
avE2 = self.E2 / self.n_cycles
avM = self.M / self.n_cycles
avM2 = self.M / self.n_cycles
if self.n_cycles>1000:
avE = self.E / (self.n_cycles-1000)
avE2 = self.E2 / (self.n_cycles-1000)
avM = self.M / (self.n_cycles-1000)
avM2 = self.M / (self.n_cycles-1000)
return avE, avE2, avM, avM2, self.n_cycles
</pre>

=== Range of temperatures ===

The experiment was run for 8x8 lattice at a temperature range of 0.25 K-5.00 K, with a point every 0.05 K. The experiments using a different number of total steps were performed and the comparison is shown below in figures 5-7.

<gallery mode="packed" heights=250 px>
Image:GE715_CMP_4.png|Figure 5: Test result for 8x8 lattice with 11000 steps
Image:GE715_CMP_5.png|Figure 6: Test result for 8x8 lattice with 21000 steps
Image:GE715_CMP_6.png|Figure 7: Test result for 8x8 lattice with 101000 steps
</gallery>

There is a small, but noticeable improvement in the smoothness of the curves obtained for energy with using more steps, thus the highest number of steps (101000 steps) was used in the further experiments regardless of the higher need for computational power.

== Lattice size effects ==

With increasing lattice size the energy and magnetisation curves become more realistic, since the long range effects are better accounted for. 7 different lattice sizes were tested and energy and magnetisation for each are shown below in figure 8-14 and 15-21 respectively.

<gallery mode="packed" heights=250 px>
Image:GE715_CMP_7.png|Figure 8: Energy result for 2x2 lattice
Image:GE715_CMP_8.png|Figure 9: Energy result for 4x4 lattice
Image:GE715_CMP_9.png|Figure 10: Energy result for 8x8 lattice
Image:GE715_CMP_10.png|Figure 11: Energy result for 16x16 lattice
Image:GE715_CMP_11.png|Figure 12: Energy result for 32x32 lattice
Image:GE715_CMP_12.png|Figure 13: Energy result for 64x64 lattice
Image:GE715_CMP_13.png|Figure 14: Energy result for 128x128 lattice
Image:GE715_CMP_14.png|Figure 15: Magnetisation result for 2x2 lattice
Image:GE715_CMP_15.png|Figure 16: Magnetisation result for 4x4 lattice
Image:GE715_CMP_16.png|Figure 17: Magnetisation result for 8x8 lattice
Image:GE715_CMP_17.png|Figure 18: Magnetisation result for 16x16 lattice
Image:GE715_CMP_18.png|Figure 19: Magnetisation result for 32x32 lattice
Image:GE715_CMP_19.png|Figure 20: Magnetisation result for 64x64 lattice
Image:GE715_CMP_20.png|Figure 21: Magnetisation result for 128x128 lattice
</gallery>

Overall comparison of energy and magnetisation per spin in each lattice was made and can be seen in figures 22 and 23 below.

<gallery mode="packed" heights=500 px>
Image:GE715_CMP_21.png|Figure 22: Energy comparison
Image:GE715_CMP_22.png|Figure 23: Magnetisation comparison
</gallery>

The code used for producing figure 23 above is shown below. Other codes are analagous. Example for defining ls2 is also shown Others are similar..

<pre>
ls2=loadtxt("2x2(100).dat")

figsize(12,8)
plot(ls2[:,0],ls2[:,3]/4,color="red",label="2x2")
plot(ls4[:,0],ls4[:,3]/16,color="orange",label="4x4")
plot(ls8[:,0],ls8[:,3]/64,color="yellow",label="8x8")
plot(ls16[:,0],ls16[:,3]/256,color="lime",label="16x16")
plot(ls32[:,0],ls32[:,3]/1024,color="blue",label="32x32")
plot(ls64[:,0],ls64[:,3]/4096,color="indigo",label="64x64")
plot(ls128[:,0],ls128[:,3]/16384,color="violet",label="128x128")
title("Magnetisation per spin")
xlabel("Temperature (K)")
ylabel("Magnetisation per spin")
legend()
savefig("Magnetisation.png")
show()
</pre>

Significant changes may be observed from 2x2 lattice up to 16x16 lattice, later the additional long range effects become less noticeable, but new effects are still present. It seems that most of the long range effects are captured by 16x16 lattice, but higher dimensions bring additional improvement.

== Heat Capacity ==

This formula holds true by definition: <math>C = \frac{\partial \left\langle E\right\rangle}{\partial T}</math>.

Variance of energy can be calculated in the following way: <math>\mathrm{Var}[X] = \left\langle E^2\right\rangle - \left\langle E\right\rangle^2</math>.

Since expected energy can also be expressed as partition function: <math>\left\langle E\right\rangle = \frac{\partial}{\partial \beta}lnZ</math>, variance may be expressed as <math>\mathrm{Var}[X] = \frac{\partial ^2}{\partial \beta ^2}lnZ = - \frac{\partial \left\langle E\right\rangle}{\partial \beta}</math>.

We can rearrange the equations given, because <math>T = \frac{1}{k_BT}</math>, so that we obtain <math>C = \frac{\mathrm{Var}[E]}{k_B T^2}</math>.

=== Heat capacity of the model ===

Heat capacities of the models were calculated using the script given below (example for 2x2 lattice).

<pre>
figsize(6,4)
plot(ls2[:,0],(ls2[:,2]-ls2[:,1]**2)/ls2[:,0],color="red",label="2x2")
xlabel("Temperature (K)")
ylabel("C")
savefig("C2.png")
show()
</pre>

Graphs of heat capacities are shown in the gallery below (figures 24-30). It can be seen that the model breaks down at low temperatures for bigger lattices, but it still gives a peak heat capacity.

<gallery mode="packed" heights=250 px>
Image:GE715_CMP_23.png|Figure 24: C result for 2x2 lattice
Image:GE715_CMP_24.png|Figure 25: C result for 4x4 lattice
Image:GE715_CMP_25.png|Figure 26: C result for 8x8 lattice
Image:GE715_CMP_26.png|Figure 27: C result for 16x16 lattice
Image:GE715_CMP_27.png|Figure 28: C result for 32x32 lattice
Image:GE715_CMP_28.png|Figure 29: C result for 64x64 lattice
Image:GE715_CMP_29.png|Figure 30: C result for 128x128 lattice
</gallery>

Comparison of heat capacity per spin is given in the figure 31 below.

<gallery mode="packed" heights=500 px>
Image:GE715_CMP_30.png|Figure 31: Heat capacity comparison
</gallery>

== Curie temperature ==

Curie temperature can be deduced from where there is a peak in heat capacity per spin. The graphs of heat capacity vs temperature is shown below (figures 32 and 33) for two different lattice sizes, once using Python script and once using C++ data provided.

<gallery mode="packed" heights=500 px>
Image:GE715_CMP_31.png|Figure 32: 16x16 lattice
Image:GE715_CMP_32.png|Figure 33: 64x64 lattice
</gallery>

=== Polynomial fit ===

The C++ data was also fitted with a polynomial, firstly over the whole range and afterwards just in the peak range. Fitting a polynomial proved to be hard as expected (since data is close to delta-like function), but it became more feasible when smaller range was used. Examples of both fits are shown in figures 34 and 35 below.

<gallery mode="packed" heights=500 px>
Image:GE715_CMP_33.png|Figure 34: Whole region fit
Image:GE715_CMP_34.png|Figure 35: Peak region fit
</gallery>

=== Finding Curie temperature ===

The last step was to find peak heat capacity for different lattice sizes in C++ data. The results are shown in the table below.

{| class="wikitable"
|-
! Lattice size
! Curie temperature
|-
|2x2
|2.48 K
|-
|4x4
|2.40 K
|-
|8x8
|2.32 K
|-
|16x16
|2.25 K
|-
|32x32
|2.30 K
|-
|64x64
|2.25 K
|}

The curie temperature for infinite lattice is known to be calculated from that of a finite lattice using equation <math>T_{C, L} = \frac{A}{L} + T_{C,\infty}</math>, where <math>T_{C, L}</math> is the Curie temperature of an <math>L\times L</math>lattice, <math>T_{C,\infty}</math> is the Curie temperature of an infinite lattice, and <math>A</math> is a constant. This was used to draw a graph shown in figure 36 below.

<gallery mode="packed" heights=500 px>
Image:GE715_CMP_35.png|Figure 36: Curie temperature for infinite lattice
</gallery>

The Curie temperature for infinite lattice was found to be <math>T=2.26 K</math>.

The code used for drawing the graph in figure 36 is shown below.

<pre>
L=array([1/2,1/4,1/8,1/16,1/32,1/64])
T=array([2.48,2.4,2.32,2.25,2.3,2.25])
X=array([1/2,0])
figsize(9,6)
plot(L,T,marker="D",linestyle=" ")
xlabel("1/L")
ylabel("T")
a=polyfit(L,T,1)
ylim(2.23,2.5)
plot(X,a[0]*X+a[1],color="lime")
title("Curie temperature")
savefig("Curie.png")
show()
</pre>

Rep:CMPY3 GE715

2017-11-22T11:40:41Z

Ge715: /* Ising Model */ Introduction

== Ising Model ==

Ising model is built on the assumption that a lattice that can be of any dimensionality contains spins which can either be up (+1 or down (-1) and interactions can be studied using this model. The model is extremely simplified, since in real systems the spins can point in any direction, but the basic properties of magnetic systems, such as magnetisation, energy levels and Curie temperature, may be studied using this simplification.

The interaction energy in Ising model is defined as <math>E=- \frac{1}{2} J \sum_i^N \sum_{j\ \in\ \mathrm{neighbours}\left(i\right)} s_i s_j</math>. The energy will be the lowest possible when all the spins point in the same direction. In that case <math>s_i s_j=1</math> for all possible <math>i,j</math> pairs. The number of pairs for each <math>i</math> depends on the number of dimensions.

<math>P=2D</math>, where <math>P</math> represents the number of pairs and <math>D</math> the number of dimensions. The lowest possible energy will thus be <math>E=- \frac{1}{2} JNP</math>, where <math>N</math> represents the total number of spins. Thus the lowest possible energy may also be simplified as <math>E=-DNJ</math>, as requested.

The multiplicity of this particular state is 2 (all spins up or all spins down), thus <math>\Omega=2</math>, entropy of the system is <math>S=k_b ln \Omega=9.6*10^{-24} JK^{-1}</math>.

=== Flip of a spin ===

If one of the spins flips, the energy increases. In the given case (<math>D=3, N=1000</math>), the <math>E_i=-3000 J</math>, while the flip of one spin increases the energy by <math>12 J</math> (6 pairs, <math>2 J</math> per pair), so that <math>E_f=-2988 J</math>. Note that <math>J</math> does not represent joules in this case.

Multiplicity has increased to <math>\Omega=2000</math>, thus <math>S=1.05*10^{-22} JK^{-1}</math>.

<math>\Delta S=1.04*10^{-22} JK^{-1}</math>.

This means that entropy has increased by more than 100 times.

=== Magnetisation in the sketch ===

The 1D system shown has the magnetisation of <math>M=1</math>. This is also true for the 2D system shown. For the model system with <math>N=1000, D=3</math> the expected magnetisation at absolute zero would be <math>M=\pm 1000</math>.

== Code for magnetisation and energy ==

Python code was written to calculate the energy and magnetisation of a configuration and tested a few times. The test results are shown in figure 1 below.

<gallery mode="packed" heights=800 px>
Image:GE715_CMP_0.png|Figure 1: Test results
</gallery>

== Accelerating the code ==

The code at the beginning was written using a loop for magnetisation and a double loop for energy. 2000 Monte Carlo steps needed was <math>M=9.32 s \pm 0.10 s</math>. After optimizing magnetisation only this has dropped to <math>M=8.88 s \pm 0.04 s</math>, while after optimizing both magnetisation and energy with using Numpy's sum, roll and multiply functions, the time was further reduced to <math>M=0.0796 s \pm 0.018 s</math>. It was later realised that only one roll in each direction is needed to produce the correct energy, so the final time for 2000 steps was successfully reduced to <math>M=0.418 s \pm 0.005 s</math>.

The whole of the code used is given below.

<pre>
import numpy as np

class IsingLattice:

E = 0.0
E2 = 0.0
M = 0.0
M2 = 0.0
energy_x=0
energy2_x=0
magnetisation_x=0
magnetisation2_x=0

n_cycles = 0

def __init__(self, n_rows, n_cols):
self.n_rows = n_rows
self.n_cols = n_cols
self.lattice = np.random.choice([-1,1], size=(n_rows, n_cols))

def energy(self):
a=np.roll(self.lattice,1,axis=0)
a1=np.multiply(a,self.lattice)
x1=np.sum(a1)
c=np.roll(self.lattice,1,axis=1)
c1=np.multiply(c,self.lattice)
x2=np.sum(c1)
x=x1+x2
energy = -x
return energy

def magnetisation(self):
"Return the total magnetisation of the current lattice configuration."
return np.sum(self.lattice)

def energy_flip(self,random_i,random_j):
"Return the total energy of the current lattice after flip configuration."
self.lattice[random_i][random_j]=-self.lattice[random_i][random_j]
a=np.roll(self.lattice,1,axis=0)
a1=np.multiply(a,self.lattice)
x1=np.sum(a1)
c=np.roll(self.lattice,1,axis=1)
c1=np.multiply(c,self.lattice)
x2=np.sum(c1)
x=x1+x2
energy = -x
return energy

def montecarlostep(self, T):
# complete this function so that it performs a single Monte Carlo step
energy = self.energy()
#the following two lines will select the coordinates of the random spin for you
random_i = np.random.choice(range(0, self.n_rows))
random_j = np.random.choice(range(0, self.n_cols))
#the following line will choose a random number in the range [0,1) for you
random_number = np.random.random()
self.lattice[random_i][random_j] = -self.lattice[random_i][random_j]
energy2 = self.energy()
dE = energy2 - energy

if energy2<energy or random_number<=np.e**(-(energy2-energy)/T):
yyy=True
else:
self.lattice[random_i][random_j] = -self.lattice[random_i][random_j]

self.n_cycles+=1

self.E += self.energy()
self.E2 += self.energy()**2
self.M += self.magnetisation()
self.M2 += self.magnetisation()**2
return self.energy(), self.magnetisation()

def statistics(self):
# complete this function so that it calculates the correct values for the averages of E, E*E (E2), M, M*M (M2), and returns them
avE = self.E / self.n_cycles
avE2 = self.E2 / self.n_cycles
avM = self.M / self.n_cycles
avM2 = self.M / self.n_cycles
return avE, avE2, avM, avM2, self.n_cycles
</pre>

== Temperature effects ==

=== Correction of averaging code ===

A few experiments were performed at different temperatures using ILfinalframe.py, to see how long it usually takes for the system to equilibrate. At low temperatures it usually happened before 1000 steps were performed, thus 1000 Monte Carlo steps was decided to be appropriate cut-off period for calculation of average values of magnetisation and energy. The figures 2-4, showing the experiments performed are shown in the gallery below.

<gallery mode="packed" heights=250 px>
Image:GE715_CMP_1.png|Figure 2: Test result @ 0.5 K
Image:GE715_CMP_2.png|Figure 3: Test result @ 1 K
Image:GE715_CMP_3.png|Figure 4: Test result @ 1.5 K
</gallery>

The definitions of two class functions were changed as shown below to accommodate this change.

<pre>
def montecarlostep(self, T):
# complete this function so that it performs a single Monte Carlo step
energy = self.energy()
#the following two lines will select the coordinates of the random spin for you
random_i = np.random.choice(range(0, self.n_rows))
random_j = np.random.choice(range(0, self.n_cols))
#the following line will choose a random number in the range [0,1) for you
random_number = np.random.random()
self.lattice[random_i][random_j] = -self.lattice[random_i][random_j]
energy2 = self.energy()
dE = energy2 - energy

if energy2<energy or random_number<=np.e**(-(energy2-energy)/T):
yyy=True
else:
self.lattice[random_i][random_j] = -self.lattice[random_i][random_j]

self.n_cycles+=1
if self.n_cycles==1000:
self.E=0
self.E2=0
self.M=0
self.M2=0
self.E += self.energy()
self.E2 += self.energy()**2
self.M += self.magnetisation()
self.M2 += self.magnetisation()**2
return self.energy(), self.magnetisation()

def statistics(self):
# complete this function so that it calculates the correct values for the averages of E, E*E (E2), M, M*M (M2), and returns them
avE = self.E / self.n_cycles
avE2 = self.E2 / self.n_cycles
avM = self.M / self.n_cycles
avM2 = self.M / self.n_cycles
if self.n_cycles>1000:
avE = self.E / (self.n_cycles-1000)
avE2 = self.E2 / (self.n_cycles-1000)
avM = self.M / (self.n_cycles-1000)
avM2 = self.M / (self.n_cycles-1000)
return avE, avE2, avM, avM2, self.n_cycles
</pre>

=== Range of temperatures ===

The experiment was run for 8x8 lattice at a temperature range of 0.25 K-5.00 K, with a point every 0.05 K. The experiments using a different number of total steps were performed and the comparison is shown below in figures 5-7.

<gallery mode="packed" heights=250 px>
Image:GE715_CMP_4.png|Figure 5: Test result for 8x8 lattice with 11000 steps
Image:GE715_CMP_5.png|Figure 6: Test result for 8x8 lattice with 21000 steps
Image:GE715_CMP_6.png|Figure 7: Test result for 8x8 lattice with 101000 steps
</gallery>

There is a small, but noticeable improvement in the smoothness of the curves obtained for energy with using more steps, thus the highest number of steps (101000 steps) was used in the further experiments regardless of the higher need for computational power.

== Lattice size effects ==

With increasing lattice size the energy and magnetisation curves become more realistic, since the long range effects are better accounted for. 7 different lattice sizes were tested and energy and magnetisation for each are shown below in figure 8-14 and 15-21 respectively.

<gallery mode="packed" heights=250 px>
Image:GE715_CMP_7.png|Figure 8: Energy result for 2x2 lattice
Image:GE715_CMP_8.png|Figure 9: Energy result for 4x4 lattice
Image:GE715_CMP_9.png|Figure 10: Energy result for 8x8 lattice
Image:GE715_CMP_10.png|Figure 11: Energy result for 16x16 lattice
Image:GE715_CMP_11.png|Figure 12: Energy result for 32x32 lattice
Image:GE715_CMP_12.png|Figure 13: Energy result for 64x64 lattice
Image:GE715_CMP_13.png|Figure 14: Energy result for 128x128 lattice
Image:GE715_CMP_14.png|Figure 15: Magnetisation result for 2x2 lattice
Image:GE715_CMP_15.png|Figure 16: Magnetisation result for 4x4 lattice
Image:GE715_CMP_16.png|Figure 17: Magnetisation result for 8x8 lattice
Image:GE715_CMP_17.png|Figure 18: Magnetisation result for 16x16 lattice
Image:GE715_CMP_18.png|Figure 19: Magnetisation result for 32x32 lattice
Image:GE715_CMP_19.png|Figure 20: Magnetisation result for 64x64 lattice
Image:GE715_CMP_20.png|Figure 21: Magnetisation result for 128x128 lattice
</gallery>

Overall comparison of energy and magnetisation per spin in each lattice was made and can be seen in figures 22 and 23 below.

<gallery mode="packed" heights=500 px>
Image:GE715_CMP_21.png|Figure 22: Energy comparison
Image:GE715_CMP_22.png|Figure 23: Magnetisation comparison
</gallery>

The code used for producing figure 23 above is shown below. Other codes are analagous. Example for defining ls2 is also shown Others are similar..

<pre>
ls2=loadtxt("2x2(100).dat")

figsize(12,8)
plot(ls2[:,0],ls2[:,3]/4,color="red",label="2x2")
plot(ls4[:,0],ls4[:,3]/16,color="orange",label="4x4")
plot(ls8[:,0],ls8[:,3]/64,color="yellow",label="8x8")
plot(ls16[:,0],ls16[:,3]/256,color="lime",label="16x16")
plot(ls32[:,0],ls32[:,3]/1024,color="blue",label="32x32")
plot(ls64[:,0],ls64[:,3]/4096,color="indigo",label="64x64")
plot(ls128[:,0],ls128[:,3]/16384,color="violet",label="128x128")
title("Magnetisation per spin")
xlabel("Temperature (K)")
ylabel("Magnetisation per spin")
legend()
savefig("Magnetisation.png")
show()
</pre>

Significant changes may be observed from 2x2 lattice up to 16x16 lattice, later the additional long range effects become less noticeable, but new effects are still present. It seems that most of the long range effects are captured by 16x16 lattice, but higher dimensions bring additional improvement.

== Heat Capacity ==

This formula holds true by definition: <math>C = \frac{\partial \left\langle E\right\rangle}{\partial T}</math>.

Variance of energy can be calculated in the following way: <math>\mathrm{Var}[X] = \left\langle E^2\right\rangle - \left\langle E\right\rangle^2</math>.

Since expected energy can also be expressed as partition function: <math>\left\langle E\right\rangle = \frac{\partial}{\partial \beta}lnZ</math>, variance may be expressed as <math>\mathrm{Var}[X] = \frac{\partial ^2}{\partial \beta ^2}lnZ = - \frac{\partial \left\langle E\right\rangle}{\partial \beta}</math>.

We can rearrange the equations given, because <math>T = \frac{1}{k_BT}</math>, so that we obtain <math>C = \frac{\mathrm{Var}[E]}{k_B T^2}</math>.

=== Heat capacity of the model ===

Heat capacities of the models were calculated using the script given below (example for 2x2 lattice).

<pre>
figsize(6,4)
plot(ls2[:,0],(ls2[:,2]-ls2[:,1]**2)/ls2[:,0],color="red",label="2x2")
xlabel("Temperature (K)")
ylabel("C")
savefig("C2.png")
show()
</pre>

Graphs of heat capacities are shown in the gallery below (figures 24-30). It can be seen that the model breaks down at low temperatures for bigger lattices, but it still gives a peak heat capacity.

<gallery mode="packed" heights=250 px>
Image:GE715_CMP_23.png|Figure 24: C result for 2x2 lattice
Image:GE715_CMP_24.png|Figure 25: C result for 4x4 lattice
Image:GE715_CMP_25.png|Figure 26: C result for 8x8 lattice
Image:GE715_CMP_26.png|Figure 27: C result for 16x16 lattice
Image:GE715_CMP_27.png|Figure 28: C result for 32x32 lattice
Image:GE715_CMP_28.png|Figure 29: C result for 64x64 lattice
Image:GE715_CMP_29.png|Figure 30: C result for 128x128 lattice
</gallery>

Comparison of heat capacity per spin is given in the figure 31 below.

<gallery mode="packed" heights=500 px>
Image:GE715_CMP_30.png|Figure 31: Heat capacity comparison
</gallery>

== Curie temperature ==

Curie temperature can be deduced from where there is a peak in heat capacity per spin. The graphs of heat capacity vs temperature is shown below (figures 32 and 33) for two different lattice sizes, once using Python script and once using C++ data provided.

<gallery mode="packed" heights=500 px>
Image:GE715_CMP_31.png|Figure 32: 16x16 lattice
Image:GE715_CMP_32.png|Figure 33: 64x64 lattice
</gallery>

=== Polynomial fit ===

The C++ data was also fitted with a polynomial, firstly over the whole range and afterwards just in the peak range. Fitting a polynomial proved to be hard as expected (since data is close to delta-like function), but it became more feasible when smaller range was used. Examples of both fits are shown in figures 34 and 35 below.

<gallery mode="packed" heights=500 px>
Image:GE715_CMP_33.png|Figure 34: Whole region fit
Image:GE715_CMP_34.png|Figure 35: Peak region fit
</gallery>

=== Finding Curie temperature ===

The last step was to find peak heat capacity for different lattice sizes in C++ data. The results are shown in the table below.

{| class="wikitable"
|-
! Lattice size
! Curie temperature
|-
|2x2
|2.48 K
|-
|4x4
|2.40 K
|-
|8x8
|2.32 K
|-
|16x16
|2.25 K
|-
|32x32
|2.30 K
|-
|64x64
|2.25 K
|}

The curie temperature for infinite lattice is known to be calculated from that of a finite lattice using equation <math>T_{C, L} = \frac{A}{L} + T_{C,\infty}</math>, where <math>T_{C, L}</math> is the Curie temperature of an <math>L\times L</math>lattice, <math>T_{C,\infty}</math> is the Curie temperature of an infinite lattice, and <math>A</math> is a constant. This was used to draw a graph shown in figure 36 below.

<gallery mode="packed" heights=500 px>
Image:GE715_CMP_35.png|Figure 36: Curie temperature for infinite lattice
</gallery>

The Curie temperature for infinite lattice was found to be <math>T=2.26 K</math>.

File:GE715 CMP 35.png

2017-11-22T10:48:12Z

Ge715:

Rep:CMPY3 GE715

2017-11-22T10:47:12Z

Ge715:

== Ising Model ==

The interaction energy in Ising model is defined as <math>E=- \frac{1}{2} J \sum_i^N \sum_{j\ \in\ \mathrm{neighbours}\left(i\right)} s_i s_j</math>. The energy will be the lowest possible when all the spins point in the same direction. In that case <math>s_i s_j=1</math> for all possible <math>i,j</math> pairs. The number of pairs for each <math>i</math> depends on the number of dimensions.

<math>P=2D</math>, where <math>P</math> represents the number of pairs and <math>D</math> the number of dimensions. The lowest possible energy will thus be <math>E=- \frac{1}{2} JNP</math>, where <math>N</math> represents the total number of spins. Thus the lowest possible energy may also be simplified as <math>E=-DNJ</math>, as requested.

The multiplicity of this particular state is 2 (all spins up or all spins down), thus <math>\Omega=2</math>, entropy of the system is <math>S=k_b ln \Omega=9.6*10^{-24} JK^{-1}</math>.

=== Flip of a spin ===

If one of the spins flips, the energy increases. In the given case (<math>D=3, N=1000</math>), the <math>E_i=-3000 J</math>, while the flip of one spin increases the energy by <math>12 J</math> (6 pairs, <math>2 J</math> per pair), so that <math>E_f=-2988 J</math>. Note that <math>J</math> does not represent joules in this case.

Multiplicity has increased to <math>\Omega=2000</math>, thus <math>S=1.05*10^{-22} JK^{-1}</math>.

<math>\Delta S=1.04*10^{-22} JK^{-1}</math>.

=== Magnetisation in the sketch ===

The 1D system shown has the magnetisation of <math>M=1</math>. This is also true for the 2D system shown. For the model system with <math>N=1000, D=3</math> the expected magnetisation at absolute zero would be <math>M=\pm 1000</math>.

== Code for magnetisation and energy ==

Python code was written to calculate the energy and magnetisation of a configuration and tested a few times. The test results are shown in figure 1 below.

<gallery mode="packed" heights=800 px>
Image:GE715_CMP_0.png|Figure 1: Test results
</gallery>

== Accelerating the code ==

The code at the beginning was written using a loop for magnetisation and a double loop for energy. 2000 Monte Carlo steps needed was <math>M=9.32 s \pm 0.10 s</math>. After optimizing magnetisation only this has dropped to <math>M=8.88 s \pm 0.04 s</math>, while after optimizing both magnetisation and energy with using Numpy's sum, roll and multiply functions, the time was further reduced to <math>M=0.0796 s \pm 0.018 s</math>. It was later realised that only one roll in each direction is needed to produce the correct energy, so the final time for 2000 steps was successfully reduced to <math>M=0.418 s \pm 0.005 s</math>.

The whole of the code used is given below.

<pre>
import numpy as np

class IsingLattice:

E = 0.0
E2 = 0.0
M = 0.0
M2 = 0.0
energy_x=0
energy2_x=0
magnetisation_x=0
magnetisation2_x=0

n_cycles = 0

def __init__(self, n_rows, n_cols):
self.n_rows = n_rows
self.n_cols = n_cols
self.lattice = np.random.choice([-1,1], size=(n_rows, n_cols))

def energy(self):
a=np.roll(self.lattice,1,axis=0)
a1=np.multiply(a,self.lattice)
x1=np.sum(a1)
c=np.roll(self.lattice,1,axis=1)
c1=np.multiply(c,self.lattice)
x2=np.sum(c1)
x=x1+x2
energy = -x
return energy

def magnetisation(self):
"Return the total magnetisation of the current lattice configuration."
return np.sum(self.lattice)

def energy_flip(self,random_i,random_j):
"Return the total energy of the current lattice after flip configuration."
self.lattice[random_i][random_j]=-self.lattice[random_i][random_j]
a=np.roll(self.lattice,1,axis=0)
a1=np.multiply(a,self.lattice)
x1=np.sum(a1)
c=np.roll(self.lattice,1,axis=1)
c1=np.multiply(c,self.lattice)
x2=np.sum(c1)
x=x1+x2
energy = -x
return energy

def montecarlostep(self, T):
# complete this function so that it performs a single Monte Carlo step
energy = self.energy()
#the following two lines will select the coordinates of the random spin for you
random_i = np.random.choice(range(0, self.n_rows))
random_j = np.random.choice(range(0, self.n_cols))
#the following line will choose a random number in the range [0,1) for you
random_number = np.random.random()
self.lattice[random_i][random_j] = -self.lattice[random_i][random_j]
energy2 = self.energy()
dE = energy2 - energy

if energy2<energy or random_number<=np.e**(-(energy2-energy)/T):
yyy=True
else:
self.lattice[random_i][random_j] = -self.lattice[random_i][random_j]

self.n_cycles+=1

self.E += self.energy()
self.E2 += self.energy()**2
self.M += self.magnetisation()
self.M2 += self.magnetisation()**2
return self.energy(), self.magnetisation()

def statistics(self):
# complete this function so that it calculates the correct values for the averages of E, E*E (E2), M, M*M (M2), and returns them
avE = self.E / self.n_cycles
avE2 = self.E2 / self.n_cycles
avM = self.M / self.n_cycles
avM2 = self.M / self.n_cycles
return avE, avE2, avM, avM2, self.n_cycles
</pre>

== Temperature effects ==

=== Correction of averaging code ===

A few experiments were performed at different temperatures using ILfinalframe.py, to see how long it usually takes for the system to equilibrate. At low temperatures it usually happened before 1000 steps were performed, thus 1000 Monte Carlo steps was decided to be appropriate cut-off period for calculation of average values of magnetisation and energy. The figures 2-4, showing the experiments performed are shown in the gallery below.

<gallery mode="packed" heights=250 px>
Image:GE715_CMP_1.png|Figure 2: Test result @ 0.5 K
Image:GE715_CMP_2.png|Figure 3: Test result @ 1 K
Image:GE715_CMP_3.png|Figure 4: Test result @ 1.5 K
</gallery>

The definitions of two class functions were changed as shown below to accommodate this change.

<pre>
def montecarlostep(self, T):
# complete this function so that it performs a single Monte Carlo step
energy = self.energy()
#the following two lines will select the coordinates of the random spin for you
random_i = np.random.choice(range(0, self.n_rows))
random_j = np.random.choice(range(0, self.n_cols))
#the following line will choose a random number in the range [0,1) for you
random_number = np.random.random()
self.lattice[random_i][random_j] = -self.lattice[random_i][random_j]
energy2 = self.energy()
dE = energy2 - energy

if energy2<energy or random_number<=np.e**(-(energy2-energy)/T):
yyy=True
else:
self.lattice[random_i][random_j] = -self.lattice[random_i][random_j]

self.n_cycles+=1
if self.n_cycles==1000:
self.E=0
self.E2=0
self.M=0
self.M2=0
self.E += self.energy()
self.E2 += self.energy()**2
self.M += self.magnetisation()
self.M2 += self.magnetisation()**2
return self.energy(), self.magnetisation()

def statistics(self):
# complete this function so that it calculates the correct values for the averages of E, E*E (E2), M, M*M (M2), and returns them
avE = self.E / self.n_cycles
avE2 = self.E2 / self.n_cycles
avM = self.M / self.n_cycles
avM2 = self.M / self.n_cycles
if self.n_cycles>1000:
avE = self.E / (self.n_cycles-1000)
avE2 = self.E2 / (self.n_cycles-1000)
avM = self.M / (self.n_cycles-1000)
avM2 = self.M / (self.n_cycles-1000)
return avE, avE2, avM, avM2, self.n_cycles
</pre>

=== Range of temperatures ===

The experiment was run for 8x8 lattice at a temperature range of 0.25 K-5.00 K, with a point every 0.05 K. The experiments using a different number of total steps were performed and the comparison is shown below in figures 5-7.

<gallery mode="packed" heights=250 px>
Image:GE715_CMP_4.png|Figure 5: Test result for 8x8 lattice with 11000 steps
Image:GE715_CMP_5.png|Figure 6: Test result for 8x8 lattice with 21000 steps
Image:GE715_CMP_6.png|Figure 7: Test result for 8x8 lattice with 101000 steps
</gallery>

There is a small, but noticeable improvement in the smoothness of the curves obtained for energy with using more steps, thus the highest number of steps (101000 steps) was used in the further experiments regardless of the higher need for computational power.

== Lattice size effects ==

With increasing lattice size the energy and magnetisation curves become more realistic, since the long range effects are better accounted for. 7 different lattice sizes were tested and energy and magnetisation for each are shown below in figure 8-14 and 15-21 respectively.

<gallery mode="packed" heights=250 px>
Image:GE715_CMP_7.png|Figure 8: Energy result for 2x2 lattice
Image:GE715_CMP_8.png|Figure 9: Energy result for 4x4 lattice
Image:GE715_CMP_9.png|Figure 10: Energy result for 8x8 lattice
Image:GE715_CMP_10.png|Figure 11: Energy result for 16x16 lattice
Image:GE715_CMP_11.png|Figure 12: Energy result for 32x32 lattice
Image:GE715_CMP_12.png|Figure 13: Energy result for 64x64 lattice
Image:GE715_CMP_13.png|Figure 14: Energy result for 128x128 lattice
Image:GE715_CMP_14.png|Figure 15: Magnetisation result for 2x2 lattice
Image:GE715_CMP_15.png|Figure 16: Magnetisation result for 4x4 lattice
Image:GE715_CMP_16.png|Figure 17: Magnetisation result for 8x8 lattice
Image:GE715_CMP_17.png|Figure 18: Magnetisation result for 16x16 lattice
Image:GE715_CMP_18.png|Figure 19: Magnetisation result for 32x32 lattice
Image:GE715_CMP_19.png|Figure 20: Magnetisation result for 64x64 lattice
Image:GE715_CMP_20.png|Figure 21: Magnetisation result for 128x128 lattice
</gallery>

Overall comparison of energy and magnetisation per spin in each lattice was made and can be seen in figures 22 and 23 below.

<gallery mode="packed" heights=500 px>
Image:GE715_CMP_21.png|Figure 22: Energy comparison
Image:GE715_CMP_22.png|Figure 23: Magnetisation comparison
</gallery>

The code used for producing figure 23 above is shown below. Other codes are analagous. Example for defining ls2 is also shown Others are similar..

<pre>
ls2=loadtxt("2x2(100).dat")

figsize(12,8)
plot(ls2[:,0],ls2[:,3]/4,color="red",label="2x2")
plot(ls4[:,0],ls4[:,3]/16,color="orange",label="4x4")
plot(ls8[:,0],ls8[:,3]/64,color="yellow",label="8x8")
plot(ls16[:,0],ls16[:,3]/256,color="lime",label="16x16")
plot(ls32[:,0],ls32[:,3]/1024,color="blue",label="32x32")
plot(ls64[:,0],ls64[:,3]/4096,color="indigo",label="64x64")
plot(ls128[:,0],ls128[:,3]/16384,color="violet",label="128x128")
title("Magnetisation per spin")
xlabel("Temperature (K)")
ylabel("Magnetisation per spin")
legend()
savefig("Magnetisation.png")
show()
</pre>

Significant changes may be observed from 2x2 lattice up to 16x16 lattice, later the additional long range effects become less noticeable, but new effects are still present. It seems that most of the long range effects are captured by 16x16 lattice, but higher dimensions bring additional improvement.

== Heat Capacity ==

This formula holds true by definition: <math>C = \frac{\partial \left\langle E\right\rangle}{\partial T}</math>.

Variance of energy can be calculated in the following way: <math>\mathrm{Var}[X] = \left\langle E^2\right\rangle - \left\langle E\right\rangle^2</math>.

Since expected energy can also be expressed as partition function: <math>\left\langle E\right\rangle = \frac{\partial}{\partial \beta}lnZ</math>, variance may be expressed as <math>\mathrm{Var}[X] = \frac{\partial ^2}{\partial \beta ^2}lnZ = - \frac{\partial \left\langle E\right\rangle}{\partial \beta}</math>.

We can rearrange the equations given, because <math>T = \frac{1}{k_BT}</math>, so that we obtain <math>C = \frac{\mathrm{Var}[E]}{k_B T^2}</math>.

=== Heat capacity of the model ===

Heat capacities of the models were calculated using the script given below (example for 2x2 lattice).

<pre>
figsize(6,4)
plot(ls2[:,0],(ls2[:,2]-ls2[:,1]**2)/ls2[:,0],color="red",label="2x2")
xlabel("Temperature (K)")
ylabel("C")
savefig("C2.png")
show()
</pre>

Graphs of heat capacities are shown in the gallery below (figures 24-30). It can be seen that the model breaks down at low temperatures for bigger lattices, but it still gives a peak heat capacity.

<gallery mode="packed" heights=250 px>
Image:GE715_CMP_23.png|Figure 24: C result for 2x2 lattice
Image:GE715_CMP_24.png|Figure 25: C result for 4x4 lattice
Image:GE715_CMP_25.png|Figure 26: C result for 8x8 lattice
Image:GE715_CMP_26.png|Figure 27: C result for 16x16 lattice
Image:GE715_CMP_27.png|Figure 28: C result for 32x32 lattice
Image:GE715_CMP_28.png|Figure 29: C result for 64x64 lattice
Image:GE715_CMP_29.png|Figure 30: C result for 128x128 lattice
</gallery>

Comparison of heat capacity per spin is given in the figure 31 below.

<gallery mode="packed" heights=500 px>
Image:GE715_CMP_30.png|Figure 31: Heat capacity comparison
</gallery>

== Curie temperature ==

Curie temperature can be deduced from where there is a peak in heat capacity per spin. The graphs of heat capacity vs temperature is shown below (figures 32 and 33) for two different lattice sizes, once using Python script and once using C++ data provided.

<gallery mode="packed" heights=500 px>
Image:GE715_CMP_31.png|Figure 32: 16x16 lattice
Image:GE715_CMP_32.png|Figure 33: 64x64 lattice
</gallery>

=== Polynomial fit ===

The C++ data was also fitted with a polynomial, firstly over the whole range and afterwards just in the peak range. Fitting a polynomial proved to be hard as expected (since data is close to delta-like function), but it became more feasible when smaller range was used. Examples of both fits are shown in figures 34 and 35 below.

<gallery mode="packed" heights=500 px>
Image:GE715_CMP_33.png|Figure 34: Whole region fit
Image:GE715_CMP_34.png|Figure 35: Peak region fit
</gallery>

=== Finding Curie temperature ===

The last step was to find peak heat capacity for different lattice sizes in C++ data. The results are shown in the table below.

{| class="wikitable"
|-
! Lattice size
! Curie temperature
|-
|2x2
|2.48 K
|-
|4x4
|2.40 K
|-
|8x8
|2.32 K
|-
|16x16
|2.25 K
|-
|32x32
|2.30 K
|-
|64x64
|2.25 K
|}

The curie temperature for infinite lattice is known to be calculated from that of a finite lattice using equation <math>T_{C, L} = \frac{A}{L} + T_{C,\infty}</math>, where <math>T_{C, L}</math> is the Curie temperature of an <math>L\times L</math>lattice, <math>T_{C,\infty}</math> is the Curie temperature of an infinite lattice, and <math>A</math> is a constant. This was used to draw a graph shown in figure 36 below.

<gallery mode="packed" heights=500 px>
Image:GE715_CMP_35.png|Figure 36: Curie temperature for infinite lattice
</gallery>

The Curie temperature for infinite lattice was found to be <math>T=2.26 K</math>.

Rep:CMPY3 GE715

2017-11-22T10:44:05Z

Ge715: /* Finding Curie temperature */

== Ising Model ==

The interaction energy in Ising model is defined as <math>E=- \frac{1}{2} J \sum_i^N \sum_{j\ \in\ \mathrm{neighbours}\left(i\right)} s_i s_j</math>. The energy will be the lowest possible when all the spins point in the same direction. In that case <math>s_i s_j=1</math> for all possible <math>i,j</math> pairs. The number of pairs for each <math>i</math> depends on the number of dimensions.

<math>P=2D</math>, where <math>P</math> represents the number of pairs and <math>D</math> the number of dimensions. The lowest possible energy will thus be <math>E=- \frac{1}{2} JNP</math>, where <math>N</math> represents the total number of spins. Thus the lowest possible energy may also be simplified as <math>E=-DNJ</math>, as requested.

The multiplicity of this particular state is 2 (all spins up or all spins down), thus <math>\Omega=2</math>, entropy of the system is <math>S=k_b ln \Omega=9.6*10^{-24} JK^{-1}</math>.

=== Flip of a spin ===

If one of the spins flips, the energy increases. In the given case (<math>D=3, N=1000</math>), the <math>E_i=-3000 J</math>, while the flip of one spin increases the energy by <math>12 J</math> (6 pairs, <math>2 J</math> per pair), so that <math>E_f=-2988 J</math>. Note that <math>J</math> does not represent joules in this case.

Multiplicity has increased to <math>\Omega=2000</math>, thus <math>S=1.05*10^{-22} JK^{-1}</math>.

<math>\Delta S=1.04*10^{-22} JK^{-1}</math>.

=== Magnetisation in the sketch ===

The 1D system shown has the magnetisation of <math>M=1</math>. This is also true for the 2D system shown. For the model system with <math>N=1000, D=3</math> the expected magnetisation at absolute zero would be <math>M=\pm 1000</math>.

== Code for magnetisation and energy ==

Python code was written to calculate the energy and magnetisation of a configuration and tested a few times. The test results are shown in figure 1 below.

<gallery mode="packed" heights=800 px>
Image:GE715_CMP_0.png|Figure 1: Test results
</gallery>

== Accelerating the code ==

The code at the beginning was written using a loop for magnetisation and a double loop for energy. 2000 Monte Carlo steps needed was <math>M=9.32 s \pm 0.10 s</math>. After optimizing magnetisation only this has dropped to <math>M=8.88 s \pm 0.04 s</math>, while after optimizing both magnetisation and energy with using Numpy's sum, roll and multiply functions, the time was further reduced to <math>M=0.0796 s \pm 0.018 s</math>. It was later realised that only one roll in each direction is needed to produce the correct energy, so the final time for 2000 steps was successfully reduced to <math>M=0.418 s \pm 0.005 s</math>.

The whole of the code used is given below.

<pre>
import numpy as np

class IsingLattice:

E = 0.0
E2 = 0.0
M = 0.0
M2 = 0.0
energy_x=0
energy2_x=0
magnetisation_x=0
magnetisation2_x=0

n_cycles = 0

def __init__(self, n_rows, n_cols):
self.n_rows = n_rows
self.n_cols = n_cols
self.lattice = np.random.choice([-1,1], size=(n_rows, n_cols))

def energy(self):
a=np.roll(self.lattice,1,axis=0)
a1=np.multiply(a,self.lattice)
x1=np.sum(a1)
c=np.roll(self.lattice,1,axis=1)
c1=np.multiply(c,self.lattice)
x2=np.sum(c1)
x=x1+x2
energy = -x
return energy

def magnetisation(self):
"Return the total magnetisation of the current lattice configuration."
return np.sum(self.lattice)

def energy_flip(self,random_i,random_j):
"Return the total energy of the current lattice after flip configuration."
self.lattice[random_i][random_j]=-self.lattice[random_i][random_j]
a=np.roll(self.lattice,1,axis=0)
a1=np.multiply(a,self.lattice)
x1=np.sum(a1)
c=np.roll(self.lattice,1,axis=1)
c1=np.multiply(c,self.lattice)
x2=np.sum(c1)
x=x1+x2
energy = -x
return energy

def montecarlostep(self, T):
# complete this function so that it performs a single Monte Carlo step
energy = self.energy()
#the following two lines will select the coordinates of the random spin for you
random_i = np.random.choice(range(0, self.n_rows))
random_j = np.random.choice(range(0, self.n_cols))
#the following line will choose a random number in the range [0,1) for you
random_number = np.random.random()
self.lattice[random_i][random_j] = -self.lattice[random_i][random_j]
energy2 = self.energy()
dE = energy2 - energy

if energy2<energy or random_number<=np.e**(-(energy2-energy)/T):
yyy=True
else:
self.lattice[random_i][random_j] = -self.lattice[random_i][random_j]

self.n_cycles+=1

self.E += self.energy()
self.E2 += self.energy()**2
self.M += self.magnetisation()
self.M2 += self.magnetisation()**2
return self.energy(), self.magnetisation()

def statistics(self):
# complete this function so that it calculates the correct values for the averages of E, E*E (E2), M, M*M (M2), and returns them
avE = self.E / self.n_cycles
avE2 = self.E2 / self.n_cycles
avM = self.M / self.n_cycles
avM2 = self.M / self.n_cycles
return avE, avE2, avM, avM2, self.n_cycles
</pre>

== Temperature effects ==

=== Correction of averaging code ===

A few experiments were performed at different temperatures using ILfinalframe.py, to see how long it usually takes for the system to equilibrate. At low temperatures it usually happened before 1000 steps were performed, thus 1000 Monte Carlo steps was decided to be appropriate cut-off period for calculation of average values of magnetisation and energy. The figures 2-4, showing the experiments performed are shown in the gallery below.

<gallery mode="packed" heights=250 px>
Image:GE715_CMP_1.png|Figure 2: Test result @ 0.5 K
Image:GE715_CMP_2.png|Figure 3: Test result @ 1 K
Image:GE715_CMP_3.png|Figure 4: Test result @ 1.5 K
</gallery>

The definitions of two class functions were changed as shown below to accommodate this change.

<pre>
def montecarlostep(self, T):
# complete this function so that it performs a single Monte Carlo step
energy = self.energy()
#the following two lines will select the coordinates of the random spin for you
random_i = np.random.choice(range(0, self.n_rows))
random_j = np.random.choice(range(0, self.n_cols))
#the following line will choose a random number in the range [0,1) for you
random_number = np.random.random()
self.lattice[random_i][random_j] = -self.lattice[random_i][random_j]
energy2 = self.energy()
dE = energy2 - energy

if energy2<energy or random_number<=np.e**(-(energy2-energy)/T):
yyy=True
else:
self.lattice[random_i][random_j] = -self.lattice[random_i][random_j]

self.n_cycles+=1
if self.n_cycles==1000:
self.E=0
self.E2=0
self.M=0
self.M2=0
self.E += self.energy()
self.E2 += self.energy()**2
self.M += self.magnetisation()
self.M2 += self.magnetisation()**2
return self.energy(), self.magnetisation()

def statistics(self):
# complete this function so that it calculates the correct values for the averages of E, E*E (E2), M, M*M (M2), and returns them
avE = self.E / self.n_cycles
avE2 = self.E2 / self.n_cycles
avM = self.M / self.n_cycles
avM2 = self.M / self.n_cycles
if self.n_cycles>1000:
avE = self.E / (self.n_cycles-1000)
avE2 = self.E2 / (self.n_cycles-1000)
avM = self.M / (self.n_cycles-1000)
avM2 = self.M / (self.n_cycles-1000)
return avE, avE2, avM, avM2, self.n_cycles
</pre>

=== Range of temperatures ===

The experiment was run for 8x8 lattice at a temperature range of 0.25 K-5.00 K, with a point every 0.05 K. The experiments using a different number of total steps were performed and the comparison is shown below in figures 5-7.

<gallery mode="packed" heights=250 px>
Image:GE715_CMP_4.png|Figure 5: Test result for 8x8 lattice with 11000 steps
Image:GE715_CMP_5.png|Figure 6: Test result for 8x8 lattice with 21000 steps
Image:GE715_CMP_6.png|Figure 7: Test result for 8x8 lattice with 101000 steps
</gallery>

There is a small, but noticeable improvement in the smoothness of the curves obtained for energy with using more steps, thus the highest number of steps (101000 steps) was used in the further experiments regardless of the higher need for computational power.

== Lattice size effects ==

With increasing lattice size the energy and magnetisation curves become more realistic, since the long range effects are better accounted for. 7 different lattice sizes were tested and energy and magnetisation for each are shown below in figure 8-14 and 15-21 respectively.

<gallery mode="packed" heights=250 px>
Image:GE715_CMP_7.png|Figure 8: Energy result for 2x2 lattice
Image:GE715_CMP_8.png|Figure 9: Energy result for 4x4 lattice
Image:GE715_CMP_9.png|Figure 10: Energy result for 8x8 lattice
Image:GE715_CMP_10.png|Figure 11: Energy result for 16x16 lattice
Image:GE715_CMP_11.png|Figure 12: Energy result for 32x32 lattice
Image:GE715_CMP_12.png|Figure 13: Energy result for 64x64 lattice
Image:GE715_CMP_13.png|Figure 14: Energy result for 128x128 lattice
Image:GE715_CMP_14.png|Figure 15: Magnetisation result for 2x2 lattice
Image:GE715_CMP_15.png|Figure 16: Magnetisation result for 4x4 lattice
Image:GE715_CMP_16.png|Figure 17: Magnetisation result for 8x8 lattice
Image:GE715_CMP_17.png|Figure 18: Magnetisation result for 16x16 lattice
Image:GE715_CMP_18.png|Figure 19: Magnetisation result for 32x32 lattice
Image:GE715_CMP_19.png|Figure 20: Magnetisation result for 64x64 lattice
Image:GE715_CMP_20.png|Figure 21: Magnetisation result for 128x128 lattice
</gallery>

Overall comparison of energy and magnetisation per spin in each lattice was made and can be seen in figures 22 and 23 below.

<gallery mode="packed" heights=500 px>
Image:GE715_CMP_21.png|Figure 22: Energy comparison
Image:GE715_CMP_22.png|Figure 23: Magnetisation comparison
</gallery>

The code used for producing figure 23 above is shown below. Other codes are analagous. Example for defining ls2 is also shown Others are similar..

<pre>
ls2=loadtxt("2x2(100).dat")

figsize(12,8)
plot(ls2[:,0],ls2[:,3]/4,color="red",label="2x2")
plot(ls4[:,0],ls4[:,3]/16,color="orange",label="4x4")
plot(ls8[:,0],ls8[:,3]/64,color="yellow",label="8x8")
plot(ls16[:,0],ls16[:,3]/256,color="lime",label="16x16")
plot(ls32[:,0],ls32[:,3]/1024,color="blue",label="32x32")
plot(ls64[:,0],ls64[:,3]/4096,color="indigo",label="64x64")
plot(ls128[:,0],ls128[:,3]/16384,color="violet",label="128x128")
title("Magnetisation per spin")
xlabel("Temperature (K)")
ylabel("Magnetisation per spin")
legend()
savefig("Magnetisation.png")
show()
</pre>

Significant changes may be observed from 2x2 lattice up to 16x16 lattice, later the additional long range effects become less noticeable, but new effects are still present. It seems that most of the long range effects are captured by 16x16 lattice, but higher dimensions bring additional improvement.

== Heat Capacity ==

This formula holds true by definition: <math>C = \frac{\partial \left\langle E\right\rangle}{\partial T}</math>.

Variance of energy can be calculated in the following way: <math>\mathrm{Var}[X] = \left\langle E^2\right\rangle - \left\langle E\right\rangle^2</math>.

Since expected energy can also be expressed as partition function: <math>\left\langle E\right\rangle = \frac{\partial}{\partial \beta}lnZ</math>, variance may be expressed as <math>\mathrm{Var}[X] = \frac{\partial ^2}{\partial \beta ^2}lnZ = - \frac{\partial \left\langle E\right\rangle}{\partial \beta}</math>.

We can rearrange the equations given, because <math>T = \frac{1}{k_BT}</math>, so that we obtain <math>C = \frac{\mathrm{Var}[E]}{k_B T^2}</math>.

=== Heat capacity of the model ===

Heat capacities of the models were calculated using the script given below (example for 2x2 lattice).

<pre>
figsize(6,4)
plot(ls2[:,0],(ls2[:,2]-ls2[:,1]**2)/ls2[:,0],color="red",label="2x2")
xlabel("Temperature (K)")
ylabel("C")
savefig("C2.png")
show()
</pre>

Graphs of heat capacities are shown in the gallery below (figures 24-30). It can be seen that the model breaks down at low temperatures for bigger lattices, but it still gives a peak heat capacity.

<gallery mode="packed" heights=250 px>
Image:GE715_CMP_23.png|Figure 24: C result for 2x2 lattice
Image:GE715_CMP_24.png|Figure 25: C result for 4x4 lattice
Image:GE715_CMP_25.png|Figure 26: C result for 8x8 lattice
Image:GE715_CMP_26.png|Figure 27: C result for 16x16 lattice
Image:GE715_CMP_27.png|Figure 28: C result for 32x32 lattice
Image:GE715_CMP_28.png|Figure 29: C result for 64x64 lattice
Image:GE715_CMP_29.png|Figure 30: C result for 128x128 lattice
</gallery>

Comparison of heat capacity per spin is given in the figure 31 below.

<gallery mode="packed" heights=500 px>
Image:GE715_CMP_30.png|Figure 31: Heat capacity comparison
</gallery>

== Curie temperature ==

Curie temperature can be deduced from where there is a peak in heat capacity per spin. The graphs of heat capacity vs temperature is shown below (figures 32 and 33) for two different lattice sizes, once using Python script and once using C++ data provided.

<gallery mode="packed" heights=500 px>
Image:GE715_CMP_31.png|Figure 32: 16x16 lattice
Image:GE715_CMP_32.png|Figure 33: 64x64 lattice
</gallery>

=== Polynomial fit ===

The C++ data was also fitted with a polynomial, firstly over the whole range and afterwards just in the peak range. Fitting a polynomial proved to be hard as expected (since data is close to delta-like function), but it became more feasible when smaller range was used. Examples of both fits are shown in figures 34 and 35 below.

<gallery mode="packed" heights=500 px>
Image:GE715_CMP_33.png|Figure 34: Whole region fit
Image:GE715_CMP_34.png|Figure 35: Peak region fit
</gallery>

=== Finding Curie temperature ===

The last step was to find peak heat capacity for different lattice sizes in C++ data. The results are shown in the table below.

{| class="wikitable"
|-
! Lattice size
! Curie temperature
|-
|2x2
|2.48 K
|-
|4x4
|2.40 K
|-
|8x8
|2.32 K
|-
|16x16
|2.25 K
|-
|32x32
|2.30 K
|-
|64x64
|2.25 K
|}

The curie temperature for infinite lattice is known to be calculated from that of a finite lattice using equation <math>T_{C, L} = \frac{A}{L} + T_{C,\infty}</math>, where <math>T_{C, L}</math> is the Curie temperature of an <math>L\times L</math>lattice, <math>T_{C,\infty}</math> is the Curie temperature of an infinite lattice, and <math>A</math> is a constant. This was used to draw a graph shown in figure X below.

The Curie temperature for infinite lattice was found to be <math>T=2.26 K</math>.

Rep:CMPY3 GE715

2017-11-21T19:13:13Z

Ge715: /* Finding Curie temperature */ Fail

Rep:CMPY3 GE715

2017-11-21T19:10:03Z

Ge715: /* Finding Curie temperature */ Formatting mistake

== Ising Model ==

The interaction energy in Ising model is defined as <math>E=- \frac{1}{2} J \sum_i^N \sum_{j\ \in\ \mathrm{neighbours}\left(i\right)} s_i s_j</math>. The energy will be the lowest possible when all the spins point in the same direction. In that case <math>s_i s_j=1</math> for all possible <math>i,j</math> pairs. The number of pairs for each <math>i</math> depends on the number of dimensions.

<math>P=2D</math>, where <math>P</math> represents the number of pairs and <math>D</math> the number of dimensions. The lowest possible energy will thus be <math>E=- \frac{1}{2} JNP</math>, where <math>N</math> represents the total number of spins. Thus the lowest possible energy may also be simplified as <math>E=-DNJ</math>, as requested.

The multiplicity of this particular state is 2 (all spins up or all spins down), thus <math>\Omega=2</math>, entropy of the system is <math>S=k_b ln \Omega=9.6*10^{-24} JK^{-1}</math>.

=== Flip of a spin ===

If one of the spins flips, the energy increases. In the given case (<math>D=3, N=1000</math>), the <math>E_i=-3000 J</math>, while the flip of one spin increases the energy by <math>12 J</math> (6 pairs, <math>2 J</math> per pair), so that <math>E_f=-2988 J</math>. Note that <math>J</math> does not represent joules in this case.

Multiplicity has increased to <math>\Omega=2000</math>, thus <math>S=1.05*10^{-22} JK^{-1}</math>.

<math>\Delta S=1.04*10^{-22} JK^{-1}</math>.

=== Magnetisation in the sketch ===

The 1D system shown has the magnetisation of <math>M=1</math>. This is also true for the 2D system shown. For the model system with <math>N=1000, D=3</math> the expected magnetisation at absolute zero would be <math>M=\pm 1000</math>.

== Code for magnetisation and energy ==

Python code was written to calculate the energy and magnetisation of a configuration and tested a few times. The test results are shown in figure 1 below.

<gallery mode="packed" heights=800 px>
Image:GE715_CMP_0.png|Figure 1: Test results
</gallery>

== Accelerating the code ==

The code at the beginning was written using a loop for magnetisation and a double loop for energy. 2000 Monte Carlo steps needed was <math>M=9.32 s \pm 0.10 s</math>. After optimizing magnetisation only this has dropped to <math>M=8.88 s \pm 0.04 s</math>, while after optimizing both magnetisation and energy with using Numpy's sum, roll and multiply functions, the time was further reduced to <math>M=0.0796 s \pm 0.018 s</math>. It was later realised that only one roll in each direction is needed to produce the correct energy, so the final time for 2000 steps was successfully reduced to <math>M=0.418 s \pm 0.005 s</math>.

The whole of the code used is given below.

<pre>
import numpy as np

class IsingLattice:

E = 0.0
E2 = 0.0
M = 0.0
M2 = 0.0
energy_x=0
energy2_x=0
magnetisation_x=0
magnetisation2_x=0

n_cycles = 0

def __init__(self, n_rows, n_cols):
self.n_rows = n_rows
self.n_cols = n_cols
self.lattice = np.random.choice([-1,1], size=(n_rows, n_cols))

def energy(self):
a=np.roll(self.lattice,1,axis=0)
a1=np.multiply(a,self.lattice)
x1=np.sum(a1)
c=np.roll(self.lattice,1,axis=1)
c1=np.multiply(c,self.lattice)
x2=np.sum(c1)
x=x1+x2
energy = -x
return energy

def magnetisation(self):
"Return the total magnetisation of the current lattice configuration."
return np.sum(self.lattice)

def energy_flip(self,random_i,random_j):
"Return the total energy of the current lattice after flip configuration."
self.lattice[random_i][random_j]=-self.lattice[random_i][random_j]
a=np.roll(self.lattice,1,axis=0)
a1=np.multiply(a,self.lattice)
x1=np.sum(a1)
c=np.roll(self.lattice,1,axis=1)
c1=np.multiply(c,self.lattice)
x2=np.sum(c1)
x=x1+x2
energy = -x
return energy

def montecarlostep(self, T):
# complete this function so that it performs a single Monte Carlo step
energy = self.energy()
#the following two lines will select the coordinates of the random spin for you
random_i = np.random.choice(range(0, self.n_rows))
random_j = np.random.choice(range(0, self.n_cols))
#the following line will choose a random number in the range [0,1) for you
random_number = np.random.random()
self.lattice[random_i][random_j] = -self.lattice[random_i][random_j]
energy2 = self.energy()
dE = energy2 - energy

if energy2<energy or random_number<=np.e**(-(energy2-energy)/T):
yyy=True
else:
self.lattice[random_i][random_j] = -self.lattice[random_i][random_j]

self.n_cycles+=1

self.E += self.energy()
self.E2 += self.energy()**2
self.M += self.magnetisation()
self.M2 += self.magnetisation()**2
return self.energy(), self.magnetisation()

def statistics(self):
# complete this function so that it calculates the correct values for the averages of E, E*E (E2), M, M*M (M2), and returns them
avE = self.E / self.n_cycles
avE2 = self.E2 / self.n_cycles
avM = self.M / self.n_cycles
avM2 = self.M / self.n_cycles
return avE, avE2, avM, avM2, self.n_cycles
</pre>

== Temperature effects ==

=== Correction of averaging code ===

A few experiments were performed at different temperatures using ILfinalframe.py, to see how long it usually takes for the system to equilibrate. At low temperatures it usually happened before 1000 steps were performed, thus 1000 Monte Carlo steps was decided to be appropriate cut-off period for calculation of average values of magnetisation and energy. The figures 2-4, showing the experiments performed are shown in the gallery below.

<gallery mode="packed" heights=250 px>
Image:GE715_CMP_1.png|Figure 2: Test result @ 0.5 K
Image:GE715_CMP_2.png|Figure 3: Test result @ 1 K
Image:GE715_CMP_3.png|Figure 4: Test result @ 1.5 K
</gallery>

The definitions of two class functions were changed as shown below to accommodate this change.

<pre>
def montecarlostep(self, T):
# complete this function so that it performs a single Monte Carlo step
energy = self.energy()
#the following two lines will select the coordinates of the random spin for you
random_i = np.random.choice(range(0, self.n_rows))
random_j = np.random.choice(range(0, self.n_cols))
#the following line will choose a random number in the range [0,1) for you
random_number = np.random.random()
self.lattice[random_i][random_j] = -self.lattice[random_i][random_j]
energy2 = self.energy()
dE = energy2 - energy

if energy2<energy or random_number<=np.e**(-(energy2-energy)/T):
yyy=True
else:
self.lattice[random_i][random_j] = -self.lattice[random_i][random_j]

self.n_cycles+=1
if self.n_cycles==1000:
self.E=0
self.E2=0
self.M=0
self.M2=0
self.E += self.energy()
self.E2 += self.energy()**2
self.M += self.magnetisation()
self.M2 += self.magnetisation()**2
return self.energy(), self.magnetisation()

def statistics(self):
# complete this function so that it calculates the correct values for the averages of E, E*E (E2), M, M*M (M2), and returns them
avE = self.E / self.n_cycles
avE2 = self.E2 / self.n_cycles
avM = self.M / self.n_cycles
avM2 = self.M / self.n_cycles
if self.n_cycles>1000:
avE = self.E / (self.n_cycles-1000)
avE2 = self.E2 / (self.n_cycles-1000)
avM = self.M / (self.n_cycles-1000)
avM2 = self.M / (self.n_cycles-1000)
return avE, avE2, avM, avM2, self.n_cycles
</pre>

=== Range of temperatures ===

The experiment was run for 8x8 lattice at a temperature range of 0.25 K-5.00 K, with a point every 0.05 K. The experiments using a different number of total steps were performed and the comparison is shown below in figures 5-7.

<gallery mode="packed" heights=250 px>
Image:GE715_CMP_4.png|Figure 5: Test result for 8x8 lattice with 11000 steps
Image:GE715_CMP_5.png|Figure 6: Test result for 8x8 lattice with 21000 steps
Image:GE715_CMP_6.png|Figure 7: Test result for 8x8 lattice with 101000 steps
</gallery>

There is a small, but noticeable improvement in the smoothness of the curves obtained for energy with using more steps, thus the highest number of steps (101000 steps) was used in the further experiments regardless of the higher need for computational power.

== Lattice size effects ==

With increasing lattice size the energy and magnetisation curves become more realistic, since the long range effects are better accounted for. 7 different lattice sizes were tested and energy and magnetisation for each are shown below in figure 8-14 and 15-21 respectively.

<gallery mode="packed" heights=250 px>
Image:GE715_CMP_7.png|Figure 8: Energy result for 2x2 lattice
Image:GE715_CMP_8.png|Figure 9: Energy result for 4x4 lattice
Image:GE715_CMP_9.png|Figure 10: Energy result for 8x8 lattice
Image:GE715_CMP_10.png|Figure 11: Energy result for 16x16 lattice
Image:GE715_CMP_11.png|Figure 12: Energy result for 32x32 lattice
Image:GE715_CMP_12.png|Figure 13: Energy result for 64x64 lattice
Image:GE715_CMP_13.png|Figure 14: Energy result for 128x128 lattice
Image:GE715_CMP_14.png|Figure 15: Magnetisation result for 2x2 lattice
Image:GE715_CMP_15.png|Figure 16: Magnetisation result for 4x4 lattice
Image:GE715_CMP_16.png|Figure 17: Magnetisation result for 8x8 lattice
Image:GE715_CMP_17.png|Figure 18: Magnetisation result for 16x16 lattice
Image:GE715_CMP_18.png|Figure 19: Magnetisation result for 32x32 lattice
Image:GE715_CMP_19.png|Figure 20: Magnetisation result for 64x64 lattice
Image:GE715_CMP_20.png|Figure 21: Magnetisation result for 128x128 lattice
</gallery>

Overall comparison of energy and magnetisation per spin in each lattice was made and can be seen in figures 22 and 23 below.

<gallery mode="packed" heights=500 px>
Image:GE715_CMP_21.png|Figure 22: Energy comparison
Image:GE715_CMP_22.png|Figure 23: Magnetisation comparison
</gallery>

The code used for producing figure 23 above is shown below. Other codes are analagous. Example for defining ls2 is also shown Others are similar..

<pre>
ls2=loadtxt("2x2(100).dat")

figsize(12,8)
plot(ls2[:,0],ls2[:,3]/4,color="red",label="2x2")
plot(ls4[:,0],ls4[:,3]/16,color="orange",label="4x4")
plot(ls8[:,0],ls8[:,3]/64,color="yellow",label="8x8")
plot(ls16[:,0],ls16[:,3]/256,color="lime",label="16x16")
plot(ls32[:,0],ls32[:,3]/1024,color="blue",label="32x32")
plot(ls64[:,0],ls64[:,3]/4096,color="indigo",label="64x64")
plot(ls128[:,0],ls128[:,3]/16384,color="violet",label="128x128")
title("Magnetisation per spin")
xlabel("Temperature (K)")
ylabel("Magnetisation per spin")
legend()
savefig("Magnetisation.png")
show()
</pre>

Significant changes may be observed from 2x2 lattice up to 16x16 lattice, later the additional long range effects become less noticeable, but new effects are still present. It seems that most of the long range effects are captured by 16x16 lattice, but higher dimensions bring additional improvement.

== Heat Capacity ==

This formula holds true by definition: <math>C = \frac{\partial \left\langle E\right\rangle}{\partial T}</math>.

Variance of energy can be calculated in the following way: <math>\mathrm{Var}[X] = \left\langle E^2\right\rangle - \left\langle E\right\rangle^2</math>.

Since expected energy can also be expressed as partition function: <math>\left\langle E\right\rangle = \frac{\partial}{\partial \beta}lnZ</math>, variance may be expressed as <math>\mathrm{Var}[X] = \frac{\partial ^2}{\partial \beta ^2}lnZ = - \frac{\partial \left\langle E\right\rangle}{\partial \beta}</math>.

We can rearrange the equations given, because <math>T = \frac{1}{k_BT}</math>, so that we obtain <math>C = \frac{\mathrm{Var}[E]}{k_B T^2}</math>.

=== Heat capacity of the model ===

Heat capacities of the models were calculated using the script given below (example for 2x2 lattice).

<pre>
figsize(6,4)
plot(ls2[:,0],(ls2[:,2]-ls2[:,1]**2)/ls2[:,0],color="red",label="2x2")
xlabel("Temperature (K)")
ylabel("C")
savefig("C2.png")
show()
</pre>

Graphs of heat capacities are shown in the gallery below (figures 24-30). It can be seen that the model breaks down at low temperatures for bigger lattices, but it still gives a peak heat capacity.

<gallery mode="packed" heights=250 px>
Image:GE715_CMP_23.png|Figure 24: C result for 2x2 lattice
Image:GE715_CMP_24.png|Figure 25: C result for 4x4 lattice
Image:GE715_CMP_25.png|Figure 26: C result for 8x8 lattice
Image:GE715_CMP_26.png|Figure 27: C result for 16x16 lattice
Image:GE715_CMP_27.png|Figure 28: C result for 32x32 lattice
Image:GE715_CMP_28.png|Figure 29: C result for 64x64 lattice
Image:GE715_CMP_29.png|Figure 30: C result for 128x128 lattice
</gallery>

Comparison of heat capacity per spin is given in the figure 31 below.

<gallery mode="packed" heights=500 px>
Image:GE715_CMP_30.png|Figure 31: Heat capacity comparison
</gallery>

== Curie temperature ==

Curie temperature can be deduced from where there is a peak in heat capacity per spin. The graphs of heat capacity vs temperature is shown below (figures 32 and 33) for two different lattice sizes, once using Python script and once using C++ data provided.

<gallery mode="packed" heights=500 px>
Image:GE715_CMP_31.png|Figure 32: 16x16 lattice
Image:GE715_CMP_32.png|Figure 33: 64x64 lattice
</gallery>

=== Polynomial fit ===

The C++ data was also fitted with a polynomial, firstly over the whole range and afterwards just in the peak range. Fitting a polynomial proved to be hard as expected (since data is close to delta-like function), but it became more feasible when smaller range was used. Examples of both fits are shown in figures 34 and 35 below.

<gallery mode="packed" heights=500 px>
Image:GE715_CMP_33.png|Figure 34: Whole region fit
Image:GE715_CMP_34.png|Figure 35: Peak region fit
</gallery>

=== Finding Curie temperature ===

The last step was to find peak heat capacity for different lattice sizes in C++ data. The results are shown in the table below.

{| class="wikitable"
|-
! Lattice size
! Curie temperature
|-
|2x2
|2.48 K
|-
|4x4
|2.40 K
|-
|8x8
|2.32 K
|-
|16x16
|2.25 K
|-
|32x32
|2.30 K
|-
|64x64
|2.25 K
|}

It becomes clear that for larger lattice sizes the calculated Curie temperature will slowly converge, thus we can conclude that for infinite lattice it will be in the region of 2.2-2.3 K.

Rep:CMPY3 GE715

2017-11-21T19:09:28Z

Ge715: Part 8

== Ising Model ==

The interaction energy in Ising model is defined as <math>E=- \frac{1}{2} J \sum_i^N \sum_{j\ \in\ \mathrm{neighbours}\left(i\right)} s_i s_j</math>. The energy will be the lowest possible when all the spins point in the same direction. In that case <math>s_i s_j=1</math> for all possible <math>i,j</math> pairs. The number of pairs for each <math>i</math> depends on the number of dimensions.

<math>P=2D</math>, where <math>P</math> represents the number of pairs and <math>D</math> the number of dimensions. The lowest possible energy will thus be <math>E=- \frac{1}{2} JNP</math>, where <math>N</math> represents the total number of spins. Thus the lowest possible energy may also be simplified as <math>E=-DNJ</math>, as requested.

The multiplicity of this particular state is 2 (all spins up or all spins down), thus <math>\Omega=2</math>, entropy of the system is <math>S=k_b ln \Omega=9.6*10^{-24} JK^{-1}</math>.

=== Flip of a spin ===

If one of the spins flips, the energy increases. In the given case (<math>D=3, N=1000</math>), the <math>E_i=-3000 J</math>, while the flip of one spin increases the energy by <math>12 J</math> (6 pairs, <math>2 J</math> per pair), so that <math>E_f=-2988 J</math>. Note that <math>J</math> does not represent joules in this case.

Multiplicity has increased to <math>\Omega=2000</math>, thus <math>S=1.05*10^{-22} JK^{-1}</math>.

<math>\Delta S=1.04*10^{-22} JK^{-1}</math>.

=== Magnetisation in the sketch ===

The 1D system shown has the magnetisation of <math>M=1</math>. This is also true for the 2D system shown. For the model system with <math>N=1000, D=3</math> the expected magnetisation at absolute zero would be <math>M=\pm 1000</math>.

== Code for magnetisation and energy ==

Python code was written to calculate the energy and magnetisation of a configuration and tested a few times. The test results are shown in figure 1 below.

<gallery mode="packed" heights=800 px>
Image:GE715_CMP_0.png|Figure 1: Test results
</gallery>

== Accelerating the code ==

The code at the beginning was written using a loop for magnetisation and a double loop for energy. 2000 Monte Carlo steps needed was <math>M=9.32 s \pm 0.10 s</math>. After optimizing magnetisation only this has dropped to <math>M=8.88 s \pm 0.04 s</math>, while after optimizing both magnetisation and energy with using Numpy's sum, roll and multiply functions, the time was further reduced to <math>M=0.0796 s \pm 0.018 s</math>. It was later realised that only one roll in each direction is needed to produce the correct energy, so the final time for 2000 steps was successfully reduced to <math>M=0.418 s \pm 0.005 s</math>.

The whole of the code used is given below.

<pre>
import numpy as np

class IsingLattice:

E = 0.0
E2 = 0.0
M = 0.0
M2 = 0.0
energy_x=0
energy2_x=0
magnetisation_x=0
magnetisation2_x=0

n_cycles = 0

def __init__(self, n_rows, n_cols):
self.n_rows = n_rows
self.n_cols = n_cols
self.lattice = np.random.choice([-1,1], size=(n_rows, n_cols))

def energy(self):
a=np.roll(self.lattice,1,axis=0)
a1=np.multiply(a,self.lattice)
x1=np.sum(a1)
c=np.roll(self.lattice,1,axis=1)
c1=np.multiply(c,self.lattice)
x2=np.sum(c1)
x=x1+x2
energy = -x
return energy

def magnetisation(self):
"Return the total magnetisation of the current lattice configuration."
return np.sum(self.lattice)

def energy_flip(self,random_i,random_j):
"Return the total energy of the current lattice after flip configuration."
self.lattice[random_i][random_j]=-self.lattice[random_i][random_j]
a=np.roll(self.lattice,1,axis=0)
a1=np.multiply(a,self.lattice)
x1=np.sum(a1)
c=np.roll(self.lattice,1,axis=1)
c1=np.multiply(c,self.lattice)
x2=np.sum(c1)
x=x1+x2
energy = -x
return energy

def montecarlostep(self, T):
# complete this function so that it performs a single Monte Carlo step
energy = self.energy()
#the following two lines will select the coordinates of the random spin for you
random_i = np.random.choice(range(0, self.n_rows))
random_j = np.random.choice(range(0, self.n_cols))
#the following line will choose a random number in the range [0,1) for you
random_number = np.random.random()
self.lattice[random_i][random_j] = -self.lattice[random_i][random_j]
energy2 = self.energy()
dE = energy2 - energy

if energy2<energy or random_number<=np.e**(-(energy2-energy)/T):
yyy=True
else:
self.lattice[random_i][random_j] = -self.lattice[random_i][random_j]

self.n_cycles+=1

self.E += self.energy()
self.E2 += self.energy()**2
self.M += self.magnetisation()
self.M2 += self.magnetisation()**2
return self.energy(), self.magnetisation()

def statistics(self):
# complete this function so that it calculates the correct values for the averages of E, E*E (E2), M, M*M (M2), and returns them
avE = self.E / self.n_cycles
avE2 = self.E2 / self.n_cycles
avM = self.M / self.n_cycles
avM2 = self.M / self.n_cycles
return avE, avE2, avM, avM2, self.n_cycles
</pre>

== Temperature effects ==

=== Correction of averaging code ===

A few experiments were performed at different temperatures using ILfinalframe.py, to see how long it usually takes for the system to equilibrate. At low temperatures it usually happened before 1000 steps were performed, thus 1000 Monte Carlo steps was decided to be appropriate cut-off period for calculation of average values of magnetisation and energy. The figures 2-4, showing the experiments performed are shown in the gallery below.

<gallery mode="packed" heights=250 px>
Image:GE715_CMP_1.png|Figure 2: Test result @ 0.5 K
Image:GE715_CMP_2.png|Figure 3: Test result @ 1 K
Image:GE715_CMP_3.png|Figure 4: Test result @ 1.5 K
</gallery>

The definitions of two class functions were changed as shown below to accommodate this change.

<pre>
def montecarlostep(self, T):
# complete this function so that it performs a single Monte Carlo step
energy = self.energy()
#the following two lines will select the coordinates of the random spin for you
random_i = np.random.choice(range(0, self.n_rows))
random_j = np.random.choice(range(0, self.n_cols))
#the following line will choose a random number in the range [0,1) for you
random_number = np.random.random()
self.lattice[random_i][random_j] = -self.lattice[random_i][random_j]
energy2 = self.energy()
dE = energy2 - energy

if energy2<energy or random_number<=np.e**(-(energy2-energy)/T):
yyy=True
else:
self.lattice[random_i][random_j] = -self.lattice[random_i][random_j]

self.n_cycles+=1
if self.n_cycles==1000:
self.E=0
self.E2=0
self.M=0
self.M2=0
self.E += self.energy()
self.E2 += self.energy()**2
self.M += self.magnetisation()
self.M2 += self.magnetisation()**2
return self.energy(), self.magnetisation()

def statistics(self):
# complete this function so that it calculates the correct values for the averages of E, E*E (E2), M, M*M (M2), and returns them
avE = self.E / self.n_cycles
avE2 = self.E2 / self.n_cycles
avM = self.M / self.n_cycles
avM2 = self.M / self.n_cycles
if self.n_cycles>1000:
avE = self.E / (self.n_cycles-1000)
avE2 = self.E2 / (self.n_cycles-1000)
avM = self.M / (self.n_cycles-1000)
avM2 = self.M / (self.n_cycles-1000)
return avE, avE2, avM, avM2, self.n_cycles
</pre>

=== Range of temperatures ===

The experiment was run for 8x8 lattice at a temperature range of 0.25 K-5.00 K, with a point every 0.05 K. The experiments using a different number of total steps were performed and the comparison is shown below in figures 5-7.

<gallery mode="packed" heights=250 px>
Image:GE715_CMP_4.png|Figure 5: Test result for 8x8 lattice with 11000 steps
Image:GE715_CMP_5.png|Figure 6: Test result for 8x8 lattice with 21000 steps
Image:GE715_CMP_6.png|Figure 7: Test result for 8x8 lattice with 101000 steps
</gallery>

There is a small, but noticeable improvement in the smoothness of the curves obtained for energy with using more steps, thus the highest number of steps (101000 steps) was used in the further experiments regardless of the higher need for computational power.

== Lattice size effects ==

With increasing lattice size the energy and magnetisation curves become more realistic, since the long range effects are better accounted for. 7 different lattice sizes were tested and energy and magnetisation for each are shown below in figure 8-14 and 15-21 respectively.

<gallery mode="packed" heights=250 px>
Image:GE715_CMP_7.png|Figure 8: Energy result for 2x2 lattice
Image:GE715_CMP_8.png|Figure 9: Energy result for 4x4 lattice
Image:GE715_CMP_9.png|Figure 10: Energy result for 8x8 lattice
Image:GE715_CMP_10.png|Figure 11: Energy result for 16x16 lattice
Image:GE715_CMP_11.png|Figure 12: Energy result for 32x32 lattice
Image:GE715_CMP_12.png|Figure 13: Energy result for 64x64 lattice
Image:GE715_CMP_13.png|Figure 14: Energy result for 128x128 lattice
Image:GE715_CMP_14.png|Figure 15: Magnetisation result for 2x2 lattice
Image:GE715_CMP_15.png|Figure 16: Magnetisation result for 4x4 lattice
Image:GE715_CMP_16.png|Figure 17: Magnetisation result for 8x8 lattice
Image:GE715_CMP_17.png|Figure 18: Magnetisation result for 16x16 lattice
Image:GE715_CMP_18.png|Figure 19: Magnetisation result for 32x32 lattice
Image:GE715_CMP_19.png|Figure 20: Magnetisation result for 64x64 lattice
Image:GE715_CMP_20.png|Figure 21: Magnetisation result for 128x128 lattice
</gallery>

Overall comparison of energy and magnetisation per spin in each lattice was made and can be seen in figures 22 and 23 below.

<gallery mode="packed" heights=500 px>
Image:GE715_CMP_21.png|Figure 22: Energy comparison
Image:GE715_CMP_22.png|Figure 23: Magnetisation comparison
</gallery>

The code used for producing figure 23 above is shown below. Other codes are analagous. Example for defining ls2 is also shown Others are similar..

<pre>
ls2=loadtxt("2x2(100).dat")

figsize(12,8)
plot(ls2[:,0],ls2[:,3]/4,color="red",label="2x2")
plot(ls4[:,0],ls4[:,3]/16,color="orange",label="4x4")
plot(ls8[:,0],ls8[:,3]/64,color="yellow",label="8x8")
plot(ls16[:,0],ls16[:,3]/256,color="lime",label="16x16")
plot(ls32[:,0],ls32[:,3]/1024,color="blue",label="32x32")
plot(ls64[:,0],ls64[:,3]/4096,color="indigo",label="64x64")
plot(ls128[:,0],ls128[:,3]/16384,color="violet",label="128x128")
title("Magnetisation per spin")
xlabel("Temperature (K)")
ylabel("Magnetisation per spin")
legend()
savefig("Magnetisation.png")
show()
</pre>

Significant changes may be observed from 2x2 lattice up to 16x16 lattice, later the additional long range effects become less noticeable, but new effects are still present. It seems that most of the long range effects are captured by 16x16 lattice, but higher dimensions bring additional improvement.

== Heat Capacity ==

This formula holds true by definition: <math>C = \frac{\partial \left\langle E\right\rangle}{\partial T}</math>.

Variance of energy can be calculated in the following way: <math>\mathrm{Var}[X] = \left\langle E^2\right\rangle - \left\langle E\right\rangle^2</math>.

Since expected energy can also be expressed as partition function: <math>\left\langle E\right\rangle = \frac{\partial}{\partial \beta}lnZ</math>, variance may be expressed as <math>\mathrm{Var}[X] = \frac{\partial ^2}{\partial \beta ^2}lnZ = - \frac{\partial \left\langle E\right\rangle}{\partial \beta}</math>.

We can rearrange the equations given, because <math>T = \frac{1}{k_BT}</math>, so that we obtain <math>C = \frac{\mathrm{Var}[E]}{k_B T^2}</math>.

=== Heat capacity of the model ===

Heat capacities of the models were calculated using the script given below (example for 2x2 lattice).

<pre>
figsize(6,4)
plot(ls2[:,0],(ls2[:,2]-ls2[:,1]**2)/ls2[:,0],color="red",label="2x2")
xlabel("Temperature (K)")
ylabel("C")
savefig("C2.png")
show()
</pre>

Graphs of heat capacities are shown in the gallery below (figures 24-30). It can be seen that the model breaks down at low temperatures for bigger lattices, but it still gives a peak heat capacity.

<gallery mode="packed" heights=250 px>
Image:GE715_CMP_23.png|Figure 24: C result for 2x2 lattice
Image:GE715_CMP_24.png|Figure 25: C result for 4x4 lattice
Image:GE715_CMP_25.png|Figure 26: C result for 8x8 lattice
Image:GE715_CMP_26.png|Figure 27: C result for 16x16 lattice
Image:GE715_CMP_27.png|Figure 28: C result for 32x32 lattice
Image:GE715_CMP_28.png|Figure 29: C result for 64x64 lattice
Image:GE715_CMP_29.png|Figure 30: C result for 128x128 lattice
</gallery>

Comparison of heat capacity per spin is given in the figure 31 below.

<gallery mode="packed" heights=500 px>
Image:GE715_CMP_30.png|Figure 31: Heat capacity comparison
</gallery>

== Curie temperature ==

Curie temperature can be deduced from where there is a peak in heat capacity per spin. The graphs of heat capacity vs temperature is shown below (figures 32 and 33) for two different lattice sizes, once using Python script and once using C++ data provided.

<gallery mode="packed" heights=500 px>
Image:GE715_CMP_31.png|Figure 32: 16x16 lattice
Image:GE715_CMP_32.png|Figure 33: 64x64 lattice
</gallery>

=== Polynomial fit ===

The C++ data was also fitted with a polynomial, firstly over the whole range and afterwards just in the peak range. Fitting a polynomial proved to be hard as expected (since data is close to delta-like function), but it became more feasible when smaller range was used. Examples of both fits are shown in figures 34 and 35 below.

<gallery mode="packed" heights=500 px>
Image:GE715_CMP_33.png|Figure 34: Whole region fit
Image:GE715_CMP_34.png|Figure 35: Peak region fit
</gallery>

=== Finding Curie temperature ===

The last step was to find peak heat capacity for different lattice sizes in C++ data. The results are shown in the table below.

{| class="wikitable"
|-
! Lattice size
! Curie temperature
|-
|2x2
|2.48 K
|-
|4x4
|2.40 K
|-
|8x8
|2.32 K
|-
|16x16
|2.25 K
|-
|32x32
|2.30 K
|-
|64x64
|2.25 K
}

It becomes clear that for larger lattice sizes the calculated Curie temperature will slowly converge, thus we can conclude that for infinite lattice it will be in the region of 2.2-2.3 K.

File:GE715 CMP 34.png

2017-11-21T18:52:49Z

Ge715:

File:GE715 CMP 33.png

2017-11-21T18:52:32Z

Ge715:

File:GE715 CMP 32.png

2017-11-21T18:52:12Z

Ge715:

File:GE715 CMP 31.png

2017-11-21T18:51:54Z

Ge715:

File:GE715 CMP 30.png

2017-11-21T18:05:53Z

Ge715:

File:GE715 CMP 29.png

2017-11-21T18:05:37Z

Ge715:

File:GE715 CMP 28.png

2017-11-21T18:05:19Z

Ge715:

File:GE715 CMP 27.png

2017-11-21T18:04:00Z

Ge715: