MODULE Shared
  USE, INTRINSIC :: ISO_Fortran_env, dp=>REAL64 !modern DOUBLE PRECISION
  IMPLICIT NONE
  INTEGER :: i_count, j, k, n_rejected, MCS, MCS_read
  INTEGER, PARAMETER :: N = 20 !total number of particles
  REAL(dp) :: v_cum, random, boltz, dE, beta, v_avg, K_avg, K_old
  REAL(dp), DIMENSION(:), ALLOCATABLE :: v_old, v_new
  REAL(dp), PARAMETER :: T = 10.0_dp
END MODULE Shared

PROGRAM Boltzmann
  USE Shared
  IMPLICIT NONE
  CALL INIT_RANDOM_SEED()
  WRITE(6, '(a29)', advance='no') "Number of Monte Carlo steps: "
  READ *, MCS_read
  MCS = MCS_read
  ALLOCATE(v_old(N + 1), v_new(N + 1))
  beta = 1.0_dp / T !Boltzmann's constant = 1 in this simulation
  v_old = T !initialize entire array
  v_cum = 0.0_dp
  initial: DO j = 1, N !calculate the initial energy of the system
    v_cum = v_cum + v_old(j) * v_old(j)
  END DO initial
  K_old = 0.5_dp * v_cum !initial total energy of the system; PE is zero (ideal gas)
  n_rejected = 0
  K_avg = 0.0_dp
  v_avg = 0.0_dp
  v_cum = 0.0_dp
  Monte_Carlo: DO k = 1, MCS
    CALL RANDOM_NUMBER(random) !select a particle at random
    i_count = INT(FLOAT(N) * random) + 1
    CALL RANDOM_NUMBER(random) !change velocity of the randomly selected particle
    v_new(i_count) = v_old(i_count) + (1.0_dp - 2.0_dp * random) * 0.5_dp
    dE = 0.5_dp * v_new(i_count)**2 - 0.5_dp * v_old(i_count)**2 !calculate the change in energy
    IF (dE <= 0.0_dp) THEN !condition for an accepted microstate
      v_old(i_count) = v_new(i_count) !update velocity
      K_old = K_old + dE !step accepted; update energy
    ELSE
      boltz = EXP(-dE * beta)
      CALL RANDOM_NUMBER(random)
      IF (boltz >= random) THEN !accept step
        v_old(i_count) = v_new(i_count)
        K_old = K_old + dE
      ELSE !reject step; keep current energy
        n_rejected = n_rejected + 1
      END IF
    END IF
    K_avg = K_avg + K_old / FLOAT(MCS)
    v_avg = v_avg + v_old(i_count) / FLOAT(MCS)
  END DO Monte_Carlo
  DEALLOCATE(v_old, v_new)
  WRITE(6, '(a8)', advance='no') "Mean KE:"
  WRITE(6, *) K_avg
  WRITE(6, '(a14)', advance='no') "Mean velocity:"
  WRITE(6, *) v_avg
  WRITE(6, '(a37)', advance='no') "Number of rejected Monte Carlo steps:"
  WRITE(6, *) n_rejected
END PROGRAM Boltzmann

SUBROUTINE INIT_RANDOM_SEED()
  USE ISO_Fortran_env, ONLY: INT64
  IMPLICIT NONE
  INTEGER, ALLOCATABLE :: seed(:)
  INTEGER :: i, n, un, istat, dt(8), pid
  INTEGER(INT64) :: t
  CALL RANDOM_SEED(size = n)
  ALLOCATE(seed(n))
  OPEN(newunit=un, file='/dev/urandom', access='stream', status='old', action='read', form='unformatted', iostat=istat)
  IF (istat == 0) THEN
    READ(un) seed
    CLOSE(un)
  ELSE
    CALL SYSTEM_CLOCK(t)
    IF (t == 0) THEN
      CALL DATE_AND_TIME(values = dt)
      t = (dt(1) - 1970) * 365_INT64 * 24 * 60 * 60 * 1000 + dt(2) * 31_INT64 * 24 * 60 * 60 * 1000 + dt(3) * 24_INT64 * 60 * 60   &
        * 1000 + dt(5) * 60 * 60 * 1000 + dt(6) * 60 * 1000 + dt(7) * 1000 + dt(8)
    END IF
    pid = GETPID()
    t = IEOR(t, INT(pid, KIND(t)))
    DO i = 1, n
      seed(i) = lcg(t)
    END DO
  END IF
  CALL RANDOM_SEED(put = seed)
  DEALLOCATE(seed)
CONTAINS
  FUNCTION lcg(s)
    INTEGER :: lcg
    INTEGER(INT64) :: s
    IF (s == 0) THEN
      s = 104729
    ELSE
      s = MOD(s, 4294967296_INT64)
    END IF
    s = MOD(s * 279470273_INT64, 4294967291_INT64)
    lcg = INT(MOD(s, INT(HUGE(0), INT64)), KIND(0))
  END FUNCTION lcg
END SUBROUTINE INIT_RANDOM_SEED