{+ file: sa_omit_map.inp +}
{+ directory: xtal_refine +}
{+ description: Make an annealed omit map with spherical omitted region. +}
{+ comment:
            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))
              (u m|Fo| e^(i phi_comb)) - (v D|Fc| e^(i phi_calc))
              (u m_obs|Fo| e^(i phi_obs )) - (v k m_obs|Fc| e^(i phi_calc))
               d(target)/dFc
              where Fc is the calculated structure factor
              where m and D are derived from sigmaa. +}
{+ authors: Axel T. Brunger, Alec Hodel and Paul D. Adams +}
{+ copyright: Yale University +}

{+ reference: A. Hodel, S.-H. Kim, and A.T. Brunger, Model Bias in 
              Macromolecular Crystal Structures, Acta Cryst. A48, 
              851-859 (1992) +}
{+ reference: L.M. Rice and A.T. Brunger, Torsion Angle Dynamics:
              Reduced Variable Conformational Sampling Enhances
              Crystallographic Structure Refinement, Proteins: Structure,
              Function, and Genetics, 19, 277-290 (1994) +}
{+ reference: R.J. Read, Improved Fourier coefficients for maps using 
              phases from partial structures with errors. Acta Cryst. 
              A42, 140-149 (1986) +}
{+ 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 
   - the selections store1 through store4 are available for general use -}

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

{* force field settings file *}
{===>} force_field_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: required *}
{* cross-validation should always be used, with the possible exception
   of a final round of refinement including all data *}
{* cross-validation is always required for the maximum likelihood targets *}
{===>} 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.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_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 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 all atoms present in model *}
{* this should include all conformations if multiple conformations are used *}
{===>} atom_select=(known and not hydrogen);

{* select atoms to be omitted from omit map refinement and calculation *}
{===>} atom_omit=(residue 20 or residue 21);

{* size of the neighborhood that is also omitted *}
{===>} sphere=3.5;

{* size of the surrounding harmonically restrained cushion *}
{===>} cushion=2;

{* select fixed atoms *}
{* note: isolated atoms and diatomic molecules are automatically 
   fixed during torsion angle dynamics, and atoms at special positions are 
   automatically fixed for all types of dynamics. So, you don't have to 
   explicitly fix them here. *}
{===>} atom_fixed=(none);

{- begin modification, ATB, 11/29/08 -}
{* atom selections for non-default rigid groups during torsion angle dynamics *}
{* note: the selections must be non-overlapping *}
{===>} atom_rigid_1=(none);
{===>} atom_rigid_2=(none);
{===>} atom_rigid_3=(none);
{===>} atom_rigid_4=(none);
{===>} atom_rigid_5=(none);
{===>} atom_rigid_6=(none);
{===>} atom_rigid_7=(none);
{===>} atom_rigid_8=(none);
{===>} atom_rigid_9=(none);
{===>} atom_rigid_10=(none);
! to add more groups add more numbered entries:
!   {===>} atom_rigid_11=(none);
!   {===>} atom_rigid_12=(none);
!   {===>} atom_rigid_13=(none);
! etc

{* select atoms to be harmonically restrained during annealing *}
{===>} atom_harm=(none);

{* harmonic restraint constant - for harmonically restrained atoms *}
{===>} k_harmonic=10;

{* select atoms in alternate conformation 1 *}
{===>} conf_1=(none);

{* select atoms in alternate conformation 2 *}
{===>} conf_2=(none);

{* select atoms in alternate conformation 3 *}
{===>} conf_3=(none);

{* select atoms in alternate conformation 4 *}
{===>} conf_4=(none);

{* additional restraints file *}
{* eg. auxiliary/dna-rna_restraints.def *}
{===>} restraints_infile="";

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

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

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

