# YASARA MACRO
# TOPIC:       5. Structure Determination
# TITLE:       Refine an NMR ensemble in explicit solvent
# REQUIRES:    YASARA Structure and NMR Structure Determination Module
# AUTHOR:      Elmar Krieger, Sander Nabuurs, Chris Spronk
# LICENSE:     GPL
# DESCRIPTION: This macro refines an NMR ensemble in explicit solvent to increase the accuracy

# To refine an ensemble in water, you need:
# - An ensemble of correctly folded structures in PDB format
#   default filenames are 'ensemble001.pdb', 'ensemble002.pdb' etc.
# - A file with distance and dihedral angle restraints,
#   default filename is 'restraints.tbl'
#
# You can either set the target by clicking on Options > Macro > Set target,
# or by uncommenting the line below and specifying it directly.
#macrotarget = 'c:\MyProject\1crn'

# If you want to change the defaults, do it here, before the defaults are included.
# (look in nmr_setdefaults to see what you can change)

# Set those defaults that have not been set already
include nmr_setdefaults

# Step 3: Trade speed for accuracy and use the best YASARA got in the final stage
# ===============================================================================
# This means explicit solvent, Particle Mesh Ewald and the YASARA force field.
RestrainPar (defaultpar)
SimSteps 5
Console Off
LogAs (logfile),append=Yes, print 'Explicit solvent refinement:'
for i=start to start+structures-1
  Clear
  # Load the ensemble member, maybe the user just started here 
  size = FileSize (ensemblefile)(i)
  if size
    LoadPDB (ensemblefile)(i)
  else
    LoadPDB (ensemblefile)(start)
  # Correction of naming deviations (Val/Leu) could be useful, but Clean may do too much
  #Clean
  LogAs (logfile),append=Yes, print '  Ensemble structure (i):'
  # Choose a nice orientation
  # (this aligns the major protein axes with the global coordinate system
  #  and makes sure the simulation cell gets as small as possible).
  NiceOriObj 1
  # Select force field
  ForceField (fofwater)
  ScaleKBP 7,7,7
  ScaleForce Angle,2
  ScaleForce Planarity,5
  Cutoff 7.86
  Boundary periodic
  Longrange Coulomb
  # Refinement at room temperature
  Temp (298+i)K
  # Create the cell
  Cell Auto,Extension=5
  
  # Check if there is a maximum allowed violation
  if violdismax!=9999 or violdihmax!=9999
    # Iteratively pull stronger and stronger until requirements are met
    RestrainPot (defaultpot)
    LoadTbl (restrainfile),1
    ScaleRest (defaultscale)
    TempCtrl Anneal
    Timestep 2,0.25
    scaledis=25
    scaledih=2
    Sim On
    do
      # Force as much as necessary
      ScaleRest Distance=(scaledis),Dihedral=(scaledih)
      Wait 10
      scaledis=scaledis+10
      scaledih=scaledih+1
      dissum,dismax,dihsum,dihmax=RestViol
      ShowMessage 'Forcing with ScaleDis=(0.00+scaledis) and ScaleDih=(0.00+scaledih), ViolDisMax=(0.00+dismax), ViolDihMax=(0.00+dihmax)'
    while scaledis<1000 and (dismax>violdismax or dihmax>violdihmax)
    DelRest
  
  # Fill with water, predict pKas, place counter ions
  Experiment Neutralization
    WaterDensity 0.997
    pH (ph)
    Speed Fast
  Experiment On
  Wait ExpEnd
  # At this point, NMR structures are still sometimes loosely packed, and water
  # molecules may have ended up in the core. So we delete all waters that are
  # relatively lonely, because nothing is worse than water refinement with
  # water inside the protein.
  DelRes HOH with <7 atoms closer 5 among HOH Atom O
  # Load the restraints
  RestrainPot (defaultpot)
  LoadTbl (restrainfile),1
  ScaleRest (defaultscale)
  # Start the Simulation
  TimeStep 2,1.0
  TempCtrl SteepDes
  Sim On
  ShowMessage 'Starting explicit solvent refinement of structure (i) with a steepest descent minimization...'
  for j=1 to 1000
    Wait 1
    if SpeedMax<3000
      break
  # Correct cis-peptide bonds and wrong isomers also by remodeling
  CorrectCis On,Proline=(correctcispro)
  CorrectIso On
  CorrectConv On
  bestenergy=1e100
  bestenergyshown='not yet determined'
  # Store current SqOffset in 'offset'
  x,x,offset=RestrainPot
  failures=0
  failuresmax=4
  while failures<failuresmax
    Style Ribbon
    # Set strong restraining forces
    RestrainPot SqOffset=(0.5+offset)
    ScaleRest (strongscale)
    # Heat
    ShowMessage 'Simulated annealing refinement in solvent, cycle (failures)/(failuresmax). Minimum violation energy is (bestenergyshown). Heating...'
    Temp (298+i)K
    TempCtrl Rescale
    Wait (watersimsteps)
    # And cool
    ShowMessage 'Simulated annealing refinement in solvent, cycle (failures)/(failuresmax). Minimum violation energy is (bestenergyshown). Cooling...'
    TempCtrl Anneal
    Wait 20
    # Back to default restraining forces
    RestrainPot (defaultpot)
    ScaleRest (defaultscale)
    Wait 20
    # See if the restraint violations are smaller
    violenergy=RestEnergy
    LogAs (logfile),append=Yes, print '    Cycle (failures)/(failuresmax): Violation energy is (0.00+violenergy) (EnergyUnit)'
    if violenergy<bestenergy
      # Yes, replace the structure. Cannot use X-PLOR format because of missing CONECTs 
      bestenergy=violenergy
      SavePDB 1,(ensemblefile)(i)
      LogAs (logfile),append=Yes, print '      Best energy so far, this is the new member of the final ensemble.'
      failures=0
    else
      failures=failures+1
    bestenergyshown=0.00+bestenergy

  # Finally, force restraints again
  if violdismax!=9999 or violdihmax!=9999
    Clear
    LoadPDB (ensemblefile)(i)
    # Correction of naming deviations (Val/Leu) could be useful, but Clean may do too much
    #Clean
    Cell Auto,Extension=5
    # Iteratively pull stronger and stronger until requirements are met
    RestrainPot (defaultpot)
    LoadTbl (restrainfile),1
    ScaleRest (defaultscale)
    TempCtrl Anneal
    Timestep 2,0.25
    scaledis=25
    scaledih=2
    Sim On
    do
      # Force as much as necessary
      ScaleRest Distance=(scaledis),Dihedral=(scaledih)
      Wait 10
      scaledis=scaledis+10
      scaledih=scaledih+1
      dissum,dismax,dihsum,dihmax=RestViol
      ShowMessage 'Forcing with ScaleDis=(0.00+scaledis) and ScaleDih=(0.00+scaledih), ViolDisMax=(0.00+dismax), ViolDihMax=(0.00+dihmax)'
    while scaledis<1000 and (dismax>violdismax or dihmax>violdihmax)
    RestrainPar (defaultpar)
    SavePDB 1,(ensemblefile)(i)