My Community

I have calculated the geometry of a Ag(111) slab + adsorbates and wanted to calculate the DOS. I used a 5x5x1 k-grid to geometrically optimize the slab. I checked the electronic properties using different MP k-grids. Surprisingly, the results changed a lot (EF = Fermi energy):

k = 5x5x1: EF = -0.4744, Total E = -208.6 eV
k = 9x9x1: EF = -0.3124, Total E = -202.5 eV
k = 13x13x1: EF = +0.2392, Total E = -190.4 eV

Changing the SIGMA parameter did not change these values. With the Fermi energy shift, also the position of all features in the DOS are shifted. I understand that the total energy could change when you change the number of k-points, but why does the Fermi energy change so much? And if it does change so much, how do I know if I can trust it or not? It is the density of states around EF in which I am interested in the end, so I should be able to calculate it more precise than with 0.7 eV precision, which is currently the EF-shift upon going from 5x5x1 to 13x13x1 k-grids.

Hi,

did you use a gamma-centred mesh?

alex

the results seem to be wrong: an increase of the k-mesh (even with different meshes) should NEVER lead to drops (!!) of the total energy by 6 (5 --> 9) and 18 eV (5-->13) if nothing else in the input is changed. I would definitely not trust these results
please carefully check your input and output files

I don't trust these results for sure, since it is counterintuitive that the larger k-grids lead to more unstable electronic configurations. However I cannot find something wrong.

My INCAR, POSCAR and POTCAR files are identical for the various calculations. Here is my INCAR file:

PREC = Normal
ENCUT = 400
IBRION = -1
ISYM = 0
ICHARG = 1
IALGO = 48
NELM = 60
NELMIN = 2
EDIFF = 1.0e-04
ISTART = 1
LWAVE = .FALSE.
LCHARG = .FALSE.
ISMEAR = 1
SIGMA = 0.2
LREAL = .FALSE.

NEDOS = 801
EMIN = -4
EMAX = 4

LORBIT = 11

The only difference in the KPOINTS is the number of k-points in the Monkhorst Pack grid. What should I be looking for in the OUTCAR or other output files to see where it goes wrong?

Post one KPOINTS file, please

Automatic Mesh
0
Monkhorst-Pack
9 9 1
0 0 0

Put "Gamma" instead of "Monkhorst-Pack". This will center your k-points at the gamma point and will retain the symmetry of the lattice, which is not guaranteed with your choice.
Convince yourself with a (small) example.

Cheers

Alex

Hello Alex,

Thanks for your reply. I thought that odd-numbered Monkhorst Pack grids were always Gamma-centered. I have visualized the k-space, given the individual k-points in the OUTCAR file, and this also shows my k-meshes are Gamma-centered.

Still I have made your suggested change and recalculated the 5x5x1 and 9x9x1 electronic configurations. The Fermi energy did not change from the Monkhorst Pack k-meshes however, so this weird thing is not due to his problem.

Does anyone know what could be wrong here?

I have tried it with different, easier systems, and found the same problem, although the differences are smaller for these smaller systems:

System: bulk Ag
One atom per unit cell
k = 5x5x5: EF = 2.6687, Total E = -2.94 eV
k = 9x9x9: EF = 3.3126, Total E = -2.84 eV
k = 13x13x13: EF = 3.1098, Total E = -2.82 eV
k = 17x17x17: EF = 3.0139, Total E = -2.82 eV
k = 21x21x21: EF = 2.9648, Total E = -2.82 eV

System: surface Ag(111)
21atoms per unit cell (7 layers, 3 atoms per layer)
k = 5x5x1: EF = -1.9366, Total E = -57.45 eV
k = 9x9x1: EF = -2.0050, Total E = -57.30 eV
k = 13x13x1: EF = -1.9700, Total E = -57.33 eV

I used exactly the same INCAR and KPOINTS file as I have posted before. Does anyone have similar problems? With which accuracy can I trust EF?

To obtain an accurate DOS it's in general required to have a dense k-point sampling. So the results for your two simple systems are not unreasonable (at least not to me at first glance). Just increase the k-point sampling further until the features in the DOS have converged.

The problem with your first input was that the total energy changed so dramatically with the k-point sampling. Assuming you have around 70 atoms you get differences corresponding to 0.1-0.15 eV/atom in the convergence in your first case compared to around 0.01 eV/atom in your surface Ag(111) calculations with 21 atoms. Therefore it seems that something has gone wrong with the input in the first case. Have you checked that all input were the same in the first case (except for the k-point sampling)? Are the positions exactly the same? Have you repeated the calculcations to see if it reoccurs? What kind of adsorbates to you add and do you add them on both surfaces?


