Browse Source

Initialise nano-scale chemically powered sphere-dimer motor

Peter Colberg 3 years ago
commit
f805ce3de5
100 changed files with 8222 additions and 0 deletions
  1. 5
    0
      .gitignore
  2. 19
    0
      LICENSE
  3. 179
    0
      doc/nano-dimer.bib
  4. 1
    0
      examples/env.sh
  5. 62
    0
      examples/many_dimer/equilibration/config.lua
  6. 157
    0
      examples/many_dimer/equilibration/many_dimer.lua
  7. 178
    0
      examples/many_dimer/many_dimer.sh
  8. 72
    0
      examples/many_dimer/production/config.lua
  9. 163
    0
      examples/many_dimer/production/many_dimer.lua
  10. 67
    0
      examples/many_dimer/production/msdd.lua
  11. 66
    0
      examples/many_dimer/production/orient.lua
  12. 61
    0
      examples/many_dimer/production/vafd.lua
  13. 62
    0
      examples/single_dimer/equilibration/config.lua
  14. 138
    0
      examples/single_dimer/equilibration/single_dimer.lua
  15. 80
    0
      examples/single_dimer/production/config.lua
  16. 35
    0
      examples/single_dimer/production/cylindf.cl
  17. 134
    0
      examples/single_dimer/production/cylindf.lua
  18. 87
    0
      examples/single_dimer/production/cylinvf.cl
  19. 170
    0
      examples/single_dimer/production/cylinvf.lua
  20. 67
    0
      examples/single_dimer/production/msdd.lua
  21. 66
    0
      examples/single_dimer/production/orient.lua
  22. 40
    0
      examples/single_dimer/production/polardf.cl
  23. 136
    0
      examples/single_dimer/production/polardf.lua
  24. 45
    0
      examples/single_dimer/production/rdf.cl
  25. 127
    0
      examples/single_dimer/production/rdf.lua
  26. 171
    0
      examples/single_dimer/production/single_dimer.lua
  27. 61
    0
      examples/single_dimer/production/vafd.lua
  28. 186
    0
      examples/single_dimer/single_dimer.sh
  29. 73
    0
      examples/single_dimer/tools/cylindrical_density.py
  30. 86
    0
      examples/single_dimer/tools/cylindrical_velocity.py
  31. 49
    0
      examples/single_dimer/tools/dimer_force.py
  32. 62
    0
      examples/single_dimer/tools/dimer_force_distribution.py
  33. 66
    0
      examples/single_dimer/tools/dimer_mean_square_displacement.py
  34. 59
    0
      examples/single_dimer/tools/dimer_orientation_autocorrelation.py
  35. 171
    0
      examples/single_dimer/tools/dimer_trajectory.py
  36. 50
    0
      examples/single_dimer/tools/dimer_velocity.py
  37. 62
    0
      examples/single_dimer/tools/dimer_velocity_autocorrelation.py
  38. 63
    0
      examples/single_dimer/tools/dimer_velocity_distribution.py
  39. 58
    0
      examples/single_dimer/tools/polar_density.py
  40. 56
    0
      examples/single_dimer/tools/radial_density.py
  41. 33
    0
      nanomotor/box.cl
  42. 41
    0
      nanomotor/box.lua
  43. 59
    0
      nanomotor/compute.lua
  44. 16
    0
      nanomotor/config.lua
  45. 173
    0
      nanomotor/correlation.lua
  46. 77
    0
      nanomotor/cumsum.cl
  47. 67
    0
      nanomotor/cumsum.lua
  48. 184
    0
      nanomotor/domain.lua
  49. 259
    0
      nanomotor/h5md.lua
  50. 76
    0
      nanomotor/hilbert.cl
  51. 111
    0
      nanomotor/integrate.lua
  52. 109
    0
      nanomotor/lj.cl
  53. 70
    0
      nanomotor/lj.lua
  54. 105
    0
      nanomotor/mpcd.cl
  55. 129
    0
      nanomotor/mpcd.lua
  56. 88
    0
      nanomotor/neigh.cl
  57. 104
    0
      nanomotor/neigh.lua
  58. 36
    0
      nanomotor/observe.lua
  59. 86
    0
      nanomotor/position.cl
  60. 82
    0
      nanomotor/position.lua
  61. 128
    0
      nanomotor/random.cl
  62. 63
    0
      nanomotor/random.lua
  63. 77
    0
      nanomotor/rattle.cl
  64. 90
    0
      nanomotor/rattle.lua
  65. 79
    0
      nanomotor/snapshot.lua
  66. 60
    0
      nanomotor/sort.cl
  67. 111
    0
      nanomotor/sort.lua
  68. 40
    0
      nanomotor/species.cl
  69. 96
    0
      nanomotor/species.lua
  70. 70
    0
      nanomotor/statistics.lua
  71. 79
    0
      nanomotor/thermo.lua
  72. 94
    0
      nanomotor/trajd.lua
  73. 90
    0
      nanomotor/trajs.lua
  74. 78
    0
      nanomotor/velocity.cl
  75. 71
    0
      nanomotor/velocity.lua
  76. 130
    0
      nanomotor/verlet.cl
  77. 158
    0
      nanomotor/verlet.lua
  78. 40
    0
      strict.lua
  79. 19
    0
      test/Makefile
  80. 30
    0
      test/cumsum.lua
  81. 160
    0
      test/h5md.lua
  82. 26
    0
      test/hilbert.cl
  83. 72
    0
      test/hilbert.cxx
  84. 98
    0
      test/hilbert.lua
  85. 28
    0
      test/hilbert.py
  86. BIN
      test/hilbert_1x3x2x4.h5
  87. BIN
      test/hilbert_2x5.h5
  88. BIN
      test/hilbert_4x2x3.h5
  89. 76
    0
      test/position.lua
  90. 91
    0
      test/random.cl
  91. 49
    0
      test/random.cxx
  92. BIN
      test/random.h5
  93. 220
    0
      test/random.lua
  94. 105
    0
      test/sort.cxx
  95. BIN
      test/sort.h5
  96. 99
    0
      test/sort.lua
  97. 34
    0
      test/sort.py
  98. 90
    0
      test/statistics.lua
  99. 76
    0
      test/velocity.lua
  100. 0
    0
      tools/center_of_mass_velocity.py

+ 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.

+ 179
- 0
doc/nano-dimer.bib View File

@@ -0,0 +1,179 @@
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{Malevanets1999,
109
+  author = {Anatoly Malevanets and Raymond Kapral},
110
+  title = {Mesoscopic model for solvent dynamics},
111
+  journal = {Journal of Chemical Physics},
112
+  year = {1999},
113
+  volume = {110},
114
+  pages = {8605--8613},
115
+  number = {17},
116
+  doi = {10.1063/1.478857}
117
+}
118
+
119
+@article{Malevanets2000,
120
+  author = {Anatoly Malevanets and Raymond Kapral},
121
+  title = {Solute molecular dynamics in a mesoscale solvent},
122
+  journal = {Journal of Chemical Physics},
123
+  year = {2000},
124
+  volume = {112},
125
+  pages = {7260--7269},
126
+  number = {16},
127
+  doi = {10.1063/1.481289}
128
+}
129
+
130
+@article{Neumann1951,
131
+  author = {John von Neumann},
132
+  title = {Various techniques used in connection with random digits},
133
+  journal = {J. Res Nat. Bur. Stand., Applied Mathematics Series},
134
+  year = {1951},
135
+  volume = {12},
136
+  pages = {36--38}
137
+}
138
+
139
+@inproceedings{Salmon2011,
140
+  author = {Salmon, John K. and Moraes, Mark A. and Dror, Ron O. and
141
+    Shaw, David E.},
142
+  title = {Parallel random numbers: as easy as 1, 2, 3},
143
+  booktitle = {Proceedings of 2011 International Conference for High
144
+    Performance Computing, Networking, Storage and Analysis},
145
+  year = {2011},
146
+  series = {SC '11},
147
+  pages = {16:1--16:12},
148
+  address = {New York, NY, USA},
149
+  publisher = {ACM},
150
+  articleno = {16},
151
+  doi = {10.1145/2063384.2063405},
152
+  isbn = {978-1-4503-0771-0},
153
+  location = {Seattle, Washington}
154
+}
155
+
156
+@inproceedings{Sengupta2007,
157
+  author = {Shubhabrata Sengupta and Mark Harris and Yao Zhang and
158
+    John D. Owens},
159
+  title = {Scan Primitives for GPU Computing},
160
+  booktitle = {Graphics Hardware 2007},
161
+  year = {2007},
162
+  pages = {97--106},
163
+  organization = {ACM},
164
+  location = {San Diego, CA}
165
+}
166
+
167
+@article{Swope1982,
168
+  author = {William C. Swope and Hans C. Andersen and Peter H. Berens and
169
+    Kent R. Wilson},
170
+  title = {A computer simulation method for the calculation of equilibrium
171
+    constants for the formation of physical clusters of molecules: Application
172
+    to small water clusters},
173
+  journal = {Journal of Chemical Physics},
174
+  year = {1982},
175
+  volume = {76},
176
+  pages = {637--649},
177
+  number = {1},
178
+  doi = {10.1063/1.442716}
179
+}

+ 1
- 0
examples/env.sh View File

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

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

@@ -0,0 +1,62 @@
1
+-- Edge lengths of simulation domain.
2
+L = {80, 80, 80}
3
+-- Periodic boundary conditions.
4
+periodic = {true, true, true}
5
+-- Number of dimer spheres.
6
+Nd = 2 * 80
7
+-- Diameters of dimer spheres.
8
+diameter = {C = 4.0, N = 8.0}
9
+-- Solvent number density.
10
+rho = 9.0
11
+-- Masses of neutrally buoyant dimer spheres.
12
+mass = {C = rho * (math.pi/6) * diameter.C^3, N = rho * (math.pi/6) * diameter.N^3}
13
+-- Number of solvent particles excluding volume of dimer spheres.
14
+Ns = math.floor(rho*L[1]*L[2]*L[3] - (Nd/2)*(mass.C + mass.N))
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
+-- Solvent temperature.
20
+temp = 1/6
21
+-- Distance of dimer spheres.
22
+bond = (diameter.C + diameter.N)/2 * 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 = 10.0, CN = 10.0, NN = 10.0}
25
+-- Potential diameters.
26
+sigma = {C = diameter.C/2, N = diameter.N/2, CC = diameter.C + 1, CN = (diameter.C+diameter.N)/2 + 1, NN = diameter.N + 1}
27
+-- Forward and backward reaction of solvent particles.
28
+reaction = {rate = 0.01}
29
+-- Verlet neighbour list skin.
30
+skin = 1.0
31
+-- Number of placeholders per dimer-sphere neighbour list.
32
+nlist_size = 5000
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 = 150
37
+-- Multi-particle collision dynamics.
38
+mpcd = {interval = 50, bin_size = 50}
39
+-- Integration time-step.
40
+timestep = 0.01
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
+-- Count solvent species.
50
+species = {interval = 1000, species = {"A", "B"}}
51
+-- OpenCL platforms.
52
+platform = 1
53
+-- OpenCL device per platform.
54
+device = 1
55
+-- Random number generator seed.
56
+seed = "0xd54e68ef741800c7"
57
+-- H5MD output filename.
58
+output = "many_dimer.h5"
59
+-- H5MD file author.
60
+author = {name = "B.J. Alder", email = "bjalder@example.org"}
61
+-- H5MD file creator.
62
+creator = {name = "many_dimer.lua", version = "1.0"}

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

