remarks Solvent refinement protocol from ARIA1.2 (Nilges and Linge), modified for XPLOR-NIH
remarks Site the following references:

set message on echo on end

{*==========================================================================*}
{*====================== SET FILENAMES AND VARIABLES =======================*}
{*==========================================================================*}

{* type of non-bonded parameters *}
{* choice: "PROLSQ" "PARMALLH6" "PARALLHDG" "OPLSX" *}
{* The water refinement uses the OPLSX parameters *}

evaluate ( $par_nonbonded = "OPLSX" )

{* parameter file(s) *}
evaluate ( $par.1 = "TOPPAR:parallhdg5.3.pro" )
evaluate ( $par.2 = "TOPPAR:parallhdg5.3.sol" )
evaluate ( $par.3 = "" )
evaluate ( $par.4 = "" )
evaluate ( $par.5 = "" )

{* solvent topology file *}
evaluate ( $solvent_topology = "TOPPAR:topallhdg5.3.sol" )

{* structure file(s) *}
evaluate ( $struct.1 = "STRUCTURES:1xxx.psf" )
evaluate ( $struct.2 = "" )
evaluate ( $struct.3 = "" )
evaluate ( $struct.4 = "" )
evaluate ( $struct.5 = "" )

{* input coordinate file(s) *}
evaluate ( $pdb.in.file.1 = "INPUTCOORDINATES:analyzed_1.pdb" )

{* output coordinate file(s) *}
evaluate ( $pdb.out.file.1 = "OUTPUTCOORDINATES:refined_1.pdb" )

{* input distance restraint file(s) *}
evaluate ( $noe.file.1 = "TABLES:1xxx_noe.tbl" )

{* Averaging for NOE restraints *}
evaluate ( $noe.ave = sum )

{* input dihedral restraint file(s) *}
evaluate ( $cdih.file.1 = "TABLES:1xxx_dihe.tbl" )

{* seed for velocity generation *}
evaluate ( $seed = 12396 )

{* set number of md steps for the heating stage *}
evaluate ( $mdsteps.heat = 200 )

{* set number of md steps for the hot md stage *}
evaluate ( $mdsteps.hot = 2000 )

{* set number of md steps for the cooling stage *}
evaluate ( $mdsteps.cool = 200 )

{* seed for velocity generation *}
evaluate ( $seed = 12396 )

{* set acceptance criteria *}
evaluate ( $accept.noe = 0.50 )
evaluate ( $accept.cdih = 5.00 )
evaluate ( $accept.coup = 1.00 )
evaluate ( $accept.sani = 0.00 )
evaluate ( $accept.vean = 5.00 )


{*==========================================================================*}
{*=========== READ THE STRUCTURE, TOPOLOGY AND PARAMETER FILES =============*}
{*==========================================================================*}

structure @$struct.1 end
topology  @$solvent_topology end
parameter @$par.1 @$par.2 end

{*==========================================================================*}
{*================== SET VALUES FOR NONBONDED PARAMETERS ===================*}
{*==========================================================================*}

parameter
  nbonds
    nbxmod=5 atom cdiel shift
    cutnb=9.5 ctofnb=8.5 ctonnb=6.5 eps=1.0 e14fac=0.4 inhibit 0.25
    wmin=0.5
    tolerance  0.5
  end
end


{*==========================================================================*}
{*============== READ THE COORDINATES AND COPY TO REFERENCE SET ============*}
{*==========================================================================*}

coor @@$pdb.in.file.1

vector do (refx = x) (all)
vector do (refy = y) (all)
vector do (refz = z) (all)

{*==========================================================================*}
{*========================= GENERATE THE WATER LAYER =======================*}
{*==========================================================================*}

vector do (segid = "PROT") (segid "    ")
@SOLVENT:generate_water.inp
vector do (segid = "    ") (segid "PROT")

{*==========================================================================*}
{*========================= READ THE EXPERIMENTAL DATA =====================*}
{*==========================================================================*}

noe reset
  nrestraints = 10000
  ceiling     = 100

  @$noe.file.1
  averaging  * $noe.ave
  potential  * soft
  scale      * 50
  sqconstant * 1.0
  sqexponent * 2
  soexponent * 1
  rswitch    * 1.0
  sqoffset   * 0.0
  asymptote  * 2.0
end

restraints dihedral reset
   @$cdih.file.1
   scale = 200
end

#{ncs constraints for symmetric dimer}

#evaluate ($kncs=0.1)
#ncs restraints

#initialize
#group
#equi (resid 6:15 or resid 24:59 or resid 46:61 or resid 66:72 or resid 78:88)
#equi (resid 306:315 or resid 324:359 or resid 346:361 or resid 366:372 or resid 378:388)
#weight = $kncs
#end
#?
#end

{*==========================================================================*}
{*============================ SET FLAGS ===================================*}
{*==========================================================================*}

flags exclude *
      include bond angle dihe impr vdw elec
              noe cdih coup oneb carb ncs dani
              vean sani prot harm
end


{*==========================================================================*}
{*========================= START THE REFINEMENT ===========================*}
{*==========================================================================*}

