Changeset 3c56da87 in sasmodels for sasmodels/models/stickyhardsphere.py


Ignore:
Timestamp:
Mar 4, 2015 10:55:38 PM (9 years ago)
Author:
Paul Kienzle <pkienzle@…>
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:
3a45c2c
Parents:
b89f519
Message:

lint cleanup

File:
1 edited

Legend:

Unmodified
Added
Removed

  • sasmodels/models/stickyhardsphere.py

    r1353f60 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 
    3 a narrow attractive well. A perturbative solution of the Percus-Yevick closure is used. 
    4 The strength of the attractive well is described in terms of "stickiness" 
    5 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)*. 
    68 
    7 The perturb (perturbation parameter), |epsilon|, should be held between 0.01 and 0.1. 
    8 It is best to hold the perturbation parameter fixed and let the "stickiness" vary to 
    9 adjust the interaction strength. The stickiness, |tau|, is defined in the equation 
    10 below and is a function of both the perturbation parameter and the interaction strength. 
    11 |tau| and |epsilon| are defined in terms of the hard sphere diameter (|sigma| = 2\*\ *R*\ ), 
    12 the width of the square well, |bigdelta| (same units as *R*), and the depth of the well, 
    13 *Uo*, in units of kT. From the definition, it is clear that smaller |tau| means stronger 
    14 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. 
    1518 
    1619.. image:: img/stickyhardsphere_228.PNG 
    … …  
    2023.. image:: img/stickyhardsphere_229.PNG 
    2124 
    22 The Percus-Yevick (PY) closure was used for this calculation, and is an adequate closure 
    23 for an attractive interparticle potential. This solution has been compared to Monte Carlo 
    24 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. 
    2529 
    26 The true particle volume fraction, |phi|, is not equal to *h*, which appears in most of 
    27 the reference. The two are related in equation (24) of the reference. The reference also 
    28 describes the relationship between this perturbation solution and the original sticky hard 
    29 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. 
    3035 
    31 NB: The calculation can go haywire for certain combinations of the input parameters, 
    32 producing unphysical solutions - in this case errors are reported to the command window and 
    33 the *S(q)* is set to -1 (so it will disappear on a log-log plot). Use tight bounds to keep 
    34 the parameters to values that you know are physical (test them) and keep nudging them until 
     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 
    3541the optimization does not hit the constraints. 
    3642 
    37 In sasview the effective radius will be calculated from the parameters used in the 
    38 form factor P(Q) that this 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. 
    3945 
    40 For 2D data: The 2D scattering intensity is calculated in the same way as 1D, where 
    41 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 
    4248 
    4349.. math:: 
    … …  
    99105 
    100106Iq = """ 
    101         double onemineps,eta; 
    102         double sig,aa,etam1,etam1sq,qa,qb,qc,radic; 
    103         double lam,lam2,test,mu,alpha,beta; 
    104         double kk,k2,k3,ds,dc,aq1,aq2,aq3,aq,bq1,bq2,bq3,bq,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; 
    105111 
    106         onemineps = 1.0-perturb; 
    107         eta = volfraction/onemineps/onemineps/onemineps; 
     112    onemineps = 1.0-perturb; 
     113    eta = volfraction/onemineps/onemineps/onemineps; 
    108114 
    109         sig = 2.0 * effect_radius; 
    110         aa = sig/onemineps; 
    111         etam1 = 1.0 - eta; 
    112         etam1sq=etam1*etam1; 
    113         //C 
    114         //C  SOLVE QUADRATIC FOR LAMBDA 
    115         //C 
    116         qa = eta/12.0; 
    117         qb = -1.0*(stickiness + eta/etam1); 
    118         qc = (1.0 + eta/2.0)/etam1sq; 
    119         radic = qb*qb - 4.0*qa*qc; 
    120         if(radic<0) { 
    121                 //if(x>0.01 && x<0.015) 
    122                 //      Print "Lambda unphysical - both roots imaginary" 
    123                 //endif 
    124                 return(-1.0); 
    125         } 
    126         //C   KEEP THE SMALLER ROOT, THE LARGER ONE IS UNPHYSICAL 
    127         lam = (-1.0*qb-sqrt(radic))/(2.0*qa); 
    128         lam2 = (-1.0*qb+sqrt(radic))/(2.0*qa); 
    129         if(lam2<lam) { 
    130                 lam = lam2; 
    131         } 
    132         test = 1.0 + 2.0*eta; 
    133         mu = lam*eta*etam1; 
    134         if(mu>test) { 
    135                 //if(x>0.01 && x<0.015) 
    136                 // Print "Lambda unphysical mu>test" 
    137                 //endif 
    138                 return(-1.0); 
    139         } 
    140         alpha = (1.0 + 2.0*eta - mu)/etam1sq; 
    141         beta = (mu - 3.0*eta)/(2.0*etam1sq); 
    142         //C 
    143         //C   CALCULATE THE STRUCTURE FACTOR 
    144         //C 
    145         kk = q*aa; 
    146         k2 = kk*kk; 
    147         k3 = kk*k2; 
    148         SINCOS(kk,ds,dc); 
    149         //ds = sin(kk); 
    150         //dc = cos(kk); 
    151         aq1 = ((ds - kk*dc)*alpha)/k3; 
    152         aq2 = (beta*(1.0-dc))/k2; 
    153         aq3 = (lam*ds)/(12.0*kk); 
    154         aq = 1.0 + 12.0*eta*(aq1+aq2-aq3); 
    155         // 
    156         bq1 = alpha*(0.5/kk - ds/k2 + (1.0 - dc)/k3); 
    157         bq2 = beta*(1.0/kk - ds/k2); 
    158         bq3 = (lam/12.0)*((1.0 - dc)/kk); 
    159         bq = 12.0*eta*(bq1+bq2-bq3); 
    160         // 
    161         sq = 1.0/(aq*aq +bq*bq); 
     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); 
    162168 
    163         return(sq); 
     169    return(sq); 
    164170""" 
    165171 
    … …  
    176182            stickiness=0.2, effect_radius_pd=0.1, effect_radius_pd_n=40) 
    177183 
    178  
    179  
Note: See TracChangeset for help on using the changeset viewer.