# -*- coding: utf-8 -*-
"""
A collection of functions that extract 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 sys
import collections as cd
import time
from operator import itemgetter


# a dictionary of van der Waals radii (most from http://periodictable.com/Properties/A/VanDerWaalsRadius.an.html)
vdwdict2 = { 1:1.20,  2:1.40,  3:1.82,  4:1.53,  5:1.80,  6:1.70,  7:1.55,  8:1.52,  9:1.47, 10:1.54,
           11:2.27, 12:1.73, 13:1.84, 14:2.10, 15:1.80, 16:1.80, 17:1.75, 18:1.88, 19:2.75, 20:2.31,
           21:2.11,                                                       28:1.63, 29:1.40, 30:1.39,
           31:1.87, 32:2.11, 33:1.85, 34:1.90, 35:1.85, 36:2.02, 37:3.03, 38:2.49, 39:2.19,
                                                        46:1.63, 47:1.72, 48:1.58, 49:1.93, 50:2.17,
           51:2.06, 52:2.06, 53:1.98, 54:2.16, 55:3.43, 56:2.58,
           
                                                                          78:1.75, 79:1.66, 80:1.55,
           81:1.96, 82:2.02, 83:2.07, 84:1.97, 85:2.02, 86:2.20, 87:3.48, 88:2.83, 
                    92:1.86,          94:2.00}


# van der Waals radii from J. Phys. Chem. A, Vol. 113, No. 19, 2009
vdwdict = { 1:1.10,  2:1.40,  3:1.81,  4:1.53,  5:1.92,  6:1.70,  7:1.55,  8:1.52,  9:1.47, 10:1.54,
           11:2.27, 12:1.73, 13:1.84, 14:2.10, 15:1.80, 16:1.80, 17:1.75, 18:1.88, 19:2.75, 20:2.31,
                                                                         
           31:1.87, 32:2.11, 33:1.85, 34:1.90, 35:1.83, 36:2.02, 37:3.03, 38:2.49, 
                                                                                   49:1.93, 50:2.17,
           51:2.06, 52:2.06, 53:1.98, 54:2.16, 55:3.43, 56:2.58,
           
                                                                          
           81:1.96, 82:2.02, 83:2.07, 84:1.97, 85:2.02, 86:2.20, 87:3.48, 88:2.83}
#@profile
def ExtractData(fpath, au=False, density=1, value_type='esp', units=None):
    """
    required arguments: path for a cube file (str)
    
    optional argument: atomic units (bool), density (float), value type (str), units(None/str)

    calls: -
    
    returns: a list of: the number of atoms, a list of the coordinates of the origin, the number of values recorded at each point, the number of voxels along x axis,
    a list of the coordinates of the increment vector for the x axis, the number of voxels along y axis, a list of the coordinates of the increment vector for the y
    axis, the number of voxels along z axis, a list of the coordinates of the increment vector for the z axis, a list of lists containing the atomic number, the nuclear
    charge and coordinates for each atom, a list of MO indices, a 3D array of all the values in the cube file, a list of all the values in the cube file (in order)

    this function extracts all the data from the cube file
    """
#    start_time=time.time() 
    
    if value_type.lower()=='esp' or value_type.lower()=='potential':
        value_type='esp'
        if au:
            f = 1  # lengths in atomic units (i.e. bohr)
            vf = 1
            units = 'a.u.'
        else:
            f = 0.529177  # length: 1 a.u. = 0.529177 angstrom
            if units is None or units.lower()=='v':
                vf = 27.21138602 # potential: 1 a.u. = 27.21138602 V
                units='V'
            elif units.lower()=='':
                pass
            else:
                vf = 27.21138602
                units='V'
                print('\nWARNING: The units you have specified were not recognised! The values will be given in '+units+'.\n')
            
    elif value_type.lower()=='dens' or value_type.lower()=='density':
        value_type='dens'
        if au:
            f = 1  # lengths in atomic units (i.e. bohr)
            vf = 1
            units = 'a.u.'
        else:
            f = 0.529177  # length: 1 a.u. = 0.529177 angstrom
            if units is None or units.lower()=='a' or units.lower()=='angstrom':
                vf = 6.74834256323  # 1 a.u. = 6.74834256323 A^(-3)
                units = '1/'+r'$\AA^3$'
            elif units.lower()=='nm':
                vf = 1000*6.74834256323  # 1 a.u. = 1000*6.74834256323 nm^(-3)
                units = '1/'+r'$nm^3$'
            elif units.lower()=='pm':
                vf = 1/10**6*6.74834256323  # 1 a.u. = 1/10**6*6.74834256323 pm^(-3)
                units = '1/'+r'$pm^3$'
            elif units.lower()=='m':
                vf = 10**30*6.74834256323  # 1 a.u. = 10**30*6.74834256323 m^(-3)
                units = '1/'+r'$m^3$'
            else:
                vf = 6.74834256323 # 1 a.u. = 6.74834256323 A^(-3)
                units = '1/'+r'$\AA^3$'
                print('\nWARNING: The units you have specified were not recognised! The values will be given in '+units+'.\n')

    with open(fpath) as file:
        lines = file.readlines()

    density=float(density)
    if density<1:
        print('\n!!! Warning: the density you chose is lower than 1. A value of 1 will be used instead. !!!')
        density=1
    elif not density.is_integer():
        print('\n!!! Warning: the density you chose is not an integer. A value of '+str(int(round(density)))+' will be used instead. !!!')
    density=int(round(density))

    # with header
    try:

        # number of atoms
        natoms = int(lines[2].split()[0])

        # origin coordinates
        o_x = f*float(lines[2].split()[1])
        o_y = f*float(lines[2].split()[2])
        o_z = f*float(lines[2].split()[3])

        # number of values recorded at each point
        try:
            nval = int(lines[2].split()[4])
        except IndexError:
            nval = 1

        # number of voxels along each axis and the increment vector for each axis in the original cube file
        n_x = int(lines[3].split()[0])
        x_x = f*float(lines[3].split()[1])
        y_x = f*float(lines[3].split()[2])
        z_x = f*float(lines[3].split()[3])

        n_y = int(lines[4].split()[0])
        x_y = f*float(lines[4].split()[1])
        y_y = f*float(lines[4].split()[2])
        z_y = f*float(lines[4].split()[3])

        n_z = int(lines[5].split()[0])
        x_z = f*float(lines[5].split()[1])
        y_z = f*float(lines[5].split()[2])
        z_z = f*float(lines[5].split()[3])
        
        # number of voxels along each axis and the increment vector for each axis after reducing the density
        d_n_x = (n_x-1)//density+1
        d_x_x = x_x*density
        d_y_x = y_x*density
        d_z_x = z_x*density

        d_n_y = (n_y-1)//density+1
        d_x_y = x_y*density
        d_y_y = y_y*density
        d_z_y = z_y*density

        d_n_z = (n_z-1)//density+1
        d_x_z = x_z*density
        d_y_z = y_z*density
        d_z_z = z_z*density

        # a list of a number of natoms lists, each list component containing the atomic number, nuclear charge and coordinates of each atom
        atoms = []
        for i in range(6, 6+abs(natoms)):
            atoms.append([int(lines[i].split()[0]), float(lines[i].split()[1]), f*float(lines[i].split()[2]), f*float(lines[i].split()[3]), f*float(lines[i].split()[4])])

    # no header
    except IndexError:
        # number of atoms
        natoms = int(lines[0].split()[0])

        # origin coordinates
        o_x = f*float(lines[0].split()[1])
        o_y = f*float(lines[0].split()[2])
        o_z = f*float(lines[0].split()[3])

        # number of values recorded at each point
        try:
            nval = int(lines[0].split()[4])
        except IndexError:
            nval = 1

        # number of voxels along each axis and the increment vector for each axis
        n_x = int(lines[1].split()[0])
        x_x = f*float(lines[1].split()[1])
        y_x = f*float(lines[1].split()[2])
        z_x = f*float(lines[1].split()[3])

        n_y = int(lines[2].split()[0])
        x_y = f*float(lines[2].split()[1])
        y_y = f*float(lines[2].split()[2])
        z_y = f*float(lines[2].split()[3])

        n_z = int(lines[3].split()[0])
        x_z = f*float(lines[3].split()[1])
        y_z = f*float(lines[3].split()[2])
        z_z = f*float(lines[3].split()[3])

        # number of voxels along each axis and the increment vector for each axis after reducing the density
        d_n_x = (n_x-1)//density+1
        d_x_x = x_x*density
        d_y_x = y_x*density
        d_z_x = z_x*density

        d_n_y = (n_y-1)//density+1
        d_x_y = x_y*density
        d_y_y = y_y*density
        d_z_y = z_y*density

        d_n_z = (n_z-1)//density+1
        d_x_z = x_z*density
        d_y_z = y_z*density
        d_z_z = z_z*density

        # a list of a number of natoms lists, each list component containing the atomic number, nuclear charge and coordinates of each atom
        atoms = []
        for i in range(4, 4+abs(natoms)):
            atoms.append([int(lines[i].split()[0]), float(lines[i].split()[1]), f*float(lines[i].split()[2]), f*float(lines[i].split()[3]), f*float(lines[i].split()[4])])
    i = i+1
    mo_list = []
    n_mos = 1
    if natoms < 0:
        n_mos = int(lines[i].split()[0])
        mo_list.extend([int(id) for id in lines[i].split()[1:]])
        k = len(lines[i].split())-1   # a counter for the number of MO ids
        i = i+1
        while k < n_mos:
            mo_list.extend([int(id) for id in lines[i].split()])
            k = k+len(lines[i].split())
            i = i+1

        natoms = abs(natoms)
    print('\n\n '+'-'*74+'\n|'+' Extracting data from the cube file '.center(74)+'|\n '+'-'*74)
    j = 0  # current position in line
    val_list = []     # a list of all the values in the cube file
    
    if nval == 1 and n_mos == 1:
        val = [[[None for m in range(d_n_z)] for l in range(d_n_y)] for k in range(d_n_x)]    # a 3d array for all the values in the cube file
        for k in range(n_x):
            for l in range(n_y):
                for m in range(n_z):
                    current_id=m+l*n_z+k*n_z*n_y+1
                    if (n_x*n_y*n_z)>=50 and current_id%((n_x*n_y*n_z)//50)==0 or (n_x*n_y*n_z)<50:
                        print('Reading data: '.ljust(21)+'#'*(current_id*50//(n_x*n_y*n_z))+' '*(54-current_id*50//(n_x*n_y*n_z)-len(str(current_id*100//(n_x*n_y*n_z))))+str(current_id*100//(n_x*n_y*n_z))+'%',end='\r')
                    if k%density==0 and l%density==0 and m%density==0:
                        val[k//density][l//density][m//density] = vf*float(lines[i].split()[j])
                        val_list.append(vf*float(lines[i].split()[j]))
                    if j == len(lines[i].split())-1:    # if reached the end of the line
                        i = i+1
                        j = 0
                    else:
                        j = j+1
    else:
        n_data = nval*n_mos
        val = [[[[] for m in range(d_n_z)] for l in range(d_n_y)] for k in range(d_n_x)]    # a 3d array for all the values in the cube file
        for k in range(n_x):
            for l in range(n_y):
                for m in range(n_z):
                    current_id=m+l*n_z+k*n_z*n_y+1
                    if (n_x*n_y*n_z)>=50 and current_id%((n_x*n_y*n_z)//50)==0 or (n_x*n_y*n_z)<50:
                        print('Reading data: '.ljust(21)+'#'*(current_id*50//(n_x*n_y*n_z))+' '*(54-current_id*50//(n_x*n_y*n_z)-len(str(current_id*100//(n_x*n_y*n_z))))+str(current_id*100//(n_x*n_y*n_z))+'%',end='\r')
                    if k%density==0 and l%density==0 and m%density==0:
                        for n in range(n_data):
                            val[k//density][l//density][m//density].append(vf*float(lines[i].split()[j]))
                            if j == len(lines[i].split())-1:    # if reached the end of the line
                                i = i+1
                                j = 0
                            else:
                                j = j+1
                        val_list.append(val[k//density][l//density][m//density])
    print('Reading data: '.ljust(21)+50*'#'+' 100%')
          
#    print("--- %s seconds ---" % (time.time() - start_time))
    return [natoms, [o_x, o_y, o_z], nval, d_n_x, [d_x_x, d_y_x, d_z_x], d_n_y, [d_x_y, d_y_y, d_z_y], d_n_z, [d_x_z, d_y_z, d_z_z], atoms, mo_list, val, val_list, units]


def Atoms(fpath=None, au=False, value_type='esp', units=None):
    """
    arguments: path for a cube file (str), atomic units (bool), value type (str), units(None/str)

    calls: ExtractData

    returns: a list of a number of natoms lists, each list component containing the atomic number, nuclear charge and coordinates of each atom
    """

    return ExtractData(fpath, au=au, value_type=value_type, units=units)[9]

def Axes(fpath=None, origin = None, n_x = None, x_vector = None, n_y = None, y_vector = None, n_z = None, z_vector = None, au=False, density=1, value_type='esp', units=None):
    """
    arguments: path for a cube file (str), origin (list), number of voxels along x axis (float), increment vector in the x direction (list), number of voxels along y axis (float), increment vector in the y direction (list),
               number of voxels along z axis (float), increment vector in the z direction (list), atomic units (bool), density (float), value type (str), units(None/str)

    calls: ExtractData (if fpath is specified)
    
    called with either [fpath, au, density, value_type, units] or [origin, n_x, x_vector, n_y, y_vector, n_z, z_vector]
    
    If specified, fpath takes priority over the other arguments and ExtractData is called (i.e. the script reads and goes through the whole cube file)

    returns: a list of three lists representing the x, y and z coordinates, respectively of all the grid points
    """
    if fpath:
        exdata = 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]
    [o_x, o_y, o_z] = origin
    [x_x, y_x, z_x] = x_vector
    [x_y, y_y, z_y] = y_vector
    [x_z, y_z, z_z] = z_vector

    x = []
    for i in range(n_x):
        x.append([o_x+i*x_x,o_y+i*y_x,o_z+i*z_x])

    y = []
    for i in range(n_y):
        y.append([o_x+i*x_y,o_y+i*y_y,o_z+i*z_y])

    z = []
    for i in range(n_z):
        z.append([o_x+i*x_z,o_y+i*y_z,o_z+i*z_z])

    return [x, y, z]


def ValuesAsMatrix(fpath=None, au=False, density=1, value_type='esp', units=None):
    """
    arguments: path for a cube file (str), atomic units (bool), density (float), value type (str), units(None/str)

    calls: ExtractData

    returns: a 3D array of all the values in the cube file
    """

    return ExtractData(fpath, au=au, value_type=value_type, density=density, units=units)[11]


def ValuesAsList(fpath=None, au=False, density=1, value_type='esp', units=None):
    """
    arguments: path for a cube file (str), atomic units (bool), density (float), value type (str), units(None/str)

    calls: ExtractData

    returns: a list of all the values in the cube file (in the order in which they appear)
    """

    return ExtractData(fpath, au=au, value_type=value_type, density=density, units=units)[12]

def ValuesAsDictionary(fpath=None, vals=None, axes_list=None, origin=None, au=False, density=1, value_type='esp', units=None):
    """
    arguments: path for a cube file (str), values (list), axes list (list), origin (list), atomic units (bool), density (float), value type (str), units(None/str)

    calls: ExtractData and Axes (if fpath is specified)
    
    called with either [fpath, au, density, value_type, units] or [vals, axes_list, origin]
    
    If specified, fpath takes priority over the other arguments and ExtractData is called (i.e. the script reads and goes through the whole cube file)

    returns: a list of an ordered dictionary of each grid point (represented through a tuple of its coordinates) and its corresponding value and
             an ordered dictionary of each grid point (represented through a tuple of its indices) and its corresponding value
    """
    
    if fpath:
        exdata = 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]
        axes_list=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)
    [x, y, z] = axes_list
        

    d_pts = cd.OrderedDict()
    d_ids = cd.OrderedDict()
    idx = 0
    
    x_id=0
    y_id=0
    z_id=0
    
    for i in x:
        for j in y:
            for k in z:
                d_pts[(i[0]+j[0]+k[0]-2*origin[0], i[1]+j[1]+k[1]-2*origin[1], i[2]+j[2]+k[2]-2*origin[2])] = vals[idx]
                d_ids[(x_id,y_id,z_id)] = vals[idx]
                idx = idx+1
                z_id += 1
            z_id = 0
            y_id += 1
        y_id=0
        x_id += 1

    return d_pts, d_ids


def Minimum(fpath=None, gridpts=None, vals=None, au=False, density=1, value_type='esp', units=None):
    """
    arguments: path for a cube file (str), grid points (dict), values (list), atomic units (bool), density (float), value type (str), units(None/str)

    calls: ExtractData, Axes and ValuesAsDictionary (if fpath is specified)
    
    called with either [fpath, au, density, value_type, units] or [gridpts, vals]
    
    If specified, fpath takes priority over the other arguments and ExtractData is called (i.e. the script reads and goes through the whole cube file)

    returns: a list the coordinates of all the grid points that have the minimum value
    """
    if fpath:
        exdata = 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]
        axes_list=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=ValuesAsDictionary(vals=vals,axes_list=axes_list, origin=origin)[0]
                
    all_pts_list = list(gridpts.items())
    m = min(vals)
    min_list = []
    for a in all_pts_list:
        if a[1] == m:
            min_list.append(a[0])

    return min_list


def Maximum(fpath=None, gridpts=None, vals=None, au=False, density=1, value_type='esp', units=None):
    """
    arguments: path for a cube file (str), grid points (dict), values (list), atomic units (bool), density (float), value type (str), units(None/str)

    calls: ExtractData, Axes and ValuesAsDictionary (if fpath is specified)
    
    called with either [fpath, au, density, value_type, units] or [gridpts, vals]
    
    If specified, fpath takes priority over the other arguments and ExtractData is called (i.e. the script reads and goes through the whole cube file)

    returns: a list the coordinates of all the grid points that have the maximum value
    """
    if fpath:
        exdata = 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]
        axes_list=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=ValuesAsDictionary(vals=vals,axes_list=axes_list, origin=origin)[0]

    all_pts_list = list(gridpts.items())
    m = max(vals)
    max_list = []
    for a in all_pts_list:
        if a[1] == m:
            max_list.append(a[0])

    return max_list


def GetVdWPoints(fpath=None, axes_list=None, gridpts=None, atoms_list=None, x_vector=None, y_vector=None, z_vector=None, origin=None, factor=1, au=False, density=1, value_type='esp', units=None, cat=False, an=False):
    """
    arguments: path for a cube file (str), axes list (list), grid points (dict), atoms list (list), increment vector in the x direction (list), origin (list), factor (float), atomic units (bool), density (float), value type (str), units(None/str)

    calls: ExtractData, Axes and ValuesAsDictionary (if fpath is specified)
    
    called with either [fpath, au, density, value_type, units, factor] or [axes_list, gridpts, atoms_list, x_vector, y_vector, z_vector, origin, factor]
    
    If specified, fpath takes priority over the other arguments and ExtractData is called (i.e. the script reads and goes through the whole cube file)

    returns: a list of three lists: the first one contains the x, y, z coordinates and the value of each point on the VdW surface, the second one contains the indices and the value of each point on the VdW surface and
    the third one contains the coordinates and values of all the points that are on or inside the VdW surface
    
    the points that reside on the VdW surface are chosen to be within plus/minus a fractional distance of the surface defined by the Van der Waals radii of the constituent atoms
    (the fractional distance represents 30% of the distance between two diagonally adjacent points in the grid)

    this function returns all the points that are close to the surface defined in terms of the Van der Waals radii of the constituent atoms
    """
#    start_time=time.time()
    
    if fpath:
        exdata = 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]
        axes_list=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=ValuesAsDictionary(vals=vals,axes_list=axes_list, origin=origin)[0]
        
    [x, y, z] = axes_list
    dgp = gridpts
    [x_x, y_x, z_x] = x_vector
    [x_y, y_y, z_y] = y_vector
    [x_z, y_z, z_z] = z_vector
    n_x=len(x)
    n_y=len(y)
    n_z=len(z)
    dist = ((x_x+x_y+x_z)**2+(y_x+y_y+y_z)**2+(z_x+z_y+z_z)**2)**0.5            
    
    half_dist = dist*0.4

    coords = []
    id_coords = []
    filled_coords = []

    if au:
        vdw_f = 0.529177  # lengths in atomic units (i.e. bohr)
    else:
        vdw_f = 1  # lengths in angstroms
    print('\n\n '+'-'*74+'\n|'+' Building the van der Waals surface '.center(74)+'|\n '+'-'*74)
    for at in atoms_list:
        vdw = vdwdict[at[0]]/vdw_f
    for i in range(n_x):
        for j in range(n_y):
            for k in range(n_z):
                current_id=k+j*n_z+i*n_z*n_y+1
                if (n_x*n_y*n_z)>=50 and current_id%((n_x*n_y*n_z)//50)==0 or (n_x*n_y*n_z)<50:
                    print('Selecting points: '.ljust(21)+'#'*(current_id*50//(n_x*n_y*n_z))+' '*(54-current_id*50//(n_x*n_y*n_z)-len(str(current_id*100//(n_x*n_y*n_z))))+str(current_id*100//(n_x*n_y*n_z))+'%',end='\r')
                is_on_vdw_surf = False
                for at in atoms_list:
                    if True:  # a condition can be added here to create a surface only around some atoms (for instance, at[0]==7 if you only want to visualise the nitrogen atoms)
                        is_cat=at[0] in [1,7] or (at[0]==6 and at[2]<0)
                        is_an=not(is_cat)
                        if (cat and is_cat) or (an and is_an) or (not cat and not an): 
                            vdw = vdwdict[at[0]]/vdw_f
                            d = ((x[i][0]+y[j][0]+z[k][0]-2*origin[0]-at[2])**2+(x[i][1]+y[j][1]+z[k][1]-2*origin[1]-at[3])**2+(x[i][2]+y[j][2]+z[k][2]-2*origin[2]-at[4])**2)**0.5   # distance from the current point to the center of the atom
                            if d < vdw*factor+half_dist:  # if within +/- the half distance between any two diagonally adjacent points in the grid
                                filled_coords.append([x[i][0]+y[j][0]+z[k][0]-2*origin[0], x[i][1]+y[j][1]+z[k][1]-2*origin[1], x[i][2]+y[j][2]+z[k][2]-2*origin[2], dgp[(x[i][0]+y[j][0]+z[k][0]-2*origin[0], x[i][1]+y[j][1]+z[k][1]-2*origin[1], x[i][2]+y[j][2]+z[k][2]-2*origin[2])]])
                                if d > vdw*factor-half_dist:  # if within +/- the half distance between any two diagonally adjacent points in the grid
                                    is_on_vdw_surf = True
                                    break
                if is_on_vdw_surf:
                    is_inside = False    # assume it is inside the vdw radius for another atom
                    for at in atoms_list:
                        vdw = vdwdict[at[0]]/vdw_f
                        d = ((x[i][0]+y[j][0]+z[k][0]-2*origin[0]-at[2])**2+(x[i][1]+y[j][1]+z[k][1]-2*origin[1]-at[3])**2+(x[i][2]+y[j][2]+z[k][2]-2*origin[2]-at[4])**2)**0.5   # distance from the current point to the center of the atom
                        if d < vdw*factor-half_dist:  # if the point is inside the vdw radius for another atom
                            is_inside = True
                            break
                    if not is_inside:  # if not inside the vdw radius for another atom
                        coords.append([x[i][0]+y[j][0]+z[k][0]-2*origin[0], x[i][1]+y[j][1]+z[k][1]-2*origin[1], x[i][2]+y[j][2]+z[k][2]-2*origin[2], dgp[(x[i][0]+y[j][0]+z[k][0]-2*origin[0], x[i][1]+y[j][1]+z[k][1]-2*origin[1], x[i][2]+y[j][2]+z[k][2]-2*origin[2])]])
                        id_coords.append([i, j, k, dgp[(x[i][0]+y[j][0]+z[k][0]-2*origin[0], x[i][1]+y[j][1]+z[k][1]-2*origin[1], x[i][2]+y[j][2]+z[k][2]-2*origin[2])]])
    print('Selecting points: '.ljust(21)+50*'#'+' 100%')
#    print("--- %s seconds ---" % (time.time() - start_time))
    return coords, id_coords, filled_coords



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)
