forked from Allen-Tildesley/examples
-
Notifications
You must be signed in to change notification settings - Fork 0
/
md_nvt_poly_lj.f90
362 lines (288 loc) · 16.3 KB
/
md_nvt_poly_lj.f90
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
! 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 <m.p.allen@warwick.ac.uk>/<m.p.allen@bristol.ac.uk> !
! and Dominic J. Tildesley <d.tildesley7@gmail.com> ("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 <http://creativecommons.org/publicdomain/zero/1.0/>. !
! !
! 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