Commit 06529fb5 authored by Jan Oliver Oelerich's avatar Jan Oliver Oelerich

Implementation of the plasmon stuff. Compiles and runs, but who knows if it is meaningful physics.

parent 831d08e5
......@@ -26,7 +26,6 @@
set(INTERNAL_FILES
libatomic/Scattering.hpp
libatomic/Plasmons.hpp
classes/Crystal.cpp classes/Crystal.hpp
classes/Params.cpp classes/Params.hpp
classes/Simulation.cpp classes/Simulation.hpp
classes/IO.cpp classes/IO.hpp
......
......@@ -149,9 +149,11 @@ void IO::initNCFile(const shared_ptr<GridManager> &gridman, const shared_ptr<ato
{
auto paramsGroup = f.defineGroup("params");
auto defocusDim = paramsGroup.defineDim("defocus", p.numberOfDefoci());
auto plasmonDim = paramsGroup.defineDim("plasmon_energies", gridman->energyLossGrid().size());
paramsGroup.defineVar<float>("defocus", vd({defocusDim}));
paramsGroup.defineVar<float>("defocus_weights", vd({defocusDim}));
paramsGroup.defineVar<float>("plasmon_energies", vd({plasmonDim}));
}
if(p.adf()) {
......@@ -169,15 +171,18 @@ void IO::initNCFile(const shared_ptr<GridManager> &gridman, const shared_ptr<ato
auto fpDim = g.defineDim("adf_phonon", p.adfAverageConfigurations() ? 1 : p.numberOfConfigurations());
auto sliceDim = g.defineDim("adf_slice", gridman->adfSliceCoords().size());
auto coordinateDim = g.defineDim("coordinate_dim", 2);
auto plasmonDim = g.defineDim("adf_plasmon_energies", gridman->energyLossGrid().size());
auto v = g.defineVar<float>("adf_intensities", vd({defocusDim, xDim, yDim, fpDim, sliceDim, angleDim}));
auto v = g.defineVar<float>("adf_intensities",
vd({defocusDim, xDim, yDim, fpDim, sliceDim, plasmonDim, angleDim}));
g.defineVar<float>("center_of_mass", vd({defocusDim, xDim, yDim, fpDim, sliceDim, coordinateDim}));
g.defineVar<double>("adf_probe_x_grid", vd({xDim}));
g.defineVar<double>("adf_probe_y_grid", vd({yDim}));
g.defineVar<double>("adf_detector_grid", vd({angleDim}));
g.defineVar<double>("adf_slice_coords", vd({sliceDim}));
v.chunking(vs({1, adf_points_x, adf_points_y, 1, gridman->adfSliceCoords().size(), p.adfNumberOfDetectorAngles()}));
v.chunking(vs({1, adf_points_x, adf_points_y, 1, gridman->adfSliceCoords().size(),
gridman->energyLossGrid().size(), p.adfNumberOfDetectorAngles()}));
if(p.outputCompress())
v.deflate(1);
......@@ -202,15 +207,18 @@ void IO::initNCFile(const shared_ptr<GridManager> &gridman, const shared_ptr<ato
auto defDim = g.defineDim("cbed_defocus", p.cbedAverageDefoci() ? 1 : p.numberOfDefoci());
auto fpDim = g.defineDim("cbed_phonon", p.cbedAverageConfigurations() ? 1 : p.numberOfConfigurations());
auto sliceDim = g.defineDim("cbed_slice", gridman->cbedSliceCoords().size());
auto plasmonDim = g.defineDim("cbed_plasmon_energies", gridman->energyLossGrid().size());
auto v = g.defineVar<float>("cbed_intensities", vd({defDim, xDim, yDim, fpDim, sliceDim, kxDim, kyDim}));
auto v = g.defineVar<float>("cbed_intensities",
vd({defDim, xDim, yDim, fpDim, sliceDim, plasmonDim, kxDim, kyDim}));
g.defineVar<double>("cbed_probe_x_grid", vd({xDim}));
g.defineVar<double>("cbed_probe_y_grid", vd({yDim}));
g.defineVar<double>("cbed_x_grid", vd({kxDim}));
g.defineVar<double>("cbed_y_grid", vd({kyDim}));
g.defineVar<double>("cbed_slice_coords", vd({sliceDim}));
v.chunking(vs({1, 1, 1, 1, gridman->cbedSliceCoords().size(), k_sampling_x, k_sampling_y}));
v.chunking(vs({1, 1, 1, 1, gridman->cbedSliceCoords().size(), gridman->energyLossGrid().size(), k_sampling_x,
k_sampling_y}));
if(p.outputCompress())
v.deflate(1);
......@@ -228,7 +236,9 @@ void IO::initBuffers(const shared_ptr<GridManager> &gridman) {
Params &p = Params::getInstance();
if(p.adf()) {
_adf_intensities_buffer.resize(gridman->adfDetectorGrid().size() * gridman->adfSliceCoords().size(), 0);
_adf_intensities_buffer.resize(gridman->adfDetectorGrid().size() *
gridman->adfSliceCoords().size() *
gridman->energyLossGrid().size(), 0);
_com_buffer.resize(2 * gridman->adfSliceCoords().size(), 0);
}
......@@ -236,7 +246,10 @@ void IO::initBuffers(const shared_ptr<GridManager> &gridman) {
unsigned int k_sampling_x = gridman->storedCbedSizeX();
unsigned int k_sampling_y = gridman->storedCbedSizeY();
// resize the write buffer
_cbed_intensities_buffer.resize(k_sampling_x * k_sampling_y * gridman->cbedSliceCoords().size(), 0);
_cbed_intensities_buffer.resize(k_sampling_x *
k_sampling_y *
gridman->cbedSliceCoords().size() *
gridman->energyLossGrid().size(), 0);
}
}
......@@ -351,6 +364,17 @@ void IO::writeParams(const shared_ptr<GridManager> &gridman, const string &conf_
g.att("enabled", p.fpEnabled());
}
{
auto g = pg.defineGroup("plasmon_scattering");
g.att("enabled", p.plasmonsEnabled());
g.att("simple_mode", p.plasmonsSimple());
g.att("max_energy", p.plasmonsMaxEnergy());
g.att("energy_grid_density", p.plasmonsEnergyGridDensity());
g.att("mean_free_path", p.plasmonMeanFreePath());
g.att("plasmon_energy", p.plasmonEnergy());
g.att("plasmon_fwhm", p.plasmonFWHM());
}
} catch(NcException &e) {
output::error("%s\n", e.what());
}
......@@ -548,6 +572,7 @@ void IO::writeAdfIntensities(unsigned int idefocus, unsigned int iconf, const sh
try {
unsigned long nangles = gridman->adfDetectorGrid().size();
unsigned int n_slices_total = (unsigned int) gridman->adfSliceCoords().size();
unsigned long nplasmons = gridman->energyLossGrid().size();
unsigned int idefocus_store = p.adfAverageDefoci() ? 0 : idefocus;
unsigned int iconf_store = p.adfAverageConfigurations() ? 0 : iconf;
......@@ -562,8 +587,8 @@ void IO::writeAdfIntensities(unsigned int idefocus, unsigned int iconf, const sh
if(iconf > iconf_store || idefocus > idefocus_store) {
v.get(vs({idefocus_store, point.adf_row, point.adf_col, iconf_store, 0, 0}),
vs({1, 1, 1, 1, n_slices_total, nangles}),
v.get(vs({idefocus_store, point.adf_row, point.adf_col, iconf_store, 0, 0, 0}),
vs({1, 1, 1, 1, n_slices_total, nplasmons, nangles}),
_adf_intensities_buffer);
for(size_t j = 0; j < _adf_intensities_buffer.size(); ++j)
......@@ -574,8 +599,8 @@ void IO::writeAdfIntensities(unsigned int idefocus, unsigned int iconf, const sh
_adf_intensities_buffer[j] = w * ibuf->value(point.adf_intensities, j);
}
v.put(vs({idefocus_store, point.adf_row, point.adf_col, iconf_store, 0, 0}),
vs({1, 1, 1, 1, n_slices_total, nangles}),
v.put(vs({idefocus_store, point.adf_row, point.adf_col, iconf_store, 0, 0, 0}),
vs({1, 1, 1, 1, n_slices_total, nplasmons, nangles}),
_adf_intensities_buffer);
......@@ -618,6 +643,7 @@ void IO::writeCBEDIntensities(unsigned int idefocus, unsigned int iconf, const s
unsigned int idefocus_store = p.cbedAverageDefoci() ? 0 : idefocus;
unsigned int iconf_store = p.cbedAverageConfigurations() ? 0 : iconf;
unsigned long nplasmons = gridman->energyLossGrid().size();
float w = 1.0;
......@@ -632,8 +658,8 @@ void IO::writeCBEDIntensities(unsigned int idefocus, unsigned int iconf, const s
if(iconf > iconf_store || idefocus > idefocus_store) {
v.get(vs({idefocus_store, point.cbed_row, point.cbed_col, iconf_store, 0, 0, 0}),
vs({1, 1, 1, 1, n_slices_total, nx, ny}),
v.get(vs({idefocus_store, point.cbed_row, point.cbed_col, iconf_store, 0, 0, 0, 0}),
vs({1, 1, 1, 1, n_slices_total, nplasmons, nx, ny}),
_cbed_intensities_buffer);
for(size_t j = 0; j < _cbed_intensities_buffer.size(); ++j)
......@@ -646,8 +672,8 @@ void IO::writeCBEDIntensities(unsigned int idefocus, unsigned int iconf, const s
}
v.put(vs({idefocus_store, point.cbed_row, point.cbed_col, iconf_store, 0, 0, 0}),
vs({1, 1, 1, 1, n_slices_total, nx, ny}),
v.put(vs({idefocus_store, point.cbed_row, point.cbed_col, iconf_store, 0, 0, 0, 0}),
vs({1, 1, 1, 1, n_slices_total, nplasmons, nx, ny}),
_cbed_intensities_buffer);
} catch(NcException &e) {
......@@ -843,6 +869,7 @@ void IO::writeGrids(const shared_ptr<GridManager> &gridman) {
{
auto g = f.group("params");
g.var("defocus").put(vs({0}), vs({gridman->defoci().size()}), gridman->defoci());
g.var("plasmon_energies").put(vs({0}), vs({gridman->energyLossGrid().size()}), gridman->energyLossGrid());
g.var("defocus_weights").put(vs({0}), vs({gridman->defocusWeights().size()}), gridman->defocusWeights());
}
......@@ -1030,6 +1057,7 @@ std::shared_ptr<stemsalabim::atomic::Cell> IO::initCrystalFromNCFile(const std::
cell->addElement(elp.elementBySymbol(at));
}
cell->setLengths(system_lengths[0], system_lengths[1], system_lengths[2]);
for(size_t i_atom = 0; i_atom < e_len[1]; i_atom++) {
cell->addAtom(make_shared<atomic::Atom>(lattice_coords[3 * i_atom + 0],
lattice_coords[3 * i_atom + 1],
......@@ -1039,7 +1067,6 @@ std::shared_ptr<stemsalabim::atomic::Cell> IO::initCrystalFromNCFile(const std::
i_atom));
}
cell->setLengths(system_lengths[0], system_lengths[1], system_lengths[2]);
return cell;
}
......
......@@ -254,6 +254,15 @@ void Params::broadcast(int rank) {
mpi_env.broadcast(_adf_average_defoci, rank);
mpi_env.broadcast(_adf_save_every_n_slices, rank);
mpi_env.broadcast(_plasmons_enabled, rank);
mpi_env.broadcast(_plasmons_simple, rank);
mpi_env.broadcast(_plasmons_max_energy, rank);
mpi_env.broadcast(_plasmons_energy_grid_density, rank);
mpi_env.broadcast(_plasmons_mean_free_path, rank);
mpi_env.broadcast(_plasmons_energy, rank);
mpi_env.broadcast(_plasmons_fwhm, rank);
}
......@@ -381,5 +390,18 @@ void Params::readParamsFromNCFile(const std::string & path) {
}
{
auto gg = g.group("plasmon_scattering");
_plasmons_enabled = gg.is("enabled");
_plasmons_simple = gg.is("simple_mode");
_plasmons_max_energy = gg.att<float>("max_energy");
_plasmons_energy_grid_density = gg.att<float>("energy_grid_density");
_plasmons_mean_free_path = gg.att<float>("mean_free_path");
_plasmons_energy = gg.att<float>("plasmon_energy");
_plasmons_fwhm = gg.att<float>("plasmon_fwhm");
}
}
......@@ -22,7 +22,6 @@
#include "Simulation.hpp"
#include <thread>
#include "../libatomic/Scattering.hpp"
#include "../utilities/http.hpp"
#include "../utilities/timings.hpp"
#include "../libatomic/Plasmons.hpp"
......@@ -92,7 +91,7 @@ void Simulation::init() {
// plasmon stuff. Needs to come after crystal and Gridmanager initialization!
if(p.plasmonsEnabled()) {
atomic::Plasmons &plasmons = atomic::Plasmons::getInstance();
plasmons.populateCache(_gridman, _crystal->electronDensity(), _crystal->atomicDensity());
plasmons.populateCache(_gridman, p.cell()->electronDensity(), p.cell()->atomicDensity());
IO::dumpWave("scattering0.dat", plasmons.scatteringFunction(0));
}
}
......@@ -305,6 +304,7 @@ void Simulation::calculatePixel(ScanPoint &point, const double defocus) {
std::vector<Wave> waves(_gridman->energyLossGrid().size());
Wave tmpwave;
if(p.plasmonsEnabled()) {
tmpwave.init(_gridman->samplingX(), _gridman->samplingY());
}
......@@ -323,12 +323,6 @@ void Simulation::calculatePixel(ScanPoint &point, const double defocus) {
// scattering
for(unsigned int energy_i = 1; energy_i < _gridman->energyLossGrid().size(); energy_i++) {
// tmpwave = waves[0];
// slice->transmit(tmpwave);
// tmpwave *= plasmons.scatteringFunction(energy_i-1);
// tmpwave *= sqrt(plasmons.beta1());
//
// waves[energy_i] += tmpwave;
// plural scattering includes single scattering with energy_j = 0
for(unsigned int energy_j = 0; energy_j < energy_i; energy_j++) {
......@@ -360,8 +354,6 @@ void Simulation::calculatePixel(ScanPoint &point, const double defocus) {
if(point.adf && _gridman->adfStoreSlice(slice->id())) {
std::fill(detector_intensities.begin(), detector_intensities.end(), 0);
// here, we iterate the k space and calculate the wave angles.
// Then, we determine the corresponding detector angle and sum the intensity up.
for(unsigned int ip = 0; ip < _gridman->energyLossGrid().size(); ip++) {
......@@ -372,34 +364,33 @@ void Simulation::calculatePixel(ScanPoint &point, const double defocus) {
int ind = _gridman->adfBinIndex(ix, iy);
if(ind >= 0) {
// check for not a number
auto intens = (float) pow(abs(waves[ip](ix, iy)), 2);
if(!::isnan(intens))
detector_intensities[ind] += intens;
if(ind >= 0) {
// check for not a number
auto intens = (float) pow(abs(waves[ip](ix, iy)), 2);
if(!::isnan(intens))
detector_intensities[ind] += intens;
}
}
}
adf_intensities.insert(adf_intensities.end(), detector_intensities.begin(), detector_intensities.end());
}
// calculate center of mass
double total_intensity = wave.integratedIntensity(), cx = 0, cy = 0;
for(unsigned int ix = 0; ix < _gridman->samplingX(); ix++) {
for(unsigned int iy = 0; iy < _gridman->samplingY(); iy++) {
auto intens = (float) pow(abs(wave(ix, iy)), 2);
cx += intens * _gridman->kx(ix);
cy += intens * _gridman->ky(iy);
}
// calculate center of mass
double total_intensity = waves[0].integratedIntensity(), cx = 0, cy = 0;
for(unsigned int ix = 0; ix < _gridman->samplingX(); ix++) {
for(unsigned int iy = 0; iy < _gridman->samplingY(); iy++) {
auto intens = (float) pow(abs(waves[0](ix, iy)), 2);
cx += intens * _gridman->kx(ix);
cy += intens * _gridman->ky(iy);
}
cx /= total_intensity;
cy /= total_intensity;
center_of_mass.push_back((float) (cx * p.wavelength() * 1e3));
center_of_mass.push_back((float) (cy * p.wavelength() * 1e3));
}
cx /= total_intensity;
cy /= total_intensity;
}
center_of_mass.push_back((float) (cx * p.wavelength() * 1e3));
center_of_mass.push_back((float) (cy * p.wavelength() * 1e3));
}
......@@ -460,8 +451,9 @@ void Simulation::calculatePixel(ScanPoint &point, const double defocus) {
}
}
for(auto &c: waves)
for(auto &c: waves) {
c.backwardFFT();
}
}
if(point.adf) {
......
......@@ -46,7 +46,7 @@ namespace stemsalabim { namespace atomic {
template<class Archive>
void serialize(Archive &ar) {
ar(_atoms, _elements, _size_x, _size_y, _size_z);
ar(_atoms, _elements, _size_x, _size_y, _size_z, _atom_density, _volume, _electron_density);
}
/*!
......@@ -91,6 +91,14 @@ namespace stemsalabim { namespace atomic {
return _elements.size();
}
double atomicDensity() const {
return _atom_density;
}
double electronDensity() const {
return _electron_density;
}
/*!
* Get a set of all elements contained in the crystal.
*/
......@@ -116,12 +124,15 @@ namespace stemsalabim { namespace atomic {
*/
void addAtom(const std::shared_ptr<Atom> &a) {
_atoms.push_back(a);
_atom_density += 1/_volume;
_electron_density += a->getElement()->valence() / _volume;
}
void setLengths(float lx, float ly, float lz) {
_size_x = lx;
_size_y = ly;
_size_z = lz;
_volume = lx*ly*lz;
}
private:
......@@ -129,6 +140,10 @@ namespace stemsalabim { namespace atomic {
std::set<std::shared_ptr<Element>> _elements;
double _size_x, _size_y, _size_z;
double _volume{0};
double _electron_density{0};
double _atom_density{0};
};
}}
#endif //STEMSALABIM_CELL_H
......@@ -52,7 +52,8 @@ namespace stemsalabim { namespace atomic {
unsigned int lx = gridman->samplingX();
unsigned int ly = gridman->samplingY();
double sdense2 = pow(1. / 3. * p.samplingDensity(), 2);
double sdense2_x = pow(1. / 3. * gridman->densityX(), 2);
double sdense2_y = pow(1. / 3. * gridman->densityY(), 2);
double energy_delta = gridman->energyLossGrid()[1];
int fine_grid_num = 10;
......@@ -71,26 +72,27 @@ namespace stemsalabim { namespace atomic {
for(unsigned int iy = 0; iy < ly; iy++) {
double theta = sqrt(pow(gridman->kx(ix), 2) + pow(gridman->ky(iy), 2)) * p.wavelength();
if(!p.bandwidthLimiting() || pow(gridman->kx(ix), 2) + pow(gridman->ky(iy), 2) < sdense2) {
if(!p.bandwidthLimiting() || pow(gridman->kx(ix), 2) / sdense2_x + pow(gridman->ky(iy), 2) / sdense2_y < 1) {
// integrate over energy
// @TODO: Replace with integration quad library
double integral = 0;
for(int fine_grid_i = 1; fine_grid_i <= fine_grid_num; fine_grid_i++) {
double energy = gridman->energyLossGrid()[eii] + (float)fine_grid_i / fine_grid_num * energy_delta;
integral += plasmonDoubleDifferentialCrossSection(energy, theta, atom_density, electron_density) * energy_delta / fine_grid_num;
}
_cache[eii](ix, iy) = integral;//(float)sqrt(integral);
_cache[eii](ix, iy) = (float)sqrt(integral);
}
}
}
_cache[eii](0, 0) = _cache[eii](0, 1);
IO::dumpWave(output::fmt("w%d_k.dat", eii), _cache[eii]);
printf("wl %e\n", gridman->kx(1) * p.wavelength());
printf("wl %e\n", gridman->kx(lx/2) * p.wavelength());
printf("wl %e\n", p.wavelength());
_cache[eii].normalize();
// IO::dumpWave(output::fmt("w%d_k.dat", eii), _cache[eii]);
// printf("wl %e\n", gridman->kx(1) * p.wavelength());
// printf("wl %e\n", gridman->kx(lx/2) * p.wavelength());
// printf("wl %e\n", p.wavelength());
_cache[eii].backwardFFT();
IO::dumpWave(output::fmt("w%d_r.dat", eii), _cache[eii]);
exit(0);
_cache[eii].normalize();
// IO::dumpWave(output::fmt("w%d_r.dat", eii), _cache[eii]);
// exit(0);
}
_cache_populated = true;
......@@ -133,7 +135,7 @@ namespace stemsalabim { namespace atomic {
// this factor is pow(3*pi^2, 2./3.)
constexpr double numerical_factor = 9.570780000627304;
// this factor is hbar^2 / (2 * m_e * e) * 10^20 (eV A^2)
constexpr double numerical_factor_2 = 3.809982080226884;
constexpr double numerical_factor_2 = 0.03809982080226884;
constexpr double fermi_factor = numerical_factor_2 * numerical_factor;
// fermi energy of free electron gas in eV
......@@ -149,7 +151,6 @@ namespace stemsalabim { namespace atomic {
double factor = 2 * _incident_energy * 1000 * (theta * theta + characteristic_angle * characteristic_angle);
double g = 3. / 5. * fermi_energy / _plasmon_energy;
g = 3;
// plasmon energy dispersion
double Ep = _plasmon_energy + g * factor;
......@@ -158,7 +159,7 @@ namespace stemsalabim { namespace atomic {
double B = 1 / (atom_density * a0 * 2 * pi * pi);
B = sin(theta) / (atom_density * a0 * pi);
B *= 2 / factor;
B *= E * _plasmon_fwhm * Ep * Ep / (pow(E * E - Ep * Ep, 2.) + pow(E * _plasmon_fwhm, 2.));
B *= E * _plasmon_fwhm * _plasmon_energy * _plasmon_energy / (pow(E * E - _plasmon_energy * _plasmon_energy - 2*g*factor*_plasmon_energy, 2.) + pow(E * _plasmon_fwhm, 2.));
return B;
}
......@@ -173,14 +174,8 @@ namespace stemsalabim { namespace atomic {
std::vector<Wave> _cache;
/// this is \f$2\pi^2 a_0 e\f$, with \f$e = 2R_ya_0\f$ (\f$R_y = 13.6 eV\f$)
constexpr static double PREFACTOR_ATOMPOT = 150.4120584;
/// pi
constexpr static double pi = 3.14159265358979323846;
/// hbar
constexpr static double hbar = 1.0545718e-34;
/// electron mass
constexpr static double m_e = 9.10938356e-31;
/// electron charge
constexpr static double e = 1.60217662e-19;
/// electron rest mass in keV. \f$m_0 c^2\f$ in keV
......
......@@ -20,7 +20,7 @@
# compilation
add_executable(ssb-test main.cpp mpi.cpp queue.cpp utilities.cpp fft.cpp buffer.cpp math.cpp nc.cpp)
add_executable(ssb-test main.cpp mpi.cpp queue.cpp utilities.cpp fft.cpp buffer.cpp math.cpp nc.cpp plasmons.cpp)
# linking ################################################
......
/*
* STEMsalabim: Magical STEM image simulations
*
* Copyright (c) 2016-2018 Jan Oliver Oelerich <jan.oliver.oelerich@physik.uni-marburg.de>
* Copyright (c) 2016-2018 Structure and Technology Research Laboratory, Philipps-Universität Marburg, Germany
*
* This program is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.
*
* This program is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with this program. If not, see <http://www.gnu.org/licenses/>.
*
*/
#include "catch.hpp"
#include <complex>
#include <cmath>
TEST_CASE("Plasmons", "[plasmons]") {
}
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment