Minor improvement reallocating arrays when N increases in mc_zvt_lj. …

…This closes #15.

GUIDE.md

@@ -296,9 +296,6 @@ The [Thol et al (2016)](https://doi.org/10.1063/1.4945000) LRC-corrected value t
 For `mc_zvt_lj` the box length was _L_=7σ; for `mc_zvt_lj_ll` _L_=10.5σ.
 Acceptance rate of creation/destruction moves is quite small, at about 0.3%.
 For other state points see below.
-These programs could be improved to use array reallocation (available in Fortran)
-to make them more resilient against large increases in the number of particles.
-For simplicity we have not included this feature.
 
 ♯ The `mc_nvt_lj_re` program was run for four temperatures, see below for details.
 

link_list_module.f90

@@ -31,7 +31,7 @@ MODULE link_list_module
   PRIVATE
 
   ! Public routines
-  PUBLIC :: initialize_list, finalize_list, make_list, check_list
+  PUBLIC :: initialize_list, extend_list, finalize_list, make_list, check_list
   PUBLIC :: move_in_list, create_in_list, destroy_in_list, c_index, neighbours
 
   ! Public (protected) data
@@ -40,6 +40,10 @@ MODULE link_list_module
   INTEGER, DIMENSION(:),     ALLOCATABLE, PROTECTED, SAVE, PUBLIC :: list   ! list(n)
   INTEGER, DIMENSION(:,:),   ALLOCATABLE, PROTECTED, SAVE, PUBLIC :: c      ! c(3,n) 3D cell index of each atom
 
+  ! Private data
+  INTEGER, DIMENSION(:),     ALLOCATABLE :: tmp_list ! for dynamic reallocation of list(n)
+  INTEGER, DIMENSION(:,:),   ALLOCATABLE :: tmp_c    ! for dynamic reallocation of c(3,n)
+
 CONTAINS
 
   SUBROUTINE initialize_list ( n, r_cut_box ) ! Routine to allocate list arrays
@@ -61,6 +65,33 @@ SUBROUTINE initialize_list ( n, r_cut_box ) ! Routine to allocate list arrays
 
   END SUBROUTINE initialize_list
 
+  SUBROUTINE extend_list ( n_new )
+    IMPLICIT NONE
+    INTEGER, INTENT(in) :: n_new ! New size of arrays
+
+    INTEGER :: n_old
+
+    n_old = SIZE ( list )
+
+    IF ( n_new <= n_old ) THEN
+       WRITE ( unit=error_unit, fmt='(a,2i15)') 'Unexpected array size ', n_old, n_new
+       STOP 'Error in extend_list'
+    END IF
+
+    ! Reallocate list array, preserving existing data
+    ALLOCATE ( tmp_list(n_new) )
+    tmp_list(1:n_old)  = list
+    tmp_list(n_old+1:) = 0
+    CALL MOVE_ALLOC ( tmp_list, list )
+
+    ! Reallocate c array, preserving existing data
+    ALLOCATE ( tmp_c(3,n_new) )
+    tmp_c(:,1:n_old)  = c
+    tmp_c(:,n_old+1:) = 0
+    CALL MOVE_ALLOC ( tmp_c, c )
+
+  END SUBROUTINE extend_list
+
   SUBROUTINE finalize_list ! Routine to deallocate list arrays
     IMPLICIT NONE
 

mc_lj_ll_module.f90

@@ -34,7 +34,7 @@ MODULE mc_module
   PRIVATE
 
   ! Public routines
-  PUBLIC :: introduction, conclusion, allocate_arrays, deallocate_arrays
+  PUBLIC :: introduction, conclusion, allocate_arrays, extend_arrays, deallocate_arrays
   PUBLIC :: box_check, potential_1, potential, force_sq
   PUBLIC :: move, create, destroy
 
@@ -45,6 +45,7 @@ MODULE mc_module
   ! Private data
   REAL,    DIMENSION(:,:), ALLOCATABLE :: f      ! Forces for force_sq calculation (3,n)
   INTEGER, DIMENSION(:),   ALLOCATABLE :: j_list ! List of j-neighbours (n)
+  REAL,    DIMENSION(:,:), ALLOCATABLE :: tmp_r  ! Temporary array for dynamic reallocation of r
 
   INTEGER, PARAMETER :: lt = -1, gt = 1 ! Options for j-range
 
@@ -120,6 +121,38 @@ SUBROUTINE allocate_arrays ( box, r_cut )
 
   END SUBROUTINE allocate_arrays
 
+  SUBROUTINE extend_arrays ( n_new )
+    USE link_list_module, ONLY : extend_list
+    IMPLICIT NONE
+    INTEGER, INTENT(in) :: n_new ! New size of arrays
+
+    INTEGER :: n_old
+
+    n_old = SIZE ( r, dim=2 )
+
+    IF ( n_new <= n_old ) THEN
+       WRITE ( unit=error_unit, fmt='(a,2i15)') 'Unexpected array size ', n_old, n_new
+       STOP 'Error in extend_arrays'
+    END IF
+
+    ! Reallocate r array, preserving existing data
+    ALLOCATE ( tmp_r(3,n_new) )
+    tmp_r(:,1:n_old)  = r
+    tmp_r(:,n_old+1:) = 0.0
+    CALL MOVE_ALLOC ( tmp_r, r )
+
+    ! Reallocate f array, no need to preserve data
+    DEALLOCATE ( f )
+    ALLOCATE ( f(3,n_new) )
+
+    ! Reallocate j_list array, no need to preserve data
+    DEALLOCATE ( j_list )
+    ALLOCATE ( j_list(n_new) )
+
+    CALL extend_list ( n_new )
+
+  END SUBROUTINE extend_arrays
+
   SUBROUTINE deallocate_arrays
     USE link_list_module, ONLY : finalize_list
     IMPLICIT NONE

mc_lj_module.f90

@@ -31,7 +31,7 @@ MODULE mc_module
   PRIVATE
 
   ! Public routines
-  PUBLIC :: introduction, conclusion, allocate_arrays, deallocate_arrays
+  PUBLIC :: introduction, conclusion, allocate_arrays, extend_arrays, deallocate_arrays
   PUBLIC :: box_check, potential_1, potential, force_sq
   PUBLIC :: move, create, destroy
 
@@ -40,7 +40,8 @@ MODULE mc_module
   REAL,    DIMENSION(:,:), ALLOCATABLE, PUBLIC :: r ! Positions (3,n)
 
   ! Private data
-  REAL, DIMENSION(:,:), ALLOCATABLE :: f ! Forces for force_sq calculation (3,n)
+  REAL, DIMENSION(:,:), ALLOCATABLE :: f     ! Forces for force_sq calculation (3,n)
+  REAL, DIMENSION(:,:), ALLOCATABLE :: tmp_r ! Temporary array for dynamic reallocation of r
 
   INTEGER, PARAMETER :: lt = -1, gt = 1 ! Options for j-range
 
@@ -113,6 +114,31 @@ SUBROUTINE allocate_arrays ( box, r_cut )
 
   END SUBROUTINE allocate_arrays
 
+  SUBROUTINE extend_arrays ( n_new )
+    IMPLICIT NONE
+    INTEGER, INTENT(in) :: n_new ! New size of arrays
+
+    INTEGER :: n_old
+
+    n_old = SIZE ( r, dim=2 )
+
+    IF ( n_new <= n_old ) THEN
+       WRITE ( unit=error_unit, fmt='(a,2i15)') 'Unexpected array size ', n_old, n_new
+       STOP 'Error in extend_arrays'
+    END IF
+
+    ! Reallocate r array, preserving existing data
+    ALLOCATE ( tmp_r(3,n_new) )
+    tmp_r(:,1:n_old)  = r
+    tmp_r(:,n_old+1:) = 0.0
+    CALL MOVE_ALLOC ( tmp_r, r )
+
+    ! Reallocate f array, no need to preserve data
+    DEALLOCATE ( f )
+    ALLOCATE ( f(3,n_new) )
+
+  END SUBROUTINE extend_arrays
+
   SUBROUTINE deallocate_arrays
     IMPLICIT NONE
 

mc_zvt_lj.f90

@@ -48,7 +48,7 @@ PROGRAM mc_zvt_lj
   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 : metropolis, random_integer, random_translate_vector
-  USE               mc_module,        ONLY : introduction, conclusion, allocate_arrays, deallocate_arrays, &
+  USE               mc_module,        ONLY : introduction, conclusion, allocate_arrays, extend_arrays, deallocate_arrays, &
        &                                     potential_1, potential, move, create, destroy, n, r, potential_type
 
   IMPLICIT NONE
@@ -116,12 +116,11 @@ PROGRAM mc_zvt_lj
   WRITE ( unit=output_unit, fmt='(a,t40,i15)'   ) 'Number of particles',   n
   WRITE ( unit=output_unit, fmt='(a,t40,f15.6)' ) 'Simulation box length', box
   WRITE ( unit=output_unit, fmt='(a,t40,f15.6)' ) 'Density',               REAL(n) / box**3
-  n = n * 2 ! Increase n for array allocation
   CALL allocate_arrays ( box, r_cut ) ! Allocate r
-  n = n / 2 ! Restore value of n
   CALL read_cnf_atoms ( cnf_prefix//inp_tag, n, box, r ) ! Second call is to get r
   r(:,:) = r(:,:) / box              ! Convert positions to box units
   r(:,:) = r(:,:) - ANINT ( r(:,:) ) ! Periodic boundaries
+  CALL extend_arrays ( 2*n )         ! Extend arrays to allow for particle creation
 
   ! Initial energy and overlap check
   total = potential ( box, r_cut ) 
@@ -191,8 +190,7 @@ PROGRAM mc_zvt_lj
               c_try = c_try + 1
 
               IF ( n+1 > SIZE(r,dim=2) ) THEN
-                 WRITE ( unit=error_unit, fmt='(a,2i5)') 'n has grown too large', n+1, SIZE(r,dim=2)
-                 STOP 'Error in mc_zvt_lj'
+                 CALL extend_arrays ( 2*n ) ! Extend arrays to allow for particle creation
               END IF
 
               CALL RANDOM_NUMBER ( ri ) ! Three uniform random numbers in range (0,1)