Calculation of "sigma_A"

[Jiang Jiansheng]

> Dear Colleagues,
>
> I need to calculate "sigma_A" as a function of resolution (Read RJ. Acta Cryst.
> (1986) A42 140-149). This plot provides a better estimate of the overall error
> in the coordinates than Luzzati plot.
> I will really appreciate it if someone has done it before could send me a
> sample of his script file, or give me a hint on how this could be easily done.
>
> Thank you in advance
> Best wishes,
> 		Igor.


There is a script file for estimating the coordinate error by using the SIGMAA
plot in X-PLOR, but  it runs on new version of X-PLOR (later than 3.1).
You may ask for a pre-release X-PLOR version directly from Axel.

Just to give you some ideas how this could be done by X-PLOR,
I am including the files following this message.

Jiansheng Jiang.

--------------cut here--------------------
! Macro: sigmaacoorderr.mac
!
! version    5-JAN-95 JSJ
!
! Function:
! 1. Estimate coordinates error by using SIGMAA,
!    return the estimated error in "#ESDERR_SIGMAA".
! 2. Produce the list file for the SIGMAA plot.
!
! Requirements:
! 1. This macro is used within xrefin.
!
! Warnings:
! Data lower than 5 angstroms resolution are not appropriate
! for estimating coordinate error, unless the solvent refinement
! had been done properly.
!
! References:
! 1. Read R.J. (1986). Acta Cryst. A42, 140-149.
!    "Improved Fourier coefficients for maps using phases from
!    partial structures with errors."
! 

macro { sigmaacoorderr }
(
   #bins=8;               { number of bins for output     }
   #f1=fobs;              { observed data                 }
   #f2=fcalc;             { calculated data               }
   #low=5.0;              { low resolution limit          }
   #high=2.0;             { high resolution limit         }
   #sel=(ampl(#f1)>0 and #low >= d >= #high);
                          { reflections selection         }
   #disp=OUTPUT;          { display output filename       }
   #mess=LONG;            { SHORT | LONG display messages }

   #ESDERR_SIGMAA=$ESDERR_SIGMAA;  { Output: symbol for coord. err. }
)


{* Declare local arrays *}
declare domain=reci type=real name=sigmaa_oa end
declare domain=reci type=real name=sigmaa_ob end
declare domain=reci type=real name=sigmaa_oc end

evaluate ($SigMacro_nbins=#bins)
evaluate ($SigMacro_low_res=#low)
evaluate ($SigMacro_high_res=#high)
evaluate ($SigMacro_display=#disp)
if ($SigMacro_display = " ") then
  evaluate ($SigMacro_display="sigmaacoorderr.list")
end if
evaluate ($SigMacro_mess=#mess)

{* check the low resolution and the number of bins *}
if ($SigMacro_low_res > 5.0) then
  display %ESDERR-WRN: below 5A resolution are not appropriate for estimating
  display %ESDERR-WRN: coord. error, unless solvent refinement has been done properly.
end if
if ($SigMacro_nbins < 4) then
  display %ESDERR-WRN: at least 4 bins for estimating coord. error.
  evaluate ($SigMacro_nbins=4)
end if

{* determine the coord. error via by the SHAPE function     *}
{* compute SIGMAA -  Read (1986), Acta Cryst. A42, 140-149. *}
{* SIGMAA = (sig_p/sig_n)^(1/2)*                            *}
{*          exp(-pi^3*<|delta_r|>^2*(sin(theta)/lambda)^2)  *}
{* where <|delta_r|>^2 is the mean square coordinate error  *}

{* compute normalized structure factor moduli *}
do (sigmaa_OA=NORM[bins=$SigMacro_nbins](amplitude(#f1))) ( #sel )
do (sigmaa_OB=NORM[bins=$SigMacro_nbins](amplitude(#f2))) ( #sel )
do (sigmaa_OC=SIGA[bins=$SigMacro_nbins](sigmaa_OA,sigmaa_OB)) ( #sel )


{* determine K scale and B value by the SHAPE function      *}
{* SHAPE(Y,X) = Y/X = K*exp(-B*(sin(theta)/lambda)^2)       *}
do (sigmaa_OA=1.0) ( #sel )
do (sigmaa_OB=SHAPE[bins=1](sigmaa_OC,sigmaa_OA)) ( #sel )
evaluate ($SigMacro_K_SIGMAA=$K_BIN1)
evaluate ($SigMacro_B_SIGMAA=$B_BIN1)

{* thus, the coord. error can be obtained from B value by   *}
{*   <|delta_r|> = sqrt(B/pi^3)                             *}
evaluate (#ESDERR_SIGMAA=sign($SigMacro_B_SIGMAA)*
                         sqrt(abs($SigMacro_B_SIGMAA)/($pi)^3))

if ($SigMacro_mess = LONG) then

  set display=$SigMacro_display end

  evaluate ($SigMacro_sigmaa=#ESDERR_SIGMAA)
  display  K= $SigMacro_K_SIGMAA  B= $SigMacro_B_SIGMAA  ESDERR_SIGMAA=  
$SigMacro_sigmaa[f8.4]

{* produce a list file for the SIGMMA plot *}

  mbins=$SigMacro_nbins
  statistic
    output=$SigMacro_display
    ( s()^2/4.0 )
    ( log(sigmaa_OC) )
    ( log($SigMacro_K_SIGMAA)-$SigMacro_B_SIGMAA*s()^2/4.0 )
    selection=( #sel )
  end


end if