search

SasView / sasmodels

Package for calculation of small angle scattering models using OpenCL.
BSD 3-Clause "New" or "Revised" License
15 stars 27 forks source link

hayter_msa S(Q) - possible limitations or numerical instabilities (Trac #859) #143

Open RichardHeenan opened 5 years ago

RichardHeenan commented 5 years ago

This requires investigation. User Paul Fitzgerald (Univ Sydney) has supplied comprehensive examples, see below, and .xls on original email, attached here. NOTE that Steve Kline in March 2015 reported to Paul Butler, that sasview does have the rescaled MSA from Hayter.

Hello SASView help,

I've occasionally run into an instability in the Hayter_MSA routine (both in SASView and the original Igor routine), which I usually get around by fudging the ionic strength and charge. However, this is not possible with my current set of data where we intentionally and precisely controlled the ionic strength while looking at the particle charge.

Here is an example of the problem (more examples in the attached spreadsheet):

In hayter_msa if I set: radius = 20Å Temp = 298K diel = 78 volfrac = 0.001 [salt] = 0.001 M

then vary the charge from 0 to infinity I get:

0-2 => all OK 3-32 => returns -1 for all values (error code -3) 33-115 => all OK 116-195 => returns -1 for all values (error code -3) 196-inf => all OK

The error comes from the sqcoeff(ir) function that handles rescaling. Note that this is for the Igor code - but not the compiled C version! I haven't been unable to find the HayterMSA code in the most recent version of the SASview source code so I don't know if it is the same function.

Can you tell me if (a) the theory is just not applicable to my systems and I need to try something else, or (b) if there is an implementation problem?

I did find the following in Hayter & Hensen's 1982 paper for rescaling (which I'm still trying to work through) but didn't know if it was applying to me or not: In the weak screening limit, the MSA solution of Hayter & Penfold [3] requires increasing computational precision because of the wide dynamic range of the coefficients. Under these conditions it is more efficient to use the simpler MSA solution for unscreened charged hard spheres [10] as a reference system, and to treat the weak screening as a perturbation [16]. The presence of the effective hard core allows an application of the 'optimized random phase approximation' (ORPA) [17] to express the structure factor S(K) of the screened system in terms of the structure factor So(K) of the unscreened reference system. In practice, S(K) and So(K) turn out to be almost indistinguishable for Ka < 0.5.

That being said, even if it turns out that the theory is not applicable, I think that there is a problem in the way that SASview deals with this problem. When the hayter_msa returns values of -1 this is an error, but SASView returns (down in the bottom left corner) "Computation completed!" instead of "error".

From a fitting point of view, SASView just sees a massive increase in chi square and moves in the opposite direction. This can be a problem especially when the solution is on the other side of the instability, which can be narrow. For example, if you try the above settings with vol frac = 0.01 then the instability is only from 4 to 6.

I also think that once the underlying problem should be documented in the SASView help file. We've had students run into this problem and being completely lost as to what to do.

Sorry if this creates problems.

Regards,

Paul

Migrated from http://trac.sasview.org/ticket/859

