Skip to content

Commit

Permalink
Revised mc_nvt_poly_lj example as discussed in #8
Browse files Browse the repository at this point in the history
  • Loading branch information
@Michael-P-Allen
Michael-P-Allen committed Sep 6, 2017
1 parent 7fe333f commit 74ea172
Show file tree
Hide file tree
Showing 3 changed files with 137 additions and 134 deletions.
72 changes: 40 additions & 32 deletions GUIDE.md
View file
Original file line number Diff line number Diff line change
Expand Up @@ -789,7 +789,7 @@ composed of Lennard-Jones interaction sites.
For simplicity the sites are taken to be identical, although the program is easily generalized.
Molecular orientations are represented by quaternions,
which are used to calculate the rotation matrix
and hence the LJ site positions.
and hence the interaction site positions.

We test this with the three-site model of orthoterphenyl, a fragile glassformer,
described in the following publications amongst others.
Expand All @@ -799,46 +799,54 @@ described in the following publications amongst others.
* S Mossa, E La Nave, HE Stanley, C Donati, F Sciortino, P Tartaglia, _Phys Rev E,_ __65,__ 041205 (2002)
* E La Nave, S Mossa, F Sciortino, P Tartaglia, _J Chem Phys,_ __120,__ 6128 (2004)

We compare with the results of Mossa et al (2002).
The sites are arranged at the vertices of an isosceles triangle with bond angle 75 degrees,
LJ parameters &epsilon;/_k_<sub>B</sub> = 600K
or &epsilon; = 5 kJ mol<sup>-1</sup> to a good approximation,
LJ parameters &epsilon; = 5.276 kJ mol<sup>-1</sup>,
&sigma;=0.483nm,
and bond lengths equal to &sigma;.
and two equal bonds of length &sigma;.
The program employs the usual reduced units based on &epsilon; and &sigma;
and in these units the potential cutoff of Mossa et al (2002) is _R_<sub>c</sub>=2.612;
the pair potential is Lennard-Jones with a shifted-force correction term, linear in _r_,
to make the potential and its derivative vanish at _r_=_R_<sub>c</sub>.
Apart from the temperatures,
default program parameters were used throughout the tests.

Tests were performed at &rho;=0.32655 which is equivalent to &rho;<sub>4</sub>=1.108g cm<sup>-3</sup>
in Mossa et al (2002).
Comparisons of potential energy (_PE_=_E_-3 _T_ converted to kJ/mol with a factor 5)
were made with the fit given by eqn (6) of that paper.

_T_ | _E_ | _P_ | _T_ (K) | _PE_ (kJ/mol) | eqn (6)
----- | ----- | ----- | ----- | ----- | -----
0.5 | -14.30(1) | 1.63(3) | 300 | -79.00(5) | -77.52
1.0 | -11.21(1) | 5.52(3) | 600 | -71.05(5) | -68.55
1.5 | -8.262(7) | 8.95(2) | 900 | -63.81(4) | -61.19
2.0 | -5.52(1) | 11.86(3) | 1200 | -57.60(5) | -54.69

Exact agreement is not expected because the potential of Mossa et al (2002) has a different
cutoff correction, but the agreement is reasonable.
A second set of tests were performed at _T_=0.63333=380K
at the specified densities &rho;<sub>1</sub>, &hellip; &rho;<sub>5</sub>.
Comparisons of potential energy (_PE_=_E_-3 _T_ converted to kJ/mol with a factor 5.276)
were made with the fit given by eqn (23) of that paper.
Note that &epsilon;/k<sub>B</sub>&asymp;635 K.

_T_ | _E_ | _P_ | _T_ (K) | _PE_ (kJ/mol) | _PE_ (kJ/mol) eqn (23)
----- | ----- | ----- | ----- | ----- | -----
0.5 | -12.993(3) | 1.773(8) | 317 | -76.47(2) | -76.945
1.0 | -9.817(9) | 5.86(2) | 635 | -67.62(5) | -67.634
1.5 | -6.84(1) | 9.38(4) | 952 | -59.83(5) | -60.011
2.0 | -4.07(1) | 12.37(4) | 1270 | -53.13(5) | -53.265

A second set of tests was performed at _T_=0.6&asymp;380K
at the specified densities &rho;<sub>1</sub>, &hellip; &rho;<sub>5</sub> of Mossa et al (2002).
Here the excess pressure (_P_(ex)=_P_-&rho;_T_ converted to MPa
with a factor 73.54 based on the values of &epsilon; and &sigma;)
with a factor 77.75 based on the values of &epsilon; and &sigma;)
is compared with the fit given by eqn (28) and the coefficients in Table III of Mossa et al (2002).
NB the volumes to insert into the equation are those of their Table I,
which are specific to their system size.

