code: def _get_I(self, a, b, size, plus_transpose=True):
          
          r_sum = np.zeros((3, 3), dtype='double', order='C')
          for r in self._rotations_cartesian:
              for i in range(3):
                  for j in range(3):
                      r_sum[i, j] += r[a, i] * r[b, j]
          if plus_transpose:
              r_sum += r_sum.T
  
          # Return None not to consume computer for diagonalization
          if (np.abs(r_sum) < 1e-10).all():
              return None
  
          # Same as np.kron(np.eye(size), r_sum), but writen as below
          # to be sure the values in memory C-congiguous with 'double'.
          I_mat = np.zeros((3 * size, 3 * size), dtype='double', order='C')
          for i in range(size):
              I_mat[(i * 3):((i + 1) * 3), (i * 3):((i + 1) * 3)] = r_sum
  
          return I_mat
docs: Return I matrix in Chaput's PRL paper.
  
          None is returned if I is zero matrix.