{
    "status": "new",
    "changetime": "2019-03-28T09:21:09",
    "_ts": "2019-03-28 09:21:09.543786+00:00",
    "description": "This requires investigation. User Paul Fitzgerald (Univ Sydney) has supplied comprehensive examples, see below, and .xls on original email, attached here. NOTE that Steve Kline in March 2015 reported to Paul Butler, that sasview does have the rescaled MSA from Hayter.\n\nHello SASView help,\n\nI've occasionally run into an instability in the Hayter_MSA routine (both in SASView and the original Igor routine), which I usually get around by fudging the ionic strength and charge.  However, this is not possible with my current set of data where we intentionally and precisely controlled the ionic strength while looking at the particle charge.   \n\nHere is an example of the problem (more examples in the attached spreadsheet):\n\nIn hayter_msa if I set: \nradius = 20\u00c5\nTemp = 298K\ndiel = 78\nvolfrac = 0.001 \n[salt] = 0.001 M\n\nthen vary the charge from 0 to infinity I get:\n\n0-2        => all OK\n3-32      => returns -1 for all values (error code -3)\n33-115  => all OK\n116-195 => returns -1 for all values (error code -3)\n196-inf  => all OK\n\nThe error comes from the sqcoeff(ir) function that handles rescaling.  Note that this is for the Igor code - but not the compiled C version!  I haven't been unable to find the HayterMSA code in the most recent version of the SASview source code so I don't know if it is the same function.\n\nCan you tell me if \n(a) the theory is just not applicable to my systems and I need to try something else, or \n(b) if there is an implementation problem?\n\nI did find the following in Hayter & Hensen's 1982 paper for rescaling (which I'm still trying to work through) but didn't know if it was applying to me or not:\nIn the weak screening limit, the MSA solution of Hayter & Penfold [3] requires increasing computational precision because of the wide dynamic range of the coefficients. Under these conditions it is more efficient to use the simpler MSA solution for unscreened charged hard spheres [10] as a reference system, and to treat the weak screening as a perturbation [16]. The presence of the effective hard core allows an application of the 'optimized random phase approximation' (ORPA) [17] to express the structure factor S(K) of the screened system in terms of the structure factor So(K) of the unscreened reference system. In practice, S(K) and So(K) turn out to be almost indistinguishable for Ka < 0.5.\n\nThat being said, even if it turns out that the theory is not applicable, I think that there is a problem in the way that SASview deals with this problem.  When the hayter_msa returns values of -1 this is an error, but SASView returns (down in the bottom left corner) \"Computation completed!\" instead of \"error\".  \n\nFrom a fitting point of view, SASView just sees a massive increase in chi square and moves in the opposite direction.  This can be a problem especially when the solution is on the other side of the instability, which can be narrow.  For example, if you try the above settings with vol frac = 0.01 then the instability is only from 4 to 6.\n\nI also think that once the underlying problem should be documented in the SASView help file.  We've had students run into this problem and being completely lost as to what to do.\n\nSorry if this creates problems.  \n\nRegards,\n\nPaul\n",
    "reporter": "richardh",
    "cc": "",
    "resolution": "",
    "workpackage": "SasModels Model Issues",
    "time": "2017-02-22T10:56:00",
    "component": "sasmodels",
    "summary": "hayter_msa S(Q) -  possible limitations or numerical instabilities",
    "priority": "minor",
    "keywords": "",
    "milestone": "SasView 4.3.0",
    "owner": "",
    "type": "task"
}
RichardHeenan commented 5 years ago

Trac update at 2017/02/22 10:58:11:

further examples for hayter_msa from Paul Fitzgerald

pkienzle commented 5 years ago

Trac update at 2017/03/01 13:55:13: pkienzle commented:

Some GPU drivers want array lengths to be a multiple of 4. Pad gMSAWave with 0 up to length 20.

RichardHeenan commented 5 years ago

Trac update at 2017/03/01 17:59:04: richardh commented:

thank you, will keep that in minds.

Quick tests shows similar issues with Penfold's routine in FISH, suspect we are up against limitations of the theory for long range interactions at low concentrations, or possibly the way the theory has been implemented. On-going discussions continue directly with Paul Fitzgerald in Sydney.

pkienzle commented 5 years ago

Trac update at 2017/08/02 07:05:57: pkienzle commented:

Random parameter generation is broken for heyter msa since too many parameter sets lead to unstable calculations.

butlerpd commented 5 years ago

Trac update at 2017/10/27 17:22:21: butler changed milestone from "SasView 4.2.0" to "SasView 4.3.0"

RichardHeenan commented 5 years ago

Trac update at 2018/08/06 09:25:34: richardh commented:

The Hayter-Penfold rmsa S(Q) has a lot of "instabilities" which may be due to the physical limitations and/or numerical approximations used in generating it. Paul Butler and I had some disscussions by emaill with Paul Fitzgerald in Sydney in Feb & March 2017 regarding issues he had seen. Jeff Penfold's view is that this S(Q) should not be pushed too far. For the record, as this is still ongoing, an edited email string is below, so start at bottom.

