get_atom_species Subroutine

subroutine get_atom_species(this, atoms_pcode, species, num_atoms, num_species, na, sp_map)

    class(mimic_communicator), intent(inout) :: this
    integer, dimension(:), intent(in) :: atoms_pcode
    !> atom species of each atom in the system (-1 for overlapped)
    integer, dimension(:,:), intent(inout), allocatable, target :: species
    integer, dimension(:), intent(inout) :: na
    integer, intent(out) :: num_atoms
    integer, intent(out) :: num_species
    integer, dimension(:,:), allocatable, intent(out) :: sp_map

    integer :: n_client, n_map, n_species, n_atom
    integer :: code_atoms
    integer, dimension(this%num_clients) :: species_pcode

    num_species = 0
    num_atoms = 0

    do n_client = 1, this%num_clients
        call this%send_command(MCL_SEND_ATOM_SPECIES, n_client)
        call mcl_receive(species(:, n_client), atoms_pcode(n_client), MCL_DATA, n_client)
        species_pcode(n_client) = maxval(species(:, n_client))
    end do ! n_client

    allocate(sp_map(maxval(species_pcode), this%num_clients))
    sp_map(:, :) = -1

    ! Treat overlaps
    do n_client = 1, this%num_clients
        do n_atom = 1, atoms_pcode(n_client)
            do n_map = 1, size(this%overlap_maps(n_client)%maps)
                if (this%overlap_maps(n_client)%maps(n_map)%id == n_atom) then
                    species(n_atom, n_client) = -1
                end if
            end do ! n_map
        end do ! n_atom

        ! Count atoms/species
        do n_species  = 1, species_pcode(n_client)
            code_atoms = size(pack(species(:, n_client), mask=(species(:, n_client) == n_species)))
            if (code_atoms > 0) then
                num_species = num_species + 1
                na(num_species) = code_atoms
                sp_map(n_species, n_client) = num_species
            end if

            if (code_atoms > num_atoms) then
                num_atoms = code_atoms
            end if
        end do ! n_species
    end do ! n_client

end subroutine get_atom_species