1
#!/usr/bin/env python3
2
3
# Molecular dynamics in 1D
4
5
# import needed modules
6
import
random
7
import
math
8
import
sys
9
10
# the main loop function
11
def
main
(
md
):
12
time
=
0.0
# initialise time
13
14
# open coordinate output files
15
cfile
=
open
(
"coords.xyz"
,
"w"
)
16
# open temperature output file
17
tfile
=
open
(
"temperature.dat"
,
"w"
)
18
tfile
.
write
(
"# time temperature\n"
)
19
# open energy output file
20
efile
=
open
(
"energy.dat"
,
"w"
)
21
efile
.
write
(
"# time energy\n"
)
22
23
# main MD loop
24
for
t
in
range
(
md
.
tsteps
):
25
(
"#---- Time = "
,
round
(
time
,
2
),
" ---- Steps = "
,
t
,
"
----"
)
# python3
26
en
=
md
.
force
()
# calculate forces
27
md
.
integrate
(
t
,
en
)
# integrate equations of motion
28
# print current coordinates to file
29
md
.
printcoords
(
time
,
cfile
)
30
# print current temperature to file
31
tfile
.
write
(
str
(
round
(
time
,
2
))+
" "
+
str
(
md
.
temp
)+
"\n"
)
32
# print current energy to file
33
efile
.
write
(
str
(
round
(
time
,
2
))+
" "
+
str
(
md
.
etot
)+
"\n"
)
34
35
time
+=
md
.
dt
# increase time by dt
36
37
md
.
statistics
(
tfile
,
efile
)
# calculate averages, SD, etc
38
39
# close output files
40
cfile
.
close
()
41
tfile
.
close
()
42
efile
.
close
()
43
44
class
MD
(
object
):
45
46
N
=
36
# number of particles (integer for loop control)
47
dN
=
36.0
#number of particles (double for doing maths)
48
L
=
36.0
# length of 1D box
49
dim
=
1.0
# dimensions
50
# initialise positions and velocites
51
def
__init__
(
self
):
52
53
# declare the global variables
54
# constants
55
self
.
a
=
self
.
L
/
self
.
dN
# lattice spacing
56
self
.
dtsteps
=
100.0
# number of time steps (double for math)
57
self
.
tsteps
=
int
(
self
.
dtsteps
)
# number of time steps
(integer for counting)
58
self
.
dt
=
0.01
# integration timestep
59
self
.
rc
=
18.0
# distance cutoff for computing LJ interactions