# 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,2,3,4

# 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

# The B-factors calculated from the root-mean-square fluctuations can be too large to fit them
# into the PDB file's B-factor column. Replace e.g. 1.0 with 0.1 to scale them down to 10%
bfactorscale=1.0

# 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
Console Off
# 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
  initiallist(i) = DuplicateObj (objlist(i))
  RemoveObj (initiallist(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)'
  Wait 1
  # 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 (initiallist(j))
    Tabulate SupAtom CA Obj (objlist(j)),CA Obj (initiallist(j))
    Tabulate SupAtom Backbone Obj (objlist(j)),Backbone Obj (initiallist(j))
    Tabulate SupAtom Element !H Obj (objlist(j)),Element !H Obj (initiallist(j)),Flip=Yes
  # Add the current atom positions to internal table to obtain RMSF and average positions
  AddPosAtom Obj (join objlist)
  RemoveObj (join initiallist)  
  # 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)

# Calculate time-average structure and set B-factors according to RMSF
AveragePosAtom Obj (join objlist)
RMSFAtom Obj (join objlist),Unit=BFactor
for obj in objlist
  if bfactorscale!=1.0
    # Scale B-factors so that they fit into the PDB format
    first,last=SpanAtom Obj (obj)
    for i=first to last
      bf=BFactorAtom (i)
      BFactorAtom (i),(bf*bfactorscale)    
  # Join to single object
  JoinObj (obj),(objlist(1))
# The time average structure has incorrect covalent geometry and should be energy minimized
SavePDB (objlist(1)),(MacroTarget)_average
HideMessage