def find_windows(elements,
                 coordinates,
                 processes=None,
                 mol_size=None,
                 adjust=1,
                 pore_opt=True,
                 increment=1.0,
                 **kwargs):
    """Return windows diameters and center of masses for a molecule."""
    # Copy the coordinates as will perform many opertaions on them
    coordinates = deepcopy(coordinates)
    # Center of our cartesian system is always at origin
    origin = np.array([0, 0, 0])
    # Initial center of mass to reverse translation at the end
    initial_com = center_of_mass(elements, coordinates)
    # Shift the cage to the origin using either the standard center of mass
    # or if pore_opt flag is True, the optimised pore center as center of mass
    if pore_opt is True:
        # Normally the pore is calculated from the COM of a molecule.
        # So, esentially the molecule's COM is the pore center.
        # To shift the molecule so that the center of the optimised pore
        # is at the origin of the system and not the center of the not
        # optimised one, we need to adjust the shift. We also have to update
        # the initial com.
        com_adjust = initial_com - opt_pore_diameter(elements, coordinates, **
                                                     kwargs)[2]
        initial_com = initial_com - com_adjust
        coordinates = shift_com(elements, coordinates, com_adjust=com_adjust)
    else:
        # Otherwise, we just shift the cage to the origin.
        coordinates = shift_com(elements, coordinates)
    # We create an array of vdw radii of elements.
    elements_vdw = np.array([[atomic_vdw_radius[x.upper()]] for x in elements])
    # We calculate maximum diameter of a molecule to determine the radius
    # of a sampling sphere neccessary to enclose the whole molecule.
    shpere_radius = max_dim(elements, coordinates)[2] / 2
    sphere_surface_area = 4 * np.pi * shpere_radius**2
    # Here we determine the number of sampling points necessary for a fine
    # sampling. Smaller molecules require more finner density of sampling
    # points on the sampling sphere's surface, whereas largen require less.
    # This formula was created so that larger molecule do not take much longer
    # to analyse, as number_sampling_points*length_of_sampling_vectors
    # results in quadratic increase of sampling time. The 250 factor was
    # specificly determined to produce close to 1 sampling point /Angstrom^2
    # for a sphere of radius ~ 24 Angstrom. We can adjust how fine is the
    # sampling by changing the adjust factor.
    number_of_points = int(np.log10(sphere_surface_area) * 250 * adjust)
    # Here I use code by Alexandre Devert for spreading points on a sphere:
    # http://blog.marmakoide.org/?p=1
    golden_angle = np.pi * (3 - np.sqrt(5))
    theta = golden_angle * np.arange(number_of_points)
    z = np.linspace(1 - 1.0 / number_of_points, 1.0 / number_of_points - 1.0,
                    number_of_points)
    radius = np.sqrt(1 - z * z)
    points = np.zeros((number_of_points, 3))
    points[:, 0] = radius * np.cos(theta) * shpere_radius
    points[:, 1] = radius * np.sin(theta) * shpere_radius
    points[:, 2] = z * shpere_radius
    # Here we will compute the eps parameter for the sklearn.cluster.DBSCAN
    # (3-dimensional spatial clustering algorithm) which is the mean distance
    # to the closest point of all points.
    values = []
    tree = KDTree(points)
    for i in points:
        dist, ind = tree.query(i.reshape(1, -1), k=10)
        values.extend(dist)
    mean_distance = np.mean(values)
    # The best eps is parametrized when adding the mean distance and it's root.
    eps = mean_distance + mean_distance**0.5
    # Here we either run the sampling points vectors analysis in serial
    # or parallel. The vectors that go through molecular pores return
    # as analysed list with the increment at vector's path with largest
    # included sphere, coordinates for this narrow channel point. vectors
    # that find molecule on theirs path are return as NoneType object.
    # Parralel analysis on user's defined number of CPUs.
    if processes:
        pool = Pool(processes=processes)
        parallel = [
            pool.apply_async(
                vector_preanalysis,
                args=(
                    point,
                    coordinates,
                    elements_vdw, ),
                kwds={'increment': increment}) for point in points
        ]
        results = [p.get() for p in parallel if p.get() is not None]
        pool.terminate()
        # Dataset is an array of sampling points coordinates.
        dataset = np.array([x[5:8] for x in results])
    else:
        results = [
            vector_preanalysis(
                point, coordinates, elements_vdw, increment=increment)
            for point in points
        ]
        results = [x for x in results if x is not None]
        dataset = np.array([x[5:8] for x in results])
    # If not a single vector was returned from the analysis it mean that
    # no molecular channels (what we call windows here) connects the
    # molecule's interior with the surroungsings (exterior space).
    # The number of windows in that case equals zero and zero is returned.
    # Otherwise we continue our search for windows.
    if len(results) == 0:
        return None
    else:
        # Perfomr DBSCAN to cluster the sampling points vectors.
        # the n_jobs will be developed later.
        # db = DBSCAN(eps=eps, n_jobs=_ncpus).fit(dataset)
        db = DBSCAN(eps=eps).fit(dataset)
        core_samples_mask = np.zeros_like(db.labels_, dtype=bool)
        core_samples_mask[db.core_sample_indices_] = True
        labels = set(db.labels_)
        # Assing cluster label to a sampling point.
        clusters = [[i, j] for i, j in zip(results, db.labels_)]
        clustered_results = {label: [] for label in labels}
        # Create a dictionary of clusters with points listed.
        [clustered_results[i[1]].append(i[0]) for i in clusters]
        # No for the sampling point vector in each cluster that had
        # the widest channel's 'neck' is assumed to pass the closest
        # to the window's center and therefore will be passed to
        # window analysis function.
        # We also pass user defined settings for window analysis.
        # Again either in serlia or in parallel.
        # Noisy points get a cluster label -1, therefore we have to exclude it.
        if processes:
            pool = Pool(processes=processes)
            parallel = [
                pool.apply_async(
                    window_analysis,
                    args=(np.array(clustered_results[cluster]), elements,
                          coordinates, elements_vdw),
                    kwds=kwargs) for cluster in clustered_results
                if cluster != -1
            ]
            window_results = [p.get() for p in parallel if p.get() is not None]
            pool.terminate()
        else:
            window_results = [
                window_analysis(
                    np.array(clustered_results[cluster]), elements,
                    coordinates, elements_vdw, **kwargs)
                for cluster in clustered_results if cluster != -1
            ]
    # The function returns two numpy arrays, one with windows diameters
    # in Angstrom, second with corresponding windows center's coordinates
    windows = np.array([result[0] for result in window_results
                        if result is not None])
    windows_coms = np.array(
        [np.add(result[1], initial_com) for result in window_results
         if result is not None])
    # Safety measures, if one of the windows is None or negative a warning
    # should be raised.
    for result in window_results:
        if result is None:
            msg_ = " ".join(
                ['Warning. One of the analysed windows has',
                 'returned as None. See manual.']
            )
            # print(msg_)
        elif result[0] < 0:
            msg_ = " ".join(
                ['Warning. One of the analysed windows has a vdW',
                 'corrected diameter smaller than 0. See manual.']
            )
            # print(msg_)
    return (windows, windows_coms)