COMCIFS / cif_core

The IUCr CIF core dictionary
14 stars 9 forks source link

Bugs in _refln.form_factor_table dREL #374

Closed rowlesmr closed 1 year ago

rowlesmr commented 1 year ago

For the scattering factor calculations, the dREL says:

    With r  as  refln

     table = Table ()
     s  =    r.sin_theta_over_lambda

    Loop t  as  atom_type  {

       If (_diffrn_radiation.probe == 'neutron') {

            f  =   t.length_neutron

       } Else If (s < 2.0 ) {

            c  =   t.Cromer_Mann_coeffs

            f  =  (c[0] + c[1] * Exp (-c[2] * s*s)
                        + c[3] * Exp (-c[4] * s*s)
                        + c[5] * Exp (-c[6] * s*s)
                        + c[7] * Exp (-c[8] * s*s))
       } Else {

            c  =   t.hi_ang_Fox_coeffs

            f  =   Exp ( c[0] + c[1]*s + c[2]*s*s + c[3]*s*s*s )
       }
        table [ t.symbol ]  =   f
    }
       _refln.form_factor_table  =   table

From my minimal understaind of dREL, I think there are two things wrong, and the altered lines should be:

#...
    Loop t  as  atom_type_scat  {
#...
            f  =   Exp ( c[0] + c[1]*s + (c[2]/10)*s*s + (c[3]/100)*s*s*s )
#...

The first means that the loop never happens, as _atom_type.length_neutron, _atom_type.Cromer_Mann_coeffs, and _atom_type.hi_ang_Fox_coeffs don't exist

The second is kind of a big deal, as it means that the the scattering factors are calculated incorrectly. Should the division take place here, or when _atom_type_scat.hi_ang_Fox_coeffs is populated?

jamesrhester commented 1 year ago

The first is not actually a bug, because atom_type_scat is a child category of atom_type, so the object parts of its datanames (the part after the period) can be accessed using the parent category name. Likewise, the contents of both categories are effectively joined on the category keys from the point of view of dREL, so the number of rows would be the number of distinct category key combinations, in this case, the number of atomic types.

I'm not so sure on the second point. The values of the hi_ang_Fox_coeffs that were used in testing are stored in templ_enum.cif if you have a handy source of them to see if they've been pre-divided. I do remember when testing my dREL interpreter using this calculation that the final F_calcs were ever so slightly different to an alternative source of F_calc, but I didn't chase that down.

rowlesmr commented 1 year ago

_atom_type_scat.hi_ang_fox_coeffs is defined as:

save_atom_type_scat.hi_ang_fox_coeffs

    _definition.id                '_atom_type_scat.hi_ang_Fox_coeffs'
    _definition.update            2012-11-30
    _description.text
;
    The set of Fox et al. coefficients for generating high angle
    X-ray scattering factors. [ c0, c1, c2, c3 ]
    Ref: International Tables for Crystallography, Vol. C
         (1991) Table 6.1.1.5
;
    _name.category_id             atom_type_scat
    _name.object_id               hi_ang_Fox_coeffs
    _type.purpose                 Number
    _type.source                  Derived
    _type.container               List
    _type.dimension               '[4]'
    _type.contents                Real
    _units.code                   none
    _method.purpose               Evaluation
    _method.expression
;
     With t  as  atom_type_scat

    _atom_type_scat.hi_ang_Fox_coeffs  =
    [t.hi_ang_Fox_c0,t.hi_ang_Fox_c1,t.hi_ang_Fox_c2,t.hi_ang_Fox_c3]
;

save_

so you're just putting the values _atom_type_scat.hi_ang_Fox_c0-3 into a list.

Picking on _atom_type_scat.hi_ang_Fox_c3, it is defined as

save_atom_type_scat.hi_ang_fox_c3

    _definition.id                '_atom_type_scat.hi_ang_Fox_c3'
    _name.category_id             atom_type_scat
    _name.object_id               hi_ang_Fox_c3

    _import.get
        [
         {'file':templ_attr.cif  'save':hi_ang_Fox_coeffs}
         {'file':templ_enum.cif  'save':hi_ang_Fox_c3}
        ]

save_

which, when you do the imports, becomes

save_atom_type_scat.hi_ang_fox_c3

    _definition.id                '_atom_type_scat.hi_ang_Fox_c3'
    _name.category_id             atom_type_scat
    _name.object_id               hi_ang_Fox_c3

    _definition.update           2012-11-29
    _description.text
;
    The set of data items used to define Fox et al.  coefficients
     for generation of high angle (s >2.0) X-ray scattering factors.

        Ref: International Tables for Crystallography, Vol. C
             (1991) Table 6.1.1.5
;
    _type.purpose                Number
    _type.source                 Assigned
    _type.container              Single
    _type.contents               Real
    _enumeration.def_index_id  '_atom_type.symbol'
    _units.code                  none

     loop_        
         _enumeration_default.index 
         _enumeration_default.value
     H     .0      D     .0      H1-   .0      He    -2.5476 Li    -.71949 
     Li1+  -.71949 Be    -0.04297 Be2+  -0.04297 B     -0.65981 C     -0.42715
     #...

save_

ITC (2004, 3rd ed.) has this table

image

where you can see that the values for a_3 are identical, but the column header says that these values are 100x their actual value. If you put them as-is into the equation

image

you get the wrong answer.

Carbon
a0   1.70560
a1   -1.56760
a2   1.18930
a3   -0.42715

s = 2

SF_C = exp(1.70560 + (-1.56760* 2) + 1.18930*4 + (-0.42715 * 8) = exp(-0.0896) = 0.914

but 

SF_C = exp(1.70560 + (-1.56760* 2) + (1.18930/10)*4 + ((-0.42715/100) * 8) = exp(-0.98805) = 0.372

From earlier in the ITC

image

Carbon, at s = 2, has a scattering factor of 0.373.

jamesrhester commented 1 year ago

I think that's pretty convincing proof that this is indeed a bug.

rowlesmr commented 1 year ago

But where is it? I think that, strictly, the values in templ_enum.cif are incorrect.