introducing scons

README.md

@@ -7,6 +7,9 @@ The intention is to clarify points made in the text,
 rather than to provide a piece of code
 suitable for direct use in a research application.
 
+## Warning
+These programs are currently in development. They are not yet ready.
+
 ## Disclaimer
 We ascribe no commercial value to the programs themselves.
 Although a few complete programs are provided,
@@ -19,5 +22,5 @@ You should always check out a routine for your particular application.
 
 ## Language
 The programs contain some explanatory comments, and
-are written, in the main, in Fortran 2003.
+are written, in the main, in Fortran 2003/2008.
 Fortran provides a standard way of inter-operating with C.

SConscript

@@ -0,0 +1,2 @@
+Import('env','sources')
+env.Program(sources)

SConstruct

@@ -0,0 +1,25 @@
+# SCons build file for the various programs
+# Run this typing 'scons'
+# 'Clean' the program by 'scons -c' (i.e., remove object, module and executable)
+
+import os, sys
+
+env_normal=Environment(ENV=os.environ)
+env_mpi=Environment(ENV=os.environ,F90='mpif90',F90LINKER='mpif90')
+
+# Assume that gfortran will be used. 
+#Tool('gfortran')(env_normal)
+#Tool('mpif90')(env_mpi)
+#MY_F90FLAGS='-O2 -finline-functions -funswitch-loops -fwhole-file'
+MY_F90FLAGS='-fdefault-real-8 -std=f2008 -Wall'
+MY_LINKFLAGS=''
+
+env_normal.Append(F90FLAGS=MY_F90FLAGS,LINKFLAGS=MY_LINKFLAGS,FORTRANMODDIRPREFIX='-J',FORTRANMODDIR = '${TARGET.dir}',F90PATH='${TARGET.dir}')
+env_mpi.Append(F90FLAGS=MY_F90FLAGS,LINKFLAGS=MY_LINKFLAGS,FORTRANMODDIRPREFIX='-J',FORTRANMODDIR = '${TARGET.dir}',F90PATH='${TARGET.dir}')
+
+variants={}
+variants['build_mc_nvt_lj']    = (['mc_nvt_lj.f90','mc_nvt_lj_module.f90','utility_module.f90'],env_normal)
+
+# Build each variant in appropriate variant directory
+for variant_dir,(sources,env) in variants.iteritems():
+    SConscript('SConscript', variant_dir=variant_dir,exports={'env':env,'sources':sources},duplicate=1)

mc_nvt_lj.f90

@@ -2,7 +2,7 @@
 ! Monte Carlo simulation, constant-NVT ensemble, Lennard-Jones atoms
 PROGRAM mc_nvt_lj
   USE utility_module,   ONLY : read_cnf_atoms, write_cnf_atoms, run_begin, run_end, blk_begin, blk_end, stp_end
-  USE mc_nvt_lj_module, ONLY : energy_i, energy, pot_lrc, vir_lrc
+  USE mc_nvt_lj_module, ONLY : energy_1, energy, pot_lrc, vir_lrc, n, r, ne
   IMPLICIT NONE
 
   ! Takes in a configuration of atoms (positions)
@@ -104,11 +104,11 @@ PROGRAM mc_nvt_lj
            zeta(1:3) = 2.0*zeta(1:3) - 1.0 ! first three now in range (-1,+1)
 
            ri(:) = r(:,i)
-           CALL  energy_i ( ri, i, j_ne_i, sigma, r_cut, pot_old, vir_old, overlap )
+           CALL  energy_1 ( ri, i, ne, sigma, r_cut, pot_old, vir_old, overlap )
            IF ( overlap ) STOP 'Overlap in current configuration'
            ri(:) = ri(:) + zeta(1:3) * dr_max ! trial move to new position
            ri(:) = ri(:) - ANINT ( ri(:) )    ! periodic boundary correction
