1!--------------------------------------------------------------------------------------------------! 2! CP2K: A general program to perform molecular dynamics simulations ! 3! Copyright (C) 2000 - 2020 CP2K developers group ! 4!--------------------------------------------------------------------------------------------------! 5 6! ************************************************************************************************** 7!> \par History 8!> 09.2004 created [tlaino] 9!> \author Teodoro Laino 10! ************************************************************************************************** 11MODULE qmmm_util 12 USE cell_types, ONLY: cell_type 13 USE cp_log_handling, ONLY: cp_logger_get_default_io_unit 14 USE cp_subsys_types, ONLY: cp_subsys_type 15 USE fist_environment_types, ONLY: fist_env_get 16 USE force_env_types, ONLY: force_env_type,& 17 use_qmmm,& 18 use_qmmmx 19 USE input_constants, ONLY: do_qmmm_wall_none,& 20 do_qmmm_wall_quadratic,& 21 do_qmmm_wall_reflective 22 USE input_section_types, ONLY: section_vals_get,& 23 section_vals_get_subs_vals,& 24 section_vals_type,& 25 section_vals_val_get 26 USE kinds, ONLY: dp 27 USE mathconstants, ONLY: pi 28 USE particle_methods, ONLY: write_fist_particle_coordinates,& 29 write_qs_particle_coordinates 30 USE particle_types, ONLY: particle_type 31 USE qmmm_types, ONLY: qmmm_env_type 32 USE qs_energy_types, ONLY: qs_energy_type 33 USE qs_environment_types, ONLY: get_qs_env 34 USE qs_kind_types, ONLY: qs_kind_type 35#include "./base/base_uses.f90" 36 37 IMPLICIT NONE 38 PRIVATE 39 40 LOGICAL, PRIVATE, PARAMETER :: debug_this_module = .FALSE. 41 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qmmm_util' 42 PUBLIC :: apply_qmmm_walls_reflective, & 43 apply_qmmm_walls, & 44 apply_qmmm_translate, & 45 apply_qmmm_wrap, & 46 apply_qmmm_unwrap, & 47 spherical_cutoff_factor 48 49CONTAINS 50 51! ************************************************************************************************** 52!> \brief Apply QM quadratic walls in order to avoid QM atoms escaping from 53!> the QM Box 54!> \param qmmm_env ... 55!> \par History 56!> 02.2008 created 57!> \author Benjamin G Levine 58! ************************************************************************************************** 59 SUBROUTINE apply_qmmm_walls(qmmm_env) 60 TYPE(qmmm_env_type), POINTER :: qmmm_env 61 62 CHARACTER(len=*), PARAMETER :: routineN = 'apply_qmmm_walls', & 63 routineP = moduleN//':'//routineN 64 65 INTEGER :: iwall_type 66 LOGICAL :: do_qmmm_force_mixing, explicit 67 TYPE(section_vals_type), POINTER :: qmmmx_section, walls_section 68 69 walls_section => section_vals_get_subs_vals(qmmm_env%qs_env%input, "QMMM%WALLS") 70 qmmmx_section => section_vals_get_subs_vals(qmmm_env%qs_env%input, "QMMM%FORCE_MIXING") 71 CALL section_vals_get(qmmmx_section, explicit=do_qmmm_force_mixing) 72 CALL section_vals_get(walls_section, explicit=explicit) 73 IF (explicit) THEN 74 CALL section_vals_val_get(walls_section, "TYPE", i_val=iwall_type) 75 SELECT CASE (iwall_type) 76 CASE (do_qmmm_wall_quadratic) 77 IF (do_qmmm_force_mixing) THEN 78 CALL cp_warn(__LOCATION__, & 79 "Quadratic walls for QM/MM are not implemented (or useful), when "// & 80 "force mixing is active. Skipping!") 81 ELSE 82 CALL apply_qmmm_walls_quadratic(qmmm_env, walls_section) 83 ENDIF 84 CASE (do_qmmm_wall_reflective) 85 ! Do nothing.. reflective walls are applied directly in the integrator 86 END SELECT 87 ENDIF 88 89 END SUBROUTINE apply_qmmm_walls 90 91! ************************************************************************************************** 92!> \brief Apply reflective QM walls in order to avoid QM atoms escaping from 93!> the QM Box 94!> \param force_env ... 95!> \par History 96!> 08.2007 created [tlaino] - Zurich University 97!> \author Teodoro Laino 98! ************************************************************************************************** 99 SUBROUTINE apply_qmmm_walls_reflective(force_env) 100 TYPE(force_env_type), POINTER :: force_env 101 102 CHARACTER(len=*), PARAMETER :: routineN = 'apply_qmmm_walls_reflective', & 103 routineP = moduleN//':'//routineN 104 105 INTEGER :: ip, iwall_type, qm_index 106 INTEGER, DIMENSION(:), POINTER :: qm_atom_index 107 LOGICAL :: explicit, is_x(2), is_y(2), is_z(2) 108 REAL(KIND=dp), DIMENSION(3) :: coord, qm_cell_diag, skin 109 REAL(KIND=dp), DIMENSION(:), POINTER :: list 110 TYPE(cell_type), POINTER :: mm_cell, qm_cell 111 TYPE(cp_subsys_type), POINTER :: subsys_mm, subsys_qm 112 TYPE(particle_type), DIMENSION(:), POINTER :: particles_mm 113 TYPE(section_vals_type), POINTER :: walls_section 114 115 NULLIFY (subsys_mm, subsys_qm, qm_atom_index, particles_mm, qm_cell, mm_cell, & 116 walls_section) 117 118 IF (force_env%in_use /= use_qmmm .AND. force_env%in_use /= use_qmmmx) RETURN 119 120 walls_section => section_vals_get_subs_vals(force_env%root_section, "FORCE_EVAL%QMMM%WALLS") 121 CALL section_vals_get(walls_section, explicit=explicit) 122 IF (explicit) THEN 123 NULLIFY (list) 124 CALL section_vals_val_get(walls_section, "WALL_SKIN", r_vals=list) 125 CALL section_vals_val_get(walls_section, "TYPE", i_val=iwall_type) 126 skin(:) = list(:) 127 ELSE 128 ![NB] 129 iwall_type = do_qmmm_wall_reflective 130 skin(:) = 0.0_dp 131 END IF 132 133 IF (force_env%in_use == use_qmmmx) THEN 134 IF (iwall_type /= do_qmmm_wall_none) & 135 CALL cp_warn(__LOCATION__, & 136 "Reflective walls for QM/MM are not implemented (or useful) when "// & 137 "force mixing is active. Skipping!") 138 RETURN 139 END IF 140 141 ! from here on we can be sure that it's conventional QM/MM 142 CPASSERT(ASSOCIATED(force_env%qmmm_env)) 143 CPASSERT(force_env%qmmm_env%ref_count > 0) 144 145 CALL fist_env_get(force_env%qmmm_env%fist_env, cell=mm_cell, subsys=subsys_mm) 146 CALL get_qs_env(force_env%qmmm_env%qs_env, cell=qm_cell, cp_subsys=subsys_qm) 147 qm_atom_index => force_env%qmmm_env%qm%qm_atom_index 148 CPASSERT(ASSOCIATED(qm_atom_index)) 149 150 qm_cell_diag = (/qm_cell%hmat(1, 1), & 151 qm_cell%hmat(2, 2), & 152 qm_cell%hmat(3, 3)/) 153 particles_mm => subsys_mm%particles%els 154 DO ip = 1, SIZE(qm_atom_index) 155 qm_index = qm_atom_index(ip) 156 coord = particles_mm(qm_index)%r 157 IF (ANY(coord < skin) .OR. ANY(coord > (qm_cell_diag - skin))) THEN 158 IF (explicit) THEN 159 IF (iwall_type == do_qmmm_wall_reflective) THEN 160 ! Apply Walls 161 is_x(1) = (coord(1) < skin(1)) 162 is_x(2) = (coord(1) > (qm_cell_diag(1) - skin(1))) 163 is_y(1) = (coord(2) < skin(2)) 164 is_y(2) = (coord(2) > (qm_cell_diag(2) - skin(2))) 165 is_z(1) = (coord(3) < skin(3)) 166 is_z(2) = (coord(3) > (qm_cell_diag(3) - skin(3))) 167 IF (ANY(is_x)) THEN 168 ! X coordinate 169 IF (is_x(1)) THEN 170 particles_mm(qm_index)%v(1) = ABS(particles_mm(qm_index)%v(1)) 171 ELSE IF (is_x(2)) THEN 172 particles_mm(qm_index)%v(1) = -ABS(particles_mm(qm_index)%v(1)) 173 END IF 174 END IF 175 IF (ANY(is_y)) THEN 176 ! Y coordinate 177 IF (is_y(1)) THEN 178 particles_mm(qm_index)%v(2) = ABS(particles_mm(qm_index)%v(2)) 179 ELSE IF (is_y(2)) THEN 180 particles_mm(qm_index)%v(2) = -ABS(particles_mm(qm_index)%v(2)) 181 END IF 182 END IF 183 IF (ANY(is_z)) THEN 184 ! Z coordinate 185 IF (is_z(1)) THEN 186 particles_mm(qm_index)%v(3) = ABS(particles_mm(qm_index)%v(3)) 187 ELSE IF (is_z(2)) THEN 188 particles_mm(qm_index)%v(3) = -ABS(particles_mm(qm_index)%v(3)) 189 END IF 190 END IF 191 ENDIF 192 ELSE 193 ! Otherwise print a warning and continue crossing cp2k's finger.. 194 CALL cp_warn(__LOCATION__, & 195 "One or few QM atoms are within the SKIN of the quantum box. Check your run "// & 196 "and you may possibly consider: the activation of the QMMM WALLS "// & 197 "around the QM box, switching ON the centering of the QM box or increase "// & 198 "the size of the QM cell. CP2K CONTINUE but results could be meaningless. ") 199 END IF 200 END IF 201 END DO 202 203 END SUBROUTINE apply_qmmm_walls_reflective 204 205! ************************************************************************************************** 206!> \brief Apply QM quadratic walls in order to avoid QM atoms escaping from 207!> the QM Box 208!> \param qmmm_env ... 209!> \param walls_section ... 210!> \par History 211!> 02.2008 created 212!> \author Benjamin G Levine 213! ************************************************************************************************** 214 SUBROUTINE apply_qmmm_walls_quadratic(qmmm_env, walls_section) 215 TYPE(qmmm_env_type), POINTER :: qmmm_env 216 TYPE(section_vals_type), POINTER :: walls_section 217 218 CHARACTER(len=*), PARAMETER :: routineN = 'apply_qmmm_walls_quadratic', & 219 routineP = moduleN//':'//routineN 220 221 INTEGER :: ip, qm_index 222 INTEGER, DIMENSION(:), POINTER :: qm_atom_index 223 LOGICAL :: is_x(2), is_y(2), is_z(2) 224 REAL(KIND=dp) :: k, wallenergy, wallforce 225 REAL(KIND=dp), DIMENSION(3) :: coord, qm_cell_diag, skin 226 REAL(KIND=dp), DIMENSION(:), POINTER :: list 227 TYPE(cell_type), POINTER :: mm_cell, qm_cell 228 TYPE(cp_subsys_type), POINTER :: subsys_mm, subsys_qm 229 TYPE(particle_type), DIMENSION(:), POINTER :: particles_mm 230 TYPE(qs_energy_type), POINTER :: energy 231 232 NULLIFY (list) 233 CALL section_vals_val_get(walls_section, "WALL_SKIN", r_vals=list) 234 CALL section_vals_val_get(walls_section, "K", r_val=k) 235 CPASSERT(ASSOCIATED(qmmm_env)) 236 CPASSERT(qmmm_env%ref_count > 0) 237 238 CALL fist_env_get(qmmm_env%fist_env, cell=mm_cell, subsys=subsys_mm) 239 CALL get_qs_env(qmmm_env%qs_env, cell=qm_cell, cp_subsys=subsys_qm) 240 241 qm_atom_index => qmmm_env%qm%qm_atom_index 242 CPASSERT(ASSOCIATED(qm_atom_index)) 243 244 skin(:) = list(:) 245 246 qm_cell_diag = (/qm_cell%hmat(1, 1), & 247 qm_cell%hmat(2, 2), & 248 qm_cell%hmat(3, 3)/) 249 particles_mm => subsys_mm%particles%els 250 wallenergy = 0.0_dp 251 DO ip = 1, SIZE(qm_atom_index) 252 qm_index = qm_atom_index(ip) 253 coord = particles_mm(qm_index)%r 254 IF (ANY(coord < skin) .OR. ANY(coord > (qm_cell_diag - skin))) THEN 255 is_x(1) = (coord(1) < skin(1)) 256 is_x(2) = (coord(1) > (qm_cell_diag(1) - skin(1))) 257 is_y(1) = (coord(2) < skin(2)) 258 is_y(2) = (coord(2) > (qm_cell_diag(2) - skin(2))) 259 is_z(1) = (coord(3) < skin(3)) 260 is_z(2) = (coord(3) > (qm_cell_diag(3) - skin(3))) 261 IF (is_x(1)) THEN 262 wallforce = 2.0_dp*k*(skin(1) - coord(1)) 263 particles_mm(qm_index)%f(1) = particles_mm(qm_index)%f(1) + & 264 wallforce 265 wallenergy = wallenergy + wallforce*(skin(1) - coord(1))*0.5_dp 266 ENDIF 267 IF (is_x(2)) THEN 268 wallforce = 2.0_dp*k*((qm_cell_diag(1) - skin(1)) - coord(1)) 269 particles_mm(qm_index)%f(1) = particles_mm(qm_index)%f(1) + & 270 wallforce 271 wallenergy = wallenergy + wallforce*((qm_cell_diag(1) - skin(1)) - & 272 coord(1))*0.5_dp 273 ENDIF 274 IF (is_y(1)) THEN 275 wallforce = 2.0_dp*k*(skin(2) - coord(2)) 276 particles_mm(qm_index)%f(2) = particles_mm(qm_index)%f(2) + & 277 wallforce 278 wallenergy = wallenergy + wallforce*(skin(2) - coord(2))*0.5_dp 279 ENDIF 280 IF (is_y(2)) THEN 281 wallforce = 2.0_dp*k*((qm_cell_diag(2) - skin(2)) - coord(2)) 282 particles_mm(qm_index)%f(2) = particles_mm(qm_index)%f(2) + & 283 wallforce 284 wallenergy = wallenergy + wallforce*((qm_cell_diag(2) - skin(2)) - & 285 coord(2))*0.5_dp 286 ENDIF 287 IF (is_z(1)) THEN 288 wallforce = 2.0_dp*k*(skin(3) - coord(3)) 289 particles_mm(qm_index)%f(3) = particles_mm(qm_index)%f(3) + & 290 wallforce 291 wallenergy = wallenergy + wallforce*(skin(3) - coord(3))*0.5_dp 292 ENDIF 293 IF (is_z(2)) THEN 294 wallforce = 2.0_dp*k*((qm_cell_diag(3) - skin(3)) - coord(3)) 295 particles_mm(qm_index)%f(3) = particles_mm(qm_index)%f(3) + & 296 wallforce 297 wallenergy = wallenergy + wallforce*((qm_cell_diag(3) - skin(3)) - & 298 coord(3))*0.5_dp 299 ENDIF 300 ENDIF 301 ENDDO 302 303 CALL get_qs_env(qs_env=qmmm_env%qs_env, energy=energy) 304 energy%total = energy%total + wallenergy 305 306 END SUBROUTINE apply_qmmm_walls_quadratic 307 308! ************************************************************************************************** 309!> \brief wrap positions (with mm periodicity) 310!> \param subsys_mm ... 311!> \param mm_cell ... 312!> \param subsys_qm ... 313!> \param qm_atom_index ... 314!> \param saved_pos ... 315! ************************************************************************************************** 316 SUBROUTINE apply_qmmm_wrap(subsys_mm, mm_cell, subsys_qm, qm_atom_index, saved_pos) 317 TYPE(cp_subsys_type), POINTER :: subsys_mm 318 TYPE(cell_type), POINTER :: mm_cell 319 TYPE(cp_subsys_type), OPTIONAL, POINTER :: subsys_qm 320 INTEGER, DIMENSION(:), OPTIONAL, POINTER :: qm_atom_index 321 REAL(dp), ALLOCATABLE :: saved_pos(:, :) 322 323 CHARACTER(len=*), PARAMETER :: routineN = 'apply_qmmm_wrap', & 324 routineP = moduleN//':'//routineN 325 326 INTEGER :: i_dim, ip 327 REAL(dp) :: r_lat(3) 328 329 ALLOCATE (saved_pos(3, subsys_mm%particles%n_els)) 330 DO ip = 1, subsys_mm%particles%n_els 331 saved_pos(1:3, ip) = subsys_mm%particles%els(ip)%r(1:3) 332 r_lat = MATMUL(mm_cell%h_inv, subsys_mm%particles%els(ip)%r) 333 DO i_dim = 1, 3 334 IF (mm_cell%perd(i_dim) /= 1) THEN 335 r_lat(i_dim) = 0.0_dp 336 END IF 337 END DO 338 subsys_mm%particles%els(ip)%r = subsys_mm%particles%els(ip)%r - MATMUL(mm_cell%hmat, FLOOR(r_lat)) 339 END DO 340 341 IF (PRESENT(subsys_qm) .AND. PRESENT(qm_atom_index)) THEN 342 DO ip = 1, SIZE(qm_atom_index) 343 subsys_qm%particles%els(ip)%r = subsys_mm%particles%els(qm_atom_index(ip))%r 344 END DO 345 END IF 346 END SUBROUTINE apply_qmmm_wrap 347 348! ************************************************************************************************** 349!> \brief ... 350!> \param subsys_mm ... 351!> \param subsys_qm ... 352!> \param qm_atom_index ... 353!> \param saved_pos ... 354! ************************************************************************************************** 355 SUBROUTINE apply_qmmm_unwrap(subsys_mm, subsys_qm, qm_atom_index, saved_pos) 356 TYPE(cp_subsys_type), POINTER :: subsys_mm 357 TYPE(cp_subsys_type), OPTIONAL, POINTER :: subsys_qm 358 INTEGER, DIMENSION(:), OPTIONAL, POINTER :: qm_atom_index 359 REAL(dp), ALLOCATABLE :: saved_pos(:, :) 360 361 CHARACTER(len=*), PARAMETER :: routineN = 'apply_qmmm_unwrap', & 362 routineP = moduleN//':'//routineN 363 364 INTEGER :: ip 365 366 DO ip = 1, subsys_mm%particles%n_els 367 subsys_mm%particles%els(ip)%r(1:3) = saved_pos(1:3, ip) 368 END DO 369 370 IF (PRESENT(subsys_qm) .AND. PRESENT(qm_atom_index)) THEN 371 DO ip = 1, SIZE(qm_atom_index) 372 subsys_qm%particles%els(ip)%r = subsys_mm%particles%els(qm_atom_index(ip))%r 373 END DO 374 END IF 375 376 DEALLOCATE (saved_pos) 377 END SUBROUTINE apply_qmmm_unwrap 378 379! ************************************************************************************************** 380!> \brief Apply translation to the full system in order to center the QM 381!> system into the QM box 382!> \param qmmm_env ... 383!> \par History 384!> 08.2007 created [tlaino] - Zurich University 385!> \author Teodoro Laino 386! ************************************************************************************************** 387 SUBROUTINE apply_qmmm_translate(qmmm_env) 388 TYPE(qmmm_env_type), POINTER :: qmmm_env 389 390 CHARACTER(len=*), PARAMETER :: routineN = 'apply_qmmm_translate', & 391 routineP = moduleN//':'//routineN 392 393 INTEGER :: bigger_ip, i_dim, ip, max_ip, min_ip, & 394 smaller_ip, tmp_ip, unit_nr 395 INTEGER, DIMENSION(:), POINTER :: qm_atom_index 396 LOGICAL, ALLOCATABLE :: avoid(:) 397 REAL(DP) :: bigger_lat_dv, center_p(3), lat_dv, lat_dv3(3), lat_min(3), lat_p(3), & 398 max_coord_lat(3), min_coord_lat(3), smaller_lat_dv 399 REAL(DP), POINTER :: charges(:) 400 REAL(KIND=dp), DIMENSION(3) :: max_coord, min_coord, transl_v 401 TYPE(cell_type), POINTER :: mm_cell, qm_cell 402 TYPE(cp_subsys_type), POINTER :: subsys_mm, subsys_qm 403 TYPE(particle_type), DIMENSION(:), POINTER :: particles_mm, particles_qm 404 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set 405 TYPE(section_vals_type), POINTER :: subsys_section 406 407 NULLIFY (subsys_mm, subsys_qm, qm_atom_index, particles_mm, particles_qm, & 408 subsys_section, qm_cell, mm_cell, qs_kind_set) 409 410 CPASSERT(ASSOCIATED(qmmm_env)) 411 CPASSERT(qmmm_env%ref_count > 0) 412 413 CALL fist_env_get(qmmm_env%fist_env, cell=mm_cell, subsys=subsys_mm) 414 CALL get_qs_env(qmmm_env%qs_env, cell=qm_cell, cp_subsys=subsys_qm) 415 qm_atom_index => qmmm_env%qm%qm_atom_index 416 CPASSERT(ASSOCIATED(qm_atom_index)) 417 418 particles_qm => subsys_qm%particles%els 419 particles_mm => subsys_mm%particles%els 420 IF (.NOT. qmmm_env%qm%center_qm_subsys0) qmmm_env%qm%do_translate = .FALSE. 421 IF (qmmm_env%qm%do_translate) THEN 422 IF (.NOT. qmmm_env%qm%center_qm_subsys_pbc_aware) THEN 423 ! naive coordinate based min-max 424 min_coord = HUGE(0.0_dp) 425 max_coord = -HUGE(0.0_dp) 426 DO ip = 1, SIZE(qm_atom_index) 427 min_coord = MIN(min_coord, particles_mm(qm_atom_index(ip))%r) 428 max_coord = MAX(max_coord, particles_mm(qm_atom_index(ip))%r) 429 END DO 430 ELSE 431 !! periodic based min max (uses complex number based mean) 432 center_p = qmmm_pbc_aware_mean(particles_mm, mm_cell, qm_atom_index) 433 ALLOCATE (avoid(SIZE(qm_atom_index))) 434 DO i_dim = 1, 3 435 IF (mm_cell%perd(i_dim) /= 1) THEN 436 ! find absolute min and max positions (along i_dim direction) in lattice coordinates 437 min_coord_lat(i_dim) = HUGE(0.0_dp) 438 max_coord_lat(i_dim) = -HUGE(0.0_dp) 439 DO ip = 1, SIZE(qm_atom_index) 440 lat_p = MATMUL(mm_cell%h_inv, particles_mm(qm_atom_index(ip))%r) 441 min_coord_lat(i_dim) = MIN(lat_p(i_dim), min_coord_lat(i_dim)) 442 max_coord_lat(i_dim) = MAX(lat_p(i_dim), max_coord_lat(i_dim)) 443 END DO 444 ELSE 445 ! find min_ip closest to (pbc-aware) mean pos 446 avoid = .FALSE. 447 min_ip = qmmm_find_closest(particles_mm, mm_cell, qm_atom_index, avoid, center_p, i_dim, 0) 448 avoid(min_ip) = .TRUE. 449 ! find max_ip closest to min_ip 450 max_ip = qmmm_find_closest(particles_mm, mm_cell, qm_atom_index, avoid, & 451 particles_mm(qm_atom_index(min_ip))%r, i_dim, 0, lat_dv) 452 avoid(max_ip) = .TRUE. 453 ! switch min and max if necessary 454 IF (lat_dv < 0.0) THEN 455 tmp_ip = min_ip 456 min_ip = max_ip 457 max_ip = tmp_ip 458 ENDIF 459 ! loop over all other atoms 460 DO WHILE (.NOT. ALL(avoid)) 461 ! find smaller below min, bigger after max 462 smaller_ip = qmmm_find_closest(particles_mm, mm_cell, qm_atom_index, & 463 avoid, particles_mm(qm_atom_index(min_ip))%r, i_dim, -1, smaller_lat_dv) 464 bigger_ip = qmmm_find_closest(particles_mm, mm_cell, qm_atom_index, & 465 avoid, particles_mm(qm_atom_index(max_ip))%r, i_dim, 1, bigger_lat_dv) 466 ! move min or max, not both 467 IF (ABS(smaller_lat_dv) < ABS(bigger_lat_dv)) THEN 468 min_ip = smaller_ip 469 avoid(min_ip) = .TRUE. 470 ELSE 471 max_ip = bigger_ip 472 avoid(max_ip) = .TRUE. 473 END IF 474 END DO 475 ! find min and max coordinates in lattice positions (i_dim ! only) 476 lat_dv3 = qmmm_lat_dv(mm_cell, particles_mm(qm_atom_index(min_ip))%r, particles_mm(qm_atom_index(max_ip))%r) 477 IF (lat_dv3(i_dim) < 0.0_dp) lat_dv3(i_dim) = lat_dv3(i_dim) + 1.0_dp 478 lat_min = MATMUL(mm_cell%h_inv, particles_mm(qm_atom_index(min_ip))%r) 479 min_coord_lat(i_dim) = lat_min(i_dim) 480 max_coord_lat(i_dim) = lat_min(i_dim) + lat_dv3(i_dim) 481 END IF ! periodic 482 END DO ! i_dim 483 ! min and max coordinates from lattice positions to Cartesian 484 min_coord = MATMUL(mm_cell%hmat, min_coord_lat) 485 max_coord = MATMUL(mm_cell%hmat, max_coord_lat) 486 DEALLOCATE (avoid) 487 END IF ! pbc aware center 488 transl_v = (max_coord + min_coord)/2.0_dp 489 490 ! 491 ! The first time we always translate all the system in order 492 ! to centre the QM system in the box. 493 ! 494 transl_v(:) = transl_v(:) - SUM(qm_cell%hmat, 2)/2.0_dp 495 496 IF (ANY(qmmm_env%qm%utrasl /= 1.0_dp)) THEN 497 transl_v = REAL(FLOOR(transl_v/qmmm_env%qm%utrasl), KIND=dp)* & 498 qmmm_env%qm%utrasl 499 END IF 500 qmmm_env%qm%transl_v = qmmm_env%qm%transl_v + transl_v 501 particles_mm => subsys_mm%particles%els 502 DO ip = 1, subsys_mm%particles%n_els 503 particles_mm(ip)%r = particles_mm(ip)%r - transl_v 504 END DO 505 IF (qmmm_env%qm%added_shells%num_mm_atoms .GT. 0) THEN 506 DO ip = 1, qmmm_env%qm%added_shells%num_mm_atoms 507 qmmm_env%qm%added_shells%added_particles(ip)%r = qmmm_env%qm%added_shells%added_particles(ip)%r - transl_v 508 qmmm_env%qm%added_shells%added_cores(ip)%r = qmmm_env%qm%added_shells%added_cores(ip)%r - transl_v 509 END DO 510 END IF 511 unit_nr = cp_logger_get_default_io_unit() 512 IF (unit_nr > 0) WRITE (unit=unit_nr, fmt='(/1X,A)') & 513 " Translating the system in order to center the QM fragment in the QM box." 514 IF (.NOT. qmmm_env%qm%center_qm_subsys) qmmm_env%qm%do_translate = .FALSE. 515 END IF 516 particles_mm => subsys_mm%particles%els 517 DO ip = 1, SIZE(qm_atom_index) 518 particles_qm(ip)%r = particles_mm(qm_atom_index(ip))%r 519 END DO 520 521 subsys_section => section_vals_get_subs_vals(qmmm_env%qs_env%input, "SUBSYS") 522 523 CALL get_qs_env(qs_env=qmmm_env%qs_env, qs_kind_set=qs_kind_set) 524 CALL write_qs_particle_coordinates(particles_qm, qs_kind_set, subsys_section, "QM/MM first QM, then MM (0 charges)") 525 ALLOCATE (charges(SIZE(particles_mm))) 526 charges = 0.0_dp 527 CALL write_fist_particle_coordinates(particles_mm, subsys_section, charges) 528 DEALLOCATE (charges) 529 530 END SUBROUTINE apply_qmmm_translate 531 532! ************************************************************************************************** 533!> \brief pbc-aware mean QM atom position 534!> \param particles_mm ... 535!> \param mm_cell ... 536!> \param qm_atom_index ... 537!> \return ... 538! ************************************************************************************************** 539 FUNCTION qmmm_pbc_aware_mean(particles_mm, mm_cell, qm_atom_index) 540 TYPE(particle_type), DIMENSION(:), POINTER :: particles_mm 541 TYPE(cell_type), POINTER :: mm_cell 542 INTEGER, DIMENSION(:), POINTER :: qm_atom_index 543 REAL(dp) :: qmmm_pbc_aware_mean(3) 544 545 COMPLEX(dp) :: mean_z(3) 546 INTEGER :: ip 547 548 mean_z = 0.0_dp 549 DO ip = 1, SIZE(qm_atom_index) 550 mean_z = mean_z + EXP(CMPLX(0.0_dp, 1.0_dp, KIND=dp)*2.0*pi* & 551 MATMUL(mm_cell%h_inv, particles_mm(qm_atom_index(ip))%r)) 552 END DO 553 mean_z = mean_z/ABS(mean_z) 554 qmmm_pbc_aware_mean = MATMUL(mm_cell%hmat, & 555 REAL(LOG(mean_z)/(CMPLX(0.0_dp, 1.0_dp, KIND=dp)*2.0_dp*pi), dp)) 556 END FUNCTION 557 558! ************************************************************************************************** 559!> \brief minimum image lattice coordinates difference vector 560!> \param mm_cell ... 561!> \param p1 ... 562!> \param p2 ... 563!> \return ... 564! ************************************************************************************************** 565 FUNCTION qmmm_lat_dv(mm_cell, p1, p2) 566 TYPE(cell_type), POINTER :: mm_cell 567 REAL(dp) :: p1(3), p2(3), qmmm_lat_dv(3) 568 569 REAL(dp) :: lat_v1(3), lat_v2(3) 570 571 lat_v1 = MATMUL(mm_cell%h_inv, p1) 572 lat_v2 = MATMUL(mm_cell%h_inv, p2) 573 574 qmmm_lat_dv = lat_v2 - lat_v1 575 qmmm_lat_dv = qmmm_lat_dv - FLOOR(qmmm_lat_dv) 576 END FUNCTION qmmm_lat_dv 577 578! ************************************************************************************************** 579!> \brief find closest QM particle, in position/negative direction 580!> if dir is 1 or -1, respectively 581!> \param particles_mm ... 582!> \param mm_cell ... 583!> \param qm_atom_index ... 584!> \param avoid ... 585!> \param p ... 586!> \param i_dim ... 587!> \param dir ... 588!> \param closest_dv ... 589!> \return ... 590! ************************************************************************************************** 591 FUNCTION qmmm_find_closest(particles_mm, mm_cell, qm_atom_index, avoid, p, i_dim, dir, closest_dv) RESULT(closest_ip) 592 TYPE(particle_type), DIMENSION(:), POINTER :: particles_mm 593 TYPE(cell_type), POINTER :: mm_cell 594 INTEGER, DIMENSION(:), POINTER :: qm_atom_index 595 LOGICAL :: avoid(:) 596 REAL(dp) :: p(3) 597 INTEGER :: i_dim, dir 598 REAL(dp), OPTIONAL :: closest_dv 599 INTEGER :: closest_ip 600 601 INTEGER :: ip, shift 602 REAL(dp) :: lat_dv3(3), lat_dv_shifted, my_closest_dv 603 604 closest_ip = -1 605 my_closest_dv = HUGE(0.0) 606 DO ip = 1, SIZE(qm_atom_index) 607 IF (avoid(ip)) CYCLE 608 lat_dv3 = qmmm_lat_dv(mm_cell, p, particles_mm(qm_atom_index(ip))%r) 609 DO shift = -1, 1 610 lat_dv_shifted = lat_dv3(i_dim) + shift*1.0_dp 611 IF (ABS(lat_dv_shifted) < ABS(my_closest_dv) .AND. (dir*lat_dv_shifted >= 0.0)) THEN 612 my_closest_dv = lat_dv_shifted 613 closest_ip = ip 614 ENDIF 615 END DO 616 END DO 617 618 IF (PRESENT(closest_dv)) THEN 619 closest_dv = my_closest_dv 620 ENDIF 621 622 END FUNCTION qmmm_find_closest 623 624! ************************************************************************************************** 625!> \brief Computes a spherical cutoff factor for the QMMM interactions 626!> \param spherical_cutoff ... 627!> \param rij ... 628!> \param factor ... 629!> \par History 630!> 08.2008 created 631!> \author Teodoro Laino 632! ************************************************************************************************** 633 SUBROUTINE spherical_cutoff_factor(spherical_cutoff, rij, factor) 634 REAL(KIND=dp), DIMENSION(2), INTENT(IN) :: spherical_cutoff 635 REAL(KIND=dp), DIMENSION(3), INTENT(IN) :: rij 636 REAL(KIND=dp), INTENT(OUT) :: factor 637 638 CHARACTER(len=*), PARAMETER :: routineN = 'spherical_cutoff_factor', & 639 routineP = moduleN//':'//routineN 640 641 REAL(KIND=dp) :: r, r0 642 643 r = SQRT(DOT_PRODUCT(rij, rij)) 644 r0 = spherical_cutoff(1) - 20.0_dp*spherical_cutoff(2) 645 factor = 0.5_dp*(1.0_dp - TANH((r - r0)/spherical_cutoff(2))) 646 647 END SUBROUTINE spherical_cutoff_factor 648 649END MODULE qmmm_util 650