{+ file: fo-fo_map.inp +}
{+ directory: xtal_refine +}
{+ description: Make an electron density map using phase information 
                from a model and 2 experimental amplitude data sets +}
{+ comment:
            choice of coefficients: 
              (u |Fo2| e^(i phi_calc)) - (v k|Fo1| e^(i phi_calc))
              where phi_calc is from the calculated structure factors.
            Fo2 will be scaled to Fo1 +}
{+ authors: Axel T. Brunger and Paul D. Adams +}
{+ copyright: Yale University +}

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

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

{* molecular topology file *}
{===>} structure_infile="pen.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="pen.pdb";

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

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

{* unit cell parameters in Angstroms and degrees *}
{+ table: rows=1 "cell" cols=6 "a" "b" "c" "alpha" "beta" "gamma" +}
{===>} a=97.37;
{===>} b=46.64;
{===>} c=65.47;
{===>} alpha=90;
{===>} beta=115.4;
{===>} 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="penicillopepsin.hkl";
{===>} reflection_infile_2="";
{===>} reflection_infile_3="";

{* reciprocal space array containing first set of observed
   amplitudes (Fo1): required *}
{* selection of reflections and data corrections will be made
   with reference to this data set *}
{===>} obs_f1="fobs";

{* reciprocal space array containing sigma values for first
   set of amplitudes (Fo1): required *}
{===>} obs_sigf1="sigma";

{* reciprocal space array containing second set of observed
   amplitudes (Fo2): required *}
{===>} obs_f2="ppu2";

{* reciprocal space array containing sigma values for second
   set of amplitudes (Fo2): required *}
{===>} obs_sigf2="ppu2_s";

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

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

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

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

{* rms outlier cutoff: applied to amplitudes *}
{* 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 map calculation *}
{===>} atom_select=(known and not hydrogen);

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

{* maps are calculated u*Fo2 - v*Fo1 *}
{* eg. 2fo2-fo1 map -> u=2 and v=1 or
        fo2-fo1 map -> u=1 and v=1 *}

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

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

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

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

{=========================== 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
         Fourier coefficients will be in: <output_root>.coeff +}
{===>} output_root="fo-fo_map";

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

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

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

