#!/bin/sh

# ============================================================
#       Outline of Heavy atom refinement and Phasing tutorial 
#            ( E.J.Dodson)
# ============================================================
#        Example based on MIR solution of Dendrotoxin from green mamba
#        Protein of 425 residues 
#        cell 73.58   38.73   23.19   90.00   90.00   90.00
#        Spacegroup  P212121
#    Two derivatives; Au I

# Steps after solving patterson for 1 site.
#    Use FFT, PEAKMAX  HAVECS, Int Tab etc..
#========================
# step 1.
#========================
#   Refine Au Site 1 using Centric differences only
#  Phasing cycle ( the last) will calculate SIR phases for all reflections.
#  Alternative:  Use a direct methods program to find a consistent set of Hg 
#     sites and go straight to Step 3.
#
# Step 2.
#========================
#   Difference Fourier to find more sites. Contour and
#   note pertubation in map about first site.
#   Check patterson peaks are consistent with 2nd site. ( do not worry 
#   too much if they are not... (gets round alternate origin problem)
#
# Step 3.
#========================
#   Refine 2 Au Sites  using Centric differences only
#  Phasing cycle ( the last) will calculate SIR phases for all reflections
#
# Step 4.
#========================
#   Difference fouriers to look for I sites.
#
# Step 5.
#========================
#  Refinement of both derivatives together.
#
# Step 6.
#========================
#  Refinement of both derivatives together including anomalous data
#    If hand correct, and there are TWO independent derivatives; 
#    anomalous occupancies will be positive if hand correct;
#                                   negative if hand incorrect.
#       If so - change positions to -x,-y,-z, ( spacegroup to 
#        enantiomorph if necessary) and go on..
#
#========================
#   Next steps not in this example:
#
# Step 7.
#========================
#  Check I pattersons.
#  Now should repeat procedure, using I phases to find Au.
#  Calculate isomorphous map...
#
# Step 8.
#========================
#  Calculate isomorphous map...
#
# Step 9.
#========================
# Do density modification. (dm/Wang procedures)
#  Calculate new map...
#
# Step 10.
#========================
# Build structure.
#
#
###   Command procedures  to carry out this procedure. 
#    ==================================================
#                          STEP 1
#    ==================================================
# step 1.
#   Refine Au Site 1 using Centric differences only
#  Phasing cycle ( the last) will calculate SIR phases for all reflections.
#----------------------------------------------------
#----------------------------------------------------
#    step1_mir_steps
#----------------------------------------------------
#----------------------------------------------------

set -e

if test ! -f $CCP4_SCR/toxd_sc.mtz; then
  echo '! run scaleit.exam first' 1>&2
  exit 1
fi

mlphare HKLIN $CCP4_SCR/toxd_sc.mtz \
       HKLOUT   $CCP4_SCR/toxd_step1.mtz        \
 << eof-phare
TITLE Refining cycles only
RESO 15.0 2.7
CENTRIC
SCALE SIGFP 1.0 0.0
CYCLE 15
THRES 2.5 0.5
ANGLE 10
PRINT AVF AVE MON 1000
LABIN FP=FTOXD3 SIGFP=SIGFTOXD3 -
      FPH1=FAU20 SIGFPH1=SIGFAU20
LABOUT ALLIN PHIB=PHI_step1 FOM=FOM_step1
FHOUT DERIV 1   ! For double difference map.
 DERIV AU20 1st Site
 DCYCLE PHASE ALL REFCYC ALL KBOV ALL
 RESO 10.0 2.7
! ISOE     1575.10 1685.89 1100.50  959.22 1019.41 1069.37  938.56  708.38
! ANOE       95.00  338.73  348.88  326.07  288.25  273.30  296.62  242.99
 ATOM   AU    0.1758  0.1016 -0.10  1.00  BFAC   25.000
 ATREF X ALL Y ALL Z ALL OCC ALL 
END
eof-phare
#

#   ==================================================
#                         STEP 2
#   ==================================================
#
# Step 2.
#   Difference Fourier to find more sites. Contour and
#   note pertubation in map about first site.
#   Check patterson peaks are consistent with 2nd.. sites.  (do not worry 
#   too much if they are not...)
#   (If you add FH=FH1 PHIH=PHIH1 you will generate a double difference map 
#      - may be clearer..)
#----------------------------------------------------
#----------------------------------------------------
#    step2_mir_steps
#----------------------------------------------------
#----------------------------------------------------
#

