diff --git a/GUIDE.md b/GUIDE.md index 3c7042c..805565a 100644 --- a/GUIDE.md +++ b/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. diff --git a/link_list_module.f90 b/link_list_module.f90 index ada81bc..ba07e45 100644 --- a/link_list_module.f90 +++ b/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 diff --git a/mc_lj_ll_module.f90 b/mc_lj_ll_module.f90 index 0556910..6f95052 100644 --- a/mc_lj_ll_module.f90 +++ b/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 diff --git a/mc_lj_module.f90 b/mc_lj_module.f90 index 2dede10..1b447dc 100644 --- a/mc_lj_module.f90 +++ b/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 diff --git a/mc_zvt_lj.f90 b/mc_zvt_lj.f90 index 8617887..544994e 100644 --- a/mc_zvt_lj.f90 +++ b/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)