{+ 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(<Fsharp^2>/<f^2>) > 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(<Fo^2>/<f^2>) 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(<Fo^2>/<f^2>) 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 <f^2> 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 <FOBS^2> over resolution bins -}
      do (wilson_OC=SAVE(wilson_OC^2)) (select=1)
      
      {- compute the ratio of log(<FOBS^2>/<f^2>) -} 
      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(<Fo^2>/<f^2>) 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(<Fo^2>/<f^2>)
   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