#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
A collection of functions that plot data from cube files

author: Oxana Andriuc

This code is under CC BY-NC-SA license: https://creativecommons.org/licenses/by-nc-sa/4.0/legalcode
"""

import matplotlib
matplotlib.use('TkAgg')
from mpl_toolkits.mplot3d import Axes3D, proj3d
import matplotlib.pyplot as plt
import numpy as np
import matplotlib.cm as cm
import sys
from pylab import NaN
import math
from matplotlib.widgets import Slider
from matplotlib.widgets import RadioButtons
from operator import itemgetter
import ReadCube as rc
import CalculateCube as cc
import gistfile1 as gf
#import time
from scipy import stats


coldict = {  1:'#cecece',   2:'#efffff',   3:'#db5bff',   4:'#ebff7a',   5:'#ffb7be',   6:'#7c7c7c',   7:'#0032e8',   8:'#ff0000',   9:'#c4fcff',  10:'#addbdd',
            11:'#d749ff',  12:'#cee25a',  13:'#d8adb1',  14:'#7b918d',  15:'#ff9400',  16:'#ffe500',  17:'#30ff52',  18:'#8cc9cc',  19:'#d235ff',  20:'#b0c43c',
            21:'#e2e2e2',  22:'#d1d1d1',  23:'#b7b7b7',  24:'#afb0d8',  25:'#ad8dc9',  26:'#8875af',  27:'#6251ff',  28:'#737dc4',  29:'#cc7b4f',  30:'#7d739b',
            31:'#b28488',  32:'#648780',  33:'#ea91ff',  34:'#ffc300',  35:'#a51000',  36:'#6eb4b7',  37:'#b722e2',  38:'#917d1b',  39:'#d2f6f7',  40:'#abd9db',
            41:'#85bcbf',  42:'#66a2a5',  43:'#4e8b8e',  44:'#3e7e82',  45:'#2c6e72',  46:'#1d5f63',  47:'#9dc1e8',  48:'#e8d99d',  49:'#9b696d',  50:'#416860',
            51:'#b25dc6',  52:'#e59900',  53:'#a600b2',  54:'#4f9b9e',  55:'#9212b7',  56:'#72610b',  57:'#95d2e5',  58:'#fffabc',  59:'#ddffe1',  60:'#bbe8c1',
            61:'#a7e2af',  62:'#91e0ba',  63:'#81e2c5',  64:'#6bdbb9',  65:'#64d6b3',  66:'#59cca8',  67:'#4fc6a1',  68:'#4fc6a1',  69:'#4fc682',  70:'#33a041',
            71:'#21842d',  72:'#5cced6',  73:'#5cb7d6',  74:'#4299b7',  75:'#2e7f9b',  76:'#20708c',  77:'#166987',  78:'#0a5a77',  79:'#ffe100',  80:'#bab9b6',
            81:'#99575c',  82:'#4b514f',  83:'#9041a3',  84:'#a56e00',  85:'#6d5941',  86:'#307c7f',  87:'#760696',  88:'#564803',  89:'#7abedd',  90:'#62c0ef',
            91:'#53b9ed',  92:'#399cce',  93:'#268bbf',  94:'#2577bf',  95:'#4436c1',  96:'#6d49d8',  97:'#9248d8',  98:'#ac48d8',  99:'#bd48d8', 100:'#cc48d8', 
           101:'#af36ab', 102:'#af3690', 103:'#b21c71', 104:'#ff9f7c', 105:'#f27f63', 106:'#cc5f47', 107:'#bc3e32', 108:'#b72a1d', 109:'#91180d'}


def PlotMesh(fpath, xlim=[None, None], ylim=[None, None], zlim=[None, None], alpha=0.1, size=0.05, lower_cutoff=None, upper_cutoff=None, axes=False, colour='rainbow', bkg=(1, 1, 1), minmax=True, minmax_range=0.01, cb_lim=None, au=False, density=1, value_type='esp', units=None):
    """
    SLOW
    
    ! ONLY WORKS FOR ORTHONORMAL GRIDS !

    required arguments: path for a cube file (str)
    
    optional arguments: x coordinate limits (list), y coordinate limits (list), z coordinate limits (list), alpha (float), size (float),
               lower cutoff value (None/float), upper cutoff value (None/float), axes (bool), colour (str), background colour (tuple/list/str), minmax (bool), minmax range (float), colour bar limits
               (list), atomic units (bool), density (int), value type (str), units(None/str)

    calls: ExtractData, Axes, ValuesAsDictionary (all from ReadCube)

    returns: -

    this function creates an interactive plot of the cube file values

    it has sliders for the lower cutoff value, upper cutoff value, alpha and size

    it has buttons for turning axes on/off, showing/hiding min and max, colour scheme, background colour

    fpath (string) = path for a cube file
    xlim (list,default=[None, None]) = a list of two values to be used as the lower and upper limit for the x values (by default all the x values from the cube file are plotted)
    ylim (list,default=[None, None]) = a list of two values to be used as the lower and upper limit for the y values (by default all the y values from the cube file are plotted)
    zlim (list,default=[None, None]) = a list of two values to be used as the lower and upper limit for the z values (by default all the z values from the cube file are plotted)
    alpha (float, default=0.1) = transparency value (between 0 and 1)
    size (float, default=0.05) = marker size
    lower_cutoff (float, default=None) = the lower limit of the values to be plotted (by default all the values from the cube file are plotted)
    upper_cutoff (float, default=None) = the upper limit of the values to be plotted (by default all the values from the cube file are plotted)
    axes (boolean, default=False) = if True, show axes; if False, hide axes
    colour (string, default='rainbow') = colour scheme ('rainbow' by default, change it to 'bwr' for blue-white-red) - full list here
    bkg (tuple/list/string, default=(1, 1, 1)) = background colour (RGB tuple/list or name - each value between 0 and 1, white by default)
    minmax (boolean, default=True) = if True, show minima and maxima; if False, hide minima and maxima
    minmax_range (float, default=0.01) = the value range to consider for the minima and maxima (i.e. ±0.01 by default)
    cb_lim (list, default=None) = limits for the colour bar (list of two floats; if not specified, the limits will be the min and max of the values; if the limits do not cover the whole range of values, a warning is printed)
    au (boolean, default=False) = if True, use atomic units (i.e. bohr for coordinates); if False, use Angstroms for coordinates and units for values
    density (integer, default=1) = an integer number n which indicates that the number of points stored along each axis is reduced by the density value n (therefore the total number of points is reduced by n3)
    value_type (string, default='esp') = can be ‘esp’/’potential’ or ‘dens’/’density’ (or any uppercase version of the aforementioned); used in conjunction with the keywords au and units in order to convert the values
                                        to the appropriate units
    units (string, default=None) = the units to be used for the value; currently supports ‘V’ for ESP and ‘A’/‘angstrom’, ’nm’, ’pm’, ’m’ for density (or any lowercase/uppercase version of these). In the case of the density,
                                    the actual units are the specified units ^(-3). The keywords value_type and au override units. If au is True, the values will be in atomic units regardless of the value passed for units.
                                    If value_type is ‘dens’, the values will be read as density values even if the units argument has a valid ESP value (such as ‘V’). The default unit for ESP is V and for density 1/Å3.
    """
    
    
    #########################################
    # Extracting data in convenient formats #
    #########################################
    
    exdata = rc.ExtractData(fpath, au=au, value_type=value_type, density=density, units=units)
    
    origin=exdata[1]
    n_x = exdata[3]
    x_vector=exdata[4]
    n_y=exdata[5]
    y_vector=exdata[6]
    n_z=exdata[7]
    z_vector=exdata[8]
    vals = exdata[12]
    val_units = exdata[13]
    axes_list=rc.Axes(origin=origin, n_x = n_x, x_vector=x_vector, n_y=n_y, y_vector=y_vector,n_z=n_z, z_vector=z_vector)
    [gridpts, d]=rc.ValuesAsDictionary(vals=vals,axes_list=axes_list, origin=origin)[0:2]

    val_array = np.asarray(vals)
    
    l = list(d.items())

    maxv = max(vals)
    minv = min(vals)

    if lower_cutoff is None:
        lower_cutoff = minv
    if upper_cutoff is None:
        upper_cutoff = maxv


    ##########################
    # Interpreting arguments #
    ##########################

    if au:
        l_units = '$a_0$'
    else:
        l_units = r'$\AA$'

    init_bkg = bkg  # copy of the initial background value

    if cb_lim is None:
        cb_lim = [min(val_array), max(val_array)]
    elif cb_lim[0] > min(val_array) or cb_lim[1] < max(val_array):
        print('\nWarning: The colour bar limits you have selected do not cover the whole range of values\n')


    ###################
    # Creating figure #
    ###################

    fig = plt.figure()
    ax = fig.add_subplot(111, projection='3d')
    ax.set_aspect('equal')
    fig.subplots_adjust(bottom=0.35)  # leave space for the sliders

    cb_ax = fig.add_axes([0.85, 0.32, 0.05, 0.6])  # axes for colour bar

    plt.gcf().text(0.02, 0.97, fpath.split('/')[-1], fontsize=7)     

    def ComputeAndPlot(alpha, size, lower_cutoff, upper_cutoff, axes, colour, minmax, bkg):
        ax.set_facecolor(bkg)   # set background colour for 3D plot

        # create colour bar
        cmap = getattr(cm, colour)
        norm = matplotlib.colors.Normalize(vmin=cb_lim[0], vmax=cb_lim[1])
        matplotlib.colorbar.ColorbarBase(cb_ax, cmap=cmap, norm=norm, label=val_units)

        #########################################################################
        # Updating values based on xlim, ylim, zlim, lower_cutoff, upper_cutoff #
        #########################################################################

        x = [i[0] for i in axes_list[0]]
        y = [i[1] for i in axes_list[1]]
        z = [i[2] for i in axes_list[2]]

        if xlim != [None, None]:
            x_low = xlim[0]
            x_up = xlim[1]
            for i in range(len(x)):
                if x[i] < x_low or x[i] > x_up:
                    x[i] = NaN    # if the x coordinate is outside the x limits, set it to NaN
        if ylim != [None, None]:
            y_low = ylim[0]
            y_up = ylim[1]
            for i in range(len(y)):
                if y[i] < y_low or y[i] > y_up:
                    y[i] = NaN    # i f the y coordinate is outside the y limits, set it to NaN
        if zlim != [None, None]:
            z_low = zlim[0]
            z_up = zlim[1]
            for i in range(len(z)):
                if z[i] < z_low or z[i] > z_up:
                    z[i] = NaN    # if the z coordinate is outside the z limits, set it to NaN

        yy, xx, zz = np.meshgrid(y, x, z)    # make coordinate matrices from coordinate vectors

        if lower_cutoff is not None:
            for i in range(len(l)):
                if l[i][1] < lower_cutoff:
                    (x_i, y_i, z_i) = l[i][0]
                    xx[x_i][y_i][z_i] = NaN   # if the value of a point is lower than the lower cutoff, set its x coordinate to NaN
        else:
            lower_cutoff = minv

        if upper_cutoff is not None:
            for i in range(len(l)):
                if l[i][1] > upper_cutoff:
                    (x_i, y_i, z_i) = l[i][0]
                    xx[x_i][y_i][z_i] = NaN   # if the value of a point is higher than the upper cutoff, set its x coordinate to NaN
        else:
            upper_cutoff = maxv


        #################
        # Plotting data #
        #################

        ax.scatter(xx, yy, zz, c=val_array, cmap=cmap, norm=norm, s=size, alpha=alpha)


        #################
        # Updating axes #
        #################
        
        if xlim != [None, None]:
            ax.set_xlim([xlim[0], xlim[1]])
        if ylim != [None, None]:
            ax.set_ylim([ylim[0], ylim[1]])
        if zlim != [None, None]:
            ax.set_zlim([zlim[0], zlim[1]])

        # lists of all the coordinates of all the points that are not excluded after applying the lower cutoff and upper cutoff
        x_pl = []
        y_pl = []
        z_pl = []
        v_pl = []
        for i in range(len(xx)):
            for j in range(len(xx[i])):
                for k in range(len(xx[i][j])):
                    if not math.isnan(xx[i][j][k]):
                        x_pl.append(xx[i][j][k])
                        y_pl.append(yy[i][j][k])
                        z_pl.append(zz[i][j][k])
                        v_pl.append(d[i,j,k])

        # Create cubic bounding box to simulate equal aspect ratio
        max_range = np.array([max(x_pl)-min(x_pl), max(y_pl)-min(y_pl), max(z_pl)-min(z_pl)]).max()
        Xb = 0.5*max_range*np.mgrid[-1:2:2, -1:2:2, -1:2:2][0].flatten() + 0.5*(max(x_pl)+min(x_pl))
        Yb = 0.5*max_range*np.mgrid[-1:2:2, -1:2:2, -1:2:2][1].flatten() + 0.5*(max(y_pl)+min(y_pl))
        Zb = 0.5*max_range*np.mgrid[-1:2:2, -1:2:2, -1:2:2][2].flatten() + 0.5*(max(z_pl)+min(z_pl))
        # Comment or uncomment following both lines to test the fake bounding box:
        for xb, yb, zb in zip(Xb, Yb, Zb):
            ax.plot([xb], [yb], [zb], 'w')

        minimum = min(v_pl)
        maximum = max(v_pl)
        min_list = []
        max_list = []
        for i in range(len(v_pl)):
            if v_pl[i] <= minimum+minmax_range:
                min_list.append([x_pl[i],y_pl[i],z_pl[i]])
            if v_pl[i] >= maximum-minmax_range:
                max_list.append([x_pl[i],y_pl[i],z_pl[i]])
        for a in min_list:
            ax.scatter([a[0]], [a[1]], [a[2]], s=20+size*2, marker='*', color='blue')

        for a in max_list:
            ax.scatter([a[0]], [a[1]], [a[2]], s=20+size*2, marker='*', color='red')

        if axes == True:  # add axis labels
            ax.set_xlabel('x/'+l_units)
            ax.set_ylabel('y/'+l_units)
            ax.set_zlabel('z/'+l_units)
            ax.grid(True)
            ax.w_xaxis.set_pane_color((0.9, 0.9, 0.9, 0.5))
            ax.w_yaxis.set_pane_color((0.9, 0.9, 0.9, 0.5))
            ax.w_zaxis.set_pane_color((0.9, 0.9, 0.9, 0.5))

            if list(bkg) == [0, 0, 0] or type(bkg) == str and bkg.lower() == 'black':
                ax.spines['bottom'].set_color([0.8, 0.8, 0.8])
                ax.spines['top'].set_color([0.8, 0.8, 0.8])
                ax.spines['right'].set_color([0.8, 0.8, 0.8])
                ax.spines['left'].set_color([0.8, 0.8, 0.8])
                ax.tick_params(axis='both', colors=[0.8, 0.8, 0.8], which='both')
                ax.yaxis.label.set_color([0.8, 0.8, 0.8])
                ax.xaxis.label.set_color([0.8, 0.8, 0.8])
                ax.zaxis.label.set_color([0.8, 0.8, 0.8])
            else:
                ax.spines['bottom'].set_color([0, 0, 0])
                ax.spines['top'].set_color([0, 0, 0])
                ax.spines['right'].set_color([0, 0, 0])
                ax.spines['left'].set_color([0, 0, 0])
                ax.tick_params(axis='both', colors=[0, 0, 0])
                ax.yaxis.label.set_color([0, 0, 0])
                ax.xaxis.label.set_color([0, 0, 0])
                ax.zaxis.label.set_color([0, 0, 0])
        else:
            ax.w_xaxis.set_pane_color((1, 1, 1, 0))
            ax.w_yaxis.set_pane_color((1, 1, 1, 0))
            ax.w_zaxis.set_pane_color((1, 1, 1, 0))
            ax.grid(False)
            ax.set_xticks([])
            ax.set_yticks([])
            ax.set_zticks([])
            ax._axis3don = False

    """
    Sliders and buttons

    for RadioButtons, the default value is the first one, so the order of the options is set so that the initial value is the first one
    """

    ComputeAndPlot(alpha, size, lower_cutoff, upper_cutoff, axes, colour, minmax, bkg)

    alpha_slider_ax = fig.add_axes([0.25, 0.15, 0.65, 0.03])
    alpha_slider = Slider(alpha_slider_ax, 'Alpha', 0.1, 1.0, valinit=alpha)

    size_slider_ax = fig.add_axes([0.25, 0.1, 0.65, 0.03])
    size_slider = Slider(size_slider_ax, 'Size', 0.1, 20.0, valinit=size)

    lower_cutoff_slider_ax = fig.add_axes([0.25, 0.2, 0.65, 0.03])
    lower_cutoff_slider = Slider(lower_cutoff_slider_ax, 'Lower cutoff', minv, maxv, valinit=lower_cutoff)

    upper_cutoff_slider_ax = fig.add_axes([0.25, 0.25, 0.65, 0.03])
    upper_cutoff_slider = Slider(upper_cutoff_slider_ax, 'Upper cutoff', minv, maxv, valinit=upper_cutoff)

    axes_ax = fig.add_axes([0.02, 0.83, 0.20, 0.13])
    if axes is False:
        axes_button = RadioButtons(axes_ax, ('axes off', 'axes on'))
    else:
        axes_button = RadioButtons(axes_ax, ('axes on', 'axes off'))

    minmax_ax = fig.add_axes([0.02, 0.69, 0.20, 0.13])
    if minmax:
        minmax_button = RadioButtons(minmax_ax, ('min/max on', 'min/max off'))
    else:
        minmax_button = RadioButtons(minmax_ax, ('min/max off', 'min/max on'))

    # The colour button includes rainbow, bwr and the value passed by the user as a parameter (if different from rainbow and bwr)
    if colour != 'rainbow' and colour != 'bwr':
        col_b_size = 0.16  # if three options, make the height of the colour button larger
    else:
        col_b_size = 0.13  # if two options, make the height of the colour button smaller
    colour_ax = fig.add_axes([0.02, 0.68-col_b_size, 0.2, col_b_size])
    if colour != 'rainbow' and colour != 'bwr':
        colour_button = RadioButtons(colour_ax, (colour, 'rainbow', 'bwr'))
    elif colour == 'rainbow':
        colour_button = RadioButtons(colour_ax, ('rainbow', 'bwr'))
    else:
        colour_button = RadioButtons(colour_ax, ('bwr', 'rainbow'))

    # The background button includes White, Black and the value passed by the user as a parameter (if different from White or Black)
    # If the background colour is specified by the user in RGB, the option displayed on the button is the closest colour name as identified using the code from gistfile1. When changing back to this 
    # option, the initial RGB colour is used, not the one defined by the closest name
    if list(bkg) != [1, 1, 1] and list(bkg) != [0, 0, 0] and not(type(bkg) == str and (bkg.lower() == 'white' or bkg.lower() == 'black')):
        bkg_b_size = 0.16  # if three options, make the height of the background button larger
    else:
        bkg_b_size = 0.13  # if two options, make the height of the background button smaller
    bkg_ax = fig.add_axes([0.02, 0.67-col_b_size-bkg_b_size, 0.2, bkg_b_size])
    if list(bkg) != [1, 1, 1] and list(bkg) != [0, 0, 0] and not(type(bkg) == str and (bkg.lower() == 'white' or bkg.lower() == 'black')):  # if different from white and black (both in RGB and colour name format)
        if type(bkg) == str:
            bkg_button = RadioButtons(bkg_ax, (bkg, 'White', 'Black'))
        else:
            bkg_button = RadioButtons(bkg_ax, (gf.ColorNames.findNearestWebColorName(bkg[0]*255, bkg[1]*255, bkg[2]*255), 'White', 'Black'))
    elif list(bkg) == [1, 1, 1] or type(bkg) == str and bkg.lower() == 'white':
        bkg_button = RadioButtons(bkg_ax, ('White', 'Black'))
    else:
        bkg_button = RadioButtons(bkg_ax, ('Black', 'White'))

    def sliders_on_changed(val):
        nonlocal alpha, size, lower_cutoff, upper_cutoff
        alpha = alpha_slider.val
        size = size_slider.val
        lower_cutoff = lower_cutoff_slider.val
        upper_cutoff = upper_cutoff_slider.val

        ax.cla()
        ComputeAndPlot(alpha, size, lower_cutoff, upper_cutoff, axes, colour, minmax, bkg)
        plt.show()

    def minmax_on_changed(val):
        nonlocal minmax
        b_dict = {'min/max on': True, 'min/max off': False}
        minmax = b_dict[val]

        ax.cla()
        ComputeAndPlot(alpha, size, lower_cutoff, upper_cutoff, axes, colour, minmax, bkg)
        plt.show()

    def axes_on_changed(val):
        nonlocal axes
        b_dict = {'axes on': True, 'axes off': False}
        axes = b_dict[val]

        ax.cla()
        ComputeAndPlot(alpha, size, lower_cutoff, upper_cutoff, axes, colour, minmax, bkg)
        plt.show()

    def colour_on_changed(val):
        nonlocal colour
        colour = colour_button.value_selected

        ax.cla()
        ComputeAndPlot(alpha, size, lower_cutoff, upper_cutoff, axes, colour, minmax, bkg)
        plt.show()

    def bkg_on_changed(val):
        nonlocal bkg
        if bkg_button.value_selected == 'White':
            bkg = [1, 1, 1]
        elif bkg_button.value_selected == 'Black':
            bkg = [0, 0, 0]
        else:
            bkg = init_bkg

        ax.cla()
        ComputeAndPlot(alpha, size, lower_cutoff, upper_cutoff, axes, colour, minmax, bkg)
        plt.show()

    alpha_slider.on_changed(sliders_on_changed)
    size_slider.on_changed(sliders_on_changed)
    lower_cutoff_slider.on_changed(sliders_on_changed)
    upper_cutoff_slider.on_changed(sliders_on_changed)
    colour_button.on_clicked(colour_on_changed)
    axes_button.on_clicked(axes_on_changed)
    minmax_button.on_clicked(minmax_on_changed)
    bkg_button.on_clicked(bkg_on_changed)

    plt.show()

#@profile
def PlotSurface(fpath, factor=1, alpha=1, size=1, axes=False, colour='rainbow', bkg=(1, 1, 1), minmax=True, minmax_range=0.05, cb_lim=None, atoms=True, au=False, curvature=False, cv_lim=None, cv_pts=5, save=False, grey_surf=False, density=1, value_type='esp', units=None):
    """
    required arguments: path for a cube file (str)
    
    optional arguments: factor (float), alpha (float), size (float), axes (bool), colour (str), background colour (tuple/list/str), minmax (bool), minmax range (float), colour bar limits (list),
                        atoms (bool), atomic units (bool), curvature (bool), curvature limits (float), curvature number of points (int), save (bool), grey surface (bool), density (int), value type (str), units (None/str)

    calls: ExtractData, Axes, ValuesAsDictionary, GetVdWPoints (from ReadCube), VdWLaplacian (from CalculateCube)

    returns: -

    this function creates an interactive plot of all the points which reside within +/- a fractional distance from the surface defined in terms of the Van der Waals radii of the atoms in the molecule
    (the fractional distance represents 30% of the distance between two diagonally adjacent points in the grid)

    it has sliders for alpha, size and zoom

    it has buttons for turning axes on/off, showing/hiding min and max, showing/hiding atoms, colour scheme, background colour

    fpath (string) = path for a cube file
    factor (float, default=1) = factor by which to multiply the VdW radii
    alpha (float, default=1) = transparency value (between 0 and 1)
    size (float, default=1) = marker size
    axes (boolean, default=False) = if True, show axes; if False, hide axes
    colour (string, default='rainbow') = colour scheme ('rainbow' by default, change it to 'bwr' for blue-white-red) - full list here
    bkg (tuple/list/string, default=(1, 1, 1)) = background colour (RGB tuple/list or name - each value between 0 and 1, white by default)
    minmax (boolean, default=False) = if True, show minima and maxima; if False, hide minima and maxima
    minmax_range (float, default=0.05) = the value range to consider for the minima and maxima (i.e. ±0.05 by default)
    cb_lim (list, default=None) = limits for the colour bar (list of two floats; if not specified, the limits will be the min and max of the values; if the limits do not cover the whole range of values, a warning is printed)
    atoms (boolean, default=True) = if True, show atoms; if False, hide atoms
    au (boolean, default=False) = if True, use atomic units (i.e. bohr for coordinates); if False, use Angstroms for coordinates and units for values
    curvature (boolean, default=False) = if True, plots curvature points (either based on cv_lim if given, or cv_pts)
    cv_lim (float, default=None) = the upper limit of the curvature for the low curvature points to be plotted (it overrides cv_pts if given)
    cv_pts (integer, default=5) = the number of low curvature points to be plotted (the points will be chosen in increasing order of their absolute curvature
    save (boolean, default=False) = if True, it saves a .txt file with a summary of the results (minimum, maximum and lowest curvature)
    grey_surf (boolean, default=False) = if True, it plots a grey surface just below the van der Waals surface
    density (integer, default=1) = an integer number n which indicates that the number of points stored along each axis is reduced by the density value n (therefore the total number of points is reduced by n3)
    value_type (string, default='esp') = can be ‘esp’/’potential’ or ‘dens’/’density’ (or any uppercase version of the aforementioned); used in conjunction with the keywords au and units in order to convert the values to the
                                        appropriate units
    units (string, default=None) = the units to be used for the value; currently supports ‘V’ for ESP and ‘A’/‘angstrom’, ’nm’, ’pm’, ’m’ for density (or any lowercase/uppercase version of these). In the case of the density,
                                    the actual units are the specified units ^(-3). The keywords value_type and au override units. If au is True, the values will be in atomic units regardless of the value passed for units.
                                    If value_type is ‘dens’, the values will be read as density values even if the units argument has a valid ESP value (such as ‘V’). The default unit for ESP is V and for density 1/Å3.
    """
#    start_time = time.time()
    
    print('\n\n   “Patience is the companion of wisdom.”\n'+'(St. Augustine)'.rjust(41))
    
    def get_abs_third_elem(iterable):
        return abs(iterable[3])

    #########################################
    # Extracting data in convenient formats #
    #########################################
    exdata = rc.ExtractData(fpath, au=au, value_type=value_type, density=density, units=units)
    origin=exdata[1]
    n_x = exdata[3]
    x_vector=exdata[4]
    n_y=exdata[5]
    y_vector=exdata[6]
    n_z=exdata[7]
    z_vector=exdata[8]
    atoms_list=exdata[9]
    val=exdata[11]  #matrix
    vals=exdata[12]  #list
    val_units=exdata[13]
    axes_list=rc.Axes(origin=origin, n_x = n_x, x_vector=x_vector, n_y=n_y, y_vector=y_vector,n_z=n_z, z_vector=z_vector)
    gridpts=rc.ValuesAsDictionary(vals=vals,axes_list=axes_list, origin=origin)[0]
    
    xdc2 = None
    ydc2 = None
    hl_x = None
    hl_y = None
    hl_z = None
    
    [coords,vdw_ind] = rc.GetVdWPoints(axes_list=axes_list,gridpts=gridpts,atoms_list=atoms_list, x_vector=x_vector, y_vector=y_vector, z_vector=z_vector, origin=origin, factor=factor, au=au)[0:2]
    coord_dict={}
    for i in range(len(coords)):
        coord_dict[(coords[i][0],coords[i][1],coords[i][2])]=coords[i][3]
    sorted_c = sorted(coords, key=itemgetter(3))  # sorted list of coordinates by value
    xs, ys, zs, vs = [list(x) for x in zip(*coords)]
    cv_vdw = sorted(cc.VdWLaplacian(axes_list=axes_list, gridpts=gridpts, atoms_list=atoms_list, x_vector=x_vector,y_vector=y_vector, z_vector=z_vector, val=val, factor=factor, vdw_pts=coords,vdw_ind=vdw_ind, au=au)[0], key=get_abs_third_elem)  # sorted curvature for points on vdw surface
    vs2 = np.asarray(vs)

    
    ##########################
    # Interpreting arguments #
    ##########################

    if au:
        l_units = '$a_0$'
    else:
        l_units = r'$\AA$'

    init_bkg = bkg  # copy of the initial background value
    
    if cb_lim is None:
        cb_lim = [min(vs2), max(vs2)]
    elif cb_lim[0] > min(vs2) or cb_lim[1] < max(vs2):
        print('\nWarning: The colour bar limits you have selected do not cover the whole range of values\n')


    ###################
    # Creating figure #
    ###################

    fig = plt.figure()
    ax = fig.add_subplot(111, projection='3d')  # 1x1 grid, 1st subplot, 3D
    ax.set_aspect('equal')   # equal aspect ratio

    fig.subplots_adjust(bottom=0.3)  # leave space for the slider

    cb_ax = fig.add_axes([0.85, 0.25, 0.05, 0.6])  # axes for colour bar

    plt.gcf().text(0.02, 0.97, fpath.split('/')[-1], fontsize=7)     
    
    ########################
    # Additional variables #
    ########################
    
    zoom = 1    
    u = np.linspace(0, 2 * np.pi, 100)
    v = np.linspace(0, np.pi, 100)


    def ComputeAndPlot(alpha, size, axes, colour, minmax, bkg, zoom, atoms, curvature):
        
        ax.set_facecolor(bkg)   # set background colour for 3D plot

        # plot atoms if True
        if atoms:
            for at in atoms_list:
                x_origin = at[2]
                y_origin = at[3]
                z_origin = at[4]
                
                if au:
                    rf = 0.529177  # lengths in atomic units (i.e. bohr)
                else:
                    rf = 1  # lengths in angstroms
                    
                r = 0.4*rc.vdwdict[at[0]]/rf

                x = x_origin + r * np.outer(np.cos(u), np.sin(v))
                y = y_origin + r * np.outer(np.sin(u), np.sin(v))
                z = z_origin + r * np.outer(np.ones(np.size(u)), np.cos(v))

                ax.plot_surface(x, y, z,  rstride=4, cstride=4, color=coldict[at[0]], linewidth=0, alpha=1)
             
        # plot surface
        if grey_surf:
            for at in atoms_list:
                x_origin = at[2]
                y_origin = at[3]
                z_origin = at[4]

                if au:
                    rf = 0.529177  # lengths in atomic units (i.e. bohr)
                else:
                    rf = 1  # lengths in angstroms

                r = 0.9*rc.vdwdict[at[0]]/rf

                x = x_origin + r * np.outer(np.cos(u), np.sin(v))
                y = y_origin + r * np.outer(np.sin(u), np.sin(v))
                z = z_origin + r * np.outer(np.ones(np.size(u)), np.cos(v))

                ax.plot_surface(x, y, z,  rstride=4, cstride=4, color=(0.6,0.6,0.6), linewidth=0, alpha=1)

        # create colour bar
        cmap = getattr(cm, colour)
        norm = matplotlib.colors.Normalize(vmin=cb_lim[0], vmax=cb_lim[1])
        matplotlib.colorbar.ColorbarBase(cb_ax, cmap=cmap, norm=norm, label = val_units)

        # plot points
        ax.scatter(xs, ys, zs, c=vs2, cmap=cmap, norm=norm, s=size, alpha=alpha)

        # Create cubic bounding box to simulate equal aspect ratio
        max_range = np.array([max(xs)-min(xs), max(ys)-min(ys), max(zs)-min(zs)]).max()
        Xb = 0.5*max_range*np.mgrid[-1:2:2, -1:2:2, -1:2:2][0].flatten() + 0.5*(max(xs)+min(xs))
        Yb = 0.5*max_range*np.mgrid[-1:2:2, -1:2:2, -1:2:2][1].flatten() + 0.5*(max(ys)+min(ys))
        Zb = 0.5*max_range*np.mgrid[-1:2:2, -1:2:2, -1:2:2][2].flatten() + 0.5*(max(zs)+min(zs))
        # Comment or uncomment following both lines to test the fake bounding box:
        for xb, yb, zb in zip(Xb, Yb, Zb):
            ax.plot([xb], [yb], [zb], 'w')

        x1, x2 = ax.get_xlim()
        y1, y2 = ax.get_ylim()
        z1, z2 = ax.get_zlim()
#        ax.set_xlim3d(x1/zoom, x2/zoom)   # with this method the grid and the panes remain the same size within the subplot frame
#        ax.set_ylim3d(y1/zoom, y2/zoom)
#        ax.set_zlim3d(z1/zoom, z2/zoom)
        ax.get_proj = lambda: np.dot(Axes3D.get_proj(ax), np.diag([zoom, zoom, zoom, 1]))   # with this method the grid and the panes extend outside of the subplot frame

        if minmax:
            i = 0
            while sorted_c[i][3] <= min(vs2)+minmax_range:
                ax.scatter([sorted_c[i][0]], [sorted_c[i][1]], [sorted_c[i][2]], s=40+size*2, marker="v", color='blue')  # plot all points with minimum value
                i = i+1
                if i>=len(sorted_c[i]):
                    break
            i = len(sorted_c)-1
            while sorted_c[i][3] >= max(vs2)-minmax_range:
                ax.scatter([sorted_c[i][0]], [sorted_c[i][1]], [sorted_c[i][2]], s=40+size*2, marker="^", color='red')  # plot all points with maximum value
                i = i-1
                if i<0:
                    break
        if curvature:
            if cv_lim:
                i = 0
                while abs(cv_vdw[i][3]) <= cv_lim:
                    
                    ax.scatter([cv_vdw[i][0]], [cv_vdw[i][1]], [cv_vdw[i][2]], s=60+size*2, marker="1", color='black')
                    i += 1
            else:
                for i in range(cv_pts):
                    ax.scatter([cv_vdw[i][0]], [cv_vdw[i][1]], [cv_vdw[i][2]], s=60+size*2, marker="1", color='black')

        if axes == True:  # add axis labels
            ax.set_xlabel('x/'+l_units)
            ax.set_ylabel('y/'+l_units)
            ax.set_zlabel('z/'+l_units)
            ax.grid(True)
            ax.w_xaxis.set_pane_color((0.9, 0.9, 0.9, 0.5))
            ax.w_yaxis.set_pane_color((0.9, 0.9, 0.9, 0.5))
            ax.w_zaxis.set_pane_color((0.9, 0.9, 0.9, 0.5))

            if list(bkg) == [0, 0, 0] or type(bkg) == str and bkg.lower() == 'black':
                ax.spines['bottom'].set_color([0.8, 0.8, 0.8])
                ax.spines['top'].set_color([0.8, 0.8, 0.8])
                ax.spines['right'].set_color([0.8, 0.8, 0.8])
                ax.spines['left'].set_color([0.8, 0.8, 0.8])
                ax.tick_params(axis='both', colors=[0.8, 0.8, 0.8], which='both')
                ax.yaxis.label.set_color([0.8, 0.8, 0.8])
                ax.xaxis.label.set_color([0.8, 0.8, 0.8])
                ax.zaxis.label.set_color([0.8, 0.8, 0.8])
            else:
                ax.spines['bottom'].set_color([0, 0, 0])
                ax.spines['top'].set_color([0, 0, 0])
                ax.spines['right'].set_color([0, 0, 0])
                ax.spines['left'].set_color([0, 0, 0])
                ax.tick_params(axis='both', colors=[0, 0, 0])
                ax.yaxis.label.set_color([0, 0, 0])
                ax.xaxis.label.set_color([0, 0, 0])
                ax.zaxis.label.set_color([0, 0, 0])
        else:
            ax.w_xaxis.set_pane_color((1, 1, 1, 0))
            ax.w_yaxis.set_pane_color((1, 1, 1, 0))
            ax.w_zaxis.set_pane_color((1, 1, 1, 0))
            ax.grid(False)
            ax.set_xticks([])
            ax.set_yticks([])
            ax.set_zticks([])
            ax._axis3don = False
        
        if hl_x and hl_y and hl_z:
            ax.scatter([hl_x],[hl_y],[hl_z],marker='o',color='yellow',s=40+size*2)
            plt.gcf().text(0.8, 0.9, "  Highlighted:".center(15)+"\n"+('%.4E' % (coord_dict[(hl_x,hl_y,hl_z)])).center(15),bbox = dict(boxstyle = 'round,pad=0.5', fc = (0.9,0.9,0.9)))     


    """
    Sliders and buttons

    for RadioButtons, the default value is the first one, so the order of the options is set so that the initial value is the first one
    """

    alpha_slider_ax = fig.add_axes([0.3, 0.11, 0.4, 0.03])
    alpha_slider = Slider(alpha_slider_ax, 'Alpha', 0, 1.0, valinit=alpha)

    size_slider_ax = fig.add_axes([0.3, 0.06, 0.4, 0.03])
    size_slider = Slider(size_slider_ax, 'Size', 0.1, 50.0, valinit=size)

    zoom_slider_ax = fig.add_axes([0.3, 0.01, 0.4, 0.03])
    zoom_slider = Slider(zoom_slider_ax, 'Zoom', 0.5, 5.0, valinit=zoom)
    
    elev_slider_ax = fig.add_axes([0.3, 0.16, 0.4, 0.03])
    elev_slider = Slider(elev_slider_ax, 'Elev', 0, 360, valinit=ax.elev)
    
    azim_slider_ax = fig.add_axes([0.3, 0.21, 0.4, 0.03])
    azim_slider = Slider(azim_slider_ax, 'Azim', 0, 360, valinit=ax.azim)

    axes_ax = fig.add_axes([0.02, 0.83, 0.20, 0.13])
    if axes is False:
        axes_button = RadioButtons(axes_ax, ('axes off', 'axes on'))
    else:
        axes_button = RadioButtons(axes_ax, ('axes on', 'axes off'))

    minmax_ax = fig.add_axes([0.02, 0.69, 0.20, 0.13])
    if minmax:
        minmax_button = RadioButtons(minmax_ax, ('min/max on', 'min/max off'))
    else:
        minmax_button = RadioButtons(minmax_ax, ('min/max off', 'min/max on'))

    atoms_ax = fig.add_axes([0.02, 0.55, 0.20, 0.13])
    if atoms:
        atoms_button = RadioButtons(atoms_ax, ('atoms on', 'atoms off'))
    else:
        atoms_button = RadioButtons(atoms_ax, ('atoms off', 'atoms on'))

    # The colour button includes rainbow, bwr and the value passed by the user as a parameter (if different from rainbow and bwr)
    if colour != 'rainbow' and colour != 'bwr':
        col_b_size = 0.16  # if three options, make the height of the colour button larger
    else:
        col_b_size = 0.13  # if two options, make the height of the colour button smaller
    colour_ax = fig.add_axes([0.02, 0.54-col_b_size, 0.2, col_b_size])
    if colour != 'rainbow' and colour != 'bwr':
        colour_button = RadioButtons(colour_ax, (colour, 'rainbow', 'bwr'))
    elif colour == 'rainbow':
        colour_button = RadioButtons(colour_ax, ('rainbow', 'bwr'))
    else:
        colour_button = RadioButtons(colour_ax, ('bwr', 'rainbow'))

    # The background button includes White, Black and the value passed by the user as a parameter (if different from White or Black)
    # If the background colour is specified by the user in RGB, the option displayed on the button is the closest colour name as identified using the code from gistfile1. When changing back to this
    # option, the initial RGB colour is used, not the one defined by the closest name
    if list(bkg) != [1, 1, 1] and list(bkg) != [0, 0, 0] and not(type(bkg) == str and (bkg.lower() == 'white' or bkg.lower() == 'black')):
        bkg_b_size = 0.16  # if three options, make the height of the background button larger
    else:
        bkg_b_size = 0.13  # if two options, make the height of the background button smaller
    bkg_ax = fig.add_axes([0.02, 0.53-col_b_size-bkg_b_size, 0.2, bkg_b_size])
    if list(bkg) != [1, 1, 1] and list(bkg) != [0, 0, 0] and not(type(bkg) == str and (bkg.lower() == 'white' or bkg.lower() == 'black')):  # if different from white and black (both in RGB and colour name format)
        if type(bkg) == str:
            bkg_button = RadioButtons(bkg_ax, (bkg, 'White', 'Black'))
        else:
            bkg_button = RadioButtons(bkg_ax, (gf.ColorNames.findNearestWebColorName(bkg[0]*255, bkg[1]*255, bkg[2]*255), 'White', 'Black'))
    elif list(bkg) == [1, 1, 1] or type(bkg) == str and bkg.lower() == 'white':
        bkg_button = RadioButtons(bkg_ax, ('White', 'Black'))
    else:
        bkg_button = RadioButtons(bkg_ax, ('Black', 'White'))

    curv_ax = fig.add_axes([0.02, 0.39-col_b_size-bkg_b_size, 0.2, 0.13])
    if curvature:
        curv_button = RadioButtons(curv_ax, ('curvature on', 'curvature off'))
    else:
        curv_button = RadioButtons(curv_ax, ('curvature off', 'curvature on'))
        

    # plot with initial parameters
    ax.cla()
    ComputeAndPlot(alpha, size, axes, colour, minmax, bkg, zoom, atoms, curvature)
    
    
    # Functions defining what happens when the values of the sliders and buttons change
    def sliders_on_changed(val):
        nonlocal alpha, size, zoom
        alpha = alpha_slider.val
        size = size_slider.val
        zoom = zoom_slider.val

        ax.cla()
        ComputeAndPlot(alpha, size, axes, colour, minmax, bkg, zoom, atoms, curvature)
        plt.show()
        
    def elev_changed(val):
        elev = elev_slider.val

        ax.cla()
        ComputeAndPlot(alpha, size, axes, colour, minmax, bkg, zoom, atoms, curvature)
        ax.view_init(elev=elev)
        plt.show()
        
    def azim_changed(val):
        azim = azim_slider.val

        ax.cla()
        ComputeAndPlot(alpha, size, axes, colour, minmax, bkg, zoom, atoms, curvature)
        ax.view_init(azim=azim)
        plt.show()
        
    def minmax_on_changed(val):
        nonlocal minmax
        b_dict = {'min/max on': True, 'min/max off': False}
        minmax = b_dict[val]

        ax.cla()
        ComputeAndPlot(alpha, size, axes, colour, minmax, bkg, zoom, atoms, curvature)
        plt.show()

    def axes_on_changed(val):
        nonlocal axes
        b_dict = {'axes on': True, 'axes off': False}
        axes = b_dict[val]

        ax.cla()
        ComputeAndPlot(alpha, size, axes, colour, minmax, bkg, zoom, atoms, curvature)
        plt.show()

    def atoms_on_changed(val):
        nonlocal atoms
        b_dict = {'atoms on': True, 'atoms off': False}
        atoms = b_dict[val]

        ax.cla()
        ComputeAndPlot(alpha, size, axes, colour, minmax, bkg, zoom, atoms, curvature)
        plt.show()

    def colour_on_changed(val):
        nonlocal colour
        colour = colour_button.value_selected

        ax.cla()
        ComputeAndPlot(alpha, size, axes, colour, minmax, bkg, zoom, atoms, curvature)
        plt.show()

    def bkg_on_changed(val):
        nonlocal bkg
        if bkg_button.value_selected == 'White':
            bkg = [1, 1, 1]
        elif bkg_button.value_selected == 'Black':
            bkg = [0, 0, 0]
        else:
            bkg = init_bkg

        ax.cla()
        ComputeAndPlot(alpha, size, axes, colour, minmax, bkg, zoom, atoms, curvature)
        plt.show()

    def curv_on_changed(val):
        nonlocal curvature
        b_dict = {'curvature on': True, 'curvature off': False}
        curvature = b_dict[val]

        ax.cla()
        ComputeAndPlot(alpha, size, axes, colour, minmax, bkg, zoom, atoms, curvature)
        plt.show()
        
    def onclick(event):
        nonlocal xdc2, ydc2, hl_x, hl_y, hl_z
        if event.dblclick:
            ax.button_pressed = None
            s = ax.format_coord(event.xdata,event.ydata)
            
            try:
                [xss, yss, zss] = s.split(",")
                x = float(xss.strip().split("=")[1])
                y = float(yss.strip().split("=")[1])
                z = float(zss.strip().split("=")[1])
                
                xdc2, ydc2, _ = proj3d.proj_transform(x,y,z, ax.get_proj())
                
                coords_2d = []
                for i in range(len(coords)):
                    x2, y2, _ = proj3d.proj_transform(coords[i][0],coords[i][1],coords[i][2], ax.get_proj())
                    coords_2d.append([coords[i][0],coords[i][1],coords[i][2],coords[i][3],x2,y2])

                xy_dc2 = np.asarray([xdc2,ydc2])  # x y double click 2d
                                                
                if xdc2 <= max(coords_2d,key=itemgetter(4))[4] and xdc2 >= min(coords_2d,key=itemgetter(4))[4] and ydc2 <= max(coords_2d,key=itemgetter(5))[5] and ydc2 >= min(coords_2d,key=itemgetter(5))[5]:
                    
                    for i in range(len(coords_2d)):
                        dist = np.sum((np.asarray([coords_2d[i][4],coords_2d[i][5]]) - xy_dc2)**2)
                        coords_2d[i].append(dist)

                    ordered_coords_2d = sorted(coords_2d,key=itemgetter(6))
                    
                    
                    hl_x = ordered_coords_2d[0][0]  # highlighted x 
                    hl_y = ordered_coords_2d[0][1]
                    hl_z = ordered_coords_2d[0][2]
                    ax.cla()
                    ComputeAndPlot(alpha, size, axes, colour, minmax, bkg, zoom, atoms, curvature)
                    plt.show()
            except:
                pass
            
        
    fig.canvas.mpl_connect('button_press_event', onclick)

    colour_button.on_clicked(colour_on_changed)
    axes_button.on_clicked(axes_on_changed)
    minmax_button.on_clicked(minmax_on_changed)
    alpha_slider.on_changed(sliders_on_changed)
    size_slider.on_changed(sliders_on_changed)
    zoom_slider.on_changed(sliders_on_changed)
    bkg_button.on_clicked(bkg_on_changed)
    atoms_button.on_clicked(atoms_on_changed)
    curv_button.on_clicked(curv_on_changed)
    elev_slider.on_changed(elev_changed)
    azim_slider.on_changed(azim_changed)

    print('\n\n '+'-'*74+'\n|'+' Summary '.center(74)+'|\n '+'-'*74)
    print('Minimum: '+'%.4E' % (sorted_c[0][3])+' at (%.3f, %.3f, %.3f)' % (sorted_c[0][0],sorted_c[0][1],sorted_c[0][2]))
    print('Maximum: '+'%.4E' % (sorted_c[len(sorted_c)-1][3])+' at (%.3f, %.3f, %.3f)' % (sorted_c[len(sorted_c)-1][0],sorted_c[len(sorted_c)-1][1],sorted_c[len(sorted_c)-1][2]))
    print('Lowest curvature: '+'%.4E' % (cv_vdw[0][3])+' at (%.3f, %.3f, %.3f)' % (cv_vdw[0][0],cv_vdw[0][1],cv_vdw[0][2])+'\n')
    
#    print("--- %s seconds ---" % (time.time() - start_time))
    plt.show()
    
    if save:
        with open(fpath[:-5]+'_summary.txt','w') as f:
            f.write('\n########## Summary ##########\n\n'+'Minimum: '+str(sorted_c[0][3])+' at (%.3f, %.3f, %.3f)' % (sorted_c[0][0],sorted_c[0][1],sorted_c[0][2])+
                    '\nMaximum: '+str(sorted_c[len(sorted_c)-1][3])+' at (%.3f, %.3f, %.3f)' % (sorted_c[len(sorted_c)-1][0],sorted_c[len(sorted_c)-1][1],sorted_c[len(sorted_c)-1][2])+
                    '\nLowest curvature: '+str(cv_vdw[0][3])+' at (%.3f, %.3f, %.3f)' % (cv_vdw[0][0],cv_vdw[0][1],cv_vdw[0][2]))

def PlotIsosurface(fpath, factor=1, alpha=1, size=4, axes=True, colour='rainbow', bkg=(1, 1, 1), minmax=False, minmax_range=0, cb_lim=None, atoms=True, au=False, curvature=False, cv_lim=None, cv_pts=5, zoom=1, density=1, iso=None, iso_range=0.01, value_type='esp', units=None):
    """
    required arguments: path for a cube file (str)
    
    optional arguments: factor (float), alpha (float), size (float), axes (bool), colour (str), background colour (tuple/list/str), minmax (bool), minmax range (float), colour bar limits (list),
                        atoms (bool), atomic units (bool), curvature (bool), curvature limits (float), curvature number of points (int), zoom (float), density (int), iso (float), iso range (float), 
                        value type (str), units (None/str)

    calls: GetVdWPoints (from ReadCube)

    returns: -

    this function creates an interactive plot of all the points which are within a specified value range from a given value

    it has sliders for alpha, size and zoom

    it has buttons for turning axes on/off, showing/hiding min and max, showing/hiding atoms, colour scheme, background colour
    
    fpath (string) = path for a cube file 
    factor (float, default=1) = factor by which to multiply the VdW radii
    alpha (float, default=1) = transparency value (between 0 and 1)
    size (float, default=4) = marker size
    axes (boolean, default=True) = if True, show axes; if False, hide axes
    colour (string, default='rainbow') = colour scheme ('rainbow' by default, change it to 'bwr' for blue-white-red) - full list here
    bkg (tuple/list/string, default=(1, 1, 1)) = background colour (RGB tuple/list or name - each value between 0 and 1, white by default)
    minmax (boolean, default=True) = if True, show minima and maxima; if False, hide minima and maxima
    minmax_range (float, default=0) = the value range to consider for the minima and maxima
    cb_lim (list, default=None) = limits for the colour bar (list of two floats; if not specified, the limits will be the min and max of the values; if the limits do not cover the whole range of values, a warning is printed)
    atoms (boolean, default=False) = if True, show atoms; if False, hide atoms
    au (boolean, default=False) = if True, use atomic units (i.e. bohr for coordinates); if False, use Angstroms for coordinates and units for values
    curvature (boolean, default=False) = if True, plots curvature points (either based on cv_lim if given, or cv_pts)
    cv_lim (float, default=None) = the upper limit of the curvature for the low curvature points to be plotted (it overrides cv_pts if given)
    cv_pts (integer, default=5) = the number of low curvature points to be plotted (the points will be chosen in increasing order of their absolute curvature
    zoom (float, default=1) = how large the plot is displayed with respect to the default size
    density (integer, default=1) = an integer number n which indicates that the number of points stored along each axis is reduced by the density value n (therefore the total number of points is reduced by n3)
    iso (float, default=None) = the value of the points to be plotted
    iso_range (float, default=0.01) = the value range to consider when selecting the points (i.e. iso±0.01 by default)
    value_type (string, default='esp') = can be ‘esp’/’potential’ or ‘dens’/’density’ (or any uppercase version of the aforementioned); used in conjunction with the keywords au and units in order to convert the values to
                                        the appropriate units
    units (string, default=None) = the units to be used for the value; currently supports ‘V’ for ESP and ‘A’/‘angstrom’, ’nm’, ’pm’, ’m’ for density (or any lowercase/uppercase version of these). In the case of the density,
                                    the actual units are the specified units ^(-3). The keywords value_type and au override units. If au is True, the values will be in atomic units regardless of the value passed for units.
                                    If value_type is ‘dens’, the values will be read as density values even if the units argument has a valid ESP value (such as ‘V’). The default unit for ESP is V and for density 1/Å3.
    """
    
    print('\n\n   “Patience is the companion of wisdom.”\n'+'(St. Augustine)'.rjust(41))
    
    #########################################
    # Extracting data in convenient formats #
    #########################################
    
    exdata = rc.ExtractData(fpath, au=au, value_type=value_type, density=density, units=units)
    origin=exdata[1]
    n_x = exdata[3]
    x_vector=exdata[4]
    n_y=exdata[5]
    y_vector=exdata[6]
    n_z=exdata[7]
    z_vector=exdata[8]
    atoms_list = exdata[9]
    vals=exdata[12]
    val_units=exdata[13]
    axes_list=rc.Axes(origin=origin, n_x = n_x, x_vector=x_vector, n_y=n_y, y_vector=y_vector,n_z=n_z, z_vector=z_vector)
    (gridpts,d_ind)=rc.ValuesAsDictionary(vals=vals,axes_list=axes_list, origin=origin)
    
    vdw_coords = rc.GetVdWPoints(axes_list=axes_list,gridpts=gridpts,atoms_list=atoms_list, x_vector=x_vector, y_vector=y_vector, z_vector=z_vector, origin=origin, factor=factor, au=au)[0]
    
    coords = [list(i)+[gridpts[i]] for i in list(gridpts.keys())]
    sorted_c = sorted(coords, key=itemgetter(3))  # sorted list of coordinates by value
    xs, ys, zs, vs = [list(j) for j in zip(*coords)]
    vs2 = np.asarray(vs)

    x_vdw, y_vdw, z_vdw, v_vdw = [list(j) for j in zip(*vdw_coords)]
    
    ##########################
    # Interpreting arguments #
    ##########################

    if iso is None:
        iso = (min(v_vdw)+max(v_vdw))/2
        
    if au:
        l_units = '$a_0$'
    else:
        l_units = r'$\AA$'

    if cb_lim is None:
        cb_lim = [min(v_vdw), max(v_vdw)]
    elif cb_lim[0] > min(vs2) or cb_lim[1] < max(vs2):
        print('\nWarning: The colour bar limits you have selected do not cover the whole range of values\n')

    ###################
    # Creating figure #
    ###################

    fig = plt.figure()
    ax = fig.add_subplot(111, projection='3d')  # 1x1 grid, 1st subplot, 3D
    ax.set_aspect('equal')   # equal aspect ratio
    fig.subplots_adjust(bottom=0.3)  # leave space for the slider
    cb_ax = fig.add_axes([0.85, 0.25, 0.05, 0.6])  # axes for colour bar
    plt.gcf().text(0.02, 0.97, fpath.split('/')[-1], fontsize=7)     


    ########################
    # Additional variables #
    ########################

    xdc2 = None
    ydc2 = None
    hl_x = None
    hl_y = None
    hl_z = None

    zoom = 1    

    init_bkg = bkg  # copy of the initial background value

    u = np.linspace(0, 2 * np.pi, 100)
    v = np.linspace(0, np.pi, 100)

    def select_pts(iso,iso_range): 
        x_p = []
        y_p = []
        z_p = []
        v_p = []
            
        for i in range(len(sorted_c)):
          
            if sorted_c[i][3]>iso+abs(iso_range):
                break
            if sorted_c[i][3]>=iso-abs(iso_range):
                x_p.append(sorted_c[i][0])
                y_p.append(sorted_c[i][1])
                z_p.append(sorted_c[i][2])
                v_p.append(sorted_c[i][3])
        return [x_p,y_p,z_p,v_p]

    [x_p,y_p,z_p,v_p]=select_pts(iso,iso_range)
    
    
    def ComputeAndPlot(alpha, size, axes, colour, minmax, bkg, zoom, atoms, curvature):
        ax.set_facecolor(bkg)   # set background colour for 3D plot
        # plot atoms if True
        if atoms:
            for at in atoms_list:
                x_origin = at[2]
                y_origin = at[3]
                z_origin = at[4]

                if au:
                    rf = 0.529177  # lengths in atomic units (i.e. bohr)
                else:
                    rf = 1  # lengths in angstroms

                r = 0.4*rc.vdwdict[at[0]]/rf

                x_box = x_origin + r * np.outer(np.cos(u), np.sin(v))
                y_box = y_origin + r * np.outer(np.sin(u), np.sin(v))
                z_box = z_origin + r * np.outer(np.ones(np.size(u)), np.cos(v))

                ax.plot_surface(x_box, y_box, z_box,  rstride=4, cstride=4, color=coldict[at[0]], linewidth=0, alpha=1)


        # create colour bar
        cmap = getattr(cm, colour)
        norm = matplotlib.colors.Normalize(vmin=cb_lim[0], vmax=cb_lim[1])
        matplotlib.colorbar.ColorbarBase(cb_ax, cmap=cmap, norm=norm, label=val_units)

        
        # plot points
        ax.scatter(x_p, y_p, z_p, c=v_p, cmap=cmap, norm=norm, s=size, alpha=alpha)
        
        # Create cubic bounding box to simulate equal aspect ratio
        max_range = np.array([max(xs)-min(xs), max(ys)-min(ys), max(zs)-min(zs)]).max()
        Xb = 0.5*max_range*np.mgrid[-1:2:2, -1:2:2, -1:2:2][0].flatten() + 0.5*(max(xs)+min(xs))
        Yb = 0.5*max_range*np.mgrid[-1:2:2, -1:2:2, -1:2:2][1].flatten() + 0.5*(max(ys)+min(ys))
        Zb = 0.5*max_range*np.mgrid[-1:2:2, -1:2:2, -1:2:2][2].flatten() + 0.5*(max(zs)+min(zs))
        # Comment or uncomment following both lines to test the fake bounding box:
        for xb, yb, zb in zip(Xb, Yb, Zb):
            ax.plot([xb], [yb], [zb], 'w')

        x1, x2 = ax.get_xlim()
        y1, y2 = ax.get_ylim()
        z1, z2 = ax.get_zlim()
#        ax.set_xlim3d(x1/zoom, x2/zoom)   # with this method the grid and the panes remain the same size within the subplot frame
#        ax.set_ylim3d(y1/zoom, y2/zoom)
#        ax.set_zlim3d(z1/zoom, z2/zoom)
        ax.get_proj = lambda: np.dot(Axes3D.get_proj(ax), np.diag([zoom, zoom, zoom, 1]))   # with this method the grid and the panes extend outside of the subplot frame

        if minmax:
            i = 0
            while sorted_c[i][3] <= min(vs2)+minmax_range:
                ax.scatter([sorted_c[i][0]], [sorted_c[i][1]], [sorted_c[i][2]], s=40+size*2, marker="v", color='blue')  # plot all points with minimum value
                i = i+1
            i = len(sorted_c)-1
            while sorted_c[i][3] >= max(vs2)-minmax_range:
                ax.scatter([sorted_c[i][0]], [sorted_c[i][1]], [sorted_c[i][2]], s=40+size*2, marker="^", color='red')  # plot all points with maximum value
                i = i-1

#        if curvature:
#            if cv_lim:
#                i = 0
#                while abs(cv_vdw[i][3]) <= cv_lim:
#                    
#                    ax.scatter([cv_vdw[i][0]], [cv_vdw[i][1]], [cv_vdw[i][2]], s=60+size*2, marker="1", color='black')
#                    i += 1
#            else:
#                for i in range(cv_pts):
#                    ax.scatter([cv_vdw[i][0]], [cv_vdw[i][1]], [cv_vdw[i][2]], s=60+size*2, marker="1", color='black')

        if axes == True:  # add axis labels
            ax.set_xlabel('x/'+l_units)
            ax.set_ylabel('y/'+l_units)
            ax.set_zlabel('z/'+l_units)
            ax.grid(True)
            ax.w_xaxis.set_pane_color((0.9, 0.9, 0.9, 0.5))
            ax.w_yaxis.set_pane_color((0.9, 0.9, 0.9, 0.5))
            ax.w_zaxis.set_pane_color((0.9, 0.9, 0.9, 0.5))

            if list(bkg) == [0, 0, 0] or type(bkg) == str and bkg.lower() == 'black':
                ax.spines['bottom'].set_color([0.8, 0.8, 0.8])
                ax.spines['top'].set_color([0.8, 0.8, 0.8])
                ax.spines['right'].set_color([0.8, 0.8, 0.8])
                ax.spines['left'].set_color([0.8, 0.8, 0.8])
                ax.tick_params(axis='both', colors=[0.8, 0.8, 0.8], which='both')
                ax.yaxis.label.set_color([0.8, 0.8, 0.8])
                ax.xaxis.label.set_color([0.8, 0.8, 0.8])
                ax.zaxis.label.set_color([0.8, 0.8, 0.8])
            else:
                ax.spines['bottom'].set_color([0, 0, 0])
                ax.spines['top'].set_color([0, 0, 0])
                ax.spines['right'].set_color([0, 0, 0])
                ax.spines['left'].set_color([0, 0, 0])
                ax.tick_params(axis='both', colors=[0, 0, 0])
                ax.yaxis.label.set_color([0, 0, 0])
                ax.xaxis.label.set_color([0, 0, 0])
                ax.zaxis.label.set_color([0, 0, 0])
        else:
            ax.w_xaxis.set_pane_color((1, 1, 1, 0))
            ax.w_yaxis.set_pane_color((1, 1, 1, 0))
            ax.w_zaxis.set_pane_color((1, 1, 1, 0))
            ax.grid(False)
            ax.set_xticks([])
            ax.set_yticks([])
            ax.set_zticks([])
            ax._axis3don = False
        
        if hl_x and hl_y and hl_z:
            ax.scatter([hl_x],[hl_y],[hl_z],marker='o',color='yellow',s=40+size*2)
            plt.gcf().text(0.8, 0.9, "  Highlighted:".center(15)+"\n"+('%.4E' % (gridpts[(hl_x,hl_y,hl_z)])).center(15),bbox = dict(boxstyle = 'round,pad=0.5', fc = (0.9,0.9,0.9)))     

            

    """
    Sliders and buttons

    for RadioButtons, the default value is the first one, so the order of the options is set so that the initial value is the first one
    """

    alpha_slider_ax = fig.add_axes([0.3, 0.09, 0.4, 0.03])
    alpha_slider = Slider(alpha_slider_ax, 'Alpha', 0, 1.0, valinit=alpha)

    size_slider_ax = fig.add_axes([0.3, 0.05, 0.4, 0.03])
    size_slider = Slider(size_slider_ax, 'Size', 0.1, 50.0, valinit=size)

    zoom_slider_ax = fig.add_axes([0.3, 0.01, 0.4, 0.03])
    zoom_slider = Slider(zoom_slider_ax, 'Zoom', 0.5, 5.0, valinit=zoom)
    
    elev_slider_ax = fig.add_axes([0.3, 0.13, 0.4, 0.03])
    elev_slider = Slider(elev_slider_ax, 'Elev', 0, 360, valinit=ax.elev)
    
    azim_slider_ax = fig.add_axes([0.3, 0.17, 0.4, 0.03])
    azim_slider = Slider(azim_slider_ax, 'Azim', 0, 360, valinit=ax.azim)
    
    iso_slider_ax = fig.add_axes([0.3, 0.21, 0.4, 0.03])
    iso_slider = Slider(iso_slider_ax, 'Value', min(v_vdw), max(v_vdw), valinit=iso)

    axes_ax = fig.add_axes([0.02, 0.83, 0.20, 0.13])
    if axes is False:
        axes_button = RadioButtons(axes_ax, ('axes off', 'axes on'))
    else:
        axes_button = RadioButtons(axes_ax, ('axes on', 'axes off'))

    minmax_ax = fig.add_axes([0.02, 0.69, 0.20, 0.13])
    if minmax:
        minmax_button = RadioButtons(minmax_ax, ('min/max on', 'min/max off'))
    else:
        minmax_button = RadioButtons(minmax_ax, ('min/max off', 'min/max on'))

    atoms_ax = fig.add_axes([0.02, 0.55, 0.20, 0.13])
    if atoms:
        atoms_button = RadioButtons(atoms_ax, ('atoms on', 'atoms off'))
    else:
        atoms_button = RadioButtons(atoms_ax, ('atoms off', 'atoms on'))

    # The colour button includes rainbow, bwr and the value passed by the user as a parameter (if different from rainbow and bwr)
    if colour != 'rainbow' and colour != 'bwr':
        col_b_size = 0.16  # if three options, make the height of the colour button larger
    else:
        col_b_size = 0.13  # if two options, make the height of the colour button smaller
    colour_ax = fig.add_axes([0.02, 0.54-col_b_size, 0.2, col_b_size])
    if colour != 'rainbow' and colour != 'bwr':
        colour_button = RadioButtons(colour_ax, (colour, 'rainbow', 'bwr'))
    elif colour == 'rainbow':
        colour_button = RadioButtons(colour_ax, ('rainbow', 'bwr'))
    else:
        colour_button = RadioButtons(colour_ax, ('bwr', 'rainbow'))

    # The background button includes White, Black and the value passed by the user as a parameter (if different from White or Black)
    # If the background colour is specified by the user in RGB, the option displayed on the button is the closest colour name as identified using the code from gistfile1. When changing back to this
    # option, the initial RGB colour is used, not the one defined by the closest name
    if list(bkg) != [1, 1, 1] and list(bkg) != [0, 0, 0] and not(type(bkg) == str and (bkg.lower() == 'white' or bkg.lower() == 'black')):
        bkg_b_size = 0.16  # if three options, make the height of the background button larger
    else:
        bkg_b_size = 0.13  # if two options, make the height of the background button smaller
    bkg_ax = fig.add_axes([0.02, 0.53-col_b_size-bkg_b_size, 0.2, bkg_b_size])
    if list(bkg) != [1, 1, 1] and list(bkg) != [0, 0, 0] and not(type(bkg) == str and (bkg.lower() == 'white' or bkg.lower() == 'black')):  # if different from white and black (both in RGB and colour name format)
        if type(bkg) == str:
            bkg_button = RadioButtons(bkg_ax, (bkg, 'White', 'Black'))
        else:
            bkg_button = RadioButtons(bkg_ax, (gf.ColorNames.findNearestWebColorName(bkg[0]*255, bkg[1]*255, bkg[2]*255), 'White', 'Black'))
    elif list(bkg) == [1, 1, 1] or type(bkg) == str and bkg.lower() == 'white':
        bkg_button = RadioButtons(bkg_ax, ('White', 'Black'))
    else:
        bkg_button = RadioButtons(bkg_ax, ('Black', 'White'))

    curv_ax = fig.add_axes([0.02, 0.39-col_b_size-bkg_b_size, 0.2, 0.13])
    if curvature:
        curv_button = RadioButtons(curv_ax, ('curvature on', 'curvature off'))
    else:
        curv_button = RadioButtons(curv_ax, ('curvature off', 'curvature on'))
        
    # plot with initial parameters
    ax.cla()
    ComputeAndPlot(alpha, size, axes, colour, minmax, bkg, zoom, atoms, curvature)
    
#    print("--- %s seconds ---" % (time.time() - start_time))
    
    # Functions defining what happens when the values of the sliders and buttons change
    def sliders_on_changed(val):
        nonlocal alpha, size, zoom
        alpha = alpha_slider.val
        size = size_slider.val
        zoom = zoom_slider.val

        ax.cla()
        ComputeAndPlot(alpha, size, axes, colour, minmax, bkg, zoom, atoms, curvature)
        plt.show()
        
    def iso_changed(val):
        nonlocal iso, x_p, y_p, z_p, v_p
        iso = iso_slider.val
        
        [x_p,y_p,z_p,v_p] = select_pts(iso,iso_range)
        
        ax.cla()
        ComputeAndPlot(alpha, size, axes, colour, minmax, bkg, zoom, atoms, curvature)
        plt.show()
        
    def angle_changed(val):
        elev = elev_slider.val
        azim = azim_slider.val

        ax.cla()
        ComputeAndPlot(alpha, size, axes, colour, minmax, bkg, zoom, atoms, curvature)
        ax.view_init(elev=elev, azim=azim)
        plt.show()
        
    def minmax_on_changed(val):
        nonlocal minmax
        b_dict = {'min/max on': True, 'min/max off': False}
        minmax = b_dict[val]

        ax.cla()
        ComputeAndPlot(alpha, size, axes, colour, minmax, bkg, zoom, atoms, curvature)
        plt.show()

    def axes_on_changed(val):
        nonlocal axes
        b_dict = {'axes on': True, 'axes off': False}
        axes = b_dict[val]

        ax.cla()
        ComputeAndPlot(alpha, size, axes, colour, minmax, bkg, zoom, atoms, curvature)
        plt.show()

    def atoms_on_changed(val):
        nonlocal atoms
        b_dict = {'atoms on': True, 'atoms off': False}
        atoms = b_dict[val]

        ax.cla()
        ComputeAndPlot(alpha, size, axes, colour, minmax, bkg, zoom, atoms, curvature)
        plt.show()

    def colour_on_changed(val):
        nonlocal colour
        colour = colour_button.value_selected

        ax.cla()
        ComputeAndPlot(alpha, size, axes, colour, minmax, bkg, zoom, atoms, curvature)
        plt.show()

    def bkg_on_changed(val):
        nonlocal bkg
        if bkg_button.value_selected == 'White':
            bkg = [1, 1, 1]
        elif bkg_button.value_selected == 'Black':
            bkg = [0, 0, 0]
        else:
            bkg = init_bkg

        ax.cla()
        ComputeAndPlot(alpha, size, axes, colour, minmax, bkg, zoom, atoms, curvature)
        plt.show()

    def curv_on_changed(val):
        nonlocal curvature
        b_dict = {'curvature on': True, 'curvature off': False}
        curvature = b_dict[val]

        ax.cla()
        ComputeAndPlot(alpha, size, axes, colour, minmax, bkg, zoom, atoms, curvature)
        plt.show()
        
    def onclick(event):
        nonlocal xdc2, ydc2, hl_x, hl_y, hl_z
        if event.dblclick:
            ax.button_pressed = None
            s = ax.format_coord(event.xdata,event.ydata)
            
            try:
                [xss, yss, zss] = s.split(",")
                x = float(xss.strip().split("=")[1])
                y = float(yss.strip().split("=")[1])
                z = float(zss.strip().split("=")[1])
                
                xdc2, ydc2, _ = proj3d.proj_transform(x,y,z, ax.get_proj())
                
                coords_2d = []
                for i in range(len(coords)):
                    x2, y2, _ = proj3d.proj_transform(coords[i][0],coords[i][1],coords[i][2], ax.get_proj())
                    coords_2d.append([coords[i][0],coords[i][1],coords[i][2],coords[i][3],x2,y2])

                xy_dc2 = np.asarray([xdc2,ydc2])  # x y double click 2d
                                                
                if xdc2 <= max(coords_2d,key=itemgetter(4))[4] and xdc2 >= min(coords_2d,key=itemgetter(4))[4] and ydc2 <= max(coords_2d,key=itemgetter(5))[5] and ydc2 >= min(coords_2d,key=itemgetter(5))[5]:
                    
                    for i in range(len(coords_2d)):
                        dist = np.sum((np.asarray([coords_2d[i][4],coords_2d[i][5]]) - xy_dc2)**2)
                        coords_2d[i].append(dist)

                    ordered_coords_2d = sorted(coords_2d,key=itemgetter(6))
                    
                    
                    hl_x = ordered_coords_2d[0][0]  # highlighted x 
                    hl_y = ordered_coords_2d[0][1]
                    hl_z = ordered_coords_2d[0][2]
                    ax.cla()
                    ComputeAndPlot(alpha, size, axes, colour, minmax, bkg, zoom, atoms, curvature)
                    plt.show()
            except:
                pass
            
        
    fig.canvas.mpl_connect('button_press_event', onclick)

    colour_button.on_clicked(colour_on_changed)
    axes_button.on_clicked(axes_on_changed)
    minmax_button.on_clicked(minmax_on_changed)
    alpha_slider.on_changed(sliders_on_changed)
    size_slider.on_changed(sliders_on_changed)
    zoom_slider.on_changed(sliders_on_changed)
    bkg_button.on_clicked(bkg_on_changed)
    atoms_button.on_clicked(atoms_on_changed)
    curv_button.on_clicked(curv_on_changed)
    elev_slider.on_changed(angle_changed)
    azim_slider.on_changed(angle_changed)
    iso_slider.on_changed(iso_changed)
    plt.show()


def PlotSlice(fpath, factor=1, alpha=1, size=4, axes=True, colour='rainbow', bkg=(1, 1, 1), minmax=False, minmax_range=0, cb_lim=None, atoms=False, au=False, curvature=False, cv_lim=None, cv_pts=5, zoom=1, a=0, b=0, c=1, x0=None, y0=None, z0=None, density=1, contour=True, value_type='esp', units=None):
    """
    arguments: path for a cube file (str), factor (float, opt), alpha (float, opt), size (float, opt), axes (bool, opt), colour (str, opt), background colour (tuple/list/str, opt), minmax (bool, opt), colour bar limits (list, opt), atoms (boolean, opt), au (boolean, opt), curvature (boolean, opt), curvature upper limit (float, opt), number of curvature points (int, opt)

    calls: GetVdWPoints (from ReadCube)

    returns: -

    this function creates an interactive plot of all the points which reside within +/- half a distance from the surface defined in terms of the Van der Waals radii of the atoms in the molecule
    (half a distance represents half the distance between two diagonally adjacent points in the grid)

    it has sliders for alpha, size and zoom

    it has buttons for turning axes on/off, showing/hiding min and max, showing/hiding atoms, colour scheme, background colour

    fpath (string) = path for a cube file 
    factor (float, default=1) = factor by which to multiply the VdW radii
    alpha (float, default=1) = transparency value (between 0 and 1)
    size (float, default=4) = marker size
    axes (boolean, default=True) = if True, show axes; if False, hide axes
    colour (string, default='rainbow') = colour scheme ('rainbow' by default, change it to 'bwr' for blue-white-red) - full list here
    bkg (tuple/list/string, default=(1, 1, 1)) = background colour (RGB tuple/list or name - each value between 0 and 1, white by default)
    minmax (boolean, default=False) = if True, show minima and maxima; if False, hide minima and maxima
    minmax_range (float, default=0) = = the value range to consider for the minima and maxima
    cb_lim (list, default=None) = limits for the colour bar (list of two floats; if not specified, the limits will be the min and max of the values; if the limits do not cover the whole range of values, a warning is printed)
    atoms (boolean, default=False) = if True, show atoms; if False, hide atoms
    au (boolean, default=False) = if True, use atomic units (i.e. bohr for coordinates); if False, use Angstroms for coordinates and units for values
    curvature (boolean, default=False) = if True, plots curvature points (either based on cv_lim if given, or cv_pts)
    cv_lim (float, default=None) = the upper limit of the curvature for the low curvature points to be plotted (it overrides cv_pts if given)
    cv_pts (integer, default=5) = the number of low curvature points to be plotted (the points will be chosen in increasing order of their absolute curvature
    zoom (float, default=1) = how large the plot is displayed with respect to the default size
    a (float, default=0) = the x component of the normal vector defining the plane
    b (float, default=0) = the y component of the normal vector defining the plane
    c (float, default=1) = the z component of the normal vector defining the plane
    x0 (float, default=None) = the x coordinate of the origin for the normal vector defining the plane
    y0 (float, default=None) = the y coordinate of the origin for the normal vector defining the plane
    z0 (float, default=None) = the z coordinate of the origin for the normal vector defining the plane
    density (integer, default=1) = an integer number n which indicates that the number of points stored along each axis is reduced by the density value n (therefore the total number of points is reduced by n3)
    contour (boolean, default=True) = if True, show van der Waals contour; if False, hide van der Waals contour
    value_type (string, default='esp') = can be ‘esp’/’potential’ or ‘dens’/’density’ (or any uppercase version of the aforementioned); used in conjunction with the keywords au and units in order to convert the values to
                                        the appropriate units
    units (string, default=None) = the units to be used for the value; currently supports ‘V’ for ESP and ‘A’/‘angstrom’, ’nm’, ’pm’, ’m’ for density (or any lowercase/uppercase version of these). In the case of the density,
                                    the actual units are the specified units ^(-3). The keywords value_type and au override units. If au is True, the values will be in atomic units regardless of the value passed for units. 
                                If value_type is ‘dens’, the values will be read as density values even if the units argument has a valid ESP value (such as ‘V’). The default unit for ESP is V and for density 1/Å3.
    """
    
#    start_time = time.time()
    
    print('\n\n   “Patience is the companion of wisdom.”\n'+'(St. Augustine)'.rjust(41))
    
    #########################################
    # Extracting data in convenient formats #
    #########################################
    
    exdata = rc.ExtractData(fpath, au=au, value_type=value_type, density=density, units=units)
    origin=exdata[1]
    n_x = exdata[3]
    x_vector=exdata[4]
    n_y=exdata[5]
    y_vector=exdata[6]
    n_z=exdata[7]
    z_vector=exdata[8]
    atoms_list = exdata[9]
    vals=exdata[12]
    val_units=exdata[13]
    axes_list=rc.Axes(origin=origin, n_x = n_x, x_vector=x_vector, n_y=n_y, y_vector=y_vector,n_z=n_z, z_vector=z_vector)
    (gridpts,d_ind)=rc.ValuesAsDictionary(vals=vals,axes_list=axes_list, origin=origin)
    
    vdw_coords = rc.GetVdWPoints(axes_list=axes_list,gridpts=gridpts,atoms_list=atoms_list, x_vector=x_vector, y_vector=y_vector, z_vector=z_vector, origin=origin, factor=factor, au=au)[0]
    
    coords = [list(i)+[gridpts[i]] for i in list(gridpts.keys())]
    
    dist = ((x_vector[0]+y_vector[0]+z_vector[0])**2+(x_vector[1]+y_vector[1]+z_vector[1])**2+(x_vector[2]+y_vector[2]+z_vector[2])**2)**0.5
    half_dist = 0.3*dist
    
    
    x_vdw, y_vdw, z_vdw, v_vdw = [list(j) for j in zip(*vdw_coords)]
    
    sorted_c = sorted(coords, key=itemgetter(3))  # sorted list of coordinates by value
    
    xs, ys, zs, vs = [list(j) for j in zip(*coords)]
    vs2 = np.asarray(vs)
    
    ##########################
    # Interpreting arguments #
    ##########################
    
    if x0 is None:
        x0 = (min(xs)+max(xs))/2
    if y0 is None:
        y0 = (min(ys)+max(ys))/2
    if z0 is None:
        z0 = (min(zs)+max(zs))/2
        
    if au:
        l_units = '$a_0$'
    else:
        l_units = r'$\AA$'

    if cb_lim is None:
        cb_lim = [min(v_vdw), max(v_vdw)]
    elif cb_lim[0] > min(vs2) or cb_lim[1] < max(vs2):
        print('\nWarning: The colour bar limits you have selected do not cover the whole range of values\n')

    ###################
    # Creating figure #
    ###################

    fig = plt.figure()
    ax = fig.add_subplot(111, projection='3d')  # 1x1 grid, 1st subplot, 3D
    ax.set_aspect('equal')   # equal aspect ratio

    fig.subplots_adjust(bottom=0.3)  # leave space for the slider

    cb_ax = fig.add_axes([0.85, 0.25, 0.05, 0.6])  # axes for colour bar
    
    plt.gcf().text(0.02, 0.97, fpath.split('/')[-1], fontsize=7)     

    ########################
    # Additional variables #
    ########################

    zoom = 1
    
    xdc2 = None
    ydc2 = None
    hl_x = None
    hl_y = None
    hl_z = None

    u = np.linspace(0, 2 * np.pi, 100)
    v = np.linspace(0, np.pi, 100)

    init_bkg = bkg  # copy of the initial background value

    
    def select_pts(a,b,c,x0,y0,z0): 
        x_p = []
        y_p = []
        z_p = []
        v_p = []
        
        x_out = []
        y_out = []
        z_out = []
        
        d = -(a*x0+b*y0+c*z0)
        
        for i in range(len(coords)):
            d_to_plane = abs(a*coords[i][0]+b*coords[i][1]+c*coords[i][2]+d)/(a**2+b**2+c**2)**0.5
            if d_to_plane < half_dist:
                x_p.append(coords[i][0])
                y_p.append(coords[i][1])
                z_p.append(coords[i][2])
                v_p.append(coords[i][3])
                
        for i in range(len(x_vdw)):
            d_to_plane = abs(a*x_vdw[i]+b*y_vdw[i]+c*z_vdw[i]+d)/(a**2+b**2+c**2)**0.5
            if d_to_plane < half_dist:
                x_out.append(x_vdw[i])
                y_out.append(y_vdw[i])
                z_out.append(z_vdw[i])                        

        return [x_p,y_p,z_p,v_p],[x_out,y_out,z_out]

    [[x_p,y_p,z_p,v_p], [x_out,y_out,z_out]]=select_pts(a,b,c,x0,y0,z0)
    
    
    def ComputeAndPlot(alpha, size, axes, colour, minmax, bkg, zoom, atoms, curvature):
        ax.set_facecolor(bkg)   # set background colour for 3D plot
        # plot atoms if True
        if atoms:
            for at in atoms_list:
                x_origin = at[2]
                y_origin = at[3]
                z_origin = at[4]

                if au:
                    rf = 0.529177  # lengths in atomic units (i.e. bohr)
                else:
                    rf = 1  # lengths in angstroms

                r = 0.4*rc.vdwdict[at[0]]/rf

                x_box = x_origin + r * np.outer(np.cos(u), np.sin(v))
                y_box = y_origin + r * np.outer(np.sin(u), np.sin(v))
                z_box = z_origin + r * np.outer(np.ones(np.size(u)), np.cos(v))

                ax.plot_surface(x_box, y_box, z_box,  rstride=4, cstride=4, color=coldict[at[0]], linewidth=0, alpha=1)


        # create colour bar
        cmap = getattr(cm, colour)
        norm = matplotlib.colors.Normalize(vmin=cb_lim[0], vmax=cb_lim[1])
        matplotlib.colorbar.ColorbarBase(cb_ax, cmap=cmap, norm=norm, label=val_units)

        
        # plot points
        ax.scatter(x_p, y_p, z_p, c=v_p, cmap=cmap, norm=norm, s=size, alpha=alpha)
        
        # plot contour
        if contour:
            ax.plot(x_out,y_out,z_out,linestyle='',marker='o',markersize=1,color=(0,0,0))
        
        # Create cubic bounding box to simulate equal aspect ratio
        max_range = np.array([max(xs)-min(xs), max(ys)-min(ys), max(zs)-min(zs)]).max()
        Xb = 0.5*max_range*np.mgrid[-1:2:2, -1:2:2, -1:2:2][0].flatten() + 0.5*(max(xs)+min(xs))
        Yb = 0.5*max_range*np.mgrid[-1:2:2, -1:2:2, -1:2:2][1].flatten() + 0.5*(max(ys)+min(ys))
        Zb = 0.5*max_range*np.mgrid[-1:2:2, -1:2:2, -1:2:2][2].flatten() + 0.5*(max(zs)+min(zs))
        # Comment or uncomment following both lines to test the fake bounding box:
        for xb, yb, zb in zip(Xb, Yb, Zb):
            ax.plot([xb], [yb], [zb], 'w')

        x1, x2 = ax.get_xlim()
        y1, y2 = ax.get_ylim()
        z1, z2 = ax.get_zlim()
#        ax.set_xlim3d(x1/zoom, x2/zoom)   # with this method the grid and the panes remain the same size within the subplot frame
#        ax.set_ylim3d(y1/zoom, y2/zoom)
#        ax.set_zlim3d(z1/zoom, z2/zoom)
        ax.get_proj = lambda: np.dot(Axes3D.get_proj(ax), np.diag([zoom, zoom, zoom, 1]))   # with this method the grid and the panes extend outside of the subplot frame

        if minmax:
            i = 0
            while sorted_c[i][3] <= min(vs2)+minmax_range:
                ax.scatter([sorted_c[i][0]], [sorted_c[i][1]], [sorted_c[i][2]], s=40+size*2, marker="v", color='blue')  # plot all points with minimum value
                i = i+1
            i = len(sorted_c)-1
            while sorted_c[i][3] >= max(vs2)-minmax_range:
                ax.scatter([sorted_c[i][0]], [sorted_c[i][1]], [sorted_c[i][2]], s=40+size*2, marker="^", color='red')  # plot all points with maximum value
                i = i-1

#        if curvature:
#            if cv_lim:
#                i = 0
#                while abs(cv_vdw[i][3]) <= cv_lim:
#                    
#                    ax.scatter([cv_vdw[i][0]], [cv_vdw[i][1]], [cv_vdw[i][2]], s=60+size*2, marker="1", color='black')
#                    i += 1
#            else:
#                for i in range(cv_pts):
#                    ax.scatter([cv_vdw[i][0]], [cv_vdw[i][1]], [cv_vdw[i][2]], s=60+size*2, marker="1", color='black')

        if axes == True:  # add axis labels
            ax.set_xlabel('x/'+l_units)
            ax.set_ylabel('y/'+l_units)
            ax.set_zlabel('z/'+l_units)
            ax.grid(True)
            ax.w_xaxis.set_pane_color((0.9, 0.9, 0.9, 0.5))
            ax.w_yaxis.set_pane_color((0.9, 0.9, 0.9, 0.5))
            ax.w_zaxis.set_pane_color((0.9, 0.9, 0.9, 0.5))

            if list(bkg) == [0, 0, 0] or type(bkg) == str and bkg.lower() == 'black':
                ax.spines['bottom'].set_color([0.8, 0.8, 0.8])
                ax.spines['top'].set_color([0.8, 0.8, 0.8])
                ax.spines['right'].set_color([0.8, 0.8, 0.8])
                ax.spines['left'].set_color([0.8, 0.8, 0.8])
                ax.tick_params(axis='both', colors=[0.8, 0.8, 0.8], which='both')
                ax.yaxis.label.set_color([0.8, 0.8, 0.8])
                ax.xaxis.label.set_color([0.8, 0.8, 0.8])
                ax.zaxis.label.set_color([0.8, 0.8, 0.8])
            else:
                ax.spines['bottom'].set_color([0, 0, 0])
                ax.spines['top'].set_color([0, 0, 0])
                ax.spines['right'].set_color([0, 0, 0])
                ax.spines['left'].set_color([0, 0, 0])
                ax.tick_params(axis='both', colors=[0, 0, 0])
                ax.yaxis.label.set_color([0, 0, 0])
                ax.xaxis.label.set_color([0, 0, 0])
                ax.zaxis.label.set_color([0, 0, 0])
        else:
            ax.w_xaxis.set_pane_color((1, 1, 1, 0))
            ax.w_yaxis.set_pane_color((1, 1, 1, 0))
            ax.w_zaxis.set_pane_color((1, 1, 1, 0))
            ax.grid(False)
            ax.set_xticks([])
            ax.set_yticks([])
            ax.set_zticks([])
            ax._axis3don = False
        
        if hl_x and hl_y and hl_z:
            ax.scatter([hl_x],[hl_y],[hl_z],marker='o',color='yellow',s=40+size*2)
            plt.gcf().text(0.8, 0.9, "  Highlighted:".center(15)+"\n"+('%.4E' % (gridpts[(hl_x,hl_y,hl_z)])).center(15),bbox = dict(boxstyle = 'round,pad=0.5', fc = (0.9,0.9,0.9)))     
            

    """
    Sliders and buttons

    for RadioButtons, the default value is the first one, so the order of the options is set so that the initial value is the first one
    """
    

    zoom_slider_ax = fig.add_axes([0.3, 0.01, 0.4, 0.03])
    zoom_slider = Slider(zoom_slider_ax, 'Zoom', 0.5, 5.0, valinit=zoom)
    
    a_slider_ax = fig.add_axes([0.3,0.16,0.15,0.03])
    a_slider = Slider(a_slider_ax,'a',0,1,valinit=a)

    b_slider_ax = fig.add_axes([0.3,0.12,0.15,0.03])
    b_slider = Slider(b_slider_ax,'b',0,1,valinit=b)

    c_slider_ax = fig.add_axes([0.3,0.08,0.15,0.03])
    c_slider = Slider(c_slider_ax,'c',0,1,valinit=c)
    
    ox_slider_ax = fig.add_axes([0.55,0.16,0.15,0.03])
    ox_slider = Slider(ox_slider_ax,'$x_0$',min(xs),max(xs),valinit=x0)
    
    oy_slider_ax = fig.add_axes([0.55,0.12,0.15,0.03])
    oy_slider = Slider(oy_slider_ax,'$y_0$',min(ys),max(ys),valinit=y0)

    oz_slider_ax = fig.add_axes([0.55,0.08,0.15,0.03])
    oz_slider = Slider(oz_slider_ax,'$z_0$',min(zs),max(zs),valinit=z0)
    
    axes_ax = fig.add_axes([0.02, 0.83, 0.20, 0.13])
    if axes is False:
        axes_button = RadioButtons(axes_ax, ('axes off', 'axes on'))
    else:
        axes_button = RadioButtons(axes_ax, ('axes on', 'axes off'))

    minmax_ax = fig.add_axes([0.02, 0.69, 0.20, 0.13])
    if minmax:
        minmax_button = RadioButtons(minmax_ax, ('min/max on', 'min/max off'))
    else:
        minmax_button = RadioButtons(minmax_ax, ('min/max off', 'min/max on'))

    atoms_ax = fig.add_axes([0.02, 0.55, 0.20, 0.13])
    if atoms:
        atoms_button = RadioButtons(atoms_ax, ('atoms on', 'atoms off'))
    else:
        atoms_button = RadioButtons(atoms_ax, ('atoms off', 'atoms on'))

    # The colour button includes rainbow, bwr and the value passed by the user as a parameter (if different from rainbow and bwr)
    if colour != 'rainbow' and colour != 'bwr':
        col_b_size = 0.16  # if three options, make the height of the colour button larger
    else:
        col_b_size = 0.13  # if two options, make the height of the colour button smaller
    colour_ax = fig.add_axes([0.02, 0.54-col_b_size, 0.2, col_b_size])
    if colour != 'rainbow' and colour != 'bwr':
        colour_button = RadioButtons(colour_ax, (colour, 'rainbow', 'bwr'))
    elif colour == 'rainbow':
        colour_button = RadioButtons(colour_ax, ('rainbow', 'bwr'))
    else:
        colour_button = RadioButtons(colour_ax, ('bwr', 'rainbow'))

    # The background button includes White, Black and the value passed by the user as a parameter (if different from White or Black)
    # If the background colour is specified by the user in RGB, the option displayed on the button is the closest colour name as identified using the code from gistfile1. When changing back to this
    # option, the initial RGB colour is used, not the one defined by the closest name
    if list(bkg) != [1, 1, 1] and list(bkg) != [0, 0, 0] and not(type(bkg) == str and (bkg.lower() == 'white' or bkg.lower() == 'black')):
        bkg_b_size = 0.16  # if three options, make the height of the background button larger
    else:
        bkg_b_size = 0.13  # if two options, make the height of the background button smaller
    bkg_ax = fig.add_axes([0.02, 0.53-col_b_size-bkg_b_size, 0.2, bkg_b_size])
    if list(bkg) != [1, 1, 1] and list(bkg) != [0, 0, 0] and not(type(bkg) == str and (bkg.lower() == 'white' or bkg.lower() == 'black')):  # if different from white and black (both in RGB and colour name format)
        if type(bkg) == str:
            bkg_button = RadioButtons(bkg_ax, (bkg, 'White', 'Black'))
        else:
            bkg_button = RadioButtons(bkg_ax, (gf.ColorNames.findNearestWebColorName(bkg[0]*255, bkg[1]*255, bkg[2]*255), 'White', 'Black'))
    elif list(bkg) == [1, 1, 1] or type(bkg) == str and bkg.lower() == 'white':
        bkg_button = RadioButtons(bkg_ax, ('White', 'Black'))
    else:
        bkg_button = RadioButtons(bkg_ax, ('Black', 'White'))

    curv_ax = fig.add_axes([0.02, 0.39-col_b_size-bkg_b_size, 0.2, 0.13])
    if curvature:
        curv_button = RadioButtons(curv_ax, ('curvature on', 'curvature off'))
    else:
        curv_button = RadioButtons(curv_ax, ('curvature off', 'curvature on'))
        
    # plot with initial parameters
    ax.cla()
    ComputeAndPlot(alpha, size, axes, colour, minmax, bkg, zoom, atoms, curvature)

    # Functions defining what happens when the values of the sliders and buttons change
    def sliders_on_changed(val):
        nonlocal zoom
        zoom = zoom_slider.val

        ax.cla()
        ComputeAndPlot(alpha, size, axes, colour, minmax, bkg, zoom, atoms, curvature)
        plt.show()
        
    def plane_sliders_on_changed(val):
        nonlocal x_p, y_p, z_p, v_p, a, b, c, x0, y0, z0, x_out, y_out, z_out
        a = a_slider.val
        b = b_slider.val
        c = c_slider.val
        x0 = ox_slider.val
        y0 = oy_slider.val
        z0 = oz_slider.val
        
        [[x_p,y_p,z_p,v_p],[x_out,y_out,z_out]] = select_pts(a,b,c,x0,y0,z0)
        
        ax.cla()
        ComputeAndPlot(alpha, size, axes, colour, minmax, bkg, zoom, atoms, curvature)
        plt.show()
        
    def minmax_on_changed(val):
        nonlocal minmax
        b_dict = {'min/max on': True, 'min/max off': False}
        minmax = b_dict[val]

        ax.cla()
        ComputeAndPlot(alpha, size, axes, colour, minmax, bkg, zoom, atoms, curvature)
        plt.show()

    def axes_on_changed(val):
        nonlocal axes
        b_dict = {'axes on': True, 'axes off': False}
        axes = b_dict[val]

        ax.cla()
        ComputeAndPlot(alpha, size, axes, colour, minmax, bkg, zoom, atoms, curvature)
        plt.show()

    def atoms_on_changed(val):
        nonlocal atoms
        b_dict = {'atoms on': True, 'atoms off': False}
        atoms = b_dict[val]

        ax.cla()
        ComputeAndPlot(alpha, size, axes, colour, minmax, bkg, zoom, atoms, curvature)
        plt.show()

    def colour_on_changed(val):
        nonlocal colour
        colour = colour_button.value_selected

        ax.cla()
        ComputeAndPlot(alpha, size, axes, colour, minmax, bkg, zoom, atoms, curvature)
        plt.show()

    def bkg_on_changed(val):
        nonlocal bkg
        if bkg_button.value_selected == 'White':
            bkg = [1, 1, 1]
        elif bkg_button.value_selected == 'Black':
            bkg = [0, 0, 0]
        else:
            bkg = init_bkg

        ax.cla()
        ComputeAndPlot(alpha, size, axes, colour, minmax, bkg, zoom, atoms, curvature)
        plt.show()

    def curv_on_changed(val):
        nonlocal curvature
        b_dict = {'curvature on': True, 'curvature off': False}
        curvature = b_dict[val]

        ax.cla()
        ComputeAndPlot(alpha, size, axes, colour, minmax, bkg, zoom, atoms, curvature)
        plt.show()
             
    def onclick(event):
        nonlocal xdc2, ydc2, hl_x, hl_y, hl_z
        if event.dblclick:
            ax.button_pressed = None
            s = ax.format_coord(event.xdata,event.ydata)
            
            try:
                [xss, yss, zss] = s.split(",")
                x = float(xss.strip().split("=")[1])
                y = float(yss.strip().split("=")[1])
                z = float(zss.strip().split("=")[1])
                
                xdc2, ydc2, _ = proj3d.proj_transform(x,y,z, ax.get_proj())
                
                coords_2d = []
                for i in range(len(coords)):
                    x2, y2, _ = proj3d.proj_transform(coords[i][0],coords[i][1],coords[i][2], ax.get_proj())
                    coords_2d.append([coords[i][0],coords[i][1],coords[i][2],coords[i][3],x2,y2])

                xy_dc2 = np.asarray([xdc2,ydc2])  # x y double click 2d
                                                
                if xdc2 <= max(coords_2d,key=itemgetter(4))[4] and xdc2 >= min(coords_2d,key=itemgetter(4))[4] and ydc2 <= max(coords_2d,key=itemgetter(5))[5] and ydc2 >= min(coords_2d,key=itemgetter(5))[5]:
                    
                    for i in range(len(coords_2d)):
                        dist = np.sum((np.asarray([coords_2d[i][4],coords_2d[i][5]]) - xy_dc2)**2)
                        coords_2d[i].append(dist)

                    ordered_coords_2d = sorted(coords_2d,key=itemgetter(6))
                    
                    
                    hl_x = ordered_coords_2d[0][0]  # highlighted x 
                    hl_y = ordered_coords_2d[0][1]
                    hl_z = ordered_coords_2d[0][2]
                    ax.cla()
                    ComputeAndPlot(alpha, size, axes, colour, minmax, bkg, zoom, atoms, curvature)
                    plt.show()
            except:
                pass
            
    fig.canvas.mpl_connect('button_press_event', onclick)
    zoom_slider.on_changed(sliders_on_changed)
    colour_button.on_clicked(colour_on_changed)
    axes_button.on_clicked(axes_on_changed)
    minmax_button.on_clicked(minmax_on_changed)
    bkg_button.on_clicked(bkg_on_changed)
    atoms_button.on_clicked(atoms_on_changed)
    curv_button.on_clicked(curv_on_changed)
    
    a_slider.on_changed(plane_sliders_on_changed)
    b_slider.on_changed(plane_sliders_on_changed)
    c_slider.on_changed(plane_sliders_on_changed)
    ox_slider.on_changed(plane_sliders_on_changed)
    oy_slider.on_changed(plane_sliders_on_changed)
    oz_slider.on_changed(plane_sliders_on_changed)

    plt.show()


def PlotHistogram(fpath, factor=1, bins=30, hist=True, kde=True, norm_freq=False, norm_hist=True, au=False, density=1, value_type='esp', units=None, colour='rainbow', cb_lim=None, xlim=None, ylim=None, save=True, cat=False, an=False, scaling_factor=1):
    """
    arguments: path for a cube file (str), factor (float, opt), number of bins (None/int, opt), histogram (bool, opt), kde (bool, opt), normalise frequency (bool, opt),
               normalise histogram (boolean, opt)

    calls: GetVdWPoints (from ReadCube)

    returns: -

    this function creates a non-interactive histogram plot of the values of all the points which reside within +/- half a distance from the surface defined in terms of the Van der Waals radii
    of the atoms in the molecule along with a smooth histogram
    (half a distance represents half the distance between two diagonally adjacent points in the grid)

    fpath (string) = path for a cube file
    factor (float, default=1) = factor by which to multiply the VdW radii
    bins (integer/None/string, default=30) = the number of bins in the histogram (see here for details)
    hist (boolean, default=True) = if True, it plots the histogram; if False, it does not plot the histogram (only the KDE)
    kde (boolean, default=True) = if True, it plots a gaussian kernel density estimate (KDE = a way to estimate the probability density function of a variable); if False, it does not plot the KDE
    norm_freq (boolean, default=False) = if True, it plots a histogram with normalised frequency (without KDE, because the KDE has a normalised area, not normalised bin heights)
    norm_hist (boolean, default=True) = if True, it normalises the histogram (i.e. integral=1) - this happens by default for plots including KDE
    au (boolean, default=False) = if True, use atomic units (i.e. bohr for coordinates); if False, use Angstroms for coordinates and units for values
    density (integer, default=1) = an integer number n which indicates that the number of points stored along each axis is reduced by the density value n (therefore the total number of points is reduced by n3)
    value_type (string, default='esp') = can be ‘esp’/’potential’ or ‘dens’/’density’ (or any uppercase version of the aforementioned); used in conjunction with the keywords au and units in order to convert the values to the appropriate units
    units (string, default=None) = the units to be used for the value; currently supports ‘V’ for ESP and ‘A’/‘angstrom’, ’nm’, ’pm’, ’m’ for density (or any lowercase/uppercase version of these). In the case of the density, 
                                    the actual units are the specified units ^(-3). The keywords value_type and au override units. If au is True, the values will be in atomic units regardless of the value passed for units. 
                                    If value_type is ‘dens’, the values will be read as density values even if the units argument has a valid ESP value (such as ‘V’). The default unit for ESP is V and for density 1/Å3.
    colour (string, default='rainbow') = colour scheme ('rainbow' by default, change it to 'bwr' for blue-white-red) - full list here
    cb_lim (list, default=None) = limits for the colour bar (list of two floats; if not specified, the limits will be the min and max of the values; if the limits do not cover the whole range of values, a warning is printed)
    xlim (list, default=None) = x axis limits
    ylim (list, default=None) = y axis limits
    save (boolean, default=True) = if True: if hist is also True, a file <original_file_name>_hist.out is created, containing the bin edges values and the frequency for each bin as shown in the histogram; if kde is also True, 
                                    a file <original_file_name>_kde.out is created, containing the range of values used for plotting the KDE and the corresponding KDE value. The file ends with two separate sections: 
                                    peak values and minimum values
    cat (boolean, default=False) = if True, it only plots the histogram and/or KDE for the points that correspond to the cation (the condition for which atoms belong to the cation has to be changed manually inside GetVdWPoints)
    an (boolean, default=False) = if True, it only plots the histogram and/or KDE for the points that correspond to the anion (the condition for which atoms belong to the cation has to be changed manually inside GetVdWPoints; 
                                   the atoms belonging to the anion are the remaining ones)    
    
    """
    
    #########################################
    # Extracting data in convenient formats #
    #########################################
    
    exdata = rc.ExtractData(fpath, au=au, value_type=value_type, density=density, units=units)
    origin=exdata[1]
    n_x = exdata[3]
    x_vector=exdata[4]
    n_y=exdata[5]
    y_vector=exdata[6]
    n_z=exdata[7]
    z_vector=exdata[8]
    atoms_list=exdata[9]
    vals=exdata[12]
    val_units=exdata[13]
    axes_list=rc.Axes(origin=origin, n_x = n_x, x_vector=x_vector, n_y=n_y, y_vector=y_vector,n_z=n_z, z_vector=z_vector)
    gridpts=rc.ValuesAsDictionary(vals=vals,axes_list=axes_list, origin=origin)[0]
    coords = rc.GetVdWPoints(axes_list=axes_list,gridpts=gridpts,atoms_list=atoms_list, x_vector=x_vector, y_vector=y_vector, z_vector=z_vector, origin=origin, factor=factor, au=au, cat=cat, an=an)[0]
    xs, ys, zs, vs = [list(x) for x in zip(*coords)]

    ##########################
    # Interpreting arguments #
    ##########################

    if value_type.lower()=='esp' or value_type.lower()=='potential':
        title = 'Electrostatic potential /'+val_units
    elif value_type.lower()=='dens' or value_type.lower()=='density':
        title = 'Density /'+val_units
    if cb_lim is None:
        cb_lim = [min(vs), max(vs)]
    elif cb_lim[0] > min(vs) or cb_lim[1] < max(vs):
        print('\nWarning: The colour bar limits you have selected do not cover the whole range of values\n')
        
    if norm_freq:   #norm_freq overrides everything else
        hist = True
        kde = False
        norm_hist = False
    elif kde:
        norm_hist = True  #kde overrides norm_hist 
    else:
        hist=True  # kde overrides hist
    
        
    if xlim is None:
        xlim=[cb_lim[0]-(max(vs)-min(vs))/10,cb_lim[1]+(max(vs)-min(vs))/10]
    xmin=xlim[0]
    xmax=xlim[1]
    
    if ylim is None:
        ylim=[0,None]
    ymin=ylim[0]
    ymax=ylim[1]
        
    ###################
    # Creating figure #
    ###################

    fig = plt.figure()
    plt.gcf().text(0.02, 0.97, fpath.split('/')[-1], fontsize=7)
    ax = fig.add_subplot(111)  # 1x1 grid, 1st subplot, 3D
    ax.set_facecolor((0.92,0.92,0.95))
    plt.grid(color=(1,1,1))
    ax.spines['bottom'].set_color([1, 1, 1])
    ax.spines['top'].set_color([1, 1, 1])
    ax.spines['right'].set_color([1, 1, 1])
    ax.spines['left'].set_color([1, 1, 1])

# This is  the colormap I'd like to use.

    cmap = getattr(cm, colour)
    norm = matplotlib.colors.Normalize(vmin=cb_lim[0], vmax=cb_lim[1])

    # Plot histogram
    if hist:
        
        if norm_freq:
            weights = np.ones_like(vs)/float(len(vs))
            n, b, patches = ax.hist(vs, bins=bins, weights=weights)
        else:
            n, b, patches = ax.hist(vs, bins=bins, normed=norm_hist)
        
        bin_centers = 0.5 * (b[:-1] + b[1:])
                
        
        # scale values to interval [0,1]
        col = bin_centers - min(bin_centers)
        col /= max(col)
    
        for c, p in zip(bin_centers, patches):
            plt.setp(p, 'facecolor', cmap(norm(c)))
    
    if kde:
        kde = stats.gaussian_kde(vs)
        xx = np.append(np.linspace(xmin,min(vs)-(max(vs)-min(vs))/2,100),np.linspace(min(vs)-(max(vs)-min(vs))/2,max(vs)+(max(vs)-min(vs))/2,1000),0)  # higher density of points in the relevant region
        xx = np.append(xx,np.linspace(max(vs)+(max(vs)-min(vs))/2,xmax,100))
        ax.plot(xx,kde(xx),color=(0.1,0.1,0.1))
                
            
    if title is None:
        ax.set_xlabel('Value')
        plt.title('Histogram')
    else:
        ax.set_xlabel(title)    
        plt.title(title.split("/")[0]+'histogram')
    ax.set_ylabel('Frequency')
    
    ax.set_axisbelow(True)
    
    if xmin is not None:
        ax.set_xlim(xmin=xmin)
    if xmax is not None:
        ax.set_xlim(xmax=xmax)
    if ymin is not None:
        ax.set_ylim(ymin=ymin)
    if ymax is not None:
        ax.set_ylim(ymax=ymax)

    if hist and save:
        if factor==1:
            if cat:
                outhist=fpath[:-5]+'_cat_hist.out'
            elif an:
                outhist=fpath[:-5]+'_an_hist.out'
            else:
                outhist=fpath[:-5]+'_hist.out'
        else:
            if cat:
                outhist=fpath[:-5]+'_f'+str(factor)+'_cat_hist.out'
            elif an:
                outhist=fpath[:-5]+'_f'+str(factor)+'_an_hist.out'
            else:
                outhist=fpath[:-5]+'_f'+str(factor)+'_hist.out'
        
        with open(outhist,'w') as oh:
            oh.write("### Histogram data ###\n")
            oh.write(("Bin edges /"+val_units).center(30)+"\t\t"+"Frequency".center(17)+"\n")
            for i in range(len(n)):    
                oh.write(str(b[i])+"\t"+str(b[i+1])+"\t\t"+str(n[i])+"\n")
                
    if kde and save:
        if factor==1:
            if cat:
                outkde=fpath[:-5]+'_cat_kde.out'
            elif an:
                outkde=fpath[:-5]+'_an_kde.out'
            else:
                outkde=fpath[:-5]+'_kde.out'
        else:
            if cat:
                outkde=fpath[:-5]+'_f'+str(factor)+'_cat_kde.out'
            elif an:
                outkde=fpath[:-5]+'_f'+str(factor)+'_an_kde.out'
            else:
                outkde=fpath[:-5]+'_f'+str(factor)+'_kde.out'
        
        with open(outkde,'w') as ok:
            ok.write("### KDE data ###\n")
            ok.write(("Value /"+val_units).center(13)+"\t\t"+"Frequency".center(17)+"\n")
            peakvals=[]
            peakfreqs=[]
            minvals=[]
            minfreqs=[]
            for i in range(len(xx)):
                ok.write(str(xx[i])+"\t\t"+str(kde(xx[i])[0]*scaling_factor)+"\n")
                if i>0 and i<len(xx)-1 and kde(xx[i])[0]>kde(xx[i-1])[0] and kde(xx[i])[0]>kde(xx[i+1])[0]:
                    peakvals.append(xx[i])
                    peakfreqs.append(kde(xx[i])[0]*scaling_factor)
                    
                if i>0 and i<len(xx)-1 and kde(xx[i])[0]<kde(xx[i-1])[0] and kde(xx[i])[0]<kde(xx[i+1])[0]:
                    minvals.append(xx[i])
                    minfreqs.append(kde(xx[i])[0]*scaling_factor)
                    
            ok.write("\n### PEAK VALUES ###\n")
            for i in range(len(peakvals)):
                ok.write(str(peakvals[i])+"\t\t"+str(peakfreqs[i])+"\n")
                
            ok.write("\n### MIN VALUES ###\n")
            for i in range(len(minvals)):
                ok.write(str(minvals[i])+"\t\t"+str(minfreqs[i])+"\n")
    plt.show()
            

if __name__ == '__main__':
    # Map command line arguments to function arguments.
    fct = str(sys.argv[1])

    comp_args = []    # compulsory arguments
    arg_dict = {}     # optional arguments
    for argument in sys.argv[2:]:
        if "=" not in str(argument):
            comp_args.append(argument)
        else:
            try:
                arg_dict[str(argument).split('=')[0]] = eval(str(argument).split('=')[1])   # eval turns a string into the right format (list, float, boolean etc.)
            except NameError:
                arg_dict[str(argument).split('=')[0]] = str(argument).split('=')[1]

    globals()[fct](*comp_args, **arg_dict)
