docs: def get_D_square_min(atoms_now, atoms_old, i_now, j_now, delta_plus_epsilon=None):
      """
      Calculate the D^2_min norm of Falk and Langer
      """
      nat = len(atoms_now)
      assert len(atoms_now) == len(atoms_old)
  
      pos_now = atoms_now.positions
      pos_old = atoms_old.positions
  
      # Compute current and old distance vectors. Note that current distance
      # vectors cannot be taken from the neighbor calculation, because neighbors
      # are calculated from the sheared cell while these distance need to come
      # from the unsheared cell. Taking the distance from the unsheared cell
      # make periodic boundary conditions (and flipping of cell) a lot easier.
      dr_now = mic(pos_now[i_now] - pos_now[j_now], atoms_now.cell)
      dr_old = mic(pos_old[i_now] - pos_old[j_now], atoms_old.cell)
  
      # Sanity check: Shape needs to be identical!
      assert dr_now.shape == dr_old.shape
  
      if delta_plus_epsilon is None:
          # Get minimum strain tensor
          delta_plus_epsilon = get_delta_plus_epsilon(nat, i_now, dr_now, dr_old)
  
      # Spread epsilon out for each neighbor index
      delta_plus_epsilon_n = delta_plus_epsilon[i_now]
  
      # Compute D^2_min
      d_sq_n = np.sum(
          (
          dr_now-
          np.sum(delta_plus_epsilon_n.reshape(-1,3,3)*dr_old.reshape(-1,1,3),
                 axis=2)
          )**2,
          axis=1)
  
      # For each atom, sum over all neighbors
      d_sq = np.bincount(i_now, weights=d_sq_n)
  
      return delta_plus_epsilon, d_sq