Merge pull request #18 from Allen-Tildesley/minor-dev

Minor dev

README.md

@@ -59,8 +59,9 @@ A few examples consist of individual routines or modules,
 rather than working programs,
 so there is no need to build them.
 
-The build process for the Fortran examples has been tested using SCons v4.2.0
+The build process for the Fortran examples has been tested using SCons v4.4.0
 (and some earlier versions back to v2.5.1 with minor changes to the `SConstruct` file).
+
 If you don't like using SCons, or can't get it to work,
 it is not difficult to compile the programs using other methods.
 Bear in mind that, with Fortran, it is usually essential to compile any
@@ -77,20 +78,8 @@ it is advisable to __compile each example in its own build directory__
 (which is what the `SConstruct` file is configured to do)
 or to delete all intermediate files before each individual compilation.
 
-We have used gfortran v11.2.0 (and earlier versions back to v6.3) for testing,
+We have used gfortran v12.1.0 (and earlier versions back to v6.3) for testing,
 but have attempted to stick to code which conforms to the Fortran 2008 standard.
-In gfortran v6, calling the intrinsic `random_seed()` function would
-generate the same sequence of random numbers every time.
-In a few examples it is important to generate different sequences each time the program is run,
-and for this purpose a subroutine `init_random_seed` was provided in the file `maths_module.f90`;
-this routine, however, included some GNU extensions.
-With gfortran v7, the random number generator was changed,
-and the intrinsic `random_seed()` function now behaves in the desired fashion.
-Therefore, the GNU-specific code has been transferred to a separate file, `gnu_v6_init_random_seed.f90`,
-which is not included in the build process,
-but may be of interest to those still using gfortran v6.
-The current `init_random_seed` routine simply calls `random_seed()`.
-__You should check the behaviour of the random number generator on your own system.__
 Note that, by default, we do not select any optimization options in compilation.
 If you are using a different compiler,
 then the compiler and linker options in the `SConstruct` file will most likely need changing.
@@ -100,6 +89,40 @@ Unfortunately, due to the enormous variety of computing platforms and compilers,
 __we cannot offer more specific advice on the build process.__
 
 The Python versions do not require building, they are simply run through the Python interpreter.
