# YASARA MACRO
# TOPIC:       3. Molecular Dynamics
# TITLE:       Running an accurate molecular dynamics simulation in water
# REQUIRES:    Dynamics
# AUTHOR:      Elmar Krieger
# LICENSE:     GPL
# DESCRIPTION: This macro sets up and runs a simulation. It can also continue a simulation that got interrupted.

# Parameter section - adjust as needed
# ====================================

# The structure to simulate must be present with a .pdb or .sce extension.
# If a .sce (=YASARA scene) file is present, the cell must have been added.
# 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'

# Simulation temperature
temperature = '298K'

# Name of solvent molecules used to measure solvent density
solvent = 'HOH'

# Solvent density in g/ml (0.997 is water at 298K)
# If you run at a significantly different temperature, this density needs to
# be adapted. Look at the PressureCtrl documentation.
# For solvents other than water, you have to create your own solvent box
# as described at the FillCellObj documentation and save it as .._solvent.sce
density = 0.997

# pH at which the simulation should be run
pH = 7.0

# NaCl concentration in mass percent (0.9% is a physiological solution)
NaCl = 0.9

# Multiple timestep: 1.25 femtoseconds for intramolecular and 2*1.25 fs for intermolecular forces
timestep = '2,1.25'

# Save simulation snapshots every 3000 simulation steps
# (with a timestep of 2*1.25=2.5 femtoseconds, that's 3000*2.5 fs = 7.5 picoseconds).
savesteps = 3000

# The format used to save the trajectories, sim or xtc
format = 'sim'

# 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

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

# Normally no change required below this point
# Temperature
Temp (temperature)

# 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 already have a 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
  # No scene with water present yet
  if solvent!='HOH'
    RaiseError "When using a solvent other than water, the cell cannot be filled automatically. Look at the documentation of the FillCellObj command"
  # Do we have a scene at all?
  scene = FileSize (MacroTarget).sce
  if scene
    LoadSce (MacroTarget)
  else
    # No scene present, assume it's a PDB or YOB file
    for type in 'pdb','yob','Error'
      size = FileSize (MacroTarget).(type)
      if size
        obj = Load(type) (MacroTarget)
        # Align object with major axes to minimize cell size
        NiceOriObj (obj)
        break
    if type=='Error'
      RaiseError "Initial structure not found. Make sure to create a project directory and place the structure there"
    Clean
    # Create the simulation cell, 2*10 A larger than the protein along each axis
    Cell Auto,Extension=10
    SaveSce (MacroTarget)
  # Fill with water, predict pKas, place counter ions
  Experiment Neutralization
    WaterDensity (density)
    pH (pH)
    NaCl (NaCl)
    pKaFile (MacroTarget).pka
    Speed Fast
  Experiment On
  Wait ExpEnd
  # Save scene with water
  SaveSce (MacroTarget)_water

# Start the Simulation
TimeStep (timestep)
Sim On
# Make sure all atoms are free to move
Free
# Alread a snapshot/trajectory present?
i=00000
filename = '(MacroTarget)(i).(format)'
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 simulation 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
  # And now start with the real simulation
  HideMessage
  Temp (temperature)
  # Reset time to 0
  Time 0
else
  # Simulation has been running before
  # Switch console off to load the snapshots quickly
  Console Off
  if format=='sim'
    # Find and load the last 'sim' snapshot
    do
      i = i+1
      found = FileSize (MacroTarget)(i).sim
    while found
    LoadSim (MacroTarget)(i-1)
  else
    # XTC format requires that the entire trajectory is read in to find the last one
    do
      i = i+1
      eof,time = LoadXTC (filename),(i)
      ShowMessage 'Searching XTC trajectory for last snapshot, showing snapshot (i) at (0+time) fs'
      Wait 1
    while !eof
    HideMessage
  Console On

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

# And finally, make sure that future snapshots are saved
Save(format) (filename),(savesteps)

Wait forever