# 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.
currobj = 1

# 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

# Consider protein residues with a Calpha plus nucleic acid residues with a C1*
atomsel='Atom CA Protein or Atom C1* NucAcid'

# First snapshot to be analyzed, increase number to ignore an equilibration period.
# (By default, md_run.mcr saves snapshots every 25ps, choosing 40 thus starts the analysis after 1 nanosecond)
firstsnapshot=0

# No change required below this point
# ===================================

# 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

# We need at least 3 Calpha atoms to superpose
calphas = CountAtom (atomsel)
if calphas<3
  RaiseError 'This macro currently requires at least three residues to analyze, because it performs a superposition and RMSD calculation'

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

# Duplicate the intial object for RMSD calculation
startobj = DuplicateObj (currobj)

# Get a selection of residues that match atomsel in the current and the start objects
for type in 'curr','start'
  reslist() = ListRes (atomsel) and Obj ((type)obj)
  (type)ressel=join reslist
RemoveObj (startobj)

# Set the summed up RMSDs for CA, BackBone and HeavyAtoms to zero
rmsds=count reslist
for i=1 to 3
  for j=1 to rmsds
    rmsd(i)sum(j)=0.

i=00000+firstsnapshot
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
  # Get per residue RMSDs between current residue selection and start residue selection
  AddObj (startobj)
  rmsd1() = SupAtom (atomsel) and Obj (currobj), (atomsel) and Obj (startobj),Unit=Res
  rmsd2() = SupAtom Res (currressel) Backbone, Res (startressel) Backbone,Unit=Res
  rmsd3() = SupAtom Res (currressel) Element !H, Res (startressel) Element !H,Unit=Res
  # Sum up atom positions to calculate RMSFs
  AddPosAtom Res (currressel)
  # 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 (startobj)
  # 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 (atomsel) and Obj (currobj)
reslist() = ListRes (atomsel) and Obj (currobj),Format='RESNAME RESNUM MOLNAME'
for j=1 to rmsds
  Tabulate '(reslist(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