{+ file: model_map_twin.inp +}
{+ directory: xtal_twin +}
{+ description: Make an electron density map using phase information 
                from a model, for data with hemihedral twinning +}
{+ 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))
               d(target)/dFc
              where Fo is the detwinned observed structure factor,
                    Fc is the calculated structure factor and
                    m and D are derived from sigmaa
            Realspace R-value calculation optional. +}
{+ authors: Axel T. Brunger and Paul D. Adams +}
{+ copyright: Yale University +}

{+ reference: T.O. Yeates, Detecting and Overcoming Crystal Twinning,
              Methods in Enzymology 276, 344-358 (1997) +}
{+ 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) +}
{+ 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 -}

{========================== IMPORTANT NOTE ===========================}

{* This input file requires that the test set has been generated with
   the CNS task file make_cv_twin.inp. This ensures that twin related
   reflections are part of the same set. *}

{- begin block parameter definition -} define(

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

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

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

{* coordinate files *}
{===>} coordinate_infile="porin.pdb";

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

{* space group *}
{* use International Table conventions with subscripts substituted
   by parenthesis *}
{===>} sg="R3";

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

{* 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="porin.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 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 and phase combined or observed maps *}
{+ 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.25;

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

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

{======================= twinning parameters =========================}

{* twinning operation *}
{===>} twin_oper="h,-h-k,-l";

{* twinning fraction *}
{===>} twin_frac=0.304;

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

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

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

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

{* specify u *}
{===>} u=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
         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" "gradient" +}
{===>} map_type="sigmaa";

{* refinement target *}
{+ list: twin_lsq: least squares residual for hemihedral twinning +}
{+ choice: "twin_lsq" +}
{===>} reftarget="twin_lsq";

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

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

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

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

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

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

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

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

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

{* calculate realspace R-values *}
{+ 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
         Fourier coefficients will be in: <output_root>.coeff +}
{===>} output_root="model_map_twin";

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

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

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

{* write a reflection file with the Fourier coefficients of the map *}
{+ list: arrays written:
         map_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

 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

 set echo=off end
 if ( &twin_frac > 0.5 ) then
   display Error: twinning fraction must be less than or equal 0.5
   abort
 end if
 set echo=on end

 {- 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=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) ( ( amplitude($STRIP%reject_obs) > $obs_lower_limit ) 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
                    ( $obs_low <= amplitude(remap[&STRIP%twin_oper]($STRIP%reject_obs)) <= $obs_high ) and
                    ( &low_res >= d >= &high_res ) )

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

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

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

   associate fcalc ( &atom_select )

   tselection=( ref_active=1 )

   cvselection=( tst_active=1 )

   method=FFT          
{- MODIFIED 8/01/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
   query name=&STRIP%obs_f domain=reciprocal end
   declare name=total domain=reciprocal type=complex end
   declare name=fmap domain=reciprocal type=complex end
 end

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

   xray
      predict
        mode=reciprocal
        to=fcalc
        selection=( all )
        atomselection=( &atom_select ) 
      end
   end

   {- BEGIN MODIFICATION 8/01/06 -}
   @CNS_XTALMODULE:scale_and_solvent_grid_search (
                             bscale=&bscale;
                             sel=( ref_active=1 );
                             sel_test=( tst_active=1 );
                             atom_select=( &atom_select );
                             bulk_sol=&bulk_sol;
                             bulk_mask=&bulk_mask_infile;
                             bulk_atoms=( &atom_select );
                             
                             sol_auto=&sol_auto;
                             sol_k=&sol_k;
                             sol_b=&sol_b;
                             sol_rad=&sol_rad;
                             sol_shrink=&sol_shrink;
  
                             fcalc=fcalc;
                             obs_f=&STRIP%obs_f;
                             obs_sigf=&STRIP%obs_sigf;
                             obs_i=$STRIP%obs_i;
                             obs_sigi=$STRIP%obs_sigi;                             
                             fpart=fbulk;

                             Baniso_11=$Baniso_11;
                             Baniso_22=$Baniso_22;
                             Baniso_33=$Baniso_33;
                             Baniso_12=$Baniso_12;
                             Baniso_13=$Baniso_13;
                             Baniso_23=$Baniso_23;
                             Biso=$Biso_model;
                             
                             sol_k_best=$sol_k_ref;
                             sol_b_best=$sol_b_ref;
			     solrad_best=$solrad_best;
			     shrink_best=$shrink_best;
                             
                             b=b;

                             low_b_flag=$low_b_flag;
                            
                             sol_output=&sol_output;
                             
                             )
 
   {- check the gridding again since the minimum B-factor may have changed -}
   show min ( b ) ( &atom_select )
   evaluate ($b_min=$result)
   @@CNS_XTALMODULE:fft_parameter_check ( 
                             d_min=&high_res; 
                             b_min=$b_min;
                             grid=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_twin (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;
                                               twin_oper=&STRIP%twin_oper;
                                               twin_frac=&twin_frac;
                                               pa=&STRIP%obs_pa;
                                               pb=&STRIP%obs_pb;
                                               pc=&STRIP%obs_pc;
                                               pd=&STRIP%obs_pd;
                                               phase=&STRIP%obs_phase;
                                               fom=&STRIP%obs_fom;
                                               sel=(ref_active=1);
                                               sel_test=(tst_active=1);
                                               statistics=true;)
       predict
          mode=dtarget(fcalc)
          to=dtarg
          selection=(ref_active=1)
          atomselection=( &atom_select ) 
       end

     end if

   end
 
 xray
   @@CNS_XTALMODULE:calculate_r_twin (fobs=&STRIP%obs_f;
                                      fcalc=fcalc;
                                      fpart=fbulk;
                                      twin_oper=&STRIP%twin_oper;
                                      twin_frac=&twin_frac;
                                      sel=(ref_active=1);
                                      sel_test=(tst_active=1);
                                      print=true;
                                      output=OUTPUT;
                                      r=$map_r;
                                      test_r=$map_free_r;)
 end

 {- detwin the observed data -}
 xray
   if ( &twin_frac < 0.5 ) then
     do (&STRIP%obs_f=fobs_orig) (all)
     do (&STRIP%obs_sigf=sigma_orig) (all)
   end if
   query name=fo_detwin domain=reciprocal end
   if ( $object_exist = false ) then
     declare name=fo_detwin domain=reciprocal type=real end
   else
     do (fo_detwin=0) (all)
   end if
   declare name=ref_det_active domain=reciprocal type=integer end
   declare name=tst_det_active domain=reciprocal type=integer end
 end

 xray

   @@CNS_XTALMODULE:data_detwin (fobs=&STRIP%obs_f;
                                 fcalc=fcalc;
                                 fpart=fbulk;
                                 twin_oper=&STRIP%twin_oper;
                                 twin_frac=&twin_frac;
                                 sel=(ref_active=1);
                                 sel_test=(tst_active=1);
                                 fobs_detwin=fo_detwin;
                                 ref_detwin=ref_det_active;
                                 tst_detwin=tst_det_active;
                                 ref_reject=$ref_detwin_reject;
                                 tst_reject=$tst_detwin_reject;
                                 total_reject=$detwin_reject;)
 end

   @CNS_XTALMODULE:scale_and_solvent_grid_search (
                             bscale=&bscale;
                             sel=( ref_det_active=1 );
                             sel_test=( tst_det_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=fo_detwin;
                             obs_sigf=&STRIP%obs_sigf;
                             obs_i=$STRIP%obs_i;
                             obs_sigi=$STRIP%obs_sigi;                             
                             fpart=fbulk;

                             Baniso_11=$Baniso_11;
                             Baniso_22=$Baniso_22;
                             Baniso_33=$Baniso_33;
                             Baniso_12=$Baniso_12;
                             Baniso_13=$Baniso_13;
                             Baniso_23=$Baniso_23;
                             Biso=$Biso_model;
                             
                             sol_k_best=$sol_k_ref;
                             sol_b_best=$sol_b_ref;
			     solrad_best=$solrad_best;
			     shrink_best=$shrink_best;
                             
                             b=b;

                             low_b_flag=$low_b_flag;
                            
                             sol_output=&sol_output;
                             
                             )


 xray
   @@CNS_XTALMODULE:calculate_r (fobs=fo_detwin;
                                 fcalc=fcalc;
                                 fpart=fbulk;
                                 sel=(ref_det_active=1);
                                 sel_test=(tst_det_active=1);
                                 print=true;
                                 output=OUTPUT;
                                 r=$detwin_r;
                                 test_r=$detwin_test_r;)

 end

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

 if ( &map_type = "unweighted" ) then

   xray

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

     do (total=fcalc+fbulk) (all)
     multiscale
       bfmin=-40 bfmax=40
       set1=fo_detwin  k1=-1 b1=0
       set2=total            b2=0
       selection=(ref_det_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=fo_detwin;
                                   fcalc=fcalc;
                                   fpart=fbulk;
                                   test=tst_det_active;
                                   sel=(ref_det_active=1);
                                   sel_test=(tst_det_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

 end if

 if ( &map_type = "gradient" ) then
 
   xray
     {- take the negative of the gradient so the map is the same sign
        as a standard difference map -}
     do (fmap=-dtarg) (ref_active=1)
     do (map_fom=0.0) (ref_active=1)
     do (map_scale=0.0) (ref_active=1)
     do (map_phase=0.0) (ref_active=1)
   end

 else

   xray
     if ( &u = &v ) then
       do (fmap= 2 ((&u   map_fom combine(amplitude(fo_detwin),map_phase)) -
                    (&v map_scale        (fcalc+fbulk))))
          (ref_det_active=1 and acentric)
       do (fmap=    (&u   map_fom combine(amplitude(fo_detwin),map_phase)) -
                    (&v map_scale        (fcalc+fbulk)))
          (ref_det_active=1 and centric)
     else
       do (fmap=(&u             map_fom combine(amplitude(fo_detwin),map_phase)) - 
                (&v           map_scale        (fcalc+fbulk)))
          (ref_det_active=1 and acentric)
       do (fmap=(max((&u-1),0)   map_fom combine(amplitude(fo_detwin),map_phase)) -
                (max((&v-1),0) map_scale        (fcalc+fbulk)))
          (ref_det_active=1 and centric)
       if ( &fill_in = true ) then
         do (fmap=(&u-&v) map_scale (fcalc+fbulk)) 
            (&low_res >= d >= &high_res and ref_det_active # 1)
       end if
     end if
   end 

 end if

 xray
   declare name=map domain=real end
   if ( &map_type = "gradient" ) then
     do (map=ft(fmap)) (ref_active=1)
   elseif ( &fill_in = true ) then
     do (map=ft(fmap)) (&low_res >= d >= &high_res)
   else
     do (map=ft(fmap)) (ref_det_active=1) 
   end if
 end

 if ( &write_coeff = true ) then
   evaluate ($coeff_out=&output_root + ".coeff")
   xray
     declare name=map_coeff domain=reciprocal type=complex end
     do (map_coeff=fmap) (all)
     write reflection
       output=$coeff_out
       if ( &map_type = "gradient" ) then
         sele=(ref_active=1)
       elseif ( &fill_in = true ) then
         sele=(&low_res >= d >= &high_res)
       else
         sele=(ref_det_active=1)
       end if
       map_coeff {- MODIFIED 8/1/06 : map_coeff only -}
     end
     undeclare name=map_coeff domain=reciprocal end
   end
 end if

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

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

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


 {- MODIFIED -}
 {--------------------------------------------------------------}
 {- begin real space local correlation coefficient calculation -}
 if ( &real_r = true ) then
    
   xray
     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

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

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

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

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

 remark twinning operator= &STRIP%twin_oper  twinning fraction= &twin_frac

 remark reflections rejected (I_detwin <= 0): $detwin_reject
 remark                          working set: $ref_detwin_reject
 if ( $tst_detwin_reject > 0 ) then
   remark                             test set: $tst_detwin_reject
 end if

 if ( $total_test > 0 ) then
   remark detwinned r= $detwin_r[f6.4] detwinned free_r= $detwin_test_r[f6.4]
 else
   remark detwinned r= $detwin_r[f6.4]
 end if

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

 if ( &obs_type = "intensity" ) then
   remark reflections with Iobs/sigma_I < &sigma_cut rejected
   remark reflections with Iobs > &obs_rms * rms(Iobs) rejected
   remark reflections with Iobs[&STRIP%twin_oper] = 0 rejected
 else
   remark reflections with |Fobs|/sigma_F < &sigma_cut rejected
   remark reflections with |Fobs| > &obs_rms * rms(Fobs) rejected
   remark reflections with |Fobs|[&STRIP%twin_oper] = 0 rejected
 end if
 xray anomalous=? end
 if ( $result = true ) then
   remark anomalous diffraction data was input
 end if
 
 {- MODIFIED 8/01/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] % )

 ! 
 ! Begin modification (8/0/106)

 if ( &bscale = "no" ) then
   remark initial B-factor correction: none
 end if

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

 if ( &write_map = true ) then

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

   if ( &map_mode = "asymmetric" ) then
     evaluate ($map_mode_string=ASYM)
   elseif ( &map_mode = "unit" ) then
     evaluate ($map_mode_string=UNIT)
   elseif ( &map_mode = "box" ) then
     evaluate ($map_mode_string=BOX)
   elseif ( &map_mode = "fract" ) then
     evaluate ($map_mode_string=FRAC)
   else
     evaluate ($map_mode_string=MOLE)
   end if

   xray
     write map
       if ( &map_format = "ezd" ) then
         type=ezd
       else
         type=cns
       end if
       automatic=false
       from=map
       output=$filename
       cushion=&map_cushion
       selection=&atom_map
       extend=$map_mode_string
       if ( &map_mode = "box" ) then
         xmin=&xmin xmax=&xmax
         ymin=&ymin ymax=&ymax
         zmin=&zmin zmax=&zmax
       end if
       if ( &map_mode = "fract" ) then
         xmin=&xmin xmax=&xmax
         ymin=&ymin ymax=&ymax
         zmin=&zmin zmax=&zmax
       end if
     end
   end
 end if

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

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

   delete sele=(segid=PEAK) end

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

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

 stop