@@ -0,0 +1,157 @@
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
+  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
+  species = require("nanomotor.species"),
25
+}
26
+
27
+local cl = require("opencl")
28
+local hdf5 = require("hdf5")
29
+local ffi = require("ffi")
30
+local syscall = require("syscall")
31
+local json = require("cjson")
32
+
33
+-- Read configuration file.
34
+local args = nm.config(arg[1])
35
+
36
+-- Seed random number generator.
37
+math.randomseed(args.seed)
38
+
39
+-- Create OpenCL context on chosen OpenCL device.
40
+local platform = cl.get_platforms()[args.platform]
41
+print(platform:get_info("version"))
42
+print(platform:get_info("name"))
43
+print(platform:get_info("vendor"))
44
+print(platform:get_info("extensions"))
45
+local device = platform:get_devices()[args.device]
46
+print(device:get_info("name"))
47
+print(device:get_info("vendor"))
48
+print(device:get_info("extensions"))
49
+local context = cl.create_context({device})
50
+local queue = context:create_command_queue(device)
51
+
52
+-- Create simulation box with given boundary conditions.
53
+local box = nm.box(args)
54
+
55
+-- Allocate zero-filled particle buffers.
56
+local dom = nm.domain(queue, box, args)
57
+
58
+-- Place sphere-dimers at random positions and orientations.
59
+do
60
+  local rd, L, bond = dom.rd, box.L, args.bond
61
+  local dist2_CC = (args.sigma.CC * args.r_cut)^2
62
+  local dist2_CN = (args.sigma.CN * args.r_cut)^2
63
+  local dist2_NN = (args.sigma.NN * args.r_cut)^2
64
+  for i = 0, dom.Nd-1, 2 do ::redo::
65
+    rd[i].x = math.random() * L[1]
66
+    rd[i].y = math.random() * L[2]
67
+    rd[i].z = math.random() * L[3]
68
+    for j = 0, i-1 do
69
+      local dx = rd[i].x - rd[j].x
70
+      local dy = rd[i].y - rd[j].y
71
+      local dz = rd[i].z - rd[j].z
72
+      dx, dy, dz = box.mindist(dx, dy, dz)
73
+      local d2 = dx*dx + dy*dy + dz*dz
74
+      if d2 < (j%2 == 0 and dist2_CC or dist2_CN) then goto redo end
75
+    end
76
+    local x, y, z = nm.random.sphere()
77
+    x = rd[i].x + bond*x
78
+    y = rd[i].y + bond*y
79
+    z = rd[i].z + bond*z
80
+    rd[i+1].x, rd[i+1].y, rd[i+1].z = box.minimage(x, y, z)
81
+    for j = 0, i-1 do
82
+      local dx = rd[i+1].x - rd[j].x
83
+      local dy = rd[i+1].y - rd[j].y
84
+      local dz = rd[i+1].z - rd[j].z
85
+      dx, dy, dz = box.mindist(dx, dy, dz)
86
+      local d2 = dx*dx + dy*dy + dz*dz
87
+      if d2 < (j%2 == 0 and dist2_CN or dist2_NN) then goto redo end
88
+    end
89
+  end
90
+  queue:enqueue_write_buffer(dom.d_rd, true, dom.rd)
91
+end
92
+
93
+-- Randomly place solvent particles excluding dimer-sphere volumes.
94
+nm.position.random(dom, args)()
95
+-- Pick solvent velocities from Maxwell-Boltzmann distribution.
96
+nm.velocity.gaussian(dom, args)()
97
+
98
+-- Create MPCDMD 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.thermo then table.insert(observables, nm.thermo(integrate, file, args.thermo)) end
119
+if args.trajd then table.insert(observables, nm.trajd(integrate, file, args.trajd)) end
120
+if args.trajs then table.insert(observables, nm.trajs(integrate, file, args.trajs)) end
121
+if args.species then table.insert(observables, nm.species(integrate, file, args.species)) end
122
+
123
+-- Periodically log step, total energy, and remaining run-time.
124
+do
125
+  local interval = 1000
126
+  local self = {step = 0}
127
+
128
+  function self.observe(nsteps)
129
+    local t1 = syscall.clock_gettime("monotonic")
130
+    while self.step <= nsteps do
131
+      coroutine.yield()
132
+      local en = integrate.compute_en()
133
+      local t2 = syscall.clock_gettime("monotonic")
134
+      local t = (t2.time - t1.time) / interval * (nsteps - self.step)
135
+      local u = "s"
136
+      if t >= 60 then t, u = t/60, "m" end
137
+      if t >= 60 then t, u = t/60, "h" end
138
+      io.write("#", self.step, " ", en.en_tot, " ", ("%.1f"):format(t), u, "\n")
139
+      io.flush()
140
+      self.step = self.step + interval
141
+      t1 = t2
142
+    end
143
+  end
144
+
145
+  table.insert(observables, self)
146
+end
147
+
148
+-- Run simulation.
149
+nm.observe(integrate, observables, args.nsteps)
150
+
151
+-- Synchronize output files to disk.
152
+file:flush()
153
+assert((syscall.fsync(ffi.cast("int *", file:get_vfd_handle())[0])))
154
+file:close()
155
+
156
+-- Rename H5MD output file after successful simulation.
157
+assert((syscall.rename(args.output..".tmp", args.output)))

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

