! ! 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 !> MiMiC main contains all API calls module mimic_main use mimic_cells use mimic_data_collect, only: mimic_communicator, unique_entries use mimic_errors use mimic_field_grids use mimic_fragments use mimic_long_range, only: compute_lr_energy, & compute_lr_forces, & compute_lr_potential, & compute_folded_tensors use mimic_particles use mimic_precision, only: dp use mimic_properties use mimic_short_range, only: compute_sr_energy_ions, & compute_sr_forces_ions, & compute_sr_forces_electrons, & compute_sr_potential_electrons use mimic_subsystems use mimic_types, only: system_type, sizes_type, map_type, maps_type,& bond_type, angle_type use mimic_utils, only: print_welcome, timer_start, timer_stop, timer_print implicit none !> communicator handling data collection type(mimic_communicator) :: communicator !> sort fragments according to distance between centroids integer, parameter :: CENTROID = 1 !> sort fragments according to distance between center-of-mass integer, parameter :: CENTER_OF_MASS = 2 !> sort fragments according to distance between center-of-charge for charged fragments and center-of-mass for neutral fragments integer, parameter :: CENTER_OF_CHARGE = 3 !> sort fragments according to minimum distance integer, parameter :: MINIMUM_DISTANCE = 4 !> sort atoms according to distance from quantum fragment centroid integer, parameter :: ATOM_WISE = 5 contains !> initialize MiMiC subroutine mimic_init(comm) !> MPI communicator integer, intent(inout) :: comm call communicator%init(comm) end subroutine mimic_init !> initialize error handling subroutine mimic_init_error_handler(process_id, error_handler) !> process id (preferrably MPI rank) integer, intent(in) :: process_id !> optional external error handler procedure(handle), pointer, optional :: error_handler call init_error_handling(process_id, error_handler) end subroutine mimic_init_error_handler !> MiMiC initialization routine subroutine mimic_handshake(paths) !> paths to client codes character(len=*), dimension(:), intent(in) :: paths call timer_start("mimic_handshake") call print_welcome if (size(paths) == 0) then !> @todo: may be downgraded to a warning as !> the development progresses call handle_error(SEVERITY_FATAL, & TYPE_INCORRECT_ARG, & "Empty paths array is provided", & __FILE__, __LINE__) endif call communicator%handshake(paths) call timer_stop end subroutine mimic_handshake !> set the number of clients (needed for non-root processes) !> @todo revise and check if it is better to do that with just !> a broadcast subroutine mimic_set_num_clients(num) !> number of clients integer, intent(in) :: num communicator%num_clients = num end subroutine mimic_set_num_clients !> initialize overlap maps using the input file data subroutine mimic_init_overlaps(overlaps) !> overlaps in the format: code_id : atom_id => code_id : atom_id !> (e.g. code_one : atom_three => code_two : atom_one) integer, dimension(:,:), intent(in) :: overlaps integer :: n_client, n_map type(maps_type), dimension(:), allocatable :: overlap_maps allocate(overlap_maps(communicator%num_clients)) do n_client = 1, communicator%num_clients do n_map = 1, size(overlaps, 2) if (overlaps(1, n_map) == n_client + 1) then call overlap_maps(n_client)%add_entry(overlaps(2, n_map), & overlaps(3, n_map), overlaps(4, n_map)) end if end do end do call communicator%init_overlaps(overlap_maps) end subroutine mimic_init_overlaps subroutine mimic_finalize call communicator%finalize() end subroutine mimic_finalize !> MiMiC destruction routine - deallocates all arrays and shuts down communicator subroutine mimic_destroy call communicator%destroy_comm end subroutine !> wrapper to request sizes of MM parts (number of atoms/species/fragments) subroutine mimic_request_sizes(sizes, system_data) !> sizes data !> @todo check if it can be made intent(out) type(sizes_type), intent(out) :: sizes !> atom species of each atom in the system are stored here (-1 for overlapped) !> @todo check if it can be made intent(out) type(system_type), intent(out) :: system_data integer :: n_code integer :: max_atom_count integer, dimension(:,:), allocatable, target :: ids integer, dimension(:), allocatable, target :: atoms_pcode_temp integer, dimension(:), allocatable :: na_temp call timer_start("mimic_request_sizes") sizes%num_atoms = 0 sizes%num_species = 0 allocate(sizes%atoms_pcode(communicator%num_clients)) allocate(sizes%multipoles_patom(communicator%num_clients)) allocate(sizes%frag_num(communicator%num_clients)) allocate(sizes%nbonds_pcode(communicator%num_clients)) allocate(sizes%nangles_pcode(communicator%num_clients)) allocate(sizes%types_length_pcode(communicator%num_clients)) allocate(sizes%types_pcode(communicator%num_clients)) allocate(atoms_pcode_temp(communicator%num_clients)) ! gather atom counts call communicator%gather_atom_count(atoms_pcode_temp) max_atom_count = maxval(atoms_pcode_temp) allocate(na_temp(sum(atoms_pcode_temp))) allocate(ids(max_atom_count, communicator%num_clients)) allocate(system_data%species(max_atom_count, communicator%num_clients)) allocate(sizes%multipoles_order(communicator%num_clients)) ids(:,:) = -1 system_data%species(:,:) = -1 na_temp = -1 ! gather number of atom types per code call communicator%gather_type_count(sizes%types_pcode) call communicator%get_atom_species(atoms_pcode_temp, system_data%species, & sizes%num_atoms, sizes%num_species, na_temp, & system_data%species_map) allocate(sizes%atoms_pspecies(sizes%num_species)) sizes%atoms_pspecies = na_temp(1:sizes%num_species) sizes%atoms_pcode = atoms_pcode_temp call communicator%get_multipole_order(sizes%multipoles_order) do n_code = 1, communicator%num_clients sizes%multipoles_patom(n_code) = (sizes%multipoles_order(n_code) + 3) * & (sizes%multipoles_order(n_code) + 2) * & (sizes%multipoles_order(n_code) + 1) / 6 end do ! gather fragment count call communicator%gather_fragment_count(sizes%frag_num) allocate(sizes%atoms_pfragment(maxval(sizes%frag_num), communicator%num_clients)) call communicator%gather_atom_fragment_count(sizes%frag_num, sizes%atoms_pfragment) ! gather number of bonds call communicator%gather_bond_count(sizes%nbonds_pcode) ! gather number of angles call communicator%gather_angle_count(sizes%nangles_pcode) call communicator%gather_bonds(system_data%bonds, sizes%nbonds_pcode) call communicator%gather_angles(system_data%angles, sizes%nangles_pcode) call timer_stop end subroutine mimic_request_sizes !> return data about atomic species, masses and multipoles in the MM part subroutine mimic_request_system_data(sizes, system_data) !> sizes data of the MM part type(sizes_type), intent(inout) :: sizes !> system data !> @todo intent(out)? type(system_type), intent(inout) :: system_data integer :: max_dim integer :: alloc_stat call timer_start("mimic_request_system_data") max_dim = maxval(sizes%atoms_pcode) allocate(system_data%multipole_values(maxval(sizes%multipoles_patom), & max_dim, communicator%num_clients), stat=alloc_stat) if (alloc_stat /= 0) then call handle_error(SEVERITY_FATAL, & TYPE_MEM, & "system_data%multipole_values", & __FILE__, __LINE__) endif ! allocate(system_data%species(max_dim, num_clients)) allocate(system_data%masses(sizes%num_species), stat=alloc_stat) if (alloc_stat /= 0) then call handle_error(SEVERITY_FATAL, & TYPE_MEM, & "system_data%masses", & __FILE__, __LINE__) endif allocate(system_data%elements(sizes%num_species), stat=alloc_stat) if (alloc_stat /= 0) then call handle_error(SEVERITY_FATAL, & TYPE_MEM, & "system_data%elements", & __FILE__, __LINE__) endif call communicator%gather_multipoles(system_data%species, & system_data%multipole_values, & sizes%atoms_pcode, sizes%multipoles_patom) call communicator%gather_masses(system_data%masses, & sizes%types_pcode, & system_data%species_map) allocate(system_data%atom_fragment_ids(maxval(sizes%atoms_pfragment),& maxval(sizes%frag_num),& communicator%num_clients), stat=alloc_stat) if (alloc_stat /= 0) then call handle_error(SEVERITY_FATAL, & TYPE_MEM, & "system_data%atom_fragment_ids", & __FILE__, __LINE__) endif call communicator%gather_atom_fragment_ids(sizes%atoms_pcode,& sizes%atoms_pfragment,& system_data%atom_fragment_ids) call communicator%gather_elements(sizes%types_pcode, & system_data%elements, & system_data%species_map) call timer_stop end subroutine mimic_request_system_data !> allocates fragments for MM part !> @todo too many arguments - consider factory @endtodo !> @todo break into several routines - this one does too much stuff now subroutine mimic_allocate_mm_struct(subsystems, factors, tau, force, num_qm_atoms,& total_species, cov_radii, sizes, system_data) !> array of subsystems associated with client codes type(subsystem_type), dimension(:), intent(inout) :: subsystems !> factors applied to forces contribution (by code) !> @todo clarify if this is also for energies, potential, etc. real(dp), dimension(:), intent(in) :: factors !> coordinate matrix of CPMD real(dp), dimension(:,:,:), intent(in) :: tau !> force matrix of CPMD real(dp), dimension(:,:,:), intent(in) :: force !> number of quantum atoms per species integer, dimension(:) :: num_qm_atoms !> total number of species integer, intent(in) :: total_species !> covalent radius of each species real(dp), dimension(:), intent(in) :: cov_radii !> sizes data type(sizes_type), intent(in) :: sizes !> system data type(system_type), intent(inout) :: system_data integer :: n_atoms, n_code, n_atom, n_species, n_fragment, n_map, n_bond, n_angle integer :: atom_index integer :: id_offset integer :: spid, atid, id integer :: sp_offset, atom_species integer :: target_code integer :: alloc_stat integer, dimension(communicator%num_clients) :: code_offset integer, dimension(total_species) :: species_counters integer, dimension(:,:), allocatable :: species_unique integer, dimension(:,:,:), allocatable :: at_sp_id integer, dimension(:), allocatable :: fragment_atoms real(dp) :: charge real(dp), dimension(3) :: origin logical :: found logical :: overlapped type(atom_type), dimension(:), allocatable :: subsystem_atoms type(multipoles_type), allocatable :: multipoles type(map_type) :: map type(bond_type), dimension(:,:), allocatable :: temp_bonds type(angle_type), dimension(:,:), allocatable :: temp_angles call timer_start("mimic_allocate_mm_struct") origin(:) = 0 allocate(temp_bonds(maxval(sizes%nbonds_pcode), communicator%num_clients), stat=alloc_stat) if (alloc_stat /= 0) then call handle_error(SEVERITY_FATAL, & TYPE_MEM, & "temp_bonds", & __FILE__, __LINE__) endif allocate(temp_angles(maxval(sizes%nangles_pcode), communicator%num_clients), stat=alloc_stat) if (alloc_stat /= 0) then call handle_error(SEVERITY_FATAL, & TYPE_MEM, & "temp_angles", & __FILE__, __LINE__) endif if (allocated(system_data%bonds)) then temp_bonds(:,:) = system_data%bonds(:,:) end if if (allocated(system_data%angles)) then temp_angles(:,:) = system_data%angles(:,:) end if call unique_entries(system_data%species, species_unique) species_counters(:) = 1 if (size(species_unique) > 0) then sp_offset = size(num_qm_atoms) do n_code = 1, communicator%num_clients code_offset(n_code) = sp_offset sp_offset = sp_offset + size(pack(species_unique(:, n_code), mask = (species_unique(:, n_code) /= -1))) end do endif allocate(at_sp_id(2, maxval(sizes%atoms_pcode, 1), communicator%num_clients + 1), stat=alloc_stat) if (alloc_stat /= 0) then call handle_error(SEVERITY_FATAL, & TYPE_MEM, & "at_sp_id", & __FILE__, __LINE__) endif at_sp_id(:,:,:) = 0 id_offset = 1 do n_species = 1, size(num_qm_atoms) do n_atom = 1, num_qm_atoms(n_species) at_sp_id(2, id_offset, 1) = n_atom at_sp_id(1, id_offset, 1) = n_species id_offset = id_offset + 1 end do end do do n_code = 1, communicator%num_clients do n_atom = 1, sizes%atoms_pcode(n_code) found = .false. atom_species = system_data%species(n_atom, n_code) if (atom_species == -1) then cycle end if do n_species = 1, size(species_unique(:, n_code)) if (species_unique(n_species, n_code) /= -1 & .and. atom_species == species_unique(n_species, n_code)) then found = .true. spid = system_data%species_map(atom_species, n_code) + code_offset(1) atid = species_counters(spid) species_counters(spid) = species_counters(spid) + 1 at_sp_id(1, n_atom, n_code + 1) = spid at_sp_id(2, n_atom, n_code + 1) = atid exit end if end do ! n_species if (.not. found) then call handle_error(SEVERITY_FATAL, & TYPE_INCORRECT_ARG, & "Unrecognized atom type is encountered", & __FILE__, __LINE__) endif end do ! n_atom end do ! n_code do n_code = 1, communicator%num_clients allocate(fragment_atoms(0:sizes%frag_num(n_code)), stat=alloc_stat) if (alloc_stat /= 0) then call handle_error(SEVERITY_FATAL, & TYPE_MEM, & "fragment_atoms", & __FILE__, __LINE__) endif allocate(subsystem_atoms(sizes%atoms_pcode(n_code)), stat=alloc_stat) if (alloc_stat /= 0) then call handle_error(SEVERITY_FATAL, & TYPE_MEM, & "subsystem_atoms", & __FILE__, __LINE__) endif fragment_atoms(0) = 0 atom_index = 1 id_offset = 0 do n_fragment = 1, sizes%frag_num(n_code) n_atoms = sizes%atoms_pfragment(n_fragment, n_code) spid = -1 do n_atom = 1, n_atoms id = system_data%atom_fragment_ids(n_atom, n_fragment, n_code) atom_species = system_data%species(id, n_code) overlapped = .false. if (atom_species == -1) then do n_map = 1, size(communicator%overlap_maps(n_code)%maps) if (id == communicator%overlap_maps(n_code)%maps(n_map)%id) then target_code = communicator%overlap_maps(n_code)%maps(n_map)%codes(1) map = communicator%overlap_maps(n_code)%maps(n_map) spid = at_sp_id(1, map%atoms(1), target_code) atid = at_sp_id(2, map%atoms(1), target_code) overlapped = .true. end if end do ! n_map else spid = at_sp_id(1, id, n_code + 1) atid = at_sp_id(2, id, n_code + 1) end if charge = system_data%multipole_values(1, id, n_code) allocate(multipoles) call multipoles%init(n_atom, origin, sizes%multipoles_order(n_code), system_data%multipole_values(1:sizes%multipoles_patom(n_code), id_offset + n_atom, n_code)) call subsystem_atoms(atom_index)%init(system_data%atom_fragment_ids(n_atom, n_fragment, n_code), & spid, n_atom, charge, cov_radii(spid), & tau(:, atid, spid), force(:, atid, spid), multipoles, overlapped) deallocate(multipoles) fragment_atoms(n_fragment) = atom_index atom_index = atom_index + 1 !! @todo Some extraodinary ugly code here - refactor ASAP do n_bond = 1, sizes%nbonds_pcode(n_code) if (system_data%bonds(n_bond, n_code)%atom_i == id) then if (spid > size(num_qm_atoms)) then temp_bonds(n_bond, n_code)%atom_i = sum(num_qm_atoms) + & sum(sizes%atoms_pspecies(1:spid - size(num_qm_atoms) - 1), n_code) + & atid else temp_bonds(n_bond, n_code)%atom_i = sum(num_qm_atoms(1:spid - 1)) + atid end if end if if (system_data%bonds(n_bond, n_code)%atom_j == id) then if (spid > size(num_qm_atoms)) then temp_bonds(n_bond, n_code)%atom_j = sum(num_qm_atoms) + & sum(sizes%atoms_pspecies(1:spid - size(num_qm_atoms) - 1), n_code) + & atid else temp_bonds(n_bond, n_code)%atom_j = sum(num_qm_atoms(1:spid - 1)) + atid end if end if end do !n_bond do n_angle = 1, sizes%nangles_pcode(n_code) if (system_data%angles(n_angle, n_code)%atom_i == id) then if (spid > size(num_qm_atoms)) then temp_angles(n_angle, n_code)%atom_i = sum(num_qm_atoms) + & sum(sizes%atoms_pspecies(1:spid - size(num_qm_atoms) - 1), n_code) + & atid else temp_angles(n_angle, n_code)%atom_i = sum(num_qm_atoms(1:spid - 1)) + atid end if end if if (system_data%angles(n_angle, n_code)%atom_j == id) then if (spid > size(num_qm_atoms)) then temp_angles(n_angle, n_code)%atom_j = sum(num_qm_atoms) + & sum(sizes%atoms_pspecies(1:spid - size(num_qm_atoms) - 1), n_code) + & atid else temp_angles(n_angle, n_code)%atom_j = sum(num_qm_atoms(1:spid - 1)) + atid end if end if if (system_data%angles(n_angle, n_code)%atom_k == id) then if (spid > size(num_qm_atoms)) then temp_angles(n_angle, n_code)%atom_k = sum(num_qm_atoms) + & sum(sizes%atoms_pspecies(1:spid - size(num_qm_atoms) - 1), n_code) + & atid else temp_angles(n_angle, n_code)%atom_k = sum(num_qm_atoms(1:spid - 1)) + atid end if end if end do ! n_angle end do ! n_atom id_offset = id_offset + n_atoms end do ! n_fragment call subsystems(n_code)%init(n_code, factors(n_code), subsystem_atoms) call subsystems(n_code)%define_fragments(sizes%frag_num(n_code), fragment_atoms) deallocate(fragment_atoms) deallocate(subsystem_atoms) end do ! n_client if (maxval(sizes%nbonds_pcode) > 0) then system_data%bonds = temp_bonds end if if (maxval(sizes%nangles_pcode) > 0) then system_data%angles = temp_angles end if call timer_stop end subroutine mimic_allocate_mm_struct !> allocates quantum fragment (at the moment one per simulation) !> @todo too many arguments - consider factory subroutine mimic_allocate_qm_struct(quantum_fragment, n_sp_atoms, charges, origin, & box, n_points, n_points_r, rho, extf, tau, force) !> quantum subsystem type(quantum_fragment_type), intent(inout) :: quantum_fragment !> number of atoms per quantum species integer, dimension(:), intent(in) :: n_sp_atoms !> number of points in the electronic grid in each direction integer, dimension(3), intent(in) :: n_points !> real number of points in the electronic grid in each direction integer, dimension(3), intent(in) :: n_points_r !> ionic charges real(dp), dimension(:), intent(in) :: charges !> point of origin of quantum box real(dp), dimension(3), intent(in) :: origin !> quantum box vectors real(dp), dimension(3,3), intent(in) :: box !> electronic density matrix real(dp), dimension(:,:,:), intent(in) :: rho !> external field on the electronic grid real(dp), dimension(:,:,:), intent(in) :: extf !> coordinate matrix of CPMD real(dp), dimension(:,:,:), intent(in) :: tau !> force matrix of CPMD real(dp), dimension(:,:,:), intent(in) :: force integer :: n_species integer :: species, atom, offset integer :: alloc_stat type(nucleus_type), dimension(:), allocatable :: nuclei type(density_type) :: density type(potential_type) :: potential type(cell_type) :: cell integer :: nuclei_count call timer_start("mimic_allocate_qm_struct") nuclei_count = 0 n_species = size(n_sp_atoms) do species = 1, n_species nuclei_count = nuclei_count + n_sp_atoms(species) end do allocate(nuclei(nuclei_count), stat=alloc_stat) if (alloc_stat /= 0) then call handle_error(SEVERITY_FATAL, & TYPE_MEM, & "nuclei", & __FILE__, __LINE__) endif offset = 0 do species = 1, n_species do atom = 1, n_sp_atoms(species) offset = offset + 1 call nuclei(offset)%init(offset, species, atom, charges(species), & tau(:, atom, species), force(:, atom, species), .false.) end do ! atom end do ! species call cell%init(1, n_points, n_points_r, origin, box) call density%init(1, rho) call potential%init(1, extf) call quantum_fragment%init(1, nuclei, cell, density, potential) call timer_stop end subroutine mimic_allocate_qm_struct !> wrapper to request coordinates of MM part subroutine mimic_collect_coordinates(subsystems) !> array of subsystems associated with client codes type(subsystem_type), dimension(:) :: subsystems call timer_start("mimic_collect_coordinates") call communicator%gather_coords(subsystems) call timer_stop end subroutine mimic_collect_coordinates subroutine mimic_send_coords(subsystems) !> array of subsystems associated with client codes type(subsystem_type), dimension(:) :: subsystems call timer_start("mimic_send_coords") call communicator%send_coords(subsystems) call communicator%issue_force_comp() call timer_stop end subroutine mimic_send_coords !> wrapper to request MM forces !> @todo separate forces and potential computation and remove logics subroutine mimic_collect_forces(subsystems) !> array of subsystems associated with client codes type(subsystem_type), dimension(:), intent(inout) :: subsystems integer :: alloc_stat call timer_start("mimic_collect_forces") call communicator%gather_forces(subsystems) call timer_stop end subroutine mimic_collect_forces !> wrapper to request MM energies subroutine mimic_collect_energies(subsystems, energy) !> array of subsystems associated with client codes type(subsystem_type), dimension(:), intent(inout) :: subsystems !> energy from external programs real(dp), intent(out) :: energy call timer_start("mimic_collect_energies") call communicator%gather_energies(energy, subsystems) call timer_stop end subroutine mimic_collect_energies !> sort fragments into short- and long-range according to distance criterion subroutine mimic_sort_fragments(subsystem, quantum_fragment, cutoff_distance, sorting_method) intrinsic :: norm2 !> array of subsystems associated with client codes type(subsystem_type), intent(inout) :: subsystem !> the quantum fragment type(quantum_fragment_type), intent(in) :: quantum_fragment !> distance cutoff real(dp), intent(in) :: cutoff_distance !> sorting method (centroid, center_of_mass, center_of_charge, minimum_distance) integer, intent(in) :: sorting_method integer :: i, j, k integer :: num_atoms integer, dimension(:), allocatable :: sr_indices integer, dimension(:), allocatable :: lr_indices integer, dimension(:), allocatable :: temporary_indices real(dp) :: distance real(dp), dimension(3) :: fragment_center real(dp), dimension(3) :: qm_fragment_center call timer_start("mimic_sort_fragments") if (allocated(subsystem%sr_atoms)) then deallocate(subsystem%sr_atoms) end if if (allocated(subsystem%lr_atoms)) then deallocate(subsystem%lr_atoms) end if if (allocated(subsystem%sr_fragments)) then deallocate(subsystem%sr_fragments) end if if (allocated(subsystem%lr_fragments)) then deallocate(subsystem%lr_fragments) end if allocate(sr_indices(0)) allocate(lr_indices(0)) select case (sorting_method) ! sort according to distance between centroids case (CENTROID) qm_fragment_center = quantum_fragment%centroid() do i = 1, subsystem%num_fragments fragment_center = subsystem%fragments(i)%centroid() distance = norm2(fragment_center - qm_fragment_center) if (distance <= cutoff_distance) then call move_alloc(sr_indices, temporary_indices) allocate(sr_indices, source=[temporary_indices, i]) else call move_alloc(lr_indices, temporary_indices) allocate(lr_indices, source=[temporary_indices, i]) end if end do ! sort according to distance between center-of-mass case (CENTER_OF_MASS) qm_fragment_center = quantum_fragment%center_of_mass() do i = 1, subsystem%num_fragments fragment_center = subsystem%fragments(i)%center_of_mass() distance = norm2(fragment_center - qm_fragment_center) if (distance <= cutoff_distance) then call move_alloc(sr_indices, temporary_indices) allocate(sr_indices, source=[temporary_indices, i]) else call move_alloc(lr_indices, temporary_indices) allocate(lr_indices, source=[temporary_indices, i]) end if end do ! sort according to distance between center-of-charge for charged fragments and center-of-mass for other fragments case (CENTER_OF_CHARGE) qm_fragment_center = quantum_fragment%center_of_charge() do i = 1, subsystem%num_fragments if (nint(subsystem%fragments(i)%charge) /= 0) then fragment_center = subsystem%fragments(i)%center_of_charge() else fragment_center = subsystem%fragments(i)%center_of_mass() end if distance = norm2(fragment_center - qm_fragment_center) if (distance <= cutoff_distance) then call move_alloc(sr_indices, temporary_indices) allocate(sr_indices, source=[temporary_indices, i]) else call move_alloc(lr_indices, temporary_indices) allocate(lr_indices, source=[temporary_indices, i]) end if end do ! sort according to minimum distance case (MINIMUM_DISTANCE) do i = 1, subsystem%num_fragments distance = compute_minimum_distance(subsystem%fragments(i)%atoms, quantum_fragment%nuclei) if (distance <= cutoff_distance) then call move_alloc(sr_indices, temporary_indices) allocate(sr_indices, source=[temporary_indices, i]) else call move_alloc(lr_indices, temporary_indices) allocate(lr_indices, source=[temporary_indices, i]) end if end do ! sort individual atoms into one short-range and one long-range fragment case (ATOM_WISE) qm_fragment_center = quantum_fragment%centroid() do i = 1, subsystem%num_atoms distance = norm2(qm_fragment_center - subsystem%atoms(i)%coordinate) if (distance <= cutoff_distance) then call move_alloc(sr_indices, temporary_indices) allocate(sr_indices, source=[temporary_indices, i]) else call move_alloc(lr_indices, temporary_indices) allocate(lr_indices, source=[temporary_indices, i]) end if end do case default call handle_error(SEVERITY_FATAL, & TYPE_INCORRECT_ARG, & "Unrecognized atom sorting method", & __FILE__, __LINE__) end select if (sorting_method == ATOM_WISE) then subsystem%num_sr_fragments = 1 allocate(subsystem%sr_fragments(subsystem%num_sr_fragments)) subsystem%num_sr_atoms = size(sr_indices) allocate(subsystem%sr_atoms(size(sr_indices))) do i = 1, subsystem%num_sr_atoms subsystem%sr_atoms(i) = subsystem%atoms(sr_indices(i)) end do call subsystem%sr_fragments(1)%init(1, subsystem%sr_atoms) subsystem%num_lr_fragments = 1 allocate(subsystem%lr_fragments(subsystem%num_lr_fragments)) subsystem%num_lr_atoms = size(lr_indices) allocate(subsystem%lr_atoms(size(lr_indices))) do i = 1, subsystem%num_lr_atoms subsystem%lr_atoms(i) = subsystem%atoms(lr_indices(i)) end do call subsystem%lr_fragments(1)%init(1, subsystem%lr_atoms) else subsystem%num_sr_fragments = size(sr_indices) allocate(subsystem%sr_fragments(subsystem%num_sr_fragments)) num_atoms = 0 do i = 1, subsystem%num_sr_fragments subsystem%sr_fragments(i) = subsystem%fragments(sr_indices(i)) num_atoms = num_atoms + subsystem%sr_fragments(i)%num_atoms end do subsystem%num_sr_atoms = num_atoms allocate(subsystem%sr_atoms(subsystem%num_sr_atoms)) k = 0 do i = 1, subsystem%num_sr_fragments j = k + 1 k = k + subsystem%sr_fragments(i)%num_atoms subsystem%sr_atoms(j:k) = subsystem%sr_fragments(i)%atoms end do subsystem%num_lr_fragments = size(lr_indices) allocate(subsystem%lr_fragments(subsystem%num_lr_fragments)) num_atoms = 0 do i = 1, subsystem%num_lr_fragments subsystem%lr_fragments(i) = subsystem%fragments(lr_indices(i)) num_atoms = num_atoms + subsystem%lr_fragments(i)%num_atoms end do subsystem%num_lr_atoms = num_atoms allocate(subsystem%lr_atoms(subsystem%num_lr_atoms)) k = 0 do i = 1, subsystem%num_lr_fragments j = k + 1 k = k + subsystem%lr_fragments(i)%num_atoms subsystem%lr_atoms(j:k) = subsystem%lr_fragments(i)%atoms end do end if call timer_stop end subroutine mimic_sort_fragments !> compute folded tensors for long-range electrostatics !> @todo currently works only for one client code subroutine mimic_compute_tensors(subsystems, & quantum_fragment, & tensors, & tensor_sums, & atom_start, & atom_end) !> array of subsystems associated with client codes type(subsystem_type), dimension(:), intent(inout) :: subsystems !> quantum subsystem type(quantum_fragment_type), intent(inout) :: quantum_fragment !> array with tensors real(dp), dimension(:,:), allocatable, intent(out) :: tensors !> array with tensor sums real(dp), dimension(:,:), allocatable, intent(out) :: tensor_sums !> first atom id in the chunk integer, dimension(:), intent(in) :: atom_start !> last atom id in the chunk integer, dimension(:), intent(in) :: atom_end real(dp), dimension(3) :: multipole_origin integer :: multipole_order integer :: i call timer_start("mimic_compute_tensors") multipole_origin = quantum_fragment%nuclear_multipoles%origin multipole_order = quantum_fragment%nuclear_multipoles%order allocate(tensors(size(quantum_fragment%nuclear_multipoles%values), & size(subsystems))) allocate(tensor_sums(size(quantum_fragment%nuclear_multipoles%values), & size(subsystems))) do i = 1, size(subsystems) call compute_folded_tensors(tensors(:,i), tensor_sums(:,i), & multipole_origin, multipole_order, & subsystems(i)%lr_atoms(atom_start(i):atom_end(i))) end do call timer_stop end subroutine mimic_compute_tensors !> compute interaction energy !> @todo currently works only for one client code subroutine mimic_compute_energy(subsystems, & quantum_fragment, & energy, & sr_atom_start, & sr_atom_end, & tensor_sums) !> array of subsystems associated with client codes type(subsystem_type), dimension(:), intent(inout) :: subsystems !> quantum subsystem type(quantum_fragment_type), intent(inout) :: quantum_fragment !> interaction energy real(dp), intent(out) :: energy !> index of the first short-range atom for this process integer, dimension(:), intent(in) :: sr_atom_start !> index of the last short-range atom for this process integer, dimension(:), intent(in) :: sr_atom_end !> Array stroing tensor sums real(dp), dimension(:,:), intent(in), allocatable :: tensor_sums integer :: i real(dp) :: temp_energy call timer_start("mimic_compute_energy") energy = 0.0_dp do i = 1, size(subsystems) temp_energy = 0.0_dp call compute_sr_energy_ions(subsystems(i)%sr_atoms, & quantum_fragment%nuclei, & temp_energy, & sr_atom_start(i), & sr_atom_end(i)) energy = energy + temp_energy if (subsystems(i)%num_lr_atoms > 0) then temp_energy = 0.0_dp call compute_lr_energy(quantum_fragment, & temp_energy, & tensor_sums(:,i)) energy = energy + temp_energy end if end do call timer_stop end subroutine mimic_compute_energy !> compute forces from electrostatic interactions !> @todo fix for multiple subsystems subroutine mimic_compute_forces(subsystems, & quantum_fragment, & x_start, & x_end, & gr_sr_atom_start, & gr_sr_atom_end, & gr_lr_atom_start, & gr_lr_atom_end, & pr_sr_atom_start, & pr_sr_atom_end, & pr_lr_atom_start, & pr_lr_atom_end, & root) !> array of subsystems associated with client codes type(subsystem_type), dimension(:), intent(inout) :: subsystems !> quantum subsystem type(quantum_fragment_type), intent(inout) :: quantum_fragment !> index of the first yz-plane for this process integer, intent(in) :: x_start !> index of the last yz-plane for this process integer, intent(in) :: x_end !> index of the first short-range atom for this process group integer, dimension(:), intent(in) :: gr_sr_atom_start !> index of the last short-range atom for this process group integer, dimension(:), intent(in) :: gr_sr_atom_end !> index of the first long-ange atom for this process group integer, dimension(:), intent(in) :: gr_lr_atom_start !> index of the last long-range atom for this process group integer, dimension(:), intent(in) :: gr_lr_atom_end !> index of the first short-range atom for this process integer, dimension(:), intent(in) :: pr_sr_atom_start !> index of the last short-range atom for this process integer, dimension(:), intent(in) :: pr_sr_atom_end !> index of the first long-ange atom for this process integer, dimension(:), intent(in) :: pr_lr_atom_start !> index of the last long-range atom for this process integer, dimension(:), intent(in) :: pr_lr_atom_end !> logical for root process logical, intent(in) :: root integer :: i call timer_start("mimic_compute_forces") do i = 1, size(subsystems) call compute_sr_forces_electrons(subsystems(i)%sr_atoms, & quantum_fragment%density%field, & quantum_fragment%cell, & x_start, x_end, & gr_sr_atom_start(i), & gr_sr_atom_end(i)) call compute_sr_forces_ions(subsystems(i)%sr_atoms, & quantum_fragment%nuclei, & pr_sr_atom_start(i), & pr_sr_atom_end(i)) if (subsystems(i)%num_lr_atoms > 0) then call compute_lr_forces(subsystems(i)%lr_atoms, & quantum_fragment, & pr_lr_atom_start(i), & pr_lr_atom_end(i)) end if end do call timer_stop end subroutine mimic_compute_forces !> compute electrostatic potential !> @todo fix for multiple subsystems subroutine mimic_compute_potential(subsystems, & quantum_fragment, & x_start, x_end, & sr_atom_start, sr_atom_end, & tensor_sums) !> array of subsystems associated with client codes type(subsystem_type), dimension(:), intent(inout) :: subsystems !> quantum subsystem type(quantum_fragment_type), intent(inout) :: quantum_fragment !> index of the first yz-plane for this process integer, intent(in) :: x_start !> index of the last yz-plane for this process integer, intent(in) :: x_end !> index of the first short-range atom for this process integer, dimension(:), intent(in) :: sr_atom_start !> index of the last short-range atom for this process integer, dimension(:), intent(in) :: sr_atom_end !> Array stroing tensor sums real(dp), dimension(:,:), intent(inout), allocatable :: tensor_sums integer :: i call timer_start("mimic_compute_potential") quantum_fragment%potential%field = 0.0_dp do i = 1, size(subsystems) call compute_sr_potential_electrons(subsystems(i)%sr_atoms, & quantum_fragment%potential%field, & quantum_fragment%cell, & x_start, x_end, & sr_atom_start(i), & sr_atom_end(i)) if (subsystems(i)%num_lr_atoms > 0) then call compute_lr_potential(quantum_fragment, & x_start, x_end, & tensor_sums(:,i)) end if end do call timer_stop end subroutine mimic_compute_potential !> compute multipoles of the quantum fragment subroutine mimic_compute_multipoles(quantum_fragment, order, x_start, x_end) !> quantum subsystem type(quantum_fragment_type), intent(inout) :: quantum_fragment !> order of the multipole expansion integer, intent(in) :: order !> index of the first yz-plane for this process integer, intent(in) :: x_start !> index of the last yz-plane for this process integer, intent(in) :: x_end real(dp), dimension(3) :: multipole_origin call timer_start("mimic_compute_multipoles") multipole_origin = compute_centroid(quantum_fragment%nuclei) call quantum_fragment%compute_nuclear_multipoles(multipole_origin, order) call quantum_fragment%compute_electronic_multipoles(multipole_origin, order, & x_start, x_end) call timer_stop end subroutine mimic_compute_multipoles !> apply MIC to the whole system with respect to a pivot point !> @todo rename because this is not MIC subroutine mimic_min_image(box, pivot, subsystems) !> vectors of the encompassing box real(dp), dimension(3,3), intent(in) :: box !> pivot point real(dp), dimension(3), intent(in) :: pivot !> classical fragments class(subsystem_type), dimension(:), intent(inout) :: subsystems real(dp), dimension(3) :: frag_center, orig_center real(dp), dimension(3) :: r_ij, stmp, s_ij real(dp), dimension(3,3) :: inv_box integer :: i, j, k type(cell_type) :: classical_box call timer_start("mimic_min_image") call classical_box%init(0, [1, 1, 1], [1, 1, 1], [0.0_dp, 0.0_dp, 0.0_dp], box) inv_box = classical_box%lattice_inverse do k = 1, size(subsystems) !$OMP PARALLEL DO PRIVATE(j, frag_center, orig_center, r_ij, stmp, & !$OMP& s_ij, i) do j = 1, subsystems(k)%num_fragments frag_center = subsystems(k)%fragments(j)%centroid() orig_center = frag_center r_ij = frag_center - pivot stmp = matmul(inv_box, r_ij) s_ij = stmp - anint(stmp) frag_center = matmul(s_ij, box) + pivot s_ij = frag_center - orig_center do i = 1, subsystems(k)%fragments(j)%num_atoms if (.not. subsystems(k)%fragments(j)%atoms(i)%overlapped) then subsystems(k)%fragments(j)%atoms(i)%coordinate = & subsystems(k)%fragments(j)%atoms(i)%coordinate + s_ij endif end do end do !$OMP END PARALLEL DO end do call timer_stop end subroutine mimic_min_image !> re-center the quantum box subroutine mimic_center_qm(r_diff, subsystem, quantum_fragment) !> resulting translation vector real(dp), dimension(3), intent(out) :: r_diff !> classical subsystem class(subsystem_type), dimension(:), intent(in) :: subsystem !> quantum fragment class(quantum_fragment_type), intent(in) :: quantum_fragment real(dp) :: max_x, max_y, max_z real(dp) :: min_x, min_y, min_z real(dp), dimension(3,3) :: quantum_box integer :: n_client, n_atom, n_qatom call timer_start("mimic_center_qm") max_x = -huge(0.0_dp) max_y = -huge(0.0_dp) max_z = -huge(0.0_dp) min_x = huge(0.0_dp) min_y = huge(0.0_dp) min_z = huge(0.0_dp) quantum_box = quantum_fragment%cell%lattice do n_qatom = 1, quantum_fragment%num_nuclei max_x = max(max_x, quantum_fragment%nuclei(n_qatom)%coordinate(1)) max_y = max(max_y, quantum_fragment%nuclei(n_qatom)%coordinate(2)) max_z = max(max_z, quantum_fragment%nuclei(n_qatom)%coordinate(3)) min_x = min(min_x, quantum_fragment%nuclei(n_qatom)%coordinate(1)) min_y = min(min_y, quantum_fragment%nuclei(n_qatom)%coordinate(2)) min_z = min(min_z, quantum_fragment%nuclei(n_qatom)%coordinate(3)) end do r_diff(1) = (quantum_box(1,1) - max_x - min_x) * 0.5_dp r_diff(2) = (quantum_box(2,2) - max_y - min_y) * 0.5_dp r_diff(3) = (quantum_box(3,3) - max_z - min_z) * 0.5_dp !$OMP PARALLEL DO PRIVATE(n_qatom) do n_qatom = 1, quantum_fragment%num_nuclei quantum_fragment%nuclei(n_qatom)%coordinate = & quantum_fragment%nuclei(n_qatom)%coordinate + r_diff end do !$OMP END PARALLEL DO do n_client = 1, size(subsystem) !$OMP PARALLEL DO PRIVATE(n_atom) do n_atom = 1, size(subsystem(n_client)%atoms) if (.not. subsystem(n_client)%atoms(n_atom)%overlapped) then subsystem(n_client)%atoms(n_atom)%coordinate = & subsystem(n_client)%atoms(n_atom)%coordinate + r_diff endif end do !$OMP END PARALLEL DO end do call timer_stop end subroutine mimic_center_qm !> Routine to count the number of interdependent constraints amoung bonds subroutine mimic_count_trcnst(nconn, conn_map, bonds, bonds_pcode) !> number of dependencies per bond integer, dimension(:), allocatable, intent(out) :: nconn !> dependency map integer, dimension(:, :), allocatable, intent(out) :: conn_map !> list of bond constraints type(bond_type), dimension(:,:), intent(in) :: bonds !> number of bonds per code integer, dimension(:), intent(in) :: bonds_pcode integer :: n_code, n_bond, n_bond2 integer :: current_maxconn = 0, current_conn = 0, id = 0, id2 = 0 integer, dimension(:, :), allocatable :: temp_conn_map allocate(nconn(sum(bonds_pcode))) do n_code = 1, size(bonds_pcode) do n_bond = 1, bonds_pcode(n_code) current_conn = 0 id = id + 1 id2 = sum(bonds_pcode(1:n_code - 1)) do n_bond2 = 1, bonds_pcode(n_code) id2 = id2 + 1 if (bonds(n_bond, n_code)%atom_i == bonds(n_bond2, n_code)%atom_i .or. & bonds(n_bond, n_code)%atom_j == bonds(n_bond2, n_code)%atom_j .or. & bonds(n_bond, n_code)%atom_i == bonds(n_bond2, n_code)%atom_j .or. & bonds(n_bond, n_code)%atom_j == bonds(n_bond2, n_code)%atom_i) then current_conn = current_conn + 1 if (current_conn > current_maxconn) then allocate(temp_conn_map(current_conn, sum(bonds_pcode))) if (allocated(conn_map)) then temp_conn_map(1:current_maxconn, :) = conn_map deallocate(conn_map) end if current_maxconn = current_conn call move_alloc(temp_conn_map, conn_map) end if conn_map(current_conn, id) = id2 end if end do ! n_bond2 nconn(id) = current_conn end do ! n_bond end do ! n_code end subroutine mimic_count_trcnst !> Routine that does translation of the whole system with a given vector subroutine mimic_translate(quantum_fragment, subsystems, shift) type(quantum_fragment_type), intent(inout) :: quantum_fragment type(subsystem_type), dimension(:), intent(inout) :: subsystems !> Vector of translation to apply real(dp), dimension(3) :: shift integer :: i_c, i_a call timer_start("mimic_translate") do i_c = 1, size(subsystems) !$OMP PARALLEL DO PRIVATE(i_a) do i_a = 1, size(subsystems(i_c)%atoms) if (.not. subsystems(i_c)%atoms(i_a)%overlapped) then subsystems(i_c)%atoms(i_a)%coordinate = subsystems(i_c)%atoms(i_a)%coordinate + shift end if end do ! i_a !$OMP END PARALLEL DO end do ! i_c !$OMP PARALLEL DO PRIVATE(i_a) do i_a = 1, quantum_fragment%num_nuclei quantum_fragment%nuclei(i_a)%coordinate = quantum_fragment%nuclei(i_a)%coordinate + shift end do !$OMP END PARALLEL DO call timer_stop end subroutine mimic_translate !> Routine that prints timings accumulated during the current run of MiMiC subroutine mimic_print_timings call timer_print end subroutine mimic_print_timings end module mimic_main