Program Listing for File bz.hpp

↰ Return to documentation for file (src/bz.hpp)

/* This file is part of brille.

Copyright © 2019,2020 Greg Tucker <greg.tucker@stfc.ac.uk>

brille is free software: you can redistribute it and/or modify it under the
terms of the GNU Affero General Public License as published by the Free
Software Foundation, either version 3 of the License, or (at your option)
any later version.

brille 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 Affero General Public License for more details.

You should have received a copy of the GNU Affero General Public License
along with brille. If not, see <https://www.gnu.org/licenses/>.            */
#ifndef BRILLE_BZ_CLASS_H_
#define BRILLE_BZ_CLASS_H_
#include <omp.h>
#include "neighbours.hpp"
#include "transform.hpp"
#include "polyhedron.hpp"
#include "phonon.hpp"
#include "approx.hpp"
namespace brille {

class BrillouinZone {
  Reciprocal lattice;
  Reciprocal outerlattice;
  Polyhedron polyhedron;
  Polyhedron ir_polyhedron;
  bArray<double> ir_wedge_normals;
  bool time_reversal;
  bool has_inversion;
  bool is_primitive;
  bool no_ir_mirroring;
public:
  BrillouinZone(const Reciprocal& lat,
                const bool toprim=true,
                const int extent=1,
                const bool tr=false,
                const bool wedge_search=true
               ):
  lattice(toprim ? lat.primitive() : lat), outerlattice(lat), time_reversal(tr)
  {
    profile_update("Start of BrillouinZone construction");
    this->is_primitive = !(this->lattice.issame(this->outerlattice));
    this->has_inversion = this->time_reversal || lat.has_space_inversion();
    this->no_ir_mirroring = true;
    double old_volume = -1.0, new_volume=0.0;
    int test_extent = extent-1;
    // initial test_extent based on spacegroup or pointgroup?
    while (!brille::approx::scalar(new_volume, old_volume)){
      old_volume = new_volume;
      this->voro_search(++test_extent);
      new_volume = this->polyhedron.get_volume();
    }
    // fallback in case voro_search fails for some reason?!?
    if (brille::approx::scalar(new_volume, 0.)){
      throw std::runtime_error("voro_search failed to produce a non-null first Brillouin zone.");
    } else {
      verbose_update("New polyhedron volume ", this->polyhedron.get_volume());
      verbose_update("First Brillouin zone found using extent ",test_extent,", a ",this->polyhedron.string_repr());
    }
    // in case we've been asked to perform a wedge search for, e.g., P1 or P-1,
    // set the irreducible wedge now as the search will do nothing.
    this->ir_polyhedron = this->polyhedron;
    if (wedge_search){
      bool success = this->wedge_brute_force()
                  || this->wedge_brute_force(false,false) // no special 2-fold or mirror handling
                  || this->wedge_brute_force(false,true) // no special 2-fold handling (but special mirror handling)
                  || this->wedge_brute_force(true, false) // no special mirror handling (maybe not useful)
                  || this->wedge_brute_force(true, true, false) // last ditch effort, handle non order(2) operations in decreasing order
                  || this->wedge_brute_force(true, false, true, false)
                  ;
      // other combinations of special_2_folds, special_mirrors,
      // and sort_by_length are possible but not necessarily useful.
      if (!success)
        throw std::runtime_error("Failed to find an irreducible Brillouin zone.");
    }
    profile_update("  End of BrillouinZone construction");
  }
  void check_if_mirroring_needed(void){
    this->no_ir_mirroring = true;
    if (!this->has_inversion){
      PointSymmetry ps = this->outerlattice.get_pointgroup_symmetry(this->time_reversal?1:0);
      double goal = this->polyhedron.get_volume() / static_cast<double>(ps.size());
      double found = this->ir_polyhedron.get_volume();
      if (brille::approx::scalar(goal, 2.0*found)){
        /*The current Polyhedron at this->ir_polyhedron has half the anticipated
        volume. We can 'fix' this the easy way by mirroring the Polyhedron and
        gluing it onto the current version; the resultant polyhedron will come
        to a point at (0,0,0) and will not be convex. Instead, try to be clever
        and look for a convex combination of the found polyhedron and its
        mirror with one of the pointgroup operations applied:*/
        std::vector<Polyhedron> all_unions;
        Polyhedron mirrored = this->ir_polyhedron.mirror();
        for (auto & r: ps.getall())
          all_unions.push_back(this->ir_polyhedron + mirrored.rotate(r));
        // The combination of a polyhedron and its rotated inverse which has
        // the least number of vertices is (hopefully) convex.
        auto min_vert_union = std::min_element(
          all_unions.begin(),all_unions.end(),
          [](const Polyhedron & a, const Polyhedron & b){
            return a.num_vertices() < b.num_vertices();
          }
        );
        // Polyhedron::operator+() needs to be fixed. As of now it does not look
        // for extraneous faces (like internal faces which are coplanar and
        // pointing in opposite directions) or coplanar external faces which
        // share an edge. Until this is fixed cross our fingers and hope that we
        // created a convex polyhedron such that the convex hull of its points
        // gives the same polyhedron back:
        Polyhedron mvu_convex_hull(min_vert_union->get_vertices());
        if (brille::approx::scalar(goal, mvu_convex_hull.get_volume())){
          // we found a polyhedron with the right volume which is convex!
          // so we can keep this as *the* ir_polyhedron
          this->ir_polyhedron = mvu_convex_hull;
          // and we need to update the found volume for the check below
          found = goal;
        }
      }
      // if found == goal no mirroring is required.
      // if found == goal/2, mirroring is required.
      this->no_ir_mirroring = brille::approx::scalar(goal, found);
    }
  }

