# YASARA MACRO
# TOPIC:       3. Molecular Dynamics
# TITLE:       Running a molecular dynamics simulation of a protein on top of a metal surface with normal or fast speed
# REQUIRES:    Dynamics
# AUTHOR:      Elmar Krieger & Kornel Ozvoldik
# LICENSE:     GPL
# DESCRIPTION: This macro sets up and runs a simulation of a protein on top of a monoelemental face-centerd cubic crystal metal surface

# Parameter section - adjust as needed, but NOTE that some changes only take
# effect if you start an entirely new simulation, not if you continue an existing one.
# ====================================================================================

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

# Extension of the cell on each side of the protein in the surface plane (=XZ plane)
# '15' means that the surface will be 30 A larger than the protein
surfextension=15

# Extension of the cell on each side of the protein along the third (water) axis (=Y-axis)
# '10' means that the cell will be 20 A higher than the solute and surface together
waterextension=10

# Flag to use a square surface. This makes sure that also elongated proteins
# on top of the metal surface can rotate freely during very long simulations.
# If only a short simulation is planned, it can be speeded up by setting
# the flag to 0, creating a rectangular surface that fits the solute more tightly.
square=1

# Surface composition: available Elements are Al, Ni, Cu, Pd, Ag, Pt, Au and Pb.
surfelement='Au'

# Number of monoatomic layers making up the metal surface
surflayers=5

# pH at which the simulation should be run, by default physiological pH 7.4.
ph=7.4

# The ion concentration as a mass fraction, here we use 0.9% NaCl (physiological solution)
ions='Na,Cl,0.9'

# Forcefield to use (this is a YASARA command, so no '=' used)
ForceField AMBER14

# Simulation temperature, which also serves as the random number seed (see Temp command).
# If you increase the temperature significantly by X%, you also need to reduce the timestep by X%
# by changing the 'tslist' that matches your speed below.
temperature='298K'

# Pressure at which the simulation should be run [bar].
pressure=1 

# Cutoff
cutoff=8

# Delay for animations, 1=maximum speed
delay=100

# The format used to save the trajectories: YASARA 'sim', GROMACS 'xtc' or AMBER 'mdcrd'.
# If you don't pick 'sim', a single *.sim restart file will be saved too, since the other
# two formats don't contain velocities, only positions.
format='sim'

# Duration of the complete simulation.
# Alternatively use e.g. duration=5000 to simulate for 5000 picoseconds
# 'if !count duration' simply checks if variable 'duration' as been defined previously (e.g. by an including macro)
if !count duration
  duration='forever'

# The simulation speed, either 'slow' (2*1 fs timestep), 'normal' (2*1.25 fs timestep) or
# 'fast' (maximize performance with 2*2.5 fs timestep and constraints)
# Do not use 'fast' if you simulate incorrect molecules (that would not be stable in reality) 
# 'if !count speed' simply checks if variable 'speed' as been defined previously (e.g. by an including macro) 
if !count speed
  speed='normal'

# The save interval for snapshots. Normally you don't need more than 500-1000 snapshots
# of your simulation, since that's the resolution limit of a typical figure in a journal.
if speed=='fast'
  # Fast speed, save simulation snapshots every 250000 fs, i.e. 250 ps.
  saveinterval=250000
else  
  # Slow or normal speed, save simulation snapshots every 100000 fs, i.e. 100 ps.
  saveinterval=100000

# Flag if the protein's surface attachment needs to be confirmed
confirmneeded=1

# Normally no change required below this point
# ============================================

RequireVersion 24.10.1

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

# Metal surface simulations are always periodic
Boundary periodic

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

# When run as a macro in text mode, add configuration details to log file
if runWithMacro and ConsoleMode
  Processors

Clear
Console off
SimSteps 1,1
# Do we already have a scene with water?
waterscene = FileSize (MacroTarget)_water.sce
if waterscene
  LoadSce (MacroTarget)_water
