Dealing with Convergence Problems - Quantum ESPRESSO

 Dealing with Convergence Problems

In the process of charge self-consistent calculation and structure relaxation, sometimes the pw.x calculation will reach the maximum number of iterations or exit with an error, which is a problem of non-convergence of the calculation.

On the one hand, the convergence problem is a physical reason, that is, the calculation does not conform to reality, including:

(1) Modeling is not in line with reality: the input atomic coordinates and element types are wrong; the calculation settings of symmetry.

(2) The approximation used is unrealistic: pseudopotentials with good performance are a necessary condition for the realization of calculations; k-space sampling method for metals; magnetic, strong spin-orbit coupling; transition metal elements with localized d and f orbitals On-site Coulomb repulsion (Hubbard U or Hybrid functional); Van der Waals force of molecular crystals, etc.; dipole electric field of surface slab model, etc. These non-convergence cases also deepen our understanding of the system, which may require corresponding methods to converge and correctly characterize the physical system.

On the other hand are the causes of numerical calculations, including matrix diagonalization errors, iterative algorithms, and structural optimization algorithms, resulting in errors or too many calculation iterations. The fact that the structure does not conform to the reality is the root cause of some numerical calculation errors: the structure that does not exist in practice will generate a charge distribution that deviates from the initial charge of the atom during calculation, which increases the difficulty of convergence. Of course, the reasons for non-convergence also include the instability of the numerical calculation algorithm itself, the contradiction between calculation accuracy and computer performance, and the lack of computer memory, hard disk space, number of hard disk files, and CPU resources.

The favorable condition for charge convergence is that the initial charge (spin) density distribution is close to the charge (spin) density distribution of the energy minimum state of the electron-ion system. The principle of achieving convergence is that the established model and the approximation used are in line with reality, and at the same time, the redundancy is appropriately reduced. If a charge does not converge in the calculation, the first thing to check is whether the structure is reasonable (for example, atoms in two overlapping positions, units of lattice constants, or even the wrong kind of element).

Sometimes something like the following occurs in the scf or nscf calculation:

     iteration #  1     ecut=    70.00 Ry     beta=0.80
     Davidson diagonalization with overlap
     c_bands:  1 eigenvalues not converged
     c_bands:  3 eigenvalues not converged
     c_bands:  1 eigenvalues not converged

This is a warning, which means that during the diagonalization iteration, the calculation of the highest 1-4 eigenvalues ​​does not converge (the threshold is controlled by diago_thr_init), and if more than 5 eigenvalues ​​do not converge (or nbnd/4), an error will be reported quit. It has a certain effect on the total energy and energy level. If it is not the last step, it can usually be ignored. The method described on this page may eliminate this warning.

Principles and Methods for Convergence of Self-Consistent Computation

Convergence of self-consistent calculation means that in the process of self-consistent calculation of charge, the total energy no longer changes significantly. The criterion is that the estimated scf accuracy in the pw.x output is less than the given value conv_thr, the default value (1.D-6) is relatively loose, and sometimes it is necessary to further improve the convergence criteria according to the calculated physical quantities.

Methods to achieve convergence:

method one:

Modify the charge mixing parameters. (1) Decrease mixing_beta. The default value is 0.7, which will be mixing_betaadjusted down to 0.3 ~ 0.1 ( even smaller 0.05 for very long slab). mixing_betaThe lower the value of , the less the new charges are mixed, and the closer the mixed charge is to the original charge, which is conducive to ensuring convergence.

(2) Increase mixing_ndim, that is, the charge of several steps back will be taken into account in mixing, the default is 8, and the memory required for calculation will increase.

mixing_ndim = 12

Method Two:

Increase the ecutwfcsum ecutrhoto increase the number of plane waves and the density of the charge and potential energy calculation grid (nr1, nr2, nr3). For supersoft pseudopotentials (US PP), the ecutrhodefault is ecutwfc4 times (4 times for modulo-conserving pseudopotentials), and it is increased to 10 times

ecutrho=10*ecutwfc

Method three:

