Implementation of wavelength scan via python

[Justin Koenig]

Hello,

as per issue 35 (https://code.google.com/archive/p/a-dda/issues/35) i know it is not possible to scan the wavelength in-software.

has anyone got an implementation in python which inputs some refractive index files (like those from https://refractiveindex.info/) and lets me use the software as usual, just over a range of wavelengths? i am just intersted in Q_sca and Q_abs.

I have only a weak grasp of python and little time to make an implementation, perhaps someone already has this. if not, I will keep this topic open for discussion and I will show my own progress in writing a code for this.

Best

Justin Koenig

[Maxim Yurkin]

Dear  Justin,

I have just added a comment to this issue, describing the existing options (that I am aware of). This should be sufficient for the simple task that you describe. But let me know, if you have further questions.

Maxim.

P.S. This answer has been forwarded to your e-mail address for your convenience. However, if you want to continue the discussion please reply to the group's e-mail. You are also advised to check the corresponding discussion thread at http://groups.google.com/group/adda-discuss and/or subscribe to it to receive automatic notifications, since other answers/comments will not be necessarily forwarded to your e-mail. 

--
You received this message because you are subscribed to the Google Groups "ADDA questions and answers" group.
To unsubscribe from this group and stop receiving emails from it, send an email to adda-discuss...@googlegroups.com.
To view this discussion visit https://groups.google.com/d/msgid/adda-discuss/f3f3c7f8-3e09-41f1-b0b2-5c5260747e99n%40googlegroups.com.

[Justin Koenig]

Dear Maxim,

Here you can find my implementation (to which i shall upon request make a github repository), which is by no means final or well-tested. It has been significantly altered from https://github.com/baptiste/adda/wiki/wrapper_primer and adapted for python. I make use of the refractiveindex.info database utilizing https://github.com/toftul/refractiveindex.

(Python 3.13)

import subprocess

import re

import numpy as np

import csv

import time

from refractiveindex import RefractiveIndexMaterial as nMat

def adda_spectrum(shape="ellipsoid",

                  wavelength=500, #wavelength of the incoming light in nm

                  radius=20, #Radius of the particle in nm

                  nPart=1.5 + 0.2j, #RI of the particles

                  euler = (0, 0, 0), #euler angles

                  granulRI = 1.4 + 0.1j, #RI of the internal granules

                  nMedium=1.333, #RI of the medium

                  dpl=None, #number of dipoles

                  test=True, #if test, do not run command

                  verbose=True, #if verbose, print command

                  **kwargs):

    if dpl is None:

        dpl = int(np.ceil(min(50, 20 * abs(nPart)))) #automatic scaling of dpl in case no argument is given

    cmd = (

        fr"C:\Users\innoFSPEC\Desktop\Amsterdam_DDA\win64_executables\adda -shape {shape} " #TODO make independent of path and implement other shapes

        f"-orient {euler[0]} {euler[1]} {euler[2]} "

        f"-lambda {wavelength * 1e-3 / nMedium} "

        f"-dpl {dpl} "

        f"-size {2 * round(int(radius) * 1e-3, 5)} "

        f"-m {np.real(nPart) / nMedium} {np.imag(nPart) / nMedium} {np.real(granulRI)/nMedium} {np.imag(granulRI)/nMedium} " #TODO: currently, granul is required or error occurs. make independent, make this row a variable and insert it as an f-string

    )

    # append additional kwargs

    if kwargs:

        extras = " ".join(f"-{k} {v}" for k, v in kwargs.items())

        cmd += extras

    if verbose:

        print("Command:", cmd)

    if test:

        Cext = 0

        Cabs = 0

        Csca = 0

        return Cext, Cabs, Csca, cmd

    # run command

    start_time = time.time()

    result = subprocess.run(cmd, shell=True, capture_output=True, text=True) #automatically waits for loop to complete

    end_time = time.time()

    elapsed_time = end_time - start_time

    print(f'Iteration done @ wavelength = {wavelength} nm. Time taken = {round(elapsed_time,1)}s ')

    output = result.stdout.splitlines()

    # extract numerical values

    def extract_vals(keyword):

        line = next((l for l in output if keyword in l), "")

        nums = re.findall(r"[-+]?\d*\.\d+|\d+", line)

        return list(map(float, nums))

    Cext = extract_vals("Cext")

    Cabs = extract_vals("Cabs")

    Csca = [Cext[i] - Cabs[i] for i in range(len(Cext))]

    return Cext[0], Cabs[0], Csca[0], cmd

# example usage

MedMat = nMat('main', 'H2O', 'Daimon-21.5C')

PartMat = nMat('main', 'SiO2', 'Arosa')

GranulMat = nMat('main', 'Ag', 'Ferrera-298K')

resultFile = open(fr'C:\Users\innoFSPEC\Desktop\Python programs\OutPut\output.csv', 'w')

writer = csv.writer(resultFile, delimiter='\t')

writer.writerow(['Wavelength /nm', 'Cext', 'Cabs', 'Csca'])

for wl in range(400, 1001, 10):

    nMed = MedMat.get_refractive_index(wl)

    nPart = PartMat.get_refractive_index(wl)

    nGranulReal = GranulMat.get_refractive_index(wl)

    nGranulIm = GranulMat.get_extinction_coefficient(wl)

    TotalRi = f'{nPart} 0 {nGranulReal} {nGranulIm}'

    Cext, Cabs, Csca, cmd = adda_spectrum(test=False, shape ="sphere", wavelength = wl, radius = 175, nPart = nPart, granulRI= complex(nGranulReal,nGranulIm), nMedium = nMed, dpl = 100, verbose = True, granul= '0.10 0.08', asym ='')

    writer.writerow([wl, Cext, Cabs, Csca])

resultFile.close()

However, with the giiven command (which reads as  "C:\Users\innoFSPEC\Desktop\Amsterdam_DDA\win64_executables\adda -shape sphere -orient 0 0 0 -lambda 0.2977493924382135 -dpl 100 -size 0.35 -m 1.094320021536666 0.0 0.06779811872466268 1.5473168450068826 -granul 0.10 0.08 -asym" (of course keeping in mind that refractive index of medium, particle, and wavelength constantly change), i've stumbled upon a problem: the result looks very much rather jagged, and i'm not sure what causes this.

The refractive indices themselves are smoothly interpolated from refractiveindex.info. Perhaps my granules are much too small? In your paper https://opg.optica.org/oe/fulltext.cfm?uri=oe-15-25-16561, the granules were 10x as large as mine, albeit the dipoles were very much lower. I'm not exactly sure where to best iterate from here without spending a lot of time. At 100 dipoles per lambda, each dipole is approx. 3nm large and the granules are 8nm large. i really hope this is not the problem, because one iteration at the lowest wavelength takes something like 15 minutes to calculate on my current PC.

Best,

Justin

[Justin Koenig]

I've just seen that i only generate 9 granules per run - they are 80 nm large, not 8 nm. let me re-run and get back to you ASAP.

[Justin Koenig]

This seems more physically acceptable now. The jitters remain - I assume because the granules are randomly generated upon each iteration. Handy but perhaps not always something I would like. Any way to set the same seed for each generation? I was unable to find anything related to that in the manual.

[Maxim Yurkin]

Hi Justin,

Thanks for sharing the script. The scanning over wavelength seems to work perfectly, while the whole following
discussion is related to the specifics of your scattering problem (that is small size relative to the wavelength and
granules). You have almost solved the problem yourself while I was writing the reply. So some of it is probably
redundant, but the rest is still useful. For the question about seed, see the end of this reply.

# dpl vs. grid

For nanoparticles, it is not recommended to use dpl, it's better to rather specify grid directly based on the complexity
of the shape. In your case, I would recommend to start from 1 voxel (dipole) per granule - on the one hand, that is
grossly inaccurate for a single granule (note also that it is generally shifted with respect to the voxel), on the other
hand, the corresponding errors will average out due to large number of granules. Still, further testing with refining
discretization (try 2, 3, 4 times smaller voxels) is recommended. However, that would be complicated due to randomness
of granule positions, as described below.

Note also that using nm as unit of all sizes is perfectly fine - this will make it easier to debug the script.

# Randomness of granule positions
Concerning the oscillations, you have correctly related them to the granules. The current script regenerates the
granules for each run, so each wavelength has slightly different particle shape. The advantage is that you get (for
free) the estimation of the variability due to the granule positions. The disadvantage is that it looks ugly and, most
probably, that is not what you expect, when asking for the spectrum.

In principle, the most robust approach is to perform multiple simulations for each wavelength (varying granule
positions) and then average. In this scenario, current implementation can even be more effective, since it can trade
some spectral resolution (which is probably good enough) for reduced number of random samples (to get the same error
estimate). But it is much more complicated to do it properly with reliable error control.

So a more straightforward (and, thus, less error-prone) approach is to use exactly the same shape for the whole
spectrum. For that use `-save_geom` in the first run (or a preliminary run with `-save_geom -prognosis`) to produce the
shape, then extract it from the run folder and use `-shape read ...` in all further runs (and using neither `-dpl ...`
nor `-grid ...`). The only drawback is that this scenario will necessary fix the size of the discretization grid for all
wavelength. This is fine for nanoparticles as described above (i.e. if -grid  was originally used).

# Additional comments about granules

For larger particles (and sufficiently large granules) the above strategy still works (if the shape is generated for the
shortest wavelength, either using -grid or -dpl), but will become progressively inefficient (computationally) when the
wavelength is increased further. In the latter case increasing the voxel size is generally recommended (can be done
keeping the errors small). When the dynamic range of wavelength is more than two, some improvement can be obtained by
using -jagged option. Specifically, for a range of [a,b] (b>2a], one can generate the shape for wavelength of 2a and use
it for [2a, b], while the same shape with `-jagged 2` for the range [a,2a). The same idea can be further refined with
larger arguments to jagged (if the dynamic range is even larger).

Finally, I have just created an issue to make the granule generator more reproducible - see
https://github.com/adda-team/adda/issues/346. When it is implemented (it is open to contributions), the above saving and
reading of shape files can be replaced by passing the random seed. It also describes a dirty fix for such reproducibility.

Maxim.

[Justin Koenig]

Helllo Maxim,

I have since created a github repository for automation and multithreading. For anyone interested in working on the further development, you may find the implementation at https://github.com/JKoenig-spec/pyADDA-spectrum

keep in mind it creates the results folders directly in the source foulder, i have yet to rectify this, and saving geometry is also not yet supported, and it is a very manual setup i.e. the source code must be written to and executed directly. I am not a coder by any means, feel free to contribute.

Best

Justin.

[Maxim Yurkin]

Thanks a lot, Justin. Could you please change AmsterdamDDA to ADDA in the description? I prefer to have a single official name of the code to avoid confusion.

Maxim.

--
You received this message because you are subscribed to the Google Groups "ADDA questions and answers" group.
To unsubscribe from this group and stop receiving emails from it, send an email to adda-discuss...@googlegroups.com.

To view this discussion visit https://groups.google.com/d/msgid/adda-discuss/2db55644-9105-4c50-90a8-1edf220b8d3dn%40googlegroups.com.