Convergence test for Ecut + k mesh points

[karianaandrea_morenosader]

Hi all,

Hope you’re doing well.

I’ve run a convergence test (Ecut and K-points) for a primitive Pt cell, using absolute energy as the property of interest. I’m using PAW_PBE pseudopotentials and the revised PB functional.

For the Ecut convergence, I used a k-mesh of 25×25×25 and observed convergence at Ecut > 650 eV. For the k-mesh convergence, I used Ecut = 650 eV and found convergence at 33×33×33. However, I’ve seen some examples where convergence is reached at slightly lower settings.

1. I’m computing the difference in absolute energy between consecutive Ecut values and checking when the difference is <1 meV/atom.
Q1: Is this a reasonable approach, or should I instead compare each Ecut value to the lowest energy obtained so far?

2. I notice some oscillations in the energy difference between consecutive Ecut values after Ecut > 650 eV of around 0.3 meV/atom.
Q2: Is this numerically expected? (Maybe the figure below is more illustrative)

3. If I want to build a slab from this bulk material, I was told to relax the primitive cell first. I first manually varied the lattice parameter and found the equilibrium value at 3.99 A. Then I did relaxation using ISIF=3 and ISIF=4
Q3: I wonder if symmetry should be turned off for any of these? I have done all my tests before with symmetry on (default ISYM).

4. ISIF=3 and =4 gave the same optimized lattice and energy, with external pressure = 0.11 kB Pullay stress = 0.00 kB.
Q4: What is a reasonable threshold value for the external pressure to indicate that the primitive cell is well converged?

Thank you very much for your response!

Below are my INCAR file settings:
SYSTEM = bulk fcc Pt primitive cell (ecut400)
ENCUT = 650 # Plane-wave cutoff energy in eV
EDIFF = 1E-7 # SCF energy convergence criterion
GGA = RP # Revised Perdew-Burke-Ernzerhof functional
LWAVE = .FALSE. # No wavefunction file output
LCHARG = .FALSE. # No charge density file output
ISMEAR = 1 # Methfessel-Paxton smearing
SIGMA = 0.1 # Smearing width in eV
PREC = Accurate
LREAL = .FALSE. # Real-space projection
KPAR = 2 # Parallelization over k-points
NCORE=1
ALGO = Normal # Electronic minimization algorithm

NELM = 200 # Maximum electronic steps
EDIFFG= -0.01
IBRION=2
NSW = 120

[ahampel]

Hi,

thank you for reaching out on the VASP forum. We are currently working on a guide for studying encut and kpoints convergence for the VASP wiki. So hopefully what I will write in the following will also be found there in the near future. Currently we are lacking such tutorial.

For the Ecut convergence, I used a k-mesh of 25×25×25 and observed convergence at Ecut > 650 eV. For the k-mesh convergence, I used Ecut = 650 eV and found convergence at 33×33×33. However, I’ve seen some examples where convergence is reached at slightly lower settings.

this is good. Always first find encut and then convergence kpoints. Make also sure to study convergence per f.u. or atom/ electron in the system. For you system it is of course only one f.u.

1. I’m computing the difference in absolute energy between consecutive Ecut values and checking when the difference is <1 meV/atom.
Q1: Is this a reasonable approach, or should I instead compare each Ecut value to the lowest energy obtained so far?

I typically use the largest encut calculated as reference value. Then I can define the "floor". Because, encut will behave asymptotically and more plane waves is always more accurate this should be the safest definition for determining the encut. I typically define converged when Etot (without entropy) change is smaller than 10 meV per atom. However, be careful with metals here as the smearing will be part of the convergence. Especially for forces and metals a higher number of kpoints and tetrahedron method (ISMEAR=-5) is advised to be used.

2. I notice some oscillations in the energy difference between consecutive Ecut values after Ecut > 650 eV of around 0.3 meV/atom.
Q2: Is this numerically expected? (Maybe the figure below is more illustrative)

