#!/bin/tcsh -f
#
#
#  Calculate anomalous difference peaks         Version 16 Mar 2005
#
#   A CCP4 map is used in the peaksearch
#   A FSFOUR map can also be calculated for use in XtalView and xcontur
#
#   Assumptions
#   (1) DANO and SIGDANO are present in the mtz file with column names DANO and SIGDANO
#   (2) The ccp4 programs fft, scaleit, mtzdump, mtz2various, sftools, peakmax will be found
#   (3) The phases program fsfour will be found
#
#   Defaults
#   (1) Require F/SIGF > 10 
#
#   Useful but not implemented yet - a -v switch to save all the temporary files
#   F/SIGF test should be an option
#

if (${#argv} < 2 || ${#argv} > 5) then

    set errno = err1
    goto errormessage            # bailing out

else
#
#
#  Set the symbols
#
#

    set mapcheck = ( ccp4 fsfour none both )

    if ( -e ${1} ) then
	set data       = "${1}"
    else
	echo "Couldn't find data file ${1}"
	exit
endif

if (${#argv} == 2) then

    set lowres     = 1000
    set highres    = 0.1

if      ( ${2} == ${mapcheck[1]} ) then
    else if ( ${2} == ${mapcheck[2]} ) then
	else if ( ${2} == ${mapcheck[3]} ) then
	    else if ( ${2} == ${mapcheck[4]} ) then
	else
	    set errno = err6        # bailing out
	    goto errormessage
    endif
    set map        = "${2}"
    set suspect    = "false"

else

    if (${#argv} == 3) then

	set notanumber = `echo "${2}" | grep -c '[^0-9\.]'`
	if ( ${notanumber} == 0 ) then   # second argument is a number (highres)
	    set lowres     = 1000
	    set highres    = ${2}
	if      ( ${3} == ${mapcheck[1]} ) then
	    else if ( ${3} == ${mapcheck[2]} ) then
		else if ( ${3} == ${mapcheck[3]} ) then
		    else if ( ${3} == ${mapcheck[4]} ) then
		else
		    set errno = err6        # bailing out
		    goto errormessage
    endif

    set map        = "${3}"
    set suspect    = "false"

else                         # second argument is a word (ccp4,fsfour,both,none)
    set lowres     = 1000
    set highres    = 0.1
	if      ( ${2} == ${mapcheck[1]} ) then
	    else if ( ${2} == ${mapcheck[2]} ) then
		else if ( ${2} == ${mapcheck[3]} ) then
		    else if ( ${2} == ${mapcheck[4]} ) then
		else
		    set errno = err6        # bailing out
		    goto errormessage
    endif

    set map        = "${2}"
    if ( ${3} == "suspect" ) then
	set suspect    = "true"
    else
	set errno = err2             # bailing out
	goto errormessage
    endif

endif

else

    if (${#argv} == 4) then

  set notanumber = `echo "${3}" | grep -c '[^0-9\.]'`
  if ( ${notanumber} == 0 ) then     # third argument is a number (highres)

    set notanumber = `echo "${2}" | grep -c '[^0-9\.]'`
    if ( ${notanumber} == 0 ) then   # second argument is also a number (lowres)
      set ratio = `awk -v x=${2} -v y=${3} 'BEGIN{ if ((x / y) <= 1) print 0; else print 1}'`
      if ( ${ratio} == 0 ) then
        set errno = err3             # bailing out with lowres less than or equal to highres
        goto errormessage
      endif
    else 
      set errno = err7               # bailing out with lowres or highres not a properly formed numb
er
      goto errormessage
    endif

    set lowres     = ${2}
    set highres    = ${3}
    if      ( ${4} == ${mapcheck[1]} ) then
     else if ( ${4} == ${mapcheck[2]} ) then
      else if ( ${4} == ${mapcheck[3]} ) then
       else if ( ${4} == ${mapcheck[4]} ) then
        else
         set errno = err6            # bailing out
         goto errormessage
    endif
    set map        = "${4}"          # since [suspect] is optional and map is not
    set suspect    = "false"
  else                               # third argument is a word (ccp4,fsfour,both,none)
    set notanumber = `echo "${2}" | grep -c '[^0-9\.]'`
    if ( ${notanumber} == 0 ) then   # second argument is a number
      set lowres     = 1000
      set highres    = ${2}
    else
      set errno = err7               # bailing out with lowres or highres not a properly formed number
      goto errormessage
    endif
    if      ( ${3} == ${mapcheck[1]} ) then
     else if ( ${3} == ${mapcheck[2]} ) then
      else if ( ${3} == ${mapcheck[3]} ) then
       else if ( ${3} == ${mapcheck[4]} ) then
        else
         set errno = err8          # bailing out with ${3} a badly formed number (highres) or map name
         goto errormessage
    endif
    set map        = "${3}"
    if ( ${4} == "suspect" ) then
      set suspect    = "true"
    else
      set errno = err4             # bailing out
      goto errormessage
    endif
  endif

else            # five arguments were found on command line

  set notanumber = `echo "${2}${3}" | grep -c '[^0-9\.]'`
  if ( ${notanumber} == 0 ) then     # second and third arguments are numbers

    set ratio = `awk -v x=${2} -v y=${3} 'BEGIN{ if ((x / y) <= 1) print 0; else print 1}'`
    if ( ${ratio} == 0 ) then
      set errno = err3             # bailing out with lowres less than or equal to highres
      goto errormessage
    endif
  else
    set errno = err9               # bailing out
    goto errormessage
  endif

  set lowres     = ${2}
  set highres    = ${3}
  if      ( ${4} == ${mapcheck[1]} ) then
   else if ( ${4} == ${mapcheck[2]} ) then
    else if ( ${4} == ${mapcheck[3]} ) then
     else if ( ${4} == ${mapcheck[4]} ) then
      else
       set errno = err6            # bailing out
       goto errormessage
  endif
  set map        = "${4}"
  if ( ${5} == "suspect" ) then
    set suspect    = "true"
  else
    set errno = err5               # bailing out
    goto errormessage
  endif

endif   # #argv == 4
endif   # #argv == 3
endif   # #argv == 2

set tmpfile = anomalous_difference$$

#
#   Which map(s) requested
#

switch ( ${map} )

case fsfour:
set fsfourmap = "${map}"
set ccp4map = "${tmpfile}"
breaksw

case ccp4:
set fsfourmap = "${tmpfile}"
set ccp4map = "${map}"
breaksw

case both:
set fsfourmap = "fsfour"
set ccp4map = "ccp4"
breaksw

case none:
set fsfourmap = "${tmpfile}"
set ccp4map = "${tmpfile}"
breaksw

default:
set errno = err6    # map keyword is not set correctly - bailing out
goto errormessage   # should never get here
breaksw

endsw

echo " "
echo "Calculating anomalous-difference maps, please be patient . . ."
echo " "
echo "A CCP4 anomalous-difference map will be searched for peaks above 2.5 sigma"
echo "and the top 200 peaks listed and output to the file peaks.pdb"
echo " "
#
#
#  Construct the PHASES parameter file
#
#
mtzdump     HKLIN ${data} << EOF >& ${tmpfile}.head
HEAD
SYMM
END
EOF
if ( $status != 0 ) then
  echo "Problem with reflection data file in MTZDUMP, exiting"
  echo "See ${tmpfile}.head for MTZDUMP error message(s)"
  exit

endif
 
 
set lattice = `grep "Lattice Type =" ${tmpfile}.head`
set cell    = `awk '/\* Cell Dimensions :/ { getline; getline; print $0 }' ${tmpfile}.head`
set ops     = `grep "Number of Symmetry Operations =" ${tmpfile}.head`
#
echo "LOGFILE=NULL"                    >  ${tmpfile}.parameters
echo "LATTICE=${lattice[${#lattice}]}" >> ${tmpfile}.parameters
echo "${cell}"                         >> ${tmpfile}.parameters
echo "${ops[${#ops}]}"                 >> ${tmpfile}.parameters
#
#
set opcount = 1
while (${opcount} <= ${ops[${#ops}]})
set operation = `grep -E "^ *Symmetry +${opcount} +" ${tmpfile}.head`
shift operation
shift operation
echo "${operation}" >> ${tmpfile}.parameters
@ opcount = ${opcount} + 1
end
#
#
#  Get the actual resolution limits
#
#
set resolution = `grep -E '^( +[0-9\.]+){2} +[(]( +[-0-9\.]+){3}' ${tmpfile}.head`

set dmax = ${resolution[4]}
set dmin = ${resolution[6]}

echo ${dmin} " " ${dmax} " " ${highres} " " ${lowres} > ${tmpfile}.limits

set dmax = `awk '{ if ($2 < $4) print $2; else print $4 }' ${tmpfile}.limits`
set dmin = `awk '{ if ($1 > $3) print $1; else print $3 }' ${tmpfile}.limits`

echo ${dmin} > ${tmpfile}.calc
#
#   Warn if the high reso limit was reset to match the limit of the data found in the mtz file
#
set message = `awk '{ if ($1 > $3) print "Setting high reso limit to " $1 " A"}' ${tmpfile}.limits`
echo "${message}"
echo " "

#   Now for FSFOUR
set spacing = `awk -v x=${dmin} 'BEGIN{ if (x < 2.25) print x / 3.0; else print 0.75 }'`
#
#
#
#  Check for unreasonable data before fft
#    Note this step requires F be present in the input mtz data file
#    First set fft defaults
#
set labin = ""
set toolarge = 130
set flabel = ""
set qlabel = ""

if ( ${suspect} == "true" ) then
set ctypes = `grep " H H H " ${tmpfile}.head`
set clabels = `grep " H K L " ${tmpfile}.head`

set counter = 1
while ( ${counter} <= ${#ctypes} )
if ( ${ctypes[${counter}]} == "F" ) then
set flabel = ${clabels[${counter}]}
@ counter ++
if ( ${ctypes[${counter}]} == "Q" ) then
set qlabel = ${clabels[${counter}]}
set labin = "F2=${flabel} SIG2=${qlabel}"
endif
break
endif
@ counter ++
end

if ( ${qlabel} == "" ) then
echo "Sorry, couldn't find columns labeled F, SIGF in the mtz file"
echo "The exclusions based on F/SIGF and unreasonable values of DANO can't be done"
echo " "
endif

#  
#  Scaleit will accept SIGF=SIGDANO, but without SIGF the F > 3*SIGF test can't be done in fft
#
if ( ${qlabel} != "" ) then

scaleit     HKLIN ${data} << EOF >& ${tmpfile}.analysis
resolution ${dmin} ${dmax}
analyze
labi FP=${flabel} SIGFP=${qlabel} FPH1=${flabel} SIGFPH1=${qlabel} DPH1=DANO SIGDPH1=SIGDANO
EOF
if ( $status != 0 ) then
  echo "Problem in SCALEIT, exiting"
  echo "See ${tmpfile}.analysis for SCALEIT error message(s)"
  exit
endif

set toolarge = `awk '/acceptable differences/ { print $NF }' ${tmpfile}.analysis`

endif   #  exit from block where F,SIGF were found

endif   #  exit from block where suspect == true

echo "Anomalous differences with |DANO| > ${toolarge[${#toolarge}]} will not be used"
echo "in the Patterson synthesis"
echo " "

#
#
#  Calculate a CCP4 map
#
#  F1MAX will be 130 or the value from scaleit
#
#
fft         HKLIN  ${data} \
            MAPOUT ${ccp4map}.map        << EOF >& ${tmpfile}.fftlog
PATTERSON
RESOLUTION  ${dmax} ${dmin}
TITLE       ${dmax} - ${dmin} A anomalous difference Patterson
EXCLUDE     F1MAX ${toolarge[${#toolarge}]}
EXCLUDE     SIG2 10
GRID        SAMPLE 3
LABIN       F1=DANO  SIG1=SIGDANO ${labin}
SCALE       F2 0.000001 0.0
END
EOF
if ( $status != 0 ) then
  echo "Problem in FFT, exiting"
  echo "See ${tmpfile}.fftlog for FFT error message(s)"
  exit
endif

set numlarge = `grep "F1 term too large" ${tmpfile}.fftlog | wc`

echo "${numlarge[1]} anomalous differences were excluded from the fft"
echo "calculation as too large"
echo " "

if ( ${qlabel} != "" ) then

sftools << eof >& ${tmpfile}.sftoolslog
read ${data}
calc col FSIGF = col ${flabel} col ${qlabel} /
select resol > ${dmin} < ${dmax}
select not centro
select col FSIGF <= 10
select col DANO present
purge yes
write ${tmpfile}.hkl format(3i5,2f10.2) col DANO SIGDANO
quit
yes
eof

if ( $status != 0 ) then
  echo "Problem in SFTOOLS, exiting"
  echo "See ${tmpfile}.sftoolslog for SFTOOLS error message(s)"
  exit
endif

set tooweak = `wc ${tmpfile}.hkl`

echo "Reflections with F < 10*SIGF will not be used"
echo "in the Patterson synthesis"
echo " "
echo "Within the resolution limits requested, the input data"
echo "includes ${tooweak[1]} such reflections"
echo " "

endif   # exit from block where F,SIGF were found
#
#
#  Search for peaks above 2.5 sigma
#
#
peakmax     MAPIN  ${ccp4map}.map \
            XYZOUT peaks.pdb << EOF
THRESH      RMS 2.5
NUMPEAKS    200
OUTPUT      BROOKHAVEN
END
EOF

if ( ${fsfourmap} !~ "fsfour" ) then     # we are done
goto cleanup
endif

#
#
#  Save the structure factors in a text file
#    Note - could use EXCLUDE keywords 
#    to try to match the terms used by fft and fsfour
#    EXCLUDE DIFF ${toolarge[${#toolarge}]}
#    and if FP were assigned on LABIN then
#    EXCLUDE SIGP 3 would be possible
#

mtz2various HKLIN  ${data} \
            HKLOUT ${tmpfile}_data.txt << EOF >& ${tmpfile}.mtz2varlog
LABIN       DP=DANO  SIGDP=SIGDANO
MISS        0.0
RESOLUTION  ${dmax} ${dmin}
OUTPUT      USER '(3I4,2F10.2,3X,"0.00")'
MONITOR     20000
END
EOF
#
#
#  FSFOUR can read the file output by mtz2various
# The awk line performs rejections and manipulates the fft coefficients as decribed below
# Note +/- ${toolarge[${#toolarge}]} could be used instead of +/- 130
# Note there is a cutoff of |DANO| > 5*SIGDANO
# Note the DANO**2 terms are weighted by SIGDANO 
#
#
awk '{ if ($4 > -130 && $4 < 130 && ($4 > 0.2 * $5 || $4 < -0.2 * $5)) \
         printf ("%4d%4d%4d %9.2f %9.2f %6.2f\n",  \
                $1, $2, $3, $4 / sqrt ($5), $5, 0) }' < ${tmpfile}_data.txt \
                                                      > ${tmpfile}.phs
#
#
#  Calculate the FSFOUR map
#
#
fsfour << EOF >& ${tmpfile}.fsfourlog
${tmpfile}.parameters
${dmax} - ${dmin} angstrom anomalous Patterson map
0 0 0 0 6 0 8 0 1 ${spacing} ${dmin}
${tmpfile}.phs
${fsfourmap}.map
EOF
if ( $status != 0 ) then
  echo "Problem in FSFOUR, exiting"
  echo "See ${tmpfile}.fsfourlog for FSFOUR error message(s)"
  exit
endif
#
#
#
#
cleanup:
# Comment out this line to see additional log files
rm -f ${tmpfile}*
#
#
endif
exit
#
# Error messages
#
errormessage:
echo " "
echo " Usage:" ${0:t} "data [[lowres] highres] map [suspect]"
echo " You must specify the mtz file and which map(s) to save"
echo " "
echo "        Calculates anomalous-difference Patterson maps"
echo "          and optionally saves map(s) for viewing"
echo " "
echo "            data        CCP4 mtz file"
echo "            lowres      low-resolution limit [defaults to the low reso limit of the data]"
echo "            highres     high-resolution limit [defaults to the high reso limit of the data]"
echo "            map         ccp4, fsfour, both, or none"
echo "            suspect     don't use suspect data [default is to use all data with |DANO| < 130]"
echo " "

switch ( ${errno} )

case err1:
breaksw

case err2:
echo " "
echo " Please check your command and try again"
echo " ${3} is not a keyword"
echo " "
breaksw

case err3:
echo " "
echo " Please check your command and try again"
echo " The low reso limit should be greater than the high reso limit"
echo " "
breaksw

case err4:
echo " "
echo " Please check your command and try again"
echo " ${4} is not a keyword"
echo " "
breaksw

case err5:
echo " "
echo " Please check your command and try again"
echo " ${5} is not a keyword"
echo " "
breaksw

case err6:
echo " "
echo " Keyword map must be ccp4, fsfour, both, or none"
echo " Specify "none" if no map should be saved"
echo " "
breaksw

case err7:
echo " "
echo " Please check your command and try again"
echo " Keyword ${2} should be a number"
echo " "
breaksw

case err8:
echo " "
echo " Please check your command and try again"
echo " ${3} is not a keyword"
echo " "
breaksw

case err9:
echo " " 
echo " Please check your command and try again"
echo " Keywords ${2} and ${3} should both be numbers"
echo " "
breaksw

case default:
breaksw

endsw

exit