Changeset c9e31e2 in sasmodels for sasmodels/models/stickyhardsphere.py


Ignore:
Timestamp:
Mar 6, 2015 10:54:22 AM (9 years ago)
Author:
richardh
Branches:
master, core_shell_microgels, costrafo411, magnetic_model, release_v0.94, release_v0.95, ticket-1257-vesicle-product, ticket_1156, ticket_1265_superball, ticket_822_more_unit_tests
Children:
9f6f2f8, ab87a12
Parents:
d60b433 (diff), 3c56da87 (diff)
Note: this is a merge changeset, the changes displayed below correspond to the merge itself.
Use the (diff) links above to see all the changes relative to each parent.
Message:

Merge branch 'master' of ​https://github.com/SasView/sasmodels

File:
1 edited

Legend:

Unmodified
Added
Removed

  • sasmodels/models/stickyhardsphere.py

    r9cb1415 r3c56da87  
    11# Note: model title and parameter table are inserted automatically 
    2 r"""This calculates the interparticle structure factor for a hard sphere fluid with a narrow attractive well. A perturbative 
    3 solution of the Percus-Yevick closure is used. The strength of the attractive well is described in terms of "stickiness" 
    4 as defined below. The returned value is a dimensionless structure factor, *S(q)*. 
     2r""" 
     3This calculates the interparticle structure factor for a hard sphere fluid 
     4with a narrow attractive well. A perturbative solution of the Percus-Yevick 
     5closure is used. The strength of the attractive well is described in terms 
     6of "stickiness" as defined below. The returned value is a dimensionless 
     7structure factor, *S(q)*. 
    58 
    6 The perturb (perturbation parameter), |epsilon|, should be held between 0.01 and 0.1. It is best to hold the 
    7 perturbation parameter fixed and let the "stickiness" vary to adjust the interaction strength. The stickiness, |tau|, 
    8 is defined in the equation below and is a function of both the perturbation parameter and the interaction strength. 
    9 |tau| and |epsilon| are defined in terms of the hard sphere diameter (|sigma| = 2\*\ *R*\ ), the width of the square 
    10 well, |bigdelta| (same units as *R*), and the depth of the well, *Uo*, in units of kT. From the definition, it is clear 
    11 that smaller |tau| means stronger attraction. 
     9The perturb (perturbation parameter), |epsilon|, should be held between 0.01 
     10and 0.1. It is best to hold the perturbation parameter fixed and let 
     11the "stickiness" vary to adjust the interaction strength. The stickiness, 
     12|tau|, is defined in the equation below and is a function of both the 
     13perturbation parameter and the interaction strength. |tau| and |epsilon| 
     14are defined in terms of the hard sphere diameter (|sigma| = 2\*\ *R*\ ), the 
     15width of the square well, |bigdelta| (same units as *R*), and the depth of 
     16the well, *Uo*, in units of kT. From the definition, it is clear that 
     17smaller |tau| means stronger attraction. 
    1218 
    1319.. image:: img/stickyhardsphere_228.PNG 
    … …  
    1723.. image:: img/stickyhardsphere_229.PNG 
    1824 
    19 The Percus-Yevick (PY) closure was used for this calculation, and is an adequate closure for an attractive interparticle 
    20 potential. This solution has been compared to Monte Carlo simulations for a square well fluid, with good agreement. 
     25The Percus-Yevick (PY) closure was used for this calculation, and is an 
     26adequate closure for an attractive interparticle potential. This solution 
     27has been compared to Monte Carlo simulations for a square well fluid, with 
     28good agreement. 
    2129 
    22 The true particle volume fraction, |phi|, is not equal to *h*, which appears in most of the reference. The two are 
    23 related in equation (24) of the reference. The reference also describes the relationship between this perturbation 
    24 solution and the original sticky hard sphere (or adhesive sphere) model by Baxter. 
     30The true particle volume fraction, |phi|, is not equal to *h*, which appears 
     31in most of the reference. The two are related in equation (24) of the 
     32reference. The reference also describes the relationship between this 
     33perturbation solution and the original sticky hard sphere (or adhesive 
     34sphere) model by Baxter. 
    2535 
    26 NB: The calculation can go haywire for certain combinations of the input parameters, producing unphysical solutions - in 
    27 this case errors are reported to the command window and the *S(q)* is set to -1 (so it will disappear on a log-log 
    28 plot). Use tight bounds to keep the parameters to values that you know are physical (test them) and keep nudging them 
    29 until the optimization does not hit the constraints. 
     36NB: The calculation can go haywire for certain combinations of the input 
     37parameters, producing unphysical solutions - in this case errors are 
     38reported to the command window and the *S(q)* is set to -1 (so it will 
     39disappear on a log-log plot). Use tight bounds to keep the parameters to 
     40values that you know are physical (test them) and keep nudging them until 
     41the optimization does not hit the constraints. 
    3042 
    31 In sasview the effective radius will be calculated from the parameters used in the form factor P(Q) that this  
    32 S(Q) is combined with. 
     43In sasview the effective radius will be calculated from the parameters 
     44used in the form factor P(Q) that this S(Q) is combined with. 
    3345 
    34 For 2D data: The 2D scattering intensity is calculated in the same way as 1D, where the *q* vector is defined as 
     46For 2D data: The 2D scattering intensity is calculated in the same way 
     47as 1D, where the *q* vector is defined as 
    3548 
    3649.. math:: 
    … …  
    5770 
    5871# TODO: refactor so that we pull in the old sansmodels.c_extensions 
    59   
    60 from numpy import pi, inf 
     72 
     73from numpy import inf 
    6174 
    6275name = "stickyhardsphere" 
    … …  
    6477description = """\ 
    6578        [Sticky hard sphere structure factor, with Percus-Yevick closure] 
    66         Interparticle structure factor S(Q)for a hard sphere fluid with  
     79        Interparticle structure factor S(Q)for a hard sphere fluid with 
    6780                a narrow attractive well. Fits are prone to deliver non-physical 
    68                 parameters, use with care and read the references in the full manual.  
    69                 In sasview the effective radius will be calculated from the  
     81                parameters, use with care and read the references in the full manual. 
     82                In sasview the effective radius will be calculated from the 
    7083                parameters used in P(Q). 
    7184""" 
     85category = "structure-factor" 
    7286 
    7387parameters = [ 
    74 #   [ "name", "units", default, [lower, upper], "type", 
    75 #     "description" ], 
    76     [ "effect_radius", "Ang", 50.0, [0, inf], "volume", 
    77       "effective radius of hard sphere" ], 
    78     [ "volfraction", "", 0.2, [0, 0.74], "", 
    79       "volume fraction of hard spheres" ], 
    80     [ "perturb", "", 0.05, [0.01, 0.1], "", 
    81       "perturbation parameter, epsilon" ], 
    82     [ "stickiness", "",  0.20, [-inf,inf], "", 
    83       "stickiness, tau" ], 
     88    #   [ "name", "units", default, [lower, upper], "type", 
     89    #     "description" ], 
     90    ["effect_radius", "Ang", 50.0, [0, inf], "volume", 
     91     "effective radius of hard sphere"], 
     92    ["volfraction", "", 0.2, [0, 0.74], "", 
     93     "volume fraction of hard spheres"], 
     94    ["perturb", "", 0.05, [0.01, 0.1], "", 
     95     "perturbation parameter, epsilon"], 
     96    ["stickiness", "", 0.20, [-inf, inf], "", 
     97     "stickiness, tau"], 
    8498    ] 
    85          
     99 
    86100# No volume normalization despite having a volume parameter 
    87101# This should perhaps be volume normalized? 
    … …  
    91105 
    92106Iq = """ 
    93         double onemineps,eta; 
    94         double sig,aa,etam1,etam1sq,qa,qb,qc,radic; 
    95         double lam,lam2,test,mu,alpha,beta; 
    96         double kk,k2,k3,ds,dc,aq1,aq2,aq3,aq,bq1,bq2,bq3,bq,sq; 
    97          
    98         onemineps = 1.0-perturb; 
    99         eta = volfraction/onemineps/onemineps/onemineps; 
    100          
    101         sig = 2.0 * effect_radius; 
    102         aa = sig/onemineps; 
    103         etam1 = 1.0 - eta; 
    104         etam1sq=etam1*etam1; 
    105         //C 
    106         //C  SOLVE QUADRATIC FOR LAMBDA 
    107         //C 
    108         qa = eta/12.0; 
    109         qb = -1.0*(stickiness + eta/etam1); 
    110         qc = (1.0 + eta/2.0)/etam1sq; 
    111         radic = qb*qb - 4.0*qa*qc; 
    112         if(radic<0) { 
    113                 //if(x>0.01 && x<0.015) 
    114                 //      Print "Lambda unphysical - both roots imaginary" 
    115                 //endif 
    116                 return(-1.0); 
    117         } 
    118         //C   KEEP THE SMALLER ROOT, THE LARGER ONE IS UNPHYSICAL 
    119         lam = (-1.0*qb-sqrt(radic))/(2.0*qa); 
    120         lam2 = (-1.0*qb+sqrt(radic))/(2.0*qa); 
    121         if(lam2<lam) { 
    122                 lam = lam2; 
    123         } 
    124         test = 1.0 + 2.0*eta; 
    125         mu = lam*eta*etam1; 
    126         if(mu>test) { 
    127                 //if(x>0.01 && x<0.015) 
    128                 // Print "Lambda unphysical mu>test" 
    129                 //endif 
    130                 return(-1.0); 
    131         } 
    132         alpha = (1.0 + 2.0*eta - mu)/etam1sq; 
    133         beta = (mu - 3.0*eta)/(2.0*etam1sq); 
    134         //C 
    135         //C   CALCULATE THE STRUCTURE FACTOR 
    136         //C 
    137         kk = q*aa; 
    138         k2 = kk*kk; 
    139         k3 = kk*k2; 
    140         SINCOS(kk,ds,dc); 
    141         //ds = sin(kk); 
    142         //dc = cos(kk); 
    143         aq1 = ((ds - kk*dc)*alpha)/k3; 
    144         aq2 = (beta*(1.0-dc))/k2; 
    145         aq3 = (lam*ds)/(12.0*kk); 
    146         aq = 1.0 + 12.0*eta*(aq1+aq2-aq3); 
    147         // 
    148         bq1 = alpha*(0.5/kk - ds/k2 + (1.0 - dc)/k3); 
    149         bq2 = beta*(1.0/kk - ds/k2); 
    150         bq3 = (lam/12.0)*((1.0 - dc)/kk); 
    151         bq = 12.0*eta*(bq1+bq2-bq3); 
    152         // 
    153         sq = 1.0/(aq*aq +bq*bq); 
    154          
    155         return(sq); 
     107    double onemineps,eta; 
     108    double sig,aa,etam1,etam1sq,qa,qb,qc,radic; 
     109    double lam,lam2,test,mu,alpha,beta; 
     110    double kk,k2,k3,ds,dc,aq1,aq2,aq3,aq,bq1,bq2,bq3,bq,sq; 
     111 
     112    onemineps = 1.0-perturb; 
     113    eta = volfraction/onemineps/onemineps/onemineps; 
     114 
     115    sig = 2.0 * effect_radius; 
     116    aa = sig/onemineps; 
     117    etam1 = 1.0 - eta; 
     118    etam1sq=etam1*etam1; 
     119    //C 
     120    //C  SOLVE QUADRATIC FOR LAMBDA 
     121    //C 
     122    qa = eta/12.0; 
     123    qb = -1.0*(stickiness + eta/etam1); 
     124    qc = (1.0 + eta/2.0)/etam1sq; 
     125    radic = qb*qb - 4.0*qa*qc; 
     126    if(radic<0) { 
     127        //if(x>0.01 && x<0.015) 
     128        //      Print "Lambda unphysical - both roots imaginary" 
     129        //endif 
     130        return(-1.0); 
     131    } 
     132    //C   KEEP THE SMALLER ROOT, THE LARGER ONE IS UNPHYSICAL 
     133    lam = (-1.0*qb-sqrt(radic))/(2.0*qa); 
     134    lam2 = (-1.0*qb+sqrt(radic))/(2.0*qa); 
     135    if(lam2<lam) { 
     136        lam = lam2; 
     137    } 
     138    test = 1.0 + 2.0*eta; 
     139    mu = lam*eta*etam1; 
     140    if(mu>test) { 
     141        //if(x>0.01 && x<0.015) 
     142        // Print "Lambda unphysical mu>test" 
     143        //endif 
     144        return(-1.0); 
     145    } 
     146    alpha = (1.0 + 2.0*eta - mu)/etam1sq; 
     147    beta = (mu - 3.0*eta)/(2.0*etam1sq); 
     148    //C 
     149    //C   CALCULATE THE STRUCTURE FACTOR 
     150    //C 
     151    kk = q*aa; 
     152    k2 = kk*kk; 
     153    k3 = kk*k2; 
     154    SINCOS(kk,ds,dc); 
     155    //ds = sin(kk); 
     156    //dc = cos(kk); 
     157    aq1 = ((ds - kk*dc)*alpha)/k3; 
     158    aq2 = (beta*(1.0-dc))/k2; 
     159    aq3 = (lam*ds)/(12.0*kk); 
     160    aq = 1.0 + 12.0*eta*(aq1+aq2-aq3); 
     161    // 
     162    bq1 = alpha*(0.5/kk - ds/k2 + (1.0 - dc)/k3); 
     163    bq2 = beta*(1.0/kk - ds/k2); 
     164    bq3 = (lam/12.0)*((1.0 - dc)/kk); 
     165    bq = 12.0*eta*(bq1+bq2-bq3); 
     166    // 
     167    sq = 1.0/(aq*aq +bq*bq); 
     168 
     169    return(sq); 
    156170""" 
     171 
    157172Iqxy = """ 
    158     // never called since no orientation or magnetic parameters. 
    159     return -1.0; 
     173    return Iq(sqrt(qx*qx+qy*qy), IQ_PARAMETERS); 
    160174    """ 
    161175 
    … …  
    165179oldname = 'StickyHSStructure' 
    166180oldpars = dict() 
    167 demo = dict(effect_radius = 200,volfraction = 0.2,perturb=0.05,stickiness=0.2,effect_radius_pd = 0.1,effect_radius_pd_n = 40) 
     181demo = dict(effect_radius=200, volfraction=0.2, perturb=0.05, 
     182            stickiness=0.2, effect_radius_pd=0.1, effect_radius_pd_n=40) 
    168183 
    169  
    170  
Note: See TracChangeset for help on using the changeset viewer.