Discrepancy in code for polymer_excl_volume

[smk78]

Reported by pkienzle:

Just went through the code for polymer_excl_volume. The result should match older versions of SasView, but this does not match the documentation, nor does it match Igor. For some reason, when Jae He added the model he scaled the gammainc term by gamma. Maybe Paul or Mat has some memory of why?

Here the code that Jae He added:

    $ git show -p ec1c5540
    2010-10-04 ec1c5540 added polymerexclVolume model and its test Jae Cho
    diff --git a/sansmodels/src/sans/models/PolymerExclVolume.py b/sansmodels/src/sans/models/PolymerExclVolume.py
    ...
    +        nu = 1.0 / mm
    +    
    +        Xx = x * x * rg * rg *(2.0 * nu + 1.0) * (2.0 * nu + 2.0) / 6.0
    +        onu = 1.0 / nu
    +        o2nu = 1.0 /(2.0 * nu)
    +        Ps =(1.0 / (nu * pow(Xx,o2nu))) * (gamma(o2nu)*gammainc(o2nu,Xx) - \
    +                        1.0 / pow(Xx,o2nu) * gamma(onu)*gammainc(onu,Xx))
    +
    +        if x == 0:
    +            Ps = 1.0
    +
    +        return (sc * Ps + bg);

Here's what is in Igor 7.12, Analysis/Models/NewModels_2010/PolymerExcludVol_v40.ipf:

        nu=1/m
        Xx=qval^2*Rg^2*(2*nu+1)*(2*nu+2)/6
        onu=1.0/(nu)
        o2nu=1.0/(2*nu)
        Ps=(1/(nu*Xx^o2nu))*(gammaInc(o2nu,Xx,0)-(1/Xx^o2nu)*gammaInc(onu,Xx,0))
        Debye=(2/Xx^2)*(exp(-Xx)-1+Xx)
        
        if(qval == 0)
                Ps = 1
        endif
        
        inten = I0*Ps + bgd 

Igor gammaInc(a, x, 0) =  gamma(a) - gammaInc(a, x) =  int_0^x e^-t t^(a-1) dt = scipy.special.gammainc(a, x).

Note: The current version looks different from Jae He's because I updated it during one of my cleanup passes. I'm pretty sure that I didn't change the result, but probably should revert it so that the equations in the code better match the equations in the documentation.

[smk78]

Hammouda reference for polymer_excl_volume model

[smk78]

SasView documentation for polymer_excl_volume model references:

H Benoit, Comptes Rendus, 245 (1957) 2244-2247

B Hammouda, SANS from Homogeneous Polymer Mixtures - A Unified Overview, Advances in Polym. Sci. 106(1993) 87-133

Second of these is attached to this ticket.

[krzywon]

From Mike Hore:

• The reason that gammainc is scaled by gamma in the SasView function is that gammainc in scipy is normalized by gamma, so that is 100% correct to do.

• Looking at the function in the version of SasView (4.1.2) that I have installed currently, the issue is in the following block of code:

  usub = (q*rg)**2 * (2.0/porod_exp + 1.0) * (2.0/porod_exp + 2.0)/6.0
  with errstate(divide='ignore', invalid='ignore'):
      upow = power(usub, -0.5*porod_exp)
      result= (porod_exp*upow *
              (gamma(0.5*porod_exp)*gammainc(0.5*porod_exp, usub) -
              upow*gamma(porod_exp)*gammainc(porod_exp, usub)))
  

• I think it should be something like the following instead:

  usub = (q*rg)**2 * (2.0/porod_exp + 1.0) * (2.0/porod_exp + 2.0)/6.0
  with errstate(divide='ignore', invalid='ignore'):
      upow1 = power(usub, -0.5*porod_exp)
      upow2 = power(usub, -porod_exp)
      result= (porod_exp*upow1 *
              (gamma(0.5*porod_exp)*gammainc(0.5*porod_exp,
              usub)upow2*gamma(porod_exp)*gammainc(porod_exp, usub)))
  

• Because U is raised to a different power for the first and second terms in the correct version of the function (½*nu for the first, and 1/nu for the second).

A second email after he sent the above:

• I see where the disagreement is from. Boualem’s paper from 1993 has the form that’s currently in SasView, but more recent publications like ​https://www.ncnr.nist.gov/staff/hammouda/publications/2017_hammouda_kim_j_mol_liq.pdf and mine ​https://www.ncnr.nist.gov/staff/hammouda/publications/2013_hore_et_al_macromolecules.pdf have the form that I’m talking about. So I guess it’s a matter of which is correct!

[smk78]

Interestingly, ​https://www.ncnr.nist.gov/staff/hammouda/publications/2013_hore_et_al_macromolecules.pdf references the polymer with excluded volume form factor with exactly the same two references that are in the SasView documentation!!!

[krzywon]

I just spoke with Boualem and we went through the math together. In the 1993 paper, [1/vX^(1/2v)] has been factored out. When multiplied through, the resulting equation is identical to the equation in the 2013 and 2017 papers. The calculation currently within SasView matches that in the 1993 paper so I think the model might be correct.

[smk78]

So can we close this ticket now, or is there still an issue?

[krzywon]

I contacted Mike Hore and he agrees, the mathematics is correct. Yes, we can close this ticket.

[smk78]

I have pushed some updates to the documentation to Master to reflect what has been learnt from this exchange.

[Paul Kienzle <pkienzle@…>]

In 547c6f0c4f0182017463ec1d667299383fe53598/sasmodels:

polymer_excl_volume: update eq. to match paper. Refs #1073

[Paul Kienzle <pkienzle@…>]

In aa900158f288ffb14cd35ce414d65dcc836750dd/sasmodels:

polymer_excl_volume: note in the code that scipy gammainc is regularized, hence the gamma scale factor. Refs #1073