# YASARA MACRO
# TOPIC:       3. Molecular Dynamics
# TITLE:       Analyzing the change of secondary structure in a molecular dynamics trajectory
# REQUIRES:    Dynamics
# AUTHOR:      Elmar Krieger
# LICENSE:     GPL
# DESCRIPTION: This macro analyzes a simulation and creates a table with per-residue secondary structure as well as a graphical plot

# Selection of residues whose secondary structure will be analyzed,
# here we take all protein residues in object 1
selection='Protein Obj 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
# 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
HideRes HOH
Move X=-50,Z=200

# Rearrange the objects so that there are no empty spots in the object list
RenumberObj All

Console Off

# Start the text table
Tabulate 0

# Get a list of CA atoms to analyze secondary structure
calist() = ListAtom CA Res (selection)
residues = count calist

# Create a text object with the residue names
lineheight = 2.5
textwidth = 4000
textobj = MakeTextObj Residues,(textwidth),(lineheight*residues)
Font Arial,Height=(0.6*lineheight),Color=Yellow,Depth=0.5,DepthCol=Red
for i=1 to residues
  PosText X=(textwidth*0.5-12),Y=(lineheight*(residues-i+1)),justify=left
  params = 'Atom (calist(i)),Format=MOLNAME RESName RESNUM'
  resid = ListRes (params)  
  Print (resid)
  Tabulate ListRes (params)
Wait 1

# Get number of next free object, which will contain the plot
plotobj = Objects+1

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)
  t = Time
  t = 0+t/1000
  Tabulate (t)
  if i%10==0
    # Label the time axis every 10th dot
    PosText X=(lineheight*i+0.5*textwidth),Y=0,justify=left
    Print (t)
  # Color residues by current secondary structure
  ColorRes (selection),SecStr
  if i
    # The plot object must not be part of the soup when loading the snapshot
    AddObj (plotobj)
  for j=1 to residues
    # Build one atom 'dot' for each residue
    resobj = BuildAtom Carbon
    # Transfer secondary structure color
    color = ColorAtom (calist(j))
    ColorObj (resobj),(color)
    # Add secondary structure to table
    Tabulate SecStrRes Atom (calist(j))
    # Shift the atom to the proper spot
    MoveAtom Obj (resobj),X=(lineheight*i),Y=(lineheight*(1.25+0.5*residues-j))
    # Join atoms to plot object, to gain speed we do this only once in a while
    if Objects>240 or j==residues
      JoinObj (plotobj)-(Objects),(plotobj),Center=No
  # Next snapshot
  ShowMessage 'Analyzing snapshot (i) at (t) picoseconds'
  Wait 1
  RemoveObj (plotobj)
  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"

AddObj (plotobj)
DelObj not (textobj) (plotobj)
Zoom
# Save table
SaveTab default,(MacroTarget)_SecStrAnalysis,Format=Text,Columns=(residues+1),NumFormat=%11.0f,'Secondary structure analysis'
HideMessage