diff --git a/metatomic-core/CMakeLists.txt b/metatomic-core/CMakeLists.txt
index 717e52e81..5cb30446a 100644
--- a/metatomic-core/CMakeLists.txt
+++ b/metatomic-core/CMakeLists.txt
@@ -5,6 +5,11 @@
 # an easier to use, idiomatic Rust API.
 cmake_minimum_required(VERSION 3.22)
 
+if (POLICY CMP0135)
+    # Use download time as timestamp when extracting files from archives.
+    cmake_policy(SET CMP0135 NEW)
+endif()
+
 # Is metatomic the main project configured by the user? Or is this being used
 # as a submodule/subdirectory?
 if (${CMAKE_CURRENT_SOURCE_DIR} STREQUAL ${CMAKE_SOURCE_DIR})
@@ -449,6 +454,17 @@ else()
     target_link_libraries(metatomic::shared INTERFACE metatensor)
 endif()
 
+include(FetchContent)
+
+# JSON library from https://github.com/nlohmann/json
+FetchContent_Declare(nlohmann_json
+    URL https://github.com/nlohmann/json/releases/download/v3.11.3/json.tar.xz
+    URL_HASH SHA256=d6c65aca6b1ed68e7a182f4757257b107ae403032760ed6ef121c9d55e81757d
+)
+FetchContent_MakeAvailable(nlohmann_json)
+
+target_link_libraries(metatomic::shared INTERFACE nlohmann_json::nlohmann_json)
+target_link_libraries(metatomic::static INTERFACE nlohmann_json::nlohmann_json)
 
 if (BUILD_SHARED_LIBS)
     add_library(metatomic ALIAS metatomic::shared)
diff --git a/metatomic-core/cmake/metatomic-config.in.cmake b/metatomic-core/cmake/metatomic-config.in.cmake
index 90fca167a..8f2fc5443 100644
--- a/metatomic-core/cmake/metatomic-config.in.cmake
+++ b/metatomic-core/cmake/metatomic-config.in.cmake
@@ -4,6 +4,7 @@ cmake_minimum_required(VERSION 3.22)
 
 include(CMakeFindDependencyMacro)
 include(FindPackageHandleStandardArgs)
+include(FetchContent)
 
 if(metatomic_FOUND)
     return()
@@ -15,6 +16,15 @@ enable_language(CXX)
 set(REQUIRED_METATENSOR_VERSION @REQUIRED_METATENSOR_VERSION@)
 find_package(metatensor ${REQUIRED_METATENSOR_VERSION} CONFIG REQUIRED)
 
+if (NOT TARGET nlohmann_json::nlohmann_json)
+    # JSON library from https://github.com/nlohmann/json
+    FetchContent_Declare(nlohmann_json
+        URL https://github.com/nlohmann/json/releases/download/v3.11.3/json.tar.xz
+        URL_HASH SHA256=d6c65aca6b1ed68e7a182f4757257b107ae403032760ed6ef121c9d55e81757d
+    )
+    FetchContent_MakeAvailable(nlohmann_json)
+endif()
+
 get_filename_component(METATOMIC_PREFIX_DIR "${CMAKE_CURRENT_LIST_DIR}/@PACKAGE_RELATIVE_PATH@" ABSOLUTE)
 
 if (WIN32)
@@ -46,7 +56,7 @@ if (@METATOMIC_INSTALL_BOTH_STATIC_SHARED@ OR @BUILD_SHARED_LIBS@)
     )
 
     target_compile_features(metatomic::shared INTERFACE cxx_std_17)
-    target_link_libraries(metatomic::shared INTERFACE metatensor)
+    target_link_libraries(metatomic::shared INTERFACE metatensor nlohmann_json::nlohmann_json)
 
     if (WIN32)
         if (NOT EXISTS ${METATOMIC_IMPLIB_LOCATION})
@@ -75,7 +85,7 @@ if (@METATOMIC_INSTALL_BOTH_STATIC_SHARED@ OR NOT @BUILD_SHARED_LIBS@)
     )
 
     target_compile_features(metatomic::static INTERFACE cxx_std_17)
-    target_link_libraries(metatomic::static INTERFACE metatensor)
+    target_link_libraries(metatomic::static INTERFACE metatensor nlohmann_json::nlohmann_json)
 endif()
 
 # Export either the shared or static library as the metatomic target
diff --git a/metatomic-core/include/metatomic.hpp b/metatomic-core/include/metatomic.hpp
index 3b5c8ac2a..9eaad1e57 100644
--- a/metatomic-core/include/metatomic.hpp
+++ b/metatomic-core/include/metatomic.hpp
@@ -1,4 +1,5 @@
 #include "metatomic/utils.hpp"    // IWYU pragma: export
+#include "metatomic/metadata.hpp" // IWYU pragma: export
 #include "metatomic/system.hpp"   // IWYU pragma: export
 #include "metatomic/model.hpp"    // IWYU pragma: export
 #include "metatomic/plugin.hpp"   // IWYU pragma: export
