# build solute in FF description from a PDB file

# author: mengen

#

# usage: ./modar step1_build_solute.inp pdb=<filename.pdb> [telltoptree=true]

#        telltoptree to print topology tree

#

if(! -var pdb) {

  echo -hsn "please tell the PDB file:"

  pdb=gettext()

  if(!-f $pdb) {

    echo -hs "no file $pdb found"

    abort

  }

}

if(! -var telltoptree) telltoptree=false

#load CHARMM force field

include "loadff.inc"

#save a copy of force filed in Modar format for reservation

#savefftop file="modar.mtf"

#saveffprm file="modar.prm"

#read molecular structure file and coordinate

tell nmodel file=$pdb

if($nmodel>1 && !-var imodel) {

  nmodel--

  getvar imodel default=0 min=0 max=$nmodel \

         msg="which model to use?"

}

if(!-var imodel) imodel=0

loadmsf fmt=pdb file=$pdb model=$imodel

#do PDB residue name to CHARMM FF residue name transform

include "do_pdb_nomination_replacement.inc"

#please put more nomination replacement here

#remove other atoms except the protein, RNA and DNA

#delete atom select="seg HET*"

#check unknown residue type and do replacement one by one if has

:tellunkr

tell unknownresiduetype

if($nresun>0) {

  echo -hs " there are $nresun unknown residue types"

  echo -hs " please replace following resname "

  echo -hs "   to a type declared by the force field"

  iunk=0

  while($iunk<$nresun){

    echo -hsn " rename $unresun[$iunk] to "

    newresn=gettext()

    if(!-z $newresn) {

       atom resn newname=$newresn select="resn $unresun[$iunk]"

    }

    iunk++

  }

  goto tellunkr #try again

}

tell unknownresiduetype

if($nresun>0) {

  echo -hs "** there are still $nresun unknown residue type for the force field" 

  echo -hs "** they are: $resun"

  echo -hs "** please check log file for more detail"

  abort

}

#build MSF according the force field topologies

build top

#put Protonation here, please check command "residue patch"

#make Disulfide bonds here if need

#Phosphorylation here if need

#put more patching or mutation here if need

#patch peptide, dna, rna terminals if has

include "do_protein_dna_rna_terminals_patch.inc"

#sort internal terms, including residues, bonds, angles ...

sortmsf

#tell mol top tree if request

if($telltoptree=="true") include "tell_top_tree.inc"

#save uninited coordinates atoms indexes for further minimization

defcls name="uninit" select="unset" outfile="nocoord.ndx"

nnocoord=$nasel

#build unknown coordinates automatically according to

#          the internal connections and FF parameters

if($nnocoord>0) {

  build unkcrd

}

#check msf

checkmsf

#check whether there are still some unknown coordinates

defcls name="testun" select="unset"

if($nasel>0)  {

  echo -hs "|  ** there are still $nasel atoms with coordinate uninitiated"

  echo -hs "|  ** please load a better initial structure file"

  abort

}

else{

  echo -hs "|  built unknown coordinates successfully"      

}

#do minimization for atoms which were with unknown coordinates

if($nnocoord>0) { #skip minimize

  nonbond type=cutoff nblcutoff=12.0 nbcutoff=10.0 swcutoff=8.0 eps=4.0

  doenergy

  atom fix select="! uninit"

  domin method=SDfuse nstep=200 nprint=10 ftol=1e-6

  domin method=SD     nstep=200 nprint=10 ftol=1e-6

}

else{

  echo -hs "|  there were no unknown coordinates, minimization is not needed"

}

#remove all fixes

atom fix select="none"

#tell solute information

tell msf_info

echo -hs "| "

echo -hs "|  Summary of solute:"

echo -hs "|    number of segments:    $nseg"

echo -hs "|    number of residues:    $nres"

echo -hs "|    number of atoms:       $natom"

echo -hs "|    number of bonds:       $nbond"

echo -hs "|    number of angles:      $nangl"

echo -hs "|    number of dihedrals:   $ndihe"

echo -hs "|    number of impropers:   $nimpr"

echo -hs "|    number of cross terms: $ncmap"

echo -hs "|    net charge:  $netcharge"

echo -hs "|    total mass:  $totmass"

#reorient

tell orient select="all"

atom reorient select="all"

#tell density

tell density

#translate to origin

atom translateto pos=(0,0,0) select="all"

#tell and save solute size

tell size select="all"

echo xext=$xext >  "solute_size.inc"

echo yext=$yext >> "solute_size.inc"

echo zext=$zext >> "solute_size.inc"

echo xmin=$xmin >> "solute_size.inc"

echo ymin=$ymin >> "solute_size.inc"

echo zmin=$zmin >> "solute_size.inc"

echo xmax=$xmax >> "solute_size.inc"

echo ymax=$ymax >> "solute_size.inc"

echo zmax=$zmax >> "solute_size.inc"

echo -hs "|    size in angstrom:"

echo -hs "|      X extension: $xext "

echo -hs "|      Y extension: $yext "

echo -hs "|      Z extension: $zext "

echo -hs "| "

#save molecular structure|topology,coordinates,pdb file respectively

savemsf file="step1_solute.msf"

savepdb file="step1_solute.pdb"

echo -hs "|  output files:"

echo -hs "|    solute_size.inc    brief size information of solute"

echo -hs "|    step1_solute.pdb   full structure in PDB format file"

echo -hs "|    step1_solute.msf   the topology file described by the Force Field"

echo -hs "| "

echo -hs "|  it's ready to go next step to prepare a water box by command:"

echo -hs "|    ./modar -nobak step2_gen_solvent_box.inp boxtype=cubic"

echo -hs "| "