!   Program to solve the two-dimensional Ising model 
!   with zero external field.
!   The coupling constant J = 1
!   Boltzmann's constant = 1, temperature has thus dimension energy
!   Metropolis sampling is used. Periodic boundary conditions.


!    This module contains all constants and declarations 
!    of variables read in by the function read_data. These
!    variables are used by many functions.

MODULE constants
  INTEGER,  PARAMETER :: dp = KIND(1.0D0)
  INTEGER, PARAMETER :: dpc = KIND((1.0D0,1.0D0))
END MODULE constants


!    F90 program library, adapted from Numerical Recipes
!    All functions have been translated to F90 from F77

MODULE F90library
  USE constants
CONTAINS

  !
  !     The function
  !           ran0()
  !     is an "Minimal" random number generator of Park and Miller
  !     (see Numerical recipe page 279). Set or reset the input value
  !     idum to any integer value (except the unlikely value MASK)
  !     to initialize the sequence; idum must not be altered between
  !     calls for sucessive deviates in a sequence.
  !     The function returns a uniform deviate between 0.0 and 1.0.
  !


  REAL(DP) FUNCTION ran0(idum)
    IMPLICIT NONE
    INTEGER :: idum,ia,im,iq,ir,mask,k
    REAL(DP) :: am
    PARAMETER (ia=16807,im=2147483647,am=1./im,iq=127773,ir=2836,mask=123459876)

    idum=ieor(idum,MASK)
    k=idum/IQ
    idum=IA*(idum-k*IQ)-IR*k
    IF  (idum < 0) idum=idum+IM
    ran0=am*idum
    idum=ieor(idum,MASK)

  END FUNCTION ran0


  !
  !      The function
  !           ran1()
  !      is an "Minimal" random number generator of Park and Miller
  !      (see Numerical recipe page 280) with Bays-Durham shuffle and
  !      added safeguards. Call with idum a negative integer to initialize;
  !      thereafter, do not alter idum between sucessive deviates in a
  !      sequence. RNMX should approximate the largest floating point value
  !      that is less than 1.
  !      The function returns a uniform deviate between 0.0 and 1.0
  !      (exclusive of end-point values).
  !



  REAL(DP) FUNCTION ran1(idum)
    IMPLICIT NONE
    INTEGER :: idum,ia,im,iq,ir,ntab,ndiv
    REAL(DP) :: am,eps,rnmx
    PARAMETER (ia=16807,im=2147483647,am=1./im,iq=127773,ir=2836, &
         ntab=32,ndiv=1+(im-1)/ntab,eps=1.2e-7,rnmx=1.-eps)
    INTEGER :: j,k,iv(ntab),iy
    SAVE iv,iy
    DATA iv /NTAB*0/, iy /0/

    IF ((idum == 0).OR.(iy==0)) THEN
       idum=MAX(-idum,1)
       DO j=ntab+8,1,-1
          k=idum/iq
          idum=ia*(idum-k*iq)-ir*k
          IF (idum < 0) idum=idum+im
          IF (j <= ntab) iv(j)=idum
       ENDDO
       iy=iv(1)
    ENDIF
    k=idum/IQ
    idum=ia*(idum-k*iq)-ir*k
    if (idum < 0) idum=idum+IM
    j=1+iy/ndiv
    iy=iv(j)
    iv(j)=idum
    ran1=MIN(am*iy,rnmx)

  END FUNCTION ran1


  REAL(DP) FUNCTION ran(idum)
    IMPLICIT NONE
    INTEGER, PARAMETER :: K4B=selected_int_kind(9)
    INTEGER(K4B), INTENT(INOUT) :: idum 
    INTEGER :: ia,im,iq,ir
    REAL(DP), SAVE :: am
    PARAMETER (ia=16807,im=2147483647,iq=127773,ir=2836)
    INTEGER(K4B), SAVE :: ix=-1, iy=-1, k 

    IF ((idum <= 0).OR.(iy < 0)) THEN
       am=NEAREST(1.,-1.)/im
       iy=IOR(IEOR(888889999,ABS(idum)),1)
       ix=IEOR(777755555,ABS(idum))
       idum=ABS(idum)+1
    ENDIF
    ix=IEOR(ix,ISHFT(ix,13))
    ix=IEOR(ix,ISHFT(ix,-17))
    ix=IEOR(ix,ISHFT(ix,5))         
    k=iy/iq
    iy=ia*(iy-k*iq)-ir*k
    IF ( iy < 0 ) iy=iy+im
    ran=am*IOR(IAND(im,IEOR(ix,iy)),1)

  END FUNCTION ran

  !
  !     The function 
  !         ran2
  !     is a long periode (> 2 x 10^18) random number generator of 
  !     L'Ecuyer and Bays-Durham shuffle and added safeguards.
  !     Call with idum a negative integer to initialize; thereafter,
  !     do not alter idum between sucessive deviates in a
  !     sequence. RNMX should approximate the largest floating point value
  !     that is less than 1.
  !     The function returns a uniform deviate between 0.0 and 1.0
  !     (exclusive of end-point values).
  !

  REAL(DP)  FUNCTION ran2(idum)
    IMPLICIT NONE
    INTEGER :: idum,im1,im2,imm1,ia1,ia2,iq1,iq2,ir1,ir2,ntab,ndiv
    REAL(DP) :: am,eps,rnmx
    PARAMETER (im1=2147483563,im2=2147483399,am=1./im1,imm1=im1-1, &
         ia1=40014,ia2=40692,iq1=53668,iq2=52774,ir1=12211,ir2=3791, &
         ntab=32,ndiv=1+imm1/ntab,eps=1.2e-7,rnmx=1.-eps)
    INTEGER :: idum2,j,k,iv(NTAB),iy
    SAVE iv,iy,idum2
    DATA idum2/123456789/, iv/NTAB*0/, iy/0/
    IF  (idum <= 0) then
       idum=MAX(-idum,1)
       idum2=idum
       DO j=NTAB+8,1,-1
          k=idum/IQ1
          idum=ia1*(idum-k*iq1)-k*ir1
          IF  (idum < 0) idum=idum+IM1
          IF (J <= NTAB) IV(J)=IDUM
       ENDDO
       IY=IV(1)
    ENDIF
    k=idum/iq1
    idum=ia1*(idum-k*iq1)-k*ir1
    if (idum < 0) idum=idum+im1
    k=idum2/iq2
    idum2=ia2*(idum2-k*iq2)-k*ir2
    IF (idum2 < 0) idum2=idum2+im2
    j=1+iy/ndiv
    iy=iv(j)-idum2
    iv(j)=idum
    IF (iy < 1)iy=iy+imm1
    ran2=MIN(am*iy,rnmx)

  END FUNCTION ran2


  !
  !     The function
  !        ran3
  !     returns a uniform random number deviate between 0.0 and 1.0. Set
  !     the idum to any negative value to initialize or reinitialize the
  !     sequence. Any large MBIG, and any small (but still large) MSEED
  !     can be substituted for the present values. 
  !

  REAL(DP) FUNCTION ran3(idum)
    IMPLICIT NONE
    INTEGER :: idum
    INTEGER :: mbig,mseed,mz
    REAL(DP) ::  fac
    PARAMETER (mbig=1000000000,mseed=161803398,mz=0,fac=1./mbig)
    INTEGER :: i,iff,ii,inext,inextp,k
    INTEGER :: mj,mk,ma(55)
    SAVE iff,inext,inextp,ma
    DATA iff /0/

    IF ( (idum < 0) .or. (iff == 0) ) THEN
       iff=1
       mj=mseed-IABS(idum)
       mj=MOD(mj,mbig)
       ma(55)=mj
       mk=1
       DO i=1,54
          ii=MOD(21*i,55)
          ma(ii)=mk
          mk=mj-mk
          IF(mk < mz)mk=mk+mbig
          mj=ma(ii)
       ENDDO
       DO k=1,4
          DO i=1,55
             ma(i)=ma(i)-ma(1+MOD(i+30,55))
             IF (ma(i) < mz)ma(i)=ma(i)+mbig
          ENDDO
       ENDDO
       inext=0
       inextp=31
       idum=1
    ENDIF
    inext=inext+1
    IF (inext == 56) inext=1
    inextp=inextp+1
    IF (inextp == 56) inextp=1
    mj=ma(inext)-ma(inextp)
    IF (mj < mz) mj=mj+mbig
    ma(inext)=mj
    ran3=mj*fac

  END FUNCTION ran3

