# YASARA MACRO
# TOPIC:       5. Structure Determination
# TITLE:       Refine an NMR ensemble in vacuo
# REQUIRES:    YASARA Structure and NMR Structure Determination Module
# AUTHOR:      Elmar Krieger, Sander Nabuurs, Chris Spronk
# LICENSE:     GPL
# DESCRIPTION: This macro refines a roughly folded structure in vacuo to create a realistic fold

# To refine an ensemble in vacuo, you need:
# - An ensemble of roughly 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 2: Proceed to a realistic structure by in vacuo simulation with NOVA
# =========================================================================
# Having created the approximate protein fold, the next step is to
# increase the resolution and let the hydrogen bonds snap into place.
RestrainPar (defaultpar)
SimSteps 10
Console Off
LogAs (logfile),append=yes,print 'In vacuo 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):'
  # This is done in vacuo with YASARA. At this stage, speed is more important
  # than accuracy, so a very small 5.24 A cutoff is used.
  ForceField YASARA
  Cutoff 5.24
  TimeStep 2,0.75
  Cell Auto,Extension=20,Scale=Yes
  # Load the restraints
  RestrainPot (defaultpot)
  LoadTbl (restrainfile),1
  ScaleRest (defaultscale)
  # Exclude non-bonded interactions for now
  Interactions Bond,Angle,Dihedral,Planarity
  Temp (500+i*10)K
  # To avoid simulation failures due to excess restraining forces,
  # brake every 200th fast atom
  Brake (Atoms/200)
  # Correct cis-peptide bonds and wrong isomers by applying a force
  CorrectCis byForce,Proline=(correctcispro)
  CorrectIso byForce
  CorrectConv On
  Sim On
  ShowMessage 'Initial steepest descent minimization...'
  TempCtrl SteepDes
  for j=1 to 1000
    Wait 1
    if SpeedMax<4000
      break
  # Shake it a bit
  ShowMessage 'Relaxing structure without non-bonded interactions...'
  TempCtrl Reassign
  Wait 20
  # Now comes the 'trick': the non-bonded interactions are slowly added,
  # which causes the protein to blow up to a realistic volume. The first
  # time, this may trap part of the peptide chain at the wrong spot,
  # that's why the cycle is repeated until not improvement is found for
  # 'failuresmax' times.
  failuresmax=4
  # To avoid loose structures, we don't use a 'radius of gyration' refinement,
  # but instead start pulling a bit earlier
  _,_,offset=RestrainPot
  RestrainPot SqOffset=(0.5+offset)

  # First everything, then only the sidechains
  for fixed in 'none','backbone'
    FixAtom (fixed)
    bestenergy=1e100
    bestenergyshown='not yet determined'
    failures=0
    while failures<failuresmax
      # Slowly add the non-bonded forces in 10 steps
      # Because the influence of non-bonded interactions is larger in the
      # beginning, more time is spent there (defined by 'waitlist')
      nbilist= 0.0, 0.0005, 0.01, 0.05, 0.1, 0.2, 0.3, 0.5, 0.7, 1.0
      waitlist=50,  30,     20,   20,    10,  10,  5,   5,   5,   30
      for k=1 to 10
        ShowMessage 'Scaling non-bonded interactions with (nbilist(k)), (k*10)% of try (failures+1)/(failuresmax) completed. Minimum violation energy is (bestenergyshown).'
        if k==1
          # First step, non-bonded interactions are disabled
          Interactions Bond,Angle,Dihedral,Planarity
        else
          # In all subsequent steps, the steepest descent is done only with
          # repulsive VdW forces, otherwise the electrostatic attraction of
          # close atoms might cause 'nuclear fusion'
          Interactions Bond,Angle,Dihedral,Planarity,VdW
          # Also, we don't want restraints at this point
          ScaleRest Distance=0,Dihedral=0
        ScaleForce NonBonded,(nbilist(k))
        # Short steepest descent to remove ugly bumps
        TempCtrl SteepDes
        for l=1 to 1000
          Wait 1
          if SpeedMax<4000
            break
        # Bring back the strong scaling
        ScaleRest (strongscale)
        if k>1
          Interactions Bond,Angle,Dihedral,Planarity,VdW,Coulomb
        # Short MD
        TempCtrl Reassign
        Wait (waitlist(k))
        # Show secondary structure changes
        Style Ribbon,Stick
      violenergy=RestEnergy
      if violenergy<bestenergy
        # Yes, replace the structure
        bestenergy=violenergy
        bestenergyshown='now (0.00+bestenergy)'
        SavePosAtom all,bestpos
        failures=0
      else
        failures=failures+1
        bestenergyshown='still (0.00+bestenergy)'
      LogAs (logfile),append=Yes, print '    Violation energy is (0.00+violenergy), best so far (0.00+bestenergy) (EnergyUnit), (failures) failures recently.'
    LoadPosAtom all,bestpos
  Free
  
  # In the next stage, we run a few heating-cooling cycles to get a structure that
  # optimially fits the experimental restraints.
  cycles=40
  # Choose the standard orientation just before switching the simulation off and on
  # (otherwise mouse rotations could lead to rounding differences - and different trajectories)
  Ori 0,0,0
  Cutoff 10.48
  # Correct cis-peptide bonds and wrong isomers also by remodeling
  CorrectCis On,Proline=(correctcispro)
  CorrectIso On
  CorrectConv On
  # Create cysteine bridges
  if cysbridgelist(1)=='Auto'
    # Automatic assignment
    LinkCys DelH=Yes
  elif cysbridgelist(1)!='None'
    # Linked cysteines have been specified explicitly
    for j=1 to count cysbridgelist step 2
      DelAtom HG Res (cysbridgelist(j)) (cysbridgelist(j+1))
      AddBond SG Res (cysbridgelist(j)),SG Res (cysbridgelist(j+1))
  # Correction of naming deviations (Val/Leu) could be useful, but Clean may do too much
  #Clean
  # Load the restraints again, because 'Clean' removed them.
  LoadTbl (restrainfile),1
  ScaleRest (defaultscale)
  # Start in vacuo refinement to arrive at a suitable candidate for explicit solvent refinement
  ForceField (fofwater)
  Cutoff 10.48
  ScaleKBP 7,7,7
  ScaleForce Angle,2
  ScaleForce Planarity,8
  Sim On
  bestenergy=1e100
  bestenergyshown='not yet determined'
  failures=0
  for j=1 to cycles
    Style Ribbon
    # Set strong restraining forces
    RestrainPot SqOffset=(0.5+offset)
    ScaleRest (strongscale)
    # Heat
    ShowMessage 'Simulated annealing refinement in vacuo, cycle (j)/(cycles). Minimum violation energy is (bestenergyshown). Heating...'
    Temp (500+i*10)K
    TempCtrl Reassign
    Wait 30
    # And cool
    ShowMessage 'Simulated annealing refinement in vacuo, cycle (j)/(cycles). 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 (j): 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 candidate for explicit solvent refinement.'
      failures=0
    else
      failures=failures+1
      if failures>7
        break
    bestenergyshown=0.00+bestenergy