mimic_data_collect.f90 Source File


Contents


Source Code

!
!    MiMiC: A Framework for Multiscale Modeling in Computational Chemistry
!    Copyright (C) 2015-2022  The MiMiC Contributors (see CONTRIBUTORS file for details).
!
!    This file is part of MiMiC.
!
!    MiMiC is free software: you can redistribute it and/or modify
!    it under the terms of the GNU Lesser General Public License as
!    published by the Free Software Foundation, either version 3 of
!    the License, or (at your option) any later version.
!
!    MiMiC is distributed in the hope that it will be useful, but
!    WITHOUT ANY WARRANTY; without even the implied warranty of
!    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
!    GNU Lesser General Public License for more details.
!
!    You should have received a copy of the GNU Lesser General Public License
!    along with this program.  If not, see <http://www.gnu.org/licenses/>.
!

!> author: V. Bolnykh
!> provides routines for interaction with client codes:
!> collection and distribution of data
module mimic_data_collect

    use mcl_fortran, only: mcl_init, mcl_handshake, mcl_send, &
                          mcl_receive, mcl_destroy
    use mcl_requests

    use mimic_fragments, only: fragment_type
    use mimic_precision, only: dp, i32
    use mimic_subsystems, only: subsystem_type
    use mimic_types, only: maps_type, bond_type, angle_type

    implicit none

    public :: mimic_communicator

    character, parameter :: delimiter = ";"

    !> communicator type providing routines to handle data interactions
    !> @todo consider removing it - possibly redundant
    type mimic_communicator
        !> number of client codes
        integer :: num_clients
        !> mapping of overlapping atoms
        type(maps_type), dimension(:), allocatable :: overlap_maps
    contains
        procedure :: init
        procedure :: handshake
        procedure :: send_command
        procedure :: init_overlaps
        procedure :: finalize
        procedure :: destroy_comm
        procedure :: get_atom_species
        procedure :: send_coords
        procedure :: gather_coords
        procedure :: gather_forces
        procedure :: gather_energies
        procedure :: gather_atom_fragment_count
        procedure :: gather_atom_fragment_ids
        procedure :: get_multipole_order
        procedure :: gather_multipoles
        procedure :: gather_masses
        procedure :: gather_bonds
        procedure :: gather_angles
        procedure :: gather_int
        procedure :: send_int
        procedure :: gather_elements
        procedure :: gather_atom_count
        procedure :: gather_fragment_count
        procedure :: gather_bond_count
        procedure :: gather_angle_count
        procedure :: gather_type_count
        procedure :: issue_force_comp
    end type mimic_communicator

contains

!> initialize communication
subroutine init(this, comm)

    class(mimic_communicator), intent(inout) :: this
    !> MPI communicator
    integer, intent(inout) :: comm

    call mcl_init(comm)

end subroutine init

!> identify and connect clients
subroutine handshake(this, paths)

    class(mimic_communicator), intent(inout) :: this
    !> array of paths to working folders of clients
    character(len=*), dimension(:), intent(in) :: paths

    integer :: n_client
    character(len=:), allocatable :: merged_paths
    integer :: n_char, i

    this%num_clients = size(paths)
    n_char = 0

    do n_client = 1, this%num_clients
        if (n_client > 1) then
            n_char = n_char + 1
        end if
        n_char = n_char + len_trim(paths(n_client))
    end do

    allocate(character(len = n_char) :: merged_paths)

    merged_paths = ""

    do n_client = 1, this%num_clients
        if (n_client > 1) then
            merged_paths = merged_paths // delimiter
        end if
        merged_paths = merged_paths // trim(paths(n_client))
    end do

    call mcl_handshake(merged_paths, delimiter, .true.)

end subroutine handshake

!> send a command via the CommLib
subroutine send_command(this, command, dest)
    implicit none
    class(mimic_communicator), intent(inout) :: this
    !> integer code of the command
    integer, intent(in) :: command
    !> destination
    integer, intent(in) :: dest

    call MCL_send(command, 1, MCL_COMMAND, dest)
end subroutine send_command

!> initialize overlap mapping
subroutine init_overlaps(this, overlap_maps)

    class(mimic_communicator), intent(inout) :: this
    !> list of overlap_maps
    type(maps_type), dimension(:) :: overlap_maps

    if (allocated(this%overlap_maps)) deallocate(this%overlap_maps)
    allocate(this%overlap_maps(size(overlap_maps)))
    this%overlap_maps = overlap_maps

