Replaced random_seed by random_init, fixes #21

README.md

@@ -96,11 +96,12 @@ They have been tested with Python 3.12.2 and NumPy 1.26.4
 
 The Fortran examples use, for simplicity,
 the built-in intrinsic subroutines
-`random_seed` and `random_number` respectively to
-initialize and generate sequences of random numbers.
-From gfortran v7 onwards,
-calling `random_seed()` generates different, non-reproducible, sequences each time,
-and the examples assume this behaviour.
+`random_init` and `random_number` respectively to
+initialize and generate different, non-reproducible, sequences of random numbers every time.
+In Fortran 2018 `random_init`  was introduced into the standard for this purpose.
+Previously, we would call `random_seed()` with no argument,
+which served the same function in recent versions of gfortran (v7 and above),
+but was not part of the standard (and hence, potentially, compiler dependent).
 Prior to gfortran v7,
 it was necessary to do something more complicated to generate different sequences each time,
 as exemplified by the routine `init_random_seed`

bd_nvt_lj.f90

@@ -77,7 +77,7 @@ PROGRAM bd_nvt_lj
   WRITE ( unit=output_unit, fmt='(a)' ) 'Particle mass m=1 throughout'
   CALL introduction
 
-  CALL RANDOM_SEED () ! Initialize random number generator
+  CALL RANDOM_INIT ( .FALSE., .TRUE. ) ! Initialize random number generator
 
   ! Set sensible default run parameters for testing
   nblock      = 10

corfun.f90

@@ -84,7 +84,7 @@ PROGRAM corfun
   WRITE( unit=output_unit, fmt='(a)' ) 'Illustrates methods for calculating time correlation functions'
   WRITE( unit=output_unit, fmt='(a)' ) 'using synthetic data from a generalized Langevin equation'
 
-  CALL RANDOM_SEED () ! Initialize random number generator
+  CALL RANDOM_INIT ( .FALSE., .TRUE. ) ! Initialize random number generator
 
   ! Example default values
   ! Agreement (to numerical precision) of direct and FFT methods is expected if origin_interval=1

diffusion_test.f90

@@ -97,7 +97,7 @@ PROGRAM diffusion_test
   ALLOCATE ( r(3,n), v(3,n), zeta(3,n) )
 
   ! Set random positions
-  CALL RANDOM_SEED ()
+  CALL RANDOM_INIT ( .FALSE., .TRUE. ) ! Initialize random number generator
   CALL RANDOM_NUMBER ( r )
   r = r - 0.5 ! Now in range (-1/2,1/2)
   r = r * box ! Now in range (-box/2,box/2)

dpd.f90

@@ -83,7 +83,7 @@ PROGRAM dpd
   WRITE ( unit=output_unit, fmt='(a)' ) 'Particle mass=1 and cutoff=1 throughout'
   CALL introduction
 
-  CALL RANDOM_SEED () ! Initialize random number generator
+  CALL RANDOM_INIT ( .FALSE., .TRUE. ) ! Initialize random number generator
 
   ! Set sensible default run parameters for testing
   nblock      = 10

error_calc.f90

@@ -114,7 +114,7 @@ PROGRAM error_calc
 
   ! Data generation
 
-  CALL RANDOM_SEED()
+  CALL RANDOM_INIT ( .FALSE., .TRUE. ) ! Initialize random number generator
 
   ! For comparison, we do n_repeat independent runs and estimate the error in run averages directly from these
   ! This is to give an empirical idea of the distribution from which the run average is sampled

hit_and_miss.f90

@@ -35,7 +35,7 @@ PROGRAM hit_and_miss
   REAL,               PARAMETER :: v_0 = PRODUCT(r_0)
   INTEGER                       :: tau, tau_shot, tau_hit
 
-  CALL RANDOM_SEED()
+  CALL RANDOM_INIT ( .FALSE., .TRUE. ) ! Initialize random number generator
   tau_hit  = 0
   tau_shot = 1000000
 

initialize.f90

@@ -68,7 +68,7 @@ PROGRAM initialize
   WRITE ( unit=output_unit, fmt='(a)' ) 'Options for molecules are "atom", "linear", "nonlinear", "chain"'
   WRITE ( unit=output_unit, fmt='(a)' ) 'Particle mass m=1 throughout'
 
-  CALL RANDOM_SEED()
+  CALL RANDOM_INIT ( .FALSE., .TRUE. ) ! Initialize random number generator
 
   ! Set default parameters
   n           = 0        ! nc takes precedence unless n is explicitly specified

