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 : !> \brief Routines used for force-mixing QM/MM calculations
10 : !> \par History
11 : !> 2.2012 created [noam]
12 : !> \author Noam Bernstein
13 : ! **************************************************************************************************
14 : MODULE qmmmx_util
15 : USE atomic_kind_types, ONLY: atomic_kind_type,&
16 : get_atomic_kind,&
17 : set_atomic_kind
18 : USE cell_types, ONLY: cell_copy,&
19 : cell_type,&
20 : pbc
21 : USE cp_log_handling, ONLY: cp_logger_get_default_io_unit,&
22 : cp_logger_get_default_unit_nr
23 : USE cp_subsys_types, ONLY: cp_subsys_get,&
24 : cp_subsys_type
25 : USE fist_neighbor_list_types, ONLY: fist_neighbor_deallocate,&
26 : fist_neighbor_type
27 : USE fist_neighbor_lists, ONLY: build_fist_neighbor_lists
28 : USE input_section_types, ONLY: &
29 : section_vals_add_values, section_vals_duplicate, section_vals_get, &
30 : section_vals_get_subs_vals, section_vals_get_subs_vals3, section_vals_release, &
31 : section_vals_remove_values, section_vals_set_subs_vals, section_vals_type, &
32 : section_vals_val_get, section_vals_val_set, section_vals_write
33 : USE kinds, ONLY: default_string_length,&
34 : dp
35 : USE memory_utilities, ONLY: reallocate
36 : USE message_passing, ONLY: mp_para_env_type
37 : USE molecule_list_types, ONLY: molecule_list_type
38 : USE molecule_types, ONLY: molecule_type
39 : USE particle_list_types, ONLY: particle_list_type
40 : USE particle_types, ONLY: particle_type
41 : USE qmmm_ff_fist, ONLY: qmmm_ff_precond_only_qm
42 : USE qmmm_types, ONLY: qmmm_env_get
43 : USE qmmm_types_low, ONLY: force_mixing_label_QM_core,&
44 : force_mixing_label_QM_core_list,&
45 : force_mixing_label_QM_dynamics,&
46 : force_mixing_label_QM_dynamics_list,&
47 : force_mixing_label_buffer,&
48 : force_mixing_label_buffer_list,&
49 : force_mixing_label_none,&
50 : force_mixing_label_termination
51 : USE qmmm_util, ONLY: apply_qmmm_translate
52 : USE qmmmx_types, ONLY: qmmmx_env_type
53 : #include "./base/base_uses.f90"
54 :
55 : IMPLICIT NONE
56 : PRIVATE
57 :
58 : LOGICAL, PRIVATE :: debug_this_module = .FALSE.
59 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qmmmx_util'
60 :
61 : PUBLIC :: setup_force_mixing_qmmm_sections, update_force_mixing_labels, &
62 : apply_qmmmx_translate
63 :
64 : CONTAINS
65 :
66 : ! **************************************************************************************************
67 : !> \brief Apply translation to the full system in order to center the QM
68 : !> system into the QM box
69 : !> \param qmmmx_env ...
70 : !> \par History
71 : !> 08.2007 created [tlaino] - Zurich University
72 : !> \author Teodoro Laino
73 : ! **************************************************************************************************
74 52 : SUBROUTINE apply_qmmmx_translate(qmmmx_env)
75 : TYPE(qmmmx_env_type), POINTER :: qmmmx_env
76 :
77 : INTEGER :: ip
78 : TYPE(cell_type), POINTER :: cell_core, cell_extended
79 : TYPE(cp_subsys_type), POINTER :: subsys_core, subsys_extended
80 52 : TYPE(particle_type), DIMENSION(:), POINTER :: particles_core, particles_extended
81 :
82 52 : NULLIFY (cell_core, cell_extended)
83 52 : NULLIFY (subsys_core, subsys_extended)
84 52 : NULLIFY (particles_core, particles_extended)
85 :
86 : ! want to center extended, and make core consistent with that
87 52 : CALL apply_qmmm_translate(qmmmx_env%ext)
88 :
89 : ! translate core fist particles
90 52 : CALL qmmm_env_get(qmmmx_env%ext, subsys=subsys_extended)
91 52 : CALL cp_subsys_get(subsys_extended, cell=cell_extended)
92 52 : CALL qmmm_env_get(qmmmx_env%core, subsys=subsys_core)
93 52 : CALL cp_subsys_get(subsys_core, cell=cell_core)
94 52 : particles_extended => subsys_extended%particles%els
95 52 : particles_core => subsys_core%particles%els
96 99838 : DO ip = 1, SIZE(particles_extended)
97 798340 : particles_core(ip)%r = particles_extended(ip)%r
98 : END DO
99 52 : CALL cell_copy(cell_extended, cell_core)
100 :
101 : ! The core QM particles will be updated the regular call
102 : ! to apply_qmmm_translate() from within qmmm_calc_energy_force()
103 :
104 52 : END SUBROUTINE apply_qmmmx_translate
105 :
106 : ! **************************************************************************************************
107 : !> \brief ...
108 : !> \param subsys ...
109 : !> \param qmmm_section ...
110 : !> \param labels_changed ...
111 : !> \par History
112 : !> 02.2012 created [noam]
113 : !> \author Noam Bernstein
114 : ! **************************************************************************************************
115 56 : SUBROUTINE update_force_mixing_labels(subsys, qmmm_section, labels_changed)
116 : TYPE(cp_subsys_type), POINTER :: subsys
117 : TYPE(section_vals_type), POINTER :: qmmm_section
118 : LOGICAL, OPTIONAL :: labels_changed
119 :
120 56 : CHARACTER(LEN=default_string_length), POINTER :: adaptive_exclude_molecules(:)
121 : INTEGER :: i_rep_section, i_rep_val, ip, max_n_qm, n_new, n_rep_exclude, n_rep_section, &
122 : n_rep_val, natoms, output_unit, QM_extended_seed_min_label_val
123 56 : INTEGER, ALLOCATABLE :: new_full_labels(:), orig_full_labels(:)
124 56 : INTEGER, POINTER :: broken_bonds(:), cur_indices(:), &
125 56 : cur_labels(:), mm_index_entry(:), &
126 56 : new_indices(:), new_labels(:)
127 : LOGICAL :: explicit, QM_extended_seed_is_core_list
128 56 : REAL(dp), ALLOCATABLE :: nearest_dist(:)
129 56 : REAL(dp), POINTER :: r_buf(:), r_core(:), r_qm(:)
130 : TYPE(cell_type), POINTER :: cell
131 : TYPE(fist_neighbor_type), POINTER :: nlist
132 : TYPE(molecule_list_type), POINTER :: molecules
133 : TYPE(molecule_type), DIMENSION(:), POINTER :: molecule_set
134 : TYPE(mp_para_env_type), POINTER :: para_env
135 : TYPE(particle_list_type), POINTER :: particles
136 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
137 : TYPE(section_vals_type), POINTER :: force_mixing_section, &
138 : non_adaptive_section, qm_kind_section, &
139 : restart_section
140 :
141 112 : output_unit = cp_logger_get_default_io_unit()
142 :
143 56 : IF (debug_this_module .AND. output_unit > 0) WRITE (output_unit, *) "BOB starting update_force_mixing_labels"
144 : ! get cur indices, labels
145 56 : force_mixing_section => section_vals_get_subs_vals3(qmmm_section, "FORCE_MIXING")
146 56 : CALL get_force_mixing_indices(force_mixing_section, cur_indices, cur_labels)
147 56 : IF (debug_this_module .AND. output_unit > 0) WRITE (output_unit, *) "BOB got cur_indices ", SIZE(cur_indices)
148 56 : IF (debug_this_module .AND. output_unit > 0) WRITE (output_unit, *) "BOB got cur_labels ", SIZE(cur_labels)
149 :
150 : ! read from input
151 : ![NB] breakable bonds will come from here, too
152 56 : NULLIFY (r_core, r_qm, r_buf, adaptive_exclude_molecules, broken_bonds)
153 56 : CALL section_vals_val_get(force_mixing_section, "R_CORE", r_vals=r_core)
154 56 : CALL section_vals_val_get(force_mixing_section, "R_QM", r_vals=r_qm)
155 : CALL section_vals_val_get(force_mixing_section, "QM_EXTENDED_SEED_IS_ONLY_CORE_LIST", &
156 56 : l_val=QM_extended_seed_is_core_list)
157 56 : CALL section_vals_val_get(force_mixing_section, "R_BUF", r_vals=r_buf)
158 56 : CALL section_vals_val_get(force_mixing_section, "MAX_N_QM", i_val=max_n_qm)
159 :
160 56 : CALL section_vals_val_get(force_mixing_section, "ADAPTIVE_EXCLUDE_MOLECULES", n_rep_val=n_rep_exclude)
161 56 : IF (n_rep_exclude > 0) THEN
162 12 : CALL section_vals_val_get(force_mixing_section, "ADAPTIVE_EXCLUDE_MOLECULES", c_vals=adaptive_exclude_molecules)
163 : END IF
164 : ![NB] need to read real list from input
165 : ! should be 2xN_bb integer arrays, with (1,:) indices of inside atoms, and (2,:) indices of outside atoms
166 : ! maybe also breakable_bond_types, with _atomic numbers_ of inside/outside atoms?
167 : ! separate lists for core/buffer?
168 :
169 : ! get particles, molecules
170 56 : NULLIFY (particles, molecules)
171 56 : CALL cp_subsys_get(subsys=subsys, particles=particles, molecules=molecules)
172 56 : particle_set => particles%els
173 56 : molecule_set => molecules%els
174 :
175 56 : natoms = SIZE(particle_set)
176 :
177 : ! initialize new indices, labels, and new_full_labels
178 56 : NULLIFY (new_indices, new_labels)
179 56 : CALL reallocate(new_indices, 1, SIZE(cur_indices))
180 56 : CALL reallocate(new_labels, 1, SIZE(cur_labels))
181 1484 : new_indices = 0
182 1484 : new_labels = force_mixing_label_none
183 168 : ALLOCATE (new_full_labels(natoms))
184 99902 : new_full_labels = force_mixing_label_none
185 :
186 : ! neighbor list for various hysteretic distance calls
187 56 : NULLIFY (cell)
188 56 : CALL cp_subsys_get(subsys, cell=cell)
189 56 : NULLIFY (nlist)
190 56 : CALL make_neighbor_list(force_mixing_section, subsys, cell, MAX(r_core(2), r_qm(2), r_buf(2)), nlist)
191 :
192 : ! create labels for core_list from QM_KIND
193 56 : NULLIFY (mm_index_entry)
194 56 : qm_kind_section => section_vals_get_subs_vals3(qmmm_section, "QM_KIND")
195 56 : CALL section_vals_get(qm_kind_section, n_repetition=n_rep_section)
196 56 : n_new = 0
197 174 : DO i_rep_section = 1, n_rep_section
198 118 : CALL section_vals_val_get(qm_kind_section, "MM_INDEX", i_rep_section=i_rep_section, n_rep_val=n_rep_val)
199 304 : DO i_rep_val = 1, n_rep_val
200 : CALL section_vals_val_get(qm_kind_section, "MM_INDEX", i_rep_section=i_rep_section, i_rep_val=i_rep_val, &
201 130 : i_vals=mm_index_entry)
202 422 : DO ip = 1, SIZE(mm_index_entry)
203 : CALL add_new_label(mm_index_entry(ip), force_mixing_label_QM_core_list, n_new, new_indices, new_labels, &
204 304 : new_full_labels, max_n_qm)
205 : END DO ! ip
206 : END DO ! i_rep_val
207 : END DO ! i_rep_section
208 :
209 56 : IF (debug_this_module .AND. output_unit > 0) THEN
210 0 : WRITE (output_unit, *) "BOB core_list new_indices ", new_indices(1:n_new)
211 0 : WRITE (output_unit, *) "BOB core_list new_labels ", new_labels(1:n_new)
212 : END IF
213 :
214 : ! create labels for non adaptive QM and buffer regions from *_NON_ADAPTIVE&QM_KIND sections
215 : non_adaptive_section => section_vals_get_subs_vals(qmmm_section, "FORCE_MIXING%QM_NON_ADAPTIVE", &
216 56 : can_return_null=.TRUE.)
217 56 : IF (ASSOCIATED(non_adaptive_section)) THEN
218 56 : qm_kind_section => section_vals_get_subs_vals3(non_adaptive_section, "QM_KIND")
219 56 : CALL section_vals_get(qm_kind_section, n_repetition=n_rep_section)
220 86 : DO i_rep_section = 1, n_rep_section
221 30 : CALL section_vals_val_get(qm_kind_section, "MM_INDEX", i_rep_section=i_rep_section, n_rep_val=n_rep_val)
222 164 : DO i_rep_val = 1, n_rep_val
223 : CALL section_vals_val_get(qm_kind_section, "MM_INDEX", i_rep_section=i_rep_section, i_rep_val=i_rep_val, &
224 78 : i_vals=mm_index_entry)
225 186 : DO ip = 1, SIZE(mm_index_entry)
226 : CALL add_new_label(mm_index_entry(ip), force_mixing_label_QM_dynamics_list, n_new, new_indices, new_labels, &
227 156 : new_full_labels, max_n_qm)
228 : END DO ! ip
229 : END DO ! i_rep_val
230 : END DO ! i_rep_section
231 : END IF
232 56 : IF (debug_this_module .AND. output_unit > 0) THEN
233 0 : WRITE (output_unit, *) "BOB core_list + non adaptive QM new_indices ", new_indices(1:n_new)
234 0 : WRITE (output_unit, *) "BOB core_list + non adaptive QM new_labels ", new_labels(1:n_new)
235 : END IF
236 : non_adaptive_section => section_vals_get_subs_vals(qmmm_section, "FORCE_MIXING%BUFFER_NON_ADAPTIVE", &
237 56 : can_return_null=.TRUE.)
238 56 : IF (ASSOCIATED(non_adaptive_section)) THEN
239 56 : qm_kind_section => section_vals_get_subs_vals3(non_adaptive_section, "QM_KIND")
240 56 : CALL section_vals_get(qm_kind_section, n_repetition=n_rep_section)
241 92 : DO i_rep_section = 1, n_rep_section
242 36 : CALL section_vals_val_get(qm_kind_section, "MM_INDEX", i_rep_section=i_rep_section, n_rep_val=n_rep_val)
243 692 : DO i_rep_val = 1, n_rep_val
244 : CALL section_vals_val_get(qm_kind_section, "MM_INDEX", i_rep_section=i_rep_section, i_rep_val=i_rep_val, &
245 600 : i_vals=mm_index_entry)
246 1236 : DO ip = 1, SIZE(mm_index_entry)
247 : CALL add_new_label(mm_index_entry(ip), force_mixing_label_buffer_list, n_new, new_indices, new_labels, &
248 1200 : new_full_labels, max_n_qm)
249 : END DO ! ip
250 : END DO ! i_rep_val
251 : END DO ! i_rep_section
252 : END IF
253 :
254 56 : IF (debug_this_module .AND. output_unit > 0) THEN
255 0 : WRITE (output_unit, *) "BOB core_list + non adaptive QM+buffer new_indices ", new_indices(1:n_new)
256 0 : WRITE (output_unit, *) "BOB core_list + non adaptive QM+buffer new_labels ", new_labels(1:n_new)
257 : END IF
258 :
259 : ! allocate and initialize full atom set labels for hysteretic loops
260 168 : ALLOCATE (nearest_dist(natoms))
261 :
262 : ! orig_full_labels is full array (natoms) with orig labels
263 112 : ALLOCATE (orig_full_labels(natoms))
264 99902 : orig_full_labels = force_mixing_label_none
265 1484 : orig_full_labels(cur_indices(:)) = cur_labels(:)
266 :
267 : ! hysteretically set QM core from QM_core_list and radii, whole molecule
268 : ![NB] need to replace all the whole molecule stuff with pad to breakable bonds. not quite done
269 : ! (need intra molecule bond info, which isn't available for QM molecules yet)
270 :
271 : ! add core using hysteretic selection(core_list, r_core) + unbreakable bonds
272 : CALL add_layer_hysteretically( &
273 : nlist, particle_set, cell, nearest_dist, &
274 : orig_full_labels, new_full_labels, n_new, new_indices, new_labels, &
275 : force_mixing_label_QM_core_list, force_mixing_label_QM_core_list, force_mixing_label_QM_core, r_core, &
276 56 : max_n_qm, adaptive_exclude_molecules, molecule_set, broken_bonds)
277 : ![NB] should actually pass this back for making link sections?
278 56 : DEALLOCATE (broken_bonds)
279 :
280 56 : IF (debug_this_module .AND. output_unit > 0) THEN
281 0 : WRITE (output_unit, *) "BOB core new_indices ", new_indices(1:n_new)
282 0 : WRITE (output_unit, *) "BOB core new_labels ", new_labels(1:n_new)
283 : END IF
284 :
285 : ![NB] need more sophisticated QM extended, buffer rules
286 :
287 : ! add QM using hysteretic selection (core_list, r_qm) + unbreakable bonds
288 56 : IF (debug_this_module .AND. output_unit > 0) THEN
289 0 : WRITE (output_unit, *) "BOB QM_extended_seed_is_core_list ", QM_extended_seed_is_core_list
290 : END IF
291 56 : IF (QM_extended_seed_is_core_list) THEN
292 0 : QM_extended_seed_min_label_val = force_mixing_label_QM_core_list
293 : ELSE ! QM region seed is all of core, not just core list + unbreakable bonds
294 56 : QM_extended_seed_min_label_val = force_mixing_label_QM_core
295 : END IF
296 : CALL add_layer_hysteretically(nlist, particle_set, cell, nearest_dist, &
297 : orig_full_labels, new_full_labels, n_new, new_indices, new_labels, &
298 : QM_extended_seed_min_label_val, force_mixing_label_QM_core_list, &
299 : force_mixing_label_QM_dynamics, r_qm, &
300 56 : max_n_qm, adaptive_exclude_molecules, molecule_set)
301 :
302 56 : IF (debug_this_module .AND. output_unit > 0) THEN
303 0 : WRITE (output_unit, *) "BOB extended new_indices ", new_indices(1:n_new)
304 0 : WRITE (output_unit, *) "BOB extended new_labels ", new_labels(1:n_new)
305 : END IF
306 :
307 : ! add buffer using hysteretic selection (>= QM extended, r_buf) + unbreakable bonds
308 : CALL add_layer_hysteretically( &
309 : nlist, particle_set, cell, nearest_dist, &
310 : orig_full_labels, new_full_labels, n_new, new_indices, new_labels, &
311 : force_mixing_label_QM_dynamics, force_mixing_label_QM_core_list, force_mixing_label_buffer, r_buf, &
312 56 : max_n_qm, adaptive_exclude_molecules, molecule_set, broken_bonds)
313 : ![NB] should actually pass this back for making link sections?
314 56 : DEALLOCATE (broken_bonds)
315 :
316 56 : IF (debug_this_module .AND. output_unit > 0) THEN
317 0 : WRITE (output_unit, *) "BOB buffer new_indices ", new_indices(1:n_new)
318 0 : WRITE (output_unit, *) "BOB buffer new_labels ", new_labels(1:n_new)
319 : END IF
320 :
321 56 : DEALLOCATE (nearest_dist)
322 :
323 65568 : IF (PRESENT(labels_changed)) labels_changed = ANY(new_full_labels /= orig_full_labels)
324 :
325 56 : DEALLOCATE (orig_full_labels)
326 56 : DEALLOCATE (new_full_labels)
327 :
328 : ! reduce new indices, labels to actually used size
329 56 : CALL reallocate(new_indices, 1, n_new)
330 56 : CALL reallocate(new_labels, 1, n_new)
331 :
332 : ! save info in input structure
333 56 : restart_section => section_vals_get_subs_vals(qmmm_section, "FORCE_MIXING%RESTART_INFO")
334 56 : CALL section_vals_get(restart_section, explicit=explicit)
335 56 : IF (explicit) CALL section_vals_remove_values(restart_section)
336 56 : CALL section_vals_val_set(restart_section, "INDICES", i_vals_ptr=new_indices)
337 56 : CALL section_vals_val_set(restart_section, "LABELS", i_vals_ptr=new_labels)
338 :
339 56 : DEALLOCATE (cur_indices, cur_labels)
340 56 : CALL fist_neighbor_deallocate(nlist)
341 :
342 : ![NB] perhap be controlled by some &PRINT section?
343 56 : CALL cp_subsys_get(subsys, para_env=para_env)
344 56 : IF (para_env%is_source() .AND. output_unit > 0) THEN
345 : WRITE (unit=output_unit, fmt='(A,A,I6,A,I5,A,I5,A,I5)') &
346 28 : "QMMM FORCE MIXING final count (not including links): ", &
347 1021 : " N_QM core_list ", COUNT(new_labels == force_mixing_label_QM_core_list), &
348 1021 : " N_QM core ", COUNT(new_labels == force_mixing_label_QM_core), &
349 28 : " N_QM extended ", COUNT(new_labels == force_mixing_label_QM_dynamics .OR. &
350 1021 : new_labels == force_mixing_label_QM_dynamics_list), &
351 28 : " N_QM buffered ", COUNT(new_labels == force_mixing_label_buffer .OR. &
352 1049 : new_labels == force_mixing_label_buffer_list)
353 : END IF
354 :
355 280 : END SUBROUTINE update_force_mixing_labels
356 :
357 : ! **************************************************************************************************
358 : !> \brief ...
359 : !> \param ip ...
360 : !> \param label ...
361 : !> \param n_new ...
362 : !> \param new_indices ...
363 : !> \param new_labels ...
364 : !> \param new_full_labels ...
365 : !> \param max_n_qm ...
366 : ! **************************************************************************************************
367 1986 : SUBROUTINE add_new_label(ip, label, n_new, new_indices, new_labels, new_full_labels, max_n_qm)
368 : INTEGER :: ip, label, n_new
369 : INTEGER, POINTER :: new_indices(:), new_labels(:)
370 : INTEGER :: new_full_labels(:), max_n_qm
371 :
372 : INTEGER :: i, old_index
373 :
374 1986 : IF (new_full_labels(ip) > force_mixing_label_none) THEN ! already marked, just change mark
375 0 : old_index = -1
376 0 : DO i = 1, n_new
377 0 : IF (new_indices(i) == ip) THEN
378 : old_index = i
379 : EXIT
380 : END IF
381 : END DO
382 0 : IF (old_index <= 0) THEN
383 : CALL cp_abort(__LOCATION__, &
384 : "add_new_label found atom with a label "// &
385 0 : "already set, but not in new_indices array")
386 : END IF
387 0 : new_labels(old_index) = label
388 : ELSE
389 1986 : n_new = n_new + 1
390 1986 : IF (n_new > max_n_qm) THEN
391 : CALL cp_abort(__LOCATION__, &
392 : "add_new_label tried to add more atoms "// &
393 0 : "than allowed by &FORCE_MIXING&MAX_N_QM!")
394 : END IF
395 1986 : IF (n_new > SIZE(new_indices)) CALL reallocate(new_indices, 1, n_new + 9)
396 1986 : IF (n_new > SIZE(new_labels)) CALL reallocate(new_labels, 1, n_new + 9)
397 1986 : new_indices(n_new) = ip
398 1986 : new_labels(n_new) = label
399 : END IF
400 1986 : new_full_labels(ip) = label
401 1986 : END SUBROUTINE add_new_label
402 :
403 : ! **************************************************************************************************
404 : !> \brief ...
405 : !> \param nlist ...
406 : !> \param particle_set ...
407 : !> \param cell ...
408 : !> \param nearest_dist ...
409 : !> \param orig_full_labels ...
410 : !> \param new_full_labels ...
411 : !> \param n_new ...
412 : !> \param new_indices ...
413 : !> \param new_labels ...
414 : !> \param seed_min_label_val ...
415 : !> \param seed_max_label_val ...
416 : !> \param set_label_val ...
417 : !> \param r_inout ...
418 : !> \param max_n_qm ...
419 : !> \param adaptive_exclude_molecules ...
420 : !> \param molecule_set ...
421 : !> \param broken_bonds ...
422 : ! **************************************************************************************************
423 168 : SUBROUTINE add_layer_hysteretically(nlist, particle_set, cell, nearest_dist, &
424 168 : orig_full_labels, new_full_labels, n_new, new_indices, new_labels, &
425 : seed_min_label_val, seed_max_label_val, set_label_val, r_inout, max_n_qm, &
426 : adaptive_exclude_molecules, molecule_set, broken_bonds)
427 : TYPE(fist_neighbor_type), POINTER :: nlist
428 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
429 : TYPE(cell_type), POINTER :: cell
430 : REAL(dp) :: nearest_dist(:)
431 : INTEGER :: orig_full_labels(:), new_full_labels(:), &
432 : n_new
433 : INTEGER, POINTER :: new_indices(:), new_labels(:)
434 : INTEGER :: seed_min_label_val, seed_max_label_val, &
435 : set_label_val
436 : REAL(dp) :: r_inout(2)
437 : INTEGER :: max_n_qm
438 : CHARACTER(len=*), POINTER :: adaptive_exclude_molecules(:)
439 : TYPE(molecule_type), DIMENSION(:), OPTIONAL, &
440 : POINTER :: molecule_set
441 : INTEGER, OPTIONAL, POINTER :: broken_bonds(:)
442 :
443 : INTEGER :: i_ind, im, im_exclude, ip, ipair, &
444 : ipairkind, j_ind, output_unit
445 : LOGICAL :: adaptive_exclude, i_in_new_seed, i_outside_new_seed, j_in_new_seed, &
446 : j_outside_new_seed, molec_in_inner, molec_in_outer
447 : REAL(dp) :: r_ij(3), r_ij_mag
448 :
449 168 : output_unit = cp_logger_get_default_unit_nr()
450 :
451 168 : IF (debug_this_module .AND. output_unit > 0) WRITE (output_unit, *) "BOB adding hysteretically seed ", &
452 0 : seed_min_label_val, seed_max_label_val, " set ", set_label_val, " r ", r_inout
453 : ! calculate nearest dist from each atom outside of new seed to nearest atom inside of new seed
454 299706 : nearest_dist = HUGE(1.0_dp)
455 : ! loop over pairs of all kinds in random order
456 4704 : DO ipairkind = 1, SIZE(nlist%neighbor_kind_pairs)
457 2621184 : DO ipair = 1, nlist%neighbor_kind_pairs(ipairkind)%npairs
458 :
459 2616480 : i_ind = nlist%neighbor_kind_pairs(ipairkind)%list(1, ipair)
460 2616480 : j_ind = nlist%neighbor_kind_pairs(ipairkind)%list(2, ipair)
461 :
462 2616480 : i_in_new_seed = (new_full_labels(i_ind) >= seed_min_label_val .AND. new_full_labels(i_ind) <= seed_max_label_val)
463 2616480 : i_outside_new_seed = (new_full_labels(i_ind) < seed_min_label_val)
464 2616480 : j_in_new_seed = (new_full_labels(j_ind) >= seed_min_label_val .AND. new_full_labels(j_ind) <= seed_max_label_val)
465 2616480 : j_outside_new_seed = (new_full_labels(j_ind) < seed_min_label_val)
466 :
467 2621016 : IF ((i_in_new_seed .AND. j_outside_new_seed) .OR. (j_in_new_seed .AND. i_outside_new_seed)) THEN
468 22936 : r_ij = pbc(particle_set(i_ind)%r - particle_set(j_ind)%r, cell)
469 22936 : r_ij_mag = SQRT(SUM(r_ij**2))
470 5734 : IF (i_in_new_seed .AND. j_outside_new_seed .AND. (r_ij_mag < nearest_dist(j_ind))) THEN
471 1694 : nearest_dist(j_ind) = r_ij_mag
472 : END IF
473 5734 : IF (j_in_new_seed .AND. i_outside_new_seed .AND. (r_ij_mag < nearest_dist(i_ind))) THEN
474 1660 : nearest_dist(i_ind) = r_ij_mag
475 : END IF
476 : END IF
477 :
478 : END DO
479 : END DO
480 :
481 : ![NB] this is whole molecule. Should be replaced with labeling of individual atoms +
482 : ! pad_to_breakable_bonds (below), but QM molecule bond information isn't available yet
483 90558 : DO im = 1, SIZE(molecule_set)
484 : ! molecule_set(im)%first_atom,molecule_set(im)%last_atom
485 90390 : IF (ASSOCIATED(adaptive_exclude_molecules)) THEN
486 89730 : adaptive_exclude = .FALSE.
487 179460 : DO im_exclude = 1, SIZE(adaptive_exclude_molecules)
488 177120 : IF (TRIM(molecule_set(im)%molecule_kind%name) == TRIM(adaptive_exclude_molecules(im_exclude)) .OR. &
489 89730 : TRIM(molecule_set(im)%molecule_kind%name) == '_QM_'//TRIM(adaptive_exclude_molecules(im_exclude))) THEN
490 2340 : adaptive_exclude = .TRUE.
491 : END IF
492 : END DO
493 89730 : IF (adaptive_exclude) CYCLE
494 : END IF
495 350854 : molec_in_inner = ANY(nearest_dist(molecule_set(im)%first_atom:molecule_set(im)%last_atom) <= r_inout(1))
496 350504 : molec_in_outer = ANY(nearest_dist(molecule_set(im)%first_atom:molecule_set(im)%last_atom) <= r_inout(2))
497 88218 : IF (molec_in_inner) THEN
498 1480 : DO ip = molecule_set(im)%first_atom, molecule_set(im)%last_atom
499 : ! labels are being rebuild from scratch, so never overwrite new label that's higher level
500 1480 : IF (new_full_labels(ip) < set_label_val) THEN
501 1110 : CALL add_new_label(ip, set_label_val, n_new, new_indices, new_labels, new_full_labels, max_n_qm)
502 : END IF
503 : END DO
504 87680 : ELSE IF (molec_in_outer) THEN
505 520 : IF (ANY(orig_full_labels(molecule_set(im)%first_atom:molecule_set(im)%last_atom) >= set_label_val)) THEN
506 32 : DO ip = molecule_set(im)%first_atom, molecule_set(im)%last_atom
507 : ! labels are being rebuild from scratch, so never overwrite new label that's higher level
508 32 : IF (new_full_labels(ip) < set_label_val) THEN
509 24 : CALL add_new_label(ip, set_label_val, n_new, new_indices, new_labels, new_full_labels, max_n_qm)
510 : END IF
511 : END DO
512 : END IF
513 : END IF
514 : END DO
515 168 : IF (PRESENT(broken_bonds)) CALL reallocate(broken_bonds, 1, 0)
516 :
517 168 : END SUBROUTINE add_layer_hysteretically
518 :
519 : ! **************************************************************************************************
520 : !> \brief ...
521 : !> \param force_mixing_section ...
522 : !> \param subsys ...
523 : !> \param cell ...
524 : !> \param r_max ...
525 : !> \param nlist ...
526 : ! **************************************************************************************************
527 56 : SUBROUTINE make_neighbor_list(force_mixing_section, subsys, cell, r_max, nlist)
528 : TYPE(section_vals_type), POINTER :: force_mixing_section
529 : TYPE(cp_subsys_type), POINTER :: subsys
530 : TYPE(cell_type), POINTER :: cell
531 : REAL(dp) :: r_max
532 : TYPE(fist_neighbor_type), POINTER :: nlist
533 :
534 : CHARACTER(LEN=default_string_length) :: kind_name
535 56 : CHARACTER(LEN=default_string_length), POINTER :: kind_name_a(:)
536 : INTEGER :: ik
537 : LOGICAL :: skip_kind
538 56 : REAL(dp), ALLOCATABLE :: r_max_a(:, :), r_minsq_a(:, :)
539 : TYPE(atomic_kind_type), POINTER :: atomic_kind
540 :
541 224 : ALLOCATE (r_max_a(SIZE(subsys%atomic_kinds%els), SIZE(subsys%atomic_kinds%els)))
542 168 : ALLOCATE (r_minsq_a(SIZE(subsys%atomic_kinds%els), SIZE(subsys%atomic_kinds%els)))
543 472612 : r_max_a = r_max
544 472612 : r_minsq_a = EPSILON(1.0_dp)
545 :
546 : ! save kind names
547 168 : ALLOCATE (kind_name_a(SIZE(subsys%atomic_kinds%els)))
548 2012 : DO ik = 1, SIZE(subsys%atomic_kinds%els)
549 1956 : atomic_kind => subsys%atomic_kinds%els(ik)
550 1956 : CALL get_atomic_kind(atomic_kind=atomic_kind, name=kind_name)
551 2012 : kind_name_a(ik) = kind_name
552 : END DO
553 :
554 : ! overwrite kind names so that none are QM, and so excluding QM-QM interactions
555 : ! (which is not what we want) will not happen
556 2012 : DO ik = 1, SIZE(subsys%atomic_kinds%els)
557 1956 : atomic_kind => subsys%atomic_kinds%els(ik)
558 1956 : CALL get_atomic_kind(atomic_kind=atomic_kind, name=kind_name)
559 : ! when atom is QM atom, kind_name is replaced with original
560 : ! mm kind name, and return status is logical .TRUE.
561 1956 : skip_kind = qmmm_ff_precond_only_qm(kind_name)
562 2012 : CALL set_atomic_kind(atomic_kind=atomic_kind, name=kind_name)
563 : END DO
564 :
565 56 : NULLIFY (nlist)
566 : CALL build_fist_neighbor_lists(subsys%atomic_kinds%els, subsys%particles%els, &
567 : cell=cell, r_max=r_max_a, r_minsq=r_minsq_a, &
568 : ei_scale14=1.0_dp, vdw_scale14=1.0_dp, nonbonded=nlist, &
569 : para_env=subsys%para_env, build_from_scratch=.TRUE., geo_check=.FALSE., &
570 56 : mm_section=force_mixing_section)
571 :
572 56 : DEALLOCATE (r_max_a, r_minsq_a)
573 :
574 : ! restore kind names
575 2012 : DO ik = 1, SIZE(subsys%atomic_kinds%els)
576 2012 : CALL set_atomic_kind(atomic_kind=atomic_kind, name=kind_name_a(ik))
577 : END DO
578 56 : DEALLOCATE (kind_name_a)
579 :
580 56 : END SUBROUTINE make_neighbor_list
581 :
582 : ! **************************************************************************************************
583 : !> \brief ...
584 : !> \param subsys ...
585 : !> \param qmmm_section ...
586 : !> \param qmmm_core_section ...
587 : !> \param qmmm_extended_section ...
588 : !> \par History
589 : !> 02.2012 created [noam]
590 : !> \author Noam Bernstein
591 : ! **************************************************************************************************
592 30 : SUBROUTINE setup_force_mixing_qmmm_sections(subsys, qmmm_section, qmmm_core_section, qmmm_extended_section)
593 : TYPE(cp_subsys_type), POINTER :: subsys
594 : TYPE(section_vals_type), POINTER :: qmmm_section, qmmm_core_section, &
595 : qmmm_extended_section
596 :
597 30 : CHARACTER(len=default_string_length), POINTER :: elem_mapping(:, :), elem_mapping_entry(:)
598 : INTEGER :: delta_charge, i_rep_section_core, i_rep_section_extended, i_rep_val_core, &
599 : i_rep_val_extended, ielem, ip, n_elements, output_unit
600 30 : INTEGER, POINTER :: cur_indices(:), cur_labels(:)
601 : LOGICAL :: mapped, new_element_core, &
602 : new_element_extended
603 30 : TYPE(particle_type), DIMENSION(:), POINTER :: particles
604 : TYPE(section_vals_type), POINTER :: buffer_non_adaptive_section, &
605 : dup_link_section, &
606 : force_mixing_section, link_section, &
607 : qm_kind_section
608 :
609 30 : NULLIFY (qmmm_core_section, qmmm_extended_section)
610 60 : output_unit = cp_logger_get_default_unit_nr()
611 :
612 : ! create new qmmm sections for core and extended
613 30 : CALL section_vals_duplicate(qmmm_section, qmmm_core_section)
614 30 : CALL section_vals_duplicate(qmmm_section, qmmm_extended_section)
615 :
616 : ! remove LINKs (specified by user for core) from extended
617 30 : link_section => section_vals_get_subs_vals(qmmm_extended_section, "LINK", can_return_null=.TRUE.)
618 30 : IF (ASSOCIATED(link_section)) THEN
619 30 : CALL section_vals_remove_values(link_section)
620 : END IF
621 : ! for LINKs to be added to extended
622 : buffer_non_adaptive_section => section_vals_get_subs_vals(qmmm_extended_section, "FORCE_MIXING%BUFFER_NON_ADAPTIVE", &
623 30 : can_return_null=.TRUE.)
624 30 : link_section => section_vals_get_subs_vals(buffer_non_adaptive_section, "LINK", can_return_null=.TRUE.)
625 30 : IF (ASSOCIATED(link_section)) THEN
626 30 : NULLIFY (dup_link_section)
627 30 : CALL section_vals_duplicate(link_section, dup_link_section)
628 30 : CALL section_vals_set_subs_vals(qmmm_extended_section, "LINK", dup_link_section)
629 30 : CALL section_vals_release(dup_link_section)
630 : END IF
631 :
632 30 : IF (debug_this_module .AND. output_unit > 0) THEN
633 0 : link_section => section_vals_get_subs_vals(qmmm_core_section, "LINK", can_return_null=.TRUE.)
634 0 : WRITE (output_unit, *) "core section has LINKs ", ASSOCIATED(link_section)
635 0 : CALL section_vals_write(link_section, unit_nr=6)
636 0 : link_section => section_vals_get_subs_vals(qmmm_extended_section, "LINK", can_return_null=.TRUE.)
637 0 : WRITE (output_unit, *) "extended section has LINKs ", ASSOCIATED(link_section)
638 0 : CALL section_vals_write(link_section, unit_nr=6)
639 : END IF
640 :
641 30 : force_mixing_section => section_vals_get_subs_vals(qmmm_section, "FORCE_MIXING")
642 :
643 : ! get QM_KIND_ELEMENT_MAPPING
644 30 : CALL section_vals_val_get(force_mixing_section, "QM_KIND_ELEMENT_MAPPING", n_rep_val=n_elements)
645 90 : ALLOCATE (elem_mapping(2, n_elements))
646 102 : DO ielem = 1, n_elements
647 72 : CALL section_vals_val_get(force_mixing_section, "QM_KIND_ELEMENT_MAPPING", i_rep_val=ielem, c_vals=elem_mapping_entry)
648 390 : elem_mapping(1:2, ielem) = elem_mapping_entry(1:2)
649 : END DO
650 :
651 : ! get CUR_INDICES, CUR_LABELS
652 30 : CALL get_force_mixing_indices(force_mixing_section, cur_indices, cur_labels)
653 30 : IF (SIZE(cur_indices) <= 0) THEN
654 0 : CPABORT("cur_indices is empty, found no QM atoms")
655 : END IF
656 :
657 30 : IF (debug_this_module .AND. output_unit > 0) THEN
658 0 : WRITE (output_unit, *) "cur_indices ", cur_indices
659 0 : WRITE (output_unit, *) "cur_labels ", cur_labels
660 : END IF
661 :
662 : ! loop through elements and atoms, and set up new QM_KIND sections
663 30 : particles => subsys%particles%els
664 :
665 1014 : DO ip = 1, SIZE(cur_indices)
666 984 : IF (cur_labels(ip) > force_mixing_label_none .AND. cur_labels(ip) < force_mixing_label_QM_core_list .AND. &
667 30 : cur_labels(ip) /= force_mixing_label_termination) THEN
668 892 : mapped = .FALSE.
669 1472 : DO ielem = 1, n_elements
670 1472 : IF (TRIM(particles(cur_indices(ip))%atomic_kind%element_symbol) == TRIM(elem_mapping(1, ielem))) THEN
671 : mapped = .TRUE.
672 : EXIT
673 : END IF
674 : END DO
675 892 : IF (.NOT. mapped) THEN
676 : CALL cp_abort(__LOCATION__, &
677 : "Force-mixing failed to find QM_KIND mapping for atom of type "// &
678 : TRIM(particles(cur_indices(ip))%atomic_kind%element_symbol)// &
679 0 : "! ")
680 : END IF
681 : END IF
682 : END DO
683 :
684 : ! pre-existing QM_KIND section specifies list of core atom
685 30 : qm_kind_section => section_vals_get_subs_vals3(qmmm_section, "QM_KIND")
686 30 : CALL section_vals_get(qm_kind_section, n_repetition=i_rep_section_core)
687 30 : IF (i_rep_section_core <= 0) THEN
688 : CALL cp_abort(__LOCATION__, &
689 : "Force-mixing QM didn't find any QM_KIND sections, "// &
690 0 : "so no core specified!")
691 : END IF
692 30 : i_rep_section_extended = i_rep_section_core
693 102 : DO ielem = 1, n_elements
694 : new_element_core = .TRUE.
695 : new_element_extended = .TRUE.
696 3522 : DO ip = 1, SIZE(cur_indices) ! particles with label
697 3420 : IF (TRIM(particles(cur_indices(ip))%atomic_kind%element_symbol) /= TRIM(elem_mapping(1, ielem))) CYCLE
698 : ! extended
699 : ! if current particle is some sort of QM atom, and not in core list
700 : ! (those the user gave explicit QM_KIND sections for), and not a
701 : ! termination atom, need to make a QM_KIND section for it
702 : IF (cur_labels(ip) > force_mixing_label_none .AND. &
703 984 : cur_labels(ip) /= force_mixing_label_QM_core_list .AND. &
704 : cur_labels(ip) /= force_mixing_label_termination) THEN
705 892 : qm_kind_section => section_vals_get_subs_vals3(qmmm_extended_section, "QM_KIND")
706 892 : IF (new_element_extended) THEN ! add new QM_KIND section for this element
707 72 : i_rep_section_extended = i_rep_section_extended + 1
708 72 : CALL section_vals_add_values(qm_kind_section)
709 : CALL section_vals_val_set(qm_kind_section, "_SECTION_PARAMETERS_", i_rep_section=i_rep_section_extended, &
710 72 : c_val=elem_mapping(2, ielem))
711 72 : i_rep_val_extended = 0
712 72 : new_element_extended = .FALSE.
713 : END IF
714 892 : i_rep_val_extended = i_rep_val_extended + 1
715 : CALL section_vals_val_set(qm_kind_section, "MM_INDEX", i_rep_section=i_rep_section_extended, &
716 892 : i_rep_val=i_rep_val_extended, i_val=cur_indices(ip))
717 : END IF ! is a non-termination QM atom
718 :
719 : ! core
720 : ! if current particle is a core QM atom, and not in core list (those the user
721 : ! gave explicit QM_KIND sections for, need to make a QM_KIND section for it
722 1056 : IF (cur_labels(ip) == force_mixing_label_QM_core) THEN
723 78 : qm_kind_section => section_vals_get_subs_vals3(qmmm_core_section, "QM_KIND")
724 78 : IF (new_element_core) THEN ! add new QM_KIND section for this element
725 48 : i_rep_section_core = i_rep_section_core + 1
726 48 : CALL section_vals_add_values(qm_kind_section)
727 : CALL section_vals_val_set(qm_kind_section, "_SECTION_PARAMETERS_", i_rep_section=i_rep_section_core, &
728 48 : c_val=elem_mapping(2, ielem))
729 48 : i_rep_val_core = 0
730 48 : new_element_core = .FALSE.
731 : END IF
732 78 : i_rep_val_core = i_rep_val_core + 1
733 : CALL section_vals_val_set(qm_kind_section, "MM_INDEX", i_rep_section=i_rep_section_core, &
734 78 : i_rep_val=i_rep_val_core, i_val=cur_indices(ip))
735 : END IF ! is a non-termination QM atom
736 :
737 : END DO ! atom index ip
738 : END DO ! element index ielem
739 :
740 30 : CALL section_vals_val_get(force_mixing_section, "EXTENDED_DELTA_CHARGE", i_val=delta_charge)
741 30 : CALL section_vals_val_set(qmmm_extended_section, "DELTA_CHARGE", i_val=delta_charge)
742 :
743 : ![NB] check
744 30 : DEALLOCATE (elem_mapping, cur_indices, cur_labels)
745 :
746 30 : IF (debug_this_module .AND. output_unit > 0) THEN
747 0 : WRITE (output_unit, *) "qmmm_core_section"
748 0 : CALL section_vals_write(qmmm_core_section, unit_nr=6)
749 0 : WRITE (output_unit, *) "qmmm_extended_section"
750 0 : CALL section_vals_write(qmmm_extended_section, unit_nr=6)
751 : END IF
752 :
753 120 : END SUBROUTINE setup_force_mixing_qmmm_sections
754 :
755 : ! **************************************************************************************************
756 : !> \brief ...
757 : !> \param force_mixing_section ...
758 : !> \param indices ...
759 : !> \param labels ...
760 : ! **************************************************************************************************
761 86 : SUBROUTINE get_force_mixing_indices(force_mixing_section, indices, labels)
762 : TYPE(section_vals_type), POINTER :: force_mixing_section
763 : INTEGER, POINTER :: indices(:), labels(:)
764 :
765 : INTEGER :: i_rep_val, n_indices, n_labels, n_reps
766 86 : INTEGER, POINTER :: indices_entry(:), labels_entry(:)
767 : LOGICAL :: explicit
768 : TYPE(section_vals_type), POINTER :: restart_section
769 :
770 86 : NULLIFY (indices, labels)
771 172 : restart_section => section_vals_get_subs_vals(force_mixing_section, "RESTART_INFO")
772 86 : CALL section_vals_get(restart_section, explicit=explicit)
773 86 : IF (.NOT. explicit) THEN ! no old indices, labels, return empty arrays
774 8 : ALLOCATE (indices(0))
775 8 : ALLOCATE (labels(0))
776 8 : RETURN
777 : END IF
778 :
779 : ![NB] maybe switch to reallocatable array
780 78 : CALL section_vals_val_get(restart_section, "INDICES", n_rep_val=n_reps)
781 78 : n_indices = 0
782 156 : DO i_rep_val = 1, n_reps
783 : CALL section_vals_val_get(restart_section, "INDICES", &
784 78 : i_rep_val=i_rep_val, i_vals=indices_entry)
785 156 : n_indices = n_indices + SIZE(indices_entry)
786 : END DO
787 234 : ALLOCATE (indices(n_indices))
788 78 : n_indices = 0
789 156 : DO i_rep_val = 1, n_reps
790 : CALL section_vals_val_get(restart_section, "INDICES", &
791 78 : i_rep_val=i_rep_val, i_vals=indices_entry)
792 4902 : indices(n_indices + 1:n_indices + SIZE(indices_entry)) = indices_entry
793 156 : n_indices = n_indices + SIZE(indices_entry)
794 : END DO
795 :
796 78 : CALL section_vals_val_get(restart_section, "LABELS", n_rep_val=n_reps)
797 78 : n_labels = 0
798 156 : DO i_rep_val = 1, n_reps
799 : CALL section_vals_val_get(restart_section, "LABELS", &
800 78 : i_rep_val=i_rep_val, i_vals=labels_entry)
801 156 : n_labels = n_labels + SIZE(labels_entry)
802 : END DO
803 234 : ALLOCATE (labels(n_labels))
804 78 : n_labels = 0
805 156 : DO i_rep_val = 1, n_reps
806 : CALL section_vals_val_get(restart_section, "LABELS", &
807 78 : i_rep_val=i_rep_val, i_vals=labels_entry)
808 4902 : labels(n_labels + 1:n_labels + SIZE(labels_entry)) = labels_entry
809 156 : n_labels = n_labels + SIZE(labels_entry)
810 : END DO
811 :
812 78 : IF (n_indices /= n_labels) THEN
813 0 : CPABORT("got unequal numbers of force_mixing indices and labels!")
814 : END IF
815 242 : END SUBROUTINE get_force_mixing_indices
816 :
817 : END MODULE qmmmx_util
|
