Electronic Convergence Issues with Nonlocal vdW-DF Functionals (in VASP)

Question

Compared to say the DFT-D3 type of dispersion corrections (set with IVDW in VASP), I've found the nonlocal vdW-DF functionals (described here in the VASP Wiki) to be much more difficult to converge to similar levels of accuracy (which I guess should be expected given the nature of the vdW correction). My monolayer calculations with the optB88-vdW functional led to some geometry optimizations with several consecutive unconverged electronic optimization steps (even with NELM = 90, i.e. 90 SCF steps), when starting from a converged structure from a database.

My first impression was that my use of EDIFF=1E-07 (VASP's electronic convergence parameter) and EDIFFG = -0.002 (VASP's geometry convergence parameter) was too fine, which seemed somewhat in line with what I thought was a rougher convergence criterion by 2DMatPedia (a 2D materials database which used the optB88-vdW functional and an ionic convergence criterion that the energy difference is below $10^{-4}$ eV for structure relaxations). However, the JARVIS database, which also contains 2D monolayers, seemed to be able to do calculations with optB88-vdW and (what I thought were) stricter convergence criteria:

"Both the internal atomic positions and the lattice constants are allowed to relax in spin-unrestricted calculations until the maximal residual Hellmann–Feynman forces on atoms are smaller than $0.001$ eV Å$^{−1}$ and energy-tolerance of $10^{-7}$ eV with accurate precision setting (PREC=Accurate)."

This leads me to my question.

Are the convergence difficulties I've faced with optB88-vdW expected (and so it is reasonable to relax my convergence criteria)? Or are there other things that I can tweak to achieve the tighter convergence criteria that would be important for later phonon and elastic constants calculations for example?

I have played around with just trying to converge static calculations tightly (EDIFF=1E-07), by changing the electronic minimization algorithm to ALGO = Normal or ALGO = ALL (I usually use ALGO = Fast), or playing around with IMIX (density mixing). However, I could not find a convenient general setting for successful convergence all the time. Different combinations have helped me achieved convergence for the 4-6 different monolayers I've tested, but I had to 'babysit' each calculation.

Answer 1

I've managed to converge my optB88-vdW geometry relaxations to within forces of 0.002 eV/Å (EDIFFG=0.002), with an electronic convergence tolerance of $10^{-6}$ eV (EDIFF=1E-6). This was achieved partly by using a finer $\mathbf{k}$-point mesh and the conjugate gradient algorithm (IBRION=2), with ALGO=Fast. I think my convergence difficulties stemmed from using the qusi-Newton algorithm (IBRION=1) with the default step-size (POTIM=0.5) which should have been tuned.

For static calculations, I found that ALGO=All has been successful most if not all of the time.

My subsequent phonon and elastic constants calculations seemed to work out fine with my "tighter" geometry relaxations. However, I did find this paper titled "Energy, Phonon, and Dynamic Stability Criteria of Two-Dimensional Materials" (ACS Appl. Mater. Interfaces 2019, 11, 28, 24876–24884) which used optB88-vdW in VASP with ionic relaxation force tolerances of 0.01 eV/Å. They managed to get nice-looking phonon dispersion curves, which perhaps indicates that a higher force tolerance may be good enough.