# YASARA MACRO
# TOPIC:       3. Molecular Dynamics
# TITLE:       Running the DHFR benchmark to measure MD performance
# REQUIRES:    Dynamics
# AUTHOR:      Elmar Krieger
# LICENSE:     GPL
# DESCRIPTION: This macro runs the DHFR benchmark for about 10 seconds and displays the results

RequireVersion 15.1.1

# Choose between dihydrofolate reductase (dhfr) and satellite tobacco mosaic virus (stmv)
system='dhfr'
#system='stmv'

# See how many threads and GPUs we have  
_,threadsmax,_,gpus = Processors
# Start with all available CPU threads and no GPU (=0), adjust as needed
threads=threadsmax
gpu=0

#gpus=0

Clear
Console off
HUD off
if !ConsoleMode
  # Allow use to interactively choose the system
  # Load Satellite Tobacco Mosaic Virus (STMV)
  LoadSce (YASARADir)/sce/stmv_water
  DelRes Water MG CL
  DelObj SimCell
  Style Ribbon
  ShowAtom NucAcid Sidechain
  SwapObj 1,2
  # Load Dihydro Folate Reductase (DHFR)
  LoadPDB 3qls,Download=Yes
  DelRes Mol A or HetGroup or Gly 805 and Obj 1,Center=Yes
  PosObj 1,X=-28, Y=3, Z=108
  PosObj 2,X=85, Y=15, Z=390
  # Rotate the proteins
  AutoRotateObj all,Y=0.5
  # Show the options
  LabelAll 'Please choose your benchmark',Height=0.5,Color=Yellow,X=0,Y=4,Z=5
  colorlist='Green','Yellow'
  for i=1 to 2
    # Create the atom to mark
    obj1 = BuildAtom Carbon
    BallObj (obj1)
    ColorObj (obj1),(colorlist(i))
    # Create the donut around it
    obj2 = ShowTorus LargeRadius=2,LargeEdges=32,SmallRadius=0.5,SmallEdges=10,Color=(colorlist(i))
    NameObj (obj2),'Option(i)'
    PosObj (join obj),((-1.5+i)*11),-3,12
  ShowMessage "Click the green atom to benchmark dihydrofolate reductase, and yellow for satellite tobacco mosaic virus (for fast systems)."
  UserInput off
  # Wait until the user chooses an option
  do
    Wait 5
    opt = MarkAtom
  while opt<Atoms-2
  # Pick the requested system, depending on which atom was clicked
  if opt==Atoms-1
    system='dhfr'
    AutoRotateObj 4,X=1,Y=2,Z=3
  else
    system='stmv'
    AutoRotateObj 6,X=-1,Y=-2,Z=-3
  Wait 100
  # Move the objects away
  AutoMoveObj 1,X=-1
  AutoMoveObj 2,X=3
  AutoMoveObj 3-6,Y=-0.1
  Wait 100
  StickRadius 23%
  UserInput on

# Prepare comparison of results
if system=='dhfr'
  # 543, AMD Ryzen Threadripper PRO 5975WX and nVIDIA RTX A6000.
  # 454, AMD Ryzen Threadripper PRO 5975WX 56 threads and RTX 4090 GPU.
  # 446,'Intel Core i9 11900K, 3.5 GHz, 16 threads, Geforce RTX 3060',
  # 304, AMD Ryzen Threadripper PRO 5975WX 56 threads and no GPU.
  resultlist=543,'AMD Ryzen Threadripper PRO 5975WX, RTX A6000',
             290,'Intel Core i9 11900K, 3.5 GHz, 16 threads, no GPU',
             290,'Intel Core i7 5960X, 3.6 GHz, 16 threads, Geforce GTX 980',
             172,'Intel Core i7 5960X, 3.6 GHz, 16 threads, no GPU',
             150,'Intel Core i7 4770, 3.4 GHz, 8 threads, Radeon R9 290X',
             125,'Intel Core i7 6700K, 4.0 GHz, 8 threads, no GPU',
             13,'Intel Atom Z3580, 2.3 GHz, 4 threads, no GPU. Android smartphone'
  plotlen=50             
  consolespeedup=1.1
else
  # 1.77,'Intel Core i7 6700K, 4.0 GHz, 8 threads, no GPU'
  # 1.25,'Intel Core i7 4770, 3.4 GHz, 8 threads, no GPU'             
  # Threadripper 3970X seems to have problems with hyperthreading... 
  resultlist=14.40,'AMD Threadripper Pro 5975WX, Geforce RTX A6000',
              6.79,'AMD Threadripper 3970X, 3.7 GHz, 26 threads, Geforce RTX 3090',
              5.02,'AMD Threadripper 3970X, no GPU',
              3.70,'Intel Core i7 5960X, 3.6 GHz, 16 threads, Gefore GTX 980',
              3.70,'AMD Ryzen 9 3900X, 3.8 GHz, 24 threads, no GPU',
              2.88,'AMD Threadripper 1950X, 3.4 GHz, 32 threads, no GPU',
              2.62,'Intel Core i7 5960X, 3.6 GHz, 16 threads, no GPU'
  plotlen=150             
  consolespeedup=1.05
oldresults=count resultlist
# Load the system to benchmark
Clear
LoadSce (YASARADir)/sce/(system)_water
#Pos Z=125


if !ConsoleMode
  # Prepare visualization of results in graphics mode
  LabelAll '(system) Benchmark',Height=1,Color=Yellow,X=0,Y=4.5,Z=5
  MakeImage Chart,1024,1024,TopCol=000080,BottomCol=8080ff
  MakeImageObj Chart,Chart,Width=(plotlen),Height=(plotlen)
  TransferObj Chart,SimCell,Local=Keep
  SwitchObj Chart,off
  # Warn about other programs
  ShowMessage "Please close all other programs, including those running in the background (e.g. MP3 music players)..."
  Wait ContinueButton
  HideMessage