mc_chain_nvt_cbmc_lj.f90

@@ -71,7 +71,7 @@ PROGRAM mc_chain_nvt_cbmc_lj
   WRITE ( unit=output_unit, fmt='(a)' ) 'Simulation uses full nonbonded potential (no cutoff)'
   CALL introduction
 
-  CALL RANDOM_SEED () ! Initialize random number generator
+  CALL RANDOM_INIT ( .FALSE., .TRUE. ) ! Initialize random number generator
 
   ! Set sensible default run parameters for testing
   nblock      = 10     ! Number of blocks

mc_chain_nvt_sw.f90

@@ -89,7 +89,7 @@ PROGRAM mc_chain_nvt_sw
   WRITE ( unit=output_unit, fmt='(a)' ) 'Monte Carlo, constant-NVT ensemble, chain molecule, square wells'
   CALL introduction
 
-  CALL RANDOM_SEED () ! Initialize random number generator
+  CALL RANDOM_INIT ( .FALSE., .TRUE. ) ! Initialize random number generator
 
   ! Set sensible default run parameters for testing
   nblock         = 10     ! Number of blocks

mc_chain_wl_sw.f90

@@ -103,7 +103,7 @@ PROGRAM mc_chain_wl_sw
   WRITE ( unit=output_unit, fmt='(a)' ) 'Monte Carlo, Wang-Landau method, chain molecule, square wells'
   CALL introduction
 
-  CALL RANDOM_SEED () ! Initialize random number generator
+  CALL RANDOM_INIT ( .FALSE., .TRUE. ) ! Initialize random number generator
 
   ! Set sensible default run parameters for testing
   nstage         = 20    ! 2**(-20) = approx 10**(-6) for smallest modification factor

mc_gibbs_lj.f90

@@ -90,7 +90,7 @@ PROGRAM mc_gibbs_lj
   WRITE( unit=output_unit, fmt='(a)' ) 'Monte Carlo, Gibbs ensemble'
   CALL introduction
 
-  CALL RANDOM_SEED () ! Initialize random number generator
+  CALL RANDOM_INIT ( .FALSE., .TRUE. ) ! Initialize random number generator
 
   ! Set sensible default run parameters for testing
   nblock      = 10

mc_npt_hs.f90

@@ -71,7 +71,7 @@ PROGRAM mc_npt_hs
   WRITE( unit=output_unit, fmt='(a)' ) 'Monte Carlo, constant-NPT'
   CALL introduction
 
-  CALL RANDOM_SEED () ! Initialize random number generator
+  CALL RANDOM_INIT ( .FALSE., .TRUE. ) ! Initialize random number generator
 
   ! Set sensible default run parameters for testing
   nblock   = 10

mc_npt_lj.f90

@@ -85,7 +85,7 @@ PROGRAM mc_npt_lj
   WRITE ( unit=output_unit, fmt='(a)' ) 'Simulation uses cut (but not shifted) potential'
   CALL introduction
 
-  CALL RANDOM_SEED () ! Initialize random number generator
+  CALL RANDOM_INIT ( .FALSE., .TRUE. ) ! Initialize random number generator
 
   ! Set sensible default run parameters for testing
   nblock      = 10

mc_npt_sc.f90

@@ -70,7 +70,7 @@ PROGRAM mc_npt_sc
   WRITE( unit=output_unit, fmt='(a)' ) 'Monte Carlo, constant-NPT, hard linear molecules'
   CALL introduction
 
-  CALL RANDOM_SEED () ! Initialize random number generator
+  CALL RANDOM_INIT ( .FALSE., .TRUE. ) ! Initialize random number generator
 
   ! Set sensible default run parameters for testing
   nblock   = 10

mc_nvt_hs.f90

@@ -68,7 +68,7 @@ PROGRAM mc_nvt_hs
   WRITE( unit=output_unit, fmt='(a)' ) 'Monte Carlo, constant-NVT'
   CALL introduction
 
-  CALL RANDOM_SEED () ! Initialize random number generator
+  CALL RANDOM_INIT ( .FALSE., .TRUE. ) ! Initialize random number generator
 
   ! Set sensible default run parameters for testing
   nblock  = 10

mc_nvt_lj.f90

