#!/usr/bin/env python3

#
# GRAIN BOUNDARY SYMMETRIC TILT ENERGY RELAXED CUBIC CRYSTAL
#
# Author:      Brandon Runnels
# Date:        June 2019
# Institution: Department of Mechanical and Aerospace Engineering
#              University of Colorado Colorado Springs
# Email:       brunnels@uccs.edu
#
# Description: This file contains a test driver for generating GB
#              energy data given an arbitrary tilt axis, tilt range,
#              crystal structure, material, and potential.
#              The purpose of this file is to interface between the
#              OpenKIM framework and the compute_gb_tilt_energy python
#              script that generates LAMMPS input and runs the LAMMPS
#              tests.
#              Please see compute_gb_tilt_energy.py for additional
#              information.
#
#
# Copyright:   2016-2019 University of Colorado Colorado Springs
#

import sys

from compute_isolated_atom_energy import compute_isolated_atom_energy
from compute_gb_tilt_energy import compute_gb_tilt_energy

#
# Import property defintion keys from test
#

model_name = sys.stdin.readline().replace("\n", "")
print("model_name=" + model_name)
a1 = [int(x) for x in sys.stdin.readline().split(" ")]
print("a1=" + str(a1))
a2 = [int(x) for x in sys.stdin.readline().split(" ")]
print("a2=" + str(a2))
a3 = [int(x) for x in sys.stdin.readline().split(" ")]
print("a3=" + str(a3))
species = sys.stdin.readline().replace("\n", "")
print("species=" + species)
theta_min = float(sys.stdin.readline().replace("\n", ""))
print("theta_min=" + str(theta_min))
theta_max = float(sys.stdin.readline().replace("\n", ""))
print("theta_max=" + str(theta_max))
max_denominator = int(sys.stdin.readline().replace("\n", ""))
print("max_denominator=" + str(max_denominator))
short_name = sys.stdin.readline().replace("\n", "")
print("short_name=" + short_name)
cutoff_distance = float(sys.stdin.readline())
print("cutoff_distance=" + str(cutoff_distance))
x_repeat = int(sys.stdin.readline().replace("\n", ""))
print("x_repeat=" + str(x_repeat))
z_repeat = int(sys.stdin.readline().replace("\n", ""))
print("z_repeat=" + str(z_repeat))
min_cell_height = float(sys.stdin.readline())
print("min_cell_height=" + str(min_cell_height))
offset_top = [float(x) for x in sys.stdin.readline().split(" ")]
print("offset_top=" + str(offset_top))
offset_bottom = [float(x) for x in sys.stdin.readline().split(" ")]
print("offset_bottom=" + str(offset_bottom))
intermediate = sys.stdin.readline()
a_cohesive_energy = [
    float(x)
    for x in intermediate.replace("[", "").replace("]", "").replace(",", "").split(" ")
]
a = a_cohesive_energy[0] * 1e10  # convert meters to angstroms
print("a=" + str(a))
cohesive_energy = (
    a_cohesive_energy[1] * 6.242e18
)  # convert from joules to electron volts
print("cohesive_energy=" + str(cohesive_energy))
mass = float(sys.stdin.readline().replace("\n", ""))
print("mass=" + str(mass))

#
# Compute energy of an isolated atom of this species for this model in eV
#
isolated_atom_energy = compute_isolated_atom_energy(
    output_dir="output", model_name=model_name, species=species
)

print("")
print("")
print("Isolated atom energy: {} eV".format(isolated_atom_energy))

#
# Compute GB energy using comput_gb_tilt_energy python script
#

thetas, relaxed_energies, minimum_distances, sigmas = compute_gb_tilt_energy(
    output_dir="output",
    species=species,
    a1=a1,
    a2=a2,
    a3=a3,
    pair_style="kim " + model_name,
    pair_coeff="* * " + species + " " + species,
    theta_min=theta_min,
    theta_max=theta_max,
    max_denominator=max_denominator,
    lattice_style=short_name,
    lattice_constant=a,
    cutoff_distance=cutoff_distance,
    x_repeat=x_repeat,
    z_repeat=z_repeat,
    min_cell_height=min_cell_height,
    #    offset_top = offset_top,
    #    offset_bottom = offset_bottom,
    num_mpi_processors=1,
    cohesive_energy=cohesive_energy,
    mass=mass,
    isolated_atom_energy=isolated_atom_energy,
    #    dump_format = 'cfg',
    #    return_positions = 'files',
    gzip=False,
    verbose=True,
)

