#!/bin/csh 
set main_version = 2
set minor_version = 11 
set edit_version = 3
#
#                   NOTE:  This script requires resolve version 2.08 or higher 
echo " "
echo "RESOLVE build script version $main_version.$minor_version.$edit_version "
echo " "
echo "Date: `date`"
echo "Working directory: `pwd`"
echo " "
#
#           RESOLVE iterative model-building script
#                       2-Feb-2007
#
#                      Tom Terwilliger
#              Los Alamos National Laboratory 
#
#
#    Input:  FP and SIGFP, experimental phases.
#            Input model  (optional)
#
#    Autobuilding:  Iteratively improve phases using statistical density
#            modification including solvent flattening, histogram-matching,
#            NCS, local pattern recognition, local fragment identification,
#            and model-building and refinement. Optionally include cycles of 
#            rebuilding (see below); accept if R decreases.
#
#    Rebuilding: Start with input model and calculate phases using
#            density modification including the model information. Optionally
#            calculate and use prime-and-switch composite omit map. 
#            Build new model into the map...iterate
#
#    Model evaluation: as rebuilding, but just evaluate the model using the
#            prime-and-switch composite omit map
#
#    Output: Optimized phases in resolve_best.mtz 
#            Partially-refined model in resolve_best.pdb
#            Log file for best cycle in resolve_best.log
#            CCP4-style map file in resolve_best.map
#
#    Major changes:  Version 2.08 allows you to build
#            outside of a region defined by a model.
#            Version 9 allows definition of a mask with a PDB file
#
#======================================================================
#     EDIT THE NEXT FEW SETS OF LINES TO MATCH YOUR DATA AND SYSTEM
#
#     PLEASE NOTE: each of the " = " below must have a SPACE on either
#       side (hklin = solve.mtz  NOT hklin=solve.mtz )
#
#======================================================================
#
#  Anything special for your location (SOLVEDIR etc);
#
setenv SOLVEDIR $PHENIX_SOLVEDIR
setenv TMPDIR .                            # large scratch files go here
#
# location of resolve and resolve libraries;  place for large scratch files
#
set resolve =  $PHENIX_SOLVEBIN/resolve
set resolve_pattern =  $PHENIX_SOLVEBIN/resolve_pattern
#
#   location of refmac5
#
set refmac5 = refmac5 # NONE if you do not have refmac5
#
#   Input mtz file with at least FP and SIGFP and PHIB
#
set hklin = perfect.mtz
set labin = 'FP=FP SIGFP=SIGFP '
set labin_cont = ''
#
#   Optional input model to be rebuilt  
#
set pdb_in =   start.pdb #  NONE if no model to rebuild
#
#   resolution, solvent_content, sequence, heavy-atom file
#
set dmin =  3.0                    #   high-res limit  
set dmax = 200                     #   low-res limit
set solvent_content = 0.40         #   solvent content (required)
set seq_file = seq.dat           #   amino acid sequence file (up to 80 char/line)
set ha_file =  NONE                  #     e.g., ha.pdb or NONE; used to get NCS
#
set refine_model =  YES        #  say NO if resolution > 3 A usually
set compare_file =  start.pdb #   NONE if no comparison file
#
#    Just build outside the region defined by input model? (Only use if input model has
#    side-chains. Not for a backbone model)
#
set mask_pdb_file = NONE   # PDB file specifying mask for density modification
#
#   criteria for ending..
#
set n_cycle_build_max = 20   # try 20
set n_cycle_build_min = 3    # try 3
set unacceptable_r = 0.4 # try 0.4 keep going if R > unacceptable_r
set n_cycle_rebuild_max = 30  # try 30 for build, 50 for rebuild
set n_cycle_rebuild_min = 20  # try 20 (Rebuild skipped if R < unacceptable_r)
set n_cycle_rebuild_omit = 10 #  every n'th cycle do omit  if 
                              # omit_on_rebuild is YES (starting with nth cycle)
#
#   misc optional parameters...
#
set use_resolve_fragments = AUTO  #  AUTO or YES or NO  (AUTO usually)
set use_resolve_pattern   =  AUTO #  AUTO or YES or NO  (AUTO usually)
#
set superquick_build =     NO #  usually NO. Use YES to get a quick answer
#
set clean_up = YES          # YES usually
set extra_verbose = NO      # only for debugging
set omit_on_rebuild = NO    # NO usually
set quick_as_possible = NO  # NO usually
set n_cycle_image_min = 3   # 3 usually use image info for at least these cycles
set cycle_start = 1         # 1 usually
set evaluate_only =  NO   # NO usually; YES ok; NEVER skip
set use_prob_mask = NO    # NO usually; yes for version 2.0 probability mask
set simple_scaling = NO   # YES to use SCALE TYPE SIMPLE in refmac5
#
#==============================================================================
#                Normally no need to edit below here...
#==============================================================================
#
# Last edits: 091206: added SIGFP=SIGFP to all statements with labin FP=FP 
#   containing also HLA. Also added simple_scaling option.
#
set skip =  NO
#
set main_version_minimum = 2
set minor_version_minimum = 08
set build_outside_model = NO  # Does not always work. refmac5 fails if resolve writes a residue
                              # with same residue number as a residue in pdb_in
#
setenv SYMOP $SOLVEDIR/symop.lib
setenv SYMINFO $SOLVEDIR/syminfo.lib
setenv CCP4_OPEN UNKNOWN
unlimit
limit coredumpsize 0
#
if ( $extra_verbose == YES ) then
 set ip = 1
else
 set ip = 0
endif
#--------------------------------------------------------------
#  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 -f test_a.dat
#--------------------------------------------------------------
#  check for all the library files we will need etc
#
if ( $resolve == NONE || $resolve_pattern == NONE ) then
 echo "Sorry, you need to define resolve and resolve_pattern..."
 exit
endif
foreach program ($resolve $resolve_pattern $refmac5)
if ( $ip ) echo "Checking for the program $program ..."
if ( -f $program ) then
 if ( $ip ) echo "OK"
 goto ok
endif
if (  $program == NONE ) then
 if ( $ip ) echo "This program is not used"
 goto ok
