diff --git a/src/metatensor/import.sh b/src/metatensor/import.sh new file mode 100755 index 0000000000..a1a80f9600 --- /dev/null +++ b/src/metatensor/import.sh @@ -0,0 +1,42 @@ +#!/usr/bin/env bash + +path=$1 + +if (($# != 1)) ; then +echo "usage: $0 /path/to/vesin" +echo +echo "All the commands to run:" +echo " cd /some/dir/" +echo " git clone https://github.com/luthaf/vesin" +echo " cd /your/plumed2/src/metatensor" +echo " ./import.sh /some/dir/vesin" +exit 0 +fi + +bash -c "$path/create-single-cpp.py > /dev/null" + +cp $path/include/vesin.h vesin.h +mv vesin-single-build.cpp vesin.cpp + +# Patch files to follow PLUMED linter +sed 's|#define VESIN_H|#define VESIN_H\n/*INDENT-OFF*/\n| + s|VESIN_H|__PLUMED_metatensor_vesin_h|g + s||| + s||| + s|extern "C" {|namespace PLMD {\nnamespace metatensor {\nnamespace vesin {\nextern "C" {| + s|} // extern "C"|} // extern "C"\n} // namespace vesin\n} // namespace metatensor\n} // namespace PLMD| + ' vesin.h > tmp +mv tmp vesin.h + +sed '1 s|^|/*INDENT-OFF*/\n#include "vesin.h"\n| + s||| + s||| + s|vesin::|PLMD::metatensor::vesin::| + s|namespace vesin {|namespace PLMD {\nnamespace metatensor {\nnamespace vesin {| + s|} // namespace vesin|} // namespace vesin\n} // namespace metatensor\n} // namespace PLMD| + s|using namespace PLMD::metatensor::vesin::cpu;|using namespace PLMD::metatensor::vesin;\nusing namespace PLMD::metatensor::vesin::cpu;| + ' vesin.cpp > tmp +mv tmp vesin.cpp + +cd .. +./header.sh metatensor diff --git a/src/metatensor/metatensor.cpp b/src/metatensor/metatensor.cpp index 081f851d0a..7052d38f08 100644 --- a/src/metatensor/metatensor.cpp +++ b/src/metatensor/metatensor.cpp @@ -672,15 +672,15 @@ metatensor_torch::TorchTensorBlock MetatensorPlumedAction::computeNeighbors( // use https://github.com/Luthaf/vesin to compute the requested neighbor // lists since we can not get these from PLUMED - VesinOptions options; + vesin::VesinOptions options; options.cutoff = cutoff; options.full = request->full_list(); options.return_shifts = true; options.return_distances = false; options.return_vectors = true; - VesinNeighborList* vesin_neighbor_list = new VesinNeighborList(); - memset(vesin_neighbor_list, 0, sizeof(VesinNeighborList)); + vesin::VesinNeighborList* vesin_neighbor_list = new vesin::VesinNeighborList(); + memset(vesin_neighbor_list, 0, sizeof(vesin::VesinNeighborList)); const char* error_message = NULL; int status = vesin_neighbors( @@ -688,7 +688,7 @@ metatensor_torch::TorchTensorBlock MetatensorPlumedAction::computeNeighbors( positions.size(), reinterpret_cast(&cell(0, 0)), !non_periodic, - VesinCPU, + vesin::VesinCPU, options, vesin_neighbor_list, &error_message diff --git a/src/metatensor/vesin_single_build.cpp b/src/metatensor/vesin.cpp similarity index 87% rename from src/metatensor/vesin_single_build.cpp rename to src/metatensor/vesin.cpp index caf259ba35..a03f950e77 100644 --- a/src/metatensor/vesin_single_build.cpp +++ b/src/metatensor/vesin.cpp @@ -1,4 +1,27 @@ +/* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ +Copyright (c) 2024 The METATENSOR code team +(see the PEOPLE-METATENSOR file at the root of this folder for a list of names) + +See https://docs.metatensor.org/latest/ for more information about the +metatensor package that this module allows you to call from PLUMED. + +This file is part of METATENSOR-PLUMED module. + +The METATENSOR-PLUMED module is free software: you can redistribute it and/or modify +it under the terms of the GNU Lesser General Public License as published by +the Free Software Foundation, either version 3 of the License, or +(at your option) any later version. + +The METATENSOR-PLUMED module is distributed in the hope that it will be useful, +but WITHOUT ANY WARRANTY; without even the implied warranty of +MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +GNU Lesser General Public License for more details. + +You should have received a copy of the GNU Lesser General Public License +along with the METATENSOR-PLUMED module. If not, see . ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ */ /*INDENT-OFF*/ +#include "vesin.h" #include #include #include @@ -24,6 +47,8 @@ #include #include +namespace PLMD { +namespace metatensor { namespace vesin { struct Vector; @@ -129,15 +154,19 @@ inline Vector operator*(Vector vector, Matrix matrix) { }; } -} +} // namespace vesin +} // namespace metatensor +} // namespace PLMD #endif +namespace PLMD { +namespace metatensor { namespace vesin { class BoundingBox { public: - BoundingBox(Matrix matrix, bool periodic): matrix_(std::move(matrix)), periodic_(periodic) { + BoundingBox(Matrix matrix, bool periodic): matrix_(matrix), periodic_(periodic) { if (periodic) { this->inverse_ = matrix_.inverse(); } else { @@ -228,11 +257,14 @@ inline CellShift operator-(CellShift a, CellShift b) { } - -} +} // namespace vesin +} // namespace metatensor +} // namespace PLMD #endif +namespace PLMD { +namespace metatensor { namespace vesin { namespace cpu { void free_neighbors(VesinNeighborList& neighbors); @@ -305,9 +337,9 @@ class GrowableNeighborList { } void set_pair(size_t index, size_t first, size_t second); - void set_shift(size_t index, vesin::CellShift shift); + void set_shift(size_t index, PLMD::metatensor::vesin::CellShift shift); void set_distance(size_t index, double distance); - void set_vector(size_t index, vesin::Vector vector); + void set_vector(size_t index, PLMD::metatensor::vesin::Vector vector); // reset length to 0, and allocate/deallocate members of // `neighbors` according to `options` @@ -317,13 +349,17 @@ class GrowableNeighborList { void grow(); }; -}} +} // namespace vesin +} // namespace metatensor +} // namespace PLMD +} // namespace cpu #endif -using namespace vesin::cpu; +using namespace PLMD::metatensor::vesin; +using namespace PLMD::metatensor::vesin::cpu; -void vesin::cpu::neighbors( +void PLMD::metatensor::vesin::cpu::neighbors( const Vector* points, size_t n_points, BoundingBox cell, @@ -420,7 +456,7 @@ static std::tuple divmod(int32_t a, size_t b) { auto quotient = a / b_32; auto remainder = a % b_32; if (remainder < 0) { - remainder += b; + remainder += b_32; quotient -= 1; } return std::make_tuple(quotient, remainder); @@ -441,7 +477,7 @@ divmod(std::array a, std::array b) { CellList::CellList(BoundingBox box, double cutoff): n_search_({0, 0, 0}), cells_shape_({0, 0, 0}), - box_(std::move(box)) + box_(box) { auto distances_between_faces = box_.distances_between_faces(); @@ -501,9 +537,9 @@ void CellList::add_point(size_t index, Vector position) { // find the cell in which this atom should go auto cell_index = std::array{ - static_cast(std::floor(fractional[0] * cells_shape_[0])), - static_cast(std::floor(fractional[1] * cells_shape_[1])), - static_cast(std::floor(fractional[2] * cells_shape_[2])), + static_cast(std::floor(fractional[0] * static_cast(cells_shape_[0]))), + static_cast(std::floor(fractional[1] * static_cast(cells_shape_[1]))), + static_cast(std::floor(fractional[2] * static_cast(cells_shape_[2]))), }; // deal with pbc by wrapping the atom inside if it was outside of the @@ -512,8 +548,8 @@ void CellList::add_point(size_t index, Vector position) { // auto (shift, cell_index) = if (box_.periodic()) { auto result = divmod(cell_index, cells_shape_); - shift = CellShift{std::move(std::get<0>(result))}; - cell_index = std::move(std::get<1>(result)); + shift = CellShift{std::get<0>(result)}; + cell_index = std::get<1>(result); } else { shift = CellShift({0, 0, 0}); cell_index = std::array{ @@ -532,7 +568,7 @@ void CellList::foreach_pair(Function callback) { for (int32_t cell_i_x=0; cell_i_x(cells_shape_[0]); cell_i_x++) { for (int32_t cell_i_y=0; cell_i_y(cells_shape_[1]); cell_i_y++) { for (int32_t cell_i_z=0; cell_i_z(cells_shape_[2]); cell_i_z++) { - auto& current_cell = this->get_cell({cell_i_x, cell_i_y, cell_i_z}); + const auto& current_cell = this->get_cell({cell_i_x, cell_i_y, cell_i_z}); // look through each neighboring cell for (int32_t delta_x=-n_search_[0]; delta_x<=n_search_[0]; delta_x++) { for (int32_t delta_y=-n_search_[1]; delta_y<=n_search_[1]; delta_y++) { @@ -590,7 +626,7 @@ void GrowableNeighborList::set_pair(size_t index, size_t first, size_t second) { this->neighbors.pairs[index][1] = second; } -void GrowableNeighborList::set_shift(size_t index, vesin::CellShift shift) { +void GrowableNeighborList::set_shift(size_t index, PLMD::metatensor::vesin::CellShift shift) { if (index >= this->capacity) { this->grow(); } @@ -608,7 +644,7 @@ void GrowableNeighborList::set_distance(size_t index, double distance) { this->neighbors.distances[index] = distance; } -void GrowableNeighborList::set_vector(size_t index, vesin::Vector vector) { +void GrowableNeighborList::set_vector(size_t index, PLMD::metatensor::vesin::Vector vector) { if (index >= this->capacity) { this->grow(); } @@ -620,30 +656,30 @@ void GrowableNeighborList::set_vector(size_t index, vesin::Vector vector) { template static scalar_t (*alloc(scalar_t (*ptr)[N], size_t size, size_t new_size))[N] { - ptr = reinterpret_cast(std::realloc(ptr, new_size * sizeof(scalar_t[N]))); + auto* new_ptr = reinterpret_cast(std::realloc(ptr, new_size * sizeof(scalar_t[N]))); - if (ptr == nullptr) { + if (new_ptr == nullptr) { throw std::bad_alloc(); } // initialize with a bit pattern that maps to NaN for double - std::memset(ptr + size, 0b11111111, (new_size - size) * sizeof(scalar_t[N])); + std::memset(new_ptr + size, 0b11111111, (new_size - size) * sizeof(scalar_t[N])); - return ptr; + return new_ptr; } template static scalar_t* alloc(scalar_t* ptr, size_t size, size_t new_size) { - ptr = reinterpret_cast(std::realloc(ptr, new_size * sizeof(scalar_t))); + auto* new_ptr = reinterpret_cast(std::realloc(ptr, new_size * sizeof(scalar_t))); - if (ptr == nullptr) { + if (new_ptr == nullptr) { throw std::bad_alloc(); } // initialize with a bit pattern that maps to NaN for double - std::memset(ptr + size, 0b11111111, (new_size - size) * sizeof(scalar_t)); + std::memset(new_ptr + size, 0b11111111, (new_size - size) * sizeof(scalar_t)); - return ptr; + return new_ptr; } void GrowableNeighborList::grow() { @@ -652,7 +688,7 @@ void GrowableNeighborList::grow() { new_size = 1; } - auto new_pairs = alloc(neighbors.pairs, neighbors.length, new_size); + auto* new_pairs = alloc(neighbors.pairs, neighbors.length, new_size); int32_t (*new_shifts)[3] = nullptr; if (options.return_shifts) { @@ -698,7 +734,7 @@ void GrowableNeighborList::reset() { this->neighbors.length = 0; // allocate/deallocate pointers as required - auto shifts = this->neighbors.shifts; + auto* shifts = this->neighbors.shifts; if (this->options.return_shifts && shifts == nullptr) { shifts = alloc(shifts, 0, capacity); } else if (!this->options.return_shifts && shifts != nullptr) { @@ -706,7 +742,7 @@ void GrowableNeighborList::reset() { shifts = nullptr; } - auto distances = this->neighbors.distances; + auto* distances = this->neighbors.distances; if (this->options.return_distances && distances == nullptr) { distances = alloc(distances, 0, capacity); } else if (!this->options.return_distances && distances != nullptr) { @@ -714,7 +750,7 @@ void GrowableNeighborList::reset() { distances = nullptr; } - auto vectors = this->neighbors.vectors; + auto* vectors = this->neighbors.vectors; if (this->options.return_vectors && vectors == nullptr) { vectors = alloc(vectors, 0, capacity); } else if (!this->options.return_vectors && vectors != nullptr) { @@ -728,7 +764,7 @@ void GrowableNeighborList::reset() { } -void vesin::cpu::free_neighbors(VesinNeighborList& neighbors) { +void PLMD::metatensor::vesin::cpu::free_neighbors(VesinNeighborList& neighbors) { assert(neighbors.device == VesinCPU); std::free(neighbors.pairs); @@ -792,16 +828,16 @@ extern "C" int vesin_neighbors( try { if (device == VesinCPU) { - auto matrix = vesin::Matrix{{{ + auto matrix = PLMD::metatensor::vesin::Matrix{{{ {{box[0][0], box[0][1], box[0][2]}}, {{box[1][0], box[1][1], box[1][2]}}, {{box[2][0], box[2][1], box[2][2]}}, }}}; - vesin::cpu::neighbors( - reinterpret_cast(points), + PLMD::metatensor::vesin::cpu::neighbors( + reinterpret_cast(points), n_points, - vesin::BoundingBox(matrix, periodic), + PLMD::metatensor::vesin::BoundingBox(matrix, periodic), options, *neighbors ); @@ -833,7 +869,7 @@ extern "C" void vesin_free(VesinNeighborList* neighbors) { if (neighbors->device == VesinUnknownDevice) { // nothing to do } else if (neighbors->device == VesinCPU) { - vesin::cpu::free_neighbors(*neighbors); + PLMD::metatensor::vesin::cpu::free_neighbors(*neighbors); } std::memset(neighbors, 0, sizeof(VesinNeighborList)); diff --git a/src/metatensor/vesin.h b/src/metatensor/vesin.h index 6efaaea523..e02231d4ce 100644 --- a/src/metatensor/vesin.h +++ b/src/metatensor/vesin.h @@ -1,20 +1,68 @@ +/* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ +Copyright (c) 2024 The METATENSOR code team +(see the PEOPLE-METATENSOR file at the root of this folder for a list of names) + +See https://docs.metatensor.org/latest/ for more information about the +metatensor package that this module allows you to call from PLUMED. + +This file is part of METATENSOR-PLUMED module. + +The METATENSOR-PLUMED module is free software: you can redistribute it and/or modify +it under the terms of the GNU Lesser General Public License as published by +the Free Software Foundation, either version 3 of the License, or +(at your option) any later version. + +The METATENSOR-PLUMED module is distributed in the hope that it will be useful, +but WITHOUT ANY WARRANTY; without even the implied warranty of +MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +GNU Lesser General Public License for more details. + +You should have received a copy of the GNU Lesser General Public License +along with the METATENSOR-PLUMED module. If not, see . ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ */ +#ifndef __PLUMED_metatensor_vesin_h +#define __PLUMED_metatensor_vesin_h /*INDENT-OFF*/ -#ifndef VESIN_H -#define VESIN_H -#include -#include + +#include +#include + +#if defined(VESIN_SHARED) + #if defined(VESIN_EXPORTS) + #if defined(__clang__) || defined(__GNUC__) || defined(__GNUG__) + #define VESIN_API __attribute__((visibility("default"))) + #elif defined(_MSC_VER) + #define VESIN_API __declspec(dllexport) + #else + #define VESIN_API + #endif + #else + #if defined(__clang__) || defined(__GNUC__) || defined(__GNUG__) + #define VESIN_API __attribute__((visibility("default"))) + #elif defined(_MSC_VER) + #define VESIN_API __declspec(dllimport) + #else + #define VESIN_API + #endif + #endif +#else + #define VESIN_API +#endif #ifdef __cplusplus +namespace PLMD { +namespace metatensor { +namespace vesin { extern "C" { #endif /// Options for a neighbor list calculation -typedef struct VesinOptions { +struct VesinOptions { /// Spherical cutoff, only pairs below this cutoff will be included double cutoff; /// Should the returned neighbor list be a full list (include both `i -> j` - /// and `j -> i` pairs) or a half list (include only `i -> j`). + /// and `j -> i` pairs) or a half list (include only `i -> j`)? bool full; // TODO: sort option? @@ -24,7 +72,7 @@ typedef struct VesinOptions { bool return_distances; /// Should the returned `VesinNeighborList` contain `vector`? bool return_vectors; -} VesinOptions; +}; /// Device on which the data can be enum VesinDevice { @@ -55,7 +103,7 @@ enum VesinDevice { /// /// Under periodic boundary conditions, two atoms can be part of multiple pairs, /// each pair having a different periodic shift. -typedef struct VesinNeighborList { +struct VESIN_API VesinNeighborList { #ifdef __cplusplus VesinNeighborList(): length(0), @@ -87,11 +135,11 @@ typedef struct VesinNeighborList { double (*vectors)[3]; // TODO: custom memory allocators? -} VesinNeighborList; +}; /// Free all allocated memory inside a `VesinNeighborList`, according the it's /// `device`. -void vesin_free(VesinNeighborList* neighbors); +void VESIN_API vesin_free(struct VesinNeighborList* neighbors); /// Compute a neighbor list. /// @@ -114,14 +162,14 @@ void vesin_free(VesinNeighborList* neighbors); /// @param error_message Pointer to a `char*` that wil be set to the error /// message if this function fails. This does not need to be freed when no /// longer needed. -int vesin_neighbors( +int VESIN_API vesin_neighbors( const double (*points)[3], size_t n_points, const double box[3][3], bool periodic, VesinDevice device, - VesinOptions options, - VesinNeighborList* neighbors, + struct VesinOptions options, + struct VesinNeighborList* neighbors, const char** error_message ); @@ -129,6 +177,9 @@ int vesin_neighbors( #ifdef __cplusplus } // extern "C" +} // namespace vesin +} // namespace metatensor +} // namespace PLMD #endif