diff --git a/README.md b/README.md index b7136f6..44c152a 100644 --- a/README.md +++ b/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. diff --git a/SConscript b/SConscript new file mode 100644 index 0000000..5ff4b5b --- /dev/null +++ b/SConscript @@ -0,0 +1,2 @@ +Import('env','sources') +env.Program(sources) diff --git a/SConstruct b/SConstruct new file mode 100644 index 0000000..718fe8e --- /dev/null +++ b/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) diff --git a/mc_nvt_lj.f90 b/mc_nvt_lj.f90 index 7cf623a..5e986b6 100644 --- a/mc_nvt_lj.f90 +++ b/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 diff --git a/mc_nvt_lj_module.f90 b/mc_nvt_lj_module.f90 index 809ef5a..7d4bd0f 100644 --- a/mc_nvt_lj_module.f90 +++ b/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 diff --git a/md_hard.f90 b/md_hard.f90 index 509cbd1..a1a3425 100644 --- a/md_hard.f90 +++ b/md_hard.f90 @@ -1,21 +1,22 @@ ! 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) @@ -23,25 +24,28 @@ PROGRAM md_hard 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 + 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)