diff --git a/metatomic-core/include/metatomic/metadata.hpp b/metatomic-core/include/metatomic/metadata.hpp
new file mode 100644
index 000000000..35a92c75d
--- /dev/null
+++ b/metatomic-core/include/metatomic/metadata.hpp
@@ -0,0 +1,192 @@
+#pragma once
+
+#include <cstdint>
+#include <cstring>
+
+#include <map>
+#include <sstream>
+#include <string>
+#include <utility>
+#include <vector>
+
+#include <nlohmann/json.hpp>
+
+namespace metatomic {
+
+/// Options for the calculation of a pair list.
+class PairListOptions final {
+public:
+    PairListOptions() = default;
+
+    PairListOptions(double cutoff_value, bool full_list_value, bool strict_value, std::vector<std::string> requestors_list = {}):
+        cutoff(cutoff_value),
+        full_list(full_list_value),
+        strict(strict_value),
+        requestors(std::move(requestors_list))
+    {}
+
+    double cutoff = 0.0;
+    bool full_list = false;
+    bool strict = false;
+    std::vector<std::string> requestors;
+
+    std::string to_json() const;
+    static PairListOptions from_json(const std::string& json);
+};
+
+/// Metadata about a specific exported model.
+class ModelMetadata final {
+public:
+    ModelMetadata() = default;
+
+    ModelMetadata(
+        std::string model_name,
+        std::string model_description,
+        std::vector<std::string> model_authors,
+        std::map<std::string, std::vector<std::string>> model_references = {},
+        std::map<std::string, std::string> extra_metadata = {}
+    ):
+        name(std::move(model_name)),
+        description(std::move(model_description)),
+        authors(std::move(model_authors)),
+        references(std::move(model_references)),
+        extra(std::move(extra_metadata))
+    {}
+
+    std::string name;
+    std::string description;
+    std::vector<std::string> authors;
+    std::map<std::string, std::vector<std::string>> references;
+    std::map<std::string, std::string> extra;
+
+    std::string to_json() const;
+    static ModelMetadata from_json(const std::string& json);
+};
+
+/// Description of a quantity used as model input or output.
+class Quantity final {
+public:
+    Quantity() = default;
+
+    Quantity(
+        std::string quantity_name,
+        std::string quantity_unit,
+        std::vector<std::string> quantity_gradients,
+        std::string quantity_sample_kind
+    ):
+        name(std::move(quantity_name)),
+        unit(std::move(quantity_unit)),
+        gradients(std::move(quantity_gradients)),
+        sample_kind(std::move(quantity_sample_kind))
+    {}
+
+    std::string name;
+    std::string unit;
+    std::vector<std::string> gradients;
+    std::string sample_kind;
+
+    std::string to_json() const;
+    static Quantity from_json(const std::string& json);
+};
+
+namespace details {
+    inline std::string double_to_hex(double value) {
+        uint64_t bits = 0;
+        static_assert(sizeof(bits) == sizeof(value), "unexpected double size");
+        std::memcpy(&bits, &value, sizeof(bits));
+
+        auto stream = std::ostringstream();
+        stream << "0x" << std::hex << bits;
+        return stream.str();
+    }
+
+    inline double hex_to_double(const std::string& value) {
+        auto bits = uint64_t(0);
+        auto stream = std::istringstream(value);
+        if (value.rfind("0x", 0) == 0 || value.rfind("0X", 0) == 0) {
+            stream.seekg(2);
+        }
+        stream >> std::hex >> bits;
+
+        auto result = 0.0;
+        static_assert(sizeof(bits) == sizeof(result), "unexpected double size");
+        std::memcpy(&result, &bits, sizeof(result));
+        return result;
+    }
+
+} // namespace details
+
+inline std::string PairListOptions::to_json() const {
+    return nlohmann::json{
+        {"type", "metatomic_pair_options"},
+        {"cutoff", details::double_to_hex(cutoff)},
+        {"full_list", full_list},
+        {"strict", strict},
+        {"requestors", requestors},
+    }.dump();
+}
+
+inline PairListOptions PairListOptions::from_json(const std::string& string) {
+    auto json = nlohmann::json::parse(string);
+    auto options = PairListOptions();
+
+    auto cutoff = std::string("0x0");
+    if (json.contains("cutoff")) {
+        cutoff = json.at("cutoff").get<std::string>();
+    }
+
+    options.cutoff = details::hex_to_double(cutoff);
+    options.full_list = json.value("full_list", false);
+    options.strict = json.value("strict", false);
+    options.requestors = json.value("requestors", std::vector<std::string>{});
+
+    return options;
+}
+
+inline std::string ModelMetadata::to_json() const {
+    return nlohmann::json{
+        {"type", "metatomic_model_metadata"},
+        {"name", name},
+        {"description", description},
+        {"authors", authors},
+        {"references", references},
+        {"extra", extra},
+    }.dump();
+}
+
+inline ModelMetadata ModelMetadata::from_json(const std::string& string) {
+    auto json = nlohmann::json::parse(string);
+    auto metadata = ModelMetadata();
+
+    metadata.name = json.value("name", "");
+    metadata.description = json.value("description", "");
+    metadata.authors = json.value("authors", std::vector<std::string>{});
+    metadata.references = json.value("references", std::map<std::string, std::vector<std::string>>{});
+    metadata.extra = json.value("extra", std::map<std::string, std::string>{});
+
+    return metadata;
+}
+
+inline std::string Quantity::to_json() const {
+    return nlohmann::json{
+        {"type", "metatomic_quantity"},
+        {"name", name},
+        {"unit", unit},
+        {"gradients", gradients},
+        {"sample_kind", sample_kind},
+    }.dump();
+}
+
+inline Quantity Quantity::from_json(const std::string& string) {
+    auto json = nlohmann::json::parse(string);
+    auto quantity = Quantity();
+
+    quantity.name = json.value("name", "");
+    quantity.unit = json.value("unit", "");
+    quantity.gradients = json.value("gradients", std::vector<std::string>{});
+    quantity.sample_kind = json.value("sample_kind", "");
+
+    return quantity;
+}
+
+} // namespace metatomic
diff --git a/metatomic-core/include/metatomic/model.hpp b/metatomic-core/include/metatomic/model.hpp
index 1cae91bdf..ee140760d 100644
--- a/metatomic-core/include/metatomic/model.hpp
+++ b/metatomic-core/include/metatomic/model.hpp
@@ -1,7 +1,410 @@
 #pragma once
 
+#include <string>
+#include <memory>
+#include <utility>
+#include <vector>
+
 #include <metatomic.h>
