# YASARA MACRO
# TOPIC:       3. Molecular Dynamics
# TITLE:       Running a steered molecular dynamics simulation in water
# REQUIRES:    Dynamics
# AUTHOR:      Elmar Krieger
# LICENSE:     GPL
# DESCRIPTION: This macro sets up and runs a steered simulation to pull a ligand from its binding site.

# 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\1bfb'

# Selection of ligand residue(s), must be adapted
# (if the ligand is a protein, you can simply use its molecule name, e.g. 'Mol B')
ligand = 'Res SGN IDS'

# Selection of receptor residues
# (if the ligand is also a protein, replace with the receptor's molecule name, e.g. 'Mol A')
receptor = 'Protein'

# Simulation temperature
temperature = '298K'

# The strength of the pulling force in picoNewton
# This must be large enough to pull off the ligand
strength = 5000.

# Time in picoseconds when the pulling starts
# (everything before can be considered equilibration)
pullstart = 1

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

Console Off

# 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

# Make sure that ligand and receptor are really present
for type in 'ligand','receptor'
  (type)atoms = CountAtom ((type))
  if !(type)atoms
    RaiseError 'This macro requires that the (type) is selected explicitly. The current (type) ((type)) was not found'

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

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

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

# Run the steered simulation
while 1
  t = Time
  t=t/1000
  if t>pullstart
    # Steer the simulation
    # Get the geometric centers of receptor and ligand
    rpos() = PosAtom (receptor)
    lpos() = PosAtom (ligand)
    # Calculate the normalized direction vector from receptor to ligand
    for i=1 to 3
      dir(i)=lpos(i)-rpos(i)
    dis=sqrt (dir1*dir1+dir2*dir2+dir3*dir3)
    ShowMessage 'Steering simulation. Current ligand - receptor distance is (0.00+dis) A'
    scale=(1./dis)*strength
    # Apply the same force, but in opposite directions on receptor and ligand
    # This requires to consider the number of atoms in receptor and ligand
    # NOTE that this recipe only works as long the complex doesn't cross a periodic boundary
    lscale=scale/ligandatoms
    rscale=-scale/receptoratoms
    ForceAtom (ligand),X=(dir1*lscale),Y=(dir2*lscale),Z=(dir3*lscale)
    ForceAtom (receptor),X=(dir1*rscale),Y=(dir2*rscale),Z=(dir3*rscale)
    # Let YASARA simulate one step per screen update/Wait
    SimSteps 1
  else
    ShowMessage 'Equilibration, steered simulation starts in (0.00+pullstart-t) picoseconds...'  
  # Proceed with one simulation step
  Wait 1