{+ file: water_pick.inp +}
{+ directory: xtal_refine +}
{+ description: pick water molecules in electron density map +}
{+ 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))
               (u m|Fo| e^(i phi_comb)) - (v D|Fc| e^(i phi_calc))
               (u m_obs|Fo| e^(i phi_obs )) - (v k|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 and Paul D. Adams +}
{+ copyright: Yale University +}

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

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

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

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

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

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

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

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

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

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

{* space group *}
{* use International Table conventions with subscripts substituted
   by parenthesis *}
{===>} sg="P2(1)2(1)2(1)";

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

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

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

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

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

{* reciprocal space array containing test set for cross-validation:
   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.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 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
         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|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";

{* refinement target - used for gradient maps and refinement of waters *}
{+ 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;

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

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

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

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

{* 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 *}
{* only peaks closer than the maximum hbond distance and further away
   than the minimum hbond distance from oxygen or nitrogen 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.mtf";

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

{===========================================================================}
{         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;
                                sgparam=$sgparam;)

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

   @CNS_XRAYLIB:scatter.lib

   binresolution &low_res &high_res
   mapresolution &high_res

   generate &low_res &high_res

   evaluate ($counter=1)
   evaluate ($done=false)
   while ( $done = false ) loop read
    if ( &exist_reflection_infile_$counter = true ) then
      if ( &BLANK%reflection_infile_$counter = false ) then
       reflection
         @@&reflection_infile_$counter
       end
      end if
   else
     evaluate ($done=true)
   end if
   evaluate ($counter=$counter+1)
  end loop read

 end

 if ( &BLANK%anom_library = false ) then
   @@&anom_library
 else
   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=ref_active domain=reciprocal type=integer end
   declare name=tst_active domain=reciprocal type=integer end

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

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

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

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

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

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

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

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

   associate fcalc ( &atom_select )

   tselection=( ref_active=1 )

   cvselection=( tst_active=1 )

   method=FFT

 {- MODIFIED 2/15/06 -}
 end
 

 show min ( b ) ( &atom_select )
 evaluate ($b_min=$result)
 @@CNS_XTALMODULE:fft_parameter_check ( 
                             d_min=&high_res; 
                             b_min=$b_min;
                             grid=&grid;
                             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)
   elseif (&ncs_type = "both" ) 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

 @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;
                             
                             )
 xray
   @@CNS_XTALMODULE:calculate_r (
                                 fobs=&STRIP%obs_f;
                                 fcalc=fcalc;
                                 fpart=fbulk;
                                 sel=( ref_active=1 );
                                 sel_test=( tst_active=1 );
                                 print=true;
                                 output=OUTPUT;
                                 r=$start_r;
                                 test_r=$start_test_r;)
 end

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

 xray
   if ( &map_type = "gradient" ) then
     @@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;)
     predict
        mode=dtarget(fcalc)
        to=dtarg
        selection=( ref_active=1 )
        atomselection=( &atom_select ) 
     end
   end if
 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=(ref_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=(ref_active=1);
                                  sel_test=(tst_active=1);
                                  fom=mod_fom;
                                  x=mod_x;
                                  pa=mod_pa;
                                  pb=mod_pb;
                                  pc=mod_pc;
                                  pd=mod_pd;
                                  dd=mod_dd; )

     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

 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=(ref_active=1);
                                  sel_test=(tst_active=1);
                                  fom=mod_fom;
                                  x=mod_x;
                                  pa=mod_pa;
                                  pb=mod_pb;
                                  pc=mod_pc;
                                  pd=mod_pd;
                                  dd=mod_dd; )

     @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)) (&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

 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=(ref_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)
     undeclare name=m domain=reciprocal end
   end

 end if

 if ( &map_type = "gradient" ) then
 
   xray
     declare name=map domain=real end

     {- BEGIN MODIFICATION -}
     {- 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 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))))
          (ref_active=1 and acentric)
       do (fmap=    (&u   map_fom combine(amplitude(&STRIP%obs_f),map_phase)) -
                    (&v map_scale        (fcalc+fbulk)))
          (ref_active=1 and centric)
     else
       do (fmap=(&u             map_fom combine(amplitude(&STRIP%obs_f),map_phase)) -
                (&v           map_scale        (fcalc+fbulk)))
          (ref_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)))
          (ref_active=1 and centric)
       if ( &fill_in = true ) then
         do (fmap=(&u-&v) map_scale (fcalc+fbulk))
            (&low_res >= d >= &high_res and ref_active # 1)
       end if
     end if
   end

   {- MODIFIED -}
   {- 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)) (ref_active=1)
     else
       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 if
   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

 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

 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

 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

 {- BEGIN MODIFICATION -}
 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 -}

 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

 {- BEGIN MODIFICATION -}
 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 -}

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

 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 )

 parameter
   nbonds wmin=1.5 end
 end

 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

   @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=&grid;
                               fft_memory=&fft_memory;
                               fft_grid=$fft_grid;   
                               fft_b_add=$fft_b_add; 
                               fft_elim=$fft_elim; 
                                        )

   xray
     @@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

   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
       
       method=lbfgs
       
       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
       
         method=lbfgs

         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 (fobs=&STRIP%obs_f;
                                fcalc=fcalc;
                                fpart=fbulk;
                                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=adjustl(format("I4",$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_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
 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 a oxygen or nitrogen were deleted
 end if
 display REMARK map resolution: &low_res - &high_res A  
 if ( $total_test > 0 ) then
   display REMARK starting r= $start_r[f6.4] free_r= $start_test_r[f6.4]
   display REMARK final    r= $full_r[f6.4] free_r= $full_test_r[f6.4]
 else
   display REMARK starting r= $start_r[f6.4]
   display REMARK final    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

 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
      display REMARK parameter file $counter  : &STRIP%parameter_infile_$counter
    end if
  else
   evaluate ($done=true)
  end if
  evaluate ($counter=$counter+1)
 end loop read
 
 if ( &BLANK%structure_infile = true ) then
   display REMARK molecular structure file: automatic

   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
        display REMARK topology file $counter  : &STRIP%topology_infile_$counter
      end if
    else
     evaluate ($done=true)
    end if
    evaluate ($counter=$counter+1)
   end loop read

   evaluate ($counter=1)
   evaluate ($done=false)
   while ( $done = false ) loop read
    if ( &exist_link_infile_$counter = true ) then
      if ( &BLANK%link_infile_$counter = false ) then
        display REMARK linkage file $counter  : &STRIP%link_infile_$counter
      end if
    else
     evaluate ($done=true)
    end if
    evaluate ($counter=$counter+1)
   end loop read

   if ( &BLANK%patch_infile = false ) then
      display REMARK custom patch file = &STRIP%patch_infile
   end if

 else
   display REMARK molecular structure file: &STRIP%structure_infile
 end if
 
 display REMARK input coordinates: &STRIP%coordinate_infile
 if ( &BLANK%anom_library = false ) then
   display REMARK anomalous f' f'' library: &STRIP%anom_library
 end if

 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
         display REMARK reflection file $counter : &STRIP%reflection_infile_$counter       
      end if
   else
     evaluate ($done=true)
   end if
   evaluate ($counter=$counter+1)
 end loop read

 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 (6/28/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
 ! End modification
 !

 {- MODIFIED 5/18/05 -}
 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 -}
 
  {- BEGIN MODIFICATION -}
 if (&bsharp # 0) then
   display REMARK B-factor sharpening applied to map: exp( Bsharp * S^2/4 ) where Bsharp = &bsharp
 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
 else
   display REMARK reflections with |Fobs|/sigma_F < &sigma_cut rejected
   display REMARK reflections with |Fobs| > &obs_rms * rms(Fobs) rejected
 end if
 xray anomalous=? end
 if ( $result = true ) then
   display REMARK anomalous diffraction data was input
 end if

 {- MODIFIED 2/15/06 -}
 display REMARK 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):          $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=true;
                            coordinate_outfile=&coordinate_outfile;
                            sgparam=$sgparam;)

 stop