Commit 6666faea authored by Jan Oliver Oelerich's avatar Jan Oliver Oelerich

A seemingly working version of the rewrite. Now comes the Cleanup.... :(

parent fbbe2c3f
......@@ -42,12 +42,15 @@ set(INTERNAL_FILES
utilities/output.hpp
utilities/algorithms.hpp
utilities/TaskQueue.hpp
utilities/timings.hpp)
utilities/timings.hpp
utilities/netcdf_utils.hpp
utilities/ndarray.hpp)
set(EXTERNAL_FILES
3rdparty/jsoncpp/jsoncpp.cpp
3rdparty/jsoncpp/json/json.h
3rdparty/jsoncpp/json/json-forwards.h utilities/netcdf_utils.hpp)
3rdparty/jsoncpp/json/json-forwards.h
3rdparty/multidimgrid/multidimgrid.hpp)
set_source_files_properties(${INTERNAL_FILES}
PROPERTIES COMPILE_FLAGS "-Wall -Wextra -pedantic")
......
#include <memory>
/*
* STEMsalabim: Magical STEM image simulations
*
......@@ -118,20 +120,15 @@ void stemsalabim::Crystal::initFromXYZ(const std::string &field_file_content) {
// set the site's occupant
if(with_slices) {
_atoms.push_back(shared_ptr<atomic::Atom>(new atomic::Atom(x,
y,
z,
msd,
p.elementBySymbol(tokens[0]),
atom_id,
slice)));
_atoms.push_back(std::make_shared<atomic::Atom>(x,
y,
z,
msd,
p.elementBySymbol(tokens[0]),
atom_id,
slice));
} else {
_atoms.push_back(shared_ptr<atomic::Atom>(new atomic::Atom(x,
y,
z,
msd,
p.elementBySymbol(tokens[0]),
atom_id)));
_atoms.push_back(std::make_shared<atomic::Atom>(x, y, z, msd, p.elementBySymbol(tokens[0]), atom_id));
}
_elements.insert(p.elementBySymbol(tokens[0]));
......@@ -140,12 +137,9 @@ void stemsalabim::Crystal::initFromXYZ(const std::string &field_file_content) {
atom_id++;
}
generateSlices();
}
void stemsalabim::Crystal::generateSlices() {
Params &p = Params::getInstance();
......@@ -159,25 +153,7 @@ void stemsalabim::Crystal::generateSlices() {
// slice everything. Now, atomic coordinates with precision as bad as two decimals should be
// taken care of.
unsigned int c = 0;
for(double z = slice_top_offset; z < _size_z; z += dz, c++)
_slices.push_back(shared_ptr<Slice>(new Slice(z, dz, c)));
}
void Crystal::assignFieldsToSlices() {
for(shared_ptr<Slice> &slic: _slices)
slic->clearFields();
size_t i = 0;
double x, y, z, field;
for(tuple<double, double, double, double> &f: _fields) {
// should be improved
std::tie(x, y, z, field) = f;
unsigned int slice_index = sliceIndex(z);
_slices[slice_index]->addField(x, y, field);
i++;
for(double z = slice_top_offset; z < _size_z; z += dz, c++) {
_slices.push_back(std::make_shared<Slice>(z, dz, c));
}
}
......@@ -114,6 +114,26 @@ namespace stemsalabim {
return _elements;
}
/*!
* Add an element to the collection of elements.
*/
void addElement(const std::shared_ptr<atomic::Element> &el) {
_elements.insert(el);
}
/*!
* Add an element to the collection of elements.
*/
void addAtom(const std::shared_ptr<atomic::Atom> &a) {
_atoms.push_back(a);
}
void setLengths(float lx, float ly, float lz) {
_size_x = lx;
_size_y = ly;
_size_z = lz;
}
/*!
* Get the number of slices of the multi-slice simulation
* @return number of slices
......@@ -138,11 +158,6 @@ namespace stemsalabim {
return _slices;
}
/*!
* Assign the field points given in the field file to their corresponding slice.
*/
void assignFieldsToSlices();
/*!
* For a given z coordinate, get the index of the slice that contains this z coordinate.
* Negative z coordinates always result in slice 0, z coordinates that are larger than the
......@@ -167,6 +182,10 @@ namespace stemsalabim {
}
/*!
* Generates the slices
*/
void generateSlices();
private:
std::vector<std::shared_ptr<atomic::Atom>> _atoms;
std::vector<std::shared_ptr<Slice>> _slices;
......@@ -174,10 +193,6 @@ namespace stemsalabim {
std::vector<std::tuple<double, double, double, double>> _fields;
double _size_x, _size_y, _size_z;
/*!
* Generates the slices
*/
void generateSlices();
};
}
......
......@@ -144,6 +144,10 @@ namespace stemsalabim {
return _current_configuration;
}
void setDisplacements(const std::vector<std::vector<std::tuple<double, double, double>>> & d) {
_displacements = d;
}
private:
std::shared_ptr<Crystal> _crystal;
std::vector<std::vector<std::tuple<double, double, double>>> _displacements;
......
This diff is collapsed.
......@@ -30,6 +30,7 @@
#include <netcdf.h>
#include "Params.hpp"
#include "../utilities/memory.hpp"
#include "../utilities/algorithms.hpp"
#include "../utilities/netcdf_utils.hpp"
#include "../utilities/Wave.hpp"
......@@ -54,12 +55,6 @@ namespace stemsalabim {
class IO {
public:
/*!
* Destructor, takes care of closing the NetCDF file and writes the stopping time to it.
* @TODO: Maybe this should be moved to its own function.
*/
~IO();
/*!
* Initialize NetCDF output file. Here, most of the hierarchical structure is written, including
* variables and dimension. Also, compression and chunking is set up and some information about the simulation
......@@ -89,7 +84,7 @@ namespace stemsalabim {
/*!
* Write all simulation parameters to the output file.
*/
void writeParams(const std::shared_ptr<GridManager> &gridman);
void writeParams(const std::shared_ptr<GridManager> &gridman, const std::string & conf_file_string = "");
/*!
* Append to the ADF intensities variable.
......@@ -157,7 +152,18 @@ namespace stemsalabim {
*/
static void dumpWave(std::string filename, const Wave &wave);
std::vector<float> readMSDs();
void initCrystalFromNCFile(const std::shared_ptr<Crystal> &crystal);
void initDisplacementsFromNCFile(const std::shared_ptr<FPConfManager> &fpman);
std::vector<std::vector<float>> readStoredPotential(const std::shared_ptr<FPConfManager> &fpman, unsigned int slice_id);
bool rewriteNCFile() const {
Params &p = Params::getInstance();
return algorithms::absoluteFilePath(p.outputFilename()) != algorithms::absoluteFilePath(p.paramsFileName());
}
void writeRuntime() const;
private:
typedef std::vector<size_t> vs;
......
......@@ -81,7 +81,6 @@ void Params::readParamsFromString(const string &prms) {
_cfg.lookupValue("probe.fwhm_defoci", _probe_fwhm_defoci);
_cfg.lookupValue("specimen.max_potential_radius", _max_potential_radius);
_cfg.lookupValue("specimen.field_file", _field_file);
if(!cliParamGiven("crystal-file"))
_cfg.lookupValue("specimen.crystal_file", _crystal_file);
......@@ -199,33 +198,22 @@ void Params::broadcast(int 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(get<0>(_cbed_scan_x), rank);
mpi_env.broadcast(get<1>(_cbed_scan_x), rank);
mpi_env.broadcast(get<0>(_cbed_scan_y), rank);
mpi_env.broadcast(get<1>(_cbed_scan_y), rank);
mpi_env.broadcast(get<0>(_cbed_size), rank);
mpi_env.broadcast(get<1>(_cbed_size), rank);
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_scan_x), rank);
mpi_env.broadcast(get<1>(_adf_scan_x), rank);
mpi_env.broadcast(get<0>(_adf_scan_y), rank);
mpi_env.broadcast(get<1>(_adf_scan_y), rank);
mpi_env.broadcast(get<0>(_adf_detector_angles), rank);
mpi_env.broadcast(get<1>(_adf_detector_angles), rank);
......@@ -240,7 +228,7 @@ void Params::broadcast(int rank) {
void Params::readParamsFromNCFile() {
NCFile f(_param_file);
NCFile f(_param_file, false, true);
auto g = f.group("params");
_init_from_ncfile = true;
......@@ -253,15 +241,10 @@ void Params::readParamsFromNCFile() {
{
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_file = _param_file;
_title = gg.att("title");
_output_compress = gg.is("output_compress");
_normalize_always = gg.is("normalize_always");
_bandwidth_limiting = gg.is("bandwidth_limiting");
......@@ -281,7 +264,7 @@ void Params::readParamsFromNCFile() {
_probe_scan_density = gg.att<double>("scan_density");
if(!cliParamGiven("defocus"))
_probe_defocus = gg.att<unsigned int>("defocus");
_probe_defocus = gg.att<double>("defocus");
_probe_num_defoci = gg.att<unsigned int>("num_defoci");
_probe_fwhm_defoci = gg.att<double>("fwhm_defoci");
......@@ -294,7 +277,7 @@ void Params::readParamsFromNCFile() {
_max_potential_radius = gg.att<double>("max_potential_radius");
if(!cliParamGiven("crystal-file"))
_crystal_file = gg.att<string>("crystal_file");
_crystal_file = gg.att("crystal_file");
}
......@@ -340,7 +323,7 @@ void Params::readParamsFromNCFile() {
_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");
auto cbed_size = gg.attList<unsigned int>("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]);
......
......@@ -106,14 +106,6 @@ namespace stemsalabim {
return _crystal_file;
}
/*!
* Return the path to the field XYZ file.
* @return path to field XYZ file.
*/
const std::string &fieldFilename() const {
return _field_file;
}
/*!
* Return path to the output file name.
* @return path to output file.
......@@ -545,27 +537,27 @@ namespace stemsalabim {
std::string _uuid{""};
double _sample_density{36.0};
double _sample_density{360};
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};
double _probe_spherical_aberration_coeff{2e4};
double _probe_defocus{-20};
double _probe_c5_coeff{5e6};
double _probe_spherical_aberration_coeff{2e3};
double _probe_defocus{-2};
unsigned int _probe_num_defoci{1};
double _probe_fwhm_defoci{60};
double _probe_fwhm_defoci{6};
double _probe_astigmatism_angle{0};
double _probe_astigmatism_ca{0};
double _probe_scan_density{4}; // in A^-1
double _probe_scan_density{40}; // in A^-1
double _beam_energy{200.0};
double _max_potential_radius{3.0};
double _max_potential_radius{0.3};
double _slice_thickness{2.715};
double _slice_thickness{0.2715};
unsigned int _number_configurations{1};
bool _fixed_slicing{true};
......@@ -573,8 +565,7 @@ namespace stemsalabim {
std::string _crystal_file{""};
std::string _param_file{""};
std::string _output_file{""};
std::string _tmp_dir{""};
std::string _field_file{""};
std::string _tmp_dir{"."};
std::string _title{""};
bool _bandwidth_limiting{true};
......
......@@ -24,6 +24,7 @@
#include <thread>
#include "../libatomic/Scattering.hpp"
#include "../utilities/timings.hpp"
#include "../3rdparty/multidimgrid/multidimgrid.hpp"
using namespace std;
using namespace chrono;
......@@ -40,6 +41,7 @@ void Simulation::init() {
_crystal = std::make_shared<Crystal>();
_gridman = std::make_shared<GridManager>(_crystal);
_fpman = std::make_shared<FPConfManager>(_crystal);
_io = std::make_shared<IO>();
if(mpi_env.isMpi())
output::print("[%s, %d] Hello from %s, MPI rank %d/%d, about to spawn %d threads.\n",
......@@ -50,68 +52,8 @@ void Simulation::init() {
mpi_env.size(),
p.numberOfThreads());
// read and communicate the content of the crystal file and initialize crystal
// object.
{
string crystal_file_content;
if(mpi_env.isMaster()) {
ifstream t(p.crystalFilename());
stringstream buffer;
buffer << t.rdbuf();
crystal_file_content = buffer.str();
}
mpi_env.broadcast(crystal_file_content, mpi::Environment::MASTER);
_crystal->initFromXYZ(crystal_file_content);
}
// read and communicate the field file content, when it is given as a parameter
if(!p.fieldFilename().empty()) {
string field_file_content;
if(mpi_env.isMaster()) {
ifstream t(p.fieldFilename());
stringstream buffer;
buffer << t.rdbuf();
field_file_content = buffer.str();
}
mpi_env.broadcast(field_file_content, mpi::Environment::MASTER);
_crystal->readFieldFile(field_file_content);
_crystal->assignFieldsToSlices();
}
// 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() *
3;
vector<double> numbers;
if(mpi_env.isMaster()) {
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
}
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);
}
_io->initCrystalFromNCFile(_crystal);
_io->initDisplacementsFromNCFile(_fpman);
_gridman->generateGrids();
......@@ -121,14 +63,18 @@ void Simulation::init() {
// initialize output file if this is master process and write some simulation information
// to it.
_io = std::make_shared<IO>();
if(mpi_env.isMaster()) {
if(mpi_env.isMaster() && _io->rewriteNCFile()) {
// As IO should be handled by the master only, we need to initialize
// the instance only when this is the master process.
// The crystal instance will be initialized here and then broadcast to
// tge slaves.
string config_file_contents;
_io->initNcFile(_gridman, _crystal);
_io->writeParams(_gridman);
{
NCFile f(p.paramsFileName(), false, true);
config_file_contents = f.group("params").att("config_file_contents");
}
_io->writeParams(_gridman, config_file_contents);
_io->writeGrids(_gridman);
}
_io->initBuffers(_gridman);
......@@ -142,13 +88,6 @@ void Simulation::init() {
_gridman->adfXGrid().size(),
_gridman->adfYGrid().size()));
// every MPI processor figures out the FFT plans itself. This should somehow be
// improved at some point
//Wave::initPlans(_gridman->samplingX(), _gridman->samplingY());
// make the atomic potentials
//Scattering::init(_crystal);
printMaster(output::fmt("Initialization took %s.\n", output::humantime(algorithms::getTimeSince(init_start))));
initBuffers();
......@@ -165,18 +104,6 @@ void Simulation::run() {
high_resolution_clock::time_point start_sim, t1;
start_sim = high_resolution_clock::now();
// in case the slicing is fixed, this can be done outside of the main loops
if(prms.isFixedSlicing()) {
_fpman->assignAtomsToSlices();
for(shared_ptr<Slice> &slice: _crystal->slices())
printMaster(output::fmt("Slice %d [%.3fnm - %.3fnm] contains %d atoms.\n",
slice->id(),
slice->z(),
(slice->z() + slice->thickness()),
slice->numberOfAtoms()));
}
SimulationStateManager state_manager(_gridman);
// loop over simulations
......@@ -197,53 +124,51 @@ void Simulation::run() {
t1 = high_resolution_clock::now();
// in each configuration, assign atoms to the corresponding slices.
if(!prms.isFixedSlicing()) {
_fpman->assignAtomsToSlices();
printMaster(st, output::fmt("Finished slicing in %s.\n", output::humantime(algorithms::getTimeSince(t1))));
for(shared_ptr<Slice> &slice: _crystal->slices())
printMaster(st,
output::fmt("Slice %d [%.3fnm - %.3fnm] contains %d atoms.\n",
slice->id(),
slice->z(),
(slice->z() + slice->thickness()),
slice->numberOfAtoms()));
}
_fpman->assignAtomsToSlices();
printMaster(st, output::fmt("Finished slicing in %s.\n", output::humantime(algorithms::getTimeSince(t1))));
for(shared_ptr<Slice> &slice: _crystal->slices())
printMaster(st,
output::fmt("Slice %d [%.3fnm - %.3fnm] contains %d atoms.\n",
slice->id(),
slice->z(),
(slice->z() + slice->thickness()),
slice->numberOfAtoms()));
// in each frozen phonon iteration, first write the crystal to the output
// file.
if(mpi_env.isMaster())
if(mpi_env.isMaster() && _io->rewriteNCFile())
_io->writeCrystal(_fpman, _crystal);
// generate transmission functions in parallel using both threads and the main thread.
// This is executed on each MPI processor.
t1 = high_resolution_clock::now();
TaskQueue slice_work;
slice_work.append(_crystal->slices());
if(!prms.storedPotentials()) {
TaskQueue slice_work;
slice_work.append(_crystal->slices());
vector<thread> threads;
for(unsigned int nt = 1; nt < prms.numberOfThreads(); nt++) {
threads.emplace_back([&] {
unsigned int slic_index;
while(slice_work.pop(slic_index)) {
_crystal->slices()[slic_index]->calculateTransmissionFunction(_fpman, _gridman);
slice_work.finish(slic_index);
}
});
}
// also do work on the main thread.
unsigned int slic_index;
while(slice_work.pop(slic_index)) {
_crystal->slices()[slic_index]->calculateTransmissionFunction(_fpman, _gridman);
slice_work.finish(slic_index);
}
vector<thread> threads;
for(unsigned int nt = 0; nt < prms.numberOfThreads(); nt++) {
threads.emplace_back([&] {
unsigned int slic_index;
while(slice_work.pop(slic_index)) {
_crystal->slices()[slic_index]->calculateTransmissionFunction(_fpman, _gridman);
slice_work.finish(slic_index);
}
});
}
slice_work.waitUntilFinished();
for(auto & t : threads)
t.join();
slice_work.waitUntilFinished();
for(auto &t : threads)
t.join();
} else {
for(auto & slic: _crystal->slices()) {
slic->calculateTransmissionFunction(_gridman, _io->readStoredPotential(_fpman, slic->id()));
}
}
printLine(st,
output::fmt("Finished calculation of transmission functions in %s.\n",
......@@ -282,6 +207,10 @@ void Simulation::run() {
};
if(mpi_env.isMaster()) {
_io->writeRuntime();
}
printMaster(output::fmt("Finished computation in %s. Find the results in %s.\n",
output::humantime(algorithms::getTimeSince(start_sim)),
prms.outputFilename()));
......@@ -592,7 +521,7 @@ void Simulation::multisliceMaster(const SimulationState &st) {
}
// join the threads
for(auto & t: threads)
for(auto &t: threads)
t.join();
......@@ -663,7 +592,7 @@ void Simulation::multisliceWorker(const SimulationState &st) {
pixel_index_queue.waitUntilFinished();
for(auto & t: threads)
for(auto &t: threads)
t.join();
threads.clear();
......@@ -748,6 +677,8 @@ void Simulation::generatePropagator() {
}
}
}
}
......
......@@ -26,21 +26,21 @@
#include "Simulation.hpp"
#include "../libatomic/Scattering.hpp"
#include "../utilities/output.hpp"
#include "../3rdparty/multidimgrid/multidimgrid.hpp"
using namespace std;
using namespace stemsalabim;
void Slice::calculateTransmissionFunction(const shared_ptr<FPConfManager> &fpman,
const shared_ptr<GridManager> &gridman) {
void Slice::calculateTransmissionFunction(const shared_ptr<GridManager> &gridman,
const std::vector<std::vector<float>> &vz) {
Params &p = Params::getInstance();
unsigned int samples_x = gridman->samplingX();
unsigned int samples_y = gridman->samplingY();
auto vz = calculatePotential(fpman, gridman);
_transmission_function.init(samples_x, samples_y, complex<float>(1.0, 0.0));
// interaction parameter, convert to 1/(eV * A) from 1/(keV * A)
double energy = p.beamEnergy();
......@@ -68,8 +68,8 @@ void Slice::calculateTransmissionFunction(const shared_ptr<FPConfManager> &fpman
}
vector<vector<float>> Slice::calculatePotential(const shared_ptr<FPConfManager> &fpman,
const shared_ptr<GridManager> &gridman) {
vector<vector<float>>
Slice::calculatePotential(const shared_ptr<FPConfManager> &fpman, const shared_ptr<GridManager> &gridman) {
Params &p = Params::getInstance();
......@@ -91,13 +91,10 @@ vector<vector<float>> Slice::calculatePotential(const shared_ptr<FPConfManager>
// Scattering calculations
vector<vector<float>> vz(samples_x, vector<float>(samples_y, 0));
atomic::Scattering & sf = atomic::Scattering::getInstance();
// initialize wave function from the field cache
_transmission_function.init(samples_x, samples_y, complex<float>(1.0, 0.0));
atomic::Scattering &sf = atomic::Scattering::getInstance();
for(auto &a: _atoms) {
auto & pot_cache = sf.getPotential(a->getElement());
auto &pot_cache = sf.getPotential(a->getElement());
// now add up all the potentials
double dx = fpman->dx(a->getId());
......
......@@ -62,16 +62,21 @@ namespace stemsalabim {
* @param fpman Reference to FPConfManager object
* @param crystal Reference to Crystal object
*/
void calculateTransmissionFunction(const std::shared_ptr<GridManager> &gridman,
const std::vector<std::vector<float>> &pot);
void calculateTransmissionFunction(const std::shared_ptr<FPConfManager> &fpman,
const std::shared_ptr<GridManager> &gridman);
const std::shared_ptr<GridManager> &gridman) {
calculateTransmissionFunction(gridman, calculatePotential(fpman, 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);
std::vector<std::vector<float>>
calculatePotential(const std::shared_ptr<FPConfManager> &fpman, const std::shared_ptr<GridManager> &gridman);
/*!
* Clear the cached transmission function.
......@@ -101,7 +106,7 @@ namespace stemsalabim {
* @param atom
*/
void addField(double x, double y, double field) {