Skip to content

Commit

Permalink
some updates to lamellar_x5 still sorting normalisations
Browse files Browse the repository at this point in the history
  • Loading branch information
@RichardHeenan
RichardHeenan committed Oct 13, 2023
1 parent 2a17ea5 commit d84ebe1
Show file tree
Hide file tree
Showing 2 changed files with 106 additions and 56 deletions.
84 changes: 58 additions & 26 deletions sasmodels/models/lamellar_x5_paracrystal_kr.c
View file
Original file line number Diff line number Diff line change
@@ -1,6 +1,15 @@
/*Lamellar_ParaCrystal - July 2019 RKH - hack previous sasview version to become Kotlarchyk & Ritzau
J.Appl.Cryst. 24(1991)753-758 with some "corrections" as per version from FISH .
SHOULD put paraCryst_kr and rlayff_kr into a library as they are also used by the original lamellar_stack_paracrystal_kr.c
/*Lamellar_ParaCrystal
July 2019 RKH - converted his FISH Fortran code for the Kotlarchyk & Ritzau form factor of sheets and structure factor
for 1d paracrystals into c++ for sasview. Ref. J.Appl.Cryst. 24(1991)753-758 with some "corrections" as per version from FISH .
SHOULD put paraCryst_kr and rlayff_kr into a library as they could also be used by other routines
12 June 2020 initial version of lamellar12345 type models, (released privately to two user groups).
Sept 2023 this x5 version has explicit amounts of 1 through 5 layers, plus a stack with an adjustable number of layers
BEWARE P(Q) equations divide by zero if sigma_t_by_t =0.0 (and may not then give results even when reset!)
This is very flexible, particularly for mixtures of unilamellar and multilamellar particles.
IF sigma_d_by_d = -1 ignore the S(Q), return only P(Q) for sheet, including (scale1 +scale2+ ...)(del rho)^2
IF sigma_t_by_t = -1 ignore the P(Q), return only AVERAGE S(Q) sum(scale_n*S(Q,n) )/sum(scale_n), or of course a
single S(Q,n) if the other scale_n are set to zero.
*/
double paraCryst_kr(double qval, double RM, double D, double GD);
double rlayff_kr(double qval, double RLBAR, double GL, double TT, double RSIG, double f1);
Expand All @@ -16,16 +25,19 @@ Iq(double qval,
double scale3,
double scale4,
double scale5,
double scaleMadj,
double Madj,
double davg,
double sigma_d_by_d,
double sld,
double solvent_sld)
{
double f1, Znq2,Znq3,Znq4,Znq5;
double f1, f2, Znq2, Znq3, Znq4, Znq5, ZnqM, sqsum, contr;
const double fp_2 = 2.0;
const double fp_3 = 3.0;
const double fp_4 = 4.0;
const double fp_5 = 5.0;
sqsum = 0.0;

/* //get the fractional part of Nlayers, to determine the "mixing" of N's
// have removed this averaging of two values of Nlayers, as the pattern changes
Expand All @@ -42,27 +54,41 @@ Iq(double qval,
//calculate the n2 contribution
Znq += (1.0-xn)*paraCryst_kr(qval,(double)n2,davg,sigma_d_by_d);
*/

Znq2 = paraCryst_kr(qval,fp_2,davg,sigma_d_by_d);
Znq3 = paraCryst_kr(qval,fp_3,davg,sigma_d_by_d);
Znq4 = paraCryst_kr(qval,fp_4,davg,sigma_d_by_d);
Znq5 = paraCryst_kr(qval,fp_5,davg,sigma_d_by_d);
//
// will one day separate the S(Q) and P(Q) here into individual models
// all these IF's will I'm told will reduce the efficiency of gpu code
if (sigma_d_by_d >= -0.999){
// normal route
if(scale2 > 0.0){sqsum += scale2*paraCryst_kr(qval,fp_2,davg,sigma_d_by_d);}
if(scale3 > 0.0){sqsum += scale3*paraCryst_kr(qval,fp_3,davg,sigma_d_by_d);}
if(scale4 > 0.0){sqsum += scale4*paraCryst_kr(qval,fp_4,davg,sigma_d_by_d);}
if(scale5 > 0.0){sqsum += scale5*paraCryst_kr(qval,fp_5,davg,sigma_d_by_d);}
if(scaleMadj > 0.0){sqsum += scaleMadj*paraCryst_kr(qval,Madj,davg,sigma_d_by_d);}
} else {
// return only scaled P(Q)
sqsum = scale2 +scale3 +scale4 +scale5 +scaleMadj;
}
// Could likely speed up a little with a bespoke paraCryst_kr for all five at once as parts of the
// calculation are repeated.
// One day separate the S(Q) and P(Q) here into individual models
// meanwhile it should be relatively simple to make other versions with shell/core/shell layers etc.
// I believe that <F> and F^2> are being properly averaged for the polydisperse layer
const double f2 = rlayff_kr(qval,thickness,sigma_t_by_t,interface_t,rsig_lorentz,f1);

if (sigma_t_by_t < -0.999){
// return only average S(Q) which should tend to 1.0 at high Q
return sqsum/(scale2 +scale3 +scale4 +scale5 +scaleMadj);
} else if( sigma_t_by_t > 1.0e-08){
contr = sld - solvent_sld;
f2 = 0.25*contr*contr*rlayff_kr(qval,thickness,sigma_t_by_t,interface_t,rsig_lorentz,f1);
} else{
// catch divide by zero when sigma_t_by_t = 0.0
// TODO: add proper P(Q) for monodisperse sheet
f2 =0.0;
}
//
const double contr = sld - solvent_sld;
// BEWARE may need to rescale f1 also for sasview 5
// do not know yet why need a rescale of 0.25 here to get same scale parameter as FISH
// need check against lamellar model!
// 12June2020 add 5 terms
const double inten = 0.25*contr*contr*(scale1 + scale2*Znq2 + scale3*Znq3 + scale4*Znq4 + scale5*Znq5)*f2;
return ( scale1 + sqsum )*f2*1.0e-04;

}