@@ -0,0 +1,178 @@
1
+#PBS -q gpu
2
+#PBS -l nodes=1:ppn=2:gpus=2
3
+#PBS -l walltime=12:00:00
4
+#PBS -m n
5
+#PBS -V
6
+set -e -o pipefail
7
+PROJECT="${PBS_O_WORKDIR}/nano-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}/equilibration/${OUTPUT}_${TASK}.h5" ]
15
+  then
16
+    mkdir -p "${PBS_O_WORKDIR}/equilibration"
17
+    cat >"${PBS_O_WORKDIR}/equilibration/${OUTPUT}_${TASK}.lua" <<EOF
18
+-- Edge lengths of simulation domain.
19
+L = {80, 80, 80}
20
+-- Periodic boundary conditions.
21
+periodic = {true, true, true}
22
+-- Number of dimer spheres.
23
+Nd = 2 * 80
24
+-- Diameters of dimer spheres.
25
+diameter = {C = 4.0, N = 8.0}
26
+-- Solvent number density.
27
+rho = 9.0
28
+-- Masses of neutrally buoyant dimer spheres.
29
+mass = {C = rho * (math.pi/6) * diameter.C^3, N = rho * (math.pi/6) * diameter.N^3}
30
+-- Number of solvent particles excluding volume of dimer spheres.
31
+Ns = math.floor(rho*L[1]*L[2]*L[3] - (Nd/2)*(mass.C + mass.N))
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
+-- Solvent temperature.
37
+temp = 1/6
38
+-- Distance of dimer spheres.
39
+bond = (diameter.C + diameter.N)/2 * r_cut
40
+-- Potential well depths of dimer spheres with solvent particles.
41
+epsilon = {CA = 1.0, CB = 1.0, NA = 1.0, NB = 0.1, CC = 10.0, CN = 10.0, NN = 10.0}
42
+-- Potential diameters.
43
+sigma = {C = diameter.C/2, N = diameter.N/2, CC = diameter.C + 1, CN = (diameter.C+diameter.N)/2 + 1, NN = diameter.N + 1}
44
+-- Forward and backward reaction of solvent particles.
45
+reaction = {rate = 0.01}
46
+-- Verlet neighbour list skin.
47
+skin = 1.0
48
+-- Number of placeholders per dimer-sphere neighbour list.
49
+nlist_size = 5000
50
+-- Number of bins per dimension for solvent particles.
51
+nbin = {L[1]/2, L[2]/2, L[3]/2}
52
+-- Number of placeholders per bin.
53
+bin_size = 150
54
+-- Multi-particle collision dynamics.
55
+mpcd = {interval = 50, bin_size = 50}
56
+-- Integration time-step.
57
+timestep = 0.01
58
+-- Number of steps.
59
+nsteps = 200000
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
+-- Count solvent species.
67
+species = {interval = 1000, species = {"A", "B"}}
68
+-- OpenCL platforms.
69
+platform = 1
70
+-- OpenCL device per platform.
71
+device = ${TASK}
72
+-- Random number generator seed.
73
+seed = "0x$(echo -ne "${SEED}\0${TASK}\0equilibration" | md5sum | cut -c-16)"
74
+-- H5MD output filename.
75
+output = "${PBS_O_WORKDIR}/equilibration/${OUTPUT}_${TASK}.h5"
76
+-- H5MD file author.
77
+author = {name = "$(git config user.name)", email = "$(git config user.email)"}
78
+-- H5MD file creator.
79
+creator = {name = "many_dimer.lua", version = "$(cd "${PROJECT}" && git describe --always)"}
80
+EOF
81
+    cd "${PROJECT}/examples/many_dimer/equilibration"
82
+    ./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" &
83
+    TASKS+=("${!}")
84
+  fi
85
+done
86
+for TASK in "${TASKS[@]}"
87
+do
88
+  wait "${TASK}"
89
+done
90
+TASKS=()
91
+for TASK in 1 2
92
+do
93
+  if [ ! -e "${PBS_O_WORKDIR}/production/${OUTPUT}_${TASK}.h5" ]
94
+  then
95
+    mkdir -p "${PBS_O_WORKDIR}/production"
96
+    cat >"${PBS_O_WORKDIR}/production/${OUTPUT}_${TASK}.lua" <<EOF
97
+-- Edge lengths of simulation domain.
98
+L = {80, 80, 80}
99
+-- Periodic boundary conditions.
100
+periodic = {true, true, true}
101
+-- Number of dimer spheres.
102
+Nd = 2 * 80
103
+-- Diameters of dimer spheres.
104
+diameter = {C = 4.0, N = 8.0}
105
+-- Solvent number density.
106
+rho = 9.0
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 dimer spheres.
110
+Ns = math.floor(rho*L[1]*L[2]*L[3] - (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 * 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 = 10.0, CN = 10.0, NN = 10.0}
119
+-- Potential diameters.
120
+sigma = {C = diameter.C/2, N = diameter.N/2, CC = diameter.C + 1, CN = (diameter.C+diameter.N)/2 + 1, NN = diameter.N + 1}
121
+-- Forward and backward reaction of solvent particles.
122
+reaction = {rate = 0.01}
123
+-- Verlet neighbour list skin.
124
+skin = 1.0
125
+-- Number of placeholders per dimer-sphere neighbour list.
126
+nlist_size = 5000
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 = 150
131
+-- Multi-particle collision dynamics.
132
+mpcd = {interval = 50, bin_size = 50}
133
+-- Integration time-step.
134
+timestep = 0.01
135
+-- Number of steps.
136
+nsteps = 10000000
137
+-- Snapshot program state.
138
+snapshot = {interval = 1000000, output = "${PBS_O_WORKDIR}/production/${OUTPUT}_${TASK}-snapshot.h5"}
139
+-- Thermodynamic variables.
140
+thermo = {interval = 10000, datatype = "double"}
141
+-- Dimer trajectory.
142
+trajd = {interval = 1000, datatype = "float", species = {"C", "N"}}
143
+-- Solvent trajectory.
144
+trajs = {interval = nsteps, datatype = "float", species = {"A", "B"}}
145
+-- Count solvent species.
146
+species = {interval = 10000, species = {"A", "B"}}
147
+-- Mean-square displacement of dimer.
148
+msdd = {interval = 100, blockcount = 5, blocksize = 100, blockshift = 10, separation = 10000}
149
+-- Velocity autocorrelation function of dimer.
150
+vafd = {interval = 100, blockcount = 5, blocksize = 100, blockshift = 10, separation = 10000}
151
+-- Orientation autocorrelation of dimer.
152
+orient = {interval = 100, blockcount = 5, blocksize = 100, blockshift = 10, separation = 10000}
153
+-- OpenCL platforms.
154
+platform = 1
155
+-- OpenCL device per platform.
156
+device = ${TASK}
157
+-- Random number generator seed.
158
+seed = "0x$(echo -ne "${SEED}\0${TASK}\0production" | md5sum | cut -c-16)"
159
+-- H5MD input filename.
160
+input = "${PBS_O_WORKDIR}/equilibration/${OUTPUT}_${TASK}.h5"
161
+-- Read trajectory at step.
162
+step = 200000
163
+-- H5MD output filename.
164
+output = "${PBS_O_WORKDIR}/production/${OUTPUT}_${TASK}.h5"
165
+-- H5MD file author.
166
+author = {name = "$(git config user.name)", email = "$(git config user.email)"}
167
+-- H5MD file creator.
168
+creator = {name = "many_dimer.lua", version = "$(cd "${PROJECT}" && git describe --always)"}
169
+EOF
170
+    cd "${PROJECT}/examples/many_dimer/production"
171
+    ./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" &
172
+    TASKS+=("${!}")
173
+  fi
174
+done
175
+for TASK in "${TASKS[@]}"
176
+do
177
+  wait "${TASK}"
178
+done

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

@@ -0,0 +1,72 @@
1
+-- Edge lengths of simulation domain.
2
+L = {80, 80, 80}
3
+-- Periodic boundary conditions.
4
+periodic = {true, true, true}
5
+-- Number of dimer spheres.
6
+Nd = 2 * 80
7
+-- Diameters of dimer spheres.
8
+diameter = {C = 4.0, N = 8.0}
9
+-- Solvent number density.
10
+rho = 9.0
11
+-- Masses of neutrally buoyant dimer spheres.
12
+mass = {C = rho * (math.pi/6) * diameter.C^3, N = rho * (math.pi/6) * diameter.N^3}
13
+-- Number of solvent particles excluding volume of dimer spheres.
14
+Ns = math.floor(rho*L[1]*L[2]*L[3] - (Nd/2)*(mass.C + mass.N))
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
+-- Distance of dimer spheres.
20
+bond = (diameter.C + diameter.N)/2 * r_cut
21
+-- Potential well depths of dimer spheres with solvent particles.
22
+epsilon = {CA = 1.0, CB = 1.0, NA = 1.0, NB = 0.1, CC = 10.0, CN = 10.0, NN = 10.0}
23
+-- Potential diameters.
24
+sigma = {C = diameter.C/2, N = diameter.N/2, CC = diameter.C + 1, CN = (diameter.C+diameter.N)/2 + 1, NN = diameter.N + 1}
25
+-- Forward and backward reaction of solvent particles.
26
+reaction = {rate = 0.01}
27
+-- Verlet neighbour list skin.
28
+skin = 1.0
29
+-- Number of placeholders per dimer-sphere neighbour list.
30
+nlist_size = 5000
31
+-- Number of bins per dimension for solvent particles.
32
+nbin = {L[1]/2, L[2]/2, L[3]/2}
33
+-- Number of placeholders per bin.
34
+bin_size = 150
35
+-- Multi-particle collision dynamics.
36
+mpcd = {interval = 50, bin_size = 50}
37
+-- Integration time-step.
38
+timestep = 0.01
39
+-- Number of steps.
40
+nsteps = 100000
41
+-- Snapshot program state.
42
+snapshot = {interval = 10000, output = "many_dimer-snapshot.h5"}
43
+-- Thermodynamic variables.
44
+thermo = {interval = 100, datatype = "double"}
45
+-- Dimer trajectory.
46
+trajd = {interval = 100, datatype = "float", species = {"C", "N"}}
47
+-- Solvent trajectory.
48
+trajs = {interval = nsteps, datatype = "float", species = {"A", "B"}}
49
+-- Count solvent species.
50
+species = {interval = 1000, species = {"A", "B"}}
51
+-- Mean-square displacement of dimer.
52
+msdd = {interval = 100, blockcount = 5, blocksize = 100, blockshift = 10, separation = 10000}
53
+-- Velocity autocorrelation function of dimer.
54
+vafd = {interval = 100, blockcount = 5, blocksize = 100, blockshift = 10, separation = 10000}
55
+-- Orientation autocorrelation of dimer.
56
+orient = {interval = 100, blockcount = 5, blocksize = 100, blockshift = 10, separation = 10000}
57
+-- OpenCL platforms.
58
+platform = 1
59
+-- OpenCL device per platform.
60
+device = 1
61
+-- Random number generator seed.
62
+seed = "0x1ad4fc482a4192e6"
63
+-- H5MD input filename.
64
+input = "../equilibration/many_dimer.h5"
65
+-- Read trajectory at step.
66
+step = 100000
67
+-- H5MD output filename.
68
+output = "many_dimer.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
+  snapshot = require("nanomotor.snapshot"),
16
+  integrate = require("nanomotor.integrate"),
17
+  observe = require("nanomotor.observe"),
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 MPCDMD 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

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

@@ -0,0 +1,67 @@
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 = {}, -- scalar
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 m = 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
+      m = m + ((dx*dx + dy*dy + dz*dz) - m) / (i+1)
46
+    end
47
+    result(m)
48
+  end
49
+
50
+  do
51
+    local space_r_cm = hdf5.create_simple_space({dom.Nd/2, 3})
52
+
53
+    function self.snapshot(group, r_cm)
54
+      local dset = group:create_dataset("r_cm", hdf5.double, space_r_cm)
55
+      dset:write(r_cm, hdf5.double, space_r_cm)
56
+    end
57
+
58
+    function self.restore(group)
59
+      local dset = group:open_dataset("r_cm")
60
+      local r_cm = vector_n(dom.Nd/2)
61
+      dset:read(r_cm, hdf5.double, space_r_cm)
62
+      return r_cm
63
+    end
64
+  end
65
+
66
+  return self
67
+end

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

@@ -0,0 +1,66 @@
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 = {}, -- scalar
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 m = 0
40
+    for i = 0, dom.Nd/2-1 do
41
+      local x = d[i].x*d0[i].x
42
+      local y = d[i].y*d0[i].y
43
+      local z = d[i].z*d0[i].z
44
+      m = m + ((x + y + z) - m) / (i+1)
45
+    end
46
+    result(m)
47
+  end
48
+
49
+  do
50
+    local space_d = hdf5.create_simple_space({dom.Nd/2, 3})
51
+
52
+    function self.snapshot(group, d)
53
+      local dset = group:create_dataset("d", hdf5.double, space_d)
54
+      dset:write(d, hdf5.double, space_d)
55
+    end
56
+
57
+    function self.restore(group)
58
+      local dset = group:open_dataset("d")
59
+      local d = vector_n(dom.Nd/2)
60
+      dset:read(d, hdf5.double, space_d)
61
+      return d
62
+    end
63
+  end
64
+
65
+  return self
66
+end

+ 61
- 0
examples/many_dimer/production/vafd.lua View File

@@ -0,0 +1,61 @@
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 = {}, -- scalar
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 m = 0
35
+    for i = 0, dom.Nd/2-1 do
36
+      local x = v[i].x*v0[i].x
37
+      local y = v[i].y*v0[i].y
38
+      local z = v[i].z*v0[i].z
39
+      m = m + ((x + y + z) - m) / (i+1)
40
+    end
41
+    result(m)
42
+  end
43
+
44
+  do
45
+    local space_v_cm = hdf5.create_simple_space({dom.Nd/2, 3})
46
+
47
+    function self.snapshot(group, v_cm)
48
+      local dset = group:create_dataset("v_cm", hdf5.double, space_v_cm)
49
+      dset:write(v_cm, hdf5.double, space_v_cm)
50
+    end
51
+
52
+    function self.restore(group)
53
+      local dset = group:open_dataset("v_cm")
54
+      local v_cm = vector_n(dom.Nd/2)
55
+      dset:read(v_cm, hdf5.double, space_v_cm)
56
+      return v_cm
57
+    end
58
+  end
59
+
60
+  return self
61
+end

+ 62
- 0
examples/single_dimer/equilibration/config.lua View File

@@ -0,0 +1,62 @@
1
+-- Edge lengths of simulation domain.
2
+L = {50, 50, 50}
3
+-- Periodic boundary conditions.
4
+periodic = {true, true, true}
5
+-- Number of dimer spheres.
6
+Nd = 2
7
+-- Diameters of dimer spheres.
8
+diameter = {C = 4.0, N = 8.0}
9
+-- Solvent number density.
10
+rho = 9.0
11
+-- Number of solvent particles.
12
+Ns = math.floor(rho*L[1]*L[2]*L[3])
13
+-- Cutoff radius of potential in units of σ.
14
+r_cut = 2^(1/6)
15
+-- Minimum number of steps between solvent particle sorts.
16
+sort = {interval = 100}
17
+-- Solvent temperature.
18
+temp = 1/6
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 * 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/2, N = diameter.N/2}
27
+-- Forward and backward reaction of solvent particles.
28
+reaction = {radius = math.max(L[1], L[2], L[3])/2}
29
+-- Verlet neighbour list skin.
30
+skin = 1.0
31
+-- Number of placeholders per dimer-sphere neighbour list.
32
+nlist_size = 5000
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 = 150
37
+-- Multi-particle collision dynamics.
38
+mpcd = {interval = 50, bin_size = 50}
39
+-- Integration time-step.
40
+timestep = 0.01
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
+-- Count solvent species.
50
+species = {interval = 1000, species = {"A", "B"}}
51
+-- OpenCL platforms.
52
+platform = 1
53
+-- OpenCL device per platform.
54
+device = 1
55
+-- Random number generator seed.
56
+seed = "0xb6002b4159803644"
57
+-- H5MD output filename.
58
+output = "single_dimer.h5"
59
+-- H5MD file author.
60
+author = {name = "B.J. Alder", email = "bjalder@example.org"}
61
+-- H5MD file creator.
62
+creator = {name = "single_dimer.lua", version = "1.0"}

