Skip to content

Adding truncated tetrahedron model pure python#718

Open
sara-mokhtari wants to merge 9 commits intomasterfrom
adding_truncated_tetrahedron_model_purePython
Open

Adding truncated tetrahedron model pure python#718
sara-mokhtari wants to merge 9 commits intomasterfrom
adding_truncated_tetrahedron_model_purePython

Conversation

@sara-mokhtari
Copy link
Copy Markdown
Contributor

Description

Adding pure python models for regular tetrahedrons (tetrahedron.py) and truncated tetrahedrons (tetrahedron_truncated.py) with orientation averaging using Fibonacci quadrature.

Fixes #713

How Has This Been Tested?

Unit tests worked.
Both of the models were compared to an another numerical calculation : for different sizes and truncatures t=0 (regular tetrahedron), t=0.33 (friauf polyhedron), t=0.5 (octahedron), good agreement was found with the Debye equation (Debye calculator: https://github.com/FrederikLizakJohansen/DebyeCalculator). The Debye equation is based on atomic positions while SasView model is based on analytical expressions. Example of comparison on a regular tetrahedron (no truncature):

regularTd_sasview_vs_debye

Since t=0.5 leads to an octahedron, the model was also compared with the previous octahedron_truncated model and showed good agreement.

Discussion points:

  • Integration: for now, the models use the Fibonacci quadrature, but they will be adapted with the most accurate integration method chosen (same issue that the prism model, see adaptive integration for the prism model #716 )
  • Singularities: the general form factor of the tetrahedron and the treatment of singularities can be found in the Fq_ortho() function. This function was inspired from the orthogonal_tetrahedra_formfactor() function in https://github.com/ch-tung/simplex_scattering/blob/main/form_factor.py. I refactored the function to improve numerical robustness in singular configurations. Previously, formulas were evaluated more directly, which could trigger division-by-zero warnings even when the final value was later corrected by singular case logic. The updated implementation uses boolean masks for each case (no singularties case, singularities on vertices, singularities on edges, singularities on faces and q getting close to 0) and evaluates each analytical expression only on its valid domain. It seems to be removing the warnings, let me know if it is the right correction, also some other optimizations could probably be done.
  • The truncated_tetrahedron model allows to create a full tetrahedron too, let me know if having the two models is considered redundant.

Review Checklist:

[if using the editor, use [x] in place of [ ] to check a box]

Documentation (check at least one)

  • There is nothing that needs documenting
  • Documentation changes are in this PR
  • There is an issue open for the documentation (link?)

Installers

  • There is a chance this will affect the installers, if so
    • Windows installer (GH artifact) has been tested (installed and worked)
    • MacOSX installer (GH artifact) has been tested (installed and worked)
    • Wheels installer (GH artifact) has been tested (installed and worked)

Licensing (untick if necessary)

  • The introduced changes comply with SasView license (BSD 3-Clause)

Copy link
Copy Markdown
Contributor

@pkienzle pkienzle left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I verified that the truncated tetrahedron model matches a regular octahedron with the same length side using

$ python -m sasmodels.compare -pars octahedron_truncated,tetrahedron_truncated background=0 t=1.0,0.5 length_a=100 "circumradius=sqrt(3)*length_a" -double -midq

The scale factor on length_a comes from conversion between length_a and edge length of the octahedron and from the conversion of twice edge length to circumradius for the tetrahedron.


.. math::

F(\mathbf{q}) = \int_V \exp(i\mathbf{q} \cdot \mathbf{r}) \, d\mathbf{r}
Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Fourier transform uses $e^{-iq \cdot r}$

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This is a good point !
Usually, I'm using the minus sign convention for Fourier transform.
But here, I decided on purpose to shift to the plus sign convention instead.
In order to be consistent with other publications.
The reference by J. Wuttke and the code for simplex decomposition written by Chi-Huan Tung which we are using for polyhedron models:
https://github.com/ch-tung/simplex_scattering

Both sign conventions are used in the SAS community, I had a look at different authors, it is about 50/50 !

Anyway, it is easy to shift from one convention to the other by replacing $\vec{q}$ by - $\vec{q}$.

The choice needs to be consistent between different models in Fqabc() computations.
And it will be important only when combining different form factor amplitudes together.
In a core/shell polyhedral particle or for a small cluster of rods for example.

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I put in a ticket #721 to check the consistency. We should at least document the convention used per model even if we don't make them consistent.

"tetrahedron scattering length density"],
["sld_solvent", "1e-6/Ang^2", 9.4, [-inf, inf], "sld",
"Solvent scattering length density"],
["circumradius", "Ang", 100, [0., inf], "volume",
Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I think "radius" is sufficiently clear, especially if the description says that it is the circumradius.

from sasmodels.special.fibonacci import fibonacci_sphere

name = "tetrahedron_truncated"
title = "Tetrahedron truncated"
Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Following the pattern of "elliptical_cylinder" and "rectangular_prism", the truncated models should be "truncated_tetrahedron" and "truncated_octahedron".

"Solvent scattering length density"],
["circumradius", "Ang", 100, [0., inf], "volume",
"Circumradius of the full tetrahedron"],
["t", "", 0, [0, 0.5], "volume",
Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Sasmodels prefers longer names, like "truncation".

Note that this definition of t is 1-t from the octahedron model. The behaviour for parameters of the same name should be similar between the models. At a minimum, the parameter description should say that zero means no truncation.

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

OK, my bad !
I have changed the definition of t, replacing it by (1-t).
I agree that it should be consistent in all truncated models.

We probably need to fix it in the octahedron model, but I need to check the figure again ...

Convention for no truncation can be either t=0 or t=1

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I find truncation in [0, 0.5] is more natural; it is the amount of truncation. I got it wrong when I first played with the octahedron model.

You could change t to a name that suggests the portion of the model to included, with portion in [0.5, 1.0]. Using truncation is better, though, since it is a truncated tetrahedron.

If octahedron is already available in a released version of sasview then we will need to add translation code in sasmodels/convert.py so that saved models can be reloaded. [The conversion code should be stored with the model, but #249 has not yet been addressed.]

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I start today to update the truncated_octahedron model (.py and .c) but from this PR (related to adaptative integration): #710
branch name is adaptive-octahedron
Please follow on this other branch ...

This in order to have the same truncature convention for all models:
t=0 for no truncation of the shape ! :)

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

None yet

Projects

None yet

Development

Successfully merging this pull request may close these issues.

Regular and truncated tetrahedron models

3 participants