// functions for the lamellar paracrystal model, hacked from FISH fortran code, July 2019, RKH
// functions for the lamellar paracrystal model, converted from my own FISH fortran code, July 2019, RKH
double
paraCryst_kr(double QQ, double RM, double D, double GD) {
// presume RMU is mu=1 = cos(angle between Q and axis of stack) in eqn (13)
Expand Down Expand Up @@ -109,6 +135,8 @@ rlayff_kr(double QQ, double RLBAR, double GL, double TT, double RSIG, double f1)
// GL = SIGMA(L)/LBAR = (Z+1)**-0.5
// 16/3/93 add extra parameter RSIG to model, as Lorentz factor,
// see Skipper et.al. J.Chem.Phys 94(1991)5751-5760
// Sept 2023, divide by RLBAR.RSIG^2 to try to get volume fraction phi scaling correct
// should also reduce strong parameter correlation between these and scale parameters.
double ZP1 = 1.0/square(GL);
const double Z = ZP1 - 1.0;
double ZM1 = Z - 1.0;
Expand All @@ -125,20 +153,24 @@ rlayff_kr(double QQ, double RLBAR, double GL, double TT, double RSIG, double f1)
// I believe that <F> and F^2> are being properly averaged for the polydisperse layer
double AL2SQ = square(AL2);
// f1 = (1. + 1.0/AL2SQ)**(-0.5*ZP1)*sqrt(AL2SQ + 1)*sin(Z*atan(1.0/AL2))*exp(0.5*Q2S2)/Z;
f1 = pow( 1.0 + 1.0/AL2SQ, -0.5*ZP1);
f1 *= sqrt(AL2SQ + 1)*sin(Z*atan(1.0/AL2))*exp(0.5*Q2S2)/Z;
f1 /= sqrt(DEN);

//
f1 = 0.0;
// start of final code for f1, comment out here as we are not using it anywhere
//f1 = pow( 1.0 + 1.0/AL2SQ, -0.5*ZP1);
//f1 *= sqrt(AL2SQ + 1)*sin(Z*atan(1.0/AL2))*exp(0.5*Q2S2)/Z;
//f1 /= sqrt(DEN);
// end of code for f1

// original fortran
// RLAYFF= (AL2**ZP1)*( (AL2**(1.0-Z)) -((AL2**2 +4.0)**(-0.5*ZM1))*
// > COS( ZM1*atan(1.0/AL) ) )*EXP(Q2S2)/ (2.0*Z*ZM1*QQ*QQ)

//rearrange to avoid overflows ( note AL can be 100, Z say 25 )
//included an extra factor of 4.0 to get beta(Q) to go to 1.0 at low
// Q and P(Q) to go to 1.0/Q**2
//included an extra factor of 4.0 to get beta(Q) to go to 1.0 at low Q and P(Q) to go to 1.0/Q**2
// f2 = 4.0*(AL2SQ)*(1. - ((1. + 4.0/(AL2SQ))**(-0.5*ZM1))*cos(ZM1*atan(1.0/AL)))*exp(Q2S2)/(2.0*Z*ZM1);

double f2 = 4.0*(AL2SQ)*(1. - pow(1. + 4.0/AL2SQ,-0.5*ZM1)*cos(ZM1*atan(1.0/AL))) * exp(Q2S2)/(2.0*Z*ZM1);
f2 /= DEN;
return f2;
return 1.0e-04*(1.0 + square(RSIG))*RLBAR*f2;
}

