allocates fragments for MM part
too many arguments - consider factory
break into several routines - this one does too much stuff now
Some extraodinary ugly code here - refactor ASAP
| Type | Intent | Optional | Attributes | Name | ||
|---|---|---|---|---|---|---|
| type(subsystem_type), | intent(inout), | dimension(:) | :: | subsystems | array of subsystems associated with client codes |
|
| real(kind=dp), | intent(in), | dimension(:) | :: | factors | factors applied to forces contribution (by code) Todoclarify if this is also for energies, potential, etc. |
|
| real(kind=dp), | intent(in), | dimension(:,:,:) | :: | tau | coordinate matrix of CPMD |
|
| real(kind=dp), | intent(in), | dimension(:,:,:) | :: | force | force matrix of CPMD |
|
| integer, | dimension(:) | :: | num_qm_atoms | number of quantum atoms per species |
||
| integer, | intent(in) | :: | total_species | total number of species |
||
| real(kind=dp), | intent(in), | dimension(:) | :: | cov_radii | covalent radius of each species |
|
| type(sizes_type), | intent(in) | :: | sizes | sizes data |
||
| type(system_type), | intent(inout) | :: | system_data | system data |
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