(git:ed6f26b)
Loading...
Searching...
No Matches
qmmm_util.F
Go to the documentation of this file.
1!--------------------------------------------------------------------------------------------------!
2! CP2K: A general program to perform molecular dynamics simulations !
3! Copyright 2000-2025 CP2K developers group <https://cp2k.org> !
4! !
5! SPDX-License-Identifier: GPL-2.0-or-later !
6!--------------------------------------------------------------------------------------------------!
7
8! **************************************************************************************************
9!> \par History
10!> 09.2004 created [tlaino]
11!> \author Teodoro Laino
12! **************************************************************************************************
13MODULE qmmm_util
14 USE cell_types, ONLY: cell_type
15 USE cp_log_handling, ONLY: cp_logger_get_default_io_unit
16 USE cp_subsys_types, ONLY: cp_subsys_type
17 USE fist_environment_types, ONLY: fist_env_get
18 USE force_env_types, ONLY: force_env_type,&
19 use_qmmm,&
20 use_qmmmx
21 USE input_constants, ONLY: do_qmmm_wall_none,&
22 do_qmmm_wall_quadratic,&
23 do_qmmm_wall_reflective
24 USE input_section_types, ONLY: section_vals_get,&
25 section_vals_get_subs_vals,&
26 section_vals_type,&
27 section_vals_val_get
28 USE kinds, ONLY: dp
29 USE mathconstants, ONLY: pi
30 USE particle_methods, ONLY: write_fist_particle_coordinates,&
31 write_qs_particle_coordinates
32 USE particle_types, ONLY: particle_type
33 USE qmmm_types, ONLY: qmmm_env_type
34 USE qs_energy_types, ONLY: qs_energy_type
35 USE qs_environment_types, ONLY: get_qs_env
36 USE qs_kind_types, ONLY: qs_kind_type
37#include "./base/base_uses.f90"
38
39 IMPLICIT NONE
40 PRIVATE
41
42 LOGICAL, PRIVATE, PARAMETER :: debug_this_module = .false.
43 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qmmm_util'
44 PUBLIC :: apply_qmmm_walls_reflective, &
45 apply_qmmm_walls, &
46 apply_qmmm_translate, &
47 apply_qmmm_wrap, &
48 apply_qmmm_unwrap, &
49 spherical_cutoff_factor
50
51CONTAINS
52
53! **************************************************************************************************
54!> \brief Apply QM quadratic walls in order to avoid QM atoms escaping from
55!> the QM Box
56!> \param qmmm_env ...
57!> \par History
58!> 02.2008 created
59!> \author Benjamin G Levine
60! **************************************************************************************************
61 SUBROUTINE apply_qmmm_walls(qmmm_env)
62 TYPE(qmmm_env_type), POINTER :: qmmm_env
63
64 INTEGER :: iwall_type
65 LOGICAL :: do_qmmm_force_mixing, explicit
66 TYPE(section_vals_type), POINTER :: qmmmx_section, walls_section
67
68 walls_section => section_vals_get_subs_vals(qmmm_env%qs_env%input, "QMMM%WALLS")
69 qmmmx_section => section_vals_get_subs_vals(qmmm_env%qs_env%input, "QMMM%FORCE_MIXING")
70 CALL section_vals_get(qmmmx_section, explicit=do_qmmm_force_mixing)
71 CALL section_vals_get(walls_section, explicit=explicit)
72 IF (explicit) THEN
73 CALL section_vals_val_get(walls_section, "TYPE", i_val=iwall_type)
74 SELECT CASE (iwall_type)
75 CASE (do_qmmm_wall_quadratic)
76 IF (do_qmmm_force_mixing) THEN
77 CALL cp_warn(__location__, &
78 "Quadratic walls for QM/MM are not implemented (or useful), when "// &
79 "force mixing is active. Skipping!")
80 ELSE
81 CALL apply_qmmm_walls_quadratic(qmmm_env, walls_section)
82 END IF
83 CASE (do_qmmm_wall_reflective)
84 ! Do nothing.. reflective walls are applied directly in the integrator
85 END SELECT
86 END IF
87
88 END SUBROUTINE apply_qmmm_walls
89
90! **************************************************************************************************
91!> \brief Apply reflective QM walls in order to avoid QM atoms escaping from
92!> the QM Box
93!> \param force_env ...
94!> \par History
95!> 08.2007 created [tlaino] - Zurich University
96!> \author Teodoro Laino
97! **************************************************************************************************
98 SUBROUTINE apply_qmmm_walls_reflective(force_env)
99 TYPE(force_env_type), POINTER :: force_env
100
101 INTEGER :: ip, iwall_type, qm_index
102 INTEGER, DIMENSION(:), POINTER :: qm_atom_index
103 LOGICAL :: explicit, is_x(2), is_y(2), is_z(2)
104 REAL(kind=dp), DIMENSION(3) :: coord, qm_cell_diag, skin
105 REAL(kind=dp), DIMENSION(:), POINTER :: list
106 TYPE(cell_type), POINTER :: mm_cell, qm_cell
107 TYPE(cp_subsys_type), POINTER :: subsys_mm, subsys_qm
108 TYPE(particle_type), DIMENSION(:), POINTER :: particles_mm
109 TYPE(section_vals_type), POINTER :: walls_section
110
111 NULLIFY (subsys_mm, subsys_qm, qm_atom_index, particles_mm, qm_cell, mm_cell, &
112 walls_section)
113
114 IF (force_env%in_use /= use_qmmm .AND. force_env%in_use /= use_qmmmx) RETURN
115
116 walls_section => section_vals_get_subs_vals(force_env%root_section, "FORCE_EVAL%QMMM%WALLS")
117 CALL section_vals_get(walls_section, explicit=explicit)
118 IF (explicit) THEN
119 NULLIFY (list)
120 CALL section_vals_val_get(walls_section, "WALL_SKIN", r_vals=list)
121 CALL section_vals_val_get(walls_section, "TYPE", i_val=iwall_type)
122 skin(:) = list(:)
123 ELSE
124 ![NB]
125 iwall_type = do_qmmm_wall_reflective
126 skin(:) = 0.0_dp
127 END IF
128
129 IF (force_env%in_use == use_qmmmx) THEN
130 IF (iwall_type /= do_qmmm_wall_none) &
131 CALL cp_warn(__location__, &
132 "Reflective walls for QM/MM are not implemented (or useful) when "// &
133 "force mixing is active. Skipping!")
134 RETURN
135 END IF
136
137 ! from here on we can be sure that it's conventional QM/MM
138 cpassert(ASSOCIATED(force_env%qmmm_env))
139
140 CALL fist_env_get(force_env%qmmm_env%fist_env, cell=mm_cell, subsys=subsys_mm)
141 CALL get_qs_env(force_env%qmmm_env%qs_env, cell=qm_cell, cp_subsys=subsys_qm)
142 qm_atom_index => force_env%qmmm_env%qm%qm_atom_index
143 cpassert(ASSOCIATED(qm_atom_index))
144
145 qm_cell_diag = (/qm_cell%hmat(1, 1), &
146 qm_cell%hmat(2, 2), &
147 qm_cell%hmat(3, 3)/)
148 particles_mm => subsys_mm%particles%els
149 DO ip = 1, SIZE(qm_atom_index)
150 qm_index = qm_atom_index(ip)
151 coord = particles_mm(qm_index)%r
152 IF (any(coord < skin) .OR. any(coord > (qm_cell_diag - skin))) THEN
153 IF (explicit) THEN
154 IF (iwall_type == do_qmmm_wall_reflective) THEN
155 ! Apply Walls
156 is_x(1) = (coord(1) < skin(1))
157 is_x(2) = (coord(1) > (qm_cell_diag(1) - skin(1)))
158 is_y(1) = (coord(2) < skin(2))
159 is_y(2) = (coord(2) > (qm_cell_diag(2) - skin(2)))
160 is_z(1) = (coord(3) < skin(3))
161 is_z(2) = (coord(3) > (qm_cell_diag(3) - skin(3)))
162 IF (any(is_x)) THEN
163 ! X coordinate
164 IF (is_x(1)) THEN
165 particles_mm(qm_index)%v(1) = abs(particles_mm(qm_index)%v(1))
166 ELSE IF (is_x(2)) THEN
167 particles_mm(qm_index)%v(1) = -abs(particles_mm(qm_index)%v(1))
168 END IF
169 END IF
170 IF (any(is_y)) THEN
171 ! Y coordinate
172 IF (is_y(1)) THEN
173 particles_mm(qm_index)%v(2) = abs(particles_mm(qm_index)%v(2))
174 ELSE IF (is_y(2)) THEN
175 particles_mm(qm_index)%v(2) = -abs(particles_mm(qm_index)%v(2))
176 END IF
177 END IF
178 IF (any(is_z)) THEN
179 ! Z coordinate
180 IF (is_z(1)) THEN
181 particles_mm(qm_index)%v(3) = abs(particles_mm(qm_index)%v(3))
182 ELSE IF (is_z(2)) THEN
183 particles_mm(qm_index)%v(3) = -abs(particles_mm(qm_index)%v(3))
184 END IF
185 END IF
186 END IF
187 ELSE
188 ! Otherwise print a warning and continue crossing cp2k's finger..
189 CALL cp_warn(__location__, &
190 "One or few QM atoms are within the SKIN of the quantum box. Check your run "// &
191 "and you may possibly consider: the activation of the QMMM WALLS "// &
192 "around the QM box, switching ON the centering of the QM box or increase "// &
193 "the size of the QM cell. CP2K CONTINUE but results could be meaningless. ")
194 END IF
195 END IF
196 END DO
197
198 END SUBROUTINE apply_qmmm_walls_reflective
199
200! **************************************************************************************************
201!> \brief Apply QM quadratic walls in order to avoid QM atoms escaping from
202!> the QM Box
203!> \param qmmm_env ...
204!> \param walls_section ...
205!> \par History
206!> 02.2008 created
207!> \author Benjamin G Levine
208! **************************************************************************************************
209 SUBROUTINE apply_qmmm_walls_quadratic(qmmm_env, walls_section)
210 TYPE(qmmm_env_type), POINTER :: qmmm_env
211 TYPE(section_vals_type), POINTER :: walls_section
212
213 INTEGER :: ip, qm_index
214 INTEGER, DIMENSION(:), POINTER :: qm_atom_index
215 LOGICAL :: is_x(2), is_y(2), is_z(2)
216 REAL(kind=dp) :: k, wallenergy, wallforce
217 REAL(kind=dp), DIMENSION(3) :: coord, qm_cell_diag, skin
218 REAL(kind=dp), DIMENSION(:), POINTER :: list
219 TYPE(cell_type), POINTER :: mm_cell, qm_cell
220 TYPE(cp_subsys_type), POINTER :: subsys_mm, subsys_qm
221 TYPE(particle_type), DIMENSION(:), POINTER :: particles_mm
222 TYPE(qs_energy_type), POINTER :: energy
223
224 NULLIFY (list)
225 CALL section_vals_val_get(walls_section, "WALL_SKIN", r_vals=list)
226 CALL section_vals_val_get(walls_section, "K", r_val=k)
227 cpassert(ASSOCIATED(qmmm_env))
228
229 CALL fist_env_get(qmmm_env%fist_env, cell=mm_cell, subsys=subsys_mm)
230 CALL get_qs_env(qmmm_env%qs_env, cell=qm_cell, cp_subsys=subsys_qm)
231
232 qm_atom_index => qmmm_env%qm%qm_atom_index
233 cpassert(ASSOCIATED(qm_atom_index))
234
235 skin(:) = list(:)
236
237 qm_cell_diag = (/qm_cell%hmat(1, 1), &
238 qm_cell%hmat(2, 2), &
239 qm_cell%hmat(3, 3)/)
240 particles_mm => subsys_mm%particles%els
241 wallenergy = 0.0_dp
242 DO ip = 1, SIZE(qm_atom_index)
243 qm_index = qm_atom_index(ip)
244 coord = particles_mm(qm_index)%r
245 IF (any(coord < skin) .OR. any(coord > (qm_cell_diag - skin))) THEN
246 is_x(1) = (coord(1) < skin(1))
247 is_x(2) = (coord(1) > (qm_cell_diag(1) - skin(1)))
248 is_y(1) = (coord(2) < skin(2))
249 is_y(2) = (coord(2) > (qm_cell_diag(2) - skin(2)))
250 is_z(1) = (coord(3) < skin(3))
251 is_z(2) = (coord(3) > (qm_cell_diag(3) - skin(3)))
252 IF (is_x(1)) THEN
253 wallforce = 2.0_dp*k*(skin(1) - coord(1))
254 particles_mm(qm_index)%f(1) = particles_mm(qm_index)%f(1) + &
255 wallforce
256 wallenergy = wallenergy + wallforce*(skin(1) - coord(1))*0.5_dp
257 END IF
258 IF (is_x(2)) THEN
259 wallforce = 2.0_dp*k*((qm_cell_diag(1) - skin(1)) - coord(1))
260 particles_mm(qm_index)%f(1) = particles_mm(qm_index)%f(1) + &
261 wallforce
262 wallenergy = wallenergy + wallforce*((qm_cell_diag(1) - skin(1)) - &
263 coord(1))*0.5_dp
264 END IF
265 IF (is_y(1)) THEN
266 wallforce = 2.0_dp*k*(skin(2) - coord(2))
267 particles_mm(qm_index)%f(2) = particles_mm(qm_index)%f(2) + &
268 wallforce
269 wallenergy = wallenergy + wallforce*(skin(2) - coord(2))*0.5_dp
270 END IF
271 IF (is_y(2)) THEN
272 wallforce = 2.0_dp*k*((qm_cell_diag(2) - skin(2)) - coord(2))
273 particles_mm(qm_index)%f(2) = particles_mm(qm_index)%f(2) + &
274 wallforce
275 wallenergy = wallenergy + wallforce*((qm_cell_diag(2) - skin(2)) - &
276 coord(2))*0.5_dp
277 END IF
278 IF (is_z(1)) THEN
279 wallforce = 2.0_dp*k*(skin(3) - coord(3))
280 particles_mm(qm_index)%f(3) = particles_mm(qm_index)%f(3) + &
281 wallforce
282 wallenergy = wallenergy + wallforce*(skin(3) - coord(3))*0.5_dp
283 END IF
284 IF (is_z(2)) THEN
285 wallforce = 2.0_dp*k*((qm_cell_diag(3) - skin(3)) - coord(3))
286 particles_mm(qm_index)%f(3) = particles_mm(qm_index)%f(3) + &
287 wallforce
288 wallenergy = wallenergy + wallforce*((qm_cell_diag(3) - skin(3)) - &
289 coord(3))*0.5_dp
290 END IF
291 END IF
292 END DO
293
294 CALL get_qs_env(qs_env=qmmm_env%qs_env, energy=energy)
295 energy%total = energy%total + wallenergy
296
297 END SUBROUTINE apply_qmmm_walls_quadratic
298
299! **************************************************************************************************
300!> \brief wrap positions (with mm periodicity)
301!> \param subsys_mm ...
302!> \param mm_cell ...
303!> \param subsys_qm ...
304!> \param qm_atom_index ...
305!> \param saved_pos ...
306! **************************************************************************************************
307 SUBROUTINE apply_qmmm_wrap(subsys_mm, mm_cell, subsys_qm, qm_atom_index, saved_pos)
308 TYPE(cp_subsys_type), POINTER :: subsys_mm
309 TYPE(cell_type), POINTER :: mm_cell
310 TYPE(cp_subsys_type), OPTIONAL, POINTER :: subsys_qm
311 INTEGER, DIMENSION(:), OPTIONAL, POINTER :: qm_atom_index
312 REAL(dp), ALLOCATABLE :: saved_pos(:, :)
313
314 INTEGER :: i_dim, ip
315 REAL(dp) :: r_lat(3)
316
317 ALLOCATE (saved_pos(3, subsys_mm%particles%n_els))
318 DO ip = 1, subsys_mm%particles%n_els
319 saved_pos(1:3, ip) = subsys_mm%particles%els(ip)%r(1:3)
320 r_lat = matmul(mm_cell%h_inv, subsys_mm%particles%els(ip)%r)
321 DO i_dim = 1, 3
322 IF (mm_cell%perd(i_dim) /= 1) THEN
323 r_lat(i_dim) = 0.0_dp
324 END IF
325 END DO
326 subsys_mm%particles%els(ip)%r = subsys_mm%particles%els(ip)%r - matmul(mm_cell%hmat, floor(r_lat))
327 END DO
328
329 IF (PRESENT(subsys_qm) .AND. PRESENT(qm_atom_index)) THEN
330 DO ip = 1, SIZE(qm_atom_index)
331 subsys_qm%particles%els(ip)%r = subsys_mm%particles%els(qm_atom_index(ip))%r
332 END DO
333 END IF
334 END SUBROUTINE apply_qmmm_wrap
335
336! **************************************************************************************************
337!> \brief ...
338!> \param subsys_mm ...
339!> \param subsys_qm ...
340!> \param qm_atom_index ...
341!> \param saved_pos ...
342! **************************************************************************************************
343 SUBROUTINE apply_qmmm_unwrap(subsys_mm, subsys_qm, qm_atom_index, saved_pos)
344 TYPE(cp_subsys_type), POINTER :: subsys_mm
345 TYPE(cp_subsys_type), OPTIONAL, POINTER :: subsys_qm
346 INTEGER, DIMENSION(:), OPTIONAL, POINTER :: qm_atom_index
347 REAL(dp), ALLOCATABLE :: saved_pos(:, :)
348
349 INTEGER :: ip
350
351 DO ip = 1, subsys_mm%particles%n_els
352 subsys_mm%particles%els(ip)%r(1:3) = saved_pos(1:3, ip)
353 END DO
354
355 IF (PRESENT(subsys_qm) .AND. PRESENT(qm_atom_index)) THEN
356 DO ip = 1, SIZE(qm_atom_index)
357 subsys_qm%particles%els(ip)%r = subsys_mm%particles%els(qm_atom_index(ip))%r
358 END DO
359 END IF
360
361 DEALLOCATE (saved_pos)
362 END SUBROUTINE apply_qmmm_unwrap
363
364! **************************************************************************************************
365!> \brief Apply translation to the full system in order to center the QM
366!> system into the QM box
367!> \param qmmm_env ...
368!> \par History
369!> 08.2007 created [tlaino] - Zurich University
370!> \author Teodoro Laino
371! **************************************************************************************************
372 SUBROUTINE apply_qmmm_translate(qmmm_env)
373 TYPE(qmmm_env_type), POINTER :: qmmm_env
374
375 INTEGER :: bigger_ip, i_dim, ip, max_ip, min_ip, &
376 smaller_ip, tmp_ip, unit_nr
377 INTEGER, DIMENSION(:), POINTER :: qm_atom_index
378 LOGICAL, ALLOCATABLE :: avoid(:)
379 REAL(dp) :: bigger_lat_dv, center_p(3), lat_dv, lat_dv3(3), lat_min(3), lat_p(3), &
380 max_coord_lat(3), min_coord_lat(3), smaller_lat_dv
381 REAL(dp), POINTER :: charges(:)
382 REAL(kind=dp), DIMENSION(3) :: max_coord, min_coord, transl_v
383 TYPE(cell_type), POINTER :: mm_cell, qm_cell
384 TYPE(cp_subsys_type), POINTER :: subsys_mm, subsys_qm
385 TYPE(particle_type), DIMENSION(:), POINTER :: particles_mm, particles_qm
386 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
387 TYPE(section_vals_type), POINTER :: subsys_section
388
389 NULLIFY (subsys_mm, subsys_qm, qm_atom_index, particles_mm, particles_qm, &
390 subsys_section, qm_cell, mm_cell, qs_kind_set)
391
392 cpassert(ASSOCIATED(qmmm_env))
393
394 CALL fist_env_get(qmmm_env%fist_env, cell=mm_cell, subsys=subsys_mm)
395 CALL get_qs_env(qmmm_env%qs_env, cell=qm_cell, cp_subsys=subsys_qm)
396 qm_atom_index => qmmm_env%qm%qm_atom_index
397 cpassert(ASSOCIATED(qm_atom_index))
398
399 particles_qm => subsys_qm%particles%els
400 particles_mm => subsys_mm%particles%els
401 IF (.NOT. qmmm_env%qm%center_qm_subsys0) qmmm_env%qm%do_translate = .false.
402 IF (qmmm_env%qm%do_translate) THEN
403 IF (.NOT. qmmm_env%qm%center_qm_subsys_pbc_aware) THEN
404 ! naive coordinate based min-max
405 min_coord = huge(0.0_dp)
406 max_coord = -huge(0.0_dp)
407 DO ip = 1, SIZE(qm_atom_index)
408 min_coord = min(min_coord, particles_mm(qm_atom_index(ip))%r)
409 max_coord = max(max_coord, particles_mm(qm_atom_index(ip))%r)
410 END DO
411 ELSE
412 !! periodic based min max (uses complex number based mean)
413 center_p = qmmm_pbc_aware_mean(particles_mm, mm_cell, qm_atom_index)
414 ALLOCATE (avoid(SIZE(qm_atom_index)))
415 DO i_dim = 1, 3
416 IF (mm_cell%perd(i_dim) /= 1) THEN
417 ! find absolute min and max positions (along i_dim direction) in lattice coordinates
418 min_coord_lat(i_dim) = huge(0.0_dp)
419 max_coord_lat(i_dim) = -huge(0.0_dp)
420 DO ip = 1, SIZE(qm_atom_index)
421 lat_p = matmul(mm_cell%h_inv, particles_mm(qm_atom_index(ip))%r)
422 min_coord_lat(i_dim) = min(lat_p(i_dim), min_coord_lat(i_dim))
423 max_coord_lat(i_dim) = max(lat_p(i_dim), max_coord_lat(i_dim))
424 END DO
425 ELSE
426 ! find min_ip closest to (pbc-aware) mean pos
427 avoid = .false.
428 min_ip = qmmm_find_closest(particles_mm, mm_cell, qm_atom_index, avoid, center_p, i_dim, 0)
429 avoid(min_ip) = .true.
430 ! find max_ip closest to min_ip
431 max_ip = qmmm_find_closest(particles_mm, mm_cell, qm_atom_index, avoid, &
432 particles_mm(qm_atom_index(min_ip))%r, i_dim, 0, lat_dv)
433 avoid(max_ip) = .true.
434 ! switch min and max if necessary
435 IF (lat_dv < 0.0) THEN
436 tmp_ip = min_ip
437 min_ip = max_ip
438 max_ip = tmp_ip
439 END IF
440 ! loop over all other atoms
441 DO WHILE (.NOT. all(avoid))
442 ! find smaller below min, bigger after max
443 smaller_ip = qmmm_find_closest(particles_mm, mm_cell, qm_atom_index, &
444 avoid, particles_mm(qm_atom_index(min_ip))%r, i_dim, -1, smaller_lat_dv)
445 bigger_ip = qmmm_find_closest(particles_mm, mm_cell, qm_atom_index, &
446 avoid, particles_mm(qm_atom_index(max_ip))%r, i_dim, 1, bigger_lat_dv)
447 ! move min or max, not both
448 IF (abs(smaller_lat_dv) < abs(bigger_lat_dv)) THEN
449 min_ip = smaller_ip
450 avoid(min_ip) = .true.
451 ELSE
452 max_ip = bigger_ip
453 avoid(max_ip) = .true.
454 END IF
455 END DO
456 ! find min and max coordinates in lattice positions (i_dim ! only)
457 lat_dv3 = qmmm_lat_dv(mm_cell, particles_mm(qm_atom_index(min_ip))%r, particles_mm(qm_atom_index(max_ip))%r)
458 IF (lat_dv3(i_dim) < 0.0_dp) lat_dv3(i_dim) = lat_dv3(i_dim) + 1.0_dp
459 lat_min = matmul(mm_cell%h_inv, particles_mm(qm_atom_index(min_ip))%r)
460 min_coord_lat(i_dim) = lat_min(i_dim)
461 max_coord_lat(i_dim) = lat_min(i_dim) + lat_dv3(i_dim)
462 END IF ! periodic
463 END DO ! i_dim
464 ! min and max coordinates from lattice positions to Cartesian
465 min_coord = matmul(mm_cell%hmat, min_coord_lat)
466 max_coord = matmul(mm_cell%hmat, max_coord_lat)
467 DEALLOCATE (avoid)
468 END IF ! pbc aware center
469 transl_v = (max_coord + min_coord)/2.0_dp
470
471 !
472 ! The first time we always translate all the system in order
473 ! to centre the QM system in the box.
474 !
475 transl_v(:) = transl_v(:) - sum(qm_cell%hmat, 2)/2.0_dp
476
477 IF (any(qmmm_env%qm%utrasl /= 1.0_dp)) THEN
478 transl_v = real(floor(transl_v/qmmm_env%qm%utrasl), kind=dp)* &
479 qmmm_env%qm%utrasl
480 END IF
481 qmmm_env%qm%transl_v = qmmm_env%qm%transl_v + transl_v
482 particles_mm => subsys_mm%particles%els
483 DO ip = 1, subsys_mm%particles%n_els
484 particles_mm(ip)%r = particles_mm(ip)%r - transl_v
485 END DO
486 IF (qmmm_env%qm%added_shells%num_mm_atoms .GT. 0) THEN
487 DO ip = 1, qmmm_env%qm%added_shells%num_mm_atoms
488 qmmm_env%qm%added_shells%added_particles(ip)%r = qmmm_env%qm%added_shells%added_particles(ip)%r - transl_v
489 qmmm_env%qm%added_shells%added_cores(ip)%r = qmmm_env%qm%added_shells%added_cores(ip)%r - transl_v
490 END DO
491 END IF
492 unit_nr = cp_logger_get_default_io_unit()
493 IF (unit_nr > 0) WRITE (unit=unit_nr, fmt='(/1X,A)') &
494 " Translating the system in order to center the QM fragment in the QM box."
495 IF (.NOT. qmmm_env%qm%center_qm_subsys) qmmm_env%qm%do_translate = .false.
496 END IF
497 particles_mm => subsys_mm%particles%els
498 DO ip = 1, SIZE(qm_atom_index)
499 particles_qm(ip)%r = particles_mm(qm_atom_index(ip))%r
500 END DO
501
502 subsys_section => section_vals_get_subs_vals(qmmm_env%qs_env%input, "SUBSYS")
503
504 CALL get_qs_env(qs_env=qmmm_env%qs_env, qs_kind_set=qs_kind_set)
505 CALL write_qs_particle_coordinates(particles_qm, qs_kind_set, subsys_section, "QM/MM first QM, then MM (0 charges)")
506 ALLOCATE (charges(SIZE(particles_mm)))
507 charges = 0.0_dp
508 CALL write_fist_particle_coordinates(particles_mm, subsys_section, charges)
509 DEALLOCATE (charges)
510
511 END SUBROUTINE apply_qmmm_translate
512
513! **************************************************************************************************
514!> \brief pbc-aware mean QM atom position
515!> \param particles_mm ...
516!> \param mm_cell ...
517!> \param qm_atom_index ...
518!> \return ...
519! **************************************************************************************************
520 FUNCTION qmmm_pbc_aware_mean(particles_mm, mm_cell, qm_atom_index)
521 TYPE(particle_type), DIMENSION(:), POINTER :: particles_mm
522 TYPE(cell_type), POINTER :: mm_cell
523 INTEGER, DIMENSION(:), POINTER :: qm_atom_index
524 REAL(dp) :: qmmm_pbc_aware_mean(3)
525
526 COMPLEX(dp) :: mean_z(3)
527 INTEGER :: ip
528
529 mean_z = 0.0_dp
530 DO ip = 1, SIZE(qm_atom_index)
531 mean_z = mean_z + exp(cmplx(0.0_dp, 1.0_dp, kind=dp)*2.0*pi* &
532 matmul(mm_cell%h_inv, particles_mm(qm_atom_index(ip))%r))
533 END DO
534 mean_z = mean_z/abs(mean_z)
535 qmmm_pbc_aware_mean = matmul(mm_cell%hmat, &
536 REAL(log(mean_z)/(cmplx(0.0_dp, 1.0_dp, kind=dp)*2.0_dp*pi), dp))
537 END FUNCTION
538
539! **************************************************************************************************
540!> \brief minimum image lattice coordinates difference vector
541!> \param mm_cell ...
542!> \param p1 ...
543!> \param p2 ...
544!> \return ...
545! **************************************************************************************************
546 FUNCTION qmmm_lat_dv(mm_cell, p1, p2)
547 TYPE(cell_type), POINTER :: mm_cell
548 REAL(dp) :: p1(3), p2(3), qmmm_lat_dv(3)
549
550 REAL(dp) :: lat_v1(3), lat_v2(3)
551
552 lat_v1 = matmul(mm_cell%h_inv, p1)
553 lat_v2 = matmul(mm_cell%h_inv, p2)
554
555 qmmm_lat_dv = lat_v2 - lat_v1
556 qmmm_lat_dv = qmmm_lat_dv - floor(qmmm_lat_dv)
557 END FUNCTION qmmm_lat_dv
558
559! **************************************************************************************************
560!> \brief find closest QM particle, in position/negative direction
561!> if dir is 1 or -1, respectively
562!> \param particles_mm ...
563!> \param mm_cell ...
564!> \param qm_atom_index ...
565!> \param avoid ...
566!> \param p ...
567!> \param i_dim ...
568!> \param dir ...
569!> \param closest_dv ...
570!> \return ...
571! **************************************************************************************************
572 FUNCTION qmmm_find_closest(particles_mm, mm_cell, qm_atom_index, avoid, p, i_dim, dir, closest_dv) RESULT(closest_ip)
573 TYPE(particle_type), DIMENSION(:), POINTER :: particles_mm
574 TYPE(cell_type), POINTER :: mm_cell
575 INTEGER, DIMENSION(:), POINTER :: qm_atom_index
576 LOGICAL :: avoid(:)
577 REAL(dp) :: p(3)
578 INTEGER :: i_dim, dir
579 REAL(dp), OPTIONAL :: closest_dv
580 INTEGER :: closest_ip
581
582 INTEGER :: ip, shift
583 REAL(dp) :: lat_dv3(3), lat_dv_shifted, my_closest_dv
584
585 closest_ip = -1
586 my_closest_dv = huge(0.0)
587 DO ip = 1, SIZE(qm_atom_index)
588 IF (avoid(ip)) cycle
589 lat_dv3 = qmmm_lat_dv(mm_cell, p, particles_mm(qm_atom_index(ip))%r)
590 DO shift = -1, 1
591 lat_dv_shifted = lat_dv3(i_dim) + shift*1.0_dp
592 IF (abs(lat_dv_shifted) < abs(my_closest_dv) .AND. (dir*lat_dv_shifted >= 0.0)) THEN
593 my_closest_dv = lat_dv_shifted
594 closest_ip = ip
595 END IF
596 END DO
597 END DO
598
599 IF (PRESENT(closest_dv)) THEN
600 closest_dv = my_closest_dv
601 END IF
602
603 END FUNCTION qmmm_find_closest
604
605! **************************************************************************************************
606!> \brief Computes a spherical cutoff factor for the QMMM interactions
607!> \param spherical_cutoff ...
608!> \param rij ...
609!> \param factor ...
610!> \par History
611!> 08.2008 created
612!> \author Teodoro Laino
613! **************************************************************************************************
614 SUBROUTINE spherical_cutoff_factor(spherical_cutoff, rij, factor)
615 REAL(kind=dp), DIMENSION(2), INTENT(IN) :: spherical_cutoff
616 REAL(kind=dp), DIMENSION(3), INTENT(IN) :: rij
617 REAL(kind=dp), INTENT(OUT) :: factor
618
619 REAL(kind=dp) :: r, r0
620
621 r = sqrt(dot_product(rij, rij))
622 r0 = spherical_cutoff(1) - 20.0_dp*spherical_cutoff(2)
623 factor = 0.5_dp*(1.0_dp - tanh((r - r0)/spherical_cutoff(2)))
624
625 END SUBROUTINE spherical_cutoff_factor
626
627END MODULE qmmm_util
cell_types
Handles all functions related to the CELL.
Definition cell_types.F:15
cp_log_handling
various routines to log and control the output. The idea is that decisions about where to log should ...
Definition cp_log_handling.F:41
cp_log_handling::cp_logger_get_default_io_unit
integer function, public cp_logger_get_default_io_unit(logger)
returns the unit nr for the ionode (-1 on all other processors) skips as well checks if the procs cal...
Definition cp_log_handling.F:515
cp_subsys_types
types that represent a subsys, i.e. a part of the system
Definition cp_subsys_types.F:16
fist_environment_types
Definition fist_environment_types.F:15
fist_environment_types::fist_env_get
subroutine, public fist_env_get(fist_env, atomic_kind_set, particle_set, ewald_pw, local_particles, local_molecules, molecule_kind_set, molecule_set, cell, cell_ref, ewald_env, fist_nonbond_env, thermo, para_env, subsys, qmmm, qmmm_env, input, shell_model, shell_model_ad, shell_particle_set, core_particle_set, multipoles, results, exclusions, efield)
Purpose: Get the FIST environment.
Definition fist_environment_types.F:140
force_env_types
Interface for the force calculations.
Definition force_env_types.F:18
force_env_types::use_qmmm
integer, parameter, public use_qmmm
Definition force_env_types.F:83
force_env_types::use_qmmmx
integer, parameter, public use_qmmmx
Definition force_env_types.F:83
input_constants
collects all constants needed in input so that they can be used without circular dependencies
Definition input_constants.F:17
input_constants::do_qmmm_wall_quadratic
integer, parameter, public do_qmmm_wall_quadratic
Definition input_constants.F:385
input_constants::do_qmmm_wall_reflective
integer, parameter, public do_qmmm_wall_reflective
Definition input_constants.F:385
input_constants::do_qmmm_wall_none
integer, parameter, public do_qmmm_wall_none
Definition input_constants.F:385
input_section_types
objects that represent the structure of input sections and the data contained in an input section
Definition input_section_types.F:15
input_section_types::section_vals_get_subs_vals
recursive type(section_vals_type) function, pointer, public section_vals_get_subs_vals(section_vals, subsection_name, i_rep_section, can_return_null)
returns the values of the requested subsection
Definition input_section_types.F:731
input_section_types::section_vals_get
subroutine, public section_vals_get(section_vals, ref_count, n_repetition, n_subs_vals_rep, section, explicit)
returns various attributes about the section_vals
Definition input_section_types.F:704
input_section_types::section_vals_val_get
subroutine, public section_vals_val_get(section_vals, keyword_name, i_rep_section, i_rep_val, n_rep_val, val, l_val, i_val, r_val, c_val, l_vals, i_vals, r_vals, c_vals, explicit)
returns the requested value
Definition input_section_types.F:1047
kinds
Defines the basic variable types.
Definition kinds.F:23
kinds::dp
integer, parameter, public dp
Definition kinds.F:34
list
An array-based list which grows on demand. When the internal array is full, a new array of twice the ...
Definition list.F:24
mathconstants
Definition of mathematical constants and functions.
Definition mathconstants.F:16
mathconstants::pi
real(kind=dp), parameter, public pi
Definition mathconstants.F:110
particle_methods
Define methods related to particle_type.
Definition particle_methods.F:14
particle_methods::write_qs_particle_coordinates
subroutine, public write_qs_particle_coordinates(particle_set, qs_kind_set, subsys_section, label)
Write the atomic coordinates to the output unit.
Definition particle_methods.F:546
particle_methods::write_fist_particle_coordinates
subroutine, public write_fist_particle_coordinates(particle_set, subsys_section, charges)
Write the atomic coordinates to the output unit.
Definition particle_methods.F:484
particle_types
Define the data structure for the particle information.
Definition particle_types.F:19
qmmm_types
Basic container type for QM/MM.
Definition qmmm_types.F:12
qmmm_util
Definition qmmm_util.F:13
qmmm_util::apply_qmmm_unwrap
subroutine, public apply_qmmm_unwrap(subsys_mm, subsys_qm, qm_atom_index, saved_pos)
...
Definition qmmm_util.F:344
qmmm_util::apply_qmmm_walls_reflective
subroutine, public apply_qmmm_walls_reflective(force_env)
Apply reflective QM walls in order to avoid QM atoms escaping from the QM Box.
Definition qmmm_util.F:99
qmmm_util::apply_qmmm_walls
subroutine, public apply_qmmm_walls(qmmm_env)
Apply QM quadratic walls in order to avoid QM atoms escaping from the QM Box.
Definition qmmm_util.F:62
qmmm_util::spherical_cutoff_factor
subroutine, public spherical_cutoff_factor(spherical_cutoff, rij, factor)
Computes a spherical cutoff factor for the QMMM interactions.
Definition qmmm_util.F:615
qmmm_util::apply_qmmm_translate
subroutine, public apply_qmmm_translate(qmmm_env)
Apply translation to the full system in order to center the QM system into the QM box.
Definition qmmm_util.F:373
qmmm_util::apply_qmmm_wrap
subroutine, public apply_qmmm_wrap(subsys_mm, mm_cell, subsys_qm, qm_atom_index, saved_pos)
wrap positions (with mm periodicity)
Definition qmmm_util.F:308
qs_energy_types
Definition qs_energy_types.F:14
qs_environment_types
Definition qs_environment_types.F:14
qs_environment_types::get_qs_env
subroutine, public get_qs_env(qs_env, atomic_kind_set, qs_kind_set, cell, super_cell, cell_ref, use_ref_cell, kpoints, dft_control, mos, sab_orb, sab_all, qmmm, qmmm_periodic, sac_ae, sac_ppl, sac_lri, sap_ppnl, sab_vdw, sab_scp, sap_oce, sab_lrc, sab_se, sab_xtbe, sab_tbe, sab_core, sab_xb, sab_xtb_pp, sab_xtb_nonbond, sab_almo, sab_kp, sab_kp_nosym, particle_set, energy, force, matrix_h, matrix_h_im, matrix_ks, matrix_ks_im, matrix_vxc, run_rtp, rtp, matrix_h_kp, matrix_h_im_kp, matrix_ks_kp, matrix_ks_im_kp, matrix_vxc_kp, kinetic_kp, matrix_s_kp, matrix_w_kp, matrix_s_ri_aux_kp, matrix_s, matrix_s_ri_aux, matrix_w, matrix_p_mp2, matrix_p_mp2_admm, rho, rho_xc, pw_env, ewald_env, ewald_pw, active_space, mpools, input, para_env, blacs_env, scf_control, rel_control, kinetic, qs_charges, vppl, rho_core, rho_nlcc, rho_nlcc_g, ks_env, ks_qmmm_env, wf_history, scf_env, local_particles, local_molecules, distribution_2d, dbcsr_dist, molecule_kind_set, molecule_set, subsys, cp_subsys, oce, local_rho_set, rho_atom_set, task_list, task_list_soft, rho0_atom_set, rho0_mpole, rhoz_set, ecoul_1c, rho0_s_rs, rho0_s_gs, do_kpoints, has_unit_metric, requires_mo_derivs, mo_derivs, mo_loc_history, nkind, natom, nelectron_total, nelectron_spin, efield, neighbor_list_id, linres_control, xas_env, virial, cp_ddapc_env, cp_ddapc_ewald, outer_scf_history, outer_scf_ihistory, x_data, et_coupling, dftb_potential, results, se_taper, se_store_int_env, se_nddo_mpole, se_nonbond_env, admm_env, lri_env, lri_density, exstate_env, ec_env, harris_env, dispersion_env, gcp_env, vee, rho_external, external_vxc, mask, mp2_env, bs_env, kg_env, wanniercentres, atprop, ls_scf_env, do_transport, transport_env, v_hartree_rspace, s_mstruct_changed, rho_changed, potential_changed, forces_up_to_date, mscfg_env, almo_scf_env, gradient_history, variable_history, embed_pot, spin_embed_pot, polar_env, mos_last_converged, eeq, rhs)
Get the QUICKSTEP environment.
Definition qs_environment_types.F:514
qs_kind_types
Define the quickstep kind type and their sub types.
Definition qs_kind_types.F:23
cell_types::cell_type
Type defining parameters related to the simulation cell.
Definition cell_types.F:55
cp_subsys_types::cp_subsys_type
represents a system: atoms, molecules, their pos,vel,...
Definition cp_subsys_types.F:89
force_env_types::force_env_type
wrapper to abstract the force evaluation of the various methods
Definition force_env_types.F:145
input_section_types::section_vals_type
stores the values of a section
Definition input_section_types.F:127
particle_types::particle_type
Definition particle_types.F:35
qmmm_types::qmmm_env_type
Definition qmmm_types.F:34
qs_energy_types::qs_energy_type
Definition qs_energy_types.F:25
qs_kind_types::qs_kind_type
Provides all information about a quickstep kind.
Definition qs_kind_types.F:171