Changeset e44432d in sasmodels for sasmodels/kernel.py


Ignore:
Timestamp:
Oct 19, 2018 3:46:26 PM (6 years ago)
Author:
Paul Kienzle <pkienzle@…>
Branches:
master, core_shell_microgels, magnetic_model, ticket-1257-vesicle-product, ticket_1156, ticket_1265_superball, ticket_822_more_unit_tests
Children:
be43e39
Parents:
2586ab72
Message:

support hollow models in structure factor calculations

File:
1 edited

Legend:

Unmodified
Added
Removed

  • sasmodels/kernel.py

    ree60aa7 re44432d  
    4343    def Iq(self, call_details, values, cutoff, magnetic): 
    4444        # type: (CallDetails, np.ndarray, np.ndarray, float, bool) -> np.ndarray 
    45         Pq, Reff = self.Pq_Reff(call_details, values, cutoff, magnetic, effective_radius_type=0) 
    46         return Pq 
     45        r""" 
     46        Returns I(q) from the polydisperse average scattering. 
     47 
     48        .. math:: 
     49 
     50            I(q) = \text{scale} \cdot P(q) + \text{background} 
     51 
     52        With the correct choice of model and contrast, setting *scale* to 
     53        the volume fraction $V_f$ of particles should match the measured 
     54        absolute scattering.  Some models (e.g., vesicle) have volume fraction 
     55        built into the model, and do not need an additional scale. 
     56        """ 
     57        _, F2, _, shell_volume, _ = self.Fq(call_details, values, cutoff, magnetic, 
     58                              effective_radius_type=0) 
     59        combined_scale = values[0]/shell_volume 
     60        background = values[1] 
     61        return combined_scale*F2 + background 
    4762    __call__ = Iq 
    4863 
    49     def Pq_Reff(self, call_details, values, cutoff, magnetic, effective_radius_type): 
     64    def Fq(self, call_details, values, cutoff, magnetic, effective_radius_type=0): 
    5065        # type: (CallDetails, np.ndarray, np.ndarray, float, bool, int) -> np.ndarray 
     66        r""" 
     67        Returns <F(q)>, <F(q)^2>, effective radius, shell volume and 
     68        form:shell volume ratio.  The <F(q)> term may be None if the 
     69        form factor does not support direct computation of $F(q)$ 
     70 
     71        $P(q) = <F^2(q)>/<V>$ is used for structure factor calculations, 
     72 
     73        .. math:: 
     74 
     75            I(q) = \text{scale} \cdot P(q) \cdot S(q) + \text{background} 
     76 
     77        For the beta approximation, this becomes 
     78 
     79        .. math:: 
     80 
     81            I(q) = \text{scale} * P (1 + <F>^2/<F^2> (S - 1)) + \text{background} 
     82                 = \text{scale}/<V> (<F^2> + <F>^2 (S - 1)) + \text{background} 
     83 
     84        $<F(q)>$ and $<F^2(q)>$ are averaged by polydispersity in shape 
     85        and orientation, with each configuration $x_k$ having form factor 
     86        $F(q, x_k)$, weight $w_k$ and volume $V_k$.  The result is: 
     87 
     88        .. math:: 
     89 
     90            P(q) = \frac{\sum w_k F^2(q, x_k) / \sum w_k}{\sum w_k V_k / \sum w_k} 
     91 
     92        The form factor itself is scaled by volume and contrast to compute the 
     93        total scattering.  This is then squared, and the volume weighted 
     94        F^2 is then normalized by volume F.  For a given density, the number 
     95        of scattering centers is assumed to scale linearly with volume.  Later 
     96        scaling the resulting $P(q)$ by the volume fraction of particles 
     97        gives the total scattering on an absolute scale. Most models 
     98        incorporate the volume fraction into the overall scale parameter.  An 
     99        exception is vesicle, which includes the volume fraction parameter in 
     100        the model itself, scaling $F$ by $\surd V_f$ so that the math for 
     101        the beta approximation works out. 
     102 
     103        By scaling $P(q)$ by total weight $\sum w_k$, there is no need to make 
     104        sure that the polydisperisity distributions normalize to one.  In 
     105        particular, any distibution values $x_k$ outside the valid domain 
     106        of $F$ will not be included, and the distribution will be implicitly 
     107        truncated.  This is controlled by the parameter limits defined in the 
     108        model (which truncate the distribution before calling the kernel) as 
     109        well as any region excluded using the *INVALID* macro defined within 
     110        the model itself. 
     111 
     112        The volume used in the polydispersity calculation is the form volume 
     113        for solid objects or the shell volume for hollow objects.  Shell 
     114        volume should be used within $F$ so that the normalizing scale 
     115        represents the volume fraction of the shell rather than the entire 
     116        form.  This corresponds to the volume fraction of shell-forming 
     117        material added to the solvent. 
     118 
     119        The calculation of $S$ requires the effective radius and the 
     120        volume fraction of the particles.  The model can have several 
     121        different ways to compute effective radius, with the 
     122        *effective_radius_type* parameter used to select amongst them.  The 
     123        volume fraction of particles should be determined from the total 
     124        volume fraction of the form, not just the shell volume fraction. 
     125        This makes a difference for hollow shapes, which need to scale 
     126        the volume fraction by the returned volume ratio when computing $S$. 
     127        For solid objects, the shell volume is set to the form volume so 
     128        this scale factor evaluates to one and so can be used for both 
     129        hollow and solid shapes. 
     130        """ 
    51131        self._call_kernel(call_details, values, cutoff, magnetic, effective_radius_type) 
    52132        #print("returned",self.q_input.q, self.result) 
    … …  
    55135        if total_weight == 0.: 
    56136            total_weight = 1. 
    57         weighted_volume = self.result[nout*self.q_input.nq + 1] 
    58         weighted_radius = self.result[nout*self.q_input.nq + 2] 
    59         effective_radius = weighted_radius/total_weight 
    60         # compute I = scale*P + background 
    61         #           = scale*(sum(w*F^2)/sum w)/(sum (w*V)/sum w) + background 
    62         #           = scale/sum (w*V) * sum(w*F^2) + background 
    63         F2 = self.result[0:nout*self.q_input.nq:nout] 
    64         scale = values[0]/(weighted_volume if weighted_volume != 0.0 else 1.0) 
    65         background = values[1] 
    66         Pq = scale*F2 + background 
    67         #print("scale",scale,background) 
    68         return Pq, effective_radius 
    69  
    70     def beta(self, call_details, values, cutoff, magnetic, effective_radius_type): 
    71         # type: (CallDetails, np.ndarray, np.ndarray, float, bool, int) -> np.ndarray 
    72         if self.dim == '2d': 
    73             raise NotImplementedError("beta not yet supported for 2D") 
    74         if not self.info.have_Fq: 
    75             raise NotImplementedError("beta not yet supported for "+self.info.id) 
    76         self._call_kernel(call_details, values, cutoff, magnetic, effective_radius_type) 
    77         total_weight = self.result[2*self.q_input.nq + 0] 
    78         if total_weight == 0.: 
    79             total_weight = 1. 
    80         weighted_volume = self.result[2*self.q_input.nq + 1] 
    81         weighted_radius = self.result[2*self.q_input.nq + 2] 
    82         volume_average = weighted_volume/total_weight 
    83         effective_radius = weighted_radius/total_weight 
    84         F2 = self.result[0:2*self.q_input.nq:2]/total_weight 
    85         F1 = self.result[1:2*self.q_input.nq:2]/total_weight 
    86         return F1, F2, volume_average, effective_radius 
     137        # Note: shell_volume == form_volume for solid objects 
     138        form_volume = self.result[nout*self.q_input.nq + 1]/total_weight 
     139        shell_volume = self.result[nout*self.q_input.nq + 2]/total_weight 
     140        effective_radius = self.result[nout*self.q_input.nq + 3]/total_weight 
     141        if shell_volume == 0.: 
     142            shell_volume = 1. 
     143        F1 = self.result[1:nout*self.q_input.nq:nout]/total_weight if nout == 2 else None 
     144        F2 = self.result[0:nout*self.q_input.nq:nout]/total_weight 
     145        return F1, F2, effective_radius, shell_volume, form_volume/shell_volume 
    87146 
    88147    def release(self): 
    … …  
    92151    def _call_kernel(self, call_details, values, cutoff, magnetic, effective_radius_type): 
    93152        # type: (CallDetails, np.ndarray, np.ndarray, float, bool, int) -> np.ndarray 
     153        """ 
     154        Call the kernel.  Subclasses defining kernels for particular execution 
     155        engines need to provide an implementation for this. 
     156        """ 
    94157        raise NotImplementedError() 
Note: See TracChangeset for help on using the changeset viewer.