{+ file: average_map.inp +}
{+ directory: xtal_util +}
{+ description: Single cycle realspace averaging of a 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(

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

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

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

{* coordinate file *}
{===>} coordinate_infile="eg1_dimer.pdb";

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

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

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

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

{* masks can be derived from the input atomic coordinates, or can be
   read in from an O format mask file. One or the other - do not mix *}

{* masks from coordinates: select atoms which form the primary protomer
   in each NCS group. In the strict NCS case only one group can be
   defined. *}

{* select atoms in NCS group 1 *}
{===>} group_1=(none);

{* select atoms in NCS group 2 *}
{===>} group_2=(none);

{* select atoms in NCS group 3 *}
{===>} group_3=(none);

{* select atoms in NCS group 4 *}
{===>} group_4=(none);

{* masks from file: give the O format mask files which form the primary
   protomer in each NCS group. In the strict NCS case only one mask
   can be defined. *}

{* NCS group 1: O format mask file *}
{===>} mask_infile_1="eg1.mask";

{* NCS group 2: O format mask file *}
{===>} mask_infile_2="";

{* NCS group 3: O format mask file *}
{===>} mask_infile_3="";

{* NCS group 4: O format mask file *}
{===>} mask_infile_4="";

{====================== 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="average_map.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

   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

 evaluate ($1=1)
 evaluate ($2=2)
 evaluate ($3=3)
 evaluate ($4=4)

 evaluate ($ngroup=1)
 evaluate ($group=1)
 evaluate ($done=false)
 while ( $done = false ) loop group
   if ( &exist_group_$group = true ) then
     show sum(1) ( &group_$group )
     if ( $result > 0 ) then
       evaluate ($ngroup=$ngroup+1)
     end if
   else
     evaluate ($done=true)
     evaluate ($ngroup=$ngroup-1)
   end if
   evaluate ($group=$group+1)
 end loop group

 if ( $ngroup <= 0 ) then
   evaluate ($ngroup=1)
   evaluate ($done=false)
   while ( $done = false ) loop group
     if ( &exist_mask_infile_$ngroup = true ) then
       if ( &BLANK%mask_infile_$ngroup = true ) then
         evaluate ($done=true)
         evaluate ($ngroup=$ngroup-1)
       else
         evaluate ($ngroup=$ngroup+1)
       end if
     else
         evaluate ($done=true)
         evaluate ($ngroup=$ngroup-1)
     end if
   end loop group
 end if

 if ($ngroup <= 0) then
   display No NCS groups found, aborting
   abort
 end if

 if (&ncs_type="strict") then
   ncs strict ? end
   evaluate ($num_op_1=$ncs)
 elseif (&ncs_type="restrain") then
   ncs restraint ? end
   evaluate ($group=1)
   while ($group <= $ngroup) loop group
     evaluate ($num_op=1)
     evaluate ($done=false)
     while ( $done = false ) loop ncsop
       if ( $exist_rot_$group_$num_op_$1_$1 # true ) then
         evaluate ($done=true)
         evaluate ($num_op=$num_op-1)
       else
         evaluate ($num_op=$num_op+1)
       end if
     end loop ncsop
     evaluate ($num_op_$group=$num_op)
     evaluate ($group=$group+1)
   end loop group
 end if

 evaluate ($group=1)
 while ($group <= $ngroup) loop mask
   evaluate ($maskname="mask" + encode($group))
   if ( &BLANK%mask_infile_$group = false ) then
     xray
       declare name=$maskname domain=real end
       read mask
         to=$maskname
         type=omask
         input=&mask_infile_$group
       end
     end
   else
     xray
       declare name=$maskname domain=real end
       mask
         average=true
         mode=vdw
         solrad=1.0
         shrink=1.0
         nshell=1
         to=$maskname
         sele=( &group_$group )
       end
     end
   end if
   evaluate ($group=$group+1)
 end loop mask

 xray
   evaluate ($group=1)
   while ($group <= $ngroup) loop average
     evaluate ($maskname="mask" + encode($group))
     average
       from=map
       evaluate ($ncsop=1)
       while ($ncsop <= $num_op_$group) loop ncsop
         group
           if ( $ncsop = 1 ) then
             mask=$maskname
           end if
           if (&ncs_type="strict") then
             matrix=
               ( $ncsop_$ncsop_$1_$1 $ncsop_$ncsop_$1_$2 $ncsop_$ncsop_$1_$3 )
               ( $ncsop_$ncsop_$2_$1 $ncsop_$ncsop_$2_$2 $ncsop_$ncsop_$2_$3 )
               ( $ncsop_$ncsop_$3_$1 $ncsop_$ncsop_$3_$2 $ncsop_$ncsop_$3_$3 )
             translation=
               ( $ncsop_$ncsop_$1_$4 $ncsop_$ncsop_$2_$4 $ncsop_$ncsop_$3_$4 )
           elseif (&ncs_type="restrain") then
             matrix=
               ( $rot_$group_$ncsop_$1_$1
                 $rot_$group_$ncsop_$1_$2
                 $rot_$group_$ncsop_$1_$3 )
               ( $rot_$group_$ncsop_$2_$1
                 $rot_$group_$ncsop_$2_$2
                 $rot_$group_$ncsop_$2_$3 )
               ( $rot_$group_$ncsop_$3_$1
                 $rot_$group_$ncsop_$3_$2
                 $rot_$group_$ncsop_$3_$3 )
             translation=
               ( $rot_$group_$ncsop_$1_$4
                 $rot_$group_$ncsop_$2_$4
                 $rot_$group_$ncsop_$3_$4 )
           end if
         end
         evaluate ($ncsop=$ncsop+1)
       end loop ncsop
     end
     display     NCS group $group :     average correlation= $av_corr[f6.4]
     evaluate ($ncsop=1)
     while ($ncsop <= $num_op_$group) loop print
       display       NCS operator $ncsop ( 1 -> $ncsop ): correlation= $av_corr_op_$ncsop[f6.4]
       evaluate ($ncsop=$ncsop+1)
     end loop print
     evaluate ($group=$group+1)
   end loop average
 end

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

 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

 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