else
  # No scene with water present yet
  # Do we have a scene with either the protein on top of a metal surface or a user-supplied initial structure that should not be modified?
  scene = FileSize (MacroTarget).sce
  if scene
    # Yes, a scene is already present. Already with surface, or just the initial structure?
    LoadSce (MacroTarget)
    # Search for the surface
    c = CountObj Surface
    if c!=1
      # No surface object, this must be a user-provided initial scene to prevent modifications like Clean/OptHyd, load later
      scene=0
  if !scene
    # No surface scene present yet
    # Load the protein, assuming it's a SCE, PDB or YOB file
    for filetype in 'sce','yob','pdb'
      size = FileSize (MacroTarget).(filetype)
      if size
        break
    if !size
      RaiseError 'Initial structure not found, expected (MacroTarget).pdb or .yob. Make sure to create a project directory and place the structure there in PDB or YOB format'
    # Load structure
    Load(filetype) (MacroTarget)
    DelObj SimCell
    NameObj all,Solute
    cleaned=(filetype!='sce')
    # In case user accidentally provided a YOb file with selected atoms
    Unselect
    Style Ribbon,Stick
    ShowAtom fixed
    # Simple default orientation as a starting point for manual adjustment
    NiceOri
    if ConsoleMode and confirmneeded
      RaiseError "The protein orientation needs to be confirmed visually. Please start YASARA in graphics mode, then continue in text mode"
    if cleaned
      # Delete long peptide bonds that bridge gaps in the structure, which tells CleanAll to add ACE/NME
      # capping groups (the structure of the missing residues could also be predicted, see LoadPDB docs).
      DelBond N,C,LenMin=5
      # Delete waters that are not involved in metal binding and not fixed, to help the calculation of binding energies 
      DelRes !fixed Water with 0 arrows to all
      # Make sure that all hydrogens are present etc. and adapt protonation states to chosen pH
      CleanAll
      pH (ph)
      # pH command recreated water hydrogens, free them
      FreeAtom Element H
      if Structure
        # Optimize the hydrogen-bonding network (more stable trajectories)
        ShowMessage "Optimizing hydrogen bonding network..."
        Wait 1
        OptHydAll
        HideMessage
    # Set cell
    Cell auto,Extension=0
    size1,size2,size3=Cell
    if square
      # Make sure that X- and Z-axis have the same length and we get a square surface
      if size1>size3
        size3=size1
      else
        size1=size3
    for i=1 to 3 step 2
      size(i)=size(i)+surfextension*2
    Cell (size)
    # Build surface and add it to the cell
    BuildSurface surfelement,surflayers,1
    # Extend cell along water axis
    size1,size2,size3=Cell
    size2=size2+waterextension*2
    Cell Y=(size2)
    # Let the solute fly in
    SwitchObj Solute,Off
    UserInput Off
    ZoomObj Simcell,Steps=(delay)
    Wait (delay)
    ShowMessage "Moving solute to suggested position..."
    solutelist()=ListObj Solute
    for i=1 to count solutelist
      posori(0000+i)()=PosOriObj (solutelist(i))
    PosOriObj Solute,Y=-40,Z=-18,Alpha=60
    SwitchObj Solute,On
    for i=1 to count solutelist
      AutoPosOriObj (solutelist(i)),(posori(0000+i)),Steps=(delay*2),Wait=No
    Wait (delay*2)
    if confirmneeded
      # Let the user tune the setup
      UserInput On
      HUD Obj
      ShowMessage "Move the solute and adjust the cell size to fine-tune the suggested setup. Click 'Continue' when done."
      Wait ContinueButton
      ShowMessage "Adjusting surface..."
      ZoomObj Simcell,Steps=(delay)
      Wait (delay)
      # User may have adjusted the cell size, rebuild and replace the surface
      pos()=PosObj Surface
      DelObj Surface
      BuildSurface surfelement,surflayers,0
      PosObj Surface,(pos)
      HideMessage
      Wait 1
    SaveSce (MacroTarget)
  # Solvate system with the surface
  # Set simulation parameters
  Cutoff (cutoff)
  Longrange Coulomb
  # Fill the cell with water, keep user-provided water fixed
  FixRes Water
  # Save pKa data in a file if MacroTarget is set
  pkacommand='Pass'
  if count MacroTarget and MacroTarget!=""
    pkacommand='pKaFile (MacroTarget)'
  Experiment Neutralization
    WaterDensity 0.998
    pH (ph)
    (pkacommand)
    Ions (ions)
    Speed Fast
  Experiment On
  Wait ExpEnd
  StickRes Water
  # Save scene with water
  SaveSce (MacroTarget)_water
  HideMessage
