#!/usr/bin/env python
'''---------------------------------------------------------------------------------------------------------------------
                     Intrinsic-extrinsic stacking fault and gamma surface calculation for FCC lattice

Description: This script computes the stacking fault properties of an FCC crystal. For more details, refer to README.txt

Author:  Subrahmanyam Pattamatta
Email:   lalithasubrahmanyam@gmail.com

Inputs:  Species    --> (required) Atomic species
         Model  --> (required) Extended ID of a KIM Model
         LatConst   --> (required) Equilibrium (zero-pressure, zero-temperature) FCC lattice constant (meters)
         Pressure   --> (optional) Hydrostatic pressure (bars).  If omitted, the pressure is taken to be zero.
                                   If the value specified is non-zero, the lattice constant specified for
                                   LatConst will be used to construct an initial lattice geometry for an NPT
                                   simulation carried out at the specified pressure and temperature of 1e-4
                                   Kelvin from which the actual lattice constant at the specified pressure is
                                   calculated.

Outputs: gamma_us,  frac_us         # Unstable stacking fault energy  (max)
         gamma_isf, frac_isf = 1.0  # Intrinsic stacking fault energy (min)
         gamma_ut,  frac_ut         # Unstable twinning fault energy  (max)
         gamma_esf, frac_esf = 2.0  # Entrinsic stacking fault energy (min)
         FracList  = []             # fractional displacements array
         SFEDList  = []             # stacking fault energy densities (eV/A^2)
         GammaSurf = []             # [frac_along_112, frac_along_110, energy (ev/A^2)] in each row

Supporting modules used:
    make_lammps_input.py  --> Makes LAMMPS input
    dump_edn.py           --> Dumps output in edn format

External resources used:
    LAMMPS executable with KIM API support
---------------------------------------------------------------------------------------------------------------------'''
from __future__ import print_function
import sys
import os
import operator
from dump_edn import *
from make_lammps_input import *
import numpy as np

# Plotting libs
import matplotlib
matplotlib.use('Agg') # Use backend for non-interactive plotting
import pylab as plt
from mpl_toolkits.mplot3d import Axes3D
from matplotlib import cm
from matplotlib.ticker import LinearLocator, FormatStrFormatter

# Function for printing to stderr
def eprint(*args, **kwargs):
    print(*args, file=sys.stderr, **kwargs)

#-----------------------------------------------------------------------------------------------------------------------
#                                                     MAIN
#-----------------------------------------------------------------------------------------------------------------------

# Program Parameter Variables
LatConst_Tol   = 10e-4

# Read the input variables from command line (see definitions in header)
Model = raw_input("Enter the extended ID of a KIM Model:\n")
if not Model:
    errmsg = "Error: No KIM Model name was specified in input. Exiting..."
    eprint(errmsg)
    sys.exit(1)

Species = raw_input("Enter a chemical element:\n")
if not Species:
    errmsg = "Error: No atomic species was specified in input. Exiting..."
    eprint(errmsg)
    sys.exit(1)

LatConst = raw_input("Enter an FCC lattice constant (meters):\n")
if not LatConst:
    errmsg = "Error: No equilibrium lattice constant was specified in input. Exiting..."
    eprint(errmsg)
    sys.exit(1)
else:
    LatConst = float(LatConst)
    if LatConst <= 0.0:
        errmsg = "Error: Lattice constant specified must be positive. Exiting..."
        eprint(errmsg)
        sys.exit(1)

    # Convert from meters to Angstroms
    LatConst = LatConst * 1.e10

try:
    Pressure = raw_input("Enter a hydrostatic pressure consistent with the lattice constant above (bars):\n")
    Pressure = float(Pressure)
    if Pressure == float(0):
        msg = ("\nInfo: Pressure specified as zero in input. Forgoing lattice constant calculation and "
               "proceeding with lattice constant specified.\n")
        print(msg)
except EOFError:
    Pressure = None
    msg = ("\nInfo: No pressure was specified in input. Taking pressure to be equal to zero and proceeding "
           "with the lattice constant specified.\n")
    print(msg)

