Dear Yujia Teng,
Thanks for taking the time to make these plots. Just for the benefit of other readers, the plots are for an InSe layer intercalated between layers of graphene and ZrTe2, not for a single graphene sheet or bulk ZrTe2.
I think there are a few critical points to consider from the get-go.
-
What do you consider as converged? You will have to define a tolerance threshold for the quantity of interest before you start your calculations. E.g., you want to converge your lattice parameter within +- 1% or to within +-0.05 Angstroms. And you might want to consider convergence reached only if N datapoints in a row fall within the predefined threshold.
-
You need to be aware of other factors that can change your results. E.g., the choice of functional (PBE vs. PBEsol, for example) might have a larger effect on the lattice parameter than running with ENCUT=ENMAX, or ENCUT=2 ENMAX. Of course, your results need to be converged with respect to all computational parameters, but there are other sources of errors as well, and just turning up a parameter to change the result by 0.5% is wasteful if another error is in the range of 5%.
-
You need to start your convergence studies with somewhat lower energies to better grasp the behaviour. For ENCUT, I would recommend starting at EMIN rather than ENMAX and then increasing the cutoff by 25 or 50 eV in each step.
-
The total energy of a system is an arbitrary number in DFT. Only energy differences (e.g., the energy difference between two phases, the energy of formation, or a surface energy) are meaningful.
-
Converging just one parameter (like the plane-wave cutoff energy) is usually not enough. Convergence behaviour might look different at a different KPOINTS mesh. It is cumbersome to converge all relevant parameters for a given problem, and they can be interdependent. But it is necessary to get trustworthy results.
Now to your concrete questions:
1. For the ZrTe2 heterostructure, the convergence is clear from both energy plots. Looks like it's 400 eV from total energy and 250 eV from energy difference. The dipole moment is not clear; it oscillates a lot but varies by about 1% of the initial value.
I would argue against your point for the total energy plots. Convergence is not clear, and the energy is also not normalized with respect to system size. You should plot the energy in e.g., meV/atom. As for the lattice parameters of the graphene system, the changes are very small, both in absolute value and as a percentage of the lattice constant. For the dipole moment of both systems: They are 0. You are mainly plotting noise.
2. The energy trend in graphene heterostructure is weird, first increase and then decrease, I think it's because the starting ENCUT is just ENMAX? And from energy, seems that graphene case doesn't converge until 800 eV, which is pretty high. Dipole moment change is also too small and can be viewed unchanged with ENCUT.
Here, I have to agree that it is a bit strange for the total energy plot. The total energy is a variational quantity with respect to the plane-wave energy cutoff. The energy should decrease with increasing number of plane waves, and the converged solution should be approached from above. Maybe the situation would change if the k-point density is increased? I already commented on the dipole moments above.
So I did the test with E-a method, which takes lots of time due to huge amount of jobs.
Note that it is better to converge with respect to the quantity you want to calculate (e.g., interlayer distance in case of a slab) than to do the Equation of state fits. This is just a common approach to illustrate the advantage of converging with respect to energy differences.
I hope this helps to clear up your confusion. Please let me know if there are still things that are unclear.
Cheers, Michael
