diff --git a/ChangeLog.md b/ChangeLog.md
index c5b22f551..d5bb22dae 100644
--- a/ChangeLog.md
+++ b/ChangeLog.md
@@ -8,6 +8,7 @@ and this project adheres to
 
 ### Added
 * Coloring plots by polytope area in `freud.locality.Voronoi`.
+* Neighbor vectors to `freud.locality.NeighborList`s.
 
 ### Changed
 * The `normalize` argument to `freud.density.RDF` is now `normalization_mode`.
diff --git a/cpp/cluster/Cluster.cc b/cpp/cluster/Cluster.cc
index be3df3170..73e9d8f75 100644
--- a/cpp/cluster/Cluster.cc
+++ b/cpp/cluster/Cluster.cc
@@ -23,9 +23,9 @@ void Cluster::compute(const freud::locality::NeighborQuery* nq, const freud::loc
         nq, nq->getPoints(), num_points, qargs, nlist,
         [&dj](const freud::locality::NeighborBond& neighbor_bond) {
             // Merge the two sets using the disjoint set
-            if (!dj.same(neighbor_bond.point_idx, neighbor_bond.query_point_idx))
+            if (!dj.same(neighbor_bond.getPointIdx(), neighbor_bond.getQueryPointIdx()))
             {
-                dj.unite(neighbor_bond.point_idx, neighbor_bond.query_point_idx);
+                dj.unite(neighbor_bond.getPointIdx(), neighbor_bond.getQueryPointIdx());
             }
         });
 
diff --git a/cpp/density/CorrelationFunction.cc b/cpp/density/CorrelationFunction.cc
index 9692f3b1a..8071ffb16 100644
--- a/cpp/density/CorrelationFunction.cc
+++ b/cpp/density/CorrelationFunction.cc
@@ -88,11 +88,11 @@ void CorrelationFunction<T>::accumulate(const freud::locality::NeighborQuery* ne
     accumulateGeneral(
         neighbor_query, query_points, n_query_points, nlist, qargs,
         [&](const freud::locality::NeighborBond& neighbor_bond) {
-            size_t value_bin = m_histogram.bin({neighbor_bond.distance});
+            size_t value_bin = m_histogram.bin({neighbor_bond.getDistance()});
             m_local_histograms.increment(value_bin);
             m_local_correlation_function.increment(
                 value_bin,
-                product(values[neighbor_bond.point_idx], query_values[neighbor_bond.query_point_idx]));
+                product(values[neighbor_bond.getPointIdx()], query_values[neighbor_bond.getQueryPointIdx()]));
         });
 }
 
