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