### Reflectometry: Fit honeycomb lattice

In this example, we want to demonstrate how to fit a more complex sample. For this purpose, we utilize the reflectometry data of an artificial magnetic honeycomb lattice published by A. Glavic et al., in this paper

The experiment was performed with polarized neutrons, but without polarization analysis. Since the magnetization of the sample was parallel to the neutron spin, there is no spin flip and we apply the scalar theory to this problem. This is primarily done to speed up computations: when the polarized computational engine is utilized the fitting procedure takes roughly three times as long.

#### Experimental data

The experimental data consists of four datasets that should be fitted simultaneously. These datasets arise from the two polarization channels for up and down polarization of the incoming beam and both of these channels are measured at two temperatures (300K and 150K).

All of this is measured on the same sample, so all parameters are assumed to be the same, except the magnetization being temperature dependent. Therefore, we introduce a scaling parameter for the magnetization as the ratio of the magnetizations at 150K and 300K: $M_{s150} = M_{150K} / M_{300K}$.

#### Magnetization model

To model a magnetic material, one can assign a magnetization vector to any material, as is demonstrated in the basic polarized reflectometry tutorial. When a non-vanishing magnetization vector is specified for at least one layer in a sample, BornAgain will automatically utilize the polarized computational engine. This leads to lower performance as the computations are more invovled.

In this example, the magnetization is (anti)parallel to the neutron spin and hence we instead parametrize the magnetic layers with an effective SLD that is the sum/difference of the nuclear and their magnetic SLD:

$$\rho_\pm = \rho_{\text{N}} \pm \rho_{\text{M}}$$

Here the $+$ is chosen for incoming neutrons with spin up and $-$ is chosen for spin down neutrons.

#### Computational model

We simulate this experiment by bulding a 6 layer model: As usual the top layer is the vacuum and the bottom layer is a silicon substrate. On top of the silicon substrate, we simulate a thin oxide layer, where we fit its roughness and thickness The SLDs of these three layers are taken from the literature and kept constant.

Then we model the lattice structure with a three-layer model: two layers to account for density fluctuations in $z$-direction and another oxide layer on top. This lattice structure is assumed to be magnetic and we fit all of their SLDs, magnetic SLDs, thicknesses and roughnesses. The magnetic SLD depends on the temperature of the dataset, according to the scaling described above, where the $M_{s150}$ parameter is fitted.

All layers are modeled without absorption, i.e. no imaginary part of the SLD. Ferthermore, we apply a resolution correction as described in this tutorial with a fixed value for $\Delta Q / Q = 0.018$. The experimental data is normalized to unity, but we still fit the intensity, as is necessary due to the resolution correction.

#### Running a computation

In order to run a computation, we define some functions to utilize a common simulation function for the two spin channels and both temperatures:

def run_Simulation_300_p( qzs, params ):
return run_simulation(qzs, params, sign=1)

def run_Simulation_300_m( qzs, params ):
return run_simulation(qzs, params, sign=-1)

def run_Simulation_150_p( qzs, params ):
return run_simulation(qzs, params, sign=1, ms150=True)

def run_Simulation_150_m( qzs, params ):
return run_simulation(qzs, params, sign=-1, ms150=True)

Here, the given arguments specify whether we want to simulate a spin-up beam (sign = 1) and whether the scaling of the magnetization should be applied (ms150=True). For the latter, true means that a dataset at 150K is simulated while false corresponds to 300K and the scaling parameter is set to unity.

All four reflectivity curves are then computed using:

q_300_p, r_300_p = qr( run_Simulation_300_p( qzs, paramsInitial ) )
q_300_m, r_300_m = qr( run_Simulation_300_m( qzs, paramsInitial ) )

q_150_p, r_150_p = qr( run_Simulation_150_p( qzs, paramsInitial ) )
q_150_m, r_150_m = qr( run_Simulation_150_m( qzs, paramsInitial ) )

We choose some sensible initial parameters and these yield the following simulation result

We have chosen the initial magnetization to be zero, hence there is only a single SLD curve for both spin directions.

#### Fitting

We fit this example by utilizing the differential evolution algorithm from Scipy. As a measure for the goodness of the fit, we use the relative difference:

$$\Delta = \sum_{j = 1}^4 \frac{1}{N_j} \sum_{i = 1}^N \left( \frac{d_{ji} - s_{ji}}{d_{ji} + s_{ji}} \right)^2$$

Here the sum over $i$ sums up the fitting error at every data point as usual and the sum over $j$ adds the contributions from all four datasets. This is implemented in the FitObjective::__call__ function and the FitObjective object holds all four datasets, and also performs the corresponding simulations. In the function run_fit_differential_evolution this is dispatched to the differential evolution algorithm.

The given uncercainty of the experimental data is not taken into account.

#### Fit Result

As usual, the fit can be run with the following command:

python3 Honeycomb_fit.py fit
On a four-core workstation, the fitting procedure takes roughly 45 minutes to complete and we obtain the following result:

As can be seen from the plot of the SLDs, the magnetization is indeed larger for the measurement at lower temperature, exactly as expected.

  1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207 208 209 210 211 212 213 214 215 216 217 218 219 220 221 222 223 224 225 226 227 228 229 230 231 232 233 234 235 236 237 238 239 240 241 242 243 244 245 246 247 248 249 250 251 252 253 254 255 256 257 258 259 260 261 262 263 264 265 266 267 268 269 270 271 272 273 274 275 276 277 278 279 280 281 282 283 284 285 286 287 288 289 290 291 292 293 294 295 296 297 298 299 300 301 302 303 304 305 306 307 308 309 310 311 312 313 314 315 316 317 318 319 320 321 322 323 324 325 326 327 328 329 330 331 332 333 334 335 336 337 338 339 340 341 342 343 344 345 346 347 348 349 350 351 352 353 354 355 356 357 358 359 360 361 362 363 364 365 366 367 368 369 370 371 372 373 374 375 376 377 378 379 380 381 382 383 384 385 386 387 388 389 390 391 392 393 394 395 396 397 398 399 400 401 402 403 404 405 406 407 408 409 410 411 412 413 414 415 416 417 418 419 420 421 422 423 424 425 426 427 428 429 430 431 432 433 434 435 436 437 438 439 440 441 442 443 444 445 446 447 448 449  #!/usr/bin/env python3 """ This example demonstrates how to fit a complex experimental setup using BornAgain. It is based on real data published in https://doi.org/10.1002/advs.201700856 by A. Glavic et al. In this example we utilize the scalar reflectometry engine to fit polarized data without spin-flip for performance reasons. """ import sys import numpy import matplotlib.pyplot as plt from scipy.optimize import differential_evolution import bornagain as ba from bornagain import angstrom from zipfile import ZipFile # number of points on which the computed result is plotted scan_size = 1500 # restrict the Q-range of the data used for fitting qmin = 0.08 qmax = 1.4 #################################################################### # Create Sample and Simulation # #################################################################### def get_sample(parameters, sign, ms150=1): m_Air = ba.MaterialBySLD("Air", 0, 0) m_PyOx = ba.MaterialBySLD("PyOx", (parameters["sld_PyOx_real"] + \ sign * ms150 * parameters["msld_PyOx"] )* 1e-6, parameters["sld_PyOx_imag"] * 1e-6) m_Py2 = ba.MaterialBySLD("Py2", ( parameters["sld_Py2_real"] + \ sign * ms150 * parameters["msld_Py2"] ) * 1e-6, parameters["sld_Py2_imag"] * 1e-6) m_Py1 = ba.MaterialBySLD("Py1", ( parameters["sld_Py1_real"] + \ sign * ms150 * parameters["msld_Py1"] ) * 1e-6, parameters["sld_Py1_imag"] * 1e-6) m_SiO2 = ba.MaterialBySLD("SiO2", parameters["sld_SiO2_real"]*1e-6, parameters["sld_SiO2_imag"]*1e-6) m_Si = ba.MaterialBySLD("Substrate", parameters["sld_Si_real"]*1e-6, parameters["sld_Si_imag"]*1e-6) l_Air = ba.Layer(m_Air) l_PyOx = ba.Layer(m_PyOx, parameters["t_PyOx"]*angstrom) l_Py2 = ba.Layer(m_Py2, parameters["t_Py2"]*angstrom) l_Py1 = ba.Layer(m_Py1, parameters["t_Py1"]*angstrom) l_SiO2 = ba.Layer(m_SiO2, parameters["t_SiO2"]*angstrom) l_Si = ba.Layer(m_Si) r_PyOx = ba.LayerRoughness() r_Py2 = ba.LayerRoughness() r_Py1 = ba.LayerRoughness() r_SiO2 = ba.LayerRoughness() r_Si = ba.LayerRoughness() r_PyOx.setSigma(parameters["r_PyOx"]*angstrom) r_Py2.setSigma(parameters["r_Py2"]*angstrom) r_Py1.setSigma(parameters["r_Py1"]*angstrom) r_SiO2.setSigma(parameters["r_SiO2"]*angstrom) r_Si.setSigma(parameters["r_Si"]*angstrom) sample = ba.MultiLayer() sample.addLayer(l_Air) sample.addLayerWithTopRoughness(l_PyOx, r_PyOx) sample.addLayerWithTopRoughness(l_Py2, r_Py2) sample.addLayerWithTopRoughness(l_Py1, r_Py1) sample.addLayerWithTopRoughness(l_SiO2, r_SiO2) sample.addLayerWithTopRoughness(l_Si, r_Si) sample.setRoughnessModel(ba.RoughnessModel.NEVOT_CROCE) return sample def get_simulation(q_axis, fitParams, sign, ms150=False): n_points = 25 n_sig = 3.0 parameters = dict(fitParams, **fixedParams) q_distr = ba.RangedDistributionGaussian(n_points, n_sig) dq = parameters["dq"]*q_axis scan = ba.QSpecScan(q_axis) scan.setAbsoluteQResolution(q_distr, dq) simulation = ba.SpecularSimulation() simulation.setScan(scan) simulation.beam().setIntensity(parameters["intensity"]) if ms150: sample = get_sample(parameters=parameters, sign=sign, ms150=parameters["ms150"]) else: sample = get_sample(parameters=parameters, sign=sign, ms150=1) simulation.setSample(sample) simulation.setBackground(ba.ConstantBackground(5e-7)) return simulation def run_simulation(q_axis, fitParams, *, sign, ms150=False): simulation = get_simulation(q_axis, fitParams, sign, ms150) simulation.runSimulation() return simulation def qr(result): """ Returns two arrays that hold the q-values as well as the reflectivity from a given simulation result """ q = numpy.array(result.result().axis(ba.Axes.QSPACE)) r = numpy.array(result.result().array(ba.Axes.QSPACE)) return q, r #################################################################### # Plot Handling # #################################################################### def plot(qs, rs, exps, shifts, labels, filename): """ Plot the simulated result together with the experimental data """ fig = plt.figure() ax = fig.add_subplot(111) for q, r, exp, shift, l in zip(qs, rs, exps, shifts, labels): ax.errorbar(exp[0], exp[1]/shift, yerr=exp[2]/shift, fmt='.', markersize=0.75, linewidth=0.5) ax.plot(q, r/shift, label=l) ax.set_yscale('log') plt.legend() plt.xlabel("Q [nm${}^{-1}$]") plt.ylabel("R") plt.tight_layout() plt.savefig(filename) def plot_sld_profile(fitParams, filename): plt.figure() parameters = dict(fitParams, **fixedParams) z_300_p, sld_300_p = ba.MaterialProfile(get_sample(parameters, 1)) z_300_m, sld_300_m = ba.MaterialProfile(get_sample(parameters, -1)) z_150_p, sld_150_p = ba.MaterialProfile( get_sample(parameters, 1, ms150=parameters["ms150"])) z_150_m, sld_150_m = ba.MaterialProfile( get_sample(parameters, -1, ms150=parameters["ms150"])) plt.figure() nsld = plt.plot(z_300_p, numpy.real(sld_300_p)*1e6, label=r"300K $+$") nsld = plt.plot(z_300_m, numpy.real(sld_300_m)*1e6, label=r"300K $-$") nsld = plt.plot(z_150_p, numpy.real(sld_150_p)*1e6, label=r"150K $+$") nsld = plt.plot(z_150_m, numpy.real(sld_150_m)*1e6, label=r"150K $-$") plt.xlabel(r"$z$ [A]") plt.ylabel(r"$\delta(z) \cdot 10^6$") plt.legend() plt.tight_layout() plt.savefig(filename) plt.close() #################################################################### # Data Handling # #################################################################### def normalizeData(data): """ Removes duplicate q values from the input data, normalizes it such that the maximum of the reflectivity is unity and rescales the q-axis to inverse nm """ # delete repeated data r0 = numpy.where(data[0] - numpy.roll(data[0], 1) == 0) data = numpy.delete(data, r0, 1) data[0] = data[0]/angstrom norm = numpy.max(data[1]) data[1] = data[1]/norm data[2] = data[2]/norm # sort by q axis so = numpy.argsort(data[0]) data = data[:, so] return data def get_Experimental_data(filename, qmin, qmax): if not hasattr(get_Experimental_data, "zipfile"): get_Experimental_data.zipfile = ZipFile("honeycomb_data.zip") input_Data = numpy.genfromtxt( get_Experimental_data.zipfile.open(filename), unpack=True, usecols=(0, 2, 3)) data = normalizeData(input_Data) minIndex = numpy.argmin(numpy.abs(data[0] - qmin)) maxIndex = numpy.argmin(numpy.abs(data[0] - qmax)) return data[:, minIndex:maxIndex + 1] #################################################################### # Fit Function # #################################################################### def relative_difference(sim, exp): result = (exp - sim)/(exp + sim) return numpy.sum(result*result)/len(sim) def create_Parameter_dictionary(parameterNames, *args): return {name: value for name, value in zip(parameterNames, *args)} class FitObjective: def __init__(self, q_axis, rdata, simulationFactory, parameterNames): if isinstance(q_axis, list) and isinstance(rdata, list) and \ isinstance(simulationFactory, list): self._q = q_axis self._r = rdata self._simulationFactory = simulationFactory elif not isinstance(q_axis, list) and not isinstance(rdata, list) \ and not isinstance(simulationFactory, list): self._q = [q_axis] self._r = [rdata] self._simulationFactory = [simulationFactory] else: raise Exception("Inconsistent parameters") self._parameterNames = parameterNames def __call__(self, *args): fitParameters = create_Parameter_dictionary(self._parameterNames, *args) print(f"FitParamters = {fitParameters}") result_metric = 0 for q, r, sim in zip(self._q, self._r, self._simulationFactory): sim_result = sim(q, fitParameters).result().array() result_metric += relative_difference(sim_result, r) return result_metric def run_fit_differential_evolution(q_axis, rdata, simulationFactory, startParams): bounds = [(par[1], par[2]) for n, par in startParams.items()] parameters = [par[0] for n, par in startParams.items()] parameterNames = [n for n, par in startParams.items()] print(f"Bounds = {bounds}") objective = FitObjective(q_axis, rdata, simulationFactory, parameterNames) chi2_initial = objective(parameters) result = differential_evolution(objective, bounds, maxiter=200, popsize=len(bounds)*10, mutation=(0.5, 1.5), disp=True, tol=1e-2) resultParameters = create_Parameter_dictionary(parameterNames, result.x) chi2_final = objective(resultParameters.values()) print(f"Initial chi2: {chi2_initial}") print(f"Final chi2: {chi2_final}") return resultParameters #################################################################### # Main Function # #################################################################### if __name__ == '__main__': fixedParams = { "sld_PyOx_imag": (0, 0, 0), "sld_Py2_imag": (0, 0, 0), "sld_Py1_imag": (0, 0, 0), "sld_SiO2_imag": (0, 0, 0), "sld_Si_imag": (0, 0, 0), "sld_SiO2_real": (3.47, 3, 4), "sld_Si_real": (2.0704, 2, 3), "dq": (0.018, 0, 0.1), } if len(sys.argv) > 1 and sys.argv[1] == "fit": # some sensible start parameters for fitting startParams = { "intensity": (1.04, 0, 3), "t_PyOx": (77, 60, 100), "t_Py2": (56, 40, 70), "t_Py1": (64, 50, 80), "t_SiO2": (16, 10, 30), "sld_PyOx_real": (1.915, 1.6, 2.2), "sld_Py2_real": (5, 3, 6), "sld_Py1_real": (4.62, 3, 6), "r_PyOx": (27, 5, 35), "r_Py2": (12, 5, 20), "r_Py1": (12, 5, 20), "r_SiO2": (17, 2, 25), "r_Si": (18, 2, 25), "msld_PyOx": (0.25, 0, 1), "msld_Py2": (0.63, 0, 1), "msld_Py1": (0.64, 0, 1), "ms150": (1, 0.9, 1.1), } fit = True else: # result from our own fitting startParams = { 'intensity': 0.9482344993285265, 't_PyOx': 74.97056190221168, 't_Py2': 61.75823766477464, 't_Py1': 54.058310970786316, 't_SiO2': 23.127048588278402, 'sld_PyOx_real': 2.199791584033569, 'sld_Py2_real': 4.980316982224387, 'sld_Py1_real': 4.612135848532186, 'r_PyOx': 31.323366207013787, 'r_Py2': 9.083768897940645, 'r_Py1': 5, 'r_SiO2': 14.43455709065263, 'r_Si': 14.948233893986075, 'msld_PyOx': 0.292684104601585, 'msld_Py2': 0.5979447434271339, 'msld_Py1': 0.56376339230351, 'ms150': 1.083311702077648 } startParams = {d: (v, ) for d, v in startParams.items()} fit = False fixedParams = {d: v[0] for d, v in fixedParams.items()} paramsInitial = {d: v[0] for d, v in startParams.items()} def run_Simulation_300_p(qzs, params): return run_simulation(qzs, params, sign=1) def run_Simulation_300_m(qzs, params): return run_simulation(qzs, params, sign=-1) def run_Simulation_150_p(qzs, params): return run_simulation(qzs, params, sign=1, ms150=True) def run_Simulation_150_m(qzs, params): return run_simulation(qzs, params, sign=-1, ms150=True) qzs = numpy.linspace(qmin, qmax, scan_size) q_300_p, r_300_p = qr(run_Simulation_300_p(qzs, paramsInitial)) q_300_m, r_300_m = qr(run_Simulation_300_m(qzs, paramsInitial)) q_150_p, r_150_p = qr(run_Simulation_150_p(qzs, paramsInitial)) q_150_m, r_150_m = qr(run_Simulation_150_m(qzs, paramsInitial)) data_300_p = get_Experimental_data("honeycomb_300_p.dat", qmin, qmax) data_300_m = get_Experimental_data("honeycomb_300_m.dat", qmin, qmax) data_150_p = get_Experimental_data("honeycomb_150_p.dat", qmin, qmax) data_150_m = get_Experimental_data("honeycomb_150_m.dat", qmin, qmax) plot_sld_profile(paramsInitial, f"Honeycomb_Fit_sld_profile_initial.pdf") plot([q_300_p, q_300_m, q_150_p, q_150_m], [r_300_p, r_300_m, r_150_p, r_150_m], [data_300_p, data_300_m, data_150_p, data_150_m], [1, 1, 10, 10], ["300K $+$", "300K $-$", "150K $+$", "150K $-$"], f"Honeycomb_Fit_reflectivity_initial.pdf") # fit and plot fit if fit: dataSimTuple = [[ data_300_p[0], data_300_m[0], data_150_p[0], data_150_m[0] ], [data_300_p[1], data_300_m[1], data_150_p[1], data_150_m[1]], [ run_Simulation_300_p, run_Simulation_300_m, run_Simulation_150_p, run_Simulation_150_m ]] fitResult = run_fit_differential_evolution(*dataSimTuple, startParams) print("Fit Result:") print(fitResult) q_300_p, r_300_p = qr(run_Simulation_300_p(qzs, fitResult)) q_300_m, r_300_m = qr(run_Simulation_300_m(qzs, fitResult)) q_150_p, r_150_p = qr(run_Simulation_150_p(qzs, fitResult)) q_150_m, r_150_m = qr(run_Simulation_150_m(qzs, fitResult)) plot([q_300_p, q_300_m, q_150_p, q_150_m], [r_300_p, r_300_m, r_150_p, r_150_m], [data_300_p, data_300_m, data_150_p, data_150_m], [1, 1, 10, 10], ["300K $+$", "300K $-$", "150K $+$", "150K $-$"], f"Honeycomb_Fit_reflectivity_fit.pdf") plot_sld_profile(fitResult, f"Honeycomb_Fit_sld_profile_fit.pdf") 
Honeycomb_fit.py Reference data