+ 138
- 0
examples/single_dimer/equilibration/single_dimer.lua View File

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

+ 80
- 0
examples/single_dimer/production/config.lua View File

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

+ 35
- 0
examples/single_dimer/production/cylindf.cl View File

@@ -0,0 +1,35 @@
1
+/*
2
+ * Cylindrical density function.
3
+ * Copyright © 2013–2015 Peter Colberg.
4
+ * Distributed under the MIT license. (See accompanying file LICENSE.)
5
+ */
6
+
7
+#ifndef CL_VERSION_1_2
8
+#pragma OPENCL EXTENSION cl_khr_fp64 : enable
9
+#endif
10
+
11
+${include "nanomotor/box.cl"}
12
+
13
+__kernel void
14
+cylindf_bin(__global const double3 *restrict d_rd,
15
+            __global const double3 *restrict d_rs,
16
+            __global const char4 *restrict d_species,
17
+            __global uint *restrict d_bin)
18
+{
19
+  const uint gid = get_global_id(0);
20
+  if (gid < ${dom.Ns}) {
21
+    const double3 z = normalize(mindist(d_rd[0] - d_rd[1]));
22
+    const double3 d = mindist(d_rs[gid] - d_rd[1]);
23
+    const double dz = dot(d, z);
24
+    if (dz > ${cutoff.z[1]} && dz < ${cutoff.z[2]}) {
25
+      const double r2 = dot(d, d) - dz*dz;
26
+      if (r2 < ${cutoff.r^2}) {
27
+        const uint bin_z = clamp(convert_int_rtn((dz - ${cutoff.z[1]}) * ${nbin.z/(cutoff.z[2]-cutoff.z[1])}), 0, ${nbin.z-1});
28
+        const uint bin_r = min(convert_int_rtn(sqrt(r2) * ${nbin.r/cutoff.r}), ${nbin.r-1});
29
+        const char s = d_species[gid].s0;
30
+        if (s == 0) atomic_inc(&d_bin[bin_z*${nbin.r*2} + bin_r*2 + 0]);
31
+        if (s == 1) atomic_inc(&d_bin[bin_z*${nbin.r*2} + bin_r*2 + 1]);
32
+      }
33
+    }
34
+  }
35
+}

+ 134
- 0
examples/single_dimer/production/cylindf.lua View File