+#include <metatensor.hpp>
+#include <nlohmann/json.hpp>
+
+#include "./metadata.hpp"
+#include "./system.hpp"
+#include "./utils.hpp"
 
 namespace metatomic {
 
+/// Abstract base class for atomistic models implemented in C++.
+class ModelInterface {
+public:
+    virtual ~ModelInterface() = default;
+
+    /// Get metadata about this model.
+    virtual ModelMetadata metadata() const = 0;
+
+    /// Get all quantities this model can compute.
+    virtual std::vector<Quantity> supported_outputs() const = 0;
+
+    /// Get all pair lists this model requires.
+    virtual std::vector<PairListOptions> requested_pair_lists() const {
+        return {};
+    }
+
+    /// Get all custom inputs this model requires.
+    virtual std::vector<Quantity> requested_inputs() const {
+        return {};
+    }
+
+    /// Execute this model.
+    virtual std::vector<metatensor::TensorMap> execute(
+        const std::vector<const System*>& systems,
+        const mts_labels_t* selected_atoms,
+        const std::vector<Quantity>& requested_outputs
+    ) = 0;
+};
+
+/// RAII wrapper around a `mta_model_t`.
+class AtomisticModel final {
+public:
+    /// Create an empty, invalid model.
+    AtomisticModel() {
+        model_ = empty_model();
+    }
+
+    /// Take ownership of a raw `mta_model_t`.
+    explicit AtomisticModel(mta_model_t model): model_(model) {}
+
+    /// Create a C API model wrapping a C++ model implementation.
+    explicit AtomisticModel(std::unique_ptr<ModelInterface> model) {
+        if (model == nullptr) {
+            throw Error("can not create a metatomic::AtomisticModel from a null ModelInterface");
+        }
+
+        model_ = empty_model();
+        model_.data = model.release();
+        model_.unload = &AtomisticModel::unload_callback;
+        model_.metadata = &AtomisticModel::metadata_callback;
+        model_.supported_outputs = &AtomisticModel::supported_outputs_callback;
+        model_.requested_pair_lists = &AtomisticModel::requested_pair_lists_callback;
+        model_.requested_inputs = &AtomisticModel::requested_inputs_callback;
+        model_.execute_inner = &AtomisticModel::execute_callback;
+    }
+
+    ~AtomisticModel() {
+        this->reset_noexcept();
+    }
+
+    AtomisticModel(const AtomisticModel&) = delete;
+    AtomisticModel& operator=(const AtomisticModel&) = delete;
+
+    AtomisticModel(AtomisticModel&& other) noexcept: AtomisticModel() {
+        *this = std::move(other);
+    }
+
+    AtomisticModel& operator=(AtomisticModel&& other) noexcept {
+        if (this != &other) {
+            this->reset_noexcept();
+            model_ = other.model_;
+            other.model_ = empty_model();
+        }
+        return *this;
+    }
+
+    /// Does this wrapper contain a model?
+    bool is_valid() const {
+        return model_.data != nullptr;
+    }
+
+    /// Unload the model.
+    void unload() {
+        if (model_.data != nullptr && model_.unload != nullptr) {
+            details::check_status(model_.unload(model_.data));
+        }
+        model_ = empty_model();
+    }
+
+    /// Get model metadata serialized as JSON.
+    std::string metadata_json() const {
+        this->check_callback(model_.metadata, "metadata");
+
+        mta_string_t metadata = nullptr;
+        details::check_status(model_.metadata(model_.data, &metadata));
+        return String(metadata).str();
+    }
+
+    /// Get model metadata.
+    ModelMetadata metadata() const {
+        return ModelMetadata::from_json(this->metadata_json());
+    }
+
+    /// Get supported outputs serialized as JSON.
+    std::string supported_outputs_json() const {
+        this->check_callback(model_.supported_outputs, "supported_outputs");
+
+        mta_string_t outputs = nullptr;
+        details::check_status(model_.supported_outputs(model_.data, &outputs));
+        return String(outputs).str();
+    }
+
+    /// Get all quantities this model can compute.
+    std::vector<Quantity> supported_outputs() const {
+        auto outputs = std::vector<Quantity>();
+        for (const auto& output: nlohmann::json::parse(this->supported_outputs_json())) {
+            outputs.push_back(Quantity::from_json(output.dump()));
+        }
+        return outputs;
+    }
+
+    /// Get all pair lists requested by this model serialized as a JSON array.
+    std::string requested_pair_lists_json() const {
+        this->check_callback(model_.requested_pair_lists, "requested_pair_lists");
+
+        mta_string_t options = nullptr;
+        details::check_status(model_.requested_pair_lists(model_.data, &options));
+        return String(options).str();
+    }
+
+    /// Get all pair lists requested by this model.
+    std::vector<PairListOptions> requested_pair_lists() const {
+        auto result = std::vector<PairListOptions>();
+        for (const auto& options: nlohmann::json::parse(this->requested_pair_lists_json())) {
+            result.push_back(PairListOptions::from_json(options.dump()));
+        }
+        return result;
+    }
+
+    /// Get all custom inputs requested by this model serialized as a JSON array.
+    std::string requested_inputs_json() const {
+        this->check_callback(model_.requested_inputs, "requested_inputs");
+
+        mta_string_t inputs = nullptr;
+        details::check_status(model_.requested_inputs(model_.data, &inputs));
+        return String(inputs).str();
+    }
+
+    /// Get all custom inputs requested by this model.
+    std::vector<Quantity> requested_inputs() const {
+        auto result = std::vector<Quantity>();
+        for (const auto& input: nlohmann::json::parse(this->requested_inputs_json())) {
+            result.push_back(Quantity::from_json(input.dump()));
+        }
+        return result;
+    }
+
+    /// Execute this model.
+    ///
+    /// The number of returned tensor maps is `requested_outputs.size()`.
+    std::vector<metatensor::TensorMap> execute(
+        const std::vector<const System*>& systems,
+        const metatensor::Labels* selected_atoms,
+        const std::vector<Quantity>& requested_outputs,
+        bool check_consistency = true
+    ) {
+        this->check_valid();
+
+        auto c_systems = std::vector<const mta_system_t*>();
+        c_systems.reserve(systems.size());
+        for (const auto* system: systems) {
+            details::check_pointer(system);
+            c_systems.push_back(system->as_mta_system_t());
+        }
+
+        auto c_requested_outputs = std::vector<const char*>();
+        auto requested_outputs_json = std::vector<std::string>();
+        requested_outputs_json.reserve(requested_outputs.size());
+        c_requested_outputs.reserve(requested_outputs.size());
+        for (const auto& output: requested_outputs) {
+            requested_outputs_json.push_back(output.to_json());
+            c_requested_outputs.push_back(requested_outputs_json.back().c_str());
+        }
+
+        auto raw_outputs = std::vector<mts_tensormap_t*>(requested_outputs.size(), nullptr);
+        details::check_status(mta_execute_model(
+            model_,
+            c_systems.data(),
+            c_systems.size(),
+            selected_atoms == nullptr ? nullptr : selected_atoms->as_mts_labels_t(),
+            c_requested_outputs.data(),
+            c_requested_outputs.size(),
+            check_consistency,
+            raw_outputs.data(),
+            raw_outputs.size()
+        ));
+
+        auto outputs = std::vector<metatensor::TensorMap>();
+        outputs.reserve(raw_outputs.size());
+
+        try {
+            for (auto*& output: raw_outputs) {
+                details::check_pointer(output);
+                outputs.emplace_back(output);
+                output = nullptr;
+            }
+        } catch (...) {
+            for (auto* output: raw_outputs) {
+                if (output != nullptr) {
+                    (void)mts_tensormap_free(output);
+                }
+            }
+            throw;
+        }
+
+        return outputs;
+    }
+
+    /// Execute this model on all atoms.
+    std::vector<metatensor::TensorMap> execute(
+        const std::vector<const System*>& systems,
+        const std::vector<Quantity>& requested_outputs,
+        bool check_consistency = true
+    ) {
+        return this->execute(systems, nullptr, requested_outputs, check_consistency);
+    }
+
+    /// Get the underlying `mta_model_t`.
+    const mta_model_t& as_mta_model_t() const & {
+        return model_;
+    }
+
+    const mta_model_t& as_mta_model_t() && = delete;
+
+    /// Release the underlying `mta_model_t` without unloading it.
+    mta_model_t release() {
+        auto model = model_;
+        model_ = empty_model();
+        return model;
+    }
+
+private:
+    static mta_model_t empty_model() {
+        mta_model_t model;
+        model.data = nullptr;
+        model.unload = nullptr;
+        model.metadata = nullptr;
+        model.supported_outputs = nullptr;
+        model.requested_pair_lists = nullptr;
+        model.requested_inputs = nullptr;
+        model.execute_inner = nullptr;
+        return model;
+    }
+
+    static ModelInterface* model_base(const void* data) {
+        details::check_pointer(data);
+        return static_cast<ModelInterface*>(const_cast<void*>(data));
+    }
+
+    static mta_status_t unload_callback(void* data) {
+        return details::catch_exceptions([&]() {
+            delete model_base(data);
+        });
+    }
+
+    static mta_status_t metadata_callback(const void* data, mta_string_t* metadata_json) {
+        return details::catch_exceptions([&]() {
+            details::check_pointer(metadata_json);
+            *metadata_json = mta_string_create(model_base(data)->metadata().to_json().c_str());
+            details::check_pointer(*metadata_json);
+        });
+    }
+
+    static mta_status_t supported_outputs_callback(const void* data, mta_string_t* outputs_json) {
+        return details::catch_exceptions([&]() {
+            details::check_pointer(outputs_json);
+            auto outputs = nlohmann::json::array();
+            for (const auto& output: model_base(data)->supported_outputs()) {
+                outputs.push_back(nlohmann::json::parse(output.to_json()));
+            }
+
+            *outputs_json = mta_string_create(outputs.dump().c_str());
+            details::check_pointer(*outputs_json);
+        });
+    }
+
+    static mta_status_t requested_pair_lists_callback(const void* data, mta_string_t* pair_options_json) {
+        return details::catch_exceptions([&]() {
+            details::check_pointer(pair_options_json);
+            auto options = nlohmann::json::array();
+            for (const auto& option: model_base(data)->requested_pair_lists()) {
+                options.push_back(nlohmann::json::parse(option.to_json()));
+            }
+
+            *pair_options_json = mta_string_create(options.dump().c_str());
+            details::check_pointer(*pair_options_json);
+        });
+    }
+
+    static mta_status_t requested_inputs_callback(const void* data, mta_string_t* inputs_json) {
+        return details::catch_exceptions([&]() {
+            details::check_pointer(inputs_json);
+            auto inputs = nlohmann::json::array();
+            for (const auto& input: model_base(data)->requested_inputs()) {
+                inputs.push_back(nlohmann::json::parse(input.to_json()));
+            }
+
+            *inputs_json = mta_string_create(inputs.dump().c_str());
+            details::check_pointer(*inputs_json);
+        });
+    }
+
+    static mta_status_t execute_callback(
+        void* data,
+        const mta_system_t* const* systems,
+        uintptr_t systems_count,
+        const mts_labels_t* selected_atoms,
+        const char* const* requested_outputs_json,
+        uintptr_t requested_outputs_count,
+        mts_tensormap_t** outputs,
+        uintptr_t outputs_count
+    ) {
+        return details::catch_exceptions([&]() {
+            if (systems_count != 0) {
+                details::check_pointer(systems);
+            }
+            if (requested_outputs_count != 0) {
+                details::check_pointer(requested_outputs_json);
+            }
+            if (outputs_count != 0) {
+                details::check_pointer(outputs);
+            }
+            if (requested_outputs_count != outputs_count) {
+                throw Error("expected one output storage slot for each requested output");
+            }
+
+            auto system_views = std::vector<System>();
+            system_views.reserve(systems_count);
+            for (uintptr_t i=0; i<systems_count; i++) {
+                details::check_pointer(systems[i]);
+                system_views.push_back(System::unsafe_view_from_ptr(systems[i]));
+            }
+
+            auto cxx_systems = std::vector<const System*>();
+            cxx_systems.reserve(system_views.size());
+            for (const auto& system: system_views) {
+                cxx_systems.push_back(&system);
+            }
+
+            auto requested_outputs = std::vector<Quantity>();
+            requested_outputs.reserve(requested_outputs_count);
+            for (uintptr_t i=0; i<requested_outputs_count; i++) {
+                details::check_pointer(requested_outputs_json[i]);
+                requested_outputs.push_back(Quantity::from_json(requested_outputs_json[i]));
+            }
+
+            auto result = model_base(data)->execute(cxx_systems, selected_atoms, requested_outputs);
+            if (result.size() != outputs_count) {
+                throw Error("model returned the wrong number of outputs");
+            }
+
+            for (uintptr_t i=0; i<outputs_count; i++) {
+                outputs[i] = mts_tensormap_copy(result[i].as_mts_tensormap_t());
+                details::check_pointer(outputs[i]);
+            }
+        });
+    }
+
+    void reset_noexcept() noexcept {
+        if (model_.data != nullptr && model_.unload != nullptr) {
+            (void)model_.unload(model_.data);
+        }
+        model_ = empty_model();
+    }
+
+    void check_valid() const {
+        if (model_.data == nullptr) {
+            throw Error("can not use an empty metatomic::AtomisticModel");
+        }
+    }
+
+    template<typename Callback>
+    void check_callback(Callback callback, const char* name) const {
+        this->check_valid();
+        if (callback == nullptr) {
+            throw Error("metatomic::AtomisticModel does not implement " + std::string(name));
+        }
+    }
+
+    mta_model_t model_;
+};
+
 } // namespace metatomic
