#!/bin/sh

#
# Procedure  for finding water peaks. 
# An atom map generated with SFALL is used to mask out the protein regions in
# both an fo-fc and 2fo-fc map.  These solvent maps are then averaged and 
# searched for peaks.  
#
# fofcmap script does the same thing but only uses fo-fc map
#
set -e

########################################################################
# Generate atom map for final coordinates
########################################################################
sfall    			\
MAPOUT $CCP4_SCR/toxd_fc.map \
XYZIN $CEXAM/toxd/toxd.pdb \
<< END-sfall 
MODE ATMMAP RESMOD
FORM NGAU 2
GRID 152 80 48
SFSG 19
symm 19
END-sfall
########################################################################
# Calculate sfs for final coordinates -
########################################################################
sfall    			\
HKLIN $CEXAM/toxd/toxd.mtz \
HKLOUT $CCP4_SCR/toxd.mtz  \
XYZIN $CEXAM/toxd/toxd.pdb \
<< END-sfall 
MODE SFCALC XYZIN HKLIN
FORM NGAU 2
RESO  10 2.0
BINS 40
RSCB  6 2.0
BADD 0
LABI FP=FTOXD3 SIGFP=SIGFTOXD3
LABO FC=FC PHIC=AC
END-sfall
#
########################################################################
# Fo-Fc map from final coordintes.  
########################################################################
#
fft \
HKLIN $CCP4_SCR/toxd.mtz  \
MAPOUT $CCP4_SCR/toxd_fo-fc.map  \
<< END-fft
RESO 10 2.0
SCALE F1 1.0 0.0
SCALE F2 1.0 0.0
GRID 152 80 48 ! must be same as used for calc Fc map
XYZLIM 0 151 0 79 0 12
BINMAPOUT
LABI F1=FTOXD3 SIG1=SIGFTOXD3 F2=FC PHI=AC
END-fft
#
########################################################################
# 2Fo -fc map from final coordintes.  
########################################################################
#
#
fft \
HKLIN $CCP4_SCR/toxd.mtz  \
MAPOUT $CCP4_SCR/toxd_2fo-fc.map  \
<< END-fft
RESO 10 2.0
SCALE F1 2.0 0.0
SCALE F2 1.0 0.0
GRID 152 80 48 ! must be same as used for calc Fc map
XYZLIM 0 151 0 79 0 12
BINMAPOUT
LABI F1=FTOXD3 SIG1=SIGFTOXD3 F2=FC PHI=AC
END-fft
#
#
########################################################################
# exclude points from 2fo-fc map which occur in fcalc map
########################################################################

overlapmap \
MAPIN1 $CCP4_SCR/toxd_fc.map  \
MAPIN2 $CCP4_SCR/toxd_2fo-fc.map  \
MAPOUT $CCP4_SCR/toxd_2fo-fc_masked.map  \
<< end-over
MAP EXCLUDE
end-over
#
########################################################################
# exclude points from fo-fc map which occur in fcalc map
########################################################################
#
overlapmap \
MAPIN1 $CCP4_SCR/toxd_fc.map  \
MAPIN2 $CCP4_SCR/toxd_fo-fc.map  \
MAPOUT $CCP4_SCR/toxd_fo-fc_masked.map  \
<< end-over
MAP EXCLUDE
end-over
#
########################################################################
# average 2 solvent maps now
########################################################################
#
overlapmap \
 MAPIN1  $CCP4_SCR/toxd_fo-fc_masked.map \
 MAPIN2  $CCP4_SCR/toxd_2fo-fc_masked.map \
 MAPOUT  $CCP4_SCR/toxd_av_masked.map \
<< end-over
MAP AVERAGE
end-over
########################################################################
# run peakmax with electron density level 0.15  - look at 2fo-fc, fo-fc 
#                                                 to get sensible level
########################################################################

#  run mapmask to cover molecule
mapmask    \
XYZIN $CEXAM/toxd/toxd.pdb \
 MAPIN  $CCP4_SCR/toxd_av_masked.map \
 MAPOUT  $CCP4_SCR/toxd_av_masked.ext \
<< END 
border 4.0
mode MAPIN
end
END
#
########################################################################
#  Some peaks may appear twice since map has been extended.
########################################################################

peakmax                 \
 XYZOUT $CCP4_SCR/toxd_peaks.pdb \
 MAPIN  $CCP4_SCR/toxd_av_masked.ext \
<< END-peakmax
NUMPEAKS 1200
THRESHOLD rms 1.5
END-peakmax
#
########################################################################
# run atpeak to check distances within 3.6A
########################################################################
#

watpeak                  \
PEAKS $CCP4_SCR/toxd_peaks.pdb \
XYZIN $CEXAM/toxd/toxd.pdb \
xyzout $CCP4_SCR/toxd_waters.pdb \
<<END-watpeak
TITLE  Toxd contacts
DISTANCE 3.6
SYMMETRY 19
!HETATOMONLY ! use if only distances from O and N wanted
END-watpeak
#