!================== IsingVanilla.f90 ==============
!--------------------------------------------------
module          IsingVanilla
 use, intrinsic          :: iso_fortran_env
 implicit none
 SAVE
 integer, parameter      :: dp = real64
 integer                 :: L
 integer                 :: N
 integer                 :: XNN, YNN
 integer ,allocatable    :: s(:)
 real(dp),dimension(0:4) :: prob
 real(dp)                :: beta
 integer                 :: nsweep,start
 real(dp)                :: r
!--------------------------------------------------
contains
!--------------------------------------------------
 subroutine     init !inherits implicit none!
  integer :: i
!----------------------
!Ask user for input:
  print *,'# Enter L, beta, nsweep, start (0 cold/1 hot):'
  read  *,         L, beta, nsweep, start
  print *,'#     ',L, beta, nsweep, start
  N = L*L; XNN = 1; YNN = L
  ALLOCATE(s(N))

  call random_init(.false.,.true.) ! random_init(REPEATABLE, IMAGE_DISTINCT)

!initialize probabilities for E_\nu > E_mu
  prob=0.0_dp
! Sum_<ik> s_i^\mu = -4, -2, 0, +2, +4
! Need:     exp(-2*beta*4), exp(- 2*beta*2) 
! and not:  exp( 2*beta*0), exp(  2*beta*2),exp(2*beta*4)
  do i=2,4,2 !i = dE/2 = (E_nu-E_mu)/2=2,4
   prob(i) = exp(-2.0_dp*beta*i)
  enddo
  !initial configuration:
  select case(start)
  case(0)!cold:
   s = 1 !all s(i) = 1
  case(1)!hot:
   do i=1,N
    call random_number(r)
    if(r .lt. 0.5_dp)then
     s(i) =  1
    else
     s(i) = -1
    endif
   enddo
  case default
   print *,'init: start= ',start,' not valid. Exiting...'
   stop
  end select
  
 end subroutine init
!--------------------------------------------------
subroutine      met()
 integer :: i,k
 integer :: nn,snn,dE

 do k=1,N
!pick a random site:
  call random_number(r)
  i = INT(N*r)+1
!snn=sum of neighboring spins:
  snn = 0
  nn=i+XNN;if(nn.gt.N)nn=nn-N;snn = snn + s(nn)
  nn=i-XNN;if(nn.lt.1)nn=nn+N;snn = snn + s(nn)
  nn=i+YNN;if(nn.gt.N)nn=nn-N;snn = snn + s(nn)
  nn=i-YNN;if(nn.lt.1)nn=nn+N;snn = snn + s(nn)
!dE=change in energy/2:
  dE=snn*s(i)
!flip:
  if(dE.le.0)then
   s(i) = -s(i) !accept
  else 
   call random_number(r)
   if(r < prob(dE)) s(i) = -s(i) !accept
  end if
 enddo !do k=1,N: end sweep
end subroutine  met
!--------------------------------------------------
subroutine      measure()
 print *, E(),M()
end subroutine  measure
!--------------------------------------------------
function        E()
 integer :: E
 integer :: en,sum,i,nn
 en = 0
 do i=1,N
!Sum of neighboring spins: only forward nn necessary in the sum
  sum = 0
  nn=i+XNN;if(nn.gt.N)nn=nn-N;sum = sum + s(nn)
  nn=i+YNN;if(nn.gt.N)nn=nn-N;sum = sum + s(nn)
  en=en+sum*s(i)
 enddo
 e = -en
end function    E
!--------------------------------------------------
function        M()
 integer :: M
 M=SUM(s)
end function    M
!--------------------------------------------------
end module      IsingVanilla
!==================-------------===================
program         ising2dProgram
 use IsingVanilla
 implicit none
 integer isweep
 call init
 do isweep = 1, nsweep
  call met
  call measure
 end do
end program ising2dProgram
!=================================================
!  ---------------------------------------------------------------------
!  Copyright by Konstantinos N. Anagnostopoulos (2004-2022)
!  Physics Dept., National Technical University,
!  konstant@mail.ntua.gr, www.physics.ntua.gr/~konstant
!  
!  This program is free software: you can redistribute it and/or modify
!  it under the terms of the GNU General Public License as published by
!  the Free Software Foundation, version 3 of the License.
!  
!  This program is distributed in the hope that it will be useful, but
!  WITHOUT ANY WARRANTY; without even the implied warranty of
!  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
!  General Public License for more details.
!  
!  You should have received a copy of the GNU General Public Liense along
!  with this program.  If not, see <http://www.gnu.org/licenses/>.
!  -----------------------------------------------------------------------
