MAINT: Migration of the codebase to C++

.gitignore

@@ -7,6 +7,7 @@
 *.bak
 .project
 .pydevproject
+.vscode
 
 # Compiled source #
 ###################
@@ -48,4 +49,10 @@ build_standalone/*.f90
 build_standalone/*.c
 build_standalone/*.h
 build_standalone/*.inc
-build_standalone/*.classes
+build_standalone/*.classes
+
+# Meson subprojects
+cpp/subprojects/packagecache
+cpp/subprojects/Catch2-*
+cpp/subprojects/eigen-*
+cpp/subprojects/pybind11-*

discover_version.py

@@ -97,9 +97,19 @@ def get_version_from_git():
     # Make version PEP 440 compliant
     if dirty:
         version = version.replace('-dirty', '')
-    version = version.replace('-', '+', 1)
-    if dirty:
-        version += '-dirty'
+
+    # Check if there's a dash (indicating commits after tag, e.g., 1.2.7-5-gabcdef)
+    if '-' in version:
+        # Replace first '-' with '+' for local version separator
+        version = version.replace('-', '+', 1)
+        # Replace remaining '-' with '.' for PEP 440 compliance
+        version = version.replace('-', '.')
+        if dirty:
+            version += '.dirty'
+    else:
+        # No commits after tag (e.g., just 1.2.7), need to add local version segment
+        if dirty:
+            version += '+dirty'
 
     return dirty, version, git_hash
 

examples/STANDALONE/C60_with_Rebo2/md.dat

@@ -36,4 +36,12 @@ Simulation {
 
   Rebo2 { };
 
+  #
+  # Trajectory output (every ps) to .NC file
+  #
+
+  OutputNC {
+    freq = "100.0";
+  };
+
 };

lib/include/atomistica/atomistica.hpp

@@ -0,0 +1,57 @@
+// ======================================================================
+// Atomistica - Interatomic potential library and molecular dynamics code
+// https://github.com/Atomistica/atomistica
+//
+// Copyright (2005-2024) Lars Pastewka <lars.pastewka@imtek.uni-freiburg.de>
+// and others. See the AUTHORS file in the top-level Atomistica directory.
+//
+// This program is free software: you can redistribute it and/or modify
+// it under the terms of the GNU General Public License as published by
+// the Free Software Foundation, either version 2 of the License, or
+// (at your option) any later version.
+//
+// This program 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 General Public License for more details.
+//
+// You should have received a copy of the GNU General Public License
+// along with this program.  If not, see <http://www.gnu.org/licenses/>.
+// ======================================================================
+
+#pragma once
+
+// Configuration and common types
+#include "config.hpp"
+
+// Core components
+#include "core/atomic_system.hpp"
+#include "core/neighbor_list.hpp"
+
+// Math utilities
+#include "math/cutoff_functions.hpp"
+#include "math/spline.hpp"
+
+// Potential base
+#include "potentials/potential_base.hpp"
+
+// Pair potentials
+#include "potentials/pair/lj.hpp"
+
+// Bond-order potentials
+#include "potentials/bop/bop_base.hpp"
+#include "potentials/bop/tersoff.hpp"
+#include "potentials/bop/brenner.hpp"
+#include "potentials/bop/kumagai.hpp"
+#include "potentials/bop/rebo2.hpp"
+
+// Coulomb potentials
+#include "potentials/coulomb/coulomb.hpp"
+#include "potentials/coulomb/pme.hpp"
+#include "potentials/coulomb/fmm.hpp"
+
+// Tight-binding
+#include "tightbinding/tightbinding.hpp"
+
+// Integrators
+#include "integrators/integrators.hpp"

lib/include/atomistica/config.hpp

@@ -0,0 +1,47 @@
+// ======================================================================
+// Atomistica - Interatomic potential library and molecular dynamics code
+// https://github.com/Atomistica/atomistica
+//
+// Copyright (2005-2024) Lars Pastewka <lars.pastewka@imtek.uni-freiburg.de>
+// and others. See the AUTHORS file in the top-level Atomistica directory.
+//
+// This program is free software: you can redistribute it and/or modify
+// it under the terms of the GNU General Public License as published by
+// the Free Software Foundation, either version 2 of the License, or
+// (at your option) any later version.
+//
+// This program 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 General Public License for more details.
+//
+// You should have received a copy of the GNU General Public License
+// along with this program.  If not, see <http://www.gnu.org/licenses/>.
+// ======================================================================
+
+#pragma once
+
+#include <Eigen/Core>
+#include <Eigen/Geometry>
+
+namespace atomistica {
+
+// Scalar type (can be changed for single precision if needed)
+using Scalar = double;
+
+// Vector and matrix types
+using Vec3 = Eigen::Matrix<Scalar, 3, 1>;
+using Mat3 = Eigen::Matrix<Scalar, 3, 3>;
+using VecX = Eigen::Matrix<Scalar, Eigen::Dynamic, 1>;
+using MatX = Eigen::Matrix<Scalar, Eigen::Dynamic, Eigen::Dynamic>;
+using MatX3 = Eigen::Matrix<Scalar, Eigen::Dynamic, 3>;  // N x 3 matrix for forces
+
+// Array types for per-atom data (row-major for cache efficiency when iterating atoms)
+using Array3X = Eigen::Array<Scalar, 3, Eigen::Dynamic>;
+using ArrayX = Eigen::Array<Scalar, Eigen::Dynamic, 1>;
+using ArrayXi = Eigen::Array<int, Eigen::Dynamic, 1>;
+
+// Constants
+constexpr Scalar PI = 3.14159265358979323846;
+
+} // namespace atomistica

lib/include/atomistica/core/atomic_system.hpp

@@ -0,0 +1,220 @@
+// ======================================================================
+// Atomistica - Interatomic potential library and molecular dynamics code
+// https://github.com/Atomistica/atomistica
+//
+// Copyright (2005-2024) Lars Pastewka <lars.pastewka@imtek.uni-freiburg.de>
+// and others. See the AUTHORS file in the top-level Atomistica directory.
+//
+// This program is free software: you can redistribute it and/or modify
+// it under the terms of the GNU General Public License as published by
+// the Free Software Foundation, either version 2 of the License, or
+// (at your option) any later version.
+//
+// This program 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 General Public License for more details.
+//
+// You should have received a copy of the GNU General Public License
+// along with this program.  If not, see <http://www.gnu.org/licenses/>.
+// ======================================================================
+
+#pragma once
+
+#include <array>
+#include <cstddef>
+#include <memory>
+#include <string>
+#include <unordered_map>
+#include <variant>
+#include <vector>
+
+#include "../config.hpp"
+
+namespace atomistica {
+
+/**
+ * @brief Type-erased property storage for per-atom data
+ *
+ * Allows storing arbitrary per-atom properties (charges, velocities, etc.)
+ * with type-safe access.
+ */
+class PropertyMap {
+public:
+    using PropertyVariant = std::variant<
+        ArrayX,      // Scalar per-atom (e.g., charge, energy)
+        Array3X,     // 3-vector per-atom (e.g., velocity, force)
+        ArrayXi      // Integer per-atom (e.g., type, molecule ID)
+    >;
+
+    template<typename T>
+    void add(const std::string& name, std::size_t size);
+
+    template<typename T>
+    T& get(const std::string& name);
+
+    template<typename T>
+    const T& get(const std::string& name) const;
+
+    bool has(const std::string& name) const;
+    void remove(const std::string& name);
+    void resize(std::size_t n);
+
+private:
+    std::unordered_map<std::string, PropertyVariant> properties_;
+};
+
+/**
+ * @brief Container for atomic system data
+ *
+ * Stores positions, atomic numbers, cell vectors, and periodic boundary conditions.
+ * Also provides extensible per-atom property storage.
+ *
+ * This replaces the Fortran particles_t type.
+ */
+class AtomicSystem {
+public:
+    AtomicSystem() = default;
+    explicit AtomicSystem(std::size_t num_atoms);
+
+    // Resize the system
+    void resize(std::size_t num_atoms);
+
+    // Number of atoms
+    std::size_t num_atoms() const { return num_atoms_; }
+
+    // Cell vectors (columns are the lattice vectors a, b, c)
+    Mat3& cell() { return cell_; }
+    const Mat3& cell() const { return cell_; }
+    void set_cell(const Mat3& cell);
+
+    // Periodic boundary conditions
+    std::array<bool, 3>& pbc() { return pbc_; }
+    const std::array<bool, 3>& pbc() const { return pbc_; }
+
+    // Positions (3 x N array, column i is position of atom i)
+    Array3X& positions() { return positions_; }
+    const Array3X& positions() const { return positions_; }
+
+    // Single atom position access - always returns a Vec3 copy for type safety
+    // Use positions().col(i) for assignment or set_position() method
+    Vec3 position(std::size_t i) const { return positions_.col(i).matrix(); }
+
+    // Set single atom position
+    void set_position(std::size_t i, const Vec3& r) {
+        positions_.col(i) = r.array();
+        ++position_revision_;
+    }
+
+    // Atomic numbers
+    ArrayXi& atomic_numbers() { return atomic_numbers_; }
+    const ArrayXi& atomic_numbers() const { return atomic_numbers_; }
+
+    // Single atom atomic number access
+    int atomic_number(std::size_t i) const { return atomic_numbers_[i]; }
+
+    // Masses
+    ArrayX& masses() { return masses_; }
+    const ArrayX& masses() const { return masses_; }
+
+    // Single atom mass access
+    Scalar mass(std::size_t i) const { return masses_[i]; }
+    void set_mass(std::size_t i, Scalar m) { masses_[i] = m; }
+
+    // Velocities
+    Array3X& velocities() { return velocities_; }
+    const Array3X& velocities() const { return velocities_; }
+
+    // Single atom velocity access - always returns a Vec3 copy for type safety
+    // Use velocities().col(i) for assignment or set_velocity() method
+    Vec3 velocity(std::size_t i) const { return velocities_.col(i).matrix(); }
+    void set_velocity(std::size_t i, const Vec3& v) { velocities_.col(i) = v.array(); }
+
+    // Add atom with position and mass
+    void add_atom(int Z, const Vec3& r, Scalar m = 1.0) {
+        std::size_t i = num_atoms_;
+        resize(num_atoms_ + 1);
+        positions_.col(i) = r.array();
+        atomic_numbers_[i] = Z;
+        masses_[i] = m;
+        velocities_.col(i).setZero();
+    }
+
+    // Forces (accumulated by potentials)
+    Array3X& forces() { return forces_; }
+    const Array3X& forces() const { return forces_; }
+
+    // Zero forces
+    void zero_forces();
+
+    // Extensible property storage
+    PropertyMap& properties() { return properties_; }
+    const PropertyMap& properties() const { return properties_; }
+
+    // Cell operations
+    Scalar volume() const;
+    Mat3 inverse_cell() const;
+
+    // Minimum image convention for distance vectors
+    Vec3 minimum_image(const Vec3& dr) const;
+
+    // Wrap position into cell
+    Vec3 wrap_position(const Vec3& r) const;
+
+    // Change tracking for lazy updates (e.g., neighbor lists)
+    int position_revision() const { return position_revision_; }
+    int cell_revision() const { return cell_revision_; }
+    void positions_changed() { ++position_revision_; }
+    void cell_changed() { ++cell_revision_; }
+
+private:
+    std::size_t num_atoms_ = 0;
+
+    Mat3 cell_ = Mat3::Identity();
+    std::array<bool, 3> pbc_ = {true, true, true};
+
+    Array3X positions_;
+    ArrayXi atomic_numbers_;
+    ArrayX masses_;
+    Array3X velocities_;
+    Array3X forces_;
+
+    PropertyMap properties_;
+
+    int position_revision_ = 0;
+    int cell_revision_ = 0;
+
+    // Cached inverse cell (lazily computed)
+    mutable Mat3 inverse_cell_;
+    mutable int inverse_cell_revision_ = -1;
+};
+
+// Template implementations
+
+template<typename T>
+void PropertyMap::add(const std::string& name, std::size_t size) {
+    if constexpr (std::is_same_v<T, ArrayX>) {
+        ArrayX arr = ArrayX::Zero(size);
+        properties_[name] = std::move(arr);
+    } else if constexpr (std::is_same_v<T, Array3X>) {
+        Array3X arr = Array3X::Zero(3, size);
+        properties_[name] = std::move(arr);
+    } else if constexpr (std::is_same_v<T, ArrayXi>) {
+        ArrayXi arr = ArrayXi::Zero(size);
+        properties_[name] = std::move(arr);
+    } else {
+        static_assert(sizeof(T) == 0, "Unsupported property type");
+    }
+}
+
+template<typename T>
+T& PropertyMap::get(const std::string& name) {
+    return std::get<T>(properties_.at(name));
+}
+
+template<typename T>
+const T& PropertyMap::get(const std::string& name) const {
+    return std::get<T>(properties_.at(name));
+}
+
+} // namespace atomistica