{+ file: composite_omit_map.inp +}
{+ directory: xtal_refine +}
{+ description: Make a composite annealed omit map +}
{+ 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 m_obs k|Fc| e^(i phi_calc))
               d(target)/dFc
              where Fc is the calculated structure factor and
              m and D are derived from sigmaa
            Averaging of structure factors from several models optional.
            Realspace R-value calculation optional. +}
{+ 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: T.N. Bhat, Calculation of an OMIT map. J Appl Crystallogr
              21, 279-281 (1988) +}
{+ reference: A. Hodel, S.-H. Kim, and A.T. Brunger, Model Bias in 
              Macromolecular Crystal Structures, Acta Cryst. A48, 
              851-859 (1992) +}
{+ reference: L.M. Rice and A.T. Brunger, Torsion Angle Dynamics:
              Reduced Variable Conformational Sampnameling Enhances
              Crystallographic Structure Refinement, Proteins: Structure,
              Function, and Genetics, 19, 277-290 (1994) +}
{+ reference: A.T. Brunger, P.D. Adams and L.M. Rice, New applications
              of simulated annealing in X-ray crystallography and
              solution NMR, Structure 5, 325-336, (1997) +}

{- Guidelines for using this file:
   - all strings must be quoted by double-quotes
   - logical variables (true/false) must not be quoted
   - do not remove any evaluate statements from the file
   - the selections store1 through store6 are available for general use -}

