"""     

```
mddf_linkedcells_parallel(trajectory::Trajectory, options::Options,samples::Samples, RNG, mddf_compute!)  
```

Computes the MDDF using linked cells, in parallel.

"""
function mddf_linkedcells_parallel(
    trajectory::Trajectory,
    options::Options,
    samples::Samples,
    RNG,
    mddf_compute!,
)

    # Simplify code by assigning some shortened names
    solute = trajectory.solute
    solvent = trajectory.solvent
    x_solute = trajectory.x_solute
    x_solvent = trajectory.x_solvent

    # Number of threads
    if options.nthreads < 0
        nspawn = Threads.nthreads() - 1
    else
        nspawn = options.nthreads - 1
    end
    if nspawn == 0
        error(" Parallel version must be executed only with more than 1 thread. ")
    end

    # Initializing the structure that carries the result per thread
    R0 = Result(trajectory, options)
    R = [Result(trajectory, options, irefatom = R0.irefatom) for i = 1:nspawn]

    # Safe passing of frame counter to the threads
    tframe = zeros(Int, nspawn)

    # Create data structures required for multithreading
    framedata = Vector{FrameData}(undef, nspawn)
    for ispawn = 1:nspawn
        framedata[ispawn] = FrameData(deepcopy(trajectory), R[1])
    end

    # Tasks
    t = Vector{Task}(undef, nspawn)
    free = ones(Bool, nspawn)

    # Skip initial frames if desired
    iframe = 0 # Counter for all frames of the input file
    while iframe < options.firstframe - 1
        nextframe!(trajectory)
        iframe += 1
    end

    # Print some information about this run
    options.silent || title(R[1], solute, solvent, nspawn)

    ndone = 0
    iframe_read = 0 # Counter for the frames that are actually being considered
    if !options.silent
        progress = Progress(R[1].nframes_read, 1)
    end
    while ndone < R[1].nframes_read

        # Launch for each free thread the computation of a frame
        while iframe_read < R[1].nframes_read && count(free) > 0
            ifree = findfirst(x -> x == true, free)
            # Skip frames that are not going to be considered 
            while (iframe + 1) % options.stride != 0
                nextframe!(trajectory)
                iframe += 1
            end
            # Actually read the data of the frame that will be considered
            nextframe!(trajectory)
            iframe += 1
            iframe_read += 1
            tframe[ifree] = iframe
            @. framedata[ifree].trajectory.x_solute = trajectory.x_solute
            @. framedata[ifree].trajectory.x_solvent = trajectory.x_solvent
            @. framedata[ifree].trajectory.sides = trajectory.sides
            # Spawn the calculations for this frame
            t[ifree] = ThreadPools.@tspawnat ifree + 1 mddf_compute!(
                tframe[ifree],
                framedata[ifree],
                options,
                RNG,
                R[ifree],
            )
            free[ifree] = false
        end

        # Wait a little bit before checking
        sleep(options.sleep)

        # Check thread status
        for ispawn = 1:nspawn
            if !free[ispawn]
                if istaskfailed(t[ispawn])
                    error(
                        " Computation of MDDF failed in thread: $ispawn",
                        fetch(t[ispawn]),
                    )
                end
                if istaskdone(t[ispawn])
                    ndone += 1
                    free[ispawn] = true
                    options.silent || next!(progress)
                    if options.GC &&
                       (Sys.free_memory() / Sys.total_memory() < options.GC_threshold)
                        GC.gc() # why we need this anyway??? There should not be so much garbage.
                    end
                end
            end
        end

    end
    closetraj(trajectory)

    # Sum up the results of all threads into the data of thread one (R1<-R1+R2)
    for ispawn = 2:nspawn
        sum!(R[1], R[ispawn])
    end

    # Setup the final data structure with final values averaged over the number of frames,
    # sampling, etc, and computes final distributions and integrals
    finalresults!(R[1], options, trajectory, samples)
    options.silent || println(bars)

    return R[1]

end
