#!/bin/csh

#
# Use DYNAMO to create a protein PDB file with protons from crystal structure:
#
#  1. Extract coords for existing crystal structure.
#
#  2. Create DYNAMO GMC from existing PDB sequence, anneal an initial structure
#     to conform to the crystal coords.
#
#  3. Remove residues which were missing in the existing PDB sequence.

if ($#argv < 2) then
   echo "Create a PDB File with Protons from an Existing PDB File:"
   echo " addH.com inPDBName outPDBName  [addHOptions] -- [pdb2gmcArgs]"
   echo "addH Options:"
   echo " -gmc       gmcDir  [addH.gmc] GMC Directory Name for Computation."
   echo " -seg       segName [A]        Segment Name for Output PDB."
   echo " -fc        fc      [500.0]    Force Constant for AC Term." 
   echo " -dHi       dHi     [0.0]      Target AC Distance, Angstroms."
   echo " -coolSteps cStep   [16000]    Number of Cooling Steps for Annealing."
   echo " -print     pStep   [100]      Annealing Steps Between Reports."
   echo "Multi-Chain Mode:"
   echo " -auto                         Determine Mode Automatically (Default)."
   echo " -noauto                       No Auto Determination."
   echo " -multi                        Use Multi-Chain PDB Input Mode."
   echo " -nomulti                      Use Single-Chain PDB Input Mode."
   echo "Restoring Original Heavy Atom Coords After Annealing:"
   echo " -restore                      Restore Original Heavy Atom Coords (Default)."
   echo " -norestore                    Keep Simulated Annealing Coords."
   echo "Other Options:"
   echo " -center                       Center PDB Coordinates."
   echo " -nocenter                     Keep PDB Coordinates (Default)."
   echo " -remark                       Keep REMARK Lines, Etc."
   echo " -hetatm                       Keep HETATM Lines."
   echo " -del                          Delete Missing Residues (Default)."
   echo " -nodel                        Replace Missing Residues with GLY."
   echo " -eval                         Evaluate Results."
   echo " -tk                           Use TK Graphics."
   exit 0
endif

set inPDBName   = $argv[1]
set outPDBName  = $argv[2]
set inPDBTitle  = (`getArgD $argv -inTitle $inPDBName`)
set adjPDBName  = (`getArgD $argv -adj     addH_adj.pdb`)
set tmpPDBName  = (`getArgD $argv -tmp     addH_tmp.pdb`)
set saPDBName   = (`getArgD $argv -sa      addH_sa.pdb`)
set gmcName     = (`getArgD $argv -gmc addH.gmc`)
set segName     = (`getArgD $argv -seg A`)
set dHi         = (`getArgD $argv -dHi 0.0`)
set fc          = (`getArgD $argv -fc  500.0`)
set coolSteps   = (`getArgD $argv -coolSteps 16000`)
set printStep   = (`getArgD $argv -print     100`)
set dynArgs     = ""
set centerFlag  = nocenter
set tkFlag      = notk
set delFlag     = del
set restoreFlag = 1
set remFlag     = 0
set hetFlag     = 0
set multiFlag   = 0
set autoFlag    = 1
set eFlag       = 0
set chainList   = NONE 

set inPDBNameOrig = $inPDBName

#
# Hidden option for server use:

if (`flagLoc $argv -saAdvanced`) then
   set restoreFlag = 0
endif

if (`flagLoc $argv -tk`) then
   set tkFlag = tk
endif
 
if (`flagLoc $argv -graph`) then
   set tkFlag = tk
endif

if (`flagLoc $argv -del`) then
   set delFlag = del
endif
 
if (`flagLoc $argv -nodel`) then
   set delFlag = nodel
endif

if (`flagLoc $argv -restore`) then
   set restoreFlag = 1
endif

if (`flagLoc $argv -norestore`) then
   set restoreFlag = 0
endif

if (`flagLoc $argv -center`) then
   set centerFlag = center 
endif

if (`flagLoc $argv -nocenter`) then
   set centerFlag = nocenter
endif

if (`flagLoc $argv -auto`) then
   set autoFlag = 1
endif

if (`flagLoc $argv -noauto`) then
   set autoFlag = 0
endif

if (`flagLoc $argv -remark`) then
   set remFlag = 1
endif

if (`flagLoc $argv -noremark`) then
   set remFlag = 0
endif

if (`flagLoc $argv -hetatm`) then
   set hetFlag = 1
endif

if (`flagLoc $argv -nohetatm`) then
   set hetFlag = 0
endif

#
# Hidden Option for addHMulti use:

if (`flagLoc $argv -addHMulti`) then
   set remFlag = 0
endif

set loc = (`flagLoc $argv --`)

if ($loc != 0) then
   @ loc++

   if ($loc <= $#argv) then
      set dynArgs = ($argv[$loc-])
      echo $dynArgs
   endif
endif

if (!(-e $inPDBName)) then
    echo "AddH Error Finding Input PDB $inPDBName"
    exit 1
endif

if ((-d $gmcName)) then
    echo "AddH: Removing Existing GMC directory $gmcName"
   /bin/rm -rf $gmcName
endif

set chainList = (`getPDBInfo.tcl -in $inPDBName -noverb -info chainNames`)

if ($autoFlag != 0) then
   if ($#chainList == 1) then
      set multiFlag = 0
   else
      set multiFlag = 1
   endif
endif

if (`flagLoc $argv -multi`) then
   set multiFlag = 1
endif

if (`flagLoc $argv -nomulti`) then
   set multiFlag = 0
endif

if ($multiFlag == 0) then
   echo "AddH: Using Single-Chain Mode."
   set multiArg = "nomulti"
else
   echo "AddH: Using Multi-Chain Mode."
   set multiArg = "multi"

   foreach chain ($chainList)
      if (`string csame $chain None`) then
         echo "AddH: Fixing Null ChainIDs in Multi-Chain Input ..."
         fixPDBNull.tcl -in $inPDBName -out $adjPDBName -chain Auto -seg Null -verb
         set inPDBName = $adjPDBName
         echo ""
         break
      endif
   end
endif

mkdir $gmcName

if ($multiFlag == 0) then
   set chainList = NONE
else
   set chainList = (`getPDBInfo.tcl -in $inPDBName -noverb -info chainNames`)
endif

if ($#chainList == 1) then
   pdb2ac.tcl -dHi $dHi -in $1 -seg $segName -fc $fc -$centerFlag -verb > $gmcName/ac.tab

   pdb2gmc.tcl $dynArgs -$tkFlag -$delFlag \
               -in $inPDBName -inTitle $inPDBTitle -out $outPDBName \
               -seg $segName -gmc $gmcName -nohelix -bb $inPDBName  \
               -bbTitle $inPDBTitle -coolSteps $coolSteps -print $printStep \
               -cis -disu -noreset
else
   pdb2ac.tcl -dHi $dHi -in $1 -seg $segName -fc $fc -$centerFlag -multi -verb > $gmcName/ac.tab

   set ii = 1

   foreach chain ($chainList)
      if ($ii == 1) then
         pdb2gmc.tcl $dynArgs -$tkFlag -nodel \
                     -in $inPDBName -inTitle $inPDBTitle -out $outPDBName \
                     -chain $chain -seg $chain -gmc $gmcName -nohelix -bb $inPDBName  \
                     -bbTitle $inPDBTitle -coolSteps $coolSteps -print $printStep \
                     -cis -nodisu -noreset -cont -multi -nosa
      else if ($ii == $#chainList) then
         pdb2gmc.tcl $dynArgs -$tkFlag -$delFlag \
                     -in $inPDBName -inTitle $inPDBTitle -out $outPDBName \
                     -chain $chain -seg $chain -gmc $gmcName -nohelix -bb $inPDBName  \
                     -bbTitle $inPDBTitle -coolSteps $coolSteps -print $printStep \
                     -cis -disu -noreset -append -multi
      else
         pdb2gmc.tcl $dynArgs -$tkFlag -nodel \
                     -in $inPDBName -inTitle $inPDBTitle -out $outPDBName \
                     -chain $chain -seg $chain -gmc $gmcName -nohelix -bb $inPDBName  \
                     -bbTitle $inPDBTitle -coolSteps $coolSteps -print $printStep \
                     -cis -nodisu -noreset -cont -append -multi -nosa
      endif

      @ ii++
   end
endif

if (`flagLoc $argv -eval`) then
   set nanInfo = (`fgrep -i " NaN " dynExt.out  | wc`)
   set xList   = (`fgrep temperature dynExt.out | sed -e 's/temperature/ /' | tail -5`)
   set lastT   = (`fgrep temperature dynExt.out | sed -e 's/temperature/ /' | tail -1`)

   echo ""

   if ($#xList != 5) then
      nmrPrintf "AddH Warning: the Simulated Annealing did not complete as expected.\n"
      nmrPrintf "     The Proton Addition procedure might not be correct."
      set eFlag = 1
   else
      set avgT     = (`getStat.tcl -noverb -stat Avg -x $xList`)
      set goodFlag = (`IMATH "$avgT < 75.0"`)
  
      if ($nanInfo[1] != 0) then
         set goodFlag = 0
      endif
 
      if ($goodFlag == 0) then
         nmrPrintf "AddH Warning: the Simulated Annealing procedure did not converge.\n"
         nmrPrintf "     Average Temperature for Last 5 Steps: %.2f Final: %.2f\n"  $avgT $lastT
         nmrPrintf "     The Proton Addition procedure might not be correct."
         set eFlag = 2
      else
         nmrPrintf "AddH: Simulated Annealing Procedure Completed Successfully.\n"
         nmrPrintf "      Average Temperature for Last 5 Steps: %.2f Final: %.2f\n" $avgT $lastT
      endif
   endif

   echo ""

   adjTab.tcl -in $outPDBName -out $saPDBName -pdb -seg2chain -noverb

   set status = 0

   pdbDiff.tcl -ref $1 -refTitle $inPDBTitle -in $saPDBName -addH -$multiArg

   if ($status != 0) then
      echo ""
      echo "AddH Warning: Atoms in the Input Structure are not in the expected geometry."
      echo "     The Proton Addition procedure might not be ideal for all atoms."
      echo ""
      set eFlag = 3
   endif
endif

if ($restoreFlag != 0) then
   echo "AddH: Restoring Heavy Atom Coords to Original Positions."
   echo "AddH: Moving Single-Bonded Atoms Accordingly ..."

   set status = 0

   updateAC.tcl -gmc addH.gmc -in $outPDBName -out $outPDBName -$multiArg

   if ($status != 0) then
      echo ""
      echo "AddH: Warning: Expected Atoms May Be Missing from the Input Structure."
      echo "AddH: The Proton Addition procedure might not be ideal for all atoms."
      echo ""
      set eFlag = 4
   endif
endif

if ($multiFlag == 0) then
   adjTab.tcl -in $outPDBName -out $outPDBName -pdb -set SEGID -val "(Null)" -noverb
else
   adjTab.tcl -in $outPDBName -out $outPDBName -pdb -seg2chain -noverb
endif

if ($remFlag != 0 || $hetFlag != 0) then
   if (-e $tmpPDBName) then
      /bin/rm -rf $tmpPDBName
   endif

   if ($remFlag != 0) then
      echo "AddH: Restoring REMARKs from $inPDBTitle ..."
      splitPDB.tcl -in $outPDBName -rem stdout -out None        >> $tmpPDBName
      splitPDB.tcl -in $inPDBName  -rem stdout -out None -clean >> $tmpPDBName
      cat $outPDBName >> $tmpPDBName
      cat $tmpPDBName >  $outPDBName
   endif

   if ($hetFlag != 0) then
      echo "AddH: Restoring HETATMs from $inPDBTitle ..."
      splitPDB.tcl -in $inPDBName -hetatm stdout -out None -clean >> $outPDBName
      adjTab.tcl -in $outPDBName -out $outPDBName -pdb -renumber -hetatm -noverb 
   endif 
endif

exit $eFlag