-           CALL  energy_i ( ri, i, j_ne_i, sigma, r_cut, pot_new, vir_new, overlap )
+           CALL  energy_1 ( ri, i, ne, sigma, r_cut, pot_new, vir_new, overlap )
 
            IF ( .NOT. overlap ) THEN ! consider non-overlapping configuration
 

mc_nvt_lj_module.f90

@@ -2,13 +2,13 @@ MODULE mc_nvt_lj_module
 
   IMPLICIT NONE
   PRIVATE
-  PUBLIC :: n, r, j_lt_i, j_ne_i, j_gt_i
-  PUBLIC :: energy_i, energy, pot_lrc, vir_lrc
+  PUBLIC :: n, r, lt, ne, gt
+  PUBLIC :: energy_1, energy, pot_lrc, vir_lrc
 
   INTEGER                              :: n ! number of atoms
   REAL,    DIMENSION(:,:), ALLOCATABLE :: r ! positions
 
-  INTEGER, PARAMETER :: j_lt_i = -1, j_ne_i = 0, j_gt_i = 1 ! j-range options
+  INTEGER, PARAMETER :: lt = -1, ne = 0, gt = 1 ! j-range options
 
 CONTAINS
 
@@ -31,7 +31,7 @@ SUBROUTINE energy ( sigma, r_cut, pot, vir, overlap )
 
     DO i = 1, n - 1
 
-       CALL energy_i ( r(:,i), i, j_gt_i, sigma, r_cut, pot_i, vir_i, overlap )
+       CALL energy_1 ( r(:,i), i, gt, sigma, r_cut, pot_i, vir_i, overlap )
 
        IF ( overlap ) RETURN
        pot = pot + pot_i
@@ -40,7 +40,7 @@ SUBROUTINE energy ( sigma, r_cut, pot, vir, overlap )
 
   END SUBROUTINE energy
 
-  SUBROUTINE energy_i ( ri, i, j_range, sigma, r_cut, pot, vir, overlap )
+  SUBROUTINE energy_1 ( ri, i, j_range, sigma, r_cut, pot, vir, overlap )
 
     REAL, DIMENSION(3), INTENT(in)  :: ri
     INTEGER,            INTENT(in)  :: i, j_range
@@ -68,13 +68,13 @@ SUBROUTINE energy_i ( ri, i, j_range, sigma, r_cut, pot, vir, overlap )
     overlap = .FALSE.
 
     SELECT CASE ( j_range )
-    CASE ( j_lt_i )
+    CASE ( lt ) ! j < i
        j1 = 1
        j2 = i-1
-    CASE ( j_gt_i )
+    CASE ( gt ) ! j > i
        j1 = i+1
        j2 = n
-    CASE ( j_ne_i )
+    CASE ( ne ) ! j /= i
        j1 = 1
        j2 = n
     END SELECT
@@ -106,7 +106,7 @@ SUBROUTINE energy_i ( ri, i, j_range, sigma, r_cut, pot, vir, overlap )
     pot = 4.0 * pot        ! LJ units (sigma = 1)
     vir = 48.0 * vir / 3.0 ! LJ units (sigma = 1)
 
-  END SUBROUTINE energy_i
+  END SUBROUTINE energy_1
 
   FUNCTION pot_lrc ( sigma, r_cut, density ) ! Long-range correction to potential per atom
     REAL             :: pot_lrc              ! Function result in LJ (sigma=1) units

md_hard.f90

@@ -1,47 +1,51 @@
 ! md_hard.f90 (also uses md_hard_module.f90 and io_module.f90)
 ! Molecular dynamics of hard spheres
 PROGRAM md_hard
+  USE utility_module
   USE md_hard_module
-  USE io_module
   IMPLICIT NONE
 
   ! Takes in a hard-sphere configuration (positions and velocities)
   ! Checks for overlaps    
-  ! Conducts molecular dynamics simulation for specified number of collisions
+  ! Conducts molecular dynamics simulation
   ! Uses no special neighbour lists
   ! ... so is restricted to small number of atoms
   ! Assumes that collisions can be predicted by looking at 
   ! nearest neighbour particles in periodic boundaries
   ! ... so is unsuitable for low densities
