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