# 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. If you have the Twinset installed, this table also contains PhiPsi, Backbone and Packing3 checks.

# 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.
temperature = '298K'
density = 0.997
pH = 7.0
ForceField Yamber3
Cutoff 7.86
Boundary periodic
Longrange Coulomb
Temp (temperature)
# Save snapshots every 2*1*12500 = 25 ps */
timestep = '2,1.0'
savesteps = 12500
# Create 20 snapshots, 20*25 = 500 ps simulation 
snapshots = 20

# 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
  
# Start the Simulation
Console Off
TimeStep (timestep)
Sim On
# 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
  # Do a short steepest descent energy minimization
  TempCtrl SteepDes
  ShowMessage "Starting refinement with a steepest descent minimization..."
  # Wait for convergence (until the maximum atom speed drops below 2200 m/s)
  Wait SpeedMax<2200
  # Continue with 500 simulated annealing steps
  ShowMessage "Continuing with 500 simulated annealing steps..."
  Temp 0
  TempCtrl Anneal
  Wait 500
  HideMessage
  # And now start the real simulation
  Temp (temperature)
  # Reset time to 0
  Time 0
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 'Waiting for (ps) ps, (0.000+t*0.1/ps)% completed'
  Wait 10
while t<ps*1000
  
SaveSim Off   

# Load the snapshots, energy minimize each snapshot and tabulate energies etc.
for i=00 to snapshots-1
  ShowMessage 'Analyzing snapshot (i+1)/(snapshots)'
  filename='(MacroTarget)_snapshot(i)'
  LoadSim (filename)
  TempCtrl Anneal
  Wait 500
  SavePDB 1,(filename)
  Tabulate (i)
  Tabulate EnergyObj 1
  if Twinset
    for check in 'PhiPsi','Backbone','Packing1'
      Tabulate CheckObj 1,(check)
Sim Off

# Save the result table
filename='(MacroTarget)_results.tab'
if Twinset
  SaveTab default,(filename),Format=Text,Columns=5,NumFormat=10.2f,'  Snapshot     Energy     PhiPsi   Backbone   Packing1'
else
  SaveTab default,(filename),Format=Text,Columns=2,NumFormat=10.2f,'  Snapshot     Energy'

ShowMessage 'The analysis results have been saved as (filename).'