end subroutine init_overlaps

!> Get the number of atoms per atomic species per code (with overlap treatment)
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

!> routine sending an exit code to each client
subroutine finalize(this)

    class(mimic_communicator), intent(inout) :: this

    integer :: n_client

    do n_client = 1, this%num_clients
        call this%send_command(MCL_EXIT, n_client)
    end do

end subroutine finalize

!> destructor for the commlib
subroutine destroy_comm(this)

    class(mimic_communicator), intent(inout) :: this

    call mcl_destroy
    deallocate(this%overlap_maps)
    this%num_clients = -1

end subroutine destroy_comm

!> get the number of fragment per each of the client code
subroutine get_fragment_count(this, frag_count)

    class(mimic_communicator), intent(inout) :: this
    !> number of fragments per code, return value
    integer, dimension(:), intent(inout) :: frag_count

    integer(kind=i32), target :: code_fragment_count
    integer :: n_code

    do n_code = 1, this%num_clients
        call this%send_command(MCL_SEND_NUM_FRAGMENTS, n_code)
        call mcl_receive(code_fragment_count, 1, MCL_DATA, n_code)

        frag_count(n_code) = code_fragment_count
    end do

end subroutine get_fragment_count

!> get the bond constraints data - atom indices and lengths
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

!> get the angle constraints data - atom indices and equilibrium values
subroutine gather_angles(this, angles, angle_number)

    class(mimic_communicator), intent(inout) :: this
    !> list of angle constraints, return value
    type(angle_type), dimension(:,:), allocatable, intent(out) :: angles
    !> number of angle constraints
    integer, dimension(:), intent(in) :: angle_number

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

    !integer, dimension(:, :), allocatable :: act_ind
    integer :: n_client, n_angle!, n_map
    integer :: offset
    !integer :: overlap_id
    !integer, dimension(:) :: temp_angle_number
    !logical :: overlapped

    allocate(indices(3 * maxval(angle_number), this%num_clients))
    allocate(angle_vals(maxval(angle_number), this%num_clients))
    allocate(angles(maxval(angle_number), this%num_clients))

    if (maxval(angle_number) == 0) return

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

        offset = 1

        call this%send_command(MCL_SEND_ANGLE_ATOMS, n_client)
        call mcl_receive(indices(:, n_client), 3 * angle_number(n_client), MCL_DATA, n_client)
        call this%send_command(MCL_SEND_ANGLE_VALS, n_client)
        call mcl_receive(angle_vals(:, n_client), angle_number(n_client), MCL_DATA, n_client)

        !! @todo check if we need to skip overlapped angles
        ! do n_angle = 1, angle_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_angle * 2, n_client) .or. &
        !             overlap_id == indices(n_angle * 2 + 1, n_client) .or. &
        !             overlap_id == indices(n_angle * 2 + 2, n_client)) then
        !             temp_angle_number(n_client) = temp_angle_number(n_client) - 1
        !             overlapped = .true.
        !             exit
        !         end if
        !     end do ! n_map
        !     if (.not. overlapped) then
        !         act_ind(offset, n_client) = n_angle
        !         offset = offset + 1
        !     endif
        ! end do ! n_angle
    end do ! n_client

    do n_client = 1, this%num_clients
        do n_angle = 1, angle_number(n_client)
            angles(n_angle, n_client)%atom_i = &
                indices(n_angle * 3 - 2, n_client)
            angles(n_angle, n_client)%atom_j = &
                indices(n_angle * 3 - 1, n_client)
            angles(n_angle, n_client)%atom_k = &
                indices(n_angle * 3, n_client)
            angles(n_angle, n_client)%angle = &
                angle_vals(n_angle, n_client)
        end do ! n_angle
    end do ! n_client

end subroutine gather_angles

