!caf Ising01.f90 -o is; echo 12 0.20 1000 1 |  cafrun -n 4 ./is
!
!
!If you need more processes than the cores of your CPU, use:
!cafrun -n 20 -oversubscribe ./is
!
!================== Ising01.f90 =======
!--------------------------------------------------
module          Ising01
 use, intrinsic          :: iso_fortran_env
 implicit none
 integer, parameter      :: dp = real64
 integer                 :: L,Lloc,Limg
 integer                 :: N,Nloc,Nimg
 integer ,allocatable    :: s(:,:)[:,:]
 real(dp),dimension(0:4) :: prob
 real(dp)                :: beta
 integer                 :: nsweep,start
 real(dp)                :: r
 integer                 :: myimg, myi, myj
 integer                 :: i,j,iloc,jloc,ip,jp
 integer                 :: fout
 character(1000)         :: fname_out
 real(dp)                :: acceptance
!--------------------------------------------------
contains
!--------------------------------------------------
 subroutine     init !inherits implicit none!
  integer :: i
!----------------------
  Nimg  = num_images()
  myimg = this_image()
  if(myimg == 1) then
   print '(A)', '# Enter L, beta, nsweep, start (0 cold/1 hot):'
   read    *  ,          L, beta, nsweep, start
  end if
  call co_broadcast(L     ,1);    ! sync all not necessary
  call co_broadcast(beta  ,1); 
  call co_broadcast(nsweep,1); 
  call co_broadcast(start ,1); 
  N = L*L; Limg = int(sqrt(Nimg+1.0e-6_dp)); Lloc = L/Limg; Nloc=Lloc*Lloc
  if(      Nimg  /= Limg**2) call co_abort('Nimg /= Limg**2')
  if(mod(L,Limg) /= 0      ) call co_abort('mod(L,Limg) /= 0')
  ALLOCATE(s(Lloc,Lloc)[Limg,*])
  
  myi = this_image(s,1) ; myj = this_image(s,2)
!----------------------
  if(myimg == 1) then
   write(fname_out,'(A,I0.3,A,F5.3,A,I0.2,A)') 'L',L,'b',beta,'i',myimg,'.out'
   open(newunit=fout,file=fname_out)
   write(fout, *)'#       L  beta  nsweep  start  myimg  nimg     (start: 0 cold/1 hot)'
   write(fout, *)'#     ',L, beta, nsweep, start, myimg, nimg
  end if
  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  iloc=1,Lloc
    do jloc=1,Lloc
     call random_number(r)
     if(r .lt. 0.5_dp)then
      s(iloc,jloc) =  1
     else
      s(iloc,jloc) = -1
     endif
    end do
   end do
  case default
   print *,'init: start= ',start,' not valid. Exiting...',myimg
   error stop
  end select
  
  acceptance = 0.0_dp
 end subroutine init
!--------------------------------------------------
!Naive implementation, generalizing the previous code: one needs a sync all at each Metropolis check, which is very expensive
subroutine      met()
 integer :: k
 integer :: nniloc,nnjloc
 integer :: nnip  ,nnjp
 integer :: snn,dE

 do k=1,Nloc
!--------------------------------------------------
!pick a random site:
  call random_number(r) ; iloc = INT(Lloc*r)+1
  call random_number(r) ; jloc = INT(Lloc*r)+1
!snn=sum of neighboring spins:
  snn = 0
!--------------------------------------------------
! East:
  nniloc    = iloc + 1 ; nnjloc = jloc
  if(nniloc > Lloc ) then
   nniloc   = 1
   nnip     = myi  + 1 ; nnjp   = myj
   if(nnip  > Limg ) nnip = 1
   snn = snn + s(nniloc,nnjloc)[nnip,nnjp]
  else
   snn = snn + s(nniloc,nnjloc)
  end if
!--------------------------------------------------
! North:
  nniloc    = iloc     ; nnjloc = jloc + 1
  if(nnjloc > Lloc ) then
   nnjloc   = 1
   nnip     = myi      ; nnjp   = myj  + 1
   if(nnjp  > Limg ) nnjp = 1
   snn = snn + s(nniloc,nnjloc)[nnip,nnjp]
  else
   snn = snn + s(nniloc,nnjloc)
  end if
