Browse Source

Initialise Ångström-scale chemically powered sphere-dimer motor

Peter Colberg 3 years ago
commit
bea321289a
100 changed files with 9182 additions and 0 deletions
  1. 5
    0
      .gitignore
  2. 19
    0
      LICENSE
  3. 157
    0
      doc/angstrom-dimer.bib
  4. 259
    0
      examples/confined_dimer/confined_dimer.sh
  5. 64
    0
      examples/confined_dimer/equilibration/config.lua
  6. 152
    0
      examples/confined_dimer/equilibration/confined_dimer.lua
  7. 74
    0
      examples/confined_dimer/production/config.lua
  8. 165
    0
      examples/confined_dimer/production/confined_dimer.lua
  9. 69
    0
      examples/confined_dimer/production/msdd.lua
  10. 65
    0
      examples/confined_dimer/production/orient.lua
  11. 60
    0
      examples/confined_dimer/production/vafd.lua
  12. 126
    0
      examples/confined_dimer/production/walldf.lua
  13. 60
    0
      examples/confined_dimer/thermostat/config.lua
  14. 139
    0
      examples/confined_dimer/thermostat/confined_dimer.lua
  15. 69
    0
      examples/confined_dimer/tools/dimer_mean_square_displacement.py
  16. 62
    0
      examples/confined_dimer/tools/dimer_orientation_autocorrelation.py
  17. 65
    0
      examples/confined_dimer/tools/dimer_velocity_autocorrelation.py
  18. 52
    0
      examples/confined_dimer/tools/wall_density.py
  19. 1
    0
      examples/env.sh
  20. 64
    0
      examples/many_dimer/equilibration/config.lua
  21. 152
    0
      examples/many_dimer/equilibration/many_dimer.lua
  22. 257
    0
      examples/many_dimer/many_dimer.sh
  23. 72
    0
      examples/many_dimer/production/config.lua
  24. 163
    0
      examples/many_dimer/production/many_dimer.lua
  25. 69
    0
      examples/many_dimer/production/msdd.lua
  26. 65
    0
      examples/many_dimer/production/orient.lua
  27. 60
    0
      examples/many_dimer/production/vafd.lua
  28. 60
    0
      examples/many_dimer/thermostat/config.lua
  29. 158
    0
      examples/many_dimer/thermostat/many_dimer.lua
  30. 171
    0
      examples/many_dimer/tools/dimer_trajectory.py
  31. 62
    0
      examples/single_dimer/equilibration/config.lua
  32. 152
    0
      examples/single_dimer/equilibration/single_dimer.lua
  33. 86
    0
      examples/single_dimer/production/config.lua
  34. 136
    0
      examples/single_dimer/production/cylindf.lua
  35. 158
    0
      examples/single_dimer/production/cylinvf.lua
  36. 113
    0
      examples/single_dimer/production/forcepdf.lua
  37. 67
    0
      examples/single_dimer/production/msdd.lua
  38. 52
    0
      examples/single_dimer/production/msds.cl
  39. 88
    0
      examples/single_dimer/production/msds.lua
  40. 66
    0
      examples/single_dimer/production/orient.lua
  41. 145
    0
      examples/single_dimer/production/polardf.lua
  42. 138
    0
      examples/single_dimer/production/rdf.lua
  43. 179
    0
      examples/single_dimer/production/single_dimer.lua
  44. 61
    0
      examples/single_dimer/production/vafd.lua
  45. 49
    0
      examples/single_dimer/production/vafs.cl
  46. 85
    0
      examples/single_dimer/production/vafs.lua
  47. 137
    0
      examples/single_dimer/production/velpdf.lua
  48. 267
    0
      examples/single_dimer/single_dimer.sh
  49. 58
    0
      examples/single_dimer/thermostat/config.lua
  50. 139
    0
      examples/single_dimer/thermostat/single_dimer.lua
  51. 73
    0
      examples/single_dimer/tools/cylindrical_density.py
  52. 88
    0
      examples/single_dimer/tools/cylindrical_velocity.py
  53. 49
    0
      examples/single_dimer/tools/dimer_force.py
  54. 62
    0
      examples/single_dimer/tools/dimer_force_distribution.py
  55. 66
    0
      examples/single_dimer/tools/dimer_mean_square_displacement.py
  56. 59
    0
      examples/single_dimer/tools/dimer_orientation_autocorrelation.py
  57. 127
    0
      examples/single_dimer/tools/dimer_solvent.py
  58. 171
    0
      examples/single_dimer/tools/dimer_trajectory.py
  59. 50
    0
      examples/single_dimer/tools/dimer_velocity.py
  60. 62
    0
      examples/single_dimer/tools/dimer_velocity_autocorrelation.py
  61. 63
    0
      examples/single_dimer/tools/dimer_velocity_distribution.py
  62. 58
    0
      examples/single_dimer/tools/polar_density.py
  63. 56
    0
      examples/single_dimer/tools/radial_density.py
  64. 51
    0
      examples/single_dimer/tools/solvent_force_distribution.py
  65. 62
    0
      examples/single_dimer/tools/solvent_mean_square_displacement.py
  66. 62
    0
      examples/single_dimer/tools/solvent_velocity_autocorrelation.py
  67. 51
    0
      examples/single_dimer/tools/solvent_velocity_distribution.py
  68. 33
    0
      nanomotor/box.cl
  69. 41
    0
      nanomotor/box.lua
  70. 59
    0
      nanomotor/compute.lua
  71. 16
    0
      nanomotor/config.lua
  72. 173
    0
      nanomotor/correlation.lua
  73. 77
    0
      nanomotor/cumsum.cl
  74. 67
    0
      nanomotor/cumsum.lua
  75. 184
    0
      nanomotor/domain.lua
  76. 259
    0
      nanomotor/h5md.lua
  77. 76
    0
      nanomotor/hilbert.cl
  78. 103
    0
      nanomotor/integrate.lua
  79. 137
    0
      nanomotor/ljd.cl
  80. 71
    0
      nanomotor/ljd.lua
  81. 78
    0
      nanomotor/ljs.cl
  82. 56
    0
      nanomotor/ljs.lua
  83. 58
    0
      nanomotor/neighd.cl
  84. 62
    0
      nanomotor/neighd.lua
  85. 79
    0
      nanomotor/neighs.cl
  86. 101
    0
      nanomotor/neighs.lua
  87. 36
    0
      nanomotor/observe.lua
  88. 99
    0
      nanomotor/position.cl
  89. 101
    0
      nanomotor/position.lua
  90. 128
    0
      nanomotor/random.cl
  91. 63
    0
      nanomotor/random.lua
  92. 77
    0
      nanomotor/rattle.cl
  93. 90
    0
      nanomotor/rattle.lua
  94. 79
    0
      nanomotor/snapshot.lua
  95. 60
    0
      nanomotor/sort.cl
  96. 111
    0
      nanomotor/sort.lua
  97. 40
    0
      nanomotor/species.cl
  98. 96
    0
      nanomotor/species.lua
  99. 70
    0
      nanomotor/statistics.lua
  100. 0
    0
      nanomotor/thermo.lua

+ 5
- 0
.gitignore View File

@@ -0,0 +1,5 @@
1
+*.h5
2
+*.h5.tmp
3
+*.log
4
+*.pdf
5
+*.mp4

+ 19
- 0
LICENSE View File

@@ -0,0 +1,19 @@
1
+Copyright © 2013–2015 Peter Colberg.
2
+
3
+Permission is hereby granted, free of charge, to any person obtaining a copy
4
+of this software and associated documentation files (the "Software"), to deal
5
+in the Software without restriction, including without limitation the rights
6
+to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
7
+copies of the Software, and to permit persons to whom the Software is
8
+furnished to do so, subject to the following conditions:
9
+
10
+The above copyright notice and this permission notice shall be included in
11
+all copies or substantial portions of the Software.
12
+
13
+THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
14
+IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
15
+FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.  IN NO EVENT SHALL THE
16
+AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
17
+LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
18
+OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
19
+THE SOFTWARE.

+ 157
- 0
doc/angstrom-dimer.bib View File

