# YASARA MACRO
# TOPIC:       5. Structure Determination
# TITLE:       Fold an ensemble of NMR structures
# REQUIRES:    YASARA Structure and NMR Structure Determination Module
# AUTHOR:      Elmar Krieger, Sander Nabuurs, Chris Spronk
# LICENSE:     GPL
# DESCRIPTION: This macro folds a structure from the stretched-out conformation using NMR restraints

# To fold a structure, you need:
# - The protein sequence, in FASTA or PDB format,
#   default filenames are 'sequence.pdb' or 'sequence.fasta'
# - 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,
# by providing it as command line argument (see docs at Essentials > The command line),
# 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 1: Fold the stretched out conformation to arrive at an approximate fold
# ============================================================================
# Especially for larger proteins, the result is far from accurate, but well
# suited to arrive at the real structure during a molecular dynamics simulation.

Clear
Console off
# Do we have a starting structure?
for type in 'yob','pdb'
  size = FileSize (MacroTarget)_start.(type)
  if size
    break
if size
  # Yes, load it
  Load(type) (MacroTarget)_start
  LogAs (logfile),append=no,print 'Folding (structures*foldmul) structures using starting structure (MacroTarget)_start.(type) and restraints `(restrainfile)`, keeping best (structures)' 
else
  # No, build a linear peptide from the sequence
  BuildMol (seqfile)
  NumberRes Obj 1,First=(firstresnum)
  LogAs (logfile),append=no,print 'Folding (structures*foldmul) structures using sequence (seqfile) and restraints `(restrainfile)`, keeping best (structures)' 
CleanAll
# Build cysteine bridges if specified by the user, these will be treated as distance restraints
if cysbridgelist(1)!='Auto' and cysbridgelist(1)!='None'
  for i=1 to count cysbridgelist step 2
    DelAtom HG Res (cysbridgelist(i)) (cysbridgelist(i+1))
    AddBond SG Res (cysbridgelist(i)),SG Res (cysbridgelist(i+1))
Style Tube
# Set the BFactors for floating assignment
BFactorAll 0
BFactorAtom (floating),25
# We need a forcefield with knowledge based potentials for sampling
ForceField YASARA
# Set restraining potential and parameters
RestrainPot (defaultpot)
RestrainPar (defaultpar)
if correctcispro=='Yes' or correctcispro=='yes'
  # Cis-prolines are not allowed. We therefore exclude the peptide bonds
  # before prolines from sampling by fixing the C and N atom:
  FixAtom N Res Pro
  FixAtom C with bond to N Res Pro
# Fold the structures from the stretched out conformation
Experiment NMRFolding
  # The object to fold
  StartObj 1
  # The file(s) with NMR restraints (extension is .tbl or .nef)
  RestrainFile (restrainfile),NameFormat=(nameformat)
  # The number of structures to generate
  Structures (structures*foldmul)
  # The filename for the structures (digits will be incremented)
  StructureFile (ensemblefile)(start)
  # The random seed, simply the number of the starting structure.
  RandomSeed (start)
# Generate the initial structures
Experiment On
Wait ExpEnd

# Keep only the best 'structures'
ForceField (fofvacuo),SetPar=Yes
end=start+structures*foldmul
for i=start to end-1
  Clear
  # Load ensemble member and apply restraints
  LoadPDB (ensemblefile)(i)
  Load(restrainformat) (restrainfile),1,NameFormat=(nameformat)
  ScaleRest Distance=(defaultscale1),Dihedral=(defaultscale2),RDC=(defaultscale3)
  ShowMessage 'Minimizing newly folded structure (i)/(end-1) to calculate energies...'
  # Add cell
  Cell Auto,Extension=10
  # Wall boundaries speed up restraint force calculations
  Boundary Wall
  # Remove worst bumps with steepest descent minimization
  TempCtrl SteepDes
  Sim On
  SimSteps 10
  for j=1 to 100
    Wait 1
    if SpeedMax<4000
      break
  # Annealing, first with small time-step to avoid simulation failures during bump removal
  Time 0
  TempCtrl Anneal
  TimeStep 2,0.5
  Wait 50
  TimeStep 2,1
  Wait 50
  # Now the most bumps are gone, the remaining ones are most likely caused by knots, remove them
  CorrectKnots on
  Wait 50
  CorrectKnots off    
  Wait ClosedBonds
  # Sum up restraint violation and force field energies (the latter scaled with 'fofscale')
  violenergy = RestEnergy
  fofenergy = Energy
  sortenergy=fofenergy*fofscale+violenergy
  LogAs (logfile),append=yes,print '  Structure (i) has energy (0.000+sortenergy) (EnergyUnit), including (0.000+violenergy) (EnergyUnit) restraint violations'
  Sim off
  if count sortenergylist<structures
    # Simply accept the initial 'structures'
    sortenergylist(count sortenergylist+1)=sortenergy
    SavePDB 1,(ensemblefile)(i)
  else
    # Get the worst structure
    for j=1 to structures
      if j==1 or sortenergylist(j)>worstenergy
        worstidx=j
        worstenergy=sortenergylist(j)
    if sortenergy<worstenergy    
      # Keep current structure, throwing out the worst
      sortenergylist(worstidx)=sortenergy
      SavePDB 1,(ensemblefile)(start+worstidx-1)
      LogAs (logfile),append=yes,print '    Replacing structure (start+worstidx-1)'
    DelFile (ensemblefile)(i).pdb