{+ file: real_space.inp +}
{+ directory: xtal_util +}
{+ description: Computes real space correlation between a
                map and residues in the input model +}
{+ comment:
            Differences can be isomorphous, dispersive, or anomalous. 
            Phases and foms can be directly read or obtained from
            Hendrickson-Lattman coefficients.  
            Weighted average of up to 4 different Fourier maps can be computed.
            Weights are given by 1/ < df ^2 > for difference maps
            and 1 / < f^2 > for native Patterson maps. When adding
            isomorphous/dispersive and anomalous maps this weighting
            scheme is approximately equal to a weighting by f' : 2f'' +}
{+ authors: Paul D. Adams +}
{+ copyright: Yale University +}

{- 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(

{======================= molecular structure =========================}

{* molecular topology file *}
{===>} structure_infile="amy.mtf";

{* parameter files *}
{===>} parameter_infile_1="CNS_TOPPAR:protein_rep.param";
{===>} parameter_infile_2="";
{===>} parameter_infile_3="";
{===>} parameter_infile_4="";
{===>} parameter_infile_5="";

{* coordinate files *}
{* structure factors from multiple coordinates sets will be averaged *}
{===>} coordinate_infile="amy.pdb";

{====================== 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=65.508;
{===>} b=72.216;
{===>} c=45.035;
{===>} alpha=90;
{===>} beta=90;
{===>} gamma=90;

{* reflection file(s) *}
{* specify non-anomalous reflection file(s) (if any)
   before anomalous reflection file(s) *}
{===>} reflection_infile_1="amy_mir.hkl";
{===>} reflection_infile_2="";
{===>} reflection_infile_3="";

{========================= general map options =======================}

{* resolution range for scaling and map calculation *}
{+ table: rows=1 "resolution" cols=2 "lowest" "highest" +}
{===>} low_res=500.;
{===>} high_res=2.8;

{* number of maps for averaging *}
{===>} number_of_maps=1;

{* use diffraction ratios to weight Fourier maps *}
{* if false weights are set to 1 *}
{+ choice: true false +}
{===>} weight_flag=true;

{* number of bins for empirical averaging weight calculation *}
{===>} bins=10;

{============================== map 1 ================================}

{* map type *}
{+ list:  "synthesis" : map = FT ( fom * |f_a| exp(i phase) )
         "dispersive" : map = FT(fom (|f_a|-|f_b|) exp(i phase))
        "isomorphous" : map = FT(fom (|f_a|-|f_b|) exp(i phase))
          "anomalous" : map = FT(fom (|f_a+|-|f_a-|) exp(i [phase-90])) +}
{+ choice: "synthesis" "dispersive" "isomorphous" "anomalous" +}
{===>} mode_1="synthesis";

{* Hendrickson-Lattman coefficient arrays for phase and fom calculation *}
{* leave these blank if centroid phases and foms are directly read *}
{+ table: rows=1 "HL coefficients" cols=4 "A" "B" "C" "D" +}
{===>} pa_1="pa";
{===>} pb_1="pb";
{===>} pc_1="pc";
{===>} pd_1="pd";
     
{* optional: name of complex structure factor array with centroid phases *}
{* you must specify this array if no Hendrickson-Lattman coefficients are 
   specified, leave blank otherwise. *}
{===>} phase_array_1="";

{* optional: name of corresponding figure of merit array *}
{* leave this blank if no fom weighting is desired *}
{===>} fom_name_1="";

{* amplitude array f_a (for all map types) *}
{===>} f_a_1="fobs";

{* corresponding sigma array *}
{* leave blank if unavailable *}
{===>} s_a_1="sigma";

{* amplitude array f_b (for "dispersive" or "isomorphous" maps only) *}
{===>} f_b_1="";

{* corresponding sigma array (for "dispersive" or "isomorphous" maps only) *}
{* leave blank if unavailable *}
{===>} s_b_1="";

{* |F|/sigma{F} amplitude cutoff *}
{* reflections with |F|/Sigma{F} < cutoff will be rejected 
   both Freidel mates are rejected if either one has |F|/Sigma{F} < cutoff *}
{===>} cut_f_1=0.0;

{* Rms outlier cutoff *}
{+ list: synthesis maps: refl. with |F| > cutoff*rms(|F|) will be rejected
        difference maps: refl. with |delta F| > cutoff*rms(|delta F|) will 
                         be rejected +}
{===>} max_df_1=100;

{* overall k-scaling for all difference maps *}
{* if "yes": f_b will be scaled to f_a *}
{+ choice: "yes" "no" +}
{===>} kscale_1="yes";

{* B-scaling for all difference maps *}
{* if other than "no": f_b will be scaled to f_a *}
{+ choice: "no" "isotropic" "anisotropic" "anisotropic_fixed_isotropic" +}
{===>} bscale_1="anisotropic";

{* scaling iterations for all difference maps *}
{* number of iterations using rejection criteria *}
{===>} nscale_iter_1=4;

{* delta |F| cutoff for all difference maps *}
{* reflections with |delta F|/Sigma{delta F} < cutoff will be rejected *}
{===>} cut_df_1=0.;

{============================== map 2 ================================}

{* map type *}
{+ list:  "synthesis" : map = FT ( fom * |f_a| exp(i phase) )
         "dispersive" : map = FT(fom (|f_a|-|f_b|) exp(i phase))
        "isomorphous" : map = FT(fom (|f_a|-|f_b|) exp(i phase))
          "anomalous" : map = FT(fom (|f_a+|-|f_a-|) exp(i [phase-90])) +}
{+ choice: "synthesis" "dispersive" "isomorphous" "anomalous" +}
{===>} mode_2="synthesis";

{* Hendrickson-Lattman coefficient arrays for phase and fom calculation *}
{* leave these blank if centroid phases and foms are directly read *}
{+ table: rows=1 "HL coefficients" cols=4 "A" "B" "C" "D" +}
{===>} pa_2="";
{===>} pb_2="";
{===>} pc_2="";
{===>} pd_2="";
     
{* optional: name of complex structure factor array with centroid phases *}
{* you must specify this array if no Hendrickson-Lattman coefficients are
   specified, leave blank otherwise. *}
{===>} phase_array_2="";

{* optional: name of corresponding figure of merit array *}
{* leave this blank if no fom weighting is desired *}
{===>} fom_name_2="";
     
{* amplitude array f_a (for all map types) *}
{===>} f_a_2="";

{* corresponding sigma array *}
{* leave blank if unavailable *}
{===>} s_a_2="";

{* amplitude array f_b (for "dispersive" or "isomorphous" maps only) *}
{===>} f_b_2="";

{* corresponding sigma array (for "dispersive" or "isomorphous" maps only) *}
{* leave blank if unavailable *}
{===>} s_b_2="";

{* |F|/sigma{F} amplitude cutoff *}
{* reflections with |F|/Sigma{F} < cutoff will be rejected 
   both Freidel mates are rejected if either one has |F|/Sigma{F} < cutoff *}
{===>} cut_f_2=0.0;

{* Rms outlier cutoff *}
{+ list: synthesis maps: refl. with |F| > cutoff*rms(|F|) will be rejected
        difference maps: refl. with |delta F| > cutoff*rms(|delta F|) will 
                         be rejected +}
{===>} max_df_2=100;

{* overall k-scaling for all difference maps *}
{* if "yes": f_b will be scaled to f_a *}
{+ choice: "yes" "no" +}
{===>} kscale_2="yes";

{* B-scaling for all difference maps *}
{* if other than "no": f_b will be scaled to f_a *}
{+ choice: "no" "isotropic" "anisotropic" "anisotropic_fixed_isotropic" +}
{===>} bscale_2="anisotropic";

{* scaling iterations for all difference maps *}
{* number of iterations using rejection criteria *}
{===>} nscale_iter_2=4;

{* delta |F| cutoff for all difference maps *}
{* reflections with |delta F|/Sigma{delta F} < cutoff will be rejected *}
{===>} cut_df_2=0.;

{============================== map 3 ================================}

{* map type *}
{+ list:  "synthesis" : map = FT ( fom * |f_a| exp(i phase) )
         "dispersive" : map = FT(fom (|f_a|-|f_b|) exp(i phase))
        "isomorphous" : map = FT(fom (|f_a|-|f_b|) exp(i phase))
          "anomalous" : map = FT(fom (|f_a+|-|f_a-|) exp(i [phase-90])) +}
{+ choice: "synthesis" "dispersive" "isomorphous" "anomalous" +}
{===>} mode_3="synthesis";

{* Hendrickson-Lattman coefficient arrays for phase and fom calculation *}
{* leave these blank if centroid phases and foms are directly read *}
{+ table: rows=1 "HL coefficients" cols=4 "A" "B" "C" "D" +}
{===>} pa_3="";
{===>} pb_3="";
{===>} pc_3="";
{===>} pd_3="";
     
{* optional: name of complex structure factor array with centroid phases *}
{* you must specify this array if no Hendrickson-Lattman coefficients are
   specified, leave blank otherwise. *}
{===>} phase_array_3="";

{* optional: name of corresponding figure of merit array *}
{* leave this blank if no fom weighting is desired *}
{===>} fom_name_3="";
     
{* amplitude array f_a (for all map types) *}
{===>} f_a_3="";

{* corresponding sigma array *}
{* leave blank if unavailable *}
{===>} s_a_3="";

{* amplitude array f_b (for "dispersive" or "isomorphous" maps only) *}
{===>} f_b_3="";

{* corresponding sigma array (for "dispersive" or "isomorphous" maps only) *}
{* leave blank if unavailable *}
{===>} s_b_3="";

{* |F|/sigma{F} amplitude cutoff *}
{* reflections with |F|/Sigma{F} < cutoff will be rejected 
   both Freidel mates are rejected if either one has |F|/Sigma{F} < cutoff *}
{===>} cut_f_3=0.0;

{* Rms outlier cutoff *}
{+ list: synthesis maps: refl. with |F| > cutoff*rms(|F|) will be rejected
        difference maps: refl. with |delta F| > cutoff*rms(|delta F|) will 
                         be rejected +}
{===>} max_df_3=100;

{* overall k-scaling for all difference maps *}
{* if "yes": f_b will be scaled to f_a *}
{+ choice: "yes" "no" +}
{===>} kscale_3="yes";

{* B-scaling for all difference maps *}
{* if other than "no": f_b will be scaled to f_a *}
{+ choice: "no" "isotropic" "anisotropic" "anisotropic_fixed_isotropic" +}
{===>} bscale_3="anisotropic";

{* scaling iterations for all difference maps *}
{* number of iterations using rejection criteria *}
{===>} nscale_iter_3=4;

{* delta |F| cutoff for all difference maps *}
{* reflections with |delta F|/Sigma{delta F} < cutoff will be rejected *}
{===>} cut_df_3=0.;

{============================== map 4 ================================}

{* map type *}
{+ list:  "synthesis" : map = FT ( fom * |f_a| exp(i phase) )
         "dispersive" : map = FT(fom (|f_a|-|f_b|) exp(i phase))
        "isomorphous" : map = FT(fom (|f_a|-|f_b|) exp(i phase))
          "anomalous" : map = FT(fom (|f_a+|-|f_a-|) exp(i [phase-90])) +}
{+ choice: "synthesis" "dispersive" "isomorphous" "anomalous" +}
{===>} mode_4="synthesis";

{* Hendrickson-Lattman coefficient arrays for phase and fom calculation *}
{* leave these blank if centroid phases and foms are directly read *}
{+ table: rows=1 "HL coefficients" cols=4 "A" "B" "C" "D" +}
{===>} pa_4="";
{===>} pb_4="";
{===>} pc_4="";
{===>} pd_4="";
     
{* optional: name of complex structure factor array with centroid phases *}
{* you must specify this array if no Hendrickson-Lattman coefficients are 
   specified, leave blank otherwise. *}
{===>} phase_array_4="";

{* optional: name of corresponding figure of merit array *}
{* leave this blank if no fom weighting is desired *}
{===>} fom_name_4="";
     
{* amplitude array f_a (for all map types) *}
{===>} f_a_4="";

{* corresponding sigma array *}
{* leave blank if unavailable *}
{===>} s_a_4="";

{* amplitude array f_b (for "dispersive" or "isomorphous" maps only) *}
{===>} f_b_4="";

{* corresponding sigma array (for "dispersive" or "isomorphous" maps only) *}
{* leave blank if unavailable *}
{===>} s_b_4="";

{* |F|/sigma{F} amplitude cutoff *}
{* reflections with |F|/Sigma{F} < cutoff will be rejected 
   both Freidel mates are rejected if either one has |F|/Sigma{F} < cutoff *}
{===>} cut_f_4=0.0;

{* Rms outlier cutoff *}
{+ list: synthesis maps: refl. with |F| > cutoff*rms(|F|) will be rejected
        difference maps: refl. with |delta F| > cutoff*rms(|delta F|) will 
                         be rejected +}
{===>} max_df_4=100;

{* overall k-scaling for all difference maps *}
{* if "yes": f_b will be scaled to f_a *}
{+ choice: "yes" "no" +}
{===>} kscale_4="yes";

{* B-scaling for all difference maps *}
{* if other than "no": f_b will be scaled to f_a *}
{+ choice: "no" "isotropic" "anisotropic" "anisotropic_fixed_isotropic" +}
{===>} bscale_4="anisotropic";

{* scaling iterations for all difference maps *}
{* number of iterations using rejection criteria *}
{===>} nscale_iter_4=4;

{* delta |F| cutoff for all difference maps *}
{* reflections with |delta F|/Sigma{delta F} < cutoff will be rejected *}
{===>} cut_df_4=0.;

{========================== atom selection ===========================}

{* select atoms to be included in real space correlation *}
{===>} atom_select=(known and not hydrogen);

{============================ map output options ======================}

{* write map file *}
{+ choice: true false +}
{===>} write_map=true;

{* scale map by dividing by the rms sigma of the map *}
{* otherwise map will be on the scale of the first data set *}
{+ choice: true false +}
{===>} map_scale=true;

{* map format *}
{+ choice: "cns" "ezd" +}
{===>} map_format="cns";

{* map grid size: dmin*grid *}
{* use grid=0.25 for better map appearance *}
{===>} grid=0.33;

{* memory allocation for FFT calculation *}
{* this will be determined automatically if a negative value is given
   otherwise the specified number of words of space will be allocated *}
{===>} fft_memory=-1;

{* extent of map *}
{+ choice: "asymmetric" "box" "unit" "molecule" "fract" +}
{===>} map_mode="asymmetric";

{* limits in orthogonal angstroms for box mode or
   fractional coordinates for fract mode *}
{+ table: rows=3 "x" "y" "z" cols=2 "minimum" "maximum" +}
{===>} xmin=0.;
{===>} xmax=0.;
{===>} ymin=0.;
{===>} ymax=0.;
{===>} zmin=0.;
{===>} zmax=0.;

{* atom selection for writing map around molecule *}
{* the map will cover the selected atoms *}
{===>} atom_map=(known and not hydrogen);

{* cushion (in Angstroms) around selected atoms *}
{===>} map_cushion=3.0;

{============================ peak picking ===========================}

{* do peak picking on map *}
{+ choice: true false +}
{===>} peak_search=true;

{* number of highest peaks *}
{===>} peak_number=30;

{* optional O-mask for restricting peak picking *}
{===>} mask_infile="";

{* set electron density in map to zero outside of mask *}
{+ choice: true false +}
{===>} zero_map=false;

{============================ output filenames =======================}

{* root for output files *}
{+ list: 
        map file will be in: <output_root>.map
        real space correlations in: <output_root>_r.list
        positive peaks in: <output_root>_positive.peaks
        negative peaks in: <output_root>_negative.peaks +}
{===>} output_root="real_space";

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

 ) {- end block parameter definition -}

checkversion 1.2

evaluate ($log_level=quiet)

set echo=off end

 structure @&structure_infile end

 coordinates @&coordinate_infile

 parameter
   if ( &BLANK%parameter_infile_1 = false ) then
      @@&parameter_infile_1
   end if
   if ( &BLANK%parameter_infile_2 = false ) then
      @@&parameter_infile_2
   end if
   if ( &BLANK%parameter_infile_3 = false ) then
      @@&parameter_infile_3
   end if
   if ( &BLANK%parameter_infile_4 = false ) then
      @@&parameter_infile_4
   end if
   if ( &BLANK%parameter_infile_5 = false ) then
      @@&parameter_infile_5
   end if
 end

 xray

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

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

   @CNS_XRAYLIB:scatter.lib

   if ( &BLANK%reflection_infile_1 = false ) then
      reflection @@&reflection_infile_1 end  
   end if                
   if ( &BLANK%reflection_infile_2 = false ) then
      reflection @@&reflection_infile_2 end
   end if               
   if ( &BLANK%reflection_infile_3 = false ) then
      reflection @@&reflection_infile_3 end
   end if               

   binresolution &low_res &high_res
   mapresolution &high_res

   fft
     grid=&grid
     if ( &fft_memory < 0 ) then
       automemory=true
     else
       memory=&fft_memory
     end if
   end

   statistics overall
     completeness
     selection=( &low_res >= d >= &high_res )
   end
   evaluate ($total_compl=$expression1)

   show sum(1) ( &low_res >= d >= &high_res )
   evaluate ($total_read=$select)
   evaluate ($total=int(1./$total_compl * $total_read))

end

xray
   declare name=df_sel domain=reciprocal  type=complex end
   declare name=df     domain=reciprocal  type=real end
   declare name=df_accum     domain=reciprocal  type=complex end
   declare name=s_df   domain=reciprocal  type=real end
   declare name=weight domain=reciprocal type=real end

   do (df_accum=0) ( all )

   if (&number_of_maps > 1) then
      buffer harvest
         display Average over &number_of_maps Fourier maps computed.  Weight=&weight_flag
      end
   end if
   
   buffer harvest
     display resolution: &low_res - &high_res  
     display sg= &STRIP%sg 
     display a= &a b= &b c= &c alpha= &alpha beta= &beta gamma= &gamma
     display map is scaled to sigma units : &map_scale
   end
      
   if ( $log_level = verbose ) then
     set echo=on end
   else
     set echo=off end
   end if

   eval ($nm=1) 
   while ($nm <= &number_of_maps) loop main
   
      if ( &BLANK%pa_$nm = false ) then
         {- Hendrickson-Lattman coefficients are specified -> compute
            phases and figure of merit arrays -}
            
         @@CNS_XTALMODULE:check_abcd (pa=&pa_$nm;
                                     pb=&pb_$nm;
                                     pc=&pc_$nm;
                                     pd=&pd_$nm;)
         if ( &BLANK%fom_name_$nm = false ) then
            display
            display *********************************************************
            display map=$nm : error, both Hendrickson-Lattman and fom
            display           specified.  Use either one or the other. 
            display *********************************************************
            display 
            abort
         elseif ( &BLANK%phase_array_$nm = false ) then
            display
            display *********************************************************
            display map=$nm : error, both Hendrickson-Lattman and centroid
            display           phases specified.  Use either one or the other. 
            display *********************************************************
            display 
            abort
         end if

         evaluate ($phase_array_$nm="phase_tem"+encode($nm))
         declare name=$STRIP%phase_array_$nm type=complex domain=reciprocal end
         
         @CNS_XTALMODULE:getfom (pa=&STRIP%pa_$nm;
                                 pb=&STRIP%pb_$nm;
                                 pc=&STRIP%pc_$nm;
                                 pd=&STRIP%pd_$nm;
                                 m=$STRIP%phase_array_$nm;
                                 phistep=5; )
                                 
         evaluate ($fom_name_$nm="fom_tem"+encode($nm))  
         declare name=$STRIP%fom_name_$nm type=real domain=reciprocal end
         do ($STRIP%fom_name_$nm=amplitude($STRIP%phase_array_$nm)) ( all )
                                 
        
      else

         evaluate ($fom_name_$nm=&fom_name_$nm) 
         if ( &BLANK%fom_name_$nm = true ) then
            evaluate ($fom_name_$nm="fom_tem"+encode($nm))
            declare name=$STRIP%fom_name_$nm type=real domain=reciprocal end
            do ($STRIP%fom_name_$nm=1) ( all )
            display
            display *********************************************************
            display map=$nm : figure of merit array not specified -> set to 1
            display *********************************************************
            display 
         end if
 
         query name=$STRIP%fom_name_$nm domain=reciprocal end
         if ( $object_exist=false ) then
            display
            display ********************************************************************
            display map=$nm : error, figure of merit array $STRIP%fom_name_$nm does not exist
            display ********************************************************************
            display 
            abort
         end if
         if ( $object_type # "REAL" ) then
            display 
            display *****************************************************************************
            display map=$nm : error, figure of merit array $STRIP%fom_name_$nm has the wrong data type
            display *****************************************************************************
            display
            abort
         end if
      
         if ( &BLANK%phase_array_$nm = true ) then
           display 
           display  ***********************************************************
           display  map=$nm : error, complex array with phases is not specified
           display  ***********************************************************
           display
           abort
         else
            evaluate ($phase_array_$nm=&phase_array_$nm) 
            query name=&STRIP%phase_array_$nm domain=reciprocal end
            if ( $object_exist=false ) then
               display
               display *************************************************************************
               display map=$nm : error, complex array with phases &phase_array_$nm does not exist
               display *************************************************************************
               display 
               abort
            end if
            if ( $object_type # "COMPLEX" ) then
               display 
               display **********************************************************************************
               display map=$nm : error, complex array with phases &phase_array_$nm has the wrong data type
               display **********************************************************************************
               display
               abort
            end if
         end if
      end if
      
      if ( &BLANK%f_a_$nm = true ) then
        display 
        display  *****************************************************
        display  map=$nm : error, f_a amplitude array is not specified
        display  *****************************************************
        display
        abort
      else
         query name=&STRIP%f_a_$nm domain=reciprocal end
         if ( $object_exist=false ) then
            display 
            display *************************************************************
            display map=$nm : error, f_a amplitude array &f_a_$nm does not exist
            display *************************************************************
            display 
            abort
         end if
         {- note: this array can be of any type -}
      end if
      
      evaluate ($s_a_$nm=&s_a_$nm)    
      if ( &BLANK%s_a_$nm = true ) then
         evaluate ($s_a_$nm="s_tema"+encode($nm)) 
         declare name=$STRIP%s_a_$nm type=real domain=reciprocal end
         do ($STRIP%s_a_$nm=0.001) ( all ) 
         display 
         display *************************************************
         display map=$nm : f_a sigma array missing -> set to 0.001
         display *************************************************
         display 
      else
         query name=$STRIP%s_a_$nm domain=reciprocal end
         if ( $object_exist=false ) then
            display
            display *********************************************************
            display map=$nm : error, f_a sigma array &s_a_$nm does not exist
            display *********************************************************
            display
            abort
         end if
         if ( $object_type # "REAL" ) then
            display 
            display **********************************************************************
            display map=$nm : error, f_a set sigma array &s_a_$nm has the wrong data type
            display **********************************************************************
            display
            abort
         end if
      end if
     
      if (&mode_$nm="synthesis") then
      
         do (df=amplitude(&STRIP%f_a_$nm)) ( all ) 
                                                     
         show rms (df) ( df > 0 and &low_res >= d >= &high_res ) 
         evaluate ($rms_df=$result)

         do (df_sel=$STRIP%fom_name_$nm * combine( amplitude(&STRIP%f_a_$nm),phase($STRIP%phase_array_$nm)) )
            ( amplitude($STRIP%phase_array_$nm)>0 and amplitude(&STRIP%f_a_$nm)> &cut_f_$nm * $STRIP%s_a_$nm
              and df < &max_df_$nm * $rms_df and  &low_res >= d >= &high_res)

 
         eval ($total_used=$select) 
         eval ($per_used=100*$total_used/$total) 
         show sum ( 1 ) ( &low_res >= d >= &high_res and amplitude(&STRIP%f_a_$nm)>0  ) 
         eval ($total_rejected=$select-$total_used) 
         eval ($per_rejected=100*($total_rejected)/$total) 
         eval ($total_unobserved=$total-$select)
         eval ($per_unobserved=100*($total_unobserved)/$total)


         buffer harvest
            display map $nm : &STRIP%mode_$nm map: 
         end
         evaluate ($remark="fom *" + &f_a_$nm + " * exp(i phase)")  
          
         buffer harvest     
            display $remark
         end
         
         if ( &BLANK%fom_name_$nm = false ) then
            buffer harvest
               display fom= &STRIP%fom_name_$nm  phase= &STRIP%phase_array_$nm
            end
         elseif ( &BLANK%phase_array_$nm = false ) then
            buffer harvest
               display fom=1 phase= &STRIP%phase_array_$nm
            end
         else
            buffer harvest
               display fom and phase obtained from HL coefficients &STRIP%pa_$nm, &STRIP%pb_$nm, &STRIP%pc_$nm, &STRIP%pd_$nm
            end
         end if
        
         buffer harvest 
            display reflections with &STRIP%f_a_$nm / $STRIP%s_a_$nm < &cut_f_$nm rejected
            display reflections with &STRIP%f_a_$nm > &max_df_$nm * rms( &STRIP%f_a_$nm ) rejected
         end

     elseif (&mode_$nm="anomalous") then

         declare name=f_a_fried domain=reciprocal type=real end
         declare name=s_a_fried domain=reciprocal type=real end
         
         do (f_a_fried=friedel(amplitude(&STRIP%f_a_$nm))) 
                                     ( friedel_pair(amplitude(&STRIP%f_a_$nm)>0))
         do (s_a_fried=friedel(amplitude(&STRIP%s_a_$nm))) 
                                     ( friedel_pair(amplitude(&STRIP%f_a_$nm)>0))
                                     

         buffer harvest
            display map $nm : &STRIP%mode_$nm difference map: 
         end
         evaluate ($remark=" fom * " + "[ " + &f_a_$nm + "+  -  " + &f_a_$nm )
         evaluate ($remark=$remark + "- ] * exp(i [ phase -90 ] )" )
         
         buffer harvest
            display $remark
         end
         if ( &BLANK%fom_name_$nm = false ) then
            buffer harvest 
               display fom= &STRIP%fom_name_$nm  phase= &STRIP%phase_array_$nm
            end
         elseif ( &BLANK%phase_array_$nm = false ) then
            buffer harvest 
               display fom=1 phase= &STRIP%phase_array_$nm
            end
         else
            buffer harvest
               display fom and phase obtained from HL coefficients &STRIP%pa_$nm, &STRIP%pb_$nm, &STRIP%pc_$nm, &STRIP%pd_$nm
            end
         end if
         
         buffer harvest
            display scaling applied to Friedel mates of &STRIP%f_a_$nm : k-scaling = &kscale_$nm   b-scaling = &bscale_$nm
         end
         
          {- intial scale factor: Friedel to wa -}
         @CNS_XTALMODULE:scalef (		    
            kscale=&kscale_$nm;	       
            bscale=&bscale_$nm;	
            sel=( friedel_pair(&low_res >= d >= &high_res 
               and amplitude(&STRIP%f_a_$nm)> &cut_f_$nm * $STRIP%s_a_$nm));  
            apply="all"; 
            fref=&STRIP%f_a_$nm;
            f=f_a_fried;
            s=s_a_fried;
            buffer_name=harvest;
            )

         buffer harvest
            concatenate  (iteration no. 0)
         end
                                     
         do (df=amplitude(&STRIP%f_a_$nm)-f_a_fried)
                                     ( friedel_pair(amplitude(&STRIP%f_a_$nm)>0)) 
                                                    
                                                     
         do (s_df=sqrt($STRIP%s_a_$nm^2 + s_a_fried^2) ) 
                                     ( friedel_pair( amplitude(&STRIP%f_a_$nm)>0))

         show rms (df) ( abs(df)>0 and  &low_res >= d >= &high_res ) 
         evaluate ($rms_df=$result)

         evaluate ($iter=0)
         
         {- iterate scale factor using rejection criteria -}
         while ($iter < &nscale_iter_$nm) loop sitr

            evaluate ($iter=$iter+1)
            
            @CNS_XTALMODULE:scalef (		    
               kscale=&kscale_$nm;	       
               bscale=&bscale_$nm;	
               sel=(friedel_pair( amplitude(&STRIP%f_a_$nm) > &cut_f_$nm * $STRIP%s_a_$nm 
                         and abs(df) > &cut_df_$nm * s_df 
                         and abs(df) < &max_df_$nm * $rms_df 
                         and &low_res >= d >= &high_res ) ) ;
               apply="all"; 
               fref=&STRIP%f_a_$nm;
               f=f_a_fried;
               s=s_a_fried;
               buffer_name=harvest;
               )

            buffer harvest
               concatenate  (iteration no. $iter)
            end
               
            do (df=amplitude(&STRIP%f_a_$nm)-f_a_fried)
                                     ( friedel_pair(amplitude(&STRIP%f_a_$nm)>0)) 
                                                    
                                                     
            do (s_df=sqrt($STRIP%s_a_$nm^2 + s_a_fried^2) ) 
                                     ( friedel_pair( amplitude(&STRIP%f_a_$nm)>0))
  
            show rms (df) ( abs(df)>0 and &low_res >= d >= &high_res) 
            evaluate ($rms_df=$result)

         end loop sitr
      
         do (df_sel= $STRIP%fom_name_$nm * combine( df , phase($STRIP%phase_array_$nm)- 90.))
            (friedel_pair( 
               amplitude($STRIP%phase_array_$nm)>0 and amplitude(&STRIP%f_a_$nm)> &cut_f_$nm * $STRIP%s_a_$nm   
               and  abs(df) > &cut_df_$nm * s_df and abs(df) < &max_df_$nm * $rms_df) 
                                        and &low_res >= d >= &high_res)


         eval ($total_used=$select) 
         eval ($per_used=100*$total_used/$total) 
         show sum ( 1 ) ( &low_res >= d >= &high_res and amplitude(&STRIP%f_a_$nm)>0  ) 
         eval ($total_rejected=$select-$total_used) 
         eval ($per_rejected=100*($total_rejected)/$total) 
         eval ($total_unobserved=$total-$select)
         eval ($per_unobserved=100*($total_unobserved)/$total) 

         buffer harvest
            display reflections with  &STRIP%f_a_$nm / $STRIP%s_a_$nm < &cut_f_$nm rejected
            display reflections with | delta &STRIP%f_a_$nm | / sigma(| delta &STRIP%f_a_$nm |) \
< &cut_df_$nm rejected
            display reflections with | delta &STRIP%f_a_$nm | > &max_df_$nm * rms(| delta \
&STRIP%f_a_$nm |) rejected
            display all centric reflections rejected
         end

         undeclare name=f_a_fried domain=reciprocal end
         undeclare name=s_a_fried domain=reciprocal end
   
       else  { dispersive or isomorphous }
       
         if ( &mode_$nm # "isomorphous" ) then
            if ( &mode_$nm # "dispersive" ) then
               display 
               display  ********************************************
               display  map=$nm : error, unknown map mode: &mode_$nm
               display  ********************************************
               display
               abort
            end if
         end if

         if ( &BLANK%f_b_$nm = true ) then
           display 
           display  *****************************************************
           display  map=$nm : error, f_b amplitude array is not specified
           display  *****************************************************
           display
           abort
         else
            query name=&STRIP%f_b_$nm domain=reciprocal end
            if ( $object_exist=false ) then
               display 
               display ************************************************************
               display map=$nm : error, f_b amplitude array &f_b_$nm does not exist
               display ************************************************************
               display
               abort
            end if
           {- note: this array can be of any type -}
         end if
         
         evaluate ($s_b_$nm=&s_b_$nm) 
         
         if ( &BLANK%s_b_$nm = true ) then
            evaluate ($s_b_$nm="s_temb"+encode($nm)) 
            declare name=$STRIP%s_b_$nm type=real domain=reciprocal end
            do ($STRIP%s_b_$nm=0.001) ( all ) 
            display 
            display *************************************************
            display map=$nm : f_b sigma array missing -> set to 0.001
            display *************************************************
            display 
         else
            query name=&STRIP%s_b_$nm domain=reciprocal end
            if ( $object_exist=false ) then
               display 
               display ********************************************************
               display map=$nm : error, f_b sigma array &s_b_$nm does not exist
               display ********************************************************
               display 
               abort
            end if
            if ( $object_type # "REAL" ) then
               display  							
               display *****************************************************************
               display map=$nm : error, f_b sigma array &s_b_$nm has the wrong data type
               display *****************************************************************
               display
               abort
            end if
         end if

         buffer harvest
            display map $nm : &STRIP%mode_$nm difference map:
         end
         evaluate ($remark=" fom * " + "[ " + &f_a_$nm + " - " + &f_b_$nm )
         evaluate ($remark=$remark + " ] * exp(i phase )" )
         buffer harvest
            display $remark
         end
         if ( &BLANK%fom_name_$nm = false ) then
            buffer harvest
               display fom= &STRIP%fom_name_$nm  phase= &STRIP%phase_array_$nm
            end
         elseif ( &BLANK%phase_array_$nm = false ) then
            buffer harvest
               display fom=1 phase= &STRIP%phase_array_$nm
            end
         else
            buffer harvest
               display fom and phase obtained from HL coefficients &STRIP%pa_$nm, &STRIP%pb_$nm, &STRIP%pc_$nm, &STRIP%pd_$nm
            end
         end if

         buffer harvest
            display scaling applied to &STRIP%f_b_$nm : k-scaling = &kscale_$nm   b-scaling = &bscale_$nm
         end
         
         {- Scale wb to wa -}
         @CNS_XTALMODULE:scalef (		    
            kscale=&kscale_$nm;	       
            bscale=&bscale_$nm;	
            sel=( &low_res >= d >= &high_res 
               and amplitude(&STRIP%f_a_$nm)> &cut_f_$nm * $STRIP%s_a_$nm 
               and amplitude(&STRIP%f_b_$nm)> &cut_f_$nm * $STRIP%s_b_$nm);  
            apply="all";
            fref=&STRIP%f_a_$nm;
            f=&STRIP%f_b_$nm;
            s=$STRIP%s_b_$nm;
            buffer_name=harvest;
            )

         buffer harvest
            concatenate  (iteration no. 0)
         end

         do (df=amplitude(&STRIP%f_a_$nm)-amplitude(&STRIP%f_b_$nm)) 
                  ( amplitude(&STRIP%f_a_$nm)>0 and amplitude(&STRIP%f_b_$nm)>0) 

         do (s_df=sqrt($STRIP%s_a_$nm^2 + $STRIP%s_b_$nm^2) ) 
                  ( amplitude(&STRIP%f_a_$nm)>0 and amplitude(&STRIP%f_b_$nm)>0)
  
         show rms (df) ( abs(df)>0 and &low_res >= d >= &high_res) 
         evaluate ($rms_df=$result)

         evaluate ($iter=0)
         
         {- iterate scale factor using rejection criteria -}
         while ($iter < &nscale_iter_$nm) loop sitr

            evaluate ($iter=$iter+1)
            
            @CNS_XTALMODULE:scalef (		    
               kscale=&kscale_$nm;	       
               bscale=&bscale_$nm;	
               sel=(  amplitude(&STRIP%f_a_$nm) > &cut_f_$nm * $STRIP%s_a_$nm 
                  and amplitude(&STRIP%f_b_$nm) > &cut_f_$nm * $STRIP%s_b_$nm 
                  and abs(df) > &cut_df_$nm * s_df 
                  and abs(df) < &max_df_$nm * $rms_df 
                  and &low_res >= d >= &high_res ) ;
               apply="all"; 
               fref=&STRIP%f_a_$nm;
               f=&STRIP%f_b_$nm;
               s=$STRIP%s_b_$nm;
               buffer_name=harvest;
               )

            buffer harvest
               concatenate  (iteration no. $iter)
            end
               
            do (df=amplitude(&STRIP%f_a_$nm) - amplitude(&STRIP%f_b_$nm) )  
                      ( amplitude(&STRIP%f_a_$nm)>0 and amplitude(&STRIP%f_b_$nm)>0) 

            do (s_df=sqrt($STRIP%s_a_$nm^2 + $STRIP%s_b_$nm^2) )  
                      ( amplitude(&STRIP%f_a_$nm)>0 and amplitude(&STRIP%f_b_$nm)>0) 

            show rms (df) ( abs(df)>0 and &low_res >= d >= &high_res ) 
            evaluate ($rms_df=$result)

         end loop sitr
         
         do (df_sel=$STRIP%fom_name_$nm * combine( df , phase($STRIP%phase_array_$nm) )  )
            ( amplitude($STRIP%phase_array_$nm)>0 and amplitude(&STRIP%f_a_$nm)> &cut_f_$nm * $STRIP%s_a_$nm
               and amplitude(&STRIP%f_b_$nm)> &cut_f_$nm * $STRIP%s_b_$nm
               and  abs(df) > &cut_df_$nm * s_df and abs(df) < &max_df_$nm * $rms_df             
               and  &low_res >= d >= &high_res)

         eval ($total_used=$select) 
         eval ($per_used=100*$total_used/$total) 
         show sum ( 1 ) ( &low_res >= d >= &high_res 
                          and ( amplitude(&STRIP%f_a_$nm)>0 and amplitude(&STRIP%f_b_$nm)>0 ) )
         eval ($total_rejected=$select-$total_used) 
         eval ($per_rejected=100*($total_rejected)/$total) 
         eval ($total_unobserved=$total-$select)
         eval ($per_unobserved=100*($total_unobserved)/$total)
         
         buffer harvest 
            display reflections with  &STRIP%f_a_$nm / $STRIP%s_a_$nm < &cut_f_$nm \
and/or  &STRIP%f_b_$nm / $STRIP%s_b_$nm < &cut_f_$nm rejected
            display reflections with | &STRIP%f_a_$nm-&STRIP%f_b_$nm | / \
sigma(| &STRIP%f_a_$nm-&STRIP%f_b_$nm |) < &cut_df_$nm rejected
            display reflections with | &STRIP%f_a_$nm-&STRIP%f_b_$nm | > \
&max_df_$nm * rms(| &STRIP%f_a_$nm-&STRIP%f_b_$nm |) rejected
         end 
     
      end if

      evaluate ($weight_flag=&weight_flag) 
      if (&number_of_maps=1) then
          evaluate ($weight_flag=false)  
      end if

      if ($weight_flag=true) then
         if ( &bins > 0 ) then
           bins=&bins
         end if
         do (weight=max(0.001,save(amplitude(df_sel)^2)))  ( amplitude(df_sel)>0 ) 
         do (df_accum=df_accum+ df_sel/ sqrt(weight) )     ( amplitude(df_sel)>0 ) 
      else
         do (df_accum=df_sel)                              ( amplitude(df_sel)>0 ) 
      end if

      buffer harvest
         display theoretical total number of refl. in resol. range: \
$total[I7] ( 100.0  % )
         display number of unobserved reflections (no entry or f=0):\
$total_unobserved[I7] ( $per_unobserved[F5.1]  % )
         display number of reflections rejected:                    \
$total_rejected[I7] ( $per_rejected[F5.1]  % )
         display total number of reflections used:                  \
$total_used[I7] ( $per_used[F5.1]  % )
      end

      evaluate ($nm=$nm+1)
   end loop main  

   set echo=on end

   buffer harvest
      to=remarks
      dump
   end
    
   declare name=map1 domain=real end
      
   do (map1=ft( df_accum )) ( amplitude(df_accum)>0) 
   
   {- put map on a sigma (rms) scale -}   
   if (&map_scale=true) then                                            
      show rms (real(map1)) ( all )
      do (map1=map1/$result) ( all )
   end if
   
   if ( &BLANK%mask_infile = false ) then
     declare name=mask domain=real end
     read mask
       to=mask
       input=&mask_infile
       type=omask
     end
   end if

   if ( &zero_map = true ) then
     if ( &BLANK%mask_infile = false ) then
       do (map1=0) ( mask > 0 )
     end if
   end if

   if ( &map_mode="box" ) then
      eval ($map_mode_string=BOX) 
   elseif ( &map_mode="molecule" ) then
      eval ($map_mode_string=MOLE)
   elseif ( &map_mode="unit" ) then
      eval ($map_mode_string=UNIT)
   elseif ( &map_mode = "fract" ) then
     evaluate ($map_mode_string=FRAC)
   else
      eval ($map_mode_string=ASYM)
   end if

   if ( &write_map = true ) then
     evaluate ($filename=&output_root + ".map") 
     write map 
       if ( &map_format = "ezd" ) then
         type=ezd
       else
         type=cns
       end if
       automatic=false
       from=map1
       extend=$map_mode_string
       if ( &map_mode = "molecule" ) then
         selection=&atom_map
         cushion=&map_cushion
       end if
       output=$filename
       if ( &map_mode = "box" ) then
         xmin=&xmin xmax=&xmax
         ymin=&ymin ymax=&ymax
         zmin=&zmin zmax=&zmax
       end if
       if ( &map_mode = "fract" ) then
         xmin=&xmin xmax=&xmax
         ymin=&ymin ymax=&ymax
         zmin=&zmin zmax=&zmax
       end if
     end  
   end if                                       
 end

 xray

   query name=fcalc domain=reciprocal end
   if ( $object_exist = false ) then
      declare name=fcalc domain=reciprocal type=complex end
   end if

   associate fcalc ( &atom_select )

   predict
     mode=reciprocal
     to=fcalc
     selection=(all)
     atomselection=( &atom_select ) 
   end

   declare name=mod_map domain=real end

   do (mod_map=ft(fcalc)) (&low_res >= d >= &high_res)

   if (&map_scale=true) then    
     show rms (real(mod_map)) ( all )
     do (mod_map=mod_map/$result) ( all )
   end if

   declare name=corr_map domain=real end
   declare name=map11    domain=real end
   declare name=map12    domain=real end
   declare name=map22    domain=real end
   declare name=prop     domain=real end
   declare name=dist     domain=real end
   declare name=mask     domain=real end
 end

 do (store9=0) (all)
 evaluate ($counter=1)
 for $id in id ( tag and &atom_select ) loop main
   do ( store9=$counter ) ( byres( id $id) )
   evaluate ($counter=$counter+1)
 end loop main

 xray
   mask
     mode=vdw
     solrad=0.1
     shrink=0.1
     nshell=1
     to=mask
     selection=( &atom_select ) 
   end

   do (prop=0) (all)
   proximity
     from=store9
     distance_map=dist
     property_map=prop
     cutoff=3.0
     selection=( &atom_select )
   end

   do (map11=gave(real(map1) * real(map1), real(prop)) ) 
                            ( real(prop) > 0 and real(mask) <= 0 )
   do (map12=gave(real(map1) * real(mod_map), real(prop)) ) 
                            ( real(prop) > 0 and real(mask) <= 0 )
   do (map22=gave(real(mod_map) * real(mod_map), real(prop)) ) 
                            ( real(prop) > 0 and real(mask) <= 0 )
   do (corr_map=real(map12)/ sqrt ( real(map11) * real(map22) ) )
                            ( real(prop) > 0 and real(mask) <= 0 )
 end

 evaluate ($display=&output_root + "_r.list")
 set display=$display end

 display    #   resid   segid   CC   R-value

 evaluate ($counter=1)
 for $id in id ( tag and &atom_select ) loop real
   xray
     show ave (real(corr_map)) (real(prop)=$counter and real(mask) <= 0)
   end
   evaluate ($corr=$result)
   evaluate ($realr=1-$corr)
   show (resid) (id $id)
   evaluate ($resid=$result)
   show (segid) (id $id)
   evaluate ($segid=$result)
   display $counter[i6]   $resid[a4]   $segid[a4] $corr[f6.3] $realr[f6.3]
   evaluate ($counter=$counter+1)
 end loop real

 xray
   undeclare name=mod_map  domain=real end
   undeclare name=corr_map domain=real end
   undeclare name=map11    domain=real end
   undeclare name=map12    domain=real end
   undeclare name=map22    domain=real end
   undeclare name=prop     domain=real end
   undeclare name=dist     domain=real end
   undeclare name=mask     domain=real end
 end

 if ( &peak_search = true ) then

    evaluate ($filename=&output_root + "_positive.peaks") 
    xray  
      peakpik
         from=map1
         mpeak=&peak_number
         if ( &BLANK%mask_infile = false ) then
           selection=( mask <= 0 ) 
         else
           selection=( all )
         end if
         atom=true
         if ( &map_mode = "molecule" ) then
           proximity=(&atom_map)
         end if
      end
    end 
    write coor output=$filename selection=(segid PEAK) end
    
    delete selection=(segid PEAK) end

    evaluate ($filename=&output_root + "_negative.peaks") 
    xray  
      do (map1=-map1) ( all ) 
      peakpik
         from=map1
         mpeak=&peak_number
         if ( &BLANK%mask_infile = false ) then
           selection=( mask <= 0 ) 
         else
           selection=( all )
         end if
         atom=true
         if ( &map_mode = "molecule" ) then
           proximity=(&atom_map)
         end if
      end
    end
    write coor output=$filename selection=(segid PEAK) end

 end if
                                              
stop