#!/bin/csh   
#
#
#    Resolve ligand identification 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 identification script version 2.08.5 of 05-Jan-2005"
echo ""
#
#   Test many ligands for fit to map and identify the best
#
set hklin = resolve_diff.mtz
set labin = "FP=FP PHIB=PHIM FOM=FOMM"
set ligand_list_file = ../ligand_list_1.dat    # file containing list of ligand PDB files
set model = PARTIAL_MODEL.pdb
set resolve = /u1/terwill/resolve/work/resolve.linux  # where is resolve
set dmin = 2.6
set dmax = 200
#
#  some parameters you can set if you want to...
#
set n_indiv_tries_min = 10    # usually 0 to 10, but set up to 20 to try harder to find soln
set n_indiv_tries_max = 10    # usually 5-10, but set up to 20 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
#
echo "Resolution limits: $dmin $dmax A"
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 "File containing list of template PDB files is $ligand_list_file"
if  ( ! -f $ligand_list_file) then
  echo "Sorry, the file $ligand_list_file 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 $model"
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
@ n1 = $n_indiv_tries_max + 1
@ nn = $n_indiv_tries_min + 1
echo "Fitting will finish when number of tries is $n1 or when"
echo "all atoms are found and the number of tries is at least $nn"
#
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
#  make sure the headers exist...
@ header_lines = `cat $model |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:
endif
#
#  make sure all the ligand files in $ligand_list_file exist
set count = 0
set keep_going = 1
set lines = `cat $ligand_list_file |wc -l`
echo "Number of files in $ligand_list_file is $lines"
set count = $lines
while ($count > 0)
 set file = `tail -${count} $ligand_list_file |head -1 `
 #echo "File: $file"
 if ( ! -f $file ) then
  echo "Sorry the file $file from $ligand_list_file does not exist?"
  exit
 endif
@ count--
end
#
echo "Evaluating all ligands in $ligand_list_file now...and placing fitted"
echo "ligand xx in resolve_ligand_xx.pdb"
echo ""
echo " LIGAND    ATOMS PLACED      CC          SCORE        TEMPLATE FILE"
set count = $lines
set count_up = 0
while ($count > 0)
@ count_up++
 set file = `tail -${count} $ligand_list_file |head -1 `
 #echo "Evaluating ligand $file"
#
$resolve<<EOD   > resolve_fit_id_$count_up.log
resolution $dmax $dmin
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 $file
ligand_resno 1 # label ligand with 
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 ( -f ligand_fit.pdb )then
 if($count_up == "1" ) then
  cat ligand_fit.pdb  > all_ligands.pdb
 else
  cat ligand_fit.pdb|$grep_type ATOM  >> all_ligands.pdb
 endif
$resolve<<EOD>TEMP.DAT
hklin $hklin
labin $labin
resolution $dmin $dmax
model ligand_fit.pdb
evaluate_model
EOD
else
echo "" > TEMP.DAT
echo "" > resolve.ligand_scores
endif
#
set correl_line = `$grep_type 'region of model' TEMP.DAT|$grep_type Map`
#if (-f TEMP.DAT ) rm -f TEMP.DAT
if ( $#correl_line > 0 ) then
 set correl = $correl_line[$#correl_line-$#correl_line]
else
 set correl = 0
endif
#
set score_line = `cat resolve.ligand_scores |$grep_type MODIFIED`
set placed_line = `cat resolve.ligand_scores |$grep_type "ATOMS PLACED"`
set poss_line = `cat resolve.ligand_scores |$grep_type TOTAL`
if ( $#score_line > 0 ) then
 set score = $score_line[1-1]
else
 set score = 0
endif
if ( $#placed_line > 0 ) then
 set placed = $placed_line[1-1]:r
else
 set placed = 0
endif
if ( $#poss_line > 0 ) then
 set poss = $poss_line[1-1]:r
else
 set poss = 0
endif
cp ligand_fit.pdb resolve_ligand_${count_up}.pdb
echo "    $count_up        ${placed}/${poss}          $correl         $score     $file"
@ count --
end
#--------------------------------------------------------------------
#--------------------------------------------------------------------
echo ""
echo "Re-evaluating all ligands in $ligand_list_file now over common region"
echo ""
echo " LIGAND   CC     TEMPLATE FILE"
set count = $lines
set count_up = 0
head -10000 all_ligands.pdb > model_2.pdb   # limit to 10000 atoms
while ($count > 0)
@ count_up++
 set file = `tail -${count} $ligand_list_file |head -1 `
 #echo "Evaluating ligand $file"
#
$resolve<<EOD>TEMP.DAT
hklin $hklin
labin $labin
resolution $dmin $dmax
model resolve_ligand_${count_up}.pdb
model_2 model_2.pdb    ! evaluate in region of model_2
evaluate_model
EOD
#
set correl_line = `$grep_type 'region of model' TEMP.DAT|$grep_type Map`
#if (-f TEMP.DAT ) rm -f TEMP.DAT
if ( $#correl_line > 0 ) then
 set correl = $correl_line[$#correl_line-$#correl_line]
else
 set correl = 0
endif
#
echo "    $count_up      $correl  $file"
@ count --
end
#
#
#