Strange results for the DOS projection in the Bader volumes

[forsdan]

Dear Wenjie Tang and Graeme Henkelman,

I tried to use the new available code for the DOS projection in the Bader volumes, but I'm puzzeled about the results. I am looking at a bulk Fe system with 4 atoms where all atoms are in the same enviroment. However, when I examine the lm-projection on each atom, the sum of the projections on the first atom is of the order of the total DOS while the projections on the remaining three atoms are identical to the projections I obtain if I comment out ICHAIN = 5.

So there is two problems when I use ICHAIN = 5: 1) The symmetry between the atoms are significantly broken, 2) The total sum of all the atomic projections are larger than the total DOS. From the Bader analysis the charge and the minimum distance are the same for each atom, and in the BDRCAR the grid points are equally distributed between the atoms. So the symmetry should be preserved.

So I am wondering if I use the program correctly. The procedure I'm using is the following:

1. Run a scf-calculation with the INCAR, POSCAR, KPOINTS, POTCAR files given below.
2. Generate AtIndex.dat with "./bader -p atom_index CHGCAR" and rename it to BDRCAR
3. Set ISTART = 1, ICHARG = 0 and Add
ICHAIN = 5
LORBIT = 12
NEDOS = 1000
EMIN = -9.5
EMAX = 15.5

to the INCAR file and restart VASP.

I also tried running with LORBIT = 10, ISYM = 0 and ISPIN = 1, and combinations as well (the above procedure was repeated in each case). The results are always the same when I use ICHAIN = 5. The first atom get too much states, while the other three get the default projections. I also ran without NEDOS, EMIN and EMAX but I got the same features.

I'm quite confident that my summations of the projections in the DOSCAR file are correct. So I wondering if there is something important that I'm missing?
I am using Vasp5.2.2.

Best regards,
/Dan


INCAR
-------

SYSTEM = Fe (Bulk)
ISTART = 0
ICHARG = 2
EDIFF = 1E-5
ENCUT = 500 eV
ENAUG = 700 eV
PREC = Accurate
LREAL = .FALSE.
ISMEAR = 1
SIGMA = 0.2
NSW = 0
NELM = 400
VOSKOWN = 1
ISPIN = 2
MAGMOM = 4*3
LMAXMIX = 4

POSCAR
---------

Fe_bulk
1.0000000000000000
2.9129299999999998 0.0000000000000000 0.0000000000000000
0.0000000000000000 2.9129299999999998 0.0000000000000000
0.0000000000000000 0.0000000000000000 5.4930599999999998
4
Direct
0.0000000000000000 0.0000000000000000 0.0000000000000000
0.5000000000000000 0.5000000000000000 0.2500000000000000
0.0000000000000000 0.0000000000000000 0.5000000000000000
0.5000000000000000 0.5000000000000000 0.7500000000000000


KPOINTS
----------
Monkhorst Pack
0
Monkhorst Pack
16 16 8
0 0 0


POTCAR
---------

PAW Fe
8.00000000000000000
parameters from PSCTR are:
VRHFIN =Fe: d7 s1
LEXCH = 91
EATOM = 595.5622 eV, 43.7726 Ry

TITEL = PAW Fe
LULTRA = F use ultrasoft PP ?
IUNSCR = 1 unscreen: 0-lin 1-nonlin 2-no
....

[Wenjie]

If possible, can you send me your BDRCAR and CHGCAR file, so I can check what happened with the calculation.

[Wenjie]

Dear Dan,

It seems that I have some problem to open your email. Is possible to send me a link so I can download the files? Thank you.

Wenjie

[forsdan]

Dear Wenjie,

Thank you for your replies. I was concerned that there would be a size or a tar-achive issue. I have now sent you a download link to the files. Hope it works properly otherwise please let me know.

Best regards,
/Dan

[forsdan]

Dear Wenjie,

I want to give an update to the above problem:

I abandoned my current system for now and instead tried the example given on the homepage. In that way I had results to compare with. I also compiled VASP 4.6.34 with the bdrpro.F, in addition to Vasp5.2.2.

For Vasp4.6.34 I get the following:
----------------------------------------

1. When I use NPAR = 1 as given in the INCAR file, I obtain the same results as in the supplied DOSCAR file.
2. When I comment out NPAR = 1, and instead use the default value I get the wrong results. The results resembles the situation for my own case where the first atom get too much states.
3. If I comment out ICHAIN = 5, I obtain the identical results for NPAR = 1 and for NPAR set to the default value, with the standard spherical projections.

An additional side point is also when I try to use LORBIT = 12 instead of LORBIT = 10 (together with NPAR = 1 and ICHAIN = 5) , I get that the columns 2 to 7 in the DOSCAR file for each atom projection with LORBIT = 12 are identical to the columns 2 to 7 obtained with LORBIT = 10. Then the remaining columns 8 to 19 for LORBIT = 12 seem to be something else, since if I include them in the summation I get much more states than the total DOS. This structure of the columns differ from the standard spherical projection scheme, which contains the lm-projections. Does this mean that I can't get the m-projections in the current implementation?

Also, if I integrate the columns 2 to 7 in DOSCAR (in the supplied tar-achive of completed files) up the Fermi-level and sum up the contributions for all the atoms, I seem to lack approximately two electrons. Now I just used a simple trapezoid integration so there are probably some numerical errors, but just to be sure: how close I can expect to get to the total DOS by using the current projection scheme? Should I always get close to 100% of the electrons?