{* 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
         combined: (u m|Fo|^exp(i phi_comb) - v D|Fc|^exp(i phi_calc))
                    model and experimental phases combined, m and D from sigmaa
         observed: (u m|Fo|^exp(i phi_obs) - v k m|Fc|^exp(i phi_calc))
                    observed phases and fom from phase probability distribution
         gradient: d(target)/dFc
                   gradient of the current crystallographic target wrt Fc
       NB. experimental phases must be supplied as a phase
           probability distribution in the Hendrickson-Lattman arrays +}
{+ choice: "sigmaa" "unweighted" "combined" "observed" "gradient" +}
{===>} map_type="sigmaa";

{* 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 format *}
{+ choice: "cns" "ezd" +}
{===>} map_format="cns";

{* b-factor sharpening (A^2), for example, -100 *}
{===>} bsharp=0;

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

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

{* select atoms around which map will be written *}
{* change if different to atoms selected for map calculation *}
{===>} atom_map=(known and not hydrogen);

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

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

{======================= annealing parameters ========================}

{* carry out annealing *}
{+ choice: true false +}
{===>} anneal=true;

{* type of molecular dynamics *}
{* it is possible that torsion angle dynamics may encounter a problem
   with the molecular topology after certain parts of the structure have
   been omitted. If this happens the program will terminate with an
   error and cartesian dynamics must be used instead *}
{+ choice: "torsion" "cartesian" +}
{===>} md_type="cartesian";

{* starting temperature *}
{===>} temperature=1000;

{* drop in temperature (K) per set of dynamics *}
{===>} cool_rate=25;

{* seed for random number generator *}
{* change to get different initial velocities *}
{===>} seed=82364;

{* torsion angle topology modification file *}
{===>} torsion_infile="CNS_TOPPAR:torsionmdmods";

{======================= minimization parameters ======================}

{* final steps of conjugate gradient minimization *}
{===>} mini_steps=100;

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

{* refinement target *}
{+ list: mlf: maximum likelihood target using amplitudes
         mli: maximum likelihood target using intensities
        mlhl: maximum likelihood target using amplitudes
              and phase probability distribution
    residual: standard crystallographic residual
      vector: vector residual
       mixed: (1-fom)*residual + fom*vector
        e2e2: correlation coefficient using normalized E^2
        e1e1: correlation coefficient using normalized E
        f2f2: correlation coefficient using F^2
        f1f1: correlation coefficient using F +}
{+ choice: "mlf" "mli" "mlhl" "residual" "vector" "mixed"
           "e2e2" "e1e1" "f2f2" "f1f1" +}
{===>} reftarget="mlf";

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

{* Wa weight for X-ray term *}
{* this will be determined automatically if a negative value is given.
   Note: wa can be very different depending on the target - if it is not
         determined automatically make sure an appropriate value is used *}
{===>} wa=-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;

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

{* root name for output files *}
{+ list:
         map file will be in: <output_root>.map
         positive peaks in: <output_root>_positive.peaks
         negative peaks in: <output_root>_negative.peaks
         Fourier coefficients will be in: <output_root>.coeff +}
{===>} output_root="sa_omit_map";

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

{* do peak picking on map *}
{* optional - use water_pick.inp to pick waters *}
{+ choice: true false +}
{===>} peak_search=true;

{* number of peaks to pick from map *}
{===>} peak_num=30;