  bool check_ir_polyhedron(void);
  bool wedge_explicit(void);
  const Reciprocal get_lattice() const { return this->outerlattice;};
  const Reciprocal get_primitive_lattice() const { return this->lattice;};
  size_t vertices_count() const { return this->get_vertices().size(0);};
  size_t faces_count() const { return this->get_points().size(0);};
  // irreducible reciprocal space wedge
  void
  set_ir_wedge_normals(const LQVec<double>& x){
    bool already_same = this->outerlattice.issame(x.get_lattice());
    LQVec<double> xp(this->outerlattice);
    PrimitiveTransform PT(this->outerlattice.get_bravais_type());
    bool transform_needed = ( PT.does_anything() && this->lattice.issame(x.get_lattice()) );
    if (!(already_same || transform_needed))
      throw std::runtime_error("ir_wedge_normals must be in the standard or primitive lattice used to define the BrillouinZone object");
    if (transform_needed)  xp = transform_from_primitive(this->outerlattice,x);
    const LQVec<double> & xref = transform_needed ? xp : x;
    this->ir_wedge_normals = xref.get_hkl();
  }
  LQVec<double> get_ir_wedge_normals() const;
  LQVec<double> get_primitive_ir_wedge_normals() const;
  template<typename... A>
  void
  set_polyhedron(const LQVec<double>& v, const LQVec<double>& p, A... args){
    bool both_same = v.get_lattice().issame(p.get_lattice());
    if (!both_same)
      throw std::runtime_error("The vertices, and points of a polyhedron must all be in the same cooridinate system");
    bool is_outer = this->outerlattice.issame(v.get_lattice());
    bool is_inner = this->lattice.issame(v.get_lattice());
    LQVec<double> vp(this->outerlattice), pp(this->outerlattice);
    PrimitiveTransform PT(this->outerlattice.get_bravais_type());
    is_inner &= PT.does_anything();
    if (!(is_outer || is_inner))
      throw std::runtime_error("The polyhedron must be described in the conventional or primitive lattice used to define the BrillouinZone object");
    if (is_inner){
      vp = transform_from_primitive(this->outerlattice, v);
      pp = transform_from_primitive(this->outerlattice, p);
    }
    const LQVec<double> & vref = is_inner ? vp : v;
    const LQVec<double> & pref = is_inner ? pp : p;
    this->polyhedron = Polyhedron(vref.get_xyz(), pref.get_xyz(), args...);
  }
  template<typename... A>
  void
  set_ir_polyhedron(const LQVec<double>& v, const LQVec<double>& p, const LQVec<double>& n, A... args){
    bool all_same = v.get_lattice().issame(p.get_lattice()) && p.get_lattice().issame(n.get_lattice());
    if (!all_same)
      throw std::runtime_error("The vertices, points, and normals of a polyhedron must all be in the same cooridinate system");
    bool is_outer = this->outerlattice.issame(v.get_lattice());
    bool is_inner = this->lattice.issame(v.get_lattice());
    LQVec<double> vp(this->outerlattice), pp(this->outerlattice), np(this->outerlattice);
    PrimitiveTransform PT(this->outerlattice.get_bravais_type());
    is_inner &= PT.does_anything();
    if (!(is_outer || is_inner))
      throw std::runtime_error("The polyhedron must be described in the conventional or primitive lattice used to define the BrillouinZone object");
    if (is_inner){
      vp = transform_from_primitive(this->outerlattice, v);
      pp = transform_from_primitive(this->outerlattice, p);
      np = transform_from_primitive(this->outerlattice, n);
    }
    const LQVec<double> & vref = is_inner ? vp : v;
    const LQVec<double> & pref = is_inner ? pp : p;
    const LQVec<double> & nref = is_inner ? np : n;
    this->ir_polyhedron = Polyhedron(vref.get_xyz(), pref.get_xyz(), nref.get_xyz(), args...);
  }
  bool set_ir_vertices(const LQVec<double>& v){
    bool is_outer = this->outerlattice.issame(v.get_lattice());
    bool is_inner = this->lattice.issame(v.get_lattice());
    LQVec<double> vp(this->outerlattice);
    PrimitiveTransform PT(this->outerlattice.get_bravais_type());
    is_inner &= PT.does_anything();
    if (!(is_outer || is_inner))
      throw std::runtime_error("The polyhedron must be described in the conventional or primitive lattice used to define the BrillouinZone object");
    if (is_inner)
      vp = transform_from_primitive(this->outerlattice, v);
    const LQVec<double> & vref = is_inner ? vp : v;
    debug_update("Generate a convex polyhedron from\n", vref.to_string());
    this->ir_polyhedron = Polyhedron(vref.get_xyz());
    debug_update("Generated a ", this->ir_polyhedron.string_repr()," from the specified vertices");
    // check that the polyhedron we've specified has the desired properties
    if (this->check_ir_polyhedron()){
      // set the wedge normals as well
      this->set_ir_wedge_normals(this->get_ir_polyhedron_wedge_normals());
      return true;
    }
    return false;
  }
  Polyhedron get_polyhedron(void) const;
  LQVec<double> get_vertices(void) const;
  LQVec<double> get_points(void) const;
  LQVec<double> get_normals(void) const;
  LQVec<double> get_half_edges(void) const;
  LQVec<double> get_primitive_vertices(void) const;
  LQVec<double> get_primitive_points(void) const;
  LQVec<double> get_primitive_normals(void) const;
  std::vector<std::vector<int>> get_faces_per_vertex(void) const;
  std::vector<std::vector<int>> get_vertices_per_face(void) const;
  Polyhedron get_ir_polyhedron(const bool true_ir=true) const;
  LQVec<double> get_ir_vertices(void) const;
  LQVec<double> get_ir_points(void) const;
  LQVec<double> get_ir_normals(void) const;
  LQVec<double> get_ir_primitive_vertices(void) const;
  LQVec<double> get_ir_primitive_points(void) const;
  LQVec<double> get_ir_primitive_normals(void) const;
  std::vector<std::vector<int>> get_ir_faces_per_vertex(void) const;
  std::vector<std::vector<int>> get_ir_vertices_per_face(void) const;
  void print() const;
  void wedge_search(const bool prefer_basis_vectors=true, const bool parallel_ok=true);
  bool wedge_brute_force(bool special_2_folds = true, bool special_mirrors = true, bool sort_by_length=true, bool sort_one_sym=true);
  void wedge_triclinic(void);
  void irreducible_vertex_search();
  void voro_search(const int extent=1);
  bool isprimitive(void) const {return this->is_primitive;};

