Commit 8b27365f authored by Jan Oliver Oelerich's avatar Jan Oliver Oelerich

Somewhere in the middle of a big refactor. Introduced a wrapper around netcdf...

Somewhere in the middle of a big refactor. Introduced a wrapper around netcdf library, and the tool ssb-mkin.
parent a438aec7
......@@ -18,7 +18,7 @@ set(PACKAGE_AUTHOR_EMAIL "jan.oliver.oelerich@physik.uni-marburg.de")
project(STEMsalabim CXX)
cmake_minimum_required(VERSION 3.3)
set(CMAKE_MODULE_PATH ${CMAKE_MODULE_PATH} "${CMAKE_SOURCE_DIR}/cmake")
set(CMAKE_CXX_STANDARD 11)
set(CMAKE_CXX_STANDARD 17)
string(TIMESTAMP DATE "%Y-%m-%dT%H:%M:%S")
......
......@@ -48,7 +48,7 @@ set(INTERNAL_FILES
set(EXTERNAL_FILES
3rdparty/jsoncpp/jsoncpp.cpp
3rdparty/jsoncpp/json/json.h
3rdparty/jsoncpp/json/json-forwards.h)
3rdparty/jsoncpp/json/json-forwards.h utilities/netcdf_utils.hpp)
set_source_files_properties(${INTERNAL_FILES}
PROPERTIES COMPILE_FLAGS "-Wall -Wextra -pedantic")
......@@ -57,6 +57,7 @@ add_library(stemsalabim_lib ${INTERNAL_FILES} ${EXTERNAL_FILES})
add_executable(stemsalabim stemsalabim.cpp)
add_executable(ssb-chk ssb-chk.cpp)
add_executable(ssb-mkin ssb-mkin.cpp)
#add_subdirectory(3rdparty/pybind11)
#pybind11_add_module(ssbpy py_stemsalabim.cpp)
......@@ -69,6 +70,7 @@ target_compile_features(stemsalabim PRIVATE cxx_range_for)
target_link_libraries(stemsalabim_lib m ${LIBS})
target_link_libraries(stemsalabim stemsalabim_lib)
target_link_libraries(ssb-chk stemsalabim_lib)
target_link_libraries(ssb-mkin stemsalabim_lib)
# MPI ####################################################
......
......@@ -75,3 +75,20 @@ void FPConfManager::assignAtomsToSlices() {
i++;
}
}
void FPConfManager::initFromMSDs(const std::vector<float> &msds) {
_displacements.resize(_total_configurations);
assert(msds.size() >= _total_configurations * _crystal->numberOfAtoms() * 3);
size_t cnt = 0;
for(size_t j = 0; j < _total_configurations; j++) {
_displacements[j].resize(_crystal->numberOfAtoms());
for(size_t i = 0; i < _crystal->numberOfAtoms(); i++) {
_displacements[j][i] = make_tuple(msds[cnt], msds[cnt+1], msds[cnt+2]);
cnt += 3;
}
}
}
......@@ -54,6 +54,9 @@ namespace stemsalabim {
*/
void init(const std::vector<double> &random_numbers);
void initFromMSDs(const std::vector<float> &msds);
/*!
* Set the current FP configuration index.
* @param iconf configuration index
......
This diff is collapsed.
......@@ -30,6 +30,7 @@
#include <netcdf.h>
#include "Params.hpp"
#include "../utilities/memory.hpp"
#include "../utilities/netcdf_utils.hpp"
#include "../utilities/Wave.hpp"
namespace stemsalabim {
......@@ -79,6 +80,12 @@ namespace stemsalabim {
*/
void writeCrystal(const std::shared_ptr<FPConfManager> &fpman, const std::shared_ptr<Crystal> &crystal);
/*!
* Writes the scattering potentials to the output file
* @param fpman FP configuration manager
*/
void writePotentials(const std::shared_ptr<FPConfManager> &fpman, const std::vector<std::vector<std::vector<float>>> & d);
/*!
* Write all simulation parameters to the output file.
*/
......@@ -97,7 +104,7 @@ namespace stemsalabim {
std::shared_ptr<memory::buffer::number_buffer<float>> &ibuf);
void writeAdfIntensities(unsigned int idefocus, unsigned int iconf,
const std::shared_ptr<GridManager> &gridman, const ScanPoint &point,
std::shared_ptr<memory::buffer::number_buffer<float>> &ibuf, int g_id, int v_id);
std::shared_ptr<memory::buffer::number_buffer<float>> &ibuf, NCVar & v);
/*!
* Append to the CBED intensities variable.
......@@ -112,7 +119,7 @@ namespace stemsalabim {
std::shared_ptr<memory::buffer::number_buffer<float>> &cbuf);
void writeCBEDIntensities(unsigned int idefocus, unsigned int iconf,
const std::shared_ptr<GridManager> &gridman, const ScanPoint &point,
std::shared_ptr<memory::buffer::number_buffer<float>> &cbuf, int g_id, int v_id);
std::shared_ptr<memory::buffer::number_buffer<float>> &cbuf, NCVar & v);
/*!
* Write temporary results to a binary file (custom format). This is done in MPI mode.
......@@ -150,6 +157,8 @@ namespace stemsalabim {
*/
static void dumpWave(std::string filename, const Wave &wave);
std::vector<float> readMSDs();
private:
typedef std::vector<size_t> vs;
......
......@@ -23,6 +23,9 @@
#include <thread>
#include <fstream>
#include <tuple>
#include "../utilities/netcdf_utils.hpp"
#ifdef __GNUC__
#pragma GCC diagnostic push
......@@ -68,6 +71,8 @@ void Params::initFromCLI(int argc, const char **argv) {
"diagnostics-dir",
"The directory to place diagnostics output files.",
{"diagnostics-dir"});
Flag stored_potentials(parser, "stored-potentials", "Work with stored potentials", {"stored-potentials"});
try {
parser.ParseCLI(argc, argv);
......@@ -91,6 +96,11 @@ void Params::initFromCLI(int argc, const char **argv) {
_skip_simulation = true;
}
if(stored_potentials) {
_cli_options.push_back(params.Name());
_stored_potentials = true;
}
if(params) {
_cli_options.push_back(params.Name());
_param_file = get(params);
......@@ -167,8 +177,6 @@ void Params::readParamsFromString(const string &prms) {
string prec_str = "single";
_cfg.lookupValue("application.precision", prec_str);
if(prec_str.compare("d") == 0 || prec_str.compare("double") == 0 || prec_str.compare("high") == 0)
_double_precision = true;
if(!cliParamGiven("title"))
_cfg.lookupValue("simulation.title", _title);
......@@ -310,3 +318,233 @@ void Params::readParamsFromString(const string &prms) {
_rgen.seed(_random_seed);
}
void Params::broadcast(int rank) {
auto &mpi_env = mpi::Environment::getInstance();
mpi_env.broadcast(_init_from_ncfile, rank);
mpi_env.broadcast(_random_seed, rank);
mpi_env.broadcast(_skip_simulation, rank);
mpi_env.broadcast(_title, rank);
mpi_env.broadcast(_output_file, rank);
mpi_env.broadcast(_tmp_dir, rank);
mpi_env.broadcast(_output_compress, rank);
mpi_env.broadcast(_normalize_always, rank);
mpi_env.broadcast(_bandwidth_limiting, rank);
mpi_env.broadcast(_probe_c5_coeff, rank);
mpi_env.broadcast(_probe_spherical_aberration_coeff, rank);
mpi_env.broadcast(_probe_astigmatism_ca, rank);
mpi_env.broadcast(_probe_astigmatism_angle, rank);
mpi_env.broadcast(_probe_min_aperture, rank);
mpi_env.broadcast(_probe_max_aperture, rank);
mpi_env.broadcast(_beam_energy, rank);
mpi_env.broadcast(_probe_scan_density, rank);
mpi_env.broadcast(_probe_defocus, rank);
mpi_env.broadcast(_probe_num_defoci, rank);
mpi_env.broadcast(_probe_fwhm_defoci, rank);
mpi_env.broadcast(_max_potential_radius, rank);
mpi_env.broadcast(_crystal_file, rank);
mpi_env.broadcast(_sample_density, rank);
mpi_env.broadcast(_slice_thickness, rank);
mpi_env.broadcast(_number_configurations, rank);
mpi_env.broadcast(_fixed_slicing, rank);
// CBED stuff
vector<double> cbed_scan_x(2);
vector<double> cbed_scan_y(2);
vector<unsigned int> cbed_size(2);
mpi_env.broadcast(_cbed_enabled, rank);
mpi_env.broadcast(cbed_scan_x, rank);
mpi_env.broadcast(cbed_scan_y, rank);
mpi_env.broadcast(cbed_size, rank);
mpi_env.broadcast(_cbed_save_every_n_slices, rank);
mpi_env.broadcast(_cbed_average_configurations, rank);
mpi_env.broadcast(_cbed_average_defoci, rank);
if(mpi_env.isSlave()) {
_cbed_scan_x = make_tuple(cbed_scan_x[0], cbed_scan_x[1]);
_cbed_scan_y = make_tuple(cbed_scan_y[0], cbed_scan_y[1]);
_cbed_size = make_tuple(cbed_size[0], cbed_size[1]);
}
// ADF stuff
vector<double> adf_scan_x(2);
vector<double> adf_scan_y(2);
mpi_env.broadcast(_adf_enabled, rank);
mpi_env.broadcast(adf_scan_x, rank);
mpi_env.broadcast(adf_scan_y, rank);
if(mpi_env.isSlave()) {
_adf_scan_x = make_tuple(adf_scan_x[0], adf_scan_x[1]);
_adf_scan_y = make_tuple(adf_scan_y[0], adf_scan_y[1]);
}
mpi_env.broadcast(get<0>(_adf_detector_angles), rank);
mpi_env.broadcast(get<1>(_adf_detector_angles), rank);
mpi_env.broadcast(get<2>(_adf_detector_angles), rank);
mpi_env.broadcast(_adf_detector_interval_exponent, rank);
mpi_env.broadcast(_adf_average_configurations, rank);
mpi_env.broadcast(_adf_average_defoci, rank);
mpi_env.broadcast(_adf_save_every_n_slices, rank);
mpi_env.broadcast(_http_reporting, rank);
mpi_env.broadcast(_http_url, rank);
mpi_env.broadcast(_http_auth_user, rank);
mpi_env.broadcast(_http_auth_pass, rank);
}
void Params::readParamsFromNCFile() {
NCFile f(_param_file);
auto g = f.group("params");
_init_from_ncfile = true;
{
auto gg = g.group("application");
_random_seed = gg.att<unsigned int>("random_seed");
_skip_simulation = gg.is("skip_simulation");
}
{
auto gg = g.group("simulation");
if(!cliParamGiven("title"))
_title = gg.att<string>("title");
if(!cliParamGiven("output-file"))
_output_file = gg.att<string>("output_file");
if(!cliParamGiven("tmp-dir"))
_tmp_dir = gg.att<string>("tmp_dir");
_output_compress = gg.is("output_compress");
_normalize_always = gg.is("normalize_always");
_bandwidth_limiting = gg.is("bandwidth_limiting");
}
{
auto gg = g.group("probe");
_probe_c5_coeff = gg.att<double>("c5");
_probe_spherical_aberration_coeff = gg.att<double>("cs");
_probe_astigmatism_ca = gg.att<double>("astigmatism_ca");
_probe_astigmatism_angle = gg.att<double>("astigmatism_angle");
_probe_min_aperture = gg.att<double>("min_apert");
_probe_max_aperture = gg.att<double>("max_apert");
_beam_energy = gg.att<double>("beam_energy");
_probe_scan_density = gg.att<double>("scan_density");
if(!cliParamGiven("defocus"))
_probe_defocus = gg.att<unsigned int>("defocus");
_probe_num_defoci = gg.att<unsigned int>("num_defoci");
_probe_fwhm_defoci = gg.att<double>("fwhm_defoci");
}
{
auto gg = g.group("specimen");
_max_potential_radius = gg.att<double>("max_potential_radius");
if(!cliParamGiven("crystal-file"))
_crystal_file = gg.att<string>("crystal_file");
}
{
auto gg = g.group("grating");
_sample_density = gg.att<double>("density");
_slice_thickness = gg.att<double>("slice_thickness");
}
{
auto gg = g.group("adf");
_adf_enabled = gg.is("enabled");
auto adf_scan_x = gg.attList<double>("x");
auto adf_scan_y = gg.attList<double>("y");
_adf_scan_x = make_tuple(adf_scan_x[0], adf_scan_x[1]);
_adf_scan_y = make_tuple(adf_scan_y[0], adf_scan_y[1]);
get<0>(_adf_detector_angles) = gg.att<double>("detector_min_angle");
get<1>(_adf_detector_angles) = gg.att<double>("detector_max_angle");
get<2>(_adf_detector_angles) = gg.att<unsigned int>("detector_num_angles");
_adf_detector_interval_exponent = gg.att<float>("detector_interval_exponent");
_adf_average_configurations = gg.is("average_configurations");
_adf_average_defoci = gg.is("average_defoci");
_adf_save_every_n_slices = gg.att<unsigned int>("save_slices_every");
// when we calculate a single defocus only, set average_defoci to false.
if(_adf_average_defoci && _probe_num_defoci == 1)
_adf_average_defoci = false;
// when defoci are averaged, configurations need to be averaged as well.
if(_adf_average_defoci)
_adf_average_configurations = true;
}
{
auto gg = g.group("cbed");
_cbed_enabled = gg.is("enabled");
auto cbed_scan_x = gg.attList<double>("x");
auto cbed_scan_y = gg.attList<double>("y");
auto cbed_size = gg.attList<unsigned int>("cbed_size");
_cbed_scan_x = make_tuple(cbed_scan_x[0], cbed_scan_x[1]);
_cbed_scan_y = make_tuple(cbed_scan_y[0], cbed_scan_y[1]);
_cbed_size = make_tuple(cbed_size[0], cbed_size[1]);
_cbed_average_configurations = gg.is("average_configurations");
_cbed_average_defoci = gg.is("average_defoci");
_cbed_save_every_n_slices = gg.att<unsigned int>("save_slices_every");
// when we calculate a single defocus only, set average_defoci to false.
if(_cbed_average_defoci && _probe_num_defoci == 1)
_cbed_average_defoci = false;
// when defoci are averaged, configurations need to be averaged as well.
if(_cbed_average_defoci)
_cbed_average_configurations = true;
}
{
auto gg = g.group("frozen_phonon");
if(!cliParamGiven("num-configurations"))
_number_configurations = gg.att<unsigned int>("number_configurations");
_fixed_slicing = gg.is("fixed_slicing");
}
// convert lengths to angstrom and densities to /angstrom
_sample_density /= 10;
_probe_scan_density /= 10;
_slice_thickness *= 10;
_max_potential_radius *= 10;
_probe_c5_coeff *= 10;
_probe_spherical_aberration_coeff *= 10;
_probe_defocus *= 10;
_probe_fwhm_defoci *= 10;
_probe_astigmatism_ca *= 10;
}
......@@ -88,6 +88,16 @@ namespace stemsalabim {
*/
void readParamsFromString(const std::string &prms);
/*!
* Read the simulation parameters from another NC file.
*/
void readParamsFromNCFile();
/*!
* broadcast parameters to other MPI procs.
*/
void broadcast(int rank);
/*!
* Return the path to the crystal XYZ file.
* @return path to crystal XYZ file.
......@@ -128,6 +138,23 @@ namespace stemsalabim {
bool outputCompress() const {
return _output_compress;
}
/*!
* Returns if potentials should be stored and read from the init file.
* @return true if stored potentials should be used.
*/
bool storedPotentials() const {
return _stored_potentials;
}
/*!
* Returns if simulation was initialized from NC file.
* @return true if init from nc file.
*/
bool initFromNCFile() const {
return _init_from_ncfile;
}
/*!
* Return the title of the simulation.
* @return title of the simulation
......@@ -281,14 +308,6 @@ namespace stemsalabim {
return _fixed_slicing;
}
/*!
* Return whether simulation runs with double precision.
* @return true if double precision
*/
bool doublePrecision() const {
return _double_precision;
}
/*!
* Return whether wave functions are bandwidth limited.
* @return true if bandwidth limiting is ON
......@@ -596,6 +615,8 @@ namespace stemsalabim {
double _sample_density{36.0};
bool _output_compress{false};
bool _init_from_ncfile{false};
double _probe_max_aperture{24.0};
double _probe_min_aperture{0.0};
double _probe_c5_coeff{5e7};
......@@ -625,7 +646,8 @@ namespace stemsalabim {
bool _bandwidth_limiting{true};
bool _normalize_always{false};
bool _double_precision{false};
bool _stored_potentials{false};
unsigned int _random_seed{0};
std::mt19937 _rgen{0};
......@@ -638,7 +660,7 @@ namespace stemsalabim {
std::tuple<double, double> _adf_scan_y{std::make_tuple(0.0, 1.0)};
unsigned int _adf_save_every_n_slices{1};
std::tuple<double, double, unsigned int> _adf_detector_angles{std::make_tuple(1.0, 300.0, 300)};
float _adf_detector_interval_exponent{1.0};
double _adf_detector_interval_exponent{1.0};
bool _adf_average_configurations{true};
bool _adf_average_defoci{true};
......
......@@ -98,7 +98,7 @@ void Simulation::init() {
// generate the required random numbers for the frozen phonon stuff.
{
if(!p.initFromNCFile()) {
const size_t number_fp_random_numbers = p.numberOfConfigurations() *
p.numberOfDefoci() *
_crystal->numberOfAtoms() *
......@@ -117,6 +117,13 @@ void Simulation::init() {
mpi_env.broadcast(numbers, mpi::Environment::MASTER);
_fpman->init(numbers);
} else {
vector<float> msds;
if(mpi_env.isMaster()) {
msds = _io->readMSDs();
}
mpi_env.broadcast(msds, mpi::Environment::MASTER);
_fpman->initFromMSDs(msds);
}
_gridman->generateGrids();
......
......@@ -56,6 +56,28 @@ namespace stemsalabim {
*/
void run();
/*!
* Print a line to stdout, with some additional information such as MPI rank, defocus and configuration etc.
* @param st The SimulationState object of the current simulation state.
* @param line the string to print.
*/
static void printLine(const SimulationState &st, const std::string &line);
/*!
* Print a line to stdout, with some additional information such as MPI rank, defocus and configuration etc.
* Line is printed only when this is the master MPI process.
* @param st The SimulationState object of the current simulation state.
* @param line the string to print.
*/
static void printMaster(const SimulationState &st, const std::string &line);
/*!
* Print a line to stdout, with some additional information such as MPI rank, defocus.
* Line is printed only when this is the master MPI process.
* @param line the string to print.
*/
static void printMaster(const std::string &line);
private:
enum STATUS_SIGNAL { START_SIMULATION, START_DEFOCUS, START_CONFIGURATION, PROGRESS, FINISH_SIMULATION };
......@@ -131,28 +153,6 @@ namespace stemsalabim {
*/
void generatePropagator();
/*!
* Print a line to stdout, with some additional information such as MPI rank, defocus and configuration etc.
* @param st The SimulationState object of the current simulation state.
* @param line the string to print.
*/
void printLine(const SimulationState &st, const std::string &line);
/*!
* Print a line to stdout, with some additional information such as MPI rank, defocus and configuration etc.
* Line is printed only when this is the master MPI process.
* @param st The SimulationState object of the current simulation state.
* @param line the string to print.
*/
void printMaster(const SimulationState &st, const std::string &line);
/*!
* Print a line to stdout, with some additional information such as MPI rank, defocus.
* Line is printed only when this is the master MPI process.
* @param line the string to print.
*/
void printMaster(const std::string &line);
/*!
* propagate and bandwidth limit a wave one slice.
* @param wave the wave to propagate.
......
......@@ -39,12 +39,8 @@ void Slice::calculateTransmissionFunction(const shared_ptr<FPConfManager> &fpman
unsigned int samples_x = gridman->samplingX();
unsigned int samples_y = gridman->samplingY();
// see autoslic.cpp line 903ff of Kirklands code
auto vz = calculatePotential(fpman, gridman);
// first of all, define some scales, as we have the real inv_scale (Angstroms) and
// the pixel inv_scale
double scale_x = 1. / gridman->scaleX();
double scale_y = 1. / gridman->scaleY();
// interaction parameter, convert to 1/(eV * A) from 1/(keV * A)
double energy = p.beamEnergy();
......@@ -56,13 +52,44 @@ void Slice::calculateTransmissionFunction(const shared_ptr<FPConfManager> &fpman
(2 * Params::ELECTRON_REST_MASS + energy) /
1000.;
// finish calculation
for(unsigned int current_x = 0; current_x < samples_x; ++current_x) {
for(unsigned int current_y = 0; current_y < samples_y; ++current_y) {
_transmission_function(current_x, current_y) = polar(1.0, sigma * vz[current_x][current_y]);
}
}
// bandwidth limiting
if(p.bandwidthLimiting()) {
_transmission_function.forwardFFT();
_transmission_function *= gridman->getBandwidthMask();
_transmission_function.backwardFFT();
}
}
vector<vector<float>> Slice::calculatePotential(const shared_ptr<FPConfManager> &fpman,
const shared_ptr<GridManager> &gridman) {
Params &p = Params::getInstance();
unsigned int samples_x = gridman->samplingX();
unsigned int samples_y = gridman->samplingY();
// see autoslic.cpp line 903ff of Kirklands code
// first of all, define some scales, as we have the real inv_scale (Angstroms) and
// the pixel inv_scale
double scale_x = 1. / gridman->scaleX();
double scale_y = 1. / gridman->scaleY();
// define a maximal radius above which an atomic potential is zero
// it is ok to have the same number for x and y, presumably, as the error
// is very small when the densities are not equal.
int pix_max = (int) ceil(p.maxPotentialRadius() * p.samplingDensity());
// Scattering calculations
vector<vector<double>> vz(samples_x, vector<double>(samples_y, 0));
vector<vector<float>> vz(samples_x, vector<float>(samples_y, 0));
atomic::Scattering & sf = atomic::Scattering::getInstance();
......@@ -118,18 +145,6 @@ void Slice::calculateTransmissionFunction(const shared_ptr<FPConfManager> &fpman
}
}
// finish calculation
for(unsigned int current_x = 0; current_x < samples_x; ++current_x) {
for(unsigned int current_y = 0; current_y < samples_y; ++current_y) {
_transmission_function(current_x, current_y) = polar(1.0, sigma * vz[current_x][current_y]);
}
}
// bandwidth limiting
if(p.bandwidthLimiting()) {
_transmission_function.forwardFFT();
_transmission_function *= gridman->getBandwidthMask();
_transmission_function.backwardFFT();
}
return vz;
}
......@@ -65,6 +65,14 @@ namespace stemsalabim {
void calculateTransmissionFunction(const std::shared_ptr<FPConfManager> &fpman,
const std::shared_ptr<GridManager> &gridman);
/*!
* Construct a coulomb potential for this slice.
* @param fpman Reference to FPConfManager object
* @param crystal Reference to Crystal object
*/
std::vector<std::vector<float>> calculatePotential(const std::shared_ptr<FPConfManager> &fpman,
const std::shared_ptr<GridManager> &gridman);
/*!
* Clear the cached transmission function.
*/
......
/*
* 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 <iostream>
#include <thread>
#include "classes/Crystal.hpp"
#include "classes/Simulation.hpp"
using namespace std;
using namespace stemsalabim;
#include "utilities/timings.hpp"
using namespace std;
using namespace chrono;
using namespace stemsalabim;
int main(int argc, const char **argv) {
auto &mpi_env = mpi::Environment::getInstance();
if(mpi_env.isSlave()) {
return 0;
}
Params &p = Params::getInstance();
p.initFromCLI(argc, argv);
std::ifstream t(p.paramsFileName());
string params_file_content = std::string((std::istreambuf_iterator<char>(t)), std::istreambuf_iterator<char>());
p.readParamsFromString(params_file_content);
high_resolution_clock::time_point init_start = high_resolution_clock::now();
auto _crystal = std::make_shared<Crystal>();
auto _gridman = std::make_shared<GridManager>(_crystal);
auto _fpman = std::make_shared<FPConfManager>(_crystal);
// read and communicate the content of the crystal file and initialize crystal
// object.
{
ifstream t(p.crystalFilename());
stringstream buffer;
buffer << t.rdbuf();
_crystal->init(buffer.str());
}
// generate the required random numbers for the frozen phonon stuff.
{
const size_t number_fp_random_numbers = p.numberOfConfigurations() *
p.numberOfDefoci() *
_crystal->numberOfAtoms() *
3;
vector<double> numbers;
numbers.resize(number_fp_random_numbers);
#pragma GCC diagnostic push
#pragma GCC diagnostic ignored "-Wmaybe-uninitialized"
normal_distribution<double> distribution(0.0, 1.0);
auto rgen = p.getRandomGenerator();
for(unsigned long i = 0; i < number_fp_random_numbers; ++i)
numbers[i] = distribution(rgen);
#pragma GCC diagnostic pop