diff --git a/src/kanapy/api.py b/src/kanapy/api.py
index 6bda1e75..4d599d5b 100644
--- a/src/kanapy/api.py
+++ b/src/kanapy/api.py
@@ -212,9 +212,9 @@ def voxelize(self, particles=None, dim=None):
             print('Volume fractions of phases in voxel structure:')
             vt = 0.
             for ip in range(self.nphases):
-                vf = 100.0*vox_count[ip]/self.mesh.nvox
+                vf = vox_count[ip]/self.mesh.nvox
                 vt += vf
-                print(f'{ip}: {self.rve.phase_names[ip]} ({vf.round(1)}%)')
+                print(f'{ip}: {self.rve.phase_names[ip]} ({(vf*100):.3f}%)')
             if not np.isclose(vt, 1.0):
                 logging.warning(f'Volume fractions of phases in voxels do not add up to 1. Value: {vt}')
 
diff --git a/src/kanapy/collisions-alt.py b/src/kanapy/collisions-alt.py
new file mode 100644
index 00000000..3b512c09
--- /dev/null
+++ b/src/kanapy/collisions-alt.py
@@ -0,0 +1,249 @@
+# -*- coding: utf-8 -*-
+import numpy as np
+
+
+def collision_routine(E1, E2, damp=0):
+    """
+    Calls the method :meth:'collide_detect' to determine whether the given two ellipsoid objects overlap using
+    the Algebraic separation condition developed by W. Wang et al. A detailed description is provided
+    therein.
+
+    Also calls the :meth:`collision_react` to evaluate the response after collision.     
+
+    :param E1: Ellipsoid :math:`i` 
+    :type E1: object :obj:`Ellipsoid`
+    :param E2: Ellipsoid :math:`j`
+    :type E2: object :obj:`Ellipsoid`
+    :param damp: Damping factor
+    :type damp: float
+
+    .. note:: 1. If both the particles to be tested for overlap are spheres, then the bounding sphere hierarchy is
+                 sufficient to determine whether they overlap.
+              2. Else, if either of them is an ellipsoid, then their coefficients, positions & rotation matrices are
+                 used to determine whether they overlap.
+    """
+
+    # call the collision detection algorithm
+    overlap_status = collide_detect(E1.get_coeffs(), E2.get_coeffs(),
+                                    E1.get_pos(), E2.get_pos(),
+                                    E1.rotation_matrix, E2.rotation_matrix)
+
+    if overlap_status:
+        collision_react(E1, E2, damp=damp)  # calculates resultant speed for E1
+        collision_react(E2, E1, damp=damp)  # calculates resultant speed for E2
+
+    return overlap_status
+
+
+def collision_react(ell1, ell2, damp=0.):
+    r"""
+    Evaluates and modifies the magnitude and direction of the ellipsoid's velocity after collision.    
+
+    :param ell1: Ellipsoid :math:`i`
+    :type ell1: object :obj:`Ellipsoid`
+    :param ell2: Ellipsoid :math:`j`
+    :type ell2: object :obj:`Ellipsoid`
+    :param damp: Damping factor
+    :type damp: float
+
+    .. note:: Consider two ellipsoids :math:`i, j` at collision. Let them occupy certain positions in space
+              defined by the position vectors :math:`\mathbf{r}^{i}, \mathbf{r}^{j}` and have certain 
+              velocities represented by :math:`\mathbf{v}^{i}, \mathbf{v}^{j}` respectively. The objective
+              is to  find the velocity vectors after collision. The elevation angle :math:`\phi` between
+              the ellipsoids is determined by,      
+
+              .. image:: /figs/elevation_ell.png                        
+                :width: 200px
+                :height: 45px
+                :align: center                 
+
+              where :math:`dx, dy, dz` are defined as the distance between the two ellipsoid centers along
+                    :math:`x, y, z` directions given by,
+
+              .. image:: /figs/dist_ell.png                        
+                  :width: 110px
+                  :height: 75px
+                  :align: center
+
+              Depending on the magnitudes of :math:`dx, dz` as projected on the :math:`x-z` plane, the angle
+              :math:`\Theta` is computed.
+              The angles :math:`\Theta` and :math:`\phi` determine the in-plane and out-of-plane directions along which
+              the ellipsoid :math:`i` would bounce back after collision. Thus, the updated velocity vector components
+              along the :math:`x, y, z` directions are determined by,
+
+              .. image:: /figs/velocities_ell.png                        
+                  :width: 180px
+                  :height: 80px
+                  :align: center                        
+
+    ell1_speed = np.linalg.norm([ell1.speedx0, ell1.speedy0, ell1.speedz0]) * (1. - damp)
+    x_diff = ell2.x - ell1.x
+    y_diff = ell2.y - ell1.y
+    z_diff = ell2.z - ell1.z
+    elevation_angle = np.arctan2(y_diff, np.linalg.norm([x_diff, z_diff]))
+    angle = np.arctan2(z_diff, x_diff)
+
+    x_speed = -ell1_speed * np.cos(angle) * np.cos(elevation_angle)
+    y_speed = -ell1_speed * np.sin(elevation_angle)
+    z_speed = -ell1_speed * np.sin(angle) * np.cos(elevation_angle)
+
+    # Assign new speeds 
+    ell1.speedx += x_speed
+    ell1.speedy += y_speed
+    ell1.speedz += z_speed
+    """
+    # check if inner polyhedra overlap if present
+    """if overlap_status and E1.inner is not None and E2.inner is not None:
+        E1.sync_poly()
+        E2.sync_poly()
+        if all(E1.inner.find_simplex(E2.inner.points) < 0) and\
+           all(E2.inner.find_simplex(E1.inner.points) < 0):
+            overlap_status = False"""
+    x_diff = ell1.x - ell2.x
+    y_diff = ell1.y - ell2.y
+    z_diff = ell1.z - ell2.z
+    r2inv = 1.0/((x_diff/ell1.a) ** 3 + (y_diff/ell1.b) ** 3 + (z_diff/ell1.c) ** 3)
+    r2inv = np.sign(r2inv) * np.minimum(np.abs(r2inv), 2.)
+
+    elevation_angle = np.arctan2(y_diff, np.linalg.norm([x_diff, z_diff]))
+    angle = np.arctan2(z_diff, x_diff)
+
+    ell1.force_x += r2inv #* np.cos(angle) * np.cos(elevation_angle)
+    ell1.force_y += r2inv #* np.sin(elevation_angle)
+    ell1.force_z += r2inv #* np.sin(angle) * np.cos(elevation_angle)
+
+    return
+
+
+def collide_detect(coef_i, coef_j, r_i, r_j, A_i, A_j):
+    r"""Implementation of Algebraic separation condition developed by
+    W. Wang et al. 2001 for overlap detection between two static ellipsoids.
+    Kudos to ChatGPT for support with translation from C++ code.
+
+    :param coef_i: Coefficients of ellipsoid :math:`i`
+    :type coef_i: numpy array
+    :param coef_j: Coefficients of ellipsoid :math:`j`
+    :type coef_j: numpy array
+    :param r_i: Position of ellipsoid :math:`i`
+    :type r_i: numpy array
+    :param r_j: Position of ellipsoid :math:`j`
+    :type r_j: numpy array
+    :param A_i: Rotation matrix of ellipsoid :math:`i`
+    :type A_i: numpy array
+    :param A_j: Rotation matrix of ellipsoid :math:`j`
+    :type A_j: numpy array
+    :returns: **True** if ellipsoids :math:`i, j` overlap, else **False**
+    :rtype: boolean
+    """
+    SMALL = 1.e-12
+
+    # Initialize Matrices A & B with zeros
+    A = np.zeros((4, 4), dtype=float)
+    B = np.zeros((4, 4), dtype=float)
+
+    A[0, 0] = 1 / (coef_i[0] ** 2)
+    A[1, 1] = 1 / (coef_i[1] ** 2)
+    A[2, 2] = 1 / (coef_i[2] ** 2)
+    A[3, 3] = -1
+
+    B[0, 0] = 1 / (coef_j[0] ** 2)
+    B[1, 1] = 1 / (coef_j[1] ** 2)
+    B[2, 2] = 1 / (coef_j[2] ** 2)
+    B[3, 3] = -1
+
+    # Rigid body transformations
+    T_i = np.zeros((4, 4), dtype=float)
+    T_j = np.zeros((4, 4), dtype=float)
+
+    T_i[:3, :3] = A_i
+    T_i[:3, 3] = r_i
+    T_i[3, 3] = 1.0
+
+    T_j[:3, :3] = A_j
+    T_j[:3, 3] = r_j
+    T_j[3, 3] = 1.0
+
+    # Copy the arrays for future operations
+    Ma = np.tile(T_i, (1, 1))
+    Mb = np.tile(T_j, (1, 1))
+
+    # aij of matrix A in det(lambda*A - Ma'*(Mb^-1)'*B*(Mb^-1)*Ma).
+    # bij of matrix b = Ma'*(Mb^-1)'*B*(Mb^-1)*Ma
+    aux = np.linalg.inv(Mb) @ Ma
+    bm = aux.T @ B @ aux
+
+    # Coefficients of the Characteristic Polynomial.
+    T0 = (-A[0, 0] * A[1, 1] * A[2, 2])
+    T1 = (A[0, 0] * A[1, 1] * bm[2, 2] + A[0, 0] * A[2, 2] * bm[1, 1] + A[1, 1] * A[2, 2] * bm[0, 0] - A[0, 0] * A[
+        1, 1] * A[2, 2] * bm[3, 3])
+    T2 = (A[0, 0] * bm[1, 2] * bm[2, 1] - A[0, 0] * bm[1, 1] * bm[2, 2] - A[1, 1] * bm[0, 0] * bm[2, 2] + A[1, 1] * bm[
+        0, 2] * bm[2, 0] -
+          A[2, 2] * bm[0, 0] * bm[1, 1] + A[2, 2] * bm[0, 1] * bm[1, 0] + A[0, 0] * A[1, 1] * bm[2, 2] * bm[3, 3] - A[
+              0, 0] * A[1, 1] * bm[2, 3] * bm[3, 2] +
+          A[0, 0] * A[2, 2] * bm[1, 1] * bm[3, 3] - A[0, 0] * A[2, 2] * bm[1, 3] * bm[3, 1] + A[1, 1] * A[2, 2] * bm[
+              0, 0] * bm[3, 3] -
+          A[1, 1] * A[2, 2] * bm[0, 3] * bm[3, 0])
+    T3 = (bm[0, 0] * bm[1, 1] * bm[2, 2] - bm[0, 0] * bm[1, 2] * bm[2, 1] - bm[0, 1] * bm[1, 0] * bm[2, 2] + bm[0, 1] *
+          bm[1, 2] * bm[2, 0] +
+          bm[0, 2] * bm[1, 0] * bm[2, 1] - bm[0, 2] * bm[1, 1] * bm[2, 0] - A[0, 0] * bm[1, 1] * bm[2, 2] * bm[3, 3] +
+          A[0, 0] * bm[1, 1] * bm[2, 3] * bm[3, 2] +
+          A[0, 0] * bm[1, 2] * bm[2, 1] * bm[3, 3] - A[0, 0] * bm[1, 2] * bm[2, 3] * bm[3, 1] - A[0, 0] * bm[1, 3] * bm[
+              2, 1] * bm[3, 2] +
+          A[0, 0] * bm[1, 3] * bm[2, 2] * bm[3, 1] - A[1, 1] * bm[0, 0] * bm[2, 2] * bm[3, 3] + A[1, 1] * bm[0, 0] * bm[
+              2, 3] * bm[3, 2] +
+          A[1, 1] * bm[0, 2] * bm[2, 0] * bm[3, 3] - A[1, 1] * bm[0, 2] * bm[2, 3] * bm[3, 0] - A[1, 1] * bm[0, 3] * bm[
+              2, 0] * bm[3, 2] +
+          A[1, 1] * bm[0, 3] * bm[2, 2] * bm[3, 0] - A[2, 2] * bm[0, 0] * bm[1, 1] * bm[3, 3] + A[2, 2] * bm[0, 0] * bm[
+              1, 3] * bm[3, 1] +
+          A[2, 2] * bm[0, 1] * bm[1, 0] * bm[3, 3] - A[2, 2] * bm[0, 1] * bm[1, 3] * bm[3, 0] - A[2, 2] * bm[0, 3] * bm[
+              1, 0] * bm[3, 1] +
+          A[2, 2] * bm[0, 3] * bm[1, 1] * bm[3, 0])
+    T4 = (bm[0, 0] * bm[1, 1] * bm[2, 2] * bm[3, 3] - bm[0, 0] * bm[1, 1] * bm[2, 3] * bm[3, 2] - bm[0, 0] * bm[1, 2] *
+          bm[2, 1] * bm[3, 3] +
+          bm[0, 0] * bm[1, 2] * bm[2, 3] * bm[3, 1] + bm[0, 0] * bm[1, 3] * bm[2, 1] * bm[3, 2] - bm[0, 0] * bm[1, 3] *
+          bm[2, 2] * bm[3, 1] -
+          bm[0, 1] * bm[1, 0] * bm[2, 2] * bm[3, 3] + bm[0, 1] * bm[1, 0] * bm[2, 3] * bm[3, 2] + bm[0, 1] * bm[1, 2] *
+          bm[2, 0] * bm[3, 3] -
+          bm[0, 1] * bm[1, 2] * bm[2, 3] * bm[3, 0] - bm[0, 1] * bm[1, 3] * bm[2, 0] * bm[3, 2] + bm[0, 1] * bm[1, 3] *
+          bm[2, 2] * bm[3, 0] +
+          bm[0, 2] * bm[1, 0] * bm[2, 1] * bm[3, 3] - bm[0, 2] * bm[1, 0] * bm[2, 3] * bm[3, 1] - bm[0, 2] * bm[1, 1] *
+          bm[2, 0] * bm[3, 3] +
+          bm[0, 2] * bm[1, 1] * bm[2, 3] * bm[3, 0] + bm[0, 2] * bm[1, 3] * bm[2, 0] * bm[3, 1] - bm[0, 2] * bm[1, 3] *
+          bm[2, 1] * bm[3, 0] -
+          bm[0, 3] * bm[1, 0] * bm[2, 1] * bm[3, 2] + bm[0, 3] * bm[1, 0] * bm[2, 2] * bm[3, 1] + bm[0, 3] * bm[1, 1] *
+          bm[2, 0] * bm[3, 2] -
+          bm[0, 3] * bm[1, 1] * bm[2, 2] * bm[3, 0] - bm[0, 3] * bm[1, 2] * bm[2, 0] * bm[3, 1] + bm[0, 3] * bm[1, 2] *
+          bm[2, 1] * bm[3, 0])
+
+    # Roots of the characteristic_polynomial (lambda0, ... , lambda4).
+    cp = np.array([T0, T1, T2, T3, T4], dtype=float)
+
+    # Solve the polynomial
+    roots = np.roots(cp)
+
+    # Find the real roots where imaginary part doesn't exist
+    real_roots = [root.real for root in roots if abs(root.imag) < SMALL]
+
+    # Count number of real negative roots
+    count_neg = sum(1 for root in real_roots if root < -SMALL)
+
+    # Sort the real roots in ascending order
+    real_roots.sort()
+
+    # Algebraic separation conditions to determine overlapping
+    if count_neg == 2:
+        if abs(real_roots[0] - real_roots[1]) > SMALL:
+            return False
+        else:
+            return True
+    else:
+        return True
+
+# Example usage:
+# coef_i = np.array([a_i, b_i, c_i])
+# coef_j = np.array([a_j, b_j, c_j])
+# r_i = np.array([x_i, y_i, z_i])
+# r_j = np.array([x_j, y_j, z_j])
+# A_i = np.array([[A11_i, A12_i, A13_i], [A21_i, A22_i, A23_i], [A31_i, A32_i, A33_i]])
+# A_j = np.array([[A11_j, A12_j, A13_j], [A21_j, A22_j, A23_j], [A31_j, A32_j, A33_j]])
+# result = collide_detect(coef_i, coef_j, r_i, r_j, A_i, A_j)
diff --git a/src/kanapy/collisions-hull.py b/src/kanapy/collisions-hull.py
new file mode 100644
index 00000000..09faf519
--- /dev/null
+++ b/src/kanapy/collisions-hull.py
@@ -0,0 +1,258 @@
+# -*- coding: utf-8 -*-
+import time
+import numpy as np
+from scipy.spatial import Delaunay
+
+
+def collision_routine(E1, E2):
+    """
+    Calls the method :meth:'collide_detect' to determine whether the given two ellipsoid objects overlap using
+    the Algebraic separation condition developed by W. Wang et al. A detailed description is provided
+    therein.
+
+    Also calls the :meth:`collision_react` to evaluate the response after collision.     
+
+    :param E1: Ellipsoid :math:`i` 
+    :type E1: object :obj:`Ellipsoid`
+    :param E2: Ellipsoid :math:`j`
+    :type E2: object :obj:`Ellipsoid`
+
+    .. note:: 1. If both the particles to be tested for overlap are spheres, then the bounding sphere hierarchy is
+                 sufficient to determine whether they overlap.
+              2. Else, if either of them is an ellipsoid, then their coefficients, positions & rotation matrices are
+                 used to determine whether they overlap.
+    """
+
+    # call the collision detection algorithm
+    start = time.time()
+    overlap_status = collide_detect(E1.get_coeffs(), E2.get_coeffs(),
+                                    E1.get_pos(), E2.get_pos(),
+                                    E1.rotation_matrix, E2.rotation_matrix)
+    end = time.time()
+    tcd = end-start
+    start = time.time()
+    if overlap_status:
+        overlap_status = collision_react(E1, E2)  # calculates resultant contact forces of both particles
+    end = time.time()
+    tcr = end-start
+
+    return overlap_status, tcd, tcr
+
+
+def collision_react(ell1, ell2):
+    r"""
+    Evaluates and modifies the magnitude and direction of the ellipsoid's velocity after collision.    
+
+    :param ell1: Ellipsoid :math:`i`
+    :type ell1: object :obj:`Ellipsoid`
+    :param ell2: Ellipsoid :math:`j`
+    :type ell2: object :obj:`Ellipsoid`
+
+    .. note:: Consider two ellipsoids :math:`i, j` at collision. Let them occupy certain positions in space
+              defined by the position vectors :math:`\mathbf{r}^{i}, \mathbf{r}^{j}` and have certain 
+              velocities represented by :math:`\mathbf{v}^{i}, \mathbf{v}^{j}` respectively. The objective
+              is to  find the velocity vectors after collision. The elevation angle :math:`\phi` between
+              the ellipsoids is determined by,      
+
+              .. image:: /figs/elevation_ell.png                        
+                :width: 200px
+                :height: 45px
+                :align: center                 
+
+              where :math:`dx, dy, dz` are defined as the distance between the two ellipsoid centers along
+                    :math:`x, y, z` directions given by,
+
+              .. image:: /figs/dist_ell.png                        
+                  :width: 110px
+                  :height: 75px
+                  :align: center
+
+              Depending on the magnitudes of :math:`dx, dz` as projected on the :math:`x-z` plane, the angle
+              :math:`\Theta` is computed.
+              The angles :math:`\Theta` and :math:`\phi` determine the in-plane and out-of-plane directions along which
+              the ellipsoid :math:`i` would bounce back after collision. Thus, the updated velocity vector components
+              along the :math:`x, y, z` directions are determined by,
+
+              .. image:: /figs/velocities_ell.png                        
+                  :width: 180px
+                  :height: 80px
+                  :align: center                        
+
+    """
+    if ell1.inner is None:
+        pts1 = ell1.surfacePointsGen(nang=180)
+        hull1 = Delaunay(pts1)
+    else:
+        ell1.sync_poly()
+        pts1 = ell1.inner.points
+        hull1 = ell1.inner
+    if ell2.inner is None:
+        pts2 = ell2.surfacePointsGen(nang=180)
+        hull2 = Delaunay(pts2)
+    else:
+        ell2.sync_poly()
+        pts2 = ell2.inner.points
+        hull2 = ell2.inner
+
+    ind12 = np.nonzero(hull1.find_simplex(pts2) >= 0)[0]
+    if len(ind12) == 0:
+        # no overlapping vertices on convex hulls
+        return False
+    ind21 = np.nonzero(hull2.find_simplex(pts1) >= 0)[0]
+    if len(ind21) == 0:
+        return False
+
+    pos_i1 = np.average(pts1[ind21], axis=0)
+    pos_i2 = np.average(pts2[ind12], axis=0)
+    ctr1 = ell1.get_pos()
+    ctr2 = ell2.get_pos()
+    dir1 = ctr1 - pos_i1
+    dir2 = ctr2 - pos_i2
+    dst1 = np.linalg.norm(dir1)
+    dst2 = np.linalg.norm(dir2)
+    fac = np.average([dst1, dst2])
+    val = np.linalg.norm(pos_i1 - pos_i2) / fac
+    val = np.minimum(5.0, val)
+
+    ell1.force_x += val * dir1[0] / dst1
+    ell1.force_y += val * dir1[1] / dst1
+    ell1.force_z += val * dir1[2] / dst1
+    ell2.force_x += val * dir2[0] / dst2
+    ell2.force_y += val * dir2[1] / dst2
+    ell2.force_z += val * dir2[2] / dst2
+
+    return True
+
+
+def collide_detect(coef_i, coef_j, r_i, r_j, A_i, A_j):
+    r"""Implementation of Algebraic separation condition developed by
+    W. Wang et al. 2001 for overlap detection between two static ellipsoids.
+    Kudos to ChatGPT for support with translation from C++ code.
+
+    :param coef_i: Coefficients of ellipsoid :math:`i`
+    :type coef_i: numpy array
+    :param coef_j: Coefficients of ellipsoid :math:`j`
+    :type coef_j: numpy array
+    :param r_i: Position of ellipsoid :math:`i`
+    :type r_i: numpy array
+    :param r_j: Position of ellipsoid :math:`j`
+    :type r_j: numpy array
+    :param A_i: Rotation matrix of ellipsoid :math:`i`
+    :type A_i: numpy array
+    :param A_j: Rotation matrix of ellipsoid :math:`j`
+    :type A_j: numpy array
+    :returns: **True** if ellipsoids :math:`i, j` overlap, else **False**
+    :rtype: boolean
+    """
+    SMALL = 1.e-12
+
+    # Initialize Matrices A & B with zeros
+    A = np.zeros((4, 4), dtype=float)
+    B = np.zeros((4, 4), dtype=float)
+
+    A[0, 0] = 1 / (coef_i[0] ** 2)
+    A[1, 1] = 1 / (coef_i[1] ** 2)
+    A[2, 2] = 1 / (coef_i[2] ** 2)
+    A[3, 3] = -1
+
+    B[0, 0] = 1 / (coef_j[0] ** 2)
+    B[1, 1] = 1 / (coef_j[1] ** 2)
+    B[2, 2] = 1 / (coef_j[2] ** 2)
+    B[3, 3] = -1
+
+    # Rigid body transformations
+    T_i = np.zeros((4, 4), dtype=float)
+    T_j = np.zeros((4, 4), dtype=float)
+
+    T_i[:3, :3] = A_i
+    T_i[:3, 3] = r_i
+    T_i[3, 3] = 1.0
+
+    T_j[:3, :3] = A_j
+    T_j[:3, 3] = r_j
+    T_j[3, 3] = 1.0
+
+    # Copy the arrays for future operations
+    Ma = np.tile(T_i, (1, 1))
+    Mb = np.tile(T_j, (1, 1))
+
+    # aij of matrix A in det(lambda*A - Ma'*(Mb^-1)'*B*(Mb^-1)*Ma).
+    # bij of matrix b = Ma'*(Mb^-1)'*B*(Mb^-1)*Ma
+    aux = np.linalg.inv(Mb) @ Ma
+    bm = aux.T @ B @ aux
+
+    # Coefficients of the Characteristic Polynomial.
+    T0 = (-A[0, 0] * A[1, 1] * A[2, 2])
+    T1 = (A[0, 0] * A[1, 1] * bm[2, 2] + A[0, 0] * A[2, 2] * bm[1, 1] + A[1, 1] * A[2, 2] * bm[0, 0] - A[0, 0] * A[
+        1, 1] * A[2, 2] * bm[3, 3])
+    T2 = (A[0, 0] * bm[1, 2] * bm[2, 1] - A[0, 0] * bm[1, 1] * bm[2, 2] - A[1, 1] * bm[0, 0] * bm[2, 2] + A[1, 1] * bm[
+        0, 2] * bm[2, 0] -
+          A[2, 2] * bm[0, 0] * bm[1, 1] + A[2, 2] * bm[0, 1] * bm[1, 0] + A[0, 0] * A[1, 1] * bm[2, 2] * bm[3, 3] - A[
+              0, 0] * A[1, 1] * bm[2, 3] * bm[3, 2] +
+          A[0, 0] * A[2, 2] * bm[1, 1] * bm[3, 3] - A[0, 0] * A[2, 2] * bm[1, 3] * bm[3, 1] + A[1, 1] * A[2, 2] * bm[
+              0, 0] * bm[3, 3] -
+          A[1, 1] * A[2, 2] * bm[0, 3] * bm[3, 0])
+    T3 = (bm[0, 0] * bm[1, 1] * bm[2, 2] - bm[0, 0] * bm[1, 2] * bm[2, 1] - bm[0, 1] * bm[1, 0] * bm[2, 2] + bm[0, 1] *
+          bm[1, 2] * bm[2, 0] +
+          bm[0, 2] * bm[1, 0] * bm[2, 1] - bm[0, 2] * bm[1, 1] * bm[2, 0] - A[0, 0] * bm[1, 1] * bm[2, 2] * bm[3, 3] +
+          A[0, 0] * bm[1, 1] * bm[2, 3] * bm[3, 2] +
+          A[0, 0] * bm[1, 2] * bm[2, 1] * bm[3, 3] - A[0, 0] * bm[1, 2] * bm[2, 3] * bm[3, 1] - A[0, 0] * bm[1, 3] * bm[
+              2, 1] * bm[3, 2] +
+          A[0, 0] * bm[1, 3] * bm[2, 2] * bm[3, 1] - A[1, 1] * bm[0, 0] * bm[2, 2] * bm[3, 3] + A[1, 1] * bm[0, 0] * bm[
+              2, 3] * bm[3, 2] +
+          A[1, 1] * bm[0, 2] * bm[2, 0] * bm[3, 3] - A[1, 1] * bm[0, 2] * bm[2, 3] * bm[3, 0] - A[1, 1] * bm[0, 3] * bm[
+              2, 0] * bm[3, 2] +
+          A[1, 1] * bm[0, 3] * bm[2, 2] * bm[3, 0] - A[2, 2] * bm[0, 0] * bm[1, 1] * bm[3, 3] + A[2, 2] * bm[0, 0] * bm[
+              1, 3] * bm[3, 1] +
+          A[2, 2] * bm[0, 1] * bm[1, 0] * bm[3, 3] - A[2, 2] * bm[0, 1] * bm[1, 3] * bm[3, 0] - A[2, 2] * bm[0, 3] * bm[
+              1, 0] * bm[3, 1] +
+          A[2, 2] * bm[0, 3] * bm[1, 1] * bm[3, 0])
+    T4 = (bm[0, 0] * bm[1, 1] * bm[2, 2] * bm[3, 3] - bm[0, 0] * bm[1, 1] * bm[2, 3] * bm[3, 2] - bm[0, 0] * bm[1, 2] *
+          bm[2, 1] * bm[3, 3] +
+          bm[0, 0] * bm[1, 2] * bm[2, 3] * bm[3, 1] + bm[0, 0] * bm[1, 3] * bm[2, 1] * bm[3, 2] - bm[0, 0] * bm[1, 3] *
+          bm[2, 2] * bm[3, 1] -
+          bm[0, 1] * bm[1, 0] * bm[2, 2] * bm[3, 3] + bm[0, 1] * bm[1, 0] * bm[2, 3] * bm[3, 2] + bm[0, 1] * bm[1, 2] *
+          bm[2, 0] * bm[3, 3] -
+          bm[0, 1] * bm[1, 2] * bm[2, 3] * bm[3, 0] - bm[0, 1] * bm[1, 3] * bm[2, 0] * bm[3, 2] + bm[0, 1] * bm[1, 3] *
+          bm[2, 2] * bm[3, 0] +
+          bm[0, 2] * bm[1, 0] * bm[2, 1] * bm[3, 3] - bm[0, 2] * bm[1, 0] * bm[2, 3] * bm[3, 1] - bm[0, 2] * bm[1, 1] *
+          bm[2, 0] * bm[3, 3] +
+          bm[0, 2] * bm[1, 1] * bm[2, 3] * bm[3, 0] + bm[0, 2] * bm[1, 3] * bm[2, 0] * bm[3, 1] - bm[0, 2] * bm[1, 3] *
+          bm[2, 1] * bm[3, 0] -
+          bm[0, 3] * bm[1, 0] * bm[2, 1] * bm[3, 2] + bm[0, 3] * bm[1, 0] * bm[2, 2] * bm[3, 1] + bm[0, 3] * bm[1, 1] *
+          bm[2, 0] * bm[3, 2] -
+          bm[0, 3] * bm[1, 1] * bm[2, 2] * bm[3, 0] - bm[0, 3] * bm[1, 2] * bm[2, 0] * bm[3, 1] + bm[0, 3] * bm[1, 2] *
+          bm[2, 1] * bm[3, 0])
+
+    # Roots of the characteristic_polynomial (lambda0, ... , lambda4).
+    cp = np.array([T0, T1, T2, T3, T4], dtype=float)
+
+    # Solve the polynomial
+    roots = np.roots(cp)
+
+    # Find the real roots where imaginary part doesn't exist
+    real_roots = [root.real for root in roots if abs(root.imag) < SMALL]
+
+    # Count number of real negative roots
+    count_neg = sum(1 for root in real_roots if root < -SMALL)
+
+    # Sort the real roots in ascending order
+    real_roots.sort()
+
+    # Algebraic separation conditions to determine overlapping
+    if count_neg == 2:
+        if abs(real_roots[0] - real_roots[1]) > SMALL:
+            return False
+        else:
+            return True
+    else:
+        return True
+
+# Example usage:
+# coef_i = np.array([a_i, b_i, c_i])
+# coef_j = np.array([a_j, b_j, c_j])
+# r_i = np.array([x_i, y_i, z_i])
+# r_j = np.array([x_j, y_j, z_j])
+# A_i = np.array([[A11_i, A12_i, A13_i], [A21_i, A22_i, A23_i], [A31_i, A32_i, A33_i]])
+# A_j = np.array([[A11_j, A12_j, A13_j], [A21_j, A22_j, A23_j], [A31_j, A32_j, A33_j]])
+# result = collide_detect(coef_i, coef_j, r_i, r_j, A_i, A_j)
diff --git a/src/kanapy/collisions.py b/src/kanapy/collisions.py
index dbb6ec65..1284ee6d 100644
--- a/src/kanapy/collisions.py
+++ b/src/kanapy/collisions.py
@@ -2,7 +2,7 @@
 import numpy as np
 
 