@@ -0,0 +1,134 @@
1
+------------------------------------------------------------------------------
2
+-- Cylindrical density function.
3
+-- Copyright © 2013–2015 Peter Colberg.
4
+-- Distributed under the MIT license. (See accompanying file LICENSE.)
5
+------------------------------------------------------------------------------
6
+
7
+local compute = require("nanomotor.compute")
8
+local hdf5 = require("hdf5")
9
+local ffi = require("ffi")
10
+
11
+-- Cache C types.
12
+local scalar = hdf5.create_space("scalar")
13
+
14
+return function(integrate, file, args)
15
+  local dom = integrate.dom
16
+  local context, device, queue = dom.context, dom.device, dom.queue
17
+  local box = dom.box
18
+  local nbin, cutoff = args.nbin, args.cutoff
19
+  local interval = args.interval
20
+
21
+  local program = compute.program(context, "cylindf.cl", {
22
+    dom = dom,
23
+    box = box,
24
+    nbin = nbin,
25
+    cutoff = cutoff,
26
+  })
27
+
28
+  local bin = ffi.typeof("struct { cl_uint A, B; }[$][$]", nbin.z, nbin.r)()
29
+  local d_bin = context:create_buffer("use_host_ptr", ffi.sizeof(bin), bin)
30
+  local count = ffi.new("uint32_t[1]")
31
+
32
+  local sample do
33
+    local kernel = program:create_kernel("cylindf_bin")
34
+    local work_size = kernel:get_work_group_info(device, "preferred_work_group_size_multiple")
35
+    local glob_size = math.ceil(dom.Ns/work_size) * work_size
36
+
37
+    function sample()
38
+      kernel:set_arg(0, dom.d_rd)
39
+      kernel:set_arg(1, dom.d_rs)
40
+      kernel:set_arg(2, dom.d_species)
41
+      kernel:set_arg(3, d_bin)
42
+      queue:enqueue_ndrange_kernel(kernel, nil, {glob_size}, {work_size})
43
+      count[0] = count[0] + 1
44
+    end
45
+  end
46
+
47
+  local group = file:create_group("structure/solvent/cylindrical_density")
48
+
49
+  local delta_r = cutoff.r / nbin.r
50
+  local r = ffi.new("float[?]", nbin.r)
51
+  for i = 0, nbin.r-1 do
52
+    r[i] = (i+0.5)*delta_r
53
+  end
54
+  do
55
+    local space_r = hdf5.create_simple_space({1, nbin.r})
56
+    local dset_r = group:create_dataset("radial", hdf5.float, space_r)
57
+    dset_r:write(r, hdf5.float)
58
+  end
59
+
60
+  local delta_z = (cutoff.z[2]-cutoff.z[1]) / nbin.z
61
+  local z = ffi.new("float[?]", nbin.z)
62
+  for i = 0, nbin.z-1 do
63
+    z[i] = (i+0.5)*delta_z + cutoff.z[1]
64
+  end
65
+  do
66
+    local space_z = hdf5.create_simple_space({nbin.z, 1})
67
+    local dset_z = group:create_dataset("axial", hdf5.float, space_z)
68
+    dset_z:write(z, hdf5.float)
69
+  end
70
+
71
+  do
72
+    local space_species = hdf5.create_simple_space({2})
73
+    local dtype_species = hdf5.c_s1:copy()
74
+    dtype_species:set_size(2)
75
+    local attr_species = group:create_attribute("species", dtype_species, space_species)
76
+    attr_species:write(ffi.new("char[2][2]", "A", "B"), dtype_species)
77
+  end
78
+
79
+  local write do
80
+    local rho = ffi.typeof("struct { float A, B; }[$][$]", nbin.z, nbin.r)()
81
+    local space_rho = hdf5.create_simple_space({nbin.z, nbin.r, 2})
82
+    local dset_rho = group:create_dataset("value", hdf5.float, space_rho)
83
+    local attr_count = group:create_attribute("count", hdf5.uint32, scalar)
84
+
85
+    function write()
86
+      local norm = 2*math.pi * count[0] * delta_z * delta_r
87
+      queue:enqueue_map_buffer(d_bin, true, "read")
88
+      for i = 0, nbin.z-1 do
89
+        for j = 0, nbin.r-1 do
90
+          rho[i][j].A = bin[i][j].A / (r[j]*norm)
91
+          rho[i][j].B = bin[i][j].B / (r[j]*norm)
92
+        end
93
+      end
94
+      queue:enqueue_unmap_mem_object(d_bin, bin)
95
+      dset_rho:write(rho, hdf5.float)
96
+      attr_count:write(count, hdf5.uint32)
97
+    end
98
+  end
99
+
100
+  local self = {step = 0}
101
+
102
+  function self.observe(nsteps)
103
+    while self.step <= nsteps do
104
+      coroutine.yield()
105
+      sample()
106
+      self.step = self.step + interval
107
+    end
108
+    write()
109
+  end
110
+
111
+  do
112
+    local space_bin = hdf5.create_simple_space({nbin.z, nbin.r, 2})
113
+
114
+    function self.snapshot(group)
115
+      local dset_bin = group:create_dataset("bin", hdf5.uint32, space_bin)
116
+      local attr_count = group:create_attribute("count", hdf5.uint32, scalar)
117
+      queue:enqueue_map_buffer(d_bin, true, "read")
118
+      dset_bin:write(bin, hdf5.uint32, space_bin)
119
+      queue:enqueue_unmap_mem_object(d_bin, bin)
120
+      attr_count:write(count, hdf5.uint32)
121
+    end
122
+
123
+    function self.restore(group)
124
+      local dset_bin = group:open_dataset("bin")
125
+      local attr_count = group:open_attribute("count")
126
+      queue:enqueue_map_buffer(d_bin, true, "write")
127
+      dset_bin:read(bin, hdf5.uint32, space_bin)
128
+      queue:enqueue_unmap_mem_object(d_bin, bin)
129
+      attr_count:read(count, hdf5.uint32)
130
+    end
131
+  end
132
+
133
+  return self
134
+end

+ 87
- 0
examples/single_dimer/production/cylinvf.cl View File

@@ -0,0 +1,87 @@
1
+/*
2
+ * Cylindrical velocity field.
3
+ * Copyright © 2013–2015 Peter Colberg.
4
+ * Distributed under the MIT license. (See accompanying file LICENSE.)
5
+ */
6
+
7
+#ifndef CL_VERSION_1_2
8
+#pragma OPENCL EXTENSION cl_khr_fp64 : enable
9
+#endif
10
+
11
+${include "nanomotor/box.cl"}
12
+
13
+|local bit = require("bit")
14
+
15
+__kernel void
16
+cylinvf_fill(__global uint *restrict d_bin_count)
17
+{
18
+  const uint gid = get_global_id(0);
19
+  if (gid < ${nbin.r*nbin.z}) d_bin_count[gid] = 0;
20
+}
21
+
22
+__kernel void
23
+cylinvf_bin(__global const double3 *restrict d_rd,
24
+            __global const double3 *restrict d_vd,
25
+            __global const double3 *restrict d_rs,
26
+            __global const double3 *restrict d_vs,
27
+            __global float4 *restrict d_bin,
28
+            __global uint *restrict d_bin_count,
29
+            __global uint *restrict d_bin_err)
30
+{
31
+  const uint gid = get_global_id(0);
32
+  if (gid < ${dom.Ns}) {
33
+    double3 z = mindist(d_rd[0] - d_rd[1]);
34
+    const double zz = dot(z, z);
35
+    const double3 r_cm = d_rd[1] + ${1/(1 + dom.mass.N/dom.mass.C)}*z;
36
+    const double3 v_cm = ${1/(1 + dom.mass.N/dom.mass.C)}*d_vd[0] + ${1/(1 + dom.mass.C/dom.mass.N)}*d_vd[1];
37
+    const double3 omega = cross(z, d_vd[0] - d_vd[1]) / zz;
38
+    z *= rsqrt(zz);
39
+    const double3 d = mindist(d_rs[gid] - d_rd[1]);
40
+    const double dz = dot(d, z);
41
+    if (dz > ${cutoff.z[1]} && dz < ${cutoff.z[2]}) {
42
+      const double r2 = dot(d, d) - dz*dz;
43
+      if (r2 < ${cutoff.r^2}) {
44
+        const double r = sqrt(r2);
45
+        const double3 n = (d - dz*z) / r;
46
+        const double3 vs = d_vs[gid];
47
+        const double3 rd = mindist(d_rs[gid] - r_cm);
48
+        const double3 vd = vs - v_cm - cross(omega, rd);
49
+        const float4 value = (float4)(dot(vs, z), dot(vs, n), dot(vd, z), dot(vd, n));
50
+        const uint bin_z = clamp(convert_int_rtn((dz - ${cutoff.z[1]})*${nbin.z / (cutoff.z[2] - cutoff.z[1])}), 0, ${nbin.z-1});
51
+        const uint bin_r = min(convert_int_rtn(r*${nbin.r/cutoff.r}), ${nbin.r-1});
52
+        const uint bin = bin_z*${nbin.r} + bin_r;
53
+        const uint count = atomic_inc(&d_bin_count[bin]);
54
+        if (count < ${bin_size}) d_bin[count + bin*${bin_size}] = value;
55
+        else atomic_max(d_bin_err, count+1);
56
+      }
57
+    }
58
+  }
59
+}
60
+
61
+|local work_size = 32
62
+
63
+__kernel void __attribute__((reqd_work_group_size(${work_size}, 1, 1)))
64
+cylinvf_sum(__global const float4 *restrict d_bin,
65
+            __global const uint *restrict d_bin_count,
66
+            __global ulong *restrict d_count,
67
+            __global float4 *restrict d_value)
68
+{
69
+  const uint lid = get_local_id(0);
70
+  const uint wid = get_group_id(0);
71
+  __local float4 l_value[${work_size}];
72
+  const uint count = min(${bin_size}u, d_bin_count[wid]);
73
+  if (count != 0) {
74
+    float4 value = 0;
75
+    for (uint i = lid; i < count; i += ${work_size}) {
76
+      value += d_bin[i + wid*${bin_size}];
77
+    }
78
+    l_value[lid] = value;
79
+    |local i = work_size
80
+    |while i > 1 do
81
+    |  i = bit.rshift(i, 1)
82
+    barrier(CLK_LOCAL_MEM_FENCE);
83
+    if (lid < ${i}) l_value[lid] += l_value[lid+${i}];
84
+    |end
85
+    if (lid == 0) d_value[wid] += (l_value[0] - d_value[wid]*count) / (d_count[wid] += count);
86
+  }
87
+}

+ 170
- 0
examples/single_dimer/production/cylinvf.lua View File

