# YASARA MACRO
# TOPIC:       3. Molecular Dynamics
# TITLE:       Analyzing a molecular dynamics trajectory per residue
# REQUIRES:    Dynamics
# AUTHOR:      Elmar Krieger
# LICENSE:     GPL
# DESCRIPTION: This macro analyzes a simulation and creates a table with average RMSFs and RMSDs from the starting structure per residue

# 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 RMSFs and RMSDs from the starting conformation will be calculated
# If this 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

# 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)
# Set the summed up RMSDs for CA, BackBone and HeavyAtoms to zero
rmsds = CountRes Protein Obj (current)
for i=1 to 3
  for j=1 to rmsds
    rmsd(i)sum(j)=0.

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)
  ShowMessage 'Analyzing snapshot (0+i)'
  Wait 1
  # Get per residue RMSDs
  AddObj (initial)
  rmsd(1)()=SupAtom CA Protein Obj (current), CA Protein Obj (initial),Unit=Res
  rmsd(2)()=SupAtom Protein Backbone Obj (current), Protein Backbone Obj (initial),Unit=Res
  rmsd(3)()=SupAtom Protein Element !H Obj (current), Protein Element !H Obj (initial),Unit=Res
  # Sum up atom positions to calculate RMSFs
  AddPosAtom Protein Obj (current)
  # Sum up RMSDs
  for j=1 to 3
    for k=1 to rmsds
      rmsd(j)sum(k)=rmsd(j)sum(k)+rmsd(j)(k)
  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"

# Calculate the average
rmsflist()=RMSFRes Protein Obj (current)
for j=1 to rmsds
  Tabulate (j),(rmsd(1)sum(j)/i),(rmsd(2)sum(j)/i),(rmsd(3)sum(j)/i),(rmsflist(j))
# Save table
SaveTab default,(MacroTarget)_analysisres,Format=Text,Columns=5,NumFormat=%12.3f,_____Residue _RMSDs[A]:CA ____Backbone __HeavyAtoms______RMSF[A]
HideMessage