For Vasp5.2.2 I get the following:
-------------------------------------

1. When I use NPAR = 1 as given in the INCAR file, the calculation either crashes or hangs (depending on the cluster I compiled it on) when it enters the Bader projection part.
2. When I use NPAR set to the default value I get the wrong results. The results are the same as for Vasp4.6.34 with the same settings.
3. If I comment out ICHAIN = 5, I obtain that the results with NPAR = 1 are identical to the results with NPAR set to the default value, with the standard spherical projections. The two settings then work as expected.


I get the same results on three different clusters. The binaries were all compiled with using MKL and ifort of different versions, together with either OpenMPI or MPICH. So my guesses are

1. The NPAR = 1 failure with Vasp5.2.2 is a compiler issue of some sort, although when I turned off the optimization the problem still persisted. I will look into it in more detail.

2. The NPAR > 1 issue is related to the memory distribution across the nodes, so when the calculation enters the algorithm it produces the incorrect results. Can you please just confirm that I have to use NPAR = 1 for the projection scheme to work or is it supposed to work for NPAR > 1 as well ?


I appreciate any comments that you have time to spare. It would be really helpful so that I know that I am doing the projections correctly.

Best regards,
/Dan

[andyfan]

Dear, Dan,

Does your problem figure out? I got similar problems to yours that the same atoms obtain different site-projected DOS.

[forsdan]

Dear andyfan,

To clarify my above post I have (still) three questions:

1. Does the algorithm require NPAR = 1 for a parallel run due to the memory distribution?
2. Are the m-projections available?
3. Will the current projection scheme account for all electrons if the atom-projections are summed up? Increasing NG(X,Y,Z)F and/or the k-point sampling doesn't seem to help in the test-case.

If I use NPAR = 1 and LORBIT = 10, I get reasonable results for my systems and an identical DOSCAR in the test case. So given that these settings should be used, the projections seems to work fine and the only question mark is point 3.

Best regards,
/Dan

[andyfan]

Dear Dan,

Using the input file of the test example, I can also obtain the identical DOSCAR to the test one, when using NPAR = 1 and LORBIT = 10. But when I performed similar calculations on a slab model, the results seemed incorrect. The partial DOS of all the surface atoms should be identical (This is indeed the truth if I do a regular DOS calculation). However, the results from the DOS projection in Bader volumes are very different for the surface atoms.

I noted that the test example is a Cu3Pt bulk. So I also perform a similar calculation on a Fe2Ni2 bulk. In that case, I can get reasonable result. So I wonder this code could not accurately deal with the slab model, where a vaccum layer is involved.

[Wenjie]

Dear Dan,

I've figured out the problem for your calculation. I've updated the dos webpage: http://theory.cm.utexas.edu/vtsttools/dos/
If you modify your main.F according to the webpage and recomplie, the code should work. Please make sure you set NPAR=1, which is very important. This time you should get identical DOS for each atom and the sum of atomic DOS should add up to the total.

Also, just a reminder, LORBIT=12 does the lm-decomposed DOS, but our split_dos and dosanalyze.pl codes don't work for this case. Please make sure you have the correct scripts to deal with lm-decomposed DOS.

Wenjie

[Wenjie]

andyfan,

Because the volumes in the vacuum are assigned to the surface atoms, the high energy states will have large contribution, which makes the DOS weird. When calculating surface atoms, we recommend the Bader code with version 0.27 to get rid of the vacuum volumes. This code should be released soon.

Wenjie

[andyfan]

Thanks, I am looking forward to the new Bader code.

[forsdan]

Dear Wenjie,

Thank you very much for your reply. The updated code indeed fixes the LORBIT = 12 issue (I have been using my own scripts the entire time) and I get identical results when I sum up the m-projections for each l-component as I get for LORBIT = 10. So that part seems to be working fine now.

However, I have to bother you again about the claim that the atomic DOS should add up to the total DOS. I have been over it several times now, but I do not get that the sum of the atomic DOS adds up to the total if I sum up the supplied DOSCAR file for the test case (I also re-ran the input-files and confirmed that I get the same DOSCAR). I also added NEDOS = 1000, PREC = Accurate and doubled the grids in each direction, i.e NG(X,Y,Z) = 40 and NG(X,Y,Z)F = 80, but the results are the same. Note that BDRCAR was regenerated in this case. I have plotted the total DOS, the sum of projections and the difference in both cases:

http://fy.chalmers.se/~forsdan/BADERDOS/files.html

where the resemblence between the difference and the total dos to some extent seems to suggest that there is states partially missing across the entire interval. The relative difference also tends to increase with the energy. So I am puzzeled about that the atomic DOS should add up to the total, when I can't even get a proper result for the supplied DOSCAR. I also got even larger differences for another of my systems (an interface). Is it with 100% certainty I should obtain exact agreement between the sum of the projections and the total dos with the current implementation? The summations are very straightforward so I am quite confident in them. Do you have any suggestion on how to sort this out?

Best regards,
/Dan

[Wenjie]

Theoretically the sum should be exactly the same as the total. However, there's always a very difference. I may still need some tests to figure out where the difference comes from. I'll let you know when I have the tests done.

[forsdan]

Dear Wenjie,

I was just wondering if you possibly have had time to look into the above issue yet. The files at the homepage are still identical to the ones I did use for the above investigations, so the differences still persists.

Yours sincerely,
/Dan Fors

[Wenjie]

I haven't figured out this issue yet. I may need more time to look into it.
Thanks for your patient.

Wenjie