-def collision_routine(E1, E2, damp=0):
+def collision_routine(E1, E2):
     """
     Calls the method :meth:'collide_detect' to determine whether the given two ellipsoid objects overlap using
     the Algebraic separation condition developed by W. Wang et al. A detailed description is provided
@@ -14,8 +14,6 @@ def collision_routine(E1, E2, damp=0):
     :type E1: object :obj:`Ellipsoid`
     :param E2: Ellipsoid :math:`j`
     :type E2: object :obj:`Ellipsoid`
-    :param damp: Damping factor
-    :type damp: float
 
     .. note:: 1. If both the particles to be tested for overlap are spheres, then the bounding sphere hierarchy is
                  sufficient to determine whether they overlap.
@@ -27,23 +25,12 @@ def collision_routine(E1, E2, damp=0):
     overlap_status = collide_detect(E1.get_coeffs(), E2.get_coeffs(),
                                     E1.get_pos(), E2.get_pos(),
                                     E1.rotation_matrix, E2.rotation_matrix)
-
-    # check if inner polyhedra overlap if present
-    """if overlap_status and E1.inner is not None and E2.inner is not None:
-        E1.sync_poly()
-        E2.sync_poly()
-        if all(E1.inner.find_simplex(E2.inner.points) < 0) and\
-           all(E2.inner.find_simplex(E1.inner.points) < 0):
-            overlap_status = False"""
-
     if overlap_status:
-        collision_react(E1, E2, damp=damp)  # calculates resultant speed for E1
-        collision_react(E2, E1, damp=damp)  # calculates resultant speed for E2
-
+        collision_react(E1, E2)  # calculates resultant contact forces of both particles
     return overlap_status
 
 
-def collision_react(ell1, ell2, damp=0.):
+def collision_react(ell1, ell2):
     r"""
     Evaluates and modifies the magnitude and direction of the ellipsoid's velocity after collision.    
 
