#!/bin/csh   
#
#
#    Resolve ligand-fitting script 
#     T. Terwilliger 05-Jan-2005
#
#    This script requires resolve version 2.08 and ccp4 version 5.0.2
#
#    Before running this script, please get map coefficients with 
#    resolve_completion.com script(output is resolve_map.mtz optimized for ligand )
#    or else input map coefficients of your own in hklin and labin
#
#
setenv SOLVEDIR /usr/local/lib/solve/
setenv SOLVETMPDIR .
setenv SYMOP $SOLVEDIR/symop.lib
setenv SYMINFO $SOLVEDIR/syminfo.lib
setenv CCP4_OPEN UNKNOWN
unlimit
#
echo "Resolve ligand-fitting script version 2.08.8 of 05-Jan-2005"
echo ""
#   get N ligands in N passes through map, 
#    excluding previous as we go
#
set hklin = resolve_diff.mtz
set labin = "FP=FP PHIB=PHIM FOM=FOMM"
set number_of_ligands = 1    # number of copies of ligand to find
set ligand_template = ligand_template.pdb    # file with 1-10 copies of ligand
set model = PARTIAL_MODEL.pdb    #  model without ligand or ""
set resolve = /u1/terwill/resolve/work/resolve.linux  # where is resolve
set dmin = NONE   # NONE to use all data
set dmax = NONE   # NONE to use all data
#
#  some parameters you can set if you want to...
#
set n_indiv_tries_min = 10    # usually 0-10, but set up to 100 to try harder to find soln
set n_indiv_tries_max = 10    # usually 5-10, but set up to 100 to try harder to find soln
set n_group_search = 3    # usually 3, but set up to 10 to try harder to find soln
set search_dist = 10    # usually 10  A; always at least 5; smaller speeds up search
set no_local_search = ""    # usually "";  "no_local_search" to force complete search
set delta_phi_ligand = 40    # usually 40 degree increments; set lower to search more
#
if ( $dmin == NONE || $dmax == NONE ) then
 set resolution_line = ""
else
 set resolution_line = "resolution $dmin $dmax"
 echo "Resolution limits: $dmin $dmax A"
endif
#
echo "Input mtz file with map: $hklin"
if  (! -f $hklin ) then
  echo "Sorry, the file $hklin does not seem to exist?"
  exit
endif
echo "Labin line for $hklin is $labin"
echo "Number of copies of ligand to find: $number_of_ligands"
echo "Template PDB file containing 1 to 10 copies of ligand: $ligand_template"
if  ( ! -f $ligand_template) then
  echo "Sorry, the file $ligand_template does not seem to exist?"
  exit
endif
if ($model != "") then
 echo "The model $model will be used to mask out part of the map"
 if  ( ! -f $model) then
  echo "Sorry, the file $model does not seem to exist?"
  exit
 endif
set model_use = "model ALL.PDB"
else
set model_use = ""
endif
#
echo "Location of resolve: $resolve"
if  ( ! -f $resolve) then
  echo "Sorry, the program $resolve does not seem to exist?"
  exit
endif

echo "Number of groups (fragments) to search for with FFT: $n_group_search"
echo "FFT search will be in increments of $delta_phi_ligand degrees"
if ($no_local_search != "") then
 echo "Entire map will be searched"
endif
echo "Fitting will finish when number of tries is $n_indiv_tries_max or when"
echo "all atoms are found and the number of tries is at least $n_indiv_tries_min"
#
echo ""
#------------------------------------------------------------------------------
#--------------------------------------------------------------
#  figure out if this machine uses grep -a or just grep for text files:
echo "A" > test_a.dat
 set test_grep = `grep -a "a" test_a.dat >& tmp.dat`
if ( $status ) then
#  there was an error...do not use grep -a
 set grep_type = "grep"
else
 set grep_type = "grep -a"
endif
set test_grep = `$grep_type "A" test_a.dat`
if ( $#test_grep != 1 ) then
 echo "Sorry, unable to set the grep command on this system...giving up"
 exit
endif
rm test_a.dat
#--------------------------------------------------------------

#------------------------------------------------------------------------------
#
if ($model != "")then
cp $model ALL.PDB
cat ALL.PDB|$grep_type 'CRYST1' > resolve_ligand.pdb
cat ALL.PDB|$grep_type 'ORIGX[1-3]' >> resolve_ligand.pdb
cat ALL.PDB|$grep_type 'SCALE[1-3]' >> resolve_ligand.pdb
#  make sure the headers exist...
@ header_lines = `cat resolve_ligand.pdb |wc -l`
if ( $header_lines >= 4 ) goto ok1
 echo "The input PDB file $model needs to have at least 4 lines of headers with "
 echo "CRYST1 SCALE1 SCALE2 SCALE3"
 echo "Yours seems to have instead $header_lines of headers"
 cat resolve_ligand.pdb
 exit
ok1:
else
 if ( -f ALL.PDB ) rm -f ALL.PDB
 if ( -f resolve_ligand.pdb ) rm -f resolve_ligand.pdb
endif
#
set cycle = 1
while ($cycle <= $number_of_ligands)
#
echo "Searching for copy $cycle of ligand..."
$resolve<<EOD  
$resolution_line
hklin $hklin
labin $labin
ha_file NONE
no_build
delta_phi_ligand $delta_phi_ligand
$no_local_search
search_dist  $search_dist
$model_use   ! exclude region defined by model
ligand_file $ligand_template
ligand_resno $cycle   # label ligand with $cycle
n_indiv_tries_min  $n_indiv_tries_min
n_indiv_tries_max  $n_indiv_tries_max
n_group_search $n_group_search  ! how many groups to search for
! some more parameters that can be set:
! search_center  3 25 6 ! center the search here
! group_search  3     ! use this group in FFT search
! fit_phi_range -180 180  ! range of torsion angles to check
! fit_phi_inc   20    ! increment for torsion angle check
EOD
#
if ($cycle == 1 ) then
  cp ligand_map.mtz ligand_map_all.mtz
  cp ligand_map.map ligand_map_all.map
endif
#
cp ligand_fit.pdb ligand_$cycle.pdb
#
#  add new ligand to ALL.PDB and to resolve_ligand.pdb
#
if ( -f ALL.PDB )then
 cat  ligand_fit.pdb |$grep_type ATOM >> ALL.PDB
else
 cp  ligand_fit.pdb ALL.PDB
endif
if ( -f resolve_ligand.pdb ) then
 cat  ligand_fit.pdb |$grep_type ATOM >> resolve_ligand.pdb
else
 cp ligand_fit.pdb resolve_ligand.pdb
endif
#
@ cycle ++
end
#
cp ligand_map_all.mtz ligand_map.mtz
cp ligand_map_all.map ligand_map.map
#
echo "Evaluating fit of ligand to map...."
$resolve<<EOD>TEMP.DAT
hklin $hklin
labin $labin
$resolution_line
model resolve_ligand.pdb
evaluate_model
EOD
$grep_type 'region of model' TEMP.DAT|$grep_type Map
rm -f TEMP.DAT
echo ""
#
echo "All done...your fitted $number_of_ligands copies of ligand are in resolve_ligand.pdb"
echo "and coefficients for a map showing just the ligand regions are in ligand_map.mtz"