ising.f90

ising.f90

program izing
    implicit none


    integer, parameter :: L = 16
    integer, parameter :: mc_steps = 500
    integer, parameter :: sim_steps = 100
    double precision, parameter :: B = 0.0
    double precision, parameter :: J_ising = 1.0 
    double precision, parameter :: T_init = 0.5
    double precision, parameter :: T_final = 4.0
    double precision, parameter :: T_step = 0.002
    integer, parameter :: T_iters = floor((T_final-T_init)/T_step)



    integer :: i, j, k
    integer :: lattice(L, L)

    ! double precision :: find_mean, find_var, find_magnetization, find_energy, T
    double precision :: temp_E_arr(sim_steps), temp_M_arr(sim_steps), T

    ! initializing the lattice
    do i = 1, L
        do j = 1, L
            if (mod((i+j),2) == 0) then
                lattice(i, j) = -1
            else
                lattice(i, j) = 1
            end if
        end do
    end do


    ! driver code

    open(file="E_vs_T.dat", unit=1)
    open(file="M_vs_T.dat", unit=2)
    open(file="Cv_vs_T.dat", unit=3)
    open(file="Chi_vs_T.dat", unit=4)

    do i = 1, T_iters

        print*, "we are at ", i, "th", " iteration of total ", T_iters, " iterations"

        T = T_init + i * T_step
        ! print*, "T1=",T

        do j = 1, mc_steps
            ! T = T_init + i * T_step
            ! print*, "T2=",T
            ! print*, i
            call one_mcs(lattice, L, T, J_ising, B)
        end do

        do k = 1, sim_steps

            ! T = T_init + i * T_step
            ! print*, i
            ! print*, "T3=",T
            call one_mcs(lattice, L, T, J_ising, B)

            temp_E_arr(k) = find_energy(B, J_ising, lattice, L)
            temp_M_arr(k) = find_magnetization(lattice, L)

            ! write(1, '(f8.4)') T, find_mean(temp_E_arr, sim_steps)
            ! print*, T
            ! print*, find_mean(temp_E_arr, sim_steps)
        end do
        write(1, *) T, find_mean(temp_E_arr, sim_steps)/(L**2)

    end do
    close(1)
    close(2)
    close(3)
    close(4)

    contains

    subroutine one_mcs(Lattice, size_lattice, temparature, J_ising, B)
        implicit none

        integer :: i, size_lattice, Lattice(size_lattice, size_lattice), p, q
        double precision :: E, M, del_E, rand, prob, pr, qr
        double precision, intent(in) :: temparature, J_ising, B

        E = find_energy(B, J_ising, Lattice, size_lattice)
        ! print*, "temparature1 is", temparature

        do i = 1, size_lattice**2

            if (mod(i, size_lattice**2) == 0) then
                ! print*, temparature
            end if

            call random_number(pr)
            call random_number(qr)
            p = floor(pr*(size_lattice+1))
            q = floor(qr*(size_lattice+1))

            del_E = 2 * J_ising * (Lattice(p, q) * (Lattice(p, mod(q + 1, size_lattice)) +&
            Lattice(p, mod(q - 1, size_lattice)) + Lattice(mod(p + 1, size_lattice), q) +&
            Lattice(mod(p - 1, size_lattice), q))) + 2*Lattice(p, q)*B


            if (del_E < 0) then
                lattice(p, q) = -lattice(p, q)
                E = E + del_E
                M = M - 2 * lattice(p, q)

            else
                call random_number(rand)
                prob = exp(-del_E/temparature)

                if (rand < prob) then
                    lattice(p, q) = -lattice(p, q)
                    E = E + del_E
                    M = M - 2 * lattice(p, q)
                end if

            end if

        end do
        ! print*, "temparature 2 is", temparature

    end subroutine one_mcs


    double precision function find_energy(B, J_ising, Lattice, size_lattice)
        implicit none

        integer :: i, j, size_lattice

        double precision :: E
        integer, intent(in) :: Lattice(size_lattice, size_lattice)
        double precision, intent(in) :: J_ising, B

        E = 0.0
        do i = 1, size_lattice
            do j = 1, size_lattice
                E = E - J_ising * (Lattice(i, j) * (Lattice(i, mod(j + 1, size_lattice)) +&
                Lattice(i, mod(j - 1, size_lattice)) +&
                Lattice(mod(i + 1, size_lattice), j) +&
                Lattice(mod(i - 1, size_lattice), j))) + Lattice(i, j)*B
            end do
        end do

        find_energy = E/2

    end function find_energy



    double precision function find_magnetization(Lattice, size_lattice)
        implicit none

        integer :: i, j

        double precision :: M = 0.0
        integer, intent(in) :: Lattice(size_lattice, size_lattice), size_lattice

        do i = 1, size_lattice
            do j = 1, size_lattice
                M = M + Lattice(i, j)
            end do
        end do

        find_magnetization = M

    end function find_magnetization



    double precision function find_mean(array, size_arr)
        implicit none

        integer :: i, size_arr
        double precision :: mean
        double precision :: some
        double precision, intent(in) :: array(size_arr)

        some = 0
        do i = 1, size_arr
            some = some + array(i)
        end do

        mean = some/size_arr

        find_mean = mean

    end function find_mean



    double precision function find_var(array, size_arr)
        implicit none

        integer :: i, size_arr
        double precision, intent(in) :: array(size_arr)
        double precision :: var
        ! double precision :: find_mean
        double precision :: sq_mean, mean_sq, sq_arr(size_arr)

        mean_sq = (find_mean(array, size_arr))**2

        do i = 1, size_arr
            sq_arr(i) = (array(i))**2
        end do

        sq_mean = find_mean(sq_arr, size_arr)

        var = sq_mean - mean_sq
        
        find_var = var
    end function find_var

end program izing