internal error in asa.F and in YLM3LOOKUP: L or LP > LMAXCG

[alpinnovianus]

Dear admin,

I encounter this bug when I am running calculations for Nd7Ce1Cu4O16 system.

--------------
!Bug!
internal error in: asa.F at line: 479
internal error in YLM3LOOKUP: L or LP > LMAXCG
L, LP, LMAXCG 7 1 6
If you are not a developer, you should not encounter this problem.
Please submit a bug report.
--------------
I attach the error files as well as the input files below.

error_inputfiles.zip

I remark that there was no such error when I ran calculations for Nd8Cu4O16.
Could you please let me know if I made a mistake in my input files for Nd7Ce1Cu4O16?

I use vasp 6.2.1.

Thank you.
Alpin

[alpinnovianus]

I did a simpler check with only Ce atoms:
POSCAR:

Ce4
1.0
0.0000000000000000 5.7713080000000003 0.0000000000000000
3.3242200000000000 0.0000000000000000 0.0000000000000000
0.0000000000000000 -1.9276719999999987 -5.4396610000000010
Ce
4
direct
0.2446390000000001 0.5000000000000000 0.7504950000000000 Ce
0.7553609999999999 0.5000000000000000 0.2495050000000000 Ce
0.7446390000000001 0.0000000000000000 0.7504950000000000 Ce
0.2553609999999999 0.0000000000000000 0.2495050000000000 Ce

and I narrowed it down that this error happens only when METAGGA is switched on.
I have found this to happen when METAGGA = SCAN and R2SCAN. I didn't check the others.

INCAR (uncommenting the metagga tag below will result in this error):

MAGMOM = 4*5.0

LREAL = .FALSE.
NPAR = 4
KPAR = 4
ISIF = 3
LWAVE = .FALSE.
LCHARG = .FALSE.

GGA = PE
#METAGGA = SCAN
LASPH = .TRUE.
LMIXTAU = .TRUE.
ADDGRID = .TRUE.
ALGO = Normal
NELM = 200
PREC = A

IBRION = 1
NSW = 200

EDIFF = 1E-05
EDIFFG = -0.008
ENCUT = 550
ENAUG = 650

ISMEAR = 0
SIGMA = 0.01

LORBIT = 11
ISPIN = 2
LMAXMIX = 6
LMAXTAU = 8

KPOINTS:

pymatgen with grid density = 608 / number of atoms
0
Gamma
4 7 4

POTCAR: the standard Ce POTCAR.

PAW_PBE Ce 23Dec2003
12.0000000000000
parameters from PSCTR are:
VRHFIN =Ce : [core= Kr 4d10]
LEXCH = PE
EATOM = 1063.0861 eV, 78.1346 Ry

TITEL = PAW_PBE Ce 23Dec2003
LULTRA = F use ultrasoft PP ?
IUNSCR = 1 unscreen: 0-lin 1-nonlin 2-no
RPACOR = 1.800 partial core radius
POMASS = 140.115; ZVAL = 12.000 mass and valenz
RCORE = 2.700 outmost cutoff radius
RWIGS = 2.500; RWIGS = 1.323 wigner-seitz radius (au A)
ENMAX = 273.042; ENMIN = 204.781 eV
RCLOC = 1.810 cutoff for local pot
LCOR = T correct aug charges
LPAW = T paw PP
EAUG = 580.196
DEXC = 0.000
RMAX = 2.749 core radius for proj-oper
RAUG = 1.300 factor for augmentation sphere
RDEP = 2.810 radius for radial grids
RDEPT = 2.162 core radius for aug-charge

Atomic configuration
14 entries
n l j E occ.
1 0 0.50 -40258.9816 2.0000
2 0 0.50 -6445.4884 2.0000
2 1 1.50 -5771.7470 6.0000
3 0 0.50 -1387.9739 2.0000
3 1 1.50 -1173.8158 6.0000
3 2 2.50 -869.9108 10.0000
4 0 0.50 -280.8706 2.0000
4 1 1.50 -210.9815 6.0000
4 2 2.50 -110.4558 10.0000
5 0 0.50 -40.5819 2.0000
6 0 0.50 -3.7444 2.0000
5 1 1.50 -23.2188 6.0000
5 2 2.50 -3.0598 1.0000
4 3 2.50 -5.3970 1.0000
Description
l E TYP RCUT TYP RCUT
0 -40.5819164 23 1.500
0 -3.7443793 23 2.300
1 -23.2188445 23 2.000
1 40.8174780 23 2.300
2 -3.0597556 23 2.500
2 13.6058260 23 2.500
3 -5.3970469 23 2.700
3 -5.8505052 23 2.700
....