!> gather element number associated to MM atom types
subroutine gather_elements(this, elements_pcode, elements, species_map)

    class(mimic_communicator), intent(inout) :: this
    !> number of elements to expect per code
    integer, dimension(:), intent(in) :: elements_pcode
    !> element numbers per code
    integer, dimension(:), intent(out) :: elements
    integer, dimension(:,:), intent(in) :: species_map

    integer :: n_client, n_species
    integer, dimension(maxval(elements_pcode), this%num_clients), target :: temp_el

    elements(:) = -1
    temp_el(:, :) = -1

    do n_client = 1, this%num_clients
        call this%send_command(MCL_SEND_ATOM_ELEMENTS, n_client)
        call mcl_receive(temp_el(:, n_client), elements_pcode(n_client), MCL_DATA, n_client)

        do n_species = 1, elements_pcode(n_client)
            if (species_map(n_species, n_client) /= - 1) then
                elements(species_map(n_species, n_client)) = temp_el(n_species, n_client)
            end if
        end do ! n_species
    end do ! n_client

end subroutine gather_elements

!> get the number of atoms per each fragment in the system
subroutine gather_atom_fragment_count(this, fragment_count, fragment_data)

    class(mimic_communicator), intent(inout) :: this
    !> number of fragments per each code
    integer, dimension(:), intent(in) :: fragment_count
    !> number of atoms per fragment per code
    integer, dimension(:, :), intent(inout) :: fragment_data

    integer :: n_code
    integer(kind=i32), dimension(:), allocatable, target :: atoms_fragment

    fragment_data(:, :) = -1

    do n_code = 1, this%num_clients
        allocate(atoms_fragment(fragment_count(n_code)))
        call this%send_command(MCL_SEND_NUM_ATOMS_IN_FRAGMENTS, n_code)
        call mcl_receive(atoms_fragment, fragment_count(n_code), MCL_DATA, n_code)
        fragment_data(1 : size(atoms_fragment), n_code) = atoms_fragment
        deallocate(atoms_fragment)
    end do

end subroutine gather_atom_fragment_count

!> get IDs of atoms per each fragment in the system
subroutine gather_atom_fragment_ids(this, atoms_pcode, fragment_data, fragment_ids)

    class(mimic_communicator), intent(inout) :: this
    !> number of fragments per each code
    integer, dimension(:), intent(in) :: atoms_pcode
    !> number of atoms per fragment per code
    integer, dimension(:, :), intent(in) :: fragment_data
    !> ids of atoms belonging to a fragment
    integer, dimension(:, :, :), intent(inout) :: fragment_ids

    integer :: n_code, n_fragment, n_atom
    integer :: offset
    integer(kind=i32), dimension(:), allocatable, target :: atoms_ids

    do n_code = 1, this%num_clients
        allocate(atoms_ids(atoms_pcode(n_code)))
        call this%send_command(MCL_SEND_ATOMS_IN_FRAGMENTS, n_code)
        call mcl_receive(atoms_ids, atoms_pcode(n_code), MCL_DATA, n_code)
        offset = 0
        do n_fragment = 1, size(fragment_data, 1)
            if (fragment_data(n_fragment, n_code) == -1) then
                exit
            end if
            do n_atom = 1, fragment_data(n_fragment, n_code)
                fragment_ids(n_atom, n_fragment, n_code) = atoms_ids(offset + n_atom)
            end do
            offset = offset + fragment_data(n_fragment, n_code)
        end do
        deallocate(atoms_ids)
    end do

end subroutine gather_atom_fragment_ids

!> get the data about atoms in the system
subroutine gather_masses(this, masses, species_pcode, species_map)

    class(mimic_communicator), intent(inout) :: this
    !> mass of each atom in the system
    real(dp), dimension(:), intent(out) :: masses
    !> number of atoms per client code
    integer, dimension(:), intent(in) :: species_pcode
    !> mapping of MM species to CPMD species
    integer, dimension(:,:), intent(in) :: species_map

    real(dp), dimension(:,:), allocatable, target :: temp_masses
    integer :: n_code, n_species
    integer(kind=i32) :: num_species

    allocate(temp_masses(maxval(species_pcode), this%num_clients))
    temp_masses(:, :) = -1.0_dp

    do n_code = 1, this%num_clients
        num_species = species_pcode(n_code)

        call this%send_command(MCL_SEND_ATOM_MASSES, n_code)
        call mcl_receive(temp_masses(1:num_species, n_code), num_species, MCL_DATA, n_code)
        do n_species = 1, species_pcode(n_code)
            if (species_map(n_species, n_code) /= - 1) then
                masses(species_map(n_species, n_code)) = temp_masses(n_species, n_code)
            end if
        end do ! n_species
    end do ! n_code

