{+ file: map_cover.inp +}
{+ directory: xtal_refine +}
{+ description: Make an electron density map masked by a model +}
{+ comment: 
            Warning: this task file should be used with caution,
                     it is intended for figure preparation, 
                     not for model buliding.
            The map is calculated from previously calculated 
            coefficients (see model_map.inp). +}
{+ 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="amy.mtf";

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

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

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

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

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

{* 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="amy_map.coeff";
{===>} 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";

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

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

{========================= map coefficients ==========================}

{* complex or real reciprocal space array with amplitude coefficients *}
{===>} map_amplitude="map_coeff";

{* complex reciprocal space array with phase coefficients *}
{===>} map_phase="map_coeff";

{* optional real reciprocal space array with FOM values for weighting
   of coefficients: will be ignored if blank *}
{* this is not required for coefficients derived from a previous map
   calculation (eg. from model_map.inp) *}
{===>} map_fom="";

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

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

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

{* select atoms which map will cover*}
{===>} atom_map=(resid 30:40);

{* map cover radius *}
{* points further than this distance from an atom will be set to zero *}
{===>} map_cover=2.0;

{* use map coefficients for unmeasured data *}
{* this can only be used if it was also used to create the map coefficients *}
{+ 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;

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

{* output map file *}
{===>} map_outfile="map_cover.map";

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

 xray

   if ( &BLANK%obs_f = true ) then
     display
     display  Error: required observed amplitude array is not specified
     display  *********************************************************
     display
     abort
   else
     query name=&STRIP%obs_f domain=reciprocal end
     if ( $object_exist = false ) then
       display
       display  Error: required observed amplitude array &obs_f does not exist
       display  **************************************************************
       display
       abort
     end if
    {- note: this array can be of any type -}
   end if

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

   evaluate ($reject_obs=&obs_f)
   evaluate ($reject_sig=&obs_sigf)
   evaluate ($obs_lower_limit=0)

   declare name=ref_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
                    ( &low_res >= d >= &high_res ) )

   show sum(1) ( ref_active=1 )
   evaluate ($total_work=$select)
   evaluate ($total_used=$total_work)

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

   tselection=( ref_active=1 )

   cvselection=( none )

   method=FFT          
   
   fft
     grid=&grid
     if ( &fft_memory < 0 ) then
       automemory=true
     else
       memory=&fft_memory
     end if
   end

   tolerance=0.0 lookup=false
   
 end                  

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

 xray
   declare name=map domain=real end
   if ( &fill_in = true ) then
     if ( &BLANK%map_fom = false ) then
       do (map=ft(&STRIP%map_fom*combine(amplitude(&STRIP%map_amplitude),
                                             phase(&STRIP%map_phase))))
          (&low_res >= d >= &high_res)
     else
       do (map=ft(combine(amplitude(&STRIP%map_amplitude),
                              phase(&STRIP%map_phase))))
          (&low_res >= d >= &high_res)
     end if
   else
     if ( &BLANK%map_fom = false ) then
       do (map=ft(&STRIP%map_fom*combine(amplitude(&STRIP%map_amplitude),
                                             phase(&STRIP%map_phase))))
          (ref_active=1)
     else
       do (map=ft(combine(amplitude(&STRIP%map_amplitude),
                              phase(&STRIP%map_phase))))
          (ref_active=1) 
     end if
   end if
 end

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

 xray
   declare name=mask domain=real end
   mask
     average=false
     mode=probe
     solrad=&map_cover
     shrink=0
     nshell=1
     to=mask
     selection=( &atom_map )
   end
   do (map=0) ( mask > 0 )
   undeclare name=mask domain=real end
 end

 evaluate ($cushion=&map_cover*1.5)

 xray
   write map
     if ( &map_format = "ezd" ) then
       type=ezd
     else
       type=cns
     end if
     automatic=false
     from=map
     output=&map_outfile
     cushion=$cushion
     selection=&atom_map
     extend=molecule
   end
 end

 stop