# 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 (save a ray-traced screenshot with 'Ball zoom' set to 2).

# 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,
# by providing it as command line argument (see docs at Essentials > The command line),
# or by uncommenting the line below and specifying it directly.
#MacroTarget 'c:\MyProject\1crn'

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

# Selection of residues whose secondary structure will be analyzed,
# here we take all protein residues in object 1
selection='Protein Obj 1'

# Set to 0 if you do not want to create a graphical plot (makes the macro run faster)
plot=1

# Maximum number of secondary structure columns in the plot
columnsmax=100

# 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

# Load scene with water or other solvent
waterscene = FileSize (MacroTarget)_water.sce
solventscene = FileSize (MacroTarget)_solvent.sce
if waterscene
  LoadSce (MacroTarget)_water
elif solventscene 
  LoadSce (MacroTarget)_solvent
else
  RaiseError 'Could not find initial scene file (MacroTarget)_water.sce. You must run a simulation with the macro md_run first'

HideRes HOH
MoveAll X=-50,Z=200
Longrange None

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

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

# Start the text table
Tabulate 'Time'

# 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 and the table header
lineheight = 2.5
textwidth = 4000
if plot
  textobj = MakeTextObj Residues,(textwidth),(lineheight*(residues+1))
  Font Arial,Height=(0.6*lineheight),Color=Yellow,Depth=0.5,DepthCol=Red
  # Get number of next free object, which will contain one column of the plot
  colobj = Objects+1
  # Build a stretch of carbon atoms, one for each residue
  for i=1 to residues
    obj = BuildAtom Carbon
    resname = NameRes Atom (calist(i))
    resnum = ListRes Atom (calist(i)),Format=RESNUM
    # Shift the atom to the proper spot
    MoveAtom Obj (obj),Y=(lineheight*(1.75+0.5*residues-i))
    NameAtom Obj (obj),CA
    NameRes Obj (obj),(resname)
    NumberRes Obj (obj),(i)
    if i==1
      colobj = obj
    else
      JoinObj (obj),(colobj),Center=No
  JoinRes Obj (colobj)
  # Orient the object away, so that one can better watch the progress
  OriObj (textobj) (colobj),Beta=-65
  RemoveObj (colobj)
# Create the list of residues on the left side
for i=1 to residues
  params = 'Atom (calist(i)),Format=MOLNAME RESName RESNUM'
  Tabulate ListRes (params)
  if plot
    PosText X=(textwidth*0.5-12),Y=(lineheight*(residues-i+2)),justify=left
    resid = ListRes (params)  
    Print (resid)
Wait 1
if plot
  # Count the number of snapshots
  i=00000+firstsnapshot
  while 1
    sim = FileSize (MacroTarget)(i).sim
    if not sim
      break
    i=i+10
  # If there are too many snapshots, don't display all of them
  step=i//columnsmax+1
columns=0

# Backwards compatibility: Starting with YASARA version 12.8.1, XTC trajectories no longer contain a number in the filename
old = FileSize (MacroTarget)00000.xtc
if old
  RenameFile (MacroTarget)00000.xtc,(MacroTarget).xtc
# Determine trajectory format
for format in 'xtc','mdcrd','sim'
  found = FileSize (MacroTarget).(format)
  if found
    break
i=00000+firstsnapshot
last=0
while !last
  # Load next snapshot from SIM, XTC or MDCRD trajectory
  if format=='sim'
    sim = FileSize (MacroTarget)(i).sim
    if not sim
      break
    LoadSim (MacroTarget)(i)
  else
    # Set last (end of file) to 1 if last snapshot loaded
    last = Load(format) (MacroTarget),(i+1)
  Sim Pause
  t = Time
  t=0+t/1000
  Tabulate (t)
  for j=1 to residues
    # Add secondary structure to table
    Tabulate SecStrRes Atom (calist(j))
  if plot 
    if !(i%step)
      x=lineheight*columns
      if columns%5==0 and x<2000
        # Label the time axis every 10th dot
        PosText X=(0.5e0*textwidth+x),Y=0,justify=center
        Print (t)
      # Color residues by current secondary structure
      ColorRes (selection),SecStr
      # Duplicate column
      AddObj (colobj)
      obj = DuplicateObj (colobj)
      # Transfer secondary structure color
      for j=1 to residues
        color = ColorAtom (calist(j))
        ColorRes (j) Obj (obj),(color)
      # Shift the atoms to the proper spot
      MoveAtom Obj (obj),X=(x)
      if !i
        plotobj=obj
      else
        AddObj (plotobj)
        JoinObj (obj),(plotobj),Center=No
      Wait 1
      RemoveObj (plotobj) (colobj)
      columns=columns+1
  else
    Wait 1
  # Next snapshot
  ShowMessage 'Analyzing snapshot (i) at (t) picoseconds'
  i=i+1

if i==firstsnapshot
  RaiseError "This macro is meant to analyze a molecular dynamics trajectory created with md_run, but none was found in this directory"

# Save table
SaveTab default,(MacroTarget)_SecStrAnalysis,Format=Text,Columns=(residues+1),NumFormat=%11.0f,'Secondary structure analysis'
if plot
  AddObj (plotobj)
  DelObj not (textobj) (plotobj)
  OriObj all,Beta=0
  ZoomAll
  CenterObj all
HideMessage

# Exit YASARA if this macro was provided as command line argument in console mode and not included from another macro
if runWithMacro and ConsoleMode and !IndentationLevel
  Exit