fft \
HKLIN $CCP4_SCR/toxd_step1.mtz \
MAPOUT $CCP4_SCR/toxd_step2.map	\
abcoeffs $CCP4_SCR/fftkw.abcoeffs \
<< 'END_fft' 
RESO 37. 2.7
TITL  Fourier map to find other sites.
GRID 128 64 40
XYZLIM 0 127 0 63 0 10
LABI F1=FAU20 SIG1=SIGFAU20 F2=FTOXD3 SIG2=SIGFTOXD3 -
PHI=PHI_step1 W=FOM_step1
END_fft
#

# Find 100 peaks above 5*rms, 
# output to PDB file "occupancy"= peak height, " Btemp"=   height/sigma
# orthogonalization code (normal)
# positive peaks only

peakmax mapin $CCP4_SCR/toxd_step2.map xyzout $CCP4_SCR/toxd_step2.frc  << END
NUMPEAKS 100 
THRESHOLD RMS 5
OUTPUT BROOKHAVEN
END

#   ==================================================
#                         STEP 2A
#   ==================================================
#
#  Get PDB file of coordinates and vectors.
#  Alternative inputs - PDB,fractional or phare_ml

# 1)   Input fractional sites 
#   Free format:   Occup   x y z b (sites.frc - typed in)
#
# 2)  Input PHARE fractional sites 
#ATOM   Au    0.026  0.000  0.287 19.75 25.74 BFAC   13.887
#
# 3) Input PDB file - usual sort of stuff CRYSTAL/SCALEi coordinates.
#
havecs				\
XYZIN $CEXAM/unix/runnable/sites.frc		\
XYZOUT $CCP4_SCR/sites.pdb		\
UVWOUT $CCP4_SCR/sitesuvw.pdb		\
<< 'END-havecs' 
19                                               ! Space Group 
phare                                            ! Input format
73.53    39.06    23.15   90.00  90.00  90.00  1 ! cell and ncode
128  64  40                          ! grid  (x,y,z)  of Difference Fourier
0 64 0 32 0 20                       ! grid limits of asymmetric unit
END-havecs

#   Plot map with Au sites marked.  Choose sensible contour levels
#
#!/bin/csh -f
npo		\
MAPIN $CCP4_SCR/toxd_step2.map	\
XYZIN1 $CCP4_SCR/sites.pdb \
PLOT $CCP4_SCR/step2.plot		\
<< 'END-npo' 
MAP SCALE 2.0
CONTRS NEG SIG 1.0 to 4.0 by 1.0
MODE BELOW 0.0 DASHED 2 1 0
SECTNS 1 10 1
INPUT BROOK
RESIDUE RED SELECT ALL
LABEL FIXED ALL 
JOIN RADII RESIDUE Au 0.01    ! stop any bonds being draw
SOLID
PLOT
END-npo

#   ==================================================
#                         STEP 3
#   ==================================================
#
# Step 3.
#========================
#   Refine 2 Au Sites  using Centric differences only
#  Phasing cycle ( the last) will calculate SIR phases for all reflections
#
#----------------------------------------------------
#----------------------------------------------------
#    step3_mir_steps 
#----------------------------------------------------
#----------------------------------------------------
#
mlphare			\
HKLIN $CCP4_SCR/toxd_sc.mtz \
HKLOUT $CCP4_SCR/toxd_step3.mtz \
DATOUT test.dat	\
<< 'END-phare'
TITLE Ref and Phasing on  2 possible Au sites
RESO 100.0 2.7
CENTRIC
ANGLE 10
THRESHOLD 2.5 0.5
PRINT AVF AVE MON 1000 
CYCLES 10
LABI FP=FTOXD3 SIGFP=SIGFTOXD3 -
FPH1=FAU20 SIGFPH1=SIGFAU20 
LABO  ALLIN PHIB=PHI_step3 FOM=FOM_step3  

DERIV   
DCYCLE PHASE ALL REFCYC ALL KBOV ALL
ISOE     100 100 100 100 100 100 100 100 100 100 
ATOM   AU    0.177  0.104 -0.114  9.92
ATREF X ALL Y ALL Z ALL OCC ALL 
ATOM   AU    0.218  0.138 -0.105  4.88
ATREF X ALL Y ALL Z ALL OCC ALL 
END
END-phare

#==================================================
#                         STEP 4
#==================================================
#
# Step 4.
#========================
#   Difference fouriers to look for I sites and Hg sites