Wait 1
# Don't keep selected atoms, LoadXTC/LoadMDCRD would load only selected ones
Unselect

# Choose timestep and activate constraints
if speed=='fast'
  # Fast simulation speed
  # Constrain bonds to hydrogens
  FixBond all,Element H
  # Constrain certain bond angles involving hydrogens
  FixHydAngle all
  # Choose a multiple timestep of 2*2.5 = 5 fs
  # For structures with severe errors, 2*2 = 4 fs is safer (tslist=2,2)
  tslist=2,2.5
else
  # Slow or normal simulation speed
  # Remove any constraints
  FreeBond all,all
  FreeAngle all,all,all
  if speed=='slow'
    # Choose a multiple timestep of 2*1.00 = 2.0 fs
    tslist=2,1.0
  else
    # Choose a multiple timestep of 2*1.25 = 2.5 fs
    tslist=2,1.25
    # With this timestep, atoms may get too fast in very rare circumstances (only
    # in a specific protein, only once every few nanoseconds). The command below
    # slows down atoms moving faster than 13000 m/s. Such a 'random collision' every
    # few nanoseconds has no more impact than the random number seed. You can comment
    # it out for most proteins, or use the smaller timestep with speed 'slow' above:
    Brake 13000
# During minimization update the pairlist every 10 steps
SimSteps Screen=10,Pairlist=10
# Calculate total timestep, we want a float, so tslist2 is on the left side
ts=tslist2*tslist1
# Snapshots are saved every 'savesteps'
savesteps=saveinterval/ts
# Set final simulation parameters
TimeStep (tslist)
Temp (temperature)
Cutoff (cutoff)
Longrange Coulomb
# Make sure all atoms are free to move
FreeAll
# Alread a snapshot/trajectory present?
i=00000
if format=='sim'
  trajectfilename='(MacroTarget)(i).sim'
else  
  trajectfilename='(MacroTarget).(format)'
  restartfilename='(MacroTarget).sim'
running = FileSize (trajectfilename)
if not running
  # Perform energy minimization
  Experiment Minimization
  Experiment On
  Wait ExpEnd
  # And now start the real simulation
  Sim On
else
  # Simulation has been running before
  ShowMessage "Simulation has been running before, loading last snapshot..."
  # 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
    i=i-1
    LoadSim (MacroTarget)(i)
    # Adjust savesteps to save snapshots in the same interval as previously
    if i>0
      t = Time
      savesteps=0+t/(ts*i)
  else
    # Do we have a restart file with atom velocities?
    found = FileSize (restartfilename)
    if found
      # Yes. First determine the savesteps if possible by loading the 2nd XTC/MDCrd snapshot
      last,t = Load(format) (trajectfilename),1
      if !last
        last,t = Load(format) (trajectfilename),2
        savesteps=0+t/ts
      # Then load the restart file
      LoadSim (restartfilename)
    else
      # No restart file found, load the last snapshot in the XTC/MDCrd trajectory
      do
        i=i+1
        last,t = Load(format) (trajectfilename),(i)
        ShowMessage 'Searching (format) trajectory for last snapshot, showing snapshot (i) at (0+t) fs'
        Sim Pause
        Wait 1
      while !last
      savesteps=0+t/(ts*(i-1))
      Sim Continue
HideMessage

# Set temperature and pressure control
TempCtrl Rescale
PressureCtrl Manometer2D,Pressure=(pressure)

# Now the simulation is running, here you can make changes to the force field

# Uncomment to fix certain atoms in space
# FixAtom Backbone Mol B

# Uncomment to add distance constraints
# AddSpring O Res Lys 80,H Res Glu 84,Len=1.9

# Uncomment to modify charges, e.g. let Trp 12 in Mol A lose an electron:
# ChargeRes Trp 12 Mol A,+1