@@ -51,8 +38,6 @@ def collision_react(ell1, ell2, damp=0.):
     :type ell1: object :obj:`Ellipsoid`
     :param ell2: Ellipsoid :math:`j`
     :type ell2: object :obj:`Ellipsoid`
-    :param damp: Damping factor
-    :type damp: float
 
     .. note:: Consider two ellipsoids :math:`i, j` at collision. Let them occupy certain positions in space
               defined by the position vectors :math:`\mathbf{r}^{i}, \mathbf{r}^{j}` and have certain 
@@ -83,23 +68,21 @@ def collision_react(ell1, ell2, damp=0.):
                   :width: 180px
                   :height: 80px
                   :align: center                        
-    """
-    ell1_speed = np.linalg.norm([ell1.speedx0, ell1.speedy0, ell1.speedz0]) * (1. - damp)
-    x_diff = ell2.x - ell1.x
-    y_diff = ell2.y - ell1.y
-    z_diff = ell2.z - ell1.z
-    elevation_angle = np.arctan2(y_diff, np.linalg.norm([x_diff, z_diff]))
-    angle = np.arctan2(z_diff, x_diff)
-
-    x_speed = -ell1_speed * np.cos(angle) * np.cos(elevation_angle)
-    y_speed = -ell1_speed * np.sin(elevation_angle)
-    z_speed = -ell1_speed * np.sin(angle) * np.cos(elevation_angle)
-
-    # Assign new speeds 
-    ell1.speedx += x_speed
-    ell1.speedy += y_speed
-    ell1.speedz += z_speed
 
+    """
+    cdir = ell1.get_pos() - ell2.get_pos()
+    dst = np.linalg.norm(cdir)
+    eqd1 = ell1.get_volume()**(1/3)
+    eqd2 = ell2.get_volume()**(1/3)
+    val = (0.5*(eqd1 + eqd2) / dst)**3
+    val = np.minimum(5.0, val)
+
+    ell1.force_x += val * cdir[0] / dst
+    ell1.force_y += val * cdir[1] / dst
+    ell1.force_z += val * cdir[2] / dst
+    ell2.force_x -= val * cdir[0] / dst
+    ell2.force_y -= val * cdir[1] / dst
+    ell2.force_z -= val * cdir[2] / dst
     return
 
 
