! Module file: wilsonplot ! ! CNS module ! ********** ! ! Authors: Jian-Sheng Jiang and Axel T. Brunger ! ! copyright Yale University ! ! Function: ! Produces a Wilson plot (message=LONG) ! Scale factor and B factor are returned in parameters &KW, &BW ! ! ! References: ! 1. Wilson, A.J.C. (1949). Acta Cryst. 2, 318-321. ! "The Probability Distribution of X-ray Intensities". ! 2. Rogers, D. (1965). In "Computing Methods in Crystallgraphy" ! (Rollett, J.S., ed.), 117-132. ! "Statistical Properties of Reciprocal Space". ! 3. Main P. (1975). In "Crystallographic Computing Techniques" ! (Ahmed, F.R., Huml, K. and Sedlacek, B. eds.), 98-103. ! "Recent Developments in the MULTAN System - ! The Use of Molecular Structure". ! ! Application to the data set: ! Fobs = Fobs * KW^(1/2)*exp(Bw/4*s^2) ! module { wilsonplot } ( &bins=8; { number of bins for average } &f1=fobs; { data set } &sel=(all); { reflection selection for plot } &wilson_low_res=10000.0; { resolution range for linear fit } &wilson_high_res=0.5; { resolution range for linear fit } &disp=OUTPUT; { display output filename } &mess=LONG; { SHORT | LONG display messages } &KW=$KW; { Output: Wilson scale factor } &BW=$BW; { Output: Wilson B factor } &atom_sel=( all ) ; { atom selection. } ) set message ? end evaluate ($message_old=$result) set echo ? end evaluate ($echo_old=$result) if ( $log_level = verbose ) then set echo=on message=normal end else set echo=off message=off end end if checkversion 1.3 evaluate ($old_bins=$bin_number) evaluate ($old_bin_res_low=$bin_resolution_low) evaluate ($old_bin_res_high=$bin_resolution_high) bins=&bins show sum (1) (&sel) if ($result>0) then show max (d) (&sel) evaluate ($Wilson_low_res=$result+0.001) show min (d) (&sel) evaluate ($Wilson_high_res=$result-0.001) if ($select=0) then display display ******************************************************* display Error: no reflections selected for Wilson plot analysis display ******************************************************* display abort end if binresolution $wilson_low_res $wilson_high_res {* Declare local arrays *} declare domain=reci type=real name=wilson_oa end declare domain=reci type=real name=wilson_ob end declare domain=reci type=real name=wilson_oc end if (&disp = " ") then evaluate (&disp="wilsonplot.list") end if {* compute f^2 i.e. the sum of squared atomic scattering forms *} declare name=ffsum domain=reci type=real end predict mode=ff2 to=ffsum selection=( &sel ) atomselection=( &atom_sel ) end {* compute average over resolution bins *} do (ffsum=ffsum*epsilon) ( &sel ) do (wilson_OB=SAVE(ffsum)) ( &sel ) {* compute average over resolution bins *} do (wilson_OC=SAVE(amplitude(&f1)^2)) ( &sel ) {* compute the ratio of / *} do (wilson_OB=wilson_OB/wilson_OC) ( &sel ) {* divide the resolution range with the equal s^2 bins *} {* the first and the last resolution bins are not fit *} evaluate ($Wilson_low_s=1.0/min(&Wilson_low_res,$wilson_low_res)^2) evaluate ($Wilson_high_s=1.0/max(&Wilson_high_res,$wilson_high_res)^2) evaluate ($Wilson_delta=($Wilson_high_s-$Wilson_low_s)/&bins) evaluate ($Wilson_low_cut=1.0/($Wilson_low_s+$Wilson_delta)^(1/2)) evaluate ($Wilson_high_cut=1.0/($Wilson_high_s-$Wilson_delta)^(1/2)) {* determine K scale and B value using the SHAPE function *} {* SHAPE(Y,X) = Y/X = K*exp(-B*(sin(theta)/lambda)^2) *} do (wilson_OA=1.0) ( &sel ) show max (d) (&sel and $Wilson_low_cut >= d >= $Wilson_high_cut) evaluate ($Wilson_low_res=$result) show min (d) (&sel and $Wilson_low_cut >= d >= $Wilson_high_cut) evaluate ($Wilson_high_res=$result) if ($select=0) then display display ******************************************************* display Error: no reflections available in selected Wilson display fit resolution range ( &wilson_low_res > d > & wilson_high-res ) display Suggestion: change Wilson resolution range or turn off display the Wilson analysis and/or scaling display ******************************************************* display abort end if bins=1 do (wilson_OC=SHAPE(wilson_OA,wilson_OB)) ( $Wilson_low_res >= d >= $Wilson_high_res and &sel ) show max (d) (&sel) evaluate ($Wilson_low_res=$result+0.001) show min (d) (&sel) evaluate ($Wilson_high_res=$result-0.001) binresolution $wilson_low_res $wilson_high_res bins=&bins evaluate ($Wilson_lnk=log($K_BIN1)) evaluate ($Wilson_KW=1.0/$K_BIN1) evaluate ($Wilson_BW=$B_BIN1/2.0) evaluate (&KW=$Wilson_KW) evaluate (&BW=$Wilson_BW) {* produce output *} set display=&disp end if (&mess = LONG) then { estimate a scale for plot data output - to the scale of 10^4 } show sum ( amplitude(&f1)^2 ) ( &sel ) evaluate ($plot_scale=$result/&bins) show sum ( ffsum ) ( &sel ) evaluate ($plot_scale=max($plot_scale,$result/&bins)) evaluate ($power_scale=nint(log($plot_scale)/log(10))-4) evaluate ($plot_scale=10^$power_scale) evaluate ($sscale=10^(-$power_scale)) display column 1: bin number display columns 2 & 3: resolution range display column 4: number of reflections in bin display column 5: <|s|^2> display column 6: <|Fobs|^2> * $sscale[F12.6] display column 7: * $sscale[F12.6] display column 8: ln(/) display column 9: linear fit to ln(/) vs. <|s|^2> for \ $Wilson_low_cut[F6.2] > d > $Wilson_high_cut[F6.2] A statistics output=&disp ( s^2 ) ( amplitude(&f1)^2/$plot_scale) ( ffsum/$plot_scale ) ( log(1.0/wilson_OB) ) ( $Wilson_lnk - $Wilson_BW*s^2/2.0 ) selection= ( &sel ) end end if display display Linear fit of ln(/) vs. in resol. range \ $Wilson_low_cut[F6.2] > d > $Wilson_high_cut[F6.2] A display produces / = (1/k)*exp(-B/2*s^2) \ with k= $Wilson_KW B= $Wilson_BW . display display In order to scale &f1 use the following expression &f1 = &f1 * k^(1/2)*exp( B /4*s^2). {* undeclare local arrays. *} undeclare domain=reci name=wilson_oa end undeclare domain=reci name=wilson_ob end undeclare domain=reci name=wilson_oc end undeclare domain=reci name=ffsum end else display %WILSONPLOT-ERROR: Zero reflections selected. end if bins=$old_bins binresolution $old_bin_res_low $old_bin_res_high set message=$message_old echo=$echo_old end