docs: def _symm_current(C):
      """To get rid of NaNs produced by _scalar2array, symmetrize operators
      where C_ijkl = C_klij"""
      nans = np.isnan(C)
      C[nans] = np.einsum('klij', C)[nans]
      return C