@@ -75,7 +75,7 @@ PROGRAM mc_nvt_lj
   WRITE ( unit=output_unit, fmt='(a)' ) 'Simulation uses cut (but not shifted) potential'
   CALL introduction
 
-  CALL RANDOM_SEED () ! Initialize random number generator
+  CALL RANDOM_INIT ( .FALSE., .TRUE. ) ! Initialize random number generator
 
   ! Set sensible default run parameters for testing
   nblock      = 10

mc_nvt_lj_re.f90

@@ -123,7 +123,7 @@ PROGRAM mc_nvt_lj_re
   WRITE( unit=output_unit, fmt='(a,t40,i15)') 'This is process rank',   m
   WRITE( unit=output_unit, fmt='(a,t40,i15)') 'Number of processes is', nproc
 
-  CALL RANDOM_SEED () ! Initialize random number generator (hopefully differently on each process)
+  CALL RANDOM_INIT ( .FALSE., .TRUE. ) ! Initialize random number generator (hopefully differently on each process)
   CALL RANDOM_NUMBER ( zeta )
   WRITE( unit=output_unit, fmt='(a,t40,f15.6)') 'Random # (different for each process?)', zeta
 

mc_nvt_poly_lj.f90

@@ -77,7 +77,7 @@ PROGRAM mc_nvt_poly_lj
   WRITE ( unit=output_unit, fmt='(a)' ) 'Monte Carlo, constant-NVT ensemble, polyatomic molecule'
   CALL introduction
 
-  CALL RANDOM_SEED () ! Initialize random number generator
+  CALL RANDOM_INIT ( .FALSE., .TRUE. ) ! Initialize random number generator
 
   ! Set sensible default run parameters for testing
   nblock      = 10

mc_nvt_sc.f90

@@ -69,7 +69,7 @@ PROGRAM mc_nvt_sc
   WRITE( unit=output_unit, fmt='(a)' ) 'Monte Carlo, constant-NVT, hard linear molecules'
   CALL introduction
 
-  CALL RANDOM_SEED () ! Initialize random number generator
+  CALL RANDOM_INIT ( .FALSE., .TRUE. ) ! Initialize random number generator
 
   ! Set sensible default run parameters for testing
   nblock  = 10

mc_zvt_lj.f90

@@ -79,7 +79,7 @@ PROGRAM mc_zvt_lj
   WRITE( unit=output_unit, fmt='(a)' ) 'Monte Carlo, constant-zVT ensemble'
   CALL introduction
 
-  CALL RANDOM_SEED () ! Initialize random number generator
+  CALL RANDOM_INIT ( .FALSE., .TRUE. ) ! Initialize random number generator
 
   ! Set sensible default run parameters for testing
   nblock      = 10

md_npt_lj.f90

@@ -138,7 +138,7 @@ PROGRAM md_npt_lj
   v(:,:) = v(:,:) - SPREAD ( vcm(:), dim = 2, ncopies = n ) ! Set COM velocity to zero
 
   ! Initial values of thermal variables
-  CALL RANDOM_SEED()
+  CALL RANDOM_INIT ( .FALSE., .TRUE. ) ! Initialize random number generator
   g    = REAL ( 3*(n-1) )
   q    = temperature * tau**2   
   q(1) = g * temperature * tau**2

md_nvt_lj.f90

@@ -124,7 +124,7 @@ PROGRAM md_nvt_lj
   v(:,:) = v(:,:) - SPREAD ( vcm(:), dim = 2, ncopies = n ) ! Set COM velocity to zero
 
   ! Initial values of thermal variables
-  CALL RANDOM_SEED()
+  CALL RANDOM_INIT ( .FALSE., .TRUE. ) ! Initialize random number generator
   g    = REAL ( 3*(n-1) )
   q    = temperature * tau**2   
   q(1) = g * temperature * tau**2

md_nvt_poly_lj.f90

@@ -79,7 +79,7 @@ PROGRAM md_nvt_poly_lj
   WRITE ( unit=output_unit, fmt='(a)' ) 'Molecular mass=1 throughout'
 
   CALL introduction
-  CALL RANDOM_SEED () ! Initialize random number generator
+  CALL RANDOM_INIT ( .FALSE., .TRUE. ) ! Initialize random number generator
 
   ! Set sensible default run parameters for testing
   nblock      = 10

mesh.f90

@@ -77,7 +77,7 @@ PROGRAM mesh
 
   ! For illustration we choose random charge positions with coordinates in range (0,1)
   ! In a real application, we would convert positions into this range
