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

Fixes for plasmon scattering

parent 88c481f5
......@@ -181,8 +181,14 @@ void IO::initNCFile(const shared_ptr<GridManager> &gridman, const shared_ptr<ato
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(),
gridman->energyLossGrid().size(), p.adfNumberOfDetectorAngles()}));
// In case of plasmon energy grid, we skip the chunking of pixels to prevent errors...
if(gridman->energyLossGrid().size() == 1) {
v.chunking(vs({1, adf_points_x, adf_points_y, 1, gridman->adfSliceCoords().size(),
gridman->energyLossGrid().size(), p.adfNumberOfDetectorAngles()}));
} else {
v.chunking(vs({1, 1, 1, 1, gridman->adfSliceCoords().size(), gridman->energyLossGrid().size(),
p.adfNumberOfDetectorAngles()}));
}
if(p.outputCompress())
v.deflate(1);
......@@ -217,8 +223,14 @@ void IO::initNCFile(const shared_ptr<GridManager> &gridman, const shared_ptr<ato
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(), gridman->energyLossGrid().size(), k_sampling_x,
k_sampling_y}));
// In case of plasmon energy grid, we skip the chunking of pixels to prevent errors...
if(gridman->energyLossGrid().size() == 1) {
v.chunking(vs({1, cbed_points_x, cbed_points_y, 1, gridman->cbedSliceCoords().size(),
gridman->energyLossGrid().size(), k_sampling_x, k_sampling_y}));
} else {
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);
......
......@@ -92,7 +92,6 @@ void Simulation::init() {
if(p.plasmonsEnabled()) {
atomic::Plasmons &plasmons = atomic::Plasmons::getInstance();
plasmons.populateCache(_gridman, p.cell()->electronDensity(), p.cell()->atomicDensity());
IO::dumpWave("scattering0.dat", plasmons.scatteringFunction(0));
}
}
......@@ -731,9 +730,9 @@ void Simulation::initBuffers() {
_gridman->energyLossGrid().size()) : 0;
_com_buffer = std::make_shared<buf_type>(prms.adf() ? (2 * _gridman->adfSliceCoords().size()) : 0,
2 * prms.workPackageSize());
_adf_intensity_buffer = std::make_shared<buf_type>(number_intensities_per_pixel, 2 * prms.workPackageSize());
_cbed_intensity_buffer = std::make_shared<buf_type>(number_cbed_per_pixel, 2 * prms.workPackageSize());
prms.workPackageSize());
_adf_intensity_buffer = std::make_shared<buf_type>(number_intensities_per_pixel, prms.workPackageSize());
_cbed_intensity_buffer = std::make_shared<buf_type>(number_cbed_per_pixel, prms.workPackageSize());
}
......
......@@ -101,7 +101,7 @@ namespace stemsalabim { namespace atomic {
}
}
//_cache[eii](0, 0) = _cache[eii](0, 1);
_cache[eii].normalize();
//_cache[eii].normalize();
_cache[eii].backwardFFT();
}
......@@ -165,16 +165,14 @@ namespace stemsalabim { namespace atomic {
double factor = 2 * _incident_energy * 1000 * (theta * theta + characteristic_angle * characteristic_angle);
double g = 3. / 5. * fermi_energy / _plasmon_energy;
// plasmon dispersion
double Ep = _plasmon_energy +
3. / 5. * fermi_energy / _plasmon_energy * 2 * _incident_energy * 1000 *
(theta * theta + characteristic_angle * characteristic_angle);
// 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 *
_plasmon_energy *
_plasmon_energy /
(pow(E * E - _plasmon_energy * _plasmon_energy - 2 * g * factor * _plasmon_energy, 2.) +
pow(E * _plasmon_fwhm, 2.));
double B = sin(theta) / (theta * theta + characteristic_angle * characteristic_angle);
B *= E * _plasmon_fwhm * Ep * Ep / (pow(E * E - Ep * Ep, 2.) + pow(E * _plasmon_fwhm, 2.));
return B;
}
......
......@@ -29,6 +29,7 @@
using namespace std;
using namespace stemsalabim;
namespace op = stemsalabim::output;
const char *getTmpDir (void) {
char *tmpdir;
......@@ -47,6 +48,8 @@ void Params::initFromCLI(int argc, const char **argv) {
ArgumentParser parser("STEMsalabim usage instructions.", "Please refer to the User documentation.");
HelpFlag help(parser, "help", "Display this help menu", {'h', "help"});
ValueFlag<string> params(parser, "params", "The parameter file path.", {'p', "params"});
ValueFlag<unsigned int> num_threads(parser, "num-threads", "The number of threads per MPI proc.", {"num-threads"});
ValueFlag<unsigned int> pkgsize(parser, "package-size", "Number of tasks in a MPI work package.", {"package-size"});
try {
parser.ParseCLI(argc, argv);
......@@ -72,6 +75,16 @@ void Params::initFromCLI(int argc, const char **argv) {
output::error("a parameter file MUST be specified!\n");
}
if(num_threads) {
_cli_options.push_back(num_threads.Name());
_num_threads = get(num_threads);
}
if(pkgsize) {
_cli_options.push_back(pkgsize.Name());
_work_package_size = get(pkgsize);
}
if(_param_file.find(".nc") == std::string::npos) {
std::string tmp(getTmpDir());
tmp += "/fileXXXXXX";
......@@ -85,6 +98,7 @@ void Params::initFromCLI(int argc, const char **argv) {
}
int main(int argc, const char **argv) {
auto &mpi_env = mpi::Environment::getInstance();
Params &p = Params::getInstance();
......@@ -94,7 +108,7 @@ int main(int argc, const char **argv) {
bool remove_output_file = false;
if(p.paramsFileName().find(".nc") == std::string::npos) {
initialize_input_file(argc, argv);
initialize_input_file(argc, argv, true);
p._param_file = p.outputFilename();
remove_output_file = true;
}
......@@ -110,6 +124,45 @@ int main(int argc, const char **argv) {
// -------------- checks
double kx_max = 1.0 / (3.0 * max(1. / gridman->scaleX(), 1. / gridman->scaleY())) * 1000;
size_t mem_per_proc = 0, result_mem_per_proc = 0;
size_t fs = sizeof(float);
// transmission functions
mem_per_proc += gridman->samplingX() * gridman->samplingY() * 2 * fs * gridman->slices().size();
// wave functions for each thread during pixel calculation
mem_per_proc += gridman->samplingX() * gridman->samplingY() * 2 * fs * (
gridman->energyLossGrid().size() + 2) * p.numberOfThreads();
// results
unsigned long number_intensities_per_pixel = p.adf() ? (gridman->adfDetectorGrid().size() *
gridman->adfSliceCoords().size() *
gridman->energyLossGrid().size()) : 0;
unsigned long number_cbed_per_pixel = p.cbed() ? (gridman->storedCbedSizeX() *
gridman->storedCbedSizeY() *
gridman->cbedSliceCoords().size() *
gridman->energyLossGrid().size()) : 0;
// buffers at least!
result_mem_per_proc += p.workPackageSize() * (number_intensities_per_pixel + number_cbed_per_pixel) * fs;
op::nakedprint("General information\n");
op::nakedprint("-------------------------------\n");
op::nakedprint(" Grid size: %d x %d\n", gridman->samplingX(), gridman->samplingY());
op::nakedprint(" Scan points: %d x %d\n", gridman->adfXGrid().size(), gridman->adfYGrid().size());
op::nakedprint(" Energy bins: %d\n", gridman->energyLossGrid().size());
op::nakedprint(" Max angle %.1f mrad\n", kx_max*p.wavelength());
op::nakedprint(" Angle resolution: %.3f mrad\n", gridman->kx(1)*p.wavelength()*1000);
op::nakedprint("\nMemory estimates (lower bounds!!)\n");
op::nakedprint("-----------------------------------\n");
op::nakedprint(" Working mem/proc: %.1f GB\n", mem_per_proc/(float)(1024*1024*1024));
op::nakedprint(" Results mem/proc: %.1f GB\n", result_mem_per_proc/(float)(1024*1024*1024));
op::nakedprint(" Total mem/proc: %.1f GB\n", (result_mem_per_proc+mem_per_proc)/(float)(1024*1024*1024));
// -------------- stop checks
if(remove_output_file)
......
......@@ -34,7 +34,7 @@
namespace stemsalabim {
void initialize_input_file(int argc, const char **argv, bool with_output = true) {
void initialize_input_file(int argc, const char **argv, bool silent = false) {
auto &mpi_env = mpi::Environment::getInstance();
Params &p = Params::getInstance();
......@@ -72,7 +72,8 @@ namespace stemsalabim {
if(p.storedPotentials()) {
std::vector<std::vector<std::vector<float>>> vz(gridman->slices().size());
output::print(output::fmt("Calculating potentials for defocus %d, configuration %d...\n"),
if(!silent)
output::print(output::fmt("Calculating potentials for defocus %d, configuration %d...\n"),
st.idefocus(),
st.iconf());
......
......@@ -98,6 +98,11 @@ namespace stemsalabim {
tfm::format(std::cout, fmt.c_str(), parameters...);
}
template<typename... ParamTypes>
void nakedprint(const std::string &str_format, const ParamTypes &... parameters) {
tfm::format(std::cout, str_format.c_str(), parameters...);
}
template<typename... ParamTypes>
std::string fmt(const std::string &str_format, ParamTypes... parameters) {
return tfm::format(str_format.c_str(), parameters...);
......
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