@@ -0,0 +1,170 @@
1
+------------------------------------------------------------------------------
2
+-- Cylindrical velocity field.
3
+-- Copyright © 2013–2015 Peter Colberg.
4
+-- Distributed under the MIT license. (See accompanying file LICENSE.)
5
+------------------------------------------------------------------------------
6
+
7
+local compute = require("nanomotor.compute")
8
+local hdf5 = require("hdf5")
9
+local ffi = require("ffi")
10
+
11
+-- Cache C types.
12
+local scalar = hdf5.create_space("scalar")
13
+
14
+return function(integrate, file, args)
15
+  local dom = integrate.dom
16
+  local context, device, queue = dom.context, dom.device, dom.queue
17
+  local box = dom.box
18
+  local cutoff, nbin, bin_size = args.cutoff, args.nbin, args.bin_size
19
+  local interval = args.interval
20
+
21
+  local program = compute.program(context, "cylinvf.cl", {
22
+    dom = dom,
23
+    box = box,
24
+    nbin = nbin,
25
+    bin_size = bin_size,
26
+    cutoff = cutoff,
27
+  })
28
+
29
+  local d_bin = context:create_buffer(nbin.z*nbin.r*bin_size*ffi.sizeof("cl_float4"))
30
+  local d_bin_count = context:create_buffer(nbin.z*nbin.r*ffi.sizeof("cl_uint"))
31
+  local bin_err = ffi.new("cl_uint[1]")
32
+  local d_bin_err = context:create_buffer("use_host_ptr", ffi.sizeof(bin_err), bin_err)
33
+
34
+  local fill do
35
+    local kernel = program:create_kernel("cylinvf_fill")
36
+    local work_size = kernel:get_work_group_info(device, "preferred_work_group_size_multiple")
37
+    local glob_size = math.ceil(nbin.z*nbin.r/work_size) * work_size
38
+
39
+    function fill()
40
+      kernel:set_arg(0, d_bin_count)
41
+      queue:enqueue_ndrange_kernel(kernel, nil, {glob_size}, {work_size})
42
+    end
43
+  end
44
+
45
+  local bin do
46
+    local kernel = program:create_kernel("cylinvf_bin")
47
+    local work_size = kernel:get_work_group_info(device, "preferred_work_group_size_multiple")
48
+    local glob_size = math.ceil(dom.Ns/work_size) * work_size
49
+
50
+    function bin()
51
+      kernel:set_arg(0, dom.d_rd)
52
+      kernel:set_arg(1, dom.d_vd)
53
+      kernel:set_arg(2, dom.d_rs)
54
+      kernel:set_arg(3, dom.d_vs)
55
+      kernel:set_arg(4, d_bin)
56
+      kernel:set_arg(5, d_bin_count)
57
+      kernel:set_arg(6, d_bin_err)
58
+      queue:enqueue_ndrange_kernel(kernel, nil, {glob_size}, {work_size})
59
+    end
60
+  end
61
+
62
+  local count = ffi.typeof("cl_ulong[$][$]", nbin.z, nbin.r)()
63
+  local value = ffi.typeof("cl_float4[$][$]", nbin.z, nbin.r)()
64
+  local d_count = context:create_buffer("use_host_ptr", ffi.sizeof(count), count)
65
+  local d_value = context:create_buffer("use_host_ptr", ffi.sizeof(value), value)
66
+
67
+  local sum do
68
+    local kernel = program:create_kernel("cylinvf_sum")
69
+    local work_size = kernel:get_work_group_info(device, "compile_work_group_size")[1]
70
+    local glob_size = nbin.z*nbin.r * work_size
71
+
72
+    function sum()
73
+      kernel:set_arg(0, d_bin)
74
+      kernel:set_arg(1, d_bin_count)
75
+      kernel:set_arg(2, d_count)
76
+      kernel:set_arg(3, d_value)
77
+      queue:enqueue_ndrange_kernel(kernel, nil, {glob_size}, {work_size})
78
+      queue:enqueue_map_buffer(d_bin_err, true, "read")
79
+      if bin_err[0] ~= 0 then error("bin overflow of "..bin_err[0].." solvent particles") end
80
+      queue:enqueue_unmap_mem_object(d_bin_err, bin_err)
81
+    end
82
+  end
83
+
84
+  local group = file:create_group("fields/solvent/cylindrical_velocity")
85
+
86
+  do
87
+    local z = ffi.new("float[?]", nbin.z)
88
+    local delta_z = (cutoff.z[2]-cutoff.z[1]) / nbin.z
89
+    for i = 0, nbin.z-1 do
90
+      z[i] = (i+0.5)*delta_z + cutoff.z[1]
91
+    end
92
+    local space_z = hdf5.create_simple_space({nbin.z, 1})
93
+    local space_cutoff_z = hdf5.create_simple_space({2})
94
+    local dset_z = group:create_dataset("axial", hdf5.float, space_z)
95
+    local attr_cutoff_z = dset_z:create_attribute("cutoff", hdf5.double, space_cutoff_z)
96
+    dset_z:write(z, hdf5.float)
97
+    attr_cutoff_z:write(ffi.new("double[2]", cutoff.z), hdf5.double)
98
+  end
99
+
100
+  do
101
+    local r = ffi.new("float[?]", nbin.r)
102
+    local delta_r = cutoff.r / nbin.r
103
+    for i = 0, nbin.r-1 do
104
+      r[i] = (i+0.5)*delta_r
105
+    end
106
+    local space_r = hdf5.create_simple_space({1, nbin.r})
107
+    local dset_r = group:create_dataset("radial", hdf5.float, space_r)
108
+    local attr_cutoff_r = dset_r:create_attribute("cutoff", hdf5.double, scalar)
109
+    dset_r:write(r, hdf5.float)
110
+    attr_cutoff_r:write(ffi.new("double[1]", cutoff.r), hdf5.double)
111
+  end
112
+
113
+  local write do
114
+    local space_count = hdf5.create_simple_space({nbin.z, nbin.r})
115
+    local space_value = hdf5.create_simple_space({nbin.z, nbin.r, 2, 2})
116
+    local dset_count = group:create_dataset("count", hdf5.uint64, space_count)
117
+    local dset_value = group:create_dataset("value", hdf5.float, space_value)
118
+
119
+    function write()
120
+      queue:enqueue_map_buffer(d_count, true, "read")
121
+      queue:enqueue_map_buffer(d_value, true, "read")
122
+      dset_count:write(count, hdf5.uint64)
123
+      dset_value:write(value, hdf5.float)
124
+      queue:enqueue_unmap_mem_object(d_count, count)
125
+      queue:enqueue_unmap_mem_object(d_value, value)
126
+    end
127
+  end
128
+
129
+  local self = {step = 0}
130
+
131
+  function self.observe(nsteps)
132
+    while self.step <= nsteps do
133
+      coroutine.yield()
134
+      fill()
135
+      bin()
136
+      sum()
137
+      self.step = self.step + interval
138
+    end
139
+    write()
140
+  end
141
+
142
+  do
143
+    local space_count = hdf5.create_simple_space({nbin.z, nbin.r})
144
+    local space_value = hdf5.create_simple_space({nbin.z, nbin.r, 4})
145
+
146
+    function self.snapshot(group)
147
+      local dset_count = group:create_dataset("count", hdf5.uint64, space_count)
148
+      local dset_value = group:create_dataset("value", hdf5.float, space_value)
149
+      queue:enqueue_map_buffer(d_count, true, "read")
150
+      queue:enqueue_map_buffer(d_value, true, "read")
151
+      dset_count:write(count, hdf5.uint64, space_count)
152
+      dset_value:write(value, hdf5.float, space_value)
153
+      queue:enqueue_unmap_mem_object(d_count, count)
154
+      queue:enqueue_unmap_mem_object(d_value, value)
155
+    end
156
+
157
+    function self.restore(group)
158
+      local dset_count = group:open_dataset("count")
159
+      local dset_value = group:open_dataset("value")
160
+      queue:enqueue_map_buffer(d_count, true, "write")
161
+      queue:enqueue_map_buffer(d_value, true, "write")
162
+      dset_count:read(count, hdf5.uint64, space_count)
163
+      dset_value:read(value, hdf5.float, space_value)
164
+      queue:enqueue_unmap_mem_object(d_count, count)
165
+      queue:enqueue_unmap_mem_object(d_value, value)
166
+    end
167
+  end
168
+
169
+  return self
170
+end

+ 67
- 0
examples/single_dimer/production/msdd.lua View File

@@ -0,0 +1,67 @@
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 = {}, -- scalar
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 m = 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
+      m = m + ((dx*dx + dy*dy + dz*dz) - m) / (i+1)
46
+    end
47
+    result(m)
48
+  end
49
+
50
+  do
51
+    local space_r_cm = hdf5.create_simple_space({dom.Nd/2, 3})
52
+
53
+    function self.snapshot(group, r_cm)
54
+      local dset = group:create_dataset("r_cm", hdf5.double, space_r_cm)
55
+      dset:write(r_cm, hdf5.double, space_r_cm)
56
+    end
57
+
58
+    function self.restore(group)
59
+      local dset = group:open_dataset("r_cm")
60
+      local r_cm = vector_n(dom.Nd/2)
61
+      dset:read(r_cm, hdf5.double, space_r_cm)
62
+      return r_cm
63
+    end
64
+  end
65
+
66
+  return self
67
+end

+ 66
- 0
examples/single_dimer/production/orient.lua View File

@@ -0,0 +1,66 @@
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 = {}, -- scalar
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 m = 0
40
+    for i = 0, dom.Nd/2-1 do
41
+      local x = d[i].x*d0[i].x
42
+      local y = d[i].y*d0[i].y
43
+      local z = d[i].z*d0[i].z
44
+      m = m + ((x + y + z) - m) / (i+1)
45
+    end
46
+    result(m)
47
+  end
48
+
49
+  do
50
+    local space_d = hdf5.create_simple_space({dom.Nd/2, 3})
51
+
52
+    function self.snapshot(group, d)
53
+      local dset = group:create_dataset("d", hdf5.double, space_d)
54
+      dset:write(d, hdf5.double, space_d)
55
+    end
56
+
57
+    function self.restore(group)
58
+      local dset = group:open_dataset("d")
59
+      local d = vector_n(dom.Nd/2)
60
+      dset:read(d, hdf5.double, space_d)
61
+      return d
62
+    end
63
+  end
64
+
65
+  return self
66
+end

