# 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,
# 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
for i=001 to structures
  # Analyze the best structure
  LoadPDB (ensemblefile)(i)
  # 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)
  LoadTbl (restrainfile),(i)
  ScaleRest (defaultscale)
  # Select force field
  ForceField (fofwater)
  Cutoff 7.86
  Boundary periodic
  Longrange Coulomb
  # Calculate energies
  ShowMessage 'Calculating energies for ensemble structure (i)...'
  Wait 1
  violenergy=RestEnergy
  disenergy=RestEnergy Distance
  dihenergy=RestEnergy Dihedral
  dissum,dismax,dihsum,dihmax=RestViol
  soluteenergy=Energy
  solvcoulombenergy,solvvdwenergy=SolvEnergy
  fofenergy=soluteenergy+solvcoulombenergy+solvvdwenergy
  totenergy=fofenergy+violenergy
  # Create a list of energies and check results, which will be sorted below
  energylist(i)=fofenergy*fofscale+violenergy
  checklist(i)='    Total energy                   (0.000+totenergy) (EnergyUnit)\n'+
               '      Restraint violation energy:  (0.000+violenergy) (EnergyUnit)\n'+
               '        Distance restraints:       (0.000+disenergy) (EnergyUnit)\n'+
               '          Maximum violation:       (0.000+dismax) A\n'+
               '        Dihedral angle restraints: (0.000+dihenergy) (EnergyUnit)\n'+
               '          Maximum violation:       (0.000+dihmax) degrees\n'+
               '      Force field energy:          (0.000+fofenergy) (EnergyUnit)\n'+
               '        Internal solute energy:    (0.000+soluteenergy) (EnergyUnit)\n'+
               '        Electrostatic solv.energy: (0.000+solvcoulombenergy) (EnergyUnit)\n'+
               '        Van der Waals solv.energy: (0.000+solvvdwenergy) (EnergyUnit)'
  if Twinset
    # Check the structure with WHAT IF
    ramchk=CheckObj (i),PhiPsi
    bbcchk=CheckObj (i),Backbone
    quachk=CheckObj (i),Packing1
    checklist(i)=checklist(i)+'\n'+
                 '    Ramachandran plot Z-score:     (0.000+ramchk)\n'+
                 '    Backbone Z-score:              (0.000+bbcchk)\n'+
                 '    Packing quality Z-score:       (0.000+quachk)'
  RemoveObj (i)
  DelObj SimCell

# Sort the ensemble, lowest energies first
AddObj all
do
  sorted=1
  for i=001 to structures-1
    if energylist(i)>energylist(i+1)
      SwapObj (i),(i+1)
      swap=energylist(i)
      energylist(i)=energylist(i+1)
      energylist(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
SupMultiAtom Backbone,Method=Theseus

# Create the final logfile output
Console Off
RecordLog (logfile),append=yes
print 'Ensemble analysis:'
for i=001 to structures
  print '  Ensemble structure (i):\n(checklist(i))'
  # Overwrite the ensemble members in the new sort order
  SavePDB (i),(ensemblefile)(i),Format=XPLOR
StopLog
Console On

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