#-------------------------------------------------------------------------------
#                        Program internal Constants
#-------------------------------------------------------------------------------
N_Layers       = 58      # No. of layers of the slip planes in the periodic cell
N_Twin_Layers  = N_Layers/2
Rigid_Grp_SIdx = 15
Rigid_Grp_EIdx = 45
Gamma_Nx_112   = 50
Gamma_Ny_110   = 50

output_dir     = "./output"                           # Output directory
if not os.path.exists(output_dir):
    os.makedirs(output_dir)

stack_inp_flnm     = output_dir + "/stack.in"         # Input file for LAMMPS
stack_log_flnm     = output_dir + "/stack.log"        # Log file for LAMMPS
stack_data_flnm    = output_dir + "/stack.dat"        # temporary file for lammps output
stack_results_flnm = output_dir + "/results.edn"  # Results file in .edn format for KIM
LAMMPS_command     = "lammps"

#-------------------------------------------------------------------------------
#      Target variables to be calculated (see definitions in header)
#-------------------------------------------------------------------------------
gamma_us  = 0.0
gamma_isf = 0.0
gamma_ut  = 0.0
gamma_esf = 0.0

frac_us  = 0.0
frac_ut  = 0.0
FracList = []
SFEDList = []
Gamma_X_112_frac = [0 + x*1.0/(Gamma_Nx_112-1) for x in range(Gamma_Nx_112)]
Gamma_Y_110_frac = [0 + y*1.0/(Gamma_Ny_110-1) for y in range(Gamma_Ny_110)]
GammaSurf = []

if not Pressure:
    #------------------------------------------------------------------------------
    #                            CASE I - ZERO PRESSURE
    # Either no pressure was specified, or a pressure of zero was specified.
    # Proceed by constructing the FCC lattice using the equilibrium lattice
    # constant given and forming the stacking faults, etc.
    #------------------------------------------------------------------------------
    Pressure = 0.0

else:
    #------------------------------------------------------------------------------
    #                         CASE II - NON-ZERO PRESSURE
    # Use the zero-temperature, zero-pressure equilibrium lattice constant
    # specified in the input to construct the initial FCC lattice for an NPT
    # simulation at 1e-4K and the specified pressure.  After 200,000 timesteps, the
    # length of the supercell along the x direction is parsed from the output and
    # divided by the number of conventional FCC cells along that direction to
    # arrive at the equilibrium lattice constant at the specified pressure.  This
    # lattice constant is then used to construct the lattice geometry for the
    # actual stacking fault calculations.
    #------------------------------------------------------------------------------
    print("Info: A non-zero pressure of %r bar was specified. Computing the corresponding "
          "FCC lattice constant...\n" % Pressure)
    with open(stack_inp_flnm, "w") as fstack:
       InpStr = compute_eq_latconst( Species, Model, LatConst, Pressure, stack_data_flnm )
       fstack.write(InpStr)
    fstack.close()

    # Run the LAMMPS script
    os.system(LAMMPS_command + " -in "+ stack_inp_flnm + " -log " + stack_log_flnm)

    # Read the LAMMPS output file for the lattice constant
    with open(stack_data_flnm) as fstack:
       linelist = fstack.readlines()
       linebuf = linelist[0].split()
       LatConst  = float(linebuf[0])
    fstack.close()

    # delete the output file
    os.system("rm " + stack_data_flnm)
    os.system("rm " + stack_inp_flnm)
    print("Info: Calculated an FCC lattice constant of %r corresponding to pressure %r\n" % (LatConst, Pressure))

#------------------------------------------------------------------------------
#                            COMPUTE GAMMMA SURFACE
#------------------------------------------------------------------------------
print("***********************************************************")
print("              COMPUTING GAMMA SURFACE                      ")
print("***********************************************************")
with open(stack_inp_flnm, "w") as fstack:
   InpStr = setup_problem(Species, Model, N_Layers, LatConst, Pressure, Rigid_Grp_SIdx, Rigid_Grp_EIdx, N_Twin_Layers)
   fstack.write(InpStr)
   InpStr = make_gammasurface_moves(stack_data_flnm,Gamma_Nx_112,Gamma_Ny_110)
   fstack.write(InpStr)

