mimic_types.f90 Source File

!
!    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/>.
!

!> contains type definitions for MiMiC
!> @todo insert sanity checks (e.g. check that field
!> dimensions fit with cell parameters) @endtodo
!> @todo perhaps split into several modules
module mimic_types

    use mimic_precision, only: dp

    implicit none

    private

    public :: map_type, maps_type
    public :: sizes_type, system_type
    public :: bond_type, angle_type


    !> map for overlapping atoms, represents
    !> many-to-one relation
    type :: map_type
        private
        !> id of the atom in the code, which should handle it
        integer, public :: id
        !> ids of codes in which this atom is present
        !> (excluding the current one)
        integer, dimension(:), allocatable, public :: codes
        !> atom ids of this atom in other codes
        integer, dimension(:), allocatable, public :: atoms
    contains
        private
        procedure, public :: init => init_map
        procedure, public :: add_entry => add_map_entry
        procedure :: copy=>copy_map
        generic, public :: assignment(=) => copy
    end type map_type

    !> collection of maps for overlap mapping (for one client)
    type :: maps_type
        private
        !> id - possibly redundant
        integer, public :: id
        !> list of maps
        type(map_type), dimension(:), allocatable, public :: maps
    contains
        private
        ! procedure, public :: init => init_maps
        procedure, public :: add_entry => add_maps_entry
        procedure :: copy => copy_maps
        generic, public :: assignment(=) => copy
    end type maps_type

    type :: bond_type
        !> first atom of the bond
        integer :: atom_i
        !> second atom of the bond
        integer :: atom_j
        !> equilibrium length of the bond
        real(dp) :: length
    contains
        procedure :: copy => copy_bond
        generic, public :: assignment(=) => copy
    end type bond_type

    type :: angle_type
        !> first atom of the angle
        integer :: atom_i
        !> second atom of the angle
        integer :: atom_j
        !> third atom of the angle
        integer :: atom_k
        !> equilibrium value of the angle
        real(dp) :: angle
    contains
        procedure :: copy => copy_angle
        generic, public :: assignment(=) => copy
    end type angle_type

    !> system sizes type, used to allocate structures
    type :: sizes_type
        !> maximum number of atoms per MM species (after overlaps)
        integer :: num_atoms = 0
        !> total number of species (after overlaps)
        integer :: num_species = 0
        !> total number of atoms per code
        integer, dimension(:), allocatable :: atoms_pcode
        !> Maximal order of multipoles per code
        integer, dimension(:), allocatable :: multipoles_order
        !> number of multipole values per atoms (computed from the order)
        integer, dimension(:), allocatable :: multipoles_patom
        !> number of fragments per code
        integer, dimension(:), allocatable :: frag_num
        !> number of atoms per species
        integer, dimension(:), allocatable :: atoms_pspecies
        !> number of atoms per fragment per code
        integer, dimension(:,:), allocatable :: atoms_pfragment
        !> number of bond constraints per code
        integer, dimension(:), allocatable :: nbonds_pcode
        !> number of angle constraints per code
        integer, dimension(:), allocatable :: nangles_pcode
        !> total length of the delimeted string with atom types
        integer, dimension(:), allocatable :: types_length_pcode
        !> total number of atom types in the system
        integer, dimension(:), allocatable :: types_pcode
        !> ID of the atom to begin electrostatic treatment
        integer :: id_start
        !> ID of the atom to finish electrostatic treatment
        integer :: id_end
    end type sizes_type

    !> type holding the data about the system, used in allocation
    !> of internal structures
    type :: system_type
        !> list of atomic species per atom and per code
        integer, dimension(:, :), allocatable :: species
        !> list of atom ids per fragment per code
        integer, dimension(:,:,:), allocatable :: atom_fragment_ids
        !> list of atomic masses per atom and per code
        real(dp), dimension(:), allocatable :: masses
        !> multipole values per atom per code
        real(dp), dimension(:,:,:), allocatable :: multipole_values
        !> list of bond constraints per code
        type(bond_type), dimension(:,:), allocatable :: bonds
        !> list of angle constraints per code
        type(angle_type), dimension(:,:), allocatable :: angles
        !> deprecated: true
        !> list of atom types per code
        character(len=:), dimension(:,:), allocatable :: atom_types
        !> integer representation of element number
        integer, dimension(:), allocatable :: elements
        !> mapping between classical and CPMD species
        integer, dimension(:,:), allocatable :: species_map
    end type system_type