endif
which $program >& tmp.dat
set test = `head -1 tmp.dat`
if ( $#test != 1 ) goto bad
if ( -f $test ) then
 if ( $ip ) echo "OK"
 goto ok
endif
bad:
 echo "Sorry, the program $program does not exist?"
 echo "Please check its definition in this script..."
 exit
endif
ok:
end
#
#
foreach lib_file ( symop.lib segments/rho.list \
 segments/fragment_1_0.pdb segments/fragment_2_0.pdb  \
 segments/segment_001_1.pdb segments/segment_001_2.pdb \
 segments/side_chains_best.dat segments/side_chains_rho.dat \
 segments/gap_library.dat \
 segments/template_30_mean_1.ezd segments/template_30_mask_1.ezd \
 segments/template_30_corr_1.ezd segments/template_30_mean_2.ezd \
 segments/template_30_mask_2.ezd segments/template_30_corr_2.ezd \
 hist.sol_91 hist.sol_92 \
 patterns/hist_cc_0.8.dat patterns/index_95.dat \
 patterns/patterns.info patterns/box001.dat )
#
if ( $ip ) echo "Checking for library file $SOLVEDIR/$lib_file..."
if ( ! -f $SOLVEDIR/$lib_file) then
 echo "Sorry, the library file $lib_file does not exist?"
 echo "Please check the definition of SOLVEDIR in your script: $SOLVEDIR"
 exit
endif
end
#--------------------------------------------------------------
#
if ( $quick_as_possible == YES ) then
 echo "Trying to do everything as quickly as possible..."
 echo "Setting superquick_build, no refinement, no rebuild, no pattern ID,"
 echo "no fragment ID, minimum cycles"
 set refine_model = NO
 set superquick_build = YES
 set use_resolve_pattern = NO
 set use_resolve_fragments = NO
endif
if ( $evaluate_only == YES ) then
 echo "Just evaluating the model $pdb_in"
 set n_cycle_rebuild_max = 0
 set n_cycle_rebuild_min = 0
endif
if ( ( $pdb_in != NONE ) && ( $evaluate_only == YES ) && ( ! -f $pdb_in ) ) then
 echo "Sorry cannot find the model $pdb_in ... required for evaluate_only"
 exit
endif
if ( $build_outside_model == NO || $pdb_in == NONE )then
  # do nothing
else
 echo "Building only around model $pdb_in (excluding region in existing model)"
endif
if ( -f HEADER.pdb ) rm -f HEADER.pdb
if ($pdb_in != NONE )then
$grep_type CRYST1 $pdb_in >>HEADER.pdb
$grep_type ORIGX1 $pdb_in >>HEADER.pdb
$grep_type ORIGX2 $pdb_in >>HEADER.pdb
$grep_type ORIGX3 $pdb_in >>HEADER.pdb
$grep_type SCALE1 $pdb_in >>HEADER.pdb
$grep_type SCALE2 $pdb_in >>HEADER.pdb
$grep_type SCALE3 $pdb_in >>HEADER.pdb
#  make sure the headers exist...
@ header_lines = `cat HEADER.pdb |wc -l`
if ( $header_lines >= 4 ) goto ok1
 echo "The input PDB file $pdb_in 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 HEADER.pdb
 exit
ok1:
echo "REMARK RESOLVE MODEL" >> HEADER.pdb
endif
#
if ( $n_cycle_build_min < 1 ) then
   echo "Sorry, please set n_cycle_build_min to at least 1"
   exit
endif
if ( $n_cycle_build_max < 1 ) then
   echo "Sorry, please set n_cycle_build_max to at least 1"
   exit
endif
if ( $n_cycle_rebuild_min < 0 ) then
   set n_cycle_rebuild_min = 0
endif
if ( $n_cycle_rebuild_max < 0 ) then
    set n_cycle_rebuild_max = 0
endif
#
#
set past_cycle = $n_cycle_build_min  
set n_cycle = $n_cycle_build_max
set n_cycle_rebuild = $n_cycle_rebuild_max
#   (set n_cycle to n_cycle_rebuild_max later)
#
set minimum_coverage = 0
#
#
#
# check for resolve and resolve_pattern
if ( -f resolve.ok ) rm -f resolve.ok
if ( -f resolve.version ) rm -f resolve.version
$resolve <<EOD  >& resolve.log_0
 quit 
EOD
if ( ! -f resolve.ok ) then
 echo "Sorry please check the location of resolve..."
 echo "    is it really $resolve?"
 exit
endif
#  make sure this sort of worked...
set test_ok = `cat resolve.ok`
set test_ok_mem = $test_ok[$#test_ok-$#test_ok]
if ( $test_ok_mem != "ALLOCATED" ) then
 echo "Sorry, resolve was not able to run properly at all"
 echo "The end of the resolve log file says..."
 echo "---------------------------------------------------------"
 tail -12 resolve.log_0
 echo "---------------------------------------------------------"
 exit
endif
if ( ! -f resolve.version ) then
 echo "Sorry, this script requires version $main_version_minimum.$minor_version_minimum or higher of resolve"
 exit
endif
@ version = `cat resolve.version|head -1`
set version_minimum = "$main_version_minimum$minor_version_minimum" 
if ( $version < $version_minimum )then
 echo "Sorry, this script requires version $main_version_minimum.$minor_version_minimum or higher of resolve"
 exit
else
 if ( $ip ) echo "Running version $version of resolve"
endif 
if ( -f resolve_pattern.ok ) rm -f resolve_pattern.ok
$resolve_pattern <<EOD >& /dev/null
end 
EOD
if ( ! -f resolve_pattern.ok ) then
 echo "Sorry please check the location of resolve_pattern... "
 echo "   is it really $resolve_pattern?"
 exit
endif
#
# make sure we have these files...
#
if ( -f $hklin ) then
   echo "Data are in the mtz file $hklin" 
else
   echo "Sorry, cannot find the mtz file $hklin" 
   exit
endif
echo "LABIN information: $labin $labin_cont"
echo " "
#
if ( $seq_file != NONE &&  -f $seq_file  ) then
  echo "Sequence information from ${seq_file}"
else if ( $seq_file != NONE ) then
  echo "Sorry, cannot find the seq_file file $seq_file"
  exit
endif
#
if ( ( $ha_file != NONE ) && ( -f $ha_file  ) ) then
  echo "Atom information for NCS from $ha_file"
else if ( $ha_file != NONE ) then
  echo "Sorry, cannot find the ha_file file $ha_file"
  exit
endif
#
if ( ( $pdb_in != NONE ) && ( ! -f $pdb_in ) ) then
 echo "Sorry, cannot find your pdb_in file $pdb_in"
 exit 
else if ( ( $cycle_start == 1 ) && ( $pdb_in != NONE ) ) then
 set use_resolve_fragments = NO
 set use_resolve_pattern   =  NO
endif
#
if ( $superquick_build == YES ) then
 set superquick_build = "superquick_build"
else
 set superquick_build = " "
endif
#
@ grep_phib = `echo $labin $labin_cont | $grep_type 'PHIB='|wc -m`
@ grep_fp = `echo $labin $labin_cont | $grep_type 'FP='|wc -m`
#
if ( $grep_fp == 0 ) then
 echo "Sorry, FP is always required in the labin line..."
 exit
endif
if ( $grep_phib == 0 && $pdb_in == NONE ) then
 echo "Sorry if you have no pdb_in file you need phases (PHIB=...)"
 exit
else if ( $pdb_in != NONE ) then
 echo "Using  phase probabilities from $hklin in rebuilding"
endif
if ( $pdb_in == NONE ) then
 set scale = 1.0
else if ( $grep_phib != 0 ) then
 set scale = 1.0
else
 set scale = 0.2
endif
#
echo " "
echo "solvent content ${solvent_content}"
echo "resolution ${dmin} ${dmax}"
if ( $refine_model == NO ) then
 echo "Model will not be refined"
else
 echo "Model will be refined"
endif
echo " "
if ( $compare_file != NONE ) then
 echo "Coordinates for comparison from $compare_file"
endif
if ( $compare_file != NONE && ! ( -f $compare_file ) ) then
  echo "Sorry, cannot find the comparison PDB file $compare_file"
  exit
endif
#
set blank = " "
#
if ( $pdb_in != NONE ) then
  set combine_resolve_fragments = $blank
  echo "Not using fragment ID"
else if ( $use_resolve_fragments == YES ) then
  set combine_resolve_fragments = "combine_map fragments.map"
  echo "Using fragment ID"
else if ( $use_resolve_fragments == AUTO ) then
  set combine_resolve_fragments = "combine_map fragments.map"
  echo "Testing using fragment ID"
else 
  set combine_resolve_fragments = $blank
  echo "Not using fragment ID"
endif 
if ( $pdb_in != NONE ) then
  set combine_resolve_pattern = $blank
  echo "Not using resolve_pattern"
else if ( $use_resolve_pattern == YES) then
  set combine_resolve_pattern = "combine_map pattern.map"
  echo "Using resolve_pattern"
else if ( $use_resolve_pattern == AUTO) then
  set combine_resolve_pattern = "combine_map pattern.map"
  echo "Testing using resolve_pattern"
else
  set combine_resolve_pattern = $blank
  echo "Not using resolve_pattern"
endif
#  make n_unacceptable_r  have 5 digits exactly
set n_unacceptable_r = $unacceptable_r:e
if ( $n_unacceptable_r < 1 ) then
 set n_unacceptable_r = 10000
endif
while ( $n_unacceptable_r < 10000 ) 
 @ n_unacceptable_r = $n_unacceptable_r * 10
end
if ( ( $pdb_in != NONE ) && ( $evaluate_only != YES ) ) then
 echo "If R-factor < 0.$n_unacceptable_r then no rebuilding will be done"
 echo "If R-factor >= 0.$n_unacceptable_r then maximum cycles will be done"
endif
if ($mask_pdb_file != NONE  && ! -f $mask_pdb_file)then
  echo "Sorry, could not find the mask pdb file $mask_pdb_file"
  exit
else if ($mask_pdb_file != NONE )then
  set mask_line = "model $mask_pdb_file"
  set use_mask_line = "use_model_mask"
  echo "Using $mask_pdb_file to define mask for density modification"
else if ($use_prob_mask == YES )then
  set mask_line = ""
  set use_mask_line = "use_prob"
  echo "Using probability-based mask (version 2.0)"
else
  set mask_line = ""
  set use_mask_line = ""
endif
if ($simple_scaling == YES)then
  set simple_scaling_line = "SCALe TYPE SIMPLE"
  echo "Using simple scaling in refmac5"
else
  set simple_scaling_line = "SCALe TYPE BULK"
  echo "Using BULK scaling in refmac5"
endif
#
#
#           SETUP:
#
# =======================================================================
# =======================================================================
# hklin_image will be our combined phase information 
#
set hklin_image = start_image.mtz
set labin_image = 'FP=FP PHIB=PHIM FOM=FOMM HLA=HLAM HLB=HLBM HLC=HLCM HLD=HLDM'
set labin_image_cont = 'FreeR_flag=FreeR_flag'
#
#  hklin_exp is the file we will use to evaluate phase quality. 
set hklin_exp = exp.mtz
set labin_exp = 'FP=FP PHIB=PHIM FOM=FOMM HLA=HLAM HLB=HLBM HLC=HLCM HLD=HLDM'
set labin_exp_cont = 'SIGFP=SIGFP FreeR_flag=FreeR_flag'
#
set build_score_best_overall = -1000
set build_score_best_overall_cycle = -1000
set r_best_overall = 99999
#
# =======================================================================
set cycle = $cycle_start
set cycle_prev = $cycle
@ cycle_prev --
# =======================================================================
if ( $cycle_start == 1 ) then
 goto cycle_one
endif
echo "Starting at cycle $cycle_start"
#
cycle_one:
# =======================================================================
set jump = NO
if ( $jump == YES ) then
set best_cc_value = 0000
set read_mask = " " 
echo "JUMPING TO MIDDLE OF FILE AT CYCLE $cycle_start"
endif
if ( $jump == YES ) goto jump_point
# =======================================================================
#
while ( $cycle <= $n_cycle )
#
set cycle_prev = $cycle
@   cycle_prev --
#
 if ( $ip ) echo " "
if ( $ip ) echo "Cycle $cycle -------------------------------------------"
#
#
if ( -f resolve.log_$cycle ) then
 if ( $ip ) echo "removing previously-existing log file resolve.log_$cycle"
 rm -f resolve.log_$cycle
endif
set build_score_best = -1000
set placed_best = 0
set built_best = 0
if ( -f work_model.pdb ) then
  rm -f work_model.pdb
endif
#
if ( $cycle > 1 ) then
 goto cycle_point
endif
#
# =======================================================================
#
#  First cycle:
#
#   remove all the files from previous runs..
if ( -f pattern.map ) then
 if ( $ip ) echo "removing existing pattern.map"
 rm -f pattern.map
endif
if ( -f fragments.map ) then
 if ( $ip ) echo "removing existing fragments.map"
 rm -f fragments.map
endif
if ( -f combine.map ) then
 if ( $ip ) echo "removing existing combine.map"
 rm -f combine.map
endif
if ( -f current_model.pdb ) then
  if ( $ip ) echo "Removing current_model.pdb"
  rm -f current_model.pdb
endif
if ( -f resolve_peaks.dat ) then
  if ( $ip ) echo "removing resolve_peaks.dat"
  rm -f resolve_peaks.dat
endif
if ( -f resolve_peaks.use) then
  if ( $ip ) echo "removing resolve_peaks.use"
  rm -f resolve_peaks.use
endif
if ( -f fragments.dat ) then
  if ( $ip ) echo "removing fragments.dat"
  rm -f fragments.dat
endif
foreach best_file (resolve_best.pdb resolve_best.log resolve_best.mtz \
    refmac_best.log )
if ( -f $best_file ) rm -f $best_file
end
#
#==============================================================================
echo "Copying phase information from $hklin to $hklin_exp"
${resolve}<<EOD >>resolve.log_0
! Copying over $hklin to $hklin_exp
hklin $hklin
labin $labin
labin $labin_cont
hklout $hklin_exp
solvent_content $solvent_content
resolution $dmin $dmax
res_start $dmin
mask_cycles 1
minor_cycles 0
ha_file NONE
no_build
EOD
#
#  make sure this sort of worked...
if ( `cat resolve.ok` != 'OK' ) then
 echo "Sorry, resolve was not able to finish even the first cycle"
 echo "The end of the resolve log file says..."
 echo "---------------------------------------------------------"
 tail -10 resolve.log_0
 echo "---------------------------------------------------------"
 exit
endif
#
#  for model start we go right on to rebuilding...cycle 2
#
if (  $pdb_in != NONE )  then
 echo " "
 echo "Going right on to rebuilding/model evaluation now..."
 break
endif
#
#  First cycle just do std density modification...2x with
#   2 sets of parameters
#
set best_dm_type = "quick"
#
if ( $quick_as_possible == YES ) then
 set types = "quick"
else
 set types = "quick"
endif
foreach build_type ( $types )
if ( $build_type == "quick" ) then
 set build_1 = "mask_cycles 3"
 set build_2 = "minor_cycles 1"
else
 set build_1 = $blank
 set build_2 = $blank
endif
echo "Testing $build_type density modification..."
#
if ( -f resolve.built ) rm -f resolve.built
$resolve <<EOD >>resolve.log_$cycle
!  Resolve_build script
!  Cycle $cycle...testing $build_type density modification
resolution $dmin $dmax
hklin $hklin
labin $labin
labin $labin_cont
solvent_content $solvent_content
ha_file $ha_file
seq_file $seq_file
$superquick_build
compare_file $compare_file
hklout dummy.mtz
$build_1
$build_2
$mask_line
$use_mask_line
EOD
#
#  make sure this sort of worked...
if ( `cat resolve.ok` != 'OK' ) then
 echo "Sorry, resolve was not able to finish even the first cycle"
 echo "The end of the resolve log file says..."
 echo "---------------------------------------------------------"
 tail -10 resolve.log_0
 echo "---------------------------------------------------------"
 exit
endif
# 
#   score this build:
#
if ( -f resolve.built ) then
    set built = `cat resolve.built|head -1` 
    set placed = `cat resolve.built | head -2|tail -1` 
    set build_score = `cat resolve.built | head -4|tail -1` 
 echo "Score for $build_type : $build_score   built: $built  placed: $placed"
else
 set build_score = 0
 set built = 0
 set placed = 0
endif
#
if ( $build_score > $build_score_best ) then
 set best_dm_type = $build_type
 set build_score_best = $build_score
 set built_best = $built
 set placed_best = $placed
 cp dummy.mtz resolve.mtz_$cycle
 cp resolve.mtz_$cycle omit.mtz
 cp resolve_peaks.dat resolve_peaks.use
 cp resolve.pdb  work_model.pdb 
endif
#
end
echo "Best density modification obtained with $best_dm_type approach"
echo $best_dm_type > resolve.type
#
echo " "
echo "Starting building at cycle $cycle"
echo " "
echo "-------------------------------------------------------------"
echo "         RESIDUES  SIDE-CHAINS"
echo " CYCLE     BUILT      PLACED        SCORE      R       Free R   BEST"
#
goto finish
#
#==============================================================================
cycle_point:
#
# Figure out if we need to try alternate ncs:
if ( -f resolve.ncs_used ) then
 set ncs_used = `cat resolve.ncs_used|head -1`
else
 set ncs_used = NO
endif
if ( ( $ncs_used == NO ) &&  ( -f ncs_test.pdb ) ) then
 set ha_file = ncs_test.pdb
 echo "Trying NCS from $ha_file based on last model built..."
endif
#
#  for first few cycles always use image info if available...
if ( ( $cycle <= $n_cycle_image_min) && ( -f combine.map ) ) then
 goto use_image
endif
#
# density-modify with model information only:
#
if ( -f resolve.built ) rm -f resolve.built
${resolve}<<EOD>>resolve.log_$cycle
! map with model information only cycle $cycle
hklin $hklin_exp
labin $labin_exp
labin $labin_exp_cont
hklout resolve.mtz_$cycle
solvent_content $solvent_content
ha_file $ha_file
seq_file $seq_file
$superquick_build
compare_file $compare_file
composite_all
composite_pdb refine.pdb_
composite_pdb_first 1 
composite_pdb_last $cycle_prev
pdb_in refine.pdb_$cycle_prev
image
mask_cycles 1
nohl
$mask_line
$use_mask_line
EOD
#
if ( -f resolve.built ) then
    set built = `cat resolve.built|head -1` 
    set placed = `cat resolve.built | head -2|tail -1` 
    set build_score = `cat resolve.built | head -4|tail -1` 
 if ( $ip ) echo "Build score with model but no image info: $build_score   built: $built  placed: $placed"
else
 set build_score = 0
 set built = 0
 set placed = 0
 echo "REMARK No model built" > resolve.pdb
 echo "REMARK No model built" 
endif
#
set best_method = "no_image"
set build_score_best = $build_score
set built_best = $built
set placed_best = $placed
cp resolve_peaks.dat resolve_peaks.use
cp resolve.pdb  work_model.pdb 
#
#  Include model information with image info and build a new model:
#
if ( ! -f combine.map ) then
 goto skip_image
endif
#
use_image:
#
#
 if ( $ip ) echo "Adding in information from combine.map to $hklin_image"
#
${resolve}<<EOD >>resolve.log_$cycle
! applying combine info cycle $cycle
resolution $dmin $dmax
hklin $hklin_exp
labin $labin_exp
labin $labin_exp_cont
solvent_content $solvent_content
pattern_phase
hklout $hklin_image
cc_map_file combine.map
ha_file NONE
no_build
EOD
#
if ( -f resolve.built ) rm -f resolve.built
${resolve}<<EOD>>resolve.log_$cycle
! Get map with model and image information cycle $cycle
hklin $hklin_image
labin $labin_image
labin $labin_image_cont
hklout dummy.mtz
solvent_content $solvent_content
ha_file $ha_file
seq_file $seq_file
$superquick_build
compare_file $compare_file
composite_all
composite_pdb refine.pdb_
composite_pdb_first 1
composite_pdb_last $cycle_prev
pdb_in refine.pdb_$cycle_prev
image
mask_cycles 1
nohl
$mask_line
$use_mask_line
EOD
#
#   score this build:
#
if ( -f resolve.built ) then
    set built = `cat resolve.built|head -1` 
    set placed = `cat resolve.built | head -2|tail -1` 
    set build_score = `cat resolve.built | head -4|tail -1` 
 if ( $ip ) echo "Score for build using image info  : $build_score   built: $built  placed: $placed"
else
 set build_score = 0
 set built = 0
 set placed = 0
endif
if ( $build_score > $build_score_best ) then
 set best_method = "image"
 set build_score_best = $build_score
 set built_best = $built
 set placed_best = $placed
 cp resolve_peaks.dat resolve_peaks.use
 cp resolve.pdb  work_model.pdb 
 cp dummy.mtz resolve.mtz_$cycle
 if ( $ip ) echo "Including image information"
else if ( $use_resolve_fragments == AUTO ) then
 echo "No longer including image information"
  set combine_resolve_pattern = $blank
  set combine_resolve_fragments = $blank
  rm -f combine.map
  goto skip_image
else
 if ( $ip ) echo "Not using image information on this cycle"
endif
#
# Also get the omit map if we are using images... 
#
if ( $ip ) echo "Getting composite ps omit map..."
if ( -f omit.mtz ) then
 rm -f omit.mtz
endif
$resolve <<EOD >>resolve.log_$cycle
!  Cycle $cycle...composite prime-and-switch omit density modification
hklin $hklin_image
labin $labin_image
labin $labin_image_cont
hklout omit.mtz
solvent_content $solvent_content
no_build
ha_file $ha_file
composite_all
composite_pdb refine.pdb_
composite_pdb_first 1
composite_pdb_last $cycle_prev
add_mask
image
omit
mask_cycles 1
nohl
$mask_line
$use_mask_line
EOD
#
skip_image:
#
finish:
#
#  Start of refine-extend cycles
#
foreach extend_cycle ( 0 1 2)
#
if ( $ip ) echo "Extend cycle $extend_cycle"
if( $extend_cycle > 0 ) then
 if ( $ip ) echo "Extending model"
#
if ( -f resolve.built ) rm -f resolve.built
${resolve}<<EOD >> resolve.log_$cycle
hklin resolve.mtz_$cycle
labin FP=FP PHIB=PHIM FOM=FOMM 
hklout dummy.mtz
solvent_content $solvent_content
resolution ${dmin} ${dmax}
extend_only
ha_file $ha_file
seq_file $seq_file
compare_file $compare_file
pdb_in resolve.pdb
no_trim
nohl
EOD
endif
if ( ( $extend_cycle > 0 ) && ( -f resolve.built ) ) then
    set built = `cat resolve.built|head -1` 
    set placed = `cat resolve.built | head -2|tail -1` 
    set build_score = `cat resolve.built | head -4|tail -1` 
 if ( $ip ) echo "Score for extension: $build_score   built: $built  placed: $placed"
else if ( $extend_cycle > 0 ) then
 set build_score = 0
 set built = 0
 set placed = 0
endif
if ( ( $extend_cycle > 0 ) && ( $build_score > $build_score_best ) ) then
 set build_score_best = $build_score
 set built_best = $built
 set placed_best = $placed
 cp resolve.pdb work_model.pdb
else if ( $extend_cycle > 0 ) then
 if ( $ip ) echo "Extension did not improve score...going on"
#  break to end of extend cycles..
 break
endif
#
#
#  decide if we refine model this cycle..
#
if ( -f resolve.coverage ) then
 set coverage1 = `cat resolve.coverage|head -3|tail -1`
 set coverage = 0
 if ( $#coverage1> 0 ) set coverage= $coverage1[$#coverage1-$#coverage1]:e
else
 set coverage = 0
endif
if ( $refine_model == NO ) then
 set refine_it = NO
else if ( $build_score_best < 2 ) then
 set refine_it = NO
else if (  $coverage < $minimum_coverage )then
 set refine_it = NO
else if ( ! -f work_model.pdb ) then
 set refine_it = NO
 echo "REMARK No model built" > work_model.pdb
 echo "REMARK No model built extend cycle $extend_cycle" 
else
 set refine_it = YES
endif
if ( $refine_it == YES ) then
 set ncyc = 20
else
 set ncyc = 0
endif
if ( ( $refmac5 != NONE ) && ( $build_score_best > 1 ) ) then
if ( $ip ) echo "Refining model $ncyc cycles and placing it in refine.pdb_$cycle"
#
cat work_model.pdb | $grep_type -v HETATM > resolve.pdb
#
# --------refine with refine.--------------
#
#  NOTE: refine against "best" map = resolve.mtz_$cycle
#
$refmac5 xyzin resolve.pdb hklin resolve.mtz_$cycle xyzout refine.pdb  \
PROTOUT $CCP4_SCR/protout.dat \
PROTCOUNTS $CCP4_SCR/counts.dat \
PROTSCR $CCP4_SCR/counts.scr \
<<EOD > refmac.log_${cycle}
MAKE HYDRogens No
MAKE CHECK 0
LABI FP=FP SIGFP=SIGFP HLA=HLAM HLB=HLBM HLC=HLCM HLD=HLDM FREE=FreeR_flag
REFI TYPE RESTrained RESOlution ${dmax} ${dmin}
#   Maximum likelihood refinement
REFI RESI MLKF
#   Overall B value refinement
REFI BREF OVER      ! Refine overall B-values
#   Tight restraints  Lowered further 051403
WEIG MATR 0.05
#   Scale down shifts at every cycle by factor 0.5
DAMP 0.5 0.5
#   No Blur (scale down reliability of) phases
REFI PHASed  BBLUrring 00.0  SCBLurring 1.
#  Babinet's bulk solvent scale parameters. But fix B value
$simple_scaling_line
SCALe LSSCale FIXBulk BBULk 200.0
#   Anisotropic scaling. It is default
SCALe LSSCale ANISotropic
#   Fix Babinet's bulkd solvent for sigmaA
SCAL MLSC FIXBulk BBULk 200.0 SCBUlk -0.05
NCYC  $ncyc
MONI MEDI
BINS 10
EOD
#----------------------------------------------------------
#
cp refine.pdb refine.pdb_$cycle
else if (  $build_score_best > 1  ) then
#   no refmac5...just copy over
 cat work_model.pdb | $grep_type -v HETATM > resolve.pdb
 cp resolve.pdb refine.pdb
 cp refine.pdb refine.pdb_$cycle
else
if ( $ip ) echo "No model built this try"
 echo "REMARK  No model built cycle $cycle " > refine.pdb_$cycle
 echo "REMARK  No model built cycle $cycle " 
endif
#
end_of_extend:
end
#
#
#  Get R and freeR value on refine.pdb_$cycle:
#
set freer_line = `$grep_type 'FREE R VALUE                     :' refine.pdb_$cycle|tail -1`
set r_line = `$grep_type 'R VALUE            (WORKING SET) :' refine.pdb_$cycle|tail -1`
if ( $#freer_line > 0 ) then
 set freer = $freer_line[$#freer_line-$#freer_line]:e
else
 set freer = 99999
endif
if ( $#r_line > 0 ) then
 set r = $r_line[$#r_line-$#r_line]:e
else
 set r = 99999
endif
if ( $ncyc == 0 )then
 set freer = $r
endif
#
if ( $ip ) echo " R/freeR for this model: 0.$r/0.$freer"
#
if ( $cycle < 10 ) then
 set b1 = " "
else
 set b1 = ""
endif
if ( $built_best < 100 ) then
 set b2 = " "
else
 set b2 = ""
endif
if ( $placed_best >= 100 ) then
 set b3 = ""
else if ( $placed_best >= 10 )then
 set b3 = " "
else
 set b3 = "  "
endif
if ( $build_score_best < 100 )then
 set b4 = " "
else
 set b4 = ""
endif
# Decide how to score. If R < unacceptable_r or r_best_overall < unacceptable_r
#   then use R. Otherwise use score
#
if ( ( $r < $n_unacceptable_r ) && ( $r < $r_best_overall ) ) then
 set this_is_best = YES
else if ( ( $r_best_overall >= $n_unacceptable_r ) && \
   ( $build_score_best > $build_score_best_overall ) ) then
 set this_is_best = YES
else
 set this_is_best = NO
endif
if ( $this_is_best == YES ) then
 if ( $ip ) echo "This is the best model overall so far. Previous best had score of $build_score_best_overall"
 if ( $ip ) echo "Copied refine.pdb_$cycle to resolve_best.pdb"
 if ( $ip ) echo "Copied resolve.log_$cycle resolve_best.log"
 if ( $ip ) echo "Copied resolve.mtz_$cycle to resolve_best.mtz"
 @ build_score_best_overall = $build_score_best
 @ build_score_best_overall_cycle = $cycle
 @ r_best_overall = $r
 cp refine.pdb_$cycle resolve_best.pdb
 cp resolve.log_$cycle resolve_best.log
 cp resolve.mtz_$cycle resolve_best.mtz
 if( -f refmac.log_$cycle ) cp refmac.log_$cycle refmac_best.log
 echo "  ${b1}$cycle       ${b2}$built_best         ${b3}$placed_best          ${b4}$build_score_best     0.$r    0.$freer    *"
else
 echo "  ${b1}$cycle       ${b2}$built_best         ${b3}$placed_best          ${b4}$build_score_best     0.$r    0.$freer    "
endif
#
#
#====================================================================
#====================================================================
#
#
set have_map = NO
if ( ( "$combine_resolve_fragments" != $blank ) && ( -f resolve_peaks.dat )  )  then
if ( $ip ) echo "Carrying out fragment ID based on resolve.mtz_$cycle"
cp resolve_peaks.use resolve_peaks.dat
#
$resolve <<EOD >> resolve.log_$cycle
hklin resolve.mtz_$cycle
labin FP=FP PHIB=PHIM FOM=FOMM
solvent_content $solvent_content
nohl
build_image
noget_peaks      ! Don't need peaks because we just got them!
mask_cycles 0
minor_cycles 0
superquick_build
ha_file NONE
$mask_line
$use_mask_line
EOD
cp dump.map fragments.map
set have_map = YES
endif
#
if ( "$combine_resolve_pattern" != $blank ) then
if ( $ip ) echo "Carrying out pattern ID based on omit.mtz"
$resolve_pattern<<EOD >> resolve.log_$cycle
hklin omit.mtz
labin FP=FP PHIB=PHIM FOM=FOMM
recover_image
cc_map_file dump.map
EOD
cp dump.map pattern.map
set have_map = YES
endif
#
#  Decide which to combine..
#
if ( $have_map == YES ) then
if ( $ip ) echo "Combining all image information into combine.map"
#
$resolve<<EOD >> resolve.log_$cycle
hklin omit.mtz
labin FP=FP PHIB=PHIM FOM=FOMM
$combine_resolve_pattern
$combine_resolve_fragments
ha_file NONE
EOD
#
else
 if ( $ip ) echo "No image information generated this cycle"
 if ( -f combine.map ) rm -f combine.map
endif
#
#
#  Decide if we should keep going or quit...
#   quit if 
#   r_best_overall < n_unacceptable_r and 
#   it has been at least past_cycle cycles since things got better
#
@ cutoff_cycle = $cycle - $past_cycle
if ( ( $build_score_best_overall_cycle <= $cutoff_cycle ) && ( $r_best_overall < $n_unacceptable_r ) ) then
 echo "Best model was built on cycle $build_score_best_overall_cycle which was"
 echo "   more than $past_cycle cycles ago..."
 @ cycle ++
 break 
endif
#
end_of_cycle:
#
#  clean up...
#
if ( $clean_up == YES  && -f resolve.log_$cycle ) then
 if ( -f resolve.log_$cycle ) rm -f resolve.log_$cycle
 if ( -f refmac.log$cycle ) rm -f refmac.log_$cycle
 if ( -f resolve.mtz_$cycle ) rm -f resolve.mtz_$cycle
endif
#===========================================================================
#
if ( $ip ) echo "END OF CYCLE $cycle"
@ cycle ++
end
#
if ( $pdb_in != NONE ) then
 echo "Copying $pdb_in to resolve_group.pdb and resolve_best.pdb for rebuilding/evaluation..."
 cp $pdb_in resolve_best.pdb
 cp $pdb_in resolve_group.pdb
 goto start_rebuild
endif
#
echo "" > resolve_group.pdb
@ cyc = $cycle - 5
if ( $cyc < 1 ) set cyc = 1
while ( $cyc < $cycle )
 echo "Adding refine.pdb_$cyc to resolve_group.pdb"
 cat refine.pdb_$cyc >>resolve_group.pdb
 @ cyc ++
end
echo "Starting PDB group model for rebuild is in resolve_group.pdb"
#
if ( -f resolve_best.pdb ) then
 echo "Best model was built on cycle $build_score_best_overall_cycle"
 echo "This model is now in resolve_best.pdb"
else
 echo "No model was built...quitting"
 exit
endif
#
#===================================================================
if (  $quick_as_possible == YES ) then 
 goto wrap_up
endif
if (  $r_best_overall < $n_unacceptable_r ) then
 goto wrap_up
endif
if ( $n_cycle_rebuild_max == 0 ) then
 goto wrap_up
endif
 echo " "
 echo "Our best R-factor (0.$r_best_overall) is > cutoff of $unacceptable_r"
 echo "so we are going to try and rebuild the model..."
 echo " "
#===================================================================
#
start_rebuild:
#
set past_cycle = $n_cycle_rebuild_min
@ cycle_start = $cycle + 1
@ n_cycle = $cycle_start + $n_cycle_rebuild - 1
#
#
#
set pdb_start = resolve_group.pdb           #   starting PDB file
#
if ( $pdb_in == NONE ) then
 set scale = 1.0
 set extend = NO
 set reuse_chain = NO
else if ( $grep_phib != 0 ) then
#   we have phases, so use scale = 1.0...but we are rebuilding model, so
#    use extend and reuse_chain to get maximum building
 set scale = 1.0
 set extend = YES
 set reuse_chain = YES
else
#  slightly different rules for rebuilding from scratch from a model
#   without any phase information:
 set scale = 0.2
 set extend = YES
 set reuse_chain = YES
endif
#
set r_max = 0.$r_best_overall # maximum free R to write anything to resolve_best.pdb
#
if ( $ip ) echo "Output density-modified phases will be in image_only_dm_image.mtz_[cycle]"
#
#
if ( $evaluate_only != YES && $extend == NO )then
 echo "Model will not be extended each cycle"
 set extend_flag = ' 0 ' 
else if ( $evaluate_only != YES ) then
 echo "Model will be extended each cycle"
 set extend_flag = ' 0 1 2 ' 
endif
#
if ( $evaluate_only != YES && $omit_on_rebuild == YES )then
 echo "Omit prime-and-switch phasing will be used every  $n_cycle_rebuild_omit cycle"
 set omit = 'omit' 
else 
 echo "Standard image-based phasing will be used"
 set omit = ' ' 
endif
if ( $ip ) echo "Scale on map weighting is $scale"
set n_r_max = $r_max:e
if ( $cycle_start > 2 ) then
 echo "Maximum R to save model: 0.$n_r_max "
endif
#
echo " "
#
@ cycle = $cycle_start - 1
@ cycle_pdb_start = $cycle
#
echo "Copying ${pdb_start} to refine.pdb_$cycle"
cp ${pdb_start} refine.pdb_$cycle
#
#======================================================================
@ cycle ++
if ( $evaluate_only != YES ) then
echo "Starting rebuilding at cycle $cycle of total of $n_cycle cycles"
#=============================================================
#=============================================================
#
echo " "
echo "-------------------------------------------------------------"
echo "         RESIDUES  SIDE-CHAINS"
echo " CYCLE     BUILT      PLACED        SCORE      R       Free R   BEST"
endif
#
set n_omit_cycle = 0
while ( $cycle <= $n_cycle )
@ n_omit_cycle ++
if ( $n_cycle_rebuild_omit == 0 ) then
#  is_omit_cycle is NO for this one; defaults for omit_use and pdb_in
#
 set is_omit_cycle = NO
 set omit_use = $omit
 set pdb_in_current_model = 'pdb_in current_model.pdb'
else if ( $n_omit_cycle == 1 || $n_omit_cycle == $n_cycle_rebuild_omit ) then
 set n_omit_cycle = 1
 set is_omit_cycle = YES
 if ( $ip ) echo "Using $omit model phasing on this cycle and building from scratch"

 set omit_use = $omit
 set pdb_in_current_model = ' '
else
 set omit_use = ' '
 set pdb_in_current_model = 'pdb_in current_model.pdb'
 set is_omit_cycle = NO
endif
#
set cycle_prev = $cycle
@ cycle_prev --
#
if ( -f resolve.log_$cycle ) then
 if ( $ip ) echo "removing previously-existing log file resolve.log_$cycle"
 rm -f resolve.log_$cycle
endif
set build_score_best = -1000
set built_best = 0
set placed_best = 0
#
if ( -f work_model.pdb ) then
  if ( $ip ) echo "removing work_model.pdb"
  rm -f work_model.pdb
endif
#
if ( $build_outside_model == YES ) then
# put all current models in all.pdb
if ( -f all.pdb ) rm -f all.pdb
set c_cycle = $cycle_pdb_start
set c_end = $cycle_prev
while ($c_cycle <= $c_end)
if ( -f refine.pdb_${c_cycle} ) cat refine.pdb_${c_cycle} >> all.pdb
@ c_cycle ++
end
set add_mask = 'add_mask'
set model_all = 'model all.pdb'
set build_outside = 'build_outside_model'
else
set add_mask = ''
set model_all = ''
set build_outside = ''
endif
# 
#  get image of pdb_start and all models up to now into image.mtz. 
#  Reflns from hklin_exp
#
if ( $ip ) echo "Getting map of ${pdb_start} and all models built into binary file dump.map"
#
${resolve}<<EOD >> resolve.log_$cycle
hklin $hklin_exp
labin $labin_exp
labin $labin_exp_cont
resolution ${dmin} ${dmax}
composite_pdb refine.pdb_
composite_pdb_start $cycle_pdb_start
composite_pdb_end $cycle_prev
composite_all
dump_composite
EOD
#
if ( $ip ) echo "Getting image of ${pdb_start} and built models using dump.map into image_only.mtz"
#
${resolve}<<EOD >> resolve.log_$cycle
hklin $hklin_exp
labin $labin_exp
labin $labin_exp_cont
resolution ${dmin} ${dmax}
no_build
pattern_phase
prior_weight 1.0
hklout image_only.mtz
read_cc_map
cc_map_file dump.map
nohl
EOD
#  make sure this sort of worked...
if ( `cat resolve.ok` != 'OK' ) then
 echo "Sorry, resolve was not able to finish on the first cycle"
 echo "The end of the resolve log file says..."
 echo "---------------------------------------------------------"
 tail -10 resolve.log_$cycle
 echo "---------------------------------------------------------"
 exit
endif
#
if ( $ip ) echo "Density-modifying image of $pdb_start and prev models to get image_only_dm.mtz"
#
${resolve}<<EOD >> resolve.log_$cycle
hklin image_only.mtz
labin FP=FP SIGFP=SIGFP PHIB=PHIM FOM=FOMM HLA=HLAM HLB=HLBM HLC=HLCM HLD=HLDM
hklout image_only_dm.mtz
solvent_content $solvent_content  
prior_weight 1.0
database 5
resolution ${dmin} ${dmax}
$add_mask      ! optionally add mask based on all.pdb (model_all)
$model_all
mask_cycles 1
no_build
nohl
EOD
#
if ( $ip ) echo "Carrying out $omit_use phasing using image_only_dm.mtz "
if ( $ip ) echo " and ${pdb_start} and building model"
#
#
#   either build using this map, or build around starting model only!
#
if ( -f resolve.built ) rm -f resolve.built
if ( -f built_model.pdb ) rm -f built_model.pdb
if ( -f resolve.pdb ) rm -f resolve.pdb
#
if ( $build_outside_model == NO || $pdb_in == NONE )then
 if ($ip) echo "Standard build"
#
${resolve}<<EOD >> resolve.log_${cycle}
!standard build
hklin $hklin_exp
labin $labin_exp
labin $labin_exp_cont
hklstart image_only_dm.mtz
labstart FP=FP PHIB=PHIM FOM=FOMM 
hklout image_only_dm_image.mtz_${cycle}
solvent_content $solvent_content  
database 5
resolution ${dmin} ${dmax}
composite_pdb refine.pdb_
composite_pdb_start $cycle_pdb_start
composite_pdb_end $cycle_prev
composite_all
prior_weight 1.0
scale_refl $scale
image
$omit_use
mask_cycles 0
$superquick_build
ha_file $ha_file
seq_file $seq_file
compare_file $compare_file
$pdb_in_current_model
nohl
$mask_line
$use_mask_line
EOD
if ( -f resolve.pdb ) cp resolve.pdb work_model.pdb
if ( -f resolve.pdb ) cp resolve.pdb built_model.pdb
cp HEADER.pdb work_model_start.pdb
#
else
#
 if ($ip) echo "Building around model $pdb_in in 2 steps"
cp HEADER.pdb work_model.pdb
cp HEADER.pdb built_model.pdb
cat $pdb_in |$grep_type ATOM >>work_model.pdb
cat $pdb_in |$grep_type HETATM >>work_model.pdb
cp work_model.pdb work_model_start.pdb
#
${resolve}<<EOD >> resolve.log_${cycle}
! building around model $pdb_in in 2 steps
hklin $hklin_exp
labin $labin_exp
labin $labin_exp_cont
hklstart image_only_dm.mtz
labstart FP=FP PHIB=PHIM FOM=FOMM 
hklout image_only_dm_image.mtz_${cycle}
solvent_content $solvent_content  
database 5
resolution ${dmin} ${dmax}
composite_pdb refine.pdb_
composite_pdb_start $cycle_pdb_start
composite_pdb_end $cycle_prev
composite_all
prior_weight 1.0
scale_refl $scale
image
$omit_use
mask_cycles 0
no_build
ha_file $ha_file
seq_file $seq_file
compare_file $compare_file
$pdb_in_current_model
$add_mask   ! add mask from model to rmsd mask
$model_all
$mask_line
$use_mask_line
nohl
EOD
#
${resolve}<<EOD >> resolve.log_${cycle}
! Building outside input model $pdb_in only
hklin  image_only_dm_image.mtz_${cycle}
labin FP=FP PHIB=PHIM FOM=FOMM
ha_file NONE
add_mask
model work_model_start.pdb
database 5
solvent_content $solvent_content
resolution ${dmin} ${dmax}
nohl
build_only
build_outside_model
ha_file $ha_file
seq_file $seq_file
compare_file $compare_file
$superquick_build
EOD
 if ( -f resolve.built)    $grep_type ATOM resolve.pdb  >>work_model.pdb 
 if ( -f resolve.built)    $grep_type ATOM resolve.pdb  >>built_model.pdb
endif
#
if ( $omit_use == omit )then
 cp resolve_composite_map.mtz  build.mtz
 if ( $ip ) echo "Copied resolve_composite_map.mtz to build.mtz"
else
 cp image_only_dm_image.mtz_${cycle} build.mtz
 if ( $ip ) echo "Copied image_only_dm_image.mtz to build.mtz"
endif
 cp image_only_dm_image.mtz_${cycle} refine.mtz
 if ( $ip ) echo "Copied image_only_dm_image.mtz to refine.mtz"
#
#
if (  -f resolve.built ) then
    set built = `cat resolve.built|head -1` 
    set placed = `cat resolve.built | head -2|tail -1` 
    set build_score = `cat resolve.built | head -4|tail -1` 
    cp resolve.pdb built_model.pdb
else 
 set build_score = 0
 set built = 0
 set placed = 0
 echo "REMARK No model built" > work_model.pdb
 if ( -f built_model.pdb ) rm -f built_model.pdb
 echo "REMARK No model built" 
endif
set build_score_best = $build_score
set built_best = $built
set placed_best = $placed
 if ( $ip ) echo "Score : $build_score   built: $built  placed: $placed"
#
#==================================================================
#
#  Start of refine-extend cycles
#
if  ( $build_outside_model == YES && $pdb_in != NONE )then
set add_mask_use = 'add_mask'
set model_in = 'model work_model_start.pdb'
set build_outside_model_use = 'build_outside_model'
else 
set add_mask_use = ''
set model_in = ''
set build_outside_model_use = ''
endif
foreach extend_cycle ( $extend_flag )  
#
if ( $ip ) echo "Extend cycle $extend_cycle"
if( -f built_model.pdb && $extend_cycle > 0 ) then
 if ( $ip ) echo "Extending model"
#
if ( -f resolve.built ) rm -f resolve.built
${resolve}<<EOD >> resolve.log_$cycle
hklin build.mtz
labin FP=FP PHIB=PHIM FOM=FOMM 
hklout dummy.mtz
solvent_content $solvent_content
resolution ${dmin} ${dmax}
extend_only
$add_mask_use
$model_in     ! building outside work_model_start only
$build_outside_model_use    ! and map to the region of this model too
ha_file $ha_file
seq_file $seq_file
compare_file $compare_file
pdb_in built_model.pdb    ! extend what we just built (only)
no_trim
nohl
$mask_line
$use_mask_line
EOD
endif
if ( ( $extend_cycle > 0 ) && ( -f resolve.built ) ) then
    set built = `cat resolve.built|head -1` 
    set placed = `cat resolve.built | head -2|tail -1` 
    set build_score = `cat resolve.built | head -4|tail -1` 
 if ( $ip ) echo "Score for extension: $build_score   built: $built  placed: $placed"
else if ( $extend_cycle > 0 ) then
 set build_score = 0
 set built = 0
 set placed = 0
endif
if ( ( $extend_cycle > 0 ) && ( $build_score > $build_score_best ) ) then
 set build_score_best = $build_score
 set built_best = $built
 set placed_best = $placed
 cp resolve.pdb built_model.pdb
 cp work_model_start.pdb work_model.pdb
 $grep_type ATOM built_model.pdb >>work_model.pdb
 if ( $ip ) echo "Adding extended model to starting work model"
else if ( $extend_cycle > 0 ) then
 if ( $ip ) echo "Extension did not improve score...going on"
#  break to end of the extend cycles...
 break
endif
#
#
#  decide if we refine model this cycle..
#
if ( -f resolve.coverage ) then
 set coverage1 = `cat resolve.coverage|head -3|tail -1`
 set coverage = 0
 if ( $#coverage1> 0 ) set coverage= $coverage1[$#coverage1-$#coverage1]:e
else
 set coverage = 0
endif
if ( $refine_model == NO ) then
 set refine_it = NO
else if ( $build_score_best < 2 ) then
 set refine_it = NO
else if ( ( $coverage < $minimum_coverage ) && ( $build_outside_model == NO) )then
 set refine_it = NO
else if ( ! -f work_model.pdb ) then
 echo "REMARK No model built" > work_model.pdb
 echo "REMARK No model built extend cycle $extend_cycle"
 set refine_it = NO
else
 set refine_it = YES
endif
if ( $refine_it == YES ) then
 set ncyc = 20
else
 set ncyc = 0
endif
if ( ( $refmac5 != NONE ) && ( $build_score_best > 1 ) ) then
if ( $ip ) echo "Refining model $ncyc cycles and placing it in refine.pdb_$cycle"
#
cat work_model.pdb | $grep_type -v HETATM > resolve.pdb
#
# --------refine with refine.--------------
#
#  NOTE: refine against "best" map = refine.mtz
#
cp resolve.pdb refine.pdb   # in case it doesn't work...
#
$refmac5 xyzin resolve.pdb hklin refine.mtz xyzout refine.pdb  \
PROTOUT $CCP4_SCR/protout.dat \
PROTCOUNTS $CCP4_SCR/counts.dat \
PROTSCR $CCP4_SCR/counts.scr \
<<EOD > refmac.log_${cycle}
MAKE HYDRogens No
MAKE CHECK 0
LABI FP=FP SIGFP=SIGFP HLA=HLAM HLB=HLBM HLC=HLCM HLD=HLDM FREE=FreeR_flag
REFI TYPE RESTrained RESOlution ${dmax} ${dmin}
#   Maximum likelihood refinement
REFI RESI MLKF
#   Overall B value refinement
REFI BREF OVER      ! Refine overall B-values
#   Tight restraints  Lowered further 051403
WEIG MATR 0.05
#   Scale down shifts at every cycle by factor 0.5
DAMP 0.5 0.5
#   No Blur (scale down reliability of) phases
REFI PHASed  BBLUrring 00.0  SCBLurring 1.
#  Babinet's bulk solvent scale parameters. But fix B value
$simple_scaling_line
SCALe LSSCale FIXBulk BBULk 200.0
#   Anisotropic scaling. It is default
SCALe LSSCale ANISotropic
#   Fix Babinet's bulkd solvent for sigmaA
SCAL MLSC FIXBulk BBULk 200.0 SCBUlk -0.05
NCYC $ncyc
MONI MEDI
BINS 10
EOD
#----------------------------------------------------------
#
cp refine.pdb refine.pdb_$cycle
else if (  $build_score_best > 1  ) then
#   no refmac5...just copy over
 cat work_model.pdb | $grep_type -v HETATM > resolve.pdb
 cp resolve.pdb refine.pdb
 cp refine.pdb refine.pdb_$cycle
else
 if ( $ip ) echo "No model built this try"
 echo "REMARK  No model built cycle $cycle " > refine.pdb_$cycle
 echo "REMARK  No model built cycle $cycle "
endif
#
end_of_extend:
end
#
#  Get R and freeR value on refine.pdb_$cycle:
set freer_line = `$grep_type 'FREE R VALUE                     :' refine.pdb_$cycle|tail -1`
set r_line = `$grep_type 'R VALUE            (WORKING SET) :' refine.pdb_$cycle|tail -1`
#
if ( $#freer_line > 0 ) then
 set freer = $freer_line[$#freer_line-$#freer_line]:e
else
 set freer = 99999
endif
if ( $#r_line > 0 ) then
 set r = $r_line[$#r_line-$#r_line]:e
else
 set r = 99999
endif
if ( $ncyc == 0 )then
 set freer = $r
endif
#
if ( $ip ) echo " R/freeR for this model: 0.$r/0.$freer"
#
if ( $cycle < 10 ) then
 set b1 = " "
else
 set b1 = ""
endif
if ( $built_best < 100 ) then
 set b2 = " "
else
 set b2 = ""
endif
if ( $placed_best >= 100 ) then
 set b3 = ""
else if ( $placed_best >= 10 )then
 set b3 = " "
else
 set b3 = "  "
endif
if ( $build_score_best < 100 )then
 set b4 = " "
else
 set b4 = ""
endif
# Decide how to score. If R < unacceptable_r or r_best_overall < unacceptable_r
#   then use R. Otherwise use score
#
if ( ( $r < $n_unacceptable_r ) && ( $r < $r_best_overall ) ) then
 set this_is_best = YES
else if ( ( $r_best_overall >= $n_unacceptable_r ) && \
   ( $r <= $n_r_max ) && \
   ( $build_score_best > $build_score_best_overall ) ) then
 set this_is_best = YES
else
 set this_is_best = NO
endif
if ( $this_is_best == YES ) then
 if ( $ip ) echo "This is the best model overall so far. Previous best had score of $build_score_best_overall"
 if ( $ip ) echo "Copied refine.pdb_$cycle to resolve_best.pdb"
 if ( $ip ) echo "Copied resolve.log_$cycle resolve_best.log"
 if ( $ip ) echo "Copied refine.mtz to resolve_best.mtz"
 @ build_score_best_overall = $build_score_best
 @ build_score_best_overall_cycle = $cycle
 @ r_best_overall = $r
 cp refine.pdb_$cycle resolve_best.pdb
 cp resolve.log_$cycle resolve_best.log
 cp refine.mtz resolve_best.mtz
 if( -f refmac.log_$cycle ) cp refmac.log_$cycle refmac_best.log
 echo "  ${b1}$cycle       ${b2}$built_best         ${b3}$placed_best          ${b4}$build_score_best     0.$r    0.$freer    *"
else
 echo "  ${b1}$cycle       ${b2}$built_best         ${b3}$placed_best          ${b4}$build_score_best     0.$r    0.$freer    "
endif
#
#==================================================================
#   
if ( $reuse_chain == YES )then
 cat  refine.pdb > current_model.pdb
 if ( $ip ) echo " Using refine.pdb fragments as part of next build cycle"
else if ( -f current_model.pdb ) then
  rm -f current_model.pdb
endif
#
#
if ( $ip ) echo " "
if ( $ip ) echo "End of cycle $cycle with score of $build_score_best"
if ( $ip ) echo " "
#
#  Decide if we should keep going or quit..
@ cutoff_cycle = $cycle - $past_cycle
if ( ( $build_score_best_overall_cycle <= $cutoff_cycle ) && ( $r_best_overall < $n_unacceptable_r ) ) then
 echo "Our best model was built on cycle $build_score_best_overall_cycle which was"
 echo "   more than $past_cycle cycles ago...calling it quits with this model,"
 echo "   now in resolve_best.pdb"
 break 
endif
#
#  clean up...
#
if ( $clean_up == YES  ) then
 if ( -f resolve.log_$cycle ) rm -f resolve.log_$cycle
 if ( -f refmac.log_$cycle ) rm -f refmac.log_$cycle
 if ( -f resolve.mtz_$cycle ) rm -f resolve.mtz_$cycle
 if ( -f refine.mtz ) rm -f refine.mtz
 if ( -f image_only_dm_image.mtz_$cycle ) rm -f image_only_dm_image.mtz_$cycle
endif
#
if ( is_omit_cycle == YES ) then
#  set this cycle as first one to consider when creating composite from PDB files
 set cycle_pdb_start = $cycle
endif
#
@ cycle ++
#
#   end of build cycle
#
end
#
wrap_up:
if ($evaluate_only == NEVER)then
  exit
endif
#
#==========================================================================
#==========================================================================
echo " "
echo "Analyzing model based on prime-and-switch composite omit map"
#
echo "First rebuilding current model to remove refinement bias..."
if ( ! -f resolve_best.pdb ) then
 echo "Sorry cannot find resolve_best.pdb?"
 exit
else
 cp resolve_best.pdb test.pdb_0
endif
if ( $skip == YES ) then
 goto evaluate_it
endif
#
${resolve}<<EOD > resolve_eval.log_1
hklin $hklin_exp
labin $labin_exp
labin $labin_exp_cont
resolution ${dmin} ${dmax}
composite_pdb test.pdb_
composite_pdb_start 0
composite_pdb_end 0
composite_all
dump_composite
EOD
#
${resolve}<<EOD >> resolve_eval.log_1
hklin $hklin_exp
labin $labin_exp
labin $labin_exp_cont
resolution ${dmin} ${dmax}
no_build
pattern_phase
prior_weight 1.0
hklout image_only.mtz
read_cc_map
cc_map_file dump.map
nohl
EOD
#
if ( $ip ) echo "Density-modifying image of resolve_best.pdb to get image_only_dm.mtz"
#
${resolve}<<EOD >> resolve_eval.log_1
hklin image_only.mtz
labin FP=FP SIGFP=SIGFP PHIB=PHIM FOM=FOMM HLA=HLAM HLB=HLBM HLC=HLCM HLD=HLDM
hklout image_only_dm.mtz
solvent_content $solvent_content  
prior_weight 1.0
database 5
resolution ${dmin} ${dmax}
no_build
nohl
EOD
#
if ( $ip ) echo "Carrying out phasing using image_only_dm.mtz "
if ( $ip ) echo " and resolve_best.pdb and rebuilding model"
#
if ( -f resolve.built ) rm -f resolve.built
${resolve}<<EOD >> resolve_eval.log_1
hklin $hklin_exp
labin $labin_exp
labin $labin_exp_cont
hklstart image_only_dm.mtz
labstart FP=FP PHIB=PHIM FOM=FOMM 
hklout image_only_dm_image.mtz
solvent_content $solvent_content  
database 5
resolution ${dmin} ${dmax}
composite_pdb test.pdb_
composite_pdb_start 0
composite_pdb_end 0
composite_all
prior_weight 1.0
scale_refl $scale
image
mask_cycles 0
$superquick_build
seq_file $seq_file
compare_file $compare_file
ha_file $ha_file
hklout test.mtz_1
nohl
$mask_line
$use_mask_line
EOD
if ( -f resolve.pdb ) then
 cp resolve.pdb test.pdb_1
else
 echo "Sorry, no model obtained on rebuilding...quitting"
 exit
endif
if (  -f resolve.built  ) then
    set built = `cat resolve.built|head -1` 
    set placed = `cat resolve.built | head -2|tail -1` 
    set build_score = `cat resolve.built | head -4|tail -1` 
else 
 set build_score = 0
 set built = 0
 set placed = 0
 echo "Sorry, no model obtained on rebuilding...quitting"
 exit
endif
#
#  report what happened in our rebuilding...
echo "Rebuilt model...$built residues built, $placed with side-chains in test.pdb_1"
#
echo "Generating composite prime-and-switch composite omit map based on rebuilt model"
#
${resolve}<<EOD >> resolve_eval.log_1
hklin $hklin_exp
labin $labin_exp
labin $labin_exp_cont
resolution ${dmin} ${dmax}
composite_pdb test.pdb_
composite_pdb_start 1
composite_pdb_end 1
composite_all
dump_composite
EOD
#
${resolve}<<EOD >> resolve_eval.log_1
hklin $hklin_exp
labin $labin_exp
labin $labin_exp_cont
resolution ${dmin} ${dmax}
no_build
pattern_phase
prior_weight 1.0
hklout image_only.mtz
read_cc_map
cc_map_file dump.map
nohl
EOD
#
if ( $ip ) echo "Density-modifying image of test.pdb_1 to get image_only_dm.mtz"
#
${resolve}<<EOD >> resolve_eval.log_1
hklin image_only.mtz
labin FP=FP SIGFP=SIGFP PHIB=PHIM FOM=FOMM HLA=HLAM HLB=HLBM HLC=HLCM HLD=HLDM
hklout image_only_dm.mtz
solvent_content $solvent_content  
prior_weight 1.0
database 5
resolution ${dmin} ${dmax}
no_build
nohl
$mask_line
$use_mask_line
EOD
#
if ( $ip ) echo "Getting composite ps omit map using image_only_dm.mtz "
if ( $ip ) echo " and test.pdb_1"
#
${resolve}<<EOD >> resolve_eval.log_1
hklin $hklin_exp
labin $labin_exp
labin $labin_exp_cont
hklstart image_only_dm.mtz
labstart FP=FP PHIB=PHIM FOM=FOMM  
hklout image_only_dm_image.mtz
solvent_content $solvent_content  
database 5
resolution ${dmin} ${dmax}
composite_pdb test.pdb_
composite_pdb_start 1
composite_pdb_end 1
composite_all
prior_weight 1.0
scale_refl $scale
image
omit
mask_cycles 0
no_build
nohl
$mask_line
$use_mask_line
EOD
cp resolve_composite_map.mtz build.mtz
if ( $ip ) echo "Composite omit is in build.mtz"
#
evaluate_it:
#
#  Get R-factor for the model resolve_best.pdb...relative to hklin_exp
#
#
$refmac5 xyzin resolve_best.pdb hklin $hklin_exp xyzout refine.pdb  \
PROTOUT $CCP4_SCR/protout.dat \
PROTCOUNTS $CCP4_SCR/counts.dat \
PROTSCR $CCP4_SCR/counts.scr \
<<EOD > refmac_eval.log_${cycle}
MAKE HYDRogens No
MAKE CHECK 0
LABI FP=FP SIGFP=SIGFP FREE=FreeR_flag
REFI TYPE RESTrained RESOlution ${dmax} ${dmin}
#   Maximum likelihood refinement
REFI RESI MLKF
#   Overall B value refinement
REFI BREF OVER      ! Refine overall B-values
#   Tight restraints  Lowered further 051403
WEIG MATR 0.05
#   Scale down shifts at every cycle by factor 0.5
DAMP 0.5 0.5
#   No Blur (scale down reliability of) phases
REFI PHASed  BBLUrring 00.0  SCBLurring 1.
#  Babinet's bulk solvent scale parameters. But fix B value
$simple_scaling_line
SCALe LSSCale FIXBulk BBULk 200.0
#   Anisotropic scaling. It is default
SCALe LSSCale ANISotropic
#   Fix Babinet's bulkd solvent for sigmaA
SCAL MLSC FIXBulk BBULk 200.0 SCBUlk -0.05
NCYC 0
MONI MEDI
BINS 10
EOD
#----------------------------------------------------------
#
#  Get R and freeR value on refine.pdb:
set freer_line = `$grep_type 'FREE R VALUE                     :' refine.pdb|tail -1`
set r_line = `$grep_type 'R VALUE            (WORKING SET) :' refine.pdb|tail -1`
#
if ( $#freer_line > 0 ) then
 set freer = $freer_line[$#freer_line-$#freer_line]:e
else
 set freer = 99999
endif
if ( $#r_line > 0 ) then
 set r = $r_line[$#r_line-$#r_line]:e
else
 set r = 99999
endif
#
echo " "
echo "======================================================================"
if ( $evaluate_only == YES ) then
 echo "Work R-factor for input model $pdb_in : 0.$r   free R: 0.$freer"
 echo "Evaluation of $pdb_in using prime-and-switch composite omit map:"
else
 echo "Work R-factor for resolve_best.pdb: 0.$r   free R: 0.$freer"
 echo "Evaluation using prime-and-switch composite omit map:"
endif
#
#
#   Evaluate the model relative to the composite omit ps map:
#
${resolve} <<EOD > evaluate.log
hklin build.mtz
labin FP=FP PHIB=PHIM FOM=FOMM
model resolve_best.pdb
evaluate_model
EOD
#
#==========================================================================
echo "See evaluate.log for more details on the analysis of this model"
echo " "
echo "All done...please carefully evaluate your model."
echo "Your best indicators are the free R factor and the match"
echo "to the prime-and-switch composite omit map in evaluate.log."
echo "Note that no waters are built in this model."
echo " "
echo "Output files are:"
if ( -f resolve_best.pdb ) then
 echo "resolve_best.pdb  ... partially-refined model"
endif
if ( -f resolve_best.log ) then
 echo "resolve_best.log  ... log file in which this model was built"
endif
if ( -f resolve_best.mtz ) then
$resolve<<EOD>final.log
hklin resolve_best.mtz
labin FP=FP PHIB=PHIM FOM=FOMM
mask_cycles 1
minor_cycles 0
no_build
ccp4_map_file resolve_best.map
EOD
 echo "resolve_best.mtz  ... FP, SIGFP, best overall phases, freeR_flag...."
 echo "resolve_best.map  ... ccp4-style map from resolve_best.mtz"
endif
if ( -f refmac_best.log ) then
 echo "refmac_best.log   ... log file for refinement of resolve_best.pdb"
endif
if ( -f build.mtz ) then
$resolve <<EOD>> final.log
hklin build.mtz
labin FP=FP PHIB=PHIM FOM=FOMM
mask_cycles 1
minor_cycles 0
no_build
ccp4_map_file resolve_composite_omit.map
EOD
#
 echo "resolve_composite_omit.map ... composite omit map "
endif
#==================================================================