finite_systems.md

   1 Using QPMS library for simulating finite systems
   2 ================================================
   3
   4 *** This tutorial is partly obsolete, the interpolators are no longer the first choice of getting the T-matrices. ***
   5
   6 The main C API for finite systems is defined in [scatsystem.h][], and the
   7 most relevant parts are wrapped into python modules. The central data structure
   8 defining the system of scatterers is [qpms_scatsys_t][],
   9 which holds information about particle positions and their T-matrices
  10 (provided by user) and about the symmetries of the system. Specifically, it
  11 keeps track about the symmetry group and how the particles transform
  12 under the symmetry operations.
  13
  14
  15 SVD of a finite symmetric system of scatterers
  16 ----------------------------------------------
  17
  18 Let's have look how thinks are done on a small python script.
  19 The following script is located in `misc/201903_finiterectlat_AaroBEC.py`.
  20
  21 ```{.py}
  22 #!/usr/bin/env python
  23 from qpms import Particle, CTMatrix, BaseSpec, FinitePointGroup, ScatteringSystem, TMatrixInterpolator, eV, hbar, c
  24 from qpms.symmetries import point_group_info
  25 import numpy as np
  26 import os
  27 import sys
  28 nm = 1e-9
  29
  30 sym = FinitePointGroup(point_group_info['D2h'])
  31 bspec = BaseSpec(lMax = 2)
  32 tmfile = '/m/phys/project/qd/Marek/tmatrix-experiments/Cylinder/AaroBEC/cylinder_50nm_lMax4_cleaned.TMatrix'
  33 #outputdatadir = '/home/necadam1/wrkdir/AaroBECfinite_new'
  34 outputdatadir = '/u/46/necadam1/unix/project/AaroBECfinite_new'
  35 os.makedirs(outputdatadir, exist_ok = True)
  36 interp = TMatrixInterpolator(tmfile, bspec, symmetrise = sym, atol = 1e-8)
  37 # There is only one t-matrix in the system for each frequency. We initialize the matrix with the lowest frequency data.
  38 # Later, we can replace it using the tmatrix[...] = interp(freq) and s.update_tmatrices NOT YET; TODO
  39
  40 omega = float(sys.argv[3]) * eV/hbar
  41 sv_threshold = float(sys.argv[4])
  42
  43 # Now place the particles and set background index.
  44 px = 571*nm; py = 621*nm
  45 n = 1.51
  46 Nx = int(sys.argv[1])
  47 Ny = int(sys.argv[2])
  48
  49 orig_x = (np.arange(Nx/2) + (0 if (Nx % 2) else .5)) * px
  50 orig_y = (np.arange(Ny/2) + (0 if (Ny % 2) else .5)) * py
  51
  52 orig_xy = np.stack(np.meshgrid(orig_x, orig_y), axis = -1)
  53
  54 tmatrix = interp(omega)
  55 particles = [Particle(orig_xy[i], tmatrix) for i in np.ndindex(orig_xy.shape[:-1])]
  56
  57
  58 ss = ScatteringSystem(particles, sym)
  59
  60
  61 k = n * omega / c
  62
  63
  64 for iri in range(ss.nirreps):
  65     mm_iri = ss.modeproblem_matrix_packed(k, iri)
  66     U, S, Vh = np.linalg.svd(mm_iri)
  67     print(iri, ss.irrep_names[iri], S[-1])
  68     starti = max(0,len(S) - np.searchsorted(S[::-1], sv_threshold, side='left')-1)
  69     np.savez(os.path.join(outputdatadir, 'Nx%d_Ny%d_%geV_ir%d.npz'%(Nx, Ny, omega/eV*hbar, iri)),
  70         S=S[starti:], omega=omega, Vh = Vh[starti:], iri=iri, Nx = Nx, Ny= Ny )
  71     # Don't forget to conjugate Vh before transforming it to the full vector!
  72 ```
  73
  74 Let's have a look at the imports.
  75
  76 ```{.py}
  77 from qpms import Particle, CTMatrix, BaseSpec, FinitePointGroup, ScatteringSystem, TMatrixInterpolator, eV, hbar, c
  78 from qpms.symmetries import point_group_info
  79 ```
  80
  81  * `Particle` is a wrapper over the C structure `qpms_particle_t`,
  82    containing information about particle position and T-matrix.
  83  * `CTMatrix` is a wrapper over the C structure `qpms_tmatrix_t`,
  84    containing a T-matrix.
  85  * `BaseSpec` is a wrapper over the C structure `qpms_vswf_set_spec_t`,
  86    defining with which subset of VSWFs we are working with and how their
  87    respective coefficients are ordered in memory. Typically, this
  88    just means having all electric and magnetic VSWFs up to a given multipole
  89    order `lMax` in the "standard" ordering, but other ways are possible.
  90    Note that different `Particle`s (or, more specifically, `CTMatrix`es)
  91    can have different `BaseSpec`s and happily coexist in the same
  92    `ScatteringSystem`.
  93    This makes sense if the system contains particles with different sizes,
  94    where the larger particles need cutoff at higher multipole orders.
  95  * `FinitePointGroup` is a wrapper over the C structure `qpms_finite_group_t`
  96    containing info about a 3D point group and its representation. Its contents
  97    are currently *not* generated using C code. Rather, it is populated using a
  98    `SVWFPointGroupInfo` instance from the
  99    `point_group_info` python dictionary, which uses sympy to generate the group
 100    and its representation from generators and some metadata.
 101  * `ScatteringSystem` is a wrapper over the C structure `qpms_scatsys_t`,
 102    mentioned earlier, containing info about the whole structure.
 103  * `TMatrixInterpolator` in a wrapper over the C structure
 104    `qpms_tmatrix_interpolator_t`  which contains tabulated T-matrices
 105    (calculated e.g. using [`scuff-tmatrix`][scuff-tmatrix])
 106    and generates frequency-interpolated T-matrices based on these.
 107  * `eV`, `hbar`, `c` are numerical constants with rather obvious meanings.
 108
 109 Let's go on:
 110 ```{.py}
 111 sym = FinitePointGroup(point_group_info['D2h'])
 112 bspec = BaseSpec(lMax = 2)
 113 tmfile = '/m/phys/project/qd/Marek/tmatrix-experiments/Cylinder/AaroBEC/cylinder_50nm_lMax4_cleaned.TMatrix'
 114 ...
 115 interp = TMatrixInterpolator(tmfile, bspec, symmetrise = sym, atol = 1e-8)
 116 ```
 117
 118 The `D2h` choice indicates that our system will have mirror symmetries along
 119 the *xy*, *xz* and *yz* axes. Using the `BaseSpec` with the standard
 120 constructor with `lMax = 2` we declare that we include all the VSWFs up to
 121 quadrupole order. Next, we create a `TMatrixInterpolator` based on a file
 122 created by `scuff-tmatrix`. We force the symmetrisation of the T-matrices with
 123 the same point group as the overall system symmetry in order to eliminate the
 124 possible asymmetries caused by the used mesh. The `atol` parameter just says
 125 that if the absolute value of a given T-matrix element is smaller than the
 126 `atol` value, it is set to zero.
 127
 128 ```{.py}
 129 omega = float(sys.argv[3]) * eV/hbar
 130 sv_threshold = float(sys.argv[4])
 131
 132 # Now place the particles and set background index.
 133 px = 571*nm; py = 621*nm
 134 n = 1.51
 135 Nx = int(sys.argv[1])
 136 Ny = int(sys.argv[2])
 137
 138 orig_x = (np.arange(Nx/2) + (0 if (Nx % 2) else .5)) * px
 139 orig_y = (np.arange(Ny/2) + (0 if (Ny % 2) else .5)) * py
 140
 141 orig_xy = np.stack(np.meshgrid(orig_x, orig_y), axis = -1)
 142
 143 tmatrix = interp(omega)
 144 particles = [Particle(orig_xy[i], tmatrix) for i in np.ndindex(orig_xy.shape[:-1])]
 145
 146 ss = ScatteringSystem(particles, sym)
 147 ```
 148 This chunk sets the light frequency and array size based on a command
 149 line argument. Then it generates a list of particles covering a quarter
 150 of a rectangular array. Finally, these particles are used to generate
 151 the final scattering system – the rest of the particles is generated
 152 automatically to satisfy the specified system symmetry.
 153
 154 ```{.py}
 155 for iri in range(ss.nirreps):
 156     mm_iri = ss.modeproblem_matrix_packed(k, iri)
 157     U, S, Vh = np.linalg.svd(mm_iri)
 158     print(iri, ss.irrep_names[iri], S[-1])
 159     starti = max(0,len(S) - np.searchsorted(S[::-1], sv_threshold, side='left')-1)
 160     np.savez(os.path.join(outputdatadir, 'Nx%d_Ny%d_%geV_ir%d.npz'%(Nx, Ny, omega/eV*hbar, iri)),
 161         S=S[starti:], omega=omega, Vh = Vh[starti:], iri=iri, Nx = Nx, Ny= Ny )
 162 ```
 163 The last part iterates over the irreducible representations of the systems.
 164 It generates scattering problem LHS (TODO ref) matrix reduced (projected)
 165 onto each irrep, and performs SVD on that reduced matrix,
 166 and saving the lowest singular values (or all singular values smaller than
 167 `sv_threshold`) together with their respective singular vectors to files.
 168
 169 The singular vectors corresponding to zero singular values represent the
 170 "modes" of the finite array.
 171
 172
 173 Analysing the results
 174 ---------------------
 175
 176 *TODO analyzing the resulting files.*
 177
 178 Examples of how the data generated above can be analysed
 179 can be seen in the jupyter notebooks from the [qpms_ipynotebooks][]
 180 repository in the `AaroBEC` directory.
 181
 182
 183
 184
 185
 186
 187 [qpms_ipynotebooks]: https://version.aalto.fi/gitlab/qd/qpms_ipynotebooks
 188 [scatsystem.h]: @ref scatsystem.h
 189 [qpms_scatsys_t]: @ref qpms_scatsys_t
 190 [scuff-tmatrix]: https://homerreid.github.io/scuff-em-documentation/applications/scuff-tmatrix/scuff-tmatrix/