+ 40
- 0
examples/single_dimer/production/polardf.cl View File

@@ -0,0 +1,40 @@
1
+/*
2
+ * Polar density function.
3
+ * Copyright © 2013–2015 Peter Colberg.
4
+ * Distributed under the MIT license. (See accompanying file LICENSE.)
5
+ */
6
+
7
+#ifndef CL_VERSION_1_2
8
+#pragma OPENCL EXTENSION cl_khr_fp64 : enable
9
+#endif
10
+
11
+${include "nanomotor/box.cl"}
12
+
13
+__kernel void
14
+polardf_bin(__global const double3 *restrict d_rd,
15
+            __global const double3 *restrict d_rs,
16
+            __global const char4 *restrict d_species,
17
+            __global uint *restrict d_bin)
18
+{
19
+  const uint gid = get_global_id(0);
20
+  if (gid < ${dom.Ns}) {
21
+    const double3 rs = d_rs[gid];
22
+    const char ss = d_species[gid].s0;
23
+    for (uint i = 0; i < 2; ++i) {
24
+      const double3 z = normalize(mindist(d_rd[i^1] - d_rd[i]));
25
+      const double3 d = mindist(rs - d_rd[i]);
26
+      const double r2 = dot(d, d);
27
+      if (r2 < ${cutoff.r^2}) {
28
+        const double r = sqrt(r2);
29
+        const double theta = acos(dot(d, z) / r);
30
+        const uint bin_r = min(convert_int_rtn(r * ${nbin.r/cutoff.r}), ${nbin.r-1});
31
+        const uint bin_theta = clamp(convert_int_rtn(theta * ${nbin.theta/math.pi}), 0, ${nbin.theta-1});
32
+        const char sd = i&1;
33
+        if (sd == 0 && ss == 0) atomic_inc(&d_bin[bin_r*${nbin.theta*4} + bin_theta*4 + 0]);
34
+        if (sd == 0 && ss == 1) atomic_inc(&d_bin[bin_r*${nbin.theta*4} + bin_theta*4 + 1]);
35
+        if (sd == 1 && ss == 0) atomic_inc(&d_bin[bin_r*${nbin.theta*4} + bin_theta*4 + 2]);
36
+        if (sd == 1 && ss == 1) atomic_inc(&d_bin[bin_r*${nbin.theta*4} + bin_theta*4 + 3]);
37
+      }
38
+    }
39
+  }
40
+}

+ 136
- 0
examples/single_dimer/production/polardf.lua View File

@@ -0,0 +1,136 @@
1
+------------------------------------------------------------------------------
2
+-- Polar density function.
3
+-- Copyright © 2013–2015 Peter Colberg.
4
+-- Distributed under the MIT license. (See accompanying file LICENSE.)
5
+------------------------------------------------------------------------------
6
+
7
+local compute = require("nanomotor.compute")
8
+local hdf5 = require("hdf5")
9
+local ffi = require("ffi")
10
+
11
+-- Cache C types.
12
+local scalar = hdf5.create_space("scalar")
13
+
14
+return function(integrate, file, args)
15
+  local dom = integrate.dom
16
+  local context, device, queue = dom.context, dom.device, dom.queue
17
+  local box = dom.box
18
+  local nbin, cutoff = args.nbin, args.cutoff
19
+  local interval = args.interval
20
+
21
+  local program = compute.program(context, "polardf.cl", {
22
+    dom = dom,
23
+    box = box,
24
+    nbin = nbin,
25
+    cutoff = cutoff,
26
+  })
27
+
28
+  local bin = ffi.typeof("struct { cl_uint CA, CB, NA, NB; }[$][$]", nbin.r, nbin.theta)()
29
+  local d_bin = context:create_buffer("use_host_ptr", ffi.sizeof(bin), bin)
30
+  local count = ffi.new("uint32_t[1]")
31
+
32
+  local sample do
33
+    local kernel = program:create_kernel("polardf_bin")
34
+    local work_size = kernel:get_work_group_info(device, "preferred_work_group_size_multiple")
35
+    local glob_size = math.ceil(dom.Ns/work_size) * work_size
36
+
37
+    function sample()
38
+      kernel:set_arg(0, dom.d_rd)
39
+      kernel:set_arg(1, dom.d_rs)
40
+      kernel:set_arg(2, dom.d_species)
41
+      kernel:set_arg(3, d_bin)
42
+      queue:enqueue_ndrange_kernel(kernel, nil, {glob_size}, {work_size})
43
+      count[0] = count[0] + 1
44
+    end
45
+  end
46
+
47
+  local group = file:create_group("structure/solvent/polar_density")
48
+
49
+  local delta_r = cutoff.r / nbin.r
50
+  local r = ffi.new("float[?]", nbin.r)
51
+  for i = 0, nbin.r-1 do
52
+    r[i] = (i+0.5)*delta_r
53
+  end
54
+  do
55
+    local space_r = hdf5.create_simple_space({nbin.r, 1})
56
+    local dset_r = group:create_dataset("radial", hdf5.float, space_r)
57
+    dset_r:write(r, hdf5.float)
58
+  end
59
+
60
+  local delta_theta = math.pi / nbin.theta
61
+  local theta = ffi.new("float[?]", nbin.theta)
62
+  for i = 0, nbin.theta-1 do
63
+    theta[i] = (i+0.5)*delta_theta
64
+  end
65
+  do
66
+    local space_theta = hdf5.create_simple_space({1, nbin.theta})
67
+    local dset_theta = group:create_dataset("polar", hdf5.float, space_theta)
68
+    dset_theta:write(theta, hdf5.float)
69
+  end
70
+
71
+  do
72
+    local space_species = hdf5.create_simple_space({2, 2})
73
+    local dtype_species = hdf5.c_s1:copy()
74
+    dtype_species:set_size(3)
75
+    local attr_species = group:create_attribute("species", dtype_species, space_species)
76
+    attr_species:write(ffi.new("char[4][3]", "CA", "CB", "NA", "NB"), dtype_species)
77
+  end
78
+
79
+  local write do
80
+    local rho = ffi.typeof("struct { float CA, CB, NA, NB; }[$][$]", nbin.r, nbin.theta)()
81
+    local space_rho = hdf5.create_simple_space({nbin.r, nbin.theta, 2, 2})
82
+    local dset_rho = group:create_dataset("value", hdf5.float, space_rho)
83
+    local attr_count = group:create_attribute("count", hdf5.uint32, scalar)
84
+
85
+    function write()
86
+      local norm = 2*math.pi * count[0] * delta_r * delta_theta
87
+      queue:enqueue_map_buffer(d_bin, true, "read")
88
+      for i = 0, nbin.r - 1 do
89
+        for j = 0, nbin.theta - 1 do
90
+          rho[i][j].CA = bin[i][j].CA / (r[i]*r[i]*math.sin(theta[j])*norm)
91
+          rho[i][j].CB = bin[i][j].CB / (r[i]*r[i]*math.sin(theta[j])*norm)
92
+          rho[i][j].NA = bin[i][j].NA / (r[i]*r[i]*math.sin(theta[j])*norm)
93
+          rho[i][j].NB = bin[i][j].NB / (r[i]*r[i]*math.sin(theta[j])*norm)
94
+        end
95
+      end
96
+      queue:enqueue_unmap_mem_object(d_bin, bin)
97
+      dset_rho:write(rho, hdf5.float)
98
+      attr_count:write(count, hdf5.uint32)
99
+    end
100
+  end
101
+
102
+  local self = {step = 0}
103
+
104
+  function self.observe(nsteps)
105
+    while self.step <= nsteps do
106
+      coroutine.yield()
107
+      sample()
108
+      self.step = self.step + interval
109
+    end
110
+    write()
111
+  end
112
+
113
+  do
114
+    local space_bin = hdf5.create_simple_space({nbin.r, nbin.theta, 4})
115
+
116
+    function self.snapshot(group)
117
+      local dset_bin = group:create_dataset("bin", hdf5.uint32, space_bin)
118
+      local attr_count = group:create_attribute("count", hdf5.uint32, scalar)
119
+      queue:enqueue_map_buffer(d_bin, true, "read")
120
+      dset_bin:write(bin, hdf5.uint32, space_bin)
121
+      queue:enqueue_unmap_mem_object(d_bin, bin)
122
+      attr_count:write(count, hdf5.uint32)
123
+    end
124
+
125
+    function self.restore(group)
126
+      local dset_bin = group:open_dataset("bin")
127
+      local attr_count = group:open_attribute("count")
128
+      queue:enqueue_map_buffer(d_bin, true, "write")
129
+      dset_bin:read(bin, hdf5.uint32, space_bin)
130
+      queue:enqueue_unmap_mem_object(d_bin, bin)
131
+      attr_count:read(count, hdf5.uint32)
132
+    end
133
+  end
134
+
135
+  return self
136
+end

+ 45
- 0
examples/single_dimer/production/rdf.cl View File

