# YASARA MACRO
# TOPIC:       3. Molecular Dynamics
# TITLE:       Refining a homology model by molecular dynamics simulation in water
# REQUIRES:    Dynamics
# AUTHOR:      Elmar Krieger
# LICENSE:     GPL
# DESCRIPTION: This macro runs a 500 ps simulation of a homology model using the protocol described in Proteins 57,678-683. It saves one PDB file every 25 picoseconds, including a table with force field energies to identify the best snapshot. In YASARA Structure, this table additionally contains Dihedrals, Packing1D and Packing3D checks (see the Check command). And in the Twinset, three more checks are included: PhiPsi, Backbone and Packing3.

# The homology model to refine must be present with a .pdb 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'

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

# Homology model refinement parameters from Proteins 57,678-683
# See md_run.mcr for an explanation of these parameters.
# This example already uses the new YAMBER3 force field in YASARA Dynamics,
# and the knowledge-based YASARA force field in YASARA Structure.
temperature='298K'
density=0.997
pH=7.0
if Structure
  ForceField YASARA2,SetPar=Yes
else
  ForceField YAMBER3,SetPar=Yes
# Save snapshots every 2*1*12500 = 25 ps */
timestep='2,1.0'
savesteps=12500
# Create 20 snapshots, 20*25 = 500 ps simulation 
snapshots=20

# Keep the solute from diffusing around and crossing periodic boundaries
CorrectDrift On

# Treat all simulation warnings as errors that stop the macro
WarnIsError On

Clear
# Do we already have a scene with water?
size = FileSize (MacroTarget)_water.sce
if not size
  # No scene present, assume it's a PDB or YOB file
  for type in 'pdb','yob'
    size = FileSize (MacroTarget).(type)
    if size
      break
  if not size
    RaiseError "Initial structure not found. Make sure to create a project directory and place the structure there"
  obj = Load(type) (MacroTarget)
  # Align model with major axes to minimize cell size
  NiceOriObj (obj)
  Clean
  # Create the simulation cell, 2*7.5 A larger than the model along each axis
  Cell Auto,Extension=7.5
  # Fill with water, predict pKas, place counter ions
  Experiment Neutralization
    WaterDensity (density)
    pH (pH)
    pKaFile (MacroTarget).pka
    Speed Fast
  Experiment On
  Wait ExpEnd
  # Save scene with water
  SaveSce (MacroTarget)_water
else
  LoadSce (MacroTarget)_water
  
# Temperature
Temp (temperature)
# Start the Simulation
Console Off
TimeStep (timestep)
# Alread a snapshot/trajectory present?
i=00
filename='(MacroTarget)_snapshot(i).sim'
running = FileSize (filename)
if not running
  # Simulation has not been running before, start now
  # Perform energy minimization
  Experiment Minimization
  Experiment On
  Wait ExpEnd
  # And now start the real simulation
  Temp (temperature)
  Sim On
else
  # Simulation has been running before
  # Find and load the last 'sim' snapshot
  do
    i=i+1
    found = FileSize (MacroTarget)_snapshot(i).sim
  while found
  LoadSim (MacroTarget)_snapshot(i-1)

# Set temperature and pressure control
TempCtrl Rescale
PressureCtrl SolventProbe,HOH,(density)

# Wait for 500ps, saving 20 snapshots
SaveSim (MacroTarget)_snapshot00,(savesteps)
SimSteps 2
do
  t=Time
  ps=snapshots*25
  ShowMessage 'Running refinement simulation, waiting for (ps) ps, (0.00+t*0.1/ps)% completed'
  Wait 10
while t<ps*1000
  
SaveSim Off
PressureCtrl Off   

# YASARA Structure checks
strchecklist()='  Dihedrals','  Packing1D','  Packing3D'
strchecks=count strchecklist

# WHATIF checks in the Twinset
wifchecklist()='     PhiPsi','   Backbone','   Packing1'
wifchecks=count wifchecklist

# Load the snapshots, energy minimize each snapshot and tabulate energies etc.
for i=00 to snapshots-1
  ShowMessage 'Analyzing snapshot (1+i)/(snapshots)'
  filename='(MacroTarget)_snapshot(i)'
  LoadSim (filename)
  TempCtrl Anneal
  Wait 500
  SavePDB 1,(filename)
  Tabulate (i)
  Tabulate EnergyObj 1
  if Structure
    for j=1 to strchecks
      resultlist1(j) = CheckObj 1,(strchecklist(j))
      Tabulate (resultlist1(j))
  if Twinset
    for j=1 to wifchecks
      resultlist2(j) = CheckObj 1,(wifchecklist(j))
      Tabulate (resultlist2(j))
  if Structure or Twinset
    Tabulate (mean resultlist)
Sim Off

# Save the result table
columns=2
header="  Snapshot     Energy"
if Structure
  columns=columns+strchecks
  header=header+fuse strchecklist
if Twinset
  columns=columns+wifchecks
  header=header+fuse wifchecklist
if Structure or Twinset
  columns=columns+1
  header=header+"    Average"
filename='(MacroTarget)_results.tab'
SaveTab default,(filename),Format=Text,Columns=(columns),NumFormat=10.2f,(header)
ShowMessage 'The analysis results have been saved as (filename).'