contains

! map_type procedures

pure subroutine init_map(this, id, codes, atoms)

    class(map_type), intent(inout) :: this
    integer, intent(in) :: id
    integer, dimension(:), intent(in) :: codes, atoms

    this%id = id
    allocate(this%atoms(size(codes)))
    allocate(this%codes(size(codes)))
    this%codes = codes
    this%atoms = atoms

end subroutine init_map

!> add an overlap entry to the map
subroutine add_map_entry(this, code_id, atom_id)

    class(map_type), intent(inout) :: this
    !> id of the code, treating the overlapped atom
    integer, intent(in) :: code_id
    !> id of an overlapped atom within its treating code
    integer, intent(in) :: atom_id

    integer, dimension(:), allocatable :: temp_codes, temp_atoms
    integer :: current_size = 0

    if (allocated(this%codes)) current_size = size(this%codes)

    allocate(temp_codes(current_size + 1))
    allocate(temp_atoms(current_size + 1))

    if (current_size /= 0) then
        temp_codes(1:current_size) = this%codes
        temp_atoms(1:current_size) = this%atoms
    end if

    call move_alloc(temp_atoms, this%atoms)
    call move_alloc(temp_codes, this%codes)

    this%codes(current_size + 1) = code_id
    this%atoms(current_size + 1) = atom_id

end subroutine add_map_entry

!> add an overlap entry to the map
subroutine copy_map(dest, source)

    class(map_type), intent(out) :: dest
    class(map_type), intent(in) :: source

    if(allocated(dest%codes)) deallocate(dest%codes)
    if(allocated(dest%atoms)) deallocate(dest%atoms)

    allocate(dest%codes, source=source%codes)
    allocate(dest%atoms, source=source%atoms)

    dest%id=source%id

end subroutine copy_map

!> add an overlap entry to the map
subroutine copy_maps(dest, source)

    class(maps_type), intent(out) :: dest
    class(maps_type), intent(in) :: source

    dest%maps = source%maps

    ! if(allocated(dest%codes)) deallocate(dest%codes)
    ! if(allocated(dest%atoms)) deallocate(dest%atoms)

    ! allocate(dest%codes, source=source%codes)
    ! allocate(dest%atoms, source=source%atoms)

    dest%id=source%id

end subroutine copy_maps

subroutine copy_bond(dest, source)
    class(bond_type), intent(out) :: dest
    class(bond_type), intent(in) :: source

    dest%atom_i = source%atom_i
    dest%atom_j = source%atom_j
    dest%length = source%length
end subroutine copy_bond

subroutine copy_angle(dest, source)
    class(angle_type), intent(out) :: dest
    class(angle_type), intent(in) :: source

    dest%atom_i = source%atom_i
    dest%atom_j = source%atom_j
    dest%atom_k = source%atom_k
    dest%angle = source%angle
end subroutine copy_angle

!> add an entry to the maps collection
subroutine add_maps_entry(this, original_atom_id, target_code_id, target_atom_id)

    class(maps_type), intent(inout) :: this
    !> id of the code which is going to treat this atom
    integer, intent(in) :: target_code_id
    !> id of an atom within the target code
    integer, intent(in) :: target_atom_id
    !> id of an atom within its original code
    integer, intent(in) :: original_atom_id

    type(map_type), dimension(:), allocatable :: temp_maps
    integer :: num_map

    if (.not. (allocated(this%maps))) then
        allocate(this%maps(1))
        this%maps(1)%id = original_atom_id
        call this%maps(1)%add_entry(target_code_id, target_atom_id)
        return
    end if

    do num_map = 1, size(this%maps)
        if (this%maps(num_map)%id == original_atom_id) then
            call this%maps(num_map)%add_entry(target_code_id, target_atom_id)
            return
        end if
    end do

    call move_alloc(this%maps, temp_maps)
    allocate(this%maps(size(temp_maps) + 1))

    this%maps(1:size(temp_maps)) = temp_maps
    do num_map = 1, size(temp_maps)
        this%maps(num_map) = temp_maps(num_map)
    end do
    this%maps(size(this%maps))%id = original_atom_id
    call this%maps(size(this%maps))%add_entry(target_code_id, target_atom_id)

end subroutine add_maps_entry

end module mimic_types