{+ file: optimize_ncsop.inp +}
{+ directory: xtal_phase +}
{+ description: Improve NCS operators for density averaging +}
{+ 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(

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

{* space group *}
{* use International Table conventions with subcripts 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;

{* reflection file *}
{===>} reflection_infile="eg1_abcd.cv";

{* 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 observed intensities: optional *}
{* may be used for reflection rejection *}
{===>} obs_i="";

{* reciprocal space array containing sigma values for intensities: optional *}
{* may be used for reflection rejection *}
{===>} obs_sigi="";

{* experimental phase probability distribution: required *}
{* required for phase combination *}
{* reciprocal space arrays with Hendrickson-Lattman coefficients A,B,C,D *}
{+ table: rows=1 "HL coefficients" cols=4 "A" "B" "C" "D" +}
{===>} obs_pa="pa";
{===>} obs_pb="pb";
{===>} obs_pc="pc";
{===>} obs_pd="pd";

{* resolution limits for improvement *}
{+ table: rows=1 "resolution" cols=2 "lowest" "highest" +}
{===>} low_res=500.0;
{===>} high_res=5.0;

{* 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 < &sigma_cut are rejected. *}
{===>} sigma_cut=0.0;

{* rms outlier cutoff: applied to amplitudes or intensities *}
{* reflections with magnitude(Obs) > &obs_rms*rms(Obs) will be rejected *}
{===>} obs_rms=10000;

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

{* NCS-restraints/constraints file *}
{* only strict NCS is allowed *}
{* see auxiliary/ncs.def *}
{===>} ncs_infile="eg1_ncs_strict.def";

{============================== mask =================================}

{* O mask file *}
{* the mask should cover the primary molecule - the other NCS masks
   will be generated from this and the NCS operators *}
{===>} ncs_mask_infile="eg1.mask";

{======================== search parameters ==========================}

{* rotation parameters *}
{* theta1, theta2 and theta3 define realspace Euler angles.
   maximum, minimum and increment size are required for each angle *}
{+ table: rows=3 "theta1" "theta2" "theta3" 
          cols=3 "minimum" "maximum" "increment" +}

{===>} theta1_min=-0.1;
{===>} theta1_max=0.1;
{===>} theta1_inc=0.1;

{===>} theta2_min=-0.1;
{===>} theta2_max=0.1;
{===>} theta2_inc=0.1;

{===>} theta3_min=-0.1;
{===>} theta3_max=0.1;
{===>} theta3_inc=0.1;

{* translation parameters *}
{* dx, dy, dz define realspace translations in orthogonal Angstroms.
   maximum, minimum and increment size are required for each translation *}
{+ table: rows=3 "dx" "dy" "dz" 
          cols=3 "minimum" "maximum" "increment" +}

{===>} dx_min=-0.5;
{===>} dx_max=0.5;
{===>} dx_inc=0.5;

{===>} dy_min=-0.5;
{===>} dy_max=0.5;
{===>} dy_inc=0.5;

{===>} dz_min=-0.5;
{===>} dz_max=0.5;
{===>} dz_inc=0.5;

{* grid spacing as a fraction of the highest resolution limit *}
{===>} grid=0.25;

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

{* phase step for integration of combined phase probability *}
{===>} np_phistep=5;

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

{* output listing file *}
{===>} list_outfile="optimize_ncs.list";

{===========================================================================}
{         things below this line do not normally need to be changed         }
{===========================================================================}

 ) {- end block parameter definition -}

 checkversion 1.3

 evaluate ($log_level=quiet)

 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

   generate &low_res &high_res

   reflection @&reflection_infile end

   declare name=fobsold domain=reciprocal type=complex end
   declare name=testold domain=reciprocal type=integer end

   do (fobsold=&STRIP%obs_f) (all)

   {- check that all data arrays are well-defined -}
   set echo=off end

   if ( &BLANK%obs_f = true ) then
      display 
      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  **************************************************************
          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  *****************************************************
      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  *************************************************************
          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

   @@CNS_XTALMODULE:check_abcd (pa=&obs_pa;
                               pb=&obs_pb;
                               pc=&obs_pc;
                               pd=&obs_pd;)

   if ( &BLANK%obs_i = false ) then
     query name=&STRIP%obs_i domain=reciprocal end
     if ( $object_exist = true ) then
       do (iobs_orig=&STRIP%obs_i) (all)
     else
        display 
        display  ********************************************
        display  Error: intensity array &obs_i does not exist
        display  ********************************************
        display
        abort
     end if
   end if
   if ( &BLANK%obs_sigi = false ) then
     query name=&STRIP%obs_sigi domain=reciprocal end
     if ( $object_exist = true ) then
       do (sigi_orig=&STRIP%obs_sigi) (all)
     else
        display 
        display  *****************************************************
        display  Error: intensity sigma array &obs_sigi does not exist
        display  *****************************************************
        display
        abort     
     end if
   end if

   set echo=on end
   
   if ( &obs_type = "intensity" ) then
     evaluate ($reject_obs=&obs_i)
     evaluate ($reject_sig=&obs_sigi)
   else
     evaluate ($reject_obs=&obs_f)
     evaluate ($reject_sig=&obs_sigf)
   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) ( ( $STRIP%reject_sig # 0 ) 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 ) )

   do (tst_active=0) (all)
   do (tst_active=1) (none)

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

   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

 end                  

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

 ncs strict ? end

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

 evaluate ($num_op=1)
 evaluate ($done=false)
 while ( $done = false ) loop ncsop
   if ( $exist_ncsop_$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

 if ( $log_level = verbose ) then
   set message=normal echo=on end
 else
   set message=off echo=off end
 end if

 xray
   declare name=m      domain=reci type=complex end  {(fom_obs, phi_obs)}
   declare name=fcomb  domain=reci type=complex end  {combined FOM*Fo,phi}
   declare name=map    domain=real           end  {density map}
   declare name=dm_fom domain=reci type=real end
   declare name=maptmp domain=real           end  {density map}

   @CNS_XTALMODULE:getfom ( pa=&STRIP%obs_pa;
                           pb=&STRIP%obs_pb;
                           pc=&STRIP%obs_pc;
                           pd=&STRIP%obs_pd;
                           m=m;
                           phistep=&np_phistep; )

   do (dm_fom=0) (all)
   do (dm_fom=amplitude(m)) 
      ( (abs(&STRIP%obs_pa) > 0 or abs(&STRIP%obs_pb) > 0 or 
         abs(&STRIP%obs_pc) > 0 or abs(&STRIP%obs_pd) > 0) )

   do (fcomb=0) (all)
   do (fcomb=combine( dm_fom amplitude(&STRIP%obs_f), phase (m) )) (dm_fom > 0)

   do (map=ft(fcomb)) (&low_res >= d >= &high_res and dm_fom > 0)

   declare name=mask domain=real end
   read mask
     to=mask
     type=omask
     input=&ncs_mask_infile
   end

 end

 set display=&list_outfile end
 display ======================== Summary of density averaging ========================
 display
 display input data ---->
 display >>> resolution: &low_res - &high_res
 display >>> sg= &STRIP%sg a= &a b= &b c= &c alpha= &alpha beta= &beta gamma= &gamma
 display >>> reflection file= &STRIP%reflection_infile
 display >>> NCS-strict file= &STRIP%ncs_infile
 display >>> number of NCS operators= $num_op
 display >>> NCS mask for primary molecule read from: &STRIP%ncs_mask_infile
 display
 display data usage ---->
 if ( &obs_type = "intensity" ) then
   display >>> reflections with Iobs/sigma_I < &sigma_cut rejected
   display >>> reflections with Iobs > &obs_rms * rms(Iobs) rejected
 else
   display >>> reflections with |Fobs|/sigma_F < &sigma_cut rejected
   display >>> reflections with |Fobs| > &obs_rms * rms(Fobs) rejected
 end if
 display >>> theoretical total number of refl. in resol. range:    $total_theor[I6] ( 100.0 % )
 display >>> number of unobserved reflections (no entry):          $unobserved[I6] ( $per_unobs[f5.1] % )
 display >>> number of reflections rejected:                       $rejected[I6] ( $per_reject[f5.1] % )
 display >>> total number of reflections used:                     $total_used[I6] ( $per_used[f5.1] % )
 display >>> number of reflections in working set:                 $total_work[I6] ( $per_work[f5.1] % )
 display
 set display=OUTPUT end
 
 evaluate ($ncsop=2)
 while ($ncsop <= $num_op) loop ncs

   set display=&list_outfile end

   display
   display =================================================================== 
   display
   display >>> initial matrix for operator $ncsop
   display     matrix= ( $ncsop_$ncsop_$1_$1 $ncsop_$ncsop_$1_$2 $ncsop_$ncsop_$1_$3 )
   display             ( $ncsop_$ncsop_$2_$1 $ncsop_$ncsop_$2_$2 $ncsop_$ncsop_$2_$3 )
   display             ( $ncsop_$ncsop_$3_$1 $ncsop_$ncsop_$3_$2 $ncsop_$ncsop_$3_$3 )
   display
   display >>> initial translation for operator $ncsop
   display     vector= ( $ncsop_$ncsop_$1_$4 $ncsop_$ncsop_$2_$4 $ncsop_$ncsop_$3_$4 )

   set display=OUTPUT end

   coord
     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 )
   end
   evaluate ($euler_$ncsop_$1=$theta1)
   evaluate ($euler_$ncsop_$2=$theta2)
   evaluate ($euler_$ncsop_$3=$theta3)
   evaluate ($ncsop=$ncsop+1)
 end loop ncs

 evaluate ($ncsop=2)
 while ($ncsop <= $num_op) loop aver

   xray
     do (maptmp=map) (all)
     average
       from=maptmp
       group
         mask=mask
         euler=( 0 0 0 )
         translation=( 0 0 0 )
       end
       group
         euler=
           ( $euler_$ncsop_$1 $euler_$ncsop_$2 $euler_$ncsop_$3 )
         translation=
           ( $ncsop_$ncsop_$1_$4 $ncsop_$ncsop_$2_$4 $ncsop_$ncsop_$3_$4 )
       end
     end
   end

   evaluate ($startcc_$ncsop=$av_corr)
   evaluate ($bestcc_$ncsop=$av_corr)
   evaluate ($bestt1_$ncsop=0.0)
   evaluate ($bestt2_$ncsop=0.0)
   evaluate ($bestt3_$ncsop=0.0)
   evaluate ($bestdx_$ncsop=0.0)
   evaluate ($bestdy_$ncsop=0.0)
   evaluate ($bestdz_$ncsop=0.0)

   set display=&list_outfile end
  
   display >>> starting correlation for 1 -> $ncsop[i3] $av_corr[f6.4]
   
   set display=OUTPUT end

   evaluate ($ncsop=$ncsop+1)

 end loop aver

 set display=&list_outfile end

 display
 display =================================================================== 
 display 
 display  NCS_op     theta1   theta2   theta3     CC
 display
   
 set display=OUTPUT end

 evaluate ($t1=&theta1_min)
 while ( $t1 <= &theta1_max ) loop t1
   evaluate ($t2=&theta2_min)
   while ( $t2 <= &theta2_max ) loop t2
     evaluate ($t3=&theta3_min)
     while ( $t3 <= &theta3_max ) loop t3

       evaluate ($ncsop=2)
       while ($ncsop <= $num_op) loop aver

         evaluate ($ncstmp_1=$euler_$ncsop_$1 + $t1)
         evaluate ($ncstmp_2=$euler_$ncsop_$2 + $t2)
         evaluate ($ncstmp_3=$euler_$ncsop_$3 + $t3)

         xray
           do (maptmp=map) (all)
           average
             from=maptmp
             group
               mask=mask
               euler=( 0 0 0 )
               translation=( 0 0 0 )
             end
             group
               euler=
                 ( $ncstmp_$1 $ncstmp_$2 $ncstmp_$3 )
               translation=
                 ( $ncsop_$ncsop_$1_$4 $ncsop_$ncsop_$2_$4 $ncsop_$ncsop_$3_$4 )
             end
           end
         end

         if ( $av_corr > $bestcc_$ncsop ) then
           evaluate ($bestcc_$ncsop=$av_corr)
           evaluate ($bestt1_$ncsop=$t1)
           evaluate ($bestt2_$ncsop=$t2)
           evaluate ($bestt3_$ncsop=$t3)
         end if

         set display=&list_outfile end
  
         display 1 -> $ncsop[i3]  $t1[f8.4] $t2[f8.4] $t3[f8.4] $av_corr[f8.4]

         set display=OUTPUT end

         evaluate ($ncsop=$ncsop+1)

       end loop aver
       evaluate ($t3=$t3+&theta3_inc)
     end loop t3
     evaluate ($t2=$t2+&theta2_inc)
   end loop t2
   evaluate ($t1=$t1+&theta1_inc)
 end loop t1

 set display=&list_outfile end
  
 display
 display =================================================================== 
 display

 display  NCS_op       dx       dy       dz       CC

 evaluate ($dx=&dx_min)
 while ( $dx <= &dx_max ) loop dx
   evaluate ($dy=&dy_min)
   while ( $dy <= &dy_max ) loop dy
     evaluate ($dz=&dz_min)
     while ( $dz <= &dz_max ) loop dz

       evaluate ($ncsop=2)
       while ($ncsop <= $num_op) loop aver

         evaluate ($ncstmp_1=$euler_$ncsop_$1 + $bestt1_$ncsop)
         evaluate ($ncstmp_2=$euler_$ncsop_$2 + $bestt2_$ncsop)
         evaluate ($ncstmp_3=$euler_$ncsop_$3 + $bestt3_$ncsop)
         evaluate ($ncstmp_4=$ncsop_$ncsop_$1_$4 + $dx)
         evaluate ($ncstmp_5=$ncsop_$ncsop_$2_$4 + $dy)
         evaluate ($ncstmp_6=$ncsop_$ncsop_$3_$4 + $dz)

         xray
           do (maptmp=map) (all)
           average
             from=maptmp
             group
               mask=mask
               euler=( 0 0 0 )
               translation=( 0 0 0 )
             end
             group
               euler=
                 ( $ncstmp_1 $ncstmp_2 $ncstmp_3 )
               translation=
                 ( $ncstmp_4 $ncstmp_5 $ncstmp_6 )
             end
           end
         end

         if ($av_corr > $bestcc_$ncsop) then
           evaluate ($bestcc_$ncsop=$av_corr)
           evaluate ($bestdx_$ncsop=$dx)
           evaluate ($bestdy_$ncsop=$dy)
           evaluate ($bestdz_$ncsop=$dz)
         end if

         set display=&list_outfile end
  
         display 1 -> $ncsop[i3]  $dx[f8.4] $dy[f8.4] $dz[f8.4] $av_corr[f8.4]

         set display=OUTPUT end

         evaluate ($ncsop=$ncsop+1)

       end loop aver
       evaluate ($dz=$dz+&dz_inc)
     end loop dz
     evaluate ($dy=$dy+&dy_inc)
   end loop dy
   evaluate ($dx=$dx+&dx_inc)
 end loop dx

 evaluate ($ncsop=2)
 while ($ncsop <= $num_op) loop aver

   set display=&list_outfile end

   evaluate ($ncstmp_1=$euler_$ncsop_$1 + $bestt1_$ncsop)
   evaluate ($ncstmp_2=$euler_$ncsop_$2 + $bestt2_$ncsop)
   evaluate ($ncstmp_3=$euler_$ncsop_$3 + $bestt3_$ncsop)

   coord
     euler=( $ncstmp_1 $ncstmp_2 $ncstmp_3 )
   end

   evaluate ($ncstmp_4=$ncsop_$ncsop_$1_$4+$bestdx_$ncsop)
   evaluate ($ncstmp_5=$ncsop_$ncsop_$2_$4+$bestdy_$ncsop)
   evaluate ($ncstmp_6=$ncsop_$ncsop_$3_$4+$bestdz_$ncsop)
  
   display
   display =================================================================== 
   display
   display >>> optimized matrix for operator $ncsop
   display >>> correlation coefficient= $bestcc_$ncsop
   display     matrix= ( $rot_1_1[f9.6] $rot_1_2[f9.6] $rot_1_3[f9.6] )
   display             ( $rot_2_1[f9.6] $rot_2_2[f9.6] $rot_2_3[f9.6] )
   display             ( $rot_3_1[f9.6] $rot_3_2[f9.6] $rot_3_3[f9.6] )
   display
   display >>> optimized translation for operator $ncsop
   display >>> correlation coefficient= $bestcc_$ncsop
   display     vector= ( $ncstmp_4[f12.6] $ncstmp_5[f12.6] $ncstmp_6[f12.6] )

   set display=OUTPUT end

   evaluate ($ncsop=$ncsop+1)

 end loop aver

 set display=&list_outfile end
  
 display
 display =================================================================== 

 stop