#!/bin/csh -f
#
##########################################################
##########################################################
#
# AVERAGING PROCEDURE USING SOLOMON
#
# Script to run several cycles of density modification, in
# this case with a flip value (solvmul) of -1.0. This 
# corresponds to adding the inverted solvent density to 
# the initial map. This is similar to adding negative noise
# to an image in order to strengthen the signal/noise ratio.
#
# The resolution of the map is extended from 3.0 to 2.7
# The variable $res defines the resolution at which the
# current refinement cycle is working. This outlines the 
# general program order that is required but obviously the 
# procedure will be different for other structures. Automatic 
# refinement over a large number of cycles is not advisable, 
# the increase in map quality is the most important factor.
#
# Note that the 'final' map from each cycle is kept i.e. 
# the ones calculated from the combined structure factors.
# Along with the MTZ file that generated it and the solomon
# log files are appended.
#
##########################################################
#


if ! (-r $CCP4_SCR/toxd_phase_mir.mtz) then
  echo '! run the mlphare procedure first'
  exit 1
endif

if (-r toxd_stats.log) then
  /bin/rm -f toxd_stats.log
endif

if (-r toxd_Solomon.log) then
  /bin/rm -f toxd_Solomon.log
endif


set solvdens=0
set EndCycle=8          # Number of refinement cycles

set res=(3.0 3.0 3.0 2.9 2.8 2.7 2.7 2.7)
#
#
set LastCycle=0
set CurrentCycle=$LastCycle
#
#
#
##########################################################
#  Calculate initial map with best phases and FOM weighted 
#  amplitudes.
##########################################################
#
#
fft HKLIN $CCP4_SCR/toxd_phase_mir.mtz \
    MAPOUT $CCP4_SCR/toxd_cycle0.map \
    << end-fft   > /dev/null
LABIN F1=FTOXD3 SIG1=SIGFTOXD3 PHI=PHI_mir W=W_mir
RESO 9999 $res[1]
GRID 88 144 80
#XYZLIM 0 1.0 0 1.0 0 1.0
end
end-fft
#
#
Begin:
#
@ CurrentCycle=$LastCycle + 1
#
##########################################################
#  Extend the map resolution for several cycles
##########################################################
#
#
##########################################################
#  1) Modify the density with Solomon
##########################################################
#
#
Solomon:
#
solomon \
mapin $CCP4_SCR/toxd_cycle${LastCycle}.map mapout $CCP4_SCR/toxd_av.map \
<< end-solomon  >> toxd_Solomon.log
slvfrc 0.51
mulsolv auto -1.0
slvdens $solvdens
radius 3
trunc 0.4 1
extrude
mapout

!Put any symmetry constraints here...

end-solomon
#
set solvdens=`grep "solvdens:" toxd_Solomon.log|tail -1|awk '{print $3}'`
#
#
##########################################################
#  2) generate sfs and phases from this modified map
##########################################################
#
#
sfall HKLIN $CCP4_SCR/toxd_phase_mir \
      MAPIN $CCP4_SCR/toxd_av.map\
      HKLOUT $CCP4_SCR/junk.mtz \
      << end-sfrkall > $CCP4_SCR/toxd_sfall.log
TITLE   Sfs from density modification
MODE SFCALC MAPIN HKLIN
GRID 88 144 80
LABIN FP=FTOXD3  SIGFP=SIGFTOXD3 F0=PHI_mir F1=W_mir F2=HLA F3=HLB -
      F4=HLC F5=HLD
RESO 9999 $res[${CurrentCycle}]
RSCB 6 $res[${CurrentCycle}]
BINS 40
SFSG 19
LABO   FC=FCmod PHIC=PHICmod
end
end-sfrkall
#
echo Cycle $CurrentCycle >> toxd_stats.log
grep "Overall Reliability index is" $CCP4_SCR/toxd_sfall.log \
						 >> toxd_stats.log
#
#
##########################################################
#  3) combine phases
##########################################################
#
sigmaa HKLIN $CCP4_SCR/junk.mtz \
       HKLOUT $CCP4_SCR/toxd_cycle${CurrentCycle}.mtz \
       << end-sigmaa > /dev/null
TITLE   Phase combination toxd
LABI FP=FTOXD3  SIGFP=SIGFTOXD3                  -
     HLA=HLA HLB=HLB HLC=HLC HLD=HLD PHIBP=PHI_mir  WP=W_mir     -
     FC=FCmod PHIC=PHICmod
LABO PHCMB=PHCMB1 WCMB=FOMCMB1 FWT=FWT PHFWT=PHFWT
RANGES 10 1000
RESO  9999 $res[$CurrentCycle]
ERROR
COMBINE PART 1 
END
end-sigmaa
#
#
##########################################################
#  4) recalculate map.
##########################################################
#
fft HKLIN $CCP4_SCR/toxd_cycle${CurrentCycle}.mtz \
    MAPOUT $CCP4_SCR/toxd_cycle${CurrentCycle}.map \
    << end-fft > /dev/null
TITLE F1 15-$res A, after ${CurrentCycle} averaging cycles
LABI F1=FWT  SIG1=FWT PHI=PHFWT
RESO 9999 $res[$CurrentCycle]
GRID 88 144 80
#XYZLIM 0 1.0 0 1.0 0 1.0
end
end-fft
#
#
@ LastCycle=$LastCycle + 1
#
#
if ( $CurrentCycle < $EndCycle ) goto Begin
#
exit 0