-  CALL RANDOM_SEED() ! same random number sequence every time
+  CALL RANDOM_INIT ( .FALSE., .TRUE. )  ! Initialize random number generator
   CALL RANDOM_NUMBER ( r )
 
   ! For illustration we choose +1 and -1 charges, alternately

qmc_pi_lj.f90

@@ -85,7 +85,7 @@ PROGRAM qmc_pi_lj
   WRITE ( unit=output_unit, fmt='(a)' ) 'Simulation uses cut (but not shifted) potential'
   CALL introduction
 
-  CALL RANDOM_SEED () ! Initialize random number generator
+  CALL RANDOM_INIT ( .FALSE., .TRUE. ) ! Initialize random number generator
 
   ! Set sensible default run parameters for testing
   nblock      = 10

qmc_pi_sho.f90

@@ -68,7 +68,7 @@ PROGRAM qmc_pi_sho
   WRITE ( unit=output_unit, fmt='(a)' ) 'Path Integral Monte Carlo simulation of a quantum oscillator'
   WRITE ( unit=output_unit, fmt='(a)' ) 'Results in atomic units'
 
-  CALL RANDOM_SEED() ! Initialize random number generator
+  CALL RANDOM_INIT ( .FALSE., .TRUE. ) ! Initialize random number generator
 
   ! Set sensible default run parameters for testing
   p           = 8

qmc_walk_sho.f90

@@ -111,7 +111,7 @@ PROGRAM qmc_walk_sho
   ALLOCATE ( x(n_max), v(n_max), replica(n_max), alive(n_max) )
   ALLOCATE ( psi(n_bin), bin(n_bin) )
 
-  CALL RANDOM_SEED()
+  CALL RANDOM_INIT ( .FALSE., .TRUE. ) ! Initialize random number generator
 
   ! Set up an initial delta distribution of walkers at origin
   n = n_target

sample_mean.f90

@@ -35,7 +35,7 @@ PROGRAM sample_mean
   REAL,               PARAMETER :: a_0 = PRODUCT(r_0)
   INTEGER                       :: tau, tau_max
 
-  CALL RANDOM_SEED()
+  CALL RANDOM_INIT ( .FALSE., .TRUE. ) ! Initialize random number generator
   tau_max = 1000000
 
   f = 0.0

smc_nvt_lj.f90

@@ -81,7 +81,7 @@ PROGRAM smc_nvt_lj
   WRITE ( unit=output_unit, fmt='(a)' ) 'Smart Monte Carlo, constant-NVT ensemble'
   CALL introduction
 
-  CALL RANDOM_SEED () ! Initialize random number generator
+  CALL RANDOM_INIT ( .FALSE., .TRUE. ) ! Initialize random number generator
 
   ! Set sensible default run parameters for testing
   nblock      = 10

t_tensor.f90

@@ -64,8 +64,7 @@ PROGRAM t_tensor
   WRITE ( unit=output_unit, fmt='(a)' ) 'Calculation of electrostatic interactions between linear molecules' 
   WRITE ( unit=output_unit, fmt='(a)' ) 'using T-tensors and Euler angles'
 
-  ! Initialize random number generator                          
-  CALL RANDOM_SEED ()
+  CALL RANDOM_INIT ( .FALSE., .TRUE. ) ! Initialize random number generator
 
   ! Default parameters
   d_min     = 0.5 ! Minimum separation

test_pot_atom.f90

@@ -46,8 +46,7 @@ PROGRAM test_pot_atom
 
   NAMELIST /nml/ delta, d_min, d_max, pot_max, ntry, npos
 
-  ! Initialize random number generator (hopefully different every time!)
-  CALL RANDOM_SEED ()
+  CALL RANDOM_INIT ( .FALSE., .TRUE. ) ! Initialize random number generator
 
   ! Default values: any of the following could be empirically adjusted
   delta   = 1.e-5 ! Small displacement

test_pot_linear.f90

@@ -49,8 +49,7 @@ PROGRAM test_pot_linear
 
   NAMELIST /nml/ delta, d_min, d_max, pot_max, ntry, npos
 
-  ! Initialize random number generator (hopefully different every time!)
-  CALL RANDOM_SEED ()
+  CALL RANDOM_INIT ( .FALSE., .TRUE. ) ! Initialize random number generator
 
   ! Default values: any of the following could be empirically adjusted
   delta   = 1.e-5 ! Small displacement