/Dan

<span class='smallblacktext'>[ Edited Mon Oct 19 2009, 06:35PM ]</span>

Thanks for your reply. My adsorbates are thiophenol (a.k.a. benzenethiol, SC6H5), and they are adsorbed on both sides of the 7 layer slab, so I have 45 atoms in total. The input parameters are exactly the same (I compared them with the linux command diff) except for the one line in KPOINTS.

I have not recalculated the original calculation with too many k-points, this might take a long time. I did recalculate the 5x5x1 electronic relaxation anew, now using ICHARG=2 and ISTART=0 to have an unbiased electronic relaxation. The total energy is however again -208.6 eV, the Fermi energy is now at -0.4756 eV. No change here. Also the calculations on the Ag-only systems were performed twice each, on two different computer systems, but that gives identical results as well. That is why I don't expect any difference in outcome if I repeat the calculation with the same input with the finer k-grids in the original calculation.

If the input were the same, you might want to have a look at the output. Can you confirm that all three calculations with the different k-point sampling were performed using the same number of bands, same number of cores and that the symmetry identification was the same? Also, did all the calculations actually converge within NELM electronic steps (default 60)? This last point is very important so that you don't compare results that haven't reached self-concistency.

Finally, do you get any types of warnings in the OUTCAR file or the stdout ?

Best regards,
/Dan


<span class='smallblacktext'>[ Edited Tue Oct 20 2009, 04:07PM ]</span>

Thanks Dan. You help me a lot to show me where to look for errors.

Alas, I still cannot find any dissimilarities. I compared the three OUTCAR files carefully. None of them contained any errors, also the output to the screen contained no errors. Up to the part where the iterations start, the OUTCAR files are identical, including same number of bands, same number of cores and same symmetry identification. After the iterations (4, 41 and 42 iterations were necessary respectively), the total energy is different, as mentioned before. After that are some other changing results:

* average potential at core: 45 values in order of -50 eV, deviations of 0.1-0.2 eV occur
* pseudopotential strength for first ion, spin component 1: largest deviation ~0.2 eV for values around -15
* Total charge: deviates ~0.05 per ion per spd-character and ~0.5 in total for the three different OUTCARS, on a value of order ~250
* Total force on cell is very different:
FORCE on cell Total X Y Z
= 2.95 2.31 -2.76 for the 5x5x1 grid,
= 4.51 -0.36 2.16 for the 9x9x1 grid and
= 15.89 9.14 9.51 for the 13x13x1 grid
* Also the forces on the atoms are totally different, they are larger for the finer k-grid.
FORCES: max atom, RMS
= 0.04 0.016 for the 5x5x1 grid,
= 0.18 0.053 for the 9x9x1 grid and
= 0.64 0.155 for the 13x13x1 grid.

Could this mean that the forces are not properly calculated for the grids with fewer k-points, and that thus the system was not in a stable ionic ground state? I still don't understand why I cannot find the proper electronic ground state, even if I am not in the ionic ground state.

Have you tried increasing EDIFF (according to the manual you should have only 3 significant digits in your calcs--> total energy to the eV level) to for example 1.0e-8. you can simply do a continuation from your current results.

Danny

As long as the you do not update the ionic coordinates the forces will not have any impact on the total energy. I have assume you only performed 1 ionic iteration? The information you can get from the stress on the cell and the forces would be if you for example seem to be in a very unfavorable geometry, since in that case either the forces would be very high and/or the stress would be high. In your case the forces don't seem to be that high and to say anything about the magnitude of the stress is hard for a surface system. Maybe it seems a bit high in the x-direction. Is the normal to the surface directed in that direction? I guess not since according to the KPOINTS file you have put the single 1 point in the z-direction. However if the surface normal would be in the x-direction and the k-point sampling would be adjusted to the z-direction this would probably cause some strange results.

Honestly I don't have any more suggestions at the moment. Can you post your POSCAR file and let me know which pseudopotentials you are using? I.e. PAW or US, PBE or PW91, how many electrons for each specific potential and which order you have put them in the POTCAR file. Maybe then I can see if I can reproduce the results if the calculations are not too time consuming. Also which VASP version are you using?

Best regards,
/Dan


<span class='smallblacktext'>[ Edited Wed Oct 21 2009, 11:11PM ]</span>