docs: def _symm_herm(C):
      """To get rid of NaNs produced by _scalar2array, symmetrize operators
      where C_ijkl = C_jilk*"""
      nans = np.isnan(C)
      C[nans] = np.einsum('jilk', C)[nans].conj()
      return C