{+ file: optimize_average.inp +}
{+ directory: xtal_refine +}
{+ description: Finds optimum combination of coordinate files which produce 
                minimum free R-value from averaged structure factors +}
{+ authors: Axel T. Brunger, Luke M. Rice and Paul D. Adams +}
{+ copyright: Yale University +}

{+ reference: A.T. Brunger, The Free R Value: a Novel Statistical
              Quantity for Assessing the Accuracy of Crystal Structures,
              Nature 355, 472-474 (1992) +}
{+ reference: G.J. Kleywegt and A.T. Brunger, Checking your 
              imagination: Applications of the free R value, 
              Structure 4, 897-904 (1996) +}
{+ reference: A.T. Brunger, P.D. Adams and L.M. Rice, New applications
              of simulated annealing in X-ray crystallography and
              solution NMR, Structure 5, 325-336 (1997) +}

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

{======================= 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_1="amy.pdb";
{===>} coordinate_infile_2="";
{===>} coordinate_infile_3="";
{===>} coordinate_infile_4="";
{===>} coordinate_infile_5="";
{===>} coordinate_infile_6="";
{===>} coordinate_infile_7="";
{===>} coordinate_infile_8="";
{===>} coordinate_infile_9="";
{===>} coordinate_infile_10="";


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

{* 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: required *}
{===>} test_set="test";

{* number for selection of test reflections: required *}
{* 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 weighting scheme for observed
   amplitudes: optional *}
{* only used for the "residual" and "vector" targets - this will
   default to a constant value of 1 if array is not present *}
{===>} obs_w="";

{* reciprocal space array containing observed intensities: optional *}
{* required for the "mli" target *}
{===>} obs_i="";

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

{* reciprocal space arrays with experimental phase probability
   distribution: optional *}
{* Hendrickson-Lattman coefficients A,B,C,D *}
{* required for the "mlhl" target *}
{+ table: rows=1 "HL coefficients" cols=4 "A" "B" "C" "D" +}
{===>} obs_pa="";
{===>} obs_pb="";
{===>} obs_pc="";
{===>} obs_pd="";

{* complex reciprocal space array containing experimental phases: optional *}
{* required for the "mixed" and "vector" targets *}
{===>} obs_phase="";

{* reciprocal space array containing experimental figures of merit: optional *}
{* required for the "mixed" target *}
{===>} obs_fom="";

{* resolution limits to be used in refinement *}
{* the full resolution range of observed data should be used in refinement.
   A bulk solvent correction should be applied to allow the use of low
   resolution terms. If no bulk solvent correction is applied, data must
   be truncated at a lower resolution limit of between 8 and 6 Angstrom. *}
{+ 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;

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

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

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

{============ 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 search *}
{+ choice: true false +}
{===>} sol_auto=true;

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

{* fixed solvent mask parameters if the automatic option is not used *}
{+ table: rows=1 "bulk solvent" cols=2 "probe radius (A)" "shrink radius (A)" +}
{===>} sol_rad=1.0;
{===>} sol_shrink=1.0;

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

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

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

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

{* output listing file *}
{===>} list_outfile="optimize_average.list";

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

 ) {- end block parameter definition -}

 checkversion 1.2

 evaluate ($log_level=quiet)

 structure @&structure_infile end

 coordinates @&coordinate_infile_1

 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

 end

 if ( &BLANK%anom_library = false ) then
   @@&anom_library
 else
   set echo=off end
   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
   set echo=on end
 end if

 evaluate ($obs_i=&obs_i)
 evaluate ($obs_sigi=&obs_sigi)
 evaluate ($obs_w=&obs_w)
 
 {- check only required arrays by setting the target to "" -}
 xray
   @@CNS_XTALMODULE:checkrefinput (
                                  reftarget=&reftarget;
                                  obs_f=&obs_f;
                                  obs_sigf=&obs_sigf;
                                  test_set=&test_set;
                                  obs_pa=&obs_pa;
                                  obs_pb=&obs_pb;
                                  obs_pc=&obs_pc;
                                  obs_pd=&obs_pd;
                                  obs_phase=&obs_phase;
                                  obs_fom=&obs_fom;
                                  obs_w=$obs_w;
                                  obs_i=$obs_i;
                                  obs_sigi=$obs_sigi;
                                  )

   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 )

   binresolution &low_res &high_res
   mapresolution &high_res

   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)
     show min (amplitude(&STRIP%obs_i)) (all)
     evaluate ($obs_lower_limit=$result-0.1)
   else
     evaluate ($reject_obs=&obs_f)
     evaluate ($reject_sig=&obs_sigf)
     evaluate ($obs_lower_limit=0)
   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) ( ( amplitude($STRIP%reject_obs) > $obs_lower_limit ) 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)
   do (tst_active=1) (ref_active=1 and &STRIP%test_set=&test_flag)

   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; 
                                      )
                            
 xray
   tolerance=0.0 lookup=false
 {- END MODIFICATION -}

 end                  

 xray
   @@CNS_XTALMODULE:checkcv (active=ref_active;
                             test=tst_active;)
 end

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

 evaluate ($struct=1)
 evaluate ($done=false)
 while ( $done = false ) loop struct
   if ( &exist_coordinate_infile_$struct = true ) then
     if ( &BLANK%coordinate_infile_$struct = true ) then
       evaluate ($done=true)
       evaluate ($struct=$struct-1)
     else
       evaluate ($struct=$struct+1)
     end if       
   else
     evaluate ($done=true)
     evaluate ($struct=$struct-1)
   end if
 end loop struct

 set display=&list_outfile end

 display >>> resolution: &low_res - &high_res A
 display >>> sg= &STRIP%sg a= &a b= &b c= &c alpha= &alpha beta= &beta gamma= &gamma
 if ( &BLANK%parameter_infile_1 = false ) then
    display >>> parameter file 1  : &STRIP%parameter_infile_1
 end if
 if ( &BLANK%parameter_infile_2 = false ) then
    display >>> parameter file 2  : &STRIP%parameter_infile_2
 end if
 if ( &BLANK%parameter_infile_3 = false ) then
    display >>> parameter file 3  : &STRIP%parameter_infile_3
 end if
 if ( &BLANK%parameter_infile_4 = false ) then
    display >>> parameter file 4  : &STRIP%parameter_infile_4
 end if
 if ( &BLANK%parameter_infile_5 = false ) then
    display >>> parameter file 5  : &STRIP%parameter_infile_5
 end if
 display >>> molecular structure file: &STRIP%structure_infile
 if ( &BLANK%anom_library = false ) then
   display >>> anomalous f' f'' library: &STRIP%anom_library
 end if
 if ( &BLANK%reflection_infile_1 = false ) then
   display >>> reflection file= &STRIP%reflection_infile_1
 end if
 if ( &BLANK%reflection_infile_2 = false ) then
   display >>> reflection file= &STRIP%reflection_infile_2
 end if
 if ( &BLANK%reflection_infile_3 = false ) then
   display >>> reflection file= &STRIP%reflection_infile_3
 end if
 if ( &BLANK%ncs_infile = false ) then
   display >>> ncs= &STRIP%ncs_type  ncs file= &STRIP%ncs_infile
 else
   display >>> ncs= none
 end if
 if ( &obs_type = "intensity" ) then
   display >>> reflections with Iobs/sigma_I < &sigma_cut rejected
   display >>> reflections with Iobs > &obs_rms * rms(Iobs) rejected
 else
   display >>> reflections with |Fobs|/sigma_F < &sigma_cut rejected
   display >>> reflections with |Fobs| > &obs_rms * rms(Fobs) rejected
 end if
 xray anomalous=? end
 if ( $result = true ) then
   display >>> anomalous diffraction data was input
 end if

 {- MODIFIED 2/15/06 -}
 display >>> fft gridding factor = $fft_grid, B factor offset = $fft_b_add A^2, Elimit = $fft_elim
 {- END MODIFICATION -}

 display >>> theoretical total number of refl. in resol. range:    $total_theor[I6] ( 100.0 % )
 display >>> number of unobserved reflections (no entry or |F|=0): $unobserved[I6] ( $per_unobs[f5.1] % )
 display >>> number of reflections rejected:                       $rejected[I6] ( $per_reject[f5.1] % )
 display >>> total number of reflections used:                     $total_used[I6] ( $per_used[f5.1] % )
 display >>> number of reflections in working set:                 $total_work[I6] ( $per_work[f5.1] % )
 display >>> number of reflections in test set:                    $total_test[I6] ( $per_test[f5.1] % )

 display ==============================================================================

 set display=OUTPUT end

 evaluate ($cycle=1)

 xray
   declare name=fcave domain=reciprocal type=complex end
   declare name=fpave domain=reciprocal type=complex end
   query name=&STRIP%obs_f domain=reciprocal end
   declare name=foave domain=reciprocal type=$object_type end
   do (fcave=0) (all)
   do (fpave=0) (all)
   do (foave=0) (all)
 end

 while ( $cycle <= $struct ) loop main

   coord init end
   coord @@&coordinate_infile_$cycle

   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
   
   set display=&list_outfile end

   display >>> structure: &STRIP%coordinate_infile_$cycle

   if ( &bscale # "no" ) then
     if ( $low_b_flag = true ) then
       display >>> warning: B-correction gave atomic B-values less than zero
       display >>>          they have been reset to zero
     end if
   end if
   ! 
   ! Begin modification (6/28/06)
   if ( &bscale = "anisotropic" ) then
    display >>> Anisotropic B-factor tensor Ucart of atomic model without isotropic component :
    display >>>   B11=$Baniso_11[f8.3] B22=$Baniso_22[f8.3] B33=$Baniso_33[f8.3]
    display >>>   B12=$Baniso_12[f8.3] B13=$Baniso_13[f8.3] B23=$Baniso_23[f8.3]
    display >>> Isotropic component added to coordinate array B: $Biso_model[f8.3]
   elseif ( &bscale = "isotropic" ) then
    display >>> B-factor applied to coordinate array B: $Biso_model[f8.3]
   else
    display >>> initial B-factor correction: none
   end if
   ! End modification
   !

   {- MODIFIED 5/18/05 -}
   if ( &bulk_sol = true ) then 
     display >>> bulk solvent: probe radius=$solrad_best, shrink value=$solrad_best
     display >>> bulk solvent: density level= $sol_k_ref e/A^3, B-factor= $sol_b_ref A^2
   else
     display >>> bulk solvent: false
   end if
   {- END MODIFICATION -}

   set display=OUTPUT end

   evaluate ($farray="fc_" + encode($cycle))
   evaluate ($fparray="fp_" + encode($cycle))
   evaluate ($foarray="fo_" + encode($cycle))

   xray
     declare name=$STRIP%farray domain=reciprocal type=complex end
     declare name=$STRIP%fparray domain=reciprocal type=complex end
     query name=&STRIP%obs_f domain=reciprocal end
     declare name=$STRIP%foarray domain=reciprocal type=$object_type end
     @@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=$r_$cycle;
                                  test_r=$test_r_$cycle;)

     do ($STRIP%farray=fcalc) (all)
     do ($STRIP%fparray=fbulk) (all)
     do ($STRIP%foarray=&STRIP%obs_f) (all)
   end

   set display=&list_outfile end

   display >>> free R-value= $test_r_$cycle R-value= $r_$cycle
   display ==============================================================================

   set display=OUTPUT end

   evaluate ($used_$cycle=false)

   evaluate ($cycle=$cycle+1)

 end loop main

 evaluate ($cycle_i=1)
 while ( $cycle_i <= $struct ) loop i
   evaluate ($cycle_j=1)
   evaluate ($best_j=1)
   while ( $cycle_j <= $struct ) loop j
     if ( $used_$cycle_j = false ) then
       if ( $test_r_$cycle_j < $best_j ) then
          evaluate ($best_j=$test_r_$cycle_j)
          evaluate ($number=$cycle_j)
       end if
     end if
     evaluate ($cycle_j=$cycle_j+1)
   end loop j
   evaluate ($order_$cycle_i=$number)
   evaluate ($used_$number=true)
   evaluate ($cycle_i=$cycle_i+1)
 end loop i

 evaluate ($cycle=1)
 
 set display=&list_outfile end

 display >>> ranking of structures by free R-value (lowest to highest)

 while ( $cycle <= $struct ) loop struct

   evaluate ($number=$order_$cycle)
   display   structure= &STRIP%coordinate_infile_$number \
