#! /bin/tcsh -f
#
#	estimate the F000 level of the 2fofc map from 
#	a pdb file and mtz file
#
#

source /local/programs/phenix-1.6.1-357/phenix_env
setenv PHENIX_OVERWRITE_ALL true
set tempfile = ${CCP4_SCR}/F000_temp$$

set pdbfile = ""
set mtzfile = ""
set reso = "2"
set ksol = ""
set Bsol = ""

# a few command-line options
foreach arg ( $* )
    if("$arg" =~ *.pdb) set pdbfile = "$arg"
    if("$arg" =~ *.mtz) set mtzfile = "$arg"
    if("$arg" =~ *A) set reso = `echo $arg | awk '{print $1+0}'`
    if("$arg" =~ ksol=*) set ksol = `echo $arg | awk -F "=" '{print $2+0}'`
    if("$arg" =~ Bsol=*) set Bsol = `echo $arg | awk -F "=" '{print $2+0}'`
end


# figure out space group from PDB file
set CELL = `awk '/^CRYST1/{print $2,$3,$4,$5,$6,$7;exit}' $pdbfile`
if( $#CELL != 6) then
    set BAD = "no cell in PDB"
    goto exit
endif
set pdbSG = `awk '/^CRYST1/{SG=substr($0,56,14);if(length(SG)==14)while(gsub(/[^ ]$/,"",SG));print SG;exit}' $pdbfile | head -1`
if("$pdbSG" == "R 32") set pdbSG = "R 3 2"
#if("$pdbSG" == "P 21") set pdbSG = "P 1 21 1"
#if("$pdbSG" == "A 2") set pdbSG = "A 1 2 1"
if("$pdbSG" == "I 21") set pdbSG = "I 1 21 1"
if("$pdbSG" == "A 1") set pdbSG = "P 1"
if("$pdbSG" == "P 21 21 2 A") set pdbSG = "P 21 21 2"
if("$pdbSG" == "P 1-") set pdbSG = "P -1"
if("$pdbSG" == "R 3 2" && $CELL[6] == 120.00) set pdbSG = "H 3 2"
if("$pdbSG" == "R 3" && $CELL[6] == 120.00) set pdbSG = "H 3"
set SG = `awk -v pdbSG="$pdbSG" -F "[\047]" 'pdbSG==$2 || pdbSG==$4{print;exit}' ${CLIBD}/symop.lib | awk '{print $4}'`
if("$SG" == "") set SG = "$pdbSG"


# need cell volume to compute F000 from map offsets
echo $CELL |\
awk 'NF==6{DTR=atan2(1,1)/45; A=cos(DTR*$4); B=cos(DTR*$5); G=cos(DTR*$6); \
 skew = 1 + 2*A*B*G - A*A - B*B - G*G ; if(skew < 0) skew = -skew;\
 printf "%.3f\n", $1*$2*$3*sqrt(skew)}' |\
cat >! ${tempfile}volume
set CELLvolume = `cat ${tempfile}volume`
rm -f ${tempfile}volume
rm -f ${tempfile}maphead.txt




ls -1 test* | awk '/edits$|cif$/{print "rm -f",$1}' | tcsh
cp $pdbfile test.pdb
echo "ready, set ..."
phenix.ready_set test.pdb >! ready_set.log
if($status) then
    set BAD = "phenix.ready_set failed"
    goto exit
endif

if("$ksol" == "" || "$Bsol" == "") then
cp test.updated.pdb phenix.pdb
set edits = `ls test* | awk '/edits$|cif$/'`
echo "go. "
rm -f phenix_refine_001.pdb phenix_refine_001_map_coeffs.mtz
rm -f refined.pdb map_coeffs.mtz
phenix.refine phenix.pdb $edits $mtzfile \
  apply_volume_scaling=True \
  strategy=rigid_body \
  xray_data.high_resolution=$reso \
  refinement.input.xray_data.r_free_flags.generate=True \
  main.number_of_macro_cycles=1 | tail
if($status || ! -e phenix_refine_001.pdb) then
    set BAD = "phenix.refine failed"
    goto exit
endif
cp phenix_refine_001.pdb refined.pdb
cp phenix_refine_001_map_coeffs.mtz map_coeffs.mtz

set ksol = `awk -F "ksol=" '/ksol=/{print $2+0}' phenix_refine_001.log | tail -1`
set Bsol = `awk -F "Bsol=" '/Bsol=/{print $2+0}' phenix_refine_001.log | tail -1`
#set phenix_reso = `awk -F "resolution:" '/resolution:/{print $2+0}' phenix_refine_001.log | tail -1`
#set finereso = `echo $phenix_reso | awk '{print $1/2}'`
#if("$reso" == "") set reso = "$finereso"

echo "refined bulk solvent paramters: $ksol $Bsol"

endif

# set all B factors to 80
cat test.updated.pdb |\
awk '/^ATOM/ || /^HETAT/{print substr($0,1,60)" 80.00"substr($0,67);next} {print}' |\
cat >! refined.pdb

# generate phenix-derived model Fs (but without F000)
echo "generating phenix models..."
echo "regular"
rm -f fmodel.mtz
phenix.fmodel refined.pdb k_sol=$ksol b_sol=$Bsol high_res=$reso \
  algorithm=direct file_name=fmodel.mtz > /dev/null

# now with bulk solvent turned off
echo "no solvent"
rm -f nobulk.mtz
phenix.fmodel refined.pdb k_sol=0 b_sol=$Bsol high_res=$reso \
  algorithm=direct file_name=nobulk.mtz > /dev/null

# sharp solvent mask (should have same Q)
echo "sharp solvent"
rm -f ss.mtz
phenix.fmodel refined.pdb k_sol=$ksol b_sol=0 high_res=$reso \
  algorithm=direct file_name=ss.mtz > /dev/null

# isotropic Bs (to compare with SFALL)
echo "isotropic Bs"
egrep -v "ANISO" refined.pdb >! iso.pdb
rm -f iso.mtz
phenix.fmodel iso.pdb k_sol=0 b_sol=0 high_res=$reso \
  algorithm=direct file_name=iso.mtz > /dev/null




# compare with sfall (no anisotropic Bs! )
echo "making sfall map for vacuum level reference"
sfall xyzin iso.pdb mapout ${tempfile}sfall.map << EOF > /dev/null
mode atmmap
SYMM $SG
resolution $reso 1000
EOF

if($?doublegrid) then
# examine this map to get the grid and axis convention for this SG
echo "" | mapdump mapin ${tempfile}sfall.map >! ${tempfile}maphead.txt
set xyzgrid = `awk '/Grid sampling on x, y, z ../{print $8,$9,$10}' ${tempfile}maphead.txt  | tail -1`
# 2x finer than default sampling...
set biggrid = `echo $xyzgrid 2 | awk '{print $1*$NF,$2*$NF,$3*$NF}'`
sfall xyzin iso.pdb mapout ${tempfile}sfall.map << EOF > /dev/null
mode atmmap
SYMM $SG
GRID $biggrid
EOF
endif
echo "xyzlim ASU" | mapmask mapin ${tempfile}sfall.map mapout sfall.map > /dev/null


# examine this map to get the grid and axis convention for this SG
echo "" | mapdump mapin sfall.map >! ${tempfile}maphead.txt
set xyzgrid = `awk '/Grid sampling on x, y, z ../{print $8,$9,$10}' ${tempfile}maphead.txt  | tail -1`
echo "using map grid: $xyzgrid"









# make the maps
echo "calculating maps of phenix models..."
fft hklin map_coeffs.mtz mapout ${tempfile}ffted.map << EOF > /dev/null
labin F1=2FOFCWT PHI=PH2FOFCWT
GRID $xyzgrid
EOF
echo "xyzlim asu" |\
 mapmask mapin ${tempfile}ffted.map mapout 2fofc.map > /dev/null
rm -f ${tempfile}ffted.map


fft hklin fmodel.mtz mapout ${tempfile}ffted.map << EOF > /dev/null
labin F1=FMODEL PHI=PHIFMODEL
GRID $xyzgrid
EOF
echo "xyzlim asu" |\
 mapmask mapin ${tempfile}ffted.map mapout fmodel.map > /dev/null
rm -f ${tempfile}ffted.map

fft hklin ss.mtz mapout ${tempfile}ffted.map << EOF > /dev/null
labin F1=FMODEL PHI=PHIFMODEL
GRID $xyzgrid
EOF
echo "xyzlim asu" |\
 mapmask mapin ${tempfile}ffted.map mapout fmodel_ss.map > /dev/null
rm -f ${tempfile}ffted.map


fft hklin nobulk.mtz mapout ${tempfile}ffted.map << EOF > /dev/null
labin F1=FMODEL PHI=PHIFMODEL
GRID $xyzgrid
EOF
echo "xyzlim asu" |\
 mapmask mapin ${tempfile}ffted.map mapout nobulk.map > /dev/null
rm -f ${tempfile}ffted.map


fft hklin iso.mtz mapout ${tempfile}ffted.map << EOF > /dev/null
labin F1=FMODEL PHI=PHIFMODEL
GRID $xyzgrid
EOF
echo "xyzlim asu" |\
 mapmask mapin ${tempfile}ffted.map mapout iso.map > /dev/null
rm -f ${tempfile}ffted.map





# subtract to isolate solvent region
echo "subtracting nobulk map from sharp-solvent map"
echo scale factor -1 | mapmask mapin nobulk.map mapout neg.map > /dev/null
echo "maps add" | mapmask mapin1 fmodel_ss.map mapin2 neg.map mapout sharpsolvent.map > /dev/null

# this is the median difference between solvent on and off (proportional to total charge)
echo "finding median with outlier rejection"
./map_vacuum_level.com sharpsolvent.map >! sharpsolvent_vac.log 
grep "vacuum level:" sharpsolvent_vac.log
set solvent_vac = `awk '/vacuum level:/{print $4}' sharpsolvent_vac.log`
set solvent_sig = `awk '/vacuum level:/{print $6}' sharpsolvent_vac.log`
cp hist0.txt sharpsolvent_hist0.txt


# subtract sfall map to measure offset
echo "subtracting sfall map from phenix map"
echo scale factor -1 | mapmask mapin sfall.map mapout neg.map > /dev/null
echo "maps add" | mapmask mapin1 iso.map mapin2 neg.map mapout baseline.map > /dev/null


# this is the median difference between phenix and sfall maps (sfall vacuum=0)
echo "finding median with outlier rejection"
./map_vacuum_level.com baseline.map >! baseline_vac.log
grep "vacuum level:" baseline_vac.log
set model_vac = `awk '/vacuum level:/{print $4}' baseline_vac.log`
set model_sig = `awk '/vacuum level:/{print $6}' baseline_vac.log`
cp hist0.txt baseline_hist0.txt


# this should be the vacuum level of the 2fofc (and fmodel) maps
set baseline = `echo $model_vac $solvent_vac | awk '{print -$1-$2}'`
set baseline_sig = `echo $model_sig  $solvent_sig | awk '{print sqrt($1*$1+$2*$2)}'`

echo scale factor 1 $baseline |\
 mapmask mapin fmodel.map mapout fmodel_vac0.map > /dev/null
echo "" | mapdump mapin fmodel_vac0.map | grep dens >! dens.txt
set mean = `awk '/Mean density/{print $NF}' dens.txt`
set F000 = `echo $CELLvolume $mean 0 | awk '{print $1*($2-$3)}'`
set SIGF000 = `echo $CELLvolume $baseline_sig | awk '{print $1*$2}'`
echo "fmodel_vac0.map should have a vacuum level of zero and F000 = $F000 +/- $SIGF000"

echo scale factor 1 $baseline |\
 mapmask mapin 2fofc.map mapout 2fofc_vac0.map > /dev/null
echo "" | mapdump mapin 2fofc_vac0.map | grep dens >! dens.txt
set mean = `awk '/Mean density/{print $NF}' dens.txt`
set F000 = `echo $CELLvolume $mean 0 | awk '{print $1*($2-$3)}'`
echo "2fofc_vac0.map should have a vacuum level of zero"


exit:
if($?BAD) then
    echo "ERROR: $BAD"
    exit 9
endif

exit

######################################################################
#
#	notes and tests
#
#
~/Develop/randompdb.com 10 10 10 90 90 90
awk '{gsub("OW1 WAT","C   IUM");print}' random.pdb > ! refme.pdb
egrep -v "ANISO|WAT|HOH|REM" decent.pdb >! refme.pdb

ano_sfall.com refme.pdb 1A

echo labin file 1 E1=FP E2=SIGFP | cad hklin1 ideal_ano.mtz hklout refme.mtz

./find_F000.com refme.pdb refme.mtz


foreach pdb ( 3lzt 4lzt 1a6g )

./find_F000.com ${pdb}.pdb ${pdb}.mtz 2A | tee ${pdb}_F000.log
cp fmodel_vac0.map ${pdb}_fmodel_vac0.map
cp 2fofc_vac0.map ${pdb}_2fofc_vac0.map

end