-----BEGIN PGP SIGNED MESSAGE-----
Hash: SHA1
P.S.: run.py reads
format_value("%.6g",ug_n.electron_density(0, b_iso))
so I thought the output of the first line states the calculated
electron density a position 0 (0,0,0) for a Carbon atom (top lines) at
the given b_iso values.
Cheers,
Tim
On 09/20/2012 02:47 PM, Tim Gruene wrote:
> Because C has 6 electrons and without thermal vibrations (T=0/B=0)
> I thought you'd catch all six of them with a box of 1A side
> length.
>
> Is this too simple thinking?
>
> Tim
>
> On 09/20/2012 02:19 PM, Ian Tickle wrote:
>> Tim, I don't follow your argument: why should the density be
>> 6A^-3 at the centre of a C atom?
>
>> -- Ian
>
>> On 20 September 2012 10:39, Tim Gruene <[email protected]>
>> wrote: tg@slartibartfast:~/tmp$ phenix.python run.py 0.001
>> 627.413 -4.01639e+06 303880 0.1 275.984
>> 275.247 435.678 0.5 92.2049 92.206
>> 93.6615 1 47.8941 47.8936 47.9421 10
>> 3.54414 3.54415 3.5439 100 0.217171
>> 0.21717 0.21714
>
>> weird numbers. A proper description would have 6e/A^3 for a C at
>> x=(0,0,0) with B=0. How are these numbers 'not inaccurate'?
>
>> Cheers, Tim
>
>> On 09/19/2012 06:47 PM, Pavel Afonine wrote:
>>>>> Hi James,
>>>>>
>>>>> using dynamic N-Gaussian approximation to form-factor
>>>>> tables as described here (pages 27-29):
>>>>>
>>>>> http://cci.lbl.gov/publications/download/iucrcompcomm_jan2004.pdf
>>>>>
>>>>>
>>>>>
>
>>>>>
and used in Phenix since 2004, avoids both: singularity at B=0 and
>>>>> inaccurate density values (compared to the raw
>>>>> forma-factor tables) for B->0.
>>>>>
>>>>> Attached is the script that proves this point. To run,
>>>>> simply "phenix.python run.py".
>>>>>
>>>>> Pavel
>>>>>
>>>>> On Sun, Sep 16, 2012 at 11:32 PM, James Holton
>>>>> <[email protected]> wrote:
>>>>>
>>>>>> Yes, the constant term in the "5-Gaussian" structure
>>>>>> factor tables does become annoying when you try to plot
>>>>>> electron density in real space, but only if you try to
>>>>>> make the B factor zero. If the B factors are ~12 (like
>>>>>> they are in 1m1n), then the electron density 2.0 A from
>>>>>> an Fe atom is not -0.2 e-/A^3, it is 0.025 e-/A^3. This
>>>>>> is only 1% of the electron density at the center of a
>>>>>> nitrogen atom with the same B factor.
>>>>>>
>>>>>> But if you do set the B factor to zero, then the
>>>>>> electron density at the center of any atom (using the
>>>>>> 5-Gaussian model) is infinity. To put it in gnuplot-ish,
>>>>>> the structure factor of Fe (in reciprocal space) can be
>>>>>> plotted with this function:
>>>>>> Fe_sf(s)=Fe_a1*exp(-Fe_b1*s*s)**+Fe_a2*exp(-Fe_b2*s*s)+Fe_a3***
>>>>>>
>>>>>>
>
>>>>>>
exp(-Fe_b3*s*s)+Fe_a4*exp(-Fe_**b4*s*s)+Fe_c
>>>>>>
>>>>>> where: Fe_c = 1.036900; Fe_a1 = 11.769500; Fe_a2 =
>>>>>> 7.357300; Fe_a3 = 3.522200; Fe_a4 = 2.304500; Fe_b1 =
>>>>>> 4.761100; Fe_b2 = 0.307200; Fe_b3 = 15.353500; Fe_b4 =
>>>>>> 76.880501; and "s" is sin(theta)/lambda
>>>>>>
>>>>>> applying a B factor is then just multiplication by
>>>>>> exp(-B*s*s)
>>>>>>
>>>>>>
>>>>>> Since the terms are all Gaussians, the inverse Fourier
>>>>>> transform can actually be done analytically, giving the
>>>>>> real-space version, or the expression for electron
>>>>>> density vs distance from the nucleus (r):
>>>>>>
>>>>>> Fe_ff(r,B) = \
>>>>>> +Fe_a1*(4*pi/(Fe_b1+B))**1.5***safexp(-4*pi**2/(Fe_b1+B)*r*r)
>>>>>>
>>>>>>
\
>>>>>> +Fe_a2*(4*pi/(Fe_b2+B))**1.5***safexp(-4*pi**2/(Fe_b2+B)*r*r)
>>>>>>
>>>>>>
\
>>>>>> +Fe_a3*(4*pi/(Fe_b3+B))**1.5***safexp(-4*pi**2/(Fe_b3+B)*r*r)
>>>>>>
>>>>>>
\
>>>>>> +Fe_a4*(4*pi/(Fe_b4+B))**1.5***safexp(-4*pi**2/(Fe_b4+B)*r*r)
>>>>>>
>>>>>>
\ +Fe_c *(4*pi/(B))**1.5*safexp(-4*pi****2/(B)*r*r);
>>>>>>
>>>>>> Where here applying a B factor requires folding it into
>>>>>> each Gaussian term. Notice how the Fe_c term blows up
>>>>>> as B->0? This is where most of the series-termination
>>>>>> effects come from. If you want the above equations for
>>>>>> other atoms, you can get them from here:
>>>>>> http://bl831.als.lbl.gov/~**jamesh/pickup/all_atomsf.**gnuplot<http://bl831.als.lbl.gov/~jamesh/pickup/all_atomsf.gnuplot>
>>>>>>
>>>>>>
>
>>>>>>
>>>>>>
> http://bl831.als.lbl.gov/~**jamesh/pickup/all_atomff.**gnuplot<http://bl831.als.lbl.gov/~jamesh/pickup/all_atomff.gnuplot>
>>>>>>
>>>>>>
>
This "infinitely sharp spike problem" seems to have led
>>>>>> some people to conclude that a zero B factor is
>>>>>> non-physical, but nothing could be further from the
>>>>>> truth! The scattering from mono-atomic gasses is an
>>>>>> excellent example of how one can observe the B=0
>>>>>> structure factor. In fact, gas scattering is how the
>>>>>> quantum mechanical self-consistent field calculations of
>>>>>> electron clouds around atoms was experimentally verified.
>>>>>> Does this mean that there really is an infinitely sharp
>>>>>> "spike" in the middle of every atom? Of course not. But
>>>>>> there is a "very" sharp spike.
>>>>>>
>>>>>> So, the problem of "infinite density" at the nucleus is
>>>>>> really just an artifact of the 5-Gaussian formalism.
>>>>>> Strictly speaking, the "5-Gaussian" structure factor
>>>>>> representation you find in ${CLIBD}/atomsf.lib (or Table
>>>>>> 6.1.1.4 in the International Tables volume C) is nothing
>>>>>> more than a curve fit to the "true" values listed in ITC
>>>>>> volume C tables 6.1.1.1 (neutral atoms) and 6.1.1.3
>>>>>> (ions). These latter tables are the Fourier transform of
>>>>>> the "true" electron density distribution around a
>>>>>> particular atom/ion obtained from quantum mechanical
>>>>>> self-consistent field calculations (like those of Cromer,
>>>>>> Mann and many others).
>>>>>>
>>>>>> The important thing to realize is that the fit was done
>>>>>> in _reciprocal_ space, and if you look carefully at
>>>>>> tables 6.1.1.1 and 6.1.1.3, you can see that even at
>>>>>> REALLY high angle (sin(theta)/lambda = 6, or 0.083 A
>>>>>> resolution) there is still significant elastic scattering
>>>>>> from the heavier atoms. The purpose of the "constant
>>>>>> term" in the 5-Gaussian representation is to try and
>>>>>> capture this high-angle "tail", and for the really heavy
>>>>>> atoms this can be more than 5 electron equivalents. In
>>>>>> real space, this is equivalent to saying that about 5
>>>>>> electrons are located within at least ~0.03 A of the
>>>>>> nucleus. That's a very short distance, but it is also
>>>>>> not zero. This is because the first few shells of
>>>>>> electrons around things like a Uranium nucleus actually
>>>>>> are very small and dense. How, then, can we have any
>>>>>> hope of modelling heavy atoms properly without using a
>>>>>> map grid sampling of 0.01A ? Easy! The B factors are
>>>>>> never zero.
>>>>>>
>>>>>> Even for a truly infinitely sharp peak (aka a single
>>>>>> electron), it doesn't take much of a B factor to spread
>>>>>> it out to a reasonable size. For example, applying a B
>>>>>> factor of 9 to a point charge will give it a
>>>>>> full-width-half max (FWHM) of 0.8 A, the same as the
>>>>>> "diameter" of a carbon atom. A carbon atom with B=12 has
>>>>>> FWHM = 1.1 A, the same as a "point" charge with B=16.
>>>>>> Carbon at B=80 and a point with B=93 both have FWHM = 2.6
>>>>>> A. As the B factor becomes larger and larger, it tends
>>>>>> to dominate the atomic shape (looks like a single
>>>>>> Gaussian). This is why it is so hard to assign atom
>>>>>> types from density alone. In fact, with B=80, a Uranium
>>>>>> atom at 1/100th occupancy is essentially
>>>>>> indistinguishable from a hydrogen atom. That is, even a
>>>>>> modest B factor pretty much "washes out" any sharp
>>>>>> features the atoms might have. Sometimes I wonder why we
>>>>>> bother with "form factors" at all, since at modest
>>>>>> resolutions all we really need is Z (the atomic number)
>>>>>> and the B factor. But, then again, I suppose it doesn't
>>>>>> hurt either.
>>>>>>
>>>>>>
>>>>>> So, what does this have to do with series termination?
>>>>>> Series termination arises in the inverse Fourier
>>>>>> transform (making a map from structure factors).
>>>>>> Technically, the "tails" of a Gaussian never reach zero,
>>>>>> so any sort of "resolution cutoff" always introduces some
>>>>>> error into the electron density calculation. That is, if
>>>>>> you create an arbitrary electron-density map, convert it
>>>>>> into structure factors and then "fft" it back, you do
>>>>>> _not_ get the same map that you started with! How much
>>>>>> do they differ? Depends on the RMS value of the
>>>>>> high-angle structure factors that have been cut off
>>>>>> (Parseval's theorem). The "infinitely sharp spike"
>>>>>> problem exacerbates this, because the B=0 structure
>>>>>> factors do not tend toward zero as fast as a Gaussian
>>>>>> with the "atomic width" would.
>>>>>>
>>>>>> So, for a given resolution, when does the B factor get
>>>>>> "too sharp"? Well, for "protein" atoms, the following B
>>>>>> factors will introduce an rms error in the electron
>>>>>> density map equal to about 5% of the peak height of the
>>>>>> atoms when the data are cut to the following resolution:
>>>>>> d B 1.0 <5 1.5 8 2.0 27 2.5 45 3.0 65 3.5 86 4.0 >99
>>>>>>
>>>>>> smaller B factors than this will introduce more than 5%
>>>>>> error at each of these resolutions. Now, of course, one
>>>>>> is often not nearly as concerned with the average error
>>>>>> in the map as you are with the error at a particular
>>>>>> point of interest, but the above numbers can serve as a
>>>>>> rough guide. If you want to see the series-termination
>>>>>> error at a particular point in the map, you will have to
>>>>>> calculate the "true" map of your model (using a program
>>>>>> like SFALL), and then run the map back and forth through
>>>>>> the Fourier transform and resolution cutoff (such as with
>>>>>> SFALL and FFT). You can then use MAPMAN or Chimera to
>>>>>> probe the electron density at the point of interest.
>>>>>>
>>>>>> But, to answer the OP's question, I would not recommend
>>>>>> trying to do fancy map interpretation to identify a
>>>>>> mystery atom. Instead, just refine the occupancy of the
>>>>>> mystery atom and see where that goes. Perhaps jiggling
>>>>>> the rest of the molecule with "kick maps" to see how
>>>>>> stable the occupancy is. Since refinement only does
>>>>>> forward-FFTs, it is formally insensitive to series
>>>>>> termination errors. It is only map calculation where
>>>>>> series termination can become a problem.
>>>>>>
>>>>>> Thanks to Garib for clearing up that last point for me!
>>>>>>
>>>>>> -James Holton MAD Scientist
>>>>>>
>>>>>>
>>>>>>
>>>>>> On 9/15/2012 3:12 AM, Tim Gruene wrote:
>>>>>>
>>>>> Dear Ian,
>>>>>
>>>>> provided that f(s) is given by the formula in the
>>>>> Cromer/Mann article, which I believe we have agreed on,
>>>>> the inset of Fig.1 of the Science article we are talking
>>>>> about is claimed to be the graph of the function g, which I
>>>>> added as pdf to this email for better readability.
>>>>>
>>>>> Irrespective of what has been plotted in any other article
>>>>> meantioned throughout this thread, this claim is
>>>>> incorrect, given a_i, b_i, c > 0.
>>>>>
>>>>> I am sure you can figure this out yourself. My argument
>>>>> was not involving mathematical programs but only
>>>>> one-dimensional calculus.
>>>>>
>>>>> Cheers, Tim
>>>>>
>>>>> On 09/14/2012 04:46 PM, Ian Tickle wrote:
>>>>>
>>>>>>>>> On 14 September 2012 15:15, Tim Gruene
>>>>>>>>> <[email protected]> wrote:
>>>>>>>>>
>>>>>>>>>> -----BEGIN PGP SIGNED MESSAGE----- Hash: SHA1
>>>>>>>>>>
>>>>>>>>>> Hello Ian,
>>>>>>>>>>
>>>>>>>>>> your article describes f(s) as sum of four
>>>>>>>>>> Gaussians, which is not the same f(s) from
>>>>>>>>>> Cromer's and Mann's paper and the one used both
>>>>>>>>>> by Niu and me. Here, f(s) contains a constant, as
>>>>>>>>>> I pointed out to in my response, which makes the
>>>>>>>>>> integral oscillate between plus and minus
>>>>>>>>>> infinity as the upper integral border (called
>>>>>>>>>> 1/dmax in the article Niu refers to) goes to
>>>>>>>>>> infinity).
>>>>>>>>>>
>>>>>>>>>> Maybe you can shed some light on why your
>>>>>>>>>> article uses a different f(s) than Cromer/Mann.
>>>>>>>>>> This explanation might be the answer to Nius
>>>>>>>>>> question, I reckon, and feed my curiosity, too.
>>>>>>>>>>
>>>>>>>>> Tim & Niu, oops yes a small slip in the paper
>>>>>>>>> there, it should have read "4 Gaussians + constant
>>>>>>>>> term": this is clear from the ITC reference given
>>>>>>>>> and the $CLIBD/atomsf.lib table referred to. In
>>>>>>>>> practice it's actually rendered as a sum of 5
>>>>>>>>> Gaussians after you multiply the f(s) and atomic
>>>>>>>>> Biso factor terms, so unless Biso = 0 (very
>>>>>>>>> unphysical!) there is actually no constant term.
>>>>>>>>> My integral for rho(r) certainly doesn't oscillate
>>>>>>>>> between plus and minus infinity as d_min -> zero.
>>>>>>>>> If yours does then I suspect that either the Biso
>>>>>>>>> term was forgotten or if not then a bug in the
>>>>>>>>> integration routine (e.g. can it handle properly
>>>>>>>>> the point at r = 0 where the standard formula for
>>>>>>>>> the density gives 0/0?). I used QUADPACK
>>>>>>>>> (http://people.sc.fsu.edu/~**jburkardt/f_src/quadpack/**quadpack.html<http://people.sc.fsu.edu/~jburkardt/f_src/quadpack/quadpack.html>
>>>>>>>>>
>>>>>>>>>
>
>>>>>>>>>
>>>>>>>>>
> )
>>>>>>>>> which seems pretty good at taking care of such
>>>>>>>>> singularities (assuming of course that the
>>>>>>>>> integral does actually converge).
>>>>>>>>>
>>>>>>>>> Cheers
>>>>>>>>>
>>>>>>>>> -- Ian
>>>>>>>>>
>>>>>>>>> - -- - --
>>>>> Dr Tim Gruene Institut fuer anorganische Chemie
>>>>> Tammannstr. 4 D-37077 Goettingen
>>>>>
>>>>> GPG Key ID = A46BEE1A
>>>>>
>>>>>>>
>>>>>>
>
>
>
>
- --
- --
Dr Tim Gruene
Institut fuer anorganische Chemie
Tammannstr. 4
D-37077 Goettingen
GPG Key ID = A46BEE1A
-----BEGIN PGP SIGNATURE-----
Version: GnuPG v1.4.12 (GNU/Linux)
Comment: Using GnuPG with Mozilla - http://enigmail.mozdev.org/
iD8DBQFQWxEvUxlJ7aRr7hoRAiH7AKDX/RtVOUdPEKoFaOrKPa3AcnsdrwCfTcHe
Mw1RupURod928hMItnfnd+I=
=Nm8A
-----END PGP SIGNATURE-----