# Run the LAMMPS script
os.system(LAMMPS_command + " -in "+ stack_inp_flnm + " -log " + stack_log_flnm)

# Read the LAMMPS output file
'''-----------------------------------------------------------------------------
File format: stack.dat
Line 1: Header
Line 1+1 to 1+NxPoints*NyPoints: [ 112_frac    110_frac    SFED ]
-----------------------------------------------------------------------------'''
with open(stack_data_flnm) as fstack:
   linelist = fstack.readlines()
   # Discard the header, index = 0
   # Read the data into arrays
   count = 1
   for yIdx in range(1, Gamma_Ny_110+1):
      temp_list_at_each_y = []
      for xIdx in range(1, Gamma_Nx_112+1):
         linebuf = linelist[count].split()
         count = count+1
         # Discard x any y coordinates
         temp_list_at_each_y.append( float(linebuf[2]) )
      GammaSurf.append(temp_list_at_each_y)

# delete the output file
os.system("rm " + stack_data_flnm)
os.system("rm " + stack_inp_flnm)

#------------------------------------------------------------------------------
#                       COMPUTE STACKING FAULT ENERGIES
#------------------------------------------------------------------------------
with open(stack_inp_flnm, "w") as fstack:
   InpStr = setup_problem(Species, Model, N_Layers, LatConst, Pressure, Rigid_Grp_SIdx, Rigid_Grp_EIdx, N_Twin_Layers)
   fstack.write(InpStr)
   InpStr = make_stack_twin_test(stack_data_flnm)
   fstack.write(InpStr)

# Run the LAMMPS script
os.system(LAMMPS_command + " -in "+ stack_inp_flnm + " -log " + stack_log_flnm)

# Read the LAMMPS output file
'''-----------------------------------------------------------------------------
File format: stack.dat
Line 1:             NPoints 1 nincr nincr
Line 1+1 to 1+1+2*nincr:    [ column1 = frac_disp, column2 = SFED ]
Line:                gamma_us
Line:                gamma_isf
Line:                gamma_ut
Line:                gamma_esf
-----------------------------------------------------------------------------'''
with open(stack_data_flnm) as fstack:
   linelist = fstack.readlines()
   linebuf = linelist[0].split()
   size_0  = int(linebuf[2])
   size_1  = int(linebuf[3])
   size_2  = int(linebuf[4])
   # Read the data into arrays
   listend = 1+size_0+size_1+size_2
   for i in range(1, listend):
      linebuf = linelist[i].split()
      FracList.append(float(linebuf[0]))
      SFEDList.append(float(linebuf[1]))
   # Store the isf and esf values
   gamma_isf = SFEDList[size_0+size_1-1]
   gamma_esf = SFEDList[size_0+size_1+size_2-1]

# delete the output file
os.system("rm " + stack_data_flnm)
os.system("rm " + stack_inp_flnm)

#------------------------------------------------------------------------------
#             Refinement to locate the unstable position - gamma_us
#------------------------------------------------------------------------------
# Locate the unstable stacking fault energy
us_rough_Idx, us_rough_val = max(enumerate(SFEDList[0:size_0+size_1-1]), key=operator.itemgetter(1))

SFrac_us = FracList[us_rough_Idx-1]
dFrac_us = FracList[us_rough_Idx] - FracList[us_rough_Idx-1]

# Make input for the refinement
with open(stack_inp_flnm, "w") as fstack:
   InpStr = setup_problem(Species, Model, N_Layers, LatConst, Pressure, Rigid_Grp_SIdx, Rigid_Grp_EIdx, N_Twin_Layers)
   fstack.write(InpStr)
   InpStr = make_refine_us(SFrac_us, dFrac_us, stack_data_flnm)
   fstack.write(InpStr)

# Run the LAMMPS script
os.system(LAMMPS_command + " -in "+ stack_inp_flnm + " -log " + stack_log_flnm)

# Read the Lammps output file
with open(stack_data_flnm) as fstack:
   linelist = fstack.readlines()
   linebuf  = linelist[1].split()
   frac_us  = float(linebuf[0])
   gamma_us = float(linebuf[1])

# delete the output file
os.system("rm " + stack_data_flnm)
os.system("rm " + stack_inp_flnm)

