Add tetrapod model implementation#705
Conversation
Co-authored-by: IndigoCarmine <79552713+IndigoCarmine@users.noreply.github.com>
…-functionality Copilot checked it using formula of the paper.
… author attribution in tetrapod.py
|
@IndigoCarmine. I notice you have left this as a draft request (not ready for review) with the last commit in late March. Just checking to see if you are waiting for or needing something? No worries if you just ran out of time. we all do. Just want to make sure you aren't abandoned. |
|
Sorry for the long silence. My workload has increased recently, and I haven't had much time to work on this. Most of the implementation is already done, but progress has been much slower than I had hoped. I apologize for the delay. I have not abandoned this work and do intend to finish it. However, I cannot promise a timeline at the moment. If someone else is interested in working on it in the meantime, I would be happy for them to build on my current work or even start from scratch. Thank you for checking in and for your patience. |
|
I have finished implementing the standard model (not the core-shell version). The calculated values appear to be consistent with those reported in the reference paper. To be honest, I am not entirely confident about the scaling, so there may still be some issues. |
|
Implemented an initial tetrapod core-shell model. The scattering calculation, documentation, and example results are available and produce reasonable output. The following features are not yet implemented:
|
pkienzle
left a comment
There was a problem hiding this comment.
Code looks okay. Math and plots match those in the paper.
$ PYTHONPATH=. python -m sasmodels.compare tetrapod -q=1e-3:10 -nq=400 background=0 radius=5,10 thickness=0 length=100 scale=100,25
$ PYTHONPATH=. python -m sasmodels.compare tetrapod -q=1e-3:10 -nq=400 background=0 radius=5 length=10,100 scale=100,10
| References | ||
| ---------- | ||
|
|
||
| #. Seoki Kyoo Seo *Korean J. Chem. Eng.* 34(2017) 1192-1198 |
There was a problem hiding this comment.
DOI:10.1007/s11814-016-0341-x
|
|
||
| static double form_volume(double L, double R, double t) { | ||
| // V = 4 * pi * (R + t)^2 * L (outer volume of 4 arms) | ||
| return 4.0 * M_PI * (R + t) * (R + t) * L; |
There was a problem hiding this comment.
An overestimate of the volume since the ends of the cylinders overlap at the origin. The effect will be significant when the arms are short and wide.
There was a problem hiding this comment.
Thank you for your reply. I agree with your point. However, the scattering model also ignores overlap effects, and any resulting artifacts appear only in the high-(q) region (provably?). Therefore, for now, I have followed the treatment used in the literature.
If this issue remains a concern, I suggest that a model such as a core–shell sphere with four arms would be more appropriate for addressing it. What do you think?
There was a problem hiding this comment.
Agreed that you don't need to do anything. The model corresponds to the published model.
You may want to note in the documentation that the approximation is only valid for length >> radius+thickness.
I'm a little concerned about the effects of volume normalization on low q values, but it may be that scattering is overestimated by as much as volume, and so errors cancel. We could use Monte Carlo sampling and the generic scattering calculator to confirm this.
| # [ "name", "units", default, [lower, upper], "type", "description"], | ||
| parameters = [ | ||
| ["length", "Ang", 1000, [0, inf], "volume", "Arm length"], | ||
| ["radius", "Ang", 50, [0, inf], "volume", "Arm outer radius"], |
There was a problem hiding this comment.
In core-shell models this is the core radius rather than the outer radius. See "thickness" in #305
| double contrast_shell = sld_shell - sld_solvent; | ||
| double total = 0.0; | ||
|
|
||
| for (int dtheta = 0; dtheta < GAUSS_N; dtheta++) { |
There was a problem hiding this comment.
Doesn't use adaptive integration.
There was a problem hiding this comment.
Thank you for the suggestion.
I tried switching from gauss76.c to gauss150.c for the orientation integration, but the fitting result did not change significantly (as shown in the figure above).
I understand that your comment refers to a more genuine adaptive integration scheme.
Could you clarify what approach you had in mind? For example, are you suggesting a specific integration scheme already used in other sasmodels models?
There was a problem hiding this comment.
Target accuracy for the models is 5 digits of precision for SANS (R=20 nm at q=1/Å) and USANS (R=20 μm at q=0.002/Å). For USANS slit smearing, which pushes q out to 0.1/Å, the target accuracy 0.2. With regard to performance, I'm trying to limit it to 2 s to evaluate 200 q points on my mac M2 chip.
I recently revised all the models to use a simple adaptive integration scheme (#658) based on qr max for the shape along the c-axis and in the a-b plane. For the outer loop I use the max of these two, but only the a-b cross section for the inner loop. I'm limiting the (θ,φ) grid to 100 000 evaluation points, with no more than 76 points in the outer loop.
Here's an example from triaxial ellipsoid:
const double qr_max_inner = fmax(q*radius_equat_minor, q*radius_equat_major);
const double qr_max = fmax(qr_max_inner, q*radius_polar);
constant double *z_outer, *w_outer;
constant double *z_inner, *w_inner;
int n_outer = gauss_weights(qr_max, ADAPTIVE_MAX_OUTER, &w_outer, &z_outer);
int n_inner = gauss_weights(qr_max_inner, n_outer, &w_inner, &z_inner);For the tetrapod, qr_max = q*sqrt(length**2 + (outer radius)**2) is probably all you need. Whether outer radius is radius or core radius + thickness depends on your parameterization.
Code in explore/check_adaptive.py compares the adaptive grid to a 5000x5000 grid on a few points and reports those that are out of tolerance.
This is not needed for this PR, but please convert this comment into a new ticket if you choose to do it later. [It is now in #705]
There was a problem hiding this comment.
@pkienzle is that change to adaptive integration in the master branch yet? or is it only in the release_1.1.0 branch?
| @@ -0,0 +1,75 @@ | |||
| // acos(-1.0/3.0)/2.0 = half of tetrahedral angle ~54.7356 deg | |||
| #define TETRAHEDRAL_HALF_ANGLE acos(-1.0 / 3.0) / 2.0 | |||
There was a problem hiding this comment.
The following can be precomputed:
cos_A = cos(TETRAHEDRAL_HALF_ANGLE) = √(1/3)
sin_A = sin(TETRAHEDRAL_HALF_ANGLE) = √(2/3)
| parameters = [ | ||
| ["length", "Ang", 1000, [0, inf], "volume", "Arm length"], | ||
| ["radius", "Ang", 50, [0, inf], "volume", "Arm outer radius"], | ||
| ["core_radius", "Ang", 50, [0, inf], "volume", "Arm core radius"], |
There was a problem hiding this comment.
You are still subtracting thickness from radius internally, so this parameter is outer radius, not the core radius. To avoid confusion with core-shell models, I suggest changing the name to radius_outer or leaving it as radius and have thickness add to the radius rather than subtract from the radius.
We are following the "mathematical" convention radius_tag rather than the natural language tag_radius for most parameter names in our models.
This pull request introduces a new tetrapod-shaped particle model.
#641
The current implementation and accompanying documentation include some AI-generated context that has not yet been thoroughly reviewed or validated. Sorry. I will revise it soon.
The issue originally proposed a CoreShell model, but for the moment I will implement the model as described in the referenced paper.
Reference:
https://github.com/user-attachments/files/19167378/s11814-016-0341-x.pdf