HomerReid / scuff-em

A comprehensive and full-featured computational physics suite for boundary-element analysis of electromagnetic scattering, fluctuation-induced phenomena (Casimir forces and radiative heat transfer), nanophotonics, RF device engineering, electrostatics, and more. Includes a core library with C++ and python APIs as well as many command-line applications.
http://www.homerreid.com/scuff-em
GNU General Public License v2.0
125 stars 50 forks source link

LDOS #100

Closed amrit-poudel closed 7 years ago

amrit-poudel commented 8 years ago

Hi Homer,

I am a bit confused by the LDOS output in Halfspace.byOmegakBloch and Halfspace.LDOS files (below).

Is the electric LDOS output in Halfspace.LDOS really the electric LDOS computed in the real space (integrated over "kx and ky". Do you integrate this internally? Since I only provided kx =0 and a few ky points, do you choose other kx and ky points yourself)? If so, why does it output so many values, even though there are only two separate distances "z" from the halfspace (I see that these values simply repeat, but I just wanted to confirm this from you)? In this simulation, I chose z = 0.1 and 1, with kx = 0 and ky spans from 0 to pi with a step size of 0.1 * pi (so there are total of 11 ky points). length scale is in micron (1e-6 m) and omega = 1 (which corresponds to 3e14 rad/s).

Please note that the formula for electric LDOS = \omega/pi/c^2 Tr [Im[G]], where G = electric dyadic Green's function. In 3D, The dimension of LDOS must be [time]/[length^3] or sec/m^3 in SI units.

What is the unit of LDOS output in Halfspace.LDOS file?

To obtain these results, I ran my simulation on a single core, just to avoid any inconsistency. Also, I am using the most recent version of scuff-em, downloaded from github.

Thank you for your help.

https://github.com/amrit-poudel/Boundary_Integral_Method/blob/master/Halfspace.LDOS https://github.com/amrit-poudel/Boundary_Integral_Method/blob/master/Halfspace.byOmegakBloch

HomerReid commented 8 years ago

Is the electric LDOS output in Halfspace.LDOS really the electric LDOS computed in the real space (integrated over "kx and ky". Do you integrate this internally?

Yes, that's right. The values reported in the .LDOS file are real-space LDOS values computed by integrating over the Brillouin zone. The sample points used by SCUFF-LDOS to evaluate this integral, and the function values at those points, are reported in the corresponding .byOmegakBloch file.

Since I only provided kx =0 and a few ky points, do you choose other kx and ky points yourself)?

I'm confused. How did you specify (kx,ky) points? If you got an .LDOS file, that means SCUFF-LDOS evaluated the Brillouin-zone integral internally, using (kx,ky) points that it selected automatically. Alternatively, if you used the --OmegakBlochFile option to give a list of (\omega, kx, ky) triples, then the code should not have produced a .LDOS file.

In the .byOmegakBlochFile you supplied above, the (kx,ky) values are changing from line to line, indicating that the code is internally selecting its own Brillouin-zone sample points.

What is the unit of LDOS output in Halfspace.LDOS file?

The components of the dyadic Green's functions (GE, GM) as reported in the .LDOS output file have units of (microns)^{-1}.

The electric and magnetic LDOS reported in the .LDOS output file are obtained by multiplying (Tr Im GE, GM) by \omega in units of c/(1 micron). So the units of these quantities are c/micron^2 = 3e14 (micron^{-3}) (rad/sec)^{-1}.

amrit-poudel commented 8 years ago

Thank you for taking your time to clarify these points. I truly appreciate it.

I ran the LDOS simulation using following arguments:

ARGS="" ARGS="${ARGS} --geometry ./ldos_inputs/Halfspace.scuffgeo" ARGS="${ARGS} --EPFile ./ldos_inputs/EPFile" ARGS="${ARGS} --OmegaFile ./ldos_inputs/OKBFile"

With these arguments, I got Halfspace.byOmegakBloch and Halfspace.LDOS files as my output. I did specify (kx,ky) values in OKBFile: https://github.com/amrit-poudel/Boundary_Integral_Method/blob/master/OKBFile

However, for some reason, the LDOS script did not seem to read my (kx, ky) values, which explains why the output .byOmegakBloch file has different kx,ky values and an additional .LDOS output file.

Could you please comment on why the OKBFile was not read by the LDOS script? Do I have to write everything as float? any tab issues, etc.?

