Minor improvements to link cells when box size varies. This closes #14.

GUIDE.md

@@ -208,6 +208,7 @@ Source                 | ρ     | _T_       | _E_ (cs)   | _P_ (cs) | _C_ (cs
 `md_nve_lj_vl`         | 0.75      | 1.0023(3) | -2.9278    | 0.992(2) |  2.24(1)  | -3.7274    | 0.391(2) |
 `md_nve_lj_ll`‡ | 0.75      | 1.0010(1) | -2.9272    | 0.992(1) |  2.28(1)  | -3.7268    | 0.391(1) |
 `md_nvt_lj_ll`‡ | 0.75      | 1.00      | -2.927(2)  | 0.994(3) |  2.3(1)   | -3.727(2)  | 0.392(3) | 2.3(1)
+`md_npt_lj_ll`§‡ | 0.7501(4) | 1.00 | -2.931(4) | 0.991(1) |           | -3.731(4)  | 0.389(1) |
 `smc_nvt_lj`♯(a) | 0.75      | 1.00      | -2.9300(5) | 0.971(2) |  2.263(5) | -3.7296(5) | 0.369(2) | 2.270(5)
 `smc_nvt_lj`♯(b) | 0.75      | 1.00      | -2.928(2)  | 0.99(1)  |  2.26(2)  | -3.728(2)  | 0.39(1)  | 2.27(2)
 `smc_nvt_lj`♯(c) | 0.75      | 1.00      | -2.930(3)  | 0.98(2)  |  2.26(3)  | -3.729(3)  | 0.38(2)  | 2.27(3)
@@ -216,7 +217,7 @@ Source                 | ρ     | _T_       | _E_ (cs)   | _P_ (cs) | _C_ (cs
 ‡ Indicates a larger system size, _N_=864, needed to make the link-list method viable. Note that
 the speedup is not enormous for this system size, corresponding to 4×4×4 cells.
 
-§ The constant-pressure simulation was run at _P_=0.99, the program default.
+§ The constant-pressure simulations were run at _P_=0.99, the program default.
 
 ♯ The `smc_nvt_lj` program was tested (a) in default, single-particle-move, mode, with δt=0.1;
 (b) in multi-particle mode, moving 100% of particles, with δt=0.02;
@@ -295,13 +296,12 @@ 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.
 
-Several of these programs could be improved to use array reallocation (available in Fortran)
-to make them more resilient against changes in box size or number of particles.
-For simplicity we have not included these features.
-
 The measured pressures _P_ (c) are systematically a little low;
 this is particularly noticeable for the constant-pressure programs,
 where they might be expected to agree with the user-defined value of _P_.

SConstruct

@@ -92,6 +92,7 @@ variants['build_md_chain_mts_lj']      = (['md_chain_mts_lj.f90','md_chain_lj_mo
 variants['build_md_chain_nve_lj']      = (['md_chain_nve_lj.f90','md_chain_lj_module.f90']+utilities,env_lapack)
 variants['build_md_lj_mts']            = (['md_lj_mts.f90','md_lj_mts_module.f90','lrc_lj_module.f90']+utnomaths,env_normal)
 variants['build_md_npt_lj']            = (['md_npt_lj.f90','md_lj_module.f90','lrc_lj_module.f90']+utilities,env_normal)
+variants['build_md_npt_lj_ll']         = (['md_npt_lj.f90','md_lj_ll_module.f90','lrc_lj_module.f90','link_list_module.f90']+utilities,env_normal)
 variants['build_md_nve_hs']            = (['md_nve_hs.f90','md_nve_hs_module.f90']+utnomaths,env_normal)
 variants['build_md_nve_lj']            = (['md_nve_lj.f90','md_lj_module.f90','lrc_lj_module.f90']+utnomaths,env_normal)
 variants['build_md_nve_lj_vl']         = (['md_nve_lj.f90','md_lj_vl_module.f90','lrc_lj_module.f90','verlet_list_module.f90']+utnomaths,env_normal)

link_list_module.f90

@@ -45,11 +45,12 @@ MODULE link_list_module
   SUBROUTINE initialize_list ( n, r_cut_box ) ! Routine to allocate list arrays
     IMPLICIT NONE
     INTEGER, INTENT(in) :: n         ! Number of particles
-    REAL,    INTENT(in) :: r_cut_box ! rcut/box, assume never changes
+    REAL,    INTENT(in) :: r_cut_box ! rcut/box
 
-    WRITE ( unit=output_unit, fmt='(a,t40,f15.6)') 'Link cells based on r_cut/box =', r_cut_box
+    WRITE ( unit=output_unit, fmt='(a,t40,f15.6)') 'Link cells based on r_cut/box', r_cut_box
 
     sc = FLOOR ( 1.0 / r_cut_box ) ! Number of cells in each dimension
+    WRITE ( unit=output_unit, fmt='(a,t40,i15)') 'Number of cells in box', sc
     IF ( sc < 3 ) THEN
        WRITE ( unit=error_unit, fmt='(a,i15)') 'System is too small to use link cells', sc
        STOP 'Error in initialize_list'
@@ -68,13 +69,25 @@ SUBROUTINE finalize_list ! Routine to deallocate list arrays
 
   END SUBROUTINE finalize_list
 
-  SUBROUTINE make_list ( n, r ) ! Routine to make list
+  SUBROUTINE make_list ( n, r_cut_box, r ) ! Routine to make list
     IMPLICIT NONE
-    INTEGER,                 INTENT(in) :: n ! Number of atoms
-    REAL,    DIMENSION(3,n), INTENT(in) :: r ! Atom coordinates
+    INTEGER,                 INTENT(in) :: n         ! Number of atoms
+    REAL,                    INTENT(in) :: r_cut_box ! rcut/box
+    REAL,    DIMENSION(3,n), INTENT(in) :: r         ! Atom coordinates
 
     INTEGER :: i
-
+
+    ! Allow for possibility of box size change
+    sc = FLOOR ( 1.0 / r_cut_box )
+    IF ( sc < 3 ) THEN
+       WRITE ( unit=error_unit, fmt='(a,i15)') 'System is too small to use link cells', sc
+       STOP 'Error in make_list'
+    END IF
+    IF ( sc > SIZE ( head, dim=1 ) ) THEN
+       DEALLOCATE ( head )
+       ALLOCATE ( head(0:sc-1,0:sc-1,0:sc-1) )
+    END IF
+
     head(:,:,:) = 0
 
     DO i = 1, n ! Loop over all atoms

mc_lj_ll_module.f90

@@ -39,7 +39,7 @@ MODULE mc_module
 
   ! Public routines
   PUBLIC :: introduction, conclusion, allocate_arrays, deallocate_arrays
-  PUBLIC :: potential_1, potential, force_sq
+  PUBLIC :: box_check, potential_1, potential, force_sq
   PUBLIC :: move, create, destroy
 
   ! Public data
@@ -134,6 +134,19 @@ SUBROUTINE deallocate_arrays
 
   END SUBROUTINE deallocate_arrays
 
+  SUBROUTINE box_check ( box, r_cut )
+    USE link_list_module, ONLY : make_list
+    IMPLICIT NONE
+    REAL, INTENT(in) :: box   ! Simulation box length
+    REAL, INTENT(in) :: r_cut ! Potential cutoff distance
+
+    REAL :: r_cut_box
+
+    r_cut_box = r_cut / box
+    CALL make_list ( n, r_cut_box, r )
+
+  END SUBROUTINE box_check
+
   FUNCTION potential ( box, r_cut ) RESULT ( total )
     USE link_list_module, ONLY : make_list
     IMPLICIT NONE
@@ -148,26 +161,23 @@ FUNCTION potential ( box, r_cut ) RESULT ( total )
     ! If this flag is .true., the values of total%pot etc should not be used
     ! Actual calculation is performed by function potential_1
 
-    ! We assume that the main program calls this function at the start,
-    ! and hence use this opportunity to call the make_list function.
-    ! We also assume that the box length remains constant throughout
-    ! The list is maintained by create_in_list, destroy_in_list, move_in_list
-    ! which are called in other routines below.
+    ! This routine uses linked lists.
+    ! We assume that the box length may vary, so we call the make_list function at the start.
+    ! The list for fixed box length is maintained by create_in_list, destroy_in_list,
+    ! and move_in_list which are called in other routines below.
 
     TYPE(potential_type) :: partial
     INTEGER              :: i
-    LOGICAL, SAVE        :: first_call = .TRUE.
-    
+    REAL                 :: r_cut_box
+
     IF ( n > SIZE(r,dim=2) ) THEN ! should never happen
        WRITE ( unit=error_unit, fmt='(a,2i15)' ) 'Array bounds error for r', n, SIZE(r,dim=2)
        STOP 'Impossible error in potential'
     END IF
 
-    IF ( first_call ) THEN
-       r(:,:) = r(:,:) - ANINT ( r(:,:) ) ! Periodic boundaries
-       CALL make_list ( n, r )
-       first_call = .FALSE.
-    END IF
+    r(:,:) = r(:,:) - ANINT ( r(:,:) ) ! Periodic boundaries
+    r_cut_box = r_cut / box
+    CALL make_list ( n, r_cut_box, r )
 
     total = potential_type ( pot=0.0, vir=0.0, lap=0.0, ovr=.FALSE. ) ! Initialize
 
@@ -210,6 +220,10 @@ FUNCTION potential_1 ( ri, i, box, r_cut, j_range ) RESULT ( partial )
     ! It is assumed that r has been divided by box
     ! Results are in LJ units where sigma = 1, epsilon = 1
 
+    ! This routine uses linked lists.
+    ! We assume that the box length may vary, but that the make_list function
+    ! has been called since the last change in box length.
+
     INTEGER               :: j, jj
     LOGICAL               :: half
     REAL                  :: r_cut_box, r_cut_box_sq, box_sq
@@ -285,14 +299,15 @@ FUNCTION potential_1 ( ri, i, box, r_cut, j_range ) RESULT ( partial )
   END FUNCTION potential_1
 
   FUNCTION force_sq ( box, r_cut ) RESULT ( fsq )
-    USE link_list_module, ONLY : c, neighbours
+    USE link_list_module, ONLY : make_list, c, neighbours
     IMPLICIT NONE
     REAL             :: fsq   ! Returns total squared force
     REAL, INTENT(in) :: box   ! Simulation box length
     REAL, INTENT(in) :: r_cut ! Potential cutoff distance
 
     ! Calculates total squared force (using array f)
     ! Uses link lists
+    ! We assume that the box length may vary, so we call the make_list function at the start.
 
     INTEGER               :: i, j, jj
     REAL                  :: r_cut_box, r_cut_box_sq, box_sq, rij_sq
@@ -302,6 +317,7 @@ FUNCTION force_sq ( box, r_cut ) RESULT ( fsq )
     r_cut_box    = r_cut / box
     r_cut_box_sq = r_cut_box ** 2
     box_sq       = box ** 2
+    CALL make_list ( n, r_cut_box, r )
 
     ! Initialize
     f = 0.0

mc_lj_module.f90

@@ -32,7 +32,7 @@ MODULE mc_module
 
   ! Public routines
   PUBLIC :: introduction, conclusion, allocate_arrays, deallocate_arrays
-  PUBLIC :: potential_1, potential, force_sq
+  PUBLIC :: box_check, potential_1, potential, force_sq
   PUBLIC :: move, create, destroy
 
   ! Public data
@@ -120,6 +120,21 @@ SUBROUTINE deallocate_arrays
 
   END SUBROUTINE deallocate_arrays
 
+  SUBROUTINE box_check ( box, r_cut )
+    IMPLICIT NONE
+    REAL, INTENT(in) :: box   ! Simulation box length
+    REAL, INTENT(in) :: r_cut ! Potential cutoff distance
+
+    REAL :: r_cut_box
+
+    r_cut_box = r_cut / box
+    IF ( r_cut_box > 0.5 ) THEN
+       WRITE ( unit=error_unit, fmt='(a,f15.6)') 'r_cut/box too large ', r_cut_box
+       STOP 'Error in box_check'
+    END IF
+
+  END SUBROUTINE box_check
+
   FUNCTION potential ( box, r_cut ) RESULT ( total )
     IMPLICIT NONE
     TYPE(potential_type) :: total ! Returns a composite of pot, vir etc
@@ -133,14 +148,23 @@ FUNCTION potential ( box, r_cut ) RESULT ( total )
     ! If this flag is .true., the values of total%pot etc should not be used
     ! Actual calculation is performed by function potential_1
 
+    ! We assume that the box length may vary, so we check r_cut/box at the start.
+
     TYPE(potential_type) :: partial ! Atomic contribution to total
     INTEGER              :: i
+    REAL                 :: r_cut_box
 
     IF ( n > SIZE(r,dim=2) ) THEN ! should never happen
        WRITE ( unit=error_unit, fmt='(a,2i15)' ) 'Array bounds error for r', n, SIZE(r,dim=2)
        STOP 'Error in potential'
     END IF
 
+    r_cut_box = r_cut / box
+    IF ( r_cut_box > 0.5 ) THEN
+       WRITE ( unit=error_unit, fmt='(a,f15.6)') 'r_cut/box too large ', r_cut_box
+       STOP 'Error in potential'
+    END IF
+
     total = potential_type ( pot=0.0, vir=0.0, lap=0.0, ovr=.FALSE. ) ! Initialize
 
     DO i = 1, n - 1
@@ -180,6 +204,9 @@ FUNCTION potential_1 ( ri, i, box, r_cut, j_range ) RESULT ( partial )
     ! It is assumed that r has been divided by box
     ! Results are in LJ units where sigma = 1, epsilon = 1
 
+    ! We assume that the box length may vary, but that r_cut/box
+    ! has been checked since the last change in box length.
+
     INTEGER              :: j, j1, j2
     REAL                 :: r_cut_box, r_cut_box_sq, box_sq
     REAL                 :: sr2, sr6, sr12, rij_sq
@@ -270,6 +297,10 @@ FUNCTION force_sq ( box, r_cut ) RESULT ( fsq )
     r_cut_box    = r_cut / box
     r_cut_box_sq = r_cut_box ** 2
     box_sq       = box ** 2
+    IF ( r_cut_box > 0.5 ) THEN
+       WRITE ( unit=error_unit, fmt='(a,f15.6)') 'r_cut/box too large ', r_cut_box
+       STOP 'Error in force_sq'
+    END IF
 
     f = 0.0 ! Initialize
 

mc_npt_lj.f90

@@ -54,7 +54,7 @@ PROGRAM mc_npt_lj
   USE               averages_module,  ONLY : run_begin, run_end, blk_begin, blk_end, blk_add
   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
+       &                                     box_check, potential_1, potential, move, n, r, potential_type
 
   IMPLICIT NONE
 
@@ -143,6 +143,7 @@ PROGRAM mc_npt_lj
 
      DO stp = 1, nstep ! Begin loop over steps
 
+        CALL box_check ( box, r_cut ) ! Prepare for single-particle moves
         moves = 0
 
         DO i = 1, n ! Begin loop over atoms

md_lj_le_module.f90

@@ -137,6 +137,10 @@ SUBROUTINE force ( box, strain, total )
     r_cut_box    = r_cut / box
     r_cut_box_sq = r_cut_box ** 2
     box_sq       = box ** 2
+    IF ( r_cut_box > 0.5 ) THEN
+       WRITE ( unit=error_unit, fmt='(a,f15.6)' ) 'r_cut/box too large ', r_cut_box
+       STOP 'Error in force'
+    END IF
 
     ! Initialize
     f     = 0.0

md_lj_ll_module.f90

@@ -158,7 +158,7 @@ SUBROUTINE force ( box, r_cut, total )
     pot_cut = sr12 - sr6 ! Without numerical factor 4
 
     r(:,:) = r(:,:) - ANINT ( r(:,:) ) ! Periodic boundary conditions in box=1 units
-    CALL make_list ( n, r )
+    CALL make_list ( n, r_cut_box, r )
 
     ! Initialize
     f     = 0.0
@@ -278,7 +278,7 @@ FUNCTION hessian ( box, r_cut ) RESULT ( hes )
     hes = 0.0
 
     r(:,:) = r(:,:) - ANINT ( r(:,:) ) ! Periodic boundary conditions in box=1 units
-    CALL make_list ( n, r )
+    CALL make_list ( n, r_cut_box, r )
 
     ! Triple loop over cells
     DO ci1 = 0, sc-1

md_lj_llle_module.f90

@@ -163,7 +163,7 @@ SUBROUTINE force ( box, strain, total )
     ! Lees-Edwards periodic boundaries
     r(1,:) = r(1,:) - ANINT ( r(2,:) ) * strain ! Extra correction in box=1 units
     r(:,:) = r(:,:) - ANINT ( r(:,:) )          ! Standard correction in box=1 units
-    CALL make_list ( n, r )
+    CALL make_list ( n, r_cut_box, r )
 
     shift = FLOOR ( strain * REAL ( sc ) ) ! Strain measured in cell lengths
 

md_lj_module.f90

@@ -135,6 +135,10 @@ SUBROUTINE force ( box, r_cut, total )
     r_cut_box    = r_cut / box
     r_cut_box_sq = r_cut_box ** 2
     box_sq       = box ** 2
+    IF ( r_cut_box > 0.5 ) THEN
+       WRITE ( unit=error_unit, fmt='(a,f15.6)' ) 'r_cut/box too large ', r_cut_box
+       STOP 'Error in force'
+    END IF
 
     ! Calculate potential at cutoff
     sr2     = 1.0 / r_cut**2 ! in sigma=1 units

python_examples/md_lj_ll_module.py

@@ -104,6 +104,7 @@ def force ( box, r_cut, r ):
 
     # Calculate cell index triplets
     sc = math.floor(box/r_cut)                          # Number of cells along box edge
+    assert sc >= 3, 'System is too small for cells'     # Guard against box being too small
     c  = np.floor((r+0.5)*sc).astype(np.int_)           # N*3 array of cell indices for all atoms
     assert np.all(c>=0) and np.all(c<sc), 'Index error' # Simplistic "guard" against roundoff
 
@@ -255,6 +256,7 @@ def hessian ( box, r_cut, r, f ):
 
     # Calculate cell index triplets
     sc = math.floor(box/r_cut)                          # Number of cells along box edge
+    assert sc >= 3, 'System is too small for cells'     # Guard against box being too small
     c  = np.floor((r+0.5)*sc).astype(np.int_)           # N*3 array of cell indices for all atoms
     assert np.all(c>=0) and np.all(c<sc), 'Index error' # Simplistic "guard" against roundoff
 

python_examples/md_lj_llle_module.py

@@ -105,6 +105,7 @@ def force ( box, strain, r ):
 
     # Calculate cell index triplets
     sc = math.floor(box/r_cut)                          # Number of cells along box edge
+    assert sc >= 3, 'System is too small for cells'     # Guard against box being too small
     c  = np.floor((r+0.5)*sc).astype(np.int_)           # N*3 array of cell indices for all atoms
     assert np.all(c>=0) and np.all(c<sc), 'Index error' # Simplistic "guard" against roundoff