end subroutine gather_masses

!> get the data about atoms in the system
subroutine gather_multipoles(this, species, multipoles, &
                              atoms_pcode, multipoles_patom)

    class(mimic_communicator), intent(inout) :: this
    !> atom species of each atom in the system (-1 for overlapped)
    integer, dimension(:,:), intent(in), target :: species
    !> multipoles of each atom in the system (0 charge for overlapped)
    real(dp), dimension(:,:,:), intent(inout), target :: multipoles
    !> number of atoms per client code
    integer, dimension(:), intent(in) :: atoms_pcode
    !> amount of multipoles per atom per code
    integer, dimension(:), intent(in) :: multipoles_patom

    integer :: offset = 1
    integer :: n_code, n_atom
    integer(kind=i32) :: n_atoms
    real(dp), dimension(:), target, allocatable :: multipoles_temp

    do n_code = 1, this%num_clients
        n_atoms = atoms_pcode(n_code)

        allocate(multipoles_temp(n_atoms * multipoles_patom(n_code)))

        call this%send_command(MCL_SEND_ATOM_MULTIPOLES, n_code)
        call mcl_receive(multipoles_temp, n_atoms * multipoles_patom(n_code), MCL_DATA, n_code)

        offset = 1

        do n_atom = 1, n_atoms
            multipoles(1:multipoles_patom(n_code), n_atom, n_code) = multipoles_temp(offset : offset + multipoles_patom(n_code) - 1)
            if (species(n_atom, n_code) == -1) then
                multipoles(:, n_atom, n_code) = 0.0_dp
            end if
            offset = offset + multipoles_patom(n_code)
        end do
        deallocate(multipoles_temp)
    end do

end subroutine gather_multipoles

!> gather coordinates of atoms within the system (dealing with overlaps)
!> and store result in the tau matrix through pointers
subroutine gather_coords(this, frag_code)

    class(mimic_communicator), intent(inout) :: this
    !> subsystem_type belonging to each code
    type(subsystem_type), dimension(:), intent(in) :: frag_code

    real(kind=dp), dimension(:), allocatable, target :: coordinates
    type(fragment_type) :: current_fragment
    integer :: code_atoms

    integer :: n_code, n_fragment, n_atom, n_map
    integer :: offset, location
    logical :: overlapped = .false.

    do n_code = 1, this%num_clients
        code_atoms = frag_code(n_code)%num_atoms

        allocate(coordinates(code_atoms * 3))
        call this%send_command(MCL_SEND_ATOM_COORDS, n_code)
        call mcl_receive(coordinates, code_atoms * 3, MCL_DATA, n_code)
        offset = 0
        do n_fragment = 1, size(frag_code(n_code)%fragments)
            current_fragment = frag_code(n_code)%fragments(n_fragment)
            do n_atom = 1, current_fragment%num_atoms
                location = (current_fragment%atoms(n_atom)%id - 1) * 3 + 1

                do n_map = 1, size(this%overlap_maps(n_code)%maps)
                    if (this%overlap_maps(n_code)%maps(n_map)%id == &
                             current_fragment%atoms(n_atom)%id) then
                        overlapped = .true.
                        exit
                    end if
                end do ! n_map

                if (overlapped) then
                    overlapped = .false.
                    ! offset = offset + 3
                    cycle
                end if

                current_fragment%atoms(n_atom)%coordinate = coordinates(location : location + 2)
            end do ! n_atom
            offset = offset + current_fragment%num_atoms * 3
        end do ! n_fragment
        deallocate(coordinates)
    end do ! n_code

end subroutine gather_coords

