[Ifeffit] Cumulant expansion fittings

[Matt Newville]

Hi Anatoly, Scott,

Of course, XAFS *is* a one-dimensional probe, not a three-dimensional one.
At least ignoring for the moment the angular dependence of multiple
scattering, XAFS is sensitive to g(r) only.  Sadly, this is sometimes
forgotten in the literature, and one sees attempts to distinguish
"sigma^2_perpendicular" and "sigma^2_parallel", which is a good sign of a
paper that is complete nonsense.

Simply stated, there is a g(r), and this distributions has moments and
cumulants.  As Anatoly points out, one can (and many have) make models for
pair potentials (giving the interaction between any two atoms) or effective
pair potentials (more precisely called potentials of mean force and giving
the net potential seen by a pair of atoms) and from these calculate g(r).
This is not easy to do in general, and, as Anatoly points out, "success" in
this business is rare for systems that are not mono-atomic close-packed
metals.  There are actually some good calculations of g(r) based on force
constants for organic molecules as well, which tends to be a different
extreme case.  That's all to make two points.  1) a general approach to
converting model potentials into g(r) is not a solved problem, and 2) the
potentials used are still parametrized models.  Given that, I don't buy
any "proof" from discussion of pair potentials about what cumulants of g(r)
can and cannot be.

The more common approach is to parametrize g(r) instead of V_eff(r), as in
Matthew Marcus' example, and also as the GNXAS approach does.  This has the
advantage of not needing to parametrize potentials.

I agree that it seems unlikely to have C3 == 0 and a non-zero C4 for a real
system.  That's mostly because, in my experience, it is very rare to have
C4 be clearly distinguishable from 0.  On the other hand, I have seen C3 be
both positive and negative, so I'd hesitate to say that it cannot be zero
while C4 is not zero in a real system.  If you're looking for a model g(r)
that will give C3==0 and a non-zero C4, then

   g(r) = exp[-(r0+delta)^2/(2sigma2) ] + exp[-(r0-delta)^2/(2sigma2) ]

will do the trick.  The attached figure shows such a distribution, looking
fairly realistic and like a flattened Gaussian.  The script to generate
this is and calculate its moments and cumulants is:

####
####  script to generate symmetric, non-Gaussian distribution
set my.r  = range(1,4,0.01)
set r0    = 2.50
set sig   = 0.05
set delr  = 0.05
set my.gr  = gauss(my.r,r0+delr,sig) + gauss(my.r,r0-delr,sig)

#  calculate moments around the mean value (r0 here)

set norm  = vsum(my.gr)
set r0    = vsum(my.gr * my.r) / norm
set my.dr = my.r - r0
set mom_1   = vsum(my.gr * my.dr   ) / norm
set mom_2   = vsum(my.gr * my.dr^2 ) / norm
set mom_3   = vsum(my.gr * my.dr^3 ) / norm
set mom_4   = vsum(my.gr * my.dr^4 ) / norm

# calculate cumulants:
set cumul_2 = mom_2 - mom_1*mom_1
set cumul_3 = mom_3 - 3*mom_2*mom_1 + 2*mom_1^3
set cumul_4 = mom_4 - 3*mom_2*mom_2 - 4*mom_3*mom_1 +
12*mom_2*mom_1*mom_1 - 6*mom_1^4

show mom_1, mom_2, mom_3, cumul_2, cumul_3, cumul_4
newplot my.r, my.gr , xlabel = 'R (\A)', ylabel='g(R)'

######

The resulting moments and cumulants are:

   mom_1          =  -0.139956462E-14
   mom_2          =       0.005000000
   mom_3          =  -0.201455968E-16
   cumul_2        =       0.005000000
   cumul_3        =   0.847872485E-18
   cumul_4        =      -0.000012500

which, to within machine precision, means C3==0 and C4=/= 0.

Now, whether such a C4 can be measured from real data, and whether it might
be better to model this g(r) as two Gaussians would be good questions.

--Matt
-------------- next part --------------
A non-text attachment was scrubbed...
Name: cumul.png
Type: image/png
Size: 5387 bytes
Desc: not available
URL: <http://millenia.cars.aps.anl.gov/pipermail/ifeffit/attachments/20090122/3e0383f7/attachment.png>