k-points and ENCUT convergence tests before or after relaxation?

Question

Let's say I want to relax a structure using VASP. K-points and ENCUT convergence tests must be conducted before or after relaxation? and which tags should be included inside the INCAR file during convergence tests?

Answer 1

I do generally follow the follwoing:

1. Keep basic tags EDIFF, EDIFF to 1E-07 and 1E-06 (or -0.005) generally. i do use defaults if i want to do a faster run just for checking. and i dont know how can we preempt or know the potential surface is flat or not, as mentioned by Rosen. As in his answer said"...However, if the potential energy surface is flat, this may lead to issues converging the structure to the local minimum because of numerically inaccurate forces....."

2. for KPOINTS i use 30 to 40 times the reciprocal lattice parameter. (Eg. if a=3, b=4, c=6, then KPOINTS be 40/a, 40/b, 40/c). Not known the origin of such thumb rule !!! :). If its slab, surface or 2D material, then the Z-direction KPOINT will be 1.

3. i avoid using any other tags, except ISMEAR, ISYM and run static calculations for ENCUT and KPOINTS convergence. Still i face convergence issues.

but I did not face such issues so far in QE (just a basic learner now) Regards

Answer 2

This is a good question. In general, it often may not matter too much, but as with all aspects of numerical convergence, the only way for you to know for sure is to investigate. The concise answer is that in most cases, you can probably feel quite comfortable using the starting geometry for your convergence tests if: 1) the structure is reasonable (e.g. it's from experiments, taken from a database, converged at a different level of theory) and 2) the lattice constants do not change substantially after structure relaxation. Regardless, for any properties of interest, it is never a bad idea to check every now and then that the settings you are using are numerically converged.

This is a matter of opinion, but here is what I would suggest if you're particularly unsure. Carry out an (initial) structure relaxation with what are likely to be some reasonably accurate settings. The general rule-of-thumb for volume relaxation is to use (at least) 1.3 times the largest ENMAX value in your pseudopotential (POTCAR) files to prevent Pulay stresses. I often just use 520 eV for the plane-wave kinetic energy cutoff (ENCUT) since that's the maximum it could possibly be using the standard VASP 5.4 pseudopotentials. As for the number of $k$-points, that's going to be a bit trickier, but you can follow the lead of the OQMD or the Materials Project and use ~1000 $k$-points per number of atoms in the cell, distributed in a way that places more $k$-points along dimensions with smaller lattice constants. This would likely give you a quite reliable geometry to start with.

With this cleaner geometry, you can carry out your convergence tests and decide what you will use for the remainder of the project when studying this system. Of course, any new settings you decide upon must be used to re-relax your structure such that it is again at a local minimum in the potential energy surface for the production-quality settings you have chosen. Since the initial geometry and final geometry here are likely to be quite similar, you can feel comfortable not running yet another convergence test, although there is never harm in checking.

As for input flags, the following are the most important:

1. Plane-wave kinetic energy cutoff (ENCUT): Higher values are better but are more computationally expensive. I would start with the default obtained using prec='Accurate' and/or 1.3 times the maximum ENMAX value in your POTCAR files. Increase ENCUT by increments of ~50 eV or so until you feel comfortable with the results.

2. The number of $k$-points (KPOINTS): Again, higher numbers are better here. In general, you want to use more $k$-points along smaller lattice constants. There are several utilities in Pymatgen and elsewhere, such as the JHU $k$-point grid server, that can be very helpful with arranging the $k$-points. The ideal number of $k$-points changes significantly with cell volume, so this is one to keep an eye on if your structure changes drastically after relaxation.

3. The numerical tolerance for converging the electronic energy (EDIFF): A smaller number is better here. For ensuring accurate geometries, this value does not have to be too small. A default of 10E-4 is often okay. However, if the potential energy surface is flat, this may lead to issues converging the structure to the local minimum because of numerically inaccurate forces, in which case you would need to decrease the value further. For other properties, you may need a lower value of EDIFF as well, which I recommend decreasing in intervals of one order of magnitude until you are satisfied. Typically, I use 1E-6 for most of my work.

4. The numerical tolerance for the forces during structure relaxation (EDIFFG): This tells VASP when to stop the structure relaxation, and values closer to zero will bring you closer to the desired point in the potential energy surface. Generally, I would recommend going no higher than 0.05 eV/Å (EDIFF=-0.05) and often recommend 0.03 eV/Å as a good starting point. Try decreasing this value in intervals of ~0.01 eV/Å until you are satisfied. You can also use a tolerance based solely on energy differences between iterations, but I generally recommend against this in practice.

The following flags are less crucial but still important to consider:

5. The size of the integration grid (NGX, NGY, NGZ): Larger values are better here. Generally speaking, I have never found this to be an issue if you use prec='Accurate', which automatically sets the values for the integration grid to quite reasonable defaults. This will effect the energy as well as several derived properties, perhaps most notably are partial atomic charges, such as those via the Bader method. Some meta-GGA functionals have been shown to be quite sensitive to the integration grid, but I have found prec='Accurate' to still yield appropriate results in this case.

6. The smearing width (SIGMA). Depending on the smearing scheme used (ISMEAR), lower values often yield better numerical precision, but potentially at the expense of a more difficult convergence of the self-consistent field. You will want to confirm the the energy before and after extrapolation to the 0 K limit (from the fictitious temperature that is dependent on SIGMA) is reasonably close to one another. The closer these energies are, the less you have to worry about an inaccurate interpolation. If there are concerns, I often suggest using Gaussian smearing (ISMEAR=0) with SIGMA set to 0.01, which you can tweak by half-orders of magnitude or so depending on what you observe.

7. The number of bands (NBANDS). Higher numbers are better, and this will most notably influence band structures. Often, the default values are a reasonable start.

8. The number of grid-points when evaluating the density of states (DOS) (NEDOS). Higher numbers are better, and generally the default of 301 is not ideal. Increasing this to 2000 or so may yield better results. This will only influence the DOS, so you can neglect this setting if you are not interested in visualizing the DOS for your work.

An important note in all of this is that whether you are numerically converged or not will depend strongly on the property of interest. A geometry is going to be much less sensitive to numerical convergence than a electronic property like a band gap. Also, if you are interested in producing numerically precise energies, then it is worth remembering that all energies are only meaningful when they are relative. For instance, if you are modeling a reaction energy, there will be both a product and a reactant. In these cases, due to fortuitous error-cancellation, energy differences will converge significantly faster than absolute energies, and you may wish to perform your convergence tests on the difference to reduce the overall computational cost of your project.

You may also find my answer to "What are good ways to reduce computing time when working with large systems in VASP?" to be helpful.

Answer 3

Andrew's answer covers a lot of other considerations but not some useful additional ones in the case of relaxing unit cell parameters in solids. If the unit cell is being optimized in addition to atomic positions, the cell stress converges more slowly with respect to energy cutoff and k-points compared to forces and total energy. A good way to check this is to converge the calculated total stress to within 0.03 kbar with your initial structure (if it's reasonably close to the relaxed one, like an experimental structure). Then when you do your relaxation (making the total stress zero or a desired set value), an additional scf calculation is performed, resetting the G-vectors to those for the new unit cell (QE does this automatically, I'm not sure about VASP). If the stress from that calculation is not also zero, you are possibly not converged with respect to cutoff.