!> gather forces acting on atoms within the system
!> and store result in the forces matrix through pointers
subroutine gather_forces(this, frag_code)

    class(mimic_communicator), intent(inout) :: this
    !> subsystem_type belonging to each code
    type(subsystem_type), dimension(:), intent(in) :: frag_code

    real(kind=dp), dimension(:), allocatable, target :: forces
    type(fragment_type) :: current_fragment
    integer :: code_atoms

    integer :: n_code, n_fragment, n_atom
    integer :: offset

    do n_code = 1, this%num_clients
        code_atoms = frag_code(n_code)%num_atoms

        allocate(forces(code_atoms * 3))
        call this%send_command(MCL_SEND_ATOM_FORCES, n_code)
        call mcl_receive(forces, code_atoms * 3, MCL_DATA, n_code)
        ! offset = 1
        do n_fragment = 1, size(frag_code(n_code)%fragments)
            current_fragment = frag_code(n_code)%fragments(n_fragment)
            do n_atom = 1, current_fragment%num_atoms
                offset = (current_fragment%atoms(n_atom)%id - 1) * 3 + 1
                current_fragment%atoms(n_atom)%force = &
                            & current_fragment%atoms(n_atom)%force &
                            & + forces(offset : offset + 2) * frag_code(n_code)%factor
                ! offset = offset + 3
            end do ! n_atom
        end do ! n_fragment
        deallocate(forces)
    end do ! n_code

end subroutine gather_forces

!> gather MM energies from clients
subroutine gather_energies(this, energy, frag_code)

    class(mimic_communicator), intent(inout) :: this
    !> resulting total MM energy
    real(dp), intent(out) :: energy
    !> mapping of fragments per code
    type(subsystem_type), dimension(:), intent(in) :: frag_code

    integer :: n_code
    real(kind=dp), target :: temp_energy

    temp_energy = 0.0_dp
    energy = 0.0_dp
    do n_code = 1, this%num_clients
        temp_energy = 0.0_dp
        call this%send_command(MCL_SEND_ENERGY, n_code)
        call mcl_receive(temp_energy, 1, MCL_DATA, n_code)
        energy = energy + frag_code(n_code)%factor * temp_energy
    end do

end subroutine gather_energies

!> get the orders of multipoles used in each client code
subroutine get_multipole_order(this, multipoles_orders)

    class(mimic_communicator), intent(inout) :: this
    !> multipoles_orders array of orders
    integer, dimension(:), intent(out), target :: multipoles_orders

    integer :: n_client

    do n_client = 1, this%num_clients
        call this%send_command(MCL_SEND_MULTIPOLE_ORDER, n_client)
        call mcl_receive(multipoles_orders(n_client), 1, MCL_DATA, n_client)
    end do

end subroutine get_multipole_order

!> send coordinates of atoms to client codes (sends a flattened array)
subroutine send_coords(this, subsystems)

    class(mimic_communicator), intent(inout) :: this
    !> subsystem_type belonging to each code
    type(subsystem_type), dimension(:), intent(in) :: subsystems

    real(kind=dp), dimension(:), allocatable, target :: coordinates
    integer :: code_atoms

    integer :: n_code, n_atom
    integer :: offset

    do n_code = 1, this%num_clients
        code_atoms = subsystems(n_code)%num_atoms

        allocate(coordinates(code_atoms * 3))

        !$OMP PARALLEL DO PRIVATE(n_atom, offset)
        do n_atom = 1, size(subsystems(n_code)%atoms)
            offset = (subsystems(n_code)%atoms(n_atom)%id - 1) * 3 + 1
            coordinates(offset : offset + 2) = subsystems(n_code)%atoms(n_atom)%coordinate
            ! offset = offset + 3
        end do ! n_atom
        !$OMP END PARALLEL DO

        call this%send_command(MCL_RECV_ATOM_COORDS, n_code)
        call mcl_send(coordinates, code_atoms * 3, MCL_DATA, n_code)
        deallocate(coordinates)
    end do ! n_code

end subroutine send_coords

!> gather integer parameter from client codes
subroutine gather_int(this, data, command)

    class(mimic_communicator), intent(inout) :: this
    !> Integer data to gather
    integer, dimension(:), target, intent(out) :: data
    !> MCL command to execute
    integer, intent(in) :: command

    integer :: n_client

    do n_client = 1, this%num_clients
        call this%send_command(command, n_client)
        call mcl_receive(data(n_client), 1, MCL_DATA, n_client)
    end do ! n_client

end subroutine gather_int

!> issue a command to compute forces and energies
subroutine issue_force_comp(this)

    class(mimic_communicator), intent(inout) :: this

    integer :: n_client

    do n_client = 1, this%num_clients
        call this%send_command(MCL_COMPUTE_FORCE_ENERGY, n_client)
    end do ! n_client

end subroutine issue_force_comp

