! How to remove ^M in source code on linux %s/^M//g (^M = Ctrl+v+m)
! Note on subroutine random_number
! Returns a single pseudorandom number or an array of pseudorandom numbers from the uniform distribution over the range [0, 1).
! The runtime-library implements the xorshift1024* random number generator(RNG). This generator has a period of 2^{1024} - 1,
! and when using multiple threads up to 2^{512} threads can each generate 2^{512} random numbers before any aliasing occurs.

module parameters
    integer, parameter :: latt_leng = 80
    integer, parameter :: EXCH_CONST = 1 ! exchange constant
    integer, dimension(-1 : latt_leng, -1 : latt_leng) :: spin_of ! (-1 ... latt_leng) X (-1 ... latt_leng)
    double precision, dimension(-2 : 2) :: accept_prob_chart
end module


program Ising_Model
    use parameters
    implicit none ! Require all variables to be explicitly declared
    integer(kind=8), parameter :: MC_max = 1E7 * latt_leng**2  
    double precision, parameter :: T = 2.269d0 ! Unit = kT/J, Tc computed with each of the physical observables
    
    call create_lattice_random 
    if(EXCH_CONST == -1) then ! Anti-ferromagnetic, EXCH_CONST = -1
        accept_prob_chart = (/ exp(-8d0 / T), exp(-4d0 / T), &
        1d0, exp(4d0 / T), exp(8d0 / T) /) ! Acceptance probability chart
    else ! Ferromagnetic, EXCH_CONST = 1
        accept_prob_chart = (/ exp(8d0 / T), exp(4d0 / T), &
        1d0, exp(-4d0 / T), exp(-8d0 / T) /)
    end if
    call flip_sites_metropolis(MC_max) ! Iterate MC simulation MC_max times per temperature
end program


subroutine create_lattice_random
    use parameters
    integer :: x, y
    double precision :: u

    do x = 0, latt_leng - 1
        do y = 0, latt_leng - 1
            call random_number(u) ! Return random real number in [0,1)
            if (u >= 0.5) then
                spin_of(x, y) = 1 ! Assign random spin values
            else
                spin_of(x, y) = -1
            end if
        end do
    end do

    do x = 0, latt_leng - 1 ! Assign spins to top and bottom line of padding sites to make it periodic
        spin_of(x, -1) = spin_of(x, latt_leng - 1)
        spin_of(x, latt_leng) = spin_of(x, 0)
    end do
    do y = 0, latt_leng -1  ! Assign spins to left and right line of padding sites to make it periodic
        spin_of(-1, y) = spin_of(latt_leng - 1, y)
        spin_of(latt_leng, y) = spin_of(0, y)
    end do
end subroutine


subroutine create_lattice_aligned
    use parameters
    integer :: x, y

    do x = 0, latt_leng - 1 ! Assign spins on each site
        do y = 0, latt_leng - 1
            spin_of(x, y) = 1 ! All up spins
        end do
    end do
    do x = 0, latt_leng - 1 ! Assign spins to top and bottom line of padding sites to make it periodic
        spin_of(x, -1) = spin_of(x, latt_leng - 1)
        spin_of(x, latt_leng) = spin_of(x, 0)
    end do
    do y = 0, latt_leng -1  ! Assign spins to left and right line of padding sites to make it periodic
        spin_of(-1, y) = spin_of(latt_leng - 1, y)
        spin_of(latt_leng, y) = spin_of(0, y)
    end do
end subroutine


