# YASARA MACRO
# TOPIC:       5. Structure Determination
# TITLE:       Analyze and merge an NMR ensemble
# REQUIRES:    YASARA Structure and NMR Structure Determination Module
# AUTHOR:      Elmar Krieger, Sander Nabuurs, Chris Spronk
# LICENSE:     GPL
# DESCRIPTION: This macro analyzes an NMR ensemble and merges the members into a single PDB file

# To analyze an ensemble, you need:
# - An ensemble of refined 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,
# 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 4: Analyze the ensemble
# ============================
Console Off
Clear
LogAs (logfile),append=Yes, print 'Explicit solvent refinement:'
# We cannot simply load all files present, the user may have restarted with different 'structures',
# and the leftovers may be for different restraints. 'i' must start with 1.
for i=1 to structures-1
  # Try with and without leading zeroes, maybe start in nmr_setdefaults was changed
  pdbfilename='(ensemblefile)(start+i).pdb'
  size = FileSize (pdbfilename)
  if !size
    RaiseError 'The ensemble member `(pdbfilename)` was not found'
  # Analyze the best structure
  LoadPDB (pdbfilename)
  # Create a new large simulation cell, so that the influence of periodic
  # boundaries on the energies is smaller.
  Cell Auto,Extension=20
  # Load the restraints
  RestrainPot (defaultpot)
  RestrainPar (defaultpar)
  Load(restrainformat) (restrainfile),(i),NameFormat=(nameformat)
  ScaleRest Distance=(defaultscale1),Dihedral=(defaultscale2),RDC=(defaultscale3)
  # Select force field for explicit solvent
  ForceField (fofsolvent)
  Cutoff 8
  Boundary periodic
  Longrange Coulomb
  # Calculate energies
  ShowMessage 'Calculating energies for ensemble structure (i)...'
  Wait 1
  violenergylist(i) = RestEnergy
  disenergylist(i) = RestEnergy Distance
  dihenergylist(i) = RestEnergy Dihedral
  rdcenergylist(i) = RestEnergy RDC
  dissum,dismax,dihsum,dihmax,rdcsum,rdcmax = RestViol
  soluteenergylist(i) = EnergyAll
  solvcoulombenergylist(i),solvvdwenergylist(i) = SolvEnergy
  fofenergylist(i)=soluteenergylist(i)+solvcoulombenergylist(i)+solvvdwenergylist(i)
  totenergylist(i)=fofenergylist(i)+violenergylist(i)
  # Calculate model quality Z-score
  ForceField YASARA2
  qualitylist(i) = CheckAll ModelQuality
  # Create a list of energies and check results, which will be sorted below
  sortlist(i)=fofenergylist(i)*fofscale+violenergylist(i)
  checklist(i)='    Total energy                   (0.000+totenergylist(i)) (EnergyUnit)\n'+
               '      Restraint violation energy:  (0.000+violenergylist(i)) (EnergyUnit)\n'+
               '        Distance restraints:       (0.000+disenergylist(i)) (EnergyUnit)\n'+
               '          Maximum violation:       (0.000+dismax) A\n'+
               '        Dihedral angle restraints: (0.000+dihenergylist(i)) (EnergyUnit)\n'+
               '          Maximum violation:       (0.000+dihmax) degrees\n'+
               '        RDC restraints:            (0.000+rdcenergylist(i)) (EnergyUnit)\n'+
               '          Maximum violation:       (0.000+rdcmax) Hz\n'+
               '      Force field energy:          (0.000+fofenergylist(i)) (EnergyUnit)\n'+
               '        Internal solute energy:    (0.000+soluteenergylist(i)) (EnergyUnit)\n'+
               '        Electrostatic solv.energy: (0.000+solvcoulombenergylist(i)) (EnergyUnit)\n'+
               '        Van der Waals solv.energy: (0.000+solvvdwenergylist(i)) (EnergyUnit)\n'+
               '    Model Quality Z-score:         (0.000+qualitylist(i))'
  if Twinset
    # Check the structure with WHAT IF
    ramchklist(i)=CheckObj (i),PhiPsi
    bbcchklist(i)=CheckObj (i),Backbone
    quachklist(i)=CheckObj (i),Packing1
    checklist(i)=checklist(i)+'\n'+
                 '    Ramachandran plot Z-score:     (0.000+ramchklist(i))\n'+
                 '    Backbone Z-score:              (0.000+bbcchklist(i))\n'+
                 '    Packing quality Z-score:       (0.000+quachklist(i))'
  RemoveObj (i)
  DelObj SimCell

# Sort the ensemble, lowest energies first
AddObj all
do
  sorted=1
  for i=1 to Objects-1
    if sortlist(i)>sortlist(i+1)
      SwapObj (i),(i+1)
      swap=sortlist(i)
      sortlist(i)=sortlist(i+1)
      sortlist(i+1)=swap
      swap=checklist(i)
      checklist(i)=checklist(i+1)
      checklist(i+1)=swap
      sorted=0
while not sorted

# Superpose the ensemble using Theseus
rmsdlist() = SupMultiAtom Backbone,Method=Theseus

# Create the final logfile output
Console Off
RecordLog (logfile),append=yes
print 'Ensemble analysis:'
for i=1 to Objects
  print '  Ensemble structure (i):\n(checklist(i))'
  print '    Backbone RMSD from average:    (0.000+rmsdlist(i)) A'
  # Load the restraints again to list the violations
  print '    List of violated restraints:'
  Load(restrainformat) (restrainfile),(i),NameFormat=(nameformat)
  ScaleRest Distance=(defaultscale1),Dihedral=(defaultscale2),RDC=(defaultscale3)
  ListRestAll Class=Violate,Sort=Yes,Spaces=6
  print ' '
  DelRestAll
  # Overwrite the ensemble members in the new sort order
  SavePDB (i),(ensemblefile)(000+i)
print '\nEnsemble average results for (Objects) members:'
print '  Total energy                   (0.000+mean totenergylist) (EnergyUnit)'
print '    Restraint violation energy:  (0.000+mean violenergylist) (EnergyUnit)'
print '      Distance restraints:       (0.000+mean disenergylist) (EnergyUnit)'
print '      Dihedral angle restraints: (0.000+mean dihenergylist) (EnergyUnit)'
print '      RDC restraints:            (0.000+mean rdcenergylist) (EnergyUnit)'
print '    Force field energy:          (0.000+mean fofenergylist) (EnergyUnit)'
print '      Internal solute energy:    (0.000+mean soluteenergylist) (EnergyUnit)'
print '      Electrostatic solv.energy: (0.000+mean solvcoulombenergylist) (EnergyUnit)'
print '      Van der Waals solv.energy: (0.000+mean solvvdwenergylist) (EnergyUnit)'
print '  Model Quality Z-score:         (0.000+mean qualitylist)'
if Twinset
  print '  Ramachandran plot Z-score:     (0.000+mean ramchklist)'
  print '  Backbone Z-score:              (0.000+mean bbcchklist)'
  print '  Packing quality Z-score:       (0.000+mean quachklist)'
print '  Backbone RMSD from average:    (0.000+mean rmsdlist) A'
StopLog
HideMessage
Console On

# Save the ensemble as one PDB file
SavePDB all,(ensemblefile)