corrected and tested Python smc_nvt_lj example

GUIDE.md

@@ -150,7 +150,6 @@ Source                 | ρ     | _T_       | _E_ (cs)   | _P_ (cs) | _C_ (cs
 ------                 | -----     | -----     | --------   | -------- | --------- | -------    | -------  | --------
 Thol et al (2015) (cs) | 0.75      | 1.00      | -2.9286    | 0.9897   |  2.2787   |            |          |          
 Thol et al (2016) (f)  | 0.75      | 1.00      |            |          |           | -3.7212    | 0.3996   | 2.2630  
-`bd_nvt_lj`            | 0.75      | 1.00      | -2.934(4)  | 0.974(7) |  2.26(8)  | -3.733(4)  | 0.373(7) | 2.27(8)
 `md_nvt_lj`            | 0.75      | 1.00      | -2.940(4)  | 0.965(6) |  2.27(12) | -3.740(4)  | 0.363(6) | 2.27(12)
 `md_npt_lj`§      | 0.7514(6) | 1.00      | -2.947(6)  | 0.995(1) |           | -3.748(7)  | 0.391(1) |
 `md_nve_lj`            | 0.75      | 1.0022(3) | -2.9289    | 0.987(2) |  2.24(1)  | -3.7284    | 0.386(2) |          
@@ -161,6 +160,7 @@ Thol et al (2016) (f)  | 0.75      | 1.00      |            |          |
 `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)
+`bd_nvt_lj`            | 0.75      | 1.00      | -2.934(4)  | 0.974(7) |  2.26(8)  | -3.733(4)  | 0.373(7) | 2.27(8)
 
 ‡ 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.

python_examples/GUIDE.md

@@ -175,7 +175,7 @@ at the start of the main module; the user may override this by editing this stat
 but the program will then run (yet another) order of magnitude more slowly,
 so it may be necessary to reduce the run length still further.
 
-### Lennard-Jones MD programs
+### Lennard-Jones MD, BD and SMC programs
 
 The first test of the MD codes is that energy, or the appropriate energy-like variable, is conserved.
 The following table uses runs of 10 blocks,
@@ -197,7 +197,8 @@ except for some small deviations at the smallest timestep.
 
 Now we compare EOS data with typical test runs of our programs.
 The results in the following table use the same run lengths
-as for the Fortran examples, and default parameters otherwise.
+as for the Fortran examples (i.e. longer than the default value in the program),
+and default parameters otherwise.
 Numbers in parentheses (here and in the following tables)
 indicate errors in the last quoted digit, estimated from block averages.
 Results without error estimates are fixed (such as the temperature or density) or conserved.
@@ -209,8 +210,16 @@ Thol et al (2016) (f)  | 0.75      | 1.00      |            |          |
 `md_nve_lj.py`         | 0.75      | 1.0027(4) | -2.9278    | 0.988(3) |  2.24(1)  | -3.7274    | 0.387(3)  |          
 `md_nvt_lj.py`         | 0.75      | 1.00      | -2.937(3)  | 0.975(4) |  2.1(1)   | -3.737(3)  | 0.374(4)  | 2.1(1)
 `md_npt_lj.py`         | 0.7509(5) | 1.00      | -2.942(5)  | 0.994(1) |           | -3.743(6)  | 0.3908(6) |
+`smc_nvt_lj`♯(a) | 0.75      | 1.00      | -2.929(1)  | 0.979(3) |  2.256(4) | -3.729(1)  | 0.378(3)  | 2.264(4)
+`smc_nvt_lj`♯(b) | 0.75      | 1.00      | -2.932(2)  | 0.95(1)  |  2.28(3)  | -3.732(2)  | 0.35(1)   | 2.29(3)
+`smc_nvt_lj`♯(c) | 0.75      | 1.00      | -2.934(3)  | 0.94(1)  |  2.24(2)  | -3.733(3)  | 0.34(1)   | 2.24(2)
 `bd_nvt_lj.py`         | 0.75      | 1.00      | -2.931(3)  | 0.976(5) |  2.2(1)   | -3.731(3)  | 0.374(5)  | 2.2(1)
 
+♯ 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;
+and (c) in multi-particle mode, moving 30% of particles, with δt=0.03.
+These values give acceptance rates in the 45% – 55% range.
+
 ### Lennard-Jones MC programs
 
 The results in the following table use the same run lengths

python_examples/smc_lj_module.py

@@ -26,7 +26,7 @@
 
 """Energy, force, and move routines for SMC, LJ potential."""
 
-fast = False # Change this to replace NumPy force evaluation with slower Python
+fast = True # Change this to replace NumPy force evaluation with slower Python
 
 class PotentialType:
     """A composite variable for interactions."""
@@ -140,6 +140,7 @@ def force_1 ( ri, box, r_cut, r ):
         rij_sq = np.sum(rij**2,axis=1)  # Squared separations
         in_range = rij_sq < r_cut_box_sq                    # Set mask for within cutoff
         rij_sq   = rij_sq * box_sq                          # Now in sigma=1 units
+        rij      = rij * box                                # Now in sigma=1 units
         sr2      = np.where ( in_range, 1.0 / rij_sq, 0.0 ) # (sigma/rij)**2, only if in range
         ovr      = sr2 > sr2_ovr                            # Set flags for any overlaps
         if np.any(ovr):
@@ -165,6 +166,7 @@ def force_1 ( ri, box, r_cut, r ):
             rij_sq = np.sum(rij**2)  # Squared separation
             if rij_sq < r_cut_box_sq:    # Check within cutoff
                 rij_sq = rij_sq * box_sq # Now in sigma=1 units
+                rij    = rij * box       # Now in sigma=1 units
                 sr2    = 1.0 / rij_sq    # (sigma/rij)**2
                 ovr    = sr2 > sr2_ovr   # Overlap if too close
                 if ovr:

smc_nvt_lj.f90

@@ -58,7 +58,7 @@ PROGRAM smc_nvt_lj
   ! Most important variables
   REAL    :: box         ! Box length
   REAL    :: temperature ! Temperature (specified)
-  INTEGER :: move_mode   ! Selects single- or multi-atom moves
+  LOGICAL :: single_atom ! Selects single- or multi-atom moves
   REAL    :: fraction    ! Fraction of atoms to move in multi-atom move
   REAL    :: dt          ! Time step (effectively determines typical displacement)
   REAL    :: r_cut       ! Potential cutoff distance
@@ -74,11 +74,8 @@ PROGRAM smc_nvt_lj
   CHARACTER(len=3), PARAMETER :: inp_tag    = 'inp'
   CHARACTER(len=3), PARAMETER :: out_tag    = 'out'
   CHARACTER(len=3)            :: sav_tag    = 'sav' ! May be overwritten with block number
-  CHARACTER(len=30)           :: mode
 
-  INTEGER, PARAMETER :: single_atom = 1, multi_atom = 2
-
-  NAMELIST /nml/ nblock, nstep, r_cut, dt, mode, temperature, fraction
+  NAMELIST /nml/ nblock, nstep, r_cut, dt, single_atom, temperature, fraction
 
   WRITE ( unit=output_unit, fmt='(a)' ) 'smc_nvt_lj'
   WRITE ( unit=output_unit, fmt='(a)' ) 'Smart Monte Carlo, constant-NVT ensemble'
@@ -90,10 +87,10 @@ PROGRAM smc_nvt_lj
   nblock      = 10
   nstep       = 10000
   r_cut       = 2.5
-  temperature = 1.0           ! Default temperature T
-  dt          = 0.1           ! Together with v_rms=sqrt(T) determines typical displacement
-  mode        = 'single-atom' ! Other option is 'multi-atom', probably requiring smaller dt
-  fraction    = 1.0           ! Only applicable in multi-atom mode
+  temperature = 1.0    ! Default temperature T
+  dt          = 0.1    ! Together with v_rms=sqrt(T) determines typical displacement
+  single_atom = .TRUE. ! .false. selects multi-atom mode, probably requiring smaller dt
+  fraction    = 1.0    ! Only applicable in multi-atom mode
 
   ! Read run parameters from namelist
   ! Comment out, or replace, this section if you don't like namelists
@@ -111,20 +108,15 @@ PROGRAM smc_nvt_lj
   WRITE ( unit=output_unit, fmt='(a,t40,f15.6)' ) 'Potential cutoff distance', r_cut
   WRITE ( unit=output_unit, fmt='(a,t40,f15.6)' ) 'Time step',                 dt
   WRITE ( unit=output_unit, fmt='(a,t40,f15.6)' ) 'Temperature',               temperature
-  IF ( INDEX( lowercase(mode), 'sing' ) /= 0 ) THEN
-     move_mode = single_atom
+  IF ( single_atom ) THEN
      WRITE ( unit=output_unit, fmt='(a,t40,a15)' ) 'Move mode is ', 'single-atom'
-  ELSE IF ( INDEX( lowercase(mode), 'mult' ) /= 0 ) THEN
-     move_mode = multi_atom
+  ELSE
      WRITE ( unit=output_unit, fmt='(a,t40,a15)'   ) 'Move mode is ', 'multi-atom'
      WRITE ( unit=output_unit, fmt='(a,t40,f15.6)' ) 'Fraction of atoms moving', fraction
      IF ( fraction < 0.0 .OR. fraction > 1.0 ) THEN
         WRITE ( unit=error_unit, fmt='(a)') 'Error: fraction out of range'
         STOP 'Error in smc_nvt_lj'
      END IF
-  ELSE
-     WRITE ( unit=error_unit, fmt='(a,a)') 'Error: unrecognized move mode ', mode
-     STOP 'Error in smc_nvt_lj'
   END IF
   v_rms = SQRT ( temperature ) ! RMS value for velocity selection
   WRITE ( unit=output_unit, fmt='(a,t40,f15.6)' ) 'Typical dr', v_rms*dt
@@ -156,9 +148,7 @@ PROGRAM smc_nvt_lj
 
      DO stp = 1, nstep ! Begin loop over steps
 
-        SELECT CASE ( move_mode )
-
-        CASE ( single_atom )
+        IF ( single_atom ) THEN ! Single-atom moves
 
            n_move = 0
            DO i = 1, n ! Loop over atoms
@@ -200,7 +190,7 @@ PROGRAM smc_nvt_lj
 
            m_ratio = REAL(n_move) / REAL(n)
 
-        CASE ( multi_atom )
+        ELSE ! Multi-atom moves
 
            CALL RANDOM_NUMBER ( zeta )                         ! Select N uniform random numbers
            move = SPREAD ( zeta < fraction, dim=1, ncopies=3 ) ! Construct mask for moving atoms
@@ -240,7 +230,7 @@ PROGRAM smc_nvt_lj
 
            END IF ! End test for overlap
 
-        END SELECT
+        END IF
 
         ! Calculate and accumulate variables for this step
         CALL blk_add ( calc_variables() )