Commit 9b2d1586 authored by Jan Oliver Oelerich's avatar Jan Oliver Oelerich

Intermediate state of plasmon scattering

parent 0e6ad42d
......@@ -66,6 +66,16 @@ cbed: {
average_defoci = true; # average the defoci in the output file
}
plasmon_scattering: {
enabled = false; # enable plasmon scattering
simple_mode = true; # No energy resolution, only E = 0 and E > 0
max_energy = 10; # max energy of the energy grid considered for plasmon energy transfer in eV
energy_grid_density = 10; # density of the energy grid in 1/eV
mean_free_path = 120; # mean free path of a plasmon in nm
plasmon_energy = 16.7; # plasmon energy in eV
plasmon_fwhm =3.7; # plasmon energy FWHM in eV
}
http_reporting: { # send POST status requests of the simulation to some HTTP endpoint.
enabled = false; # whether HTTP reporting is enabled.
url = ""; # The URL to POST to.
......
......@@ -25,6 +25,7 @@
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
......
......@@ -27,6 +27,7 @@
#include "../utilities/output.hpp"
#include "../libatomic/Scattering.hpp"
#include "Slice.hpp"
#include "../libatomic/Plasmons.hpp"
using namespace std;
using namespace stemsalabim;
......@@ -82,6 +83,7 @@ void stemsalabim::Crystal::init(const std::string &crystal_file_content) {
_size_z = std::stod(tokens[2]) * 10;
}
double volume = _size_x * _size_y * _size_z;
atomic::Scattering &sf = atomic::Scattering::getInstance();
......@@ -134,10 +136,15 @@ void stemsalabim::Crystal::init(const std::string &crystal_file_content) {
atom_id)));
}
_elements.insert(p.elementBySymbol(tokens[0]));
sf.initPotential(p.elementBySymbol(tokens[0]));
_atom_density += 1/volume;
_electron_density += p.elementBySymbol(tokens[0])->valence() / volume;
atom_id++;
}
......
......@@ -100,6 +100,14 @@ namespace stemsalabim {
return _atoms.size();
}
double atomicDensity() const {
return _atom_density;
}
double electronDensity() const {
return _electron_density;
}
/*!
* Get the number of different elements that appear in the crystal
* @return number of different elements
......@@ -176,6 +184,9 @@ namespace stemsalabim {
std::vector<std::tuple<double, double, double, double>> _fields;
double _size_x, _size_y, _size_z;
double _electron_density{0};
double _atom_density{0};
/*!
* Generates the slices
*/
......
......@@ -139,11 +139,11 @@ void GridManager::generateGrids() {
double k2 = sqrt(pow(kx(ix), 2) + pow(ky(iy), 2)) * p.wavelength();
int index = algorithms::getIndexOfAdaptiveSpace(k2 * 1e3,
get<0>(p.adfDetectorAngles()),
get<1>(p.adfDetectorAngles()),
get<2>(p.adfDetectorAngles()),
p.adfDetectorIntervalExponent(),
true);
get<0>(p.adfDetectorAngles()),
get<1>(p.adfDetectorAngles()),
get<2>(p.adfDetectorAngles()),
p.adfDetectorIntervalExponent(),
true);
if(index >= 0 && index < (long) _detector_grid.size())
_adf_bin_index[ix][iy] = index;
else
......@@ -176,12 +176,27 @@ void GridManager::generateGrids() {
d /= running_sum;
}
// energy loss grid
if(p.plasmonsEnabled() && p.plasmonsSimple()) {
_energy_loss_grid.resize(2,0);
_energy_loss_grid[1] = p.plasmonsMaxEnergy();
} else if(p.plasmonsEnabled()) {
auto num_plasmon_energy_bins = (unsigned int) ceil(p.plasmonsMaxEnergy() * p.plasmonsEnergyGridDensity()) + 1;
for(unsigned int ibin = 0; ibin < num_plasmon_energy_bins; ++ibin) {
_energy_loss_grid[ibin] = ibin / (double) num_plasmon_energy_bins * p.plasmonsMaxEnergy();
}
} else {
_energy_loss_grid.resize(1,0);
}
if(p.bandwidthLimiting()) {
_bandwidth_limit_mask.init(samplingX(), samplingY());
_bandwidth_limit_mask.setIsKSpace(true);
for(unsigned int ix = 0; ix < samplingX(); ix++) {
for(unsigned int iy = 0; iy < samplingY(); iy++) {
if(pow(kx(ix), 2) + pow(ky(iy), 2) > pow(1.0 / (3.0 * max(1./scaleX(), 1./scaleY())), 2)) {
if(pow(kx(ix), 2) + pow(ky(iy), 2) > pow(1.0 / (3.0 * max(1. / scaleX(), 1. / scaleY())), 2)) {
_bandwidth_limit_mask(ix, iy) = 0.0;
} else {
_bandwidth_limit_mask(ix, iy) = 1.0;
......
......@@ -481,6 +481,15 @@ namespace stemsalabim {
return _defocus_weights;
}
/*!
* Get the energy loss grid. Each entry of the vector corresponds to the lower bound of the
* energy loss bins, i.e., 0-th entry is energy loss between 0 and the energy of index 1.
* @return vector with energy loss bins
*/
const std::vector<double> &energyLossGrid() const {
return _energy_loss_grid;
}
/*!
* Get the size of the stored CBED, taking into account bandwidth limiting
* and rescaling in X direction.
......@@ -539,6 +548,8 @@ namespace stemsalabim {
std::vector<double> _cbed_x_grid;
std::vector<double> _cbed_y_grid;
std::vector<double> _energy_loss_grid;
Wave _bandwidth_limit_mask;
};
......
This diff is collapsed.
......@@ -296,6 +296,19 @@ void Params::readParamsFromString(const string &prms) {
}
}
_cfg.lookupValue("plasmon_scattering.enabled", _plasmons_enabled);
_cfg.lookupValue("plasmon_scattering.simple_mode", _plasmons_simple);
_cfg.lookupValue("plasmon_scattering.max_energy", _plasmons_max_energy);
_cfg.lookupValue("plasmon_scattering.energy_grid_density", _plasmons_energy_grid_density);
_cfg.lookupValue("plasmon_scattering.mean_free_path", _plasmons_mean_free_path);
_cfg.lookupValue("plasmon_scattering.plasmon_energy", _plasmons_energy);
_cfg.lookupValue("plasmon_scattering.plasmon_fwhm", _plasmons_fwhm);
if(_cfg.exists("plasmon_scattering.mean_free_path"))
_plasmons_mean_free_path *= 10;
if(_tmp_dir.empty()) {
// use the path of the output file directory for tmp
int found = _output_file.find_last_of("/\\");
......
......@@ -572,6 +572,62 @@ namespace stemsalabim {
return _cbed_save_every_n_slices;
}
/*!
* Return true if plasmon scattering is enabled.
* @return true if enabled
*/
bool plasmonsEnabled() const {
return _plasmons_enabled;
}
/*!
* Return true when the simple plasmon mode is activated
* @return true if enabled
*/
bool plasmonsSimple() const {
return _plasmons_simple;
}
/*!
* Return the maximal energy for the plasmon energy grid
* @return plasmon maximal energy
*/
float plasmonsMaxEnergy() const {
return _plasmons_max_energy;
}
/*!
* Return the density in 1/eV for the plasmon energy grid.
* @return density in 1/eV
*/
float plasmonsEnergyGridDensity() const {
return _plasmons_energy_grid_density;
}
/*!
* Plasmon mean free path in nm
* @return mean free path in nm
*/
float plasmonMeanFreePath() const {
return _plasmons_mean_free_path;
}
/*!
* Plasmon energy in eV
* @return plasmon energy in eV
*/
float plasmonEnergy() const {
return _plasmons_energy;
}
/*!
* Plasmon energy full with half maximum in eV
* @return plasmon energy FWHM in eV
*/
float plasmonFWHM() const {
return _plasmons_fwhm;
}
/// electron rest mass in keV. \f$m_0 c^2\f$ in keV
constexpr static double ELECTRON_REST_MASS = 510.99906;
......@@ -651,6 +707,15 @@ namespace stemsalabim {
bool _cbed_average_configurations{true};
bool _cbed_average_defoci{true};
// plasmon scattering
bool _plasmons_enabled{false};
bool _plasmons_simple{false};
float _plasmons_max_energy{10};
float _plasmons_energy_grid_density{10};
float _plasmons_mean_free_path{110};
float _plasmons_energy{16.6};
float _plasmons_fwhm{3.7};
bool _http_reporting{false};
std::string _http_url{""};
std::string _http_auth_user{""};
......
This diff is collapsed.
......@@ -80,6 +80,14 @@ namespace stemsalabim {
wave *= _transmission_function;
}
/*!
* Transmit a Complex2D wave through the transmission function.
* @param wave The electronic wave to transmit
*/
const Wave & transmissionFunction() const {
return _transmission_function;
}
/*!
* Add an atom to this slice.
* @param atom
......
......@@ -103,6 +103,14 @@ namespace stemsalabim { namespace atomic {
return _atomic_weight;
}
/*!
* Get the atomic valence of the element.
* @return atomic valence
*/
int valence() const {
return _valence;
}
/*!
* Set atomic weight of the element.
* @param atomic_weight atomic weight
......@@ -127,6 +135,14 @@ namespace stemsalabim { namespace atomic {
_atomic_radius = atomic_radius;
}
/*!
* Set atomic valence of the element.
* @param valence atomic valence
*/
void setValence(const int valence) {
_valence = valence;
}
/*!
* Set the Doyle Turner scattering parameters of the element.
* They are encoded as a sequence in a vector of doubles, where
......@@ -187,6 +203,7 @@ namespace stemsalabim { namespace atomic {
unsigned _atomic_number;
float _atomic_weight;
float _atomic_radius;
int _valence;
float _lattice_constant;
std::vector<double> _scattering_params_dt;
std::string _lattice_structure;
......@@ -248,6 +265,7 @@ namespace stemsalabim { namespace atomic {
e->setLatticeStructure(val.get("lattice_structure", "").asString());
e->setAtomicWeight(val.get("atomic_weight", 0.f).asFloat());
e->setValence(val.get("valence_electrons", 0).asInt());
try {
// convert pm to angstrom
......
/*
* 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/>.
*
*/
#ifndef STEMSALABIM_PLASMONS_HPP
#define STEMSALABIM_PLASMONS_HPP
#include <memory>
#include <mutex>
#include <algorithm>
#include <cmath>
#include "../classes/GridManager.hpp"
#include "../classes/IO.hpp"
namespace stemsalabim { namespace atomic {
class Plasmons {
public:
static Plasmons &
getInstance() {
static Plasmons instance;
return instance;
}
void populateCache(const std::shared_ptr<GridManager> & gridman, double electron_density, double atom_density) {
if(_cache_populated)
return;
Params &p = Params::getInstance();
if(!p.plasmonsEnabled())
return;
unsigned int lx = gridman->samplingX();
unsigned int ly = gridman->samplingY();
double sdense2 = pow(1. / 3. * p.samplingDensity(), 2);
double energy_delta = gridman->energyLossGrid()[1];
int fine_grid_num = 10;
if(p.plasmonsSimple())
fine_grid_num = 1000;
_cache.resize(gridman->energyLossGrid().size()-1);
//unsigned int energy_grid_index = 0;
for(unsigned int eii = 0; eii < _cache.size(); ++eii) {
_cache[eii].init(lx, ly);
_cache[eii].setIsKSpace(true);
for(unsigned int ix = 0; ix < lx; ix++) {
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) {
// integrate over energy
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](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();
_cache[eii].backwardFFT();
IO::dumpWave(output::fmt("w%d_r.dat", eii), _cache[eii]);
exit(0);
}
_cache_populated = true;
}
const Wave & scatteringFunction(unsigned int energy_loss_bin_index) const {
return _cache[energy_loss_bin_index];
}
double beta0() const {
return _beta_0;
}
double beta1() const {
return _beta_1;
}
private:
explicit Plasmons() {
Params &p = Params::getInstance();
_plasmon_energy = p.plasmonEnergy();
_plasmon_fwhm = p.plasmonFWHM();
_incident_energy = p.beamEnergy();
_beta_0 = exp(-p.sliceThickness() / p.plasmonMeanFreePath());
_beta_1 = p.sliceThickness() / p.plasmonMeanFreePath() * _beta_0;
}
/*!
* Calculates the double differential cross section for plasmon scattering with an incident
* electron wave by using the free electron gas model and the Lindhardt approximation.
* @param E energy in eV
* @param theta scattering angle in rad
* @return double differential cross section in angstroms^2/eV
*/
double plasmonDoubleDifferentialCrossSection(double E, double theta, double atom_density, double electron_density) {
// calculate fermi energy in eV
// 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 fermi_factor = numerical_factor_2 * numerical_factor;
// fermi energy of free electron gas in eV
double fermi_energy = pow(electron_density, 2. / 3.) * fermi_factor;
// characteristic angle of plasmons in rad
double characteristic_angle = (_incident_energy + ELECTRON_REST_MASS) /
(2 * ELECTRON_REST_MASS + _incident_energy) *
E /
_incident_energy *
1e-3;
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;
// double differential cross section of plasmon scattering in angstroms^2/eV
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.));
return B;
}
double _incident_energy;
double _plasmon_energy;
double _plasmon_fwhm;
double _beta_0;
double _beta_1;
bool _cache_populated{false};
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
constexpr static double ELECTRON_REST_MASS = 510.99906;
/// bohr radius
constexpr static double a0 = 5.2917721067e-11;
};
}}
#endif //STEMSALABIM_PLASMONS_HPP
This diff is collapsed.
......@@ -315,6 +315,29 @@ namespace stemsalabim {
return *this;
}
/*!
* Multiply with a real number, element wise.
* @param m the number
* @return this
*/
Wave &operator*=(const double xf) {
for(unsigned int i = 0; i < _lx * _ly; i++)
_data[i] *= xf;
return *this;
}
/*!
* Assign a real number to all elements of this function. The imaginary
* parts will be zero.
* @param m the number
* @return this
*/
Wave &operator=(const double xf) {
for(unsigned int i = 0; i < _lx * _ly; i++)
_data[i] = xf;
return *this;
}
/*!
* Initialize the class, i.e., allocate memory.
* @param lx width
......
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