Line data Source code
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 : ! **************************************************************************************************
13 : MODULE 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 :
52 : CONTAINS
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 11406 : 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 3802 : walls_section => section_vals_get_subs_vals(qmmm_env%qs_env%input, "QMMM%WALLS")
70 3802 : qmmmx_section => section_vals_get_subs_vals(qmmm_env%qs_env%input, "QMMM%FORCE_MIXING")
71 3802 : CALL section_vals_get(qmmmx_section, explicit=do_qmmm_force_mixing)
72 3802 : CALL section_vals_get(walls_section, explicit=explicit)
73 3802 : IF (explicit) THEN
74 404 : CALL section_vals_val_get(walls_section, "TYPE", i_val=iwall_type)
75 202 : SELECT CASE (iwall_type)
76 : CASE (do_qmmm_wall_quadratic)
77 404 : IF (do_qmmm_force_mixing) THEN
78 : CALL cp_warn(__LOCATION__, &
79 : "Quadratic walls for QM/MM are not implemented (or useful), when "// &
80 0 : "force mixing is active. Skipping!")
81 : ELSE
82 202 : 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 3802 : 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 42259 : SUBROUTINE apply_qmmm_walls_reflective(force_env)
100 : TYPE(force_env_type), POINTER :: force_env
101 :
102 : INTEGER :: ip, iwall_type, qm_index
103 40941 : 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 40941 : 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 40941 : TYPE(particle_type), DIMENSION(:), POINTER :: particles_mm
110 : TYPE(section_vals_type), POINTER :: walls_section
111 :
112 40941 : NULLIFY (subsys_mm, subsys_qm, qm_atom_index, particles_mm, qm_cell, mm_cell, &
113 40941 : walls_section)
114 :
115 39671 : IF (force_env%in_use /= use_qmmm .AND. force_env%in_use /= use_qmmmx) RETURN
116 :
117 1318 : walls_section => section_vals_get_subs_vals(force_env%root_section, "FORCE_EVAL%QMMM%WALLS")
118 1318 : CALL section_vals_get(walls_section, explicit=explicit)
119 1318 : IF (explicit) THEN
120 400 : NULLIFY (list)
121 400 : CALL section_vals_val_get(walls_section, "WALL_SKIN", r_vals=list)
122 400 : CALL section_vals_val_get(walls_section, "TYPE", i_val=iwall_type)
123 1600 : skin(:) = list(:)
124 : ELSE
125 : ![NB]
126 918 : iwall_type = do_qmmm_wall_reflective
127 918 : skin(:) = 0.0_dp
128 : END IF
129 :
130 1318 : IF (force_env%in_use == use_qmmmx) THEN
131 48 : 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 48 : "force mixing is active. Skipping!")
135 : END IF
136 48 : RETURN
137 : END IF
138 :
139 : ! from here on we can be sure that it's conventional QM/MM
140 1270 : CPASSERT(ASSOCIATED(force_env%qmmm_env))
141 :
142 1270 : CALL fist_env_get(force_env%qmmm_env%fist_env, cell=mm_cell, subsys=subsys_mm)
143 1270 : CALL get_qs_env(force_env%qmmm_env%qs_env, cell=qm_cell, cp_subsys=subsys_qm)
144 1270 : qm_atom_index => force_env%qmmm_env%qm%qm_atom_index
145 1270 : CPASSERT(ASSOCIATED(qm_atom_index))
146 :
147 : qm_cell_diag = [qm_cell%hmat(1, 1), &
148 : qm_cell%hmat(2, 2), &
149 5080 : qm_cell%hmat(3, 3)]
150 1270 : particles_mm => subsys_mm%particles%els
151 7120 : DO ip = 1, SIZE(qm_atom_index)
152 5850 : qm_index = qm_atom_index(ip)
153 23400 : coord = particles_mm(qm_index)%r
154 48034 : IF (ANY(coord < skin) .OR. ANY(coord > (qm_cell_diag - skin))) THEN
155 12 : IF (explicit) THEN
156 12 : IF (iwall_type == do_qmmm_wall_reflective) THEN
157 : ! Apply Walls
158 2 : is_x(1) = (coord(1) < skin(1))
159 2 : is_x(2) = (coord(1) > (qm_cell_diag(1) - skin(1)))
160 2 : is_y(1) = (coord(2) < skin(2))
161 2 : is_y(2) = (coord(2) > (qm_cell_diag(2) - skin(2)))
162 2 : is_z(1) = (coord(3) < skin(3))
163 2 : is_z(2) = (coord(3) > (qm_cell_diag(3) - skin(3)))
164 2 : IF (ANY(is_x)) THEN
165 : ! X coordinate
166 2 : IF (is_x(1)) THEN
167 2 : particles_mm(qm_index)%v(1) = ABS(particles_mm(qm_index)%v(1))
168 0 : ELSE IF (is_x(2)) THEN
169 0 : particles_mm(qm_index)%v(1) = -ABS(particles_mm(qm_index)%v(1))
170 : END IF
171 : END IF
172 6 : IF (ANY(is_y)) THEN
173 : ! Y coordinate
174 0 : IF (is_y(1)) THEN
175 0 : particles_mm(qm_index)%v(2) = ABS(particles_mm(qm_index)%v(2))
176 0 : ELSE IF (is_y(2)) THEN
177 0 : particles_mm(qm_index)%v(2) = -ABS(particles_mm(qm_index)%v(2))
178 : END IF
179 : END IF
180 6 : IF (ANY(is_z)) THEN
181 : ! Z coordinate
182 0 : IF (is_z(1)) THEN
183 0 : particles_mm(qm_index)%v(3) = ABS(particles_mm(qm_index)%v(3))
184 0 : ELSE IF (is_z(2)) THEN
185 0 : 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 0 : "the size of the QM cell. CP2K CONTINUE but results could be meaningless. ")
196 : END IF
197 : END IF
198 : END DO
199 :
200 40941 : 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 404 : 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 202 : 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 202 : 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 202 : TYPE(particle_type), DIMENSION(:), POINTER :: particles_mm
224 : TYPE(qs_energy_type), POINTER :: energy
225 :
226 202 : NULLIFY (list)
227 202 : CALL section_vals_val_get(walls_section, "WALL_SKIN", r_vals=list)
228 202 : CALL section_vals_val_get(walls_section, "K", r_val=k)
229 202 : CPASSERT(ASSOCIATED(qmmm_env))
230 :
231 202 : CALL fist_env_get(qmmm_env%fist_env, cell=mm_cell, subsys=subsys_mm)
232 202 : CALL get_qs_env(qmmm_env%qs_env, cell=qm_cell, cp_subsys=subsys_qm)
233 :
234 202 : qm_atom_index => qmmm_env%qm%qm_atom_index
235 202 : CPASSERT(ASSOCIATED(qm_atom_index))
236 :
237 808 : skin(:) = list(:)
238 :
239 : qm_cell_diag = [qm_cell%hmat(1, 1), &
240 : qm_cell%hmat(2, 2), &
241 808 : qm_cell%hmat(3, 3)]
242 202 : particles_mm => subsys_mm%particles%els
243 202 : wallenergy = 0.0_dp
244 808 : DO ip = 1, SIZE(qm_atom_index)
245 606 : qm_index = qm_atom_index(ip)
246 2424 : coord = particles_mm(qm_index)%r
247 5014 : IF (ANY(coord < skin) .OR. ANY(coord > (qm_cell_diag - skin))) THEN
248 12 : is_x(1) = (coord(1) < skin(1))
249 12 : is_x(2) = (coord(1) > (qm_cell_diag(1) - skin(1)))
250 12 : is_y(1) = (coord(2) < skin(2))
251 12 : is_y(2) = (coord(2) > (qm_cell_diag(2) - skin(2)))
252 12 : is_z(1) = (coord(3) < skin(3))
253 12 : is_z(2) = (coord(3) > (qm_cell_diag(3) - skin(3)))
254 12 : IF (is_x(1)) THEN
255 12 : wallforce = 2.0_dp*k*(skin(1) - coord(1))
256 : particles_mm(qm_index)%f(1) = particles_mm(qm_index)%f(1) + &
257 12 : wallforce
258 12 : wallenergy = wallenergy + wallforce*(skin(1) - coord(1))*0.5_dp
259 : END IF
260 12 : IF (is_x(2)) THEN
261 0 : 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 0 : wallforce
264 : wallenergy = wallenergy + wallforce*((qm_cell_diag(1) - skin(1)) - &
265 0 : coord(1))*0.5_dp
266 : END IF
267 12 : IF (is_y(1)) THEN
268 0 : wallforce = 2.0_dp*k*(skin(2) - coord(2))
269 : particles_mm(qm_index)%f(2) = particles_mm(qm_index)%f(2) + &
270 0 : wallforce
271 0 : wallenergy = wallenergy + wallforce*(skin(2) - coord(2))*0.5_dp
272 : END IF
273 12 : IF (is_y(2)) THEN
274 0 : 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 0 : wallforce
277 : wallenergy = wallenergy + wallforce*((qm_cell_diag(2) - skin(2)) - &
278 0 : coord(2))*0.5_dp
279 : END IF
280 12 : IF (is_z(1)) THEN
281 0 : wallforce = 2.0_dp*k*(skin(3) - coord(3))
282 : particles_mm(qm_index)%f(3) = particles_mm(qm_index)%f(3) + &
283 0 : wallforce
284 0 : wallenergy = wallenergy + wallforce*(skin(3) - coord(3))*0.5_dp
285 : END IF
286 12 : IF (is_z(2)) THEN
287 0 : 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 0 : wallforce
290 : wallenergy = wallenergy + wallforce*((qm_cell_diag(3) - skin(3)) - &
291 0 : coord(3))*0.5_dp
292 : END IF
293 : END IF
294 : END DO
295 :
296 202 : CALL get_qs_env(qs_env=qmmm_env%qs_env, energy=energy)
297 202 : energy%total = energy%total + wallenergy
298 :
299 202 : 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 104 : 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 312 : ALLOCATE (saved_pos(3, subsys_mm%particles%n_els))
320 199676 : DO ip = 1, subsys_mm%particles%n_els
321 798288 : saved_pos(1:3, ip) = subsys_mm%particles%els(ip)%r(1:3)
322 2594436 : r_lat = MATMUL(mm_cell%h_inv, subsys_mm%particles%els(ip)%r)
323 798288 : DO i_dim = 1, 3
324 798288 : IF (mm_cell%perd(i_dim) /= 1) THEN
325 0 : r_lat(i_dim) = 0.0_dp
326 : END IF
327 : END DO
328 3791972 : subsys_mm%particles%els(ip)%r = subsys_mm%particles%els(ip)%r - MATMUL(mm_cell%hmat, FLOOR(r_lat))
329 : END DO
330 :
331 104 : IF (PRESENT(subsys_qm) .AND. PRESENT(qm_atom_index)) THEN
332 2444 : DO ip = 1, SIZE(qm_atom_index)
333 18824 : subsys_qm%particles%els(ip)%r = subsys_mm%particles%els(qm_atom_index(ip))%r
334 : END DO
335 : END IF
336 104 : 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 104 : 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 199676 : DO ip = 1, subsys_mm%particles%n_els
354 798392 : subsys_mm%particles%els(ip)%r(1:3) = saved_pos(1:3, ip)
355 : END DO
356 :
357 104 : IF (PRESENT(subsys_qm) .AND. PRESENT(qm_atom_index)) THEN
358 2444 : DO ip = 1, SIZE(qm_atom_index)
359 18824 : subsys_qm%particles%els(ip)%r = subsys_mm%particles%els(qm_atom_index(ip))%r
360 : END DO
361 : END IF
362 :
363 104 : DEALLOCATE (saved_pos)
364 104 : 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 3914 : 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 3914 : 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 3914 : 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 3914 : TYPE(particle_type), DIMENSION(:), POINTER :: particles_mm, particles_qm
388 3914 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
389 : TYPE(section_vals_type), POINTER :: subsys_section
390 :
391 3914 : NULLIFY (subsys_mm, subsys_qm, qm_atom_index, particles_mm, particles_qm, &
392 3914 : subsys_section, qm_cell, mm_cell, qs_kind_set)
393 :
394 0 : CPASSERT(ASSOCIATED(qmmm_env))
395 :
396 3914 : CALL fist_env_get(qmmm_env%fist_env, cell=mm_cell, subsys=subsys_mm)
397 3914 : CALL get_qs_env(qmmm_env%qs_env, cell=qm_cell, cp_subsys=subsys_qm)
398 3914 : qm_atom_index => qmmm_env%qm%qm_atom_index
399 3914 : CPASSERT(ASSOCIATED(qm_atom_index))
400 :
401 3914 : particles_qm => subsys_qm%particles%els
402 3914 : particles_mm => subsys_mm%particles%els
403 3914 : IF (.NOT. qmmm_env%qm%center_qm_subsys0) qmmm_env%qm%do_translate = .FALSE.
404 3914 : IF (qmmm_env%qm%do_translate) THEN
405 972 : IF (.NOT. qmmm_env%qm%center_qm_subsys_pbc_aware) THEN
406 : ! naive coordinate based min-max
407 3792 : min_coord = HUGE(0.0_dp)
408 3792 : max_coord = -HUGE(0.0_dp)
409 7996 : DO ip = 1, SIZE(qm_atom_index)
410 28192 : min_coord = MIN(min_coord, particles_mm(qm_atom_index(ip))%r)
411 29140 : 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 24 : center_p = qmmm_pbc_aware_mean(particles_mm, mm_cell, qm_atom_index)
416 72 : ALLOCATE (avoid(SIZE(qm_atom_index)))
417 96 : DO i_dim = 1, 3
418 96 : IF (mm_cell%perd(i_dim) /= 1) THEN
419 : ! find absolute min and max positions (along i_dim direction) in lattice coordinates
420 0 : min_coord_lat(i_dim) = HUGE(0.0_dp)
421 0 : max_coord_lat(i_dim) = -HUGE(0.0_dp)
422 0 : DO ip = 1, SIZE(qm_atom_index)
423 0 : lat_p = MATMUL(mm_cell%h_inv, particles_mm(qm_atom_index(ip))%r)
424 0 : min_coord_lat(i_dim) = MIN(lat_p(i_dim), min_coord_lat(i_dim))
425 0 : 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 72 : avoid = .FALSE.
430 72 : min_ip = qmmm_find_closest(particles_mm, mm_cell, qm_atom_index, avoid, center_p, i_dim, 0)
431 72 : 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 72 : particles_mm(qm_atom_index(min_ip))%r, i_dim, 0, lat_dv)
435 72 : avoid(max_ip) = .TRUE.
436 : ! switch min and max if necessary
437 72 : IF (lat_dv < 0.0) THEN
438 0 : tmp_ip = min_ip
439 0 : min_ip = max_ip
440 0 : max_ip = tmp_ip
441 : END IF
442 : ! loop over all other atoms
443 2726 : 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 612 : 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 612 : avoid, particles_mm(qm_atom_index(max_ip))%r, i_dim, 1, bigger_lat_dv)
449 : ! move min or max, not both
450 684 : IF (ABS(smaller_lat_dv) < ABS(bigger_lat_dv)) THEN
451 180 : min_ip = smaller_ip
452 180 : avoid(min_ip) = .TRUE.
453 : ELSE
454 432 : max_ip = bigger_ip
455 432 : avoid(max_ip) = .TRUE.
456 : END IF
457 : END DO
458 : ! find min and max coordinates in lattice positions (i_dim ! only)
459 72 : lat_dv3 = qmmm_lat_dv(mm_cell, particles_mm(qm_atom_index(min_ip))%r, particles_mm(qm_atom_index(max_ip))%r)
460 72 : IF (lat_dv3(i_dim) < 0.0_dp) lat_dv3(i_dim) = lat_dv3(i_dim) + 1.0_dp
461 936 : lat_min = MATMUL(mm_cell%h_inv, particles_mm(qm_atom_index(min_ip))%r)
462 72 : min_coord_lat(i_dim) = lat_min(i_dim)
463 72 : 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 312 : min_coord = MATMUL(mm_cell%hmat, min_coord_lat)
468 312 : max_coord = MATMUL(mm_cell%hmat, max_coord_lat)
469 24 : DEALLOCATE (avoid)
470 : END IF ! pbc aware center
471 3888 : 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 12636 : transl_v(:) = transl_v(:) - SUM(qm_cell%hmat, 2)/2.0_dp
478 :
479 3726 : IF (ANY(qmmm_env%qm%utrasl /= 1.0_dp)) THEN
480 : transl_v = REAL(FLOOR(transl_v/qmmm_env%qm%utrasl), KIND=dp)* &
481 216 : qmmm_env%qm%utrasl
482 : END IF
483 3888 : qmmm_env%qm%transl_v = qmmm_env%qm%transl_v + transl_v
484 972 : particles_mm => subsys_mm%particles%els
485 1150822 : DO ip = 1, subsys_mm%particles%n_els
486 4600372 : particles_mm(ip)%r = particles_mm(ip)%r - transl_v
487 : END DO
488 972 : IF (qmmm_env%qm%added_shells%num_mm_atoms > 0) THEN
489 0 : DO ip = 1, qmmm_env%qm%added_shells%num_mm_atoms
490 0 : qmmm_env%qm%added_shells%added_particles(ip)%r = qmmm_env%qm%added_shells%added_particles(ip)%r - transl_v
491 0 : 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 972 : unit_nr = cp_logger_get_default_io_unit()
495 972 : IF (unit_nr > 0) WRITE (unit=unit_nr, fmt='(/1X,A)') &
496 490 : " Translating the system in order to center the QM fragment in the QM box."
497 972 : IF (.NOT. qmmm_env%qm%center_qm_subsys) qmmm_env%qm%do_translate = .FALSE.
498 : END IF
499 3914 : particles_mm => subsys_mm%particles%els
500 22520 : DO ip = 1, SIZE(qm_atom_index)
501 152762 : particles_qm(ip)%r = particles_mm(qm_atom_index(ip))%r
502 : END DO
503 :
504 3914 : subsys_section => section_vals_get_subs_vals(qmmm_env%qs_env%input, "SUBSYS")
505 :
506 3914 : CALL get_qs_env(qs_env=qmmm_env%qs_env, qs_kind_set=qs_kind_set)
507 3914 : CALL write_qs_particle_coordinates(particles_qm, qs_kind_set, subsys_section, "QM/MM first QM, then MM (0 charges)")
508 11742 : ALLOCATE (charges(SIZE(particles_mm)))
509 1377072 : charges = 0.0_dp
510 3914 : CALL write_fist_particle_coordinates(particles_mm, subsys_section, charges)
511 3914 : DEALLOCATE (charges)
512 :
513 7828 : 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 24 : 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 24 : mean_z = 0.0_dp
532 276 : DO ip = 1, SIZE(qm_atom_index)
533 : mean_z = mean_z + EXP(gaussi*2.0*pi* &
534 4056 : MATMUL(mm_cell%h_inv, particles_mm(qm_atom_index(ip))%r))
535 : END DO
536 96 : mean_z = mean_z/ABS(mean_z)
537 : qmmm_pbc_aware_mean = MATMUL(mm_cell%hmat, &
538 456 : 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 7488 : 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 97344 : lat_v1 = MATMUL(mm_cell%h_inv, p1)
555 97344 : lat_v2 = MATMUL(mm_cell%h_inv, p2)
556 :
557 29952 : qmmm_lat_dv = lat_v2 - lat_v1
558 29952 : 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 1368 : 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 1368 : closest_ip = -1
588 1368 : my_closest_dv = HUGE(0.0)
589 16056 : DO ip = 1, SIZE(qm_atom_index)
590 14688 : IF (avoid(ip)) CYCLE
591 7416 : lat_dv3 = qmmm_lat_dv(mm_cell, p, particles_mm(qm_atom_index(ip))%r)
592 31032 : DO shift = -1, 1
593 22248 : lat_dv_shifted = lat_dv3(i_dim) + shift*1.0_dp
594 36936 : IF (ABS(lat_dv_shifted) < ABS(my_closest_dv) .AND. (dir*lat_dv_shifted >= 0.0)) THEN
595 2330 : my_closest_dv = lat_dv_shifted
596 2330 : closest_ip = ip
597 : END IF
598 : END DO
599 : END DO
600 :
601 1368 : IF (PRESENT(closest_dv)) THEN
602 1296 : closest_dv = my_closest_dv
603 : END IF
604 :
605 1368 : 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 1845816 : 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 7383264 : r = NORM2(rij)
624 1845816 : r0 = spherical_cutoff(1) - 20.0_dp*spherical_cutoff(2)
625 1845816 : factor = 0.5_dp*(1.0_dp - TANH((r - r0)/spherical_cutoff(2)))
626 :
627 1845816 : END SUBROUTINE spherical_cutoff_factor
628 :
629 : END MODULE qmmm_util
|