{* write a reflection file with the Fourier coefficients of the map *}
{+ list: arrays written:
         map_fom:   FOM weight applied to observed data
         map_phase: phases used for observed data
         map_scale: scale factor applied to calculated data
         map_coeff: Fourier map coefficients - map=ft(map_coeff) +}
{+ choice: true false +}
{===>} write_coeff=true;

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

 set message=normal echo=on end

 xray

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

   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
   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)
   else
     evaluate ($reject_obs=&obs_f)
     evaluate ($reject_sig=&obs_sigf)
   end if

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

   do (all_active=0) ( all )
   do (all_active=1) ( &low_res >= d >= &high_res )

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

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

   show sum (1) ( map_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) (map_active=1 and &STRIP%test_set=&test_flag)
   end if

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

   do (map_active=0) ( all )
   do (map_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 ) )

   show sum (1) ( map_active=1 and tst_active=0 )
   evaluate ($total_work=$select)
   show sum (1) ( map_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=( map_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=&grid;
                             fft_memory=&fft_memory;
                             fft_grid=$fft_grid;   
                             fft_b_add=$fft_b_add; 
                             fft_elim=$fft_elim; 
                                      )
                            
 xray
 {- END MODIFICATION -}

   if ( &wa >= 0 ) then
      wa=&wa
   end if
      
 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

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

 if ( &BLANK%restraints_infile = false ) then
     @&restraints_infile
 end if

 do (store7=0) (all)

 evaluate ($nalt=1)
 evaluate ($alt=1)
 evaluate ($done=false)
 while ( $done = false ) loop nalt
   if ( &exist_conf_$alt = true ) then
     show sum(1) ( &conf_$alt )
     if ( $result > 0 ) then
       evaluate ($nalt=$nalt+1)
     end if
   else
     evaluate ($done=true)
     evaluate ($nalt=$nalt-1)
   end if
   evaluate ($alt=$alt+1)
 end loop nalt

 evaluate ($alt=1)
 while ( $alt <= $nalt ) loop alt
   do (store7=$alt) ( &conf_$alt )
   evaluate ($alt=$alt+1)
 end loop alt

 ident ( store8 ) ( byresidue ( &atom_omit saround &sphere ) )
 
 {- BEGIN MODIFICATION ATB 5/14/08 -}
 show sum ( 1 ) ( store8 )
  if ($select=0) then
  display 
  display  *****************************************************************
  display  Error: selection for omitted atoms is zero.  Check atom selection 
  display         for "select atoms to be omitted from map calculation"
  display  *****************************************************************
  display
  abort
 end if

 evaluate ($atom_omit_select=$select)
 {- END MODIFICATION -}
 
 ident ( store9 )  ( byresidue ( ( store8 ) saround &cushion and
                               not store8 ) )

 show sum ( 1 ) ( store9 )
 evaluate ($cushion_select=$select)
			       
 igroup
   interaction ( &atom_select and not(attr store7 > 0))
               ( &atom_select and not(attr store7 > 0))
   evaluate ($alt=1)
   while ( $alt <= $nalt ) loop alcs
     interaction ( &atom_select and
                   ( attr store7 = $alt or attr store7 = 0 ))
                 ( &atom_select and ( attr store7 = $alt ))
     evaluate ($alt=$alt+1)
   end loop alcs
 end
 
 {- check isolated atoms and atoms at special positions and add to
    list of fixed atoms if needed - store7 will be used -}
 
 if (&anneal=true) then
    evaluate ($mode=&md_type)
 else
    evaluate ($mode="minimization")
 end if

 @CNS_XTALMODULE:setupfixed (
                           mode=$mode;
                           atom_select=&atom_select;
                           atom_fixed=&atom_fixed;
                           atom_total_fixed=store7;
                           atom_multiplicity=rmsd;
                           mset=$mset;
                           )
                           
 fix selection=( store7 ) end

 xray
   associate fcalc ( &atom_select and not store8 )
 end

 fastnb grid end

 flags                                       
    include xref                   
   ?                                        
 end      
 if ( &BLANK%force_field_infile = true ) then
    flags                                       
       exclude elec pele include vdw pvdw                 
      ?                                        
    end      
 else
     @&force_field_infile
 end if                                   

 set seed=&seed end

 if ( &md_type = "torsion" ) then
   evaluate ($start_temp=&temperature)
   evaluate ($time_step=0.004)
   evaluate ($md_steps=6)
   evaluate ($fbeta=200)
 end if

 if ( &md_type = "cartesian" ) then
   evaluate ($start_temp=&temperature)
   evaluate ($time_step=0.0005)
   evaluate ($md_steps=50)
   evaluate ($fbeta=100)
 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
 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

 xray
    predict
      mode=reciprocal
      to=fcalc
      selection=( all_active=1 )
      atomselection=( &atom_select and not store8 )
    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=&grid;
                             fft_memory=&fft_memory;
                             fft_grid=$fft_grid;   
                             fft_b_add=$fft_b_add; 
                             fft_elim=$fft_elim; 
                                      )
 {- END MODIFICATION -}
 
 xray
  
   tselection=(ref_active=1)

   @@CNS_XTALMODULE:refinementtarget (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;
                                     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

 if ( &anneal = true ) then

   xray
     tolerance=0.2 lookup=true
   end

   do ( harm=0 ) ( all )
   do ( harmonic=&k_harmonic ) ( &atom_harm )
   do ( harmonic=20 ) ( store9 or store8 )
   do ( refx=x ) ( all )
   do ( refy=y ) ( all )
   do ( refz=z ) ( all )
   flags include harm end
  
   if ( &wa < 0 ) then
     @@CNS_XTALMODULE:getweight (
                                selected=&atom_select;
                                fixed=(store7);
                               )
   end if
 
   parameter
     nbonds
        repel ?  evaluate ($repel_old=$result)
        rcon ?   evaluate ($rcon_old=$result)
        if ($repel_old =1  ) then
           repel=1.  rcon=100.
        else
           repel=.75 rcon=50.
        end if
     end
   end

   {* modification: moved mass modification before velocity assignment, ATB 7/8/2008 *}
   do (store6=mass) ( all )
   do (mass=max(10,min(30,mass))) ( all )

   do (fbeta=$fbeta) ( ( &atom_select ) and not (store7) )

   do (vx=maxwell($start_temp)) (&atom_select and not (store7))
   do (vy=maxwell($start_temp)) (&atom_select and not (store7))
   do (vz=maxwell($start_temp)) (&atom_select and not (store7))

 if ( &md_type = "torsion" ) then
   dynamics torsion
     topology
     
       maxlength=-1
       maxchain=-1
       maxtree=-1
       kdihmax = 95.
       
         {- begin modification, ATB, 11/29/08 -}
         @CNS_XTALMODULE:messagesilent(&message_old_silent=$message_old_silent; &echo_old_silent=$echo_old_silent; )     
         evaluate ($atr_count=1)
         evaluate ($atr_done=false)
         while ( $atr_done = false ) loop atrl
          if ( &exist_atom_rigid_$atr_count = true ) then
           fix group ( &atom_rigid_$atr_count ) 
           evaluate ($atr_count=$atr_count+1)
          else
           evaluate ($atr_done=true)
          end if
         end loop atrl
         @CNS_XTALMODULE:messagerevert(&message_old_silent=$message_old_silent; &echo_old_silent=$echo_old_silent; )     
         
         if ( &BLANK%torsion_infile = false ) then
            @&torsion_infile
         else         
            @CNS_TOPPAR:torsionmdmods
         end if
       
     end
     nstep=0
     cmremove=true
   end
 end if


 evaluate ( $curr_temp = &temperature )
 while ( $curr_temp > 0.0 ) loop cool
   if ( &md_type = "torsion" ) then
     dynamics torsion
       timestep=$time_step
       nstep=$md_steps
       nprint=5
       cmremove=false
       vscaling=true
       temperature=$curr_temp
     end
   end if
   if ( &md_type = "cartesian" ) then
     dynamics cartesian
       if ($curr_temp=&temperature) then
          cmremove=true
       else
          cmremove=false
       end if
       timestep=$time_step
       nstep=$md_steps
       nprint=10
       vscaling=true
       temperature=$curr_temp
     end
   end if
   evaluate ( $curr_temp = $curr_temp - &cool_rate )
 end loop cool

 do (mass=store6) ( all )

 if ( &md_type = "torsion" ) then
   dynamics torsion
     nstep = 0
     cmremove=false
     topology
       reset
     end
   end
 end if

 parameter
     nbonds
       repel=$repel_old rcon=$rcon_old
     end
   end
 end if

 xray
   tolerance=0.0 lookup=false
 end
 
 if ( &wa < 0 ) then
    @@CNS_XTALMODULE:getweight (
                               selected=&atom_select;
                               fixed=(store7);
                              )
 end if
   
 if ( &mini_steps > 0 ) then
   minimize lbfgs
     nstep=&mini_steps
     drop=10.0 nprint=10
   end
 end if
 
 xray

   tselection=(map_active=1)

   if ( &map_type = "gradient" ) then
     predict
       mode=dtarget(fcalc)
       to=dtarg
       selection=( map_active=1 )
       atomselection=( &atom_select and not store8 )
     end
   else
     predict
        mode=reciprocal
        to=fcalc
        selection=( all_active=1 )
        atomselection=( &atom_select and not store8 )
     end
   end if
 end

 xray
  @@CNS_XTALMODULE:calculate_r (fobs=&STRIP%obs_f;
                                fcalc=fcalc;
                                fpart=fbulk;
                                sel=(map_active=1);
                                sel_test=(tst_active=1);
                                print=true;
                                output=OUTPUT;
                                r=$map_r;
                                test_r=$map_free_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=&STRIP%obs_f  k1=-1 b1=0
       set2=total          b2=0
       selection=(map_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=&STRIP%obs_f;
                                  fcalc=fcalc;
                                  fpart=fbulk;
                                  test=tst_active;
                                  sel=(map_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; )

     do (map_fom=mod_fom) (all)
     do (map_scale=distribute(mod_dd)) (all_active=1)

     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

 elseif ( &map_type = "combined" ) then

   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

     @CNS_XTALMODULE:fomsigmaacv ( sig_sigacv=0.07;
                                  mbins=&target_bins;
                                  statistics=true;
                                  fobs=&STRIP%obs_f;
                                  fcalc=fcalc;
                                  fpart=fbulk;
                                  test=tst_active;
                                  sel=(map_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; )

     @CNS_XTALMODULE:combineprobability ( messages="off";
                                          addname="model phases";
                                          pa=mod_pa;
                                          pb=mod_pb;
                                          pc=mod_pc;
                                          pd=mod_pd;
                                          w=1;
                                          addname="experimental phases";
                                          adda=&STRIP%obs_pa;
                                          addb=&STRIP%obs_pb;
                                          addc=&STRIP%obs_pc;
                                          addd=&STRIP%obs_pd;
                                          addw=1;)

     @CNS_XTALMODULE:getfom ( pa=mod_pa;
                              pb=mod_pb;
                              pc=mod_pc;
                              pd=mod_pd;
                              m=m;
                              phistep=5; )

     do (map_phase=phase(m)) (all)
     do (map_fom=amplitude(m)) (all)
     do (map_scale=distribute(mod_dd)) (all_active=1)

     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

 elseif ( &map_type = "observed" ) then

   xray

     do (total=fcalc+fbulk) (all)
     multiscale
       bfmin=-40 bfmax=40
       set1=&STRIP%obs_f  k1=-1 b1=0
       set2=total          b2=0
       selection=(map_active=1)
     end
     do (map_scale=$k2) (all)

     declare name=m          domain=reciprocal type=complex end

     @CNS_XTALMODULE:getfom ( pa=&STRIP%obs_pa;
                              pb=&STRIP%obs_pb;
                              pc=&STRIP%obs_pc;
                              pd=&STRIP%obs_pd;
                              m=m;
                              phistep=5; )

     do (map_phase=phase(m)) (all)
     do (map_fom=amplitude(m)) (all)
     do (map_scale=map_scale*map_fom) (all)
     undeclare name=m domain=reciprocal end
   end

 end if
 
 if ( &map_type = "gradient" ) then
 
   xray

     declare name=map domain=real end
     
     {- BEGIN MODIFICATION -}
     if (&bsharp # 0) then 
        {- b-factor sharpening -} 
        evaluate ($bsharp= (-1) * &bsharp) 
        {- take the negative of the gradient so the map is the same sign
        as a standard difference map -}
        do (fmap=exp( $bsharp * s^2/4) * (-1) * dtarg) ( all )
 
        do (map=ft(fmap)) (map_active=1)
     end if
     {- END MODIFICATION -}
     
   end

 else

   xray
     if ( &u = &v ) then
       do (fmap= 2 ((&u   map_fom combine(amplitude(&STRIP%obs_f),map_phase)) -
                    (&v map_scale        (fcalc+fbulk))))
          (map_active=1 and acentric)
       do (fmap=    (&u   map_fom combine(amplitude(&STRIP%obs_f),map_phase)) -
                    (&v map_scale        (fcalc+fbulk)))
          (map_active=1 and centric)
     else
       do (fmap=(&u             map_fom combine(amplitude(&STRIP%obs_f),map_phase)) -
                (&v           map_scale        (fcalc+fbulk)))
          (map_active=1 and acentric)
       do (fmap=(max((&u-1),0)   map_fom combine(amplitude(&STRIP%obs_f),map_phase)) -
                (max((&v-1),0) map_scale        (fcalc+fbulk)))
          (map_active=1 and centric)
       if ( &fill_in = true ) then
         do (fmap=(&u-&v) map_scale (fcalc+fbulk))
            ( all_active=1 and map_active # 1)
       end if
     end if
   end
   
   {- BEGIN MODIFICATION -}
   {- b-factor sharpening -} 
   evaluate ($bsharp= (-1) * &bsharp) 
   xray
      do (fmap=exp( $bsharp * s^2/4) * fmap) ( all )
   end
   {- END MODIFICATION -}

   xray
     declare name=map domain=real end
     if ( &u = &v ) then
       do (map=ft(fmap)) ( map_active=1 )
     else
       if ( &fill_in = true ) then
         do (map=ft(fmap)) ( all_active=1 )
       else
         do (map=ft(fmap)) ( map_active=1 )
       end if
     end if
   end
  
 end if

 if ( &write_coeff = true ) then
   evaluate ($coeff_out=&output_root + ".coeff")
   xray
     declare name=map_coeff domain=reciprocal type=complex end
     do (map_coeff=fmap) (all)

    set display=$coeff_out end
   end

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

   xray
     write reflection
       output=$coeff_out
       if ( &map_type = "gradient" ) then
         sele=(map_active=1)
       elseif ( &fill_in = true ) then
         sele=(&low_res >= d >= &high_res)
       else
         sele=(map_active=1)
       end if
       map_fom map_phase map_scale map_coeff
     end
     
     set display=OUTPUT end
     
     undeclare name=map_coeff domain=reciprocal end
   end
 end if

 xray
   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

 set remarks=reset end
 set remarks=accumulate end

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

 evaluate ($remark="")
 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 + " omit map")

 if ( $total_test > 0 ) then
   remark $remark  r= $map_r[f6.4] free_r= $map_free_r[f6.4]
 else
   remark $remark  r= $map_r[f6.4]
 end if
 
 {- BEGIN MODIFICATION -}
 remark number of omitted atoms: $atom_omit_select   number of atoms in cushion: $cushion_select
 {- END MODIFICATION -}

 if ( &obs_type = "intensity" ) then
   remark reflections with Iobs/sigma_I < &sigma_cut rejected
   remark reflections with Iobs > &obs_rms * rms(Iobs) rejected
 else
   remark reflections with |Fobs|/sigma_F < &sigma_cut rejected
   remark reflections with |Fobs| > &obs_rms * rms(Fobs) rejected
 end if
 xray anomalous=? end
 if ( $result = true ) then
   remark anomalous diffraction data was input
 end if
 
 {- MODIFIED 2/15/06 -}
 remark fft gridding factor = $fft_grid, B factor offset = $fft_b_add A^2, Elimit = $fft_elim
 {- END MODIFICATION -}
 
 remark theoretical total number of refl. in resol. range:    $total_theor[I6] ( 100.0 % )
 remark number of unobserved reflections (no entry):          $unobserved[I6] ( $per_unobs[f5.1] % )
 remark number of reflections rejected:                       $rejected[I6] ( $per_reject[f5.1] % )
 remark total number of reflections used:                     $total_used[I6] ( $per_used[f5.1] % )
 remark number of reflections in working set:                 $total_work[I6] ( $per_work[f5.1] % )
 remark number of reflections in test set:                    $total_test[I6] ( $per_test[f5.1] % )

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

 ! 
 ! Begin modification (6/28/06)
 if ( &bscale = "anisotropic" ) then
   REMARK Anisotropic B-factor tensor Ucart of atomic model without isotropic component :
   REMARK   B11=$Baniso_11[f8.3] B22=$Baniso_22[f8.3] B33=$Baniso_33[f8.3]
   REMARK   B12=$Baniso_12[f8.3] B13=$Baniso_13[f8.3] B23=$Baniso_23[f8.3]
   REMARK Isotropic component added to coordinate array B: $Biso_model[f8.3]
 elseif ( &bscale = "isotropic" ) then
   REMARK B-factor applied to coordinate array B: $Biso_model[f8.3]
 else
   REMARK initial B-factor correction: none
 end if
 ! End modification
 !

 {- BEGIN MODIFICATION -}
 if (&bsharp # 0) then
   remark B-factor sharpening applied to map: exp( Bsharp * S^2/4 ) where Bsharp = &bsharp
 end if
 {- END MODIFICATION -}

 remark a= &a b= &b c= &c alpha= &alpha beta= &beta gamma= &gamma sg= &STRIP%sg

 if (&map_scale=true) then
   xray
     show rms (real(map)) ( all )
     do (map=map/$result) ( all )
     remarks map has been scaled by 1/rms (rms= $result[F9.3] )
   end
 end if

 if ( &write_map = true ) then

   evaluate ($map_outfile=&output_root + ".map")

   if ( &map_mode = "asymmetric" ) then
     evaluate ($map_mode_string=ASYM)
   elseif ( &map_mode = "unit" ) then
     evaluate ($map_mode_string=UNIT)
   elseif ( &map_mode = "box" ) then
     evaluate ($map_mode_string=BOX)
   elseif ( &map_mode = "fract" ) then
     evaluate ($map_mode_string=FRAC)
   else
     evaluate ($map_mode_string=MOLE)
   end if
   xray
     write map
       if ( &map_format = "ezd" ) then
         type=ezd
       else
         type=cns
       end if
       automatic=false
       from=map
       output=$map_outfile
       cushion=&map_cushion
       selection=&atom_map
       extend=$map_mode_string
       if ( &map_mode = "box" ) then
         xmin=&xmin xmax=&xmax
         ymin=&ymin ymax=&ymax
         zmin=&zmin zmax=&zmax
       end if
       if ( &map_mode = "fract" ) then
         xmin=&xmin xmax=&xmax
         ymin=&ymin ymax=&ymax
         zmin=&zmin zmax=&zmax
       end if
     end
   end
 end if

 if ( &peak_search = true ) then
 
   show sum ( x ) ( segid PEAK) 
   if ($SELECT > 0) then
      display
      display  WARNING: there are atoms that have the reserved segid PEAK. 
      display           They will be deleted prior to the peak search.
      display
   end if
   delete sele=(segid=PEAK) end
 
   evaluate ($filename=&output_root + "_positive.peaks")
   xray
     peakpik
       from=map
       mpeak=&peak_num
       selection=( all )
       atom=true
       proximity=(&atom_map)
     end
   end

   write coor output=$filename selection=(segid=PEAK) end

   delete sele=(segid=PEAK) end

   evaluate ($filename=&output_root + "_negative.peaks")
   xray
     do (map=-map) ( all )
     peakpik
       from=map
       mpeak=&peak_num
       selection=( all )
       atom=true
       proximity=(&atom_map)
     end
   end

   write coor output=$filename selection=(segid=PEAK) end
 end if

 stop