78 changes: 48 additions & 30 deletions sasmodels/models/lamellar_x5_paracrystal_kr.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -9,12 +9,13 @@
allow for some curvature or waviness of the layers.
This x5 version has explicit amounts of N = 1,2,3,4 or 5 layers,
so "scale" is complicated to interpret in terms of volume
fractions. Thus a distribution in N may be designed, though beware
that the best fit may be poorly determined, so do not adjust too
many parameters at once! In practice this model has proven to be
very effective in systems with a mixture of
unilammellar and multilammellar vesicles with a few layers.
plus an additional stack with an adjustable number of M\_adj layers.
Thus a distribution in N may be designed and/or the number of layers
fit. Beware that the best fit may be poorly determined, so do not
adjust too many parameters at once! In practice this model has proven
to be very effective in systems with a mixture of
unilammellar and multilammellar vesicles.
The vesicles must be sufficiently large (or flexible) that the
interference in Q space from opposite faces, across their diameter, is
Expand All @@ -26,12 +27,16 @@
This is a first attempt at the equations, reports of any issues
with them or the calculations in the code will be gratefully received.
Further explanations of the scale parameters will follow. It will
be best to leave the overall scale parameter included by sasview set
to 1.0 or some other convenient value.
Definition
----------
In the equations below,
- *scale* for each of $Ni = 1, 2, 3, 4, 5$ layer stacks determines
- *scale* for each of $Ni = 1, 2, 3, 4, 5$ and then $Madj$ layer stacks determines
the amount of each in the final intensity,
- *sld* $-$ *sld_solvent* is the contrast $\Delta \rho$,
Expand All @@ -44,10 +49,11 @@
- *interface_t* is the width of the interfaces, note $\sigma_{interface} = interface\_t/\sqrt{2\pi}$,
- *The number of layers* in each stack is $Ni = 1, 2, 3, 4, 5$,
- *The number of layers* in each stack is $Ni = 1, 2, 3, 4, 5$, then $Madj >= 2$,
- *d_spacing* is the average distance between adjacent layers
$d = \langle D \rangle$,
$d = \langle D \rangle$, its initial value will need to be close to $2\pi/Qpeak$
for a fit to work,
- *sigma_d* is the standard deviation of the Gaussian
layer distance distribution from $\sigma_D / \langle D \rangle$, and
Expand All @@ -63,17 +69,21 @@
.. math::
I(q) = \frac{\Delta\rho^2P_{layer}(q)}{4(1+\text{rsig_lorentz}^2q^2/2)} \sum_{i=1}^{5} \left[ scale_i.S(q,i)\right]
I(q) = \frac{\Delta\rho^2P_{layer}(q)}{4(1+rsig{\_}lorentz^2q^2/2)}( scale\_N1 + \sum_{i=2}^{5} \left[ scale\_Ni.S(q,i)\right] )
The form factor of the layer follows eqn (17) in the reference (but has been
re-arranged in the c++ code to help avoid computational errors). There is some
uncertainty regarding the leading factor of 4, which now seems to cancel with
the 4 in the numerator of the previous equation, see notes in c++ code.
This will then be multiplied by the leading sasview scale parameter and a
flat background added. The form factor of the layer, $P\_{layer}(q)$ follows
eqn (17) in the reference (but has been re-arranged in the c++ code to help
avoid computational errors).
There is some uncertainty regarding the leading factor of 4, which now seems
to cancel with the 4 in the numerator of the previous equation, see notes
in c++ code.
.. math::
P(q) = \frac{4(2\alpha^{z+1}(2\alpha^{(1-z)} -
P\_{layer}(q) = \frac{4(2\alpha^{z+1}(2\alpha^{(1-z)} -
((2\alpha)^2+4)^{-(z-1)/2})}{2z(z-1)}
\cos((z-1)tan^{-1}(1/\alpha))\exp{(-q^2\sigma_{interface}^2/2)}
Expand Down Expand Up @@ -101,15 +111,15 @@
.. math::
I_c(q,N) = -2F((1+F^2)\cos(qd)-2F-F^N\cos((N+1)qd)+
2F^{(N+1)}\cos(Nqd)-F^{(N+2)}\cos((N-1)qd))/(1-2F\cos(qd)+F^2)^2
I_c(q,N) = -2F&\biggl\{(1+F^2)\cos(qd)-2F-F^N\cos((N+1)qd) \\
& +2F^{(N+1)}\cos(Nqd)-F^{(N+2)}\cos((N-1)qd) \biggr\}/(1-2F\cos(qd)+F^2)^2
where
.. math::
Z(q,N) = \frac{(1-F^2)}{(1-2F\cos(qd)+F^2)}
Z(q) = \frac{(1-F^2)}{(1-2F\cos(qd)+F^2)}
and
Expand Down Expand Up @@ -143,8 +153,8 @@
name = "lamellar_x5_paracrystal_kr"
title = "Random lamellar sheet plus stacks with paracrystal structure factors"
description = """\
[Random lamellar sheet plus stacks with paracrystal structure factors]
work in progress
[lamellar sheet plus stacks with paracrystal structure factors]
"""
category = "shape:lamellae"
have_Fq = False
Expand All @@ -169,6 +179,10 @@
"relative amount of N = 4, four layers"],
["scale_N5", "", 0.0, [0.0, inf], "",
"relative amount of N = 5, five layers"],
["scale_Madj", "", 0.0, [0.0, inf], "",
"relative amount of Madj layers"],
["Madj", "", 2.0, [2.0, 50], "",
"adjustable number of layers"],
["d_spacing", "Ang", 75., [0.0, inf], "",
"lamellar spacing of paracrystal stack"],
["sigma_d_by_d", "", 0.05, [0.0, inf], "",
Expand All @@ -179,52 +193,56 @@
"Solvent scattering length density"]
]


source = ["lamellar_x5_paracrystal_kr.c"]

form_volume = """
return 1.0;
"""

#TODO test this random function
def random():
total_thickness = 10**np.random.uniform(2, 4.7)
scale_N1 = np.random.randint(0, 1)
scale_N2 = np.random.randint(0, 1)
scale_N3 = np.random.randint(0, 1)
scale_N4 = np.random.randint(0, 1)
scale_N5 = np.random.randint(0, 1)
scale_Madj = np.random.randint(0, 1)
Madj = np.random.randint(6, 50)
thickness = np.random.uniform(10, 1000)
d_spacing = thickness*np.random.uniform(1.2,5.0)
# Let polydispersity peak around 15%; 95% < 0.4; max=100%
sigma_d = np.random.beta(1.5, 7)
sigma_d_by_d = np.random.beta(1.5, 7)
pars = dict(
thickness=thickness,
scale_N1=scale_N1,
scale_N2=scale_N2,
scale_N3=scale_N3,
scale_N4=scale_N4,
scale_N5=scale_N5,
scale_Madj=scale_Madj,
Madj=Madj,
d_spacing=d_spacing,
sigma_d=sigma_d,
sigma_d=sigma_d_by_d,
)
return pars

# ER defaults to 0.0
# VR defaults to 1.0

#TODO check the demo plot
demo = dict(scale=1, background=0,
thickness=33, scale_N1=0.5, scale_N2=0.0, scale_N3=0.5, scale_N4=0.0, d_spacing=250, sigma_d=0.2,
sld=1.0, sld_solvent=6.34,
thickness_pd=0.2, thickness_pd_n=40, plus_monolayer=0.5)
thickness=33, scale_N1=0.5, scale_N2=0.0, scale_N3=0.5, scale_N4=0.0, scale_Madj=0.75, Madj=40,
d_spacing=250, sigma_d_by_d=0.2,
sld=1.0, sld_solvent=6.34)

# needs new unit tests here
#TODO needs new unit tests here
'''
# this from a different model - edit it!
tests = [
[{'scale': 1.0, 'background': 0.0, 'thickness': 33., 'Nlayers': 20.0,
'd_spacing': 250., 'sigma_d': 0.2, 'sld': 1.0,
'sld_solvent': 6.34, 'thickness_pd': 0.0, 'thickness_pd_n': 40},
[0.001, 0.215268], [21829.3, 0.00487686]],
]
'''
# ADDED by: RKH ON: June 11, 2020


0 comments on commit d84ebe1

Please sign in to comment.