#!/usr/bin/env python3
  
import sys
import getopt
import pyhma 
 
try:
  opts, args = getopt.getopt(sys.argv[1:],'rv',['pressure_qh=', 'steps_eq=', 'steps_tot=', 'blocksize=', 'force_tol=', 'meV', 'raw_files', 'verbose'])
except:
  print('Usage: pyhma --pressure_qh=quasiharmonic pressure (GPa) --steps_eq=equilibaration steps --blocksize=block size [--steps_tot=total steps] [--force_tol=force tolerance] [--raw_files|-r] [--meV] [--verbose|-v] vasprun-1.xml vasprun-2.xml ...\n')
  raise
    
filenames = args

pressure_qh = 0      # *required*
steps_eq    = 0      # *required*
blocksize   = 0      # *required*

force_tol   = 0.001  # optional
raw_files   = False  # optional
verbose     = False  # optional
steps_tot   = None   # optional (default is total MD steps in vasprun.xml)
meV         = False  # optional

for opt, val in opts:
  if opt == '--pressure_qh':
    pressure_qh = float(val)
  elif opt == '--steps_eq':   
    steps_eq = int(val)
  elif opt == '--steps_tot':
    steps_tot = int(val) 
  elif opt == '--blocksize':
    blocksize = int(val)
  elif opt == '--force_tol':   
    force_tol = float(val)
  elif opt == '--meV':
    meV = True
  elif opt == '--raw_files' or opt == '-r': 
    raw_files = True
  elif opt == '--verbose' or opt == '-v': 
    verbose = True

if len(args) == 0 or pressure_qh == 0 or steps_eq == 0 or blocksize == 0:
  print('Usage: pyhma --pressure_qh=quasiharmonic pressure (GPa) --steps_eq=equilibaration steps --blocksize=block size [--steps_tot=total steps] [--force_tol=force tolerance] [--raw_files|-r] [--meV] [--verbose|-v] vasprun-1.xml vasprun-2.xml ...\n')
  sys.exit(1)

# Read MD simulation data from vasprun.xml files
data = pyhma.read(filenames, force_tol=force_tol, raw_files=raw_files, verbose=verbose) # a dictionary of data

# Creat simulation object
proc = pyhma.Processor(data, pressure_qh=pressure_qh, meV=meV)

# Compute anharmonic energy and pressure (Conv and HMA) at each step
proc.process(verbose=verbose, steps_tot=steps_tot)
# Get statistics using block averaging method
stats = proc.get_stats(steps_eq=steps_eq, blocksize=blocksize, verbose=verbose)
proc.print_stats(stats)

