ELPH_MODE = TRANSPORT: "Sub-Space-Matrix is not hermitian in DAV" vasp 6.5.1

[ghaithjamal_nadhreen]

Dear VASP developers,

I am running an electron-phonon transport study (MRTA) on with VASP 6.5.1 for 6 Re-containing compounds (CrReAs, MoReAs, MoReBi, MoReSb, NbReTe, WReAs), the run crashes deterministically at the same point, and I have run out of input-side workarounds. I would be grateful for advice.

## System and workflow

- Half-Heusler XYZ, F-43m, 3-atom primitive cell, ISYM = 0
- Example: MoReAs, PAW_PBE Mo_sv / Re / As, NELECT = 26, ENCUT = 520 eV, PREC = Accurate, LREAL = .FALSE.
- Electron-phonon potential from phelel (phelel_params.hdf5), then 'ELPH_MODE = TRANSPORT', 'ELPH_SCATTERING_APPROX = MRTA_TAU', 'ELPH_ISMEAR = -15', KPOINTS_ELPH 12x12x12 Γ-centred, 107 carrier densities, 13 temperatures
- 1 node, 192 MPI ranks (AMD Genoa, Cray MPICH, Intel classic 2023.1 build), OMP_NUM_THREADS = 1
- Input files are attached.

## What happens

1. The SCF converges normally.
2. The KPOINTS_OPT loop over the 1728 KPOINTS_ELPH points completes without any warning (with 24 k-point groups, 'k-point batch [ 36–37 / 868 ]' is the last batch).
3. Right after the last 'k-point batch' line, before anything else is printed, the run emits ~10^6 lines of

Code: Select all

   WARNING: Sub-Space-Matrix is not hermitian in DAV            1
    -5.735061726052421E-004

and the printed deviation grows up to ~7e2. Then one group of ranks calls 'MPI_Abort(MPI_COMM_WORLD, 1)', within ~1.5 minutes of the start. No other error message is printed.

In runs that succeed, the line after the last 'k-point batch' is 'direction 1', 'direction 2', 'direction 3', followed by the interatomic force constants and the q-point loop. So the failure happens in the step that prints 'direction 1..3'. OUTCAR shows 'elph_velocity_mode = 2' and 'elph_lr = 1', so I suspect the band-velocity (linear-response / commutator) step, but I cannot confirm this from the output.

## What I have tested (all on MoReAs unless stated otherwise; each changes one thing)

| Change | Result |
|---|---|
| Band window ELPH_SELFEN_BAND_START/STOP (several windows, chosen from the band structure, edges not splitting degenerate bands); NBANDS 24 vs 48 | fails (also on VReTe in its failing layout) |
| ALGO = All (only affects the SCF) | fails (VReTe) |
| VASP 6.6.1 | fails (VReTe, in the layout that fails there) |
| EDIFF 1E-8 → 1E-6 | fails (NbReTe at the same point; MoReAs failed earlier, in the SCF, in that layout, see note) |
| Parallelisation: KPAR 24 / NCORE 1 (1 core per k-point in el-ph), KPAR 1 / NCORE 8, KPAR 24 / NCORE 8 | fails at the same point in all three for NbReTe; for MoReAs in KPAR 24 / NCORE 1 and KPAR 24 / NCORE 8 (KPAR 1 / NCORE 8 fails earlier, in the SCF, see note) |
| ELPH_WF_COMM_OPT = 1 (known issue 68 workaround) | fails, exactly the same number of warning lines (1,034,080) as without it |
| ELPH_WF_REDISTRIBUTE = .FALSE. | fails, again exactly the same number of warning lines |
| KPOINTS_ELPH 4x4x4, 6x6x6, 8x8x8 | passes (velocity step, q-point loop, transport output all fine) |
| KPOINTS_ELPH 11x11x11, 12x12x12, 13x13x13 | fails at the same point |
| KPOINTS_ELPH 12x12x12 with shift 0.5 0.5 0.5 | run ends after ~36 s (output not kept) |

The number of warning lines is identical across the parallelisation and communication variants with the same mesh and NBANDS. The failure looks deterministic and numerical rather than an MPI issue.

Note on parallelisation: with KPAR 1 / NCORE 8, the SCF itself diverges for CrReAs, MoReAs, MoReBi and WReAs ("not hermitian in DAV", energy → −1e9 eV after nearly converging), while KPAR 24 / NCORE 8 and KPAR 24 / NCORE 1 converge the SCF normally. This may be a separate issue.

## Materials that do work

- VReTe and CrReSb fail with KPAR 24 / NCORE 1 but pass at 12x12x12 with KPAR 1 / NCORE 8.
- CrReBi, TaReTe, WReBi and WReSb pass with the original KPAR 24 / NCORE 1 settings.

## Questions

1. Is this a known problem in the velocity step of the electron-phonon transport driver (elph_velocity_mode = 2)?
2. Is there an INCAR setting to change how the band velocities are computed for ELPH_TRANSPORT (e.g. a different velocity mode), or to make this step more robust?
3. Is there anything else we should try, short of using a coarser mesh? We would like to keep 12x12x12 for consistency with the rest of the dataset.

I can provide full OUTCAR, std.out, std.err and phelel_params.hdf5 if needed.

Thank you very much,
Ghaith Nadhreen
KAUST

[bernhard_schmiedmayer]

Dear Ghaith Nadhreen,

thank you for bringing this to our attention and for providing the input files. This made it easy to reproduce the issue.

The crash is not related to the electron-phonon transport driver, as it occurs independently of it. I was able to reproduce the failure, using your POSCAR, INCAR and POTCAR files without any ELPH tags. Your calculation fails during the blocked Davidson diagonalisation of the dense KPOINTS_ELPH mesh. The velocity step is only reached afterwards.

VASP does detect the breakdown and raises a fatal error: "EDDAV: Call to ZHEGV failed. Returncode = 10 2 16". Unfortunately, you never saw this error message because it was raised in only one k-point group (MPI rank), and our error routine only writes on the master rank. The job is then terminated with MPI_Abort, so the message is lost.

The calculation most likely crashes because EDIFF is used as the stopping criterion for two different energies: First in the SCF loop, and then on the eigenvalues of the bands. The latter increases linearly with the number of k-points. As EDIFF is not updated, this criterion becomes increasingly strict as the number of k-points increases, eventually becoming unreachable. This is the point at which the DAV optimizer fails. However, I have not yet found a reliable workaround (loosening EDIFF further is not a reliable fix), so there may be more to it than that. I will continue to look into it and let you know when I have found a solution.

Best,
Bernhard