elif runWithMacro
  # When run as a macro in text mode, add configuration details to log file
  Processors
  
# Choose force field with default parameters (PME, 8 A cutoff)
ForceField AMBER14,SetPar=Yes
TempCtrl Rescale
# Enable pressure control
PressureCtrl WaterProbe,Density=0.997
# Choose a timestep of 5 fs
TimeStep 2,2.5
# Apply constraints to reach the timestep above
FixBond Element H,all
FixHydAngle all
bestgpu=0
bestnspday=0
prevthreads=threads
# Typical speedup in console mode as compared to graphics mode
consolespeedup=1.1
# Try various combinations of CPU threads and GPUs
while gpu<=gpus
  Processors CPUThreads=(threads),GPU=(gpu)
  # Determine the simulation steps per screen update (spu) and the duration of the simulation (in fs)
  if system=='dhfr'
    spu=100
    if gpu
      duration=12500
    else
      duration=5000
  elif gpu
    spu=25
    duration=1000
  else 
    spu=10
    duration=500
  if gpu
    message='Running benchmark with (threads) CPU threads and GPU (gpu) of (gpus)'
    SimSteps Screen=(spu),Pairlist=25
  else
    message='Running benchmark with (threads) of (threadsmax) CPU threads and no GPU'
    SimSteps Screen=(spu),Pairlist=10
  Wait 1
  # Run MD
  Sim On
  updates=duration/(spu*5)
  start = SystemTime
  for i=1 to updates
    ShowMessage '(message), (i*100/updates)% complete'
    Wait 1
  end = SystemTime
  Energy all
  Sim Off
  HideMessage
  # Convert from fs/ms to ns/day
  nspday=3.600e0*24*duration/(end-start)
  if !ConsoleMode
    # Add a bit when running in graphics mode to account for graphics delay
    nspday=nspday*consolespeedup
  # Remember best result
  if nspday>bestnspday
    bestnspday=nspday
    bestgpu=gpu
  # Add to results
  results=count resultlist
  text='Your computer with (threads) of (threadsmax) CPU threads and '
  if gpu
    text=text+'GPU (gpu) of (gpus)'
  else
    text=text+'no GPU'
  # Round as needed
  if system=='dhfr'
    resultlist(results+1)=0+nspday
  else
    resultlist(results+1)=0.00+nspday
  resultlist(results+2)=text
  # Choose next method
  if gpu
    # Try next available GPU
    gpu=gpu+1
  else
    # Try with fewer threads, sometimes this helps
    if nspday<bestnspday
      # Reducing threads didn't help, continue with GPU benchmark
      threads=prevthreads
      gpu=1
    else
      prevthreads=threads  
      if threads>32
        threads=threads-4
      elif threads>16 or threads>threadsmax
        threads=threads-2
      else
        threads=threads-1
      if threads==0
        threads=1
        gpu=1
# Set final result
Processors CPUThreads=(threads),GPU=(bestgpu)
# Show a diagnostic message
message=''
if threads<threadsmax-1
  message=message+'Best results were achieved by using only (threads) of (threadsmax) available threads, which is unusual. '+
                  'Either there were still programs running, or this operating system has issues with high performance computing. '
if bestgpu==0 and gpus
  if gpus==1
    message=message+'The available GPU was'
  else
    message=message+'The (gpus) available GPUs were'
  message=message+' unfortunately not fast enough to accelerate simulations, see www.yasara.org/gpu for recommended models. '
if message!=''
  ShowMessage (message)
# Go back to interactive mode with frequent screen updates
SimSteps 10,10
if !ConsoleMode
  # Show result in graphics mode
  LabelAll '(0.00+bestnspday/consolespeedup) ns/day in graphics mode',Height=0.6,Color=Yellow,X=0,Y=-3.5,Z=5
  LabelAll '(0.00+bestnspday) ns/day in console mode [estimated]',Height=0.6,Color=Yellow,X=0,Y=-4.5,Z=5
  # Generate a bar plot for comparison
  FillRect Color=ffffff
  FillRect 16,16,992,992,Color=None
  threads = Processors
  perfmax=max resultlist
  # Draw the bars
  bars=count resultlist/2
  barheight=(830-64)/(bars*2-1)
  colorlist='Black','White'
  for i=1 to bars
    posy=64+(i-1)*2*barheight
    FillRect 64,(posy),((1024-128)*resultlist(i*2-1)/perfmax),(barheight),Color=(300*i/bars)
    for j=1 to 2
      Font Arial,Height=(barheight*0.5),Color=(colorlist(j))
      PosText (72+j*3),(posy+barheight*0.25),Justify=Left
      Print '(resultlist(i*2-1)) ns/day'
    Font Arial,Height=20,Color=(colorlist((i>oldresults/2)+1))
    PosText 72,(posy+barheight+4),Justify=Left
    Print (resultlist(i*2))
  PosText 72,916
  print '(system) Benchmark, PME, 8.0 A VdW cutoff,\n'
  print 'correct atom masses, reproducible trajectory.\n'
  MoveMesh Chart,Z=(-plotlen)
  SwitchObj Chart,on
  if system=='dhfr'
    Sim On
else
  Print 'DHFR benchmark results:'
  Print '-----------------------'
  for i=1 to count resultlist step 2
    Print '(resultlist(i)) ns/day: (resultlist(i+1))'
  # Exit YASARA if this macro was provided as command line argument
  if runWithMacro
    Energy all
    Exit