01/2021: updated to pyGDM v1.1+
In a first example, we will search for the geometry of a gold nano-structure which leads to maximum scattering at a specific wavelength and polarization.
Load the modules¶
from pyGDM2 import core from pyGDM2 import propagators from pyGDM2 import structures from pyGDM2 import materials from pyGDM2 import fields from pyGDM2.EO.core import run_eo from pyGDM2.EO.problems import ProblemScat from pyGDM2.EO.models import RectangularAntenna
Setup (1): pyGDM simulation¶
In a first step, we configure the pyGDM simulation which we will use. This is the same as for any other simulation, except that we will leave the geometry blank and use an empty list instead. The structure geometry will be defined below by the EO model and will be varied during the optimization.
## ---------- Setup structure mesh = 'cube' step = 15 material = materials.gold() # material: gold ## --- Empty dummy-geometry, will be replaced on run-time by EO trial geometries geometry =  struct = structures.struct(step, geometry, material) ## ---------- Setup incident field field_generator = fields.planewave # planwave excitation kwargs = dict(theta = [0.0]) # target lin. polarization angle wavelengths =  # target wavelength efield = fields.efield(field_generator, wavelengths=wavelengths, kwargs=kwargs) ## ---------- environment dyads = propagators.DyadsQuasistatic123(n1=1, n2=1) ## ---------- Simulation initialization sim = core.simulation(struct, efield, dyads)
structure initialization - automatic mesh detection: cube structure initialization - consistency check: 0/0 dipoles valid
/home/hans/.local/lib/python3.8/site-packages/pyGDM2/tools.py:817: UserWarning: Empty structure. Setting mesh to 'cubic'. warnings.warn("Empty structure. Setting mesh to 'cubic'.") /home/hans/.local/lib/python3.8/site-packages/pyGDM2/structures.py:183: UserWarning: Emtpy structure geometry. warnings.warn("Emtpy structure geometry.")
Note: We can safely ignore the warning about the empty structure: The geometry will be handled by the EO model.
Setup (2): Evolutionary optimization setup¶
In the second step we now configure the evolutionary optimization: We will setup the structure model and the optimization problem as well as the optimization algorithm.
The model will determine the structure geometry as function of certain free parameters which will be subjet to the optimization. For demonstration, we will choose a very simple model: A planar cuboid of rectangular footprint, where the free parameters are only its width and length (the height shall be fixed).
The problem will specify the optimization objective, hence a certain optical property to be maximized (or minimized). In our example, we want to maximize the scattering from the plasmonic nano-structure.
As algorithm we use the “jde” differential evolution algorithm, implemented in pyGMO/paGMO.
## --- structure model: Rectangular planar antenna of fixed height limits_W = [2, 20] # units of "step" limits_L = [2, 20] # units of "step" limits_pos = [-1, 1] # units of nm --> effectively no shift of the struct. allowed height = 3 # units of "step" model = RectangularAntenna(sim, limits_W, limits_L, limits_pos, height) ## --- optimization problem: Scattering opt_target = 'Qscat' # 'Qscat' --> scat. efficiency problem = ProblemScat(model, opt_target=opt_target) ## --- filename to save results results_filename = 'eo_Qscat.eo' ## --- size of population population = 25 # Nr of individuals ## --- stop criteria max_time = 60 # seconds max_iter = 20 # max. iterations max_nonsuccess = 5 # max. consecutive iterations without improvement ## --- other config generations = 1 # generations to evolve between status reports plot_interval = 1 # plot each N improvements save_all_generations = False ## Use algorithm "sade" (jDE variant, a self-adaptive form of differential evolution) import pygmo as pg algorithm = pg.sade algorithm_kwargs = dict() # optional kwargs passed to the algorithm
Rectangular Antenna optimziation model: Note that this simple model is rather intended for testing and demonstration purposes.
Run the optimization¶
Now let’s run this optimization:
eo_dict = run_eo(problem, population=population, algorithm=algorithm, plot_interval=plot_interval, generations=generations, max_time=max_time, max_iter=max_iter, max_nonsuccess=max_nonsuccess, filename=results_filename)
---------------------------------------------- Starting new optimization ----------------------------------------------
/home/hans/.local/lib/python3.8/site-packages/pyGDM2/tools.py:835: UserWarning: Mesh not detected, falling back to 'cubic'. warnings.warn("Mesh not detected, falling back to 'cubic'.") /home/hans/.local/lib/python3.8/site-packages/numba/core/dispatcher.py:237: UserWarning: Numba extension module 'numba_scipy' failed to load due to 'ValueError(No function '__pyx_fuse_0pdtr' found in __pyx_capi__ of 'scipy.special.cython_special')'. entrypoints.init_all() /home/hans/.local/lib/python3.8/site-packages/pyGDM2/tools.py:835: UserWarning: Mesh not detected, falling back to 'cubic'. warnings.warn("Mesh not detected, falling back to 'cubic'.")
iter # 1, time: 27.2s, progress # 1, f_evals: 50
/home/hans/.local/lib/python3.8/site-packages/pyGDM2/EO/models.py:155: UserWarning: 'models.BaseModel.plot_structure' not re-implemented! Using `pyGDM2.visu.structure`. warnings.warn("'models.BaseModel.plot_structure' not re-implemented! Using `pyGDM2.visu.structure`.") /home/hans/.local/lib/python3.8/site-packages/pyGDM2/visu.py:49: UserWarning: 3D data. Falling back to XY projection... warnings.warn("3D data. Falling back to XY projection...")
/home/hans/.local/lib/python3.8/site-packages/pyGDM2/tools.py:835: UserWarning: Mesh not detected, falling back to 'cubic'. warnings.warn("Mesh not detected, falling back to 'cubic'.")
- champion fitness: [-14.484] iter # 2, time: 49.9s, progress # 2, f_evals: 75
/home/hans/.local/lib/python3.8/site-packages/pyGDM2/EO/models.py:155: UserWarning: 'models.BaseModel.plot_structure' not re-implemented! Using `pyGDM2.visu.structure`. warnings.warn("'models.BaseModel.plot_structure' not re-implemented! Using `pyGDM2.visu.structure`.")
- champion fitness: [-26.052] iter # 3, time: 63.0s, progress # 2, f_evals: 100(non-success: 1) -------- timelimit reached
The output will look something like this.
Note: In this run, the optimum structure was already found in the third generation, it may even happen to be in the randomized initial population, since there are only 2 free parameters which are furthermore very constrained. For more complex structure models, the probability of this to happen will of course be very low.
Post-processing of optimum solution¶
Now we want to calculate the scattering spectrum for the optimum solution. We will load the simulation from the files, generated by do_eo. Then we create a new simulation with a spectrum of wavelengths.
## --- load additional modules from pyGDM2.EO.tools import get_best_candidate from pyGDM2 import linear from pyGDM2 import tools from pyGDM2 import visu import copy import numpy as np import matplotlib.pyplot as plt #============================================================================== # load the final fittest candidate (=optimum geometry) #============================================================================== ## --- optimization results file results_filename = 'eo_Qscat.eo' sim = get_best_candidate(results_filename, iteration=-1, verbose=True)
Best candidate after 3 iterations (with 1 improvements): fitness = ['-26.052'] Testing: recalculating fitness...Done. Everything OK.
Note: We loaded the best candidate from iteration Nr “-1”., which is the last iteration of the evolution (“-2” would be second last and so on; positive numbers starting from “0” can be used as well, this is just python indexing.).
#============================================================================== # setup new simulation to calculate spectrum #============================================================================== ## --- structure struct = copy.deepcopy(sim.struct) ## --- incident field field_generator = fields.planewave # planwave excitation wavelengths = np.arange(600, 1410, 30) # spectrum kwargs = dict(theta = [0.0, 90.0]) # 0 / 90 deg polarizations efield = fields.efield(field_generator, wavelengths=wavelengths, kwargs=kwargs) ## --- env. dyads = propagators.DyadsQuasistatic123(n1=1, n2=1) ## --- simulation sim_spectrum = core.simulation(struct, efield, dyads) #============================================================================== # run simulation for the spectrum #============================================================================== core.scatter(sim_spectrum, verbose=False) ## --- calculate the spectrum for X and Y polarization wl, spec_ext0 = tools.calculate_spectrum(sim_spectrum, 0, linear.extinct) asca0 = spec_ext0.T wl, spec_ext90 = tools.calculate_spectrum(sim_spectrum, 1, linear.extinct) asca90 = spec_ext90.T geom_cs = tools.get_geometric_cross_section(sim_spectrum)
Plot the scattering spectrum for the best solution¶
## --- plot plt.figure(figsize=(10,5)) ## --- spectra plt.subplot2grid((1,5), (0,0), colspan=3) plt.title("scattering spectrum") plt.plot(wl, asca0/geom_cs, label="0deg") plt.plot(wl, asca90/geom_cs, label="90deg") plt.legend(loc='best', fontsize=10) plt.xlabel("wavelength (nm)") plt.ylabel("Q_scat") ## --- structure plt.subplot2grid((1,5), (0,3), colspan=2, aspect="equal") plt.title('structure geometry') visu.structure(sim_spectrum, show=False) plt.show()
Indeed, a structure was found which has a strong plasmon resonance at the target wavelength of 1000nm.