Hi Tim, I'm not sure I understand your argument either. Anyway, I hope this Ralf's paper (and references therein) will make it more clear:
http://cci.lbl.gov/~rwgk/my_papers/CCN_2011_01_electron_density.pdf All the best, Pavel On Thu, Sep 20, 2012 at 5:19 AM, Ian Tickle <ianj...@gmail.com> 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 <t...@shelx.uni-ac.gwdg.de> wrote: > > -----BEGIN PGP SIGNED MESSAGE----- > > Hash: SHA1 > > > > 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 <jmhol...@lbl.gov> > >> 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 > >>>>>> <t...@shelx.uni-ac.gwdg.de> 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/ > > > > iD8DBQFQWuRdUxlJ7aRr7hoRAix4AJ9u/dT5RQuuL9wY0+2BTQre9TdsywCeJ5Jg > > uWoFQGip/L4ZTbGQvaMAIHQ= > > =Gdhk > > -----END PGP SIGNATURE----- >