diff --git a/src/kanapy/entities.py b/src/kanapy/entities.py
index 95481184..23371876 100644
--- a/src/kanapy/entities.py
+++ b/src/kanapy/entities.py
@@ -2,7 +2,7 @@
 import logging
 import numpy as np
 import random
-from kanapy.collisions import collision_routine
+from kanapy.collisions import collision_react, collision_routine
 from scipy.spatial import Delaunay
 
 
@@ -89,14 +89,14 @@ def __init__(self, iden, x, y, z, a, b, c, quat, phasenum=0, dup=None, points=No
         self.x = x
         self.y = y
         self.z = z
+        self.xold = x
+        self.yold = y
+        self.zold = z
         self.a = a
         self.b = b
         self.c = c
         self.quat = quat
         self.oria, self.orib, self.oric = a, b, c  # Store the original size of the particle
-        self.speedx0 = 0.
-        self.speedy0 = 0.
-        self.speedz0 = 0.
         self.speedx = 0.
         self.speedy = 0.
         self.speedz = 0.
@@ -277,20 +277,24 @@ def sync_poly(self, scale=1.6):
         if any(self.inner.find_simplex(self.surface_points) >= 0):
             logging.error(f'Polyhedron too large for ellipsoid {self.id}. Reduce scale.')"""
 
-
-    def move(self):
+    def move(self, dt):
         """
-        Moves the ellipsoid by updating its position vector according to Eulerian integration method 
+        Moves the ellipsoid by updating its position vector according to the Verlet integration method
 
         .. note:: The :class:`~Cuboid` object of the ellipsoid has to be updated everytime it moves
         """
-        self.x += self.speedx
-        self.y += self.speedy
-        self.z += self.speedz
-
-        self.speedx += self.force_x
-        self.speedy += self.force_y
-        self.speedz += self.force_z
+        xx = 2.0 * self.x - self.xold + self.force_x * dt * dt
+        yy = 2.0 * self.y - self.yold + self.force_y * dt * dt
+        zz = 2.0 * self.z - self.zold + self.force_z * dt * dt
+        self.speedx = (xx - self.xold) / (2.0 * dt)
+        self.speedy = (yy - self.yold) / (2.0 * dt)
+        self.speedz = (zz - self.zold) / (2.0 * dt)
+        self.xold = self.x
+        self.yold = self.y
+        self.zold = self.z
+        self.x = xx
+        self.y = yy
+        self.z = zz
         self.set_cub()
 
     def gravity_effect(self, value):
@@ -713,32 +717,32 @@ def wallCollision(self, sim_box, periodicity):
                 diff = sim_box.left - self.bbox_xmin
                 # move the ellipsoid in opposite direction after bouncing
                 self.x += diff
-                self.speedx *= -1
+                self.xold = self.x
 
             if self.bbox_ymin < sim_box.top:
                 diff = sim_box.top - self.bbox_ymin
                 self.y += diff
-                self.speedy *= -1
+                self.yold = self.y
 
             if self.bbox_zmin < sim_box.front:
                 diff = sim_box.front - self.bbox_zmin
                 self.z += diff
-                self.speedz *= -1
+                self.zold = self.z
 
             if self.bbox_xmax > sim_box.right:
                 diff = self.bbox_xmax - sim_box.right
                 self.x -= diff
-                self.speedx *= -1
+                self.xold = self.x
 
             if self.bbox_ymax > sim_box.bottom:
                 diff = self.bbox_ymax - sim_box.bottom
                 self.y -= diff
-                self.speedy *= -1
+                self.yold = self.y
 
             if self.bbox_zmax > sim_box.back:
                 diff = self.bbox_zmax - sim_box.back
                 self.z -= diff
-                self.speedz *= -1
+                self.zold = self.z
         return duplicates
 
 
@@ -856,15 +860,22 @@ def make_neighborlist(self):
         Finds the neighborlist for each particle
         """
         for particle in self.particles:
+            particle.neighborlist = set()
             for branch in particle.branches:
                 particle.neighborlist.update(branch.particles)
 
-    def collisionsTest(self, damp=0.):
+    def collisionsTest(self):
         """
         Tests for collision between all ellipsoids in the particle list of octree         
         """
         self.make_neighborlist()
         ncoll = 0
+        count = 0
+        colp = []
+        ppar = []
+        timecd = 0.
+        timecr = 0.
+        third = 1.0 / 3.0
         for E1 in self.particles:
             for E2 in E1.neighborlist:
                 id1 = E1.id if E1.duplicate is None else (E1.duplicate + len(self.particles))
@@ -872,13 +883,21 @@ def collisionsTest(self, damp=0.):
                 if id2 > id1:
                     # Distance between the centers of ellipsoids
                     dist = np.linalg.norm(np.subtract(E1.get_pos(), E2.get_pos()))
+                    psize = 0.5 * (E1.get_volume() ** third + E2.get_volume() ** third)
                     # If the bounding spheres collide then check for collision
-                    if dist <= (E1.a + E2.a):
+                    if dist <= psize:
                         # Check if ellipsoids overlap and update their speeds accordingly
-                        if collision_routine(E1, E2, damp=damp):
+                        count += 1
+                        ppar.append([dist, psize])
+                        if collision_routine(E1, E2):
+                            colp.append([dist, psize])
                             ncoll += 1
                             E1.ncollision += 1
                             E2.ncollision += 1
+        # print(f'Checked {count} pairs with {ncoll} collisions.')
+        # print(f'collision time: {timecd}, force time: {timecr}')
+        # print(f'Pair params:', ppar)
+        # print(f'Collisions:', colp)
         return ncoll
 
     def update(self):
diff --git a/src/kanapy/packing.py b/src/kanapy/packing.py
index bca821cf..036e1777 100644
--- a/src/kanapy/packing.py
+++ b/src/kanapy/packing.py
@@ -5,7 +5,7 @@
 
 from kanapy.input_output import write_dump
 from kanapy.entities import Ellipsoid, Cuboid, Octree
-from kanapy.collisions import collide_detect
+from kanapy.collisions import collision_routine, collide_detect
 
 
 def particle_generator(particle_data, sim_box, periodic, poly):
@@ -30,7 +30,7 @@ def particle_generator(particle_data, sim_box, periodic, poly):
     id_ctr = 0
     vell = 0.
     # introduce scaling factor to reduce particle overlap for non-periodic box
-    sf = 0.5 # if periodic else 0.45
+    sf = 0.5  # if periodic else 0.45
     for particle in particle_data:
         num_particles = particle['Number']  # Total number of ellipsoids
         for n in range(num_particles):
@@ -71,16 +71,7 @@ def particle_generator(particle_data, sim_box, periodic, poly):
                                   phasenum=particle['Phase'], points=poly)
             ellipsoid.color = (random.randint(0, 255), random.randint(
                 0, 255), random.randint(0, 255))