#------------------------------------------------------------------------------
#             Refinement to locate the unstable position - gamma_ut
#------------------------------------------------------------------------------
# Locate the unstable stacking fault energy
ut_rough_Idx, ut_rough_val = max(enumerate(SFEDList[size_0+size_1:size_0+size_1+size_2-1]), key=operator.itemgetter(1))

SFrac_ut = FracList[size_0+size_1+ut_rough_Idx-1]
dFrac_ut = FracList[size_0+size_1+ut_rough_Idx] - FracList[size_0+size_1+ut_rough_Idx-1]

# Make input for the refinement
with open(stack_inp_flnm, "w") as fstack:
   InpStr = setup_problem(Species, Model, N_Layers, LatConst, Pressure, Rigid_Grp_SIdx, Rigid_Grp_EIdx, N_Twin_Layers)
   fstack.write(InpStr)
   InpStr = make_refine_ut(SFrac_ut, dFrac_ut, stack_data_flnm)
   fstack.write(InpStr)

# Run the LAMMPS script
os.system(LAMMPS_command + " -in "+ stack_inp_flnm + " -log " + stack_log_flnm)

# Read the Lammps output file
with open(stack_data_flnm) as fstack:
   linelist = fstack.readlines()
   linebuf  = linelist[1].split()
   frac_ut  = float(linebuf[0])
   gamma_ut = float(linebuf[1])

# delete the output and log files
os.system("rm " + stack_data_flnm)
os.system("rm " + stack_inp_flnm)
os.system("rm " + stack_log_flnm)
os.system("rm kim.log")

#------------------------------------------------------------------------------
#                    PRINT FINAL OUTPUTS TO KIM EDN FORMAT
#------------------------------------------------------------------------------

# Convert pressure to match relevant KIM Property Definitions
CauchyStress = [-Pressure, -Pressure, -Pressure, 0.0, 0.0, 0.0]

