Opened 3 years ago

Closed 3 years ago

#593 closed defect (fixed)

suspicious code in sasview 3.x rpa.cc

Reported by: pkienzle Owned by: butler
Priority: critical Milestone: SasView 4.1.0
Component: SasView Keywords:
Cc: Work Package: SasView Bug Fixing

Description

Initialization of Xd, Xd looks suspicious in rpa.cc line 187:

Xa=q*q*ba*ba*Na/6.0;
Xb=q*q*bb*bb*Nb/6.0;
Xc=q*q*bc*bb*Nc/6.0;
Xd=q*q*bd*bb*Nd/6.0;

The remainder of the initialization code is quite regular in its symmetries, so having Xa=f(ba*ba), Xb=f(bb*bb), but Xc=f(bc*bb) and Xd=f(bd*bb) looks like an error.

Also, initializing the unused portions of a 4x4 matrix as in 222-293 looks prone to causing ill-condition solution matrices. Why not 1 on the diagonal and zero of the diagonal?

Link to the source:

https://github.com/SasView/sasview/blob/792e6be8e13f67657de88e5fb9a9afb91a750c30/src/sas/models/c_extension/c_models/rpa.cpp

Change History (11)

comment:1 Changed 3 years ago by smk78

I don't think the definitions of Xc & Xd are suspicious: I think they are wrong!

The RPA10 model allows for up to 4 polymers (i=A/B/C/D). Xi is the reduced Rg in the evaluation of the Debye function, where the actual Rg is calculated from the segment length bi and number of segments Ni.

So Xi=q*q*bi*bi*Ni/6.0

comment:2 Changed 3 years ago by butler

Hummm… but there are also cross terms I guess since it is a mixture of homopolymer and diblocks? so verifying the entire function is important but the current reference seems rather difficult to use for this purpose. I asked Adrian Rennie if he had some suggestions and he offered the following suggestions/input which may be helpful in this regard.


Part of the problem with the current documentation in SASview is that the 10 models under this heading were probably not ever really described together in one paper and so if one wants an overview one may have to refer to a text book or perhaps Boualem's SAS document. The reference given to Ackasu et al is really quite complicated as it gives the dynamic structure factor. Actually quite a few of the early papers from de Gennes et al. using the random phase approximation were also concerned with dynamics.

For diblock copolymer melts I used the equations given in the following paper:

  1. Olvera de la Cruz, I. C. Sanchez, Microphase separation in block copolymer/homopolymer blends Macromolecules 20, (1987), 440-443.

http://dx.doi.org/10.1021/ma00168a040

For polymer solutions the older reference is helpful:

  1. Jannink, P. G. de Gennes, Quasielastic Scattering by Semidilute Polymer Solutions J. Chem. Phys. 48, (1968), 2260-2265. http://dx.doi.org/10.1063/1.1669422

Polymer mixtures (homopolymers) were described by Brochard and de Gennes:

  1. Brochard, P. G. de Gennes, Dynamics of compatible polymer mixtures Physica A 118, (1983), 289-299.

http://dx.doi.org/10.1016/0378-4371(83)90195-4

There are a number of papers by Ludwik Leibler that give the equations for phase copolymers and discuss the separation:

Ludwik Leibler Theory of Microphase Separation in Block Copolymers Macromolecules 13, (1980), 1602-1617.

http://dx.doi.org/10.1021/ma60078a047

Unfortunately I do not know of a place where all the various models and equations for mixtures, solutions and melts with hompolymers and copolymers with different architecture are set out in a uniform notation.

comment:3 Changed 3 years ago by pkienzle

The Xc, Xd portion is fixed.

The matrix initialization issue still stands.

comment:4 Changed 3 years ago by butler

  • Priority changed from minor to critical

comment:5 Changed 3 years ago by butler

  • Owner set to butler
  • Status changed from new to accepted

comment:6 Changed 3 years ago by butler

  • Milestone changed from SasView 4.0.0 to SasView 4.1.0

This should be only an instability problem. Since we don't have time to figure out what the right answer is on this we could pull it from the release but suggest that we keep. If this is less stable than it should be and we can fix the maths we can release it to 4.1 and upload a more stable version to the marketplace. At any rate this ticket is now moving to 4.1

comment:7 Changed 3 years ago by butler

Have added the test case published in IGOR which Boualem asserts is correct to the 1d test folder of SasView. It is named rpa_igor.txt and is two column data. It is for case 0, 100 points between q=0.001 and 0.5.
values are -0.0004 for Kcd (K34), b3=b4=5, L3 is 1e-12 and L4 is 0, N3=N4=1000, v3=v4=100, Phi3=Phi4=0.5 scale =1 and background does NOT exist in IGOR rpa so set to 0

NOTE: sasmodels is NOT handling units correctly at the moment. L is listed in fm but clearly is no e-12 is applied. Thus at present units are really m.

NOTE2: According to Steve Kline notes only the segment lengths, b, and xi parameters, Kxy, are variable parameters (fittable). the vol fractions, scattering lengths, and specific volumes are supposed to be input by the user and used as is (fixed). SasView is currently treating them all as variable.

NOTE3: the sum of the Phi must be 1 always - as currently set up the defaults are such that that is only true for cases where all 4 polymers are present (so NOT true for case 0 where the sum of the Phis = 0.5)

FINALLY — SasView does match at high q but deviates completely (goes crazy) at low q. Deviation is evident starting around 0.1 1/A

comment:8 Changed 3 years ago by butler

In 4329539369ac6b5b77841eb856439cff61b188a2/sasmodels:

add comments to rpa c code to help code readability. Addresses ticket
#593

comment:9 Changed 3 years ago by butler

In bb73096ada08f02e4b87abe96b34d91f5c188a01/sasmodels:

Hide volume fraction parameter of component D since it is treated as the
matrix polymer and thus must always be 1-phia-phib-phid (and so
caculcate in c code). Relates to ticket #593.

comment:10 Changed 3 years ago by butler

In acfb09411425445889900e65e4acf7e2ca31b0f0/sasmodels:

scale intensity returned from rpa.c code to be close to correct
(currently using 1e-24). However th exact power of 10 will depend on
the units of Li (the scattering lengths) that are assumed in the
computation (we are now telling users it is in fm but clearly the code
is assuming something else — need to verify with rpa authors. Once done
this should complete work on #593. The documentation could use a lot of
TLC but that is a different ticket.

comment:11 Changed 3 years ago by butler

  • Resolution set to fixed
  • Status changed from accepted to closed

In 20c856a138b7cfbe5d8eb03e1655371ef25b4b06/sasmodels:

Final rpa model corrections (correct scale and units) and verified with
one of original authors on paper. Fixes #593.

while documentation was updated minimally to be usable it could use more
thorough work as can many which is the subject of other tickets. Also
annoying is that the example plot is useless. It is however
autogenerated from "default parameters" which is not helpful in this
particular case. Will submit new ticket for that.

Note: See TracTickets for help on using tickets.