{+ file: water_delete.inp +}
{+ directory: xtal_refine +}
{+ description: delete water molecules by simple analysis of 
                electron density map based on real space local
                correlation coefficient and peak height +}
{+ comments:
             Water molecules can be deleted on the basis of realspace
             correlation coefficient and/or maximum density level at
             the water position.
             The water molecules selected for analysis are excluded from 
             the map calculation.
             Choice of map coefficients:
               (u m|Fo| e^(i phi_calc)) - (v D|Fc| e^(i phi_calc))
               (u  |Fo| e^(i phi_calc)) - (v k|Fc| e^(i phi_calc))
               (u m|Fo| e^(i phi_comb)) - (v D|Fc| e^(i phi_calc))
               (u m_obs|Fo| e^(i phi_obs )) - (v k|Fc| e^(i phi_calc))
                d(target)/dFc
               where Fc is the calculated structure factor
               where m and D are derived from sigmaa +}
{+ authors: Axel T. Brunger and Paul D. Adams +}
{+ copyright: Yale University +}

{+ reference: R.J. Read, Improved Fourier coefficients for maps using
              phases from partial structures with errors.
              Acta Cryst. A42, 140-149 (1986) +}
{+ reference: G.J. Kleywegt and A.T. Brunger, Checking your imagination:
              Applications of the free R value, Structure 4,
              897-904 (1996) +}

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