@@ -0,0 +1,45 @@
1
+/*
2
+ * Radial density function.
3
+ * Copyright © 2013–2015 Peter Colberg.
4
+ * Distributed under the MIT license. (See accompanying file LICENSE.)
5
+ */
6
+
7
+#ifndef CL_VERSION_1_2
8
+#pragma OPENCL EXTENSION cl_khr_fp64 : enable
9
+#endif
10
+
11
+${include "nanomotor/box.cl"}
12
+
13
+__kernel void
14
+rdf_bin(__global const double3 *restrict d_rd,
15
+        __global const double3 *restrict d_rs,
16
+        __global const char4 *restrict d_species,
17
+        __global uint *restrict d_bin)
18
+{
19
+  const uint gid = get_global_id(0);
20
+  if (gid < ${dom.Ns}) {
21
+    const double3 rs = d_rs[gid];
22
+    const char ss = d_species[gid].s0;
23
+    for (uint i = 0; i < 2; ++i) {
24
+      const double3 d = mindist(rs - d_rd[i]);
25
+      const double d2 = dot(d, d);
26
+      if (d2 < ${cutoff^2}) {
27
+        const uint bin = min(convert_int_rtn(sqrt(d2) * ${nbin/cutoff}), ${nbin-1});
28
+        const char sd = i&1;
29
+        const double3 z = mindist(d_rd[i] - d_rd[i^1]);
30
+        if (dot(d, z) > 0) {
31
+          if (sd == 0 && ss == 0) atomic_inc(&d_bin[8*bin + 0]);
32
+          if (sd == 0 && ss == 1) atomic_inc(&d_bin[8*bin + 2]);
33
+          if (sd == 1 && ss == 0) atomic_inc(&d_bin[8*bin + 4]);
34
+          if (sd == 1 && ss == 1) atomic_inc(&d_bin[8*bin + 6]);
35
+        }
36
+        else {
37
+          if (sd == 0 && ss == 0) atomic_inc(&d_bin[8*bin + 1]);
38
+          if (sd == 0 && ss == 1) atomic_inc(&d_bin[8*bin + 3]);
39
+          if (sd == 1 && ss == 0) atomic_inc(&d_bin[8*bin + 5]);
40
+          if (sd == 1 && ss == 1) atomic_inc(&d_bin[8*bin + 7]);
41
+        }
42
+      }
43
+    }
44
+  }
45
+}

+ 127
- 0
examples/single_dimer/production/rdf.lua View File

@@ -0,0 +1,127 @@
1
+------------------------------------------------------------------------------
2
+-- Radial density function.
3
+-- Copyright © 2013–2015 Peter Colberg.
4
+-- Distributed under the MIT license. (See accompanying file LICENSE.)
5
+------------------------------------------------------------------------------
6
+
7
+local compute = require("nanomotor.compute")
8
+local hdf5 = require("hdf5")
9
+local ffi = require("ffi")
10
+
11
+-- Cache C types.
12
+local scalar = hdf5.create_space("scalar")
13
+
14
+return function(integrate, file, args)
15
+  local dom = integrate.dom
16
+  local context, device, queue = dom.context, dom.device, dom.queue
17
+  local box = dom.box
18
+  local nbin, cutoff = args.nbin, args.cutoff
19
+  local interval = args.interval
20
+
21
+  local program = compute.program(context, "rdf.cl", {
22
+    dom = dom,
23
+    box = box,
24
+    nbin = nbin,
25
+    cutoff = cutoff,
26
+  })
27
+
28
+  local bin = ffi.new("struct { cl_uint CA1, CA2, CB1, CB2, NA1, NA2, NB1, NB2; }[?]", nbin)
29
+  local d_bin = context:create_buffer("use_host_ptr", ffi.sizeof(bin), bin)
30
+  local count = ffi.new("uint32_t[1]")
31
+
32
+  local sample do
33
+    local kernel = program:create_kernel("rdf_bin")
34
+    local work_size = kernel:get_work_group_info(device, "preferred_work_group_size_multiple")
35
+    local glob_size = math.ceil(dom.Ns/work_size) * work_size
36
+
37
+    function sample()
38
+      kernel:set_arg(0, dom.d_rd)
39
+      kernel:set_arg(1, dom.d_rs)
40
+      kernel:set_arg(2, dom.d_species)
41
+      kernel:set_arg(3, d_bin)
42
+      queue:enqueue_ndrange_kernel(kernel, nil, {glob_size}, {work_size})
43
+      count[0] = count[0] + 1
44
+    end
45
+  end
46
+
47
+  local group = file:create_group("structure/solvent/radial_density")
48
+
49
+  local r = ffi.new("float[?]", nbin)
50
+  local delta_r = cutoff / nbin
51
+  for i = 0, nbin-1 do
52
+    r[i] = (i+0.5)*delta_r
53
+  end
54
+  do
55
+    local space_r = hdf5.create_simple_space({nbin})
56
+    local dset_r = group:create_dataset("radial", hdf5.float, space_r)
57
+    dset_r:write(r, hdf5.float)
58
+  end
59
+
60
+  do
61
+    local space_species = hdf5.create_simple_space({2, 2})
62
+    local dtype_species = hdf5.c_s1:copy()
63
+    dtype_species:set_size(3)
64
+    local attr_species = group:create_attribute("species", dtype_species, space_species)
65
+    attr_species:write(ffi.new("char[4][3]", "CA", "CB", "NA", "NB"), dtype_species)
66
+  end
67
+
68
+  local write do
69
+    local rho = ffi.new("struct { float CA1, CA2, CB1, CB2, NA1, NA2, NB1, NB2; }[?]", nbin)
70
+    local space_rho = hdf5.create_simple_space({nbin, 2, 2, 2})
71
+    local dset_rho = group:create_dataset("value", hdf5.float, space_rho)
72
+    local attr_count = group:create_attribute("count", hdf5.uint32, scalar)
73
+
74
+    function write()
75
+      local norm = 2*math.pi * count[0] * delta_r
76
+      queue:enqueue_map_buffer(d_bin, true, "read")
77
+      for i = 0, nbin-1 do
78
+        rho[i].CA1 = bin[i].CA1 / (r[i]*r[i]*norm)
79
+        rho[i].CA2 = bin[i].CA2 / (r[i]*r[i]*norm)
80
+        rho[i].CB1 = bin[i].CB1 / (r[i]*r[i]*norm)
81
+        rho[i].CB2 = bin[i].CB2 / (r[i]*r[i]*norm)
82
+        rho[i].NA1 = bin[i].NA1 / (r[i]*r[i]*norm)
83
+        rho[i].NA2 = bin[i].NA2 / (r[i]*r[i]*norm)
84
+        rho[i].NB1 = bin[i].NB1 / (r[i]*r[i]*norm)
85
+        rho[i].NB2 = bin[i].NB2 / (r[i]*r[i]*norm)
86
+      end
87
+      queue:enqueue_unmap_mem_object(d_bin, bin)
88
+      dset_rho:write(rho, hdf5.float)
89
+      attr_count:write(count, hdf5.uint32)
90
+    end
91
+  end
92
+
93
+  local self = {step = 0}
94
+
95
+  function self.observe(nsteps)
96
+    while self.step <= nsteps do
97
+      coroutine.yield()
98
+      sample()
99
+      self.step = self.step + interval
100
+    end
101
+    write()
102
+  end
103
+
104
+  do
105
+    local space_bin = hdf5.create_simple_space({nbin, 8})
106
+
107
+    function self.snapshot(group)
108
+      local dset_bin = group:create_dataset("bin", hdf5.uint32, space_bin)
109
+      local attr_count = group:create_attribute("count", hdf5.uint32, scalar)
110
+      queue:enqueue_map_buffer(d_bin, true, "read")
111
+      dset_bin:write(bin, hdf5.uint32, space_bin)
112
+      queue:enqueue_unmap_mem_object(d_bin, bin)
113
+      attr_count:write(count, hdf5.uint32)
114
+    end
115
+
116
+    function self.restore(group)
117
+      local dset_bin = group:open_dataset("bin")
118
+      local attr_count = group:open_attribute("count")
119
+      queue:enqueue_map_buffer(d_bin, true, "write")
120
+      dset_bin:read(bin, hdf5.uint32, space_bin)
121
+      queue:enqueue_unmap_mem_object(d_bin, bin)
122
+      attr_count:read(count, hdf5.uint32)
123
+    end
124
+  end
125
+
126
+  return self
127
+end

+ 171
- 0
examples/single_dimer/production/single_dimer.lua View File

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

+ 61
- 0
examples/single_dimer/production/vafd.lua View File

@@ -0,0 +1,61 @@
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 = {}, -- scalar
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 m = 0
35
+    for i = 0, dom.Nd/2-1 do
36
+      local x = v[i].x*v0[i].x
37
+      local y = v[i].y*v0[i].y
38
+      local z = v[i].z*v0[i].z
39
+      m = m + ((x + y + z) - m) / (i+1)
40
+    end
41
+    result(m)
42
+  end
43
+
44
+  do
45
+    local space_v_cm = hdf5.create_simple_space({dom.Nd/2, 3})
46
+
47
+    function self.snapshot(group, v_cm)
48
+      local dset = group:create_dataset("v_cm", hdf5.double, space_v_cm)
49
+      dset:write(v_cm, hdf5.double, space_v_cm)
50
+    end
51
+
52
+    function self.restore(group)
53
+      local dset = group:open_dataset("v_cm")
54
+      local v_cm = vector_n(dom.Nd/2)
55
+      dset:read(v_cm, hdf5.double, space_v_cm)
56
+      return v_cm
57
+    end
58
+  end
59
+
60
+  return self
61
+end

+ 186
- 0
examples/single_dimer/single_dimer.sh View File

@@ -0,0 +1,186 @@
1
+#PBS -q gpu
2
+#PBS -l nodes=1:ppn=2:gpus=2