{+ file: water_pick_twin.inp +}
{+ directory: xtal_twin +}
{+ description: pick water molecules in electron density map 
                using phase information from a model, for data
                with hemihedral twinning +}
{+ comments:
             choice of coefficients:
               (u m|Fo| e^(i phi_calc)) - (v D|Fc| e^(i phi_calc))
               (u  |Fo| e^(i phi_calc)) - (v k|Fc| e^(i phi_calc))
                d(target)/dFc
             where Fo is the detwinned observed structure factor,
                   Fc is the calculated structure factor and
                   m and D are derived from sigmaa +}
{+ authors: Axel T. Brunger and Paul D. Adams +}
{+ copyright: Yale University +}

{+ reference: T.O. Yeates, Detecting and Overcoming Crystal Twinning,
              Methods in Enzymology 276, 344-358 (1997) +}
{+ reference: R.J. Read, Improved Fourier coefficients for maps using
              phases from partial structures with errors.
              Acta Cryst. A42, 140-149 (1986) +}
{+ reference: G.J. Kleywegt and A.T. Brunger, Checking your imagination:
              Applications of the free R value, Structure 4,
              897-904 (1996) +}

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

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

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

{====================== crystallographic data ========================}

{* space group *}
{* use International Table conventions with subscripts substituted
   by parenthesis *}
{===>} sg="R3";

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

