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