Id | &rho; | _E_ | _P_ | _P_(ex) (MPa) | eqn (28)
----- | ----- | ----- | ----- | ----- | -----
1 | 0.30533 | -12.737(7) | 0.35(2) | 12(2) | 19.077
2 | 0.31240 | -13.025(7) | 0.99(2) | 58(2) | 60.143
3 | 0.31918 | -13.293(8) | 1.66(2) | 107(2) | 112.798
4 | 0.32655 | -13.50(1) | 2.60(2) | 176(1) | 177.222
5 | 0.33451 | -13.65(1) | 3.95(2) | 274(1) | 253.510

Although not perfect at the ends of the range, the agreement is not bad;
once more, the difference in cutoff correction should be borne in mind.
Note that Mossa et al (2002) used a different value
&epsilon;=5.276 kJ/mol which, with their cutoff term, gives a well depth 4.985 kJ/mol.
In addition their eqn (29) with coefficients in Table V is a fit to their potential energy,
which we calculate from the simulation as described above.

Id | &rho; | _E_ | _P_ | _P_(ex) (MPa) | _P_(ex) (MPa) eqn (28) | _PE_ (kJ/mol) | _PE_ (kJ/mol) eqn (29)
----- | ----- | ----- | ----- | ----- | ----- | ------ | ------
1 | 0.30533 | -11.625(4) | 0.45(1) | 20.7(7) | 19.077 | -70.83(2) | -70.818
2 | 0.31240 | -11.914(6) | 1.01(2) | 64(2) | 60.143 | -72.36(3) | -72.289
3 | 0.31918 | -12.143(5) | 1.74(1) | 120(2) | 112.798 | -73.56(3) | -73.601
4 | 0.32655 | -12.400(4) | 2.53(1) | 181(1) | 177.222 | -74.92(2) | -74.886
5 | 0.33451 | -12.487(4) | 3.84(1) | 283(1) | 253.510 | -75.38(2) | -75.825

In making these comparisons,
our default run length (10 blocks of 20000 sweeps each) should be borne in mind,
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 10<sup>7</sup> MD timesteps) at the lowest temperatures.

## DPD program
For the `dpd` example, we recommend generating an initial configuration
Expand Down
72 changes: 41 additions & 31 deletions mc_nvt_poly_lj.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -39,14 +39,14 @@ PROGRAM mc_nvt_poly_lj
! 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 mc_module
! The model (including the cutoff distance) is defined in mc_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 : metropolis, random_rotate_quaternion, random_translate_vector
USE maths_module, ONLY : metropolis, random_rotate_quaternion, random_translate_vector, q_to_a
USE mc_module, ONLY : introduction, conclusion, allocate_arrays, deallocate_arrays, &
& potential_1, potential, n, r, e, potential_type
& potential_1, potential, n, na, db, r, e, d, potential_type

IMPLICIT NONE

Expand All @@ -55,35 +55,34 @@ PROGRAM mc_nvt_poly_lj
REAL :: dr_max ! Maximum MC translational displacement
REAL :: de_max ! Maximum MC rotational displacement
REAL :: temperature ! Specified temperature
REAL :: r_cut ! Potential cutoff distance

! Composite interaction = pot & vir & ovr variables
TYPE(potential_type) :: total, partial_old, partial_new

INTEGER :: blk, stp, i, nstep, nblock, moves, ioerr
REAL :: delta, m_ratio
REAL, DIMENSION(3) :: ri
REAL, DIMENSION(0:3) :: ei
INTEGER :: blk, stp, i, a, nstep, nblock, moves, ioerr
REAL :: delta, m_ratio
REAL, DIMENSION(3) :: ri
REAL, DIMENSION(0:3) :: ei
REAL, DIMENSION(3,3) :: ai
REAL, DIMENSION(na,3) :: di

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, temperature, r_cut, dr_max, de_max
NAMELIST /nml/ nblock, nstep, temperature, dr_max, de_max

WRITE ( unit=output_unit, fmt='(a)' ) 'mc_nvt_poly_lj'
WRITE ( unit=output_unit, fmt='(a)' ) 'Monte Carlo, constant-NVT ensemble, polyatomic molecule'
WRITE ( unit=output_unit, fmt='(a)' ) 'Simulation uses cut-and-shifted potential'
CALL introduction

CALL RANDOM_SEED () ! Initialize random number generator

! Set sensible default run parameters for testing
nblock = 10
nstep = 10000
nstep = 20000
temperature = 1.0
r_cut = 2.5
dr_max = 0.05
de_max = 0.05

Expand All @@ -101,7 +100,6 @@ PROGRAM mc_nvt_poly_lj
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)' ) 'Temperature', temperature
WRITE ( unit=output_unit, fmt='(a,t40,f15.6)' ) 'Potential cutoff distance', r_cut
WRITE ( unit=output_unit, fmt='(a,t40,f15.6)' ) 'Maximum r displacement', dr_max
WRITE ( unit=output_unit, fmt='(a,t40,f15.6)' ) 'Maximum e displacement', de_max

