diff --git a/CONTENTS.md b/CONTENTS.md
index a0dd375..313c5c9 100644
--- a/CONTENTS.md
+++ b/CONTENTS.md
@@ -26,6 +26,11 @@ but are mentioned in the online [GUIDE](GUIDE.md):
* [eos_lj.f90](eos_lj.f90) with [eos_lj_module.f90](eos_lj_module.f90)
* [eos_hs.f90](eos_hs.f90)
+An additional molecular dynamics program using quaternions
+[md_nvt_poly_lj.f90](md_nvt_poly_lj.f90) and [md_poly_lj_module.f90](md_poly_lj_module.f90)
+is described in the online [GUIDE](GUIDE.md),
+but does not appear in the text.
+
### 1.1 Calculation of T tensors
[t_tensor.f90](t_tensor.f90).
@@ -100,6 +105,9 @@ We also supply a constant-pressure version [mc_npt_hs.f90](mc_npt_hs.f90).
### 4.7 Monte Carlo program using quaternions
[mc_nvt_poly_lj.f90](mc_nvt_poly_lj.f90) and [mc_poly_lj_module.f90](mc_poly_lj_module.f90).
+We also supply a molecular dynamics program
+[md_nvt_poly_lj.f90](md_nvt_poly_lj.f90) and [md_poly_lj_module.f90](md_poly_lj_module.f90)
+to simulate the same model.
### 4.8 Monte Carlo of hard spherocylinders
[mc_nvt_sc.f90](mc_nvt_sc.f90) and [mc_sc_module.f90](mc_sc_module.f90).
diff --git a/GUIDE.md b/GUIDE.md
index 073965b..c9a7730 100644
--- a/GUIDE.md
+++ b/GUIDE.md
@@ -876,6 +876,7 @@ _T_ | _E_ | _P_ | _T_ (K) | _PE_ (kJ/mol) | _PE_ (kJ/mol) eqn (23
A second set of tests was performed at _T_=0.6≈380K
at the specified densities ρ1, … ρ5 of Mossa et al (2002).
+A set of starting configurations is provided in the [Data repository](https://github.com/Allen-Tildesley/data).
Here the excess pressure (_P_(ex)=_P_-ρ_T_ converted to MPa
with a factor 77.75 based on the values of ε and σ)
is compared with the fit given by eqn (28) and the coefficients in Table III of Mossa et al (2002).
@@ -898,6 +899,71 @@ since this system can show sluggish behaviour.
The MD simulations of Mossa et al (2002) are reported to extend to several hundred nanoseconds
(of order 107 MD timesteps) at the lowest temperatures.
+For comparison we provide a molecular dynamics code `md_nvt_poly_lj` for the same model.
+The program takes the molecular mass _M_ to be unity.
+Mossa et al (2002) ascribe a notional mass of 78u to each of the three LJ sites,
+so _M_≈3.9×10-25kg.
+Combined with the above values of ε and σ,
+this gives a time scale (_M_/ε)1/2σ ≈ 3.22 ps.
+The timestep of δt=0.01 ps used by Mossa et al (2002)
+corresponds to the default value in the program `dt=0.003` in these units.
+By default, the program simulates the constant-_NVE_ ensemble,
+but there is an option to simulate at constant _NVT_ by velocity randomization (Andersen thermostat).
+If the latter option is selected,
+the program will read configurations in the same format as `mc_nvt_poly_lj` (positions and quaternions only),
+selecting random initial velocities and angular momenta,
+which can be convenient.
+
+By default the program calculates the inertia tensor from the LJ site bond vectors,
+assuming equal masses.
+For simplicity it is assumed that the bond vectors are defined such that
+the principal axes of the inertia tensor coincide with
+the xyz axes of the molecular coordinate system,
+with the centre of mass at the origin;
+it is always possible to arrange this.
+In general,
+the three principal moments of inertia will all be different,
+so the molecule is an asymmetric top.
+The MD algorithm for rotation is a symplectic one
+in which a `kick` propagator advances the space-fixed angular momenta,
+using the torque on each molecule,
+and a succession of `drift` steps implement free rotation about each of the principal axes.
+This is described in the text, section 3.3; see
+
+* A Dullweber, B Leimkuhler, R McLachlan, _J Chem Phys,_ __107,__ 5840 (1997),
+* TF Miller, M Eleftheriou, P Pattnaik, A Ndirango, D Newns, GJ Martyna, _J Chem Phys,_ __116,__ 8649 (2002).
+
+The results below are for test runs in both constant-_NVE_ and constant-_NVT_ ensembles,
+at (approximately) the same state points as those given above.
+All runs were 10×20000 steps in length and used program defaults,
+except for `t_interval=1` and the specified temperature in the _NVT_ case.
+For constant-_NVE_ runs we report RMS energy fluctuations,
+and _T_ is the average translational temperature.
+
+ ρ | _T_ | _E_ | _P_ | _E_(RMS)
+ ----- | ----- | ----- | ----- | -----
+ 0.32655 | 0.5 | -12.984(5) | 1.81(1) |
+ 0.32655 | 0.5082(4) | -12.9838 | 1.755(6) | 1.23×10-8
+ 0.32655 | 1.0 | -9.80(2) | 5.87(4) |
+ 0.32655 | 1.004(1) | -9.8006 | 5.896(5) | 9.88×10-8
+ 0.32655 | 1.5 | -6.83(1) | 9.35(3) |
+ 0.32655 | 1.506(1) | -6.8326 | 9.378(6) | 3.67×10-7
+ 0.32655 | 2.0 | -4.05(1) | 12.42(4) |
+ 0.32655 | 2.007(1) | -4.0507 | 12.405(4) | 9.39×10-7
+
+ ρ | _T_ | _E_ | _P_ | _E_ (RMS)
+ ----- | ----- | ----- | ----- | -----
+ 0.30533 | 0.6 | -11.613(5) | 0.45(2) |
+ 0.30533 | 0.6013(8) | -11.6131 | 0.485(4) | 1.31×10-8
+ 0.31240 | 0.6 | -11.892(6) | 1.04(2) |
+ 0.31240 | 0.6039(8) | -11.8923 | 1.058(5) | 1.52×10-8
+ 0.31918 | 0.6 | -12.147(4) | 1.75(1) |
+ 0.31918 | 0.603(1) | -12.1465 | 1.710(9) | 1.73×10-8
+ 0.32655 | 0.6 | -12.362(3) | 2.61(1) |
+ 0.32655 | 0.604(2) | -12.3616 | 2.58(1) | 2.05×10-8
+ 0.33451 | 0.6 | -12.453(7) | 3.96(2) |
+ 0.33451 | 0.612(1) | -12.4532 | 3.87(1) | 2.53×10-8
+
## DPD program
For the `dpd` example, we recommend generating an initial configuration
using the `initialize` program, with namelist input similar to the following
diff --git a/SConstruct b/SConstruct
index 41b5083..9cf1f8b 100644
--- a/SConstruct
+++ b/SConstruct
@@ -98,6 +98,7 @@ variants['build_md_nvt_lj'] = (['md_nvt_lj.f90','md_lj_module.f90','l
variants['build_md_nvt_lj_ll'] = (['md_nvt_lj.f90','md_lj_ll_module.f90','lrc_lj_module.f90','link_list_module.f90']+utilities,env_normal)
variants['build_md_nvt_lj_le'] = (['md_nvt_lj_le.f90','md_lj_le_module.f90']+utnomaths,env_normal)
variants['build_md_nvt_lj_llle'] = (['md_nvt_lj_le.f90','md_lj_llle_module.f90','link_list_module.f90']+utnomaths,env_normal)
+variants['build_md_nvt_poly_lj'] = (['md_nvt_poly_lj.f90','md_poly_lj_module.f90']+utilities,env_normal)
variants['build_mesh'] = (['mesh.f90','mesh_module.f90'],env_normal)
variants['build_pair_distribution'] = (['pair_distribution.f90','config_io_module.f90'],env_normal)
variants['build_qmc_pi_sho'] = (['qmc_pi_sho.f90','averages_module.f90','maths_module.f90'],env_normal)
diff --git a/config_io_module.f90 b/config_io_module.f90
index 457374b..425202b 100644
--- a/config_io_module.f90
+++ b/config_io_module.f90
@@ -183,11 +183,11 @@ SUBROUTINE read_cnf_mols ( filename, n, box, r, e, v, w ) ! Read in molecular co
REAL, DIMENSION(:,:), OPTIONAL, INTENT(out) :: r ! Molecular positions (3,n)
REAL, DIMENSION(:,:), OPTIONAL, INTENT(out) :: e ! Orientation vectors (3,n) or quaternions (4,n)
REAL, DIMENSION(:,:), OPTIONAL, INTENT(out) :: v ! Molecular velocities (3,n)
- REAL, DIMENSION(:,:), OPTIONAL, INTENT(out) :: w ! Angular velocities (3,n)
+ REAL, DIMENSION(:,:), OPTIONAL, INTENT(out) :: w ! Angular velocities/momenta (3,n)
! The first call of this routine is just to get n and box, used to allocate arrays
! The second call attempts to read in the molecular positions, orientations
- ! and optionally velocities, angular velocities
+ ! and optionally velocities, angular velocities/momenta
INTEGER :: cnf_unit, ioerr, i
@@ -252,7 +252,7 @@ SUBROUTINE read_cnf_mols ( filename, n, box, r, e, v, w ) ! Read in molecular co
STOP 'Error in read_cnf_mols'
END IF
- ! Read positions, orientation vectors or quaternions, velocities, angular velocities
+ ! Read positions, orientation vectors or quaternions, velocities, angular velocities/momenta
DO i = 1, n
READ ( unit=cnf_unit, fmt=*, iostat=ioerr ) r(:,i), e(:,i), v(:,i), w(:,i)
IF ( ioerr /= 0 ) THEN
@@ -292,7 +292,7 @@ SUBROUTINE write_cnf_mols ( filename, n, box, r, e, v, w ) ! Write out molecular
REAL, DIMENSION(:,:), INTENT(in) :: r ! Molecular positions (3,n)
REAL, DIMENSION(:,:), INTENT(in) :: e ! Orientation vectors (3,n) or quaternions (0:3,n)
REAL, DIMENSION(:,:), OPTIONAL, INTENT(in) :: v ! Molecular velocities (3,n)
- REAL, DIMENSION(:,:), OPTIONAL, INTENT(in) :: w ! Angular velocities (3,n)
+ REAL, DIMENSION(:,:), OPTIONAL, INTENT(in) :: w ! Angular velocities/momenta (3,n)
INTEGER :: cnf_unit, ioerr, i
@@ -330,7 +330,7 @@ SUBROUTINE write_cnf_mols ( filename, n, box, r, e, v, w ) ! Write out molecular
STOP 'Error in write_cnf_mols'
END IF
- ! Write positions, orientation vectors or quaternions, velocities, angular velocities
+ ! Write positions, orientation vectors or quaternions, velocities, angular velocities/momenta
DO i = 1, n
WRITE ( unit=cnf_unit, fmt='(*(f15.10))') r(:,i), e(:,i), v(:,i), w(:,i)
END DO
diff --git a/mc_poly_lj_module.f90 b/mc_poly_lj_module.f90
index f0b985f..561d3f1 100644
--- a/mc_poly_lj_module.f90
+++ b/mc_poly_lj_module.f90
@@ -61,7 +61,7 @@ MODULE mc_module
! Public derived type
TYPE, PUBLIC :: potential_type ! A composite variable for interactions comprising
- REAL :: pot ! the potential energy cut and shifted to 0 at r_cut, and
+ REAL :: pot ! the potential energy
REAL :: vir ! the virial and
LOGICAL :: ovr ! a flag indicating overlap (i.e. pot too high to use)
CONTAINS
diff --git a/md_nvt_poly_lj.f90 b/md_nvt_poly_lj.f90
new file mode 100644
index 0000000..f069a45
--- /dev/null
+++ b/md_nvt_poly_lj.f90
@@ -0,0 +1,362 @@
+! md_nvt_poly_lj.f90
+! Molecular dynamics, NVT ensemble (NVE option), polyatomic molecules
+PROGRAM md_nvt_poly_lj
+
+ !------------------------------------------------------------------------------------------------!
+ ! This software was written in 2016/17 !
+ ! by Michael P. Allen / !
+ ! and Dominic J. Tildesley ("the authors"), !
+ ! to accompany the book "Computer Simulation of Liquids", second edition, 2017 ("the text"), !
+ ! published by Oxford University Press ("the publishers"). !
+ ! !
+ ! LICENCE !
+ ! Creative Commons CC0 Public Domain Dedication. !
+ ! To the extent possible under law, the authors have dedicated all copyright and related !
+ ! and neighboring rights to this software to the PUBLIC domain worldwide. !
+ ! This software is distributed without any warranty. !
+ ! You should have received a copy of the CC0 Public Domain Dedication along with this software. !
+ ! If not, see . !
+ ! !
+ ! DISCLAIMER !
+ ! The authors and publishers make no warranties about the software, and disclaim liability !
+ ! for all uses of the software, to the fullest extent permitted by applicable law. !
+ ! The authors and publishers do not recommend use of this software for any purpose. !
+ ! It is made freely available, solely to clarify points made in the text. When using or citing !
+ ! the software, you should not imply endorsement by the authors or publishers. !
+ !------------------------------------------------------------------------------------------------!
+
+ ! Takes in a configuration of atoms (positions, quaternions, velocities and angular momenta)
+ ! Cubic periodic boundary conditions
+ ! Conducts molecular dynamics with optional velocity thermalization
+ ! Uses no special neighbour lists
+ ! The rotational algorithm is described in the text, section 3.3. See A Dullweber, B Leimkuhler,
+ ! R McLachlan, J Chem Phys 107, 5840 (1997) and TF Miller, M Eleftheriou, P Pattnaik,
+ ! A Ndirango, D Newns, GJ Martyna, J Chem Phys 116, 8649 (2002).
+
+ ! Reads several variables and options from standard input using a namelist nml
+ ! Leave namelist empty to accept supplied defaults
+
+ ! Positions r are divided by box length after reading in and we assume mass=1 throughout
+ ! However, input configuration, output configuration, most calculations, and all results
+ ! are given in simulation units defined by the model
+ ! For example, for Lennard-Jones, sigma = 1, epsilon = 1
+
+ ! Despite the program name, there is nothing here specific to Lennard-Jones
+ ! The model is defined in md_module
+
+ USE, INTRINSIC :: iso_fortran_env, ONLY : input_unit, output_unit, error_unit, iostat_end, iostat_eor
+ USE config_io_module, ONLY : read_cnf_mols, write_cnf_mols
+ USE averages_module, ONLY : run_begin, run_end, blk_begin, blk_end, blk_add
+ USE maths_module, ONLY : q_to_a
+ USE md_module, ONLY : introduction, conclusion, allocate_arrays, deallocate_arrays, &
+ & force, r, e, d, v, ell, f, tau, n, na, db, inertia, potential_type
+
+ IMPLICIT NONE
+
+ ! Most important variables
+ REAL :: box ! Box length
+ REAL :: dt ! Time step
+ REAL :: temperature ! Specified temperature
+ LOGICAL :: nvt ! Flag to indicate NVT ensemble
+
+ ! Composite interaction = pot & vir & ovr variables
+ TYPE(potential_type) :: total
+
+ INTEGER :: blk, stp, nstep, nblock, t_interval, ioerr, i, a
+ REAL :: norm
+ REAL, DIMENSION(3) :: vcm
+ REAL, DIMENSION(3,3) :: ai
+
+ CHARACTER(len=4), PARAMETER :: cnf_prefix = 'cnf.'
+ 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
+
+ NAMELIST /nml/ nblock, nstep, dt, temperature, t_interval
+
+ WRITE ( unit=output_unit, fmt='(a)' ) 'md_nvt_poly_lj'
+ WRITE ( unit=output_unit, fmt='(a)' ) 'Molecular dynamics, constant-NVT/NVE ensemble'
+ WRITE ( unit=output_unit, fmt='(a)' ) 'Molecular mass=1 throughout'
+
+ CALL introduction
+ CALL RANDOM_SEED () ! Initialize random number generator
+
+ ! Set sensible default run parameters for testing
+ nblock = 10
+ nstep = 20000
+ dt = 0.003
+ temperature = 1.0
+ t_interval = 0
+
+ ! Read run parameters from namelist
+ ! Comment out, or replace, this section if you don't like namelists
+ READ ( unit=input_unit, nml=nml, iostat=ioerr )
+ IF ( ioerr /= 0 ) THEN
+ WRITE ( unit=error_unit, fmt='(a,i15)') 'Error reading namelist nml from standard input', ioerr
+ IF ( ioerr == iostat_eor ) WRITE ( unit=error_unit, fmt='(a)') 'End of record'
+ IF ( ioerr == iostat_end ) WRITE ( unit=error_unit, fmt='(a)') 'End of file'
+ STOP 'Error in md_poly_lj'
+ END IF
+
+ ! Write out run parameters
+ nvt = t_interval > 0 .AND. t_interval < nstep
+ WRITE ( unit=output_unit, fmt='(a,t40,i15)' ) 'Number of blocks', nblock
+ WRITE ( unit=output_unit, fmt='(a,t40,i15)' ) 'Number of steps per block', nstep
+ WRITE ( unit=output_unit, fmt='(a,t40,f15.6)' ) 'Time step', dt
+ IF ( nvt ) THEN
+ WRITE ( unit=output_unit, fmt='(a)' ) 'NVT ensemble'
+ WRITE ( unit=output_unit, fmt='(a,t40,i15)' ) 'Thermalization interval', t_interval
+ WRITE ( unit=output_unit, fmt='(a,t40,f15.6)' ) 'Temperature', temperature
+ ELSE
+ WRITE ( unit=output_unit, fmt='(a)' ) 'NVE ensemble'
+ t_interval = nstep+1
+ END IF
+
+ ! Read in initial configuration and allocate necessary arrays
+ CALL read_cnf_mols ( cnf_prefix//inp_tag, n, box ) ! First call is just to get n and box
+ WRITE ( unit=output_unit, fmt='(a,t40,i15)' ) 'Number of particles', n
+ WRITE ( unit=output_unit, fmt='(a,t40,f15.6)' ) 'Simulation box length', box
+ WRITE ( unit=output_unit, fmt='(a,t40,f15.6)' ) 'Density', REAL(n) / box**3
+ CALL allocate_arrays ( box )
+ IF ( nvt ) THEN
+ CALL read_cnf_mols ( cnf_prefix//inp_tag, n, box, r, e ) ! Second call gets r, e
+ CALL ran_velocities
+ ELSE
+ CALL read_cnf_mols ( cnf_prefix//inp_tag, n, box, r, e, v, ell ) ! Second call gets r, e, v and ell
+ vcm(:) = SUM ( v(:,:), dim=2 ) / REAL(n) ! Centre-of mass velocity
+ v(:,:) = v(:,:) - SPREAD ( vcm(:), dim = 2, ncopies = n ) ! Set COM velocity to zero
+ END IF
+ r(:,:) = r(:,:) / box ! Convert positions to box units
+ r(:,:) = r(:,:) - ANINT ( r(:,:) ) ! Periodic boundaries
+
+ ! Calculate all bond vectors
+ DO i = 1, n ! Loop over molecules
+ norm = SQRT(SUM(e(:,i)**2))
+ e(:,i) = e(:,i) / norm ! Ensure normalized quaternions
+ ai = q_to_a ( e(:,i) ) ! Rotation matrix for i
+ DO a = 1, na ! Loop over all atoms
+ d(:,a,i) = MATMUL ( db(:,a), ai ) ! NB: equivalent to ai_T*db, ai_T=transpose of ai
+ END DO ! End loop over all atoms
+ END DO ! End loop over molecules
+
+ ! Initial forces, potential, etc plus overlap check
+ CALL force ( box, total )
+ IF ( total%ovr ) THEN
+ WRITE ( unit=error_unit, fmt='(a)') 'Overlap in initial configuration'
+ STOP 'Error in md_poly_lj'
+ END IF
+
+ ! Initialize arrays for averaging and write column headings
+ CALL run_begin ( calc_variables() )
+
+ DO blk = 1, nblock ! Begin loop over blocks
+
+ CALL blk_begin
+
+ DO stp = 1, nstep ! Begin loop over steps
+
+ IF ( nvt .AND. MOD(stp,t_interval) == 0 ) THEN
+ CALL ran_velocities
+ END IF
+
+ CALL kick_propagator ( 0.5 * dt ) ! Half-kick step
+
+ CALL drift_translate ( dt ) ! Drift step for positions
+
+ ! Succession of drift steps for rotation about body-fixed axes
+ ! Depending on the values of the moments of inertia, a different nested
+ ! sequence of axes may produce better or worse energy conservation
+ CALL drift_rotate ( 1, 0.5*dt )
+ CALL drift_rotate ( 2, 0.5*dt )
+ CALL drift_rotate ( 3, dt )
+ CALL drift_rotate ( 2, 0.5*dt )
+ CALL drift_rotate ( 1, 0.5*dt )
+
+ ! Calculate all bond vectors
+ DO i = 1, n ! Loop over molecules
+ norm = SQRT(SUM(e(:,i)**2))
+ e(:,i) = e(:,i) / norm ! Guard against cumulative roundoff
+ ai = q_to_a ( e(:,i) ) ! Rotation matrix for i
+ DO a = 1, na ! Loop over all atoms
+ d(:,a,i) = MATMUL ( db(:,a), ai ) ! NB: equivalent to ai_T*db, ai_T=transpose of ai
+ END DO ! End loop over all atoms
+ END DO ! End loop over molecules
+
+ CALL force ( box, total ) ! Force evaluation
+ IF ( total%ovr ) THEN
+ WRITE ( unit=error_unit, fmt='(a)') 'Overlap in configuration'
+ STOP 'Error in md_poly_lj'
+ END IF
+
+ CALL kick_propagator ( 0.5 * dt ) ! Half-kick step
+
+ ! Calculate and accumulate variables for this step
+ CALL blk_add ( calc_variables() )
+
+ END DO ! End loop over steps
+
+ CALL blk_end ( blk ) ! Output block averages
+ IF ( nblock < 1000 ) WRITE(sav_tag,'(i3.3)') blk ! Number configuration by block
+ CALL write_cnf_mols ( cnf_prefix//sav_tag, n, box, r*box, e, v, ell ) ! Save configuration
+
+ END DO ! End loop over blocks
+
+ CALL run_end ( calc_variables() ) ! Output run averages
+
+ CALL force ( box, total )
+ IF ( total%ovr ) THEN ! should never happen
+ WRITE ( unit=error_unit, fmt='(a)') 'Overlap in final configuration'
+ STOP 'Error in md_poly_lj'
+ END IF
+
+ CALL write_cnf_mols ( cnf_prefix//out_tag, n, box, r*box, e, v, ell ) ! Write out final configuration
+
+ CALL deallocate_arrays
+ CALL conclusion
+
+CONTAINS
+
+ SUBROUTINE kick_propagator ( t )
+ IMPLICIT NONE
+ REAL, INTENT(in) :: t ! Time over which to propagate (typically dt/2)
+
+ ! Advances velocities and body-fixed angular momenta
+ v(:,:) = v(:,:) + t * f(:,:) ! Linear momenta are equivalent to linear velocities
+ ell(:,:) = ell(:,:) + t * tau(:,:) ! Space-fixed angular momenta ell and torques tau
+
+ END SUBROUTINE kick_propagator
+
+ SUBROUTINE drift_translate ( t )
+ IMPLICIT NONE
+ REAL, INTENT(in) :: t ! Time over which to propagate (typically dt)
+
+ ! Advances positions
+
+ r(:,:) = r(:,:) + t * v(:,:) / box ! Drift step (positions in box=1 units)
+ r(:,:) = r(:,:) - ANINT ( r(:,:) ) ! Periodic boundaries
+
+ END SUBROUTINE drift_translate
+
+ SUBROUTINE drift_rotate ( xyz, t )
+ USE maths_module, ONLY : q_to_a, rotate_quaternion
+ IMPLICIT NONE
+ INTEGER, INTENT(in) :: xyz ! Body-fixed axis about which to rotate
+ REAL, INTENT(in) :: t ! Time over which to propagate (typically dt or dt/2)
+
+ ! Advances quaternion orientations about a specified axis
+
+ INTEGER :: i
+ REAL :: w_mag
+ REAL, DIMENSION(3) :: w_hat
+ REAL, DIMENSION(3,3) :: ai
+
+ DO i = 1, n
+ ai = q_to_a ( e(:,i) ) ! Rotation matrix
+ w_hat = ai(xyz,:) ! Space-fixed axis about which to rotate
+ w_mag = DOT_PRODUCT ( ell(:,i), w_hat ) / inertia(xyz) ! Angular velocity about this axis
+ e(:,i) = rotate_quaternion ( w_mag*t, w_hat, e(:,i) ) ! Rotate by specified angle
+ END DO
+
+ END SUBROUTINE drift_rotate
+
+ SUBROUTINE ran_velocities
+ USE maths_module, ONLY : random_normals
+
+ INTEGER :: i
+ REAL :: factor
+ REAL, DIMENSION(3) :: v_cm, factors
+ REAL, DIMENSION(3,3) :: ai
+
+ CALL random_normals ( 0.0, 1.0, v ) ! Random velocities
+ v_cm(:) = SUM ( v(:,:), dim=2 ) / REAL ( n ) ! Compute centre of mass velocity
+ v(:,:) = v(:,:) - SPREAD ( v_cm(:), dim=2, ncopies=n ) ! Set net momentum to zero
+ factor = SUM ( v**2 ) / REAL(3*n-3) ! Estimate of kinetic temperature
+ factor = SQRT ( temperature / factor ) ! Necessary correction factor
+ v = factor * v ! Make correction
+
+ CALL random_normals ( 0.0, 1.0, ell ) ! Random body-fixed angular momenta
+ factors = SUM ( ell**2, dim=2 ) / (REAL(n)*inertia(:)) ! Estimate of kinetic temperatures
+ factors = SQRT ( temperature / factors ) ! Necessary correction factors
+ ell = SPREAD(factors,dim=2,ncopies=n) * ell ! Make corrections
+
+ ! Convert to space-fixed angular momenta
+ DO i = 1, n ! Loop over molecules
+ ai = q_to_a ( e(:,i) ) ! Rotation matrix for i
+ ell(:,i) = MATMUL ( ell(:,i), ai ) ! NB: equivalent to ell_s = ai_T*ell_b, ai_T=transpose of ai
+ END DO ! End loop over molecules
+
+ END SUBROUTINE ran_velocities
+
+ FUNCTION calc_variables ( ) RESULT ( variables )
+ USE averages_module, ONLY : variable_type, msd
+ USE maths_module, ONLY : q_to_a
+ IMPLICIT NONE
+ TYPE(variable_type), DIMENSION(5) :: variables ! The 5 variables listed below
+
+ ! This routine calculates all variables of interest and (optionally) writes them out
+ ! They are collected together in the variables array, for use in the main program
+
+ TYPE(variable_type) :: e_sf, p_sf, t_t, t_r, conserved_msd
+ REAL :: vol, rho, kin_t, kin_r, eng, tmp_t, tmp_r
+ INTEGER :: i
+ REAL, DIMENSION(3) :: ell_i
+ REAL, DIMENSION(3,3) :: ai
+
+ ! Preliminary calculations
+ vol = box**3 ! Volume
+ rho = REAL(n) / vol ! Density
+ kin_t = 0.5*SUM(v**2) ! Translational kinetic energy
+
+ kin_r = 0.0
+ DO i = 1, n ! Loop over molecules
+ ai = q_to_a ( e(:,i) ) ! Rotation matrix for i
+ ell_i = MATMUL ( ai, ell(:,i) ) ! Get body-fixed angular momentum
+ kin_r = kin_r + SUM((ell_i**2)/inertia)
+ END DO ! End loop over molecules
+ kin_r = 0.5*kin_r ! Rotational kinetic energy
+
+ tmp_t = 2.0 * kin_t / REAL(3*n-3) ! Remove three degrees of freedom for momentum conservation
+ tmp_r = 2.0 * kin_r / REAL(3*n) ! 3N degrees of rotational freedom
+ eng = kin_t + kin_r + total%pot ! Total energy for simulated system
+
+ ! Variables of interest, of type variable_type, containing three components:
+ ! %val: the instantaneous value
+ ! %nam: used for headings
+ ! %method: indicating averaging method
+ ! If not set below, %method adopts its default value of avg
+ ! The %nam and some other components need only be defined once, at the start of the program,
+ ! but for clarity and readability we assign all the values together below
+
+ ! Internal energy (shifted-force potential) per atom
+ ! Total translational and rotational KE plus total PE divided by N
+ IF ( nvt ) THEN
+ e_sf = variable_type ( nam = 'E/N shifted force', val = 3.0*temperature+total%pot/REAL(n) )
+ ELSE
+ e_sf = variable_type ( nam = 'E/N shifted force', val = eng/REAL(n) )
+ END IF
+
+ ! Pressure (shifted-force potential)
+ ! Ideal gas contribution plus total virial divided by V
+ IF ( nvt ) THEN
+ p_sf = variable_type ( nam = 'P shifted force', val = rho*temperature + total%vir/vol )
+ ELSE
+ p_sf = variable_type ( nam = 'P shifted force', val = rho*tmp_t + total%vir/vol )
+ END IF
+
+ ! Kinetic translational temperature
+ t_t = variable_type ( nam = 'T translational', val = tmp_t )
+
+ ! Kinetic rotational temperature
+ t_r = variable_type ( nam = 'T rotational', val = tmp_r )
+
+ ! Mean-squared deviation of conserved energy per atom
+ conserved_msd = variable_type ( nam = 'Conserved MSD', val = eng/REAL(n), &
+ & method = msd, e_format = .TRUE., instant = .FALSE. )
+
+ ! Collect together for averaging
+ variables = [ e_sf, p_sf, t_t, t_r, conserved_msd ]
+
+ END FUNCTION calc_variables
+
+END PROGRAM md_nvt_poly_lj
+
diff --git a/md_poly_lj_module.f90 b/md_poly_lj_module.f90
new file mode 100644
index 0000000..6d0586c
--- /dev/null
+++ b/md_poly_lj_module.f90
@@ -0,0 +1,271 @@
+! md_poly_lj_module.f90
+! Force routine for MD simulation, polyatomic molecule, LJ atoms
+MODULE md_module
+
+ !------------------------------------------------------------------------------------------------!
+ ! This software was written in 2016/17 !
+ ! by Michael P. Allen / !
+ ! and Dominic J. Tildesley ("the authors"), !
+ ! to accompany the book "Computer Simulation of Liquids", second edition, 2017 ("the text"), !
+ ! published by Oxford University Press ("the publishers"). !
+ ! !
+ ! LICENCE !
+ ! Creative Commons CC0 Public Domain Dedication. !
+ ! To the extent possible under law, the authors have dedicated all copyright and related !
+ ! and neighboring rights to this software to the PUBLIC domain worldwide. !
+ ! This software is distributed without any warranty. !
+ ! You should have received a copy of the CC0 Public Domain Dedication along with this software. !
+ ! If not, see . !
+ ! !
+ ! DISCLAIMER !
+ ! The authors and publishers make no warranties about the software, and disclaim liability !
+ ! for all uses of the software, to the fullest extent permitted by applicable law. !
+ ! The authors and publishers do not recommend use of this software for any purpose. !
+ ! It is made freely available, solely to clarify points made in the text. When using or citing !
+ ! the software, you should not imply endorsement by the authors or publishers. !
+ !------------------------------------------------------------------------------------------------!
+
+ USE, INTRINSIC :: iso_fortran_env, ONLY : output_unit, error_unit
+
+ IMPLICIT NONE
+ PRIVATE
+
+ ! Public routines
+ PUBLIC :: introduction, conclusion, allocate_arrays, deallocate_arrays
+ PUBLIC :: force
+
+ ! Public data
+ INTEGER, PUBLIC :: n ! number of molecules
+ REAL, DIMENSION(:,:), ALLOCATABLE, PUBLIC :: r ! centre of mass positions (3,n)
+ REAL, DIMENSION(:,:), ALLOCATABLE, PUBLIC :: e ! quaternions (0:3,n)
+ REAL, DIMENSION(:,:,:), ALLOCATABLE, PUBLIC :: d ! bond vectors (0:3,na,n)
+ REAL, DIMENSION(:,:), ALLOCATABLE, PUBLIC :: v ! centre of mass velocities (3,n)
+ REAL, DIMENSION(:,:), ALLOCATABLE, PUBLIC :: ell ! angular momenta (3,n)
+ REAL, DIMENSION(:,:), ALLOCATABLE, PUBLIC :: f ! centre of mass forces (3,n)
+ REAL, DIMENSION(:,:), ALLOCATABLE, PUBLIC :: tau ! torques (3,n)
+
+ ! Bond vectors in body-fixed frame (na and db are public)
+ ! Isosceles triangle, 3 sites, with unit bond length and bond angle alpha, which we set to 75 degrees here
+ REAL, PARAMETER :: pi = 4.0*ATAN(1.0)
+ REAL, PARAMETER :: alpha = 75.0 * pi / 180.0, alpha2 = alpha / 2.0
+ INTEGER, PARAMETER, PUBLIC :: na = 3
+ REAL, DIMENSION(3,na), PARAMETER, PUBLIC :: db = RESHAPE ( [ &
+ & -SIN(alpha2), 0.0, -COS(alpha2)/3.0, &
+ & 0.0, 0.0, 2.0*COS(alpha2)/3.0, &
+ & SIN(alpha2), 0.0, -COS(alpha2)/3.0 ], [3,na] )
+
+ ! Atomic masses and (public) moments of inertia
+ REAL, DIMENSION(na), PARAMETER :: m = [ 1.0/3.0, 1.0/3.0, 1.0/3.0 ] ! Masses add up to 1.0
+ REAL, DIMENSION(3), PUBLIC :: inertia
+
+ ! Cutoff distance and force-shift parameters (all private) chosen as per the reference:
+ ! S Mossa, E La Nave, HE Stanley, C Donati, F Sciortino, P Tartaglia, Phys Rev E, 65, 041205 (2002)
+ REAL, PARAMETER :: r_cut = 2.612 ! in sigma=1 units, where r_cut = 1.2616 nm, sigma = 0.483 nm
+ REAL, PARAMETER :: sr_cut = 1.0/r_cut, sr_cut6 = sr_cut**6, sr_cut12 = sr_cut6**2
+ REAL, PARAMETER :: lambda1 = 4.0*(7.0*sr_cut6-13.0*sr_cut12)
+ REAL, PARAMETER :: lambda2 = -24.0*(sr_cut6-2.0*sr_cut12)*sr_cut
+
+ ! Public derived type
+ TYPE, PUBLIC :: potential_type ! A composite variable for interactions comprising
+ REAL :: pot ! the potential energy
+ REAL :: vir ! the virial and
+ LOGICAL :: ovr ! a flag indicating overlap (i.e. pot too high to use)
+ CONTAINS
+ PROCEDURE :: add_potential_type
+ GENERIC :: OPERATOR(+) => add_potential_type
+ END TYPE potential_type
+
+CONTAINS
+
+ FUNCTION add_potential_type ( a, b ) RESULT (c)
+ IMPLICIT NONE
+ TYPE(potential_type) :: c ! Result is the sum of
+ CLASS(potential_type), INTENT(in) :: a, b ! the two inputs
+ c%pot = a%pot + b%pot
+ c%vir = a%vir + b%vir
+ c%ovr = a%ovr .OR. b%ovr
+ END FUNCTION add_potential_type
+
+ SUBROUTINE introduction
+ IMPLICIT NONE
+
+ INTEGER :: i
+ REAL :: diameter
+ REAL, DIMENSION(3) :: com
+ REAL, PARAMETER :: tol = 1.0e-9
+
+ WRITE ( unit=output_unit, fmt='(a)' ) 'Lennard-Jones potential'
+ WRITE ( unit=output_unit, fmt='(a)' ) 'Cut-and-force-shifted'
+ WRITE ( unit=output_unit, fmt='(a)' ) 'Diameter, sigma = 1'
+ WRITE ( unit=output_unit, fmt='(a)' ) 'Well depth, epsilon = 1'
+
+ WRITE ( unit=output_unit, fmt='(a,t40,i15)' ) 'Number of atoms per molecule', na
+ DO i = 1, na ! Loop over atoms
+ WRITE ( unit=output_unit, fmt='(a,i1,t40,3f15.6)' ) 'Body-fixed atom vector ', i, db(:,i)
+ END DO ! End loop over atoms
+
+ diameter = 2.0 * SQRT ( MAXVAL ( SUM(db**2,dim=1) ) )
+ WRITE ( unit=output_unit, fmt='(a,t40,f15.6)' ) 'Molecular diameter', diameter
+ WRITE ( unit=output_unit, fmt='(a,t40,f15.6)' ) 'r_cut', r_cut
+ WRITE ( unit=output_unit, fmt='(a,t40,f15.6)' ) 'Force-shift lambda1', lambda1
+ WRITE ( unit=output_unit, fmt='(a,t40,f15.6)' ) 'Force-shift lambda2', lambda2
+
+ ! The following section sets the diagonal moments of inertia "realistically"
+ ! based on the values of atomic masses and bond vectors (above), with some checking.
+ ! However, there is nothing to stop the user replacing this section with a statement setting
+ ! the values of inertia(1:3). The masses m are not passed back to the calling program.
+ ! It might be advantageous, for instance, to artificially increase the values in inertia.
+
+ ! Ensure that the db bonds, xyz molecular axes, and masses are chosen such that
+ ! the total mass is 1 and the centre-of-mass is at the origin
+ IF ( ABS ( SUM(m) - 1.0 ) > tol ) THEN
+ WRITE ( unit=error_unit, fmt='(a,f15.6)') 'Molecular mass is not 1.0 ', SUM(m)
+ STOP 'Error in introduction'
+ END IF
+ com = SUM ( SPREAD(m,dim=1,ncopies=3)*db, dim = 2 )
+ IF ( ANY ( ABS(com) > tol ) ) THEN
+ WRITE ( unit=error_unit, fmt='(a,3f15.6)') 'Molecular centre-of-mass error', com
+ STOP 'Error in introduction'
+ END IF
+
+ ! Ensure that the db bonds, xyz molecular axes, and masses are chosen such that
+ ! the off-diagonal elements of the inertia tensor are zero
+ inertia(1) = -SUM ( m*db(2,:)*db(3,:) )
+ inertia(2) = -SUM ( m*db(3,:)*db(1,:) )
+ inertia(3) = -SUM ( m*db(1,:)*db(2,:) )
+ IF ( ANY ( ABS(inertia) > tol ) ) THEN
+ WRITE ( unit=error_unit, fmt='(a,3f15.6)') 'Molecular off-diagonal inertia error', inertia
+ STOP 'Error in introduction'
+ END IF
+
+ ! Calculate the diagonal elements of the inertia tensor
+ inertia(1) = SUM ( m*db(2,:)**2 ) + SUM ( m*db(3,:)**2 )
+ inertia(2) = SUM ( m*db(3,:)**2 ) + SUM ( m*db(1,:)**2 )
+ inertia(3) = SUM ( m*db(1,:)**2 ) + SUM ( m*db(2,:)**2 )
+ WRITE ( unit=output_unit, fmt='(a,t40,3f15.6)' ) 'Inertia Ixx, Iyy, Izz', inertia
+
+ END SUBROUTINE introduction
+
+ SUBROUTINE conclusion
+ IMPLICIT NONE
+
+ WRITE ( unit=output_unit, fmt='(a)') 'Program ends'
+
+ END SUBROUTINE conclusion
+
+ SUBROUTINE allocate_arrays ( box )
+ IMPLICIT NONE
+ REAL, INTENT(in) :: box ! simulation box length
+
+ REAL :: rm_cut_box, diameter
+
+ ALLOCATE ( r(3,n), e(0:3,n), d(3,na,n) )
+ ALLOCATE ( v(3,n), ell(3,n), f(3,n), tau(3,n) )
+
+ diameter = 2.0 * SQRT ( MAXVAL ( SUM(db**2,dim=1) ) )
+ rm_cut_box = ( r_cut+diameter ) / box
+ IF ( rm_cut_box > 0.5 ) THEN
+ WRITE ( unit=error_unit, fmt='(a,f15.6)') 'rm_cut/box too large ', rm_cut_box
+ STOP 'Error in allocate_arrays'
+ END IF
+
+ END SUBROUTINE allocate_arrays
+
+ SUBROUTINE deallocate_arrays
+ IMPLICIT NONE
+
+ DEALLOCATE ( r, e, d )
+ DEALLOCATE ( v, ell, f, tau )
+
+ END SUBROUTINE deallocate_arrays
+
+ SUBROUTINE force ( box, total )
+ USE maths_module, ONLY : cross_product
+ IMPLICIT NONE
+ REAL, INTENT(in) :: box ! Simulation box length
+ TYPE(potential_type), INTENT(out) :: total ! Composite of pot, vir, etc
+
+ ! total%pot is the nonbonded cut-and-shifted potential energy for whole system
+ ! total%vir is the corresponding virial
+ ! total%ovr is a warning flag that there is an overlap
+ ! This routine also calculates forces and torques
+ ! and stores them in the arrays f, tau
+
+ ! It is assumed that r has been divided by box
+ ! The d array of space-fixed bond vectors for each molecule should be computed already
+ ! Results are in LJ units where sigma = 1, epsilon = 1
+ ! Note that this is the force-shifted LJ potential with a linear smoothing term
+ ! S Mossa, E La Nave, HE Stanley, C Donati, F Sciortino, P Tartaglia, Phys Rev E, 65, 041205 (2002)
+
+ INTEGER :: i, j, a, b
+ REAL :: diameter, rm_cut_box, rm_cut_box_sq, r_cut_sq
+ REAL :: sr2, sr6, sr12, rij_sq, rab_sq, virab, rmag
+ REAL, DIMENSION(3) :: rij, rab, fab
+ REAL, PARAMETER :: sr2_ovr = 1.77 ! overlap threshold (pot > 100)
+ TYPE(potential_type) :: pair
+
+ diameter = 2.0 * SQRT ( MAXVAL ( SUM(db**2,dim=1) ) )
+ rm_cut_box = ( r_cut + diameter ) / box ! Molecular cutoff in box=1 units
+ rm_cut_box_sq = rm_cut_box**2 ! squared
+ r_cut_sq = r_cut**2 ! Potential cutoff squared in sigma=1 units
+
+ ! Initialize
+ f = 0.0
+ tau = 0.0
+ total = potential_type ( pot=0.0, vir=0.0, ovr=.FALSE. )
+
+ DO i = 1, n - 1 ! Begin outer loop over molecules
+
+ DO j = i + 1, n ! Begin inner loop over molecules
+
+ rij(:) = r(:,i) - r(:,j) ! Centre-centre separation vector
+ rij(:) = rij(:) - ANINT ( rij(:) ) ! Periodic boundaries in box=1 units
+ rij_sq = SUM ( rij**2 ) ! Squared centre-centre separation in box=1 units
+
+ IF ( rij_sq < rm_cut_box_sq ) THEN ! Test within molecular cutoff
+
+ rij = rij * box ! Now in sigma = 1 units
+
+ ! Double loop over atoms on both molecules
+ DO a = 1, na
+ DO b = 1, na
+
+ rab = rij + d(:,a,i) - d(:,b,j) ! Atom-atom vector, sigma=1 units
+ rab_sq = SUM ( rab**2 ) ! Squared atom-atom separation, sigma=1 units
+
+ IF ( rab_sq < r_cut_sq ) THEN ! Test within potential cutoff
+
+ sr2 = 1.0 / rab_sq ! (sigma/rab)**2
+ pair%ovr = sr2 > sr2_ovr ! Overlap if too close
+
+ rmag = SQRT(rab_sq)
+ sr6 = sr2**3
+ sr12 = sr6**2
+ pair%pot = 4.0*(sr12-sr6) + lambda1 + lambda2*rmag ! LJ atom-atom pair potential (force-shifted)
+ virab = 24.0*(2.0*sr12-sr6 ) - lambda2*rmag ! LJ atom-atom pair virial
+ fab = rab * virab * sr2 ! LJ atom-atom pair force
+ pair%vir = DOT_PRODUCT ( rij, fab ) ! Contribution to molecular virial
+
+ total = total + pair
+ f(:,i) = f(:,i) + fab
+ f(:,j) = f(:,j) - fab
+ tau(:,i) = tau(:,i) + cross_product ( d(:,a,i), fab )
+ tau(:,j) = tau(:,j) - cross_product ( d(:,b,j), fab )
+
+ END IF ! End test within potential cutoff
+
+ END DO
+ END DO
+ ! End double loop over atoms on both molecules
+
+ END IF ! End test within molecular cutoff
+
+ END DO ! End inner loop over molecules
+ END DO ! End outer loop over molecules
+
+ ! Include numerical factors
+ total%vir = total%vir / 3.0 ! Divide virial by 3
+
+ END SUBROUTINE force
+
+END MODULE md_module
diff --git a/python_examples/CONTENTS.md b/python_examples/CONTENTS.md
index 71de1d7..e303cdb 100644
--- a/python_examples/CONTENTS.md
+++ b/python_examples/CONTENTS.md
@@ -26,6 +26,11 @@ but are mentioned in the online [GUIDE](GUIDE.md):
* [eos_lj.py](eos_lj.py) with [eos_lj_module.py](eos_lj_module.py)
* [eos_hs.py](eos_hs.py)
+An additional molecular dynamics program using quaternions
+[md_nvt_poly_lj.py](md_nvt_poly_lj.py) and [md_poly_lj_module.py](md_poly_lj_module.py)
+is described in the online [GUIDE](GUIDE.md),
+but does not appear in the text.
+
### 1.1 Calculation of T tensors
[t_tensor.py](t_tensor.py).
@@ -99,6 +104,9 @@ We also supply a constant-pressure version [mc_npt_hs.py](mc_npt_hs.py).
### 4.7 Monte Carlo program using quaternions
[mc_nvt_poly_lj.py](mc_nvt_poly_lj.py) and [mc_poly_lj_module.py](mc_poly_lj_module.py).
+We also supply a molecular dynamics program
+[md_nvt_poly_lj.py](md_nvt_poly_lj.py) and [md_poly_lj_module.py](md_poly_lj_module.py)
+to simulate the same model.
### 4.8 Monte Carlo of hard spherocylinders
[mc_nvt_sc.py](mc_nvt_sc.py) and [mc_sc_module.py](mc_sc_module.py).
diff --git a/python_examples/GUIDE.md b/python_examples/GUIDE.md
index db7a214..d103504 100644
--- a/python_examples/GUIDE.md
+++ b/python_examples/GUIDE.md
@@ -655,7 +655,8 @@ When comparing results with the MC program, several points should be remembered.
2. While we use _k_spring=10000 to highlight the multiple timestep method,
it is quite likely that energy flow between bond vibrations and other degrees of freedom will be inefficient,
due to the timescale separation.
-3. The constant-_NVE_ and constant-_NVT_ ensembles are expected to yield different behaviour around the collapse transition. We do not focus on this region, in the Python tests.
+3. The constant-_NVE_ and constant-_NVT_ ensembles are expected to yield different behaviour around the collapse transition.
+We do not focus on this region, in the Python tests.
4. Molecular dynamics is not expected to thoroughly explore the energy landscape at low temperatures,
giving instead (typically) quasi-harmonic vibrations in a single basin.
We do not focus on this region, in the Python tests.
@@ -771,6 +772,7 @@ _T_ | _E_ | _P_ | _T_ (K) | _PE_ (kJ/mol) | _PE_ (kJ/mol) eqn (23
A second set of tests was performed at _T_=0.6≈380K
at the specified densities ρ1, … ρ5 of Mossa et al (2002).
+A set of starting configurations is provided in the [Data repository](https://github.com/Allen-Tildesley/data).
Here the excess pressure (_P_(ex)=_P_-ρ_T_ converted to MPa
with a factor 77.75 based on the values of ε and σ)
is compared with the fit given by eqn (28) and the coefficients in Table III of Mossa et al (2002).
@@ -794,6 +796,74 @@ The MD simulations of Mossa et al (2002) are reported to extend to several hundr
(of order 107 MD timesteps) at the lowest temperatures.
It should be noted that this Python code is quite slow compared to the simpler atomic LJ examples.
+For comparison we provide a molecular dynamics code `md_nvt_poly_lj.py` for the same model.
+The program takes the molecular mass _M_ to be unity.
+Mossa et al (2002) ascribe a notional mass of 78u to each of the three LJ sites,
+so _M_≈3.9×10-25kg.
+Combined with the above values of ε and σ,
+this gives a time scale (_M_/ε)1/2σ ≈ 3.22 ps.
+The timestep of δt=0.01 ps used by Mossa et al (2002)
+corresponds to the default value in the program `dt=0.003` in these units.
+By default, the program simulates the constant-_NVE_ ensemble,
+but there is an option to simulate at constant _NVT_ by velocity randomization (Andersen thermostat).
+If the latter option is selected,
+the program will read configurations in the same format as `mc_nvt_poly_lj.py` (positions and quaternions only),
+selecting random initial velocities and angular momenta,
+which can be convenient.
+
+By default the program calculates the inertia tensor from the LJ site bond vectors,
+assuming equal masses.
+For simplicity it is assumed that the bond vectors are defined such that
+the principal axes of the inertia tensor coincide with
+the xyz axes of the molecular coordinate system,
+with the centre of mass at the origin;
+it is always possible to arrange this.
+In general,
+the three principal moments of inertia will all be different,
+so the molecule is an asymmetric top.
+The MD algorithm for rotation is a symplectic one
+in which a `kick` propagator advances the space-fixed angular momenta,
+using the torque on each molecule,
+and a succession of `drift` steps implement free rotation about each of the principal axes.
+This is described in the text, section 3.3; see
+
+* A Dullweber, B Leimkuhler, R McLachlan, _J Chem Phys,_ __107,__ 5840 (1997),
+* TF Miller, M Eleftheriou, P Pattnaik, A Ndirango, D Newns, GJ Martyna, _J Chem Phys,_ __116,__ 8649 (2002).
+
+The results below are for test runs in both constant-_NVE_ and constant-_NVT_ ensembles,
+at (approximately) the same state points as those given above.
+All runs were 10×5000 steps in length and used program defaults,
+except for `t_interval=1` and the specified temperature in the _NVT_ case.
+Because of the slow execution of the Python code,
+these runs are significantly shorter than the comparable Fortran examples,
+and very much shorter than the runs of Mossa et al (2002).
+For constant-_NVE_ runs we report RMS energy fluctuations,
+and _T_ is the average translational temperature.
+
+ ρ | _T_ | _E_ | _P_ | _E_(RMS)
+ ----- | ----- | ----- | ----- | -----
+ 0.32655 | 0.5 | -13.035(7) | 1.66(2) |
+ 0.32655 | 0.5025(9) | -13.0355 | 1.606(9) | 1.14×10-8
+ 0.32655 | 1.0 | -9.74(1) | 6.03(3) |
+ 0.32655 | 1.013(1) | -9.7372 | 5.975(7) | 1.01×10-7
+ 0.32655 | 1.5 | -6.81(1) | 9.41(4) |
+ 0.32655 | 1.508(1) | -6.8136 | 9.407(5) | 3.64×10-7
+ 0.32655 | 2.0 | -4.16(3) | 12.26(7) |
+ 0.32655 | 1.986(2) | -4.1570 | 12.291(8) | 9.04×10-7
+
+ ρ | _T_ | _E_ | _P_ | _E_ (RMS)
+ ----- | ----- | ----- | ----- | -----
+ 0.30533 | 0.6 | -11.624(6) | 0.48(2) |
+ 0.30533 | 0.6018(7) | -11.6239 | 0.453(6) | 1.33×10-8
+ 0.31240 | 0.6 | -11.902(9) | 1.13(2) |
+ 0.31240 | 0.604(1) | -11.9018 | 1.049(7) | 1.49×10-8
+ 0.31918 | 0.6 | -12.206(7) | 1.59(2) |
+ 0.31918 | 0.596(1) | -12.2065 | 1.63(1) | 1.64×10-8
+ 0.32655 | 0.6 | -12.379(8) | 2.50(2) |
+ 0.32655 | 0.599(1) | -12.3793 | 2.58(1) | 1.93×10-8
+ 0.33451 | 0.6 | -12.541(4) | 3.66(1) |
+ 0.33451 | 0.601(1) | -12.5411 | 3.70(1) | 2.40×10-8
+
## DPD program
For the `dpd.py` example, we recommend generating an initial configuration
using the `initialize.py` program, with JSON input similar to the following
diff --git a/python_examples/config_io_module.py b/python_examples/config_io_module.py
index be6b67f..274cf8f 100644
--- a/python_examples/config_io_module.py
+++ b/python_examples/config_io_module.py
@@ -64,7 +64,7 @@ def read_cnf_mols ( filename, with_v=False, quaternions=False ):
if with_v:
assert cols >= cols_re+6, "{:d}{}{:d}".format(cols,' cols less than',cols_re+6)
v = revw[:,cols_re :cols_re+3].astype(np.float_) # Velocity array
- w = revw[:,cols_re+3:cols_re+6].astype(np.float_) # Angular velocity array
+ w = revw[:,cols_re+3:cols_re+6].astype(np.float_) # Angular velocity/momentum array
return n, box, r, e, v, w
else:
return n, box, r, e
@@ -116,7 +116,7 @@ def write_cnf_mols ( filename, nn, box, *args ):
re=np.concatenate((args[0],args[1]),axis=1) # Just r and e
np.savetxt(filename,re,header=my_header,fmt='%15.10f',comments='')
- else: # Positions, orientations, velocities and angular velocities
+ else: # Positions, orientations, velocities and angular velocities/momenta
n, d = args[2].shape
assert n==nn, "v shape mismatch {:d}{:d}".format(n,nn)
assert d==3, "v shape mismatch {:d}".format(d)
diff --git a/python_examples/mc_nvt_poly_lj.py b/python_examples/mc_nvt_poly_lj.py
index c018cce..d712a43 100755
--- a/python_examples/mc_nvt_poly_lj.py
+++ b/python_examples/mc_nvt_poly_lj.py
@@ -140,11 +140,9 @@ def calc_variables ( ):
# Calculate all bond vectors
d = np.empty ( (n,na,3), dtype=np.float_ )
-di = np.empty ( (na,3), dtype=np.float_ )
for i, ei in enumerate(e):
ai = q_to_a ( ei ) # Rotation matrix for i
- for a in range(na):
- d[i,a,:] = np.dot ( db[a,:], ai ) # NB: equivalent to ai_T*db, ai_T=transpose of ai
+ d[i,:,:] = np.dot ( db, ai ) # NB: equivalent to ai_T*db, ai_T=transpose of ai
# Initial energy and overlap check
total = potential ( box, r, d )
@@ -172,8 +170,7 @@ def calc_variables ( ):
ri = ri - np.rint ( ri ) # Periodic boundary correction
ei = random_rotate_quaternion ( de_max, e[i,:] ) # Trial rotation
ai = q_to_a ( ei ) # Rotation matrix for i
- for a in range(na):
- di[a,:] = np.dot ( db[a,:], ai ) # NB: equivalent to ai_T*db, ai_T=transpose of ai
+ di = np.dot ( db, ai ) # NB: equivalent to ai_T*db, ai_T=transpose of ai
partial_new = potential_1 ( ri, di, box, rj, dj ) # New atom potential, virial etc
if not partial_new.ovr: # Test for non-overlapping configuration
diff --git a/python_examples/mc_poly_lj_module.py b/python_examples/mc_poly_lj_module.py
index 1abdb4e..f1aa48f 100644
--- a/python_examples/mc_poly_lj_module.py
+++ b/python_examples/mc_poly_lj_module.py
@@ -54,7 +54,7 @@ class PotentialType:
"""A composite variable for interactions."""
def __init__(self, pot, vir, ovr):
- self.pot = pot # the potential energy cut (but not shifted) at r_cut
+ self.pot = pot # the potential energy
self.vir = vir # the virial
self.ovr = ovr # a flag indicating overlap (i.e. pot too high to use)
diff --git a/python_examples/md_nvt_poly_lj.py b/python_examples/md_nvt_poly_lj.py
new file mode 100755
index 0000000..f4f2490
--- /dev/null
+++ b/python_examples/md_nvt_poly_lj.py
@@ -0,0 +1,305 @@
+#!/usr/bin/env python3
+# md_nvt_poly_lj.py
+
+#------------------------------------------------------------------------------------------------#
+# This software was written in 2016/17 #
+# by Michael P. Allen / #
+# and Dominic J. Tildesley ("the authors"), #
+# to accompany the book "Computer Simulation of Liquids", second edition, 2017 ("the text"), #
+# published by Oxford University Press ("the publishers"). #
+# #
+# LICENCE #
+# Creative Commons CC0 Public Domain Dedication. #
+# To the extent possible under law, the authors have dedicated all copyright and related #
+# and neighboring rights to this software to the PUBLIC domain worldwide. #
+# This software is distributed without any warranty. #
+# You should have received a copy of the CC0 Public Domain Dedication along with this software. #
+# If not, see . #
+# #
+# DISCLAIMER #
+# The authors and publishers make no warranties about the software, and disclaim liability #
+# for all uses of the software, to the fullest extent permitted by applicable law. #
+# The authors and publishers do not recommend use of this software for any purpose. #
+# It is made freely available, solely to clarify points made in the text. When using or citing #
+# the software, you should not imply endorsement by the authors or publishers. #
+#------------------------------------------------------------------------------------------------#
+
+"""Molecular dynamics, NVT ensemble (NVE option), polyatomic molecules."""
+
+def kick_propagator ( t ):
+ """Advances velocities and space-fixed angular momenta."""
+
+ # t is time over which to propagate, typically dt/2
+ # v and ell must be declared global to be updated in main program
+
+ global v, ell
+ import numpy as np
+
+ v = v + t * f # Linear momenta are equivalent to velocities
+ ell = ell + t * tau # Space-fixed angular momenta ell and torques tau
+
+def drift_translate ( t ):
+ """Advances positions."""
+
+ # t is time over which to propagate, typically dt
+ # r must be declared global to be updated in main program
+
+ global r
+ import numpy as np
+
+ r = r + t * v / box # Drift step (positions in box=1 units)
+ r = r - np.rint ( r ) # Periodic boundaries
+
+def drift_rotate ( xyz, t ):
+ """Advances quaternion orientations about a specified axis."""
+
+ # t is time over which to propagate, typically dt or dt/2
+ # e must be declared global to be updated in main program
+
+ global e
+ import numpy as np
+ from maths_module import rotate_quaternion
+
+ for i, ei in enumerate(e):
+ ai = q_to_a ( ei ) # Rotation matrix for i
+ w_hat = ai[xyz,:] # Space-fixed axis about which to rotate
+ w_mag = np.dot ( ell[i,:], w_hat ) / inertia[xyz] # Angular velocity about this axis
+ e[i,:] = rotate_quaternion ( w_mag*t, w_hat, ei ) # Rotate by specified angle
+
+def ran_velocities ( temperature, inertia ):
+ """Returns random velocities and space-fixed angular momenta."""
+
+ import numpy as np
+ from maths_module import q_to_a
+
+ v = np.random.randn ( n, 3 ) # Random velocities
+ v_cm = np.average ( v, axis=0 ) # Compute centre of mass velocity
+ v = v - v_cm # Set net momentum to zero
+ factor = np.sum ( v**2 ) / (3*n-3) # Estimate of kinetic temperature
+ factor = np.sqrt ( temperature / factor ) # Necessary correction factor
+ v = v * factor # Make correction
+
+ ell = np.random.randn ( n, 3 ) # Random body-fixed angular momenta
+ factors = np.sum ( ell**2, axis=0 ) / ( n*inertia ) # Estimate of kinetic temperatures
+ factors = np.sqrt ( temperature / factors ) # Necessary correction factors
+ ell = factors * ell # Make corrections
+
+ # Convert to space-fixed angular momenta
+ for i, ei in enumerate(e):
+ ai = q_to_a ( ei ) # Rotation matrix for i
+ ell[i,:] = np.dot ( ell[i,:], ai ) # NB: equivalent to ell_s = ai_T*ell_b, ai_T=transpose of ai
+
+ return v, ell
+
+def calc_variables ( ):
+ """Calculates all variables of interest.
+
+ They are collected and returned as a list, for use in the main program."""
+
+ import numpy as np
+ import math
+ from averages_module import msd, VariableType
+ from maths_module import q_to_a
+
+ # Preliminary calculations (n,r,v,total are taken from the calling program)
+ vol = box**3 # Volume
+ rho = n / vol # Density
+ kin_t = 0.5*np.sum(v**2) # Translational kinetic energy
+ kin_r = 0.0
+ for i, ei in enumerate(e):
+ ai = q_to_a ( ei ) # Rotation matrix for i
+ ell_i = np.dot ( ai, ell[i,:] ) # Get body-fixed angular momentum
+ kin_r = kin_r + np.sum((ell_i**2)/inertia) # Increment kinetic energy
+ kin_r = kin_r / 2 # Rotational kinetic energy
+
+ tmp_t = 2.0 * kin_t / (3*n-3) # Remove three degrees of freedom for momentum conservation
+ tmp_r = 2.0 * kin_r / (3*n) # 3N degrees of rotational freedom
+ eng = kin_t + kin_r + total.pot # Total energy for simulated system
+
+ # Variables of interest, of class VariableType, containing three attributes:
+ # .val: the instantaneous value
+ # .nam: used for headings
+ # .method: indicating averaging method
+ # If not set below, .method adopts its default value of avg
+ # The .nam and some other attributes need only be defined once, at the start of the program,
+ # but for clarity and readability we assign all the values together below
+
+ # Internal energy (shifted-force potential) per atom
+ # Total translational and rotational KE plus total PE divided by N
+ if nvt:
+ e_sf = VariableType ( nam = 'E/N shifted force', val = 3.0*temperature + total.pot/n )
+ else:
+ e_sf = VariableType ( nam = 'E/N shifted force', val = eng/n )
+
+ # Pressure (shifted-force potential)
+ # Ideal gas contribution plus total virial divided by V
+ if nvt:
+ p_sf = VariableType ( nam = 'P shifted force', val = rho*temperature + total.vir/vol )
+ else:
+ p_sf = VariableType ( nam = 'P shifted force', val = rho*tmp_t + total.vir/vol )
+
+ # Kinetic translational temperature
+ t_t = VariableType ( nam = 'T translational', val = tmp_t )
+
+ # Kinetic rotational temperature
+ t_r = VariableType ( nam = 'T rotational', val = tmp_r )
+
+ # Mean-squared deviation of conserved energy per atom
+ conserved_msd = VariableType ( nam = 'Conserved MSD', val = eng/n,
+ method = msd, e_format = True, instant = False )
+
+ # Collect together into a list for averaging
+ return [ e_sf, p_sf, t_t, t_r, conserved_msd ]
+
+# Takes in a configuration of atoms (positions, quaternions, velocities and angular momenta)
+# Cubic periodic boundary conditions
+# Conducts molecular dynamics with optional velocity thermalization
+# Uses no special neighbour lists
+# The rotational algorithm is described in the text, section 3.3. See A Dullweber, B Leimkuhler,
+# R McLachlan, J Chem Phys 107, 5840 (1997) and TF Miller, M Eleftheriou, P Pattnaik,
+# A Ndirango, D Newns, GJ Martyna, J Chem Phys 116, 8649 (2002).
+
+# Reads several variables and options from standard input using JSON format
+# Leave input empty "{}" to accept supplied defaults
+
+# Positions r are divided by box length after reading in and we assume mass=1 throughout
+# However, input configuration, output configuration, most calculations, and all results
+# are given in simulation units defined by the model
+# For example, for Lennard-Jones, sigma = 1, epsilon = 1
+
+# Despite the program name, there is nothing here specific to Lennard-Jones
+# The model is defined in md_poly_lj_module
+
+import json
+import sys
+import numpy as np
+import math
+from config_io_module import read_cnf_mols, write_cnf_mols
+from averages_module import run_begin, run_end, blk_begin, blk_end, blk_add
+from maths_module import q_to_a
+from md_poly_lj_module import introduction, conclusion, force, na, db, inertia, PotentialType
+
+cnf_prefix = 'cnf.'
+inp_tag = 'inp'
+out_tag = 'out'
+sav_tag = 'sav'
+
+print('md_nvt_poly_lj')
+print('Molecular dynamics, constant-NVT/NVE ensemble')
+print('Molecular mass=1 throughout')
+
+# Read parameters in JSON format
+try:
+ nml = json.load(sys.stdin)
+except json.JSONDecodeError:
+ print('Exiting on Invalid JSON format')
+ sys.exit()
+
+# Set default values, check keys and typecheck values
+defaults = {"nblock":10, "nstep":5000, "dt":0.003, "temperature":1.0, "t_interval":0 }
+for key, val in nml.items():
+ if key in defaults:
+ assert type(val) == type(defaults[key]), key+" has the wrong type"
+ else:
+ print('Warning', key, 'not in ',list(defaults.keys()))
+
+# Set parameters to input values or defaults
+nblock = nml["nblock"] if "nblock" in nml else defaults["nblock"]
+nstep = nml["nstep"] if "nstep" in nml else defaults["nstep"]
+dt = nml["dt"] if "dt" in nml else defaults["dt"]
+temperature = nml["temperature"] if "temperature" in nml else defaults["temperature"]
+t_interval = nml["t_interval"] if "t_interval" in nml else defaults["t_interval"]
+
+introduction()
+np.random.seed()
+
+# Write out parameters
+nvt = t_interval > 0 and t_interval < nstep
+print( "{:40}{:15d} ".format('Number of blocks', nblock) )
+print( "{:40}{:15d} ".format('Number of steps per block', nstep) )
+print( "{:40}{:15.6f}".format('Time step', dt) )
+if nvt:
+ print( "{:40} ".format('NVT ensemble' ) )
+ print( "{:40}{:15d} ".format('Thermalization interval', t_interval ) )
+ print( "{:40}{:15.6f}".format('Temperature', temperature ) )
+else:
+ print( "{:40} ".format('NVE ensemble' ) )
+ t_interval = nstep+1
+
+# Read in initial configuration
+if nvt:
+ n, box, r, e = read_cnf_mols ( cnf_prefix+inp_tag, quaternions=True )
+ v, ell = ran_velocities ( temperature, inertia )
+else:
+ n, box, r, e, v, ell = read_cnf_mols ( cnf_prefix+inp_tag, with_v=True, quaternions=True )
+ vcm = np.sum ( v, axis=0 ) / n # Centre-of mass velocity
+ v = v - vcm # Set COM velocity to zero
+
+print( "{:40}{:15d} ".format('Number of particles', n) )
+print( "{:40}{:15.6f}".format('Box length', box) )
+print( "{:40}{:15.6f}".format('Density', n/box**3) )
+r = r / box # Convert positions to box units
+r = r - np.rint ( r ) # Periodic boundaries
+
+# Calculate all bond vectors
+d = np.empty ( (n,na,3), dtype=np.float_ )
+norm = np.sqrt ( np.sum(e**2,axis=1) ) # Quaternion norms
+e = e / norm[:,np.newaxis] # Ensure normalized quaternions
+for i, ei in enumerate(e):
+ ai = q_to_a ( ei ) # Rotation matrix for i
+ d[i,:,:] = np.dot ( db, ai ) # NB: equivalent to ai_T*db, ai_T=transpose of ai
+
+# Initial forces, potential, etc plus overlap check
+total, f, tau = force ( box, r, d )
+assert not total.ovr, 'Overlap in initial configuration'
+
+# Initialize arrays for averaging and write column headings
+run_begin ( calc_variables() )
+
+for blk in range(1,nblock+1): # Loop over blocks
+
+ blk_begin()
+
+ for stp in range(nstep): # Loop over steps
+
+ if nvt and stp%t_interval == 0:
+ v, ell = ran_velocities ( temperature, inertia )
+
+ kick_propagator ( 0.5 * dt ) # Half-kick step
+
+ drift_translate ( dt ) # Drift step for positions
+
+ # Succession of drift steps for rotation about body-fixed axes
+ # Depending on the values of the moments of inertia, a different nested
+ # sequence of axes may produce better or worse energy conservation
+ drift_rotate ( 0, 0.5*dt )
+ drift_rotate ( 1, 0.5*dt )
+ drift_rotate ( 2, dt )
+ drift_rotate ( 1, 0.5*dt )
+ drift_rotate ( 0, 0.5*dt )
+
+ # Calculate all bond vectors
+ norm = np.sqrt ( np.sum(e**2,axis=1) ) # Quaternion norms
+ e = e / norm[:,np.newaxis] # Ensure normalized quaternions
+ for i, ei in enumerate(e):
+ ai = q_to_a ( ei ) # Rotation matrix for i
+ d[i,:,:] = np.dot ( db, ai ) # NB: equivalent to ai_T*db, ai_T=transpose of ai
+
+ total, f, tau = force ( box, r, d ) # Force evaluation
+ assert not total.ovr, 'Overlap in configuration'
+
+ kick_propagator ( 0.5 * dt ) # Half-kick step
+
+ blk_add ( calc_variables() )
+
+ blk_end(blk) # Output block averages
+ sav_tag = str(blk).zfill(3) if blk<1000 else 'sav' # Number configuration by block
+ write_cnf_mols ( cnf_prefix+sav_tag, n, box, r*box, e, v, ell ) # Save configuration
+
+run_end ( calc_variables() )
+
+total, f, tau = force ( box, r, d ) # Force evaluation
+assert not total.ovr, 'Overlap in final configuration'
+
+write_cnf_mols ( cnf_prefix+out_tag, n, box, r*box, e, v, ell ) # Save configuration
+conclusion()
diff --git a/python_examples/md_poly_lj_module.py b/python_examples/md_poly_lj_module.py
new file mode 100644
index 0000000..94837f0
--- /dev/null
+++ b/python_examples/md_poly_lj_module.py
@@ -0,0 +1,209 @@
+#!/usr/bin/env python3
+# md_poly_lj_module.py
+
+#------------------------------------------------------------------------------------------------#
+# This software was written in 2016/17 #
+# by Michael P. Allen / #
+# and Dominic J. Tildesley ("the authors"), #
+# to accompany the book "Computer Simulation of Liquids", second edition, 2017 ("the text"), #
+# published by Oxford University Press ("the publishers"). #
+# #
+# LICENCE #
+# Creative Commons CC0 Public Domain Dedication. #
+# To the extent possible under law, the authors have dedicated all copyright and related #
+# and neighboring rights to this software to the PUBLIC domain worldwide. #
+# This software is distributed without any warranty. #
+# You should have received a copy of the CC0 Public Domain Dedication along with this software. #
+# If not, see . #
+# #
+# DISCLAIMER #
+# The authors and publishers make no warranties about the software, and disclaim liability #
+# for all uses of the software, to the fullest extent permitted by applicable law. #
+# The authors and publishers do not recommend use of this software for any purpose. #
+# It is made freely available, solely to clarify points made in the text. When using or citing #
+# the software, you should not imply endorsement by the authors or publishers. #
+#------------------------------------------------------------------------------------------------#
+
+"""Force routine for MD simulation, polyatomic molecules, LJ atoms."""
+
+import numpy as np
+
+fast = True # Change this to replace NumPy force evaluation with slower Python
+
+# Bond vectors in body-fixed frame
+# Isosceles triangle, 3 sites, with unit bond length and bond angle alpha
+# which we set to 75 degrees here
+alpha = 75.0 * np.pi / 180.0
+alpha2 = alpha / 2.0
+na = 3
+db = np.array([[-np.sin(alpha2), 0.0, -np.cos(alpha2)/3.0],
+ [0.0, 0.0, 2.0*np.cos(alpha2)/3.0],
+ [np.sin(alpha2), 0.0, -np.cos(alpha2)/3.0]], dtype=np.float_)
+diameter = 2.0 * np.sqrt ( np.max ( np.sum(db**2,axis=1) ) ) # Molecular diameter
+
+# Cutoff distance and force-shift parameters (all private) chosen as per the reference:
+# S Mossa, E La Nave, HE Stanley, C Donati, F Sciortino, P Tartaglia, Phys Rev E, 65, 041205 (2002)
+r_cut = 2.612 # in sigma=1 units, where r_cut = 1.2616 nm, sigma = 0.483 nm
+sr_cut = 1.0/r_cut
+sr_cut6 = sr_cut**6
+sr_cut12 = sr_cut6**2
+lambda1 = 4.0*(7.0*sr_cut6-13.0*sr_cut12)
+lambda2 = -24.0*(sr_cut6-2.0*sr_cut12)*sr_cut
+
+m = np.array([ 1.0/3.0, 1.0/3.0, 1.0/3.0 ],dtype=np.float_) # Masses add up to 1.0
+
+# The following section sets the diagonal moments of inertia "realistically"
+# based on the values of atomic masses and bond vectors (above), with some checking
+# that the total mass is 1, the COM is at the origin, and the inertia tensor is diagonal.
+# However, there is nothing to stop the user replacing this section with a statement just
+# setting the values of inertia[:]. The masses m are not used by the calling program.
+# It might be advantageous, for instance, to artificially increase the values in inertia.
+
+# Ensure that the db bonds, xyz molecular axes, and masses are chosen such that
+# the total mass is 1 and the centre-of-mass is at the origin
+assert np.isclose(np.sum(m),1.0), 'M is not 1.0 {}'.format(np.sum(m))
+com = np.sum ( m[:,np.newaxis]*db, axis = 0 )
+assert np.all ( np.isclose(com,0.0) ), 'COM error {} {} {}'.format(*com)
+
+# Ensure that the db bonds, xyz molecular axes, and masses are chosen such that
+# the off-diagonal elements of the inertia tensor are zero
+inertia = -np.einsum('ij,ik->jk',m[:,np.newaxis]*db,db)
+offdiag = inertia[np.triu_indices(3,1)]
+assert np.all ( np.isclose(offdiag,0.0) ), 'Inertia not diagonal {} {} {}'.format(*offdiag)
+
+# Calculate the diagonal elements of the inertia tensor
+inertia = np.sum(m[:,np.newaxis]*db**2) + np.diagonal ( inertia )
+
+class PotentialType:
+ """A composite variable for interactions."""
+
+ def __init__(self, pot, vir, ovr):
+ self.pot = pot # the potential energy
+ self.vir = vir # the virial
+ self.ovr = ovr # a flag indicating overlap (i.e. pot too high to use)
+
+ def __add__(self, other):
+ pot = self.pot + other.pot
+ vir = self.vir + other.vir
+ ovr = self.ovr or other.ovr
+
+ return PotentialType(pot,vir,ovr)
+
+def introduction():
+ """Prints out introductory statements at start of run."""
+
+ print('Lennard-Jones potential')
+ print('Cut-and-force-shifted')
+ print('Diameter, sigma = 1')
+ print('Well depth, epsilon = 1')
+ if fast:
+ print('Fast NumPy force routine')
+ else:
+ print('Slow Python force routine')
+
+ print( "{:40}{:15d}".format('Number of atoms per molecule', na) )
+ for i, b in enumerate(db):
+ print( "{}{:2d}{:15.6f}{:15.6f}{:15.6f}".format('Body-fixed atom vector',i,*b))
+ print( "{:40}{:15.6f}".format('Molecular diameter', diameter) )
+
+ print( "{:40}{:15.6f}".format('r_cut', r_cut) )
+ print( "{:40}{:15.6f}".format('Force-shift lambda1', lambda1) )
+ print( "{:40}{:15.6f}".format('Force-shift lambda2', lambda2) )
+
+ print( "{:40}{:15.6f}{:15.6f}{:15.6f}".format('Inertia Ixx, Iyy, Izz', *inertia) )
+
+def conclusion():
+ """Prints out concluding statements at end of run."""
+
+ print('Program ends')
+
+def force ( box, r, d ):
+ """Takes in box, and r & d arrays, and calculates forces, torques and potentials etc."""
+
+ import numpy as np
+
+ # It is assumed that positions are in units where box = 1
+ # Forces are calculated in units where sigma = 1 and epsilon = 1
+ # Note that this is the force-shifted LJ potential with a linear smoothing term
+ # S Mossa, E La Nave, HE Stanley, C Donati, F Sciortino, P Tartaglia, Phys Rev E, 65, 041205 (2002)
+
+ n, ndim = r.shape
+ assert ndim==3, 'Dimension error for r'
+ nn, nna, ndim = d.shape
+ assert nna==na and ndim==3, 'Dimension error for d'
+ assert n==nn, 'Dimension mismatch for r and d'
+
+ sr2_ovr = 1.77 # Overlap threshold (pot > 100)
+ rm_cut_box = ( r_cut + diameter ) / box # Molecular cutoff in box=1 units
+ rm_cut_box_sq = rm_cut_box**2 # squared
+ assert rm_cut_box<0.5, 'rm_cut/box too large'
+ r_cut_sq = r_cut ** 2
+
+ # Initialize
+ f = np.zeros_like(r)
+ tau = np.zeros_like(r)
+ total = PotentialType ( pot=0.0, vir=0.0, ovr=False )
+
+ if fast:
+ for i in range(n-1):
+ rij = r[i,:]-r[i+1:,:] # Separation vectors for j>i
+ rij = rij - np.rint(rij) # Periodic boundary conditions in box=1 units
+ rij = rij * box # Now in sigma=1 units
+ for a in range(na):
+ for b in range(na):
+ rab = rij + d[i,a,:] - d[i+1:,b,:] # All atom-atom vectors for given a and b
+ rab_sq = np.sum(rab**2,axis=1) # Squared separations
+ in_range = rab_sq < r_cut_sq # Set flags for within cutoff
+ sr2 = 1.0 / rab_sq # (sigma/rab)**2
+ ovr = sr2 > sr2_ovr # Set flags for any overlaps
+ rmag = np.sqrt(rab_sq)
+ sr6 = sr2 ** 3
+ sr12 = sr6 ** 2
+ pot = np.where ( in_range,
+ 4.0*(sr12-sr6) + lambda1 + lambda2*rmag, 0.0 ) # force-shifted pair potentials
+ virab = np.where ( in_range,
+ 24.0*(2.0*sr12-sr6) - lambda2*rmag, 0.0 ) # pair virials
+ fab = virab * sr2
+ fab = rab * fab[:,np.newaxis] # atom-atom pair forces
+
+ total = total + PotentialType ( pot=np.sum(pot), vir=np.sum(rij*fab), ovr=np.any(ovr) )
+ fia = np.sum(fab,axis=0)
+ f[i,:] = f[i,:] + fia
+ f[i+1:,:] = f[i+1:,:] - fab
+ tau[i,:] = tau[i,:] + np.cross ( d[i,a,:], fia )
+ tau[i+1:,:] = tau[i+1:,:] - np.cross ( d[i+1:,b,:], fab )
+
+ else:
+ for i in range(n-1): # Outer loop
+ for j in range(i+1,n): # Inner loop
+ rij = r[i,:]-r[j,:] # Separation vector
+ rij = rij - np.rint(rij) # Periodic boundary conditions in box=1 units
+ rij_sq = np.sum(rij**2) # Squared separation
+
+ if rij_sq < rm_cut_box_sq: # Check within cutoff
+ rij = rij * box # Now in sigma=1 units
+ for a in range(na):
+ for b in range(na):
+ rab = rij + d[i,a,:] - d[j,b,:] # Atom-atom vector for given a and b
+ rab_sq = np.sum(rab**2) # Squared separation
+ if rab_sq < r_cut_sq: # Test within potential cutoff
+ sr2 = 1.0 / rab_sq # (sigma/rab)**2
+ ovr = sr2 > sr2_ovr # Set flag for overlap
+ rmag = np.sqrt(rab_sq)
+ sr6 = sr2 ** 3
+ sr12 = sr6 ** 2
+ pot = 4.0*(sr12-sr6) + lambda1 + lambda2*rmag # force-shifted pair potential
+ virab = 24.0*(2.0*sr12-sr6) - lambda2*rmag # pair virial
+ fab = virab * sr2
+ fab = rab * fab # atom-atom pair force
+
+ total = total + PotentialType ( pot=pot, vir=np.sum(rij*fab), ovr=ovr )
+ f[i,:] = f[i,:] + fab
+ f[j,:] = f[j,:] - fab
+ tau[i,:] = tau[i,:] + np.cross ( d[i,a,:], fab )
+ tau[j,:] = tau[j,:] - np.cross ( d[j,b,:], fab )
+
+ # Multiply results by numerical factors
+ total.vir = total.vir / 3.0 # Divide virial by 3
+
+ return total, f, tau