diff --git a/metatomic-core/include/metatomic/plugin.hpp b/metatomic-core/include/metatomic/plugin.hpp
index 1cae91bdf..f5aee0692 100644
--- a/metatomic-core/include/metatomic/plugin.hpp
+++ b/metatomic-core/include/metatomic/plugin.hpp
@@ -1,7 +1,148 @@
 #pragma once
 
+#include <string>
+#include <type_traits>
+#include <utility>
+
 #include <metatomic.h>
 
+#include "./model.hpp"
+#include "./utils.hpp"
+
 namespace metatomic {
 
+/// Abstract base class for metatomic plugins implemented in C++.
+class Plugin {
+public:
+    virtual ~Plugin() = default;
+
+    /// Name used to identify this plugin.
+    virtual std::string name() const = 0;
+
+    /// Load a model from `load_from`, using the provided JSON options.
+    virtual AtomisticModel load_model(
+        const std::string& load_from,
+        const std::string& options_json = "{}"
+    ) = 0;
+};
+
+/// Handle to a plugin registered in metatomic's global plugin registry.
+class PluginHandle final {
+public:
+    explicit PluginHandle(std::string name): name_(std::move(name)) {}
+
+    /// Name used to identify this plugin.
+    const std::string& name() const {
+        return name_;
+    }
+
+    /// Load a model from `load_from`, using the provided JSON options.
+    AtomisticModel load_model(
+        const std::string& load_from,
+        const std::string& options_json = "{}"
+    ) const {
+        auto model = mta_model_t{};
+        details::check_status(mta_load_model(
+            name_.c_str(),
+            load_from.c_str(),
+            options_json.c_str(),
+            &model
+        ));
+
+        return AtomisticModel(model);
+    }
+
+private:
+    std::string name_;
+};
+
+namespace details {
+    template<typename PluginT>
+    struct PluginRegistration {
+        static PluginT* plugin;
+        static std::string name;
+
+        static mta_status_t load_model(
+            const char* load_from,
+            const char* options_json,
+            mta_model_t* model
+        ) {
+            return details::catch_exceptions([&]() {
+                details::check_pointer(plugin);
+                details::check_pointer(model);
+
+                auto loaded = plugin->load_model(
+                    load_from == nullptr ? "" : load_from,
+                    options_json == nullptr ? "{}" : options_json
+                );
+
+                *model = loaded.release();
+            });
+        }
+    };
+
+    template<typename PluginT>
+    PluginT* PluginRegistration<PluginT>::plugin = nullptr;
+
+    template<typename PluginT>
+    std::string PluginRegistration<PluginT>::name;
+} // namespace details
+
+/// Register a C++ plugin.
+///
+/// Due to the current C plugin ABI, this stores one plugin instance per concrete
+/// C++ plugin type. The registered object must outlive all model-loading calls.
+template<typename PluginT>
+void register_plugin(PluginT& plugin) {
+    static_assert(
+        std::is_base_of<Plugin, PluginT>::value,
+        "register_plugin expects a class derived from metatomic::Plugin"
+    );
+
+    details::PluginRegistration<PluginT>::plugin = &plugin;
+    details::PluginRegistration<PluginT>::name = plugin.name();
+
+    auto c_plugin = mta_plugin_t{
+        details::PluginRegistration<PluginT>::name.c_str(),
+        &details::PluginRegistration<PluginT>::load_model,
+    };
+
+    mta_register_plugin(c_plugin);
+}
+
+/// Load a plugin dynamic library from the given path.
+inline void load_plugin(const std::string& path) {
+    details::check_status(mta_load_plugin(path.c_str()));
+}
+
+/// Get a handle to a plugin in metatomic's global plugin registry.
+inline PluginHandle plugin(const std::string& name) {
+    return PluginHandle(name);
+}
+
+/// Load a model using the given plugin.
+inline AtomisticModel load_model(
+    const std::string& plugin_name,
+    const std::string& load_from,
+    const std::string& options_json = "{}"
+) {
+    return plugin(plugin_name).load_model(load_from, options_json);
+}
+
+/// Load a model, letting metatomic pick the plugin.
+inline AtomisticModel load_model(
+    const std::string& load_from,
+    const std::string& options_json = "{}"
+) {
+    auto model = mta_model_t{};
+    details::check_status(mta_load_model(
+        nullptr,
+        load_from.c_str(),
+        options_json.c_str(),
+        &model
+    ));
+
+    return AtomisticModel(model);
+}
+
 } // namespace metatomic