#----------------------------------------------------
#    step4a_mir_steps
#----------------------------------------------------
#----------------------------------------------------
#    In this case it was important to EXCLUDE ridiculous differences.
#
#!/bin/csh -f
fft \
HKLIN $CCP4_SCR/toxd_step3.mtz \
MAPOUT $CCP4_SCR/step4.map	\
abcoeffs $CCP4_SCR/fftkw.abcoeffs \
<< 'END_fft' 
RESO 10 2.7
EXCLUDE DIFF  1300
FFTSPG 19 
TITL  Fourier map  to find I sites.
grid 128 64 40  
xyzlim 0 127 0 63 0 10 
LABI F1=FI100 SIG1=SIGFI100 F2=FTOXD3 SIG2=SIGFTOXD3 -
PHI=PHI_step3 W=FOM_step3
END_fft
#
#
peakmax    \
MAPIN $CCP4_SCR/step4.map	\
XYZOUT $CCP4_SCR/siteI.frc  \
<< 'END' 
NUMPEAKS 100 
THRESHOLD RMS 4.0
END
#
#----------------------------------------------------
#----------------------------------------------------
#    step4b_mir_steps
#----------------------------------------------------
#----------------------------------------------------
#
#
#!/bin/csh -f
fft \
HKLIN $CCP4_SCR/toxd_step3.mtz \
MAPOUT $CCP4_SCR/step4b.map	\
abcoeffs $CCP4_SCR/fftkw.abcoeffs \
<< 'END_fft' 
RESO 10 2.8
EXCLUDE DIFF  3000
FFTSPG 19 
TITL  Fourier map  to find Hg sites.
grid 128 64 40  
xyzlim 0 127 0 63 0 10 
LABI F1=FMM11 SIG1=SIGFMM11 F2=FTOXD3 SIG2=SIGFTOXD3 -
PHI=PHI_step3 W=FOM_step3
END_fft
#
#
peakmax    \
MAPIN $CCP4_SCR/step4b.map	\
XYZOUT $CCP4_SCR/siteHg.frc  \
<< 'END' 
NUMPEAKS 100 
THRESHOLD RMS 8
END

#==================================================
#                         STEP 5
#==================================================
#
# Step 5.
#========================
#  Refinement of all derivatives together.
#  Include refinement of anomalous occupancies.
#----------------------------------------------------
#----------------------------------------------------
#    step5_mir_steps
#----------------------------------------------------
#----------------------------------------------------

mlphare			\
HKLIN $CCP4_SCR/toxd_sc.mtz \
HKLOUT $CCP4_SCR/toxd_step5.mtz \
DATOUT $CCP4_SCR/test.dat	\
<< 'END-phare'
TITLE Ref and Phasing
RESO 100.0 2.7
#CENTRIC
ANGLE 10
THRESHOLD 2.5 0.5
PRINT AVF AVE MON 1000 
CYCLES 10
LABI FP=FTOXD3 SIGFP=SIGFTOXD3 -
FPH1=FAU20 SIGFPH1=SIGFAU20  DPH1=ANAU20 SIGDPH1=SIGANAU20 -
FPH2=FI100 SIGFPH2=SIGFI100 -
FPH3=FMM11 SIGFPH3=SIGFMM11 
LABO  ALLIN PHIB=PHI_step5 FOM=FOM_step5  

DERIV     Gold
DCYCLE PHASE ALL REFCYC ALL KBOV ALL
#ISOE      746.82  548.14  407.42  365.98  379.86  329.45  346.91  345.99
ATOM AU    0.177  0.104 -0.114  9.92 BFAC   25.000
ATREF X ALL Y ALL Z ALL OCC ALL AOCC ALL
ATOM  AU    0.218  0.138 -0.105  4.88 BFAC   25.000
ATREF X ALL Y ALL Z ALL OCC ALL AOCC ALL

DERIV     Iodine
DCYCLE PHASE ALL REFCYC ALL KBOV ALL
#ISOE     100 100 100 100 100 100 100 100 100 100 
ATOM  I     0.491  0.370  0.487  8.40 BFAC   25.000
ATREF X ALL Y ALL Z ALL OCC ALL 

