{+ file: bsharp.inp +} {+ directory: xtal_util +} {+ description: determine optimum value of Bsharp. +} {+ comment: determines optimum value for B-sharpening: it determines the smallest Bsharp value such that ln(/) > 0 for all values where Fsharp = exp ( - Bsharp * S^2 / 4 ) +} {+ authors: Axel T. Brunger, Byron DeLaBarre +} {+ copyright: Yale University +} {+ reference: Wilson, A. J. C. The probability distibution of X-ray intensities, Acta. Cryst. (1949) vol 2, 318 +} {+ reference: B.DeLaBarre and A.T. Brunger, Nucleotide dependent motion and mechanism of action of p97/VCP, J. Mol. Biol. (2005) vol 347, 437 +} {- Guidelines for using this file: - all strings must be quoted by double-quotes - logical variables (true/false) are not quoted - do not remove any evaluate statements from the file -} {- begin block parameter definition -} define( {====================== crystallographic data ========================} {* space group *} {* use International Table conventions with subscripts substituted by parenthesis *} {===>} sg="P4(1)2(1)2"; {* unit cell parameters in Angstroms and degrees *} {+ table: rows=1 "cell" cols=6 "a" "b" "c" "alpha" "beta" "gamma" +} {===>} a=101.4; {===>} b=101.4; {===>} c=199.5; {===>} alpha=90; {===>} beta=90; {===>} gamma=90; {* reflection file *} {===>} reflection_infile="eg1_abcd.cv"; {* specify amplitude array to be analysed *} {===>} f="fobs"; {==================== ln(/) plot options ==================} {* resolution range *} {+ table: rows=1 "resolution" cols=2 "lowest" "highest" +} {===>} low_res=500.0; {===>} high_res=4.0; {* reflection selection in addition to resolution selection *} {===>} ref_sele=(all); {* number of resolution shells for analysis *} {===>} bins=20; {* approximate number of amino acid residues in the asymmetric unit *} {===>} aa_residues=740; {* give any additional atoms of the specified chemical type and the number of these atoms in the asymmetric unit *} {+ table: rows=9 cols=2 "chemical type" "number of atoms in ASU" +} {===>} atom_type_1="C"; {===>} atom_number_1=0; {===>} atom_type_2="N"; {===>} atom_number_2=0; {===>} atom_type_3="O"; {===>} atom_number_3=0; {===>} atom_type_4="S"; {===>} atom_number_4=10; {===>} atom_type_5=""; {===>} atom_number_5=0; {===>} atom_type_6=""; {===>} atom_number_6=0; {===>} atom_type_7=""; {===>} atom_number_7=0; {===>} atom_type_8=""; {===>} atom_number_8=0; {===>} atom_type_9=""; {===>} atom_number_9=0; {* step in B-factor sharpening value *} {===>} decrement=-2.5; {=========================== output files ============================} {* output file *} {+ list: output files will be: +} {===>} output="bsharp.list"; {===========================================================================} { things below this line do not normally need to be changed } {===========================================================================} ) {- end block parameter definition -} checkversion 1.2 evaluate ($log_level=quiet) xray a=&a b=&b c=&c alpha=&alpha beta=&beta gamma=&gamma @CNS_XTALLIB:spacegroup.lib (sg=&sg;) if ( &BLANK%reflection_infile = false ) then reflection @@&reflection_infile end end if set echo=off end if ( &BLANK%f = false ) then query name=&STRIP%f domain=reciprocal end if ( $object_exist = false ) then display display ************************************************** display Error : data set amplitude array &f does not exist display ************************************************** display abort end if {- this object can be of any type -} end if set echo=on end binresolution &low_res &high_res bins=&bins declare name=diff domain=reciprocal type=real end end topology {- create a dummy residue for isolated point atoms -} residue scat atom=scat type scat mass 10. end end end if (&aa_residues>0) then {- approximate number of atoms: 5 C + 3 N + 1 O -} evaluate ($atom_number=&aa_residues * 5) {- BEGIN MODIFICATION -} if ($atom_number > 9999999 ) then display ERROR: number of residues too large (resulting in more than 999,9999 atoms) abort else evaluate ($counter=0) evaluate ($atom_nn=$atom_number) while ($atom_nn > 0) loop aac evaluate ($atom_ii=min(9999,$atom_nn)) evaluate ($atom_nn=$atom_nn-$atom_ii) evaluate ($counter=$counter+1) evaluate ($seg="C" + encode ( $counter )) segment name= $seg mole number=$atom_ii name=scat end end do (chemical="C") ( segid $seg ) end loop aac end if evaluate ($atom_number=&aa_residues * 3) if ($atom_number > 9999999 ) then display ERROR: number of residues too large (resulting in more than 999,9999 atoms) abort else evaluate ($counter=0) evaluate ($atom_nn=$atom_number) while ($atom_nn > 0) loop aan evaluate ($atom_ii=min(9999,$atom_nn)) evaluate ($atom_nn=$atom_nn-$atom_ii) evaluate ($counter=$counter+1) evaluate ($seg="N" + encode ( $counter )) segment name=$seg mole number=$atom_ii name=scat end end do (chemical="N") ( segid $seg ) end loop aan end if evaluate ($atom_number=&aa_residues ) if ($atom_number > 9999999 ) then display ERROR: number of residues too large (resulting in more than 999,9999 atoms) abort else evaluate ($counter=0) evaluate ($atom_nn=$atom_number) while ($atom_nn > 0) loop aao evaluate ($atom_ii=min(9999,$atom_nn)) evaluate ($atom_nn=$atom_nn-$atom_ii) evaluate ($counter=$counter+1) evaluate ($seg="O" + encode ( $counter )) segment name=$seg mole number=$atom_ii name=scat end end do (chemical="O") ( segid $seg ) end loop aao end if end if evaluate ($ii=1) while (&exist_&atom_type_$ii=true) loop atm if ( &BLANK%atom_type_$ii = false ) then if (&atom_number_$ii > 9999) then display ERROR: number of atoms too large (> 9999) for &atom_type_$ii abort elseif (&atom_number_$ii > 0) then evaluate ($atom_type_$ii=capitalize(&atom_type_$ii)) segment name=$atom_type_$ii mole number=&atom_number_$ii name=scat end end do (chemical=$atom_type_$ii) ( segid $atom_type_$ii ) end if end if evaluate ($ii=$ii+1) end loop atm {- END MODIFICATION -} {- coordinates are not required. However, we must initialize the coordinates in order to avoid warning messages -} do (x=0.) ( all ) do (y=0.) ( all ) do (z=0.) ( all ) xray @CNS_XRAYLIB:scatter.lib show sum ( amplitude(&STRIP%f) ) ( amplitude(&STRIP%f)>0 and &low_res >= d >= &high_res and &ref_sele ) if ($select=0) then display display ************************************************* display Error: no reflections selected display in overall analysis range display ( &low_res >= d >= &high_res ) display and overall selection ( &ref_sele ) display display ************************************************** abort end if {- determine minimum B-factor applied to amplitudes such that ln(/) becomes positive everywhere -} {- 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 declare name=ffsum domain=reci type=real end declare name=select domain=reci type=integer end do ( select=0 ) ( all ) do ( select=1 ) ( amplitude(&STRIP%f)>0 and &low_res >= d >= &high_res and &ref_sele) {- compute f^2 i.e. the sum of squared atomic scattering forms -} predict mode=ff2 to=ffsum selection=( select=1 ) atomselection=( all ) end {- compute average over resolution bins -} do (ffsum=ffsum*epsilon) (select=1) do (wilson_OB=SAVE(ffsum)) (select=1) evaluate ($bsharp=10) evaluate ($min=-1) evaluate ($trial=0) while ($min < 0 ) loop main evaluate ($bsharp=$bsharp+&decrement) do (wilson_OC=exp( - $bsharp * s^2/4) * amplitude(&STRIP%f)) (select=1) {- compute average over resolution bins -} do (wilson_OC=SAVE(wilson_OC^2)) (select=1) {- compute the ratio of log(/) -} do (wilson_OC=log(max(0.00001,wilson_OC/wilson_OB))) (select=1) show min ( wilson_OC ) (select=1) evaluate ($min=$result) evaluate ($trial=$trial+1) if ($trial > 10000) then display maximum number of trials exceeded display check data display check decrement value abort end if end loop main set display=&output end display optimum value for Bsharp = $bsharp display display display plot of ln(/) for Fo = exp ( $bsharp s^2 / 4 ) * &STRIP%f display 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: ln(/) statistics output=&output ( s^2 ) ( wilson_OC ) selection= (select=1) end close &output end {- 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 end stop