1!--------------------------------------------------------------------------------------------------! 2! CP2K: A general program to perform molecular dynamics simulations ! 3! Copyright (C) 2000 - 2020 CP2K developers group ! 4!--------------------------------------------------------------------------------------------------! 5 6! ************************************************************************************************** 7!> \brief Set of routines handling the localization for molecular properties 8! ************************************************************************************************** 9MODULE qs_loc_molecules 10 USE cell_types, ONLY: pbc 11 USE cp_log_handling, ONLY: cp_get_default_logger,& 12 cp_logger_type 13 USE cp_para_types, ONLY: cp_para_env_type 14 USE distribution_1d_types, ONLY: distribution_1d_type 15 USE kinds, ONLY: dp 16 USE memory_utilities, ONLY: reallocate 17 USE message_passing, ONLY: mp_max,& 18 mp_minloc 19 USE molecule_kind_types, ONLY: get_molecule_kind,& 20 molecule_kind_type 21 USE molecule_types, ONLY: molecule_type 22 USE particle_types, ONLY: particle_type 23 USE qs_loc_types, ONLY: qs_loc_env_new_type 24#include "./base/base_uses.f90" 25 26 IMPLICIT NONE 27 28 PRIVATE 29 30 ! *** Public *** 31 PUBLIC :: wfc_to_molecule 32 33 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_loc_molecules' 34 35CONTAINS 36 37! ************************************************************************************************** 38!> \brief maps wfc's to molecules and also prints molecular dipoles 39!> \param qs_loc_env ... 40!> \param center ... 41!> \param molecule_set ... 42!> \param ispin ... 43!> \param nspins ... 44! ************************************************************************************************** 45 SUBROUTINE wfc_to_molecule(qs_loc_env, center, molecule_set, ispin, nspins) 46 TYPE(qs_loc_env_new_type), INTENT(IN) :: qs_loc_env 47 REAL(KIND=dp), INTENT(IN) :: center(:, :) 48 TYPE(molecule_type), POINTER :: molecule_set(:) 49 INTEGER, INTENT(IN) :: ispin, nspins 50 51 CHARACTER(len=*), PARAMETER :: routineN = 'wfc_to_molecule', & 52 routineP = moduleN//':'//routineN 53 54 INTEGER :: counter, first_atom, i, iatom, ikind, imol, imol_now, istate, k, local_location, & 55 natom, natom_loc, natom_max, nkind, nmol, nstate 56 INTEGER, POINTER :: wfc_to_atom_map(:) 57 REAL(KIND=dp) :: dr(3), mydist(2), ria(3) 58 REAL(KIND=dp), POINTER :: distance(:), r(:, :) 59 TYPE(cp_logger_type), POINTER :: logger 60 TYPE(cp_para_env_type), POINTER :: para_env 61 TYPE(distribution_1d_type), POINTER :: local_molecules 62 TYPE(molecule_kind_type), POINTER :: molecule_kind 63 TYPE(particle_type), POINTER :: particle_set(:) 64 65 logger => cp_get_default_logger() 66 67 particle_set => qs_loc_env%particle_set 68 para_env => qs_loc_env%para_env 69 local_molecules => qs_loc_env%local_molecules 70 nstate = SIZE(center, 2) 71 ALLOCATE (wfc_to_atom_map(nstate)) 72 !--------------------------------------------------------------------------- 73 !--------------------------------------------------------------------------- 74 nkind = SIZE(local_molecules%n_el) 75 natom = 0 76 natom_max = 0 77 DO ikind = 1, nkind 78 nmol = SIZE(local_molecules%list(ikind)%array) 79 DO imol = 1, nmol 80 i = local_molecules%list(ikind)%array(imol) 81 molecule_kind => molecule_set(i)%molecule_kind 82 CALL get_molecule_kind(molecule_kind=molecule_kind, natom=natom) 83 natom_max = natom_max + natom 84 IF (.NOT. ASSOCIATED(molecule_set(i)%lmi)) THEN 85 ALLOCATE (molecule_set(i)%lmi(nspins)) 86 DO k = 1, nspins 87 NULLIFY (molecule_set(i)%lmi(k)%states) 88 END DO 89 ENDIF 90 molecule_set(i)%lmi(ispin)%nstates = 0 91 IF (ASSOCIATED(molecule_set(i)%lmi(ispin)%states)) THEN 92 DEALLOCATE (molecule_set(i)%lmi(ispin)%states) 93 END IF 94 END DO 95 END DO 96 natom_loc = natom_max 97 natom = natom_max 98 99 CALL mp_max(natom_max, para_env%group) 100 101 ALLOCATE (r(3, natom_max)) 102 103 ALLOCATE (distance(natom_max)) 104 105 !Zero all the stuff 106 r(:, :) = 0.0_dp 107 distance(:) = 1.E10_dp 108 109 !--------------------------------------------------------------------------- 110 !--------------------------------------------------------------------------- 111 counter = 0 112 nkind = SIZE(local_molecules%n_el) 113 DO ikind = 1, nkind 114 nmol = SIZE(local_molecules%list(ikind)%array) 115 DO imol = 1, nmol 116 i = local_molecules%list(ikind)%array(imol) 117 molecule_kind => molecule_set(i)%molecule_kind 118 first_atom = molecule_set(i)%first_atom 119 CALL get_molecule_kind(molecule_kind=molecule_kind, natom=natom) 120 121 DO iatom = 1, natom 122 counter = counter + 1 123 r(:, counter) = particle_set(first_atom + iatom - 1)%r(:) 124 END DO 125 END DO 126 END DO 127 128 !--------------------------------------------------------------------------- 129 !--------------------------------------------------------------------------- 130 DO istate = 1, nstate 131 distance(:) = 1.E10_dp 132 DO iatom = 1, natom_loc 133 dr(1) = r(1, iatom) - center(1, istate) 134 dr(2) = r(2, iatom) - center(2, istate) 135 dr(3) = r(3, iatom) - center(3, istate) 136 ria = pbc(dr, qs_loc_env%cell) 137 distance(iatom) = SQRT(DOT_PRODUCT(ria, ria)) 138 END DO 139 140 !combine distance() from all procs 141 local_location = MAX(1, MINLOC(distance, DIM=1)) 142 143 mydist(1) = distance(local_location) 144 mydist(2) = para_env%mepos 145 146 CALL mp_minloc(mydist, para_env%group) 147 148 IF (mydist(2) == para_env%mepos) THEN 149 wfc_to_atom_map(istate) = local_location 150 ELSE 151 wfc_to_atom_map(istate) = 0 152 END IF 153 END DO 154 !--------------------------------------------------------------------------- 155 !--------------------------------------------------------------------------- 156 IF (natom_loc /= 0) THEN 157 DO istate = 1, nstate 158 iatom = wfc_to_atom_map(istate) 159 IF (iatom /= 0) THEN 160 counter = 0 161 nkind = SIZE(local_molecules%n_el) 162 DO ikind = 1, nkind 163 nmol = SIZE(local_molecules%list(ikind)%array) 164 DO imol = 1, nmol 165 imol_now = local_molecules%list(ikind)%array(imol) 166 molecule_kind => molecule_set(imol_now)%molecule_kind 167 CALL get_molecule_kind(molecule_kind=molecule_kind, natom=natom) 168 counter = counter + natom 169 IF (counter >= iatom) EXIT 170 END DO 171 IF (counter >= iatom) EXIT 172 END DO 173 i = molecule_set(imol_now)%lmi(ispin)%nstates 174 i = i + 1 175 molecule_set(imol_now)%lmi(ispin)%nstates = i 176 CALL reallocate(molecule_set(imol_now)%lmi(ispin)%states, 1, i) 177 molecule_set(imol_now)%lmi(ispin)%states(i) = istate 178 END IF 179 END DO 180 END IF 181 182 DEALLOCATE (distance) 183 DEALLOCATE (r) 184 DEALLOCATE (wfc_to_atom_map) 185 186 END SUBROUTINE wfc_to_molecule 187 !------------------------------------------------------------------------------ 188 189END MODULE qs_loc_molecules 190 191