{+ file: model_phase.inp +}
{+ directory: xtal_util +}
{+ description: Generate calculated amplitudes, phases, Hendrickson-Lattman 
                coefficients, FOMs from model(s) +}
{+ comments:
            The output amplitudes and phases include the contribution from bulk solvent. +}
{+ authors: Axel T. Brunger and Paul D. Adams +}
{+ copyright: Yale University +}

{- Guidelines for using this file:
   - all strings must be quoted by double-quotes
   - logical variables (true/false) must not be quoted
   - do not remove any evaluate statements from the file -}

{- begin block parameter definition -} define(

{============================ coordinates ============================}

{* coordinate file *}
{===>} coordinate_infile="amy.pdb";

{==================== molecular information ==========================}

{* topology files *}
{===>} topology_infile_1="CNS_TOPPAR:protein.top";
{===>} topology_infile_2="CNS_TOPPAR:dna-rna.top";
{===>} topology_infile_3="CNS_TOPPAR:water.top";
{===>} topology_infile_4="CNS_TOPPAR:ion.top";
{===>} topology_infile_5="CNS_TOPPAR:carbohydrate.top";
{===>} topology_infile_6="";
{===>} topology_infile_7="";
{===>} topology_infile_8="";

{* linkage files for linear, continuous polymers (protein, DNA, RNA) *}
{===>} link_infile_1="CNS_TOPPAR:protein.link";
{===>} link_infile_2="CNS_TOPPAR:dna-rna-pho.link";
{===>} link_infile_3="";

{* parameter files *}
{===>} parameter_infile_1="CNS_TOPPAR:protein_rep.param";
{===>} parameter_infile_2="CNS_TOPPAR:dna-rna_rep.param";
{===>} parameter_infile_3="CNS_TOPPAR:water_rep.param";
{===>} parameter_infile_4="CNS_TOPPAR:ion.param";
{===>} parameter_infile_5="CNS_TOPPAR:carbohydrate.param";
{===>} parameter_infile_6="";
{===>} parameter_infile_7="";
{===>} parameter_infile_8="";

{* molecular topology file: optional (leave blank for auto generation) *}
{* 
   Auto generation of the molecular topology from the coordinates should only 
   be used if:
   (1) Each distinct protein, DNA, or RNA chain must have a separate segid 
       (or chainid if the chainid is non-blank). 
   (2) Each contiguous protein, RNA, or RNA chain must not be disrupted by 
       other types of residues or ligands.  Rather, these other residues 
       should be listed after protein, RNA/DNA chains. 
   (3) Disulphides are automatically detected based on distances between the sulfur atoms
      (must be less than 3 A apart).
   (4) Broken protein/RNA/DNA chains without terminii must be more than 2.5 A apart to be recognized as such.
   (5) N-linked glycan links are automatically recognized if the bonded atoms are less than 2.5 A apart.
   (6) Automatic generation cannot be used with alternate conformations. 
   For ligands, the user must make suitable topology and parameter files.
   For non-standard covalent linkages, the custom patch file should be used.
   Alternatively, the generate.inp or generate_easy.inp task files
   can be used to generated the mtf prior to running this task file.
    *}
{===>} structure_infile="amy.mtf";

{* for auto generation: extra linkages and modifications by custom patches *}
{===>} patch_infile="";

{====================== 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;

{* anomalous f' f'' library file *}
{* If a file is not specified, no anomalous contribution will be included *}
{+ choice: "CNS_XRAYLIB:anom_cu.lib" "CNS_XRAYLIB:anom_mo.lib" "" user_file +}
{===>} anom_library="";

{* reflection files *}
{* specify non-anomalous reflection files before anomalous reflection files. *}
{* files must contain unique array names otherwise errors will occur *}
{===>} reflection_infile_1="amy.cv";
{===>} reflection_infile_2="";
{===>} reflection_infile_3="";
{===>} reflection_infile_4="";

{* reciprocal space array containing observed amplitudes: required *}
{===>} obs_f="fobs";

{* reciprocal space array containing sigma values for amplitudes: required *}
{===>} obs_sigf="sigma";

{* reciprocal space array containing test set for cross-validation: optional *}
{* required for the calculation of cross-validated sigmaA values *}
{===>} test_set="test";

{* flag for selection of test reflections: required for cross-validation *}
{* ie. reflections with the test set array equal to this number will be
       used for cross-validation, all other reflections form the working set *}
{===>} test_flag=1;

{* reciprocal space array containing observed intensities: optional *}
{===>} obs_i="";

{* reciprocal space array containing sigma values for intensities: optional *}
{===>} obs_sigi="";

{* output: complex reciprocal space array with model amplitudes and
           centroid phases from model, includes bulk solvent if used.
           note: do not use the reserved name "fcalc" *}
{===>} out_f="f_model";

{* output: real reciprocal space array with 
           figures-of-merit of model phases derived from sigmaA *}
{===>} out_fom="fom_model";

{* output: phase probability distribution *}
{* reciprocal space arrays with Hendrickson-Lattman 
   coefficients A,B,C,D from model derived from sigmaA *}
{+ table: rows=1 "HL coefficients" cols=4 "A" "B" "C" "D" +}
{===>} pa_out="pa_model";
{===>} pb_out="pb_model";
{===>} pc_out="pc_model";
{===>} pd_out="pd_model";

{* resolution limits *}
{+ table: rows=1 "resolution" cols=2 "lowest" "highest" +}
{===>} low_res=500.0;
{===>} high_res=2.0;

{* apply rejection criteria to amplitudes or intensities *}
{+ choice: "amplitude" "intensity" +}
{===>} obs_type="amplitude";

{* Observed data cutoff criteria: applied to amplitudes or intensities *}
{* reflections with magnitude(Obs)/sigma < cutoff are rejected. *}
{===>} sigma_cut=0.0;

{* rms outlier cutoff: applied to amplitudes or intensities *}
{* reflections with magnitude(Obs) > cutoff*rms(Obs) will be rejected *}
{===>} obs_rms=10000;

{=================== non-crystallographic symmetry ===================}

{* NCS-restraints/constraints file *}
{* see auxiliary/ncs.def *}
{===>} ncs_file="";

{============ overall B-factor and bulk solvent corrections ==========}

{* overall B-factor correction *}
{+ choice: "no" "isotropic" "anisotropic" +}
{===>} bscale="anisotropic";

{* bulk solvent correction *}
{* a mask is required around the molecule(s). The region
   outside this mask is the solvent region *}
{+ choice: true false +}
{===>} bulk_sol=true;

{* bulk solvent mask file *}
{* mask will be read from O type mask file if a name is given
   otherwise calculated from coordinates of selected atoms *}
{===>} bulk_mask_infile="";

{* automatic bulk solvent parameter optimization for e-density level sol_k (e/A^3) and B-factor sol_b (A^2) *}
{+ choice: true false +}
{===>} sol_auto=true;

{* fixed solvent parameters (used if the automatic option is turned off) *}
{+ table: rows=1 "bulk solvent" cols=2 "e-density level sol_k (e/A^3)" "B-factor sol_b (A^2) " +}
{===>} sol_k=0.3;
{===>} sol_b=50.0;

{* optional file with a listing of the results of the automatic bulk solvent optimization *}
{===>} sol_output="";

{* solvent mask parameters *}
{+ table: rows=1 "bulk solvent" cols=2 "probe radius (A) (usually set to 1)" "shrink radius (A) (usually set to 1)" +}
{===>} sol_rad=1.0;
{===>} sol_shrink=1.0;

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

{* select atoms to be included in calculating phases *}
{===>} atom_select=(known and not hydrogen);

{=================== phase calculation parameters ====================}

{* number of bins for sigmaa calculation *}
{* this will be determined automatically if a negative value is given
   otherwise the specified number of bins will be used *}
{===>} target_bins=-1;

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

{* calculate FOM using sigmaa *}
{* FOMs of 1 will be assigned if sigmaa is not used *}
{+ choice: true false +}
{===>} use_sigmaa=true;

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

{* output reflection file *}
{* note: only reflections will be written in specified resolution range 
         that satisfy data selection criteria *}
{===>} reflection_outfile="model_phase.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)

 if ( $log_level = verbose ) then
   set message=normal echo=on end
 else
   set message=off echo=off end
 end if

 if ( &BLANK%structure_infile = true ) then
 
    {- read topology files -}
    topology
     evaluate ($counter=1)
     evaluate ($done=false)
     while ( $done = false ) loop read
      if ( &exist_topology_infile_$counter = true ) then
         if ( &BLANK%topology_infile_$counter = false ) then
            @@&topology_infile_$counter
         end if
      else
        evaluate ($done=true)
      end if
      evaluate ($counter=$counter+1)
     end loop read
    end
    
    @CNS_XTALMODULE:mtfautogenerate (
                                  coordinate_infile=&coordinate_infile;
                                  convert=true;
                                  separate=true;
                                  atom_delete=(not known);
                                  hydrogen_flag=true;
                                  break_cutoff=2.5;
                                  disulphide_dist=3.0;
                                  carbo_dist=2.5;
                                  patch_infile=&patch_infile;
                                  O5_becomes="O";
                                 )

 else

   structure @&structure_infile end
   coordinates @&coordinate_infile

 end if

 {- read parameter files -}
 parameter
  evaluate ($counter=1)
  evaluate ($done=false)
  while ( $done = false ) loop read
   if ( &exist_parameter_infile_$counter = true ) then
      if ( &BLANK%parameter_infile_$counter = false ) then
         @@&parameter_infile_$counter
      end if
   else
     evaluate ($done=true)
   end if
   evaluate ($counter=$counter+1)
  end loop read
 end

 xray

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

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

   @CNS_XRAYLIB:scatter.lib

   binresolution &low_res &high_res
   mapresolution &high_res

   generate &low_res &high_res

   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

 end

 if ( &BLANK%anom_library = false ) then
   @@&anom_library
 else
   xray anomalous=? end
   if ( $result = true ) then
     display Warning: no anomalous library has been specified
     display          no anomalous contribution will used in refinement
   end if
 end if

 evaluate ($obs_i=&obs_i)
 evaluate ($obs_sigi=&obs_sigi)

 {- BEGIN MODIFICATION -}
 evaluate ($array_name=&STRIP%out_f)
 if ($array_name = "FCALC") then
    display 
    display  error: do not use reserved name "fcalc" for output array
    display  aborting script
    display 
    abort
 end if
 {- END MODIFICATION -}

 xray

   if ( &BLANK%obs_f = true ) then
      display 
      display  *********************************************************
      display  Error: required observed amplitude array is not specified
      display  *********************************************************
      display
      abort
   else
      query name=&STRIP%obs_f domain=reciprocal end
      if ( $object_exist = false ) then
          display 
          display  **************************************************************
          display  Error: required observed amplitude array &obs_f does not exist
          display  **************************************************************
          display
          abort
       end if
    {- note: this array can be of any type -}
    end if

    if ( &BLANK%obs_sigf = true ) then
      display 
      display  *****************************************************
      display  Error: required observed sigma array is not specified
      display  *****************************************************
      display
      abort
    else
       query name=&STRIP%obs_sigf domain=reciprocal end
       if ( $object_exist = false ) then
          display 
          display  *************************************************************
          display  Error: required observed sigma array &obs_sigf does not exist
          display  *************************************************************
          display
          abort
       end if
       if ( $object_type # "REAL" ) then
          display 
          display  **********************************************************************
          display  Error: required observed sigma array &obs_sigf has the wrong data type
          display  **********************************************************************
          display
          abort
       end if
    end if

    if ( &BLANK%test_set = false ) then
      query name=&STRIP%test_set domain=reciprocal end
      if ( $object_exist = false ) then
        display 
        display  **********************************************
        display  Error: test set array &test_set does not exist
        display  **********************************************
        display
        abort
      end if
      {- test set array can be integer or real -}
      if ( $object_type = "COMPLEX" ) then
        display 
        display  *******************************************************
        display  Error: test set array &test_set has the wrong data type
        display  *******************************************************
        display
        abort
      end if
   end if
   
   query name=fcalc domain=reciprocal end
   if ( $object_exist = false ) then
      declare name=fcalc domain=reciprocal type=complex end
   end if
   
   declare name=fbulk domain=reciprocal type=complex end
   do (fbulk=0) ( all )

   query name=&STRIP%obs_f domain=reciprocal end
   declare name=fobs_orig domain=reciprocal type=$object_type end
   declare name=sigma_orig domain=reciprocal type=real end

   do (fobs_orig=&STRIP%obs_f) (all)
   do (sigma_orig=&STRIP%obs_sigf) (all)

   if ( &BLANK%obs_i = false ) then
     query name=&STRIP%obs_i domain=reciprocal end
     declare name=iobs_orig domain=reciprocal type=$object_type end
     declare name=sigi_orig domain=reciprocal type=real end
     do (iobs_orig=&STRIP%obs_i) (all)
     do (sigi_orig=&STRIP%obs_sigi) (all)
   end if

   if ( &obs_type = "intensity" ) then
     if ( &BLANK%obs_i = true ) then
       display  Error: observed intensity array is undefined
       display         aborting script
       abort
     end if
     evaluate ($reject_obs=&obs_i)
     evaluate ($reject_sig=&obs_sigi)
   else
     evaluate ($reject_obs=&obs_f)
     evaluate ($reject_sig=&obs_sigf)
   end if

   declare name=ref_active domain=reciprocal type=integer end
   declare name=tst_active domain=reciprocal type=integer end

   do (ref_active=0) ( all )
   do (ref_active=1) ( ( $STRIP%reject_sig # 0 ) and
                      ( &low_res >= d >= &high_res ) )

   statistics overall
     completeness
     selection=( ref_active=1 )
   end
   evaluate ($total_compl=$expression1)

   show sum(1) ( ref_active=1 )
   evaluate ($total_read=$select)
   evaluate ($total_theor=int(1./$total_compl * $total_read))

   show rms (amplitude($STRIP%reject_obs)) ( ref_active=1 )
   evaluate ($obs_high=$result*&obs_rms)
   show min (amplitude($STRIP%reject_obs)) ( ref_active=1 )
   evaluate ($obs_low=$result)

   do (ref_active=0) ( all )
   do (ref_active=1)
      ( ( amplitude($STRIP%reject_obs) > &sigma_cut*$STRIP%reject_sig ) and
        ( $STRIP%reject_sig # 0 ) and
        ( $obs_low <= amplitude($STRIP%reject_obs) <= $obs_high ) and
        ( &low_res >= d >= &high_res ) )

   do (tst_active=0) (all)
   if ( &BLANK%test_set = false ) then
     do (tst_active=1) (ref_active=1 and &STRIP%test_set=&test_flag)
   end if

   show sum(1) ( ref_active=1 and tst_active=0 )
   evaluate ($total_work=$select)
   show sum(1) ( ref_active=1 and tst_active=1 )
   evaluate ($total_test=$select)
   evaluate ($total_used=$total_work+$total_test)

   evaluate ($unobserved=$total_theor-$total_read)
   evaluate ($rejected=$total_read-$total_used)
   evaluate ($per_unobs=100*($unobserved/$total_theor))
   evaluate ($per_reject=100*($rejected/$total_theor))
   evaluate ($per_used=100*($total_used/$total_theor))
   evaluate ($per_work=100*($total_work/$total_theor))
   evaluate ($per_test=100*($total_test/$total_theor))

   associate fcalc ( &atom_select )

   tselection=( ref_active=1 )

   cvselection=( tst_active=1 )

   method=FFT          
   
 {- MODIFIED 2/15/06 -}
 end
 

 show min ( b ) ( &atom_select )
 evaluate ($b_min=$result)
 @@CNS_XTALMODULE:fft_parameter_check ( 
                             d_min=&high_res; 
                             b_min=$b_min;
                             grid=auto;
                             fft_memory=&fft_memory;
                             fft_grid=$fft_grid;   
                             fft_b_add=$fft_b_add; 
                             fft_elim=$fft_elim; 
                                      )
                            

 {- END MODIFICATION -}
             

 if ( &BLANK%ncs_file = false ) then
    inline @&ncs_file
 end if

 xray 
   do (&STRIP%obs_f=fobs_orig) (all)
   do (&STRIP%obs_sigf=sigma_orig) (all)
   if ( &BLANK%obs_i = false ) then
     do (&STRIP%obs_i=iobs_orig) (all)
     do (&STRIP%obs_sigi=sigi_orig) (all)
   end if
 end

 xray
    predict
      mode=reciprocal
      to=fcalc
      selection=(ref_active=1)
      atomselection=( &atom_select ) 
    end
 end
                             
 {- BEGIN MODIFICATION -}
 @CNS_XTALMODULE:scale_and_solvent_grid_search (
                             bscale=&bscale;
                             sel=( ref_active=1 );
                             sel_test=( tst_active=1 );
                             atom_select=( &atom_select );
                             bulk_sol=&bulk_sol;
                             bulk_mask=&bulk_mask_infile;
                             bulk_atoms=( &atom_select );
                             
                             sol_auto=&sol_auto;
                             sol_k=&sol_k;
                             sol_b=&sol_b;
                             sol_rad=&sol_rad;
                             sol_shrink=&sol_shrink;
  
                             fcalc=fcalc;
                             obs_f=&STRIP%obs_f;
                             obs_sigf=&STRIP%obs_sigf;
                             obs_i=$STRIP%obs_i;
                             obs_sigi=$STRIP%obs_sigi;                             
                             fpart=fbulk;
                             
!
! Begin modification (6/28/06)                             
                             Baniso_11=$Baniso_11;
                             Baniso_22=$Baniso_22;
                             Baniso_33=$Baniso_33;
                             Baniso_12=$Baniso_12;
                             Baniso_13=$Baniso_13;
                             Baniso_23=$Baniso_23;
                             Biso=$Biso_model;
! End modification
! 
                             
                             sol_k_best=$sol_k_ref;
                             sol_b_best=$sol_b_ref;
			     solrad_best=$solrad_best;
			     shrink_best=$shrink_best;
                             
                             b=b;

                             low_b_flag=$low_b_flag;
                            
                             sol_output=&sol_output;
                             
                             )
 xray
   @@CNS_XTALMODULE:calculate_r (
                                 fobs=&STRIP%obs_f;
                                 fcalc=fcalc;
                                 fpart=fbulk;
                                 sel=( ref_active=1 );
                                 sel_test=( tst_active=1 );
                                 print=true;
                                 output=OUTPUT;
                                 r=$start_r;
                                 test_r=$start_test_r;)
 end

 {- check the gridding again since the minimum B-factor may have changed -}
 show min ( b ) ( &atom_select )
 evaluate ($b_min=$result)
 @@CNS_XTALMODULE:fft_parameter_check ( 
                             d_min=&high_res; 
                             b_min=$b_min;
                             grid=auto;
                             fft_memory=&fft_memory;
                             fft_grid=$fft_grid;   
                             fft_b_add=$fft_b_add; 
                             fft_elim=$fft_elim; 
                                      )
{- END MODIFICATION -}                         

 xray
   declare name=m          domain=reciprocal type=complex end
   declare name=mod_fom    domain=reciprocal type=real end
   declare name=mod_x      domain=reciprocal type=real end
   declare name=mod_pa     domain=reciprocal type=real end
   declare name=mod_pb     domain=reciprocal type=real end
   declare name=mod_pc     domain=reciprocal type=real end
   declare name=mod_pd     domain=reciprocal type=real end
   declare name=mod_dd     domain=reciprocal type=real end

   do (mod_pa=0) (all)
   do (mod_pb=0) (all)
   do (mod_pc=0) (all)
   do (mod_pd=0) (all)

   if ( &use_sigmaa = true ) then

     @CNS_XTALMODULE:fomsigmaacv ( sig_sigacv=0.07;
                                   mbins=&target_bins;
                                   statistics=true;
                                   fobs=&STRIP%obs_f;
                                   fcalc=fcalc;
                                   fpart=fbulk;
                                   test=tst_active;
                                   sel=(ref_active=1);
                                   sel_test=(tst_active=1);
                                   fom=mod_fom;
                                   x=mod_x;
                                   pa=mod_pa;
                                   pb=mod_pb;
                                   pc=mod_pc;
                                   pd=mod_pd; 
                                   dd=mod_dd; )

   else

     do (mod_pa=10*cos(phase(fcalc+fbulk))) ( ref_active=1 and centric )
     do (mod_pb=10*sin(phase(fcalc+fbulk))) ( ref_active=1 and centric )
     do (mod_pa=1000*cos(phase(fcalc+fbulk))) ( ref_active=1 and acentric )
     do (mod_pb=1000*sin(phase(fcalc+fbulk))) ( ref_active=1 and acentric )
     do (mod_pc=0) ( ref_active=1 )
     do (mod_pd=0) ( ref_active=1 )

     do (mod_fom=1) ( ref_active=1 )

   end if

 end

 xray

   if ( &BLANK%pa_out = true ) then
     display 
     display  *********************************************************
     display  Error: output Hendrickson-Lattman array A not specified
     display  *********************************************************
     display
     abort
   elseif ( &BLANK%pb_out = true ) then
     display 
     display  *********************************************************
     display  Error: output Hendrickson-Lattman array B not specified
     display  *********************************************************
     display
     abort
   elseif ( &BLANK%pc_out = true ) then
     display 
     display  *********************************************************
     display  Error: output Hendrickson-Lattman array C not specified
     display  *********************************************************
     display
     abort
   elseif ( &BLANK%pd_out = true ) then
     display 
     display  *********************************************************
     display  Error: output Hendrickson-Lattman array D not specified
     display  *********************************************************
     display
     abort
   elseif ( &BLANK%out_fom = true ) then
     display 
     display  *************************************************
     display  Error: output figure of merit array not specified
     display  *************************************************
     display
     abort
   elseif ( &BLANK%out_f = true ) then
     display 
     display  **********************************************************
     display  Error: output complex structure factor array not specified
     display  **********************************************************
     display
     abort
   end if
   	
   {- the output arrays must not already exist -}

   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

   declare name=&STRIP%out_fom domain=reciprocal type=real end
   declare name=&STRIP%out_f domain=reciprocal type=complex end

   {- reinstall original (unscaled) amplitude, sigma, and intensity arrays -}
   do (&STRIP%obs_f=fobs_orig) (all)
   do (&STRIP%obs_sigf=sigma_orig) (all)
   if ( &BLANK%obs_i = false ) then
     do (&STRIP%obs_i=iobs_orig) (all)
     do (&STRIP%obs_sigi=sigi_orig) (all)
   end if
 end

 xray
    do (&STRIP%pa_out=mod_pa) (all)
    do (&STRIP%pb_out=mod_pb) (all)
    do (&STRIP%pc_out=mod_pc) (all)
    do (&STRIP%pd_out=mod_pd) (all)
     
    do (&STRIP%out_f=fcalc+fbulk) (all)
      
    do (&STRIP%out_fom=mod_fom) (all)
 end

 xray
 
   undeclare name=fobs_orig  domain=reciprocal end
   undeclare name=sigma_orig domain=reciprocal end
   undeclare name=iobs_orig  domain=reciprocal end
   undeclare name=sigi_orig  domain=reciprocal end
   undeclare name=mod_pa     domain=reciprocal end
   undeclare name=mod_pb     domain=reciprocal end
   undeclare name=mod_pc     domain=reciprocal end
   undeclare name=mod_pd     domain=reciprocal end
   undeclare name=mod_fom    domain=reciprocal end
   undeclare name=mod_x      domain=reciprocal end
   undeclare name=mod_dd     domain=reciprocal end
   undeclare name=ref_active domain=reciprocal end
   undeclare name=tst_active domain=reciprocal end
   undeclare name=m          domain=reciprocal end
   undeclare name=fcalc      domain=reciprocal end    
    
    
   declare ? end
   
   {- modification: only selected reflections will be written -}

   set display=&reflection_outfile end

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

  xray
   
   if ( &merge_inout = true ) then  
      write reflection
        output=&reflection_outfile
        sele=( $STRIP%reject_sig # 0 and
                       ( &low_res >= d >= &high_res ) )
      end
   else
      write reflection
        &STRIP%out_f 
        &STRIP%out_fom 
        &STRIP%pa_out &STRIP%pb_out &STRIP%pc_out &STRIP%pd_out 
        output=&reflection_outfile
        sele=( $STRIP%reject_sig # 0 and
                       ( &low_res >= d >= &high_res ) )
      end
    end if
    
    set display=OUTPUT end

 end
 
 stop