diff --git a/src/diffpy/srreal/ObjCrystStructureAdapter.cpp b/src/diffpy/srreal/ObjCrystStructureAdapter.cpp index 2b5e566..8190e81 100644 --- a/src/diffpy/srreal/ObjCrystStructureAdapter.cpp +++ b/src/diffpy/srreal/ObjCrystStructureAdapter.cpp @@ -86,7 +86,9 @@ fetchSymmetryOperations(const ObjCryst::SpaceGroup& spacegroup) assert(nbtran * last <= nbsym); for (int nt = 0; nt < nbtran; ++nt) { - const REAL* pt = sgtrans[nt].tr; + // Keep this compatible with ObjCryst builds where REAL is either + // float or double. + const auto* pt = sgtrans[nt].tr; R3::Vector sgt(pt[0], pt[1], pt[2]); for (int i = 0; i < last; ++i) { diff --git a/src/diffpy/srreal/PDFCalculator.cpp b/src/diffpy/srreal/PDFCalculator.cpp index 47c7189..e2df550 100644 --- a/src/diffpy/srreal/PDFCalculator.cpp +++ b/src/diffpy/srreal/PDFCalculator.cpp @@ -187,10 +187,7 @@ QuantityType PDFCalculator::getExtendedPDF() const QuantityType PDFCalculator::getExtendedRDF() const { QuantityType rdf(this->countExtendedPoints()); - const double& totocc = mstructure_cache.totaloccupancy; - double sfavg = this->sfAverage(); - double rdf_scale = (totocc * sfavg == 0.0) ? 0.0 : - 1.0 / (totocc * sfavg * sfavg); + double rdf_scale = this->getRDFScale(); QuantityType::iterator iirdf = rdf.begin(); QuantityType::const_iterator iival, iival_last; iival = this->value().begin() + @@ -491,9 +488,7 @@ void PDFCalculator::resetValue() // when applicable, configure linear baseline if (this->getBaseline()->type() == "linear") { - double partialpdfscale = - (0.0 == mstructure_cache.totaloccupancy) ? 0.0 : - mstructure_cache.activeoccupancy / mstructure_cache.totaloccupancy; + double partialpdfscale = this->getPartialPDFScale(); double pnumdensity = partialpdfscale * mstructure->numberDensity(); PDFBaseline& bl = *(this->getBaseline()); bl.setDoubleAttr("slope", -4 * M_PI * pnumdensity); @@ -677,6 +672,24 @@ double PDFCalculator::sfAverage() const } +double PDFCalculator::getPartialPDFScale() const +{ + const double totocc = mstructure_cache.totaloccupancy; + return (totocc == 0.0) ? 0.0 : + (mstructure_cache.activeoccupancy / totocc); +} + + +double PDFCalculator::getRDFScale() const +{ + const double& totocc = mstructure_cache.totaloccupancy; + double sfavg = this->sfAverage(); + double rv = (totocc * sfavg == 0.0) ? 0.0 : + 1.0 / (totocc * sfavg * sfavg); + return rv; +} + + void PDFCalculator::cacheStructureData() { int cntsites = this->countSites(); diff --git a/src/diffpy/srreal/PDFCalculator.hpp b/src/diffpy/srreal/PDFCalculator.hpp index 79220e0..5d5eba4 100644 --- a/src/diffpy/srreal/PDFCalculator.hpp +++ b/src/diffpy/srreal/PDFCalculator.hpp @@ -114,6 +114,14 @@ class PDFCalculator : // support for PQEvaluatorOptimized virtual void stashPartialValue(); virtual void restorePartialValue(); + /// activeoccupancy / totaloccupancy used by baseline background term + double getPartialPDFScale() const; + /// RDF scale from the cached total occupancy and average scattering factor + double getRDFScale() const; + /// effective scattering factor at a given site scaled by occupancy + const double& sfSite(int) const; + /// average scattering factor + double sfAverage() const; private: @@ -145,10 +153,6 @@ class PDFCalculator : void cutRipplePoints(QuantityType& y) const; // structure factors - fast lookup by site index - /// effective scattering factor at a given site scaled by occupancy - const double& sfSite(int) const; - /// average scattering factor - double sfAverage() const; void cacheStructureData(); void cacheRlimitsData(); diff --git a/src/diffpy/srreal/R3linalg.cpp b/src/diffpy/srreal/R3linalg.cpp index 19a70cc..337acd4 100644 --- a/src/diffpy/srreal/R3linalg.cpp +++ b/src/diffpy/srreal/R3linalg.cpp @@ -16,6 +16,8 @@ * *****************************************************************************/ +#include +#include #include #include #include @@ -82,6 +84,88 @@ const Matrix& inverse(const Matrix& A) } +void eigen_solve_3x3(const Matrix& A, Vector& w, Matrix& V) +{ + V = identity(); + Matrix m = A; + + const int max_iter = 50; + const double eps = 1e-10; + + for (int iter = 0; iter < max_iter; ++iter) + { + double max_off_diag = 0.0; + int p = 0; + int q = 1; + + for (int i = 0; i < Ndim; ++i) + { + for (int j = i + 1; j < Ndim; ++j) + { + if (std::abs(m(i, j)) > max_off_diag) + { + max_off_diag = std::abs(m(i, j)); + p = i; + q = j; + } + } + } + + if (max_off_diag < eps) break; + + double phi = 0.5 * std::atan2( + 2.0 * m(p, q), m(q, q) - m(p, p)); + double c = std::cos(phi); + double s = std::sin(phi); + + double m_pp = m(p, p); + double m_qq = m(q, q); + double m_pq = m(p, q); + + m(p, p) = c * c * m_pp - 2.0 * s * c * m_pq + s * s * m_qq; + m(q, q) = s * s * m_pp + 2.0 * s * c * m_pq + c * c * m_qq; + m(p, q) = 0.0; + m(q, p) = 0.0; + + for (int i = 0; i < Ndim; ++i) + { + if (i == p || i == q) continue; + double m_ip = m(i, p); + double m_iq = m(i, q); + m(i, p) = c * m_ip - s * m_iq; + m(p, i) = m(i, p); + m(i, q) = s * m_ip + c * m_iq; + m(q, i) = m(i, q); + } + + for (int i = 0; i < Ndim; ++i) + { + double v_ip = V(i, p); + double v_iq = V(i, q); + V(i, p) = c * v_ip - s * v_iq; + V(i, q) = s * v_ip + c * v_iq; + } + } + + w[0] = m(0, 0); + w[1] = m(1, 1); + w[2] = m(2, 2); + + for (int i = 0; i < Ndim - 1; ++i) + { + for (int j = 0; j < Ndim - 1 - i; ++j) + { + if (w[j] <= w[j + 1]) continue; + std::swap(w[j], w[j + 1]); + for (int k = 0; k < Ndim; ++k) + { + std::swap(V(k, j), V(k, j + 1)); + } + } + } +} + + size_t hash_value(const Vector& v) { return boost::hash_range(v.begin(), v.end()); diff --git a/src/diffpy/srreal/R3linalg.hpp b/src/diffpy/srreal/R3linalg.hpp index bffd5c7..9188834 100644 --- a/src/diffpy/srreal/R3linalg.hpp +++ b/src/diffpy/srreal/R3linalg.hpp @@ -150,6 +150,8 @@ const Matrix& identity(); const Matrix& zeromatrix(); double determinant(const Matrix& A); const Matrix& inverse(const Matrix& A); +void eigen_solve_3x3(const Matrix& A, Vector& eigenvalues, + Matrix& eigenvectors); const Vector& floor(const Vector&); template double norm(const V&); diff --git a/src/diffpy/srreal/ThreeDPDFCalculator.cpp b/src/diffpy/srreal/ThreeDPDFCalculator.cpp new file mode 100644 index 0000000..28230c7 --- /dev/null +++ b/src/diffpy/srreal/ThreeDPDFCalculator.cpp @@ -0,0 +1,1239 @@ +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#ifndef M_PI +#define M_PI 3.14159265358979323846 +#endif + +using namespace std; + +namespace diffpy { +namespace srreal { + +// Static constants +constexpr double ThreeDPDFCalculator::DEFAULT_RMAX_3D; +constexpr double ThreeDPDFCalculator::DEFAULT_GRID_STEP; + +// Constructor --------------------------------------------------------------- + +ThreeDPDFCalculator::ThreeDPDFCalculator() : + mdr(DEFAULT_GRID_STEP), + mnbins(0), + maccumblocksize(32), + mapplyrho0background3d(true), + musecqwindow3d(true), + mcalculationmode3d(0), + mhistogramweightmode3d(0), + menablennDelta3d(false), + mnnDelta3d(0.0), + mnnDelta3dUpperBound(0.0), + mnnDeltaPositiveEta3d(0.9), + mdistanceDelta1_3d(0.0), + mdistanceDelta2_3d(0.0), + mdeltaPairA3d("*"), + mdeltaPairB3d("*"), + mdeltaShellIndex3d(1), + mdeltaShellTolerance3d(0.05), + mdeltaKeyTolerance3d(1.0e-6), + museadpscaleSensitivity3d(false), + madpScale3d(1.0), + mrho0backgroundscale3d(1.0) +{ + this->registerDoubleAttribute("enable_nn_delta3d", + this, + &ThreeDPDFCalculator::getEnableNNDelta3DAttr, + &ThreeDPDFCalculator::setEnableNNDelta3DAttr); + this->registerDoubleAttribute("nn_delta3d", + this, + &ThreeDPDFCalculator::getNNDelta3D, + &ThreeDPDFCalculator::setNNDelta3D); + this->registerDoubleAttribute("nn_delta3d_upper_bound", + this, + &ThreeDPDFCalculator::getNNDelta3DUpperBound); + this->registerDoubleAttribute("nn_delta_positive_eta3d", + this, + &ThreeDPDFCalculator::getNNDeltaPositiveEta3D, + &ThreeDPDFCalculator::setNNDeltaPositiveEta3D); + this->registerDoubleAttribute("delta1_3d", + this, + &ThreeDPDFCalculator::getDistanceDelta1_3D, + &ThreeDPDFCalculator::setDistanceDelta1_3D); + this->registerDoubleAttribute("delta2_3d", + this, + &ThreeDPDFCalculator::getDistanceDelta2_3D, + &ThreeDPDFCalculator::setDistanceDelta2_3D); + this->registerDoubleAttribute("delta_shell_index3d", + this, + &ThreeDPDFCalculator::getDeltaShellIndex3DAttr, + &ThreeDPDFCalculator::setDeltaShellIndex3DAttr); + this->registerDoubleAttribute("delta_shell_tolerance3d", + this, + &ThreeDPDFCalculator::getDeltaShellTolerance3D, + &ThreeDPDFCalculator::setDeltaShellTolerance3D); + this->registerDoubleAttribute("delta_key_tolerance3d", + this, + &ThreeDPDFCalculator::getDeltaKeyTolerance3D, + &ThreeDPDFCalculator::setDeltaKeyTolerance3D); + this->registerDoubleAttribute("use_adp_scale_sensitivity3d", + this, + &ThreeDPDFCalculator::getUseADPScaleSensitivity3DAttr, + &ThreeDPDFCalculator::setUseADPScaleSensitivity3DAttr); + this->registerDoubleAttribute("adp_scale3d", + this, + &ThreeDPDFCalculator::getADPScale3D, + &ThreeDPDFCalculator::setADPScale3D); + this->registerDoubleAttribute("rho0_background_scale3d", + this, + &ThreeDPDFCalculator::getRho0BackgroundScale3D, + &ThreeDPDFCalculator::setRho0BackgroundScale3D); + this->registerDoubleAttribute("calculation_mode3d", + this, + &ThreeDPDFCalculator::getCalculationMode3DAttr, + &ThreeDPDFCalculator::setCalculationMode3DAttr); + this->registerDoubleAttribute("histogram_weight_mode3d", + this, + &ThreeDPDFCalculator::getHistogramWeightMode3DAttr, + &ThreeDPDFCalculator::setHistogramWeightMode3DAttr); + + this->setRmax(DEFAULT_RMAX_3D); + this->setRstep(mdr); + this->setQmax(12.0); +} + +// Public Methods ------------------------------------------------------------ + +void ThreeDPDFCalculator::setGridStep(double dr) +{ + if (dr <= 0) throw std::invalid_argument("Grid step must be positive."); + if (dr != mdr) + { + mdr = dr; + this->setRstep(dr); + this->resetValue(); // triggers re-allocation + } +} + +double ThreeDPDFCalculator::getGridStep() const +{ + return mdr; +} + +void ThreeDPDFCalculator::setAccumBlockSize(int bs) +{ + if (bs <= 0) throw std::invalid_argument("Accumulation block size must be positive."); + maccumblocksize = bs; +} + +int ThreeDPDFCalculator::getAccumBlockSize() const +{ + return maccumblocksize; +} + +void ThreeDPDFCalculator::setApplyRho0Background3D(bool v) +{ + mapplyrho0background3d = v; +} + +bool ThreeDPDFCalculator::getApplyRho0Background3D() const +{ + return mapplyrho0background3d; +} + +void ThreeDPDFCalculator::setUseCQWindow3D(bool v) +{ + musecqwindow3d = v; +} + +bool ThreeDPDFCalculator::getUseCQWindow3D() const +{ + return musecqwindow3d; +} + +void ThreeDPDFCalculator::setCalculationMode3D(int mode) +{ + if (mode != 0 && mode != 1) + throw std::invalid_argument("3D calculation mode must be 0 (ADP PDF) or 1 (vector histogram)."); + if (mcalculationmode3d != mode) mticker.click(); + mcalculationmode3d = mode; +} + +int ThreeDPDFCalculator::getCalculationMode3D() const +{ + return mcalculationmode3d; +} + +double ThreeDPDFCalculator::getCalculationMode3DAttr() const +{ + return static_cast(mcalculationmode3d); +} + +void ThreeDPDFCalculator::setCalculationMode3DAttr(double v) +{ + const int mode = static_cast(std::lround(v)); + if (std::fabs(v - static_cast(mode)) > 1.0e-8) + throw std::invalid_argument("3D calculation mode must be an integer value."); + this->setCalculationMode3D(mode); +} + +void ThreeDPDFCalculator::setHistogramWeightMode3D(int mode) +{ + if (mode != 0 && mode != 1) + throw std::invalid_argument("3D histogram weight mode must be 0 (scattering) or 1 (count)."); + if (mhistogramweightmode3d != mode) mticker.click(); + mhistogramweightmode3d = mode; +} + +int ThreeDPDFCalculator::getHistogramWeightMode3D() const +{ + return mhistogramweightmode3d; +} + +double ThreeDPDFCalculator::getHistogramWeightMode3DAttr() const +{ + return static_cast(mhistogramweightmode3d); +} + +void ThreeDPDFCalculator::setHistogramWeightMode3DAttr(double v) +{ + const int mode = static_cast(std::lround(v)); + if (std::fabs(v - static_cast(mode)) > 1.0e-8) + throw std::invalid_argument("3D histogram weight mode must be an integer value."); + this->setHistogramWeightMode3D(mode); +} + +void ThreeDPDFCalculator::setEnableNNDelta3D(bool v) +{ + if (menablennDelta3d != v) mticker.click(); + menablennDelta3d = v; +} + +bool ThreeDPDFCalculator::getEnableNNDelta3D() const +{ + return menablennDelta3d; +} + +double ThreeDPDFCalculator::getEnableNNDelta3DAttr() const +{ + return menablennDelta3d ? 1.0 : 0.0; +} + +void ThreeDPDFCalculator::setEnableNNDelta3DAttr(double v) +{ + this->setEnableNNDelta3D(v != 0.0); +} + +void ThreeDPDFCalculator::setNNDelta3D(double v) +{ + if (v < 0.0) throw std::invalid_argument("3D NN delta must be non-negative."); + const double bound = this->currentNNDeltaUpperBound(); + if (!mdeltaShellRecords.empty()) + { + if (bound <= 0.0 && v > 0.0) + throw std::invalid_argument("3D NN delta cannot be positive without a positive-definite shell bound."); + if (bound > 0.0 && v >= bound) + throw std::invalid_argument("3D NN delta exceeds the positive-definite shell bound."); + } + if (mnnDelta3d != v) mticker.click(); + mnnDelta3d = v; +} + +const double& ThreeDPDFCalculator::getNNDelta3D() const +{ + return mnnDelta3d; +} + +double ThreeDPDFCalculator::getNNDelta3DUpperBound() const +{ + return this->currentNNDeltaUpperBound(); +} + +const double& ThreeDPDFCalculator::getNNDeltaPositiveEta3D() const +{ + return mnnDeltaPositiveEta3d; +} + +void ThreeDPDFCalculator::setNNDeltaPositiveEta3D(double v) +{ + if (v <= 0.0 || v >= 1.0) + throw std::invalid_argument("Positive-definite eta must be between 0 and 1."); + const double minproj = (mnnDeltaPositiveEta3d > 0.0) ? + (mnnDelta3dUpperBound / mnnDeltaPositiveEta3d) : 0.0; + if (mnnDeltaPositiveEta3d != v) mticker.click(); + mnnDeltaPositiveEta3d = v; + if (minproj > 0.0) mnnDelta3dUpperBound = mnnDeltaPositiveEta3d * minproj; + this->validateNNDelta3D(); + this->validateDistanceDecayDelta3D(); +} + +void ThreeDPDFCalculator::setDistanceDelta1_3D(double v) +{ + if (v < 0.0) throw std::invalid_argument("3D distance-delta delta1 must be non-negative."); + if (mdistanceDelta1_3d != v) mticker.click(); + mdistanceDelta1_3d = v; + this->validateDistanceDecayDelta3D(); +} + +const double& ThreeDPDFCalculator::getDistanceDelta1_3D() const +{ + return mdistanceDelta1_3d; +} + +void ThreeDPDFCalculator::setDistanceDelta2_3D(double v) +{ + if (v < 0.0) throw std::invalid_argument("3D distance-delta delta2 must be non-negative."); + if (mdistanceDelta2_3d != v) mticker.click(); + mdistanceDelta2_3d = v; + this->validateDistanceDecayDelta3D(); +} + +const double& ThreeDPDFCalculator::getDistanceDelta2_3D() const +{ + return mdistanceDelta2_3d; +} + +void ThreeDPDFCalculator::setDeltaPairTypes3D(const std::string& a, const std::string& b) +{ + if (a.empty() || b.empty()) + throw std::invalid_argument("Delta pair atom types must be non-empty."); + if (mdeltaPairA3d != a || mdeltaPairB3d != b) mticker.click(); + mdeltaPairA3d = a; + mdeltaPairB3d = b; +} + +const std::string& ThreeDPDFCalculator::getDeltaPairA3D() const +{ + return mdeltaPairA3d; +} + +const std::string& ThreeDPDFCalculator::getDeltaPairB3D() const +{ + return mdeltaPairB3d; +} + +void ThreeDPDFCalculator::setDeltaShellIndex3D(int v) +{ + if (v <= 0) throw std::invalid_argument("Delta shell index must be positive."); + if (mdeltaShellIndex3d != v) mticker.click(); + mdeltaShellIndex3d = v; +} + +int ThreeDPDFCalculator::getDeltaShellIndex3D() const +{ + return mdeltaShellIndex3d; +} + +double ThreeDPDFCalculator::getDeltaShellIndex3DAttr() const +{ + return static_cast(mdeltaShellIndex3d); +} + +void ThreeDPDFCalculator::setDeltaShellIndex3DAttr(double v) +{ + this->setDeltaShellIndex3D(static_cast(std::floor(v + 0.5))); +} + +const double& ThreeDPDFCalculator::getDeltaShellTolerance3D() const +{ + return mdeltaShellTolerance3d; +} + +void ThreeDPDFCalculator::setDeltaShellTolerance3D(double v) +{ + if (v <= 0.0) throw std::invalid_argument("Delta shell tolerance must be positive."); + if (mdeltaShellTolerance3d != v) mticker.click(); + mdeltaShellTolerance3d = v; +} + +const double& ThreeDPDFCalculator::getDeltaKeyTolerance3D() const +{ + return mdeltaKeyTolerance3d; +} + +void ThreeDPDFCalculator::setDeltaKeyTolerance3D(double v) +{ + if (v <= 0.0) throw std::invalid_argument("Delta key tolerance must be positive."); + if (mdeltaKeyTolerance3d != v) mticker.click(); + mdeltaKeyTolerance3d = v; +} + +int ThreeDPDFCalculator::getDeltaEligiblePairCount3D() const +{ + return static_cast(mdeltaEligibleBondKeys.size()); +} + +void ThreeDPDFCalculator::setUseADPScaleSensitivity3D(bool v) +{ + if (museadpscaleSensitivity3d != v) mticker.click(); + museadpscaleSensitivity3d = v; + this->validateNNDelta3D(); +} + +bool ThreeDPDFCalculator::getUseADPScaleSensitivity3D() const +{ + return museadpscaleSensitivity3d; +} + +double ThreeDPDFCalculator::getUseADPScaleSensitivity3DAttr() const +{ + return museadpscaleSensitivity3d ? 1.0 : 0.0; +} + +void ThreeDPDFCalculator::setUseADPScaleSensitivity3DAttr(double v) +{ + this->setUseADPScaleSensitivity3D(v != 0.0); +} + +void ThreeDPDFCalculator::setADPScale3D(double v) +{ + if (v <= 0.0) throw std::invalid_argument("3D ADP scale must be positive."); + if (madpScale3d != v) mticker.click(); + madpScale3d = v; + this->validateNNDelta3D(); +} + +const double& ThreeDPDFCalculator::getADPScale3D() const +{ + return madpScale3d; +} + +void ThreeDPDFCalculator::setRho0BackgroundScale3D(double v) +{ + if (v < 0.0) throw std::invalid_argument("3D rho0 background scale must be non-negative."); + if (mrho0backgroundscale3d != v) mticker.click(); + mrho0backgroundscale3d = v; +} + +const double& ThreeDPDFCalculator::getRho0BackgroundScale3D() const +{ + return mrho0backgroundscale3d; +} + +QuantityType ThreeDPDFCalculator::getThreeDPDF() const +{ + QuantityType result; + std::vector grid = mgrid3d; + if (!this->isVectorHistogramMode3D()) + { + this->applyPostProcessing3D(grid); + } + + size_t nonzero_count = 0; + for (double val : grid) { + if (val != 0.0) ++nonzero_count; + } + result.reserve(nonzero_count * 4); + + for (size_t i = 0; i < grid.size(); ++i) + { + if (grid[i] != 0.0) + { + double x, y, z; + indexToCoord(i, x, y, z); + + result.push_back(x); + result.push_back(y); + result.push_back(z); + result.push_back(grid[i]); + } + } + return result; +} + +QuantityType ThreeDPDFCalculator::getRadialHistogram3D() const +{ + QuantityType result; + result.reserve(mradialhistogram3d.size() * 2); + for (size_t i = 0; i < mradialhistogram3d.size(); ++i) + { + result.push_back(static_cast(i) * mdr); + result.push_back(mradialhistogram3d[i]); + } + return result; +} + +void ThreeDPDFCalculator::exportGrid3DBinary(const std::string& path, bool usefloat32, bool applypost) const +{ + std::vector grid = mgrid3d; + if (applypost) + { + this->applyPostProcessing3D(grid); + } + + std::ofstream ofs(path.c_str(), std::ios::binary | std::ios::trunc); + if (!ofs) throw std::runtime_error("Failed to open output file: " + path); + + const size_t nxy = static_cast(mnbins) * mnbins; + if (usefloat32) + { + std::vector row(nxy); + for (int iz = 0; iz < mnbins; ++iz) + { + const size_t base = static_cast(iz) * nxy; + for (size_t i = 0; i < nxy; ++i) row[i] = static_cast(grid[base + i]); + ofs.write(reinterpret_cast(row.data()), static_cast(nxy * sizeof(float))); + } + } + else + { + for (int iz = 0; iz < mnbins; ++iz) + { + const size_t base = static_cast(iz) * nxy; + ofs.write(reinterpret_cast(&grid[base]), static_cast(nxy * sizeof(double))); + } + } + if (!ofs) throw std::runtime_error("Failed while writing output file: " + path); +} + +void ThreeDPDFCalculator::applyPostProcessing3D(std::vector& grid) const +{ + if (musecqwindow3d) applyQWindow3D(grid); + + const double rdf_scale = this->getRDFScale(); + if (rdf_scale != 1.0) + { + for (double& val : grid) val *= rdf_scale; + } + + if (mapplyrho0background3d) + { + const double rho0_bg = computeRho0Background(); + if (rho0_bg != 0.0) + { + for (double& val : grid) val -= rho0_bg; + } + } + + double qdamp = 0.0; + try + { + qdamp = this->getEnvelopeByType("qresolution")->getDoubleAttr("qdamp"); + } + catch (...) + { + qdamp = 0.0; + } + if (qdamp > 0.0) + { + for (size_t i = 0; i < grid.size(); ++i) + { + if (grid[i] == 0.0) continue; + double x, y, z; + indexToCoord(i, x, y, z); + const double r = sqrt(x * x + y * y + z * z); + grid[i] *= exp(-0.5 * (r * qdamp) * (r * qdamp)); + } + } +} + +double ThreeDPDFCalculator::computeRho0Background() const +{ + const StructureAdapterPtr& structure = this->getStructure(); + if (!structure) return 0.0; + const double partialpdfscale = this->getPartialPDFScale(); + return mrho0backgroundscale3d * partialpdfscale * structure->numberDensity(); +} + +void ThreeDPDFCalculator::applyQWindow3D(std::vector& grid) const +{ + if (grid.empty()) return; + + const double qmin = this->getQmin(); + const double qmax = this->getQmax(); + if (qmin <= 0.0 && qmax <= 0.0) return; + + const int n = mnbins; + int npad = 1; + while (npad < n) npad <<= 1; + + const size_t n3 = static_cast(npad) * npad * npad; + std::vector data(2 * n3, 0.0); + + const int center = n / 2; + for (int iz = 0; iz < n; ++iz) + { + const int sz = (iz - center + npad) % npad; + for (int iy = 0; iy < n; ++iy) + { + const int sy = (iy - center + npad) % npad; + const size_t base_src = static_cast(iz * n + iy) * n; + const size_t base_dst = static_cast(sz * npad + sy) * npad; + for (int ix = 0; ix < n; ++ix) + { + const int sx = (ix - center + npad) % npad; + const size_t src = base_src + ix; + const size_t dst = base_dst + sx; + data[2 * dst] = grid[src]; + } + } + } + + for (int iz = 0; iz < npad; ++iz) + { + for (int iy = 0; iy < npad; ++iy) + { + double* row = &data[2 * ((iz * npad + iy) * npad)]; + gsl_fft_complex_radix2_forward(row, 1, npad); + } + } + + for (int iz = 0; iz < npad; ++iz) + { + for (int ix = 0; ix < npad; ++ix) + { + double* col = &data[2 * (iz * npad * npad + ix)]; + gsl_fft_complex_radix2_forward(col, npad, npad); + } + } + + for (int iy = 0; iy < npad; ++iy) + { + for (int ix = 0; ix < npad; ++ix) + { + double* line = &data[2 * (iy * npad + ix)]; + gsl_fft_complex_radix2_forward(line, npad * npad, npad); + } + } + + const double qstep = (npad > 0) ? (2.0 * M_PI / (npad * mdr)) : 0.0; + for (int iz = 0; iz < npad; ++iz) + { + const int kz = (iz <= npad / 2) ? iz : iz - npad; + const double qz = kz * qstep; + for (int iy = 0; iy < npad; ++iy) + { + const int ky = (iy <= npad / 2) ? iy : iy - npad; + const double qy = ky * qstep; + for (int ix = 0; ix < npad; ++ix) + { + const int kx = (ix <= npad / 2) ? ix : ix - npad; + const double qx = kx * qstep; + const double q = sqrt(qx * qx + qy * qy + qz * qz); + if ((qmax > 0.0 && q > qmax) || (qmin > 0.0 && q < qmin)) + { + const size_t idx = (static_cast(iz) * npad + iy) * npad + ix; + data[2 * idx] = 0.0; + data[2 * idx + 1] = 0.0; + } + } + } + } + + for (int iy = 0; iy < npad; ++iy) + { + for (int ix = 0; ix < npad; ++ix) + { + double* line = &data[2 * (iy * npad + ix)]; + gsl_fft_complex_radix2_inverse(line, npad * npad, npad); + } + } + + for (int iz = 0; iz < npad; ++iz) + { + for (int ix = 0; ix < npad; ++ix) + { + double* col = &data[2 * (iz * npad * npad + ix)]; + gsl_fft_complex_radix2_inverse(col, npad, npad); + } + } + + for (int iz = 0; iz < npad; ++iz) + { + for (int iy = 0; iy < npad; ++iy) + { + double* row = &data[2 * ((iz * npad + iy) * npad)]; + gsl_fft_complex_radix2_inverse(row, 1, npad); + } + } + + for (int iz = 0; iz < n; ++iz) + { + const int sz = (iz - center + npad) % npad; + for (int iy = 0; iy < n; ++iy) + { + const int sy = (iy - center + npad) % npad; + const size_t base_src = static_cast(sz * npad + sy) * npad; + const size_t base_dst = static_cast(iz * n + iy) * n; + for (int ix = 0; ix < n; ++ix) + { + const int sx = (ix - center + npad) % npad; + const size_t src = base_src + sx; + const size_t dst = base_dst + ix; + grid[dst] = data[2 * src]; + } + } + } +} + +// Protected Methods --------------------------------------------------------- + +void ThreeDPDFCalculator::resetValue() +{ + // Ensure odd number of bins so origin is centered. + mnbins = static_cast(2 * ceil(this->getRmax() / mdr)) + 1; + size_t total = static_cast(mnbins) * mnbins * mnbins; + mgrid3d.assign(total, 0.0); + const size_t nr = static_cast(ceil(this->getRmax() / mdr)) + 1; + mradialhistogram3d.assign(nr, 0.0); + + PDFCalculator::resetValue(); + if (!this->isVectorHistogramMode3D()) + { + this->buildDeltaEligibleShellTable(); + this->validateNNDelta3D(); + this->validateDistanceDecayDelta3D(); + } + if (mevaluator) mevaluator->setFlag(USEFULLSUM, true); +} + +void ThreeDPDFCalculator::addPairContribution(const BaseBondGenerator& bnds, int summationscale) +{ + if (bnds.distance() == 0.0) return; + + int i0 = bnds.site0(); + int i1 = bnds.site1(); + int cntsites = this->countSites(); + + if (i0 >= cntsites || i1 >= cntsites) + return; + + const R3::Vector& rvec = bnds.r01(); + const double pairscale = bnds.multiplicity() * static_cast(summationscale); + double sfprod = this->sfSite(i0) * this->sfSite(i1) * pairscale; + + if (this->isVectorHistogramMode3D()) + { + const double weight = (mhistogramweightmode3d == 1) ? pairscale : sfprod; + addVectorHistogramToGrid(rvec, bnds.distance(), weight); + return; + } + + const R3::Matrix& U_i = bnds.Ucartesian0(); + const R3::Matrix& U_j = bnds.Ucartesian1(); + R3::Matrix Sigma = this->effectivePairCovariance(bnds, U_i, U_j); + + addAnisotropicGaussianToGrid(rvec, Sigma, sfprod); +} + +void ThreeDPDFCalculator::addVectorHistogramToGrid(const R3::Vector& r_ij, double distance, double weight) +{ + if (weight == 0.0) return; + + const size_t idx = this->coordToIndex(r_ij[0], r_ij[1], r_ij[2]); + if (idx < mgrid3d.size()) + { +#ifdef _OPENMP +#pragma omp atomic +#endif + mgrid3d[idx] += weight; + } + + if (distance >= 0.0 && !mradialhistogram3d.empty()) + { + const int ir = static_cast(std::lround(distance / mdr)); + if (ir >= 0 && ir < static_cast(mradialhistogram3d.size())) + { +#ifdef _OPENMP +#pragma omp atomic +#endif + mradialhistogram3d[static_cast(ir)] += weight; + } + } +} + +void ThreeDPDFCalculator::addAnisotropicGaussianToGrid(const R3::Vector& r_ij, const R3::Matrix& Sigma, double sfprod) +{ + // Eigenvalue decomposition + R3::Vector eigenvalues; + R3::Matrix eigenvectors; + + R3::eigen_solve_3x3(Sigma, eigenvalues, eigenvectors); + + // Check for positive definiteness (eigenvalues are sorted ascending) + if (eigenvalues[0] <= 1e-8) { + return; + } + + // 3. Invert the covariance matrix + R3::Matrix Sigma_inv = R3::inverse(Sigma); + double det_Sigma = R3::determinant(Sigma); + + // 4. Normalization factor + const double two_pi = 2.0 * M_PI; + double norm_factor = sfprod / (pow(two_pi, 1.5) * sqrt(det_Sigma)); + + // 5. Determine sampling bounding box + double max_sigma = 0.0; + for(int k=0; k<3; ++k) { + max_sigma = std::max(max_sigma, sqrt(eigenvalues[k])); + } + + // 4.0 sigma cutoff + double cutoff_radius = 4.0 * max_sigma; + + // Safety clamps + if (cutoff_radius < mdr) cutoff_radius = mdr; + if (cutoff_radius > 10.0) cutoff_radius = 10.0; + + // 6. Iterate over the local grid indices + double halfspan = (mnbins / 2) * mdr; + + auto get_index_range = [&](double center_val) -> std::pair { + double min_val = center_val - cutoff_radius; + double max_val = center_val + cutoff_radius; + + // Node-centered grid: points at k * md + int start = static_cast(ceil((min_val + halfspan) / mdr)); + int end = static_cast(floor((max_val + halfspan) / mdr)); + + start = std::max(0, start); + end = std::min(mnbins - 1, end); + + return {start, end}; + }; + + std::pair xr = get_index_range(r_ij[0]); + std::pair yr = get_index_range(r_ij[1]); + std::pair zr = get_index_range(r_ij[2]); + + double cutoff_sq = 16.0; + + const int bs = std::max(1, maccumblocksize); + const int nbz = (zr.second - zr.first + bs) / bs; + const int nby = (yr.second - yr.first + bs) / bs; + const int nbx = (xr.second - xr.first + bs) / bs; + +#ifdef _OPENMP +#pragma omp parallel +#endif + { + std::vector local; + +#ifdef _OPENMP +#pragma omp for collapse(3) schedule(dynamic, 1) +#endif + for (int tbz = 0; tbz < nbz; ++tbz) + { + for (int tby = 0; tby < nby; ++tby) + { + for (int tbx = 0; tbx < nbx; ++tbx) + { + const int bz = zr.first + tbz * bs; + const int by = yr.first + tby * bs; + const int bx = xr.first + tbx * bs; + + const int izhi = std::min(zr.second, bz + bs - 1); + const int iyhi = std::min(yr.second, by + bs - 1); + const int ixhi = std::min(xr.second, bx + bs - 1); + + const int tz = izhi - bz + 1; + const int ty = iyhi - by + 1; + const int tx = ixhi - bx + 1; + const size_t nloc = static_cast(tz) * ty * tx; + local.assign(nloc, 0.0); + + bool any = false; + for (int iz = bz; iz <= izhi; ++iz) + { + const double z_grid = iz * mdr - halfspan; + const double dz = z_grid - r_ij[2]; + for (int iy = by; iy <= iyhi; ++iy) + { + const double y_grid = iy * mdr - halfspan; + const double dy = y_grid - r_ij[1]; + for (int ix = bx; ix <= ixhi; ++ix) + { + const double x_grid = ix * mdr - halfspan; + const double dx = x_grid - r_ij[0]; + R3::Vector delta(dx, dy, dz); + R3::Vector tmp = R3::mxvecproduct(Sigma_inv, delta); + const double mahalanobis_sq = R3::dot(delta, tmp); + if (mahalanobis_sq > cutoff_sq) continue; + + const double val = norm_factor * exp(-0.5 * mahalanobis_sq); + const int lz = iz - bz; + const int ly = iy - by; + const int lx = ix - bx; + const size_t lidx = (static_cast(lz) * ty + ly) * tx + lx; + local[lidx] += val; + any = true; + } + } + } + + if (!any) continue; + +#ifdef _OPENMP +#pragma omp critical(threedpdf_tile_reduce) +#endif + { + for (int lz = 0; lz < tz; ++lz) + { + const int iz = bz + lz; + for (int ly = 0; ly < ty; ++ly) + { + const int iy = by + ly; + for (int lx = 0; lx < tx; ++lx) + { + const size_t lidx = (static_cast(lz) * ty + ly) * tx + lx; + const double v = local[lidx]; + if (v == 0.0) continue; + const int ix = bx + lx; + const size_t gidx = (static_cast(iz) * mnbins + iy) * mnbins + ix; + mgrid3d[gidx] += v; + } + } + } + } + } + } + } + } +} + +// Private Helpers ----------------------------------------------------------- + +bool ThreeDPDFCalculator::DeltaBondKey::operator==(const DeltaBondKey& other) const +{ + return site0 == other.site0 && site1 == other.site1 && + rx == other.rx && ry == other.ry && rz == other.rz; +} + +size_t ThreeDPDFCalculator::DeltaBondKeyHash::operator()(const DeltaBondKey& key) const +{ + size_t seed = 0; + seed ^= std::hash()(key.site0) + 0x9e3779b9 + (seed << 6) + (seed >> 2); + seed ^= std::hash()(key.site1) + 0x9e3779b9 + (seed << 6) + (seed >> 2); + seed ^= std::hash()(key.rx) + 0x9e3779b9 + (seed << 6) + (seed >> 2); + seed ^= std::hash()(key.ry) + 0x9e3779b9 + (seed << 6) + (seed >> 2); + seed ^= std::hash()(key.rz) + 0x9e3779b9 + (seed << 6) + (seed >> 2); + return seed; +} + +bool ThreeDPDFCalculator::matchesDeltaPairTypes( + const std::string& atom0, const std::string& atom1) const +{ + const bool allpairs = (mdeltaPairA3d == "*" && mdeltaPairB3d == "*"); + if (allpairs) return true; + const bool forward = + (mdeltaPairA3d == "*" || atom0 == mdeltaPairA3d) && + (mdeltaPairB3d == "*" || atom1 == mdeltaPairB3d); + const bool reverse = + (mdeltaPairA3d == "*" || atom1 == mdeltaPairA3d) && + (mdeltaPairB3d == "*" || atom0 == mdeltaPairB3d); + if (forward || reverse) return true; + return (atom0 == mdeltaPairA3d && atom1 == mdeltaPairB3d) || + (atom0 == mdeltaPairB3d && atom1 == mdeltaPairA3d); +} + +ThreeDPDFCalculator::DeltaBondKey +ThreeDPDFCalculator::makeDeltaBondKey(const BaseBondGenerator& bnds) const +{ + const R3::Vector& r = bnds.r01(); + DeltaBondKey key; + key.site0 = bnds.site0(); + key.site1 = bnds.site1(); + key.rx = static_cast(std::lround(r[0] / mdeltaKeyTolerance3d)); + key.ry = static_cast(std::lround(r[1] / mdeltaKeyTolerance3d)); + key.rz = static_cast(std::lround(r[2] / mdeltaKeyTolerance3d)); + return key; +} + +bool ThreeDPDFCalculator::isDeltaEligiblePair(const BaseBondGenerator& bnds) const +{ + if (!menablennDelta3d || mnnDelta3d == 0.0) return false; + return mdeltaEligibleBondKeys.find(this->makeDeltaBondKey(bnds)) != + mdeltaEligibleBondKeys.end(); +} + +double ThreeDPDFCalculator::projectedSigmaAlongBond( + const R3::Matrix& Sigma, + const R3::Vector& r, + double distance) const +{ + if (distance <= 0.0) return 0.0; + R3::Vector e(r[0] / distance, r[1] / distance, r[2] / distance); + R3::Vector tmp = R3::mxvecproduct(Sigma, e); + return R3::dot(e, tmp); +} + +double ThreeDPDFCalculator::positiveDeltaBoundAlongBond( + const R3::Matrix& Sigma, + const R3::Vector& r, + double distance) const +{ + if (distance <= 0.0) return 0.0; + R3::Vector e(r[0] / distance, r[1] / distance, r[2] / distance); + + R3::Vector eigenvalues; + R3::Matrix eigenvectors; + R3::eigen_solve_3x3(Sigma, eigenvalues, eigenvectors); + if (eigenvalues[0] <= 1.0e-10) + return 0.0; + + try + { + R3::Matrix Sigma_inv = R3::inverse(Sigma); + R3::Vector tmp = R3::mxvecproduct(Sigma_inv, e); + const double denom = R3::dot(e, tmp); + return (std::isfinite(denom) && denom > 0.0) ? (1.0 / denom) : 0.0; + } + catch (...) + { + return 0.0; + } +} + +bool ThreeDPDFCalculator::isDistanceDecayDeltaActive() const +{ + return menablennDelta3d && + (mdistanceDelta1_3d > 0.0 || mdistanceDelta2_3d > 0.0); +} + +double ThreeDPDFCalculator::distanceDecayDeltaFraction(double distance) const +{ + if (distance <= 0.0) return 0.0; + return mdistanceDelta1_3d / distance + + mdistanceDelta2_3d / (distance * distance); +} + +double ThreeDPDFCalculator::currentNNDeltaUpperBound() const +{ + if (mnnDelta3dUpperBound <= 0.0) return 0.0; + return museadpscaleSensitivity3d ? + (mnnDelta3dUpperBound * madpScale3d) : mnnDelta3dUpperBound; +} + +bool ThreeDPDFCalculator::isVectorHistogramMode3D() const +{ + return mcalculationmode3d == 1; +} + +void ThreeDPDFCalculator::validateNNDelta3D() const +{ + if (this->isVectorHistogramMode3D()) return; + if (mnnDelta3d < 0.0) + throw std::invalid_argument("3D NN delta must be non-negative."); + if (mnnDelta3d == 0.0 || mdeltaShellRecords.empty()) return; + const double bound = this->currentNNDeltaUpperBound(); + if (bound <= 0.0) + throw std::invalid_argument("3D NN delta cannot be positive without a positive-definite shell bound."); + if (mnnDelta3d >= bound) + throw std::invalid_argument("3D NN delta exceeds the positive-definite shell bound."); +} + +void ThreeDPDFCalculator::validateDistanceDecayDelta3D() const +{ + if (this->isVectorHistogramMode3D()) return; + if (mdistanceDelta1_3d < 0.0 || mdistanceDelta2_3d < 0.0) + throw std::invalid_argument("3D distance-decay delta parameters must be non-negative."); + if (!this->isDistanceDecayDeltaActive() || mdeltaShellRecords.empty()) return; + + for (const DeltaShellRecord& record : mdeltaShellRecords) + { + if (record.distance <= 0.0) continue; + const double f = this->distanceDecayDeltaFraction(record.distance); + if (f <= 0.0) continue; + if (record.projectedSigma <= 0.0 || record.positiveDeltaBound <= 0.0) + continue; + const double delta = record.projectedSigma * f; + const double bound = mnnDeltaPositiveEta3d * record.positiveDeltaBound; + if (bound <= 0.0 || delta >= bound) + throw std::invalid_argument("3D distance-decay delta violates the pair positive-definite covariance bound."); + } +} + +void ThreeDPDFCalculator::buildDeltaEligibleShellTable() +{ + mdeltaEligibleBondKeys.clear(); + mdeltaShellRecords.clear(); + mnnDelta3dUpperBound = 0.0; + + const StructureAdapterPtr& structure = this->getStructure(); + if (!structure) return; + + BaseBondGeneratorPtr bnds = structure->createBondGenerator(); + this->configureBondGenerator(*bnds); + const int cntsites = structure->countSites(); + + for (int i0 = 0; i0 < cntsites; ++i0) + { + bnds->selectAnchorSite(i0); + bnds->selectSiteRange(0, cntsites); + for (bnds->rewind(); !bnds->finished(); bnds->next()) + { + if (bnds->distance() == 0.0) continue; + const std::string& atom0 = structure->siteAtomType(bnds->site0()); + const std::string& atom1 = structure->siteAtomType(bnds->site1()); + if (!this->matchesDeltaPairTypes(atom0, atom1)) continue; + + const R3::Vector& r = bnds->r01(); + R3::Matrix Sigma = bnds->Ucartesian0() + bnds->Ucartesian1(); + if (museadpscaleSensitivity3d) + { + for (int row = 0; row < 3; ++row) + for (int col = 0; col < 3; ++col) + Sigma(row, col) *= madpScale3d; + } + const double projected = this->projectedSigmaAlongBond( + Sigma, r, bnds->distance()); + const double positive_bound = this->positiveDeltaBoundAlongBond( + Sigma, r, bnds->distance()); + + DeltaShellRecord record; + record.site0 = bnds->site0(); + record.site1 = bnds->site1(); + record.r01x = r[0]; + record.r01y = r[1]; + record.r01z = r[2]; + record.distance = bnds->distance(); + record.projectedSigma = projected; + record.positiveDeltaBound = positive_bound; + mdeltaShellRecords.push_back(record); + } + } + + std::vector order(mdeltaShellRecords.size()); + for (size_t i = 0; i < order.size(); ++i) order[i] = i; + std::sort(order.begin(), order.end(), + [&](size_t a, size_t b) { + return mdeltaShellRecords[a].distance < mdeltaShellRecords[b].distance; + }); + + int shell = 0; + double last_distance = 0.0; + bool first = true; + double min_positive_bound = std::numeric_limits::infinity(); + for (size_t idx : order) + { + const double d = mdeltaShellRecords[idx].distance; + if (first || d - last_distance > mdeltaShellTolerance3d) + { + ++shell; + first = false; + } + last_distance = d; + if (shell == mdeltaShellIndex3d) + { + if (mdeltaShellRecords[idx].positiveDeltaBound < min_positive_bound) + min_positive_bound = mdeltaShellRecords[idx].positiveDeltaBound; + DeltaBondKey key; + key.site0 = mdeltaShellRecords[idx].site0; + key.site1 = mdeltaShellRecords[idx].site1; + key.rx = static_cast(std::lround( + mdeltaShellRecords[idx].r01x / mdeltaKeyTolerance3d)); + key.ry = static_cast(std::lround( + mdeltaShellRecords[idx].r01y / mdeltaKeyTolerance3d)); + key.rz = static_cast(std::lround( + mdeltaShellRecords[idx].r01z / mdeltaKeyTolerance3d)); + mdeltaEligibleBondKeys.insert(key); + } + } + + if (min_positive_bound < std::numeric_limits::infinity() && + min_positive_bound > 0.0) + { + mnnDelta3dUpperBound = mnnDeltaPositiveEta3d * min_positive_bound; + } +} + +R3::Matrix ThreeDPDFCalculator::effectivePairCovariance( + const BaseBondGenerator& bnds, + const R3::Matrix& U_i, + const R3::Matrix& U_j) const +{ + R3::Matrix Sigma = U_i + U_j; + if (museadpscaleSensitivity3d) + { + for (int row = 0; row < 3; ++row) + for (int col = 0; col < 3; ++col) + Sigma(row, col) *= madpScale3d; + } + + const double d = bnds.distance(); + if (d <= 0.0) return Sigma; + + R3::Vector e(bnds.r01()[0] / d, bnds.r01()[1] / d, bnds.r01()[2] / d); + double delta_to_subtract = 0.0; + if (this->isDistanceDecayDeltaActive()) + { + const StructureAdapterPtr& structure = this->getStructure(); + if (!structure) return Sigma; + const std::string& atom0 = structure->siteAtomType(bnds.site0()); + const std::string& atom1 = structure->siteAtomType(bnds.site1()); + if (!this->matchesDeltaPairTypes(atom0, atom1)) return Sigma; + + const double fraction = this->distanceDecayDeltaFraction(d); + if (fraction <= 0.0) return Sigma; + const double projected = this->projectedSigmaAlongBond(Sigma, bnds.r01(), d); + delta_to_subtract = projected * fraction; + } + else + { + if (!this->isDeltaEligiblePair(bnds)) return Sigma; + delta_to_subtract = mnnDelta3d; + } + + if (delta_to_subtract <= 0.0) return Sigma; + + const double positive_bound = this->positiveDeltaBoundAlongBond(Sigma, bnds.r01(), d); + if (positive_bound <= 0.0) return Sigma; + + const double local_bound = mnnDeltaPositiveEta3d * positive_bound; + if (local_bound <= 0.0 || delta_to_subtract >= local_bound) + throw std::invalid_argument("3D delta violates the pair positive-definite covariance bound."); + + for (int row = 0; row < 3; ++row) + { + for (int col = 0; col < 3; ++col) + { + Sigma(row, col) -= delta_to_subtract * e[row] * e[col]; + } + } + return Sigma; +} + +size_t ThreeDPDFCalculator::coordToIndex(double x, double y, double z) const +{ + double halfspan = (mnbins / 2) * mdr; + if (fabs(x) > halfspan || fabs(y) > halfspan || fabs(z) > halfspan) + return mgrid3d.size(); // out of bounds + + int ix = static_cast(std::lround((x + halfspan) / mdr)); + int iy = static_cast(std::lround((y + halfspan) / mdr)); + int iz = static_cast(std::lround((z + halfspan) / mdr)); + + ix = std::max(0, std::min(ix, mnbins - 1)); + iy = std::max(0, std::min(iy, mnbins - 1)); + iz = std::max(0, std::min(iz, mnbins - 1)); + + return static_cast((iz * mnbins + iy) * mnbins + ix); +} + +void ThreeDPDFCalculator::indexToCoord(size_t idx, double& x, double& y, double& z) const +{ + size_t iz = idx / (static_cast(mnbins) * mnbins); + size_t rem = idx % (static_cast(mnbins) * mnbins); + size_t iy = rem / mnbins; + size_t ix = rem % mnbins; + + double halfspan = (mnbins / 2) * mdr; + x = ix * mdr - halfspan; + y = iy * mdr - halfspan; + z = iz * mdr - halfspan; +} + +} // namespace srreal +} // namespace diffpy + +#include +DIFFPY_INSTANTIATE_SERIALIZATION(diffpy::srreal::ThreeDPDFCalculator) diff --git a/src/diffpy/srreal/ThreeDPDFCalculator.hpp b/src/diffpy/srreal/ThreeDPDFCalculator.hpp new file mode 100644 index 0000000..425fcf0 --- /dev/null +++ b/src/diffpy/srreal/ThreeDPDFCalculator.hpp @@ -0,0 +1,233 @@ +#ifndef THREEDPDFCALCULATOR_HPP_INCLUDED +#define THREEDPDFCALCULATOR_HPP_INCLUDED + +#include +#include +#include +#include +#include +#include +#include +#include + +namespace diffpy { +namespace srreal { + +class ThreeDPDFCalculator : public PDFCalculator +{ +public: + // constructor + ThreeDPDFCalculator(); + + // Public interface to retrieve 3D PDF data + // Returns a flat vector of [x0, y0, z0, G0, x1, y1, z1, G1, ...] + QuantityType getThreeDPDF() const; + + // Export dense 3D grid (nz, ny, nx) directly to binary file. + // usefloat32=true writes float32, otherwise float64. + // applypost=true applies q-window/rdf_scale/rho0/qdamp before export. + void exportGrid3DBinary(const std::string& path, bool usefloat32 = true, bool applypost = true) const; + QuantityType getRadialHistogram3D() const; + + // Grid configuration + void setGridStep(double dr); + double getGridStep() const; + + // Blocked accumulation configuration + void setAccumBlockSize(int bs); + int getAccumBlockSize() const; + + void setApplyRho0Background3D(bool); + bool getApplyRho0Background3D() const; + + void setUseCQWindow3D(bool v); + bool getUseCQWindow3D() const; + + // Nearest-neighbor correlated-motion correction for anisotropic 3D PDF models. + // Boolean switches are also exposed as double attributes with 0/1 + // values because libdiffpy's generic attribute system is double-only. + void setEnableNNDelta3D(bool); + bool getEnableNNDelta3D() const; + double getEnableNNDelta3DAttr() const; + void setEnableNNDelta3DAttr(double); + + void setNNDelta3D(double); + const double& getNNDelta3D() const; + double getNNDelta3DUpperBound() const; + const double& getNNDeltaPositiveEta3D() const; + void setNNDeltaPositiveEta3D(double); + void setDistanceDelta1_3D(double); + const double& getDistanceDelta1_3D() const; + void setDistanceDelta2_3D(double); + const double& getDistanceDelta2_3D() const; + + void setDeltaPairTypes3D(const std::string&, const std::string&); + const std::string& getDeltaPairA3D() const; + const std::string& getDeltaPairB3D() const; + void setDeltaShellIndex3D(int); + int getDeltaShellIndex3D() const; + double getDeltaShellIndex3DAttr() const; + void setDeltaShellIndex3DAttr(double); + const double& getDeltaShellTolerance3D() const; + void setDeltaShellTolerance3D(double); + const double& getDeltaKeyTolerance3D() const; + void setDeltaKeyTolerance3D(double); + int getDeltaEligiblePairCount3D() const; + + void setUseADPScaleSensitivity3D(bool); + bool getUseADPScaleSensitivity3D() const; + double getUseADPScaleSensitivity3DAttr() const; + void setUseADPScaleSensitivity3DAttr(double); + void setADPScale3D(double); + const double& getADPScale3D() const; + + void setRho0BackgroundScale3D(double); + const double& getRho0BackgroundScale3D() const; + + void setCalculationMode3D(int); + int getCalculationMode3D() const; + double getCalculationMode3DAttr() const; + void setCalculationMode3DAttr(double); + + void setHistogramWeightMode3D(int); + int getHistogramWeightMode3D() const; + double getHistogramWeightMode3DAttr() const; + void setHistogramWeightMode3DAttr(double); + +protected: + // Override PairQuantity virtual methods + virtual void resetValue() override; + virtual void addPairContribution(const BaseBondGenerator& bnds, int) override; + +private: + struct DeltaBondKey + { + int site0; + int site1; + long rx; + long ry; + long rz; + bool operator==(const DeltaBondKey&) const; + }; + + struct DeltaBondKeyHash + { + size_t operator()(const DeltaBondKey&) const; + }; + + struct DeltaShellRecord + { + int site0; + int site1; + double r01x; + double r01y; + double r01z; + double distance; + double projectedSigma; + double positiveDeltaBound; + }; + + // Helper to add anisotropic Gaussian to the grid + void addAnisotropicGaussianToGrid(const R3::Vector& r_ij, const R3::Matrix& Sigma, double sfprod); + void addVectorHistogramToGrid(const R3::Vector& r_ij, double distance, double weight); + + // Helper: map 3D position to linear index + size_t coordToIndex(double x, double y, double z) const; + void indexToCoord(size_t idx, double& x, double& y, double& z) const; + + void applyPostProcessing3D(std::vector& grid) const; + void applyQWindow3D(std::vector& grid) const; + double computeRho0Background() const; + void buildDeltaEligibleShellTable(); + bool isDeltaEligiblePair(const BaseBondGenerator&) const; + bool matchesDeltaPairTypes(const std::string&, const std::string&) const; + DeltaBondKey makeDeltaBondKey(const BaseBondGenerator&) const; + R3::Matrix effectivePairCovariance( + const BaseBondGenerator&, + const R3::Matrix&, + const R3::Matrix&) const; + double projectedSigmaAlongBond( + const R3::Matrix&, + const R3::Vector&, + double distance) const; + double positiveDeltaBoundAlongBond( + const R3::Matrix&, + const R3::Vector&, + double distance) const; + bool isDistanceDecayDeltaActive() const; + double distanceDecayDeltaFraction(double distance) const; + void validateDistanceDecayDelta3D() const; + void validateNNDelta3D() const; + double currentNNDeltaUpperBound() const; + bool isVectorHistogramMode3D() const; + + // Data members + double mdr; // grid spacing (assumes cubic grid centered at origin) + int mnbins; // number of bins per dimension (odd number, center at 0) + std::vector mgrid3d; // flattened 3D histogram (size = mnbins^3) + std::vector mradialhistogram3d; + int maccumblocksize; // block size for tiled accumulation loops + + bool mapplyrho0background3d; + bool musecqwindow3d; + int mcalculationmode3d; + int mhistogramweightmode3d; + + bool menablennDelta3d; + double mnnDelta3d; + double mnnDelta3dUpperBound; + double mnnDeltaPositiveEta3d; + double mdistanceDelta1_3d; + double mdistanceDelta2_3d; + std::string mdeltaPairA3d; + std::string mdeltaPairB3d; + int mdeltaShellIndex3d; + double mdeltaShellTolerance3d; + double mdeltaKeyTolerance3d; + bool museadpscaleSensitivity3d; + double madpScale3d; + double mrho0backgroundscale3d; + std::unordered_set mdeltaEligibleBondKeys; + std::vector mdeltaShellRecords; + + // Internal constants + static constexpr double DEFAULT_RMAX_3D = 10.0; // Angstrom + static constexpr double DEFAULT_GRID_STEP = 0.1; // Angstrom + + // serialization + friend class boost::serialization::access; + template + void serialize(Archive& ar, const unsigned int version) + { + using boost::serialization::base_object; + ar & base_object(*this); + ar & mdr; + ar & mnbins; + ar & mgrid3d; + ar & mradialhistogram3d; + ar & maccumblocksize; + ar & mapplyrho0background3d; + ar & musecqwindow3d; + ar & mcalculationmode3d; + ar & mhistogramweightmode3d; + ar & menablennDelta3d; + ar & mnnDelta3d; + ar & mnnDelta3dUpperBound; + ar & mnnDeltaPositiveEta3d; + ar & mdistanceDelta1_3d; + ar & mdistanceDelta2_3d; + ar & mdeltaPairA3d; + ar & mdeltaPairB3d; + ar & mdeltaShellIndex3d; + ar & mdeltaShellTolerance3d; + ar & mdeltaKeyTolerance3d; + ar & museadpscaleSensitivity3d; + ar & madpScale3d; + ar & mrho0backgroundscale3d; + } +}; + +} // namespace srreal +} // namespace diffpy + +#endif // THREEDPDFCALCULATOR_HPP_INCLUDED diff --git a/src/tests/TestObjCrystStructureAdapter.hpp b/src/tests/TestObjCrystStructureAdapter.hpp index cd4bf90..fbe8af3 100644 --- a/src/tests/TestObjCrystStructureAdapter.hpp +++ b/src/tests/TestObjCrystStructureAdapter.hpp @@ -19,6 +19,7 @@ * *****************************************************************************/ +#include #include #include @@ -144,6 +145,21 @@ class TestObjCrystStructureAdapter : public CxxTest::TestSuite } + void test_real_translation_vectors() + { + // Ni.cif is F-centered. Its four equivalent atoms require + // fetchSymmetryOperations to read ObjCryst translation vectors, + // including builds where ObjCryst::REAL is double. + CrystalStructureAdapterPtr ni = + std::dynamic_pointer_cast(m_ni); + TS_ASSERT(ni); + TS_ASSERT_EQUALS(4, ni->siteMultiplicity(0)); + CrystalStructureAdapter::AtomVector equivalent = + ni->getEquivalentAtoms(0); + TS_ASSERT_EQUALS(4u, equivalent.size()); + } + + void test_numberDensity() { const double eps = 1.0e-7; diff --git a/src/tests/TestR3linalg.hpp b/src/tests/TestR3linalg.hpp index 8f779df..e84ee17 100644 --- a/src/tests/TestR3linalg.hpp +++ b/src/tests/TestR3linalg.hpp @@ -134,6 +134,53 @@ class TestR3linalg : public CxxTest::TestSuite } + void test_eigen_solve_3x3() + { + R3::Matrix A( + 3.0, 0.2, 0.0, + 0.2, 2.0, 0.1, + 0.0, 0.1, 1.0); + R3::Vector eigenvalues; + R3::Matrix eigenvectors; + R3::eigen_solve_3x3(A, eigenvalues, eigenvectors); + + TS_ASSERT(eigenvalues[0] <= eigenvalues[1]); + TS_ASSERT(eigenvalues[1] <= eigenvalues[2]); + + EpsilonEqual eigclose(1.0e-8); + for (int i = 0; i < R3::Ndim; ++i) + { + R3::Vector v( + eigenvectors(0, i), + eigenvectors(1, i), + eigenvectors(2, i)); + R3::Vector Av = R3::mxvecproduct(A, v); + R3::Vector lv( + eigenvalues[i] * v[0], + eigenvalues[i] * v[1], + eigenvalues[i] * v[2]); + TS_ASSERT(eigclose(Av, lv)); + TS_ASSERT_DELTA(1.0, R3::norm(v), 1.0e-8); + } + + for (int i = 0; i < R3::Ndim; ++i) + { + for (int j = i + 1; j < R3::Ndim; ++j) + { + R3::Vector vi( + eigenvectors(0, i), + eigenvectors(1, i), + eigenvectors(2, i)); + R3::Vector vj( + eigenvectors(0, j), + eigenvectors(1, j), + eigenvectors(2, j)); + TS_ASSERT_DELTA(0.0, R3::dot(vi, vj), 1.0e-8); + } + } + } + + }; // class TestR3linalg // End of file diff --git a/src/tests/TestThreeDPDFCalculator.hpp b/src/tests/TestThreeDPDFCalculator.hpp new file mode 100644 index 0000000..e4836d8 --- /dev/null +++ b/src/tests/TestThreeDPDFCalculator.hpp @@ -0,0 +1,541 @@ +/***************************************************************************** +* +* class TestThreeDPDFCalculator -- unit tests for ThreeDPDFCalculator +* +*****************************************************************************/ + +#include + +#include +#include +#include +#include +#include +#include +#include +#include + +#include +#include +#include +#include + +using namespace std; +using namespace diffpy::srreal; + +// Local Helpers ------------------------------------------------------------- + +namespace { + +AtomicStructureAdapterPtr +makeDimer(const string& atomtype, double uiso) +{ + AtomicStructureAdapterPtr stru = + std::make_shared(); + Atom atom; + atom.atomtype = atomtype; + atom.xyz_cartn = R3::Vector(0.0, 0.0, 0.0); + atom.uij_cartn = R3::zeromatrix(); + atom.uij_cartn(0, 0) = atom.uij_cartn(1, 1) = + atom.uij_cartn(2, 2) = uiso; + stru->append(atom); + atom.xyz_cartn[0] = 1.0; + stru->append(atom); + return stru; +} + + +PeriodicStructureAdapterPtr +makePeriodicDimer() +{ + PeriodicStructureAdapterPtr stru = + std::make_shared(); + stru->setLatPar(10.0, 10.0, 10.0, 90.0, 90.0, 90.0); + AtomicStructureAdapterPtr atoms = makeDimer("Ni", 0.004); + stru->append(atoms->at(0)); + stru->append(atoms->at(1)); + return stru; +} + + +void +configureSmallStandardCalculation(ThreeDPDFCalculator& calc) +{ + calc.setRmax(2.0); + calc.setGridStep(0.25); + calc.setUseCQWindow3D(false); + calc.setApplyRho0Background3D(false); +} + + +bool +allSparseValuesFinite(const QuantityType& data) +{ + if (data.size() % 4) return false; + for (size_t i = 0; i < data.size(); ++i) + { + if (!std::isfinite(data[i])) return false; + } + return true; +} + + +bool +findSparseValue( + const QuantityType& data, + double x, + double y, + double z, + double eps, + double& value) +{ + for (size_t i = 0; i + 3 < data.size(); i += 4) + { + if (std::fabs(data[i] - x) <= eps && + std::fabs(data[i + 1] - y) <= eps && + std::fabs(data[i + 2] - z) <= eps) + { + value = data[i + 3]; + return true; + } + } + return false; +} + + +bool +quantitiesClose( + const QuantityType& a, + const QuantityType& b, + double eps) +{ + if (a.size() != b.size()) return false; + for (size_t i = 0; i < a.size(); ++i) + { + if (std::fabs(a[i] - b[i]) > eps) return false; + } + return true; +} + + +size_t +binaryFileSize(const string& path) +{ + ifstream stream(path.c_str(), ios::binary | ios::ate); + if (!stream) return 0; + return static_cast(stream.tellg()); +} + +} // namespace + +////////////////////////////////////////////////////////////////////////////// +// class TestThreeDPDFCalculator +////////////////////////////////////////////////////////////////////////////// + +class TestThreeDPDFCalculator : public CxxTest::TestSuite +{ + private: + + std::shared_ptr mpdfc; + AtomicStructureAdapterPtr mstru2; + double meps; + const string mfloat32path = "threedpdf_test_float32.bin"; + const string mfloat64path = "threedpdf_test_float64.bin"; + + public: + + void setUp() + { + meps = diffpy::mathutils::SQRT_DOUBLE_EPS; + mpdfc.reset(new ThreeDPDFCalculator); + mstru2 = makeDimer("Ni", 0.004); + std::remove(mfloat32path.c_str()); + std::remove(mfloat64path.c_str()); + } + + + void tearDown() + { + std::remove(mfloat32path.c_str()); + std::remove(mfloat64path.c_str()); + } + + + void test_defaults_and_attributes() + { + TS_ASSERT_DELTA(0.1, mpdfc->getGridStep(), meps); + TS_ASSERT_EQUALS(32, mpdfc->getAccumBlockSize()); + TS_ASSERT_EQUALS(true, mpdfc->getApplyRho0Background3D()); + TS_ASSERT_EQUALS(true, mpdfc->getUseCQWindow3D()); + TS_ASSERT_EQUALS(false, mpdfc->getEnableNNDelta3D()); + TS_ASSERT_EQUALS(0, mpdfc->getCalculationMode3D()); + TS_ASSERT_EQUALS(0, mpdfc->getHistogramWeightMode3D()); + TS_ASSERT_DELTA(0.0, mpdfc->getNNDelta3D(), meps); + TS_ASSERT_DELTA(0.9, mpdfc->getNNDeltaPositiveEta3D(), meps); + TS_ASSERT_DELTA(0.0, mpdfc->getDistanceDelta1_3D(), meps); + TS_ASSERT_DELTA(0.0, mpdfc->getDistanceDelta2_3D(), meps); + TS_ASSERT_EQUALS(string("*"), mpdfc->getDeltaPairA3D()); + TS_ASSERT_EQUALS(string("*"), mpdfc->getDeltaPairB3D()); + TS_ASSERT_EQUALS(1, mpdfc->getDeltaShellIndex3D()); + TS_ASSERT_DELTA( + 0.05, mpdfc->getDeltaShellTolerance3D(), meps); + TS_ASSERT_DELTA( + 1.0e-6, mpdfc->getDeltaKeyTolerance3D(), meps); + TS_ASSERT_EQUALS( + false, mpdfc->getUseADPScaleSensitivity3D()); + TS_ASSERT_DELTA(1.0, mpdfc->getADPScale3D(), meps); + TS_ASSERT_DELTA( + 1.0, mpdfc->getRho0BackgroundScale3D(), meps); + + TS_ASSERT_DELTA(0.0, + mpdfc->getDoubleAttr("enable_nn_delta3d"), meps); + TS_ASSERT_DELTA(0.0, + mpdfc->getDoubleAttr("nn_delta3d"), meps); + TS_ASSERT_DELTA(0.9, + mpdfc->getDoubleAttr("nn_delta_positive_eta3d"), meps); + TS_ASSERT_DELTA(0.0, + mpdfc->getDoubleAttr("delta1_3d"), meps); + TS_ASSERT_DELTA(0.0, + mpdfc->getDoubleAttr("delta2_3d"), meps); + TS_ASSERT_DELTA(1.0, + mpdfc->getDoubleAttr("delta_shell_index3d"), meps); + TS_ASSERT_DELTA(1.0, + mpdfc->getDoubleAttr("adp_scale3d"), meps); + + mpdfc->setDoubleAttr("calculation_mode3d", 1.0); + mpdfc->setDoubleAttr("histogram_weight_mode3d", 1.0); + mpdfc->setDoubleAttr("enable_nn_delta3d", 1.0); + mpdfc->setDoubleAttr("use_adp_scale_sensitivity3d", 1.0); + TS_ASSERT_EQUALS(1, mpdfc->getCalculationMode3D()); + TS_ASSERT_EQUALS(1, mpdfc->getHistogramWeightMode3D()); + TS_ASSERT_EQUALS(true, mpdfc->getEnableNNDelta3D()); + TS_ASSERT_EQUALS( + true, mpdfc->getUseADPScaleSensitivity3D()); + } + + + void test_invalid_configuration() + { + TS_ASSERT_THROWS(mpdfc->setGridStep(0.0), invalid_argument); + TS_ASSERT_THROWS(mpdfc->setGridStep(-0.1), invalid_argument); + TS_ASSERT_THROWS( + mpdfc->setAccumBlockSize(0), invalid_argument); + TS_ASSERT_THROWS( + mpdfc->setCalculationMode3D(2), invalid_argument); + TS_ASSERT_THROWS( + mpdfc->setHistogramWeightMode3D(-1), + invalid_argument); + TS_ASSERT_THROWS( + mpdfc->setDoubleAttr("calculation_mode3d", 0.5), + invalid_argument); + TS_ASSERT_THROWS( + mpdfc->setDoubleAttr("histogram_weight_mode3d", 0.5), + invalid_argument); + TS_ASSERT_THROWS( + mpdfc->setNNDelta3D(-0.01), invalid_argument); + TS_ASSERT_THROWS( + mpdfc->setNNDeltaPositiveEta3D(0.0), + invalid_argument); + TS_ASSERT_THROWS( + mpdfc->setNNDeltaPositiveEta3D(1.0), + invalid_argument); + TS_ASSERT_THROWS( + mpdfc->setDistanceDelta1_3D(-0.01), + invalid_argument); + TS_ASSERT_THROWS( + mpdfc->setDistanceDelta2_3D(-0.01), + invalid_argument); + TS_ASSERT_THROWS( + mpdfc->setDeltaPairTypes3D("", "Ni"), + invalid_argument); + TS_ASSERT_THROWS( + mpdfc->setDeltaShellIndex3D(0), invalid_argument); + TS_ASSERT_THROWS( + mpdfc->setDeltaShellTolerance3D(0.0), + invalid_argument); + TS_ASSERT_THROWS( + mpdfc->setDeltaKeyTolerance3D(0.0), + invalid_argument); + TS_ASSERT_THROWS( + mpdfc->setADPScale3D(0.0), invalid_argument); + TS_ASSERT_THROWS( + mpdfc->setRho0BackgroundScale3D(-1.0), + invalid_argument); + } + + + void test_standard_grid_mode() + { + configureSmallStandardCalculation(*mpdfc); + mpdfc->eval(mstru2); + + QuantityType result = mpdfc->getThreeDPDF(); + TS_ASSERT(!result.empty()); + TS_ASSERT_EQUALS(0u, result.size() % 4); + TS_ASSERT(allSparseValuesFinite(result)); + + double positive = 0.0; + double negative = 0.0; + TS_ASSERT(findSparseValue( + result, 1.0, 0.0, 0.0, meps, positive)); + TS_ASSERT(findSparseValue( + result, -1.0, 0.0, 0.0, meps, negative)); + TS_ASSERT(positive > 0.0); + TS_ASSERT(negative > 0.0); + TS_ASSERT_DELTA(positive, negative, meps); + } + + + void test_vector_histogram_and_pre_rename_regression() + { + mpdfc->setRmax(2.0); + mpdfc->setGridStep(0.5); + mpdfc->setCalculationMode3D(1); + mpdfc->setHistogramWeightMode3D(1); + mpdfc->eval(mstru2); + + QuantityType result = mpdfc->getThreeDPDF(); + TS_ASSERT_EQUALS(8u, result.size()); + double positive = 0.0; + double negative = 0.0; + TS_ASSERT(findSparseValue( + result, 1.0, 0.0, 0.0, meps, positive)); + TS_ASSERT(findSparseValue( + result, -1.0, 0.0, 0.0, meps, negative)); + TS_ASSERT_DELTA(1.0, positive, meps); + TS_ASSERT_DELTA(1.0, negative, meps); + + QuantityType radial = mpdfc->getRadialHistogram3D(); + size_t nr = static_cast( + std::ceil(mpdfc->getRmax() / + mpdfc->getGridStep())) + 1; + TS_ASSERT_EQUALS(2 * nr, radial.size()); + TS_ASSERT_DELTA(2.0, radial[5], meps); + } + + + void test_histogram_mode_does_not_use_adps() + { + AtomicStructureAdapterPtr largeadp = + makeDimer("Ni", 0.4); + mpdfc->setRmax(2.0); + mpdfc->setGridStep(0.5); + mpdfc->setCalculationMode3D(1); + mpdfc->setHistogramWeightMode3D(1); + + mpdfc->eval(mstru2); + QuantityType small = mpdfc->getThreeDPDF(); + mpdfc->eval(largeadp); + QuantityType large = mpdfc->getThreeDPDF(); + TS_ASSERT(quantitiesClose(small, large, 0.0)); + } + + + void test_q_window_changes_grid() + { + configureSmallStandardCalculation(*mpdfc); + mpdfc->eval(mstru2); + QuantityType unwindowed = mpdfc->getThreeDPDF(); + + mpdfc->setUseCQWindow3D(true); + mpdfc->setQmin(0.0); + mpdfc->setQmax(4.0); + QuantityType windowed = mpdfc->getThreeDPDF(); + + TS_ASSERT(!windowed.empty()); + TS_ASSERT(allSparseValuesFinite(windowed)); + TS_ASSERT(!quantitiesClose(unwindowed, windowed, 1.0e-10)); + } + + + void test_rho0_background() + { + PeriodicStructureAdapterPtr periodic = makePeriodicDimer(); + configureSmallStandardCalculation(*mpdfc); + mpdfc->eval(periodic); + QuantityType raw = mpdfc->getThreeDPDF(); + TS_ASSERT(!raw.empty()); + + mpdfc->setApplyRho0Background3D(true); + mpdfc->setRho0BackgroundScale3D(0.75); + QuantityType corrected = mpdfc->getThreeDPDF(); + + double correctedvalue = 0.0; + TS_ASSERT(findSparseValue( + corrected, raw[0], raw[1], raw[2], + meps, correctedvalue)); + const double expected = + 0.75 * periodic->numberDensity(); + TS_ASSERT_DELTA( + expected, raw[3] - correctedvalue, 1.0e-10); + } + + + void test_delta_pair_filter_and_shell_covariance() + { + configureSmallStandardCalculation(*mpdfc); + mpdfc->setEnableNNDelta3D(true); + mpdfc->setDeltaPairTypes3D("O", "O"); + mpdfc->eval(mstru2); + TS_ASSERT_EQUALS(0, mpdfc->getDeltaEligiblePairCount3D()); + + mpdfc->setDeltaPairTypes3D("Ni", "Ni"); + mpdfc->eval(mstru2); + TS_ASSERT( + mpdfc->getDeltaEligiblePairCount3D() > 0); + const double bound = + mpdfc->getNNDelta3DUpperBound(); + TS_ASSERT(bound > 0.0); + QuantityType baseline = mpdfc->getThreeDPDF(); + + mpdfc->setNNDelta3D(0.5 * bound); + mpdfc->eval(mstru2); + QuantityType correlated = mpdfc->getThreeDPDF(); + TS_ASSERT(!quantitiesClose( + baseline, correlated, 1.0e-10)); + TS_ASSERT_THROWS( + mpdfc->setNNDelta3D(bound), invalid_argument); + } + + + void test_distance_delta_covariance() + { + configureSmallStandardCalculation(*mpdfc); + mpdfc->setEnableNNDelta3D(true); + mpdfc->setDeltaPairTypes3D("Ni", "Ni"); + mpdfc->eval(mstru2); + QuantityType baseline = mpdfc->getThreeDPDF(); + + mpdfc->setDistanceDelta1_3D(0.1); + mpdfc->setDistanceDelta2_3D(0.05); + mpdfc->eval(mstru2); + QuantityType correlated = mpdfc->getThreeDPDF(); + TS_ASSERT(!quantitiesClose( + baseline, correlated, 1.0e-10)); + } + + + void test_rdf_normalization_cancels_atom_scattering_scale() + { + AtomicStructureAdapterPtr carbon = + makeDimer("C", 0.004); + configureSmallStandardCalculation(*mpdfc); + mpdfc->eval(mstru2); + QuantityType nickelresult = mpdfc->getThreeDPDF(); + mpdfc->eval(carbon); + QuantityType carbonresult = mpdfc->getThreeDPDF(); + + TS_ASSERT(quantitiesClose( + nickelresult, carbonresult, 1.0e-10)); + } + + + void test_binary_export() + { + mpdfc->setRmax(1.0); + mpdfc->setGridStep(0.5); + mpdfc->setUseCQWindow3D(false); + mpdfc->setApplyRho0Background3D(false); + mpdfc->eval(mstru2); + + mpdfc->exportGrid3DBinary( + mfloat32path, true, false); + mpdfc->exportGrid3DBinary( + mfloat64path, false, false); + const size_t gridcells = 5u * 5u * 5u; + TS_ASSERT_EQUALS( + gridcells * sizeof(float), + binaryFileSize(mfloat32path)); + TS_ASSERT_EQUALS( + gridcells * sizeof(double), + binaryFileSize(mfloat64path)); + TS_ASSERT_THROWS( + mpdfc->exportGrid3DBinary( + "libdiffpy_missing_test_directory/grid.bin"), + runtime_error); + } + + + void test_serialization_round_trip() + { + mpdfc->setRmax(1.0); + mpdfc->setGridStep(0.25); + mpdfc->setAccumBlockSize(7); + mpdfc->setApplyRho0Background3D(false); + mpdfc->setUseCQWindow3D(false); + mpdfc->setCalculationMode3D(1); + mpdfc->setHistogramWeightMode3D(1); + mpdfc->setEnableNNDelta3D(true); + mpdfc->setNNDelta3D(0.02); + mpdfc->setNNDeltaPositiveEta3D(0.8); + mpdfc->setDistanceDelta1_3D(0.03); + mpdfc->setDistanceDelta2_3D(0.04); + mpdfc->setDeltaPairTypes3D("Ni", "Ni"); + mpdfc->setDeltaShellIndex3D(2); + mpdfc->setDeltaShellTolerance3D(0.06); + mpdfc->setDeltaKeyTolerance3D(1.0e-5); + mpdfc->setUseADPScaleSensitivity3D(true); + mpdfc->setADPScale3D(1.2); + mpdfc->setRho0BackgroundScale3D(0.5); + + stringstream storage(ios::in | ios::out | ios::binary); + diffpy::serialization::oarchive oa( + storage, ios::binary); + oa << mpdfc; + diffpy::serialization::iarchive ia( + storage, ios::binary); + std::shared_ptr restored; + ia >> restored; + + TS_ASSERT_DIFFERS(restored.get(), mpdfc.get()); + TS_ASSERT_DELTA( + 0.25, restored->getGridStep(), meps); + TS_ASSERT_EQUALS(7, restored->getAccumBlockSize()); + TS_ASSERT_EQUALS( + false, restored->getApplyRho0Background3D()); + TS_ASSERT_EQUALS( + false, restored->getUseCQWindow3D()); + TS_ASSERT_EQUALS( + 1, restored->getCalculationMode3D()); + TS_ASSERT_EQUALS( + 1, restored->getHistogramWeightMode3D()); + TS_ASSERT_EQUALS( + true, restored->getEnableNNDelta3D()); + TS_ASSERT_DELTA( + 0.02, restored->getNNDelta3D(), meps); + TS_ASSERT_DELTA( + 0.8, restored->getNNDeltaPositiveEta3D(), meps); + TS_ASSERT_DELTA( + 0.03, restored->getDistanceDelta1_3D(), meps); + TS_ASSERT_DELTA( + 0.04, restored->getDistanceDelta2_3D(), meps); + TS_ASSERT_EQUALS( + string("Ni"), restored->getDeltaPairA3D()); + TS_ASSERT_EQUALS( + string("Ni"), restored->getDeltaPairB3D()); + TS_ASSERT_EQUALS( + 2, restored->getDeltaShellIndex3D()); + TS_ASSERT_DELTA( + 0.06, restored->getDeltaShellTolerance3D(), meps); + TS_ASSERT_DELTA( + 1.0e-5, restored->getDeltaKeyTolerance3D(), meps); + TS_ASSERT_EQUALS( + true, restored->getUseADPScaleSensitivity3D()); + TS_ASSERT_DELTA( + 1.2, restored->getADPScale3D(), meps); + TS_ASSERT_DELTA( + 0.5, restored->getRho0BackgroundScale3D(), meps); + + restored->eval(mstru2); + QuantityType restoredresult = + restored->getThreeDPDF(); + mpdfc->eval(mstru2); + QuantityType originalresult = + mpdfc->getThreeDPDF(); + TS_ASSERT(quantitiesClose( + originalresult, restoredresult, 0.0)); + } + +}; // class TestThreeDPDFCalculator + +// End of file