#! /usr/bin/env python

import numpy as np
from scipy import linalg as LA
from subprocess import call
from scipy.constants import physical_constants as phc

# constants
atmas = phc['atomic mass constant'][0]
eV2J  = phc['electron volt-joule relationship'][0]
Bc    = phc['Boltzmann constant'][0]
hbar  = phc['Planck constant'][0] / (2*np.pi)
H2J   = phc['Hartree energy'][0]
A2B   = phc['Angstrom star'][0] / phc['Bohr radius'][0]
H2eV  = phc['Hartree energy in eV'][0]

# *****************************************
# Inputs
# *****************************************

#### Cell #######################################
Cell = np.array([[     8.6930937287162351,    0.0000000000000000,    0.0000000000000000],
     [4.3465468643581175,    7.5284400065474486,    0.0000000000000000],
     [0.0000000000000000,    0.0000000000000000,   22.2281052000000017]]) * A2B
Cell_ = np.linalg.inv(Cell)
#################################################

Surface = [[0.0000000000000000,  0.0000000000000000,  0.0000000000000000],
  [0.9999999995943938,  0.3333333330434343,  0.0000000000000000],
  [0.0000000003391207,  0.6666666660868685,  0.0000000000000000],
  [0.3333333333825621,  0.0000000000000000,  0.0000000000000000],
  [0.3333333329769488,  0.3333333330434343,  0.0000000000000000],
  [0.3333333337216828,  0.6666666660868685,  0.0000000000000000],
  [0.6666666667651171,  0.0000000000000000,  0.0000000000000000],
  [0.6666666663595109,  0.3333333330434343,  0.0000000000000000],
  [0.6666666671042378,  0.6666666660868685,  0.0000000000000000],
  [0.1111111108302509,  0.1111111105717129,  0.8902105340044884],
  [0.1111111109108336,  0.4444444449434428,  0.8902105340044884],
  [0.1111111105052274,  0.7777777779868771,  0.8902105340044884],
  [0.4444444442128130,  0.1111111105717129,  0.8902105340044884],
  [0.4444444442933957,  0.4444444449434428,  0.8902105340044884],
  [0.4444444438877895,  0.7777777779868771,  0.8902105340044884],
  [0.7777777775953680,  0.1111111105717129,  0.8902105340044884],
  [0.7777777776759507,  0.4444444449434428,  0.8902105340044884],
  [0.7777777772703445,  0.7777777779868771,  0.8902105340044884],
  [0.2222222221466978,  0.2222222224717214,  0.7844846748340899],
  [0.2222222217410916,  0.5555555555151557,  0.7844846748340899],
  [0.2222222224858186,  0.8888888885585899,  0.7844846748340899],
  [0.5555555555292599,  0.2222222224717214,  0.7844846748340899],
  [0.5555555551236466,  0.5555555555151557,  0.7844846748340899],
  [0.5555555558683807,  0.8888888885585899,  0.7844846748340899],
  [0.8888888889118149,  0.2222222224717214,  0.7844846748340899],
  [0.8888888885062087,  0.5555555555151557,  0.7844846748340899],
  [0.8888888892509428,  0.8888888885585899,  0.7844846748340899],
  [0.0000000000000000,  0.0000000000000000,  0.6748213518442441],
  [0.9999999995943938,  0.3333333330434343,  0.6748213518442441],
  [0.0000000003391207,  0.6666666660868685,  0.6748213518442441],
  [0.3333333333825621,  0.0000000000000000,  0.6748213518442441],
  [0.3333333329769488,  0.3333333330434343,  0.6748213518442441],
  [0.3333333337216828,  0.6666666660868685,  0.6748213518442441],
  [0.6666666667651171,  0.0000000000000000,  0.6748213518442441],
  [0.6666666663595109,  0.3333333330434343,  0.6748213518442441],
  [0.6666666671042378,  0.6666666660868685,  0.6748213518442441]]

Surface_b = np.dot( Surface, Cell )

#-----------------------------------------------------------------------------------------------------
def POSCAR_analyze(CONT):

	CONTlines = open(CONT).readlines()
	pos = Surface_b
	for i in range(45, 47):
		pos = np.append( pos, [[ float(x) * A2B for x in CONTlines[i].split()[0:3] ]], axis=0 )

	return pos

#-----------------------------------------------------------------------------------------------------
#-----------------------------------------------------------------------------------------------------
def states_analyze(CONT):

	CONTlines = open(CONT).readlines()
	pos = Surface_b
	for i in range(2, 0, -1):
		pos = np.append( pos, [[ float(x) * A2B for x in CONTlines[-i].split()[1:4] ]], axis=0 )

	return pos

#-----------------------------------------------------------------------------------------------------

######################### end of inputs #############################
def writeinput(pos):
    runner = '/home/nick/source-rev854/RuNNer.serial.x'

    Type = 36 * ["Au"] + ["H"] + ["Cl"]

    g=open('input.data', 'w')
    g.write( "begin \n" )
    g.write( "lattice   16.427566      0.00000000      0.00000000 \n" )
    g.write( "lattice    8.213783     14.226690      0.00000000 \n" )
    g.write( "lattice    0.00000000      0.00000000     42.005031 \n" )
    for i  in range( len( pos ) ):
        g.write( "atom {0:15.8f} {1:15.8f} {2:15.8f} {3:10s} 0.0  0.0  0.0  0.0  0.0\n" .format(pos[i,0], pos[i,1], pos[i,2], Type[i])  )
    g.write( "energy 0.0\n")
    g.write( "charge 0.0\n")
    g.write( "end\n")
    g.close()
    # calculate energy and forces by RuNNer
    myf = open("mytmp", "w")
    call([runner],stdout=myf)
    call(["rm", "nnatoms.out", "nnforces.out", "output.data", "debug.out","mytmp", "input.data"])
    data=np.loadtxt("energy.out",dtype=float)*H2eV
    call(["rm", "energy.out"])

    return data

####################################################################

#POSCAR = open('barrier/CONTCAR_asymp', 'r').readlines()
POS = POSCAR_analyze( 'POSCAR' )
Pid0 = writeinput(POS)
POS = states_analyze( 'states.xyz' )
Pid1 = writeinput(POS)

profile = open('ET.profile', 'r').readlines()
comdat = open('analysis_rcom.dat', 'r').readlines()
Pdist0 = float( profile[2].split()[2] )
Pdist1 = float( profile[-1].split()[2] )
Pau0 = Pdist0 - Pid0
Pau1 = Pdist1 - Pid1
Kau0 = float( profile[2].split()[4] )
Kau1 = float( profile[-1].split()[4] )
ET = ( Pau1 + Kau1 ) - ( Pau0 + Kau0 )

output = open('PostAnalysis.dat', 'a')
output.write( 'Energy transfer (Efinal - Einitial): {0:f}\n'.format( ET ) )