{+ file: ncs_average_map.inp +}
{+ directory: xtal_refine +}
{+ description: Make an annealed omit map followed by a single cycle of
                electron density averaging +}
{+ 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 is the calculated structure factor and
              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, 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: 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
   - the selections store1 through store4 are available for general use -}

{- begin block parameter definition -} define(

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

{* coordinate file *}
{===>} coordinate_infile="eg1_dimer.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="eg1_dimer.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="P4(1)2(1)2";

{* unit cell parameters in Angstroms and degrees *}
{+ table: rows=1 "cell" cols=6 "a" "b" "c" "alpha" "beta" "gamma" +}
{===>} a=101.4;
{===>} b=101.4;
{===>} c=199.5;
{===>} 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="eg1_abcd.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=4.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="ncs.def";

{============ 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 *}
{* this should include all conformations if multiple conformations are used *}
{===>} atom_select=(known and not hydrogen);

{* 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=(segid BBBB or segid DDDD);

{* select atoms to be omitted from map calculation *}
{===>} atom_omit=(residue 85:90);

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

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

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

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

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

{=========================== NCS masks ===============================}

{* masks can be derived from the input atomic coordinates, or can be
   read in from an O format mask file. One or the other - do not mix *}

{* masks from coordinates: select atoms which form the primary protomer 
   in each NCS group. In the strict NCS case only one group can be 
   defined. *}

{* select atoms in NCS group 1 *}
{===>} ncs_group_1=(none);

{* select atoms in NCS group 2 *}
{===>} ncs_group_2=(none);

{* select atoms in NCS group 3 *}
{===>} ncs_group_3=(none);

{* select atoms in NCS group 4 *}
{===>} ncs_group_4=(none);

{* masks from file: give the O format mask files which form the primary 
   protomer in each NCS group. In the strict NCS case only one mask
   can be defined. *}

{* NCS group 1: O format mask file *}
{===>} mask_infile_1="eg1.mask";

{* NCS group 2: O format mask file *}
{===>} mask_infile_2="";

{* NCS group 3: O format mask file *}
{===>} mask_infile_3="";

{* NCS group 4: O format mask file *}
{===>} mask_infile_4="";

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

{* map grid size: dmin*grid *}
{* 0.25 is reasonable for map averaging *}
{===>} grid=0.25;

{* 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 *}
{+ choice: "torsion" "cartesian" +}
{===>} md_type="torsion";

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

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

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

{======================= 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 file name for output files *}
{+ list: files created:
         output NCS averaged map file in: <output_root>.map
         positive peaks in: <output_root>_positive.peaks
         negative peaks in: <output_root>_negative.peaks
         listing file in: <output_root>.list
         Fourier coefficients will be in: <output_root>.coeff +}
{===>} output_root="ncs_average_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

 {- MODIFICATION: removed "refine_low_res" parameter -}
 
 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 )
   
   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 

   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
   evaluate ($test_hl=true)
 elseif ( &map_type = "combined" ) then
   evaluate ($test_hl=true)
 else
   evaluate ($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 ) )

 ident ( store9 )  ( byresidue ( ( store8 ) saround &cushion and 
                               not store8 ) )

 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
   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=auto;
                             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 ( &mini_steps > 0 ) then
   if ( &wa < 0 ) then
     @@CNS_XTALMODULE:getweight (selected=&atom_select;
                                 fixed=(store7);)
   end if
   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
   @@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
     {- take the negative of the gradient so the map is the same sign
        as a standard difference map -}
     do (map=ft(-dtarg)) (map_active=1)
   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

   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

 evaluate ($1=1)
 evaluate ($2=2)
 evaluate ($3=3)
 evaluate ($4=4)

 evaluate ($ngroup=1)
 evaluate ($group=1)
 evaluate ($done=false)
 while ( $done = false ) loop group
   if ( &exist_ncs_group_$group = true ) then
     show sum(1) ( &ncs_group_$group )
     if ( $result > 0 ) then
       evaluate ($ngroup=$ngroup+1)
     end if
   else
     evaluate ($done=true)
     evaluate ($ngroup=$ngroup-1)
   end if
   evaluate ($group=$group+1)
 end loop group

 if ( $ngroup <= 0 ) then
   evaluate ($ngroup=1)
   evaluate ($done=false)
   while ( $done = false ) loop group
     if ( &exist_mask_infile_$ngroup = true ) then
       if ( &BLANK%mask_infile_$ngroup = true ) then
         evaluate ($done=true)
         evaluate ($ngroup=$ngroup-1)
       else
         evaluate ($ngroup=$ngroup+1)
       end if
     else
         evaluate ($done=true)
         evaluate ($ngroup=$ngroup-1)
     end if
   end loop group
 end if

 if ($ngroup <= 0) then
   display No NCS groups found, aborting
   abort
 end if

 if (&ncs_type="strict") then
   ncs strict ? end
   evaluate ($num_op_1=$ncs)
 elseif (&ncs_type="restrain") then
   ncs restraint ? end
   evaluate ($group=1)
   while ($group <= $ngroup) loop group
     evaluate ($num_op=1)
     evaluate ($done=false)
     while ( $done = false ) loop ncsop
       if ( $exist_rot_$group_$num_op_$1_$1 # true ) then
         evaluate ($done=true)
         evaluate ($num_op=$num_op-1)
       else
         evaluate ($num_op=$num_op+1)
       end if
     end loop ncsop
     evaluate ($num_op_$group=$num_op)
     evaluate ($group=$group+1)
   end loop group
 {- modification, ATB, 12/20/08 -}
 elseif (&ncs_type="both") then
   display this script file currently does not support both NCS restraints and constraints at the same time
   abort
 else
   display unknown NCS type
   abort
 end if

 evaluate ($group=1)
 while ($group <= $ngroup) loop mask
   evaluate ($maskname="mask" + encode($group))
   if ( &BLANK%mask_infile_$group = false ) then
     xray
       declare name=$maskname domain=real end
       read mask
         to=$maskname
         type=omask
         input=&mask_infile_$group
       end
     end
   else
!
! modification ATB 5/29/08.  For strict NCS and automatic mask generation
! we need to compute the mask for the primary molecule only.  However, the
! mask statement computes the mask over all molecules (primary + strict).  
! We work around this problem by temporarily reseting the NCS operators
!
     if (&ncs_type="strict") then
      ncs strict init end
      xray
        declare name=$maskname domain=real end
        mask
         average=true
         mode=vdw
         solrad=1.0
         shrink=1.0
         nshell=1
         to=$maskname
         sele=( &ncs_group_$group )
        end
      end

     if ( &BLANK%ncs_infile = false ) then
      inline @&ncs_infile
     end if
     ncs strict ? end
     
     else
      xray
        declare name=$maskname domain=real end
        mask
         average=true
         mode=vdw
         solrad=1.0
         shrink=1.0
         nshell=1
         to=$maskname
         sele=( &ncs_group_$group )
        end
      end
     end if
   end if
   evaluate ($group=$group+1)
 end loop mask 

 evaluate ($display=&output_root + ".list")
 set display=$display end
 display ================ Summary of single-pass NCS density averaging ================
 display
 display input data ---->
 display >>> resolution: &low_res - &high_res A
 display >>> sg= &STRIP%sg a= &a b= &b c= &c alpha= &alpha beta= &beta gamma= &gamma

 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 >>> reflection file $counter : &STRIP%reflection_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_parameter_infile_$counter = true ) then
    if ( &BLANK%parameter_infile_$counter = false ) then
      display >>> 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 >>> 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 >>> 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 >>> 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 >>> custom patch file = &STRIP%patch_infile
   end if

 else
   display >>> molecular structure file: &STRIP%structure_infile
 end if
 
 display >>> input coordinates: &STRIP%coordinate_infile
 if ( &BLANK%anom_library = false ) then
   display >>> anomalous f' f'' library: &STRIP%anom_library
 end if
 if ( &BLANK%restraints_infile = false ) then
   display >>> additional restraints file: &STRIP%restraints_infile
 end if
 display >>> ncs= &STRIP%ncs_type  ncs file= &STRIP%ncs_infile
 display
 display data usage ---->
 if ( &obs_type = "intensity" ) then
   display >>> reflections with Iobs/sigma_I < &sigma_cut rejected
   display >>> reflections with Iobs > &obs_rms * rms(Iobs) rejected
 else
   display >>> reflections with |Fobs|/sigma_F < &sigma_cut rejected
   display >>> reflections with |Fobs| > &obs_rms * rms(Fobs) rejected
 end if
 
 {- MODIFIED 2/15/06 -}
 display >>> fft gridding factor = $fft_grid, B factor offset = $fft_b_add A^2, Elimit = $fft_elim
 {- END MODIFICATION -}
 
 display >>> theoretical total number of refl. in resol. range:    $total_theor[I6] ( 100.0 % )
 display >>> number of unobserved reflections (no entry):          $unobserved[I6] ( $per_unobs[f5.1] % )
 display >>> number of reflections rejected:                       $rejected[I6] ( $per_reject[f5.1] % )
 display >>> total number of reflections used:                     $total_used[I6] ( $per_used[f5.1] % )
 display >>> number of reflections in working set:                 $total_work[I6] ( $per_work[f5.1] % )
 display
 display NCS information ---->
 display >>> number of NCS groups= $ngroup
 evaluate ($group=1)
 while ($group <= $ngroup) loop group
   display >>> NCS group $group : number of NCS operators= $num_op_$group
   evaluate ($group=$group+1)
 end loop group

 xray
   evaluate ($group=1)
   while ($group <= $ngroup) loop average
     evaluate ($maskname="mask" + encode($group))
     average
       from=map
       evaluate ($ncsop=1)
       while ($ncsop <= $num_op_$group) loop ncsop
         group
           if ( $ncsop = 1 ) then
             mask=$maskname
           end if
           if (&ncs_type="strict") then
             matrix=
               ( $ncsop_$ncsop_$1_$1 $ncsop_$ncsop_$1_$2 $ncsop_$ncsop_$1_$3 )
               ( $ncsop_$ncsop_$2_$1 $ncsop_$ncsop_$2_$2 $ncsop_$ncsop_$2_$3 )
               ( $ncsop_$ncsop_$3_$1 $ncsop_$ncsop_$3_$2 $ncsop_$ncsop_$3_$3 )
             translation=
               ( $ncsop_$ncsop_$1_$4 $ncsop_$ncsop_$2_$4 $ncsop_$ncsop_$3_$4 )
           elseif (&ncs_type="restrain") then
             matrix=
               ( $rot_$group_$ncsop_$1_$1 
                 $rot_$group_$ncsop_$1_$2 
                 $rot_$group_$ncsop_$1_$3 )
               ( $rot_$group_$ncsop_$2_$1 
                 $rot_$group_$ncsop_$2_$2 
                 $rot_$group_$ncsop_$2_$3 )
               ( $rot_$group_$ncsop_$3_$1 
                 $rot_$group_$ncsop_$3_$2 
                 $rot_$group_$ncsop_$3_$3 )
             translation=
               ( $rot_$group_$ncsop_$1_$4 
                 $rot_$group_$ncsop_$2_$4 
                 $rot_$group_$ncsop_$3_$4 )
           end if
         end
         evaluate ($ncsop=$ncsop+1)
       end loop ncsop
     end
     display     NCS group $group :     average correlation= $av_corr[f6.4]
     evaluate ($ncsop=1)
     while ($ncsop <= $num_op_$group) loop print
       display       NCS operator $ncsop ( 1 -> $ncsop ): correlation= $av_corr_op_$ncsop[f6.4]
       evaluate ($ncsop=$ncsop+1)
     end loop print
     evaluate ($group=$group+1)     
   end loop average
 end

 set remarks=reset end
 set remarks=accumulate end

 evaluate ($remark="")
 evaluate ($remark="NCS averaged ")

 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="NCS averaged")
 if ( &map_type = "sigmaa" ) then
   if ( $test_exist = true ) then
     evaluate ($remark=$remark + " cross-val.")
   end if
   evaluate ($remark=$remark + " sigmaa")
 end if
 if ( &map_type = "combined" ) then
   if ( $test_exist = true ) then
     evaluate ($remark=$remark + " cross-val.")
   end if
   evaluate ($remark=$remark + " sigmaa")
 end if
 evaluate ($remark=$remark + " omit map")

 remark $remark

 if ( &map_type = "unweighted" ) then
   evaluate ($remark="(" + encode(&u) +  " |Fo| - " +
                           encode(&v) + " k|Fc|)e^(i phi_calc)")
 elseif ( &map_type = "sigmaa" ) then
   evaluate ($remark="(" + encode(&u) + " m|Fo| - " +
                           encode(&v) + " D|Fc|)e^(i phi_calc)")
 elseif ( &map_type = "combined" ) then
   evaluate ($remark="(" + encode(&u) + " m|Fo|)e^(i phi_comb) - " +
                     "(" + encode(&v) + " D|Fc|)e^(i phi_calc)")
 elseif ( &map_type = "observed" ) then
   evaluate ($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="( d(" + &reftarget + ")/dFc )")
 end if

 remark $remark

 if ( $total_test > 0 ) then
   remark r= $map_r[f6.4] free_r= $map_free_r[f6.4]
 else
   remark r= $map_r[f6.4]
 end if

 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
 !

 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

   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

   evaluate ($output_map=&output_root + ".map")
   xray
     write map
       if ( &map_format = "ezd" ) then
         type=ezd
       else
         type=cns
       end if
       automatic=false
       from=map
       output=$output_map
       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