-  ! Box is taken to be of unit length.                       
-  ! However, input configuration, output configuration and all
-  ! results are given in units sigma = 1, kT = 1, mass = 1
+  ! Box is taken to be of unit length during the dynamics.
+  ! However, input configuration, output configuration,
+  ! most calculations, and all results 
+  ! are given in units sigma = 1, mass = 1
 
   REAL                        :: sigma   ! atomic diameter (in units where box=1)
   REAL                        :: box     ! box length (in units where sigma=1)
   REAL                        :: density ! reduced density n*sigma**3/box**3
   CHARACTER(len=7), PARAMETER :: prefix = 'md_hard'
   CHARACTER(len=7), PARAMETER :: cnfinp = '.cnfinp', cnfout = '.cnfout'
 
-  INTEGER            :: i, j, k, ncoll, coll
-  REAL               :: tij, t, rate, kinetic_energy
+  INTEGER            :: i, j, k, ncoll, coll, nblock
+  REAL               :: tij, t, rate, kin
   REAL               :: virial, pvnkt1, virial_avg, temperature, tbc, sigma_sq
   REAL, DIMENSION(3) :: total_momentum
 
-  NAMELIST /run_parameters/ ncoll
+  NAMELIST /run_parameters/ nblock, ncoll
 
   WRITE(*,'(''md_hard'')')
   WRITE(*,'(''Molecular dynamics of hard spheres'')')
-  WRITE(*,'(''Results in units kT = sigma = 1'')')
+  WRITE(*,'(''Results in units sigma = 1, mass = 1'')')
 
+  ! Set sensible defaults for testing
+  nblock = 10
+  ncoll  = 10000
   READ(*,nml=run_parameters)
