VASP 6.5.1: recurrent "internal error in rot.F / EDWAV: the gradient is not orthogonal"

[andres_bort1]

Dear VASP developers,

I repeatedly hit the internal error below in a 287–297-atom DFT+U perovskite supercell (VASP 6.5.1). The error: internal error in: rot.F (...) EDWAV: internal error, the gradient is not orthogonal (...). It has killed 32 jobs so far (about 30% of cases) of SQS generated bulk supercells. It occurs either after 1 to 56 successful ionic steps. A rerun from the crashed job's CONTCAR with identical inputs gets past the crash point most of the times. It happens both with cold starts (ISTART=0, ICHARG=2: 13 crashes) or with warm-started relaxations (ISTART=1, ICHARG=0: 19 crashes), both with Gamma-point calculations (but with vasp_std bin - I have not tried with vasp_gam) and with 4 irreducible KPOINTS.

I have encountered a (maybe) related anomaly in some of the jobs that crashed: a wrong electron count on the FINAL SCF iteration of an ionic step. In relaxations of this cell I see occasional number of electron lines off NELECT by 0.003–0.12 e. Every instance checked is on the last SCF iteration of an ionic step, i.e. the iteration whose density gives the forces. But some of the cases that showed this anomaly eventually do converged after restoring the proper electron count.

What I have tried (None of these removed the crash):
NELMDL (0 and −5); three mixing/algorithm combinations (ALGO=Normal and ALGO=All, default and conservative AMIX/BMIX/AMIX_MAG BMIX_MAG); LREAL=.FALSE. (failed earlier); LDAU=.FALSE.; ISPIN=1. I also ruled out NBANDS exhaustion (ample headroom).

I have not tried: IWAVPR changes, ALGO=Damped/Fast for the relaxations, ADDGRID, LMAXMIX=6, or a different MPI/compiler build. I would rather ask you before changing anything that could affect results across an ongoing campaign.

System and settings (identical in all crashed jobs)

• Sr0.95Ti0.3Fe0.6Cu0.1O3−x perovskite, 3×4×5 supercell (60 B sites: 18 Ti, 36 Fe, 6 Cu; 57 Sr; 170–180 O; 287–297 atoms);
  fixed lattice about 11.76 × 15.72 × 19.54 Å .
• PBE+U, Dudarev (LDAUTYPE=2): U(Ti)=4.91, U(Fe)=4.5, U(Cu)=6.47 eV. Collinear AFM A-type order.
• POTCARs: PAW_PBE O, Ti_sv, Fe, Cu, Sr_sv.
• ENCUT=600, PREC=Accurate, LREAL=Auto, ISYM=0, ISMEAR=0, SIGMA=0.05, NBANDS=1344, ALGO=All (IALGO 58), ISEARCH=1,
  AMIX=0.05, BMIX=0.01, AMIX_MAG=0.1, BMIX_MAG=0.01, LASPH=.TRUE., LMAXMIX=4, EDIFF=1E-6, NELM=900,
  IBRION=2 (POTIM 0.5), EDIFFG=−0.03.
• Parallelisation: NCORE=8, 64 MPI ranks per node, 1 thread per rank, KPAR=1 (Γ) or 2/4 (k-mesh).

Build:

• vasp.6.5.1.
• Compiler and libraries:
  • GNU Fortran 11.5.0 (Red Hat 11.5.0-2) via Intel MPI's mpif90;
  • -O2 -march=native;
  • MKL 2021.4.0 (gf_lp64, sequential, ScaLAPACK, BLACS-IntelMPI);
  • CPP options: -DMPI -DMPI_BLOCK=8000 -Duse_collective -DscaLAPACK -DCACHE_SIZE=4000 -Davoidalloc -Dvasp6 -Dtbdyn
    -Dfock_dblbuf -DVASP_HDF5 -DVASP2WANNIER90.
  • makefile.include is attached.
  • The binary was built against Intel MPI 2021.16 (library path embedded in the binary).
  • At run time the environment loads oneapi/mpi/latest, which has pointed to Intel MPI 2021.18.1
• The binary contains no AVX-512 instructions. It was built with -march=native on an AMD EPYC 7262

I attach two bundles with complete inputs and outputs for 7 crashed jobs and 2 successful comparison jobs, the full crash list, and the build and hardware details.