{* write a reflection file with the Fourier coefficients of the map *}
{+ list: arrays written:
         map_phase: phases calculated from atomic model
         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

 xray

   set echo=off end

   if ( &BLANK%obs_f1 = true ) then
     display
     display  *********************************************************
     display  Error: required observed amplitude array is not specified
     display  *********************************************************
     display
     abort
   else
     query name=&STRIP%obs_f1 domain=reciprocal end
     if ( $object_exist = false ) then
       display
       display  ***************************************************************
       display  Error: required observed amplitude array &obs_f1 does not exist
       display  ***************************************************************
       display
       abort
     end if
   end if

   if ( &BLANK%obs_sigf1 = true ) then
     display
     display  *****************************************************
     display  Error: required observed sigma array is not specified
     display  *****************************************************
     display
     abort
   else
     query name=&STRIP%obs_sigf1 domain=reciprocal end
     if ( $object_exist = false ) then
       display
       display  **************************************************************
       display  Error: required observed sigma array &obs_sigf1 does not exist
       display  **************************************************************
       display
       abort
     end if
     if ( $object_type # "REAL" ) then
       display
       display  ***********************************************************************
       display  Error: required observed sigma array &obs_sigf1 has the wrong data type
       display  ***********************************************************************
       display
       abort
     end if
   end if

   if ( &BLANK%obs_f2 = true ) then
     display
     display  *********************************************************
     display  Error: required observed amplitude array is not specified
     display  *********************************************************
     display
     abort
   else
     query name=&STRIP%obs_f2 domain=reciprocal end
     if ( $object_exist = false ) then
       display
       display  ***************************************************************
       display  Error: required observed amplitude array &obs_f2 does not exist
       display  ***************************************************************
       display
       abort
     end if
   end if

   if ( &BLANK%obs_sigf2 = true ) then
     display
     display  *****************************************************
     display  Error: required observed sigma array is not specified
     display  *****************************************************
     display
     abort
   else
     query name=&STRIP%obs_sigf2 domain=reciprocal end
     if ( $object_exist = false ) then
       display
       display  **************************************************************
       display  Error: required observed sigma array &obs_sigf2 does not exist
       display  **************************************************************
       display
       abort
     end if
     if ( $object_type # "REAL" ) then
       display
       display  ***********************************************************************
       display  Error: required observed sigma array &obs_sigf2 has the wrong data type
       display  ***********************************************************************
       display
       abort
     end if
   end if

   if ( &BLANK%test_set = false ) then
     query name=&STRIP%test_set domain=reciprocal end
     if ( $object_exist = false ) then
       display
       display  **********************************************
       display  Error: test set array &test_set does not exist
       display  **********************************************
       display
       abort
     end if
     {- test set array can be integer or real -}
     if ( $object_type = "COMPLEX" ) then
       display
       display  *******************************************************
       display  Error: test set array &test_set has the wrong data type
       display  *******************************************************
       display
       abort
     end if
   end if

   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 )

   evaluate ($obs_lower_limit=0)

   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%obs_f1) > $obs_lower_limit ) and
                       ( amplitude(&STRIP%obs_f2) > $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%obs_f1)) ( ref_active=1 )
   evaluate ($obs_high_f1=$result*&obs_rms)
   show min (amplitude(&STRIP%obs_f1)) ( ref_active=1 )
   evaluate ($obs_low_f1=$result)

   show rms (amplitude(&STRIP%obs_f2)) ( ref_active=1 )
   evaluate ($obs_high_f2=$result*&obs_rms)
   show min (amplitude(&STRIP%obs_f2)) ( ref_active=1 )
   evaluate ($obs_low_f2=$result)

   do (ref_active=0) ( all )
   do (ref_active=1)
                  ( ( amplitude(&STRIP%obs_f1) >= &sigma_cut*&STRIP%obs_sigf1 ) and
                    ( amplitude(&STRIP%obs_f2) >= &sigma_cut*&STRIP%obs_sigf2 ) and
                    ( &STRIP%obs_sigf1 # 0 ) and
                    ( &STRIP%obs_sigf2 # 0 ) and
                    ( $obs_low_f1 <= amplitude(&STRIP%obs_f1) <= $obs_high_f1 ) and
                    ( $obs_low_f2 <= amplitude(&STRIP%obs_f2) <= $obs_high_f2 ) and
                    ( &low_res >= d >= &high_res ) )

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

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

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

   associate fcalc ( &atom_select )

   tselection=( ref_active=1 )

   cvselection=( tst_active=1 )

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

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

   tolerance=0.0 lookup=false
   
 end                  

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

 xray
   declare name=total domain=reciprocal type=complex end
   declare name=fmap domain=reciprocal type=complex end
 end

 xray
    predict
      mode=reciprocal
      to=fcalc
      selection=(ref_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_f1;
                             obs_sigf=&STRIP%obs_sigf1;
                             obs_i="";
                             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_f1;
                                 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
   @@CNS_XTALMODULE:calculate_r (fobs=&STRIP%obs_f1;
                                 fcalc=fcalc;
                                 fpart=fbulk;
                                 sel=(ref_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

 xray

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

   if ( &bscale = "no" ) then
     @CNS_XTALMODULE:scalef (kscale="yes";
                             bscale="no";
                             fref=&STRIP%obs_f1;
                             f=&STRIP%obs_f2;
                             s=&STRIP%obs_sigf2;
                             sel=(ref_active=1);)
   elseif ( &bscale = "isotropic" ) then
     @CNS_XTALMODULE:scalef (kscale="yes";
                             bscale="isotropic";
                             fref=&STRIP%obs_f1;
                             f=&STRIP%obs_f2;
                             s=&STRIP%obs_sigf2;
                             sel=(ref_active=1);)
   elseif ( &bscale = "anisotropic" ) then
     @CNS_XTALMODULE:scalef (kscale="yes";
                             bscale="anisotropic";
                             fref=&STRIP%obs_f1;
                             f=&STRIP%obs_f2;
                             s=&STRIP%obs_sigf2;
                             sel=(ref_active=1);)
   end if

   do (map_scale=1.0) (all)

   do (map_fom=1.0) (all)

 end

 xray
   if ( &u = &v ) then
     do (fmap= 2 ((&u   map_fom combine(amplitude(&STRIP%obs_f2),map_phase)) -
                  (&v map_scale combine(amplitude(&STRIP%obs_f1),map_phase))))
        (ref_active=1 and acentric)
     do (fmap=    (&u   map_fom combine(amplitude(&STRIP%obs_f2),map_phase)) -
                  (&v map_scale combine(amplitude(&STRIP%obs_f1),map_phase)))
        (ref_active=1 and centric)
   else
     do (fmap=(&u             map_fom combine(amplitude(&STRIP%obs_f2),map_phase)) - 
              (&v           map_scale combine(amplitude(&STRIP%obs_f1),map_phase)))
        (ref_active=1 and acentric)
     do (fmap=(max((&u-1),0)   map_fom combine(amplitude(&STRIP%obs_f2),map_phase)) -
              (max((&v-1),0) map_scale combine(amplitude(&STRIP%obs_f1),map_phase)))
        (ref_active=1 and centric)
   end if
 end 


 xray
 
   {- BEGIN MODIFICATION -}
   if (&bsharp # 0) then 
      {- b-factor sharpening -} 
      evaluate ($bsharp= (-1) * &bsharp) 
      do (fmap=exp( $bsharp * s^2/4) * fmap) ( all )
   end if
   {- END MODIFICATION -}
           
   declare name=map domain=real end
   do (map=ft(fmap)) (ref_active=1) 
 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
       sele=(ref_active=1)
       map_phase map_coeff
     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=total domain=reciprocal end
   undeclare name=fmap  domain=reciprocal end
 end

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

 set remarks=reset end
 set remarks=accumulate end

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

 evaluate ($remark="")
 evaluate ($remark=$remark + "(" + encode(&u) +  " |Fo2| - " + 
                                   encode(&v) + " k|Fo1|)e^(i phi_calc)")
 evaluate ($remark=$remark + " map")

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

 if ( &obs_type = "intensity" ) then
   remark reflections with Iobs/sigma_I < &sigma_cut rejected
   remark reflections with Iobs > &obs_rms * rms(Iobs) rejected
 else
   remark reflections with |Fobs|/sigma_F < &sigma_cut rejected
   remark reflections with |Fobs| > &obs_rms * rms(Fobs) rejected
 end if
 xray anomalous=? end
 if ( $result = true ) then
   remark anomalous diffraction data was input
 end if
 
 {- MODIFIED 2/15/06 -}
 remark fft gridding factor = $fft_grid, B factor offset = $fft_b_add A^2, Elimit = $fft_elim
 {- END MODIFICATION -}
 
 remark theoretical total number of refl. in resol. range:    $total_theor[I6] ( 100.0 % )
 remark number of unobserved reflections (no entry 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 -}

 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