R= $r_$number free-R= $test_r_$number

   evaluate ($cycle=$cycle+1)

 end loop struct

 display ==============================================================================
 display >>> free R-value and R-value after averaging structure factors
 evaluate ($cycle=1)

 set display=OUTPUT end

 evaluate ($best_free=1)
 
 while ( $cycle <= $struct ) loop struct

   evaluate ($add=1)
   xray
     do (fcave=0) (all)
     do (fpave=0) (all)
     do (foave=0) (all)
   end
   while ( $add <= $cycle ) loop average

     evaluate ($number=$order_$add)
     xray
       evaluate ($fcarray="fc_" + encode($number))
       evaluate ($fparray="fp_" + encode($number))
       evaluate ($foarray="fo_" + encode($number))
       do (fcave=fcave+$STRIP%fcarray) (all)
       do (fpave=fpave+$STRIP%fparray) (all)
       do (foave=foave+$STRIP%foarray) (all)
     end

     evaluate ($add=$add+1)

   end loop average

   xray
     do (fcave=fcave/$cycle) (all)
     do (fpave=fpave/$cycle) (all)
     do (foave=foave/$cycle) (all)
   end

   xray
     @@CNS_XTALMODULE:calculate_r (fobs=foave;
                                  fcalc=fcave;
                                  fpart=fpave;
                                  sel=(ref_active=1);
                                  sel_test=(tst_active=1);
                                  print=true;
                                  output=OUTPUT;
                                  r=$r;
                                  test_r=$test_r;)
   end
   
   if ( $test_r < $best_free ) then
     evaluate ($best_free=$test_r)
     evaluate ($best_num=$cycle)
   end if

   set display=&list_outfile end

   display   top $cycle structures  R-value= $r free R-value= $test_r

   set display=OUTPUT end

   evaluate ($cycle=$cycle+1)

 end loop struct
 
 set display=&list_outfile end

 display ==============================================================================
 display
 display ----> lowest free R-value from averaging top $best_num structures
 display 
 display >>> these structures are

 evaluate ($cycle=1)
 
 while ( $cycle <= $best_num ) loop struct

   evaluate ($number=$order_$cycle)
   display   &STRIP%coordinate_infile_$number

   evaluate ($cycle=$cycle+1)

 end loop struct

 display ==============================================================================

 set display=OUTPUT end

 stop