Questions:

1. Is this a known issue in 6.5.1 for ALGO=All with large spin-polarised DFT+U cells? Is there a fix or patch? Is the occasional wrong electron count on the final SCF iteration a known symptom of the same problem?
2. In the relaxation case, is there a recommended setting that avoids the failure without changing the converged result? For example IWAVPR, ALGO=Damped, or a smaller trial step.
3. Could the build/run MPI mismatch (built 2021.16, run 2021.18.1) or -march=native built on Zen 2 and run on Zen 4 plausibly cause this? Should we rebuild?

Best regards,

[MS edit: changed a few lines that were misinterpreted as LaTeX]

[martin.schlipf]

Thank you for this bug report, I will investigate.

[martin.schlipf]

This is indeed a recurring issue but it was not on our list of known issues yet. I have added it now. Within you fill find a link to this thread that proposes a solution. The fix has been merged into our working branch and will be shipped with the next version. You can add the following lines before the first EDWAV call in rot.F

Code: Select all

IF (ABS(INFO%NELMDL) > 0 .AND. N == ABS(INFO%NELMDL)+1) THEN
   CALL ORTHCH(WDES,W, INFO%LOVERL, LMDIM,CQIJ)
ENDIF 

There may be a second issue that causes the warm-relaxation crashes. You mentioned that ALGO = Normal did result in the same crash but that should not be possible because it does not travel the same path in the code. Perhaps it crashes too, but with a different error? Anyway, it is not recommended for complex magnetic structures because it often struggles to converge. To address this second issue you could try

• You could try to change IWAVPR that would alter the orbitals in the next ionic step.
• Switching off scalapack LSCALAPACK = .FALSE. (no recompilation needed), if that helps we need to look into the SCALAPACK calls.

Perhaps the fix above is to restricted. It will only apply to the first ionic step. You could try the more general

Code: Select all

IF (N == ABS(INFO%NELMDL)+1) THEN
   CALL ORTHCH(WDES,W, INFO%LOVERL, LMDIM,CQIJ)
ENDIF

That would orthonormalize before every ionic step. Please report back if that helps in your case.

[andres_bort1]

Dear Martin,

I really appreciate your quick response. And many thanks again for the patch suggestion. We tested the general variant, which orthonormalises before every ionic step, and it fixed the crash for us.

What we did: using ALGO = All, Γ-only with vasp_std (which we chose for the relaxations given the size of the cells and our computational resources), IBRION = 2, each warm-started from a converged WAVECAR. As N, INFO and CQIJ are not in scope in rot.F, we put the call in electron_all.F (ELMIN_ALL), just before the first CALL EDWAV:

IF (N == ABS(INFO%NELMDL)+1) THEN
CALL ORTHCH(WDES,W, INFO%LOVERL, LMDIM,CQIJ)
ENDIF

We also printed a marker line to OUTCAR, which confirmed that the call ran on every ionic step.

We reran 8 relaxations that had all crashed with rot.F:822 / EDWAV on the stock binary.
- Inputs and seed WAVECAR were byte-identical; only the binary changed.
- On the stock binary these runs crashed after 1 to 56 ionic steps.
- Across all our warm relaxations, the stock binary crashed 19 times in 1498 ionic steps, about 1.3 % per step.

Result with the patch:
- 0 rot.F/EDWAV crashes in 404 ionic steps.
- All 8 relaxations converged (EDIFFG = −0.03 eV/Å).

We confirmed that the patch does not change results:
- A warm static with identical inputs and seed gave TOTENs 0.003 meV apart. Per-ion moments agreed to within 0.003 μB.
- For the 4 relaxations whose stock reruns had also converged, final energies agree to within 2 meV (−1.1 to +1.9 meV), and all Fe moment signs match.

Would you recommend this general variant for production use, or is there a reason to prefer the first-ionic-step-only version? We're happy to share inputs if they're useful.

Best,

[martin.schlipf]

If this fixes the issue for you, you can use this variant for production use. The solution for the first iteration, will be released with the next version of VASP.

The main reason not to use it, is that it is often the signature of a different underlying issue (e.g. atoms in incorrect positions). By circumventing the early crash, you may still end up with convergence issues later on but then you have already used up more compute time.