END MODULE F90library
!
! definition of MPI constants
!
MODULE mpi_constants
  INCLUDE 'mpif.h'
  INTEGER, PUBLIC :: no_intervalls, iam, myloop_begin, myloop_end, &
                     numprocs, ierror

END MODULE mpi_constants

!
!   Main program starts here, standard Metropolis algorithm
!

PROGRAM ising2dim
  USE constants
  USE f90library
  USE mpi_constants
  INTEGER :: idum, n_spins, mcs, de, t, cycles, tstep
  INTEGER, ALLOCATABLE, DIMENSION(:,:) :: spin_matrix
  REAL(dp) :: w(-8:8), average(5), total_average(5)
  REAL(DP) :: initial_temp, final_temp, E, M, temp_step, temperature

  OPEN(UNIT=5,FILE='input.dat')
  CALL read_input(n_spins, mcs, initial_temp, final_temp, temp_step)
  ALLOCATE(spin_matrix(n_spins,n_spins))
  ! Find temperature step
  tstep = INT((final_temp-initial_temp)/temp_step); temperature = initial_temp
  ! random starting point
  idum = -1 
  OPEN(UNIT=6,FILE='mc_calc.dat')
  ! We assume that we start with the lowest temperature and that all
  ! spins are ordered. These spins are then updated with the latest value
  ! for every temperature in the loop over temperatures. The next temp step
  ! receives as input the previous spin matrix
  spin_matrix = 1
  ! initialize mpi
  !
  CALL MPI_INIT(ierror)
  CALL MPI_COMM_SIZE( MPI_COMM_WORLD, numprocs, ierror )
  CALL MPI_COMM_RANK( MPI_COMM_WORLD, iam, ierror )
  ! determine number of intervalls (sub matrices) which are summed 
  ! all processes
  ! myloop_begin gives the starting point on process iam
  ! myloop_end gives the end point for summation on process iam
  !
  no_intervalls = mcs/numprocs
  myloop_begin = iam*no_intervalls + 1
  myloop_end = (iam+1)*no_intervalls
  IF ( iam == numprocs-1 .and. myloop_end < mcs ) myloop_end = mcs
  DO t = 1, tstep
     !   initialise energy, magnetization and array for energy changes 
     e = 0.0_dp; m = 0.0_dp; w = 0.0_dp
     DO de =-8, 8, 4 
        w(de) = EXP(-de/temperature)
     ENDDO
     ! initialise array for expectation values
     total_average = 0.0_dp
     average = 0.0_dp
     CALL initialize(n_spins, temperature, spin_matrix, e, m)
     ! start Monte Carlo computation
     DO cycles = myloop_begin, myloop_end
        CALL Metropolis(n_spins, idum, spin_matrix, e, m, w)
        ! update expectation values
        average(1) = average(1) + e;     average(2) = average(2) + e*e
        average(3) = average(3)+m;     average(4) = average(4)+m*m 
        average(5) = average(5)+ABS(m)
     ENDDO
     ! collect all results
     CALL  mpi_reduce(average(:), total_average(:), 5, MPI_DOUBLE_PRECISION, MPI_SUM, 0, &
             MPI_COMM_WORLD, ierror)
     ! print results
     CALL output(n_spins, mcs, temperature, total_average)
     ! update temp
     temperature = temperature + temp_step
  ENDDO
  DEALLOCATE(spin_matrix)
  CALL MPI_FINALIZE(ierror)

