{+ file: hlcoeff_blur.inp +}
{+ directory: xtal_util +}
{+ description: "Blur" Hendrickson-Lattman coefficients for use
                 in refinement with the MLHL target +}
{+ comment:
            The phase probability distribution is "blurred"
            by application of a scale factor (S) and a B-factor (B):
              HLA_new = S * e^(-B*s*s) * HLA_old
              HLB_new = S * e^(-B*s*s) * HLB_old
              HLC_new = S * e^(-B*s*s) * HLC_old
              HLD_new = S * e^(-B*s*s) * HLD_old
            This is performed to compensate for overestimation of phase
            accuracy which most often occurs after density modification
            or when probability distributions are derived from an
            atomic model.
            Warning: this should only be used when there is a good
                     reason to believe the phases are biased. MAD or MIR
                     phase probability distributions should usually
                     not be modified in this way. +}
{+ authors: Axel T. Brunger and Paul D. Adams +}
{+ copyright: Yale University +}

{+ reference: N.S. Pannu, G.N. Murshudov, E.J. Dodson, and R.J. Read.
              Incorporation of prior phase information strengthens
              maximum likelihood structural refinement,
              Acta Cryst. D54, in press (1998) +}

{- 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="P2(1)2(1)2(1)";

{* unit cell parameters in Angstroms and degrees *}
{+ table: rows=1 "cell" cols=6 "a" "b" "c" "alpha" "beta" "gamma" +}
{===>} a=61.76;
{===>} b=40.73;
{===>} c=26.74;
{===>} alpha=90;
{===>} beta=90;
{===>} gamma=90;

{* reflection file(s) *}
{* if anomalous and non-anomalous data are to be combined
   read non-anomalous data first *}
{===>} reflection_infile_1="amy_mir.hkl";
{===>} reflection_infile_2="";
{===>} reflection_infile_3="";
{===>} reflection_infile_4="";
{===>} reflection_infile_5="";

{============================ input arrays ===========================}

{* input Hendrickson-Lattman coefficients *}
{+ table: rows=1 "HL coefficients" cols=4 "A" "B" "C" "D" +}
{===>} pa_in="pa";
{===>} pb_in="pb";
{===>} pc_in="pc";
{===>} pd_in="pd";

{=========================== output arrays ===========================}

{* output Hendrickson-Lattman coefficients *}
{+ table: rows=1 "HL coefficients" cols=4 "A" "B" "C" "D" +}
{===>} pa_out="pa_blur";
{===>} pb_out="pb_blur";
{===>} pc_out="pc_blur";
{===>} pd_out="pd_blur";

{========================= blurring parameters =======================}

{* Blurring scale factor *}
{===>} blur_scale=0.75;

{* Blurring B-factor *}
{===>} blur_bfactor=30;

{============================ output files ===========================}

{* output reflection file *}
{===>} reflection_outfile="hlcoeff_blur.hkl";

{* merge input data arrays with output arrays, otherwise only output
   arrays will be written *}
{+ choice: true false +}
{===>} merge_inout=false;

{===========================================================================}
{         things below this line do not normally need to be changed         }
{===========================================================================}

 ) {- end block parameter definition -}

 checkversion 1.3

 evaluate ($log_level=quiet)

 xray   
                                  
   @CNS_XTALLIB:spacegroup.lib (sg=&sg;sgparam=$sgparam;)

   a=&a b=&b c=&c  alpha=&alpha beta=&beta gamma=&gamma      

   evaluate ($counter=1)
   evaluate ($done=false)
   while ( $done = false ) loop read
    if ( &exist_reflection_infile_$counter = true ) then
      if ( &BLANK%reflection_infile_$counter = false ) then
       reflection
         @@&reflection_infile_$counter
       end
      end if
    else
     evaluate ($done=true)
    end if
    evaluate ($counter=$counter+1)
   end loop read

   declare name=&STRIP%pa_out domain=reciprocal type=real end
   declare name=&STRIP%pb_out domain=reciprocal type=real end
   declare name=&STRIP%pc_out domain=reciprocal type=real end
   declare name=&STRIP%pd_out domain=reciprocal type=real end

   group type=hl
     object=&STRIP%pa_out
     object=&STRIP%pb_out
     object=&STRIP%pc_out
     object=&STRIP%pd_out
   end

   @@CNS_XTALMODULE:check_abcd (pa=&pa_in;
                                pb=&pb_in;
                                pc=&pc_in;
                                pd=&pd_in;)

   do (&STRIP%pa_out=&blur_scale * exp(-&blur_bfactor*s*s/4.0) * &STRIP%pa_in)
      (all)
   do (&STRIP%pb_out=&blur_scale * exp(-&blur_bfactor*s*s/4.0) * &STRIP%pb_in)
      (all)
   do (&STRIP%pc_out=&blur_scale * exp(-&blur_bfactor*s*s/4.0) * &STRIP%pc_in)
      (all)
   do (&STRIP%pd_out=&blur_scale * exp(-&blur_bfactor*s*s/4.0) * &STRIP%pd_in)
      (all)

 end
 
 set display=&STRIP%reflection_outfile end

 @CNS_XTALMODULE:write_hkl_header (sg=&STRIP%sg;
                                   sgparam=$sgparam;)

 xray

   if ( &merge_inout = true ) then
      write reflection
        output=&STRIP%reflection_outfile
        sele=(all)
      end
   else
      write reflection
        &STRIP%pa_out &STRIP%pb_out &STRIP%pc_out &STRIP%pd_out
        output=&STRIP%reflection_outfile
        sele=(all)
      end
   end if
   
   set display=OUTPUT end

 end

 stop