diff --git a/metatomic-core/include/metatomic/system.hpp b/metatomic-core/include/metatomic/system.hpp
index 1cae91bdf..15d136b83 100644
--- a/metatomic-core/include/metatomic/system.hpp
+++ b/metatomic-core/include/metatomic/system.hpp
@@ -1,7 +1,267 @@
 #pragma once
 
+#include <string>
+#include <vector>
+
 #include <metatomic.h>
+#include <metatensor.hpp>
+#include <nlohmann/json.hpp>
+
+#include "./metadata.hpp"
+#include "./utils.hpp"
 
 namespace metatomic {
 
+/// A System contains all the information about an atomistic system, and should
+/// be used as input of atomistic models.
+class System final {
+public:
+    /// Create a new `System` from DLPack tensors.
+    ///
+    /// Ownership of all DLPack tensors is transferred to the C API.
+    System(
+        const std::string& length_unit,
+        DLManagedTensorVersioned* types,
+        DLManagedTensorVersioned* positions,
+        DLManagedTensorVersioned* cell,
+        DLManagedTensorVersioned* pbc
+    ): system_(nullptr), is_view_(false) {
+        details::check_status(mta_system_create(
+            length_unit.c_str(),
+            types,
+            positions,
+            cell,
+            pbc,
+            &system_
+        ));
+        details::check_pointer(system_);
+    }
+
+    ~System() {
+        if (system_ != nullptr && !is_view_) {
+            (void)mta_system_free(system_);
+        }
+    }
+
+    System(const System&) = delete;
+    System& operator=(const System&) = delete;
+
+    System(System&& other) noexcept: system_(nullptr), is_view_(true) {
+        *this = std::move(other);
+    }
+
+    System& operator=(System&& other) noexcept {
+        if (system_ != nullptr && !is_view_) {
+            (void)mta_system_free(system_);
+        }
+
+        system_ = other.system_;
+        is_view_ = other.is_view_;
+        other.system_ = nullptr;
+        other.is_view_ = true;
+        return *this;
+    }
+
+    /// Get the number of particles in this system.
+    size_t size() const {
+        uintptr_t result = 0;
+        details::check_status(mta_system_size(system_, &result));
+        return static_cast<size_t>(result);
+    }
+
+    /// Get the length unit used by positions and cell.
+    std::string length_unit() const {
+        mta_string_t length_unit = nullptr;
+        details::check_status(mta_system_get_length_unit(system_, &length_unit));
+        return String(length_unit).str();
+    }
+
+    /// Get particle types for all particles in the system.
+    DLPackTensor types() const {
+        return this->data(MTA_SYSTEM_DATA_TYPES);
+    }
+
+    /// Get the positions for all particles in the system.
+    DLPackTensor positions() const {
+        return this->data(MTA_SYSTEM_DATA_POSITIONS);
+    }
+
+    /// Get the unit cell/bounding box of the system.
+    DLPackTensor cell() const {
+        return this->data(MTA_SYSTEM_DATA_CELL);
+    }
+
+    /// Get the periodic boundary conditions for the system.
+    DLPackTensor pbc() const {
+        return this->data(MTA_SYSTEM_DATA_PBC);
+    }
+
+    /// Add a new pair list in this system.
+    ///
+    /// Ownership of `pairs` is transferred to the C API.
+    void add_pairs(const PairListOptions& options, mts_block_t* pairs) {
+        this->add_pairs(options.to_json(), pairs);
+    }
+
+    /// Add a new pair list in this system.
+    ///
+    /// Ownership of `pairs` is transferred to the C API.
+    void add_pairs(const std::string& options_json, mts_block_t* pairs) {
+        details::check_status(mta_system_add_pairs(system_, options_json.c_str(), pairs));
+    }
+
+    /// Add a new pair list in this system.
+    ///
+    /// Ownership of `pairs` is transferred to the C API.
+    void set_pairs(const PairListOptions& options, mts_block_t* pairs) {
+        this->add_pairs(options, pairs);
+    }
+
+    /// Add a new pair list in this system.
+    ///
+    /// Ownership of `pairs` is transferred to the C API.
+    void set_pairs(const std::string& options_json, mts_block_t* pairs) {
+        this->add_pairs(options_json, pairs);
+    }
+
+    /// Retrieve a previously stored pair list with the given options.
+    const mts_block_t* pairs_raw(const PairListOptions& options) const {
+        return this->pairs_raw(options.to_json());
+    }
+
+    /// Retrieve a previously stored pair list with the given options.
+    const mts_block_t* pairs_raw(const std::string& options_json) const {
+        const mts_block_t* pairs = nullptr;
+        details::check_status(mta_system_get_pairs(system_, options_json.c_str(), &pairs));
+        details::check_pointer(pairs);
+        return pairs;
+    }
+
+    /// Retrieve a previously stored pair list with the given options as a
+    /// non-owning metatensor view.
+    metatensor::TensorBlock pairs(const PairListOptions& options) const {
+        return this->pairs(options.to_json());
+    }
+
+    /// Retrieve a previously stored pair list with the given options as a
+    /// non-owning metatensor view.
+    metatensor::TensorBlock pairs(const std::string& options_json) const {
+        return metatensor::TensorBlock::unsafe_view_from_ptr(
+            const_cast<mts_block_t*>(this->pairs_raw(options_json))
+        );
+    }
+
+    /// Get the options for all pair lists registered with this `System`,
+    /// serialized as a JSON array.
+    std::string known_pairs_json() const {
+        mta_string_t pairs_options = nullptr;
+        details::check_status(mta_system_known_pairs(system_, &pairs_options));
+        return String(pairs_options).str();
+    }
+
+    /// Get the options for all pair lists registered with this `System`.
+    std::vector<PairListOptions> known_pairs() const {
+        auto result = std::vector<PairListOptions>();
+        for (const auto& options: nlohmann::json::parse(this->known_pairs_json())) {
+            result.push_back(PairListOptions::from_json(options.dump()));
+        }
+        return result;
+    }
+
+    /// Get the options for all pair lists registered with this `System`,
+    /// each one serialized as JSON.
+    std::vector<std::string> pairs_options() const {
+        auto result = std::vector<std::string>();
+        for (const auto& options: nlohmann::json::parse(this->known_pairs_json())) {
+            result.push_back(options.dump());
+        }
+        return result;
+    }
+
+    /// Add custom data to this system.
+    ///
+    /// Ownership of `data` is transferred to the C API.
+    void add_data(const std::string& name, mts_tensormap_t* data) {
+        details::check_status(mta_system_add_custom_data(system_, name.c_str(), data));
+    }
+
+    /// Add custom data to this system.
+    ///
+    /// Ownership of `data` is transferred to the C API.
+    void set_data(const std::string& name, mts_tensormap_t* data) {
+        this->add_data(name, data);
+    }
+
+    /// Retrieve custom data stored in this system.
+    ///
+    /// The returned pointer is borrowed and owned by this `System`.
+    const mts_tensormap_t* data_raw(const std::string& name) const {
+        const mts_tensormap_t* data = nullptr;
+        details::check_status(mta_system_get_custom_data(system_, name.c_str(), &data));
+        details::check_pointer(data);
+        return data;
+    }
+
+    /// Get the names of all custom data registered with this `System`.
+    std::string known_data_json() const {
+        mta_string_t names = nullptr;
+        details::check_status(mta_system_known_custom_data(system_, &names));
+        return String(names).str();
+    }
+
+    /// Get the names of all custom data registered with this `System`.
+    std::vector<std::string> known_data() const {
+        return nlohmann::json::parse(this->known_data_json()).get<std::vector<std::string>>();
+    }
+
+    /// Get the names of all custom data registered with this `System`.
+    std::vector<std::string> data_names() const {
+        return this->known_data();
+    }
+
+    /// Get the underlying `mta_system_t` pointer.
+    mta_system_t* as_mta_system_t() & {
+        details::check_pointer(system_);
+        return system_;
+    }
+
+    /// Get the underlying `mta_system_t` pointer.
+    const mta_system_t* as_mta_system_t() const & {
+        details::check_pointer(system_);
+        return system_;
+    }
+
+    mta_system_t* as_mta_system_t() && = delete;
+
+    /// Take ownership of a raw `mta_system_t*`.
+    static System unsafe_from_ptr(mta_system_t* system) {
+        return System(system, false);
+    }
+
+    /// Create a non-owning view of a raw `mta_system_t*`.
+    static System unsafe_view_from_ptr(const mta_system_t* system) {
+        return System(const_cast<mta_system_t*>(system), true);
+    }
+
+    /// Release the raw `mta_system_t*` without freeing it.
+    mta_system_t* release() {
+        auto* system = system_;
+        system_ = nullptr;
+        return system;
+    }
+
+private:
+    explicit System(mta_system_t* system, bool is_view): system_(system), is_view_(is_view) {}
+
+    DLPackTensor data(mta_system_data_kind request) const {
+        DLManagedTensorVersioned* data = nullptr;
+        details::check_status(mta_system_get_data(system_, request, &data));
+        details::check_pointer(data);
+        return DLPackTensor(data);
+    }
+
+    mta_system_t* system_;
+    bool is_view_;
+};
+
 } // namespace metatomic