set seed $seed end

! We loop untill we have an accepted structure, maximum trials=3
evaluate ($end_count = 3)
evaluate ($count = 0)

while ($count < $end_count ) loop main

  evaluate ($accept = 0)

  ! 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 ($kbonds = 1000)
  evaluate ($kangle = 50)
  evaluate ($kimpro = 5)
  evaluate ($kchira = 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
  constraints fix (not resn tip3) end
  minimize powell nstep=40 drop=100 end

  ! release protein and restrain harmonically
  constraints fix (not all) end
  vector do (refx=x) (all)
  vector do (refy=y) (all)
  vector do (refz=z) (all)
  restraints harmonic
     exponent = 2
  end
  vector do (harmonic = 0)  (all)
  vector do (harmonic = 10) (not name h*)
  vector do (harmonic = 20.0)(resname ANI and name OO)
  vector do (harmonic = 0.0) (resname ANI and name Z )
  vector do (harmonic = 0.0) (resname ANI and name X )
  vector do (harmonic = 0.0) (resname ANI and name Y )

  constraints
    interaction (not resname ANI) (not resname ANI)
    interaction ( resname ANI) ( resname ANI)
  end

  minimize powell nstep=40 drop=10 end
  vector do (refx=x) (not resname ANI)
  vector do (refy=y) (not resname ANI)
  vector do (refz=z) (not resname ANI)
  minimize powell nstep=40 drop=10 end
  vector do (refx=x) (not resname ANI)
  vector do (refy=y) (not resname ANI)
  vector do (refz=z) (not resname ANI)

  vector do (mass =50) (all)
  vector do (mass=1000) (resname ani)
  vector do (fbeta = 0) (all)
  vector do (fbeta = 20. {1/ps} ) (not resn ani)
  evaluate ($kharm = 50)
  ! heat to 500 K
  for $bath in (100 200 300 400 500 ) loop heat
     vector do (harm = $kharm) (not name h* and not resname ANI)
     vector do (vx=maxwell($bath)) (all)
     vector do (vy=maxwell($bath)) (all)
     vector do (vz=maxwell($bath)) (all)
     dynamics verlet
        {time step was 3 fs}
        nstep=$mdsteps.heat timest=0.003{ps}
        tbath=$bath  tcoupling = true
        iasvel=current
        nprint=50
     end
     evaluate ($kharm = max(0, $kharm - 4))
     vector do (refx=x) (not resname ANI)
     vector do (refy=y) (not resname ANI)
     vector do (refz=z) (not resname ANI)
  end loop heat


  
  ! refinement at high T:
  constraints
    interaction (not resname ANI) (not resname ANI) weights * 1 dihed 3 end
    interaction ( resname ANI) ( resname ANI) weights * 1 end
  end

  vector do (harm = 0)  (not resname ANI)
  vector do (vx=maxwell($bath)) (all)
  vector do (vy=maxwell($bath)) (all)
  vector do (vz=maxwell($bath)) (all)
  dynamics verlet
     nstep=$mdsteps.hot timest=0.004 {ps}
     tbath=$bath  tcoupling = true
     iasvel=current
     nprint=50
     !trajectory=1xxx_hot.dat nsavc=5
  end

  constraints
    interaction (not resname ANI) (not resname ANI) weights * 1 dihed 5 end
    interaction ( resname ANI) ( resname ANI) weights * 1  end
  end


  ! cool
  evaluate ($bath = 500)
  while ($bath >= 25) loop cool

     evaluate ($kbonds    = max(250,$kbonds / 1.1))
     evaluate ($kangle    = min(250,$kangle * 1.1)) ! equi-ed after 17 steps
     evaluate ($kimpro    = min(300,$kimpro * 1.4)) ! 12     
     evaluate ($kchira    = min(1000,$kchira * 1.4))! 16
     evaluate ($komega    = min(100,$komega * 1.4)) !  9

     parameter
        bond     (not resn tip3 and not name H*)(not resn tip3 and not name H*) $kbonds  TOKEN
        angle    (not resn tip3 and not name H*)(not resn tip3 and not name H*)(not resn tip3 and not name H*) $kangle  TOKEN
        improper (all)(all)(all)(all) $kimpro  TOKEN TOKEN
        !VAL: stereo CB
        improper (name HB and resn VAL)(name CA and resn VAL)(name CG1 and resn VAL)(name CG2 and resn VAL) $kchira TOKEN TOKEN
        !THR: stereo CB
        improper (name HB and resn THR)(name CA and resn THR)(name OG1 and resn THR)(name CG2 and resn THR) $kchira TOKEN TOKEN
        !LEU: stereo CG
        improper (name HG and resn LEU)(name CB and resn LEU)(name CD1 and resn LEU)(name CD2 and resn LEU) $kchira TOKEN TOKEN
        !ILE: chirality CB
        improper (name HB and resn ILE)(name CA and resn ILE)(name CG2 and resn ILE)(name CG1 and resn ILE) $kchira TOKEN TOKEN
        !chirality CA
        improper (name HA)(name N)(name C)(name CB) $kchira TOKEN TOKEN

        improper (name O)  (name C) (name N) (name CA) $komega TOKEN TOKEN
        improper (name HN) (name N) (name C) (name CA) $komega TOKEN TOKEN
        improper (name CA) (name C) (name N) (name CA) $komega TOKEN TOKEN
        improper (name CD) (name N) (name C) (name CA) $komega TOKEN TOKEN
     end

     vector do (vx=maxwell($bath)) (all)
     vector do (vy=maxwell($bath)) (all)
     vector do (vz=maxwell($bath)) (all)
     dynamics verlet
        !was a time step of 0.004
        nstep=$mdsteps.cool timest=0.004{ps}
        tbath=$bath  tcoupling = true
        iasvel=current
        nprint=50
     end

     evaluate ($bath = $bath - 25)
  end loop cool


  !final minimization:
  mini powell nstep 200 end

  {*==========================================================================*}
  {*======================= CHECK RESTRAINT VIOLATIONS =======================*}
  {*==========================================================================*}

  print threshold=$accept.noe noe
  evaluate ($v_noe = $violations)
  print threshold=0.5 noe
  evaluate ($v_noe_0.5 = $violations)
  print threshold=0.4 noe
  evaluate ($v_noe_0.4 = $violations)
  print threshold=0.3 noe
  evaluate ($v_noe_0.3 = $violations)
  print threshold=0.2 noe
  evaluate ($v_noe_0.2 = $violations)
  print threshold=0.1 noe
  evaluate ($v_noe_0.1 = $violations)
  evaluate ($rms_noe = $result)

  print threshold=$accept.cdih cdih
  evaluate ($rms_cdih=$result)
  evaluate ($v_cdih = $violations)

  coupl print thres = $accept.coup class c1 end
  evaluate ($rms_coup = $result)
  evaluate ($v_coup = $violations)

  sani print threshold = $accept.sani class rdc1 end
  evaluate ($rms_sani = $result)
  evaluate ($v_sani = $violations)

  vean print threshold = $accept.vean class vea1 end
  evaluate( $rms_vean = $result)
  evaluate( $v_vean = $violations)

  print thres=0.05 bonds
  evaluate ($rms_bonds=$result)
  evaluate ($v_bonds = $violations)
  print thres=5. angles
  evaluate ($rms_angles=$result)
  evaluate ($v_angles = $violations)
  print thres=5. impropers
  evaluate ($rms_impropers=$result)
  evaluate ($v_impropers = $violations)

  if ($v_noe > 0) then  evaluate ( $accept=$accept + 1 ) end if
  if ($v_cdih > 0) then  evaluate ( $accept=$accept + 1 ) end if
  if ($v_coup > 0) then  evaluate ( $accept=$accept + 1 ) end if
  if ($v_sani > 0) then  evaluate ( $accept=$accept + 1 ) end if
  if ($v_vean > 0) then  evaluate ( $accept=$accept + 1 ) end if

  if ($accept = 0 ) then
    evaluate ( $label = "ACCEPTED" )
    exit main
  else
    evaluate ( $label = "NOT ACCEPTED" )
    evaluate ( $count = $count + 1 )
  end if

end loop main

{*==========================================================================*}
{*================ COMPUTE RMS DIFFERENCE BETWEEN OBSERVED =================*}
{*========================= AND MODEL DISTANCES. ===========================*}
{*==========================================================================*}

constraints interaction
   (not resname TIP* and not resname ANI)
   (not resname TIP* and not resname ANI)
end

energy end

remarks Structure $label
remarks E-overall:                             $ener
remarks E-NOE_restraints:                      $noe
remarks E-CDIH_restraints:                     $cdih
remarks E-COUP_restraints:                     $coup
remarks E-SANI_restraints:                     $sani
remarks E-VEAN_restraints:                     $vean
remarks RMS-NOE_restraints:                    $rms_noe
remarks RMS-CDIH_restraints:                   $rms_cdih
remarks RMS-COUP_restraints:                   $rms_coup
remarks RMS-SANI_restraints:                   $rms_sani
remarks RMS-VEAN_restraints:                   $rms_vean
remarks NOE Acceptance criterium:              $accept.noe
remarks NOE violations > $accept.noe:          $v_noe
remarks All NOE viol. (>0.5,0.4,0.3,0.2,0.1A): $v_noe_0.5 $v_noe_0.4 $v_noe_0.3 $v_noe_0.2 $v_noe_0.1
remarks CDIH Acceptance criterium:             $accept.cdih
remarks CDIH violations > $accept.cdih:        $v_cdih
remarks COUP Acceptance criterium:             $accept.coup
remarks COUP violations > $accept.coup:        $v_coup
remarks SANI Acceptance criterium:             $accept.sani
remarks SANI violations > $accept.sani:        $v_sani
remarks VEAN Acceptance criterium:             $accept.vean
remarks VEAN violations > $accept.vean:        $v_vean

write coordinates sele= (not resn TIP3) output =$pdb.out.file.1 end
!write coordinates sele= (all) output =$pdb.out.file.1 end

stop