1!--------------------------------------------------------------------------------------------------! 2! CP2K: A general program to perform molecular dynamics simulations ! 3! Copyright (C) 2000 - 2019 CP2K developers group ! 4!--------------------------------------------------------------------------------------------------! 5 6! ************************************************************************************************** 7!> \par History 8!> Harald Forbert (Dec-2000): Changes for multiple linked lists 9!> linklist_internal_data_type 10!> 07.02.2005: using real coordinates for r_last_update; cleaned (MK) 11!> \author CJM,MK 12! ************************************************************************************************** 13MODULE fist_neighbor_list_control 14 15 USE atomic_kind_types, ONLY: atomic_kind_type,& 16 get_atomic_kind_set 17 USE cell_types, ONLY: cell_clone,& 18 cell_create,& 19 cell_release,& 20 cell_type 21 USE cp_log_handling, ONLY: cp_get_default_logger,& 22 cp_logger_type 23 USE cp_output_handling, ONLY: cp_print_key_finished_output,& 24 cp_print_key_unit_nr 25 USE cp_para_types, ONLY: cp_para_env_type 26 USE distribution_1d_types, ONLY: distribution_1d_type 27 USE exclusion_types, ONLY: exclusion_type 28 USE fist_neighbor_list_types, ONLY: fist_neighbor_type 29 USE fist_neighbor_lists, ONLY: build_fist_neighbor_lists 30 USE fist_nonbond_env_types, ONLY: fist_nonbond_env_get,& 31 fist_nonbond_env_set,& 32 fist_nonbond_env_type,& 33 pos_type 34 USE input_section_types, ONLY: section_vals_type,& 35 section_vals_val_get 36 USE kinds, ONLY: dp 37 USE message_passing, ONLY: mp_max 38 USE pair_potential_types, ONLY: pair_potential_pp_type,& 39 siepmann_type,& 40 tersoff_type 41 USE particle_types, ONLY: particle_type 42#include "./base/base_uses.f90" 43 44 IMPLICIT NONE 45 46 PRIVATE 47 48 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'fist_neighbor_list_control' 49 50 PUBLIC :: list_control 51 52!*** 53 54CONTAINS 55 56! to decide whether the neighbor list is to be updated or not 57! based on a displacement criterion; 58! if any particle has moved by 0.5*verlet_skin from the previous 59! list update, then the list routine is called. 60 61! ************************************************************************************************** 62!> \brief ... 63!> \param atomic_kind_set ... 64!> \param particle_set ... 65!> \param local_particles ... 66!> \param cell ... 67!> \param fist_nonbond_env ... 68!> \param para_env ... 69!> \param mm_section ... 70!> \param shell_particle_set ... 71!> \param core_particle_set ... 72!> \param force_update ... 73!> \param exclusions ... 74! ************************************************************************************************** 75 SUBROUTINE list_control(atomic_kind_set, particle_set, local_particles, & 76 cell, fist_nonbond_env, para_env, mm_section, shell_particle_set, & 77 core_particle_set, force_update, exclusions) 78 79 TYPE(atomic_kind_type), POINTER :: atomic_kind_set(:) 80 TYPE(particle_type), POINTER :: particle_set(:) 81 TYPE(distribution_1d_type), POINTER :: local_particles 82 TYPE(cell_type), POINTER :: cell 83 TYPE(fist_nonbond_env_type), POINTER :: fist_nonbond_env 84 TYPE(cp_para_env_type), POINTER :: para_env 85 TYPE(section_vals_type), POINTER :: mm_section 86 TYPE(particle_type), OPTIONAL, POINTER :: shell_particle_set(:), & 87 core_particle_set(:) 88 LOGICAL, INTENT(IN), OPTIONAL :: force_update 89 TYPE(exclusion_type), DIMENSION(:), OPTIONAL :: exclusions 90 91 CHARACTER(LEN=*), PARAMETER :: routineN = 'list_control', routineP = moduleN//':'//routineN 92 93 INTEGER :: counter, handle, ikind, iparticle, iparticle_kind, iparticle_local, ishell, & 94 jkind, last_update, nparticle, nparticle_kind, nparticle_local, nshell, num_update, & 95 output_unit 96 LOGICAL :: build_from_scratch, geo_check, & 97 shell_adiabatic, shell_present, & 98 update_neighbor_lists 99 LOGICAL, DIMENSION(:, :), POINTER :: full_nl 100 REAL(KIND=dp) :: aup, dr2, dr2_max, ei_scale14, lup, & 101 vdw_scale14, verlet_skin 102 REAL(KIND=dp), DIMENSION(3) :: dr, rab, s, s2r 103 REAL(KIND=dp), DIMENSION(:, :), POINTER :: rlist_cut, rlist_lowsq 104 TYPE(cell_type), POINTER :: cell_last_update 105 TYPE(cp_logger_type), POINTER :: logger 106 TYPE(fist_neighbor_type), POINTER :: nonbonded 107 TYPE(pair_potential_pp_type), POINTER :: potparm 108 TYPE(pos_type), DIMENSION(:), POINTER :: r_last_update, r_last_update_pbc, & 109 rcore_last_update_pbc, & 110 rshell_last_update_pbc 111 112 CALL timeset(routineN, handle) 113 NULLIFY (logger) 114 logger => cp_get_default_logger() 115 116 ! *** Assigning local pointers *** 117 CALL fist_nonbond_env_get(fist_nonbond_env, & 118 nonbonded=nonbonded, & 119 rlist_cut=rlist_cut, & 120 rlist_lowsq=rlist_lowsq, & 121 aup=aup, & 122 lup=lup, & 123 ei_scale14=ei_scale14, & 124 vdw_scale14=vdw_scale14, & 125 counter=counter, & 126 r_last_update=r_last_update, & 127 r_last_update_pbc=r_last_update_pbc, & 128 rshell_last_update_pbc=rshell_last_update_pbc, & 129 rcore_last_update_pbc=rcore_last_update_pbc, & 130 cell_last_update=cell_last_update, & 131 num_update=num_update, & 132 potparm=potparm, & 133 last_update=last_update) 134 135 nparticle = SIZE(particle_set) 136 nparticle_kind = SIZE(atomic_kind_set) 137 nshell = 0 138 CALL get_atomic_kind_set(atomic_kind_set=atomic_kind_set, & 139 shell_present=shell_present, shell_adiabatic=shell_adiabatic) 140 IF (shell_present) THEN 141 nshell = SIZE(shell_particle_set) 142 END IF 143 144 ! *** Check, if the neighbor lists have to be built or updated *** 145 update_neighbor_lists = .FALSE. 146 CALL section_vals_val_get(mm_section, "NEIGHBOR_LISTS%NEIGHBOR_LISTS_FROM_SCRATCH", & 147 l_val=build_from_scratch) 148 CALL section_vals_val_get(mm_section, "NEIGHBOR_LISTS%GEO_CHECK", & 149 l_val=geo_check) 150 IF (ASSOCIATED(r_last_update)) THEN 151 ! Determine the maximum of the squared displacement, compared to 152 ! r_last_update. 153 CALL section_vals_val_get(mm_section, "NEIGHBOR_LISTS%VERLET_SKIN", & 154 r_val=verlet_skin) 155 dr2_max = 0.0_dp 156 DO iparticle_kind = 1, nparticle_kind 157 nparticle_local = local_particles%n_el(iparticle_kind) 158 DO iparticle_local = 1, nparticle_local 159 iparticle = local_particles%list(iparticle_kind)%array(iparticle_local) 160 s2r = r_last_update(iparticle)%r 161 s = particle_set(iparticle)%r(:) 162 dr(:) = s2r - s 163 dr2 = dr(1)*dr(1) + dr(2)*dr(2) + dr(3)*dr(3) 164 dr2_max = MAX(dr2_max, dr2) 165 END DO 166 END DO 167 168 CALL mp_max(dr2_max, para_env%group) 169 170 ! If the maximum distplacement is too large, ... 171 IF (dr2_max > 0.25_dp*verlet_skin**2 .OR. build_from_scratch) THEN 172 DO iparticle = 1, nparticle 173 r_last_update(iparticle)%r = particle_set(iparticle)%r(:) 174 END DO 175 update_neighbor_lists = .TRUE. 176 END IF 177 ELSE 178 ! There is no r_last_update to compare with. Neighbor lists from scratch. 179 ALLOCATE (r_last_update(nparticle)) 180 DO iparticle = 1, nparticle 181 r_last_update(iparticle)%r = particle_set(iparticle)%r(:) 182 END DO 183 184 update_neighbor_lists = .TRUE. 185 build_from_scratch = .TRUE. 186 END IF 187 ! Force Update 188 IF (PRESENT(force_update)) THEN 189 IF (force_update) update_neighbor_lists = .TRUE. 190 END IF 191 192 ! Allocate the r_last_update_pbc, rshell_last_update_pbc, rcore_last_update_pbc 193 IF (.NOT. ASSOCIATED(r_last_update_pbc)) THEN 194 ALLOCATE (r_last_update_pbc(nparticle)) 195 END IF 196 IF (shell_present .AND. .NOT. ASSOCIATED(rshell_last_update_pbc)) THEN 197 ALLOCATE (rshell_last_update_pbc(nshell)) 198 END IF 199 IF (shell_present .AND. .NOT. ASSOCIATED(rcore_last_update_pbc)) THEN 200 ALLOCATE (rcore_last_update_pbc(nshell)) 201 END IF 202 203 ! update the neighbor lists 204 IF (update_neighbor_lists) THEN 205 ! determine which pairs of atom kinds need full neighbor lists. Full 206 ! means that atom a is in the neighbor list of atom b and vice versa. 207 ALLOCATE (full_nl(nparticle_kind, nparticle_kind)) 208 IF (ASSOCIATED(potparm)) THEN 209 DO ikind = 1, nparticle_kind 210 DO jkind = ikind, nparticle_kind 211 full_nl(ikind, jkind) = .FALSE. 212 IF (ANY(potparm%pot(ikind, jkind)%pot%type == tersoff_type)) THEN 213 full_nl(ikind, jkind) = .TRUE. 214 END IF 215 IF (ANY(potparm%pot(ikind, jkind)%pot%type == siepmann_type)) THEN 216 full_nl(ikind, jkind) = .TRUE. 217 END IF 218 full_nl(jkind, ikind) = full_nl(ikind, jkind) 219 END DO 220 END DO 221 ELSE 222 full_nl = .FALSE. 223 END IF 224 CALL build_fist_neighbor_lists(atomic_kind_set, particle_set, & 225 local_particles, cell, rlist_cut, rlist_lowsq, ei_scale14, & 226 vdw_scale14, nonbonded, para_env, & 227 build_from_scratch=build_from_scratch, geo_check=geo_check, & 228 mm_section=mm_section, full_nl=full_nl, & 229 exclusions=exclusions) 230 231 CALL cell_release(cell_last_update) 232 CALL cell_create(cell_last_update) 233 CALL cell_clone(cell, cell_last_update) 234 235 IF (counter > 0) THEN 236 num_update = num_update + 1 237 lup = counter + 1 - last_update 238 last_update = counter + 1 239 aup = aup + (lup - aup)/REAL(num_update, KIND=dp) 240 ELSE 241 num_update = 0 242 lup = 0 243 last_update = 1 244 aup = 0.0_dp 245 END IF 246 247 CALL fist_nonbond_env_set(fist_nonbond_env, & 248 lup=lup, & 249 aup=aup, & 250 r_last_update=r_last_update, & 251 r_last_update_pbc=r_last_update_pbc, & 252 rshell_last_update_pbc=rshell_last_update_pbc, & 253 rcore_last_update_pbc=rcore_last_update_pbc, & 254 nonbonded=nonbonded, & 255 num_update=num_update, & 256 last_update=last_update, & 257 cell_last_update=cell_last_update) 258 259 output_unit = cp_print_key_unit_nr(logger, mm_section, "PRINT%NEIGHBOR_LISTS", & 260 extension=".mmLog") 261 IF (output_unit > 0) THEN 262 WRITE (UNIT=output_unit, & 263 FMT="(/,T2,A,/,T52,A,/,A,T31,A,T49,2(1X,F15.2),/,T2,A,/)") & 264 REPEAT("*", 79), "INSTANTANEOUS AVERAGES", & 265 " LIST UPDATES[steps]", "= ", lup, aup, REPEAT("*", 79) 266 END IF 267 CALL cp_print_key_finished_output(output_unit, logger, mm_section, & 268 "PRINT%NEIGHBOR_LISTS") 269 DEALLOCATE (full_nl) 270 END IF 271 272 ! Store particle positions after the last update, translated to the 273 ! primitive cell, in r_last_update_pbc. 274 DO iparticle = 1, nparticle 275 rab = r_last_update(iparticle)%r 276 IF (cell%orthorhombic) THEN 277 rab(1) = -cell%hmat(1, 1)*cell%perd(1)*ANINT(cell_last_update%h_inv(1, 1)*rab(1)) 278 rab(2) = -cell%hmat(2, 2)*cell%perd(2)*ANINT(cell_last_update%h_inv(2, 2)*rab(2)) 279 rab(3) = -cell%hmat(3, 3)*cell%perd(3)*ANINT(cell_last_update%h_inv(3, 3)*rab(3)) 280 ELSE 281 s(1) = cell_last_update%h_inv(1, 1)*rab(1) + cell_last_update%h_inv(1, 2)*rab(2) + & 282 cell_last_update%h_inv(1, 3)*rab(3) 283 s(2) = cell_last_update%h_inv(2, 1)*rab(1) + cell_last_update%h_inv(2, 2)*rab(2) + & 284 cell_last_update%h_inv(2, 3)*rab(3) 285 s(3) = cell_last_update%h_inv(3, 1)*rab(1) + cell_last_update%h_inv(3, 2)*rab(2) + & 286 cell_last_update%h_inv(3, 3)*rab(3) 287 s(1) = -cell%perd(1)*ANINT(s(1)) 288 s(2) = -cell%perd(2)*ANINT(s(2)) 289 s(3) = -cell%perd(3)*ANINT(s(3)) 290 rab(1) = +cell%hmat(1, 1)*s(1) + cell%hmat(1, 2)*s(2) + cell%hmat(1, 3)*s(3) 291 rab(2) = +cell%hmat(2, 1)*s(1) + cell%hmat(2, 2)*s(2) + cell%hmat(2, 3)*s(3) 292 rab(3) = +cell%hmat(3, 1)*s(1) + cell%hmat(3, 2)*s(2) + cell%hmat(3, 3)*s(3) 293 END IF 294 r_last_update_pbc(iparticle)%r = particle_set(iparticle)%r + rab 295 ! Use the same translation for core and shell. 296 ishell = particle_set(iparticle)%shell_index 297 IF (ishell /= 0) THEN 298 rshell_last_update_pbc(ishell)%r = rab + shell_particle_set(ishell)%r(:) 299 IF (shell_adiabatic) THEN 300 rcore_last_update_pbc(ishell)%r = rab + core_particle_set(ishell)%r(:) 301 ELSE 302 rcore_last_update_pbc(ishell)%r = r_last_update_pbc(iparticle)%r(:) 303 END IF 304 END IF 305 END DO 306 307 counter = counter + 1 308 CALL fist_nonbond_env_set(fist_nonbond_env, counter=counter) 309 CALL timestop(handle) 310 311 END SUBROUTINE list_control 312 313END MODULE fist_neighbor_list_control 314