ESPResSo
Extensible Simulation Package for Research on Soft Matter Systems
Loading...
Searching...
No Matches
h5md_core.cpp
Go to the documentation of this file.
1/*
2 * Copyright (C) 2010-2022 The ESPResSo project
3 * Copyright (C) 2002,2003,2004,2005,2006,2007,2008,2009,2010
4 * Max-Planck-Institute for Polymer Research, Theory Group
5 *
6 * This file is part of ESPResSo.
7 *
8 * ESPResSo is free software: you can redistribute it and/or modify
9 * it under the terms of the GNU General Public License as published by
10 * the Free Software Foundation, either version 3 of the License, or
11 * (at your option) any later version.
12 *
13 * ESPResSo is distributed in the hope that it will be useful,
14 * but WITHOUT ANY WARRANTY; without even the implied warranty of
15 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
16 * GNU General Public License for more details.
17 *
18 * You should have received a copy of the GNU General Public License
19 * along with this program. If not, see <http://www.gnu.org/licenses/>.
20 */
21
22#include "h5md_core.hpp"
23
24#include "BoxGeometry.hpp"
25#include "Particle.hpp"
28
29#include "config/version.hpp"
30
31#include <utils/Vector.hpp>
32
33#include <boost/mpi/collectives.hpp>
34
35#include <mpi.h>
36
37#include <algorithm>
38#include <cstddef>
39#include <fstream>
40#include <functional>
41#include <iterator>
42#include <stdexcept>
43#include <string>
44#include <vector>
45
46namespace Writer {
47namespace H5md {
48
49using MultiArray3i = boost::multi_array<int, 3>;
53
54static void backup_file(const std::string &from, const std::string &to) {
55 /*
56 * If the file itself *and* a backup file exists, something must
57 * have gone wrong.
58 */
59 boost::filesystem::path pfrom(from), pto(to);
60 auto constexpr option_fail_if_exists = boost::filesystem::copy_options::none;
61 try {
62 boost::filesystem::copy_file(pfrom, pto, option_fail_if_exists);
63 } catch (const boost::filesystem::filesystem_error &) {
64 throw left_backupfile();
65 }
66}
67
68template <typename extent_type>
69static void extend_dataset(h5xx::dataset &dataset,
70 extent_type const &change_extent) {
71 auto const rank = static_cast<h5xx::dataspace>(dataset).rank();
72 auto extents = static_cast<h5xx::dataspace>(dataset).extents();
73 /* Extend the dataset for another timestep */
74 for (int i = 0; i < rank; i++) {
75 extents[i] += change_extent[i];
76 }
77 H5Dset_extent(dataset.hid(), extents.data()); // extend all dims is collective
78}
79
80template <typename value_type, typename extent_type>
81static void write_dataset(value_type const &data, h5xx::dataset &dataset,
82 extent_type const &change_extent,
83 extent_type const &offset, extent_type const &count) {
84 extend_dataset(dataset, change_extent);
85 /* write the data to the dataset. */
86 h5xx::write_dataset(dataset, data, h5xx::slice(offset, count));
87}
88
89static void write_script(std::string const &target,
90 boost::filesystem::path const &script_path) {
91 if (!script_path.empty()) {
92 std::ifstream scriptfile(script_path.string());
93 std::string buffer((std::istreambuf_iterator<char>(scriptfile)),
94 std::istreambuf_iterator<char>());
95 auto file = h5xx::file(target, h5xx::file::out);
96 auto const group = h5xx::group(file, "parameters/files");
97 h5xx::write_attribute(group, "script", buffer);
98 file.close();
99 }
100}
101
102/* Initialize the file-related variables after parameters have been set. */
103void File::init_file(std::string const &file_path) {
104 m_backup_filename = file_path + ".bak";
105 if (m_script_path.empty()) {
106 m_absolute_script_path = boost::filesystem::path();
107 } else {
108 boost::filesystem::path script_path(m_script_path);
109 m_absolute_script_path = boost::filesystem::canonical(script_path);
110 }
111 auto const file_exists = boost::filesystem::exists(file_path);
112 auto const backup_file_exists = boost::filesystem::exists(m_backup_filename);
113 /* Perform a barrier synchronization. Otherwise one process might already
114 * create the file while another still checks for its existence. */
115 m_comm.barrier();
116 if (file_exists) {
117 if (m_h5md_specification.is_compliant(file_path)) {
118 /*
119 * If the file exists and has a valid H5MD structure, let's create a
120 * backup of it. This has the advantage, that the new file can
121 * just be deleted if the simulation crashes at some point and we
122 * still have a valid trajectory backed up, from which we can restart.
123 */
124 if (m_comm.rank() == 0)
125 backup_file(file_path, m_backup_filename);
126 load_file(file_path);
127 } else {
128 throw incompatible_h5mdfile();
129 }
130 } else {
131 if (backup_file_exists)
132 throw left_backupfile();
133 create_file(file_path);
134 }
135}
136
137void File::load_datasets() {
138 for (auto const &d : m_h5md_specification.get_datasets()) {
139 if (d.is_link)
140 continue;
141 datasets[d.path()] = h5xx::dataset(m_h5md_file, d.path());
142 }
143}
144
145void File::create_groups() {
146 h5xx::group group(m_h5md_file);
147 for (auto const &d : m_h5md_specification.get_datasets()) {
148 h5xx::group new_group(group, d.group);
149 }
150}
151
152static std::vector<hsize_t> create_dims(hsize_t rank, hsize_t data_dim) {
153 switch (rank) {
154 case 3:
155 return std::vector<hsize_t>{0, 0, data_dim};
156 case 2:
157 return std::vector<hsize_t>{0, data_dim};
158 case 1:
159 return std::vector<hsize_t>{data_dim};
160 default:
161 throw std::runtime_error(
162 "H5MD Error: datasets with this dimension are not implemented\n");
163 }
164}
165
166static std::vector<hsize_t> create_chunk_dims(hsize_t rank, hsize_t data_dim) {
167 hsize_t chunk_size = (rank > 1) ? 1000 : 1;
168 switch (rank) {
169 case 3:
170 return {1, chunk_size, data_dim};
171 case 2:
172 return {1, chunk_size};
173 case 1:
174 return {chunk_size};
175 default:
176 throw std::runtime_error(
177 "H5MD Error: datasets with this dimension are not implemented\n");
178 }
179}
180
181void File::create_datasets() {
182 namespace hps = h5xx::policy::storage;
183 for (const auto &d : m_h5md_specification.get_datasets()) {
184 if (d.is_link)
185 continue;
186 auto maxdims = std::vector<hsize_t>(d.rank, H5S_UNLIMITED);
187 auto dataspace = h5xx::dataspace(create_dims(d.rank, d.data_dim), maxdims);
188 auto storage = hps::chunked(create_chunk_dims(d.rank, d.data_dim))
189 .set(hps::fill_value(-10));
190 datasets[d.path()] = h5xx::dataset(m_h5md_file, d.path(), d.type, dataspace,
191 storage, H5P_DEFAULT, H5P_DEFAULT);
192 }
193}
194
195void File::load_file(const std::string &file_path) {
196 m_h5md_file = h5xx::file(file_path, m_comm, MPI_INFO_NULL, h5xx::file::out);
197 load_datasets();
198}
199
200static void write_attributes(h5xx::file &h5md_file) {
201 auto h5md_group = h5xx::group(h5md_file, "h5md");
202 h5xx::write_attribute(h5md_group, "version",
203 boost::array<hsize_t, 2>{{1, 1}});
204 auto h5md_creator_group = h5xx::group(h5md_group, "creator");
205 h5xx::write_attribute(h5md_creator_group, "name", "ESPResSo");
206 h5xx::write_attribute(h5md_creator_group, "version", ESPRESSO_VERSION);
207 auto h5md_author_group = h5xx::group(h5md_group, "author");
208 h5xx::write_attribute(h5md_author_group, "name", "N/A");
209 auto group = h5xx::group(h5md_file, "particles/atoms/box");
210 h5xx::write_attribute(group, "dimension", 3);
211 h5xx::write_attribute(group, "boundary", "periodic");
212}
213
214void File::write_units() {
215 if (!mass_unit().empty() and (m_fields & H5MD_OUT_MASS)) {
216 h5xx::write_attribute(datasets["particles/atoms/mass/value"], "unit",
217 mass_unit());
218 }
219 if (!charge_unit().empty() and (m_fields & H5MD_OUT_CHARGE)) {
220 h5xx::write_attribute(datasets["particles/atoms/charge/value"], "unit",
221 charge_unit());
222 }
223 if (!length_unit().empty() and (m_fields & H5MD_OUT_BOX_L)) {
224 h5xx::write_attribute(datasets["particles/atoms/position/value"], "unit",
225 length_unit());
226 h5xx::write_attribute(datasets["particles/atoms/box/edges/value"], "unit",
227 length_unit());
228 }
229 if (!length_unit().empty() and (m_fields & H5MD_OUT_LE_OFF)) {
230 h5xx::write_attribute(datasets["particles/atoms/lees_edwards/offset/value"],
231 "unit", length_unit());
232 }
233 if (!velocity_unit().empty() and (m_fields & H5MD_OUT_VEL)) {
234 h5xx::write_attribute(datasets["particles/atoms/velocity/value"], "unit",
235 velocity_unit());
236 }
237 if (!force_unit().empty() and (m_fields & H5MD_OUT_FORCE)) {
238 h5xx::write_attribute(datasets["particles/atoms/force/value"], "unit",
239 force_unit());
240 }
241 if (!time_unit().empty()) {
242 h5xx::write_attribute(datasets["particles/atoms/id/time"], "unit",
243 time_unit());
244 }
245}
246
247void File::create_hard_links() {
248 std::string path_step = "particles/atoms/id/step";
249 std::string path_time = "particles/atoms/id/time";
250 for (auto &ds : m_h5md_specification.get_datasets()) {
251 if (ds.is_link) {
252 char const *from = nullptr;
253 if (ds.name == "step") {
254 from = path_step.c_str();
255 } else if (ds.name == "time") {
256 from = path_time.c_str();
257 }
258 assert(from != nullptr);
259 if (H5Lcreate_hard(m_h5md_file.hid(), from, m_h5md_file.hid(),
260 ds.path().c_str(), H5P_DEFAULT, H5P_DEFAULT) < 0) {
261 throw std::runtime_error("Error creating hard link for " + ds.path());
262 }
263 }
264 }
265}
266
267void File::create_file(const std::string &file_path) {
268 if (m_comm.rank() == 0)
269 write_script(file_path, m_absolute_script_path);
270 m_comm.barrier();
271 m_h5md_file = h5xx::file(file_path, m_comm, MPI_INFO_NULL, h5xx::file::out);
272 create_groups();
273 create_datasets();
274 write_attributes(m_h5md_file);
275 write_units();
276 create_hard_links();
277}
278
280 if (m_comm.rank() == 0)
281 boost::filesystem::remove(m_backup_filename);
282}
283
284namespace detail {
285
286template <std::size_t rank> struct slice_info {};
287
288template <> struct slice_info<3> {
289 static auto extent(hsize_t n_part_diff) {
290 return Vector3hs{1, n_part_diff, 0};
291 }
292 static constexpr auto count() { return Vector3hs{1, 1, 3}; }
293 static auto offset(hsize_t n_time_steps, hsize_t prefix) {
294 return Vector3hs{n_time_steps, prefix, 0};
295 }
296};
297
298template <> struct slice_info<2> {
299 static auto extent(hsize_t n_part_diff) { return Vector2hs{1, n_part_diff}; }
300 static constexpr auto count() { return Vector2hs{1, 1}; }
301 static auto offset(hsize_t n_time_steps, hsize_t prefix) {
302 return Vector2hs{n_time_steps, prefix};
303 }
304};
305
306} // namespace detail
307
308template <std::size_t dim, typename Op>
309void write_td_particle_property(hsize_t prefix, hsize_t n_part_global,
310 ParticleRange const &particles,
311 h5xx::dataset &dataset, Op op) {
312 auto const old_extents = static_cast<h5xx::dataspace>(dataset).extents();
313 auto const extent_particle_number =
314 std::max(n_part_global, old_extents[1]) - old_extents[1];
315 extend_dataset(dataset,
316 detail::slice_info<dim>::extent(extent_particle_number));
317 auto const count = detail::slice_info<dim>::count();
318 auto offset = detail::slice_info<dim>::offset(old_extents[0], prefix);
319 for (auto const &p : particles) {
320 h5xx::write_dataset(dataset, op(p), h5xx::slice(offset, count));
321 // advance in the particle dimension
322 offset[1] += 1;
323 }
324}
325
326static void write_box(BoxGeometry const &box_geo, h5xx::dataset &dataset) {
327 auto const extents = static_cast<h5xx::dataspace>(dataset).extents();
328 extend_dataset(dataset, Vector2hs{1, 0});
329 h5xx::write_dataset(dataset, box_geo.length(),
330 h5xx::slice(Vector2hs{extents[0], 0}, Vector2hs{1, 3}));
331}
332
333static void write_le_off(LeesEdwardsBC const &lebc, h5xx::dataset &dataset) {
334 auto const extents = static_cast<h5xx::dataspace>(dataset).extents();
335 extend_dataset(dataset, Vector2hs{1, 0});
336 h5xx::write_dataset(dataset, Utils::Vector<double, 1>{lebc.pos_offset},
337 h5xx::slice(Vector2hs{extents[0], 0}, Vector2hs{1, 1}));
338}
339
340static void write_le_dir(LeesEdwardsBC const &lebc, h5xx::dataset &dataset) {
341 auto const shear_direction = static_cast<int>(lebc.shear_direction);
342 auto const extents = static_cast<h5xx::dataspace>(dataset).extents();
343 extend_dataset(dataset, Vector2hs{1, 0});
344 h5xx::write_dataset(dataset, Utils::Vector<int, 1>{shear_direction},
345 h5xx::slice(Vector2hs{extents[0], 0}, Vector2hs{1, 1}));
346}
347
348static void write_le_normal(LeesEdwardsBC const &lebc, h5xx::dataset &dataset) {
349 auto const shear_plane_normal = static_cast<int>(lebc.shear_plane_normal);
350 auto const extents = static_cast<h5xx::dataspace>(dataset).extents();
351 extend_dataset(dataset, Vector2hs{1, 0});
352 h5xx::write_dataset(dataset, Utils::Vector<int, 1>{shear_plane_normal},
353 h5xx::slice(Vector2hs{extents[0], 0}, Vector2hs{1, 1}));
354}
355
356void File::write(const ParticleRange &particles, double time, int step,
357 BoxGeometry const &box_geo) {
358 if (m_fields & H5MD_OUT_BOX_L) {
359 write_box(box_geo, datasets["particles/atoms/box/edges/value"]);
360 }
361 auto const &lebc = box_geo.lees_edwards_bc();
362 if (m_fields & H5MD_OUT_LE_OFF) {
363 write_le_off(lebc, datasets["particles/atoms/lees_edwards/offset/value"]);
364 }
365 if (m_fields & H5MD_OUT_LE_DIR) {
366 write_le_dir(lebc,
367 datasets["particles/atoms/lees_edwards/direction/value"]);
368 }
369 if (m_fields & H5MD_OUT_LE_NORMAL) {
370 write_le_normal(lebc,
371 datasets["particles/atoms/lees_edwards/normal/value"]);
372 }
373
374 auto const n_part_local = static_cast<int>(particles.size());
375 // calculate count and offset
376 int prefix = 0;
377 // calculate prefix for write of the current process
378 BOOST_MPI_CHECK_RESULT(MPI_Exscan,
379 (&n_part_local, &prefix, 1, MPI_INT, MPI_SUM, m_comm));
380
381 auto const n_part_global =
382 boost::mpi::all_reduce(m_comm, n_part_local, std::plus<int>());
383
384 write_td_particle_property<2>(
385 prefix, n_part_global, particles, datasets["particles/atoms/id/value"],
386 [](auto const &p) { return Utils::Vector<int, 1>{p.id()}; });
387
388 {
389 h5xx::dataset &dataset = datasets["particles/atoms/id/value"];
390 auto const extents = static_cast<h5xx::dataspace>(dataset).extents();
392 datasets["particles/atoms/id/time"], Vector1hs{1},
393 Vector1hs{extents[0]}, Vector1hs{1});
395 datasets["particles/atoms/id/step"], Vector1hs{1},
396 Vector1hs{extents[0]}, Vector1hs{1});
397 }
398
399 if (m_fields & H5MD_OUT_TYPE) {
400 write_td_particle_property<2>(
401 prefix, n_part_global, particles,
402 datasets["particles/atoms/species/value"],
403 [](auto const &p) { return Utils::Vector<int, 1>{p.type()}; });
404 }
405 if (m_fields & H5MD_OUT_MASS) {
406 write_td_particle_property<2>(
407 prefix, n_part_global, particles,
408 datasets["particles/atoms/mass/value"],
409 [](auto const &p) { return Utils::Vector<double, 1>{p.mass()}; });
410 }
411 if (m_fields & H5MD_OUT_POS) {
412 write_td_particle_property<3>(
413 prefix, n_part_global, particles,
414 datasets["particles/atoms/position/value"],
415 [&](auto const &p) { return box_geo.folded_position(p.pos()); });
416 }
417 if (m_fields & H5MD_OUT_IMG) {
418 write_td_particle_property<3>(
419 prefix, n_part_global, particles,
420 datasets["particles/atoms/image/value"], [&](auto const &p) {
421 return box_geo.folded_image_box(p.pos(), p.image_box());
422 });
423 }
424 if (m_fields & H5MD_OUT_VEL) {
425 write_td_particle_property<3>(prefix, n_part_global, particles,
426 datasets["particles/atoms/velocity/value"],
427 [](auto const &p) { return p.v(); });
428 }
429 if (m_fields & H5MD_OUT_FORCE) {
430 write_td_particle_property<3>(prefix, n_part_global, particles,
431 datasets["particles/atoms/force/value"],
432 [](auto const &p) { return p.force(); });
433 }
434 if (m_fields & H5MD_OUT_CHARGE) {
435 write_td_particle_property<2>(
436 prefix, n_part_global, particles,
437 datasets["particles/atoms/charge/value"],
438 [](auto const &p) { return Utils::Vector<double, 1>{p.q()}; });
439 }
440 if (m_fields & H5MD_OUT_BONDS) {
441 write_connectivity(particles);
442 }
443}
444
445void File::write_connectivity(const ParticleRange &particles) {
446 MultiArray3i bond(boost::extents[0][0][0]);
447 for (auto const &p : particles) {
448 auto nbonds_local = static_cast<decltype(bond)::index>(bond.shape()[1]);
449 for (auto const b : p.bonds()) {
450 auto const partner_ids = b.partner_ids();
451 if (partner_ids.size() == 1u) {
452 bond.resize(boost::extents[1][nbonds_local + 1][2]);
453 bond[0][nbonds_local][0] = p.id();
454 bond[0][nbonds_local][1] = partner_ids[0];
455 nbonds_local++;
456 }
457 }
458 }
459
460 auto const n_bonds_local = static_cast<int>(bond.shape()[1]);
461 int prefix_bonds = 0;
462 BOOST_MPI_CHECK_RESULT(
463 MPI_Exscan, (&n_bonds_local, &prefix_bonds, 1, MPI_INT, MPI_SUM, m_comm));
464 auto const n_bonds_total =
465 boost::mpi::all_reduce(m_comm, n_bonds_local, std::plus<int>());
466 auto const extents =
467 static_cast<h5xx::dataspace>(datasets["connectivity/atoms/value"])
468 .extents();
469 Vector3hs offset_bonds = {extents[0], static_cast<hsize_t>(prefix_bonds), 0};
470 Vector3hs count_bonds = {1, static_cast<hsize_t>(n_bonds_local), 2};
471 auto const n_bond_diff =
472 std::max(static_cast<hsize_t>(n_bonds_total), extents[1]) - extents[1];
473 Vector3hs change_extent_bonds = {1, static_cast<hsize_t>(n_bond_diff), 0};
474 write_dataset(bond, datasets["connectivity/atoms/value"], change_extent_bonds,
475 offset_bonds, count_bonds);
476}
477
478void File::flush() { m_h5md_file.flush(); }
479
480} /* namespace H5md */
481} /* namespace Writer */
Vector implementation and trait types for boost qvm interoperability.
auto folded_position(Utils::Vector3d const &pos) const
Calculate coordinates folded to primary simulation box.
Utils::Vector3d const & length() const
Box length.
LeesEdwardsBC const & lees_edwards_bc() const
auto folded_image_box(Utils::Vector3d const &pos, Utils::Vector3i const &image_box) const
Calculate image box of coordinates folded to primary simulation box.
A range of particles.
base_type::size_type size() const
void write(const ParticleRange &particles, double time, int step, BoxGeometry const &geometry)
Write data to the hdf5 file.
auto const & length_unit() const
Retrieve the set length unit.
auto const & time_unit() const
Retrieve the set time unit.
auto file_path() const
Retrieve the path to the hdf5 file.
void close()
Method to perform the renaming of the temporary file from "filename" + ".bak" to "filename".
auto const & force_unit() const
Retrieve the set force unit.
auto const & mass_unit() const
Retrieve the set mass unit.
auto const & charge_unit() const
Retrieve the set charge unit.
auto const & velocity_unit() const
Retrieve the set velocity unit.
auto const & script_path() const
Retrieve the path to the simulation script.
void flush()
Method to enforce flushing the buffer to disk.
static void write_script(std::string const &target, boost::filesystem::path const &script_path)
Definition h5md_core.cpp:89
boost::multi_array< int, 3 > MultiArray3i
Definition h5md_core.cpp:49
Utils::Vector< hsize_t, 3 > Vector3hs
Definition h5md_core.cpp:52
static void write_le_dir(LeesEdwardsBC const &lebc, h5xx::dataset &dataset)
static void write_le_normal(LeesEdwardsBC const &lebc, h5xx::dataset &dataset)
static void backup_file(const std::string &from, const std::string &to)
Definition h5md_core.cpp:54
static void write_dataset(value_type const &data, h5xx::dataset &dataset, extent_type const &change_extent, extent_type const &offset, extent_type const &count)
Definition h5md_core.cpp:81
static std::vector< hsize_t > create_dims(hsize_t rank, hsize_t data_dim)
static void write_box(BoxGeometry const &box_geo, h5xx::dataset &dataset)
Utils::Vector< hsize_t, 2 > Vector2hs
Definition h5md_core.cpp:51
static void write_attributes(h5xx::file &h5md_file)
static void extend_dataset(h5xx::dataset &dataset, extent_type const &change_extent)
Definition h5md_core.cpp:69
static void write_le_off(LeesEdwardsBC const &lebc, h5xx::dataset &dataset)
static std::vector< hsize_t > create_chunk_dims(hsize_t rank, hsize_t data_dim)
void write_td_particle_property(hsize_t prefix, hsize_t n_part_global, ParticleRange const &particles, h5xx::dataset &dataset, Op op)
unsigned int shear_plane_normal
unsigned int shear_direction
bool is_compliant(std::string const &filename) const