For metal systems (or narrow band gaps, semi-metals, etc.): (1) Add some empty bands: increase nbnd, the specific value is based on the total number of electrons in the system, the default is at least 4 empty bands, increase to enough empty bands, so that the uppermost The occupancy probability of the empty band tends to zero; (2) increase the k-point grid density to eliminate the influence of the half-full band on the convergence of the total energy; (3) at the same time gradually increase the broadening until the total energy converges.

occupations = 'smearing', 
smearing='marzari-vanderbilt', 
degauss=0.01

For insulators and semiconductors,

occupations = 'fixed', 

or use a tiny degauss

occupations = 'smearing', 
smearing='gauss', 
degauss=1.0d-9,

Method four:

For systems with errors in diagonalization, after checking whether the structure is reasonable (for example, two atoms in the same position), try to use diagonalization='cg', the default is 'david', and set at the same time startingwfc='random'. Decrease the diagonalization threshold at the initial iteration diago_thr_init=1.0d-6(scf defaults to 1.0d-2, this option is also the nscf calculation threshold).

Method five:

Do you really need such a large vacuum size, in a slab model (slab), the thicker the vacuum and slab, the more problems you have with convergence. Have you checked the convergence of the surface energy (not just the total energy)? Start with a thin slab and gradually build up.

Principles and Methods for Realizing Structural Relaxation and Convergence

The criterion for structural relaxation convergence is that the total energy change of two adjacent ion steps is less than etot_conv_thr(1.0D-4), and the stress of each ion is less than forc_conv_thr(1.0D-3), and for vc-relax, the cell pressure is also required to be less than press_conv_thr( 0.5). By default, only nstepone ion step is performed (default 50). When using bfgs to optimize the structure, the ion step scf convergence threshold conv_thr will gradually decrease until conv_thr/upscale, upscalewhich is 100.0 by default. These convergence defaults are relatively loose, and it is necessary to further improve the convergence criteria if necessary. It is recommended to use the default ion_dynamics='bfgs', cell_dynamics='bfgs', and the 'damp'series method converges very slowly, and is only suitable for the initial close to the optimal structure.

If the charge in the first steps of the relax calculation converges normally, but an error is reported at a certain step, that is, the structure optimization does not converge, refer to the following method.

method one:

Do vc-relax on the primitive cell to get the theoretical lattice constant. The supercell is built with the theoretical lattice constant, and only ion is optimized for the supercell, that is calculation='relax', no vc-relax is performed.

Method Two:

If relax stops after several steps, this is sometimes because the band gap disappears after the system is optimized. To check the band gap, you can use the script here (gappw.sh relax.out). At this time, the self-consistent calculation is performed again according to the metal.

occupations = 'smearing', 
smearing='gauss', 
degauss=1.0d-2,

Method three:

If the relax stops after reaching the maximum number of iteration steps nstep, this is also a normal situation for a system with a large number of atoms (about 64). When the relaxation convergence standard is not too high, the structure of the last step can be used as the initial structure to continue the optimization. and increase nstep (default 50).

Method four:

Sharpening knives does not cut firewood by mistake, and calculates power with high precision. Appropriately increase the cut-off energy ecutwfc, ecutrhoand decrease it conv_thr(1.0d-8 to 1.0d-9), that is, to improve the convergence accuracy of scf to calculate a more accurate force. Decrease press_conv_thr to 0.1 during vc-relax calculation.

Method five:

When the calculation is initialized, the symmetry of the system will be automatically searched, and the amount of calculation will be reduced by the symmetry found. If the symmetry in the model is realistic, it should be retained, but there are also some cases, the initial symmetry is too high, and the symmetry of the system is spontaneous Defects (eg, ferroelectric, degenerate ground state Jahn-Teller effect), need to turn off the symmetry nosym=.true., so that the system can relax to a structure with low symmetry.

Method six:

Try not to perform lattice optimization on all atoms in a large system, for example, in a surface slab model, for internal slab or substrate atom fixation is not optimized.

Original post : https://yyyu200.github.io/DFTbook/blogs/2019/06/27/converge/