Re: [ccp4bb] Series termination effect calculation.

2012-09-20 Thread Tim Gruene 
-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  
>> 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 1000.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 
>  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*pi2/(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

Re: [ccp4bb] Series termination effect calculation.

2012-09-20 Thread Tim Gruene 
-BEGIN PGP SIGNED MESSAGE-
Hash: SHA1

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 
> 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 1000.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
  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*pi2/(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 (io

Re: [ccp4bb] Series termination effect calculation.

2012-09-20 Thread Pavel Afonine 
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  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  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
> > 1000.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 
> >> 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*pi2/(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

Re: [ccp4bb] Series termination effect calculation.

2012-09-20 Thread Ian Tickle 
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  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
> 1000.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 
>> 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*pi2/(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

Re: [ccp4bb] Series termination effect calculation.

2012-09-20 Thread Tim Gruene 
-BEGIN PGP SIGNED MESSAGE-
Hash: SHA1

Sorry, 6e^-/A^3 (or -6e/A^3 for charge density people) this should
have said.


On 09/20/2012 11:39 AM, Tim Gruene 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 1000.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
>>  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*pi2/(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(

Re: [ccp4bb] Series termination effect calculation.

2012-09-20 Thread Tim Gruene 
-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
1000.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 
> 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*pi2/(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 "

Re: [ccp4bb] Series termination effect calculation.

2012-09-20 Thread Pavel Afonine 
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  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*pi2/(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, a

Re: [ccp4bb] Series termination effect calculation.

2012-09-19 Thread Tim Gruene 
-BEGIN PGP SIGNED MESSAGE-
Hash: SHA1

Dear James,

Mann only fitted their data to sin\lambda/theta < 1.5, and up to there
the fit is pretty good. 30 years ago the computational means must have
been very different and what takes 5s now would have taken minutes or
hours then.

I am going to do the re-fitting, compare with a couple of test
structures where we have observed problems which might well be related
to this discussion. Fortunately with shelxl the testing is a piece of
cake and does not even require a change of the source code.

The results should be published within the year. Maybe you get an
updated atomsf.lib then with ccp4 6.1.4 ;-)

Best wishes,
Tim

On 09/19/2012 02:16 AM, James Holton wrote:
> 
> That's really interesting!  Since the fits then and now were both 
> least-squares, I wonder how Cromer & Mann could have gotten it so
> far off?  Looking at the residuals, I see that although that of
> nitrogen oscillates badly, even the worst outlier is still within
> 0.01 electrons of the Hartree-Fock values.  Perhaps 1% of an
> electron was their convergence limit?
> 
> Either way, I think it would be valuable to have a "re-fit" of the
> Table 6.1.1.1/3 values without the "c" term.  Then we can go all
> the way to B=0 without worrying about singularities.
> 
> For example, I attach here a plot of the electron density at the
> center of a nitrogen atom vs B factor (in real space).  The red
> curve is the result of a 20-Gaussian fit to the data for nitrogen
> in table 6.1.1.1 all the way out to sin(theta)/lambda = 6 (although
> 7 Gaussians is more than enough).  This "true" curve approaches
> 1000 e-/A^3 as B approaches zero, but the 5-Gaussian model using
> the Cromer-Mann coefficients form 6.1.1.4 (blue curve) starts to
> deviate when B becomes less than one, and actually goes negative
> for B < 0.1.   A simpler model (without the c term, but re-fit) is
> the green line.  Very much like what Tim suggested.
> 
> Not exactly a problem for typical macromolecular refinement, but 
> still...  I wonder what would happen if I edited my
> ${CLIBD}/atomsf.lib ?
> 
> -James Holton MAD Scientist
> 
> On 9/18/2012 6:32 AM, Tim Gruene wrote: Hello Oliver,
> 
> when you fit the values from ICA Tab 6.1.1.1 with gnuplot, the
> values of C and N become much more comparable. c(C) = 0.017 and
> especially c(N) = 0.025 > 0!!! for C: Final set of parameters
> Asymptotic Standard Error ===
> ==
> 
> a1  = 0.604126 +/- 0.02326  (3.85%) a2
> = 2.63343  +/- 0.03321  (1.261%) a3  =
> 1.52123  +/- 0.03528  (2.319%) a4  = 1.2211
> +/- 0.02225  (1.822%) b1  = 0.185807 +/-
> 0.00629  (3.385%) b2  = 14.6332  +/- 0.1355
> (0.9263%) b3  = 41.6948  +/- 0.5345
> (1.282%) b4  = 0.717984 +/- 0.01251
> (1.743%) c   = 0.0171359+/- 0.002045
> (11.93%)
> 
> for N: Final set of parametersAsymptotic Standard
> Error ===
> ==
> 
> a1  = 0.723788 +/- 0.04334  (5.988%) a2
> = 3.24589  +/- 0.04074  (1.255%) a3  =
> 1.90049  +/- 0.04422  (2.327%) a4  =
> 1.10071  +/- 0.0413   (3.752%) b1  =
> 0.157345 +/- 0.007552 (4.8%) b2  = 10.106
> +/- 0.1041   (1.03%) b3  = 30.0211  +/-
> 0.3946   (1.314%) b4  = 0.567116 +/-
> 0.01914  (3.376%) c   = 0.0252303+/-
> 0.003284 (13.01%)
> 
> In 1967, Mann only calculated to sin \theta/lambda = 0, ... 1.5,
> and their tabulated values do indeed fit decently within that
> range, but not out to 6A.
> 
> I thought this was notworthy, and I am curious which values for
> these constants refinement programs use nowadays. Maybe George,
> Garib, Pavel, and Gerard may want to comment?
> 
> Cheers, Tim
> 
> On 09/18/2012 10:11 AM, Oliver Einsle wrote:
 Hi there,
 
 I was just pointed to this thread and should comment on the 
 discussion, as actually made the plots for this paper. James
 has clarified the issue much better than I could have, and
 indeed the calculations will fail for larger Bragg angles if
 you do not assume a reasonable B-factor (I used B=10 for the
 plots).
 
 Doug Rees has pointed out at the time that for large theta
 the c-term of the Cromer/Mann approximation becomes dominant,
 and this is where chaos comes in, as the Cromer/Mann
 parameters are only derived from a fit to the actual
 HF-calculation. They are numbers without physical meaning,
 which becomes particularly obvious if you compare the
 parameters for C and N:
 
 
 C:   2.3100  20.8439   1.0200  10.2075   1.5886  0.5687
 0.8650 51.6512 0.2156 N:  12.2126  0.0057   3.1322  9.8933
 2.

Re: [ccp4bb] Series termination effect calculation.

2012-09-18 Thread James Holton 


That's really interesting!  Since the fits then and now were both 
least-squares, I wonder how Cromer & Mann could have gotten it so far 
off?  Looking at the residuals, I see that although that of nitrogen 
oscillates badly, even the worst outlier is still within 0.01 electrons 
of the Hartree-Fock values.  Perhaps 1% of an electron was their 
convergence limit?


Either way, I think it would be valuable to have a "re-fit" of the Table 
6.1.1.1/3 values without the "c" term.  Then we can go all the way to 
B=0 without worrying about singularities.


For example, I attach here a plot of the electron density at the center 
of a nitrogen atom vs B factor (in real space).  The red curve is the 
result of a 20-Gaussian fit to the data for nitrogen in table 6.1.1.1 
all the way out to sin(theta)/lambda = 6 (although 7 Gaussians is more 
than enough).  This "true" curve approaches 1000 e-/A^3 as B approaches 
zero, but the 5-Gaussian model using the Cromer-Mann coefficients form 
6.1.1.4 (blue curve) starts to deviate when B becomes less than one, and 
actually goes negative for B < 0.1.   A simpler model (without the c 
term, but re-fit) is the green line.  Very much like what Tim suggested.


Not exactly a problem for typical macromolecular refinement, but 
still...  I wonder what would happen if I edited my ${CLIBD}/atomsf.lib ?


-James Holton
MAD Scientist

On 9/18/2012 6:32 AM, Tim Gruene wrote:

-BEGIN PGP SIGNED MESSAGE-
Hash: SHA1

Hello Oliver,

when you fit the values from ICA Tab 6.1.1.1 with gnuplot, the values
of C and N become much more comparable. c(C) = 0.017 and especially
c(N) = 0.025 > 0!!!
for C:
Final set of parametersAsymptotic Standard Error
=====

a1  = 0.604126 +/- 0.02326  (3.85%)
a2  = 2.63343  +/- 0.03321  (1.261%)
a3  = 1.52123  +/- 0.03528  (2.319%)
a4  = 1.2211   +/- 0.02225  (1.822%)
b1  = 0.185807 +/- 0.00629  (3.385%)
b2  = 14.6332  +/- 0.1355   (0.9263%)
b3  = 41.6948  +/- 0.5345   (1.282%)
b4  = 0.717984 +/- 0.01251  (1.743%)
c   = 0.0171359+/- 0.002045 (11.93%)

for N:
Final set of parametersAsymptotic Standard Error
=====

a1  = 0.723788 +/- 0.04334  (5.988%)
a2  = 3.24589  +/- 0.04074  (1.255%)
a3  = 1.90049  +/- 0.04422  (2.327%)
a4  = 1.10071  +/- 0.0413   (3.752%)
b1  = 0.157345 +/- 0.007552 (4.8%)
b2  = 10.106   +/- 0.1041   (1.03%)
b3  = 30.0211  +/- 0.3946   (1.314%)
b4  = 0.567116 +/- 0.01914  (3.376%)
c   = 0.0252303+/- 0.003284 (13.01%)

In 1967, Mann only calculated to sin \theta/lambda = 0, ... 1.5, and
their tabulated values do indeed fit decently within that range, but
not out to 6A.

I thought this was notworthy, and I am curious which values for these
constants refinement programs use nowadays. Maybe George, Garib,
Pavel, and Gerard may want to comment?

Cheers,
Tim

On 09/18/2012 10:11 AM, Oliver Einsle wrote:

Hi there,

I was just pointed to this thread and should comment on the
discussion, as actually made the plots for this paper. James has
clarified the issue much better than I could have, and indeed the
calculations will fail for larger Bragg angles if you do not assume
a reasonable B-factor (I used B=10 for the plots).

Doug Rees has pointed out at the time that for large theta the
c-term of the Cromer/Mann approximation becomes dominant, and this
is where chaos comes in, as the Cromer/Mann parameters are only
derived from a fit to the actual HF-calculation. They are numbers
without physical meaning, which becomes particularly obvious if you
compare the parameters for C and N:


C:   2.3100  20.8439   1.0200  10.2075   1.5886  0.5687  0.8650
51.6512 0.2156 N:  12.2126  0.0057   3.1322  9.8933   2.0125
28.9975  1.1663  0.5826 -11.5290

The scattering factors for these are reasonably similar, but the
c-values are entirely different. The B-factor dampens this out and
this is an essential point.



For clarity: I made the plots using Waterloo Maple with the
following code:

restart; SF :=Matrix(17,9,readdata("scatter.dat",float,9));

biso := 10; e:=  1; AFF  :=
(e)->(SF[e,1]*exp(-SF[e,2]*s^2)+SF[e,3]*exp(-SF[e,4]*s^2)
+SF[e,5]*exp(-SF[e,6]*s^2)+SF[e,7]*exp(-SF[e,8]*s^2)
+SF[e,9])*exp(-biso*s^2/4);

H:=  AFF(1); C:=  AFF(2); N:=  AFF(3); Ox   :=
AFF(4); S:=  AFF(5); Fe   :=  AFF(6); Fe2  :=  AFF(7); Fe3  :=
AFF(8); Cu   :=  AFF(9); Cu1  :=  AFF(10); Cu2  :=  AFF(11); Mo
:=  AFF(12); Mo4  :=  AFF(13); Mo5  :=  AFF(14); Mo6  :=  AFF(15);

// Plot scattering factors

plot

[ccp4bb] AW: [ccp4bb] Series termination effect calculation.

2012-09-18 Thread Oliver Einsle 
Good point! With 9 parameters there must be a zillion combinations that produce 
a decent fit.  But even if they are chosen to be similar they have no physical 
meaning, right?

Cheers, Oliver.


-Ursprüngliche Nachricht-
Von: Tim Gruene [mailto:t...@shelx.uni-ac.gwdg.de]
Gesendet: Dienstag, 18. September 2012 15:32
An: Oliver Einsle
Cc: CCP4BB@JISCMAIL.AC.UK
Betreff: Re: [ccp4bb] Series termination effect calculation.

-BEGIN PGP SIGNED MESSAGE-
Hash: SHA1

Hello Oliver,

when you fit the values from ICA Tab 6.1.1.1 with gnuplot, the values of C and 
N become much more comparable. c(C) = 0.017 and especially
c(N) = 0.025 > 0!!!
for C:
Final set of parametersAsymptotic Standard Error
=====

a1  = 0.604126 +/- 0.02326  (3.85%)
a2  = 2.63343  +/- 0.03321  (1.261%)
a3  = 1.52123  +/- 0.03528  (2.319%)
a4  = 1.2211   +/- 0.02225  (1.822%)
b1  = 0.185807 +/- 0.00629  (3.385%)
b2  = 14.6332  +/- 0.1355   (0.9263%)
b3  = 41.6948  +/- 0.5345   (1.282%)
b4  = 0.717984 +/- 0.01251  (1.743%)
c   = 0.0171359+/- 0.002045 (11.93%)

for N:
Final set of parametersAsymptotic Standard Error
=====

a1  = 0.723788 +/- 0.04334  (5.988%)
a2  = 3.24589  +/- 0.04074  (1.255%)
a3  = 1.90049  +/- 0.04422  (2.327%)
a4  = 1.10071  +/- 0.0413   (3.752%)
b1  = 0.157345 +/- 0.007552 (4.8%)
b2  = 10.106   +/- 0.1041   (1.03%)
b3  = 30.0211  +/- 0.3946   (1.314%)
b4  = 0.567116 +/- 0.01914  (3.376%)
c   = 0.0252303+/- 0.003284 (13.01%)

In 1967, Mann only calculated to sin \theta/lambda = 0, ... 1.5, and their 
tabulated values do indeed fit decently within that range, but not out to 6A.

I thought this was notworthy, and I am curious which values for these constants 
refinement programs use nowadays. Maybe George, Garib, Pavel, and Gerard may 
want to comment?

Cheers,
Tim

On 09/18/2012 10:11 AM, Oliver Einsle wrote:
> Hi there,
>
> I was just pointed to this thread and should comment on the
> discussion, as actually made the plots for this paper. James has
> clarified the issue much better than I could have, and indeed the
> calculations will fail for larger Bragg angles if you do not assume a
> reasonable B-factor (I used B=10 for the plots).
>
> Doug Rees has pointed out at the time that for large theta the c-term
> of the Cromer/Mann approximation becomes dominant, and this is where
> chaos comes in, as the Cromer/Mann parameters are only derived from a
> fit to the actual HF-calculation. They are numbers without physical
> meaning, which becomes particularly obvious if you compare the
> parameters for C and N:
>
>
> C:   2.3100  20.8439   1.0200  10.2075   1.5886  0.5687  0.8650
> 51.6512 0.2156 N:  12.2126  0.0057   3.1322  9.8933   2.0125
> 28.9975  1.1663  0.5826 -11.5290
>
> The scattering factors for these are reasonably similar, but the
> c-values are entirely different. The B-factor dampens this out and
> this is an essential point.
>
>
>
> For clarity: I made the plots using Waterloo Maple with the following
> code:
>
> restart; SF :=Matrix(17,9,readdata("scatter.dat",float,9));
>
> biso := 10; e:=  1; AFF  :=
> (e)->(SF[e,1]*exp(-SF[e,2]*s^2)+SF[e,3]*exp(-SF[e,4]*s^2)
> +SF[e,5]*exp(-SF[e,6]*s^2)+SF[e,7]*exp(-SF[e,8]*s^2)
> +SF[e,9])*exp(-biso*s^2/4);
>
> H:=  AFF(1); C:=  AFF(2); N:=  AFF(3); Ox   :=
> AFF(4); S:=  AFF(5); Fe   :=  AFF(6); Fe2  :=  AFF(7); Fe3  :=
> AFF(8); Cu   :=  AFF(9); Cu1  :=  AFF(10); Cu2  :=  AFF(11); Mo
> :=  AFF(12); Mo4  :=  AFF(13); Mo5  :=  AFF(14); Mo6  :=  AFF(15);
>
> // Plot scattering factors
>
> plot([C,N,Fe,S], s=0..1);
>
>
> // Figure 1:
>
> rho0 := (r) ->  Int((4*Pi*s^2)*Fe2*sin(2*Pi*s*r)/(2*Pi*s*r),
> s=0..1/dmax); dmax := 1.0; plot (rho0, -5..5);
>
>
> // Figure 1 (inset): Electron Density Profile
>
> rho := (r,f)
> ->(Int((4*Pi*s^2)*f*sin(2*Pi*s*r)/(2*Pi*s*r),s=0..1/dmax));
> cofactor:= 9*rho(3.3,S) + 6*rho(2.0,Fe2) + 1*rho(3.49,Mo6) +
> 1*rho(3.51,Fe3); plot(cofactor, dmax=0.5..3.5);
>
>
> The file scatter.dat is simply a collection of some form factors,
> courtesy of atomsf.lib (see attachment).
>
>
>
> Cheers,
>
> Oliver.
>
>
>
> Am 9/17/12 11:24 AM schrieb "Tim Gruene" unter
> :
>
> Dear James et al.,
>
> so to summarise, the

Re: [ccp4bb] Series termination effect calculation.

2012-09-18 Thread Tim Gruene 
-BEGIN PGP SIGNED MESSAGE-
Hash: SHA1

Hello Oliver,

when you fit the values from ICA Tab 6.1.1.1 with gnuplot, the values
of C and N become much more comparable. c(C) = 0.017 and especially
c(N) = 0.025 > 0!!!
for C:
Final set of parametersAsymptotic Standard Error
=====

a1  = 0.604126 +/- 0.02326  (3.85%)
a2  = 2.63343  +/- 0.03321  (1.261%)
a3  = 1.52123  +/- 0.03528  (2.319%)
a4  = 1.2211   +/- 0.02225  (1.822%)
b1  = 0.185807 +/- 0.00629  (3.385%)
b2  = 14.6332  +/- 0.1355   (0.9263%)
b3  = 41.6948  +/- 0.5345   (1.282%)
b4  = 0.717984 +/- 0.01251  (1.743%)
c   = 0.0171359+/- 0.002045 (11.93%)

for N:
Final set of parametersAsymptotic Standard Error
=====

a1  = 0.723788 +/- 0.04334  (5.988%)
a2  = 3.24589  +/- 0.04074  (1.255%)
a3  = 1.90049  +/- 0.04422  (2.327%)
a4  = 1.10071  +/- 0.0413   (3.752%)
b1  = 0.157345 +/- 0.007552 (4.8%)
b2  = 10.106   +/- 0.1041   (1.03%)
b3  = 30.0211  +/- 0.3946   (1.314%)
b4  = 0.567116 +/- 0.01914  (3.376%)
c   = 0.0252303+/- 0.003284 (13.01%)

In 1967, Mann only calculated to sin \theta/lambda = 0, ... 1.5, and
their tabulated values do indeed fit decently within that range, but
not out to 6A.

I thought this was notworthy, and I am curious which values for these
constants refinement programs use nowadays. Maybe George, Garib,
Pavel, and Gerard may want to comment?

Cheers,
Tim

On 09/18/2012 10:11 AM, Oliver Einsle wrote:
> Hi there,
> 
> I was just pointed to this thread and should comment on the
> discussion, as actually made the plots for this paper. James has
> clarified the issue much better than I could have, and indeed the
> calculations will fail for larger Bragg angles if you do not assume
> a reasonable B-factor (I used B=10 for the plots).
> 
> Doug Rees has pointed out at the time that for large theta the
> c-term of the Cromer/Mann approximation becomes dominant, and this
> is where chaos comes in, as the Cromer/Mann parameters are only
> derived from a fit to the actual HF-calculation. They are numbers
> without physical meaning, which becomes particularly obvious if you
> compare the parameters for C and N:
> 
> 
> C:   2.3100  20.8439   1.0200  10.2075   1.5886  0.5687  0.8650
> 51.6512 0.2156 N:  12.2126  0.0057   3.1322  9.8933   2.0125
> 28.9975  1.1663  0.5826 -11.5290
> 
> The scattering factors for these are reasonably similar, but the
> c-values are entirely different. The B-factor dampens this out and
> this is an essential point.
> 
> 
> 
> For clarity: I made the plots using Waterloo Maple with the
> following code:
> 
> restart; SF :=Matrix(17,9,readdata("scatter.dat",float,9));
> 
> biso := 10; e:=  1; AFF  :=
> (e)->(SF[e,1]*exp(-SF[e,2]*s^2)+SF[e,3]*exp(-SF[e,4]*s^2) 
> +SF[e,5]*exp(-SF[e,6]*s^2)+SF[e,7]*exp(-SF[e,8]*s^2) 
> +SF[e,9])*exp(-biso*s^2/4);
> 
> H:=  AFF(1); C:=  AFF(2); N:=  AFF(3); Ox   :=
> AFF(4); S:=  AFF(5); Fe   :=  AFF(6); Fe2  :=  AFF(7); Fe3  :=
> AFF(8); Cu   :=  AFF(9); Cu1  :=  AFF(10); Cu2  :=  AFF(11); Mo
> :=  AFF(12); Mo4  :=  AFF(13); Mo5  :=  AFF(14); Mo6  :=  AFF(15);
> 
> // Plot scattering factors
> 
> plot([C,N,Fe,S], s=0..1);
> 
> 
> // Figure 1:
> 
> rho0 := (r) ->  Int((4*Pi*s^2)*Fe2*sin(2*Pi*s*r)/(2*Pi*s*r),
> s=0..1/dmax); dmax := 1.0; plot (rho0, -5..5);
> 
> 
> // Figure 1 (inset): Electron Density Profile
> 
> rho := (r,f)
> ->(Int((4*Pi*s^2)*f*sin(2*Pi*s*r)/(2*Pi*s*r),s=0..1/dmax)); 
> cofactor:= 9*rho(3.3,S) + 6*rho(2.0,Fe2) + 1*rho(3.49,Mo6) + 
> 1*rho(3.51,Fe3); plot(cofactor, dmax=0.5..3.5);
> 
> 
> The file scatter.dat is simply a collection of some form factors,
> courtesy of atomsf.lib (see attachment).
> 
> 
> 
> Cheers,
> 
> Oliver.
> 
> 
> 
> Am 9/17/12 11:24 AM schrieb "Tim Gruene" unter
> :
> 
> 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 e

Re: [ccp4bb] Series termination effect calculation.

2012-09-18 Thread Oliver Einsle 
Hi there,

I was just pointed to this thread and should comment on the discussion, as
actually made the plots for this paper. James has clarified the issue much
better than I could have, and indeed the calculations will fail for larger
Bragg angles if you do not assume a reasonable B-factor (I used B=10 for
the plots).

Doug Rees has pointed out at the time that for large theta the c-term of
the Cromer/Mann approximation becomes dominant, and this is where chaos
comes in, as the Cromer/Mann parameters are only derived from a fit to the
actual HF-calculation. They are numbers without physical meaning, which
becomes particularly obvious if you compare the parameters for C and N:


C:   2.3100  20.8439   1.0200  10.2075   1.5886  0.5687  0.8650  51.6512
 0.2156
N:  12.2126  0.0057   3.1322  9.8933   2.0125  28.9975  1.1663  0.5826
-11.5290

The scattering factors for these are reasonably similar, but the c-values
are entirely different. The B-factor dampens this out and this is an
essential point.



For clarity: I made the plots using Waterloo Maple with the following code:

restart;
SF :=Matrix(17,9,readdata("scatter.dat",float,9));

biso := 10;
e:=  1;
AFF  := (e)->(SF[e,1]*exp(-SF[e,2]*s^2)+SF[e,3]*exp(-SF[e,4]*s^2)
+SF[e,5]*exp(-SF[e,6]*s^2)+SF[e,7]*exp(-SF[e,8]*s^2)
+SF[e,9])*exp(-biso*s^2/4);

H:=  AFF(1);
C:=  AFF(2);
N:=  AFF(3);
Ox   :=  AFF(4);
S:=  AFF(5);
Fe   :=  AFF(6);
Fe2  :=  AFF(7);
Fe3  :=  AFF(8);
Cu   :=  AFF(9);
Cu1  :=  AFF(10);
Cu2  :=  AFF(11);
Mo   :=  AFF(12);
Mo4  :=  AFF(13);
Mo5  :=  AFF(14);
Mo6  :=  AFF(15);

// Plot scattering factors
 
plot([C,N,Fe,S], s=0..1);


// Figure 1:

rho0 := (r) ->  Int((4*Pi*s^2)*Fe2*sin(2*Pi*s*r)/(2*Pi*s*r), s=0..1/dmax);
dmax := 1.0;
plot (rho0, -5..5);


// Figure 1 (inset): Electron Density Profile

rho := (r,f) ->(Int((4*Pi*s^2)*f*sin(2*Pi*s*r)/(2*Pi*s*r),s=0..1/dmax));
cofactor:= 9*rho(3.3,S) + 6*rho(2.0,Fe2) + 1*rho(3.49,Mo6) +
1*rho(3.51,Fe3);
plot(cofactor, dmax=0.5..3.5);


The file scatter.dat is simply a collection of some form factors, courtesy
of atomsf.lib (see attachment).



Cheers,

Oliver.



Am 9/17/12 11:24 AM schrieb "Tim Gruene" unter :

>-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 ho

Re: [ccp4bb] Series termination effect calculation.

2012-09-17 Thread Ian Tickle 
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  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!  Th

Re: [ccp4bb] Series termination effect calculation.

2012-09-17 Thread Tim Gruene 
-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

Re: [ccp4bb] Series termination effect calculation.

2012-09-17 Thread DUMAS Philippe (UDS) 

Le Lundi 17 Septembre 2012 08:32 CEST, James Holton  a écrit

Hello
May I add a few words after the thorough comments by James.
I lmay be easier to consider series termination in real space as follows.

The effect of series termination in 3D on rho(r) is of convoluting the exact 
rho(r) with the "approximation" of a delta function resulting from the limit in 
resolution. This "approximation" in 3D is given exactly by the function G[X] = 
3*[Sin(X) - X*Cos(X)]/X^3, where X = 2*Pi*r/d (r in Angstrom and d the 
resolution, also in Angstrom). This is the function appearing in the rotation 
function (for exactly the same reason of truncating the resolution).
If you consider that the iron atom is punctual (i.e. its Fourier transform 
would be merely constant), then the approximation resulting  from series 
termination is just given by  G[X] (apart for a scaling factor). And if you 
convolute the exact and ideal rho(r) with G[X], you will obtain the exact form 
of rho[r] affected by series termination. Note that, considering the Gaussian 
approximation of the structure factors, this would amount to convolute 
gaussians with G[X] (see James comments).
I join a figure corresponding to the simplification of a punctual iron atom. I 
only put on this figure the curves corresponding to the limits in resolution, 
1.3, 2 an 2.5 Angstrom because at a resolution of 1 Angstrom, the iron atom is 
definitely not punctual.
I used the same color codes as in Fig. 1 of the paper. One can see that the 
ripples on my approximate figure are essentially the same as in Fig. 1 of the 
paper. Of course, it cannot reproduce the features of rho[r] for r-->0 since 
the iron aton is definitely not punctual.

Practical comment. It is quite useful to consider the following  rule of thumb: 
the first minimum of G[X] appears at a distance equal to  0.92*d (d = 
resolution) and the first maximum  at 1.45*d. Therefore, if one suspects that 
series terminaiton effects might cause a spurious through, or peak, it may be 
enough to recalculate the e.d. map at different resolutions to check whether 
these features are moving or not.

Philippe Dumas

PS: it is instructive to make a comparison with the Airy function in astronomy. 
Airy calculated this function to take into account the distorsion brought by 
the limlited optical resolution of a telescope to a punctual image of a star. 
Nothing else than our problem, with an iron atom replacing a star...
Plus ça change, plus c'est la même chose.



> 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

Re: [ccp4bb] Series termination effect calculation.

2012-09-16 Thread James Holton 
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

Re: [ccp4bb] Series termination effect calculation.

2012-09-15 Thread Tim Gruene 
-BEGIN PGP SIGNED MESSAGE-
Hash: SHA1

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 
> 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

-BEGIN PGP SIGNATURE-
Version: GnuPG v1.4.12 (GNU/Linux)
Comment: Using GnuPG with Mozilla - http://enigmail.mozdev.org/

iD8DBQFQVFSBUxlJ7aRr7hoRAoPYAKDNQu84ozIz5Mn/qmRKiLxXPw/zPgCgwd75
KUHsKzaSdi9mL5kzZBeOqUI=
=mbnY
-END PGP SIGNATURE-


integral-crop.pdf
Description: Adobe PDF document

Re: [ccp4bb] Series termination effect calculation.

2012-09-14 Thread Ian Tickle 
On 14 September 2012 15:15, Tim Gruene  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

Re: [ccp4bb] Series termination effect calculation.

2012-09-14 Thread Tim Gruene 
-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.

Cheers,
Tim

On 09/14/2012 02:39 PM, Ian Tickle wrote:
> On 14 September 2012 13:05, Tim Gruene 
> wrote:
>> -BEGIN PGP SIGNED MESSAGE- Hash: SHA1
>> 
>> Dear Niu,
>> 
>> as far as I can tell, all your parameters are correct and the 
>> scattering term for f(s) you use is also correct. f(s)
>> furthermore matches very closely those tabulated in the Intl.
>> Tables C Tab. 6.1.1.1.
>> 
>> My reproduction of the mentioned formula with r=2.0 using MAXIMA
>> also shows quite a different graph.
>> 
>> The graph does not make much sense: as d_max -> 0 1/d_max ->
>> infinity and the integrand goes to infinity,
> 
> d_min surely, i.e. minimum d-spacing?
> 
>> because f(s) contains a constant positive term. Hence the
>> integral should oscilatingly approach infinity and not stabilise
>> as the upper integral limit approaches infinity.
> 
> Tim, exactly so, in fact like this:
> 
> http://journals.iucr.org/d/issues/2012/04/00/dz5235/dz5235.pdf
> 
> See eqn 2 and Fig 11(b).  Note that although rho(r) itself does
> indeed tend to zero as d_min -> 0 as expected, the volume integral
> of rho(r) (i.e. the calculated number of electrons) does not
> (unless d_min -> 0) !
> 
> Cheers
> 
> -- Ian
> 

- -- 
- --
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/

iD8DBQFQUzwQUxlJ7aRr7hoRAuQWAJ40ufwGUZrPmt2cLpXdccR2MMrboACgiOaG
CTfi2MtYiBpy5oxY1kcyFJ4=
=fIwH
-END PGP SIGNATURE-

Re: [ccp4bb] Series termination effect calculation.

2012-09-14 Thread Ian Tickle 
On 14 September 2012 13:05, Tim Gruene  wrote:
> -BEGIN PGP SIGNED MESSAGE-
> Hash: SHA1
>
> Dear Niu,
>
> as far as I can tell, all your parameters are correct and the
> scattering term for f(s) you use is also correct. f(s) furthermore
> matches very closely those tabulated in the Intl. Tables C Tab. 6.1.1.1.
>
> My reproduction of the mentioned formula with r=2.0 using MAXIMA also
> shows quite a different graph.
>
> The graph does not make much sense:
> as d_max -> 0 1/d_max -> infinity and the integrand goes to infinity,

d_min surely, i.e. minimum d-spacing?

> because f(s) contains a constant positive term. Hence the integral
> should oscilatingly approach infinity and not stabilise as the upper
> integral limit approaches infinity.

Tim, exactly so, in fact like this:

http://journals.iucr.org/d/issues/2012/04/00/dz5235/dz5235.pdf

See eqn 2 and Fig 11(b).  Note that although rho(r) itself does indeed
tend to zero as d_min -> 0 as expected, the volume integral of rho(r)
(i.e. the calculated number of electrons) does not (unless d_min -> 0)
!

Cheers

-- Ian

Re: [ccp4bb] Series termination effect calculation.

2012-09-14 Thread Tim Gruene 
-BEGIN PGP SIGNED MESSAGE-
Hash: SHA1

Dear Niu,

as far as I can tell, all your parameters are correct and the
scattering term for f(s) you use is also correct. f(s) furthermore
matches very closely those tabulated in the Intl. Tables C Tab. 6.1.1.1.

My reproduction of the mentioned formula with r=2.0 using MAXIMA also
shows quite a different graph.

The graph does not make much sense:
as d_max -> 0 1/d_max -> infinity and the integrand goes to infinity,
because f(s) contains a constant positive term. Hence the integral
should oscilatingly approach infinity and not stabilise as the upper
integral limit approaches infinity.

Best wishes,
Tim

On 09/13/2012 11:09 PM, Niu Tou wrote:
> Dear Colleagues,
> 
> I am trying to repeat a series termination effect calculation
> displayed as figure 2 in a publihsed paper 
> (http://www.ncbi.nlm.nih.gov/pubmed/12215645). Formula (1) was used
> to implement this calculation. Since f(s) is not defined in detail
> in this paper, I used formula and parameters listed in another 
> paper (http://scripts.iucr.org/cgi-bin/paper?a05896) to calculate
> it.
> 
> However, the result I got is not consistent with figure 2 of the
> first paper. I am not sure if the formulas I used are right or not.
> Or if there is any problem in the MatLab code, which I list below:
> 
> ###
> 
> clear all;clc;format compact;format long;
> 
> 
> 
> % matrix of a, b, c coefficients:
> 
> % rows: Fe, S, Fe1, Mo
> 
> % columns: A1; B1; A2; B2; A3; B3; A4; B4; C
> 
> fM = ...
> 
> [11.9185 4.87394 7.04848 0.34023 3.34326 15.9330 2.27228 79.0339
> 1.40818;...
> 
> 7.18742 1.43280 5.88671 0.02865 5.15858 22.1101 1.64403 55.4651 
> -3.87732;...
> 
> 11.9185 4.87394 7.04848 0.34023 3.34326 15.9330 2.27228 79.0339
> 1.40818;...
> 
> 19.3885 0.97877 11.8308 10.0885 3.75919 31.9738 1.46772 117.932
> 5.55047];
> 
> 
> 
> %%% store radius data:
> 
> % distance from: origin
> 
> % columns: Fe, S, Fe, Mo
> 
> R_el = [2.0 3.3 3.5 3.5];
> 
> RHO_t = zeros(4,400);
> 
> for numel = 1:4
> 
> EL = numel;
> 
> RHO = zeros(1,400);
> 
> dmax = zeros(1,400);
> 
> for iter = 1:400
> 
> dmax(iter) = iter/100; % in angstroms
> 
> % numerical integration
> 
> int_fun = @(s) 4*pi*(s.^2).* ...
> 
> (fM(EL,1).*exp(-fM(EL,2).*(s.^2)*0.25) + ...
> 
> fM(EL,3).*exp(-fM(EL,4).*(s.^2)*0.25) + ...
> 
> fM(EL,5).*exp(-fM(EL,6).*(s.^2)*0.25) + ...
> 
> fM(EL,7).*exp(-fM(EL,8).*(s.^2)*0.25) + fM(EL,9)).* ...
> 
> sin(2*pi*s*R_el(EL))./(2*pi*s*R_el(EL));
> 
> 
> 
> RHO(iter) = quad(int_fun,0,1/dmax(iter));
> 
> clc;display(iter);display(numel);
> 
> end
> 
> RHO_t(numel,:) = RHO;
> 
> end
> 
> 
> 
> RHO_t(1,:)= 6*RHO_t(1,:);
> 
> RHO_t(2,:)= 9*RHO_t(2,:);
> 
> 
> 
> figure;
> 
> axis([0.5 3.5 -10 10]); hold on;
> 
> plot(dmax,RHO_t(1,:),...
> 
> dmax,RHO_t(2,:),...
> 
> dmax,RHO_t(3,:),...
> 
> dmax,RHO_t(4,:),...
> 
> dmax,sum(RHO_t,1));
> 
> title('Electron Density Profile');
> 
> legend('Fe','S','Fe1','Mo','Sum');
> 
> xlabel('d_m_a_x'); ylabel('Rho(r)');
> 
> set(gca,'XDir','reverse');
> 
> ##
> 
> 
> 
> Any suggestions will be appreciated. Thanks!
> 
> 
> 
> Niu
> 

- -- 
- --
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/

iD8DBQFQUx2lUxlJ7aRr7hoRAi6IAJ9cU5etz4nb8Y7ti9zSjjL5P3ptJgCgtsl4
RM8Zn+l70MLPXoO1MpGKYUI=
=Lq9Z
-END PGP SIGNATURE-

Re: [ccp4bb] Series termination effect calculation.

2012-09-13 Thread Pavel Afonine 
Hi,
pointers listed here may be of help:

1) CCP4 Newsletterhttp://www.ccp4.ac.uk/newsletters/newsletter42/content.html
On the Fourier series truncation peaks at subatomic resolution
Anne Bochow, Alexandre Urzhumtsev


2) https://www.phenix-online.org/presentations/latest/pavel_maps.pdf

3) Central Ligand in the FeMo-Cofactor Nitrogenase MoFe-Protein at 1.16 Å
Resolution: A.

Oliver Einsle, et al. Science, 1696 (2002) 297

4) Page 267 Figure 4:

On the possibility of the observation of valence electron density for
individual
bonds in proteins in conventional difference maps

P. V. Afonine, V. Y. Lunin, N. Muzet and A. Urzhumtsev
Acta Cryst. (2004). D60, 260-274


Pavel

On Thu, Sep 13, 2012 at 2:09 PM, Niu Tou  wrote:

> Dear Colleagues,
>
> I am trying to repeat a series termination effect calculation displayed as
> figure 2 in a publihsed paper (http://www.ncbi.nlm.nih.gov/pubmed/12215645). 
> Formula
> (1) was used to implement this calculation. Since f(s) is not defined in
> detail in this paper, I used formula and parameters listed in another
> paper (http://scripts.iucr.org/cgi-bin/paper?a05896) to calculate it.
>
> However, the result I got is not consistent with figure 2 of the first
> paper. I am not sure if the formulas I used are right or not. Or if there
> is any problem in the MatLab code, which I list below:
>
> ###
>
> clear all;clc;format compact;format long;
>
>
>
> % matrix of a, b, c coefficients:
>
> % rows: Fe, S, Fe1, Mo
>
> % columns: A1; B1; A2; B2; A3; B3; A4; B4; C
>
> fM = ...
>
> [11.9185 4.87394 7.04848 0.34023 3.34326 15.9330 2.27228 79.0339
> 1.40818;...
>
>  7.18742 1.43280 5.88671 0.02865 5.15858 22.1101 1.64403 55.4651
> -3.87732;...
>
>  11.9185 4.87394 7.04848 0.34023 3.34326 15.9330 2.27228 79.0339
> 1.40818;...
>
>  19.3885 0.97877 11.8308 10.0885 3.75919 31.9738 1.46772 117.932 5.55047];
>
>
>
> %%% store radius data:
>
> % distance from: origin
>
> % columns: Fe, S, Fe, Mo
>
> R_el = [2.0 3.3 3.5 3.5];
>
> RHO_t = zeros(4,400);
>
>  for numel = 1:4
>
>  EL = numel;
>
>  RHO = zeros(1,400);
>
>  dmax = zeros(1,400);
>
>  for iter = 1:400
>
> dmax(iter) = iter/100; % in angstroms
>
> % numerical integration
>
>  int_fun = @(s) 4*pi*(s.^2).* ...
>
> (fM(EL,1).*exp(-fM(EL,2).*(s.^2)*0.25) + ...
>
>  fM(EL,3).*exp(-fM(EL,4).*(s.^2)*0.25) + ...
>
>  fM(EL,5).*exp(-fM(EL,6).*(s.^2)*0.25) + ...
>
>  fM(EL,7).*exp(-fM(EL,8).*(s.^2)*0.25) + fM(EL,9)).* ...
>
>  sin(2*pi*s*R_el(EL))./(2*pi*s*R_el(EL));
>
>
>
>  RHO(iter) = quad(int_fun,0,1/dmax(iter));
>
> clc;display(iter);display(numel);
>
>  end
>
>  RHO_t(numel,:) = RHO;
>
>  end
>
>
>
> RHO_t(1,:)= 6*RHO_t(1,:);
>
> RHO_t(2,:)= 9*RHO_t(2,:);
>
>
>
>  figure;
>
>  axis([0.5 3.5 -10 10]); hold on;
>
>  plot(dmax,RHO_t(1,:),...
>
>   dmax,RHO_t(2,:),...
>
>   dmax,RHO_t(3,:),...
>
>   dmax,RHO_t(4,:),...
>
>   dmax,sum(RHO_t,1));
>
>   title('Electron Density Profile');
>
>   legend('Fe','S','Fe1','Mo','Sum');
>
>   xlabel('d_m_a_x'); ylabel('Rho(r)');
>
>   set(gca,'XDir','reverse');
>
> ##
>
>
>
> Any suggestions will be appreciated. Thanks!
>
>
>
> Niu
>

[ccp4bb] Series termination effect calculation.

2012-09-13 Thread Niu Tou 
Dear Colleagues,

I am trying to repeat a series termination effect calculation displayed as
figure 2 in a publihsed paper
(http://www.ncbi.nlm.nih.gov/pubmed/12215645). Formula
(1) was used to implement this calculation. Since f(s) is not defined in
detail in this paper, I used formula and parameters listed in another
paper (http://scripts.iucr.org/cgi-bin/paper?a05896) to calculate it.

However, the result I got is not consistent with figure 2 of the first
paper. I am not sure if the formulas I used are right or not. Or if there
is any problem in the MatLab code, which I list below:

###

clear all;clc;format compact;format long;



% matrix of a, b, c coefficients:

% rows: Fe, S, Fe1, Mo

% columns: A1; B1; A2; B2; A3; B3; A4; B4; C

fM = ...

[11.9185 4.87394 7.04848 0.34023 3.34326 15.9330 2.27228 79.0339 1.40818;...

 7.18742 1.43280 5.88671 0.02865 5.15858 22.1101 1.64403 55.4651
-3.87732;...

 11.9185 4.87394 7.04848 0.34023 3.34326 15.9330 2.27228 79.0339 1.40818;...

 19.3885 0.97877 11.8308 10.0885 3.75919 31.9738 1.46772 117.932 5.55047];



%%% store radius data:

% distance from: origin

% columns: Fe, S, Fe, Mo

R_el = [2.0 3.3 3.5 3.5];

RHO_t = zeros(4,400);

 for numel = 1:4

 EL = numel;

 RHO = zeros(1,400);

 dmax = zeros(1,400);

 for iter = 1:400

dmax(iter) = iter/100; % in angstroms

% numerical integration

 int_fun = @(s) 4*pi*(s.^2).* ...

(fM(EL,1).*exp(-fM(EL,2).*(s.^2)*0.25) + ...

 fM(EL,3).*exp(-fM(EL,4).*(s.^2)*0.25) + ...

 fM(EL,5).*exp(-fM(EL,6).*(s.^2)*0.25) + ...

 fM(EL,7).*exp(-fM(EL,8).*(s.^2)*0.25) + fM(EL,9)).* ...

 sin(2*pi*s*R_el(EL))./(2*pi*s*R_el(EL));



 RHO(iter) = quad(int_fun,0,1/dmax(iter));

clc;display(iter);display(numel);

 end

 RHO_t(numel,:) = RHO;

 end



RHO_t(1,:)= 6*RHO_t(1,:);

RHO_t(2,:)= 9*RHO_t(2,:);



 figure;

 axis([0.5 3.5 -10 10]); hold on;

 plot(dmax,RHO_t(1,:),...

  dmax,RHO_t(2,:),...

  dmax,RHO_t(3,:),...

  dmax,RHO_t(4,:),...

  dmax,sum(RHO_t,1));

  title('Electron Density Profile');

  legend('Fe','S','Fe1','Mo','Sum');

  xlabel('d_m_a_x'); ylabel('Rho(r)');

  set(gca,'XDir','reverse');

##



Any suggestions will be appreciated. Thanks!



Niu