!$Revision: 2.6 $
!$Date: 2002/07/23 16:19:27 $
!$RCSfile: re_h2o.inp,v $

{* ==========================================================================  *}
remarks  re_h2o.inp (derived from sa_l.inp)
remarks  Author: Michael Nilges, Jens Linge 26.8.98
remarks  seperated from ARIA1.2 and modified by Jinwon Jung 2003.2.25
{* ==========================================================================  *}

{*  ============ things you want to change ================= *}
evaluate ($seed = 498037388 )                {* seed for random generator        *}
evaluate ($count = 4 )                 {* starting structure number        *}
evaluate ($strucfile = "sgr4_h2o.mtf")  {* molec. topology file             *}
evaluate ($Noeavg = sum)                   {* Your choice of center, sum       *}
evaluate ($xtimes = 1)                {* apply user selected scale        *}
evaluate ($miparm = "PARAM19")             {* choice of OPLSX|PARAM19|PARMALLH6
                                                        PROLSQ|CONTACT         *}
evaluate ($heatingc = "200")
evaluate ($hotcycle = "1000")
evaluate ($coolcycl = "100")

{* ============= Things you _may_ want to change ============ *}
evaluate ($rootname = "cnsPDB/sa_cns_")
evaluate ($outname = "refinedPDB/resa_")
evaluate ($outwatr = "refinedPDB_w/resa_")
evaluate ($finomega = 100)
evaluate ($wpdb = "TOPOWAT:boxtyp20.pdb")    {* For generate_water.cns         *}

evaluate ($HaveNoe = "yes")                   ! change to "no" if no NOE data
evaluate ($HaveHbond = "no")
evaluate ($HaveDih = "yes")                   ! change to "no" if no DIH data
evaluate ($noe_rstrs = "sgr4_noe.tbl")    ! File holding the NOE data
evaluate ($hbn_rstrs = "sgr4_hbond.tbl")
evaluate ($dih_rstrs = "sgr4_dihe.tbl")   ! File holding DIHEDRAL data

{******************* not yet developed
evaluate ($HaveJcoup = "no")
evaluate ($filejcoup = "sgr4_Jhnha.tbl")
evaluate ($HaveRDC = "no")
evaluate ($filesani = "sgr4_sani.tbl")    ! file with RDC restraints
evaluate ($axispdb = "axis.pdb")             ! file with axis COORDS
evaluate ($HaveCaCbShifts = "no")
evaluate ($fileshifts_CaCb = "sgr4_shiftsCACB.tbl")
**********************}

