diff --git a/IsingReport.pdf b/IsingReport.pdf new file mode 100755 index 0000000..0b07829 Binary files /dev/null and b/IsingReport.pdf differ diff --git a/create_mvsSteps_Plot.p b/create_mvsSteps_Plot.p new file mode 100755 index 0000000..32f76cf --- /dev/null +++ b/create_mvsSteps_Plot.p @@ -0,0 +1,25 @@ +reset + + +set xlabel "Millions of Steps" +set ylabel "Magnetization, m" + +set key right top +set key outside + +plot "Results1.0" using 1:2 with line title "kB T/J=1.0" +replot "Results2.0" using 1:2 with line title "kB T/J=2.0" +replot "Results2.1" using 1:2 with line title "kB T/J=2.1" +replot "Results2.2" using 1:2 with line title "kB T/J=2.2" +replot "Results2.3" using 1:2 with line title "kB T/J=2.3" +replot "Results2.4" using 1:2 with line title "kB T/J=2.4" +replot "Results2.5" using 1:2 with line title "kB T/J=2.5" +replot "Results2.6" using 1:2 with line title "kB T/J=2.6" +replot "Results2.7" using 1:2 with line title "kB T/J=2.7" +replot "Results2.8" using 1:2 with line title "kB T/J=2.8" +replot "Results2.9" using 1:2 with line title "kB T/J=2.9" +replot "Results3.0" using 1:2 with line title "kB T/J=3.0" +replot "Results4.0" using 1:2 with line title "kB T/J=4.0" +set terminal png size 600,300 enhanced font ",20" +set output "m_vs_Steps_plot.png" +replot diff --git a/create_mvsT_Plot.p b/create_mvsT_Plot.p new file mode 100755 index 0000000..87ce1a7 --- /dev/null +++ b/create_mvsT_Plot.p @@ -0,0 +1,22 @@ +reset + +set xlabel "k_B*T/J" +set ylabel "Magnetization, m" + +unset key + +plot "avg_m_values" lt rgb "black" with linespoints, "avg_m_values" using 1:2:3 with yerrorbars +set terminal png size 400,300 enhanced font ",20" + +set term png +set output "m_vs_T_plot.png" +replot + + + + + + + + + diff --git a/isingCode.f95 b/isingCode.f95 new file mode 100755 index 0000000..7f245bd --- /dev/null +++ b/isingCode.f95 @@ -0,0 +1,503 @@ +program isingtest + +IMPLICIT NONE + +!Definition of all varables, arrays and strings + +integer :: i, j, flipSpins, k, Tstep, everyHundred +integer :: size, i_swap_value, j_swap_value +integer :: avgCount, mCount +integer, parameter :: out_unit=20 +integer, parameter :: out_unit2=10 +integer, dimension(100,100) :: spins + +real, dimension(3,7000000) :: m_T +real, dimension(20) :: m_bins +real :: e_new, e_old, delta_e +real :: m, rand, kB, beta, g, T +real :: avgM, m_avg, stdev, diff +real :: diffSq, diffSqSum + +character(len=10) :: out_file +character(len=7) :: out_file_base +character(len=3) :: out_file_no + +! Dimension of square array +size=100 + +! Boltzmann Constant +kB=1 + +out_file_base="Results" + +open (unit=out_unit2,file="avg_m_values",action="write",status="replace") + +! Loops through spin flip iterations for T=1.0 and T-2.0 +do Tstep=1,2 + + k=1 + + T=Tstep + beta=1/(kB*T) + +! Creates new outfile for each T with m values +! for each one hundred flips + write(out_file_no,"(F3.1)") T + out_file=out_file_base//out_file_no + open (unit=out_unit,file=out_file,action="write",status="replace") + +!Initiate array to have every value equal to positive 1, spin up + do i=1,size + do j=1,size + spins(i,j)=1 + end do + end do + +! Resets variables + m_Avg=0 + everyHundred=1 + avgCount=1 + mCount=1 + avgM=0 + diffSqSum=0 + +! Iterate 5 million times, flipping a random spin each time + do flipSpins=1,5000000 + + rand = get_random_number() + +! Selects random numbers and uses to select i,j coors to +! determine which spin will be flipped + rand = get_random_number() + i_swap_value = get_random_whole_number(rand,size) + rand = get_random_number() + j_swap_value = get_random_whole_number(rand,size) + +! Flip Selected spin + spins(i_swap_value, j_swap_value)=spins(i_swap_value, j_swap_value)*(-1) + + e_new=0 + e_old=0 + +! Calculate energy contribution of flipped spin in new and old state while +! accounting for periodic boundary conditions + if (i_swap_value>1) then + e_new=e_new+(spins(i_swap_value,j_swap_value)*spins((i_swap_value-1),j_swap_value)) + e_old=e_old+(spins(i_swap_value,j_swap_value)*spins((i_swap_value-1),j_swap_value)*(-1)) + else + e_new=e_new+(spins(i_swap_value,j_swap_value)*spins((size),j_swap_value)) + e_old=e_old+(spins(i_swap_value,j_swap_value)*spins((size),j_swap_value)*(-1)) + end if + if (j_swap_value>1) then + e_new=e_new+(spins(i_swap_value,j_swap_value)*spins(i_swap_value,(j_swap_value-1))) + e_old=e_old+(spins(i_swap_value,j_swap_value)*spins(i_swap_value,(j_swap_value-1))*(-1)) + else + e_new=e_new+(spins(i_swap_value,j_swap_value)*spins(i_swap_value,(size))) + e_old=e_old+(spins(i_swap_value,j_swap_value)*spins(i_swap_value,(size))*(-1)) + end if + if (i_swap_value0) then + rand=get_random_number() + g=exp((-beta*delta_e)) + if(g<=rand) then + spins(i_swap_value,j_swap_value)=spins(i_swap_value,j_swap_value)*(-1) + end if + end if + +! Calculate and record m value every 100 steps + if (everyHundred == 100) then + m=0 + do i=1,size + do j=1,size + m=m+spins(i,j) + end do + end do + + m=m/(size*size) + if (m<0) then + m=(-1)*m + end if + m_T(2,k)=m + m_T(1,k)=flipSpins + m_T(1,k)=m_T(1,k)/1000000.0 + write (out_unit,*) m_T(1,k),m_T(2,k) + +! If the iteration is 3 million or above use to calculate average + if (flipSpins>2999999) then + m_avg=m_avg+m +! Calculate average m every 100,000 spin flip iterations + if (avgCount==1000) then + m_avg=m_avg/avgCount + avgCount=0 + m_bins(mCount)=m_avg + mCount=mCount+1 + m_avg=0 + end if + avgCount=avgCount+1 + end if + + end if + + k=k+1 + if (everyHundred==100) then + everyHundred=1 + else + everyHundred=everyHundred+1 + end if + + !Iterations end + end do + + mCount=mCount-1 + do i=1,mCount + avgM=avgM+m_bins(i) + end do + avgM=avgM/mCount + + do i=1,mCount + diff=avgM-m_bins(i) + diffSq=diff**2 + diffSqSum=diffSqSum+diffSq + end do + stdev=sqrt((diffSqSum/mCount)) + + write (out_unit2,*) T, avgM, stdev + close(out_unit) + +!Temp end +end do + +!Loops through spin flip iterations for T=2.1-2.9 +do Tstep=1,9 + + k=1 + + T=2+(Tstep/10.0) + beta=1/(kB*T) + write(out_file_no,"(F3.1)") T + out_file=out_file_base//out_file_no + open (unit=out_unit,file=out_file,action="write",status="replace") + do i=1,size + do j=1,size + spins(i,j)=1 + end do + end do + + m_Avg=0 + everyHundred=1 + avgCount=1 + mCount=1 + avgM=0 + diffSqSum=0 + do flipSpins=1,5000000 + + rand = get_random_number() + + rand = get_random_number() + i_swap_value = get_random_whole_number(rand,size) + + rand = get_random_number() + j_swap_value = get_random_whole_number(rand,size) + + spins(i_swap_value, j_swap_value)=spins(i_swap_value, j_swap_value)*(-1) + + e_new=0 + e_old=0 + + if (i_swap_value>1) then + e_new=e_new+(spins(i_swap_value,j_swap_value)*spins((i_swap_value-1),j_swap_value)) + e_old=e_old+(spins(i_swap_value,j_swap_value)*spins((i_swap_value-1),j_swap_value)*(-1)) + else + e_new=e_new+(spins(i_swap_value,j_swap_value)*spins((size),j_swap_value)) + e_old=e_old+(spins(i_swap_value,j_swap_value)*spins((size),j_swap_value)*(-1)) + end if + if (j_swap_value>1) then + e_new=e_new+(spins(i_swap_value,j_swap_value)*spins(i_swap_value,(j_swap_value-1))) + e_old=e_old+(spins(i_swap_value,j_swap_value)*spins(i_swap_value,(j_swap_value-1))*(-1)) + else + e_new=e_new+(spins(i_swap_value,j_swap_value)*spins(i_swap_value,(size))) + e_old=e_old+(spins(i_swap_value,j_swap_value)*spins(i_swap_value,(size))*(-1)) + end if + if (i_swap_value0) then + rand=get_random_number() + g=exp((-beta*delta_e)) + if(g<=rand) then + spins(i_swap_value,j_swap_value)=spins(i_swap_value,j_swap_value)*(-1) + end if + end if + + if (everyHundred == 100) then + m=0 + do i=1,size + do j=1,size + m=m+spins(i,j) + end do + end do + + m=m/(size*size) + if (m<0) then + m=(-1)*m + end if + m_T(2,k)=m + m_T(1,k)=flipSpins + m_T(1,k)= m_T(1,k)/1000000 + write (out_unit,*) m_T(1,k),m_T(2,k) + + if (flipSpins>2999999) then + m_avg=m_avg+m + if (avgCount==1000) then + m_avg=m_avg/avgCount + avgCount=0 + m_bins(mCount)=m_avg + mCount=mCount+1 + m_avg=0 + end if + avgCount=avgCount+1 + end if + + end if + + k=k+1 + if (everyHundred==100) then + everyHundred=1 + else + everyHundred=everyHundred+1 + end if + + !Iterations end + end do + mCount=mCount-1 + do i=1,mCount + avgM=avgM+m_bins(i) + end do + avgM=avgM/mCount + + do i=1,mCount + diff=avgM-m_bins(i) + diffSq=diff**2 + diffSqSum=diffSqSum+diffSq + end do + stdev=sqrt((diffSqSum/mCount)) + + write (out_unit2,*) T, avgM, stdev + close(out_unit) + +!Temp end +end do + +!Loops through spin flip iterations for T=3.0 and T-4.0 +do Tstep=3,4 + + k=1 + + T=Tstep + beta=1/(kB*T) + write(out_file_no,"(F3.1)") T + out_file=out_file_base//out_file_no + open (unit=out_unit,file=out_file,action="write",status="replace") + do i=1,size + do j=1,size + spins(i,j)=1 + end do + end do + + m_Avg=0 + everyHundred=1 + avgCount=1 + mCount=1 + avgM=0 + diffSqSum=0 + do flipSpins=1,5000000 + + rand = get_random_number() + + rand = get_random_number() + i_swap_value = get_random_whole_number(rand,size) + + rand = get_random_number() + j_swap_value = get_random_whole_number(rand,size) + + spins(i_swap_value, j_swap_value)=spins(i_swap_value, j_swap_value)*(-1) + + e_new=0 + e_old=0 + + if (i_swap_value>1) then + e_new=e_new+(spins(i_swap_value,j_swap_value)*spins((i_swap_value-1),j_swap_value)) + e_old=e_old+(spins(i_swap_value,j_swap_value)*spins((i_swap_value-1),j_swap_value)*(-1)) + else + e_new=e_new+(spins(i_swap_value,j_swap_value)*spins((size),j_swap_value)) + e_old=e_old+(spins(i_swap_value,j_swap_value)*spins((size),j_swap_value)*(-1)) + end if + if (j_swap_value>1) then + e_new=e_new+(spins(i_swap_value,j_swap_value)*spins(i_swap_value,(j_swap_value-1))) + e_old=e_old+(spins(i_swap_value,j_swap_value)*spins(i_swap_value,(j_swap_value-1))*(-1)) + else + e_new=e_new+(spins(i_swap_value,j_swap_value)*spins(i_swap_value,(size))) + e_old=e_old+(spins(i_swap_value,j_swap_value)*spins(i_swap_value,(size))*(-1)) + end if + if (i_swap_value0) then + rand=get_random_number() + g=exp((-beta*delta_e)) + if(g<=rand) then + spins(i_swap_value,j_swap_value)=spins(i_swap_value,j_swap_value)*(-1) + end if + end if + + if (everyHundred == 100) then + m=0 + do i=1,size + do j=1,size + m=m+spins(i,j) + end do + end do + + m=m/(size*size) + if (m<0) then + m=(-1)*m + end if + m_T(2,k)=m + m_T(1,k)=flipSpins + m_T(1,k)= m_T(1,k)/1000000 + write (out_unit,*) m_T(1,k),m_T(2,k) + + if (flipSpins>2999999) then + m_avg=m_avg+m + if (avgCount==1000) then + m_avg=m_avg/avgCount + avgCount=0 + m_bins(mCount)=m_avg + mCount=mCount+1 + m_avg=0 + end if + avgCount=avgCount+1 + end if + + end if + + k=k+1 + if (everyHundred==100) then + everyHundred=1 + else + everyHundred=everyHundred+1 + end if + + !Iterations end + end do + mCount=mCount-1 + do i=1,mCount + avgM=avgM+m_bins(i) + end do + avgM=avgM/mCount + + do i=1,mCount + diff=avgM-m_bins(i) + diffSq=diff**2 + diffSqSum=diffSqSum+diffSq + end do + stdev=sqrt((diffSqSum/mCount)) + + write (out_unit2,*) T, avgM, stdev + close(out_unit) + +!Temp end +end do + +close(out_unit2) + + +contains + + function get_random_number() + integer :: i_seed + integer, dimension(:), ALLOCATABLE :: Set_seed + integer, dimension(1:8) :: dateSet_seed + real :: get_random_number, r + + CALL RANDOM_SEED(size=i_seed) + ALLOCATE(Set_seed(1:i_seed)) + CALL RANDOM_SEED(get=Set_seed) + CALL DATE_AND_TIME(values=dateSet_seed) + Set_seed(i_seed)=dateSet_seed(8) + !Set_seed(2)=dateSet_seed(8) + Set_seed(1)=dateSet_seed(8)*dateSet_seed(6) + CALL RANDOM_SEED(put=Set_seed) + DEALLOCATE(Set_seed) + + CALL RANDOM_NUMBER(r) + get_random_number = r + end function get_random_number + + function get_random_whole_number(y,x) + integer :: get_random_whole_number, x + real :: y + get_random_whole_number = (x*y)+1 + if (get_random_whole_number==(x+1)) then + get_random_whole_number = (x*y)+1 + end if + end function + + +end program isingtest \ No newline at end of file diff --git a/isingRunandPlotBashScript.txt b/isingRunandPlotBashScript.txt new file mode 100755 index 0000000..7fcf5f3 --- /dev/null +++ b/isingRunandPlotBashScript.txt @@ -0,0 +1,9 @@ +#!/bin/bash -login + +gfortran isingCode.f95 -o test + +./test + +gnuplot> load create_mvsSteps_Plot.p + +gnuplot> load create_mvsT_Plot.p