Introduce site and hopping registries

They track site and hopping families introduced both by the
lattice and generators.

cppcore/CMakeLists.txt

@@ -61,6 +61,7 @@ add_library(cppcore
     include/system/CompressedSublattices.hpp
     include/system/Foundation.hpp
     include/system/HoppingBlocks.hpp
+    include/system/Registry.hpp
     include/system/Shape.hpp
     include/system/StructureModifiers.hpp
     include/system/Symmetry.hpp
@@ -88,6 +89,7 @@ add_library(cppcore
     src/system/CompressedSublattices.cpp
     src/system/Foundation.cpp
     src/system/HoppingBlocks.cpp
+    src/system/Registry.cpp
     src/system/Shape.cpp
     src/system/StructureModifiers.cpp
     src/system/Symmetry.cpp

cppcore/include/Lattice.hpp

@@ -1,24 +1,10 @@
 #pragma once
-#include "numeric/dense.hpp"
-#include "detail/opaque_alias.hpp"
+#include "system/Registry.hpp"
 
 #include <vector>
-#include <string>
-#include <unordered_map>
 
 namespace cpb {
 
-// TODO: replace with proper string_view
-using string_view = std::string const&;
-
-/// Sublattice and hopping ID data types
-using SubID = detail::OpaqueIntegerAlias<class SubIDTag>;
-using SubAliasID = detail::OpaqueIntegerAlias<class SubAliasIDTag>;
-using HopID = detail::OpaqueIntegerAlias<class HopIDTag>;
-
-/// Map from friendly sublattice/hopping name to numeric ID
-using NameMap = std::unordered_map<std::string, storage_idx_t>;
-
 class OptimizedUnitCell;
 
 /**
@@ -118,41 +104,20 @@ class Lattice {
     /// Get a single vector -- Expects: index < lattice.ndim()
     Cartesian vector(size_t index) const { return vectors[index]; }
 
-    /// Access sublattice information by name or ID
-    Sublattice const& sublattice(string_view name) const;
-    Sublattice const& sublattice(SubID id) const;
-    Sublattice const& operator[](string_view name) const { return sublattice(name); }
-    Sublattice const& operator[](SubID id) const { return sublattice(id); }
-    /// Assess hopping family information by name or ID
-    HoppingFamily const& hopping_family(string_view name) const;
-    HoppingFamily const& hopping_family(HopID id) const;
-    HoppingFamily const& operator()(string_view name) const { return hopping_family(name); }
-    HoppingFamily const& operator()(HopID id) const { return hopping_family(id); }
-
-    /// Return name for this ID
-    string_view sublattice_name(SubID id) const;
-    /// Return name for this ID
-    string_view hopping_family_name(HopID id) const;
-
     /// Get the maximum possible number of hoppings from any site of this lattice
     int max_hoppings() const;
     /// Does at least one sublattice have non-zero onsite energy on the main diagonal?
     bool has_diagonal_terms() const;
     /// Does at least one sublattice have non-zero onsite energy (including off-diagonal terms)?
     bool has_onsite_energy() const;
-    /// Does at least one sublattice have multiple orbitals?
-    bool has_multiple_orbitals() const;
-    /// Is at least one hopping complex?
-    bool has_complex_hoppings() const;
 
 public: // optimized access
     /// Return the optimized version of the lattice structure -- see `OptimizedLatticeStructure`
     OptimizedUnitCell optimized_unit_cell() const;
-
-    /// Mapping from friendly sublattice names to their unique IDs
-    NameMap sub_name_map() const;
-    /// Mapping from friendly hopping names to their unique IDs
-    NameMap hop_name_map() const;
+    /// Site information for the system builder
+    SiteRegistry site_registry() const;
+    /// Hopping information for the system builder
+    HoppingRegistry hopping_registry() const;
 
 public: // utilities
     /// Calculate the spatial position of a unit cell or a sublattice site if specified
@@ -167,6 +132,14 @@ class Lattice {
     Vector3f translate_coordinates(Cartesian position) const;
 
 private:
+    /// Access sublattice information by name or ID
+    Sublattice const& sublattice(string_view name) const;
+    Sublattice const& sublattice(SubID id) const;
+
+    /// Assess hopping family information by name or ID
+    HoppingFamily const& hopping_family(string_view name) const;
+    HoppingFamily const& hopping_family(HopID id) const;
+
     SubID make_unique_sublattice_id(string_view name);
 
 private:

cppcore/include/Model.hpp

@@ -1,4 +1,5 @@
 #pragma once
+#include "system/Registry.hpp"
 #include "system/System.hpp"
 #include "system/Shape.hpp"
 #include "system/Symmetry.hpp"
@@ -17,11 +18,11 @@ namespace cpb {
 
 class Model {
 public:
-    Model(Lattice const& lattice) : lattice(lattice) {}
+    Model(Lattice const& lattice);
 
     /// The arguments can any type accepted by `Model::add()`
     template<class... Args>
-    Model(Lattice const& lattice, Args&&... args) : lattice(lattice) {
+    Model(Lattice const& lattice, Args&&... args) : Model(lattice) {
         detail::eval_ordered({(add(std::forward<Args>(args)), 0)...});
     }
 
@@ -42,13 +43,17 @@ class Model {
     void set_wave_vector(Cartesian const& k);
 
 public:
+    /// Are any of the onsite or hopping energies given as matrices instead of scalars?
+    bool is_multiorbital() const;
     /// Uses double precision values in the Hamiltonian matrix?
     bool is_double() const;
     /// Uses complex values in the Hamiltonian matrix?
     bool is_complex() const;
 
 public: // get parameters
     Lattice const& get_lattice() const { return lattice; }
+    SiteRegistry const& get_site_registry() const { return site_registry; }
+    HoppingRegistry const& get_hopping_registry() const { return hopping_registry; }
     Primitive const& get_primitive() const { return primitive; }
     Shape const& get_shape() const { return shape; }
     TranslationalSymmetry const& get_symmetry() const { return symmetry; }
@@ -87,6 +92,9 @@ class Model {
 
 private:
     Lattice lattice;
+    SiteRegistry site_registry;
+    HoppingRegistry hopping_registry;
+
     Primitive primitive;
     Shape shape;
     TranslationalSymmetry symmetry;

cppcore/include/hamiltonian/Hamiltonian.hpp

@@ -72,16 +72,16 @@ class Hamiltonian {
 namespace detail {
 
 template<class scalar_t>
-void build_main(SparseMatrixX<scalar_t>& matrix, System const& system,
+void build_main(SparseMatrixX<scalar_t>& matrix, System const& system, Lattice const& lattice,
                 HamiltonianModifiers const& modifiers, bool simple_build) {
     auto const size = system.hamiltonian_size();
     matrix.resize(size, size);
 
     if (simple_build) {
         // Fast path: No generators were used (only unit cell replication + modifiers)
         // so we can easily predict the maximum number of non-zero values per row.
-        auto const has_diagonal = system.lattice.has_diagonal_terms() || !modifiers.onsite.empty();
-        auto const num_per_row = system.lattice.max_hoppings() + has_diagonal;
+        auto const has_diagonal = lattice.has_diagonal_terms() || !modifiers.onsite.empty();
+        auto const num_per_row = lattice.max_hoppings() + has_diagonal;
         matrix.reserve(ArrayXi::Constant(size, num_per_row));
 
         modifiers.apply_to_onsite<scalar_t>(system, [&](idx_t i, idx_t j, scalar_t onsite) {
@@ -161,11 +161,11 @@ inline bool is(Hamiltonian const& h) {
 }
 
 template<class scalar_t>
-Hamiltonian make(System const& system, HamiltonianModifiers const& modifiers,
-                 Cartesian k_vector, bool simple_build) {
+Hamiltonian make(System const& system, Lattice const& lattice,
+                 HamiltonianModifiers const& modifiers, Cartesian k_vector, bool simple_build) {
     auto matrix = std::make_shared<SparseMatrixX<scalar_t>>();
 
-    detail::build_main(*matrix, system, modifiers, simple_build);
+    detail::build_main(*matrix, system, lattice, modifiers, simple_build);
     detail::build_periodic(*matrix, system, modifiers, k_vector);
 
     matrix->makeCompressed();

cppcore/include/hamiltonian/HamiltonianModifiers.hpp

@@ -105,8 +105,7 @@ namespace detail {
 
 template<class scalar_t, class Fn>
 void HamiltonianModifiers::apply_to_onsite(System const& system, Fn lambda) const {
-    auto const& lattice = system.lattice;
-    auto const has_intrinsic_onsite = lattice.has_onsite_energy();
+    auto const has_intrinsic_onsite = system.site_registry.has_nonzero_energy();
     if (!has_intrinsic_onsite && onsite.empty()) {
         return;
     }
@@ -118,8 +117,9 @@ void HamiltonianModifiers::apply_to_onsite(System const& system, Fn lambda) cons
 
         if (has_intrinsic_onsite) {
             // Intrinsic lattice onsite energy -- just replicate the value at each site
-            auto const& site_family = lattice[SubID(sub.alias_id())];
-            auto const intrinsic_energy = num::force_cast<scalar_t>(site_family.energy);
+            auto const intrinsic_energy = num::force_cast<scalar_t>(
+                system.site_registry.energy(sub.id())
+            );
 
             auto start = idx_t{0};
             for (auto const& value : intrinsic_energy) {
@@ -133,7 +133,7 @@ void HamiltonianModifiers::apply_to_onsite(System const& system, Fn lambda) cons
             auto onsite_ref = (norb == 1) ? arrayref(onsite_energy.data(), nsites)
                                           : arrayref(onsite_energy.data(), norb, norb, nsites);
             auto const position_ref = system.positions.segment(sub.sys_start(), nsites);
-            auto const sub_name = lattice.sublattice_name(SubID(sub.alias_id()));
+            auto const sub_name = system.site_registry.name(sub.id());
 
             for (auto const& modifier : onsite) {
                 modifier.apply(onsite_ref, position_ref, sub_name);
@@ -240,14 +240,15 @@ template<class scalar_t, class SystemOrBoundary, class Fn>
 void HamiltonianModifiers::apply_to_hoppings_impl(System const& system,
                                                   SystemOrBoundary const& system_or_boundary,
                                                   Fn lambda) const {
-    auto const& lattice = system.lattice;
+    auto const& hopping_registry = system.hopping_registry;
 
     // Fast path: Modifiers don't need to be applied and the single-orbital model
     // allows direct mapping between sites and Hamiltonian matrix elements.
-    if (hopping.empty() && !lattice.has_multiple_orbitals()) {
+    if (hopping.empty() && !hopping_registry.has_multiple_orbitals()) {
         for (auto const& block : system_or_boundary.hopping_blocks) {
-            auto const& hopping_family = lattice(block.family_id());
-            auto const energy = num::force_cast<scalar_t>(hopping_family.energy);
+            auto const energy = num::force_cast<scalar_t>(
+                hopping_registry.energy(block.family_id())
+            );
 
             auto const value = energy(0, 0); // single orbital
             for (auto const& coo : block.coordinates()) {
@@ -263,11 +264,11 @@ void HamiltonianModifiers::apply_to_hoppings_impl(System const& system,
     for (auto const& block : system_or_boundary.hopping_blocks) {
         if (block.size() == 0) { continue; }
 
-        auto const& hopping_family = lattice(block.family_id());
-        auto const hopping_name = lattice.hopping_family_name(block.family_id());
-        auto const index_translator = IndexTranslator(system, hopping_family.energy);
+        auto const& hopping_energy = hopping_registry.energy(block.family_id());
+        auto const hopping_name = hopping_registry.name(block.family_id());
+        auto const index_translator = IndexTranslator(system, hopping_energy);
 
-        auto buffer = HoppingBuffer<scalar_t>(hopping_family.energy, block.size());
+        auto buffer = HoppingBuffer<scalar_t>(hopping_energy, block.size());
         for (auto const coo_slice : sliced(block.coordinates(), buffer.size)) {
             auto size = idx_t{0};
             for (auto const& coo : coo_slice) {

cppcore/include/leads/HamiltonianPair.hpp

@@ -5,22 +5,24 @@ namespace cpb { namespace leads {
 
 namespace detail {
     template<class scalar_t>
-    Hamiltonian make_h0(System const& lead_system, HamiltonianModifiers const& modifiers) {
+    Hamiltonian make_h0(System const& lead_system, Lattice const& lattice,
+                        HamiltonianModifiers const& modifiers) {
         auto h0 = std::make_shared<SparseMatrixX<scalar_t>>();
-        cpb::detail::build_main(*h0, lead_system, modifiers, /*simple_build*/true);
+        cpb::detail::build_main(*h0, lead_system, lattice, modifiers, /*simple_build*/true);
         h0->makeCompressed();
         cpb::detail::throw_if_invalid(*h0);
         return h0;
     }
 
     template<class scalar_t>
-    Hamiltonian make_h1(System const& system, HamiltonianModifiers const& modifiers) {
+    Hamiltonian make_h1(System const& system, Lattice const& lattice,
+                        HamiltonianModifiers const& modifiers) {
         auto h1 = std::make_shared<SparseMatrixX<scalar_t>>();
         auto& matrix = *h1;
 
         auto const size = system.hamiltonian_size();
         matrix.resize(size, size);
-        matrix.reserve(ArrayXi::Constant(size, system.lattice.max_hoppings()));
+        matrix.reserve(ArrayXi::Constant(size, lattice.max_hoppings()));
 
         modifiers.apply_to_hoppings<scalar_t>(system, 0, [&](idx_t i, idx_t j, scalar_t hopping) {
             matrix.insert(i, j) = hopping;
@@ -40,23 +42,23 @@ struct HamiltonianPair {
     Hamiltonian h0; ///< hoppings within the unit cell
     Hamiltonian h1; ///< hoppings between unit cells
 
-    HamiltonianPair(System const& lead_system, HamiltonianModifiers const& modifiers,
-                    bool is_double, bool is_complex) {
+    HamiltonianPair(System const& lead_system, Lattice const& lattice,
+                    HamiltonianModifiers const& modifiers, bool is_double, bool is_complex) {
         if (is_double) {
             if (is_complex) {
-                h0 = detail::make_h0<std::complex<double>>(lead_system, modifiers);
-                h1 = detail::make_h1<std::complex<double>>(lead_system, modifiers);
+                h0 = detail::make_h0<std::complex<double>>(lead_system, lattice, modifiers);
+                h1 = detail::make_h1<std::complex<double>>(lead_system, lattice, modifiers);
             } else {
-                h0 = detail::make_h0<double>(lead_system, modifiers);
-                h1 = detail::make_h1<double>(lead_system, modifiers);
+                h0 = detail::make_h0<double>(lead_system, lattice, modifiers);
+                h1 = detail::make_h1<double>(lead_system, lattice, modifiers);
             }
         } else {
             if (is_complex) {
-                h0 = detail::make_h0<std::complex<float>>(lead_system, modifiers);
-                h1 = detail::make_h1<std::complex<float>>(lead_system, modifiers);
+                h0 = detail::make_h0<std::complex<float>>(lead_system, lattice, modifiers);
+                h1 = detail::make_h1<std::complex<float>>(lead_system, lattice, modifiers);
             } else {
-                h0 = detail::make_h0<float>(lead_system, modifiers);
-                h1 = detail::make_h1<float>(lead_system, modifiers);
+                h0 = detail::make_h0<float>(lead_system, lattice, modifiers);
+                h1 = detail::make_h1<float>(lead_system, lattice, modifiers);
             }
         }
     }

cppcore/include/leads/Leads.hpp

@@ -49,7 +49,8 @@ class Leads {
     void make_structure(Foundation const& foundation);
 
     /// Create a Hamiltonian pair for each lead
-    void make_hamiltonian(HamiltonianModifiers const& modifiers, bool is_double, bool is_complex);
+    void make_hamiltonian(Lattice const& lattice, HamiltonianModifiers const& modifiers,
+                          bool is_double, bool is_complex);
 
     /// Clear any existing structural data, implies clearing Hamiltonian
     void clear_structure();

cppcore/include/system/CompressedSublattices.hpp

@@ -15,7 +15,7 @@ namespace cpb {
 class CompressedSublattices {
 public:
     struct Element {
-        SubAliasID alias_id; ///< the alias ID of each sublattice (unique among elements)
+        SiteID id; ///< the alias ID of each sublattice (unique among elements)
         storage_idx_t num_sites; ///< number of sublattice sites in the final system
         storage_idx_t num_orbitals; ///< number of orbitals on this sublattice
     };
@@ -31,7 +31,7 @@ class CompressedSublattices {
         It(std::vector<Element>::const_iterator it) : it(it) {}
 
         /// Directly correspond to fields in `Element`
-        SubAliasID alias_id() const { return it->alias_id; }
+        SiteID id() const { return it->id; }
         idx_t num_sites() const { return it->num_sites; }
         idx_t num_orbitals() const { return it->num_orbitals; }
 
@@ -70,7 +70,7 @@ class CompressedSublattices {
                           ArrayXi const& orbital_counts);
 
     /// Start a new sublattice block or increment the site count for the existing block
-    void add(SubAliasID id, idx_t norb);
+    void add(SiteID id, idx_t norb);
 
     /// Remove sites for which `keep == false`
     void filter(VectorX<bool> const& keep);