{* 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="porin.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 for calculation of cross-validated sigmaA values *}
{===>} test_set="test";

{* number 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 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 for data included in map calculation *}
{* all data available should be included in the map calculation *}
{+ table: rows=1 "resolution" cols=2 "lowest" "highest" +}
{===>} low_res=500.0;
{===>} high_res=2.25;

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

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

{======================= twinning parameters =========================}

{* twinning operation *}
{===>} twin_oper="h,-h-k,-l";

{* twinning fraction *}
{===>} twin_frac=0.304;

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

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

{==================== map generation parameters ======================}

{* maps are calculated u*Fo - v*Fc *}
{* eg. 2fo-fc map -> u=2 and v=1 or
        fo-fc map -> u=1 and v=1 *}

{* specify u *}
{===>} u=1;

{* specify v *}
{===>} v=1;

{* type of map *}
{+ list:   sigmaa: (u m|Fo| - v D|Fc|)^exp(i phi_calc)
                    m and D calculated from sigmaa
       unweighted: (u |Fo| - v k|Fc|)^exp(i phi_calc)
                    no figure-of-merit weighting
         gradient: d(target)/dFc
                   gradient of the current crystallographic target wrt Fc +}
{+ choice: "sigmaa" "unweighted" "gradient" +}
{===>} map_type="sigmaa";

{* refinement target *}
{+ list: twin_lsq: least squares residual for hemihedral twinning +}
{+ choice: "twin_lsq" +}
{===>} reftarget="twin_lsq";

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

{* use model amplitudes for unmeasured data *}
{* this will not be applied to gradient or difference maps *}
{+ choice: true false +}
{===>} fill_in=false;

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

{* map grid size: dmin*grid *}
{===>} grid=0.25;

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

{===================== water picking parameters ======================}

{* minimum peak height that will be allowed *}
{+ list: if map scaling is true: in sigma units
         if map scaling is false: in electron density (observed) units +}
{===>} peak=3.0;

{* maximum number of peaks searched in map *}
{===>} npeaks=400;

{* residue name for picked waters *}
{+ choice: "TIP" "TIP3" "WAT" "HOH" +}
{===>} wat_resname="TIP";

{* segid for newly picked waters (must be 4 characters or less) *}
{* do not use the segid PEAK - this is needed elsewhere in 
   the task file *}
{===>} wat_segid="S";

{* first residue number for waters *}
{* if waters are already present in the segid selected above then
   numbering will begin after the last residue number in the segment *}
{===>} wat_number=1;

{=================== water refinement parameters =====================}

{* refine position and B-factors for newly picked waters *}
{+ choice: true false +}
{===>} refine=true;

{* number of steps of conjugate gradient minimization for coordinates *}
{===>} minimize_nstep=20;

{* number of steps of conjugate gradient minimization for B-factors *}
{===>} bfactor_nstep=20;

{* B-factor limits *}
{+ table: rows=1 "B-factor" cols=2 "minimum" "maximum" +}
{===>} bmin=1;
{===>} bmax=200;

{==================== water deleting parameters ======================}

{* selection criteria *}

{* minimum distance between water and any atom *}
{===>} min=2.6;
{* maximum distance between water and any atom *}
{===>} max=4.0;

{* minimum hydrogen bonding distance between water and O or N *}
{===>} hmin=2.0;
{* maximum hydrogen bonding distance between water and O or N *}
{===>} hmax=3.2;

{* distance based peak deletion criteria *}
{+ list: if "none", all peaks greater than the minimum peak height and 
                    closer than the maximum distance are kept 
         if "distance", peaks greater than the minimum distance and 
                        less than maximum distance from any atom are kept 
         if "hbond", the above distance criterion is applied except that
                     peaks further than the minimum hydrogen bonding distance
                     from oxygen or nitrogen atoms are also kept +}
{+ choice: "none" "distance" "hbond" +}
{===>} peak_select="hbond";

{* hydrogen bonding based peak deletion criteria *}
{* delete all peaks not forming hydrogen bonds with the starting model.
   Only peaks closer than the maximum hbond distance and further away 
   than the minimum hbond distance from oxygen or nitrogen atoms will
   be kept *}
{+ choice: true false +}
{===>} peak_hbond=true;

{* cutoff distance (Angstroms) to identify atoms on special positions *}
{* any atom which is within this distance of itself after application
   of symmetry will be considered to be on a special position *}
{===>} special_dist=0.25;

{* keep waters on special positions *}
{* peaks which coincide with special positions can be 
   retained or deleted *}
{+ choice: true false +}
{===>} peak_special=false;

{* maximum B-factor (A^2) for picked waters *}
{* if coordinate and B-factor refinement are applied, then waters
   with a B-factor greater than this value will be deleted *}
{===>} bfactor_max=50;

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

{* merge picked waters with input coordinates *}
{* if true then output molecular structure and coordinate files
   will contain original and new atoms *}
{+ choice: true false +}
{===>} merge_coord=true;

{* output molecular structure file *}
{===>} structure_outfile="water_pick_twin.mtf";

{* output coordinate file *}
{===>} coordinate_outfile="water_pick_twin.pdb";

{* format output coordinates for use in o *}
{* if false then the default CNS output coordinate format will be used *}
{+ choice: true false +}
{===>} pdb_o_format=true;

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

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

   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

 {- copy define parameters of optional arrays into symbols so 
    we can redefine them -}

 evaluate ($obs_i=&obs_i)
 evaluate ($obs_sigi=&obs_sigi)
 evaluate ($obs_w=&obs_w)
 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 )

   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
                    ( $obs_low <= amplitude(remap[&STRIP%twin_oper]($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 8/01/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
 {- END MODIFICATION -}

   tolerance=0.0 lookup=false

 end

 if ( &map_type = "observed" ) then
   evalaute ($test_hl=true)
 elseif ( &map_type = "combined" ) then
   evalaute ($test_hl=true)
 else
   evalaute ($test_hl=false)
 end if

 if ( $test_hl = true ) then
   xray
     @@CNS_XTALMODULE:check_abcd (pa=&obs_pa;
                                 pb=&obs_pb;
                                 pc=&obs_pc;
                                 pd=&obs_pd;)
   end
 end if

 evaluate ($strict_ncs=false)
 if ( &BLANK%ncs_infile = false ) then
   inline @&ncs_infile
   if ( &ncs_type = "strict" ) then
     evaluate ($strict_ncs=true)
   end if
 end if

 xray
   declare name=dtarg domain=reciprocal type=complex end
   declare name=total domain=reciprocal type=complex end
   declare name=fmap domain=reciprocal type=complex end
   predict
     mode=reciprocal
     to=fcalc
     selection=( &low_res >= d >= &high_res )
     atomselection=( &atom_select ) 
   end
 end

 {- BEGIN MODIFICATION 8/01/06 -}
 @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;

                             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;
                             
                             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;
                             
                             )
 
 {- 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
   if ( &map_type = "gradient" ) then
       @@CNS_XTALMODULE:refinementtarget_twin (target=&reftarget;
                                               sig_sigacv=0.07;
                                               mbins=&target_bins;
                                               fobs=&STRIP%obs_f;
                                               sigma=&STRIP%obs_sigf;
                                               weight=$STRIP%obs_w;
                                               iobs=$STRIP%obs_i;
                                               sigi=$STRIP%obs_sigi;
                                               test=tst_active;
                                               fcalc=fcalc;
                                               fpart=fbulk;
                                               twin_oper=&STRIP%twin_oper;
                                               twin_frac=&twin_frac;
                                               pa=&STRIP%obs_pa;
                                               pb=&STRIP%obs_pb;
                                               pc=&STRIP%obs_pc;
                                               pd=&STRIP%obs_pd;
                                               phase=&STRIP%obs_phase;
                                               fom=&STRIP%obs_fom;
                                               sel=(ref_active=1);
                                               sel_test=(tst_active=1);
                                               statistics=true;)
     predict
        mode=dtarget(fcalc)
        to=dtarg
        selection=( &low_res >= d >= &high_res )
        atomselection=( &atom_select ) 
     end
   end if
   @@CNS_XTALMODULE:calculate_r_twin (fobs=&STRIP%obs_f;
                                      fcalc=fcalc;
                                      fpart=fbulk;
                                      twin_oper=&STRIP%twin_oper;
                                      twin_frac=&twin_frac;
                                      sel=(ref_active=1);
                                      sel_test=(tst_active=1);
                                      print=true;
                                      output=OUTPUT;
                                      r=$start_r;
                                      test_r=$start_test_r;)
 end

 {- detwin the observed data -}
 xray
   if ( &twin_frac < 0.5 ) then
     do (&STRIP%obs_f=fobs_orig) (all)
     do (&STRIP%obs_sigf=sigma_orig) (all)
   end if
   query name=fo_detwin domain=reciprocal end
   if ( $object_exist = false ) then
     declare name=fo_detwin domain=reciprocal type=real end
   else
     do (fo_detwin=0) (all)
   end if
   declare name=ref_det_active domain=reciprocal type=integer end
   declare name=tst_det_active domain=reciprocal type=integer end
 end

 xray

   @@CNS_XTALMODULE:data_detwin (fobs=&STRIP%obs_f;
                                 fcalc=fcalc;
                                 fpart=fbulk;
                                 twin_oper=&STRIP%twin_oper;
                                 twin_frac=&twin_frac;
                                 sel=(ref_active=1);
                                 sel_test=(tst_active=1);
                                 fobs_detwin=fo_detwin;
                                 ref_detwin=ref_det_active;
                                 tst_detwin=tst_det_active;
                                 ref_reject=$ref_detwin_reject;
                                 tst_reject=$tst_detwin_reject;
                                 total_reject=$detwin_reject;)
 end

 {- BEGIN MODIFICATION 8/01/06 -}
 @CNS_XTALMODULE:scale_and_solvent_grid_search (
                             bscale=&bscale;
                             sel=( ref_det_active=1 );
                             sel_test=( tst_det_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=fo_detwin;
                             obs_sigf=&STRIP%obs_sigf;
                             obs_i=$STRIP%obs_i;
                             obs_sigi=$STRIP%obs_sigi;                             
                             fpart=fbulk;

                             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;
                             
                             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;
                             
                             )
 
{- END MODIFICATION -}

 xray
   @@CNS_XTALMODULE:calculate_r (fobs=fo_detwin;
                                 fcalc=fcalc;
                                 fpart=fbulk;
                                 sel=(ref_det_active=1);
                                 sel_test=(tst_det_active=1);
                                 print=true;
                                 output=OUTPUT;
                                 r=$detwin_r;
                                 test_r=$detwin_test_r;)

 end

 xray
   declare name=map_phase domain=reciprocal type=real end 
   declare name=map_fom   domain=reciprocal type=real end 
   declare name=map_scale domain=reciprocal type=real end 
 end

 if ( &map_type = "unweighted" ) then

   xray

     do (map_phase=phase(fcalc+fbulk)) (all)

     do (total=fcalc+fbulk) (all)
     multiscale
       bfmin=-40 bfmax=40
       set1=fo_detwin  k1=-1 b1=0
       set2=total            b2=0
       selection=(ref_det_active=1)
     end
     do (map_scale=$k2) (all)

     do (map_fom=1.0) (all)

   end

 elseif ( &map_type = "sigmaa" ) then

   xray

     do (map_phase=phase(fcalc+fbulk)) (all)

     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

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

     do (map_fom=mod_fom) (all)
     do (map_scale=distribute(mod_dd)) (&low_res >= d >= &high_res)

     undeclare name=m          domain=reciprocal end
     undeclare name=mod_fom    domain=reciprocal end
     undeclare name=mod_x      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_dd     domain=reciprocal end

   end

 end if

 if ( &map_type = "gradient" ) then
 
   xray
     declare name=map domain=real end
     {- take the negative of the gradient so the map is the same sign
        as a standard difference map -}
     do (map=ft(-dtarg)) (ref_active=1)
   end

 else

   xray
     if ( &u = &v ) then
       do (fmap= 2 ((&u   map_fom combine(amplitude(fo_detwin),map_phase)) -
                    (&v map_scale        (fcalc+fbulk))))
          (ref_det_active=1 and acentric)
       do (fmap=    (&u   map_fom combine(amplitude(fo_detwin),map_phase)) -
                    (&v map_scale        (fcalc+fbulk)))
          (ref_det_active=1 and centric)
     else
       do (fmap=(&u             map_fom combine(amplitude(fo_detwin),map_phase)) -
                (&v           map_scale        (fcalc+fbulk)))
          (ref_det_active=1 and acentric)
       do (fmap=(max((&u-1),0)   map_fom combine(amplitude(fo_detwin),map_phase)) -
                (max((&v-1),0) map_scale        (fcalc+fbulk)))
          (ref_det_active=1 and centric)
       if ( &fill_in = true ) then
         do (fmap=(&u-&v) map_scale (fcalc+fbulk))
            (&low_res >= d >= &high_res and ref_det_active # 1)
       end if
     end if
   end

   xray
     declare name=map domain=real end
     if ( &fill_in = true ) then
       do (map=ft(fmap)) (&low_res >= d >= &high_res)
     else
       do (map=ft(fmap)) (ref_active=1)
     end if
   end

 end if

 xray
   undeclare name=fo_detwin domain=reciprocal end
   undeclare name=map_phase domain=reciprocal end 
   undeclare name=map_fom   domain=reciprocal end 
   undeclare name=map_scale domain=reciprocal end 
 end

 xray
   undeclare name=dtarg domain=reciprocal end
   undeclare name=total domain=reciprocal end
   undeclare name=fmap domain=reciprocal end
 end

 if ( $strict_ncs = true ) then
   ncs strict initialize end
   xray
     declare name=mask_ncs domain=real end
     mask
       mode=vdw
       solrad=&max
       shrink=0.
       to=mask_ncs
       sele=( &atom_select )
     end
   end
   inline @&ncs_infile
 end if

 xray
   if (&map_scale=true) then
      show rms (real(map)) ( all )
      if ( $result > 0 ) then
        do (map=map/$result) ( all )
      end if
   end if

   declare name=mask1 domain=real end
   mask mode=vdw solrad=0.5 shrink=0.7 to=mask1 
        sele=( &atom_select ) end

   declare name=mask2 domain=real end
   mask mode=vdw  solrad=&max shrink=0. to=mask2 sele=( &atom_select ) end

   peakpik
     from=map
     mpeak=&npeaks
     if ( $strict_ncs = true ) then
       selection=( mask1=1 and mask2=0 and mask_ncs=0 )
     else
       selection=( mask1=1 and mask2=0 )
     end if
     proximity=( &atom_select )
     atom=true
   end

 end
 
 if ( $strict_ncs = true ) then
   xray
     undeclare name=mask_ncs domain=real end
   end
 end if

 do (mass=16) ( segid PEAK ) 

 delete sele=( segid PEAK and ( attr B < &peak ) or ( attr B = &peak ) ) end

 delete sele=( segid PEAK and not ( 
             ( &atom_select and not segid PEAK) saround &max ) )
 end

 if ( &peak_select = "distance" ) then
   delete sele= ( segid PEAK and ( 
                ( &atom_select and not segid PEAK) saround &min ) )
   end
 elseif ( &peak_select = "hbond" ) then
   delete sele= ( segid PEAK and (
                ( &atom_select and not segid PEAK ) saround &hmin ) )
   end
   delete sele=      ( segid PEAK and 
                 ( ( ( &atom_select and not segid PEAK ) saround &min ) and
              not  ( ( ( &atom_select and not segid PEAK ) and 
                       ( name O* or name N* ) ) saround &min ) ) )
   end
 end if

 if ( &peak_hbond = true ) then
   delete sele=( segid PEAK and not ( 
               ( &atom_select and not segid PEAK) saround &hmax ) )
   end
   delete sele=      ( segid PEAK and 
                 ( ( ( &atom_select and not segid PEAK ) saround &hmax ) and
              not  ( ( ( &atom_select and not segid PEAK ) and 
                       ( name O* or name N* ) ) saround &hmax ) ) )
   end
 end if

 parameter nbonds special_position=&special_dist end end

 if ( &peak_special = false ) then
   xray
     @@CNS_XRAYLIB:scatter.lib
     special selection=( segid PEAK ) to=rmsd end
   end
   delete sele=( segid PEAK and ( attribute rmsd > 1 ) ) end
 end if

 show sum(1) ( segid PEAK )
 if ( $result = 0 ) then
   display Aborting: No waters picked or all removed by deletion criteria
   stop
 end if

 flags include pvdw end

 parameter
   nbond wmin=&max end
   nonbonded PEAK 0.1591   2.8509   0.1591   2.8509
 end

 show min(decode(resid)) ( segid PEAK )
 evaluate ($minid=$result)
 show max(decode(resid)) ( segid PEAK )
 evaluate ($maxid=$result)

 evaluate ($counter=$maxid) 

 while ( $counter >= $minid ) loop id

   show sum(1) (segid PEAK and resid $counter)

   if ( $select = 1 ) then
  
     igroup interaction
       ( segid PEAK and resid $counter ) ( segid PEAK )
     end

     distance from=( segid PEAK and resid $counter ) 
                to=( segid PEAK ) 
       disposition=rmsd cuton=0.0 cutoff = &max
     end
 
     show min(rmsd) ( segid PEAK and ( attribute rmsd > 0 ) )

     if ( $select > 0 ) then
       if ( $result < &hmin ) then
         delete 
           sele=( segid PEAK and resid $counter )
         end
       end if
     end if

   end if

   evaluate ($counter=$counter-1)

 end loop id

 {- BEGIN MODIFICATION 8/01/06 -}
 show sum(1) ( segid PEAK )
 if ( $result = 0 ) then
   display ==============================================================
   display Warning: No waters picked or all removed by deletion criteria
   display          No output files written
   display ==============================================================
   stop
 end if
 {- END MODIFICATION -}

 show sum(1) ( segid PEAK )
 if ( $result = 0 ) then
   display Aborting: All waters removed by deletion criteria
   stop
 end if

 flags exclude pvdw end

 do (resname=&wat_resname) ( segid PEAK)
 if ( &wat_resname = "HOH" ) then
   evaluate ($wat_oname="O")
 else
   evaluate ($wat_oname="OH2")
 end if
 do (name=$wat_oname) ( segid PEAK)
 do (chemical="OT")   ( segid PEAK)

 evaluate ($wat_segid=capitalize(&wat_segid))

 show max(decode(resid)) (segid $wat_segid)
 if ( $result > 0 ) then
   evaluate ($maxsolv=$result+1)
 else
   evaluate ($maxsolv=&wat_number)
 end if

 show ave(b) ( &atom_select and not ( segid PEAK ) )

 do (b=$result) ( segid PEAK )

 evaluate ($read_wat=true)
 for $counter in ( 1 2 3 4 5 ) loop wat
   if ( &parameter_infile_$counter = "CNS_TOPPAR:water_rep.param" ) then
     evaluate ($read_wat=false)
   end if
 end loop wat

 if ( $read_wat = true ) then
   parameter
     @@CNS_TOPPAR:water_rep.param
   end
 end if

 xray
   @CNS_XRAYLIB:scatter.lib
 end

 if ( &refine = true ) then

   xray
     associate fcalc ( &atom_select or (segid PEAK) )
   end

   xray
     predict
       mode=reciprocal
       to=fcalc
       selection=(ref_active=1)
       atomselection=( &atom_select or (segid PEAK) )
     end
   end

   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

   {- BEGIN MODIFICATION 8/01/06 -}
   @CNS_XTALMODULE:scale_and_solvent_grid_search (
                             bscale=&bscale;
                             sel=( ref_active=1 );
                             sel_test=( tst_active=1 );
                             atom_select=( &atom_select or (segid PEAK) );
                             bulk_sol=&bulk_sol;
                             bulk_mask=&bulk_mask_infile;
                             bulk_atoms=( &atom_select or (segid PEAK) );
                             
                             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;

                             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;
                             
                             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;
                             
                             )
  {- END MODIFICATION -}

   xray
     @@CNS_XTALMODULE:refinementtarget_twin (target=&reftarget;
                                             sig_sigacv=0.07;
                                             mbins=&target_bins;
                                             fobs=&STRIP%obs_f;
                                             sigma=&STRIP%obs_sigf;
                                             weight=$STRIP%obs_w;
                                             iobs=$STRIP%obs_i;
                                             sigi=$STRIP%obs_sigi;
                                             test=tst_active;
                                             twin_oper=&STRIP%twin_oper;
                                             twin_frac=&twin_frac;
                                             fcalc=fcalc;
                                             fpart=fbulk;
                                             pa=&STRIP%obs_pa;
                                             pb=&STRIP%obs_pb;
                                             pc=&STRIP%obs_pc;
                                             pd=&STRIP%obs_pd;
                                             phase=&STRIP%obs_phase;
                                             fom=&STRIP%obs_fom;
                                             sel=(ref_active=1);
                                             sel_test=(tst_active=1);
                                             statistics=true;)
   end

   fix selection=( not ( segid PEAK ) ) end

   igroup interaction
     ( &atom_select or (segid PEAK) ) 
     ( segid PEAK )
   end

   flags exclude elec include xref pvdw end

   @@CNS_XTALMODULE:getweight (selected=(&atom_select or (segid PEAK));
                               fixed=(not(segid PEAK));)

   flags exclude * include xref vdw pvdw end

   minimize lbfgs
     nstep=&minimize_nstep
     drop=10
   end

   flags exclude * include xref end   

   @@CNS_XTALMODULE:one_term_wa (wa=$wa;)

   xray
     optimize bfactors

       bmin=&bmin
       bmax=&bmax

       nstep=&bfactor_nstep
       drop=10.0

       rweight=0
     end
   end

   do ( q = 0 ) ( segid PEAK and ( attribute b > &bfactor_max ) )

   fix selection=( not ( segid PEAK ) or (segid PEAK and (attr q = 0) ) ) end

   igroup interaction
     ( &atom_select or (segid PEAK and (attr q > 0) ) ) 
     ( segid PEAK and (attr q > 0) )
   end

   if ( $select > 0 ) then

     @@CNS_XTALMODULE:one_term_wa (wa=$wa;)

     xray
       optimize bfactors

         bmin=&bmin
         bmax=&bmax

         nstep=&bfactor_nstep
         drop=1.0

         rweight=0
       end
     end
   end if

 end if

 xray
   predict
     mode=reciprocal
     to=fcalc
     selection=( &low_res >= d >= &high_res )
     atomselection=( &atom_select or (segid PEAK)) 
   end
   @@CNS_XTALMODULE:calculate_r_twin (fobs=&STRIP%obs_f;
                                      fcalc=fcalc;
                                      fpart=fbulk;
                                      twin_oper=&STRIP%twin_oper;
                                      twin_frac=&twin_frac;
                                      sel=(ref_active=1);
                                      sel_test=(tst_active=1);
                                      print=true;
                                      output=OUTPUT;
                                      r=$full_r;
                                      test_r=$full_test_r;)
 end

 if ( &refine = true ) then
   delete sele=( segid PEAK and ( attribute q = 0 ) ) end
 end if

 show sum(1) ( segid PEAK )
 if ( $result = 0 ) then
   display Aborting: All waters removed by deletion criteria
   stop
 end if

 if ( &merge_coord = false ) then
   delete sele=( not ( segid PEAK and name $wat_oname ) ) end
 end if

 show sum(1) (name $wat_oname and segid PEAK)
 evaluate ($new_water=$result)

 evaluate ($counter=$maxsolv)
 for $id in id (segid PEAK and name $wat_oname) loop water
   do (resid=encode($counter)) ((id $id) and segid PEAK)
   do (segid=$wat_segid) ((id $id) and segid PEAK)
   evaluate ($counter=$counter+1)   
 end loop water

 set display=&coordinate_outfile end

 xray
   show sum (1) (tst_det_active=1)
   if ( $result > 0 ) then
     evaluate ($test_exist=true)
   else
     evaluate ($test_exist=false)
   end if
 end

 evaluate ($remark="in ")

 if ( &map_type = "unweighted" ) then
   evaluate ($remark=$remark + "(" + encode(&u) +  " |Fo| - " +
                                     encode(&v) + " k|Fc|)e^(i phi_calc)")
 elseif ( &map_type = "sigmaa" ) then
   evaluate ($remark=$remark + "(" + encode(&u) + " m|Fo| - " +
                                     encode(&v) + " D|Fc|)e^(i phi_calc)")
   if ( $test_exist = true ) then
     evaluate ($remark=$remark + " cross-val.")
   end if
   evaluate ($remark=$remark + " sigmaa")
 elseif ( &map_type = "combined" ) then
   evaluate ($remark=$remark + "(" + encode(&u) + " m|Fo|)e^(i phi_comb) - " +
                               "(" + encode(&v) + " D|Fc|)e^(i phi_calc)")
   if ( $test_exist = true ) then
     evaluate ($remark=$remark + " cross-val.")
   end if
   evaluate ($remark=$remark + " sigmaa")
 elseif ( &map_type = "observed" ) then
   evaluate ($remark=$remark + "(" + encode(&u) +  " m|Fo|)e^(i phi_obs) - " +
                                 "(" + encode(&v) + " k m|Fc|)e^(i phi_calc))")
 elseif ( &map_type = "gradient" ) then
   evaluate ($remark=$remark + "( d(" + &reftarget + ")/dFc )")
 end if

 evaluate ($remark=$remark + " map")

 display REMARK coordinates from water picking from twinned data
 display REMARK twinning operator= &STRIP%twin_oper  twinning fraction= &twin_frac
 display REMARK $new_water waters picked at level greater than &peak
 display REMARK $remark
 display REMARK peak selection criteria: &STRIP%peak_select
 if ( &peak_select = "distance" ) then
   display REMARK peaks closer than &min A or further than &max A were deleted
 elseif ( &peak_select = "hbond" ) then
   display REMARK peaks closer than &min A or further than &max A were deleted
   display REMARK but peaks &hmin A from oxygen or nitrogen were kept
 else
   display REMARK peaks further than &max A were deleted
 end if
 if ( &peak_hbond = true ) then
   display REMARK peaks further than &hmax A from oxygen or nitrogen were deleted
 end if
 display REMARK map resolution: &low_res - &high_res A  
 if ( $total_test > 0 ) then
   display REMARK starting twinned r= $start_r[f6.4] twinned free_r= $start_test_r[f6.4]
   display REMARK final    twinned r= $full_r[f6.4] twinned free_r= $full_test_r[f6.4]
 else
   display REMARK starting twinned r= $start_r[f6.4]
   display REMARK final    twinned r= $full_r[f6.4]
 end if
 display REMARK sg= &STRIP%sg a= &a b= &b c= &c alpha= &alpha beta= &beta gamma= &gamma
 if ( &BLANK%parameter_infile_1 = false ) then
    display REMARK parameter file 1  : &STRIP%parameter_infile_1
 end if
 if ( &BLANK%parameter_infile_2 = false ) then
    display REMARK parameter file 2  : &STRIP%parameter_infile_2
 end if
 if ( &BLANK%parameter_infile_3 = false ) then
    display REMARK parameter file 3  : &STRIP%parameter_infile_3
 end if
 if ( &BLANK%parameter_infile_4 = false ) then
    display REMARK parameter file 4  : &STRIP%parameter_infile_4
 end if
 if ( &BLANK%parameter_infile_5 = false ) then
    display REMARK parameter file 5  : &STRIP%parameter_infile_5
 end if
 display REMARK molecular structure file: &STRIP%structure_infile
 display REMARK input coordinates: &STRIP%coordinate_infile
 if ( &BLANK%anom_library = false ) then
   display REMARK anomalous f' f'' library: &STRIP%anom_library
 end if
 if ( &BLANK%reflection_infile_1 = false ) then
   display REMARK reflection file= &STRIP%reflection_infile_1
 end if
 if ( &BLANK%reflection_infile_2 = false ) then
   display REMARK reflection file= &STRIP%reflection_infile_2
 end if
 if ( &BLANK%reflection_infile_3 = false ) then
   display REMARK reflection file= &STRIP%reflection_infile_3
 end if
 if ( &BLANK%ncs_infile = false  ) then
   display REMARK ncs= &STRIP%ncs_type  ncs file= &STRIP%ncs_infile
 else
   display REMARK ncs= none
 end if
 if ( &bscale # "no" ) then
   if ( $low_b_flag = true ) then
     display REMARK warning: B-correction gave atomic B-values less than zero
     display REMARK          they have been reset to zero
   end if
 end if
 
 ! 
 ! Begin modification (8/01/06)
 if ( &bscale = "anisotropic" ) then
   display REMARK Anisotropic B-factor tensor Ucart of atomic model without isotropic component :
   display REMARK   B11=$Baniso_11[f8.3] B22=$Baniso_22[f8.3] B33=$Baniso_33[f8.3]
   display REMARK   B12=$Baniso_12[f8.3] B13=$Baniso_13[f8.3] B23=$Baniso_23[f8.3]
   display REMARK Isotropic component added to coordinate array B: $Biso_model[f8.3]
 elseif ( &bscale = "isotropic" ) then
   display REMARK B-factor applied to coordinate array B: $Biso_model[f8.3]
 else
   display REMARK initial B-factor correction: none
 end if

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

 if ( &obs_type = "intensity" ) then
   display REMARK reflections with Iobs/sigma_I < &sigma_cut rejected
   display REMARK reflections with Iobs > &obs_rms * rms(Iobs) rejected
   display REMARK reflections with Iobs[&STRIP%twin_oper] = 0 rejected
 else
   display REMARK reflections with |Fobs|/sigma_F < &sigma_cut rejected
   display REMARK reflections with |Fobs| > &obs_rms * rms(Fobs) rejected
   display REMARK reflections with |Fobs|[&STRIP%twin_oper] = 0 rejected
 end if
 xray anomalous=? end
 if ( $result = true ) then
   display REMARK anomalous diffraction data was input
 end if
 
 {- MODIFIED 8/01/06 -}
 display REMARK fft gridding factor = $fft_grid, B factor offset = $fft_b_add A^2, Elimit = $fft_elim
 {- END MODIFICATION -}

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

 remark

 write structure output=&structure_outfile end

 @CNS_XTALMODULE:write_pdb (pdb_o_format=&pdb_o_format;
                            coordinate_outfile=&coordinate_outfile;
                            sgparam=$sgparam;)

 stop