# Dielectric Constant Calculation

from numpy import *
import os


global HomeDir
HomeDir = os.getcwd()
print(HomeDir)
field_intensity = 0.1028
epsilon_0 = 8.854187817
save = open("DielectricConstants", "at")


def pbc(length, boundary):
    if length > boundary/2.00:
        length -= boundary
    return length

def read_length(system_name, mode):

    source = open(HomeDir + "/" + mode + "/" + system_name + "-P/CONTCAR", 'r')
    # PBC_on = 0
    for i in range(5):
        line = source.readline()
    size_z = float(line.split()[2])
    for i in range(4):
        line = source.readline()
    z_max = 0.0
    z_min = 1.0
    while True:
        line = (source.readline()).strip("\n")
        line = line.split()
        if len(line) < 3:
            break
        if z_max < pbc(float(line[2]), 1):
            z_max = pbc(float(line[2]), 1)
        if z_min > pbc(float(line[2]), 1):
            z_min = pbc(float(line[2]), 1)
    coordMax = size_z*float(z_max)
    coordMin = size_z*float(z_min)
    if coordMin < 0:
        # PBC_on = 1
        coordMin += size_z

    return coordMin, coordMax

def read_cell_parameters(system_name, mode):
    # mode is "Optical" or "Static"

    file = HomeDir + "/" + mode + "/" + system_name + "-P/LOCPOT"
    source0 = open(file, 'r')

    source0.readline()
    latt_parameter = float(source0.readline())
    line = source0.readline().split()
    a = [float(line[0])*latt_parameter, float(line[1])*latt_parameter, float(line[2])*latt_parameter]
    line = source0.readline().split()
    b = [float(line[0])*latt_parameter, float(line[1])*latt_parameter, float(line[2])*latt_parameter]
    line = source0.readline().split()
    c = [float(line[0])*latt_parameter, float(line[1])*latt_parameter, float(line[2])*latt_parameter]
    bxc = [b[1]*c[2]-b[2]*c[1], b[2]*c[0]-b[0]*c[2], b[0]*c[1]-b[1]*c[0]]
    axb = [a[1]*b[2]-a[2]*b[1], a[2]*b[0]-a[0]*b[2], a[0]*b[1]-a[1]*b[0]]
    volume = a[0]*bxc[0]+a[1]*bxc[1]+a[2]*bxc[2]
    surface = sqrt(pow(axb[0], 2)+pow(axb[1], 2)+pow(axb[2], 2))

    while True:
        line = source0.readline().split()
        if len(line) < 1:
            break

    line = source0.readline().split()
    nx = int(line[0])
    ny = int(line[1])
    nz = int(line[2])
    total_n = nx*ny*nz
    elemental_surface = surface/(nx*ny)

    parameters = {'vectors': [a, b, c], 'volume': [volume], 'surface': surface, "grid_size": [nx, ny, nz],
                  'number_of_points': total_n, 'surface_element': elemental_surface}
    return parameters


def subtract_potentials(system_name, mode):

    pos = open(HomeDir + "/" + mode + "/" + system_name + "-P/LOCPOT", 'r')
    neg = open(HomeDir + "/" + mode + "/" + system_name + "-N/LOCPOT", 'r')

    while True:
        line = pos.readline().split()
        if len(line) < 1:
            break
    pos.readline()

    while True:
        line = neg.readline().split()
        if len(line) < 1:
            break
    neg.readline()

    parameters = read_cell_parameters(system_name, mode)
    total_n = parameters['number_of_points']
    nx = parameters["grid_size"][0]
    ny = parameters["grid_size"][1]
    nz = parameters["grid_size"][2]
    points_field = zeros(total_n)

    n = 0
    while True:
        line = pos.readline().split()
        if n == total_n:
            break
        else:
            for i in line:
                points_field[n] = float(i)
                n += 1

    n = 0
    while True:
        line = neg.readline().split()
        if n == total_n:
            break
        else:
            for i in line:
                points_field[n] -= float(i)
                n += 1

    potential_difference = zeros((nx, ny, nz))
    n = 0
    for k in range(nz):
        for j in range(ny):
            for i in range(nx):
                potential_difference[i][j][k] = points_field[n]
                n += 1

    pos.close()
    neg.close()

    a = parameters['vectors'][0]
    b = parameters['vectors'][1]
    c = parameters['vectors'][2]
    grid = [[[[0, 0, 0] for x in range(nz)] for y in range(ny)] for z in range(nx)]
    for i in range(nx):
        for j in range(ny):
            for k in range(nz):
                grid[i][j][k] = [i*a[0]/nx, j*b[1]/ny, k*c[2]/nz]
                # print(i,j,k,grid[i][j][k])

    return potential_difference, grid