+They have been tested with Python 3.10.6
+(and before that, versions back to 3.6.0).
+
+## Random number generators
+
+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.
+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`
+contained in the file `gnu_v6_init_random_seed.f90`;
+but this file is now retained only for historical reference
+and is not included in the build process.
+
+The Python examples use, for simplicity,
+NumPy convenience routines such as
+`random.seed()` and `random.rand()`
+for the same purpose.
+Since NumPy v1.17.0,
+this random number generator is classified as "legacy":
+a different, and more flexible, approach is now provided and recommended
+in the NumPy documentation.
+Nonetheless, the legacy generator continues to be supported within NumPy,
+for backwards compatibility,
+and so we continue to use it in these examples.
+
+We do not recommend the above choices (using built-in and/or legacy random number generators) for production work.
+In any case, quite generally,
+__you should check the behaviour of the random number generator on your own system.__
 
 ## Reporting errors
 If you spot an error in these program files, or the accompanying documentation,

cluster.f90

@@ -114,7 +114,7 @@ PROGRAM cluster
 
      IF ( ALL ( done > 0 ) ) EXIT ! Loop until all done
 
-     i = MINLOC ( done, dim = 1 ) ! Find first zero (FINDLOC is not implemented in gfortran at the time of writing)
+     i = FINDLOC ( done, 0, dim = 1 ) ! Find first zero (FINDLOC was implemented in gfortran v9)
      cluster_id = cluster_id + 1
      WRITE ( unit=output_unit, fmt='(a,i5,a)', advance='no' ) 'Cluster ', cluster_id, ' = '
      j = i

error_calc.f90

@@ -32,7 +32,7 @@ PROGRAM error_calc
   ! and AD Baczewski and SD Bond J Chem Phys 139 044107 (2013)
 
   USE, INTRINSIC :: iso_fortran_env, ONLY : input_unit, output_unit, error_unit, iostat_end, iostat_eor
-  USE               maths_module,    ONLY : random_normal, init_random_seed
+  USE               maths_module,    ONLY : random_normal
 
   IMPLICIT NONE
 
@@ -105,7 +105,7 @@ PROGRAM error_calc
      d = x * ( d1 + x * ( d2 + x * ( d3 + x * d4 ) ) )
   END IF
   b      = SQRT ( b )
-  b      = b * sqrt ( kappa/2.0 ) ! alpha in B&B paper  
+  b      = b * SQRT ( kappa/2.0 ) ! alpha in B&B paper  
   stddev = SQRT(2.0*variance)     ! NB stddev of random forces, not data
 
   ! For this process, the results of interest can be calculated exactly
@@ -123,7 +123,7 @@ PROGRAM error_calc
 
   ! Data generation
 
-  CALL init_random_seed
+  CALL RANDOM_SEED()
 
   ! 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

ewald_module.f90

@@ -210,8 +210,8 @@ FUNCTION pot_k_pm_ewald ( nk, n, r, q, kappa ) RESULT ( pot )
 
     INTEGER                             :: ix, iy, iz
     REAL                                :: kx, ky, kz, fx, fxy, fxyz, g, k_sq, dr
-    real, DIMENSION(:,:,:), ALLOCATABLE :: rho_sq ! Square modulus of charge density (0:sc-1,0:sc-1,0:sc-1)
-    real, dimension(:),     allocatable :: kmesh  ! k-values in wraparound convention (0:sc-1)
+    REAL, DIMENSION(:,:,:), ALLOCATABLE :: rho_sq ! Square modulus of charge density (0:sc-1,0:sc-1,0:sc-1)
+    REAL, DIMENSION(:),     ALLOCATABLE :: kmesh  ! k-values in wraparound convention (0:sc-1)
 
     REAL, PARAMETER    :: pi = 4.0*ATAN(1.0), twopi = 2.0*pi, fourpi = 4.0*pi, dk = twopi
 
@@ -220,7 +220,7 @@ FUNCTION pot_k_pm_ewald ( nk, n, r, q, kappa ) RESULT ( pot )
     ! Use nk to determine mesh dimension
     sc = 2*nk
     ALLOCATE ( fft_inp(0:sc-1,0:sc-1,0:sc-1), fft_out(0:sc-1,0:sc-1,0:sc-1) )
-    allocate ( rho_sq(0:sc-1,0:sc-1,0:sc-1), kmesh(0:sc-1) )
+    ALLOCATE ( rho_sq(0:sc-1,0:sc-1,0:sc-1), kmesh(0:sc-1) )
 
     ! Assign charge density to complex array ready for FFT
     ! Assume r in unit box with range (-0.5,0.5)
@@ -231,10 +231,10 @@ FUNCTION pot_k_pm_ewald ( nk, n, r, q, kappa ) RESULT ( pot )
     CALL fftw_destroy_plan ( fft_plan )                                                      ! Release plan
 
     fft_out = fft_out / REAL(sc**3)           ! Incorporate scaling by number of grid points
-    rho_sq  = real ( fft_out*conjg(fft_out) ) ! Convert to square modulus of charge density
+    rho_sq  = REAL ( fft_out*CONJG(fft_out) ) ! Convert to square modulus of charge density
 
-    kmesh      = [ ( real(ix)*dk, ix = 0, sc-1 ) ] ! Set up k-components
-    kmesh(nk:) = kmesh(nk:) - real(sc)*dk          ! in wraparound convention
+    kmesh      = [ ( REAL(ix)*dk, ix = 0, sc-1 ) ] ! Set up k-components
+    kmesh(nk:) = kmesh(nk:) - REAL(sc)*dk          ! in wraparound convention
     pot = 0.0
 
     ! Triple loop over xyz grid points (uses wraparound indexing)

gnu_v6_init_random_seed.f90

@@ -30,8 +30,8 @@
 ! For this version, calling RANDOM_SEED() initializes the random number generator
 ! with the same random seed to a default state, which may result in the same sequence
 ! being generated every time. The routines below are intended to generate different
-! sequences on different calls. If you have this compiler version, you may like to replace the
-! routine init_random_seed in maths_module.f90 with the following two routines.
+! sequences on different calls. If you have this compiler version, you may like to replace any
+! call of RANDOM_SEED() with the following two routines.
 
 ! In gfortran v7 the random number generator was changed.
 ! Now, calling RANDOM_SEED() initializes the random number generator with random data

maths_module.f90

@@ -31,7 +31,7 @@ MODULE maths_module
   PRIVATE
 
   ! Public random number routines
-  PUBLIC :: init_random_seed, random_integer, random_normal, random_normals, pick
+  PUBLIC :: random_integer, random_normal, random_normals, pick
   PUBLIC :: random_vector
   PUBLIC :: random_vector_1, random_vector_2, random_vector_3
   PUBLIC :: random_perpendicular_vector
@@ -84,27 +84,6 @@ MODULE maths_module
 
   ! Routines associated with random number generation
 
-  SUBROUTINE init_random_seed
-    IMPLICIT NONE
-
-    ! Initializes random number generator differently each time
-
-    ! Prior to gfortran v7, calling the intrinsic RANDOM_SEED() routine initializes the
-    ! random number generator with the same random seed to a default state, 
-    ! which may result in the same sequence being generated every time.
-    ! If you have gfortran v6, you may like to replace this routine init_random_seed 
-    ! with the contents of file gnu_v6_init_random_seed.f90
-
-    ! In gfortran v7 the random number generator was changed.
-    ! Now, calling RANDOM_SEED() initializes the random number generator with random data
-    ! retrieved from the operating system. Various other compilers behave the same way.
-    ! We assume that this applies here.
-
-    ! YOU SHOULD INVESTIGATE THE BEHAVIOUR FOR YOUR OWN COMPILER AND MACHINE IMPLEMENTATION 
-
-    CALL RANDOM_SEED()
-  END SUBROUTINE init_random_seed
-
   FUNCTION random_integer ( k1, k2 ) RESULT ( k )
     IMPLICIT NONE
     INTEGER             :: k      ! Returns a uniformly distributed random integer

mc_nvt_lj_re.f90

@@ -62,7 +62,7 @@ PROGRAM mc_nvt_lj_re
   USE, INTRINSIC :: iso_fortran_env,  ONLY : input_unit, output_unit, error_unit, iostat_end, iostat_eor
   USE               config_io_module, ONLY : read_cnf_atoms, write_cnf_atoms
   USE               averages_module,  ONLY : run_begin, run_end, blk_begin, blk_end, blk_add
-  USE               maths_module,     ONLY : init_random_seed, metropolis, random_translate_vector
+  USE               maths_module,     ONLY : metropolis, random_translate_vector
   USE               mc_module,        ONLY : introduction, conclusion, allocate_arrays, deallocate_arrays, &
        &                                     potential_1, potential, move, n, r, potential_type
   USE               mpi
@@ -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 init_random_seed () ! Initialize random number generator (hopefully differently on each process)
+  CALL RANDOM_SEED () ! 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
 

python_examples/GUIDE.md

@@ -40,10 +40,10 @@ the main aim of these examples is to illustrate ideas in the text,
 not to provide programs for practical use.
 The reader should feel free to experiment with ways to make the programs run faster!
 
-In the past few years, the community has been making the transition from Python 2 to Python 3.
+In the past few years, the community has made the transition from Python 2 to Python 3.
 There are some incompatibilities between the two,
 and since a choice had to be made,
-we have settled on __Python 3__ for these examples.
+we settled on __Python 3__ for these examples.
 We indicate this by the string
 ```
 #!/usr/bin/env python3