Moreover, for this simple check on Ce atoms only, I can use METAGGA = SCAN without seeing this error if I change the POTCAR from the standard Ce (tetravalent) to the Ce_3 (trivalent, f in core) POTCAR.

The problem is since I want to study Ce-substituted Nd-cuprate as mentioned in the first post:
I believe I need to use potential for tetravalent Ce in my metaGGA calculations in order to obtain electron doped conditions for this cuprate material: Nd valency: +3 --> Ce valency= +4

Do you have any suggestions for getting around this bug/error, please?

[alpinnovianus]

also, I have checked the Ce_h PBE POTCAR too in the Ce-only structure test. They also gave the same error when computed with METAGGA=SCAN.

[ferenc_karsai]

You have specified "LMAXTAU=8" in your calculation.
This value is not supported, because the maximum l value in asa.F is 6.
I also talked to Georg Kresse and he said values above 6 make no difference to the results, that's why the upper limit is hardcoded.

By setting "LMAXTAU=6" your problem should be solved.

[alpinnovianus]

You have specified "LMAXTAU=8" in your calculation.
This value is not supported, because the maximum l value in asa.F is 6.
I also talked to Georg Kresse and he said values above 6 make no difference to the results, that's why the upper limit is hardcoded.

By setting "LMAXTAU=6" your problem should be solved.

thank you for your reply.

Can you elaborate why, please?

The choice of LMAXTAU=8 is because of this wiki content. I reckon since Ce has f-elements it should be 2*3+2 = 8.
the last wiki sentence is warning that f-elements may need value larger than 6.

https://www.vasp.at/wiki/index.php/LMAXTAU

The PAW one-center expansion of the density has component up to and including L=2*lmax, where lmax is the l-quantum number of the partial waves on the POTCAR file, with the highest angular moment. If the PAW one-center expansion of the density has component up to L, then the one-center expansion of the kinetic energy density has components up to L+2.

This means that as a rule of thumb, for s-elements: LMAXTAU=2, for p: LMAXTAU=4, and for d: LMAXTAU=6. If you are willing to live with the computational costs, the default for LMAXTAU should be safe in all cases, except those involving f-elements.

[ferenc_karsai]

Ok I've talked to Georg Kresse again.

He said that he has tested that if you don't use kinetic energy density like in meta GGAs LMAX=6 enough. That is why it is hard-coded in asa2.F.
On second thought, for meta GGAs he is not shure. So we will test this carefully and change the default in an update of VASP.

Mind it's not enough to change asa2.F but also the calling routines.

So for the moment you have to use LMAXTAU=6 (possibly 7 works too, please try it).

[alpinnovianus]

On second thought, for meta GGAs he is not shure. So we will test this carefully and change the default in an update of VASP.

May I know if you guys have made any changes related to this in vasp 6.4.0? (or figured out whether it is safe to use LMAXTAU=6 only?)

[matthewkuner]

I am also interested in this--do you have any updates @ferenc ?

[ferenc_karsai]

I've brought this up in a development meeting. There is also an issue about this in our git repository. Unfortunately I fear that nobody had time until now to work on that problem. I will bring that up again, especially since already more than one of you is interested in that feature.

[ferenc_karsai]

Ok, I have changed asa.F to support lmax=8 and we have looked at the convergence for Ce that you have posted. The energy difference between a calculation using LMAXTAU=6 and LMAXTAU=8 is in the sub-sub meV region per atom for this particular system (LMAXTAU=8 has the lower energy). So for this system we conclude that going beyond LMAXTAU=6 makes no significant differences. We don't know if this is a general rule or particular only for few systems, since it would require much more test systems.

We looked into the code and haven't found anything that would prohibit a proper functioning of LMAXTAU=8. But to be entirel shure one would need very extensive investigations. One thing we saw is that the calculation using LMAXTAU=8 required more steps to converge than LMAXTAU=6, but this could mean anything.
So until now we can't say that it is safe to employ the changes in the main branch of VASP and it is unlikely that it will be included in the next patch.

Nevertheless if you want to try the fix, you need two changes:
asa.F:
change
INTEGER, PARAMETER :: NCG=13000 ! LMAX= 6 (f electrons)
to
INTEGER, PARAMETER :: NCG=45000 ! LMAX= 8 (f electrons)

and

main.F:
change
LMAX_TABLE=6; CALL YLM3ST_(LMAX_TABLE)
to
LMAX_TABLE=8; CALL YLM3ST_(LMAX_TABLE)


Good luck!


PS: In the INCAR for your Ce example you set "MAGMOM=4*5.0". Since the system has very small magnetization (around 0.07) we achieved much faster convergence starting from "MAGMOM=4*0.0".