  template<class T>
  std::vector<bool>
  isinside(const LQVec<T>& p) const {
    bool isouter = this->outerlattice.issame(p.get_lattice());
    bool isinner = this->lattice.issame(p.get_lattice());
    if (!(isouter||isinner))
      throw std::runtime_error("Q points must be in the standard or primitive lattice");
    std::vector<bool> out(p.size(0), true);
    LQVec<double> points, normals;
    if (isouter){
      points = this->get_points();
      normals = this->get_normals();
    } else {
      points = this->get_primitive_points();
      normals = this->get_primitive_normals();
    }
    for (size_t i=0; i<p.size(0); ++i)
      out[i] = dot(normals, p.view(i)-points).all(brille::cmp::le, 0.);
    return out;
  }
  template<class T>
  std::vector<bool>
  isinside_wedge(const LQVec<T> &p, const bool constructing=false) const {
    bool isouter = this->outerlattice.issame(p.get_lattice());
    bool isinner = this->lattice.issame(p.get_lattice());
    if (!(isouter||isinner)){
      std::string msg = "Q points provided to BrillouinZone::isinside_wedge ";
      msg += "must be in the standard or primitive lattice ";
      msg += "used to define the BrillouinZone object";
      throw std::runtime_error(msg);
    }
    std::vector<bool> out(p.size(0), true);
    LQVec<double> normals;
    if (isouter)
      normals = this->get_ir_wedge_normals();
    else
      normals = this->get_primitive_ir_wedge_normals();
    if (normals.size(0)){ // with no normals *all* points are "inside" the wedge
      // If a pointgroup has inversion symmetry then for every point, p, there is
      // an equivalent point, -p. This indicates that a point p is already in
      // the irreducible wedge if it has n̂ᵢ⋅p ≥ 0 for all irredudible-bounding-
      // plane normals, ̂nᵢ, *or* the opposite -- n̂ᵢ⋅ p ≤ 0 for all ̂nᵢ.
      // Since we are interested in enforcing a smallest-possible irreducible
      // Brillouin zone, we want to exclude the all n̂ᵢ⋅ p ≤ 0 half of the
      // reciprocal wedge precisely because they are equivalent.
      // It is only in the case where a pointgroup *does not have* space-inversion
      // symmetry and time-reversal symmetry is to be excluded that we can allow
      // n̂ᵢ⋅ p ≤ 0 to be a valid in-wedge solution.
      // The Array all method has a special switched execution path
      // for checking whether all values are ≤ or ≥ a value simultaneously

      // when constructing the irreducible Brillouin zone we need to only consider
      // the ≥0 case so that we end up with a convex polyhedron. The ir_polyhedron
      // accessor method mirrors the half-polyhedron in this case, so when
      // identifying whether a point is inside of the irreducible Brillouin zone
      // we must allow for the ≤0 case as well.
      brille::cmp c = constructing||this->no_ir_mirroring ? brille::cmp::ge : brille::cmp::le_ge;
      for (size_t i=0; i<p.size(0); ++i)
        out[i] = dot(normals, p.view(i)).all(c, 0.);
    }
    return out;
  }
  bool
  moveinto(const LQVec<double>& Q, LQVec<double>& q, LQVec<int>& tau, const int threads=0) const
  {
    profile_update("BrillouinZone::moveinto called with ",threads," threads");
    omp_set_num_threads( (threads > 0) ? threads : omp_get_max_threads() );
    bool already_same = this->lattice.issame(Q.get_lattice());
    LQVec<double> Qprim(this->lattice);
    LQVec<double> qprim(this->lattice);
    LQVec<int> tauprim(this->lattice);
    PrimitiveTransform PT(this->outerlattice.get_bravais_type());
    bool transform_needed = ( PT.does_anything() && this->outerlattice.issame(Q.get_lattice()) );
    if (!(already_same || transform_needed)){
      std::string msg = "Q points provided to BrillouinZone::moveinto must be ";
      msg += "in the standard or primitive lattice used to define ";
      msg += "the BrillouinZone object";
      throw std::runtime_error(msg);
    }
    if (transform_needed)  Qprim = transform_to_primitive(this->outerlattice,Q);
    const LQVec<double> & Qsl = transform_needed ? Qprim : Q;
    LQVec<double> & qsl = transform_needed ? qprim : q;
    LQVec<int> & tausl = transform_needed? tauprim : tau;

    // the face centre points and normals in the primitive lattice:
    auto points = this->get_primitive_points();
    auto normals = this->get_primitive_normals();
    normals = normals/norm(normals); // ensure they're normalised
    auto taus = (2.0*points).round();
    auto taulen = norm(taus);
    size_t max_count = taus.size(0);
    // ensure that qsl and tausl can hold each qi and taui
    qsl.resize(Qsl.size(0));
    tausl.resize(Qsl.size(0));
    long long snQ = brille::utils::u2s<long long, size_t>(Qsl.size(0));
  #pragma omp parallel for default(none)\
  shared(Qsl, tausl, qsl, points, normals, taus, taulen, snQ, max_count)\
  schedule(dynamic)
    for (long long si=0; si<snQ; si++){
      size_t i = brille::utils::s2u<size_t, long long>(si);
      auto taui = Qsl.view(i).round();
      auto qi = Qsl.view(i) - taui;
      auto last_shift = taui;
      size_t count{0};
      while (count++ < max_count && dot(normals, qi-points).any(brille::cmp::gt,0.)){
        auto qi_dot_normals = dot(qi , normals);
        auto Nhkl = (qi_dot_normals/taulen).round().to_std();
        auto qidn = qi_dot_normals.to_std();
        if (std::any_of(Nhkl.begin(), Nhkl.end(), [](int a){return a > 0;})){
          int maxnm{0};
          size_t maxat{0};
          for (size_t j=0; j<Nhkl.size(); ++j)
                                                           // protect against oscillating by ±τ
          if (Nhkl[j]>0 && Nhkl[j]>=maxnm && (0==maxnm || (norm(taus.view(j)+last_shift).all(brille::cmp::gt, 0.) && qidn[j]>qidn[maxat]))){
            maxnm = Nhkl[maxat=j];
          }
          qi -= taus.view(maxat) * static_cast<double>(maxnm); // ensure we subtract LQVec<double>
          taui += taus.view(maxat) * maxnm; // but add LQVec<int>
          last_shift = taus.view(maxat) * maxnm;
        }
      }
      qsl.set(i, qi);
      tausl.set(i, taui);
    }
    if (transform_needed){ // then we need to transform back q and tau
      q   = transform_from_primitive(this->outerlattice,qsl);
      tau = transform_from_primitive(this->outerlattice,tausl);
    }
    auto allinside = this->isinside(q);
    if (std::count(allinside.begin(), allinside.end(), false) > 0){
      for (size_t i=0; i<Q.size(0); ++i) if (!allinside[i]){
        info_update("Q  =",Q.to_string(i)  ," tau  =",tau.to_string(i)  ," q  =",q.to_string(i));
        info_update("Qsl=",Qsl.to_string(i)," tausl=",tausl.to_string(i)," qsl=",qsl.to_string(i),"\n");
      }
      throw std::runtime_error("Not all points inside Brillouin zone");
      // return false;
    }
    return true; // otherwise an error has been thrown
  }
  bool
  ir_moveinto(
    const LQVec<double>& Q, LQVec<double>& q, LQVec<int>& tau,
    std::vector<size_t>& Ridx, std::vector<size_t>& invRidx, const int threads=0)
  const {
    profile_update("BrillouinZone::ir_moveinto called with ",threads," threads");
    omp_set_num_threads( (threads > 0) ? threads : omp_get_max_threads() );
    /* The Pointgroup symmetry information comes from, effectively, spglib which
    has all rotation matrices defined in the conventional unit cell -- which is
    our `outerlattice`. Consequently we must work in the outerlattice here.  */
    if (!this->outerlattice.issame(Q.get_lattice()))
      throw std::runtime_error("Q points provided to ir_moveinto must be in the standard lattice used to define the BrillouinZone object");
    // get the PointSymmetry object, containing all operations
    PointSymmetry psym = this->get_pointgroup_symmetry();
    // ensure q, tau, and Rm can hold one for each Q.
    size_t nQ = Q.size(0);
    auto Qshape = Q.shape();
    q.resize(Qshape);
    tau.resize(Qshape);
    Ridx.resize(nQ);
    invRidx.resize(nQ);
    // find q₁ₛₜ in the first Brillouin zone and τ ∈ [reciprocal lattice vectors]
    // such that Q = q₁ₛₜ + τ
    this->moveinto(Q, q, tau, threads);
    // by chance some first Bz points are likely already in the IR-Bz:
    std::vector<bool> in_ir = this->isinside_wedge(q);
    //LQVec<double> qj(Q.get_lattice(), 1u);
    auto lat = Q.get_lattice();
    // OpenMP 2 (VS) doesn't like unsigned loop counters
    size_t n_outside{0};
    long long snQ = brille::utils::u2s<long long, size_t>(nQ);
    #pragma omp parallel for default(none) shared(psym, Ridx, invRidx, q, in_ir, lat, snQ) reduction(+:n_outside) schedule(dynamic)
    for (long long si=0; si<snQ; ++si){
      size_t i = brille::utils::s2u<size_t, long long>(si);
      // any q already in the irreducible zone need no rotation → identity, but we need to find the index of E
      bool outside=!in_ir[i];
      if (outside){
        // for others find the jᵗʰ operation which moves qᵢ into the irreducible zone
        LQVec<double> qj(lat, 1u); // a place to hold the multiplication result
        for (size_t j=0; j<psym.size(); ++j) if (outside) {
          // The point symmetry matrices relate *real space* vectors! We must use
          // their transposes' to rotate reciprocal space vectors.
          brille::utils::multiply_matrix_vector(qj.ptr(0), transpose(psym.get(j)).data(), q.ptr(i));
          if ( this->isinside_wedge(qj)[0] ){ /* store the result */
            // and (Rⱼᵀ)⁻¹ ∈ G, such that Qᵢ = (Rⱼᵀ)⁻¹⋅qᵢᵣ + τᵢ.
            q.set(i, qj); // keep Rⱼᵀ⋅qᵢ as qᵢᵣ
            invRidx[i] = j; // Rⱼ *is* the inverse of what we want for ouput
            Ridx[i] = psym.get_inverse_index(j); // find the index of Rⱼ⁻¹
            outside = false;
          }
        }
      } else {
        invRidx[i] = Ridx[i] = psym.find_index({1,0,0, 0,1,0, 0,0,1});
      }
      if (outside) {
        ++n_outside;
        in_ir[i] = false;
      }
    }
    if (n_outside > 0) for (size_t i=0; i<nQ; ++i) if (!in_ir[i]){
      std::string msg = "Q = " + Q.to_string(i);
      msg += " is outside of the irreducible BrillouinZone ";
      msg += " : tau = " + tau.to_string(i) + " , q = " + q.to_string(i);
      throw std::runtime_error(msg);
      return false;
    }
    return true; // otherwise we hit the runtime error above
  }
  bool
  ir_moveinto_wedge(const LQVec<double>& Q, LQVec<double>& q, std::vector<std::array<int,9>>& R, const int threads=0) const {
    omp_set_num_threads( (threads > 0) ? threads : omp_get_max_threads() );
    /* The Pointgroup symmetry information comes from, effectively, spglib which
    has all rotation matrices defined in the conventional unit cell -- which is
    our `outerlattice`. Consequently we must work in the outerlattice here.  */
    if (!this->outerlattice.issame(Q.get_lattice()))
      throw std::runtime_error("Q points provided to ir_moveinto must be in the standard lattice used to define the BrillouinZone object");
    // get the PointSymmetry object, containing all operations
    PointSymmetry psym = this->outerlattice.get_pointgroup_symmetry(this->time_reversal);
    // ensure q and R can hold one for each Q.
    size_t nQ = Q.size(0);
    auto Qshape = Q.shape();
    q.resize(Qshape);
    R.resize(nQ);
    // by chance some first Bz points are likely already in the IR-wedge:
    std::vector<bool> in_ir = this->isinside_wedge(Q);
    //LQVec<double> qj(Q.get_lattice(), 1u);
    auto lat = Q.get_lattice();
    // OpenMP 2 (VS) doesn't like unsigned loop counters
    size_t n_outside{0};
    long long snQ = brille::utils::u2s<long long, size_t>(nQ);
    #pragma omp parallel for default(none) shared(psym, R, q, Q, in_ir, lat, snQ) reduction(+:n_outside) schedule(dynamic)
    for (long long si=0; si<snQ; ++si){
      size_t i = brille::utils::s2u<size_t, long long>(si);
      // any q already in the irreducible zone need no rotation → identity
      bool outside=!in_ir[i];
      if (outside){
        // for others find the jᵗʰ operation which moves qᵢ into the irreducible zone
        LQVec<double> qj(lat, 1u); // a place to hold the multiplication result
        for (size_t j=0; j<psym.size(); ++j) if (outside) {
          // The point symmetry matrices relate *real space* vectors! We must use
          // their transposes' to rotate reciprocal space vectors.
          brille::utils::multiply_matrix_vector(qj.ptr(0), transpose(psym.get(j)).data(), Q.ptr(i));
          if ( this->isinside_wedge(qj)[0] ){ /* store the result */
            q.set(i, qj); // keep Rⱼᵀ⋅Qᵢ as qᵢᵣ
            R[i] = transpose(psym.get_inverse(j)); // and (Rⱼᵀ)⁻¹ ∈ G, such that Q = (Rⱼᵀ)⁻¹⋅qᵢᵣ
            outside = false;
          }
        }
      } else {
        q.set(i, Q.view(i));
        R[i] = {1,0,0, 0,1,0, 0,0,1};
      }
      if (outside) {
        ++n_outside;
        in_ir[i] = false;
      }
    }
    if (n_outside > 0) for (size_t i=0; i<nQ; ++i) if (!in_ir[i]){
      std::string msg = "Q = " + Q.to_string(i);
      msg += " is outside of the irreducible reciprocal space wedge ";
      msg += " , irQ = " + q.to_string(i);
      throw std::runtime_error(msg);
      return false;
    }
    return true; // otherwise we hit the runtime error above
  }