diff --git a/cpp/density/LocalDensity.cc b/cpp/density/LocalDensity.cc
index 68cdb1ee7..e66aa9043 100644
--- a/cpp/density/LocalDensity.cc
+++ b/cpp/density/LocalDensity.cc
@@ -43,7 +43,7 @@ void LocalDensity::compute(const freud::locality::NeighborQuery* neighbor_query,
             for (freud::locality::NeighborBond nb = ppiter->next(); !ppiter->end(); nb = ppiter->next())
             {
                 // count particles that are fully in the r_max sphere
-                if (nb.distance < (m_r_max - m_diameter / float(2.0)))
+                if (nb.getDistance() < (m_r_max - m_diameter / float(2.0)))
                 {
                     num_neighbors += float(1.0);
                 }
@@ -54,7 +54,7 @@ void LocalDensity::compute(const freud::locality::NeighborQuery* neighbor_query,
                     // lots of them. It smooths out the neighbor count distributions and avoids noisy spikes
                     // that obscure data
                     num_neighbors
-                        += float(1.0) + (m_r_max - (nb.distance + m_diameter / float(2.0))) / m_diameter;
+                        += float(1.0) + (m_r_max - (nb.getDistance() + m_diameter / float(2.0))) / m_diameter;
                 }
                 m_num_neighbors_array[i] = num_neighbors;
                 if (m_box.is2D())
diff --git a/cpp/density/RDF.cc b/cpp/density/RDF.cc
index ab09fc5dd..cf59e8e7d 100644
--- a/cpp/density/RDF.cc
+++ b/cpp/density/RDF.cc
@@ -88,7 +88,7 @@ void RDF::accumulate(const freud::locality::NeighborQuery* neighbor_query, const
 {
     accumulateGeneral(neighbor_query, query_points, n_query_points, nlist, qargs,
                       [&](const freud::locality::NeighborBond& neighbor_bond) {
-                          m_local_histograms(neighbor_bond.distance);
+                          m_local_histograms(neighbor_bond.getDistance());
                       });
 }
 
diff --git a/cpp/environment/BondOrder.cc b/cpp/environment/BondOrder.cc
index fbe88681b..aed715074 100644
--- a/cpp/environment/BondOrder.cc
+++ b/cpp/environment/BondOrder.cc
@@ -85,9 +85,9 @@ void BondOrder::accumulate(const locality::NeighborQuery* neighbor_query, quat<f
 {
     accumulateGeneral(neighbor_query, query_points, n_query_points, nlist, qargs,
                       [&](const freud::locality::NeighborBond& neighbor_bond) {
-                          const quat<float>& ref_q(orientations[neighbor_bond.point_idx]);
-                          vec3<float> v(bondVector(neighbor_bond, neighbor_query, query_points));
-                          const quat<float>& q = query_orientations[neighbor_bond.query_point_idx];
+                          const quat<float>& ref_q(orientations[neighbor_bond.getPointIdx()]);
+                          vec3<float> v(neighbor_bond.getVector());
+                          const quat<float>& q(query_orientations[neighbor_bond.getQueryPointIdx()]);
                           if (m_mode == obcd)
                           {
                               // give bond directions of neighboring particles rotated by the matrix
diff --git a/cpp/environment/LocalBondProjection.cc b/cpp/environment/LocalBondProjection.cc
index 7d99435c8..3577b53d5 100644
--- a/cpp/environment/LocalBondProjection.cc
+++ b/cpp/environment/LocalBondProjection.cc
@@ -71,11 +71,12 @@ void LocalBondProjection::compute(const locality::NeighborQuery* nq, const quat<
                 const size_t j(m_nlist.getNeighbors()(bond, 1));
 
                 // compute bond vector between the two particles
-                vec3<float> local_bond(bondVector(locality::NeighborBond(i, j), nq, query_points));
+                vec3<float> local_bond(m_nlist.getVectors()[bond]);
+
                 // rotate bond vector into the local frame of particle p
                 local_bond = rotate(conj(orientations[j]), local_bond);
                 // store the length of this local bond
-                float local_bond_len = std::sqrt(dot(local_bond, local_bond));
+                float local_bond_len = m_nlist.getDistances()[bond];
 
                 for (unsigned int k = 0; k < n_proj; k++)
                 {
diff --git a/cpp/environment/LocalDescriptors.cc b/cpp/environment/LocalDescriptors.cc
index dce2b4c3d..13fead52d 100644
--- a/cpp/environment/LocalDescriptors.cc
+++ b/cpp/environment/LocalDescriptors.cc
@@ -52,8 +52,7 @@ void LocalDescriptors::compute(const locality::NeighborQuery* nq, const vec3<flo
                      && m_nlist.getNeighbors()(bond_copy, 0) == i && neighbor_count < max_num_neighbors;
                      ++bond_copy, ++neighbor_count)
                 {
-                    const size_t j(m_nlist.getNeighbors()(bond_copy, 1));
-                    const vec3<float> r_ij(bondVector(locality::NeighborBond(i, j), nq, query_points));
+                    const vec3<float> r_ij(m_nlist.getVectors()[bond_copy]);
                     const float r_sq(dot(r_ij, r_ij));
 
                     for (size_t ii(0); ii < 3; ++ii)
@@ -105,14 +104,11 @@ void LocalDescriptors::compute(const locality::NeighborQuery* nq, const vec3<flo
                  ++bond, ++neighbor_count)
             {
                 const unsigned int sphCount(bond * getSphWidth());
-                const size_t j(m_nlist.getNeighbors()(bond, 1));
-                const vec3<float> r_ij(bondVector(locality::NeighborBond(i, j), nq, query_points));
-                const float r_sq(dot(r_ij, r_ij));
+                const vec3<float> r_ij(m_nlist.getVectors()[bond]);
+                const float magR(m_nlist.getDistances()[bond]);
                 const vec3<float> bond_ij(dot(rotation_0, r_ij), dot(rotation_1, r_ij),
                                           dot(rotation_2, r_ij));
 
-                const float magR(std::sqrt(r_sq));
-
                 // Wrap theta into [0, 2*pi]
                 float theta(std::atan2(bond_ij.y, bond_ij.x));
                 theta = util::modulusPositive(theta, constants::TWO_PI);
diff --git a/cpp/environment/MatchEnv.cc b/cpp/environment/MatchEnv.cc
index e10be058f..3d51bc38b 100644
--- a/cpp/environment/MatchEnv.cc
+++ b/cpp/environment/MatchEnv.cc
@@ -7,7 +7,6 @@
 
 #include "MatchEnv.h"
 
-#include "NeighborBond.h"
 #include "NeighborComputeFunctional.h"
 
 namespace freud { namespace environment {
@@ -488,8 +487,7 @@ MatchEnv::~MatchEnv() = default;
 
 EnvironmentCluster::~EnvironmentCluster() = default;
 
-Environment MatchEnv::buildEnv(const freud::locality::NeighborQuery* nq,
-                               const freud::locality::NeighborList* nlist, size_t num_bonds, size_t& bond,
+Environment MatchEnv::buildEnv(const freud::locality::NeighborList* nlist, size_t num_bonds, size_t& bond,
                                unsigned int i, unsigned int env_ind)
 {
     Environment ei = Environment();
@@ -502,7 +500,7 @@ Environment MatchEnv::buildEnv(const freud::locality::NeighborQuery* nq,
         const size_t j(nlist->getNeighbors()(bond, 1));
         if (i != j)
         {
-            vec3<float> delta(bondVector(locality::NeighborBond(i, j), nq, nq->getPoints()));
+            vec3<float> delta(nlist->getVectors()[bond]);
             ei.addVec(delta);
         }
     }
@@ -540,7 +538,7 @@ void EnvironmentCluster::compute(const freud::locality::NeighborQuery* nq,
     // if you don't do this, things will get screwy.
     for (unsigned int i = 0; i < Np; i++)
     {
-        Environment ei = buildEnv(nq, &env_nlist, env_num_bonds, env_bond, i, i);
+        Environment ei = buildEnv(&env_nlist, env_num_bonds, env_bond, i, i);
         dj.s.push_back(ei);
         dj.m_max_num_neigh = std::max(dj.m_max_num_neigh, ei.num_vecs);
         ;
@@ -726,7 +724,7 @@ void EnvironmentMotifMatch::compute(const freud::locality::NeighborQuery* nq,
     for (unsigned int i = 0; i < Np; i++)
     {
         unsigned int dummy = i + 1;
-        Environment ei = buildEnv(nq, &nlist, num_bonds, bond, i, dummy);
+        Environment ei = buildEnv(&nlist, num_bonds, bond, i, dummy);
         dj.s.push_back(ei);
 
         // if the environment matches e0, merge it into the e0 environment set
@@ -803,7 +801,7 @@ void EnvironmentRMSDMinimizer::compute(const freud::locality::NeighborQuery* nq,
     for (unsigned int i = 0; i < Np; i++)
     {
         unsigned int dummy = i + 1;
-        Environment ei = buildEnv(nq, &nlist, num_bonds, bond, i, dummy);
+        Environment ei = buildEnv(&nlist, num_bonds, bond, i, dummy);
         dj.s.push_back(ei);
 
         // if the environment matches e0, merge it into the e0 environment set
diff --git a/cpp/environment/MatchEnv.h b/cpp/environment/MatchEnv.h
index 544a77fa5..148971550 100644
--- a/cpp/environment/MatchEnv.h
+++ b/cpp/environment/MatchEnv.h
@@ -205,8 +205,7 @@ class MatchEnv
 
     //! Construct and return a local environment surrounding the particle indexed by i. Set the environment
     //! index to env_ind.
-    static Environment buildEnv(const freud::locality::NeighborQuery* nq,
-                                const freud::locality::NeighborList* nlist, size_t num_bonds, size_t& bond,
+    static Environment buildEnv(const freud::locality::NeighborList* nlist, size_t num_bonds, size_t& bond,
                                 unsigned int i, unsigned int env_ind);
 
     //! Returns the entire Np by m_num_neighbors by 3 matrix of all environments for all particles
diff --git a/cpp/locality/AABBQuery.cc b/cpp/locality/AABBQuery.cc
index bb7b3a2f6..9938975a8 100644
--- a/cpp/locality/AABBQuery.cc
+++ b/cpp/locality/AABBQuery.cc
@@ -183,7 +183,7 @@ NeighborBond AABBQueryBallIterator::next()
                         // Check ii exclusion before including the pair.
                         if (r_sq < r_max_sq && r_sq >= r_min_sq)
                         {
-                            return NeighborBond(m_query_point_idx, j, std::sqrt(r_sq));
+                            return NeighborBond(m_query_point_idx, j, std::sqrt(r_sq), 1, r_ij);
                         }
                     }
                 }
@@ -230,7 +230,7 @@ NeighborBond AABBQueryIterator::next()
             // r_min filtering because we're querying beyond the normally safe
             // bounds, so we have to do it in this class.
             m_current_neighbors.clear();
-            m_all_distances.clear();
+            m_all_bonds_minimum_distance.clear();
             m_query_points_below_r_min.clear();
             std::shared_ptr<NeighborQueryPerPointIterator> ball_it = std::make_shared<AABBQueryBallIterator>(
                 static_cast<const AABBQuery*>(m_neighbor_query), m_query_point, m_query_point_idx,
@@ -243,26 +243,28 @@ NeighborBond AABBQueryIterator::next()
                     continue;
                 }
 
-                if (!m_exclude_ii || m_query_point_idx != nb.point_idx)
+                if (!m_exclude_ii || m_query_point_idx != nb.getPointIdx())
                 {
-                    nb.query_point_idx = m_query_point_idx;
+                    nb.setQueryPointIdx(m_query_point_idx);
                     // If we've expanded our search radius beyond safe
                     // distance, use the map instead of the vector.
                     if (m_search_extended)
                     {
-                        if ((m_all_distances.count(nb.point_idx) == 0)
-                            || m_all_distances[nb.point_idx] > nb.distance)
+                        const unsigned int nb_point_idx = nb.getPointIdx();
+                        const float nb_distance = nb.getDistance();
+                        if ((m_all_bonds_minimum_distance.count(nb_point_idx) == 0)
+                            || m_all_bonds_minimum_distance[nb_point_idx].getDistance() > nb_distance)
                         {
-                            m_all_distances[nb.point_idx] = nb.distance;
-                            if (nb.distance < m_r_min)
+                            m_all_bonds_minimum_distance[nb_point_idx] = nb;
+                            if (nb_distance < m_r_min)
                             {
-                                m_query_points_below_r_min.insert(nb.point_idx);
+                                m_query_points_below_r_min.insert(nb_point_idx);
                             }
                         }
                     }
                     else
                     {
-                        if (nb.distance >= m_r_min)
+                        if (nb.getDistance() >= m_r_min)
                         {
                             m_current_neighbors.emplace_back(nb);
                         }
@@ -281,18 +283,18 @@ NeighborBond AABBQueryIterator::next()
             }
 
             if ((m_r_cur >= m_r_max) || (m_r_cur >= max_plane_distance)
-                || ((m_all_distances.size() - m_query_points_below_r_min.size()) >= m_num_neighbors))
+                || ((m_all_bonds_minimum_distance.size() - m_query_points_below_r_min.size())
+                    >= m_num_neighbors))
             {
                 // Once this condition is reached, either we found enough
                 // neighbors beyond the normal min_plane_distance
                 // condition or we conclude that there are not enough
                 // neighbors left in the system.
-                for (const auto& bond_distance : m_all_distances)
+                for (const auto& minimum_distance_bond : m_all_bonds_minimum_distance)
                 {
-                    if (bond_distance.second >= m_r_min)
+                    if (minimum_distance_bond.second.getDistance() >= m_r_min)
                     {
-                        m_current_neighbors.emplace_back(m_query_point_idx, bond_distance.first,
-                                                         bond_distance.second);
+                        m_current_neighbors.emplace_back(minimum_distance_bond.second);
                     }
                 }
                 std::sort(m_current_neighbors.begin(), m_current_neighbors.end());
@@ -320,7 +322,7 @@ NeighborBond AABBQueryIterator::next()
     while ((m_count < m_num_neighbors) && (m_count < m_current_neighbors.size()))
     {
         m_count++;
-        if (m_current_neighbors[m_count - 1].distance > m_r_max)
+        if (m_current_neighbors[m_count - 1].getDistance() > m_r_max)
         {
             m_finished = true;
             return ITERATOR_TERMINATOR;
diff --git a/cpp/locality/AABBQuery.h b/cpp/locality/AABBQuery.h
index b4c9764fd..c35ffe11c 100644
--- a/cpp/locality/AABBQuery.h
+++ b/cpp/locality/AABBQuery.h
@@ -147,7 +147,7 @@ class AABBQueryIterator : public AABBIterator
                       float r_min, float scale, bool exclude_ii)
         : AABBIterator(neighbor_query, query_point, query_point_idx, r_max, r_min, exclude_ii), m_count(0),
           m_num_neighbors(num_neighbors), m_search_extended(false), m_r_cur(r_guess), m_scale(scale),
-          m_all_distances(), m_query_points_below_r_min()
+          m_all_bonds_minimum_distance(), m_query_points_below_r_min()
     {
         updateImageVectors(0);
     }
@@ -167,8 +167,9 @@ class AABBQueryIterator : public AABBIterator
     float
         m_r_cur; //!< Current search ball cutoff distance in use for the current particle (expands as needed).
     float m_scale; //!< The amount to scale m_r by when the current ball is too small.
-    std::map<unsigned int, float> m_all_distances; //!< Hash map of minimum distances found for a given point,
-                                                   //!< used when searching beyond maximum safe AABB distance.
+    std::map<unsigned int, NeighborBond>
+        m_all_bonds_minimum_distance; //!< Hash map of minimum distances found for a given point,
+                                      //!< used when searching beyond maximum safe AABB distance.
     std::unordered_set<unsigned int> m_query_points_below_r_min; //!< The set of query_points that were too
                                                                  //!< close based on the r_min threshold.
 };
diff --git a/cpp/locality/LinkCell.cc b/cpp/locality/LinkCell.cc
index adef1ad23..446a3761b 100644
--- a/cpp/locality/LinkCell.cc
+++ b/cpp/locality/LinkCell.cc
@@ -516,7 +516,7 @@ NeighborBond LinkCellQueryBallIterator::next()
 
             if (r_sq < r_max_sq && r_sq >= r_min_sq)
             {
-                return NeighborBond(m_query_point_idx, j, std::sqrt(r_sq));
+                return NeighborBond(m_query_point_idx, j, std::sqrt(r_sq), 1, r_ij);
             }
         }
 
@@ -604,7 +604,7 @@ NeighborBond LinkCellQueryIterator::next()
                     const float r_sq(dot(r_ij, r_ij));
                     if (r_sq < r_max_sq && r_sq >= r_min_sq)
                     {
-                        m_current_neighbors.emplace_back(m_query_point_idx, j, std::sqrt(r_sq));
+                        m_current_neighbors.emplace_back(m_query_point_idx, j, std::sqrt(r_sq), 1, r_ij);
                     }
                 }
             }
@@ -639,7 +639,7 @@ NeighborBond LinkCellQueryIterator::next()
             // closest possible neighbor in the new shell.
             std::sort(m_current_neighbors.begin(), m_current_neighbors.end());
             if ((m_current_neighbors.size() >= m_num_neighbors)
-                && (m_current_neighbors[m_num_neighbors - 1].distance
+                && (m_current_neighbors[m_num_neighbors - 1].getDistance()
                     < static_cast<float>(m_neigh_cell_iter.getRange() - 1) * m_linkcell->getCellWidth()))
             {
                 break;
@@ -650,7 +650,7 @@ NeighborBond LinkCellQueryIterator::next()
     while ((m_count < m_num_neighbors) && (m_count < m_current_neighbors.size()))
     {
         m_count++;
-        if (m_current_neighbors[m_count - 1].distance > m_r_max)
+        if (m_current_neighbors[m_count - 1].getDistance() > m_r_max)
         {
             m_finished = true;
             return ITERATOR_TERMINATOR;
diff --git a/cpp/locality/NeighborBond.h b/cpp/locality/NeighborBond.h
index 93037c195..26b74fa3b 100644
--- a/cpp/locality/NeighborBond.h
+++ b/cpp/locality/NeighborBond.h
@@ -1,6 +1,11 @@
+// Copyright (c) 2010-2020 The Regents of the University of Michigan
+// This file is from the freud project, released under the BSD 3-Clause License.
+
 #ifndef NEIGHBOR_BOND_H
 #define NEIGHBOR_BOND_H
 
+#include "VectorMath.h"
+
 namespace freud { namespace locality {
 
 //! Simple data structure encoding neighboring points.
@@ -9,20 +14,38 @@ namespace freud { namespace locality {
  *  class defines the less than operator according to distance, making it
  *  possible to sort.
  */
-struct NeighborBond
+class NeighborBond
 {
+public:
     // For now, id = query_point_idx and ref_id = point_idx (into the NeighborQuery).
     constexpr NeighborBond() = default;
 
-    constexpr NeighborBond(unsigned int query_point_idx, unsigned int point_idx, float d = 0, float w = 1)
-        : query_point_idx(query_point_idx), point_idx(point_idx), distance(d), weight(w)
+    /*
+     * This constructor should only be used in cases where the correct distance/vector
+     * pair is known before instantiating this object to save time and eliminate
+     * redundancy. Common cases include NeighbQuery objects which have to compute
+     * the distance to know if a neighbor pair fits the query arguments and NeighborList
+     * objects which have already computed the distances corresponding to their vectors.
+     * */
+    constexpr NeighborBond(unsigned int query_point_idx, unsigned int point_idx, float d, float w,
+                           const vec3<float>& v)
+        : query_point_idx(query_point_idx), point_idx(point_idx), distance(d), weight(w), vector(v)
+    {}
+
+    /**
+     * This constructor should be used for the majority of cases when instantiating
+     * a NeighborBond object.
+     */
+    NeighborBond(unsigned int query_point_idx, unsigned int point_idx, float w, const vec3<float>& v)
+        : query_point_idx(query_point_idx), point_idx(point_idx), distance(std::sqrt(dot(v, v))), weight(w),
+          vector(v)
     {}
 
     //! Equality checks both query_point_idx and distance.
     bool operator==(const NeighborBond& other) const
     {
         return (query_point_idx == other.query_point_idx) && (point_idx == other.point_idx)
-            && (distance == other.distance);
+            && (distance == other.distance) && (vector == other.vector);
     }
 
     //! Not equals checks inverse of equality.
@@ -87,10 +110,58 @@ struct NeighborBond
         return weight < n.weight;
     }
 
+    unsigned int getQueryPointIdx() const
+    {
+        return query_point_idx;
+    }
+
+    void setQueryPointIdx(unsigned int idx)
+    {
+        query_point_idx = idx;
+    }
+
+    unsigned int getPointIdx() const
+    {
+        return point_idx;
+    }
+
+    void setPointIdx(unsigned int new_pidx)
+    {
+        point_idx = new_pidx;
+    }
+
+    float getWeight() const
+    {
+        return weight;
+    }
+
+    void setWeight(float w)
+    {
+        weight = w;
+    }
+
+    const vec3<float>& getVector() const
+    {
+        return vector;
+    }
+
+    void setVector(const vec3<float>& v)
+    {
+        vector = v;
+        distance = std::sqrt(dot(v, v));
+    }
+
+    float getDistance() const
+    {
+        return distance;
+    }
+
+private:
     unsigned int query_point_idx {0}; //! The query point index.
     unsigned int point_idx {0};       //! The reference point index.
     float distance {0};               //! The distance between the points.
     float weight {0};                 //! The weight of this bond.
+    vec3<float> vector;               //! The directed vector from query point to reference point.
 };
 
 }; }; // end namespace freud::locality
diff --git a/cpp/locality/NeighborComputeFunctional.h b/cpp/locality/NeighborComputeFunctional.h
index 97e8103f3..e768994e6 100644
--- a/cpp/locality/NeighborComputeFunctional.h
+++ b/cpp/locality/NeighborComputeFunctional.h
@@ -24,17 +24,6 @@ NeighborList makeDefaultNlist(const NeighborQuery* nq, const NeighborList* nlist
                               const vec3<float>* query_points, unsigned int num_query_points,
                               locality::QueryArgs qargs);
 
-//! Compute the vector corresponding to a NeighborBond.
-/*! The primary purpose of this function is to standardize the directionality
- * of the delta vector, which is defined as pointing from the query_point to
- * the point (point - query_point), wrapped into the box.
- */
-inline vec3<float> bondVector(const NeighborBond& nb, const NeighborQuery* nq,
-                              const vec3<float>* query_points)
-{
-    return nq->getBox().wrap((*nq)[nb.point_idx] - query_points[nb.query_point_idx]);
-}
-
 //! Implementation of per-point finding logic for NeighborList objects.
 /*! This class provides a concrete implementation of the per-point neighbor
  *  finding interface specified by the NeighborPerPointIterator. In particular,
@@ -68,9 +57,10 @@ class NeighborListPerPointIterator : public NeighborPerPointIterator
 
         NeighborBond nb = NeighborBond(
             m_nlist->getNeighbors()(m_current_index, 0), m_nlist->getNeighbors()(m_current_index, 1),
-            m_nlist->getDistances()[m_current_index], m_nlist->getWeights()[m_current_index]);
+            m_nlist->getDistances()[m_current_index], m_nlist->getWeights()[m_current_index],
+            m_nlist->getVectors()[m_current_index]);
         ++m_current_index;
-        m_returned_point_index = nb.query_point_idx;
+        m_returned_point_index = nb.getQueryPointIdx();
         return nb;
     }
 
@@ -187,7 +177,8 @@ void loopOverNeighbors(const NeighborQuery* neighbor_query, const vec3<float>* q
                 for (size_t bond = begin; bond != end; ++bond)
                 {
                     const NeighborBond nb(nlist->getNeighbors()(bond, 0), nlist->getNeighbors()(bond, 1),
-                                          nlist->getDistances()[bond], nlist->getWeights()[bond]);
+                                          nlist->getDistances()[bond], nlist->getWeights()[bond],
+                                          nlist->getVectors()[bond]);
                     cf(nb);
                 }
             },
diff --git a/cpp/locality/NeighborList.cc b/cpp/locality/NeighborList.cc
index 1101b83fd..f7cc424ba 100644
--- a/cpp/locality/NeighborList.cc
+++ b/cpp/locality/NeighborList.cc
@@ -4,18 +4,19 @@
 #include <algorithm>
 #include <stdexcept>
 
+#include "ManagedArray.h"
 #include "NeighborList.h"
 
 namespace freud { namespace locality {
 
 NeighborList::NeighborList()
-    : m_num_query_points(0), m_num_points(0), m_neighbors({0, 2}), m_distances(0), m_weights(0),
+    : m_num_query_points(0), m_num_points(0), m_neighbors({0, 2}), m_distances(0), m_weights(0), m_vectors(0),
       m_segments_counts_updated(false)
 {}
 
 NeighborList::NeighborList(unsigned int num_bonds)
     : m_num_query_points(0), m_num_points(0), m_neighbors({num_bonds, 2}), m_distances(num_bonds),
-      m_weights(num_bonds), m_segments_counts_updated(false)
+      m_weights(num_bonds), m_vectors(num_bonds), m_segments_counts_updated(false)
 {}
 
 NeighborList::NeighborList(const NeighborList& other)
@@ -27,9 +28,9 @@ NeighborList::NeighborList(const NeighborList& other)
 
 NeighborList::NeighborList(unsigned int num_bonds, const unsigned int* query_point_index,
                            unsigned int num_query_points, const unsigned int* point_index,
-                           unsigned int num_points, const float* distances, const float* weights)
+                           unsigned int num_points, const vec3<float>* vectors, const float* weights)
     : m_num_query_points(num_query_points), m_num_points(num_points), m_neighbors({num_bonds, 2}),
-      m_distances(num_bonds), m_weights(num_bonds), m_segments_counts_updated(false)
+      m_distances(num_bonds), m_weights(num_bonds), m_vectors(num_bonds), m_segments_counts_updated(false)
 {
     unsigned int last_index(0);
     for (unsigned int i = 0; i < num_bonds; i++)
@@ -48,10 +49,7 @@ NeighborList::NeighborList(unsigned int num_bonds, const unsigned int* query_poi
         {
             throw std::invalid_argument("NeighborList point_index values must be less than num_points.");
         }
-        m_neighbors(i, 0) = index;
-        m_neighbors(i, 1) = point_index[i];
-        m_weights[i] = weights[i];
-        m_distances[i] = distances[i];
+        setNeighborEntry(i, NeighborBond(index, point_index[i], weights[i], vectors[i]));
         last_index = index;
     }
 }
@@ -132,6 +130,7 @@ template<typename Iterator> unsigned int NeighborList::filter(Iterator begin)
     auto new_neighbors = util::ManagedArray<unsigned int>({new_size, 2});
     auto new_distances = util::ManagedArray<float>(new_size);
     auto new_weights = util::ManagedArray<float>(new_size);
+    auto new_vectors = util::ManagedArray<vec3<float>>(new_size);
 
     auto current_element = begin;
     unsigned int num_good(0);
@@ -141,8 +140,9 @@ template<typename Iterator> unsigned int NeighborList::filter(Iterator begin)
         {
             new_neighbors(num_good, 0) = m_neighbors(i, 0);
             new_neighbors(num_good, 1) = m_neighbors(i, 1);
-            new_weights[num_good] = m_weights[i];
             new_distances[num_good] = m_distances[i];
+            new_weights[num_good] = m_weights[i];
+            new_vectors[num_good] = m_vectors[i];
             ++num_good;
         }
         ++current_element;
@@ -151,6 +151,7 @@ template<typename Iterator> unsigned int NeighborList::filter(Iterator begin)
     m_neighbors = new_neighbors;
     m_distances = new_distances;
     m_weights = new_weights;
+    m_vectors = new_vectors;
     m_segments_counts_updated = false;
     return old_size - new_size;
 }
@@ -199,6 +200,7 @@ void NeighborList::resize(unsigned int num_bonds)
     auto new_neighbors = util::ManagedArray<unsigned int>({num_bonds, 2});
     auto new_distances = util::ManagedArray<float>(num_bonds);
     auto new_weights = util::ManagedArray<float>(num_bonds);
+    auto new_vectors = util::ManagedArray<vec3<float>>(num_bonds);
 
     // On shrinking resizes, keep existing data.
     if (num_bonds <= getNumBonds())
@@ -209,12 +211,14 @@ void NeighborList::resize(unsigned int num_bonds)
             new_neighbors(i, 1) = m_neighbors(i, 1);
             new_distances[i] = m_distances[i];
             new_weights[i] = m_weights[i];
+            new_vectors[i] = m_vectors[i];
         }
     }
 
     m_neighbors = new_neighbors;
     m_distances = new_distances;
     m_weights = new_weights;
+    m_vectors = new_vectors;
     m_segments_counts_updated = false;
 }
 
@@ -222,8 +226,9 @@ void NeighborList::copy(const NeighborList& other)
 {
     setNumBonds(other.getNumBonds(), other.getNumQueryPoints(), other.getNumPoints());
     m_neighbors = other.m_neighbors.copy();
-    m_weights = other.m_weights.copy();
     m_distances = other.m_distances.copy();
+    m_weights = other.m_weights.copy();
+    m_vectors = other.m_vectors.copy();
     m_segments_counts_updated = false;
 }
 
diff --git a/cpp/locality/NeighborList.h b/cpp/locality/NeighborList.h
index 17dde7b2c..f0789d235 100644
--- a/cpp/locality/NeighborList.h
+++ b/cpp/locality/NeighborList.h
@@ -35,7 +35,7 @@ class NeighborList
     NeighborList(const NeighborList& other);
     //! Construct from arrays
     NeighborList(unsigned int num_bonds, const unsigned int* query_point_index, unsigned int num_query_points,
-                 const unsigned int* point_index, unsigned int num_points, const float* distances,
+                 const unsigned int* point_index, unsigned int num_points, const vec3<float>* vectors,
                  const float* weights);
 
     //! Return the number of bonds stored in this NeighborList
@@ -50,49 +50,30 @@ class NeighborList
     //! Update the arrays of neighbor counts and segments
     void updateSegmentCounts() const;
 
-    //! Access the neighbors array for reading and writing
-    util::ManagedArray<unsigned int>& getNeighbors()
-    {
-        return m_neighbors;
-    }
-    //! Access the distances array for reading and writing
-    util::ManagedArray<float>& getDistances()
-    {
-        return m_distances;
-    }
-    //! Access the weights array for reading and writing
-    util::ManagedArray<float>& getWeights()
-    {
-        return m_weights;
-    }
-    //! Access the counts array for reading
-    util::ManagedArray<unsigned int>& getCounts()
-    {
-        updateSegmentCounts();
-        return m_counts;
-    }
-    //! Access the segments array for reading
-    util::ManagedArray<unsigned int>& getSegments()
-    {
-        updateSegmentCounts();
-        return m_segments;
-    }
-
     //! Access the neighbors array for reading
     const util::ManagedArray<unsigned int>& getNeighbors() const
     {
         return m_neighbors;
     }
+
     //! Access the distances array for reading
     const util::ManagedArray<float>& getDistances() const
     {
         return m_distances;
     }
+
     //! Access the weights array for reading
     const util::ManagedArray<float>& getWeights() const
     {
         return m_weights;
     }
+
+    //! Access the vectors array for reading
+    const util::ManagedArray<vec3<float>>& getVectors() const
+    {
+        return m_vectors;
+    }
+
     //! Access the counts array for reading
     const util::ManagedArray<unsigned int>& getCounts() const
     {
@@ -106,6 +87,18 @@ class NeighborList
         return m_segments;
     }
 
+    /**
+     * Set the values for the neighbor index to be that of the given neighborbond
+     */
+    void setNeighborEntry(size_t neighbor_index, const NeighborBond& nb)
+    {
+        m_neighbors(neighbor_index, 0) = nb.getQueryPointIdx();
+        m_neighbors(neighbor_index, 1) = nb.getPointIdx();
+        m_vectors[neighbor_index] = nb.getVector();
+        m_distances[neighbor_index] = nb.getDistance();
+        m_weights[neighbor_index] = nb.getWeight();
+    }
+
     //! Remove bonds in this object based on an array of boolean values. The
     //  array must be at least as long as the number of neighbor bonds.
     //  Returns the number of bonds removed.
@@ -140,6 +133,8 @@ class NeighborList
     util::ManagedArray<float> m_distances;
     //! Neighbor list per-bond weight array
     util::ManagedArray<float> m_weights;
+    //!< Directed vectors per-bond array
+    util::ManagedArray<vec3<float>> m_vectors;
 
     //! Track whether segments and counts are up to date
     mutable bool m_segments_counts_updated;
diff --git a/cpp/locality/NeighborQuery.h b/cpp/locality/NeighborQuery.h
index 0d3c094ba..64e41fc7a 100644
--- a/cpp/locality/NeighborQuery.h
+++ b/cpp/locality/NeighborQuery.h
@@ -39,7 +39,7 @@ constexpr float DEFAULT_R_GUESS(-1.0);                    //!< Default guess que
 constexpr float DEFAULT_SCALE(-1.0);      //!< Default scaling parameter for AABB nearest neighbor queries.
 constexpr bool DEFAULT_EXCLUDE_II(false); //!< Default for whether or not to include self-neighbors.
 constexpr auto ITERATOR_TERMINATOR
-    = NeighborBond(-1, -1, 0); //!< The object returned when iteration is complete.
+    = NeighborBond(-1, -1, 0, 0, vec3<float>()); //!< The object returned when iteration is complete.
 
 //! POD class to hold information about generic queries.
 /*! This class provides a standard method for specifying the type of query to
@@ -395,7 +395,8 @@ class NeighborQueryIterator
                     // If we're excluding ii bonds, we have to check before adding.
                     if (nb != ITERATOR_TERMINATOR)
                     {
-                        local_bonds.emplace_back(nb.query_point_idx, nb.point_idx, nb.distance);
+                        local_bonds.emplace_back(nb.getQueryPointIdx(), nb.getPointIdx(), nb.getWeight(),
+                                                 nb.getVector());
                     }
                 }
             }
@@ -420,10 +421,7 @@ class NeighborQueryIterator
         util::forLoopWrapper(0, num_bonds, [&](size_t begin, size_t end) {
             for (size_t bond = begin; bond < end; ++bond)
             {
-                nl->getNeighbors()(bond, 0) = linear_bonds[bond].query_point_idx;
-                nl->getNeighbors()(bond, 1) = linear_bonds[bond].point_idx;
-                nl->getDistances()[bond] = linear_bonds[bond].distance;
-                nl->getWeights()[bond] = float(1.0);
+                nl->setNeighborEntry(bond, linear_bonds[bond]);
             }
         });
 
diff --git a/cpp/locality/Voronoi.cc b/cpp/locality/Voronoi.cc
index 35feb907f..347772bf1 100644
--- a/cpp/locality/Voronoi.cc
+++ b/cpp/locality/Voronoi.cc
@@ -1,6 +1,7 @@
 // Copyright (c) 2010-2020 The Regents of the University of Michigan
 // This file is from the freud project, released under the BSD 3-Clause License.
 
+#include <algorithm>
 #include <cmath>
 #include <iterator>
 #include <tbb/parallel_sort.h>
@@ -124,8 +125,10 @@ void Voronoi::compute(const freud::locality::NeighborQuery* nq)
 
             // Compute cell neighbors
             size_t neighbor_counter(0);
+            size_t face_vertices_index(0);
             for (auto neighbor_iterator = neighbors.begin(); neighbor_iterator != neighbors.end();
-                 neighbor_iterator++, neighbor_counter++)
+                 ++neighbor_iterator, ++neighbor_counter,
+                      face_vertices_index += face_vertices[face_vertices_index] + 1)
             {
                 // Get the normal to the current face
                 const vec3<double> normal(normals[3 * neighbor_counter], normals[3 * neighbor_counter + 1],
@@ -148,13 +151,41 @@ void Voronoi::compute(const freud::locality::NeighborQuery* nq)
                 // Fetch neighbor information
                 const int point_id = *neighbor_iterator;
                 const float weight(face_areas[neighbor_counter]);
-                const vec3<double> point_system_coords((*nq)[point_id]);
 
-                // Compute the distance from query_point to point.
-                const vec3<float> rij = box.wrap(point_system_coords - query_point_system_coords);
-                const float distance(std::sqrt(dot(rij, rij)));
-
-                bonds.emplace_back(query_point_id, point_id, distance, weight);
+                // Find a vertex on the current face: this leverages the
+                // structure of face_vertices, which has a count of the
+                // number of vertices for a face followed by the
+                // corresponding vertex ids for that face. We use this
+                // structure later when incrementing face_vertices_index.
+                // face_vertices_index always points to the "vertex
+                // counter" element of face_vertices for the current face.
+
+                // Get the id of the vertex on this face that is most parallel to the normal
+                const auto normal_length = std::sqrt(dot(normal, normal));
+                auto cosine_vertex_to_normal = [&](const auto& vertex_id_on_face) {
+                    const vec3<double> rv(vertices[3 * vertex_id_on_face],
+                                          vertices[3 * vertex_id_on_face + 1],
+                                          vertices[3 * vertex_id_on_face + 2]);
+                    const vec3<double> riv(rv - query_point);
+                    return dot(riv, normal) / std::sqrt(dot(riv, riv)) / normal_length;
+                };
+
+                const int vertex_id_on_face = *std::max_element(
+                    &(face_vertices[face_vertices_index + 1]),
+                    &(face_vertices[face_vertices_index + face_vertices[face_vertices_index]]),
+                    [&](const auto& a, const auto& b) {
+                        return cosine_vertex_to_normal(a) < cosine_vertex_to_normal(b);
+                    });
+
+                // Project the vertex vector onto the face normal to get a
+                // vector from query_point to the face, then double it to
+                // get the vector to the neighbor particle.
+                const vec3<double> rv(vertices[3 * vertex_id_on_face], vertices[3 * vertex_id_on_face + 1],
+                                      vertices[3 * vertex_id_on_face + 2]);
+                const vec3<double> riv(rv - query_point);
+                const vec3<float> vector(2.0 * dot(riv, normal) * normal);
+
+                bonds.emplace_back(query_point_id, point_id, weight, vector);
             }
 
         } while (voronoi_loop.inc());
@@ -172,10 +203,7 @@ void Voronoi::compute(const freud::locality::NeighborQuery* nq)
     util::forLoopWrapper(0, num_bonds, [&](size_t begin, size_t end) {
         for (size_t bond = begin; bond != end; ++bond)
         {
-            m_neighbor_list->getNeighbors()(bond, 0) = bonds[bond].query_point_idx;
-            m_neighbor_list->getNeighbors()(bond, 1) = bonds[bond].point_idx;
-            m_neighbor_list->getDistances()[bond] = bonds[bond].distance;
-            m_neighbor_list->getWeights()[bond] = bonds[bond].weight;
+            m_neighbor_list->setNeighborEntry(bond, bonds[bond]);
         }
     });
 }
diff --git a/cpp/order/HexaticTranslational.cc b/cpp/order/HexaticTranslational.cc
index 8f706158c..ba477a921 100644
--- a/cpp/order/HexaticTranslational.cc
+++ b/cpp/order/HexaticTranslational.cc
@@ -28,8 +28,8 @@ void HexaticTranslational<T>::computeGeneral(Func func, const freud::locality::N
             for (freud::locality::NeighborBond nb = ppiter->next(); !ppiter->end(); nb = ppiter->next())
             {
                 // Compute vector from query_point to point
-                const vec3<float> delta = box.wrap((*points)[nb.point_idx] - ref);
-                const float weight(m_weighted ? nb.weight : 1.0);
+                const vec3<float> delta = box.wrap((*points)[nb.getPointIdx()] - ref);
+                const float weight(m_weighted ? nb.getWeight() : 1.0);
 
                 // Compute psi for this vector
                 m_psi_array[i] += weight * func(delta);
diff --git a/cpp/order/Steinhardt.cc b/cpp/order/Steinhardt.cc
index c403d433f..5891d527e 100644
--- a/cpp/order/Steinhardt.cc
+++ b/cpp/order/Steinhardt.cc
@@ -130,8 +130,8 @@ void Steinhardt::baseCompute(const freud::locality::NeighborList* nlist,
 
             for (freud::locality::NeighborBond nb = ppiter->next(); !ppiter->end(); nb = ppiter->next())
             {
-                const vec3<float> delta = points->getBox().wrap((*points)[nb.point_idx] - ref);
-                const float weight(m_weighted ? nb.weight : float(1.0));
+                const vec3<float> delta = points->getBox().wrap((*points)[nb.getPointIdx()] - ref);
+                const float weight(m_weighted ? nb.getWeight() : float(1.0));
 
                 // phi is usually in range 0..2Pi, but
                 // it only appears in Ylm as exp(im\phi),
@@ -142,11 +142,11 @@ void Steinhardt::baseCompute(const freud::locality::NeighborList* nlist,
                 // aligned along z, otherwise due to floating point error we
                 // could get delta.z/nb.distance = -1-eps, which is outside the
                 // valid range of std::acos.
-                float theta = std::acos(util::clamp(delta.z / nb.distance, -1, 1)); // 0..Pi
+                float theta = std::acos(util::clamp(delta.z / nb.getDistance(), -1, 1)); // 0..Pi
 
                 // If the points are directly on top of each other,
                 // theta should be zero instead of nan.
-                if (nb.distance == float(0))
+                if (nb.getDistance() == float(0))
                 {
                     theta = 0;
                 }
@@ -225,7 +225,7 @@ void Steinhardt::computeAve(const freud::locality::NeighborList* nlist,
                     auto& qlmiAve = m_qlmiAve[l_index];
                     auto& qlmi = m_qlmi[l_index];
                     const auto ave_index = qlmiAve.getIndex({i, 0});
-                    const auto nb_index = qlmi.getIndex({nb.point_idx, 0});
+                    const auto nb_index = qlmi.getIndex({nb.getPointIdx(), 0});
                     for (size_t k = 0; k < m_num_ms[l_index]; ++k)
                     {
                         // Adding all the qlm of the neighbors. We use the
diff --git a/cpp/pmft/PMFTR12.cc b/cpp/pmft/PMFTR12.cc
index 3a2cbd55e..119e43b6c 100644
--- a/cpp/pmft/PMFTR12.cc
+++ b/cpp/pmft/PMFTR12.cc
@@ -82,16 +82,17 @@ void PMFTR12::accumulate(const locality::NeighborQuery* neighbor_query, const fl
     neighbor_query->getBox().enforce2D();
     accumulateGeneral(neighbor_query, query_points, n_query_points, nlist, qargs,
                       [&](const freud::locality::NeighborBond& neighbor_bond) {
-                          vec3<float> delta(bondVector(neighbor_bond, neighbor_query, query_points));
+                          const vec3<float>& delta(neighbor_bond.getVector());
                           // calculate angles
-                          float d_theta1 = std::atan2(delta.y, delta.x);
-                          float d_theta2 = std::atan2(-delta.y, -delta.x);
-                          float t1 = orientations[neighbor_bond.point_idx] - d_theta1;
-                          float t2 = query_orientations[neighbor_bond.query_point_idx] - d_theta2;
+                          const float d_theta1 = std::atan2(delta.y, delta.x);
+                          const float d_theta2 = std::atan2(-delta.y, -delta.x);
                           // make sure that t1, t2 are bounded between 0 and 2PI
-                          t1 = util::modulusPositive(t1, constants::TWO_PI);
-                          t2 = util::modulusPositive(t2, constants::TWO_PI);
-                          m_local_histograms(neighbor_bond.distance, t1, t2);
+                          const float t1 = util::modulusPositive(
+                              orientations[neighbor_bond.getPointIdx()] - d_theta1, constants::TWO_PI);
+                          const float t2 = util::modulusPositive(
+                              query_orientations[neighbor_bond.getQueryPointIdx()] - d_theta2,
+                              constants::TWO_PI);
+                          m_local_histograms(neighbor_bond.getDistance(), t1, t2);
                       });
 }
 
diff --git a/cpp/pmft/PMFTXY.cc b/cpp/pmft/PMFTXY.cc
index 93ce56dd5..1a3c0157f 100644
--- a/cpp/pmft/PMFTXY.cc
+++ b/cpp/pmft/PMFTXY.cc
@@ -62,13 +62,13 @@ void PMFTXY::accumulate(const locality::NeighborQuery* neighbor_query, const flo
     neighbor_query->getBox().enforce2D();
     accumulateGeneral(neighbor_query, query_points, n_query_points, nlist, qargs,
                       [&](const freud::locality::NeighborBond& neighbor_bond) {
-                          vec3<float> delta(bondVector(neighbor_bond, neighbor_query, query_points));
+                          const vec3<float>& delta(neighbor_bond.getVector());
 
                           // rotate interparticle vector
-                          vec2<float> myVec(delta.x, delta.y);
-                          rotmat2<float> myMat
-                              = rotmat2<float>::fromAngle(-query_orientations[neighbor_bond.query_point_idx]);
-                          vec2<float> rotVec = myMat * myVec;
+                          const vec2<float> myVec(delta.x, delta.y);
+                          const rotmat2<float> myMat(rotmat2<float>::fromAngle(
+                              -query_orientations[neighbor_bond.getQueryPointIdx()]));
+                          const vec2<float> rotVec = myMat * myVec;
 
                           m_local_histograms(rotVec.x, rotVec.y);
                       });
diff --git a/cpp/pmft/PMFTXYT.cc b/cpp/pmft/PMFTXYT.cc
index 35c8eb528..93f5a67ec 100644
--- a/cpp/pmft/PMFTXYT.cc
+++ b/cpp/pmft/PMFTXYT.cc
@@ -70,18 +70,18 @@ void PMFTXYT::accumulate(const locality::NeighborQuery* neighbor_query, const fl
     neighbor_query->getBox().enforce2D();
     accumulateGeneral(neighbor_query, query_points, n_query_points, nlist, qargs,
                       [&](const freud::locality::NeighborBond& neighbor_bond) {
-                          vec3<float> delta(bondVector(neighbor_bond, neighbor_query, query_points));
+                          const vec3<float>& delta(neighbor_bond.getVector());
 
                           // rotate interparticle vector
-                          vec2<float> myVec(delta.x, delta.y);
-                          rotmat2<float> myMat
-                              = rotmat2<float>::fromAngle(-query_orientations[neighbor_bond.query_point_idx]);
-                          vec2<float> rotVec = myMat * myVec;
+                          const vec2<float> myVec(delta.x, delta.y);
+                          const rotmat2<float> myMat(rotmat2<float>::fromAngle(
+                              -query_orientations[neighbor_bond.getQueryPointIdx()]));
+                          const vec2<float> rotVec = myMat * myVec;
                           // calculate angle
-                          float d_theta = std::atan2(-delta.y, -delta.x);
-                          float t = orientations[neighbor_bond.point_idx] - d_theta;
+                          const float d_theta = std::atan2(-delta.y, -delta.x);
                           // make sure that t is bounded between 0 and 2PI
-                          t = util::modulusPositive(t, constants::TWO_PI);
+                          const float t = util::modulusPositive(
+                              orientations[neighbor_bond.getPointIdx()] - d_theta, constants::TWO_PI);
                           m_local_histograms(rotVec.x, rotVec.y, t);
                       });
 }
diff --git a/cpp/pmft/PMFTXYZ.cc b/cpp/pmft/PMFTXYZ.cc
index 5efd975da..6701e3392 100644
--- a/cpp/pmft/PMFTXYZ.cc
+++ b/cpp/pmft/PMFTXYZ.cc
@@ -116,9 +116,10 @@ void PMFTXYZ::accumulate(const locality::NeighborQuery* neighbor_query, const qu
     accumulateGeneral(neighbor_query, query_points, n_query_points, nlist, qargs,
                       [&](const freud::locality::NeighborBond& neighbor_bond) {
                           // create the reference point quaternion
-                          quat<float> query_orientation(query_orientations[neighbor_bond.query_point_idx]);
+                          const quat<float> query_orientation(
+                              query_orientations[neighbor_bond.getQueryPointIdx()]);
                           // make sure that the particles are wrapped into the box
-                          vec3<float> delta(bondVector(neighbor_bond, neighbor_query, query_points));
+                          const vec3<float>& delta(neighbor_bond.getVector());
 
                           for (unsigned int k = 0; k < num_equiv_orientations; k++)
                           {
diff --git a/cpp/util/VectorMath.h b/cpp/util/VectorMath.h
index 9b1d53f39..2bd1e3832 100644
--- a/cpp/util/VectorMath.h
+++ b/cpp/util/VectorMath.h
@@ -32,16 +32,16 @@ template<class Real> struct vec3
         \param _y y-component
         \param _z z-component
     */
-    vec3(const Real& _x, const Real& _y, const Real& _z) : x(_x), y(_y), z(_z) {}
+    constexpr vec3(const Real& _x, const Real& _y, const Real& _z) : x(_x), y(_y), z(_z) {}
 
     //! Implicit cast from vec3<double> to the current Real
-    vec3(const vec3<double>& a) : x(a.x), y(a.y), z(a.z) {}
+    constexpr vec3(const vec3<double>& a) : x(a.x), y(a.y), z(a.z) {}
 
     //! Implicit cast from vec3<float> to the current Real
-    vec3(const vec3<float>& a) : x(a.x), y(a.y), z(a.z) {}
+    constexpr vec3(const vec3<float>& a) : x(a.x), y(a.y), z(a.z) {}
 
     //! Default construct a 0 vector
-    vec3() = default;
+    constexpr vec3() = default;
 
     //! Swap with another vector
     void swap(vec3<Real>& v)
@@ -63,7 +63,7 @@ template<class Real> struct vec3
     Addition is component wise.
     \returns The vector (a.x+b.x, a.y+b.y, a.z+b.z).
 */
-template<class Real> inline vec3<Real> operator+(const vec3<Real>& a, const vec3<Real>& b)
+template<class Real> constexpr vec3<Real> operator+(const vec3<Real>& a, const vec3<Real>& b)
 {
     return vec3<Real>(a.x + b.x, a.y + b.y, a.z + b.z);
 }
@@ -75,7 +75,7 @@ template<class Real> inline vec3<Real> operator+(const vec3<Real>& a, const vec3
     Subtraction is component wise.
     \returns The vector (a.x-b.x, a.y-b.y, a.z-b.z).
 */
-template<class Real> inline vec3<Real> operator-(const vec3<Real>& a, const vec3<Real>& b)
+template<class Real> constexpr vec3<Real> operator-(const vec3<Real>& a, const vec3<Real>& b)
 {
     return vec3<Real>(a.x - b.x, a.y - b.y, a.z - b.z);
 }
@@ -87,7 +87,7 @@ template<class Real> inline vec3<Real> operator-(const vec3<Real>& a, const vec3
     Multiplication is component wise.
     \returns The vector (a.x*b.x, a.y*b.y, a.z*b.z).
 */
-template<class Real> inline vec3<Real> operator*(const vec3<Real>& a, const vec3<Real>& b)
+template<class Real> constexpr vec3<Real> operator*(const vec3<Real>& a, const vec3<Real>& b)
 {
     return vec3<Real>(a.x * b.x, a.y * b.y, a.z * b.z);
 }
@@ -99,7 +99,7 @@ template<class Real> inline vec3<Real> operator*(const vec3<Real>& a, const vec3
     Division is component wise.
     \returns The vector (a.x/b.x, a.y/b.y, a.z/b.z).
 */
-template<class Real> inline vec3<Real> operator/(const vec3<Real>& a, const vec3<Real>& b)
+template<class Real> constexpr vec3<Real> operator/(const vec3<Real>& a, const vec3<Real>& b)
 {
     return vec3<Real>(a.x / b.x, a.y / b.y, a.z / b.z);
 }
@@ -110,7 +110,7 @@ template<class Real> inline vec3<Real> operator/(const vec3<Real>& a, const vec3
     Negation is component wise.
     \returns The vector (-a.x, -a.y, -a.z).
 */
-template<class Real> inline vec3<Real> operator-(const vec3<Real>& a)
+template<class Real> constexpr vec3<Real> operator-(const vec3<Real>& a)
 {
     return vec3<Real>(-a.x, -a.y, -a.z);
 }
@@ -122,7 +122,7 @@ template<class Real> inline vec3<Real> operator-(const vec3<Real>& a)
     Addition is component wise.
     \returns The vector (a.x += b.x, a.y += b.y, a.z += b.z).
 */
-template<class Real> inline vec3<Real>& operator+=(vec3<Real>& a, const vec3<Real>& b)
+template<class Real> constexpr vec3<Real>& operator+=(vec3<Real>& a, const vec3<Real>& b)
 {
     a.x += b.x;
     a.y += b.y;
@@ -137,7 +137,7 @@ template<class Real> inline vec3<Real>& operator+=(vec3<Real>& a, const vec3<Rea
     Subtraction is component wise.
     \returns The vector (a.x -= b.x, a.y -= b.y, a.z -= b.z).
 */
-template<class Real> inline vec3<Real>& operator-=(vec3<Real>& a, const vec3<Real>& b)
+template<class Real> constexpr vec3<Real>& operator-=(vec3<Real>& a, const vec3<Real>& b)
 {
     a.x -= b.x;
     a.y -= b.y;
@@ -152,7 +152,7 @@ template<class Real> inline vec3<Real>& operator-=(vec3<Real>& a, const vec3<Rea
     Multiplication is component wise.
     \returns The vector (a.x *= b.x, a.y *= b.y, a.z *= b.z).
 */
-template<class Real> inline vec3<Real>& operator*=(vec3<Real>& a, const vec3<Real>& b)
+template<class Real> constexpr vec3<Real>& operator*=(vec3<Real>& a, const vec3<Real>& b)
 {
     a.x *= b.x;
     a.y *= b.y;
@@ -167,7 +167,7 @@ template<class Real> inline vec3<Real>& operator*=(vec3<Real>& a, const vec3<Rea
     Division is component wise.
     \returns The vector (a.x /= b.x, a.y /= b.y, a.z /= b.z).
 */
-template<class Real> inline vec3<Real>& operator/=(vec3<Real>& a, const vec3<Real>& b)
+template<class Real> constexpr vec3<Real>& operator/=(vec3<Real>& a, const vec3<Real>& b)
 {
     a.x /= b.x;
     a.y /= b.y;
@@ -182,7 +182,7 @@ template<class Real> inline vec3<Real>& operator/=(vec3<Real>& a, const vec3<Rea
     Multiplication is component wise.
     \returns The vector (a.x*b, a.y*b, a.z*b).
 */
-template<class Real> inline vec3<Real> operator*(const vec3<Real>& a, const Real& b)
+template<class Real> constexpr vec3<Real> operator*(const vec3<Real>& a, const Real& b)
 {
     return vec3<Real>(a.x * b, a.y * b, a.z * b);
 }
@@ -194,7 +194,7 @@ template<class Real> inline vec3<Real> operator*(const vec3<Real>& a, const Real
     Multiplication is component wise.
     \returns The vector (a.x*b, a.y*b, a.z*b).
 */
-template<class Real> inline vec3<Real> operator*(const Real& b, const vec3<Real>& a)
+template<class Real> constexpr vec3<Real> operator*(const Real& b, const vec3<Real>& a)
 {
     return vec3<Real>(a.x * b, a.y * b, a.z * b);
 }
@@ -206,7 +206,7 @@ template<class Real> inline vec3<Real> operator*(const Real& b, const vec3<Real>
     Division is component wise.
     \returns The vector (a.x/b, a.y/b, a.z/b).
 */
-template<class Real> inline vec3<Real> operator/(const vec3<Real>& a, const Real& b)
+template<class Real> constexpr vec3<Real> operator/(const vec3<Real>& a, const Real& b)
 {
     Real q = Real(1.0) / b;
     return a * q;
@@ -219,7 +219,7 @@ template<class Real> inline vec3<Real> operator/(const vec3<Real>& a, const Real
     Multiplication is component wise.
     \returns The vector (a.x *= b, a.y *= b, a.z *= b).
 */
-template<class Real> inline vec3<Real>& operator*=(vec3<Real>& a, const Real& b)
+template<class Real> constexpr vec3<Real>& operator*=(vec3<Real>& a, const Real& b)
 {
     a.x *= b;
     a.y *= b;
@@ -234,7 +234,7 @@ template<class Real> inline vec3<Real>& operator*=(vec3<Real>& a, const Real& b)
     Division is component wise.
     \returns The vector (a.x /= b, a.y /= b, a.z /= b).
 */
-template<class Real> inline vec3<Real>& operator/=(vec3<Real>& a, const Real& b)
+template<class Real> constexpr vec3<Real>& operator/=(vec3<Real>& a, const Real& b)
 {
     a.x /= b;
     a.y /= b;
@@ -247,7 +247,7 @@ template<class Real> inline vec3<Real>& operator/=(vec3<Real>& a, const Real& b)
     \param b Second vector
     \returns true if the two vectors are identically equal, false if they are not
 */
-template<class Real> inline bool operator==(const vec3<Real>& a, const vec3<Real>& b)
+template<class Real> constexpr bool operator==(const vec3<Real>& a, const vec3<Real>& b)
 {
     return (a.x == b.x) && (a.y == b.y) && (a.z == b.z);
 }
@@ -257,7 +257,7 @@ template<class Real> inline bool operator==(const vec3<Real>& a, const vec3<Real
     \param b Second vector
     \returns true if the two vectors are not identically equal, and false if they are
 */
-template<class Real> inline bool operator!=(const vec3<Real>& a, const vec3<Real>& b)
+template<class Real> constexpr bool operator!=(const vec3<Real>& a, const vec3<Real>& b)
 {
     return (a.x != b.x) || (a.y != b.y) || (a.z != b.z);
 }
@@ -268,7 +268,7 @@ template<class Real> inline bool operator!=(const vec3<Real>& a, const vec3<Real
 
     \returns the dot product a.x*b.x + a.y*b.y + a.z*b.z.
 */
-template<class Real> inline Real dot(const vec3<Real>& a, const vec3<Real>& b)
+template<class Real> constexpr Real dot(const vec3<Real>& a, const vec3<Real>& b)
 {
     return (a.x * b.x + a.y * b.y + a.z * b.z);
 }
@@ -279,7 +279,7 @@ template<class Real> inline Real dot(const vec3<Real>& a, const vec3<Real>& b)
 
     \returns the cross product (a.y * b.z - a.z * b.y, a.z * b.x - a.x * b.z, a.x * b.y - a.y * b.x).
 */
-template<class Real> inline vec3<Real> cross(const vec3<Real>& a, const vec3<Real>& b)
+template<class Real> constexpr vec3<Real> cross(const vec3<Real>& a, const vec3<Real>& b)
 {
     return vec3<Real>(a.y * b.z - a.z * b.y, a.z * b.x - a.x * b.z, a.x * b.y - a.y * b.x);
 }
@@ -303,16 +303,16 @@ template<class Real> struct vec2
     /*! \param _x x-component
         \param _y y-component
     */
-    vec2(const Real& _x, const Real& _y) : x(_x), y(_y) {}
+    constexpr vec2(const Real& _x, const Real& _y) : x(_x), y(_y) {}
 
     //! Default construct a 0 vector
-    vec2() : x(0), y(0) {}
+    constexpr vec2() : x(0), y(0) {}
 
     //! Implicit cast from vec2<double> to the current Real
-    vec2(const vec2<double>& a) : x(a.x), y(a.y) {}
+    constexpr vec2(const vec2<double>& a) : x(a.x), y(a.y) {}
 
     //! Implicit cast from vec2<float> to the current Real
-    vec2(const vec2<float>& a) : x(a.x), y(a.y) {}
+    constexpr vec2(const vec2<float>& a) : x(a.x), y(a.y) {}
 
     //! Swap with another vector
     void swap(vec2<Real>& v)
@@ -338,7 +338,7 @@ template<class Real> struct vec2
     Addition is component wise.
     \returns The vector (a.x+b.x, a.y+b.y).
 */
-template<class Real> inline vec2<Real> operator+(const vec2<Real>& a, const vec2<Real>& b)
+template<class Real> constexpr vec2<Real> operator+(const vec2<Real>& a, const vec2<Real>& b)
 {
     return vec2<Real>(a.x + b.x, a.y + b.y);
 }
@@ -350,7 +350,7 @@ template<class Real> inline vec2<Real> operator+(const vec2<Real>& a, const vec2
     Subtraction is component wise.
     \returns The vector (a.x-b.x, a.y-b.y).
 */
-template<class Real> inline vec2<Real> operator-(const vec2<Real>& a, const vec2<Real>& b)
+template<class Real> constexpr vec2<Real> operator-(const vec2<Real>& a, const vec2<Real>& b)
 {
     return vec2<Real>(a.x - b.x, a.y - b.y);
 }
@@ -362,7 +362,7 @@ template<class Real> inline vec2<Real> operator-(const vec2<Real>& a, const vec2
     Multiplication is component wise.
     \returns The vector (a.x*b.x, a.y*b.y).
 */
-template<class Real> inline vec2<Real> operator*(const vec2<Real>& a, const vec2<Real>& b)
+template<class Real> constexpr vec2<Real> operator*(const vec2<Real>& a, const vec2<Real>& b)
 {
     return vec2<Real>(a.x * b.x, a.y * b.y);
 }
@@ -374,7 +374,7 @@ template<class Real> inline vec2<Real> operator*(const vec2<Real>& a, const vec2
     Division is component wise.
     \returns The vector (a.x/b.x, a.y/b.y).
 */
-template<class Real> inline vec2<Real> operator/(const vec2<Real>& a, const vec2<Real>& b)
+template<class Real> constexpr vec2<Real> operator/(const vec2<Real>& a, const vec2<Real>& b)
 {
     return vec2<Real>(a.x / b.x, a.y / b.y);
 }
@@ -385,7 +385,7 @@ template<class Real> inline vec2<Real> operator/(const vec2<Real>& a, const vec2
     Negation is component wise.
     \returns The vector (-a.x, -a.y).
 */
-template<class Real> inline vec2<Real> operator-(const vec2<Real>& a)
+template<class Real> constexpr vec2<Real> operator-(const vec2<Real>& a)
 {
     return vec2<Real>(-a.x, -a.y);
 }
@@ -397,7 +397,7 @@ template<class Real> inline vec2<Real> operator-(const vec2<Real>& a)
     Addition is component wise.
     \returns The vector (a.x += b.x, a.y += b.y).
 */
-template<class Real> inline vec2<Real>& operator+=(vec2<Real>& a, const vec2<Real>& b)
+template<class Real> constexpr vec2<Real>& operator+=(vec2<Real>& a, const vec2<Real>& b)
 {
     a.x += b.x;
     a.y += b.y;
@@ -411,7 +411,7 @@ template<class Real> inline vec2<Real>& operator+=(vec2<Real>& a, const vec2<Rea
     Subtraction is component wise.
     \returns The vector (a.x -= b.x, a.y -= b.y).
 */
-template<class Real> inline vec2<Real>& operator-=(vec2<Real>& a, const vec2<Real>& b)
+template<class Real> constexpr vec2<Real>& operator-=(vec2<Real>& a, const vec2<Real>& b)
 {
     a.x -= b.x;
     a.y -= b.y;
@@ -425,7 +425,7 @@ template<class Real> inline vec2<Real>& operator-=(vec2<Real>& a, const vec2<Rea
     Multiplication is component wise.
     \returns The vector (a.x *= b.x, a.y *= b.y).
 */
-template<class Real> inline vec2<Real>& operator*=(vec2<Real>& a, const vec2<Real>& b)
+template<class Real> constexpr vec2<Real>& operator*=(vec2<Real>& a, const vec2<Real>& b)
 {
     a.x *= b.x;
     a.y *= b.y;
@@ -439,7 +439,7 @@ template<class Real> inline vec2<Real>& operator*=(vec2<Real>& a, const vec2<Rea
     Division is component wise.
     \returns The vector (a.x /= b.x, a.y /= b.y).
 */
-template<class Real> inline vec2<Real>& operator/=(vec2<Real>& a, const vec2<Real>& b)
+template<class Real> constexpr vec2<Real>& operator/=(vec2<Real>& a, const vec2<Real>& b)
 {
     a.x /= b.x;
     a.y /= b.y;
@@ -453,7 +453,7 @@ template<class Real> inline vec2<Real>& operator/=(vec2<Real>& a, const vec2<Rea
     Multiplication is component wise.
     \returns The vector (a.x*b, a.y*b).
 */
-template<class Real> inline vec2<Real> operator*(const vec2<Real>& a, const Real& b)
+template<class Real> constexpr vec2<Real> operator*(const vec2<Real>& a, const Real& b)
 {
     return vec2<Real>(a.x * b, a.y * b);
 }
@@ -465,7 +465,7 @@ template<class Real> inline vec2<Real> operator*(const vec2<Real>& a, const Real
     Multiplication is component wise.
     \returns The vector (a.x*b, a.y*b).
 */
-template<class Real> inline vec2<Real> operator*(const Real& b, const vec2<Real>& a)
+template<class Real> constexpr vec2<Real> operator*(const Real& b, const vec2<Real>& a)
 {
     return vec2<Real>(a.x * b, a.y * b);
 }
@@ -477,7 +477,7 @@ template<class Real> inline vec2<Real> operator*(const Real& b, const vec2<Real>
     Division is component wise.
     \returns The vector (a.x/b, a.y/b, a.z/b).
 */
-template<class Real> inline vec2<Real> operator/(const vec2<Real>& a, const Real& b)
+template<class Real> constexpr vec2<Real> operator/(const vec2<Real>& a, const Real& b)
 {
     Real q = Real(1.0) / b;
     return a * q;
@@ -490,7 +490,7 @@ template<class Real> inline vec2<Real> operator/(const vec2<Real>& a, const Real
     Multiplication is component wise.
     \returns The vector (a.x *= b, a.y *= b).
 */
-template<class Real> inline vec2<Real>& operator*=(vec2<Real>& a, const Real& b)
+template<class Real> constexpr vec2<Real>& operator*=(vec2<Real>& a, const Real& b)
 {
     a.x *= b;
     a.y *= b;
@@ -504,7 +504,7 @@ template<class Real> inline vec2<Real>& operator*=(vec2<Real>& a, const Real& b)
     Division is component wise.
     \returns The vector (a.x /= b, a.y /= b).
 */
-template<class Real> inline vec2<Real>& operator/=(vec2<Real>& a, const Real& b)
+template<class Real> constexpr vec2<Real>& operator/=(vec2<Real>& a, const Real& b)
 {
     a.x /= b;
     a.y /= b;
@@ -516,7 +516,7 @@ template<class Real> inline vec2<Real>& operator/=(vec2<Real>& a, const Real& b)
     \param b Second vector
     \returns true if the two vectors are identically equal, false if they are not
 */
-template<class Real> inline bool operator==(const vec2<Real>& a, const vec2<Real>& b)
+template<class Real> constexpr bool operator==(const vec2<Real>& a, const vec2<Real>& b)
 {
     return (a.x == b.x) && (a.y == b.y);
 }
@@ -526,7 +526,7 @@ template<class Real> inline bool operator==(const vec2<Real>& a, const vec2<Real
     \param b Second vector
     \returns true if the two vectors are not identically equal, false if they are
 */
-template<class Real> inline bool operator!=(const vec2<Real>& a, const vec2<Real>& b)
+template<class Real> constexpr bool operator!=(const vec2<Real>& a, const vec2<Real>& b)
 {
     return (a.x != b.x) || (a.y != b.y);
 }
@@ -537,7 +537,7 @@ template<class Real> inline bool operator!=(const vec2<Real>& a, const vec2<Real
 
     \returns the dot product a.x*b.x + a.y*b.y.
 */
-template<class Real> inline Real dot(const vec2<Real>& a, const vec2<Real>& b)
+template<class Real> constexpr Real dot(const vec2<Real>& a, const vec2<Real>& b)
 {
     return (a.x * b.x + a.y * b.y);
 }
@@ -546,7 +546,7 @@ template<class Real> inline Real dot(const vec2<Real>& a, const vec2<Real>& b)
 /*! \param a vector
     \returns a vector perpendicular to *a* (a.y, -a.x)
 */
-template<class Real> inline vec2<Real> perp(const vec2<Real>& a)
+template<class Real> constexpr vec2<Real> perp(const vec2<Real>& a)
 {
     return vec2<Real>(-a.y, a.x);
 }
@@ -556,7 +556,7 @@ template<class Real> inline vec2<Real> perp(const vec2<Real>& a)
     \param b second vector
     \returns the perpendicular dot product of a and b
 */
-template<class Real> inline Real perpdot(const vec2<Real>& a, const vec2<Real>& b)
+template<class Real> constexpr Real perpdot(const vec2<Real>& a, const vec2<Real>& b)
 {
     return dot(perp(a), b);
 }
@@ -598,19 +598,19 @@ template<class Real> struct quat
     /*! \param _s scalar component
         \param _v vector component
     */
-    quat(const Real& _s, const vec3<Real>& _v) : s(_s), v(_v) {}
+    constexpr quat(const Real& _s, const vec3<Real>& _v) : s(_s), v(_v) {}
 
     //! Implicit cast from quat<double> to the current Real
-    quat(const quat<double>& a) : s(a.s), v(a.v) {}
+    constexpr quat(const quat<double>& a) : s(a.s), v(a.v) {}
 
     //! Implicit cast from quat<float> to the current Real
-    quat(const quat<float>& a) : s(a.s), v(a.v) {}
+    constexpr quat(const quat<float>& a) : s(a.s), v(a.v) {}
 
     //! Default construct a unit quaternion
-    quat() : s(1), v(vec3<Real>(0, 0, 0)) {}
+    constexpr quat() : s(1), v(vec3<Real>(0, 0, 0)) {}
 
     //! Construct a quaternion from a rotation matrix
-    explicit quat(const rotmat3<Real>& r);
+    constexpr explicit quat(const rotmat3<Real>& r);
 
     //! Construct a quat from an axis and an angle.
     /*! \param axis angle to represent
@@ -636,7 +636,7 @@ template<class Real> struct quat
     Multiplication is component wise.
     \returns The quaternion (b*a.s, b*a.v).
 */
-template<class Real> inline quat<Real> operator*(const Real& b, const quat<Real>& a)
+template<class Real> constexpr quat<Real> operator*(const Real& b, const quat<Real>& a)
 {
     return quat<Real>(b * a.s, b * a.v);
 }
@@ -648,7 +648,7 @@ template<class Real> inline quat<Real> operator*(const Real& b, const quat<Real>
     Multiplication is component wise.
     \returns The quaternion (a.s*b, a.v*b).
 */
-template<class Real> inline quat<Real> operator*(const quat<Real>& a, const Real& b)
+template<class Real> constexpr quat<Real> operator*(const quat<Real>& a, const Real& b)
 {
     return quat<Real>(a.s * b, a.v * b);
 }
@@ -660,7 +660,7 @@ template<class Real> inline quat<Real> operator*(const quat<Real>& a, const Real
     Addition is component wise.
     \returns The quaternion (a.s + b.s, a.v+b.v).
 */
-template<class Real> inline quat<Real> operator+(const quat<Real>& a, const quat<Real>& b)
+template<class Real> constexpr quat<Real> operator+(const quat<Real>& a, const quat<Real>& b)
 {
     return quat<Real>(a.s + b.s, a.v + b.v);
 }
@@ -672,7 +672,7 @@ template<class Real> inline quat<Real> operator+(const quat<Real>& a, const quat
     Addition is component wise.
     \returns The quaternion (a.s += b.s, a.v += b.v).
 */
-template<class Real> inline quat<Real>& operator+=(quat<Real>& a, const quat<Real>& b)
+template<class Real> constexpr quat<Real>& operator+=(quat<Real>& a, const quat<Real>& b)
 {
     a.s += b.s;
     a.v += b.v;
@@ -686,7 +686,7 @@ template<class Real> inline quat<Real>& operator+=(quat<Real>& a, const quat<Rea
     Subtraction is component wise.
     \returns The quaternion (a.s - b.s, a.v - b.v).
 */
-template<class Real> inline quat<Real> operator-(const quat<Real>& a, const quat<Real>& b)
+template<class Real> constexpr quat<Real> operator-(const quat<Real>& a, const quat<Real>& b)
 {
     return quat<Real>(a.s - b.s, a.v - b.v);
 }
@@ -698,7 +698,7 @@ template<class Real> inline quat<Real> operator-(const quat<Real>& a, const quat
     Subtraction is component wise.
     \returns The quaternion (a.s -= b.s, a.v -= b.v).
 */
-template<class Real> inline quat<Real>& operator-=(quat<Real>& a, const quat<Real>& b)
+template<class Real> constexpr quat<Real>& operator-=(quat<Real>& a, const quat<Real>& b)
 {
     a.s -= b.s;
     a.v -= b.v;
@@ -717,7 +717,7 @@ template<class Real> inline quat<Real>& operator-=(quat<Real>& a, const quat<Rea
     \returns The quaternion
     (a.s * b.s − dot(a.v, b.v), a.s*b.v + b.s * a.v + cross(a.v, b.v)).
 */
-template<class Real> inline quat<Real> operator*(const quat<Real>& a, const quat<Real>& b)
+template<class Real> constexpr quat<Real> operator*(const quat<Real>& a, const quat<Real>& b)
 {
     return quat<Real>(a.s * b.s - dot(a.v, b.v), a.s * b.v + b.s * a.v + cross(a.v, b.v));
 }
@@ -730,7 +730,7 @@ template<class Real> inline quat<Real> operator*(const quat<Real>& a, const quat
 
     \returns The quaternion (a.s * b.s − dot(a.v, b.v), a.s*b.v + b.s * a.v + cross(a.v, b.v)).
 */
-template<class Real> inline quat<Real> operator*(const vec3<Real>& a, const quat<Real>& b)
+template<class Real> constexpr quat<Real> operator*(const vec3<Real>& a, const quat<Real>& b)
 {
     return quat<Real>(0, a) * b;
 }
@@ -743,7 +743,7 @@ template<class Real> inline quat<Real> operator*(const vec3<Real>& a, const quat
 
     \returns The quaternion (a.s * b.s − dot(a.v, b.v), a.s*b.v + b.s * a.v + cross(a.v, b.v)).
 */
-template<class Real> inline quat<Real> operator*(const quat<Real>& a, const vec3<Real>& b)
+template<class Real> constexpr quat<Real> operator*(const quat<Real>& a, const vec3<Real>& b)
 {
     return a * quat<Real>(0, b);
 }
@@ -753,7 +753,7 @@ template<class Real> inline quat<Real> operator*(const quat<Real>& a, const vec3
 
     \returns the norm of the quaternion, squared. (a.s*a.s + dot(a.v,a.v))
 */
-template<class Real> inline Real norm2(const quat<Real>& a)
+template<class Real> constexpr Real norm2(const quat<Real>& a)
 {
     return (a.s * a.s + dot(a.v, a.v));
 }
@@ -763,7 +763,7 @@ template<class Real> inline Real norm2(const quat<Real>& a)
 
     \returns the conjugate of the quaternion. (a.s, -a.v)
 */
-template<class Real> inline quat<Real> conj(const quat<Real>& a)
+template<class Real> constexpr quat<Real> conj(const quat<Real>& a)
 {
     return quat<Real>(a.s, -a.v);
 }
@@ -772,7 +772,7 @@ template<class Real> inline quat<Real> conj(const quat<Real>& a)
 /*! \note The rotation matrix must have positive determinant
  http://www.euclideanspace.com/maths/geometry/rotations/conversions/matrixToQuaternion/
  */
-template<class Real> inline quat<Real>::quat(const rotmat3<Real>& r)
+template<class Real> constexpr quat<Real>::quat(const rotmat3<Real>& r)
 {
     Real tr = r.row0.x + r.row1.y + r.row2.z;
 
@@ -808,7 +808,7 @@ template<class Real> inline quat<Real>::quat(const rotmat3<Real>& r)
 
     \returns the vector rotated by the quaternion, equivalent to the vector component of a*b*conj(a);
 */
-template<class Real> inline vec3<Real> rotate(const quat<Real>& a, const vec3<Real>& b)
+template<class Real> constexpr vec3<Real> rotate(const quat<Real>& a, const vec3<Real>& b)
 {
     // quat<Real> result = a*b*conj(a);
     // return result.v;
@@ -831,7 +831,7 @@ template<class Real> inline vec3<Real> rotate(const quat<Real>& a, const vec3<Re
 
     *b* is promoted to a 3d vector with z=0 for the rotation.
 */
-template<class Real> inline vec2<Real> rotate(const quat<Real>& a, const vec2<Real>& b)
+template<class Real> constexpr vec2<Real> rotate(const quat<Real>& a, const vec2<Real>& b)
 {
     vec3<Real> b3(b.x, b.y, Real(0.0));
     // b3 = (a*b3*conj(a)).v;
@@ -853,7 +853,7 @@ template<class Real> inline vec2<Real> rotate(const quat<Real>& a, const vec2<Re
 
     \returns the dot product a.s*b.s+a.v.x*b.v.x + a.v.y*b.v.y + a.v.z*b.v.z.
 */
-template<class Real> inline Real dot(const quat<Real>& a, const quat<Real>& b)
+template<class Real> constexpr Real dot(const quat<Real>& a, const quat<Real>& b)
 {
     return (a.s * b.s + dot(a.v, b.v));
 }
@@ -932,7 +932,7 @@ template<class Real> struct rotmat2
 
     Multiplication is matrix multiplication, where the vector is represented as a column vector.
 */
-template<class Real> inline vec2<Real> operator*(const rotmat2<Real>& A, const vec2<Real>& b)
+template<class Real> constexpr vec2<Real> operator*(const rotmat2<Real>& A, const vec2<Real>& b)
 {
     return vec2<Real>(dot(A.row0, b), dot(A.row1, b));
 }
@@ -944,7 +944,7 @@ template<class Real> inline vec2<Real> operator*(const rotmat2<Real>& A, const v
     A rotation matrix has an inverse equal to its transpose. There may be times where an algorithm needs to
    undo a rotation, so the transpose method is provided.
 */
-template<class Real> inline rotmat2<Real> transpose(const rotmat2<Real>& A)
+template<class Real> constexpr rotmat2<Real> transpose(const rotmat2<Real>& A)
 {
     return rotmat2<Real>(vec2<Real>(A.row0.x, A.row1.x), vec2<Real>(A.row0.y, A.row1.y));
 }
@@ -1036,7 +1036,7 @@ template<class Real> struct rotmat3
 
     Multiplication is matrix multiplication, where the vector is represented as a column vector.
 */
-template<class Real> inline vec3<Real> operator*(const rotmat3<Real>& A, const vec3<Real>& b)
+template<class Real> constexpr vec3<Real> operator*(const rotmat3<Real>& A, const vec3<Real>& b)
 {
     return vec3<Real>(dot(A.row0, b), dot(A.row1, b), dot(A.row2, b));
 }
@@ -1046,7 +1046,7 @@ template<class Real> inline vec3<Real> operator*(const rotmat3<Real>& A, const v
     \param B matrix
     \returns A*b
 */
-template<class Real> inline rotmat3<Real> operator*(const rotmat3<Real>& A, const rotmat3<Real>& B)
+template<class Real> constexpr rotmat3<Real> operator*(const rotmat3<Real>& A, const rotmat3<Real>& B)
 {
     rotmat3<Real> r;
     rotmat3<Real> B_t = transpose(B);
@@ -1069,7 +1069,7 @@ template<class Real> inline rotmat3<Real> operator*(const rotmat3<Real>& A, cons
     A rotation matrix has an inverse equal to its transpose. There may be times where an algorithm needs to
    undo a rotation, so the transpose method is provided.
 */
-template<class Real> inline rotmat3<Real> transpose(const rotmat3<Real>& A)
+template<class Real> constexpr rotmat3<Real> transpose(const rotmat3<Real>& A)
 {
     return rotmat3<Real>(vec3<Real>(A.row0.x, A.row1.x, A.row2.x), vec3<Real>(A.row0.y, A.row1.y, A.row2.y),
                          vec3<Real>(A.row0.z, A.row1.z, A.row2.z));
@@ -1083,7 +1083,7 @@ template<class Real> inline rotmat3<Real> transpose(const rotmat3<Real>& A)
     \returns the projection of a onto b
     \note projection() can be applied to 2d or 3d vectors
 */
-template<class Vec> inline Vec project(const Vec& a, const Vec& b)
+template<class Vec> constexpr Vec project(const Vec& a, const Vec& b)
 {
     return dot(a, b) / dot(b, b) * b;
 }
diff --git a/doc/source/gettingstarted/migration.rst b/doc/source/gettingstarted/migration.rst
index 44808781b..3ab1d9997 100644
--- a/doc/source/gettingstarted/migration.rst
+++ b/doc/source/gettingstarted/migration.rst
@@ -26,3 +26,9 @@ Overview of API Changes
     * - Color voronoi diagram by number of cell sides.
       - ``voro.plot(..., color_by_sides=True)``
       - ``voro.plot(..., color_by='sides')``
+    * - Get vectors corresponding to neighbor pairs.
+      - ``points[nlist.point_indices] - points[nlist.query_point_indices]``
+      - ``nlist.vectors``
+    * - Create a custom neighborlist from arrays.
+      - ``freud.locality.NeighborList.from_arrays(..., distances=...)``
+      - ``freud.locality.NeighborList.from_arrays(..., vectors=...)``
diff --git a/doc/source/reference/credits.rst b/doc/source/reference/credits.rst
index 4deaefe3c..890668642 100644
--- a/doc/source/reference/credits.rst
+++ b/doc/source/reference/credits.rst
@@ -124,6 +124,7 @@ Bradley Dice - **Lead developer**
 * Fixed ``Box.contains`` to run in linear time, ``O(num_points)``.
 * Contributed code, design, documentation, and testing for ``StaticStructureFactorDirect`` class.
 * Fixed doctests to run with pytest.
+* Work with Tommy Waltmann on adding neighbor vectors to ``freud.locality.NeighborList``.
 
 Eric Harper, University of Michigan - **Former lead developer**
 
@@ -334,6 +335,7 @@ Tommy Waltmann
 * Add support for compilation with the C++17 standard.
 * Update and test the ``normalization_mode`` argument to ``freud.density.RDF`` class.
 * Extending plotting options for the ``Voronoi`` module.
+* Work with Bradley Dice on adding neighbor vectors to ``freud.locality.NeighborList``.
 
 Maya Martirossyan
 
diff --git a/freud/_locality.pxd b/freud/_locality.pxd
index f28828fd6..ae654754d 100644
--- a/freud/_locality.pxd
+++ b/freud/_locality.pxd
@@ -13,10 +13,11 @@ from freud.util cimport vec3
 
 cdef extern from "NeighborBond.h" namespace "freud::locality":
     cdef cppclass NeighborBond:
-        unsigned int query_point_idx
-        unsigned int point_idx
-        float distance
-        float weight
+        unsigned int getQueryPointIdx() const
+        unsigned int getPointIdx() const
+        float getDistance() const
+        float getWeight() const
+        vec3[float] getVector() const
         bool operator==(const NeighborBond &) const
         bool operator!=(const NeighborBond &) const
         bool operator<(const NeighborBond &) const
@@ -72,12 +73,13 @@ cdef extern from "NeighborList.h" namespace "freud::locality":
         NeighborList()
         NeighborList(unsigned int)
         NeighborList(unsigned int, const unsigned int*, unsigned int,
-                     const unsigned int*, unsigned int, const float*,
+                     const unsigned int*, unsigned int, const vec3[float]*,
                      const float*) except +
 
         freud.util.ManagedArray[unsigned int] &getNeighbors()
         freud.util.ManagedArray[float] &getDistances()
         freud.util.ManagedArray[float] &getWeights()
+        freud.util.ManagedArray[vec3[float]] &getVectors()
         freud.util.ManagedArray[float] &getSegments()
         freud.util.ManagedArray[float] &getCounts()
 
diff --git a/freud/locality.pyx b/freud/locality.pyx
index c8109d161..7ce3bf2a5 100644
--- a/freud/locality.pyx
+++ b/freud/locality.pyx
@@ -198,7 +198,8 @@ cdef class NeighborQueryResult:
 
         npoint = dereference(iterator).next()
         while npoint != ITERATOR_TERMINATOR:
-            yield (npoint.query_point_idx, npoint.point_idx, npoint.distance)
+            yield (npoint.getQueryPointIdx(), npoint.getPointIdx(),
+                   npoint.getDistance())
             npoint = dereference(iterator).next()
 
         raise StopIteration
@@ -473,8 +474,7 @@ cdef class NeighborList:
        nlist = aq.query(positions, {'r_max': 3}).toNeighborList()
 
        # Get all vectors from central particles to their neighbors
-       rijs = (positions[nlist.point_indices] -
-              positions[nlist.query_point_indices])
+       rijs = nlist.vectors
        rijs = box.wrap(rijs)
 
     The NeighborList can be indexed to access bond particle indices. Example::
@@ -485,7 +485,7 @@ cdef class NeighborList:
 
     @classmethod
     def from_arrays(cls, num_query_points, num_points, query_point_indices,
-                    point_indices, distances, weights=None):
+                    point_indices, vectors, weights=None):
         r"""Create a NeighborList from a set of bond information arrays.
 
         Example::
@@ -495,15 +495,14 @@ cdef class NeighborList:
             box = freud.box.Box(2, 3, 4, 0, 0, 0)
             query_points = np.array([[0, 0, 0], [0, 0, 1]])
             points = np.array([[0, 0, -1], [0.5, -1, 0]])
-            query_point_indices = np.array([0, 0, 1])
-            point_indices = np.array([0, 1, 1])
-            distances = box.compute_distances(
-                query_points[query_point_indices], points[point_indices])
             num_query_points = len(query_points)
             num_points = len(points)
+            query_point_indices = np.array([0, 0, 1])
+            point_indices = np.array([0, 1, 1])
+            vectors = box.wrap(points[point_indices] - query_points[query_point_indices])
             nlist = freud.locality.NeighborList.from_arrays(
                 num_query_points, num_points, query_point_indices,
-                point_indices, distances)
+                point_indices, vectors)
 
 
         Args:
@@ -518,10 +517,10 @@ cdef class NeighborList:
             point_indices (:class:`np.ndarray`):
                 Array of integers corresponding to indices in the set of
                 points.
-            distances (:class:`np.ndarray`):
-                Array of distances between corresponding query points and
+            vectors (:math:`\left(N_{bonds}, 3\right)` :class:`numpy.ndarray`):
+                Array of bond vectors from query points to corresponding
                 points.
-            weights (:class:`np.ndarray`, optional):
+            weights (:math:`\left(N_{bonds} \right)` :class:`np.ndarray`, optional):
                 Array of per-bond weights (if :code:`None` is given, use a
                 value of 1 for each weight) (Default value = :code:`None`).
         """  # noqa 501
@@ -530,8 +529,8 @@ cdef class NeighborList:
         point_indices = freud.util._convert_array(
             point_indices, shape=query_point_indices.shape, dtype=np.uint32)
 
-        distances = freud.util._convert_array(
-            distances, shape=query_point_indices.shape)
+        vectors = freud.util._convert_array(
+            vectors, shape=(len(query_point_indices), 3))
 
         if weights is None:
             weights = np.ones(query_point_indices.shape, dtype=np.float32)
@@ -541,7 +540,7 @@ cdef class NeighborList:
         cdef const unsigned int[::1] l_query_point_indices = \
             query_point_indices
         cdef const unsigned int[::1] l_point_indices = point_indices
-        cdef const float[::1] l_distances = distances
+        cdef const float[:, ::1] l_vectors = vectors
         cdef const float[::1] l_weights = weights
         cdef unsigned int l_num_bonds = l_query_point_indices.shape[0]
         cdef unsigned int l_num_query_points = num_query_points
@@ -551,7 +550,8 @@ cdef class NeighborList:
         result = cls()
         result.thisptr = new freud._locality.NeighborList(
             l_num_bonds, &l_query_point_indices[0], l_num_query_points,
-            &l_point_indices[0], l_num_points, &l_distances[0], &l_weights[0])
+            &l_point_indices[0], l_num_points, <vec3[float]*> &l_vectors[0, 0],
+            &l_weights[0])
 
         return result
 
@@ -632,6 +632,14 @@ cdef class NeighborList:
             &self.thisptr.getDistances(),
             freud.util.arr_type_t.FLOAT)
 
+    @property
+    def vectors(self):
+        r"""(:math:`N_{bonds}`, 3) :class:`np.ndarray`: The vectors for each
+        bond."""
+        return freud.util.make_managed_numpy_array(
+            &self.thisptr.getVectors(),
+            freud.util.arr_type_t.FLOAT, 3)
+
     @property
     def segments(self):
         """(:math:`N_{query\\_points}`) :class:`np.ndarray`: A segment array
diff --git a/setup.cfg b/setup.cfg
index b94e92fce..0af5b0a6b 100644
--- a/setup.cfg
+++ b/setup.cfg
@@ -10,7 +10,7 @@ filename = *.py,*.pyx,*.pxi,*.pxd
 force-check = True
 exclude = .eggs,*.egg,build,extern,doc/source/gettingstarted/examples
 select = E,F,W
-ignore = E203,E225,E226,E227,E741,E999,W503,W504
+ignore = E203,E225,E226,E227,E402,E741,E999,W503,W504
 per-file-ignores = 
 	freud/__init__.py: F401
 	freud/*.pxd: E402
diff --git a/tests/test_environment_BondOrder.py b/tests/test_environment_BondOrder.py
index b15bd0a4b..a55531745 100644
--- a/tests/test_environment_BondOrder.py
+++ b/tests/test_environment_BondOrder.py
@@ -113,7 +113,7 @@ def test_points_ne_query_points(self):
         r_max = 1.6
 
         num_neighbors = 12
-        n_bins_theta = 30
+        n_bins_theta = 35
         n_bins_phi = 2
         test_set = util.make_raw_query_nlist_test_set(
             box, points, query_points, "nearest", r_max, num_neighbors, False
diff --git a/tests/test_environment_LocalBondProjection.py b/tests/test_environment_LocalBondProjection.py
index 2f58a3f53..ff92fea19 100644
--- a/tests/test_environment_LocalBondProjection.py
+++ b/tests/test_environment_LocalBondProjection.py
@@ -86,10 +86,10 @@ def test_compute(self):
         ang = freud.environment.LocalBondProjection()
         ang.compute((box, points), ors, proj_vecs, neighbors=query_args)
 
-        dnlist = freud.locality.AABBQuery(box, points).query(
+        nlist = freud.locality.AABBQuery(box, points).query(
             points, dict(num_neighbors=num_neighbors, r_guess=r_guess, exclude_ii=True)
         )
-        bonds = [(i[0], i[1]) for i in dnlist]
+        bonds = [(i[0], i[1]) for i in nlist]
 
         # We will look at the bond between [1, 0, 0] as query_point
         # and [0, 0, 0] as point
diff --git a/tests/test_environment_LocalDescriptors.py b/tests/test_environment_LocalDescriptors.py
index 9930d3194..beb484b08 100644
--- a/tests/test_environment_LocalDescriptors.py
+++ b/tests/test_environment_LocalDescriptors.py
@@ -239,11 +239,12 @@ def test_ql(self, unit_cell):
             # in cases where there is no symmetry. Since simple cubic
             # should have a 0 ql value in many cases, we need to set high
             # tolerances for those specific cases.
-            atol = 1e-3 if unit_cell == freud.data.UnitCell.sc else 1e-6
+            atol = 1e-3 if unit_cell == freud.data.UnitCell.sc else 1e-5
             npt.assert_allclose(
                 steinhardt.particle_order,
                 ql[:, L],
                 atol=atol,
+                rtol=1e-5,
                 err_msg=f"Failed for {unit_cell.__name__}, L = {L}",
             )
 
@@ -276,7 +277,7 @@ def test_ql_weighted(self, unit_cell):
             len(points),
             nl.query_point_indices,
             nl.point_indices,
-            nl.distances,
+            nl.vectors,
             np.random.rand(len(nl.weights)),
         )
 
@@ -368,7 +369,7 @@ def test_ld(self):
                     scipy_val = sph_harm(m, l, phi, theta)
                     ld_val = (-1) ** abs(m) * ld.sph[idx, count]
                     assert np.isclose(scipy_val, ld_val, atol=atol), (
-                        "Failed for l={}, m={}, x={}, y = {}" "\ntheta={}, phi={}"
+                        "Failed for l={}, m={}, x={}, y={}, theta={}, phi={}"
                     ).format(l, m, scipy_val, ld_val, theta, phi)
                     count += 1
 
@@ -377,7 +378,7 @@ def test_ld(self):
                     scipy_val = sph_harm(m, l, phi, theta)
                     ld_val = ld.sph[idx, count]
                     assert np.isclose(scipy_val, ld_val, atol=atol), (
-                        "Failed for l={}, m={}, x={}, y = {}" "\ntheta={}, phi={}"
+                        "Failed for l={}, m={}, x={}, y={}, theta={}, phi={}"
                     ).format(l, m, scipy_val, ld_val, theta, phi)
                     count += 1
 
@@ -399,15 +400,14 @@ def test_query_point_ne_points(self):
         # again later anyway.
         lc = freud.locality.AABBQuery(box, points)
         nl = lc.query(
-            points, dict(exclude_ii=True, num_neighbors=num_neighbors)
+            query_points, dict(exclude_ii=False, num_neighbors=num_neighbors)
         ).toNeighborList()
 
         ld = freud.environment.LocalDescriptors(l_max, mode="global")
         ld.compute((box, points), query_points, neighbors=nl)
 
         # Loop over the sphs and compute them explicitly.
-        for idx, (i, j) in enumerate(nl):
-            bond = box.wrap(points[j] - query_points[i])
+        for idx, bond in enumerate(nl.vectors):
             r = np.linalg.norm(bond)
             theta = np.arccos(bond[2] / r)
             phi = np.arctan2(bond[1], bond[0])
@@ -423,7 +423,7 @@ def test_query_point_ne_points(self):
                     scipy_val = sph_harm(m, l, phi, theta)
                     ld_val = (-1) ** abs(m) * ld.sph[idx, count]
                     assert np.isclose(scipy_val, ld_val, atol=atol), (
-                        "Failed for l={}, m={}, x={}, y = {}" "\ntheta={}, phi={}"
+                        "Failed for l={}, m={}, x={}, y={}, theta={}, phi={}"
                     ).format(l, m, scipy_val, ld_val, theta, phi)
                     count += 1
 
@@ -432,6 +432,6 @@ def test_query_point_ne_points(self):
                     scipy_val = sph_harm(m, l, phi, theta)
                     ld_val = ld.sph[idx, count]
                     assert np.isclose(scipy_val, ld_val, atol=atol), (
-                        "Failed for l={}, m={}, x={}, y = {}" "\ntheta={}, phi={}"
+                        "Failed for l={}, m={}, x={}, y={}, theta={}, phi={}"
                     ).format(l, m, scipy_val, ld_val, theta, phi)
                     count += 1
diff --git a/tests/test_environment_MatchEnv.py b/tests/test_environment_MatchEnv.py
index c5b6935d8..5b6be31c1 100644
--- a/tests/test_environment_MatchEnv.py
+++ b/tests/test_environment_MatchEnv.py
@@ -490,4 +490,4 @@ def test_api(self):
         query_args = dict(r_guess=r_max, num_neighbors=num_neighbors)
         match.compute((box, points), motif, neighbors=query_args)
         assert np.all(match.rmsds[:-1] > 0)
-        assert match.rmsds[-1] == 0
+        assert np.isclose(match.rmsds[-1], 0, atol=1e-6)
diff --git a/tests/test_locality_NeighborList.py b/tests/test_locality_NeighborList.py
index 3ea2960b3..7d65dd4e8 100644
--- a/tests/test_locality_NeighborList.py
+++ b/tests/test_locality_NeighborList.py
@@ -138,76 +138,80 @@ def test_segments(self):
     def test_from_arrays(self):
         query_point_indices = [0, 0, 1, 2, 3]
         point_indices = [1, 2, 3, 0, 0]
-        distances = np.ones(len(query_point_indices))
+        vectors = np.ones((len(query_point_indices), 3))
 
         # implicit weights
         nlist = freud.locality.NeighborList.from_arrays(
-            4, 4, query_point_indices, point_indices, distances
+            4, 4, query_point_indices, point_indices, vectors
         )
         assert np.allclose(nlist.weights, 1)
 
         # explicit weights
         weights = np.ones(len(query_point_indices)) * 4.0
         nlist = freud.locality.NeighborList.from_arrays(
-            4, 4, query_point_indices, point_indices, distances, weights
+            4, 4, query_point_indices, point_indices, vectors, weights
         )
         assert np.allclose(nlist.weights, 4)
 
         # copy of existing nlist by arrays
         weights = np.random.rand(len(query_point_indices))
         nlist = freud.locality.NeighborList.from_arrays(
-            4, 4, query_point_indices, point_indices, distances, weights
+            4, 4, query_point_indices, point_indices, vectors, weights
         )
         nlist2 = freud.locality.NeighborList.from_arrays(
             4,
             4,
             nlist.query_point_indices,
             nlist.point_indices,
-            nlist.distances,
+            nlist.vectors,
             nlist.weights,
         )
         npt.assert_equal(nlist.query_point_indices, nlist2.query_point_indices)
         npt.assert_equal(nlist.point_indices, nlist2.point_indices)
         npt.assert_equal(nlist.distances, nlist2.distances)
         npt.assert_equal(nlist.weights, nlist2.weights)
+        npt.assert_equal(nlist.vectors, nlist2.vectors)
         npt.assert_equal(nlist.neighbor_counts, nlist2.neighbor_counts)
         npt.assert_equal(nlist.segments, nlist2.segments)
 
         # too few reference particles
         with pytest.raises(ValueError):
             freud.locality.NeighborList.from_arrays(
-                3, 4, query_point_indices, point_indices, distances
+                3, 4, query_point_indices, point_indices, vectors
             )
 
         # too few target particles
         with pytest.raises(ValueError):
             freud.locality.NeighborList.from_arrays(
-                4, 3, query_point_indices, point_indices, distances
+                4, 3, query_point_indices, point_indices, vectors
             )
 
         # query particles not sorted
         with pytest.raises(ValueError):
             freud.locality.NeighborList.from_arrays(
-                4, 4, point_indices, query_point_indices, distances
+                4, 4, point_indices, query_point_indices, vectors
             )
 
         # mismatched array sizes
         with pytest.raises(ValueError):
+            # wrong number of query point indices
             freud.locality.NeighborList.from_arrays(
-                4, 4, query_point_indices[:-1], point_indices, distances
+                4, 4, query_point_indices[:-1], point_indices, vectors
             )
         with pytest.raises(ValueError):
+            # wrong number of point indices
             freud.locality.NeighborList.from_arrays(
-                4, 4, query_point_indices, point_indices[:-1], distances
+                4, 4, query_point_indices, point_indices[:-1], vectors
             )
         with pytest.raises(ValueError):
+            # wrong number of vectors
             freud.locality.NeighborList.from_arrays(
-                4, 4, query_point_indices, point_indices, distances[:-1]
+                4, 4, query_point_indices, point_indices, vectors[:-1]
             )
         with pytest.raises(ValueError):
-            weights = np.ones((len(query_point_indices) - 1,))
+            # wrong number of weights
             freud.locality.NeighborList.from_arrays(
-                4, 4, query_point_indices, point_indices, distances, weights
+                4, 4, query_point_indices, point_indices, vectors, weights[:-1]
             )
 
     def test_indexing_empty(self):
@@ -282,11 +286,11 @@ def test_ordering_distance(self):
     def test_num_points(self):
         query_point_indices = [0, 0, 1, 2, 3]
         point_indices = [1, 2, 3, 0, 0]
-        distances = np.ones(len(query_point_indices))
+        vectors = np.ones((len(query_point_indices), 3))
 
         # test num_query_points and num_points when built from arrays
         nlist = freud.locality.NeighborList.from_arrays(
-            42, 99, query_point_indices, point_indices, distances
+            42, 99, query_point_indices, point_indices, vectors
         )
         assert nlist.num_query_points == 42
         assert nlist.num_points == 99
diff --git a/tests/test_locality_Voronoi.py b/tests/test_locality_Voronoi.py
index 1104cbe83..944dfcaa4 100644
--- a/tests/test_locality_Voronoi.py
+++ b/tests/test_locality_Voronoi.py
@@ -11,6 +11,21 @@
 
 
 class TestVoronoi:
+    def _check_vectors_and_distances(self, vor, sys):
+        """Assert the neighbor vectors/distances have the right length/direction."""
+        box = sys[0]
+        points = sys[1]
+
+        # Verify the neighbor vectors
+        wrapped_points = box.wrap(
+            points[vor.nlist.query_point_indices] + vor.nlist.vectors
+        )
+        npt.assert_allclose(wrapped_points, points[vor.nlist.point_indices], atol=1e-5)
+
+        # Verify the neighbor distances
+        vector_lengths = np.linalg.norm(vor.nlist.vectors, axis=-1)
+        npt.assert_allclose(vector_lengths, vor.nlist.distances)
+
     def test_random_2d(self):
         # Test that voronoi tessellations of random systems have the same
         # number of points and polytopes
@@ -25,14 +40,7 @@ def test_random_2d(self):
         npt.assert_equal(len(vor.volumes), len(points))
         npt.assert_almost_equal(np.sum(vor.volumes), box.volume)
 
-        # Verify the neighbor distances
-        wrapped_distances = np.linalg.norm(
-            box.wrap(
-                points[vor.nlist.point_indices] - points[vor.nlist.query_point_indices]
-            ),
-            axis=-1,
-        )
-        npt.assert_allclose(wrapped_distances, vor.nlist.distances)
+        self._check_vectors_and_distances(vor, (box, points))
 
         # Ensure every point has neighbors
         assert np.all(vor.nlist.neighbor_counts > 0)
@@ -62,14 +70,7 @@ def test_random_3d(self):
         npt.assert_equal(len(vor.volumes), len(points))
         npt.assert_almost_equal(np.sum(vor.volumes), box.volume)
 
-        # Verify the neighbor distances
-        wrapped_distances = np.linalg.norm(
-            box.wrap(
-                points[vor.nlist.point_indices] - points[vor.nlist.query_point_indices]
-            ),
-            axis=-1,
-        )
-        npt.assert_allclose(wrapped_distances, vor.nlist.distances)
+        self._check_vectors_and_distances(vor, (box, points))
 
         # Ensure every point has neighbors
         assert np.all(vor.nlist.neighbor_counts > 0)
@@ -116,14 +117,7 @@ def test_voronoi_tess_2d(self):
             vor.nlist.weights[vor.nlist.query_point_indices == 4], 1
         )
 
-        # Verify the neighbor distances
-        wrapped_distances = np.linalg.norm(
-            box.wrap(
-                points[vor.nlist.point_indices] - points[vor.nlist.query_point_indices]
-            ),
-            axis=-1,
-        )
-        npt.assert_allclose(wrapped_distances, vor.nlist.distances)
+        self._check_vectors_and_distances(vor, (box, points))
 
         # Double the points (still inside the box) and test again
         points *= 2
@@ -139,14 +133,7 @@ def test_voronoi_tess_2d(self):
             vor.nlist.weights[vor.nlist.query_point_indices == 4], 2
         )
 
-        # Verify the neighbor distances
-        wrapped_distances = np.linalg.norm(
-            box.wrap(
-                points[vor.nlist.point_indices] - points[vor.nlist.query_point_indices]
-            ),
-            axis=-1,
-        )
-        npt.assert_allclose(wrapped_distances, vor.nlist.distances)
+        self._check_vectors_and_distances(vor, (box, points))
 
     def test_voronoi_tess_3d(self):
         # Test that the voronoi polytope works for a 3D system
@@ -211,14 +198,7 @@ def test_voronoi_tess_3d(self):
             vor.nlist.weights[vor.nlist.query_point_indices == 13], 1
         )
 
-        # Verify the neighbor distances
-        wrapped_distances = np.linalg.norm(
-            box.wrap(
-                points[vor.nlist.point_indices] - points[vor.nlist.query_point_indices]
-            ),
-            axis=-1,
-        )
-        npt.assert_allclose(wrapped_distances, vor.nlist.distances)
+        self._check_vectors_and_distances(vor, (box, points))
 
         # Double the points (still inside the box) and test again
         points *= 2
@@ -243,14 +223,7 @@ def test_voronoi_tess_3d(self):
             vor.nlist.weights[vor.nlist.query_point_indices == 13], 4
         )
 
-        # Verify the neighbor distances
-        wrapped_distances = np.linalg.norm(
-            box.wrap(
-                points[vor.nlist.point_indices] - points[vor.nlist.query_point_indices]
-            ),
-            axis=-1,
-        )
-        npt.assert_allclose(wrapped_distances, vor.nlist.distances)
+        self._check_vectors_and_distances(vor, (box, points))
 
     @pytest.mark.parametrize(
         "func, neighbors",
@@ -285,14 +258,7 @@ def test_voronoi_neighbors_wrapped(self, func, neighbors):
         npt.assert_equal(counts, neighbors)
         npt.assert_almost_equal(np.sum(vor.volumes), box.volume)
 
-        # Verify the neighbor distances
-        wrapped_distances = np.linalg.norm(
-            box.wrap(
-                points[vor.nlist.point_indices] - points[vor.nlist.query_point_indices]
-            ),
-            axis=-1,
-        )
-        npt.assert_allclose(wrapped_distances, vor.nlist.distances)
+        self._check_vectors_and_distances(vor, (box, points))
 
     def test_voronoi_weights_fcc(self):
         # Test that voronoi neighbor weights are computed properly for 3D FCC
@@ -323,14 +289,27 @@ def test_voronoi_weights_fcc(self):
             atol=1e-5,
         )
 
-        # Verify the neighbor distances
-        wrapped_distances = np.linalg.norm(
-            box.wrap(
-                points[vor.nlist.point_indices] - points[vor.nlist.query_point_indices]
-            ),
-            axis=-1,
+        self._check_vectors_and_distances(vor, (box, points))
+
+    def test_diamond(self):
+        box = freud.box.Box(
+            Lx=0.9075878262519836,
+            Ly=1.0908596515655518,
+            Lz=0.7142124772071838,
+            xy=0.6126953363418579,
+            xz=0.17466391623020172,
+            yz=0.16529279947280884,
         )
-        npt.assert_allclose(wrapped_distances, vor.nlist.distances)
+        points = np.array(
+            [
+                [0.6287098, 0.39470837, 0.33198223],
+                [-0.6287098, -0.39470834, -0.33198214],
+            ]
+        )
+        vor = freud.locality.Voronoi()
+        vor.compute((box, points))
+        nlist = vor.nlist
+        assert not np.any(np.all(np.isclose(nlist.vectors, [0, 0, 0]), axis=-1))
 
     def test_repr(self):
         vor = freud.locality.Voronoi()
diff --git a/tests/test_order_Hexatic.py b/tests/test_order_Hexatic.py
index 486767d23..f6e495490 100644
--- a/tests/test_order_Hexatic.py
+++ b/tests/test_order_Hexatic.py
@@ -171,7 +171,6 @@ def test_normalization(self, k):
         query_point_indices = np.array([0, 0, 0, 0])
         point_indices = np.array([1, 2, 3, 4])
         rijs = box.wrap(points[point_indices] - points[query_point_indices])
-        distances = np.linalg.norm(rijs, axis=-1)
         thetas = np.arctan2(rijs[:, 1], rijs[:, 0])
         weights = np.array([1, 0.7, 0.3, 0])
         nlist = freud.NeighborList.from_arrays(
@@ -179,7 +178,7 @@ def test_normalization(self, k):
             len(points),
             query_point_indices,
             point_indices,
-            distances,
+            rijs,
             weights,
         )
 
diff --git a/tests/test_order_Steinhardt.py b/tests/test_order_Steinhardt.py
index 28bfc2096..42e528a78 100644
--- a/tests/test_order_Steinhardt.py
+++ b/tests/test_order_Steinhardt.py
@@ -232,7 +232,7 @@ def test_weighted(self, wt):
                 len(positions),
                 nlist.query_point_indices,
                 nlist.point_indices,
-                nlist.distances,
+                nlist.vectors,
                 weights,
             )
 
@@ -301,14 +301,14 @@ def test_rotational_invariance(self, seed):
                 [0, 1, 1],
             ]
         )
-        query_point_indices = np.zeros(len(positions) - 1)
+        query_point_indices = np.zeros(len(positions) - 1, dtype=int)
         point_indices = np.arange(1, len(positions))
         nlist = freud.locality.NeighborList.from_arrays(
             len(positions),
             len(positions),
             query_point_indices,
             point_indices,
-            np.full(len(query_point_indices), np.sqrt(2)),
+            box.wrap(positions[point_indices] - positions[query_point_indices]),
         )
 
         q6 = freud.order.Steinhardt(6)
diff --git a/tests/validation/test_Minkowski_Structure_Metrics.py b/tests/validation/test_Minkowski_Structure_Metrics.py
index 6d153a00d..3270b6e08 100644
--- a/tests/validation/test_Minkowski_Structure_Metrics.py
+++ b/tests/validation/test_Minkowski_Structure_Metrics.py
@@ -48,12 +48,12 @@ def test_minkowski_structure_metrics(self, structure):
             # Test q'l
             comp = freud.order.Steinhardt(sph_l, weighted=True)
             comp.compute(snap, neighbors=voro.nlist)
-            npt.assert_allclose(comp.order, expected_ql[:, sph_l], atol=1e-5)
+            npt.assert_allclose(comp.order, expected_ql[:, sph_l], atol=1e-4)
 
             # Test average q'l
             comp = freud.order.Steinhardt(sph_l, average=True, weighted=True)
             comp.compute(snap, neighbors=voro.nlist)
-            npt.assert_allclose(comp.order, expected_avql[:, sph_l], atol=1e-5)
+            npt.assert_allclose(comp.order, expected_avql[:, sph_l], atol=1e-4)
 
             # w'2 tests fail for unknown (probably numerical) reasons.
             if sph_l != 2: