(git:58e3e09)
qmmm_util.F
Go to the documentation of this file.
1 !--------------------------------------------------------------------------------------------------!
2 ! CP2K: A general program to perform molecular dynamics simulations !
3 ! Copyright 2000-2024 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 ! **************************************************************************************************
13 MODULE qmmm_util
14  USE cell_types, ONLY: cell_type
16  USE cp_subsys_types, ONLY: cp_subsys_type
18  USE force_env_types, ONLY: force_env_type,&
19  use_qmmm,&
20  use_qmmmx
26  section_vals_type,&
28  USE kinds, ONLY: dp
29  USE mathconstants, ONLY: pi
32  USE particle_types, ONLY: particle_type
33  USE qmmm_types, ONLY: qmmm_env_type
34  USE qs_energy_types, ONLY: qs_energy_type
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, &
50 
51 CONTAINS
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)
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
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
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 
627 END MODULE qmmm_util
Handles all functions related to the CELL.
Definition: cell_types.F:15
various routines to log and control the output. The idea is that decisions about where to log should ...
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...
types that represent a subsys, i.e. a part of the system
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.
Interface for the force calculations.
integer, parameter, public use_qmmm
integer, parameter, public use_qmmmx
collects all constants needed in input so that they can be used without circular dependencies
integer, parameter, public do_qmmm_wall_quadratic
integer, parameter, public do_qmmm_wall_reflective
integer, parameter, public do_qmmm_wall_none
objects that represent the structure of input sections and the data contained in an input section
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
subroutine, public section_vals_get(section_vals, ref_count, n_repetition, n_subs_vals_rep, section, explicit)
returns various attributes about the section_vals
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
Defines the basic variable types.
Definition: kinds.F:23
integer, parameter, public dp
Definition: kinds.F:34
An array-based list which grows on demand. When the internal array is full, a new array of twice the ...
Definition: list.F:24
Definition of mathematical constants and functions.
Definition: mathconstants.F:16
real(kind=dp), parameter, public pi
Define methods related to particle_type.
subroutine, public write_qs_particle_coordinates(particle_set, qs_kind_set, subsys_section, label)
Write the atomic coordinates to the output unit.
subroutine, public write_fist_particle_coordinates(particle_set, subsys_section, charges)
Write the atomic coordinates to the output unit.
Define the data structure for the particle information.
Basic container type for QM/MM.
Definition: qmmm_types.F:12
subroutine, public apply_qmmm_unwrap(subsys_mm, subsys_qm, qm_atom_index, saved_pos)
...
Definition: qmmm_util.F:344
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
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
subroutine, public spherical_cutoff_factor(spherical_cutoff, rij, factor)
Computes a spherical cutoff factor for the QMMM interactions.
Definition: qmmm_util.F:615
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
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
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_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, 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, rhs)
Get the QUICKSTEP environment.
Define the quickstep kind type and their sub types.
Definition: qs_kind_types.F:23