HomerReid commented 8 years ago

You used the option --OmegaFile, which is used to specify only the temporal frequencies \omega. With this option, SCUFF-LDOS only looks at the first column of the input file and ignores all subsequent columns; for each \omega value, the code evaluates the full Brillouin-zone integral using automatically-chosen sample points.

If you want to bypass the internal Brillouin-zone integrator and request values of the Brillouin-zone integrand at specific k-points, you want to use the option --OmegakBlochFile. From what you have described, you should be able to change --OmegaFile to --OmegakBlochFile and leave everything else unchanged in the example you quoted.

amrit-poudel commented 8 years ago

Aside: I looked more carefully into .LDOS file (attached in my previous message). The LDOS values are repeating after every two lines. Since I only supplied two different values of z (which are 0.1 and 1, also clearly seen in.LDOS file), I do not expect .LDOS file to have more than 2 lines of data. However, it generates so many lines of repeated data. Could this be a (harmless) bug?

amrit-poudel commented 8 years ago

Thank you very much for your quick reply! That explains why my (kx, ky) values were not read. I will re-run my calculation with these corrections. Thanks.

amrit-poudel commented 8 years ago

Just wanted to confirm that ----OmegakBlochFile option works fine in terms of reading w, kx, ky triple values.

That said, I decided to do the test by using --OmegaFile option so that I don't have to supply kx, ky points. Here are the outputs from this test run:

https://github.com/amrit-poudel/Boundary_Integral_Method/blob/master/Halfspace.LDOS https://github.com/amrit-poudel/Boundary_Integral_Method/blob/master/Halfspace.byOmegakBloch

I have three questions:

(1) Why is LDOS negative in .LDOS file?

(2) Why does .LDOS file have trailing zeros on the far right ? I counted them and there are 38 trailing zeros.

(3) I supplied 4 observation points (different distances from the metallic half space) in EPFile. Therefore, I expect only 4 rows (and 43 columns) in .LDOS file. The reason why I expect only 4 rows is because this is a real space LDOS obtained by integrating out kx and ky. However, there are more than 4 rows in the .LDOS output. Could you please clarify why?

Thanks.

HomerReid commented 8 years ago

(1) Why is LDOS negative in .LDOS file?

Because this is only the scattering contribution to the LDOS. To get the total LDOS you have to add in the vacuum contributions as well.

(2) Why does .LDOS file have trailing zeros on the far right ? I counted them and there are 38 trailing zeros.

The value in column N+38 is the estimated error incurred by the numerical Brillouin-zone cubature in computing the output quantity in column N. This is only nonzero if you chose an integration algorithm that provides error estimates (the default does not). For example, if you say --BZIMethod CC --BZIOrder 0 you will get nonzero error estimates.

(3) I supplied 4 observation points (different distances from the metallic half space) in EPFile. Therefore, I expect only 4 rows (and 43 columns) in .LDOS file. The reason why I expect only 4 rows is because this is a real space LDOS obtained by integrating out kx and ky. However, there are more than 4 rows in the .LDOS output. Could you please clarify why?

My guess is that you prepared a input file to be passed to the --OmegakBlochFile option, but instead passed it to the --OmegaFile option. Presumably your file consists of many lines with the same \omega value but different (kx,ky) values. But when you specify --OmegaFile the code ignores everything but the first column on each line. If there are multiple lines specifying the same \omega value, the LDOS calculation at that frequency will be repeated.

The --OmegakBlochFile and --OmegaFile options are really designed to accept different types of input files:

amrit-poudel commented 8 years ago

Thank you for clarifying. Everything you said makes perfect sense and completely clears my confusion. I will report back with the comparison between numerical and analytical LDOS results from scuffem.

amrit-poudel commented 8 years ago

Here are the outputs from my test runs (from scuff LDOS executable):

https://github.com/amrit-poudel/Boundary_Integral_Method/blob/master/Halfspace.Aluminum.LDOS https://github.com/amrit-poudel/Boundary_Integral_Method/blob/master/Halfspace.LDOS

In both cases, I chose N = 20 in Square.geo file, \omega = 1.0 in OFile (file contains single \omega value) and z = 0.1, 1, 10, 100. I chose N= 20 to get accurate result all the way down to z = 0.1. However, that did not turn out to be the case. A plot comparing the two is attached.

ldos_vs_z_1

I could not get these two to agree even though I took N=20.