# And finally, make sure that future snapshots are saved
Save(format) (trajectfilename),(savesteps)
if format!='sim'
  # We save an XTC/MDCrd trajectory plus a single Sim restart file
  SaveSim (restartfilename),(savesteps),Number=no

# Keep the protein from diffusing around and crossing periodic boundaries
CorrectDrift on

# Update the pairlist every 10 (CPU) or 25 (GPU) steps
_,_,gpu = Processors
if gpu
  SimSteps Screen=25,Pairlist=25
else
  SimSteps Screen=10,Pairlist=10

if duration=='forever'
  Console On
  if ConsoleMode
    # In the console, we need to wait forever to avoid a prompt for user input
    Wait forever
else
  Console Off
  measurements=0
  # Wait for given number of picoseconds
  do
    # Tabulate properties you want to monitor during the simulation,
    # e.g. the speeds and velocity vectors of atoms 4, 5 and 7:
    #Tabulate SpeedAtom 4 5 7
    # Or apply a pulling force, in this example downwards
    #AccelRes GLI 1027,Y=-10000
    if ConsoleMode
      # In console model, display the ns/day
      Console On
    # Wait for one screen update
    Wait 1
    measurements=measurements+1
    t = Time
  while t<1000.*duration+1
  # Did we create a table with measurements?
  vallist() = Tab Default
  if count vallist
    # Yes, save the table
    SaveTab default,(MacroTarget)_duringsim,Format=Text,Columns=(count vallist/measurements),Header='Insert your own header here'
  Sim Off
# 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

# BUILD SURFACE
# =============
# Build face-centered mono'element' crystal with monoatomic 'layers',
# 'placed' at the bottom of the current cell.
def BuildSurface element,layers,placed
  # Available elements and their optimized VdW radii
  elementlist='Al','Ni','Cu','Pd','Ag','Pt','Au','Pb'
  radiuslist=1.4625,1.2760,1.3080,1.4095,1.4775,1.4225,1.4755,1.7825
  radius=0
  for i=1 to count elementlist
    # Get radius of selected element
    if element==elementlist(i)
      radius=radiuslist(i)
      break
  if !radius
    RaiseError '(element) is not one of the available surface composition elements (elementlist). Please adapt the `surfelement` variable at the beginning of the macro'
  size1,size2,size3=Cell
  # Side lenght of crystal unit cube
  unitlen=sqrt 8*radius
  # Determine if we have an odd number of layers
  odd=layers%2
  # Get number of atoms in each dimension
  for i=1 to 3 step 2
    gridlen(i)=0+size(i)/unitlen
  gridlen2=(layers+odd)/2
  # Build the face-centered crystal from four primitive cubic grids
  for i=1 to 2
    grid(i)=BuildGrid (element),(gridlen),(unitlen)
  MoveAtom Obj (grid2),(unitlen/2),0,(unitlen/2)
  if layers>1
    for i=3 to 4
      grid(i)=BuildGrid (element),(gridlen1),(gridlen2-odd),(gridlen3),(unitlen)
    if odd
      MoveAtom Obj (grid3) (grid4),0,(-unitlen/2)
    MoveAtom Obj (grid3),(unitlen/2),(unitlen/2)
    MoveAtom Obj (grid4),0,(unitlen/2),(unitlen/2)
    JoinObj (grid3) (grid4),(grid2)
  # Join grids and rename
  NameObj (grid1),Surface
  JoinObj (grid2),Surface
  # Rename residues and atoms to get correct force field parameters
  SplitAtom Obj Surface
  NumberRes Obj Surface
  JoinMol Obj Surface
  NameRes Obj Surface,(element)0
  NameAtom Element Au,Au
  # Adjust cell size
  Cell X=(unitlen*(gridlen1)),Z=(unitlen*(gridlen3))
  if placed
    # Put surface at the bottom of the cell and extend the cell
    Cell (unitlen*(gridlen1)),(size2+unitlen*layers/2),(unitlen*(gridlen3))
    pos()=PosObj SimCell
    PosObj Surface,(pos)
    MoveObj SimCell Surface,Y=(-unitlen*layers/4)
    MoveObj Surface,Y=(-size2/2)