Open amrit-poudel opened 7 years ago
Almost certainly this is due to inaccuracies in BZ integration. It seems to be extremely tricky to devise a universal BZ integration strategy that works well for all materials, frequencies, and evaluation points.
To check this, I would try plotting the LDOS integrand vs. Bloch vector [columns 8 (electric LDOS) or 9 (magnetic LDOS) vs. columns 6,7 (kx,ky) in the .byOmegakBloch
file] and looking at the shape of the integrand. My guess is that looking at the integrand will suggest why the numerical cubature is inaccurate.
Also, I have found that the BZ integration proceeds much more robustly for the case of a perfect metallic half space (see the PECPlate
geometries in the HalfSpaceLDOS
example), so I would encourage you to investigate that case too for comparison. Comparing results for PEC and non-PEC problems may shed light on optimal BZ integration strategies.
I would also suggest looking at higher frequencies. For a geometry with lengthscales on the order of 10s or 100s of nanometers, \omega=1 (in SCUFF units) is a very low frequency. What this means is that the pole in the BZ integrand at |k| = \omega/c lies very close to the origin of the Brillouin zone. (Indeed, for a square lattice with lattice constant L=10 nm, the Brillouin zone extends out to |k| = \pm 314 inverse microns, so the pole at \omega/c=1 inverse micron lies less than 1% of a reciprocal-lattice constant from the origin. This results in a very sharply-peaked integrand that is hard to evaluate by numerical cubature (although the --BZIMethod Polar2
strategy should be helpful).
Finally, N=20 is probably too high-resolution meshing unless you are at very high frequency or very small distance from the surface. This may be slowing down your calculations unnecessarily. For lattice constant L=10 nm I got good results all the way up to Omega=50 or so with just N=4.
Of course, for the purpose of testing BZ integration strategies you should just use --HalfSpace Aluminum
, which will bypass numerical SCUFF calculations altogether and use analytical formulas instead. (In this case the .scuffgeo
file is consulted only for the purpose of reading off the lattice constant, but the mesh files are not referenced).
For analytical calculations for a PEC half-space you can say either --HalfSpace PEC
or --GroundPlane
; these describe the same geometry but are computed in different ways, which I have tested are in agreement in all cases I considered.
Thank you for the update. I just wanted to add that the above calculation was done for L = 1 (scuffem units) which corresponds to 1 micron and \omega = 1 scuffem units which corresponds to 3e14 rad/s. I am not sure why this is a low frequency.
As for N=20, the reason I took such a fine mesh is because I was interested in accuracy of the calculation at z = 0.1 (scuffem units) which corresponds to 0.1 micron.
I just wanted to report back on the LDOS simulation I did after your most recent fixes of BZ integration.
I can now confirm that numerical and analytical LDOS computed by SCUFFEM are consistent with each other for the parameters I chose ( N = 20, \omega =1, Aluminum half space). Please find the attached plot below.
My only concern now is that theoretically (and physically) I expect LDOS of a metallic half space to fall of as 1/z^3 near the surface and 1/z^2 further away from the surface, where z = distance from the surface of the metal. However, this does not seem to be the case for numerical or theoretical results obtained from SCUFFEM.
Could you please clarify why LDOS computed by SCUFFEM does not behave like 1/z^3 for metallic half space?
I tried with finer mesh (N = 40) also, but I did not see any improvement on 1/z^3 behavior.
Thanks!![ldos_vs_z_1_report](https://cloud.githubusercontent.com/assets/12847183/17535192/b6eedc62-5e52-11e6-94b8-3b30c9556c76.jpg)