with open(stack_results_flnm, "w") as fstack:

   # Global bracket
   fstack.write("[\n")

   #----------------------------------------------------------------------------
   #                  gamma-surface-relaxed-fcc-crystal-npt
   #----------------------------------------------------------------------------
   # Printing Header
   header = ''' {
   "property-id" "tag:staff@noreply.openkim.org,2015-05-26:property/gamma-surface-relaxed-fcc-crystal-npt"
   "instance-id" 1'''
   fstack.write("%s\n" %(header))
   dump_string_array( "species", Species, fstack )
   dump_dbl_scalar( "a", LatConst, fstack, "angstrom" )
   dump_dbl_vector( "cauchy-stress", CauchyStress, fstack, "bar" )
   dump_dbl_vector( "fault-plane-shift-fraction-112", Gamma_X_112_frac, fstack )
   dump_dbl_vector( "fault-plane-shift-fraction-110", Gamma_Y_110_frac, fstack )
   dump_dbl_matrix( "gamma-surface", GammaSurf, fstack, "eV/A^2" )
   fstack.write(" }\n\n")

   #----------------------------------------------------------------------------
   #                  unstable-stacking-energy-fcc-crystal
   #----------------------------------------------------------------------------
   # Printing Header
   header = ''' {
   "property-id" "tag:staff@noreply.openkim.org,2015-05-26:property/unstable-stacking-fault-relaxed-energy-fcc-crystal-npt"
   "instance-id" 2'''
   fstack.write("%s\n" %(header))
   dump_string_array( "species", Species, fstack )
   dump_dbl_scalar( "a", LatConst, fstack, "angstrom" )
   dump_dbl_vector( "cauchy-stress", CauchyStress, fstack, "bar" )
   dump_dbl_scalar( "unstable-stacking-energy", gamma_us, fstack, "eV/A^2" )
   dump_dbl_scalar( "unstable-slip-fraction", frac_us, fstack )
   fstack.write(" }\n\n")

   #----------------------------------------------------------------------------
   #                intrinsic-stacking-fault-energy-fcc-crystal
   #----------------------------------------------------------------------------
   # Printing Header
   header = ''' {
   "property-id" "tag:staff@noreply.openkim.org,2015-05-26:property/intrinsic-stacking-fault-relaxed-energy-fcc-crystal-npt"
   "instance-id" 3'''
   fstack.write("%s\n" %(header))
   dump_string_array( "species", Species, fstack )
   dump_dbl_scalar( "a", LatConst, fstack, "angstrom" )
   dump_dbl_vector( "cauchy-stress", CauchyStress, fstack, "bar" )
   dump_dbl_scalar( "intrinsic-stacking-fault-energy", gamma_isf, fstack, "eV/A^2" )
   fstack.write(" }\n\n")

   #----------------------------------------------------------------------------
   #                   unstable-twinning-energy-fcc-crystal
   #----------------------------------------------------------------------------
   # Printing Header
   header = ''' {
   "property-id" "tag:staff@noreply.openkim.org,2015-05-26:property/unstable-twinning-fault-relaxed-energy-fcc-crystal-npt"
   "instance-id" 4'''
   fstack.write("%s\n" %(header))
   dump_string_array( "species", Species, fstack )
   dump_dbl_scalar( "a", LatConst, fstack, "angstrom" )
   dump_dbl_vector( "cauchy-stress", CauchyStress, fstack, "bar" )
   dump_dbl_scalar( "unstable-twinning-energy", gamma_ut, fstack, "eV/A^2" )
   dump_dbl_scalar( "unstable-slip-fraction", frac_ut, fstack )
   fstack.write(" }\n\n")

   #----------------------------------------------------------------------------
   #                 extrinsic-stacking-fault-energy-fcc-crystal
   #----------------------------------------------------------------------------
   # Printing Header
   header = ''' {
   "property-id" "tag:staff@noreply.openkim.org,2015-05-26:property/extrinsic-stacking-fault-relaxed-energy-fcc-crystal-npt"
   "instance-id" 5'''
   fstack.write("%s\n" %(header))
   dump_string_array( "species", Species, fstack )
   dump_dbl_scalar( "a", LatConst, fstack, "angstrom" )
   dump_dbl_vector( "cauchy-stress", CauchyStress, fstack, "bar" )
   dump_dbl_scalar( "extrinsic-stacking-fault-energy", gamma_esf, fstack, "eV/A^2" )
   fstack.write(" }\n\n")
   #----------------------------------------------------------------------------
   #                     stacking-energy-curve-fcc-crystal
   #----------------------------------------------------------------------------
   # Printing Header
   header = ''' {
   "property-id" "tag:staff@noreply.openkim.org,2015-05-26:property/stacking-fault-relaxed-energy-curve-fcc-crystal-npt"
   "instance-id" 6'''
   fstack.write("%s\n" %(header))
   dump_string_array( "species", Species, fstack )
   dump_dbl_scalar( "a", LatConst, fstack, "angstrom" )
   dump_dbl_vector( "cauchy-stress", CauchyStress, fstack, "bar" )
   dump_dbl_vector( "fault-plane-shift-fraction", FracList, fstack )
   dump_dbl_vector( "fault-plane-energy", SFEDList, fstack, "eV/A^2" )
   fstack.write(" }\n\n")

   # Global bracket
   fstack.write("]")

#------------------------------------------------------------------------------
#        Print stacking energy curve data to file in asymptote format
#------------------------------------------------------------------------------
with open(stack_data_flnm, "w") as fstack:
   fstack.write("Properties = fractional_displacement    eV/A^2\n")
   fstack.write("%.12f %.12f \n" %(frac_us, gamma_us))
   fstack.write("%.12f %.12f \n" %(1.0, gamma_isf))
   fstack.write("%.12f %.12f \n" %(frac_ut, gamma_ut, ))
   fstack.write("%.12f %.12f \n" %(2.0, gamma_esf))
   fstack.write("%d \n" %(size_0+size_1+size_2))
   for i in range(0, size_0+size_1+size_2):
      fstack.write("%.12f %.12f \n" %(FracList[i], SFEDList[i]))