Expand All @@ -110,13 +108,21 @@ PROGRAM mc_nvt_poly_lj
WRITE ( unit=output_unit, fmt='(a,t40,i15)' ) 'Number of molecules', n
WRITE ( unit=output_unit, fmt='(a,t40,f15.6)' ) 'Simulation box length', box
WRITE ( unit=output_unit, fmt='(a,t40,f15.6)' ) 'Density of molecules', REAL(n) / box**3
CALL allocate_arrays ( box, r_cut )
CALL allocate_arrays ( box )
CALL read_cnf_mols ( cnf_prefix//inp_tag, n, box, r, e ) ! Second call is to get r and e
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
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 energy and overlap check
total = potential ( box, r_cut )
total = potential ( box )
IF ( total%ovr ) THEN
WRITE ( unit=error_unit, fmt='(a)') 'Overlap in initial configuration'
STOP 'Error in mc_nvt_poly_lj'
Expand All @@ -136,7 +142,7 @@ PROGRAM mc_nvt_poly_lj

DO i = 1, n ! Begin loop over atoms

partial_old = potential_1 ( r(:,i), e(:,i), i, box, r_cut ) ! Old molecule potential, virial, etc
partial_old = potential_1 ( r(:,i), d(:,:,i), i, box ) ! Old molecule potential, virial, etc

IF ( partial_old%ovr ) THEN ! should never happen
WRITE ( unit=error_unit, fmt='(a)') 'Overlap in current configuration'
Expand All @@ -146,19 +152,24 @@ PROGRAM mc_nvt_poly_lj
ri = random_translate_vector ( dr_max/box, r(:,i) ) ! Trial move to new position (in box=1 units)
ri = ri - ANINT ( ri ) ! Periodic boundary correction
ei = random_rotate_quaternion ( de_max, e(:,i) ) ! Trial rotation
ai = q_to_a ( ei ) ! Rotation matrix for i
DO a = 1, na ! Loop over all atoms
di(:,a) = MATMUL ( db(:,a), ai ) ! NB: equivalent to ai_T*db, ai_T=transpose of ai
END DO ! End loop over all atoms

partial_new = potential_1 ( ri, ei, i, box, r_cut ) ! New molecule potential, virial etc
partial_new = potential_1 ( ri, di, i, box ) ! New molecule potential, virial etc

IF ( .NOT. partial_new%ovr ) THEN ! Test for non-overlapping configuration

delta = partial_new%pot - partial_old%pot ! Use cut-and-shifted potential
delta = delta / temperature ! Divide by temperature

IF ( metropolis ( delta ) ) THEN ! Accept Metropolis test
total = total + partial_new - partial_old ! Update potential energy
r(:,i) = ri ! Update position
e(:,i) = ei ! Update quaternion
moves = moves + 1 ! Increment move counter
total = total + partial_new - partial_old ! Update potential energy
r(:,i) = ri ! Update position
e(:,i) = ei ! Update quaternion
d(:,:,i) = di ! Update bond vectors
moves = moves + 1 ! Increment move counter
END IF ! End accept Metropolis test

END IF ! End test for non-overlapping configuration
Expand Down Expand Up @@ -195,12 +206,11 @@ FUNCTION calc_variables ( ) RESULT ( variables )
! 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

! In this example we simulate using the cut-and-shifted potential only
! The values of < p_s >, < e_s > and density should be consistent (for this potential)
! In this example we simulate using the shifted-force potential only
! The values of < p_sf >, < e_sf > and density should be consistent (for this potential)
! There are no long-range or delta corrections
! The value of the cut (but not shifted) potential is not used, in this example

TYPE(variable_type) :: m_r, e_s, p_s
TYPE(variable_type) :: m_r, e_sf, p_sf
REAL :: vol, rho

! Preliminary calculations
Expand All @@ -218,16 +228,16 @@ FUNCTION calc_variables ( ) RESULT ( variables )
! Move acceptance ratio
m_r = variable_type ( nam = 'Move ratio', val = m_ratio, instant = .FALSE. )

! Internal energy per molecule (cut-and-shifted potential)
! Ideal gas contribution (assuming nonlinear molecules) plus total (cut-and-shifted) PE divided by N
e_s = variable_type ( nam = 'E/N cut&shift', val = 3.0*temperature + total%pot/REAL(n) )
! Internal energy per molecule (shifted-force potential)
! Ideal gas contribution (assuming nonlinear molecules) plus total PE divided by N
e_sf = variable_type ( nam = 'E/N shifted force', val = 3.0*temperature + total%pot/REAL(n) )

! Pressure (cut-and-shifted potential)
! Pressure (shifted-force potential)
! Ideal gas contribution plus total virial divided by V
p_s = variable_type ( nam = 'P cut&shift', val = rho*temperature + total%vir/vol )
p_sf = variable_type ( nam = 'P shifted force', val = rho*temperature + total%vir/vol )

! Collect together for averaging
variables = [ m_r, e_s, p_s ]
variables = [ m_r, e_sf, p_sf ]

END FUNCTION calc_variables

Expand Down
Loading

0 comments on commit 74ea172

Please sign in to comment.