# YASARA MACRO
# TOPIC:       3. Molecular Dynamics
# TITLE:       Analyzing a molecular dynamics trajectory of multiple objects
# REQUIRES:    Dynamics
# AUTHOR:      Elmar Krieger
# LICENSE:     GPL
# DESCRIPTION: This macro analyzes a simulation and creates a table with energies and RMSDs from the starting structures for multiple objects

# The structure to analyze must be present with a .sce extension.
# You can either set the target structure by clicking on Options > Macro > Set target,
# or by uncommenting the line below and specifying it directly.
#MacroTarget = 'c:\MyProject\1crn'

# Total number of objects for RMSD calculations, ObjList contains the individual object numbers
# This will calculate RMSDs for 2 objects: object 1 and object 3
# If one of the objects is an oligomer, check the documentation of the 'Sup' command at 'analyzing a simulation' to avoid pitfalls.
ObjList = 1,3

# Forcefield to use (these are all YASARA commands, so no '=' used)
ForceField Amber99

# Cutoff
Cutoff 7.86

# Cell boundary
Boundary periodic

# Use longrange coulomb forces (particle-mesh Ewald)
Longrange Coulomb

# Do we have a target?
if MacroTarget==''
  RaiseError "This macro requires a target. Either edit the macro file or click Options > Macro > Set target to choose a target structure"

Clear
# Do we have a scene with water?
scene = FileSize (MacroTarget)_water.sce
if not scene
  RaiseError 'Could not find initial scene file (MacroTarget)_water.sce'

# Load the scene
LoadSce (MacroTarget)_water
# Duplicate the objects
supobjects = count ObjList
for i=1 to supobjects
  InitialObj(i) = DuplicateObj (ObjList(i))
  RemoveObj (InitialObj(i))

i=00000
while 1
  # See if next snapshot is present
  sim = FileSize (MacroTarget)(i).sim
  if not sim
    break
  # Yes, load it
  LoadSim (MacroTarget)(i)
  ShowMessage 'Analyzing snapshot (0+i)'
  # Add time in picoseconds to table
  simtime = Time
  Tabulate (simtime/1000)
  # Add energy to table, first the total energy, then all components individually
  Tabulate Energy
  Tabulate Energy All
  # Add CA, backbone and heavy atom RMSDs to table
  for j=1 to supobjects
    AddObj (InitialObj(j))
    rmsd=SupAtom CA Obj (InitialObj(j)),CA Obj (ObjList(j))
    Tabulate (rmsd)
    rmsd=SupAtom Backbone Obj (InitialObj(j)),Backbone Obj (ObjList(j))
    Tabulate (rmsd)
    rmsd=SupAtom Element !H Obj (InitialObj(j)),Element !H Obj (ObjList(j)),Flip=Yes
    Tabulate (rmsd)
    RemoveObj (InitialObj(j))
  # Next snapshot
  i=i+1

if !i
  RaiseError "This macro is meant to analyze a molecular dynamics trajectory created with md_run, but none was found in this directory"

# Save table
header='____Time[ps] Energy:Total ________Bond _______Angle ____Dihedral ___Planarity _____Coulomb _________VdW '
for i=1 to supobjects
  header=header+'RMSDs-(000+ObjList(i)):CA ____Backbone __HeavyAtoms '
SaveTab default,(MacroTarget)_Analysis,Format=Text,Columns=(8+supobjects*3),NumFormat=%12.3f,(header)