From 96e1a7ebc9819c4ce4c42e319a2389af05880428 Mon Sep 17 00:00:00 2001 From: Sparks29032 Date: Fri, 6 Mar 2026 12:22:33 -0800 Subject: [PATCH 1/4] Need this for develop mode building --- SConstruct | 2 +- src/diffpy/srreal/ObjCrystStructureAdapter.cpp | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/SConstruct b/SConstruct index 5179bba5..8bf9920a 100644 --- a/SConstruct +++ b/SConstruct @@ -66,7 +66,7 @@ vars.Add(PathVariable( vars.Add(EnumVariable( 'build', 'compiler settings', - 'fast', allowed_values=('fast', 'debug', 'coverage'))) + 'fast', allowed_values=('fast', 'debug', 'develop', 'coverage'))) vars.Add(EnumVariable( 'tool', 'C++ compiler toolkit to be used', diff --git a/src/diffpy/srreal/ObjCrystStructureAdapter.cpp b/src/diffpy/srreal/ObjCrystStructureAdapter.cpp index f74395ee..f6beea70 100644 --- a/src/diffpy/srreal/ObjCrystStructureAdapter.cpp +++ b/src/diffpy/srreal/ObjCrystStructureAdapter.cpp @@ -86,7 +86,7 @@ fetchSymmetryOperations(const ObjCryst::SpaceGroup& spacegroup) assert(nbtran * last <= nbsym); for (int nt = 0; nt < nbtran; ++nt) { - const double* pt = sgtrans[nt].tr; + const float* pt = sgtrans[nt].tr; R3::Vector sgt(pt[0], pt[1], pt[2]); for (int i = 0; i < last; ++i) { From 0f3e7eef8f2f8aca9694d2733b868bb011f8d13c Mon Sep 17 00:00:00 2001 From: joeseaer <2969111832@qq.com> Date: Fri, 20 Mar 2026 09:41:37 +0800 Subject: [PATCH 2/4] docs: add development instructions for compiling source code --- README.md | 1 + docs/development_instructions.md | 242 +++++++++++++++++++++++++++++++ 2 files changed, 243 insertions(+) create mode 100644 docs/development_instructions.md diff --git a/README.md b/README.md index c227b2b8..4ec67f65 100644 --- a/README.md +++ b/README.md @@ -149,6 +149,7 @@ construction environment can be further customized in a `sconscript.local` script. The library integrity can be verified by executing unit tests with `scons -j4 test` (requires the CxxTest framework). +[Instructions for compiling and modifying source code in development mode](docs/development_instructions.md) ## CONTACTS diff --git a/docs/development_instructions.md b/docs/development_instructions.md new file mode 100644 index 00000000..c609b246 --- /dev/null +++ b/docs/development_instructions.md @@ -0,0 +1,242 @@ +# Compile and Modify diffpy Source Code + +--- + +## I. Prerequisites + +### 1) Install System Dependencies + +Open a terminal and run the following commands in order: + +```bash +# Update package list +sudo apt update + +# Install build toolchain, CMake, Git, Python development headers, etc. +sudo apt install -y build-essential cmake git python3-dev + +# Install SCons (used to build libdiffpy) +sudo apt install -y scons + +# Install GSL (GNU Scientific Library) and ObjCryst (X-ray crystallography library) +sudo apt install -y libgsl-dev libobjcryst-dev +``` + +--- + +### 2) Create a Conda Virtual Environment + +> If you haven't installed Miniconda, please download and install it first + +```bash +# Create your environment (e.g., named your_env_name) with Python 3.11 +conda create -n your_env_name python=3.11 -y + +# Activate the environment (shell prompt will show (your_env_name)) +conda activate your_env_name + +# Install Python dependencies +conda install -y numpy + +# Install additional scientific libraries +conda install -c conda-forge boost-cpp gsl pyobjcryst +``` + +All subsequent operations should be performed within the (your_env_name) environment! + +--- + +## II. Get the Source Code + +Assume our working directory is `~/workspace` (adjust to your actual setup): + +```bash +# Switch to the working directory +cd ~/workspace + +# Clone the two core repositories +git clone https://github.com/diffpy/libdiffpy.git diffpy.libdiffpy +git clone https://github.com/diffpy/diffpy.srreal.git diffpy.srreal +``` + +> The directory structure should be: `~/workspace/` +> +> ``` +> ~/workspace/ +> ├── diffpy.libdiffpy/ ← C++ core library +> └── diffpy.srreal/ ← Python interface layer +> ``` + +--- + +## III. Initial Build (Developer Mode) + +> We use `--inplace` so after modifying code you only need to recompile, no `pip install` required. + +--- + +### Step 1: Build libdiffpy (C++ Library) + +```bash +# Enter the libdiffpy directory +cd ~/workspace/diffpy.libdiffpy + +# Clean old builds (optional but recommended) +rm -rf build/ + +# Build the C++ library with SCons +scons +``` + +Success indicators: + +- `build/fast-x86_64/libdiffpy.so` is generated +- No errors in terminal output; warnings can be ignored + +Verify it was generated: + +```bash +ls -l build/fast-x86_64/libdiffpy.so +``` + +> If you see `ObjCryst not found` or `gsl-config not found`, make sure you ran the `apt install` commands in Step 1. + +--- + +### Step 2: Set Environment Variables + +To allow srreal to find the newly built `libdiffpy.so`, set the following variables and replace the paths with your actual paths: + +```bash +# Set the libdiffpy build path +export DIFFPY_LIBDIFFPY_BUILD="$HOME/workspace/diffpy.libdiffpy/build/fast-x86_64" + +# Tell the linker and runtime where to find .so files +export LD_LIBRARY_PATH="$DIFFPY_LIBDIFFPY_BUILD:$LD_LIBRARY_PATH" +export LIBRARY_PATH="$DIFFPY_LIBDIFFPY_BUILD:$LIBRARY_PATH" + +# Let Python import srreal source directly (developer mode) +export PYTHONPATH="$HOME/workspace/diffpy.srreal/src:$PYTHONPATH" +``` + +--- + +### Step 3: Build diffpy.srreal (Python Extension) + +```bash +# Enter the srreal directory +cd ~/workspace/diffpy.srreal + +# Clean old builds +rm -rf build/ + +# Build the Python extension (generate .so in place) +python setup.py build_ext --inplace +``` + +Success indicators: + +- `diffpy/srreal/srreal_ext.cpython-*.so` is generated +- No `cannot find -ldiffpy` errors + +--- + +### Step 4: Verify the Installation + +```bash +# Test whether it can be imported +python -c "from diffpy.srreal.pdfcalculator import PDFCalculator; print('srreal Success！')" +``` + +If you see `srreal Success！`, everything is working! + +--- + +## IV. Daily Development: Modify Code + Rebuild + +### Scenario A: Only Modified diffpy.srreal Python or C++ Code + +```bash +# 1. Activate your Conda environment +conda activate your_env_name + +# 2. Enter the srreal directory +cd ~/workspace/diffpy.srreal + +# 3. Clean and rebuild +rm -rf build/ +python setup.py build_ext --inplace +``` + +No need to rebuild libdiffpy + +--- + +### Scenario B: Modified diffpy.libdiffpy C++ Code + +```bash +# 1. Activate your Conda environment +conda activate your_env_name + +# 2. Rebuild libdiffpy +cd ~/workspace/diffpy.libdiffpy +rm -rf build/ +scons + +# 3. Update environment variables (ensure they point to the latest .so) +export DIFFPY_LIBDIFFPY_BUILD="$HOME/workspace/diffpy.libdiffpy/build/fast-x86_64" +export LD_LIBRARY_PATH="$DIFFPY_LIBDIFFPY_BUILD:$LD_LIBRARY_PATH" + +# 4. Rebuild srreal (since it depends on libdiffpy) +cd ~/workspace/diffpy.srreal +rm -rf build/ +python setup.py build_ext --inplace +``` + +Key point: if you change the underlying C++ library, you must rebuild the Python extension + +--- + +## V. Common Issues and Fixes + +### Issue 1: ModuleNotFoundError: No module named 'numpy' + +→ Fix: run `conda install numpy` in your Conda environment + +### Issue 2: cannot find -ldiffpy + +→ Fix: + +1. Confirm `libdiffpy.so` has been generated +2. Confirm `LD_LIBRARY_PATH` and `LIBRARY_PATH` are set correctly + +### Issue 3: gsl-config: not found + +→ Fix: run `sudo apt install libgsl-dev` + +### Issue 4: Checking for C++ library ObjCryst... no + +→ Fix: run `sudo apt install libobjcryst-dev` + +If you encounter other issues, you can search relevant forums or contact qiaohai@tongji.edu.cn + +--- + +## VI. Appendix: One-Click Environment Variable Setup + +You can append the following to `~/.bashrc` so you don't have to set environment variables every time: + +```bash +# diffpy development environment +export DIFFPY_WORKSPACE="$HOME/workspace" +export DIFFPY_LIBDIFFPY_BUILD="$DIFFPY_WORKSPACE/diffpy.libdiffpy/build/fast-x86_64" +export LD_LIBRARY_PATH="$DIFFPY_LIBDIFFPY_BUILD:$LD_LIBRARY_PATH" +export LIBRARY_PATH="$DIFFPY_LIBDIFFPY_BUILD:$LIBRARY_PATH" +export PYTHONPATH="$DIFFPY_WORKSPACE/diffpy.srreal/src:$PYTHONPATH" +``` + +Then run: + +```bash +source ~/.bashrc +``` From ee254fa0c7c077780170f3b14b6ab80e03b41368 Mon Sep 17 00:00:00 2001 From: joeseaer <2969111832@qq.com> Date: Sun, 5 Jul 2026 14:30:56 +0800 Subject: [PATCH 3/4] ENH:Add 3D PDF calculator support --- .../srreal/ObjCrystStructureAdapter.cpp | 2 +- src/diffpy/srreal/PDF3DCalculator.cpp | 1308 +++++++++++++++++ src/diffpy/srreal/PDF3DCalculator.hpp | 257 ++++ src/diffpy/srreal/PDFCalculator.cpp | 8 + src/diffpy/srreal/PDFCalculator.hpp | 2 + src/diffpy/srreal/R3linalg.cpp | 84 ++ src/diffpy/srreal/R3linalg.hpp | 2 + 7 files changed, 1662 insertions(+), 1 deletion(-) create mode 100644 src/diffpy/srreal/PDF3DCalculator.cpp create mode 100644 src/diffpy/srreal/PDF3DCalculator.hpp diff --git a/src/diffpy/srreal/ObjCrystStructureAdapter.cpp b/src/diffpy/srreal/ObjCrystStructureAdapter.cpp index f6beea70..f74395ee 100644 --- a/src/diffpy/srreal/ObjCrystStructureAdapter.cpp +++ b/src/diffpy/srreal/ObjCrystStructureAdapter.cpp @@ -86,7 +86,7 @@ fetchSymmetryOperations(const ObjCryst::SpaceGroup& spacegroup) assert(nbtran * last <= nbsym); for (int nt = 0; nt < nbtran; ++nt) { - const float* pt = sgtrans[nt].tr; + const double* 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/PDF3DCalculator.cpp b/src/diffpy/srreal/PDF3DCalculator.cpp new file mode 100644 index 00000000..e839ab35 --- /dev/null +++ b/src/diffpy/srreal/PDF3DCalculator.cpp @@ -0,0 +1,1308 @@ +/************************************************************************************* +* +* libdiffpy by DANSE Diffraction group +* Simon J. L. Billinge +* (c) 2009 The Trustees of Columbia University +* in the City of New York. All rights reserved. +* +* File coded by: Hai Qiao +* +* See AUTHORS.txt for a list of people who contributed. +* See LICENSE_DANSE.txt for license information. +* +/************************************************************************************* +* +* class PDF3DCalculator -- Implementation of PDF3DCalculator +* +*************************************************************************************/ + +#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 PDF3DCalculator::DEFAULT_RMAX_3D; +constexpr double PDF3DCalculator::DEFAULT_GRID_STEP; + +// Constructor --------------------------------------------------------------- + +PDF3DCalculator::PDF3DCalculator() : + 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), + mtotaloccupancy(0.0), + msfaverage(0.0) +{ + this->registerDoubleAttribute("enable_nn_delta3d", + this, + &PDF3DCalculator::getEnableNNDelta3DAttr, + &PDF3DCalculator::setEnableNNDelta3DAttr); + this->registerDoubleAttribute("nn_delta3d", + this, + &PDF3DCalculator::getNNDelta3D, + &PDF3DCalculator::setNNDelta3D); + this->registerDoubleAttribute("nn_delta3d_upper_bound", + this, + &PDF3DCalculator::getNNDelta3DUpperBound); + this->registerDoubleAttribute("nn_delta_positive_eta3d", + this, + &PDF3DCalculator::getNNDeltaPositiveEta3D, + &PDF3DCalculator::setNNDeltaPositiveEta3D); + this->registerDoubleAttribute("delta1_3d", + this, + &PDF3DCalculator::getDistanceDelta1_3D, + &PDF3DCalculator::setDistanceDelta1_3D); + this->registerDoubleAttribute("delta2_3d", + this, + &PDF3DCalculator::getDistanceDelta2_3D, + &PDF3DCalculator::setDistanceDelta2_3D); + this->registerDoubleAttribute("delta_shell_index3d", + this, + &PDF3DCalculator::getDeltaShellIndex3DAttr, + &PDF3DCalculator::setDeltaShellIndex3DAttr); + this->registerDoubleAttribute("delta_shell_tolerance3d", + this, + &PDF3DCalculator::getDeltaShellTolerance3D, + &PDF3DCalculator::setDeltaShellTolerance3D); + this->registerDoubleAttribute("delta_key_tolerance3d", + this, + &PDF3DCalculator::getDeltaKeyTolerance3D, + &PDF3DCalculator::setDeltaKeyTolerance3D); + this->registerDoubleAttribute("use_adp_scale_sensitivity3d", + this, + &PDF3DCalculator::getUseADPScaleSensitivity3DAttr, + &PDF3DCalculator::setUseADPScaleSensitivity3DAttr); + this->registerDoubleAttribute("adp_scale3d", + this, + &PDF3DCalculator::getADPScale3D, + &PDF3DCalculator::setADPScale3D); + this->registerDoubleAttribute("rho0_background_scale3d", + this, + &PDF3DCalculator::getRho0BackgroundScale3D, + &PDF3DCalculator::setRho0BackgroundScale3D); + this->registerDoubleAttribute("calculation_mode3d", + this, + &PDF3DCalculator::getCalculationMode3DAttr, + &PDF3DCalculator::setCalculationMode3DAttr); + this->registerDoubleAttribute("histogram_weight_mode3d", + this, + &PDF3DCalculator::getHistogramWeightMode3DAttr, + &PDF3DCalculator::setHistogramWeightMode3DAttr); + + this->setRmax(DEFAULT_RMAX_3D); + this->setRstep(mdr); + this->setQmax(12.0); +} + +// Public Methods ------------------------------------------------------------ + +void PDF3DCalculator::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 PDF3DCalculator::getGridStep() const +{ + return mdr; +} + +void PDF3DCalculator::setAccumBlockSize(int bs) +{ + if (bs <= 0) throw std::invalid_argument("Accumulation block size must be positive."); + maccumblocksize = bs; +} + +int PDF3DCalculator::getAccumBlockSize() const +{ + return maccumblocksize; +} + +void PDF3DCalculator::setApplyRho0Background3D(bool v) +{ + mapplyrho0background3d = v; +} + +bool PDF3DCalculator::getApplyRho0Background3D() const +{ + return mapplyrho0background3d; +} + +void PDF3DCalculator::setUseCQWindow3D(bool v) +{ + musecqwindow3d = v; +} + +bool PDF3DCalculator::getUseCQWindow3D() const +{ + return musecqwindow3d; +} + +void PDF3DCalculator::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 PDF3DCalculator::getCalculationMode3D() const +{ + return mcalculationmode3d; +} + +double PDF3DCalculator::getCalculationMode3DAttr() const +{ + return static_cast(mcalculationmode3d); +} + +void PDF3DCalculator::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 PDF3DCalculator::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 PDF3DCalculator::getHistogramWeightMode3D() const +{ + return mhistogramweightmode3d; +} + +double PDF3DCalculator::getHistogramWeightMode3DAttr() const +{ + return static_cast(mhistogramweightmode3d); +} + +void PDF3DCalculator::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 PDF3DCalculator::setEnableNNDelta3D(bool v) +{ + if (menablennDelta3d != v) mticker.click(); + menablennDelta3d = v; +} + +bool PDF3DCalculator::getEnableNNDelta3D() const +{ + return menablennDelta3d; +} + +double PDF3DCalculator::getEnableNNDelta3DAttr() const +{ + return menablennDelta3d ? 1.0 : 0.0; +} + +void PDF3DCalculator::setEnableNNDelta3DAttr(double v) +{ + this->setEnableNNDelta3D(v != 0.0); +} + +void PDF3DCalculator::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& PDF3DCalculator::getNNDelta3D() const +{ + return mnnDelta3d; +} + +double PDF3DCalculator::getNNDelta3DUpperBound() const +{ + return this->currentNNDeltaUpperBound(); +} + +const double& PDF3DCalculator::getNNDeltaPositiveEta3D() const +{ + return mnnDeltaPositiveEta3d; +} + +void PDF3DCalculator::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 PDF3DCalculator::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& PDF3DCalculator::getDistanceDelta1_3D() const +{ + return mdistanceDelta1_3d; +} + +void PDF3DCalculator::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& PDF3DCalculator::getDistanceDelta2_3D() const +{ + return mdistanceDelta2_3d; +} + +void PDF3DCalculator::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& PDF3DCalculator::getDeltaPairA3D() const +{ + return mdeltaPairA3d; +} + +const std::string& PDF3DCalculator::getDeltaPairB3D() const +{ + return mdeltaPairB3d; +} + +void PDF3DCalculator::setDeltaShellIndex3D(int v) +{ + if (v <= 0) throw std::invalid_argument("Delta shell index must be positive."); + if (mdeltaShellIndex3d != v) mticker.click(); + mdeltaShellIndex3d = v; +} + +int PDF3DCalculator::getDeltaShellIndex3D() const +{ + return mdeltaShellIndex3d; +} + +double PDF3DCalculator::getDeltaShellIndex3DAttr() const +{ + return static_cast(mdeltaShellIndex3d); +} + +void PDF3DCalculator::setDeltaShellIndex3DAttr(double v) +{ + this->setDeltaShellIndex3D(static_cast(std::floor(v + 0.5))); +} + +const double& PDF3DCalculator::getDeltaShellTolerance3D() const +{ + return mdeltaShellTolerance3d; +} + +void PDF3DCalculator::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& PDF3DCalculator::getDeltaKeyTolerance3D() const +{ + return mdeltaKeyTolerance3d; +} + +void PDF3DCalculator::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 PDF3DCalculator::getDeltaEligiblePairCount3D() const +{ + return static_cast(mdeltaEligibleBondKeys.size()); +} + +void PDF3DCalculator::setUseADPScaleSensitivity3D(bool v) +{ + if (museadpscaleSensitivity3d != v) mticker.click(); + museadpscaleSensitivity3d = v; + this->validateNNDelta3D(); +} + +bool PDF3DCalculator::getUseADPScaleSensitivity3D() const +{ + return museadpscaleSensitivity3d; +} + +double PDF3DCalculator::getUseADPScaleSensitivity3DAttr() const +{ + return museadpscaleSensitivity3d ? 1.0 : 0.0; +} + +void PDF3DCalculator::setUseADPScaleSensitivity3DAttr(double v) +{ + this->setUseADPScaleSensitivity3D(v != 0.0); +} + +void PDF3DCalculator::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& PDF3DCalculator::getADPScale3D() const +{ + return madpScale3d; +} + +void PDF3DCalculator::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& PDF3DCalculator::getRho0BackgroundScale3D() const +{ + return mrho0backgroundscale3d; +} + +QuantityType PDF3DCalculator::get3DPDF() const +{ + QuantityType result; + std::vector grid = mgrid3d; + if (!this->isVectorHistogramMode3D()) + { + if (musecqwindow3d) applyQWindow3D(grid); + + const double rdf_scale = (mtotaloccupancy * msfaverage == 0.0) ? 0.0 : + 1.0 / (mtotaloccupancy * msfaverage * msfaverage); + 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)); + } + } + } + + 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 PDF3DCalculator::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 PDF3DCalculator::exportGrid3DBinary(const std::string& path, bool usefloat32, bool applypost) const +{ + std::vector grid = mgrid3d; + if (applypost) + { + if (musecqwindow3d) applyQWindow3D(grid); + + const double rdf_scale = (mtotaloccupancy * msfaverage == 0.0) ? 0.0 : + 1.0 / (mtotaloccupancy * msfaverage * msfaverage); + 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)); + } + } + } + + 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); +} + +double PDF3DCalculator::computeRho0Background() const +{ + const StructureAdapterPtr& structure = this->getStructure(); + if (!structure) return 0.0; + const double partialpdfscale = this->getPartialPDFScale(); + return mrho0backgroundscale3d * partialpdfscale * structure->numberDensity(); +} + +void PDF3DCalculator::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 PDF3DCalculator::resetValue() +{ + const StructureAdapterPtr& structure = this->getStructure(); + if (!structure) + { + msfCache.clear(); + mtotaloccupancy = 0.0; + msfaverage = 0.0; + } + else + { + int nsite = this->countSites(); + msfCache.resize(nsite); + const auto& sftable = this->getScatteringFactorTable(); + + mtotaloccupancy = structure->totalOccupancy(); + double totsf = 0.0; + for (int i = 0; i < nsite; ++i) + { + std::string atom_type = structure->siteAtomType(i); + const double sf = sftable->lookup(atom_type); + const double occ = structure->siteOccupancy(i); + const double mult = structure->siteMultiplicity(i); + msfCache[i] = sf * occ; + totsf += msfCache[i] * mult; + } + msfaverage = (mtotaloccupancy == 0.0) ? 0.0 : (totsf / mtotaloccupancy); + } + + // 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 PDF3DCalculator::addPairContribution(const BaseBondGenerator& bnds, int summationscale) +{ + if (bnds.distance() == 0.0) return; + + int i0 = bnds.site0(); + int i1 = bnds.site1(); + + if (i0 >= static_cast(msfCache.size()) || i1 >= static_cast(msfCache.size())) + return; + + const R3::Vector& rvec = bnds.r01(); + const double pairscale = bnds.multiplicity() * static_cast(summationscale); + double sfprod = msfCache[i0] * msfCache[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 PDF3DCalculator::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 PDF3DCalculator::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 * mdr + 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(pdf3d_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 PDF3DCalculator::DeltaBondKey::operator==(const DeltaBondKey& other) const +{ + return site0 == other.site0 && site1 == other.site1 && + rx == other.rx && ry == other.ry && rz == other.rz; +} + +size_t PDF3DCalculator::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 PDF3DCalculator::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); +} + +PDF3DCalculator::DeltaBondKey +PDF3DCalculator::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 PDF3DCalculator::isDeltaEligiblePair(const BaseBondGenerator& bnds) const +{ + if (!menablennDelta3d || mnnDelta3d == 0.0) return false; + return mdeltaEligibleBondKeys.find(this->makeDeltaBondKey(bnds)) != + mdeltaEligibleBondKeys.end(); +} + +double PDF3DCalculator::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 PDF3DCalculator::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 PDF3DCalculator::isDistanceDecayDeltaActive() const +{ + return menablennDelta3d && + (mdistanceDelta1_3d > 0.0 || mdistanceDelta2_3d > 0.0); +} + +double PDF3DCalculator::distanceDecayDeltaFraction(double distance) const +{ + if (distance <= 0.0) return 0.0; + return mdistanceDelta1_3d / distance + + mdistanceDelta2_3d / (distance * distance); +} + +double PDF3DCalculator::currentNNDeltaUpperBound() const +{ + if (mnnDelta3dUpperBound <= 0.0) return 0.0; + return museadpscaleSensitivity3d ? + (mnnDelta3dUpperBound * madpScale3d) : mnnDelta3dUpperBound; +} + +bool PDF3DCalculator::isVectorHistogramMode3D() const +{ + return mcalculationmode3d == 1; +} + +void PDF3DCalculator::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 PDF3DCalculator::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 PDF3DCalculator::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 PDF3DCalculator::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 PDF3DCalculator::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 PDF3DCalculator::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::PDF3DCalculator) diff --git a/src/diffpy/srreal/PDF3DCalculator.hpp b/src/diffpy/srreal/PDF3DCalculator.hpp new file mode 100644 index 00000000..c9c22515 --- /dev/null +++ b/src/diffpy/srreal/PDF3DCalculator.hpp @@ -0,0 +1,257 @@ +/************************************************************************************* +* +* libdiffpy by DANSE Diffraction group +* Simon J. L. Billinge +* (c) 2009 The Trustees of Columbia University +* in the City of New York. All rights reserved. +* +* File coded by: Hai Qiao +* +* See AUTHORS.txt for a list of people who contributed. +* See LICENSE_DANSE.txt for license information. +* +/************************************************************************************* +* +* class PDF3DCalculator -- 3D real-space pair distribution function calculator +* +*************************************************************************************/ + +#ifndef PDF3DCALCULATOR_HPP_INCLUDED +#define PDF3DCALCULATOR_HPP_INCLUDED + +#include +#include +#include +#include +#include +#include +#include +#include + +namespace diffpy { +namespace srreal { + +class PDF3DCalculator : public PDFCalculator +{ +public: + // constructor + PDF3DCalculator(); + + // Public interface to retrieve 3D PDF data + // Returns a flat vector of [x0, y0, z0, G0, x1, y1, z1, G1, ...] + QuantityType get3DPDF() 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 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; + + std::vector msfCache; + double mtotaloccupancy; + double msfaverage; + + // 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; + ar & msfCache; + ar & mtotaloccupancy; + ar & msfaverage; + } +}; + +} // namespace srreal +} // namespace diffpy + +#endif // PDF3DCALCULATOR_HPP_INCLUDED diff --git a/src/diffpy/srreal/PDFCalculator.cpp b/src/diffpy/srreal/PDFCalculator.cpp index 47c7189a..adb3cb4a 100644 --- a/src/diffpy/srreal/PDFCalculator.cpp +++ b/src/diffpy/srreal/PDFCalculator.cpp @@ -677,6 +677,14 @@ double PDFCalculator::sfAverage() const } +double PDFCalculator::getPartialPDFScale() const +{ + const double totocc = mstructure_cache.totaloccupancy; + return (totocc == 0.0) ? 0.0 : + (mstructure_cache.activeoccupancy / totocc); +} + + void PDFCalculator::cacheStructureData() { int cntsites = this->countSites(); diff --git a/src/diffpy/srreal/PDFCalculator.hpp b/src/diffpy/srreal/PDFCalculator.hpp index 79220e0a..bcce3850 100644 --- a/src/diffpy/srreal/PDFCalculator.hpp +++ b/src/diffpy/srreal/PDFCalculator.hpp @@ -114,6 +114,8 @@ class PDFCalculator : // support for PQEvaluatorOptimized virtual void stashPartialValue(); virtual void restorePartialValue(); + /// activeoccupancy / totaloccupancy used by baseline background term + double getPartialPDFScale() const; private: diff --git a/src/diffpy/srreal/R3linalg.cpp b/src/diffpy/srreal/R3linalg.cpp index 19a70cc7..337acd4e 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 bffd5c7c..91888346 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&); From 39555d4777dc71b9cea38fe185effa1fe5e04029 Mon Sep 17 00:00:00 2001 From: joeseaer <2969111832@qq.com> Date: Thu, 9 Jul 2026 01:29:43 +0800 Subject: [PATCH 4/4] MAINT: address PDF3DCalculator review feedback --- src/diffpy/srreal/PDF3DCalculator.cpp | 165 +++++++---------------- src/diffpy/srreal/PDF3DCalculator.hpp | 26 +--- src/diffpy/srreal/PDFCalculator.cpp | 19 ++- src/diffpy/srreal/PDFCalculator.hpp | 10 +- src/tests/TestPDF3DCalculator.hpp | 180 ++++++++++++++++++++++++++ src/tests/TestR3linalg.hpp | 47 +++++++ 6 files changed, 294 insertions(+), 153 deletions(-) create mode 100644 src/tests/TestPDF3DCalculator.hpp diff --git a/src/diffpy/srreal/PDF3DCalculator.cpp b/src/diffpy/srreal/PDF3DCalculator.cpp index e839ab35..e3bbca03 100644 --- a/src/diffpy/srreal/PDF3DCalculator.cpp +++ b/src/diffpy/srreal/PDF3DCalculator.cpp @@ -1,21 +1,3 @@ -/************************************************************************************* -* -* libdiffpy by DANSE Diffraction group -* Simon J. L. Billinge -* (c) 2009 The Trustees of Columbia University -* in the City of New York. All rights reserved. -* -* File coded by: Hai Qiao -* -* See AUTHORS.txt for a list of people who contributed. -* See LICENSE_DANSE.txt for license information. -* -/************************************************************************************* -* -* class PDF3DCalculator -- Implementation of PDF3DCalculator -* -*************************************************************************************/ - #include #include #include @@ -65,9 +47,7 @@ PDF3DCalculator::PDF3DCalculator() : mdeltaKeyTolerance3d(1.0e-6), museadpscaleSensitivity3d(false), madpScale3d(1.0), - mrho0backgroundscale3d(1.0), - mtotaloccupancy(0.0), - msfaverage(0.0) + mrho0backgroundscale3d(1.0) { this->registerDoubleAttribute("enable_nn_delta3d", this, @@ -444,44 +424,7 @@ QuantityType PDF3DCalculator::get3DPDF() const std::vector grid = mgrid3d; if (!this->isVectorHistogramMode3D()) { - if (musecqwindow3d) applyQWindow3D(grid); - - const double rdf_scale = (mtotaloccupancy * msfaverage == 0.0) ? 0.0 : - 1.0 / (mtotaloccupancy * msfaverage * msfaverage); - 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)); - } - } + this->applyPostProcessing3D(grid); } size_t nonzero_count = 0; @@ -523,34 +466,7 @@ void PDF3DCalculator::exportGrid3DBinary(const std::string& path, bool usefloat3 std::vector grid = mgrid3d; if (applypost) { - if (musecqwindow3d) applyQWindow3D(grid); - - const double rdf_scale = (mtotaloccupancy * msfaverage == 0.0) ? 0.0 : - 1.0 / (mtotaloccupancy * msfaverage * msfaverage); - 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)); - } - } + this->applyPostProcessing3D(grid); } std::ofstream ofs(path.c_str(), std::ios::binary | std::ios::trunc); @@ -578,6 +494,47 @@ void PDF3DCalculator::exportGrid3DBinary(const std::string& path, bool usefloat3 if (!ofs) throw std::runtime_error("Failed while writing output file: " + path); } +void PDF3DCalculator::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 PDF3DCalculator::computeRho0Background() const { const StructureAdapterPtr& structure = this->getStructure(); @@ -721,33 +678,6 @@ void PDF3DCalculator::applyQWindow3D(std::vector& grid) const void PDF3DCalculator::resetValue() { - const StructureAdapterPtr& structure = this->getStructure(); - if (!structure) - { - msfCache.clear(); - mtotaloccupancy = 0.0; - msfaverage = 0.0; - } - else - { - int nsite = this->countSites(); - msfCache.resize(nsite); - const auto& sftable = this->getScatteringFactorTable(); - - mtotaloccupancy = structure->totalOccupancy(); - double totsf = 0.0; - for (int i = 0; i < nsite; ++i) - { - std::string atom_type = structure->siteAtomType(i); - const double sf = sftable->lookup(atom_type); - const double occ = structure->siteOccupancy(i); - const double mult = structure->siteMultiplicity(i); - msfCache[i] = sf * occ; - totsf += msfCache[i] * mult; - } - msfaverage = (mtotaloccupancy == 0.0) ? 0.0 : (totsf / mtotaloccupancy); - } - // 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; @@ -771,13 +701,14 @@ void PDF3DCalculator::addPairContribution(const BaseBondGenerator& bnds, int sum int i0 = bnds.site0(); int i1 = bnds.site1(); + int cntsites = this->countSites(); - if (i0 >= static_cast(msfCache.size()) || i1 >= static_cast(msfCache.size())) + if (i0 >= cntsites || i1 >= cntsites) return; const R3::Vector& rvec = bnds.r01(); const double pairscale = bnds.multiplicity() * static_cast(summationscale); - double sfprod = msfCache[i0] * msfCache[i1] * pairscale; + double sfprod = this->sfSite(i0) * this->sfSite(i1) * pairscale; if (this->isVectorHistogramMode3D()) { @@ -860,7 +791,7 @@ void PDF3DCalculator::addAnisotropicGaussianToGrid(const R3::Vector& r_ij, const double min_val = center_val - cutoff_radius; double max_val = center_val + cutoff_radius; - // Node-centered grid: points at k * mdr + // 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)); diff --git a/src/diffpy/srreal/PDF3DCalculator.hpp b/src/diffpy/srreal/PDF3DCalculator.hpp index c9c22515..b14abc97 100644 --- a/src/diffpy/srreal/PDF3DCalculator.hpp +++ b/src/diffpy/srreal/PDF3DCalculator.hpp @@ -1,21 +1,3 @@ -/************************************************************************************* -* -* libdiffpy by DANSE Diffraction group -* Simon J. L. Billinge -* (c) 2009 The Trustees of Columbia University -* in the City of New York. All rights reserved. -* -* File coded by: Hai Qiao -* -* See AUTHORS.txt for a list of people who contributed. -* See LICENSE_DANSE.txt for license information. -* -/************************************************************************************* -* -* class PDF3DCalculator -- 3D real-space pair distribution function calculator -* -*************************************************************************************/ - #ifndef PDF3DCALCULATOR_HPP_INCLUDED #define PDF3DCALCULATOR_HPP_INCLUDED @@ -153,6 +135,7 @@ class PDF3DCalculator : public PDFCalculator 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(); @@ -207,10 +190,6 @@ class PDF3DCalculator : public PDFCalculator std::unordered_set mdeltaEligibleBondKeys; std::vector mdeltaShellRecords; - std::vector msfCache; - double mtotaloccupancy; - double msfaverage; - // Internal constants static constexpr double DEFAULT_RMAX_3D = 10.0; // Angstrom static constexpr double DEFAULT_GRID_STEP = 0.1; // Angstrom @@ -245,9 +224,6 @@ class PDF3DCalculator : public PDFCalculator ar & museadpscaleSensitivity3d; ar & madpScale3d; ar & mrho0backgroundscale3d; - ar & msfCache; - ar & mtotaloccupancy; - ar & msfaverage; } }; diff --git a/src/diffpy/srreal/PDFCalculator.cpp b/src/diffpy/srreal/PDFCalculator.cpp index adb3cb4a..e2df550f 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); @@ -685,6 +680,16 @@ double PDFCalculator::getPartialPDFScale() const } +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 bcce3850..5d5eba4b 100644 --- a/src/diffpy/srreal/PDFCalculator.hpp +++ b/src/diffpy/srreal/PDFCalculator.hpp @@ -116,6 +116,12 @@ class PDFCalculator : 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: @@ -147,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/tests/TestPDF3DCalculator.hpp b/src/tests/TestPDF3DCalculator.hpp new file mode 100644 index 00000000..dfc58d7e --- /dev/null +++ b/src/tests/TestPDF3DCalculator.hpp @@ -0,0 +1,180 @@ +/***************************************************************************** +* +* class TestPDF3DCalculator -- unit tests for PDF3DCalculator class +* +*****************************************************************************/ + +#include + +#include +#include +#include +#include + +#include +#include +#include + +using namespace std; +using namespace diffpy::srreal; + +class TestPDF3DCalculator : public CxxTest::TestSuite +{ + private: + + boost::shared_ptr mpdfc; + AtomicStructureAdapterPtr mstru2; + double meps; + + public: + + void setUp() + { + meps = diffpy::mathutils::SQRT_DOUBLE_EPS; + mpdfc.reset(new PDF3DCalculator); + + mstru2 = boost::make_shared(); + Atom ai; + ai.atomtype = "Ni"; + ai.xyz_cartn[0] = 0.0; + ai.xyz_cartn[1] = 0.0; + ai.xyz_cartn[2] = 0.0; + ai.uij_cartn = R3::identity(); + ai.uij_cartn(0, 0) = ai.uij_cartn(1, 1) = + ai.uij_cartn(2, 2) = 0.004; + mstru2->append(ai); + ai.xyz_cartn[0] = 1.0; + mstru2->append(ai); + } + + + 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->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); + TS_ASSERT_EQUALS(1, mpdfc->getCalculationMode3D()); + TS_ASSERT_EQUALS(1, mpdfc->getHistogramWeightMode3D()); + TS_ASSERT_EQUALS(true, mpdfc->getEnableNNDelta3D()); + } + + + void test_invalid_configuration() + { + TS_ASSERT_THROWS(mpdfc->setGridStep(0.0), invalid_argument); + TS_ASSERT_THROWS(mpdfc->setAccumBlockSize(0), invalid_argument); + TS_ASSERT_THROWS(mpdfc->setCalculationMode3D(2), invalid_argument); + TS_ASSERT_THROWS(mpdfc->setHistogramWeightMode3D(2), + invalid_argument); + TS_ASSERT_THROWS( + mpdfc->setDoubleAttr("calculation_mode3d", 0.5), + invalid_argument); + TS_ASSERT_THROWS(mpdfc->setADPScale3D(0.0), invalid_argument); + TS_ASSERT_THROWS(mpdfc->setRho0BackgroundScale3D(-1.0), + invalid_argument); + } + + + void test_vector_histogram_eval() + { + mpdfc->setRmax(2.0); + mpdfc->setGridStep(0.5); + mpdfc->setCalculationMode3D(1); + mpdfc->setHistogramWeightMode3D(1); + mpdfc->eval(mstru2); + + QuantityType pdf3d = mpdfc->get3DPDF(); + TS_ASSERT(!pdf3d.empty()); + TS_ASSERT_EQUALS(0u, pdf3d.size() % 4); + + QuantityType rh = mpdfc->getRadialHistogram3D(); + size_t nr = static_cast( + std::ceil(mpdfc->getRmax() / mpdfc->getGridStep())) + 1; + TS_ASSERT_EQUALS(2 * nr, rh.size()); + + bool found_distance_one = false; + for (size_t i = 0; i < rh.size(); i += 2) + { + if (std::fabs(rh[i] - 1.0) < meps && rh[i + 1] != 0.0) + { + found_distance_one = true; + break; + } + } + TS_ASSERT(found_distance_one); + } + + + void test_serialization() + { + 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); + boost::shared_ptr pdfc1; + ia >> pdfc1; + + TS_ASSERT_DIFFERS(pdfc1.get(), mpdfc.get()); + TS_ASSERT_DELTA(0.25, pdfc1->getGridStep(), meps); + TS_ASSERT_EQUALS(7, pdfc1->getAccumBlockSize()); + TS_ASSERT_EQUALS(false, pdfc1->getApplyRho0Background3D()); + TS_ASSERT_EQUALS(false, pdfc1->getUseCQWindow3D()); + TS_ASSERT_EQUALS(1, pdfc1->getCalculationMode3D()); + TS_ASSERT_EQUALS(1, pdfc1->getHistogramWeightMode3D()); + TS_ASSERT_EQUALS(true, pdfc1->getEnableNNDelta3D()); + TS_ASSERT_DELTA(0.02, pdfc1->getNNDelta3D(), meps); + TS_ASSERT_DELTA(0.8, pdfc1->getNNDeltaPositiveEta3D(), meps); + TS_ASSERT_DELTA(0.03, pdfc1->getDistanceDelta1_3D(), meps); + TS_ASSERT_DELTA(0.04, pdfc1->getDistanceDelta2_3D(), meps); + TS_ASSERT_EQUALS(string("Ni"), pdfc1->getDeltaPairA3D()); + TS_ASSERT_EQUALS(string("Ni"), pdfc1->getDeltaPairB3D()); + TS_ASSERT_EQUALS(2, pdfc1->getDeltaShellIndex3D()); + TS_ASSERT_DELTA(0.06, pdfc1->getDeltaShellTolerance3D(), meps); + TS_ASSERT_DELTA(1.0e-5, pdfc1->getDeltaKeyTolerance3D(), meps); + TS_ASSERT_EQUALS(true, pdfc1->getUseADPScaleSensitivity3D()); + TS_ASSERT_DELTA(1.2, pdfc1->getADPScale3D(), meps); + TS_ASSERT_DELTA(0.5, pdfc1->getRho0BackgroundScale3D(), meps); + } + +}; // class TestPDF3DCalculator + +// End of file diff --git a/src/tests/TestR3linalg.hpp b/src/tests/TestR3linalg.hpp index 8f779dfb..e84ee175 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