@@ -0,0 +1,157 @@
1
+@article{Andersen1983,
2
+  author = {Hans C. Andersen},
3
+  title = {Rattle: A ``velocity'' version of the shake algorithm for molecular
4
+    dynamics calculations},
5
+  journal = {Journal of Computational Physics},
6
+  year = {1983},
7
+  volume = {52},
8
+  pages = {24--34},
9
+  number = {1},
10
+  doi = {10.1016/0021-9991(83)90014-1},
11
+}
12
+
13
+@techreport{Blelloch1990,
14
+  author = {Guy E. Blelloch},
15
+  title = {Prefix Sums and Their Applications},
16
+  institution = {School of Computer Science, Carnegie Mellon University},
17
+  year = {1990},
18
+  number = {CMU-CS-90-190},
19
+  url = {http://www.cs.cmu.edu/~scandal/papers/CMU-CS-90-190.html}
20
+}
21
+
22
+@article{Colberg2011,
23
+  author = {Peter H. Colberg and Felix H{\"o}f{}ling},
24
+  title = {Highly accelerated simulations of glassy dynamics using {GPU}s:
25
+    Caveats on limited floating-point precision},
26
+  journal = {Computer Physics Communications},
27
+  year = {2011},
28
+  volume = {182},
29
+  pages = {1120--1129},
30
+  number = {5},
31
+  doi = {10.1016/j.cpc.2011.01.009}
32
+}
33
+
34
+@article{Cook1957,
35
+  author = {J. M. Cook},
36
+  title = {Rational Formulae for the Production of a Spherically Symmetric
37
+    Probability Distribution},
38
+  journal = {Mathematics of Computation},
39
+  year = {1957},
40
+  volume = {11},
41
+  pages = {81--82},
42
+  doi = {10.1090/S0025-5718-1957-0690630-7}
43
+}
44
+
45
+@misc{Ferguson2010,
46
+  author = {Niels Ferguson and Stefan Lucks and Bruce Schneier and Doug Whiting
47
+    and Mihir Bellare and Tadayoshi Kohno and Jon Callas and Jesse Walker},
48
+  title = {The Skein Hash Function Family},
49
+  year = {2010},
50
+  url = {http://www.skein-hash.info/sites/default/files/skein1.3.pdf}
51
+}
52
+
53
+@techreport{Hamilton2006,
54
+  author = {Chris Hamilton},
55
+  title = {Compact Hilbert Indices},
56
+  institution = {Dalhousie University, Faculty of Computer Science},
57
+  year = {2006},
58
+  number = {Technical Report CS-2006-07},
59
+  url = {http://www.cs.dal.ca/research/techreports/cs-2006-07}
60
+}
61
+
62
+@article{Hamilton2008,
63
+  author = {Chris H. Hamilton and Andrew Rau-Chaplin},
64
+  title = {Compact Hilbert indices: Space-filling curves for domains with
65
+    unequal side lengths},
66
+  journal = {Information Processing Letters},
67
+  year = {2008},
68
+  volume = {105},
69
+  pages = {155--163},
70
+  number = {5},
71
+  doi = {10.1016/j.ipl.2007.08.034}
72
+}
73
+
74
+@article{Hilbert1891,
75
+  author = {Hilbert, David},
76
+  title = {Ueber die stetige Abbildung einer Linie auf ein Fl{\"a}chenst{\"u}ck},
77
+  journal = {Mathematische Annalen},
78
+  year = {1891},
79
+  volume = {38},
80
+  pages = {459-460},
81
+  doi = {10.1007/BF01199431},
82
+  issue = {3},
83
+  language = {German}
84
+}
85
+
86
+@article{Marsaglia1964,
87
+  author = {Marsaglia, G. and Bray, T.},
88
+  title = {A Convenient Method for Generating Normal Variables},
89
+  journal = {SIAM Review},
90
+  year = {1964},
91
+  volume = {6},
92
+  pages = {260-264},
93
+  number = {3},
94
+  doi = {10.1137/1006063}
95
+}
96
+
97
+@article{Marsaglia1972,
98
+  author = {George Marsaglia},
99
+  title = {Choosing a Point from the Surface of a Sphere},
100
+  journal = {Annals of Mathematical Statistics},
101
+  year = {1972},
102
+  volume = {43},
103
+  pages = {645--646},
104
+  number = {2},
105
+  doi = {10.1214/aoms/1177692644}
106
+}
107
+
108
+@article{Neumann1951,
109
+  author = {John von Neumann},
110
+  title = {Various techniques used in connection with random digits},
111
+  journal = {J. Res Nat. Bur. Stand., Applied Mathematics Series},
112
+  year = {1951},
113
+  volume = {12},
114
+  pages = {36--38}
115
+}
116
+
117
+@inproceedings{Salmon2011,
118
+  author = {Salmon, John K. and Moraes, Mark A. and Dror, Ron O. and
119
+    Shaw, David E.},
120
+  title = {Parallel random numbers: as easy as 1, 2, 3},
121
+  booktitle = {Proceedings of 2011 International Conference for High
122
+    Performance Computing, Networking, Storage and Analysis},
123
+  year = {2011},
124
+  series = {SC '11},
125
+  pages = {16:1--16:12},
126
+  address = {New York, NY, USA},
127
+  publisher = {ACM},
128
+  articleno = {16},
129
+  doi = {10.1145/2063384.2063405},
130
+  isbn = {978-1-4503-0771-0},
131
+  location = {Seattle, Washington}
132
+}
133
+
134
+@inproceedings{Sengupta2007,
135
+  author = {Shubhabrata Sengupta and Mark Harris and Yao Zhang and
136
+    John D. Owens},
137
+  title = {Scan Primitives for GPU Computing},
138
+  booktitle = {Graphics Hardware 2007},
139
+  year = {2007},
140
+  pages = {97--106},
141
+  organization = {ACM},
142
+  location = {San Diego, CA}
143
+}
144
+
145
+@article{Swope1982,
146
+  author = {William C. Swope and Hans C. Andersen and Peter H. Berens and
147
+    Kent R. Wilson},
148
+  title = {A computer simulation method for the calculation of equilibrium
149
+    constants for the formation of physical clusters of molecules: Application
150
+    to small water clusters},
151
+  journal = {Journal of Chemical Physics},
152
+  year = {1982},
153
+  volume = {76},
154
+  pages = {637--649},
155
+  number = {1},
156
+  doi = {10.1063/1.442716}
157
+}

+ 259
- 0
examples/confined_dimer/confined_dimer.sh View File

@@ -0,0 +1,259 @@
1
+#PBS -q gpu
2
+#PBS -l nodes=1:ppn=2:gpus=2
3
+#PBS -l walltime=3:00:00
4
+#PBS -m n
5
+#PBS -V
6
+set -e -o pipefail
7
+PROJECT="${PBS_O_WORKDIR}/angstrom-dimer"
8
+cd "${PROJECT}" && . examples/env.sh
9
+OUTPUT="${OUTPUT:-${PBS_JOBNAME%%.sh*}_${PBS_JOBID%%[^0-9]*}}"
10
+SEED="$(date +%s)\0${PBS_O_HOST}\0${PBS_JOBID}"
11
+TASKS=()
12
+for TASK in 1 2
13
+do
14
+  if [ ! -e "${PBS_O_WORKDIR}/thermostat/${OUTPUT}_${TASK}.h5" ]
15
+  then
16
+    mkdir -p "${PBS_O_WORKDIR}/thermostat"
17
+    cat >"${PBS_O_WORKDIR}/thermostat/${OUTPUT}_${TASK}.lua" <<EOF
18
+-- Edge lengths of simulation domain.
19
+L = {50, 50, 20}
20
+-- Periodic boundary conditions.
21
+periodic = {true, true, false}
22
+-- Diameters of dimer spheres.
23
+diameter = {C = 4.0, N = 8.0}
24
+-- Lennard-Jones wall potential.
25
+wall = {dimer = {epsilon = 5.0, r_cut = 3^(1/6), sigma = {C = (diameter.C + 1)/2, N = (diameter.N + 1)/2}}, solvent = {epsilon = 1.0, r_cut = 3^(1/6), sigma = 1.0}}
26
+-- Number of dimer spheres.
27
+Nd = 2
28
+-- Solvent number density.
29
+rho = 0.8
30
+-- Number of solvent particles excluding volume of walls.
31
+Ns = math.floor(rho * L[1] * L[2] * (L[3] - 1.0))
32
+-- Cutoff radius of potential in units of σ.
33
+r_cut = 2^(1/6)
34
+-- Minimum number of steps between solvent particle sorts.
35
+sort = {interval = 100}
36
+-- Lattice constant and region for solvent particles.
37
+lattice = {a = 1.6, slab = {{0, 0, 1}, {L[1], L[2], L[3]-1}}}
38
+-- Solvent temperature.
39
+temp = 1
40
+-- Masses of neutrally buoyant dimer spheres.
41
+mass = {C = rho * (math.pi/6) * diameter.C^3, N = rho * (math.pi/6) * diameter.N^3}
42
+-- Distance of dimer spheres.
43
+bond = ((diameter.C + diameter.N)/2 + 1) * r_cut
44
+-- Potential well depths of dimer spheres with solvent particles.
45
+epsilon = {CA = 1.0, CB = 1.0, NA = 1.0, NB = 0.1}
46
+-- Potential diameters.
47
+sigma = {C = (diameter.C + 1)/2, N = (diameter.N + 1)/2}
48
+-- Verlet neighbour list skin.
49
+skin = 0.5
50
+-- Number of placeholders per neighbour list.
51
+nlist_size = {ss = 25, ds = 500}
52
+-- Number of bins per dimension for solvent particles.
53
+nbin = {L[1]/2, L[2]/2, L[3]/2}
54
+-- Number of placeholders per bin.
55
+bin_size = 20
56
+-- Integration time-step.
57
+timestep = 0.001
58
+-- Number of steps.
59
+nsteps = 20000
60
+-- Thermodynamic variables.
61
+thermo = {interval = 1000, datatype = "double"}
62
+-- Dimer trajectory.
63
+trajd = {interval = 1000, datatype = "double", species = {"C", "N"}}
64
+-- Solvent trajectory.
65
+trajs = {interval = nsteps, datatype = "double", species = {"A", "B"}}
66
+-- OpenCL platforms.
67
+platform = 1
68
+-- OpenCL device per platform.
69
+device = ${TASK}
70
+-- Random number generator seed.
71
+seed = "0x$(echo -ne "${SEED}\0${TASK}\0thermostat" | md5sum | cut -c-16)"
72
+-- H5MD output filename.
73
+output = "${PBS_O_WORKDIR}/thermostat/${OUTPUT}_${TASK}.h5"
74
+-- H5MD file author.
75
+author = {name = "$(git config user.name)", email = "$(git config user.email)"}
76
+-- H5MD file creator.
77
+creator = {name = "confined_dimer.lua", version = "$(cd "${PROJECT}" && git describe --always)"}
78
+EOF
79
+    cd "${PROJECT}/examples/confined_dimer/thermostat"
80
+    ./confined_dimer.lua "${PBS_O_WORKDIR}/thermostat/${OUTPUT}_${TASK}.lua" | awk '{ print strftime("%Y-%m-%d %H:%M:%S"), $0; fflush(); }' >>"${PBS_O_WORKDIR}/thermostat/${OUTPUT}_${TASK}.log" &
81
+    TASKS+=("${!}")
82
+  fi
83
+done
84
+for TASK in "${TASKS[@]}"
85
+do
86
+  wait "${TASK}"
87
+done
88
+TASKS=()
89
+for TASK in 1 2
90
+do
91
+  if [ ! -e "${PBS_O_WORKDIR}/equilibration/${OUTPUT}_${TASK}.h5" ]
92
+  then
93
+    mkdir -p "${PBS_O_WORKDIR}/equilibration"
94
+    cat >"${PBS_O_WORKDIR}/equilibration/${OUTPUT}_${TASK}.lua" <<EOF
95
+-- Edge lengths of simulation domain.
96
+L = {50, 50, 20}
97
+-- Periodic boundary conditions.
98
+periodic = {true, true, false}
99
+-- Diameters of dimer spheres.
100
+diameter = {C = 4.0, N = 8.0}
101
+-- Lennard-Jones wall potential.
102
+wall = {dimer = {epsilon = 5.0, r_cut = 3^(1/6), sigma = {C = (diameter.C + 1)/2, N = (diameter.N + 1)/2}}, solvent = {epsilon = 1.0, r_cut = 3^(1/6), sigma = 1.0}}
103
+-- Number of dimer spheres.
104
+Nd = 2
105
+-- Solvent number density.
106
+rho = 0.8
107
+-- Number of solvent particles excluding volume of walls.
108
+Ns = math.floor(rho * L[1] * L[2] * (L[3] - 1.0))
109
+-- Cutoff radius of potential in units of σ.
110
+r_cut = 2^(1/6)
111
+-- Minimum number of steps between solvent particle sorts.
112
+sort = {interval = 100}
113
+-- Masses of neutrally buoyant dimer spheres.
114
+mass = {C = rho * (math.pi/6) * diameter.C^3, N = rho * (math.pi/6) * diameter.N^3}
115
+-- Distance of dimer spheres.
116
+bond = ((diameter.C + diameter.N)/2 + 1) * r_cut
117
+-- Potential well depths of dimer spheres with solvent particles.
118
+epsilon = {CA = 1.0, CB = 1.0, NA = 1.0, NB = 0.1}
119
+-- Potential diameters.
120
+sigma = {C = (diameter.C + 1)/2, N = (diameter.N + 1)/2}
121
+-- Forward and backward reaction of solvent particles.
122
+reaction = {rate = 0.001}
123
+-- Verlet neighbour list skin.
124
+skin = 0.5
125
+-- Number of placeholders per neighbour list.
126
+nlist_size = {ss = 25, ds = 550}
127
+-- Number of bins per dimension for solvent particles.
128
+nbin = {L[1]/2, L[2]/2, L[3]/2}
129
+-- Number of placeholders per bin.
130
+bin_size = 20
131
+-- Integration time-step.
132
+timestep = 0.001
133
+-- Number of steps.
134
+nsteps = 2000000
135
+-- Thermodynamic variables.
136
+thermo = {interval = 1000, datatype = "double"}
137
+-- Dimer trajectory.
138
+trajd = {interval = 1000, datatype = "double", species = {"C", "N"}}
139
+-- Solvent trajectory.
140
+trajs = {interval = nsteps, datatype = "double", species = {"A", "B"}}
141
+-- Count solvent species.
142
+species = {interval = 1000, species = {"A", "B"}}
143
+-- OpenCL platforms.
144
+platform = 1
145
+-- OpenCL device per platform.
146
+device = ${TASK}
147
+-- Random number generator seed.
148
+seed = "0x$(echo -ne "${SEED}\0${TASK}\0equilibration" | md5sum | cut -c-16)"
149
+-- H5MD input filename.
150
+input = "${PBS_O_WORKDIR}/thermostat/${OUTPUT}_${TASK}.h5"
151
+-- Read trajectory at step.
152
+step = 20000
153
+-- H5MD output filename.
154
+output = "${PBS_O_WORKDIR}/equilibration/${OUTPUT}_${TASK}.h5"
155
+-- H5MD file author.
156
+author = {name = "$(git config user.name)", email = "$(git config user.email)"}
157
+-- H5MD file creator.
158
+creator = {name = "confined_dimer.lua", version = "$(cd "${PROJECT}" && git describe --always)"}
159
+EOF
160
+    cd "${PROJECT}/examples/confined_dimer/equilibration"
161
+    ./confined_dimer.lua "${PBS_O_WORKDIR}/equilibration/${OUTPUT}_${TASK}.lua" | awk '{ print strftime("%Y-%m-%d %H:%M:%S"), $0; fflush(); }' >>"${PBS_O_WORKDIR}/equilibration/${OUTPUT}_${TASK}.log" &
162
+    TASKS+=("${!}")
163
+  fi
164
+done
165
+for TASK in "${TASKS[@]}"
166
+do
167
+  wait "${TASK}"
168
+done
169
+TASKS=()
170
+for TASK in 1 2
171
+do
172
+  if [ ! -e "${PBS_O_WORKDIR}/production/${OUTPUT}_${TASK}.h5" ]
173
+  then
174
+    mkdir -p "${PBS_O_WORKDIR}/production"
175
+    cat >"${PBS_O_WORKDIR}/production/${OUTPUT}_${TASK}.lua" <<EOF
176
+-- Edge lengths of simulation domain.
177
+L = {50, 50, 20}
178
+-- Periodic boundary conditions.
179
+periodic = {true, true, false}
180
+-- Diameters of dimer spheres.
181
+diameter = {C = 4.0, N = 8.0}
182
+-- Lennard-Jones wall potential.
183
+wall = {dimer = {epsilon = 5.0, r_cut = 3^(1/6), sigma = {C = (diameter.C + 1)/2, N = (diameter.N + 1)/2}}, solvent = {epsilon = 1.0, r_cut = 3^(1/6), sigma = 1.0}}
184
+-- Number of dimer spheres.
185
+Nd = 2
186
+-- Solvent number density.
187
+rho = 0.8
188
+-- Number of solvent particles excluding volume of walls.
189
+Ns = math.floor(rho * L[1] * L[2] * (L[3] - 1.0))
190
+-- Cutoff radius of potential in units of σ.
191
+r_cut = 2^(1/6)
192
+-- Minimum number of steps between solvent particle sorts.
193
+sort = {interval = 100}
194
+-- Masses of neutrally buoyant dimer spheres.
195
+mass = {C = rho * (math.pi/6) * diameter.C^3, N = rho * (math.pi/6) * diameter.N^3}
196
+-- Distance of dimer spheres.
197
+bond = ((diameter.C + diameter.N)/2 + 1) * r_cut
198
+-- Potential well depths of dimer spheres with solvent particles.
199
+epsilon = {CA = 1.0, CB = 1.0, NA = 1.0, NB = 0.1}
200
+-- Potential diameters.
201
+sigma = {C = (diameter.C + 1)/2, N = (diameter.N + 1)/2}
202
+-- Forward and backward reaction of solvent particles.
203
+reaction = {rate = 0.001}
204
+-- Verlet neighbour list skin.
205
+skin = 0.5
206
+-- Number of placeholders per neighbour list.
207
+nlist_size = {ss = 25, ds = 550}
208
+-- Number of bins per dimension for solvent particles.
209
+nbin = {L[1]/2, L[2]/2, L[3]/2}
210
+-- Number of placeholders per bin.
211
+bin_size = 20
212
+-- Integration time-step.
213
+timestep = 0.001
214
+-- Number of steps.
215
+nsteps = 10000000
216
+-- Thermodynamic variables.
217
+thermo = {interval = 10000, datatype = "double"}
218
+-- Dimer trajectory.
219
+trajd = {interval = 1000, datatype = "float", species = {"C", "N"}}
220
+-- Solvent trajectory.
221
+trajs = {interval = nsteps, datatype = "float", species = {"A", "B"}}
222
+-- Count solvent species.
223
+species = {interval = 10000, species = {"A", "B"}}
224
+-- Mean-square displacement of dimer.
225
+msdd = {interval = 10, blockcount = 6, blocksize = 100, blockshift = 10, separation = 10000}
226
+-- Velocity autocorrelation function of dimer.
227
+vafd = {interval = 10, blockcount = 6, blocksize = 100, blockshift = 10, separation = 10000}
228
+-- Orientation autocorrelation of dimer.
229
+orient = {interval = 10, blockcount = 6, blocksize = 100, blockshift = 10, separation = 10000}
230
+-- Wall density function.
231
+walldf = {interval = 1000, nbin = 1000, cutoff = (diameter.C + diameter.N)/2 + bond}
232
+-- OpenCL platforms.
233
+platform = 1
234
+-- OpenCL device per platform.
235
+device = ${TASK}
236
+-- Random number generator seed.
237
+seed = "0x$(echo -ne "${SEED}\0${TASK}\0production" | md5sum | cut -c-16)"
238
+-- H5MD input filename.
239
+input = "${PBS_O_WORKDIR}/equilibration/${OUTPUT}_${TASK}.h5"
240
+-- Read trajectory at step.
241
+step = 2000000
242
+-- H5MD output filename.
243
+output = "${PBS_O_WORKDIR}/production/${OUTPUT}_${TASK}.h5"
244
+-- Snapshot program state.
245
+snapshot = {interval = 1000000, output = "${PBS_O_WORKDIR}/production/${OUTPUT}_${TASK}-snapshot.h5"}
246
+-- H5MD file author.
247
+author = {name = "$(git config user.name)", email = "$(git config user.email)"}
248
+-- H5MD file creator.
249
+creator = {name = "confined_dimer.lua", version = "$(cd "${PROJECT}" && git describe --always)"}
250
+EOF
251
+    cd "${PROJECT}/examples/confined_dimer/production"
252
+    ./confined_dimer.lua "${PBS_O_WORKDIR}/production/${OUTPUT}_${TASK}.lua" | awk '{ print strftime("%Y-%m-%d %H:%M:%S"), $0; fflush(); }' >>"${PBS_O_WORKDIR}/production/${OUTPUT}_${TASK}.log" &
253
+    TASKS+=("${!}")
254
+  fi
255
+done
256
+for TASK in "${TASKS[@]}"
257
+do
258
+  wait "${TASK}"
259
+done

+ 64
- 0
examples/confined_dimer/equilibration/config.lua View File

@@ -0,0 +1,64 @@
1
+-- Edge lengths of simulation domain.
2
+L = {50, 50, 20}
3
+-- Periodic boundary conditions.
4
+periodic = {true, true, false}
5
+-- Diameters of dimer spheres.
6
+diameter = {C = 4.0, N = 8.0}
7
+-- Lennard-Jones wall potential.
8
+wall = {dimer = {epsilon = 5.0, r_cut = 3^(1/6), sigma = {C = (diameter.C + 1)/2, N = (diameter.N + 1)/2}}, solvent = {epsilon = 1.0, r_cut = 3^(1/6), sigma = 1.0}}
9
+-- Number of dimer spheres.
10
+Nd = 2
11
+-- Solvent number density.
12
+rho = 0.8
13
+-- Number of solvent particles excluding volume of walls.
14
+Ns = math.floor(rho * L[1] * L[2] * (L[3] - 1.0))
15
+-- Cutoff radius of potential in units of σ.
16
+r_cut = 2^(1/6)
17
+-- Minimum number of steps between solvent particle sorts.
18
+sort = {interval = 100}
19
+-- Masses of neutrally buoyant dimer spheres.
20
+mass = {C = rho * (math.pi/6) * diameter.C^3, N = rho * (math.pi/6) * diameter.N^3}
21
+-- Distance of dimer spheres.
22
+bond = ((diameter.C + diameter.N)/2 + 1) * r_cut
23
+-- Potential well depths of dimer spheres with solvent particles.
24
+epsilon = {CA = 1.0, CB = 1.0, NA = 1.0, NB = 0.1}
25
+-- Potential diameters.
26
+sigma = {C = (diameter.C + 1)/2, N = (diameter.N + 1)/2}
27
+-- Forward and backward reaction of solvent particles.
28
+reaction = {rate = 0.001}
29
+-- Verlet neighbour list skin.
30
+skin = 0.5
31
+-- Number of placeholders per neighbour list.
32
+nlist_size = {ss = 25, ds = 550}
33
+-- Number of bins per dimension for solvent particles.
34
+nbin = {L[1]/2, L[2]/2, L[3]/2}
35
+-- Number of placeholders per bin.
36
+bin_size = 20
37
+-- Integration time-step.
38
+timestep = 0.001
39
+-- Number of steps.
40
+nsteps = 100000
41
+-- Thermodynamic variables.
42
+thermo = {interval = 100, datatype = "double"}
43
+-- Dimer trajectory.
44
+trajd = {interval = 100, datatype = "double", species = {"C", "N"}}
45
+-- Solvent trajectory.
46
+trajs = {interval = nsteps, datatype = "double", species = {"A", "B"}}
47
+-- Count solvent species.
48
+species = {interval = 1000, species = {"A", "B"}}
49
+-- OpenCL platforms.
50
+platform = 1
51
+-- OpenCL device per platform.
52
+device = 1
53
+-- Random number generator seed.
54
+seed = "0x5ed7aa41fd2cd724"
55
+-- H5MD input filename.
56
+input = "../thermostat/confined_dimer.h5"
57
+-- Read trajectory at step.
58
+step = 100000
59
+-- H5MD output filename.
60
+output = "confined_dimer.h5"
61
+-- H5MD file author.
62
+author = {name = "B.J. Alder", email = "bjalder@example.org"}
63
+-- H5MD file creator.
64
+creator = {name = "confined_dimer.lua", version = "1.0"}

+ 152
- 0
examples/confined_dimer/equilibration/confined_dimer.lua View File

@@ -0,0 +1,152 @@
1
+#!/usr/bin/env luajit
2
+------------------------------------------------------------------------------
3
+-- Equilibration phase of confined-dimer simulation.
4
+-- Copyright © 2013–2015 Peter Colberg.
5
+-- Distributed under the MIT license. (See accompanying file LICENSE.)
6
+------------------------------------------------------------------------------
7
+
8
+-- Check uses of undeclared global variables.
9
+require("strict")
10
+
11
+local nm = {
12
+  config = require("nanomotor.config"),
13
+  box = require("nanomotor.box"),
14
+  domain = require("nanomotor.domain"),
15
+  integrate = require("nanomotor.integrate"),
16
+  observe = require("nanomotor.observe"),
17
+  h5md = require("nanomotor.h5md"),
18
+  trajd = require("nanomotor.trajd"),
19
+  trajs = require("nanomotor.trajs"),
20
+  thermo = require("nanomotor.thermo"),
21
+  species = require("nanomotor.species"),
22
+}
23
+
24
+local cl = require("opencl")
25
+local hdf5 = require("hdf5")
26
+local ffi = require("ffi")
27
+local syscall = require("syscall")
28
+local json = require("cjson")
29
+
30
+-- Read configuration file.
31
+local args = nm.config(arg[1])
32
+
33
+-- Seed random number generator.
34
+math.randomseed(args.seed)
35
+
36
+-- Create OpenCL context on chosen OpenCL device.
37
+local platform = cl.get_platforms()[args.platform]
38
+print(platform:get_info("version"))
39
+print(platform:get_info("name"))
40
+print(platform:get_info("vendor"))
41
+print(platform:get_info("extensions"))
42
+local device = platform:get_devices()[args.device]
43
+print(device:get_info("name"))
44
+print(device:get_info("vendor"))
45
+print(device:get_info("extensions"))
46
+local context = cl.create_context({device})
47
+local queue = context:create_command_queue(device)
48
+
49
+-- Create simulation box with given boundary conditions.
50
+local box = nm.box(args)
51
+
52
+-- Allocate zero-filled particle buffers.
53
+local dom = nm.domain(queue, box, args)
54
+
55
+-- Read particle coordinates from H5MD input file.
56
+local file = hdf5.open_file(args.input)
57
+do
58
+  local group = file:open_group("particles/dimer")
59
+  local obs_rd = nm.h5md.open_observable(group:open_group("position"))
60
+  local obs_vd = nm.h5md.open_observable(group:open_group("velocity"))
61
+  local space_rd = hdf5.create_simple_space({dom.Nd, 4})
62
+  local space_vd = hdf5.create_simple_space({dom.Nd, 4})
63
+  space_rd:select_hyperslab("set", {0, 0}, nil, {dom.Nd, 3})
64
+  space_vd:select_hyperslab("set", {0, 0}, nil, {dom.Nd, 3})
65
+  obs_rd:read(args.step, dom.rd, hdf5.double, space_rd)
66
+  obs_vd:read(args.step, dom.vd, hdf5.double, space_vd)
67
+  queue:enqueue_write_buffer(dom.d_rd, true, dom.rd)
68
+  queue:enqueue_write_buffer(dom.d_vd, true, dom.vd)
69
+end
70
+do
71
+  local group = file:open_group("particles/solvent")
72
+  local obs_rs = nm.h5md.open_observable(group:open_group("position"))
73
+  local obs_vs = nm.h5md.open_observable(group:open_group("velocity"))
74
+  local obs_species = nm.h5md.open_observable(group:open_group("species"))
75
+  local space_rs = hdf5.create_simple_space({dom.Ns, 4})
76
+  local space_vs = hdf5.create_simple_space({dom.Ns, 4})
77
+  local space_species = hdf5.create_simple_space({dom.Ns, 4})
78
+  space_rs:select_hyperslab("set", {0, 0}, nil, {dom.Ns, 3})
79
+  space_vs:select_hyperslab("set", {0, 0}, nil, {dom.Ns, 3})
80
+  space_species:select_hyperslab("set", {0, 0}, nil, {dom.Ns, 1})
81
+  local dtype_species = hdf5.int8:enum_create()
82
+  dtype_species:enum_insert("A", ffi.new("cl_char[1]", 0))
83
+  dtype_species:enum_insert("B", ffi.new("cl_char[1]", 1))
84
+  obs_rs:read(args.step, dom.rs, hdf5.double, space_rs)
85
+  obs_vs:read(args.step, dom.vs, hdf5.double, space_vs)
86
+  obs_species:read(args.step, dom.species, dtype_species, space_species)
87
+  queue:enqueue_write_buffer(dom.d_rs, true, dom.rs)
88
+  queue:enqueue_write_buffer(dom.d_vs, true, dom.vs)
89
+  queue:enqueue_write_buffer(dom.d_species, true, dom.species)
90
+end
91
+file:close()
92
+
93
+-- Create MD integrator.
94
+local integrate = nm.integrate(dom, args)
95
+
96
+-- Create H5MD output file with temporary filename.
97
+local file = hdf5.create_file(args.output..".tmp")
98
+-- Write H5MD metadata group.
99
+nm.h5md.write_meta(file, args)
100
+
101
+-- Serialize simulation parameters to JSON string attribute.
102
+do
103
+  local space = hdf5.create_space("scalar")
104
+  local dtype = hdf5.c_s1:copy()
105
+  local s = json.encode(args)
106
+  dtype:set_size(#s+1)
107
+  local attr = file:create_attribute("parameters", dtype, space)
108
+  attr:write(ffi.cast("const char *", s), dtype)
109
+end
110
+
111
+-- Create observables.
112
+local observables = {}
113
+if args.thermo then table.insert(observables, nm.thermo(integrate, file, args.thermo)) end
114
+if args.trajd then table.insert(observables, nm.trajd(integrate, file, args.trajd)) end
115
+if args.trajs then table.insert(observables, nm.trajs(integrate, file, args.trajs)) end
116
+if args.species then table.insert(observables, nm.species(integrate, file, args.species)) end
117
+
118
+-- Periodically log step, total energy, and remaining run-time.
119
+do
120
+  local interval = 1000
121
+  local self = {step = 0}
122
+
123
+  function self.observe(nsteps)
124
+    local t1 = syscall.clock_gettime("monotonic")
125
+    while self.step <= nsteps do
126
+      coroutine.yield()
127
+      local en = integrate.compute_en()
128
+      local t2 = syscall.clock_gettime("monotonic")
129
+      local t = (t2.time - t1.time) / interval * (nsteps - self.step)
130
+      local u = "s"
131
+      if t >= 60 then t, u = t/60, "m" end
132
+      if t >= 60 then t, u = t/60, "h" end
133
+      io.write("#", self.step, " ", en.en_tot, " ", ("%.1f"):format(t), u, "\n")
134
+      io.flush()
135
+      self.step = self.step + interval
136
+      t1 = t2
137
+    end
138
+  end
139
+
140
+  table.insert(observables, self)
141
+end
142
+
143
+-- Run simulation.
144
+nm.observe(integrate, observables, args.nsteps)
145
+
146
+-- Synchronize output files to disk.
147
+file:flush()
148
+assert((syscall.fsync(ffi.cast("int *", file:get_vfd_handle())[0])))
149
+file:close()
150
+
151
+-- Rename H5MD output file after successful simulation.
152
+assert((syscall.rename(args.output..".tmp", args.output)))

+ 74
- 0
examples/confined_dimer/production/config.lua View File

@@ -0,0 +1,74 @@
1
+-- Edge lengths of simulation domain.
2
+L = {50, 50, 20}
3
+-- Periodic boundary conditions.
4
+periodic = {true, true, false}
5
+-- Diameters of dimer spheres.
6
+diameter = {C = 4.0, N = 8.0}
7
+-- Lennard-Jones wall potential.
8
+wall = {dimer = {epsilon = 5.0, r_cut = 3^(1/6), sigma = {C = (diameter.C + 1)/2, N = (diameter.N + 1)/2}}, solvent = {epsilon = 1.0, r_cut = 3^(1/6), sigma = 1.0}}
9
+-- Number of dimer spheres.
10
+Nd = 2
11
+-- Solvent number density.
12
+rho = 0.8
13
+-- Number of solvent particles excluding volume of walls.
14
+Ns = math.floor(rho * L[1] * L[2] * (L[3] - 1.0))
15
+-- Cutoff radius of potential in units of σ.
16
+r_cut = 2^(1/6)
17
+-- Minimum number of steps between solvent particle sorts.
18
+sort = {interval = 100}
19
+-- Masses of neutrally buoyant dimer spheres.
20
+mass = {C = rho * (math.pi/6) * diameter.C^3, N = rho * (math.pi/6) * diameter.N^3}
21
+-- Distance of dimer spheres.
22
+bond = ((diameter.C + diameter.N)/2 + 1) * r_cut
23
+-- Potential well depths of dimer spheres with solvent particles.
24
+epsilon = {CA = 1.0, CB = 1.0, NA = 1.0, NB = 0.1}
25
+-- Potential diameters.
26
+sigma = {C = (diameter.C + 1)/2, N = (diameter.N + 1)/2}
27
+-- Forward and backward reaction of solvent particles.
28
+reaction = {rate = 0.001}
29
+-- Verlet neighbour list skin.
30
+skin = 0.5
31
+-- Number of placeholders per neighbour list.
32
+nlist_size = {ss = 25, ds = 550}
33
+-- Number of bins per dimension for solvent particles.
34
+nbin = {L[1]/2, L[2]/2, L[3]/2}
35
+-- Number of placeholders per bin.
36
+bin_size = 20
37
+-- Integration time-step.
38
+timestep = 0.001
39
+-- Number of steps.
40
+nsteps = 100000
41
+-- Thermodynamic variables.
42
+thermo = {interval = 100, datatype = "double"}
43
+-- Dimer trajectory.
44
+trajd = {interval = 100, datatype = "float", species = {"C", "N"}}
45
+-- Solvent trajectory.
46
+trajs = {interval = 10000, datatype = "float", species = {"A", "B"}}
47
+-- Count solvent species.
48
+species = {interval = 1000, species = {"A", "B"}}
49
+-- Mean-square displacement of dimer.
50
+msdd = {interval = 10, blockcount = 6, blocksize = 100, blockshift = 10, separation = 1000}
51
+-- Velocity autocorrelation function of dimer.
52
+vafd = {interval = 10, blockcount = 6, blocksize = 100, blockshift = 10, separation = 1000}
53
+-- Orientation autocorrelation of dimer.
54
+orient = {interval = 10, blockcount = 6, blocksize = 100, blockshift = 10, separation = 1000}
55
+-- Wall density function.
56
+walldf = {interval = 1000, nbin = 1000, cutoff = (diameter.C + diameter.N) / 2 + bond}
57
+-- OpenCL platforms.
58
+platform = 1
59
+-- OpenCL device per platform.
60
+device = 1
61
+-- Random number generator seed.
62
+seed = "0xcb22f5cf04a094bf"
63
+-- H5MD input filename.
64
+input = "../equilibration/confined_dimer.h5"
65
+-- Read trajectory at step.
66
+step = 100000
67
+-- H5MD output filename.
68
+output = "confined_dimer.h5"
69
+-- Snapshot program state.
70
+snapshot = {interval = 10000, output = "confined_dimer-snapshot.h5"}
71
+-- H5MD file author.
72
+author = {name = "B.J. Alder", email = "bjalder@example.org"}
73
+-- H5MD file creator.
74
+creator = {name = "confined_dimer.lua", version = "1.0"}

+ 165
- 0
examples/confined_dimer/production/confined_dimer.lua View File

@@ -0,0 +1,165 @@
1
+#!/usr/bin/env luajit
2
+------------------------------------------------------------------------------
3
+-- Production phase of confined-dimer simulation.
4
+-- Copyright © 2013–2015 Peter Colberg.
5
+-- Distributed under the MIT license. (See accompanying file LICENSE.)
6
+------------------------------------------------------------------------------
7
+
8
+-- Check uses of undeclared global variables.
9
+require("strict")
10
+
11
+local nm = {
12
+  config = require("nanomotor.config"),
13
+  box = require("nanomotor.box"),
14
+  domain = require("nanomotor.domain"),
15
+  integrate = require("nanomotor.integrate"),
16
+  observe = require("nanomotor.observe"),
17
+  snapshot = require("nanomotor.snapshot"),
18
+  h5md = require("nanomotor.h5md"),
19
+  trajd = require("nanomotor.trajd"),
20
+  trajs = require("nanomotor.trajs"),
21
+  thermo = require("nanomotor.thermo"),
22
+  species = require("nanomotor.species"),
23
+  corr = require("nanomotor.correlation"),
24
+  msdd = require("msdd"),
25
+  vafd = require("vafd"),
26
+  orient = require("orient"),
27
+  walldf = require("walldf"),
28
+}
29
+
30
+local cl = require("opencl")
31
+local hdf5 = require("hdf5")
32
+local ffi = require("ffi")
33
+local syscall = require("syscall")
34
+local json = require("cjson")
35
+
36
+-- Read configuration file.
37
+local args = nm.config(arg[1])
38
+
39
+-- Seed random number generator.
40
+math.randomseed(args.seed)
41
+
42
+-- Create OpenCL context on chosen OpenCL device.
43
+local platform = cl.get_platforms()[args.platform]
44
+print(platform:get_info("version"))
45
+print(platform:get_info("name"))
46
+print(platform:get_info("vendor"))
47
+print(platform:get_info("extensions"))
48
+local device = platform:get_devices()[args.device]
49
+print(device:get_info("name"))
50
+print(device:get_info("vendor"))
51
+print(device:get_info("extensions"))
52
+local context = cl.create_context({device})
53
+local queue = context:create_command_queue(device)
54
+
55
+-- Create simulation box with given boundary conditions.
56
+local box = nm.box(args)
57
+
58
+-- Allocate zero-filled particle buffers.
59
+local dom = nm.domain(queue, box, args)
60
+
61
+-- Read particle coordinates from H5MD input file.
62
+local file = hdf5.open_file(args.input)
63
+do
64
+  local group = file:open_group("particles/dimer")
65
+  local obs_rd = nm.h5md.open_observable(group:open_group("position"))
66
+  local obs_vd = nm.h5md.open_observable(group:open_group("velocity"))
67
+  local space_rd = hdf5.create_simple_space({dom.Nd, 4})
68
+  local space_vd = hdf5.create_simple_space({dom.Nd, 4})
69
+  space_rd:select_hyperslab("set", {0, 0}, nil, {dom.Nd, 3})
70
+  space_vd:select_hyperslab("set", {0, 0}, nil, {dom.Nd, 3})
71
+  obs_rd:read(args.step, dom.rd, hdf5.double, space_rd)
72
+  obs_vd:read(args.step, dom.vd, hdf5.double, space_vd)
73
+  queue:enqueue_write_buffer(dom.d_rd, true, dom.rd)
74
+  queue:enqueue_write_buffer(dom.d_vd, true, dom.vd)
75
+end
76
+do
77
+  local group = file:open_group("particles/solvent")
78
+  local obs_rs = nm.h5md.open_observable(group:open_group("position"))
79
+  local obs_vs = nm.h5md.open_observable(group:open_group("velocity"))
80
+  local obs_species = nm.h5md.open_observable(group:open_group("species"))
81
+  local space_rs = hdf5.create_simple_space({dom.Ns, 4})
82
+  local space_vs = hdf5.create_simple_space({dom.Ns, 4})
83
+  local space_species = hdf5.create_simple_space({dom.Ns, 4})
84
+  space_rs:select_hyperslab("set", {0, 0}, nil, {dom.Ns, 3})
85
+  space_vs:select_hyperslab("set", {0, 0}, nil, {dom.Ns, 3})
86
+  space_species:select_hyperslab("set", {0, 0}, nil, {dom.Ns, 1})
87
+  local dtype_species = hdf5.int8:enum_create()
88
+  dtype_species:enum_insert("A", ffi.new("cl_char[1]", 0))
89
+  dtype_species:enum_insert("B", ffi.new("cl_char[1]", 1))
90
+  obs_rs:read(args.step, dom.rs, hdf5.double, space_rs)
91
+  obs_vs:read(args.step, dom.vs, hdf5.double, space_vs)
92
+  obs_species:read(args.step, dom.species, dtype_species, space_species)
93
+  queue:enqueue_write_buffer(dom.d_rs, true, dom.rs)
94
+  queue:enqueue_write_buffer(dom.d_vs, true, dom.vs)
95
+  queue:enqueue_write_buffer(dom.d_species, true, dom.species)
96
+end
97
+file:close()
98
+
99
+-- Create MD integrator.
100
+local integrate = nm.integrate(dom, args)
101
+
102
+-- Create H5MD output file with temporary filename.
103
+local file = hdf5.create_file(args.output..".tmp")
104
+-- Write H5MD metadata group.
105
+nm.h5md.write_meta(file, args)
106
+
107
+-- Serialize simulation parameters to JSON string attribute.
108
+do
109
+  local space = hdf5.create_space("scalar")
110
+  local dtype = hdf5.c_s1:copy()
111
+  local s = json.encode(args)
112
+  dtype:set_size(#s+1)
113
+  local attr = file:create_attribute("parameters", dtype, space)
114
+  attr:write(ffi.cast("const char *", s), dtype)
115
+end
116
+
117
+-- Create observables.
118
+local observables = {}
119
+if args.snapshot then table.insert(observables, nm.snapshot(dom, observables, args.snapshot)) end
120
+if args.thermo then table.insert(observables, nm.thermo(integrate, file, args.thermo)) end
121
+if args.trajd then table.insert(observables, nm.trajd(integrate, file, args.trajd)) end
122
+if args.trajs then table.insert(observables, nm.trajs(integrate, file, args.trajs)) end
123
+if args.species then table.insert(observables, nm.species(integrate, file, args.species)) end
124
+if args.msdd then table.insert(observables, nm.corr(integrate, nm.msdd(dom, file), args.msdd)) end
125
+if args.vafd then table.insert(observables, nm.corr(integrate, nm.vafd(dom, file), args.vafd)) end
126
+if args.orient then table.insert(observables, nm.corr(integrate, nm.orient(dom, file), args.orient)) end
127
+if args.walldf then table.insert(observables, nm.walldf(integrate, file, args.walldf)) end
128
+
129
+-- Periodically log step, total energy, and remaining run-time.
130
+do
131
+  local interval = 10000
132
+  local self = {step = 0}
133
+
134
+  function self.observe(nsteps)
135
+    local t1 = syscall.clock_gettime("monotonic")
136
+    while self.step <= nsteps do
137
+      coroutine.yield()
138
+      local en = integrate.compute_en()
139
+      local t2 = syscall.clock_gettime("monotonic")
140
+      local t = (t2.time - t1.time) / interval * (nsteps - self.step)
141
+      local u = "s"
142
+      if t >= 60 then t, u = t/60, "m" end
143
+      if t >= 60 then t, u = t/60, "h" end
144
+      io.write("#", self.step, " ", en.en_tot, " ", ("%.1f"):format(t), u, "\n")
145
+      io.flush()
146
+      self.step = self.step + interval
147
+      t1 = t2
148
+    end
149
+  end
150
+
151
+  table.insert(observables, self)
152
+end
153
+
154
+-- Run simulation.
155
+nm.observe(integrate, observables, args.nsteps)
156
+
157
+-- Synchronize output files to disk.
158
+file:flush()
159
+assert((syscall.fsync(ffi.cast("int *", file:get_vfd_handle())[0])))
160
+file:close()
161
+
162
+-- Rename H5MD output file after successful simulation.
163
+assert((syscall.rename(args.output..".tmp", args.output)))
164
+-- Remove snapshot file.
165
+if args.snapshot then assert((syscall.unlink(args.snapshot.output))) end

+ 69
- 0
examples/confined_dimer/production/msdd.lua View File

@@ -0,0 +1,69 @@
1
+------------------------------------------------------------------------------
2
+-- Mean-square displacement of dimer.
3
+-- Copyright © 2013–2015 Peter Colberg.
4
+-- Distributed under the MIT license. (See accompanying file LICENSE.)
5
+------------------------------------------------------------------------------
6
+
7
+local hdf5 = require("hdf5")
8
+local ffi = require("ffi")
9
+
10
+-- Cache C types.
11
+local vector_n = ffi.typeof("struct { double x, y, z; }[?]")
12
+
13
+return function(dom, file)
14
+  local box = dom.box
15
+
16
+  local self = {
17
+    group = file:create_group("dynamics/dimer/mean_square_displacement"),
18
+    dims = {3},
19
+  }
20
+
21
+  function self.sample()
22
+    local r_cm = vector_n(dom.Nd/2)
23
+    dom.queue:enqueue_read_buffer(dom.d_rd, true, dom.rd)
24
+    dom.queue:enqueue_read_buffer(dom.d_imd, true, dom.imd)
25
+    local rd, imd, L = dom.rd, dom.imd, box.L
26
+    local m_N = 1 / (1 + dom.mass.C/dom.mass.N)
27
+    for i = 0, dom.Nd/2-1 do
28
+      local dx = rd[2*i+1].x - rd[2*i].x
29
+      local dy = rd[2*i+1].y - rd[2*i].y
30
+      local dz = rd[2*i+1].z - rd[2*i].z
31
+      dx, dy, dz = box.mindist(dx, dy, dz)
32
+      r_cm[i].x = rd[2*i].x + L[1]*imd[2*i].x + m_N*dx
33
+      r_cm[i].y = rd[2*i].y + L[2]*imd[2*i].y + m_N*dy
34
+      r_cm[i].z = rd[2*i].z + L[3]*imd[2*i].z + m_N*dz
35
+    end
36
+    return r_cm
37
+  end
38
+
39
+  function self.correlate(r, r0, result)
40
+    local x, y, z = 0, 0, 0
41
+    for i = 0, dom.Nd/2-1 do
42
+      local dx = r[i].x - r0[i].x
43
+      local dy = r[i].y - r0[i].y
44
+      local dz = r[i].z - r0[i].z
45
+      x = x + (dx*dx - x) / (i+1)
46
+      y = y + (dy*dy - y) / (i+1)
47
+      z = z + (dz*dz - z) / (i+1)
48
+    end
49
+    result[0](x); result[1](y); result[2](z)
50
+  end
51
+
52
+  do
53
+    local space_r_cm = hdf5.create_simple_space({dom.Nd/2, 3})
54
+
55
+    function self.snapshot(group, r_cm)
56
+      local dset = group:create_dataset("r_cm", hdf5.double, space_r_cm)
57
+      dset:write(r_cm, hdf5.double, space_r_cm)
58
+    end
59
+
60
+    function self.restore(group)
61
+      local dset = group:open_dataset("r_cm")
62
+      local r_cm = vector_n(dom.Nd/2)
63
+      dset:read(r_cm, hdf5.double, space_r_cm)
64
+      return r_cm
65
+    end
66
+  end
67
+
68
+  return self
69
+end

+ 65
- 0
examples/confined_dimer/production/orient.lua View File

@@ -0,0 +1,65 @@
1
+------------------------------------------------------------------------------
2
+-- Orientation autocorrelation of dimer.
3
+-- Copyright © 2013–2015 Peter Colberg.
4
+-- Distributed under the MIT license. (See accompanying file LICENSE.)
5
+------------------------------------------------------------------------------
6
+
7
+local hdf5 = require("hdf5")
8
+local ffi = require("ffi")
9
+
10
+-- Cache library functions.
11
+local sqrt = math.sqrt
12
+-- Cache C types.
13
+local vector_n = ffi.typeof("struct { double x, y, z; }[?]")
14
+
15
+return function(dom, file)
16
+  local box = dom.box
17
+
18
+  local self = {
19
+    group = file:create_group("dynamics/dimer/orientation_autocorrelation"),
20
+    dims = {3},
21
+  }
22
+
23
+  function self.sample()
24
+    local d = vector_n(dom.Nd/2)
25
+    dom.queue:enqueue_read_buffer(dom.d_rd, true, dom.rd)
26
+    local rd = dom.rd
27
+    for i = 0, dom.Nd/2-1 do
28
+      local dx = rd[2*i+1].x - rd[2*i].x
29
+      local dy = rd[2*i+1].y - rd[2*i].y
30
+      local dz = rd[2*i+1].z - rd[2*i].z
31
+      dx, dy, dz = box.mindist(dx, dy, dz)
32
+      local dn = 1 / sqrt(dx*dx + dy*dy + dz*dz)
33
+      d[i].x, d[i].y, d[i].z = dx*dn, dy*dn, dz*dn
34
+    end
35
+    return d
36
+  end
37
+
38
+  function self.correlate(d, d0, result)
39
+    local x, y, z = 0, 0, 0
40
+    for i = 0, dom.Nd/2-1 do
41
+      x = x + (d[i].x*d0[i].x - x) / (i+1)
42
+      y = y + (d[i].y*d0[i].y - y) / (i+1)
43
+      z = z + (d[i].z*d0[i].z - z) / (i+1)
44
+    end
45
+    result[0](x); result[1](y); result[2](z)
46
+  end
47
+
48
+  do
49
+    local space_d = hdf5.create_simple_space({dom.Nd/2, 3})
50
+
51
+    function self.snapshot(group, d)
52
+      local dset = group:create_dataset("d", hdf5.double, space_d)
53
+      dset:write(d, hdf5.double, space_d)
54
+    end
55
+
56
+    function self.restore(group)
57
+      local dset = group:open_dataset("d")
58
+      local d = vector_n(dom.Nd/2)
59
+      dset:read(d, hdf5.double, space_d)
60
+      return d
61
+    end
62
+  end
63
+
64
+  return self
65
+end

+ 60
- 0
examples/confined_dimer/production/vafd.lua View File

@@ -0,0 +1,60 @@
1
+------------------------------------------------------------------------------
2
+-- Velocity autocorrelation function of dimer.
3
+-- Copyright © 2013–2015 Peter Colberg.
4
+-- Distributed under the MIT license. (See accompanying file LICENSE.)
5
+------------------------------------------------------------------------------
6
+
7
+local hdf5 = require("hdf5")
8
+local ffi = require("ffi")
9
+
10
+-- Cache C types.
11
+local vector_n = ffi.typeof("struct { double x, y, z; }[?]")
12
+
13
+return function(dom, file)
14
+  local self = {
15
+    group = file:create_group("dynamics/dimer/velocity_autocorrelation"),
16
+    dims = {3},
17
+  }
18
+
19
+  function self.sample()
20
+    local v_cm = vector_n(dom.Nd/2)
21
+    dom.queue:enqueue_read_buffer(dom.d_vd, true, dom.vd)
22
+    local vd = dom.vd
23
+    local m_C = 1 / (1 + dom.mass.N/dom.mass.C)
24
+    local m_N = 1 / (1 + dom.mass.C/dom.mass.N)
25
+    for i = 0, dom.Nd/2-1 do
26
+      v_cm[i].x = m_C*vd[2*i].x + m_N*vd[2*i+1].x
27
+      v_cm[i].y = m_C*vd[2*i].y + m_N*vd[2*i+1].y
28
+      v_cm[i].z = m_C*vd[2*i].z + m_N*vd[2*i+1].z
29
+    end
30
+    return v_cm
31
+  end
32
+
33
+  function self.correlate(v, v0, result)
34
+    local x, y, z = 0, 0, 0
35
+    for i = 0, dom.Nd/2-1 do
36
+      x = x + (v[i].x*v0[i].x - x) / (i+1)
37
+      y = y + (v[i].y*v0[i].y - y) / (i+1)
38
+      z = z + (v[i].z*v0[i].z - z) / (i+1)
39
+    end
40
+    result[0](x); result[1](y); result[2](z)
41
+  end
42
+
43
+  do
44
+    local space_v_cm = hdf5.create_simple_space({dom.Nd/2, 3})
45
+
46
+    function self.snapshot(group, v_cm)
47
+      local dset = group:create_dataset("v_cm", hdf5.double, space_v_cm)
48
+      dset:write(v_cm, hdf5.double, space_v_cm)
49
+    end
50
+
51
+    function self.restore(group)
52
+      local dset = group:open_dataset("v_cm")
53
+      local v_cm = vector_n(dom.Nd/2)
54
+      dset:read(v_cm, hdf5.double, space_v_cm)
55
+      return v_cm
56
+    end
57
+  end
58
+
59
+  return self
60
+end

+ 126
- 0
examples/confined_dimer/production/walldf.lua View File

@@ -0,0 +1,126 @@
1
+------------------------------------------------------------------------------
2
+-- Wall density function.
3
+-- Copyright © 2013–2015 Peter Colberg.
4
+-- Distributed under the MIT license. (See accompanying file LICENSE.)
5
+------------------------------------------------------------------------------
6
+
7
+local hdf5 = require("hdf5")
8
+local ffi = require("ffi")
9
+
10
+-- Cache library functions.
11
+local floor, max, min = math.floor, math.max, math.min
12
+-- Cache C types.
13
+local scalar = hdf5.create_space("scalar")
14
+
15
+return function(integrate, file, args)
16
+  local dom = integrate.dom
17
+  local box = dom.box
18
+  local L = box.L
19
+  local nbin, cutoff = args.nbin, args.cutoff
20
+  local interval = args.interval
21
+
22
+  local bin = ffi.new("struct { uint64_t F, N; }[?]", nbin)
23
+  local count = ffi.new("uint64_t[1]")
24
+
25
+  local function sample()
26
+    dom.queue:enqueue_read_buffer(dom.d_rd, true, dom.rd)
27
+    dom.queue:enqueue_read_buffer(dom.d_rs, true, dom.rs)
28
+    dom.queue:enqueue_read_buffer(dom.d_species, true, dom.species)
29
+    local rd, rs, species = dom.rd, dom.rs, dom.species
30
+    local m_N = 1 / (1 + dom.mass.C/dom.mass.N)
31
+    local zx = rd[1].x - rd[0].x
32
+    local zy = rd[1].y - rd[0].y
33
+    local zz = rd[1].z - rd[0].z
34
+    zx, zy, zz = box.mindist(zx, zy, zz)
35
+    local rx = rd[0].x + m_N*zx
36
+    local ry = rd[0].y + m_N*zy
37
+    local rz = rd[0].z + m_N*zz
38
+    for i = 0, dom.Ns-1 do
39
+      local bin_z = max(0, min(nbin-1, floor(rs[i].z*(nbin/L[3]))))
40
+      local dx = rs[i].x - rx
41
+      local dy = rs[i].y - ry
42
+      dx, dy = box.mindist(dx, dy)
43
+      if dx*dx + dy*dy > cutoff^2 then
44
+        bin[bin_z].F = bin[bin_z].F + 1
45
+      else
46
+        bin[bin_z].N = bin[bin_z].N + 1
47
+      end
48
+    end
49
+    count[0] = count[0] + 1
50
+  end
51
+
52
+  local group = file:create_group("structure/solvent/wall_density")
53
+
54
+  local delta_z = box.L[3] / nbin
55
+  do
56
+    local z = ffi.new("float[?]", nbin)
57
+    for i = 0, nbin-1 do
58
+      z[i] = (i+0.5)*delta_z
59
+    end
60
+    local space_z = hdf5.create_simple_space({nbin})
61
+    local dset_z = group:create_dataset("distance", hdf5.float, space_z)
62
+    dset_z:write(z, hdf5.float)
63
+  end
64
+
65
+  do
66
+    local dtype_region = hdf5.c_s1:copy()
67
+    dtype_region:set_size(5)
68
+    local space_region = hdf5.create_simple_space({2})
69
+    local attr_region = group:create_attribute("region", dtype_region, space_region)
70
+    attr_region:write(ffi.new("char[2][5]", "far", "near"), dtype_region)
71
+  end
72
+
73
+  do
74
+    local attr_cutoff = group:create_attribute("cutoff", hdf5.double, scalar)
75
+    attr_cutoff:write(ffi.new("double[1]", cutoff), hdf5.double)
76
+  end
77
+
78
+  local write do
79
+    local rho = ffi.new("struct { float F, N; }[?]", nbin)
80
+    local space_rho = hdf5.create_simple_space({nbin, 2})
81
+    local dset_rho = group:create_dataset("value", hdf5.float, space_rho)
82
+    local attr_count = group:create_attribute("count", hdf5.uint64, scalar)
83
+
84
+    function write()
85
+      local norm_F = tonumber(count[0]) * delta_z * (L[1]*L[2] - math.pi*cutoff*cutoff)
86
+      local norm_N = tonumber(count[0]) * delta_z * (math.pi*cutoff*cutoff)
87
+      for i = 0, nbin-1 do
88
+        rho[i].F = tonumber(bin[i].F) / norm_F
89
+        rho[i].N = tonumber(bin[i].N) / norm_N
90
+      end
91
+      dset_rho:write(rho, hdf5.float)
92
+      attr_count:write(count, hdf5.uint64)
93
+    end
94
+  end
95
+
96
+  local self = {step = 0}
97
+
98
+  function self.observe(nsteps)
99
+    while self.step <= nsteps do
100
+      coroutine.yield()
101
+      sample()
102
+      self.step = self.step + interval
103
+    end
104
+    write()
105
+  end
106
+
107
+  do
108
+    local space_bin = hdf5.create_simple_space({nbin, 2})
109
+
110
+    function self.snapshot(group)
111
+      local dset_bin = group:create_dataset("bin", hdf5.uint64, space_bin)
112
+      local attr_count = group:create_attribute("count", hdf5.uint64, scalar)
113
+      dset_bin:write(bin, hdf5.uint64, space_bin)
114
+      attr_count:write(count, hdf5.uint64)
115
+    end
116
+
117
+    function self.restore(group)
118
+      local dset_bin = group:open_dataset("bin")
119
+      local attr_count = group:open_attribute("count")
120
+      dset_bin:read(bin, hdf5.uint64, space_bin)
121
+      attr_count:read(count, hdf5.uint64)
122
+    end
123
+  end
124
+
125
+  return self
126
+end

+ 60
- 0
examples/confined_dimer/thermostat/config.lua View File

@@ -0,0 +1,60 @@
1
+-- Edge lengths of simulation domain.
2
+L = {50, 50, 20}
3
+-- Periodic boundary conditions.
4
+periodic = {true, true, false}
5
+-- Diameters of dimer spheres.
6
+diameter = {C = 4.0, N = 8.0}
7
+-- Lennard-Jones wall potential.
8
+wall = {dimer = {epsilon = 5.0, r_cut = 3^(1/6), sigma = {C = (diameter.C + 1)/2, N = (diameter.N + 1)/2}}, solvent = {epsilon = 1.0, r_cut = 3^(1/6), sigma = 1.0}}
9
+-- Number of dimer spheres.
10
+Nd = 2
11
+-- Solvent number density.
12
+rho = 0.8
13
+-- Number of solvent particles excluding volume of walls.
14
+Ns = math.floor(rho * L[1] * L[2] * (L[3] - 1.0))
15
+-- Cutoff radius of potential in units of σ.
16
+r_cut = 2^(1/6)
17
+-- Minimum number of steps between solvent particle sorts.
18
+sort = {interval = 100}
19
+-- Lattice constant and region for solvent particles.
20
+lattice = {a = 1.6, slab = {{0, 0, 1}, {L[1], L[2], L[3]-1}}}
21
+-- Solvent temperature.
22
+temp = 1
23
+-- Masses of neutrally buoyant dimer spheres.
24
+mass = {C = rho * (math.pi/6) * diameter.C^3, N = rho * (math.pi/6) * diameter.N^3}
25
+-- Distance of dimer spheres.
26
+bond = ((diameter.C + diameter.N)/2 + 1) * r_cut
27
+-- Potential well depths of dimer spheres with solvent particles.
28
+epsilon = {CA = 1.0, CB = 1.0, NA = 1.0, NB = 0.1}
29
+-- Potential diameters.
30
+sigma = {C = (diameter.C + 1)/2, N = (diameter.N + 1)/2}
31
+-- Verlet neighbour list skin.
32
+skin = 0.5
33
+-- Number of placeholders per neighbour list.
34
+nlist_size = {ss = 25, ds = 500}
35
+-- Number of bins per dimension for solvent particles.
36
+nbin = {L[1]/2, L[2]/2, L[3]/2}
37
+-- Number of placeholders per bin.
38
+bin_size = 20
39
+-- Integration time-step.
40
+timestep = 0.001
41
+-- Number of steps.
42
+nsteps = 100000
43
+-- Thermodynamic variables.
44
+thermo = {interval = 100, datatype = "double"}
45
+-- Dimer trajectory.
46
+trajd = {interval = 100, datatype = "double", species = {"C", "N"}}
47
+-- Solvent trajectory.
48
+trajs = {interval = nsteps, datatype = "double", species = {"A", "B"}}
49
+-- OpenCL platforms.
50
+platform = 1
51
+-- OpenCL device per platform.
52
+device = 1
53
+-- Random number generator seed.
54
+seed = "0xa66800f15c3fa1a"
55
+-- H5MD output filename.
56
+output = "confined_dimer.h5"
57
+-- H5MD file author.
58
+author = {name = "B.J. Alder", email = "bjalder@example.org"}
59
+-- H5MD file creator.
60
+creator = {name = "confined_dimer.lua", version = "1.0"}

+ 139
- 0
examples/confined_dimer/thermostat/confined_dimer.lua View File

@@ -0,0 +1,139 @@
1
+#!/usr/bin/env luajit
2
+------------------------------------------------------------------------------
3
+-- Thermostat phase of confined-dimer simulation.
4
+-- Copyright © 2013–2015 Peter Colberg.
5
+-- Distributed under the MIT license. (See accompanying file LICENSE.)
6
+------------------------------------------------------------------------------
7
+
8
+-- Check uses of undeclared global variables.
9
+require("strict")
10
+
11
+local nm = {
12
+  config = require("nanomotor.config"),
13
+  box = require("nanomotor.box"),
14
+  domain = require("nanomotor.domain"),
15
+  position = require("nanomotor.position"),
16
+  velocity = require("nanomotor.velocity"),
17
+  integrate = require("nanomotor.integrate"),
18
+  observe = require("nanomotor.observe"),
19
+  random = require("nanomotor.random"),
20
+  h5md = require("nanomotor.h5md"),
21
+  trajd = require("nanomotor.trajd"),
22
+  trajs = require("nanomotor.trajs"),
23
+  thermo = require("nanomotor.thermo"),
24
+}
25
+
26
+local cl = require("opencl")
27
+local hdf5 = require("hdf5")
28
+local ffi = require("ffi")
29
+local syscall = require("syscall")
30
+local json = require("cjson")
31
+
32
+-- Read configuration file.
33
+local args = nm.config(arg[1])
34
+
35
+-- Seed random number generator.
36
+math.randomseed(args.seed)
37
+
38
+-- Create OpenCL context on chosen OpenCL device.
39
+local platform = cl.get_platforms()[args.platform]
40
+print(platform:get_info("version"))
41
+print(platform:get_info("name"))
42
+print(platform:get_info("vendor"))
43
+print(platform:get_info("extensions"))
44
+local device = platform:get_devices()[args.device]
45
+print(device:get_info("name"))
46
+print(device:get_info("vendor"))
47
+print(device:get_info("extensions"))
48
+local context = cl.create_context({device})
49
+local queue = context:create_command_queue(device)
50
+
51
+-- Create simulation box with given boundary conditions.
52
+local box = nm.box(args)
53
+
54
+-- Allocate zero-filled particle buffers.
55
+local dom = nm.domain(queue, box, args)
56
+
57
+-- Place sphere-dimer at random position and orientation.
58
+do
59
+  local rd, L, bond = dom.rd, box.L, args.bond
60
+  for i = 0, dom.Nd-1, 2 do
61
+    rd[i].x = math.random()*L[1]
62
+    rd[i].y = math.random()*L[2]
63
+    rd[i].z = 0.5*L[3]
64
+    local x, y = nm.random.circle()
65
+    x = rd[i].x + bond*x
66
+    y = rd[i].y + bond*y
67
+    local z = rd[i].z
68
+    rd[i+1].x, rd[i+1].y, rd[i+1].z = box.minimage(x, y, z)
69
+  end
70
+  queue:enqueue_write_buffer(dom.d_rd, true, dom.rd)
71
+end
72
+
73
+-- Place solvent particles on fcc lattice excluding dimer-sphere volumes.
74
+nm.position.lattice(dom, args)()
75
+-- Pick solvent velocities from Maxwell-Boltzmann distribution.
76
+local velocity = nm.velocity.gaussian(dom, args)
77
+
78
+-- Create MD integrator.
79
+local integrate = nm.integrate(dom, args)
80
+
81
+-- Create H5MD output file with temporary filename.
82
+local file = hdf5.create_file(args.output..".tmp")
83
+-- Write H5MD metadata group.
84
+nm.h5md.write_meta(file, args)
85
+
86
+-- Serialize simulation parameters to JSON string attribute.
87
+do
88
+  local space = hdf5.create_space("scalar")
89
+  local dtype = hdf5.c_s1:copy()
90
+  local s = json.encode(args)
91
+  dtype:set_size(#s+1)
92
+  local attr = file:create_attribute("parameters", dtype, space)
93
+  attr:write(ffi.cast("const char *", s), dtype)
94
+end
95
+
96
+-- Create observables.
97
+local observables = {}
98
+if args.thermo then table.insert(observables, nm.thermo(integrate, file, args.thermo)) end
99
+if args.trajd then table.insert(observables, nm.trajd(integrate, file, args.trajd)) end
100
+if args.trajs then table.insert(observables, nm.trajs(integrate, file, args.trajs)) end
101
+
102
+-- Periodically log step, total energy, and remaining run-time.
103
+do
104
+  local interval = 1000
105
+  local self = {step = 0}
106
+
107
+  function self.observe(nsteps)
108
+    local t1 = syscall.clock_gettime("monotonic")
109
+    while self.step <= nsteps do
110
+      ffi.fill(dom.vd, dom.d_vd:get_info("size"))
111
+      queue:enqueue_write_buffer(dom.d_vd, true, dom.vd)
112
+      velocity()
113
+      coroutine.yield()
114
+      local en = integrate.compute_en()
115
+      local t2 = syscall.clock_gettime("monotonic")
116
+      local t = (t2.time - t1.time) / interval * (nsteps - self.step)
117
+      local u = "s"
118
+      if t >= 60 then t, u = t/60, "m" end
119
+      if t >= 60 then t, u = t/60, "h" end
120
+      io.write("#", self.step, " ", en.en_tot, " ", ("%.1f"):format(t), u, "\n")
121
+      io.flush()
122
+      self.step = self.step + interval
123
+      t1 = t2
124
+    end
125
+  end
126
+
127
+  table.insert(observables, self)
128
+end
129
+
130
+-- Run simulation.
131
+nm.observe(integrate, observables, args.nsteps)
132
+
133
+-- Synchronize output files to disk.
134
+file:flush()
135
+assert((syscall.fsync(ffi.cast("int *", file:get_vfd_handle())[0])))
136
+file:close()
137
+
138
+-- Rename H5MD output file after successful simulation.
139
+assert((syscall.rename(args.output..".tmp", args.output)))

+ 69
- 0
examples/confined_dimer/tools/dimer_mean_square_displacement.py View File

@@ -0,0 +1,69 @@
1
+#!/usr/bin/env python
2
+# -*- coding: utf-8 -*-
3
+#
4
+# Plot mean-square displacement of dimer.
5
+# Copyright © 2013–2015 Peter Colberg.
6
+# Distributed under the MIT license. (See accompanying file LICENSE.)
7
+#
8
+
9
+import matplotlib as mpl
10
+mpl.use("pgf")
11
+import matplotlib.pyplot as plt
12
+import numpy as np
13
+import h5py as h5
14
+import argparse
15
+
16
+parser = argparse.ArgumentParser()
17
+parser.add_argument("output", help="plot filename")
18
+parser.add_argument("input", nargs="+", help="input filenames")
19
+parser.add_argument("--figsize", type=float, nargs=2, default=(6., 4.), help="figure size")
20
+args = parser.parse_args()
21
+
22
+mpl.rc("pgf", texsystem="xelatex", preamble=(r"\usepackage{amsmath,bm,txfonts}",))
23
+mpl.rc("font", family="serif", serif=(), size=11)
24
+mpl.rc("legend", fontsize=11, numpoints=1)
25
+
26
+with h5.File(args.input[0], "r") as f:
27
+  g = f["dynamics/dimer/mean_square_displacement"]
28
+  time = g["time"][:]
29
+  index = g["count"][:] != 0
30
+  mean = np.zeros(shape=g["value"].shape)
31
+  error = np.zeros(shape=g["value"].shape)
32
+
33
+count = 0
34
+for fn in args.input:
35
+  with h5.File(fn, "r") as f:
36
+    g = f["dynamics/dimer/mean_square_displacement"]
37
+    assert np.all(time == g["time"])
38
+    value = g["value"][:]
39
+    count += 1
40
+    delta = value - mean
41
+    mean += delta / count
42
+    error += delta * (value - mean)
43
+
44
+error = np.sqrt(error / (count - 1) / count)
45
+
46
+time = np.reshape(time[:, 1:], (-1,))
47
+index = np.reshape(index[:, 1:], (-1, 3))
48
+mean = np.reshape(mean[:, 1:], (-1, 3))
49
+error = np.reshape(error[:, 1:], (-1, 3))
50
+
51
+order = time.argsort()
52
+time = time[order]
53
+index = index[order]
54
+mean = mean[order]
55
+error = error[order]
56
+
57
+fig = plt.figure(figsize=args.figsize)
58
+ax = fig.add_subplot(1, 1, 1)
59
+ax.errorbar(time[index[:, 0]], mean[index[:, 0], 0], error[index[:, 0], 0], errorevery=10, label="x")
60
+ax.errorbar(time[index[:, 1]], mean[index[:, 1], 1], error[index[:, 1], 1], errorevery=10, label="y")
61
+ax.errorbar(time[index[:, 2]], mean[index[:, 2], 2], error[index[:, 2], 2], errorevery=10, label="z")
62
+ax.set_xscale("log")
63
+ax.set_yscale("log")
64
+ax.minorticks_on()
65
+ax.set_xlabel(r"$t$")
66
+ax.set_ylabel(r"$\Delta L^2(t)$")
67
+ax.legend(loc="best", fancybox=True)
68
+fig.tight_layout(pad=0.1)
69
+fig.savefig(args.output)

+ 62
- 0
examples/confined_dimer/tools/dimer_orientation_autocorrelation.py View File

@@ -0,0 +1,62 @@
1
+#!/usr/bin/env python
2
+# -*- coding: utf-8 -*-
3
+#
4
+# Plot orientation autocorrelation of dimer.
5
+# Copyright © 2013–2015 Peter Colberg.
6
+# Distributed under the MIT license. (See accompanying file LICENSE.)
7
+#
8
+
9
+import matplotlib as mpl
10
+mpl.use("pgf")
11
+import matplotlib.pyplot as plt
12
+import numpy as np
13
+import h5py as h5
14
+import argparse
15
+
16
+parser = argparse.ArgumentParser()
17
+parser.add_argument("output", help="plot filename")
18
+parser.add_argument("input", nargs="+", help="input filenames")
19
+parser.add_argument("--figsize", type=float, nargs=2, default=(6., 4.), help="figure size")
20
+parser.add_argument("--block", type=int, default=4, help="block level")
21
+args = parser.parse_args()
22
+
23
+mpl.rc("pgf", texsystem="xelatex", preamble=(r"\usepackage{amsmath,bm,txfonts}",))
24
+mpl.rc("font", family="serif", serif=(), size=11)
25
+mpl.rc("legend", fontsize=11, numpoints=1)
26
+
27
+with h5.File(args.input[0], "r") as f:
28
+  g = f["dynamics/dimer/orientation_autocorrelation"]
29
+  time = g["time"][:]
30
+  index = g["count"][:] != 0
31
+  mean = np.zeros(shape=g["value"].shape)
32
+  error = np.zeros(shape=g["value"].shape)
33
+
34
+count = 0
35
+for fn in args.input:
36
+  with h5.File(fn, "r") as f:
37
+    g = f["dynamics/dimer/orientation_autocorrelation"]
38
+    assert np.all(time == g["time"])
39
+    value = g["value"][:]
40
+    count += 1
41
+    delta = value - mean
42
+    mean += delta / count
43
+    error += delta * (value - mean)
44
+
45
+error = np.sqrt(error / (count - 1) / count)
46
+
47
+time = time[args.block]
48
+index = index[args.block]
49
+mean = mean[args.block]
50
+error = error[args.block]
51
+
52
+fig = plt.figure(figsize=args.figsize)
53
+ax = fig.add_subplot(1, 1, 1)
54
+ax.errorbar(time[index[:, 0]], mean[index[:, 0], 0], error[index[:, 0], 0], errorevery=10, label="x")
55
+ax.errorbar(time[index[:, 1]], mean[index[:, 1], 1], error[index[:, 1], 1], errorevery=10, label="y")
56
+ax.errorbar(time[index[:, 2]], mean[index[:, 2], 2], error[index[:, 2], 2], errorevery=10, label="z")
57
+ax.minorticks_on()
58
+ax.set_xlabel(r"$t$")
59
+ax.set_ylabel(r"$\langle\hat{\bm{z}}(t)\cdot\hat{\bm{z}}_0\rangle$")
60
+ax.legend(loc="best", fancybox=True)
61
+fig.tight_layout(pad=0.1)
62
+fig.savefig(args.output)

+ 65
- 0
examples/confined_dimer/tools/dimer_velocity_autocorrelation.py View File

@@ -0,0 +1,65 @@
1
+#!/usr/bin/env python
2
+# -*- coding: utf-8 -*-
3
+#
4
+# Plot velocity autocorrelation function of dimer.
5
+# Copyright © 2013–2015 Peter Colberg.
6
+# Distributed under the MIT license. (See accompanying file LICENSE.)
7
+#
8
+
9
+import matplotlib as mpl
10
+mpl.use("pgf")
11
+import matplotlib.pyplot as plt
12
+import numpy as np
13
+import h5py as h5
14
+import argparse
15
+
16
+parser = argparse.ArgumentParser()
17
+parser.add_argument("output", help="plot filename")
18
+parser.add_argument("input", nargs="+", help="input filenames")
19
+parser.add_argument("--figsize", type=float, nargs=2, default=(6., 4.), help="figure size")
20
+parser.add_argument("--block", type=int, default=1, help="block level")
21
+args = parser.parse_args()
22
+
23
+mpl.rc("pgf", texsystem="xelatex", preamble=(r"\usepackage{amsmath,bm,txfonts}",))
24
+mpl.rc("font", family="serif", serif=(), size=11)
25
+mpl.rc("legend", fontsize=11, numpoints=1)
26
+
27
+with h5.File(args.input[0], "r") as f:
28
+  g = f["dynamics/dimer/velocity_autocorrelation"]
29
+  time = g["time"][:]
30
+  index = g["count"][:] != 0
31
+  mean = np.zeros(shape=g["value"].shape)
32
+  error = np.zeros(shape=g["value"].shape)
33
+
34
+count = 0
35
+for fn in args.input:
36
+  with h5.File(fn, "r") as f:
37
+    g = f["dynamics/dimer/velocity_autocorrelation"]
38
+    assert np.all(time == g["time"])
39
+    value = g["value"][:]
40
+    count += 1
41
+    delta = value - mean
42
+    mean += delta / count
43
+    error += delta * (value - mean)
44
+
45
+error = np.sqrt(error / (count - 1) / count)
46
+
47
+time = time[args.block]
48
+index = index[args.block]
49
+mean = mean[args.block]
50
+error = error[args.block]
51
+
52
+mean[:] /= np.sum(mean[0, :])
53
+error[:] /= np.sum(mean[0, :])
54
+
55
+fig = plt.figure(figsize=args.figsize)
56
+ax = fig.add_subplot(1, 1, 1)
57
+ax.errorbar(time[index[:, 0]], mean[index[:, 0], 0], error[index[:, 0], 0], errorevery=10, label="x")
58
+ax.errorbar(time[index[:, 1]], mean[index[:, 1], 1], error[index[:, 1], 1], errorevery=10, label="y")
59
+ax.errorbar(time[index[:, 2]], mean[index[:, 2], 2], error[index[:, 2], 2], errorevery=10, label="z")
60
+ax.minorticks_on()
61
+ax.set_xlabel(r"$t$")
62
+ax.set_ylabel(r"$\langle\hat{\bm{V}}(t)\cdot\hat{\bm{V}}_0\rangle$")
63
+ax.legend(loc="best", fancybox=True)
64
+fig.tight_layout(pad=0.1)
65
+fig.savefig(args.output)

+ 52
- 0
examples/confined_dimer/tools/wall_density.py View File

@@ -0,0 +1,52 @@
1
+#!/usr/bin/env python
2
+# -*- coding: utf-8 -*-
3
+#
4
+# Plot wall density function.
5
+# Copyright © 2013–2015 Peter Colberg.
6
+# Distributed under the MIT license.
7
+#
8
+
9
+import matplotlib as mpl
10
+mpl.use("pgf")
11
+import matplotlib.pyplot as plt
12
+import numpy as np
13
+import h5py as h5
14
+import argparse
15
+
16
+parser = argparse.ArgumentParser()
17
+parser.add_argument("output", help="plot filename")
18
+parser.add_argument("input", nargs="+", help="input filenames")
19
+parser.add_argument("--figsize", type=float, nargs=2, default=(6., 3.), help="figure size")
20
+args = parser.parse_args()
21
+
22
+mpl.rc("pgf", texsystem="xelatex", preamble=(r"\usepackage{amsmath,bm,txfonts}",))
23
+mpl.rc("font", family="serif", serif=(), size=11)
24
+mpl.rc("legend", fontsize=11, numpoints=1)
25
+
26
+with h5.File(args.input[0], "r") as f:
27
+  g = f["structure/solvent/wall_density"]
28
+  Z = g["distance"][:]
29
+  value = np.zeros(shape=g["value"].shape)
30
+  region = g.attrs["region"]
31
+  L = f["particles/dimer/box/edges"][:]
32
+
33
+count = 0
34
+for fn in args.input:
35
+  with h5.File(fn, "r") as f:
36
+    g = f["structure/solvent/wall_density"]
37
+    assert np.all(Z == g["distance"])
38
+    count += 1
39
+    value += (g["value"][:] - value) / count
40
+
41
+fig = plt.figure(figsize=args.figsize)
42
+ax = fig.add_subplot(1, 1, 1)
43
+ax.plot(Z, value[:, 0], label=region[0])
44
+ax.plot(Z, value[:, 1], label=region[1], ls=":")
45
+ax.minorticks_on()
46
+ax.set_xlim(0, L[2])
47
+ax.set_ylim(bottom=0)
48
+ax.set_xlabel(r"$Z$")
49
+ax.set_ylabel(r"$\varrho_s(Z)$")
50
+ax.legend(loc="best", fancybox=True)
51
+fig.tight_layout(pad=0.1)
52
+fig.savefig(args.output)

+ 1
- 0
examples/env.sh View File

@@ -0,0 +1 @@
1
+export LUA_PATH="$PWD/?.lua;${LUA_PATH:-;}"

+ 64
- 0
examples/many_dimer/equilibration/config.lua View File

@@ -0,0 +1,64 @@
1
+-- Edge lengths of simulation domain.
2
+L = {60, 60, 30}
3
+-- Periodic boundary conditions.
4
+periodic = {true, true, false}
5
+-- Lennard-Jones wall potential.
6
+wall = {dimer = {epsilon = 5.0, sigma = {C = L[3]/2, N = L[3]/2}}, solvent = {epsilon = 1.0, r_cut = 3^(1/6), sigma = 1.0}}
7
+-- Number of dimer spheres.
8
+Nd = 20
9
+-- Diameters of dimer spheres.
10
+diameter = {C = 4.0, N = 8.0}
11
+-- Solvent number density.
12
+rho = 0.8
13
+-- Masses of neutrally buoyant dimer spheres.
14
+mass = {C = rho * (math.pi/6) * diameter.C^3, N = rho * (math.pi/6) * diameter.N^3}
15
+-- Number of solvent particles excluding volume of walls and dimer spheres.
16
+Ns = math.floor(rho * L[1] * L[2] * (L[3] - 1.0) - (Nd / 2) * (mass.C + mass.N))
17
+-- Cutoff radius of potential in units of σ.
18
+r_cut = 2^(1/6)
19
+-- Minimum number of steps between solvent particle sorts.
20
+sort = {interval = 100}
21
+-- Distance of dimer spheres.
22
+bond = ((diameter.C + diameter.N)/2 + 1) * r_cut
23
+-- Potential well depths of dimer spheres with solvent particles.
24
+epsilon = {CA = 1.0, CB = 1.0, NA = 1.0, NB = 0.1, CC = 5.0, CN = 5.0, NN = 5.0}
25
+-- Potential diameters.
26
+sigma = {C = (diameter.C + 1)/2, N = (diameter.N + 1)/2, CC = diameter.C + 3, CN = (diameter.C + diameter.N)/2 + 3, NN = diameter.N + 3}
27
+-- Forward and backward reaction of solvent particles.
28
+reaction = {rate = 0.001}
29
+-- Verlet neighbour list skin.
30
+skin = 0.5
31
+-- Number of placeholders per neighbour list.
32
+nlist_size = {ss = 25, ds = 1000}
33
+-- Number of bins per dimension for solvent particles.
34
+nbin = {L[1]/2, L[2]/2, L[3]/2}
35
+-- Number of placeholders per bin.
36
+bin_size = 20
37
+-- Integration time-step.
38
+timestep = 0.001
39
+-- Number of steps.
40
+nsteps = 100000
41
+-- Thermodynamic variables.
42
+thermo = {interval = 100, datatype = "double"}
43
+-- Dimer trajectory.
44
+trajd = {interval = 100, datatype = "double", species = {"C", "N"}}
45
+-- Solvent trajectory.
46
+trajs = {interval = nsteps, datatype = "double", species = {"A", "B"}}
47
+-- Count solvent species.
48
+species = {interval = 1000, species = {"A", "B"}}
49
+-- OpenCL platforms.
50
+platform = 1
51
+-- OpenCL device per platform.
52
+device = 1
53
+-- Random number generator seed.
54
+seed = "0xd946f32b8853b7a1"
55
+-- H5MD input filename.
56
+input = "../thermostat/many_dimer.h5"
57
+-- Read trajectory at step.
58
+step = 100000
59
+-- H5MD output filename.
60
+output = "many_dimer.h5"
61
+-- H5MD file author.
62
+author = {name = "B.J. Alder", email = "bjalder@example.org"}
63
+-- H5MD file creator.
64
+creator = {name = "many_dimer.lua", version = "1.0"}

+ 152
- 0
examples/many_dimer/equilibration/many_dimer.lua View File

@@ -0,0 +1,152 @@
1
+#!/usr/bin/env luajit
2
+------------------------------------------------------------------------------
3
+-- Equilibration phase of many-dimer simulation.
4
+-- Copyright © 2013–2015 Peter Colberg.
5
+-- Distributed under the MIT license. (See accompanying file LICENSE.)
6
+------------------------------------------------------------------------------
7
+
8
+-- Check uses of undeclared global variables.
9
+require("strict")
10
+
11
+local nm = {
12
+  config = require("nanomotor.config"),
13
+  box = require("nanomotor.box"),
14
+  domain = require("nanomotor.domain"),
15
+  integrate = require("nanomotor.integrate"),
16
+  observe = require("nanomotor.observe"),
17
+  h5md = require("nanomotor.h5md"),
18
+  trajd = require("nanomotor.trajd"),
19
+  trajs = require("nanomotor.trajs"),
20
+  thermo = require("nanomotor.thermo"),
21
+  species = require("nanomotor.species"),
22
+}
23
+
24
+local cl = require("opencl")
25
+local hdf5 = require("hdf5")
26
+local ffi = require("ffi")
27
+local syscall = require("syscall")
28
+local json = require("cjson")
29
+
30
+-- Read configuration file.
31
+local args = nm.config(arg[1])
32
+
33
+-- Seed random number generator.
34
+math.randomseed(args.seed)
35
+
36
+-- Create OpenCL context on chosen OpenCL device.
37
+local platform = cl.get_platforms()[args.platform]
38
+print(platform:get_info("version"))
39
+print(platform:get_info("name"))
40
+print(platform:get_info("vendor"))
41
+print(platform:get_info("extensions"))
42
+local device = platform:get_devices()[args.device]
43
+print(device:get_info("name"))
44
+print(device:get_info("vendor"))
45
+print(device:get_info("extensions"))
46
+local context = cl.create_context({device})
47
+local queue = context:create_command_queue(device)
48
+
49
+-- Create simulation box with given boundary conditions.
50
+local box = nm.box(args)
51
+
52
+-- Allocate zero-filled particle buffers.
53
+local dom = nm.domain(queue, box, args)
54
+
55
+-- Read particle coordinates from H5MD input file.
56
+local file = hdf5.open_file(args.input)
57
+do
58
+  local group = file:open_group("particles/dimer")
59
+  local obs_rd = nm.h5md.open_observable(group:open_group("position"))
60
+  local obs_vd = nm.h5md.open_observable(group:open_group("velocity"))
61
+  local space_rd = hdf5.create_simple_space({dom.Nd, 4})
62
+  local space_vd = hdf5.create_simple_space({dom.Nd, 4})
63
+  space_rd:select_hyperslab("set", {0, 0}, nil, {dom.Nd, 3})
64
+  space_vd:select_hyperslab("set", {0, 0}, nil, {dom.Nd, 3})
65
+  obs_rd:read(args.step, dom.rd, hdf5.double, space_rd)
66
+  obs_vd:read(args.step, dom.vd, hdf5.double, space_vd)
67
+  queue:enqueue_write_buffer(dom.d_rd, true, dom.rd)
68
+  queue:enqueue_write_buffer(dom.d_vd, true, dom.vd)
69
+end
70
+do
71
+  local group = file:open_group("particles/solvent")
72
+  local obs_rs = nm.h5md.open_observable(group:open_group("position"))
73
+  local obs_vs = nm.h5md.open_observable(group:open_group("velocity"))
74
+  local obs_species = nm.h5md.open_observable(group:open_group("species"))
75
+  local space_rs = hdf5.create_simple_space({dom.Ns, 4})
76
+  local space_vs = hdf5.create_simple_space({dom.Ns, 4})
77
+  local space_species = hdf5.create_simple_space({dom.Ns, 4})
78
+  space_rs:select_hyperslab("set", {0, 0}, nil, {dom.Ns, 3})
79
+  space_vs:select_hyperslab("set", {0, 0}, nil, {dom.Ns, 3})
80
+  space_species:select_hyperslab("set", {0, 0}, nil, {dom.Ns, 1})
81
+  local dtype_species = hdf5.int8:enum_create()
82
+  dtype_species:enum_insert("A", ffi.new("cl_char[1]", 0))
83
+  dtype_species:enum_insert("B", ffi.new("cl_char[1]", 1))
84
+  obs_rs:read(args.step, dom.rs, hdf5.double, space_rs)
85
+  obs_vs:read(args.step, dom.vs, hdf5.double, space_vs)
86
+  obs_species:read(args.step, dom.species, dtype_species, space_species)
87
+  queue:enqueue_write_buffer(dom.d_rs, true, dom.rs)
88
+  queue:enqueue_write_buffer(dom.d_vs, true, dom.vs)
89
+  queue:enqueue_write_buffer(dom.d_species, true, dom.species)
90
+end
91
+file:close()
92
+
93
+-- Create MD integrator.
94
+local integrate = nm.integrate(dom, args)
95
+
96
+-- Create H5MD output file with temporary filename.
97
+local file = hdf5.create_file(args.output..".tmp")
98
+-- Write H5MD metadata group.
99
+nm.h5md.write_meta(file, args)
100
+
101
+-- Serialize simulation parameters to JSON string attribute.
102
+do
103
+  local space = hdf5.create_space("scalar")
104
+  local dtype = hdf5.c_s1:copy()
105
+  local s = json.encode(args)
106
+  dtype:set_size(#s+1)
107
+  local attr = file:create_attribute("parameters", dtype, space)
108
+  attr:write(ffi.cast("const char *", s), dtype)
109
+end
110
+
111
+-- Create observables.
112
+local observables = {}
113
+if args.thermo then table.insert(observables, nm.thermo(integrate, file, args.thermo)) end
114
+if args.trajd then table.insert(observables, nm.trajd(integrate, file, args.trajd)) end
115
+if args.trajs then table.insert(observables, nm.trajs(integrate, file, args.trajs)) end
116
+if args.species then table.insert(observables, nm.species(integrate, file, args.species)) end
117
+
118
+-- Periodically log step, total energy, and remaining run-time.
119
+do
120
+  local interval = 1000
121
+  local self = {step = 0}
122
+
123
+  function self.observe(nsteps)
124
+    local t1 = syscall.clock_gettime("monotonic")
125
+    while self.step <= nsteps do
126
+      coroutine.yield()
127
+      local en = integrate.compute_en()
128
+      local t2 = syscall.clock_gettime("monotonic")
129
+      local t = (t2.time - t1.time) / interval * (nsteps - self.step)
130
+      local u = "s"
131
+      if t >= 60 then t, u = t/60, "m" end
132
+      if t >= 60 then t, u = t/60, "h" end
133
+      io.write("#", self.step, " ", en.en_tot, " ", ("%.1f"):format(t), u, "\n")
134
+      io.flush()
135
+      self.step = self.step + interval
136
+      t1 = t2
137
+    end
138
+  end
139
+
140
+  table.insert(observables, self)
141
+end
142
+
143
+-- Run simulation.
144
+nm.observe(integrate, observables, args.nsteps)
145
+
146
+-- Synchronize output files to disk.
147
+file:flush()
148
+assert((syscall.fsync(ffi.cast("int *", file:get_vfd_handle())[0])))
149
+file:close()
150
+
151
+-- Rename H5MD output file after successful simulation.
152
+assert((syscall.rename(args.output..".tmp", args.output)))

+ 257
- 0
examples/many_dimer/many_dimer.sh View File

@@ -0,0 +1,257 @@
1
+#PBS -q gpu
2
+#PBS -l nodes=1:ppn=2:gpus=2
3
+#PBS -l walltime=6:00:00
4
+#PBS -m n
5
+#PBS -V
6
+set -e -o pipefail
7
+PROJECT="${PBS_O_WORKDIR}/angstrom-dimer"
8
+cd "${PROJECT}" && . examples/env.sh
9
+OUTPUT="${OUTPUT:-${PBS_JOBNAME%%.sh*}_${PBS_JOBID%%[^0-9]*}}"
10
+SEED="$(date +%s)\0${PBS_O_HOST}\0${PBS_JOBID}"
11
+TASKS=()
12
+for TASK in 1 2
13
+do
14
+  if [ ! -e "${PBS_O_WORKDIR}/thermostat/${OUTPUT}_${TASK}.h5" ]
15
+  then
16
+    mkdir -p "${PBS_O_WORKDIR}/thermostat"
17
+    cat >"${PBS_O_WORKDIR}/thermostat/${OUTPUT}_${TASK}.lua" <<EOF
18
+-- Edge lengths of simulation domain.
19
+L = {60, 60, 30}
20
+-- Periodic boundary conditions.
21
+periodic = {true, true, false}
22
+-- Lennard-Jones wall potential.
23
+wall = {dimer = {epsilon = 5.0, sigma = {C = L[3]/2, N = L[3]/2}}, solvent = {epsilon = 1.0, r_cut = 3^(1/6), sigma = 1.0}}
24
+-- Number of dimer spheres.
25
+Nd = 20
26
+-- Diameters of dimer spheres.
27
+diameter = {C = 4.0, N = 8.0}
28
+-- Solvent number density.
29
+rho = 0.8
30
+-- Masses of neutrally buoyant dimer spheres.
31
+mass = {C = rho * (math.pi/6) * diameter.C^3, N = rho * (math.pi/6) * diameter.N^3}
32
+-- Number of solvent particles excluding volume of walls and dimer spheres.
33
+Ns = math.floor(rho * L[1] * L[2] * (L[3] - 1.0) - (Nd / 2) * (mass.C + mass.N))
34
+-- Cutoff radius of potential in units of σ.
35
+r_cut = 2^(1/6)
36
+-- Minimum number of steps between solvent particle sorts.
37
+sort = {interval = 100}
38
+-- Lattice constant and region for solvent particles.
39
+lattice = {a = 1.6, slab = {{0, 0, 1}, {L[1], L[2], L[3]-1}}}
40
+-- Solvent temperature.
41
+temp = 1
42
+-- Distance of dimer spheres.
43
+bond = ((diameter.C + diameter.N)/2 + 1) * r_cut
44
+-- Potential well depths of dimer spheres with solvent particles.
45
+epsilon = {CA = 1.0, CB = 1.0, NA = 1.0, NB = 0.1, CC = 5.0, CN = 5.0, NN = 5.0}
46
+-- Potential diameters.
47
+sigma = {C = (diameter.C + 1)/2, N = (diameter.N + 1)/2, CC = diameter.C + 3, CN = (diameter.C + diameter.N)/2 + 3, NN = diameter.N + 3}
48
+-- Verlet neighbour list skin.
49
+skin = 0.5
50
+-- Number of placeholders per neighbour list.
51
+nlist_size = {ss = 25, ds = 1000}
52
+-- Number of bins per dimension for solvent particles.
53
+nbin = {L[1]/2, L[2]/2, L[3]/2}
54
+-- Number of placeholders per bin.
55
+bin_size = 20
56
+-- Integration time-step.
57
+timestep = 0.001
58
+-- Number of steps.
59
+nsteps = 20000
60
+-- Thermodynamic variables.
61
+thermo = {interval = 1000, datatype = "double"}
62
+-- Dimer trajectory.
63
+trajd = {interval = 1000, datatype = "double", species = {"C", "N"}}
64
+-- Solvent trajectory.
65
+trajs = {interval = nsteps, datatype = "double", species = {"A", "B"}}
66
+-- OpenCL platforms.
67
+platform = 1
68
+-- OpenCL device per platform.
69
+device = ${TASK}
70
+-- Random number generator seed.
71
+seed = "0x$(echo -ne "${SEED}\0${TASK}\0thermostat" | md5sum | cut -c-16)"
72
+-- H5MD output filename.
73
+output = "${PBS_O_WORKDIR}/thermostat/${OUTPUT}_${TASK}.h5"
74
+-- H5MD file author.
75
+author = {name = "$(git config user.name)", email = "$(git config user.email)"}
76
+-- H5MD file creator.
77
+creator = {name = "many_dimer.lua", version = "$(cd "${PROJECT}" && git describe --always)"}
78
+EOF
79
+    cd "${PROJECT}/examples/many_dimer/thermostat"
80
+    ./many_dimer.lua "${PBS_O_WORKDIR}/thermostat/${OUTPUT}_${TASK}.lua" | awk '{ print strftime("%Y-%m-%d %H:%M:%S"), $0; fflush(); }' >>"${PBS_O_WORKDIR}/thermostat/${OUTPUT}_${TASK}.log" &
81
+    TASKS+=("${!}")
82
+  fi
83
+done
84
+for TASK in "${TASKS[@]}"
85
+do
86
+  wait "${TASK}"
87
+done
88
+TASKS=()
89
+for TASK in 1 2
90
+do
91
+  if [ ! -e "${PBS_O_WORKDIR}/equilibration/${OUTPUT}_${TASK}.h5" ]
92
+  then
93
+    mkdir -p "${PBS_O_WORKDIR}/equilibration"
94
+    cat >"${PBS_O_WORKDIR}/equilibration/${OUTPUT}_${TASK}.lua" <<EOF
95
+-- Edge lengths of simulation domain.
96
+L = {60, 60, 30}
97
+-- Periodic boundary conditions.
98
+periodic = {true, true, false}
99
+-- Lennard-Jones wall potential.
100
+wall = {dimer = {epsilon = 5.0, sigma = {C = L[3]/2, N = L[3]/2}}, solvent = {epsilon = 1.0, r_cut = 3^(1/6), sigma = 1.0}}
101
+-- Number of dimer spheres.
102
+Nd = 20
103
+-- Diameters of dimer spheres.
104
+diameter = {C = 4.0, N = 8.0}
105
+-- Solvent number density.
106
+rho = 0.8
107
+-- Masses of neutrally buoyant dimer spheres.
108
+mass = {C = rho * (math.pi/6) * diameter.C^3, N = rho * (math.pi/6) * diameter.N^3}
109
+-- Number of solvent particles excluding volume of walls and dimer spheres.
110
+Ns = math.floor(rho * L[1] * L[2] * (L[3] - 1.0) - (Nd / 2) * (mass.C + mass.N))
111
+-- Cutoff radius of potential in units of σ.
112
+r_cut = 2^(1/6)
113
+-- Minimum number of steps between solvent particle sorts.
114
+sort = {interval = 100}
115
+-- Distance of dimer spheres.
116
+bond = ((diameter.C + diameter.N)/2 + 1) * r_cut
117
+-- Potential well depths of dimer spheres with solvent particles.
118
+epsilon = {CA = 1.0, CB = 1.0, NA = 1.0, NB = 0.1, CC = 5.0, CN = 5.0, NN = 5.0}
119
+-- Potential diameters.
120
+sigma = {C = (diameter.C + 1)/2, N = (diameter.N + 1)/2, CC = diameter.C + 3, CN = (diameter.C + diameter.N)/2 + 3, NN = diameter.N + 3}
121
+-- Forward and backward reaction of solvent particles.
122
+reaction = {rate = 0.001}
123
+-- Verlet neighbour list skin.
124
+skin = 0.5
125
+-- Number of placeholders per neighbour list.
126
+nlist_size = {ss = 25, ds = 1000}
127
+-- Number of bins per dimension for solvent particles.
128
+nbin = {L[1]/2, L[2]/2, L[3]/2}
129
+-- Number of placeholders per bin.
130
+bin_size = 20
131
+-- Integration time-step.
132
+timestep = 0.001
133
+-- Number of steps.
134
+nsteps = 2000000
135
+-- Thermodynamic variables.
136
+thermo = {interval = 1000, datatype = "double"}
137
+-- Dimer trajectory.
138
+trajd = {interval = 1000, datatype = "double", species = {"C", "N"}}
139
+-- Solvent trajectory.
140
+trajs = {interval = nsteps, datatype = "double", species = {"A", "B"}}
141
+-- Count solvent species.
142
+species = {interval = 1000, species = {"A", "B"}}
143
+-- OpenCL platforms.
144
+platform = 1
145
+-- OpenCL device per platform.
146
+device = ${TASK}
147
+-- Random number generator seed.
148
+seed = "0x$(echo -ne "${SEED}\0${TASK}\0equilibration" | md5sum | cut -c-16)"
149
+-- H5MD input filename.
150
+input = "${PBS_O_WORKDIR}/thermostat/${OUTPUT}_${TASK}.h5"
151
+-- Read trajectory at step.
152
+step = 20000
153
+-- H5MD output filename.
154
+output = "${PBS_O_WORKDIR}/equilibration/${OUTPUT}_${TASK}.h5"
155
+-- H5MD file author.
156
+author = {name = "$(git config user.name)", email = "$(git config user.email)"}
157
+-- H5MD file creator.
158
+creator = {name = "many_dimer.lua", version = "$(cd "${PROJECT}" && git describe --always)"}
159
+EOF
160
+    cd "${PROJECT}/examples/many_dimer/equilibration"
161
+    ./many_dimer.lua "${PBS_O_WORKDIR}/equilibration/${OUTPUT}_${TASK}.lua" | awk '{ print strftime("%Y-%m-%d %H:%M:%S"), $0; fflush(); }' >>"${PBS_O_WORKDIR}/equilibration/${OUTPUT}_${TASK}.log" &
162
+    TASKS+=("${!}")
163
+  fi
164
+done
165
+for TASK in "${TASKS[@]}"
166
+do
167
+  wait "${TASK}"
168
+done
169
+TASKS=()
170
+for TASK in 1 2
171
+do
172
+  if [ ! -e "${PBS_O_WORKDIR}/production/${OUTPUT}_${TASK}.h5" ]
173
+  then
174
+    mkdir -p "${PBS_O_WORKDIR}/production"
175
+    cat >"${PBS_O_WORKDIR}/production/${OUTPUT}_${TASK}.lua" <<EOF
176
+-- Edge lengths of simulation domain.
177
+L = {60, 60, 30}
178
+-- Periodic boundary conditions.
179
+periodic = {true, true, false}
180
+-- Lennard-Jones wall potential.
181
+wall = {dimer = {epsilon = 5.0, sigma = {C = L[3]/2, N = L[3]/2}}, solvent = {epsilon = 1.0, r_cut = 3^(1/6), sigma = 1.0}}
182
+-- Number of dimer spheres.
183
+Nd = 20
184
+-- Diameters of dimer spheres.
185
+diameter = {C = 4.0, N = 8.0}
186
+-- Solvent number density.
187
+rho = 0.8
188
+-- Masses of neutrally buoyant dimer spheres.
189
+mass = {C = rho * (math.pi/6) * diameter.C^3, N = rho * (math.pi/6) * diameter.N^3}
190
+-- Number of solvent particles excluding volume of walls and dimer spheres.
191
+Ns = math.floor(rho * L[1] * L[2] * (L[3] - 1.0) - (Nd / 2) * (mass.C + mass.N))
192
+-- Cutoff radius of potential in units of σ.
193
+r_cut = 2^(1/6)
194
+-- Minimum number of steps between solvent particle sorts.
195
+sort = {interval = 100}
196
+-- Distance of dimer spheres.
197
+bond = ((diameter.C + diameter.N)/2 + 1) * r_cut
198
+-- Potential well depths of dimer spheres with solvent particles.
199
+epsilon = {CA = 1.0, CB = 1.0, NA = 1.0, NB = 0.1, CC = 5.0, CN = 5.0, NN = 5.0}
200
+-- Potential diameters.
201
+sigma = {C = (diameter.C + 1)/2, N = (diameter.N + 1)/2, CC = diameter.C + 3, CN = (diameter.C + diameter.N)/2 + 3, NN = diameter.N + 3}
202
+-- Forward and backward reaction of solvent particles.
203
+reaction = {rate = 0.001}
204
+-- Verlet neighbour list skin.
205
+skin = 0.5
206
+-- Number of placeholders per neighbour list.
207
+nlist_size = {ss = 25, ds = 1000}
208
+-- Number of bins per dimension for solvent particles.
209
+nbin = {L[1]/2, L[2]/2, L[3]/2}
210
+-- Number of placeholders per bin.
211
+bin_size = 20
212
+-- Integration time-step.
213
+timestep = 0.001
214
+-- Number of steps.
215
+nsteps = 10000000
216
+-- Thermodynamic variables.
217
+thermo = {interval = 10000, datatype = "double"}
218
+-- Dimer trajectory.
219
+trajd = {interval = 1000, datatype = "float", species = {"C", "N"}}
220
+-- Solvent trajectory.
221
+trajs = {interval = nsteps, datatype = "float", species = {"A", "B"}}
222
+-- Count solvent species.
223
+species = {interval = 10000, species = {"A", "B"}}
224
+-- Mean-square displacement of dimer.
225
+msdd = {interval = 10, blockcount = 6, blocksize = 100, blockshift = 10, separation = 10000}
226
+-- Velocity autocorrelation function of dimer.
227
+vafd = {interval = 10, blockcount = 6, blocksize = 100, blockshift = 10, separation = 10000}
228
+-- Orientation autocorrelation of dimer.
229
+orient = {interval = 10, blockcount = 6, blocksize = 100, blockshift = 10, separation = 10000}
230
+-- OpenCL platforms.
231
+platform = 1
232
+-- OpenCL device per platform.
233
+device = ${TASK}
234
+-- Random number generator seed.
235
+seed = "0x$(echo -ne "${SEED}\0${TASK}\0production" | md5sum | cut -c-16)"
236
+-- H5MD input filename.
237
+input = "${PBS_O_WORKDIR}/equilibration/${OUTPUT}_${TASK}.h5"
238
+-- Read trajectory at step.
239
+step = 2000000
240
+-- H5MD output filename.
241
+output = "${PBS_O_WORKDIR}/production/${OUTPUT}_${TASK}.h5"
242
+-- Snapshot program state.
243
+snapshot = {interval = 1000000, output = "${PBS_O_WORKDIR}/production/${OUTPUT}_${TASK}-snapshot.h5"}
244
+-- H5MD file author.
245
+author = {name = "$(git config user.name)", email = "$(git config user.email)"}
246
+-- H5MD file creator.
247
+creator = {name = "many_dimer.lua", version = "$(cd "${PROJECT}" && git describe --always)"}
248
+EOF
249
+    cd "${PROJECT}/examples/many_dimer/production"
250
+    ./many_dimer.lua "${PBS_O_WORKDIR}/production/${OUTPUT}_${TASK}.lua" | awk '{ print strftime("%Y-%m-%d %H:%M:%S"), $0; fflush(); }' >>"${PBS_O_WORKDIR}/production/${OUTPUT}_${TASK}.log" &
251
+    TASKS+=("${!}")
252
+  fi
253
+done
254
+for TASK in "${TASKS[@]}"
255
+do
256
+  wait "${TASK}"
257
+done

+ 72
- 0
examples/many_dimer/production/config.lua View File

@@ -0,0 +1,72 @@
1
+-- Edge lengths of simulation domain.
2
+L = {60, 60, 30}
3
+-- Periodic boundary conditions.
4
+periodic = {true, true, false}
5
+-- Lennard-Jones wall potential.
6
+wall = {dimer = {epsilon = 5.0, sigma = {C = L[3]/2, N = L[3]/2}}, solvent = {epsilon = 1.0, r_cut = 3^(1/6), sigma = 1.0}}
7
+-- Number of dimer spheres.
8
+Nd = 20
9
+-- Diameters of dimer spheres.
10
+diameter = {C = 4.0, N = 8.0}
11
+-- Solvent number density.
12
+rho = 0.8
13
+-- Masses of neutrally buoyant dimer spheres.
14
+mass = {C = rho * (math.pi/6) * diameter.C^3, N = rho * (math.pi/6) * diameter.N^3}
15
+-- Number of solvent particles excluding volume of walls and dimer spheres.
16
+Ns = math.floor(rho * L[1] * L[2] * (L[3] - 1.0) - (Nd / 2) * (mass.C + mass.N))
17
+-- Cutoff radius of potential in units of σ.
18
+r_cut = 2^(1/6)
19
+-- Minimum number of steps between solvent particle sorts.
20
+sort = {interval = 100}
21
+-- Distance of dimer spheres.
22
+bond = ((diameter.C + diameter.N)/2 + 1) * r_cut
23
+-- Potential well depths of dimer spheres with solvent particles.
24
+epsilon = {CA = 1.0, CB = 1.0, NA = 1.0, NB = 0.1, CC = 5.0, CN = 5.0, NN = 5.0}
25
+-- Potential diameters.
26
+sigma = {C = (diameter.C + 1)/2, N = (diameter.N + 1)/2, CC = diameter.C + 3, CN = (diameter.C + diameter.N)/2 + 3, NN = diameter.N + 3}
27
+-- Forward and backward reaction of solvent particles.
28
+reaction = {rate = 0.001}
29
+-- Verlet neighbour list skin.
30
+skin = 0.5
31
+-- Number of placeholders per neighbour list.
32
+nlist_size = {ss = 25, ds = 1000}
33
+-- Number of bins per dimension for solvent particles.
34
+nbin = {L[1]/2, L[2]/2, L[3]/2}
35
+-- Number of placeholders per bin.
36
+bin_size = 20
37
+-- Integration time-step.
38
+timestep = 0.001
39
+-- Number of steps.
40
+nsteps = 100000
41
+-- Thermodynamic variables.
42
+thermo = {interval = 100, datatype = "double"}
43
+-- Dimer trajectory.
44
+trajd = {interval = 100, datatype = "float", species = {"C", "N"}}
45
+-- Solvent trajectory.
46
+trajs = {interval = 10000, datatype = "float", species = {"A", "B"}}
47
+-- Count solvent species.
48
+species = {interval = 1000, species = {"A", "B"}}
49
+-- Mean-square displacement of dimer.
50
+msdd = {interval = 10, blockcount = 6, blocksize = 100, blockshift = 10, separation = 1000}
51
+-- Velocity autocorrelation function of dimer.
52
+vafd = {interval = 10, blockcount = 6, blocksize = 100, blockshift = 10, separation = 1000}
53
+-- Orientation autocorrelation of dimer.
54
+orient = {interval = 10, blockcount = 6, blocksize = 100, blockshift = 10, separation = 1000}
55
+-- OpenCL platforms.
56
+platform = 1
57
+-- OpenCL device per platform.
58
+device = 1
59
+-- Random number generator seed.
60
+seed = "0xc0ede07b69816a9e"
61
+-- H5MD input filename.
62
+input = "../equilibration/many_dimer.h5"
63
+-- Read trajectory at step.
64
+step = 100000
65
+-- H5MD output filename.
66
+output = "many_dimer.h5"
67
+-- Snapshot program state.
68
+snapshot = {interval = 10000, output = "many_dimer-snapshot.h5"}
69
+-- H5MD file author.
70
+author = {name = "B.J. Alder", email = "bjalder@example.org"}
71
+-- H5MD file creator.
72
+creator = {name = "many_dimer.lua", version = "1.0"}

+ 163
- 0
examples/many_dimer/production/many_dimer.lua View File

@@ -0,0 +1,163 @@
1
+#!/usr/bin/env luajit
2
+------------------------------------------------------------------------------
3
+-- Production phase of many-dimer simulation.
4
+-- Copyright © 2013–2015 Peter Colberg.
5
+-- Distributed under the MIT license. (See accompanying file LICENSE.)
6
+------------------------------------------------------------------------------
7
+
8
+-- Check uses of undeclared global variables.
9
+require("strict")
10
+
11
+local nm = {
12
+  config = require("nanomotor.config"),
13
+  box = require("nanomotor.box"),
14
+  domain = require("nanomotor.domain"),
15
+  integrate = require("nanomotor.integrate"),
16
+  observe = require("nanomotor.observe"),
17
+  snapshot = require("nanomotor.snapshot"),
18
+  h5md = require("nanomotor.h5md"),
19
+  trajd = require("nanomotor.trajd"),
20
+  trajs = require("nanomotor.trajs"),
21
+  thermo = require("nanomotor.thermo"),
22
+  species = require("nanomotor.species"),
23
+  corr = require("nanomotor.correlation"),
24
+  msdd = require("msdd"),
25
+  vafd = require("vafd"),
26
+  orient = require("orient"),
27
+}
28
+
29
+local cl = require("opencl")
30
+local hdf5 = require("hdf5")
31
+local ffi = require("ffi")
32
+local syscall = require("syscall")
33
+local json = require("cjson")
34
+
35
+-- Read configuration file.
36
+local args = nm.config(arg[1])
37
+
38
+-- Seed random number generator.
39
+math.randomseed(args.seed)
40
+
41
+-- Create OpenCL context on chosen OpenCL device.
42
+local platform = cl.get_platforms()[args.platform]
43
+print(platform:get_info("version"))
44
+print(platform:get_info("name"))
45
+print(platform:get_info("vendor"))
46
+print(platform:get_info("extensions"))
47
+local device = platform:get_devices()[args.device]
48
+print(device:get_info("name"))
49
+print(device:get_info("vendor"))
50
+print(device:get_info("extensions"))
51
+local context = cl.create_context({device})
52
+local queue = context:create_command_queue(device)
53
+
54
+-- Create simulation box with given boundary conditions.
55
+local box = nm.box(args)
56
+
57
+-- Allocate zero-filled particle buffers.
58
+local dom = nm.domain(queue, box, args)
59
+
60
+-- Read particle coordinates from H5MD input file.
61
+local file = hdf5.open_file(args.input)
62
+do
63
+  local group = file:open_group("particles/dimer")
64
+  local obs_rd = nm.h5md.open_observable(group:open_group("position"))
65
+  local obs_vd = nm.h5md.open_observable(group:open_group("velocity"))
66
+  local space_rd = hdf5.create_simple_space({dom.Nd, 4})
67
+  local space_vd = hdf5.create_simple_space({dom.Nd, 4})
68
+  space_rd:select_hyperslab("set", {0, 0}, nil, {dom.Nd, 3})
69
+  space_vd:select_hyperslab("set", {0, 0}, nil, {dom.Nd, 3})
70
+  obs_rd:read(args.step, dom.rd, hdf5.double, space_rd)
71
+  obs_vd:read(args.step, dom.vd, hdf5.double, space_vd)
72
+  queue:enqueue_write_buffer(dom.d_rd, true, dom.rd)
73
+  queue:enqueue_write_buffer(dom.d_vd, true, dom.vd)
74
+end
75
+do
76
+  local group = file:open_group("particles/solvent")
77
+  local obs_rs = nm.h5md.open_observable(group:open_group("position"))
78
+  local obs_vs = nm.h5md.open_observable(group:open_group("velocity"))
79
+  local obs_species = nm.h5md.open_observable(group:open_group("species"))
80
+  local space_rs = hdf5.create_simple_space({dom.Ns, 4})
81
+  local space_vs = hdf5.create_simple_space({dom.Ns, 4})
82
+  local space_species = hdf5.create_simple_space({dom.Ns, 4})
83
+  space_rs:select_hyperslab("set", {0, 0}, nil, {dom.Ns, 3})
84
+  space_vs:select_hyperslab("set", {0, 0}, nil, {dom.Ns, 3})
85
+  space_species:select_hyperslab("set", {0, 0}, nil, {dom.Ns, 1})
86
+  local dtype_species = hdf5.int8:enum_create()
87
+  dtype_species:enum_insert("A", ffi.new("cl_char[1]", 0))
88
+  dtype_species:enum_insert("B", ffi.new("cl_char[1]", 1))
89
+  obs_rs:read(args.step, dom.rs, hdf5.double, space_rs)
90
+  obs_vs:read(args.step, dom.vs, hdf5.double, space_vs)
91
+  obs_species:read(args.step, dom.species, dtype_species, space_species)
92
+  queue:enqueue_write_buffer(dom.d_rs, true, dom.rs)
93
+  queue:enqueue_write_buffer(dom.d_vs, true, dom.vs)
94
+  queue:enqueue_write_buffer(dom.d_species, true, dom.species)
95
+end
96
+file:close()
97
+
98
+-- Create MD integrator.
99
+local integrate = nm.integrate(dom, args)
100
+
101
+-- Create H5MD output file with temporary filename.
102
+local file = hdf5.create_file(args.output..".tmp")
103
+-- Write H5MD metadata group.
104
+nm.h5md.write_meta(file, args)
105
+
106
+-- Serialize simulation parameters to JSON string attribute.
107
+do
108
+  local space = hdf5.create_space("scalar")
109
+  local dtype = hdf5.c_s1:copy()
110
+  local s = json.encode(args)
111
+  dtype:set_size(#s+1)
112
+  local attr = file:create_attribute("parameters", dtype, space)
113
+  attr:write(ffi.cast("const char *", s), dtype)
114
+end
115
+
116
+-- Create observables.
117
+local observables = {}
118
+if args.snapshot then table.insert(observables, nm.snapshot(dom, observables, args.snapshot)) end
119
+if args.thermo then table.insert(observables, nm.thermo(integrate, file, args.thermo)) end
120
+if args.trajd then table.insert(observables, nm.trajd(integrate, file, args.trajd)) end
121
+if args.trajs then table.insert(observables, nm.trajs(integrate, file, args.trajs)) end
122
+if args.species then table.insert(observables, nm.species(integrate, file, args.species)) end
123
+if args.msdd then table.insert(observables, nm.corr(integrate, nm.msdd(dom, file), args.msdd)) end
124
+if args.vafd then table.insert(observables, nm.corr(integrate, nm.vafd(dom, file), args.vafd)) end
125
+if args.orient then table.insert(observables, nm.corr(integrate, nm.orient(dom, file), args.orient)) end
126
+
127
+-- Periodically log step, total energy, and remaining run-time.
128
+do
129
+  local interval = 10000
130
+  local self = {step = 0}
131
+
132
+  function self.observe(nsteps)
133
+    local t1 = syscall.clock_gettime("monotonic")
134
+    while self.step <= nsteps do
135
+      coroutine.yield()
136
+      local en = integrate.compute_en()
137
+      local t2 = syscall.clock_gettime("monotonic")
138
+      local t = (t2.time - t1.time) / interval * (nsteps - self.step)
139
+      local u = "s"
140
+      if t >= 60 then t, u = t/60, "m" end
141
+      if t >= 60 then t, u = t/60, "h" end
142
+      io.write("#", self.step, " ", en.en_tot, " ", ("%.1f"):format(t), u, "\n")
143
+      io.flush()
144
+      self.step = self.step + interval
145
+      t1 = t2
146
+    end
147
+  end
148
+
149
+  table.insert(observables, self)
150
+end
151
+
152
+-- Run simulation.
153
+nm.observe(integrate, observables, args.nsteps)
154
+
155
+-- Synchronize output files to disk.
156
+file:flush()
157
+assert((syscall.fsync(ffi.cast("int *", file:get_vfd_handle())[0])))
158
+file:close()
159
+
160
+-- Rename H5MD output file after successful simulation.
161
+assert((syscall.rename(args.output..".tmp", args.output)))
162
+-- Remove snapshot file.
163
+if args.snapshot then assert((syscall.unlink(args.snapshot.output))) end

+ 69
- 0
examples/many_dimer/production/msdd.lua View File

@@ -0,0 +1,69 @@
1
+------------------------------------------------------------------------------
2
+-- Mean-square displacement of dimer.
3
+-- Copyright © 2013–2015 Peter Colberg.
4
+-- Distributed under the MIT license. (See accompanying file LICENSE.)
5
+------------------------------------------------------------------------------
6
+
7
+local hdf5 = require("hdf5")
8
+local ffi = require("ffi")
9
+
10
+-- Cache C types.
11
+local vector_n = ffi.typeof("struct { double x, y, z; }[?]")
12
+
13
+return function(dom, file)
14
+  local box = dom.box
15
+
16
+  local self = {
17
+    group = file:create_group("dynamics/dimer/mean_square_displacement"),
18
+    dims = {3},
19
+  }
20
+
21
+  function self.sample()
22
+    local r_cm = vector_n(dom.Nd/2)
23
+    dom.queue:enqueue_read_buffer(dom.d_rd, true, dom.rd)
24
+    dom.queue:enqueue_read_buffer(dom.d_imd, true, dom.imd)
25
+    local rd, imd, L = dom.rd, dom.imd, box.L
26
+    local m_N = 1 / (1 + dom.mass.C/dom.mass.N)
27
+    for i = 0, dom.Nd/2-1 do
28
+      local dx = rd[2*i+1].x - rd[2*i].x
29
+      local dy = rd[2*i+1].y - rd[2*i].y
30
+      local dz = rd[2*i+1].z - rd[2*i].z
31
+      dx, dy, dz = box.mindist(dx, dy, dz)
32
+      r_cm[i].x = rd[2*i].x + L[1]*imd[2*i].x + m_N*dx
33
+      r_cm[i].y = rd[2*i].y + L[2]*imd[2*i].y + m_N*dy
34
+      r_cm[i].z = rd[2*i].z + L[3]*imd[2*i].z + m_N*dz
35
+    end
36
+    return r_cm
37
+  end
38
+
39
+  function self.correlate(r, r0, result)
40
+    local x, y, z = 0, 0, 0
41
+    for i = 0, dom.Nd/2-1 do
42
+      local dx = r[i].x - r0[i].x
43
+      local dy = r[i].y - r0[i].y
44
+      local dz = r[i].z - r0[i].z
45
+      x = x + (dx*dx - x) / (i+1)
46
+      y = y + (dy*dy - y) / (i+1)
47
+      z = z + (dz*dz - z) / (i+1)
48
+    end
49
+    result[0](x); result[1](y); result[2](z)
50
+  end
51
+
52
+  do
53
+    local space_r_cm = hdf5.create_simple_space({dom.Nd/2, 3})
54
+
55
+    function self.snapshot(group, r_cm)
56
+      local dset = group:create_dataset("r_cm", hdf5.double, space_r_cm)
57
+      dset:write(r_cm, hdf5.double, space_r_cm)
58
+    end
59
+
60
+    function self.restore(group)
61
+      local dset = group:open_dataset("r_cm")
62
+      local r_cm = vector_n(dom.Nd/2)
63
+      dset:read(r_cm, hdf5.double, space_r_cm)
64
+      return r_cm
65
+    end
66
+  end
67
+
68
+  return self
69
+end

+ 65
- 0
examples/many_dimer/production/orient.lua View File

@@ -0,0 +1,65 @@