#
# Export results to EDN file
#

edn_out = open("output/results.edn", "w")

edn_out.write("{\n")
edn_out.write(
    '  "property-id" "tag:brunnels@noreply.openkim.org,2016-02-18:property/grain-boundary-symmetric-tilt-energy-relaxed-relation-cubic-crystal" \n'
)
edn_out.write('  "instance-id" 1 \n')
# a
edn_out.write('  "a" { \n')
edn_out.write('    "source-value" ' + str(a) + "\n")
edn_out.write('    "source-unit" "angstrom"\n')
edn_out.write("  }\n")
# basis-atom-coordinates TODO
edn_out.write('  "basis-atom-coordinates" { \n')
if short_name == "fcc":
    edn_out.write(
        '    "source-value" [ [0 0 0] [0 0.5 0.5] [0.5 0 0.5] [0.5 0.5 0] ] \n'
    )
elif short_name == "bcc":
    edn_out.write('    "source-value" [ [0 0 0] [0.5 0.5 0.5] ] \n')
elif short_name == "sc":
    edn_out.write('    "source-value" [ [0 0 0] ] \n')
else:
    raise ("short_name is not of type fcc, bcc, or sc")
edn_out.write("  }\n")
# interface-offset (ignored)
# minimum-atom-separation
edn_out.write('  "minimum-atom-separation" {\n')
edn_out.write('    "source-value" [ ')
for minimum_distance in minimum_distances:
    edn_out.write(str(minimum_distance) + " ")
edn_out.write("] \n")
edn_out.write('    "source-unit" "angstrom"\n')
edn_out.write("  }\n")
# relaxed-grain-boundary-energy
edn_out.write('  "relaxed-grain-boundary-energy" {\n')
edn_out.write('    "source-value" [ ')
for energy in relaxed_energies:
    edn_out.write(str(energy) + " ")
edn_out.write("] \n")
edn_out.write('    "source-unit" "J/m^2"\n')
edn_out.write("  }\n")
# relaxed-interface-positions
# edn_out.write('  "relaxed-interface-positions" { \n')
# edn_out.write('    "source-value" [\n')
# for f in files:
#    edn_out.write('           "'+f.split('output/')[1]+'"\n' )
# edn_out.write('                   ]\n')
# edn_out.write('  }\n')
# short-name
edn_out.write('  "short-name" {\n')
edn_out.write('    "source-value" [ "' + short_name + '" ] \n')
edn_out.write("  }\n")
# sigma
edn_out.write('  "sigma" {\n')
edn_out.write('    "source-value" [ ')
for sigma in sigmas:
    edn_out.write(str(sigma) + " ")
edn_out.write("] \n")
edn_out.write("  }\n")
# space-group (ignored)
# species
edn_out.write('  "species" {\n')
edn_out.write(
    '    "source-value" ["'
    + species
    + '" "'
    + species
    + '" "'
    + species
    + '" "'
    + species
    + '" ]\n'
)
edn_out.write("  }\n")
# tilt-angle
edn_out.write('  "tilt-angle" {\n')
edn_out.write('    "source-value" [ ')
for theta in thetas:
    edn_out.write(str(theta) + " ")
edn_out.write("] \n")
edn_out.write('    "source-unit" "degrees"\n')
edn_out.write("  }\n")
# tilt-axis
edn_out.write('  "tilt-axis" {\n')
edn_out.write(
    '    "source-value" [' + str(a3[0]) + " " + str(a3[1]) + " " + str(a3[2]) + "]\n"
)
edn_out.write("  }\n")
edn_out.write("}\n")
edn_out.close()
# wyckoff-coordinates (ignored)
# wyckoff-multiplicity-and-letter (ignored)
# wyckoff-species

exit(0)