{* =============== things yo don't want to change ================ *}
{* -- from here down you need to change nothing unless you know -- *}
{* -- what you are doing .......                                -- *}
{* ====  protocol starts ==== *}

set seed $seed end
set abort normal end
topology
    @@TOPOWAT:topallhdg5.3.pro
    @@TOPOWAT:topallhdg5.3.sol
    @@TOPOWAT:ion.top
end

! read coordinate and copy to reference coordinate set

! the water refinement uses a full Lenard-Jones potential:
evaluate ($par_nonbonded = $miparm)
!!evaluate ($par_nonbonded = $toppar.par_nonbonded)

!read the parameter files:
PARAMETER
      @@TOPOWAT:parallhdg5.3.pro
      @@TOPOWAT:parallhdg5.3.sol
      @@TOPOWAT:ion.param

      nbonds   {* load a few defaults, main load is done below *}
          nbxmod=5 tolerance 0.5 wmin=1.5
      end
end

{* -- adjust NBONds depending on param being used -- *}
{*
   -- depending on the choice of paramaters to be used here the 
   -- most advised selection for nonbonded parameters is done.
*}
     if ( $miparm = "OPLSX" ) then
        parameter
           nbonds
              cutnb=9.5  ctofnb=8.5   ctonnb=6.5 wmin=1.5
              nbxmod=5   atom cdiel shift
              repel=0.0  tolerance 0.5 
              eps=1.0    e14fac=0.4   inhibit 0.25
           end
        end {* this to end parameter *}
     {---------------- PROLSQ ---------------------}
     elseif ( $miparm = "PROLSQ" ) then
        parameter
           nbonds
              cutnb=9.5  ctofnb=8.5   ctonnb=6.5 wmin=1.5
              nbxmod=5   
              repel=1.0  tolerance 0.5 
	      rexponent=4 irexponent=1 rconst=16.0
           end
        end
     {-----------------PARAM19---------------------}
     elseif ( $miparm = "PARAM19" ) then
        parameter
           nbonds
              cutnb=9.5  ctofnb=8.5   ctonnb=6.5 wmin=1.5
              nbxmod=5   atom cdiel shift
              repel=0.0  tolerance 0.5 
           end
        end
     {----------------PARMALLH6--------------------}
     elseif ( $miparm = "PARMALLH6" ) then
        parameter
           nbonds
              cutnb=9.5  ctofnb=8.5   ctonnb=6.5 wmin=1.5
              nbxmod=5   
              repel=0.8  tolerance 0.5 
	      rexponent=2 irexponent=2 rconst=5.0 
           end
        end
     {----------------   CONTACT  -----------------}
     elseif ( $miparm = "CONTACT" ) then
        parameter
           nbonds
              cutnb=7.0  ctofnb=6.0   ctonnb=5.5 wmin=1.5
              nbxmod=5   
              repel=1.0  tolerance 0.5 
	      rexponent=4 irexponent=1 rconst=16.0
           end
        end
     end if

{* before it was a cycle, now it is done only once *}

     evaluate ($pdbfile = $rootname+encode($count)+".pdb")
     structure reset end 
     structure @@$strucfile end

     evaluate ($HaveCis = "yes")    {* -- Apply possible CIS peptide patches -- *}
     if ( $HaveCis = "yes" ) then
         !CISpep info
         patch CISP                 
              reference=NIL=(resid 42) 
         end                          
     end if

     evaluate ($HaveHisd = "no")    {* -- Apply possible HISD peptide patches -- *}
     if ( $HaveHisd = "yes" ) then
         !HISDpep info
     end if

     evaluate ($HaveHise = "no")    {* -- Apply possible HISE peptide patches -- *}
     if ( $HaveHise = "yes" ) then
         !HISEpep info
     end if

     evaluate ($HaveDisu = "no")        {* -- Getting ready for S-S bridges -- *}
     if ( $HaveDisu = "yes" ) then
         !SSBridge  info
     end if

     coor init end
     coor @@$pdbfile
     do (refx = x) (all)
     do (refy = y) (all)
     do (refz = z) (all)

     ! generate water layer
     do (segid = "PROT") (segid "    ")
     @TOPOWAT:generate_water.cns
     do (segid = "    ") (segid "PROT")

     flags include bonds angle vdw impro cdih dihe noe elec end

     {* -- write the coordinates with water to see what's going on -- *}
     evaluate ($waternam0 = $outwatr+encode($count)+"_waterIni.pdb")
     write coordinates output = $waternam0 end 

     noe                                 {* -- NOE section -- *}   
          reset 
          ceiling 1000
          nrestraints 10000
     end

     if ( $HaveNoe = "yes" ) then        {* -- load only if present -- *}
         noe class dist @@$noe_rstrs end
     end if
 
     if ( $HaveHbond = "yes" ) then      {* -- load only if present -- *} 
         noe class hbond @@$hbn_rstrs end 
     end if 
         
     noe
          averaging  * $Noeavg
          potential  * soft
          scale      * $xtimes
          sqconstant * 1.0
          sqexponent * 2
          soexponent * 1
          rswitch    * 1.0
          sqoffset   * 0.0
          asymptote  * 2.0
          msoexponent * 1
          masymptote  * -0.1
          mrswitch    * 1.0
          avexpo hbond 20
     end

     if ( $HaveDisu = "yes" ) then
         !SSBondsAsNoe                   {* -- Add SSbonds as noe restraints *** }
     end if
      
     if ( $HaveDih = "yes" ) then 
          restraints dihedral            {* -- Read dihedral constraints     -- *}
               @@$dih_rstrs
          end
     end if

     !!!Values taken from run.cns
     evaluate ($k_unamb= 50 )
     evaluate ($k_amb= 50 )
     evaluate ($k_hbond= 50 )

     {* -- evaluate scales values -- *}
     evaluate ($scalambi =  50 * $xtimes)
     evaluate ($scaldist =  50 * $xtimes)
     evaluate ($scalhbnd = 100 * $xtimes)
     evaluate ($scaldihd = 200 * $xtimes)

     if ( $HaveNoe = "yes" ) then 
         noe                                {* -- scale and configure NOE term -- *}
            rswitch ambi 0.5 
            rswitch dist 0.5
            mrswitch ambi 0.5
            mrswitch dist 0.5
            asym ambi 0.1
            asym dist 0.1
            masym ambi -0.1
            masym dist -0.1
            scale ambi  $scalambi 
            scale dist  $scaldist
	 end 
   end if 

   if ( $HaveHbond = "yes" ) then         {* -- scale and configure for HBONS -- *} 
         noe
            rswitch hbon 0.5
            mrswitch hbon 0.5
            asym hbon 0.1
            masym hbon -0.1
            scale hbond $scalhbnd 
         end
    end if

    if ( $HaveDih = "yes" ) then 
         restraints dihedral               {* -- scale dihedral term -- *}
              scale=$scaldihd
         end
    end if 

    ! since we do not use SHAKe, increase the water bond angle energy constant
    parameter 
            angle (resn tip3) (resn tip3) (resn tip3) 500 TOKEN 
    end

    ! reduce improper and angle force constant for some atoms
    evaluate ($kangle = 50)
    evaluate ($kimpro = 5)
    evaluate ($komega = 5)
    parameter
          angle    (not resn tip3)(not resn tip3)(not resn tip3) $kangle  TOKEN
          improper (all)(all)(all)(all) $kimpro  TOKEN TOKEN
    end

    ! fix the protein for initial minimization
    fix sele = (not resn tip3) end
    minimize powell nstep=40 drop=100 end

    ! release protein and restrain harmonically
    fix sele = (not all) end
    do (refx=x) (all)
    do (refy=y) (all)
    do (refz=z) (all)
    restraints harmonic 
          exponent = 2
    end
    do (harm = 0)  (all)
    do (harm = 10) (not name h*)
    igroup
         interaction (all) (all)
    end

    minimize powell nstep=40 drop=10 end
    do (refx=x) (all)
    do (refy=y) (all)
    do (refz=z) (all)

    minimize powell nstep=40 drop=10 end
    do (refx=x) (all)
    do (refy=y) (all)
    do (refz=z) (all)
    do (mass = 100) (all)
    do (fbeta = 0) (all)
    do (fbeta = 20. {1/ps} ) (all)                
    evaluate ($kharm = 50) 

    {* -- heat to 500 K  -- *}
    for $bath in (100 200 300 400 500) loop heat
          do (harm = $kharm) (not name h* )
          do (vx=maxwell($bath)) (all)
          do (vy=maxwell($bath)) (all)
          do (vz=maxwell($bath)) (all)  
          dynamics cartesian
               nstep=$heatingc timest=0.003 {ps}       
               temperature=$bath  tcoupling = true
               nprint=50 
          end 
          evaluate ($kharm = max(0, $kharm - 4))
          do (refx=x) (all)
          do (refy=y) (all)
          do (refz=z) (all)
    end loop heat

    {* -- refinement at high T: -- *}
    igroup
         interaction (all) (all) weights * 1 dihed 2 end
    end

    do (harm = 0)  (all)
    do (vx=maxwell($bath)) (all)
    do (vy=maxwell($bath)) (all)
    do (vz=maxwell($bath)) (all)  
    dynamics cartesian
            nstep=$hotcycle timest=0.004{ps}      
            temperature=$bath  tcoupling = true
            nprint=50 
    end 

    igroup
           interaction (all) (all) weights * 1 dihed 3 end
    end

    {* -- cool -- *}
    evaluate ($bath = 500)
    while ($bath ge 25) loop cool
           evaluate ($kangle = min(500,$kangle * 1.2))
           evaluate ($kimpro = min(500,$kimpro * 1.5))
           evaluate ($komega = min($finomega,$komega * 1.5))
           parameter
               angle    (not resn tip3)(not resn tip3)(not resn tip3) $kangle  TOKEN
               improper (all)(all)(all)(all) $kimpro  TOKEN TOKEN
               improper (name CA) (name C) (name N) (name CA) $komega TOKEN TOKEN
           end
           do (vx=maxwell($bath)) (all)
           do (vy=maxwell($bath)) (all)
           do (vz=maxwell($bath)) (all)
           dynamics cartesian
                nstep=$coolcycl timest=0.004{ps}      
                temperature=$bath  tcoupling = true                       
                nprint=50 
           end 
           evaluate ($bath = $bath - 25)
    end loop cool

    {* -- final minimization: -- *}
    mini powell nstep 200 end
    igroup interaction 
           (not resname TIP* ) 
           (not resname TIP* ) 
    end

    {* -- prepare names for output files -- *}
    evaluate ($filename= $outname+encode($count)+".pdb")
    evaluate ($waternam= $outwatr+encode($count)+"_waterEnd.pdb")
    evaluate ($viofile= $outname+encode($count)+".vio")
      
    set display_file=$viofile end
    set print_file=$viofile end
    energy end
    @TOPOWAT:print_coorheader.cns
    display overall,bonds,angles,improper,dihe,vdw,elec,noe,cdih,coup,sani,vean
    display energies: $ener, $bond, $angl, $impr, $dihe, $vdw, $elec, $noe, $cdih, $coup, $sani, $vean
    do (q=1) (all)
    write coordinates sele= (not resn TIP3) output = $filename end

    {* -- write the coordinates with water to see what's going on -- *}
    write coordinates output = $waternam end 

    delete sele = (resn TIP3) end

stop