#------------------------------------------------------------------------------
#       Plot stacking energy curve to png and svg using matplotlib
#------------------------------------------------------------------------------
plt.figure(1, (8,6))
plt.plot(FracList, SFEDList, 'ro', markersize=2)
plt.xlim(FracList[0], FracList[-1])
plt.grid('on', linestyle='--', alpha=0.5)
plt.xlabel(r"$\frac{s\,_{[112]}}{a/\sqrt{6}}$", fontsize=15)
plt.ylabel("Stacking fault energy (eV/$\mathrm{\AA}^2$)")
plt.savefig(os.path.join(output_dir, 'stacking-fault-relaxed-energy-curve-fcc-' + Species + '-' + Model + '.png'), dpi=300)
plt.savefig(os.path.join(output_dir, 'stacking-fault-relaxed-energy-curve-fcc-' + Species + '-' + Model + '.svg'))

#------------------------------------------------------------------------------
#         Plot gamma surface to png and svg using matplotlib
#------------------------------------------------------------------------------
# Convert data to numpy for matplotlib
Gamma_X_112_frac, Gamma_Y_110_frac = np.asarray(Gamma_X_112_frac), np.asarray(Gamma_Y_110_frac)
GammaSurf = np.array(GammaSurf)
Gamma_X_112_frac_grid, Gamma_Y_110_frac_grid = np.meshgrid(Gamma_X_112_frac, Gamma_Y_110_frac)

label112 = r"$\frac{s\,_{[112]}}{\sqrt{6}a/2}$"
label110 = r"$\frac{s\,_{[\mathrm{\overline{1}}10]}}{\sqrt{2}a/2}$"
energy_label = r"$\gamma$ (eV/$\mathrm{\AA}^2$)"
labelfontsize = 15
labelpadding3d = 20

plt.figure(figsize=plt.figaspect(0.5)*1.5)
ax_3d = plt.gca(projection='3d')
gamma_surf = ax_3d.plot_surface(Gamma_X_112_frac_grid, Gamma_Y_110_frac_grid, GammaSurf, cmap=cm.bone, linewidth=0, antialiased=True, vmin=0.0)
ax_3d.set_xlabel(label112, labelpad=labelpadding3d, fontsize=labelfontsize)
ax_3d.set_ylabel(label110, labelpad=labelpadding3d, fontsize=labelfontsize)
ax_3d.set_zlabel(energy_label, labelpad=10)

ax_3d.set_xlim(0,1)
ax_3d.set_ylim(0,1)
ax_3d.tick_params(pad=7)

ax_3d.zaxis.set_major_formatter(FormatStrFormatter('%.05f'))
plt.colorbar(gamma_surf, shrink=0.65, aspect=10, pad=0.035)

# Render at one view
ax_3d.view_init(45,-135)
plt.savefig(os.path.join(output_dir, 'gamma-surface-relaxed-fcc-' + Species + '-' + Model + '-view1.png'), bbox_inches='tight', dpi=300)
plt.savefig(os.path.join(output_dir, 'gamma-surface-relaxed-fcc-' + Species + '-' + Model + '-view1.svg'), bbox_inches='tight')

# Rotate and render again in a second image
ax_3d.view_init(45,135)
plt.savefig(os.path.join(output_dir, 'gamma-surface-relaxed-fcc-' + Species + '-' + Model + '-view2.png'), bbox_inches='tight', dpi=300)
plt.savefig(os.path.join(output_dir, 'gamma-surface-relaxed-fcc-' + Species + '-' + Model + '-view2.svg'), bbox_inches='tight')

# Finally, draw the 2d projection of the gamma surface
plt.figure()
ax_2d = plt.gca()
projected_gamma_surf = ax_2d.pcolor(Gamma_X_112_frac_grid, Gamma_Y_110_frac_grid, GammaSurf, cmap=cm.bone)
ax_2d.set_xlabel(label112, fontsize=labelfontsize)
ax_2d.set_ylabel(label110, fontsize=labelfontsize)
plt.colorbar(projected_gamma_surf, shrink=1, aspect=10, label=energy_label)
plt.savefig(os.path.join(output_dir, 'gamma-surface-relaxed-fcc-' + Species + '-' + Model + '-projected.png'), bbox_inches='tight', dpi=300)
plt.savefig(os.path.join(output_dir, 'gamma-surface-relaxed-fcc-' + Species + '-' + Model + '-projected.svg'), bbox_inches='tight')