diff --git a/metatomic-core/include/metatomic/utils.hpp b/metatomic-core/include/metatomic/utils.hpp
index 1cae91bdf..056e18224 100644
--- a/metatomic-core/include/metatomic/utils.hpp
+++ b/metatomic-core/include/metatomic/utils.hpp
@@ -1,7 +1,251 @@
 #pragma once
 
+#include <cstring>
+#include <cstdio>
+
+#include <exception>
+#include <stdexcept>
+#include <string>
+#include <utility>
+
 #include <metatomic.h>
 
 namespace metatomic {
 
+/// Exception class used for all errors in metatomic.
+class Error: public std::runtime_error {
+public:
+    /// Create a new Error with the given `message`.
+    explicit Error(const std::string& message): std::runtime_error(message) {}
+};
+
+/// RAII wrapper around a `DLManagedTensorVersioned*`.
+///
+/// This owns the DLPack managed tensor object, and calls its deleter when the
+/// wrapper is destroyed.
+class DLPackTensor final {
+public:
+    /// Create an empty wrapper.
+    DLPackTensor(): tensor_(nullptr) {}
+
+    /// Take ownership of an existing DLPack managed tensor.
+    explicit DLPackTensor(DLManagedTensorVersioned* tensor): tensor_(tensor) {}
+
+    ~DLPackTensor() {
+        if (tensor_ != nullptr && tensor_->deleter != nullptr) {
+            tensor_->deleter(tensor_);
+        }
+    }
+
+    DLPackTensor(const DLPackTensor&) = delete;
+    DLPackTensor& operator=(const DLPackTensor&) = delete;
+
+    DLPackTensor(DLPackTensor&& other) noexcept: DLPackTensor() {
+        *this = std::move(other);
+    }
+
+    DLPackTensor& operator=(DLPackTensor&& other) noexcept {
+        if (tensor_ != nullptr && tensor_->deleter != nullptr) {
+            tensor_->deleter(tensor_);
+        }
+
+        tensor_ = other.tensor_;
+        other.tensor_ = nullptr;
+        return *this;
+    }
+
+    /// Check if this wrapper contains a tensor.
+    explicit operator bool() const {
+        return tensor_ != nullptr;
+    }
+
+    /// Get the underlying DLPack managed tensor.
+    DLManagedTensorVersioned* get() const {
+        return tensor_;
+    }
+
+    /// Get the underlying DLPack managed tensor.
+    DLManagedTensorVersioned* as_dlpack() const {
+        return tensor_;
+    }
+
+    /// Release the DLPack managed tensor without calling its deleter.
+    DLManagedTensorVersioned* release() {
+        auto* tensor = tensor_;
+        tensor_ = nullptr;
+        return tensor;
+    }
+
+private:
+    DLManagedTensorVersioned* tensor_;
+};
+
+namespace details {
+    /// Check if a return status from the C API indicates an error, and throw a
+    /// `metatomic::Error` with the last error message if this is the case.
+    inline void check_status(mta_status_t status) {
+        if (status == MTA_SUCCESS) {
+            return;
+        }
+
+        const char* message = nullptr;
+        const char* origin = nullptr;
+        void* data = nullptr;
+        (void)mta_last_error(&message, &origin, &data);
+
+        if (origin != nullptr && std::strcmp(origin, "C++ exception") == 0 && data != nullptr) {
+            std::rethrow_exception(*static_cast<std::exception_ptr*>(data));
+        }
+
+        throw Error(message == nullptr ? "unknown error" : message);
+    }
+
+    /// Call the given `function`, catching any C++ exception and translating it
+    /// to a metatomic error code.
+    ///
+    /// This is required to prevent callbacks unwinding through the C API.
+    template<typename Function, typename ...Args>
+    inline mta_status_t catch_exceptions(Function function, Args ...args) {
+        try {
+            function(std::move(args)...);
+            return MTA_SUCCESS;
+        } catch (...) {
+            auto* exception_ptr = new std::exception_ptr(std::current_exception());
+
+            const char* message = nullptr;
+            try {
+                std::rethrow_exception(*exception_ptr);
+            } catch (const std::exception& e) {
+                message = e.what();
+            } catch (...) {
+                message = "C++ code threw an exception that was not a std::exception";
+            }
+
+            auto status = mta_set_last_error(
+                message,
+                "C++ exception",
+                exception_ptr,
+                [](void* ptr) { delete static_cast<std::exception_ptr*>(ptr); }
+            );
+
+            if (status != MTA_SUCCESS) {
+                std::fprintf(
+                    stderr,
+                    "INTERNAL ERROR: unable to set last error after C++ callback failure (status: %d). ",
+                    static_cast<int>(status)
+                );
+                if (message != nullptr) {
+                    std::fprintf(stderr, "C++ error was: %s\n", message);
+                } else {
+                    std::fprintf(stderr, "Unknown C++ error\n");
+                }
+                delete exception_ptr;
+            }
+
+            return MTA_ERROR_OTHER;
+        }
+    }
+
+    /// Check if a pointer allocated by the C API is null.
+    inline void check_pointer(const void* pointer) {
+        if (pointer != nullptr) {
+            return;
+        }
+
+        const char* message = nullptr;
+        const char* origin = nullptr;
+        void* data = nullptr;
+        (void)mta_last_error(&message, &origin, &data);
+
+        if (origin != nullptr && std::strcmp(origin, "C++ exception") == 0 && data != nullptr) {
+            std::rethrow_exception(*static_cast<std::exception_ptr*>(data));
+        }
+
+        throw Error(message == nullptr ? "received a null pointer from the metatomic C API" : message);
+    }
+} // namespace details
+
+/// RAII wrapper for `mta_string_t`.
+class String final {
+public:
+    /// Create an empty string wrapper.
+    String(): string_(nullptr) {}
+
+    /// Create a new string managed by metatomic.
+    explicit String(const std::string& string): string_(mta_string_create(string.c_str())) {
+        details::check_pointer(string_);
+    }
+
+    /// Take ownership of an existing `mta_string_t`.
+    explicit String(mta_string_t string): string_(string) {}
+
+    ~String() {
+        if (string_ != nullptr) {
+            mta_string_free(string_);
+        }
+    }
+
+    String(const String&) = delete;
+    String& operator=(const String&) = delete;
+
+    String(String&& other) noexcept: String() {
+        *this = std::move(other);
+    }
+
+    String& operator=(String&& other) noexcept {
+        if (string_ != nullptr) {
+            mta_string_free(string_);
+        }
+
+        string_ = other.string_;
+        other.string_ = nullptr;
+        return *this;
+    }
+
+    /// Get the underlying C pointer.
+    mta_string_t as_mta_string_t() const {
+        return string_;
+    }
+
+    /// View this string as a null-terminated C string.
+    const char* c_str() const {
+        if (string_ == nullptr) {
+            return "";
+        }
+
+        auto* ptr = mta_string_view(string_);
+        details::check_pointer(ptr);
+        return ptr;
+    }
+
+    /// Copy this string into a C++ string.
+    std::string str() const {
+        return std::string(this->c_str());
+    }
+
+private:
+    mta_string_t string_;
+};
+
+/// Get the runtime version of metatomic as a string.
+inline std::string version() {
+    auto* raw = mta_version();
+    details::check_pointer(raw);
+    return std::string(raw);
+}
+
+/// Get the conversion factor from `from_unit` to `to_unit`.
+inline double unit_conversion_factor(const std::string& from_unit, const std::string& to_unit) {
+    double conversion = 0.0;
+    details::check_status(mta_unit_conversion_factor(from_unit.c_str(), to_unit.c_str(), &conversion));
+    return conversion;
+}
+
+/// Format model metadata JSON for display.
+inline std::string format_metadata(const std::string& metadata) {
+    mta_string_t printed = nullptr;
+    details::check_status(mta_format_metadata(metadata.c_str(), &printed));
+    return String(printed).str();
+}
+
 } // namespace metatomic
diff --git a/metatomic-core/tests/utils/mod.rs b/metatomic-core/tests/utils/mod.rs
index e12e6897c..78540ba96 100644
--- a/metatomic-core/tests/utils/mod.rs
+++ b/metatomic-core/tests/utils/mod.rs
@@ -148,6 +148,7 @@ pub fn create_python_venv(build_dir: PathBuf) -> PathBuf {
         let mut cmd = Command::new(find_python());
         cmd.arg("-m");
         cmd.arg("venv");
+        cmd.arg("--clear");
         cmd.arg(&build_dir);
 
         run_command(cmd, "python to create virtualenv with `venv`");