END PROGRAM ising2dim
!
! read in input data
!
SUBROUTINE  read_input(n_spins, mcs, initial_temp, final_temp, temp_step)
  USE constants
  IMPLICIT NONE
  INTEGER, INTENT(inout) :: n_spins, mcs
  REAL(dp), INTENT(inout) :: initial_temp, final_temp, temp_step

  READ(5,*) ! number of cycles
  READ(5,*) mcs
  READ(5,*) ! size of spin lattice
  READ(5,*) n_spins
  READ(5,*) ! Initial temperature with dimension energy
  READ(5,*) initial_temp
  READ(5,*) !Final temperature with dimension energy
  READ(5,*) final_temp
  READ(5,*) !Temperature step with dimension energy
  READ(5,*) temp_step

END SUBROUTINE read_input
!
! function to initialise energy, spin matrix and magnetization
!
SUBROUTINE  initialize(n_spins, temperature, spin_matrix, E, M)
  USE constants
  IMPLICIT NONE
  INTEGER, INTENT(in) :: n_spins
  INTEGER, INTENT(inout) :: spin_matrix(n_spins, n_spins)
  REAL(dp), INTENT(in) :: temperature
  REAL(dp), INTENT(inout) :: E, M
  INTEGER :: x, y, right, left, up, down

  ! setup initial energy and magnetization
  DO y =1, n_spins
     DO x= 1, n_spins
        right = x+1 ; IF(x == n_spins  ) right = 1 
        left = x-1 ; IF(x == 1  ) left = n_spins 
        up = y+1 ; IF(y == n_spins  ) up = 1 
        down = y-1 ; IF(y == 1  ) down = n_spins 
        e= e-spin_matrix(x,y)*(spin_matrix(right,y)+&
             spin_matrix(left,y)+spin_matrix(x,up)+ &
             spin_matrix(x,down) )
        m = m + spin_matrix(x,y) 
     ENDDO
  ENDDO
  e = e*0.5_dp

