finite_systems.md

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