!> gather integer parameter from client codes
subroutine send_int(this, data, command)

    class(mimic_communicator), intent(inout) :: this
    !> Integer data to gather
    integer, dimension(:), target, intent(in) :: data
    !> MCL command to execute
    integer, intent(in) :: command

    integer :: n_client

    do n_client = 1, this%num_clients
        call this%send_command(command, n_client)
        call mcl_send(data(n_client), 1, MCL_DATA, n_client)
    end do ! n_client

end subroutine send_int

!> gather atom count from all clients
subroutine gather_atom_count(this, atoms_pcode)

    class(mimic_communicator), intent(inout) :: this
    !> Integer data to gather
    integer, dimension(:), target, intent(out) :: atoms_pcode

    call this%gather_int(atoms_pcode, MCL_SEND_NUM_ATOMS)

end subroutine gather_atom_count

!> gather fragment count from all clients
subroutine gather_fragment_count(this, frags_pcode)

    class(mimic_communicator), intent(inout) :: this
    !> Integer data to gather
    integer, dimension(:), target, intent(out) :: frags_pcode

    call this%gather_int(frags_pcode, MCL_SEND_NUM_FRAGMENTS)

end subroutine gather_fragment_count

!> gather bond count from all clients
subroutine gather_bond_count(this, bonds_pcode)

    class(mimic_communicator), intent(inout) :: this
    !> Integer data to gather
    integer, dimension(:), target, intent(out) :: bonds_pcode

    call this%gather_int(bonds_pcode, MCL_SEND_NUM_BONDS)

end subroutine gather_bond_count

!> gather angle count from all clients
subroutine gather_angle_count(this, angles_pcode)

    class(mimic_communicator), intent(inout) :: this
    !> Integer data to gather
    integer, dimension(:), target, intent(out) :: angles_pcode

    call this%gather_int(angles_pcode, MCL_SEND_NUM_ANGLES)

end subroutine gather_angle_count

!> gather atom types count from all clients
subroutine gather_type_count(this, types_pcode)

    class(mimic_communicator), intent(inout) :: this
    !> Integer data to gather
    integer, dimension(:), target, intent(out) :: types_pcode

    call this%gather_int(types_pcode, MCL_SEND_NUM_ATOM_TYPES)

end subroutine gather_type_count

!> search for unique occurences in the array
subroutine unique_entries(input, output)

    !> original array
    integer, dimension(:,:), intent(in) :: input
    !> resulting array, containing unique entries from the original one
    integer, dimension(:,:), intent(out), allocatable :: output

    integer, dimension(:,:), allocatable :: temp

    integer :: i, j, n_code
    integer :: out_size, entry_num
    logical :: found

    out_size = 0
    found = .false.
    if (allocated(output)) deallocate(output)
    allocate(output(0,0))

    do n_code = 1, size(input, 2)
        do i = 1, size(input, 1)
            if (input(i, n_code) == -1) then
                cycle
            end if

            if (out_size /= 0) then
                do j = 1, out_size
                    if (input(i, n_code) == output(j, n_code)) then
                        found = .true.
                        exit
                    end if
                end do
            end if


            if (.not. found) then
                out_size = out_size + 1
                if (out_size > size(output, 1)) then
                    entry_num = out_size
                else
                    entry_num = size(output, 1)
                end if
                allocate(temp(entry_num, size(input, 2)))
                temp(:,:) = -1
                if (size(output, 1) > 0) then
                    temp(1:size(output, 1), 1:size(input, 2)) = output
                end if
                call move_alloc(temp, output)
                output(out_size, n_code) = input(i, n_code)
            end if
            found = .false.
        end do
        out_size = 0
    end do

end subroutine unique_entries

subroutine move_alloc_str(source, dest)

    character(len=:), dimension(:, :), allocatable, &
        intent(inout) :: source
    character(len=:), dimension(:, :), allocatable, &
        intent(inout) :: dest

    integer :: i, j

    if (allocated(dest)) deallocate(dest)
    if (.not. allocated(source)) return

    allocate(character(len(source)) :: dest(size(source, 1), size(source, 2)))

    do i = 1, size(source, 1)
        do j = 1, size(source, 2)
            dest(i, j)(1:len(source)) = source (i, j)(1:len(source))
        end do
    end do

    deallocate(source)

end subroutine move_alloc_str

end module mimic_data_collect