END SUBROUTINE  initialize
!
!  Here we perform the sampling
!
SUBROUTINE metropolis(n_spins, idum, spin_matrix, e, m, w)
  USE constants
  USE f90library
  IMPLICIT NONE
  INTEGER, INTENT(in) :: n_spins, idum
  INTEGER :: x, y, ix, iy,  deltae, right, left, up, down
  INTEGER, INTENT(inout) :: spin_matrix(n_spins, n_spins)
  REAL(dp) :: e, m, w(-8:8)

  ! loop over all spins
  DO y =1, n_spins
     DO x= 1, n_spins
        ix = INT(ran1(idum)*n_spins)+1
        iy = INT(ran1(idum)*n_spins)+1
        right = ix+1 ; IF(ix == n_spins  ) right = 1 
        left = ix-1 ; IF(ix == 1  ) left = n_spins 
        up = iy+1 ; IF(iy == n_spins  ) up = 1 
        down = iy-1 ; IF(iy == 1  ) down = n_spins 
        deltae = 2*spin_matrix(ix,iy)*(spin_matrix(right,iy)+&
             spin_matrix(left,iy)+spin_matrix(ix,up)+ &
             spin_matrix(ix,down) )
        IF ( ran1(idum) <= w(deltae) ) THEN
           spin_matrix(ix,iy) = -spin_matrix(ix,iy)  !flip one spin and accept new spin config
           m = m+2*spin_matrix(ix,iy)
           e = e+deltaE
        ENDIF
     ENDDO
  ENDDO
END SUBROUTINE metropolis
!
!  Write out <E>, <M>, Heat capacity and susceptibility.
!
SUBROUTINE output(n_spins, mcs, temperature, average)
  USE constants
  USE mpi_constants
  IMPLICIT NONE
  INTEGER, INTENT(in) :: n_spins, mcs
  REAL(dp), INTENT(in) :: temperature, average(5)
  REAL(dp) :: norm, Eaverage, E2average, Maverage, M2average, Mabsaverage, &
       Evariance, Mvariance

  norm = 1.0_dp/mcs  ! divided by total number of cycles 
  Eaverage = average(1)*norm
  E2average = average(2)*norm
  Maverage = average(3)*norm
  M2average = average(4)*norm
  Mabsaverage = average(5)*norm
  ! all expectation values are per spin, divide by 1/n_spins/n_spins
  Evariance = (E2average- Eaverage*Eaverage)/n_spins/n_spins
  Mvariance = (M2average - Mabsaverage*Mabsaverage)/n_spins/n_spins
  IF (iam == 0 ) &
  WRITE(6,'(6F12.6)')temperature, Eaverage/n_spins/n_spins, Evariance/temperature/temperature, &
       Maverage/n_spins/n_spins, Mvariance/temperature, Mabsaverage/n_spins/n_spins

END SUBROUTINE output









