# 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'

# List of objects for RMSD calculations
# If one of the objects is an oligomer, check the documentation of the 'Sup' command at 'analyzing a simulation' to avoid pitfalls.
objlist() = ListObj Protein

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

# 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

ShowMessage "Preparing analysis, please wait..."
Wait 1

# 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)
  Sim Pause
  simtime = Time
  ShowMessage 'Analyzing snapshot (0+i) at (0+(simtime/1000)) ps'
  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))
    rmsd1(j) = SupAtom CA Obj (objlist(j)),CA Obj (initiallist(j))
    rmsd2(j) = SupAtom Backbone Obj (objlist(j)),Backbone Obj (initiallist(j))
    rmsd3(j) = SupAtom Element !H Obj (objlist(j)),Element !H Obj (initiallist(j)),Flip=Yes
    Tabulate (rmsd1(j)),(rmsd2(j)),(rmsd3(j))
  # Add the average RMSDs
  Tabulate (mean rmsd1),(mean rmsd2),(mean rmsd3)
  # 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[(EnergyUnit)]_____Bond _______Angle ____Dihedral ___Planarity _____Coulomb _________VdW '
for i=0 to supobjects
  if i<supobjects
    id=000+objlist(i+1)
  else
    id='AVR'
  header=header+'RMSDs-(id):CA ____Backbone __HeavyAtoms '
SaveTab default,(MacroTarget)_analysis,Format=Text,Columns=(8+(supobjects+1)*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