gather_bonds Subroutine

subroutine gather_bonds(this, bonds, bonds_number)

    class(mimic_communicator), intent(inout) :: this
    !> list of constrained bonds in the simulation, return value
    type(bond_type), dimension(:,:), allocatable, intent(out) :: bonds
    !> number of bond constraints in the system.
    !> WARNING!!! this routine updates this value
    integer, dimension(:), intent(inout) :: bonds_number

    integer, dimension(:, :), target, allocatable :: indices
    real(dp), dimension(:, :), target, allocatable :: lengths

    integer, dimension(:, :), allocatable :: act_ind
    integer :: n_client, n_bond, n_map
    integer :: offset
    integer :: overlap_id
    integer, dimension(:), allocatable :: temp_bonds_number
    logical :: overlapped

    allocate(indices(2 * maxval(bonds_number), this%num_clients))
    allocate(act_ind(maxval(bonds_number), this%num_clients))
    allocate(lengths(maxval(bonds_number), this%num_clients))
    allocate(temp_bonds_number(this%num_clients))

    if (maxval(bonds_number) == 0) return

    temp_bonds_number = bonds_number
    do n_client = 1, this%num_clients
        if (bonds_number(n_client) == 0) cycle

        offset = 1

        call this%send_command(MCL_SEND_BOND_ATOMS, n_client)
        call mcl_receive(indices(:, n_client), 2 * bonds_number(n_client), MCL_DATA, n_client)
        call this%send_command(MCL_SEND_BOND_LENGTHS, n_client)
        call mcl_receive(lengths(:, n_client), bonds_number(n_client), MCL_DATA, n_client)

        do n_bond = 1, bonds_number(n_client)
            overlapped = .false.
            do n_map = 1, size(this%overlap_maps(n_client)%maps)
                overlap_id = this%overlap_maps(n_client)%maps(n_map)%id
                if (overlap_id == indices(n_bond * 2 - 1, n_client) .or. &
                    overlap_id == indices(n_bond * 2, n_client)) then
                    temp_bonds_number(n_client) = temp_bonds_number(n_client) - 1
                    overlapped = .true.
                    exit
                end if
            end do ! n_map
            if (.not. overlapped) then
                act_ind(offset, n_client) = n_bond
                offset = offset + 1
            endif
        end do ! n_bond
    end do ! n_client
    bonds_number = temp_bonds_number
    allocate(bonds(maxval(bonds_number), this%num_clients))

    do n_client = 1, this%num_clients
        do n_bond = 1, bonds_number(n_client)
            bonds(n_bond, n_client)%atom_i = &
                indices(act_ind(n_bond, n_client) * 2 - 1, n_client)
            bonds(n_bond, n_client)%atom_j = &
                indices(act_ind(n_bond, n_client) * 2, n_client)
            bonds(n_bond, n_client)%length = &
                lengths(act_ind(n_bond, n_client), n_client)
        end do ! n_bond
    end do ! n_client

end subroutine gather_bonds