!--------------------------------------------------
! West:
  nniloc    = iloc - 1 ; nnjloc = jloc
  if(nniloc < 1 ) then
   nniloc   = Lloc
   nnip     = myi  - 1 ; nnjp   = myj
   if(nnip  < 1 ) nnip = Limg
   snn = snn + s(nniloc,nnjloc)[nnip,nnjp]
  else
   snn = snn + s(nniloc,nnjloc)
  end if
!--------------------------------------------------
! South:
  nniloc    = iloc     ; nnjloc = jloc - 1
  if(nnjloc < 1 ) then
   nnjloc   = Lloc
   nnip     = myi      ; nnjp   = myj  - 1
   if(nnjp  < 1 ) nnjp = Limg
   snn = snn + s(nniloc,nnjloc)[nnip,nnjp]
  else
   snn = snn + s(nniloc,nnjloc)
  end if
!--------------------------------------------------
!dE=change in energy/2:
  dE=snn*s(iloc,jloc)
!Oooops, necessary but really expensive! Without it, images may even be at different k-values and one 
!effectively simulates a "time-skewed" version of the Ising model.
  sync all                       
!Metropolis test:
  if(dE.le.0)then
   s(iloc,jloc)  = -s(iloc,jloc) !accept
   acceptance    = acceptance + 1.0_dp
  else 
   call random_number(r)
   if(r < prob(dE)) then
    s(iloc,jloc) = -s(iloc,jloc) !accept
    acceptance   = acceptance + 1.0_dp
   end if
  end if
  sync all !same: all images must have the spin updates synced
 enddo !do k=1,N: end sweep
end subroutine  met
!--------------------------------------------------
subroutine      measure()
 integer ene, mag
 ene = E(); mag = M() ! all images must execute the functions, each computes its own portion of ene+mag
 if(myimg == 1) write(fout,*) ene, mag
end subroutine  measure
!--------------------------------------------------
function        E()
 integer :: E
 integer :: en    ,snn
 integer :: nni   ,nnj
 integer :: nniloc,nnjloc
 integer :: nnip  ,nnjp
 en = 0
 do iloc=1,Lloc
  do jloc=1,Lloc
!Sum of neighboring spins: only forward nn necessary in the sum
   snn = 0
!--------------------------------------------------
! North:
   nniloc    = iloc + 1 ; nnjloc = jloc
   if(nniloc > Lloc ) then
    nniloc   = 1
    nnip     = myi  + 1 ; nnjp   = myj
    if(nnip  > Limg ) nnip = 1
    snn = snn + s(nniloc,nnjloc)[nnip,nnjp]
   else
    snn = snn + s(nniloc,nnjloc)
   end if
!--------------------------------------------------
! East:
   nniloc    = iloc     ; nnjloc = jloc + 1
   if(nnjloc > Lloc ) then
    nnjloc   = 1
    nnip     = myi      ; nnjp   = myj  + 1
    if(nnjp  > Limg ) nnjp = 1
    snn = snn + s(nniloc,nnjloc)[nnip,nnjp]
   else
    snn = snn + s(nniloc,nnjloc)
   end if
!--------------------------------------------------
   en=en+snn*s(iloc,jloc)
  end do
 end do
 call co_sum(en)
 e = -en
end function    E
!--------------------------------------------------
function        M()
 integer :: M,mag
 mag=SUM(s)
 call co_sum(mag)
 M  = mag
end function    M
!--------------------------------------------------
 subroutine     endsim
  call co_sum(acceptance)
  if(myimg == 1) write(fout,'(A,F15.6)') '# acceptance= ',acceptance/N/nsweep
 end subroutine endsim
!--------------------------------------------------
 subroutine         co_abort(errmes) 

  character(*), intent(in) :: errmes

  write(error_unit,'(A,I3)') 'Error: ' //trim(errmes) //' Image: ',myimg

  error stop 1

 end subroutine     co_abort
!--------------------------------------------------
end module      Ising01
!==================-------------===================
program         ising2dProgram
 use Ising01
 implicit none
 integer isweep
 call  init
 do isweep = 1, nsweep
  call met
  call measure
 end do
 call  endsim
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/>.
!  -----------------------------------------------------------------------