@@ -56,7 +56,6 @@ or just `python sample_mean.py`,
 depending on your particular installation of Python.
 (Of course, in most cases, it will also be necessary to input data to the program.)
 The examples _will not work_ with Python 2!
-They have been tested with Python 3.7.5 (and before that, 3.6.0).
 For an introduction to the differences between Python 2 and Python 3,
 see the [What's New in Python 3.0](https://docs.python.org/3/whatsnew/3.0.html "What's New in Python 3.0") page.
 The most obvious changes are

python_examples/README.md

@@ -30,5 +30,6 @@ and are written in Python3 using the NumPy and SciPy libraries.
 The [User Guide](./GUIDE.md) contains some comments on the Python language,
 some notes to assist in running the programs,
 and some typical results.
-We have used Python 3.9.5 (and some earlier versions back to 3.6.0) for testing.
 The Python versions do not require building, they are simply run through the Python interpreter.
+They have been tested with Python 3.10.6
+(and before that, versions back to 3.6.0).

t_tensor.f90

@@ -44,7 +44,7 @@ PROGRAM t_tensor
 
   USE, INTRINSIC :: iso_fortran_env, ONLY : input_unit, output_unit, error_unit, iostat_end, iostat_eor
 
-  USE maths_module,    ONLY : init_random_seed, random_vector, outer_product, cross_product
+  USE maths_module, ONLY : random_vector, outer_product, cross_product
 
   IMPLICIT NONE
 
@@ -65,7 +65,7 @@ PROGRAM t_tensor
   WRITE ( unit=output_unit, fmt='(a)' ) 'using T-tensors and Euler angles'
 
   ! Initialize random number generator                          
-  CALL init_random_seed
+  CALL RANDOM_SEED ()
 
   ! Default parameters
   d_min     = 0.5 ! Minimum separation

test_pot_atom.f90

@@ -27,7 +27,7 @@ PROGRAM test_pot_atom
 
   USE, INTRINSIC :: iso_fortran_env, ONLY : input_unit, output_unit, error_unit, iostat_end, iostat_eor
   USE               test_pot_module, ONLY : n, force
-  USE               maths_module,    ONLY : init_random_seed, cross_product
+  USE               maths_module,    ONLY : cross_product
 
   IMPLICIT NONE
 
@@ -47,7 +47,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 init_random_seed
+  CALL RANDOM_SEED ()
 
   ! Default values: any of the following could be empirically adjusted
   delta   = 1.e-5 ! Small displacement

test_pot_linear.f90

@@ -27,7 +27,7 @@ PROGRAM test_pot_linear
 
   USE, INTRINSIC :: iso_fortran_env, ONLY : input_unit, output_unit, error_unit, iostat_end, iostat_eor
   USE               test_pot_module, ONLY : n, force
-  USE               maths_module,    ONLY : init_random_seed, rotate_vector, cross_product
+  USE               maths_module,    ONLY : rotate_vector, cross_product
 
   IMPLICIT NONE
 
@@ -50,7 +50,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 init_random_seed
+  CALL RANDOM_SEED ()
 
   ! Default values: any of the following could be empirically adjusted
   delta   = 1.e-5 ! Small displacement