def potential_difference(potential, properties, z_a, z_b):

    length = properties['vectors'][2][2]
    nx = properties['grid_size'][0]
    ny = properties['grid_size'][1]
    nz = properties['grid_size'][2]
    element = length/nz

    # if PBC_on:
    #     k_a = int((z_a+length)/element)
    #     k_b = int(z_b/element)
    # else:
    #     k_a = int(z_a/element)
    #     k_b = int(z_b/element)
    k_a = int(z_a/element)
    k_b = int(z_b/element)

    aver = 0
    potential_difference = zeros((nx, ny))
    for i in range(0, nx):
        for j in range(0, ny):
            # linear interpolation between the grid points
            potential_z_a = (potential[i, j, k_a+1] - potential[i, j, k_a])/element # the slope between the two points
            potential_z_a = potential_z_a*(z_a-(k_a*element)) + potential[i, j, k_a]
            potential_z_b = (potential[i, j, k_b+1] - potential[i, j, k_b])/element # the slope between the two points
            potential_z_b = potential_z_b*(z_b-(k_b*element)) + potential[i, j, k_b]
            potential_difference[i][j] = potential_z_a - potential_z_b
            aver += potential_difference[i][j]
    aver = aver/(nx*ny)

    return potential_difference, aver

def local_dielectric(potential, properties, field_intensity=0.1028):

    length = properties['vectors'][2][2]
    nx = properties['grid_size'][0]
    ny = properties['grid_size'][1]
    nz = properties['grid_size'][2]
    element = length/nz

    dielectric_z = zeros(nz)
    dielectric_z_1D = zeros(nz)
    for k in range(0, nz-1):
        aver_pot = 0
        for i in range(0, nx):
            for j in range(0, ny):
                potential_difference = potential[i][j][k] - potential[i][j][k+1]
                aver_pot += potential_difference
                dielectric_x_y = element*field_intensity/potential_difference
                dielectric_z[k] += dielectric_x_y
        dielectric_z[k] = dielectric_z[k]/(nx*ny)
        aver_pot = aver_pot/(nx*ny)
        dielectric_z_1D[k] = element*field_intensity/aver_pot

    # return dielectric_z, dielectric_z_1D
    return dielectric_z_1D
    # Enter your system name(s) here
SYSTEM_NAME_DICTIONARY = {'Arth-meta1', 'Arth-meta2', 'Arth-para', 'Phen-meta1', 'Phen-meta2', 'Phen-para'}
for SYSTEM_NAME in SYSTEM_NAME_DICTIONARY:
    # sys_name = input("Enter the system name: ")
    # sys_mode = input("Enter the mode (Static/Optical): ")
    sys_name = SYSTEM_NAME
    sys_mode = "Static"
    z_a, z_b = read_length(sys_name, sys_mode)
    print(z_a, z_b)
    properties = read_cell_parameters(sys_name, sys_mode)
    print(properties)

    potential_subtracted, grid = subtract_potentials(sys_name, sys_mode)
    potential_diff, av_potential_diff = potential_difference(potential_subtracted, properties, z_a, z_b)

    print("\nComputation of the medium dielectric constant")
    aver = 0
    #If the dielectric constant < 0, that is because this line below.
    length = abs(z_b - z_a)
    epsilon_min = [0.0, 0.0, 500.0]
    epsilon_max = [0.0, 0.0, 0.0]
    for i in range(0,len(potential_diff)):
        for j in range(0,len(potential_diff[0])):
            dielectric_x_y = length*field_intensity/potential_diff[i][j]
            if dielectric_x_y > epsilon_max[2]:
                epsilon_max = [grid[i][j][0][0], grid[i][j][0][1], dielectric_x_y]
            elif dielectric_x_y < epsilon_min[2]:
                epsilon_min = [grid[i][j][0][0], grid[i][j][0][1], dielectric_x_y]
            aver += dielectric_x_y
            print("%.3f\t%.3f\t%.4f" % (grid[i][j][0][0],grid[i][j][0][1], dielectric_x_y))
    dielectric_constant = aver/(len(potential_diff)*len(potential_diff[0]))
    dielectric_constant_1D = length*field_intensity/av_potential_diff
    capacitance_per_surface = epsilon_0*dielectric_constant/length
    print("-------------------------------------------------------------------")
    print("\nDielectric constant of the medium between %.3f and %.3fA : %.4f" % (z_a,z_b,dielectric_constant))
    print("1D averaged potential dielectric between %.3f and %.3fA : %.4f" % (z_a,z_b,dielectric_constant_1D))
    print("Minimal dielectric is %.4f at (%.3f,%.3f) ; Maximal dielectric is %.4f at (%.3f,%.3f)" % (epsilon_min[2],epsilon_min[0],epsilon_min[1],epsilon_max[2],epsilon_max[0],epsilon_max[1]))