  PointSymmetry get_pointgroup_symmetry() const{
    return this->outerlattice.get_pointgroup_symmetry(this->time_reversal);
  }
  int add_time_reversal() const {
    return this->time_reversal ? 1 : 0;
  }
private:
  template<class T, class R>
  bool
  wedge_normal_check(const LQVec<T>& n, LQVec<R>& normals, size_t& num){
    std::string msg = "Considering " + n.to_string(0) + "... ";
    if (norm(n).all(brille::cmp::eq, 0.0)){
      debug_update(msg, "rejected; zero-length");
      return false;
    }
    if (num==0){
      debug_update(msg, "accepted; first normal");
      normals.set(num, n.view(0));
      num=num+1;
      return true;
    }
    // num > 0 for all following views
    if (norm(cross(normals.view(0,num), n)).any(brille::cmp::eq, 0.)){
      debug_update(msg, "rejected; already present");
      return false;
    }
    normals.set(num,  n.view(0));
    if (this->ir_wedge_is_ok(normals.view(0,num+1))){
      debug_update(msg, "accepted");
      num=num+1;
      return true;
    }
    normals.set(num, -n.extract(0));
    if (this->ir_wedge_is_ok(normals.view(0,num+1))){
      debug_update(msg, "accepted (*-1)");
      num=num+1;
      return true;
    }
    debug_update(msg, "rejected; addition causes null wedge");
    return false;
  }
  template<class T, class R, class U>
  bool
  wedge_normal_check(const LQVec<T>& n0, const LQVec<R>& n1, LQVec<U>& normals, size_t& num){
    std::string msg = "Considering " + n0.to_string(0)+ " and " + n1.to_string(0) + "... ";
    bool p0=false, p1=false;
    if (num>0){
      p0 = norm(cross(normals.view(0,num), n0)).any(brille::cmp::eq,0.);
      p1 = norm(cross(normals.view(0,num), n1)).any(brille::cmp::eq,0.);
    }
    if (p0 && p1){
      debug_update(msg, "rejected; already present");
      return false;
    }
    if (num>0 && dot(normals.view(0,num)/norm(n0),n0/norm(n0)).any(brille::cmp::eq,1.)){
      normals.set(num,  n1.view(0));
      if (this->ir_wedge_is_ok(normals.view(0,num+1))){
        debug_update(msg, "n1 accepted (n0 present)");
        num=num+1;
        return true;
      }
      debug_update(msg, "n1 rejected (n0 present); addition causes null wedge");
      return false;
    }
    if (num>0 && dot(normals.view(0,num)/norm(n1),n1/norm(n1)).any(brille::cmp::eq,1.)){
      normals.set(num,  n0.view(0));
      if (this->ir_wedge_is_ok(normals.view(0,num+1))){
        debug_update(msg, "n0 accepted (n1 present)");
        num=num+1;
        return true;
      }
      debug_update(msg, "n0 rejected (n1 present); addition causes null wedge");
      return false;
    }
    if (num>0 && (p0 || p1)){
      debug_update_if(p0, msg, "-n0 present; addition causes null wedge)");
      debug_update_if(p1, msg, "-n1 present; addition causes null wedge)");
      return false;
    }
    normals.set(num,   n0.view(0));
    normals.set(num+1, n1.view(0));
    if (this->ir_wedge_is_ok(normals.view(0,num+2))){
      debug_update(msg, "n0 & n1 accepted");
      num=num+2;
      return true;
    }
    debug_update(msg, "n0 & n1 rejected; adding both would cause null wedge");
    return false;
  }
  void
  shrink_and_prune_outside(const size_t cnt, LQVec<double>& vrt, bArray<int>& ijk) const {
    verbose_update("shrinking to ",cnt);
    if(vrt.size(0) && ijk.size(0)){
      vrt.resize(cnt);
      ijk.resize(cnt);
      if (cnt){ // isinside has a problem with vrt.size()==0
        auto isin = this->isinside(vrt);
        verbose_update("and retaining ", std::count(isin.begin(), isin.end(), true), " inside vertices");
        vrt = vrt.extract(isin);
        ijk = ijk.extract(isin);
      }
    }
  }
  template<class T>
  bool
  ir_wedge_is_ok(const LQVec<T>& normals){
    this->set_ir_wedge_normals(normals); // assigns this->ir_wedge_normals
    this->irreducible_vertex_search(); // assigns this->ir_polyhedron
    return !brille::approx::scalar(this->ir_polyhedron.get_volume(), 0.0);
  }
  LQVec<double> get_ir_polyhedron_wedge_normals(void) const;
};

bool
intersect_at(const LQVec<double>& ni, const LQVec<double>& pi,
             const LQVec<double>& nj, const LQVec<double>& pj,
             const LQVec<double>& nk, const LQVec<double>& pk,
             LQVec<double>& intersect, const int idx);
bool
intersect_at(const LQVec<double>& ni, const LQVec<double>& pi,
             const LQVec<double>& nj, const LQVec<double>& pj,
             const LQVec<double>& nk,
             LQVec<double>& intersect, const int idx);
bool
intersect_at(const LQVec<double>& ni, const LQVec<double>& pi,
             const LQVec<double>& nj,
             const LQVec<double>& nk,
             LQVec<double>& intersect, const int idx);
bool
intersect_at(const LQVec<double>& ni,
             const LQVec<double>& nj,
             const LQVec<double>& nk,
             LQVec<double>& intersect, const int idx);
double
normals_matrix_determinant(const LQVec<double>& a, const LQVec<double>&b, const LQVec<double>& c);

} // end namespace brille
#endif