-
-            # Define random speed values along the 3 axes x, y & z
-            ellipsoid.speedx0 = np.random.uniform(low=-c / 20., high=c / 20.)
-            ellipsoid.speedy0 = np.random.uniform(low=-c / 20., high=c / 20.)
-            ellipsoid.speedz0 = np.random.uniform(low=-c / 20., high=c / 20.)
-            ellipsoid.speedx = ellipsoid.speedx0
-            ellipsoid.speedy = ellipsoid.speedy0
-            ellipsoid.speedz = ellipsoid.speedz0
             vell += ellipsoid.get_volume()
-
             Ellipsoids.append(ellipsoid)  # adds ellipsoid to list
         id_ctr += num_particles
         print(f'    Total volume of generated ellipsoids: {vell}')
@@ -122,47 +113,94 @@ def particle_grow(sim_box, Ellipsoids, periodicity, nsteps,
               By default, periodic images are written to the output file, 
               but this option can be disabled within the function.         
     """
+
+    def t_step(N):
+        return K * N ** m
+
+    def stop_part(ell):
+        ell.speedx = 0.
+        ell.speedy = 0.
+        ell.speedz = 0.
+        ell.force_x = 0.
+        ell.force_y = 0.
+        ell.force_z = 0
+        ell.x *= 0.9
+        ell.y *= 0.9
+        ell.z *= 0.9
+        ell.oldx = ell.x
+        ell.oldy = ell.y
+        ell.oldz = ell.z
+        return
+
     if fill_factor is None:
         fill_factor = 0.65  # 65% should be largest packing density of ellipsoids
     # Reduce the volume of the particles to (1/nsteps)th of its original value
     end_step = int(fill_factor * nsteps) - 1  # grow particles only to given volume fraction
-    fac = nsteps**(-1/3)
+    m = -1 / 2.5
+    K = 5
+    Niter = 1
+    time = 0
+    while time < end_step:
+        time += t_step(Niter)
+        Niter += 1
+    fac = nsteps ** (-1 / 3)
     vorig = 0.
     ve = 0.
     for ell in Ellipsoids:
         ell.a, ell.b, ell.c = ell.oria * fac, ell.orib * fac, ell.oric * fac
         ve += ell.get_volume()
-        vorig += ell.oria*ell.orib*ell.oric*np.pi*4/3
+        vorig += ell.oria * ell.orib * ell.oric * np.pi * 4 / 3
     dv = vorig / nsteps
-    print(f'Volume of simulation box: {sim_box.w*sim_box.h*sim_box.d}')
+    print(f'Volume of simulation box: {sim_box.w * sim_box.h * sim_box.d}')
     print(f'Volume of unscaled particles: {vorig}')
-    print(f'Initial volume of scaled ellipsoids: {ve}, targeted final volume: {ve+end_step*dv}')
+    print(f'Initial volume of scaled ellipsoids: {ve}, targeted final volume: {ve + end_step * dv}')
     print(f'Volume increment per time step: {dv}')
 
     # Simulation loop for particle growth and interaction steps
-    for i in tqdm(range(end_step)):
-
+    damping = 0.4
+    ncol = 0
+    ndump = np.maximum(int(Niter / 1000), 1)
+    time = 0
+    for i in tqdm(range(1, Niter)):
+        dt = t_step(i)
+        time += dt
         # Initialize Octree and test for collision between ellipsoids
-        for ellipsoid in Ellipsoids:
-            ellipsoid.branches = []
-            if periodicity:
-                ellipsoid.speedx = 0.
-                ellipsoid.speedy = 0.
-                ellipsoid.speedz = 0.
-            ellipsoid.ncollision = 0
+
+        ekin = 0.
+        for ell in Ellipsoids:
+            # apply damping force opposite to current movement
+            ell.force_x = -damping * ell.speedx / dt
+            ell.force_y = -damping * ell.speedy / dt
+            ell.force_z = -damping * ell.speedz / dt
+            ekin += ell.speedx ** 2 + ell.speedy ** 2 + ell.speedz ** 2
+            ell.ncollision = 0
+            ell.branches = []
 
         tree = Octree(0, Cuboid(sim_box.left, sim_box.top, sim_box.right,
-                                sim_box.bottom, sim_box.front, sim_box.back),
-                      Ellipsoids)
+                      sim_box.bottom, sim_box.front, sim_box.back), Ellipsoids)
+        if (not np.isclose(k_att, 0.0)) or (not np.isclose(k_rep, 0.0)):
+            calculateForce(Ellipsoids, sim_box, periodicity,
+                           k_rep=k_rep, k_att=k_att)
         tree.update()
-
-        if periodicity:
-            for ellipsoid in Ellipsoids:
-                if ellipsoid.ncollision == 0:
-                    ellipsoid.speedx = ellipsoid.speedx0
-                    ellipsoid.speedy = ellipsoid.speedy0
-                    ellipsoid.speedz = ellipsoid.speedz0
-        if dump:
+        #print('Tree updated', i, time)
+        nc = tree.collisionsTest()
+        #print('Collision Test done', nc)
+        ncol += nc
+        """if i % 100 == 0:
+            print(f'Total time {time:.1f}/{end_step} | iteration {i}/{Niter} | '
+                  f'collisions in last period: {ncol} | time step: {dt:.5f} | '
+                  f'kinetic energy: {ekin}')
+            ncol = 0
+            for ell in Ellipsoids:
+                speed = np.sqrt(ell.speedx ** 2 + ell.speedy ** 2 + ell.speedz ** 2)
+                if speed > 2.:
+                    print('Fast particle:', ell.id, speed, ell.x, ell.y, ell.z)
+                    # stop fast particle and move it closer to the center
+                    stop_part(ell)
+                    if ell.duplicate is not None:
+                        stop_part([ell.duplicate])"""
+
+        if dump and (i - 1) % ndump == 0:
             # Dump the ellipsoid information to be read by OVITO 
             # (Includes duplicates at periodic boundaries)
             write_dump(Ellipsoids, sim_box)
@@ -172,28 +210,22 @@ def particle_grow(sim_box, Ellipsoids, periodicity, nsteps,
         inter_ell = [ell for ell in Ellipsoids if not isinstance(ell.id, str)]
         Ellipsoids = inter_ell
 
-        if (not np.isclose(k_att, 0.0)) or (not np.isclose(k_rep, 0.0)):
-            calculateForce(Ellipsoids, sim_box, periodicity,
-                           k_rep=k_rep, k_att=k_att)
-
         dups = []
-        ekin = 0.
         # Loop over the ellipsoids: move, set Bbox, & 
         # check for wall collision / PBC
-        sc_fac = ((i + 1) / nsteps)**(1/3)  # scale factor for semi-axes during growth
+        sc_fac = (time / nsteps) ** (1 / 3)  # scale factor for semi-axes during growth
+        #print('Loop over ellipsoids:', i, time)
         for ellipsoid in Ellipsoids:
-            ekin += np.linalg.norm([ellipsoid.speedx, ellipsoid.speedy,
-                                    ellipsoid.speedz])
-            # Move the ellipsoid according to collision status
-            ellipsoid.move()
             # grow the ellipsoid
             ellipsoid.growth(sc_fac)
+            # Move the ellipsoid according to collision status
+            ellipsoid.move(dt)
             # Check for wall collision or create duplicates
             ell_dups = ellipsoid.wallCollision(sim_box, periodicity)
             dups.extend(ell_dups)
             # Update the BBox of the ellipsoid
             ellipsoid.set_cub()
-
+        #print('Done:', i, time)
         # Update the actual list with duplicates
         Ellipsoids.extend(dups)
 
@@ -228,10 +260,7 @@ def calculateForce(Ellipsoids, sim_box, periodicity, k_rep=0.0, k_att=0.0):
     d_half = sim_box.d / 2
 
     for ell in Ellipsoids:
-        ell.force_x = 0.
-        ell.force_y = 0.
-        ell.force_z = 0.
-        for ell_n in Ellipsoids:
+        for ell_n in ell.neighborlist:
             if ell.id != ell_n.id:
                 dx = ell.x - ell_n.x
                 dy = ell.y - ell_n.y
@@ -259,9 +288,9 @@ def calculateForce(Ellipsoids, sim_box, periodicity, k_rep=0.0, k_att=0.0):
 
                 # add repulsive or attractive force for dual phase systems
                 if ell.phasenum == ell_n.phasenum:
-                    Force = -k_rep * r2inv
+                    Force = k_rep * r2inv
                 else:
-                    Force = k_att * r2inv
+                    Force = -k_att * r2inv
 
                 ell.force_x += Force * dx / r
                 ell.force_y += Force * dy / r
@@ -296,38 +325,36 @@ def packingRoutine(particle_data, periodic, nsteps, sim_box,
     print('    Particle packing by growth simulation')
 
     particles, simbox = particle_grow(sim_box, Particles, periodic,
-                                      nsteps,
-                                      k_rep=k_rep, k_att=k_att, fill_factor=fill_factor,
-                                      dump=save_files)
+                                            nsteps,
+                                            k_rep=k_rep, k_att=k_att, fill_factor=fill_factor,
+                                            dump=save_files)
 
     # statistical evaluation of collisions
     if particles is not None:
         # check if particles are overlapping after growth
         ncoll = 0
-        ekin0 = 0.
         ekin = 0.
         for E1 in particles:
             E1.ncollision = 0
-            ekin0 += np.linalg.norm([E1.speedx0, E1.speedy0, E1.speedz0])
-            ekin += np.linalg.norm([E1.speedx, E1.speedy, E1.speedz])
+            ekin += E1.speedx ** 2 + E1.speedy ** 2 + E1.speedz ** 2
             for E2 in particles:
-                if E1.id != E2.id:
+                id1 = E1.id if E1.duplicate is None else (E1.duplicate + len(particles))
+                id2 = E2.id if E2.duplicate is None else (E2.duplicate + len(particles))
+                if id2 > id1:
                     # Distance between the centers of ellipsoids
                     dist = np.linalg.norm(np.subtract(E1.get_pos(),
                                                       E2.get_pos()))
                     # If the bounding spheres collide then check for collision
-                    if dist <= (E1.a + E2.a):
+                    if dist <= np.max([E1.a, E1.b, E1.c]) + np.max([E2.a, E2.b, E2.c]):
                         # Check if ellipsoids overlap and update their speeds
                         # accordingly
-                        if collide_detect(E1.get_coeffs(), E2.get_coeffs(),
-                                          E1.get_pos(), E2.get_pos(),
-                                          E1.rotation_matrix, E2.rotation_matrix):
+                        if collision_routine(E1, E2):
                             E1.ncollision += 1
+                            E2.ncollision += 1
                             ncoll += 1
         print('Completed particle packing')
         print(f'{ncoll} overlapping particles detected after packing')
         print(f'Kinetic energy of particles after packing: {ekin}')
-        print(f'Initial kinetic energy: {ekin0}')
         print('')
 
     return particles, simbox
diff --git a/src/kanapy/voxelization.py b/src/kanapy/voxelization.py
index ca1357ec..4845aac8 100644
--- a/src/kanapy/voxelization.py
+++ b/src/kanapy/voxelization.py
@@ -64,7 +64,6 @@ def assign_voxels_to_ellipsoid(cooDict, Ellipsoids, voxel_dict, vf_target=None):
     while vf_cur < vf_target:
         # call the growth value for the ellipsoids
         scale = next(growth)
-        print('Voxelization loop', scale, vf_cur, vf_target)
         for ellipsoid in Ellipsoids:
             # scale ellipsoid dimensions by the growth factor and generate surface points
             ellipsoid.a, ellipsoid.b, ellipsoid.c = scale * ellipsoid.a, scale * ellipsoid.b, scale * ellipsoid.c
@@ -282,35 +281,6 @@ def reassign_shared_voxels(cooDict, Ellipsoids, voxel_dict):
                             clo_ells[small_loc].inside_voxels.append(vox)
                         assigned_voxel.append(vox)
                         shared_voxels.remove(vox)
-
-            """ells = list(ells)                                     # convert to list
-            vox_coord = cooDict[vox]                              # Get the voxel position        
-            ells_pos = np.array([el.get_pos() for el in ells])    # Get the ellipsoids positions
-                    
-            # Distance b/w points along three axes
-            XDiff = vox_coord[0] - ells_pos[:, 0]  # 'x'-axis
-            YDiff = vox_coord[1] - ells_pos[:, 1]  # 'y'-axis
-            ZDiff = vox_coord[2] - ells_pos[:, 2]  # 'z'-axis
-    
-            # Find the distance from the 1st ellipsoid
-            dist = np.sqrt((XDiff**2)+(YDiff**2)+(ZDiff**2))
-       
-            clo_loc = np.where(dist == dist.min())[0]             # closest ellipsoid index        
-            clo_ells = [ells[loc] for loc in clo_loc]             # closest ellipsoids
-        
-            # If '1' closest ellipsoid: assign voxel to it        
-            if len(clo_ells) == 1:
-                clo_ells[0].inside_voxels.append(vox)
-            # Else: Determine the smallest and assign to it
-            else:
-                #clo_vol = np.array([ce.get_volume() for ce in clo_ells])    # Determine the volumes
-                #small_loc = np.where(clo_vol == clo_vol.min())[0]     # Smallest ellipsoid index
-                small_ells = [ells[loc] for loc in clo_loc]           # Smallest ellipsoids
-            
-                # assign to the smallest one regardless how many are of the same volume
-                small_ells[0].inside_voxels.append(vox)
-        
-            assigned_voxel.append(vox)"""
         ncyc += 1
     return
 
@@ -356,6 +326,7 @@ def ell2vox():
                 # sys.exit(0)
 
     def poly2vox():
+        assigned_vox = set()
         for pa in Ellipsoids:
             if pa.duplicate is not None:
                 gid = pa.duplicate
@@ -367,10 +338,23 @@ def poly2vox():
             # search for voxel centers inside the polygon
             for iv, ctr in mesh.vox_center_dict.items():
                 if pa.inner.find_simplex(ctr) >= 0:
-                    if type(iv) is not int:
-                        print(iv, ctr)
-                        iv = int(iv)
-                    vox.append(iv)
+                    if iv in assigned_vox:
+                        # voxel has already been assigned to different grain.
+                        # check nearest center and re-assign if necessary
+                        dst1 = np.linalg.norm(ctr - pa.get_pos())
+                        for igr, vlist in mesh.grain_dict.items():
+                            if iv in vlist:
+                                dst2 = np.linalg.norm(ctr - Ellipsoids[igr-1].get_pos())
+                                if dst1 < dst2:
+                                    vox.append(iv)
+                                    pa.inside_voxels.append(iv)
+                                    mesh.grain_dict[igr].remove(iv)
+                                    Ellipsoids[igr-1].inside_voxels.remove(iv)
+                                break
+                    else:
+                        assigned_vox.add(iv)
+                        vox.append(iv)
+                        pa.inside_voxels.append(iv)
             if gid in mesh.grain_dict.keys():
                 for iv in vox:
                     mesh.grain_dict[gid].append(iv)
diff --git a/tests/test_collide_detect_react.py b/tests/test_collide_detect_react.py
deleted file mode 100644
index 8dc0b051..00000000
--- a/tests/test_collide_detect_react.py
+++ /dev/null
@@ -1,293 +0,0 @@
-#!/usr/bin/env python
-# -*- coding: utf-8 -*-
-
-import pytest
-import numpy as np
-from unittest import mock
-
-import kanapy
-from kanapy.collisions import *
-from kanapy.entities import Ellipsoid
-
-
-@pytest.fixture
-def ellip(mocker):
-
-    # Define attributes to mocker object
-    a, b, c = 2.0, 1.5, 1.5                                     # Semi axes lengths    
-    vec_a = np.array([a*np.cos(90), a*np.sin(90), 0.0])         # Tilt vector wrt (+ve) x axis    
-    cross_a = np.cross(np.array([1, 0, 0]), vec_a)              # Find the quaternion axis    
-    norm_cross_a = np.linalg.norm(cross_a, 2)                   # norm of the vector (Magnitude)    
-    quat_axis = cross_a/norm_cross_a                            # normalize the quaternion axis
-
-    # Find the quaternion components
-    qx, qy, qz = quat_axis * np.sin(90/2)
-    qw = np.cos(90/2)
-    quat = np.array([qw, qx, qy, qz])
-
-    # Generate rotation matrix
-    Nq = qw*qw + qx*qx + qy*qy + qz*qz
-
-    s = 2.0/Nq
-    X = qx*s
-    Y = qy*s
-    Z = qz*s
-    wX = qw*X
-    wY = qw*Y
-    wZ = qw*Z
-    xX = qx*X
-    xY = qx*Y
-    xZ = qx*Z
-    yY = qy*Y
-    yZ = qy*Z
-    zZ = qz*Z
-
-    rotation_matrix = np.array([[1.0-(yY+zZ), xY-wZ, xZ+wY],
-                                [xY+wZ, 1.0-(xX+zZ), yZ-wX],
-                                [xZ-wY, yZ+wX, 1.0-(xX+yY)]])
-
-    # Rotation matrix has to be transposed as OVITO uses the transposed matrix for visualization.
-    rotation_matrix = rotation_matrix.T
-
-    # Points on the outer surface of Ellipsoid
-    u = np.linspace(0, 2 * np.pi, 20)
-    v = np.linspace(0, np.pi, 20)
-
-    # Cartesian coordinates that correspond to the spherical angles:
-    xval = a * np.outer(np.cos(u), np.sin(v))
-    yval = b * np.outer(np.sin(u), np.sin(v))
-    zval = c * np.outer(np.ones_like(u), np.cos(v))
-
-    # combine the three 2D arrays element wise
-    stacked_xyz = np.stack((xval.ravel(), yval.ravel(), zval.ravel()), axis=1)
-
-    # Define the mocker objects
-    el1 = mocker.MagicMock()
-    el2 = mocker.MagicMock()
-
-    # assign attributes to the mocker objects
-    el1.rotation_matrix = rotation_matrix
-    el1.surface_points = stacked_xyz.dot(rotation_matrix)
-    el1.id = 1
-    el1.a, el1.b, el1.c = 2.0, 1.5, 1.5
-    el1.x, el1.y, el1.z = 1, 0.5, 0.75
-    el1.speedx, el1.speedy, el1.speedz = -0.02, 0.075, -0.05
-    el1.quat = quat
-    el1.inside_voxels = []
-    el1.get_pos.return_value = np.array([1, 0.5, 0.75])
-    el1.get_coeffs.return_value = np.array([2.0, 1.5, 1.5])
-    el1.get_volume.return_value = (4/3)*np.pi*a*b*c
-
-    el2.rotation_matrix = rotation_matrix
-    el2.surface_points = stacked_xyz.dot(rotation_matrix)
-    el2.id = 2
-    el2.a, el2.b, el2.c = 2.0, 1.5, 1.5
-    el2.x, el2.y, el2.z = 1.9, 1.68, 2.6
-    el2.speedx, el2.speedy, el2.speedz = 0.5, -0.025, -0.36
-    el2.quat = quat
-    el2.inside_voxels = []
-    el2.get_pos.return_value = np.array([1.9, 1.68, 2.6])
-    el2.get_coeffs.return_value = np.array([2.0, 1.5, 1.5])
-    el2.get_volume.return_value = (4/3)*np.pi*a*b*c
-
-    return [el1, el2]
-
-
-def test_collideDetect(ellip):
-
-    status = collide_detect(ellip[0].get_coeffs(), ellip[1].get_coeffs(),
-                           ellip[0].get_pos(), ellip[1].get_pos(),
-                           ellip[0].rotation_matrix, ellip[1].rotation_matrix)
-
-    assert status == True
-
-
-def test_collision_react():
-
-    # Initialize the Ellipsoids
-    el1 = Ellipsoid(1, 5, 5, 5, 2.0, 2.0, 2.0, np.array(
-        [0.52532199, 0., -0., 0.85090352]))
-    el2 = Ellipsoid(2, 5.5, 5.5, 5.5, 2.0, 2.0, 2.0,
-                    np.array([0.52532199, 0., -0., 0.85090352]))
-
-    # el1.speedx, el1.speedy, el1.speedz = 0.1, 0.075, -0.05
-    # el2.speedx, el2.speedy, el2.speedz = -0.1, -0.025, -0.36
-    
-    el1.speedx0, el1.speedy0, el1.speedz0 = 0.1, 0.075, -0.05
-    el2.speedx0, el2.speedy0, el2.speedz0 = -0.1, -0.025, -0.36
-        
-    # Test different conditions
-    # Condition: xdiff > 0 && zdiff > 0          
-    collision_react(el1, el2)
-    collision_react(el2, el1)
-    
-    assert round(el1.speedx, 6) == -0.077728
-    assert round(el1.speedy, 6) == -0.077728
-    assert round(el1.speedz, 6) == -0.077728
-    
-    assert round(el2.speedx, 6) == 0.216198
-    assert round(el2.speedy, 6) == 0.216198
-    assert round(el2.speedz, 6) == 0.216198    
-
-    # Condition: xdiff > 0 && zdiff < 0
-    el2.z = 4.5 
-    el1.speedx, el1.speedy, el1.speedz = 0.0, 0.0, 0.0
-    el2.speedx, el2.speedy, el2.speedz = 0.0, 0.0, 0.0
-    
-    collision_react(el1, el2)
-    collision_react(el2, el1)   
-
-    assert round(el1.speedx, 6) == -0.077728
-    assert round(el1.speedy, 6) == -0.077728
-    assert round(el1.speedz, 6) == 0.077728
-    
-    assert round(el2.speedx, 6) == 0.216198
-    assert round(el2.speedy, 6) == 0.216198
-    assert round(el2.speedz, 6) == -0.216198          
-
-    # Condition: xdiff < 0 && zdiff > 0
-    el2.x = 4.5
-    el2.z  = 5.5  
-    el1.speedx, el1.speedy, el1.speedz = 0.0, 0.0, 0.0
-    el2.speedx, el2.speedy, el2.speedz = 0.0, 0.0, 0.0
-
-    collision_react(el1, el2)
-    collision_react(el2, el1)   
-
-    assert round(el1.speedx, 6) == 0.077728
-    assert round(el1.speedy, 6) == -0.077728
-    assert round(el1.speedz, 6) == -0.077728
-    
-    assert round(el2.speedx, 6) == -0.216198
-    assert round(el2.speedy, 6) == 0.216198
-    assert round(el2.speedz, 6) == 0.216198          
-
-    # Condition: xdiff < 0 && zdiff < 0
-    el2.x = 4.5
-    el2.z  = 4.5  
-    el1.speedx, el1.speedy, el1.speedz = 0.0, 0.0, 0.0
-    el2.speedx, el2.speedy, el2.speedz = 0.0, 0.0, 0.0
-
-    collision_react(el1, el2)
-    collision_react(el2, el1)   
-
-    assert round(el1.speedx, 6) == 0.077728
-    assert round(el1.speedy, 6) == -0.077728
-    assert round(el1.speedz, 6) == 0.077728
-    
-    assert round(el2.speedx, 6) == -0.216198
-    assert round(el2.speedy, 6) == 0.216198
-    assert round(el2.speedz, 6) == -0.216198          
-
-    # Condition: xdiff = 0 && zdiff = 0
-    el2.x = 5.0
-    el2.z = 5.0  
-    el1.speedx, el1.speedy, el1.speedz = 0.0, 0.0, 0.0
-    el2.speedx, el2.speedy, el2.speedz = 0.0, 0.0, 0.0
-
-    collision_react(el1, el2)
-    collision_react(el2, el1)   
-
-    assert round(el1.speedx, 6) == 0.0
-    assert round(el1.speedy, 6) == -0.134629
-    assert round(el1.speedz, 6) == 0.0
-    
-    assert round(el2.speedx, 6) == 0.0
-    assert round(el2.speedy, 6) == 0.374466
-    assert round(el2.speedz, 6) == 0.0          
-
-    # Condition: xdiff = 0 && zdiff > 0
-    el2.x = 5.0
-    el2.z  = 5.5  
-    el1.speedx, el1.speedy, el1.speedz = 0.0, 0.0, 0.0
-    el2.speedx, el2.speedy, el2.speedz = 0.0, 0.0, 0.0
-
-    collision_react(el1, el2)
-    collision_react(el2, el1)   
-
-    assert round(el1.speedx, 6) == -0.0
-    assert round(el1.speedy, 6) == -0.095197
-    assert round(el1.speedz, 6) == -0.095197
-    
-    assert round(el2.speedx, 6) == 0.0
-    assert round(el2.speedy, 6) == 0.264788
-    assert round(el2.speedz, 6) == 0.264788          
-
-    # Condition: xdiff = 0 && zdiff < 0    
-    el2.x = 5.0
-    el2.z  = 4.5  
-    el1.speedx, el1.speedy, el1.speedz = 0.0, 0.0, 0.0
-    el2.speedx, el2.speedy, el2.speedz = 0.0, 0.0, 0.0
-
-    collision_react(el1, el2)
-    collision_react(el2, el1)   
-
-    assert round(el1.speedx, 6) == -0.0
-    assert round(el1.speedy, 6) == -0.095197
-    assert round(el1.speedz, 6) == 0.095197
-    
-    assert round(el2.speedx, 6) == 0.0
-    assert round(el2.speedy, 6) == 0.264788
-    assert round(el2.speedz, 6) == -0.264788          
-
-    # Condition: xdiff < 0 && zdiff = 0    
-    el2.x = 4.5
-    el2.z  = 5.0  
-    el1.speedx, el1.speedy, el1.speedz = 0.0, 0.0, 0.0
-    el2.speedx, el2.speedy, el2.speedz = 0.0, 0.0, 0.0
-
-    collision_react(el1, el2)
-    collision_react(el2, el1)   
-    
-    assert round(el1.speedx, 6) == 0.095197
-    assert round(el1.speedy, 6) == -0.095197
-    assert round(el1.speedz, 6) == 0.0
-    
-    assert round(el2.speedx, 6) == -0.264788
-    assert round(el2.speedy, 6) == 0.264788
-    assert round(el2.speedz, 6) == 0.0          
-
-
-    # Condition: xdiff > 0 && zdiff = 0    
-    el2.x = 5.5
-    el2.z  = 5.0  
-    el1.speedx, el1.speedy, el1.speedz = 0.0, 0.0, 0.0
-    el2.speedx, el2.speedy, el2.speedz = 0.0, 0.0, 0.0
-
-    collision_react(el1, el2)
-    collision_react(el2, el1)   
-
-    assert round(el1.speedx, 6) == -0.095197
-    assert round(el1.speedy, 6) == -0.095197
-    assert round(el1.speedz, 6) == 0.0
-    
-    assert round(el2.speedx, 6) == 0.264788
-    assert round(el2.speedy, 6) == 0.264788
-    assert round(el2.speedz, 6) == 0.0          
-                            
-
-def test_collision_routine_sphere():
-
-    # Initialize the Ellipsoids
-    el1 = Ellipsoid(1, 5, 5, 5, 2.0, 2.0, 2.0, np.array(
-        [0.52532199, 0., -0., 0.85090352]))
-    el2 = Ellipsoid(2, 5.5, 5.5, 5.5, 2.0, 2.0, 2.0,
-                    np.array([0.52532199, 0., -0., 0.85090352]))
-
-    el1.speedx, el1.speedy, el1.speedz = 0.1, 0.075, -0.05
-    el2.speedx, el2.speedy, el2.speedz = -0.1, -0.025, -0.36
-
-    # Test if the 'collision_react' function is called twice
-    with mock.patch('kanapy.collisions.collision_react') as mocked_method:
-        collision_routine(el1, el2)
-        assert mocked_method.call_count == 2
-
-
-def test_collision_routine_ellipsoid(ellip):
-
-    # Test if the 'collision_react' function is called twice
-    with mock.patch('kanapy.collisions.collision_react') as mocked_method:
-        collision_routine(ellip[0], ellip[1])
-        assert mocked_method.call_count == 2
-        
-        
diff --git a/tests/test_entities.py b/tests/test_entities.py
index 302e8beb..419b433f 100644
--- a/tests/test_entities.py
+++ b/tests/test_entities.py
@@ -209,13 +209,13 @@ def test_move(self, rot_surf):
 
         # Initialize the Ellipsoid
         self.ell = Ellipsoid(1, 1, 0.5, 0.75, 2.0, 1.5, 1.5, rot_surf[0])
-        self.ell.speedx, self.ell.speedy, self.ell.speedz = -0.02, 0.075, -0.05
+        self.ell.force_x, self.ell.force_y, self.ell.force_z = -0.02, 0.075, -0.05
 
-        self.ell.move()
+        self.ell.move(1.0)
 
-        assert self.ell.x == 1 + (-0.02)
-        assert self.ell.y == 0.5 + (0.075)
-        assert self.ell.z == 0.75 + (-0.05)
+        assert np.isclose(self.ell.x, 0.98)
+        assert np.isclose(self.ell.y, 0.575)
+        assert np.isclose(self.ell.z, 0.7)
 
     def test_gravity(self, rot_surf):
 
@@ -257,20 +257,19 @@ def test_wallCollision_Periodic(self, rot_surf, position, duplicates):
         assert len(dups) == duplicates
                                                           
                 
-    @pytest.mark.parametrize("position, speeds", [(np.array([0.2, 0.5, 0.75]), [0.02, -0.075, 0.05]), 
-                                                  (np.array([9.8, 9.5, 9.25]), [0.02, -0.075, 0.05])])
-    def test_wallCollision_Regular(self, rot_surf, position, speeds):
+    @pytest.mark.parametrize("position", [(np.array([0.1, 5.0, 0.75]))])
+    def test_wallCollision_Regular(self, rot_surf, position):
 
         sbox = Simulation_Box((10, 10, 10))
         self.ell = Ellipsoid(1, position[0], position[1], position[2], 2.0, 1.5, 1.5, rot_surf[0])
-        self.ell.speedx, self.ell.speedy, self.ell.speedz = -0.02, 0.075, -0.05
 
         # Test for periodicity == FALSE
-        self.ell.wallCollision(sbox, False)
-        assert self.ell.speedx == speeds[0]
-        assert self.ell.speedy == speeds[1]
-        assert self.ell.speedz == speeds[2]
-
+        dup = self.ell.wallCollision(sbox, False)
+        assert dup == []
+        assert np.isclose(self.ell.x, 1.6072736670942065)
+        assert np.isclose(self.ell.y, 5.0)
+        assert np.isclose(self.ell.z, 1.5)
+        assert np.isclose(self.ell.xold, self.ell.x)
         assert isinstance(self.ell.get_cub(), Cuboid)
 
 
@@ -323,7 +322,7 @@ def test_subdivide_particles(self, rot_surf):
 
     def test_collisions(self, rot_surf):
 
-        cbox = Cuboid(0, 0, 10, 10, 0, 10)              # Initialize the Cuboid
+        cbox = Cuboid(0, 0, 10, 10, 0, 10)  # Initialize the Cuboid
 
         # Initialize the Ellipsoids
         ell1 = Ellipsoid(1, 1, 0.5, 0.75, 2.0, 1.5, 1.5, rot_surf[0])
@@ -335,15 +334,14 @@ def test_collisions(self, rot_surf):
         self.tree = Octree(0, cbox, particles=[ell1, ell2])
         self.tree.update()
         #self.tree.subdivide_particles()
-        self.tree.collisionsTest()
-
-        assert round(ell1.speedx, 6) == -0.035037
-        assert round(ell1.speedy, 6) == -0.045938
-        assert round(ell1.speedz, 6) == -0.072021
-
-        assert round(ell2.speedx, 6) == 0.233994
-        assert round(ell2.speedy, 6) == 0.306793
-        assert round(ell2.speedz, 6) == 0.480988
+        nc = self.tree.collisionsTest()
+
+        assert nc == 1
+        assert np.isclose(ell1.force_x, -0.5361842103706771)
+        assert np.isclose(ell2.force_x, 0.5361842103706771)
+        assert np.isclose(ell1.force_y + ell2.force_y, 0.0)
+        assert np.isclose(ell1.force_z, -1.1021564324286142)
+        assert np.isclose(ell2.force_z, 1.1021564324286142)
 
     def test_update1(self, mocker, rot_surf):
 
diff --git a/tests/test_packing.py b/tests/test_packing.py
index 8ebcd9ac..0c53642f 100644
--- a/tests/test_packing.py
+++ b/tests/test_packing.py
@@ -104,11 +104,12 @@ def test_particle_grow(rot_surf):
     ells = [ell1, ell2]
     sb = Simulation_Box((10, 10, 10))
 
-    particle_grow(sb, ells, True, 10, dump=True)
-    
-    assert np.isclose(ell1.x, 0.9)
-    assert np.isclose(ell1.y, 0.875)
-    assert np.isclose(ell2.z, 0.8)
+    particles, simbox = particle_grow(sb, ells, True, 10, dump=True)
+
+    assert np.isclose(ell1.a, 1.5874010519681996)
+    assert np.isclose(ell1.oria, 2.0)
+    assert np.isclose(ell2.get_volume(), 9.424777960769381)
+    assert particles[0] == ell1
 
     
 def test_packingRoutine():