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