I will appreciate any help on how I can get accurate result.

Please note that I divided both numerical and analytical LDOS by c*a^2, where c = 3e8 m/s and a = 1e-6 m. This conversion to SI units is irrelevant (since i divide both numerical and analytical by the same factor) for the shake of comparison between numerical and analytical results.

HomerReid commented 8 years ago

The long story short here is that the problems are associated with the Brillouin-zone integration.

I am pretty sure that SCUFF-LDOS is computing the right values for each individual kBloch point, i.e. each sample of the integrand. But the integrand is badly behaved as a function of kBloch, with crazy death-defying blowups and plunges at points within the Brillouin zone. You can get a feel for this by plotting the data in the .byOmegakBlochFile -- plot column 8 (electric LDOS integrand) against columns 6 and 7 (x,y components of kBloch) and you will see wild variations of the integrand as a function of kBloch, which make the integral very difficult to evaluate numerically.

I am working on ways to get around this and will post some updates to the online documentation regarding Brillouin-zone integration shortly.

amrit-poudel commented 8 years ago

thanks for looking into this. I appreciate it.

For completeness, I have also included .byOmegakBloch files here:

https://github.com/amrit-poudel/Boundary_Integral_Method/blob/master/Halfspace.byOmegakBloch

https://github.com/amrit-poudel/Boundary_Integral_Method/blob/master/Halfspace.Aluminum.byOmegakBloch

amrit-poudel commented 8 years ago

Hi Homer, I was just wondering whether you got a chance to update this. Thanks.

I have found that adaptive quadrature based on a Gauss-Kronrod pair implemented in Matlab's "quadgk" function is better suited for k integration of the Green's function. I have not tested other algorithms.

HomerReid commented 8 years ago

Thanks for checking in---I am working on an LDOS tutorial, and it will be ready by the middle of next week. If not, please pester and bug me (seriously).

I have recently updated the available options for Brillouin-zone integration and symmetry factors, and I have updated the documentation on that:

http://homerreid.github.io/scuff-em-documentation/reference/BrillouinZoneIntegration/

The long story short is that it seems to be too hard for the code to figure out automatically how to sample the Brillouin-zone integrand. There is too much variation in the shape of the integrand depending on the parameters of the problem (material, frequency, evaluation points, etc.)

Instead, the way forward will probably be for users to be involved in designing the right integration scheme. This will probably involve users specifying a list of kBloch points (i.e. an --OmegakBlochFile) and asking SCUFF-LDOS to output the integrand at each point, then using my separate SCUFF-INTEGRATE utility code to integrate the data. Then, as necessary, the user can go back and ask SCUFF-LDOS for integrand samples over a finer grid in some narrow region of the Brillouin zone, which can then be combined with the samples previously obtained to get a more accurate estimate of the integral.

I am thinking about how to make this process most user-friendly, and I would appreciate suggestions from users.

HomerReid commented 7 years ago

OK, the tutorial is up:

http://homerreid.github.io/scuff-em-documentation/examples/HalfSpaceLDOS/HalfSpaceLDOS/

Everything about this has taken much longer than I expected and there is still more to do, but there's certainly enough there to get started.

I think we should close this issue at this point. If you agree, please close it. If you have specific questions regarding the new LDOS tutorial or LDOS calculations in general, open new issues for those.

amrit-poudel commented 7 years ago

Thank you very much for the update. I will try to reproduce your results and close this issue.

Just a quick note: On your tutorial, under the "Green's function and LDOS" section \Delta \rho^{E,M} does not seem to have the correct unit. I think you are missing a factor of 1/c^2, assuming G_E has units of 1/meter.

amrit-poudel commented 7 years ago

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?

Thanks!

ldos_vs_z_1_report

HomerReid commented 7 years ago

Thanks for looking into this. I have a response, but let me ask to you re-post your comment as a new issue. It is really a different question from the rest of the content in this thread, and it is a question that deserves its own thread, which can generate useful fodder for discussion going forward, particularly since (as of recently) we are in something of a new era regarding the accuracy of LDOS calculations in SCUFF. In contrast, I think the original content of this thread has been addressed adequately and this this issue should now be considered closed.

I don't want to be overly pedantic and particular about precise classification of issues, but I want to make sure we maintain a sense of progress, so that old issues are seen to have been resolved even if they birth new issues in their wake.

amrit-poudel commented 7 years ago

Sure, I will open a new issue. Thanks.