The original Hayter and Penfold MSA paper is MOLECULAR PHYSICS, 1981, VOL. 42, No. 1, 109-118    (http://dx.doi.org/10.1080/00268978100100091)

But, the paper describing rescaling was by Jean-Pierre Hansen and Hayter a year later MOLECULAR PHYSICS, 1982, VOL. 46, No. 3, 651-656    (http://dx.doi.org/10.1080/00268978200101471)

Can anyone explain the following from the original paper, just under equation 18 on page 113:

"Note that in general, however, the MSA fails at low density; letting n -> 0 in (4) and using (5 a) yields g(x) -> 1 - U(x)/kT for x > 1. Since U(x) is generally larger than thermal energies for small interparticle separations, g(x) will generally be negative (and hence unphysical) near the particle at very low densities."

Failure at low density sounds a lot like the problem that I am having.  Although I think that I've found a work around for my system, but I'll keep working on it because we run into it reasonably often.

Regards,

Paul F

PS FYI I also tried the Hayter MSA in SASFit (from the PSI).  It too failed in the same way.  However, their code is suspiciously similar to the SASView code (except for the copyright notice from the PSI).

'''Sent:''' 01 March 2017 01:47 '''To:''' Heenan, Richard (STFC,RAL,ISIS) '''Cc:''' [mailto:butlerpd@udel.edu]

Dear Richard,

Thanks for that.  I had also tried changing the maximum number of iterations (itm = 400) but it didn't work for me either. 

I'll have a go at the modified code when I get a chance.  Maybe next week.

I've also been trying to read up on the theory a bit more, but I'm only part way through "Introduction to Statistical Mechanics of Fluids", which is a book chapter written mainly by John Hayter (I think).  Correct me if I am wrong but S(q) is calculated from the pair distribution function, g(r)?  And g(r) is calculated from the pair potential, v(r), using the Ornstein-Zernike equation?  But the Ornstein-Zernike equation can only be solved with a "closure relation" (still not sure what that means yet) such as the MSA.  But the MSA only works for a volume fraction greater than 0.2 because otherwise g(r) is unphysically negative at particle contact.  So Hanson and Hayter (1982) fudged a solution by increasing (rescaling) the diameter until they got g(r) = 0 at contact?  Which is the step causing me trouble.  Does that sound right?

I'll keep reading and let you know once I get bored of that and just go ask the Stat Mech guys down the hall.

Regards, Paul F

On 1/03/2017 5:00 AM, [mailto:richard.heenan@stfc.ac.uk] wrote:

Dear Paul (cc Paul B)

 As a temporary fix our own developer code is now returning NAN instead of -1, -2 etc for errors from the Hayter_msa routine (that really ought to be called Hayter_rmsa).  [now in released code, see #142 ] We are still working on finding a proper way for Sasview to flag errors,  and also to return some “derived values” such as in this case the Debye length. The NAN return gives an error message from the plotting routine, but at least it flags up that something is wrong.

As to whether the instabilities you see are due to the limitations of the method or to computational issues I am not knowledgeable enough to be able to work this out. You will notice some hard wired parameters such as itm=40 and acc=5e-06. In particular itm=40 seems to control the maximum number of iterations in a Newton loop which returns -1 if exceeded.

You are welcome to try changing these to see what happens. Having said which since a custom S(Q) is not actually recognised yet, you would have to add a new S(Q) model in situ [ rather than use the plugin model directory is !c:/users/yourname/.sasview/ ].  Sasview will recompile it each time you start the program if.

You can put your own version, say  hayter_msa2.py and .c, as attached, into the same directory as the usual models, note that I have edited lines

name = "hayter_msa2"

source = !["hayter_msa2.c"]

inside the hayter_msa2.py file. This version currently has itm=80 instead of 40, alas such a simple change did not make any difference to your numerical example with R=20, vol frac=0.001, salt=0.001 etc  Don’t know if we can try quad precision. Suspect that understanding the theory better is the way to go however.

Richard

On 22-Feb-17 9:30 PM, [mailto:richard.heenan@stfc.ac.uk] wrote:

Hi Paul,

Your careful diagnosis and test examples are extremely helpful.  According to Steve Kline (see email copied below), both Sasview and Igor have the “rescaled” MSA. My own, very old, FISH code also has Jeff Penfold’s Fortran RMSA version (likely the same as Igor acquired from John Hayter), [oops I recently re-discovered that I fudged one of the solvent parameters, so my implementation with only charge, inverse Debye length, and radius as parameters is only correct for water at 25C]. I will now try some of your good/bad parameter sets in the FISH code to see if the same happens there. It would be useful to work out whether the theory is falling down, or computational issues are hitting us, in which case we may be able to force improved precision and/or refactor calculations so they work better.

The “-1” error return certainly needs to be trapped in Sasview. In fact we were discussing this sort of thing at our team video meeting yesterday. Perhaps some of the computational experts on the help list could comment. Should the routine be returning say NAN ? [ticket was #142.]

Richard

 

'''From:''' Paul Butler [mailto:butlerpd@udel.edu] '''Sent:''' 10 March 2015 20:45 '''To:''' Heenan, Richard (STFC,RAL,ISIS) '''Subject:''' MSA vs RMSA

I noticed your working on the MSA routine about some instabilities.  I also saw you mentioning that RMSA would be better than MSA.  It turns out that recently Karen Edler had asked me about this very thing (and I couldn't remember so had to go to Steve Kline who has dealt with that question many times and had checked) and got verification that it was in fact the RMSAnot the MSA -- I edited the documentation last week to that effect.  Actually Karen was asking about the IGOR code but what !SasView currently has is from the IGOR code.  You may have more stable code to put in eventually but just wanted to let you know what I know in case people are using it and ask. [No, the code Richard has in FISH does much the same.]

The sad irony is that while I had to ask Steve Kline, it turns out that I am the one who originally "coded" that  model in IGOR.. nearly 15 years ago? I say "coded" since it was not so much coded as converting John Hayter's FORTRAN code (from him directly as I recall?) into IGOR (meaning dealing with block statements and such). 

 Cheers

 Paul B

 On Tue, Feb 21, 2017 at 8:31 PM, Paul !FitzGerald [mailto:paul.fitzgerald@sydney.edu.au] wrote:

Hello SASView help,

I've occasionally run into an instability in the Hayter_MSA routine (both in SASView and the original Igor routine), which I usually get around by fudging the ionic strength and charge.  However, this is not possible with my current set of data where we intentionally and precisely controlled the ionic strength while looking at the particle charge.  

Here is an example of the problem (more examples in the attached spreadsheet):

In hayter_msa if I set: radius = 20Å Temp = 298K diel = 78 volfrac = 0.001 [salt] = 0.001 M

then vary the charge from 0 to infinity I get:

0-2        => all OK 3-32      => returns -1 for all values (error code -3) 33-115  => all OK 116-195 => returns -1 for all values (error code -3) 196-inf  => all OK

The error comes from the sqcoeff(ir) function that handles rescaling.  Note that this is for the Igor code - but not the compiled C version!  I haven't been unable to find the HayterMSA code in the most recent version of the SASview source code so I don't know if it is the same function.

Can you tell me if (a) the theory is just not applicable to my systems and I need to try something else, or (b) if there is an implementation problem?

I did find the following in [http://dx.doi.org/10.1080/00268978200101471 Hayter & Hensen's 1982 paper for rescaling] (which I'm still trying to work through) but didn't know if it was applying to me or not:

In the weak screening limit, the MSA solution of Hayter & Penfold ![3] requires increasing computational precision because of the wide dynamic range of the coefficients. Under these conditions it is more efficient to use the simpler MSA solution for unscreened charged hard spheres ![10] as a reference system, and to treat the weak screening as a perturbation ![16]. The presence of the effective hard core allows an application of the 'optimized random phase approximation' (ORPA) ![17] to express the structure factor S(K) of the screened system in terms of the structure factor So(K) of the unscreened reference system. In practice, S(K) and So(K) turn out to be almost indistinguishable for Ka < 0.5.

That being said, even if it turns out that the theory is not applicable, I think that there is a problem in the way that SASview deals with this problem.  When the hayter_msa returns values of -1 this is an error, but SASView returns (down in the bottom left corner) "Computation completed!" instead of "error". 

From a fitting point of view, SASView just sees a massive increase in chi square and moves in the opposite direction.  This can be a problem especially when the solution is on the other side of the instability, which can be narrow.  For example, if you try the above settings with vol frac = 0.01 then the instability is only from 4 to 6.

I also think that once the underlying problem should be documented in the SASView help file.  We've had students run into this problem and being completely lost as to what to do.

Sorry if this creates problems. 

Regards, '''Paul !FitzGerald''' 

butlerpd commented 5 years ago

Trac update at 2019/03/26 19:58:46: butler changed workpackage from "SasView Bug Fixing" to "SasModels Model Issues"

smk78 commented 5 years ago

Trac update at 2019/03/28 09:21:09: smk78 commented:

This ticket and SasView/sasview#1193 are now referenced in the model source code, and the model help docs note the instabilities caused by very low or very high charge values.

Is there actually any more that we need to (or can) do on this? If so, it should be specified here. If not the ticket should be closed.