My Community

Dear VASP team and community,

I am currently working to use the new XANES functionality in VASP 6.2.1 to calculate core-loss spectra for a project. I am technically interested in ELNES, but they should essentially be the same.

First, however, I am trying to reproduce some of the results in the following paper:
W. Zhou et al, "Direct Determination of the Chemical Bonding of Individual Impurities in Graphene," Phys Rev Lett 109, 206803 (2012).
In particular, I am attempting to calculate the spectra for the defect in which an Si atom fills a carbon divacancy (Si-4C), this would be the L2/3 edge of silicon and the composite theoretical spectra is shown in Fig. 3(b) of the paper. In the SI, the breakdown into the in-plane and out-of-plane contribution from the dielectric function are shown in S2.

Now, I've used a slightly larger than this paper. However, when I calculate the L2/3 edge for the silicon, I obtain a spectra that is not anisotropic. That is the core dielectric tensor has the same value for the xx, yy, and zz at each energy. In the attached plot the two lines are exactly on top of each other - the plot was re-broadened from the raw output, but this is not the issue. I am not entirely sure what the issue here is as I expect the spectra to not be isotropic between the in-plane and out-of-plane components. The peaks obtained via VASP match approximately to those obtained before but their intensities are quite different which suggest to me that some sort of averaging effect is being introduced. If instead I calculate the carbon K-edge for pristine graphene (also re-broadened in postprocessing), I obtain an anisotropic tensor as one would expect. I can then average them based on the experimental beam parameters. It is not clear to me if I am missing a setting in the INCAR or if something is happening inside the code though. Any help would be much appreciated. I have included my input files that were used to perform the calculation as well as some images. Please let me know if additional information is needed to help.

Thank you for your help,
Andrew

Dear Andrew,

Thanks for the detailed report.
Indeed something seems to be wrong. We will look into it
I have one request and a suggestion:

1. Do you think you could try to reproduce this issue on a system with fewer atoms?
Same defect but a smaller graphene supercell.
This would help us a lot to track down the problem.

2. Could you try to run the calculation with ISYM=-1 and see if the result changes?
(This could be done on the system with fewer atoms)

Henrique,

Thank you for looking into this to help resolve the issue. I apologize for the delay in providing the additional calculation due to the American holiday.

I re-ran the calculation using a 5x5 graphene supercell with the same defect of a silicon atom sitting inside a carbon divacancy (but unrelaxed - the only changes would be light in-plane relaxations). For me, this calculation can run on a single node (48 cores) in under 5 minutes. I also set ISYM = -1 as requested. The spectra below is the result (unlike before, no additional broadening was employed - only the setting in VASP):

As shown, the issue still remains with ISYM = -1.

I have provided the input/output files in the attached tar file.

Thanks again for your help,
Andrew

In your calculation you were using the normal Si potentials. These potentials only have s and p projectors. Since you excite from Si 2p states, the p->p transiotions are forbidden by the transition rules. So only p->s transitions are left which are always isotropic. You would need p->d transitions. These can be obtained if you use the Si GW potentials. It is generally a good rule to use the GW for XAS calculations.

Another problem in your calculation is that Si in C is a metal. The current code only works properly for insulators. We will fix this in the upcoming VASP 6.3 release.
If you need the fix now, you would have to do the following changes in the core_con_mat.F and recompile:
-) Remove the following line:
IF (W%FERWE(I_CONBANDS,KPOINT,ISPIN) >= 0.5_q) cycle

-) Change the following lines:
EPSILON_HELP = PREFACTOR/FREQUENCY**2.0_q * &
WDES%WTKPT(KPOINT)*DELTAFUN*&
MECC1*CONJG(MECC2)*CH_AMPLIFICATION
to:
EPSILON_HELP = (1.0_q-W%FERWE(I_CONBANDS,KPOINT,ISPIN)) * &
PREFACTOR/FREQUENCY**2.0_q * &
WDES%WTKPT(KPOINT)*DELTAFUN * &
MECC1*CONJG(MECC2)*CH_AMPLIFICATION