DERIV     Mercury 
DCYCLE PHASE ALL REFCYC ALL KBOV ALL
#ISOE     100 100 100 100 100 100 100 100 100 100 
ATOM  HG     0.178  0.297  0.103  12.0 BFAC   25.000
ATREF X ALL Y ALL Z ALL OCC ALL 
ATOM  HG     0.678  0.401  0.124  2.0 BFAC   25.000
ATREF X ALL Y ALL Z ALL OCC ALL 

END
END-phare
#

#   ==================================================
#                         STEP 6
#   ==================================================
#  Refinement of both derivatives together.
#  Include refinement of anomalous occupancies on both hands.
#----------------------------------------------------
#----------------------------------------------------
#    step6_mir_steps
#----------------------------------------------------
#----------------------------------------------------

mlphare			\
HKLIN $CCP4_SCR/toxd_sc.mtz \
HKLOUT $CCP4_SCR/toxd_step6.mtz \
DATOUT $CCP4_SCR/test.dat	\
<< 'END-phare'
TITLE Ref and Phasing
RESO 100.0 2.7
ANGLE 10
THRESHOLD 2.5 0.5
PRINT AVF AVE MON 1000 
CYCLES 10
LABI FP=FTOXD3 SIGFP=SIGFTOXD3 -
FPH1=FAU20 SIGFPH1=SIGFAU20  -
FPH2=FI100 SIGFPH2=SIGFI100 
LABO  ALLIN PHIB=PHI_step5 FOM=FOM_step5  

DERIV     Gold
DCYCLE PHASE ALL REFCYC ALL KBOV ALL
ISOE      746.82  548.14  407.42  365.98  379.86  329.45  346.91  345.99
ANOE      50 50 50 50 50 50 50 50 50 50 50 50 50 50
ATOM AU    0.177  0.104 -0.114  9.92 BFAC   25.000
ATREF X ALL Y ALL Z ALL OCC ALL 
ATOM  AU    0.218  0.138 -0.105  4.88 BFAC   25.000
ATREF X ALL Y ALL Z ALL OCC ALL 

DERIV     Iodine
DCYCLE PHASE ALL REFCYC ALL KBOV ALL
ISOE     100 100 100 100 100 100 100 100 100 100 
ANOE      50 50 50 50 50 50 50 50 50 50 50 50 50 50
ATOM  I     0.491  0.370  0.487  8.40 BFAC   25.000
ATREF X ALL Y ALL Z ALL OCC ALL 
END
END-phare

#----------------------------------------------------
#----------------------------------------------------
#       step6a_mlphare.com
#----------------------------------------------------
#----------------------------------------------------

# only necessary for certain spacegroups -- dummy step here
mtzutils			\
HKLIN1 $CCP4_SCR/toxd_sc.mtz \
HKLOUT $CCP4_SCR/scaled.mtz \
<< 'END-phare'
onefile
symm P212121
end
END-phare
#
#
mlphare			\
HKLIN $CCP4_SCR/scaled.mtz \
HKLOUT $CCP4_SCR/toxd_step6.mtz \
DATOUT $CCP4_SCR/test.dat	\
<< 'END-phare'
TITLE Ref and Phasing
RESO 100.0 2.7
ANGLE 10
THRESHOLD 2.5 0.5
PRINT AVF AVE MON 1000 
CYCLES 10
LABI FP=FTOXD3 SIGFP=SIGFTOXD3 -
FPH1=FAU20 SIGFPH1=SIGFAU20  -
FPH2=FI100 SIGFPH2=SIGFI100 
LABO  ALLIN PHIB=PHI_step5 FOM=FOM_step5  

DERIV     Gold
DCYCLE PHASE ALL REFCYC ALL KBOV ALL
ISOE      746.82  548.14  407.42  365.98  379.86  329.45  346.91  345.99
ANOE      50 50 50 50 50 50 50 50 50 50 50 50 50 50
ATOM AU    0.177  0.104 -0.114  9.92 BFAC   25.000
ATREF X ALL Y ALL Z ALL OCC ALL 
ATOM  AU    0.218  0.138 -0.105  4.88 BFAC   25.000
ATREF X ALL Y ALL Z ALL OCC ALL 

DERIV     Iodine
DCYCLE PHASE ALL REFCYC ALL KBOV ALL
ISOE     100 100 100 100 100 100 100 100 100 100 
ANOE      50 50 50 50 50 50 50 50 50 50 50 50 50 50
ATOM  I     0.491  0.370  0.487  8.40 BFAC   25.000
ATREF X ALL Y ALL Z ALL OCC ALL 
END
END-phare
#