This should not happen. Either the kpoint mesh is not accurate enough here, the smearing is not good enough (ISMEAR=-5) should be safest for metals), or your number of PW is really crazy high that some numerics gets messed up. But I doubt that it is the last of these options. I looked at your INCAR and please try ISMEAR=-5 . Also please try to use the "Pt_pv" PAW_PBE POTCAR. This one will be more accurate especially for cell relaxations. The normal Pt is fine but the Pt_pv has more PAW projectors and is closer to an all electron result. But of course convergence for encut and kpoints can be reached with both.

3. If I want to build a slab from this bulk material, I was told to relax the primitive cell first. I first manually varied the lattice parameter and found the equilibrium value at 3.99 A. Then I did relaxation using ISIF=3 and ISIF=4
Q3: I wonder if symmetry should be turned off for any of these? I have done all my tests before with symmetry on (default ISYM).

This is almost more a science question. Typically it is best practice to first relax with symmetry. If you have no idea whether the used symmetry is correct or maybe the material wants to go into a different structure. You can calculate a phonon spectrum on the relaxed spectrum to see if there are imaginary phonon modes. Then it is typically more advisable to construct a unit cell with this symmetry. If you remove all symmetry all together it could happen that you will just get garbage in the relaxation. Again for your simple unit cell it might just work. But I recommend the phonon approach to find instable phonon modes and follow those.

4. ISIF=3 and =4 gave the same optimized lattice and energy, with external pressure = 0.11 kB Pullay stress = 0.00 kB.
Q4: What is a reasonable threshold value for the external pressure to indicate that the primitive cell is well converged?

I tried your files and was able to basically remove all pressure that is printed. Try to change these parameters:

Code: Select all

ISMEAR=-5
EDIFF= 1E-10
EDIFFG=-0.0001

This took only 3 ionic steps for me and changed the volume slightly. That said an external pressure below 1 kB is not bad, but for this cell one can get two more digits to 0 without too much computational effort. This computes forces and stress down to a precision of 10^-4 eV/A . I would consider this the gold standard when reviewing a paper.

I hope these steps help a bit. Give this a try and I am happy to discuss further.

Best,
Alex

[karianaandrea_morenosader]

Hi Alex,

Thanks a lot for the detailed response.
Q1: Sounds great! I changed the way I was calculating the energy difference to use the largest calculated ecut (850 eV) as my reference.
Q2: I re-ran the ecut convergence tests using ISMEAR = -5 and Pt_pv and was able to get rid of the oscillations (see attached figure). I am now getting energy differences of less than 10 meV after 550 eV. Thanks!
Q3: I also think that for my simple cell it might work to remove symmetry, but I will do a vibrational frequency analysis.
Q4: I ran with your recommended settings, and ISIF = 3 removed all external pressure as well. However, I read in the VASP manual that ISMEAR = -5 is not recommended for metal relaxations, and I wonder if the slight distortion in the volume with ISMEAR = -5 should be fine—like having 0.5002892380305481 instead of exactly 0.5 for my fcc cell

3.990000000000000
0.5002892380305481 0.5002892380305481 0.0000000000000000
-0.0000000000000000 0.5002892380305481 0.5002892380305481
0.5002892380305481 -0.0000000000000000 0.5002892380305481

I also did a manual swap of lattice parameters to values closer to those from the relaxation. I computed total energies with ISMEAR = -5 (without relaxation) for these lattice parameters and found one that gave me almost zero external pressure and Pulay stress, while the cell shape remained with 0.5.

3.992295
0.5 0.5 0.0
0.0 0.5 0.5
0.5 0.0 0.5

Total 0.00004 0.00004 0.00004 0.00000 -0.00000 0.00000
in kB 0.00365 0.00365 0.00365 0.00000 -0.00000 0.00000
external pressure = 0.00 kB Pulay stress = 0.00 kB

Do you think it should be fine to have 0.5002892380305481 from the relaxation with ISMEAR = -5 instead of exactly 0.5?