-  WRITE(*,'(''Collisions required'',t40,i15)'  ) ncoll
+  WRITE(*,'(''Number of blocks'',              t40,i15)'  ) nblock
+  WRITE(*,'(''Number of collisions per block'',t40,i15)'  ) ncoll
 
   CALL read_cnf_atoms ( prefix//cnfinp, n, box )
-  WRITE(*,'(''Number of particles'',t40,i15)') n
+  WRITE(*,'(''Number of particles'', t40,i15  )') n
   WRITE(*,'(''Box (in sigma units)'',t40,f15.5)') box
   sigma = 1.0
-  WRITE(*,'(''Sigma (in sigma units)'',t40,f15.5)') sigma
   density = REAL (n) * ( sigma / box ) ** 3
   WRITE(*,'(''Reduced density'',t40,f15.5)') density
 
@@ -51,12 +55,12 @@ PROGRAM md_hard
   total_momentum = SUM(v,dim=2)
   total_momentum = total_momentum / REAL(n)
   WRITE(*,'(''Net momentum'',t40,3f15.5)') total_momentum
-  v = v - SPREAD(total_momentum,dim=1,ncopies=3)
-  kinetic_energy = 0.5 * SUM ( v**2 )
-  temperature = 2.0 * kinetic_energy / REAL ( 3*(n-1) )
-  WRITE(*,'(''Initial temperature (sigma units)'',t40,f15.5)') temperature
+  v   = v - SPREAD(total_momentum,dim=1,ncopies=3)
+  kin = 0.5 * SUM ( v**2 )
+  temperature = 2.0 * kin / REAL ( 3*(n-1) )
+  WRITE(*,'(''Temperature (sigma units)'',t40,f15.5)') temperature
 
-  ! Convert to box units
+  ! Convert to box units (time units are unaffected)
   r(:,:) = r(:,:) / box
   v(:,:) = v(:,:) / box
   r(:,:) = r(:,:) - ANINT ( r(:,:) ) ! Periodic boundaries
@@ -72,7 +76,7 @@ PROGRAM md_hard
   partner(:) = n
 
   DO i = 1, n
-     CALL update ( i, i+1, n, sigma_sq ) ! initial search for collision partners >i
+     CALL update ( i, gt, sigma_sq ) ! initial search for collision partners >i
   END DO
 
   virial_avg = 0.0 ! zero virial accumulator
@@ -98,12 +102,12 @@ PROGRAM md_hard
 
      DO k = 1, n
         IF ( ( k == i ) .OR. ( partner(k) == i ) .OR. ( k == j ) .OR. ( partner(k) == j ) ) THEN
-           CALL update ( k, k+1, n, sigma_sq ) ! search for partners >k
+           CALL update ( k, gt, sigma_sq ) ! search for partners >k
         ENDIF
      END DO
 
-     CALL update ( i, 1, i-1, sigma_sq ) ! search for partners <i
-     CALL update ( j, 1, j-1, sigma_sq ) ! search for partners <j
+     CALL update ( i, lt, sigma_sq ) ! search for partners <i
+     CALL update ( j, lt, sigma_sq ) ! search for partners <j
 
   END DO ! End of main loop
 
@@ -114,8 +118,8 @@ PROGRAM md_hard
      WRITE(*,'(''Particle overlap in final configuration'')')
   ENDIF
 
-  kinetic_energy = 0.5 * SUM ( v**2 ) ! in box units
-  temperature = 2.0 * kinetic_energy / REAL ( 3*(n-1) ) ! in box units
+  kin = 0.5 * SUM ( v**2 ) ! in box units
+  temperature = 2.0 * kin / REAL ( 3*(n-1) ) ! in box units
   pvnkt1 = virial_avg / REAL ( n ) / 3.0 / t / temperature ! dimensionless
   t = t * SQRT ( temperature ) / sigma ! in sigma units
   rate = REAL ( ncoll ) / t ! in sigma units

md_hard_module.f90

@@ -4,28 +4,44 @@ MODULE md_hard_module
 
   IMPLICIT NONE
   PRIVATE
-  PUBLIC :: n, r, v, coltime, partner
+  PUBLIC :: n, r, v, coltime, partner, lt, ne, gt
   PUBLIC :: update, overlap, collide
 
   INTEGER                              :: n       ! number of atoms
   REAL,    DIMENSION(:,:), ALLOCATABLE :: r, v    ! positions, velocities (3,n)
   REAL,    DIMENSION(:),   ALLOCATABLE :: coltime ! time to next collision (n)
   INTEGER, DIMENSION(:),   ALLOCATABLE :: partner ! collision partner (n)
 
+  INTEGER, PARAMETER :: lt = -1, ne = 0, gt = 1 ! j-range options
+
 CONTAINS
 
-  SUBROUTINE update ( i, j1, j2, sigma_sq ) ! updates collision details for atom i
-    INTEGER, INTENT(in) :: i, j1, j2
+  SUBROUTINE update ( i, j_range, sigma_sq ) ! updates collision details for atom i
+    INTEGER, INTENT(in) :: i, j_range
     REAL,    INTENT(in) :: sigma_sq
 
-    INTEGER            :: j
+    INTEGER            :: j, j1, j2
     REAL, DIMENSION(3) :: rij, vij
     REAL               :: rijsq, vijsq, bij, tij, discr
 
+    SELECT CASE ( j_range )
+    CASE ( lt ) ! j < i
+       j1 = 1
+       j2 = i-1
+    CASE ( gt ) ! j > i
+       j1 = i+1
+       j2 = n
+    CASE ( ne ) ! j /= i
+       j1 = 1
+       j2 = n
+    END SELECT
+
     coltime(i) = HUGE(1.0)
 
     DO j = j1, j2
 
+       IF ( i == j ) CYCLE
+
        rij(:) = r(:,i) - r(:,j)
        rij(:) = rij(:) - ANINT ( rij(:) )
        vij(:) = v(:,i) - v(:,j)