{- begin block parameter definition -} define(

{======================= molecular structure =========================}

{* molecular topology file *}
{===>} structure_infile="amy.mtf";

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

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

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

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

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

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

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

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

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

{* reciprocal space array containing test set for cross-validation: required *}
{* cross-validation should always be used, with the possible exception
   of a final round of refinement including all data *}
{* cross-validation is always required for the maximum likelihood targets *}
{===>} test_set="test";

{* number for selection of test reflections: required for cross-validation *}
{* ie. reflections with the test set array equal to this number will be
       used for cross-validation, all other reflections form the working set *}
{===>} test_flag=1;

{* reciprocal space array containing weighting scheme for observed
   amplitudes: optional *}
{* only used for the "residual" and "vector" targets - this will
   default to a constant value of 1 if array is not present *}
{===>} obs_w="";

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

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

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

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

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

{* resolution limits for data included in map calculation *}
{* all data available should be included in the map calculation *}
{+ table: rows=1 "resolution" cols=2 "lowest" "highest" +}
{===>} low_res=500.0;
{===>} high_res=2.0;

{* apply rejection criteria to amplitudes or intensities *}
{+ choice: "amplitude" "intensity" +}
{===>} obs_type="amplitude";

{* Observed data cutoff criteria: applied to amplitudes or intensities *}
{* reflections with magnitude(Obs)/sigma < cutoff are rejected. *}
{===>} sigma_cut=0.0;

{* rms outlier cutoff: applied to amplitudes or intensities *}
{* reflections with magnitude(Obs) > cutoff*rms(Obs) will be rejected *}
{===>} obs_rms=10000;

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

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

{============ overall B-factor and bulk solvent corrections ==========}

{* overall B-factor correction *}
{+ choice: "no" "isotropic" "anisotropic" +}
{===>} bscale="anisotropic";

{* bulk solvent correction *}
{* a mask is required around the molecule(s). The region
   outside this mask is the solvent region *}
{+ choice: true false +}
{===>} bulk_sol=true;

{* bulk solvent mask file *}
{* mask will be read from O type mask file if a name is given
   otherwise calculated from coordinates of selected atoms *}
{===>} bulk_mask_infile="";

{* automatic bulk solvent parameter search *}
{+ choice: true false +}
{===>} sol_auto=true;

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

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

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

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

{* select atoms to be included in structure factor calculation *}
{* this should include all conformations if multiple conformations are used *}
{===>} atom_select=(known and not hydrogen);

{* select atoms around which map will be calculated *}
{* use this to limit the map calculation to a part of the molecule *}
{===>} atom_map=(known and not hydrogen);

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

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

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

{* select atoms to be treated as rigid groups during torsion angle dynamics *}
{===>} atom_rigid=(none);

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

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

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

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

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

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

{* maximum percentage of model to be omitted *}
{* at most this should be 10 percent *}
{===>} percent_omit=5.0;

{* maps are calculated u*Fo - v*Fc *}
{* eg. 2fo-fc map -> u=2 and v=1 *}
{* note: a fo-fc map is not useful since the whole model is omitted *}

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

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

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

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

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

{* map format *}
{* choice: "cns" "ezd" *}
{===>} map_format="cns";

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

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

{======================= 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 MD parameters *}
{* increase these values if the program terminates with the message
   that one of these parameters is exceeded *}

{* maximum unbranched chain length *}
{* increase for long stretches of polyalanine *}
{===>} torsion_maxlength=50;
{* maximum number of distinct bodies *}
{===>} torsion_maxtree=15;
{* maximum number of chains (increase for large molecules) *}
{===>} torsion_maxchain=1000;
{* maximum number of bonds to an atom *}
{===>} torsion_maxbond=6;

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

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

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

{========================= realspace R-value =========================}

{* calculate realspace R-values *}
{* this will be computed for all atoms selected for map calculation *}
{+ choice: true false +}
{===>} real_r=true;

{* select residues for which realspace R-values will be calculated *}
{===>} atom_real=(known and not hydrogen);

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

{* root name for output files *}
{+ list:
         map file will be in: <output_root>.map
         positive peaks in: <output_root>_positive.peaks
         negative peaks in: <output_root>_negative.peaks
         realspace R-values in: <output_root>_r.list
                         and in <output_root>_r.pdb
         progress of omit refinements: <output_root>.list
         Fourier coefficients will be in: <output_root>.coeff +}
{===>} output_root="composite_omit_map";

{* 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: array written:
         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.2

 {- MODIFICATION: removed "refine_low_res" parameter -}
 {- MODIFICATION: removed multiple coordinate entries 
    due to map resets by solvent mask calculation -}
 
 evaluate ($log_level=quiet)

 structure @&structure_infile end

 coordinates @&coordinate_infile

 parameter
   if ( &BLANK%parameter_infile_1 = false ) then
      @@&parameter_infile_1
   end if
   if ( &BLANK%parameter_infile_2 = false ) then
      @@&parameter_infile_2
   end if
   if ( &BLANK%parameter_infile_3 = false ) then
      @@&parameter_infile_3
   end if
   if ( &BLANK%parameter_infile_4 = false ) then
      @@&parameter_infile_4
   end if
   if ( &BLANK%parameter_infile_5 = false ) then
      @@&parameter_infile_5
   end if
 end

 xray

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

   a=&a b=&b c=&c  alpha=&alpha beta=&beta gamma=&gamma
  
   @CNS_XRAYLIB:scatter.lib

   binresolution &low_res &high_res
   mapresolution &high_res

   generate &low_res &high_res

   if ( &BLANK%reflection_infile_1 = false ) then
     reflection @@&reflection_infile_1 end
   end if
   if ( &BLANK%reflection_infile_2 = false ) then
     reflection @@&reflection_infile_2 end
   end if
   if ( &BLANK%reflection_infile_3 = false ) then
     reflection @@&reflection_infile_3 end
   end if

 end

 if ( &BLANK%anom_library = false ) then
   @@&anom_library
 else
   set echo=off end
   xray anomalous=? end
   if ( $result = true ) then
     display Warning: no anomalous library has been specified
     display          no anomalous contribution will used in refinement
   end if
   set echo=on end
 end if

 {- copy define parameters of optional arrays into symbols so 
    we can redefine them -}
    
 evaluate ($obs_i=&obs_i)
 evaluate ($obs_sigi=&obs_sigi)
 evaluate ($obs_w=&obs_w)
 xray
   @@CNS_XTALMODULE:checkrefinput (
                                  reftarget=&reftarget;
                                  obs_f=&obs_f;
                                  obs_sigf=&obs_sigf;
                                  test_set=&test_set;
                                  obs_pa=&obs_pa;
                                  obs_pb=&obs_pb;
                                  obs_pc=&obs_pc;
                                  obs_pd=&obs_pd;
                                  obs_phase=&obs_phase;
                                  obs_fom=&obs_fom;
                                  obs_w=$obs_w;
                                  obs_i=$obs_i;
                                  obs_sigi=$obs_sigi;
                                  )

   query name=fcalc domain=reciprocal end
   if ( $object_exist = false ) then
      declare name=fcalc domain=reciprocal type=complex end
   end if
   declare name=fbulk domain=reciprocal type=complex end
   do (fbulk=0) ( all )

   query name=&STRIP%obs_f domain=reciprocal end
   declare name=fobs_orig domain=reciprocal type=$object_type end
   declare name=sigma_orig domain=reciprocal type=real end

   do (fobs_orig=&STRIP%obs_f) (all)
   do (sigma_orig=&STRIP%obs_sigf) (all)

   if ( &BLANK%obs_i = false ) then
     query name=&STRIP%obs_i domain=reciprocal end
     declare name=iobs_orig domain=reciprocal type=$object_type end
     declare name=sigi_orig domain=reciprocal type=real end
     do (iobs_orig=&STRIP%obs_i) (all)
     do (sigi_orig=&STRIP%obs_sigi) (all)
   end if

   if ( &obs_type = "intensity" ) then
     if ( &BLANK%obs_i = true ) then
       display  Error: observed intensity array is undefined
       display         aborting script
       abort
     end if
     evaluate ($reject_obs=&obs_i)
     evaluate ($reject_sig=&obs_sigi)
     show min (amplitude(&STRIP%obs_i)) (all)
     evaluate ($obs_lower_limit=$result-0.1)
   else
     evaluate ($reject_obs=&obs_f)
     evaluate ($reject_sig=&obs_sigf)
     evaluate ($obs_lower_limit=0)
   end if

   declare name=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) ( ( amplitude($STRIP%reject_obs) > $obs_lower_limit ) and
                       ( &low_res >= d >= &high_res ) )

   do (map_active=0) ( all )
   do (map_active=1) ( ( amplitude($STRIP%reject_obs) > $obs_lower_limit ) and
                       ( &low_res >= d >= &high_res ) )

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

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

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

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

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

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

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

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

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

   associate fcalc ( &atom_select )

   tselection=( map_active=1 )

   cvselection=( tst_active=1 )

   method=FFT          
   
 {- MODIFIED 2/15/06 -}
 end
 

 show min ( b ) ( &atom_select )
 evaluate ($b_min=$result)
 @@CNS_XTALMODULE:fft_parameter_check ( 
                             d_min=&high_res; 
                             b_min=$b_min;
                             grid=&grid;
                             fft_memory=&fft_memory;
                             fft_grid=$fft_grid;   
                             fft_b_add=$fft_b_add; 
                             fft_elim=$fft_elim; 
                                      )
                            
 xray
 {- END MODIFICATION -}

   if ( &wa >= 0 ) then
      wa=&wa
   end if
      
 end                  

 if ( &map_type = "observed" ) then
   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

 fastnb grid end

 flags
    exclude elec include pvdw xref
    ?
 end
 
 set seed=&seed end

 do (store8=0) ( &atom_select )

 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 (store8=$alt) ( &conf_$alt )
   evaluate ($alt=$alt+1)
 end loop alt

 igroup
   interaction ( &atom_select and not(attr store8 > 0))
               ( &atom_select and not(attr store8 > 0))
   evaluate ($alt=1)
   while ( $alt <= $nalt ) loop alcs
     interaction ( &atom_select and ( attr store8 = $alt or attr store8 = 0 ))
                 ( &atom_select and ( attr store8 = $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 - store8 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=store8;
                           atom_multiplicity=rmsd;
                           mset=$mset;
                           )
 
 fix selection=( store8 ) end

 xray
   associate fcalc ( &atom_select )
 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

 if ( &anneal = true ) then
   if ( &md_type = "torsion" ) then
     dyna tors
       topology
       
         maxlength=&torsion_maxlength
         maxchain=&torsion_maxchain
         maxtree=&torsion_maxtree
         maxbond=&torsion_maxbond
         kdihmax = 95.
         
         @CNS_TOPPAR:torsionmdmods
         
         fix group ( &atom_rigid )
         
       end
        nstep=0
        cmremove=true
     end
   end if
 end if

 xray
   do (&STRIP%obs_f=fobs_orig) (all)
   do (&STRIP%obs_sigf=sigma_orig) (all)
   if ( &BLANK%obs_i = false ) then
     do (&STRIP%obs_i=iobs_orig) (all)
     do (&STRIP%obs_sigi=sigi_orig) (all)
   end if
 end

 xray
   predict
     mode=reciprocal
     to=fcalc
     selection=( all_active=1 )
     atomselection=( &atom_select ) 
   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 -}
 
 
 {- BEGIN MODIFICATION - MOVED AFTER SOLVENT MASK COMPUTATION -}
 xray
   declare name=map_total  domain=real end
   declare name=omit_mask  domain=real end
   declare name=map        domain=real end  
   declare name=fmap domain=reciprocal type=complex end    
   declare name=dtarg domain=reciprocal type=complex end
   
   do (map_total=0) ( all )
 end
{- 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

 do ( harm=0 ) ( all )
 do ( harmonic=&k_harmonic ) ( &atom_harm )
 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=(store8);
                                )
 end if

 { get the box that covers the molecule in fractional coordinates }
 
 coor fract end

 show sum (1) (&atom_select)
 eval ($natom=$result)

 show sum (1) (&atom_select and not hydrogen)

 evaluate ($volume=$result*9.0)
 evaluate ($cube_length=($volume*(&percent_omit/100))**(1/3))

 evaluate ($cushion=0.1*$cube_length)

 { wrap the molecule into the primary unit cell }
 do (x=mod(x+10,1)) ( all )
 do (y=mod(y+10,1)) ( all )
 do (z=mod(z+10,1)) ( all )

 evaluate ($map_cushion=2.5)
 
 show max ( x ) ( &atom_map and &atom_select )
 eval ($xmax=$result+($map_cushion/&a))
 show min ( x ) ( &atom_map and &atom_select )
 eval ($xmin=$result-($map_cushion/&a))
 show max ( y ) ( &atom_map and &atom_select )
 eval ($ymax=$result+($map_cushion/&b))
 show min ( y ) ( &atom_map and &atom_select )
 eval ($ymin=$result-($map_cushion/&b))
 show max ( z ) ( &atom_map and &atom_select )
 eval ($zmax=$result+($map_cushion/&c))
 show min ( z ) ( &atom_map and &atom_select )
 eval ($zmin=$result-($map_cushion/&c))
 
 evaluate ($display=&output_root + ".list")
 set display=$display end
 display ========================================================
 display maximum percentage of model to omit= &percent_omit %
 display cube length= $cube_length Angstroms
 display cushion= $cushion Angstroms
 display limits of overall map box in fractional space:
 display    xmin= $xmin[f7.4]  xmax=$xmax[f7.4]
 display    ymin= $ymin[f7.4]  ymax=$ymax[f7.4]
 display    zmin= $zmin[f7.4]  zmax=$zmax[f7.4]
 display ========================================================
 set display=OUTPUT end

 do ( x=refx ) ( all )
 do ( y=refy ) ( all )
 do ( z=refz ) ( all )
 
 if ( $log_level = verbose ) then
   set message=normal echo=on end
 else
   set message=off echo=off end
 end if

 eval ($total=0)
 eval ($anneal_total=0)

 eval ($x=$xmin)
 while ($x < $xmax) loop xx
   eval ($x1=$x-$cushion/&a)
   eval ($x2=$x+($cube_length+$cushion)/&a)
   eval ($xstop=$x+$cube_length/&a)
 
   eval ($y=$ymin)
   while ($y < $ymax) loop yy
     eval ($y1=$y-$cushion/&b)
     eval ($y2=$y+($cube_length+$cushion)/&b)
     eval ($ystop=$y+$cube_length/&b)
 
     eval ($z=$zmin)
     while ($z < $zmax) loop zz
       eval ($z1=$z-$cushion/&c)
       eval ($z2=$z+($cube_length+$cushion)/&c)
       eval ($zstop=$z+$cube_length/&c)

       {- retrieve coordinates -}

       do ( x=refx ) ( all )
       do ( y=refy ) ( all )
       do ( z=refz ) ( all )

       {- check if there are primary molecule atoms in the box+cushion -}
       show sum (1) ( ( fbox $x1 $x2 $y1 $y2 $z1 $z2 ) and 
                        &atom_select )

       if ($select > 0) then

         show sum (1) ( ( fbox $x $xstop $y $ystop $z $zstop ) and 
                          &atom_select )
         eval ($natom_prim=$select)
 
         show sum (1) ( ( sfbox $x $xstop $y $ystop $z $zstop ) and 
                          &atom_select )
         eval ($natom_box=$select)
 
         eval ($anneal_total=$anneal_total+1)

         {- we need to omit all symmetry and lattice translated atoms 
            that are in the box+cushion -}
         {- de-select omitted atoms -}
         xray
           associate fcalc 
               (not(sfbox $x1 $x2 $y1 $y2 $z1 $z2 ) and &atom_select)    
           eval ($natom_box_cushion=$natom-$select)
           eval ($percent=$natom_box_cushion/$natom * 100)
         end
 
         {- harmonically restrain omitted atoms -}
         do ( harmonic=0) ( all )
         do ( harmonic=&k_harmonic ) ( &atom_harm )
         do ( harmonic=20) (sfbox $x1 $x2 $y1 $y2 $z1 $z2 ) 
 
         if ( &anneal = true ) then
 
           xray
             tolerance=0.2 lookup=true
           end
 
           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
 
           do (fbeta=$fbeta) ( ( &atom_select ) and not store8 )

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

           do (store7=mass) ( all )
           do (mass=max(10,min(30,mass))) ( all )

           if ( &md_type = "torsion" ) then
             dynamics torsion
               reassign=true
               timestep=$time_step
               nstep=$md_steps
               nprint=5
               cmremove=false
               vscaling=true  
               temperature=&temperature
             end
           end if

           evaluate ( $curr_temp = &temperature )
           while ( $curr_temp > 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=store7) ( all )
 
           parameter
             nbonds
               repel=$repel_old rcon=$rcon_old
             end
           end
 
         end if
 
         xray
           tolerance=0.0 lookup=false
         end
  
         if (&mini_steps > 0 ) then 
           minimize lbfgs
             nstep=&mini_steps
             drop=10.0 nprint=5
           end
         end if

       else

         evaluate ($natom_prim=0)

         show sum (1) ( ( sfbox $x $xstop $y $ystop $z $zstop ) and 
                          &atom_select )
         eval ($natom_box=$select)

         show sum (1) ( ( sfbox $x1 $x2 $y1 $y2 $z1 $z2 ) and 
                          &atom_select )
         evaluate ($natom_box_cushion=$select)
 
       end if

       xray

         tselection=(map_active=1)

         if ( &map_type = "gradient" ) then
           predict
               mode=dtarget(fcalc)
               to=dtarg
               selection=( map_active=1 )
               atomselection=(not(sfbox $x1 $x2 $y1 $y2 $z1 $z2 ) and 
                                  &atom_select)
           end
         else
            predict
               mode=reciprocal
               to=fcalc
               selection=( all_active=1 )
               atomselection=(not(sfbox $x1 $x2 $y1 $y2 $z1 $z2 ) and 
                                  &atom_select)
            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)
      
           declare name=total domain=reciprocal type=complex end

           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)

           undeclare name=total domain=reciprocal end
      
           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="expt. 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
      
           declare name=total domain=reciprocal type=complex end

           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)

           undeclare name=total domain=reciprocal end
      
           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
                 
           {- BEGIN MODIFICATION -}
           if (&bsharp # 0) then 
              {- b-factor sharpening -} 
               evaluate ($bsharp= (-1) * &bsharp) 
              {- take the negative of the gradient so the map is the same sign
                 as a standard difference map -}
              do (fmap=exp( $bsharp * s^2/4) * (-1) * dtarg) ( all )
 
              do (map=ft(fmap)) (map_active=1)
           end if
           {- END MODIFICATION -}
           
         end
      
       else
      
         xray
           if ( &u = &v ) then
             do (fmap= 2 ((&u map_fom combine(amplitude(&STRIP%obs_f),map_phase)) -
                          (&v map_scale        (fcalc+fbulk))))
                (map_active=1 and acentric)
             do (fmap=    (&u map_fom combine(amplitude(&STRIP%obs_f),map_phase)) -
                          (&v map_scale        (fcalc+fbulk)))
                (map_active=1 and centric)
           else
             do (fmap=(&u    map_fom combine(amplitude(&STRIP%obs_f),map_phase)) - 
                      (&v    map_scale        (fcalc+fbulk)))
                (map_active=1 and acentric)
             do (fmap=(max((&u-1),0) map_fom 
                       combine(amplitude(&STRIP%obs_f),map_phase)) -
                      (max((&v-1),0) map_scale        (fcalc+fbulk)))
                (map_active=1 and centric)
             if ( &fill_in = true ) then
               do (fmap=(&u-&v) map_scale (fcalc+fbulk)) 
                  (all_active=1 and map_active # 1)
             end if
           end if
         end 

        {- BEGIN MODIFICATION -}
        {- b-factor sharpening -} 
        evaluate ($bsharp= (-1) * &bsharp) 
        xray
           do (fmap=exp( $bsharp * s^2/4) * fmap) ( all )
        end
        {- END MODIFICATION -}
       
         xray
           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
      
       xray
         undeclare name=map_phase domain=reciprocal end 
         undeclare name=map_fom   domain=reciprocal end 
         undeclare name=map_scale domain=reciprocal end 
       end
      
       xray
         mask
           mode=fbox
           to=omit_mask
           xmin=$x xmax=$xstop
           ymin=$y ymax=$ystop
           zmin=$z zmax=$zstop
           selection=( &atom_select ) 
         end
         do (map_total=map)     (real(omit_mask)<=0)
       end

       evaluate ($display=&output_root + ".list")
       set display=$display end

       display ========================================================
       display   map box limits (fractional): x=$x[F8.4] to $xstop[F8.4]
       display                                y=$y[F8.4] to $ystop[F8.4]
       display                                z=$z[F8.4] to $zstop[F8.4]
       display   number of primary atoms in box=$natom_prim
       display   number of primary+symmetry atoms in box=$natom_box
       display   number of primary+symmetry atoms in box+cushion=$natom_box_cushion
       if ( $total_test > 0 ) then 
         display   R-free=$map_free_r, R=$map_r
       else
         display   R=$map_r
       end if

       if ( $natom_prim > 0 ) then
         display   percentage of molecule omitted (with cushion)=$percent %
         display   omitted atoms (including symmetry equivalents):
         for $atom in id ( sfbox $x1 $x2 $y1 $y2 $z1 $z2 ) loop atom
           show (segid) (id $atom)
           evaluate ($segid=$result)
           show (resid) (id $atom)
           evaluate ($resid=$result)
           show (resname) (id $atom)
           evaluate ($resname=$result)
           show (name) (id $atom)
           evaluate ($name=$result)
           display     $resname[a4] $resid[a4] $name[a4] $segid[a4]
         end loop atom
       end if

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

       set display=OUTPUT end

       eval ($total=$total+1)
 
       eval ($z=$z+$cube_length/&c)
     end loop zz

     eval ($y=$y+$cube_length/&b)
   end loop yy

   eval ($x=$x+$cube_length/&a)
 end loop xx

 evaluate ($display=&output_root + ".list")
 set display=$display end
 
 display
 display -----------------------------------------
 display    total number of map boxes=      $total
 display    total number of annealed boxes= $anneal_total
 display -----------------------------------------
 display

 set display=OUTPUT end
  
 xray
   undeclare name=dtarg domain=reciprocal end
   undeclare name=fmap domain=reciprocal end
 end

 set messages=normal end
 set echo=off end
 
 if (&map_scale=true) then    
   xray                                        
     show rms (real(map_total)) ( all )
     if ( $result > 0 ) then
        do (map_total=map_total/$result) ( all )
     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=ft(map_total)) (all)
     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_coeff
     end
     undeclare name=map_coeff domain=reciprocal end
   end
 end if
 
 {- retrieve coordinates -}
 do ( x=refx ) ( all )
 do ( y=refy ) ( all )
 do ( z=refz ) ( all )

 {- MODIFIED -}
 {--------------------------------------------------------------}
 {- begin real space local correlation coefficient calculation -}
 if ( &real_r = true ) then
    
   xray
     associate fcalc ( &atom_select )
      
     predict
        mode=reciprocal
        to=fcalc
        selection=( all_active=1 )
        atomselection=( &atom_select ) 
     end
     
     declare name=mod_map domain=real end

     do (mod_map=ft(fcalc+fbulk)) (&low_res >= d >= &high_res)     
     
     if (&map_scale=true) then    
       show rms (real(mod_map)) ( all )
       do (mod_map=mod_map/$result) ( all )
     end if

     declare name=corr_map domain=real end
     declare name=map11    domain=real end
     declare name=map12    domain=real end
     declare name=map22    domain=real end
     declare name=map11_ave domain=real end {- ADDED -}
     declare name=map22_ave domain=real end {- ADDED -}
     declare name=prop     domain=real end 
     declare name=dist     domain=real end
     declare name=mask     domain=real end
   end

   do (store9=0) (all)
   evaluate ($counter=1)
   for $id in id ( tag and byresidue ( &atom_real and &atom_select)  ) loop main
     do ( store9=$counter ) ( byres( id $id) )
     evaluate ($counter=$counter+1)
   end loop main

   xray
     mask
       mode=vdw
       solrad=0.1
       shrink=0.1
       nshell=1
       to=mask
       selection=( &atom_real and &atom_select ) 
     end

     do (prop=0) (all)
     proximity
       from=store9
       distance_map=dist
       property_map=prop
       cutoff=3.0
       selection=( &atom_real and &atom_select )
     end

{- ADDED / MODIFIED -}
     do (map11_ave=gave(real(map_total), real(prop))) 
                              ( real(prop) > 0 and real(mask) <= 0 )
     do (map22_ave=gave(real(mod_map), real(prop))) 
                              ( real(prop) > 0 and real(mask) <= 0 )            
                                          
     do (map11=gave((real(map_total)-map11_ave )* (real(map_total)-map11_ave), real(prop)) ) 
                              ( real(prop) > 0 and real(mask) <= 0 )
     do (map12=gave((real(map_total)-map11_ave )* (real(mod_map)-map22_ave ), real(prop)) ) 
                              ( real(prop) > 0 and real(mask) <= 0 )
     do (map22=gave((real(mod_map)-map22_ave )* (real(mod_map)-map22_ave), real(prop)) ) 
                              ( real(prop) > 0 and real(mask) <= 0 )
     do (corr_map=real(map12)/ max( 0.0001, sqrt ( real(map11) * real(map22) ) ))
                              ( real(prop) > 0 and real(mask) <= 0 )
   end

   evaluate ($display=&output_root + "_r.list")
   set display=$display end

   display    #   resid   segid   CC   R-value

   evaluate ($realr_ave=0)
   evaluate ($realr_ave_2=0) 
   evaluate ($realr_num=0)
   
   do (store9=0) ( all )
   evaluate ($counter=1)
   for $id in id ( tag and byresidue (&atom_real and &atom_select) ) loop real
     xray
       show ave (real(corr_map)) (real(prop)=$counter and real(mask) <= 0)
     end
     evaluate ($corr=$result)
     evaluate ($realr=1-$corr)
     evaluate ($realr_ave=$realr_ave + $realr)
     evaluate ($realr_ave_2=$realr_ave_2 + $realr^2)
     evaluate ($realr_num=$realr_num+1)
     
     do (store9=$realr) ( byresidue ( id $id ) and &atom_real and &atom_select ) 
     show (resid) (id $id)
     evaluate ($resid=$result)
     show (segid) (id $id)
     evaluate ($segid=$result)
     display $counter[i6]   $resid[a4]   $segid[a4] $corr[f6.3] $realr[f6.3]
     evaluate ($counter=$counter+1)
   end loop real
   
   if ($realr_num>0) then
      evaluate ($realr_ave=$realr_ave/$realr_num) 
      evaluate ($realr_ave_2=$realr_ave_2/$realr_num)
      evaluate ($realr_sigma=sqrt ( $realr_ave_2 - $realr_ave^2) ) 
      display  average R-value = $realr_ave  sigma = $realr_sigma
   end if
   

   evaluate ($filename=&output_root + "_r.pdb")
   
   do ( b=recall9) ( all )
   set remarks=reset end
   remarks Real space R value ( = 1 - local correlation coefficient ) is in B-factor array
   write coor 
      output=$filename 
      selection=( byresidue ( &atom_real and &atom_select ) ) 
      format=PDBO
   end

   xray
     undeclare name=mod_map  domain=real end
     undeclare name=corr_map domain=real end
     undeclare name=map11    domain=real end
     undeclare name=map12    domain=real end
     undeclare name=map22    domain=real end
     undeclare name=map11_ave    domain=real end {- ADDED -}
     undeclare name=map22_ave    domain=real end {- ADDED -}
     undeclare name=prop     domain=real end
     undeclare name=dist     domain=real end
     undeclare name=mask     domain=real end
   end
 end if
 {- end real space local correlation coefficient calculation -}
 {------------------------------------------------------------}
 {- END MODIFICATION -}
 
 set remarks=reset end
 set remarks=accumulate end

 evaluate ($remark="")

 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=$remark + encode(&u) + "fo-" + encode(&v) + "fc ")
 if ( &map_type = "sigmaa" ) then
   evaluate ($remark=$remark + "sigmaa")
   if ( $test_exist = true ) then
     evaluate ($remark=$remark + "cv")
   end if
 end if
 evaluate ($remark=$remark + " complete omit map")

 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 or |F|=0): $unobserved[I6] ( $per_unobs[f5.1] % )
 remark number of reflections rejected:                       $rejected[I6] ( $per_reject[f5.1] % )
 remark total number of reflections used:                     $total_used[I6] ( $per_used[f5.1] % )
 remark number of reflections in working set:                 $total_work[I6] ( $per_work[f5.1] % )
 remark number of reflections in test set:                    $total_test[I6] ( $per_test[f5.1] % )

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

 ! 
 ! Begin modification (6/28/06)
 if ( &bscale = "anisotropic" ) then
   REMARK Anisotropic B-factor tensor Ucart of atomic model without isotropic component :
   REMARK   B11=$Baniso_11[f8.3] B22=$Baniso_22[f8.3] B33=$Baniso_33[f8.3]
   REMARK   B12=$Baniso_12[f8.3] B13=$Baniso_13[f8.3] B23=$Baniso_23[f8.3]
   REMARK Isotropic component added to coordinate array B: $Biso_model[f8.3]
 elseif ( &bscale = "isotropic" ) then
   REMARK B-factor applied to coordinate array B: $Biso_model[f8.3]
 else
   REMARK initial B-factor correction: none
 end if
 ! End modification
 !
 
 {- BEGIN MODIFICATION -}
 if (&bsharp # 0) then
   remark B-factor sharpening applied to map: exp( Bsharp * S^2/4 ) where Bsharp = &bsharp
 end if
 {- END MODIFICATION -}

 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_total
     output=$output_map
     extent=molecule
     cushion=3.0
     selection=( &atom_map and &atom_select )
   end
 end

 if ( &peak_search = true ) then
   evaluate ($filename=&output_root + "_positive.peaks")
   xray
     peakpik
       from=map_total
       mpeak=&peak_num
       selection=( all )
       atom=true
       proximity= (&atom_map and &atom_select )
     end
   end

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

   delete sele=(segid=PEAK) end

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

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

 {- MODIFICATION 04/03/06 -}
 xray
    undeclare name=map_total     domain=real end
    undeclare name=omit_mask     domain=real end 
    undeclare name=map           domain=real end  
 end
 {- END MODIFICATION -} 
 
 stop