{+ file: mask_map.inp +}
{+ directory: xtal_util +}
{+ description: Mask out regions of an electron density map +}
{+ 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
   - the selections store1 through store4 are available for general use -}

{- begin block parameter definition -} define(

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

{* coordinate file *}
{===>} coordinate_infile="eg1_dimer.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="eg1_dimer.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="P4(1)2(1)2";

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

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

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

{=========================== input map ===============================}

{* input electron density map *}
{===>} map_infile="eg1_dimer.map";

{* resolution limits used in map calculation *}
{* this must be the same as the limits used to calculate the map *}
{+ table: rows=1 "resolution" cols=2 "lowest" "highest" +}
{===>} low_res=500.0;
{===>} high_res=4.0;

{* map grid size: dmin*grid *}
{* this must be the same as the grid used to calculate the map *}
{===>} grid=0.33;

{* memory allocation for FFT calculation *}
{* this must be the same as the FFT memory allocation used
   to calculate the map *}
{===>} fft_memory=-1;

{======================== masking parameters =========================}

{* select atoms to be included in mask calculation *}
{* a mask will be made around all selected atoms, including
   atoms generated by any strict NCS operators defined above *}
{===>} atom_mask=(known and not hydrogen);

{* optional O-mask for masking out map *}
{* if blank the mask will be calculated using the current model *}
{===>} mask_infile="";

{* zero density in macromolecule or solvent region *}
{+ choice: "macromolecule" "solvent" +}
{===>} mask_region="macromolecule";

{====================== output map parameters ========================}

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

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

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

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

{* atom selection for writing map around molecule *}
{* the map will cover the selected atoms *}
{===>} atom_map=(known and not hydrogen);

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

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

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

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

   binresolution &low_res &high_res
   mapresolution &high_res

   method=FFT

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

 end

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

 xray
   declare name=map domain=real end
   read map
     to=map
     input=&map_infile
   end
 end

 if ( &BLANK%mask_infile = false ) then
   xray
     declare name=mask domain=real end
     read mask
       to=mask
       type=omask
       input=&mask_infile
     end
   end
 else
   xray
     declare name=mask domain=real end
     mask
       average=true
       mode=vdw
       solrad=1.0
       shrink=1.0
       nshell=1
       to=mask
       sele=( &atom_mask )
     end
   end
 end if

 xray
   if ( &mask_region = "macromolecule" ) then
     do (map=0.0) ( real(mask) <= 0 )
   else
     do (map=0.0) ( real(mask) > 0 )
   end if
 end

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

 remark a= &a b= &b c= &c alpha= &alpha beta= &beta gamma= &gamma sg= &STRIP%sg

 xray
   write map
     auto=false
     from=map
     output=&map_outfile
     extent=$map_mode_string
     if ( &map_format = "ezd" ) then
       type=ezd
     else
       type=cns
     end if
     if ( &map_mode = "molecule" ) then
       selection=&atom_map
       cushion=&map_cushion
     end if
     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

 stop