{- begin block parameter definition -} define(

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

{* coordinate file *}
{===>} coordinate_infile="water_pick.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="water_pick.mtf";

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


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

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

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

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

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

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

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

{* reciprocal space array containing test set for cross-validation:
   required for calculation of cross-validated sigmaA values *}
{===>} test_set="test";

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

{* automatic bulk solvent parameter optimization for e-density level sol_k (e/A^3) and B-factor sol_b (A^2) *}
{+ choice: true false +}
{===>} sol_auto=true;

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

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

{* solvent mask parameters *}
{+ table: rows=1 "bulk solvent" cols=2 "probe radius (A) (usually set to 1)" "shrink radius (A) (usually set to 1)" +}
{===>} sol_rad=1.0;
{===>} sol_shrink=1.0;

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

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

{* select water molecules to be analysed for possible deletion *}
{* these atoms will not be included in the map calculation *}
{===>} water_select=(segid S);

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

{* minimum realspace correlation at a water position *}
{* waters with a realspace correlation less than this will be deleted *}
{===>} min_correl=0.50;

{* minimum peak height in sigma units at a water position *}
{* waters with a peak height less than this will be deleted *}
{===>} min_peak=3.0;

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

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

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

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

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

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

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

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

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

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

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

{* memory allocation for FFT calculation *}
{* this will be determined automatically if a negative value is given
   otherwise the specified number of words will be allocated *}
{===>} fft_memory=-1;

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

{* root name for output files *}
{+ list:
         listing file will be: <output_root>.list
         output molecular structure file will be: <output_root>.mtf
         output coordinate file will be: <output_root>.pdb
         water mask will be: <output_root>_water.mask +}
{===>} output_root="water_delete";

{===========================================================================}
{         things below this line do not normally need to be changed         }
{===========================================================================}

 ) {- end block parameter definition -}

 checkversion 1.3

 evaluate ($log_level=quiet)

 if ( $log_level = verbose ) then
   set message=normal echo=on end
 else
   set message=off echo=off end
 end if

 if ( &BLANK%structure_infile = true ) then
 
    {- read topology files -}
    topology
     evaluate ($counter=1)
     evaluate ($done=false)
     while ( $done = false ) loop read
      if ( &exist_topology_infile_$counter = true ) then
         if ( &BLANK%topology_infile_$counter = false ) then
            @@&topology_infile_$counter
         end if
      else
        evaluate ($done=true)
      end if
      evaluate ($counter=$counter+1)
     end loop read
    end
    
    @CNS_XTALMODULE:mtfautogenerate (
                                  coordinate_infile=&coordinate_infile;
                                  convert=true;
                                  separate=true;
                                  atom_delete=(not known);
                                  hydrogen_flag=true;
                                  break_cutoff=2.5;
                                  disulphide_dist=3.0;
                                  carbo_dist=2.5;
                                  patch_infile=&patch_infile;
                                  O5_becomes="O";
                                 )

 else

   structure @&structure_infile end
   coordinates @&coordinate_infile

 end if

 {- read parameter files -}
 parameter
  evaluate ($counter=1)
  evaluate ($done=false)
  while ( $done = false ) loop read
   if ( &exist_parameter_infile_$counter = true ) then
      if ( &BLANK%parameter_infile_$counter = false ) then
         @@&parameter_infile_$counter
      end if
   else
     evaluate ($done=true)
   end if
   evaluate ($counter=$counter+1)
  end loop read
 end

 set message=normal echo=on end

 xray

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

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

   @CNS_XRAYLIB:scatter.lib

   binresolution &low_res &high_res
   mapresolution &high_res

   generate &low_res &high_res

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

 end

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

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

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

   if ( &obs_type = "intensity" ) then
     if ( &BLANK%obs_i = true ) then
       display  Error: observed intensity array is undefined
       display         aborting script
       abort
     end if
     evaluate ($reject_obs=&obs_i)
     evaluate ($reject_sig=&obs_sigi)
   else
     evaluate ($reject_obs=&obs_f)
     evaluate ($reject_sig=&obs_sigf)
   end if

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

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

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

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

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

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

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

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

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

   associate fcalc ( &atom_select and not ( &water_select ) )

   tselection=( ref_active=1 )

   cvselection=( tst_active=1 )

   method=FFT

 {- MODIFIED 2/15/06 -}
 end
 

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

   tolerance=0.0 lookup=false

 end

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

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

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

 xray
   declare name=dtarg domain=reciprocal type=complex end
   declare name=total domain=reciprocal type=complex end
   declare name=fmap domain=reciprocal type=complex end
   predict
     mode=reciprocal
     to=fcalc
     selection=( &low_res >= d >= &high_res )
     atomselection=( &atom_select and not ( &water_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 -}

 xray
   if ( &map_type = "gradient" ) then
     @@CNS_XTALMODULE:refinementtarget (target=&reftarget;
                                       sig_sigacv=0.07;
                                       mbins=&target_bins;
                                       fobs=&STRIP%obs_f;
                                       sigma=&STRIP%obs_sigf;
                                       weight=$STRIP%obs_w;
                                       iobs=$STRIP%obs_i;
                                       sigi=$STRIP%obs_sigi;
                                       test=tst_active;
                                       fcalc=fcalc;
                                       fpart=fbulk;
                                       pa=&STRIP%obs_pa;
                                       pb=&STRIP%obs_pb;
                                       pc=&STRIP%obs_pc;
                                       pd=&STRIP%obs_pd;
                                       phase=&STRIP%obs_phase;
                                       fom=&STRIP%obs_fom;
                                       sel=(ref_active=1);
                                       sel_test=(tst_active=1);
                                       statistics=true;)
     predict
        mode=dtarget(fcalc)
        to=dtarg
        selection=( ref_active=1 )
        atomselection=( &atom_select and not ( &water_select ) ) 
     end
   end if
 end

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

 if ( &map_type = "unweighted" ) then

   xray

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

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

     do (map_fom=1.0) (all)

   end

 elseif ( &map_type = "sigmaa" ) then

   xray

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

     declare name=m          domain=reciprocal type=complex end
     declare name=mod_fom    domain=reciprocal type=real end
     declare name=mod_x      domain=reciprocal type=real end
     declare name=mod_pa     domain=reciprocal type=real end
     declare name=mod_pb     domain=reciprocal type=real end
     declare name=mod_pc     domain=reciprocal type=real end
     declare name=mod_pd     domain=reciprocal type=real end
     declare name=mod_dd     domain=reciprocal type=real end

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

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

     undeclare name=m          domain=reciprocal end
     undeclare name=mod_fom    domain=reciprocal end
     undeclare name=mod_x      domain=reciprocal end
     undeclare name=mod_pa     domain=reciprocal end
     undeclare name=mod_pb     domain=reciprocal end
     undeclare name=mod_pc     domain=reciprocal end
     undeclare name=mod_pd     domain=reciprocal end
     undeclare name=mod_dd     domain=reciprocal end

   end

 elseif ( &map_type = "combined" ) then

   xray
     declare name=m          domain=reciprocal type=complex end
     declare name=mod_fom    domain=reciprocal type=real end
     declare name=mod_x      domain=reciprocal type=real end
     declare name=mod_pa     domain=reciprocal type=real end
     declare name=mod_pb     domain=reciprocal type=real end
     declare name=mod_pc     domain=reciprocal type=real end
     declare name=mod_pd     domain=reciprocal type=real end
     declare name=mod_dd     domain=reciprocal type=real end

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

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

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

     do (map_phase=phase(m)) (all)
     do (map_fom=amplitude(m)) (all)
     do (map_scale=distribute(mod_dd)) (&low_res >= d >= &high_res)

     undeclare name=m          domain=reciprocal end
     undeclare name=mod_fom    domain=reciprocal end
     undeclare name=mod_x      domain=reciprocal end
     undeclare name=mod_pa     domain=reciprocal end
     undeclare name=mod_pb     domain=reciprocal end
     undeclare name=mod_pc     domain=reciprocal end
     undeclare name=mod_pd     domain=reciprocal end
     undeclare name=mod_dd     domain=reciprocal end
   end

 elseif ( &map_type = "observed" ) then

   xray

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

     declare name=m          domain=reciprocal type=complex end

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

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

 end if

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

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

 else

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

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

 end if

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

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

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

 xray
   declare name=mod_map domain=real end

   predict
     mode=real
     to=mod_map
     atomselection=(&atom_select and &water_select)
   end

   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 &water_select ) loop main
   do ( store9=$counter ) ( byres( id $id) )
   evaluate ($counter=$counter+1)
 end loop main

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

   do (prop=0) (all)

   proximity
     from=store9
     distance_map=dist
     property_map=prop
     cutoff=3.0
     selection=( &atom_select and &water_select )
   end

   {- BEGIN MODIFICAITON -}
   do (map11_ave=gave(real(map), 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)-map11_ave )* (real(map)-map11_ave), real(prop)) ) 
                            ( real(prop) > 0 and real(mask) <= 0 )
   do (map12=gave((real(map)-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 MODIFICATION -}

 end

 evaluate ($water_del=0)
 ident (store9) (none)

 evaluate ($display_outfile=&output_root + ".list")

 set display=$display_outfile end
 display
 display ==============================================================
 display  waters deleted if:
 display                    realspace correlation < &min_correl
 display                 or peak height < &min_peak
 display ==============================================================
 display
 set display=OUTPUT end

 evaluate ($counter=1)
 for $id in id ( tag and &water_select ) loop water
   xray
     show ave (real(corr_map)) (real(prop)=$counter and real(mask) <= 0)
   end
   evaluate ($corr=$result)
   xray
     show max (real(map)) (real(prop)=$counter and real(mask) <= 0)
   end
   evaluate ($peak=$result)
   show (resname) (id $id)
   evaluate ($resname=$result)
   show (resid) (id $id)
   evaluate ($resid=$result)
   show (segid) (id $id)
   evaluate ($segid=$result)
   evaluate ($del=false)
   if ( $corr < &min_correl ) then
     evaluate ($del=true)
   end if
   if ( $peak < &min_peak ) then
     evaluate ($del=true)
   end if
   if ( $del = true ) then
     ident (store9) (store9 or byres(id $id))
     set display=$display_outfile end
     display deleting: water $resname $resid[a4] $segid[a4]  \
correlation=$corr[f6.3] peak=$peak[f6.3]
     set display=OUTPUT end
     evaluate ($water_del=$water_del+1)
   else 
     set display=$display_outfile end
     display  keeping: water $resname $resid[a4] $segid[a4]  \
correlation=$corr[f6.3] peak=$peak[f6.3]
     set display=OUTPUT end
   end if
   evaluate ($counter=$counter+1)
 end loop water

 delete sele=(store9) end

 xray

   @CNS_XRAYLIB:scatter.lib

   predict
     mode=reciprocal
     to=fcalc
     selection=( &low_res >= d >= &high_res )
     atomselection=( &atom_select )
   end

   @@CNS_XTALMODULE:calculate_r (fobs=&STRIP%obs_f;
                                 fcalc=fcalc;
                                 fpart=fbulk;
                                 sel=(ref_active=1);
                                 sel_test=(tst_active=1);
                                 print=true;
                                 output=OUTPUT;
                                 r=$full_r;
                                 test_r=$full_test_r;)
 end

 evaluate ($mask_outfile=&output_root + "_water.mask")

 xray
   write mask
     from=mask
     output=$mask_outfile
     cushion=3.0
     extent=molecule
   end
 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

 evaluate ($coordinate_outfile=&output_root + ".pdb")

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

 evaluate ($remark="in ")

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

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

 set display=$coordinate_outfile end

 display REMARK coordinates from water deletion
 display REMARK $water_del waters deleted with
 display REMARK   real space correlation < &min_correl or peak height < &min_peak
 display REMARK $remark
 display REMARK map resolution: &low_res - &high_res A  
 if ( $total_test > 0 ) then
   display REMARK final r= $full_r[f6.4] free_r= $full_test_r[f6.4]
 else
   display REMARK final r= $full_r[f6.4]
 end if
 display REMARK sg= &STRIP%sg a= &a b= &b c= &c alpha= &alpha beta= &beta gamma= &gamma

 evaluate ($counter=1)
 evaluate ($done=false)
 while ( $done = false ) loop read
  if ( &exist_parameter_infile_$counter = true ) then
    if ( &BLANK%parameter_infile_$counter = false ) then
      display REMARK parameter file $counter  : &STRIP%parameter_infile_$counter
    end if
  else
   evaluate ($done=true)
  end if
  evaluate ($counter=$counter+1)
 end loop read
 
 if ( &BLANK%structure_infile = true ) then
   display REMARK molecular structure file: automatic

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

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

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

 else
   display REMARK molecular structure file: &STRIP%structure_infile
 end if

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

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

 if ( &BLANK%ncs_infile = false  ) then
   display REMARK ncs= &STRIP%ncs_type  ncs file= &STRIP%ncs_infile
 else
   display REMARK ncs= none
 end if
 if ( &bscale # "no" ) then
   if ( $low_b_flag = true ) then
     display REMARK warning: B-correction gave atomic B-values less than zero
     display REMARK          they have been reset to zero
   end if
 end if
 
 ! 
 ! Begin modification (6/28/06)
 if ( &bscale = "anisotropic" ) then
   display REMARK Anisotropic B-factor tensor Ucart of atomic model without isotropic component :
   display REMARK   B11=$Baniso_11[f8.3] B22=$Baniso_22[f8.3] B33=$Baniso_33[f8.3]
   display REMARK   B12=$Baniso_12[f8.3] B13=$Baniso_13[f8.3] B23=$Baniso_23[f8.3]
   display REMARK Isotropic component added to coordinate array B: $Biso_model[f8.3]
 elseif ( &bscale = "isotropic" ) then
   display REMARK B-factor applied to coordinate array B: $Biso_model[f8.3]
 else
   display REMARK initial B-factor correction: none
 end if
 ! End modification
 !
 
 {- MODIFIED 5/18/05 -}
 if ( &bulk_sol = true ) then 
   display REMARK bulk solvent: probe radius=$solrad_best, shrink value=$solrad_best
   display REMARK bulk solvent: density level= $sol_k_ref e/A^3, B-factor= $sol_b_ref A^2
 else
   display REMARK bulk solvent: false
 end if
 {- END MODIFICATION -}
 
  {- BEGIN MODIFICATION -}
 if (&bsharp # 0) then
   display REMARK B-factor sharpening applied to map: exp( Bsharp * S^2/4 ) where Bsharp = &bsharp
 end if
 {- END MODIFICATION -}

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

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

 remark

 evaluate ($structure_outfile=&output_root + ".mtf")

 write structure output=$structure_outfile end

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

 stop