Thanks a lot for all your feedback, I really appreciate it!

[ahampel]

Hi,

Do you think it should be fine to have 0.5002892380305481 from the relaxation with ISMEAR = -5 instead of exactly 0.5?

I think it makes sense to move this shift in all lattice vectors to the scaling factor on top. Ideally, like here they are all the same. Actually, I did not consider the fact that ISMEAR=-5 is not well suited for highly accurate forces. So yes, please try to avoid using this for metals and use instead ISMEAR=0 or larger with small SIGMA values (depending of the k-mesh).

I think your final encut convergence plots looks all good to me! We hope to have a tutorial on the wiki to reflect these tips we discussed here up soon.

Best regards,
Alex

[karianaandrea_morenosader]

Hi Alex,

Thanks a lot! One last thing regarding the phonon calculation: I followed your recommendations and got the relaxed cell with zero external pressure. Then, to obtain the phonon modes, I performed DFPT using IBRION=7. I did not observe any negative phonon modes, and the values I am getting are very small. However, I am not sure whether these can be considered essentially zero (within numerical noise)?

Thank you very much!

INTERNAL STRAIN TENSOR FOR ION 1 for displacements in x,y,z (eV/Angst):
X Y Z XY YZ ZX
--------------------------------------------------------------------------------
x 0.00550 0.00129 0.00129 0.00000 0.00000 -0.00000
y 0.00127 0.00554 0.00127 -0.00000 0.00000 0.00000
z 0.00022 0.00022 0.00596 0.00000 0.00000 -0.00000

--------------------------------------------------------------------------------------------------------

SECOND DERIVATIVES (NOT SYMMETRIZED)
------------------------------------
1X 1Y 1Z
1X -0.002050 0.000000 -0.000000
1Y 0.000000 -0.002265 -0.000000
1Z 0.000000 0.000000 -0.001983


Eigenvectors and eigenvalues of the dynamical matrix
----------------------------------------------------


1 f = 0.053275 THz 0.334736 2PiTHz 1.777058 cm-1 0.220327 meV
X Y Z dx dy dz
0.000000 0.000000 0.000000 0.000000 1.000000 -0.000000

2 f = 0.050678 THz 0.318421 2PiTHz 1.690447 cm-1 0.209589 meV
X Y Z dx dy dz
0.000000 0.000000 0.000000 1.000000 -0.000000 0.000000

3 f = 0.049837 THz 0.313136 2PiTHz 1.662389 cm-1 0.206110 meV
X Y Z dx dy dz
0.000000 0.000000 0.000000 0.000000 0.000000 1.000000

Eigenvectors after division by SQRT(mass)

Eigenvectors and eigenvalues of the dynamical matrix
----------------------------------------------------


1 f = 0.053275 THz 0.334736 2PiTHz 1.777058 cm-1 0.220327 meV
X Y Z dx dy dz
0.000000 0.000000 0.000000 0.000000 0.071597 -0.000000

2 f = 0.050678 THz 0.318421 2PiTHz 1.690447 cm-1 0.209589 meV
X Y Z dx dy dz
0.000000 0.000000 0.000000 0.071597 -0.000000 0.000000

3 f = 0.049837 THz 0.313136 2PiTHz 1.662389 cm-1 0.206110 meV
X Y Z dx dy dz
0.000000 0.000000 0.000000 0.000000 0.000000 0.071597

[ahampel]

Hi,

I think those numbers are well converged. I assume you only calculated phonons in the primitive unit cell? I should have been more precise. Only phonon modes beyond Gamma will show us breaking of space group symmetry: https://www.vasp.at/wiki/index.php/Comp ... on_and_DOS . So you would need to do the IBRION=7 calculations within a super cell that is large enough to host the phonon mode that break the space group symmetry. Or use this tutorial as guideline: https://vasp.at/tutorials/latest/phonon/part2/ .

Alternatively, if this is a bit too cumbersome it would be of course an alternative to turn off symmetries and relax the structure.

Best,
Alex