extensive mae convergence tests on hcp Co

[Kevin Diekmann]

Dear Elk users and developers.

I am currently doing my master's thesis under supervision of Markus Meinert on mae. I have run quite extensive convergence tests in elk-3.3.17 for hcp Co in LSDA and GGA regarding all the relevant parameters defined by highq = .true. and more. Co was chosen for the convergence tests because of the simple hcp lattice and the large mae of about 90 micro eV with easy axis in the z direction.

The main insights are:

• The scf cycle for task 28 does not converge well. Better preconverge without spin orbit coupling with task 0 and then run task 29.
• LSDA yields the wrong easy axis for hcp Co but correct order of magnitude for super-ultra-high-quality parameters. This seems to be an intrinsic problem of the LSDA because the wrong easy axis is reproduced by vasp and in literature.
  It turned out, improvements by use of DFT+U are necessary for better outcomes.
• GGA tends to give the correct easy axis, but the scf cycle does not converge for lmaxvr > 15. This may either be a bug or numerical problem and is especially problematic because lmaxvr needs to be really large for convergence. Otherwise the preferred easy axis tends to switch and the mae is not continuous at all as a function of lmaxvr (in neither LSDA nor GGA).
• Of all the parameters in highq only lmaxapw, gmaxvr and nempty show nice convergent behaviour for a reasonable high value of these parameters.
• Very dense k-integration is needed. Checking for convergence with just ~5 values might trick you into thinking it converged but your values might just by accident lie next to each other. E.g. radkpt = 130 & 150 might yield the same value, but 140 is way off.
• The mae is nicely proportional to socscf^2 as expected. Thus in principle it should be possible to extrapolate to socscf = 1 with a linearized fit for problems with really small mae.
• There seems to be a minor bug in the mae subroutine, where it imports reducebf, but later sets it to .75 in any case.
• Just using the written value in MAE.OUT might fool you because it does not show the direction. Always check in MAE_INFO.OUT what the determined easy axis is as well.
• Regarding computational time, the first spin direction of task 29 is by far most computationally expensive.

All in all it does not seam to be possible to calculate mae with higher accuracy than about 40 micro eV within this framework and reasonable computational time because the value varies by change of nearly every inputparameter.

---

Details

I have run all tests with highq = .true. and swidth = 0.005 as the default values on hcp Co. Then varied all parameters in the highq list from highq to uhighq and partly beyond. This has been done for both GGA and LSDA. A value of the mae > 0 means that the preferred easy axis is into (the correct) z-direction. mae < 0 means into x-direction.

It turned out the scf cycle does not converge well in all cases by running task 28. First running task 0 and afterwards task 29 did converge in all cases except for the lmaxvr issue in GGA.

Input files were changed between tasks as follows.

---

tasks
0

spinpol
.true.

bfieldc
0.0 0.0 -2.0

reducebf
0.7

xctype
3

mixtype
3

highq
.true.

lmaxvr
8

lmaxapw
10

lorbcnd
.true.

nempty
10

swidth
0.005

epsengy
1e-7

avec
0.866025403 -0.5 0.0
0.0 1.0 0.0
0.0 0.0 1.0

scale1
4.71711

scale2
4.71711

scale3
7.63073

sppath
'elk_dir/species/'

atoms
1
'Co.in'
2
0.3333333333 0.6666666667 0.25 0.0 0.0 0.0
0.6666666667 0.3333333333 0.75 0.0 0.0 0.0

---

tasks
29

npmae
2

!all other parameters exactly as above

---

GGA

LDA

I have also tried run long tests on radkpt convergence for radkpt up to 150 for 5 values of from swidth = 0.001 ... 0.02, as well as autoswidth = .true. It does not converge at all for any of the temperatures, but the trend is obvious that the mae is going down for too big swidth. Using autoswidth = .true. yields a more systematic convergence but is still not converged until radkpt = 200. In any case, kpoint convergence and choice of swidth have to be done with a lot of care.

Anyway, kind of well converged paramters were chosen as follows named khighq.

        | uhighq  |  vhighq |   highq    | khighq

------------|---------|---------|------------|--------
rgkmax= | 10 | 9 | 8 | 12
gmaxvr= | 36 | 24 | 20 | 30
trimvg= | .true. | .true.| .true. |
lmaxapw= | 16 | 12 | 10 | 18
lmaxvr= | 12 | 9 | 8 | 18
lmaxinr= | 5 | 4 | 4 | 5
lmaxmat= | 12 | 12 | 8 | 13
fracinr= | 0.001 | 0.005 | 0.005 | 0.001
nrmtscf= | 4 | 2 | 1.5 | 6
nxlo= | 3 | 3 | 2 | 3
radkpt= | 120 | 90 | 50 | >= 150
autokpt= | .true. | .true.| .true. |
vkloff(:)= | 0 | 0 | 0 |
nempty= | 60 | 20 | 10 | 30
epspot= | 1d-7 | 1.d-8 | 1e-7 | 1e-7
epsengy= | 1e-5 | 1.d-6 | 1e-5 | 1e-8
epsforce= | 1d-4 | 1.d-5 | 1e-4 |
autolinengy=| .true. | .true.| .true. |
mixtype= | 3 | 3 | 3 |

I ran the tasks with khighq twice with radkpt = 150 & 145 to see if the result is anywhere near the literature. It yielded an mae of -68 micro eV and -27 micro eV respectively. These are then at least both in the same order of magnitude and in the same direction as other LSDA results of hcp Co [1]. Note that this is not the correct easy axis compared to experiment.
It also shows a kmesh as generated by radkpt = 150 does not give a converged result (as expected). This corresponds to a k-grid of 37x37x20 points. [1] argues a sampling of the order of 100x100x100 would be needed for a well converged result, but this was not possible within these tests as one of these khighq runs takes about 4 days on a newly bought 32 core workstation. This result (same order of magnitude and easy axis for similar number of k-points) is reproduced in vasp.

I hope this is of interest to some people. Feel free to contact me under kevin.diekmann@phyik.uni-bielefeld.de

[1] http://journals.aps.org/prb/abstract/10.1103/PhysRevB.41.11919

[JW]

Great report!

I have a question when comparing the MAE with GGA and LDA, did you use the same geometry?