diff --git a/GUIDE.md b/GUIDE.md index afba81f..3c7042c 100644 --- a/GUIDE.md +++ b/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_. diff --git a/SConstruct b/SConstruct index cec543d..d5306a4 100644 --- a/SConstruct +++ b/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) diff --git a/link_list_module.f90 b/link_list_module.f90 index 49564b1..ada81bc 100644 --- a/link_list_module.f90 +++ b/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 diff --git a/mc_lj_ll_module.f90 b/mc_lj_ll_module.f90 index ed373c7..60f627d 100644 --- a/mc_lj_ll_module.f90 +++ b/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,7 +299,7 @@ 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 @@ -293,6 +307,7 @@ FUNCTION force_sq ( box, r_cut ) RESULT ( fsq ) ! 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 diff --git a/mc_lj_module.f90 b/mc_lj_module.f90 index 0b65d13..2dede10 100644 --- a/mc_lj_module.f90 +++ b/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 diff --git a/mc_npt_lj.f90 b/mc_npt_lj.f90 index c8475ff..9ca332b 100644 --- a/mc_npt_lj.f90 +++ b/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 diff --git a/md_lj_le_module.f90 b/md_lj_le_module.f90 index 0caab32..b5acfa3 100644 --- a/md_lj_le_module.f90 +++ b/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 diff --git a/md_lj_ll_module.f90 b/md_lj_ll_module.f90 index 9c85256..f6527b4 100644 --- a/md_lj_ll_module.f90 +++ b/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 diff --git a/md_lj_llle_module.f90 b/md_lj_llle_module.f90 index 4149e7a..f8eddf6 100644 --- a/md_lj_llle_module.f90 +++ b/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 diff --git a/md_lj_module.f90 b/md_lj_module.f90 index e9cd6dd..712f331 100644 --- a/md_lj_module.f90 +++ b/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 diff --git a/python_examples/md_lj_ll_module.py b/python_examples/md_lj_ll_module.py index 1a30de5..4fac5ee 100644 --- a/python_examples/md_lj_ll_module.py +++ b/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= 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= 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