subroutine flip_sites_metropolis(MC_max)
    use parameters
    integer :: first_open = 1
    integer :: x, y, sum_spin_nbr, sum_nbr_spin
    integer :: out_unit1 = 10, out_unit2 = 20 !,out_unit3 = 30  
    integer(kind=8), intent(in) :: MC_max
    integer(kind=8) :: MC_step
    double precision :: accept_prob, tot_en, tot_mag, tot_en2, tot_E, tot_M, tot_E2
    double precision :: u, v, r

    n_sample = 0 ! Initialize variable to be used in averaging samples
    sum_tot_E = 0d0
    sum_tot_M = 0d0
    sum_tot_E2 = 0d0
    call compute_tot_E_M_E2(tot_en, tot_mag, tot_en2) ! Compute initial values of physical observables
    tot_E = tot_en
    tot_M = tot_mag
    tot_E2 = tot_en2
    do MC_step = 0, MC_max
        call random_number(u) ! Return random real number in [0,1)
        call random_number(v)
        x = FLOOR(u * latt_leng) ! Pick random site, [0, latt_leng - 1]
        y = FLOOR(v * latt_leng)
        sum_nbr_spin = spin_of(x - 1, y) + spin_of(x + 1, y) + spin_of(x, y + 1) + spin_of(x, y - 1)
        accept_prob = accept_prob_chart((spin_of(x, y) * sum_nbr_spin)/2) ! Sigma(SS) = spin_of(x, y) * sum_surrounding = -4, -2 , 0, 2, 4
        call random_number(r) ![0, 1)
        if (r < accept_prob) then
            spin_of(x, y) = - spin_of(x, y) ! If random number < accept_prob => flip site
            
            if (x == 0) then ! Update padding sites
                spin_of(latt_leng, y) = spin_of(x, y) 
            else if (x == latt_leng - 1) then
                spin_of(-1, y) = spin_of(x, y)
            end if 
            if (y == 0) then
                spin_of(x, latt_leng) = spin_of(x, y)
            else if (y == latt_leng - 1) then
                spin_of(x, -1) = spin_of(x, y)
            end if
            
            tot_M = tot_M + 2*spin_of(x, y) ! Update physical observables
            sum_spin_nbr = spin_of(x + 1, y) + spin_of(x - 1, y) + spin_of(x, y + 1) + spin_of(x, y - 1)
            tot_E = tot_E + (- EXCH_CONST) * spin_of(x, y) * sum_spin_nbr
            tot_E2 = tot_E**2
        end if
        
        if (MOD(MC_step, latt_leng**2) == 0) then
            if(first_open == 1) then
                open (out_unit1, FILE = 'output_ising_model_fixed_T_E_change_plot.txt', action = "write", status = 'replace')
                write (out_unit1,*) MC_step / latt_leng**2 , " ", tot_E
                close(out_unit1)
                open (out_unit2, FILE = 'output_ising_model_fixed_T_M_change_plot.txt', action = "write", status = 'replace')
                write (out_unit2,*) MC_step / latt_leng**2, " ", tot_M
                close(out_unit2)
                first_open = -1
            else 
                open (out_unit1, FILE = 'output_ising_model_fixed_T_E_change_plot.txt', action = "write", position = 'append')
                write (out_unit1,*) MC_step / latt_leng**2, " ", tot_E
                close(out_unit1)
                open (out_unit2, FILE = 'output_ising_model_fixed_T_M_change_plot.txt', action = "write", position = 'append')
                write (out_unit2,*) MC_step / latt_leng**2, " ", tot_M
                close(out_unit2)
            end if
        end if
    end do
end subroutine


subroutine compute_tot_E_M_E2(tot_E, tot_M, tot_E2)
    use parameters
    integer :: x, y, sum_spin_nbr
    double precision, intent(out) :: tot_E, tot_M, tot_E2
    
    tot_E = 0d0
    tot_M = 0d0
    tot_E2 = 0d0
    do x = 0, latt_leng - 1
        do y = 0, latt_leng - 1
            tot_M = tot_M + spin_of(x, y)
            sum_spin_nbr = spin_of(x + 1, y) + spin_of(x - 1, y) + spin_of(x, y + 1) + spin_of(x, y - 1)
            tot_E = tot_E + (- EXCH_CONST) * spin_of(x, y) * sum_spin_nbr
        end do
    end do
    tot_E = tot_E / 2
    tot_E2 = tot_E**2
end subroutine
