# YASARA MACRO
# TOPIC:       3. Molecular Dynamics
# TITLE:       Analyzing a molecular dynamics trajectory
# REQUIRES:    Dynamics
# AUTHOR:      Elmar Krieger
# LICENSE:     GPL
# DESCRIPTION: This macro analyzes a simulation and creates a table with energies and RMSDs from the starting structure, as well as the time average structure with B-factors calculated from the root mean square fluctuations.

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

# Number of the object whose RMSDs from the starting conformation will be calculated
# If the protein is an oligomer, check the documentation of the 'Sup' command at 'analyzing a simulation' to avoid pitfalls.
current = 1

# 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 intial object for RMSD calculation
initial = DuplicateObj (current)
RemoveObj (initial)

Sim Pause
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)
  # Add time in picoseconds to table
  simtime = Time
  ShowMessage 'Analyzing snapshot (0+i) at (0+simtime) fs'
  Wait 1
  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
  AddObj (initial)
  rmsd=SupAtom CA Obj (current),CA Obj (initial)
  Tabulate (rmsd)
  rmsd=SupAtom Backbone Obj (current),Backbone Obj (initial)
  Tabulate (rmsd)
  rmsd=SupAtom Element !H Obj (current),Element !H Obj (initial)
  Tabulate (rmsd)
  # Add the current atom positions to internal table to obtain RMSF and average positions
  AddPosAtom Obj (current)
  RemoveObj (initial)
  # 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 main analysis table
eunit=EnergyUnit
if eunit=='kJ/mol'
  eunit=' kJ/mol '
SaveTab default,(MacroTarget)_analysis,Format=Text,Columns=11,NumFormat=12.3f,____Time[ps] Energy[(eunit)]_____Bond _______Angle ____Dihedral ___Planarity _____Coulomb _________VdW _RMSDs[A]:CA ____Backbone __HeavyAtoms

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

"""
# Example: Create an additional RMSF file, in case B-factors are too large for the PDB format
DelTab default
first,last = SpanAtom Obj (current)
rmsflist() = RMSFAtom Obj (current)
for i=first to last
  res = ListAtom (i),Format='RESNAME,RESNUM'
  Tabulate (i),(res),(rmsflist(i-first+1))
SaveTab default,(MacroTarget)_rmsf,Format=Text,Columns=4,NumFormat=8.2f,"RMSF Table"  
"""