Commit 52ff0c6d authored by Jan Oliver Oelerich's avatar Jan Oliver Oelerich

Moved some code to libatomic and changed calculation of scattering factors to...

Moved some code to libatomic and changed calculation of scattering factors to real space, as this is much cleaner code and doesn't require so much fiddling with ffts.
parent 1cf6638e
......@@ -24,7 +24,7 @@
# compilation
set(INTERNAL_FILES
classes/StructureFactors.hpp
libatomic/Scattering.hpp
classes/Crystal.cpp classes/Crystal.hpp
classes/Params.cpp classes/Params.hpp
classes/Simulation.cpp classes/Simulation.hpp
......@@ -32,11 +32,11 @@ set(INTERNAL_FILES
classes/GridManager.cpp classes/GridManager.hpp
classes/FPConfManager.cpp classes/FPConfManager.hpp
classes/Slice.cpp classes/Slice.hpp
classes/Atom.hpp
classes/Element.hpp
libatomic/Atom.hpp
libatomic/Element.hpp
classes/SimulationState.hpp
utilities/Wave.hpp
utilities/elements_json.hpp
libatomic/elements_json.hpp
utilities/mpi.hpp
utilities/memory.hpp
utilities/output.hpp
......
......@@ -21,10 +21,11 @@
#include "Crystal.hpp"
#include <memory>
#include <regex>
#include "../utilities/output.hpp"
#include "StructureFactors.hpp"
#include "../libatomic/Scattering.hpp"
#include "Slice.hpp"
using namespace std;
......@@ -34,7 +35,7 @@ using namespace stemsalabim;
void stemsalabim::Crystal::init(const std::string &crystal_file_content) {
// this order is very important!
ElementProvider p = ElementProvider::getInstance();
atomic::ElementProvider p = atomic::ElementProvider::getInstance();
// some useful variable declarations
istringstream iss(crystal_file_content);
......@@ -106,7 +107,7 @@ void stemsalabim::Crystal::init(const std::string &crystal_file_content) {
msd = std::stod(tokens[4]) * 100;
// set the site's occupant
_atoms.push_back(shared_ptr<Atom>(new 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]));
atom_id++;
......
......@@ -30,8 +30,8 @@
#include <cmath>
#include <stdlib.h>
#include "Atom.hpp"
#include "Element.hpp"
#include "../libatomic/Atom.hpp"
#include "../libatomic/Element.hpp"
#include "../utilities/algorithms.hpp"
#include "Params.hpp"
......@@ -64,7 +64,7 @@ namespace stemsalabim {
* Returns a reference to the vector of pointers to all Atom objects.
* @return atoms
*/
const std::vector<std::shared_ptr<Atom>> &getAtoms() const {
const std::vector<std::shared_ptr<atomic::Atom>> &getAtoms() const {
return _atoms;
};
......@@ -112,7 +112,7 @@ namespace stemsalabim {
* Get a set of all elements contained in the crystal.
* @return set of elements in the crystal
*/
const std::set<std::shared_ptr<Element>> &elements() const {
const std::set<std::shared_ptr<atomic::Element>> &elements() const {
return _elements;
}
......@@ -170,9 +170,9 @@ namespace stemsalabim {
private:
std::vector<std::shared_ptr<Atom>> _atoms;
std::vector<std::shared_ptr<atomic::Atom>> _atoms;
std::vector<std::shared_ptr<Slice>> _slices;
std::set<std::shared_ptr<Element>> _elements;
std::set<std::shared_ptr<atomic::Element>> _elements;
std::vector<std::tuple<double, double, double, double>> _fields;
double _size_x, _size_y, _size_z;
......
......@@ -61,7 +61,7 @@ void FPConfManager::assignAtomsToSlices() {
slic->clearAtoms();
size_t i = 0;
for(shared_ptr<Atom> a: _crystal->getAtoms()) {
for(shared_ptr<atomic::Atom> a: _crystal->getAtoms()) {
// should be improved
double z = a->getZ() + dz(i);
unsigned int slice_index = _crystal->sliceIndex(z);
......
......@@ -175,6 +175,20 @@ void GridManager::generateGrids() {
for(auto &d: _defocus_weights)
d /= running_sum;
}
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)) {
_bandwidth_limit_mask(ix, iy) = 0.0;
} else {
_bandwidth_limit_mask(ix, iy) = 1.0;
}
}
}
}
}
......
......@@ -24,6 +24,7 @@
#include "Params.hpp"
#include "Crystal.hpp"
#include "../utilities/Wave.hpp"
#include "../utilities/memory.hpp"
namespace stemsalabim {
......@@ -350,6 +351,11 @@ namespace stemsalabim {
return (unsigned int) algorithms::round_even(Params::getInstance().samplingDensity() * _crystal->sizeY());
}
const Wave & getBandwidthMask() const {
return _bandwidth_limit_mask;
}
/*!
* Calculates and returns the k grid values that are actually stored in the output file in kx direction.
* Units: mrad
......@@ -501,6 +507,8 @@ namespace stemsalabim {
std::vector<double> _adf_y_grid;
std::vector<double> _cbed_x_grid;
std::vector<double> _cbed_y_grid;
Wave _bandwidth_limit_mask;
};
}
......
......@@ -179,9 +179,7 @@ putVara(int ncid, int var_id, const vector<size_t> &index, const vector<size_t>
}
void
IO::initNcFile(const shared_ptr<GridManager> &gridman, const shared_ptr<Crystal> &crystal) {
void IO::initNcFile(const shared_ptr<GridManager> &gridman, const shared_ptr<Crystal> &crystal) {
Params &p = Params::getInstance();
auto &mpi_env = mpi::Environment::getInstance();
......@@ -475,7 +473,7 @@ IO::initNcFile(const shared_ptr<GridManager> &gridman, const shared_ptr<Crystal>
// insert atom types
unsigned long i = 0;
for(shared_ptr<Element> el : crystal->elements()) {
for(shared_ptr<atomic::Element> el : crystal->elements()) {
CHK(nc_put_vara_text(amber_id,
atom_types_var,
......@@ -502,10 +500,8 @@ IO::initNcFile(const shared_ptr<GridManager> &gridman, const shared_ptr<Crystal>
}
void
IO::initBuffers(const shared_ptr<GridManager> &gridman, const shared_ptr<Crystal> &crystal) {
void IO::initBuffers(const shared_ptr<GridManager> &gridman) {
Params &p = Params::getInstance();
auto &mpi_env = mpi::Environment::getInstance();
if(p.adf()) {
_adf_intensities_buffer.resize(gridman->adfDetectorGrid().size() * gridman->adfSliceCoords().size());
......@@ -647,8 +643,7 @@ void IO::writeParams(const shared_ptr<GridManager> &gridman) {
}
void
IO::writeCrystal(const shared_ptr<FPConfManager> &fpman, const shared_ptr<Crystal> &crystal) {
void IO::writeCrystal(const shared_ptr<FPConfManager> &fpman, const shared_ptr<Crystal> &crystal) {
int id;
auto &mpi_env = mpi::Environment::getInstance();
......@@ -737,9 +732,8 @@ IO::writeCrystal(const shared_ptr<FPConfManager> &fpman, const shared_ptr<Crysta
}
void IO::writeAdfIntensities(unsigned int idefocus, unsigned int iconf,
const shared_ptr<GridManager> &gridman, const ScanPoint &point,
std::shared_ptr<memory::buffer::number_buffer<float>> &ibuf) {
void IO::writeAdfIntensities(unsigned int idefocus, unsigned int iconf, const shared_ptr<GridManager> &gridman,
const ScanPoint &point, std::shared_ptr<memory::buffer::number_buffer<float>> &ibuf) {
unique_lock<mutex> qlck(_io_lock);
Params &p = Params::getInstance();
......@@ -821,9 +815,8 @@ void IO::writeAdfIntensities(unsigned int idefocus, unsigned int iconf,
}
void IO::writeCBEDIntensities(unsigned int idefocus, unsigned int iconf,
const shared_ptr<GridManager> &gridman, const ScanPoint &point,
std::shared_ptr<memory::buffer::number_buffer<float>> &cbuf) {
void IO::writeCBEDIntensities(unsigned int idefocus, unsigned int iconf, const shared_ptr<GridManager> &gridman,
const ScanPoint &point, std::shared_ptr<memory::buffer::number_buffer<float>> &cbuf) {
unique_lock<mutex> qlck(_io_lock);
Params &p = Params::getInstance();
......@@ -914,8 +907,7 @@ string getTempFileName() {
return output::fmt("%s/%s_%d.bin", p.tmpDir(), p.uuid(), mpi_rank);
}
void IO::writeTemporaryResult(unsigned int idefocus, unsigned int iconf,
ScanPoint &point,
void IO::writeTemporaryResult(unsigned int idefocus, unsigned int iconf, ScanPoint &point,
shared_ptr<memory::buffer::number_buffer<float>> &ibuf,
shared_ptr<memory::buffer::number_buffer<float>> &cbuf) {
......@@ -932,7 +924,7 @@ void IO::writeTemporaryResult(unsigned int idefocus, unsigned int iconf,
// float[] adf_intensities
// float[] cbed_intensities
if (tmpfile.fail()) {
if(tmpfile.fail()) {
output::error("Couldn't open temporary binary %s file for writing.\n", getTempFileName());
}
......@@ -946,12 +938,14 @@ void IO::writeTemporaryResult(unsigned int idefocus, unsigned int iconf,
tmpfile.write(reinterpret_cast<const char *>(&point.index), sizeof(point.index));
if(point.hasAdfIntensities()) {
tmpfile.write(reinterpret_cast<const char *>(ibuf->ptr(point.adf_intensities)), point.adf_intensities.byte_size());
tmpfile.write(reinterpret_cast<const char *>(ibuf->ptr(point.adf_intensities)),
point.adf_intensities.byte_size());
point.clearAdfIntensities(ibuf);
}
if(point.hasCbedIntensities()) {
tmpfile.write(reinterpret_cast<const char *>(cbuf->ptr(point.cbed_intensities)), point.cbed_intensities.byte_size());
tmpfile.write(reinterpret_cast<const char *>(cbuf->ptr(point.cbed_intensities)),
point.cbed_intensities.byte_size());
point.clearCBEDIntensities(cbuf);
}
......@@ -981,7 +975,7 @@ void IO::copyFromTemporaryFile(const shared_ptr<GridManager> &gridman,
unsigned int iconf;
unsigned int pix_idx;
while((int)tmpfile.tellg() < length) {
while((int) tmpfile.tellg() < length) {
tmpfile.read(reinterpret_cast<char *>(&idefocus), sizeof(idefocus));
tmpfile.read(reinterpret_cast<char *>(&iconf), sizeof(iconf));
tmpfile.read(reinterpret_cast<char *>(&pix_idx), sizeof(pix_idx));
......@@ -989,7 +983,7 @@ void IO::copyFromTemporaryFile(const shared_ptr<GridManager> &gridman,
ScanPoint p = gridman->scanPoints()[pix_idx];
if(p.adf) {
auto & entry = ibuf->add_empty();
auto &entry = ibuf->add_empty();
tmpfile.read(reinterpret_cast<char *>(ibuf->ptr(entry)), entry.byte_size());
p.storeAdfIntensities(entry);
writeAdfIntensities(idefocus, iconf, gridman, p, ibuf);
......@@ -997,7 +991,7 @@ void IO::copyFromTemporaryFile(const shared_ptr<GridManager> &gridman,
}
if(p.cbed) {
auto & entry = cbuf->add_empty();
auto &entry = cbuf->add_empty();
tmpfile.read(reinterpret_cast<char *>(cbuf->ptr(entry)), entry.byte_size());
p.storeCBEDIntensities(entry);
writeCBEDIntensities(idefocus, iconf, gridman, p, cbuf);
......
......@@ -67,7 +67,7 @@ namespace stemsalabim {
*/
void initNcFile(const std::shared_ptr<GridManager> &gridman, const std::shared_ptr<Crystal> &crystal);
void initBuffers(const std::shared_ptr<GridManager> &gridman, const std::shared_ptr<Crystal> &crystal);
void initBuffers(const std::shared_ptr<GridManager> &gridman);
/*!
* Appends the atomic positions of the current FP configuration to the AMBER group.
......
......@@ -21,7 +21,7 @@
#include "Simulation.hpp"
#include <memory>
#include "StructureFactors.hpp"
#include "../libatomic/Scattering.hpp"
#include "../utilities/http.hpp"
using namespace std;
......@@ -142,7 +142,7 @@ void Simulation::init() {
_io->writeParams(_gridman);
_io->writeGrids(_gridman);
}
_io->initBuffers(_gridman, _crystal);
_io->initBuffers(_gridman);
printMaster(output::fmt("Calculating on a real and k space lattice of %dx%d points.\n",
_gridman->samplingX(),
......@@ -161,7 +161,7 @@ void Simulation::init() {
//Wave::initPlans(_gridman->samplingX(), _gridman->samplingY());
// make the atomic potentials
//StructureFactors::init(_crystal);
//Scattering::init(_crystal);
if(mpi_env.isMpi())
PRINT_DIAGNOSTICS("finished initialization\n");
......@@ -719,8 +719,7 @@ void Simulation::initBuffers() {
number_cbed_per_pixel * sizeof(float);
_serialization_buffer.resize(prms.workPackageSize() * max_size_work_package + /*Add 10MB*/10 * 1024 * 1024);
_adf_intensity_buffer = std::make_shared<buf_type>(new buf_type(number_intensities_per_pixel,
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());
}
......
......@@ -19,10 +19,12 @@
*
*/
#include <chrono>
#include "Slice.hpp"
#include "Simulation.hpp"
#include "StructureFactors.hpp"
#include "../libatomic/Scattering.hpp"
#include "../utilities/output.hpp"
using namespace std;
......@@ -59,165 +61,16 @@ void Slice::calculateTransmissionFunction(const shared_ptr<FPConfManager> &fpman
// is very small when the densities are not equal.
int pix_max = (int) ceil(p.maxPotentialRadius() * p.samplingDensity());
// StructureFactors calculations
// Scattering calculations
vector<vector<double>> vz(samples_x, vector<double>(samples_y, 0));
// initialize fields
if(_field_cache.empty() && !p.fieldFilename().empty()) {
_field_cache.resize(samples_x, vector<double>(samples_y, 0));
double x, y, field;
int index_x, index_y;
int x1 = -1, y1 = -1, y1_initial = -1;
vector<vector<int>> grid_points_filled(samples_x, vector<int>(samples_y, 0));
for(tuple<double, double, double> &f : _fields) {
tie(x, y, field) = f;
index_x = min(max((int) round(x / scale_x), 0), (int) samples_x - 1);
index_y = min(max((int) round(y / scale_y), 0), (int) samples_y - 1);
// moving average of the field strength in a grid point of one slice.
_field_cache[index_x][index_y] = (_field_cache[index_x][index_y] * grid_points_filled[index_x][index_y] +
field * _thickness) / (grid_points_filled[index_x][index_y] + 1);
grid_points_filled[index_x][index_y]++;
if(x1 < 0 || x1 > index_x)
x1 = index_x;
if(y1 < 0 || y1 > index_y)
y1 = y1_initial = index_y;
}
if(x1 < 0 || y1 < 0)
output::error("The given field has no points in slice [%.2f, %.2f + %.2f].\n", _z, _z, _thickness);
while(true) {
int x2 = x1;
int num_missing_points_x = -1;
do {
num_missing_points_x++;
x2++;
if(x2 >= (int) samples_x)
x2 -= samples_x;
} while(!grid_points_filled[x2][y1]);
if(num_missing_points_x == 0)
break;
y1 = y1_initial;
while(true) {
int y2 = y1;
int num_missing_points_y = -1;
do {
num_missing_points_y++;
y2++;
if(y2 >= (int) samples_y)
y2 -= samples_y;
} while(!grid_points_filled[x1][y2]);
if(num_missing_points_y == 0)
break;
int x2x1, y2y1, x2x, y2y, yy1, xx1;
double q11 = _field_cache[x1][y1];
double q21 = _field_cache[x2][y1];
double q12 = _field_cache[x1][y2];
double q22 = _field_cache[x2][y2];
int i = x1;
do {
int j = y1;
do {
if(!grid_points_filled[i][j]) {
x2x1 = x2 - x1;
y2y1 = y2 - y1;
x2x = x2 - i;
y2y = y2 - j;
xx1 = i - x1;
yy1 = j - y1;
if(x2x1 < 0)
x2x1 += samples_x;
if(y2y1 < 0)
y2y1 += samples_y;
if(x2x < 0)
x2x += samples_x;
if(y2y < 0)
y2y += samples_y;
if(xx1 < 0)
xx1 += samples_x;
if(yy1 < 0)
yy1 += samples_y;
// taken from https://helloacm.com/cc-function-to-compute-the-bilinear-interpolation/
_field_cache[i][j] = 1.0 /
(x2x1 * y2y1) *
(q11 * x2x * y2y +
q21 * xx1 * y2y +
q12 * x2x * yy1 +
q22 * xx1 * yy1);
}
j++;
if(j >= (int) samples_y)
j -= samples_y;
} while(j != y2);
i++;
if(i >= (int) samples_x)
i -= samples_x;
} while(i != x2);
if(y1 < y2)
y1 = y2;
else
break;
}
if(x1 < x2)
x1 = x2;
else
break;
}
vz = _field_cache;
} else if(!_field_cache.empty()) {
vz = _field_cache;
}
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));
// iterate atoms in slice
Wave potential_cache;
// We use a different grating for the calculation of atomic potentials, which is finer than the
// one of the rest of the simulation. This allows us to have distances that suffer less from
// discretization effects. However, as the atomic potential has a divergence at 0 distance, the
// density of the grating influences the potential at the position of the atoms. For that, we
// introduce a minimum radius at which the potential is evaluated, similar to how Kirkland does
// it in his book.
double min_distance = 0.25 * scale_x;
double gscale = StructureFactors::gratingScale();
for(auto &a: _atoms) {
potential_cache = StructureFactors::cache(a->getElement());
// now add up all the potentials
double dx = fpman->dx(a->getId());
double dy = fpman->dy(a->getId());
......@@ -253,21 +106,11 @@ void Slice::calculateTransmissionFunction(const shared_ptr<FPConfManager> &fpman
while(y >= (int) samples_y)
y -= samples_y;
// TODO: occupation as present in kirkland's code is missing.
// the cutoff here is just to ensure circular symmetry
dist_x = abs(current_x * scale_x - (a->getX() + dx));
dist_y = abs(current_y * scale_y - (a->getY() + dy));
distance = (double) sqrt(pow(dist_x, 2) + pow(dist_y, 2));
if(distance < min_distance) {
dist_x = min_distance;
dist_y = 0;
}
dist_x = current_x * scale_x - (a->getX() + dx);
dist_y = current_y * scale_y - (a->getY() + dy);
distance = sqrt(pow(dist_x, 2) + pow(dist_y, 2));
if(distance < p.maxPotentialRadius()) {
vz[x][y] += potential_cache((unsigned int) round(dist_x * gscale),
(unsigned int) round(dist_y * gscale)).real();
}
vz[x][y] += sf.getPotential(distance, a->getElement());
}
}
}
......@@ -282,13 +125,7 @@ void Slice::calculateTransmissionFunction(const shared_ptr<FPConfManager> &fpman
// bandwidth limiting
if(p.bandwidthLimiting()) {
_transmission_function.forwardFFT();
for(unsigned int ix = 0; ix < samples_x; ix++) {
for(unsigned int iy = 0; iy < samples_y; iy++) {
if(pow(gridman->kx(ix), 2) + pow(gridman->ky(iy), 2) > pow(1.0 / (3.0 * max(scale_x, scale_y)), 2)) {
_transmission_function(ix, iy) *= 0.0;
}
}
}
_transmission_function *= gridman->getBandwidthMask();
_transmission_function.backwardFFT();
}
......
......@@ -26,6 +26,7 @@
#include <memory>
#include "../utilities/Wave.hpp"
#include "../libatomic/Atom.hpp"
namespace stemsalabim {
......@@ -34,9 +35,6 @@ namespace stemsalabim {
class FPConfManager;
class Atom;
class GridManager;
/*!
......@@ -86,7 +84,7 @@ namespace stemsalabim {
* Add an atom to this slice.
* @param atom
*/
void addAtom(std::shared_ptr<Atom> &atom) {
void addAtom(std::shared_ptr<atomic::Atom> &atom) {
_atoms.push_back(atom);
}
......@@ -149,7 +147,7 @@ namespace stemsalabim {
* Get a vector of pointers to the Atom objects belonging to this Slice.
* @return Vector of pointers to atoms.
*/
const std::vector<std::shared_ptr<Atom>> &atoms() const {
const std::vector<std::shared_ptr<atomic::Atom>> &atoms() const {
return _atoms;
}
......@@ -158,7 +156,7 @@ namespace stemsalabim {
double _z;
double _thickness;
unsigned int _id;
std::vector<std::shared_ptr<Atom>> _atoms;
std::vector<std::shared_ptr<atomic::Atom>> _atoms;
std::vector<std::tuple<double, double, double>> _fields;
Wave _transmission_function;
std::vector<std::vector<double>> _field_cache;
......
/*
* STEMsalabim: Magical STEM image simulations
*
* Copyright (c) 2016-2018 Jan Oliver Oelerich <[email protected]>
* 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 "StructureFactors.hpp"
#include "Crystal.hpp"
using namespace std;
using namespace stemsalabim;
std::map<std::string, Complex2D> StructureFactors::_s_potential_cache;
double StructureFactors::_s_grating_scale;
complex<double> StructureFactors::feKirkland(const double k, const shared_ptr<Element> &element) {
double sum = 0;
vector<double> scat_prms = element->scatteringParamsDT();
double q2 = k * k;
sum += scat_prms[0] / (q2 + scat_prms[1]);
sum += scat_prms[2] / (q2 + scat_prms[3]);
sum += scat_prms[4] / (q2 + scat_prms[5]);
sum += scat_prms[6] * exp(-scat_prms[7] * q2);
sum += scat_prms[8] * exp(-scat_prms[9] * q2);
sum += scat_prms[10] * exp(-scat_prms[11] * q2);
double factor = PREFACTOR_ATOMPOT / pi;
return complex<double>(sum * factor, 0);
}
void StructureFactors::init(const shared_ptr<Crystal> &crystal) {
Params &p = Params::getInstance();
Complex2D potential_cache;
// We use a different grating for the calculation of atomic potentials, which is finer than the
// one of the rest of the simulation. This allows us to have distances that suffer less from
// discretization effects. However, as the atomic potential has a divergence at 0 distance, the
// density of the grating influences the potential at the position of the atoms. For that, we
// introduce a minimum radius at which the potential is evaluated, similar to how Kirkland does
// it in his book.
// TODO: it would be even more favourable to use an adaptive grid here.
int grating_atompot = 8 * algorithms::round_even(p.samplingDensity() * p.maxPotentialRadius());
double length_atompot = 2 * p.maxPotentialRadius();
_s_grating_scale = grating_atompot / length_atompot;
potential_cache.initPlans((unsigned int) grating_atompot, (unsigned int) grating_atompot, true);
double inv_scale = grating_atompot * grating_atompot / (length_atompot * length_atompot);
for(shared_ptr<Element> e: crystal->elements()) {