Tim, Correct, but take care that 's' in Niu's program is 2.sin(theta)/lambda whereas in exp(-Bs^2) it is just sin(theta)/lambda (as it is in the usual expression for f(s)).
Cheers -- Ian On 17 September 2012 10:24, Tim Gruene <t...@shelx.uni-ac.gwdg.de> wrote: > -----BEGIN PGP SIGNED MESSAGE----- > Hash: SHA1 > > Dear James et al., > > so to summarise, the answer to Niu's question is that he must add a > factor of e^(-Bs^2) to the formula of Cromer/Mann and then adjust the > value of B until it matches the inset. Given that you claim > rho=0.025e/A^3 (I assume for 1/dmax approx. 0) for B=12 and the inset > shows a value of about 0.6, a somewhat higher B-value should work. > > Cheers, > Tim > > On 09/17/2012 08:32 AM, James Holton 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_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) >>>>> 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/ > > iD8DBQFQVuxEUxlJ7aRr7hoRAoUcAKD6v7hHaQqigX3HLZy/rHJ97zPoeACeM5Cp > XWoDIX07jHlKPC7eDcQAfq0= > =lTl4 > -----END PGP SIGNATURE-----
